基于正則化思想的tilt-Euler法在邊緣深度反演中的應用

摘要：地質體邊緣深度在重、磁位場數(shù)據(jù)半定量解釋中起著至關重要的作用。由于重、磁異常及其各階導數(shù)均滿足歐拉齊次方程，tilt-Euler法在邊緣深度反演方面?zhèn)涫芮嗖A。然而，當重、磁異常的總水平導數(shù)或者總梯度模等于0時，傾斜角的一階導數(shù)無法計算，導致傾斜角不能滿足歐拉方程，tilt-Euler法無法使用。為了解決此問題，本文基于正則化思想，對傾斜角的一階導數(shù)進行修改，使得重、磁異常的總水平導數(shù)或者總梯度模等于0時，傾斜角的一階導數(shù)依然可以計算，修改后的傾斜角導數(shù)依然滿足歐拉方程，稱改進的方法為rtilt-Euler法；同時利用識別精度更高的歸一化總水平導數(shù)垂向導數(shù)（NVDR-THDR）邊緣識別方法對反演結果進行約束，剔除偏離邊緣位置的壞點。理論模型試驗結果表明，改進后的方法消除了重、磁異常總水平導數(shù)或者總梯度模很小或者等于0時，傾斜角導數(shù)無法計算以及反演結果不穩(wěn)定的問題。將該方法應用到澳大利亞奧林匹克壩氧化鐵銅金礦床邊緣深度反演中，反演結果顯示氧化鐵銅金礦床邊緣深度主要集中在0～100 m和100～200 m這兩個深度段內，與沉積物剖面顯示的礦床邊緣深度0～200 m相符，證明了該方法的有效性。

關鍵詞：地質體邊緣深度；重磁位場；正則化；rtilt-Euler；NVDR-THDR

doi：10.13278/j.cnki.jjuese.20230035

中圖分類號：P631

文獻標志碼：A

0引言

重磁位場勘探方法以不同巖、礦石之間物理性質（密度、磁性等）的差異為基礎，當不同地質體邊界處存在明顯的密度或磁性差異時，可以根據(jù)重、磁異常來推斷地質體邊界的空間位置2]。20世紀80年代以來，隨著位場資料處理技術的發(fā)展，出現(xiàn)了一些地質體邊緣深度反演新方法，主要包括歐拉（Euler）反褶積法、沃納（Werner）反褶積法、局部波數(shù)法、解析信號振幅法、tilt-depth法和曲率屬性反演法等。其中：沃納反褶積法、局部波數(shù)法、解析信號振幅法、tilt-depth法和曲率屬性反演法都是從二維模型推廣到三維模型中，因此在三維模型中的應用受到了理論基礎的限制；而歐拉反褶積法是基于歐拉齊次方程提出的，位場數(shù)據(jù)及其各階導數(shù)均滿足歐拉齊次方程，因此該方法在邊緣深度反演方面?zhèn)涫芮嗖A。

