CN108549116A - NRIET blending融合预报方法 - Google Patents

Abstract

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

Description

NRIET blending融合预报方法

技术领域

本发明涉及一种NRIET blending融合预报方法，属于气象预报系统。

背景技术

近年来气象预报格点化、精细化逐渐成为发展趋势。精细化格点预报主要依靠数值预报，但数值模式在起报阶段往往表现较差，因而短临预报作为有效补充手段必不可少。

大风、雷暴、强降水等强对流天气是造成航运中断、飞机延误的主要原因之一，会给运输、民航等行业造成巨大损失。而利用短临预报技术有助于提前调整运输路线和航班，减少社会成本和经济损失。

利用短临外推的方法仅能解决0-2小时的预报问题，数值模式的表现可在3-6小时内达到稳定，因此需要融合预报方法，结合二者优势，形成覆盖0-6小时的预报能力。近年来，融合短临预报以其兼具短临外推和数值预报的优势逐渐成为短临预报的发展方向，具有很高的应用前景和研究价值。融合预报主要用于解决雷暴、降水和大风的0-6小时短临预报。

目前国外主流的融合短临预报系统包括：

美国Niwot-human短临系统

美国CoSPA短临系统

英国NIMROD短临系统

英国、澳大利亚STEPS短临系统

法国SIGOONS短临系统

日本VSRF短临系统

这些系统预报时效基本为0-6小时，预报对象主要为降水和灾害性天气，预报要素包括降水、风、雷达回波等。

其中，NIMROD系统采用模式与外推系统线性加权平均的方式进行融合预报，外推预报融合权重随时间指数下降，而数值模式权重不变。S-STEPS系统在外推和模式之外，考虑了噪音的影响，同时外推预报融合权重根据预报技巧实时更新。Niwot-human系统中将模式风场和倾向以一定比例加入到外推中，从而使得外推预报除了位置变化外，具备一定的强度预报信息。CoSPA是美国航空局最新研制的短临预报系统，主要目的是解决雷暴的短临预报问题。其融合权重在线性变化的基础上，还考虑了雷达回波的强度信息。

当前，国内的融合短临预报发展落后于国外，主要表现在缺乏先进的外推算法及融合预报算法，目前业务化的融合短临系统主要有两个：

中国气象局GRAPES-SWIFT短临系统

香港RAPIDS短临系统

其中GRAPES-SWIFT系统是由中国气象局研制，主要服务于气象预报业务领域。其主要融合方式是通过模式和外推的平流风场进行融合以预报雷暴的移动，同时模式给出风暴的环境潜势以预报雷暴的强度。其主要的预报变量为风和雷达回波。

RAPIDS短临预报系统是香港机场使用的短临系统。其采用外推预报和模式预报加权方式进行融合，权重随时间呈指数变化，其预报对象为雷达回波和降水。

除GRAPES-SWIFT、RAPIDS系统外，国内其他各省、市气象局、民航气象系统、水利等部门及行业尚无成熟的短临预报系统，但精细化、格点化的0-6小时短临预报对其业务开展有重要意义。因此，近年来国内短临预报技术的研究受到较多关注，但融合预报技术仍属于较为前沿的领域。

发明内容

针对以上问题本发明提供了一种NRIET blending融合预报方法，该融合预报的目的是结合短临外推与数值预报的优势，形成覆盖0-6小时的最优预报结果，其主要面临的技术难题包括：

1.如何构建短临外推预报？

根据连续时次雷达回波产品，生成雷达回波移动矢量，假定在很短时间步长内，天气系统的移动速度和方向不变，在此基础上进行外推预报。

2.如何搭建区域快速循环同化数值模式？

建立0-6小时的短临预报依赖于区域数值模式。为了解决数值模式起报阶段表现不佳的问题，需要利用资料同化方法改进数值模式初始场。

3.外推与数值模式输出场的时空分辨率如何匹配？

短临外推基于雷达回波，水平分辨率为1km，而数值模式的分辨率各不相同，需要进行统一。短临外推的时间分辨率可以达到10分钟，而数值预报的输出频次不固定，需要考虑如何进行统一。

4.短临外推和数值模式的融合权重如何确定？

短临外推和数值预报的融合权重是决定融合预报效果的关键，因此需要给出一套具体算法，解决融合权重的问题。

5.融合效果如何评判？

基于融合的预报场，需要给出评估融合预报效果的方法，并对融合预报技术进行检验。

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

