空間—波數域二維強磁場及其梯度張量數值模擬研究

戴世坤，冉應強*，張瑩，陳輕蕊，凌嘉宣，賈金榮

(1.中南大學有色金屬成礦預測與地質環境監測教育部重點實驗室，湖南長沙 410083； 2.中南大學地球科學與信息物理學院，湖南長沙 410083)

0 引言

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

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

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

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

1 方法原理

1.1 控制方程

強磁性體二維磁位方程[25]為

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

(1)

關于磁介質受到外部磁場的作用，可用磁化強度M衡量磁化程度，它與磁場強度H之間的關系為

(2)

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

(3)

1.2 邊界條件

式(3)的通解為

(4)

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

(5)

(6)

1.3 磁場和磁場梯度張量

對磁位求導可得到相應的磁異常場及其梯度張量，空間域磁異常場Ba、磁異常場梯度張量T和總強度磁異常ΔT[26]的計算公式分別為

(7)

(8)

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

(9)

(10)

(11)

1.4 空間—波數域常微分方程的求解

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

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

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

(12)

(13)

(14)

由于δu≠0，故

Ku=P

(15)

1.5 緊算子

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

(16)

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

1.6 算法流程

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

圖1 空間—波數混合域二維強磁場及梯度張量數值模擬迭代算法流程

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

2 算例分析

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

2.1 正確性驗證

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

圖2 圓柱體模型示意圖

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

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

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

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

圖6為該模型的磁異常水平分量Bax和垂直分量Baz的相對誤差隨迭代次數的變化曲線?？梢钥闯?，經9次迭代后，計算結果即可滿足迭代終止條件。

圖6和圖7表明本文迭代算法適應于不同磁化率的模型，收斂穩定。

圖6 磁異常分量相對誤差隨迭代次數變化曲線

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

2.2 算法效率分析

圖8 棱柱體模型示意圖

參照文獻[31]中利用有限體積法計算得到的數值解，衡量COMSOL軟件和本文算法的精度。引入相對均方根誤差rrms[33]統計整條觀測線上總強度磁異常ΔT的相對誤差

(17)

式中：N表示觀測總點數；Xi表示COMSOL或者本文算法計算得到的第i個觀測點的總強度磁異常；Yi表示有限體積法計算得到的第i個觀測點的總強度磁異常。

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

表1 COMSOL計算時間與占用內存統計

網格數rrms%計算時間s占用內存MB200×100(不擴邊)110.901480.5600×300(擴邊2倍)9.57130887.41000×500(擴邊4倍)2.926182387.1

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

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

對于圖8模型，改變模型網格剖分個數，分析本文算法計算時間和占用內存隨網格節點數的變化，結果見圖10和表2。對比表1與表2，可見本文算法計算時間更短、占用內存更少，占用內存和計算時間與網格節點數均呈現近似線性正相關變化趨勢。

圖10 棱柱體模型本文算法計算時間和占用內存隨網格剖分節點數的變化

表2 棱柱體模型本文算法計算時間與占用內存統計

網格數計算時間/s占用內存/MB200×1000.7512.0400×2002.9443.1600×3005.96108.8800×4009.50183.81000×50014.39280.2

2.3 適應性驗證

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

圖11 起伏地形模型示意圖

x∈[-2000,2000]

(18)

基于該模型研究起伏地形條件下插值平面個數對數值解計算精度的影響。起伏地形在z方向的區域為[-320 m，320 m]，該方向采用均勻網格剖分，將插值平面個數分別設為5、9、13、17、21、25、29和33。采用式(17)計算整條起伏觀測線上的總強度磁異常ΔT、磁異常場Ba及梯度張量T的解析解。

圖12為該起伏地形模型的總強度磁異常、磁異常場及梯度張量的解析解和數值解的相對均方根誤差隨插值平面個數的變化曲線?？梢钥闯?，總強度磁異常、磁異常場與其梯度張量的rrms均隨插值平面個數的增加而降低。采用5個平面進行插值時，起伏地形上總強度磁異常、磁異常場與其梯度張量的rrms均大于1%；當采用25個及以上不同高度平面進行插值時，起伏地形上總強度磁異常、磁異常場與其梯度張量的rrms均低于1%，約為0.9%，說明采用25個插值平面即可達到滿意的計算精度。

圖12 rrms隨插值平面數的變化曲線

另外，采用25個插值平面時，耗時約2.330 s，占用內存52.8 MB，表明本文算法計算速度快、占用內存小，適用于復雜地形下二度體磁異常場及其梯度張量的計算，計算效率較高。

采用25個高度平面進行插值，基于本文方法得到的總強度磁異常、磁異常場及其梯度張量分別見圖13、圖14和圖15?？梢钥闯?，解析解和數值解的吻合度很高，表明了本文算法對復雜地形模型的計算精度較高。

圖13 起伏地形模型總強度磁異常解析解與數值解對比

圖14 起伏地形模型磁異常場解析解和數值解對比

(a)Bax； (b)Baz

圖15 起伏地形模型磁異常場梯度張量解析解和數值解對比

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

3 結論

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

附錄A 三維強磁場變分問題求解

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

(A-1)

式中

(A-2)

令

(A-3)

式中：j、m分別為單元兩端節點坐標；p為單元中點坐標；積分區域q為單元范圍?？傻?/p>

(A-4)

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

(A-5)

其中

(A-6)

式中l為單元長度。

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

(A-7)

其中

(A-8)

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

(A-9)

利用形函數可將上式中的Jaz表示為

Jaz=JazjNj+JazpNp+JazmNm

(A-10)

其中

szq=(Jazj,Jazp,Jazm)T

(A-11)

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

(A-12)

其中

(A-13)

式(A-1)中，z=zmin、z=zmax分別為第一個和最后一個節點，可將其分別擴展為

(A-14)

其中

(A-15)

(A-16)

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

F(u)=uTKu-2uTP

(A-17)

附錄B二維圓柱體強磁場及其梯度張量解析解

Ha=

(B-1)

(B-2)

沿走向無限長圓柱體的磁場梯度張量

(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分別為背景磁場H0的x分量和z分量。