基于非結構網格的三種CSEM有限元三維正演系統分析

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

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

0 引言

可控源電磁法(Controlled Source Electromagnetic Method，CSEM)具有抗干擾能力強、工作效率高、測量精度高等諸多優點，是研究中淺部電性結構的主要地球物理方法之一，廣泛應用于固體礦產、水文、石油天然氣普查、地熱田勘探、環境地質調查及環境與工程地球物理勘查等領域[1-3]。目前，實際勘探任務大多以復雜三維地質構造為主，傳統的一維、二維假設條件下的資料解釋技術不能滿足目前的地質解釋要求，迫切需要開發三維可控源電磁資料解釋技術。正演是反演的基礎，因此國內外學者都致力于任意復雜地形的地電結構CSEM三維正演模擬研究，尋求快速、高精度的三維正演方法。

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

1 基本理論

1.1 CSEM偏微分方程

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

(1)

(2)

(3)

(4)

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

(5)

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

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

電場E和磁感應強度B可用矢量位A和標量位Φ分別表示為

(6)

(7)

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

(8)

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

(9)

聯立式(8)和式(9)，可得基于雙旋度結構的A-Φ耦合勢可控源電磁方程

(10)

(11)

最后形成系統方程

(12)

式(12)的求解系統采用節點基函數對矢量位A和標量位Φ進行離散，可獲得拉普拉斯結構的A-Φ可控源電磁法求解方程。另外，本文通過擴邊加載截斷邊界條件，實現電磁場在無窮遠處消失為零。

1.2 有限元分析

以式(5)為控制方程，應用Galerkin法推導有限元方程。設余量為

(13)

令re在計算區域T滿足

?TN·redT=0

(14)

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

(15)

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

(16)

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

同理，對式(10)應用Galerkin法進行推導并化簡，最終的有限元積分方程為

(17)

(18)

式中L表示有限元中的節點基函數。

同理，對式(12)應用Galerkin法進行公式推導及化簡，得到最終的有限元積分表達式為

(19)

(20)

(21)

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

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

(22)

式中矢量基函數Ni可用節點基函數[29]表示為

(23)

式中：li為棱邊i的長度；Lj1與Lj2表示該棱邊兩個節點的基函數，具體標號如圖1所示。

圖1 四面體單元棱邊和節點關系1～4表示節點編號，①～⑥表示棱邊編號

(2)對于雙旋度結構的A-Φ耦合勢CSEM方程，矢量位A和標量位Φ可分別采用矢量基函數和節點基函數進行描述

(24)

(3)拉普拉斯結構的A-Φ耦合勢CSEM方程，矢量位A和標量位Φ可采用節點基函數表示為

(25)

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

1.3 場源加載

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

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

(26)

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

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

(27)

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

垂直磁偶源M的表達式為

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

(28)

(29)

(30)

由上式可得

(31)

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

2 算例分析

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

2.1 模型1：均勻半空間

2.1.1 算法驗證

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

圖2 垂直磁偶源網格剖分示意圖

計算3Hz的垂直磁偶源電磁響應。整個求解區域被剖分成403789個四面體單元，471329條邊和66521個節點。沿x方向等間距設置觀測點，點距為120m，共50個觀測點。

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

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

圖3 3Hz磁場Hz實部(上)和虛部(下)響應曲線

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

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

圖5 100Hz磁場Hz的實部(上)和虛部(下)響應曲線

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

2.1.2 求解系統特點分析

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

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

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

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

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

2.2 模型2：低阻模型

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

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

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

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

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

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

表2 三種求解系統相關技術性能對比

3 結論

本文采用非結構化網格對求解區域進行離散，實現了三種求解系統的CSEM正演計算。建立均勻半空間模型和低阻異常體地電模型，從計算精度、效率以及內存消耗三個方面分析基于兩種A-Φ耦合勢方程和電場方程正演問題的求解特性。得出以下幾點認識。

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

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

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

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