非結構化網格與散度校正技術結合的三維CSEM矢量有限元正演

湯文武 鄧居智 黃清華

(①北京大學地球與空間科學學院，北京 100871；②東華理工大學核資源與環境國家重點實驗室，江西南昌 330013)

0 引言

電磁勘探方法分支眾多，在礦產資源勘查、油氣勘探、地質環境調查等方面得到廣泛應用[1-3]。可控源電磁法(CSEM)作為其中很重要的一個分支，其應用效果一定程度上依賴于高精度的三維正反演解釋技術，而三維正演效率和求解精度的提高可為三維高效反演奠定基礎[4]。目前，用于三維電磁正演的主流數值模擬方法包括積分方程法[5-7]、有限差分法[8-11]、有限體積法[12-14]和有限單元法[15-16]。

有限單元法因其模擬網格的多樣性，可精確模擬地形起伏、不規則異常體形態而備受關注。國內外諸多學者采用非結構網格開展了三維電磁有限元正演模擬：Puzyrev等[17]基于二次耦合勢和非結構化四面體網格開展了三維電磁正演模擬研究，通過引入一種TAI(剛度矩陣的近似截斷逆矩陣)預條件矩陣加速了迭代求解過程；Ansari等[18]基于二次耦合勢的雙旋度方程開展了非結構化四面體網格的三維電磁正演，通過ILUT 預條件的QMR(Quasi-Minimum Residual)方法求解線性方程組；Cai等[19]基于總場和非結構化網格開展了三維電磁矢量有限元正演，通過直接求解器MUMPS對線性方程組進行求解，這種方法適用于多源電磁問題；劉長生等[20]、楊軍等[21]、殷長春等[22]、劉穎等[23]、葉益信等[24]及趙寧等[25]也開展了基于非結構化網格的三維電磁正演研究。

基于非結構化網格對電場直接求解時，其有限元正演方程的求解絕大多數是基于直接法。究其原因，一方面是非結構化網格相比結構化網格得到的線性方程組條件數更大，另一方面對于較低頻率，基于電場方程迭代求解時電流密度散度不嚴格為零。因此，利用迭代法求解存在收斂慢甚至不收斂的問題；而基于直接法求解時，對內存需求大，在個人PC上不易實現大規模正演計算。當采用結構化網格進行正演時，往往可通過散度校正改善低頻正演的迭代收斂速度[26-29],但目前未見散度校正技術在非結構化網格正演中的應用。秦策等[30]利用分塊對角預條件并結合輔助空間預條件實現了高效率迭代解法，可將迭代求解次數降低到幾十以內，提高了正演計算效率。但該方法實現較復雜，需借助開源程序庫實現輔助空間預條件算法。

本文結合非結構化四面體網格和散度校正技術開展三維CSEM正演模擬研究，為非結構三維正演的迭代求解提供一條路徑。首先介紹基于二次電場的正演模擬方法及散度校正技術，之后分別通過一個層狀介質模型和兩個三維地電模型的數值模擬驗證算法的有效性及高精度。

1 方法原理

1.1 三維電磁正演邊值問題

對于頻率小于1×105Hz的電磁問題，從麥克斯韋方程出發(準靜態極限)，在頻率域采用負諧變形式(e-iω t)，則電磁場滿足

(1)

(2)

式中：E表示電場矢量;H為磁場矢量；B表示磁場強度矢量;ω表示角頻率；μ0表示真空中的磁導率；σ表示電導率參數；J外表示外電流密度。對式(1)兩邊求旋度，得到

××E=iωμ0×H

(3)

將式(2)代入式(3)，整理即可得到關于電場E的雙旋度方程

××E-iωμ0σE=iωμ0J外

(4)

為了避免場源的奇異性，兼顧網格劃分的便利性，采用一次場與二次場分離辦法，即總場E分解為一次場Ep及二次場Es

E=Ep+Es

(5)

對于基于參考模型電導率分布σp的一次場Ep，其雙旋度方程為

××Ep-iωμ0σpEp=iωμ0J外

(6)

參考模型一般選為均勻半空間或水平層狀介質模型，Ep可預先通過快速漢克爾變換計算得到。用式(4)減去式(6)，并結合式(5)，整理可得關于二次電場的雙旋度方程

××Es-iωμ0σEs=iωμ0ΔσEp

(7)

式中Δσ=σ-σp，表示異常電導率的分布。當研究區域足夠大，以至于邊界上的二次場可忽略時，可采用狄利克雷邊界條件，得到邊值問題[31]

(8)

