王 寧， 肖 曉， 張逗逗

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083；2.中南大學(xué) 有色金屬成礦預(yù)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室, 長沙 410083)

起伏地形條件下二維大地電磁各向異性正演

王 寧1,2， 肖 曉1,2， 張逗逗1,2

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083；2.中南大學(xué) 有色金屬成礦預(yù)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室, 長沙 410083)

在大地電磁資料處理和解釋中，大地電磁的各向異性正演一直是國內(nèi)、外研究的前沿課題。這里首先從麥克斯韋方程組出發(fā)，得出各向異性介質(zhì)二維大地電磁場(chǎng)的邊值問題以及等價(jià)的變分問題，進(jìn)而對(duì)計(jì)算區(qū)域進(jìn)行完全非結(jié)構(gòu)化三角形剖分，在單元內(nèi)采用線性插值，將變分方程轉(zhuǎn)化成線性方程組，求解出有限單元法數(shù)值解，獲得各向異性介質(zhì)下的電磁場(chǎng)值。討論了電性主軸與笛卡爾坐標(biāo)軸之間的夾角以及地形的起伏變化對(duì)視電阻率和相位值的影響，結(jié)果表明：起伏地形條件下，TE和TM兩種極化模式的視電阻率和相位值都會(huì)受到影響，且TM模式受地形影響較大；由于介質(zhì)的各向異性，視電阻率和相位斷面圖出現(xiàn)明顯不同于各向同性的形態(tài)，導(dǎo)致橫向分界面模糊。

大地電磁； 有限單元法； 各向異性； 起伏地形

0 引言

大地電磁測(cè)深法(Magnetotelluric Sounding Method)是利用天然電磁場(chǎng)源，研究地殼及上地幔電性結(jié)構(gòu)的一種重要的地球物理勘探方法[1-2]。但由于地球構(gòu)造應(yīng)力場(chǎng)、地球形變帶、巖石裂隙、空隙水及地質(zhì)沉積等因素造成的地球各向異性問題[3]，對(duì)大地電磁勘探數(shù)據(jù)準(zhǔn)確合理地解釋造成極大的困難，因而各向異性研究備受國內(nèi)、外地球物理學(xué)家的關(guān)注。O'Brien 等[4]推導(dǎo)了層狀各向異性電磁場(chǎng)計(jì)算的遞推公式，由此開啟了大地電磁各向異性研究的大門 ； Pek等[5]推出了一種聯(lián)合計(jì)算一般層狀各向異性介質(zhì)中大地電磁阻抗張量和參數(shù)靈敏度矩陣的算法。在大地電磁二維各向異性正演問題中，有限單元法、有限差分法和特殊情況下解析解法是主要的研究方法。有限差分方面， Saraf 等[6]用有限差分法計(jì)算了均勻半空間中嵌入一個(gè)水平和垂直方向電導(dǎo)率不同的異常體時(shí)的H偏振模式； Pek等[7]推導(dǎo)出了全張量電導(dǎo)率的電磁場(chǎng)偏微分方程，并且應(yīng)用有限差分法對(duì)任意各向異性介質(zhì)進(jìn)行了正演計(jì)算；霍光譜等[8]利用Pek[7]的求解方法研究了起伏地形條件下，各向異性對(duì)大地電磁實(shí)測(cè)資料的影響。但傳統(tǒng)有限差分采用矩形網(wǎng)格，這將導(dǎo)致該方法在擬合地形及復(fù)雜地質(zhì)構(gòu)時(shí)不能完成很好的模擬。Reddy 等[9]率先使用伽遼金有限元法對(duì)二維各向異性問題中的偏微分方程做了數(shù)值計(jì)算；徐世浙等[10]用有限元法計(jì)算了二維各向異性地電斷面的大地電磁場(chǎng)；李予國[11]在前人的基礎(chǔ)上給出了任意各向異性地質(zhì)體的感應(yīng)電磁場(chǎng)的有限元離散公式；熊彬等[12-13]用有限元法模擬了二維各向異性MT響應(yīng)，并且討論了層面傾角和頻率對(duì)低電阻率的影響；秦林江等[14-15]采用Fourier變化的方法獲得了對(duì)角各向異性無限深斷裂的MT響應(yīng)函數(shù)的解析解，并對(duì)解析解計(jì)算過程中的精度問題進(jìn)行了討論，同時(shí)對(duì)實(shí)測(cè)電磁資料中遇到的兩支視電阻率在無窮遠(yuǎn)處不一致的現(xiàn)象給出了解釋；田紅軍等[16]討論了起伏地形和各向異性對(duì)大地電磁測(cè)深中TM極化模式響應(yīng)的影響規(guī)律。

