黃遠(yuǎn)生，王彥國，羅瀟

(1.東華理工大學(xué) 地球物理與測控技術(shù)學(xué)院,江西 南昌 330013；2.海南水文地質(zhì)工程地質(zhì)勘察院，海南 ?？?571100；3.核工業(yè)230研究所,湖南 長沙 410011)

0 引言

磁異常快速反演一直是磁法數(shù)據(jù)處理與解釋的研究重點與熱點，該類方法能夠在無先驗信息約束下快速獲取場源的位置與幾何參數(shù)。目前較為常用的快速反演方法有歐拉反褶積法、解析信號法、tilt-depth法等。歐拉反褶積是Peters[1]提出的，Thompson[2]推導(dǎo)了二維歐拉反褶積，Reid 等[3]將其推廣至三維。張量歐拉反褶積[4]擴(kuò)展了歐拉方程個數(shù)，提高了反演解收斂性；Huang 等[5]證明了位場解析信號同樣滿足歐拉齊次方程；AN-EUL 法[6]能夠快速估算場源深度與構(gòu)造指數(shù)；Tilt-Euler 法[7]無需已知場源構(gòu)造指數(shù)，避免了因構(gòu)造指數(shù)選取不當(dāng)導(dǎo)致反演解發(fā)散的問題；張量局部波數(shù)的歐拉反褶積法[8]進(jìn)一步擴(kuò)展了方程組個數(shù)，獲得了良好的應(yīng)用效果。不過歐拉反褶積及其改進(jìn)方法的反演解中存在大量的虛假解，需要建立有效的篩選機(jī)制，另外背景場、噪聲干擾及窗口大小選擇也會對反演結(jié)果產(chǎn)生一定影響。二維磁解析信號振幅不受磁化方向影響[9]，Macleod[10]推導(dǎo)出了二維磁解析信號振幅的通用表達(dá)式；Salem 等[11]在解析信號及其水平總梯度基礎(chǔ)上采用線性最小二乘法計算磁源深度及構(gòu)造指數(shù)；Salem 等[12]又根據(jù)解析信號及磁異常垂向?qū)?shù)解析信號關(guān)系進(jìn)行磁源參數(shù)估計；Ma 和Du[13]在解析信號振幅的解析信號與解析信號比值基礎(chǔ)上提出了改進(jìn)解析信號來實現(xiàn)單一場源深度及構(gòu)造指數(shù)的估算；Cooper[14-15]同樣采用不同階次解析信號來估計巖脈及臺階的深度；Cooper 和Whitehead[16]利用不同階次解析信號比值進(jìn)行場源深度估算，不過改進(jìn)方法不需要已知構(gòu)造指數(shù)；Cooper[17]又采用了解析信號對數(shù)實現(xiàn)磁源深度估計工作。Wang等[18]采用解析信號及其倒數(shù)進(jìn)行磁源參數(shù)反演計算，對深部場源反演效果較好。但三維解析信號在一定程度上受磁化角度影響，在三維磁異常反演中使用較小。Tilt-depth 法可用于快速估算場源上頂深度，該方法是由Salem 等[19]在Tilt 梯度及磁場通用梯度公式進(jìn)行理論推導(dǎo)而來；Fairhead 等[20]提出了基于化極與化赤相結(jié)合的tilt-depth 法，還可以實現(xiàn)場源磁化率估計；張恒磊等[21]提出了基于二階導(dǎo)數(shù)的磁源邊界與頂部深度快速反演方法，有效地消除了區(qū)域場對反演結(jié)果的影響；Wang 等[22]提出了改進(jìn)tilt-depth 法來估計磁源的上頂與下底深度，并采用多特征點聯(lián)合計算來提高反演解的可靠性；Cooper[23]在Hilbert 變換基礎(chǔ)上推導(dǎo)出了巖脈模型的垂直磁化磁位的tilt-depth 公式，提高了方法實用性及穩(wěn)定性；曹偉平等[24]對tilt-depth 法進(jìn)行了深入研究，指出了該方法并不適用于埋深大、水平尺度小的磁源深度反演計算。由于tilt-depth 法及改進(jìn)算法均需要在化極或化赤基礎(chǔ)上完成計算，因此并不能直接用于斜磁化磁源深度的估算。

