CN114488350A - 一种基于雷达图信息处理的山地城市短时降雨预报方法 - Google Patents

Abstract

本发明公开了一种基于雷达图信息处理的山地城市短时降雨预报方法，其特征在于，根据云层雷达信息图提取云团信息，并计算云团移动方向和速度，判断需评估时间间隔后云团是否移动至评估区域，并进行降雨判断。本发明能够更好地识别云团并实现追踪，能够避免跟踪丢失和误差累计的问题，对于移速快、分散、复杂的云团仍有较好预测效果，特别适合对于山地城市的短时降雨(尤其是短时强降雨)预报。

Description

一种基于雷达图信息处理的山地城市短时降雨预报方法

技术领域

本发明属于降雨预报领域，具体涉及一种基于雷达图信息处理的山地城市短时降雨预报方法。

背景技术

降雨事件对人类的生活生产起着重要的作用，因此在较短时间内(0-1小时)准确预测出降雨量，以便对恶劣天气进行准确的评估，减少对人民生命财产安全造成的危害，同时也对农业、工业生产起到应急预报的作用。随着气候的变化和城镇化建设加快，各种气象现象频发，大气情况变得越来越复杂，国内外目前缺乏能够实时显示内涝情况的模拟系统，短时、准确的降雨预报方法便是实时内涝模拟系统的关键所在。尤其是在山地城市，由于地势影响原因，风云变幻快速，在强对流天气的山地城市降雨往往存在短时，量大，危害性高等特点，由于山地城市降雨一般持续时间较短，若能准确的对山地城市未来一小时左右的的降雨事件并作出预报并与实时内涝模拟系统所结合，能及时的预报出内涝发生时间与地点，提前做好防备措施，将大大减少所引起的损失。

当前，最常用的降雨量预测方法是雷达回波外推法和时间序列模型方法，由于降水往往具有周期性、突发性、随机性等因素，时间序列模型很可能受此影响，并产生较大误差并影响预测结果。预测需要大量的气象数据，训练的模型需跟着预测区域的变换而重新训练，这种方法不适合山地城市短时降雨预测的情况。

另外，对于降雨预测，还可以采用实时获取多普勒天气雷达的反射率图像，通过跟踪反射率图像的变化趋势，利用雷达反射率强度因子与降雨强度之间的关系，进行降雨事件的预测。其中通过跟踪反射率图像的变化趋势，通常会采用人工神经网络、TITAN、SCIT、光流法等方法。但这些方法仍然存在如下的缺点，包括：TITAN算法只适用于移动缓慢的云，因为存在跟踪丢失的问题；SCIT算法对运动复杂的云团具有较弱的预测性，但可以预测分散的云团；光流法不仅存在上述问题，还容易积累误差；神经网络法在缺乏足够的数据时难以进行训练。

故现有方法均不是特别适用于山地城市的短时降雨预测。

发明内容

针对上述现有技术的不足，本发明所要解决的技术问题是：怎样提供一种能够更好地识别云团并实现追踪，以更加精确地完成山地城市短时降雨预报的基于雷达图信息处理的山地城市短时降雨预报方法，能够避免跟踪丢失和误差累计的问题，对于移速快、分散、复杂的云团仍有较好预测效果，使其特别适合对于山地城市的短时降雨(尤其是短时强降雨)预报，其中短时可指数分钟到数小时范围内的时间。

为了解决上述技术问题，本发明采用了如下的技术方案：

一种基于雷达图信息处理的山地城市短时降雨预报方法，其特征在于，根据云层雷达信息图提取云团信息，并计算云团移动方向和速度，判断需评估时间间隔后云团是否移动至评估区域，再进行降雨判断。

具体地说，本发明包括以下步骤：

步骤a，获取包括预测区域及其周边地区的原图，原图包括未含有云层雷达信息的原始模板图(即和原始雷达图规格一致的但仅含有地形信息的图)，以及含有云层雷达信息的原始雷达图；采用计算机识别获得原始模板图和原始雷达图的RGB色彩值，并计算获得以RGB色彩值为表征的云层雷达信息图；

步骤b，将云层雷达信息图进行二值化处理，形成二值化图像；

步骤c，将二值化图像进行边界提取，划清云团与云团之间的界限，获得云团边界图像；

步骤d，计算获得各云团边界的质心位置；