歐拉反褶積法是Thompson在1982年基于歐拉齊次方程首次提出的。1990年，Reid等對歐拉反褶積法做了進一步研究，將歐拉反褶積法推廣到三維空間，提出了基于歐拉方程的邊緣識別方法。1998年，Keating針對不規(guī)則重力數(shù)據(jù)分布導致解發(fā)散的問題，提出用正比于測點精度的誤差函數(shù)對歐拉反褶積法的結果進行加權，消除了部分發(fā)散的解。2003年，Salem等提出ANEUL法，將總梯度模和歐拉反褶積兩種方法結合，識別地磁或航磁異常體的邊界，但這種方法受高頻干擾較大。2006年，范美寧對歐拉反褶積進行了比較系統(tǒng)的研究，用不同類型的數(shù)據(jù)作為原始數(shù)據(jù)進行反演，使反演結果的收斂性有較大的改善。后續(xù)很多學者針對歐拉反褶積法解發(fā)散的問題做了大量的研究13]，改善了其反演精度和穩(wěn)定性。由于最開始歐拉反褶積法在使用時必須提前已知構造指數(shù)，而構造指數(shù)是與場源形狀有關的參數(shù)，當已知場源先驗信息較少時，構造指數(shù)的獲取非常困難，尤其是存在多個場源時，受到異常疊加效應的影響，構造指數(shù)幾乎是不可能得到的。前人關于構造指數(shù)的選取也做了大量的研究16]，反演效果有所改善，但這仍然制約著該方法在位場數(shù)據(jù)中的應用。因此許多學者基于歐拉齊次方程嘗試使用不同的數(shù)據(jù)源類型來消除反演結果對構造指數(shù)的依賴22]。Salem等將傾斜角（tilt angle）法和歐拉反褶積法結合起來，消除了反演結果對構造指數(shù)的依賴，稱之為tilt-Euler法。該方法在邊緣位置探測和邊緣深度反演方面得到了廣泛的應用26]。2019年，Liang等通過對tilt-Euler法的研究，發(fā)現(xiàn)當位場數(shù)據(jù)的總水平導數(shù)等于0時，傾斜角的導數(shù)無法計算，存在奇點，因此對該方法進行了改進，重新定義了傾斜角的表達式，解決了該問題。然而本文在研究過程中發(fā)現(xiàn)當總梯度模等于0時，傾斜角的導數(shù)依然無法計算，存在奇點，并且當總梯度模很小時，傾斜角的導數(shù)非常不穩(wěn)定，這勢必也會引起反演結果的不穩(wěn)定。歐拉反褶積法存在的另一個問題是解的篩選，由于實際資料存在噪聲，因此反演結果中必然存在許多虛假解。針對解的篩選問題，Salem等利用傾斜角總水平導數(shù)的峰值位置對反演數(shù)據(jù)范圍進行約束，但當異常的總梯度模等于0時，傾斜角的總水平導數(shù)仍然無法計算，存在奇點。

本文基于正則化思想，對傾斜角的導數(shù)進行修改，消除奇點，并且證明修改后傾斜角各個方向的導數(shù)依然滿足tilt-Euler方程；同時利用識別精度更高的歸一化總水平導數(shù)垂向導數(shù)（NVDR-THDR）邊緣識別方法30]對反演數(shù)據(jù)范圍進行約束，使得反演的結果集中在邊緣位置附近；利用模型試驗證明改進后的方法解決總梯度模等于0時奇點問題以及總梯度模很小時反演結果不穩(wěn)定性問題的有效性。最后，將該方法應用于澳大利亞奧林匹克壩氧化鐵銅金礦床的邊緣深度反演，為礦產(chǎn)資源儲量估算提供參考資料。

1基本原理

1.1tilt-Euler法

1.2基于正則化思想改進的tilt-Euler法（rtilt-Euler）

2模型試算

本次模型試算選取3個重力異常模型，模型參數(shù)見表1。通過對比tilt-Euler法和rtilt-Euler法的反演結果，驗證本文提出的正則化方法的有效性。各模型重力異常及邊界位置如圖2所示。

2.1模型A——單一直立六面體模型

圖3為模型A的邊緣深度反演結果及其分段直方圖。由邊緣深度反演結果可知，對于單一直立六面體來說，tilt-Euler法（圖3a）和rtilt-Euler法（圖3c）的反演結果基本相同。這是由于單一直立六面體在邊界處不存在奇點，同時也證明了改進方法的正確性。從反演結果的分段直方圖上可以看出，改進方法在深度上更靠近直立六面體邊緣深度的頂面位置，眾值集中在100～120 m的范圍內（圖3d），而原來的方法集中在110～130 m之間（圖3b）。由圖3e和圖3f可知，如果以NVDR-THDR對邊緣位置進行約束，其反演效果會更好，剔除了邊界外的其他壞點。但由圖3b、圖3d和圖3f可知，即使是使用NVDR-THDR對邊緣位置進行約束，邊界上的反演結果值雖然大部分靠近模型的頂面深度，但也有小部分點的深度偏離頂面深度較遠，最深達到了120 m。這實際上是因為tilt-Euler法最開始是用來反演孤立球體的中心深度的，當將tilt-Euler法推廣到邊緣深度反演中時，對模型是有要求的，理論上來說要求臺階模型在走向方向無限延伸，在垂直于走向的方向滿足一側無限延伸，同時臺階的下底深度要遠遠大于頂面埋深。如果滿足這些條件，那么反演的結果將不存在額外的結果，所有反演解都將集中在地質體邊緣頂面深度的位置。但在實際應用中，這些條件是不可能滿足的，為了模擬真實的地質情況，本文設計的模型A也不滿足這些條件。

