歸一化總水平導數(shù)法在位場數(shù)據(jù)解釋中的應(yīng)用

李麗麗，黃大年，韓立國

吉林大學地球探測科學與技術(shù)學院，長春 130061

1 引言

邊界識別是位場（重、磁）數(shù)據(jù)解釋中的主要任務(wù)之一，依據(jù)邊界識別結(jié)果能有效地識別出場源體的水平位置（Miller and Singh，1994；Blakely，1995）.重磁異常水平導數(shù)的極值與垂直導數(shù)的零值對應(yīng)地質(zhì)體的邊界（Verduzco et al.，2005；Cooper and Cowan，2006；2008；Wang et al.，2009），為此人們利用此特性構(gòu)建邊界識別濾波器來進行地質(zhì)體邊界的識別（Wijns et al.，2005；Ma，2013）.總水平導數(shù)法是一種進行地質(zhì)體邊界識別的常用方法（Nabighian，1972；Roest et al.，1992；馬國慶等，2012），其極大值點與地質(zhì)體的邊界相對應(yīng)，但該方法并不能給出較深地質(zhì)體的邊界.Ma和Li（2012）提出正則化總水平導數(shù)法，其通過對局部水平導數(shù)進行正則化處理使淺部與深部地質(zhì)體的邊界同時顯示出來，但該方法的輸出結(jié)果依賴于窗口尺寸的大小.王萬銀等（2009）給出了總水平導數(shù)的空間變化特征，從結(jié)果中可知利用總水平導數(shù)所識別出的邊界隨著地質(zhì)埋深的增加與理論邊界的誤差越大，且總水平導數(shù)法不能給出地質(zhì)體的深度信息.

為了提供地質(zhì)體的水平位置和深度信息，本文提出歸一化總水平導數(shù)法，通過空間歸一化總水平導數(shù)的極大值點可獲得異常體的水平位置和深度，并推導出基于歸一化總水平導數(shù)的歐拉反褶積法來估算出地下地質(zhì)體的空間位置，兩種方法計算結(jié)果的相互驗證可有效地提高反演結(jié)果的準確性.通過理論模型試驗和實際數(shù)據(jù)應(yīng)用驗證了歸一化總水平導數(shù)法的應(yīng)用效果.

2 歸一化總水平導數(shù)法

總水平導數(shù)（Total Horizontal Derivative，THD）是進行地質(zhì)體邊界識別的常用方法，其表達式為：

其中，f為原始重力或磁異常.總水平導數(shù)僅能給出較淺地質(zhì)體的邊界，而對于較深地質(zhì)體的邊界分辨率較差，且總水平導數(shù)極值點所識別的邊界與地質(zhì)體埋深也存在關(guān)系，隨著埋深的增加所識別邊界與理論邊界的誤差也增大.歸一化總梯度法在地質(zhì)體中心點處取得最大值（肖鵬飛等，2006），依據(jù)此理論提出空間歸一化總水平導數(shù)法，可知歸一化總水平導數(shù)到達場源體邊界時表現(xiàn)為最大值，根據(jù)其最值點可判斷異常體的水平位置和深度，且在異常體位置時埋深對總水平導數(shù)法所識別邊界的影響可忽略不計，因此歸一化總水平導數(shù)法所獲得結(jié)果是最為準確的.本文采用幾何平均對總水平導數(shù)進行空間歸一化處理，空間歸一化總水平導數(shù)（Normalized Total Horizontal Derivative，NTHD）的表達式為：

本文還推導出基于歸一化總水平導數(shù)的歐拉反褶積法來估算地質(zhì)體的空間位置.常規(guī)歐拉反褶積方程（Thompson，1982；Dewangan et al.，2007）的表達式為

其中，x，y，z為已知的觀測點坐標；x0，y0，z0為待求的場源體中心坐標；分別為位場異常和其在x，y，z方向上的導數(shù)；B為未知背景場值；N為構(gòu)造指數(shù)，不同地質(zhì)構(gòu)造對應(yīng)特定的構(gòu)造指數(shù)（Reid et al.，1990；Barbosa et al.，1999）.表1列出了不同地質(zhì)構(gòu)造所對應(yīng)的構(gòu)造指數(shù).

