張力文，聶俊麗

（貴州大學國土資源部喀斯特環(huán)境與地質(zhì)災害重點實驗室，貴州貴陽 550025）

1 引言

探地雷達（Ground Penetrating Radar，GPR）以其快速無損、探測精度高等特點，被廣泛用于工程檢測［1，2］、環(huán)境保護［3，4］、文物考古［5～7］、地質(zhì)災害［8，9］、反恐安檢［10，11］、水文水利［12］等領(lǐng)域，是淺表地球物理探測技術(shù)中重要的手段. 在探測過程中，探地雷達會不斷接收到周圍物體界面反射的電磁波信號，這些信號不僅包括天線下方的反射信號，還包括天線四周的反射信號，四周的反射信號會對雷達解釋造成一定的干擾，尤其是由墻體、傾斜地層等面狀界面產(chǎn)生的空間干擾，這些干擾在雷達剖面表現(xiàn)為具有一定傾向的線性干擾，其能量強、范圍廣. 因此，壓制這種線性干擾對探地雷達數(shù)據(jù)的解釋至關(guān)重要.

針對數(shù)據(jù)中的線性信號，地震數(shù)據(jù)處理中常采用f-k變換進行處理.f-k變換實質(zhì)上是一種二維傅里葉變換，主要是利用時間-空間域和頻率-波數(shù)域2 個域之間的關(guān)系［13］. 如提取頻散曲線等面波分析方法依據(jù)面波在f-k域上振幅能量最大的特點，在對應的f-v域中提取與頻率對應的相速度值，從而得到面波記錄中的頻散曲線［14，15］；波場分離方法依據(jù)有效波與干擾波的視速度差異，利用f-k變換分離不同波組［16，17］. VSP 數(shù)據(jù)處理中，Treitel 等采用f-k扇形濾波器，保留濾波范圍內(nèi)的信號，區(qū)分上行信號和下行信號［18］，但存在混波、振幅畸變等缺陷［19～21］. 對此，朱海龍等人［22］提出了一種高保真的VSP波場分離方法以降低混波、振幅畸變，但需要對實際采集的數(shù)據(jù)進行偏移、坐標旋轉(zhuǎn)等偏振濾波處理得到虛擬的上行波射線坐標，將復雜波場轉(zhuǎn)換成簡單波場. 此外，由于f-k變換方法屬于多道處理的范疇，其處理結(jié)果受空間采樣點的分布、深度采樣間隔等采集參數(shù)的影響，使得f-k變換方法在使用中有著很大的限制［23］.

鑒于探地雷達的原理和地震勘探相似，且電磁波類似于彈性波，因此一些學者嘗試將f-k濾波方法應用于探地雷達數(shù)據(jù)處理中.Hayashi等人［24］利用f-k濾波壓制數(shù)據(jù)中的直達波；余海忠等人［25］利用f-k濾波壓制地面干擾波；翟波等人［26］介紹了f-k濾波的方法原理與常用的f-k濾波器. 在實際探測中，探地雷達仍有別于地震勘探，如地震勘探采集方式一般為一發(fā)多收的多偏移距方式，探地雷達則采用自發(fā)自收的共偏移距方式.采集方式的不同使得探地雷達剖面與地震剖面存在差異，且不具有時-空變特性，因而f-k變換方法在探地雷達中的使用限制小. 但采用常規(guī)f-k濾波方法壓制線性干擾時，仍會存在濾波效果不理想的問題.

為此，本文提出了一種改進的f-k濾波方法壓制探地雷達中的線性干擾. 本文將通過對模擬數(shù)據(jù)與實測數(shù)據(jù)中不同類型的線性干擾在f-k域的分布及頻譜特征進行分析，針對不同區(qū)域內(nèi)線性干擾的頻譜特征分別設(shè)計特定的濾波范圍進行濾波，達到壓制數(shù)據(jù)中線性干擾的目的，并與常規(guī)f-k濾波方法進行比較，為工程應用提供一定的參考.

2 基本方法原理

2.1 常規(guī)的f-k濾波方法

