一种基于WRF模式的台风动力平衡约束下的变分同化方法

摘要

本发明涉及一种基于WRF模式的台风动力平衡约束下的变分同化方法，包含如下步骤：i、根据WRF模式的垂直σ坐标下大气原始运动方程得到动力平衡方程；ii、利用弱约束条件法，将动力平衡方程表示成弱约束罚项J

说明书

技术领域

本发明涉及台风在动力约束条件下的变分同化方法，尤其是一种基于WRF模式的台风 动力平衡约束下的变分同化方法，属于大气科学研究领域。

背景技术

台风是热带地区的强灾害性天气系统，每年都给沿海地区带来巨大损失。因此，提高台 风预报的准确性，关乎生命和社会安全。数值模式作为当今天气预报的重要手段，能够对台 风路径和强度进行客观预报。由于大气数值预报属于初值问题，不仅要求模式能够精确的反 映天气过程的演变规律，还要求初始时刻的大气状态是精确的。

对于改善模式初始场，资料同化是当今的主流方法，它的作用是将预报场和多源观测资 料以最优化的方式相结合，从而得到对数值模式初始状态具有最优估计的分析场。在所有同 化方法中，三维变分(3DVAR)方法(Parrish and Derber1992)具备节约计算资源的重要优 势，在许多国家的业务预报中被广泛应用。在台风数值预报的应用和研究方面，3DVAR方法 也取得了许多重要成果，包括同化常规资料(Wang et al.2011)、飞机观测(chou et al.2011)、 人造涡旋(Xiao et al.2000)、卫星资料(Zou and Xiao et al.2000)和雷达资料(Zhao et al.2009) 等。然而由于缺少数值模式的动力约束，3DVAR方法的缺陷在于同化过程中非直接观测量的 调整往往不够合理，容易引起同化分析场中大气各状态量的不平衡，比如气压场和风场的结 构配置不满足物理平衡关系。因此，在现有的3DVAR同化理论基础上，加入合适的动力或 物理约束，用来改善分析场的协调性，是提高台风初始场准确程度的有效途径。

变分同化方法基于最大似然估计理论和最小二乘法，通过对所定义的目标泛函(代价函 数)求极值来估计最优分析场。而约束条件下的代价函数极值问题与简单多元函数的极值问 题类似，Sasaki(1970a，b)提出了2种方法来处理约束条件下的变分同化问题，即强约束法 和弱约束法。此后的一些研究探讨了弱约束条件下的变分同化问题。Gao et al.(1999)将连续 方程加入ARPS模式的3DVAR同化系统，用来约束水平风场的连续性；Liang et al.(2006)基 于MM5模式的3DVAR同化系统，将大气状态量的时间倾向作为弱约束，提出模式约束的 3DVAR同化方案，旨在让分析场尽可能的满足模式平衡；Ge et al.(2012)将气压诊断约束方 程引入ARPS-3DVAR同化系统，用来模拟对流和风暴系统。这些研究都表明加入适合的物理 或动力约束，对改善资料同化效果具有正面作用。

考虑到对于台风预报而言，初始状态的风场和气压场平衡很重要，因此迫切需要研究适 用于台风同化的约束平衡条件，并将其应用于变分同化中。在同化各类观测资料时，能够有 效保证台风分析场的平衡性，从而提高台风数值预报能力。

发明内容

本发明针对现有技术的不足，基于气象领域广泛应用的WRF(Weather Research and  Forecasting)中尺度天气模式的变分同化系统WRFDA(WRF data assimilation)，提出了适用 于台风的动力约束变分同化方案，将表征台风动力平衡的约束方程加入到原始的WRFDA变 分同化系统中，旨在改善各类观测资料在台风同化中的效果，提高台风预报能力。

为实现以上的技术目的，本发明将采取以下的技术方案：

①将WRF模式的垂直σ坐标下大气原始运动方程进行定常和无粘处理，同时对原始运动方 程中的垂直平流项忽略，得到非线性动力平衡约束矢量方程

②将步骤①中的矢量方程采用弱约束方法，得到变分同化中的动力约束罚项 $J_c=\frac{1}{2}G{\left(x\right)}^T{\Gamma }^{-1}G\left(x\right),$所述动力权重Г^-1取值为10³；

③将步骤②中的动力约束罚项J_c加入WRFDA原始的变分代价函数中，得到弱约束条件下的 变分同化代价函数J＝J_o+β₁J_b+β₂J_e+J_c；

④将上述步骤③中的弱约束条件下的变分同化代价函数采用共轭梯度法求解代价函数的极小 值，得到同化分析场。

上述方案的具体方案如下：

