1 引言

冰雹灾害是由强对流天气系统引起的一种剧烈的气象灾害, 它突发性强, 常伴随着狂风、强降水等灾害性天气过程, 易给农业、建筑、通讯、电力、交通以及人民的生命财产安全带来较大的损失。由于冰雹过程具有突发性、局地性强以及生命周期短的特点, 其预报效果长期以来都不理想, 哪怕只是提前半小时的临近预报, 其进展也十分缓慢, 其中最大的困难就在冰雹回波的识别上^[1]。因此, 国内外学者对冰雹的识别做了大量的研究, 在粒子的微观形态方面, Hyang等^[2]对模糊逻辑分类算法进行了一些修改和改进, 发现了液态水凝物和固态水凝物之间的区别, 提高了分类的准确率。Ryzhkvo等^[3]研究了融化冰雹的衰减, 并提出了一种在冰雹存在下进行极化衰减订正的新技术。Zhang等^[4]对朴素贝叶斯的出色分类性能进出了新颖的解释。Manzato等^[5]研究了1992年4—9月的意大利北部Friuli Venezia Giulia地区平原上收集的冰雹数据, 使用线性和非线性方法进行多元预测分析, 并通过神经网络集成方法, 增加了冰雹预测的准确率。美国国家灾害风暴实验室研究了冰雹的雷达和环境属性, 在此基础上Marzban等^[6]利用贝叶斯神经网络方法开发用于预测冰雹尺寸和预测冰雹大小分类的两个神经网络。Marzano等^[7]还通过贝叶斯方法, 使用C-Band双偏振雷达对水凝物相态进行分类, 区分水、冰和混合相的10种不同水凝物类别。

相对于国外, 国内在冰雹识别的宏观层面进行的研究较多, 对冰雹云的识别方法有很多, 如宏观识别法(即目测法)、气象要素变化识别法、卫星动力识别法、闪电定位系统识别法^[8]、雷达识别法等。在早期, 刘黎平等^[9-11]利用扩展边界条件法计算了不同相态的模型雹块对C波段雷达波的散射和衰减, 并在此基础上建立了C波段双线偏振雷达识别雹区的方法。曹俊武等^[12]利用美国KOUN雷达的一次强对流天气过程的观测资料, 讨论分析了双偏振雷达雹区识别, 给出一个模糊逻辑算法冰雹识别模式, 该模式不仅可以反映出实际的冰雹区位置, 还可以对其分类, 识别的结果比较符合实际的天气过程。苏德斌等^[13-14]利用北京地区X波段双线偏振雷达资料, 建立了X波段双线性偏振雷达冰雹识别参量H_DR, 当H_DR > 0时表示有冰雹, H_DR < 0时表示无冰雹。该判别与地面降雹情况基本一致, 但是存在一定的虚警率。

双偏振雷达已成为天气雷达的发展趋势, 其在降水估测、粒子的相态识别等方面具有明显的优势^[15-18]。至2019年, 广东省新建或者是改建的9部新一代天气雷达, 全部为S波段双偏振雷达。本文利用广东省S波段双偏振雷达2019年观测到的冰雹数据, 研究贝叶斯方法对华南地区的冰雹识别效果。

2 冰雹数据集的构建 2.1 资料来源

为了加强天气监测, 全国各地正逐步推进新一代天气雷达的偏振升级改造。CINRAD\SA雷达升级过程中, 吸收了美国WSR-88D偏振升级的研究成果, 采用双发双收的技术体制, 充分应用S波段单偏振雷达在线自动定标及补偿功能, 保证对双线偏振各参量的探测精度, 探测参量包括反射率(Z)、速度(V)、谱宽(W)、差分反射率(Z_DR)、差分传播相位常数(K_DP)、差分传播相移(Φ_DP)、偏振雷达相关系数(CC)。且依据反射率、差分反射率、零滞后相关系数等双线偏振探测量, 采用模糊逻辑算法, 双线偏振雷达可以初步探测目标的降水类型(雨、大雨、冰雹、雪、非气象目标物等)。

