数据同化（Data Assimilation ， DA），又称资料同化，是一种源自于大气领域的天气预报方法。它的主要思想是在模型动力框架内，将观测信息不断融合到动力系统中，从而自动调整模型轨迹改进系统状态估计，并且减小误差的预报系统。数据同化的基础是估计理论和控制论^［1］。传统的数据同化方法主要分为变分类方法和Kalman滤波类方法。在实际应用中，Kalman滤波存在特定的条件限制，一般存在假定的条件与实际状态估计之间的不匹配^［2］，如通常假定模型误差和观测误差为零均值高斯白噪声，在实际水文模型和陆面模型中，误差统计特性大都不满足该特性；同时量化和测量模型误差的统计特性难以确定^{［3］［4］}。针对Kalman滤波对模型误差和观测误差的特殊要求，Zames等^［5］提出了H∞滤波，该方法隶属于鲁棒类滤波，能够保证在外界误差条件较差的情况下的估计性能，如模型模拟中存在严重人为扰动，引入未知的误差，人为的输入造成模型预测的非平稳和不可预测的有偏误差。在鲁棒滤波中将H∞范数作为目标函数优化设计，目标函数中包含模型的初始条件，量测噪声和模型噪声，使得鲁棒滤波对初始条件的不确定性不敏感且能有效纠正系统偏差，保证误差估计有界。Wang和Cai^［6］系统比较了Kalman滤波与H∞滤波（H-infinity Filter，HF）的优劣，当外界条件变化不可预测时，HF能够处理模型模拟中的未知误差，减少误差对估计性能的影响且HF不需要对误差的统计特性做相应的假设，能够保证其估计特性性能处于一定的有效范围内。为了应用于顺序数据同化和解决高维同化系统问题，Luo和Hoteit^［7］将集合的思想应用于时间局地化的H∞滤波（Time-Local H-infinite Filter，TLHF）。Triantafyllou等^［8］将基于集合鲁棒卡尔曼滤波应用到生态资料同化中，提高了滤波的精度及鲁棒性。从观测角度优化数据同化方法，摆玉龙等^［9］提出放大观测误差协方差矩阵的集合鲁棒滤波，利用分析误差协方差矩阵的逆阵为半正定的条件，放大转移矩阵的特征值，间接地放大分析误差协方差矩阵，避免直接对分析误差协方差矩阵进行复杂的奇异值分解（Singular Value Decomposition， SVD），降低了计算复杂度，提高滤波效率^［10］。然而，上述研究中的鲁棒类滤波方法都将观测误差假定为单位阵且观测误差不具备相关性，显然不能满足实际同化系统中的需要。

观测误差协方差在数据同化方法中起着重要作用，它决定赋予模型预测和解决方案中观测值的权重，观测误差矩阵会影响分析值的准确性。数据同化中，为了减少计算复杂度，根据经验通常将观测误差协方差定义为单位对角阵，然而研究表明，在一些特定情况下，观测误差是相关的^{［11，12］}且与时间和状态存在依赖性^［13-16］。观测误差通常包括两部分，仪器误差和代表性误差。其中，仪器误差可以根据仪器的使用来校准^［17］。代表性误差的处理一直以来是数据同化研究中心的难点。根据2014年欧洲航天局研讨会有关代表性误差术语的讨论，Janjic等^［18］给出代表性误差更为统一的定义，将过程中无法解决的尺度问题，模型正演误差或是观测算子带来的误差和预处理数据的误差统一定义为代表性误差。Janjic和Cohn^［13］表示由于尺度不匹配造成的观测误差是相关的且与时间状态有依赖性。同时，Swewart^［11］和Desroziers等^［19］认为，在同化中考虑这些相关的观测误差时，可以得到更准确的分析值^［20］。Healy等^［21］也证明粗略的估计观测误差协方差矩阵能提高估计性能。

在观测误差估计与处理方面，2005年， Desroziers等^［19］介绍了一种用来估计观测误差协方差的方法，称为DBCP诊断，该方法主要利用观测值与背景值和分析值之间的偏差（新息统计）来估计观测误差协方差。2009年，Anderson^［22］和Li等^［25］利用新息统计方法估计自适应预测协方差矩阵的膨胀系数，来改善预测协方差低估等问题，Bocquet等^［23］将该方法用于迭代集合Kalman滤波，Gharamti^［24］使用该方法来改善自适应滤波，提高滤波精度，Li等^［25］和Miyoshi等^［26］将该方法应用于集合转换卡尔曼滤波得到观测误差协方差的估计，在Kalman滤波中，观测误差协方差的粗略估计能够提高状态估计的精度^{［11，21，26］}。因此，文中考虑将观测误差协方差估计方法与鲁棒滤波方法相结合，使得观测误差协方差与最近的时间和状态相关，来提高状态滤波的精度。

