基于泊松方程的空間波數(shù)混合域二度體磁異常數(shù)值模擬

趙東東，王旭龍，周印明，戴世坤

1.中國地質(zhì)調(diào)查局南京地質(zhì)調(diào)查中心，南京 210016 2.中南大學(xué)地球科學(xué)與信息物理學(xué)院，長沙 410083 3.有色金屬成礦預(yù)測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室(中南大學(xué))，長沙 410083 4.湖南科技大學(xué)信息與電氣工程學(xué)院，湖南 湘潭 411201

0 引言

磁法勘探是地球物理勘探中發(fā)展比較成熟的方法之一，已經(jīng)廣泛應(yīng)用于直接尋找磁鐵礦及其共生礦床、固體礦產(chǎn)、石油天然氣構(gòu)造的普查和不同比例尺的地質(zhì)填圖以及深部、區(qū)域、全球構(gòu)造的研究等[1-3]。在實際野外生產(chǎn)中，通常將走向方向尺度遠比垂直走向方向尺度大的地質(zhì)體近似用沿走向無限延伸的二度體代替。

目前，磁異常數(shù)值模擬方法一般主要分為空間域方法和頻率域方法?？臻g域方法通過異常場解析式直接求取空間任意點的精確異常場，其優(yōu)點是原理簡單、結(jié)果精確，缺點是解析式比較復(fù)雜，推導(dǎo)過程繁瑣，且當(dāng)計算大量位置點異常場數(shù)據(jù)時，速度較慢。有關(guān)二度體的磁異常數(shù)值模擬研究最早就是在空間域開始的，Talwani等[4]首次給出了多邊形截面二度體磁異常解析表達式；在此基礎(chǔ)上，Talwani[5]在計算機上實現(xiàn)了近似二維棱柱狀多邊形數(shù)值模擬；Won等[6]根據(jù)泊松關(guān)系導(dǎo)出了多邊形截面二度體磁異常的解析式，并利用Fortran77實現(xiàn)了相關(guān)算法，優(yōu)勢在于可以模擬地下任意觀測點的磁場；Shuey等[7]、Rasmussen等[8]以及Cady[9]在前人研究基礎(chǔ)上對原來的算法進行校正，使其能模擬實際地質(zhì)體所產(chǎn)生的磁異常；賈真[10]系統(tǒng)地推導(dǎo)了均勻多邊形截面二度體磁異常及磁異常梯度正演計算公式，并重點探討了正演公式中所包含的奇異性；Jeshvaghani等[11]采用自適應(yīng)有限單元法實現(xiàn)了基于拉普拉斯方程的二度體磁異常正演；Kravchinsky等[12]在Talwani[5]解析公式的基礎(chǔ)上編寫了一套基于Matlab的二度體磁異常模擬軟件，并提高了數(shù)值精度。

頻率域方法根據(jù)場源產(chǎn)生異常場的傅里葉變換解析表達式，通過數(shù)值方法計算該異常頻譜的反傅里葉變換得到異常場。相對于空間域方法而言，其優(yōu)點是表達式簡潔，計算效率較高。在頻率域，Bhattacharyya等[13]推導(dǎo)了任意二度體的磁場頻率域表達式，并進行了反演研究；Pedersen[14]給出了任意二度體、二度半體(以三角形組合為其表面的模型)重磁異常頻率域表達式，并指出該公式相對空間域解析公式更為簡潔；吳宣志[15]將任意二維物體磁場的波譜解析表達式推廣至可用于任意變磁化強度模型；柴玉璞[16]將偏移抽樣方法發(fā)展為一整套完整的理論體系，有效提升了傅里葉變換的數(shù)值精度，使得頻率域方法逐漸在處理復(fù)雜模型正演問題中得以廣泛應(yīng)用；吳樂園[17]在偏移抽樣理論的基礎(chǔ)上提出Gauss-FFT(fast Fourier transform)理論方法技術(shù)，有效提升了變換精度，降低了水平方向邊界條件帶來的影響，并將其應(yīng)用于頻率域二度體磁異常正演，取得了不錯的效果。

