霍光譜, 胡祥云, 黃一凡, 韓波

1 河南省地質(zhì)調(diào)查院， 鄭州 4500002 中國(guó)地質(zhì)大學(xué)地球物理與空間信息學(xué)院，地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 430074

?

帶地形的大地電磁各向異性二維模擬及實(shí)例對(duì)比分析

霍光譜1, 胡祥云2*, 黃一凡2, 韓波2

1 河南省地質(zhì)調(diào)查院， 鄭州 4500002 中國(guó)地質(zhì)大學(xué)地球物理與空間信息學(xué)院，地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室， 武漢 430074

帶地形的大地電磁二維正演數(shù)值模擬多數(shù)基于電性各向同性理論，由于地球內(nèi)部電性各向異性現(xiàn)象的普遍存在，基于電性各向異性理論研究地形起伏情況下大地電磁二維正演數(shù)值模擬就顯得非常迫切.本文首先由麥克斯韋方程出發(fā)，引入張量電導(dǎo)率，求得一組關(guān)于平行走向的電場(chǎng)分量Ex和磁場(chǎng)分量Hx的二階偏微分方程，使用有限差分法求解出Ex和Hx的近似解，并以此求得其他場(chǎng)分量；其次，引入地形因素，改變變量在網(wǎng)格節(jié)點(diǎn)中的排列方式，選擇交錯(cuò)排列方式從而給有限差分系數(shù)矩陣的最大帶寬分配合理的存儲(chǔ)空間；最后，使用Weaver的方法解決TM模式下，在地-空分界面垂直于構(gòu)造走向的一些區(qū)域存在不同電導(dǎo)率的問(wèn)題.通過(guò)對(duì)帶地形的二維電性各向異性結(jié)構(gòu)做正演模擬，研究地形因素對(duì)大地電磁響應(yīng)的影響；以電性各向異性理論為基礎(chǔ)，將地形因素引入對(duì)實(shí)測(cè)大地電磁資料的處理中，通過(guò)做二維正演擬合和未引入地形因素的結(jié)果做對(duì)比，說(shuō)明電性各向異性現(xiàn)象的普遍存在，認(rèn)識(shí)地形因素對(duì)觀測(cè)大地電磁場(chǎng)的影響，為今后分析解釋實(shí)測(cè)大地電磁資料包含地形因素和電性各向異性情況提供理論基礎(chǔ)和技術(shù)指導(dǎo).

地形； 電性各向異性； 大地電磁； 有限差分； 張量電導(dǎo)率

1 引言

20世紀(jì)50年代初Tikhonov（1950）和Cagnird（1953）分別提出利用天然電磁場(chǎng)進(jìn)行勘探的方法，從70年代開始大地電磁法進(jìn)入飛速發(fā)展時(shí)期，到目前已經(jīng)成為國(guó)內(nèi)外學(xué)者（Weiss and Newman，2002；趙國(guó)澤等，2004；Wannamaker et al，2005；魏文博等，2006,2010；胡祥云等，2013）研究地球內(nèi)部電性結(jié)構(gòu)最常用的方法之一.Ku等（1973）發(fā)現(xiàn)地形對(duì)大地電磁場(chǎng)的影響，20世紀(jì)80年代開始，不斷有學(xué)者將地形因素引入大地電磁數(shù)值模擬中， Wannamaker等（1986）用有限元法對(duì)起伏地形下二維大地電磁響應(yīng)進(jìn)行數(shù)值模擬.Chouteau和Bouehard（1988）研究了二維大地電磁響應(yīng)的地形校正問(wèn)題.徐世浙等（1985,1992,1997）用有限元和邊界元法計(jì)算起伏地形下水平層狀介質(zhì)的大地電磁響應(yīng)，并對(duì)畸變后的大地電磁響應(yīng)做地形改正，這些研究成果均是基于電性各向同性理論.

20世紀(jì)60年代末，Parkhomenko（1967），Keller（1968）和Hill（1972）分別討論了巖石中的微各向異性問(wèn)題，與此同時(shí)，電性各向異性現(xiàn)象的理論研究也逐漸發(fā)展起來(lái).O′Brienh和 Morrison（1967）研究層狀各向異性介質(zhì)中電磁場(chǎng)的遞推公式；Reddy和Rankin （1971）研究?jī)A斜各向異性對(duì)大地電磁場(chǎng)的影響.Reddy和Rankin（1975）用有限元法對(duì)偏微分方程做數(shù)值計(jì)算，較早對(duì)二維電性各向異性問(wèn)題進(jìn)行了理論研究.雖然在其后的20多年間（直到2000之后才有學(xué)者將其理論改進(jìn)后用于大地電磁實(shí)測(cè)資料處理）未能廣泛運(yùn)用在大地電磁實(shí)測(cè)資料處理中，但其所做研究在理論意義上的影響卻尤為深遠(yuǎn)，不僅讓世人重視電性各向異性現(xiàn)象，更為今后的二維電性各向異性研究奠定了理論基礎(chǔ).Pek和Toh（1997）在Reddy和Rankin（1975）研究的基礎(chǔ)上，并將前人（Haak，1972；Brewitt-Taylor and Wearver，1976；Cerv et al.，1978）在二維電性各向同性介質(zhì)中使用的有限差分法引入求解二維電性各向異性問(wèn)題中，用有限差分法對(duì)偏微分方程做數(shù)值模擬計(jì)算，是正演計(jì)算二維電性各向異性介質(zhì)大地電磁響應(yīng)方法中影響最為深遠(yuǎn)和廣泛的一種求解方法.Li（2002，2008）細(xì)述用有限元法計(jì)算二維電性各向異性結(jié)構(gòu)感應(yīng)電磁場(chǎng)，并于2008年用自適應(yīng)網(wǎng)格剖分法代替以往的網(wǎng)格剖分法，是一次成功的改進(jìn)，取得顯著的效果.霍光譜等（2012）通過(guò)引入權(quán)因子成功改進(jìn)馬奎特反演理論，將基于電性各向同性理論的反演方法擴(kuò)展應(yīng)用到電性各向異性理論.

