[发明专利]NRIET blending融合预报方法

摘要

本发明公开了一种NRIET blending融合预报方法，包括以下内容：利用多尺度光流法计算相邻时刻雷达组网的移动矢量，基于移动矢量，利用半拉格朗日方案进行外推计算，得到外推回波预报场。利用区域快速循环同化模式，同化雷达资料，并进行数值预报，预报时长为6小时，得到反演的雷达回波。对模拟雷达回波进行插值，得到10分钟一次，空间分辨率1km的预报雷达回波。利用融合预报方法，对短临外推及数值预报结果进行融合，得到融合预报回波场。利用CSI检验评分，得到回报预报的评分结果；可以基于评分结果，调整融合参数。

主权项

1.一种NRIET blending融合预报方法，其特征在于，包括以下内容：步骤1、利用多尺度光流法获得外推风场，结合半拉格朗日外推算法，生成回波和降水外推产品；1.1.短临外推方法：光流法可以根据不同时刻的回波图像获得雷达回波的移动矢量，结合外推计算，得到雷达回波的外推预报；图像中所有像素点的光流就构成了图像的光流场，而光流法的核心正是从连续的图像序列中计算光流场；设图像的亮度在短时间内移动过程中不发生变化，即：I(x+uΔt，y+vΔt，t+Δt)＝I(x，y，t)    (1)其中I是图像亮度，x，y为像素点的空间坐标位置，u，v为对应移动方向的移动速度，t为当前时间，Δt为时间间隔；得到：Ixu+Iyv+It＝0    (2)Ix，Iy为图像亮度的空间微分，It为图像亮度的时间微分，即：由于(2)式中有u、v两个变量，因此需要额外的条件约束，才能得到完整的光流场；采用的办法是通过最小二乘法，设计算出的风场误差平方和最小时为最优解；光流法假定的前提条件是微小的运动，当运动速度较快时，移动的速度与移动位置的偏导存在较大偏差；为了解决这个问题，需要采用多尺度法，将图像分为不同分辨率，从最粗的分辨率开始估算移动位置，作为下一层的初始位置，依次向下搜索，这种方法被称为金字塔分层；我们采用Gunnar Farneback光流法和3层金字塔分层技术来获得雷达回波的移动矢量，其中每层分辨率的比例为2；1.2.半拉格朗日外推法利用多尺度光流法得到雷达回波的移动矢量后，根据移动矢量进行外推预报；采用的外推方法为线性外推法，即雷达回波未来移动的距离等于当前时刻移动矢量乘以预报时间；因此，对于单个像素点而言，在预报时间内仅能沿直线运动，而实际中回波的移动路径往往具有一定的曲率，利用线性外推无法有效的反映这一特征，从而导致外推预报的失败；考虑到线性外推的不足，我们提出利用半拉格朗日方法进行外推；外推算法的基本表达式如下：F(t0+τ，x)＝F(t0，x‑a)    (4)其中F为回波强度，t为时间，Δt为时间间隔，x为空间位置，a为空间间隔，即当前位置未来时刻的回波是有当前时候特定位置上的回波平移得到；不同外推算法的意义在于如何建立二者之间的联系；考虑二维情况，当雷达回波从(x‑2a，y‑2β)移动到(x，y)，有：其中a,β为Δt时间内，雷达回波在x，y方向上的移动距离，因此：a＝ΔtU(x‑a，y‑β，t)   (6)β＝ΔtU(x‑a，y‑β，t)   (7)其中u(x‑a，y‑β，t),v(x‑a，y‑β，t)分别为x，y方向上x‑a位置、y‑β位置、t时刻的移动速度，考虑雷达回波的移动速度随着空间在变化，因此可以将整个外推预报时段分成更小的时间段，采用迭代的方法来获得每个时间步长的移动距离，k为迭代次数：ak+1＝ΔtU(x‑ak，y‑βk，t)   (8)βk+1＝ΔtV(x‑ak，y‑βk，t)   (9)最后将各个步长的移动距离累加，即可得到总的移动距离；利用短临外推算法，得到0‑2小时特定高度的雷达回波外推预报结果，空间分辨率为1km，时间分辨率为10分钟；步骤2.利用区域数值模式，得到0‑6小时数值预报结果；基于数值预报场，得到反演的相同高度雷达回波，要求时间分辨率也为10分钟；中尺度数值模式采用WRF模式；同化模块采用格点统计插值，以下简称GSI；GSI同化系统除了能够同化常规观测资料之外，特别地是在同化各种非常规观测资料方面的能力表现尤为突出，能为数值模式提供接近真实情况的初始场，从而进一步提升数值预报模式的准确度。随着该系统的不断完善升级和计算资源的优化，GSI将来将具备四维变分同化和观测敏感性分析的能力，同时与集合预报系统相结合，可以实现集合‑变分混合同化技术；背景场wrfinput可由两种方式获得，一种是经过WPS前处理模块得到，又称模式的冷启动过程；另一种是由上一次模式的预报结果提供，这种方式称为循环模式Cycling，即模式的热启动过程；对于GSI同化系统，观测场资料须存为二进制通用气象数据格式，即BUFR，统计控制数据是背景误差协方差和观测误差协方差及与同化有关的参数信息；以上三部分的输入数据进入GSI同化模块，并由用户作调节设置，得到分析场及各种诊断结果；最后利用分析场和经由WPS模块得到的边界条件wrfbdy更新模式区域的边界条件；最后将同化后得到分析场作为预报的初始场文件wrfinput，连同更新的边界条件wrfbdy，输入WRF预报模块进行模式的积分运算；考虑到计算资源和数据时效性的问题，采用的同化方案中只更新WRF预报模块的初始场，边界条件采用NCEP全球预报系统GFS数值产品提供的边界场，不进行实时更新；步骤3.对区域数值预报的雷达回波进行空间插值，利用双线性插值，得到与短临外推相同分辨率的预报产品；利用数值预报的输出产品：雨水、雪、霰的混合比，计算等效的雷达回波；其中假定水成物粒子均为球形，粒子浓度与尺度的关系假定为e指数变化，其公式为：n(D)＝n0e‑λD,其中n(D)为粒子浓度，D为粒子直径，n0及λ为谱参数，n0可取8*103,λ取值4.1R‑0.21，其中R为降雨率；外推预报的空间分辨率取决于雷达组网回波的分辨率；对于S及C波段天气雷达，组网分辨率为1km；利用双线性插值方法，调整区域数值预报的雷达回波产品分辨率；双线性插值流程如下：令f为各个格点的值，数值模式格点为Q11，Q12，Q21，Q21，对应坐标为(x1,y1)、(x1,y2)、(x2,y1)、(x2,y2)，对应的值为f(Q11),f(Q12),f(Q21),f(Q22)，目标位置为P，坐标为(x,y)；先对X方向进行插值，得到R1，R2点的值f(R1)，f(R2)：然后在Y方向插值，得到：其中f(P)即为观测格点上的模式结果；按照上述方法，对数值模式模拟雷达回波进行插值处理，使模式模拟雷达回波分辨率与观测保持一致；步骤4.计算融合权重：对短临外推和模式计算回波进行融合处理，融合方法为加权平均；由于短临外推预报的评分随时间呈指数下降，因此构建融合权重时也考虑采用指数型；公式为：其中w为权重，t为时间，α和β分别代表0、6小时的预报评分，需要利用一段时间的评测来确定具体数值，t为时间，单位为小时；考虑在融合预报中，短临外推和数值模式有可能会在同一位置上产生较大的差异，此时简单的考虑时间权重，可能会损失一部分有用信息，因此考虑引入回波强度的显著性指标：其中N₁和N₂代表归一化的短临外推雷达回波和数值预报雷达回波强度；为累积概率密度函数，当N₁(x，y)‑N₂(x，y)趋近于1或者‑1，相应的r取值为1或0，代表短临外推与数值预报的差异越大；引入显著性指标后，联合时间权重w(t)，构建新的权重函数：该函数统筹兼顾时间权重w和显著性权重r，从而构建二维的权重函数；综上，最终的融合算法表达式为:S＝ws*I1+(1‑ws)*I2，其中I1代表短临外推的雷达回波，I2代表数值预报的雷达回波，ws为I1的融合权重，(1‑ws)为I2的融合权重，则S为融合预报的雷达回波；步骤5.利用回波分档CSI评分，检验评估融合预报效果；CSI计算方法为：其中N_r为预报准确的点，N_e为误报，N_m为漏报；利用评分结果，对融合参数进行评估，从而动态调整α和β等参数。