電阻率各向異性介質(zhì)大地電磁二維非結(jié)構(gòu)有限元數(shù)值模擬

惠 鑫， 吳小平*

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

0 引言

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

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

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

1 有限元理論分析

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

(1)

(2)

(3)

(4)

其中：

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

D=σ22σ33-σ23σ32

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

1.1 插值函數(shù)

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

i=1,…,6

(6)

其中:Ei和Hi分別是全局坐標(biāo)系yoz中三角形單元的第i個(gè)節(jié)點(diǎn)的電場(chǎng)和磁場(chǎng);二次差值的插值函數(shù)Ni采用文獻(xiàn) [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 邊界條件

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

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

1.3 單元分析

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

(9)

(10)

公式(9)和公式(10)中的Γe為單元e的邊界，ey和ez是方向單位向量。在內(nèi)邊界Γe上，電場(chǎng)的切向分量Et和磁場(chǎng)的切向分量Ht連續(xù)。在求積分的過(guò)程中，沿著每一個(gè)邊界都進(jìn)行了兩次積分，積分的方向相反，因而，所有的內(nèi)部邊界積分之和為零。又因?yàn)樵谕膺吔缟衔覀儾捎肈irichlet邊界條件，故電場(chǎng)的變分δEx和磁場(chǎng)的變分δHx為零，故所有線性積分之和等于零，公式(9)和公式(10)可以進(jìn)一步化簡(jiǎn)為式(11)與式(12)。

(11)

▽?duì)腍xdΩ=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àn)樗鼈儾辉僖蕾?lài)于y和z)，可得到式(13)和式(14)。

(13)

(14)

KU=0

(15)

其中，

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

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

圖4 矩陣稀疏Fig.4 Sparse matrix

1.5 視電阻率的計(jì)算

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

(16)

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

(17)

2 算法驗(yàn)證

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

(18)

2.1 一維模型驗(yàn)證

一維層狀各向異性模型的大地電磁響應(yīng)，存在解析解，因此設(shè)了一個(gè)三層的電阻率參數(shù)為“高低高”的模型和一個(gè)三層的電阻率參數(shù)為“低高低”的模型。兩個(gè)模型的第一層和第二層深度均為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°。

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

2.2 二維模型驗(yàn)證

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

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

為了進(jìn)一步檢測(cè)本文有限元程序的正確性，筆者參考Pe等[4]設(shè)計(jì)的一個(gè)較為復(fù)雜的模型。模型如圖10所示，具有電阻率水平各向異性的水平地層位于一個(gè)二維異常塊體下方，該二維體出露至地表，且其電阻率亦是水平各向異性的，主軸電阻率依次為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，即二維異常體的各向異性水平主軸與水平地層的各向異性主軸相互垂直。該模型的外邊界設(shè)置為230 km，上下邊界頂點(diǎn)的剖分尺寸為80 km，測(cè)線兩側(cè)頂點(diǎn)剖分尺寸為2 km，異常體剖分尺寸為200 m，網(wǎng)格剖分結(jié)果如圖11所示，異常體剖分放大圖如圖圖12所示。MT2DFEANI模擬的結(jié)果和文獻(xiàn)[4]的有限差分模擬結(jié)果對(duì)比如圖13所示，模擬周期為T(mén)=30 s。從圖13的對(duì)比結(jié)果可以看出，在復(fù)雜各向異性的介質(zhì)中，本文的算法模擬結(jié)果和有限差分的模擬結(jié)果吻合一致。視電阻率ρxx、ρxy、ρyx、ρyy和φxx、φxy、φyx、φyy最大相對(duì)誤差均控制在2.5%以?xún)?nèi)，驗(yàn)證了MT2DFEANI模擬結(jié)果的正確性。

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

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

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

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

圖11 Pek and Verner設(shè)計(jì)模型的網(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 各向異性影響研究

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

3.1 傾斜各向異性

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

(19)

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

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

(20)

(21)

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

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

3.2 主軸各向異性

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

(22)

圖15 傾斜各向異性模型的MT響應(yīng)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)

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

3.3 水平各向異性

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

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

(24)

(25)

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

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

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

4 結(jié)論

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

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

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

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

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