基于電性各向異性理論帶地形的大地電磁數(shù)值模擬在最近20年間才逐漸發(fā)展起來(lái)，用于實(shí)際資料處理的更為少見.Pek和Toh（2000）將公式擴(kuò)展應(yīng)用在大地含地形和海洋測(cè)深中.Key和Ovall（2011）使用有限元法采用自適應(yīng)非結(jié)構(gòu)性網(wǎng)格剖分方式，系統(tǒng)研究二維大地電磁及海洋可控源電磁法，是最近較為系統(tǒng)而完善的方法.

2 帶地形的大地電磁電性各向異性二維正演理論

本文采用有限差分法參照Pek等（1997,2000,2002）的方法研究帶地形的大地電磁二維電性各向異性正演問(wèn)題.

假設(shè)在笛卡爾坐標(biāo)系中，有一個(gè)二維電性結(jié)構(gòu)，它的構(gòu)造走向平行于x軸，y垂直于x軸，z軸垂直于xy平面，且正向向下.模型由一個(gè)有限區(qū)域組成，在其中的電性各向異性區(qū)域是二維有限結(jié)構(gòu).假設(shè)地表為水平，對(duì)應(yīng)z=0.在地表上空（z<0）為絕緣空氣層.來(lái)自z→-∞的平面波垂直向下傳播，ω=2π/T，T為周期，單位為秒（s）.

諧變電磁場(chǎng)的麥克斯韋方程組為

（1）

（2）

上式為描述電磁場(chǎng)分布的微分方程.e-iω t為時(shí)諧因子，σ為張量電導(dǎo)率.

張量電導(dǎo)率σ為三階對(duì)稱矩陣，可分解為由3個(gè)電導(dǎo)率σ1,σ2,σ3和一個(gè)旋轉(zhuǎn)矩陣，旋轉(zhuǎn)矩陣又可通過(guò)3個(gè)Euler′s空間旋轉(zhuǎn)元素依次旋轉(zhuǎn)角度為αL,αD,αS求得，角度αS,αD,αL分別定義為各向異性介質(zhì)的走向，傾角，傾向（圖1所示）.

（3）

其中R為旋轉(zhuǎn)矩陣元素：

（4）

通過(guò)公式（3）可以表示空間中任何位置的電性各向異性結(jié)構(gòu).

由2-D性質(zhì)?/?x=0，將公式簡(jiǎn)化為只含Ex和Hx的一組二階偏微分方程（胡祥云等，2013）.

電場(chǎng)沿走向時(shí)（即TE模式）：

（5）

磁場(chǎng)沿走向時(shí)（即TM模式）：

圖1 各向異性基本參數(shù)示意圖（依次運(yùn)用三個(gè)Euler′s旋轉(zhuǎn)元素αS,αD,αL可將導(dǎo)電面變化到空間任何位置）

（6）

公式（5）—（6）中，

A=（σxyσyz-σxzσyy）/D,B=（σxyσyz-σxzσzz）/D,

用有限差分法求解公式（5）—（6），求得網(wǎng)格節(jié)點(diǎn)（j,k）上TE模式的近似差分方程，這個(gè)近似差分方程代表公式（5）的有限差分形式：

（7）

公式（7）中，

（p,q）∈{（j-1,k-1）,（j-1,k+1）,（j+1,

k-1）,（j+1,k+1）}

同樣可以求得在網(wǎng)格節(jié)點(diǎn)（j,k）上TM模式下的線性差分方程：

（8）

在公式（8）中，

（p,q）∈{（j-1,k-1）,（j-1,k+1）,（j+1,

k-1）,（j+1,k+1）}.

用有限差分法對(duì)公式（5）—（6）做數(shù)值近似后，須將線性代數(shù)公式（7）—（8）安排在一個(gè)系統(tǒng)，以便進(jìn)一步處理.因?yàn)榘珽x和Hx兩個(gè)變量，且這兩個(gè)變量多數(shù)時(shí)候并不同時(shí)出現(xiàn)在同一個(gè)網(wǎng)格節(jié)點(diǎn)，TE模式的方程需要在全部網(wǎng)格節(jié)點(diǎn)做近似計(jì)算，而TM模式的方程只需在導(dǎo)電底層內(nèi)部做近似計(jì)算.

引入地形因素時(shí)，對(duì)應(yīng)的有限差分網(wǎng)格可通過(guò)起伏地形對(duì)應(yīng)的階梯函數(shù)來(lái)近似.并不影響用有限體積法在中心網(wǎng)格節(jié)點(diǎn)的任何位置對(duì)公式（5）—（6）做近似計(jì)算.將變量排入一個(gè)系統(tǒng)有兩種方式：順序排列和交錯(cuò)排列.本文選取交錯(cuò)排列方法.

