亚星游戏官网-www.yaxin868.com

专利名称：基于水平集的雷暴预测方法
技术领域：
本发明涉及的是ー种图像处理技术领域，具体是ー种基于水平集的雷暴预测方法。
背景技术：
经过对现有技术的检索发现，目前雷暴预测技术可以分成两类数值模型仿真和跟踪外推预测。数值模型仿真需要对大气物理过程的精确理解、完善反映物理过程的方程组以及高精度的数值计算，此种方法不在本专利领域讨论范围内。跟踪外推预测是根据以 往和当前时刻雷暴的发展情况外推下ー时刻的发展情況。但是现有技术不论使用何种外推方法，由于不能精确反映雷暴単体形状的变化，雷暴単体位置的预测精度较低。水平集方法把N维问题看做N+1维处理，使得2D形状的分裂与合并反映在3D上只是曲面的升降。其在图像跟踪领域中的应用开始于Osher和Sethian,经过Paragios、Mansouri等人的发展，已经比较成熟。如Mansouri等在《IEEE Trans. Pattern Analysisand Machine Intelligence》(2002,24(7) :947-961.)上发表的“Region trackingvialevel set PDEs without motion computation”。此方法目前广泛用于图像分割和图像跟踪领域，在图像预测特别是形变较大的雷达回波预测中未见。

发明内容
本发明针对现有技术存在的上述不足，提供一种基于水平集的雷暴预测方法，能够实现在目标尺度和外形发生较大变化条件下的精确图像预测。本发明是通过以下技术方案实现的，首先将3D雷达体扫回波数据投影到2D作为图像看待，根据3D雷达体扫回波数据计算相关特征，作为水平集发展的驱动カ推动雷暴单体形状演化，变形演化后的雷暴単体再根据跟踪外推平移到适当的位置。本发明的技术解决方案如下—种基于水平集的雷暴预测方法，其特点在于，该方法包括如下步骤第一歩、数据质量控制为了提高3D雷达体扫回波数据的质量，減少干扰，需要对其进行预处理，包括资料的缺测填补、滤波，地物、海浪杂波的识别，云和降水衰减订正、径向速度的处理和ー些特殊的气象处理算法等。量测数据预处理技术是对雷达数据进行正确处理的前提条件，有效的量测数据预处理可以降低雷达数据处理的计算量、提升后续算法的精度。第二步、坐标系转化处理在雷达系统中，目标量测所在的坐标系与数据处理所在坐标系经常是不一致的，此时就需要通过坐标转换技术将所有的数据信息格式统一到同一坐标系中，坐标系的选择直接影响数据精度和计算量的大小。第三步、雷暴识别处理3D雷达体扫回波数据在检测、跟踪、预报方面的算法复杂度都偏高，而对雷暴最感兴趣的是其在地球表面的覆盖情况，故将数据压缩到2D作为图像，方法是平面上的每ー个点都取所在位置所有仰角层中的最大值。
第四步、雷暴跟踪处理将当前时刻的雷暴目标与上一时刻的目标相关联，以WSR-88D型雷达为例，其一次体扫间隔大约为6分钟，即比较单前时刻的雷暴目标与6分钟前的雷暴目标。本发明采用随机粒子方法处理合并、分裂、新生和消散，简化了跟踪过程，易于处理分裂合并这ー传统难题。第五步、雷暴预测处理通过将雷暴単体历史位置的运动轨迹外推得到雷暴単体的未来位置。在该未来位置上，以雷暴単体本身的结构特征作为水平集演化的驱动カ之一。雷暴単体的大小和形状则根据计算出的驱动カ由水平集演化来确定。驱动力的计算使用经过数据质量控制后的3D数据，对雷暴目标周围一定距离的点进行发展趋势判断，继续加强的趋势作为正方向的驱动，减弱的趋势作为负方向的驱动。此步是本发明的核心，属于从未被公开过的新技术。所述的雷暴单体是指回波強度超过第一指定阈值Tz，且连续面积超过第二指定阈值Tv的区域。当Tz为40-50dBz时，所述的雷暴单体为对流单体；当Tz为30_40dBZ时，所述的雷暴单体为对流风暴；当Tz为25-30dBz时，所述的雷暴单体为中尺度对流单体群；当Tz为15-25dBz时，所述的雷暴单体为降雪带。实验仿真数据采用某地基多普勒气象雷达观测数据，评价方式采用预测结果与实际结果对比。与现有技术相比，本发明的有益效果是提高了雷暴単体位置预测的预测精确度，降低了误差。


