基于歸一化磁源強度垂向差分的磁源參數快速估計方法

黃遠生，王彥國，羅瀟

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

0 引言

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

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

1 基本原理

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

(1)

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

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

(2)

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

(3)

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

(4)

Beiki等[26]給出了歸一化磁源強度的統一表達式：

(5)

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

當Δz→0，對于兩個不同高度z1,z2，有：

(6)

(7)

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

(8)

(9)

(10)

當x→x0,y→y0時：

(11)

則磁源深度z0可表示為：

(12)

將z0表達式帶入式(8)或(9)，可求得構造指數N：

(13)

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

2 模型試驗

2.1 單一球體

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

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

2.2 組合模型

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

圖2 組合模型正演磁異常

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

表1 模型體參數

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

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

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

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

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

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

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

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

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

3 實例應用

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

圖8 內蒙古M區地質

圖9 內蒙古M地區磁數據處理(z1=-0.5 km，z2=-0.75 km)

4 結論

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