由于海气耦合系统的“记忆”主要存储在海洋特别是上层海洋中，尽可能准确地确定海洋的初始状态对于准确的短期气候预测至关重要(Jin et al., 2008)。本文介绍新完成的一个全球海洋资料同化系统(ZFL_GODAS)的设计方案和初步检验评估。

与澳大利亚国家气象局的海洋资料同化业务系统BODAS系统(Oke et al., 2008)类似，ZFL_ GODAS系统也是一个基于集合最优插值(Ensemble Optimal Interpolation，EnOI)方法的同化系统。EnOI方法是集合卡尔曼滤波(Ensemble Kalman Filter，EnKF)方法的一个简化版本，以历史模式状态样本(静态样本)取代EnKF中由多模式向前积分所得到的动态样本，从而大大减小了计算量，而且能较好地规避EnKF中潜在的滤波发散风险。与3DVAR等同化方法相比，EnOI有下列优势：第一，EnOI方法使用一个由海洋模式自由积分得到的静态样本来估计背景场误差协方差。这样的基于集合样本的背景场误差协方差具有多变量协变、各向异性的特征，且能反映海洋物理过程固有的空间尺度特征。第二，EnOI没有引入显式平衡约束假定(Burgers et al., 2002)，这种假定并不是处处成立的。第三，EnOI方法能够以一种简单直接的方式，同时同化各种不同类型的观测资料，而其他很多同化方法在同化海面高度异常(sea level anomaly，SLA)和温盐廓线资料时，往往需要一些专门的处理步骤(Cooper and Haines， 1996；Troccoli and Haines， 1999；Yan et al., 2004；Cummings，2005；Martin et al., 2007)。

ZFL_GODAS系统是一个短期气候数值预测业务系统的子系统，其主要功能是为大洋环流模式提供初始场，作为气候预测的基础。系统是在Fu et al.(2009a)工作的基础上，开发完成的。Fu et al.(2009a)开展了对卫星高度计资料的同化试验，尽管只同化高度计资料一种资料，系统对海面高度进行有效修正的同时，对海表温度(sea surface temperature，SST)、次表层温盐廓线及次表层的流速都有较明显的改善；他还做了EnOI与3DVAR同化方案的对比试验，结果表明EnOI同化结果要明显优于3DVAR方案(Fu et al., 2009b；Fu and Zhu， 2011)。

在其工作基础上，我们对该系统做了大量修改和完善：第一，对于系统的静态样本生成方式做了重要改进，静态样本是决定同化系统性能的最重要因素。经过试验比较，并参照前人经验，我们放弃Fu et al.(2009a)所采用的SVD方式，改为从模式积分结果中直接提取“快照”的方式取样；从模式月平均结果中取样，改为从模式输出的瞬时结果中取样；考虑样本的季节代表性，改为分月样本，每个月同化时都使用不同的静态样本，从而能够得到更准确的背景场误差协方差矩阵，同化效果更优。第二，为了建立能够业务运行的同化和预报系统，同时又能减少单个节点上内存的占用，我们突破一系列技术难题，开发了一套特色鲜明、负载均衡、高效的并行化同化程序。第三，经过大量的敏感性试验，优化了系统中的若干参数。第四，扩展了同化资料的种类，使其能够同化SST、SLA以及各类现场温盐廓线观测资料。第五，对资料的预处理、质量控制模块以及后处理模块进行了优化和完善，改进了同化模块与海洋模式进行数据交互的方式，优化了运行流程等。

ZFL_GODAS系统所采用的海洋模式为中国科学院大气物理研究所(IAP)大气科学和地球流体力学数值模拟国家重点实验室(LASG)开发的气候系统海洋模式LICOM1.0(刘海龙等，2004)。

