CN103529492A - 基于多普勒雷达反射率图像的风暴体位置及形态预测方法 - Google Patents

Abstract

本发明公开了一种基于多普勒雷达反射率图像的风暴体位置及形态预测方法，涉及气象学领域，所述方法包括：利用前后两时刻相同反射率强度阈值的预测层，采用双向膨胀预测算法得到F_ti+1时刻风暴体的形态预测结果；以25dBZ预测层在F_ti+1时刻风暴体的形态预测结果为基底，依次把其他6层的F_ti+1时刻风暴体的形态预测结果叠加，并将t_i时刻的保留层进行最终叠加，获取F_ti+1时刻风暴体的形态预测图像；根据位置变化和方向变化对风暴体进行方位预测，获取F_ti+1时刻风暴体的方向角度γ和位置坐标C，结合F_ti+1时刻风暴体的形态预测图像得到最终预测图。本发明实现了对风暴体的方位与结构形态的预测，能够更加准确地判断出强对流天气的类型。

Description

基于多普勒雷达反射率图像的风暴体位置及形态预测方法

技术领域

本发明涉及气象领域，特别涉及一种基于多普勒雷达反射率图像的风暴体位置及形态预测方法。

背景技术

新一代天气雷达可以给出降水粒子或降雹粒子的实时回波强度（即反射率因子Z）、径向速度（V）和速度谱宽（W），并以多仰角图像的方式提供天气信息。利用雷达资料监测和临近预报强对流天气的一个主要方法是基于对雷达回波准确跟踪及合理预测的。目前国际上最常用的跟踪法有质心跟踪法和交叉相关跟踪法。以质心跟踪法为基础发展起来的代表算法有雷暴的识别、跟踪、分析和临近预报(TITAN)、风暴体识别和跟踪算法(SCIT)；交叉相关跟踪算法是利用矩形网格在搜索半径内寻找最大相关系数以跟踪雷达回波(TREC)。

以上述算法为基础发展出多种临近预报系统，例如，美国气象开发实验室MDL（Meteorological Development Lab）开发的SCAN（The System for Convectional Analysis andNow-casting）预报系统，用于风暴以及强风暴0-3h预报；美国NSSL（National Severe StormsLaboratory）开发的WDSS-II（Warning Decision Support System-Integrated Information）采用SCIT、冰雹监测、中尺度气旋检测、风暴多尺度识别等算法来监测和预警强对流天气,并进行0-1h强对流天气预报；美国国家大气科学研究中心NCAR（the National Center for AtmosphericResearch）开发的临近预报系统ANC（Auto-Nowcaster）采用TITAN、TREC算法跟踪风暴，并综合雷达、卫星、探空、闪电等多种传感器的资料进行风暴的0-2h临近预报。

发明人在实现本发明的过程中，发现现有技术中至少存在以下不足：所有算法都倾向于在风暴体识别与跟踪的基础上外推出风暴位置，而不强调风暴体内部的形态及结构预测，而形态及结构信息恰恰对判断强对流天气的类型十分有用。如果能够根据风暴的演变痕迹合理地预测出其回波图像的结构和形态，将会更有助于强对流风暴的临近预报。

发明内容

本发明提供了一种基于多普勒雷达反射率图像的风暴体位置及形态预测方法，本发明实现了对风暴体的方位与结构形态的预测，能够更加准确地判断出强对流天气的类型，详见下文描述：

一种基于多普勒雷达反射率图像的风暴体位置及形态预测方法，所述方法包括以下步骤：

（1）利用风暴体识别技术，对当前时刻的多普勒雷达反射率全局图进行图像分割得到风暴体关注区；

（2）获取风暴体在上一时刻时所处的位置，根据当前时刻和上一时刻的位置变化Δd，通过椭圆拟合获取风暴体的方向变化Δθ；

（3）旋转上一时刻风暴体使之与当前时刻风暴体方向相同，对上一时刻和当前时刻风暴体反射率图像进行分层，分别获取7对以不同反射率强度为阈值的预测层和1对保留层；

（4）利用前后两时刻相同反射率强度阈值的预测层，采用双向膨胀预测算法得到F_ti+1时刻风暴体的形态预测结果；

