海洋学报

方法	积分窗口长度	集合成员个数	观测采样间隔	参数标准值
解析4DEnVar	500	500	10	10	28	8/3
传统4DVar	500	0	10	10	28	8/3

方法

积分窗口
长度

集合成员
个数

观测采样
间隔

参数标准值

$\sigma$

解析4DEnVar

500

8/3

传统4DVar

500

8/3

方法	积分窗口长度	集合成员个数	观测采样间隔	参数标准值
解析4DEnVar	500	500	10	10	28	8/3
传统4DVar	500	0	10	10	28	8/3

方法

积分窗口
长度

集合成员
个数

观测采样
间隔

参数标准值

$\sigma$

解析4DEnVar

500

8/3

传统4DVar

500

8/3

方法	积分窗口长度	集合成员数量	观测采样间隔	参数标准值
解析4DEnVar	100	500	5	10	20	10	28	8/3
200	500	5	10	20	10	28	8/3
传统4DVar	300	0	5	10	20	10	28	8/3
400	0	5	10	20	10	28	8/3

方法

积分窗口长度

集合成员数量

观测采样间隔

参数标准值

$\sigma$

解析4DEnVar

100

500

8/3

200

500

8/3

传统4DVar

300

8/3

400

8/3

方法	积分窗口长度	集合成员数量	观测采样间隔	参数标准值
解析4DEnVar	100	500	5	10	20	10	28	8/3
200	500	5	10	20	10	28	8/3
传统4DVar	300	0	5	10	20	10	28	8/3
400	0	5	10	20	10	28	8/3

方法

积分窗口长度

集合成员数量

观测采样间隔

参数标准值

$\sigma$

解析4DEnVar

100

500

8/3

200

500

8/3

传统4DVar

300

8/3

400

8/3

方法	积分窗口长度	集合成员个数	观测采样间隔	参数标准值
解析4DEnVar	500	10	10	10	28	8/3
500	100	10	10	28	8/3
500	200	10	10	28	8/3
500	300	10	10	28	8/3

方法

积分窗口
长度

集合成员
个数

观测采样
间隔

参数标准值

$\sigma $

解析4DEnVar

500

8/3

500

100

8/3

500

200

8/3

500

300

8/3

方法	积分窗口长度	集合成员个数	观测采样间隔	参数标准值
解析4DEnVar	500	10	10	10	28	8/3
500	100	10	10	28	8/3
500	200	10	10	28	8/3
500	300	10	10	28	8/3

方法

积分窗口
长度

集合成员
个数

观测采样
间隔

参数标准值

$\sigma $

解析4DEnVar

500

8/3

500

100

8/3

500

200

8/3

500

300

8/3

方法	积分窗口长度	集合成员个数	观测采样间隔	参数标准值选取范围
解析4DEnVar	500	500	10	（9, 11）	（25.2, 30.8）	（2.4, 2.93）

方法

积分窗口长度

集合成员个数

观测采样间隔

参数标准值选取范围

$\sigma $

解析4DEnVar

500

（9, 11）

（25.2, 30.8）

（2.4, 2.93）

方法	积分窗口长度	集合成员个数	观测采样间隔	参数标准值选取范围
解析4DEnVar	500	500	10	（9, 11）	（25.2, 30.8）	（2.4, 2.93）

方法

积分窗口长度

集合成员个数

观测采样间隔

参数标准值选取范围

$\sigma $

解析4DEnVar

500

（9, 11）

（25.2, 30.8）

（2.4, 2.93）

解析四维集合变分参数优化方法研究

PDF下载

贾彬鹤 ¹ , 李威 ¹^,^* , 梁康壮 ¹^,^*

海洋学报 | 论文 2021,43(10): 61-69

收起

海洋学报 | 论文 2021, 43(10): 61-69

解析四维集合变分参数优化方法研究

全屏

贾彬鹤¹, 李威¹^,^*, 梁康壮¹^,^*

作者信息

¹天津大学海洋科学与技术学院，天津 300072

贾彬鹤（1996-），男，山西省运城市人，主要从事海洋数据同化研究。E-mail：775341237@qq.com

通讯作者:

李威，男，教授，主要从事业务化海洋学研究。E-mail：liwei_nmdis@163.com

梁康壮，男，博士，主要从事业务化海洋学研究。E-mail：liang_kz@tju.edu.cn

Research on the optimization method of analytical four dimensional ensemble variational data assimilation

Binhe Jia¹, Wei Li¹^,^*, Kangzhuang Liang¹^,^*

Affiliations

¹School of Marine Science and Technology, Tianjin University, Tianjin 300072, China

出版时间: 2021-10-25 doi: 10.12284/hyxb2021129

文章导航

摘要

收起

传统的四维变分数据同化方法在同化观测资料的同时可以对数值模式参数进行优化，然而传统的四维变分方法需要针对不同的数值模式编写特有的伴随模式，因此算法的可移植性差，同时计算时耗费大量资源。本文提出了一种新的基于解析四维集合变分的参数优化方法，该方法以迭代搜索得到的模式参数为基准展开扰动并构建样本集合，由此显式地计算协方差矩阵，并得到代价函数极小值的解析解，从而避免了伴随模式的使用。基于Lorenz-63模型对该方法进行单参数和多参数数值试验和优化效果检验，并在不同的同化时间窗口长度和观测采样间隔情况下，采用传统四维变分方法与之进行对比，结果显示，新方法表现出与传统四维变分相同的优化性能，都能有效收敛到真值，而新方法不需要计算伴随模式，可移植性好。本文还测试了不同的集合成员个数和模式参数真值的情况下新方法的同化效果，结果表明，新方法对集合样本个数及模型参数真值不敏感，采用较少的集合样本即可完成数据同化。

