1. 引言

资料同化是一种将观测资料融合于数值模式的分析技术，目的是通过有效利用一切信息(包括观测数据、背景场、模式以及相应的误差统计等)为非线性动力预报模型提供最优初始值 [1] [2] [3] [4] 。资料同化是一种具有普适性的通用技术，在海洋预报、气象预报、地震预测、电离层建模和流体等领域有广泛应用，因此对其方法的研究具有重要价值 [5] - [12] 。目前先进同化方法主要分为两大类：一类是变分资料同化方法 [3] [4] [10] [11] [12] ，另一类是顺序资料同化方法 [5] [7] [8] [9] 。顺序同化方法主要包括卡尔曼滤波(KF)、变形卡尔曼滤波、集合卡尔曼滤波(EnKF)和粒子滤波等。变分同化方法主要包括三维和四维变分资料同化方法(3D/4D-Var)，后者是前者在时间维上的扩展且同化效果更好。随着高性能计算机能力的提高和精细化数值预报需求增大，混合资料同化技术开始出现，例如：集合变分资料同化方法、集合最优插值同化方法、集合卡尔曼滤波和变分同化的混合方法，等等。虽然4D-Var是目前国际上气象和海洋业务预报中最先进和应用最成功的同化方法，但有效求解4D-Var问题时需要引入切线性和伴随模式，利用自动微分工具或手工编码方法仍然无法开发出完美无缺的伴随模式 [13] [14] [15] 。因此，众多的科研工作者一直都在探索和研究新资料同化技术 [2] - [12] 。

资料同化问题是典型的反问题，和正问题主要研究解的性质和数值求解方法等不同，反问题是通过试验或运行中的观测资料反求模型的未知参数：模式初始值、模型参数和模式误差等，从而使模型预测尽量准确或接近观测资料。因为观测量与未知参数之间常常不存在显式的直接关系，同时由于观测不准确、不充分和系统非线性等特征，所以导致反问题求解经常是不适定的。即解不一定存在、即使解存在也不唯一、在解存在唯一条件下也不稳定(即解不连续依赖于观测数据) [1] 。因此非线性动力系统资料同化问题的求解必须采用特殊方法 [4] - [9] 。本文在贝叶斯理论的基础上，提出基于马尔科夫链蒙特卡罗(Markov Chain Monte Carlo，简称MCMC)算法 [16] 来定量计算资料同化中初始值和模式误差的概率密度分布。综合利用贝叶斯方法和MCMC算法求解资料同化问题，具有以下优点：1) 能方便地将各种先验信息融合到资料同化问题的求解过程中，减小问题不确定性；2) 与确定性算法不同，反问题的不适定性不再是MCMC算法要考虑的问题，且计算获得的是全局最可能解，而在变分资料同化方法中，如果背景场准确性不够高，那么最优化算法可能陷入目标函数局部极小值；3) 能对定义在高维空间且复杂的概率分布密度函数进行数值计算，而确定性方法无法解决此类问题；4) MCMC算法通过构造Markov链在初始状态和模式误差所构成的参数空间进行重要性采样，最后获得的初始状态样本序列之间是平等关系，因此能提供给集合预报系统作为初始场集合；而变分资料同化极小化计算是一个逐步寻优过程，后面的迭代值要优于初始或中间值，即最终的分析场在理论上是唯一最优的。

本文首先利用贝叶斯公式推导了非线性动力系统需要估计的初始状态和模式误差分布规律的后验概率密度函数，参数估计值被认为是对应一维边缘后验分布的数学期望。接着采用自适应Metropolis算法 [17] 以后验概率分布为极限不变分布来构造Markov链，即对未知参数进行重要性抽样，并截取收敛后的样本序列计算数学期望，从而得到初始状态和模式误差的估计值；利用初始条件和模式误差样本集合定量计算了每个未知参数的一维边缘后验密度分布以及两个参数之间的二维边缘后验分布，后者定量描述了初始状态值和模式误差之间的相互关系。最后进行了计算机模拟，数值试验结果表明：贝叶斯MCMC方法能有效地估计非线性动力模型的初始条件，即具有较好的同化效果。

2. 贝叶斯理论

首先考虑如下非线性动力系统的资料同化问题：