表1 不同類型地質(zhì)體對應(yīng)的構(gòu)造指數(shù)Table 1 Corresponding structural indexes of geological bodies of different types

在利用式（3）進行場源體位置計算時需要給定構(gòu)造指數(shù)的大小，因此式（3）反演結(jié)果的準確度取決于給定的構(gòu)造指數(shù)與實際構(gòu)造指數(shù)的差距，但由于實際測量中區(qū)域構(gòu)造指數(shù)是難以獲取且不唯一的，為反演結(jié)果帶來了很大的不確定性，將構(gòu)造指數(shù)與位置參數(shù)一起作為未知數(shù)進行反演將是有效的解決手段（Fitzgerald et al.，2004；Salem and Ravat，2003）.計算式（3）在x和y方向上的導數(shù)

由于背景異常的變化較為平緩，因此其導數(shù)項較小可以忽略（Silva and Barbosa，2003）.式（4）和（5）可改寫為：

計算歸一化總水平導數(shù)在x，y和z方向的導數(shù)

根據(jù)式（12）可以得到地質(zhì)體的空間位置及其構(gòu)造指數(shù)，該方法相對常規(guī)歐拉反褶積法具有不受構(gòu)造指數(shù)選取誤差的影響，且有效地綜合了不同方向的水平導數(shù)，能獲得更多的構(gòu)造信息.對于式（12）反演得到的結(jié)果采用如下的篩選準則來獲得更加準確的結(jié)果.

（1）計算點（xi，yi）與利用該點值求解得到的坐標x0，y0之間的距離小于窗口長度，因為當窗口在場源附近時解的準確性較高；

（2）利用一個較小窗口濾除一些孤立的解；

（3）濾除獨立的求解結(jié)果后，針對已經(jīng)集中的解進行篩選，當求解結(jié)果與集中解均值的距離大于窗口時該解無效；

（4）集中解的個數(shù)小于一定數(shù)量時直接濾除，因為當解的個數(shù)小于一定數(shù)量時不認為是有效異常.

經(jīng)過上述的篩選后就已經(jīng)獲得比較準確的場源位置坐標x0，y0，z0.

3 理論模型試驗

利用長方體產(chǎn)生的剖面重力異常來驗證歸一化總水平導數(shù)法的應(yīng)用效果（圖1）.圖1a為中心位置分別為50m和130m，上、下頂埋深均為12m和15m，與圍巖的密度差為1g／cm3，寬度為10m的兩個板狀體所產(chǎn)生的重力異常.

大多數(shù)歸一化方法均采用算術(shù)平均作為歸一化函數(shù).圖1b為基于算術(shù)平均的歸一化總水平導數(shù)結(jié)果，其中白框代表地質(zhì)體的真實位置，可以看出其最大值對應(yīng)的地質(zhì)體的埋深為13.3m，也能準確地標識出地質(zhì)體的位置.采用公式（2）計算空間歸一化總水平導數(shù)，其結(jié)果如圖1c所示.從計算結(jié)果中可以看出歸一化總水平導數(shù)最大值的深度為13.4m，與地質(zhì)體的真實深度13.5m之間誤差很小，反演結(jié)果能準確地給出異常體的位置和深度信息，且得到的結(jié)果相對算術(shù)平均歸一化總水平導數(shù)結(jié)果更加聚焦，發(fā)散程序較小.圖1d為利用歸一化總水平導數(shù)歐拉反褶積法計算得到的結(jié)果，窗口大小為7.反演結(jié)果顯示異常體的頂面深度為12.2m，與異常體的真實空間位置相一致.通過該模型試驗可知空間歸一化總水平導數(shù)和歸一化總水平導數(shù)歐拉反褶積法能準確地完成異常的反演工作.

圖1 （a）長方體產(chǎn)生的重力異常；（b）算術(shù)平均歸一化總水平導數(shù)；（c）基于幾何平均的歸一化總水平導數(shù)；（d）歸一化總水平導數(shù)歐拉反褶積法反演結(jié)果Fig.1 （a）Synthetic gravity anomalies generated by a rectangular prism；（b）Normalized total horizontal derivative based on arithmetic mean；（c）Normalized total horizontal derivative based on geometric mean；（d）Inversion results computed by the Euler deconvolution of normalized total horizontal derivative