步骤e，在往前一个时间间隔的时间和当前时间为两个时间点，分别重复步骤a-d，获得两个时间点各自的的云团边界图像和质心位置信息，获得两个时间点各自云团边界图像之间对应关系，根据质心移动变化距离计算当前时间点各云团的质心位移速度并获得质心位移方向；

步骤f，以步骤e获得的各云团质心位移速度和质心位移方向，按照当前时间点各云团的图像位置，计算当前时间点的云层雷达信息图中各云团在需评估时间间隔后所移动到达的位置，判断是否落入预测区域范围；

步骤g，再根据f步骤获得的落入到预测区域范围内的云团信息叠加情况，判断是否降雨以及雨量大小。

这样，本方法直接对预测所在地周边云团信息识别处理，获取其运动速度和方向，进而评估若干时间后周边各云团移动情况，判断是否到达待预测区域，进而实现短时降雨预测，具有步骤简单、快速、高效且可靠性高的特点。特别适合用于山地城市的短时强降雨预测情况。

进一步地，步骤a中，原始模板图以及含有云层雷达信息的原始雷达图均来自气象局发布图片。

这样，保证信息来源可靠性，保证后续计算准确性。

进一步地，步骤a中，计算公式为：

在上述公式(1)中，x，y为雷达图坐标，f₁(x,y)为原始雷达图的RGB色彩强度值，f₂(x,y)为原始模板图RGB色彩强度值，g(x,y)为处理后的强度值。

这样是因为气象局发布的原始雷达图中，云团位置的地形信息已被色彩强烈的云团雷达信息覆盖；故色彩强度相同的区域表明未含有云团信息，故计算为零。而色彩强度不同的地方，因为地形信息已被覆盖(注意是覆盖而不是叠加)，故直接采用原始雷达图中的色彩强度信息，这样即可快速、准确地获得以RGB色彩值为表征的云层雷达信息图。

具体地说，实施时，步骤a中，获取原图的像素，通常整图为512×672像素点比例，其中左边512×512部分为有效图像部分，右边小部分为图像说明介绍。其中有效图像部分的每个像素点包含一个RGB色彩值，可用二维函数f(x,y)来表示，其中x,y是空间坐标，f为坐标点上的强度值，采用MATLAB软件读取图片后为512×672×3阶矩阵，其中雷达基本反射率信息位于原图左边512×512像素之中，利用两张图的二维函数f(x,y)作模板逻辑运算处理，目的在于将时间、仰角、反射强度等包含在其余像素中的信息修剪掉，仅保留云团雷达反射率信息，为后续运算提供依据，保证后续计算过程的精确度和可靠性。

进一步地，步骤a中，所述周边地区范围，是指以预测区域为中心向四周扩展到至少包含一个已达到降雨强度云团所在的区域，或者至少包含两个未达到降雨强度云团所在的区域。

这样是因为短时降雨有可能是已降雨云团移动到位后导致降雨，或者多个未降雨云团碰撞后导致降雨，故这样可以更好地兼顾了多种实际情况，更好地保证了预测精度。

进一步地，在所述步骤b中二值化处理过程为：

在上述公式(2)中，A为雷达反射率强度大于5dBZ的色彩强度值集合。

这样，在步骤b中将所述模板逻辑运算处理结果图(即云层雷达信息图)进行正向二值化处理，得到二值化图像，二值化处理的目的在于图中还包含有地名、地貌等噪音信息，该部分信息与研究对象的云层重合，易干扰算法，需将其分离，仅保留反射率信息，并且处理过程中去除掉反射率强度0-5dBZ的代表颜色，保留其余颜色，原因在于雷达反射率强度小于5dBZ对降雨几乎没有贡献。这样可以更好地区分云团边界，使其有利于后续计算。

进一步地，步骤c中将所述二值化图像使用Bwperim-imfill算子提取细节更加完美的边界图像G(x,y)，bwperim-imfill组合算子步骤为先使用bwperim算子对图像进行边界提取，再使用imfill算子处理提取的边界，对处理结果再次使用bwperim算子计算。

