基于構(gòu)造系統(tǒng)函數(shù)的大地電磁時(shí)間序列模擬方法

張寶強(qiáng)　裴建新*②　王　啟

(①中國海洋大學(xué)海洋地球科學(xué)學(xué)院，山東青島 266100； ②海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東青島 266100； ③廣州海洋地質(zhì)調(diào)查局，廣東廣州 510760)

1　引言

大地電磁測(cè)深(MT)法是一種在地質(zhì)探測(cè)和礦產(chǎn)普查中用途非常廣泛的地球物理方法，通過分析視電阻率曲線和相位曲線的特征獲取地下結(jié)構(gòu)電性分布。MT觀測(cè)數(shù)據(jù)的阻抗計(jì)算以及對(duì)阻抗特性產(chǎn)生嚴(yán)重影響的噪聲處理是研究的重要內(nèi)容之一[1-4]。近些年，國內(nèi)外學(xué)者在MT阻抗估算和噪聲處理研究上做出了很多努力，頻率域模擬電磁場(chǎng)數(shù)據(jù)已用于驗(yàn)證和評(píng)估MT阻抗估計(jì)方法[5-7]和提取被污染信號(hào)的預(yù)處理技術(shù)[8,9]。與此同時(shí)，時(shí)間域數(shù)據(jù)處理方法的發(fā)展[10-13]也推動(dòng)了MT數(shù)據(jù)處理方法的進(jìn)步。如今，國內(nèi)外學(xué)者使用不同的方法模擬MT時(shí)間序列； Larsen等[6]使用實(shí)測(cè)磁場(chǎng)時(shí)間序列、通過失真函數(shù)修改估計(jì)的一維傳遞函數(shù)，在測(cè)試區(qū)域中生成電場(chǎng)時(shí)間序列； Chen[14]在測(cè)試經(jīng)驗(yàn)?zāi)B(tài)分解算法時(shí)，采用的模擬數(shù)據(jù)是由實(shí)測(cè)磁場(chǎng)和一維地球模型的阻抗合成的頻率域電場(chǎng)數(shù)據(jù)。鑒于實(shí)測(cè)MT數(shù)據(jù)均含有噪聲，且通常地下電性結(jié)構(gòu)是未知的，所以實(shí)測(cè)數(shù)據(jù)在測(cè)試時(shí)間域阻抗估計(jì)方法和信號(hào)處理技術(shù)可靠性方面具有一定的局限性。Lamarque[9]和Goubau等[5]分別使用簡單的高斯隨機(jī)序列和均勻分布在(-1,1)的隨機(jī)序列作為無噪聲的MT時(shí)間序列及磁場(chǎng)分量的實(shí)部和虛部。在此基礎(chǔ)上，有學(xué)者利用二維阻抗張量生成了電場(chǎng)數(shù)據(jù)； Yee[7]和Filipo等[8]利用數(shù)字遞歸濾波器設(shè)計(jì)了隨頻率增大而衰減的指數(shù)函數(shù)描繪天然磁場(chǎng)的主要頻譜特征，將濾波器系數(shù)和偽隨機(jī)序列進(jìn)行卷積計(jì)算獲得磁場(chǎng)時(shí)間序列，并用模擬的磁場(chǎng)時(shí)間序列和模型阻抗的脈沖響應(yīng)進(jìn)行卷積計(jì)算獲得無噪的電場(chǎng)時(shí)間序列；Loddo等[15]提出了一種合成寬頻帶電磁場(chǎng)時(shí)間序列的方法，研究了場(chǎng)的極化特征，其思想是將構(gòu)造的系統(tǒng)函數(shù)作為頻率域磁場(chǎng)，建立層狀模型阻抗的頻率域積分解析式，使用傅里葉積分變換算法獲得磁場(chǎng)和阻抗的時(shí)間域脈沖響應(yīng)，磁場(chǎng)脈沖響應(yīng)與一個(gè)隨機(jī)數(shù)序列做褶積運(yùn)算生成合成磁場(chǎng)時(shí)間序列，然后與阻抗脈沖響應(yīng)做褶積運(yùn)算獲得合成電場(chǎng)時(shí)間序列； 俄羅斯學(xué)者Varentsov等[16]使用合成數(shù)據(jù)庫(COMDAT Project)比較不同的MT數(shù)據(jù)處理方法，其中，由選定的隨機(jī)函數(shù)乘以簡單的場(chǎng)得到磁場(chǎng)頻譜，并結(jié)合任意一維標(biāo)量阻抗計(jì)算電場(chǎng)頻譜； 嚴(yán)家斌[17]、蔡劍華等[18]、陳知富等[19]在研究中使用隨機(jī)數(shù)作為MT各分量的時(shí)間序列；凌振寶等[20]將多個(gè)單頻率的正弦波線性疊加模擬電磁場(chǎng)數(shù)據(jù)； 湯井田等[2]使用實(shí)測(cè)磁場(chǎng)數(shù)據(jù)和簡單二維電阻率模型阻抗計(jì)算合成電場(chǎng)數(shù)據(jù)。