式中：Ω表示包含地下介質和空氣層的長方體研究區域；Γ表示其外邊界；n為兩個相鄰單元界面的法向單位向量。

1.2 基于伽遼金有限元法的邊值問題求解

為了適應任意地形起伏及復雜地質目標模擬的需要，采用伽遼金有限元法對式(8)邊值問題進行求解。圖1所示為研究區域示意圖，研究區域包含空氣層及地下介質層。同時，為了在電導率分界面上滿足電場法向分量不連續的邊界條件，采用基于棱邊的矢量有限元，并采用四面體網格對研究區域進行剖分(圖2)。

圖1 研究區域示意圖

圖2 四面體剖分單元示意圖

設對應各條棱邊的形函數為Ni，那么一次電場及二次電場可分別表示為

(9)

(10)

式中Epi和Esi分別代表第i條棱邊上的一次場和二次場值，i=1～6。本文選用矢量基函數

Ni=(Li1Li2-Li2Li1)li

(11)

式中：i1和i2分別表示第i條棱邊的兩個節點編號；li代表該棱邊的長度；L為對應于各節點的形函數，L為L的梯度[32]。根據伽遼金有限元方法，可在式(8)兩邊同時乘以形函數Ni并對整個研究區域Ω積分

(12)

用LHS表示上式左端項，可分解為兩項，即

(13)

根據矢量恒等式

-A·(×B)+B·(×A)=·(A×B)

式(13)右邊的第一項可作如下變換

×Es·(×Ni)]dΩ

s×Es×Ni·ds

(14)

式中S表示面積矢量。根據二次電場和二次磁場的關系，當邊界離異常區域足夠遠時，滿足

s×Es×Ni·ds=siωμ0Hs×Ni·ds=0

(15)

結合式(13)～式(15)，式(12)可表示為

(16)

式(16)可表示為線性方程組

KEs=MEp

(17)

式中：K代表式(17)左邊作用于Es的剛度矩陣；M代表作用于Ep的矩陣。K和M的矩陣元素分別為

式(17)中的剛度矩陣K為大型、稀疏、復對稱矩陣。為了降低三維正演模擬的內存需求，采用Krylov子空間迭代方法QMR對式(17)中的線性方程組進行求解，并基于不完全Cholesky分解得到預條件矩陣[33]。

1.3 散度校正技術

在迭代求解正演線性方程組時，盡管采用了基于不完全Cholesky分解的預條件矩陣以加速迭代求解，但對于頻率較低或非結構化網格質量較差時，迭代過程無法經較少迭代次數收斂到預期殘差，甚至可能不收斂。常用的改善措施是在迭代求解過程中對電流密度的散度進行校正。散度校正最早由Smith[26]在三維大地電磁有限差分正演中提出，爾后普遍應用于結構化網格三維電磁正演迭代求解[26-29]，但在非結構化網格正演中的應用尚未見報道。以下將以非結構化網格正演為例探討散度校正方法。

對式(8)中的第1個方程兩邊求散度并略去常系數項，可得

·(σEs+ΔσEp)=0

(19)

上式即為二次電場需滿足的散度條件。在剖分單元內，由于矢量元插值函數的散度為零，式(19)自動滿足。但由于矢量元插值函數定義在棱邊上，其方向是棱邊切向方向，故僅可保證切向分量的連續性，無法保證法向分量的躍變數值，法向躍變數值需要依靠電導率項約束。因此，對于低頻正演，由于雙旋度方程中的電導率項比重很小，故迭代初期在不同剖分單元的界面上有可能不滿足散度條件。當在界面上的散度條件不滿足時，通過這些過電荷計算對應的靜態電勢，再對迭代得到的近似電場進行校正，則校正之后的電場滿足散度條件。在迭代求解有限元方程時，交替施加散度校正，與未進行散度校正的場相比，近似電場會更快地收斂到真實解，從而達到加快迭代求解的目的。

通過求解偏微分方程

·(σφ)=0 ∈Ω

(20)

得到靜態電勢φ。圖3為電導率分界面上的散度校正示意圖，結合該分界面上的內邊界條件

圖3 電導率分界面上的散度校正示意圖

n·(σiφi-σjφj)|Sij=

n·(σiEsi-σjEsj+ΔσiEpi-ΔσjEpj) ∈Γ

(21)

式中：Sij為單元i與單元j的界面；φj表示單元j的靜態電勢。同時，界面兩邊的一次場及二次場必須同時取在棱邊上，保持一致性。取外邊界條件為勢函數在邊界上為零。式(20)可通過傳統的節點有限單元法求解。設剖分區域內有M個節點，則電勢可表示為