（5）以25dBZ预测层在F_ti+1时刻风暴体的形态预测结果为基底，依次把其他6层的F_ti+1时刻风暴体的形态预测结果叠加，并将t_i时刻的保留层进行最终叠加，获取F_ti+1时刻风暴体的形态预测图像；

（6）根据位置变化和方向变化对风暴体进行方位预测，获取F_ti+1时刻风暴体的方向角度γ和位置坐标C，结合F_ti+1时刻风暴体的形态预测图像得到最终预测图；

（7）将下一个时刻视为当前时刻，重复步骤（1）至步骤（6），实现对风暴体的持续预测，获取风暴体在更远时刻的预测图像。

所述获取风暴体在上一时刻时所处的位置，根据当前时刻和上一时刻的位置变化Δd，通过椭圆拟合获取风暴体的方向变化Δθ的步骤具体为：

对于在风暴体提取中得到的当前时刻t_i的一个风暴体Ω，用其外包矩形区域的左上角坐标标记Ω于t_i时刻的位置；在上一时刻反射率图中匹配风暴体Ω的所在区域，并将匹配到的区域外包矩形左上角坐标记为Ω于t_i-1时刻的位置，进而得到风暴体Ω的位置变化信息Δd；对于t_i时刻和t_i-1时刻得到的风暴体40dBZ以上核区进行椭圆拟合，通过两时刻椭圆方向角变化得到风暴体Ω的旋转角度Δθ。

所述利用前后两时刻相同反射率强度阈值的预测层，采用双向膨胀预测算法得到F_ti+1时刻风暴体的形态预测结果的步骤具体为：

设p为图像中任意一点，定义当前时刻风暴体一个预测层T_i中风暴区为F_ti，对应前一时刻同阈值下的预测层风暴区为F_ti-1，并将F_ti-1未包含且F_ti包含的区域称为新生区，记为

F_ti-1包含且F_ti未包含的区域称为消散区，记为 $P_d=\left\{p\right|p\∈F_{\mathrm{ti}-1}\}\∩\{p|p\∉F_{\mathrm{ti}}\};$

1）对于同阈值下前后两个时刻的风暴体单色预测层F_ti和F_ti-1，获得新生区P_e和消散区P_d；

2）对新生区P_e和消散区P_d进行膨胀操作分别获得膨胀后新生区S_e和膨胀后消散区S_d；

3）根据 $F_{\mathrm{ti}+1}=\left(\right\{p|p\∈F_{\mathrm{ti}}\}\∪\left\{p\right|p\∈S_e\left\}\right)\∩\left\{p\right|p\∉S_d\}$ 的定义，获得风暴体预测层的预测结果F_ti+1。

所述根据位置变化和方向变化对风暴体进行方位预测，获取F_ti+1时刻风暴体的方向角度γ和位置坐标的步骤具体为：

设点A为风暴体于t_i-1时刻外包矩形区域的左上角坐标，α为此时刻风暴体的方向角，点B为风暴体于t_i时刻外包矩形区域左上角坐标，β为此时刻风暴体的方向角；由

Δθ=β-α，即将该风暴体方位外推至t_i+1时刻，获得坐标C=B+Δd，方向角γ=β+Δθ。

本发明提供的技术方案的有益效果是：本发明在风暴体方位外推的基础上，深入到云体的内部，关注了风暴体内部变化的趋势，同时考虑延续风暴体的旋转变化特点，填补了临近预报中内部结构形态不能有效预测的一项空白，此方法同时能够适应风暴体分裂、合并等客观演变现象，有益于加强风暴体的预报准确性。

附图说明

图1为本发明提供的一种基于多普勒雷达反射率图像的风暴体的位置及形态预测方法的流程图；

图2为本发明提供的风暴体提取示意图；

图3为本发明提供的风暴体分层分解处理示意图；

图4为本发明提供的双相膨胀算法示意图；

图5为本发明提供的膨胀失真效果示意图；

图6为本发明提供的结构元素自适应选择示例图；

图7为本发明提供的分层预测叠加原理示意图；

图8为本发明提供的线性外推模型示例图；

图9为本发明提供的外推6～18分钟的测试示例的示意图；

