刘磊， 蒋星亮， 李志锦， 费建芳， 吴新荣， 王辉赞，叶芳， 张福颖， 史文丽

1. 国防科技大学 气象海洋学院, 长沙 410073;

2. 加州大学洛杉矶分校 区域地球系统科学与工程联合研究所, 洛杉矶 90095-7228;

3. 复旦大学 大气与海洋科学系&大气科学研究院, 上海 200438;

4. 南京信息工程大学, 南京 210044;

5. 自然资源部海洋环境信息保障技术重点实验室, 天津 300171;

6. 陆军炮兵防空兵学院南京校区, 南京 210000

1 引 言

海面高度是气候变化最重要的指标之一，反映了气候系统对人为强迫以及自然变化的响应（Abraham等，2022）。20世纪初开始，通过对实况潮汐观测发现全球海平面高度GMSL（Global Mean Sea Level） 持 续 以 每 年1.7±0.3 mm 的 速 度 抬 升（Jevrejeva等，2008；Woppelmann等，2009；Church等，2011；Stocker等，2013）；从90年代开始，高精度的卫星高度计探测手段在海表高度场测量中的使用，为海洋分析和预报提供了全天候、高分辨率的海表高度场观测数据（Lellouche 等，2013；Oke 等，2015）。

过去二十多年，全球范围内很多研究团队都在关注于卫星高度计资料的处理并提供全球或区域尺度的海表高度时间序列产品（Greenslade 等，1997；Le Traon 和Dibarboure，1999）。与卫星高度计的沿轨数据相比，在时间和空间尺度上具有连续性的多卫星融合产品能够直观的表现海表地转流和海洋动力过程，具有很高的实际应用价值（Koblinsky 等，1992；Ubelmann 等，2016）。卫星高度计融合产品的质量主要受以下3 个方面的影响：（1）制作产品所选择的卫星种类：因不同种类卫星的轨道高度、周期差异，导致获取的海表高度时空精度也存在不同（Pascual 等，2006；Dibarboure 等，2011）；（2）高度计标准的选择：包括高度计L2P资料处理中的仪器参数、地球物理参数、环境校正、平均海面的构造等（Dibarboure等，2011；Ablain 等，2015）；（3）融合处理资料的方法：不同融合方法在背景场的选择、背景误差相关系数尺度的差异以及观测误差的处理等（Ducet 等，2000；Dussurget 等，2011；Griffin 和Cahill，2012；Escudier等，2013）。

针对现有卫星高度计融合产品的研发现状，本文系统地介绍了当前卫星高度计融合产品的种类，对比分析产品质量，并针对卫星高度计的融合方法，着重介绍卫星高度计融合产品质量的提升途径，进一步强调融合方法与有效分辨率的关系，为提升卫星高度计融合产品质量，提高融合产品的有效分辨率，扩大卫星高度计融合产品的应用范围提供参考。

2 卫星高度计融合产品现状

1991 年和1992 年发射的ERS-1 和TOPEX/Poseidon（TP）开启了卫星高度计高精度测量的历史。随着卫星自身高度准确定位技术取得突破，TP 卫星的定位准确度从米量级提高到厘米量级，最终对海面高度测量的精度达到了几厘米，从而使得卫星高度计资料在海洋中尺度涡等研究中得到广泛应用（Ubelmann 等，2016；Ablain 等，2017）。表1 给出了国内外高度计卫星的发展情况（http：//marine.copernicus.eu/documents/QUID/CMEMS-SLQUID-008-032-062.pdf［2020-07-03］）。

表1 高度计卫星的概况Table 1 Basic information about altimeter satellites