用于冰雹分析的雷达数据来自广东省S波段双偏振雷达(S-POL)2019年观测数据。根据广东2019年地面冰雹观测记录, 如降雹时间、地点、降雹直径等(表 1), 采取人机交互的方法, 挑选出降雹点上空的雷达基数据, 对雷达数据进行质量控制^[19]与处理分析, 得出冰雹过程中的各个参量的分布特征。表 1中零度层高度为清远站08或20时(北京时, 下同)探空资料得到。

2.2 S波段偏振雷达冰雹特征

在我国南方地区, 热力、水汽条件比较好, 对流旺盛, 过冷水充裕, 往往在孤立单体, 或者多单体对流云中, 能够发展出冰雹。根据2019年广东S-POL观测数据, 冰雹在雷达回波上表现为: PPI图像上中心回波强度在60 dBZ以上, 中心区域面积较小, 云体发展很高, 在较高的仰角层上仍能看到很强的回波; 同时该地区冰雹的Z_DR值与相关系数CC值变化较大, Z_DR值分布在-2~4 dB区间, CC值在0.98附近。在剖面图上有较明显的回波墙, 但弱回波区, 悬垂回波等回波特征没有超级单体中形成的冰雹那么明显, 超级单体形成的雹灾通常具有较大的有界弱回波区^[20-21]。

以2019年4月20日11:50左右, 江门恩平恩城街道观测到10~20 mm的冰雹为例, 图 1、图 2是广东阳江S-POL雷达在11:54时刻, 仰角0.5 °的PPI和强回波区的剖面图像中的Z、Z_DR、CC参量。在图 1中, 冰雹天气具有明显的局地性特点, 强中心范围较小, 回波强度比较强, 中心强度普遍达到了50 dBZ左右, 最强的区域在55 dBZ以上。Z_DR值变化较大, 在冰雹区域外围, 值分布在0 dB左右, 到回波中心Z_DR值在1~2 dB, 回波最强的区域Z_DR值达到4 dB以上, 范围较小, 可能是区域外围是小雨区, 而中心区域内为雨夹雹。CC观测显示, 整个冰雹区域相关系数均在0.98左右, 在回波中心内CC值较小的部分值也达到了0.97, 与纯的冰雹CC值存在一定差异。图 2中显示, 此次天气过程云体发展高度非常高, 其回波结构具有高悬的强反射率因子核心, 45 dBZ回波超过10 km, 最大回波强度在65 dBZ以上, 强回波略微倾斜, 但是悬垂回波不明显, 在融化层(当日零度层高度为4 859 m)以上, 强回波对应Z_DR值为0 dB左右。在剖面图中椭圆标注处, 位于强回波区, 且Z_DR和CC较小, 有比较典型的冰雹偏振特征, 为纯冰雹区域。同时在融化层以下, Z_DR剖面图标注区域有Z_DR柱存在, 表明该区域为大粒区, 或者是冰雹融化区域。冰雹在下落的过程中快速融化, 在其表面形成水膜, 当水膜厚度达到冰雹直径的30%时, 冰雹的介电性质将与雨滴相近, 形成大雨滴的偏振特性。在CC的剖面图中, 整个强回波区的CC值偏小, 特别是在2~4 km高度, 有一个明显的低值区, 这也说明冰雹在长达2 km的降落过程中, 不断融化, 冰雹、水膜冰雹、雨滴并存, 导致这个高度区间出现较小的CC值。

2.3 构建冰雹数据集

本文旨在利用贝叶斯方法对华南冰雹的识别效果进行改进, 在利用贝叶斯方法进行冰雹识别的过程中, 需要构建冰雹的数据集。

