基于非結(jié)構(gòu)網(wǎng)格的三種CSEM有限元三維正演系統(tǒng)分析

周 峰 張志勇 陳 煌 湯井田 鄧居智 李 勇

(①東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，江西南昌 330013；②東華理工大學(xué)地球物理與測(cè)控技術(shù)學(xué)院，江西南昌 330013；③中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南長(zhǎng)沙 410083；④中國(guó)地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所自然資源部地球物理電磁法探測(cè)技術(shù)重點(diǎn)實(shí)驗(yàn)室，河北廊坊 065000)

0 引言

可控源電磁法(Controlled Source Electromagnetic Method，CSEM)具有抗干擾能力強(qiáng)、工作效率高、測(cè)量精度高等諸多優(yōu)點(diǎn)，是研究中淺部電性結(jié)構(gòu)的主要地球物理方法之一，廣泛應(yīng)用于固體礦產(chǎn)、水文、石油天然氣普查、地?zé)崽锟碧?、環(huán)境地質(zhì)調(diào)查及環(huán)境與工程地球物理勘查等領(lǐng)域[1-3]。目前，實(shí)際勘探任務(wù)大多以復(fù)雜三維地質(zhì)構(gòu)造為主，傳統(tǒng)的一維、二維假設(shè)條件下的資料解釋技術(shù)不能滿足目前的地質(zhì)解釋要求，迫切需要開發(fā)三維可控源電磁資料解釋技術(shù)。正演是反演的基礎(chǔ)，因此國(guó)內(nèi)外學(xué)者都致力于任意復(fù)雜地形的地電結(jié)構(gòu)CSEM三維正演模擬研究，尋求快速、高精度的三維正演方法。

針對(duì)CSEM三維正演問(wèn)題，學(xué)者們做了大量研究工作。首先針對(duì)場(chǎng)源奇異性導(dǎo)致正演精度損失的問(wèn)題，將CSEM正演分解為總場(chǎng)算法和二次場(chǎng)算法。二次場(chǎng)算法能夠去除場(chǎng)源奇異性，但無(wú)法實(shí)現(xiàn)任意復(fù)雜地形條件下地電模型的正演模擬[4-9]。基于總場(chǎng)算法的求解方案早期主要是以偽delta函數(shù)刻畫場(chǎng)源問(wèn)題，一定程度上降低了場(chǎng)源奇異性，但無(wú)法從根本上解決場(chǎng)源的奇異性問(wèn)題[10-12]，甚至?xí)觿?chǎng)源加載的難度。為此，得益于近年來(lái)非結(jié)構(gòu)化網(wǎng)格離散技術(shù)的發(fā)展，廣泛應(yīng)用局部加密技術(shù)降低場(chǎng)源奇異性，結(jié)合delta函數(shù)場(chǎng)源積分能夠避免場(chǎng)源加載的困難[13-17]。隨后，基于自適應(yīng)的有限元正演技術(shù)被廣泛應(yīng)用于CSEM三維正演模擬，進(jìn)一步彌補(bǔ)了場(chǎng)源奇異性帶來(lái)的精度損失[14,18]。在實(shí)現(xiàn)大規(guī)模的電磁快速正演求解方面，早期基于節(jié)點(diǎn)有限元結(jié)合散度矯正的技術(shù)被應(yīng)用，雖然彌補(bǔ)了電場(chǎng)方程迭代求解不收斂現(xiàn)象，但同時(shí)降低了正演求解的精度[7]。之后，基于電場(chǎng)雙旋度方程的矢量有限元直接求解技術(shù)被提出，雖大幅度提高了正演方程的求解精度[13,16,19]，但無(wú)法避免雙旋度結(jié)構(gòu)帶來(lái)的空解問(wèn)題，使常規(guī)Krylov子空間迭代求解困難，導(dǎo)致直接求解的三維反演算法對(duì)內(nèi)存的需求更高，限制了該算法的應(yīng)用。為了使電場(chǎng)雙旋度方程可用于迭代計(jì)算，Schwarzbach等[14]開展了E-Φ方程的CSEM正演模擬(E和Φ分別表示電場(chǎng)和電標(biāo)量位)，取得了一定效果，但亦無(wú)法從根本上解決迭代收斂慢等問(wèn)題。之后，有學(xué)者提出基于勢(shì)位(如A-Φ勢(shì)[13,20-25]，A表示磁矢量位)的三維電磁正演求解系統(tǒng)，在一定程度上擴(kuò)展了三維電磁求解系統(tǒng)的種類，計(jì)算精度和效率得到一定的改善。