图10为本发明提供的与TREC算法对比数据统计图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地详细描述。

就预报而言，人们更关心的问题是怎样在多普勒雷达图像上准确预测强对流风暴体的移动和变化趋势，问题包括风暴体方位变化的外推和结构变化的预测等。冰雹、飓线、龙卷风和雷暴等中小尺度灾害性天气具有发生演变快、移动异常迅速和破坏力巨大等特点，对其进行准确地分析和预测是气象预报中亟待解决的难点问题，同时也是灾害性天气预测中迫切需要解决的问题。为了提高风暴体预测的准确性，降低误报率，本发明实施例提供了一种基于多普勒雷达反射率图像的风暴体位置及形态预测方法，参见图1，详见下文描述：

101：利用风暴体识别技术，对当前时刻的多普勒雷达反射率全局图进行图像分割得到风暴体关注区；

其中，该步骤具体为：在多普勒雷达的反射率图中，一般采用16种伪彩色分别表示16个等级的反射率强度，用以反映探测目标内部降水粒子的尺度和密度。一般来讲，在反射率图像中，风暴体总是由一个或几个以高反射率为中心、由中心向外反射率强度逐级降低的子区域构成，这些子区域内部的高反射率部分称为核区。以此为依据从整幅反射率图中提取出风暴体所在的矩形区域作为风暴体关注区，参见图2a和图2b，以图2b为例进行说明：矩形框所包围的区域为所识别出的风暴体关注区，此关注区具有典型的单体结构特点，即有一个高反射率区域中心、由中心区域向外反射率强度逐级降低。

102：获取风暴体在上一时刻时所处的位置，根据当前时刻和上一时刻的位置变化Δd，通过椭圆拟合获取风暴体的方向变化Δθ；

其中，该步骤具体为：对于在风暴体提取中得到的当前时刻t_i的一个风暴体Ω，用其外包矩形区域的左上角坐标标记Ω于t_i时刻的位置；利用交叉相关跟踪算法在上一时刻反射率图中匹配风暴体Ω的所在区域，并将匹配到的区域外包矩形左上角坐标记为Ω于t_i-1时刻的位置，进而得到风暴体Ω的位置变化信息Δd；对于t_i时刻和t_i-1时刻得到的风暴体40dBZ以上核区进行椭圆拟合，通过两时刻椭圆方向角变化得到风暴体Ω的旋转角度Δθ（例如：当前时刻t_i的椭圆方向角为β，t_i-1时刻的椭圆方向角为α，旋转角度Δθ为β-α），用于对风暴体Ω的方位预测。

103：旋转上一时刻风暴体使之与当前时刻风暴体方向相同，对上一时刻和当前时刻风暴体反射率图像进行分层，分别获取7对以不同反射率强度为阈值的预测层和1对保留层；

其中，该步骤具体为：从-5dBZ、0dBZ、5dBZ、10dBZ、15dBZ、20dBZ、25dBZ、30dBZ、35dBZ、40dBZ、45dBZ、50dBZ、55dBZ、60dBZ、大于等于65dBZ等15个反射率量级上选择25dBZ到55dBZ七个强度作为分层分割阈值，再保留60dBZ及以上区域。即将风暴图像分为7个单色预测层和1个高反射率保留层。

每一个单色预测层T_k（k=25,30,35,40,45,50,55）的生成方法如下：对于多普勒雷达反射率图像中的像素，若反射率值R≥R_kdBZ，则令其在T_k中的取值R=R_kdBZ，否则令R为背景色，其中，R_k=25,30,35,40,45,50,55dBZ，参见图3。

高反射率区域的保留层的获得方法是令多普勒雷达反射率图像中反射率值小于60dBZ的像素置成黑色，参见图3。

104：利用前后两时刻相同反射率强度阈值的预测层，采用双向膨胀预测算法得到F_ti+1时刻风暴体的形态预测结果；

其中，该步骤具体为：

