长沙机场阵列天气雷达地物识别算法

魏万益 马舒庆 杨 玲* 甄小琼4) 吕寺炜

1)(成都信息工程大学电子工程学院， 成都 610225)2)(中国气象局大气探测重点开放实验室， 成都 610225)3)(中国气象局气象探测中心， 北京 100081)4)(中国科学院大气物理研究所， 北京 100029)5)(雷象科技(北京)有限公司， 北京 100089)

引 言

地物杂波识别和去除是天气雷达数据质量控制研究中的一项基本工作。地物会影响雷达回波，掩盖天气回波，甚至直接影响雷达产品的可信度。

为了消除地物杂波，国内天气雷达专家们提出了多种去除地物杂波的方法，这些方法可分为以下3类：①通过改变雷达的设计和布局位置减少地物杂波，②在信号处理端对雷达信号进行处理，③针对雷达基数据进行后期处理。Smith[1]和Mann等[2]通过选择雷达的安放位置和采用更短的波长减少杂波及其相对强度。第2类方法通常是对收集到的所有数据采用时域[2-4]或频域[5-7]滤波器去除速度接近0的回波数据。Li等[8]在频谱上提取特征参数作为朴素贝叶斯算法的输入，以判断是否存在地物杂波并加以去除。针对雷达基数据的处理方法主要是从基数据中提取二维或三维地物识别特征参数。Steiner等[9]使用回波强度三维结构作为决策树算法输入参数构建自动化杂波抑制算法，发现最有效的3个参数是雷达回波的垂直范围、回波强度的水平变化及回波强度的垂直梯度。Zhang等[10]使用类似的回波强度垂直梯度去除地物杂波，并在计算该参数时考虑仰角抬高对水平距离的影响。Lakshmanan等[11-12]采用神经网络方法，将回波强度的三维结构、速度和谱宽以及其他特征作为神经网络的输入数据去除地物杂波。此外，Kessinger等[13]、Ferrer等[14]、Cho等[15]、刘黎平等[16]和江源等[17]在雷达基数据中提取各种特征，使用模糊逻辑算法识别并去除地物杂波。李丰等[18-19]采用类似的模糊逻辑算法识别非降水回波，并对比分析S波段和C波段多普勒天气雷达各种特征的分布。

Ruiz等[20]使用相控阵天气雷达探测数据，提取回波强度与时间的相关性纹理特征等参数，并将其用于朴素贝叶斯算法，有效识别和去除地物杂波。但该算法运算效率较低，不能满足快速扫描相控阵天气雷达的产品生成需求。文中同时指出，与传统抛物线天线天气雷达相比，日本大阪大学的相控阵天气雷达具有较强的旁瓣和较宽的波束，会导致雷达数据中存在更多的地物杂波，对地物杂波的识别和去除算法要求更高。

为了进一步研究小尺度强对流天气系统，在国内外网络化和相控阵技术发展的背景下[21-23]，2015年中国气象局气象探测中心设计并联合有关厂家开始研制阵列天气雷达，它结合网络化雷达和分布式相控阵技术，具有高时空分辨率的特点。本文采用基于模糊逻辑的地物识别算法，参考日本大阪大学的相控阵天气雷达地物杂波去除算法，提出了一个新的参数——回波强度时间变化量(time variability of reflectivity factor,TVR)，并在不同天气情况下验证了地物杂波去除算法的有效性。

1 阵列天气雷达简介

阵列天气雷达是一种分布式相控阵天气雷达，由控制处理中心和至少3个相控阵接收发射子阵(简称收发子阵)构成，通过增加收发子阵数量扩大探测区域[24]。每3个相邻的收发子阵为一组进行协同扫描，保证3个相邻收发子阵在空间同一点的数据时差小于2 s。