台风动力平衡约束下的WRF变分同化方案，具体包括以下步骤：i、根据控制大气运动 的基本方程推导适用于台风的动力平衡方程；ii、根据变分同化原理，并采用Sasaki(1970a， b)提出的弱约束条件法，将动力平衡方程表示成变分同化代价函数中的弱约束罚项J_e，并加 入WRFDA系统的代价函数中，得到新的变分代价函数；iii、采用量纲分析法估算约束罚项J_e在代价函数中的权重(即动力权重)大小；iv、根据变分同化原理，在原始的WRFDA同化 系统的基础上，开发与动力平衡方程相对应的切线性算子和伴随算子，计算约束罚项J_e的梯 度，得到新的代价函数梯度，完成改进的WRFDA同化系统。

所述步骤i中，根据WRF模式的垂直σ坐标下的大气水平运动方程来推导表征台风动力 特征的平衡方程。σ坐标下的大气原始运动方程为，

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+\omega \frac{\partial u}{\partial \sigma }-\mathrm{fv}+\frac{\partial \phi }{\partial x}+\frac{1}{\rho }\frac{\partial p}{\partial x}+F_x=0$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+\omega \frac{\partial v}{\partial \sigma }+\mathrm{fu}+\frac{\partial \phi }{\partial y}+\frac{1}{\rho }\frac{\partial p}{\partial y}+F_y=0$

对原始运动方程做定常和无粘近似，同时不考虑垂直平流项的作用，简化后的大气运动方程 为，

$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}-\mathrm{fv}+\frac{\partial \phi }{\partial x}+\frac{1}{\rho }\frac{\partial p}{\partial x}=0$

$u=\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+\mathrm{fu}+\frac{\partial \phi }{\partial y}+\frac{1}{\rho }\frac{\partial p}{\partial y}=0$

将上式写成矢量形式，得到非线性动力平衡约束方程，

上式右端第一项为非线性平流项，第二项为科氏力项，第三和第四项为气压梯度项。其中，是水平风速矢量(风分量u和v)，f是科氏参数，φ是位势高度，ρ是大气密度，p是气压。 所述步骤ii中，采用Sasaki(1970a，b)提出的弱约束法，变分同化中的动力约束罚项表示为，

$J_c=\frac{1}{2}G{\left(x\right)}^T{\Gamma }^{-1}G\left(x\right)$

其中x表示大气各状态量，G(x)代表表征台风动力特征的平衡方程；而Γ^-1为罚值系数(或称 动力权重)，其数值的大小决定了约束条件在同化中所起的作用。

WRFDA原始的变分代价函数为，

J＝J_o+β₁J_b+β₂J_e

其中，J_o表示与观测相关的惩罚项，J_b和J_e表示与背景场相关的惩罚项。J_b与静止背景误差 协方差有关，J_e与集合背景误差协方差有关。

在此基础上，将动力平衡方程G(x)以弱约束罚项J_c的形式加入代价函数，实现弱约束条 件下的变分同化方案。新的代价函数为，

J＝J_o+β₁J_b+β₂J_e+J_c

所述步骤iii中，为了让观测罚项和约束罚项在最小化迭代过程中起到相当的作用，通过 量纲分析法来估计动力权重Г^-1的大小。在动力方程G(x)中，平流项和气压梯度项为较重要 的贡献项，由此二项的量纲分析来确定方程G(x)量纲。

所述步骤iv中，由于WRFDA变分同化采用共轭梯度法(Navon and Legler1987)求解代 价函数的极小值，在最小化过程的每一次迭代步都需要计算代价函数的梯度，

$▿J=▿J_o+{\beta }_1▿J_b+{\beta }_2▿J_e+▿J_c$

其中约束罚项J_c的梯度表示为，

$▿J_c=G^T{\Gamma }^{-1}[{\mathrm{Gx}}^{\prime }+G\left(x_b\right)]$

其中x′表示大气各状态量的切线性增量，x_b表示背景场，G为动力平衡方程G的切线性算子， G^T为其对应的伴随算子。

切线性算子G为，

其中上横线表示背景态，上标表示增量。

通过上述计算约束方程的切线性算子和伴随算子G^T，求得约束项梯度与WRFDA变分 同化系统现有的观测罚项和背景罚项的梯度相加，求得新的代价函数梯度采用共轭梯度 法求解代价函数的极小值，得到同化分析场。

根据以上的技术方案，可以实现以下的有益效果：

1、本发明基于弱约束条件下的变分同化理论，在WRFDA同化系统的框架上，开发了台 风动力平衡约束下的变分同化方法。与现有技术相比，这种方法可以在台风资料同化中有效 的平衡分析场中的气压和风场结构，满足台风数值预报对初始条件各状态量的平衡性和准确 性要求。