圖2a是水平地形，對(duì)應(yīng)的存儲(chǔ)矩陣是2c；2b是起伏地形，對(duì)應(yīng)的存儲(chǔ)矩陣是2d.由圖2b可以看到，在網(wǎng)格中，不同的列，Hx的個(gè)數(shù)是可變的，所以圖2d中有限差分系數(shù)矩陣的帶寬是可變的：帶寬的最寬部分由地形的海拔高度決定；帶寬的最窄部分由地形的盆地深度決定.這也是起伏地形相對(duì)于水平地形不同的地方，帶地形情況下，就必須考慮對(duì)公式做相應(yīng)的有效修改，即：在用高斯消元法對(duì)有限差分方程組求解時(shí)，須為有限差分系數(shù)矩陣的最大帶寬分配相應(yīng)的存儲(chǔ)空間.

與水平地表相比，含起伏地形時(shí)，TM模式在地-空分界面（y方向上）存在一些不同電導(dǎo)率的區(qū)域，如：σj,kσj+1,k+1≠σj+1,kσj,k+1.在這些區(qū)域的節(jié)點(diǎn)上，沒有連續(xù)的電流和連續(xù)的切向電流.為解決這個(gè)問(wèn)題，須將Weaver等（1985,1986）的光滑方法應(yīng)用在這些網(wǎng)格節(jié)點(diǎn)的局部導(dǎo)電區(qū)域.

Weaver等（1985,1986）認(rèn)識(shí)到用矩形網(wǎng)格剖分模型時(shí)，有時(shí)網(wǎng)格節(jié)點(diǎn)四周的4個(gè)單元格會(huì)存在不同的電導(dǎo)率，并給出解決此問(wèn)題的方法：

1） 對(duì)選定的網(wǎng)格節(jié)點(diǎn)周圍的局部電導(dǎo)率重新塑造為光滑的電導(dǎo)率分布，這并不影響用有限差分法近似求解偏微分方程.

2） 使用同樣的近似方法（有限差分法）求原始偏微分方程的近似值.

在用有限差分方法對(duì)二階偏微分方程（5）—（6）做近似計(jì)算后，可以得到有限差分線性方程組，并求出整個(gè)模型中各個(gè)節(jié)點(diǎn)上的Ex和Hx的近似值.進(jìn)而可以求出磁場(chǎng)分量Hy和Hz.

求解過(guò)程為：

1） 將場(chǎng)分量Ex在中心節(jié)點(diǎn)的兩個(gè)負(fù)方向上，分別在各自的網(wǎng)格步長(zhǎng)中用泰勒公式展開到二階.

2） 用二階導(dǎo)數(shù)替換這些展開項(xiàng)，它們等于原來(lái)的二階方程.

3） 在中心節(jié)點(diǎn)的兩邊求出合適的電導(dǎo)率平均值.

4） 通過(guò)從先前步驟中減去適當(dāng)比例方程，消除Ex剩下的二階導(dǎo)數(shù).

在電性各向異性情況中，磁場(chǎng)分量Hx的影響也同樣必須重新考慮.它會(huì)從本質(zhì)上影響計(jì)算電導(dǎo)率平均值的過(guò)程.

在積分單元上對(duì)電導(dǎo)率或電阻率求平均值是可以重復(fù)多次進(jìn)行的，但之間的差別很微小，在電性各向異性結(jié)構(gòu)中，電流密度的各個(gè)分量不是對(duì)應(yīng)的電場(chǎng)分量乘以標(biāo)量電導(dǎo)率，而是所有電場(chǎng)分量乘以一個(gè)電導(dǎo)率張量.在分界面的邊界條件，須對(duì)各個(gè)區(qū)域的不同電導(dǎo)率張量進(jìn)行合適的修改.

下面將表明通過(guò)使用有限體積近似法，Weaver等（1986）的求導(dǎo)公式可用于電性各向異性情況.例如：Ex在z方向上的導(dǎo)數(shù)?Ex（yj,zk）/?z.對(duì)公式

圖2 變量在網(wǎng)格節(jié)點(diǎn)中的交錯(cuò)排列方式

（5）用體積積分法分別計(jì)算積分單元Gj,k的兩個(gè)子單元，即：積分區(qū)域的上半部分（z>zk）和下半部分（z

（9）

式中，

（10）

公式（10）前兩項(xiàng)只包含電場(chǎng)分量Ex，這與Weaver等（1986）用于計(jì)算電性各向同性情況下改進(jìn)的求導(dǎo)公式一致.公式（10）中的第一項(xiàng)代表一個(gè)導(dǎo)數(shù)值，它由網(wǎng)格節(jié)點(diǎn)（j,k-1）,（j,k）,（j,k+1）對(duì)Ex做拋物差值計(jì)算求出.第二項(xiàng)是一個(gè)修改量，它是原始偏微分方程的結(jié)果，代表了二階偏導(dǎo)數(shù)?2Ex/?z2的值.這個(gè)修改量與中心節(jié)點(diǎn)上下平均電導(dǎo)率的差成正比.對(duì)于存在較大電導(dǎo)率差值時(shí)，忽視這個(gè)修改量會(huì)造成導(dǎo)出的場(chǎng)值精度快速下降，即使場(chǎng)分量Ex非常準(zhǔn)確并有著很高的精度.

