自适应阈值方法去除风廓线雷达地物杂波

张福贵 范 潇 何建新2)

1)(成都信息工程大学，成都 610225)2)(中国气象局大气探测重点开放实验室，成都 610225)



自适应阈值方法去除风廓线雷达地物杂波

张福贵1)2)*范 潇1)何建新1)2)

1)(成都信息工程大学，成都 610225)2)(中国气象局大气探测重点开放实验室，成都 610225)

风廓线雷达探测过程中电磁波传输会受到各类杂波的干扰，其中，地物是主要来源。从功率谱数据上看，地物杂波主要集中在零频附近，且幅度较高，不加以抑制会影响气象回波的识别。针对目前常用的小波阈值滤波法在处理近零频回波被杂波覆盖时效果不佳的情况，该文结合风廓线雷达特点，提出一种根据小波分解高频系数自适应确定阈值的方法，并通过模拟数据与风廓线雷达实测数据进行检验，结果表明：即便信号靠近零频，且被杂波覆盖，该方法也能快速准确识别信号回波。同时，该算法原理简单、计算量小、易于实现，在实际应用中能够增加谱峰识别准确率，可为改善风廓线雷达产品质量提供参考。

风廓线雷达； 地物杂波；小波分析；自适应阈值方法

引 言

风廓线雷达(wind profile radar, WPR)作为探测大气风场的遥感设备，同时能够观测大气湿度与检测降水[1-4]、探测大气稳定度及湍流强度[5]，相较于传统的测风仪器与测风方法，风廓线雷达无需人员定时操作，能够连续实时对大气风场进行监测，且有较高的时空分辨率[6]。它主要采用相控阵天线技术，通过发射2对正交的斜波束与1个垂直波束对大气脉动湍流进行探测，然后由5个波束合成风的水平及垂直分量，提供局部范围的三维风场信息。

近年来，随着雷达相关技术的日益完善，风廓线雷达数据逐步被应用于气象服务及科学研究中。1999年美国风暴预测中心通过风廓线雷达数据分析，及时提高了奥克拉荷马-堪萨斯龙卷风的预警级别；2001年美国在一次雪暴预报中将快速更新循环(the Rapid Update Cycle, RUC)天气预报系统和风廓线雷达数据相结合，提高了预报质量[7]；2004年Kitamura等[8]利用风廓线雷达数据，在雨滴谱方面进行了研究；2012年Wang等[9]利用塔中油田区域的风廓线雷达，对当地一次降水成因及形成高度进行分析。

虽然风廓线雷达数据被广泛应用，但如何提高风廓线雷达的数据质量仍是人们研究的热点。每部雷达的使用场合都有各自独特的地形与环境，这为风廓线雷达探测带来了干扰，其中,最明显的是地物杂波干扰，地物杂波主要是树木建筑引起的旁瓣回波和雷达发射接收回路泄漏的直流信号，它干扰雷达回波功率谱的谱峰识别，谱峰识别正确与否影响着后期产品的准确性和可靠性。研究人员采用各种手段力求消除或减弱这种影响，硬件方面如添加杂波隔离网[10]，软件方面如引入“直流抑制”算法[11]。本文以小波分析为主要手段，提出一种自适应阈值方法，通过理论分析、模拟数据测试及实际应用论证该方法的合理性。

1 基于小波分析的风廓线雷达去地物杂波方法

1.1 小波去地物杂波可行性分析

大气回波通常是非平稳信号，这就限制了傅里叶分析的作用。短时傅里叶变换，通常被称为Gabor变换，使用矩形窗将非平稳信号时域隔开，窗内函数被视为平稳信号，一定程度上解决了这一问题，但窗口一经固定，信号的时频分辨率也被固定。在实际分析过程中，低频对应时域变化慢，分析需要窄频窗，高频对应时域变化快，需要宽频窗，于是小波变换被广泛应用于分析此类非平稳信号[12]。

