戴世坤，冉應(yīng)強(qiáng)*，張瑩，陳輕蕊，凌嘉宣，賈金榮

(1.中南大學(xué)有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，湖南長沙 410083； 2.中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南長沙 410083)

0 引言

近年來，強(qiáng)磁性體的退磁效應(yīng)逐漸成為磁場(chǎng)正、反演領(lǐng)域的研究熱點(diǎn)。當(dāng)磁性體磁化率低于0.1 SI時(shí)，可忽略退磁效應(yīng)[1]；當(dāng)磁性體磁化率高于0.1 SI時(shí)，退磁效應(yīng)會(huì)顯著影響磁異常的幅值、形態(tài)等，進(jìn)而影響磁測(cè)數(shù)據(jù)的處理和解釋結(jié)果，所以在強(qiáng)磁性體磁場(chǎng)正、反演過程中必須考慮退磁效應(yīng)的影響。考慮到三維模型的計(jì)算量太大，且實(shí)際勘探過程中，當(dāng)?shù)刭|(zhì)體沿走向的尺度遠(yuǎn)大于垂直走向的尺度時(shí)(如斷層、接觸帶等)，這類地質(zhì)體可用走向方向無限延伸的二度體近似，既能滿足勘探要求，也能大大提高計(jì)算效率[2]。

目前，對(duì)于二度體磁異常的研究多是基于弱磁情況，可忽略退磁效應(yīng)。Bhattacharyya等[3-4]、Peder-sen[5]、吳宣志[6]給出了頻率域不同二度體模型的弱磁異常場(chǎng)表達(dá)式；Talwani等[7]、Shuey等[8]、Plouff[9]、Singh等[10]、Won等[11]、賈真[12]提出了空間域二度體數(shù)值模擬算法；對(duì)于頻率域弱磁正演計(jì)算方法，柴玉璞[13]基于偏移抽樣提高了反傅里葉變換數(shù)值精度；Tontinif等[14]研究了快速傅里葉變換擴(kuò)邊法與誤差的關(guān)系；Wu等[15]在偏移采樣的基礎(chǔ)上提出了高斯傅里葉變換方法，有效提高了弱磁場(chǎng)數(shù)值模擬的精度和效率。

考慮退磁效應(yīng)的情況下，強(qiáng)磁場(chǎng)數(shù)值模擬方法主要有以下三類：積分方程法、有限體積法和有限單元法。Eskola等[16]基于積分方程提出了考慮退磁效應(yīng)的靜磁方程，其解的準(zhǔn)確性與代數(shù)方程組的規(guī)模有關(guān)；Ouyang等[17]基于強(qiáng)磁場(chǎng)積分方程，采用棱柱體剖分方法實(shí)現(xiàn)了強(qiáng)磁性體的正演計(jì)算，采用Gauss-FFT方法減小了截?cái)嘈?yīng)；Kostrov[18]提出了一種采用三角單元的體積分方法求解考慮退磁效應(yīng)的二維磁場(chǎng)正演問題；付文祥[19]基于體積分方法提出一種采用三角棱柱單元擬合二度體的方法，并給出了考慮退磁影響的無限長水平圓柱體的總磁場(chǎng)強(qiáng)度異常的理論解析表達(dá)式；歐洋等[20]利用有限體積法實(shí)現(xiàn)了同時(shí)考慮退磁和剩磁的正演模擬；王書惠[21]在不引入均勻磁化設(shè)定、考慮退磁效應(yīng)影響的情況下，利用有限元計(jì)算了強(qiáng)磁性復(fù)雜形體的磁異常；劉雙等[22]利用FlexPDE軟件實(shí)現(xiàn)了二度圓柱體、板狀體的有限元強(qiáng)磁場(chǎng)正演，總結(jié)了退磁作用對(duì)磁異常幅值影響的規(guī)律；另外，劉鑫磊[23]計(jì)算了具有各向異性、非均勻磁化率分布的任意形態(tài)地質(zhì)體的磁場(chǎng)，計(jì)算過程考慮了退磁效應(yīng)，但計(jì)算結(jié)果可能會(huì)不連續(xù)；Purssm等[24]采用迭代方法推導(dǎo)了復(fù)雜地形下高磁化率地質(zhì)體的磁場(chǎng)表達(dá)式，但計(jì)算結(jié)果誤差較大。

