一种基于物理方程的降水外推预报方法与流程

1.本发明属于降水预报技术领域，特别是涉及一种基于物理方程的降水外推预报方法。


背景技术：

2.降水预报是气象相关部门最重要的业务服务之一。降水预报的好坏不仅与人们的日常生活息息相关，而且影响着各级政府开展的防灾减灾工作。由于降水生效在时空上的不确定性、降水系统演变在物理过程上的复杂性，当前降水预报的准确率相比温度、气压、风等其它气象要素依然偏低。因此，如何有效提高降水预报水平一直是气象领域的重点和难点问题。
3.一般地，针对不同时效，降水预报的主要方法有所差别。对于时效超过6小时的短期预报而言，降水预报主要依托数值预报模式开展，通过模式数值积分大气动力学方程组直接给出未来降水的预测。而对于0~6小时的短临预报而言，降水预报主要依靠统计外推技术实现。这一技术的基本思想是通过比较连续时次的降水量空间分布图像，识别不同的降水系统并追踪其移动的水平矢量，进而外推得到未来时次的降水分布。当前，在降水短临预报领域主流的统计外推技术主要包括光流（optical flow）法、追踪（trec）法等。然而，不论以上哪种方法都是在降水没有生效的假设上应用的，即只考虑降水系统的水平移动、忽略新降水系统的生成和现有降水系统的消亡。这不仅会导致新生降水系统预报不出来，而且会引起现有降水系统一直移动、无法消散的问题。


技术实现要素：