因此, 根据表 1中记录的冰雹发生时间与地点, 查找并输出对应时空的雹云各个仰角的雷达数据。考虑到广东地区纯粹的冰雹很少, 基本上都是雨夹雹的特点, 输出的数据集限定为: Z值在45 dBZ以上, Z_DR值在(-1, 4)附近, CC集中在(0.85~0.98), 剖面图上具有明显的冰雹回波墙结构, 但不一定有穹窿状的弱回波区(图 3)。图 3为广东地区的三次雹云单体剖面图(对应表 1的序号分别为21、22、30), 在图 3a1、3b1、3c1中标注区域为冰雹区域, 从偏振量上, 冰雹区域强度很大, Z_DR、CC值也明显偏小; 在图 3a2、3b2、3c2中, CC剖面上有三个小值区, 标注区左边40 km处判断为信噪比小造成, 标注区为真正的冰雹区域, 而在标注区右边的小值区是因为回波经过冰雹区域衰减, 导致后面Z、CC都在一定程度上减小。虽然S波段衰减较小, 但是大冰雹对其仍有明显影响; 图 3a3、3b3、3c3中标注区域冰雹特征更为明显, 强回波, 弱Z_DR、CC值。

此种方式共确定133 263组(Z、Z_DR、CC)冰雹数据, 作为正样本; 另外, 从其他没有冰雹记录的体扫资料中, 随机取反射率强度大于45 dBZ的降水数据作为非冰雹数据, 共选取了140 000组(Z、Z_DR、CC)数据, 作为负样本, 正、负样本数据的比例为1:1.05。将以上的冰雹和非冰雹数据作为训练数据, 利用贝叶斯方法建立冰雹识别模型。

3 基于贝叶斯方法改进冰雹识别步骤 3.1 贝叶斯原理

贝叶斯定理是概率论中的一个定理, 描述在已知的一些条件下, 某事件发生的机率。通常, 事件A在事件B已发生的条件下发生的概率, 与事件B在事件A已发生的条件下发生的概率是不一样的。然而, 两者是有确定的关系的, 贝叶斯定理就是这种关系的陈述。如下所示:

$ P(\mathrm{~A} \mid \mathrm{B})=\frac{P(\mathrm{~B} \mid \mathrm{A}) P(\mathrm{~A})}{P(\mathrm{~B})} $

(1)

其中A以及B为随机事件, 且P(B)不为零。P(A|B)是指在事件B发生的情况下事件A发生的概率。在贝叶斯定理中, 每个名词都有约定俗成的名称:

P(A|B)是已知B发生后, A的条件概率。也由于得知B的取值而被称作A的后验概率。

P(A)是A的先验概率(或称边缘概率)。之所以称作“先验”是因为它不考虑任何B方面的因素。

P(B|A)是已知A发生后, B的条件概率。也由于得知A的取值而被称作B的后验概率。

P(B)是B的先验概率。

对于变量有两个以上的情况, 贝叶斯定理亦成立。例如:

$ P(\mathrm{~A} \mid \mathrm{B}, \mathrm{C})=\frac{P(\mathrm{~A}) P(\mathrm{~B} \mid \mathrm{A}) P(\mathrm{C} \mid \mathrm{A}, \mathrm{B})}{P(\mathrm{~B}) P(\mathrm{C} \mid \mathrm{B})} $

(2)

3.2 贝叶斯识别冰雹原理

如果事先知道冰雹区域内Z、Z_DR和CC概率分布, 则可以利用贝叶斯方法开展冰雹识别。基于贝叶斯的冰雹识别算法描述如下: 将雷达回波分成两类, C=(H, N), 其中C表示类别, 冰雹区域回波用H表示, 非冰雹回波用N表示, 三个偏振参量Z、Z_DR和CC结合起来形成雷达观测差别向量 $y=(\hat Z, {\hat Z_{{\rm{DR}}}}, {\rm{\hat CC)}}$, 要判断观测向量y是否属于H, 只有在p(H|y) > p(N|y)时, 属于H, p是概率密度, 根据贝叶斯理论:

$ p\left(C_{i} \mid y\right)=\frac{p\left(y \mid C_{i}\right) p\left(C_{i}\right)}{p(y)} $