現(xiàn)有的強(qiáng)磁場(chǎng)數(shù)值模擬方法大多基于有限單元法、有限體積法等，這些方法最終都?xì)w結(jié)為大型線性方程組的求解，對(duì)于復(fù)雜條件下的大規(guī)模二維強(qiáng)磁場(chǎng)數(shù)值的模擬存在計(jì)算效率低的問題。此外，由重磁位場(chǎng)的研究結(jié)論可知，與總磁場(chǎng)相比，梯度張量對(duì)異常體邊界的識(shí)別能力更強(qiáng)、分辨率更高，而現(xiàn)有關(guān)于強(qiáng)磁場(chǎng)梯度張量的數(shù)值模擬研究較少。因此，本文提出一種適用于大規(guī)模復(fù)雜條件(復(fù)雜磁化率分布和起伏地形)的二度體強(qiáng)磁場(chǎng)及其張量梯度的正演方法——空間—波數(shù)混合域二維強(qiáng)磁場(chǎng)及梯度張量數(shù)值模擬方法。該方法利用一維傅里葉變換將空間域二維強(qiáng)磁性體磁位滿足的二維偏微分方程轉(zhuǎn)換為空間—波數(shù)混合域磁位滿足的一維常微分方程，利用有限單元法對(duì)每個(gè)波數(shù)下的空間—波數(shù)混合域磁位進(jìn)行計(jì)算，并利用追趕法求解定帶寬線性方程組，得到空間—波數(shù)混合域磁位，利用位場(chǎng)與張量梯度之間的關(guān)系得到相應(yīng)的空間—波數(shù)域磁場(chǎng)分量，再經(jīng)過反傅里葉變換求得空間域場(chǎng)。該算法通過迭代對(duì)強(qiáng)磁場(chǎng)進(jìn)行數(shù)值逼近，確保了算法穩(wěn)定和收斂。

1 方法原理

1.1 控制方程

強(qiáng)磁性體二維磁位方程[25]為

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

(1)

關(guān)于磁介質(zhì)受到外部磁場(chǎng)的作用，可用磁化強(qiáng)度M衡量磁化程度，它與磁場(chǎng)強(qiáng)度H之間的關(guān)系為

(2)

對(duì)式(1)沿x方向進(jìn)行一維傅里葉變換，得到空間—波數(shù)混合域磁位方程

(3)

1.2 邊界條件

式(3)的通解為

(4)

式中a和b為常數(shù)。在笛卡爾坐標(biāo)系中，令z軸向下為正，模型的上邊界z=zmin，下邊界z=zmax，根據(jù)上行波和下行波相關(guān)理論，上、下邊界條件分別為

(5)

(6)

1.3 磁場(chǎng)和磁場(chǎng)梯度張量

對(duì)磁位求導(dǎo)可得到相應(yīng)的磁異常場(chǎng)及其梯度張量，空間域磁異常場(chǎng)Ba、磁異常場(chǎng)梯度張量T和總強(qiáng)度磁異常ΔT[26]的計(jì)算公式分別為

(7)

(8)

ΔT=‖B0+Ba‖-‖B0‖

(9)

(10)

(11)

1.4 空間—波數(shù)域常微分方程的求解

聯(lián)立式(3)、式(5)、式(6)，可求解空間—波數(shù)混合域磁位的邊值問題。

對(duì)于二維強(qiáng)磁性體，磁位在空間—波數(shù)混合域滿足邊值問題，可采用基于二次插值的一維有限單元法對(duì)其進(jìn)行求解。該方法保留垂向方向?yàn)榭臻g域，可根據(jù)實(shí)際需求靈活改變垂向網(wǎng)格密度，兼顧了計(jì)算精度和計(jì)算效率；利用有限元法得到的五對(duì)角線性方程組，可利用追趕法[27]實(shí)現(xiàn)方程組的高效、高精度求解。