关键词

解析四维集合变分 / 参数优化 / 传统四维变分 / Lorenz-63

Abstract

收起

The traditional four-dimensional variational data assimilation method can optimize the parameters of the numerical model while assimilating the observation data. However, the traditional four-dimensional variational method needs to compile special adjoint models for different numerical models, so the portability of the traditional four-dimensional variational method is poor and a lot of resources are consumed in the calculation. In this paper, a new parameter optimization method based on the analytic four-dimensional ensemble variation is proposed, which expands the perturbation and constructs the ensemble based on the model parameters obtained by iterative search, and then explicitly calculates the covariance matrix, and obtains the analytic solution of the minimum value of the cost function, so as to avoid the use of adjoint model. Using Lorenz-63 model, single-parameter and multi-parameter numerical tests and optimization effect tests were carried out on the analytic four-dimensional ensemble variation method, and in the case of different assimilation time window length and observation sampling interval, the traditional four-dimensional variational method was used to compare with the new method, the results show that the new method has the same optimization performance as the traditional four-dimensional variational method, and it can converge to the truth value effectively, and the new method does not need to calculate adjoint mode, so it has good portability. This paper also test the assimilation effect of the new method with different ensemble members and true values of model parameters, and the results show that the new method is insensitive to the number of ensemble members and the true values of model parameters, and the data assimilation can be completed with fewer ensemble members.

Key words

analytical four dimensional ensemble variational data assimilation / parameter optimization / traditional four dimensional variational data assimilation / Lorenz-63

引用本文

贾彬鹤, 李威, 梁康壮. 解析四维集合变分参数优化方法研究. 海洋学报, 2021 , 43 (10) : 61 -69 . DOI: 10.12284/hyxb2021129

Binhe Jia, Wei Li, Kangzhuang Liang. Research on the optimization method of analytical four dimensional ensemble variational data assimilation[J]. Haiyang Xuebao, 2021 , 43 (10) : 61 -69 . DOI: 10.12284/hyxb2021129

正文

收起

1 引言

收起

进行数值预报时，预报的准确率会随着预报时间的增长而迅速下降。造成这种现象的原因主要有3点：一是模式方程对物理过程的参数化方案有一定误差；二是数值预报模式的初始场有误差；三是数值预报的计算方法有一定的误差^[1]。

为提高预报准确率，如何改进模式物理参数是需要解决的问题。数据同化方法可以从观测数据中提取有效信息，并修正数值模式误差。四维变分数据同化（Four Dimensional Variational Data Assimilation，4DVar）方法能够把各个不同时刻的观测资料纳入统一的分析预报系统中，通过相应的切线性模式和伴随模式优化初值条件，并使各要素自然满足动力热力条件^[2-3]。由于其较好的同化能力，众多学者也采用4DVar方法进行了模式参数优化的研究^[4-7]。例如，Inazu等^[8]将进化算法应用于区域海潮模型，优化了模型的边界条件和物理参数。Wang等^[9]提出通过遗传算法和神经网络交替组合的方式来纠正系统参数误差。王云峰等^[10]提出了一种新的利用观测资料来同时优化模式初始场和物理参数的扩展四维变分同化方法，并以Ekman边界层模式和Lorenz模型为例进行了数值试验。魏敏^[1]采用4DVar方法，利用Lorenz模型和Hamilton方程，结合观测数据对该模型的参数进行了优化。然而，4DVar方法需要采用伴随模式，对于原数值模式需要编写与之相对应的特定的伴随模式，并且伴随复杂程度高、代码编辑量大，使得这种同化方法可移植性很差。另一方面，4DVar采用的静态背景误差协方差矩阵不具有流依赖特性，限制了四维变分的同化能力。

序列同化方法如集合卡尔曼滤波（Ensemble Kalman Filter，EnKF）能在避免伴随的同时产生流依赖的背景误差协方差矩阵。为结合两种方法的优点，可以在变分框架下引入集合状态变量，从而产生四维集合变分（Four Dimensional Ensemble Variational Data Assimilation，4DEnVar）。最初，Liu等^[11] 将四维变分增量优化方法与集合方法结合提出了该方法，但将其命名为集合四维变分数据同化方法（Ensemble Four Dimensional Variational Data Assimilation，En4DVar），并采用一维浅水模型对该方法进行了验证。随后，Liu等^[12]设计了局地化方法，采用高等天气研究与预报（ARW-WRF）模型，设计了观测系统模拟试验，并检验了它们在实际数据同化中的性能，证实了4DEnVar 局地化技术可以有效地减轻采样误差对分析的影响。根据世界气象组织（WMO）会议建议，Liu和Xiao^[13]将En4DVar更名为4DEnVar，以便将4DEnVar和使用伴随模型的En4DVar区分开来。4DEnVar 方法在业务化系统中也有所应用。Arbogast等^[14]对四维轨迹的输入和存储等问题，提出了一种具有分布式输入和存储扰动的四维集合变分集成并行实现办法。Yang和Mémin^[15]探索了一个动力学公式，该公式允许在一个大规模的流动模型中同化高分辨率的数据，并将该原理结合4DEnVar技术来估计随机浅水模型的流动初始条件和非均匀时变子网格参数。Song和Kang^[16]采用混合四维集合变分方法（hybrid-4DEnVar），改进了北半球500 hPa夏季月位势高度预报。杨雨轩^[17]利用WRF模型研究了POD-4DEnVar和传统三维变分方法（Three Dimensional Variational Data Assimilation，3DVar）对华南暴雨的模拟效果，并研究了物理集合方案对4DEnVar的影响。

