基于奇異值分解及波數濾波的三維地質雷達噪聲抑制技術研究

徐利劍，胡松濤，2，鄒晨陽，2

（1. 江西省水利科學院，南昌 江西，330029；2. 江西省水工安全技術研究中心，南昌 江西，330029）

0引言

地下管線作為保障城市運行的重要基礎設施，被譽為城市生命線。隨著城市建設不斷進行，密集雜亂的地下管網因年代久遠、施工破壞等因素導致滲漏、淤堵或破裂等問題頻發[1，2]，污染城市環境，擾亂運行秩序。開展地下管線探測工作，對指導地下管線施工、運維和修復等具有重要的科學及現實意義。在眾多探測方法中，地質雷達（GPR）以其快速高效、無損直觀及較高分辨率的特點，在地下管網調查中得到了廣泛關注[3，4]。而三維地質雷達（3D GPR）作為一項近年來發展迅速的高新技術，采用天線陣列收發模式，實現了大量雷達數據的同時采集，覆蓋面更廣，能更為清晰準確的定位地下目標體[5，6]。

在實際應用中，由于城市地球物理環境的復雜多變，加之陣列天線耦合及各種雜亂回波等因素，3D GPR接收信號不可避免地包含了雜波和噪聲，限制了采集數據的質量及后期解譯的精度。如何抑制雜波噪聲干擾，有效提取地下管線回波信號成為重要課題。歷年來國內外學者提出了多種降噪方法，Donoho 等[7]利用基于二階小波的軟閾值法對信號進行優化平滑去噪，該方法簡單有效，但當信噪比較低時，選擇合適的閾值就變得非常困難。傅里葉變換和小波變換亦被廣泛應用于信號降噪[8]，Huang 等[9]將二維連續小波變換引入到GPR 復信號分析中，在強干擾條件下成功提取了異常體信息。Neelamani 等[10]指出小波處理占用內存過多，計算效率低下。他們轉用曲波變換降低了GPR 相干噪聲及背景噪聲，但是該方法只對點狀目標有效。王超等[11]基于希爾伯特-黃變換對信號噪聲進行抑制。Santos 等[12]將時間反轉（TR）技術與各種GPR 預處理方法相結合，將地下目標體與其他雜波噪聲等分離開來。

近年來，基于子空間的經驗模態分解（EMD）[13]，奇異值分解（SVD）[14]、主成分分析（PCA）[15]或獨立成分分析（ICA）[16]等方法，被廣泛應用于GPR 噪聲抑制。其主要思想是將雷達信號分解為目標和噪聲，去除噪聲部分后重構信號以達到突顯目標體的目的。這些方法的主要問題在于實際探測信號中的噪聲子空間難以精確劃分，去除噪聲子空間的同時可能損失目標有效信號，導致重構信號缺失部分重要信息。此外，前人的研究成果基本都集中于二維GPR 信號處理，鮮少涉及對三維GPR 信號的噪聲抑制研究。

本文針對三維地質雷達探測某城區地下管線所得回波深度切片，提出首先基于SVD 對信號進行分解得到近似（Gross）和細節（Detail）兩個部分，對兩部分信號同時進行基于二維高斯低通的波數濾波（Wavenumber filtering）后再重新組合成為重構GPR 信號，實現在保證信號完整性的基礎上對噪聲成分的有效抑制。

1 方法概述

1.1 SVD 基本原理

三維GPR 深度切片C-scan 由m×n 的二維矩陣表示，記為：

式中：xi，j是第i 行第j 列的元素。對A 矩陣做SVD可得：

式中：U 和A 是大小分別為m×m 和n×n 的單位正交矩陣，U 為左奇異矩陣，V 為右奇異矩陣。D 是大小為m×n 的對角矩陣：

式中：σr為矩陣A 的奇異值，并且σ1≥σ2≥…≥σr，r 代表矩陣A 的秩。這樣矩陣A 可以改寫為：