2.2模型B——相同埋深的兩個組合直立六面體模型法（圖4c）在非相鄰邊界上反演結果基本相同，深度基本位于100～140 m之間（圖4b和圖4d）；在相鄰邊界處，rtilt-Euler法反演結果的連續(xù)性優(yōu)于tilt-Euler法。這是因為在形體的相鄰邊界處可能出現(xiàn)奇點，而rtilt-Euler法消除了奇點，提高了解的穩(wěn)定性，反演的深度位于200～360 m之間（圖4d）。如果以NVDR-THDR對邊緣位置進行約束，其反演結果（圖4e）會集中在模型體的邊緣位置，尤其是較窄的邊能夠得到很好的改善。另外，由圖4b、圖4d和圖4f可知，相較于模型A，模型B的反演結果具有兩個深度峰值，深度較淺的峰值（100～140 m）為直立六面體非相鄰的兩個邊的頂面埋深，深度較大的峰值（200～360 m）為相鄰兩個邊的反演結果，該峰值明顯大于模型邊緣的頂面埋深。這是由于當將tilt-Euler法推廣到邊緣深度反演中時，除了要求臺階模型在走向方向無限延伸、在垂直于走向的方向滿足一側無限延伸、臺階的下底深度遠遠大于頂面埋深外，還要求模型為孤立形體。所以當?shù)刭|體不是孤立形體時，受相鄰形體的疊加異常影響，反演的結果會偏深。

2.3模型C——不同埋深的多個組合直立六面體模型

圖5為模型C的邊緣深度反演結果。由邊緣深度反演結果可知，對于不同埋深的組合直立六面體來說，rtilt-Euler法的反演結果（圖5c）明顯優(yōu)于tilt-Euler法（圖5a）。這是由于當模型增多、埋深變化時，傾斜角導數(shù)奇點數(shù)量增多，tilt-Euler法反演的結果變得極為不穩(wěn)定。從深度的反演結果上來看，tilt-Euler法（圖5a）僅能刻畫埋深最淺形體遠離其他形體的那一條邊界，剩下的形體邊界均無法刻畫。rtilt-Euler法反演結果（圖5c）能夠清晰刻畫埋深較淺形體的邊界位置，對于埋深較大的形體能夠刻畫部分邊界；但受淺部形體異常的影響，隨著形體頂面埋深的增加，反演結果的精度越來越差，因此在實際使用中應對異常做適當?shù)姆蛛x。如果以NVDR-THDR對邊緣位置進行約束，rtilt-Euler法反演結果會集中在模型體的邊緣位置。另外，相較于模型A和模型B，模型C的反演結果較為復雜。tilt-Euler法反演的結果雖然深度峰值較多（圖5b），但結合圖5a可知，除了最淺的深度峰值外，其他的都是假峰值，與模型的邊緣位置并不對應。改進后的rtilt-Euler法反演結果具有兩個明顯的深度峰值（圖5d），較淺的峰值（60～80 m）對應的是埋深最小的模型，較深的峰值（180～300 m）對應另外兩個埋深較大模型。使用NVDR-THDR對邊緣位置進行約束后，其反演結果（圖5f）具有四個深度峰值：40～80 m、120～160 m、180～220 m、220～300 m。40～80 m的深度峰值對應埋深最淺的直立六面體的邊緣；120～160m和180～220 m的兩個深度峰值分別對應第二個直立六面體的窄邊和寬邊；220～300 m的深度峰值對應于埋深最大的直立六面體的邊緣。

3實際資料處理

澳大利亞奧林匹克壩氧化鐵銅金礦床（IOCG）是澳大利亞價值最大的礦產(chǎn)出口源，礦產(chǎn)主要包括鐵、金、銅和鈾。澳大利亞有兩個主要的IOCG省：Mount Isa Block和Gawler Craton。奧林匹克壩a、c、e分別為tilt-Euler、rtilt-Euler以及使用NVDR-THDR約束的rtilt-Euler的反演結果；b、d、f分別為tilt-Euler、rtilt-Euler以及使用NVDR-THDR約束的rtilt-Euler反演結果對應的分段直方圖。