这样，bwperim-imfill组合算子在Bwperim算子基础上进行填充腐蚀操作，该方法在多次腐蚀过程中将距离较近的云层合并到一起，极大程度的降低原图噪音干扰带来的影响，恢复了原有云团样貌，划清了云团与云团之间的界限，在Bwperim-imfill算子提取结果G(x,y)之上实现风暴云团的追踪，G(x,y)经过imfill算子处理之后，所有云团边界必定处于闭合状态，边界部分在图上体现为白色，g(x,y)值为1，其余为黑色，g(x,y)值为0。这样，可以更好地实现云团边界区分。其中使用到了Bwperim-imfill组合算子，Bwperim是形态学边界提取算子，Imfill是形态学填充算子，分别使用膨胀运算和腐蚀运算提取更加完美的边界图像G(x,y)。但几种算子各自自身均属于现有技术，本申请的贡献在于其组合应用，故各自的具体计算过程不在此详述。

进一步地，步骤d中，计算获得各云团边界的质心位置具体包括以下步骤：

1)利用追踪算法来获取距离预测区域最近云团边界满足下式的点，公式为：

在上述公式(3)中，d₈为棋盘距离，D为欧式距离，G为经算法提取后的边界图像，p为预测区域中心点坐标，Q为预测区域最近云团边界某点，X₀，Y₀为预测区域最近云团边界某点坐标；

2)利用以下(自定义)算法获取云团质心坐标，公式为：

在上述公式(4)中，

为云团质心，n为某个云团边界点个数，X_k，Y_k为某个云团所有边界点的坐标。

上述步骤d的具体计算步骤原理为：根据云团边界图像实现风暴云团的追踪，追踪得到离预测区域最近的边界点，得到最近的边界点后，利用自定义算法求得云团的质心。由步骤c得到的图形一定是云团与云团之间划清了界限，每个云团都是独立闭合的，且云团边界值G(x,y)为1，其余值为0。

这样，在步骤1)中，以预测区域中心点p为棋盘的中心，按照图像中1个像素点为1个棋盘距离，以棋盘距离由小到大开始获取G(x,y)＝1的点，若找到某棋盘距离发现有G(x,y)＝1的点，若这样的点不止一个点，则比较这些点与预测区域中心点的欧式距离，取最小值为预测区域最近云团边界点Q(X₀,Y₀)。

在步骤2)中寻找离预测区域最近云团全部边界坐标，作如下步骤，上述过程中求得了预测区域最近云团边界某点Q(X₀,Y₀)，以这个点为棋盘的中心，棋盘距离等于1的点有8个，这8个点至少有两个点满足G(x,y)＝1，设第k个点为Q_k(X_k,Y_k)(k＝1,2,3,…,8)，以Q点8领域中随机一点为初始起点，按顺时针方向对8个点逐一判断，若满足：G(X_k,Y_k)＝1，则将该点坐标储存至集合C中，并又以该点为棋盘的中心，又重复上述操作，直至同时满足X_k＝X₀,Y_k＝Y₀，这时整个云团G(x,y)＝1的点都遍历了一遍，即找到了云团所有边界点并储存在C中，利用这些点按照质心公式计算出云团质心，同理也可得到所有云团的质心。

故采用上述自定义算法，能够以计算机能识别和执行的方式，方便快捷地完成对云团质心的准确计算，运算快捷且可靠性高。

进一步地，所述步骤e中，所述时间间隔为气象局发布气象信息最短时间间隔。通常为6分钟。

这样，可以更好地提高预测精度，实现短时降雨预测。

进一步地，所述步骤e中，根据前后时间点云团数量差异，按照以下方式计算各云团的质心位移速度：

1)云团前后时刻数量相等时，云团速度：

2)云团前后时刻数量不相等时，云团速度：

上述公式(5)和公式(6)中，T为t时刻与下一时刻图片的时间间隔，一般为6min；r为t时刻云团数量，r′为t+T时刻云团数量；

分别为t时刻图片的第m个云团质心的横坐标与纵坐标，

为t+T时刻图片的第m个云团质心的横坐标与纵坐标，m＝1，2，3，…，r；

这样，当前后时间点云团数量相同时，分别根据各云团之间对应关系，根据云团位置变化计算出各云团的质心移动速度，作为当前时间各云团即将要位移的速度和方向。而当云团数量变换时，表明有云团散开或聚合，故此时在前后各自时间点上，计算出各自所有云团的共同质心，再根据共同质心的位置变化计算出质心移动速度和方向，作为当前时间各云团即将要移动的速度和方向。这样就保证了各种情况下均可以比较准确地获得当前云团所要移动的速度和方向，以方便后续计算。

进一步地，云团前后时刻数量相等时，按照前后时刻的图片中，以同一方向对各云团顺序排序并按照排序建立云团对应关系。