然而，由于在面向?qū)嶋H應(yīng)用的大規(guī)模、高效、高精度、精細化三維反演成像中，網(wǎng)格單元剖分?jǐn)?shù)目達數(shù)千萬甚至數(shù)億，導(dǎo)致計算量和存儲量巨大，使得常規(guī)從積分方程的空間域和頻率域很難高效、精細模擬任意復(fù)雜模型；故在介于空間域和頻率域之間的空間波數(shù)混合域，從微分方程角度出發(fā)提出一種可將三維或二維問題分解為多個一維小問題的數(shù)值求解方法。該方法已成功應(yīng)用于二度體重力異常和三度體位場正演模擬[18-19]，并且在計算精度和效率方面表現(xiàn)出一定的優(yōu)勢。在此基礎(chǔ)上，本文將基于泊松方程的空間波數(shù)混合域方法[18]用于磁異常二維數(shù)值模擬。該方法通過一維傅里葉變換把二維磁位問題轉(zhuǎn)化為多個一維問題，由此大大減少了計算量，不同波數(shù)之間常微分方程相互獨立，具有高度并行性；保留垂向為空間域，有利于適應(yīng)復(fù)雜地形的模擬，兼顧了計算精度與計算效率；采用有限單元法求解不同波數(shù)滿足的常微分方程，充分利用追趕法求解對角線性方程組的高效性以進一步提高計算效率[20]。在模型算例中，分別設(shè)計突變介質(zhì)模型、連續(xù)介質(zhì)模型和復(fù)雜地形模型用于驗證空間波數(shù)混合域二度體磁異常正演數(shù)值模擬方法的計算精度與計算效率。

1 基本理論

1.1 控制方程

弱磁性體滿足的二維泊松方程[21]為(詳見附錄A)

?2U(x,z)=?·M(x,z)。

(1)

式中：U(x,z)為磁位；M(x,z)為磁化強度矢量。對式(1)沿x方向做一維傅里葉變換，則有

(2)

1.2 邊界條件

在無磁源區(qū)域，空間波數(shù)混合域磁位滿足的常微分方程(式(2))通解為

(3)

式中，A和B為通解的系數(shù)，表示任意常數(shù)。規(guī)定垂向坐標(biāo)軸向下為正，上邊界z=zmin，下邊界z=zmax，如圖1所示。則上、下邊界條件為：

(4)

1.3 磁場和磁梯度張量計算

空間域磁感應(yīng)強度、磁梯度張量與磁位的關(guān)系如下：

圖1 邊界條件示意圖

(5)

(6)

式中：B為空間域磁感應(yīng)強度；T為空間域磁梯度張量。對式(5)和式(6)進行一維傅里葉變換可以得到：

(7)

(8)

1.4 零波數(shù)處理

特別地，當(dāng)k=0時，空間波數(shù)混合域磁位滿足的常微分方程為

(9)

式(9)相當(dāng)于水平層狀磁性介質(zhì)產(chǎn)生的磁場，即空間波數(shù)混合域的水平磁場為零，垂直磁場為磁化強度的垂直分量，具體形式如下：

(10)

(11)

(12)

2 常微分方程求解

式(2)與引力位滿足的空間波數(shù)混合域常微分方程類似[8]，本文采用基于二次插值的一維有限單元法進行數(shù)值模擬。該方法一方面能實現(xiàn)復(fù)雜模型和復(fù)雜地形的磁異常高精度模擬，另一方面利用追趕法求解有限單元法形成的對角線性方程組(五對角)具有計算速度快、占用內(nèi)存小的特征。

聯(lián)立式(2)和式(4)，即是空間波數(shù)混合域二度體磁位滿足的一維邊值問題。基于變分原理構(gòu)造泛函[20]，可得到與邊值問題等價的變分問題(詳見附錄B)：

(13)

(14)

其中：

和列陣Xe、Ze詳見附錄C。由于δuT≠0，故

Ku=P。

(15)

對于該五對角定帶寬線性方程組，采用追趕法實現(xiàn)高效、高精度求解，通過磁場、磁梯度張量與磁位之間的關(guān)系(式(7)和式(8))易得到空間波數(shù)混合域的磁場和磁梯度張量，采用Gauss-FFT對空間波數(shù)混合域磁位、磁場及磁梯度張量進行反變換，得到空間域的磁位、磁場和磁梯度張量[17-18]。

3 數(shù)值試驗

