吳 楠 黃忠來

1)貴州師范大學(xué)機(jī)械與電氣工程學(xué)院，貴陽 550001

2)廈門大學(xué)信息科學(xué)與技術(shù)學(xué)院，廈門361005

1 引言

層狀介質(zhì)是探地雷達(dá)實(shí)際應(yīng)用中經(jīng)常遇到的介質(zhì)分布類型，當(dāng)?shù)叵履繕?biāo)層為相對厚層時，層的頂、底反射波可以在時域被區(qū)分，因此這類目標(biāo)層的厚度、深度以及電性參數(shù)可利用反射面反射波的振幅進(jìn)行估算［1，2］。但是，當(dāng)?shù)叵麓嬖诤穸刃∮谡{(diào)諧厚度的薄層時，頂、底反射波會在時域混疊在一起。這時，利用時域反演方法得到的結(jié)果將不再準(zhǔn)確。因此，頻率域反演方法被引入處理地下含有薄層的情況。頻率域反演方法需要建立地下介質(zhì)模型，構(gòu)造目標(biāo)函數(shù)，并利用線性或非線性反演方法進(jìn)行迭代計算［3，4］。當(dāng)?shù)叵履繕?biāo)層的個數(shù)較多時，如果要同時反演目標(biāo)層的幾何與電性參數(shù)，則需要反演的參數(shù)數(shù)量眾多，導(dǎo)致計算效率低下，且反演的結(jié)果容易收斂在局部極小值。通過采用接近真實(shí)值的迭代初值，或減少反演迭代的參數(shù)個數(shù)，都可以加快反演算法的收斂速度，提高反演結(jié)果的準(zhǔn)確度。在實(shí)際資料處理過程中，先驗(yàn)信息的獲取有時非常困難。雖然文獻(xiàn)［5－8］提出的多步驟的譜反演算法，可以提高反演的效率和準(zhǔn)確度。但是，算法中需要對地下反射面的分布有初步的劃分，而薄層識別則是難點(diǎn)。

為解決這一問題，本文結(jié)合譜反演算法，提出了一種利用倒頻譜技術(shù)的頻率域聯(lián)合反演算法。

圖1 地下包含N 個反射面的層狀介質(zhì)模型Fig.1 A layered media model containing N underground layers

2 算法原理

2.1 模型及目標(biāo)函數(shù)

對于圖1 中所示的水平層狀模型，設(shè)地下第i個反射面的反射系數(shù)為ri，若不考慮各反射面之間的多次波，則該反射面反射波的幅度Ai可以表示為［7］

其中Ainc是探地雷達(dá)入射波幅度為反射面的廣義反射系數(shù)，αk是電磁波在地下第k 層介質(zhì)中傳播時的衰減常數(shù)，dk是電磁波在第k 層介質(zhì)中的單向傳播距離。地表接收到的雷達(dá)回波s(t)可以表示成雷達(dá)子波和反射系數(shù)序列的褶積，

其中w(t)為雷達(dá)發(fā)射波，r(t)為地下介質(zhì)反射系數(shù)序列，n(t)為加性噪聲。若地下有N 個反射面，將反射系數(shù)序列中每相鄰兩個反射系數(shù)組成一個反射系數(shù)對，并對其做奇偶分解，如圖2 所示。因此反演算法的目標(biāo)函數(shù)可以寫為［7］:

由于反射面上的廣義反射系數(shù)是電性參數(shù)的函數(shù)，所以在反演時反射系數(shù)的奇偶分量可以用含有電性參數(shù)的表達(dá)式替換，這樣不僅可以反演出目標(biāo)層的物理屬性，還可以反演出所含介質(zhì)的電性參數(shù)。

2.2 基于倒頻譜的薄層識別

考慮地下只含有兩個反射系數(shù)，且目標(biāo)層為薄層情況。設(shè)上反射面的廣義反射系數(shù)為，下反射面的為;上反射面距離地表的時間距離為t1，下反射面的為t2。由于發(fā)射波的能量有限，且有噪聲存在，所以需要用窗口限制分析的帶寬。帶限的反射系數(shù)序列頻譜可以表示為:

其中rp(f)為反射系數(shù)序列頻譜，rRe(f)和rIm(f)分別是rp(f)的實(shí)部和虛部。注意到實(shí)部和虛部的相位相差π/2，且倒頻率相同項(xiàng)的符號相反，幅值相同。win(f)為窗函數(shù)，nwRe(f)和nwIm(f)分別是n(f)/w(f)項(xiàng)的實(shí)部和虛部。將的實(shí)部和虛部分別表示為和，即

序列頻譜的實(shí)部和虛部都是頻率域函數(shù)，它們包含相同的倒頻率成分，此時倒頻率即是反射系數(shù)的時間位置參數(shù)t1和t2，因此，在實(shí)部和虛部倒頻譜的t1和t2位置處會出現(xiàn)兩個鐘形尖峰。對于薄層來說，兩個鐘形尖峰會重合在一起，尖峰所對應(yīng)的倒頻率值tp滿足t1≤tp≤t2。由于反射系數(shù)頻譜實(shí)部和頻譜虛部對應(yīng)倒頻率成份的符號相反，相位相差π/2，因此如果移動虛部(實(shí)部)的某一倒頻率成分的相位，使之和實(shí)部(虛部)的相應(yīng)倒頻率成分相位一致，那么兩者相加后，這一倒頻率能量必然會消失。如果原先的頻譜實(shí)(虛)部中只存在一種倒頻率成分，那么兩者之和的頻譜中必然沒有能量峰值;如果原先頻譜實(shí)(虛)部中存在兩種倒頻率成分，那么兩者之和的頻譜中會有能量殘留。因此，根據(jù)能量是否殘留，就可以判斷倒頻率成分的數(shù)量，即反射面的個數(shù)。

在改變頻譜實(shí)(虛)部中某一倒頻譜成分的相位時，還需要考慮分析頻帶帶寬的影響。例如，我們的目的是改變反射系數(shù)頻譜虛部包含倒頻率t1項(xiàng)的相位，即將虛部中的sin(2πt1f)項(xiàng)轉(zhuǎn)化為實(shí)部中的cos(2πt1f)項(xiàng)。所以將反射系數(shù)序列頻譜虛部的頻譜乘以，再將結(jié)果進(jìn)行逆傅立葉變換，這里，相位因子τ=(4n +1)/4t1。為了避免在分析頻譜中引入較大噪聲，我們在限制分析頻帶后再進(jìn)行相位移動。分別利用窗口函數(shù)win(f)和winτ(f)= win(f－τ)對反射系數(shù)序列頻譜進(jìn)行截取。利用winτ(f)進(jìn)行分析頻帶限制的反射系數(shù)序列頻譜記為，

圖3 目標(biāo)為單反射面時的結(jié)果Fig.3 result of the target as a reflection interface

圖4 目標(biāo)為薄層時的結(jié)果Fig.4 result of the targ thin layer

3 實(shí)際資料處理

實(shí)際GPR 數(shù)據(jù)來自公路檢測資料，公路表面為瀝青，且瀝青面層分為3 層，厚度分別是0.04 m，0.06 m 和0.1 m。在瀝青面層之下還鋪了一層防水層。如圖5(a)所示，在8 ns 處的即為防水層。天線發(fā)射頻率為2 GHz，每道采集1 024 個數(shù)據(jù)，采集10 ns。我們以第15 道數(shù)據(jù)為例，來檢驗(yàn)算法的處理效果，數(shù)據(jù)的時域波形如圖5(b)所示。其實(shí)部倒頻譜如圖5(c)中虛線所示，可以看到，實(shí)部倒頻譜中出現(xiàn)了4 個峰值，其中第一個代表地表;第二個代表瀝青面上層和中間層的分界面反射波，峰值倒頻率為tp1;第三個代表中間層和下面層分界面的反射波，峰值倒頻率為tp2(圖5(d));而最后一個峰值則代表防水層，峰值倒頻率為tp3(圖5(e))。

圖5 探地雷達(dá)資料剖面圖及波形圖Fig.5 GPR data profile and waveform of the data

