Argo数据同化方法及其网格化产品研究进展

王丹阳，苏 涵，张春玲,2

(1.上海海洋大学 海洋科学学院，上海 201306;2.上海海洋大学极地研究中心，上海 201306)

引 言

Argo(Array for Real-time Geostrophic Oceanography)俗称“全球Argo海洋观测网”,类似于陆地上的探空气球，是由3000个自动剖面浮标组成一个庞大的全球海洋观测网，旨在准实时、大范围、高分辨率获取全球海洋次表层观测资料[1]，进而利用这些数据帮助我们更好地认知海洋内部的长期变化，提高气候预报精度，从而有效抵御全球日益严重的气候灾害给人类造成的威胁。2013年，国际Argo计划提出了从“核心Argo”向“全球Argo”拓展的设想，计划建成由4000个浮标组成的覆盖水域更深厚、涉及领域更宽广、观测时域更长远的全球Argo实时海洋观测网，并计划在2025年之前建成一个由4700个自动剖面浮标组成的全球(包含有冰覆盖的南北极海域和重要边缘海区域)、全海深(0～6000 m)、多学科(包括物理海洋和生物地球化学等10多个海洋环境要素)的实时海洋观测网。生物地球化学Argo(BGC-Argo)和深海Argo子计划应运而生，目前BGC-Argo已经成为全球Argo实时海洋观测系统中的重要组成部分。

自2000年以来，国际Argo计划各成员国已经在全球海洋陆续投放了约1.7万个自动剖面浮标，截至2020年9月底，在全球海洋上维持正常工作的浮标有3904个。这些浮标所获得的观测剖面数量已逾200万条，其中包含生物地球化学(BGC)环境要素的逾20万条。由于Argo的标准任务是一个“悬停漂移-剖面测量”的过程[2],即，浮标首先在1000 dbar深度层上悬停，并随洋流漂移，每10天下沉到2000 dbar，然后在上升到海面过程中收集海水的温、盐度数据，浮出水面后通过卫星将观测数据传回岸上，然后再次下潜到预设的漂移层上，继续下一个观测任务。因此，其观测剖面数据在空间上是散点观测，具有空间不均匀性，为了克服Argo剖面资料被广泛应用于气候预测和海洋、大气科学研究所受到的限制，国内外许多学者借助于一些常用的数据同化方法，如逐步订正[3]、最优插值(Optimal Interpolation, OI)[4]、集合卡尔曼滤波(Ensemble Kalman Filter，EnKF)[5]、集合最优插值(Ensemble Optimal Interpolation, EnOI)[6],以及三维变分(Three Dimensional Variational，3-D Var)[7]等等，研制了包含多个海洋要素、多种分辨率的网格化产品。

本文对目前适用于Argo剖面资料的常用数据同化方法、国内外基于Argo观测资料构建的网格化数据产品进行了比较系统地回顾和总结，对比分析了各种同化方法的适用条件和计算效率，以及各种Argo网格化产品的特点，并针对目前Argo资料的观测现状，给出了未来Argo数据同化的研究展望。

1 方法分析

1.1 早期方法

1949年Panofsky开创性地提出了多项式拟合(Polynomial Fitting，PF)[8]，多项式拟合是用一个多项式展开去拟合包含数个分析格点的一小块分析区域中的所有观测点，得到观测数据的客观分析场。展开系数用最小二乘拟合确定。但此方法的区域多项式拟合并不稳定，当资料缺测时更是如此,而且会导致分析在拟合的各个区域之间不连续。1954年Gilchrist和Cressman在多项式拟合的基础上，利用频谱分析，对气象数据所对应的频率进行选择，通过选择不同的频率成分来过滤错误信息，从而给出了函数滤波(Function Filtering，FF)[9]的简单函数模型。

而Gilchrist在1954年提出了理想逐步订正法(Successive Correction，SC)的原型[9]，后由Bergthorson[10](1955)对其进行理论论证，并由Cressman[3]发展成熟，形成基于迭代算法的逐步订正法。这种方法首先要求给出网格点的初始值，然后从每一个观测中减去对该观测点的估计值得到观测增量，通过将分析格点周围影响区域内的观测增量进行加权组合得到分析增量，再将分析增量加到背景场上得到最终的分析场，并进行逐步迭代，直到分析值达到某种预期的精度[11]。起初，Cressman给出了与距离平方成反比的二次权重函数，后来，Barnes[12]提出采用高斯型的权重函数,但该方案收敛较慢，一般需要迭代3～4次，之后，Barnes[13]对这一方案又作了改进，采用收敛因子来加快收敛速度。