Le Traon等（2003）最早融合了GEOSAT后续卫星、T/P 和ERS-2 这3 颗卫星高度计资料。2001 年12 月，替代T/P 高度计卫星的Jason-1（简称J1；Jason-2 简称J2；Jason-3 简称J3）高度计卫星升空，出现了4颗卫星同时在轨运行时期。众多科学家（Fu 等，2003；Le Traon 和Dibarboure，2004；Pascual 等，2006）从研究角度提出多颗卫星高度计同时在轨时，得到的融合产品质量显著提高。但上述研究处于试验分析阶段，并没有形成多源卫星高度计融合产品。

1997 年AVISO （Archiving，Validation，and Interpretation of Satellite Oceanographic data） 的 高度计数据融合系统DUACS （Data Unification and Altimeter Combination System）开始制作发布在国际上影响最为广泛的卫星高度计产品，包括延迟几小时到一天的准实时NRT（Near Real Time）产品和延迟若干月发布的后处理产品DT（Delayed Time）。为了尽快获取多源高度计数据，NRT 产品采用非中央处理时间窗口，仅处理窗口时间前六周的高度计数据，主要应用在业务预报领域；DT产品利用所有在轨运行的卫星高度计数据，得到经交叉校准的高精度长时间序列的卫星高度计融合产品，主要应用在科学研究领域。

2015 年5 月起该系统由欧洲哥白尼计划（European Copernicus Program）负责，DT2018融合产品分别由哥白尼海事局CMEMS（Copernicus Marine Service）和哥白尼气候服务机构C3S（Copernicus Climate Change Service）生产和发布。其中，CMEMS产品主要针对全球或区域尺度范围内海洋模式和海洋环流分析中的中尺度信号反演，采用所有高度计卫星数据（图1），具有更为精确的海表高度精度（Taburet 等，2019）。C3S 产品主要用来监测海表高度的长期变化，分析气候响应因子，因此选用了具有稳定性能的两个高度计卫星数据（T/P，J1-3 系列与CMEMS 产品中的另一颗卫星融合），上述产品的具体差异可以参阅Taburet 等（2019）。表2给出了卫星高度计融合产品（全球）的概况。

图1 DT2018采用的卫星高度计数据时间线（Taburet 等，2019）Fig. 1 Timeline of modern radar altimetry missions mapped by DT2018 （Taburet et al.， 2019）

DUACS 每隔4 年（会适当调整时间）会依据最新的处理方法更新DT产品，本文仅介绍DT2010之后的版本。当前最新的DT2018 产品由CMEMS发布，采用了累计76 个使命年（commission years）共12 种不同的高度计数据（Taburet 等，2019），制作了从1993 年1 月1 日开始超过了26 年时间序列的高度计融合产品（图1）。产品要素包括绝对动力高度ADT（Absolute Dynamic Topography）融合场和沿轨ADT 数据，同时也发布了地转流的绝对融合场和异常融合场。ADT 场由平均动力高度MDT（Mean Dynamic Topography）和海表高度异常SLA（Sea Level Anomaly）组成，我们将在第3 节中介绍ADT 场的划分原因。为提高产品质量，DT2018 主要从以下3 个方面进行了改进：（1）增加了DT2014 未考虑的卫星种类，如在部分地区增加了2016 年3 月—2017 年2 月的HY-2A 卫星资料（林明森 等，2019）；（2）采用了新的高度计标准和地理修正，原始数据的选择进一步优化，最优插值参数区分了全球和区域范围；（3）相对DT2010和DT2014 的SLA 计算中MSS 统计场分别采用7 年和20 年平均，DT2018 进一步提升到了25 年平均，提高了SLA的计算准确性（Taburet 等，2019）。

表 2 卫星高度计融合产品（全球）概况，当前在轨运行卫星标粗Table 2 Introduction about merged satellite altimeter products （global）， and the bold are running satellites at present

