基于奇异值分解及波数滤波的三维地质雷达噪声抑制技术研究

徐利剑，胡松涛，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 信号的信噪比并增加图像平滑度和清晰度，更加适用于城市环境下地下管线的识别。