CN110502849B - 一种应用于四维变分同化系统的扰动模式构建方法 - Google Patents

Abstract

本发明涉及气象数据处理和预测技术领域，提出一种应用于四维变分同化系统的扰动模式构建方法，包括以下步骤：从一组大气扰动方程出发，得到以扰动气压Π′ⁿ⁺¹为变量的亥姆霍兹方程，从n＝0时刻开始，依次求解亥姆霍兹方程得到n+1时刻的扰动变量集合x′ⁿ⁺¹；判断当前时刻是否到达观测时刻t，若否，则令求解得到的n+1时刻的扰动变量x′ⁿ⁺¹替代扰动变量x′ⁿ，然后继续求解其亥姆霍兹方程；若是，则将到达观测时刻t的扰动变量x′^t输入四维变分同化系统，计算其观测余差及极小化结果，判断极小化结果是否满足预设的计算精度，若否，则回到n＝0时刻对初始扰动变量集合x′⁰中各个变量进行调整，然后重新依次求解其亥姆霍兹方程；若是，则结束扰动模式积分过程。

Description

一种应用于四维变分同化系统的扰动模式构建方法

技术领域

本发明涉及气象数据处理和预测技术领域，更具体地，涉及一种应用于四维变分同化系统的扰动模式构建方法。

背景技术

现代天气预报是以数值天气预报为基础开展的，而数值天气预报则是根据当前气象观测数据，利用高性能计算机数值求解大气运动数学物理方程组，进而计算出未来指定日期内的天气。数值天气预报系统包括数值模式和资料同化系统，其中，资料同化系统针对各种观测数据进行资料分析处理，形成数值模式所必需的网格化初始大气要素场；数值模式是根据资料同化系统的初始场求解大气运动方程组的数学计算方法或程序。

GRAPES系统(Global/Regional Assimilation PrEdiction System)是中国气象局自主研发的新一代全球/区域通用数值天气预报系统。为提高区域GRAPES模式系统的预报准确率和区域观测资料的利用率需要构建可应用于四维变分同化系统的切线性模式。而为构建与GRAPES模式物理过程相互协调合理的切线性模式，需要发展基于GRAPES基本运动方程组的扰动模式。

发明内容

本发明为提高区域模式尤其是区域GRAPES数值模式的预报准确率和观测资料利用率，提供一种应用于四维变分同化系统的扰动模式的构建方法。

本发明的技术方案如下：

一种应用于四维变分同化系统的扰动模式构建方法，包括以下步骤：

S1：从GRAPES系统相关的一组大气扰动方程出发，由扰动纬向风速u′、扰动经向风速v′、扰动垂直风速w′、扰动比湿q′和扰动位温θ′组成扰动变量集合x′，并构造以n+1时刻扰动气压Π′ⁿ⁺¹为变量的亥姆霍兹方程；其中，扰动变量集合x′为包含所有待求解扰动变量的集合，即x′＝{u′,v′,w′,Π′,q′,θ′}；

S2：求解所述亥姆霍兹方程，由n时刻扰动气压Π′ⁿ得到n+1时刻的扰动气压Π′ⁿ⁺¹，再由n+1时刻的扰动气压Π′ⁿ⁺¹分别计算n+1时刻的扰动纬向风速u′ⁿ⁺¹、扰动经向风速v′ⁿ⁺¹、扰动垂直风速w′ⁿ⁺¹、扰动比湿q′ⁿ⁺¹、扰动位温θ′ⁿ⁺¹，得到n+1时刻的扰动变量集合x′ⁿ⁺¹；

S3：从n＝0时刻开始，对扰动模式进行连续积分：将S2步骤计算得到的n+1时刻的扰动变量集合x′ⁿ⁺¹更替为扰动变量集合x′ⁿ，并设置n＝n+1，然后跳转执行S2步骤，直到当前计算的时刻n为观测时刻t；

S4：将S3步骤中得到的扰动模式预报结果接入四维变分同化系统，计算其观测余差R，然后将所述观测余差R输入极小化过程，判断极小化过程的计算结果是否满足预设的计算精度，若否，则回到n＝0时刻对初始扰动变量集合x′⁰中各个变量进行调整，跳转执行S2步骤；若是，则将该初始扰动变量集合x′⁰作为最优分析增量进行输出，并结束扰动模式积分过程，完成扰动模式的构建。

