基于主成份分析的探地雷达信号去噪研究

2019-11-19谭方玉

物探化探计算技术 2019年5期

刘雷, 谭方玉

(1.重庆市勘测院，重庆 401121； 2.贵州省有色金属和核工业地质勘查局物化探总队，都匀 558000)

0 前言

探地雷达(Ground penetrating radar, GPR)是利用高频宽带电磁(1Mhz-10GHz)波，对地下结构和特性或者物体内部不可见目标体进行探测定位的一种浅地表地球物理方法。凭借其具有操作简单、分辨率高、抗干扰能力强、适应性强、解释成果准确及效率高等优点，在交通建设、地质勘探、水利工程、城市建设、地质灾害、农业地质探测、考古、隧道、军事等许多领域都表现出广阔的应用前景[1]。在实际探测过程中，浅部地下介质复杂多变，具有频散特性，电磁波在其中传播衰减速度快，散射非常强烈，且存在各种人为设施的杂乱回波以及仪器本身或系统噪声等多种因素的影响，使得得到的GPR剖面中含有各种各样的杂波和噪声，降低了雷达信号的质量。因此，研究高信噪比的探地雷达数据去噪方法显得尤为重要[2]。实际噪声在实践中是不可避免的，必须在预处理数据的同时对噪声处理。对GPR数据进行去噪已经得到了广泛地研究，这些方法涉及将数据转换为不同的域，然后利用波场在域内的特性进行去噪(如傅立叶变换、加窗傅立叶变换[3]、小波变换[4]、小波-KL变换[5]、S变换[6]、希尔伯特-黄变换[7]、经验模态分解[8]以及独立成份分析[9]等)。这些方法也存在着各自难以克服的缺点。以小波变换为例，近年来小波变换被广泛应用于信号处理领域，其中以小波阈值去噪应用最为广泛，虽然这种方法具有多分辨率的特点，但其不足的地方在于具体小波阈值函数、软硬阈值、尺度系数的选择对去噪结果影响较大，同时对于非连续信号，采用小波阈值去噪后容易产生伪吉布斯现象[8]，使得小波变换在探地雷达数据处理的应用和推广中受到了限制。因此，有必要引入一种快速便捷的方法来去除探地雷达信号中的噪声干扰，以提高雷达数据的信噪比。

地下介质通常具有一定规律，实测的GPR数据中由异常体或者界面所引起的数据信号一般具有相关性，而有噪声所对应的成分无明显相关性[10]。因此可以找到一个域使得其相关性降序排列，达到信噪分离的目的。主成分分析[11](Principal component analysis, PCA)是一种常用的多元统计数据分析方法，该方法能够将原始数据中相关性较强的多个变量转化成含有原有变量最大信息的彼此独立的少数几个变量。PCA 方法可以将数据中最“主要”的成分和结构找出来，去除噪声和冗余成分，使原始复杂数据的维数降低，揭示出在复杂数据下所隐藏的简单结构，采用这些主要成分重构信号，可以达到去噪的效果。

笔者将主成分分析用于探地雷达信号的去噪处理中，为探地雷达信号去噪提供一种快速便捷的方法，以提高实际生产中GPR数据的信噪比。采用PCA确定二维剖面数据中最大变化的模式，将PCA分析之后的特征值降序排列，采用L曲线信噪分离准则选取K值，并采用前K个主成份将信号重构，在保留原始信号的同时消除大部分噪声。

1 主成分分析(PCA)

主成分分析是一种统计方法，通过正交变换将一组可能存在相关性的变量转换为一组线性不相关的变量，转换后的这组变量叫主成分，其常用于降维和特征提[12-14]。

(1)

其中

(2)

(3)

令Sx为(xi)的协方差矩阵。因此tr(Sy) =tr(ASxAT)通过使用拉格朗日乘数并取导数，可以得到

Sxuk=λkuk

(4)

这意味着uk是Sx的特征向量。此时xi可以表示为：

(5)

因此xi可以近似为：

(6)

PCA的算法步骤：

1) 将原始数据按列组成N行M列矩阵x。

2) 将x的每一行进行零均值化，即减去这一行的均值。

3) 求出协方差矩阵Sx。

4) 求出协方差矩阵的特征值λk及对应的特征向量uk。

5) 将特征向量按对应特征值大小从上到下按行排列成矩阵，取前K行组成矩阵A。

6)y=Ax即为降维到K维后的数据。

2 GRP数据PCA去噪流程

实际测量中，噪声的类型是多种多样的，除开地层中不同介质产生的背景干扰噪声，还有仪器本底噪声，直流漂移的干扰和各种外部因素干扰[15]。从野外采集到的探地雷达原始数据，其中既包含有用信息，也包含各种干扰噪声，有时有用信息会被噪声所掩盖，此时需要经过数据处理，才能得到有助于解释的数据和图像。需要压制噪声，增强信号，提高资料信噪比，以便从数据中提取速度、振幅、频率、相位等特征信息，帮助解释人员对资料进行地质解释。