CuUAuAg礦床是太古宙—元古代Gawler Craton邊緣的幾個氧化鐵銅金鈾（IOCGU）礦床中最大的礦床，如圖6所示，F(xiàn)e礦床的平均深度為0～1 000 m，最深部分達1 500 m。奧林匹克壩與高重力異常有關，通常由大量赤鐵礦角礫巖引起。因此，基于重力異常可以反演該礦床的邊緣深度，對礦床儲量評估提供參考。

研究區(qū)的布格重力異常（圖7a）來自于SARIG（sarig.pir.sa.gov.au）。為了得到奧林匹克壩氧化鐵銅金礦床引起的重力異常，本文使用最小曲率分離方法對布格重力異常進行位場分離，得到研究區(qū)的剩余布格重力異常（圖7b），認為該重力異常是由奧林匹克壩氧化鐵銅金礦床引起的。同時使用NVDRTHDR方法識別了該礦床的邊緣位置。以NVDRTHDR方法識別結果作為約束，使用rtilt-

Euler法反演礦體的邊緣深度（圖8a）。反演結果顯示氧化鐵銅金礦床邊緣深度主要集中在0～200 m之間，反演結果的分段直方圖（圖8b）也顯示深度主要集中在0～100 m和100～200 m這兩個深度段內。圖6中沉積物剖面受起始深度限制，顯示的礦床頂面邊緣深度在200 m左右；但根據(jù)剖面圖中沉積物的形態(tài)，可以推測其在200 m以淺明顯還有延伸。因此rtilt-Euler法反演礦體的邊緣深度與實際情況相符，證明了該方法的有效性和正確性。

4結論與展望

根據(jù)模型試算和實際資料處理結果，本文得出以下幾點結論：

1）對于單一直立六面體來說，tilt-Euler法和rtilt-Euler法的反演結果基本相同；這是由于單一直立六面體在邊界處不存在奇點，同時也證明了改進后的方法的正確性。為了消除奇點，引入的正則化因子通常情況下取1即可。

2）對于不同埋深的組合直立六面體來說，rtilt-Euler法的反演結果明顯優(yōu)于tilt-Euler法。從深度的反演結果上來看，tilt-Euler法僅能刻畫埋深最淺形體遠離其他形體的那一條邊界，剩下的形體邊界均無法刻畫；rtilt-Euler法反演結果能夠清晰刻畫埋深較淺形體的邊界位置，對于埋深較大的形體也能刻畫部分邊界。另外，受淺部形體異常的影響，隨著形體頂面埋深的增加，反演結果的精度越來越差，因此在實際應用中應對異常做適當?shù)姆蛛x。

3）使用NVDR-THDR邊緣識別方法對反演結果進行約束，無論是單一模型還是組合模型，邊緣深度的反演結果都位于形體邊界附近，剔除了形體邊界以外的壞點，反演效果得到了改善。

模型試算和實際資料處理表明rtilt-Euler法反演邊緣深度的結果優(yōu)于tilt-Euler法，證明了該方法消除奇點的有效性。實際上邊緣深度反演還存在許多問題，如數(shù)據(jù)源的適應性、方法的穩(wěn)定性、解的篩選問題等，本文只是解決了tilt-Euler法的奇點問題和解的篩選問題，而數(shù)據(jù)源的適應性和方法的穩(wěn)定性問題本文并沒有討論，但這兩個問題對邊緣深度反演非常重要。因此在未來的研究工作中，應從理論模型的正演公式出發(fā)，研究不同輸入數(shù)據(jù)的適應性和影響反演方法穩(wěn)定性的因素，進而對邊緣深度反演方法進行改進和創(chuàng)新。

參考文獻（References）：

1.趙強， 管彥武， 盧鵬羽， 等. 基于歸一化特征值的重磁數(shù)據(jù)邊界識別方法及其在西藏古堆地熱勘探中的應用. 吉林大學學報（地球科學版）， 2023， 53（2）： 578588.

Zhao Qiang， Guan Yanwu， Lu Pengyu， et al. Edge Detection Method of Gravity and Magnetic Data Based on Normalized Eigenvalue and Its Application of Geothermal Exploration in Gudui， Tibet. "Journal of Jilin University （Earth Science Edition）， 2023， 53（2）： 578588.