收发子阵是阵列天气雷达的前端，由于缺少独立的控制和数据处理部分，不是完整的雷达，故被称为收发子阵而非相控阵雷达。收发子阵由相控阵天线阵列、收发模块阵列、信号处理阵列、方位旋转伺服、同步和通信模块等部分组成。水平方向采用机械扫描方式覆盖0°～360°，垂直方向发射4个宽度为22.5°的波束覆盖0°～90°，接收到回波信号后，通过数字波束得到16个平均波束宽度为1.6°的波束。2018年3月在长沙机场布设的阵列天气雷达，天线旋转速度为30 °·s-1，完成1个体扫仅需要12 s，约为现在业务中使用的多普勒天气雷达常用体扫模式1个体扫时间的1/30。每个收发子阵的最大探测距离为20.28 km，径向分辨率为30 m。

2 基于模糊逻辑的地物识别算法

2.1 数据预处理

由于雷达基数据中含有部分孤立的点状或线状回波(图1a)，本文采用Zhang等[10]提出的孤立点滤波算法对雷达基数据进行预处理。具体算法如下：以空间中点X为中心位置，在方位、径向和俯仰上5×5×1个距离库中有效回波点的百分比为PX，如果PX小于预设阈值(默认为75%)，则认为该点为孤立点并将其剔除。图1为孤立点滤波前后的回波对比。由图1b可知，通过孤立点滤波可以消除大部分孤立回波点。

图1 孤立点滤波前(a)、滤波后(b)回波强度Fig.1 Intensity echo of isolated point before filtering(a) and after filtering(b)

2.2 识别算法介绍

阵列天气雷达的时间分辨率空间分辨率均高于常规多普勒天气雷达。本文在已有地物特征参数基础上加入表征回波强度随时间变化的参数,回波强度时间变化量(TVR)，作为模糊逻辑算法的输入参数来识别地物杂波[25-27]。其他特征参数包括回波强度纹理(TDBZ)、回波强度垂直梯度(GDBZ)、径向库间回波强度变化程度(SPIN)、径向速度平均值(MDVE)和速度谱宽平均值(MDSW)。这些参数的定义见式(1)～(7)[16]：

(1)

(2)

(3)

(4)

(5)

(6)

(7)

其中，NA，NR表示在方位和径向定义的计算范围；i，j为范围内的索引；Z，V和SW分别为回波强度、径向速度和速度谱宽；Zlow和Zup为对应的本层和上层PPI的回波强度；Zthresh为径向上相邻距离库强度变化的阈值；ZN和ZL分别表示当前时刻和上一时刻的强度。在实际特征值计算过程中，设定NA=5，NR=5，Zthresh=2 dBZ，ZN和ZL的时间差为1 min。

除少数特殊地物外，大部分地物位置固定不动且形状大小不随时间改变，所以地物杂波的径向速度主要分布在0附近，谱宽相对较小，回波强度时间变化量(TVR)接近0。相比降水回波，地物杂波大部分出现在低仰角数据中，强度在空间上变化大，且随着仰角的抬高，强度迅速减小。因此地物杂波的强度纹理(TDBZ)和径向库间强度变化程度(SPIN)较大，强度的垂直梯度(GDBZ)为负值且绝对值比较大。

本文收集了2019年3—7月不同天气情况下长沙机场阵列天气雷达基数据，根据地物杂波和天气回波的特征预先确定回波类型。判断依据：地物主要出现在较低的仰角，且随着仰角的抬高强度迅速减小，地物的径向速度接近零，位置不随时间改变；对于降水回波，可根据回波形状、垂直结构、演变和移动等方面确定[16]。使用部分回波数据，统计得到地物和降水各个参数的概率分布情况(图2)，其余回波数据验证，检验算法的识别效果。

由图2可知，地物的径向速度平均值(MDVE)集中分布在0附近，与降水的区别较大(图2a)。地物和降水回波的速度谱宽平均值(MDSW)重合面积较大，区别较小，未达到预期效果，因此，在识别算法中剔除该参数(图2b)。与降水回波相比，地物的回波强度纹理(TDBZ)分布更加离散，符合地物杂波空间分布差异较大的特征(图2c)。地物杂波的强度垂直梯度(GDBZ)存在少数大于0的部分，与预期的效果存在差异，但与降水回波明显不同，不影响作为识别地物的输入参数(图2d)。地物和降水回波的强度径向库间变化程度(SPIN)分布的波峰位置区别明显，地物杂波的波峰位置更靠右，即在计算区域内，地物杂波强度波动超过2 dB的数据点更多(图2e)。回波强度时间变化量(TVR)在地物和降水回波中有明显差异,地物杂波主要分布在0～2 dB范围内，而降水回波的分布则相对离散(图2f)。综上所述，地物和降水回波的大部分参数存在较明显差异，但也有部分重合区域，如在回波强度纹理(TDBZ)靠近0的区域也有部分地物杂波，径向速度平均值(MDVE)较小的区域也有降水回波。因此，使用单一特征参数难以识别地物杂波，需综合使用多种特征参数才能较好地区分地物和降水回波。

