基于矢量位和標量位的廣域電磁法三維有限元數(shù)值模擬

周印明 王金海 胡曉穎 何展翔 熊 彬

(①中南大學地球科學與信息物理學院，湖南長沙 410083； ②中南大學有色金屬成礦預測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室，湖南長沙 410083； ③東方地球物理公司綜合物化探處，河北涿州 072751； ④青海省第三地質(zhì)勘查院，青海西寧 810029； ⑤南方科技大學前沿與交叉科學研究院,廣東深圳 518055； ⑥桂林理工大學地球科學學院，廣西桂林 541006)

0 引言

Goldstein[1]提出以人工場源代替天然場源的可控源音頻大地電磁法(CSAMT)，有效克服了大地電磁法(MT)和音頻大地電磁法(AMT)天然場源的隨機性和信號強度弱等缺陷。該方法采用接地導線和不接地回線為場源，在其一側(cè)或者兩側(cè)60°張角的扇形區(qū)域內(nèi)測量相互正交的電磁場切向分量，沿用MT中計算卡尼亞視電阻率的公式，且保留了AMT的數(shù)據(jù)解釋方法。與MT相比，CSAMT采用人工場源，可彌補場源的隨機性和信號微弱等不足。同時，該方法在工作效率及縱、橫向分辨率等方面均有明顯改進。因此，該方法經(jīng)過40多年的長足發(fā)展，在金屬礦產(chǎn)、水文、工程、環(huán)境等領域均得到廣泛應用。但是，CSAMT也存在不足：①沿用MT法采用的卡尼亞視電阻率公式，只適用于“遠區(qū)”數(shù)據(jù)，對于“非遠區(qū)”數(shù)據(jù)，該公式便不再適用，這明顯背離了采用人工源增強信號強度、提高觀測精度的初衷，并且實際測量數(shù)據(jù)中很容易進入“過渡區(qū)”或“近區(qū)”數(shù)據(jù)；②CSAMT依然遵從卡尼亞公式計算視電阻率，該公式相對比較簡單，人為地舍棄了代表“非遠區(qū)”數(shù)據(jù)特點的高次項，引入了不必要的人為誤差；③由于建場的局限性，探測深度受發(fā)射功率的限制。

對于CSAMT，何繼善[2]利用人工場源克服了場源的隨機性，利用磁偶源頻率測深法(MELOS)“非遠區(qū)”觀測優(yōu)勢，提出了廣域電磁法。該方法摒棄了CSAMT方法“遠區(qū)”信號微弱的觀測方式，拓展了觀測范圍；也摒棄了MELOS校正計算視電阻率的方法，保留了計算公式中的高次項。在理論上提出一種適用于全域的視電阻率計算公式，從根本上突破了CSAMT法“遠區(qū)”理論的束縛，有效擴展了人工源電磁法的觀測范圍，提高了觀測精度和野外工作效率[3]。

目前廣域電磁法的資料處理與解釋技術以一維和二維為主。但是，隨著勘探難度的日益增加和方法技術應用的不斷深入，人們對廣域電磁法的勘探效果提出了更高的要求，因而以大范圍觀測為代表的三維精細勘探將逐漸發(fā)展成為廣域電磁法的研究趨勢和熱點。其中，廣域電磁法三維數(shù)值模擬作為反演成像與定量解釋的基礎，尤為重要。有關頻率域電磁法三維正演方法的研究已經(jīng)超過了30年[4-5]，并且在各個方面均取得了很大的進展。從數(shù)值方法的角度而言，常用的電磁三維數(shù)值模擬方法主要有積分方程法[6-10]、有限差分法[11-13]、有限體積法[14-19]、有限單元法[20-24]等。盡管這些方法采用的數(shù)值算法有所不同，但總體研究目標都是提高算法的計算精度與計算效率。目前專門針對廣域電磁法三維數(shù)值模擬研究主要包括：①李帝銓等[25]采用積分方程法實現(xiàn)了廣域電磁法三維數(shù)值模擬，并分析了典型三維地電模型的電磁響應；②彭榮華等[19]采用交錯網(wǎng)格有限體積法實現(xiàn)了基于二次耦合勢的廣域電磁法三維正演。