经上述改进，采用独立沿轨数据评估后发现，DT2018 产品中尺度部分的均方根误差（RMSE）相对DT2014 减小了3%—4%，尤其在沿岸区域更为显著，达到了10%，地转流（Geostrophic current）的均方根误差减小了5%到10%。DT2018 产品相对DT2014 有了显著提升，但仍然存在以下不足：沿轨资料处理采用的低通滤波以及格点产品制作过程中使用的最优插值方法，使海表高度场上波长在200 km 到65 km 范围内的波动能量损失60%，即3/5 的中尺度信号缺失（Pujol 等，2016）。中纬度地区的水平分辨尺度为1.7°，平均200 km，使得产品无法有效分辨更小尺度的中尺度信号（Chelton 等，2011和2014）。

针对现有卫星高度计融合产品中对海表高度场中尺度信号缺失的现状，本文着重从融合方法上分析原因。

3 高度计数据获取与处理

在DUACS 制作DT 系列高度计产品的流程中，首先选定计算SLA 的高度计标准，包括仪器订正，地理参数、环境校正、平均海面MSS（Mean Sea Surface，）等。随着研究深入，标准不断进行调整（见表1），之后继续完成（Taburet 等，2019）：高度计数据获取与均质化处理、输入数据质量控制、多星交叉校准、沿轨SLA 数据生成，最后生成融合产品并进行最终质量控制等多个步骤（图2）。

图2 DUACS系统处理流程（Taburet 等，2019）Fig. 2 DUACS system processing（Taburet et al.， 2019）

在高度计数据获取环节，DT 产品获取自1993 年起若干年份的不同卫星种类的可用高度计数据，随着新的高度计卫星升空，DT 也相应的补充了部分相关数据（详见表2），然后参照各研究机构和研究团队选定高度计标准，对数据进行均质化处理（Taburet 等，2019）。

在输入数据质量控制中，为了检测出高度计观测中的错误虚假值，通过改变阈值进行质量控制（Ablain 等，2010）。如在DT2014 产品的数据处理中发现，高纬地区以及离岸20 km区域MSS分辨率不高，导致非重复轨道卫星和新到达该地区的卫星的沿轨SLA 出现异常值（Pujol 等，2016）。为避免该问题出现在DT2014 产品中，在这些关键区域严格筛选了错误检测标准（Pujol 等，2016）。DT2018 中，MSS 数据质量在沿岸减小的问题得到解决，因此，在海岸附近尽可能地保留了有效测量数据（Taburet 等，2019）。

在数据获取环节中，从不同高度计卫星获取的高度场已经均质化处理，但由于仪器误差、轨道残差以及处理误差等存在，仍需要多星交叉校准（Ablain 等，2015；Fu 和Haines，2013）。该环节对于气候信号具有很高的重要性，需要确保TP、J1 和J2 卫星数据中MSS 一致，之后进行多卫星间的轨道误差校正。郭靖等（2013）对Haiyang-2A（HY-2A）轨道进行校正，定轨精度可达1—2 cm。DT 系列（包括DT2010、DT2014 和DT2018）产品制作中的校正方法保持一致，具体可参阅Le Traon和Ogor（1998）。为了减小相邻轨道间地理位置相近引起的观测误差相关性，采用长波误差减小算法，引入长波相关误差ELW（Pujol 等，2016）。

海表高度SSH（Sea Surface Height）为卫星相对地球重力中心轨道高度减去高度计距离海面的高度，由大地水准面高度（Geoid）和ADT 组成（图3）。其中，Geoid 是等重力位势面，如果海洋处于静止状态—即仅在重力场下没有洋流，它将与海洋表面相对应（Nerem 等，1994；Rapp 等，1996）。ADT 是与地转流等动力结构密切相关的海洋高度场。当海洋受到风、热交换、降水以及其他动力因素的影响时，海洋表面会偏离Geoid 面（Kelly和Gille，1990；Qiu，1994）。因此，Geoid面的变化提供了海洋动力学的信息。过去，由于准确获取地球Geoid 面存在较大困难，为了得到具有动力意义的高度场信息，AVISO 利用高度计卫星轨道的重复特性，通过SSH 减去经多年时间平均计算的MSS 后（如DT2018 采用25 年平均），得到SLA（Taburet 等，2019）。当高度计卫星为非重复轨道时，无法采用同轨道的时间平均计算MSS，SLA的计算采用格点化的MSS作为参考面（Ophaug等，2021）。近年来，随着GRACE和GOCE等重力卫星的升空（郑伟 等，2010；冯伟 等，2012），获取较高测量精度的Geoid 面成为可能，MSS 减去Geoid 后可获得MDT，SLA 与MDT 相加，即可得到ADT。