参见图4（a），设p为图像中任意一点，定义当前时刻风暴体一个预测层T_i中风暴区为F_ti，对应前一时刻同阈值下的预测层风暴区为F_ti-1，并将F_ti-1未包含且F_ti包含的区域称为新生区（打斜线的区域），记为 $P_e=\left\{p\right|p\∈F_{\mathrm{ti}}\}\∩\{p|p\∉F_{\mathrm{ti}-1}\};$ F_ti-1包含且F_ti未包含的区域称为消散区（灰色区域），记为 $P_d=\left\{p\right|p\∈F_{\mathrm{ti}-1}\}\∩\{p|p\∉F_{\mathrm{ti}}\}.$

在风暴体演变过程中，前后时刻发生的变化会在很大程度上延续到下一时刻，即新生区会在下一时刻继续向周边生长，消散区周边继续消散。为此，首先检出所有新生区和消散区，通过对新生区的膨胀操作，使新生区继续长大；通过对消散区的膨胀操作是风暴体在出现消散的地方继续回缩，即所谓双向膨胀策略，膨胀出的新生区反应了反射率在局部区域继续扩大的倾向，膨胀出的消散区反应了反射率在某一局部区域继续减弱的倾向。

参见图4（b），膨胀后新生区记作S_e（虚线所围的区域）；膨胀后消散区记作S_d（灰色点线所围的区域）。单色预测层预测结果为当前时刻图层经过新生区膨胀，并减掉膨胀后消散区的剩余区域。记为 $F_{\mathrm{ti}+1}=\left(\right\{p|p\∈F_{\mathrm{ti}}\}\∪\left\{p\right|p\∈S_e\left\}\right)\∩\left\{p\right|p\∉S_d\}.$ 双向膨胀的具体操作流程为

1）对于同阈值下前后两个时刻的风暴体单色预测层F_ti和F_ti-1，获得新生区P_e和消散区P_d。

2）对新生区P_e和消散区P_d进行膨胀操作分别获得膨胀后新生区S_e和膨胀后消散区S_d，膨胀过程根据联通域面积选择适当结构元素进行自适应膨胀。

3）参见图4（c），根据 $F_{\mathrm{ti}+1}=\left(\right\{p|p\∈F_{\mathrm{ti}}\}\∪\left\{p\right|p\∈S_e\left\}\right)\∩\left\{p\right|p\∉S_d\}$ 的定义，获得风暴体预测层的预测结果F_ti+1。

对新生区和消亡区进行膨胀操作时，膨胀结构元素的选取对于外推结果有着重要影响。本实施例中针对相对稳定的强对流风暴体模型选取各向同性结构元素，就图像膨胀处理中常用的3×3矩形结构元素会使膨胀结果出现凹形缺口缺失的情况，采用十字形结构元素，它能够在更好的保持原有结构基础上进行膨胀，参见图5。

结构元素确定之后，被膨胀的区域将向外同尺度扩展，这往往与实际情况不符；另外，当膨胀次数一定时，对于新生或消散区域延续的程度将取决于结构元素的尺度。为此，对于大小不同的新生区和消散区应当采用不同方法进行自适应性的膨胀处理。具体如下：

1）对于面积较小联通区域，其实际边界不会太过复杂，采用小十字结构元素进行膨胀；

2）对于面积较大的联通区域采用自适应的结构元素选择方案，即小十字结构元素和大十字结构元素适当结合的方法进行膨胀。首先对一个联通区域边缘取凸外包，把与凸外包重合的边缘标记为“粗膨胀区”，未与凸外包重合的边缘标记为“细膨胀区”。“粗膨胀区”标记的边缘采用大十字结构元素进行膨胀操作，以保证膨胀的尺度要求；“细膨胀区”标记的边缘采用小十字结构元素进行膨胀操作，以保证联通域边界处细节部分的形态，参见图6。

105：以25dBZ预测层在F_ti+1时刻风暴体的形态预测结果为基底，依次把其他6层的F_ti+1时刻风暴体的形态预测结果叠加，并将t_i时刻的保留层进行最终叠加，获取F_ti+1时刻风暴体的形态预测图像；

其中，该步骤具体为：以25dBZ单色预测层T₂₅在F_ti+1时刻预测结果为基底，依次用较高反射率单色预测层的预测结果对其进行覆盖，经6次覆盖后，再用t_i时刻大于等于60dBZ的区域（保留层）进行最终覆盖，得到风暴体整体预测图像，参见图7。