歸一化磁源強(qiáng)度是在磁梯度張量特征值基礎(chǔ)上提出的，該方法不(少)受磁化方向影響，是三維磁異常解釋的常用工具。Wilson[25]首先給出了磁偶極子下歸一化磁源強(qiáng)度的表達(dá)式，Beiki等[26]給出了歸一化磁源強(qiáng)度的通用表達(dá)形式，建立了歸一化磁源強(qiáng)度與場源位置、構(gòu)造指數(shù)的關(guān)系，同時結(jié)合歐拉反褶積實現(xiàn)了三維磁源的位置及構(gòu)造指數(shù)反演；Pikington 和Beiki[27]在歸一化磁源強(qiáng)度基礎(chǔ)上實現(xiàn)了剩磁情況下的磁源空間位置及磁化率反演；Guo等[28]將歸一化磁源強(qiáng)度與相關(guān)成像結(jié)合來反映地下磁性體的三維空間展布狀態(tài)；饒椿鋒等[29]將歸一化磁源強(qiáng)度與正則化共軛梯度法結(jié)合來反演磁源的磁化率。本文在歸一化磁源強(qiáng)度垂向差分基礎(chǔ)上，提出了一種場源快速反演方法，并進(jìn)行了模型分析與實例應(yīng)用。

1 基本原理

對于體積為V，磁化強(qiáng)度為M的磁性體，其磁位U的表達(dá)式[30]為：

(1)

式中：|r-r0|是觀測點到場源點的距離。

磁梯度張量[31]可以表示為：

(2)

式中：μ0為真空磁導(dǎo)率；Bx、By、Bz是磁場B的x、y、z這3個方向分量。對梯度張量Γ進(jìn)行對角化處理，得：

(3)

式中：l1、l2和l3為Γ的特征向量；λ1、λ2、λ3為對應(yīng)的特征值。Wilson[25]給出磁偶極子的歸一化磁源強(qiáng)度μ的定義，并推導(dǎo)出梯度張量特征值表達(dá)式：

(4)

Beiki等[26]給出了歸一化磁源強(qiáng)度的統(tǒng)一表達(dá)式：

(5)

式中：α是與磁化強(qiáng)度有關(guān)的物理量；(x,y,z)和(x0,y0,z0)分別是觀測點與場源點坐標(biāo)；N為構(gòu)造指數(shù)，與磁源幾何形狀有關(guān)，當(dāng)N=0，1，2，3時，分別對應(yīng)臺階、巖脈、圓柱體及球體。

當(dāng)Δz→0，對于兩個不同高度z1,z2，有：

(6)

(7)

分別將式(6)、(7)對式(5)做比值，得：

(8)

(9)

(10)

當(dāng)x→x0,y→y0時：

(11)

則磁源深度z0可表示為：

(12)

將z0表達(dá)式帶入式(8)或(9)，可求得構(gòu)造指數(shù)N：

(13)

顯然，估計場源深度z0和構(gòu)造指數(shù)N，需要事先給定場源的水平位置(x0,y0)，Δz取值以及z1、z2的大小。由于歸一化磁源強(qiáng)度的垂向差分極大值對應(yīng)著場源水平位置，因此場源的水平坐標(biāo)(x0,y0)則看作是歸一化磁源強(qiáng)度垂向差分(Δμ1)極大值位置；Δz不易過大，一般取0.05～0.2倍點距，否則式(6)、(7)不滿足，反演誤差較大，本文取0.1倍點距；z1、z2的大小需要根據(jù)磁異常的信噪比決定，當(dāng)噪聲含量較小時，z1、z2可較小。另外，當(dāng)異常較為復(fù)雜時，|z2-z1|≤5Δx且z1和z2兩者取值均不易過大，以避免延拓過大導(dǎo)致不同場源異常相互影響而使反演結(jié)果誤差增大。需要指出的是，由于反演結(jié)果僅在場源位置(x0,y0)有效，而遠(yuǎn)離場源位置的反演結(jié)果完全無參考意義，因此在反演圖中可僅保留歸一化磁源強(qiáng)度垂向差分極大值附近的結(jié)果。