本文從麥克斯韋方程組出發(fā)，推導了基于庫侖規(guī)范的矢量位和標量位控制方程。該方程具有非常完備的物理意義，可以將大地電磁法(MT)、直流電阻率法(DC)和可控源電磁法(CSEM)統(tǒng)一到一個耦合方程組。通過改變場源和頻率，基于庫侖規(guī)范的矢量位和標量位控制方程可以適用于MT、DC和CSEM三維數(shù)值模擬。廣域電磁法屬于頻率域可控源電磁法(CSEM)中的一種方法，其基本理論完全相同?；诖?，本文以庫侖規(guī)范的矢量位和標量位控制方程為基礎，對廣域電磁法三維有限元數(shù)值模擬進行了研究。在模型算例中，設計棱柱體模型，對比了傳統(tǒng)積分方程算法與本文算法的計算結(jié)果，驗證了本文方法的正確性和計算精度。在此基礎上，對比分析了廣域電磁法與CSAMT對典型三維目標體的探測能力，結(jié)果表明在相同的條件下，廣域電磁法能夠更準確地反映地下目標體信息，具有更高的分辨能力。

1 理論與方法

1.1 基于庫侖規(guī)范的矢量位和標量位控制方程

假定時間諧變因子為e-iωt，頻率域麥克斯韋方程組可以表示為

×E=iωμ0H

(1)

×H=Jp+(σ-iωε0)E

(2)

式中：E為電場強度；H為磁場強度；ω是角頻率；μ0是真空磁導率；σ為電導率；Jp為場源的電流密度；ε0為真空中的介電常數(shù)，對于廣域電磁法頻段而言，ε0可以忽略。

對式(1)兩邊取旋度，并將其帶入式(2)可得關于電場強度E的二階亥姆霍茲方程

××E-iωμ0σE=iωμ0σJp

(3)

采用數(shù)值方法，直接求解式(3)可得電場強度，既而得到磁場的分布。直接求解電磁場存在兩個明顯問題[26]：一是由于只利用了場矢量的旋度，對散度沒有規(guī)范，存在“偽解現(xiàn)象”；二是電磁場分量在界面上不連續(xù)，與節(jié)點有限元在求解區(qū)域內(nèi)連續(xù)這一基本假設相矛盾，這給求解帶來很大的困難。一些學者利用矢量有限單元法[27-29]有效解決了以上問題，但同時也增加了方程組的復雜度。為了避免上述問題，先求取電磁場的勢函數(shù)，進而通過差分格式間接求解電磁場[20,30]。另外，這種間接求解方法還可以有效減小方程組系數(shù)矩陣的條件數(shù)，顯著提高收斂速度[11]。

將磁矢量位和電標量位與電磁場的基本關系式[20]

B=×A

(4)

E=iωA-Φ

(5)

代入式(2)，可得磁矢量位和電標量位的雙旋度方程

×(×A)=μ0Jp+μ0σ(iωA-Φ)

(6)

式中：A為矢量位；B為磁感應強度；Φ為標量位。

把矢量恒等式

×(×A)=(·A)-2A

代入式(6)，得到

(·A)-2A=μ0Jp+μ0σ(iωA-Φ)

(7)

為了保證勢函數(shù)的唯一性，把庫侖規(guī)范

·A=0

(8)

代入式(7)，得到

2A+iωμσA-μ0σΦ=-μ0Jp

(9)

式中μ為相對磁導率。

·(×H)=·J=0

(10)

-·Js=

(11)

將式(5)代入式(11)可得

Φ)-iω·Jp

(12)

·(fA)=f·A+A·f

并將式(8)代入式(12)，可得