公式（10）的后兩項(xiàng)包含了磁場(chǎng)分量Hx，表示了在電性各向異性情況下，H極化模式的影響.這些項(xiàng)由電導(dǎo)率張量的非對(duì)角元素由集合參數(shù)A和B及它們的平均值B-和B+所決定.

同樣的，重復(fù)以上過(guò)程，可以容易求得積分區(qū)域Gj,k的左右積分子區(qū)域IE（Gj,k,yj-）和IE（Gj,k,yj+），并得到?Ex（j,k）/?y的近似值.

隨后可以計(jì)算由麥克斯韋方程組導(dǎo)出的電場(chǎng)分量Ey和Ez，將有限體積積分用于公式（6）中，求出?Hx/?y和?Hx/?z.

（11）

3 帶地形的二維電性各向異性模型計(jì)算與分析

在數(shù)據(jù)采集過(guò)程中，數(shù)據(jù)質(zhì)量是基礎(chǔ)是關(guān)鍵，大地電磁方法應(yīng)用成功的與否在很大程度上取決于數(shù)據(jù)采集質(zhì)量.我國(guó)地域遼闊，地形復(fù)雜多變，工業(yè)生產(chǎn)、人文因素引起的強(qiáng)干擾區(qū)較為普遍.

大地電磁法多數(shù)是研究大尺度構(gòu)造，雖然在多數(shù)情況下，地形的起伏相對(duì)于測(cè)線長(zhǎng)度變化較小，但也存在一些情況，例如在構(gòu)造帶附近，地形會(huì)有較為劇烈的變化，此時(shí)就要求依據(jù)實(shí)際情況建立含有起伏地形的地電模型用于處理解釋.對(duì)含地形因素的二維電性各向異性地電模型做正演研究，可以為分析研究存在二維電性各向異性介質(zhì)情況下，地形因素對(duì)大地電磁場(chǎng)產(chǎn)生的影響提供理論支撐.

為說(shuō)明地形的影響，首先對(duì)均勻半空間包含地壘和地塹兩種地形的情況分別做正演計(jì)算.均勻半空間ρ=100Ωm（圖3所示）.

圖3 含地形的均勻半空間

正演計(jì)算時(shí)將模型剖分為69×47的剖分單元（垂向47個(gè)剖分單元），測(cè)點(diǎn)數(shù)為69個(gè)（即在地表的網(wǎng)格節(jié)點(diǎn)處），采用21個(gè)周期點(diǎn)（0.1, 0.3, 0.5, 0.7, 1, 3, 5, 7, 10, 30,50,70,100,300,500,700,1000,1300,1500,1700,2000）.計(jì)算結(jié)果如圖4、圖5所示.

圖4、圖5分別是均勻半空間中地壘地形和地塹地形的大地電磁響應(yīng)結(jié)果.可以直觀看出：無(wú)論是地壘還是地塹，地形的影響主要表現(xiàn)在地形起伏區(qū)域，突出的畸變出現(xiàn)在地形劇烈變化的交匯處（平地與山坡交界及山坡與山頂交界處）.在地壘地形的起伏區(qū)域，TE模式的視電阻率值和阻抗相位值呈現(xiàn)出小-大-小的變化，TM模式的視電阻率值和阻抗相位值呈現(xiàn)出大-小-大的變化.而在地塹地形的起伏區(qū)域，TE模式和TM模式中視電阻率值和阻抗相位值的變化正好與地壘情況相反.此外，無(wú)論是地壘還是地塹，地形的影響對(duì)于視電阻率值來(lái)說(shuō)主要集中在TE模式的中高頻段和TM模式的低頻段，而對(duì)于阻抗相位來(lái)說(shuō)影響幾乎都是在中高頻段.最后，地形的影響程度和地形傾角有關(guān).傾角越大，影響就越大；反之，則越小.

圖4 含地壘地形的均勻半空間的大地電磁響應(yīng)擬斷面圖

圖5 含地塹地形的均勻半空間的大地電磁響應(yīng)擬斷面圖

網(wǎng)格剖分情況（圖7）及周期點(diǎn)與均勻半空間地壘相同，計(jì)算結(jié)果如圖8所示.

圖7 網(wǎng)格剖分圖

將圖8的結(jié)果同之前（胡祥云等，2013）未含地形的正演結(jié)果對(duì)比可知：含地壘情況下TE模式的視電阻率擬斷面圖中低阻異常區(qū)域在水平方向上的延展有所擴(kuò)大，圖8a中低阻區(qū)域展示了二維電性各向異性介質(zhì)的存在范圍，受地形影響在33～35km，lg（f）=0的低阻區(qū)域有明顯凹陷；阻抗相位擬斷面圖的形態(tài)基本一致，細(xì)微區(qū)別在于含地壘情況下阻抗相位值最大值和最小值區(qū)域水平方向上的延伸稍有增加，在33～35km，lg（f）=-1.5的阻抗相位值有所增大；受地形影響TE模式下的視電阻率最小值有所增大，阻抗相位值幾乎無(wú)變化.TM模式的視電阻率值和阻抗相位值受地形的影響較為明顯，在山腳處（24km和44km處）視電阻率值和阻抗相位值均發(fā)生了畸變，表現(xiàn)為視電阻率值跳躍增大，阻抗相位值跳躍降低；在山頂處（33～35km）視電阻率值有明顯降低，阻抗相位值有明顯增大；受地形影響TM模式下的視電阻率最小值有所降低，阻抗相位最大值有所增大.