风廓线雷达的主要探测目标是晴空湍流，探测过程中，常伴有噪声与杂波的干扰[13]，式(1)给出了风廓线雷达实际回波信号的组成成分[14]：

A(t)=S(t)+C(t)+N(t)。

(1)

式(1)中，A(t)为风廓线雷达探测到的实际回波，S(t)为晴空湍流净回波，C(t)为杂波，N(t)为噪声。杂波又分为地物杂波与间歇性杂波，本文重点讨论地物杂波，它主要出现在零频附近。噪声N(t)造成原始回波上下小幅度波动，S(t)形成的谱峰是真实的大气回波峰，地物杂波C(t)的强度通常比S(t)强几个量级，因此，湍流峰常被杂波峰覆盖。S(t)，C(t)，N(t)在时域有不同表现，这为利用小波分析方法区分三者提供了可能。

在一个很长的相关时间内，地物杂波的变化相较于湍流回波慢一些，即地物杂波的周期较长。噪声的I路信号(inphase signal，同相信号)、Q路信号(quadrature phase signal，90°相移信号)表现为平均值为0的高斯随机变量，不具相关性。湍流也是平均值为0的高斯随机分布的变量，但具有时间相关性。根据相关性特征，进行多次积累即可区分噪声与湍流信号。而引入小波分析后，则能进一步区分湍流信号与地物杂波。

1.2 小波去地物杂波传统方法

定义时域连续信f(t)，式(2)为其的小波变换：

(2)

图1反映了多尺度分析将信号逐级分解的流程，首先信号通过滤波器被分解为高、低频两路信号，其中,G(ω),H(ω) 分别表示低通、高通滤波器[12]，接着对低频部分进一步分解，直至分解到能够识别临近零频的大气回波频率成分。假设将信号进行i尺度分解，则小波系数由低频系数段ai和高频系数段[di,di-1,…,d1]组成。

图1 信号逐级分解流程Fig.1 The flow chart of signal decomposition

Daubechies 20(db20)小波具备正交性、紧支撑性和近似对称性等特征[15]，是去地物杂波中常用的小波基。Jordan等[16]利用db20小波结合地物杂波特征，Ai等[17]选择db20小波并结合信号正则性，均对地物杂波进行分离，取得显著效果。

小波分解后通过设定阈值对小波系数进行调整，最后重构信号[18-19]。地物杂波主要受低频系数ai控制，传统的设定阈值方法是取高频系数d1段中最大模值为阈值，因其不受地物杂波影响。Jordan等[16]收集分析大量数据，设定了干燥模式和潮湿模式，采用后半段小波系数的最大值及其倍数作为阈值。丁敏等[20]根据后半段小波系数求得阈值，重构信号取得了较好的去地物杂波效果。

1.3 自适应阈值去地物杂波方法

传统阈值方法重点关注高频系数的后半段，即使天气情况发生变化，阈值也只是在原有基础上乘相应的系数，这种方法在风速较大、大气湍流谱峰远离零频时能有较好的发挥；当风速较小时，上述方法因忽略了高频系数的前半段，往往会抑制真实大气回波，本文提出一个改进的方法，利用高频系数整段数据，根据小波系数模值大小自适应选择阈值，这种方法在处理一般情况时，与传统方法效果基本一致；但在应对谱峰靠近零频的情况时较传统方法有更好的表现。

自适应阈值方法步骤如下：

①将信号进行i尺度小波变换，得到长度为N的低频系数ak和长度为[Nj,Nj-1,…,N2,N1]的高频系数[dj,sj,dj-1,sj-1,…,d2,s2,d1，s1]，其中，k∈[1,N],sj∈[1,Nj],sj-1∈[1,Nj-1],s2∈[1,N2],s1∈[1,N1]，根据式(3)对低频系数进行去直流操作：

(3)