式中：uiviT是大小為的m×n 單秩矩陣，代表了矩陣A 的第i 個特征像元。

由于奇異值是遞減排列的，前幾個奇異值對原始GPR 信號貢獻更大，后面較小的奇異值則可視為不重要的信息或隨機噪聲。在得到特征像元和奇異值后，可利用公式（4）來重構GPR 信號。

1.2 Wavenumber filtering 基本原理

基于高斯低通的波數濾波是一種非線性濾波器，可以達到平滑降噪的效果。其基本原理是：由GPR 陣列天線中不同位置天線接收的回波深度切片均可表達為空間的函數g（x，y），對函數在x，y 空間方向做二維傅里葉變換得到波數域表達：

式中：k 為空間波數；f 為電磁波頻率；v 為電磁波速。

經SVD 算法分離出的噪聲信號和目標信號頻譜成分十分接近甚至重合，但電磁波速度不同則波數特征不同。當進行地下管線探測時，各向異性的土壤等地下介質是噪聲的主要來源，由于土壤中電磁波速一般小于地下管線的反射回波，在介質中就表現為高波數的偽聲，可通過高斯低通濾波器進行壓制。再經二維傅里葉逆變換，即可獲得清晰有效的探測信號。

需要注意的是天線陣列的間距布置△x 必須滿足空間采樣定理，即：

式中：km為信號的最大波數。否則容易出現假頻。

2 工程案例應用

2.1 測區概況

本次野外測量在江西省某探測試驗場區內進行。測區路面為瀝青，下鋪水泥碎石穩定基層，再下即為土壤，成份以粉質粘土、壤土、砂壤土為主，局部夾有薄層粉細砂及雜填土。下伏基巖主要以白奎系鈣質泥巖、泥質粉砂巖、灰礫巖、砂礫巖透鏡體和第三系紫紅色泥巖、泥質粉砂巖為主。

試驗段在深度約1.8m 設三處測點，分別埋設了1根110mm PVC 管，1 根75mm PVC 管，以及石棉瓦及雜物若干。兩處PVC 管均是水平布置，第三處的石棉瓦及雜物為傾斜布置，傾角約20°左右。

2.2 3D GPR 測量

測量使用某公司三維GPR，陣列天線中心頻率為200MHz，天線采用一發雙收，道間距約7.5cm。采集時窗設置為100ns，采樣點數1024。對試驗段全部掃描完成后我們截取了深度1.832m 處的反射回波切片，如圖1 所示。

圖1 3D GPR 深度切片原始圖像

從圖1 可以看出，在橫坐標17～18m 處大致可見PVC 管的反射，在橫坐標21～23m 之間隱約可見另一根PVC 管的反射。而在橫坐標35m 附近的石棉瓦加雜物的反射基本不可見。根本原因還是因為地下介質復雜多變，土壤又具頻散特性，電磁波傳播衰減快，且散射繞射強烈，使得有效反射信號被噪聲等干擾信號掩蓋。

2.3 算法處理應用

為抑制各種噪聲和雜波等對于地下管線反射信號的干擾，有效提取出圖像中目標信息，利用本文上述算法對原始GPR 深度切片進行降噪去雜波處理。

首先對原始圖像做SVD 處理。SVD 方法難點在于：城市環境中GPR 數據噪聲干擾較大，難以選取合適的奇異值區段對信號進行分解，人為選取較為困難且有其主觀性。陳高祥等[17]基于直達波和奇異值個數的線性相關性提出了一種自適應SVD 算法，但只限于2D GPR B-scan 中直達波的去除，且方法過程較為繁雜?；诖?，我們給出了一個更為簡化的奇異值數目閾值定量判據：

式中：r 為原始GPR 數據中奇異值的總個數；k 為選取的奇異值個數。