σ·Φ+Φ·σ-iωA·σ=·Jp

(13)

式中f表示任意標量位。

綜合式(9)和式(13)可得基于庫侖規(guī)范的矢量位和標量位滿足的耦合偏微分方程

(14)

將式(14)展開，寫成分量形式

(15)

由式(15)可以看出：

(1)該控制方程由4個方程表達式組成，矢量位A的三個分量形式Ax、Ay和Az之間本身并無直接聯(lián)系，而是通過標量位Φ的一階偏導耦合在一起的。

(2)該控制方程相對比較完備，物理意義明確，可以通過改變場源和頻率用于多種電磁勘探方法數(shù)值模擬：

1.2 邊界條件

目前可控源電磁法三維數(shù)值模擬中常用的邊界條件主要包括Dirichlet邊界條件(令邊界處切向電場為零)[20]和Neumann邊界條件(令邊界處切向電場的空間一階導數(shù)為零)[31-32]。這兩種邊界條件都是通過擴邊來降低邊界反射的影響。Streich[13]對上述兩種邊界條件進行過對比，認為Neumann邊界條件會使線性方程組系統(tǒng)不對稱。如果以截斷邊界上的近似場值作為邊界條件，那么截斷邊界造成的邊界反射就能得到相應的壓制。在可控源電磁法中，總場為一次場和二次場的疊加，目標體產(chǎn)生的二次場由一次場激發(fā)，其數(shù)值遠遠小于一次場。當截斷邊界離計算區(qū)域較遠時，可以忽略二次場對邊界的影響，進而利用邊界上一次場近似代替總場作為邊界條件。

在均勻介質(zhì)中，矢量位Ax、Ay和Az滿足微分方程

(16)

以式(16)中第一式為例，其基本解為

Ax=Ae-kz+Bekz

(17)

由于上界面只有上行波，從而可以得到上邊界條件

(18)

同理，可以得到下邊界條件、左邊界條件、右邊界條件、前邊界條件和后邊界條件分別為

(19)

(20)

(21)

(22)

(23)

式中下標“min”和“max”分別表示極小值和極大值。

式(18)～式(23)構(gòu)成了矢量位在截斷邊界的邊界條件?？梢钥闯觯噶课挥砷L導線源產(chǎn)生的電磁感應場構(gòu)成，在邊界上可以當成平面波處理，與場源的大小和位置無關。

在截斷邊界上，標量位Φ的微分方程直接簡化為直流電阻率法三維控制方程，故可以將直流電阻率法中的混合邊界條件[33-34]作為標量位Φ在截斷邊界上的邊界條件，具體形式如下

(24)

(25)

2 有限單元法

基于庫侖規(guī)范的矢量位和標量位控制方程(式(15))及邊界條件(式(18)～式(25))構(gòu)成了廣域電磁法三維數(shù)值模擬的邊值問題。針對該邊值問題，采用加權(quán)余量法將矢量位和標量位滿足的邊值問題轉(zhuǎn)換為變分問題，并采用格林公式進行降階處理，可得

(26)

式中：Ni表示第i個節(jié)點的線性插值函數(shù)；nx、ny、nz表示邊界上的方向法向分量；v和s分別表示體積分單元和面積分單元。

本文采用圖1所示的離散方案進行網(wǎng)格剖分。首先將三維模型進行六面體單元剖分，然后在六面體中進行四面體單元的剖分[35-36]，每個六面體單元均被剖分成六個任意形狀的四面體單元。該剖分方法具有很好的數(shù)值精度對稱性[37]。在每個四面體單元內(nèi)，采用線性插值，其計算表達式參見文獻[34]。對變分問題式逐項進行單元分析，單元積分過程用到的積分表參見文獻[34]。得到擴展矩陣或陣列

圖1 六面體單元及其四面體單元剖分示意圖

As=b

(27)

式中：A是大型稀疏復矩陣；s是解向量：b是源向量。