图3 高度计测高原理（引用自AVISO手册：http：//marine.copernicus.eu/documents/PUM/CMEMS-SL-PUM-008-032-062.pdf）Fig. 3 Altimetry principle （Cited from AVISO：http：//marine.copernicus.eu/documents/PUM/CMEMS-SL-PUM-008-032-062.pdf）

SLA受到观测噪音的影响，在沿轨SLA数据的制作中，采用Lanczos 低通沿轨滤波减小观测噪音，并随纬度不同而调整截断波长（Le Traon 和Dibarboure，1999）。DT2018 沿轨产品的生成方法与DT2014 保持一致，只是采用了新的高度计标准和新的卫星高度计数据等。通过信噪比分析（Dufau 等，2016），沿轨滤波数据的截断尺度平均为65 km，表明沿轨滤波数据可有效分辨65 km 以上波长尺度的海洋结构。滤波方法降低了沿轨产品的分辨率，为了在模式和预报系统里的应用而不降低分辨率，从DT2014 开始增加了沿轨非滤波产品，保留了沿轨1 Hz的短波信号。

生成沿轨滤波数据和沿轨非滤波数据后，即进入数据融合环节。采用融合方法对多源卫星高度计沿轨数据进行插值，在规则网格上产生多源卫星高度计融合产品（刘磊 等，2020）。不同的融合方法对产品质量是否有影响，在具体的资料融合中该如何选择？本文给出具体分析和介绍。

4 融合方法与背景场

DUACS 发布的DT 系列融合数据采用基于高斯—马尔科夫理论的最优插值（Bretherton 等，1976）：

图4 DT2018产品中纬向相关长度尺度（Pujol 等，2016）Fig. 4 Zonal correlation scales in DT2018（Pujol et al.，2016）

Bretherton 等（1976）的最优插值估计广泛的应用在海表温度和海表高度的资料融合中（He等，2003；Guan 和Kawamura，2004；Barron 和Kara，2006；Kawai 等，2006）。在国内，徐广珺等（2013）利用最优插值反距离加权法融合了HY-2A、J1/2 以及ENVISAT 卫星高度计数据，验证了HY-2A 在资料融合中可替代ENVISAT 数据；崔伟（2016）采用Bretherton 等（1976）的最优插值方法融合了ENVISAT、J1、J2以及HY-2A等4颗同时在轨运行的卫星高度计测高数据，得到融合产品质量与DUACS产品基本一致。

高度计卫星获取的SSH 场为准地转平衡动力场，可以由线性的最优插值方法将大部分的中尺度动力信息融合网格化（Ducet 等，2000），但非线性的小尺度信息，不能被最优插值方法识别并融合。为减小插值误差，Ubelmann 等（2015 和2016）提出了动力插值方法，加入准地转位涡方程，并应用在高度计沿轨资料SLA 的网格化中，估计场由最优插值方法表示为：

式中，xa是估计场，yo是观测场，x为分析场，εN为测量误差，N为观测数量，B为背景误差协方差矩阵，Rs为观测误差协方差矩阵，Hs为观测算子，将网格点投射到观测点，T为矩阵转置。由于应用在SLA 融合中，因此式（4）没有列出背景场。为了考虑SSH的非线性传播，引入准地转位涡方程：