本文提出一種MT時(shí)間序列模擬方法，不僅可以避免直接使用實(shí)測(cè)數(shù)據(jù)作為合成數(shù)據(jù)而引入干擾信息的不確定性，還可以改善僅使用隨機(jī)數(shù)序列及正弦波疊加序列作為合成數(shù)據(jù)而難以體現(xiàn)MT基本頻譜特性的問題，該方法支持不同的采樣率，得到的模擬數(shù)據(jù)可用于進(jìn)一步的數(shù)據(jù)處理方法研究。

2　方法原理

Loddo等[15]用構(gòu)造的系統(tǒng)函數(shù)作為磁場(chǎng)單邊譜刻畫MT頻譜。在此基礎(chǔ)上，本文提出了基于構(gòu)造系統(tǒng)函數(shù)、通過模型阻抗實(shí)現(xiàn)包含地下地電結(jié)構(gòu)信息的電磁場(chǎng)數(shù)據(jù)的模擬方法。

2.1　電磁場(chǎng)頻譜的求取

在確定模擬時(shí)間序列的采樣頻率Fs和序列長度N之后，輸入系統(tǒng)函數(shù)；為使模擬的電磁場(chǎng)頻譜和后續(xù)模擬的時(shí)間序列中帶有一定限度的波動(dòng)和跳躍、且隨機(jī)數(shù)的DFT序列的量級(jí)在整個(gè)頻段保持在1附近而不改變系統(tǒng)函數(shù)刻畫的磁場(chǎng)頻譜的量級(jí)[7,8,15,16,21,22]，需引入一均值為0、方差為1、長度為N的隨機(jī)序列，對(duì)其進(jìn)行歸一化FFT計(jì)算[22]，然后在頻率域與系統(tǒng)函數(shù)乘積，該計(jì)算結(jié)果就作為模擬的頻率域磁場(chǎng)。

為了使模擬數(shù)據(jù)包含地下的地電結(jié)構(gòu)信息，建立一維水平層狀電阻率模型，使用正頻率序列計(jì)算模型阻抗，利用頻率域電磁場(chǎng)和阻抗關(guān)系式計(jì)算模擬的頻率域電場(chǎng)。

2.2　電磁場(chǎng)時(shí)間序列的計(jì)算

傅里葉分析建立了函數(shù)時(shí)間域與頻率域之間的聯(lián)系，由模擬的頻率域磁場(chǎng)和電場(chǎng)分別做IFFT計(jì)算得到磁場(chǎng)和電場(chǎng)時(shí)間序列；根據(jù)實(shí)信號(hào)DFT的正負(fù)頻率復(fù)共軛對(duì)稱性質(zhì)[21-23]，在引入隨機(jī)數(shù)序列之前對(duì)系統(tǒng)函數(shù)做復(fù)共軛對(duì)稱處理，擴(kuò)展到負(fù)頻率使其成為雙邊對(duì)稱譜。復(fù)共軛對(duì)稱處理公式為

