强对流单体结构及结构特征可视化时空剖析方法与流程

本发明涉及气象学领域，特别涉及一种基于光流法的单体跟踪算法和一种基于移动立方体算法的单体三维重建算法。

背景技术：

冰雹、强降水、雷雨大风等强对流天气，具有空间尺度小，生命周期短，突发性强，破坏力强等特点。同时，不同类型的强对流单体，其内部结构和粒子分布具备差异。例如冰雹云强回波的中心值、强回波核(45-65dbz)高度均高于强降水回波，易出现悬垂结构^[1-4]，体现为冰雹云一般于降雹前在高、强、垂三个方面有别于单纯性短时强降水；且降雹前vil值会不同程度持续跃增，降雹时vil值明显大于短时强降水^[4]。

现有雷达产品存在以下缺点和不足：

(1)雷达图像现采取的ppi显示方式是在一系列固定仰角上扫描360°进行采样的，即在某一仰角，雷达天线绕垂直轴z进行360°扫描^[5]，直观性不高。

(2)组合反射率图只显示了最高反射率，无法观察反射率在高度上的分布。

(3)剖面不易制作，人机交互体验差。

(4)不提供强对流单体的结构化特征，如描述强上升气流造成的回波悬挂结构的悬垂度，描述单体内核体反射率强度的高度分布、梯度分布等；更谈不上对它们演变规律的展示。

[参考文献]

[1]李德俊，唐仁茂，熊守全等；强冰雹和短时强降水天气雷达特征及临近预警[j]，气象，2011，37(4)：474-480。

[2]井喜，屠妮妮，井宇等；毛乌素沙漠东部边缘中β尺度暴雨和冰雹综合分析[j]，灾害学，2010，25(1)：73-77。

[3]李湘，张腾飞，胡娟等；云南冰雹灾害的多普勒雷达特征统计及预警指标[j]，灾害学，2015(2015年03)：88-93，98。

[4]张崇莉，向明堃，赖云华等；滇西北高原冰雹、短时强降水的多普勒雷达回波特征比较[j]，暴雨灾害，2011，30(1)：64-69。

[5]俞小鼎等；多普勒天气雷达原理与业务应用[m]，北京：气象出版社，2009。

[6]李聪.强冰雹自动识别技术与防雹作业决策方法研究[d].天津：天津大学,2014：22-31

[7]王萍，潘跃.基于显著性特征的大冰雹识别模型[j].物理学报,2013,62(6):515-524

[8]amburnsa,wolfpl.vildensityasahailindicator.weatherandforecasting,1997,12(3):473-478

[9]greenedr,clarkra.anindicatorofexplosivedevelopmentinseverestorms[c]7thconferenceofseverelocalstorms,missouri.1971.

[10]barronjl,fleetdj,beaucheminss.performanceofopticalflowtechniques.internationaljournalofcomputervision,1994,12(1):43-77

[11]lorensenwe,clinehe.marchingcubes:ahighresolution3dsurfaceconstructionalgorithm[c]acmsiggraphcomputergraphics.acm,1987,21(4):163-169.

[12]wangping,licong,zhangyuan.anadaptivesegmentationarithmeticadaptedtointertwinedirregularconvectivestormimages[c]，icml2013,vol.2,pp896-898

技术实现要素：

针对上述现有技术，本发明提供一种强对流单体结构及结构特征可视化时空剖析方法，可以用一组结构化特征的变化曲线表现强对流单体的特征变化趋势和变化规律；可以通过透明度的设置为用户展现对流单体更好的三维透视性，通过透明度参数的设置为用户展现特定反射率值域的三维构成与分布；并改进了垂直剖面的人机交互体验；通过交互式剖线角参量的设置，为用户展现反射率强度在不同方位的垂直面上的分布；增加了水平剖面的人机交互体验，通过交互式高度参量的设置，为用户展现在不同高度的水平面上的分布。

为了解决上述技术问题，本发明提出的一种强对流单体结构及结构特征可视化时空剖析方法，包括获得强对流单体特征的变化趋势，获得强对流单体三维图像及垂直和水平的剖面；

一、获得强对流单体特征的变化趋势的步骤如下：

步骤1、选取强对流单体特征

以雷达基数据为原始数据，求取获得有效厚度、高回波比、悬垂度、加权核高、峰度、组合反射率图核区梯度高值、垂直累计液态水含量和累计液态水密度；

同时，通过对雷达基数据进行插值得到规则三维场数据，根据规则三维场数据得到各反射率的高度分布，并用一组柱形图表示；