这样，是因为本方法用于山地城市短时降雨预报，前后时间间隔大小是按照气象局发布消息最短时间间隔计算，通常为6分钟，时间非常短暂，故各云团移动变化情况不大。所以采用此方法可以方便快捷且相对可靠地建立其前后云团的对应关系，以方便后续计算。

进一步地，所述步骤e中，各云团移动方向表征为：

在上述公式(7)中，α为云团运动方向与直角坐标系x轴夹角；dx，dy分别表示t+T时刻至t时刻质心在x，y轴上的增量。

这样，以角度大小的方式表征各云团的方向，方便后续以坐标变化的方式判断云团是否到达预测区域。

进一步地，步骤f中，使用坐标变换的方式判断评估时间间隔后，各云团是否落入预测区域范围，坐标变换的方式为以坐标系原点为中心，按照云团移动方向夹角，逆时针旋转进行变换，坐标变换公式为：

在上述公式(8)中，(x′_i，y′_i)为旋转后各点坐标，(x′₀，y′₀)为预测区域中心点坐标，第m个云团运动方向与直角坐标系X轴夹角为α_m。

这样步骤f中，求出云团移动方向后采用坐标系变换判断是否能到达预测区域，并判断降雨事件是否发生，这样做的目的在于简化计算，避免云团自身外形不规则对计算带来的困扰，如果把坐标系旋转到云团移动方向α_m，只需要判断预测区域旋转后的纵坐标y′₀是否在y′_i，min，y′_i，max之间，若不在，则取欧式距离由小到大对各云团边界依次判断，直至满足该条件，若全不满足，则表明各云团在评估时间内均不会落入到预测区域，故该评估时间时不出现降雨。其中y′_i，min，y′_i，max分别为坐标变换后云团在y′轴上的最小最大值，能直观的感受出云团是否能到达预测区域。其中降雨时间为预测区域中心点与距离预测区域最近点M的欧式距离除以所述步骤e的移动速度。

进一步地，步骤g中，采用未拟合的经典Z-R关系式模型进行降雨强度预测。

这样，根据预测区域到达云团的叠加的反射率强度dBZ来计算雷达反射率因子Z，再通过典型的Z-R关系式来预测降雨强度，为成熟的现有计算方式，通常预测最大降雨强度误差为10％，误差较低。

故相对于现有技术，本发明的技术特点在于：针对传统降雨预报方法误差大，需要大量研究区域气象资料等问题，提出了新的预报模型，利用Bwperim-imfill算子对云团边界进行检测，能有效识别云团边界，结合TITAN等传统方法优势，提出新算法实现云团的追踪，避免跟踪丢失、累计误差的问题，对于移速快、分散、复杂的云团仍有较好预测效果。

综上所述，本发明能够更好地识别云团并实现追踪，能够避免跟踪丢失和误差累计的问题，对于移速快、分散、复杂的云团仍有较好预测效果，特别适合对于山地城市的短时降雨(尤其是短时强降雨)预报。

附图说明

图1为具体实施方式中步骤a中获得的原始雷达图，为满足专利图片格式要求，图中色彩转换为了灰度。

图2为步骤a中获得的云层雷达信息图，为满足专利图片格式要求，图中色彩转换为了灰度。

图3为步骤b获得的二值化图像。

图4为步骤c中先使用bwperim算子对图像进行边界提取后的结果示意图。

图5为步骤c中使用bwperim-imfill组合算子对图像处理后的结果示意图。

图6为步骤f中坐标系旋转变换的示意图。

具体实施方式

下面结合具体实施方式对本发明作进一步的详细说明。

具体实施方式：一种基于雷达图信息处理的山地城市短时降雨预报方法，其特点在于，根据云层雷达信息图提取云团信息，并计算云团移动方向和速度，判断需评估时间间隔后云团是否移动至评估区域，再进行降雨判断。

本发明包括以下步骤：

步骤a，获取包括预测区域及其周边地区的原图，原图包括未含有云层雷达信息的原始模板图(即和原始雷达图规格一致的但仅含有地形信息的图)，以及含有云层雷达信息的原始雷达图(参见图1)；采用计算机识别获得原始模板图和原始雷达图的RGB色彩值，并计算获得以RGB色彩值为表征的云层雷达信息图(参见图2)；

步骤b，将云层雷达信息图进行二值化处理，形成二值化图像(参见图3)；