目前國內(nèi)各向異性的正演模擬中大多使用的是結(jié)構(gòu)化網(wǎng)格[12-16]，但是該網(wǎng)格在處理復(fù)雜構(gòu)造和地形問題時(shí)很難完成精確地模擬。筆者采用非結(jié)構(gòu)化網(wǎng)格[18-21]有限單元法，對(duì)大地電磁中常見的電性主軸其中一個(gè)方向與走向一致的二維各向異性體結(jié)構(gòu)[10,12-16]進(jìn)行正演模擬，研究視電阻率和相位值隨其他兩個(gè)電性主軸與笛卡爾坐標(biāo)軸夾角α的變化規(guī)律，以及地形對(duì)視電阻率和相位值的影響。

1 變分問題

1.1 基本理論

由電磁場(chǎng)理論可知，電磁波在介質(zhì)中的傳播遵從麥克斯韋方程組[1-2,17]：

▽×E=iωμH

(1)

▽×H=(σ-iωμε)E

(2)

式中：E、H分別是電場(chǎng)強(qiáng)度矢量和磁場(chǎng)強(qiáng)度矢量；ω=2πf是角頻率，對(duì)于大地電磁法所涉及的頻率，巖石中σ>>ωε，因而可以忽略介電常數(shù)的各向異性，筆者只考慮電導(dǎo)率的各向異性，因而電導(dǎo)率σ是張量?？諝庖约暗叵陆橘|(zhì)的導(dǎo)磁率和介電常數(shù)分別取μ=μ0=4π×10-7，ε=ε0=8.85×10-12。對(duì)于各向異性介質(zhì)，可以找到一個(gè)坐標(biāo)系ikm(圖1)，使張量電導(dǎo)率化為對(duì)角矩陣形式：

(3)

其中：σi、σk、σm分別為電性主軸方向上的電導(dǎo)率。

圖1 各向異性主軸與笛卡爾坐標(biāo)系Fig.1 Anisotropic principal direction and Cartesian coordinate system

圖2 邊界條件Fig.2 The boundary conditions

我們考慮的各向異性模型的電性主軸i與笛卡爾坐標(biāo)系中的x方向一致，而互相正交的電性主軸k、m與y、z坐標(biāo)軸成一角度α。按照張量在不同坐標(biāo)系中的變換公式，在xyz坐標(biāo)系中，電導(dǎo)率張量σ可表示為：

σ=Aσ′AT

其中:A為坐標(biāo)變換張量：

(4)

TE極化模式：

(5)

TM極化模式：

(6)

式(5)和式(6)可以統(tǒng)一寫成式(7)。

▽·(τ▽u)+λu=0

(7)

式中：▽為二維哈密頓算符。

1.2 邊界條件

二維各向異性介質(zhì)大地電磁場(chǎng)兩種極化模式下滿足的邊值問題為：

(8)

1.3 變分問題

與邊值問題等價(jià)的變分問題為式(9)～式(11)。

(9)

u|AB=1

(10)

δF(u)=0

(11)

2 有限單元法

筆者采用開源程序triangle將求解區(qū)域Ω自動(dòng)剖分成三角形單元，對(duì)于三角單元e，采用線性基函數(shù)，u表示為式(12)。

u=Niui+Njuj+Nmum

(12)

其中：

式中(yk,zk)(k=i,j,m)為三角形單元頂點(diǎn)的坐標(biāo)。

最終得到的線性方程組為式(13)。

Ku=0

(13)

采用乘大數(shù)法，加載式(8)中AB邊界條件，形成右端項(xiàng)，得到大型稀疏方程組。對(duì)方程組采用Eigen中的直接求解器PardisoLU來求解線性方程組。

為了計(jì)算視電阻率還需要知道地面測(cè)點(diǎn)上的Ey和Hy，由式(5)、式(6)有：

(14)

(15)

采用四點(diǎn)差分求得相應(yīng)的磁場(chǎng)Hy和電場(chǎng)Ex，即可求出視電阻率和相位值：

TE模式：

TM模式：

3 正演結(jié)果分析

