吳海波，張平松，胡雄武，秦 鎮(zhèn)

(1.中國礦業(yè)大學(xué)(北京)煤炭資源與安全開采國家重點(diǎn)實(shí)驗(yàn)室，北京 100083;2.安徽理工大學(xué)地球與環(huán)境學(xué)院，安徽 淮南 232001)

探地雷達(dá)作為一種新興的工程物探手段，現(xiàn)廣泛用于城市工程、水利以及交通工程檢測(cè)中的多個(gè)環(huán)節(jié)，包括道路路基檢測(cè)、地下管線探測(cè)、隧道超前探測(cè)與襯砌檢測(cè)，以及土體、壩體檢測(cè)等。但在探地雷達(dá)數(shù)據(jù)采集過程中，受到環(huán)境、人為以及儀器等因素制約，探地雷達(dá)數(shù)據(jù)不可避免的會(huì)受到噪聲干擾。噪聲的存在不僅限制了有效信號(hào)的精確識(shí)別，同時(shí)會(huì)產(chǎn)生“虛假”異常，對(duì)探地雷達(dá)的數(shù)據(jù)處理和解釋工作造成嚴(yán)重的影響[1-3]。

目前，頻域?yàn)V波是去除探地雷達(dá)數(shù)據(jù)中噪聲干擾的有效方法[4]。帶通濾波器、巴特沃斯濾波器等常被用于探地雷達(dá)數(shù)據(jù)的濾波處理[5]。除此之外，一些非線性的時(shí)頻分析及變換方法，如小波閾值、經(jīng)驗(yàn)?zāi)７纸?、主成分分析?奇異值分解等也在探地雷達(dá)數(shù)據(jù)處理中得以應(yīng)用[6-9]。但總體而言，現(xiàn)階段多數(shù)濾波方法仍較多的關(guān)注單一頻率因素，造成部分情況下濾波的效果不理想。為了克服單一因素濾波的局限性，部分學(xué)者嘗試將地震數(shù)據(jù)勘探中的f-k濾波方法應(yīng)用到探地雷達(dá)數(shù)據(jù)處理中來，但這部分研究更多的是關(guān)注直達(dá)波與有效回波的分解問題[10-12]。

為此，本文將重點(diǎn)聚焦于探地雷達(dá)數(shù)據(jù)中噪聲的濾除，尋求通過分析含噪聲正演記錄與實(shí)采數(shù)據(jù)的f-k譜特征，并設(shè)計(jì)具有針對(duì)性的扇形濾波器，通過f-k濾波實(shí)現(xiàn)探地雷達(dá)數(shù)據(jù)中噪聲的濾除，保留并突出被噪聲“掩蓋”的有效回波。

1 f-k濾波方法

1.1 二維傅里葉變換

f-k濾波的理論基礎(chǔ)為二維傅里葉變換。設(shè)雷達(dá)信號(hào)為y(t，x)。t為時(shí)間變量，而x為空間變量，則y(t)表示單道記錄；y(t，x)表示多道記錄。定義y(t，x)的二維正反傅里葉變換分別為

(1)

式中：Y(f，k)稱為y(t，x)的頻率—波數(shù)譜(f-k譜)，相應(yīng)的變換稱為f-k變換。

通常采集的探地雷達(dá)數(shù)據(jù)為離散數(shù)據(jù)，則離散二維正反傅里葉變換可表示為

(2)

1.2 濾波器設(shè)計(jì)

f-k域的二維濾波方程表示為

(3)

式中：Y(f，k)可由二維傅氏變換得到，H(f，k)為設(shè)計(jì)的濾波器。

探地雷達(dá)反射信號(hào)與地震反射波信號(hào)在波的類型和頻率上存在顯著的差異。但針對(duì)f-k濾波器設(shè)計(jì)，只需考慮兩種信號(hào)在頻率方面的差異，即地震信號(hào)的頻段通常為100Hz以內(nèi)，而探地雷達(dá)信號(hào)的頻率達(dá)50MHz以上。因此，根據(jù)探地雷達(dá)數(shù)據(jù)中有效回波和噪聲在f-k平面上的分布特征，設(shè)計(jì)扇形濾波器進(jìn)行濾波如圖1所示。

圖1 扇形濾波器示例

設(shè)計(jì)的濾波器還必須滿足如下條件

(4)

其中，有效區(qū)為圖1中灰色部分，及有效回波主要的f-k譜能量分布區(qū)。

1.3 探地雷達(dá)數(shù)據(jù)的f-k濾波流程

依據(jù)離散傅里葉變換的原理和扇形濾波器設(shè)計(jì)的方法，本文建立探地雷達(dá)數(shù)據(jù)f-k濾波的流程如圖2所示，具體可分為三步：

1)讀取探地雷達(dá)剖面數(shù)據(jù)，利用二維離散傅里葉變換計(jì)算數(shù)據(jù)的f-k譜；

2)根據(jù)數(shù)據(jù)的f-k譜特征，設(shè)計(jì)合適的扇形濾波器，并與探地雷達(dá)數(shù)據(jù)的f-k譜相乘，得到濾波后的f-k譜；

