湯文武 鄧居智 黃清華

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

0 引言

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

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

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

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

1 方法原理

1.1 三維電磁正演邊值問題

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

(1)

(2)

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

××E=iωμ0×H

(3)

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

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

(4)

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

E=Ep+Es

(5)

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

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

(6)

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

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

(7)

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

(8)

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

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

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

圖1 研究區(qū)域示意圖

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

設(shè)對應(yīng)各條棱邊的形函數(shù)為Ni，那么一次電場及二次電場可分別表示為

(9)

(10)

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

Ni=(Li1Li2-Li2Li1)li

(11)

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

(12)

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

(13)

根據(jù)矢量恒等式

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

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

×Es·(×Ni)]dΩ

s×Es×Ni·ds

(14)

式中S表示面積矢量。根據(jù)二次電場和二次磁場的關(guān)系，當邊界離異常區(qū)域足夠遠時，滿足

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

(15)

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

(16)

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

KEs=MEp

(17)

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

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

1.3 散度校正技術(shù)

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

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

·(σEs+ΔσEp)=0

(19)

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

通過求解偏微分方程

·(σφ)=0 ∈Ω

(20)

得到靜態(tài)電勢φ。圖3為電導(dǎo)率分界面上的散度校正示意圖，結(jié)合該分界面上的內(nèi)邊界條件

圖3 電導(dǎo)率分界面上的散度校正示意圖

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

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

(21)

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

(22)

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

(23)

(24)

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

Dp=t

(25)

式中：

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

(26)

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

2 數(shù)值模擬實驗

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

2.1 水平層狀介質(zhì)模型

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

圖4 水平三層介質(zhì)模型yOz面示意圖

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

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

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

2.2 三維低阻體模型

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

圖7 三維地電模型示意圖

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

圖8 三維地電模型網(wǎng)格剖分切面圖

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

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

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

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

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

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

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

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

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

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

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

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

2.3 三維油氣監(jiān)測模型

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

圖16 三維水驅(qū)油氣監(jiān)測模型yOz斷面圖

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

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

3 結(jié)論

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