一种基于改进f-k 滤波的编码雷达信号去噪方法

卢丹平 沈绍祥 李玉喜 赵博 方广有

（1.中国科学院空天信息创新研究院, 北京 100094；2.中国科学院电磁辐射与探测技术重点实验室, 北京 100190；3.中国科学院大学电子电气与通信工程学院, 北京 100049）

0 引 言

超宽带雷达以其高距离分辨率、强穿透力在地质勘探、道路质量检测、穿墙成像、军事等领域得到广泛应用[1-2].探地雷达是超宽带雷达在非导电分层媒质下目标探测与定位的一种综合应用，探测深度与分辨率是其性能的两大衡量指标[3].

由于地下介质的衰减特性，探地雷达回波信号具有低信噪比和大动态范围的特点，这就要求脉冲压缩信号具有低旁瓣特性以避免深层目标的微弱信号淹没在强信号的旁瓣里[4].伪随机编码信号具有自相关特性好、抗干扰能力强以及大时宽带宽积的特点[5]，所以伪随机编码体制的探地雷达能够较好地解决上述问题.但编码雷达在探测过程中仍不可避免地受到噪声的干扰，这对回波数据的解释带来一定困难.雷达回波信号中的噪声包括相干噪声和非相干噪声，相干噪声主要来自剧烈变化的电性界面的散射和绕射，非相干噪声主要是探地雷达系统内部的随机噪声[6].将频率-波数(f-k)滤波应用于探地雷达，有针对性地设计滤波器对噪声有良好的滤除效果，并且能保留并突出有效回波.

滤波是去除干扰的主要手段，利用信号与噪声频谱分布的差异，去除一定频率范围的信号，从而达到去噪的目的[7].探地雷达图像降噪算法包括傅里叶变换滤波、奇异值分解(singular value decomposition，SVD)滤波[8]、小波变换滤波[9]以及f-k滤波等.本文主要针对f-k滤波开展研究，f-k滤波是在频率-波数域(f-k域)中进行的二维过程，通过沿空间采样入射的视速度差异实现f-k谱信号区与干扰区的分离.2008 年，Aziz 等人[10]提出了一种基于f-k滤波器的高性能地震仪降噪算法，通过并行处理提高了f-k算法的执行效率和性能.2009 年，Hayashi 等人[11]设计了两种f-k滤波器抑制双基地探地雷达的直达波，根据发射天线和接收天线的位置信息设计滤波器，直达波分量的抑制比降到-40 dB 以下，实现了对目标的可靠性探测.2012 年，Che 等人[12]提出了一种基于改进矩阵铅笔算法和f-k滤波的去噪方法，采用改进的矩阵铅笔法构造合适的f-k滤波器参数，有效抑制了噪声信号.2016 年，Liu 等人[13]设计了一种用于光纤水听器拖拽阵列的f-k变换降噪方法，结合多尺度滤波和拉东变换，信噪比提高了4～9 dB，结果表明该方法可作为波束形成、目标识别、目标跟踪的预处理方法.2021 年，Yang 等人[14]提出了一种新的fk联合带通角加权模板，利用角信息与频率信息，逐步提高图像的对比度，减弱了噪声对图像质量的影响.f-k滤波可以有效去除噪声信号，但前人所提fk滤波算法无法有效滤除信号高频段的噪声，对探地雷达的数据解释带来困难.

本文提出一种分段式f-k滤波方法，根据接收到信号的能量分布情况，将频段划分为高频段和低频段，分别对高低频段设置不同视速度阈值的f-k滤波器，通过提高高频段的视速度阈值达到抑制高频段噪声的目的.通过仿真实验验证了该方法的有效性，并使用该方法进行了雷达实测数据实验验证，即先对采集的数据进行去除直流、自动增益等预处理；再分别采用传统f-k滤波方法与分段式f-k滤波方法去噪；最后通过脉冲压缩得到雷达图像后，对比原图与两种方法去噪后的图像熵[15]大小，判断噪声滤除效果.结果表明，提出的分段式f-k滤波方法的去噪效果优于传统方法.

1 f-k 滤波方法

1.1 f-k 滤波理论

f-k滤波本质是一种二维信号变换[16].雷达回波是一种时间和空间的二维函数，假设为d(t,x)，其中t为时间变量，x为空间变量[17].对雷达回波d(t,x)做二维傅里叶变换，即可得d(t,x)在频率-波数域的频谱信息，其二维傅里叶变换及对应的二维傅里叶反变换为[18]