Re[f(-ω)]=Re[f(ω)]

Im[f(-ω)]=-Im[f(ω)]

(1)

式中f(ω)代表正頻率ω的頻域磁場(chǎng)h(ω)或阻抗序列z(ω)。同樣地，將正演的阻抗序列做復(fù)共軛對(duì)稱處理擴(kuò)展為雙邊對(duì)稱譜，最后將電場(chǎng)和磁場(chǎng)的雙邊對(duì)稱譜利用IFFT算法即可獲得模擬的實(shí)數(shù)信號(hào)。

基于以上思路，模擬MT時(shí)間序列的基本流程如圖1所示。

圖1　模擬MT時(shí)間序列的基本流程

3　一維磁場(chǎng)數(shù)據(jù)模擬

3.1　頻率域磁場(chǎng)數(shù)據(jù)模擬

模擬MT數(shù)據(jù)之前需要分析天然電磁場(chǎng)的基本頻譜特性。天然電磁場(chǎng)的產(chǎn)生有多種原因，高于6Hz的電磁場(chǎng)源主要來自于高能電磁現(xiàn)象，如閃電、磁暴等氣象活動(dòng)，低于6Hz的電磁場(chǎng)來源于地球磁場(chǎng)和地球電離層。不同頻帶的天然電磁場(chǎng)的基本特征已有大量學(xué)者做了詳細(xì)闡述[24-26]。磁場(chǎng)水平分量的平均振幅值通常隨著周期的增長而增大，在104～10-4s內(nèi)可以達(dá)到10nT； 此外，在1～6Hz及1k～2kHz之間分別存在著一個(gè)極小值，這兩個(gè)極小值代表了低頻和高頻天然電磁場(chǎng)場(chǎng)源機(jī)理的不同。圖2所示為一個(gè)通過線性時(shí)不變系統(tǒng)模擬的、經(jīng)平滑處理的水平磁場(chǎng)分量的平均頻譜，可描述天然磁場(chǎng)的主要頻譜特點(diǎn)。

圖2　線性時(shí)不變系統(tǒng)模擬的水平磁場(chǎng)頻譜[15]

一般地磁測(cè)量、大地電磁、音頻大地電磁方法僅研究有限頻帶范圍內(nèi)的電磁場(chǎng)，因此本研究也僅在有限的頻帶范圍內(nèi)模擬MT時(shí)間序列數(shù)據(jù)。Yee[7]和Filipo等[8]提出，可以直接以Z域中極點(diǎn)和零點(diǎn)的方式編寫系統(tǒng)的離散傳遞函數(shù)來近似天然場(chǎng)的頻率特性; Oppenheim等[27]認(rèn)為還有一種替代方案，即通過雙線性變換將傳遞函數(shù)從連續(xù)域變換到離散域。在模擬MT時(shí)間序列時(shí)，這些方法都是適用的。通過使用連續(xù)域中的極點(diǎn)和零點(diǎn)形式，一個(gè)系統(tǒng)函數(shù)可以表示為[28]

(2)

式中：s=σ+iω是拉普拉斯復(fù)數(shù)變量，σ是系統(tǒng)阻尼因子； 復(fù)數(shù)s01,s02,…,s0m分別代表系統(tǒng)函數(shù)的m個(gè)零點(diǎn)；sp1,sp2,…,spn分別代表系統(tǒng)函數(shù)的n個(gè)極點(diǎn)；K是比例系數(shù)。