1.利用多尺度光流法获得外推风场，结合半拉格朗日外推算法，生成回波和降水外推产品。

1.1.短临外推方法：

光流法可以根据不同时刻的回波图像获得雷达回波的移动矢量，结合外推计算，得到雷达回波的外推预报，因而在临近预报中具有重要意义。光流法是计算机视觉领域中的重要方法。光流的概念由Gibson于1950年首先提出，它是指由于被观测的目标和传感器之间的相对运动，而在序列图像中产生的瞬时位移场，体现了图像亮度模式的表观运动。图像中所有像素点的光流就构成了图像的光流场，而光流法的核心正是从连续的图像序列中计算光流场。

假定图像的亮度在短时间内移动过程中不发生变化，即：

I(x+uΔt，y+vΔt，t+Δt)＝I(x，y，t) (1)

其中I是图像亮度，x，y为像素点的空间坐标位置，u，v为对应移动方向的移动速度，t为当前时间，Δt为时间间隔。

得到：

I_xu+I_yv+I_t＝0 (2)

I_x，I_y为图像亮度的空间微分，I_t为图像亮度的时间微分，即：

由于(2)式中有u、v两个变量，因此需要额外的条件约束，才能得到完整的光流场。通常采用的办法是通过最小二乘法，假定计算出的风场误差平方和最小时为最优解。

光流法假定的前提条件是微小的运动，当运动速度较快时，移动的速度与移动位置的偏导存在较大偏差。为了解决这个问题，需要采用多尺度法，将图像分为不同分辨率，从最粗的分辨率开始估算移动位置，作为下一层的初始位置，依次向下搜索，这种方法被称为金字塔分层。

我们采用Gunnar Farneback光流法和3层金字塔分层技术来获得雷达回波的移动矢量，其中每层分辨率的比例为2。

1.2.半拉格朗日外推法

利用多尺度光流法得到雷达回波的移动矢量后，可以根据移动矢量进行外推预报。通常采用的外推方法为线性外推法，即雷达回波未来移动的距离等于当前时刻移动矢量乘以预报时间。因此，对于单个像素点而言，在预报时间内仅能沿直线运动，而实际中回波的移动路径往往具有一定的曲率，利用线性外推无法有效的反映这一特征，从而导致外推预报的失败。

考虑到线性外推的不足，我们提出利用半拉格朗日方法进行外推。外推算法的基本表达式如下：