根据观测误差协方差对时间和状态间的依赖性，使用观测误差估计技术来代替经验假定的对角单位阵。利用前一时刻观测值分别与背景值和分析值的偏差来估计观测误差协方差，得到新的观测误差协方差与上一时刻的背景状态值相关，新的观测误差协方差将用于下一次同化。在实验中，观测误差协方差随着同化进程更新，使用估计得到的观测误差协方差来代替单位对角阵，减小观测误差带来的干扰，提高滤波精度。

数据同化方法将某时刻的观测值yi∈RNp与模型状态xnf∈RNm预测值相结合完成预报过程，其中Np和Nm分别是观测矩阵和模型状态变量的维数。通过对观测值与背景值取不同的权重值来提供精确地分析值xia，得到的状态分析值通过模型演进成为下一同化时刻的背景值。以下为预报方程和观测方程：

其中：xi+1∈RNm表示第i+1时刻m维状态变量的预测值；Mi+1,i表示转移算子，将第i时刻的m维状态变量映射到第i+1时刻；Hi是观测算子，将m维的状态变量映射到p维的观测空间；为了方便计算，假定系统是线性的，则Mi+1,i和Hi是两个矩阵，ui和vi表示预报误差向量和观测误差向量。

Kalman滤波是将当前观测值与背景值（上一时刻的预测值）结合得到当前时刻预测值的一种方法。传统的Kalman滤波方法假定预先知道模型和观测误差的统计特性，已知模型结构、参数和观测中涉及的不确定性，对观测误差进行量化。为了改进这种限制，Luo和Hoteit^［7］将HF中的目标函数用于Kalman滤波中，提出一种新的滤波方法时间局地化滤波。

在HF中^［7］存在目标函数满足下面的不等式：

其中： 上标a和b分别表示分析和背景变量，xi-xib(Δib)-12表示背景误差能量，||ui||Qi-12，||vi||Ri-12分别表示模型误差，观测误差的能量。Δib表示背景协方差矩阵，λi是第i次循环的局地化性能水平，Qi和Ri分别为系统噪声协方差观测误差协方差，Si表示权重矩阵，由设计者选择^［7］。式（5）是最小最大准则，inf表示目标函数的下界，sup表示目标函数的上界^［7］。λ*保证存在一个λ值满足最小最大准则^{［6］［26］}。公式（4）中，当λi为零时，公式右边趋于无穷，所以估计误差的能量小于无穷，说明该算法不能保证估计误差能量有界，鲁棒性不强；在λi不为零时，此时为鲁棒滤波，能够保证滤波在受到扰动时误差增长能量有界，说明鲁棒滤波鲁棒性较好。

在线性系统中，TLHF的滤波过程如下：

预报：

其中：xib是i时刻的背景值，Δib表示背景误差协方差矩阵。

滤波：

其中：符号Δi表示不确定矩阵，类似于Kalman滤波中的协方差矩阵Pi，Gi是增益矩阵，在Kalman滤波中用Ki表示^［7］，Hi表示线性观测算子，xia和Δia分别为i时刻的分析值和分析误差协方差矩阵。类比于HF，存在一个目标限制函数：

其中：满足式（11）表示矩阵是半正定矩阵，其取值与Ri，λi，Si等变量取值有关，在公式（5）中λi的取值依赖于初始条件和观测系统^［7］，目标函数Jx,iHF对应变量xi，ui和vi，所以λi的取值的与Ri，Qi，Si有关， Ri，Qi，Si通常由设计者自由选定^［27］。

Desroziers等^［19］提出检验状态空间的各种新息值，随后这种估计方法应用于数据同化领域。该方法主要基于观测、预报、分析之间的不同新息值统计，分别叫做背景新息值：db=y-H(xf)，分析新息值：da=y-H(xa)。当系统满足线性和高斯条件时，此时系统进行无偏预测和观测，对背景新息值和分析新息值求期望得到观测误差协方差估计值，公式如下：

数据同化方法中，观测误差协方差主要分为两部分：一是与时间状态无关的观测误差协方差，主要指仪器误差，用εnI表示，对应的方差为δD2，得到RI=δD2εnI；另一部分与时间状态相关的观测误差协方差，主要由观测算子引起的误差和代表性误差组成，用εH表示，对应的方差为δC2，得到RH=δC2εH，所以此时观测误差公式为^［28］：