3)利用二維離散傅里葉反變換將濾波后的f-k譜轉(zhuǎn)換成濾波后的探地雷達(dá)剖面輸出顯示。

圖2 探地雷達(dá)數(shù)據(jù)f-k濾波的流程

2 濾波效果分析

2.1 含噪聲正演記錄的濾波效果分析

圖3(a)所示為常見探地雷達(dá)數(shù)值模型的正演記錄，圖3(b)為加入高斯白噪聲后的正演記錄(信噪比為10，峰值信噪比為11.77)，所示噪聲對(duì)有效回波產(chǎn)生了干擾，“掩蓋”了部分有效回波的細(xì)節(jié)，影響有效回波的準(zhǔn)確識(shí)別，如圖3(b)箭頭所示。

圖4所示為正演記錄f-k譜的正周期和主周期，可以看出，有效回波的f-k譜能量主要集中在f軸附近的扇形區(qū)；圖5所示含噪聲記錄的f-k譜中，有效回波的能量部分被噪聲的能量“掩蓋”；正演記錄f-k譜原先沒有能量的位置受到噪聲的干擾，也存在一定值的能量。

(a)原始正演記錄 (b)含噪聲正演記錄圖3 正演記錄與含噪聲正演記錄

(a)正周期 (b)主周期圖5 含噪聲正演記錄的f-k譜

依據(jù)圖1設(shè)計(jì)扇形濾波器，從含噪聲記錄的f-k譜著手進(jìn)行濾波，濾波后的f-k譜如圖6所示，扇形濾波窗以外的f-k譜能量為0，與原紀(jì)錄一致，扇形濾波窗以內(nèi)的位置多為有效回波的f-k譜能量，得以保留。通過離散的二維傅里葉反變換得到濾波后的記錄如圖7所示，圖中的記錄相對(duì)于圖3(b)記錄在細(xì)節(jié)上有了很大的改善(如圖7中箭頭所示)，濾波后的記錄峰值信噪比為22.84，相對(duì)于含噪聲正演記錄的峰值信噪比顯著提高。

圖6 濾波后的f-k譜

圖7 f-k濾波后的正演記錄

2.2 實(shí)采探地雷達(dá)剖面的濾波效果分析

在某高校的地質(zhì)工程實(shí)驗(yàn)場(chǎng)地采集的探地雷達(dá)剖面經(jīng)過球面擴(kuò)散補(bǔ)償后如圖8(a)所示，在擴(kuò)散補(bǔ)償?shù)倪^程中有效回波得到了增強(qiáng)，但噪聲干擾也隨之增強(qiáng)，造成圖像下部的有效回波信號(hào)被噪聲部分“掩蓋”，有效信號(hào)識(shí)別難度加大，很難基于該記錄進(jìn)行有效的解釋。不僅如此，從對(duì)應(yīng)的f-k譜(見圖8(b))同樣可以看出， 在頻率0～2GHz、 波數(shù)0～120的范圍內(nèi)噪聲的f-k能量與有效回波的能量部分混疊在一起。

(a)實(shí)采探地雷達(dá)剖面

(b)f-k譜圖8 實(shí)采探地雷達(dá)剖面及其f-k譜

設(shè)計(jì)如圖1所示的扇形濾波器進(jìn)行實(shí)采數(shù)據(jù)的f-k濾波，濾波后的探地雷達(dá)剖面及其與原記錄的殘差分別如圖9和圖10所示。圖9中，濾波后剖面上部的有效回波基本不受影響，殘差普遍小于0.2，而剖面中下部原先受噪聲影響嚴(yán)重的位置，噪聲得以較好的濾除，殘差普遍在0.5左右，部分有效回波得以突出，能有效的進(jìn)行識(shí)別，如圖9中箭頭所示。

圖9 f-k濾波處理后的探地雷達(dá)剖面

圖10 f-k濾波后的殘差

3 結(jié)論

為了減小噪聲對(duì)探地雷達(dá)數(shù)據(jù)中有效回波的干擾，提高有效回波的識(shí)別與解釋精度，本文借助于地震數(shù)據(jù)處理中的f-k濾波方法進(jìn)行探地雷達(dá)數(shù)據(jù)的濾波處理，得到如下結(jié)論：

(1)探地雷達(dá)有效回波的f-k譜能量集中在頻率軸附近的扇形區(qū)域，借助于扇形濾波器能有效濾除正演記錄中的噪聲；利用含噪聲正演記錄測(cè)試表明，經(jīng)過f-k濾波后，地震記錄的峰值信噪比從11.77增加到22.84。

(2)實(shí)采探地雷達(dá)剖面的f-k譜中有效回波能量與噪聲能量存在部分混疊， 通過設(shè)計(jì)針對(duì)性的扇形濾波器， 并對(duì)比分析濾波后的剖面與殘差， 發(fā)現(xiàn)f-k濾波方法對(duì)記錄剖面中下部的噪聲能有效的濾除，殘差值達(dá)0.5左右；剖面上部有效回波基本不受影響，而原先下中部被噪聲“掩蓋”的有效回波能得以突出。因此，可以認(rèn)為f-k濾波方法在壓制探地雷達(dá)數(shù)據(jù)中的噪聲，提高資料信噪比方面效果顯著。