以上这些早期的资料同化方法没有充分利用模式和观测资料的误差统计信息，也没有利用模式的时空演变信息，并且缺乏强有力的理论基础，都属于经验分析方法，在实际数值预报，特别是在海洋科学研究中并没有得到广泛应用。

1.2 现代方法

60年代初，最优插值法(Optimal Interpolation，OI)的提出，使得资料同化方法有了基于统计估计理论的基础。目前的数据同化方法根据其理论可分为两类，一类是基于统计估计理论的，如最优插值、卡尔曼滤波、扩展卡尔曼滤波、集合卡尔曼滤波、集合最优插值等；另外一类是基于最优控制理论的，如三维变分，强约束四维变分和弱约束四维变分等。这里主要总结回顾广泛用于Argo资料同化的几种方法的基本原理，主要包括：最优插值、集合卡尔曼滤波、集合最优插值、三维变分和四维变分等。

最优插值(OI)的分析场是背景场与由权重矩阵加权的修正量之和[4]采用最小二乘方法求得最佳线性无偏估计方程中的最优权重矩阵。此方法的基本假设是,对于一个模型变量,在确定它的增量时,只有几个观测值是重要的。因此，OI易于编码，计算量相对较小，但其分析结果并非全局最优，并且，OI是针对线性系统发展起来的，难以确保大小尺度分析的一致性。集合卡尔曼滤波方法[5]是在卡尔曼滤波[14]的基础上发展起来的，其基于蒙特卡罗算法，结合了卡尔曼滤波和集合预报的优点,用有限的集合样本来估算误差协方差矩阵的不确定性。EnKF概念简单,不需要作线性假设,也无需求解模式的切线性及其伴随,适合于并行计算等优点,是一个目前比较流行的方法。但EnKF计算量比OI要大得多，并且会因为样本集合离散度不够而产生样本误差问题，系统的非线性，及通常利用扰动观测法获取样本初值，使得这种样本误差问题更为明显。此后，Evensen[6]将集合思想吸收到最优插值同化技术中,提出了集合最优插值数据同化法：格点的分析值在一个固定的模式向量样本集合(如长时间序列的模式积分)空间内进行计算，模式统计误差不随时间变化,从而减少计算量。EnOI能够保持准动力一致性，避免假设均匀和各向同性等，但其在误差计算时仍沿用EnKF集合预报的方式，以获得较传统OI方法更优的分析值，由于模式误差不随模式积分时间改变，EnOI较EnKF得到一个次优解。

最初出现的变分同化方法是三维变分(3D-Var)，3D-Var以极大似然估计理论为基础，通过求解一个目标函数(也称代价函数)的极小值(一般利用目标函数的梯度求其极小值),产生一个分析时刻的综合考虑背景场和观测值的大气或海洋真实状态的最大似然估计，并且给出背景场和观测场各自相应的精度。3D-VAR的分析解为全局最优，可以处理观测算子是非线性的情况，并能同化各种不同来源的观测资料。但此方法是在某一时刻进行的分析,前一时刻的同化结果可作为后一时刻模式运行的初始场。但在使用时,无法用后面时刻的资料来订正前面的结果,同化的解在时间上不连续。为弥补3D-Var的这一缺陷，LeDimet等人于二十世纪八十年代提出了4D-Var[15]。4D-Var是在3D-Var的基础上，增加了时间变量的同化，在时间窗口内利用完整的动力模式作为强约束，自动调整模式误差，以便得到更精确的同化结果。并将某一时间段上的观测数据均纳入到同化系统，背景场误差协方差隐式发展，误差信息随动力模式而向前传播，这是4D-Var的主要优势。由于4D-Var需要求解伴随模式，并且代价函数求解通常采用最速下降法、共轭梯度法及准—牛顿迭代法等迭代计算，计算量特别大。

1.3 主要方法对比分析