地下介质中会包含多种的天然或人为的杂质，导致接收的回波中包含很多干扰信号，这些干扰统称为背景随机噪声。探地雷达勘探中,记录雷达波时, 为了保持更多的波的特征, 通常采用宽频带进行记录。因此, 在宽频带范围内记录了各种反射波的同时，也记录了各种噪声。假设探地雷达信号是长度为N的时间域离散信号s(n)，被噪声e(n)污染，所测得的含噪数据表示为[16]式(7)。

f(n)=s(n)+e(n)

(7)

去噪的主要任务是尽可能地将实际信号与噪声信号分离开。通常出现的噪声都服从或近似服从高斯分布，因此假设e(n)为高斯噪声。一般利用数学变换解决信号去噪问题，才能够时域转换到频域加以解决，如傅里叶数字滤波和小波变换。笔者研究了PCA在时间域和频率域下去除噪声的有效性：在时间域直接将数据进行PCA去噪；对于频率域而言，首先将每道单道数据进行的离散傅立叶变换，对非负频率数据进行PCA去噪，负频率数据可以通过恢复的非负频率片的复共轭得到，然后使用逆离散傅立叶变换将完整频率数据变换到时域与原始信号进行比较。时间域与频率域的实验流程如下。

1)时域处理流程:①在模拟数据上加入高斯白噪声;②在时间域数据上进行PCA去噪;③对降噪质量进行评估对比。

2)频域处理流程:①在模拟数据上加入高斯白噪声;②对时间维度进行离散傅里叶变换;③在非负频率数据上执行PCA;④对时间维度进行逆离散傅里叶变换;⑤对降噪质量进行评估对比。

将二维时间域和频率域数据进行PCA去噪，其特征值表示为每个相应的特征向量的方差，使用前K个特征向量重构数据，重构数据中只包含最大变化模式的数据，达到保留原始信号的同时消除大部分噪声。而不同K值的选取会对去噪效果产生影响。为了计算最佳模式数K，笔者采用均方根误差(RMSE)以及信噪比(SNR)两个性能指标对降噪质量进行评估对比。

均方根误差(RMSE)：去噪信号与原信号的均方误差为式(8)。

(8)

信噪比(SNR)：原始信号能量与噪声能量的比值为式(9)。

(9)

信噪比越高，则去噪效果越好。

3 合成GRP数据PCA去噪实验

为了研究PCA去噪效果，分析最佳K值的取值方法，首先对合成的理论GPR记录进行处理。此处建立了两个数值模拟模型数据。

建立图1所示的层状模型。模拟区域为1.0 m×1.0 m，模型相对介电常数从上至下依次为3、5、8、5，电导率依次为3 mS/m、5 mS/m、8 mS/m、5 mS/m。模型离散网格为201×201，网格间距为0.005 m，CPML吸收层为10层。将主频为900 MHz布莱克曼-哈里斯脉冲置于地表，模拟时间间隔为0.01 ns，模拟时窗为25 ns，采样时间间隔为0.04 ns。我们采用gprMax软件[17]对该模型进行正演数值模拟，该GPR剖面共181道，道间距为0.005 m，收发距为0.1 m。

图2显示了模拟的GPR剖面和加入高斯白噪声的GPR剖面。采用PCA去噪方法，分别在时间域和频率域对噪声数据进行去噪处理，为了计算PCA去噪中的最佳数量K，在去噪过程中保留了重构信号和原始信号之间的均方误差(RMSE)和信噪比(SNR)用于衡量恢复信号的准确性(图3)。

图1 层状模型Fig.1 Layered model(a)相对介电常数;(b)导电率

图2 高斯白噪声GPR剖面图Fig.2 GPR profile of white gaussian noise(a)层状模型的原始剖面图;(b)含噪剖面

图3 不同域去噪处理后均方误差(RMSE)和信噪比(SNR)对比图Fig.3 Comparison diagram of mean square error (RMSE) and signal to noise ratio (SNR) after denoising in different domains(a)RMSE随 k值的变化图；(b)SNR随 k值的变化图

由图3可以发现，两种不同域去噪处理均存在一个合适的K值使得均方根误差最小且信噪比最大。选择低于最佳值的K会使得信号失真，而选择高于最佳值的K会使得恢复信号中含有更多的噪音。由于合成噪声具有随机频率成分，在时间和频率范围内的去噪效果都相似的。时域和频域中均方根误差与信噪比与K值的关系满足相同的规律，RMSE分别在K=5和K=4取得最小值。

