基于二维统计分析的降水粒子分类方法

降水粒子分类研究对于灾害天气的识别、预警、数值预报以及人工影响天气的作业指挥和效果评估[1]都有十分重要的意义。与传统多普勒气象雷达相比，双偏振气象雷达可以交替发射水平和垂直偏振波，并接收2个偏振方向的回波信号，除获取反射率强度、径向速度、速度谱宽外，还可探测差分反射率因子、差分传播相移、相关系数等参量，具有识别粒子相态功能，是近年降水粒子分类研究领域的研究热点[2]。

模糊逻辑方法是最经典的降水粒子分类方法，是一种解决不精确和不完全信息问题的有效方法, 最大的特点是可以利用专家经验进行分类，能够比较自然地处理模糊的概念, 其主要的优点是可以处理不确定的信息。20世纪90年代初，Straka等[3]将模糊逻辑应用于双偏振雷达的降水粒子相态识别研究中，与经典布尔逻辑方法相比，模糊逻辑方法以隶属函数代替经典布尔逻辑的固定阈值，同时减少了判断步骤。Zrnic等[4]将隶属函数由二维矩阵变成了5个极化参量的矩阵，同时可以识别9种降水粒子。曹俊武等[5]将模糊逻辑方法应用于国内的双偏振雷达。Park等[6]改进了模糊逻辑相态识别算法(Hydrometeor Classification Algorithm, HCA)，将隶属函数修改为非对称的梯形，并引入了置信向量来平衡数据误差。

近几年，基于模糊逻辑的双偏振气象雷达方法在国内应用开始兴起。文献[5]中利用模糊逻辑算法对上海的S波段双极化气象雷达实测数据进行降水分类的业务分析。文献[7]中利用模糊逻辑算法对X波段双极化气象雷达探测降水数据进行分类，分析研究了各降水粒子在雷暴演变过程中的相态特性。汪舵等[8]建立了一种基于相态识别的S波段双线偏振雷达最优化定量降水估测算法。文献[9]以美国HCA方法为基础，建立了适合广东S波段双偏振气象雷达的分类方法。武崇等[10]对一次超级单体进行讨论，用HCA方法识别出雨及冰雹混合物区域。同年，他通过对飑线过程的分析[9]，对HCA算法的参数和阈值进行优化，得到适用于华南珠海的偏振算法，该研究为发展适合中国降水特征的 HCA 算法奠定基础。王洪等[11]将Park[6]的算法简化并应用于珠澳偏振雷达中，通过S波段雷达对一次超级单体进行分析，有效识别出冰雹区域，由于仅对一次观测结果得出结论，因此对各偏振量之间的定量关系也需要进一步研究。杨章等[12]基于HCA-LIQ算法对广州9次降雨评估，由于该算法低估了中小颗粒的沉淀，过高估计了中等和大颗粒的沉淀，导致对大雨的分类不准确，需要研究极大雨来进一步验证并改进算法。冯亮等[13]利用双偏振气象雷达数据建立了基于偏振参量、温度、纹理参数的模糊逻辑降水粒子识别算法。文献[14]基于深度学习和模糊逻辑算法对降水粒子进行分类，该方法提高了水凝物识别的准确性。王金虎等[15]利用模糊逻辑算法和神经网络算法对云雷达联合微波辐射计探测数据进行降水粒子分类研究，结果表明模糊逻辑算法对退极化比极为敏感。发展至今，模糊逻辑方法仍存在一定的局限性，例如怎样确定隶属函数参数，以及如何提高降水粒子分类的准确率。

本文使用二维统计分析的方法确定钟型隶属函数的参数。首先将数据进行二维统计，分析得到钟型隶属函数的参数中心以及宽度的取值，然后取隶属函数图像的临界值作为钟型隶属函数斜率。再根据不同极化参数对降水粒子的贡献程度以及分类准确率调节各极化参数的权重系数。最后使用调节好的参数以及权重系数对一次雷达数据进行降水粒子分类。该方法能够较为准确的调节隶属函数参数，实现降水粒子分类。