函數(shù)的解析性質(zhì)完全由極點(diǎn)和零點(diǎn)在s復(fù)數(shù)平面上的分布決定。為了使系統(tǒng)穩(wěn)定，極點(diǎn)必須在s復(fù)平面的左半部分，并且任何落在復(fù)軸iω上的極點(diǎn)必須是單極點(diǎn)；此外，當(dāng)變量s虛部為0(即為實(shí)數(shù))時(shí)，對(duì)應(yīng)的函數(shù)值也是實(shí)數(shù)，要求極點(diǎn)和零點(diǎn)必須是實(shí)數(shù)或者復(fù)共軛對(duì)。單極點(diǎn)和單零點(diǎn)、復(fù)共軛對(duì)極點(diǎn)和復(fù)共軛對(duì)零點(diǎn)的分子和分母可以分解為[15]

H(s)=

(3)

式中：i、j、l、k分別表示單零點(diǎn)、單極點(diǎn)、復(fù)共軛零點(diǎn)和復(fù)共軛極點(diǎn)；K1為比例系數(shù)；ζ1和χ1代表每個(gè)臨界頻率的阻尼系數(shù)。

對(duì)于圖2所示的天然磁場(chǎng)的頻譜，坡度向下變化(能量降低)可以由式(3)的極點(diǎn)刻畫，而向上變化(能量升高)可以由式(3)的零點(diǎn)刻畫，這樣四個(gè)極點(diǎn)值和四個(gè)零點(diǎn)值就基本表達(dá)了天然磁場(chǎng)頻譜的總體變化趨勢(shì)，單極點(diǎn)和單零點(diǎn)都選擇實(shí)數(shù)，如表1所示。對(duì)于復(fù)共軛對(duì)，選擇系數(shù)ζ1=1、χ1=1，就相當(dāng)于二重極點(diǎn)和二重零點(diǎn)都選為實(shí)數(shù)。為了使系統(tǒng)穩(wěn)定，所有極點(diǎn)頻率都取負(fù)值[15]，為了使磁場(chǎng)頻譜曲線低頻處的漸近線為10nT，選擇比例系數(shù)K1=10。

表1　模擬磁場(chǎng)使用的極點(diǎn)和零點(diǎn)的頻率絕對(duì)值和特性

滿足以上條件后，式(3)轉(zhuǎn)化為

(4)

式中：zi(i=1,2,3,4)為零點(diǎn)；pi(i=1,2,3,4)為極點(diǎn)。

3.2　時(shí)間域磁場(chǎng)數(shù)據(jù)模擬

歸一化的DFT序列做乘積，作為模擬的頻率域磁場(chǎng)。這樣可以使得磁場(chǎng)曲線能夠平穩(wěn)變化，更接近實(shí)際情況。得到頻率域磁場(chǎng)數(shù)據(jù)之后，由IFFT便可獲得磁場(chǎng)時(shí)間序列。

為了驗(yàn)證本文方法，模擬一段采樣頻率為2Hz、采樣點(diǎn)數(shù)為10000的電磁場(chǎng)時(shí)間序列數(shù)據(jù)，其頻率范圍為0.0002～1Hz，參與計(jì)算的正頻率點(diǎn)數(shù)為5000個(gè)，磁場(chǎng)振幅譜由式(4)構(gòu)造，如圖3所示。圖4為模擬的磁場(chǎng)頻譜和時(shí)間序列，時(shí)間序列記錄長度為5000s，總體上模擬的磁場(chǎng)信號(hào)隨時(shí)間變化平穩(wěn)，具體數(shù)據(jù)統(tǒng)計(jì)參量：均值為-9.7868×10-4nT，均方差為1.5348×10-5nT，最小值為-9.7×10-3nT，最大值為1.05×10-2nT，能量為0.1630nT2。

圖3　系統(tǒng)函數(shù)構(gòu)造的磁場(chǎng)振幅譜

圖4　模擬的磁場(chǎng)頻譜和時(shí)間序列，頻率范圍為0.0002～1Hz，時(shí)間長度為5000s

4　一維電場(chǎng)數(shù)據(jù)模擬