實際數(shù)據(jù)中噪聲是不可避免的干擾因素.圖2a為在圖1a所示異常中加入幅度為異常最大值5%噪聲后異常，利用本文提出的歸一化總水平導數(shù)法對圖2a所示異常進行解釋（圖2）.圖2b為基于算術(shù)平均的歸一化總水平導數(shù)結(jié)果，可以看出最大值對應(yīng)的地質(zhì)體埋深為12.9m，水平位置與地質(zhì)體的真實位置相一致，能準確地完成異常的反演工作.圖2c為利用公式（2）計算得到的空間歸一化總水平導數(shù)，根據(jù)其最大值判斷出異常的埋深為13.1m，水平位置與理論位置吻合較好，但由于受到噪聲的影響，埋深與理論深度之間的誤差加大.圖2d為歸一化總水平導數(shù)歐拉反褶積法計算結(jié)果，窗口大小為7，反演得到的異常體的頂面深度11.6m，反演得到的異常體的水平位置與真實位置吻合較好.總體來講，歸一化總水平導數(shù)法具有一定的抗噪聲干擾能力.

下面驗證一下本文方法在二維情況下的應(yīng)用效果.圖3a為埋深分別為11m和13m的正方體所產(chǎn)生的重力異常.圖3b為圖3a所示異常的總水平導數(shù)，從圖中可以看出總水平導數(shù)能清晰地識別出較淺地質(zhì)體的邊界，而所得到的較深地質(zhì)體的邊界則比較模糊，且相對理論邊界要大一些.圖3c為空間歸一化總水平導數(shù)結(jié)果，通過空間歸一化總水平導數(shù)空間不同位置的幅值特征可判斷出地質(zhì)體的埋深.為了更加清晰地顯示出地質(zhì)體的邊界，圖3d給出了3、5、7、9、11、13m 深度上歸一化總水平導數(shù)的極大值點的位置，通過該結(jié)果可以看出隨著下延深度的增加歸一化總水平導數(shù)所識別出的異常體的水平位置與真實邊界更加接近，通過最大值的分布可判定出正方體的埋深分別為11m和13m.利用歸一化總水平導數(shù)歐拉反褶積法估算異常體的位置參數(shù)，窗口大小為7×7.圖3e為利用該方法計算得到的異常體的位置分布，圖3f為該方法反演結(jié)果的三維顯示，反演得到的異常體的平均深度分別為10.7m和12.6m，結(jié)果與理論埋深相一致.

圖2 （a）含噪重力異常；（b）算術(shù)平均歸一化總水平導數(shù)；（c）式（2）計算得到的空間歸一化總水平導數(shù)；（d）歸一化總水平導數(shù)歐拉反褶積法反演結(jié)果Fig.2 （a）Noise－bearing gravity anomalies；（b）Normalized total horizontal derivative based on arithmetic mean；（c）Spatially normalized total horizontal derivative computed by Eq.（2）；（d）Inversion results computed by the Euler deconvolution of normalized total horizontal derivative

為了驗證方法的適用性，將其應(yīng)用于磁異常的解釋.圖4a為埋深分別為10m和12m板狀體在磁化傾角為70°時產(chǎn)生的磁異常.圖4b為埋深分別為10m和12m板狀體在垂直磁化時候產(chǎn)生的磁異常.圖4c為圖4a所示異常的空間歸一化總水平導數(shù)，根據(jù)其最大值依舊能判斷出地質(zhì)體的埋深，但是由于傾斜磁化的影響所獲得水平位置與真實水平位置存在一定的差距.圖4d為圖4a所示異常的空間歸一化總水平導數(shù)，標識出的異常體的水平位置與實際水平位置吻合較好.圖4e為圖4a所示異常歸一化總水平導數(shù)歐拉反褶積法得到的結(jié)果，異常體的深度分別為8.2m和10.1m；圖4f為圖4b所示異常歸一化總水平導數(shù)歐拉反褶積法得到的結(jié)果，異常體的平均深度為9.1m和11.1m.通過試驗可知，由于傾斜磁化的干擾造成反演結(jié)果與理論結(jié)果之間的差距較大，因此為了獲得更加準確的結(jié)果，在利用本文方法進行磁異常解釋前需要對數(shù)據(jù)進行化磁極處理.此外還可以發(fā)現(xiàn)即使在垂直磁化條件下反演結(jié)果與理論值之間還是存在一定的差距，這是由于異常體距離較近，受鄰近異常體干擾較大造成的.

