Integration of Aeromagnetic Data and Landsat Imagery for Structural Analysis Purposes: A Case Study in the Southern Part of Jordan ()
1. Introduction
In this study, several digital image-processing techniques were used to enhance the ability of structural interpretation of Landsat ETM+ images. These techniques were integrated with field geophysical exploration using Aeromagnetic Techniques.
Aerial magnetic or aeromagnetic survey applications are very well known in a wide range of geological and exploration studies and they play a very significance role in demarcating the lithological contacts, local and regional geological structures, faults, dykes, lineaments, exploration of mineral ores, hydrocarbon reservoirs and their associated structures, depth to basement, crustal and tectonic studies.
Numerous geophysical studies based on the analysis and processing of potential field have been carried out in south Jordan, e.g. [1] investigated the crustal structures in Jordan using potential field data [2]. Confirmed the sinistral strike-slip motion on the Dead Sea Rift based on aeromagnetic data surveys.
The principle objective of this study is to delineate and mapped the aeromagnetic lineaments and compare them with the surface lineaments deduced from analysis of Landsat Imagery data and surface geological information, demarcate the subsurface faults and depths of magnetic sources bodies. It can also be used to identify the important trends and structures in the magnetic anomaly field.
Landsat ETM+ imageries were downloaded from USGS Global Visualization Viewer website (http://glovis. usgs.gov/). The downloaded images were acquired in August 2003; PCI Geomatics version 9.1 was used for digital image processing and lineaments extraction from the satellite imageries.
The area of study forms the southern part of Jordan Figure 1 with coordinates 13,000E - 270,000E, −140,000N - −63,000N according to Palestine Grid Coordinate system 1924 (PGS). It covers an area of about 11,000 km2. The elevation of study area ranges from 250 m below mean sea level (bmsl) at the most western part of the study area to more than 1800 above mean sea level (amsl) at the eastern and northern parts of the study area. Figure 2 shows the digital elevation model (DEM) for the study area.
2. Geological and Structural Setting
Miocene to recent northward motion of the Arabian plate, which caused the opening of the Red Sea and Gulf of Aden, has controlled much of the fault pattern of Jordan Figure 3. The following fault systems can be noticed in the study area:
1) Ras En-Naqab Fault System (RNFS), which has a NW-SE trend. It most likely has formed as a consequence of extension related to the opening of the Red Sea during the Miocene. Tensional stress evidenced by the NW-SE basaltic dyke intrusions of Ras En-Naqab Fault System has formed an escarpment extending from Ras En-Naqab area in the northwest to Batn Al-Ghul area in the southeast near the Saudi Arabian border with more than 100 km length [3]. This fault zone includes many elongated grabens and horsts in addition to tilted fault blocks in the same trend of the main fault and with variable widths and lengths Figure 3 ([3,4]).
2) Wadi Araba-Dead Sea fault and associated faults. This major fault is distinctive on the Landsat images Figure 2 and extends outside the study area; with a left
Figure 1. Location map of the study area.
Figure 3. Structural pattern of Jordan (modified [4]).
lateral movement of about 110 km.
3) Al Quweira Fault Zone (QFZ). These faults pass about 4 km west of Petra, (Finan area) and extend for several hundreds of kilometers into Saudi Arabia in the south. These faults are particularly obvious in the AlQuweira area in Jordan Figure 3 [3]. East of the AlQuweira fault zone (QFZ), many N-S faults exhibit sinistral horizontal movement, other faults in south Jordan trend NW-SE and appear to have dip-slip and reverse movements [3]. The displacement along this fault was determined to be 40 km sinistral movement [3]. The study area is characterized by its high diversity in terms of the exposed rock units. Five major types of rock units and formations can be traced; these are Figure 4:
1) The Basement Complex which crops out in the southwestern part of the study area and involves different igneous rock units of Precambrian age. The main rock type is granite (Basement Complex), with varying chemical and mineralogical constituents, such as, granodiorite, basic and acidic intrusive igneous rocks and dolerite. The granite basement is interlaced with many igneous dykes with variable constituents and orientations [5].
2) Disi Sandstone of Paleozoic age, which consist of white to grey cross-bedded sandstone. These rocks were deposited predominantly as large-scale sub-aqueous dunes in high velocity, high discharged braided river with marine incursions [5].
3) Khrayim group of Middle Ordovician-Upper Ordovician age, consisting of silty clay with clay concretion, micaceous fine-grained sandstone with hummocky crossstratification and ripple marks in its lower parts. In the upper part it is composed of highly burrowed micaceous sandstone, thick, massive lensoids of limestone with low angle trough cross-bedding and ripples cross-lamentations inter-bedded with fine-grained sandstone and siltstone. This group covers the most southeastern part of the study area [5].
4) Mesozoic Sandstone, consists mainly of the Lower Cretaceous Kurnub Sandstone Group, this Group rests with a regional angular unconformity on the Ordovician rocks along Ras En-Naqab-Batn Al-Ghoul escarpment. It is characterized by its white-varicolored, fine medium-
Figure 4. Simplified geological map of the study area [4].
grained, moderate to well sorted cross-bedded sandstone with silty clay and some soil horizons. The depositional environment was predominantly fluvial to marginal marine [5].
5) Volcanic rocks represented by NW-SE trending Miocene-Pliocene basaltic dykes found in the southern part of the study area, locally intruded along the QuweiraMudawwera fault system and containing blocks of Upper Cretaceous limestone and chert [5].
3. Aeromagnetic Study
3.1. Data Sources and Methodology
The aeromagnetic data were provided by the Natural Resources Authority (NRA) in Jordan. They are part from a countrywide airborne magnetic survey which was carried out in 1979-1980 by Geoterrex limited [6]. The data were available in the form of a final dataset corrected for the time-magnetic variations [7]. These data were digitized, plotted and interpolated at a grid points using Surfer Version 8.0, and a total aeromagnetic map was prepared Figure 5. High and low pass filtering techniques were applied on the dataset for enhancement and removing undesired signals. The total aeromagnetic intensity datasets were incorporated into ArcGIS 9.3 environment for further subsequent analysis.
3.1.1. Total Aeromagnetic Intensity Map
The total magnetic field, scale 1:500,000 in the study area Figure 5 ranges from 42,700 to 43,570 nT. It shows two large dominant magnetic anomalies. The first is elongated feature with NE-SW orientation located in the western part of the study area, having a strike 015 - 030
degrees, and a length of more 110 km. It was interpreted as a major fault because its boundaries are sharply terminated, and it is coincides with the Quwiera fault. The second large magnetic anomaly is a NW-SE elongated feature, having a strike of about 110 - 130 degrees, and a length is more than 100 Km, and is associated with the Ras En-Naqab escarpment and Mudawwara Fault. The aeromagnetic intensity map clearly reflects the underlying basic geology of the study area. The basement that crop out in the study area produces high frequency, low amplitude magnetic anomalies, which are superimposed upon larger long wavelength anomalies of deeper crustal origin [1].
3.1.2. Reduction to the Magnetic Pole
The total aeromagnetic field map was reduced to the pole (RTP) using the Hansen and Paw bosky method [8], in which the magnetic field is transformed in such a way as if all the sources are magnetized vertically [9], in order to locate the observed magnetic anomalies over the magnetic source bodies causing these anomalies. The pole reduction is a useful aid to interpretation, because asymmetric anomalies become symmetric [10]. The reduction operation assumes that the total magnetic intensity to be 44,222 nT, inclination 47˚, and declination 2.64˚ as deduced from the IGRF model for the corresponding epoch 1980. The resulted reduction to the magnetic pole map Figure 6 is characterized by intense magnetic anomalies constituting several positive and negative peaks. The positive peaks reflect very strong magnetic anomalies with different amplitude, and may be attributed to the occurrence of subsurface intrusions of high magnetic contents with high magnetic susceptibility at different
Figure 6. Reduced to the magnetic pole map.
depths. It shows several dominant trends generally, aligned along definite axis that can be used to define the magnetic provinces. The most dominant trend is the NW-SE and the, NE-SW directions. In addition, several minor small high and low magnetic features were found to trend in E-W direction.
3.1.3. Horizontal Gradient
The enhancement of magnetic anomalies associated with the faults and other structural discontinuities was achieved by the application of band pass filter, and calculating the horizontal derivative of the aeromagnetic datasets. The horizontal gradient (HG) can be defined as follows, if T(x,y) is the magnetic field and the horizontal derivatives of the field in x and y directions are (δT/δx, and δT/δy) respectively, then the Horizontal Gradient HG (x,y) is given by [11]:
(1)
The horizontal gradient (HG) method is considered as the simplest method approach to estimate the contact locations (e.g. faults) [11]. The horizontal gradient map aided in defining the locations of the faults which in turn were related to fault trend and fractures in the area. The HG method was applied to the aeromagnetic dataset at different azimuths using Geosoft Package (V.6.4.2). The results were incorporated into GIS environment and presented in Figures 7 and 8. It displays detailed information about the structures, contacts and the tectonic setting of the study area.
3.1.4. Vertical Gradient
The vertical derivative (gradient) operator is often applied
Figure 7. Horizontal zero degree magnetic field map.
Figure 8. Horizontal 90˚ magnetic field map.
to ground or airborne magnetic data to sharpen the edges of anomalies and other significant features [10]. The algorithm is given by Gunn [12]:
(2)
where A (u,v) is the amplitude present at those frequentcies, and n is the order of the derivative. The higher the order of derivative, the noisier the map will become [10]. The first vertical gradient of aeromagnetic data was calculated using Geosoft package V.6.4.2, incorporated into the GIS Environment, and presented in Figure 9.
3.1.5. Analytical Signal
The analytical signal of the magnetic anomaly field can be effectively used to map the edges of 3-D and 2-D bodies [13]. This method is based on the use of the spatial derivatives of magnetic anomalies (F(x,y)) computed along three orthogonal directions [14]. By this operation the total absolute value of analytical signal
of a given magnetic anomaly field can be obtained by summing the three square differentials [15]:
(3)
The analytical signal displays the maxima over 2-D source body edges even when the direction of magnetization is not vertical [15]. In this study, the analytical signal of the aeromagnetic data was calculated using Geosoft package V.6.4.2 and the result is incorporated into GIS environment and presented in Figure 10. The
Figure 9. Vertical gradient magnetic field map.
Figure 10. Analytical signal magnetic field map.
structural features, aeromagnetic lineaments, and their different alignments in the study area were deduced from the different geophysical maps Figure 11.
3.1.6. Euler Deconvolution
Euler deconvolution [16] is a method to estimate the depth of subsurface magnetic anomalies and can be applied to any homogeneous field, such as the analytical signal of magnetic data [17]. It is particularly good at delineating the subsurface contacts. It has been observed that the depth estimates from magnetic data are more
Figure 11. Aeromagnetic lineaments map and location of cross section profiles.
accurate if the pole-reduced magnetic field is used [18], than using the magnetic data themselves.
In this study, the reduced to the magnetic pole map was used to map the subsurface structures and contacts using EULDEP software [17]. This software is based on calculating the horizontal and vertical gradients of the magnetic field (RTP in this case) along profile data [17]. The method is an automatic depth estimates one, and is designed to provide computer assisted-analysis on large volumes of magnetic data [19]. In addition, to making use of the property of the rapid decay of intensity of magnetization with distance from the source. The intensity of magnetization or magnetic field
as a result of the distribution of magnetic poles can be written as follow:
(4)
where r is
, M is proportional to magnetization and N is the structural index with values ranging from 0 to 3 [17].
In this Study, different aeromagnetic profiles with different orientations Figure 11 distributed over the study area and crossing the major magnetic anomalies have been constructed processed with different structural indices Table 1 and interpreted in terms of sub-surface structures Figure 12.
The dataset for each profile has been digitized with the aid of GIS techniques. The profile data were then processed using the EULEDEP software. In addition, the surface fault traces from surface detailed geological maps were mapped and delineated using GIS techniques along each profile and were incorporated along with the results
Table 1. Symbols used to plot structural indices (S.I) and subsurface faults shown in (Figure 11).
of the processed profile data (Figure 12, Red and Green dashed lines). The results show that, several subsurface faults have been inferred from the modeled aeromagnetic profiles data which were not documented in surface geological maps (Black dashed lines, Figure 12). Some of the delineated geological faults were not detected by modeled aeromagnetic profile data (Red dashed lines, Figure 12). But, most of the delineated geological faults were found coincident with the interpreted aeromagnetic faults (Green dashed lines, Figure 12(e) and Figure 12). Moreover, the modeled aeromagnetic profiles data showed that the inferred faults may extend downward for several kilometers.
3D Euler solution for the aeromagnetic data of the study area, at Structural Index (S.I = 1) have been produced using the GEOSOFT packages and the result is presented in Figure 13.
3.2. Landsat Interpretation
Four scenes of Landsat ETM+ images were downloaded from USGS website that covers the study area. These images were converted into PIX format to become suitable for PCI Geomatics 9.1 software. Mosaic and subset operations were performed to clip the images for the study area Figure 14.
Available structural maps were also used in this study. However, automated lineament extraction techniques had to be used to increase the details of existing data of structural maps. The main advantages of automated lineament extraction over the manual lineament extraction are; the ability to uniformly approach the different images; processing operations are performed in a short
Figure 12. Sub-surface structural interpretation and magnetic causative anomalies.
Figure 13. 3D Euler solution from the analysis the aeromagnetic data over the study area, with structural index (S.I = 1) combined with geological and inferred aeromagnetic structures.
Figure 14. Land sat ETM+ imagery for the study area.
time, and the ability to extract lineaments which are not recognized by human eyes. This was done by applying LINE module of PCI Geomatica 9.1 software. This module extracts linear features from an image and records the polylines in vector segments by using six parameters. These parameters are (PCI Geomatics 9.1 Users Manual):
1) RADI (Filter radius): This parameter is used in the first step of the first stage of the process for the “Canny edge detection”. It specifies the radius of the edge detection filter (in pixels), and roughly determines the smallest-detail level in the input image to be detected.
2) GTHR (Gradient threshold): This parameter is used in the second stage of the process for the “Canny edge detection”. It specifies the threshold for the minimum gradient level for an edge pixel to obtain a binary image.
3) LTHR (Length threshold): This parameter is used in the third stage of the process. It specifies the minimum length of curve (in pixels) to be considered as lineament or for further consideration (e.g., linking with other curves).
4) FTHR (Line fitting error threshold): This parameter is used in the second step of the third stage. It specifies the maximum error (in pixels) allowed in fitting a polyline to a pixel curve. Low FTHR values give better fitting but also shorter segments in polyline.
5) ATHR (Angular difference threshold): This parameter is used in the last step of the third stage of the process. It specifies the maximum angle (in degrees) between segments of a polyline. Otherwise, it is segmented into two or more vectors. It is also the maximum angle between two vectors for them to be linked.
6) DTHR (Linking distance threshold): This parameter is used in the last step of the third stage of the process. It specifies the minimum distance (in pixels) between the end points of two vectors for them to be linked. The used parameters and their values appear in Table 2.
Figure 15 shows the results obtained by applying the LINE module.