4.1　一維層狀模型阻抗計(jì)算

電場(chǎng)可利用關(guān)系式E(ω)=Z(ω)·H(ω)計(jì)算，所以必須獲得地下介質(zhì)的阻抗信息。本文引入一維水平層狀模型，如圖5所示，設(shè)大地由n層水平層狀介質(zhì)所組成，各層的電阻率為ρ1,ρ2,…,ρn, 厚度為h1,h2,…,hn-1。可利用同一層中兩個(gè)不同深度處阻抗值之間的關(guān)系及界面上阻抗連續(xù)條件推導(dǎo)地面阻抗[29]。

圖5　n層介質(zhì)模型

由于層狀介質(zhì)中各層的電阻率ρ和傳播系數(shù)k不同，任一層介質(zhì)中的波動(dòng)方程為

(5)

假定第m層的頂面深度為zm，底面深度為zm+1，則厚度hm=zm+1-zm。第m層的特征阻抗為Z0m=-iωμ/km，于是可得阻抗的遞推公式

(6)

需注意，在計(jì)算模型阻抗的過程中，參與計(jì)算的頻率點(diǎn)與計(jì)算磁場(chǎng)單邊譜的頻率點(diǎn)需保持一致，頻率范圍為(0，F(xiàn)s/2)，為使模擬的電場(chǎng)頻譜為對(duì)稱譜，需將阻抗擴(kuò)展到負(fù)頻率部分使其對(duì)稱。圖6為設(shè)計(jì)的一個(gè)一維水平四層模型，使用相同的頻率點(diǎn)計(jì)算該模型的大地阻抗的實(shí)部和虛部(圖7)。

4.2　時(shí)間域電場(chǎng)數(shù)據(jù)模擬

獲得模擬磁場(chǎng)數(shù)據(jù)和模型的阻抗之后，在頻率域根據(jù)電磁場(chǎng)關(guān)系式E(ω)=Z(ω)·H(ω)計(jì)算電場(chǎng)的頻譜序列，對(duì)電場(chǎng)對(duì)稱頻譜做離散反傅里葉變換，就可得到模擬的電場(chǎng)時(shí)間序列。下列幾點(diǎn)需注意： ①模擬電場(chǎng)的零頻率值無法根據(jù)關(guān)系式得到，本文取模擬的磁場(chǎng)零頻率對(duì)應(yīng)的DFT值作為電場(chǎng)零頻率的DFT值； ②需考慮采樣長度N的情況，尤其是當(dāng)N為偶數(shù)時(shí)，奈奎斯特頻率處的DFT值必須是實(shí)數(shù)。圖8是由模擬磁場(chǎng)和阻抗計(jì)算的模擬電場(chǎng)頻譜和時(shí)間序列，其時(shí)間長度為5000s，頻率范圍為0.0002～1Hz。總體來看模擬的電場(chǎng)信號(hào)也是隨時(shí)間平穩(wěn)變化的，如圖8d所示，其數(shù)據(jù)統(tǒng)計(jì)參量： 均值為-1.2292V/m，均方差為7.2409×10-5V/m，最小值為-1.2504V/m，最大值為-1.2055V/m，能量為1.5111×104(V/m)2。

圖6　一維水平層狀模型

圖7　模型阻抗的實(shí)部和虛部

通過模擬的電磁場(chǎng)時(shí)間序列，可進(jìn)一步計(jì)算得到模型的視電阻率和相位曲線，結(jié)果如圖9所示。

4.3　模擬數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)對(duì)比