F(t₀+τ，x)＝F(t₀，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为迭代次数：

a^k+1＝ΔtU(x-a^k，y-β^k，t) (8)

β^k+1＝ΔtV(x-a^k，y-β^k，t) (9)

最后将各个步长的移动距离累加，即可得到总的移动距离。

利用短临外推算法，得到0-2小时特定高度的雷达回波外推预报结果，空间分辨率为1km，时间分辨率为10分钟。

2.利用区域数值模式，得到0-6小时数值预报结果。基于数值预报场，得到反演的相同高度雷达回波，要求时间分辨率也为10分钟；

中尺度数值模式采用WRF模式(Weather Research and Forecasting Model)。同化模块采用格点统计插值(Gridpoint Statistical Interpolation，以下简称GSI)

业务应用及科研结果均证实GSI同化系统除了能够同化常规观测资料之外，特别地是在同化各种非常规观测资料(包括来自N0AA15、16、17、18、19及METOP、AQUA等卫星的辐射率观测资料、AIRS多通道大气红外辖射资料、雷达径向风资料及GPSRO无线电掩星观测资料等)方面的能力表现尤为突出，能为数值模式提供接近真实情况的初始场，从而进一步提升数值预报模式的准确度。随着该系统的不断完善升级和计算资源的优化，GSI将来将具备四维变分同化和观测敏感性分析的能力，同时与集合预报系统相结合，可以实现集合-变分混合同化技术。

背景场(wrfinput)可由两种方式获得，一种是经过WPS前处理模块得到，又称模式的冷启动过程；另一种是由上一次模式的预报结果提供，这种方式称为循环模式(Cycling)，即模式的热启动过程。对于GSI同化系统，观测场资料须存为二进制通用气象数据格式，即BUFR(Binary Universal Form for the Representation of meteorologicaldata)，统计控制数据主要是背景误差协方差和观测误差协方差及与同化有关的参数信息。以上三部分的输入数据进入GSI同化模块，并由用户作调节设置，得到分析场及各种诊断结果。最后利用分析场和经由WPS模块得到的边界条件(wrfbdy)更新模式区域的边界条件。最后将同化后得到分析场作为预报的初始场文件(wrfinput)，连同更新的边界条件(wrfbdy)，输入WRF预报模块进行模式的积分运算。考虑到计算资源和数据时效性的问题，本项目中采用的同化方案中只更新WRF预报模块的初始场，边界条件采用NCEP全球预报系统(GFS)数值产品提供的边界场，不进行实时更新。

3.对区域数值预报的雷达回波进行空间插值，利用双线性插值，得到与短临外推相同分辨率的预报产品；

利用数值预报的输出产品：雨水、雪、霰的混合比，可以计算等效的雷达回波。其中假定水成物粒子均为球形，粒子浓度与尺度的关系假定为e指数变化，其公式为：n(D)＝n₀e^-λD,其中n(D)为粒子浓度，D为粒子直径，n₀及λ为谱参数，n₀可取8*10³,λ取值4.1R^-0.21，其中R为降雨率。

外推预报的空间分辨率取决于雷达组网回波的分辨率。对于S及C波段天气雷达，通常的组网分辨率约为1km。利用双线性插值方法，可以调整区域数值预报的雷达回波产品分辨率。

双线性插值流程如下：

令f为各个格点的值，数值模式格点为Q₁₁，Q₁₂，Q₂₁，Q₂₁，对应坐标为(x₁,y₁)、(x₁,y₂)、(x₂,y₁)、(x₂,y₂)，对应的值为f(Q₁₁),f(Q₁₂),f(Q₂₁),f(Q₂₂)，目标位置为P，坐标为(x,y)。先对X方向进行插值，得到R₁，R₂点的值f(R₁)，f(R₂)：

然后在Y方向插值，得到：

其中f(P)即为观测格点上的模式结果。按照上述方法，对数值模式模拟雷达回波进行插值处理，使模式模拟雷达回波分辨率与观测保持一致。

4计算融合权重：

对短临外推和模式计算回波进行融合处理，融合方法为加权平均。

由于短临外推预报的评分随时间呈指数下降，因此构建融合权重时也考虑采用指数型。公式为：其中w为权重，t为时间，α和β分别代表0、6小时的预报评分，需要利用一段时间的评测来确定具体数值，t为时间，单位为小时。

考虑在融合预报中，短临外推和数值模式有可能会在同一位置上产生较大的差异，此时简单的考虑时间权重，可能会损失一部分有用信息，因此考虑引入回波强度的显著性指标：其中x,y为空间位置，N₁和N₂代表归一化的短临外推雷达回波和数值预报雷达回波强度。为累积概率密度函数，当N₁(x，y)-N₂(x，y)趋近于1或者-1，相应的r取值为1或0，代表短临外推与数值预报的差异越大。引入显著性指标后，联合时间权重w(t)，可以构建新的权重函数：该函数可以统筹兼顾时间权重w和显著性权重r，从而构建二维的权重函数。

综上，最终的融合算法表达式为:S＝w_s*I₁+(1-w_s)*I₂，其中I₁代表短临外推的雷达回波，I₂代表数值预报的雷达回波，w_s为I₁的融合权重，(1-w_s)为I₂的融合权重，则S为融合预报的雷达回波。

5.利用回波分档CSI评分，检验评估融合预报效果。CSI计算方法为：其中N_r为预报准确的点，N_e为误报，N_m为漏报。利用评分结果，对融合参数进行评估，从而动态调整α和β等参数。

本发明相对于最接近的现有技术而言，有以下有益效果：利用多尺度光流法和半拉格朗日短临外推算法可以构建精细化、客观的短临预报雷达回波，时效为0-2小时。利用WRF和GSI模式可以进行区域快速循环同化预报，时效为0-6小时，其中预报技巧在3-6小时达到平稳。利用融合算法，可以在0-2小时时段偏重短临外推，2-6小时偏向数值预报，从而形成覆盖0-6小时的最优预报结果。从整体来看，预报效果要优于数值预报及短临外推至少3-5％。

附图说明

图1为NRIET blending融合预报方法流程图。

具体实施方式

下面结合附图1对本发明作具体说明。

系统需要的数据主要包括雷达组网数据及GFS全球背景场数据。其中雷达组网数据用来作为短临外推初始场及同化资料进入数值模式，GFS提供数值模式初始场及边界条件。

如图1所示，本发明提供了一种NRIET blending融合预报方法包括以下内容：

1.利用多尺度光流法基于相邻时刻雷达组网数据，选取3km高度，进行移动矢量的计算；基于计算得到的移动矢量，利用半拉格朗日方法进行外推计算，外推时间为2小时，间隔10分钟输出一次，空间分辨率为1km；

2.基于WRF，将GFS全球背景场降尺度为模式初始场及边界条件，利用GSI同化系统，将雷达资料同化进模式初始场；利用同化初始场进行预报，预报时间为6小时，空间分辨率为3km；基于数值模式输出结果，进一步反演得到模拟雷达回波；

3.经过空间双线性插值，得到3km高度对应观测回波空间位置的模拟回波；

4.利用Blending融合算法，对短临外推回波及数值模式模拟回波进行融合，得到最优融合预报结果；

5.对融合结果进行评分检验，并继续检验结果，动态调整融合参数。

得到融合结果在0-6小时内平均高于短临外推与数值预报，其中0-1小时与外推评分相近，平均比数值预报评分高约30％；2-6小时与数值预报评分相近，评分远高于短临外推；1-2小时数值预报与短临预报评分相近，此时融合预报评分比二者高约10％。

以上所述仅为本发明的优选实施例而已，并不限制于本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (2)

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为时间间隔；

得到：

I_xu+I_yv+I_t＝0 (2)

I_x，I_y为图像亮度的空间微分，I_t为图像亮度的时间微分，即：

由于(2)式中有u、v两个变量，因此需要额外的条件约束，才能得到完整的光流场；采用的办法是通过最小二乘法，设计算出的风场误差平方和最小时为最优解；

光流法假定的前提条件是微小的运动，当运动速度较快时，移动的速度与移动位置的偏导存在较大偏差；为了解决这个问题，需要采用多尺度法，将图像分为不同分辨率，从最粗的分辨率开始估算移动位置，作为下一层的初始位置，依次向下搜索，这种方法被称为金字塔分层；

我们采用Gunnar Farneback光流法和3层金字塔分层技术来获得雷达回波的移动矢量，其中每层分辨率的比例为2；

1.2.半拉格朗日外推法

利用多尺度光流法得到雷达回波的移动矢量后，根据移动矢量进行外推预报；采用的外推方法为线性外推法，即雷达回波未来移动的距离等于当前时刻移动矢量乘以预报时间；因此，对于单个像素点而言，在预报时间内仅能沿直线运动，而实际中回波的移动路径往往具有一定的曲率，利用线性外推无法有效的反映这一特征，从而导致外推预报的失败；

考虑到线性外推的不足，我们提出利用半拉格朗日方法进行外推；外推算法的基本表达式如下：

F(t₀+τ，x)＝F(t₀，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为迭代次数：

a^k+1＝ΔtU(x-a^k，y-β^k，t) (8)

β^k+1＝ΔtV(x-a^k，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)＝n₀e^-λD,其中n(D)为粒子浓度，D为粒子直径，n₀及λ为谱参数，n₀可取8*10³,λ取值4.1R^-0.21，其中R为降雨率；

外推预报的空间分辨率取决于雷达组网回波的分辨率；对于S及C波段天气雷达，组网分辨率为1km；利用双线性插值方法，调整区域数值预报的雷达回波产品分辨率；

双线性插值流程如下：

令f为各个格点的值，数值模式格点为Q₁₁，Q₁₂，Q₂₁，Q₂₁，对应坐标为(x₁,y₁)、(x₁,y₂)、(x₂,y₁)、(x₂,y₂)，对应的值为f(Q₁₁),f(Q₁₂),f(Q₂₁),f(Q₂₂)，目标位置为P，坐标为(x,y)；先对X方向进行插值，得到R₁，R₂点的值f(R₁)，f(R₂)：

然后在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＝w_s*I₁+(1-w_s)*I₂，其中I₁代表短临外推的雷达回波，I₂代表数值预报的雷达回波，w_s为I₁的融合权重，(1-w_s)为I₂的融合权重，则S为融合预报的雷达回波；

步骤5.利用回波分档CSI评分，检验评估融合预报效果；CSI计算方法为：其中N_r为预报准确的点，N_e为误报，N_m为漏报；利用评分结果，对融合参数进行评估，从而动态调整α和β等参数。

2.根据权利要求1所述的NRIET blending融合预报方法，其特征在于，所述的非常规观测资料包括来自N0AA15、16、17、18、19及METOP、AQUA的辐射率观测资料、AIRS多通道大气红外辖射资料、雷达径向风资料及GPSRO无线电掩星观测资料。