帶有低通濾波的廣義S變換在探地雷達層位識別中的應(yīng)用

張先武，高云澤，方廣有＊

1 中國科學(xué)院電子學(xué)研究所，北京 100190

2 中國科學(xué)院大學(xué)，北京 100049

1 引 言

在淺層地球物理勘探中，探地雷達以其無損、高效、高分辨率等優(yōu)點正越來越受到人們的重視[1-5］.作為探地雷達的常見勘探目標，如冰川、公路路基等，通?？杀灰暈閷訝罱橘|(zhì)[6-9］.準確識別這些層狀介質(zhì)的層位信息（深度、厚度等）對后續(xù)處理和解釋探地雷達數(shù)據(jù)至關(guān)重要[10］.電磁波在這些層狀介質(zhì)中傳播時，會產(chǎn)生衰減和頻散現(xiàn)象[11-15］，造成探地雷達接收到的回波信號為時變信號.利用回波信號的時變特性可以為識別層狀介質(zhì)的層位信息提供幫助.

對于時變信號，人們提出了一系列信號分析和處理方法，如短時Fourier變換（STFT）[16］、小波變換（WT）[17］、S變換[18］等.同 WT、STFT 等時間-頻率域信號處理方法相比，S變換具有很多獨特優(yōu)點，如時間-頻率域的分辨率與頻率相關(guān)；S變換譜與Fourier變換譜保持直接聯(lián)系；S變換的基本小波不用滿足容許性條件等[19-21］.

探地雷達勘探中，雷達信號穿過由層厚小于調(diào)諧厚度的薄層構(gòu)成的層序列時，反射的雷達信號的頻率會增高[22-24］.依據(jù)此特點，利用S變換可以對這種薄層序列進行有效識別[25］.為了獲得層的厚度、深度等信息，需要S變換具有較高的時間分辨率[26-27］.而S變換的時窗函數(shù)固定，不能根據(jù)需要來調(diào)節(jié)時間分辨率.對此，本文對S變換的時窗函數(shù)加以改進，使其可調(diào).針對變換中提高時間分辨率會降低頻率分辨率的問題，本文將低通濾波函數(shù)融入到時窗函數(shù)中，來調(diào)節(jié)變換中的頻率分辨率.通過算例，文中對改進后的廣義S變換與S變換進行了對比分析，同時將廣義S變換應(yīng)用到實測探地雷達資料的層位識別中，提高了S變換的層位識別能力.

2 基本原理

2.1 S變換

Stockwell給出函數(shù)h（t）的S變換定義如下[18］：

式中，t，τ表示時間，f表示頻率，均為實數(shù).

定義S變換的時窗函數(shù)為G（t，f），

S變換的時窗函數(shù)G（t，f）要滿足以下條件[28］，即

在滿足（3）式的條件下，可以得到：

（4）式中，H（f）為函數(shù)h（t）的Fourier變換.由此可以得到S變換的反變換為：

為了提高計算效率，通常不是直接按照（1）式在時間域?qū)崿F(xiàn)S變換，而是在頻率域?qū)崿F(xiàn).頻率域S正變換公式為：

2.2 廣義S變換

為了使廣義S變換的時窗函數(shù)可調(diào)，將低通濾波函數(shù)融入其中，同時在時窗函數(shù)中引入調(diào)節(jié)參數(shù).通過改變調(diào)節(jié)參數(shù)的大小來調(diào)節(jié)廣義S變換的時間分辨率和頻率分辨率.

定義廣義S變換的時窗函數(shù)為W（t，f，λG，λL），W （t，f，λG，λL）=

式中，λG，λL表示的是不同的調(diào)節(jié)參數(shù).λG主要用來調(diào)節(jié)時窗的寬度，λG越大，時窗越窄，廣義S變換的時間分辨率越高，但是同時會降低變換的頻率分辨率.此時再通過調(diào)節(jié)λL，對廣義S變換的頻率分辨率進行調(diào)節(jié)，λL越小，廣義S變換的頻率分辨率越高.

廣義S變換的時窗函數(shù) W （t ，f，λG，λL）滿足S變換的時窗函數(shù)條件，即

定義廣義S變換如下：

在滿足（8）式的條件下，可以得到：

（10）式中，H（f）為函數(shù)h（t）的 Fourier變換.由此可以得到廣義S變換的反變換為：

同S變換類似，為了提高計算效率，可以在頻率域?qū)崿F(xiàn)廣義S變換.變換公式為：

（12）式中L（ζ，f）為：為了保證廣義S正、反變換完全可逆.當(dāng)f=0時，廣義S正變換滿足：

在計算機上實現(xiàn)廣義S變換時，采用（12）和（13）式的離散形式.

3 對比分析

為了研究廣義S變換分析時變信號的特點，合成一個時變信號h（t）.h（t）的解析式為：

信號h（t）主要包含三個頻率成分，分別為100Hz、200Hz、300Hz.各個頻率成分起止時間對應(yīng)分別為：0～449ms、274～723ms、550～999ms.計算機上處理時，需要先對上述信號h（t）進行采樣，本文采用的采樣時間間隔為1ms，總的采樣時間為999ms.圖1為信號h（t）采樣后的示意圖.

圖1 合成時變信號采樣后示意圖Fig.1 Sketch of sampling results of synthetic time-varying signal

3.1 合成信號的S變換

為了對比分析廣義S變換和S變換，對信號h（t）進行S變換，圖2為S變換結(jié)果.圖2中，信號的頻率成分以及各個頻率成分的起止時間能很好地與合成信號相對應(yīng).隨著信號頻率的增大，S變換的頻率分辨率降低，時間分辨率增高，這同小波變換相似.

圖2 合成時變信號S變換結(jié)果Fig.2 S transform results of synthetic time-varying signal

3.2 合成信號的廣義S變換

3.2.1 調(diào)節(jié)參數(shù)λG對廣義S變換的影響

為了分析調(diào)節(jié)參數(shù)λG對廣義S變換的影響，先固定參數(shù)λL（文中固定λL為10），再調(diào)節(jié)λG.對不同的λG，分別對信號h（t）進行廣義S變換.圖3顯示了λG對廣義S變換的影響.從圖3中可以看到，在λL不變的情況下，隨著λG的增大，廣義S變換的時間分辨率在增大，λG=1時，廣義S變換結(jié)果與S變換結(jié)果相同，當(dāng)λG＞1時，廣義S變換的時間分辨率比S變換的高.但是，隨著λG的增大，廣義S變換的頻率分辨率在減小.λG=2時，由于頻率分辨率過低，在550～723ms時間段內(nèi)，200Hz和300Hz頻率成分相互重疊，此時，在廣義S變換后的時間-頻率域內(nèi)已經(jīng)無法將它們識別出來.

3.2.2 調(diào)節(jié)參數(shù)λL對廣義S變換的影響

固定λG（文中固定λG為2），調(diào)節(jié)λL，對不同的λL，分別對信號h（t）進行廣義S變換，變換結(jié)果為圖4.

對比圖4a和圖3d，兩者變換結(jié)果相同，在550～723ms時間段內(nèi)，無法將200Hz和300Hz頻率成分識別出來.隨著λL的減小，廣義S變換的頻率分辨率在增大，λL≤0.2時，200Hz和300Hz頻率成分逐漸分離，此時在時間-頻率域內(nèi)已經(jīng)可以將它們識別出來.由圖4a和圖3d所對應(yīng)的變換結(jié)果相同，表明在λG不變的情況下，λL增大到某一數(shù)值后，廣義S變換結(jié)果不發(fā)生改變.當(dāng)λG=1，λL=+∞時，對任何信號，廣義S變換與S變換結(jié)果相同.

4 廣義S變換識別實測探地雷達資料的層位

圖5為某一湖區(qū)實測探地雷達資料.測量時，天線中心頻率選用150MHz，采用的時間窗口為500ns，采樣點數(shù)為512.總的記錄道數(shù)為141道.本文對其中的第90道（圖6）進行分析，分別對其進行S變換和廣義S變換，圖7為S變換結(jié)果.

對實測探地雷達數(shù)據(jù)進行廣義S變換時，需要確定合適的調(diào)節(jié)參數(shù)λG和λL.選取λG和λL時，可以先確定λL（一般可令λL為1），再根據(jù)變換結(jié)果來調(diào)節(jié)λG，λG越大，廣義S變換的時間分辨率越高.調(diào)節(jié)λG時，可先取λG為1，再逐漸增大λG，對比分析不同λG對應(yīng)的廣義S變換結(jié)果.當(dāng)變換結(jié)果的時間分辨率隨著λG的增大而變化不大時，選取此時的λG.

在選取λG的過程中，為了獲得較高的時間分辨率，逐漸增大λG，但是同時會降低廣義S變換的頻率分辨率.為了對廣義S變換的頻率分辨率進行一定的補償，在確定調(diào)節(jié)參數(shù)λG后，再適當(dāng)減小λL，當(dāng)變換結(jié)果中的不同頻率成分能較好地分開時，選取此時的λL.λL不宜過小，否則會影響廣義S變換的時間分辨率.圖8為第90道數(shù)據(jù)的廣義S變換結(jié)果，變換時選取的調(diào)節(jié)參數(shù)分別為λG=2、λL=0.5和λG=5、λL=0.5.

圖5 實測探地雷達資料Fig.5 Measured ground penetrating radar data

圖6 實測探地雷達資料第300道記錄Fig.6 Record of trace 300in ground penetrating radar data

從S變換結(jié)果（圖7）中可以看到，105～135采樣點間的回波信號中心頻率在200MHz左右，相對于直達波的中心頻率（130MHz左右）有明顯提升.據(jù)此可以推斷105～135采樣點間的回波信號對應(yīng)的為湖底淤泥沉積序列的反射回波信號，且淤泥沉積層的單層厚度小于為湖底淤泥的相對介電常數(shù)）.

圖7 第300道記錄的S變換結(jié)果Fig.7 S transform result of the record of trace 300

對比分析圖7、圖8，105～135采樣點間的回波信號在廣義S變換結(jié)果（圖8）中，可以被更精細地劃分為兩部分，前一部分信號的中心頻率（相對于直達波的中心頻率）有明顯提升，而后一部分信號的中心頻率變化不大.從而我們可以將105～135采樣點間的回波信號對應(yīng)的地下介質(zhì)劃分為兩部分，前一部分對應(yīng)的是湖底淤泥沉積序列，后一部分對應(yīng)的是沉積序列底部的獨立反射層.

為了更好地驗證單道（第90道）雷達記錄廣義S變換分析結(jié)果，本文對實測的141道雷達數(shù)據(jù)分別進行S變換和廣義S變換（變換時，取λG=5、λL=0.5）.再對變換結(jié)果進行頻率切片分析（文中選用的頻率為150MHz），圖9和圖10分別為S變換和廣義S變換頻率切片分析結(jié)果.從圖10（A區(qū)域）中可以明顯地看出，105～135采樣點間地下介質(zhì)可被劃分為湖底淤泥沉積序列（A1）和沉積序列底部的獨立反射層（A2）兩部分，而在圖9的相同位置，無法劃分出沉積序列底部的獨立反射層.頻率切片分析結(jié)果也表面廣義S變換提高了S變換的層位劃分能力.

對比分析圖7、圖8，探地雷達回波信號經(jīng)過廣義S變換后，回波信號在時頻平面內(nèi)的能量分布沿時間方向上變窄，時間分辨率在增大，這一點也可以從頻率切片分析結(jié)果中得到印證.對比分析圖9和圖10的B、C區(qū)域，圖10中層位的反射回波信號清晰，能量分布沿時間方向上細窄，可以更精細地確定層位的厚度和深度信息.

5 結(jié) 論

本文在S變換的基礎(chǔ)上，對S變換加以改進和推廣，提出了一種廣義S變換.通過引入調(diào)節(jié)參數(shù)，同時將低通濾波函數(shù)融入到廣義S變換的時窗函數(shù)中，來調(diào)節(jié)廣義S變換的時間分辨率和頻率分辨率.利用廣義S變換分析時變信號時，可以根據(jù)需要來選擇合適的調(diào)節(jié)參數(shù).理論信號分析和實測數(shù)據(jù)計算表面該變換具有很強的適用性.

廣義S變換可以在頻率域?qū)崿F(xiàn)，有較高的計算效率.實際應(yīng)用中，廣義S變換的離散采樣條件與S變換相同.在實測探地雷達資料的層位識別中，同S變換相比，廣義S變換在確定層位深度、厚度和劃分層位等方面的能力得到提高.

（References）

[1]Annan A P.Ground penetrating radar：Principles，procedures，＆applications，Sensors ＆Software Inc.Technical Paper，2003.

[2]陳義群，肖柏勛.論探地雷達現(xiàn)狀與發(fā)展.工程地球物理學(xué)報，2005，2（2）：149-155.Chen Y Q，Xiao B X.On the status quo and development of Ground Penetrating Radar.Chinese Journal of Engineering Geophysics （in Chinese），2005，2（2）：149-155.

[3]李嘉，郭成超，王復(fù)明等.探地雷達應(yīng)用概述.地球物理學(xué)進展，2007，22（2）：629-637.Li J，Guo C C，Wang F M，et al.The summary of the surface ground penetrating radar applied in subsurface investigation.Progress in Geophysics（in Chinese），2007，22（2）：629-637.

[4]Slob E，Sato M，Olhoeft G.Surface and borehole groundpenetrating-radar developments.Geophysics，2010，75（5）：A103-A120.

[5]李華，魯光銀，何現(xiàn)啟等.探地雷達的發(fā)展歷程及其前景探討.地球物理學(xué)進展，2010，25（4）：1492-1502.Li H，Lu G Y，He X Q，et al.The progress of the GPR and discussion on its future development.Progress in Geophysics（in Chinese），2010，25（4）：1492-1502.

[6]武震，張世強，劉時銀等.祁連山老虎溝12號冰川冰內(nèi)結(jié)構(gòu)特征分析.地球科學(xué)進展，2011，26（6）：631-641.Wu Z，Zhang S Q，Liu S Y，et al.Structural characteristics of the No.12Glacier in Laohugou Valley，Qilian Mountain based on the ground penetrating radar combined with FDTD simulation.Advances in Earth Science （in Chinese），2011，26（6）：631-641.

[7]楊天春，呂紹林，王齊仁.探地雷達檢測道路厚度結(jié)構(gòu)的應(yīng)用現(xiàn)狀及進展.物探與化探，2003，27（1）：79-82.Yang T C，LüS L，Wang Q R.The application and development of GPR in detection of pavement thickness and highway structures.Geophysical and Geochemical Exploration（in Chinese），2003，27（1）：79-82.

[8]馮德山，戴前偉.高速公路路面厚度探地雷達檢測.地球物理學(xué)進展，2008，23（1）：289-294.Feng D S，Dai Q W.Application of ground penetrating radar in the survey of the pavement thickness in highway.Progress in Geophysics（in Chinese），2008，23（1）：289-294.

[9]秦瑤，陳潔，方廣有等.基于探地雷達頻譜反演法的薄層識別技術(shù)研究.電子與信息學(xué)報，2010，32（11）：2760-2763.Qin Y，Chen J，F(xiàn)ang G Y，et al.Research on thin-layer recognition technique based on the spectrum inversion method of Ground Penetrating Radar.Journal of Electronics &Information Technology （in Chinese），2010，32（11）：2760-2763.

[10]王惠濂，李大心，鄧世坤等.探地雷達在建筑地基探測中的應(yīng)用.地球科學(xué) （中國地質(zhì)大學(xué)學(xué)報），1993，18（3）：323-328.Wang H L，Li D X，Deng S K，et al.Application of Ground Penetrating Radar to construction foundation exploration.Earth Science（Journal of China University of Geosciences）（in Chinese），1993，18（3）：323-328.

[11]孫洪星，康永華.有耗介質(zhì)探地雷達波傳播衰減特性的研究.工程地質(zhì)學(xué)報，1999，7（4）：344-348.Sun H X，Kang Y H.Study on GPR wave propagation attenuation behaviors in heterogeneous medium.Journal of Engineering Geology （in Chinese），1999，7（4）：344-348.

[12]底青云，許琨，王妙月.衰減雷達波有限元偏移.地球物理學(xué)報，2000，43（2）：257-262.Di Q Y，Xu K，Wang M Y.The attenuated radar wave migration with finite element method.Chinese Journal of Geophysics （in Chinese），2000，43（2）：257-262.

[13]Irving J D，Knight R J.Removal of wavelet dispersion from ground-penetrating radar data.Geophysics，2003，68（3）：960-970.

[14]劉四新，曾昭發(fā).頻散介質(zhì)中地質(zhì)雷達波傳播的數(shù)值模擬.地球物理學(xué)報，2007，50（1）：320-326.Liu S X，Zeng Z F.Numerical simulation for Ground Penetrating Radar wave propagation in the dispersive medium.Chinese Journal of Geophysics（in Chinese），2007，50（1）：320-326.

[15]John H B.Frequency-dependent attenuation analysis of ground-penetrating radar data.Geophysics，2007，72（3）：7-16.

[16]Gabor D.Theory of communication.Journal of the IEEE，1946，93：429-457.

[17]Gressman A，Morlet J.Decomposition of hardy functions into square integrable wavelets of constant shape.SIAM J.Math.Anal.，1984，15（4）：723-736.

[18]Stockwell R G，Mansinha L，Lowe R P.Localization of the complex spectrum：the S transform.IEEE Transactions on Signal Processing，1996，44（4）：998-1001.

[19]高靜懷，陳文超，李幼銘等.廣義S變換與薄互層地震響應(yīng)分析.地球物理學(xué)報，2003，46（4）：526-532.Gao J H，Chen W C，Li Y M，et al.Generalized S transform and seismic response analysis of thin interbeds.Chinese Journal of Geophysics（in Chinese），2003，46（4）：526-532.

[20]劉喜武，年靜波，劉洪.基于廣義S變換的吸收衰減補償方法.石油物探，2006，45（1）：9-14.Liu X W，Nian J B，Liu H.Generalized S-transform based compensation for stratigraphic absorption of seismic attenuation.Geophysical Prospecting for Petroleum（in Chinese），2006，45（1）：9-14.

[21]Pinnegar C R，Mansinha L.Time-local Fourier analysis with a scalable，phase-modulated analyzing function：the S-transform with a complex window.Signal Processing，2004，84（7）：1167-1176.

[22]Widess M B.How thin is a thin bed.Geophysics，1973，38（6）：1176-1180.

[23]Knapp R W.Vertical resolution of thick beds，thin beds and thin-bed cyclothems.Geophysics，1990，55（9）：1183-1190.

[24]Guha S，Kruse S E，Wright E E，et al.Spectral analysis of ground penetrating radar response to thin sedimentary layers.Geophysical Research Letters，2005，32：L23304，doi：10.1029／2005GL023933.

[25]Guha S，Kruse S，Wang P.Joint time-frequency analysis of GPR data over layered sequences.The Leading Edge，2008，27（11）：1454-1460.

[26]劉喜武，劉洪，李幼銘等.基于廣義S變換研究地震地層特征.地球物理學(xué)進展，2006，21（2）：440-451.Liu X W，Liu H，Li Y M，et al.Study on characteristics of seismic stratigraphy by generalized S-transform.Progress in Geophysics（in Chinese），2006，21（2）：440-451.

[27]陳學(xué)華，賀振華，黃德濟.時頻域高分辨地震層序識別.吉林大學(xué)學(xué)報 （地球科學(xué)版），2008，38（1）：152-155.Chen X H，He Z H，Huang D J.Seismic sequence identifying with high resolution in time-frequency domain.Journal of Jilin University （Earth Science Edition） （in Chinese），2008，38（1）：152-155.

[28]Pinnegar C R，Mansinha L.The S-transform with windows of arbitrary and varying shape.Geophysics，2003，68（1）：381-385.