实验中，使用SOAR（second-order-autoregressive correlation function）函数构造RH，计算公式为：

其中：r表示两个观测点之间的距离，L表示相关函数的长度标量，通常取值为6^［28］。

估计得到的观测误差协方差用于下一步滤波过程中的Kalman增益，每次同化时Kalman增益值随时间变化，每完成同化一次，观测误差协方差更新一次，迭代计算观测误差协方差函数。在目标函数中使用最新观测误差协方差，对目标函数进行判断，使其分析误差协方差满足条件限制，即为正定矩阵。同时得到满足该情况下的λi取值。观测误差的取值会影响λi的取值，能得到较为精确的λi的取值范围。

类似于EnKF方法，将集合的思想运用于EnTLHF（Ensemble Time-Local H-infinity， EnTLHF），提高滤波的鲁棒性和滤波精度^［7］。xib={xi,jb:xi,jb=Mi,i-1(xi-1a,j),j=1,⋯n}是i时刻的n个背景集合成员，是第i-1时刻的分析集合xi-1a=xi-1,ja,j=1⋯n的预报，转移算子Mi,i-1是非线性的。

EnTLHF的预报过程：

其中：x¯ib和Δib分别是背景和估计所对应的不确定矩阵。

滤波过程：

目标限制：

其中：Im表示m维的单位阵，其中Xib表示预报集合协方差矩阵平方根，x¯ia分析状态值。

在EnTLHF滤波中，观测误差通常由设计者定义，设计者根据经验将观测误差假定为对角阵，对角阵满足距离越远，观测误差的相关性越小，这种假设在观测数较小时有效^［17］。然而最新研究表明，观测误差协方差与状态变量相关且与时间变量有一定的依赖性，观测误差的估计与假定会影响预测误差协方差的大小，进而影响观测值在同化时的权重，对同化的精确性和同化的分析值产生影响。所以研究中考虑观测误差协方差估计技术，将观测误差估计技术与鲁棒滤波技术结合研究新的鲁棒滤波方法的鲁棒性。

将观测误差估计技术与鲁棒滤波方法相结合，得到带有观测误差协方差估计的鲁棒滤波方法（Ensemble Time-Local H-infinity with observation error covariance estimation， EnTLHFR）。实验中首先对集合数，背景误差协方差矩阵和观测误差协方差矩阵进行初始化，定义观测误差估计的样本数为同化步数的一半。然后进行EnTLHF算法预报和滤波，每完成一次滤波，使用背景新息值和分析新息值估计该时刻的观测误差协方差，更新的观测误差协方差将用于下一时刻的滤波步骤，代替上一时刻的观测误差协方差。依次迭代计算Ri，每一次观测误差协方差都由上一时刻的分析值和背景值计算得到，大小与上一时刻的状态值相关，舍弃了无关的状态变量仅保留最近时刻的相关误差，能够得到更加准确的估计^［28］。相关算法如表1所示。

公式（22）中，当λ>0时，-λiSi<0，所以该方法放大估计误差协方差。下面通过EnTLHF的放大项进行几种特殊的变换，得到与已经存在的协方差放大理论相似的理论关系，有如下情形：

①当λiSi=c[(Δib)-1+(Hi)T(Ri)-1Hi]=c(pia)-1,0<c<1时：

这种条件下的放大称为分析协方差方法，标记为EnTLHFR-IANA^［6］，c称为性能水平系数（Performance Level Coefficient，PLC）^{［7］［27］} ，此时λi对滤波效果的影响由性能水平系数c决定，所以在实验中考虑性能水平系数c的变化对滤波效果的影响。

②当λiSi=c(Hi)T(Ri)-1Hi,0<c<1时：

这种方法是基于观测构建的协方差放大法，等同于协方差放大技术，标记为EnTLHFR-IR。

③当Si=Im时：

在ETKF（Ensemble Transform Kalman Filter，ETKF）中，对分析误差协方差矩阵进行特征值分解，得到：

Pia=(XibCi)-T(Γi+I)(XibCi)-1，Ci和Γ表示对矩阵(HXb)TR-1HXb特征值分解对应的特征向量和特征值。权重矩阵Γ=diag(σi,1,σi,2,...,σi,n-1)，是包含特征值σi,j的对角阵，当j<l时，σi,j≥σi,l，为了满足目标函数的限制，即满足σi,j-(λi-1)≥0，令(λi-1)=cσi,n-1,0<c<1，得到：

