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

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

0 引言

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

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

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

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

1 基本理論

1.1 控制方程

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

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

(1)

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

(2)

1.2 邊界條件

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

(3)

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

(4)

1.3 磁場(chǎng)和磁梯度張量計(jì)算

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

圖1 邊界條件示意圖

(5)

(6)

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

(7)

(8)

1.4 零波數(shù)處理

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

(9)

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

(10)

(11)

(12)

2 常微分方程求解

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

聯(lián)立式(2)和式(4)，即是空間波數(shù)混合域二度體磁位滿足的一維邊值問(wèn)題?；谧兎衷順?gòu)造泛函[20]，可得到與邊值問(wèn)題等價(jià)的變分問(wèn)題(詳見附錄B)：

(13)

(14)

其中：

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

Ku=P。

(15)

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

3 數(shù)值試驗(yàn)

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

3.1 突變介質(zhì)模型

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

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

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

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

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

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

(16)

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

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

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

3.3 計(jì)算效率

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

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

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

3.4 起伏地形模型

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

地面以上100 m高度的磁場(chǎng)及磁梯度張量模擬結(jié)果見圖7。可以看出：在山脊處，磁場(chǎng)Bz及磁梯度分量Tzz和Txz相對(duì)水平地形模型變化較快，正異常幅值增加，磁場(chǎng)Bx相對(duì)水平地形模型負(fù)異常幅值有所增加；在山谷處，磁場(chǎng)Bz和磁梯度分量Tzz和Txz相對(duì)水平地形模型變化減緩，但負(fù)異常幅值略有增強(qiáng)，磁場(chǎng)Bx相對(duì)水平地形模型變化明顯減緩，進(jìn)一步驗(yàn)證了本文算法在起伏地形條件下的適用性和可靠性。

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

a. 水平地形磁場(chǎng)；b. 水平地形磁梯度；c. 起伏地形磁場(chǎng)；d. 起伏地形磁梯度。

4 結(jié)論

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

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

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

?×H=0；

(A1)

?·B=0。

(A2)

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

B=H+M。

(A3)

式中：M是磁化強(qiáng)度，單位為emu，由感應(yīng)磁化強(qiáng)度Ji和剩余磁化強(qiáng)度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由正常磁場(chǎng)強(qiáng)度Ho和異常磁場(chǎng)強(qiáng)度Ha組成，即可表示為H=Ho+Ha，令

Ha=-?U。

(A8)

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

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

(A9)

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

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

(A10)

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

?2U=?·M。

(A11)

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

附錄B

針對(duì)邊值問(wèn)題，構(gòu)造一個(gè)泛函：

(B1)

其中：

(B2)

則其變分為

(B3)

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

(B4)

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

(B5)

所以

(B6)

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

(B7)

附錄C

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

(C1)

令

(C2)

(C3)

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

(C4)

其中，

(C5)

式中，l為單元長(zhǎng)度。式(C1)中第二項(xiàng)單元積分為

(C6)

其中，

(C7)

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

(C8)

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

(C9)

其中：

(C10)

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

(C11)

其中，

(C12)

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

(C13)

其中，

(C14)

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

(C15)

其中，

(C16)

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

(C17)