設(shè)計截面為矩形的常磁化率和變磁化率二度體模型[21]，用于驗證本文算法的正確性和可靠性。常磁化率和變磁化率模型空間展布及模型剖分參數(shù)均一致(圖2)，模擬范圍：x方向為-500～500 m，z方向為0～500 m，網(wǎng)格剖分規(guī)模為201×101，水平網(wǎng)格和垂向網(wǎng)格間距均為5 m；異常體x方向延伸范圍為-100～100 m，z方向為150～300 m。本文算法均采用Fortran95語言編寫，模型算例均在配置為CPU-i7-4790、主頻為3.30 GHz、物理內(nèi)存為32.00 GB的工作站上串行測試得到。

3.1 突變介質(zhì)模型

圖2a所示的突變介質(zhì)模型中異常體的磁化率為0.01；給定地球正常磁場強度為45 000 nT，磁傾角為45°，磁偏角為0°。通過對比本文算法結(jié)果與解析結(jié)果[18]驗證該算法的正確性和計算精度。圖3為突變介質(zhì)模型磁位、磁場、磁梯度張量數(shù)值解和解析解以及地面節(jié)點的相對誤差。

從圖3可以看出：突變介質(zhì)模型磁位、磁場和磁梯度張量的數(shù)值解與解析解吻合較好，磁位最大相對誤差絕對值僅為0.42%，磁場分量Bx和Bz最大相對誤差均絕對值為0.48%；磁梯度張量的最大相對誤差均絕對值為0.45%。磁位、磁場以及磁梯度張量的相對誤差絕對值均小于0.50%，僅在零值線附近相對誤差較大，驗證了本文算法的正確性和可靠性，且精度較高。

a. 突變介質(zhì)模型；b. 連續(xù)介質(zhì)模型。κ.磁化率。

3.2 連續(xù)介質(zhì)模型

設(shè)計磁化率隨深度變化的連續(xù)介質(zhì)模型(圖2b)，用于進一步驗證本文算法的適用性和可靠性。給定地球正常磁場強度為45 000 nT，磁傾角為90°，磁偏角為0°。異常體磁化率在深度方向上有二次變化規(guī)律：

κ=-(z-150)(z-300)×0.00001，
z∈[150,300]。

(16)

圖4為磁化率連續(xù)變化的二度體模型數(shù)值解與解析解及相對誤差?？梢钥闯觯哼B續(xù)介質(zhì)模型磁位的數(shù)值解與解析解最大相對誤差絕對值為0.004%，磁場Bx的數(shù)值解與解析解最大相對誤差絕對值小于0.01%，與突變介質(zhì)模型相比精度更高。這是由于連續(xù)介質(zhì)模型在傅里葉變換過程中引入的誤差較小。

左.解析解；中.數(shù)值解；右.相對誤差。

a. U數(shù)值解與解析解；b. Bx數(shù)值解與解析解；c. U相對誤差；d. Bx相對誤差。

3.3 計算效率

在任意復(fù)雜模型剖分滿足模擬精度的前提下，網(wǎng)格剖分規(guī)模越大，計算成本也越高；這是由于算法的計算效率主要由剖分網(wǎng)格規(guī)模決定的。圖5給出了不同網(wǎng)格剖分規(guī)模的模擬時間，可以看出：當(dāng)網(wǎng)格

圖5 不同網(wǎng)格剖分規(guī)模的模擬時間

剖分為1 001×1 001時，模擬時間約為0.38 s，當(dāng)網(wǎng)格剖分規(guī)模為2 501×2 501時，模擬時間為4.18 s；并且隨著網(wǎng)格剖分的指數(shù)增長，計算時間僅呈線性增長，而不是指數(shù)增長?？梢姳疚乃惴ú粌H計算效率較高，而且網(wǎng)格剖分規(guī)模越大，計算效率優(yōu)勢越明顯。

3.4 起伏地形模型

設(shè)計水平地形模型(圖6a)和起伏地形模型(圖6b)，用于進一步驗證本文算法在復(fù)雜起伏地形條件下的適用性和可靠性。模擬范圍：x方向為-250～250 m，z方向為-100～200 m，網(wǎng)格剖分規(guī)模為501×301。在水平地形模型(圖6a)中，異常體截面為矩形，沿x方向延伸范圍為-100～100 m，沿z方向延伸范圍為0～100 m；起伏地形模型(圖6b)由山谷和山脊構(gòu)成，異常范圍與水平模型一致。異常體磁化率均為0.02，地球正常磁場強度為45 000 nT，磁傾角為45°，磁偏角為0°。利用本文算法計算觀測線高度為地面以上100 m的磁場及磁梯度張量，計算時間為0.21 s。

