袁園, 黃大年, 余青露
1 吉林大學地球探測科學與技術學院, 長春 130026 2 國家海洋局第二海洋研究所, 杭州 310012 3 國家海洋局海底科學重點實驗室, 杭州 310012 4 中國石化石油物探技術研究院, 南京 211103
?
利用加強水平方向總水平導數(shù)對位場全張量數(shù)據(jù)進行邊界識別
袁園1,2,3, 黃大年1*, 余青露4
1 吉林大學地球探測科學與技術學院, 長春 130026 2 國家海洋局第二海洋研究所, 杭州 310012 3 國家海洋局海底科學重點實驗室, 杭州 310012 4 中國石化石油物探技術研究院, 南京 211103
位場全張量梯度數(shù)據(jù)以其信息量大、含有更高頻的信號成分,能更好地描述小的異常特征等優(yōu)點在地球物理領域中得到廣泛應用.邊界檢測是位場解釋中不可缺少的任務,需要新的邊界探測器來處理位場梯度張量數(shù)據(jù).為了充分利用位場梯度張量數(shù)據(jù)的多信息成分,本文定義了方向總水平導數(shù)和加強方向總水平導數(shù),并利用其定義新的邊界檢測器.為了能同時顯示不同振幅大小異常的邊界,本文對其進行了歸一化處理.通過模型試驗,證明了歸一化方法能更加清晰準確地顯示淺部和深部的地質(zhì)體邊界信息.最后將該邊界檢測方法用于加拿大圣喬治灣實際測得全張量重力梯度數(shù)據(jù)和中國朱日和地區(qū)的磁異常數(shù)據(jù)中,并得到了較好的邊界檢測結(jié)果.
邊界檢測; 位場梯度張量; 加強方向總水平導數(shù); 歸一化
邊界檢測是位場數(shù)據(jù)解釋中的重要組成部分.目前,有許多方法用于邊界檢測和邊界加強.水平導數(shù)和垂向?qū)?shù)經(jīng)常用于加強這些潛在特征.垂向?qū)?shù)通常被用于描述密度體或磁性體的邊界(Evjen, 1936).后來人們發(fā)現(xiàn)總水平導數(shù)和解析信號振幅的最大值直接描述密度和磁化強度邊界的檢測(Cordell and Grauch, 1985; Cordell, 1979; Nabighian, 1984).Hsu等(1996)提出加強解析信號方法進行到高階導數(shù)以增加該方法的求解能力.Cooper和Cowan (2004)使用不同階的垂向?qū)?shù)進行增強地球物理異常的細節(jié).
Miller和Singh(1994)使用垂向?qū)?shù)與總水平導數(shù)的比來進行邊界檢測.它是第一個均衡濾波器,能均衡不同振幅強度的異常信息.Rajagopalan和Milligan (1995)使用自動增益控制技術均衡不同大小的振幅.然而,這種方法與窗口大小的選擇相關.Verduzco等(2004)用傾斜角的總水平導數(shù)(THDR)來進行邊界檢測.Wijns等(2005)提出了Theta圖濾波器進行邊界探測,它是用解析信號振幅來對總水平導數(shù)進行歸一化處理.這些均衡邊界探測器雖然能同時顯示不同振幅大小的異常的邊界信息,但識別出來的邊界分辨率較低.馬國慶等(2012)利用高階導數(shù)來定義均衡濾波器,提高了對地質(zhì)體邊界的分辨率.但是,當測量的數(shù)據(jù)中同時包含正負異常體的情況下,這些均衡探測器可能會引入錯誤的邊界信息.
近十幾年,位場梯度張量數(shù)據(jù)被廣泛使用.位場梯度張量數(shù)據(jù)可以通過測量或者數(shù)值計算得到,它提供了9個張量成分,包含了比傳統(tǒng)位場數(shù)據(jù)更高頻的信號成分,因此,對于梯度張量數(shù)據(jù)的解釋能得到更高分辨率和更加詳細的地質(zhì)構(gòu)造.這要求新的方法來解釋梯度張量數(shù)據(jù),尤其是對小地質(zhì)構(gòu)造的邊界檢測.Cooper 和Cowan (2006)以及Oru?和Keskinsezer(2008)提出了方向傾斜角來探測邊界,但是只有垂直方向的傾斜角,即傳統(tǒng)的傾斜角方法,能得到較好的檢測結(jié)果.Beiki(2010)和Mikhailov等(2007)定義了方向解析信號來處理重力梯度張量數(shù)據(jù).Beiki(2010)利用方向解析信號進行邊界檢測,比傳統(tǒng)的解析信號方法能得到更好的結(jié)果,但是該方法也不能同時描述大振幅和小振幅異常體.近年來,一些學者利用重力梯度張量數(shù)據(jù)的特征值來進行邊界識別(Oru? et al., 2013; Sertcelik and Kafadar, 2012; Zhou et al., 2013).
Marson和Klingele(1993)指出總水平導數(shù)比解析信號有更高的分辨率.因此,本文定義了方向總水平導數(shù)和加強方向總水平導數(shù),并根據(jù)其定義不同的邊界探測器.為了同時顯示不同振幅大小異常的邊界,給出了歸一化方法.
位場梯度張量數(shù)據(jù)為位場U在空間上的二階偏導數(shù),張量矩陣的形式:
(1)
Beiki(2010)和Mikhailov等(2007)定義了方向解析信號,其表達式為
(2)
(3)
(4)
其中,下標x,y和z表示方向.
這里,我們在x,y和z方向上定義方向總水平導數(shù),其表達式分別為
(5)
(6)
(7)
其中,下標x,y和z表示方向.THDx的最大值描繪地質(zhì)體的N-S向邊界;THDy的最大值描繪地質(zhì)體的E-W向邊界;垂直方向上的總水平導數(shù)THDz與傳統(tǒng)使用的總水平導數(shù)相同,振幅的最大值描繪地質(zhì)體的邊界.
圖1比較分析了方向解析信號振幅和方向總水平導數(shù)的邊界識別能力.圖1a為一垂直棱柱體產(chǎn)生的重力異常.圖1b—1d為x,y和z方向上的方向解析信號振幅;圖1e—1g為x,y和z方向上的方向總水平導數(shù).通過比較可知,方向總水平導數(shù)相對于方向解析信號振幅對于邊界識別有著更高的分辨率.因此,根據(jù)水平方向的總水平導數(shù)來定義邊界探測器有著更高的分辨率.利用THDx和THDy定義地質(zhì)體邊界探測器:
(8)
EDT的最大值定義地質(zhì)體的邊界.
Hsu等(1996)應用高階導數(shù)定義加強解析信號振幅來提高邊界的分辨率.為了提高邊界檢測的分辨率,我們使用高階導數(shù)來定義方向總水平導數(shù),叫做加強方向總水平導數(shù).
x,y和z方向上的加強方向總水平導數(shù)分別為
(9)
(10)
(11)
其中,下標x,y和z表示方向.上標n表示階數(shù),大于等于1.當n=1時,與方向總水平導數(shù)相同.
同樣,加強方向總水平導數(shù)ETHDx的最大值描繪地質(zhì)體的N-S向邊界;ETHDy的最大值描繪地質(zhì)體的E-W向邊界.因此,我們可以根據(jù)加強方向總水平導數(shù)來定義地質(zhì)體邊界探測器:
(12)
EEDTn的最大值定義地質(zhì)體的邊界.通常情況下,我們只取n=1,2.因為位場數(shù)據(jù)高階導數(shù)的計算將放大噪聲的影響.為了識別出來的邊界更加準確,我們應避免高階導數(shù)的計算.圖1h—1j為x,y和z方向上的二階加強方向總水平導數(shù).可以看出,相對于方向總水平導數(shù)有著更高的分辨率.這里,我們用EEDT來表示由二階加強方向總水平導數(shù)定義的邊界探測器EEDT2.在EEDT中,需要計算二階加強方向總水平的導數(shù),然而,垂向?qū)?shù)的直接計算會增加噪聲的影響.我們用如下表達式來計算二階ETHDx和ETHDy導數(shù):
(13)
(14)
(15)
邊界探測器EDT和EEDT都為非均衡濾波器,不能很好地同時顯示不同振幅大小的異常的邊界.為了避免這一缺點,對其采用了歸一化方法.
本文采用位場數(shù)據(jù)的一階垂向?qū)?shù)來對EDT進行歸一化處理,表達式為
(16)
這里引入常數(shù)p,一般大小在0.001到0.1.常數(shù)p的引入是為了避免產(chǎn)生額外錯誤的邊界.大的p值將降低均衡效果,小的p值將加強弱振幅異常的邊界.
邊界探測器EEDT的歸一化表達式為NEEDT=
(17)
同樣,參數(shù)p與(13)式相同.這里,用如下表達式來計算Uzzz:
(18)
為了驗證本文方法的可行性,我們選用兩種常用的方法進行對比.它們是傾角總水平導數(shù)(THDR)和Theta圖(Verduzcoetal., 2004;Wijnsetal., 2005).
建立如下重力異常模型一:模型中包含4個垂直棱柱體,各個棱柱體的上頂埋深分別為1km,1km,2km,2km.圖2顯示了模型體的平面圖和3維圖.各個棱柱體的厚度都為1km,剩余密度為200kg·m-3.圖3為合成的重力異常模型.重力異常模型的梯度張量數(shù)據(jù)可以通過數(shù)值計算和實際測量得到.這里,合成模型的重力梯度張量值通過理論公式計算得到(Forsberg, 1984)(圖4a—4c).圖4d—4i分別為EDT,EEDT,NEDT,NEEDT,THDR和Theta圖的邊界結(jié)果.可以看出,EDT和EEDT為非均衡濾波器,對于深部異常體的邊界檢測能力較弱,顯示結(jié)果不清晰.NEDT和NEEDT相對于兩種常用的邊界探測器THDR和Theta圖法更加清晰和準確.
為了進一步驗證方法的有效性,建立更加復雜的重力異常模型二,其中包含正負異常,更加符合實際地質(zhì)情況.模型體的大小和位置與圖2相同.唯一不同的是棱柱體1和3的剩余密度為200kg·m-3;棱柱體2和4的剩余密度為-200kg·m-3.合成的重力異常見圖5.圖6為計算的梯度張量和相應的邊界檢測結(jié)果.可以看出,當同時包含正負異常時,傳統(tǒng)的邊界探測器THDR和Theta圖雖然能很好地識別地質(zhì)體的邊界,但是引入了一些額外的錯誤邊界.然而,本文提出的均衡探測器能很好地識別異常體的邊界,沒有引入錯誤的邊界信息,能更好地完成邊界識別工作.其中,NEDT中的p值為0.1;NEEDT中的p值為0.08.
圖2 模型體平面圖(a)和三維圖(b)Fig.2 Plan view (a) and 3D view (b) of synthetic model
圖3 合成重力異常圖,棱柱體的剩余密度為200 kg·m-3Fig.3 Synthetic gravity anomaly model. The contrasted densities of prisms are 200 kg·m-3
為了驗證方法的穩(wěn)定性,對圖6中模型二的各個重力梯度值增加5%的高斯白噪聲,圖7a—7c為加入噪聲后的重力梯度值.圖7d—7i為噪聲數(shù)據(jù)的邊界檢測結(jié)果.其中,NEDT中的p值為0.1,NEEDT中的p值為0.0.001.可以看出,EEDT相對于EDT受噪聲的影響較大,主要是由于EEDT利用了位場數(shù)據(jù)的三階導數(shù),放大了噪聲的影響;由于噪聲的干擾,THDR已經(jīng)無法給出地質(zhì)體的邊界,而NEEDT和Theta圖法受到噪聲的影響較大.相對而言,NEDT抗噪能力強,更加穩(wěn)定.
圖5 合成的重力異常,棱柱體1和3的剩余密度為 200 kg·m-3,棱柱體2和4的剩余密度為-200 kg·m-3Fig.5 Synthetic gravity anomaly model. The contrasted densities of prisms 1 and 3 are 200 kg·m-3, the contrasted densities of prisms 2 and 4 are -200 kg·m-3
圖6 同圖4,但為模型二的重力梯度和不同方法的邊界識別結(jié)果Fig.6 Same as Fig.4, but for gravity gradients and the edges detected by different methods of model two
圖7 同圖4,但為加入噪聲后的模型二的重力梯度和不同方法的邊界識別結(jié)果Fig.7 Same as Fig.4, but for gravity gradients and the edges detected by different methods of model two add 5% Gaussian noise
圖8 同圖4,但為加拿大圣喬治灣地區(qū)全張量重力梯度值邊界檢測結(jié)果Fig.8 Same as Fig. 4, but for edges of full tensor gravity gradients anomalies in St. Georges Bay, Canada
圖9 朱日和地區(qū)磁異常邊界檢測結(jié)果 (a)化極后的磁異常數(shù)據(jù);(b)異常的EDT結(jié)果;(c)異常的EEDT結(jié)果;(d)異常的NEDT結(jié)果; (e)異常的NEEDT結(jié)果;(f)異常的THDR結(jié)果;(g)異常的Theta圖結(jié)果.Fig.9 Edge results of magnetic anomalies in the Zhurihe area (a) Reduction to pole of magnetic data; (b) EDT of the data in (a); (c) EEDT of the data in (a); (d) NEDT of the data in (a); (e) NEEDT of the data in (a); (f) THDR of the data in (a); (g) Theta map of the data in (a).
利用本文提出的方法對磁異常數(shù)據(jù)進行處理前,需要對磁異常數(shù)據(jù)進行化極處理,因為磁數(shù)據(jù)及其導數(shù)均受磁化方向的干擾.因此采用化極后的數(shù)據(jù)進行解釋所獲得的結(jié)果將更加準確(Li, 2006).
為了驗證方法在實際中的應用效果,分別對加拿大圣喬治灣的全張量重力梯度張量數(shù)據(jù)和中國朱日和地區(qū)的磁異常數(shù)據(jù)進行邊界檢測.
圖8a—8c為Air-FTG全張量重力梯度儀在圣喬治灣測得的重力梯度張量數(shù)據(jù),飛行線距為500 m,聯(lián)絡線距為5000 m,飛行的高度為100 m;圖8d—8i為識別的邊界結(jié)果.可以看出,非均衡邊界探測器EDT和EEDT的識別能力較差,識別的邊界結(jié)果不夠清晰,且不能同時顯示不同振幅大小的異常的邊界.常用的均衡探測器THDR的顯示能力較差,不能顯示地質(zhì)體邊界.Theta圖法可以很好地顯示地質(zhì)體的邊界,但是識別出來的邊界并不十分清晰.而NEDT和NEEDT識別出來的邊界相對于Theta圖的結(jié)果更加清晰、準確,尤其是NEEDT識別出來的邊界更清晰準確.其中,NEDT中的p值為0.005,NEEDT中的p值為0.001.
圖9a為朱日和地區(qū)化極后的磁異常數(shù)據(jù),勘探區(qū)域面積為73 km×117 km,采樣間隔為20 m,飛行高度為1500 m,飛行線距為500 m.勘探區(qū)域為新元古代超層序,主要由大陸沉積物組成,除一些富含鐵的沙堰堤壩外基本不含任何鐵磁性物質(zhì).磁異常數(shù)據(jù)主要由SE-NW向的堤壩產(chǎn)生的近線性異常和鐵異常構(gòu)成.圖中的黑色點線為已知的構(gòu)造.相應的梯度值通過數(shù)值計算得到(Minkus and Hinojosa, 2001).圖9b—9g為相應的邊界檢測結(jié)果.可以看出,NEDT和NEEDT的邊界識別效果明顯好于其他方法的,能更加清晰地顯示邊界信息.其中,NEDT和NEEDT中的p值為0.005.
本文針對位場梯度張量數(shù)據(jù)包含多個信息成分的優(yōu)勢,提出方向總水平導數(shù)和加強方向總水平導數(shù).根據(jù)方向總水平導數(shù)和加強方向總水平導數(shù)建立了兩個新的邊界探測器,并對其進行歸一化處理,以便能同時檢測不同振幅大小異常的邊界.當異常數(shù)據(jù)同時包含正負異常時,本文提出的新探測器能很好地避免產(chǎn)生錯誤的邊界信息,而傳統(tǒng)的THDR和Theta圖將產(chǎn)生額外的錯誤邊界信息.通過模型試驗和實際數(shù)據(jù),證明了基于方向總水平導數(shù)和加強方向總水平導數(shù)定義的邊界探測器能很好地識別地質(zhì)體的邊界,且識別出來的邊界更加清晰、準確,避免了錯誤邊界信息的引入.
致謝 感謝加拿大自然資源部(Natural Resources Canada)提供圣喬治灣(St. Georges Bay)的Air-FTG全張量重力梯度測量數(shù)據(jù),并允許發(fā)表數(shù)據(jù)處理結(jié)果.感謝為評審本文所付出努力的專家.
Beiki M. 2010. Analytic signals of gravity gradient tensor and their application to estimate source location.Geophysics, 75(6): I59-I74.
Cooper G R J, Cowan D R. 2004. Filtering using variable order vertical derivatives.Computer&Geosciences, 30(5): 455-459.
Cooper G R J, Cowan D R. 2006. Enhancing potential field data using filters based on the local phase.Computer&Geosciences, 32(10): 1585-1591.
Cordell L. 1979. Gravimetric expression of Graben faulting in Santa Fe country and the Espanola Basin, New Mexico. New Mexico Geol. Soc. Guidebook, 30thField Conf., 1979: 59-64.
Cordell L, Grauch V J S. 1985. Mapping basement magnetization zones from aeromagnetic data in the San Juan basin, New Mexico. ∥Hinzc W J ed. The Utility of Regional Gravity and Magnetic Anomaly. Society of Exploration Geophysics, 181-197.
Evjen H M. 1936. The place of the vertical gradient in gravitational interpretations.Geophysics, 1(1): 127-136.
Forsberg R. 1984. A study of terrain reductions, density anomalies and geophysical inversion methods in gravity field modelling. Dept. of Geodetic Science and Surveying, Report 355, OhioState University.
Hsu S H, Sibuet J C, Shyu C T. 1996. High-resolution detection of geologic boundaries from potential-field anomalies: An enhanced analytic signal technique.Geophysics, 61(2): 373-386. Li X. 2006. Understanding 3D analytic signal amplitude.Geophysics, 71(2): L13-L16. Ma G Q, Huang D N, Yu P, et al. 2012. Application of improved balancing filters to edge identification of potential field data.ChineseJ.Geophys. (in Chinese), 55(12): 4288-4295, doi: 10.6038/j.issn.0001-5733.2012.12.040.
Marson I, Klingele E E. 1993. Advantages of using the vertical gradient of gravity for 3-D interpretation.Geophysics, 58(11): 1588-1595. Mikhailov V, Pajot G, Diament M, et al. 2007. Tensor deconvolution: A method to locate equivalent sources from full tensor gravity data.Geophysics, 72(5): 161-169.
Miller H G, Singh V. 1994. Potential field tilt—a new concept for location of potential field sources.JournalofAppliedGeophysics, 32(2-3): 213-217. Minkus K L, Hinojosa J H. 2001. The complete gravity gradient tensor derived from the vertical component of gravity: a Fourier transform technique.JournalofAppliedGeophysics, 46(3): 159-174.
Nabighian M N. 1984. Toward a three-dimensional automatic interpretation of potential field data via generalized Hilbert transforms: Fundamental relations.Geophysics, 49(6): 780-786.
Oru? B, Keskinsezer A. 2008. Structural setting of the northeastern
Biga Peninsula (Turkey) from tilt derivatives of gravity gradient tensors and magnitude of horizontal gravity components.PureandAppliedGeophysics, 165(9-10): 1913-1927. Oru? B, Sert?elik I, Kafadar ?, et al. 2013. Structural interpretation of the Erzurum Basin, eastern Turkey, using curvature gravity gradient tensor and gravity inversion of basement relief.JournalofAppliedGeophysics, 88: 105-113. Rajagopalan S, Milligan P. 1994. Image enhancement of aeromagnetic data using automatic gain control.ExplorationGeophysics, 25(4): 173-178.
Sertcelik I, Kafadar O. 2012. Application of edge detection to potential field data using eigenvalue analysis of structure tensor.JournalofAppliedGeophysics, 84: 86-94.
Verduzco B, Fairhead J D, Green C M, et al. 2004. New insights into magnetic derivatives for structural mapping.TheLeadingEdge, 23(2): 116-119.
Wijns C, Perez C, Kowalczyk P. 2005. Theta map: edge detection in magnetic data.Geophysics, 70(4): L39-L43.
Zhou W N, Du X J, Li J Y. 2013. The limitation of curvature gravity gradient tensor for edge detection and a method for overcoming it.JournalofAppliedGeophysics, 98: 237-242.
附中文參考文獻
馬國慶, 黃大年, 于平等. 2012. 改進的均衡濾波器在位場數(shù)據(jù)邊界識別中的應用. 地球物理學報, 55(12): 4288-4295, doi: 10.6038/j.issn.0001-5733.2012.12.040.
(本文編輯 胡素芳)
Using enhanced directional total horizontal derivatives to detect the edges of potential-field full tensor data
YUAN Yuan1,2,3, HUANG Da-Nian1*, YU Qing-Lu4
1CollegeofGeoExplorationScienceandTechnology,JilinUniversity,Changchun130026,China2SecondInstituteofOceanography,StateOceanicAdministration,Hangzhou310012,China3KeyLaboratoryofSubmarineGeoscience,StateOceanicAdministration,Hangzhou310012,China4SinopecGeophysicalResearchInstitute,Nanjing211103,China
With the ongoing development of the full potential field gradient techniques, more and more gradient tensor data have been widely used in geophysical exploration for its large amount of information and containing higher frequency signals than potential field data. The full potential field gradient technique can simultaneously measure six gradient components. Each component has its own geophysical meaning. The high frequency gradient tensor data can be used to delineate small scale anomalies.Edge detection is required for interpretation of potential field data, and has been used in exploration technology for discovery of mineral resources. The main geological edges are fault lines and the borders of geological or rock bodies of different density, magnetic nature, etc. On account of the higher-frequency signals, their interpretation enables high-resolution and detailed investigation of geological structures. Development of new methods is required to enable interpretation of these data, especially for edge detection for small geological structures.In order to make use of multiple component information, we define directional total horizontal derivatives and enhanced directional total horizontal derivatives and use them to define new edge detectors. However, the new defined edge detectors can balance the detected edge signals of different amplitude anomalies. In order to balance the edge signal amplitude, we present a normalization method, where we divide the new edge detectors by the vertical potential field gradient data. Besides, we introduce a constant parameter in the denominator of the normalization method, which can effectively avoid bringing some additional false edges when real geological bodies contain positive and negative anomalies simultaneously.These methods have been tested with synthetic data to verify that the new methods can delineate the edges of different amplitude anomalies clearly. The results show that new defined directional total horizontal derivatives and enhanced directional total horizontal derivatives have higher resolution than directional analytic signal. The normalized edge detectors can display the edges of large and small amplitude anomalies simultaneously, and avoid introducing additional false edges. To further test the stability, we demonstrate the new edge detectors with the model data corrupted with 5% Gaussian noise. Finally, we apply these methods to real full gravity gradient tensor data in St. Georges Bay, Canada and magnetic anomalies in Zhurihe area, China, and get good edge results.
Edge detection; Potential field gradient tensor; Enhanced directional total horizontal derivative; Normalization
10.6038/cjg20150730.
國家高技術研究發(fā)展計劃(863計劃)課題(2014AA06A613),國家重點基礎研究發(fā)展計劃項目(973項目)西南印度洋洋中脊熱液成礦過程與硫化物礦區(qū)預測的第五課題“硫化物礦區(qū)特征和找礦標志”(2012CB417305),大洋十二五項目西南印度洋脊合同區(qū)多金屬硫化物資源評價聯(lián)合資助.
袁園,男,1988年生,助理研究員,主要從事位場數(shù)據(jù)處理及解釋方面的研究.E-mail:yuanyuan12@mails.jlu.edu.cn
*通訊作者 黃大年,男,1958年生,教授,博士生導師,主要從事移動平臺探測數(shù)據(jù)處理與解釋及一體化軟件平臺開發(fā). E-mail:dnhuang@jlu.edu.cn
10.6038/cjg20150730
P631
2015-02-22,2015-06-16收修定稿
袁園, 黃大年, 余青露. 2015. 利用加強水平方向總水平導數(shù)對位場全張量數(shù)據(jù)進行邊界識別.地球物理學報,58(7):2556-2565,
Yuan Y, Huang D N, Yu Q L. 2015. Using enhanced directional total horizontal derivatives to detect the edges of potential-field full tensor data.ChineseJ.Geophys. (in Chinese),58(7):2556-2565,doi:10.6038/cjg20150730.