基于MUSIC谱估计的探地雷达回波处理方法

(长安大学 信息工程学院，陕西 西安 710064)

探地雷达(Ground Penetrating Radar，GPR，也称地质雷达)系统是一种通过频段位于106～109Hz的无线电波来探测地下结构的无损探测仪器。探地雷达可以探测金属及非金属埋藏物，现在已经广泛应用于考古、矿产资源勘探、灾害地质勘察、路面结构检测等众多领域。在实际应用中，探地雷达接收到的波形常常会产生畸变且信号的信噪比较弱[1]。而由于被测物体具有强反射性、被测媒质内的大尺寸颗粒反射及不均匀的物体电导率、雷达设备天线间的耦合等因素导致探地雷达接收到的信号中含有大量的噪声、直达波等各种散射体的杂波分量。这些信号由于传播路径短且衰减少，因此在到达时间上相互重叠，导致实际所测的目标信号常常被淹没，因此需要对采集到的探地雷达回波信号进行滤波处理，提高所测数据的精度[2-5]。

探地雷达回波信号的处理通常从时域和频域两个方面进行：① 时域处理方法，对接收到的信号直接提取并分离出有效的速度、振幅等信号值,例如数字滤器的处理方法；② 频域处理方法，对接收到的雷达回波信号进行频谱分析。经典常用的算法有小波分析法、FFT算法等。其中，FFT算法对数据样本的长度要求较高，而小波分析在特定情况下(已知噪声的频率范围且信号和噪声的频带相互分离时)非常有效。但在实际应用中对广泛存在的白噪声处理效果较差[5-8]。因此，笔者在雷达信号预处理后引入谱分析中信号子空间与噪声子空间的概念，将MUSIC算法与最小二乘法相结合来进行探地雷达杂波及噪声处理。

1 探地雷达数据处理方法

1.1 直流偏移量处理

在探地雷达系统中，有时雷达剖面上的数据含有直流漂移量，这时的数据全是正的或全是负的，亦或正负半周出现不对称情况。此时需要对数据进行处理，若直接利用接收到的回波信号，肯定会产生误差，这对后续波形的利用产生影响，因此需要消除或压制[9]直流成分。通常直流波去除采用消除水平同相轴的方法，该方法假设直流偏移量是一个常量，因此在探地雷达测量过程中,将每道回波数据求平均值即可得到直流偏移量，再将每条原始回波数据依次减去直流偏移量即可得到移除直流偏移量的回波数据。

假设回波信号为X(n)，消除水平同相轴去除直流偏移的方法可以表示为

(1)

式中，X(n)为未经处理的原始信号;n为采样点数量;X′(n)为去除直流偏移量后得到的回波信号。

1.2 均值法滤波

均值滤波是一种线性处理技术，是一种最简单且最常用的数字滤波方法，对抑制高斯噪声有较好的效果。该方法是在某时刻连续多次进行信号采样，对得到的采样值进行算术平均处理，得到的值将作为此刻的最终信号值，其中连续采样的次数视具体情况而定。

假设在测量过程中探地雷达获得N段采样数据，在此记为xi(i=1,2,…,n)。根据最小二乘法的原理可知：存在一个数值y，使得各采样值间的偏差平方和为最小值，则y就是最接近真实值的数据，则有

(2)

依据一元函数求极值的原理可以得到

(3)

在随机噪声干扰的环境下，每次采样接收的信号均可以表示为

xi(n)=si(n)+ei(n)

(4)

式中，si(n)为目标的回波信号；ei(n)为噪声。

探地雷达对同一目标多次采样时，N次采样值的算术平均值为

sk(n)+ek(n)]

(5)

当采样时间很短时，采样信号可被看作近似相等，因此在N次采样后，得到的信号s(n)的分量可表示为

(6)

(7)

假定采样过程中，噪声环境不发生变化，则噪声ei(n)为相互独立且概率密度相同的随机变量。这样根据噪声ei(n)的性质，经过平均值滤波后的噪声大小为