步骤c，将二值化图像进行边界提取，划清云团与云团之间的界限，获得云团边界图像；

步骤d，计算获得各云团边界的质心位置；

步骤e，在往前一个时间间隔的时间和当前时间为两个时间点，分别重复步骤a-d，获得两个时间点各自的的云团边界图像和质心位置信息，获得两个时间点各自云团边界图像之间对应关系，根据质心移动变化距离计算当前时间点各云团的质心位移速度并获得质心位移方向；

步骤f，以步骤e获得的各云团质心位移速度和质心位移方向，按照当前时间点各云团的图像位置，计算当前时间点的云层雷达信息图中各云团在需评估时间间隔后所移动到达的位置，判断是否落入预测区域范围；

步骤g，再根据f步骤获得的落入到预测区域范围内的云团信息叠加情况，判断是否降雨以及雨量大小。

这样，本方法直接对预测所在地周边云团信息识别处理，获取其运动速度和方向，进而评估若干时间后周边各云团移动情况，判断是否到达待预测区域，进而实现短时降雨预测，具有步骤简单、快速、高效且可靠性高的特点。特别适合用于山地城市的短时强降雨预测情况。

本实施方式中，步骤a中，原始模板图以及含有云层雷达信息的原始雷达图均来自气象局发布图片。

这样，保证信息来源可靠性，保证后续计算准确性。

本实施方式中，步骤a中，计算公式为：

在上述公式(1)中，x，y为雷达图坐标，f₁(x,y)为原始雷达图的RGB色彩强度值，f₂(x,y)为原始模板图RGB色彩强度值，g(x,y)为处理后的强度值。

这样是因为气象局发布的原始雷达图中，云团位置的地形信息已被色彩强烈的云团雷达信息覆盖；故色彩强度相同的区域表明未含有云团信息，故计算为零。而色彩强度不同的地方，因为地形信息已被覆盖(注意是覆盖而不是叠加)，故直接采用原始雷达图中的色彩强度信息，这样即可快速、准确地获得以RGB色彩值为表征的云层雷达信息图。

具体地说，实施时，步骤a中，获取原图的像素，通常整图为512×672像素点比例，其中左边512×512部分为有效图像部分，右边小部分为图像说明介绍。其中有效图像部分的每个像素点包含一个RGB色彩值，可用二维函数f(x,y)来表示，其中x,y是空间坐标，f为坐标点上的强度值，采用MATLAB软件读取图片后为512×672×3阶矩阵，其中雷达基本反射率信息位于原图左边512×512像素之中，利用两张图的二维函数f(x,y)作模板逻辑运算处理，目的在于将时间、仰角、反射强度等包含在其余像素中的信息修剪掉，仅保留云团雷达反射率信息，为后续运算提供依据，保证后续计算过程的精确度和可靠性。

本实施方式中，步骤a中，所述周边地区范围，是指以预测区域为中心向四周扩展到至少包含一个已达到降雨强度云团所在的区域，或者至少包含两个未达到降雨强度云团所在的区域。

这样是因为短时降雨有可能是已降雨云团移动到位后导致降雨，或者多个未降雨云团碰撞后导致降雨，故这样可以更好地兼顾了多种实际情况，更好地保证了预测精度。

本实施方式中，在所述步骤b中二值化处理过程为：

在上述公式(2)中，A为雷达反射率强度大于5dBZ的色彩强度值集合。

这样，在步骤b中将所述模板逻辑运算处理结果图(即云层雷达信息图)进行正向二值化处理，得到二值化图像，二值化处理的目的在于图中还包含有地名、地貌等噪音信息，该部分信息与研究对象的云层重合，易干扰算法，需将其分离，仅保留反射率信息，并且处理过程中去除掉反射率强度0-5dBZ的代表颜色，保留其余颜色，原因在于雷达反射率强度小于5dBZ对降雨几乎没有贡献。这样可以更好地区分云团边界，使其有利于后续计算。

本实施方式中，步骤c中将所述二值化图像使用Bwperim-imfill算子提取细节更加完美的边界图像G(x,y)，bwperim-imfill组合算子步骤为先使用bwperim算子对图像进行边界提取(结果参见图4)，再使用imfill算子处理提取的边界，对处理结果再次使用bwperim算子计算(结构参见图5)。

