重力测量地形校正的误差来源（测绘学报李杨）

重力测量地形校正的误差来源（测绘学报李杨）关键词：全球海平面变化 海平面指纹 比容海平面变化 联合反演 时变重力 卫星测高 摘要：理解全球海平面变化具有十分重要的意义，它间接反映了地球系统中气候性相关因素的变化。本文基于一组海平面指纹和比容经验正交函数，联合时变重力数据和卫星测高数据反演了2002年4月至2020年2月的全球海平面变化，将全球海平面变化分解成南极冰盖融化、格陵兰冰盖融化、陆地冰川融化、陆地水储量变化、冰川均衡调整和海水比容效应这6个分量的贡献。联合反演结果显示，全球平均比容海平面变化为1.08±0.05mm/a，与相关文献的结果相吻合。研究发现，联合测高数据和时变重力数据的反演方法能够一定程度上减弱GRACE Follow-On卫星时期海水质量变化被低估的现象。本文利用联合反演的结果研究了区域海平面变化，在大部分近海区域反演效果较好，这表明该方法可用于区域海平面变化的研究。2. 福州大学空间数据挖掘与信息共享教育部

本文内容来源于《测绘学报》2022年第8期（审图号GS京（2022）0490号）

联合时变重力数据与测高数据反演全球海平面变化及其分量贡献

李杨¹

郭金运¹

孙玉² 袁佳佳¹ 常晓涛³ 张红日¹

1. 山东科技大学测绘与空间信息学院，山东 青岛 266590;

2. 福州大学空间数据挖掘与信息共享教育部重点实验室，福建 福州 350108;

3. 自然资源部国土卫星遥感应用中心，北京 100048

基金项目：国家自然科学基金(41774001；41374009)

摘要：理解全球海平面变化具有十分重要的意义，它间接反映了地球系统中气候性相关因素的变化。本文基于一组海平面指纹和比容经验正交函数，联合时变重力数据和卫星测高数据反演了2002年4月至2020年2月的全球海平面变化，将全球海平面变化分解成南极冰盖融化、格陵兰冰盖融化、陆地冰川融化、陆地水储量变化、冰川均衡调整和海水比容效应这6个分量的贡献。联合反演结果显示，全球平均比容海平面变化为1.08±0.05mm/a，与相关文献的结果相吻合。研究发现，联合测高数据和时变重力数据的反演方法能够一定程度上减弱GRACE Follow-On卫星时期海水质量变化被低估的现象。本文利用联合反演的结果研究了区域海平面变化，在大部分近海区域反演效果较好，这表明该方法可用于区域海平面变化的研究。

关键词：全球海平面变化 海平面指纹 比容海平面变化 联合反演 时变重力 卫星测高

引文格式：李杨 郭金运 孙玉 等. 联合时变重力数据与测高数据反演全球海平面变化及其分量贡献[J]. 测绘学报，2022，51(8)：1768-1778. DOI: 10.11947/j.AGCS.2022.20210169

LI Yang GUO Jinyun SUN Yu et al. Inversion of global sea level change and its component contributions by combining time-varying gravity data and altimetry data[J]. Acta Geodaetica et Cartographica Sinica 2022 51(8): 1768-1778. DOI: 10.11947/j.AGCS.2022.20210169

阅读全文：http://xb.chinasmp.com/article/2022/1001-1595/20220810.htm

引 言

全球变暖使地球气候系统产生了一系列变化，例如地表气温升高、海洋热含量增加、海洋浮冰融化和两极冰盖融化。在过去的50年间，温室气体排放使地球气候系统中堆积了大量的热量，海洋吸收了其中93%的热量，剩余7%的热量被土地、冰川、冰盖和大气吸收^[1]。全球平均海平面(global mean sea level，GMSL)变化主要受海水体积热膨胀、南极冰盖融化、格陵兰冰盖融化、陆地冰川融化和陆地水储量(terrestrial water storage，TWS)变化的影响。验潮站数据显示，19世纪中叶以后，全球平均海平面以1.2~1.9mm/a的速率上升^[2-4]。在GMSL上升和局部因素的作用下，沿海地区正受到海平面上升的威胁^[5-7]。

目前，GMSL上升已成为地球气候变化的一个重要响应，研究GMSL有助于更好地观测和解释地球气候变化，发现当今地球观测系统的不足^[1]。政府间气候变化专门委员会在第5次评估报告中指出^[8]：1993年至2010年间使用卫星测高数据估计的GMSL变化为3.2±0.4mm/a，而使用其他数据估算的GMSL变化仅为2.8±0.5mm/a，其中ANT、GRE和GLA的贡献占比为50%，TWS减少导致的海平面上升占比为13%，海水温盐度变化引起的GMSL变化占比为37%，称为比容(Steric)海平面变化。卫星测高数据计算得到的海平面上升速率与其他数据估算的速率存在差异，这种差异可能反映出现今地球观测系统不能完全捕捉到各因素对海平面变化的贡献。