为了简化计算，本文根据地物和降水特征分布的范围、交点和波峰位置确定采用梯形隶属函数(图3)。根据隶属函数，计算每个点各参数的判据M。然后，对每个点的所有判据(N个)进行累加平均(式(8))，可得到每个点地物的判断值EM(0～1)。该值越大，代表该点是地物的可能性越大，反之，代表降水回波的可能性越大。最后将得到的判断值EM与预先设定的阈值进行比较，大于阈值判断为地物，反之判断为降水。

(8)

图2 地物和降水回波各特征参数的概率分布(a)MDVE，(b)MDSW，(c)TDBZ，(d)GDBZ，(e)SPIN，(f)TVRFig.2 Probability distributions of ground clutter and precipitation echo(a)MDVE，(b)MDSW，(c)TDBZ，(d)GDBZ，(e)SPIN，(f)TVR

图3 地物识别的隶属函数(a)MDVE，(b)TDBZ，(c)GDBZ，(d)SPIN，(e)TVRFig.3 Membership functions of ground clutter detection(a)MDVE，(b)TDBZ，(c)GDBZ，(d)SPIN，(e)TVR

3 地物识别效果

不同阈值的回波强度时间变化量(TVR)参数对识别效果影响见表1。随着阈值的提高，地物识别准确率和降水误判率降低，在阈值不大于0.55时，TVR对改善地物识别准确率和降水识别误判率有贡献。当阈值大于0.55后，TVR不再对识别算法产生影响。

虽然阈值为0.4时地物识别的准确率最高，但降水识别的误判率也最大。综合考虑地物识别的准确率和降水识别误判率，最终选取阈值为0.45，并用于分析3个子阵在无降水、混合性降水和对流性降水天气情况的地物识别算法效果。

表1 地物识别准确率和降水识别误判率Table 1 Accuracy of ground clutter detection and error rate of precipitation detection

3.1 无降水情况地物识别效果

为了检验算法在无降水情况下的地物识别效果，采用2019年3—7月无降水的雷达数据，分析了3个子阵在无降水情况下算法识别效果(表2)。结果表明：地物识别算法在无降水情况下表现良好，地物识别的准确率高达96%。图4为子阵1在2019年7月31日10：21地物识别前后的对比，经过地物识别算法，大部分地物杂波可被有效识别并去除。

表2 无降水情况下地物识别算法效果Table 2 Accuracy of ground clutter detection algorithm under no precipitation condition

图4 2019年7月31日10：21子阵1在1.4°仰角的回波强度(相邻距离圈间距为5 km)(a)地物识别前，(b)地物识别后Fig.4 The echo intensity of subarray 1 at 1.4° elevation angle at 1027 BT 31 Jul 2019(the distance between adjacent rang rings is 5 km)(a)before ground clutter detection,(b)after ground clutter detection

3.2 混合性降水情况地物识别效果

分析混合性降水情况下地物识别算法效果(表3)。在此类降水情况下，各个子阵地物杂波识别的准确率均达到91%以上，而降水回波识别的误判率最高仅10%。图5为子阵1在2019年6月21日15：17地物识别前后的对比。由图5可知，降水识别误判率较高主要由降水回波边缘部分的误判引起。相对于无降水情况，地物识别准确率略有下降，从速度对比可以看到，与降水回波重合的地物杂波部分未被准确识别。

表3 混合性降水算法识别效果Table 3 Accuracy of ground clutter detection and error rate of precipitation detection under mixed prcipitation condition