(3)

其中C_i=H, N。P(y)=K是观测差别因子y的概率, 假定它对H和N分类概率是一样的, 即P(H)=P(N)= 1/2, 因此P(C_i|y) 和P(y|C_i)P(C_i)成正比, 则公式(3)转换成:

$ P\left(C_{i} \mid y\right)=\frac{1}{2 K} P\left(y \mid C_{i}\right) $

(4)

基于简单的贝叶斯判断中分类之间独立的假设, 条件概率密度可以分解为:

$ P\left(y \mid C_{i}\right)=P\left(\hat{Z}, \hat{Z}_{D R}, \hat{C} C \mid C_{i}\right)=P\left(\hat{Z} \mid C_{i}\right) P\left(\hat{Z}_{D R} \mid C_{i}\right) P\left(\hat{C} C \mid C_{i}\right) $

(5)

利用贝叶斯方法识别冰雹具体流程如图 4所示:

3.3 先验概率分布的获取

对表 1中冰雹记录区域上空的所有S-POL的Z、Z_DR、CC分布作散点图(图 5), 散点的透明度表示了粒子的集中程度。在冰雹区域内, 反射率Z从低层至高层, 分布都较为均匀, 表明在整个云体当中, 反射率都比较大; Z值主要分布在45~65 dBZ区间, 其中65 dBZ以上的云体部分主要在10 km以下, 表明其冰雹核心区域是位于中低层。冰雹区域中Z_DR集中分布在0 dB附近, 具有一定的波动, 在6 km以下的区域其波动比较明显, 说明冰雹在高层更加纯粹, 基本维持着单一的相态; 在随高度降低的过程中, 由于气温变化较大, 地面气温比较高, 导致冰雹不断的融化, 形成了水膜、小冰雹完全融化形成雨夹雹, 影响了Z_DR的分布。从CC分布图中可以看到, 相关系数越接近1, 分布越集中, 与Z、Z_DR分布较为类似, 在10 km以上的云体中比较集中, 而在高度下降的过程中, 分布趋于散乱, 在融化层附近(零度层基本位于4~5 km), 粒子的相态发生了变化, 由于雨水和冰相粒子的混合, 导致了相关系数的变化, CC值快速减小, 小于0.9的区域变大; 在0~2 km间, CC的波动更加明显, 最小值在0.1附近。

利用构建的冰雹数据集, 分析得出Z、Z_DR和CC独立概率分布(图 6)。图中蓝色线表示冰雹, 红色线表示非冰雹, 从图中可以看出, 冰雹和非冰雹的Z值差异在于, 冰雹的分布集中在50 dBZ以上的区域, 概率的峰值在0.04附近; 而非冰雹的分布趋向于Z值较小的区域, 在45 dBZ左右的概率为0.07~0.08。而在Z_DR分布中, 相对于非冰雹的Z_DR分布, 冰雹的分布更为分散, 跨度较大, 概率的峰值区在1 dB附近, 为0.015, 而非冰雹的Z_DR分布比较集中, 概率的大值在1~2 dB左右, 为0.035。冰雹区域的CC值比较集中在0.98, 峰值区概率在0.10左右; 而非冰雹的CC值更加偏向0.99~1间, 概率的峰值达到了0.2以上。

4 个例分析