这样，bwperim-imfill组合算子在Bwperim算子基础上进行填充腐蚀操作，该方法在多次腐蚀过程中将距离较近的云层合并到一起，极大程度的降低原图噪音干扰带来的影响，恢复了原有云团样貌，划清了云团与云团之间的界限，在Bwperim-imfill算子提取结果G(x,y)之上实现风暴云团的追踪，G(x,y)经过imfill算子处理之后，所有云团边界必定处于闭合状态，边界部分在图上体现为白色，g(x,y)值为1，其余为黑色，g(x,y)值为0。这样，可以更好地实现云团边界区分。其中使用到了Bwperim-imfill组合算子，Bwperim是形态学边界提取算子，Imfill是形态学填充算子，分别使用膨胀运算和腐蚀运算提取更加完美的边界图像G(x,y)。但几种算子各自自身均属于现有技术，本申请的贡献在于其组合应用，故各自的具体计算过程不在此详述。

本实施方式中，步骤d中，计算获得各云团边界的质心位置具体包括以下步骤：

1)利用追踪算法来获取距离预测区域最近云团边界满足下式的点，公式为：

在上述公式(3)中，d₈为棋盘距离，D为欧式距离，G为经算法提取后的边界图像，p为预测区域中心点坐标，Q为预测区域最近云团边界某点，X₀，Y₀为预测区域最近云团边界某点坐标；

2)利用以下(自定义)算法获取云团质心坐标，公式为：

在上述公式(4)中，

为云团质心，n为某个云团边界点个数，X_k，Y_k为某个云团所有边界点的坐标。

上述步骤d的具体计算步骤原理为：根据云团边界图像实现风暴云团的追踪，追踪得到离预测区域最近的边界点，得到最近的边界点后，利用自定义算法求得云团的质心。由步骤c得到的图形一定是云团与云团之间划清了界限，每个云团都是独立闭合的，且云团边界值G(x,y)为1，其余值为0。

这样，在步骤1)中，以预测区域中心点p为棋盘的中心，按照图像中1个像素点为1个棋盘距离，以棋盘距离由小到大开始获取G(x,y)＝1的点，若找到某棋盘距离发现有G(x,y)＝1的点，若这样的点不止一个点，则比较这些点与预测区域中心点的欧式距离，取最小值为预测区域最近云团边界点Q(X₀,Y₀)。

在步骤2)中寻找离预测区域最近云团全部边界坐标，作如下步骤，上述过程中求得了预测区域最近云团边界某点Q(X₀,Y₀)，以这个点为棋盘的中心，棋盘距离等于1的点有8个，这8个点至少有两个点满足G(x,y)＝1，设第k个点为Q_k(X_k,Y_k)(k＝1,2,3,…,8)，以Q点8领域中随机一点为初始起点，按顺时针方向对8个点逐一判断，若满足：G(X_k,Y_k)＝1，则将该点坐标储存至集合C中，并又以该点为棋盘的中心，又重复上述操作，直至同时满足X_k＝X₀,Y_k＝Y₀，这时整个云团G(x,y)＝1的点都遍历了一遍，即找到了云团所有边界点并储存在C中，利用这些点按照质心公式计算出云团质心，同理也可得到所有云团的质心。

故采用上述自定义算法，能够以计算机能识别和执行的方式，方便快捷地完成对云团质心的准确计算，运算快捷且可靠性高。

本实施方式中，所述步骤e中，所述时间间隔为气象局发布气象信息最短时间间隔。通常为6分钟。

这样，可以更好地提高预测精度，实现短时降雨预测。

本实施方式中，所述步骤e中，根据前后时间点云团数量差异，按照以下方式计算各云团的质心位移速度：

1)云团前后时刻数量相等时，云团速度：

2)云团前后时刻数量不相等时，云团速度：

上述公式(5)和公式(6)中，T为t时刻与下一时刻图片的时间间隔，一般为6min；r为t时刻云团数量，r′为t+T时刻云团数量；

分别为t时刻图片的第m个云团质心的横坐标与纵坐标，

为t+T时刻图片的第m个云团质心的横坐标与纵坐标，m＝1，2，3，…，r。

这样，当前后时间点云团数量相同时，分别根据各云团之间对应关系，根据云团位置变化计算出各云团的质心移动速度，作为当前时间各云团即将要位移的速度和方向。而当云团数量变换时，表明有云团散开或聚合，故此时在前后各自时间点上，计算出各自所有云团的共同质心，再根据共同质心的位置变化计算出质心移动速度和方向，作为当前时间各云团即将要移动的速度和方向。这样就保证了各种情况下均可以比较准确地获得当前云团所要移动的速度和方向，以方便后续计算。

