郭楚楓,張世暉,劉天佑

(中國(guó)地質(zhì)大學(xué)(武漢) 地球物理與空間信息學(xué)院，湖北 武漢 430074)

0 引言

有限單元法是地球物理數(shù)值正演模擬方法的一種，以變分問題為基礎(chǔ)，根據(jù)地球物理中的偏微分方程和邊界條件，用數(shù)值方法近似計(jì)算地球物理場(chǎng)值，適用于復(fù)雜物性分布和復(fù)雜邊界形狀的地球物理計(jì)算[1]。自1971年Coggon首次系統(tǒng)地闡述了有限單元法在電、電磁法中的應(yīng)用以來[2]，有限單元法已經(jīng)廣泛應(yīng)用于復(fù)雜形態(tài)、復(fù)雜地質(zhì)條件的地球物理正演中。

在電磁法領(lǐng)域，基于非結(jié)構(gòu)化網(wǎng)格的矢量有限單元法，具有很強(qiáng)的適應(yīng)性，是目前國(guó)內(nèi)外研究的熱點(diǎn)[3-5]。隨著研究的深入，各向同性假設(shè)模型已難以滿足解釋的需要，基于有限元的各向異性研究受到了廣泛關(guān)注[6-8]。在重力勘探中，有限單元法常用來進(jìn)行梯度張量及重力矢量的正演模擬計(jì)算[9-12]，Cai和Wang[11]實(shí)現(xiàn)了利用有限元法計(jì)算非均勻復(fù)雜形體重力異常的快速算法，朱自強(qiáng)等[13]則利用有限單元法進(jìn)行重力異常的地形校正。在磁場(chǎng)正演模擬中，有限單元法的優(yōu)勢(shì)在于不需要引入均勻磁化設(shè)定，能同時(shí)考慮磁各向異性、退磁效應(yīng)、剩磁三種因素的影響[14]，所以在充分考慮退磁作用、進(jìn)行高精度磁測(cè)資料反演解釋時(shí)，常用有限單元法進(jìn)行正演模擬[15-16]。

上述的有限元計(jì)算存在相同的問題，即采用傳統(tǒng)截?cái)噙吔?，在一個(gè)相對(duì)較大的區(qū)域內(nèi)，將無窮的邊界問題近似為有限區(qū)域進(jìn)行[17]。而有限單元法是一種區(qū)域型方法，必須在全區(qū)域進(jìn)行剖分，對(duì)于地球物理正演中的無界區(qū)域問題，在邊界處產(chǎn)生畸變會(huì)影響對(duì)異常的解釋。一般的解決方法是進(jìn)行擴(kuò)邊，即設(shè)定更大的求解區(qū)域以減小傳統(tǒng)截?cái)噙吔绲挠绊憽Ｈ缡Y甫玉等[9]進(jìn)行三度體重力矢量有限元正演計(jì)算時(shí)提出，設(shè)置求解域邊界長(zhǎng)度應(yīng)大于場(chǎng)源體長(zhǎng)度的7.5倍以達(dá)到理想效果；Ren和Kalscheuer等[5]進(jìn)行三維有限元大地電磁正演時(shí)，設(shè)置了邊長(zhǎng)為20 km的求解域，以保證長(zhǎng)度100 m的板狀體模型正演能達(dá)到模擬精度。擴(kuò)邊往往造成有限元計(jì)算區(qū)域大，節(jié)點(diǎn)數(shù)多，進(jìn)一步降低了有限元正演的計(jì)算效率[17]，較合理的策略是將計(jì)算區(qū)域限定在較小的目標(biāo)區(qū)域范圍內(nèi)，減少單元與節(jié)點(diǎn)數(shù)量，以節(jié)約計(jì)算資源，提高正演效率。利用無限單元取代傳統(tǒng)的邊界條件，通過坐標(biāo)映射將“無限單元”沿某一方向無限擴(kuò)展，最后結(jié)合有限單元法對(duì)內(nèi)部區(qū)域進(jìn)行正演模擬，提高正演效率。