2 模型試驗

2.1 單一球體

為了驗證本文算法的正確性，設(shè)計了一個單一球體模型，模型參數(shù)為：半徑為0.5 km，埋深為1 km，磁化傾角和磁化偏角均為60°，磁化強(qiáng)度為1 A/m。圖1a為單一球體模型的磁異常，圖1b與1c是歸一化磁源強(qiáng)度及其垂向差分結(jié)果，可以看出，歸一化磁源強(qiáng)度(圖1b)極大值與球體質(zhì)心相吻合，等值線成同心圓分布，異常特征完全不受磁化方向影響，而歸一化磁源強(qiáng)度垂向差分(圖1c)相對于歸一化磁源強(qiáng)度，異常更加匯聚。圖1d和1e是本文算法獲得的磁源深度及構(gòu)造指數(shù)，可以看出，深度及構(gòu)造指數(shù)反演圖均在球體質(zhì)心位置存在極小值，極小值分別為1.02 km及3.1，均與理論值1 km和3非常接近，這反映了本文算法的正確性。

圖1 單一球體磁異常反演結(jié)果(z1=0 km, z2=0.2 km)

2.2 組合模型

為了驗證本文方法對復(fù)雜情況的處理能力，設(shè)計了一個球體、長方體、薄板、巖脈構(gòu)成的組合模型，各個模型體參數(shù)見表1。圖2是組合模型產(chǎn)生的磁異常，可以看出，歸一化磁源強(qiáng)度(圖3a)能夠利用極大值很好地反映出球體質(zhì)心、棱柱體及薄板邊界和巖脈位置；歸一化磁源強(qiáng)度垂向差分(圖3b)對各模型體刻畫的更加精細(xì)，異常分辨率獲得了明顯的提高。圖3c、d是本文方法反演得到的場源埋深及構(gòu)造指數(shù)。在球體①質(zhì)心的反演深度為1.52 km，構(gòu)造指數(shù)反演值為3.0；棱柱體②邊界位置上的深度反演解在0.5～0.7 km之間，平均值為0.61 km，構(gòu)造指數(shù)反演解集中在0.2～0.5之間，平均值為0.34；薄板體③邊界上的深度反演值在0.47～0.65 km之間，平均值為0.51 km，構(gòu)造指數(shù)在0.55～1之間，平均值為0.77；巖脈④的反演深度集中在1.0～1.15 km，平均值為1.08 km，構(gòu)造指數(shù)在1.05～1.25之間，平均值為1.16。從反演結(jié)果可以看出，由于異常疊加影響，反演結(jié)果均與理論值存在偏差，但較接近于理論值，且深度反演的平均相對誤差僅為8.3%，反演結(jié)果仍較可靠。

圖2 組合模型正演磁異常

圖3 組合模型磁異常的反演結(jié)果(延拓點距z1=0 km，z2=-0.2 km)

表1 模型體參數(shù)

為了了解噪聲干擾對本文算法的影響情況，對疊加異常(圖2)添加了1%的隨機(jī)噪聲，見圖4。圖5、6和7分別是選取了不同延拓高度的歸一化磁源強(qiáng)度及其垂向差分，和深度及構(gòu)造指數(shù)反演結(jié)果。

圖4 含1%隨機(jī)噪聲的組合模型磁異常

圖5 含噪組合模型磁異常的反演結(jié)果(z1=0 km, z2=-0.2 km)

