吳鴻卿，嚴(yán)家斌 ，王 慧，李大雁

(1．中南大學(xué) 地球科學(xué)與信息物理學(xué)院, 湖南 長沙 410083;2．中南大學(xué) 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，湖南 長沙 410083;3.核工業(yè)二三零研究所，湖南 長沙 410007)

1 引 言

大地電磁測(cè)深法(MT)自20世紀(jì)50年代問世以來，以其工作效率高、裝置輕便、勘探深度大、不受高阻層影響等優(yōu)點(diǎn)而受到人們的關(guān)注[1]，在礦產(chǎn)勘查、地下水與地?zé)峥碧?、油氣普查、地殼及巖石圈深部結(jié)構(gòu)探測(cè)等領(lǐng)域起著十分重要的作用。由于大地電磁信號(hào)源本身的不穩(wěn)定性、人類活動(dòng)及其電磁信號(hào)干擾和測(cè)量系統(tǒng)自身噪聲的存在，使得大地電磁信號(hào)具有非平穩(wěn)性、非線形、非最小相位等特性[2]。大地電磁測(cè)深方法的數(shù)據(jù)處理就是一個(gè)獲得準(zhǔn)確阻抗的過程，也是一個(gè)不斷去噪的過程[3]。Sims提出了最小二乘法估計(jì)阻抗方法，該方法在數(shù)據(jù)噪聲滿足獨(dú)立高斯分布規(guī)律時(shí)，最小二乘法估計(jì)為最優(yōu)無偏估計(jì)[4]。當(dāng)測(cè)量的數(shù)據(jù)存在異常點(diǎn)，在使用最小二乘法進(jìn)行阻抗估計(jì)時(shí)，結(jié)果往往會(huì)產(chǎn)生嚴(yán)重的偏差，研究人員提出了各種形式的Robust阻抗估計(jì)法[5]。Egber利用M回歸估計(jì)降低飛點(diǎn)的權(quán)，使其對(duì)阻抗函數(shù)的估計(jì)影響減小[6]。M回歸估計(jì)僅考慮了電道上存在了噪聲，但是磁道上也存在著噪聲，湯井田等提出了有界影響估計(jì)來解決這一問題[7]。Smirnov在重復(fù)中值估計(jì)的基礎(chǔ)上將其應(yīng)用到大地電磁阻抗估計(jì)中，保證了數(shù)據(jù)不受異常值和大誤差的影響[8]。張剛等將基于重復(fù)中位數(shù)估計(jì)的Robust法應(yīng)用長周期大地電磁阻抗估算中，取得了較好的效果[9]。對(duì)于復(fù)雜的海洋環(huán)境，常規(guī)的阻抗估計(jì)方法難以得出可靠的結(jié)果，柳建新等提出了根據(jù)相關(guān)系數(shù)的改進(jìn)Robust阻抗估算方法[10]。Chave等根據(jù)數(shù)據(jù)自身的噪聲分布進(jìn)行加權(quán)，提出了基于Stable分布的最大似然估計(jì)法，得出了可靠的阻抗估計(jì)值[11]。穩(wěn)健性估計(jì)方法是基于對(duì)未受干擾或者受干擾較小的數(shù)據(jù)對(duì)異常數(shù)據(jù)進(jìn)行降權(quán)，但是當(dāng)異常數(shù)據(jù)幅度很大或電道與磁道均存在噪聲時(shí)，就會(huì)使得阻抗結(jié)果產(chǎn)生偏倚。電道和磁道數(shù)據(jù)中都存在噪聲，Gambal等提出了遠(yuǎn)參考道陣列觀測(cè)方法，解決了磁場(chǎng)輸入信號(hào)受噪聲干擾的問題[12]。Larsen等[13]提出了穩(wěn)健光滑電磁的轉(zhuǎn)換函數(shù)，在遠(yuǎn)參考點(diǎn)存在的條件下，把Robust用到遠(yuǎn)參考估計(jì)中。Lesniak等把Kalman濾波方法應(yīng)用到大地電磁數(shù)據(jù)處理中，消除了不相關(guān)噪聲的影響[14]。湯井田等[15]，蔡劍華等[16]將HHT算法應(yīng)用到大地電磁信號(hào)分析中，短時(shí)Fourier變換、Wing-Ville分布和小波變換時(shí)頻分析方法，HHT方法更具適應(yīng)性,能更好地揭示電磁場(chǎng)的分布規(guī)律[16]。大地電磁信號(hào)分析，頻率域阻抗估計(jì)是基于觀測(cè)電磁信號(hào)的穩(wěn)定性，而對(duì)于非穩(wěn)定信號(hào),時(shí)域比頻率域估算阻抗具有更好的穩(wěn)定性和可靠性。時(shí)域上能很好表現(xiàn)出噪聲的局部特征,有利于信號(hào)的篩選，在時(shí)域上估算阻抗能表現(xiàn)出對(duì)噪聲的不敏感性,阻抗的穩(wěn)定估算只需短時(shí)穩(wěn)定信號(hào)。Yee等[17]用自適應(yīng)參數(shù)時(shí)域方法對(duì)阻抗張量進(jìn)行重構(gòu)，將觀測(cè)到的MT數(shù)據(jù)直接用于估計(jì)。Spagnolini[18]提出了一種基于實(shí)測(cè)電場(chǎng)和磁場(chǎng)數(shù)據(jù)的阻抗張量自適應(yīng)時(shí)域估計(jì)方法，自適應(yīng)估計(jì)在時(shí)域內(nèi)最小均方解，使實(shí)測(cè)數(shù)據(jù)與模擬數(shù)據(jù)的均方誤差最小，發(fā)現(xiàn)時(shí)域估計(jì)只需要短時(shí)間的平穩(wěn)性。本文從頻率域傳遞函數(shù)出發(fā)，利用傅里葉變換，將其轉(zhuǎn)化為時(shí)間域傳遞函數(shù)，得到時(shí)間域下的層狀介質(zhì)視電阻率表達(dá)式及時(shí)間域測(cè)深曲線，比較分析了時(shí)域阻抗估計(jì)的特點(diǎn)。