儀器采集的數(shù)據(jù)是時間和空間的二維函數(shù)d(t,x)，通過二維正、反傅里葉變換可得到其頻率-波數(shù)（f-k）域頻譜信息和時-空（t-x）域記錄為

其中，f表示頻率；k表示波數(shù)；D(f,k)為d(t,x)的頻率-波數(shù)譜. 式（1）表明d(t,x)是由不同頻率和波數(shù)的平面簡諧波組成，頻率和波數(shù)的關(guān)系為

其中，由于儀器總是沿測線進行探測的，故v為波的視速度. 而在f-k域圖像中，縱坐標代表頻率f，橫坐標代表波數(shù)k，斜率代表視速度v，故具有一定視速度的波組在f-k域頻譜特征的表現(xiàn)呈一條譜值射線，并且視速度越大的波組其譜斜率越大. 因此根據(jù)波組的視速度設(shè)計的合適濾波器進行處理可達到壓制干擾波的目的.

f-k濾波器的性質(zhì)由時間-空間特性h(t,x)或頻率-波數(shù)特性H(f,k)所確定. 設(shè)y(t,x)為d(t,x)的f-k濾波輸出結(jié)果，時間域中f-k濾波結(jié)果由輸入信號d(t,x)與濾波因子h(t,x)的二維褶積運算實現(xiàn)；頻率域中則由輸入信號的頻率-波數(shù)譜D(f,k)與濾波器的頻率-波數(shù)特性H(f,k)相乘完成，即

其中，Y(f,k)和H(f,k)是y(t,x)和h(t,x)的傅里葉變換結(jié)果.

由于儀器采集的數(shù)據(jù)是時間序列長度有限的N道離散記錄，故離散化后的f-k濾波公式為

其中，n為原始道號；m為結(jié)果道號. 由式（4）可知，測線上任一點處，f-k濾波的結(jié)果可由N道探地雷達數(shù)據(jù)道通過一維濾波結(jié)果相加而得. 可見，f-k濾波本身也是一種空間域的二維濾波.

在地震勘探數(shù)據(jù)處理中，利用實信號傅里葉變換結(jié)果具有對稱性的性質(zhì)，即負波數(shù)信號是正波數(shù)信號的復共軛，f-k濾波方法只對采集數(shù)據(jù)的f-k域正波數(shù)信號進行處理后，再通過對稱性即可得到處理后的全頻域信號.

2.2 改進的f-k濾波方法

實際探測中，探地雷達常用采集方式為自發(fā)自收的共偏移距方式，其偏移距可視為零，則反射時距曲線關(guān)系式為

其中，V為波速；H是天線到面狀界面的垂直距離. 采集過程中，隨著天線的移動，H會發(fā)生呈正或負的線性變化，其反射時間V也隨之呈正或負的線性變化. 因此，探地雷達剖面中會出現(xiàn)傾向不同、視速度方向不同的2種線性干擾（見圖1），其中向下傾的線性干擾（下行干擾）具有正視速度，向上傾的線性干擾（上行干擾）具有負視速度，其視速度v為

圖1 探地雷達剖面中的傾斜線性干擾

其中，l為天線移動距離；Δt為相同波組的反射時間變量. 聯(lián)合式（2）可知，探地雷達剖面中具有正、負視速度方向的線性干擾經(jīng)f-k變換后，會位于f-k域的不同區(qū)域中，且根據(jù)其視速度的不同，所代表的f-k域頻譜特征也不同.

常規(guī)的f-k濾波方法在處理過程中未考慮到探地雷達實際數(shù)據(jù)中線性干擾具有正、負視速度，且視速度值不同，導致其在f-k域分布及頻譜上具有差異，在濾波時僅限于處理正波數(shù)域內(nèi)的信號，因而對線性干擾的壓制效果不理想.

為此，本文對f-k濾波方法進行改進，在濾波過程中對f-k域進行全域分析，對實際數(shù)據(jù)在f-k域中正、負波數(shù)域內(nèi)的信號均進行處理，而不僅限于分析處理正波數(shù)信號. 具體方法步驟如下：