②计算高频系数各段的模值平均值，确定最大平均值与最小平均值的位置，若最大平均值或最小平均值出现d1,s1段，则取阈值为d1,s1段的模值平均值，否则按式(4)取最大模值平均值法求得阈值λ：

(4)

③按照式(5)求得缩放比例z：

1

(5)

④通过式(6)对ak值进行截取：

(6)

其中，1

完成上述步骤获得滤波后的低频系数，通过小波逆变换重构信号，地物杂波即被抑制。

2 模拟数据测试

为测试本文算法的作用效果，将首先采用模拟数据进行验证。模拟环境如图2所示，模拟数据通过信号源直接输入接收机获取，修改信号频率，可以获得不同条件的模拟回波数据。本文算法在风廓线雷达数据处理阶段进行。同时,时钟信号也由信号源产生，保证整个数据模拟过程的相干性。

本文使用的模拟数据是与风廓线雷达数字中频偏移1.5 Hz的信号和偏移50 Hz的信号，旨在模拟回波靠近零频与远离零频的情况。根据多普勒频移与多普勒速度关系可得1.5 Hz频偏对应0.18 m·s-1，50 Hz频偏对应5.8 m·s-1。

图2 数据模拟环境Fig.2 The figure of data simulation environment

图3a为传统阈值方法处理1.5 Hz信号的小波系数图，图3b为自适应阈值方法处理后的小波系数图，可知高频系数模值最大值在di段，说明零频附近存在信号，确定阈值进行滤波后低频段小波系数与di相当；传统阈值方法将高频系数d1的最大模值作为阈值系数，导致低频段小波系数接近于0，与di段系数相差较大，易错过真实回波。

图4a为1.5 Hz模拟数据的原始功率谱，零频处存在功率较强的尖峰，即地物杂波，它覆盖了真实的回波信号，增加了识别难度。图4b采用零频剔除法，杂波宽度取3，去除中心3点并用临近值插补，该方法在目前的雷达运行中使用较多，图中(图4b中右上角为-1.5～1.5 m·s-1区间的功率谱放大图，图4c和4d与之相同)该方法有较好表现，能够有效识别近零频信号，但该方法需要人为选择合适的杂波宽度，宽度过宽，会忽略近零频信号，宽度过窄，则起不到去除地物杂波的作用；图4c采用了传统阈值方法，处理后零频处出现凹槽，零频两侧由于凹槽形成两个波峰；图4d采用了自适应阈值方法，选择db20小波基进行6尺度分解，6尺度是一个经验值，可以看到，直流信号得到了有效抑制，右侧1.5 Hz 频偏信号回波变得清晰。说明本文提出的自适应阈值方法较传统阈值方法在识别近零频回波信号时有良好的表现，且不需要人为设置固定参数或更改相应模式。

图3 1.5 Hz信号小波系数图 (a)传统阈值方法，(b)自适应阈值方法Fig.3 Wavelet coefficients of 1.5 Hz signal (a)traditional threshold processing，(b)self-adapting threshold processing

图4 1.5 Hz信号功率谱图 (a)原始功率谱，(b)零频中心3点剔除法处理后的功率谱，(c)传统阈值方法处理后的功率谱，(d)自适应阈值方法处理后的功率谱Fig.4 Power spectrum of 1.5 Hz signal (a)original power spectrum,(b)power spectrum after zero-frequency elimination of 3 points，(c)power spectrum after traditional threshold processing，(d)power spectrum after self-adapting threshold processing

续图4

图5a为50 Hz模拟数据的功率谱，图5b为采用自适应阈值方法处理后的功率谱图，图中显示自适应阈值方法在风速较大的情况下,也能够准确抑制零频，突显真实回波信号。

图5 50 Hz模拟数据功率谱图 (a)原始功率谱,(b)自适应阈值方法处理后的功率谱Fig.5 Power spectrum of 50 Hz signal (a)original power spectrum，(b)power spectrum after self-adapting threshold processing

