惠 鑫， 吳小平*

(1.中國科學技術大學 地球和空間科學學院 地震與地球內(nèi)部物理實驗室，合肥 230026 2. 蒙城地球物理國家野外科學觀測研究站，蒙城 233500)

0 引言

介質(zhì)的電性各向異性現(xiàn)象普遍存在[1]，在地電模型和地球動力學解釋中發(fā)揮著重要作用。早在20世紀70年代，文獻[2]和文獻[3]就已經(jīng)開始了對各向異性大地電磁(MT)的研究。在各向異性MT數(shù)值模擬領域， Pek等[4]公布了二維各向異性MT有限差分公式；Li[5]實現(xiàn)了二維各向異性介質(zhì)的MT有限元數(shù)值模擬，其采用線性插值基函數(shù)，數(shù)值精度在電磁場變化劇烈的地區(qū)受限，且其數(shù)值模擬結(jié)果沒有得到各向異性模型解析解的驗證；QIN L[6]給出了二維各向異性斷層模型的MT解析解，為各向異性二維MT的數(shù)值模擬提供了比對結(jié)果。依托MT有限差分正演模擬，胡祥云和霍光譜等[7-8]在實際大地電磁數(shù)據(jù)解釋中引入了各向異性。

在MT的應用方面，20世紀80年代，Cox[9]和Ranganayaki等[10]應用MT展開對海岸效應的研究，并給出了該模型同性介質(zhì)TM模式觀測的特征頻率和特征距離的經(jīng)驗公式；Constable等[11]在的TE模式觀測中也發(fā)現(xiàn)了海岸效應特征異常。

但是這些研究均基于各向同性介質(zhì)模型，與各向異性普遍存在的地球模型有一定出入，因此有必要研究各向異性介質(zhì)下的海岸效應現(xiàn)象。筆者從麥克斯韋方程組出發(fā)，推導了任意各向異性介質(zhì)二維大地電磁的非結(jié)構(gòu)二階有限元公式，實現(xiàn)了二維大地電磁有限元正演模擬。相比較一階插值，二階插值有限元模擬結(jié)果更精確，而非結(jié)構(gòu)有限元的優(yōu)勢在于，可以模擬任意復雜地形的MT觀測結(jié)果。

1 有限元理論分析

任意各向異性介質(zhì)下的麥克斯韋方程組如式(1)所示。

(1)

(2)

(3)

(4)

其中：

C=σ11+σ12B+σ13A

D=σ22σ33-σ23σ32

求解公式(2)和公式(3)，可以得到水平電磁場量Ex和Hx。

1.1 插值函數(shù)

圖1是二維非結(jié)構(gòu)網(wǎng)格剖分圖，采用Gmsh軟件進行剖分。圖2是非結(jié)構(gòu)網(wǎng)格中的某個三角形單元，編號為該單元內(nèi)插值節(jié)點的內(nèi)部編號。圖中的坐標系y軸為走向方向，x軸垂直于走向方向，z軸正方向垂直向下。我們在三角形單元內(nèi)進行二次差值，即假定電場Ex和磁場Hx是y和z的二次函數(shù)，可以近似為：

i=1,…,6

(6)

其中:Ei和Hi分別是全局坐標系yoz中三角形單元的第i個節(jié)點的電場和磁場;二次差值的插值函數(shù)Ni采用文獻 [12]中的定義，Ni公式為式(7)。

Ni(y,z)=2(Li-1)Lii=1,2,3

N4(y,z)=4L1L2

N5(y,z)=4L2L3

N6(y,z)=4L3L1

(7)

Li公式為式(8)。

c1=y2z3-z2y3;a1=z2-z3;b1=y3-y2

c2=y3z1-z3y1;a2=z3-z1;b2=y1-y3

(8)

1.2 邊界條件

筆者利用有限元方法求解的是一個關于電磁場的邊值問題，內(nèi)部邊界的電場切向分量Et和磁場切向分量Ht連續(xù)。外部邊界條件選用Dirichlet邊界條件。

左右兩側(cè)邊界加載時，將模型的邊界設置離異常體足夠遠，大概異常體規(guī)模10倍左右的距離，采用一維程序計算邊界場值并加載。上邊界(下邊界)場值的加載，則將左右兩側(cè)的上邊界(下邊界)值讀取出來，以y軸的距離為插值變量進行線性插值，從而計算出上邊界(下邊界)的場值。

1.3 單元分析

將公式(2)乘以電場的任意變分δEx，公式(3)乘以電場的任意變分δHx，并在模擬區(qū)域Ω上求積分。求解區(qū)域Ω由ne個三角形單元組成，單元編號記為e=1、2、…、ne，電阻率參數(shù)賦值到每一個三角形單元。利用矢量計算公式化簡可得到如下積分公式:

(9)

(10)

公式(9)和公式(10)中的Γe為單元e的邊界，ey和ez是方向單位向量。在內(nèi)邊界Γe上，電場的切向分量Et和磁場的切向分量Ht連續(xù)。在求積分的過程中，沿著每一個邊界都進行了兩次積分，積分的方向相反，因而，所有的內(nèi)部邊界積分之和為零。又因為在外邊界上我們采用Dirichlet邊界條件，故電場的變分δEx和磁場的變分δHx為零，故所有線性積分之和等于零，公式(9)和公式(10)可以進一步化簡為式(11)與式(12)。

(11)

▽δHxdΩ=0

(12)

圖1 二維非結(jié)構(gòu)單元Fig.1 2D unstructured element

圖2 三角形單元Fig.2 Triangle element

1.4 剛度矩陣

將插值基函數(shù)公式(6)代入式(11)和式(12)，并將單元矢量Exe、Hxe、δExe和δHxe移到積分外面(因為它們不再依賴于y和z)，可得到式(13)和式(14)。

(13)

(14)

KU=0

(15)

其中，

本文中調(diào)用Fortran中的pardiso庫函數(shù)，運用直接法求解系數(shù)矩陣。

圖3 網(wǎng)格剖分實例Fig.3 Grid example

圖4 矩陣稀疏Fig.4 Sparse matrix

1.5 視電阻率的計算

相對于各向同性介質(zhì)而言，各向異性介質(zhì)中的電場和磁場一般不能解耦，則阻抗張量Zxx和Zyy不為零，所以需要在兩種不同的極化方式下求解這4個阻抗，其公式如下：

(16)

E和H下標中的1和2代表兩種不同的極化方式。公式(16)中的輔助場Ey和Hy依據(jù)公式求得。根據(jù)公式(16)求得的阻抗值，就可以進一步求解視電阻率和相位，二者的計算公式為式(17)。

(17)

2 算法驗證

為了驗證本文任意各向異性介質(zhì)二維大地電磁非結(jié)構(gòu)有限元正演模擬算法(MT2DFEANI)的正確性，模擬了一維和二維測試模型的大地電磁響應，并依據(jù)公式(18)，計算模擬數(shù)據(jù)Dc與對比數(shù)據(jù)Dr的相對誤差r。

(18)

2.1 一維模型驗證

一維層狀各向異性模型的大地電磁響應，存在解析解，因此設了一個三層的電阻率參數(shù)為“高低高”的模型和一個三層的電阻率參數(shù)為“低高低”的模型。兩個模型的第一層和第二層深度均為3 km，其中“高低高”模型(D模型)的第一層和第三層的電阻率為1 000 Ω·m，第二層為電阻率各向異性層，其主軸電阻率一次為參數(shù)30 Ω·m、100 Ω·m、30 Ω·m，水平旋轉(zhuǎn)角度αs為30°。低高低模型(H模型)的第一層和第三層的電阻率為100 Ω·m，第二層為電阻率各向異性層，其主軸電阻率一次為參數(shù)300 Ω·m、1 000 Ω·m、300 Ω·m，水平旋轉(zhuǎn)角度αs為30°。

對比圖5和圖 6，可以得到三點結(jié)論：①模擬結(jié)果對比中，高頻段的模擬結(jié)果出現(xiàn)振蕩現(xiàn)象，而在低頻段模擬結(jié)果很好，最大相對誤差在2%以內(nèi)；②在周期T=10 s之后，本文算法的模擬結(jié)果相當準確，相對誤差接近零；③相較一次插值有限元模擬結(jié)果，二次插值有限元模擬的結(jié)果能夠在低頻段出現(xiàn)較小的震蕩現(xiàn)象，且相對較快的收斂到解析解結(jié)果，證明了本文二次插值程序的優(yōu)勢性。

2.2 二維模型驗證

圖7是d'Erceville等[13]提出的斷層模型，介質(zhì)1的主軸電阻率依次為20 Ω·m、20 Ω·m、100 Ω·m，介質(zhì)2的主軸電阻率依次為1 000 Ω·m、1 000 Ω·m、20 Ω·m。該斷層模型，TE模式觀測結(jié)果只與主軸電阻率σ11有關，即該各向異性斷層模型的TE模式觀測結(jié)果與以σ11為電導率的各向同性介質(zhì)斷層模型觀測結(jié)果相同，應用MT2DFEANI計算了該模型在周期T=10 s時的ρyx和φyx。斷層模型的邊界設置為180 km，上下邊界頂點的剖分尺寸為80 km，測線兩側(cè)頂點處和異常體剖分尺寸為500 m，網(wǎng)格剖分結(jié)果如圖8所示。本文模擬結(jié)果與文獻[6]的解析解結(jié)果對比如圖9所示。