(8)

1.3 MUSIC功率谱估计

探地雷达接收到的信号包括被测物体回波、噪声和强杂波信号。通常雷达回波可以用时间的实数振荡函数来模拟，而雷达系统中的所接收回波中的强杂波可以近似看成服从高斯分布的有色噪声，在本文中将杂波信号假定为经过滤波后的高斯白噪声[10]。因此，雷达回波信号X(n)可看作K个复正弦信号与高斯白噪声信号的组合，即

(9)

此时，探地雷达信号的滤波问题变为由已测得到的回波信号X(n)恢复出回波信号中的未知量Ai、ωi。下面采用多重信号分类的MUSIC算法来进行ωi估算。

MUSIC(Multiple Signal Classification)谱估计法，又称为多重信号分类算法。该算法是Schmidt于1979年提出的一种高分辨率空间谱频率估计方法[11]。此算法依据矩阵特征值分解理论，将信号分为信号子空间和噪声子空间，利用信号与噪声的正交性，从自相关矩阵求出最小的特征值及相对应的特征向量并构建空间谱函数，再通过谱峰搜索估计出未知的频率。

依据MUSIC算法的工作原理，探地雷达信号X(n)的自相关函数可表示为

(10)

假设接收到的雷达信号可表示为N个采样值的组合，则信号矢量可表示为

ei=[1,ejωi,…，ej(N-1)ωi]T

(11)

此时式(10)就可以写成如下形式：

(12)

式中，I为N×N的单位矩阵。

Rx=Rs+Rw

(13)

将矩阵Rs特征分解得

(14)

式中，λi为特征向量υi所对应的特征值。由于特征向量υi之间是正交的，矩阵Rs的秩为M,所以其特征值中必有N-M个为0。此时式(14)可以写成如下形式:

(15)

对于噪声相关矩阵也可以表示为

(16)

式(10)可以表示为

(17)

(18)

实际做法是构造函数，构造的函数为

(19)

设当ω=ωi时功率谱为无穷大，此时相关矩阵估计值为ωi，为PMUSIC(ω)取到最大值时所对应的值。由于相关矩阵是估计值，必然存在误差，因此PMUSIC(ωi)为有限值，但会呈现出很尖的峰值，此时峰值对应的频率就是复正弦信号的频率。

1.4 计算幅值与相角

信号的最佳频率ωi和正弦信号个数M可以由MUSIC谱方法准确地估算出来，但是通过该方法无法检测出正弦波的幅值和相角，因此需要结合其他方法进行信号参数计算。主要采用求解过程简单、计算量较小的最小二乘法[12]来估算振幅值和相角。

假设已知信号记为X(t)，振幅和初相分别为Ai和φi，最佳组成频率成分为fi，其中fi=2π/ωi,噪声为w(t)。此时信号的数学关系表达式为

(20)

由式(20)可知，雷达信号经过时域n点采样后可以表达为

(21)

通常情况下，当n>>2m时，式(19)为可表示为一组线性超定方程组，将此方程写为矩阵形式：

A=

当噪声较小时，式(21)可写成如下形式：

B=AX

(22)

取n=2m,矩阵A即为常数阵，将式(22)两边同乘以A-1,则可以求解出X，即

X=A-1B

(23)

但在实际应用中常取n>2m,此时矩阵A不为方阵不存在逆矩阵，可以用A的伪逆矩阵A+来求得未知数X的解，式(22)的解可表示成如下形式：

X=A+B

(24)

求出X中的所有元素后，可以用式(25)和式(26)求出振幅和相位：

(25)

(26)

2 仿真与实测数据分析

实验仿真主要分为两个部分，一部分用于验证MUSIC算法的可行性，另一部分为实测探地雷达回波处理过程。

2.1 基于模拟数据的MUSIC算法验证

假设信号的采样率为Fs=1000 Hz，信号长度为N=2048，表达式为

利用MUSIC谱估计方法， 加入0.01倍信号噪声后的测试结果如表1所示。