1 基本理論

1.1 CSEM偏微分方程

在準(zhǔn)靜態(tài)條件下，可控源電磁法滿足Maxwell方程組，取時(shí)間因子e-iωt，其頻率域表達(dá)式為

(1)

(2)

(3)

(4)

為消除磁場(chǎng)參數(shù)，將式(1)代入式(2)，可得

(5)

上式即是三維CSEM需滿足的電場(chǎng)方程。

有限元技術(shù)的基函數(shù)經(jīng)歷了節(jié)點(diǎn)基函數(shù)到矢量基函數(shù)的發(fā)展。為了消除偽解并遵循電場(chǎng)法向不連續(xù)性特征，基于矢量基函數(shù)的有限元技術(shù)求解電場(chǎng)方程成為目前主流方法。但是，可控源電磁法中的梯度電場(chǎng)是由地表和地下電導(dǎo)率不連續(xù)界面上的積累電荷產(chǎn)生的，導(dǎo)致基于電場(chǎng)方程的矢量有限元法很難實(shí)現(xiàn)梯度電場(chǎng)的模擬，而該梯度電場(chǎng)對(duì)于可控源電磁法非常重要，這一矛盾被稱為矢量有限元的空解難題。因此，有學(xué)者將目光聚焦到A-Φ耦合勢(shì)方程，進(jìn)行了較系統(tǒng)的研究。

電場(chǎng)E和磁感應(yīng)強(qiáng)度B可用矢量位A和標(biāo)量位Φ分別表示為

(6)

(7)

將式(6)和式(7)代入式(5)可得A-Φ耦合勢(shì)方程

(8)

根據(jù)電荷守恒以及歐姆定律，可得輔助方程

(9)

聯(lián)立式(8)和式(9)，可得基于雙旋度結(jié)構(gòu)的A-Φ耦合勢(shì)可控源電磁方程

(10)

(11)

最后形成系統(tǒng)方程

(12)

式(12)的求解系統(tǒng)采用節(jié)點(diǎn)基函數(shù)對(duì)矢量位A和標(biāo)量位Φ進(jìn)行離散，可獲得拉普拉斯結(jié)構(gòu)的A-Φ可控源電磁法求解方程。另外，本文通過(guò)擴(kuò)邊加載截?cái)噙吔鐥l件，實(shí)現(xiàn)電磁場(chǎng)在無(wú)窮遠(yuǎn)處消失為零。

1.2 有限元分析

以式(5)為控制方程，應(yīng)用Galerkin法推導(dǎo)有限元方程。設(shè)余量為

(13)

令re在計(jì)算區(qū)域T滿足

?TN·redT=0

(14)

式中N為矢量基函數(shù)。將式(13)代入式(14)，可得

(15)

對(duì)式(15)應(yīng)用第一矢量格林定理，可將其簡(jiǎn)化為

(16)

式中：Γ=?T表示求解區(qū)域外邊界；n為外邊界單位法向量。另外，擴(kuò)邊處理使得求解區(qū)域足夠大，因此式(16)中的面積分項(xiàng)可忽略不計(jì)。

同理，對(duì)式(10)應(yīng)用Galerkin法進(jìn)行推導(dǎo)并化簡(jiǎn)，最終的有限元積分方程為

(17)

(18)

式中L表示有限元中的節(jié)點(diǎn)基函數(shù)。

同理，對(duì)式(12)應(yīng)用Galerkin法進(jìn)行公式推導(dǎo)及化簡(jiǎn)，得到最終的有限元積分表達(dá)式為