3.1 程序驗(yàn)證

筆者采用c語言編寫了相應(yīng)的正演程序。為了驗(yàn)證程序的精度和可靠性， 以標(biāo)準(zhǔn)模型COMMEMI 2D-1作為驗(yàn)證模型(圖3)，選用頻率f=10 Hz，與開源程序MT2D(https:// sourceforge.net/projects/mt2d/)計(jì)算結(jié)果進(jìn)行對(duì)比，結(jié)果如圖4所示。TE模式下，二者視電阻率相對(duì)誤差在0.22%以內(nèi)，相位差在0.35°以內(nèi)；TM模式下，二者視電阻率相對(duì)誤差在0.34%以內(nèi)，相位差在0.39°以內(nèi)。

圖3 2D-1模型Fig.3 The 2D-1 model

圖4 程序結(jié)果驗(yàn)證Fig.4 Comparison of the results with MT2D(a)視電阻率對(duì)比圖；(b)相位對(duì)比圖

3.2 不同傾角對(duì)MT響應(yīng)的影響

圖5 二維各向異性模型Fig.5 The two-dimensional anisotropic model

從圖6可以看出，TE模式視電阻率和相位均不隨夾角α而變。但TM模式視電阻率和相位值隨著夾角α變化有很大的差異：在α=0°和α=90°時(shí)，視電阻率和相位值左右兩端對(duì)稱；在α=30°和α=60°時(shí)，視電阻率和相位值左右兩端不對(duì)稱。在無窮遠(yuǎn)處，不同夾角下的視電阻率值均趨于均勻半空間的值，而相位值趨于45°。

3.3 在起伏地形下不同夾角的對(duì)比模型

由圖6可知，α對(duì) TE模式的視電阻率和相位值都沒有影響，故圖8只給出α=0°時(shí)三種模型的視電阻率和相位擬斷面圖。圖8(a)是無地形時(shí)視電阻率擬斷面圖，圖中顯示有一個(gè)低阻的閉合；圖8(b)由于地壘模型的存在，使得地壘頂部覆蓋范圍內(nèi)視電阻率值有所增加；圖8(c)由于地塹模型的存在，使得地塹底部覆蓋范圍內(nèi)電阻率值有所減小。從圖8(d)～圖8(f)中可以看出，高頻段相位高于低頻段，地壘模型減小了高低頻段的相位差值，地塹增大了這種差異。

圖9是TM模式下的視電阻率和相位擬斷面圖。從圖9(a)中可以看出：當(dāng)α=0°和α=90°時(shí)，異常體橫向分界面明顯，是左右對(duì)稱的低阻異常體；當(dāng)α=30°和α=60°時(shí)，異常體橫向分界面不再明顯，在左側(cè)出現(xiàn)低阻的異常，而右側(cè)，顯示高阻的異常；在地壘頂部覆蓋范圍內(nèi)，視電阻率明顯減小，在地壘與平面交界處有明顯的高阻異常，在地塹底部覆蓋范圍內(nèi)，視電阻率明顯增大，而在地塹與平地面交界處有明顯的低阻異常。從圖9(b)中可以看出：α=0°和α=90°時(shí)，相位橫向分界面明顯，且是左右對(duì)稱；當(dāng)α=30°和α=60°時(shí)，在左側(cè)顯示高相位，而右側(cè)顯示低相位。地形對(duì)相位也有一定的影響，但是影響的數(shù)值比較小。

圖6 不同夾角對(duì)應(yīng)的視電阻率和相位值(頻率為1 Hz)Fig.6 Corresponding the values of apparent resistivity and phase to different angles(with the frequency of 1 Hz)(a)TE模式下視電阻率；(b)TM模式下視電阻率；(c)TE模式下相位；(d)TM模式下相位

圖7 帶有起伏地形的二維各向異性模型Fig.7 2D anisotropic model with topography

圖8 TE模式視電阻率和相位擬斷面圖Fig.8 Apparent resistivity and phase cross-section diagram of TE(a)水平地形視電阻率擬斷面圖;(b)地壘地形視電阻率擬斷面圖;(c)地塹地形視電阻率擬斷面圖;(d)水平地形相位擬斷面圖;(e)地壘地形相位擬斷面圖;(f)地塹地形相位擬斷面圖

圖9 TM模式視電阻率和相位擬斷面圖Fig.9 Apparent resistivity and phase cross-section diagram of TM(a)視電阻率擬斷面圖；(b)相位擬斷面圖

4 結(jié)論

1)當(dāng)α改變時(shí)，TE極化模式的視電阻率和相位保持不變，TM極化模式視電阻率和相位發(fā)生明顯變化：在α=0°和α=90°時(shí)，視電阻率和相位擬斷面圖是對(duì)稱的，且橫向分界面明顯；在α≠0°且α≠90°時(shí)，視電阻率和相位形態(tài)不再嚴(yán)格對(duì)稱，甚至出現(xiàn)很大的差異；橫向分界面也變得比較模糊。