公式（9）中分子可表征Gross 信號的總能量，分母可表征Detail 信號的總能量。選取不同的k 值，兩者之間比值不同；當比值近似收斂，即可確定k 值。經對比分析計算，本次試驗我們選取了前三奇異值作為閾值點分離出近似（Gross）信號和細節（Detail）信號，分別如圖2b 和圖2c 所示。

圖2 SVD-Wavenumber 濾波算法處理3D GPR 深度切片

從圖2b 可以看出，SVD 算法能較好的將地下目標的反射信號從干擾信號中分離出來，即：橫坐標約18m處應為110mm PVC 管的反射；橫坐標約22m 處應為75mm PVC 管的反射；橫坐標35m 至36m 處應為石棉瓦及雜物的反射，反射回波的形態也較好地反映出該處目標體的傾斜形態，但仍存在部分干擾。從圖2c 可以看出，Detail 信號較為混雜，不能確認除噪聲信號外是否還包含地下目標體的有效信號。

然后對Gross 信號和Detail 信號采取基于高斯低通的波數濾波，處理后所得圖像如圖2d 和圖2e 所示。與圖2b 和圖2c 進行對比后可以發現，波數濾波后的Gross 圖像基本壓制住了干擾信號，更加突出了目標管線的反射信號，源自于地下不均勻土壤等介質的背景噪聲和散射雜波等高波數的偽聲得到了更嚴厲的抑制，使圖像更為清晰平滑。

最后，基于波數濾波后的Gross 信號和Detail 信號，我們根據公式（4）重構了3D GPR 深度切片信號，所得結果見圖2f。與圖2a 相比，重構后的深度切片目標體反射信號更為清晰，可準確獲得地下管線等目標體的形態、橫向位置及分布情況等重要信息，達到了對圖像平滑去噪的目的。

3 討論

為分析本文方法對3D GPR 信號的去噪效果，采用峰值信噪比指標（PSNR）和圖像熵值指標（Q）做定量評價。PSNR 表征圖像最大信號和噪聲之間的比值，由公式（10）計算[18]。PSNR 越大，表明圖像中信息損失越少，目標信號越明顯，處理效果越好。

式中：MSE 為均方誤差；x[m，n]代表經算法處理后的雷達信號；y[m，n]代表原始雷達信號；k 為每像素的比特數。

Q 表征圖像的信息量，由公式（12）計算[19]。Q 值越小，則圖像越清晰有序，噪聲抑制效果越好。

通過matlab 程序分別計算了原始GPR 圖像（圖2a）、SVD 分離出的近似圖像（圖2b）、波數濾波后的近似圖像（圖2d）及最終重構的GPR 圖像（圖2f）的PSNR和Q 值，所得結果記錄于表1 中。

表1 不同圖像PSNR 和Q 值對比

分析表1 中的數據可以發現，經SVD-波數濾波后再重構的信號PSNR 得到明顯提高，Q 值得到降低，說明該方法增強了地下目標體的信號同時完整保存了有效信號。值得注意的是，圖2f 的Q 值比圖2d 的Q 值高，這是因為重構圖像（圖2f）中包含了波數濾波后近似信號（圖2d）和細節信號（圖2e）中的全部有效信息。綜合分析圖像效果和圖像信噪比及熵值，本文所提出的SVD-波數濾波在3D GPR 噪聲抑制上都取得了良好效果。

4 結論

復雜城市環境中3D GPR 實測地下管線所得原始信號存在大量噪聲及雜波干擾，難以有效提取目標信號。本研究提出了基于SVD 和波數濾波相結合的降噪技術，并給出了SVD 選取奇異值閾值的定量判別方法。結合了SVD 對數據體的切割和重構能力及波數濾波對噪聲雜波的去除能力，在壓制噪聲、增強目標信號的同時可減少目標信息的丟失。研究結果表明，該方法可有效提升GPR 信號的信噪比并增加圖像平滑度和清晰度，更加適用于城市環境下地下管線的識別。