(19)

(20)

(21)

下面分別敘述如何選擇上述三類方程有限元基函數(shù)。

(1)對(duì)電場(chǎng)方程采用矢量有限元進(jìn)行離散。對(duì)于每一個(gè)四面體單元e，任意點(diǎn)(x,y,z)處的電場(chǎng)可由各條棱邊的電場(chǎng)Ei和矢量基函數(shù)Ni表示為

(22)

式中矢量基函數(shù)Ni可用節(jié)點(diǎn)基函數(shù)[29]表示為

(23)

式中：li為棱邊i的長(zhǎng)度；Lj1與Lj2表示該棱邊兩個(gè)節(jié)點(diǎn)的基函數(shù)，具體標(biāo)號(hào)如圖1所示。

圖1 四面體單元棱邊和節(jié)點(diǎn)關(guān)系1～4表示節(jié)點(diǎn)編號(hào)，①～⑥表示棱邊編號(hào)

(2)對(duì)于雙旋度結(jié)構(gòu)的A-Φ耦合勢(shì)CSEM方程，矢量位A和標(biāo)量位Φ可分別采用矢量基函數(shù)和節(jié)點(diǎn)基函數(shù)進(jìn)行描述

(24)

(3)拉普拉斯結(jié)構(gòu)的A-Φ耦合勢(shì)CSEM方程，矢量位A和標(biāo)量位Φ可采用節(jié)點(diǎn)基函數(shù)表示為

(25)

因此，通過(guò)以上不同的離散形式可得到三種不同的CSEM求解系統(tǒng)，即電場(chǎng)方程、雙旋度結(jié)構(gòu)的A-Φ耦合勢(shì)和拉普拉斯結(jié)構(gòu)的A-Φ耦合勢(shì)，本文結(jié)合穩(wěn)健的直接求解器PARDISO和Krylov子空間算法對(duì)這三種求解系統(tǒng)進(jìn)行對(duì)比分析[27,30]。

1.3 場(chǎng)源加載

對(duì)于CSEM三維正演來(lái)說(shuō)，源項(xiàng)加載不準(zhǔn)確和場(chǎng)源的奇異性會(huì)嚴(yán)重影響正演的求解精度。針對(duì)這兩種問(wèn)題，常用的場(chǎng)源處理技術(shù)主要包含兩類：①二次場(chǎng)算法，其目的是解決源的奇異性問(wèn)題，但需要計(jì)算具有解析表達(dá)式的一次場(chǎng)，因而對(duì)復(fù)雜地電模型模擬的局限性較大；②基于總場(chǎng)算法并結(jié)合偶極源積分實(shí)現(xiàn)場(chǎng)源加載，具有較好的適用性，但仍無(wú)法避免場(chǎng)源奇異性問(wèn)題。為此，本文采用偶極源積分結(jié)合非結(jié)構(gòu)化網(wǎng)格局部細(xì)化技術(shù)，不僅能有效降低場(chǎng)源奇異性，還能實(shí)現(xiàn)場(chǎng)源的統(tǒng)一表示形式。

對(duì)于電場(chǎng)方程的三維CSEM有限元系統(tǒng)，場(chǎng)源加載可表示為

(26)

對(duì)于x方向電性源來(lái)說(shuō)，源線段置于四面體棱邊上，因此有

?T(Nxix+Nyiy+Nziz)·[Θ(xi+1)-Θ(xi)]δ(y-y0)δ(z-z0)ixdT

(27)

式中：Θ表示Heaviside函數(shù)；δ是delta函數(shù)；Nx、Ny、Nz分別表示矢量基函數(shù)N的三個(gè)分量；ix、iy、iz表示笛卡爾坐標(biāo)系下沿坐標(biāo)軸方向的單位向量；(x0,y0,z0)表示偶極源的中心坐標(biāo)。

垂直磁偶源M的表達(dá)式為

M=mδ(x-x0)δ(y-y0)δ(z-z0)iz

(28)

(29)

(30)

由上式可得

(31)