無限元的概念最早由Ungless在1973年給出，并實(shí)現(xiàn)了衰減型無限單元[18]。無限元大體可以分為3類：Bettess映射無限元[19-20]、Astley波包映射無限元[21-22](mapped wave envelope elements)和Burnett多級(jí)展開無限元[23-24]，其中Astley波包映射無限元應(yīng)用最廣。無限元的基本思路是通過坐標(biāo)變換，采用和位場(chǎng)不同的插值函數(shù)，通過坐標(biāo)映射實(shí)現(xiàn)無限域積分，并在映射后的局部坐標(biāo)里建立位場(chǎng)插值形函數(shù)，使位場(chǎng)在有限元網(wǎng)格和無限元網(wǎng)格結(jié)合處保持一致性，而在無限遠(yuǎn)處衰減為0[25]。目前，無限元廣泛應(yīng)用于聲學(xué)、巖土力學(xué)等方面[26-27]，在地球物理學(xué)中則主要應(yīng)用于地震、測(cè)井和電法領(lǐng)域：Fu和Wu[28]應(yīng)用無限元處理了彈性波的吸收邊界條件；朱軍等[29]則利用無限元模擬測(cè)井中面臨的無窮邊界問題；湯井田等[30]應(yīng)用無限元進(jìn)行了三維直流電阻率的數(shù)值模擬；張林成等[17]應(yīng)用無限元進(jìn)行了基于二次場(chǎng)的可控源電磁法的三維數(shù)值模擬。目前，混合有限元—無限元在磁場(chǎng)正演模擬中應(yīng)用較少。為提高有限元三維磁場(chǎng)正演模擬的效率，本文基于COMSOL Multiphysics軟件，在求解域外部邊界上設(shè)置無限元，綜合考慮剩磁、退磁及地表起伏，通過一系列的數(shù)值模擬，與解析算法和傳統(tǒng)有限元方法進(jìn)行對(duì)比，一方面驗(yàn)證了方法的有效性，另一方面證明本算法在保證精度的同時(shí)能夠大大提高正演效率。

1 方法理論

1.1 有限單元法原理

對(duì)于地球物理中磁場(chǎng)的正演問題，其邊界條件是依據(jù)麥克斯韋方程組給出的：

(1)

即在無電流源的磁場(chǎng)空間中磁場(chǎng)強(qiáng)度H和磁感應(yīng)強(qiáng)度B，應(yīng)該滿足：

H和B的關(guān)系為：

B=μ0(H+M)，

其中：μ0為真空中的磁導(dǎo)率，M為磁化強(qiáng)度，由感應(yīng)磁化強(qiáng)度Mi和剩余磁化強(qiáng)度Mr組成

M=Mi+Mr，

而感應(yīng)磁化強(qiáng)度Mi與磁場(chǎng)強(qiáng)度H成正比

Mi=χH，

其中：χ為磁化率，由以上關(guān)系式可以推出

(2)

圖1 非均勻介質(zhì)分布(改自徐世浙，1994)Fig.1 Distribution of inhomogeneous medium(modified from Xu Shizhe，1994)

帶入式(2)中，得

(3)

這是包括感磁和剩磁的磁標(biāo)位基本方程。

圖1為磁場(chǎng)正演中非均勻介質(zhì)分布的簡(jiǎn)單模型，Γ0為兩介質(zhì)共有邊界，Γ1、Γ2為兩介質(zhì)的外部邊界，?！逓閰^(qū)域外部邊界，Vm1、Vm2、μr1、μr2、Mr1、Mr2分別代表界面兩側(cè)介質(zhì)中的磁位、相對(duì)磁導(dǎo)率和剩余磁化強(qiáng)度。徐世浙[1]已經(jīng)證明，磁位邊值問題中的內(nèi)部邊界條件屬于自然邊界條件，使用有限單元法解偏微分方程(3)時(shí)，無需討論磁位的內(nèi)部邊界條件。因此有限單元法具有解決非均勻介質(zhì)的磁場(chǎng)正演能力。

給定Vm的外邊界條件，將計(jì)算區(qū)域Ω取足夠大，使得區(qū)域的邊界?！夼c磁性體足夠遠(yuǎn)，這時(shí)磁性體引起的磁位在邊界上近似為零。但在磁法勘探中存在著正常地磁場(chǎng)T，?！奚系拇盼皇怯蒚決定的，所以?！奚嫌械谝活愡吔鐥l件：