在含地塹地形的層狀介質(zhì)中放置一個(gè)二維電性各向異性介質(zhì)，對(duì)地電模型做正演計(jì)算，并與之前（胡祥云等，2013）未含地形的正演結(jié)果做對(duì)比分析，介質(zhì)1、2、3參數(shù)同地壘地電模型.地塹的頂部寬度為20km（從24km到44km），底部為2km，海拔高度-3km.如圖9所示.

網(wǎng)格剖分情況及周期點(diǎn)與均勻半空間地塹情況相同，計(jì)算結(jié)果如圖10所示.

圖8 含地壘地形的層狀電性各向異性地電模型的大地電磁響應(yīng)擬斷面圖

圖9 含地塹地形的層狀電性各向異性地電模型

將圖10的結(jié)果同之前未含地形的正演結(jié)果（胡祥云等，2013）對(duì)比可知：含地塹情況下TE模式的視電阻率擬斷面圖中低阻異常區(qū)域在水平方向上的延展有所擴(kuò)大，圖1a中低阻區(qū)域展示了二維電性各向異性介質(zhì)的存在范圍，受地形影響在33～35km，lg（f）=0的低阻區(qū)域有明顯凸起；阻抗相位擬斷面圖的形態(tài)基本一致，細(xì)微區(qū)別在于含地塹情況下阻抗相位值最大值和最小值區(qū)域水平方向上的延伸稍有增加，在33～35km，lg（f）=-1位置阻抗相位值有所減小.TM模式的視電阻率值和阻抗相位值受地形的影響較為明顯，在山腳處（24km和44km處）視電阻率值和阻抗相位值均發(fā)生了畸變，表現(xiàn)為視電阻率值跳躍減小，阻抗相位值跳躍增大；在山頂處（33～35km）視電阻率值表現(xiàn)為增大-減小-增大的跳躍式畸變，阻抗相位值表現(xiàn)為減小-增大-減小的跳躍式畸變.

總體來(lái)說(shuō)，地形因素對(duì)于大地電磁響應(yīng)的影響主要存在于地形的突變位置（平地與山腳交匯處及山坡與山頂交匯處），而在平穩(wěn)過(guò)渡位置（山坡）對(duì)于大地電磁響應(yīng)的影響較為平緩.其次，不論是地壘地形還是地塹地形，對(duì)于TE模式下的視電阻率值和阻抗相位值影響都較??；對(duì)于TM模式的影響則較大.此外，地塹地形的影響要強(qiáng)于地壘地形的影響，最后，地形的影響程度和地形傾角有直接關(guān)系，地形傾角越大，影響越強(qiáng)烈；反之則越小.

4 大地電磁實(shí)測(cè)資料解釋及電性結(jié)構(gòu)對(duì)比分析

選取的大地電磁測(cè)線位于貴州某地，測(cè)線位于貴州山區(qū)，地形起伏較大，加入地形的正演擬合結(jié)果展示了地形對(duì)大地電磁響應(yīng)的影響.

圖11是測(cè)線位置圖，測(cè)線由西南向東北首先穿過(guò)山谷，隨后翻過(guò)山脊進(jìn)入山坡平緩區(qū)域.測(cè)線長(zhǎng)測(cè)線的測(cè)點(diǎn)距為40m，對(duì)于長(zhǎng)度為1400m的測(cè)線來(lái)說(shuō)，數(shù)據(jù)量是十分充足的.圖12中的4幅擬斷面圖都不同程度的出現(xiàn)了許多閉合圓圈，即使通過(guò)數(shù)據(jù)插值成圖，也不能使這些區(qū)域的高/低值平滑.說(shuō)明實(shí)測(cè)資料中在同一周期點(diǎn)，水平方向上的數(shù)據(jù)值跳動(dòng)十分劇烈.本文通過(guò)構(gòu)建一個(gè)含地形的二維電性各向異性地電模型，并對(duì)其做正演計(jì)算，根據(jù)正演擬合的結(jié)果分析該區(qū)域的電性結(jié)構(gòu)并研究地形對(duì)大地電磁響應(yīng)的影響.

圖10 含地塹地形的層狀電性各向異性地電模型的大地電磁響應(yīng)擬斷面圖

圖11 測(cè)線位置圖

1400m，地形起伏高差240m，測(cè)點(diǎn)距40m.測(cè)線的視電阻率擬斷面圖和阻抗相位擬斷面圖如圖12所示.