地面以上100 m高度的磁場及磁梯度張量模擬結(jié)果見圖7?？梢钥闯觯涸谏郊固?，磁場Bz及磁梯度分量Tzz和Txz相對水平地形模型變化較快，正異常幅值增加，磁場Bx相對水平地形模型負(fù)異常幅值有所增加；在山谷處，磁場Bz和磁梯度分量Tzz和Txz相對水平地形模型變化減緩，但負(fù)異常幅值略有增強，磁場Bx相對水平地形模型變化明顯減緩，進一步驗證了本文算法在起伏地形條件下的適用性和可靠性。

圖6 水平地形模型(a)和起伏地形模型(b)

a. 水平地形磁場；b. 水平地形磁梯度；c. 起伏地形磁場；d. 起伏地形磁梯度。

4 結(jié)論

1)基于泊松方程的空間波數(shù)混合域數(shù)值模擬，主要利用了沿水平方向傅里葉變換方法和追趕法求解定帶寬線性方程組的獨特優(yōu)勢，為二度體磁異?？焖儆嬎闾峁┝艘环N新的研究思路和方案。

2)突變介質(zhì)模型、連續(xù)介質(zhì)模型和起伏地形模型的模擬結(jié)果表明本文算法數(shù)值模擬精度高、計算速度快，具有一定的實用價值。

3)該方法可能由于傅里葉變換的截斷效應(yīng)或強制周期性在介質(zhì)突變邊界引起較大的數(shù)值誤差，有待在下一步工作中繼續(xù)開展深入研究。

?×H=0；

(A1)

?·B=0。

(A2)

在高斯單位制(CGSM)中，B(Gs)和H(Gs)的關(guān)系為

B=H+M。

(A3)

式中：M是磁化強度，單位為emu，由感應(yīng)磁化強度Ji和剩余磁化強度Jr組成。

M=Ji+Jr。

(A4)

且Ji與H成正比：

Ji=κH。

(A5)

式中，κ是磁化率。將式(A4)、式(A5)代入式(A3)可得

B=(1+κ)H+Jr=μH+Jr。

(A6)

式中，μ=1+κ是磁導(dǎo)率。將式(A6)代入式(A2)可得

?·[(1+κ)H+Jr]=0。

(A7)

由于H由正常磁場強度Ho和異常磁場強度Ha組成，即可表示為H=Ho+Ha，令

Ha=-?U。

(A8)

式中，U為磁位。將式(A8)代入式(A7)：

?·[(1+κ)?U]=?·Jr+?·[(1+κ)Ho]。

(A9)

又由于?·Ho=0，且在只考慮弱磁性體的條件下，將式(A9)移項化簡，有

?2U=?·(κHo+Jr)。

(A10)

式(A10)即為同時考慮感應(yīng)磁化強度和剩余磁化強度的磁位基本方程。對于弱磁性體，式(A10)可以寫為

?2U=?·M。

(A11)

在實際磁法勘探中，一般采用CGSM單位制,而CGSM單位制中磁場強度以奧斯特(Oe)為單位，其分?jǐn)?shù)單位是伽馬(γ)，1 γ=10-5Oe；在真空中μ=1，故B和H的量綱相同，它們的單位Gs與Oe大小相等，1 γ=10-5Oe=10-5Gs=10-9T=1 nT。

附錄B

針對邊值問題，構(gòu)造一個泛函：

(B1)

其中：

(B2)

則其變分為

(B3)

將式(2)代入式(B3)可得

(B4)

將式(4)代入式(B4)得

(B5)

所以

(B6)

即空間波數(shù)混合域磁位的邊值問題與下列變分問題等價：

(B7)

附錄C

將整個區(qū)域的單元積分分解為各單元的積分之和，則

(C1)

令

(C2)

(C3)

式(C1)中第一項單元積分為

(C4)

其中，

(C5)

式中，l為單元長度。式(C1)中第二項單元積分為

(C6)

其中，

(C7)

式(C1)中第三項單元積分為

(C8)

利用形函數(shù)將jax、jaz表示為：

(C9)

其中：

(C10)

則式(C8)中第一個單元積分為

(C11)

其中，

(C12)

式(C8)中第2個單元積分為

(C13)

其中，

(C14)

式(C1)中，z=zmin、z=zmax分別為第一個和最后一個節(jié)點，可將其分別擴展成：

(C15)

其中，

(C16)

綜上，可將式(C1)表示為

(C17)