基于變分原理[25]，得到與邊值問題等價(jià)的變分問題

(12)

(13)

(14)

由于δu≠0，故

Ku=P

(15)

1.5 緊算子

求解空間—波數(shù)混合域磁位的邊值問題時(shí)，磁化強(qiáng)度M的表達(dá)式(式2)中存在未知項(xiàng)Ha，若直接求解，不能得到五對(duì)角線性方程組。本文采用迭代法對(duì)真解進(jìn)行逐次逼近。在電磁法數(shù)值模擬中，Torres-Verdín等[28]提出了一種電磁場(chǎng)擴(kuò)展的Bron近似法，這是一種基于內(nèi)部場(chǎng)的新型近似法，通過將背景電場(chǎng)投影到散射張量上對(duì)電磁場(chǎng)進(jìn)行近似表示，該方法適用于異常體體積較大、電阻率差異明顯及寬頻的情形。Zhdanov等[29]構(gòu)建了一個(gè)收斂的Bron級(jí)數(shù)，Gao[30]進(jìn)一步修改得到新的收斂Bron近似。Ouyang等[17]基于前人研究，根據(jù)電磁場(chǎng)緊算子推導(dǎo)過程構(gòu)造了適用于強(qiáng)磁場(chǎng)穩(wěn)定收斂的迭代格式，其迭代公式為

(16)

式中：i表示迭代次數(shù);H(i)和H(i+1)分別表示第i次和第i+1次迭代得到的總場(chǎng)。

1.6 算法流程

本文提出的適用于大規(guī)模復(fù)雜條件(復(fù)雜磁化率分布和起伏地形)的二度體強(qiáng)磁場(chǎng)及其張量梯度的正演方法，即空間—波數(shù)混合域二維強(qiáng)磁場(chǎng)及梯度張量數(shù)值模擬方法流程見圖1，具體步驟如下：

圖1 空間—波數(shù)混合域二維強(qiáng)磁場(chǎng)及梯度張量數(shù)值模擬迭代算法流程

(1)輸入背景場(chǎng)磁場(chǎng)強(qiáng)度H0，即第一次迭代時(shí)總磁場(chǎng)強(qiáng)度的初始值H(0)；

2 算例分析

首先設(shè)計(jì)一個(gè)二維圓柱體模型和一個(gè)棱柱體模型分別驗(yàn)證本文算法的正確性和高效性，然后設(shè)計(jì)一個(gè)帶起伏地形的、包含一個(gè)順磁性異常體和一個(gè)強(qiáng)磁性異常體的模型，驗(yàn)證正演算法對(duì)復(fù)雜條件的適應(yīng)性。算例采用串行計(jì)算，算法中的傅里葉變換通過Gauss-FFT實(shí)現(xiàn)。算法基于Fortran95語言編程，電腦配置為：Intel Core i3-4150 CPU，主頻為3.50 GHz，內(nèi)存為12 GB。

2.1 正確性驗(yàn)證

設(shè)計(jì)圖2所示二維圓柱體模型，圓柱體沿y方向無限延伸，頂部埋深為300 m。模型的研究區(qū)域?yàn)椋簒=[-2000 m，2000 m]，z=[0，1000 m]；網(wǎng)格數(shù)為800(x)×400(y)，水平采樣間隔為5 m，垂向采樣間隔為2.5 m；觀測(cè)面為z=0。圓柱體磁異常場(chǎng)的x分量Bax、z分量Baz、梯度張量水平分量?Bax/?x和垂向分量?Bax/?z的解析表達(dá)式詳見附錄B。該模型的總強(qiáng)度磁異常ΔT、磁異常場(chǎng)Ba及其梯度張量T的解析解和數(shù)值解及觀測(cè)面上各點(diǎn)的相對(duì)誤差分別見圖3、圖4和圖5。

圖2 圓柱體模型示意圖