4DEnVar可以构造多个能反映出背景误差协方差分布特征的样本形成集合，为变分同化系统提供流依赖背景误差协方差估计。然而，需要指明的是，Liu等^[11]将集合成员扰动关于集合均值展开来构造的协方差矩阵，反映了集合的统计性质，而非误差的传播性质。在线性模型下，这两种协方差矩阵等价，但对于非线性模型则不然。此外，目前4DEnVar方法研究集中在初始场的优化，尚未有进行参数优化的研究。为此，本文提出了专门进行模式参数优化的解析四维集合变分数据同化（Analytical Four Dimensional Ensemble Variational Data Assimilation，A-4DEnVar）方法，该方法在4DVar框架下，将集合样本扰动关于迭代生成的模式参数处展开，有效刻画了误差在非线性模型下的传播方式，避免了伴随模式的使用。在此基础上，进一步求得代价函数梯度为0的解析解，能迅速得到代价函数极小值，完成模式参数优化。本文采用Lorenz-63模型^[18]通过4DVar与A-4DEnVar 方法进行参数优化孪生试验和敏感性试验，证实算法的有效性。

2 解析四维集合变分参数优化方法

收起

假设动力系统为

(1) ${{\boldsymbol{X}}_i} = {{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}\left[ {{{\boldsymbol{X}}_{i - 1}},{\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}} \right],$

式中，${{\boldsymbol{X}}_{i - 1}}$为第$\left( {i - 1} \right)$时刻的状态变量；${{\boldsymbol{X}}_i}$为第$i$时刻的状态变量；${{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}$为从时刻$\left( {i - 1} \right)$到时刻$i$的演化算符；${{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}$中包含了有偏差的模式参数${\boldsymbol{\varLambda }}$。按照式（1）可以递推得到如下方程

(2) ${{\boldsymbol{X}}_i} = {{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}\left[ { \cdots {{\boldsymbol{M}}_{0 \to 1}}\left( {{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{0 \to 1}}} \right) \cdots ,{{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}} \right].$

假设初始场${{\boldsymbol{X}}_0}$与驱动场${{\boldsymbol{F}}_i}$都是准确的，唯一不准确的是模式参数${\boldsymbol{\varLambda }}$，在${\boldsymbol{\varLambda }}$上叠加一个小扰动${\mathbf{\tilde \lambda }}$得到另外的一个模式参数为${\mathbf{\lambda }} \equiv {\boldsymbol{\varLambda }} + {\mathbf{\tilde \lambda }}$，可以预期，只要积分步长不是很长，受新的模式参数${\mathbf{\lambda }}$影响，各个时刻状态变量可以表示为原状态变量${{\boldsymbol{X}}_i}$与扰动量${{\boldsymbol{\tilde x}}_i}$相加的形式${{\boldsymbol{X}_i} + {{\boldsymbol{\tilde x}}}_i}$，则可以得到扰动量所满足的方程

(3) ${{\boldsymbol{\tilde x}}_i} \approx {\left. {{{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}} \right|_{{{\boldsymbol{X}}_{i - 1}},{\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}}}_{\left( {\boldsymbol{x}} \right)}{{\boldsymbol{\tilde x}}_{i - 1}} + {\left. {{{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}} \right|_{{{\boldsymbol{X}}_{i - 1}},{\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}}}_{\left( {\mathbf{\lambda }} \right)}{\mathbf{\tilde \lambda }},$

式中，${\left. {{{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}} \right|_{{{\boldsymbol{X}}_{i - 1}},{\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}}}_{\left( {\boldsymbol{x}} \right)}$是${{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}\left[ \bullet \right]$算符在${{\boldsymbol{X}}_{i - 1}}$、${\boldsymbol{\varLambda }}$和${{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}$附近对${{\boldsymbol{X}}_{i - 1}}$的一阶偏导数；${\left. {{{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}} \right|_{{{\boldsymbol{X}}_{i - 1}},{\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}}}_{\left( {\mathbf{\lambda }} \right)}$是${{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}\left[ \bullet \right]$算符在${{\boldsymbol{X}}_{i - 1}}$、${\boldsymbol{\varLambda }}$和${{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}$对${\boldsymbol{\varLambda }}$的一阶偏导数。需要注意的是，${\left. {{{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}} \right|_{{{\boldsymbol{X}}_{i - 1}},{\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}}}_{\left( {\mathbf{x}} \right)}$和${\left. {{{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}} \right|_{{{\boldsymbol{X}}_{i - 1}},{\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}}}_{\left( {\mathbf{\lambda }} \right)}$的矩阵元素的数值虽然取决于状态变量${{\boldsymbol{X}}_{i - 1}}$、模式参数${\boldsymbol{\varLambda }}$和驱动场${{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}$的数值，但根据递推关系，其归根到底都取决于以${\boldsymbol{\varLambda }}$为模式参数和以${{\boldsymbol{X}}_0}$为初始场在${\boldsymbol{F}}$驱动下积分得到的状态变量值，因此，如果假设$n < m$，则可以定义

(4) ${\left. {{{\boldsymbol{M}}_{n \to m}}} \right|_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}} \equiv {\left. {{{\boldsymbol{M}}_{\left( {m - 1} \right) \to m}}} \right|_{{{\boldsymbol{X}}_{m - 1}},{\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{\left( {m - 1} \right) \to m}}}}_{\left( {\boldsymbol{x}} \right)} \cdots {\left. {{{\boldsymbol{M}}_{n \to \left( {n + 1} \right)}}} \right|_{{{\boldsymbol{X}}_n},{\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{\left( {n - 1} \right) \to n}}}}_{\left( {\boldsymbol{x}} \right)},$

和

(5) ${\left. {{{\boldsymbol{\varGamma }}_{n \to \left( {n + 1} \right)}}} \right|_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}} \equiv {\left. {{{\boldsymbol{M}}_{n \to \left( {n + 1} \right)}}} \right|_{{{\boldsymbol{X}}_n},{\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{\left( {n - 1} \right) \to n}}}}_{\left( {\mathbf{\lambda }} \right)},$

于是演化公式可改为

(6) $\begin{split} {{\boldsymbol{\tilde x}}_i} \approx & \left[ {{\left. {{{\boldsymbol{M}}_{1 \to i}}} \right|}_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}}{{\left. {{{\boldsymbol{\varGamma }}_{0 \to 1}}} \right|}_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}} + \cdots + {{\left. {{{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}} \right|}_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}}{{\left. {{{\boldsymbol{\varGamma }}_{\left( {i - 2} \right) \to \left( {i - 1} \right)}}} \right|}_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}} +\right.\\& \qquad\qquad\qquad\left.{{\left. {{{\boldsymbol{\varGamma }}_{\left( {i - 1} \right) \to i}}} \right|}_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}} \right]{\mathbf{\tilde \lambda }}.\\[-12pt]\end{split}$

注意到，此处假设初始场和外界强迫场都是准确的，即初始场和外界强迫场都不扰动，在这种情况下，状态变量的扰动都来自于模式参数的扰动。

定义广义背景场误差协方差矩阵为

(7) $\begin{split} {\left. {{{\boldsymbol{D}}_i}} \right|_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}} \equiv & {\left. {{{\boldsymbol{M}}_{1 \to i}}} \right|_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}}{\left. {{{\boldsymbol{\varGamma }}_{0 \to 1}}} \right|_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}} + \cdots + {\left. {{{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}} \right|_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}}\\& {\left. {{{\boldsymbol{\varGamma }}_{\left( {i - 2} \right) \to \left( {i - 1} \right)}}} \right|_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}} + {\left. {{{\boldsymbol{\varGamma }}_{\left( {i - 1} \right) \to i}}} \right|_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}},\end{split}$

式中，矩阵实际上包含了切线性演化算符，其转置矩阵则包含了伴随算符信息，于是，

(8) ${{\boldsymbol{\tilde x}}_i} \approx {\left. {{{\boldsymbol{D}}_i}} \right|_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}}{\mathbf{\tilde \lambda }}.$

如果只通过扰动参数构造集合，则可以在最小二乘意义下，显式地求解出式（7）矩阵的具体形式，即

(9) ${\left. {{{\boldsymbol{D}}_i}} \right|_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}} = E\left\{ {{{{\boldsymbol{\tilde x}}}_i}{{{\mathbf{\tilde \lambda }}}^{\rm{T}}}} \right\}E{\left\{ {{\mathbf{\tilde \lambda }}{{{\mathbf{\tilde \lambda }}}^{\rm{T}}}} \right\}^{ - 1}}.$

如果只优化模式参数，则四维变分的目标函数为

(10) $J\left( {\boldsymbol{\varLambda }} \right) = \frac{1}{2}\sum\limits_i {{{\left( {{{\boldsymbol{H}}_i}{{\boldsymbol{X}}_i} - {{\boldsymbol{Y}}_i}} \right)}^{\rm{T}}}{\boldsymbol{R}}_i^{ - 1}\left( {{{\boldsymbol{H}}_i}{{\boldsymbol{X}}_i} - {{\boldsymbol{Y}}_i}} \right)} ,$

(11) ${{\boldsymbol{X}}_i} = {{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}\left[ { \cdots {{\boldsymbol{M}}_{0 \to 1}}\left( {{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{0 \to 1}}} \right) \cdots ,{{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}} \right],$

式中，${{\boldsymbol{Y}}_i}$为第$i$时刻的观测；${{\boldsymbol{H}}_i}$为第$i$时刻的观测投影算符；${{\boldsymbol{R}}_i}$为第$i$时刻的观测误差协方差矩阵。

传统四维变分对模式参数进行最优估计的时候，就是将目标函数对模式参数${\boldsymbol{\varLambda }}$求梯度，并将目标函数值和梯度值代入到最优化算法中，通过线性搜索和逐步迭代得到最优的模式参数。仿照这一过程，在当前模式参数${\boldsymbol{\varLambda }}$附近进行扰动展开，推导${\boldsymbol{\varLambda }}$附近使得目标函数取极小值的最优扰动量，设扰动量为$\delta {\boldsymbol{\varLambda }}$，将${{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}\left[ { \cdots {{\boldsymbol{M}}_{0 \to 1}}\left( {{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }} + \delta {\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{0 \to 1}}} \right) \cdots ,{{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}} \right]$在当前模式参数${\boldsymbol{\varLambda }}$附近进行展开，由于扰动量是小量，忽略高阶项后，可以得到

(12) $\begin{split}&{{\boldsymbol{M}}_{\left( {i - 1} \right) \to i}}\left[ { \cdots {{\boldsymbol{M}}_{0 \to 1}}\left( {{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }} + \delta {\boldsymbol{\varLambda }},{{\boldsymbol{F}}_{0 \to 1}}} \right) \cdots ,{{\boldsymbol{F}}_{\left( {i - 1} \right) \to i}}} \right] \\&\qquad\qquad\quad\approx {{\boldsymbol{X}}_i} + {\left. {{{\boldsymbol{D}}_i}} \right|_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}}\delta {\boldsymbol{\varLambda }},\end{split}$

将式（12）代入式（10），则四维变分目标函数为

(13) $\begin{split} J\left( {{\boldsymbol{\varLambda }} + \delta {\boldsymbol{\varLambda }}} \right) =& \frac{1}{2}\sum\limits_i {{\left[ {{{\boldsymbol{H}}_i}{{\left. {{{\boldsymbol{D}}_i}} \right|}_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}}\delta {\boldsymbol{\varLambda }} - \left( {{{\boldsymbol{Y}}_i} - {{\boldsymbol{H}}_i}{{\boldsymbol{X}}_i}} \right)} \right]}^{\rm{T}}}\times\\&{\boldsymbol{R}}_i^{ - 1}\left[ {{{\boldsymbol{H}}_i}{{\left. {{{\boldsymbol{D}}_i}} \right|}_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}}\delta {\boldsymbol{\varLambda }} - \left( {{{\boldsymbol{Y}}_i} - {{\boldsymbol{H}}_i}{{\boldsymbol{X}}_i}} \right)} \right] ,\end{split}$

接着求出式（13）中目标函数相对于扰动量$\delta {\boldsymbol{\varLambda }}$的梯度，

(14) $\nabla J\left( {\delta {\boldsymbol{\varLambda }}} \right) = \sum\limits_i {{{\left. {{\boldsymbol{D}}_i^{\rm{T}}} \right|}_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}}{\boldsymbol{H}}_i^{\rm{T}}{\boldsymbol{R}}_i^{ - 1}\left[ {{{\boldsymbol{H}}_i}{{\left. {{{\boldsymbol{D}}_i}} \right|}_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}}\delta {\boldsymbol{\varLambda }} - \left( {{{\boldsymbol{Y}}_i} - {{\boldsymbol{H}}_i}{{\boldsymbol{X}}_i}} \right)} \right]} ,$