4.针对以上问题，本发明旨在提供一种基于物理方程的降水外推预报方法，从动力学方程出发，在方程中直接包含降水的平移和生效的贡献，最终发展更合理有效的降水短临预报技术。
5.本发明提供一种基于物理方程的降水外推预报方法，具体步骤如下，其特征在于：步骤一：获取当前时刻t0和其前两个时次t
‑1、t
‑2的降水量r数据r0、r
‑1、r
‑2；步骤二：对r0、r
‑1、r
‑2分别进行二维平滑处理；步骤二中的平滑公式为：
ꢀꢀ
（1）其中z为包含n个样本一维空间序列，平滑时，先在x方向进行一维平滑，再在y方向
进行一维平滑；步骤三：应用带常数项的多元线性回归方法，估计出降水移动速度项和降水生消强迫项，其中u表示x方向移动速度、v表示y方向移动速度、f表示降水生消强迫；步骤三中的回归方程为：（2）其中i和j分别为x方向和y方向的坐标，r
‑
1.5
、r
‑
0.5
分别表示t
‑1和t
‑2、t0和t
‑1的平均降水场，即：（3）其中和分别表示两个时次平均的x方向移动速度和y方向移动速度；步骤四：利用拉格朗日向前积分方案对r0进行外推，得到未来时刻t1的降水平移预报场rs1；步骤四中的外推公式为：（4）其中：（5）由于外推得到的降水平移预报场所处格点不一定正好位于原网格点i，j上，所以利用kriging插值方法将外推得到的降水平移预报场插值到与观测降水场相同的网格上；步骤五：应用不带常数项的多元线性回归方法，估计出强迫项对应的x方向移动速度和y方向移动速度；步骤五中的回归方程为：（6）其中f
‑1表示平均强迫场，即：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（7）步骤六：利用拉格朗日向前积分方案对f
‑
0.5
进行外推，得到未来时刻t1的降水生消预报场rf1；外推公式为：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（8）由于外推得到的降水生消预报场所处格点不一定正好位于原网格点i，j上，所以利用kriging插值方法将外推得到的降水生消预报场插值到与观测降水场相同的网格上；步骤七：将降水平移预报场和降水生消预报场进行线性叠加，得到总的降水预报场r1。
6.作为本发明进一步改进，所述步骤六中的叠加公式为：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（9）由于降水量不可能为负，因此最后将所有小于零的r1赋值为零。
7.本发明提供了一种基于物理方程的降水外推预报方法，从动力学方程出发，在方程中直接包含降水的平移和生效的贡献，因此不仅考虑降水系统的水平移动、忽略新降水系统的生成和现有降水系统的消亡，还考虑降水的生效贡献，从而进行更合理有效的降水短临预报。
附图说明
8.图1是技术架构和执行流程图；图2是实施例原始降水观测场；图3是实施例平滑后的降水场；图4是实施例降水移动速度场；图5是实施例降水生效强迫场；图6是实施例降水平移预报场；图7是实施例降水强迫移动速度场；图8是实施例降水生效预报场；图9是实施例降水总预报场。
具体实施方式
9.下面结合附图与具体实施方式对本发明作进一步详细描述：本发明旨在提供一种基于物理方程的降水外推预报方法，从动力学方程出发，在方程中直接包含降水的平移和生效的贡献，最终发展更合理有效的降水短临预报技术。
10.本发明工作流程图如图1所示，具体为：第一步，获取当前时刻（t0）和其前两个时次（t
‑1、t
‑2）的降水量（r）数据。三者（r0、
r
‑1、r
‑2）均为二维空间格点场，三者的经向（x方向）和纬向（y方向）格点数分别相同，网格间距一致。
11.第二步，对r0、r
‑1、r
‑2分别进行二维平滑处理，去除孤立、偶发降水点，减少降水随机观测误差。平滑公式为： （10）其中z为包含n个样本一维空间序列。平滑时，先在x方向进行一维平滑，再在y方向进行一维平滑。
12.第三步，对物理要素的动力学方程进行差分化，对其中的参数在局部区域内利用统计方法求解。降水变化的原始动力学方程为：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（11）其中t为时间，f为外强迫，即r的个别变化率等于f。
13.将公式（11）方程等号的左边进行分解，可以进一步得到：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（12）其中u和v分别表示降水在x方向和y方向的平移速度，分别表示r在x方向和y方向的梯度，即将r的个别变化率分解为r的局地变化率和水平平流。
14.将公式（12）写在每一个网格点上，并且用差分代替微分，可以得到：（13）其中i和j分别为x方向和y方向的坐标，r
‑
1.5
、r
‑
0.5
分别表示t
‑1和t
‑2、t0和t
‑1的平均降水场，即：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（14）以每一个网格点（i，j）为中心，计算公式（13）中与r有关的各项，然后在9
×
9的区
域内应用带常数项的多元线性回归方法，客观估计得到；。
15.第四步，利用拉格朗日向前积分方案对r0进行外推，得到未来时刻（t1）的降水平移预报场（rs1）。外推公式为：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（15）其中：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（16）由于外推得到的降水平移预报场所处格点不一定正好位于原网格点（i，j）上，所以利用kriging插值方法将外推得到的降水平移预报场插值到与观测降水场相同的网格上。
16.第五步，考虑对降水生消的强迫不随时间发生变化，即：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（17）与公式（12）类似，将方程等号的左边进行分解可以得到：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（18）其中，为与公式（12）相区别，f在x方向和y方向的风速分别用uf和vf表示。
17.将公式（18）写在每一个网格点上，并且用差分代替微分，可以得到：（19）其中f
‑1表示平均强迫场，即：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（20）以每一个网格点（i，j）为中心，计算公式（19）中与f有关的各项，然后在9
×
9的区域内应用不带常数项的多元线性回归方法，客观估计得到和。
18.第六步，利用拉格朗日向前积分方案对f
‑
0.5
进行外推，得到未来时刻（t1）的降水生消预报场（rf1）。外推公式为：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（21）由于外推得到的降水生消预报场所处格点不一定正好位于原网格点（i，j）上，所以利用kriging插值方法将外推得到的降水生消预报场插值到与观测
降水场相同的网格上。
19.第七步，将降水平移预报场和降水生消预报场在每一个网格点上进行线性叠加，得到总的降水预报场（r1），即：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（22）由于降水量不可能为负，因此最后将所有小于零的r1赋值为零。
20.作为本发明具体实施例：对2021年4月30日15：00世界时的降水进行预报，预报区域：30
°
~40
°
n、115
°
~130
°
e。
21.第一步：获取当前时刻t0（2021年4月30日12：00世界时）和其前两个时次t
‑1（2021年4月30日13：00世界时）、t
‑2（2021年4月30日14：00世界时）的降水量（r）数据（r0、r
‑1、r
‑2；图2）；第二步：对r0、r
‑1、r
‑2分别进行二维平滑处理（图3）；第三步：应用带常数项的多元线性回归方法，估计出降水移动速度项（图4）和降水生消强迫项（图5）；第四步：利用拉格朗日向前积分方案对r0进行外推，得到未来时刻（t1）的降水平移预报场（rs1；图6）；第五步：应用不带常数项的多元线性回归方法，估计出强迫项对应的移动速度和（图7）；第六步：利用拉格朗日向前积分方案对f
‑
0.5
进行外推，得到未来时刻（t1）的降水生消预报场（rf1；图8）；第七步：将降水平移预报场和降水生消预报场进行线性叠加，得到总的降水预报场（图9）。
22.以上所述，仅是本发明的较佳实施个例而已，并非是对本发明作任何其他形式的限制，而是依据本发明的技术实质所作的任何修改或等同变化，仍然属于本发明所要求保护的范围。