1 基于二维统计分析的模糊逻辑降水粒子分类方法

模糊逻辑方法是最经典的分类方法，主要分为3个部分：模糊化、规则计算和退模糊。模糊化主要是对隶属函数参数的确定，而隶属函数参数和权重系数直接影响规则计算的准确率，从而影响退模糊后的降水粒子分类准确率，退模糊中引入融化层信息调节降水粒子分类结果。隶属函数参数一般是使用专家经验值进行确定[16]，不同隶属函数的参数也不同，本文使用二维统计分析方法对钟型隶属函数参数进行确定，使用模糊逻辑方法实现降水粒子分类，流程如图1所示。在模糊化过程中，首先对带标签的不同降水粒子各极化参数的数据进行处理，然后对处理后的数据进行二维统计分析，对钟型隶属函数的中心、宽度进行确定，再取临界值作为斜率的取值；在规则计算过程中，调节不同降水粒子各极化参量的权重系数，得到参数调整的结果；最后将调节好的参数使用模糊逻辑方法对一次雷达数据进行降水粒子分类。

1.1 数据处理

本文使用4个极化参数进行识别、分类：水平反射率因子ZH，差分反射率因子ZDR，互相关系数ρhv和差分相移率KDP，其中差分相移率KDP取对数形式LKDP。

1) 水平/垂直反射率因子ZH/ZV(单位：dBz)

双偏振气象雷达发射水平和垂直偏振方向电磁波，通过雷达的回波数据可以计算得出水平反射率因子ZH和垂直反射率因子ZV。

2) 差分反射率因子ZDR(单位：dBz)

差分反射率ZDR表示为ZH和ZV的差，计算方式为

ZDR反映了降水粒子下落过程中的形变程度，形变程度与降水粒子的空间取向和轴比(粒子的长轴与短轴的比值)有关，粒子形状越接近圆球形时ZDR的值越接近于0。

3) 互相关系数ρhv

ρhv是水平和垂直偏振波的后向散射特征的相关性。定义为水平与垂直偏振回波信号的互相关系数。

4) 差分相移率KDP(单位:(°)/km)

KDP不受雷达定标、衰减和波束遮挡等影响[17]，在液态的降水粒子中一般为正值，在垂直取向的柱状冰晶处会出现负值。由于KDP参量较大的变化范围不利于降水粒子的相态识别，因此在KDP的基础上引入了LKDP，即使用对数坐标替代线性坐标，具体如下：

1.2 二维统计分析方法确定钟型隶属函数参数

1) 钟型隶属函数

不同降水粒子类型各极化参量的范围通常以适当的隶属函数为特征，本文使用钟型隶属函数，一方面，钟形隶属度函数具有扁平宽阔的区域，它对应于降水粒子的参数取值范围，而不是单一的数值；另一方面，钟形隶属度函数有一个很宽的过渡区，可以改善模糊逻辑识别算法在过渡区的可靠性。钟型隶属函数表达式如式(3)所示：

其中，i=1,2,…,6，代表6种分类结果，本文将双偏振雷达回波识别为6种降水粒子分类：1)冰晶(IC);2)干雪(DS)；3)湿雪(WS)；4)雨(RA)；5)霰(GR)；6)冰雹(HA)。j=1,2,3,4，代表有4种用于相态识别的雷达极化参量，xj代表4种极化参量，分别是：1)ZH；2)ZDR；3)ρhv；4)LKDP。mi_j是隶属函数的中心，ai_j是隶属函数宽度的一半，bi_j是隶属函数的斜率，每个参数的意义如图2所示，横坐标代表极化参数，纵坐标代表隶属度，其中，不同的分类结果的各极化参数有不同的隶属函数范围。

2) 中心值mi_j及宽度ai_j的确定方法