本技术方案中，由于扰动变量集合x′ⁿ中需要包含能够完整反映大气变化的基本变量，即需要保证扰动模式在物理意义上具备完备性，因此本技术方案中的扰动变量集合x′ⁿ包括扰动纬向风速u′、扰动经向风速v′、扰动垂直风速w′、扰动气压Π′ⁿ⁺¹、扰动比湿q′和扰动位温θ′。在扰动模式构建过程中，通过从GRAPES系统相关的一组大气扰动方程出发，构造以n+1时刻扰动气压Π′ⁿ⁺¹为变量的亥姆霍兹方程，然后通过对上述亥姆霍兹方程求解，实现对扰动模式的连续积分，计算预报观测时刻t的扰动变量集合x′^t，再输入四维变分同化系统中进行计算精度预测，当计算精度满足需求，即完成扰动模式的构建。本技术方案通过在四维变分同化系统引入了扰动模式，替代现有技术中四维变分同化系统的切线性模式，在实际应用时能够节省现有的从非线性模式代码到切线性代码再到伴随模式代码的复杂编码过程，同时使数值模式系统和资料同化系统之间的物理联系更加紧密。

优选地，S1步骤中，亥姆霍兹方程中以地形追随高度坐标

代替自然高度坐标z，其中，/>

λ为纬度，/>

为经度，z为自然高度。

优选地，S1步骤中的具体步骤如下：

S1.1：根据地形追随高度坐标下的球面GRAPES非线性模式大气原始方程组，构造扰动变量集合x′的扰动预报方程组，其中，扰动预报方程组的具体公式如下：

其中，

/>

其中，t表示时间，x表示纬向坐标，y表示经向坐标；A_x、A_y、

表示变量A＝u′,v′,w′,Π′,q′,θ′的空间梯度分量；C_p表示定压比热，f为科氏参数，L表示凝结潜热，P_w′为扰动降水率，R_d为干气体常数，/>

为/>

坐标的垂直速度；Z_sx、Z_sy、/>

分别表示模式面空间梯度的三个分量；L_πx、L_πy、/>

分别表示地形追随坐标下气压梯度的三个分量；D₃为散度；/>

为干空气热力常数；

S1.2：采用两个时间层半隐式半拉格朗日时间差分方案对所述扰动预报方程组进行时间离散差分计算，并通过调节合适的隐式权重α，得到最佳预报效果，其中所述两个时间层半隐式半拉格朗日时间差分方案对应的具体公式如下：

其中，a表示到达点，d表示上游点，Δt为时间步长；A′代表扰动变量，且A′＝u′,v′,w′,Π′,q′,θ′；L_A代表每个扰动方程右端所有各项的总和；

表示隐式计算，/>

表示显式计算；α为隐式权重，β为显式权重，且满足α+β＝1；

S1.3：上游点d在三维空间中的位置由GRAPES系统预报的基本态风场(u,v,w)通过计算大气质点的后向轨迹来确定，其中u为纬向风速，v为经向风速，w为垂直风速；

S1.4：在上游点计算n时刻的已知项，把n时刻的扰动平流项和位温扰动项均写在网格点上；本步骤中，与快波(重力波、声波)有关的气压梯度力、科氏力和散度项半隐式计算，与平流有关的慢过程全显式计算，故此把扰动平流项和位温扰动项均在n时刻网格点上计算；

S1.5：对所述扰动预报方程组通过时间和空间差分离散化后，再通过代数消元得到只含变量Π′的亥姆霍兹方程，其具体公式如下：

B₁Π_i,j,k+B₂Π_i-1,j,k+B₃Π_i+1,j,k+B₄Π_i,j-1,k+B₅Π_i,j+1,k+B₆Π_i+1,j+1,k+B₇Π_i+1,j-1,k+B₈Π_i-1,j-1,k+B₉Π_i-1,j+1,k+B₁₀Π_i,j,k-1+B₁₁Π_i-1,j,k-1+B₁₂Π_i+1,j,k-1+B₁₃Π_i,j-1,k-1+B₁₄Π_i,j+1,k-1+B₁₅Π_i,j,k+1+B₁₆Π_i-1,j,k+1+B₁₇Π_i+1,j,k+1+B₁₈Π_i,j-1,k+1+B₁₉Π_i,j+1,k+1＝B₀

