電阻率各向異性介質大地電磁二維非結構有限元數值模擬

惠 鑫， 吳小平*

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

0 引言

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

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

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

1 有限元理論分析

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

(1)

(2)

(3)

(4)

其中：

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

D=σ22σ33-σ23σ32

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

1.1 插值函數

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

i=1,…,6

(6)

其中:Ei和Hi分別是全局坐標系yoz中三角形單元的第i個節點的電場和磁場;二次差值的插值函數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 邊界條件

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

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

1.3 單元分析

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

(9)

(10)

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

(11)

▽δHxdΩ=0

(12)

圖1 二維非結構單元Fig.1 2D unstructured element

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

1.4 剛度矩陣

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

(13)

(14)

KU=0

(15)

其中，

本文中調用Fortran中的pardiso庫函數，運用直接法求解系數矩陣。

圖3 網格剖分實例Fig.3 Grid example

圖4 矩陣稀疏Fig.4 Sparse matrix

1.5 視電阻率的計算

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

(16)

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

(17)

2 算法驗證

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

(18)

2.1 一維模型驗證

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

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

2.2 二維模型驗證

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

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

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

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

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

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

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

圖8 斷層模型的非結構網格剖分Fig.8 Grid for the fault model

圖9 斷層模型的解析解和本文有限元結果的對比圖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設計模型的網格剖分圖Fig.11 Grid for model by Pek and Vernerr

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

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

3 各向異性影響研究

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

3.1 傾斜各向異性

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

(19)

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

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

(20)

(21)

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

圖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介質的視電阻率ρyx和相位φyx Fig.16 Apparent and phase results for the TIH and TIV media

(23)

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

3.3 水平各向異性

水平各向異性是x軸和y軸在xoy平面內繞z軸旋轉一定角度α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均發生變化。從圖19的阻抗張量圖中可以得出三個結論，一是隨著觀測點逐漸遠離異常體，阻抗Zxx和Zyy趨于零，在異常體中心位置出現最大值，二是在αs從0°增長到90°的過程中，Zxx和Zyy先增大后減小，且在αs=0°和αs=90°是均為零，三是αs從0°增長到90°的過程中，Zxy和Zyx在異常體附近減小，而在αs=90°時，Zxy和Zyx曲線出現較大變化，因為此時σ11和σ22發生互換。

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

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

4 結論

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

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

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

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

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