将这种放大称为转移矩阵放大，标记为EnTLHFR-IT，类似于Luo和Hoteit（2001）介绍的放大转移矩阵的特征值^{［7］［9］}。实验中，每同化一次，估计观测误差协方差得到更新，更新后的观测误差协方差用于下一时刻的滤波，同时判断目标函数是否为正定矩阵，得到满足条件的性能水平系数取值。

Lorenz-96模型是由微分方程表示的二阶非线性动力系统，即：

其中xNx=x0,x-1=xNx-1,x-2=xNx-2（Nx=40），F是驱动参数，时间间隔为0.05 s。利用经典四阶Runge-Kutta得到Lorenz-96方程的解。实验中迭代总步长为1 500步，为了避免暂态影响，舍去前500步进行同化实验。其中每4个时间步长进行观测，观测系统用yi=Hi(xi)+vi表示，其中vi服从N(0,δD2I)。

实验中采用RMSE检验同化方法误差性能，公式为：

其中：m是系统维数，Xt是真实状态值，||•||2是欧几里得范数。

实验中为了更加准确评价滤波性能指标，引进另一个性能指标，即同化随时间变化的平均集合散度（Average Ensemble Spread，AES）^［29］：

其中：σi,k2是xi,k的方差。

性能水平系数的取值与协方差放大的效果相关，同时性能水平系数的大小反映性能水平γi对滤波效果的影响，性能水平系数对目标限制函数存在约束条件，不同的性能水平系数同化效果不同。图1表示带有观测误差估计的鲁棒滤波方法的RMSE随着性能水平系数的变化图。设置参数：F=5，集合数为30，初始观测误差为RT。

结果显示：①在分析协方差放大条件下，带有观测误差估计的鲁棒滤波性能水平系数最大值为0.5，当取值为0.5时，估计误差均值出现极大值，说明此时滤波不稳定；②当大于0.5时，滤波发散。滤波中，考虑观测误差的相关性并进行其值估计，能够限制目标函数进而影响性能水平系数的取值。

图2表示在性能水平系数c=0.3，集合数为20，初始观测误差协方差为单位阵和RT，强迫参数F分别为6，8，9时，带有观测误差估计的鲁棒滤波和鲁棒滤波的均方根误差均值（RMSE）变化图。

从图中可以看出：①当初始观测误差为单位阵时，即在同化中只考虑仪器误差，随着同化步长的演进，鲁棒滤波方法的估计误差没有明显的减小。且随着强迫参数的的增加，即模型模拟误差的增大，鲁棒滤波估计误差逐渐增大，带有观测误差估计的鲁棒滤波方法的估计误差略微小于鲁棒滤波的估计误差，此时带有观测误差估计的鲁棒滤波方法的精确性效果不明显。②强迫参数较大时，模型非线性程度较强时，此时两种滤波方法波动较大，说明两种方法稳定性相差不大。③当初始观测误差为RT时，即考虑仪器误差和代表性误差，在整个强迫参数变化时，带有观测误差估计的鲁棒滤波估计误差均小于鲁棒滤波，比较于原鲁棒滤波分析精度更高。总体说明带有观测误差估计的鲁棒滤波方法对不同程度的模型误差具有较好的鲁棒性，证明观测误差协方差的估计可以改善滤波性能。

观测数目在同化方法扮演非常重要的作用，有效合理的观测数目能够使估计值更加接近真值，观测数目越多，越能表现当前状态的真实性。图3和图4分别表示性能水平系数c=0.3，集合数为20，强迫参数F=6，观测误差协方差为0.1倍单位矩阵和观测误差协方差为RT时，改变同化间隔步长和状态量观测数目的条件下（分为4种情况：每一个状态变量都有观测，每一步都进行同化；每隔一个状态变量有观测，每一步都进行同化；每隔一个状态变量有观测，每隔3步进行同化；每个状态变量都有观测，每隔3步进行同化） EnTLHFR-IR和EnTLHF-IR两种方法的AES（Average Ensemble Spread）随着同化时间的变化。