由于M为非线性算子，无法由y0反向线性拟合得到xtref，所以给出一个初猜场，初猜误差为准线性：

基于相同的数理基础，变分方法推导的式（24）与最优插值式（1）具有相同的形式。但是，二维变分方法相对最优插值具有以下计算优点：（1）在整个网格点上进行最小化，计算效率高，可以解决高分辨率条件下格点数巨大的求解问题；（2）背景和观测误差协方差矩阵灵活多变，可更方便的简化处理，同时有利于引进动力约束；（3）依据不同的观测算子Hs，方便融合不同性质的观测资料。

在卫星高度计的融合方法中，依据选定的背景场，构建与该背景场对应的背景误差协方差矩阵。二维变分方法计算背景误差相关系数时，以高斯函数构建背景误差协方差矩阵，用最小二乘法进行高斯分布拟合得到的背景误差相关系数尺度为约50 km（图5），远小于DT2018。融合方法中背景误差相关系数尺度是影响融合产品质量的重要因素，而背景误差相关系数尺度又决定于背景场的选择；背景场的选取决定了背景误差协方差，从而影响了背景误差协方差相关系数尺度（Li 等，2015 和2016）。因此，DT2018 的背景误差相关系数远大于选用前一天的融合产品结果作为背景场的二维变分方法（图5），其中，背景误差以相应的融合场与沿轨滤波资料相减得到。红实线为2DVAR 计算得到的相关系数分布；红虚线为用最小二乘法拟合高斯函数分布（相关系数尺度为50 km）；蓝实线为AVISO 的DT2018 融合产品25 年气候平均背景场计算得到的相关系数分布（刘磊 等，2020）。

图5 背景误差相关系数Fig. 5 Background error correlation coefficients

5 有效分辨率的提升途径

融合产品的有效分辨率是融合产品质量的重要指标，即为可分辨海面涡旋的最小尺度，定义为融合产品与沿轨资料能量谱密度的比值为1/2 的尺度（Chelton 等，2007）。在多源卫星高度计资料的融合中，如何在现有在轨运行卫星高度计数量维持不变的前提下，进一步提升有效分辨率，是众多学者关注的热点问题。

Chelton 和Schlax（2003）、Chelton 等（2011 和2014）通过计算SSH 的融合误差或融合产品与沿轨数据的能量谱密度的比值，估计了早期SLA 融合产品的有效分辨率。分析结果认为中纬度的有效分辨率在2°—6°之间，主要受融合产品中的高度计数量以及采样轨道的影响。

Ballarotta 等（2019）详细介绍了DUACS 的DT系列融合产品有效分辨率的变化，并得出DT2014产品在强变化区域的有效分辨率相对DT2010 产品提升2%，如湾流、黑潮和南极绕流海域。有效分辨率的提升主要在于改进了高度计卫星探测技术、大气指数修正、潮汐指数修正、校正方法以及融合方法中采用更小的相关系数尺度等（Pujol 等，2016；Taburet 等，2019）。在赤道附近有效分辨率降低可能与该区域更长的相关系数尺度有关（Pujol 等，2016）。DT2018 的有效分辨率总体上与DT2014 一致，但在强变化区域有2%—10%的提升，改进的原因主要在于融合方法中观测值的优化和选用更长统计周期的背景场，由此得到更为准确的背景误差（Ballarotta 等，2019）。另外，融合产品中高度计卫星选用的数量也对有效分辨率有一定影响。在DUACS 产品中，融合高度计卫星的数量由两颗提升到3颗时，有效分辨率平均增加5% （Ballarotta 等，2019）。Ballarotta 等（2019）认为，DUACS 产品有效分辨率随卫星数量的增加而增幅不明显的原因在于采用的最优插值方法中，误差协方差参数没有变化，对融合场施加了较强的滤波效果。