2 時(shí)域大地電磁測(cè)深阻抗估計(jì)

水平層狀介質(zhì)時(shí)域大地電磁測(cè)深表達(dá)式是從頻率域傳遞函數(shù)出發(fā)并對(duì)其進(jìn)行求解，之后利用傅里葉變換，將其轉(zhuǎn)化為時(shí)間域傳遞函數(shù)，得到時(shí)間域?qū)訝罱橘|(zhì)視電阻率公式。

2.1 時(shí)間域電磁場(chǎng)傳遞函數(shù)

時(shí)間域中聯(lián)系電磁場(chǎng)分量的傳遞函數(shù)[19]記為D(t)，定義為：

(1)

式(1)中,*表示褶積;t表示時(shí)間，E(t)是時(shí)間域電場(chǎng)分量lH(t)是時(shí)間域磁場(chǎng)分量。對(duì)于D(t)的求解，需要從與其對(duì)應(yīng)的頻率域的傳遞函數(shù)出發(fā)，其表達(dá)式為：

E(ω)=D(ω)·[-iωH(ω)]

(2)

式(2)中,D(ω)為頻率域傳遞函數(shù);ω表示頻率;E(ω)頻率域電場(chǎng)分量;H(ω)是頻率域磁場(chǎng)分量。E(ω)由大地電磁場(chǎng)的阻抗關(guān)系可表示為[19]：

(3)

其中,Z(ω)是阻抗，將式(2)和式(3)相比得到：

(4)

其中μ是導(dǎo)磁率;σ1是導(dǎo)電率。式中Rn(ω)[19]為頻率域下的N層介質(zhì)的阻抗函數(shù):

(5)

(6)

利用傅氏變換公式：

(7)

可求出D(ω)的傅氏反變換為：

(8)