（1）采集探地雷達數(shù)據(jù)經(jīng)過傅里葉變換后，得到其f-k域頻譜；

（2）對f-k域全頻域內(nèi)的頻譜特征進行分析，根據(jù)正、負速度線性干擾的f-k頻譜特征及分布的差異，在正、負波數(shù)域內(nèi)分別設(shè)計不同濾波區(qū)域，框選正、負波數(shù)域內(nèi)線性干擾所在的區(qū)域；

（3）將選定的濾波區(qū)域內(nèi)譜值置零，得到處理后的f-k譜；

（4）將處理后的f-k譜通過二維傅里葉反變換轉(zhuǎn)換到時空域輸出顯示.

3 數(shù)值模擬及f-k濾波結(jié)果對比

本研究建立了如圖2所示的數(shù)值模型，模型大小為1.4 m×1.1 m，網(wǎng)格大小為0.002 m×0.002 m×0.002 m. 模型分為兩層介質(zhì)：上面是1.4 m×0.1 m的空氣層，相對介電常數(shù)為1；下面是粘土層，范圍1.4 m×1.0 m，相對介電常數(shù)為9. 在粘土層中部有一埋深為0.66 m、半徑為0.04 m 的圓形空洞，其中心位置為（0.7 m，0.3 m）. 模型兩端對稱放置了兩塊0.2 m×1.0 m三角金屬物體，傾斜面角度約為78.7°. 雷達天線的激勵源采用點電偶極子源，激勵波形為Ricker子波，頻率為800 MHz，收發(fā)天線相距0.1 m，天線的移動步長為0.002 m. 為反映側(cè)面物體在探地雷達剖面中產(chǎn)生的干擾情況，將測線的起點布置于左側(cè)金屬物體右側(cè)（如圖2所示），與左側(cè)金屬物底端相距0.02 m，測線終點布置于右側(cè)金屬物體左側(cè)，如圖2所示，與右側(cè)金屬物底端同樣相距0.02 m，共采集430道數(shù)據(jù). 模型中雷達測線的正下方只存在圓形空洞這一異常體.

圖2 數(shù)值模型

將模型的各參數(shù)編寫入“.in”文件中，采用GprMax3.0軟件進行數(shù)值模擬計算，得到如圖3（a）所示的探地雷達剖面圖. 在圖3（a）中，雙曲線信號為圓形空洞在雷達剖面產(chǎn)生的有效信號，與之相交的上、下行干擾則是由兩側(cè)金屬物體產(chǎn)生的空間干擾波. 其中，下行干擾由左側(cè)金屬物體產(chǎn)生，上行干擾由右側(cè)金屬物體產(chǎn)生，且兩者視速度值大小相等，約為0.05 m/ns；雙曲線信號則不具有固定的視速度，其視速度隨曲線切線改變而改變，但總體上視速度值均大于上、下行干擾的視速度值.

圖3（b）是對數(shù)值模擬數(shù)據(jù)進行f-k變換后得到的f-k域頻譜特征圖，明顯可見正、負波數(shù)域內(nèi)頻譜特征一致. 在圖3（b）中，上、下行干擾在f-k域中均表現(xiàn)為以坐標原點為起點的譜值射線，其中下行干擾的f-k頻譜特征位于正波數(shù)域內(nèi)，上位行干擾的f-k頻譜特征則于負波數(shù)域內(nèi)，兩者關(guān)于f軸對稱. 同時，根據(jù)上、下行干擾的頻譜射線得到兩者的視速度值均在0.05 m/ns 左右，這與雷達剖面的視速度值基本一致. 而雙曲線信號由于視速度不定，且視速度值均大于上、下行干擾，因而f-k域中正、負波數(shù)域內(nèi)均有分布，且譜值特征更靠近f軸并關(guān)于f軸對稱.