通过模拟数据的测试，验证了本文提出的地物杂波抑制方法能够成功抑制地物杂波、提高气象信号的识别能力，即使在气象信号被地物杂波信号覆盖的情况下，也有明显效果。

3 外场试验

本文探测数据来自TWP3-M型移动式边界层风廓线雷达，该雷达由北京敏视达公司生产，安置在成都信息工程大学(CUIT)大气观测场(30°34′47″N，103°58′48″E，海拔为450 m)。

图6a为2014年5月17日12：50(北京时，下同)垂直波束低模第13个距离库的I路信号和Q路信号数据。图6c为原始I路信号和Q路信号直接进行6尺度小波分解后的小波系数示意图，由图6a和图6c可以看出，该时段数据明显有强地物杂波混入，且地物杂波成分集中在低频部分。图6d是经过自适应阈值方法处理后的小波系数，低频部分的小波系数值减小了1个量级，图6b为其重构的I路信号和Q路信号，杂波得到抑制。

图6 2014年5月17日12：50成都信息工程大学风廓线雷达I路和Q路信号图及小波系数(a)原始I路和Q路信号，(b)重构I路和Q路信号，(c)原始小波系数，(d)自适应阈值方法处理后的小波系数Fig.6 The time series of I component and Q component and wavelet coefficients of WPR at CUIT at 1250 BT 17 May 2014 (a)original time series of I component and Q component，(b)reconstructed time series of I component and Q component，(c)original wavelet coefficients，(d)wavelet coefficients after self-adapting threshold processing

图7a为5月17日12：50未经处理的功率谱图，零频附近出现一个窄尖脉冲，即地物杂波，由该图基本无法识别出真正的大气回波谱峰。图7b采用自适应阈值去除地物杂波方法，可以看出，零频信号得到了明显的抑制，其中，粗实线为高斯拟合的谱数据，回波信号得到较好的识别，信杂比提高了约30 dB。

图7 2014年5月17日12：50成都信息工程大学风廓线雷达分功率谱(a)原始功率谱，(b)自适应阈值处理后的功率谱Fig.7 Power spectrum of WPR at CUIT at 1250 BT 17 May 2014(a)original power spectrum，(b)power spectrum after self-adapting threshold processing

2014年4月25日上午出现大风现象，09：30垂直波束低模第13个距离库的谱数据如图8a所示，根据气象回波特征及当时天气情况，10 m·s-1附近的信号为真实湍流峰，距零速度区较远，零速度信号强度大于气象回波信号，且零频附近杂波信号谱宽较宽。图8b为自适应阈值方法去杂波后的谱信号，粗实线为高斯拟合后的数据，与原始功率谱相比，信杂比提高了35 dB。图8c采用了零频中心3点剔除方法处理，零频信号幅度略有下降，仍存在较强干扰，这是由于该处发生雷达频谱泄露造成的信号零频展宽现象，固定模式的3点剔除法在处理该状况时效果不理想。而本文方法在处理该情形时效果仍然显著。

图8 2014年4月25日09：30成都信息工程大学风廓线雷达功率谱 (a)原始功率谱， (b)自适应阈值方法处理后的功率谱，(c)零频中心3点剔除方法处理后的功率谱图 Fig.8 Power spectrumof WPR at CUIT at 0930 BT 25 Apr 2014 (a)original power spectrum，(b)power spectrum after self-adapting threshold processing，(c)power spectrum after zero-frequency elimination of 3 points

图9a为2014年5月17日12：00低空连续高度的功率谱图，“+”号为自动识别的谱峰位置，图中连续4个距离库高度上均存在能量较强的地物杂波干扰，幅度远大于真实气象回波谱，在谱峰识别方面增加了难度，图9b是经过自适应阈值方法处理后得到的效果图，地物杂波得到抑制，有助于自动识别算法识别真实的气象回波。