所述组合反射率图核区梯度高值用于量化单体核区的变化程度，其求取方法如下：

在组合反射率图中进行单体分割后，将单体中反射率r>40dbz的区域定义为单体核区，按照下述公式求取单体核区的中心p，

公式(1)中，pr>40dbz为所有反射率大于40dbz的点的坐标，rr>40dbz为对应反射率；

以单体核区的中心p为圆心作圆覆盖单体，并将该圆等分为八个扇形，选取三个相邻扇形使之核区面积之和最小，并求取这三个扇形的平均梯度g：

公式(2)中，∑g为三个扇形区域所有像素点的梯度之和，像素点的梯度只考虑大小，不考虑方向，num为像素点总数，该平均梯度g即为组合反射率图核区梯度高值；

步骤2、通过光流法实现前后体扫时刻的单体追踪：

利用光流法求得光流数：将前后体扫时刻的两帧组合反射率图进行单体分割得到前体扫时刻单体和后体扫时刻单体；然后输入该两帧组合反射率图，并通过光流法计算出光流场数据；遍历光流场数据点，对于横向速度或者纵向速度大于30的光流场数据点赋值为0；通过光流场数据计算得到前体扫时刻所有像素在后体扫时刻的流入位置；对比前后体扫时刻各单体的范围，计算出光流数；

光流数记为sum(a→b)，用以表示由前体扫时刻单体a流入后体扫时刻单体b的像素数；用光流数量化前体扫时刻单体a和后体扫时刻单体b的关联程度，光流数越高，前体扫时刻单体a和后体扫时刻单体b的相关联程度就越高，反之越低；

结合前后体扫时刻单体的面积s，综合考虑前后体扫时刻单体出现分裂或合并的情况，制定如下三个条件，若满足其一则认为前体扫时刻单体a和后体扫时刻单体b存在关联，

条件一、sum(a→b)>s(a)×50％；

条件二、sum(a→b)>s(b)×50％；

条件三、设sum(a→b)＝k，对前体扫时刻的l个单体流入后体扫时刻单体b的所有像素求和，记为sum，k/sum>50％，且k/s(b)>33％；

步骤3.根据步骤1选取的强对流单体特征和步骤2实现的前后体扫时刻的单体追踪，用一组曲线表示强对流单体特征的变化趋势。

进一步讲，本发明中，通过结构可视化剖析获得强对流单体三维图像及垂直和水平的剖面，步骤如下：

步骤1、构建强对流单体三维图像：首先，读取雷达基数据获得组合反射率图，并进行单体分割，插值获得强对流单体的规则三维场数据；然后，通过移动立方体算法提取等值面；最终，绘制显示等值面；

步骤2、单体垂直及水平剖面显示：根据步骤1中获得的规则三维场数据，插值获得垂直和水平剖面。

与现有技术相比，本发明的有益效果是：

(1)本发明通过基于光流法的单体跟踪算法实现单体跟踪，用一组结构化特征的变化曲线表现强对流单体的特征变化趋势和变化规律。

(2)本发明通过移动立方块算法实现单体的等值面提取后，采用面绘制的方法绘制单体的三维结构，通过透明度的设置为用户展现对流单体更好的三维透视性，通过透明度参数的设置为用户展现特定反射率值域的三维构成与分布。

(3)本发明通过连续调整角度，为用户展现对流单体垂直剖面上的构成及强度分布，改进了垂直剖面的人机交互体验；通过高度参量的设置，为用户展现任意高度上对流单体的水平构成和强度分布，增加了水平剖面的人机交互体验。

附图说明

图1是对前体扫时刻的组合反射率图进行单体分割，单体分割结果如图所示，矩形框为单体的最小外包矩形，矩形框左上角标号为对该单体的编号；

图2是对后体扫时刻的组合反射率图进行单体分割，单体分割结果如图所示，矩形框为单体的最小外包矩形，矩形框左上角标号为对该单体的编号；

图3(a)至图3(e)是部分单体特征的变化趋势；横轴表示时间，纵轴表示特征的量值，虚线箭头是对当前时刻的标识；其中：图3(a)是悬垂度变化趋势，图3(b)是有效厚度，图3(c)是加权核高，图3(d)是高回波比，图3(e)是各反射率在高度上的分布示意图；

图4表示单体中55dbz的等值面提取结果；

图5表示各反射率等值面提取结果，包括：30dbz，35dbz，40dbz，45dbz，50dbz，55dbz的等值面提取结果；