將場(chǎng)源置于四面體單元內(nèi)，上式中的面積分項(xiàng)(公式右邊第二項(xiàng))等于零，因此積分式只剩下右邊第一項(xiàng)。若采用電性源的形式加載磁性源，場(chǎng)源的積分公式與電性源一樣，在此不再敘述。

2 算例分析

為了分析上述三種求解系統(tǒng)對(duì)三維CSEM的模擬效果，設(shè)計(jì)了兩個(gè)模型進(jìn)行正演計(jì)算。一個(gè)模型是存在解析解的均勻半空間，另外一個(gè)是無(wú)解析解的低阻地電模型。測(cè)試本文開發(fā)算法的計(jì)算機(jī)配置為：內(nèi)存32G，intel(R)CoreTMi7-9700 CPU@ 3.00GHz。

2.1 模型1：均勻半空間

2.1.1 算法驗(yàn)證

對(duì)均勻半空間模型加載一個(gè)正方形的垂直磁偶源，磁偶源是大小為1.0m×1.0m的線框，位于求解區(qū)域的中心。地下背景電阻率為100Ω·m，空氣電阻率為108Ω·m。求解區(qū)域大小為[-30km,30km]3。為了降低場(chǎng)源的奇異性，對(duì)場(chǎng)源進(jìn)行局部加密，被離散為27個(gè)小的線段，如圖2所示。

圖2 垂直磁偶源網(wǎng)格剖分示意圖

計(jì)算3Hz的垂直磁偶源電磁響應(yīng)。整個(gè)求解區(qū)域被剖分成403789個(gè)四面體單元，471329條邊和66521個(gè)節(jié)點(diǎn)。沿x方向等間距設(shè)置觀測(cè)點(diǎn)，點(diǎn)距為120m，共50個(gè)觀測(cè)點(diǎn)。

采用Krylov子空間的GMRES算法結(jié)合ILU(0)預(yù)條件因子分別對(duì)前述兩種結(jié)構(gòu)的A-Φ耦合勢(shì)有限元方程進(jìn)行求解，采用PARDISO直接求解器對(duì)電場(chǎng)方程進(jìn)行求解，得到磁場(chǎng)垂直分量Hz的實(shí)部和虛部，并與解析解進(jìn)行對(duì)比(圖3)。兩種結(jié)構(gòu)的A-Φ耦合勢(shì)方程組求解的收斂曲線見(jiàn)圖4。

從圖3可知，三種求解系統(tǒng)的結(jié)果與解析解都高度吻合，誤差基本都小于1.5%。從圖4的收斂曲線可知，雙旋度結(jié)構(gòu)的A-Φ耦合勢(shì)方程組求解需要迭代3000次，相對(duì)殘差為2.31×10-15，而拉普拉斯結(jié)構(gòu)的A-Φ耦合勢(shì)方程組求解僅需200次，相對(duì)殘差為1.9×10-15。

圖3 3Hz磁場(chǎng)Hz實(shí)部(上)和虛部(下)響應(yīng)曲線

圖4 3Hz磁場(chǎng)Hz的收斂曲線

本文還計(jì)算了100Hz時(shí)垂直磁偶源的電磁響應(yīng)。求解區(qū)域被離散為814255個(gè)四面體單元、675519條邊單元和138736個(gè)節(jié)點(diǎn)，并沿x方向等間距設(shè)置觀測(cè)點(diǎn)，間距為20m，共50個(gè)觀測(cè)點(diǎn)。圖5展示了100Hz時(shí)的電磁響應(yīng)。兩種結(jié)構(gòu)的A-Φ耦合勢(shì)收斂曲線如圖6所示。從圖6可知，雙旋度結(jié)構(gòu)的A-Φ耦合勢(shì)迭代3000次未達(dá)到預(yù)設(shè)的求解精度，而拉普拉斯結(jié)構(gòu)的A-Φ耦合勢(shì)經(jīng)迭代2400次即達(dá)到預(yù)設(shè)精度。雖然前者得到的垂直磁場(chǎng)Hz與解析解吻合較好，但求解的迭代次數(shù)較多。對(duì)比分析表明，本文開發(fā)的三種CSEM求解方程算法準(zhǔn)確，結(jié)果可靠。