2.魯寶亮， 蘇子旺， 趙志剛， 等. 基于衛(wèi)星重力異常的南海南部中生界反演. 吉林大學學報（地球科學版）， 2023， 53（1）： 261273.

Lu Baoliang， Su Ziwang， Zhao Zhigang， et al. Mesozoic Distribution Derived from Satellite Gravity Anomaly in the Southern South China Sea. Journal of Jilin University （Earth Science Edition）， 2023， 53（1）： 261273.

Thompson D T. EULDPH： A New Technique for Making Computer Assisted Depth Estimates from Magnetic Data. Geophysics， 1982， 47（1）： 3137.

Reid A B， Allsop J M， Granser H， et al. Magnetic Interpretation in Three Dimensions Using Euler Deconvolution. Geophysics， 1990， 55（1）： 8091.

Keating P B. Weighted Euler Deconvolution of Gravity Data. Geophysics， 1998， 63（5）： 15951603.

Salem A，Ravat D. A Combined Analytic Signal and Euler Method （An-Eul） for Automatic Interpretation of Magnetic Data. Geophysics， 2003， 68（6）： 19521961.

3.范美寧. 歐拉反褶積方法的研究與應用. 長春：吉林大學， 2006.

Fan Meining. The Study and Application of Euler Deconvolution Method. Changchun： Jilin University， 2006.

Fedi M， Florio G， Quarta T A M. Multi-Ridge Analysis of Potential Fields： Geometric Method and Reduced Euler Deconvolution. Geophysics， 2009， 74（4）： 5365.

4.馬國慶， 孟令順， 李麗麗. 龍門山及鄰區(qū)斷裂分布及地震前后斷裂形態(tài)差異. 吉林大學學報（地球科學版）， 2012， 42（2）： 519525， 553.

Ma Guoqing， Meng Lingshun， Li Lili. Fault Distribution of Longmenshan and Adjacent Regions and Fault Morphological Differences Before and After the Earthquake. Journal of Jilin University （Earth Science Edition）， 2012， 42（2）： 519525， 553.

Reid A B， Ebbing J， Webb S J. Avoidable Euler Errors：The Use and Abuse of Euler Deconvolution Applied to Potential Fields. Geophysical Prospecting， 2015， 62（5）： 11621168.

Zhou W， Nan Z， Li J. Self-Constrained Euler Deconvolution Using Potential Field Data of Different Altitudes. Pure and Applied Geophysics， 2016， 173（6）： 20732085.

5.姚剛， 董向欣， 李麗麗，等. 東海陸架盆地構造劃分的高精度重磁解釋技術.吉林大學學報（地球科學版）， 2018， 48（2）： 517524.

Yao Gang， Dong Xiangxin， Li Lili， et al. High-Precision Gravity and Magnetic Interpretation Technique in Tectonic Division of East China Sea Shelf Basin. Journal of Jilin University （Earth Science Edition）， 2018， 48（2）： 517524.

6.喬中坤， 馬國慶， 于平， 等. 基于歐拉反褶積法的無人機航磁找礦應用. 吉林大學學報（地球科學版）， 2021， 51（2）： 552560.

Qiao Zhongkun， Ma Guoqing， Yu Ping， et al. Application of UAV Aeromagnetic Prospecting Based on Euler Deconvolution. Journal of Jilin University （Earth Science Edition）， 2021， 51（2）： 552560.

Barbosa V C F， Silva J B C， Medeiros W E. Stability Analysis and Improvement of Structural Index Estimation in Euler Deconvolution. Geophysics， 1999， 64 （1）： 4860.

Fedi M. DEXP： A Fast Method to Determine the Depth and the Structural Index of Potential Fields Sources. Geophysics， 2007， 72（1）： I1I11.

Melo F F， Barbosa V C F. Correct Structural Index in Euler Deconvolution Via Base-Level Estimates. Geophysics， 2018， 83 （6）： J87J98.

Hsu S K. Imaging Magnetic Sources Using Euler’s Equation. Geophysical Prospecting， 2010， 50（1）： 1525.

Salem A，Ravat D， Smith R， et al. Interpretation of Magnetic Data Using an Enhanced Local Wavenumber （ELW） Method. Geophysics， 2005， 70（2）： L7L12.