為了檢驗(yàn)?zāi)M效果，選取一段某海域?qū)崪y(cè)的海洋MT時(shí)間序列數(shù)據(jù)，為TE模式的水平電場(chǎng)x分量和磁場(chǎng)y分量，采樣頻率為10Hz，時(shí)間長度為15min，計(jì)算其振幅譜，然后用本文提出的方法模擬電磁場(chǎng)數(shù)據(jù)，如圖10所示。對(duì)比實(shí)測(cè)數(shù)據(jù)與模擬數(shù)據(jù)振幅譜可以看到，模擬的電磁場(chǎng)振幅譜基本體現(xiàn)了大地電磁場(chǎng)在低、中頻段隨著頻率的增大，能量逐漸衰減的總體趨勢(shì)，即信號(hào)頻譜呈下降趨勢(shì)，而中頻區(qū)段之后頻譜緩慢爬升的一般特征[29,30]。

圖8　模擬的電場(chǎng)頻譜和時(shí)間序列

圖9　模擬數(shù)據(jù)計(jì)算的視電阻率(左)和相位(右)曲線

圖10　實(shí)測(cè)海洋MT數(shù)據(jù)(左)與模擬MT數(shù)據(jù)(右)的時(shí)間序列和振幅譜

5　二維水平非均勻介質(zhì)電磁場(chǎng)數(shù)據(jù)合成

在實(shí)現(xiàn)一維MT時(shí)間序列模擬的基礎(chǔ)上，進(jìn)一步合成得到二維水平非均勻介質(zhì)電磁場(chǎng)數(shù)據(jù)。假設(shè)地下介質(zhì)是非均勻的，并且測(cè)量軸與介質(zhì)的電性主軸重合，則任何一個(gè)橢圓偏振波都可以看成是兩個(gè)線性偏振波的合成，這兩個(gè)線性偏振波就如同分別在兩個(gè)不同電阻率的均勻介質(zhì)中傳播。那么對(duì)應(yīng)的場(chǎng)可以分解為分別平行和垂直于構(gòu)造走向[29]的兩個(gè)部分。在這種情況下，每個(gè)電場(chǎng)分量和與之垂直的磁場(chǎng)分量有關(guān)，相應(yīng)的電磁場(chǎng)阻抗關(guān)系為

(7)

式中:ZTE=Zxy;ZTM=Zyx;Zxy≠Zyx。

湯井田等[2]和Loddo等[15]對(duì)給定的兩個(gè)地層模型進(jìn)行一維正演，得到兩組隨頻率變化的阻抗值，將兩個(gè)阻抗值分別作為一個(gè)矩陣的非對(duì)角元素，合成的矩陣即可用來模擬地下二維結(jié)構(gòu)。利用該方法，本文在一維模擬的基礎(chǔ)上通過分量合成的方法獲得二維水平非均勻介質(zhì)電磁場(chǎng)數(shù)據(jù)，建立兩個(gè)不同的一維水平層狀電阻率模型(表2)，將獲得的兩個(gè)阻抗響應(yīng)分別作為Zxy和Zyx。在模擬兩個(gè)磁場(chǎng)分量的時(shí)候，均使用式(4)構(gòu)造磁場(chǎng)單邊頻譜，但是兩個(gè)分量采用的隨機(jī)序列是不同的，以此來獲得兩個(gè)不同的磁場(chǎng)分量Hx和Hy，并由式(7)計(jì)算得到電場(chǎng)分量Ex和Ey。

表2　兩個(gè)一維水平層狀模型的電性參數(shù)

設(shè)定采樣頻率為10Hz，時(shí)間長度為1h，采樣點(diǎn)數(shù)為36000，計(jì)算得到該模型的電磁場(chǎng)水平分量如圖11所示。在上述水平層狀電阻率模型下模擬的水平電場(chǎng)和磁場(chǎng)是變化緩慢的平穩(wěn)信號(hào)。此外，相比磁場(chǎng)，電場(chǎng)變化相對(duì)更加平穩(wěn)，體現(xiàn)了大地電磁信號(hào)隨時(shí)間變化的基本特征，由模擬時(shí)間序列計(jì)算TE模式和TM模式的模型視電阻率和相位曲線(圖12)，可見視電阻率和相位曲線光滑連續(xù)，反映了地下模型介質(zhì)的電阻率和相位特征。由于模擬的每一個(gè)采樣點(diǎn)的電場(chǎng)和磁場(chǎng)都滿足阻抗關(guān)系式 ，所以模擬信號(hào)不包含任何噪聲，這是因?yàn)槿魏伟陔妶?chǎng)或磁場(chǎng)中的相關(guān)或不相關(guān)噪聲都不會(huì)滿足該阻抗關(guān)系式。