Vm=-T·r， ∈?！?；

(4)

式中：r是坐標(biāo)原點(diǎn)至邊界點(diǎn)的矢徑。轉(zhuǎn)為第二類邊界條件

(5)

針對(duì)基本方程(3)及邊界條件(5)，可以構(gòu)造泛函

(6)

由于內(nèi)部邊界條件為自然邊界條件，考慮?！薜倪吔鐥l件式(5)，以及自由空間中μr=1，Mr=0，式(6)的變分為

因此，無電流靜磁場(chǎng)的磁位的邊值問題與下列的變分問題相應(yīng)：

δF(Vm)=0。

(7)

1.2 有限元—無限元耦合算法

有限元—無限元耦合算法將整個(gè)求解區(qū)域分為有限元區(qū)域和無限元區(qū)域，用無限元區(qū)域代替?zhèn)鹘y(tǒng)的外邊界，在兩個(gè)區(qū)域分別采取有限元分析和無限元分析，通過總體剛度矩陣組裝將兩者結(jié)合起來從而進(jìn)行數(shù)值求解。

圖2a中劃分了有限元和無限元計(jì)算區(qū)域，圖中有限元區(qū)域?yàn)槟繕?biāo)區(qū)域，包含場(chǎng)源、目標(biāo)體及測(cè)點(diǎn)等；無限元區(qū)域?yàn)橛邢拊獏^(qū)域邊界至無窮遠(yuǎn)，為邊界計(jì)算區(qū)域。無限元分析就是在該區(qū)域的某一方向上，采用無限元映射和形函數(shù)，將整體坐標(biāo)映射到局部坐標(biāo)上來，其原理與有限元分析相同。

使用有限單元法計(jì)算時(shí)，采用四面體單元網(wǎng)格剖分，線性插值單元中的插值函數(shù)為：

(8)

式中：Ni為插值形函數(shù)，與四面體自然坐標(biāo)形同；Vmi為節(jié)點(diǎn)上的磁標(biāo)位值。

圖2 三維無限元映射(改自湯井田等，2010)Fig.2 3-D infinite element mapping(modified from Tang et al.，2010)

無限元單元分析時(shí)，采用六節(jié)點(diǎn)，圖2a給出了三維無限元映射示意，圖中節(jié)點(diǎn)1、2、3、4、5、6為無限元基本要素，P點(diǎn)為坐標(biāo)原點(diǎn)，節(jié)點(diǎn)1、2、3為有限元邊界上某單元的三節(jié)點(diǎn)，P點(diǎn)到節(jié)點(diǎn)4、5、6的距離分別為P點(diǎn)到節(jié)點(diǎn)1、2、3距離的兩倍。無限單元最外圍的3個(gè)節(jié)點(diǎn)為無窮遠(yuǎn)。無限元分析就是通過無限元映射將無窮坐標(biāo)映射到圖中的局部坐標(biāo)上。

圖2b中ζ方向表示無限元映射的方向，在ξ-η平面內(nèi)，無限單元與有限單元采用相同的映射形式和形函數(shù)。無限元的坐標(biāo)映射為

式中：Li為三角形面積坐標(biāo)。結(jié)合圖2有：

(9)

結(jié)合ζ方向的坐標(biāo)映射關(guān)系式，可以得到六節(jié)點(diǎn)無限單元的結(jié)點(diǎn)坐標(biāo)映射函數(shù)：

(10)

其中：

經(jīng)過上述映射之后，磁位借助無限單元延伸至無限遠(yuǎn)處，并且在無限遠(yuǎn)處衰減為零。對(duì)于無限元形函數(shù)，其具體的表達(dá)式如下：

(11)

該形函數(shù)是Astley映射無限元理論[21]中采用的形函數(shù)中的系數(shù)(1-ζ)/2與二階Lagrange插值多項(xiàng)式的乘積。