2)起伏地形條件下，TE和TM兩種極化模式的視電阻率和相位都會(huì)受到影響，但地形對(duì)TM模式的影響較大，TE模式則受地形影響微弱。

3)對(duì)于TE極化模式，當(dāng)電性主軸其中一個(gè)方向與走向一致時(shí)可以不考慮主軸夾角變化的影響，但對(duì)于TM極化模式，當(dāng)起伏地形和各向異性同時(shí)存在時(shí)，二者的影響都應(yīng)該考慮。

[1] TIKHONOV A N. On determining electric characteristics of the deep layers of the Earth's crust[J]. Sov.math.dokl, 1950, 73：295-297.

[2] CAGNIARD L. Basic theory of the magneto-telluric method of geophysical prospecting[J]. Geophysics, 2002, 18(3):605-635.

[3] 霍光譜, 胡祥云, 劉敏. 各向異性介質(zhì)中大地電磁正演研究綜述[J]. 地球物理學(xué)進(jìn)展, 2011, 26(6): 1976-1982.

HUO G P,HU X Y,LIU M. Review of the forward modeling of magnetotelluric in the anisotropy medium research[J]. Progress in Geophysics, 2011, 26(6): 1976-1982.( In Chinese)

[4] O'BRIEN D P,MORRISON H F. Electromagnetic fields in an n-layer anisotropic half-space[J].Geophysics,1967,32(4):668-677.

[5] PEK J, SANTOS F A M. Magnetotelluric impedances and parametric sensitivities for 1-D anisotropic layered media[J]. Computers & geosciences, 2002, 28(8): 939-950.

[6] SARAF P D, NEGI J G, ERV V. Magnetotelluric response of a laterally inhomogeneous anisotropic inclusion[J]. Physics of the earth and planetary interiors, 1986, 43(3): 196-198.

[7] PEK J, VERNER T. Finite-difference modelling of magnetotelluric fields in two-dimensional anisotropic media[J]. Geophysical Journal International, 1997, 128(3): 505-521.

[8] 霍光譜,胡祥云,黃一凡,等. 帶地形的大地電磁各向異性二維模擬及實(shí)例對(duì)比分析[J]. 地球物理學(xué)報(bào),2015(12):4696-4708.

HUO G P, HU X Y, HUANG Y F,et al. MT modeling for 2D anisotropic conductivity structure with topography and examples of comparative analyse [J]. Chinese Journal of Geophysics,2015(12):4696-4708. ( In Chinese)

[9] REDDY I K, RANKIN D. Magnetotelluric response of laterally inhomogeneous and anisotropic media[J]. Geophysics, 1975, 40(6): 1035-1045.

[10] 徐世浙, 趙生凱. 二維各向異性地電斷面大地電磁場(chǎng)的有限元法解法[J]. 地震學(xué)報(bào), 1985, 7(1): 80-89.

XU S Z, ZHAO S K. Solution of magnetotelluric field equations for a two-dimensional, anisotropic geoelectric section by the finite element method[J].Acta Seismologica Sinica, 1985, 7(1): 80-89. ( In Chinese)

[11] LI Y. A finite-element algorithm for electromagnetic induction in two-dimensional anisotropic conductivity structures[J]. Geophysical Journal International, 2002, 148(3): 389-401.

[12] 熊彬, 羅天涯, 李長偉, 等. 用有限單元法模擬各向異性介質(zhì)中二維電性異常體的大地電磁響應(yīng)[J]. 中山大學(xué)學(xué)報(bào): 自然科學(xué)版, 2013, 52(4): 143-148.

XIONG B, LUO T Y, LI C W, et al. Modelling of magnetotelluric responses of 2-D electrical anomalous body in anisotropic media using finite element method[J]. Acta Scientiarum Naturalium Universitatis Sunyatseni, 2013, 52(4): 143-148. ( In Chinese)