进入本世纪以来，海洋资料同化技术取得了快速的发展，从早期比较简单的客观分析法，发展到现在能够同化大量非常规资料的四维变分和集合卡尔曼滤波等比较流行的方法。如上所述，每种同化方法均有其优缺点(表1)：逐步订正的优点在于计算量小，算法简单，计算速度快；缺点在于没有理论基础，属于经验分析方法；OI的优点在于计算量小、易于编码且可以得到较合理的分析，缺点在于它是局部最优、不随时间变化且主要针对线性系统；EnKF的优点在于它不需要求解模式伴随且适用于并行计算，缺点在于计算量比OI要大得多且样本集合离散度不够；EnOI的优点在于计算量小、能保持准动力一致性，缺点在于模式误差不随积分时间变化且仍采用集合预报的方式；3D-Var的优点在于它可以处理非线性算子并且可以进行全局分析，缺点在于计算量大、解在时间上不连续；4D-Var的优点在于它可以同化多时刻观测资料、可在目标函数上加上其他约束项，缺点在于计算量比3D-VAR大得多且难以得到伴随方程的离散形式。这使得每种方法在实际应用中都有其一定的限制和使用范围，而且不同方法对于资料质量的依赖性也不同，逐步订正、最优插值及三维变分都是基于观测误差恒定的假设，因此对于观测精度较低的数据，其同化结果精度也会随之降低。而卡尔曼滤波方法则可以通过滤波的方法，将异常观测数据加以平滑，同化结果对观测资料质量的依赖度较低。

表1 各种同化方法比较分析

2 Argo网格化产品研制现状

目前Argo剖面浮标获取的观测资料已经成为海洋气候模式中的重要数据来源，但是由于Argo剖面浮标具有“随波逐流”的特性，其观测剖面在位置和时间上有很大的随意性，这使得Argo观测资料的应用范围受到了一定的限制。随着数据同化技术的不断发展以及Argo观测剖面数量的持续快速增长，国内外学者利用各种同化手段，有效地将各种类型的海洋观测资料进行融合，并推出了一系列网格化数据集。

2.1 国内网格化产品