根据带标签的数据，对不同降水粒子的各极化参数数据进行统计，分析确定每一类隶属函数的参数取值。即需要分析4个极化参数：水平反射率因子ZH，差分反射率因子ZDR，互相关系数ρhv和差分相移率LKDP。本文将ZH与其他极化参数进行二维统计分析，得到ZH_ZDR、ZH_ρhv以及ZH_LKDP的频次图，然后基于点与钟型隶属函数区域之间的相对位置，使数据大部分集中在钟型隶属函数包含的区域内，以封闭区域(如图3中黑色实线所包含的区域)代表钟型隶属函数，封闭区域内降水粒子的隶属度都大于等于0.5，封闭区域外，降水粒子的隶属度逐渐趋于0。根据封闭区域的边界值

通过公式(4)、(5)得到钟型隶属函数的中心mi_j和宽度ai_j。

图3是IC的ZH_ZDR二维数据统计(频次图)结果，以图3横坐标为例，封闭区域有两个边界值

从而得到IC水平反射率因子ZH的宽度a1_1=5 dBz和中心值m1_1=14 dBz。同理可得不同降水粒子各极化参数的中心和宽度。

3) 斜率bi_j的确定方法

斜率的取值一般是0～20范围内的数字，斜率太小，不同降水粒子的隶属函数重叠区域比较大；斜率太大，钟型隶属函数所包含的区域几乎不会发生变化，因此，本文分别将HA和RA、DS和IC、WS和GR的隶属函数对比，先画出斜率值为0.5时的钟型隶属函数图像，然后不断以0.5～1的间隔增大斜率的取值，得到多个隶属函数图像，经过对斜率取值的调整使钟型隶属函数所包含区域几乎不会发生变化，确定临界值，取临界值作为斜率bi_j。

图4代表HA和RA的ZH的钟型隶属函数，实线代表HA的隶属函数，虚线代表RA的隶属函数，不同颜色代表不同斜率的取值。由图可得，当HA斜率为2，RA斜率为4时，钟型隶属函数之间重叠区域较大，当HA的斜率大于10，RA的斜率大于12时，钟型隶属函数区域不会发生变化，所以选取HA的ZH钟型隶属函数斜率为b6_1=10，RA的ZH钟型隶属函数斜率为b4_1=12。

1.3 权重系数Wij确定方法

权重系数一般由综合理论值与专家的经验确定[18]，本文以文献[6]给出的权重系数为参考，如表1所示，权重系数的取值一般是0,0.2,0.4,0.6，0.8，1.0。根据极化参数对降水粒子识别的贡献程度高低以及分类结果准确率调节权重系数大小，提高识别能力更强的极化参数的权重能有效提高准确率。本文根据极化参数的贡献程度来调节权重系数大小，并使用混淆矩阵来判别不同分类结果的准确率，取准确率最高的权重系数作为最终的调整结果。

混淆矩阵如图5所示，横坐标代表修改参数后的分类结果，纵坐标代表数据的标签，对角线代表不同类别的准确率。图5(a)是未调节权重系数的分类结果准确率，识别效果不好的IC与DS主要在ZH隶属函数区域内更容易区分，在ZDR、ρhv、LKDP隶属函数区域重叠，因此，试验将IC的ρhv权重系数改为0.8，0.6，0.2，将DS的ρhv改为0.8，1.0，0.4，统计混淆矩阵的准确率，将准确率最高的结果对应的权重系数作为权重系数调节后的结果，确定将IC的ρhv权重改为0.6，DS的ρhv改为1.0。同理，HA和RA主要依靠LKDP隶属函数区域进行识别，考虑到HA的LKDP的权重已最高，仅修改RA的权重，试验将RA的ZDR改为1.0，0.6，0.4，0.2，统计混淆矩阵的准确率，确定将RA的ZDR权重系数改为0.6。WS和GR在ZDR隶属函数区域内更容易区分，考虑到权重系数已经最大，试验将GR的ρhv权重改为0.8，0.6，0.2，统计混淆矩阵准确率，确定将GR的ρhv权重改为0.6。由图5(b)可知，修改权重系数后，除雨外的其他降水粒子准确率均有所提高，说明调整权重系数后，结果更加准确。

1.4 双偏振气象雷达模糊逻辑降水粒子分类方法

模糊逻辑方法本质是简单的规则运算，其主要是将各极化参数对降水粒子的判别信息进行规则运算，得到分类结果，主要分为3个步骤进行分类：模糊化、规则计算和退模糊。