其中，B₀～B₁₉为系数，B₁～B₁₉由GRAPES系统预报的基本态决定，B₀与n时刻扰动态在上游点的取值有关；Π_i,j,k表示GRAPES系统网格点上的扰动气压，i,j,k表示GRAPES系统中网格点的三维序号；

S1.6：采用广义共轭余差法对所述亥姆霍兹方程进行迭代求解，得到下一时刻的扰动气压Π′ⁿ⁺¹，然后根据所述扰动气压Π′ⁿ⁺¹计算得到n+1时刻的其他五个扰动变量u′ⁿ⁺¹,v′ⁿ⁺¹,w′ⁿ⁺¹,q′ⁿ⁺¹,θ′ⁿ⁺¹。

本优选方案中，在GRAPES系统网格点上计算到达点a，通过拉格朗日后向轨迹计算得到上游点d，且一般默认到达点a在n+1时刻，上游点d在n时刻。

优选地，S1.2步骤中，对扰动预报方程组右端项L_A进行计算时，其空间离散差分计算中垂直方向采用Charney-Phillips跳点，水平方向采用ArakawaC网格；进行垂直差分离散计算时，采用非均匀分层的二阶精度差分方案进行计算，能够有效提高计算精度。

优选地，S1.2步骤中的隐式权重α取值为0.55，能够确保稳定积分。

优选地，S1.3步骤中上游点d位置上变量的取值通过双线性插值法由三维网格点计算。

优选地，S2步骤中，求解所述亥姆霍兹方程的过程中，采用有预条件的广义共轭余差法进行迭代求解，其求解计算精度设置为10^-32。

优选地，S4步骤中，对初始扰动变量集合x'⁰中各个变量进行调整时，对其极小化计算过程采用有限内存BFGS方法进行迭代求解。

优选地，S4步骤中，观测余差R的计算公式如下：

其中，x′^t为观测时刻t对应的扰动变量集合，x^t为GRAPES非线性模式积分到观测时刻t的基本态变量集合，y^t表示观测时刻t的观测量；H(·)为非线性观测算子，

为线性观测算子。

与现有技术相比，本发明技术方案的有益效果是：通过在GRAPES基本原始方程组的基础上引入了扰动模式，替代现有技术中四维变分同化系统的切线性模式，能够提高扰动变量的线性预报，有效强化了数值模式系统和资料同化系统之间的物理联系；能够简化现有的从非线性模式代码到切线性代码再到伴随模式代码的复杂编码过程，有效提高代码编写的效率；有效提高数值天气预报的精度，降低同化系统的计算复杂度。

附图说明

图1为本实施例的应用于四维变分同化系统的扰动模式构建方法的流程图。

图2为本实施例进行11步积分后第20层的水平扰动风场和由非线性模式预测的扰动风场对比图。

其中，图2(a)为积分两组非线性模式的差值场示意图；图2(b)为从扰动初值出发经过扰动模式积分11步后得到的结果示意图。

图3为本实施例的预报值与非线性模式预测的扰动风场的对比图。

图4为本实施例进行11步积分后扰动量的垂直剖面图和非线性模式预测的扰动量垂直剖面图的对比图。

具体实施方式

附图仅用于示例性说明，不能理解为对本专利的限制；

对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。

下面结合附图和实施例对本发明的技术方案做进一步的说明。

如图1所示，为本实施例的应用于四维变分同化系统的扰动模式构建方法的流程图。

本实施例提出一种应用于四维变分同化系统的扰动模式构建方法，包括以下步骤：

S1：从GRAPES系统相关的一组大气扰动方程出发，由扰动纬向风速u′、扰动经向风速v′、扰动垂直风速w′、扰动气压Π′ⁿ⁺¹、扰动比湿q′和扰动位温θ′组成扰动变量集合x′，并构造以n+1时刻扰动气压Π′ⁿ⁺¹为变量的亥姆霍兹方程；其中，扰动变量集合x′为包含所有待求解扰动变量的集合，即x′＝{u′,v′,w′,Π′,q′,θ′}；其具体步骤如下：

S1.1：根据地形追随高度坐标下的球面GRAPES非线性模式大气原始方程组构造扰动变量集合x′的预报方程组，其中，扰动预报方程组的具体公式如下：

其中，

其中，t表示时间，x表示纬向坐标，y表示经向坐标；A_x、A_y、

表示变量A＝u′,v′,w′,Π′,q′,θ′的空间梯度分量；C_p表示定压比热，f为科氏参数，L表示凝结潜热，P_w′为扰动降水率，R_d为干气体常数，/>