圖11　模擬二維MT水平分量時(shí)間序列

圖12　模擬數(shù)據(jù)計(jì)算的視電阻率(上)和相位(下)曲線

6　結(jié)束語

本文實(shí)現(xiàn)了一種模擬MT時(shí)間序列的方法，該方法不僅在模擬電磁場(chǎng)時(shí)結(jié)合了天然電磁場(chǎng)的頻譜特性，而且模擬的電磁場(chǎng)數(shù)據(jù)包含了地下介質(zhì)的視電阻率和相位信息，可以任意選取采樣頻率和序列長度，適用于一維電阻率模型和二維水平非均勻介質(zhì)的電阻率模型。由該方法獲得的模擬數(shù)據(jù)與實(shí)測(cè)的海洋MT數(shù)據(jù)比較發(fā)現(xiàn)，模擬的電磁場(chǎng)振幅譜基本體現(xiàn)了大地電磁場(chǎng)在低、中頻段隨著頻率的增加，能量衰減的總體趨勢(shì)。該方法模擬的的數(shù)據(jù)并未考慮噪聲的影響，因此下一步的研究應(yīng)包含噪聲和更復(fù)雜的電阻率模型，使其適用范圍更廣。

[1]張全勝,王家映.大地電磁測(cè)深資料的去噪方法.石油地球物理勘探,2004，39(增刊1):17-23.

Zhang Quansheng,Wang Jiaying.An method of noise elimination for magnetotelluric sounding data.OGP,2004，39(S1):17-23.

[2]湯井田,張弛,肖曉等.大地電磁阻抗估計(jì)方法對(duì)比.中國有色金屬學(xué)報(bào),2013，23(9):2351-2358.

Tang Jingtian,Zhang Chi,Xiao Xiao et al.Comparison of methods for magnetotelluric impedance estimation.The Chinese Journal of Nonferrous Metals,2013，23(9):2351-2358.

[3]蔡劍華,王先春,胡惟文.基于經(jīng)驗(yàn)?zāi)B(tài)分解與小波閾值的MT信號(hào)去噪方法.石油地球物理勘探,2013,48(2):303-307.

Cai Jianhua,Wang Xianchun,Hu Weiwen.De-noising of MT signal based on empirical mode decomposition and wavelet threshold method.OGP,2013,48(2):303-307.

[4]張翔,嚴(yán)良俊,蘇朱劉等.山區(qū)MT資料處理新技術(shù)及效果.石油地球物理勘探,2007,42(4):454-456.

Zhang Xiang,Yan Liangjun,Su Zhuliu et al.New technology and effects of mountainous MT data processing.OGP,2007,42(4):454-456.

[5]Goubau W M,Gamble T D,Clarke J.Magnetotelluric data analysis:removal of bias.Geophysics，1978,43(6):1157-1166.

[6]Larsen J C,Mackie R L,Manzella A et al.Robust smooth magnetotelluric transfer functions.Geophysical Journal International,2007,124(3):801-819.

[7]Yee E.The reconstruction of the magnetotelluric impedance tensor:An adaptive parametric time-domain approach.Geophysics,1988,53(3):1080-1087.

[8]Filipo W A S,Hohmann G W.Computer simulation of low-frequency electromagnetic data acquisition.Geophysics,1982,48(9):1219-1232.

[9]Lamarque G.Improvement of MT data processing using stationary and coherence tests.Geophysical Prospecting,1999,47(6):819-840.