106：根据位置变化和方向变化对风暴体进行方位预测，获取F_ti+1时刻风暴体的方向角度γ和位置坐标C，结合F_ti+1时刻风暴体的形态预测图像得到最终预测图；

其中，该步骤具体为：

风暴体回溯过程中可以得到其相邻时刻位置以及方向角的变化，根据线性公式y=kx+b，可外推得到下一时刻风暴体的坐标位置及其旋转角度，以此从宏观上确定其方位。

设点A为风暴体于t_i-1时刻外包矩形区域的左上角坐标，α为此时刻风暴体的方向角。点B为风暴体于t_i时刻外包矩形区域左上角坐标，β为此时刻风暴体的方向角。由

Δθ=β-α，即将该风暴体方位外推至t_i+1时刻，获得坐标C=B+Δd，方向角γ=β+Δθ，参见图8。

107：将下一个时刻视为当前时刻，重复步骤101-步骤106，实现对风暴体的持续预测，获取风暴体在更远时刻的预测图像。

其中，该步骤具体为：为了获得更远时刻的外推图像，将此算法在时间上递进式使用，即先用T_i-1和T_i推出当前时刻6分钟之后的图像T_i+1，再用T_i和T_i+1推出当前时刻12分钟之后的图像T_i+2，以此类推。一个递进式预测的结果示例可参见图9。

下面以一个具体的试验来验证本发明实施例提供的一种基于多普勒雷达反射率图像的风暴体位置及形态预测方法的可行性，详见下文描述：

以雷达基数据为数据源，取来自强对流过程中的150个风暴体反射率图像做样本，对本方法进行差异性检验，并作出预测图与真实图间的相似度评估，评估指标设计如下。

悬垂度：

其中q_m为第m层回波对m-1层回波的局部悬垂，N为回波层数；

有效厚度：ΔH=max{H_R>45dBZ}-H_-20℃，即大于45dBZ的风暴体核的最大高度与-20度层高度H_-20℃之差；

高回波比：

其中N_≥50为反射率值大于等于50dBZ的像素点数，N_≥30为反射率值大于等于30dBZ的像素点数。

峰度：

其中Y_i是每个像素点的像素值，

是基于反射率值分布直方图的均值，s是标准差。

随机选取来自40个强对流风暴过程的150个风暴体进行预测并展开基于上述4个特征的统计检验。

对于各项指标，设预测数据与真实数据均来自服从正态分布的总体，并假设其均值相等，使用式下所示的t分布展开置信水平为(1-α)的显著性差异检验：

$t=\frac{(\stackrel{\‾}{x}-\stackrel{\‾}{y})}{\sqrt{\frac{(n_1-1)S_1^2+(n_2-1)S_2^2}{n_1+n_2-2}\·(\frac{1}{n_1}+\frac{1}{n_2})}}~t(n_1+n_2-2)$

式中，

分别为预测数据和真实数据的均值,

分别为对应的方差，n₁和n₂为预测数据和真实数据的样本个数。

取显著性水平α=0.05查表得t_0.05（150）=1.94。对于各个指标t分布统计结果如表1所示。

表1

由此可知，风暴体上述四个特征t<t_0.05（150），则接受假设，即预测图像特征与真实图像特征没有显著差异。

选择经典的TREC算法对风暴体提取后的矩形区域进行预测，得到预测图像并统计其与真实图像的相似度。本算法与TREC算法的统计结果参见图10。本算法与传统TREC算法预测准确性都会随时间下降，但就整体的相似度评估效果而言，本算法较TREC算法的准确度提高了15个百分点。

综上所述，本方法在对其进行分层分解的基础上，提出了旋转双向膨胀算法，实现了保留风暴体内部降水粒子分布等细节信息的预测，解决了保留风暴体内强度分布规律和变化趋势的问题，此方法同时能够适应风暴体分裂、合并等客观演变现象；并且经过实验验证通过本方法得到的预测结果具有较高的准确度。

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种基于多普勒雷达反射率图像的风暴体位置及形态预测方法，其特征在于，所述方法包括以下步骤：

（1）利用风暴体识别技术，对当前时刻的多普勒雷达反射率全局图进行图像分割得到风暴体关注区；