圖9的結(jié)果可以看出，MT2DFEANI的結(jié)果與解析解的結(jié)果吻合一致。經(jīng)計算，解析解與MT2DFEANI模擬的ρyx和φyx的最大相對誤差分別為2.12%、2.07%，這驗證了本文MT2DFEANI模擬的準確性。

為了進一步檢測本文有限元程序的正確性，筆者參考Pe等[4]設計的一個較為復雜的模型。模型如圖10所示，具有電阻率水平各向異性的水平地層位于一個二維異常塊體下方，該二維體出露至地表，且其電阻率亦是水平各向異性的，主軸電阻率依次為30 Ω·m、100 Ω·m、30 Ω·m，其各向異性主軸x′與走向方向x軸之間的夾角為30°，而水平地層的各向異性主軸x′與走向方向x軸之間的夾角為120°，主軸

圖5 高低高模型模擬結(jié)果Fig.5 Results for 1D D-model

圖6 低高低模型模擬結(jié)果Fig.6 Results for 1D H-model

電阻率依次為10 Ω·m、100 Ω·m、10 Ω·m，即二維異常體的各向異性水平主軸與水平地層的各向異性主軸相互垂直。該模型的外邊界設置為230 km，上下邊界頂點的剖分尺寸為80 km，測線兩側(cè)頂點剖分尺寸為2 km，異常體剖分尺寸為200 m，網(wǎng)格剖分結(jié)果如圖11所示，異常體剖分放大圖如圖圖12所示。MT2DFEANI模擬的結(jié)果和文獻[4]的有限差分模擬結(jié)果對比如圖13所示，模擬周期為T=30 s。從圖13的對比結(jié)果可以看出，在復雜各向異性的介質(zhì)中，本文的算法模擬結(jié)果和有限差分的模擬結(jié)果吻合一致。視電阻率ρxx、ρxy、ρyx、ρyy和φxx、φxy、φyx、φyy最大相對誤差均控制在2.5%以內(nèi)，驗證了MT2DFEANI模擬結(jié)果的正確性。

圖7 各向異性半空間內(nèi)無限延伸的斷層模型Fig.7 Anisotropic semi-infinite fault model

圖8 斷層模型的非結(jié)構(gòu)網(wǎng)格剖分Fig.8 Grid for the fault model

圖9 斷層模型的解析解和本文有限元結(jié)果的對比圖Fig.9 Comparison between analytic results and the MT2DFEANI results for the fault model

圖10 Pek and Verner設計的二維各向異性模型Fig.10 Two dimensional anisotropic model designed by Pek and Verner

圖11 Pek and Verner設計模型的網(wǎng)格剖分圖Fig.11 Grid for model by Pek and Vernerr

圖12 二維各向民性模型Fig.12 Two dimensional anisotropic model

圖13 模型的有限差分(FD)和本文 有限元(FE)模擬結(jié)果Fig.13 Results comparison between finite- difference and finite element algorithm

3 各向異性影響研究

模型為含有一個二維各向異性異常體的半無限空間介質(zhì)，異常體埋深2 km，半無限空間介質(zhì)的電阻率為1 000 Ω·m。

3.1 傾斜各向異性

傾斜各向異性介質(zhì)是y,z軸在yoz平面內(nèi)繞x軸旋轉(zhuǎn)αd角度得到，其電阻率張量與主軸電阻率張量和αd的關系式表達如下：

(19)

將表達式(19)代入式(2)和式(3)，電場和磁場耦合在一起的方程組可以化簡為：

▽2Ex+iωμ0σ11Ex=0

(20)

(21)

由式(20)可以看出，在傾斜各向異性情況下，TE模式已經(jīng)解耦，且其求解方式與電阻率為σ11的各向同性介質(zhì)求解方式完全相同。式(21)表達的TM模式也已經(jīng)解耦，但是其仍然較為復雜，不能用各向同性介質(zhì)求解方式求解。傾斜各向異性的模型如圖14所示，各向異性異常體的主軸電阻率ρx'、ρy'、ρz'，依次為300 Ω·m、10 Ω·m、300 Ω·m。為討論αd對大地電磁響應的影響程度，我們設置αd依次為0°、30°、45°、60°、90°，該模型在T=10 s下的大地電磁響應如圖15所示。從圖15我們可以得到幾個結(jié)論：①ρxy和φxy不受αd的影響，僅與σ11有關，這與偏微分方程的表達一致；②ρyx和φyx會隨著αd變化，且αd越大，ρyx和φyx的變化也越大，其曲線也不再是關于中心對稱；③當αd=90°時，介質(zhì)轉(zhuǎn)變成垂直各向異性介質(zhì)，ρyx和φyx與垂直各向異性介質(zhì)響應一致，曲線又恢復了中心對稱形式。