该模式是在LASG第三代大洋环流模式基础上发展起来的，是一个非刚盖近似的完全原始方程模式，垂直坐标采用η坐标(Zhang et al., 1996)。模式设计的目的就是提供能够模拟大尺度风生环流和热盐环流的数值海洋模式，用以作为耦合气候模式的海洋分模式。LICOM1.0模式的水平分辨率为0.5°(纬度)×0.5°(经度)，为准全球版本，纬度范围为90°N～78.5°S。垂直方向分为30层，其中海面以下0～1000 m深度有20层。模式地形采用DBDB5海洋深度资料，分辨率为(1/12)°(纬度)×(1/12)°(经度)，使用时平均到(1/3)°(纬度)×(1/3)°(经度)的网格上。在处理次网格过程时，采用了湍流黏性/扩散假设。水平方向湍流黏性项采用常系数的Laplace形式，湍流扩散使用了常系数的Laplace形式、热带海洋的P−P垂直扩散方案(Pacanowski and Philander， 1981)以及Gent and McWilliams(1990)用来描写中尺度涡旋对于大尺度位温、盐度等示踪物作用的参数化方案(简称GM90)。

模式首先从静止状态开始spin-up运行，用于初始化的温盐初始场来自Levitus94气候温盐资料集。spin-up运行使用的强迫场资料，包括风应力、净短波辐射、非短波热通量、耦合系数采用由德国马克斯－普朗克气象研究所(MPI)整理的海洋模式比较计划(OMIP)的月平均强迫场。spin-up运行960年后，模式达到稳定状态，之后模式改由欧洲中期数值天气预报中心的ERA40逐日风应力驱动。用于恢复的月平均海表温盐资料来自美国国家海洋资料中心(NODC)发布的《世界海洋图集1998》(WOA98，Levitus et al., 1998)。

3 EnOI同化方案

对海面高度h₀、温度T、盐度S，以及水平流速(u，v)的分析，通过求解如下分析方程得到(Evensen，2003)：

对于EnOI方法而言，静态样本的选取非常重要，它在相当程度上决定了最终分析的效果。好的样本应该是有足够代表性的，并且独立性较好的，其平均值与分析时刻的“真实”状态场越接近越好，而不同样本之间的离散度则要与状态变量的自然变率大小相当。事实上，“真实”状态场是未知的，如果选取的样本平均值与气候平均值很接近，而样本离散度与状态场中各变量的气候变率大体相当，则这个样本就是比较好的样本。海洋的状态变化存在显著的季节循环，如果忽略季节循环，同化中都使用相同的样本，则样本的代表性就不够好，样本均值与分析时刻的海洋状态可能相差较大，样本离散度也往往比该季节(或月份)的实际变率大得多，从而导致分析结果恶化。理论上，样本数越多，由有限个样本统计得到的背景场误差协方差将越趋近于实际的背景场误差协方差，分析结果越理想。但实际上，样本数的增加，会造成计算机内存需求和计算开销的大幅增加，因而实际应用尤其是业务应用中，样本数都不能过大。研究表明，100个左右的样本数是适当的(Mitchell et al., 2002)。

基于上述考虑，本文分析方案中，静态样本取为96个，且采用分月样本，每个月进行分析时使用不同的样本。静态样本由控制试验1970～2001年32年的模式积分结果中提取1152个瞬时模式状态场，将其按不同月份分别保存到12个静态样本数据文件中，每个文件包含96个样本。

背景场误差协方差在每个观测点周围，使用局地相关函数ρ进行了水平方向的局地化处理。ρ的定义采用Gaspari and Cohn(1999)提出的准高斯函数形式。研究表明，EnKF或EnOI中，通过局地化可以减小取样误差的影响(Hamill et al., 2001; Oke and Schiller， 2007)。本文试验中，影响半径取为180 km。

尽管EnOI与EnKF相比，计算代价大大减小，但由于系统计算范围为全球且分辨率较高，计算量巨大，串行程序远不能满足业务需要，对其进行并行化是必要的。针对EnOI同化程序的特点，我们开发了一套特色鲜明、负载均衡、高效的并行化同化程序。并行化采用MPI方案来实现，使用了多种节省内存的技术；使用并行IO，节省了大量通信时间；进行了巧妙的独特的计算量分割技术，采用相邻纬线交错分配的并行方案，最大限度保证了负载的均衡。

