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

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

0 引言

可控源電磁法(Controlled Source Electromagnetic Method，CSEM)具有抗干擾能力強(qiáng)、工作效率高、測量精度高等諸多優(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ǔ)，因此國內(nèi)外學(xué)者都致力于任意復(fù)雜地形的地電結(jié)構(gòu)CSEM三維正演模擬研究，尋求快速、高精度的三維正演方法。

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

1 基本理論

1.1 CSEM偏微分方程

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

(1)

(2)

(3)

(4)

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

(5)

上式即是三維CSEM需滿足的電場方程。

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

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

(6)

(7)

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

(8)

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

(9)

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

(10)

(11)

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

(12)

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

1.2 有限元分析

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

(13)

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

?TN·redT=0

(14)

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

(15)

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

(16)

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

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

(17)

(18)

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

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

(19)

(20)

(21)

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

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

(22)

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

(23)

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

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

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

(24)

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

(25)

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

1.3 場源加載

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

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

(26)

對于x方向電性源來說，源線段置于四面體棱邊上，因此有

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

(27)

式中：Θ表示Heaviside函數(shù)；δ是delta函數(shù)；Nx、Ny、Nz分別表示矢量基函數(shù)N的三個分量；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)

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

2 算例分析

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

2.1 模型1：均勻半空間

2.1.1 算法驗(yàn)證

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

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

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

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

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

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

圖4 3Hz磁場Hz的收斂曲線

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

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

圖6 100Hz磁場Hz的收斂曲線

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

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

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

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

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

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

2.2 模型2：低阻模型

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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

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