由于负荷自吸引效应，地表或地球内部的质量迁移过程会引起大地水准面变化和地表径向位移，大地水准面变化和地表位移呈现出的不均匀空间特征通常称为海平面指纹(sea level fingerprint，SLF)或静力平衡指纹^[9-11]。随着重力恢复与气候试验卫星(gravity recovery and climate experiment，GRACE)的发射，SLF在海平面变化、地心运动和地球扁率变化等研究中得到了广泛的应用^[12-14]。文献[12 15]提出了一种联合反演方法用于解释全球海平面变化，该方法首先定义一组SLF和比容经验正交函数(empirical orthogonal function，EOF)，SLF和比容EOF分别表示不同的质量迁移过程和海水比容对全球海平面变化造成的影响。基于这组预定义的SLF和比容EOF，联合GRACE数据和卫星测高海平面异常(sea level anomaly，SLA)数据反演GMSL变化，并根据SLF和比容EOF将GMSL变化分解成ANT、GRE、GLA、TWS和Steric的贡献。与传统方法相比^[16]，该反演方法可使用SLF将海水质量变化引起的海平面变化进一步细分成冰盖融化、冰川融化、海陆水循环的贡献；不过分依赖先验的地球物理模型，可以同时估计冰川均衡调整(glacial isostatic adjustment，GIA)对绝对海平面变化的贡献；无须预先扣除GRACE数据和卫星测高数据中的GIA信号；不需要对GRACE数据进行滤波，可以最大限度地保留GRACE数据的原始信号^[12]。

文献[12]实际反演得到的全球平均比容海平面(global mean steric sea level，GMSSL)变化仅为0.17±0.08mm/a，其联合反演结果中存在1.03±0.10mm/a的残差海平面变化无法解释，反演效果不佳。文献[15]利用有限元/体积海冰-海洋模型(finite-element/volume sea ice-ocean model，FESOM)重新计算了一组比容EOF，对GMSL变化进行了再次反演。其反演得的GMSSL变化为1.38±0.16mm/a，该结果与直接使用海洋温盐数据产品计算的GMSSL变化相差较大，其估计的GMSSL变化占据GMSL变化的一半以上，目前普遍认为海水质量对GMSL变化的贡献最大^{[15 17]}。

本研究在文献[12 15]的基础上，使用EN4.2.1海洋温盐产品^[18]计算了一组比容EOF，并根据相关辅助数据解算了一组SLF。基于这组比容EOF和SLF，联合GRACE数据、GRACE-FO数据和卫星测高SLA数据，反演了2002年4月至2020年2月的全球海平面变化，将海平面变化分解成ANT、GRE、GLA、TWS、GIA、Steric的贡献。以期解决文献[12 15]中估算的GMSSL变化与前人相关文献不一致的问题，使用联合反演方法获得合理的结果。

1 数据与方法

1.1 SLF和比容EOF

本文定义SLF的方式与文献[15]定义SLF的方式一致，假设某个地理区域发生质量变化，其造成的大地水准面变化、地表径向位移和相对海平面变化称为SLF。本研究使用的SLF用于区分不同的质量迁移过程，使用所有的SLF同时拟合GRACE时变重力数据或GRACE后继卫星的(GRACE Follow-On，GRACE-FO)时变重力数据可以得到每个SLF在时间域的变化。SLF不在于数量的多少，而在于使用的SLF能否恢复重力场的变化。

本文将地表质量迁移过程分为两类，一类是冰盖和冰川的融化，另一类是TWS变化。研究中把第一类质量迁移过程所在的位置划分成一些地理区域，假设每个地理区域发生单位质量变化，然后对每个地理区域求解一个SLF。根据文献[19—20]对南极冰盖汇水区的定义，将南极冰盖划分成了25个地理区域；根据格陵兰冰盖汇水区的定义并以2000m的高程为界限，将格陵兰冰盖划分成15个地理区域；陆地冰川的位置由全球陆地冰雪空间观测项目提供^[21]，划分了14个陆地冰川区域。TWS变化较为复杂，难以手动划分地理区域，为了计算TWS相关的SLF，本文对Water GAP全球水文模型(water GAP global hydrology model，WGHM)进行经验正交函数分析提取了前60个EOF^[22-23]，这些EOF的累计方差贡献率在90%以上。在提取EOF之前剔除了WGHM水文模型中的冰川区域，对每个EOF求解了一个SLF。水文模型不能捕捉到人为因素引起的TWS变化，为了弥补水文模型和GRACE/GRACE-FO数据的不一致，本文额外引入了3个SLF捕捉里海、北印度平原和维多利亚湖附近区域的地下水减少信号，一共解算了63个TWS相关的SLF。此外，研究中额外定义了一组SLF以考虑GIA对海平面变化的影响：基于ICE-6G冰川历史模型^[24]，解算了2个Fennoscandia区域的SLF，3个Laurentide区域的SLF；基于南极的IJ05冰川历史模型^[25]，解算了1个Antarctica的SLF。使用这6个SLF共同表示GIA对全球大地水准面变化和地表径向位移的影响。