主成分分析去噪可以看做是将观测到的信号划分成信号子空间和噪声子空间两部分。信号子空间是其中的主要成分，可以由前n个最大的特征值对应的特征向量计算得到。相应的，其余的特征值对应的是噪声子空间，也可以由特征值对应的特征向量计算得到。这样得到的信号空间中可以看成是不含噪声的。主成分分析的关键问题是如何确定需要保留的主成分的个数K。在主成分分析中，确定合理的K值是一个比较难解决的问题，其将信号分成了两个子空间，一般来说，信号的信噪比越低，合理的K值越难确定。

然而对于实际数据而言，我们并不能通过均方根误差曲线以及信噪比曲线来选取最佳的K值。通常在主成分分析的过程中，选取主要成分时使用的方法是计算各成分贡献率以及前个成分的累计贡献率[13]。只要大致规定一个累计贡献率的百分比(一般为 85%以上即可)便可以决定选择几个主成分。

这里我们采用L曲线法[18]来进行最佳K值的选取。在基于L曲线的信噪分离准则[19]中，对特征值进行降序排列，以特征值索引作为横坐标，特征值作为纵坐标。那么所获得的特征值曲线可以近似为L曲线。一般来讲，含噪数据的L曲线可明显分为两个部分：①对应有效信号成分的特征值幅度较大，衰减较快；②而对应噪声成分的特征值幅度要小得多，衰减较慢且数值变化比较平稳。根据这个特性，我们可以大致确定 L曲线的拐点位置。将拐点位置对应的K值作为保留成分的依据，该值表示作为PCA分析中的信噪分离的基准点。

图4 不同域特征值的变化曲线对比图Fig.4 Comparison graph of change curves of eigenvalues in different domains(a) 时间域特征值变化曲线；(b)频率域特征值变化曲线

图5 RMSE最小值点与L曲线拐点对应K值的去噪结果Fig.5 Denoising results of the minimum value point of RMSE and the corresponding value of K at the inflection point of L curve(a) 时间域K=5时PCA去噪结果；(b) 时间域K=8(L曲线法)时PCA去噪结果；(c)频率域K=4时PCA去噪结果；(d) 频率域K=7(L曲线法)时PCA去噪结果

图4为时间域和频率域中特征值的变化曲线，其拐点分别为K=8和K=7。图5为RMSE最小值点以及L曲线拐点对应K值的去噪结果，由图5可以发现，L曲线拐点处K值的去噪效果与RMSE最小点K值的去噪效果基本相当，仅在细微地方有些差别，均能在保留原始信号的基础上去除大部分噪声。

为了进一观察基于 L 曲线的信噪分离准则的PCA去噪效果，图6显示了第90道的原始记录、噪声记录以及不同域下不同K值的PCA去噪结果。可以发现RMSE对应的K值在3个反射波处(10 ns ～15 ns)出现了一些波形失真，而基于L曲线信噪分离准则选取K值PCA去噪曲线此处波形与原始信号基本相同

从上述分析可以看到，基于L曲线的信噪分离准则PCA去噪方法可以在保持有用数据不丢失，原幅值不失真的基础上，使噪声得到了较好地抑制，提高数据的信噪比。

4 实测数据处理

图7(a)为探测某隧道地下管线的实测探地雷达剖面，检测过程中采用美国GSSI公司生产的SIR-3000型探地雷达仪进行检测，按照探测目的及要求，选用900 MHz天线。

图6 第90道的原始记录、噪声记录以及不同域下不同K值的PCA去噪结果Fig.6 The 90th data of original records and noise records and PCA denoising results with different K values in different domains(a)原始单道记录；(b)含噪单道记录；(c)时间域K=5时PCA去噪结果；(d)时间域K=8(L曲线法)时PCA去噪结果；(e)频率域K=4时PCA去噪结果；(f)频率域K=7(L曲线法)时PCA去噪结果

图7 实测探地雷达剖面Fig.7 Measured GPR profile of actual data(a)实测探地雷达信号剖面; (b)PCA去噪剖面

图8 实测数据PCA变换特征值变化曲线Fig.8 Actual data PCA eigenvalues value change curve

图9 实测数据第100道数据与PCA其消噪图Fig.9 The 100th data before and after PCA denoising of actual data

共采集400道数据，每道512个采样点，剖面水平距离为4 m，时窗长度为20 ns。采用之前讨论的L曲线信噪分离准则的PCA去噪方法在时间域对该实测数据进行去噪处理。图8为原始数据PCA的归一化特征值曲线，根据L曲线信噪分离原则，其拐点位于K=27的位置。图7(b)为取前27个主成份重构得到的PCA去噪数据剖面图，图中去除了深部(采样点250～500)的雪花状的随机噪声并使得有效信息得到保留，深部数据剖面图变得清晰，特别是位于第125道400采样点的位置的双曲线反射异常以及280道250采样点位置的弱双曲线反射异常，进一步证了该方法可行性。图9为第125道的单道数据去噪前后的对比，从图9可以看出，该方法在有效抑制噪声的同时，能大幅地保留原始信号信息。本文方法能够很好地去除高频噪声，并且可以很好应用于实际生产中。