$\{\begin{array}{l}\stackrel{\dot{}}{x}\left(t\right)=F\left[x\left(t\right),t\right]+w,t\in \left[0,T\right]\\ x|{}_{t=0}=x_0,\end{array}$

其中， $x={\left(x_1,x_2,\cdots ,x_n\right)}^{\text{T}}\in R^n$表示随时间变化的状态向量， $w={\left(w_1,w_2,\cdots ,w_n\right)}^{\text{T}}\in R^n$表示不随时间变化的常数型模式误差向量， $x_0$是初始状态向量，同时在 $\left[0,T\right]$时间段内分布有一系列观测量 $y\in R^m$。资料同化的目标是利用非线性动力预报模式(1)提取和融合观测信息，从而有效地估计由初始状态和模式误差参数构成的未知参数向量 $m={\left(x_0,w\right)}^{\text{T}}$。

贝叶斯推理方法的基本思想是：认为整体分布中未知参数 $m$本身就是随机向量，在对 $m$进行统计推断时，除了使用样本信息外，还需要对 $m$设置一个先验分布 $p\left(m\right)$。先验分布是在采集观测前就已存在的关于未知参数信息的概率表示。在获取观测资料以后，通过利用观测信息和贝叶斯公式可将未知参数的先验分布改进为后验分布。根据贝叶斯理论，参数的先验分布、观测样本信息和后验分布具有如下关系：

$p\left(m|y\right)=\frac{p\left(y|m\right)p\left(m\right)}{p\left(y\right)}$(2)

其中 $m^{\text{T}}=\left(m_1,m_2,\cdots ,m_N\right)$表示维数为N的未知随机向量， $y^{\text{T}}=\left(d_1,d_2,\cdots ,d_M\right)$是包含M个观测量的列向量。 $p\left(m\right)$表示未知随机向量的先验概率密度， $p\left(y|m\right)$是条件概率密度， $p\left(m|y\right)$是后验概率密度。因为观测量已经给出，所以 $p\left(y\right)$是一个与 $m$无关的常数，于是(2)式可以写成：

$p\left(m|y\right)\propto p\left(y|m\right)p\left(m\right)$(3)

(3)式是进行贝叶斯推理的基础，各项的物理或统计意义解释如下：先验概率分布 $p\left(m\right)$包含了观测获取前就存在的未知随机向量先验信息，它是在进行贝叶斯推理时的一种必要信息。从反问题的角度，先验信息的引入增加了信息容量、减小了不确定性范围，从而可以部分克服反问题不适定性。先验信息主要来源于历史观测资料、经验和主观判断等，例如：变分资料同化中的背景场及其误差协方差矩阵就是一种典型先验信息。条件概率密度 $p\left(y|m\right)$又称似然(Likelihood)函数，也可以写作 $L\left(y|m\right)$，包含了观测信息，即在有观测条件下未知参数的似然度信息。通俗地说，其表示了所估计的模式初始状态和模式误差与观测量之间的拟合程度，值越大表明拟合程度越好，反之越差。通过融合先验信息和观测信息后，就得到反映未知随机向量 $m$整体信息的后验概率密度函数 $p\left(m|y\right)$，它定义在初始状态和模式误差的整个解空间，表示问题的“完全”解。

在资料同化中，通过贝叶斯公式导出初始状态和模式误差等未知随机向量的后验分布表达式后，理论上可以获得未知参数的统计特征，例如：边缘分布、均值和方差等。但在实际应用中除了极简单情形，后验概率密度函数都无明确的数学表达式。另外，采用一般的数值积分方法(例如，Monte Carlo方法)也存在极大困难：必须对整个未知向量空间进行随机抽样以获得代表性样本，计算量将随向量维数增加而呈指数增长。上述原因使得直接利用贝叶斯方法无法解决复杂的实际问题，但是随着马尔可夫链蒙特卡洛算法的新发展，该情形不断被改善。本文结合贝叶斯理论和MCMC算法对资料同化问题中的初始状态和模式误差参数进行估计。

3. MCMC算法

通过上面的分析易知，虽然利用贝叶斯公式可以导出资料同化中初始状态和模式误差等未知参数的后验概率密度函数公式，但仍然无法求解资料同化问题。另外，考虑到观测误差和模式误差等不确定性因素的影响，单一“最可能”解存在巨大的局限性；而利用后验分布函数对解进行整体推断更具有合理性。换言之，研究重点不应该是通过最优化获得单一解，而是对解空间中最可能区域进行重要性抽样(importance sampling) [16] 后，然后基于参数样本集合计算解的估计值和置信区间。马尔可夫链蒙特卡洛算法是一种直接模拟后验概率密度函数的方法 [16] ，通过自动搜索概率大值区域而对未知向量进行随机抽样，然后由所得抽样序列对每个未知参数进行各种整体性推断。

MCMC (Markov chain Monte Carlo)方法是现代统计计算中最重要的算法之一，通过MCMC方法可以得到复杂模型中众多物理参数的有效范围。MCMC算法基于马尔科夫链的采样机制，可以对定义在高维随机向量空间 $M$上无明确数学表达式的概率分布 $\pi$进行有效抽样，基本思想是产生大量服从分布 $\pi$的随机向量序列 $\left\{m_{\text{1}},m_{\text{2}},\cdots ,m_I\right\}$，其中I为抽样数。如果向量序列满足马尔可夫性质：向量 $m_{i+1}$的产生仅依赖于前一个向量 $m_i$，而与过去时刻 $i-1,i-2,\cdots ,1$的状态向量 $m_{i-1},m_{i-2},\cdots ,m_1$都无关，则该向量序列称为马尔可夫链。马尔可夫性质的另一种表述是：若抽样算法当前访问的是 $m_j$点，则下一步访问另一点 $m_i$的概率只依赖于 $m_j$，而与先前访问的点无关。马尔科夫性质意味着抽样算法完全可由转移概率矩阵 $\pi$描述，矩阵元素 ${\pi }_{ij}$表示算法在当前访问 $m_j$的条件下接着将要访问 $m_i$的条件概率。按照构造Markov链所用转移概率矩阵的不同，MCMC方法的主要抽样算法有：Gibbs抽样器算法、Metropolis-Hastings算法和自适应Metropolis算法 [17] 。

本文采用自适应Metropolis算法以非线性动力预报系统初始状态和模式误差参数的后验分布为不变极限分布来构造Markov链。与传统的Metropolis-Hastings算法相比，自适应Metropolis算法不需要预先确定参数的推荐分布，而是由后验参数的协方差矩阵来估算参数分布 [17] ，从而大大提高计算效率。后验参数的协方差矩阵在每一次迭代后需要自适应地调整。如此，第i步参数的推荐分布就是均值为 ${\bar{m}}_i$、协方差矩阵为 $C_i$的多元正态分布 [17] 。协方差的计算公式如下面(4)式所示，在初始 $i_0$次迭代中，协方差矩阵 $C_i$取固定值 $C_0$，之后进行自适应更新。 ${\hat{m}}_i$是未知模式参数向量 $m$中某个元素在第i次迭代的估计值。

$C_i=\{\begin{array}{l}{C}_0,i\le i_0\\ s_{d}Cov\left({\hat{m}}_0,\cdots ,{\hat{m}}_{i-1}\right)+s_d\epsilon I_d,i\ge i_0\end{array}$(4)

其中， $\epsilon ={10}^{-6}$，其引入是为了确保 $C_i$不成为奇异矩阵； $s_d$是比例因子，依赖于未知随机向量空间的维数d，目的是保证接受率在一个合适范围内，在本文中 $s_d={\left(2.4\right)}^2/d$。 $I_d$为d维单位矩阵。当进行第 $i+1$次迭代时，由公式(4)可导出协方差的计算公式(5)。

$C_{i+1}=\frac{i-1}{i}{C}_i+\frac{s_d}{i}\left(i{\bar{m}}_{i-1}{\bar{m}}_{i-1}^{\text{T}}-\left(i+1\right){\bar{m}}_i{\bar{m}}_i^{\text{T}}+m_i{m}_i^{\text{T}}+\epsilon I_d\right)$(5)

其中， ${\bar{m}}_{i-1}$和 ${\bar{m}}_i$是前面 $i-1$次和i次迭代的参数均值。自适应Metropolis算法的计算流程 [17] 如下：

1) 设定 $i=0$，对所有未知参数变量进行初始化；

2) 随机量的生成和接受，构造Markov链：