ρ1為第一層電阻率，其中

其中

且有

Ln-1=Kn-1,n

圖1 一維模型分層后的地層參數(shù)Fig.1 Stratigraphic parameter map after layering of one-dimensional model

(10)

2.2 時(shí)域電阻率估計(jì)

由式(8)得到：

(11)

得到時(shí)間域歸一化視電阻率為：

ρa(bǔ)=ρ1×S2(t)=4π2tD2(t)

(12)

3 模型分析

建立模型對(duì)時(shí)間域視電阻率估計(jì)的方法進(jìn)行驗(yàn)證，并分析子層綜合參數(shù)以及劃分的子層數(shù)對(duì)時(shí)域電阻率的影響。

3.1 模型建立與方法驗(yàn)證

給定三層地質(zhì)模型ρ1=900 Ω·m,ρ2=400 Ω·m,ρ3=900 Ω·m,h1=900 m,h2=600 m,h3=,令各子層的綜合參數(shù)則對(duì)于第一層分為30個(gè)子層，每層層厚30 m；第二層分為30個(gè)子層，每層厚20 m，如圖2所示。分析時(shí)間范圍10-5～102s，計(jì)算得到的時(shí)域電阻率如圖3所示。從圖中可以看到時(shí)域曲線與頻率曲線一致，首支曲線和尾支曲線與模型的第一層和第三層電阻率相同，證明了該算法的可靠性。

圖3 H型三層介質(zhì)時(shí)間域視電阻率曲線Fig.3 H-type three-layer dielectric time domain apparent resistivity curve

3.2 模型分析

3.2.1 二層G型模型

由圖4可以看出，當(dāng)子層綜合參數(shù)c給定的時(shí)候，所取的子層數(shù)a不同，所得到的曲線不一樣，只要滿足a×c=350的條件，所得到視電阻率曲線(圖4)接近理論曲線(圖5b)。當(dāng)子層綜合參數(shù)取1，所取子層數(shù)為350時(shí)，所得到的視電阻率曲線最接近理論曲線(圖4d)。

圖4 不同子層參數(shù)的G型時(shí)間域視電阻率曲線Fig.4 G-type time domain apparent resistivity curve of different sublayer parameters注：c為子層綜合參數(shù)，a為選取的總的層數(shù)

圖5 不同子層參數(shù)的G型時(shí)間域視電阻率曲線對(duì)比和G型理論曲線Fig.5 Comparison of G-type apparent resistivity curves in time domain with different sub-layer parameters and G-type theoretical curves

3.2.2 三層H型模型

由圖6可以看出，在相同的H型模型下，當(dāng)子層綜合參數(shù)c給定的時(shí)候，所取的子層數(shù)a不同，所得到的理論曲線不一樣，只要滿足a×c=330的條件，所得到視電阻率曲線(圖6)就接近理論曲線(圖7b)。當(dāng)子層綜合參數(shù)取1時(shí)，所取子層數(shù)為330時(shí)，所得到的視電阻率曲線最接近理論曲線(圖6d)。

圖6 不同子層參數(shù)的H型時(shí)間域視電阻率曲線Fig.6 H-type time domain apparent resistivity curve of different sublayer parameters

圖7 不同子層參數(shù)的H型時(shí)域視電阻率曲線對(duì)比及H型理論曲線Fig.7 Comparison of H-type apparent resistivity curves in time domain with different sub-layer parameters and H-type theoretical curves

4 大地電磁測(cè)深曲線的時(shí)域與頻域比較

4.1 G型曲線

模型參數(shù)取ρ1=400 Ω·m，ρ2=900 Ω·m，h1=1 000 m。