为/>

坐标的垂直速度；Z_sx、Z_sy、/>

分别表示模式面空间梯度的三个分量；L_πx、L_πy、/>

分别表示地形追随坐标下气压梯度的三个分量；D₃为散度；/>

为干空气热力常数；

S1.2：采用两个时间层半隐式半拉格朗日时间差分方案对所述扰动预报方程组进行时间离散差分计算，并通过调节合适的隐式权重α得到最佳预报效果，其中，两个时间层半隐式半拉格朗日时间差分方案对应的具体公式如下：

其中，a表示到达点，d表示上游点，本实施例中，到达点a即为n+1时刻，上游点d即为n时刻；Δt为时间步长；A′代表扰动变量，且A′＝u′,v′,w′,Π′,q′,θ′；L_A代表每个扰动方程右端所有各项的总和；

表示隐式计算，/>

表示显式计算；α为隐式权重，β为显式权重，且满足α+β＝1；本实施例中，α＝0.55，β＝0.45；/>

其中，在对扰动预报方程组右端项L_A进行空间差分计算时，垂直方向采用Charney-Phillips跳点，水平方向采用ArakawaC网格；对扰动预报方程组进行垂直差分离散计算时，采用非均匀分层的二阶精度差分方案进行计算，用于提高其计算精度；

S1.3：上游点d在三维空间中的位置由GRAPES系统预报的基本态风场(u,v,w)通过计算大气质点的后向轨迹来确定，其中u为纬向风速，v为经向风速，w为垂直风速，且本步骤中，上游点位置上变量的取值通过双线性插值法由三维网格点计算；

S1.4：在上游点计算n时刻的已知项，把n时刻的扰动平流项和位温扰动项均写在网格点上；本实施例中，与快波(重力波、声波)有关的气压梯度力、科氏力和散度项半隐式计算，与平流有关的慢过程全显式计算，故此把扰动平流项和位温扰动项均在n时刻网格点上计算，例如，对于

方程的右端项中的平流项和位温扰动项均写在n时刻网格点上，得到以下公式X′_u：

其中，

为上游点处的扰动纬向风速；β_c＝βΔt，β为显式权重，且β＝0.45；L_θ＝C_pθ，C_p为定压比热，θ是基本态位温，θ′ⁿ为扰动位温，Δt为时间步长，/>

表示扰动风场对纬向基本态风场的平流，变量上标n表示n时刻。

S1.5：对所述扰动预报方程组通过时间和空间差分离散化后，再通过代数消元得到只含变量Π′的亥姆霍兹方程，其具体公式如下：

B₁Π_i,j,k+B₂Π_i-1,j,k+B₃Π_i+1,j,k+B₄Π_i,j-1,k+B₅Π_i,j+1,k+B₆Π_i+1,j+1,k+B₇Π_i+1,j-1,k+B₈Π_i-1,j-1,k+B₉Π_i-1,j+1,k+B₁₀Π_i,j,k-1+B₁₁Π_i-1,j,k-1+B₁₂Π_i+1,j,k-1+B₁₃Π_i,j-1,k-1+B₁₄Π_i,j+1,k-1+B₁₅Π_i,j,k+1+B₁₆Π_i-1,j,k+1+B₁₇Π_i+1,j,k+1+B₁₈Π_i,j-1,k+1+B₁₉Π_i,j+1,k+1＝B₀

其中，B₀～B₁₉为系数，B₁～B₁₉由GRAPES系统预报的基本态决定，B₀与n时刻扰动态在上游点的取值有关；Π_i,j,k表示GRAPES系统网格点上的扰动气压，i,j,k表示GRAPES系统中网格点的三维序号；

S1.6：采用广义共轭余差法对所述亥姆霍兹方程进行迭代求解，得到下一时刻的扰动气压Π′ⁿ⁺¹，然后根据所述扰动气压Π′ⁿ⁺¹计算得到n+1时刻的其他五个扰动变量u′ⁿ⁺¹,v′ⁿ⁺¹,w′ⁿ⁺¹,q′ⁿ⁺¹,θ′ⁿ⁺¹。