图9 2014年5月17日12：00成都信息工程大学风廓线雷达功率谱随高度分布(a)原始功率谱随高度分布，(b)自适应阈值方法处理后功率谱随高度分布Fig.9 Spectral distribution of WPR at CUIT with height at 1200 BT 17 May 2014(a)original spectral distribution with height，(b)spectral distribution with height after self-adapting threshold processing

4 小 结

本文提出一种基于小波变换的自适应阈值去除地物杂波方法，通过理论分析、模拟数据及实测数据等多方面的检验论证，得到以下结论：

1) 与传统阈值方法、零频中心点剔除法进行对比，自适应阈值方法去除地物杂波效果显著，并证明db20小波在去除地物杂波应用中具有可靠性与有效性。

2) 自适应阈值方法在处理近零频信号、远零频信号以及由频谱泄露导致的零频信号展宽等情况时，均能准确地从地物杂波中分离出气象回波，且不需要设定参数以及选择模式。同时，根据对实测数据处理前后进行比较，使用该方法后信杂比可提高30 dB左右。

本文讨论的是单一地物杂波情形，没有将间歇性杂波纳入考虑范围内，因此,面对含有多类型杂波干扰的情况，自适应阈值方法还有待验证与完善；同时，本文的小波尺度分解系数由经验值确定，如何进行最优值选取，也需要进一步分析与推导。

目前，风廓线雷达数据处理中的积累与谱变换多已集成在信号处理板卡上，通过硬件设备实现，确定小波族后的小波变换算法也能采用同样的方式实现。同时，本文方法中的阈值是由信号小波分解后的高频系数关系确定，原理简单，计算量小，实现过程方便，应用于实际系统及业务中，能够简化气象目标识别过程，有效提升风廓线雷达产品数据的准确率。

[1] 阮征,吴志根.风廓线仪探测降水云体结构方法的研究.应用气象学报,2002,13(3):330-338.

[2] 何平,朱小燕,阮征,等.风廓线雷达探测降水过程的初步研究.应用气象学报,2009,20(4):465-470.

[3] 郭虎,段丽,杨波,等.0679香山局地大暴雨的中小尺度天气分析.应用气象学报,2008,19(3):265-275.

[4] 王莎,阮征,葛润生.风廓线雷达探测大气返回信号谱的仿真模拟.应用气象学报,2012,23(1):20-29.

[5] 王欣,卞林根,彭浩,等.风廓线仪系统探测试验与应用.应用气象学报,2006,16(5):693-698.

[6] 何平.相控阵风廓线雷达.北京:气象出版社,2006:11-12.

[7] Benjamin S G,Schwartz B E,Koch S E,et al.The value of wind profiler data in US weather forecasting.BullAmerMeteorSoc,2004,85(12):1871-1886.

[8] Kitamura Y,Nakagawa K,Sekizawa S,et al.Vertical Profile of Raindrop Size Distribution by Using 400 MHz Wind Profiler in Stratiform Rainfall∥32nd Conference on Radar Meteorology.2005:P7R.

[9] Wang M Z,Wei W S,He Q,et al.Application of wind profiler data to rainfall analyses in Tazhong Oilfield region, Xinjiang, China.JournalofAridLand,2012,4(4):369-377.

[10] Russell C A,Jordan J R.Portable Clutter Fence for UHF Wind Profiling Radar∥7th Symposium on Meteorological Observations and Instrumentations.1991:152.

[11] Strauch R G,Merritt D A,Moran K P,et al.The Colorado wind-profiling network.JAtmosOceanTechnol,1984,1(1):37-49.

[12] 刘涛,曾祥利,曾军.实用小波分析入门.北京:国防工业出版社,2006:73-75.

[13] 何平,李柏,吴蕾,等.确定风廓线雷达功率谱噪声功率方法.应用气象学报,2013,24(3):297-303.