本文使用英国气象局(Met Office)发布的EN4.2.1格网化海洋温盐产品并借助开源程序Gibbs-SeaWater海洋学工具箱计算比容海平面变化^[18]，然后对比容海平面变化进行经验正交函数分析提取比容EOF以参数化海水比容变化对全球海平面的影响。EN4.2.1海洋温盐产品融合了全球海洋数据库13(world ocean database 13，WOD13)、全球温度和盐度剖面项目(global temperature and salinity profile project，GTSPP)、Argo计划和北极概要流域观测项目(arctic synoptic basin-wide observations，ASBO)，提供了1900年至今的海洋三维温盐场信息，相对于仅包含Argo数据的温盐产品，EN4.2.1海洋温盐产品包含的比容海平面变化信息更多。研究中使用了1990年1月至2020年12月的EN4.2.1海洋温盐产品，EN4.2.1海洋温盐产品并未用作观测数据，而是用于提取比容EOF，一共提取了239个比容EOF，其累计方差贡献率在99%以上。使用1990年至2020年的比容海平面变化提取比容EOF能使比容EOF包含的比容海平面变化空间特征更多，例如1997年至1998年的强El Nino事件对比容海平面变化的影响也被囊括在内。

1.2 卫星观测数据

采用美国德克萨斯大学空间研究中心(Center for Space Research，CSR)发布的GRACE/GRACE-FO RL06版本的月时变重力场模型，该模型为截断至60阶次的球谐系数产品，已经扣除了非潮汐大气和海洋信号^[26]。研究中使用了2002年4月至2020年2月共182个月的时变重力场模型，未替换GRACE/GRACE-FO数据中缺失的C₁₀、C₁₁、S₁₁系数和不准确的C₂₀系数，已有研究表明基于海平面指纹的反演方法可以估计这些系数的变化^[13-14]。2016年8月以后，受卫星单加速度计工作模式影响，GRACE/GRACE-FO数据的C₃₀系数不准确，得益于2012年激光相对论卫星(laser relativity satellite，LARES)的发射，卫星激光测距(satellite laser ranging，SLR)数据提供了高精度的C₃₀系数，本文采用戈达德航天中心(Goddard Space Flight Center，GSFC)解算的SLRC₃₀系数替换了2016年8月以后GRACE/GRACE-FO产品的C₃₀系数^[27]，该C₃₀系数与原始GRACE/GRACE-FO的C₃₀系数一致性好，已经被CSR用于改正其生产的Mascon产品。根据这182个月的时变重力场模型计算了研究时间段的平均重力场，将每个月的时变重力场模型都扣除这个平均重力场再乘以平均地球半径以表示大地水准面变化，作为联合反演的输入数据。

本文采用AVISO(archiving validation and interpretation of satellite oceanographic data)提供的2002年4月至2020年2月的格网化月均SLA数据。该数据融合了多颗卫星测高数据，进行了必要的地球物理改正及大气改正，包括电离层改正、干湿对流层校正、固体潮、海洋潮汐、海况偏差和逆气压改正等，该数据的格网大小为0.25°×0.25°。由于时变重力数据在某些月份缺失，为了保证联合反演中时变重力数据和卫星测高数据的参考时间一致，本研究剔除了缺失时变重力数据当月的格网SLA数据，然后计算研究时间段内SLA数据的均值，并从每月的SLA数据中扣除该均值。