(22)

式中：vk為節點k的插值函數；φk為節點k的電勢。一次場及二次場可通過各條棱邊上的電場表示為

(23)

(24)

式中：ωj為矢量插值函數；N為剖分區域內總的棱邊數。求解偏微分方程(式(20)、式(21))可等效于求解線性方程組

Dp=t

(25)

式中：

將復線性方程組(式(24))的實部及虛部分開，可得到兩組實線性方程組，可分別應用預處理雙共軛梯度法求解。然后，通過下式進行散度校正

(26)

以上所述的散度校正過程可以在迭代法求解過程中交替進行，這樣可加快迭代法的收斂速度，以更短的計算時間求得滿足所需精度的解。

2 數值模擬實驗

為了驗證所提算法的有效性，設計了三個地電模型。首先通過水平層狀介質模型驗證散度校正的有效性及其對正演精度的影響。在此基礎上，對一個三維低阻體模型進行模擬，并對比基于二次耦合勢的正演結果，進一步驗證所提算法的有效性。最后，通過一個油氣模型證實該算法在油氣勘探數值模擬中的適用性。

2.1 水平層狀介質模型

圖4為一個水平層狀介質模型的yOz面示意圖，該模型為三層介質模型。采用右手坐標系統，z軸垂直向下，x、y方向在水平面內正交。將沿x方向的單位電偶極置于原點激發諧變電磁場信號。

圖4 水平三層介質模型yOz面示意圖

為了驗證本文算法的有效性，采用同樣的迭代求解方法對未施加散度校正與施加了散度校正后的迭代收斂情況及數值解進行對比。利用開源的TetGen工具[34]對該模型進行剖分,采用統一的剖分網格，分別計算1、4、8Hz下的電磁場。圖5為正演線性方程組的迭代收斂曲線圖。其中圖5b為迭代過程中每迭代100次施加一次散度校正得到的結果。由圖可知，當不施加散度校正時，在3000次迭代后各個頻率下的正演線性方程組無法收斂到預期殘差(1×10-6)；而施加散度校正后，盡管在初期迭代殘差有所上升，但之后殘差迅速下降，總體上使得各個頻率下的正演線性方程組迭代過程明顯得到加速，且成功收斂到預期殘差。圖中還給出了各個頻率下的正演計算時間t，正演計算時間基本與迭代次數成線性關系，說明散度校正過程幾乎不會額外增加正演計算時間。總體而言，施加了散度校正后，由于迭代收斂加快，正演計算時間也有所減少，且計算頻率越高，加速越快。圖6展示了頻率為1Hz時本文方法計算得到的二次電場數值解與解析解的對比。由圖可知，未施加散度校正時正演得到的數值解與解析解差別很大，無法滿足精度要求，而施加散度校正之后得到的數值解與解析解吻合較好，精度明顯提高。該數值實驗表明散度校正可加速非結構化有限元正演，且相比未施加散度校正時數值解精度更高。

圖5 水平層狀介質模型線性方程組迭代求解收斂曲線對比

圖6 水平層狀介質模型二次電場分量數值解與解析解對比(1Hz)

2.2 三維低阻體模型

為了進一步驗證本文算法對三維地電模型模擬的有效性，設計了一個低阻體埋藏于均勻半空間的三維地電模型(圖7)。模型采用右手坐標系，Idl為中心點位于坐標原點、沿x方向的單位電偶源。低阻體體積為800m(x)×800m(y)×200m(z)，關于x軸對稱。由于三維模型不存在解析解，以下將本文算法的數值解與經過精度驗證的基于二次耦合勢正演計算結果[35]對比。

圖7 三維地電模型示意圖

采用四面體網格對該模型進行離散，分別計算8、32、128Hz時的電磁場。圖8所示為該模型的部分網格剖分切面圖。藍色區域代表良導體所在區域，采用較密的網格，四面體網格棱邊長不大于40m。白色線條代表位于地面長為2km的沿x方向的測線，點距為100m，坐標范圍為-1000～1000m。對測線附近區域采用最密的網格，網格棱邊長不大于20m。從觀測區域往外網格逐步加大。圖9所示為各頻率下正演線性方程組迭代收斂曲線，迭代過程中每100次迭代施加1次散度校正。由圖可知，未施加散度校正時線性方程組的相對殘差降至1×10-4后無法繼續下降；而施加散度校正后各個頻率下都能在3000次內收斂到預期殘差，且總體正演計算時間有所減少，尤其是頻率為128Hz時，正演加速比達1.77。

圖8 三維地電模型網格剖分切面圖