本步骤中，根据地形追随高度坐标下的球面GRAPES非线性模式大气原始方程组构造初始扰动变量x′ⁿ的扰动预报方程组，具体地，将每个初始扰动变量x分解为基本态与扰动变量x′ⁿ，然后经过线性化处理得到扰动变量x′ⁿ的预报方程组，然后采用两个时间层半隐式半拉格朗日时间差分方案对所述扰动预报方程组进行时间和空间离散差分计算以及代数消元，得到只含变量Π′的亥姆霍兹方程，应用于扰动模式的连续积分中。

S2：求解所述亥姆霍兹方程，由n时刻扰动气压Π′ⁿ得到n+1时刻的扰动气压Π′ⁿ⁺¹，再由n+1时刻的扰动气压Π′ⁿ⁺¹分别计算n+1时刻的扰动纬向风速u′ⁿ⁺¹、扰动经向风速v′ⁿ⁺¹、扰动垂直风速w′ⁿ⁺¹、扰动比湿q′ⁿ⁺¹、扰动位温θ′ⁿ⁺¹，得到n+1时刻的扰动变量集合x′ⁿ⁺¹。

本步骤在求解所述亥姆霍兹方程的过程中，采用有预条件的广义共轭余差法进行迭代求解，其求解计算精度设置为10^-32

S3：对扰动模式进行连续积分：将S2步骤计算得到的n+1时刻的扰动变量集合x′ⁿ⁺¹更替为扰动变量集合x^′n，并设置n＝n+1，然后跳转执行S2步骤，直到当前计算的时刻n为观测时刻t；即从最初扰动开始，按照时间先后依次连续执行S2。

S4：将S3步骤中得到的扰动模式预报结果接入四维变分同化系统，计算其观测余差，然后将所述观测余差输入极小化过程，判断极小化过程的计算结果是否满足预设的计算精度，若否，则回到n＝0时刻，对初始扰动变量集合x′⁰中各个变量进行调整，跳转执行S2步骤；若是，则将该初始扰动变量集合x′⁰作为最优分析增量进行输出，并结束扰动模式积分过程，完成扰动模式的构建。

本步骤中，从n时刻扰动x′ⁿ出发，求解其亥姆霍兹方程，得到n+1时刻的扰动x′ⁿ⁺¹，然后以n+1时刻的扰动量作为初值，重复上述步骤，得到目标观测时刻的扰动量，再输入资料同化系统中进行分析，得到其分析增量，若达到预设的计算精度则完成数值预测，若还未达到预设的计算精度要求，则需要对初始扰动变量进行调整，再重复上述步骤。

本步骤中，观测余差R的计算公式如下：

其中，x′^t为观测时刻t对应的扰动变量集合，x^t为GRAPES非线性模式积分到观测时刻t的基本态变量集合，y^t表示观测时刻t的观测量；H(·)为非线性观测算子，

为线性观测算子。

本步骤中，在对初始扰动变量集合x′⁰中各个变量进行调整时，对其极小化计算过程采用有限内存BFGS方法进行迭代求解。

在具体实施过程中，为了测试本实施例的应用于四维变分同化系统的扰动模式构建方法，采用Frobonius范数对本实施例的合理性进行检测。首先利用非线性模式产生检验所需要的扰动场，具体地，将非线性模式GRAPES的物理过程全部关闭(干模式)，设计两组试验，第一组试验直接以全球模式分析场x₀作为非线性模式初始场进行时间积分，第二组试验在全球分析场上叠加一个小扰动δx₀后进行时间积分，然后将这两组试验的差值用于对扰动模式积分结果进行检验。采用Frobonius范数进行正确性检验，将第一组试验结果表示为NLM(x₀)，将第二组试验结果表示为NLM(x₀+γδx₀)，其中γ为正数，用于控制初始扰动的大小。

假设R_0→iγδx₀表示为非线性模式的切线性模式从初值γδx₀出发的一次积分，得到以下公式：

NLM(x₀+γδx₀)＝NLM(x₀)+R_0→iγδx₀+O(γ²)

其中，i表示积分步数，O(γ²)表示以γ为变量的泰勒展开的二阶项。

假设PFM(γδx₀)为初值γδx₀出发的一次扰动模式积分，其中本实施例的扰动模式表示为PFM(γδx₀)＝M_0→iγδx₀，则得到如下Frobonius范数的计算公式：

其中，当γ参数趋近于0时，如果扰动模式的积分结果与非线性模式的余差接近，则F(γ)应趋近1，进而说明扰动模式M_0→i与切线性模式R_0→i相当。

