0 引 言
通常情况下，地物杂波是指雷达波束在正常传播情况下探测到的地物回波。地物杂波存在于雷达的低仰角扫描数据或雷达较近的范围，对于任一特定仰角，典型的地物杂波污染从一个体扫到下一个体扫很少有变化，并且大多数时间都会出现。这严重影响雷达的数据估算，噪声则会影响弱信号的检测，从而导致参量估算的问题，因此，在信号处理系统中，滤除地物杂波能够更好地进行气象观测。
在实际项目中可选择全程滤波和动态杂波图滤波。通过CMD(Clutter Mitigation Decision)算法识别杂波图，根据杂波图对杂波位置进行滤波，减少计算量，实现实时处理。在实际项目中设计了多个滤波器以适应复杂环境条件。
本文描述了IIR(Infinite Impulse Response) 椭圆滤波器、自适应高斯频域滤波器(GMAP)和双高斯滤波器(BGMAP)的原理，对GMAP滤波器和BGMAP滤波器进行了详细的分析和比较，并使用实际天气雷达数据对IIR椭圆滤波器和GMAP滤波器进行了测试和比较，结果表明：BGMAP滤波器在分析数据时效果比GMAP滤波器效果好，但耗费时间长，不能满足实时处理要求；GMAP滤波器滤波效果优于IIR滤波器，能满足实时处理要求。
1 IIR椭圆滤波器
通常情况下，采用的是IIR椭圆滤波器[1]进行地物滤波，因为不需要很大的存储器和计算量，在工程上比较容易实现。本文选择的是3阶极点和3阶零点的IIR滤波器，滤波器系数按照重复频率 Fs= 300∶50∶2 000等35个频率，凹口宽度按0.2 m/s步进进行滤波器设置，以满足不同重复频率、不同凹口宽度要求的应用需要。
IIR滤波器存在的主要问题是：
1) IIR滤波器有暂态响应时间，滤波器凹口宽度越窄，输出暂态响应时间就越长；反之就比较短。暂态响应时间长意味着在一个波束宽度内需要有比较多的回波脉冲，否则难以达到所设计的稳态响应。
2) 滤波器凹口宽度如果设置凹口过宽，会对气象回波造成损失；如果设置凹口过窄，滤波后地物剩余仍很多。
3) 会对速度较小的气象回波或者折叠到零频附近的气象回波造成影响。
2 CMD算法
GMAP滤波和BGMAP滤波算法拟合需要一定时间，应当配合动态杂波图进行滤波，动态杂波图使用CMD杂波识别方法[2-3]，处理流程如下：
1) 检查信噪比SNR，如果SNR<3 dB,该距离库为噪声，不作任何处理。
2) 计算3个特征量TDBZ，SPIN，CPA，其中TDBZ选择以当前距离库为中心的9个距离库，SPIN选择11个库。
3) CPA 5点中值滤波。
4) 通过隶属函数将3个特征量映射到0-1区间，如式(7)、式(8)、式(9)所示。
5) 使用模糊逻辑组合SPIN和TDBZ。
6) 计算地物概率CP，如式(10)所示。
7) 将地物概率超过0.5的距离库标识为地物。
8) 对地物标识CF进行平滑和填充。
9) 对标识的距离库进行滤波。
10) 对滤波后的距离库重新计算谱数据，包括反射率、速度和谱宽。
反射率纹理(TDBZ)是相邻距离库的反射率的差平方均值，测量相邻距离反射率的变化，反映了反射率的平滑程度。天气信号是大范围的平滑信号，按式(1)计算：
(1)
杂波相位阵列校准值CPA(Clutter Phase Alignment)用于判断由一组雷达脉冲回波的相位(I和Q)的稳定性。CPA定义所有回波IQ时序的矢量和的幅度与所有回波IQ幅度的计算和的比值。
(2)
式中，xi=Ii+jQi。
SPIN反映了反射率因子在径向梯度的变化，表示一定距离内反射率符号变化的频率。当符号相反，dBZi距离点SPIN数值的增加需要满足
sign{dBZi+1-dBZi}=-sign{dBZi-dBZi-1}
(3)
(4)
时：
MSPINi=1
(5)
(6)
式中，spin_thres表示反射率因子门限，典型值为5 dBz。最后，将SPIN数值除以计算单元数，再乘以100。
(7)
(8)
(9)
(10)
3 GMAP滤波算法
当前国内的主流滤地物杂波算法有IIR滤波及自适应谱滤波，这两种方案对零频附近存在天气信号时，都会造成谱矩估算的偏差。采用有限积累点数进行谱矩估计，相当于对信号加矩形窗，当地物杂波强度较强时，矩形窗对旁瓣的抑制小(-13 dB)，从杂波中泄露出的旁瓣功率会把弱的天气信号淹没。
GMAP算法[4-5]是在2004年由SIGMET公司的两位工程师Siggia和Passarelli提出的，该算法有以下优点：
1) 在零频附近存在天气信号与杂波叠加的情况下，可在滤除杂波的基础上，基本恢复天气信号，使得谱矩估计更为准确；
2) 在对地物杂波信号的谱宽估计更为准确，可以更好且自适应地估量出剔除地物杂波时所需凹口的宽度(传统的IIR滤波需手动切换选择凹口大小)；
3) 可以根据不同的杂信比(CSR)，自适应选择所加窗函数，从而能对旁瓣达到很好的抑制，又不至于使得主瓣过宽影响气象目标的谱宽估计。
4 BGMAP滤波算法
BGMAP算法在频域内去除地物回波，其基本假设是气象回波和地物回波的功率谱为高斯分布[5]。雷达信号的功率谱由3个部分组成：气象回波、地物杂波和噪声，其模型表示如下：
S(v)=
(11)
式中：第一项为气象回波模型，Pw为气象回波功率，vrw为平均多普勒速度，σvw为谱宽；第二项为地物杂波模型，Pc为地物杂波功率，vrc为平均多普勒速度，σvc为谱宽，并且σvc<σvw；第三项为噪声模型，Pn为噪声功率，vN为Nyquist速度。
如果假定功率谱的每条谱线是指数分布的，则有
(12)
式中，
表示根据第m根谱线的计算值，Sm表示期望值。由贝叶斯定理[6-7]取
最大值从而求得Sm。求最大值需要求价值函数的最小值，即
(13)
因此，需要求出满足高斯模型的回波功率P、平均多普勒速度vr、谱宽σ使得式(13)最小。
在BGMAP中，使用标准的非线性拟合算法来拟合功率谱。考虑到气象回波频谱不一定都满足高斯分布，因此重建的频谱尽可能保留原来的频谱而只去掉地物的部分。偏振雷达的参量不能够直接从BGMAP的结果中计算出来。重新构建不含地物回波的IQ数据按照式(18)～式(22)进行处理。BGMAP算法流程图[8]如图1所示。
图1 BGMAP流程图
1) 选择窗函数、FFT求频谱以及估算噪声
首先分别对IQ数据加海明窗，经过FFT变换计算功率谱。如式(18)所示，原始频谱的幅度为A(v)、相位为φ(v)。利用Hildebrand和Sekhon于1974年提出的方法计算噪声功率。通过式(15)计算噪声信号均值的平方，式(16)计算信号平方均值减去信号均值的平方。通过式(17)求出噪声比值R，从而确定噪声线，低于噪声线的为噪声功率谱，高于噪声线的为信号和杂波谱。
(14)
var(Sn)=〈Sn〉2
(15)
var(Sn)=〈Sn2〉-〈Sn〉2
(16)
(17)
2) 拟合高斯气象谱和杂波谱
按照双高斯模型去掉原始频谱的地物成分，得到新的幅度A′(v)和相位φ′(v)。当气象回波的速度小于地物回波的速度时，气象回波的相位会被地物回波干扰，导致与相位相关的参量估算问题，使用随机相位来代替被地物干扰的那些相位，如式(21)所示。这种方式可以有效减少地物干扰。当新的幅度和相位确定后，通过傅里叶逆变换(IFFT)得到新的IQ数据。
FFT[s(n)]=Fs(v)=A(v)ejφ(v)
(18)
F′s(v)=A′(v)ejφ′(v)
(19)
(20)
(21)
s′(n)=IFFT[F′s(v)]
(22)
根据双高斯拟合的结果，vc的值对应地物功率小于气象回波和噪声功率之和的位置。
BGMAP的实时性会受到两个问题影响：首先，在收敛之前非线性的递归运算需要大量的循环；其次，如果频谱不能用双高斯模型来表示，递归可能不收敛，使用单高斯气象模型拟合。单高斯模型和双高斯模型的主要区别就在于是否对存在地物干扰的那些频点进行拟合，如式(23)所示：
S′(v)=swexp
(23)
单高斯拟合方法基于S′(v)，该频谱只包含气象回波和噪声成分。
不含地物回波的IQ数据重建过程与双高斯滤波的重建过程类似，只是新的幅度A′(v)按照式(24)计算：
(24)
3) 检验杂信比(CSR)是否满足窗函数标准
如果CSR>40 dB，加布莱克曼窗重复BGMAP滤波器的滤波过程和动态噪声计算；如果CSR>20 dB，加布莱克曼窗重复BGMAP滤波器的滤波过程；如果CSR>25 dB，就取布莱克曼窗的结果；如果CSR<2.5 dB，加矩形窗重复BGMAP滤波器的滤波过程；如果 CSR<1 dB，就取矩形窗的结果；其他情况输出海明窗的结果。
5 滤波效果分析
图2 GMAP分析结果
图3 BGMAP分析结果
分析气象信号，对比GMAP滤波器和BGMAP滤波器效果。在有气象和强地物的信号时，GMAP滤波器识别如图2所示，BGMAP滤波器识别如图3所示。GMAP滤波器能够识别出杂波点并拟合出气象信号，并计算出速度和谱宽。BGMAP滤波器能拟合出杂波信号和气象信号，并计算出速度和谱宽等信息。二者计算出的速度和谱宽相差不大，但时间差200倍左右。
只有气象信号时，GMAP分析结果如图4所示，BGMAP分析结果如图5所示，二者虽然错误地识别了杂波信号，但都能拟合出气象信号，并计算出速度和谱宽，时间同样相差200倍左右。
图4 只有气象信号时GMAP分析结果
图5 只有气象信号时BGMAP分析结果
图6是对GMAP和BGMAP进行1 000个距离库运算统计的时间图，红线显示BGMAP与GMAP耗费时间比值，可看出BGMAP计算花费的时间是GMAP的2个数量级，在工程应用中实时性差。
本文的实际回波数据是从安徽四创电子股份有限公司的多普勒天气雷达获取。
图7显示的是未滤波的多普勒四屏图，显示了滤波前后反射率因子、滤波后速度V和谱宽W。图8显示了使用IIR滤波器后的反射率因子速度谱宽图。从图中可以看出，在零速线上，天气信号与杂波叠加的情况下，IIR滤波气象回波损失较大。图9中GMAP滤波器能识别气象信号和杂波信号，恢复气象信号，滤波效果优于IIR滤波器。
图6 GMAP与BGMAP花费时间统计
图7 未滤波的强度速度谱宽图
图8 IIR滤波的强度速度谱宽图
图9 GMAP算法强度速度谱宽图
6 结束语
本文对IIR滤波器、GMAP滤波器和BGMAP滤波器进行了分析和处理，结果表明：在非实时分析时，BGMAP的效果要优于GMAP，能更好地识别出气象信号和杂波信号；BGMAP滤波器和GMAP滤波器的效果优于IIR滤波器，由于BGMAP滤波器计算时间要比GMAP计算时间多很多，在工程实践中，不适合使用BGMAP，会造成资源不够用。
鉴于本文只是对双高斯地物滤波器算法的初步研究，还未通过长期实际回波验证，在今后的项目中应作进一步的分析和研究。
[1] 张海蓉. IIR数字滤波器的设计与实现[J]. 贵州师范学院学报, 2016, 32(12):8-11.
[2] HUBBERT J C, DIXON M, ELLIS S M. Weather Radar Ground Clutter. Part II: Real-Time Identification and Filtering[J]. Journal of Atmospheric and Oceanic Technology, 2009, 26(7):1181-1197.
[3] 李腾伟. X波段天气雷达地杂波滤波技术[D]. 长沙:国防科学技术大学, 2016.
[4] 何建新,王旭,刘艳. 自适应高斯频域滤波器在天气雷达中的应用[J].气象, 2010, 36(6):117-121.
[5] 孙召平,张持岸,张建云. 一种基于高斯模型的自适应地物杂波滤波器算法[J]. 太赫兹科学与电子信息学报, 2013, 11(2):250-253.
[6] 吴子博. 数学之美：神奇的贝叶斯方法[J]. 中国新通信, 2018, 20(8):194-196.
[7] 胡胜利,王琦敏. 基于属性加权的朴素贝叶斯分类算法改进[J]. 电脑知识与技术, 2017, 13(29):257-258.
[8] LI Yinguang. Bayesian Approaches to Detect and Mitigate Ground Clutter Mixed with Weather Signals[D]. Norman, Oklahoma: The University of Oklahoma, 2013.