目前，求解大型稀疏線性方程組的方法主要有迭代法[38]和直接法[39]。

迭代解法的優(yōu)點是算法相對簡單，編程容易實現(xiàn)，耗費計算資源較少，當稀疏矩陣的條件數(shù)不大(良態(tài))時收斂較快。缺點是當稀疏矩陣的條件數(shù)過大(病態(tài))時收斂性很差或者發(fā)散，且在有限迭代次數(shù)內(nèi)求解精度無法保證。采用最多的迭代解法是復線性方程組穩(wěn)定雙共軛梯度法[26,31,37,40]。由于地—空和地下介質(zhì)內(nèi)部電阻率變化較大，以及網(wǎng)格的非均勻剖分，導致系數(shù)矩陣出現(xiàn)嚴重的病態(tài)或者條件數(shù)很大，從而使得迭代解法收斂很慢甚至發(fā)散。因此需要采用預處理技術對系數(shù)矩陣進行修正，降低其條件數(shù)。應用較多的預處理方法包括對角線系數(shù)、不完全Cholesky分解、SSOR及不完全LU分解法。預處理算法本身并不復雜，但是由于CSEM三維有限元數(shù)值模擬中系數(shù)矩陣是稀疏的，因此在存儲的過程中通常只存入非零元素，通過編碼的方式記錄該非零元素對應的行和列位置。

直接解法的優(yōu)點是求解精度較高，多源情形下可實現(xiàn)快速回代求解，當稀疏矩陣的條件數(shù)很大(病態(tài)或嚴重病態(tài))時也可得到比較穩(wěn)定的解。缺點是當稀疏矩陣維數(shù)很大時，需要占用大量的內(nèi)存。近年來隨著計算機技術的飛速發(fā)展和大型稀疏矩陣分解算法的不斷優(yōu)化[41-42]，利用直接解法求解電磁法三維數(shù)值模擬中大型復線性稀疏方程組的案例逐漸增多[18,24,36,43-44]。本文調(diào)用Pardiso_64位并行求解器進行求解，得到矢量位和標量位計算結(jié)果。

3 電磁場分量的計算

根據(jù)電磁場與矢量位和標量位之間的關系式(式(4)和式(5))可以得到電磁場分量表達式

(28)

式中的求導運算可采用五點差分法。

4 廣域視電阻率的計算

在可控源音頻大地電磁法中，主要通過電場和磁場的正交分量計算卡尼亞視電阻率[1]

(29)

式中：Ex和Ey為地面電場水平分量；Hx和Hy為地面磁場水平分量。

卡尼亞視電阻率具有計算方便、簡單的特點，式(29)在“遠區(qū)”能比較客觀地反應地電斷面的變化，但在“過渡區(qū)”和“近區(qū)”計算的卡尼亞視電阻率會發(fā)生嚴重畸變。

為了突破卡尼亞視電阻率在“過渡區(qū)”和“近區(qū)”的局限性，何繼善[2]根據(jù)電磁場表達式的特點，定義了廣域視電阻率，將電磁測深的范圍擴大到包括“遠區(qū)”在內(nèi)的廣大區(qū)域。電偶極源激勵的電場分量Ex對應的廣域視電阻率為

(30)

5 數(shù)值試驗

設計一個均勻半空間模型，驗證本文算法的正確性及可靠性。在此基礎上，設計長方體模型，利用本文正演算法對比分析廣域電磁法與CSAMT對典型三維目標體的探測能力。算法代碼采用Fortran 95語言編寫，計算平臺為Intel(R) Xeon(R) CPU3.1G，128GB RAM，16CPUs。

5.1 均勻半空間模型