表1 仿真信号的参数估计值

由表1中列出的测试结果可以发现，原始信号与加入噪声后的仿真结果基本一致，只有幅值和相位结果存在少许差别，这主要是因为噪声的存在使得结果产生了一些误差。而在处理实际测量数据时，会采用第1节中的均值法和去直流偏移量进行回波信号预处理，这样可以有效地降低噪声，而这种近似带来的误差是可以接受的。

2.2 探地雷达实测数据分析

如图1所示，以探地雷达探测层状结构层为例来说明雷达回波数据的处理过程。先采用消除水平同轴法与均值滤波法对回波信号进行数据预处理，再采用MUSIC谱估计与反褶积方法进行去噪声处理，以期达到更好的效果。

图1 探地雷达回波信号处理过程

采用的雷达设备为WGPR系列无线雷达，雷达中心频率为1 GHz。为了尽可能降低外部环境干扰和人为干扰对测量结果的影响，实验中的沙箱尺寸为150 cm×50 cm×50 cm，将探地雷达悬挂在支架上，探地雷达距离沙子的距离为1.45 m，沙箱中沙子厚度为10 cm。图2为探地雷达测量层状结构示意图。

图2 探地雷达测量层状结构示意图

雷达采样后得到的实时图形如图3所示，该段数据为10 cm砂层的雷达剖面图，由450条采样数据构成，采样点为 8192。其中，图3中的右侧波形即为采样的单道回波信号，在实际数据处理时将提取每道回波的采样值进行数据预处理。

图3 厚度为10 cm沙层的采样数据

(1) 直流波去除和均值法滤波。

去直流波后利用均值法滤波将450条雷达回波进行数据处理后，得到结果如图4所示。

图4 均值法滤波

(2) MUSIC谱估计去噪声。

对图4中进行滤波后的信号运用MUSIC法进行功率谱估计，得到图5中的功率谱曲线及图6中的滤波结果，其中红线代表0.2×max(20lg|PMUSIC(ω)|)的阈值线。

(3) 回波信号处理有效性验证。

利用重构波形计算沙层的厚度来验证本文中回波信号处理过程的有效性。

图5 MUSIC方法进行谱估计的功率谱图

图6 MUSIC与均值法滤波比较图

已知结构层厚度的通用公式为

(27)

式中，c为电磁波在真空中的速度；Δt为探地雷达电磁波在该层的传播时间;d为结构层厚度；ε为该结构层材料的介电常数值。

本实验中通过介电常数测量仪来确定沙层的实际介电常数。取多次测量平均值为实际沙层的介电常数，大小为εsand=2.96。

分别对均值法滤波与高斯重构回波数据进行时间计算，计算得出沙层的传播时间Δt约为108 ns，代入式(27)中，通过计算得到砂层厚度为9.493 cm，与沙层厚度实际值的绝对误差为0.507，相对误差为5.07%。而利用本文中方法进行高斯信号回波重构，计算得出雷达波在沙层中的传播时间Δt为113 ns，通过计算得到沙层厚度为9.865 cm，与沙层厚度实际值的绝对误差仅为0.135，相对误差为1.35 %，误差在5%之内，相较于均值法滤波计算结果，精度得到了显著提高，说明本文中的回波处理方法有效，适用于实际测量。

3 结束语

针对探地雷达提出一种简单有效的处理回波信号方法，对雷达接收机接收到的回波信号进行了时域和频域联合处理，先将回波信号进行去直流波、均值滤波等回波后预处理，利用MUSIC谱估计方法估计各信号分量的频率信息，再利用最小二乘法估计其幅值和相位，最终通过反褶积计算进行波形重建，从而达到去除噪声精确估计实用参数的目的。通过理论仿真和实测探地雷达数据处理，表明了该方法能很好地消除雷达回波数据的噪声和杂波，能有效提高测量数据的精度，适合实际应用。

图7 均值法滤波与重构回波时间计算