4 實際數(shù)據(jù)應(yīng)用

將歸一化總水平導數(shù)方法用于實際磁測數(shù)據(jù)的解釋，由于水平導數(shù)受磁化方向的干擾而使識別結(jié)果會出現(xiàn)較大誤差（Ma，2013；王萬銀，2010），因此在對其進行計算之前首先對磁異常進行磁極計算.

圖5a為朱日和地區(qū)化極后磁異常，測量點距為250m.圖5b為空間歸一化總水平導數(shù)結(jié)果.可以看出異常體的走勢及異常體的大致埋深，深度在400～600m之間，且根據(jù)不同深度極大值點的變化可判斷出異常體的延深方向為東南.圖5c為歸一化總水平導數(shù)歐拉反褶積計算結(jié)果，窗口大小為5×5.可以看出異常體的深度大多集中在400～600m之間，與前一種方法得到的結(jié)果一致.圖5d為反演得到的異常體的構(gòu)造指數(shù)，構(gòu)造指數(shù)范圍為0.42～1.71，其均值為1.12，地下地質(zhì)體接近于脈狀分布.由于異常體的不規(guī)則形與異常體之間的相互干擾造成構(gòu)造指數(shù)并不是一個整數(shù)，很多作者也證明了這一點（Salem et al.，2008；Fedi et al.，2012）.

5 結(jié)論

本文提出歸一化總水平導數(shù)法，其通過兩種方式來完成異常的解釋工作，一提出空間歸一化總水平導數(shù)，通過其最大值點可獲得異常體的空間位置，二推導出基于歸一化總水平導數(shù)的歐拉反褶積法來完成異常體位置參數(shù)的反演，且避免了構(gòu)造指數(shù)選取誤差為反演結(jié)果帶來的不確定性.理論模型試驗證明空間歸一化總水平導數(shù)的最大值能準確地獲得異常體的空間位置，基于歸一化總水平導數(shù)的歐拉反褶積法能同時完成場源體位置和構(gòu)造指數(shù)的計算，與理論值之間的誤差較小，且歸一化總水平導數(shù)法具有一定的抗噪聲能力.此外試驗結(jié)果表明在利用歸一化總水平導數(shù)法進行磁異常解釋時，事先對磁異常進行化磁極處理可獲得更加準確的結(jié)果.將本文方法應(yīng)用于實際航磁數(shù)據(jù)的解釋，獲得了巖脈的大致延伸趨勢及深度.

圖3 （a）原始重力異常；（b）圖3a所示異常的總水平導數(shù)；（c）空間歸一化總水平導數(shù)；（d）空間歸一化總水平導數(shù)不同深度極大值點的三維顯示；（e）歸一化總水平導數(shù)歐拉反褶積法得到的異常體的位置；（d）歸一化總水平導數(shù)歐拉反褶積法反演結(jié)果的三維顯示Fig.3 （a）Original gravity anomalies；（b）Total horizontal derivative of data in （a）；（c）Spatially normalized total horizontal derivative；（d）3Dview of the maximum values of spatially normalized total horizontal derivative at different depths；（e）Locations of sources computed by Euler deconvolution of normalized total horizontal derivative；（f）3D view of results computed by Euler deconvolution of normalized total horizontal derivative