圖14中TE和TM模式的正演結(jié)果受地形的影響都較大.正演計(jì)算的結(jié)果和實(shí)測(cè)數(shù)據(jù)有很好的擬合.圖14a中，在lg（f）=1Hz，水平方向的200～400m，800～1400m區(qū)域出現(xiàn)了對(duì)應(yīng)的高阻區(qū)，在X=1000m和1200m出現(xiàn)的高阻閉合圓圈也有很好的對(duì)應(yīng)，并且高阻最大值也保持一致.圖14c和圖12c的形態(tài)基本保持一致，在阻抗相位的低值區(qū)域（lg（f）=2-1Hz，X=100～300m及800～1400m）都表現(xiàn)出高阻介質(zhì)（5、6、7）在深度方向上的延展范圍.圖14b中，在lg（f）=0.5Hz，水平方向的100～300m，1000～1200m區(qū)域，出現(xiàn)高阻閉合圓圈，與圖12b中的情況一致，高阻最大值也保持一致.圖14d和圖12d形態(tài)稍有差異，主要表現(xiàn)在水平方向600～800m范圍內(nèi).圖14d的低值區(qū)域在水平方向上很好的呈現(xiàn)出5、6、7三個(gè)電性各向異性介質(zhì)的延展范圍.根據(jù)正演擬合的結(jié)果，推斷地形的突變點(diǎn)（X=200、700、1200m）及5、6、7三個(gè)介質(zhì)的變化點(diǎn)（X=240、320、700、900、1140）是造成數(shù)據(jù)在同一周期點(diǎn)呈鋸齒狀起伏的主要原因.總體來(lái)說(shuō)，正演擬合的結(jié)果基本上與實(shí)測(cè)數(shù)據(jù)一致，不但在位置上一一對(duì)應(yīng)出高阻閉合圈的范圍，而且在高阻區(qū)域的最大值也能夠保持一致.

圖12 帶地形的二維大地電磁響應(yīng)擬斷面圖

圖13 帶地形的二維電性各向異性地電模型

實(shí)測(cè)資料的正演擬合結(jié)果表明：在大尺度下地形起伏較緩時(shí)可以忽略地形因素對(duì)大地電磁響應(yīng)造成的影響，而在地形起伏較大區(qū)域（如山區(qū)和構(gòu)造帶區(qū)域），地形變化會(huì)對(duì)大地電磁響應(yīng)造成顯著的影響，使得視電阻率值及阻抗相位值在地形突變位置產(chǎn)生畸變，此時(shí)就必須重視地形變化帶來(lái)的影響，必要時(shí)需對(duì)地形造成的畸變進(jìn)行校正.

5 結(jié)論

本文基于Pek（1997,2000）的研究理論，系統(tǒng)而詳實(shí)的闡述了含地形的二維電性各向異性正演理論；并通過(guò)理論研究說(shuō)明地形因素對(duì)大地電磁響應(yīng)的影響，實(shí)測(cè)資料的對(duì)比研究說(shuō)明地形因素對(duì)大地電磁響應(yīng)的影響程度.

基于水平地形二維大地電磁各向異性研究成果，引入地形因素，改變變量在網(wǎng)格節(jié)點(diǎn)中的排列方式，選擇交錯(cuò)排列方式從而給有限差分系數(shù)矩陣的最大帶寬分配合理的存儲(chǔ)空間，隨后，使用Weaver等（1985,1986）的方法解決TM模式下，在地-空分界面Y方向上的一些區(qū)域存在不同電導(dǎo)率的問(wèn)題，對(duì)帶地形的二維大地電磁電性各向異性正演問(wèn)題進(jìn)行系統(tǒng)研究.正演計(jì)算地壘及地塹情況下的地電模型，研究地形對(duì)大地電磁響應(yīng)的影響.

圖14 帶地形的二維電性各向異性地電模型的大地電磁響應(yīng)擬斷面圖

以本文的研究成果為基礎(chǔ)，將帶地形的二維大地電磁電性各向異性理論引入對(duì)實(shí)測(cè)大地電磁資料的處理解釋中，說(shuō)明電性各向異性現(xiàn)象的普遍存在，同時(shí)驗(yàn)證理論的正確性及程序的實(shí)用性，并且對(duì)山區(qū)地形起伏劇烈情況下地形變化對(duì)大地電磁響應(yīng)的影響程度深入分析，為今后處理解釋大地電磁資料中包含地形因素的電性各向異性現(xiàn)象提供理論依據(jù)和技術(shù)指導(dǎo)，并開拓了對(duì)大地電磁實(shí)測(cè)資料處理的思路和方法.

致謝 感謝捷克科學(xué)院地球物理研究所JosefPek教授對(duì)本文的幫助.

Brewitt-Taylor C R, Weaver J T. 1976. On the finite difference solution of two-dimensional induction problems.Geophys.J.Int., 47（2）: 375-396.

Cagniard L. 1953. Basic theory of the magneto-telluric method of geophysical prospecting.Geophysics, 18（3）: 605-635.

Chouteau M, Bouehard K. 1988. Two-dimensional terrain correction in magnetotelluric surveys.Geophysics, 53（6）: 854-862.

Haak V. 1972. Magnetotellurics: Determination of the transfer functions in regions with lateral changes of the electrical conductivity.Z.Geophys. （in German）, 38: 85-102.

Hill D G. 1972. A laboratory investigation of electrical anisotropy in Precambrian rocks.Geophysics, 37（6）: 1022-1038.

Hu X Y, Huo G P, Gao R, et al. 2013. The magnetotelluric anisotropic two-dimensional simulation and case analysis.ChineseJ.Geophys. （in Chinese）, 56（12）: 4268-4277, doi: 10.6038/cjg20131229.

Huo G P, Hu X Y, Fang H, et al. 2012. Magnetotelluric joint inversion for anisotropic conductivities in layered media.ActaPhys.Sin. （in Chinese）, 62（12）: 129101.