本实施方式中，云团前后时刻数量相等时，按照前后时刻的图片中，以同一方向对各云团顺序排序并按照排序建立云团对应关系。

这样，是因为本方法用于山地城市短时降雨预报，前后时间间隔大小是按照气象局发布消息最短时间间隔计算，通常为6分钟，时间非常短暂，故各云团移动变化情况不大。所以采用此方法可以方便快捷且相对可靠地建立其前后云团的对应关系，以方便后续计算。

本实施方式中，所述步骤e中，各云团移动方向表征为：

在上述公式(7)中，α为云团运动方向与直角坐标系x轴夹角；dx，dy分别表示t+T时刻至t时刻质心在x，y轴上的增量。

这样，以角度大小的方式表征各云团的方向，方便后续以坐标变化的方式判断云团是否到达预测区域。

本实施方式中，步骤f中，使用坐标变换的方式判断评估时间间隔后，各云团是否落入预测区域范围，坐标变换的方式为以坐标系原点为中心，按照云团移动方向夹角，逆时针旋转进行变换(参见图6)，坐标变换公式为：

在上述公式(8)中，(x′_i，y′_i)为旋转后各点坐标，(x′₀，y′₀)为预测区域中心点坐标，第m个云团运动方向与直角坐标系X轴夹角为α_m。

这样步骤f中，求出云团移动方向后采用坐标系变换判断是否能到达预测区域，并判断降雨事件是否发生，这样做的目的在于简化计算，避免云团自身外形不规则对计算带来的困扰，如果把坐标系旋转到云团移动方向α_m，只需要判断预测区域旋转后的纵坐标y′₀是否在y′_i，min，y′_i，max之间，若不在，则取欧式距离由小到大对各云团边界依次判断，直至满足该条件，若全不满足，则表明各云团在评估时间内均不会落入到预测区域，故该评估时间时不出现降雨。其中y′_i，min，y′_i，max分别为坐标变换后云团在y′轴上的最小最大值，能直观的感受出云团是否能到达预测区域。其中降雨时间为预测区域中心点与距离预测区域最近点M的欧式距离除以所述步骤e的移动速度。

本实施方式中，步骤g中，采用未拟合的经典Z-R关系式模型进行降雨强度预测。

这样，根据预测区域到达云团的叠加的反射率强度dBZ来计算雷达反射率因子Z，再通过典型的Z-R关系式来预测降雨强度，为成熟的现有计算方式，通常预测最大降雨强度误差为10％，误差较低。

故相对于现有技术，本发明的技术特点在于：针对传统降雨预报方法误差大，需要大量研究区域气象资料等问题，提出了新的预报模型，利用Bwperim-imfill算子对云团边界进行检测，能有效识别云团边界，结合TITAN等传统方法优势，提出新算法实现云团的追踪，避免跟踪丢失、累计误差的问题，对于移速快、分散、复杂的云团仍有较好预测效果。

Claims (10)

1.一种基于雷达图信息处理的山地城市短时降雨预报方法，其特征在于，根据云层雷达信息图提取云团信息，并计算云团移动方向和速度，判断需评估时间间隔后云团是否移动至评估区域，再进行降雨判断。

2.如权利要求1所述的基于雷达图信息处理的山地城市短时降雨预报方法，其特征在于，包括以下步骤：

步骤a，获取包括预测区域及其周边地区的原图，原图包括未含有云层雷达信息的原始模板图，以及含有云层雷达信息的原始雷达图；采用计算机识别获得原始模板图和原始雷达图的RGB色彩值，并计算获得以RGB色彩值为表征的云层雷达信息图；

步骤b，将云层雷达信息图进行二值化处理，形成二值化图像；

步骤c，将二值化图像进行边界提取，划清云团与云团之间的界限，获得云团边界图像；

步骤d，计算获得各云团边界的质心位置；

步骤e，在往前一个时间间隔的时间和当前时间为两个时间点，分别重复步骤a-d，获得两个时间点各自的的云团边界图像和质心位置信息，获得两个时间点各自云团边界图像之间对应关系，根据质心移动变化距离计算当前时间点各云团的质心位移速度并获得质心位移方向；