a) 利用公式(4)计算协方差矩阵 $C_i$；

b) 产生服从高斯正态分布的推荐参数值 $m^{\ast }~N\left(m_i,C_i\right)$；

c)计算接受概率 $\alpha =\min \left\{1,\frac{p\left(d|m^{\ast }\right)p\left(m^{\ast }\right)}{p\left(d|m_i\right)p\left(m_i\right)}\right\}$；

d) 产生服从均匀分布的随机数 $u~U\left(0,1\right)$；

e) 若 $u<\alpha$，则接受 $m_{i+1}=m^{\ast }$，否则 $m_{i+1}=m_i$。

3) 重复上面的步骤(a)~(e)，直到产生预先指定数量的样本为止。

4. 资料同化问题求解

下面以典型的海洋生物种群动力学模型为例，说明利用贝叶斯MCMC方法对资料同化问题进行求解的有效性。该非线性动力系统由如下的常微分方程组表示：

$\{\begin{array}{l}\text{d}{N}_1/\text{d}t=a{N}_1-c{N}_1{N}_2+w_1,\\ \text{d}{N}_2/\text{d}t=b{N}_2+d{N}_1{N}_2-e{N}_2^2+w_2,\\ \left(N_1,N_2\right)|{}_{t=0}=\left(N_1\left(0\right),N_2\left(0\right)\right),\end{array}$(6)

系统(6)表示了海洋生物二种群动力学模型，主要模拟种群之间捕食与被捕食的相互作用，刻画的是两个种群数量密度 $N_1$、 $N_2$随时间t的演化，其中 $\left(a,b,c,d,e\right)$为已知的模型参数。数值试验中采用四阶龙格-库塔算法求解微分方程组，步长设置为h。观测量仿真过程如下：首先设置模式状态初值，然后任其演化至 $Th$时刻，从而得到系统在离散时间序列 $\left\{0,h,2h,\cdots ,Th\right\}$上的状态量序列 $\left(N_1\left(t_i\right),N_2\left(t_i\right)\right)$， $i=0,1,2,\cdots ,T$。状态量可直接作为观测量 $\left(y_1\left(t_i\right),y_2\left(t_i\right)\right)$，也可加入随机噪声得到更真实的模拟观测资料。本文资料同化的目标是在已知部分观测资料的条件下，采用贝叶斯MCMC方法估计非线性动力系统(6)的初始状态和模式误差参数。

下面利用第1部分中的贝叶斯公式推导系统(6)中未知参数的后验概率密度函数公式。由于系统(6)中影响模式状态量的四个参数 $\left(N_1\left(0\right),N_2\left(0\right),w_1,w_2\right)$都是未知的，即需要对多参数进行联合反演，且只使用部分观测资料。根据贝叶斯公式(3)，且不考虑分母表示的常数项，则未知参数的后验分布 $p\left(m|y\right)=p\left[N_1\left(0\right),N_2\left(0\right),w_1,w_2|y\right]$可以通过下式进行计算：

$\begin{array}{l}p\left[N_1\left(0\right),N_2\left(0\right),w_1,w_2|y\right]\\ =p\left[y|N_1\left(0\right),N_2\left(0\right),w_1,w_2\right]p\left[N_1\left(0\right),N_2\left(0\right),w_1,w_2\right]\end{array}$(7)

公式(7)已经假设 $y=M\left(x_0,w\right)+\omega$，其中 $M\left(x_0,w\right)$表示系统(6)的离散数值模式， $\omega$为包含观测误差的独立分布随机噪声。假定随机噪声 $\omega$服从均值为零、标准偏差为 ${\sigma }_o$的正态分布，即 $\omega ~N\left(0,{\sigma }_o^2\right)$。同时假设由数值模式M引进的离散等模拟误差包含在模式误差向量 $w$中。通过上述假设，可以得到以下形式的似然函数：

$L\left(y|m\right)=L\left(y|x_0,w\right)=\frac{1}{{\left(2\text{π}{\sigma }_o^2\right)}^{n/2}}\exp \left[-{\Vert y-M\left(x_0,w\right)\Vert }^2/\left(2{\sigma }_o^2\right)\right]$(8)

式中n表示观测数量， $\Vert \cdot \Vert$表示欧几里得范数。

在贝叶斯推理中，认为未知向量 $m$是随机向量，首先需要设定未知参数的先验分布。在本文中设定所有先验分布都是独立的均匀分布，初始状态和模式误差参数分别满足均匀分布： $U\left[N_1\left(0\right)\right]$、 $U\left[N_2\left(0\right)\right]$、 $U\left(w_1\right)$和 $U\left(w_2\right)$，则总的先验分布表示为：

$p\left(m\right)=U\left[N_1\left(0\right)\right]U\left[N_2\left(0\right)\right]U\left(w_1\right)U\left(w_2\right)$(9)

均匀分布是一种最简单的先验分布，虽然只能指定未知参数变化的上下界，但是可以缩小对参数随机抽样的目标区域，有利于提高参数的估计精度。一般来说，可以通过经验知识和历史统计确定更复杂和准确的先验分布。自适应Metropolis算法的一个优点是对于 $m$的任何先验分布都能够收敛于目标分布。由式(7)、(8)和(9)可得，在给定观测条件下资料同化问题未知参数的后验概率密度函数表示为：

$p\left(x_0,w|y\right)=\frac{1}{{\left(2\text{π}{\sigma }_o^2\right)}^{n/2}}\exp \left[-\frac{{\Vert y-M\left(x_0,w\right)\Vert }^2}{2{\sigma }_o^2}\right]U\left[N_1\left(0\right)\right]U\left[N_2\left(0\right)\right]U\left(w_1\right)U\left(w_2\right)$(10)

5. 数值试验结果

在数值试验中，首先设定海洋生物种群动力模式初始状态量 $\left(N_1\left(0\right),N_2\left(0\right)\right)$、模式误差参数和时间积分区间 $\left[0,4\right]$，取积分步长 $h=0.01$和利用四阶龙格-库塔算法对模式进行积分，就可以得到系统在离散时间序列 $0h,1h,2h,\cdots ,400h$上的标准状态值。试验中只抽取 $0h,20h,40h,\cdots ,400h$时刻的状态量作为 $n=21$个观测数据，同时在所选状态量上叠加高斯型观测噪声 $N\left(0,{\sigma }_o\right)$。噪声的均值和标准偏差分别取为0和 ${\sigma }_o$。未知参数 $N_1\left(0\right)$、 $N_2\left(0\right)$、 $w_1$和 $w_2$的先验分布分别取如下形式：

$U\left[N_1\left(0\right)\right]=\{\begin{array}{l}1/30-15<N_1\left(0\right)<15\\ 0otherwise\end{array}$，(11)

$U\left[N_2\left(0\right)\right]=\{\begin{array}{l}1/30-15<N_2\left(0\right)<15\\ 0otherwise\end{array}$，(12)

$U\left(w_1\right)=\{\begin{array}{l}1/20-10<w_1<10\\ 0otherwise\end{array}$，(13)

$U\left(w_2\right)=\{\begin{array}{l}1/20-10<w_2<10\\ 0otherwise\end{array}$。(14)

针对(10)式表示的后验概率密度函数，利用MCMC方法中的自适应Metropolis算法按照其计算流程中的步骤(a)~(e)构造每个未知参数的Markov链。未知参数初始值取对应先验均匀分布中的随机值，为了提高准确性，采用两次重要性采样，首先利用样本数为1000的Markov链获得一个比较准确的参数估计值。然后执行第二次重要性采样，进一步提高参数估计值。构造Markov链的过程实际上是在由先验分布界定的区间内对未知参数进行随机最优搜索。抽样过程中未知参数的每次更新，都要对非线性动力系统(6)式进行一次数值积分，以便计算后验概率密度函数的大小。图1~图3显示了初始状态和第一个模式误差参数的Markov链，可以看出经过第一次重要性采样后马尔科夫链基本上达到收敛；取Markov链中序号2000~5000之间的样本值计算数学期望，从而得到各个参数的估计值。

表1给出了在不同观测噪声水平条件下，新同化方法对海洋生物种群非线性动力系统(6)式初始状态和模式误差的估计值。从表中可知，在不存在观测误差时，未知初始状态的估计精度最高，初始状态量可以精确到小数点后第3位，表明了基于贝叶斯MCMC的资料同化方法估计非线性物理系统未知初始状态的有效性。另外，随着观测误差增加，资料同化结果的精度下降，但迭代估计结果仍然收敛到真实值附近；即使当观测误差的标准偏差为 ${\sigma }_o=0.2$时，虽然模式误差参数的估计误差较大，但是初始状态估计值仍然较接近真实值，可以精确到小数点后1位，说明基于贝叶斯MCMC的资料同化方法具有较强的抗噪声性能。图1、图2和图3分别给出了初始状态 $N_1\left(0\right)$、 $N_2\left(0\right)$和第一个模式误差参数 $w_1$的Markov链。从图中可以看出，经过第一阶段马尔科夫链1000次迭代后，第二阶段初始状态的马尔科夫链变化幅度较小，基本达到收敛。取2000~5000步的样本序列并采用后验均值方法计算参数估计值，分别得到 $N_1\left(0\right)=8.0010$和 $N_2\left(0\right)=1.0023$。分析图中所示的试验结果可得如下结论：采用贝叶斯MCMC方法估计非线性动力模式初始状态值具有较高的准确性和稳定性。对于模式误差参数，从图3中可知，相对于初始状态的Markov链，变化幅度增大，但基于2000~5000步Markov链值计算的后验均值非常接近真实值(见表1)。模式误差参数 $w_2$的Markov链情形与 $w_1$类似，受篇幅所限不再给出。采用贝叶斯MCMC方法估计出初始状态和模式误差值后，代入海洋生物种群非线性动力系统(6)的数值模式后可以计算出状态量 $\left(N_1,N_2\right)$的预报轨迹。图4和图5分别显示了 $N_1$和 $N_2$的预报轨迹、观测值和离散时间点上的真实值，从图中可知，预报轨迹与观测结果、真实值吻合得很好，验证了用贝叶斯MCMC同化方法对多个未知参数进行联合估计的正确性。

图6展示了在观测误差标准偏差取 ${\sigma }_o=0.1$时初始状态和模式误差等4个自由参数的边缘后验分布。从图中可以看出所有参数的一维后验分布(对角线位置)都是类高斯分布，参数样本分布在一定的置信区间，没有出现发散情形，而且频率曲线的最高峰位置对应横坐标上参数真实值，这种情况说明观测资料对初始状态和模式误差具有很强的限制和约束作用，也就是贝叶斯MCMC资料同化算法能很好地将观测信息传递给未知自由参数，在逐次重要性采样中不断调整未知参数值，最终收敛到真实值附近。一维后验分布给出了资料同化问题的“完整解”，即各个参数的概率分布数值，通过后验均值方法能计算出准确的初始状态和模式误差(如表1第三行所示)。图6下三角每个二维后验分布子图表示了不同参数之间的相关关系，从图中可知：任意两个参数之间都具有较大的相关性，即初始状态和模式误差存在显著的相关作用，说明在资料同化过程中模式误差对初始状态值有重大的影响；任意两个未知参数之间的相关形态都不相同，即两个参数之间的相互影响和作用不相同；显示的相关函数呈现出近似椭圆状，即是类二维高斯分布，说明观测信息也对二维后验分布具有很强的限制和约束作用。

6. 结束语

本文在贝叶斯理论框架下，提出了一种基于Bayesian MCMC算法 [16] [17] 从观测信息中同时估计非线性模型初始状态和模式误差等未知变量的新同化方法。主要结论如下：1) 贝叶斯MCMC方法在不需要切线性/伴随模式的情况下，能准确地估计出模式初始状态和模式误差参数，验证了新方法的有效性。2) 与确定性同化(例如：变分资料同化)方法只能获得单一最优解不同，贝叶斯MCMC同化方法不但能获得作为分析场的后验均值，而且能定量计算出估计参数的一/二维后验分布。其中一维后验分布给出了同化结果的分布范围和取值频率，而二维后验分布定量刻画了不同参数之间的相关关系。3) 新同化方法获得的收敛初始状态样本序列之间是平等关系，因此能提供给集合预报系统作为初始场集合。综上所述，与基于伴随模式的变分资料同化等确定性方法 [6] [7] [8] [9] 一样，Bayesian MCMC方法具有求解资料同化问题的能力：能准确地估计出初始状态和模式误差，同时具有较好的抗噪声性能。

基金项目

国家自然科学基金资助项目(41475094; 41105063; 41375105)。

参考文献