本系统能够同化卫星高度计资料、卫星SST资料，以及Argo、XBT、TAO等各种不同来源的现场温盐廓线资料。本文的评估试验中同时同化了上述3种资料，同化间隔都是每隔7 d同化一次。

其中，SST资料使用的是NOAA/NCDC发布的海表温度的逐日高分辨率融合分析资料(Reynolds et al., 2007)，资料覆盖范围为全球，水平分辨率0.25°(纬度)×0.25°(经度)，时间分辨率1 d。同化的海面高度资料使用的是法国SSALTO/ DUACS延时的多任务融合SLA格点产品。该资料融合了Jason-1、Jason-2、T/P、Envisat、GFO、RS1/2和Geosat等卫星高度计产品。资料采用(1/3)°的Mercator 投影网格系统，总共1080×915个网 格，资料覆盖范围从82°S～82°N。时间分辨率为7 d，资料起始时间为1993年1月。本文试验中，同化时对SST和SLA资料进行了稀疏化处理，采用简单的均匀跳点取样方式，资料稀疏到1°(纬度)×1°(经度)。

同化的现场温盐廓线资料使用的是英国气象局哈得莱中心(Hadley Center)的经过质量控制的EN3次表层海洋温盐观测资料集。用于生成该资料集的资料来源包括WOD05、GTSPP、Argo以及ASBO项目，收集处理的资料类型包括机械式温深仪(MBT)、投弃式温深仪(XBT)、电导—温度—深度测量仪(CTD)及其它水道测量廓线、锚定浮标阵列(TRITON/TAO、PIRATA等)、漂流浮标(主要是Argo)等。对这些原始资料进行了复杂、严格的质控，最后生成可直接用于资料同化的高质量的数据集。

除了资料本身提供的质量控制外，实际同化时还进行了如下的质量控制(朱江等，2007)：

(1)观测资料不能出现在由模式定义的陆地上。

(2)值域控制：海温观测数值介于－1～32 °C，盐度数值介于30～40 psu，海面高度异常介于－2～2 m。

(3)所用资料无重复记录。

(4)如果某观测与背景场在该点的插值之差，超过事先设定的阈值(温度、盐度、海面高度分别设为5 °C、4 psu、1 m)，则剔除该观测。

本节对ZFL_GODAS系统的同化性能进行评估。首先，LICOM模式由静止状态启动，在气候月平均大气强迫场驱动下，spin-up积分960年，模式达到稳定状态。以spin-up运行960年的模式结果为初值，用ERA40逐日风应力资料驱动，不加海洋资料同化，从1961年1月积分到2002年8月，这个试验称为控制试验。

以控制试验1996年7月1日的预报结果作为初值，模式向前积分过程中，系统分别同化各种温盐廓线资料、卫星SST资料、卫星高度计资料，由1996年7月至2001年12月。同化间隔都是每隔7 d同化一次。这5.5年的试验称为同化试验，取 1997年1月至2001年12月的结果进行分析。

首先，来考查海表温度和海面高度的年际变率。图 1和图 2分别给出了由1997～2001年的60个月的月平均资料计算的海表温度和海面高度的变率(标准差)的全球分布，计算标准差时使用的资料都是月平均资料，因此月以下时间尺度的高频变率一定程度上被滤除了。

图 1(a)控制试验、(b)同化试验、(c)卫星观测资料计算得到的1997～2001年的海表温度变率(标准差)的全球分布 Fig. 1 Sea surface temperature(SST)variability(st and ard deviation)during 1997−2001 from(a)control run，(b)assimilation run， and (c)satellite observation

图 2同图 1，但为海面高度变率 Fig.2 Same as Fig.1, but for sea level anomaly (SLA) variability