[13] 羅天涯, 熊彬, 李長偉, 等. 各向異性介質(zhì)中大地電磁場(chǎng)的異常特征[J]. 物探化探計(jì)算技術(shù), 2013, 35(6): 645-650.

LUO T Y, XIONG B, LI C W, et al. The abnormal charcteristic of magnetotelluric in anisotrpy media [J].Computing Techniques for Geophysical and Geochemical Exploration, 2013, 35(6): 645-650. ( In Chinese)

[14] QIN L, YANG C, CHEN K. Quasi-analytic solution of 2-D magnetotelluric fields on an axially anisotropic infinite fault[J]. Geophysical Journal International, 2013, 192(1): 67-74.

[15] QIN L J, YANG C F, CHEN K. Analytic solution to the magnetotelluric response over anisotropic medium and its discussion[J]. Science China Earth Sciences, 2013, 56(9): 1607-1615.

[16] 田紅軍, 佟鐵鋼. 帶地形的二維各向異性大地電磁測(cè)深數(shù)值模擬[J]. 工程地球物理學(xué)報(bào), 2015, 12(2): 139-144.

TIAN H J,TONG T G. The forward numerical simulation of 2D anisotropic magnetotelluric sounding[J].Chinese Journal of Engineering Geophysics,2015, 12(2): 139-144. (In Chinese)

[17] 徐世浙.地球物理中的有限單元法[M].北京：科學(xué)出版社，1994.

XU S Z. The finite element method in geophysics [M].Beijing:The Science press,1994.( In Chinese)

[18] SHEWCHUK J R. Triangle: Engineering a 2D quality mesh generator and Delaunay triangulator[C]// Selected papers from the Workshop on Applied Computational Geormetry, Towards Geometric Engineering. Springer-Verlag, 1996:203-222.

[19] REN Z, TANG J. 3D direct current resistivity modeling with unstructured mesh by adaptive finite-element method[J]. Geophysics, 2010, 75(1):H7.

[20] REN Z, KALSCHEUER T, GREENHALGH S, et al. A goal-oriented adaptive finite-element approach for plane wave 3-D electromagnetic modelling[J]. Geophysical Journal International, 2012, 194(2):700-718.

[21] 原源, 湯井田, 任政勇,等. 基于非結(jié)構(gòu)化網(wǎng)格的任意復(fù)雜2D RMT有限元模擬[J]. 地球物理學(xué)報(bào), 2015(12):4685-4695.

YUAN Y，TANG J T，REN Z Y，et al. Two-dimensional complicated radio-magnetotelluric finite-element modeling using unstructured grid[J]. Chinese Journal of Geophysics, 2015(12):4685-4695. ( In Chinese)

Theforwardmodelingoftwo-dimensionalanisotropicmagnetotelluricwithtopography

WANG Ning, XIAO Xiao, ZHANG Doudou

(1.Central South University School of Geosciences and Info-Physics, Changsha 410083,China;2.Central South University Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education, Changsha 410083,China)

In the interpretation of magnetotelluric data processing, the phenomenon of anisotropy and topography are the important factors that cannot be ignored. The forward modeling of anisotropic magnetotelluric has always been the forefront subject at home and abroad. In this paper, we given the boundary value problems and its corresponding variational problem of magnetotelluric in 2D anisotropic geoelectric section from the Maxwell's equations. A triangular grid is used to discretize the calculation area, and linear interpolation is used at each triangular to compute the numerical solution with finite element method obtained the apparent resistivity and phase values. In the subsequent section we discuss the effect of the angle between the electrical spindle and Cartesian coordinate axis and topography.

magnetotelluric; FEM; anisotropic; topography

2016-10-20 改回日期： 2016-11-09

國家自然科學(xué)基金(41174105)；國家高技術(shù)研究發(fā)展計(jì)劃(2014AA06A602)

王寧(1990-)，男，碩士，主要從事電磁場(chǎng)理論和應(yīng)用研究， E-mail：w_ning2016@163.com。

肖曉(1981-)，男，博士后，副教授，碩士生導(dǎo)師，主要從事電磁場(chǎng)理論和應(yīng)用研究， E-mail：csuxiaox@csu.edu.cn。

1001-1749(2017)06-0711-08

P 631.2

A

10.3969/j.issn.1001-1749.2017.06.01