李宏,许建平,2
(1.国家海洋局第二海洋研究所,浙江 杭州 310012;2. 卫星海洋环境动力学国家重点实验室,浙江 杭州 310012)
资料同化技术的发展及其在海洋科学中的应用
李宏1,许建平1,2
(1.国家海洋局第二海洋研究所,浙江 杭州 310012;2. 卫星海洋环境动力学国家重点实验室,浙江 杭州 310012)
回顾了资料同化技术,特别是基于最优控制和统计估计这两大理论基础发展起来的几种资料同化方法的研究进展,以及这些方法在海洋科学研究中的应用现状。可以看到,由于海洋观测资料(如地转海洋学实时观测阵(Array for Real-time Geostrophic Oceanography,Argo)、热带大气海洋阵列(Tropical Atmosphere/Ocean array, TAO)、抛弃式温度深度仪(Expendable Bathythermograph,XBT)、温盐深仪(Conductance-Temperature-Depth, CTD)、海洋表层温度(Sea Surface Temperaturem, SST)、海面高度等(Sea Surface Height,SSH),存在种类繁多,且量大、均一性差等特点,资料同化技术面临种种挑战;同时,随着海洋科学的发展,无论是海洋数值预报精度的提高、海洋数据产品的开发,还是对物理海洋现象的细致描述和深化认识,也都离不开资料同化技术的进一步发展。
资料同化技术;发展;海洋科学;应用
资料同化(Data Assimilation,也叫数据同化)最初来源于为数值天气预报提供必要的初值,现在已经发展成为能够有效利用大量多源非常规资料的一种新颖技术手段,它不仅可以为海洋数值预报模式提供初始场,还可以构造海洋再分析资料集,为海洋观测计划和数值预报模式物理量及参数等提供设计依据[1]。近十年来,资料同化技术取得了快速的发展,从早期比较简单的客观分析法(Objective Analysis,OA)发展到现在能够同化大量非常规资料[2]的四维变分(4 Dimensional Variation,4D-Var)和集合卡尔曼滤波(Ensemble Kalman Filter,EnKF)等。如今,随着海洋观测技术和方法的不断改进,观测资料(如卫星遥感、雷达、船舶报、验潮站、XBT、TAO、Argo等)在大量递增。但这些资料存在空间分布不均匀(南半球资料偏少,北半球资料多,近海偏多,深海大洋偏少)、时间分布不一致、资料来源的不同性,以及各种误差信息的不统一性等特点。因此,如何将这些不同类别的资料合理地融合进入海洋数值模式,并能客观地揭示海洋中存在的各种物理海洋现象及其形成机制,从而能有效改进数值预报的结果,是资料同化技术面临的问题。
本文首先回顾了资料同化技术的发展过程,并对不同方法的优缺点做了简要阐述,随后重点介绍了资料同化技术在海洋科学领域中的应用,并展望了其未来发展的方向。
早期的资料同化方法,也称客观分析法(Objective Analysis,OA),如1949年由Panosky开创性提出的多项式拟合(Polynomial Fitting,PF)[3]、1954年Gilchrist和Cressman发展的函数滤波法(Function Filtering,FF)等[4],以及最初由 Gilchrist提出的理想逐步订正法(Successive Correction,SC)原型[4],后(1955年)由 Bergthorson对其进行理论论证[5]并由 Cressman(1959年)发展成熟的基于迭代算法的逐步订正法[6]等,这些其实都是经验分析方法,没有充分利用模式和观测资料的误差统计信息,也没有利用模式的时空演变信息,并且缺乏强有力的理论基础。因此,在实际数值预报特别是在海洋科学研究中并没有得到广泛应用。直到20世纪60年代初,最优插值法(Optimal Interpolation,OI)的提出,资料同化方法才有了基于统计估计理论的基础。目前的资料同化方法根据其理论原理可分为两类,一类是基于统计估计理论的,如最优插值、卡尔曼滤波(Kalman Filter,KF)、扩展卡尔曼滤波(Extend Kalman Filter,EKF)、集合卡尔曼滤波(EnKF)等;另外一类是基于最优控制的(也有称变分的),如三维变分(3 Dimensional Variation,3D-Var)、四维变分(4D-Var)等。本文主要回顾基于最优控制和统计估计这两大理论基础发展起来的资料同化技术。
1.1.1 OI OI[8]首先由Gandin在1963年提出,它是一种新的基于统计理论基础的均方差最小线性插值法,该方法考虑了模式和观测数据的误差统计信息,并加入必要的权重,能够较好的刻画出实际大气(海洋)状态,因而在 20世纪八、九十年代业务化数值预报当中占主流地位。
OI的一个基本假定是:在确定每个模式变量的分析增量时,仅有几个观测值是重要的。基于这一假设,OI就易于编码并且计算量相对较小,这是它的主要优点[9]。但OI的缺点也是显而易见的,由于所用协方差矩阵是固定的,不随时间变化,这就限制了它不能将动力模式和观测信息很好地融合在一起;且一般的 OI是单变量分析,这可能造成物理量的不协调[10]。OI通常选择分析格点附近的观测资料来做局部分析,这可以减小计算量,但分析结果并非全局最优,分析在空间上不协调[11];OI是针对线性系统发展起来的,难以处理观测算子非线性的情况;当对模式状态的不同部分采取不同的观测值时,会使分析场产生虚假的噪音[9];另外,OI无法确保大小尺度分析的一致性[12]。
1.1.2 KF KF算法最初由Kalman 1960年[13]引进用于离散时间下的线性系统滤波,且由 Kalman和Bucy[14]扩展用于连续时间系统滤波,类似于OI,KF也是基于统计估计理论发展起来的。在系统是线性、误差是白噪音和高斯型的情况下,KF以分析误差的最小方差为标准提供分析最优值。区别于后面的4D-Var,KF显式发展背景场误差协方差,因而不需要伴随算子,这是其一大优点;另外,KF可以直接提供分析误差协方差矩阵,这是 4D-Var不具备的优势。但由于其高昂的计算代价而难以应用于实际数据同化当中。
实际上,大气和海洋预报模式大多数是高维非线性系统,KF算法对此无能为力。于是,有学者针对非线性系统提出了被后人所称的“EKF”。
1.1.3 EKF EKF基于切线性假设(仅保留模式一阶导数项),对一般的弱非线性问题是一种很好的近似,但对于强非线性问题,这一假设本身就偏离实际很远,简化后的方程恰恰去掉了原始方程中最重要的部分,仅保留二三阶导数项的“闭合”技术,对非线性模式本身是一个很困难的问题,处理不好直接导致滤波发散[15]。因而,EKF难以应用于强非线性系统,同时,与EK一样计算量巨大。所以,EKF在实际应用中也是难以发挥作用。
为了解决(E)KF的计算代价问题,许多学者做了大量探索性工作[16-20],试图寻求一种次优方案来代替EKF,具有代表性的有降秩平方根卡尔曼滤波(Reduced Rank Square Root,RRSQRT)[17]和奇异进化扩展卡尔曼滤波(Singular Evolution Extended Kalman Filter,SEEK)[18],这两种方法均是通过对误差协方差矩阵进行特征值分解,并将其投影于主特征向量空间,由此来降低协方差矩阵的秩,以此来降低计算量和存储量。在所有的改进方案中,以Evensen 在二十世纪九十年代提出的基于集合思想的EnKF[16]最引人注目。
1.1.4 EnKF EnKF是一种基于蒙特卡罗算法的集合方法,它用有限的集合样本来估算误差协方差矩阵的不确定性[16,19-20],这样计算量明显减小。在系统为线性,且当样本数量趋于无穷时,EnKF和KF是等价的[21]。EnKF利用集合扰动的方法构造初始场易于编码实施,且用集合方法来估计背景场的误差协方差,这相对4D-Var来说,就不需要切线性假设和伴随算子的构造,也无需模式的反向积分[20-21],这是其主要优势;另外,EnKF自动提供分析误差协方差矩阵,这也是变分同化法所不具备的优点。
当然,EnKF也还存在诸多问题。首先卡尔曼滤波是基于误差的无偏估计及概论密度的高斯分布假设,这两个假设在实际当中并不一定是可靠的[9]。其次,由于EnKF是通过选取有限的样本来构造背景场误差协方差,这势必使得样本集合离散度不够(样本量有限),产生样本误差问题[22-25]。并且在实际操作中,系统的非线性性,及通常利用扰动观测法获取样本初值[24]使得这种样本误差问题更为明显。随后发展的平方根法(或者确定性方法)来获取样本初值[24-25]能有效避免传统样本误差问题,另外局地化方法能抑制由有限集合样本数引起的虚假远程观测相关性[26-27],但它也能导致初始条件的动力不平衡性[2,28]。再者,如 Houtekamer等[29]指出的那样,EnKF仅用常规意义下的一组集合实施同化会产生“近交(inbreeding)”问题,提出用两组集合实施同化:根据每组集合进行短期预报得出的协方差矩阵来计算的Kalman权重函数互相交换使用,这种方案可以避免“近交”问题,后被广泛采用[29-30]。
EnKF 的计算量相对4D-Var来说,依然很大,但其集合思想为资料同化技术的发展提供了一种新的思路,一些学者已将集合吸收到其他较为简单的同化技术当中,如Envesen等[20]最初提到集合最优插值(Ensemble Optimal Interpolation,EnOI)方法,就是EnKF的简化版。与EnKF一样,EnOI采用集合样本的方法来生成背景协方差矩阵,但EnOI用静态样本,只用一个积分模式,而EnKF用动态样本,需要多个积分模式,因此EnOI计算量远小于EnKF[20]。目前,EnOI也取得了较快的发展,感兴趣的读者可以参考相关文献[31-32]。
1.2.1 3D-Var 最初出现的变分同化技术是3D-Var,3D-Var基于最优控制理论而发展,通过分析预报值与观测值之间的距离最小化来得到海洋或大气状态的最优估计量。
相对于 OI,3D-Var可以做多变量同化分析,且在三维空间中进行全局分析,分析解为全局最优;也可以处理观测算子非线性的情况,这样可以同化各种不同来源的观测资料,包括常规的和非常规的、非同步的等,如 XBT、TAO、Argo[33]以及卫星观测资料[34]等;另外可以在代价函数上加入额外的平衡约束项,这样能抑制分析场带来的重力噪音[35]。
然而,由于三维变分无时间变量,因此动力模式不能对其进行约束,其获得的初值在时间上是不连续的,也难以保障与模式协调[9];另外,模式在同化时间窗口内被认为是静止的,而且使时间窗口内的任何观察数据都被认为是同一时刻的观测值,这无疑会使得同化结果与实际产生某些偏差[36],这是它的主要缺点。为弥补 3D-Var的这一缺陷,F.X.LeDimet等人于 20世纪 80年代提出了4D-Var[37]。
1.2.2 4D-Var 4D-Var[37]是将 3D-Var进一步扩展成为包含时间变量的同化分析。4D-Var在时间窗口内利用完整的动力模式作为强约束,能自动调整模式误差,使得同化结果更可靠。而且,规定在某一时间段上的观测数据均可纳入到同化系统。背景场误差协方差隐式发展,误差信息随动力模式而向前传播,这些是4D-Var的主要优势。由于4D-Var需要求解伴随模式,并且代价函数求解通常采用最速下降法、共轭梯度法及准—牛顿迭代法等迭代计算,所以,计算量特别大。
针对4D-Var计算量大的缺点,许多学者提出了新的改进方法,如早期由 Courtier提出的增量法[38],利用转换算子将原来模式高分辨率的增量场转换为低分辨率场,在低维空间中进行计算,最后利用逆转换算子将迭代获得的低维增量转换到原来的空间增量中,但增量法无法确保结果的收敛性[9]。国内有学者提出基于奇异值分解(Singular Value Decomposition,SVD)技术的显式四维变分同化法[39]及基于本征正交分解(Proper Orthogonal Decomposition,POD)函数技术的显式四维变分法[40]等,这两种方法在某些地方是相似的,都是通过对协方差矩阵进行分解重构来缩减计算量,且均无需伴随模式的求解及切线性假设,能够有效减少计算量,但其稳定性仍需进一步测试。
除此以外,4D-Var(特别是强约束4D-Var)在同化时间窗口内隐含了“完美”模式这一假设,当模式误差大的时候,这一假设本身就不成立。此外,实施伴随算子的编码本身也是一件相当繁重和复杂的工作,物理过程参数化也会引起目标泛函产生不连续问题[35],对同化时间窗口长度的确定没有形成较统一的方法[41]。
海洋科学的发展,特别是物理海洋学研究和业务化海洋预测预报技术的发展,长期以来一直受到观测资料不足的制约。然而,这一状况在近十年中得到了根本改善。随着海洋观测技术的发展,特别是TOPEX/Poseidon、Jason-1和NOAA等带有探测海面高度和海面温度等海洋要素的卫星成功发射,以及全球 Argo实时海洋观测网的建立等,人们开始有能力获取广阔洋面上大量的、高分辨率的实时海洋观测数据。与此同时资料同化技术在海洋科学研究中也取得了长足的进步,并在业务化海洋预报的构建、海洋再分析资料集的制作和深化对物理海洋现象的认识等方面取得了许多可喜的应用成果。
长期以来,海洋预报的准确度由于受到观测区域的覆盖有限,且资料观测周期短、均一性差等因素的影响,尽管计算机容量和计算速度已经是今非昔比,但预报结果依然不尽人意。究其原因,除了上面所述的缺陷外,还与模式本身的误差有关。于是人们发现,资料同化技术能将各种不同类别和时间段的观测资料不断地融于数值模式,并将短期分析预报结果作为模式预报的初值,以此将观测与模式的结果不断融合成为一个最优值,最终减小误差,提高预报精度。为此,许多国家首先针对海洋预报建立了专门的研究计划,利用资料同化技术并结合海洋模式开发出相应的海洋预报系统。
澳大利亚的海洋模型分析预报系统[42](Ocean Model,Analysis and Prediction System,即OceanMAPS),就是利用中尺度海洋模式4(Modular Ocean Model 4,即MOM4)并结合集EnOI开发而成。其中,观测误差协方差考虑了仪器测量误差与代表性误差,背景误差协方差利用集合样本生成,并假设非同质及各向异性[42],这一误差方案能很好地代表不同时空范围内的海洋变化状态,同化获得的分析增量能够反映出不同区域的冷暖涡状况,对围绕澳大利亚海域的环流特征及中尺度现象等能够获得准确的预报结果。该系统同化资料有 SST(AMSR-E),T/S 剖面和SSH(Jason 和 Envisat)等。
英国的预报海洋同化模型系统[43](Forecast Ocean Assimilation Model,即FOAM)采用了NEMO模式和较为简单的同化技术,即分析订正与增量分析更新方案(ACS和IAU),这类似于最优插值法。利用的同化资料主要有:SST(OSTIA)、T/S剖面、SSH(Jason、Envisat、GFO及SSALTO/DUACS),以及海冰浓度。该系统假设预报误差来源于2个方面[43],一个是模式外强迫误差(风,热通量,淡水通量强迫等),另一个是模式自身内部动力过程引起的误差,温盐剖面资料的同化对模式温盐场(尤其盐度场)的改进较大,进而对密度场能做出快速的调整。FOAM对温盐的预报最大误差出现在上层海洋(主要在混合层),均方根误差分别为1.3oC和0.3,而在 1 000 m以下,误差分别降到0.5oC和0.1[43]。
美国海军的复杂坐标海洋模式/海军耦合海洋资料同化系统[44](Hybrid Coordinate Ocean Model/Navy Coupled Ocean Data Assimilation,即HYCOM/NCODA)采用多变量OI同化技术,将背景误差协方差分离为背景误差方差与相关矩阵,考虑垂向和水平相关,并以流-依赖(flow-dependent)部分来修正这种关系,变量间的相关性采用自相关模型函数来刻画。利用的同化资料主要有:SST(卫星和现场),海冰浓度、SSH(Jason、Envisat、GFO及NAVO ADFC)和T/S剖面等。该系统目前已经业务化,对全球海洋SST的预报均方根误差可以控制在1oC以内,海面高度预报误差则在3 cm以内[44]。
日本气象研究所开发的多变量海洋变分估计系统[45](The Meteorological Research Institue(MRI)multivariate ocean Variational estimation,即MOVE/MRI.COM)采用了3D-Var,且同化中的背景协方差矩阵由气候态温盐垂直耦合 EOF来分解[45]获取。利用的同化资料主要有:SST(MGD),SSH(Jason和Ensivat)和T/S剖面。该系统对北太平洋中尺度涡的分布特征有较为清晰的描述,对黑潮输运量的模拟预报与由高度计资料估算的输运量相关系数最高可达0.59[45]。
意大利地中海预报系统(Mediterranean Forecast System,即 MFS)[46]采用 3D-Var和 NEMO模式开发,具有较高的分辨率(水平 1/16o,垂向71层)。利用的同化资料主要有 SST(AVHRR),SSH(Jason、Envisat、GFO 及 NAVO ADFC)和T/S剖面。MFS是较早实现业务化海洋预报系统之一,并经过不断改进[47],目前该系统基于新的背景场误差协方差模型而发展起来的,考虑了复杂的海岸边界条件及海底地形变化,同化试验表明海表高度误差随着预报时间的推进不断减小并趋于稳定,同化后,海表高度误差能降低到5 cm以下,而在未同化观测资料时,此误差是递增的。
中国国家海洋环境预报中心(The National Marine Environment Forcast Centre of China,即NMEFC)开发的海洋系统[48],其背景协方差矩阵构造与HYCOM/NCODA系统类似,均将水平和垂直方向的相关性分离[48],但 NMEFC采用 3D-Var同化技术,利用非线性的温盐关系作为平衡约束并以此构建协方差矩阵。目前,该中心不仅利用NMEFC实时发布热带太平洋温、盐度场数据产品,还对ENSO进行实时监测。利用的同化资料主要有:SST(Reynolds)、SSH(Jason)和T/S剖面剖面。
综上所述,各国的海洋预报系统尽管利用的数值模式不尽相同,但同化技术几乎均采用了较为先进、也是比较成熟的OI、EnOI和3D-Var等方法;且利用的同化资料除了由卫星观测的 SST和 SSH外,大都采用了温、盐度剖面资料,以及海冰浓度资料等。为此,各预报系统的预报精度都有了不同程度的提高。
资料同化技术除了被广泛应用于构建完美的海洋预报系统外,人们还利用它能有效的将各种类型的海洋观测资料融入海洋模式中,生成时空分布更加完善的分析资料的特点,广泛用来制作海洋再分析资料集,以便充分利用通过现有观测技术所能获得的全部信息,尽可能真实地揭示海洋的真实状态。
一种被称为简单海洋资料同化(Simple Ocean Data Assimilation)的再分析资料集,即SODA,应该是资料同化技术在海洋科学中应用的成功典范,曾被各国海洋和大气科学领域中的专家、学者广泛运用。早在1999年,Garton等[49,50]利用MOM2模式并结合 OI技术,利用当时能获得的所有观测资料,如MBT、XBT、CTD、SST、岸站资料、水文资料和SSH等,从而构建了海洋科学领域首个再分析资料集,即包含了1950年至1995年期间全球海洋上层海水温度、盐度和海面高度及流场等要素的月平均再分析资料集。
当时,该资料集的水平分辨率为0.5o×0.5o,垂向分40层。同化考虑了背景场误差的空间相关性,其预报误差是基于统计估计理论来计算的,并假设其不随时间变化而仅与空间尺度相关[49]。初步的检验结果表明,SODA资料集反映的EI Nino及太平洋-北美(Pacific-North American)异常信号与观测资料提供的信号比较相符。但早期SODA资料集的不足在于[50],同化缺乏海表盐度资料,限制了对混合层盐度年际变化的正确估计;其次,SODA资料集对温跃层及一些典型水团的变化形成过程刻画不充分,这与应用的OI同化技术的自身缺陷有关,因为 OI所用的协方差矩阵是固定的,且不随时间变化,这就使得它不能将动力模式和观测信息很好地融合在一起,这势必会引起误差;再者SODA资料集对地形变化不敏感,难以捕捉近海环流特征及中尺度涡现象;此外,SODA资料获得的海洋要素在南半球和印度洋存在较大误差,而这与观测资料的稀疏存在一定关系。SODA资料集在不断更新中,Carton等[51]2008年利用类似的方法制作了平均分辨率为0.25o×0.4o,垂向40层的再分析资料,分为5天平均和月平均,资料年限为1958-2001年,少量的 Argo温盐剖面资料加入到同化中。另外,SODA资料集的几个最新版本已在其网站(http://www.atmos.umd.edu/~ocean/)发布。
不同于SODA资料集,由日本科学家利用其开发的多变量海洋变分估计与综合海洋模式、预报、分析和合成系统[52](Multivariate Ocean Variational Estimation &Comprehensive Ocean Modeling, Prediction, Analysis and Synthesis System in the Kuroshio region)生成的再分析资料集,即 MOVE & COMPASS-K,这是一个专门应用于黑潮区域的再分析资料集,利用了 3D-Var和Nuding同化技术,并同化了现场温、盐度观测值(来自于 GTSPP、XBT、TAO、Argo等观测手段)和海表高度(SSH)、SST等来自于卫星观测的资料。MOVE & COMPASS-K 再分析资料集的水平分辨率在北太平洋(23oN–45oN 和 120oE–180oE)海域为 1/4o×1/4o;其他区域为 0.5o×0.5o,垂向为 21层,资料年限从 1993-2001年。可见,虽然应用了全球Argo实时海洋观测网提供的温、盐度剖面资料,但由于该观测网那时还刚刚开始建设,相信不会有太多的 Argo资料应用其中。而由 MOVE &COMPASS-K资料集[51]估算的黑潮最大(最小)流速与实测ADCP值均方根误差为0.749 m/s(0.271 m/s),并能再现黑潮在日本南部海域的流轴结构。不过,对琉球群岛流系的体积输运量估算则存在 3.1 Sv(1Sv=106m3/s)的偏差。
中国在制作海洋再分析资料集方面起步较晚。2007年,国家海洋环境预报中心利用业已建立的热带太平洋温度盐度同化业务化系统,开始在其网站(http://dell1500sc.nmefc.gov.cn/qxzl/qxzl.asp)上发布热带太平洋温度盐度再分析产品,使中国成为继美国、日本、法国、意大利和加拿大等发达国家之后具有发布热带太平洋温度盐度再分析产品的国家之一。该系统能对海洋卫星遥感、Argo及TAO等观测资料进行同化分析处理,实现了热带太平洋海域温、盐度多元资料同化再分析的业务化运行。该产品每月(25日左右)发布,主要提供热带太平洋区域(120oE-70oW,30oN-30oS,水平分辨率为经向2o、纬向1o;垂向10~630 m)的温度、盐度月平均同化再分析场。2009年,国家海洋信息中心发布了另一个海洋再分析资料集,简称 CORA(China Ocean Reanalysis)[53]。该资料集借助普林斯顿广义坐标系统海洋模式(Princeton Ocean Model with generalized coordinate system,即POMgcs)制作完成。考虑到海洋观测数据的时空分布不均匀性,为了有效提取海洋观测数据中的多尺度信息,并结合海洋要素的变化特征及其观测数据的时空分布特点,研发了多重网格三维变分海洋数据同化方法,用其同化温、盐度现场观测和卫星测温以及卫星测高数据。CORA再分析产品时段是从1986年1月至2008年12月共23年;海区范围为99oE~148oE、10oS~52oN,即同时包括渤海、黄海、东海和南海及其邻近海域;时间分辨率为历年、各月,空间水平网格分辨率为0.5o×0.5o、垂向为25层。
此外,美国的Stammer等[54]和澳大利亚的Oke等[31]还分别利用麻省理工一般环流模式(Massachusetts Institute of Technology General Circulation Model,即MITgcm)和澳大利亚海洋预报模式(Ocean Forecasting Australia Model,即OFAM),并运用4D-Var和EnOI同化技术,分别制作完成了 1992-1997年间的海洋环流及气候估计-1(Estimating the Circulation and Climate of the Ocean-1,即 ECCO-1)资料集和 1992-2004年间的BRAN(Bluelink ReAnalysis)再分析资料集,而他们利用的同化资料不仅有卫星遥感观测资料(如SSH、SLA等)、长期锚碇站观测资料(如TAO)和短期临时站观测资料(如XBT,CTD等),还有长期漂移浮标观测的资料(如ALACE,Argo等)。由此可见,再分析资料集是一种将不同时间段内各种类别零散的观测资料,通过资料同化技术并结合数值模式进行综合分析而生成的网格化数据集,它能再现海面高度、三维温度、盐度和流场的时空演变过程,为系统深入研究海洋对气候变化的响应、大洋环流与温盐结构、海洋过程与现象,以及监测方案的设计等提供更加丰富的信息[53]。
准确地预报海洋,不仅需要大量的现场观测资料、一个完美的数值预报模式,更需要对引起海洋变化的各种现象,特别是物理海洋现象的全面而又系统的认识。资料同化技术可以帮助人们从大量的观测资料中提取出更多、更有用的信息,从而深化人们对各种物理海洋现象的认识。
早期,Derber等[55]的工作仅限于同化温度资料,包括海表及垂直剖面温度资料,他先对温度资料进行严格的质量控制,随后借助 3D-Var同化工具,利用观测资料连续订正模式解,其每步订正均插入前后15 d(也即30 d的同化窗口)的观测值,有效弥补了观测资料不足的问题;最后将同化得出的SST与气候中心业务化SST产品进行了比较,发现在大尺度特征上两者表现出一致性,但相比过于平滑的业务化SST,同化后的SST在捕捉小尺度变化特征上表现出一定的优势。此外,同化温度垂直剖面资料有效地修正了温跃层和混合层的深度,并改进了对赤道潜流的强度及其垂向结构的描述。
Stammer等[54]利用MITgcm模式和4D-Var同化技术来研究全球海洋环流现象,同化利用了温盐剖面、SST及卫星高度计资料等。同化有效修正了湾流延伸体的空间结构,由于模式分辨率(水平2o×2o,垂向32层)相对较低,并且用于同化的次表层观测资料不足,同化对次表层海洋运动(比如典型水团的形成和运动)特征的刻画效果欠佳。
Masuda等[56]则利用4D-Var和MOM模式构造一个精度稍好(水平1o×1o,垂直36层)的同化系统来改进模式对北太平洋动力场的估计。首先用气候态风场和温盐场强迫并采用温盐张弛边界条件使模式积分 60年,以最后一年结果为初始场,随后利用4D-Var将WOA98气候态月平均和Argo温盐剖面、OISST及TOPEX/Poseidon卫星10天平均海面动力高度异常值同化进入模式,得益于分辨率的提高及同化资料的增加,同化的效果是明显的:同化强化了北太平洋中层水(North Pacific Intermediate Water,NPIW)的南侵,削减了南太平洋热带水高盐区的影响范围,这些与观测事实更为相符,并调整了净海气热通量在热带太平洋中的分布结构,使其与卫星观测结果趋于一致;另外,还对黑潮过于北侵及亲潮显著性北退的分布特征得到了修正。
Ishikawa等[57]的工作更令人鼓舞,它用4D-Var构建西北太平洋高精度(水平 1/6o×1/8o,垂直 78层)同化监测系统,该系统同化的资料包括现场观测的常规温、盐度资料和卫星反演的SST、SSH数据。对西北太平洋的两大重要流系(黑潮和亲潮)进行了深入研究,发现同化手段不仅优化了模式对二者输运量的模拟结果,还修正了两大流系交界处的锋面结构以及中尺度涡的空间分布特征。
朱江等[10]曾用 OI发展了一个海温短期数值预报系统。他首先利用最优插值质量控制法对船舶报SST资料进行质量控制,然后将资料插值到模式格点进行同化,类似于Derber等[55]的方案,在每次用OI 进行SST同化过程中,加入了前后相差6 h的相邻资料作为同时资料来解决资料稀释问题。同化使得绝大部分海域的SST得到修正,观测密集区的误差可以降低到 0.6oC,另外同化后黑潮(暖)、黄海暖流、朝鲜西岸沿岸流(冷)等主要流系的分布得到了更清楚的反映。近些年来,朱江等[33]还利用研发的 OVALS(Ocean Variational Analysis System)同化系统,并结合一个大洋环流模式,探讨了热带太平洋海域的三维温、盐场。在海面高度资料同化中引入了一个新的基于三维变分的同化方案,并考虑了背景场误差的垂直相关性和非线性的温-盐关系[58]。其中非线性的温盐关系是通过将Derber等[59]的线性平衡约束方案推广到非线性[60]的情况来考虑,从而可以通过同化高度计资料来直接调整模式的温、盐度场。由热带太平洋的同化试验结果表明,经过同化后的温度和盐度场较非同化试验产生了显著的改善。如利用独立的月平均温度和盐度的 TAO 观测资料进行检验, 温度场的均方根误差从 2.35℃减小到了 0.63℃,盐度场总的均方根误差从0.78 psu减小到了0.34 psu。尤其在南、北纬的副热带海域和东、西沿岸地带的盐度场修正得更大,更接近实际观测。高山等[34]利用该同化系统研究了西北太平洋中尺度涡的特征,他们将卫星高度计资料同化反演成为温盐“伪观测”数据, 然后再进行常规温盐同化,得到温盐分析场。由高度计反演温盐场的同化方法,极大地提高了中尺度涡模拟的精确性,且模拟的中尺度涡能够迅速响应温盐场的变化,并对自身流场做出快速调整。
王东晓等[61]则尝试了卫星高度计资料的同化,采用简单的牛顿松弛逼近同化技术,仅提取高度计资料的海表信息,未考虑垂向投影技术,同样取得了不错的效果。同化结果有效地修正了南海大尺度环流特征,南海天气尺度涡旋亦被同化所“唤醒”。随着同化技术的发展,韩桂军等[62]采用基于最优控制理论的伴随同化技术建立潮汐模型,在同化验潮站水位和T/P测高资料后,发现模型开边界条件中M2分潮调和常数得到了优化,说明同化对模型的校正是十分有效的。Han等[63]还采用相同的方法进一步拓展了上述工作,随着观测资料的逐步融合,水深、底摩擦系数及潮汐开边界条件等未知参数得到有效校正,发现经过参数优化后的模式对水位的模拟精度显著提高了,对中国近海 14个潮位站模拟的水位最大均方根误差由原先的75 cm降到45 cm。另外,近些年来,国内关于资料同化在南海的发展及应用方面也取得了一些有益的成果,有兴趣的读者可以参考文献[64,65]。
资料同化在海洋科学领域中的应用相当广泛,除了上面介绍的几种应用外,在为海洋模式提供初始场[53],确定海洋模式未知参数[62,63]和反演模式外强迫场(如风场,海气界面热通量)[56,66]等方面,也有许多学者进行了尝试,并取得了丰富的应用研究成果。因限于篇幅,这里不再赘述。
综上所述,进入 21世纪以来,无论是资料同化技术本身还是其在海洋科学领域中的应用都得到了快速的发展,但也还存在一些有待解决的“瓶颈”问题。就资料同化技术的发展而言,如今变分同化和集合同化(分别以4D-Var和EnKF为代表)已经发展成为两大主流方向。4D-Var在理论研究上相对较为完善,但其计算量庞大依然是人们重点关注的问题;同样,EnKF的计算量也相当大,并且仍有理论问题尚未完全解决,如集合样本生成和滤波发散等问题。不过,EnKF以其简便易行的算法和易于操作的优点受到越来越多人的关注,有学者[67]通过同化雷达数据来比较 2种方法(EnKF和4D-Var)的优劣。初步的结果表明,随着同化时间的推进,EnKF的效果要优于 4D-Var;且 EnKF对模式误差较 4D-Var来说更为敏感[41]。这些研究结果或许揭示了这样一个事实,即EnKF较4D-Var来说,具有更大的发展潜力。
另就资料同化技术在海洋科学中的应用发展而言,早期的海洋资料同化仅考虑了温度的同化调整而忽略了盐度的变化,从而导致密度场被严重恶化,同化后的结果甚至比没有同化任何观测资料时还要差[68,69]。为了弥补盐度观测资料的不足,解决同化中只有温度没有盐度的问题,需要把盐度和温度关联起来,才能有效减小同化中产生的虚假信息[69]。但已有的研究表明,将盐度资料同化进入模式当中,即便是考虑了适当的温盐约束关系,得到的盐度场也可能会与实际值存在偏差。因此,对盐度的同化尤应引起重视。这里不但存在温盐多变量约束问题,还可能存在海洋混合和热盐环流等这些难以刻画的现象[70]。国内有学者[33,53,71]对Argo温盐多变量同化做了研究,并通过同化高度计资料来直接调整温盐度场,同化表明 Argo资料能和其他资料形成很好的互补,在改进盐度场方面具有十分重要的作用。也有学者[70]针对 Argo温盐剖面的一致性尝试了将沿着等温线的盐度观测值 S(T)同化进入模式,而非常规盐度观测S(z)同化,结果表明采用这种盐度同化方式能有效保持海洋水团性质,因而温盐场,乃至密度场和流场在同化中可以得到一致性校正。因此,盐度资料的同化仍是未来的研究重点之一。
此外,卫星遥感数据和现场观测数据的多源同化问题。由于卫星遥感反演的是一种非常规的观测数据,且主要反映海洋表面的信息,其特点是时空覆盖率高;而现场观测资料(如Argo,XBT,TAO,CTD等)主要反映的是海洋表层以下的变化信息,其时空覆盖率相对较低,如何将它们有效融入数值模式,形成各种观测资料的互补,更好地揭示海洋中存在的多种尺度的物理海洋现象,并能摸清它们分布的特征和变化的规律,从而到达海洋预测预报精度的目标,这是资料同化技术在海洋科学领域中广泛应用所要解决的问题。虽然,目前国内有学者[33,53,71]对卫星高度计和其他常规资料结合同化方面做了一些有益的尝试,但离业务化应用尚有差距。此外,由遥感反演的其他海表信息如叶绿素、溶解氧、PH值等,如何在海洋科学研究中同化这类资料,也是需要进一步探讨的问题。
随着全球 Argo实时海洋观测网的进一步完善和长期(至少 10年)维持,海表盐度观测卫星的即将升空,以及计算机性能的不断提高,资料同化这一先进技术必将在海洋科学研究中得到愈来愈广泛的应用,并有望取得长足的进步,促进海洋科学的快速发展。
[1] 万利颖. 集合同化方法在太平洋海洋高度计资料同化中的应用研究 [D]. 北京: 中国科学院大气物理研究所, 2006.
[2] Lorenc A C. The potential of the ensemble Kalman filter for NWPA comparison with 4D-Var [J]. Quart. J. Roy. Meteor. Soc. , 2003b,129: 3183-3203.
[3] Panofsky H A. Objective weather-map analysis [J]. J. Appl. Meteor.,1949, 6: 386-392.
[4] Gilchrist B, Cressman G P. An experiment in Objective analysis [J].Tellus. , 1954, 6(4): 309-318.
[5] Bergthorson P, Döös B. Numerical weather map analysis [J].Tellus. , 1955, 7, 329-340.
[6] Cressman G P. An operational objective analysis system [J]. Mon.Wea. Rev. , 1959, 87: 367-374.
[7] Hoke J, Anthes R. The initialization of numerical models by a dynamic relexation technique [J]. Mon. Wea. Rev. , 1976, 104:1551-1556.
[8] Gandin L S. Objective Analysis of Meteorological Fields,Gidrometoro-logicheskoe Izdatelstvo, Leningrad(in Russian),English translation by Israel program for scientific translations,Jerusalem, 1963, 242.
[9] Bouttier F, Courtier P. Data assimilation concepts and methods March 1999, ECMWF培训教程.Available at: http: //twister. ou.edu/OBAN2006/Assim_concepts. pdf.
[10] 朱江, 徐启春, 王赐霞, 等. 海温数值预报资料同化试验Ⅰ. 客观分析的最优插值法试验 [J]. 海洋学报, 1996, 17(6): 9-20.
[11] Kalnay E. Atmospheric Modeling, Data Assimilation and Pridctability [M]. Cambridge: Cambridge University Press , 2003.
[12] Lorenc A. A Global Three-Dimensional Multivariate Statistical Interpolation Scheme [J]. Mon. Wea. Rev. , 1981, 109(4): 701-720.
[13] Kalman R E. A new approach to linear filter and prediction problems [J]. Transaction of the ASME, Journal of Basic Engineering, 1960, 82D: 34-45.
[14] Kalman R E, Bucy R. New results in linear filtering and predication theory [J]. Transaction of the ASME, Journal of Basic Engineering,1961, 83: 95-108.
[15] 高山红, 吴增茂, 谢红琴. Kalman滤波在气象数据同化中的发展与应用 [J]. 地球科学进展, 2000, 15(5): 571-575.
[16] Evensen G. Sequential data assimilation with a Nonlinearquasigeostrophic model using Monte Carlo methods to forecast error statistics [J]. J. Geophys. Res, 1994, 10: 143-162.
[17] Heemink A W, Verlaan M, Segers A J. Variance reduced ensemble Kalman Filtering [J]. Mon. Wea. Rev. , 2001, 129: 1718-1728.
[18] Pham D T, Verron J, Roubaud M C. Singular evolutive extended Kalman filter with EOF initialization for data assimilation in oceanography [J]. J. Mar. Syst. , 1998, 16: 323-340.
[19] Burgers G, Leeuwen P J, Evensen G. Analysis scheme in the ensemble Kalman filter [J]. Mon. Wea. Rev. , 1998, 126:1719-1724.
[20] Evensen G. The ensemble Kalman filter: theoretical formulation and practical implementation [J]. Ocen. Dyn. , 2003, 53: 343-367.
[21] Houtekamer P L, Mitchell H L. Ensemble Kalman filtering. [J]. Q. J.R. Meteorol. Soc. , 2005, 131: 3269-3289.
[22] Anderson J L. An ensemble adjustment Kalman filter for data assimilation [J]. Mon. Wea. Rev. , 2001, 129: 2 884-2 903.
[23] Bishop C H, Etherton B J, Majumdar S J. Adaptive sampling with the ensemble transform Kalman filter. Part I: Theoretical aspects [J].Mon. Wea. Rev. , 2001, 129: 420-436.
[24] Whitaker J S, Hamill T M. Ensemble Data Assimilation without perturbed observations [J]. Mon. Wea Rev., 2002, 130(7):1913-1924.
[25] Tippett M K, Anderson J L, Bishop C H, et al. Ensemble square root filters [J]. Mon. Wea. Rev. , 2003, 131: 1485-1490.
[26] Ott E, Hunt B R, Szunyogh I, et al. A local ensemble Kalman filter for atmospheric data assimilation [J]. Tellus., 2004(a), 56A:415-428.
[27] Hunt B R, Kostelich E J, Szunyogh I. Efficient data assimilation for spatiotemporal chaos: a localensemble transform Kalman filter [J].Physica. D. , 2007, 230: 112-126.
[28] Houtekamer P L, Mitchell H L, Pellerin G, et al. Atmospheric data assimilation with an ensemble Kalman filter: results with real observation [J]. Mon. Wea. Rev. , 2005, 133: 604-620.
[29] Houtekamer P L, Mitchell H L. Data assimilation using an ensemble Kalman filter technique [J]. Mon. Wea. Rev. , 1998, 126: 796-811
[30] Houtekamer P L, Mitchell H L. A sequential ensemble Kalman filter for atmospheric data assimilation [J]. Mon. Wea. Rev., 2001, 29:123-134.
[31] Oke P R, Schiller A, Griffin A, et al. Ensemble data assimilation for an eddy-resolving ocean model of Australian region [J]. Q. J. R.Meteorol. Soc. , 2005, 131: 3301-3311.
[32] Oke P R, Brassington G B, Griffin D A, et al. The Bluelink ocean data assimilation system (BODAS) [J]. Ocean. Modelling. , 2008,21: 46-70.
[33] 朱江, 周广庆, 闫长香, 等. 一个三维变分海洋资料同化系统的设计和初步应用 [J]. 中国科学D辑, 2007, 37(2): 261-271.
[34] 高山, 王凡, 李明悝, 等. 中尺度涡的高度计资料同化模拟 [J].中国科学D辑, 2007, 37(12): 1669-1678.
[35] Tsuyuki T, Miyoshi T. Recent progress of data assimilation methods in meteorology [J]. J. Meteor. Soc. Japan. , 2007, 85B: 331-361.
[36] 王斌, 赵颖. 一种新的资料同化方法 [J]. 气象学报, 2005, 63(5):694-701.
[37] Le Dimet F X, Talagrand O. Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects [J].Tellus. , 1986, 38A: 97-100.
[38] Courtier P. A strategy for operational implementation of 4D-Var,using an incremental Approach [J]. Q. J. R. Meteorol. Soc. , 1994,120: 1367-1387.
[39] 邱崇践, 张蕾, 邵爱梅. 一种显式四维变分资料同化方法 [J].中国科学D辑, 2007, 37(5), 698-704.
[40] 田向军, 谢正辉. 基于本征正交分解的显式四维变分同化方法:理论与验证 [J]. 中科科学D辑, 2009, 39(4): 529-536.
[41] Kalnay E, Li H, Miyoshi T, et al. 4-D-Var or ensemble Kalman filter? [J]. Tellus. , 2007, 59A: 758-773.
[42] Brassington G B, Pugh C S, Schulz E, et al. BLUElink>Development of operational oceanography and servicing in Australia[J]. Journal of Research and Practice in information Technology. ,2007, 39: 151-164
[43] Martin M J, Hines A, Bell M J. Data assimilation in the FOAM operational short-range ocean forecasting system: a description of the scheme and its impact [J]. Q. J. R. Meteorol. Soc. , 2007, 133:981–995
[44] Metzger E J, Hurlburt H E, Wallcraft A J, et al. Validation Test Report for the Global Ocean Prediction System V3.0–1/12°HYCOM/NCODA: Phase I .Available at: http: //www. dtic.mil/cgi-bin/GetTRDoc?Location=U2&doc=GetTRDoc. pdf&AD=ADA518693. orecasting system: first phase
[45] Norihisa U, Shiro I, Yosuke F, et al. Meteorological Research Institute multivariateocean variational estimation (MOVE) system:Some early results [J]. Advances in space Research, 2006, 37:806-822.
[46] Pinardi N. The Mediterranean ocean forecasting system: first phase of implementation (1998–2001) [J]. AnnalesGeophysicae. , 2003, 21:3-20.
[47] Dobricic S, Pinardi N. An oceanographic three-dimensional variational data assimilation scheme [J]. Ocean. Modelling. , 2008,22: 89-105
[48] Dombrowsky L, Bertino L, Brassington G, et al. GODAE systems in operation. Available at: http: //www.earth-prints. org/bitstream/2122/4779/1/20081214% 20Godae%20Paper. pdf.
[49] Carton J A, Chepurin G, Cao X. A Simple Ocean Data Assimilation Analysis of the Global Upper Ocean 1950–95. Part I: Methodology[J]. J. Phys. Oceanogr. , 2000, 30: 294-309.
[50] Carton J A, Chepurin G, Cao X. A Simple Ocean Data Assimilation Analysis of the Global Upper Ocean 1950–95. Part II: Results [J]. J.Phys. Oceanogr. , 2000, 30: 311-326.
[51] Carton J A, Giese B S. A Reanalysis of Ocean Climate Using Simple Ocean Data Assimilation(SODA) [J]. Mon. Wea. Rev. ,2008, 136: 2999-3017.
[52] Masafumi K, Tsurane K, Hiroshi I, et al. Operational Data Assimilation System for the Kuroshio South of Japan: Reanalysis and Validation [J]. Journal of Oceanography. , 2004, 60: 303-312.
[53] 韩桂军: 中国近海及邻近海域海洋再分析技术报告[R/OL]. 国家海洋信息中心, 天津. 2009[2010-10-21]. http://www.cora.net.cn/CORA技术报告.pdf.
[54] Stammer D C, Wunsch R, Giering C, et al. Global ocean circulation during 1992–1997, estimated from ocean bservations and a general circulation model [J]. J. Geophys. Res. , 2002, 107(C9), 3118,DOI10. 1029/2001JC000888.
[55] Derber J, Rosati A. A global oceanic data assimilation system [J]. J.Phys. Oceanogr. , 1989, 19: 1333-1347.
[56] Masuda S, Awaji T, Sugiura N, et al. Improved estimates of the dynamical state of the North Pacific ocean from a 4 dimensional variational data assimilation [J]. Geophys. Res. Lett. , 30(16), 1868,doi: 10. 1029/2003GL017604, 20032003.
[57] Ishikawa Y. High-resolution synthetic monitoring by a 4-dimension--al variational dataassimilation system in the northwestern North Pacific [J]. J. Mar. Syst., 2009, 78(2): 237-248, doi: 10. 1016/j.jmarsys. 2009. 02. 016.
[58] Chang X Y, Jiang Z, Rongfeng L, et al. Roles of vertical correlations of background error and T-Srelations in estimation of temperature andsalinity profiles from sea surface dynamicheight [J].J. Gepphys. Res. , 2004, 109, C08010, doi: 10. 1029/2003JC002224.
[59] Derber J, Bouttier F. A reformulation of the background error covariance in ECMWF global data assimilation system [J]. Tellus. ,1999, 51A: 195–221.
[60] Jiang Z, Changxiang Y. Nonlinear balance constrains in 3DVAR data assimilation [J]. Science in China(D). , 2006, 49(3): 331-336.
[61] 王东晓, 施平, 杨昆, 等. 南海TOPEX海面高度资料的混合同化试验 [J]. 海洋与湖沼, 2001, 32(1): 101-108.
[62] 韩桂军, 方国洪, 马继瑞, 等. 利用伴随法优化非线性潮汐模型的开边界条件Ⅱ: 黄海、东海潮汐资料的同化试验 [J]. 海洋学报, 2001, 23(2): 25-31.
[63] GuiJun H, Wei L, ZhongJie H, et al. Assimilated Tidal Results of Tide Gauge and TOPEX/POSEIDON Data over the China Seas Using a Variational Adjpint Approach with a Nonlinear Numerical Model. [J]. Advances In Atmospheric Sciences., 2006, 23(3):449-460.
[64] Ye Q S, Jiang Z, Dong X W, et al. Performace of four Sea surfacetemperature assimilation schemes in South China Sea[J].Continetal Shelf Research. , 2009, 29: 1489-1510.
[65] Xian J X, Dong X W, Jian J X. The assimilation experiment in the southwestern of Shouth China sea in summer 2000 [J]. Chinese Science Bulletin. , 2006, 51(z2): 31-37.
[66] Tziperman E, Thacker W C. An optimal control/adjoint equations approach to studying the oceanic general circulation [J]. J. Phys.Oceanogr. , 1989, 19: 1471-1485.
[67] Caya A, Sun J, Snyder C. A Comparison between the 4DVAR and the Ensemble Kalman Filter Techniques for Radar Data Assimilation [J]. Mon. Wea. Rev. , 2005, 133: 3081-3094.
[68] Troccoli A, Haines K. Use of the temperature–salinity relation in a data assimilation context [J]. J. Atmos. Oceanic. Technol. , 1999, 16:2011-2025.
[69] Troccoli A, Coauthors. Salinity adjustments in the presence of temperature data assimilation [J]. Mon. Wea. Rev. , 2002, 130:89-102.
[70] Gregory C S, Keith H. Evaluation of the S(T) assimilation method with the Argo dataset [J]. Q. J. R. Meteorol. Soc. , 2009, 135:739-756.
[71] 朱江, 闫长香, 万莉颖. Argo资料在太平洋海洋资料同化中的应用. [C]//许建平主编, Argo应用研究论文集, 北京: 海洋出版社,2006: 151-162.
Development of data assimilation and its application in ocean science
LI Hong1, XU Jian-ping1,2
(1. Second Institute of Oceanography, State Oceanic Administration, HangZhou 310012, China;
2. State Key Laboratory of Satellite Ocean Environment Dynamics, HangZhou 310012,China)
Development of data assimilation techniques, especially those based on two major theories of optimal control and statistical estimation, and their application status in the study of ocean sciences are reviewed in this paper.It can be seen that owing to the grate variety, huge amount and poor uniformity of the ocean observation data such as Argo(Array for Real-time Geostrophic Oceanography),TAO(Tropical Atmosphere/Ocean array), XBT(Expendable Bathythermograph), CTD (Conductance-Temperature-Depth), SST(Sea Surface Temperature)and sea surface height etc., and data assimilation techniques are facing various challenges. With the development of ocean science, the improvement of ocean numerical prediction accuracy, development of physical data products, and detailed description or deepened understanding of physical phenomenon can not be achieved without further development of data assimilation techniques.
data assimilation techniques; development; ocean science; application.
P731; O232
A
1001-6932(2011)04-0463-10
2010-10-21;
2011-03-10
国家重点基础研究发展计划课题(2007CB816001);国家海洋局第二海洋研究所基本科研业务费专项资助(1477-50);国家海洋局青年海洋科学基金资助项目(1059-50)。
李宏(1986-),男,硕士研究生,主要从事物理海洋学资料分析研究。电子邮箱:slvester_hong@163.com。
许建平,研究员,博导,电子邮箱:13805744970@139.com。