圖14 傾斜各向異性模型Fig.14 Dipping anisotropic model

3.2 主軸各向異性

▽2Ex+iωμ0σ11Ex=0

(22)

圖15 傾斜各向異性模型的MT響應Fig.15 Modeling results for dipping anisotropic model

圖16 TIH和TIV介質(zhì)的視電阻率ρyx和相位φyx Fig.16 Apparent and phase results for the TIH and TIV media

(23)

對于TIH介質(zhì)，隨著ρ2的增大，ρyx和φyx逐漸趨于平緩，直到當ρ2=1 000 Ω·m時，此時ρyx和φyx變成一條直線，即TIH介質(zhì)中當ρ2與圍巖介質(zhì)電阻率相同時，此時無法檢測到異常電阻率體。TIV介質(zhì)的ρyx和φyx也隨著ρ2的增大趨于平緩，但是相對于TIH介質(zhì)的曲線變化微乎其微。根據(jù)以上分析得出以下幾個結(jié)論：①在主軸各向異性情況下，相較于σ33而言，ρyx和φyx受σ22的影響更大；②ρyx和φyx不受σ11的影響，這一點由式(23)可知；③ρxy和φxy只與σ11有關。

3.3 水平各向異性

水平各向異性是x軸和y軸在xoy平面內(nèi)繞z軸旋轉(zhuǎn)一定角度αs，其電阻率張量與主軸電阻率張量和αs的關系式表達如下：

將其電阻率張量的表達式代入式(2)和式(3)，方程組可以化簡為：

(24)

(25)

由式(23)和式(24)可以看到Ex和Hx依然耦合在一起，此時的電場Ex不僅和σ11、σ12、σ22有關，且阻抗張量中的Zxx和Zyy為非零元。水平各向異性的模型如圖17所示，各向異性異常體的主軸電阻率依次為300 Ω·m、10 Ω·m、300 Ω·m，為討論αs對大地電磁響應的影響程度，設置αs依次為0°、30°、45°、60°、90°。該模型在T=10 s下的大地電磁響應如圖18所示，圖19是其阻抗張量圖。從圖18可以看出，隨著αs的變化，ρxy、φxy、ρyx和φyx均發(fā)生變化。從圖19的阻抗張量圖中可以得出三個結(jié)論，一是隨著觀測點逐漸遠離異常體，阻抗Zxx和Zyy趨于零，在異常體中心位置出現(xiàn)最大值，二是在αs從0°增長到90°的過程中，Zxx和Zyy先增大后減小，且在αs=0°和αs=90°是均為零，三是αs從0°增長到90°的過程中，Zxy和Zyx在異常體附近減小，而在αs=90°時，Zxy和Zyx曲線出現(xiàn)較大變化，因為此時σ11和σ22發(fā)生互換。

圖17 水平各向異性模型Fig.17 Horizontal anisotropic model

圖18 水平各向異性模型的大地電磁響應Fig.18 Modeling results for horizontal anisotropic model

4 結(jié)論

通過本文算法的推導、驗證及各向異性參數(shù)的分析研究，可以得到幾點結(jié)論。

1)通過兩個一維及兩個二維各向異性模型算例，驗證了本文的二階插值非結(jié)構(gòu)有限元算法的正確性。通過一維模型算例中，與一次插值非結(jié)構(gòu)有限元結(jié)果的對比，驗證了本文二次插值非結(jié)構(gòu)有限元算法的優(yōu)勢性。

圖19 水平各向異性模型的阻抗張量圖Fig.19 Impedance results for horizontal anisotropic model

2)對于較為特殊的傾斜各向異性介質(zhì)和主軸各向異性介質(zhì)而言，其電場和磁場已經(jīng)完全解耦，其阻抗分量中的Zxx和Zyy仍然為零，這與各向同性介質(zhì)相同，但是其Zxy≠Zyx。其TE模式只依賴于σ11。而其TM模式雖已經(jīng)解耦，但是由于各向異性電阻率張量的存在，其磁場的傳播與多個電阻率分量有關，使得其TM模式求解不再像各向同性介質(zhì)一樣簡單，也是其視電阻率和相位有著復雜的變化。

電阻率各向異性現(xiàn)象非常復雜，但卻符合真實的地質(zhì)模型，因此在后續(xù)工作中，將進一步討論電阻率各向異性對大地電磁觀測的影響。