步骤f，以步骤e获得的各云团质心位移速度和质心位移方向，按照当前时间点各云团的图像位置，计算当前时间点的云层雷达信息图中各云团在需评估时间间隔后所移动到达的位置，判断是否落入预测区域范围；

步骤g，再根据f步骤获得的落入到预测区域范围内的云团信息叠加情况，判断是否降雨以及雨量大小。

3.如权利要求2所述的基于雷达图信息处理的山地城市短时降雨预报方法，其特征在于，步骤a中，原始模板图以及含有云层雷达信息的原始雷达图均来自气象局发布图片。

4.如权利要求2所述的基于雷达图信息处理的山地城市短时降雨预报方法，其特征在于，步骤a中，计算公式为：

在上述公式(1)中，x，y为雷达图坐标，f₁(x,y)为原始雷达图的RGB色彩强度值，f₂(x,y)为原始模板图RGB色彩强度值，g(x,y)为处理后的强度值。

5.如权利要求2所述的基于雷达图信息处理的山地城市短时降雨预报方法，其特征在于，步骤a中，所述周边地区范围，是指以预测区域为中心向四周扩展到至少包含一个已达到降雨强度云团所在的区域，或者至少包含两个未达到降雨强度云团所在的区域。

6.如权利要求2所述的基于雷达图信息处理的山地城市短时降雨预报方法，其特征在于，在所述步骤b中二值化处理过程为：

在上述公式(2)中，A为雷达反射率强度大于5dBZ的色彩强度值集合。

7.如权利要求2所述的基于雷达图信息处理的山地城市短时降雨预报方法，其特征在于步骤c中将所述二值化图像使用Bwperim-imfill算子提取细节更加完美的边界图像G(x,y)，bwperim-imfill组合算子步骤为先使用bwperim算子对图像进行边界提取，再使用imfill算子处理提取的边界，对处理结果再次使用bwperim算子计算。。

8.如权利要求2所述的基于雷达图信息处理的山地城市短时降雨预报方法，其特征在于：步骤d中，计算获得各云团边界的质心位置具体包括以下步骤：

1)利用追踪算法来获取距离预测区域最近云团边界满足下式的点，公式为：

在上述公式(3)中，d₈为棋盘距离，D为欧式距离，G为经算法提取后的边界图像，p为预测区域中心点坐标，Q为预测区域最近云团边界某点，X₀，Y₀为预测区域最近云团边界某点坐标；

2)利用以下(自定义)算法获取云团质心坐标，公式为：

在上述公式(4)中，

为云团质心，n为某个云团边界点个数，X_k，Y_k为某个云团所有边界点的坐标。

9.如权利要求2所述的基于雷达图信息处理的山地城市短时降雨预报方法，其特征在于，所述步骤e中，所述时间间隔为气象局发布气象信息最短时间间隔；

所述步骤e中，根据前后时间点云团数量差异，按照以下方式计算各云团的质心位移速度：

1)云团前后时刻数量相等时，云团速度：

2)云团前后时刻数量不相等时，云团速度：

上述公式(5)和公式(6)中，T为t时刻与下一时刻图片的时间间隔，一般为6min；r为t时刻云团数量，r′为t+T时刻云团数量；

分别为t时刻图片的第m个云团质心的横坐标与纵坐标，

为t+T时刻图片的第m个云团质心的横坐标与纵坐标，m＝1，2，3，…，r；

云团前后时刻数量相等时，按照前后时刻的图片中，以同一方向对各云团顺序排序并按照排序建立云团对应关系。

10.如权利要求9所述的基于雷达图信息处理的山地城市短时降雨预报方法，其特征在于，所述步骤e中，各云团移动方向表征为：

在上述公式(7)中，α为云团运动方向与直角坐标系x轴夹角；dx，dy分别表示t+T时刻至t时刻质心在x，y轴上的增量；

步骤f中，使用坐标变换的方式判断评估时间间隔后，各云团是否落入预测区域范围，坐标变换的方式为以坐标系原点为中心，按照云团移动方向夹角，逆时针旋转进行变换，坐标变换公式为：

在上述公式(8)中，(x′_i，y′_i)为旋转后各点坐标，(x′₀，y′₀)为预测区域中心点坐标，第m个云团运动方向与直角坐标系X轴夹角为α_m；

步骤g中，采用未拟合的经典Z-R关系式模型进行降雨强度预测。

Images