Huo G P, Hu X Y, Liu M. 2011. Review of the forward modeling of magnetotelluric in the anisotropy medium research.ProgressinGeophys. （in Chinese）, 26（6）: 1976-1982,doi: 10.3969/j.issn.1004-2903.2011.06.011.

Jin G W, Zhao G Z, Xu C F, et al. 1998. The affection and correction on magnetotelluric response data for inclination two-dimension terrain.SeismologyandGeology（in Chinese）, 20（4）: 454-458.

Keller G V. 1968. Electrical prospecting for oil.∥ Quarterly of the Colorado School of Mines, v. 63, no. 2. Golden: Colorado School of Mines.Key K, Ovall J. 2011. A parallel goal-oriented adaptive finite element method for 2.5-D electromagnetic modelling.Geophys.J.Int., 186（1）: 137-154.

Ku C C, Hsieh M S, Lim S H. 1973. The topographic effect in electromagnetic fields.CanadianJ.EarthSci., 10（5）: 645-656.

Li Y G. 2002. A finite-element algorithm for electromagnetic induction in two-dimensional anisotropic conductivity structures.Geophys.J.Int., 148（3）: 389-401.

Li Y G, Pek J. 2008. Adaptive finite element modelling of two-dimensional magnetotelluric fields in general anisotropic media.Geophys.J.Int., 175（3）: 942-954.

O’Brien D P, Morrison H F. 1967. Electromagnetic fields in an N-layer anisotropic half-space.Geophysics, 32（4）: 668-677. Parkhomenko E I. 1967. Electrical Properties of Rocks. US: Springer.

Pek J, Verner T. 1997. Finite-difference modelling of magnetotelluric fields in two-dimensional anisotropic media.Geophys.J.Int., 128（3）: 505-521.

Pek J, Toh H. 2000. Numerical modeling of MT fields in 2-D anisotropic structures with topography and bathymetry considered.∥ Protokolluber das 18, Kolloquium “Electromagnetische Tiefenforschung”. Altenberg, Germany, Hordt, A., and J. Stoll, eds., DGG, 190-199.

Pek J, Santos F A M. 2002. Magnetotelluric impedances and parametric sensitivities for 1-D anisotropic layered media.ComputersandGeosciences, 28（8）: 939-950.

Reddy I K, Rankin D. 1971. Magnetotelluric effect of dipping anisotropies.GeophysicalProspecting, 19（1）: 84-97.

Reddy I K, Rankin D. 1975. Magnetotelluric response of laterally inhomogeneous and anisotropic media.Geophysics, 40（6）: 1035-1045.

Tikhonov A N. 1950. On determining electrical characteristics of the deep layers of the earth′s crust.Deki.Akud.Nuck., 73: 295-297. Wang X B, Li Y N, Gao Y C. 1999. Two dimensional topographic responses in magneto telluric sounding and it′s correction methods.ComputingTechniquesforGeophysicalandGeochemicalExploration（in Chinese）, 21（4）: 327-332.

Wannamaker P E, Stodt J A, Rijo L. 1986. Two-dimensional topographic responses in magnetotellurics modeled using finite elements.Geophysics, 51（11）: 2131-2144.

Wannamaker P E. 2005. Anisotropy versus heterogeneity in continental solid Earth electromagnetic studies: fundamental response characteristics and implications for physicochemical state.Surv.Geophys., 26（6）: 733-765.

Weaver J T, Le Quang B V, Fischer G. 1985. A comparison of analytic and numerical results for a two-dimensional control model in electromagnetic induction-I. B-polarization calculations.Gephys.J.Int., 82（2）: 263-277.

Weaver J T, Le Quang B V, Fischer G. 1986. A comparison of analytical and numerical results for a 2-D control model in .

J.Int., 87（3）: 917-948.

Wei W B, Jin S, Ye G F, et al. 2006. Conductivity structure of crust and upper mantle beneath the northern Tibetan Plateau: Results of super-wide band magnetotelluric sounding.ChineseJ.Geophys. （in Chinese）, 49（4）: 1215-1225.

Wei W B, Jin S, Ye G F, et al. 2010. On the conductive structure of Chinese continental lithosphere-experiment on “standard monitoring network” of continental EM parameters （SinoProbe-01）.ActaGeologicaSinica（in Chinese）, 84（6）: 788-800.

Weiss C J, Newman G A. 2002. Electromagnetic induction in a fully 3-D anisotropic earth.Geophysics, 67（4）: 1104-1114.

Xu S Z, Zhao S K. 1985. The topographic effects on magnetotelluric response.NorthwesternSeismologicalJournal（in Chinese）, 7（4）: 69-78.

Xu S Z, Wang Q Y, Wang J. 1992. Modeling 2-D terrain effect on MT by the boundary element method.ActaGeophysicaSinica（in Chinese）, 35（3）: 380-388.

Xu S Z, Li Y G, Liu B. 1997. The method and efficiency of 2-D terrain correction for Hx-polarization of MT.ChineseJ.Geophys. （in Chinese）, 40（6）: 842-846.

Zhao G Z, Tang J, Zhan Y, et al. 2004. Study on the crustal electrical structure and block deformation in Northeast margin of Tibetan Plateau.ScienceinChina（Ser.D） （in Chinese）, 34（10）: 908-918.

附中文參考文獻(xiàn)

胡祥云, 霍光譜, 高銳等. 2013. 大地電磁各向異性二維模擬及實(shí)例分析. 地球物理學(xué)報(bào), 56（12）: 4268-4277, doi: 10.6038/cjg20131229.

