1. 南京信息工程大学气象灾害预报预警与评估协同创新中心，江苏 南京 210044;
2. 广东省气象台，广东 广州 510640
收稿日期：2017-05-22；修订日期：2018-03-26
基金项目：国家自然科学基金项目(41475034);国家重点基础研究发展计划(973计划)项目(2013CB430101)共同资助
THE ALGORITHM AND VERIFICATION OF GROUND CLUTTER IDENTIFICATION BASED ON FUZZY LOGIC
1. Collaborative Innovation Center on Forecast and Evaluation of Meteorological Disasters, NUIST, Nanjing 210044, China;
2. Guangdong Meteorological Observatory, Guangzhou 510640, China
1 引言
新一代天气雷达不仅可以探测气象目标，同时也可以探测到非气象目标[1]。非气象目标主要包括地物、海浪以及生物(鸟、昆虫等)[2]，由这些目标形成的回波不利于业务人员对天气形势的理解和判断，特别是地物回波会对雷达数据质量造成严重的影响。因此，对于雷达数据应用，尤其对定量降水估计(QPE)而言，剔除地物目标的影响必不可少。目前，对于地物回波的剔除可以从3个方面进行：(1)雷达安装(选址、硬件等)；(2)数据处理；(3)与其他数据源进行对比[3]。国内外科研人员在这方面做了很多研究，并提出多种方法。Steiner等[3]提出使用雷达反射率因子的空间变化以及回波顶高作为判据的决策树方法来识别非气象回波。Zhang等[2]使用参考平面法并结合雷达反射率因子纹理以及垂直梯度对雷达数据进行质量控制。Cho等[4]使用动态隶属度函数的模糊逻辑方法剔除地物回波以及超折射回波。刘黎平等[5]提出分步式的模糊逻辑方法对超折射地物回波进行识别。庄薇等[6]采用距离加权法对应用于灾害性天气短时临近预报系统(SWAN)和新一代天气雷达业务系统(ROSE)中的地物回波识别算法进行了改进。
本文采用的方法基于模糊逻辑，且类似于Cho等[4]提出的方法，通过对不同反射率因子范围的雷达回波进行统计[7-8]，客观地得到隶属度函数以及权重。然而，本文提出的方法针对于降水强度量级的回波，即反射率因子不小于15 dBz，而对于非降水强度回波不进行处理，从而在保留对短临预报具有指示作用且强度较弱的特征回波(如晴空湍流回波以及阵风锋回波)的前提下[9-10]，最大程度消除地物回波对QPE、数据同化等雷达数据定量化应用的影响[11]。
2 数据
本文采用的数据来自于广州新一代天气雷达CINRAD/SA，雷达主要性能参数如表 1所示[12]。输出数据包括空间分辨率为1 km的反射率因子(Z)、250 m的径向速度(V)和250 m的谱宽(W)。体扫模式采用VCP21，即6分钟完成9层仰角扫描，对于最低的2个仰角采用高PRF和低PRF各扫描一次[1]。因此Z与多普勒数据(V和W)不是一一对应，且在计算反射率因子垂直梯度时还要对比上下两层的反射率因子大小，为此，本文采用刘黎平等[5]提出的方法，将雷达数据进行径向处理，按照严格的1 °间隔顺次排列径向数据，以满足同一距离库3种数据的对应，以及不同仰角反射率因子的对应。雷达站位于广州市番禺区，海拔高度为179 m。
表 1
表 1 广州新一代天气雷达主要性能参数[12]
| 物理量 |
性能参数 |
| 工作频率/MHz |
2 885 |
| 发射机峰值功率/kW |
650~800 |
| 脉冲重复频率/Hz |
318~1 304 |
| 波束宽度/° |
0.99 |
| 天线直径/m |
8.54 |
| 天线增益/dB |
45 |
| 动态范围/dB |
93 |
| 极化方向 |
水平 |
|
表 1 广州新一代天气雷达主要性能参数[12]
|
3 方法
3.1 特征参量
本文使用3个特征参量作为模糊逻辑输入：反射率因子平均径向纹理(TDBZ)、反射率因子垂直梯度(VGZ)和径向速度绝对值(Vabs)：
|
(1) |
|
(2) |
|
(3) |
式中，TDBZ可表征Z在某仰角径向方向的平滑程度，而降水回波明显比地物回波平滑，所以TDBZ可用来对二者进行区分，较大的TDBZ更有可能为地物回波。TDBZ在本文中的形式不同于其原始形式[13]，本文通过用对TDBZ取平方根的方式来缩小TDBZ的取值范围，从而便于对隶属度函数和权重进行统计。式(1)中i、j、k分别表示距离、方位角和仰角的索引。Nbeams和Ngates分别为计算TDBZ时沿方位角和径向所使用的距离库的个数。N = Nbeams × Ngates，代表使用距离库的总数。VGZ可表征三维反射率因子场在垂直方向的连续程度，由于地物回波只有在低层仰角才能被雷达探测到，因此地物回波通常具有较大的VGZ。本文中使用的VGZ类似于Cho等[4]，但利用高度代替仰角作为分母，从而可以抑制远距离处由波束展宽造成的影响[2]。式(2)中Z和H分别表示目标距离库的反射率因子值及其所在高度，Zup和Hup分别表示目标距离库的上层仰角相同位置处的反射率因子值及其所在高度。Vabs可以用来判断目标是否固定，对于降水目标来说，径向速度通常为非零值，因此具有较大径向速度则更有可能为降水目标。然而，单独使用上述3个特征参量中的任何一个时，都不能对非气象回波进行准确识别。因此通过模糊逻辑方法将它们以合适的权重进行组合，可以弥补各参量单独识别时的不足。
3.2 算法流程
算法流程图如图 1所示。首先通过输入的体扫数据计算特征参量TDBZ、VGZ和Vabs；将计算得到的特征参量代入预先得到的隶属度函数中计算隶属度值(MF(Z))。每一个特征参量有3个适用于不同回波强度范围(15~25 dBz、25~35 dBz和大于35 dBz)的隶属度函数及权重。通过对MF(Z)和对应权重(W(Z))的乘积进行累加得到总隶属度值(MFtot(Z))：
|
(4) |
将MFtot(Z)与预设的阈值(Thresh)进行比较，若MFtot(Z)>Thresh，则认为是地物回波。当距离较远时，降水回波与地物回波的区分趋于模糊，因此对于大于75 km且被认为是地物回波的距离库需要进行延展测试，即与目标距离库最近的8个距离库中大于50%同样被识别为地物回波或者无回波时，目标距离库才认为是地物回波。
在以上步骤执行完毕后，算法已完成对地物回波的识别，将被识别为地物回波的距离库设置为无回波即可剔除地物回波。但是，当地物回波嵌入在降水回波之内时，仅仅将地物回波剔除并不能消除其对雷达数据质量造成的影响。简单地对地物回波进行剔除不仅会造成回波不连续，还会造成对降水量的低估[14]。因此，利用上层仰角反射率因子数据对剔除地物回波后形成的“空洞”进行填补，可进一步提高QPE等雷达数据定量应用的精度。今后可考虑使用垂直反射率廓线(VPR)算法进行订正[15]。
本文根据降水回波在垂直方向上的连续性，对因剔除地物回波而造成的“空洞”进行填补，从而进一步还原真实的降水回波。(1)对地物回波上层相同位置处的距离库进行检测，如果存在回波，则使用上层距离库对“空洞”进行填补。(2)考虑到回波水平连续性，本算法对填补后的反射率因子场进行平滑处理。平滑采取“9点平均”方式，即取目标距离库及其周围最近的8个距离库的平均值作为平滑结果并赋值给该距离库。平滑处理的范围为被识别为地物回波的距离库及其周围最近的8个距离库。经过平滑处理后的数据为本算法的输出数据。
4 隶属度函数及其权重
4.1 雷达回波特性统计
本文选取200次体扫的降水数据以及200次体扫的晴空数据用于统计由特征参量频率分布所表示的雷达回波特性，其中降水数据主要包括层状云降水、混合型降水、孤立单体和飑线过程。通过对多个晴空条件下体扫数据进行平均来获得地物回波分布图，并根据此分布图将雷达回波分为G和GF两部分，分别代表地物回波区和非地物回波区。地物回波特性由200次晴空体扫中的G区域数据获得，降水回波特性由200次降水体扫中的GF区域数据获得。图 2~4分别为3个特征参量(VGZ、TDBZ和Vabs)在3个反射率因子范围(15~25 dBz、25~35 dBz和大于35 dBz)的频率分布来代表降水回波。随着反射率因子增加，降水回波的VGZ分布除略微变窄之外并无明显变化，而地物回波的变化较明显，表现为分布逐渐向VGZ增加的方向移动，与降水回波分布的重叠面积逐渐减小。因此，VGZ对于地物回波的识别效果随着反射率因子增加而得到改善。与VGZ不同，地物回波的TDBZ分布不随反射率因子增加而明显改变，且降水回波的分布随反射率因子的增加而逐渐展宽，表明TDBZ对于地物回波的识别能力随反射率因子增大而减弱，与Zhang等[2]所得结论相同。降水回波的Vabs分布随反射率因子增加而均匀展开，而地物回波的Vabs分布则随反射率因子的增加逐渐变窄，表明反射率因子越大，Vabs对于地物的识别效果越好。
4.2 隶属度函数与权重
利用统计得到的降水和地物的回波特性可以构建识别地物的隶属度函数。对于某个给定特征参量(如VGZ)以及反射率因子范围，可以得到地物回波和降水回波的频率分布FG和FP，则地物回波的隶属度函数为：
|
(5) |
式中，Z代表图 2~4中不同的反射率因子范围，下同。经式(5)计算得到的隶属度函数如图 5所示，值得注意的是隶属度函数随反射率因子范围的变化而发生改变。例如，当VGZ=14 dB/km时，反射率因子间隔取15~25 dBz对应隶属度函数值接近于1，反射率因子间隔取25~35 dBz对应隶属度函数值大致为0.5，反射率因子间隔取大于35 dBz对应隶属度函数值为0。因此，针对不同的反射率因子间隔应采用相应的隶属度函数有助于提高识别效果。
由于在某一特定反射率因子间隔内不同特征参量的识别效果不同，因此，对于不同特征参量需要分配不同的权重。本文采取Cho等[4]的方法，通过计算地物和降水回波特性频率分布(图 2~4)中的重叠面积得到的相应权重，计算公式为：
|
(6) |
|
(7) |
|
(8) |
|
(9) |
式中，A代表对应反射率因子范围和特征参量的重叠面积，S为对应反射率因子范围各特征参量重叠面积倒数之和。
5 效果检验
本文采用两种方法对识别算法效果进行检验：(1)选取典型个例，将本文提出的算法处理结果与原始数据以及文献[5]所提出的算法(简写为R5)进行对比；(2)采用统计方式，对比算法处理前、后雷达回波强度的频率分布。
5.1 个例检验法
选取4个典型个例(图 6~9)，分别为层状云降水、混合型降水、孤立单体和飑线，其中图a、b分别代表原始数据和本文提出算法处理结果，图c、d分别为本算法未经填补之前的结果和R5的处理结果。如图 6a、6b、7a、7b所示，在大范围降水条件下，本算法可对地物回波进行准确识别，并通过填补剔除地物回波所产生的“空洞”使得降水回波连续。如图 8a、8b、9a、9b所示，在强对流天气条件下，本算法对强地物回波和强降水回波具有很好的辨识度，在保留降水回波的前提下，最大程度地对地物回波进行剔除。通过对比图 6~9的图c、d可以发现，在降水回波边缘处理以及孤立强对流单体识别等方面，本文提出的算法优于R5。主要原因在于本算法通过对Z进行分段处理，从而在降低弱回波处误判的同时提高在强回波处的识别度。另外，本文提出的算法采用统计得到的隶属度函数和权重，较R5更适合当地的特定气候环境。
5.2 统计检验法
本文采用广州新一代天气雷达2011年全年体扫数据进行统计分析。统计各距离库反射率因子取值范围在15~25 dBz、25~35 dBz和大于35 dBz的频率，计算公式为：
|
(10) |
式中，x、y代表目标距离库所在位置(第x根径向，第y个距离库)，N为目标距离库所在位置处反射率因子在Z范围内出现的次数，Tot为用于统计的总体扫数。将得到的频率采用dB形式表示：
|
(11) |
统计结果如图 10~12所示，其中图a和b分别代表原始数据和经过算法处理之后的数据，辐辏状区域为部分阻挡导致。经算法处理后，不同反射率因子范围均有显著改善，具体表现为频率偏高或偏低区域明显减少，以及地物区域频率取值与周围区域趋于近似。另外，随反射率因子增加算法对于数据的改善效果愈发明显。但随距离增加算法处理能力不断减弱，表明本文所提出的算法对于近距离强回波有更好的效果。
6 结论
本文采用模糊逻辑方法对地物回波进行识别，并根据反射率因子范围采用不同的隶属度函数与权重，且二者是通过典型个例的回波特性进行统计得到，从而保证了客观性且更加适应当地气候条件。该方法针对降水强度量级的回波，即反射率因子不小于15 dBz，对于非降水强度回波则不进行处理，从而保留对短临预报具有指示作用且强度较弱的特征回波，如晴空湍流回波以及阵风锋回波。当地物回波嵌入在降水回波之内时，简单地对地物回波进行剔除不仅会造成回波不连续，同时还会造成对降水量的低估。因此，利用上层仰角反射率因子数据对剔除地物回波后形成的“空洞”进行填补，可进一步提高QPE等雷达数据定量应用的精度。今后可考虑使用垂直反射率廓线(VPR)算法进行订正。
本文采用两种方法对识别算法效果进行检验。个例检验法通过将本文提出的算法处理结果与原始数据和R5的处理结果进行对比，表明算法对于层状云降水、混合型降水、孤立单体以及飑线条件下的地物均有很好的识别效果，在降水回波边缘处理以及孤立强对流单体识别等方面，本文提出的算法R5更有优势。统计检验法通过对广州雷达2011年数据进行统计并对比算法处理前、后雷达回波强度的频率分布对算法进行检验，得到经算法处理后地物回波区的频率分布与周围区域趋于近似，而且随反射率因子的增加效果愈明显，但算法识别效果随距离的增加而减弱，表明本文所提出的算法对于近距离强回波有更好的效果。
2018, Vol. 34
