張恒磊，胡祥云，劉天佑

1 中國地質(zhì)大學（武漢）地球物理與空間信息學院，武漢 430074

2 中國地質(zhì)大學（武漢）構(gòu)造與油氣資源教育部重點實驗室，武漢 430074

1 引 言

磁法勘探作為地球物理勘探的一個重要分支，其發(fā)展最早，理論成熟，廣泛應(yīng)用于油氣與固體礦產(chǎn)勘探、環(huán)境與工程地球物理勘查、區(qū)域與深部構(gòu)造研究等領(lǐng)域.對磁異常進行反演是地質(zhì)-地球物理解釋中的一項重要工作，通過它可以獲取場源位置、埋深等參數(shù)，并與地質(zhì)資料相結(jié)合，可以了解更加豐富、全面的地下地質(zhì)信息，提高解釋質(zhì)量.

在磁異常反演[1］中，除了一些相對簡單的經(jīng)驗切線法，大多需要經(jīng)過反復(fù)的正反演計算，同時需要提供物性參數(shù).為了約束反演的多解性，常常需要充分結(jié)合地質(zhì)信息才能獲得滿意的結(jié)果，其過程費時費力.尤其在地質(zhì)工作初期，物性參數(shù)缺乏、地質(zhì)因素不甚明了，很難實現(xiàn)磁異常的有效反演.2007年，Salem等[2］在Tilt梯度分析場源邊界的基礎(chǔ)上，重新分析了Tilt梯度的含義，并通過理論推導(dǎo)，擴展了Tilt梯度的物理意義[3］，即 Tilt＿depth反演方法.研究指出，在簡單理想模型數(shù)值實驗中，Tilt＿depth方法無需提供先驗信息，即可以精確的反演出場源的上頂埋深，他們利用該方法對坦桑尼亞東南部的航磁異常進行分析，推斷解釋的基底構(gòu)造與地質(zhì)填圖結(jié)果相吻合，證明了方法的簡潔有效.2008年，Salem等[4］撰文討論了基于磁異常梯度的解釋方法，并把Tilt＿depth方法應(yīng)用于納米比亞中北部的航磁異常處理中，得到了客觀的解釋結(jié)果；2009年，Lee等[5］討論了Tilt＿depth方法的應(yīng)用前提，隨后Lahti和Karinen于2010年撰文討論了磁異常幅值大小對深度反演的影響并進一步分析了Tilt＿depth方法的理論條件[6］.2011年，F(xiàn)airhead等[7］對Tilt方法的應(yīng)用效果進行了討論，并成功應(yīng)用于沉積盆地的結(jié)構(gòu)劃分與深度反演.目前，Tilt＿depth方法得到了廣泛的應(yīng)用[8-11］；在國內(nèi)研究領(lǐng)域，也有不少學者對Tilt＿depth方法進行跟蹤報道[12-15］，體現(xiàn)了這種方法的優(yōu)越性與廣泛適用的特點.

本文指出了Tilt＿depth方法存在的問題，在此基礎(chǔ)上，推導(dǎo)了基于垂向二階導(dǎo)數(shù)的磁異常非參數(shù)反演，以期得到更好的邊界識別效果以及更精確的深度反演.同時，筆者討論了本文方法的穩(wěn)定計算問題，提出采用向上延拓的處理方法代替?zhèn)鹘y(tǒng)的低通濾波消除高頻干擾，實現(xiàn)穩(wěn)定計算.

2 方法原理

2.1 Tilt梯度與Tilt＿depth

1994年，Miller和Singh首次給出了Tilt梯度的數(shù)學定義[2］，并用于位場異常邊界分析.它一定程度上克服了常規(guī)位場導(dǎo)數(shù)對深部異常、弱異常反應(yīng)欠佳的弊端.Tilt梯度值對于場源的深度不敏感，計算結(jié)果不受場源埋深情況制約，因此它能夠很好的探測出具有不同埋深的復(fù)雜場源體的邊界.Tilt梯度的公式是根據(jù)二維解析信號發(fā)展而來的，它定義為位場異常的垂直導(dǎo)數(shù)與水平導(dǎo)數(shù)比值的反正切，表達式為：

式中，f表示位場異常，?f／?z指位場異常的垂向?qū)?shù)，其值在場源上方為正，在場源外側(cè)為負，在場源邊界位置附近其值為零；?f／?h指位場異常的水平導(dǎo)數(shù)，對于平面異常

2007年，Salem等用二度直立板狀體垂直磁化時的磁異常公式，導(dǎo)出場源頂深與異常的兩個一階導(dǎo)數(shù)之比（垂向一階導(dǎo)數(shù)與水平一階導(dǎo)數(shù)）的關(guān)系式[3］：

式中x指空間坐標，z0指二度直立板狀體的上頂深度.

分析可知，這個關(guān)系對二度直立板狀體垂直磁化情形是精確的，其它情形都是近似的.對剖面數(shù)據(jù)，由兩個一階導(dǎo)數(shù)之比的反正切繪制剖面圖，圖上兩點（-45°和45°）之間的水平距離為場源頂深z0的2倍，0°值點對應(yīng)場源邊界；對平面數(shù)據(jù)，由兩個一階導(dǎo)數(shù)之比的反正切繪制平面等值線圖，由該圖上兩條等值線（-45°和45°）之間的距離的一半確定場源頂深，由0°等值線確定場源邊界.

2.2 V2D＿depth非參數(shù)反演

Salem等[2］基于垂向一階導(dǎo)數(shù)的Tilt方法不利于區(qū)分疊加異常，同時它識別的異常邊界與實際情況差別較大.本文把前人得到的一階導(dǎo)數(shù)關(guān)系式推廣至二階形式，導(dǎo)出場源頂深與異常的兩個二階導(dǎo)數(shù)之比（垂向二階導(dǎo)數(shù)與水平二階導(dǎo)數(shù)）的關(guān)系式，即V2D＿depth方法.

二度棱柱體模型產(chǎn)生的磁異常計算公式為[16］：ΔT（P）=

圖1 二度模型導(dǎo)數(shù)計算坐標示意圖Fig.1 Notation used to derive the magnetic anomaly due to a two-dimensional

其中，k表示磁化率，F(xiàn)表示地球磁場值的大小，c=1-cos2i·sin2A（A為測線磁方位角，i為地磁傾角），φ=2I-d-90，I表示有效磁化傾角 （tanI=tani／cos A），其它參數(shù)詳見圖1中說明.

1972年，Nabighian推導(dǎo)了（3）式的垂向梯度公式，如下：

對于垂直物性邊界的地質(zhì)體，d=90°，對于化極磁異常，I=90°，故有φ=0°.（4）式可以進一步表示為

進一步推導(dǎo)可得如下關(guān)系式：

上述公式中Tzz表示ΔT磁異常的垂向二階導(dǎo)數(shù)，Tzh表示ΔT磁異常的垂向一階導(dǎo)數(shù)Tz的水平導(dǎo)數(shù)，對于平面磁異常，

根據(jù)式（6）和式（8）即可以得到如下比值關(guān)系：

由此可以確定基于垂向二階導(dǎo)數(shù)的非參數(shù)反演算法，即

2.3 分析討論與計算流程

在位場反演中，有一些方法利用了理論上場與場源特征參數(shù)之間的某種特殊關(guān)系，直接利用位場異常的特征點進行反演，即所謂的直接法，比如生產(chǎn)實踐中常用的經(jīng)驗切線法等.與這些簡單的方法相似，Tilt＿depth方法與本文方法反演過程無需輸入?yún)?shù)，筆者稱這種無需先驗信息即無需輸入?yún)?shù)的反演方法為“非參數(shù)反演方法”.

分析Tilt＿depth與本文方法，有以下幾方面性質(zhì).

（1）Salem 等[3］指出 Tilt＿depth方法容易受疊加異常干擾，±45°等值線間距在場源四周有疏有密并且差異很大，直接影響埋深參數(shù)的反演，限制了方法的使用.

（2）在文獻[17-18］中，筆者利用類似“廣義的水平二階導(dǎo)數(shù)”類方法計算重磁源邊界，并通過理論與實踐證明了方法的有效性.聯(lián)系位場拉普拉斯方程Txx+Tyy+Tzz=0，因此本文的公式（10）得到的邊界與文獻[17-18］有相似的分析結(jié)果，它比Tilt＿depth方法分析的邊界更接近真實情況.

（3）V2D＿depth方法也是通過計算 V2D＿Depth對應(yīng)的兩條等值線（特征等值線）之間的空間距離來反演深度.比如，x=±z0，x=±z0／2或者x=±z0／4兩條等值線，分別對應(yīng)±arctan（1）、±arctan（4／5）和±arctan（8／17）等值線，將它們之間的空間距離分別乘以0.5、1、2即可以得到反演的上頂深度.

（4）基于導(dǎo)數(shù)的計算方法對噪聲干擾異常敏感，本文方法也不例外.本文中筆者建議采用向上延拓的方式避免高頻干擾引起的計算振蕩現(xiàn)象.

基于以上分析，本文算法的計算流程如下：

（1）將獲得的磁異常進行化到地磁極計算；

（2）對化極磁異常向上延拓消除高頻干擾——上延的高度視噪聲水平靈活設(shè)定；

（3）計算異常的垂向二階導(dǎo)數(shù)、垂向一階導(dǎo)數(shù)的水平導(dǎo)數(shù)；

（4）按照公式（10）將上述導(dǎo)數(shù)進行比值計算，并取反正切，得到本文的θ（V2D＿Depth）值；

（5）計算θ（V2D＿Depth）特征等值線之間距，減去相應(yīng)的延拓高度，得到地質(zhì)體的上頂深度.

3 理論模型分析

為對本文方法進行驗證，我們做了模型數(shù)值實驗：首先設(shè)計了一個組合模型（模型位置如圖2黑色線），4個形體均為向下無限延深，磁化傾角I取90°，其它參數(shù)如表1所示.正演得到的磁異常如圖2所示.

圖2 理論模型正演ΔT磁異常Fig.2 Forward magnetic anomaly by theoretical model

表1 模型參數(shù)Table 1 The model parameters

為了獲取四個場源的水平位置及埋深，我們采用Tilt＿depth和V2D＿depth方法進行對比處理，結(jié)果如圖3所示.可以看出，Salem提出的Tilt＿depth方法是以垂向一階導(dǎo)數(shù)零值線識別邊界，與真實邊界存在一定誤差.此外，它所反演的深度結(jié)果受旁側(cè)異常影響很大，疊加異?？赡軐?dǎo)致多個不同的場源深度，影響結(jié)果的可靠性.例如對地質(zhì)體B（如圖3a所示），其±45°等值線間距（圖中黑色箭頭線）之一半（即場源埋深）最大值2.25km、最小值約0.79km，嚴重偏離了真實深度（2.5km）.

圖3b顯示的是V2D＿Depth非參數(shù)反演結(jié)果，較之Tilt＿depth有顯著改善.不僅識別的邊界信息更清晰，所反演的深度也更接近真實深度（表1）.其特征等值線之間距（即場源埋深）相對均衡，避免了Tilt＿depth方法受疊加異常的影響在場源四周的反演深度不同（圖3a），表明本方法受旁側(cè)疊加異常的影響小.

為了驗證方法的適應(yīng)性，我們設(shè)計了向下有限延深、傾斜磁化情況時的模型，如圖4所示.其中A、B、C三個棱柱體下延深度分別為2km、2km、4km，其上頂深度分別為1km，0.5km，1.5km.在磁化傾角45°，偏角0°時正演得到的ΔT磁異常如圖4所示.在利用本文方法處理前，將該ΔT磁異常作化極處理，同時為了考慮計算穩(wěn)定性的問題，對該異常加入了0.5%的隨機干擾，處理結(jié)果見圖5、6.

高階導(dǎo)數(shù)的計算對高頻干擾非常敏感，詳見文獻[18］的討論.在本文的處理中，采用向上延拓的方式來消除高頻振蕩的計算問題：即在延拓后的異常基礎(chǔ)上計算的反演深度減去相應(yīng)的延拓高度.圖5表示利用原平面化極磁異常與5個不同高度的延拓異常的深度反演結(jié)果，對應(yīng)的V2D＿depth計算結(jié)果見圖6.

圖3 正演磁異常非參數(shù)反演結(jié)果（a）Tilt＿depth；（b）V2D＿depth；黑色實線表示場源邊界，黑色虛線表示計算的邊界.Fig.3 Non-parameter inversion results via Tilt＿depth and V2D＿depth corresponding to figure 2

表2 模型反演結(jié)果對比分析Table 2 Comparison of the inversion results corresponding to different methods

圖6a中，受高頻干擾的影響，由原平面磁異常計算的V2D＿depth結(jié)果出現(xiàn)振蕩現(xiàn)象，經(jīng)過不同高度的延拓計算之后，得到的V2D＿depth結(jié)果（圖6b—f）消除了噪聲干擾的影響，突出了異常邊界與上頂深度，從反演的結(jié)果看（圖5），運用5個高度的向上延拓結(jié)果計算的深度大致反映出了三個場源的深度.

為了討論非理想模型條件下本文方法的普適性，設(shè)計了如圖7所示的上頂面傾斜的直立棱柱體模型.地質(zhì)體橫縱長度為20km，其邊界在地面的投影如圖中黑色線所示.下底深度50km，上頂是一由北向南傾斜的鍥型傾斜面，其北側(cè)深度5km，南側(cè)深度10km.垂直磁化條件下正演磁異常如圖7a.圖7b是對應(yīng)的非參數(shù)反演結(jié)果，其中h1和h2分別表示采用±arctan（8／17）等值線間距計算的深度：分別為6km和12km.該結(jié)果表明V2D＿depth方法可以突出上頂面傾斜時的傾斜面傾向.

4 實例應(yīng)用

項目主要目的是了解工區(qū)內(nèi)的火山巖分布，偵查可能存在的石炭系火山巖圈閉，了解石炭系地層分布情況.圖8是工區(qū)構(gòu)造分區(qū)圖，可以看出，工區(qū)跨兩個一級構(gòu)造單元：烏倫古坳陷和陸梁隆起，四個二級構(gòu)造單元：青格里底山凸起、索索泉凹陷、三個泉凸起以及陸南凸起.本區(qū)有2個主要磁性層：一個是分布在石炭系中的火成巖所形成的局部磁性層；另一個是深部古老結(jié)晶基底的磁性層.從磁異常特征分析，結(jié)晶基底引起區(qū)域強磁異常，而石炭系火成巖引起局部弱磁異常.

該區(qū)布設(shè)了1∶5萬磁法面積測量工作，ΔT化極磁異常如圖9所示.可以看出，工區(qū)磁異常總體特征為北西走向，總體上反映了海西期褶皺基底的起伏，其中東北、西南角異常高，中部存在北西向的低緩異常帶.除南部有一些局部異常外，總體形態(tài)十分簡單，幅值大約在-650～450nT變化.由于石炭系地層埋深較大，火成巖引起的局部異常為弱異常，其幅值只有十幾至幾十納特；而前寒武強磁性結(jié)晶基底，它引起的區(qū)域異?？蛇_上百納特，強大的區(qū)域背景場掩蓋了石炭系火成巖產(chǎn)生的弱磁異常.

根據(jù)地質(zhì)認識，索索泉凹陷南部Ⅱ2構(gòu)造單元是石炭系地層火成巖局部構(gòu)造發(fā)育區(qū).但受埋藏深度大、加之較強的背景場干擾等影響，并沒有引起可以識別的磁異常.圖10a是采用Tilt＿depth方法獲得局部異常，一定程度上突出了弱異常，但是局部場源細節(jié)模糊，尤其是受疊加異常的影響，Tilt＿depth方法難以得到滿意的深度反演結(jié)果，其±45°等值線并沒有突出相應(yīng)的信息，零值線所反應(yīng)的異常邊界也相對粗糙.

圖10b是采用本文方法處理的結(jié)果：為了消除高頻干擾，我們對磁異常上延了1km高度.該結(jié)果突出了局部異常信息：其中測區(qū)北東角（F1斷裂以北）以及測區(qū)西南部（F4、F2斷裂以南）的高磁異常可能反應(yīng)了基底隆起的影響.從圖9可以看出，對本研究區(qū)有重要意義的Ⅱ2構(gòu)造單元并沒有明顯的磁異常顯示，經(jīng)過處理后（圖10b），原本由火山巖引起的微弱磁異常得到體現(xiàn)，如圖中的F5斷裂與F6斷裂所挾持的部位.另外，磁異常V2D＿depth非參數(shù)反演也得出了深度信息，例如圖10b中A異常反演的深度為1990m，B異常反演的深度為2267m（如圖11所示），該反演深度與地震勘探、鉆孔資料揭露的石炭系地層上頂深度相符合1）中石化西部新區(qū)勘探指揮部，準噶爾盆地青格里底山區(qū)塊重、磁力、電法MT勘探重磁項目成果報告，2008.：鉆井Zk01在井深1981～2894m鉆遇石炭系火山沉積相凝灰?guī)r，反應(yīng)了處理結(jié)果的可靠性.因為鉆孔Zk01鉆遇的沉積相凝灰?guī)r相對其它火成巖具有低磁性特征，因此在圖10b中沒有顯示局部強磁異常.

圖10 工區(qū)磁異常非參數(shù)反演結(jié)果（a）Tilt＿depth方法結(jié)果；（b）V2D＿depth方法結(jié)果，圖中紅色線表示斷裂，黑色線表示零值線.Fig.10 Result of non-parametric fast inversion corresponding to Figure 9

圖11 V2D＿depth方法反演異常深度（圖10b局部放大）Fig.11 Depth inversion of V2D＿depth

對圖中的 A異常或B異常，±arctan（8／17）（≈±25.2°）兩條等值線間距隨邊界位置變化，則由此間距確定的局部場源深度也隨之變化，如A異常的北邊界（或B異常的東邊界）處確定的局部場源深度小得多.究其原因，主要是實際問題中的局部場源頂面一般非水平且非規(guī)則曲面（對比圖7上頂為傾斜面的模型示例）.此外，旁側(cè)異常的疊加也可能是引起這種場源周圍特征等值線間距不均勻的原因.

5 分析與結(jié)論

在平面磁異常處理中，大多傳統(tǒng)方法處理結(jié)果單一，并且對疊加異常計算效果欠佳，尤其是難以獲取受背景場覆蓋的局部弱異常信息.本文根據(jù)“化極磁異常在垂直物性邊界正上方具有梯度極大”的特點進行異常識別，在此基礎(chǔ)上根據(jù)Nabighian提出的通用磁場梯度公式，推導(dǎo)了基于磁異常高階導(dǎo)數(shù)的非參數(shù)反演算法，可以同時獲得場源邊界以及深度反演結(jié)果.為了避免導(dǎo)數(shù)計算引起的計算振蕩，建議采用向上延拓的處理方法，消除高頻振蕩.理論模型與實際應(yīng)用表明，本文方法可以獲得較好的結(jié)果，同時通過一組上頂面非水平的模型試驗說明了該方法可以反映傾斜面的不同深度，以此指導(dǎo)實際應(yīng)用中模型上頂面傾向的判斷.較之Tilt＿depth方法，本文方法具有以下特點：

（1）V2D＿depth可以同時獲得場源邊界與埋深信息，且能夠?qū)崿F(xiàn)微弱異常的分析.

（2）相比Tilt＿depth方法以垂向一階導(dǎo)數(shù)零值點確定場源邊界，V2D＿depth方法能獲得更好的分析效果.

（3）Tilt＿depth方法的應(yīng)用瓶頸是它對復(fù)雜的疊加異常分析效果欠佳：疊加異常導(dǎo)致Tilt＿depth在場源四周的特征等值線間距存在很大差別，這意味著深度參數(shù)計算的不確定性.而V2D＿depth方法在一定程度上克服了這種缺點，模型試驗及實際數(shù)據(jù)處理表明本文方法受疊加異常的影響小，并且獲得的場源邊界及深度參數(shù)都更接近于客觀情況.

（4）關(guān)于數(shù)值奇異的問題.Tilt方法在異常極大值處水平導(dǎo)數(shù)值為零，即分母為零計算產(chǎn)生數(shù)值奇異；而本文方法，如公式（9），產(chǎn)生奇異的前提是垂向二階導(dǎo)數(shù)與垂向?qū)?shù)的水平導(dǎo)數(shù)同時為零，這種情況在異常區(qū)是不可能的.此外，即使出現(xiàn)分母為零的情況，因為分子部分也是零，所以比值結(jié)果為1.分析之，即本文方法不會像Tilt＿depth方法出現(xiàn)數(shù)值奇異的問題.

通過以上分析，本文提出的磁異常非參數(shù)快速反演方法，計算簡單實用，對不同強度的異常、疊加異常都能夠獲得理想的分析結(jié)果，并通過理論模型及實踐應(yīng)用證明了方法的有效性.

致 謝 本文得到了評審專家中肯的建議，并為論文的完善提供了大量有價值的修改意見，筆者向兩位匿名評審專家致以謝意！

（References）

[1]管志寧.地磁場與磁力勘探.北京：地質(zhì)出版社，2005.Guan Z N.Geomagnetic Field and Magnetic Exploration（in Chinese）.Beijing：Geological Publishing House，2005.

[2]Salem A，Williams S，F(xiàn)airhead J，et al.Tilt-depth method：A simple depth estimation method using first-order magnetic derivatives.The Leading Edge，2007，26（12）：1502-1505.

[3]Miller H G，Singh V.Potential filed tilt—a new concept for location of potential field sources.Journal of Applied Geophysics，1994，32（2-3）：213-217.

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

[5]Lee M，Morris W，Ugalde H.Amplitude effects of magnetic signals on depth estimations.11th SAGA biennial technical meeting and exhibition，Swaziland，2009，16-18，276-279.

[6]Lahti I，Karinen T.Tilt derivative multiscale edges of magnetic data.The Leading Edge，2010，29（1）：24-29.

[7]Fairhead J D，Salem A，Cascone L，et al.New developments of the magnetic tilt-depth method to improve structural mapping of sedimentary basins.Geophysical Prospecting，2011，59（6）：1072-1086.

[8]Oru?B，Selim H H.Interpretation of magnetic data in the Sinop area of Mid Black Sea，Turkey，using tilt derivative，Euler deconvolution，and discrete wavelet transform.Journal of Applied Geophysics，2011，74（4）：194-204.

[9]Valléel M A，Smith R S，Keating P.Metalliferous mining geophysics—State of the art after a decade in the new millennium.Geophysics，2011，76（4）：W31-W50.

[10]White J C，Beamish D.Magnetic structural information obtained from the HiRES airborne survey of the Isle of Wight.Proceedings of the Geologists′Association，2011，122（5）：781-786.

[11]Smith R S，Thurston J B，Salem A，et al.A grid implementation of the SLUTH algorithm for visualising the depth and structural index of magnetic sources.Computers &Geosciences，2012，44：100-108.

[12]鄭偉軍.Tilt導(dǎo)數(shù)方法研究及其在重磁數(shù)據(jù)處理中的應(yīng)用[碩士論文］.北京：中國地質(zhì)大學，2010.Zheng W J.Research on Tilt derivative and application in processing gravity and magnetic data [Master thesis］（in Chinese）.Beijing：China University of Geosciences，2010.

[13]王想，李桐林.Tilt梯度及其水平導(dǎo)數(shù)提取重磁源邊界位置.地球物理學進展，2004，9（3）：625-630.Wang X，Li T L.Locating the boundaries of magnetic or gravity sources with tdrand tdr-thdrmethods.Progress in Geophysics（in Chinese），2004，9（3）：625-630.

[14]郭華，熊盛青，于長春.Tilt-depth方法的應(yīng)用.物探與化探，2009，33（5）：583-586.Guo H，Qing S X，Yu C C.The application of tilt-depth technology.Geophysical and Geochemical Exploration（in Chinese），2009，33（5）：583-586.

[15]郭華，于長春，吳燕岡.改進的斜導(dǎo)數(shù)方法及應(yīng)用.物探與化探，2009，33（2）：212-216.Guo H，Yu C C，Wu Y G.An improved tilt derivative method and its application.Geophysical and Geochemical Exploration （in Chinese），2009，33（2）：212-216.

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

[17]張恒磊.現(xiàn)代信號處理的位場濾波方法[博士論文］.武漢：中國地質(zhì)大學，2011.Zhang H L.Potential field data filtering via modern signal processing[Ph.D.thesis］（in Chinese）.Wuhan：China University of Geosciences，2011.

[18]張恒磊，劉天佑，楊宇山.各向異性標準化方差計算重磁源邊界.地球物理學報，2011，54（7）：1921-1927.Zhang H L，Liu T Y，Yang Y S.Calculation of gravity and magnetic source boundary based on anisotropy normalized variance.Chinese J.Geophys.（in Chinese），2011，54（7）：1921-1927.