由图 1可见，LICOM模式模拟的SST 5年变率与观测比较一致，抓住了SST变率全球水平分布的主要特征：第一，低纬和高纬变率较小，中纬度变率较大；第二，全球SST变率最大的区域在30°N～60°N之间的北太平洋西部和北大西洋西部，其中西北太平洋变率最强；第三，从秘鲁沿岸向西延伸到赤道中东太平洋的“冷舌”区是热带海表温度变率最大的地方，变率超过2.5 °C。但总体上，模式模拟的SST变率整体偏小，特别是在西北太平洋、湾流区、东赤道太平洋等上述SST变率较大的区域，显著偏小。这可能与模式分辨率较低，不能分辨这些区域的中尺度涡有关。相比较而言，经过海洋资料同化，模式模拟的SST变率与观测更为接近，特别是赤道中东太平洋“冷舌”区、南大洋等区域，同化后更有明显改善。

海面高度与海水的温度、盐度及海流都有关系，其变率分布比SST复杂。观测SLA变率全球水平分布的主要特征：第一，在南北半球的西边界流附近以及南极绕极流附近，变率最大，最大超过0.4 m；第二，赤道中东太平洋、西太平洋的赤道南北两侧、印度洋赤道以南，海面高度变率也较大(＞0.1 m)；第三，热带大西洋、东南太平洋等处，SLA变率很小，小于0.05 m。SLA变率较大区域与全球主要流系对应较好。由图 2可见，LICOM模式模拟的海面高度5年变率整体偏低，形势与观测较为接近。特别是西边界流及其延伸体附近，显著偏低，这主要是由于在这些地方存在大量不断产生、移动、消亡的中尺度涡，导致变率很大，而LICOM模式分辨率较低，无法分辨这些中尺度涡所造成的。经过海洋资料同化，模式模拟的海面高度(sea surface height，SSH)变率也有了明显改善，与观测更为接近，特别是热带的几块大值区改善特别显著。在一些模式无法很好地分辨的中尺度涡活跃的地区，比如西边界流区和南极绕极流区，同化后变率也都有所改善，但仍比观测的变率低0.1～0.2 m。

　　从上面的对比分析可以看到，通过SST、SLA以及温盐廓线资料的同化，得到的海表温度以及海面高度的标准差在全球大部分海域更接近观测的变率。关于SST或SSH标准差的比较，Brankart et al.(2003)、Castruccio et al.(2008)、Fu et al.(2009a)等的研究中也有类似的结果。

赤道月平均海表温度异常(sea surface temperature anomaly，SSTA)随时间变化的经度－时间剖面图在图 3中给出。1997～2001年，赤道SSTA最突出的特征是1997～1998年的强El Niño过程，日界线以东直到南美沿岸长达10～16个月的持续强烈增温；其次是1998年下半年到1999年初，以及1999～2000年秋冬季节的两次La Niña过程。控制试验结果也抓住了这3个主要过程，但是不论是强度、范围、起止时间，与观测都存在一些差别，如对1997～1998年El Niño的模拟明显偏弱，最强增温仅3 °C，持续时间也较短，中太平洋还分裂为两次较明显的异常增温；而经过资料同化以后，模拟结果大大改善了，1997～1998年El Niño的强度、范围、起止时间，都与观测非常一致，最高增温大于5 °C，只是最强烈增温比实测范围更大一些。两次La Niña过程与之后的发展，以及印度洋和大西洋的发展演变特征，同化后也都有明显改善。与Fu et al.(2009a)的结果相比，也有明显改善，以1997～1998年El Niño期间的最强增温为例，无论是最强增温范围、出现时间、增温强度都更加接近观测，表明我们对系统所做的优化效果显著。

图 3（a）控制试验、（b）同化试验以及（c）由观测得到的1997～2001年赤道月平均海表温度异常随时间的变化 Fig.3 The evolution of sea surface temperature anomaly (SSTA) along the equator during 1997−2001 from (a) control run, (b) assimilation run, and (c) observation