圖3 圓柱體模型總強(qiáng)度磁異常ΔT解析解和數(shù)值解(a)及其相對(duì)誤差(b)

圖4 圓柱體模型磁異常場(chǎng)Bax(上)和Baz(下)的解析解和數(shù)值解(a)及其相對(duì)誤差曲線(b)

圖5 圓柱體模型磁異常場(chǎng)梯度張量?Bax/?x(上)和?Bax/?z(下)的解析解和數(shù)值解(a)及其相對(duì)誤差(b)

由圖3～圖5可以看出，ΔT的相對(duì)誤差小于0.5%，Bax和Baz的相對(duì)誤差小于0.5%，?Bax/?x和?Bax/?z的相對(duì)誤差小于0.1%，即ΔT、Bax、Baz、?Bax/?x和?Bax/?z的解析解和數(shù)值解基本吻合，僅在極少數(shù)零值點(diǎn)附近相對(duì)誤差較大。該圓柱體模型的正演結(jié)果驗(yàn)證了本算法的正確性。

圖6為該模型的磁異常水平分量Bax和垂直分量Baz的相對(duì)誤差隨迭代次數(shù)的變化曲線?？梢钥闯觯?jīng)9次迭代后，計(jì)算結(jié)果即可滿足迭代終止條件。

圖6和圖7表明本文迭代算法適應(yīng)于不同磁化率的模型，收斂穩(wěn)定。

圖6 磁異常分量相對(duì)誤差隨迭代次數(shù)變化曲線

圖7 迭代次數(shù)隨圓柱體磁化率的變化曲線

2.2 算法效率分析

圖8 棱柱體模型示意圖

參照文獻(xiàn)[31]中利用有限體積法計(jì)算得到的數(shù)值解，衡量COMSOL軟件和本文算法的精度。引入相對(duì)均方根誤差rrms[33]統(tǒng)計(jì)整條觀測(cè)線上總強(qiáng)度磁異常ΔT的相對(duì)誤差

(17)

式中：N表示觀測(cè)總點(diǎn)數(shù)；Xi表示COMSOL或者本文算法計(jì)算得到的第i個(gè)觀測(cè)點(diǎn)的總強(qiáng)度磁異常；Yi表示有限體積法計(jì)算得到的第i個(gè)觀測(cè)點(diǎn)的總強(qiáng)度磁異常。

當(dāng)網(wǎng)格數(shù)為200×100時(shí)，采用本文算法的rrms為2.84%，計(jì)算時(shí)間0.75 s，占用內(nèi)存12.0 MB。表1為不同網(wǎng)格剖分方案下基于COMSOL軟件的計(jì)算精度、耗時(shí)和內(nèi)存占用?；贑OMSOL計(jì)算時(shí)需對(duì)研究范圍進(jìn)行擴(kuò)邊處理，當(dāng)網(wǎng)格擴(kuò)邊4倍時(shí)，rrms為2.92%，計(jì)算時(shí)間618 s，占用內(nèi)存2387.1 MB。因此，與COMSOL軟件相比，在獲得相似精度的情況下，本文算法耗時(shí)更短，占用內(nèi)存更少，計(jì)算效率更高。

表1 COMSOL計(jì)算時(shí)間與占用內(nèi)存統(tǒng)計(jì)

網(wǎng)格數(shù)rrms%計(jì)算時(shí)間s占用內(nèi)存MB200×100(不擴(kuò)邊)110.901480.5600×300(擴(kuò)邊2倍)9.57130887.41000×500(擴(kuò)邊4倍)2.926182387.1

有限體積法、COMSOL軟件(擴(kuò)邊4倍)和本文算法的總強(qiáng)度磁異常見圖9a，以有限體積法的數(shù)值解為參照，采用COMSOL軟件和本文方法計(jì)算結(jié)果的相對(duì)誤差見圖9b?？梢钥闯觯捎肅OMSOL軟件和空間—波數(shù)域方法計(jì)算得到的總強(qiáng)度磁異常與有限體積法的計(jì)算結(jié)果圖形均大致吻合。