三個峰值倒頻率可以得到三個不同的相位因子，τ1、τ2和τ3，根據(jù)這三個因子對反射系數(shù)頻譜虛部進(jìn)行相位移動，分別得到和。計算它們的頻譜，結(jié)果如圖5(b)、(c)、(d)中的實(shí)線所示。從圖5(b)和(c)可以看到，由于這兩個峰值所指示的都是單一反射面，所以tp1 和tp2 處的能量在進(jìn)行相加后都消失了。而tp3 處是一個薄層，所以在和的頻譜中有明顯能量殘留。根據(jù)倒頻譜結(jié)果，我們用三個窗口劃分?jǐn)?shù)據(jù)，并對每個窗口數(shù)據(jù)進(jìn)行反演，反演采用隨機(jī)爬山算法。第一個窗口包含了地表回波和第一個反射面的回波。第二個窗口只包含第二個反射面的回波，第三個窗口包含薄層。圖6(a)和(b)分別是數(shù)據(jù)的時域剖面圖和反射面的反演結(jié)果。薄層的平均反演厚度為0.01m，大約為波長的1/5。從上倒下各層的介電常數(shù)分別為3.7、5.4、12、9.3 和13.2。電導(dǎo)率分別是0.013、0.02、0.018 和0.01。

圖6 時域剖面圖(a)和層面反演結(jié)果(b)Fig.6 GPR data profile(a)and inversion result of the layer interfaces(b)

4 結(jié)語

由于倒頻譜上的峰值指示了地下反射面的位置，因此，利用兩者對薄層進(jìn)行辨識，以解決譜反演算法中的層面識別問題。首先，將反射系數(shù)序列頻譜的虛部進(jìn)行相位移動，將結(jié)果與未進(jìn)行相位移動的頻譜實(shí)部進(jìn)行相加，隨后計算兩者之和的倒頻譜，通過觀察倒頻譜上峰值能量的變化，可以判斷峰值指示的是單反射面還是薄層。若和的倒頻譜在峰值處能量消失，則該處為單反射面;若峰值處能量有較大殘留，則該處為薄層。該方法可以為譜反演方法提供層面位置和數(shù)量的先驗(yàn)信息，從而提高反演算法的計算效率和結(jié)果準(zhǔn)確度。

1 AL-Qadi I L and Lahouar S.Measuring layer thicknesses with GPR-theory to practice［J］.Construction and Building Materials，2005，19:763－772.

2 Kao Chienping，et al.Measurement of layer thickness and permittivity using a new multilayer model from GPR data［J］.IEEE Transactions on Geoscience and Remote Sensing，2007，45(8):2 463－2 470.

3 Qin Yao，et al.Research on thin-layer recognition technique based on the spectrum inversion method of ground penetrating radar［J］.Journal of Electronics and Information Technology，2010，32(11):2 760－2 763.

4 Julien M，et al.A generalized frequency domain reflectometry modeling technique for soil electrical properties determination［J］.Vadose Zone Journal，2010，9(4):1 063－1 072.

5 Puryear C I and Castagna J P.Layer-thickness determination and stratigraphic interpretation using spectral inversion:Theory and application［J］.Geophysics，2008，73(2):37－48.

6 黃忠來，張建中.探地雷達(dá)薄層信號的譜反演算法［J］.大地測量與地球動力學(xué)，2011，(4):154－159.(Huang Zhonglai and Zhang Jianzhong.A spectral inversion algorithm for GPR signals of thin layers［J］.Journal of Geodesy and Geodynamics，2011，(4):154－159)

7 黃忠來，張建中.利用探地雷達(dá)頻譜反演層狀介質(zhì)幾何與電性參數(shù)［J］.地球物理學(xué)報，2013，56(4):1 381－1 391.(Huang Zhonglai and Zhang Jianzhong.An inversion method for inverting geometric and electric parameters of layered media using spectrum of GPR signal［J］.Chinese Journal of Geophysics，2013，56(4):1 381－1 391.

8 Giannopoulos A.Modeling ground penetrating radar by GprMax［J］.Construction and Building Materials，2005，19:755－762.