對(duì)有限元—無限元進(jìn)行單元分析后，對(duì)方程組進(jìn)行總體合成，再解相應(yīng)的大型線性方程組即可得到各個(gè)節(jié)點(diǎn)的磁位。得到各節(jié)點(diǎn)的磁位后，通過磁場(chǎng)強(qiáng)度與磁位的關(guān)系，求得磁場(chǎng)的各個(gè)分量：

(12)

2 模型試算

為說明有限元—無限元算法的優(yōu)勢(shì)，綜合考慮剩磁和退磁的影響，設(shè)計(jì)了兩種模型：球體模型和組合體模型。使用COMSOL Multiphysics軟件進(jìn)行建模，以下數(shù)值模擬平臺(tái)均為Intel(R) Xeon(R) CPU3.91G，64GB RAM，16CPUs。

2.1 球體模型

在磁場(chǎng)正演模擬中，忽略退磁和剩磁會(huì)對(duì)正演結(jié)果產(chǎn)生較大的影響，忽略剩磁的相對(duì)誤差與Q(科尼斯布格比)相關(guān)，而忽略退磁的相對(duì)誤差與磁化率相關(guān)，若Q=0.5，忽略剩磁的相對(duì)誤差可達(dá)到33%以上[31]。相比于一般的正演方法，有限單元法的優(yōu)勢(shì)在于適用于復(fù)雜條件、復(fù)雜形狀的地球物理計(jì)算。為突出有限單元法的優(yōu)勢(shì)，驗(yàn)證基于有限元—無限元耦合的磁場(chǎng)正演算法的準(zhǔn)確性，設(shè)計(jì)了如圖3所示的模型。由于球體的退磁系數(shù)為常數(shù)，可以利用退磁改正公式計(jì)算磁異常，將傳統(tǒng)有限元方法，有限元—無限元耦合方法和解析解進(jìn)行對(duì)比。

圖3 球體模型及網(wǎng)格剖分示意(藍(lán)色為球體模型，紅色直線代表測(cè)線，黃色為觀測(cè)平面范圍)Fig.3 Diagram of the Sphere model and mesh generation(the blue region is the sphere, the red line is observation line, the yellow plane is the range of observation plane)

在相同的求解域設(shè)置和相同網(wǎng)格剖分條件下，分別使用有限元—無限元耦合算法和傳統(tǒng)有限元方法進(jìn)行數(shù)值模擬。圖4a為經(jīng)過退磁改正后的球體ΔT磁異常特征圖，圖4b為有限元—無限元耦合算法的數(shù)值模擬結(jié)果，圖4c為傳統(tǒng)有限單元法數(shù)值模擬結(jié)果。圖4d、圖4e為平面誤差分布。圖4b中基于有限元—無限元耦合算法的結(jié)果與圖4a中解析解的異常形態(tài)特征及幅值基本一致，圖4d中，有限元—無限元耦合算法在邊界處幾乎沒有產(chǎn)生誤差。傳統(tǒng)有限元數(shù)值模擬結(jié)果與解析解差別較大，圖4c和圖4e中，求解域邊長(zhǎng)為200 m的傳統(tǒng)有限元數(shù)值模擬得到的ΔT磁異常幅值與解析解相比，最大偏差達(dá)到837 nT，而在邊界處異常形態(tài)差別較大。此例說明，傳統(tǒng)有限單元法進(jìn)行數(shù)值模擬時(shí)在邊界處會(huì)產(chǎn)生較大的誤差，嚴(yán)重影響了正演模擬的精度。而在相同的計(jì)算范圍內(nèi)，基于有限元—無限元耦合的算法的數(shù)值模擬結(jié)果非常優(yōu)秀，特別在邊界附近達(dá)到了很高的精度，明顯優(yōu)于傳統(tǒng)有限元算法的結(jié)果。

擴(kuò)大求解域的范圍，利用傳統(tǒng)有限單元法進(jìn)行數(shù)值模擬，再與有限元—無限元耦合算法進(jìn)行對(duì)比。求解域一的長(zhǎng)寬高分別為300、300、300 m，求解域二為400、400、400 m。圖5為相同磁場(chǎng)參數(shù)的球體模型，在相同網(wǎng)格剖分條件下，使用不同方法計(jì)算得到的ΔT磁異常曲線。曲線1為解析解，曲線2為基于有限元—無限元的計(jì)算結(jié)果，曲線3為求解域邊界長(zhǎng)度為200 m的傳統(tǒng)有限元計(jì)算結(jié)果，曲線4為求解域邊界長(zhǎng)度為300 m的傳統(tǒng)有限元計(jì)算結(jié)果，曲線5為求解域邊界長(zhǎng)度為400 m的傳統(tǒng)有限元計(jì)算結(jié)果。