本文没有恢复GRACE/GRACE-FO时变重力数据的非潮汐大气和海洋信号，如果恢复该信号，那么在联合反演中，需要额外引入一组基函数表示非潮汐大气和海洋信号。为保持变重力数据和卫星测高数据的信号一致，使用AOD1B(atmosphere and ocean de-aliasing level-1B)模型扣除了SLA数据中的非潮汐海洋信号^{[26 28]}。扣除AOD1B模型之前，每个月的AOD1B模型需扣除其在海洋上的空间均值。非潮汐海洋信号对GMSL变化没有贡献，在研究GMSL变化时可以不恢复该信号，但是该信号对区域海平面变化有影响，在研究区域海平面变化时，必须恢复该信号。考虑到时变重力数据和卫星测高数据的空间分辨率差异，本研究将SLA数据重采样成1°×1°的格网数据，作为联合反演的输入数据。

1.3 联合反演方法

本文所使用的方法称为联合反演方法，其基本假设是造成全球海平面变化的信号在空间上可以较好地分离并建模，并且各类信号造成的空间特征彼此线性无关，不同信号造成的空间特征由SLF或者比容EOF表示。

时变重力数据表示的大地水准面变化是由质量迁移过程导致的，本文认为大地水准面变化是SLF的线性组合，每月可以在频域建立一组时变重力观测方程

(1)

式中，b(t)_grace是GRACE/GRACE-FO数据表示的大地水准面变化；设计矩阵G由SLF组成，每个列向量是一个SLF，表示相应的质量迁移过程对大地水准面造成的影响，设计矩阵G和观测向量b(t)_grace中的元素均为球谐系数且相同行系数的阶次一致；解向量x(t)_SLF表示SLF在时间域的变化。

本文没有填补GRACE/GRACE-FO数据的C₁₀、C₁₁和S₁₁系数，而是剔除了观测方程中的C₁₀、C₁₁和S₁₁系数，在联合反演中估计它们的变化^[14]。考虑到GRACE/GRACE-FO数据的C₂₀系数不准确，还额外剔除了观测方程中的C₂₀系数，C₂₀系数也可以基于海平面指纹反演得到^[13]。尝试了利用联合反演方法估计C₃₀系数的变化，但结果并不理想，因此在数据预处理阶段使用了SLR的C₃₀系数代替了2016年8月以后GRACE/GRACE-FO的C₃₀系数。

绝对海平面会受质量迁移引起的大地水准面变化而变化，还会受到海水温盐度、海洋环流、大气等影响。在文中非潮汐海洋信号已经从SLA数据中扣除，因此，本文认为绝对海平面变化可由预定义的SLF和比容EOF线性表示，每月可在空间域建立一组卫星测高观测方程

(2)

式中，观测向量b(t)_altimetry是SLA数据表示的绝对海平面变化；矩阵B由SLF组成，每个列向量是一个SLF，矩阵C由比容EOF组成，每个列向量是一个比容EOF，设计矩阵A由矩阵B和矩阵C组成；x(t)_EOF表示比容EOF在时间域的变化，解向量x(t)由x(t)_SLF和x(t)_EOF组成，表示所有基函数在时间域的变化。该观测方程中使用的SLF和比容EOF均由格网数据表示，矩阵B、矩阵C和观测向量b(t)_altimetry中相同行元素的经纬度一致。对于格网大小为1°×1°的SLA数据，观测向量b(t)_altimetry有64800个元素，需要注意的是，卫星测高观测方程中应仅含有海洋上的有效格网点值，本研究剔除了陆地上的格网点和格网值为空值的格网点，未剔除近海区域的测高数据。剔除无效数据之后，每月的观测向量b(t)_altimetry约有33000个元素。

为了联立解算时变解，需要给出观测值的权阵，然后根据最小二乘原理求解。设时变重力观测方程中观测值的权阵是P_grace，卫星测高观测方程中观测值的权阵是P_altimetry，其中P_grace是对角矩阵，每个对角线元素与GRACE/GRACE-FO球谐系数的方差成反比，假设权阵P_altimetry是单位矩阵。GRACE/GRACE-FO数据和卫星测高SLA数据是异类数据，无法直接确定它们之间的权重，本研究使用方差分量估计方法确定其权重。在方差分量估计之前，为了保持时变重力观测方程的设计矩阵和卫星测高观测方程的设计矩阵包含的列向量数量一致，需要对时变重力观测方程的设计矩阵G右侧添加239列零向量扩充成矩阵F(239是比容EOF的数量)。每个月的时变重力数据和卫星测高数据之间的相对权重由方差分量σ_grace²和σ_altimetry²确定，方差分量σ_grace²和σ_altimetry²可迭代计算得到，具体算法见文献[29]。

σ_grace²和σ_altimetry²确定后即可联立时变重力观测方程和卫星测高观测方程求最小二乘解

(3)

式中

(4)

(5)

解算得到最小二乘解

(t)后，利用

(t)对SLF和比容EOF进行线性组合即可恢复观测向量或求解某一信号的贡献。例如，重建重力场低阶项系数的变化、全球绝对海平面变化、大地水准面变化等。