圖4 （a）磁傾角為70°時板狀體產(chǎn)生的磁異常；（b）垂直磁化時板狀體產(chǎn)生的磁異常；（c）圖4a所示數(shù)據(jù)的空間歸一化總水平導數(shù)；（d）圖4b所示數(shù)據(jù)的空間歸一化總水平導數(shù)；（e）圖4a所示數(shù)據(jù)的歸一化總水平導數(shù)歐拉反褶積法反演結(jié)果；（d）圖4b所示數(shù)據(jù)的歸一化總水平導數(shù)歐拉反褶積法反演結(jié)果Fig.4 （a）Synthetic magnetic anomalies of a prism with inclination of 70°；（b）Synthetic magnetic anomalies of a prism with inclination of 90°；（c）Spatially normalized total horizontal derivative of data in（a）；（d）Spatially normalized total horizontal derivative of data in （b）； （e）Inversion results of data in （a） computed by Euler deconvolution of normalized total horizontal derivative；（f）Inversion results of data in（b）computed by Euler deconvolution of normalized total horizontal derivative

附錄 位場數(shù)據(jù)向下延拓的穩(wěn)定算法

Fedi和Florio采用泰勒展開式來完成數(shù)據(jù)的向下延拓，其表達式為：

圖5 （a）朱日和地區(qū)化極后磁異常；（b）空間歸一化總水平導數(shù)；（c）歸一化總水平導數(shù)歐拉反褶積法得到的異常體的位置參數(shù)；（d）歸一化總水平導數(shù)歐拉反褶積法得到的異常體的構(gòu)造指數(shù)Fig.5 （a）Reduction to the north pole of magnetic anomalies in the Zhurihe area；（b）Spatially normalized total horizontal derivative of data in（a）；（c）Locations of data in（a）computed by Euler deconvolution of normalized total horizontal derivative；（d）Structural index estimated by Euler deconvolution of normalized total horizontal derivative

其中，T（x，y，h）是觀測面h上的異常，h為延拓高度，T（x，y，0）是觀測面上異常.該方法是通過垂直導數(shù)來進行向下延拓操作，但垂直導數(shù)的計算會明顯增大噪聲的干擾，且當點距較大時，垂直導數(shù)的計算是不穩(wěn)定的.

向上延拓是一種穩(wěn)定的運算，且受噪聲干擾較小，因此本文提出利用向上延拓運算來完成異常的向下延拓.向上延拓運算的泰勒展開式可表示為：

其中，T（x，y，－h(huán)）為觀測面－h(huán)上的異常.三階或更高階垂直導數(shù)值較小，且計算結(jié)果穩(wěn)定性差，因此設(shè)定展開式的項數(shù)為3，并將公式（A1）和（A2）相加后可得到：

由于垂直導數(shù)的計算會增大噪聲的干擾，為此采用Laplace方程來計算異常的二階垂直導數(shù)：

水平導數(shù)是空間域中進行，不會增大噪聲的影響.因此利用公式（A3）進行異常的向下延拓工作只需要計算異常的水平導數(shù)和向上延拓結(jié)果，可有效地增強向下延拓的穩(wěn)定性.但由于該方法所保留的泰勒展開式的項數(shù)較少，因此對于復雜異常所獲得的結(jié)果精度較低.為了獲得更加準確的結(jié)果，給出一種迭代的方法來計算向下延拓結(jié)果.

將公式（A3）計算得到的下界面異常T（x，y，h）向上延拓高度h得 到異常T1（x，y，0），理論上T（x，y，0）與T1（x，y，0）是相等的，但由于T（x，y，h）僅是一個近似值，因此兩者之間會存在一定的差距，本文通過下面的迭代過程來消除這一差距.計算觀測面上原始異常T（x，y，0）與T1（x，y，0）的差

利用公式（A3）將 ΔT1（x，y，0）向下延拓高度h，可以得到

其中，ΔT1（x，y，h）為異常 ΔT1（x，y，0）下延h后異常，ΔT1（x，y，－h(huán)）為異常 ΔT1（x，y，0）上延h后異常.ΔT1（x，y，h）作下界面異常的一個修正量，則下界面異?？筛膶憺?/p>

重復式（A5）—（A7）的計算過程，直至ΔTm（x，y，0）的均方差小于給定值，因此觀測面h上的最終異?？蓪憺椋?/p>

Barbosa V C F，Silva J B C，Medeiros W E.1999.Stability analysis and improvement of structural index estimation in Euler deconvolution.Geophysics，64（1）：48－60.

