基于气象雷达图信息处理的短时降雨预报方法.pdf

水 文JOURNAL OF CHINA HYDROLOGY第43卷第3期2023年6月Vol.43 No.3Jun.,2023DOI：10.19797/ki.1000-0852.20220041收稿日期：2021-01-24网络首发日期：2023-02-17网络首发地址：https:/ 河海学院，重庆400074；2.重庆交通大学 环境水利工程重庆市工程实验室，重庆400074）摘要：气象雷达图显示的风暴云团具有极强的不规则性，导致传统降雨预报方法普遍存在累积误差大、云团追踪易丢失、气象资料依赖性强等问题。采用正向二值化方法对原始雷达图进行预处理，利用Bwperim-Imfill算子可实现雷达图中分散且不规则云团边缘的有效检测，在此基础上提出的风暴云团追踪模型可识别云团移动方向，估测到达研究区域的时间（即降雨时间）。结果表明，在气象资料严重匮乏条件下，利用风暴云团追踪模型并根据未拟合的经典Z-R关系式可实现单站雷达覆盖区域内的降雨预报。最大降雨强度预报误差为10%，降雨发生时间预报误差为3 min，降雨结束时间预报误差为16 min，纳什效率系数（NSE）高达0.71，均方根误差（RMSE）为12.05mm/h，满足模型应用精度要求。关键词：雷达图；风暴云团；边缘检测；云团追踪；Z-R关系式；降雨强度中图分类号：TV11；P456；P338文献标识码：A文章编号：1000-0852（2023）03-0076-060引言目前，通常采用雨量站和气象雷达两种方法对降雨过程进行预报1。雨量站提供的气象资料虽可相对准确地测量近地表点位的降雨量，但其无法捕捉降雨的空间变异性，也难以预测长时间序列的降雨过程，进而影响径流模型的模拟精度2。相反，气象雷达不仅具有流域尺度降雨实时探测功能，还可捕捉降雨的空间变异性，从而实现短时降雨预报，甚至可提前预报降雨过程3。通常使用雷达回波外推法、时间序列模型法和神经网络模型法等方法实现降雨量的预测4-6。其中，时间序列模型法对历史气象资料依赖性高，且由于降雨具有周期性、突变性和随机性等特点7，导致预测结果易产生较大误差；而神经网络模型法的训练同样依赖于气象资料样本数，难以适用于气象资料匮乏情况，且模型具有区域针对性。雷达回波外推法主要通过追踪反射率的图像变化趋势，利用雷达反射率强度因子与降雨强度间的关系实现降雨预报8。该方法主要分为风暴云团边界识别和云团追踪两个步骤，且云团边界识别是追踪云团的关键步骤9。基于图像信息的边缘检测算法可实现云团边界识别目标，但由于风暴云团具有极强的不规则性，而传统算法识别效果有限，可能造成云团边界识别不闭合，易导致云团追踪时发生目标丢失或错误追踪等现象。因此，选择正确的边缘检测算法尤为关键。目前常用TITAN、SCIT、光流法等算法进行云团追踪，其中TITAN算法存在跟踪丢失问题，仅适用移速缓慢的云团；SCIT算法则适用于分散云团的跟踪，难以预测复杂云团的运动情况；而光流法不仅存在跟踪丢失，难以跟踪复杂云团外，还存在易产生累积误差问题10-11。基于此，本文考察了模版逻辑运算、大津法（OTSU）和正向二值化法对于雷达图的预处理效果；探讨了Bwperim-Imfill算子的云团边界检测效果；结合TITAN等传统方法，提出了云团追踪模型；最后在实现云团追踪的基础上结合经典Z-R关系式预测降雨强度、降雨发生和结束时刻，并采用重庆市秀山县的一程麒铭等：基于气象雷达图信息处理的短时降雨预报方法第3期场实测降雨数据进行了验证，以期提出一种适用于气象资料严重匮乏条件下利用气象雷达图信息实现单站雷达覆盖区域内降雨预报的方法。1研究方法1.1数据来源与预处理所用雷达图来源于中央气象台（http：/ dBZ对降雨几乎没有贡献12，故去除反射率强度低于5 dBZ的代表颜色，并将其余14种色彩RGB值提取至集合A中，二值化后图像g()x,y可由式（1）表示。g()x,y=1，f1()x,yA0，f1()x,yA（1）采用模版逻辑运算、OTSU法和正向二值化法分别对原始图进行处理，可得到图1所示的分离效果。结果表明，采用最大类间方差法确定OTSU法分离阈值t为0.62，但图像处理效果不佳，难以去除地名等噪音信息，且以0.01为阈值间隔，遍历0-1阈值取值区域，分离结果无法达到预期目标；模版逻辑运算与正向二值化的处理效果基本一致，但前者具有单一性，模板受研究区域限制，不具备普适性。而正向二值化法则适用于所有区域的单站雷达图，且处理后图片矩阵降为512512，大幅降低了原始图片数据量，可提高后续云团追踪算法的运算速度。1.2云团边缘检测方法在尝试利用Sobel、Prewitt、Canny以及传统形态学4种边缘检测手段13-17对云团边界进行处理，均不能获得较好的边缘处理效果后，本文采用Bwperim-Imfill组合算子进行检测。g()x,y为二值化处理后图像，将其进行腐蚀与膨胀运算（具体步骤详见文献18）。形态学边界提取采用Bwperim算子19，形态学填充采用Imfill算子20。在进行检测运算时依次采用Bwperim算子、Imfill算子、Bwperim算子，重复该步骤，直至满足Imfill算子收敛条件，并将处理后最终图像用G()x,y表示。1.3降雨预测1.3.1云团追踪方法本文尝试在Bwperim-Imfill算子提取结果G()x,y基础上实现风暴云团的追踪。具体为G()x,y经Imfill算子处理后，所有云团边界必定处于闭合状态，边界部分在图上体现为白色。设研究区域中心点p坐标为()x0,y0，雷达图获取时间为t，中心点与图像G()x,y上任意一点q()xi,yj的棋盘距离为d8，欧拉距离为D，追踪步骤如下：若二值图g()x0,y0值为0，则表明t时刻研究区域未降雨。利用棋盘距离公式寻找研究区域最近云团边界某点Q()X0,Y0，该点满足以下条件：Mind8()p,Q+D()p,Q（2）s.t.G()X0,Y0=1X0,Y0=1,2,3,512（3）为寻找距离研究区域最近云团全部边界坐标，定义算法如下：满足Q()X0,Y0点棋盘距离d8()q,Q=1的点共 8个，设第k个点为Qk()Xk,Yk()k=1,2,3,8，以Q点8领域中随机一点为初始起点，按顺时针方向对8个点注：（a）原图；（b）模板逻辑运算；（c）OTSU二值化运算；（d）正向二值化运算。图1不同方法下风暴云图图像预处理效果Fig.1 Image pre-processing results of storm cloud images with different methods（a）（b）（c）（d）基本反射率雷达站名：铜仁观测日期：2021-05-29观测时间：00:26:19BJT.98JT.9BJ.99BJT数据范围：125 km显示仰角：1.570dBZ6560555045403530252015105077第43卷水 文逐一判断，若满足：G()Xk,Yk=1（4）则将该点坐标储存至集合C中，并以Qk点代替Q点，重复上述步骤直至同时满足Xk=X0,Yk=Y0，集合C中坐标不重复点数量至少4个。利用上述方法可找到图中所有云团C1,C2,C3,Cr，则第m个云团质心计算方法为：x=1nk=1nXk（5）y=1nk=1nYk（6）式中：n为集合C中坐标点个数，n=4,5,6。1.3.2降雨时间预测方法假设在一个时间间隔（6 min）内所有云团的运动速度与方向相等，由于云团可能出现分离、合并、消散的情况，导致难以追踪。若t时刻的云团数量与下一时刻不相等，则说明云团已发生分离、合并或消散，该情况下的云团移动速度V可用式（7）表示。V=1T1rm=1r-xt,m-1rm=1r-xt+T,m2+1rm=1r-yt,m-1rm=1r-yt+T,m2（7）若t时刻的云团数量与下一时刻相等，云团移动速度V可用式（8）表示。V=1rm=1r1T()-xt,m-xt+T,m2+()-yt,m-yt+T,m2（8）式中：T为t时刻与下一时刻图片的时间间隔，一般为6 min；r为t时刻云团数量，r为下一时刻云团数量；-xt,m，-yt,m为 t 时 刻 图 片 的 第 m 个 云 团 质 心，m=1,2,3,r。云团移动方向可用下式表示：=arctan-yt+T,m-yt,m-xt+T,m-xt,mdx0arctan-yt+T,m-yt,m-xt+T,m-xt,m+dx0-2dx=0,dy0（9）式中：为云团运动方向与直角坐标系x轴夹角；dx,dy分别表示t+T时刻至t时刻质心在x,y轴上的增量。在得到m个云团的移动速度Vm与方向m后，为判断是否能到达研究区域，可将坐标系以原点为中心，逆时针旋转，坐标系变换公式为：xiyi=cosmsinm-sinmcosmxiyi（10）旋转后各点坐标为()xi,yi，研究区域中心点坐标为()x0,y0，若云团能到达研究区域，则降雨能发生，该云团在yi轴方向上的投影值如图2。()x0,yi,max()x0,y0()x0,yi,minxyyxOamam质心图2坐标系旋转过程Fig.2Rotation process of coordinate system且应满足：yi,miny0yi,max（11）式中：yi,min,yi,max分别为坐标变换后云团在y轴上的最小最大值。距离研究区域最近点M的横坐标为xi,max，降雨时间Train可用式（12）表示。Train=1VD()p,M（12）1.3.3降雨强度预测方法雷达反射率因子Z与降雨强度R之间存在的联系通常采用Z-R关系式进行表示21-22，由于缺少足够气象资料进行拟合，本文采用典型Z-R关系式（即Marshall-Palmer公式）23，见式（13）。Z=200R1.6（13）2结果与讨论2.1边缘检测结果在MATLAB 中编程实现Prewitt、Sobel和Canny经典算子，以及传统形态学方法对图1预处理后的图形进行边界提取，结果如图3所示。由图3可知，一阶、二阶算子以及传统形态学手段处理后仍然存在云团界限不分明、内部点线杂乱、云团高度零散、边界较粗等问题。实际上，原图中原本集中的云层在二值化处理的过程中由于区域边界线噪音干扰，形成了较多较为零散的云层，上述几种算法均无法降低该噪音的干扰，而Bwperim-Imfill组合算子在Bwperim算子基础上进行填充腐蚀操作，Imfill算子可实现在多次腐蚀过程中将原本属同一云团而78程麒铭等：基于气象雷达图信息处理的短时降雨预报方法第3期误分开的云团重新合并，降低噪音干扰，恢复原有云团样貌，提取清晰边缘界限。2.2云团追踪结果在MATLAB中采用2021年9月6日秀山县一场完整的实测降雨数据对算法进行验证，其中，坐标系原点为图上左下角，研究区域中心点位于图中坐标（220,394），研究区域为中心点d85之内所有区域。经过边缘检测算法提取后的每个云团均为闭合区域，其边界所围成像素点个数即为面积，在实际情况下，面积较小的云团很可能消散或被合并，易干扰预测结果。本文将云团面积阈值取20，剔除低于该阈值的云团，最终处理效果如图4。由于算法识别到云团数量不相等，故采用整体质心法进行计算，移动速度采用该时刻与T时刻速度均值，V=1.67 像素点/min，=19.5，运动方向上云团距研究区域中心点距离为23.3像素点。坐标变换后，yi,min,yi,max所在点满足yi,miny0yi,max，则降雨能发生，在图5中以红色圆圈表示，研究区域以红色矩形表示，注：（a）Prewitt算子；（b）Sobel算子；（c）Canny算子；（d）传统形态学方法；（e）Bwperim算子；（f）Bwperim-Imfill组合算子。图3不同算法下的风暴云图边缘检测结果Fig.3 Edge detection results of storm cloud image using differentalgorithms（a）（b）（c）（d）（e）（f）（a）（b）基本反射率雷达站名：铜仁观测日期：2021-09-06观测时间：18:31:30BJT数据范围：125 km显示仰角：1.570dBZ65605550454035302520151050基本反射率雷达站名：铜仁观测日期：2021-09-06观测时间：18:37:17BJT数据范围：125 km显示仰角：1.570dBZ65605550454035302520151050注：（a）风暴云团质心为（148.3,309.7）；（b）风暴云团质心为（149.9,309.2）。图4雷达图中风暴云团追踪结果Fig.4Results of storm cloud tracking in radar map79第43卷水 文离研究区域最近云团边界点以*号表示。图5降雨发生检测结果Fig.5Result of rainfall detecting2.3降雨时间与降雨强度预测结果根据云团追踪算法，得到降雨时间Train=14min，此时预测降雨时间在图4（b）时刻的14min之后，雷达图刷新频率为6min/张，意味着此过程中还有2张图片。由于云团时刻发生变化，若直接采用该时间作为最终预测结果，可能误差较大，故采用迭代递推方式，每6 min进行一次预测，逐渐逼近降雨过程，结果如表1所示。预测19：00时刻风暴云团刚好到达研究区域，认为此刻已经降雨，降雨前6 min与12 min预测的降雨发生时间均值为19：00。根据秀山县雨量计实时反馈结果可知，降雨实际发生在19：00附近，且19：00至19：05实际降雨强度为54 mm/h，19：05至19：10实际降雨强度为 66 mm/h。降雨预测结果表明，降雨发生前6 min的预报结果更接近于实际测量值，故在发生降雨后，继续采用算法进行迭代递推，可较好预测出降雨过程（见图6），且降雨预测模型的纳什效率系数（NSE）为0.71，均方根误差（RMSE）为12.05 mm/h，可满足模型应用精度要求。70605040302010019:00 19:1019:20 19:3019:40 19:50 20:0020:10 20:2020:30 20:4070503010-10-30-50实测降雨强度预测降雨强度绝对误差降雨强度/（mmh-1）降雨强度/（mmh-1）-70图6降雨过程预测与实测结果Fig.6Prediction and measured results of rainfall process由图6知，绝对误差最大值发生在19：37附近，可能是因为研究区域在雷达图中的选取无法完全吻合，导致本已移出研究区域的云团也被计算在内，产生一定误差。实测与预测的最大降雨强度均发生在19：10附近，实测最大降雨强度为 66 mm/h，而预测值为59.4 mm/h，误差仅为10%，预测20：34降雨结束，而实测20：20降雨结束，两者相差16 min。实际上，本研究在研究区域设置的雨量计最低检测雨量为0.5 mm，实测降雨过程在20：20至21：00累计雨量并无变化，而在21：00雨量计的雨量增加0.5 mm，表明在该时刻内有降雨发生，平均降雨强度在6 mm/h以内，与预测降雨强度基本符合。3结论与展望针对传统降雨预报方法误差大，且需大量气象资料等问题，研究提出了一种基于雷达图信息处理的短时降雨预报方法。结果证实，在气象数据匮乏情况下，采用经典Z-R关系式模型可实现降雨过程预测，预测最大降雨强度误差仅为10%，预测降雨开始时间和最大降雨强度发生时间与实测基本一致，受雨量计检测限限制预测降雨结束时间与实测数据相差16 min。本预测方法的NSE高达0.71，RMSE为12.05mm/h，可满足应用精度要求，且适用于单站雷达覆盖所有区域。限于研究条件与时间，本文收集到的雷达图与实测数据有限，仅验证一场降雨可能存在偶然误差，未雷达图显示原始时刻18：37；18：4318：43；18：4918：49；18：5418：54；19：0018：31；18：37移动速度（像素/min）1.010.981.041.061.67Train/min1489014预测降雨发生时刻18：5718：5719：0319：0018：51预测降雨强度/（mmh-1）Train+041.134.253.532.235.2Train+554.739.852.751.639.8Train+1045.158.753.159.356.3表1基于风暴云团追踪的降雨预测结果Table1Rainfall prediction results based on storm cloud tracking80程麒铭等：基于气象雷达图信息处理的短时降雨预报方法第3期来将对不同雨型进行预测。同时采用的Z-R关系仅为经典公式，若研究区域气象资料充足，可利用历史气象资料拟合出相应的Z-R关系式，也可采用BP神经网络、遗传算法等拟合方法，进一步提高预报精度。针对研究区域在雷达图中无法完全吻合的问题可在收集足够资料后采用卷积神经网络进行解决，进一步提高预报精度。参考文献：1COLE S J,MOORE R J.Hydrological modelling using raingauge-and radar-based estimators of areal rainfall J.Journal of Hydrology,2008,358：159-181.2OCHOA-RODRIGUEZ S,WANG L P,GIRES A,et al.Impact of spatial and temporal resolution of rainfall inputs on urban hydrodynamicmodelling outputs：A multi-catchment investigation J.Journal ofHydrology,2015,531：389-407.3LOOPER J P,VIEUX B E.An assessment of distributed flash floodforecasting accuracy using radar and rain gauge input for a physics-based distributed hydrologic model J.Journal of Hydrology,2012,412：114-132.4OCHOA RODRIGUEZ S,WANG L P,WILLEMS P,et al.A Reviewof Radar Rain Gauge Data Merging Methods and Their Potential forUrban Hydrological Applications J.Water Resources Research,2019,55（8）：6356-6391.5戴铁丕,张培昌,魏鸣.713测雨雷达测定区域降水量初探 J.南京气象学院学报,1987,10（1）：87-94.6张培昌,顾松山,王振会,等.气象万千探本索源南京信息工程大学“大气探测学科”发展历程回顾与展望 J.大气科学学报,2021,44（2）：165-177.7赵然杭,甘甜,逄晓腾,等.基于时间序列分解的降雨数据挖掘与预测 EB/OL.http：/ M.北京：气象出版社,2000.9刘畅.基于多普勒雷达的风暴云团的结构及形态外推方法研究 D.天津：天津大学,2012.10庄旭东,胡胜,陈荣,等.“雨燕”中风暴算法与新一代雷达SCIT产品的对比分析 J.热带气象学报,2011,27（3）：299-306.11黄兴友,马玉蓉,胡苏蔓.基于深度学习的天气雷达回波序列外推及效果分析 J.气象学报,2021,79（5）：817-827.12吴亚昊,刘黎平,周筠珺,等.雨滴谱的变化对降水估测的影响研究 J.高原气象,2016,35（1）：220-230.13虎玲,常霞,纪峰.图像边缘检测方法研究新进展 J.现代电子技术,2018,41（23）：40-45.14张红霞,王灿,刘鑫,等.图像边缘检测算法研究新进展 J.计算机工程与应用,2018,54（14）：11-18.15付永庆,王咏胜.一种基于数学形态学的灰度图像边缘检测算法J.哈尔滨工程大学学报,2005,26（5）：685-687.16欧先锋,贾振红,郝军.多结构抗噪膨胀-腐蚀型数字灰度图像边缘检测的研究 J.激光杂志,2009,30（1）：40-41.17邓仕超,黄寅.二值图像膨胀腐蚀的快速算法 J.计算机工程与应用,2017,53（5）：207-211.18蔡利梅.MATLAB图像处理：理论、算法与实例分析 M.北京：清华大学出版社,2020.19涂盼盼,何毅斌,熊文豪,等.一种基于数学形态学的齿轮边缘检测方法研究 J.机电工程,2019,36（8）：839-841+850.20SOILLE,PIERRE.Morphological Image Analysis：Principles and Applications M.Berlin：Springer-Verlag,1999.21李玲.基于雷达数据的短临降雨预报方法和应用研究 D.成都：成都理工大学,2020.22张鹏,刘西川,周则明,等.基于实测雨滴谱数据的微波链路和天气雷达降水估计关系研究 J.气象,2021,47（7）：843-853.23姚铖.基于Z-R关系优化的雷达降雨估测及水文应用研究 D.大连：大连理工大学,2021.Short-Term Rainfall Forecasting Method Based on Information Processing of Weather Radar MapCHENG Qiming1，CHEN Yao1,2，LIU Fei1,2，MAO Jun1，LIU Linhao1，TANG Jie1（1.School of River and Ocean Engineering,Chongqing Jiaotong University,Chongqing400074,China;2.Engineering Laboratory of Environmental Hydraulic Engineering of Chongqing MunicipalDevelopment and Reform Commission,Chongqing Jiaotong University,Chongqing400074,China）Abstract：The storm clouds shown on weather radar map are very irregular,which leads to the problems of large cumulative error,easy loss of cloud tracking and strong dependence of meteorological data when using traditional rainfall forecasting methods.Therefore,the forward binarization method was used to preprocess the original radar map.The Bwperim-Imfill operator could beused to effectively detect the edges of scattered and irregular clouds in the radar map.Based on this,the storm cloud trackingmodel proposed in this study can identify the direction of cloud movement,and then estimate the time to arrive at the study area（i.e.,the time of rainfall）.The results showed that,under the condition of severe shortage of meteorological data,rainfall prediction can be achieved by using the storm cloud tracking model coupling unfitted classical Z-R relation.The prediction error ofmaximum rainfall intensity was only 10%,the prediction error of rainfall occurrence time was 3 min,and the prediction error ofrainfall end time was 16 min.The Nash efficiency coefficient（NSE）was as high as 0.71,root mean square error（RMSE）was12.05mm/h,which meets the application accuracy requirements of the model.Keywords：radar map；storm clouds；edge detection；cloud tracking；Z-R formula；rainfall intensity81