式中：f为频率；k为波数；D(f,k)为d(t,x)的f-k谱.

离散的二维傅里叶正、反变换可表示为

设f-k滤波的输出结果为y(t,x) ，y(t,x)由雷达回波d(t,x)与 时间-空间特性h(t,x)进行二维褶积运算得到，二维滤波方程可表示为

式中，Y(f,k)、H(f,k)由y(t,x)、h(t,x)二维傅里叶变换得到，而h(t,x)和H(f,k)决定了滤波器的性质.

在f-k域中，雷达信号d(t,x)沿 测线以视速度v传播.针对侧面墙体或物体，探测过程中的视速度模型如图1 所示.

图1 雷达回波视速度模型Fig.1 Radar echo apparent velocity model

假设收发天线间的距离为S，发射天线发出的电磁波经媒质表面到达目标，目标反射后电磁波到达接收天线的时间为t1；电磁波经侧面墙体或物体反射到接收天线的时间为t2.当收发天线中心位于x0处时，经地下目标反射的电磁波所用时间t1为

式中：tc为电磁波经媒质表面反射回接收天线的时间；v′为 电磁波在媒质中的传播速度；x为天线位于目1标正上方的位置；h为目标位置深度.天线的高度通常为厘米级，传播延迟可以忽略不计，则式(4)空气中的传播时间tc可忽略.

经侧面物体反射的电磁波所用时间t2为

式中，c为空气中电磁波波速.

则目标和侧面物体的视速度v1和v2分别为：

如图2 所示，f-k域图像中纵坐标代表频率f，横坐标代表波数k，故斜率的大小代表视速度v.因此，滤波器设计可以直接根据视速度的不同，分离有效信号与干扰信号[21]，即保留视速度高于v∗的有效信号区，滤除视速度低于v∗的混叠区和干扰区.

图2 f-k 滤波器设计Fig.2 Design of f-k filter

1.2 分段式 f-k滤波方法

在传统f-k滤波过程中，滤波器呈上宽下窄的扇形结构，高频段f-k谱的保留范围大于低频段.但一般情况下，能量集中于低频段，若采用传统方法，高低频段视速度阈值相同，视速度大于阈值的高频段噪声无法被滤除.此外，一般情况下，有效信号相对于全波数域的波数分布范围较窄且高频段能量低，若采用上述滤波器，在保证低频段信号被有效保留的情况下，其高频段仍存在一部分未被滤除的噪声.基于此提出一种分频段工作的f-k滤波器，设定的高频段视速度阈值大于低频段视速度值阈值，改进的滤波器模型如图3 所示，仅保留绿色的有效信号区.在保留低频段内有效信息的同时，去除高频段的干扰信息.

图3 分段式f-k 滤波器设计Fig.3 Design of sub-band f-k filter

分段式f-k滤波的步骤具体可分为： 1)获取探地雷达回波数据，进行二维傅里叶变换得到f-k谱；2)根据数据的f-k谱，确定高低频段的视速度阈值和，再根据不同阈值设计滤波器，得到处理后的fk谱；3)进行二维傅里叶反变换，得到滤波后的回波数据[22].其中，步骤2 最关键，滤波器设计得合理与否将直接影响成像质量.

2 数值模拟

为研究分段式f-k滤波方法对噪声的滤除效果，建立数值模型开展仿真实验.

模型为长5 m、深2.5 m 的土壤层，土壤的相对介电常数 εr1≈3；土壤深1 m 处分别倾斜和水平放置一块长0.5 m、厚0.06 m 的大理石板，其水平中心位置为(1.8 m，1 m)和(2.6 m,1 m)，相对介电常数均为εr2≈6.

发射天线中心频率设为800 MHz，发射天线和接收天线间距为0.4 m.采样频率设为5.12 GHz，图4(a)所示为该模型经GprMax3.0 生成的B-Scan 图，两条双曲线表示两块大理石板产生的有效回波信号.图4(b)所示为加入高斯白噪声后的B-Scan 图，加入噪声后的回波峰值信噪比为13.1 dB，有效回波信号被噪声严重干扰，无法准确分辨目标物体.

图4 B-Scan 图Fig.4 B-Scan image

对加噪声后的回波信号进行f-k变换，通过设计合理的滤波器滤除噪声，加入噪声后回波的f-k谱如图5 所示，颜色越亮，表示该区域能量越强.将信号的f-k谱往频率轴和波数轴投影，分析能量在f-k域分别随频率和波数的变化，结果如图6 所示.可以看出：在频率域，能量集中在-1.28～1.28 GHz 的频段内；在波数域，能量集中在-15～15 m-1范围内.根据图6(a)能量随频率的分布情况划分频段，|f|≤1.28 GHz 为低频段，1.28 GHz<|f|<2.56 GHz 为高频段.