由圖8可見，頻率域和時(shí)間域兩條曲線在10-3s以前即首支曲線重合，與第一層的電阻率400 Ω·m一致；隨著時(shí)間的增加或頻率的降低，視電阻率逐漸增大，頻率域和時(shí)間域兩條曲線的尾支漸近線，與第二層真實(shí)視電阻率相近。時(shí)間域曲線在10-3s開始受第二層電阻率的影響，視電阻率開始增大；而頻率域曲線要在3×10-3s 才受第二層的影響；同時(shí)時(shí)間域曲線在1 s時(shí)電阻率達(dá)到第二層的真電阻，而頻率域曲線10 s還沒達(dá)到第二層的真電阻率，表明時(shí)間域電阻率相對(duì)于頻率域有更大的探測(cè)深度；頻率域曲線大約在3×10-3～5×10-2s出現(xiàn)了向下凹，出現(xiàn)假的低阻異常，這可能是由于界面反射波和入射波干涉引起的[20,21]，而在時(shí)間域曲線中并沒有出現(xiàn)。

圖8 G型二層介質(zhì)時(shí)頻域視電阻率對(duì)比Fig.8 Comparison of apparent resistivity of frequency domain in G-type two-layer medium

4.2 D型曲線

模型參數(shù)取ρ1=1 600 Ω·m，ρ2=400 Ω·m，h1=2 000 m。

由圖9可見，頻率域和時(shí)間域兩條曲線在10-3s以前即首支曲線重合，與第一層的電阻率1 600 Ω·m一致；隨著時(shí)間的增加或頻率的降低，視電阻率逐漸減小，頻率域和時(shí)間域兩條曲線的支漸近線，與第二層真實(shí)視電阻率相近。時(shí)間域曲線在10-3s開始受第二層電阻率的影響，頻率域曲線要在3×10-3s才受第二層的影響；同時(shí)時(shí)間域曲線在10-1s時(shí)電阻率達(dá)到第二層的真電阻率，而頻率域曲線100 s還沒達(dá)到第二層的真電阻率，表明時(shí)間域電阻率相對(duì)于頻率域有更大的探測(cè)深度；頻率域曲線大約在3×10-3～8×10-2s時(shí)間段出現(xiàn)向上隆起，出現(xiàn)假的低阻異常，同樣在時(shí)間域曲線中并沒有出現(xiàn)。

4.3 H型曲線

模型參數(shù)取ρ1=900 Ω·m，ρ2=400 Ω·m，ρ3=900 Ω·m，h1=400 m，h2=600 m。

圖9 D型二層介質(zhì)時(shí)頻域視電阻率對(duì)比Fig.9 Comparison of apparent resistivity of frequency domain in D-type two-layer medium

圖10 H型三層介質(zhì)時(shí)頻域視電阻率對(duì)比Fig.10 Comparison of time-frequency domain apparent resistivity of H-type three-layer medium

圖11 A型三層介質(zhì)時(shí)頻域視電阻率對(duì)比Fig.11 Comparison of the apparent resistivity of the frequency domain in the A-type three-layer medium

由圖10可見，頻率域和時(shí)間域兩條曲線在3×10-4s以前即首支曲線重合，與第一層的電阻率900 Ω·m一致；隨著時(shí)間的增加或頻率的降低，視電阻率逐漸減小，到達(dá)一個(gè)極小值后，又逐漸增加，頻率域和時(shí)間域曲線都存在尾支漸近線，與第三層真實(shí)視電阻率相近，時(shí)間域曲線在3×10-4s開始受第二層電阻率的影響，頻率域曲線要在2×10-3s才受第二層的影響；同時(shí)時(shí)間域曲線在1 s時(shí)電阻率達(dá)到第三層的真電阻率，而頻率域曲線100 s還沒達(dá)到第三層的真電阻率，表明時(shí)間域電阻率相對(duì)于頻率域有更大的探測(cè)深度；時(shí)間域曲線在2×10-3s左右達(dá)到極大值，而頻率域曲線在0.2 s左右才達(dá)到極小值且大于時(shí)間域曲線的極小值，表明時(shí)間域電阻率相對(duì)于頻率域電阻率更逼近中間層電阻率值；頻率域曲線大約在2×10-3～7×10-3s時(shí)間段出現(xiàn)向上隆起，出現(xiàn)假的低阻異常，同樣在時(shí)間域曲線中并沒有出現(xiàn)。