設定均勻半空間模型參數(shù)：電阻率為10000Ω·m，有限長線源沿x方向布設，起點為坐標原點，長度為500m，發(fā)射電流為1A，發(fā)射頻率為1Hz。模擬計算區(qū)域大小為10km(x)×10km(y)×6km(z)，x、y、z方向網(wǎng)格距均為100m，網(wǎng)格節(jié)點數(shù)為101×101×61，水平方向各取5個節(jié)點作擴邊處理。圖2為解析解[46]、本文算法計算的數(shù)值解及相對誤差。

從圖2a～圖2c可以看出，電場分量Ex的數(shù)值解與解析解吻合度高，y=0時的最大相對誤差為1.4%。從圖2d～圖2f可以看出，磁場分量Hy的數(shù)值解與解析解吻合度高，y=0時的最大相對誤差為1.2%，這是因為場源具有奇異性，場源附近的相對誤差相對較大，其他區(qū)域的數(shù)值精度較高，但都在誤差允許的范圍內(nèi)，從而驗證了本文算法的正確性及高精度。

5.2 典型三維目標體模型

設計圖3所示的三維低阻模型[19]，分別進行可控源音頻大地電磁法(CSAMT)和廣域電磁法(WFEM)模擬計算，分析二者的分辨能力。

圖3 三維低阻模型示意圖

以x方向的電偶極子作為發(fā)射源(Tx)，發(fā)射頻點數(shù)為12，頻率在0.1～500Hz范圍內(nèi)呈對數(shù)等間隔分布。目標體剖分單元網(wǎng)格尺寸為100m×100m×50m，網(wǎng)格數(shù)為51×81×51，其中空氣層深度方向網(wǎng)格剖分為11層。

圖4為主測線y=0上卡尼亞視電阻率與廣域視電阻率擬斷面圖。圖5是頻率分別為1.02、2.21、4.80和10.40 Hz時卡尼亞視電阻率和廣域視電阻率響應切片?？梢钥闯觯孩倏醽喴曤娮杪?圖4a和圖5上)不能較好地反映低阻異常體的空間位置。當發(fā)射頻率較低時，測線上卡尼亞視電阻率已經(jīng)遠大于100Ω·m；隨著頻率的降低，卡尼亞視電阻率與真實電阻率偏離增大，這是由于頻率越低，趨膚深度越大，y=0主測線位置已經(jīng)不滿足“遠區(qū)”觀測條件，即超出卡尼亞視電阻率公式的適用范圍。②廣域視電阻率(圖4b和圖5下)能清晰地反映低阻異常體的空間位置，不受“遠區(qū)”測量條件的限制；頻率的增大或降低都不影響廣域視電阻率的異常幅值，即使在頻率很低時也能很好地反映低阻異常體。所以，廣域電磁法具有更大的有效觀測范圍，可以提高野外施工效率。

圖4 主測線 y=0上卡尼亞視電阻率(a)與廣域視電阻率(b)擬斷面圖

圖5 不同頻率時卡尼亞視電阻率(上)與廣域視電阻率(下)切片對比(a)1.02Hz; (b)2.21Hz; (c)4.80Hz;(d)10.40Hz

6 結(jié)論

本文從頻率域麥克斯韋方程出發(fā)，推導了庫侖條件下的矢量位和標量位耦合方程，并結(jié)合邊界條件給出了電磁場矢量位和標量位滿足的邊值問題。利用有限單元法對矢量位和標量位控制方程進行離散化，實現(xiàn)了廣域電磁法的三維數(shù)值模擬，通過與均勻半空間模型的地面解析解對比，驗證了本文三維數(shù)值模擬算法的正確性和高精度。在此基礎上，設計了一個三維低阻模型，對比分析了廣域電磁法與CSAMT對典型三維目標體的探測能力。

研究結(jié)果表明：廣域電磁法提出了一種適用于全域的視電阻率計算公式，從根本上突破了CSAMT“遠區(qū)”測量的束縛，在非“遠區(qū)”依然能夠有效反映地下三維目標體的電性特征和分布范圍，有效地擴展了人工源電磁法的觀測范圍，提高了觀測精度和野外工作效率。