图5 加入噪声后回波的f-k 谱Fig.5 The f-k spectrum of the echo after adding noise

图6 不同域能量分布情况Fig.6 The energy distribution

为了进一步分析信号能量随波数分布的情况，更合理地设计滤波器，图7 给出了低频段与高频段信号能量随波数的变化情况.在有效信号区内，低频段信号能量的波数分布在-15～15 m-1，高频段信号能量的波数分布在-8～8 m-1，其波数分布范围小于低频段，且有效信号区内低频段的能量高于高频段.由于两块目标大理石板放置距离较近带来的干扰，高频段信号能量在波数-20～-15 m-1以及15～20 m-1范围有两处尖峰.图中能量较前两者低的其余区域为噪声干扰区，此区域中高频段的能量高于低频段.

图7 不同频段能量随波数分布情况Fig.7 Distribution of energy with wavenumber in low and high frequency bands

综上，与低频段相比，高频段在有效信号区范围内的能量较弱，但混叠区和噪声干扰区的能量却较强，所以滤除高频段的噪声是很有必要的.

滤波器的设计原则是在保留有效信号的基础上尽可能多地滤除噪声.如果视速度阈值v∗选择过高，则低频段有效信号的一部分会被滤除；如果视速度阈值v∗选择过低，则高频段的很大部分噪声无法被滤除.传统方法在较完整保留低频段有效信号的情况下，高频段的噪声滤除效果欠佳.基于此，本文提出一种分频段工作的f-k滤波器.

根据能量分布情况和式(8)，将传统f-k滤波方法的视速度阈值v∗设为0.09 m/ns，将分段式f-k滤波方法的低频段视速度阈值设为0.09 m/ns，高频段视速度阈值设为0.32 m/ns.高于视速度阈值的信号区fk谱保留，低于视速度阈值的信号区f-k谱置零.基于传统f-k滤波方法和分段式f-k滤波方法处理的fk谱如图8(a)和(b)所示.

图8 传统与分段式方法处理仿真数据得到的f-k 谱及B-Scan 图Fig.8 The f-k spectrum and B-Scan image obtained after processing the simulation data by traditional and frequency-division method

将传统f-k方法与分段式f-k方法滤波后的f-k谱通过离散二维傅里叶反变换得到去噪后的BScan 图，结果如图8 (c)和(d)所示.与加入噪声后的图像相比，传统f-k方法处理后的图像已经有了很大改善，基本可以观测到两个目标，经该方法处理后的峰值信噪比为25.3 dB；但分段式f-k方法对图像细节的改善效果更好，经该方法处理后的峰值信噪比为33.2 dB.

根据仿真实验结果，传统f-k滤波方法对探地雷达图像有良好的去噪效果，经传统方法处理后的峰值信噪比提升了12.2 dB.分段式f-k滤波方法在传统方法的基础上去除了高频段的噪声，经该方法处理后的峰值信噪比较传统方法处理后又提升7.9 dB，进一步优化了去噪效果.

3 实验验证与效果分析

3.1 数据采集与预处理

为验证数值模拟结果与实际雷达探测结果的一致性，开展实测实验.伪随机编码超宽带雷达采用格雷互补码，码元平衡直发的信号源设计方案[23].通过对格雷码元信号的每一位进行曼彻斯特编码[24]获得中心频率为0.8 GHz 的码元信号，发射信号频率为0.1～1.4 GHz.发射端功率设置为27 dBm，收发天线采用Botiew 天线的HH 极化形式，采样频率为5.12 GHz.该信号的时频域波形如图9 所示.

图9 发射信号时频域图Fig.9 Time-frequency domain diagram of the transmitted signal

实景测试实验在一个人工测试池进行，该池大小为7 m×3 m×2.5 m，填充材料的相对介电常数εr1≈3.在测试池底部放置金属板，并于池深1 m 处分别倾斜和水平放置一块花岗岩板，两块花岗岩板间隔0.3 m，厚度均为0.06 m，花岗岩板相对介电常数εr2≈6.实验预设的具体模型如图10 所示.

图10 预埋目标示意图Fig.10 Schematic diagram of pre-buried targets