圖4 不同方法正演的球體平面ΔT磁異常及其平面誤差分布Fig.4 The plane total-field anomaly of the sphere using different methods

圖5 不同方法正演的球體磁異常曲線Fig.5 The total-field anomaly curve of the sphere using different methods

在進(jìn)行三維有限元數(shù)值模擬時(shí)，由于受到6個(gè)方向的截?cái)噙吔缬绊懀褂脗鹘y(tǒng)有限單元法進(jìn)行三維正演時(shí)，若要提高精度，需要在空間的6個(gè)方向?qū)δＰ瓦M(jìn)行延展。圖5和圖6中，對(duì)比解析解(曲線1)和傳統(tǒng)有限元(曲線3、曲線4和曲線5)的結(jié)果，隨著求解區(qū)域的逐漸增大，傳統(tǒng)有限單元法在觀測(cè)平面邊界處產(chǎn)生的誤差逐漸減小，正演精度也在逐漸提高。當(dāng)求解域邊界長(zhǎng)度達(dá)到400 m時(shí)，即對(duì)應(yīng)曲線5的結(jié)果，傳統(tǒng)有限元的數(shù)值模擬結(jié)果和解析計(jì)算結(jié)果已經(jīng)十分接近，但是仍然存在約85 nT的平均絕對(duì)誤差，且隨著求解區(qū)域的擴(kuò)大，傳統(tǒng)有限單元法的精度增加有限，圖5中曲線5和曲線1也不重合。圖5中，解析解對(duì)應(yīng)的曲線1和有限元—無限元數(shù)值模擬對(duì)應(yīng)的曲線2幾乎完全重合，基于有限元—無限元算法所產(chǎn)生的絕對(duì)誤差在0左右，說明了基于有限元—無限元算法的準(zhǔn)確性。對(duì)比有限元—無限元(曲線2)和求解域邊長(zhǎng)為200 m的傳統(tǒng)有限元(曲線3)的數(shù)值模擬結(jié)果，在相同的區(qū)域內(nèi)，基于有限元—無限元的算法(曲線2)，在邊界處所產(chǎn)生的絕對(duì)誤差基本為0，基本避免了截?cái)嗾`差，大大改善了邊界處的計(jì)算結(jié)果，再次說明基于有限元—無限元算法相對(duì)于傳統(tǒng)有限元算法，在邊界處能夠得到更高的精度。

圖6 絕對(duì)誤差分布曲線Fig.6 The absolute error distribution curve

不同有限元算法模型的計(jì)算效率和誤差對(duì)比如表1所示?；谟邢拊獰o限元的算法需要對(duì)外部設(shè)置的無限元區(qū)域進(jìn)行網(wǎng)格剖分，所以相對(duì)初始模型會(huì)生成更多的網(wǎng)格節(jié)點(diǎn)。由表1可知，當(dāng)求解域邊界長(zhǎng)度設(shè)置為400 m時(shí)，傳統(tǒng)有限元算法的平均相對(duì)誤差為8.4%，均方根誤差為86.1 nT，最大絕對(duì)誤差為120.9 nT，相對(duì)于初始模型的計(jì)算精度有了很大的提升。但是其設(shè)置邊界長(zhǎng)度過大，求解區(qū)域的大小增大了近12倍，計(jì)算所需要的時(shí)間為240 s，占用了超過20 GB的內(nèi)存，若觀測(cè)平面范圍增大，在保證精度的情況下，求解域會(huì)進(jìn)一步增大，計(jì)算量也會(huì)進(jìn)一步增加，嚴(yán)重影響了有限元數(shù)值模擬的正演效率，大大增加了運(yùn)算成本。而基于有限元—無限元算法的平均相對(duì)誤差僅有0.78%，是使用擴(kuò)邊算法精度的10倍，占用內(nèi)存為4.85 GB，僅為擴(kuò)邊算法的 1/5，而運(yùn)算時(shí)間則提高了2.4倍。所以，綜合考慮剩磁和退磁的球體模型實(shí)驗(yàn)結(jié)果說明：基于有限元—無限元的算法，在相同的計(jì)算區(qū)域時(shí)，邊界處能夠達(dá)到更高的精度；同時(shí)也說明在相同的精度下，基于有限元—無限元算法與傳統(tǒng)有限元算法相比能夠在相對(duì)更小的區(qū)域內(nèi)達(dá)到精度要求，使用的節(jié)點(diǎn)數(shù)更少，占用內(nèi)存更少，大大降低了運(yùn)算成本，提高了正演模擬的精度。

表1 不同算法模型計(jì)算效率及誤差對(duì)比Table 1 Comparison about efficiency and relative error of different model

2.2 組合體模型

圖7為多個(gè)地質(zhì)體賦存的較復(fù)雜組合體模型，地磁場(chǎng)大小T0=50 000 nT，地磁傾角I0=45°，偏角A0=0°；球體中心坐標(biāo)為(60，0，-50)，半徑為30 m，剩余磁化強(qiáng)度大小為3.98 A/m，剩磁傾角Ir1=50°，偏角Ar1=60°，球體模型磁化率k1=0.4 SI；圓柱體沿y方向無限延伸，埋深為50 m，截面半徑為15 m，剩余磁化強(qiáng)度大小為7.96 A/m，剩磁傾角Ir2=70°，偏角Ar2=50°，磁化率k2=0.5 SI；棱柱體中心坐標(biāo)為(-60,0,-60)，長(zhǎng)100 m，寬40 m，高40 m，方位為西偏北10°，磁化率k3=0.1 SI，剩余磁化強(qiáng)度為0。球體的退磁系數(shù)N1=1/3，而對(duì)于任意橫截面形狀的無限長(zhǎng)柱體，當(dāng)與長(zhǎng)度垂直方向磁化時(shí)，退磁系數(shù)N2=1/2[32]，對(duì)球體和無限延伸圓柱體進(jìn)行退磁改正，棱柱體磁化率較小且沒有剩磁，則不進(jìn)行退磁改正。

圖7所示，模型求解域邊界長(zhǎng)度為200 m，觀測(cè)平面設(shè)置在z=0平面，為200 m×200 m的方形區(qū)域，x和y的坐標(biāo)范圍均為-100～100 m。測(cè)線位于觀測(cè)平面上，測(cè)點(diǎn)分布在y=-100～100 m，x=0的直線上。在相同的網(wǎng)格剖分情況下，分別使用基于有限元—無限元的算法和傳統(tǒng)有限單元法進(jìn)行計(jì)算。不同有限元算法數(shù)值模擬結(jié)果及平面誤差分布如圖8所示，計(jì)算效率和測(cè)線中的誤差對(duì)比如表2所示。

圖8中，兩種有限單元法的數(shù)值模擬結(jié)果都能較好地展示z=0平面的ΔT磁異常形態(tài)特征，但是數(shù)值模擬的結(jié)果有較大的差異。傳統(tǒng)三維有限元數(shù)值模擬共受到6個(gè)方向截?cái)噙吔绲挠绊?，?dāng)求解域邊界長(zhǎng)度為600 m時(shí)，即求解域邊界長(zhǎng)度為觀測(cè)區(qū)域邊長(zhǎng)的3倍，相對(duì)于初始模型求解范圍增大了26倍，使平均相對(duì)誤差降低到7.62%，但是其內(nèi)存占用增加了近7倍，計(jì)算耗時(shí)增加了近9倍，圖8d中，在z=0平面上仍存在近60 nT的絕對(duì)誤差，說明傳統(tǒng)算法提升計(jì)算精度能力差。而基于有限元—無限元算法，在邊界處的誤差為10 nT左右，平均相對(duì)誤差為1.68%，占用內(nèi)存僅為3.67 G。表明基于有限元—無限元算法能在占用更少運(yùn)算資源的情況下，達(dá)到遠(yuǎn)高于傳統(tǒng)有限單元法的計(jì)算精度，更適用于復(fù)雜模型體的三維磁場(chǎng)正演計(jì)算，進(jìn)一步說明了基于有限元—無限元算法的優(yōu)越性。

圖7 組合形體模型及網(wǎng)格剖分示意(紅色直線代表測(cè)線，黃色為觀測(cè)平面范圍)Fig.7 Diagram of combined model and mesh generation(the red line is observation line, the yellow plane is the range of observation plane)

表2 組合形體計(jì)算效率及誤差對(duì)比Table 2 Comparison about efficiency and relative error of different model

圖8 有限元數(shù)值模擬結(jié)果及平面誤差分布Fig.8 The map of the Finite element numerical simulation results and plane error distribution

在組合體模型的基礎(chǔ)上，將觀測(cè)面設(shè)置為起伏地表，如圖9所示，利用有限元—無限元方法進(jìn)行正演數(shù)值模擬，網(wǎng)格剖分參數(shù)和表2中一致。

圖10為模擬結(jié)果，圖10b中有限元—無限元耦合算法的結(jié)果與圖10a中的解析解的異常形態(tài)及幅值基本一致。數(shù)值模擬產(chǎn)生的最大絕對(duì)誤差為49.87 nT，平均相對(duì)誤差為3.04%，相對(duì)于觀測(cè)面為z=0平面時(shí)誤差有所增加，但仍處于較低的水平。說明基于有限元—無限元耦合算法適用于起伏地形情況，表明了有限元—無限元耦合算法的靈活性。

3 結(jié)論及討論

1) 在進(jìn)行磁場(chǎng)三維正演模擬時(shí)，有限單元法會(huì)受到6個(gè)方向的截?cái)噙吔缧?yīng)的影響，傳統(tǒng)算法在擴(kuò)大邊界時(shí)會(huì)成倍的增加求解區(qū)域，占用大量運(yùn)算資源，而計(jì)算精度的提升有限。在進(jìn)行復(fù)雜條件下、復(fù)雜形體的正演模擬時(shí)，應(yīng)該充分使用基于有限元—無限元的算法，以充分利用計(jì)算資源，提高正演效率。

圖9 起伏地表模型(黃色為觀測(cè)面范圍)Fig.9 Diagram of rugged surface model (The yellow plane is the range of observation plane)

圖10 起伏地形下有限元數(shù)值模擬結(jié)果及平面誤差分布Fig.10 The map of the Finite element numerical simulation results and plane error distribution in rugged surface

2) 通過孤立球體模型和組合模型的數(shù)值模擬，在相同測(cè)區(qū)的計(jì)算范圍內(nèi)，有限元—無限元算法與傳統(tǒng)的有限元算法相比，能夠達(dá)到更高的計(jì)算精度。尤其是距離邊界較近時(shí)，傳統(tǒng)有限元方法誤差較大，而在本文設(shè)置的模型中，基于有限元—無限元的算法能夠達(dá)到2%以內(nèi)的相對(duì)誤差。

3) 傳統(tǒng)的有限元方法為達(dá)到較高的精度，需要進(jìn)行擴(kuò)邊處理，設(shè)置較大的求解區(qū)域范圍，且隨著求解域范圍的擴(kuò)大，精度提升有限。而基于有限元—無限元算法，只需與計(jì)算區(qū)域相當(dāng)?shù)木W(wǎng)格剖分范圍就可滿足精度要求。相對(duì)于傳統(tǒng)的有限元算法，有限元—無限元算法能在占用更少運(yùn)算資源的情況下，達(dá)到遠(yuǎn)高于傳統(tǒng)有限單元法的計(jì)算精度，正演效率更高。

4) 在起伏地形下，基于有限元—無限元耦合算法仍能高精度的完成數(shù)值模擬，進(jìn)一步說明了有限元—無限元算法的靈活性，同時(shí)也表明該算法在復(fù)雜地質(zhì)條件下有較好的應(yīng)用前景。

致謝：感謝中國(guó)地質(zhì)大學(xué)(武漢)馬火林副教授、朱丹博士后對(duì)本文提出的指導(dǎo)意見。