图 4分别给出了同化和不同化两种情形下模式预报SST相对于观测(再分析资料)的均方根误差，以及Fu et al.(2009a)试验中的对应结果。我们看到，在全球大部分海区，LICOM模拟的SST的均方根误差(root-mean-square error，RMSE)都小于1 °C，西北太平洋、西北大西洋、赤道东太平洋、南大洋等海域误差较大，特别是东中国海沿岸地区、日本东北部近岸海域、湾流区等处，均方根误差达到5～7 °C以上。经过海洋资料同化以后，效果得到很大改善，上述几大海域绝大部分都降至1 °C以下，但包括中国东部沿岸、美国东北部沿岸等在内的几小块，RMSE仍然较大，但大误差区范围大幅缩减，RMSE的值平均减小2～3 °C。与Fu et al.(2009a)的结果相比，本文同化后的改善要大得多，这主要是由于本文试验将SST资料直接进行了同化，而Fu et al.(2009a)是通过同化SLA资料来间接的修正SST；另外，如引言中所述，我们的静态样本也从几个方面大大改善。

图 4 （a）控制试验和（b）同化试验预报的SST的RMSE Fig.4 The root-mean-square error (RMSE) of SST between (a) control run and (b) assimilation run with observation

我们还比较了同化前后的SST模拟结果与观测之间的相关系数，经过同化以后也有非常显著的改善，对于SLA也作了类似比较，限于篇幅，这里不再给出具体的分析。

为了更加客观地检验同化系统的性能，事先在TRITON/TAO锚定浮标温盐廓线资料中抽取了几个点的资料，没有参加同化，专门用于检验同化结果。下面考查控制试验和同化试验的结果，与这些完全独立的观测资料的对比情况。

图 5给出了4个赤道TAO浮标站点上控制试验和同化试验中的模式温度的均方根误差。可以看出，同化系统对温度的同化效果非常好，4个点在有观测的450 m以上，经过同化后的RMSE均比不同化的结果大幅减小，在200～350 m深度改善最大，最多减小1.2 °C。我们还比较了这4个站点5年的平均温度偏差(模式平均温度对TAO浮标观测温度平均值的偏差)廓线(图略)，同样看出，同化系统对温度的同化效果非常好，4个点在有观测的450 m以上，平均温度偏差更接近0。最大的改善出现在200～300 m深度，平均温度偏差减小达1 °C以上。

图 5 模式预报TAO浮标站点(a)(0°，147°E)、(b)(0°，180°)、(c)(0°，170°W)、(d)(0°，155°W)温度的均方根误差 Fig. 5 The RMSE of temperature with respect to four TAO station points at(a)(0°，147°E)，(b)(0°，180°)，(c)(0°，170°W)， and (d)(0°，155°W)

由图 6可见，同化系统对盐度的同化也不错，50～450 m之间，都有明显改善，将模式结果向观测拉近了。最大的改善在100～150 m深度，平均盐度偏差减小0.1 psu左右。但50 m以上的表层不太理想，同化后的表层盐度偏高。这部分是由于表层的盐度向气候平均盐度恢复的结果。

受资料所限，以上同TRITON/TAO锚定浮标温盐廓线资料的对比仅限于热带太平洋的少数几个站点。这不足以反映全球大洋同化之后的总体效果。但就这几个点而言，ZFL_GODAS同化系统效果是比较好的，特别是温度的同化效果更好，这可能与温度观测资料数量比盐度观测多得多有关。

图 6 1997～2001年TAO浮标站点(a)(2°N，156°E)、(b)(0°，156°E)、(c)(2°S，156°E)的5年平均盐度廓线 Fig. 6 The mean salinity profiles over the period 1997–2001 from control run，the assimilation， and the TAO mooring data at three TAO station points:(a)(2°N，156°E)，(b)(0°，156°E)， and (c)(2°S，156°E)