图5 2019年6月21日15：17子阵1在1.4°仰角的回波强度和径向速度(相邻距离圈间距为5 km)(a)识别前的回波强度，(b)识别后的回波强度，(c)识别前的径向速度，(d)识别后的径向速度Fig.5 The echo intensity and radial velocity of subarray 1 at 1.4°elevation angle at 1517 BT 21 Jun 2019(the distance between adjacent range rings is 5 km)(a)echo intensity before ground clutter detection,(b)echo intensity after ground clutter detection,(c)radial velocity before ground clutter detection,(d)radial velocity after ground clutter detection

3.3 对流性降水情况地物识别效果

分析对流性降水情况下地物识别算法效果见表4。由于子阵1和子阵3的个例中地物杂波和降水回波无重合部分，地物识别的准确率与无降水情况下差别不大，分别达到92%和94%。子阵2的大部分个例中降水和地物杂波重合区域中存在不能被准确识别部分，准确率仅91%。图6和图7分别是2019年7月21日13：57子阵2的1.4°和2.8°仰角地物识别前后的对比，2.8°仰角地物明显减少，降水区域被更完整保留。因此，该类型降水情况下，算法识别处理后，大部分降水回波被较好保留，仅降水回波边缘部分和西南方向一小块孤立的降水回波被错误识别。

表4 对流性降水算法识别效果Table 4 Accuracy of ground clutter detection and error rate of precipitation detection under convective precipitation condition

图6 2019年7月21日13：57子阵2在1.4°仰角的回波强度和径向速度(相邻距离圈间距为5 km)(a)识别前的回波强度，(b)识别后的回波强度，(c)识别前的径向速度，(d)识别后的径向速度Fig.6 The echo intensity and radial velocity of subarray 2 at 1.4°elevation angle at 1357 BT 21 Jul 2019(the distance between adjacent range rings is 5 km)(a)echo intensity before ground clutter detection,(b)echo intensity after ground clutter detection,(c)radial velocity before ground clutter detection,(d)radial velocity after ground clutter detection

图7 2019年7月21日13：57子阵2在2.8°仰角的回波强度和径向速度(相邻距离圈间距为5 km)(a)识别前的回波强度，(b)识别后的回波强度，(c)识别前的径向速度，(d)识别后的径向速度Fig.7 The echo intensity and radial velocity of subarray 2 at 2.8° elevation angle at 1357 BT 21 Jul 2019(the distance between adjacent range rings is 5 km)(a)echo intensity before ground clutter detection,(b)echo intensity after ground clutter detection,(c)radial velocity before ground clutter detection,(d)radial velocity after ground clutter detection

续图7

4 小 结

本文利用阵列天气雷达高时空分辨率的特点，提出回波强度时间变化量(TVR)这一新参数。采用2019年3—7月长沙黄花机场X波段阵列天气雷达数据，分析地物识别算法中6个特征参数的统计特征，改进基于模糊逻辑的地物识别算法，并利用雷达基数据分析TVR对算法的贡献和不同天气情况下识别算法的效果，得到以下结论：

1) 降水回波和地物杂波回波强度纹理(TDBZ)、回波强度垂直梯度(GDBZ)、径向库间回波强度变化程度(SPIN)、回波强度时间变化量(TVR)和径向速度平均值(MDVE)5个参数特征分布差异较大，但速度谱宽平均值(MDSW)差异较小，因此，只采用前5个参数作为模糊逻辑算法的输入特征。

2) 对回波强度时间变化量(TVR)的贡献程度进行定量评价，通过选择最优的阈值，采用该参数可将地物识别准确率提高4%，降水识别误判率降低2%。

3) 在不同天气情况下分别对地物识别算法的表现进行验证。无降水时，地物识别的准确率为94%；有降水时，地物识别准确率为92%，降水识别误判率稳定在10%左右。

本算法在阵列天气雷达地物识别中取得了较好的效果，为生成更加准确的气象产品提供了保障。但还存在一些问题：降水回波的边缘部分可能被误判，且当降水回波与地物回波混合在一起时，去除地物回波也损失了降水回波信息。在后续的研究中，将考虑采用多子阵联合的地物去除方法，进一步提高识别准确率，以及补偿与地物一起去除的降水回波。