圖3（c）是利用常規(guī)f-k濾波方法進行處理后得到的雷達剖面，f-k濾波器采用扇形濾波器，選取的濾波參數(shù)為2 000 MHz的頻率上限和0.06 m/ns的視速度下限，主要保留f-k域內(nèi)視速度大于0.06 m/ns 的信號. 經(jīng)濾波處理后得到圖3（c）所示的雷達剖面圖，如圖所示，雷達剖面中的上、下行干擾均得到了有效壓制，雙曲線反射雖得到保留，但曲線兩側(cè)形態(tài)發(fā)生改變且能量有所衰減.

根據(jù)模擬數(shù)據(jù)f-k域頻譜特征的分析，將f-k域內(nèi)雙曲線信號頻譜特征外的上、下行干擾的頻譜特征進行框選，得到了如圖3（b）虛線區(qū)域所示的濾波區(qū)域.由于正、負波數(shù)域內(nèi)的線性干擾頻譜特征相同，因而該濾波區(qū)域關(guān)于f軸對稱，該濾波區(qū)域內(nèi)視速度上限為0.65 m/ns.

圖3（d）是利用改進的f-k濾波方法，采用圖3（b）所示的濾波區(qū)域，對線性干擾的f-k頻譜特征進行處理后得到的雷達剖面圖. 如圖3（d）所示，兩側(cè)金屬物體產(chǎn)生的線性干擾得到了有效壓制，圓形空洞產(chǎn)生的雙曲線仍有清晰的反映. 將圖3（c）與圖3（d）進行對比，發(fā)現(xiàn)常規(guī)的f-k濾波方法和改進的f-k濾波方法均能有效壓制模擬數(shù)據(jù)中的傾斜線性干擾，但改進的f-k濾波方法保留有效信號的效果更好.

圖3 數(shù)值模擬結(jié)果及f-k濾波對比

4 數(shù)值模擬及f-k濾波結(jié)果對比

測區(qū)為貴州某小型封閉巖溶洼地（見圖4），洼地形如“碗”狀，底部呈不規(guī)則的圓形，直徑約46 m，洼地四周側(cè)壁較陡且傾斜角度不一致，目測均大于30°. 洼地內(nèi)表層為粘土層，其下為灰?guī)r層. 采用SIR-20探地雷達及100 MHz 天線在洼地內(nèi)部進行巖溶探測，在所采集到的東西向與南北向雷達剖面中均出現(xiàn)了由洼地四周側(cè)壁產(chǎn)生的強烈空間干擾. 本文選取了其中沿洼地南北向長測線得到的雷達剖面（見圖5）進行研究.

圖4 巖溶洼地示意圖

如圖5（a）所示，洼地四周側(cè)壁產(chǎn)生的空間干擾在雷達剖面表現(xiàn)均為具有不同反射強度、不同傾向、不同視速度的傾斜線性干擾. 在該剖面中，明顯可見一條能量強的下行干擾，大致由雷達剖面（4 m，45 ns）處起始，到剖面中部（26 m，174 ns）處截止，視速度約為0.17 m/ns，推測該干擾是由測線起始處、洼地南面的側(cè)壁產(chǎn)生. 此外，在圖5（a）明顯可見多條不同的上行干擾，其中，剖面底部（17 m，250 ns）和（23 m，210 ns）處有兩條長度短、能量較強的干擾，其視速度相近，約為0.15 m/ns，推測可能是測線17 m、23 m 附近兩側(cè)洼地側(cè)壁產(chǎn)生的；剖面（30 m，160 ns）至（41 m，57 ns）處有一能量強的上行干擾，其視速約為0.11 m/ns，推測可能是由測線終止處、洼地北部側(cè)壁產(chǎn)生；剖面右下部可見能量弱的多條上行干擾，其視速度大體近似，在0.16 m/ns 左右，推測可能是由位于洼地南面的測線兩側(cè)邊坡側(cè)壁產(chǎn)生的.