在具体实施过程中，采用2015年5月16日00时ECMWF分析场作为初值进行积分。作为对照的非线性模式GRAPES与本实施例的扰动模式GRAPE_PF的水平分辨率均取12km，垂直层为67层，模式顶高度为30km，时间积分步长为60s。模式范围为96-110.04°E，16-29.8°N。利用三维变分系统3Dvar中的单站探空分析功能，在模式面约6km高度，约在模式25层的东经102.7°北纬22.78°位置上进行气压要素的单点分析，得到水平方向尺度半径约120km符合高斯分布特征的初始高压扰动Π′。

如图2所示，为本实施例进行11步积分后第20层的水平扰动风场和由非线性模式预报的扰动风场对比图。其中，图2(a)表示积分两组非线性模式的差值场，代表扰动的实际演变；图2(b)表示从扰动初值出发经过扰动模式积分11步后得到的结果，代表扰动的预报值。

如图3所示，为本实施例的预报值与非线性模式预测值的对比图。

由图可知，本实施例能够很好地模拟出扰动风场的变化。

如图4所示，为本实施例进行11步积分后扰动量的垂直剖面图和非线性模式预测的扰动量垂直剖面图的对比图。图中，本实施例得到的各个扰动变量的数量级和正负数值分布与非线性模式得到的各个扰动变量的数量级和正负数值分布高度一致，说明本实施例的扰动模式所包含的物理方案和计算结果是正确的，即由于初始高压扰动激发的重力波涟漪，引起了垂直方向的运动增量，并进而引起位温和湿度的垂直重新分配。

附图中描述位置关系的用语仅用于示例性说明，不能理解为对本专利的限制；

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (7)

1.一种应用于四维变分同化系统的扰动模式构建方法，其特征在于，包括以下步骤：

S1：从GRAPES系统相关的一组大气扰动方程出发，由扰动纬向风速u′、扰动经向风速v′、扰动垂直风速w′、扰动气压Π′ⁿ⁺¹、扰动比湿q′和扰动位温θ′组成扰动变量集合x′，并构造以n+1时刻扰动气压Π′ⁿ⁺¹为变量的亥姆霍兹方程；其中，扰动变量集合x′为包含所有待求解扰动变量的集合，即x′＝{u′,v′,w′,Π′,q′,θ′}；所述亥姆霍兹方程中以地形追随高度坐标

代替自然高度坐标z，其中，/>

λ为纬度，/>

为经度，z为自然高度；

其中，构造以n+1时刻扰动气压Π′ⁿ⁺¹为变量的亥姆霍兹方程的步骤包括：

S1.1：根据球面GRAPES非线性模式大气原始方程组，构造扰动变量集合x′的扰动预报方程组；所述扰动预报方程组的具体公式如下：

其中，

其中，t表示时间，x表示纬向坐标，y表示经向坐标；A_x、A_y、

表示变量A＝u′,v′,w′,Π′,q′,θ′的空间梯度分量；C_p表示定压比热，f为科氏参数，L表示凝结潜热，P_w′为扰动降水率，R_d为干气体常数，/>

为/>

坐标的垂直速度；Z_sx、Z_sy、/>

分别表示模式面空间梯度的三个分量；L_πx、L_πy、/>

分别表示地形追随坐标下气压梯度的三个分量；D₃为散度；/>

为干空气热力常数；

S1.2：采用两个时间层半隐式半拉格朗日时间差分方案对所述扰动预报方程组进行时间差分计算，并通过调节合适的隐式权重α，得到最佳预报效果，其中所述两个时间层半隐式半拉格朗日时间差分方案对应的具体公式如下：

其中，a表示到达点，d表示上游点，Δt为时间步长；A′代表扰动变量，且A′＝u′,v′,w′,Π′,q′,θ′；L_A代表每个扰动方程右端所有各项的总和；

表示隐式计算，/>

表示显式计算，α为隐式权重，β为显式权重，且满足α+β＝1；

S1.3：上游点d在三维空间中的位置由GRAPES系统预报的基本态风场(u,v,w)通过计算大气质点的后向轨迹来确定，其中u为纬向风速，v为经向风速，w为垂直风速；

S1.4：在上游点计算n时刻的已知项，把n时刻的扰动平流项和位温扰动项均写在网格点上；

S1.5：对所述扰动预报方程组通过时间和空间差分离散化后，再通过代数消元得到只含变量Π′的亥姆霍兹方程，其具体公式如下：

