Interpretation of magnetic data in the Sinop area of Mid Black Sea

Interpretation of magnetic data in the Sinop area of Mid Black Sea,Turkey,using tilt derivative,Euler deconvolution,and discrete wavelet transformB.Oruça ,⁎,H.H.Selim ba Kocaeli University,Engineering Faculty,Department of Geophysical Engineering,Umuttepe Campus,41380İzmit/Kocaeli,Turkey bKocaeli University,Engineering Faculty,Department of Geological Engineering,Umuttepe Campus,41380İzmit/Kocaeli,Turkeya b s t r a c ta r t i c l e i n f o Article history:Received 8March 2011Accepted 25May 2011Available online 12June 2011Keywords:Magnetic anomalies Tilt derivativeEuler deconvolutionDiscrete wavelet transform Structural lineamentsThe identi fication of lineaments in magnetic anomaly data is an important step in interpretation.In this paper,three methods for locating magnetic lineaments are applied to total field magnetic anomaly data of the Sinop area of mid Black Sea,northern Turkey.Based upon the variation in magnetic lineaments,structural map of the study area was obtained.The tilt derivative (TDR),Euler deconvolution (ED)and discrete wavelet transform (DWT)have been of great utility in the interpretation of magnetic anomaly data.These three methods were applied to enhance magnetic anomalies due to anomalous sources in part of the basement complex of study area.The TDR obtained from the magnetic anomaly reduced to pole (RTP)was designed to map geologic contacts.The source boundaries and depths are determined from the zero contours,and the half distance between ±π/4contours or the distance between zero and +π/4or −π/4contour of TDR,respectively.The ED technique is used to estimate the source depth at the semi-in finite contact location.The DWT technique was carried out to infer the substructure orientation.The direction of maxima and minima coef ficients of the horizontal,vertical and diagonal decompositions of the DWT works in imaging the horizontal locations of the anomalous sources.The results allow us to image and characterize subsurface lineaments for the survey area.©2011Elsevier B.V.All rights reserved.Contents 1.Introduction ..............................................................1942.Geological setting............................................................1953.Magnetic data .............................................................1953.1.Data enhancement........................................................1953.1.1.TDR image .......................................................1963.1.2.Euler deconvolution solutions ..............................................1973.1.3.Two-dimensional DWT parison of ED with TDR andDWT images.......................................2003.2.Delineation of magnetic features .................................................2014.Conclusions ..............................................................202Acknowledgments ..............................................................203References ........................... (203)1.IntroductionAn important objective in the interpretation of potential field data is to improve the resolution of observed data.In magnetic prospecting,inareas of limited exposure,delineating lateral change in magnetic susceptibilities provides information not only on lithological changes but also on structural trends.Especially,mapping the edges of causative bodies is fundamental to the application of potential field data to geological mapping.The edge detection techniques are used to distinguish between different sizes and different depths of the geological discontinuous.The derivatives of magnetic data are used to enhance the edges of anomalies and improve signi ficantly the visibilityJournal of Applied Geophysics 74(2011)194–204⁎Corresponding author.Fax:+902623352812.E-mail addresses:bulent.oruc@.tr (B.Oruç),haluk.selim@.tr (H.H.Selim).0926-9851/$–see front matter ©2011Elsevier B.V.All rights reserved.doi:10.1016/j.jappgeo.2011.05.007Contents lists available at ScienceDirectJournal of Applied Geophysicsj o u r n a l h o me p a g e :w w w.e l sev i e r.c om /l oc a te /j a p p ge oof such features.Many techniques for mapping have been developed to delineate structural features from magnetic data.Cordell and Grauch (1985)suggested a method for the location of the horizontal extents of the sources from the maxima of the horizontal gradient of the pseudogravity computed from the magnetic anomalies.Finding the maxima is then efficiently done with the curve-fitting approach of Blakely and Simpson(1986).The wavelet transform has been introduced into the geophysical field(Foufoula-Georgiou and Kumar,1994)as a powerful analysis technique.In contrast to the large number of studies into wavelet methods for seismic applications,there appears to be a relatively small amount of published research into wavelet methods for potentialfields.Moreau et al.(1997)applied the DWT analysis to the case of gravity or magnetic anomalies,using wavelets built on the Poisson kernel.Fedi and Quarta(1998)used the wavelet transform for the separation of regional and residualfields.Hornby et al.(1999) have used the wavelet transform to the analysis of gravity data for the characterization of geological boundaries.Ridsdill-Smith and Dentith (1999)has separated aeromagnetic anomalies using wavelet matched filters.Moreau et al.(1999)have developed an efficient technique for identification of sources of potentialfields.The DWT can be employed to transform the potentialfield data in the space scale domain(Fedi and Quarta,1998;Moreau et al.,1997,1999)without any a priori information.Li and Oldenburg(2003)have demonstrated the application of the wavelet domain to be inversion of magnetic data. Albora et al.(2004)have applied the DWT to estimate the boundaries of an archeological site.In the last few years there have been a number of methods proposed to help normalize the signatures in images of magnetic data so that weak,small amplitude anomalies can be amplified relative to stronger,larger amplitude anomalies.Verduzco et al.(2004)have developed tilt derivative from gravity or magneticfield anomaly map and suggested using the HGM of the tilt derivative as an edge detector. Keating and Pilkington(2004)have showed that the HGM of the tilt derivative is related to the local frequency.Wijns et al.(2005)use the analytic signal amplitude to normalize the total horizontal derivative. Cooper and Cowan(2006)have used tilt angles which act as normalized horizontal derivatives to detect the edges of anomalous sources.Salem et al.(2008)have developed a new method for interpretation of gridded magnetic data based on the Tilt derivative, without specifying prior information about the nature of the source.The main goal of this work is to analyze magnetic anomaly data to identify possible lineaments in the Sinop basin area,mid Black Sea, Turkey.Our pattern recognition criteria are based on the TDR,ED,and DWT decompositions obtained from RTP anomaly data.All three methods agree closely in determining the horizontal locations of contacts;thefirst two also give similar results for source depths. 2.Geological settingThe Black Sea is a large marginal sea located within the Alpine orogenic belt,surrounded by compressive tectonic belts,the Pontides orogeny in the south,Caucasus in the northeast and Crimean Range in the north(Çiftçi et al.,2003).It is located on the westernflank of the active Arabia–Eurasia collision and north of the North Anatolian Fault (NAF)that permits the tectonic escape of Anatolia(Rangin et al., 2002).The region was under the effect of compressive forces in the direction of North Northeastern–South Southwestern between Upper Cretaceous–Upper Miocene and at the end of the Miocene,this activity was suppressed by the activity of North Anatolian Fault(NAF) system(Özhan,2004).As shown in Fig.1a,the Black Sea basin comprises western and eastern Black Sea subbasins,which are separated by a regional high,the mid Black Sea Ridge,which is divided into two parts,the Andrusov Ridge in the north and the Archangelsky Ridge in the south(Dondurur andÇiftçi,2008).Gökaşan(1996),Özhan(2004),and Yılmaz et al.(1997)have reported the geological evolution of Black Sea in land and sea areas. The Sinop Basin is located between the Archangelsky Ridge and the Turkish coastline(Fig.1a)and has been affected by late Miocene normal faults along the Turkish margin and the Archangelsky Ridge (Rangin et al.,2002).The overall compressional deformation of the Sinop Basin was superimposed onto the footwall block forming the Archangelsky Ridge(Çiftçi et al.,2003).In Sinop peninsula,Miocene aged sandy limestones overly Upper Crateceous agglomerates (Fig.1b).3.Magnetic dataThe total intensity magneticfield dataset provided by the Turkish Navy Department of Navigation,Hydrography and Oceanography were recorded using the cesium-vapor marine magnetometer with sensitivity at0.004nT.The acquisition of magnetic data in the southern part of the Sinop basin was carried out along a north–south-striking profile perpendicular to the main tectonic structures. The magnetic data was acquired to allow structural interpretations to be carried out in regions where little or any data had not previously existed.Fig.2a shows the totalfield anomaly map covering a portion of the Sinop basin(onshore partially).In the study area,the magnetic signature is generally smooth,with anomalies of different wave-lengths and amplitudes,denoting the different magnetic contents of the Sinop basin infill.A correlation is seen between linear and circular patterns on the anomaly map.The anomaly map shows high amplitude and high wavenumber anomalies in the region.The high gradient contours suggest fault like structures in the central part of the study area in WNW–ESE direction(Fig.2a).This major trend was observed along the border of prominent high-amplitude linear anomalies with contrasting magnetic properties.3.1.Data enhancementThe totalfield anomaly map was processed to prepare for analysis and interpretation.The reason for this pre-processing is that the major features in the data which may be defined by the edges and minor features that are not easily seen in the original image are to detect.Before applying the TDR,ED and DWT methods,the totalfield anomaly data were converted to RTP using a magnetic inclination of 60°and a declination of5.5°.Thus,the RTP reduces the effect of the Earth's ambient magneticfield and provides a more accurate determination of the position of anomalous sources.It is therefore understood that the total magnetization direction is equivalent to that of the current inducingfiled.As is well known,this assumption states that the direction of magnetization in the anomalous sources lies in the same direction as the local magneticfield.Unfortunately,we had to ignore the remanent magnetization because of a lack of relevant information.Thus the real data example demonstrated in this section could correspond to the case of zero remanent magnetization. However,in the case of the presence of strong remanent magnetiza-tion,this will adversely affect the interpretation of the magnetic anomaly data and lead to erroneous size or shape.The RTP signature of the Sinop Basin is marked by a rugged relief with positive and negative anomalies of several wavelengths and amplitudes,with values ranging between−900and1400nT (Fig.2b).Northward there is a low NE–SW gradient which appears on both magnetic anomaly maps(Fig.2a and b).The RTP anomaly,in the central portion of the study area,shows an important gradient zone elongated in approximate E–W direction,characterized by much longer wavelength,most likely reflecting variations in lateral basement magnetization contrasts.In the south of study area,other important positive anomalies in the NE–SW and E–W directions are caused probably by the volcanic structures.195B.Oruç,H.H.Selim/Journal of Applied Geophysics74(2011)194–2043.1.1.TDR imageThe magnetic images are used to indicate areas of considerable magnetic contrast and to visualize features such as faults and dykes, which are depicted as lineaments.Because the amplitude of magnetic anomalies depends on magneticfield strength and depth of anomalous sources,lower amplitude,anomalies may be suppressed at the expense of higher amplitudes.For this reason,the edge-detection filters are presented for delineating linear features without diminish-ing the long-wavelength information.One of the most important filters is TDR,an edge-detectionfilter.The tilt derivative or tilt angle wasfirst proposed by Miller and Singh(1994)as a tool for locating magnetic sources on magnetic profile data.The horizontal gradient magnitude(HGM)is given by the square root of the sum of the squares of the horizontal derivatives of the potentialfield f:HGM=ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi∂f∂x2+∂f∂y2s:ð1ÞFig.1.(a)Simplified tectonic map showing major tectonic elements of the Black Sea(after Robinson et al.,1996;Kazmin et al.,2000).(b)Geological sketch map of the Sinop area (Varol,2004)and the location of the study area.196 B.Oruç,H.H.Selim/Journal of Applied Geophysics74(2011)194–204where ∂f /∂x and ∂f /∂y are the first derivatives of the field f in the x and y directions.Verduzco et al.(2004)generalized the concept so that it could be applied to grid data by de fining a generalized tilt derivative asTDR =tan −1∂f ∂z HGM 0B B @1C C A ;ð2Þwhere ∂f /∂z is the first vertical derivative of the field f .The TDR values are restricted to values between −π/2and +π/ler and Singh (1994)have shown that the TDR crosses through zero at or near the edge of a vertical-sided source and is negative outside the source region.The TDR is therefore very effective in allowing anomalies to be traced out along strike.Verduzco et al.(2004)have shown that this derivative also performs an automatic-gain-control (AGC)filter which tends to equalize the response from both weak and strong anomalies.Salem et al.(2007)have shown that half-distance between ±π/4contours provides an estimate of the source depth for vertical contacts or the distance between zero and +π/4or −π/4contour obtained from the TDR corresponds to the depth to the top of the vertical contact model.Fig.3shows that the TDR map facilitates the recognition of the horizontal location and extent of edges of anomalous sources assuming vertical contact model.The zero contours estimate the horizontal location of abrupt lateral changes in susceptibility.The TDR map accentuates short wavelength and reveals the presence of magnetic lineaments.As shown in Fig.3,the structural elements are enhancedsimply by inspecting the alignment of anomalies and by distinguishing the noticeable abrupt change between positive and negative anomalies in Fig.2b,particularly at the location of sharp gradients.It is interesting to note that the zero contours of TDR images produce the fully elongated zone by nonclosed contours in an approximate E –W trending.Fig.3shows also the high frequency images of TDR from RTP anomaly data.This area is ornamented with zero contours of TDR,are characterized by four zones.The half distance between ±π/4contours is used to estimate the depth of the edge of the magnetized sources.Accordingly,the TDR indicated the anomalous sources began at a depth of 0.2km and extended to 1km with an average depth estimated at 0.6km (Fig.3).Thus,one can easily trace the lineaments corresponding to a zone along the southern edge of Sinop basin.In this gradient zone,the east part of basin edge is deepest edge and is obtained shallow depths towards the west.3.1.2.Euler deconvolution solutionsIn order to estimate the positions of structural lineaments in the study area from total field data,Euler deconvolution method was employed.The reason for this is that the method requires no prior information about the source magnetization direction,and its results are not affected by the presence of remanence (Ravat,1996).The method primarily responds to the gradients in the data and effectively traces the edges and de fines the depths of the source bodies.The 3D Euler equation is written as (Reid et al.,1990;Thompson,1982)x −x 0ðÞ∂f ∂x +y −y 0ðÞ∂f ∂y +z −z 0ðÞ∂f ∂z=N B −f ðÞð3Þwhere x ,y and z are the coordinates of a point of observation,x 0,y 0,and z 0are the coordinates of the source location,and B is a base level.The structural index (SI)N ,de fines the anomaly attenuation rate at the observation point and depends on the geometry of the source.The advantages of the derivative mode of the Euler method are actually derived from the general advantages of using the derivatives in the potential field applications.The combination,using the Euler method on derivatives,is a powerful tool to characterize shallow sources.As long as the signal-to-noise ratio is suf ficiently high,the method can be applied to higher-order derivatives to locate shallow sources.Accordingly,Hsu (2002)gave the general formula for Euler's equation as∂∂n f n x −x 0ðÞ+∂∂n f n y −y 0ðÞ+∂∂n f n z −z 0ðÞ=−N ∂n fn ;ð4ÞFig.3.TDR image obtained from Fig.2b.Dashed lines show the 0radian contour of the TDR.Solid lines are contours of TDR for +π/4and −π/4radians.Fig.2.(a)Total field anomaly map of survey area,(b)RTP anomaly data computed fromtotal field anomaly.197B.Oruç,H.H.Selim /Journal of Applied Geophysics 74(2011)194–204where n is the order of the gradient used.When the ED can be applied to the first vertical gradient (n =1)of the magnetic field (T =f ),Eq.(4)provides a solution:x −x 0ðÞ∂2T ∂x ∂z +y −y 0ðÞ∂2T ∂y ∂z +z −z 0ðÞ∂2T ∂z2=−N ∂T∂z :ð5ÞThe source points that are calculated as solutions by ED arepositioned at the estimated border of the susceptibility inhomogene-ities.Thus,the ED relies on the derivatives of the magnetic data,the resulting depth estimates relate primarily to the areas of basement heterogeneities identi fied as distinct sources of the field.The first vertical gradient of magnetic data is calculated by using the fast Fourier transform (FFT)technique (Gunn,1975).The horizontal derivatives and vertical derivative of the first vertical gradient,necessary for the calculation of Eq.(5),are also been calculated using the FFT technique.The horizontal source locations from ED solutions can be used for delineation of structural and lithological trends.A location where these solutions tend to cluster is considered to be the most likely location of the source.In our study,we seek to locate the edges of anomalous sources assuming magnetic contacts which may delineate sedimentary basin.As s well known,two important parameters of Euler deconvolution are the choice of window and N .The N value must be assumed as a priori information.Reid et al.(1990)and Thompson (1982)showed that the optimum structural index usually yields the tightest clustering of the solutions.For the first derivative field of magnetic data,the N values are 1for a contact,2for a semi-in finite thin dike (Stavrev,1997).However,Schmidt (2006)has showed the N value could also be solved and accuracy of the ED to the magnetic gradient tensor is improved.These are theoretical model-derived estimates of field fall-off rates and should be integer values.In practice,the Euler solutions provide non-integer values of N ,which are rounded to the nearest integer values for subsequent analysis since real geological bodies are more complex than simple model bodies,and thus the application of ED to real data is always an approximation.Ravat (1996)discussed the effect of the size of the moving window to estimate the source location using the ED.Generally,the choice of the window size should be selected to be large enough to incorporate substantial variation of the total field anomaly and their gradients (Ravat,1996).Therefore,the window size has to be adapted to the structural dimensions of the anomalous sources to obtain optimal results.We have adopted a window size of 5×5with a 0.78km grid cell size in the ED to the size of the observed features to mostly encompass the effects of edges of causative bodies.Fig.4a shows the good clustering solutions computed using N =1.2with depth uncertainty set to a maximum of 20%using the first vertical derivative field.The N value indicates that anomalous sources are close to contact-like structures.Thus the causative bodies are mainly contacts,but also the contacts itself show a slight contrast and N is slightly shifted towards the thin dyke value.Euler's depth estimates imaged well the source edges in central,NW and SE sectors of the study area (Fig.4).Note that the ED that yielded best results to resolve source depths,that mapped the source edges are as shown in Fig.3.The depth values from the ED are positioned at the boundaries of contrasting susceptibility and range from 0.27to 0.7km.3.1.3.Two-dimensional DWTIn this study,we use orthogonal wavelet transform,since it allows an input data to be decomposed into a set of independent coef ficients,corresponding to each orthogonal basis.The basic idea of wavelet transform (WT)is to analyze different frequencies of a signal using different scales.High frequencies are analyzed using low scales while low frequencies are analyzed in high scales and thus enable us to analyze both local and global features.In wavelet transform,all of thebasis functions,which are called wavelets,are derived from scaling and translation of a single function,called mother wavelet.Many types of mother wavelets and associated wavelets exist (Daubechies,1992).Mallat (1989)introduced an ef ficient algorithm known as DWT.The DWT of a signal f(x)may be achieved in different scales of the frequency domain by means of an orthogonal family of wavelet functions.Each wavelet,ψa ,b ,is de fined by the scaling and shifting of the mother wavelet ψas follows:ψa ;b x ðÞ=1ffiffiffiap ψx −bð6Þwhere a and b are scaling (or level)and shifting (or location)parameters,respectively.For certain mother functions,the set of wavelets ψa ,b forms a smooth and orthogonal basis.The wavelet transform is used to de fine the wavelets de fined in Eq.(6)as basis functions:w a ;b ðÞ=1ffiffiffia p ∫∞−∞f x ðÞψx −ba dx :ð7ÞEq.(7)is also called wavelet decomposition of the f (x )by means ofthe set of wavelets ψa ,b .Therefore,the set of all wavelet coef ficients w (a ,b )gives the wavelet domain representation of the f (x ).If the parameters a ,b in Eq.(7)take discrete values,this transform is called a discrete wavelet transform (DWT),otherwise it is called a continuous wavelet transform (CWT).The two-dimensional DWT can be implemented by performing the one-dimensional splitting algorithm to the horizontal (rows)and vertical lines (columns)of an image,successively.However,one-dimensional WT is performed sequentially level by level.After the first level wavelet transform operation,the input image can be divided into four:approximation,horizontal detail,vertical detail and diagonal detail where the size of each part is reduced by the factor of two compared to the input image as depicted by Fig.5a.The horizontal detail coef ficients,vertical detail coef ficients and diagonal detail coef ficients contain the horizontal,vertical and diagonal components of the input data.When the second level WT is applied,the approximation part of the first level will be furtherdecomposedFig.4.(a)Euler deconvolution solutions obtained from the first vertical derivative calculated by Fig.2b.The solutions are superimposed with first vertical derivative map.Euler's depth estimates are classi fied in two different classes represented by circles based on their depths.(b)Comparison of tilt images and ED results.The dashed lines represent the zero contours of the TDR.198 B.Oruç,H.H.Selim /Journal of Applied Geophysics 74(2011)194–204into four components as shown in Fig.4b.For the higher level,the iteration is done in the same way until the desired level is reached.Fig.4c shows the schematic representation of this process at level j +1.Accordingly,the two-dimensional DWT is illustrated from level j to level j +1,where the wavelet representation at level j +1consists of four sub-images.The first sub-image cA j +1is obtained by applying the horizontal low-pass filter and the vertical low-pass filter successively (approximation detail coef ficients).The second sub-image cD j +1h (horizontal detail coef ficients)is obtained by applying the horizontal low-pass filter followed by the vertical high-pass filter.The third sub-image cD j +1v (vertical detail coef ficients)is obtained by applying the horizontal high-pass filter followed by the vertical low-pass filter.Finally,the fourth sub-image cD j +1d (diagonal detail coef ficients)is obtained by applying the horizontal and vertical high-pass filters successively.The same process can be applied to every wavelet level.3.1.3.1.Theoretical example.In this section,the robustness of the DWT technique used for the edge enhancement is tested with wavelet decomposition of total field anomaly map (Fig.6)caused by two vertical-sided prisms at a depth to the top of 1.5km (labeled 1)and 2km (labeled 2).Fig.6shows the theoretical total field anomaly map due to the two vertical sided prisms in the case of induced magnetization with inclination and declination of 90°and 0°,respectively.The theoretical magnetic anomalies of the vertical sided prisms are calculated using the formula given by Rao and Babu (1993)on a regular grid with a spacing of 0.1km.Fig.7a and b shows that the DWT leads to a decomposition of approximation coef ficients at levels 1and 2in four components:the approximation coef ficients (A 1and A 2),and the details in three orientations (horizontal —H 1and H 2,vertical —V 1and V 2,and diagonal —D 1and D 2),respectively.It should be noted that A 1and A 2coef ficients well re flect the theoretical anomaly map in Fig.6.Thus,one can conclude that theoretical anomaly was successfully decomposed in three orientations.The resolution of the coef ficients (H 1,H 2,V 1,V 2,D 1and D 2)in three orientations is similar in that their maxima and minimawork well in imaging the edges of the prisms (Fig.7a and b).It is interesting to note that the maxima and minima of H 1,H 2,V 1and V 2locate the positions of all the edges of the tilted prism while imaging N –S and E –W edges of the other one.The D 1and D 2have imaged the corners of the prism models,as expected.As shown in Fig.7a and b,the advantage of wavelet transform is that it preserves most of the information of the original image in the coarse image.Consequently,we concluded that image decomposition process keeps maxima and minima coef ficients of the decomposed image in the first level and second level of decompo-sition using the Haarwavelets.Fig.5.The concept of wavelet transforms.(a)The first level wavelet transform.(b)The second level wavelet transform.(c)Two-dimensional wavelet transform leads to a decomposition of approximation coef ficients at level j in four components:the approximation at level j +1,and the details in three orientations (horizontal,vertical,and diagonal).The c and A ,D represent the wavelet coef ficients,approximation coef ficients and detail coef ficients,respectively.The arrows represent sub-sampling,a decimation of the coef ficients by a factor of2.Fig.6.Total field magnetic anomaly of vertical prism models at the depths of 1.5km (prism 1)and 2km (prism 2)for magnetization vector at 90°and magnetization strength of 470A/m.199B.Oruç,H.H.Selim /Journal of Applied Geophysics 74(2011)194–204。