圖9 棱柱體模型不同方法得到的總強(qiáng)度磁異常ΔT(a)及相對(duì)誤差(b)

對(duì)于圖8模型，改變模型網(wǎng)格剖分個(gè)數(shù)，分析本文算法計(jì)算時(shí)間和占用內(nèi)存隨網(wǎng)格節(jié)點(diǎn)數(shù)的變化，結(jié)果見圖10和表2。對(duì)比表1與表2，可見本文算法計(jì)算時(shí)間更短、占用內(nèi)存更少，占用內(nèi)存和計(jì)算時(shí)間與網(wǎng)格節(jié)點(diǎn)數(shù)均呈現(xiàn)近似線性正相關(guān)變化趨勢(shì)。

圖10 棱柱體模型本文算法計(jì)算時(shí)間和占用內(nèi)存隨網(wǎng)格剖分節(jié)點(diǎn)數(shù)的變化

表2 棱柱體模型本文算法計(jì)算時(shí)間與占用內(nèi)存統(tǒng)計(jì)

網(wǎng)格數(shù)計(jì)算時(shí)間/s占用內(nèi)存/MB200×1000.7512.0400×2002.9443.1600×3005.96108.8800×4009.50183.81000×50014.39280.2

2.3 適應(yīng)性驗(yàn)證

為了驗(yàn)證本文正演算法對(duì)復(fù)雜模型的適應(yīng)性，設(shè)計(jì)一個(gè)起伏地形模型(圖11)，模型同時(shí)包含順磁性異常體和強(qiáng)磁性異常體。模型的背景磁場(chǎng)T0=50000 nT，磁傾角α=45°，磁偏角β=5°。研究區(qū)域：x=[-2000 m，2000 m]，網(wǎng)格數(shù)為400；z=[-320 m，820 m]，其中z方向區(qū)間[320 m，820 m]為起伏地形最低點(diǎn)以下的計(jì)算區(qū)域，網(wǎng)格數(shù)為100。起伏地形高程滿足函數(shù)

圖11 起伏地形模型示意圖

x∈[-2000,2000]

(18)

基于該模型研究起伏地形條件下插值平面?zhèn)€數(shù)對(duì)數(shù)值解計(jì)算精度的影響。起伏地形在z方向的區(qū)域?yàn)閇-320 m，320 m]，該方向采用均勻網(wǎng)格剖分，將插值平面?zhèn)€數(shù)分別設(shè)為5、9、13、17、21、25、29和33。采用式(17)計(jì)算整條起伏觀測(cè)線上的總強(qiáng)度磁異常ΔT、磁異常場(chǎng)Ba及梯度張量T的解析解。

圖12為該起伏地形模型的總強(qiáng)度磁異常、磁異常場(chǎng)及梯度張量的解析解和數(shù)值解的相對(duì)均方根誤差隨插值平面?zhèn)€數(shù)的變化曲線。可以看出，總強(qiáng)度磁異常、磁異常場(chǎng)與其梯度張量的rrms均隨插值平面?zhèn)€數(shù)的增加而降低。采用5個(gè)平面進(jìn)行插值時(shí)，起伏地形上總強(qiáng)度磁異常、磁異常場(chǎng)與其梯度張量的rrms均大于1%；當(dāng)采用25個(gè)及以上不同高度平面進(jìn)行插值時(shí)，起伏地形上總強(qiáng)度磁異常、磁異常場(chǎng)與其梯度張量的rrms均低于1%，約為0.9%，說明采用25個(gè)插值平面即可達(dá)到滿意的計(jì)算精度。

圖12 rrms隨插值平面數(shù)的變化曲線

另外，采用25個(gè)插值平面時(shí)，耗時(shí)約2.330 s，占用內(nèi)存52.8 MB，表明本文算法計(jì)算速度快、占用內(nèi)存小，適用于復(fù)雜地形下二度體磁異常場(chǎng)及其梯度張量的計(jì)算，計(jì)算效率較高。