4.4 A型曲線

模型參數(shù)取ρ1=400 Ω·m，ρ2=1 600 Ω·m，ρ3=2 400 Ω·m，h1=1 200 m，h2=2 400 m。

由圖11可知，頻率域和時(shí)間域兩條曲線在10-3s以前即首支曲線重合，與第一層的電阻率400 Ω·m一致；隨著時(shí)間的增加或頻率的降低，視電阻率逐漸增加，頻率域和時(shí)間域兩條曲線都存在尾支漸近線，與第三層真實(shí)視電阻率相近，時(shí)間域曲線在10-3s開始受第二層電阻率的影響，頻率域曲線要10-2s才受第二層的影響；同時(shí)時(shí)間域曲線在1 s時(shí)電阻率達(dá)到第二層的真電阻率，而頻率域曲線100 s還沒達(dá)到第二層的真電阻率，表明時(shí)間域電阻率相對(duì)于頻率域有更大的探測(cè)深度；頻率域曲線大約在10-2～10-1s時(shí)間段出現(xiàn)向下凹，出現(xiàn)假的低阻異常，同樣在時(shí)間域曲線中并沒有出現(xiàn)。

4.5 Q型曲線

模型參數(shù)取ρ1=1 600 Ω·m，ρ2=900 Ω·m，ρ3=400 Ω·m，h1=2 400 m，h2=1 800 m。

由圖12可見，頻率域和時(shí)間域兩條曲線在10-3s以前即首支曲線重合，與第一層的電阻率1 600 Ω·m一致；隨著時(shí)間的增加或頻率的降低，視電阻率逐漸減小，頻率域和時(shí)間域兩條曲線都存在尾支漸近線，與第三層真實(shí)視電阻率相近，時(shí)間域曲線在10-3s開始受第二層電阻率的影響，頻率域曲線要5×10-2s才受第二層的影響；同時(shí)時(shí)間域曲線在1 s時(shí)電阻率達(dá)到第二層的真電阻率，而頻率域曲線100 s還沒達(dá)到第二層的真電阻率，表明時(shí)間域電阻率相對(duì)于頻率域有更大的探測(cè)深度；頻率域曲線大約在5×10-2s～5×10-1s時(shí)間段出現(xiàn)向上隆起，出現(xiàn)假的低阻異常，同樣在時(shí)間域曲線中并沒有出現(xiàn)。

圖12 Q型三層介質(zhì)時(shí)頻域視電阻率對(duì)比Fig.12 Comparison of the apparent resistivity of the Q-type three-layer medium

4.6 K型曲線

模型參數(shù)取ρ1=400 Ω·m，ρ2=1 600 Ω·m，ρ3=400 Ω·m，h1=1 200 m，h2=2 400 m。

圖13 K型三層介質(zhì)時(shí)頻域視電阻率對(duì)比Fig.13 Comparison of apparent resistivity of frequency domain in K-type three-layer medium

由圖13，頻率域和時(shí)間域兩條曲線在10-3s以前即首支曲線重合，與第一層的電阻率400 Ω·m一致；隨著時(shí)間的增加或頻率的降低，視電阻率逐漸增加，到達(dá)一個(gè)極大值后，又逐漸減小，頻率域和時(shí)間域兩條曲線都存在尾支漸近線，與第三層真實(shí)視電阻率相近，時(shí)間域曲線在10-3s開始受第二層電阻率的影響，頻率域曲線要5×10-2s才受第二層的影響；同時(shí)時(shí)間域曲線在1 s時(shí)電阻率達(dá)到第二層的真電阻率，而頻率域曲線100 s還沒達(dá)到第二層的真電阻率，表明時(shí)間域電阻率相對(duì)于頻率域有更大的探測(cè)深度；時(shí)間域曲線在2×10-2s左右達(dá)到極大值，而頻率域曲線在2 s左右才達(dá)到極大值且小于時(shí)間域曲線的極大值，表明時(shí)間域電阻率相對(duì)于頻率域電阻率更逼近中間層電阻率值；頻率域曲線大約在5×10-2～5×10-1s時(shí)間段出現(xiàn)向下凹，出現(xiàn)假的低阻異常，同樣在時(shí)間域曲線中并沒有出現(xiàn)。