令式（14）梯度为0，即$\nabla J\left( {\delta {\boldsymbol{\varLambda }}} \right) = 0$，则可以求出$\delta {\boldsymbol{\varLambda }}$的解析表达式

(15) $\begin{split} \delta {\boldsymbol{\varLambda }} =& {\left( {\sum\limits_i {{{\left. {{\boldsymbol{D}}_i^{\rm{T}}} \right|}_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}}{\boldsymbol{H}}_i^{\rm{T}}{\boldsymbol{R}}_i^{ - 1}{{\boldsymbol{H}}_i}{{\left. {{{\boldsymbol{D}}_i}} \right|}_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}}} } \right)^{ - 1}}\times\\& \sum\limits_j {{{\left. {{\boldsymbol{D}}_j^{\rm{T}}} \right|}_{{{\boldsymbol{X}}_0},{\boldsymbol{\varLambda }},{\boldsymbol{F}}}}{\boldsymbol{H}}_j^{\rm{T}}{\boldsymbol{R}}_j^{ - 1}\left( {{{\boldsymbol{Y}}_j} - {{\boldsymbol{H}}_j}{{\boldsymbol{X}}_j}} \right)} .\end{split}$

需要指出的是，这一最优增量是在${\boldsymbol{\varLambda }}$为初猜场的情况下得到的，对于非线性动力系统，只要稍微有一点改动，上述最优增量的值就必然会发生变化。因此，为了让这一方法具有一定的稳定性，需要仿照传统四维变分的做法，引入线性搜索过程，即在相空间的$\delta {\boldsymbol{\varLambda }}$方向上，以一定的步长$\beta $线性搜索（例如采用二分法）使得目标函数在该方向上取极小的最优值$\beta \cdot \delta {\boldsymbol{\varLambda }}$作为订正量的最优取值，将${\boldsymbol{\varLambda }} + \beta \cdot \delta {\boldsymbol{\varLambda }}$作为下一次迭代的猜测值代入目标函数，继续计算解析解，循环往复上述过程从而获得最优的模式参数。