從圖5可以看出，未延拓情況下的歸一化磁源強(qiáng)度(圖5a)雖然存在明顯的異常波動，但仍可以清晰地展示出所有模型體的水平位置；垂向差分(圖5b)受噪聲干擾影響明顯，僅能識別出薄板的邊界位置，對其他模型反映極為模糊；反演得到的深度(圖5c)和構(gòu)造指數(shù)(圖5d)也僅在薄板邊界上有異常展示，但反演結(jié)果波動大，與真值也存在較大偏差。

從圖6中可以看出，向上延拓3倍點距的歸一化磁源強(qiáng)度(圖6a)異常較為圓滑，垂向差分(圖6b)的穩(wěn)定性獲得較大程度上的提升，能夠較好地展示所有模型體的水平位置。圖6c和6d是場源深度與構(gòu)造指數(shù)反演結(jié)果，表2也給出了無噪聲和含噪聲不同延拓高度時的反演結(jié)果。在該延拓情況下，球體質(zhì)心埋深反演結(jié)果為1.46 km，構(gòu)造指數(shù)為2.88；棱柱體深度反演解在0.4～1.1 km之間，均值為0.74 km，構(gòu)造指數(shù)反演解在0～1.2之間，均值為0.63；薄板的深度反演值在0.4～0.8 km之間，均值為0.51 km，構(gòu)造指數(shù)反演值在0.4～1.5之間，均值為0.76；巖脈深度反演解在0.8～1.4 km之間，均值為1.09 km，構(gòu)造指數(shù)反演值在0.7～1.6之間，均值為1.18。對比無噪聲時的反演結(jié)果可知，所有模型體的深度反演精度都有所降低，尤其棱柱體深度的相對誤差高達(dá)48%，不過其他模型體的深度反演誤差均小于10%，且深度反演的平均相對誤差為 15.4%，因此，認(rèn)為反演結(jié)果仍具有一定的可靠性。

表2 無噪及含噪組合磁異常不同高度上的反演結(jié)果

圖6 含噪組合模型磁異常的反演結(jié)果(z1=-0.3 km, z2=-0.5 km)

從圖7可以看出，向上延拓5倍點距的歸一化磁源強(qiáng)度(圖7a)及其垂向差分(圖7b)異常更加圓滑。圖7c和7d是場源深度與構(gòu)造指數(shù)反演結(jié)果，球體質(zhì)心埋深反演結(jié)果為1.54 km，構(gòu)造指數(shù)為 3.05；棱柱體深度反演解在0.5～1.4 km之間，均值為0.94 km，構(gòu)造指數(shù)反演解在0.3～1.5之間，均值為0.80；薄板的深度反演值在0.3～0.7 km之間，均值為0.45 km，構(gòu)造指數(shù)反演值在0.3～1之間，均值為0.61；巖脈深度反演解在0.8～1.3 km之間，均值為1.09 km，構(gòu)造指數(shù)反演值在0.5～1.5之間，均值為1.18。對比z1=-0.3 km，z2=-0.5 km時的反演結(jié)果可知，幾乎所有的模型體深度反演精度進(jìn)一步下降，尤其棱柱體深度的相對誤差由48%升至88%，深度反演的平均相對誤差為27.4%。由此可見，延拓高度較大時，雖然歸一化磁源強(qiáng)度及其垂向差分異常更加光滑，但是反演精度卻有所下降，主要在于延拓高度較大時，不同場源的異常相互影響更嚴(yán)重。

圖7 含噪組合模型磁異常的反演結(jié)果(z1=-0.5 km, z2=-0.7 km)

3 實例應(yīng)用