2 结果与讨论

2.1 反演结果分析

根据联合反演的结果重建了全球平均绝对海平面(global mean absolute sea level，GMASL)变化。为验证联合反演的结果，本文还根据观测向量b(t)_altimetry计算了GMASL变化，并根据联合反演中观测向量b(t)_altimetry的残差计算了全球平均的残差时间序列，如图 1所示。GRACE/GRACE-FO数据在某些月份缺失，无法进行联合反演，因此本文展示的时间序列中含有缺失值。这些时间序列在2017年7月—2018年5月连续断开11个月，本文未对该时段的缺失值做处理，仅对其他的缺失值进行了线性插值补全。

图 1 全球平均绝对海平面变化Fig. 1 Global mean absolute sea level change

图选项

由图 1可知，残差时间序列具有明显的周期特征，拟合得知其周年振幅为1.29±0.21mm。联合反演和卫星测高数据得到的GMASL变化的线性趋势分别为3.01±0.09mm/a和3.29±0.10mm/a，残差序列的线性趋势为0.28±0.02mm/a。从线性趋势看，当前SLF和比容EOF配置下的联合反演可以恢复超过90%的GMASL变化，仅有一小部分海平面变化无法解释。

联合反演中观测向量b(t)_altimetry残差的标准差如图 2(a)所示。由图可知，联合反演可以重建大空间尺度的海平面变化，但是在一些海域反演的效果不理想。研究发现联合反演残差的标准差与SLA数据的标准差在空间上存在相关性，如图 2(b)所示。这表明联合反演的效果受SLA数据的可变性影响，在SLA变化大的海域，联合反演的效果不好，SLA的高可变性主要与洋流和洋流导致的涡流等海洋现象有关。受日本暖流、千岛寒流和北太平洋暖流的影响，联合反演在日本东部海域的反演效果较差；受墨西哥湾暖流的影响，联合反演在美国东部海域的反演效果较差；受东澳大利亚暖流的影响，联合反演在澳大利亚东部海域的反演效果较差。卫星测高数据在数据预处理阶段使用AOD1B模型扣除了非潮汐海洋信号，由于AOD1B模型本身存在误差且其空间分辨率较低，这样的操作可能影响反演的结果。如果在联合反演中加回AOD1B模型，对AOD1B模型或其他海洋模型提取EOF以表示非潮汐海洋信号，同时估计其时变，有望能减少海洋现象显著海域的反演残差。

图 2 联合反演中残差标准差的空间分布与卫星测高数据标准差的空间分布Fig. 2 The spatial distribution of standard deviation of residuals in joint inversion and the spatial distribution of standard deviation of satellite altimetry data

图选项

为进一步验证方法的可靠性，本文将联合反演得到的C₁₀、C₁₁、S₁₁和C₂₀系数与仅使用式(1)进行反演得到的结果^[13]、GSFC发布的SLR结果^[27]和文献[30]得到的结果(GRACE-OBP方法)进行了对比。C₁₀、C₁₁和S₁₁系数的时间序列如图 3所示，从上到下依次为C₁₀、C₁₁和S₁₁，各时间序列均未加回AOD1B模型。C₂₀系数的时间序列如图 4所示，虚线表示相应的时间序列进行12个月的滑动平均滤波得到的结果，对时间序列进行了上下平移以便辨识，各时间序列均未加回AOD1B模型。

图 3 使用不同方法得到的C₁₀、C₁₁和S₁₁时间序列Fig. 3 Time series ofC₁₀ C₁₁andS₁₁derived from different approaches

图选项

图 4 C₂₀时间序列Fig. 4 Time series ofC₂₀

图选项

由图 3可知，联合反演得到的3个一阶项系数与GRACE-OBP方法的一阶项系数在周年振幅方面比较一致，3种方法的C₁₀和C₁₁系数的线性趋势基本一致，基于海平面指纹反演得到的S₁₁系数的线性趋势均大于GRACE-OBP的结果。C₁₀、C₁₁和S₁₁系数的线性趋势、周年振幅和周年相位见表 1。由图 4可知，联合反演、GRACE-OBP方法和仅使用GRACE/GRACE-FO进行指纹法反演得到的C₂₀系数较为一致，SLR的C₂₀时间序列振幅小于其他3种时间序列的振幅，年际变化基本一致。C₂₀时间序列的线性趋势、周年振幅和周年相位，见表 2。由拟合结果可知，联合反演得到的C₂₀时间序列在线性趋势和周年相位两个方面与GRACE-OBP方法和SLR的结果一致，SLR的C₂₀周年振幅小于联合反演和GRACE-OBP方法得到的周年振幅。