3 数值试验

收起

1963年美国麻省理工学院的气象学家Lorenz在对天气预报动力学模型进行数值计算时，发现了一个由非线性微分方程组所描述的著名的 Lorenz方程^[18]，该方程是一个3个变量模型，是对流非周期模式的简化，同时也是混沌理论的原型个例^[10]。由于Lorenz方程具有较强的非线性，且模型较为简洁，所以常采用Lorenz-63模型来检验数据同化算法的性能^[19-20]。

Lorenz方程为

(16) $\left\{ \begin{aligned}& \frac{{{{\rm{d}}{{x}}}}}{{{{\rm{d}}{{t}}}}} = {{\sigma }}({{y}} - {{x}}),\\& \frac{{{{\rm{d}}{{y}}}}}{{{{\rm{d}}{{t}}}}} = {{rx}} - {{y}} - {{xz}},\\& \frac{{{{\rm{d}}z}}}{{{{\rm{d}}{{t}}}}} = {{xy}} - {{bz}},\end{aligned} \right.$

式中，${\rm{\sigma }} > 0$，${{b}} > 0 $；参数${\rm{\sigma }}$、${{r}}$、${{b}}$、${{t}}$分别表示Prantl数、Rayleigh数、对流尺度联系的参数和时间； ${{x}}$是对流强度；${{y}}$是最大温度差；${{z}}$是对流引起的层变化^[21]。模型标准参数一般为${\rm{\sigma }} = 10,\;{{r}} = 28,\;{{b}} = 8/3$，从而产生混沌解^[22]。

本文采用Lorenz-63模型进行数值试验。首先，对标准参数模型采用四阶荣格库塔（Runge-Kutta）方法积分生成真实场（积分步长设置为${{\rm{d}}t} = 0.01$）；随后在一个同化时间窗口内以一定的采样间隔进行观测采集，并在观测数据上加入方差为1的白噪声，生成带噪声的观测数据。具体优化步骤为：根据式（15）计算每次优化迭代之后的参数最优扰动方向；确定最优扰动步长，将计算出的最优扰动叠加到当前参数值上得到下一次优化迭代的参数初值；经过多次迭代，得到收敛的模型参数值，作为最终参数优化结果。本文设置单参数、双参数和三参数试验来检测新方法的模型参数优化能力，并与传统4DVar方法进行模型参数优化对比试验；此外，通过设置不同的模型参数真值、集合成员个数、同化时间窗口长度及观测采样间隔等敏感性试验，考察新方法的参数估计能力。

3.1 新方法参数优化结果

本节试验积分窗口取500步，集合成员个数取500个，将${\rm{\sigma }} = 10,\;{{r}} = 28,\;{{b}} = 8/3$这些特殊参数值作为本次试验的标准参数值，取背景场状态变量为${{{x}}_0} = 1, {{{y}}_0} = 2,\;{{{z}}_0} = 3$对模式进行积分构造真实场，以10步为采样间隔采集观测并加入白噪声构造观测场。假设初猜的参数值为${{\rm{\sigma }}_1} = 9,\;{{{r}}_1} = 25,\;{{{b}}_1} = 5$，并采用相同的背景场状态变量向前积分500步。传统4DVar方法采用了同样的参数设置，编辑伴随模式计算代价函数梯度并采用同样的线性搜索方法优化迭代计算。表1为本试验的主要参数配置。

本节设计7种试验方案，包括单参数试验、双参数试验以及三参数试验用以检验新方法的模型参数优化能力。结果表明，针对不同个数的参数，新方法均能收敛到真值，对三参数优化后的模型分析场与真实场如图1所示（其余结果未列出）。由图1可知，A-4DEnVar参数优化方法和传统的4DVar参数优化方法均可达到一个理想的优化效果，分析场拟合真实场程度高。各试验中代价函数随迭代变化情况如图2所示。由图2可知，A-4DEnVar在优化迭代过程中代价函数值在100步左右可以实现收敛，收敛结果较为稳定，因而具有良好的模型参数优化能力。

3.2 同化时间窗口长度和观测采样间隔影响

为了测试同化时间窗口长度和观测采样间隔对新方法在模型参数优化方面能力的影响，本节采用与3.1节相同的初始场和初猜参数，选取4组不同的同化时间窗口长度进行试验，并计算分析场与真实场的均方根误差（RMSE）。同化时间窗口长度为：100步、200步、300步和400步，在这4个时间窗口下再分别取5步、10步和20步的观测采样间隔（Observation Sampling Interval，OSI）。为方便对比对比，传统4DVar方法采用了同样的参数设置进行试验。

图3给出了改变积分时间窗口长度和观测采样间隔后的模型变量的平均RMSE变化。试验结果表明，针对不同的同化窗口和采样间隔，A-4DEnVar方法均可实现模型参数的理想优化，效果与传统4DVar方法相当；新方法的模型参数优化速度较传统的4DVar方法快，这是由于A-4DEnVar方法采用了解析解而不是梯度下降方向作为搜索方向。图4描绘了不同的同化时间窗口和观测采样间隔下优化迭代收敛后的RMSE。结果表明，采用新方法与传统4DVar方法均能使分析场RMSE低于观测误差标准差，可实现模型参数的收敛，同一种方法在同一采样间隔的情况下，随着积分窗口长度的增长，最终收敛的平均RMSE会增大；同一种方法在同一积分窗口长度的情况下，随着采样间隔的增大，最终收敛的平均RMSE会增大。

3.3 集合成员数量影响

为了测试集合成员数量的大小对新方法参数优化能力的影响，本节采用与3.1节相同的初始场和初猜参数及同化时间窗口和参数标准值，选取4组不同的集合成员数量进行试验。4组集合成员数量分别是：10、100、200和300。试验参数如表3所示。

图5展示了不同集合成员数量情况下模型分析场与真实场的拟合程度。由图5可知，在不同的集合成员数量下，新方法均可实现参数的理想收敛，模型分析场拟合真实场程度高，这说明新方法对集合成员数量不敏感；本节最小集合成员个数取的是10，若再取更小的集合成员个数（例如3或4），则需要更多的迭代次数来实现模型参数的收敛。

3.4 参数标准值影响

为了测试不同参数标准值对提出方法的影响，本节采用与3.1节相同的初始场和初猜参数及同化时间窗口和集合成员个数，在模式原标准参数上加上一定的扰动倍数，形成3组不同的参数标准值进行组合试验^[10]。3种参数所属区间分别是：${\rm{\sigma }} \in (0.9 \times 10,1.1 \times 10)$，$ {{r}}\in (0.9\times 28,1.1\times 28)$和$ {{b}}\in (0.9\times 8/3,1.1\times 8/3)$，每个区间生成10个任意值，一共1 000组试验。试验参数如表4所示。

图6展示了改变模型参数标准值后的各组代价函数值随着迭代次数的变化，横坐标表示组别，纵坐标表示代价函数值。图7展示了改变模型参数标准值后平均代价函数值随着迭代次数的变化，横坐标表示迭代次数，纵坐标表示代价函数值。由图6和图7可知，对选取的1 000组参数标准值，在迭代50步之内，代价函数值即可收敛至最小值，此时模型参数收敛至标准值，由此可知新方法优化能力较强，且对参数的选取不敏感。

4 结论

收起

本文将模型参数视为一种特殊的控制变量，在传统4DVar框架下，以迭代得到的模式参数为基准展开集合扰动，计算误差协方差矩阵。在此基础上，使用代价函数最小值的解析解构造线性搜索最优化方法，提出了A-4DEnVar参数优化方法，避免了传统4DVar方法中的伴随模式的使用。为验证有效性，采用传统的4DVar方法和新的A-4DEnVar参数优化方法对Lorenz-63模型进行参数优化对比试验，针对不同个数的参数，新方法同化效果与传统方法4DVar相当。针对不同的同化时间窗口、观测采样间隔、集合样本个数及不同标准参数的敏感性试验结果显示，新方法能达到一个准确的收敛效果，在较长的同化时间窗口和观测采样间隔时也可以实现理想的模型参数优化效果，并且该方法对集合样本个数和标准参数不敏感，采用较少的集合样本即可达到同化效果，节约了工作量，这在实际模型的参数优化工作中，将具有重要意义。

基金

收起

国家自然科学基金（41876014）

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

魏敏. 四维变分方法在微分方程参数优化中的应用[D]. 南京: 南京师范大学, 2013.

Wei Min. Applications of the four-dimensional variation methods on parameter optimization of differential equations[D]. Nanjing: Nanjing Normal University, 2013.

Lewis J M, Derber J C. The use of adjoint equations to solve a variational adjustment problem with advective constraints[J]. Tellus A: Dynamic Meteorology and Oceanography, 1985, 37(4): 309−322.

Le Dimet F X, Talagrand O. Variational algorithms for analysis and assimilation of meteorological observations: Theoretical aspects[J]. Tellus A: Dynamic Meteorology and Oceanography, 1986, 38(2): 97−110.

Chu Kekuan, Tan Zhemin, Ming Xue. Impact of 4DVAR assimilation of rainfall data on the simulation of mesoscale precipitation systems in a Mei-Yu heavy rainfall event[J]. Advances in Atmospheric Sciences, 2007, 24(2): 281−300.

Zhao Juan, Wang Bin, Liu Juanjuan. Impact of analysis-time tuning on the performance of the DRP-4DVar approach[J]. Advances in Atmospheric Sciences, 2011, 28(1): 207−216.

Wang Yunfeng, Wang Bin, Fei Jianfang, et al. The effects of assimilating satellite brightness temperature and bogus data on the simulation of typhoon Kalmaegi (2008)[J]. Acta Meteorologica Sinica, 2013, 27(3): 415−434.

Zhong Jian, Huang Sixun, Fei Jianfang, et al. Application of tikhonov regularization method to wind retrieval from scatterometer data II: Cyclone wind retrieval with consideration of rain[J]. Chinese Physics B, 2011, 20(6): 064301.

Inazu D, Higuchi T, Nakamura K. Optimization of boundary condition and physical parameter in an ocean tide model using an evolutionary algorithm[J]. Theoretical and Applied Mechanics Japan, 2010, 58: 101−112.

Wang Tingting, Li Wenlong, Chen Zhanghui, et al. Correcting the systematic error of the density functional theory calculation: The alternate combination approach of genetic algorithm and neural network[J]. Chinese Physics B, 2010, 19(7): 076401.

王云峰, 顾成明, 张晓辉, 等. 优化模式物理参数的扩展四维变分同化方法[J]. 物理学报, 2014, 63(24): 12−19.

Wang Yunfeng, Gu Chengming, Zhang Xiaohui, et al. Expanded four-dimensional variatiaonal data assimilation method to optimize model physical parameters[J]. Acta Physica Sinica, 2014, 63(24): 12−19.

Liu Chengsi, Xiao Qingnong, Wang Bin. An ensemble-based four-dimensional variational data assimilation scheme. Part I: Technical formulation and preliminary test[J]. Monthly Weather Review, 2008, 136(9): 3363−3373.

Liu Chengsi, Xiao Qingnong, Wang Bin. An ensemble-based four-dimensional variational data assimilation scheme. Part II: Observing system simulation experiments with advanced research WRF (ARW)[J]. Monthly Weather Review, 2009, 137(5): 1687−1704.

Liu Chengsi, Xiao Qingnong. An ensemble-based four-dimensional variational data assimilation scheme. Part III: Antarctic applications with advanced research WRF using real data[J]. Monthly Weather Review, 2013, 141(8): 2721−2739.

Arbogast É, Desroziers G, Berre L. A parallel implementation of a 4DEnVar ensemble[J]. Quarterly Journal of the Royal Meteorological Society, 2017, 143(706): 2073−2083.

Yang Yin, Mémin E. High-resolution data assimilation through stochastic subgrid tensor and parameter estimation from 4DEnVar[J]. Tellus A: Dynamic Meteorology and Oceanography, 2017, 69(1): 1308772.

Song H J, Kang J H. Effects of the wind-mass balance constraint on ensemble forecasts in the hybrid-4DEnVar[J]. Quarterly Journal of the Royal Meteorological Society, 2019, 145(719): 434−449.

杨雨轩. 基于华南冬季暴雨的雷达资料四维集合变分同化技术研究[D]. 长沙: 国防科技大学, 2017.

Yang Yuxuan. Technical research of four-dimensional ensemble variational assimilation of doppler radar data based on awinter heavy rainstorm in south China[D]. Changsha: National University of Defense Technology, 2017.

Lee K S, Bang S H, Chang K S. Feedback-assisted iterative learning control based on an inverse process model[J]. Journal of Process Control, 1994, 4(2): 77−89.

杜川利, 黄向宇, 俞小鼎. 变分同化方法在Lorenz系统中的简单应用研究[J]. 气象, 2005, 31(2): 23−26.

Du Chuanli, Huang Xiangyu, Yu Xiaoding. Simple application of variational four-dimensional assimilation in Lorenz system[J]. Meteorological Monthly, 2005, 31(2): 23−26.

Hall M C G. Application of adjoint sensitivity theory to an atmospheric general circulation model[J]. Journal of the Atmospheric Sciences, 1986, 43(22): 2644−2652.

李建平, 丑纪范. 非线性大气动力学的进展[J]. 大气科学, 2003, 27(4): 653−673.

Li Jianping, Chou Jifan. Advances in nonlinear atmospheric dynamics[J]. Chinese Journal of Atmospheric Sciences, 2003, 27(4): 653−673.

郜吉东, 丑纪范. 数值模式初值的敏感性程度对四维同化的影响——基于Lorenz系统的研究[J]. 气象学报, 1995, 53(4): 471−479.

Gao Jidong, Chou Jifan. The effects of the model sensitivity to initial condition upon the variational four-diensional assimilation−the study based on Lorenz model[J]. Acta Meteorologica Sinica, 1995, 53(4): 471−479.

2021年第43卷第10期

PDF下载

158

引用本文

BibTeX

文章信息

doi: 10.12284/hyxb2021129

接收时间：2020-05-11
首发时间：2026-02-26
出版时间：2021-10-25

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2020-05-11
修回日期：2020-11-23

基金

国家自然科学基金（41876014）

作者信息

¹天津大学海洋科学与技术学院，天津 300072

通讯作者:

李威，男，教授，主要从事业务化海洋学研究。E-mail：liwei_nmdis@163.com

梁康壮，男，博士，主要从事业务化海洋学研究。E-mail：liang_kz@tju.edu.cn

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/hyxb/CN/10.12284/hyxb2021129

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT

方法	积分窗口长度	集合成员个数	观测采样间隔	参数标准值
解析4DEnVar	500	500	10	10	28	8/3
传统4DVar	500	0	10	10	28	8/3

方法

积分窗口
长度

集合成员
个数

观测采样
间隔

参数标准值

$\sigma$

解析4DEnVar

500

8/3

传统4DVar

500

8/3

方法	积分窗口长度	集合成员数量	观测采样间隔	参数标准值
解析4DEnVar	100	500	5	10	20	10	28	8/3
200	500	5	10	20	10	28	8/3
传统4DVar	300	0	5	10	20	10	28	8/3
400	0	5	10	20	10	28	8/3

方法

积分窗口长度

集合成员数量

观测采样间隔

参数标准值

$\sigma$

解析4DEnVar

100

500

8/3

200

500

8/3

传统4DVar

300

8/3

400

8/3

方法	积分窗口长度	集合成员个数	观测采样间隔	参数标准值
解析4DEnVar	500	10	10	10	28	8/3
500	100	10	10	28	8/3
500	200	10	10	28	8/3
500	300	10	10	28	8/3

方法

积分窗口
长度

集合成员
个数

观测采样
间隔

参数标准值

$\sigma $

解析4DEnVar

500

8/3

500

100

8/3

500

200

8/3

500

300

8/3

方法	积分窗口长度	集合成员个数	观测采样间隔	参数标准值选取范围
解析4DEnVar	500	500	10	（9, 11）	（25.2, 30.8）	（2.4, 2.93）

方法

积分窗口长度

集合成员个数

观测采样间隔

参数标准值选取范围

$\sigma $

解析4DEnVar

500

（9, 11）

（25.2, 30.8）

（2.4, 2.93）