由于海流的观测非常稀少，因此ZFL_GODAS同化系统没有对海流资料进行同化，下面考查通过对SST、SLA和温盐廓线资料的同化，对次表层海流的影响。

从图 7可以看出，总体上，同化系统对纬向流速的同化效果也较好，模式对纬向流速随深度变化的趋势、最大流速所在深度模拟得比较准确，经过同化后，结果基本上都进一步向观测靠拢。在赤道潜流(Equator Under Current，EUC)的流轴附近，模式预报结果存在系统性的低估，经过同化以后都有所修正，流速增加2～7 cm/s。需要指出的是，与很多学者(Bellucci et al., 2007)的工作相比，我们的结果是比较好的。很多前人的工作进行资料同化以后，温度改善明显，而纬向流速廓线却变坏了，最大流速位置偏差较大。这主要是由于各变量单独同化，破坏了变量之间固有的物理约束引起不协调导致的。而我们的同化方案是多变量一起同化，一起构建历史样本，计算背景场误差协方差，因此同化一种变量时，对其他变量能够协调地一起修改。

图 7 1997～2001年赤道上的4个TAO浮标站点(a)(0°，165°E)、(b)(0°，170°W)、(c)(0°，140°W)、(d)(0°，110°W)的平均纬向流速廓线 Fig. 7 Mean zonal velocity profiles over the period 1997－2001 from control run，the assimilation， and the TAO mooring data at four points:(a)(0°，165°E)，(b)(0°，170°W)，(c)(0°，140°W)， and (d)(0°,110°W)

本文介绍了新完成的一个全球海洋资料同化系统ZFL_GODAS，该系统基于EnOI方法设计实现，使用中国科学院大气物理研究所开发的LICOM气候系统海洋模式。不同于许多其他的资料同化方案，系统能够以一种简单直接的方式，无需专门的处理步骤，同时同化SST、SLA以及各种来源的现场温度、盐度廓线资料等等，且对一种资料同化时，所有模式状态变量都能得到合理协调的调整。该系统目前的主要应用是作为一个短期气候预测业务系统的一个子系统，为该气候预测系统的海洋模式分量提供初始场。

本文介绍了应用ZFL_GODAS同化系统所做的5年同化试验，分析其同化效果。同化试验中对SST、SLA以及各种不同来源的温盐廓线现场观测资料都进行了同化。其中，SST和SLA资料是以卫星观测资料为基础，由专业数据中心处理发布的格点产品；而温盐廓线资料使用的是英国Hadley中心发布的EN3次表层海洋温盐观测资料集，资料时空分布不均，数量也少得多。

我们给出了一系列控制试验和同化试验结果以及观测资料的定量分析比较，用于对比的观测资料有部分(较小一部分)参加同化从而并非完全“独立”的卫星资料，也有事先预留出来没有参加同化的部分TRITON/TAO锚定浮标资料。通过比较海表温度和海面高度的年际变率，SSTA随时间变化，SST、SLA以及次表层温盐预报产品的RMSE，5年平均温度偏差廓线、平均盐度廓线、平均纬向流速廓线等，我们发现建立的系统工作正常，同化效果较好，经过同化以后，各变量都更加接近观测，误差更小，与观测场的相关性(图略)更好。

经过资料同化以后的模式预报月平均产品，全球95%以上的海域，SST的RMSE小于1 °C；SLA的结果略差，这可能与模式结果和卫星产品使用的作为基准的平均海面高度不同有很大关系，但即便如此，全球绝大部分海域的SLA的RMSE也小于0.1 m。而对于次表层温盐廓线，与独立观测资料的抽样比较表明，温度的最大RMSE出现在150 m深度附近，为1.3～1.8 °C；盐度的最大误差则出现在近表层，最大RMSE达0.45 psu。

最后，我们可以初步得出结论，ZFL_GODAS系统同化性能优良，可以为短期气候预测系统提供较为理想的海洋初始场，也可以为物理海洋学的研究提供有效的再分析资料。