[10]Ernst T.A comparison of two methods of the transfer function calculation using the least-square criterion in time and frequency domain.Publications of Institute of Geophysics Polish Academy of Sciences,1981,143:13-24.

[11]Kunetz G.Processing and interpretation of magnetotelluric soundings.Geophysics,1972,37(6):1005-1021.

[12]Mcmechan G A,Barrodale I.Processing electromagnetic data in the time domain.Geophysical Journal International,1985,81(1):277-293.

[13]Spagnolini U.Time-domain estimation of MT impe-dance tensor.Geophysics,1994,59(5):712-721.

[14]Chen J.Using empirical mode decomposition to process magnetotelluric data.Christian-Albrechts-Universit?t,2013.

[15]Loddo M,Schiavone D,Siniscalchi A.Generation of synthetic wide-band electromagnetic time series.Geophysics,2009,45(2):289-301.

[16]Varentsov I M,Sokolova E Y.Generation of synthetic magnetotelluric data.Izvestiya Physics of the Solid

Earth,1995,30(6):554-562.

[17]嚴(yán)家斌.大地電磁信號(hào)處理理論及方法研究[學(xué)位論文].湖南長沙：中南大學(xué),2003.

Yan Jiabin.The Study on Theory and Method of Magnetotelluric Signal Processing[D].Central South University,Changsha,Hu’nan,2003.

[18]蔡劍華,李晉.基于頻率域小波去噪的大地電磁信號(hào)工頻干擾處理.地質(zhì)與勘探,2015,51(2):353-359.

Cai Jianhua,Li Jin.Suppression of power line interfe-rence on MT signals based on the frequency domain wavelet method.Geology and Exploration,2015,51(2):353-359.

[19]陳知富,鄧居智,陳輝等.基于小波變換的大地電磁信號(hào)去噪研究.工程地球物理學(xué)報(bào),2012,9(6):732-737.

Chen Zhifu,Deng Juzhi,Chen Hui et al.De-noising research of magnetotelluric signal based on wavelet transform.Chinese Journal of Engineering Geophy-sics,2012,9(6)：732-737.

[20]凌振寶,王沛元,萬云霞等.強(qiáng)人文干擾環(huán)境的電磁數(shù)據(jù)小波去噪方法研究.地球物理學(xué)報(bào),2016,59(9):3436-3447.

Ling Zhenbao,Wang Peiyuan,Wan Yunxia et al.A combined wavelet transform algorithm used for de-noising magnetotellurics data in the strong human noise.Chinese Journal of Geophysics,2016,59(9)：3436-3447.

[21]Bracewell R,Kahn P B.The Fourier Transform and Its Applications.McGraw-Hill,2000.

[22]Rao K R,Kim D N,Hwang J J.Fast Fourier Transform - Algorithms and Applications.Springer Netherlands,2010.

[23]程乾生.數(shù)字信號(hào)處理.北京：北京大學(xué)出版社,2010.

Chen Qiansheng.Digital Signal Processing.Peking University Press,Beijing,2010.

[24]Labson V F,Becker A,Morrison H F et al.Geophysical exploration with audiofrequency natural magnetic fields.Geophysics,1985,50(4):656-664.

[25]Serson P H.Instrumentation for induction studies on land.Physics of the Earth & Planetary Interiors,1973,7(3):313-322.

[26]Filloux J H.Techniques and instrumentation for study of natural electromagnetic induction at sea.Physics of the Earth & Planetary Interiors,1973,7(3):323-338.

[27]Oppenheim A V,Schafer R W.Digital Signal Processing.Prentice-Hall,New Jersey，1975.

[28]Seshu S,Balabanian N.Linear Network Analysis.Wiley-Blackwell,New Jersey,1959.

[29]李金銘.地電場(chǎng)與電法勘探.北京：地質(zhì)出版社,2005.

[30]劉國棟,陳樂壽.大地電磁測(cè)深研究.北京：地震出版社,1984.