模糊化是将极化参数对应于不同降水粒子的特定范围，本文使用二维统计分析方法确定钟型隶属函数参数，并使用钟型隶属函数表示不同降水粒子各极化参数的范围。

规则计算是模糊逻辑最重要的一步，使用隶属函数的加权求和实现分类，在此过程中，主要确定权重系数的取值，本文通过混淆矩阵代表准确率来调节权重系数，使规则计算的结果更加准确，提高降水粒子分类准确率。函数Ai表示来自所有分类的总贡献，公式如式(6)：

其中P(i)(xj)是钟型隶属函数，表示第j类极化参数xj的范围; n=4，代表4个极化参数；Wij表示某一类的权重系数。

最后是退模糊，比较不同分类的Ai并找到具有最大值的类，作为最终输出结果C，其中C=i，代表6种降水粒子。在实际中，此过程应考虑一些物理约束。对于降水粒子分类，最常见的约束来自融化层。根据融化层顶部高度Ht和底部高度Hb，计算波束中心和-3 dB波束宽度与融化层的相对位置，得到雷达波束在融化层上下界Hb、Ht的投影距离Rtt、Rt、Rb、Rbb。

根据融化层的边界以及融化层的温度信息，对降水粒子分类结果做一定的限制[18]：雪只能出现在融化层底部及以上的高度；纯雨仅能出现在融化层顶部及以下的高度。根据上述原理，在不同距离范围内可以存在的降水粒子如表2所示。

本文通过二维统计分析方法对模糊逻辑钟型隶属函数参数修改，调节权重系数，并引入融化层信息对一次雷达数据进行分类，得到降水粒子分类结果。

2 算法流程与步骤

基于二维统计分析的模糊逻辑双偏振气象雷达降水粒子分类方法具体实现步骤如下：

步骤1 对数据进行处理，得到本文所需的4种极化参量；

步骤2 根据数据进行二维统计分析，确定钟型隶属函数参数；

步骤3 根据极化参数对降水粒子识别的贡献程度高低以及分类结果准确率调节权重系数大小；

步骤4 利用调节后的参数，并引入附加信息使用模糊逻辑方法对一次雷达数据进行降水粒子分类。

3 实验结果及分析

3.1 参数调整结果

1)中心mi_j以及宽度ai_j的取值

对不同分类结果数据进行二维统计分析，得到不同极化参数的频次图，再根据图6～图11 的边界值计算得到不同降水粒子各极化参数的隶属函数参数：中心mi_j和宽度ai_j。

2) 斜率bi_j的确定

不同降水粒子各极化参数的钟型隶属函数图如图12～图14所示，斜率太小，不同降水粒子的隶属函数重叠区域比较大；斜率越大，钟型隶属函数所包含的区域几乎不会变化，所以取临界值作为最终的斜率bi_j。

3) 钟型隶属函数参数取值

通过数据分析确定参数，得到模糊逻辑钟型隶属函数参数取值如表3所示。

4) 权重系数取值

权重系数是为了显示不同极化参数对不同降水粒子的贡献程度，通过混淆矩阵表示准确率来调整权重系数，调整后的权重系数如表4所示。

3.2 降水粒子分类结果

选取KTLX雷达2019年6月9日当地时间14：55数据，利用调节参数后的模糊逻辑方法进行降水粒子分类。图15(b)是调节参数后模糊逻辑的分类结果，同时与NOAA提供的降水粒子分类结果(图15(a))作为对比和参照。从图15(a)和图15(b)的对比中可以看出：调节钟型隶属函数参数后，模糊逻辑方法与NOAA提供的分类结果相似，仅存在一定程度的误识别，对于冰雹来说，识别区域要比NOAA提供的区域大，把少量的雨分为了冰雹。总体来说，与NOAA提供的结果相似度比较高。

4 结束语

