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

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

0 引言

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

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

1 基本原理

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

(1)

式中：|r-r0|是觀測點(diǎn)到場源點(diǎn)的距離。

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

(2)

式中：μ0為真空磁導(dǎo)率；Bx、By、Bz是磁場B的x、y、z這3個(gè)方向分量。對梯度張量Γ進(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)分別是觀測點(diǎn)與場源點(diǎn)坐標(biāo)；N為構(gòu)造指數(shù)，與磁源幾何形狀有關(guān)，當(dāng)N=0，1，2，3時(shí)，分別對應(yīng)臺(tái)階、巖脈、圓柱體及球體。

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

(6)

(7)

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

(8)

(9)

(10)

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

(11)

則磁源深度z0可表示為：

(12)

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

(13)

顯然，估計(jì)場源深度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倍點(diǎn)距，否則式(6)、(7)不滿足，反演誤差較大，本文取0.1倍點(diǎn)距；z1、z2的大小需要根據(jù)磁異常的信噪比決定，當(dāng)噪聲含量較小時(shí)，z1、z2可較小。另外，當(dāng)異常較為復(fù)雜時(shí)，|z2-z1|≤5Δx且z1和z2兩者取值均不易過大，以避免延拓過大導(dǎo)致不同場源異常相互影響而使反演結(jié)果誤差增大。需要指出的是，由于反演結(jié)果僅在場源位置(x0,y0)有效，而遠(yuǎn)離場源位置的反演結(jié)果完全無參考意義，因此在反演圖中可僅保留歸一化磁源強(qiáng)度垂向差分極大值附近的結(jié)果。

2 模型試驗(yàn)

2.1 單一球體

為了驗(yàn)證本文算法的正確性，設(shè)計(jì)了一個(gè)單一球體模型，模型參數(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 組合模型

為了驗(yàn)證本文方法對復(fù)雜情況的處理能力，設(shè)計(jì)了一個(gè)球體、長方體、薄板、巖脈構(gòu)成的組合模型，各個(gè)模型體參數(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é)果(延拓點(diǎn)距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)雖然存在明顯的異常波動(dòng)，但仍可以清晰地展示出所有模型體的水平位置；垂向差分(圖5b)受噪聲干擾影響明顯，僅能識別出薄板的邊界位置，對其他模型反映極為模糊；反演得到的深度(圖5c)和構(gòu)造指數(shù)(圖5d)也僅在薄板邊界上有異常展示，但反演結(jié)果波動(dòng)大，與真值也存在較大偏差。

從圖6中可以看出，向上延拓3倍點(diǎn)距的歸一化磁源強(qiáng)度(圖6a)異常較為圓滑，垂向差分(圖6b)的穩(wěn)定性獲得較大程度上的提升，能夠較好地展示所有模型體的水平位置。圖6c和6d是場源深度與構(gòu)造指數(shù)反演結(jié)果，表2也給出了無噪聲和含噪聲不同延拓高度時(shí)的反演結(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。對比無噪聲時(shí)的反演結(jié)果可知，所有模型體的深度反演精度都有所降低，尤其棱柱體深度的相對誤差高達(dá)48%，不過其他模型體的深度反演誤差均小于10%，且深度反演的平均相對誤差為 15.4%，因此，認(rèn)為反演結(jié)果仍具有一定的可靠性。

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

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

從圖7可以看出，向上延拓5倍點(diǎn)距的歸一化磁源強(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時(shí)的反演結(jié)果可知，幾乎所有的模型體深度反演精度進(jìn)一步下降，尤其棱柱體深度的相對誤差由48%升至88%，深度反演的平均相對誤差為27.4%。由此可見，延拓高度較大時(shí)，雖然歸一化磁源強(qiáng)度及其垂向差分異常更加光滑，但是反演精度卻有所下降，主要在于延拓高度較大時(shí)，不同場源的異常相互影響更嚴(yán)重。

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

3 實(shí)例應(yīng)用

為了驗(yàn)證本文方法的實(shí)用性，選取了內(nèi)蒙古M地區(qū)的地面磁異常進(jìn)行實(shí)驗(yà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ū)東部還存在一個(gè)明顯的高磁異常圈閉區(qū)，地表出露的閃長巖上方基本都是以高磁異常顯示。圖9b是z1=-0.5 km時(shí)的歸一化磁源強(qiáng)度，可以看出，歸一化磁源強(qiáng)度異常形態(tài)更加簡潔，高異常分布區(qū)與磁異常中的高異常區(qū)基本吻合，其主要原因是該地區(qū)處于高緯度地區(qū)，高磁異常基本與磁性體位置相對應(yīng)。圖9c是z1=-0.5 km時(shí)的歸一化磁源強(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ì)算，并根據(jù)不同高度上差分關(guān)系，提出了一種磁源參數(shù)快速反演方法。模型分析表明了歸一化磁源強(qiáng)度垂向差分提高了原方法的分辨率，能更清晰地展示異常間關(guān)系；基于歸一化磁源強(qiáng)度垂向差分的反演方法能夠有效地反演出磁源的深度與構(gòu)造指數(shù)。在內(nèi)蒙古M區(qū)實(shí)例應(yīng)用中，歸一化磁源強(qiáng)度垂向差分識別出了一條近東西走向的高磁異常帶和一個(gè)范圍較大的高磁異常圈閉區(qū)，反演算法計(jì)算結(jié)果表明了東西走向的高磁異常帶對應(yīng)的磁源埋深較淺，形狀偏向于巖脈，而高磁異常圈閉區(qū)對應(yīng)的磁源埋深較大。本文提出的磁異?？焖俜囱莘椒軌蚝芎玫毓烙?jì)場源位置信息與幾何形狀信息，方法原理簡單，易于實(shí)現(xiàn)，為三維斜磁化磁異常解釋提供了新的研究思路。