表 1 C₁₀、C₁₁和S₁₁系数的线性趋势、周年振幅和周年相位Tab. 1 Linear trend annual amplitude and annual phase of time series ofC₁₀ C₁₁andS₁₁

系数	方法	线性趋势(10-10)	周年振幅(10-10)	周年相位/d
C10	joint inversion	-0.44±0.01	2.56±0.10	81.54±0.04
	GRACE inversion	-0.43±0.02	2.88±0.12	79.37±0.04
	GRACE-OBP	-0.51±0.02	2.46±0.14	82.52±0.06
C11	joint inversion	-0.09±0.01	1.40±0.08	103.80±0.06
	GRACE inversion	-0.06±0.01	1.45±0.10	99.08±0.06
	GRACE-OBP	-0.08±0.02	1.26±0.06	102.67±0.08
S11	joint inversion	0.14±0.01	1.45±0.10	297.30±0.06
	GRACE inversion	0.19±0.02	1.46±0.10	299.62±0.07
	GRACE-OBP	0.08±0.02	1.42±0.06	286.04±0.08
注：表中的误差估计指最小二乘线性拟合得到的形式误差，95%置信区间。

表选项

表 2 C₂₀时间序列的线性趋势、周年振幅和周年相位Tab. 2 Linear trend annual amplitude and annual phase of time series ofC₂₀

方法	线性趋势	周年振幅(10-10)	周年相位/d
joint inversion	-0.16±0.01	0.82±0.06	75±3
GRACE inversion	-0.15±0.01	0.84±0.07	77±3
GRACE-OBP	-0.16±0.01	0.80±0.09	78±4
SLR GSFC	-0.17±0.01	0.47±0.08	73±6
注：表中的误差估计指最小二乘线性拟合得到的形式误差，95%置信区间。

表选项

2.2 全球平均海平面变化

根据联合反演的结果计算了ANT、GRE、GLA、TWS、Steric对GMSL和GMASL的贡献。GMSL变化及其分量的贡献如图 5所示，划分成6个分量的贡献，拟合并去除了时间序列中的年周期项和半年周期项，进行了6个月的滑动平均滤波，对时间序列进行了上下平移以便辨识。各时间序列的线性趋势、周年振幅和周年相位，见表 3。由图 5可知，冰盖、冰川融化对GMSL变化的贡献较为稳定，TWS和Steric对GMSL变化的贡献波动较大。2004—2014年TWS时间序列的线性趋势为0.05±0.08mm/a，这表明TWS对GMSL变化几乎没有长期贡献，2014年后TWS时间序列出现上升的转折点，这暗示2014年后全球陆地水储量的减少。据表 3可知，冰盖和冰川融化是当今GMSL上升的主要原因，它们对GMSL贡献的线性趋势达到1.68mm/a，受全球气候变暖影响，海水温盐度变化导致的海水体积膨胀是GMSL上升的次要原因。比容作用和TWS变化是主导GMSL周期性变化的主要因素，其周年振幅分别达到了8.66±0.34mm和4.01±0.44mm。残差时间序列的相位为61±3d，与GMSSL的相位(75±4d)接近，海水质量变化在联合反演中被GRACE/GRACE-FO数据约束，从这两方面看残差属于比容信号的可能性较大。

图 5 全球平均海平面变化Fig. 5 Global mean sea level change

图选项

表 3 拟合得到的线性趋势、周年振幅和周年相位Tab. 3 Estimated annual amplitudes phases and linear trends

分量		线性趋势	周年项
		mm/a	振幅/mm	相位/d
ANT	GMASL	0.42±0.01	0.35±0.11	50±2
	GMSL	0.46±0.01	0.38±0.12	50±2
GRE	GMASL	0.65±0.01	0.46±0.07	282±10
	GMSL	0.68±0.01	0.49±0.07	282±10
GLA	GMASL	0.61±0.01	1.00±0.07	291±7
	GMSL	0.64±0.01	1.05±0.07	291±7
TWS	GMASL	0.31±0.05	8.66±0.34	279±2
	GMSL	0.33±0.06	9.21±0.35	279±2
GIA	GMASL	-0.06±0.00	-	-
Steric		1.08±0.05	4.01±0.44	75±4
residuals		0.28±0.02	1.29±0.21	61±3
total GMASL		3.29±0.09	5.73±0.19	305±8
total GMSL		3.47±0.09	6.26±0.23	305±8
注：表中的误差估计指最小二乘线性拟合得到的形式误差，95%置信区间。

表选项