圖5 100Hz磁場(chǎng)Hz的實(shí)部(上)和虛部(下)響應(yīng)曲線

圖6 100Hz磁場(chǎng)Hz的收斂曲線

2.1.2 求解系統(tǒng)特點(diǎn)分析

為了分析三種CSEM求解方程的特點(diǎn)，仍然采用上節(jié)的均勻半空間地電模型。沿x方向、在地面放置一個(gè)電偶源，源長(zhǎng)度為1m，發(fā)射電流為1A，發(fā)射頻率為10Hz。求解區(qū)域的大小為[-30km,30km]3，分別采用前述三種求解方程對(duì)該模型進(jìn)行求解，地面上各場(chǎng)分量等值線圖見(jiàn)圖7～圖9。

由圖7～圖9可以看出，基于雙旋度結(jié)構(gòu)的A-Φ耦合勢(shì)計(jì)算的矢量和標(biāo)量位無(wú)明顯規(guī)律，基于拉普拉斯結(jié)構(gòu)的A-Φ耦合勢(shì)計(jì)算的矢量和標(biāo)量位特征符合預(yù)期。通過(guò)兩種A-Φ耦合勢(shì)計(jì)算的電場(chǎng)Ex等值線與基于電場(chǎng)方程計(jì)算的電場(chǎng)Ex等值線具有高度一致性，表明矢量基函數(shù)離散的矢量位A違背了連續(xù)性條件，導(dǎo)致矢量位A切向不連續(xù)，使得計(jì)算得到的矢量和標(biāo)量位無(wú)明顯規(guī)律，影響其收斂性。

圖7 均勻半空間模型基于雙旋度結(jié)構(gòu)的A-Φ耦合勢(shì)矢量位A各分量及電場(chǎng)分量Ex平面等值線圖

圖8 均勻半空間模型基于拉普拉斯結(jié)構(gòu)的A-Φ耦合勢(shì)矢量位A各分量及電場(chǎng)分量Ex平面等值線圖

圖9 均勻半空間模型基于電場(chǎng)方程求解的電場(chǎng)E和磁場(chǎng)H各分量平面等值線分布

2.2 模型2：低阻模型

圖10 低阻模型剖面(左)及模型剖分方案(右)

圖11 低阻模型基于三種求解系統(tǒng)得到的3Hz電磁場(chǎng)分量實(shí)部(左)和虛部(右)對(duì)比

圖12 低阻模型基于三種求解系統(tǒng)得到的10Hz電磁場(chǎng)分量實(shí)部(左)和虛部(右)對(duì)比

3Hz所對(duì)應(yīng)的基于三種求解系統(tǒng)的不同求解器以及預(yù)條件因子的收斂性能統(tǒng)計(jì)結(jié)果見(jiàn)表1。由表可知，基于電場(chǎng)方程的CSEM矢量有限元方程在常規(guī)迭代求解器和預(yù)條件因子下均不能得到滿意的計(jì)算結(jié)果。在同等精度、迭代求解器條件下，基于拉普拉斯結(jié)構(gòu)的A-Φ耦合勢(shì)CSEM求解方程收斂性總體上優(yōu)于基于電場(chǎng)方程和雙旋度結(jié)構(gòu)的A-Φ耦合勢(shì)方程。對(duì)于同一種求解方程，GMRES和BICGSTAB迭代求解器的求解效果相差不大。除此之外，GMRES和BICGSTAB迭代求解器結(jié)合SOR(超松弛)預(yù)條件因子對(duì)兩種結(jié)構(gòu)的A-Φ耦合勢(shì)CSEM求解方程進(jìn)行求解，收斂性差異不明顯；而GMRES和BICGSTAB迭代求解器結(jié)合其他預(yù)條件因子進(jìn)行求解的收斂性明顯變差，所需迭代次數(shù)明顯增加。這說(shuō)明選擇合適的預(yù)條件因子可提高方程組的求解效率。