Li 等（2015）指出：融合产品的有效分辨率决定于背景误差相关系数尺度，与观测信息由观测点插值到规则网格点的信息量密切相关。当背景误差相关系数尺度给定时，无论观测资料密度如何变化，融合产品的有效分辨率都不高于两倍相关系数尺度。因此，可以从以下几个方面提升融合产品的有效分辨率：（1）提升高度计观测密度；（2）考虑同轨观测误差相关性；（3）依据选取的背景场准确估算背景误差；（4）构建更为合理的背景误差相关系数表示函数，考虑不同区域的背景误差相关系数尺度变化。最终使资料融合中背景误差相关系数尺度与观测密度相匹配。在一定的观测密度分布下，背景误差相关系数尺度小，融合产品的有效分辨率高，相反，背景误差相关系数尺度大，融合产品有效分辨率低（Daley，1991）。从理论上讲，如果选用的背景场合理，包含更多的大尺度信号，从而将更少的大尺度信号归于背景误差中，可显著降低背景误差协方差相关系数尺度，有效提升多源卫星高度计融合产品有效分辨率。依据上述原则选用前一日融合场作为背景场的2DVAR 方法，可显著提高融合产品的有效分辨率，见图6：AVISO （蓝） 和2DVAR（红）融合产品被插值到图6（a）3（Jason-3）和图6（b）S3A 轨道（黑）后的平均沿轨能量谱密度对比；AVISO（蓝）、2DVAR（红）与图6（c）J3 以及图6（d）S3A 平均沿轨能量谱密度的比值。由于卫星沿轨滤波数据滤除了65 km 以下波动，没有画出小于65 km尺度部分（刘磊 等，2020）。

图6 平均沿轨能量谱密度比值Fig. 6 Average energy spectral density and its ratio

6 结 论

卫星高度计融合产品广泛应用在海洋环境监测、海洋中尺度系统研究以及海洋数值预报等业务领域，具有极高的科学和社会价值，为人类更好认识海洋、研究海洋提供了新的工具。卫星高度计融合产品的研发涉及数据采集、数据均一化处理、数据质量控制、多卫星交叉校准、沿轨SLA数据生成、多源卫星融合以及最终的产品质量控制等多个流程（Pujol 等，2016），是众多科学家跨专业通力合作的成果。如何立足现有观测密度不变的条件下，最大程度地提取有效观测信息，以改进融合产品的有效分辨率和进一步提升产品质量，是海洋学家们研究的热点和当前亟需解决的难点问题，需对当前卫星高度计产品的研究现状进行梳理。

本文以国际上广泛应用的AVISO 产品为主线，详细介绍了DUACS制作的DT系列融合产品。从融合方法入手，分析了不同融合方法对产品质量的影响，进一步阐述了融合方法中背景场的选择对融合产品有效分辨率的决定作用。结合近期已发布数据的Sentinel-3A（S3A）和J3 卫星，以及CryoSat，SARAL/AltiKa 和HY-2A 等卫星，卫星高度计沿轨观测密度显著提升，而DUACS 制作的DT2018 产品的有效分辨率相较DT2014 和DT2010，并没有显著提升，主要原因在于背景场依次选用了7 年、20 年和25 年平均的MSS 统计场，将更多的大尺度信号包含于背景误差中，导致由此统计得到的背景误差协方差相关系数尺度更长，对观测资料的滤波效果更为显著，有效分辨率无法有效提升。

随着Sentinel-3B（S3B）进入运行阶段，HY-2B和HY-2C卫星也分别于2018年10月和2020年9月成功发射（Guo，2018），结合即将于2021 年发射的SWOT 高度计卫星（徐永生 等，2017），卫星高度计观测数据的精度和密度会跨上一个新的台阶。如何在观测中考虑同轨观测误差相关性，如何准确估算背景误差，如何建立更为合理的背景误差相关系数表示函数，为进一步提升卫星高度计融合产品质量和有效分辨尺度提供有效手段及方法改进，是科学家们继续探索的方向。