本文除了对1990—2020年的比容海平面变化进行经验正交函数分析提取比容EOF外，还对2000—2020年的比容海平面变化提取了比容EOF，分析了不同比容EOF对联合反演的影响，两组比容EOF的累计方差贡献率均为99%。2条反演得到的GMSSL时间序列与使用EN4.2.1海洋温盐产品计算的GMSSL时间序列如图 6所示。

图 6 联合反演和EN4.2.1海洋温盐产品计算得到的全球平均比容海平面变化Fig. 6 6Global mean steric sea level change derived from joint inversion and EN4.2.1 product

图选项

由图 6可知，相较于使用EN4.2.1海洋温盐产品计算得的GMSSL变化，利用不同比容EOF进行联合反演得到的GMSSL在年际变化和季节性变化方面更加接近。EN4.2.1海洋温盐产品计算得到的GMSSL在2002—2006年的年际变化明显，而同期联合反演解算得到的GMSSL不存在明显的年际变化，Argo计划早期浮标剖面数据少，联合反演方法对GMSSL的估计可能更加合理。

拟合得到了图 6中未经滤波时间序列的线性趋势、周年振幅和周年相位，结果见表 4，两组联合反演得到的GMSSL变化的周年振幅均大于使用EN4.2.1海洋温盐产品计算得到的振幅。受浅海温盐剖面数据少、缺少深层海水温盐信息等缺点，实测海洋温盐数据可能低估比容海平面变化^[31]，本文认为联合反演的比容EOF具备捕捉深海的温盐度变化信息的能力，也能弥补早期剖面数据不足的缺点，因此得到的GMSSL变化应大于直接计算的结果。然而，使用2000—2020年的比容EOF进行联合反演得到的GMSSL变化的线性趋势仅为0.87±0.04mm/a，小于使用EN4.2.1海洋温盐产品计算的结果，比容EOF并没有从卫星测高数据中捕捉到足够的比容海平面变化信息，这说明2000—2020年的比容EOF包含的比容信号过少，使用1990—2020年的比容EOF得到的结果更合理一些。

表 4 全球平均比容海平面变化的线性趋势、周年振幅和周年相位Tab. 4 Linear trend annual amplitude and annual phase of global mean steric sea level change

GMSSL	线性趋势/(mm/a)	周年振幅/mm	周年相位/d
EOF(1990—2020年)	1.08±0.05	4.01±0.44	75±4
EOF(2000—2020年)	0.87±0.04	3.18±0.35	82±5
EN4.2.1	0.97±0.06	2.33±0.60	49±1
注：表中的误差估计指最小二乘线性拟合得到的形式误差，95%置信区间。

表选项

使用1990—2020年的比容EOF反演得到的GMSSL变化的线性趋势为1.08±0.05mm/a，大约占GMSL变化的三分之一，相比于文献[12]和文献[15]研究估计的GMSSL，本文结果与相关文献的记载的结果更符合^[32-33]。联合反演方法对GMSSL变化的估计很大程度上取决于先验比容EOF包含比容信号的多少。文献[12]首次反演得到的GMSSL变化仅为0.17±0.08mm/a，本文认为是因为其使用ishii海洋温盐产品计算0~700m的比容海平面变化导致的，这使得比容EOF包含的比容变化信号过少，直接导致了比容EOF无法从测高数据中捕捉到完整的比容海平面变化。文献[15]使用FESOM模型提取了新的比容EOF用于反演，得到的GMSSL变化为1.38±0.16mm/a，超过了海水质量变化的贡献(1.06mm/a)，本文认为这是由于其解算的比容EOF包含一部分噪声信息，噪声信息可能与SLF存在一定的相关性，导致反演时错误地估计了海水比容效应和海水质量变化。

在GMSL变化研究中通常使用先验的GIA模型扣除GRACE数据和卫星测高数据中的GIA信号，或扣除一个大约-0.3mm/a的GIA改正值将GMASL转换成GMSL^[34-35]。本联合反演得到的GIA对GMASL的贡献仅为-0.06mm/a，一方面这是因为研究中使用的南极地区的GIA SLF是基于IJ-05冰川历史模型解算得到，相比于ICE系列模型，IJ-05模型估计的南极GIA信号对现今GMASL贡献的量级较低^[36]，另一方面，是因为本研究解算得到的南极区域GIA信号量级比先验的IJ-05模型量级小。这一结果与文献[15]的研究类似，其估计的GIA贡献为-0.1mm/a，虽然使用的SLF不同，但是在联合反演中解算得到的南极GIA信号都较低。