2、本发明所述方法，充分考虑了台风的动力平衡过程，不同于简单的梯度风平衡关系， 除离心力项外还包含了其它非线性平流过程的作用，更接近完整的大气运动方程，提供了能 够正确反映台风质量场和动量场平衡的同化方法。

3、本发明所述方法相比现有的技术，在改善台风同化效果的基础上，有助于台风强度预 报水平的提高。

附图说明

图1是改进的变分同化方案流程图；

图2是改进的变分同化方案迭代收敛性检验图；

图3是WRFDA原始同化方案与改进的同化方案在2010年7月21日18时同化雷达 T-TREC风场资料前后的台风“灿都”结果(风场和气压场)比较图，其中同化前背景场(a， d)、原始方案同化分析场(b，e)以及改进方案同化分析场(c，f)，图中黑色圆点表示最佳 路径观测的台风中心位置；

图4是基于WRFDA原始同化方案与改进的同化方案对台风强度预报的结果比较图。

具体实施方式

1.台风动力平衡方程的推导

由于WRF模式采用地形跟随σ坐标，本发明根据σ坐标下的大气水平运动方程来推导表征 台风动力特征的平衡方程。σ坐标下的大气水平运动方程为，

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+\omega \frac{\partial u}{\partial \sigma }-\mathrm{fv}+\frac{\partial \phi }{\partial x}+\frac{1}{\rho }\frac{\partial p}{\partial x}+F_x=0$

$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+\omega \frac{\partial v}{\partial \sigma }+\mathrm{fu}+\frac{\partial \phi }{\partial y}+\frac{1}{\rho }\frac{\partial p}{\partial y}+F_y=0$

对原始大气运动方程做定常和无粘近似，同时不考虑垂直平流项的作用，简化后的大气运动 方程为，

$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}-\mathrm{fv}+\frac{\partial \phi }{\partial x}+\frac{1}{\rho }\frac{\partial p}{\partial x}=0$

$u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+\mathrm{fu}+\frac{\partial \phi }{\partial y}+\frac{1}{\rho }\frac{\partial p}{\partial y}=0$

将上式写成矢量形式，得到非线性动力平衡约束方程G，

其中，是水平风速矢量(风分量u和v)，f是科氏参数，φ是位势高度，ρ是大气密度，p 是气压。

非线性动力平衡方程G的第一项为非线性平流项，第二项为科氏力项，第三和第四项为气压 梯度项，其本质上包含了表征台风动力特征的梯度风平衡关系。将非线性平流项做矢量变换，

并且将其转换至柱坐标系只考虑切向风v_T，

上式表明非线性平流项中包含了梯度风平衡关系中的离心力项因此，非线性动力平衡约 束方程G实际上体现了比梯度风平衡方程更为复杂的动力平衡关系。

2.动力平衡约束下的变分代价函数

变分同化原理是通过对所定义的目标泛函(或称代价函数)求极值来求得最优分析场。 原始的WRFDA变分代价函数为，

J＝J_o+β₁J_b+β₂J_e

其中，J_o表示与观测相关的惩罚项，J_b和J_e表示与背景场相关的惩罚项，J_b与静止背景误差 协方差有关，J_e与集合背景误差协方差有关。

在此基础上，将动力平衡方程以弱约束罚项J_c的形式加入代价函数，用来实现弱约束条 件下的变分同化方案，约束分析场中的风场和气压场。新的代价函数为，

J＝J_o+β₁J_b+β₂J_e+J_c

上式右端第四项J_c为动力约束罚项，表示为，

$J_c=\frac{1}{2}G{\left(x\right)}^T{\Gamma }^{-1}G\left(x\right)$

其中x表示大气各状态量，G(x)代表表征台风动力特征的平衡方程；而Г^-1为罚值系数(或称 动力权重)，其数值的大小决定了约束条件在同化中所起的作用。

3.动力权重值的估计

为了使得观测罚项和约束罚项在最小化迭代过程中起到相当的作用，通过量纲分析法来 估计动力权重Г^-1的大小。在动力方程G(x)中，平流项和气压梯度项为较重要的贡献项，由 此二项的量纲分析来确定方程G(x)的量纲。

对于台风天气系统，水平特征尺度L～10⁵m，水平风速平方的特征尺度U²～10³m²s^-2，气 压梯度的特征尺度Δp～10³Pa。因此，平流项的量级为气压梯度项量级也为 这意味着，如果动力权重Г^-1设为10⁴，则动力约束罚项J_c的量级大约为1，与观 测罚项一致。由于动力约束方程G(x)是近似的大气运动方程，具有无粘和定常等假设，为了 将该约束设为较弱的约束，本发明采用10³作为默认的动力权重值，使得动力约束罚项在同 化中的作用等于或略小于观测罚项。