工程模型雷达支持的垂直分辨率理论最高为5 cm，在目标埋深位置支持的水平分辨率[25]为23 cm.将其放置于测试池上方，从北向南以缓慢步行的速度移动采集数据，空间采样率为0.01 m，实验场景如图11 所示.通过沙坑实验，获得700 条通道数据.接收条件的差异导致雷达记录道与道之间的能量不均衡，影响剖面上同相轴的连续性.为改善剖面质量，须进行平均道处理.雷达回波通常浅层能量强，深层能量弱，需要对回波进行自动增益处理.

图11 实验现场照片Fig.11 Photo of the experiment site

3.2 f-k滤波处理及效果分析

将预处理后的数据进行二维傅里叶变换，二维谱分布于波数-200～200 m-1、频率-2.56～2.56 GHz，放大后得到的f-k谱如图12 所示.

图12 预处理后回波信号f-k 谱Fig.12 The f-k spectrum of echo signal after preprocessing

在频率域，能量集中于-1.28～1.28 GHz 的低频段内；在波数域，能量集中于-12～12 m-1范围内.根据高低频段能量分布情况和式(8)，将传统f-k滤波方法的视速度阈值设为0.11 m/ns，分段式f-k滤波方法的低频段视速度阈值v设为0.11 m/ns，高频段视速度阈值设为0.42 m/ns.基于传统f-k滤波方法和分段式f-k滤波方法处理的f-k谱如图13(a) 和(b)所示.

图13 传统与分段式方法处理实测数据后得到的f-k 谱Fig.13 The f-k spectrum obtained after processing the measured data by traditional and frequency-division method

分段式f-k滤波方法滤除的高频段噪声区域大于传统方法，在保留低频段内有效信息的同时进一步去除了高频段的干扰信息.采用以上两种滤波器进行对比实验，滤波后通过脉冲压缩成像，并分析滤波效果， 进行区域截取与深度反演后，原始图像与两种不同f-k滤波器处理后的图像如图14 所示.图14(a)原始图像可以清晰显示放置于1 m 处的异常反射体，即预先放置的花岗岩板.雷达回波中的噪声掩盖了一部分有效信号，使得系统实际水平分辨率降低，在水平距离2～4 m 内只能观察到一个异常反射体，无法区分两块花岗岩板.此外，测试池的底部金属板可以被清晰观察到.

图14 原始图像和两种方法进行区域截取与深度反演后结果对比Fig.14 Image comparison before and after processing by 2 methods

经f-k滤波后的图像如图14(b)所示，在有效回波基本不受影响的情况下，滤除了大部分噪声，实现了两个目标的初步分离.分段式f-k滤波方法进一步滤除了回波信号中高频段的噪声，减少了干扰信息，图像更平滑，对于目标的观测更有利.如图14(c)所示，两个目标在图像中的分离情况更彻底，分段式设计的方法进一步提高了水平分辨率，且能够一定程度上滤除目标物体下方沙粒产生的绕射波.

探地雷达主要使用信噪比来评价噪声抑制效果，该方法需要结合目标特征区域的判断，且目标特征区域的选择受经验和选择方法影响较大[11].本文选取图像熵的评价方法来判断噪声抑制情况.在雷达B-scan 中，图像熵表示图像的混乱程度，图像熵越大，目标信号与背景信号差异越小，目标提取越困难；图像熵越小，目标信号与背景信号差异越大，目标提取越容易.根据文献[14]中提到的方法进行图像熵的计算，即：

式中：m为通道数；n为时间采样点数.

图像熵的计算结果如表1 所示.

表1 图像熵计算结果Tab.1 The result of calculating the image entropy

结合表1 和图13 可以看出，图像经f-k滤波处理后减小了8.6%的图像熵，提高了信噪比，增大了目标信号与背景信号的差异，使得两个目标能被分开观测到.改进的分段式f-k滤波方法较原图减小了14.6%的图像熵，图像熵在传统方法基础上进一步减小6%，进一步增大了目标信号与背景信号的差异，使得两个目标的分离情况更乐观.

4 结 论

本文针对传统f-k滤波对高频段的噪声滤除效果欠佳的问题，提出一种分段式f-k滤波方法，对高低频段设置不同的视速度阈值，优化高频段的滤波效果.根据伪随机编码体制下探地雷达实测数据的处理结果，该方法在传统方法基础上进一步提高了图像的质量及水平分辨率，有利于目标的观测.采用图像熵的杂波抑制评价标准，结果显示经该方法处理后的图像熵较原图减小了14.6%，较原方法减小了6%，有效提高了信噪比.