圖5（b）是對實測數(shù)據(jù)進行f-k變換后得到的f-k域頻譜特征圖，明顯可見正、負波數(shù)域內(nèi)頻譜特征存在差異. 在圖5（b）中，明顯可見兩條以原點為起點、譜值較大的譜值射線，其中f-k域中正波數(shù)域內(nèi)為下行干擾的反映，負波數(shù)域內(nèi)則為能量強的上行干擾在f-k域的反映. 同時，由于下行干擾視速度值大于上行干擾的視速度值，因而下行干擾的譜值射線斜率更大，相對上行干擾的譜值射線斜率更靠近f軸. 圖5（b）中還可見正、負波數(shù)域內(nèi)還存在一些不以坐標原點為起點、譜值較小的譜值射線，為實測數(shù)據(jù)中其他能量較弱的上、下行干擾在f-k域中的反映. 此外，實測數(shù)據(jù)中地下異常體反射的有效信號受各種噪聲的干擾，在f-k域中的表現(xiàn)并不明顯，而線性干擾的頻譜特征相對更突出. 在巖溶洼地中，溶洞、破碎帶等異常體的反射信號的視速度值通常大于線性干擾的視速度值，因此在f-k域中的頻譜特征相對線性干擾的頻譜特征會更靠近f軸.

圖5（c）是利用常規(guī)f-k濾波進行處理后的雷達剖面，f-k濾波器采用扇形濾波器，濾波參數(shù)中頻率上限為200 MHz，視速度下限為0.4 m/ns，主要保留f-k域內(nèi)視速度大于0.4 m/ns的信號. 經(jīng)濾波處理后，剖面中線性干擾的能量得到了衰減，可見剖面中部出現(xiàn)溶洞反映，但剖面中的線性干擾仍清晰可見，且在剖面（15 m，110 ns）和（37 m，100 ns）范圍產(chǎn)生了新的干擾.

根據(jù)實測數(shù)據(jù)f-k域頻譜特征的分析，在保留f-k域中f軸附近頻譜信號的情況下，對f-k域中上、下行干擾代表的頻譜特征進行框選，選擇了如圖5（b）所示的濾波區(qū)域. 由于f-k域中正、負波數(shù)域內(nèi)頻譜特征不同，故正、負波數(shù)域內(nèi)的濾波區(qū)域也有所不同，其中正波數(shù)域內(nèi)濾波區(qū)域視速度上限為0.42 m/ns，負波數(shù)域內(nèi)濾波區(qū)域視速度上限為0.32 m/ns.

圖5（d）是利用改進的f-k濾波方法，采用圖5（b）所示的濾波區(qū)域，對線性干擾的f-k頻譜特征進行處理后得到的雷達剖面. 如圖5（d）所示，經(jīng)濾波處理后，雷達剖面中的線性干擾得到了有效壓制，且未產(chǎn)生其他強烈的干擾，中部的溶洞反映也清晰可見. 將圖5（c）與圖5（d）對比可看出，相對于常規(guī)f-k濾波方法，改進的f-k濾波方法對線性干擾的壓制效果更好.

5 結(jié)論

常規(guī)的f-k濾波方法利用FFT 變換的復共軛特征，只對正波數(shù)信號進行處理即可得到全頻域內(nèi)的信號，忽視了實際數(shù)據(jù)中與正波數(shù)信號不同的負波數(shù)信號，因此雖能有效壓制模擬數(shù)據(jù)中的線性干擾，但對實際數(shù)據(jù)中線性干擾的壓制效果不理想. 對此，本文對f-k濾波方法進行改進，提出一種基于f-k域全域頻譜特征的濾波方法. 在對雷達數(shù)據(jù)進行f-k變換后進行f-k全域分析發(fā)現(xiàn)，不同傾向線性干擾其視速度不同，因而在f-k域的頻譜特征及分布區(qū)域存在差異. 基于這種差異，在f-k域正、負波數(shù)域內(nèi)設(shè)計不同的濾波區(qū)域框選出線性干擾所在的頻譜區(qū)域進行濾波，區(qū)別壓制正、負波數(shù)域內(nèi)線性干擾的頻譜，達到壓制雷達數(shù)據(jù)中線性干擾的目的. 通過模擬數(shù)據(jù)及實測數(shù)據(jù)的分析處理，與常規(guī)f-k濾波方法相比，改進的f-k濾波方法對線性干擾的壓制效果更顯著，對探地雷達剖面的處理效果更好.