4.动力平衡约束下的变分同化方法求解过程

WRFDA变分同化采用共轭梯度法(Navon and Legler1987)求解代价函数的极小值，在 最小化过程的每一次迭代步都需要计算代价函数各项的梯度，

$▿J=▿J_o+{\beta }_1▿J_b+{\beta }_2▿J_e+▿J_c$

其中约束罚项J_c的梯度表示为，

$▿J_c=G^T{\Gamma }^{-1}[{\mathrm{Gx}}^{\prime }+G\left(x_b\right)]$

其中x′表示大气各状态量的切线性增量，x_b表示背景场，G为动力平衡方程G的切线性算子， G^T为其对应的伴随算子。切线性算子G表示为，

式中上横线表示背景态，上标表示增量。

在最小化过程的每个迭代步中，首先完成从控制变量空间v向模式变量空间x′的转换 (Barker et al.2003)，模式变量空间的物理量x′包含水平风速u和v，温度t，混合比q和地 表气压ps的增量。对于切线性方程G，需要另外诊断位势高度φ′和气压p′的增量(Skamarock  et al.2008)。

气压增量p′由地表干气压增量μ_d′和水汽混合比增量q′诊断，

${\mathrm{dp}}^{\prime }={{\mu }_d}^{\prime }(1+\bar{q})\mathrm{d\eta }+q^{\prime }\overline{{\mu }_d}\mathrm{d\eta }$

气压增量p′的三维结构通过将上式从模式顶层至底层进行垂直积分而得。

其中地表干气压μ_d′通过地表气压p_s′除掉湿度的影响计算而得，

${{\mu }_d}^{\prime }=\frac{{p_s}^{\prime }-\overline{{\mu }_d}{\int }_0^1{q}^{\prime }\mathrm{d\eta }}{{\int }_0^1(1+\bar{q})\mathrm{d\eta }}$

对于位势高度φ′的诊断，利用静力平衡关系以及地表干气压增量μ_d′和干空气比容增量计算而得，

${\mathrm{d\phi }}^{\prime }=(-\overline{{\mu }_d}{\alpha }_d^{\prime }-\overline{{\alpha }_d}{\mu }_d^{\prime })\mathrm{d\eta }$

位势高度增量φ′的三维结构可通过将上式从模式底层至顶层进行垂直积分计算而得。

通过上述计算完成约束方程的切线性算子，之后计算其对应的伴随算子G^T，求得约束项 梯度与WRFDA变分同化系统现有的观测罚项和背景罚项的梯度相加，求得新的代价 函数梯度并完成模式变量空间x′向控制变量空间v的转换。通过共轭梯度法进行最优化 迭代搜索，完成一次迭代。最小化过程的结束条件为代价函数梯度小于初始代价函数梯度的 1/100。改进的同化方案流程示意图见图1，灰色方框内为新加入的动力平衡约束模块。

5.改进同化方案的迭代收敛性检验

图2验证了加入动力约束罚项后的代价函数最小化迭代过程，总的代价函数(黑线)以 及观测罚项(红线)和约束罚项(蓝线)的代价函数都下降较快(图2a)，与其对应的梯度 也下降正常(图2b)。

6.同化结果

应用原始的同化方案和改进的同化方案，分别对2010年03号台风“灿都”进行同化试 验。图3表明同化之前背景场的环流很弱(图3a)，海平面最低气压为994hPa(图3d)，远 高于观测值的970hPa，且台风中心与观测相差较远。同化雷达资料T-TREC风场后，原始同 化方案和改进同化方案对应实验的风场环流都有显著改善(图3b，c)。然而，原始同化方案 的海平面最低气压只有988hPa(图3e)，与较强的台风涡旋风场不匹配；改进同化方案的海 平面最低气压下降到983hPa，与观测值最为接近，同时海平面最低气压中心也与观测位置十 分吻合(图3f)。这表明，在原始的变分同化方案中加入台风动力平衡约束，有助于平衡分析 场中的气压和风场结构，使得分析场更为准确。

以背景场、原始同化方案的分析场和改进同化方案的分析场分别作为预报初始场，对台 风“灿都”进行24小时预报。图4表明，改进同化方案对应的预报实验在原始同化方案的基 础上，进一步的改善了台风强度预报水平。原始同化方案的海平面最低气压和地面最大风速 平均预报误差分别为8.3hPa和4.4ms^-1，；改进同化方案对应的预报误差分别下降至6.9hPa 和4.1ms^-1。表明台风动力约束下的同化方案在改善台风同化效果的基础上，有助于台风强 度预报水平的提高。