选取2019年4月11日16时左右广州增城和4月12日09时左右江门的两次冰雹过程, 作为冰雹识别效果的检验, 这两次冰雹个例的数据没有加入建模的训练集。增城地区冰雹过程的PPI和剖面图像如图 7、图 8所示, 在图 7冰雹PPI图像中, 回波的区域较小, 强度很强, 中心的回波强度普遍在50 dBZ以上, 最强的区域达到了65 dBZ左右; 在图 7b中Z_DR值波动也比较大, 平均处于1.5 dB左右, 部分地区在0.5 dB左右, 而有的区域却达到了4 dB以上; 相关系数CC的值较小, 分布在0.90附近, 在回波中心边缘部分也有值处于0.98附近。在图 8剖面图像中, 可以看到冰雹的结构特征, 标注区域偏振量体现出了明显的反射率大值区与Z_DR、CC的弱值区。在冰雹中心区域, 回波的强度从低层至高层都比较强, 基本在50 dBZ以上, 在2~8 km的区域, 有65 dBZ以上的大值中心, 为冰雹的集中区域; 在图 8b中, Z_DR值在4 km以上的高空都比较小, 在0.5 dB左右, 而在4 km以下的区域中, Z_DR值出现一个大值区, 为冰雹经过融化层(当日零度层为4 183 m), 表面融化导致。在图 8c相关系数图中也体现出这一特征, 在融化层以上的区域, CC值较大, 均在0.96以上, 在融化层后, CC值快速减小, 部分区域处于0.90附近, 体现出了粒子相态颇为驳杂的特点。

图 9、图 10分别是2019年4月12日江门台山北陡镇附近一次冰雹过程的PPI和相应的剖面图, 与增城的过程相似, PPI图中回波强度较大, 普遍在50 dBZ以上, Z_DR值从外围至中心缓慢增大, 云体周边值在0 dB, 中心最大值达到了4.0 dB以上, 相关系数图中显示回波中心CC值在0.97左右。剖面图中, 可以看到云体发展较高, 回波最强在65 dBZ以上, 但没有增城过程的明显, 云体稍稍倾斜, 悬垂不明显; 当日融化层位于4 km左右(零度层高度为4 149 m), 在椭圆标注区右侧, 经过融化层时, Z_DR值突然变大, 在云体的下部达到了2.5 dB以上, 同时CC值减小到0.9附近, 然后在云体下部增加达到0.98左右, 为冰雹区域; 而标注区域, 出现CC的小值区, 认为是处于下沉气流区域, 回波墙边缘, 有少量的冰雹越过主上升气流, 在这个区域并没有液态水凝物, 或者是存在少量的雨滴, 强的反射率主要是由于少量冰雹引起的。

为了验证贝叶斯方法识别冰雹的效果, 选取2019年4月11日16:06和4月12日11:54的雷达数据进行冰雹识别, 并使用美国WSR-88D的冰雹识别算法与贝叶斯方法的识别结果进行综合分析。WSR-88D是对风暴系列算法识别出的每一风暴单体分析、计算其结构而进行冰雹探测的^{[2, 22]}。

两种方法具体的识别结果如图 11所示, 图 11为4月11日的识别结果, 其中第1列是16:06的雷达反射率强度图, 第2列是WSR-88D的识别结果, 第4列是贝叶斯方法的识别结果, 从上至下分别是仰角依次为0.50、1.50、2.50、3.50、4.50、6.00、10.00;图中第2列和第3列中的红色区域是识别出冰雹的区域。从反射率强度图中可以看到, 该次冰雹天气过程中心强度达到了60 dBZ以上, 并且从低层至高层都维持着很高的强度; 且回波强中心的范围较为集中, 面积相对于整片雷达回波来说较小, 反映了冰雹的局地性特点。在第2列WSR-88D的识别结果中可以看到, 算法识别出冰雹的区域基本与反射率强度图中的强回波中心位置重合, 表明该算法能有效识别出冰雹, 对于位置和区域的判断较为准确。在低层, 预测的冰雹范围和反射率的强中心范围大小不匹配, 可能原因是冰雹在下落过程中, 逐渐融化, 表面形成水膜, 导致Z_DR等偏振特性变化, 在低层的识别的冰雹范围减小; 在随高度的变化过程中, WSR-88D的识别算法预测出来的冰雹范围与60 dBZ强回波中心的范围趋于一致。在第3列的贝叶斯方法识别效果图中, 其预测冰雹的位置与强回波中心保持一致, 范围与50 dBZ以上的回波区域大小相近; 从低层到高层, 贝叶斯方法预测冰雹范围的形状也与强回波区域的形状相似。对比贝叶斯方法和WSR-88D冰雹识别算法对该次冰雹过程的识别效果, 可以看到贝叶斯方法识别出的区域相对较大, 而且与WSR-88D冰雹识别算法不同的是, 贝叶斯方法在低层识别的区域就比较大, 造成这一明显区别的原因是WSR-88D识别的是纯的冰雹, 而在华南地区, 由于中低层大气温度较高, 冰雹在下落的过程中不断融化, 高空的冰雹, 到了近地面, 小冰雹就变成了雨滴, 大冰雹外包水膜, 从而形成雨夹雹; 相对而言, 贝叶斯方法识别的是雨夹雹, 所以从整层来说, 贝叶斯方法识别出的范围比WSR-88D冰雹识别算法识别出的要大; 此外在构建冰雹数据集时, Z_H、Z_DR、CC范围的设定, 无法完全排除液态水的存在, 这在一定程度上也导致了贝叶斯方法识别出的范围较大。