4.7 識(shí)別薄層能力的分析

對(duì)于三層介質(zhì)來說，就識(shí)別中間薄層電阻率的能力進(jìn)行分析。通過改變第二層的厚度，將得到的時(shí)間域和頻率域視電阻率曲線進(jìn)行分析：

1)模型參數(shù)取ρ1=1 600 Ω·m，ρ2=400 Ω·m，ρ3=3 000 Ω·m，h1=1 000 m，h2取1 000 m、500 m、100 m和50 m。

圖14中(a)、(b)、(c)、(d)分別是第二層厚度取1 000 m、500 m、100 m和50 m的時(shí)頻域視電阻率對(duì)比圖。從圖14中可以看出，在中間層厚度為1 000 m和500 m時(shí)，時(shí)間域曲線中反映的中間層視電阻率值相較于頻率域曲線的中間層視電阻率值更接近中間層真電阻率；當(dāng)中間層厚度為100 m和50 m這樣的薄層的時(shí)候，頻率域曲線中反映的中間層視電阻率值相較于時(shí)間域曲線的中間層視電阻率值更接近中間層真電阻率，由此可以得出，時(shí)間域視電阻率識(shí)別中間薄層的能力可能沒有頻率域視電阻率的強(qiáng)。

圖14 H型三層介質(zhì)不同中間層厚度的時(shí)頻域視電阻率對(duì)比Fig.14 Comparison of time-frequency domain apparent resistivity of different intermediate layer thicknesses of H-type three-layer medium

圖15 K型三層介質(zhì)不同中間層厚度的時(shí)頻域視電阻率對(duì)比Fig.15 Comparison of time-frequency domain apparent resistivity of different intermediate layer thicknesses of K-type three-layer medium

2)模型參數(shù)取ρ1=1 600 Ω·m ，ρ2=6 400 Ω·m，ρ3=1 000 Ω·m，h1=800 m，h2取800 m、400 m、160 m和80 m。

圖15中(a)、(b)、(c)、(d)分別是第二層厚度取800 m、400 m、160 m和80 m的時(shí)頻域視電阻率對(duì)比圖。從圖5中可以看出，在中間層厚度為800 m和400 m時(shí)，時(shí)間域曲線中反映的中間層視電阻率值相較于頻率域曲線的中間層視電阻率值更接近中間層真電阻率；當(dāng)中間層厚度為160 m和80 m這樣的薄層的時(shí)候，頻率域曲線中反映的中間層視電阻率值相較于時(shí)間域曲線的中間層視電阻率值更接近中間層真電阻率，從中發(fā)現(xiàn)時(shí)間域視電阻率識(shí)別中間薄層的能力可能沒有頻率域視電阻率的強(qiáng)。

5 結(jié) 論

1)通過模型分析表明，在所分層數(shù)相同的情況下，子層綜合參數(shù)越大，即分的層數(shù)越多，計(jì)算得到的視電阻率越精確。

2)在相同頻率下，時(shí)間域視電阻率比頻率域視電阻率探測(cè)深度更大；由于界面反射波和入射波干涉，頻率域視電阻率曲線在高頻段會(huì)出現(xiàn)假異?，F(xiàn)象，而時(shí)間域視電阻率曲線并未出現(xiàn)；可以看出時(shí)間域比頻率域估算阻抗或電阻率更準(zhǔn)確。

3)通過改變第二層厚度，對(duì)所建立的K型和H型三層介質(zhì)模型進(jìn)行分析，發(fā)現(xiàn)時(shí)間域視電阻率識(shí)別中間薄層的能力可能沒有頻率域視電阻率的強(qiáng)。