（2）获取风暴体在上一时刻时所处的位置，根据当前时刻和上一时刻的位置变化Δd，通过椭圆拟合获取风暴体的方向变化Δθ；

（3）旋转上一时刻风暴体使之与当前时刻风暴体方向相同，对上一时刻和当前时刻风暴体反射率图像进行分层，分别获取7对以不同反射率强度为阈值的预测层和1对保留层；

（4）利用前后两时刻相同反射率强度阈值的预测层，采用双向膨胀预测算法得到F_ti+1时刻风暴体的形态预测结果；

（5）以25dBZ预测层在F_ti+1时刻风暴体的形态预测结果为基底，依次把其他6层的F_ti+1时刻风暴体的形态预测结果叠加，并将t_i时刻的保留层进行最终叠加，获取F_ti+1时刻风暴体的形态预测图像；

（6）根据位置变化和方向变化对风暴体进行方位预测，获取F_ti+1时刻风暴体的方向角度γ和位置坐标C，结合F_ti+1时刻风暴体的形态预测图像得到最终预测图；

（7）将下一个时刻视为当前时刻，重复步骤（1）至步骤（6），实现对风暴体的持续预测，获取风暴体在更远时刻的预测图像。

2.根据权利要求1所述的一种基于多普勒雷达反射率图像的风暴体位置及形态预测方法，其特征在于，所述获取风暴体在上一时刻时所处的位置，根据当前时刻和上一时刻的位置变化Δd，通过椭圆拟合获取风暴体的方向变化Δθ的步骤具体为：

对于在风暴体提取中得到的当前时刻t_i的一个风暴体Ω，用其外包矩形区域的左上角坐标标记Ω于t_i时刻的位置；在上一时刻反射率图中匹配风暴体Ω的所在区域，并将匹配到的区域外包矩形左上角坐标记为Ω于t_i-1时刻的位置，进而得到风暴体Ω的位置变化信息Δd；对于t_i时刻和t_i-1时刻得到的风暴体40dBZ以上核区进行椭圆拟合，通过两时刻椭圆方向角变化得到风暴体Ω的旋转角度Δθ。

3.根据权利要求1所述的一种基于多普勒雷达反射率图像的风暴体位置及形态预测方法，其特征在于，所述利用前后两时刻相同反射率强度阈值的预测层，采用双向膨胀预测算法得到F_ti+1时刻风暴体的形态预测结果的步骤具体为：

设p为图像中任意一点，定义当前时刻风暴体一个预测层T_i中风暴区为F_ti，对应前一时刻同阈值下的预测层风暴区为F_ti-1，并将F_ti-1未包含且F_ti包含的区域称为新生区，记为F_ti-1包含且F_ti未包含的区域称为消散区，记为 $P_d=\left\{p\right|p\&Element;F_{\mathrm{ti}-1}\}\&cap;\{p|p\&NotElement;F_{\mathrm{ti}}\};$

1）对于同阈值下前后两个时刻的风暴体单色预测层F_ti和F_ti-1，获得新生区P_e和消散区P_d；

2）对新生区P_e和消散区P_d进行膨胀操作分别获得膨胀后新生区S_e和膨胀后消散区S_d；

3）根据 $F_{\mathrm{ti}+1}=\left(\right\{p|p\&Element;F_{\mathrm{ti}}\}\&cup;\left\{p\right|p\&Element;S_e\left\}\right)\&cap;\left\{p\right|p\&NotElement;S_d\}$ 的定义，获得风暴体预测层的预测结果F_ti+1。

4.根据权利要求1所述的一种基于多普勒雷达反射率图像的风暴体位置及形态预测方法，其特征在于，所述根据位置变化和方向变化对风暴体进行方位预测，获取F_ti+1时刻风暴体的方向角度γ和位置坐标的步骤具体为：

设点A为风暴体于t_i-1时刻外包矩形区域的左上角坐标，α为此时刻风暴体的方向角，点B为风暴体于t_i时刻外包矩形区域左上角坐标，β为此时刻风暴体的方向角；由

Δθ=β-α，即将该风暴体方位外推至t_i+1时刻，获得坐标C=B+Δd，方向角γ=β+Δθ。

Images