張寶強　裴建新*②　王　啟

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

1　引言

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

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

2　方法原理

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

2.1　電磁場頻譜的求取

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

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

2.2　電磁場時間序列的計算

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

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

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

(1)

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

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

圖1　模擬MT時間序列的基本流程

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

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

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

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

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

(2)

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

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

H(s)=

(3)

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

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

表1　模擬磁場使用的極點和零點的頻率絕對值和特性

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

(4)

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

3.2　時間域磁場數(shù)據(jù)模擬

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

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

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

圖4　模擬的磁場頻譜和時間序列，頻率范圍為0.0002～1Hz，時間長度為5000s

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

4.1　一維層狀模型阻抗計算

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

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

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

(5)

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

(6)

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

4.2　時間域電場數(shù)據(jù)模擬

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

圖6　一維水平層狀模型

圖7　模型阻抗的實部和虛部

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

4.3　模擬數(shù)據(jù)與實測數(shù)據(jù)對比

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

圖8　模擬的電場頻譜和時間序列

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

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

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

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

(7)

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

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

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

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

圖11　模擬二維MT水平分量時間序列

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

6　結(jié)束語

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

[1]張全勝,王家映.大地電磁測深資料的去噪方法.石油地球物理勘探,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]湯井田,張弛,肖曉等.大地電磁阻抗估計方法對比.中國有色金屬學(xué)報,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)驗?zāi)B(tài)分解與小波閾值的MT信號去噪方法.石油地球物理勘探,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)家斌.大地電磁信號處理理論及方法研究[學(xué)位論文].湖南長沙：中南大學(xué),2003.

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

[18]蔡劍華,李晉.基于頻率域小波去噪的大地電磁信號工頻干擾處理.地質(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]陳知富,鄧居智,陳輝等.基于小波變換的大地電磁信號去噪研究.工程地球物理學(xué)報,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]凌振寶,王沛元,萬云霞等.強人文干擾環(huán)境的電磁數(shù)據(jù)小波去噪方法研究.地球物理學(xué)報,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ù)字信號處理.北京：北京大學(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]李金銘.地電場與電法勘探.北京：地質(zhì)出版社,2005.

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