图 12为4月12日的识别结果, 从反射率强度图中可以看到, 该次冰雹天气过程中心强度达到了55 dBZ以上, 并且从低层至高层都维持着很高的强度, 在高层, 云体中心的回波强度达到65 dBZ以上; 且回波强中心的范围较为集中, 面积相对于整片雷达回波来说较小, 反映了冰雹的局地性特点。在第2列WSR-88D的识别结果中可以看到, 算法识别出冰雹的位置与回波的位置接近, 且识别出的冰雹区域呈现出线状, 从南向北排列, 面积较小。在图 12b2中, WSR-88D算法识别出右上角也有冰雹出现而贝叶斯方法没有识别。在高层, 图 12b6、12b7中, 因回波出现分裂的情况, WSR-88D识别结果也出现了南北两个区域, 两个区域与回波位置、形状基本保持一致。在第3列贝叶斯方法的识别结果中, 和WSR-88D识别结果相比, 贝叶斯方法识别的结果区域更加集中, 且范围较大, 与反射率图中回波中心的位置和区域较为吻合, 从低层至高层都有此特征; 在第六层和第七层, 贝叶斯方法的识别结果仍是一个区域, 与65 dBZ以上的回波中心保持一致。

结合两次个例的预测结果, 与WSR-88D识别算法相比, 由于贝叶斯方法识别的是雨夹雹, 识别的区域更大更集中, 但基本的位置区域两种方法预测都较为准确。

对于华南地区的冰雹识别来说, 由于其明显的地域特征, WSR-88D冰雹识别算法虽然效果不差, 但是仍然有不足之处。在该区域, 纯冰雹的现象是比较难以出现的, 而外包水膜的冰雹同样能形成较大的雹灾, 因此使用贝叶斯方法来识别冰雹是有必要的。

5 讨论

本文通过对广东地区S波段双线偏振雷达2019年的冰雹观测数据进行分析, 首先在地面观测记录中有冰雹过程, 其次是根据观测记录的时间和地点挑选出雷达数据, 如果雷达数据的PPI图像具有冰雹特征, 并且剖面图中有明显的回波墙等冰雹结构, 将此过程认定为真实的冰雹过程, 可以用于构建冰雹数据集。在此基础上给出S波段偏振雷达在冰雹区域和非冰雹区域反射率Z、差分反射率Z_DR、和相关系数CC的先验独立概率密度分布。利用贝叶斯方法, 对2019年4月11日16:06和4月12日11:54的冰雹体扫数据, 比较和分析了(Z、Z_DR、CC)独立概率密度分布条件下冰雹的识别效果, 并得出独立概率密度分布能有效地识别雨夹雹的天气现象。

个例分析表明, 贝叶斯方法能有效地识别出冰雹, 且识别的准确率较高。在与广东现在的冰雹识别算法——WSR-88D冰雹识别算法的对比分析过程中, 贝叶斯方法识别的位置、范围大小与实际情况较为吻合, 且识别出雨夹雹天气现象, 是有利于提高冰雹的识别效果的。