采用25個(gè)高度平面進(jìn)行插值，基于本文方法得到的總強(qiáng)度磁異常、磁異常場(chǎng)及其梯度張量分別見圖13、圖14和圖15。可以看出，解析解和數(shù)值解的吻合度很高，表明了本文算法對(duì)復(fù)雜地形模型的計(jì)算精度較高。

圖13 起伏地形模型總強(qiáng)度磁異常解析解與數(shù)值解對(duì)比

圖14 起伏地形模型磁異常場(chǎng)解析解和數(shù)值解對(duì)比

(a)Bax； (b)Baz

圖15 起伏地形模型磁異常場(chǎng)梯度張量解析解和數(shù)值解對(duì)比

(a)?Bax/?x； (b)?Bax/?z

3 結(jié)論

針對(duì)強(qiáng)磁性體的退磁效應(yīng)不能忽略的問題，本文提出了一種空間—波數(shù)混合域二度體磁異常場(chǎng)及其梯度張量正演方法。該算法充分利用了傅里葉變換的快速性、追趕法求解的高效性和迭代算法的穩(wěn)定性的優(yōu)勢(shì)，實(shí)現(xiàn)了二維強(qiáng)磁性體磁異常場(chǎng)及其梯度張量的高效、高精度數(shù)值模擬。設(shè)計(jì)了水平地形圓柱體模型進(jìn)行數(shù)值試驗(yàn)，通過對(duì)比解析解與數(shù)值解驗(yàn)證了本算法的正確性和較高的精度；以有限體積法數(shù)值解為參照，對(duì)比了本文算法和COMSOL軟件的計(jì)算效率，證明了本算法的高效性；利用起伏地形下同時(shí)包含順磁性異常體和強(qiáng)磁性異常體的復(fù)雜模型驗(yàn)證了本算法對(duì)復(fù)雜地質(zhì)條件的適應(yīng)性。另外，在本文算法的計(jì)算環(huán)節(jié)中，正、反傅里葉變換和求解常微分方程這兩個(gè)環(huán)節(jié)均有很好的并行性，采用并行方式可進(jìn)一步提高計(jì)算效率。本文為地下二度體磁異常場(chǎng)及其張量梯度的正演提供了一種高效、高精度算法，對(duì)于二度體磁異常的反演和人機(jī)交互正、反演解釋具有重要意義。

附錄A 三維強(qiáng)磁場(chǎng)變分問題求解

求解變分問題(式(12))時(shí)，可將其寫為

(A-1)

式中

(A-2)

令

(A-3)

式中：j、m分別為單元兩端節(jié)點(diǎn)坐標(biāo)；p為單元中點(diǎn)坐標(biāo)；積分區(qū)域q為單元范圍?？傻?/p>

(A-4)

式(A-1)中第一項(xiàng)單元積分為

(A-5)

其中

(A-6)

式中l(wèi)為單元長度。

式(A-1)中第二項(xiàng)單元積分為

(A-7)

其中

(A-8)

式(A-1)中第三項(xiàng)單元積分為

(A-9)

利用形函數(shù)可將上式中的Jaz表示為

Jaz=JazjNj+JazpNp+JazmNm

(A-10)

其中

szq=(Jazj,Jazp,Jazm)T

(A-11)

則式(A-9)中右端項(xiàng)可表示為

(A-12)

其中

(A-13)

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

(A-14)

其中

(A-15)

(A-16)

綜上所述，式(A-1)可表示為

F(u)=uTKu-2uTP

(A-17)

附錄B二維圓柱體強(qiáng)磁場(chǎng)及其梯度張量解析解

Ha=

(B-1)

(B-2)

沿走向無限長圓柱體的磁場(chǎng)梯度張量

(B-3)

其中

H0z(z-z0)[(z-z0)2-3(x-x0)2]]

(B-4)

H0z(x-x0)[(x-x0)2-3(z-z0)2]]

(B-5)

式中H0x、H0z分別為背景磁場(chǎng)H0的x分量和z分量。