B₁Π_i,j,k+B₂Π_i-1,j,k+B₃Π_i+1,j,k+B₄Π_i,j-1,k+B₅Π_i,j+1,k+B₆Π_i+1,j+1,k+B₇Π_i+1,j-1,k+B₈Π_i-1,j-1,k+B₉Π_i-1,j+1,k+B₁₀Π_i,j,k-1+B₁₁Π_i-1,j,k-1+B₁₂Π_i+1,j,k-1+B₁₃Π_i,j-1,k-1+B₁₄Π_i,j+1,k-1+B₁₅Π_i,j,k+1+B₁₆Π_i-1,j,k+1+B₁₇Π_i+1,j,k+1+B₁₈Π_i,j-1,k+1+B₁₉Π_i,j+1,k+1＝B₀

其中，B₀～B₁₉为系数，B₁～B₁₉由GRAPES系统预报的基本态决定，B₀与n时刻扰动态在上游点的取值有关；Π_i,j,k表示GRAPES系统所定义的网格点上的扰动气压；i,j,k表示GRAPES系统中网格点的三维序号；

S1.6：采用广义共轭余差法对所述亥姆霍兹方程进行迭代求解，得到下一时刻的扰动气压Π′ⁿ⁺¹，然后根据所述扰动气压Π′ⁿ⁺¹计算得到n+1时刻的其他五个扰动变量u′ⁿ⁺¹,v′ⁿ⁺¹,w′ⁿ⁺¹,q′ⁿ⁺¹,θ′ⁿ⁺¹；

S2：求解所述亥姆霍兹方程，由n时刻扰动气压Π′ⁿ得到n+1时刻的扰动气压Π′ⁿ⁺¹，再由n+1时刻的扰动气压Π′ⁿ⁺¹分别计算n+1时刻的扰动纬向风速u′ⁿ⁺¹、扰动经向风速v′ⁿ⁺¹、扰动垂直风速w′ⁿ⁺¹、扰动比湿q′ⁿ⁺¹、扰动位温θ′ⁿ⁺¹，得到n+1时刻的扰动变量集合x′ⁿ⁺¹；

S3：从n＝0时刻开始，对扰动模式进行连续积分：将S2步骤计算得到的n+1时刻的扰动变量集合x′ⁿ⁺¹更替为扰动变量集合x′ⁿ，并设置n＝n+1，然后跳转执行S2步骤，直到当前计算的时刻n为观测时刻t；

S4：将S3步骤中得到的扰动模式积分结果接入四维变分同化系统，计算其观测余差R，然后将所述观测余差R输入极小化过程，判断极小化过程的计算结果是否满足预设的计算精度，若否，则回到n＝0时刻对初始扰动变量集合x′⁰中各个变量进行调整，跳转执行S2步骤；若是，则将该初始扰动变量集合x′⁰作为最优分析增量进行输出，并结束扰动模式积分过程，完成扰动模式的构建。

2.根据权利要求1所述的扰动模式构建方法，其特征在于：所述S1.2步骤中，对所述扰动预报方程组右端项L_A进行计算时，其空间离散差分计算中垂直方向采用Charney-Phillips跳点，水平方向采用Arakawa C网格；进行垂直差分离散计算时，采用非均匀分层的二阶精度差分方案进行计算。

3.根据权利要求1所述的扰动模式构建方法，其特征在于：所述S1.2步骤中的隐式权重α取值为0.55。

4.根据权利要求1所述的扰动模式构建方法，其特征在于：所述S1.3步骤中，上游点d位置上变量的取值通过双线性插值法由三维网格点计算。

5.根据权利要求1所述的扰动模式构建方法，其特征在于：所述S2步骤中，求解所述亥姆霍兹方程的过程中，采用有预条件的广义共轭余差法进行迭代求解，其求解计算精度设置为10^-32。

6.根据权利要求1所述的扰动模式构建方法，其特征在于：所述S4步骤中，对初始扰动变量集合x′⁰中各个变量进行调整时，对其极小化计算过程采用有限内存BFGS方法进行迭代求解。

7.根据权利要求1所述的扰动模式构建方法，其特征在于：所述S4步骤中，观测余差R的计算公式如下：

其中，x′^t为观测时刻t对应的扰动变量集合，x^t为GRAPES非线性模式积分到观测时刻t的基本态变量集合，y^t表示观测时刻t的观测量；H(·)为非线性观测算子，

为线性观测算子。/>

Images