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

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

0 引言

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

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

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

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

1 理論與方法

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

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

×E=iωμ0H

(1)

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

(2)

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

對(duì)式(1)兩邊取旋度，并將其帶入式(2)可得關(guān)于電場(chǎng)強(qiáng)度E的二階亥姆霍茲方程

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

(3)

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

將磁矢量位和電標(biāo)量位與電磁場(chǎng)的基本關(guān)系式[20]

B=×A

(4)

E=iωA-Φ

(5)

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

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

(6)

式中：A為矢量位；B為磁感應(yīng)強(qiáng)度；Φ為標(biāo)量位。

把矢量恒等式

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

代入式(6)，得到

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

(7)

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

·A=0

(8)

代入式(7)，得到

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

(9)

式中μ為相對(duì)磁導(dǎo)率。

·(×H)=·J=0

(10)

-·Js=

(11)

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

Φ)-iω·Jp

(12)

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

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

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

(13)

式中f表示任意標(biāo)量位。

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

(14)

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

(15)

由式(15)可以看出：

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

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

1.2 邊界條件

目前可控源電磁法三維數(shù)值模擬中常用的邊界條件主要包括Dirichlet邊界條件(令邊界處切向電場(chǎng)為零)[20]和Neumann邊界條件(令邊界處切向電場(chǎng)的空間一階導(dǎo)數(shù)為零)[31-32]。這兩種邊界條件都是通過(guò)擴(kuò)邊來(lái)降低邊界反射的影響。Streich[13]對(duì)上述兩種邊界條件進(jìn)行過(guò)對(duì)比，認(rèn)為Neumann邊界條件會(huì)使線性方程組系統(tǒng)不對(duì)稱。如果以截?cái)噙吔缟系慕茍?chǎng)值作為邊界條件，那么截?cái)噙吔缭斐傻倪吔绶瓷渚湍艿玫较鄳?yīng)的壓制。在可控源電磁法中，總場(chǎng)為一次場(chǎng)和二次場(chǎng)的疊加，目標(biāo)體產(chǎn)生的二次場(chǎng)由一次場(chǎng)激發(fā)，其數(shù)值遠(yuǎn)遠(yuǎn)小于一次場(chǎng)。當(dāng)截?cái)噙吔珉x計(jì)算區(qū)域較遠(yuǎn)時(shí)，可以忽略二次場(chǎng)對(duì)邊界的影響，進(jìn)而利用邊界上一次場(chǎng)近似代替總場(chǎng)作為邊界條件。

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

(16)

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

Ax=Ae-kz+Bekz

(17)

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

(18)

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

(19)

(20)

(21)

(22)

(23)

式中下標(biāo)“min”和“max”分別表示極小值和極大值。

式(18)～式(23)構(gòu)成了矢量位在截?cái)噙吔绲倪吔鐥l件?？梢钥闯?，矢量位由長(zhǎng)導(dǎo)線源產(chǎn)生的電磁感應(yīng)場(chǎng)構(gòu)成，在邊界上可以當(dāng)成平面波處理，與場(chǎng)源的大小和位置無(wú)關(guān)。

在截?cái)噙吔缟?，?biāo)量位Φ的微分方程直接簡(jiǎn)化為直流電阻率法三維控制方程，故可以將直流電阻率法中的混合邊界條件[33-34]作為標(biāo)量位Φ在截?cái)噙吔缟系倪吔鐥l件，具體形式如下

(24)

(25)

2 有限單元法

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

(26)

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

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

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

As=b

(27)

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

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

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

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

3 電磁場(chǎng)分量的計(jì)算

根據(jù)電磁場(chǎng)與矢量位和標(biāo)量位之間的關(guān)系式(式(4)和式(5))可以得到電磁場(chǎng)分量表達(dá)式

(28)

式中的求導(dǎo)運(yùn)算可采用五點(diǎn)差分法。

4 廣域視電阻率的計(jì)算

在可控源音頻大地電磁法中，主要通過(guò)電場(chǎng)和磁場(chǎng)的正交分量計(jì)算卡尼亞視電阻率[1]

(29)

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

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

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

(30)

5 數(shù)值試驗(yàn)

設(shè)計(jì)一個(gè)均勻半空間模型，驗(yàn)證本文算法的正確性及可靠性。在此基礎(chǔ)上，設(shè)計(jì)長(zhǎng)方體模型，利用本文正演算法對(duì)比分析廣域電磁法與CSAMT對(duì)典型三維目標(biāo)體的探測(cè)能力。算法代碼采用Fortran 95語(yǔ)言編寫，計(jì)算平臺(tái)為Intel(R) Xeon(R) CPU3.1G，128GB RAM，16CPUs。

5.1 均勻半空間模型

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

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

5.2 典型三維目標(biāo)體模型

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

圖3 三維低阻模型示意圖

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

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

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

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

6 結(jié)論

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

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