利用联合反演结果计算了全球平均海水质量(global mean ocean mass，GMOM)变化，作为对照，本文依据式(1)仅使用GRACE/GRACE-FO数据反演了GMOM变化，还间接地基于测高数据和海洋温盐数据计算了GMOM变化，结果如图 7所示。2006—2016年间，测高数据减比容数据估算的GMOM变化与反演得到的GMOM变化基本一致，残差都较小。2016年8月后，受卫星单加速度计工作模式的影响，仅使用GRACE/GRACE-FO时变重力数据的反演可能低估了GMOM变化，这一现象与文献[37]吻合。值得注意的是，融合了卫星测高数据的联合反演在一定程度上减弱了GMOM变化被低估的现象。在GRACE-FO数据时期，将仅使用时变重力数据反演的GMOM变化与基于测高数据和比容数据估算的GMOM变化相减，其残差的RMS为13.60mm，而联合反演估计的GMOM变化与基于测高数据和比容数据估算的GMOM变化相减，其残差的RMS仅为5.78mm。

图 7 全球平均海水质量变化(以等效海平面变化表示)Fig. 7 Global mean ocean mass change(expressed as equivalent sea level change)

图选项

2.3 区域海平面变化

本文使用的联合反演方法理论上可以分解海洋上任意位置的海平面变化，上文已经展示了在近海区域联合反演残差的标准差没有明显增加的特点，这表明该方法不仅适用于GMSL变化的研究，还可用于区域海平面变化的研究^[38]。本文试验恢复了非潮汐海洋信号，选取了20个海域计算区域的相对海平面变化，并将海平面变化分解成ANT、GRE、GLA、TWS、GIA、Steric、Other的贡献，Other指的是联合反演的残差和AOD1B模型的贡献，结果如图 8所示。

图 8 各区域的相对海平面变化Fig. 8 Relative sea level change in some regions

图选项

不同于GMSL变化，区域相对海平面变化还会受GIA影响，近海的相对海平面变化通常还受地质构造、地震、陆地沉降等影响^[39]。本文并未顾及这些区域信号，但可在联合反演的后处理中加上这些信号。由图 8可知，虽然海水比容效应是GMSL变化的次要因素，但在很多近海区域，如墨西哥西部海域、印度西部海域、塔斯曼海等，比容变化的贡献更大，超过了海水质量变化的贡献。

本文得到的中国南海相对海平面变化为4.33±1.81mm/a，其中海水比容变化的贡献占1.38mm/a，ANT、GRE和GLA冰融化的贡献分别为0.47、0.79和0.69mm/a，TWS的贡献为0.34mm/a。中国东海区域的反演效果不佳，这可能与中国近海Argo剖面数据较少有关。在北海、欧洲西部海域和地中海联合反演的质量最差，残差变化的线性趋势分别为2.48、2.12和1.40mm/a，澳大利亚东部海域和塔斯曼海的反演效果最好，其残差趋势仅为0.02mm/a和0.12mm/a。由图 8中的20个海域反演情况来看，仅有欧洲的3个海域和中国东海反演效果较差，相比于GMSL变化0.28mm/a的残差，联合反演中其他16个海域的Other贡献在可接受的范围内，表明本研究结果可以得到对GMSL变化和区域海平面变化的一致性估计，可适用于近海的海平面变化研究。

3 结论

本文基于一组预定义的SLF和比容EOF，联合2002年4月—2020年2月的时变重力数据和卫星测高数据反演了全球海平面变化，在最小二乘意义上给出了ANT、GRE、GLA、TWS、GIA和Steric对GMSL和GMASL变化的贡献。

结果显示，本文解算的SLF和比容EOF可以解释90%以上的GMASL，2002—2020年GMASL和GMSL的线性趋势分别为3.29±0.09mm/a和3.47±0.09mm/a。联合反演得到的GMSSL的线性趋势为1.08±0.05mm/a，与相关文献中的结果一致，比文献[12]和文献[15]估计的GMSSL变化合理。研究发现，结合卫星测高数据和时变重力数据的联合反演方法能够减弱GRACE-FO卫星时期海水质量变化被低估的现象。本文分析了联合反演中卫星测高观测方程中观测向量残差的标准差，发现联合反演在海洋动力现象显著的海域效果较差，在近海区域反演效果并没有降低。基于此发现，利用该联合反演方法分解了20个区域的相对海平面变化，将该联合反演方法用于区域海平面变化的研究。

作者简介

第一作者简介：李杨(1995—)，男，硕士生，研究方向为大地测量学。E-mail：isyang.li@outlook.com

通信作者：郭金运 E-mail：jinyunguo1@126.com

初审：张艳玲复审：宋启凡终审：金 君

资讯