Blakely R J.1995.Potential Theory in Gravity and Magnetic Applications.Cambridge：Cambridge University Press.

Cooper G R J，Cowan D R.2006.Enhancing potential field data using filters based on the local phase.Computers＆Geosciences，32（10）：1585－1591.

Cooper G R J，Cowan D R.2008.Edge enhancement of potentialfield data using normalized statistics.Geophysics，73（3）：H1－H4.

Dewangan P，Ramprasad T，Ramana M V，et al.2007.Automatic interpretation of magnetic data using Euler Deconvolution with nonlinear background.PureAppl.Geophys.，164（11）：2359－2372.

Fedi M，F(xiàn)lorio G，Cascone L.2012.Multiscale analysis of potential fields by a ridge consistency criterion：the reconstruction of the Bishop basement.Geophys.J.Int.，188（1）：103－114.

Fitzgerald D，Reid A，Mcinerney P.2004.New discrimination techniques for Euler deconvolution.Computers＆Geosciences，30（5）：461－469.

Ma G，Li L L.2012.Edge detection in potential fields with the normalized total horizontal derivative.Computers＆Geosciences，41：83－87.

Ma G Q，Liu C，Huang D N，et al.2013.A stable iterative downward continuation of potential field data.J.Appl.Geophys.，98：205－211.

Ma G Q.2013.Edge detection of potential field data using improved local phase filter.Explor.Geophys.，44（1）：36－41.

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.

Miller H G，Singh V.1994.Potential field tilt—A new concept for location of potential field sources.J.Appl.Geophys.，32（2－3）：213－217.

Nabighian M N.1972.The analytic signal of two－dimensional magnetic bodies with polygonal cross－section：its properties and use for automated anomaly interpretation.Geophysics，37（3）：507－517.

Reid A B，Allsop J M，Granser H，et al.1990.Magnetic interpretation in three dimensions using Euler deconvolution.Geophysics，55（1）：80－91.

Roest W R，Verhoef J，Pilkington M.1992.Magnetic interpretation using the 3－D analytic signal.Geophysics，57（1）：116－125.

Salem A，Ravat D.2003.A combined analytic signal and Euler method （AN－EUL）for automatic interpretation of magnetic data.Geophysics，68（6）：1952－1961.

Salem A，Williams S，F(xiàn)airhead D，et al.2008.Interpretation of magnetic data using tilt－angle derivatives.Geophysics，73（1）：L1－L10.

Silva J B C，Barbosa V C F.2003.3DEuler deconvolution：Theoretical basis for automatically selecting good solutions.Geophysics，68（6）：1962－1968.

Thompson D T.1982.EULDPH：A new technique for making computer－assisted depth estimates from magnetic data.Geophysics，47（1）：31－37.

Verduzco B，F(xiàn)airhead J D，Green C M.2004.New insights into magnetic derivatives for structural mapping.TheLeading Edge，23（2）：116－119.

Wang W Y，Pan Y，Qiu Z Y.2009.A new edge recognition technology based on the normalized vertical derivative of the total horizontal derivative for potential field data.Appl.Geophys.，6（3）：226－233.

Wang W Y.2010.Spatial variation law of the extreme value positions of total horizontal derivative for potential field data.ChineseJ.Geophys.（in Chinese），53（9）：2257－2270.

Wijns C，Perez C，Kowalczyk P.2005.Theta map：Edge detection in magnetic data.Geophysics，70（4）：L39－L43.

Xiao P F，Li M，Xu S Z，et al.2006.Stable solution of normalized total gravity gradient.OilGeophysicalProspecting，41（5）：596－598.

附中文參考文獻

馬國慶，黃大年，于平等.2012.改進的均衡濾波器在位場數(shù)據(jù)邊界識別中的應(yīng)用.地球物理學報，55（12）：4288－4295.

王萬銀.2010.位場總水平導數(shù)極值位置空間變化規(guī)律研究.地球物理學報，53（9）：2257－2270.

肖鵬飛，李明，徐世浙等.2006.重力歸一化總梯度的穩(wěn)定解法.石油地球物理勘探，41（5）：596－598.