针对模糊逻辑钟型隶属函数参数确定问题，本文使用二维统计分析的方法来确定钟型隶属函数的参数。首先将双偏振气象雷达的极化参数进行二维统计分析，得到钟型隶属函数参数：中心以及宽度，然后将临界值作为钟型隶属函数的斜率。再根据不同极化参数对降水粒子的贡献程度以及准确率调节降水粒子每一类极化参数的权重系数。使用调节好的参数以及权重系数对一次雷达数据进行降水粒子分类，实验结果表明，该方法能够较为有效地实现降水粒子分类，具有很好的研究潜质和应用价值。

[1] 楼小凤,师宇,李集明.云降水和人工影响天气催化数值模式的发展及应用[J].气象科技进展,2016,6(3):75-82.

[2] 刘黎平, 胡志群, 吴翀. 双线偏振雷达和相控阵天气雷达技术的发展和应用[J]. 气象科技进展, 2016, 6(3):28-33.

[3] STRAKA J M, ZRNIC D S. Algorithm to Deduce Hydrometeor Types and Contents from Multiparameter Radar Data[C]∥ 26th International Conference on Radar Meterology, Norman, Oklahoma:Amer Meteror Soc,1993:513-516.

[4] ZRNIC D S, RYZHKOV A V. Polarimetry for Weather Surveillance Radars[J]. Bulletin of the American Meteorological Society, 1999,80(3):389-406.

[5] 曹俊武，刘黎平，葛润生．模糊逻辑法在双线偏振雷达识别降水粒子相态中的研究[J].大气科学，2005,29(5):827-836．

[6] PARK H S, RYZHKOV A V, ZRNIC D S, et al. The Hydrometeor Classification Algorithm for the Polarimetric WSR-88D: Description an Application to an MCS[J]. Wea Forecasting, 2009,24(3):730-748.

[7] 谷娟，周筠，王君，等. 一次雷暴过程中水凝物相态演变特征[J].成都信息工程学院学报，2015,30(5):505-511.

[8] 汪舵，刘黎平，吴翀. 基于相态识别的S波段双线偏振雷达最优化定量降水估测方法研究[J].气象，2017,43(9):1041-1051.

[9] WU C, LIU L P, WEI M, et al. Statistics-Based Optimization of the Polarimetric Radar Hydrometeor Classification Algorithm and Its Application for a Squall Line in South China [J]. Advances in Atmospheric Sciences, 2018,35(3):296-316.

[10] WU C，LIU L P，LIU X，et al. Advances in Chinese Dual-Polarization and Phased-Array Weather Radars： Observational Analysis of a Supercell in Southern China[J].Journal of Atmospheric and Oceanic Technology，2018,35(9):1785-1806.

[11] 王洪，吴乃庚，万齐林，等.一次华南超级单体风暴的 S 波段偏振雷达观测分析[J].气象学报，2018,76(1):92-103.

[12] ZHANG Y，LIU L P，WEN H，et al. Evaluation of the Polarimetric-Radar Quantitative Precipitation Estimates of an Extremely Heavy Rainfall Event and Nine Common Rainfall Events in Guangzhou[J].Atmosphere，2018,9(9):330.

[13] 冯亮,肖辉,孙跃. X波段双偏振雷达水凝物粒子相态识别研究[J]. 气候与环境研究，2018,23(3):366-386.

[14] 冉元波. 基于模糊逻辑与深度学习的水凝物及云识别研究[D].成都：成都信息工程大学,2019.

[15] 王金虎,谢槟泽,刘光普,等.云雷达联合微波辐射计的水凝物识别研究[J].电波科学学报，2020，35(6):924-934.

[16] LIM S, CHANDRASEKAR V, BRINGI V N. Hydrometeor Classification System Using Duai-Polarization Radar Measurements: Model Improvements and in Situ Verification[J]. IEEE Trans on Geoscience and Remote Sensing,2005,43(4);792-801.

[17] 孙丝雨，沈永海，霍苗，等.双线偏振雷达在一次强降雹过程中的初步应用[J]. 暴雨灾害，2013,32(3):249-255.

[18] 吴翀. 双偏振雷达的资料质量分析，相态识別及组网应用[D]. 南京：南京信息工程大学，2018.