利用改進(jìn)的f-k濾波方法壓制探地雷達(dá)中的線性干擾

張力文，聶俊麗

（貴州大學(xué)國(guó)土資源部喀斯特環(huán)境與地質(zhì)災(zāi)害重點(diǎn)實(shí)驗(yàn)室，貴州貴陽(yáng) 550025）

1 引言

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

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

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

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

2 基本方法原理

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

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

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

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

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

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

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

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

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

2.2 改進(jìn)的f-k濾波方法

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

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

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

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

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

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

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

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

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

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

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

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

圖2 數(shù)值模型

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

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

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

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

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

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

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

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

圖4 巖溶洼地示意圖

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

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

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

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

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

5 結(jié)論

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