為了驗證本文方法的實用性，選取了內(nèi)蒙古M地區(qū)的地面磁異常進(jìn)行實驗。圖8是研究區(qū)的地質(zhì)圖，可以看出，該地區(qū)地表分布的是沉積巖和火山巖，具有一定磁性的晚奧陶世閃長巖主要分布在中西部(東2～10 km，北3～6 km)、中部(東23～28 km，北3～7 km)及中北部(東29～35 km，北5～9 km)3塊，其他類型巖石基本無磁性。圖9a是研究區(qū)地面磁異常，測網(wǎng)密度為250 m(東)×50 m(北)，圖中可以看出，測區(qū)中部存在一條近EW走向轉(zhuǎn)EN走向的條帶狀磁異常，在測區(qū)東部還存在一個明顯的高磁異常圈閉區(qū)，地表出露的閃長巖上方基本都是以高磁異常顯示。圖9b是z1=-0.5 km時的歸一化磁源強(qiáng)度，可以看出，歸一化磁源強(qiáng)度異常形態(tài)更加簡潔，高異常分布區(qū)與磁異常中的高異常區(qū)基本吻合，其主要原因是該地區(qū)處于高緯度地區(qū)，高磁異?；九c磁性體位置相對應(yīng)。圖9c是z1=-0.5 km時的歸一化磁源強(qiáng)度垂向差分結(jié)果，可以看出，該圖與歸一化磁源強(qiáng)度基本一致，不過異常更加收斂，另外東側(cè)的高異常分布區(qū)較模糊，這可能反映了該區(qū)磁源埋深較大。圖9d和9e是利用本文算法得到的磁源深度z0及構(gòu)造指數(shù)N分布圖，可以看出，地表出露的閃長巖分布區(qū)，反演深度基本都是0 km，構(gòu)造指數(shù)多分布在1附近，即可當(dāng)作巖脈看待；中部條帶狀磁異常上的深度反演結(jié)果基本都在0～1 km之間，即隱伏的高磁性體(推測為閃長巖)埋深相對較淺，構(gòu)造指數(shù)反演結(jié)果則由西往東有逐步增加的趨勢，基本從1變化到3，這可能是西側(cè)磁性體寬度小而東側(cè)偏大，或東側(cè)磁源從淺往深規(guī)模逐漸擴(kuò)大所致；東部高磁異常圈閉區(qū)的深度反演結(jié)果表明該磁性體從西南往東北埋深逐漸減小，而構(gòu)造指數(shù)從西南往東北同樣逐漸減小，這可能反映了磁性體形態(tài)上從穹窿狀往脈狀的一種轉(zhuǎn)變。

圖8 內(nèi)蒙古M區(qū)地質(zhì)

圖9 內(nèi)蒙古M地區(qū)磁數(shù)據(jù)處理(z1=-0.5 km，z2=-0.75 km)

4 結(jié)論

歸一化磁源強(qiáng)度是三維磁數(shù)據(jù)解釋常用的方法，主要在于該方法不(少)受磁化方向影響。本文在歸一化磁源強(qiáng)度基礎(chǔ)上，進(jìn)行了垂向差分計算，并根據(jù)不同高度上差分關(guān)系，提出了一種磁源參數(shù)快速反演方法。模型分析表明了歸一化磁源強(qiáng)度垂向差分提高了原方法的分辨率，能更清晰地展示異常間關(guān)系；基于歸一化磁源強(qiáng)度垂向差分的反演方法能夠有效地反演出磁源的深度與構(gòu)造指數(shù)。在內(nèi)蒙古M區(qū)實例應(yīng)用中，歸一化磁源強(qiáng)度垂向差分識別出了一條近東西走向的高磁異常帶和一個范圍較大的高磁異常圈閉區(qū)，反演算法計算結(jié)果表明了東西走向的高磁異常帶對應(yīng)的磁源埋深較淺，形狀偏向于巖脈，而高磁異常圈閉區(qū)對應(yīng)的磁源埋深較大。本文提出的磁異?？焖俜囱莘椒軌蚝芎玫毓烙媹鲈次恢眯畔⑴c幾何形狀信息，方法原理簡單，易于實現(xiàn)，為三維斜磁化磁異常解釋提供了新的研究思路。