[14] Sifuzzaman M,Islam M R,Ali M Z.Application of wavelet transform and its advantages compared to Fourier transform.JournalofPhysicalSciences,2009,13:121-134.

[15] Ingrid Daubechies.小波十讲.李建平,杨万年,译.北京:国防工业出版社,2005:29-31.

[16] Jordan J R,Lataitis R J,Carter D A.Removing ground and intermittent clutter contamination from wind profiler signals using wavelet transforms.JAtmosOceanTechnol,1997,14(6):1280-1297.

[17] Ai W,Huang Y,Hu M,et al.Ground Clutter Removing for Wind Profiler Radar Signal Using Adaptive Wavelet Threshold∥2010 International Conference on Measuring Technology and Mechatronics Automation (ICMTMA),IEEE,2010:370-373.

[18] Justen L A,Teschke G,Lehmann V.Wavelet-based Methods for Clutter Removal from Radar Wind Profiler Data∥Photonics Technologies for Robotics,Automation,and Manufacturing.International Society for Optics and Photonics,2004:157-168.

[19] Lehmann V,Teschke G.Wavelet based methods for improved wind profiler signal processing.AnnaleGeophysicae,1999,19(8):825-836.

[20] 丁敏,卜祥元.提升小波去地杂波方法.现代电子技术,2006,29(11):32-34.

A Modified Method of Removing Ground Clutter from Wind Profiler Radar Based on Adaptive Wavelet Threshold

Zhang Fugui1)2)Fan Xiao1)He Jianxin1)2)

1)(ChengduUniversityofInformationTechnology,Chengdu610225)2)(KeyLaboratoryofAtmosphericSoundingofCMA,Chengdu610225)

Wind profiler radar (WPR) can be used to retrieve real-time atmospheric wind field data of high resolution. Backscattered echo caused by irregularities of atmospheric refractive index is received by radar antenna and wind velocities is calculated with Doppler frequency shifting speed formula. It is widely used in fields of very short-term weather forecasting, airport operations and public protection, air pollution monitoring, wind field analyses and forecasts of toxic plume trajectories resulting from chemical or nuclear incidents. As a result of being widely used in different situations, WPR is always sited near the city with a large population and complicated geographical environment. Transmission of electromagnetic wave during WPR detecting period is often interfered by various clutters that contaminate WPR data introduce bias in moments and wind velocity estimation. Of all clutters, ground clutter is the primary source because it happens more often than the others. Ground clutter is radar return from more or less stationary targets such as trees, buildings near the cited place. How to eliminate the influence of ground clutter is a most concerned aspect. Ground clutter mainly concentrates around the zero-frequency and it has high amplitude on the power spectrum. The most frequently used methods, such as traditional wavelet threshold processing and zero-frequency elimination of 3 points, both have the ability to separate the meteorology echo from the ground clutter when the turbulent peak is away from the zero-frequency and not covered with ground clutter peak. However, when the near zero-frequency echo is taken into consideration, both of the traditional methods meet their limitation. Based on the wavelet high frequency coefficients, a method of determining threshold adaptively is proposed and the validation of the method is done by using of simulated data and WPR measured data. The corresponding power spectrum before and after self-adapting wavelet threshold processing are compared. Results show that this method performs well even when the signal is close to the zero-frequency and covered completely. Meanwhile, the method has some important features, such as simple theory, small amount of calculation and easy to implement. Cases analysis shows that self-adapting threshold processing can increase the accuracy of peak identification, also provide approach and basis for improving the WPR products.

wind profiler radar; ground clutter; wavelets analysis; self-adapting threshold processing

10.11898/1001-7313.20150409

国家公益性行业(气象)科研专项(GYHY200906039)

张福贵,范潇,何建新. 自适应阈值方法去除风廓线雷达地物杂波. 应用气象学报，2015，26(4)：472-481.

2014-11-13收到， 2015-03-18收到再改稿。

* email: zfg@cuit.edu.cn