图I为本发明数据质量控制的原理图。图2a为坐标系转换示意图。图2b为插值过程的示意图。图3为ニ维数据中识别雷暴的例子。图4a为雷暴单体形状示意图。图4b为单体采样示意图。图4c为采样点扩散示意图。图5a为时次140的雷暴单体分布情况。图5b为时次141的雷暴单体分布情况。图6为窄带示意图。图7为140-149时次的真实結果。图8为140-149时次演化力的黑白显示。图9为预测結果。图10为不使用水平集预测的结果，所有的雷暴区域保持形状大小不变平移。图11为使用水平集前后预测精确度的对比。处于上方的线为命中率，处于下方的 线为打偏率。图12为命中率减去打偏率的差值。
具体实施例方式下面结合附图和实施例，对本发明做详细说明，但不应以此限制本发明的保护范围。一种基于水平集的雷暴预测方法，该方法包括如下步骤第一歩、数据质量控制对3D雷达体扫回波数据进行量测数据预处理。图I为本发明数据质量控制的原理图，如图所示，包括噪声滤除、水平反射率纹理计算、垂直反射率纹理计算及非降水回波滤除的预处理过程。第二步、坐标系转化处理，使量测数据所在的坐标系与数据处理所在坐标系一致。对于体扫中的任何一个笛卡尔坐标点，考虑地球曲率的影响，換算其在球坐标系中的位置，设此时其仰角为α，径向距离为r，方位角为Θ，这个位置几乎肯定不在原有球坐标球坐标系的记录中。寻找最靠近α的仰角Ct1;^ α彡Ct2，最靠近Θ的方位角Θ彡Θ 2，径向距离不存在插值(落在哪ー个库是可以直接对应的)。在h仰角层内在θ” θ2方位角之间插值出ー个值，在α2仰角层内在θ” Θ 2方位角之间插值出ー个值，最后再将这两个插值在h、α2仰角层之间插出最终值，插值时两个已知点的权重与距离欲插值点的距离有夫。图2a为坐标系转换示意图，其中，箭头代表雷达的扫描方向，小矩形代表在球坐标系下某一点的高度，大矩形代表考虑地球曲率之后的笛卡尔坐标高度。图2b为插值过程的示意图。第三步、雷暴识别处理，通过对平面上的每ー个点取所在位置所有层中的最大值，将3D数据压缩到2D。雷暴单体是指回波強度超过第一指定阈值Tz，且连续区域面积超过第二指定阈值Tv的区域。显然Tz决定了雷暴单体的类型，如对于对流单体，Tz为40-50dBz，对于对流风暴Tz为30-40dBz，对于中尺度对流单体群，Tz为25-30dBz，对于降雪带，Tz为15_25dBz。定义Tz的大小将直接决定感兴趣区域的大小，40dBz的感兴趣区域将会比30dBz更加紧密和稀少，Tz越�。�感兴趣区域的合并就越多。Tv的作用是防止杂波和噪声，并控制可能出现的目标数量。本实施例中，第一指定阈值Tz为35dBz,第二指定阈值Tv为50km2。首先在横向寻找大于第一指定阈值Tz的点，并将其标上相同序号，直到出现小于第一指定阈值Tz的点，停止标序号。如果接着又发现大于第一指定阈值Tz的点，则序号值加一，继续标相同序号，直到此横向结束，下一横向开始，以此类推。当所有横向捜索完毕，进行纵向捜索，将相互连结的区域统一成一个整体。最后判断各个整体的大�。�面积大于Tv的保留，其他的删除。图3为ニ维数据中识别雷暴的例子，如图所示，雷暴単体4、雷暴単体 5由于面积太�。�将会被删除。对于大单体，随机选取100个点作为采样点。对上ー时刻每ー个雷暴单体采样，如果其含有数据点数小于100个，则所有点都作为采样点，第四步、雷暴跟踪处理，将当前时刻的雷暴目标与上一时刻的雷暴目标相关联。图4a为雷暴単体形状示意图，图4b为单体采样示意图，图4c为采样点扩散图示意图。首先，利用最大相关方法计算当前时刻和上一时刻的雷暴目标图，得出所有雷暴単体的整体移动速度，以此作为粒子随机运动的均值，方差则暂时设为4(如果运动剧烈就要大些)。由此得到扩散后的采样点位置，注意每个采样点来源于上ー时刻哪ー个单体是已知的。
然后统计落入当前时刻所有雷暴単体内的采样点，将得到ー个ニ维矩阵，记录着
当前时刻雷暴単体与上一时刻雷暴単体的关系，横排的列数代表当前时刻的雷暴単体数
量，竖排的列数代表上 一时刻的雷暴单体数量，如 该矩阵表不上ー时
刻的1、3、4号雷暴単体分别与当前时刻的1、3、4号雷暴単体ー脉相承，当前时刻的2号雷暴単体由上ー时刻I号雷暴単体和2号雷暴単体合并而成。某一横排全O代表该雷暴単体新生，某一竖排全O代表该雷暴単体消失。图5所示的情况比较复杂，图5a为时次140的雷暴单体分布情况，图5b为时次141的雷暴単体分布情況。141的雷暴単体的I号区域是140的雷暴単体的I号和2号区域合并得到的，141的2号、4号区域和140的3号区域ー脉相承，141的3号、5号都是140
的4号区域分裂的。矩阵为 5雷暴整体的运动方向和速度可以通过互相
关计算决定，雷暴単体运动方向和速度可以通过中心法确定。中心法运算量虽�。�但由于雷暴単体形状变化多祥，误差较大，故雷暴単体在已经跟踪关联成功的情况下，也可以使用相关法确定运动速度和方向。新生単体因为没有关联，按照背景整体速度移动。第五步、雷暴预测处理，将雷暴単体历史位置的运动轨迹外推得到雷暴単体的未来位置，在该未来位置上，通过水平集演化确定雷暴単体的大小和形状。Mansouri把跟踪问题定义为在已知上一帧目标的情况下，求当前帧目标形状的最大可能，这个可能性受到演化速度项、弾性项和梯度项的影响。演化速度项是当前帧像素点属于目标的概率，弾性项确保曲线光滑，梯度项使得曲线向边界靠拢。本发明与之不同，采用的水平集算法从演变而来。Mansouri处理的问题是跟踪，故其演化速度项是通过对当前帧和上一帧的处理得到的，而本方法是预测形状变化，演化速度项是通过上ー帧单独处理得到的，所以必须加入目标平移速度，而且对于气象目标，梯度项也不再有必要。演化速度项应该反映雷暴的动态演化趋势，所以要从3D信息中得到，定义在每个数据点垂直方向6-12层的数据中寻找大于Tz的点，如果找到3个以上，就认为这一点将要保持或发展为雷暴点，反之则认为此点将要保持或发展为非雷暴点。算法流程如下I)根据雷暴定义得到目标的位置和形状；2)在目标边界周围扩展一圈窄带；3)目标在窄带内演化,演化动力是弹性力和演化速度项。图6为窄带示意图，中间黑色线为目标边界，外圈线为窄带外边界，内圈线为窄带内边界。将预测的结果与下一时刻的真实结果相比，得到两个指标。ー是命中率，即预测结果与真实结果共同拥有的点数与真实结果总点数的比值；另ー是偏离率，即属于预测结果而不是真实结果的点数与真实结果总点数的比值。下列灰度图7是第140-149时次的真实结果，灰度不同代表不同的雷暴区域。
水平集的參数之一，第140-149时次演化力的黑白显示，如图8所示。图9为预测結果。图10为不使用水平集预测的结果，所有的雷暴区域保持形状大小不变平移。图11为使用水平集前后预测精确度的对比，其中，处于上方的线为命中率，处于下方的线为打偏率。图12为命中率减去打偏率的差值，左图的波动区间大小为O. 49，方差为O. 026，右图的波动区间大小为O. 65，方差为O. 0338。从上图可以看到，水平集的确提高了预测精确度，并且降低了误差。注意到有时预测的雷暴単体位置偏离实际情况，这是由于通过最大相关计算的単体移动速度不是百分之 百准确的，如果可以得到更精准的移动速度，那么此处水平及带来的改善还会更可观。精确度方面，由于水平集的演化目前只考虑了回波反射率ー项因素，如果考虑了气旋等因素(有气旋的区域加强雷暴)结果会更好。在本实施例中，由于引入了基于水平集的雷暴単体形状跟踪与预测方法，提高了预测的准确率，降低了误差。
权利要求
1.一种基于水平集的雷暴预测方法，其特征在于，该方法包括如下步骤 第一步、数据质量控制，对3D雷达体扫回波数据进行量测数据预处理； 第二步、坐标系转化处理，使量测数据所在的坐标系与数据处理所在坐标系一致； 第三步、雷暴识别处理，将3D雷达体扫回波数据压缩到2D ； 第四步、雷暴跟踪处理，将当前时刻的雷暴目标与上一时刻的雷暴目标相关联； 第五步、雷暴预测处理，将雷暴单体历史位置的运动轨迹外推得到雷暴单体的未来位置，在该未来位置上，通过水平集演化确定雷暴单体的大小和形状。
2.根据权利要求I所述的基于水平集的雷暴预测方法，其特征在于，所述的雷暴单体是指回波强度超过第一指定阈值Tz，且连续面积超过第二指定阈值Tv的区域。
3.根据权利要求2所述的基于水平集的雷暴预测方法，其特征在于，当Tz为40-50dBz时，所述的雷暴单体为对流单体；当Tz为30-40dBz时，所述的雷暴单体为对流风暴；当Tz为25-30dBz时，所述的雷暴单体为中尺度对流单体群；当Tz为15-25dBz时，所述的雷暴单体为降雪带。
4.根据权利要求1-3任一所述的基于水平集的雷暴预测方法，其特征在于，在所述的雷暴单体的未来位置上，以雷暴单体本身的结构特征作为水平集演化的驱动力之一。
全文摘要
一种基于水平集的雷暴预测方法，包括第一步、数据质量控制，对3D雷达体扫回波数据进行量测数据预处理；第二步、坐标系转化处理，使量测数据所在的坐标系与数据处理所在坐标系一致；第三步、雷暴识别处理，将3D雷达体扫回波数据压缩到2D；第四步、雷暴跟踪处理，将当前时刻的雷暴目标与上一时刻的雷暴目标相关联；第五步、雷暴预测处理，将雷暴单体历史位置的运动轨迹外推得到雷暴单体的未来位置，在该未来位置上，通过水平集演化确定雷暴单体的大小和形状。本发明提高了雷暴单体位置预测的预测精确度，降低了误差。
文档编号G01S13/95GK102662173SQ20121012239
公开日2012年9月12日 申请日期2012年4月24日 优先权日2012年4月24日
发明者宋金泽, 李元祥, 胡琦 申请人:上海交通大学