表1 3Hz時(shí)三種求解系統(tǒng)、不同求解器及預(yù)條件因子的求解收斂性能對(duì)比

表2是不同求解方程/求解器組合下CSEM求解內(nèi)存消耗、求解時(shí)間、迭代次數(shù)及相對(duì)殘差的統(tǒng)計(jì)。對(duì)于方程組未知數(shù)在百萬(wàn)級(jí)以內(nèi)的情況來(lái)說(shuō)，直接求解器求解電場(chǎng)方程的優(yōu)勢(shì)明顯，小型計(jì)算機(jī)即可滿足內(nèi)存需求；若需求解的未知數(shù)超過(guò)100萬(wàn)，計(jì)算機(jī)內(nèi)存需求則超過(guò)15GB；當(dāng)需求解的未知數(shù)繼續(xù)增加，利用直接求解器求解線性方程組的計(jì)算機(jī)內(nèi)存需求則難以滿足。對(duì)于A-Φ耦合勢(shì)求解方程，當(dāng)未知數(shù)超過(guò)100萬(wàn)，內(nèi)存需求不超過(guò)6GB，即A-Φ耦合勢(shì)方程的迭代解法在內(nèi)存需求方面的優(yōu)勢(shì)非常明顯。但是，在時(shí)間消耗方面，對(duì)于未知數(shù)不超過(guò)100萬(wàn)的方程組求解，直接求解器相比迭代求解器優(yōu)勢(shì)更明顯。

表2 三種求解系統(tǒng)相關(guān)技術(shù)性能對(duì)比

3 結(jié)論

本文采用非結(jié)構(gòu)化網(wǎng)格對(duì)求解區(qū)域進(jìn)行離散，實(shí)現(xiàn)了三種求解系統(tǒng)的CSEM正演計(jì)算。建立均勻半空間模型和低阻異常體地電模型，從計(jì)算精度、效率以及內(nèi)存消耗三個(gè)方面分析基于兩種A-Φ耦合勢(shì)方程和電場(chǎng)方程正演問(wèn)題的求解特性。得出以下幾點(diǎn)認(rèn)識(shí)。

(1)本文開發(fā)的三種正演求解算法程序準(zhǔn)確且計(jì)算精度較高，采用的統(tǒng)一場(chǎng)源處理技術(shù)能夠較好地?cái)U(kuò)展場(chǎng)源的適用性，同時(shí)結(jié)合局部網(wǎng)格加密技術(shù)可有效降低場(chǎng)源奇異性問(wèn)題。

(2)在同等條件下，基于非結(jié)構(gòu)網(wǎng)格的矢量有限元電場(chǎng)方程不適合于采用常規(guī)的Krylov子空間算法進(jìn)行方程組的求解；相比于雙旋度結(jié)構(gòu)，基于拉普拉斯結(jié)構(gòu)的A-Φ耦合勢(shì)更適合于電磁場(chǎng)的迭代求解，其適用性更強(qiáng)。

(3)相比于A-Φ耦合勢(shì)方程，基于電場(chǎng)方程的直接解法在方程組未知數(shù)不高于100萬(wàn)時(shí)，在內(nèi)存需求和求解效率方面具有明顯的優(yōu)勢(shì)；但是，若方程組未知數(shù)超過(guò)100萬(wàn)級(jí)，相比直接解法，A-Φ耦合勢(shì)方程迭代算法在內(nèi)存需求方面優(yōu)勢(shì)更明顯。另外，若網(wǎng)格密度增加，常規(guī)的Krylov子空間迭代算法求解大型方程組所需要的時(shí)間成本會(huì)逐漸增加。

總之，本文系統(tǒng)對(duì)比了三種頻率域CSEM正演求解系統(tǒng)，并對(duì)正演模擬的內(nèi)存需求及效率進(jìn)行了探討。對(duì)于小尺度三維CSEM正演問(wèn)題，電場(chǎng)方程的直接解法能滿足需求；對(duì)于超大尺度三維CSEM正演問(wèn)題，需要借助迭代算法開展線性方程組求解。