从图3和图4可以看出：①在观测状态数和同化步数变化的4种情况下，当观测误差协方差为单位阵和RT时，EnTLHFR-IR的AES在不同观测条件下都小于EnTLHF-IR的值，说明带有观测误差协方差估计的协方差放大鲁棒滤波的效果更佳；②在观测误差为单位矩阵和RT时，图（a）和图（b）表示每一步都进行同化，图（a）每个状态变量都有观测数，而图（b）只有一半的状态变量有观测数，从图中看出，每个状态变量都有观测的AES小于只有一半状态变量有观测数的AES，带有观测误差协方差估计的鲁棒滤波的AES小于原鲁棒滤波方法的AES，说明带有观测误差估计的鲁棒滤波对观测数目变化具有较好的鲁棒性。③图中（a）和（d），（b）和（c）的观测状态数一样，同化步长分别为1和3，随着同化间隔步数增加鲁棒滤波方法的AES增大，当同化间隔步长增加时，即先验信息的缺失，带有观测误差估计的鲁棒滤波方法的AES波动小于原鲁棒滤波，说明此该滤波收敛程度较优。④图（c）表示只有一半的状态变量有观测数且每隔3步同化一次，这种情况下，带有观测误差估计的鲁棒滤波方法与原鲁棒滤波方法的AES最大。总体来讲，EnTLHFR-IR的AES值偏小，说明在状态变量有无观测及同化间隔步长变化时，带有观测误差估计的鲁棒滤波方法的收敛程度小于原鲁棒滤波方法。

图5表示矩阵放大实验。实验中，观测误差协方差为0.1倍单位矩阵，集合数为20，强迫参数分别为6、8、9，性能水平系数不同时，带有观测误差协方差估计的鲁棒滤波方法与放大矩阵鲁棒滤波方法的RMSE随性能水平系数变化图。

结果显示：①当强迫参数分别为6、8、9，c不为0时，采用放大矩阵的鲁棒滤波的RMSE小于当c等于零的滤波方法（ETKF（Ensemble Transform Kalman Filter）方法），说明鲁棒滤波方法对模型误差具有更优的鲁棒性。②在性能水平系数和强迫参数相同时，带有观测误差估计的鲁棒滤波方法的RMSE值均小于原鲁棒滤波方法。③带有观测误差估计的鲁棒滤波方法在性能水平系数大于0.6时，此时该滤波方法的RMSE会出现突增，说明在该种情况下，即当性能水平系数大于0.6时，该鲁棒滤波方法不收敛，会出现滤波发散。从图（d）中可以看出，当性能水平系数相同时，随着强迫参数F增大，带有观测误差估计的鲁棒滤波方法的估计误差也会增大，当F取值较大时，此时估计误差偏大，进一步验证了带有观测误差估计鲁棒滤波方法放大转移矩阵的有效性和鲁棒性。说明带有观测误差估计的鲁棒滤波方法估计鲁棒稳健性更强。

将观测误差协方差估计技术与鲁棒滤波技术相结合，分析几种不同角度下带有观测误差估计的鲁棒滤波理论。新方法中，使用观测误差估计技术估计观测误差协方差，在滤波步骤中，使用估计得到的观测误差代替固定的单位对角阵，估计得到的观测误差协方差随状态量的变化而变化。同时，使用估计技术得到的观测误差协方差与最近的时间相关，提高状态分析值的精确性；减小估计误差，提高滤波精度，鲁棒性更好。

利用非线性Lorenz-96系统，将带有观测误差估计技术的鲁棒滤波与鲁棒滤波方法进行比较。首先，在分析放大实验中，通过固定性能水平系数和集合数，改变强迫项参数，考察系统模型对两种方法的鲁棒稳健性影响。结果表明，虽然两种方法的均方根误差随着强迫参数的增大而增大，但带有观测误差协方差估计技术的鲁棒滤波的均方根值始终略小于鲁棒滤波。当强迫参数增大，观测误差为RT，即考虑相关观测误差时，带有观测误差估计的鲁棒滤波的均方根明显小于鲁棒滤波，说明带有观测误差估计的鲁棒滤波对模型误差稳健性更强，提高了滤波精度。其次，通过观测角度来探究两种方法，在其他参数相同的情况下，改变观测数和同化步长来探究两种方法的收敛性，实验结果表明，4种情况下带有观测误差估计的鲁棒滤波的AES均小于鲁棒滤波，并且在整个同化过程中波动较小，即观测误差估计技术降低了误差扰动，能够提高滤波精度，新方法鲁棒性能较好。最后，利用转移矩阵放大法来观察新方法的鲁棒性能，实验结果表明，当性能水平系数不同时，新方法的均方根均小于鲁棒滤波，新方法中性能水平系数的最大值为0.6，当大于0.6时，会发生滤波发散，新方法的应用提高了状态估计的准确性，改善了滤波性能，稳健性更好。

通过数值实验验证了带有观测误差估计的鲁棒滤波的性能，考虑观测误差的相关性和对时间状态的依赖性，使用估计得到观测误差协方差代替对角阵，得到的新方法鲁棒性更好。后续工作将重点研究估计观测误差协方差的样本数的选择。