图6表示各反射率的等值面的绘制结果；

图7是单体整体三维图像加上水平剖面的结果。

具体实施方式

下面结合附图和具体实施例对本发明技术方案作进一步详细描述，所描述的具体实施例仅对本发明进行解释说明，并不用以限制本发明。

本发明提出的一种强对流单体结构及结构特征可视化时空剖析方法，主要包括获得强对流单体三维图像及垂直和水平的剖面、获得强对流单体特征的变化趋势。

一、获得强对流单体特征的变化趋势的步骤如下：

步骤1、选取强对流单体特征

以雷达基数据为原始数据，求取获得有效厚度、高回波比、悬垂度、加权核高、峰度、组合反射率图核区梯度高值、垂直累计液态水含量和累计液态水密度；其中，选取有效厚度^[7]、高回波比^[6]等特征表示单体的“高”、“强”特点，用悬垂度表示单体的强悬垂回波结构^[6]，用加权核高表示水成物固态率^[6]，用峰度表示回波分布的紧密程度^[7]，用组合反射率图核区梯度高值量化核区的变化程度，也选取常见的冰雹判别指标垂直累计液态水含量^[9]和累计液态水密度^[8]；

同时，通过对雷达基数据进行插值得到规则三维场数据，根据规则三维场数据得到各反射率的高度分布，并用一组柱形图表示；

所选特征的具体定义如下：

(1)有效厚度

定义单体核反射率值大于45dbz的顶高和-20℃层高度之差h-20为有效厚度：

δh＝max{hr>45dbz}-h-20

(2)高回波比

强冰雹反射率图通常出现大面积的高强反射率区域。求取反射率值在50dbz以上的区域面积sr>50占30dbz以上区域面积sr>30的比例ρ，并将其命名为高回波比；

(3)悬垂度

强对流云体中，强烈的上升气流位于弱回波区，其上方出现强悬垂回波，悬挂区域与弱回波区域呈正相关关系，本文将0℃层之上的弱回波体积定义为回波的悬垂度；

(4)垂直累计液态水含量，累计液态水密度与加权核高

垂直累计液态水含量vil通过将雷达反射率因子z转换成等价的液态水，积分得到垂直柱体内液态水总量：

vil＝∫mdhm＝3.44x10^-3z^4/7

设∫mdh的柱底面积为s，则进一步可以计算出累计液态水密度vild：

考虑到云中水成物粒子在零度层以上会越来越多的转变成固态，加入衰减系数k，形成加权垂直累积液态水含量密度：

式中，h0为回波底，h1是回波顶高；k取为高度的函数，且：

若h0<h<h0℃，则k＝1；

若h0℃<h<h-20℃，则k由1线性衰减到0.1；

若h1>h>h-20℃，则k由0.1线性衰减0。

设的柱底面积为s，得体积∑(h1-h0)*s，则加权垂直累积液态水含量密度为：

只考虑反射率大于45dbz的单体部分，且不区分核区内的反射率值的差异，则加权垂直累积液态水含量密度简化为下式，并将该特征命名为加权核高；

(5)峰度

峰度是一种突出直方图厚重尾部的四阶统计量，在一定程度上反映出单体低层反射率图所呈现的冰雹单体核偏向高梯度一侧、反射率分布紧致的特点：

其中，是均值，s是标准差；

(6)组合反射率图核区梯度高值

在组合反射率图中进行单体分割后，将单体中反射率大于40dbz的区域定义为单体核区，按照下述公式求取单体核区的中心p，

公式(1)中，pr>40dbz为所有反射率大于40dbz的点的坐标，rr>40dbz为对应反射率；

以单体核区的中心p作为圆心作圆覆盖单体，并将该圆等分为八个扇形，选取三个相邻扇形使之核区面积之和最小，并求取这三个扇形的平均梯度g：

公式(2)中，∑g为三个扇形区域所有像素点的梯度之和，像素点的梯度只考虑大小，不考虑方向，num为像素点总数，该平均梯度g即为组合反射率图核区梯度高值；

步骤2、通过光流法实现前后体扫时刻的单体追踪：

利用光流法^[10]求得光流数：将前后体扫时刻的两帧组合反射率图进行单体分割^[12]得到前体扫时刻单体和后体扫时刻单体；然后输入该两帧组合反射率图，并通过光流法计算出光流场数据；遍历光流场数据点，对于横向速度或者纵向速度大于30的光流场数据点赋值为0；通过光流场数据计算得到前体扫时刻所有像素在后体扫时刻的流入位置；对比前后体扫时刻各单体的范围，计算出光流数；