霍光譜, 胡祥云, 方慧等. 2012. 層狀各向異性介質(zhì)大地電磁聯(lián)合反演研究. 物理學(xué)報(bào), 61（12）: 129101.

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

魏文博, 金勝, 葉高峰等. 2006. 藏北高原地殼及上地幔導(dǎo)電性結(jié)構(gòu)-超寬頻帶大地電磁測(cè)深研究結(jié)果. 地球物理學(xué)報(bào), 49（4）: 1215-1225.

魏文博, 金勝, 葉高峰等. 2010. 中國(guó)大陸巖石圈導(dǎo)電性結(jié)構(gòu)研究-大陸電磁參數(shù)“標(biāo)準(zhǔn)網(wǎng)”試驗(yàn)（SinoProbe-01）. 地質(zhì)學(xué)報(bào), 84（6）: 788-800.

徐世浙, 趙生凱. 1985. 地形對(duì)大地電磁勘探的影響. 西北地震學(xué)報(bào), 7（4）: 69-78.

徐世浙, 王慶乙, 王軍. 1992. 用邊界單元法模擬二維地形對(duì)大地電磁場(chǎng)的影響. 地球物理學(xué)報(bào), 35（3）: 380-388.

徐世浙, 李予國(guó), 劉斌. 1997. 大地電磁Hx型波二維地形改正的方法與效果. 地球物理學(xué)報(bào), 40（6）: 842-846.

趙國(guó)澤, 湯吉, 詹艷等. 2004. 青藏高原東北緣地殼電性結(jié)構(gòu)和地塊變形關(guān)系的研究. 中國(guó)科學(xué)（D輯）, 34（10）: 908-918.

（本文編輯 汪海英）

MT modeling for two-dimensional anisotropic conductivity structure with topography and examples of comparative analyses

HUO Guang-Pu1, HU Xiang-Yun2*, HUANG Yi-Fan2, HAN Bo2

1HenanInstituteofGeologicalSurvey,Zhengzhou450000,China2HubeiSubsurfaceMulti-scaleImagingKeyLaboratory,InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan430074,China

The magnetotelluric （MT） method is a technique for probing electrical conductivity structure of the Earth, from the near-surface down to upper mantle. MT observations can be significantly influenced by topography. Most of existing forward modeling algorithms for 2-D MT with topography are based on the electrical isotropic theory. However, it has been well established that the electrical anisotropy is widely present in the earth interior. Since the electrical anisotropy in crust and upper mantle is the connection between electric structure and the geological structure, it is vital to account for anisotropic effects while modeling 2-D MT fields with topography.We present a 2-D MT modeling approach for anisotropic media with topography. The solution is based on a finite-difference （FD） discretization of the second-order Maxwell′s equation in terms of electric fields parallel to strike for TE mode and magnetic fields parallel to strike for TM mode. The topography effect is accounted by changing the way in which the sampling electromagnetic fields are arranged, i.e. the sampling fields are in a staggered-order to make sure that the maximum bandwidth of the FD coefficient matrix can be assigned moderate memory spaces. The Weaver's approach is used to deal with the problem that the conductivity of some area on the air-earth interface may vary in the direction perpendicular to strike. Once the primary field is solved, the dual field can be obtained very easily by applying a discrete differential operator to the primary field.Two synthetic models for modeling the horst structure and graben structure, respectively, are used to evaluate the effects of topography. The responses of models with topography are compared with that of models without topography. It is found that most significant differences occur in the regions with sharp topography boundaries, such as the boundary between the foot of a mountain and the mountainslope and the boundary between the mountain slope and the mountaintop, while the minor differences appear in flat regions. The topography has much greater impact on the TM responses than TE responses, either for the horst model or the graben model. The graben structure can have more effects than the horst structure. The sharper the slope, the greater the influence is.We demonstrate the validity of the algorithm for 2-D MT anisotropic forward modeling with topography by numerical tests on both synthetic models and real datasets. The effects of topography are assessed by analyzing the modeling results of the horst model and the graben model. Finally, the impacts of the topography on MT responses for both large scale models and sharp slope models are evaluated by fitting the observed MT data through forward modeling.

Topography； Electrical anisotropy; Magnetotelluric; Finite difference; Tensor conductivity

國(guó)家深部探測(cè)專項(xiàng)第3項(xiàng)目（SinoProbe-03）,“十二五”國(guó)家科技支撐計(jì)劃課題（2011BAB04B01），國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃（2013CB733203）和國(guó)家自然科學(xué)基金（41274077、41474055）聯(lián)合資助.

霍光譜,男,1983年生,博士,主要從事電磁法正反演研究.E-mail：huoguangpu1218@163.com

*通訊作者 胡祥云,男,1966年生,教授,博士生導(dǎo)師,主要從事電磁法理論應(yīng)用教學(xué)與研究工作.E-mail:xyhu@cug.edu.cn

10.6038/cjg20151230.

10.6038/cjg20151230

P631

2015-05-15，2015-10-30收修定稿

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

Huo G P, Hu X Y, Huang Y F, et al. 2015. MT modeling for two-dimensional anisotropic conductivity structure with topography and examples of comparative analyses.ChineseJ.Geophys. （in Chinese）,58（12）:4696-4708,doi:10.6038/cjg20151230.