圖9 三維地電模型不同頻率時的迭代收斂曲線

圖10是計算得到的8Hz時x方向的二次電場Exs在地面位置的平面圖。由圖可知，本文算法(基于二次電場)與基于二次耦合勢得到的電場分布形態保持一致，具有很高的吻合度。

圖10 基于二次耦合勢(a)及二次電場(b)計算的8Hz時Exs數值解實部(上)及虛部(下)

進一步地，對基于二次電場和基于二次耦合勢計算得到的Exs和Eys進行對比，以證實本文方法的可靠性。

圖11和圖12分別為8Hz和32Hz時y=3km測線上二次電場分量Exs對比，其中圖11給出了相對誤差曲線。由圖可見，總體上基于二次電場與基于二次耦合勢計算得到的Exs數值解吻合較好。另外，從圖11可知，除了在x=-500m處相對誤差偏大(達17%)之外，實部分量相對誤差都小于5.0%，虛部分量相對誤差都小于1.4%。進一步分析圖11a可知，在x=-500m附近實部分量剛好在零值附近變換符號，場值非常小，因此其絕對誤差足夠小，但相對誤差偏大，這是符合預期的。

圖11 8Hz時y=3km剖面分別基于二次電場和二次耦合勢計算得到的Exs數值解對比

圖12 32Hz時y=3km剖面分別基于二次電場和二次耦合勢計算得到的Exs數值解對比

圖13為y=5km測線上基于二次電場和基于二次耦合勢計算得到的二次電場Exs電場分布。由圖可知，電場實部和虛部分量在2km長的測線范圍內未改變符號，此時基于二次電場與基于二次耦合勢的數值解實部和虛部吻合很好，相對誤差都在5%以下。

圖13 16Hz時y=5km剖面分別基于二次電場和二次耦合勢計算得到的Exs數值解對比

圖14和圖15分別為8Hz、32Hz時二次電場Eys數值解對比。總體而言，除了在峰值位置差異稍大外，基于二次電場和基于二次耦合勢計算得到的Eys數值解吻合較好。因此，該三維地電模型的數值解對比進一步說明了本文算法的可靠性。

圖14 8Hz時y=3km剖面分別基于二次電場和二次耦合勢計算得到的Eys數值解對比

圖15 32Hz時y=3km剖面分別基于二次電場和二次耦合勢計算得到的Eys數值解對比

2.3 三維油氣監測模型

為了分析可控源電磁法在油氣監測方面應用的可能性，設計了一個水驅油氣監測模型(圖16)。單位電偶源置于坐標原點。該水平三層介質模型包含一個2000m×2000m×150m規模的長方體異常區域，一邊是低阻的水，另一邊是高阻油氣，異常區域關于x軸對稱。該模型可模擬從左邊注入水驅動油氣的過程。通過設置不同的水/油氣的體積占比，分析油氣占比分別為100%、50%及0時電場分量Ex和Ey分布，計算頻率為1Hz，結果如圖17所示。由圖可知，電場分量是關于x軸對稱的，這與模型本身關于x軸對稱是一致的，同時，在遠離電偶源區域，場的變化較為明顯。為了更清晰展示電場振幅隨偏移距的變化特征，繪制了圖18。由圖可知，油氣含量從100%降至50%時，電場振幅變化較大，基本在偏移距y>1km時異常明顯。當油氣含量從50%降至0時，電場振幅的變化相對弱一些，直至偏移距y>3km時異常才變得明顯。這說明當水驅動油氣發生變化時，地面的電場分量變化是可探測的，有規律的。該模型實例說明本文算法可用于可控源油氣監測模型的模擬。

圖16 三維水驅油氣監測模型yOz斷面圖

圖17 三維油氣監測模型電場水平分量Ex(上)和Ey(下)振幅

圖18 三維油氣監測模型x=-1.75km(上)和x=-1.00km(下)剖面電場分量Ex(a)和Ey(b)振幅曲線

3 結論

本文開展了基于二次電場的三維可控源電磁有限元正演研究，通過結合非結構化網格和散度校正技術，為基于非結構網格的三維電磁正演線性方程組的迭代求解提供了技術支撐。層狀介質模型和兩個三維地電模型的數值模擬實驗表明，散度校正技術適用于非結構化網格的三維正演，可加速復雜三維模型電磁響應的正演計算。通過在迭代過程中交替施加散度校正，可加快迭代求解收斂速度，在較少的迭代次數收斂到預期殘差，從而達到減少正演計算時間和節約內存需求的目的。相比未施加散度校正的計算方法，本文算法能進一步提高正演精度。