光流数记为sum(a→b)，用以表示由前体扫时刻单体a流入后体扫时刻单体b的像素数；用光流数量化前体扫时刻单体a和后体扫时刻单体b的关联程度，光流数越高，前体扫时刻单体a和后体扫时刻单体b的相关联程度就越高，反之越低；

结合前后体扫时刻单体的面积s，综合考虑前后体扫时刻单体出现分裂或合并的情况，制定如下三个条件，若满足其一则认为前体扫时刻单体a和后体扫时刻单体b存在关联，

条件一、sum(a→b)>s(a)×50％；

条件二、sum(a→b)>s(b)×50％；

条件三、设sum(a→b)＝k，对前体扫时刻的l个单体流入后体扫时刻单体b的所有像素求和，记为sum，k/sum>50％，且k/s(b)>33％；

消亡单体和新生单体均不会满足上述任一条件；

若前体扫时刻单体a既不分裂也不合并，且大部分被后体扫时刻单体b继承，此时对应条件一；

若后体扫时刻单体b大部分继承于前体扫时刻单体a，此时对应条件二；

若后体扫时刻单体b大部分来源于单体合并，合并部分的主要贡献者是前体扫时刻单体a，此时对应条件三；

步骤3.根据步骤1选取的强对流单体特征和步骤2实现的前后体扫时刻的单体追踪，用一组曲线表示强对流单体特征的变化趋势。

二、通过结构可视化剖析获得强对流单体三维图像及垂直和水平的剖面，步骤如下：

步骤1、构建强对流单体三维图像：首先，读取雷达基数据获得组合反射率图，并进行单体分割^[12]，插值获得强对流单体的规则三维场数据；然后，通过移动立方体算法^[11]提取等值面；最终，绘制显示等值面；

步骤2、单体纵向及横向显示：根据步骤1中获得的规则三维场数据，插值获得垂直和水平剖面。

综上，本发明提出的一种强对流单体结构及结构特征可视化时空剖析方法，其中，通过光流法实现前后体扫时刻的单体追踪，是在实现单体跟踪的基础上，选取一组强对流单体的结构化特征，用折线图的方式体现它们各自的变化趋势及变化规律，可为气象工作者的预报与科研提供方便。

研究材料：数据来源是天津市塘沽雷达站2015年8月18日的雷达数据。

图1和图2分别是前体扫时刻和后体扫时刻的单体分割结果，单体标号如图所示，根据光流法得到图1，图2的光流数如表1所示，前后体扫时刻单体面积如表2所示：

表1光流数表

表2单体面积表

根据上表得到sum(0→0)＝290>s(后体扫时刻0号单体)×50％＝339×50％＝169.5，故后体扫时刻0号单体与前体扫时刻0号单体存在联系，同理可判断后体扫时刻1，2号单体与前体扫时刻满0号单体存在联系，后体扫时刻3号单体来与前体扫时刻1号单体存在联系。以上结果表明前体扫时刻0号单体分裂为后体扫时刻的0，1，2号单体，后体扫时刻3号单体主要继承了前体扫时刻的1号单体。

在实现单体跟踪的基础上，计算每个体扫时刻所有单体的特征，用图表的形式展现单体的特征变化趋势，部分特征的变化趋势如图3(a)至图3(e)所示，横轴表示时间，纵轴表示特征的量值，虚线箭头是对当前时刻的标识；其中：图3(a)是悬垂度变化趋势，图3(b)是有效厚度，图3(c)是加权核高，图3(d)是高回波比，图3(e)是各反射率在高度上的分布示意图；

单体三维重建的实现主要分为以下三个步骤：

(1)在读取雷达基数据→单体分割→插值成规则三维场数据。

(2)通过移动立方体算法提取特定等值面，图4为55dbz等值面提取结果，可见其提取结果是由多个三角面构成的网格面。

(3)图5为30dbz，35dbz，40dbz，45dbz，50dbz，55dbz等值面提取结果，在提取等值面后，将相应的网格数据绘制为面，并设置合适的透明度。由于反射率高的等值面一般被反射率低的等值面包含，故反射率高的等值面的透明度低，反射率低的等值面反之，以此带来较好的透视性，图6为最后绘制结果，图7是单体三维图像加上水平剖面的效果图。

尽管上面结合附图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保护之内。