Salem A， Williams S， Fairhead D， et al. Interpretation of Magnetic Data Using Tilt-Angle Derivatives. Geophysics， 2008， 73（1）： L1L10.

Ma G. Combination of Horizontal Gradient Ratio and Euler （HGR-EUL） Methods for the Interpretation of Potential Field Data. Geophysics， 2013， 78（5）：J53J60.

Cooper G R J， Whitehead R C. Determining the Distance to Magnetic Sources. Geophysics， 2016， 81（2）： J25J34.

Cooper G R J. Determining the Depth and Location of Potential Field Sources Without Specifying the Structural Index. Arabian Journal of Geosciences， 2017， 10（19）： 438445.

Rahman M M， Ullah S E. Constrained Interpretation of Aeromagnetic Data Using Tilt-Angle Derivatives from North-Western Part of Bangladesh. Global Advanced Research Journal of Engineering， Technology and Innovation， 2013，2 （5）： 196204.

Khan S D， Fathy M S， Abdelazeem M. Remote Sensing and Geophysical Investigations of Moghra Lake in the Qattara Depression， Western Desert， Egypt. Geomorphology， 2014， 207（1）： 1022.

Yandjimain J， Ndougsa-Mbarga T， Meying A， et al. Combination of Tilt-Angle and Euler Deconvolution Approaches to Determine Structural Features from Aeromagnetic Data Modeling over Akonolinga-Loum Area （Centre-East， Cameroon）. International Journal of Geosciences， 2017， 8（7）： 925947.

7.羅瀟， 王彥國， 葛坤朋， 等. 基于方向tilt-Euler的三維磁數(shù)據(jù)快速反演. 地球物理學報， 2021， 64（6）： 21272140.

Luo Xiao， Wang Yanguo， Ge Kunpeng， et al. Fast Inversion of 3D Magnetic Data Based on the Directional Tilt-Euler Method. Chinese Journal of Geophysics， 2021， 64（6）： 21272140.

Liang H，Henglei Z， Sence S， et al. An Improved Tilt-Euler Deconvolution and Its Application on a FePolymetallic Deposit. Ore Geology Reviews， 2019， 114： 103114103122.

Wang W， Pan Y，Qiu Z. A New Edge Recognition Technology Based on the Normalized Vertical Derivative of the Total Horizontal Derivative for Potential Field Data. Applied Geophysics， 2009， 6（3）： 226233.

8.何濤， 王萬銀， 黃金明， 等. 正則化方法在比值類位場邊緣識別方法中的研究. 物探與化探， 2019， 43（2）： 308319.

He Tao， WangWanyin， Huang Jinming， et al. The Research of the Regularization Method in the Ratio Methods of Edge Recognition by Potential Field. Geophysical and Geochemical Exploration， 2019， 43（2）： 308319.

Zhu Y， Wang W， Farquharson C G， et al. Normalized Vertical Derivatives in the Edge Enhancement of Maximum-Edge-Recognition Methods in Potential Fields. Geophysics， 2021， 86（4）： G23G34.

Miller H G， Singh V. Potential Field Tilt： A New Concept for Location of Potential Field Sources. Journal of Applied Geophysics， 1994， 32（2/3）： 213217.

Verduzco B， Fairhead J D， Green C M， et al. The Meter Reader-New Insights into Magnetic Derivatives for Structural Mapping. The Leading Edge， 2004， 23（2）： 116119.

Austin J， Foss C. Rich， Attractive and Extremely Dense： A Geophysical Review of Australian IOCGs// ASEG Extended Abstracts. Brisbane： ASEG， 2012： 14.

Kirchenbaur M， Maas R， Ehrig K， et al. Uranium and Sm Isotope Studies of the Supergiant Olympic Dam CuAuUAg Deposit， South Australia. Geochimica et Cosmochimica Acta， 2016， 180： 1532.

9.紀曉琳， 王萬銀， 邱之云. 最小曲率位場分離方法研究. 地球物理學報， 2015， 58（3）： 10421058.

Ji Xiaolin， Wang Wanyin， Qiu Zhiyun. The Research to the Minimum Curvature Technique for Potential Field Data Separation. Chinese Journal of Geophysics， 2015， 58（3）： 10421058.