近年来，国内学者陆续研发了Argo网格化数据产品和再分析数据集，通过中国Argo实时资料中心网站(http://www.argo.org.cn/)公开发布，供广大用户使用。其中，中国Argo实时资料中心自主研发的网格化产品涵盖了多个要素，针对不同科学研究需求类型，并不定期的更新：《全球海洋Argo网格数据集(BOA_Argo)》是利用Barnes逐步订正法[3,11,13]，完全基于Argo观测剖面资料本身，首先融合形成多年平均气候态温、盐度分布场，进而构建季节平均和月平均初始场，然后对2004年1月以来全球海域的Argo温、盐度剖面资料进行客观分析而构建的。该数据集水平分辨率为1°×1°，其垂向分辩率在0～1975 dbar水深范围内分为58个标准层，其除了温度和盐度两个基本要素外，还包括混合层的相关参数，并且每年更新，目前最新版本的计算结果截止到2019年12月[16]。利用基于梯度依赖相关尺度的最优插值客观分析方法[17-18]，融合太平洋海域(120°E～70°W,60°S～60°N)Argo温、盐度剖面资料构建的《Argo三维网格资料(GDCSM_Argo)》，同样也是立足于Argo观测资料本身研发的1°×1°、逐年逐月的Argo三维网格资料产品。此数据集目前已更新过两个版本，最新结果截止到2017年12月。《热带太平洋海域Argo衍生数据(热、盐含量)产品》和《西太平洋海域Argo衍生数据(混合层、温跃层)产品》则是分别基于不同的温、盐度网格数据产品，采用最大角度法[19]和梯度法[20-21]计算的逐年逐月，空间分辨率1°×1°的混合层和温跃层的各物理参数，这两个数据集也在不定期更新中。

国内许多涉海科研院所也研发制作了多个不同类型的数据集，如，2008年，中国科学院大气物理研究所基于卡尔曼滤波，估算了1999—2010年期间全球海洋的Argo浮标漂流轨迹[22]，制作了全球海洋表层流资料集。对该数据集的时间分辨率分年平均和月平均两种，对应的空间分辨率分别为1°×1°和2°×2°；国家海洋信息中心为了给出更为客观和准确的分析场，使用三维空间的多重网格，通过三维变分数据同化方法[7，23]对2005年1月—2009年12月期间全球的Argo温、盐度剖面资料进行同化分析，于2011年构建了主要包含温度和盐度的逐年逐月的全球Argo网格化产品，该产品水平分辨率为1°×1°，垂向在0～2000 m水深范围内分为26个不等间隔的标准层；2017年，国家海洋预报中心研发的《西太平洋海域Argo资料同化再分析数据集(ROSWPOA)》是基于HYCOM海洋模式[24]，采用集合最优插值同化方案(EnOI)[6]，同化了2005年1月—2015年12月期间西太平洋海域(120°E～180°E，30°S～40°N)的Argo剖面资料、卫星海平面异常和海面温度制作完成的。该数据集的水平分辨率为1/4°×1/4°,垂向在5～1000 m水深范围内分为22个标准层；中科院大气物理所也使用EnOI方法制作了全球1°×1°的月平均温、盐度产品[25-26]，并且为了准确评估OHC，温度和盐度分别采用了WOD数据库1960—2015年和1960—2017年期间的观测资料。

2.2 国外网格化产品

法国Coriolis数据中心为综合构建Argo数据集设计开发了ISAS(In Situ Analysis System)温盐分析系统[27]。该系统以三维最优插值法为理论基础，融合了Argo、XBT、CTD等温盐观测资料，从而进行全球0～2000 m上层海洋的业务化预报，水平分辨率为0.5°×0.5°。目前，该系统已发展到ISAS_V15版本[28]；英国气象局的FOAM(Forecasting Ocean Assimilation Model)业务化海洋预报系统[29]基于最优插值法理论，使用分析订正方法进行四维迭代分析，在北大西洋进行了两组同化试验，同化资料包括实地测量的海表温度、Argo温盐剖面数据、卫星遥感反演的海表异常及海冰资料，实验结果表明在该系统中同化Argo资料对于提高温盐度预报精度是必不可少的，该数据集更新频率为每月一次，目前已经更新到EN.4.2.1版本；日本海洋地球科技厅(Japan Agency for Marine-Earth Science and Technology, JAMSTEC)利用二维最优插值方法，同化Argo浮标、TRITON锚碇浮标和船载CTD仪等观测资料，制作了空间分辨率为1°×1°的全球海洋(包括白令海峡，但除去边缘海)0～2000 dbar的月平均温、盐度数据集[30]；美国夏威夷大学国际太平洋研究中心(International Pacific Research Center, IPRC)使用2005年以来的Argo剖面浮标资料及Aviso卫星高度计资料，利用变分插值技术生成了垂向27层(0～2000 m)，全球2005—2009年气候态及逐年的年平均、季平均和月平均网格化温度、盐度及动力高度资料；印度国家海洋信息中心(Indian National Centre for Ocean Information Services, INCOIS)的客观分析系统是Argo数据处理系统的一部分，用来分析处理印度洋海域2002年以来的Argo浮标温、盐度观测资料[31]。该系统采用的客观分析方法是最优插值的简化形式，即不引入初始场，分析值仅依赖于观测数据，观测值经高斯权重函数加权平均插值到格点上，整个分析过程类似于目标函数的单一迭代，从而退化为最优客观分析过程；而美国Scripps海洋研究所的Roemmich[32]不加入任何其他观测资料，通过线性插值将2004—2008年五年的近35万条Argo浮标剖面观测数据垂向插值到2.5～1975 dbar，构建了空间分辨率为1°×1°，垂向58层的气候态温、盐度分布场，并计算了月平均异常值，给出了完全基于Argo数据的现代上层海洋的一个基本描述，为现行Argo数据与过去数据集作对比提供了基线，表明了Argo数据足以衡量大尺度波动异常，并证实了Argo数据集与其他相关海洋观测资料的相容性。Roemmich等的研究被视为朝着整合Argo资料和其他对气候变化有重要影响的海洋数据迈进的起步阶段。该数据集目前也在持续更新中，更新的频率为每月一次。

这些数据集针对不同的科研问题研制而成，其时间分辨率、空间分辨率、包含的环境要素、涵盖的时间范围和空间范围以及采用的数据同化(分析)方法等，均具有各自的特点。总体而言，大部分数据集仅采用Argo观测剖面研制，且空间分辨率、时间分辨率多数达到1°×1°和月平均及以上。这也再一次证明，Argo剖面资料，尤其是核心Argo计划观测得到的温、盐度剖面数量已足以满足高分辨率网格数据产品制作的需求。这些数据集的制作，70%以上采用了最优插值法(OI)数据同化方法，这或许与OI的编码简单、理论充分、分析结果合理以及计算效率高等优点息息相关。但纵观所有Argo数据产品，同化结果的环境要素仍以温度、盐度及其衍生的热含量、跃层参数等变量为主，生物地球化学观测要素并未囊括其中，这或许是受制于BGC-Argo观测数据的数量与质量问题。尽管如此，目前这些基于Argo观测剖面研发的网格化产品，将Argo散点观测资料进行时空规则均匀化，极大地扩展了Argo资料的应用范围。

3 小结与展望

海洋科学的发展，特别是物理海洋学研究和业务化海洋预测预报技术的发展，长期以来一直受到观测资料不足的制约。随着海洋观测技术的发展，特别是海洋卫星成功发射和全球Argo实时海洋观测网的建立，人们开始有能力获取广阔海洋上大量的、高分辨率的实时海洋观测数据。与此同时，资料同化技术在海洋科学研究中的广泛应用，可以帮助人们从大量的观测资料中提取出更多、更有用的信息，海洋再分析资料集的制作等方面取得的一系列应用成果，有利于深化人们对各种物理海洋现象的认识。

随着国际Argo计划的成功运行和不断扩展，观测要素在核心Argo观测网所积累的温度、盐度、压力(深度)剖面的基础上不断扩充，如BGC-Argo可以观测一系列的生物地球化学参数，包括溶解氧含量、叶绿素a浓度、硝酸盐浓度、pH值、悬浮颗粒物和下行辐照等[33]，这为Argo剖面资料的网格化研究提出了新的要求，主要体现在以下两个方面。

3.1 Argo子计划的剖面观测资料需要进行网格化才能更有效的得到广泛应用

2016年10月正式启动的BGC-Argo计划是第一个有能力监测全球生物地球化学过程的项目[34]，由于BGC-Argo计划的正式实施距今天只有5年的时间，目前BGC-Argo溶解氧观测剖面在全球海域空间分布极不均匀，但在其集中分布的南大洋、印度洋、北大西洋、西北太平洋以及地中海五个试验区，数据密度已能够满足空间分辨率为1°的网格化的基本要求。同时，溶解氧作为海洋生态系统和生物地球化学中的一个基本参数，在海洋新陈代谢[35]、海气通量[36-37]、海气交换过程[38]、净群落生产力[39-41]以及最低溶解氧区(oxygenminimumzones,OMZS)变异[42-44]等问题的研究中都发挥着重要的作用。但相对于散点剖面资料，很多热点科学问题的深入研究更依赖于空间分布均匀的网格数据[45]。而其他生物地球化学要素观测剖面也日益增长，截止到2020年10月底，硝酸盐、叶绿素a、pH值、悬浮颗粒物及下行辐照的观测剖面总数均已超过3万条，且叶绿素a已接近10万条。

3.2 针对新的观测要素的数据同化方法有待于改进发展

不同类型的观测数据，尤其是新的传感器所获取的观测数据，难以采用同一种校正方法，而且校正系统之间也可能存在系统性误差，即使采用同一种方法进行质量检验的数据，其所对应的校正误差也不尽相同，这将给数据网格化过程带来了一定的困难。例如，目前鲜有的包括溶解氧的全球海洋客观分析数据集WOA18，由于其采用的客观分析方法无法合理给出不同系统误差的权重[44]，该溶解氧数据集制作时仅仅采用了利用化学法测得的溶解氧数据。而如前文所述，法国Coriolis数据中心、美国Scripps海洋研究所、日本海洋科技中心和中国Argo实时资料中心等，采用不同的客观分析方法研制的Argo网格化数据产品中，目前也仍主要基于“Core-Argo”观测网所积累的温度、电导率(盐度)和压力(深度)剖面数据构建的网格数据集，尚未包含由“BGC-Argo”观测网所积累的生物地球化学要素，其中一个关键所在即为发展一种合适的数据同化方法。未来，长时间序列的全球海洋剖面数据将主要依靠Argo、Glider等自动平台的组网观测，如何将这些数据有效融合，并改进或突破传统方法，利用机器学习等人工智能方法，构建或重构高时空分辨率的网格化产品，将或是数据同化的未来发展方向。