f-k濾波在探地雷達數(shù)據(jù)去噪中的應用

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

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

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

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

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

1 f-k濾波方法

1.1 二維傅里葉變換

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

(1)

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

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

(2)

1.2 濾波器設計

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

(3)

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

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

圖1 扇形濾波器示例

設計的濾波器還必須滿足如下條件

(4)

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

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

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

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

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

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

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

2 濾波效果分析

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

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

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

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

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

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

圖6 濾波后的f-k譜

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

2.2 實采探地雷達剖面的濾波效果分析

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

(a)實采探地雷達剖面

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

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

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

圖10 f-k濾波后的殘差

3 結論

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

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

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