海平面在现今气候变化中扮演重要角色，比如其变化可以预测地表温度变化(Peyser et al., 2016).研究表明，自20世纪以来，全球平均海平面以1.3~1.9 mm·a^-1的速率上升(Calafat et al., 2014; Hay et al., 2015)，进入20世纪90年代后，海平面有加速上升的趋势(Watson, 2016).海平面加速上升会导致严重的后果，比如全球海岸线的城市面临被淹没的危险，会对全球经济造成极大破坏，直接威胁人类的生存和发展(Nicholls and Cazenave, 2010).因此，确定海平面上升的速率有着重要的科学价值和现实意义.

自20世纪90年代以来，卫星测高显示全球平均海平面的上升速率大约为3.1 mm·a^-1(Cazenave et al., 2014).值得注意的是海平面上升具有明显的空间差异，比如西太平洋，特别是菲律宾海域的上升速率达到了15 mm·a^-1，而美洲沿岸的太平洋区域却呈现出3~5 mm·a^-1的下降趋势(Nicholls and Cazenave, 2010).中国近海的海平面变化略高于全球平均值，其中南海的海平面上升速率约为4.2 mm·a^-1(郭金运等，2015).需要指出的是，以往的研究在确定海平面上升速率时是通过最小二乘拟合获取线性速率，该速率值表示研究时间段内的海平面平均速率变化.事实上，作为一种非稳态变化，海平面上升的速率是随着时间变化的，比如利用滑动窗口拟合的线性速率表明，全球平均海平面的速率在2010年后有显著的上升(Yi et al., 2015).线性速率的估计不仅受到时间段的影响，而且也受到海平面年际变化的影响(Cazenave et al., 2014).比如，滑动窗口拟合的结果表明，全球平均海平面的速率从2002年的3.5 mm·a^-1下降至2008年的2.5 mm·a^-1左右，而在经过海平面质量的年际变化改正后，这种速率下降的变化基本消失(Cazenave et al., 2014).这表明改正或者分离海平面的年际变化，对于认识海平面上升的速率十分重要.

本文使用集成经验模态分解(Ensemble Empirical Mode Decomposition, EEMD)(Wu and Huang, 2009)来分离得到海平面变化的非线性趋势部分，然后计算其瞬时速率变化.作为数据自适应的方法，集成经验模态分解无需假定海平面是以线性速率(或者伴随加速率)变化，而是将海平面变化时间序列分解为本征模态函数，特别适用于非稳态的时间序列.事实上，海平面变化就是一种非稳态的时间变化.与以往最小二乘拟合或者滑动窗口拟合得到的速率不同，我们将计算该非线性趋势的瞬时速率，即每个时间节点处的速率值.本文将分析2002—2014年全球平均海平面上升的瞬时速率，并讨论其各个贡献成分的瞬时速率，比如格陵兰岛和南极冰川消融导致的海平面上升，以及比容海平面变化等.

本文使用的数据是由Rietbroek et al.(2016)提供的海平面时间序列(2002年4月至2014年6月).该数据同化吸收了卫星测高的沿轨数据、海洋模型的比容数据以及GRACE(Gravity Recovery and Climate Experiment)卫星(Tapley et al., 2004)的观测数据，强制海平面方程闭合，得到了总体海平面变化，比容海平面变化，格陵兰岛、南极和山地冰川消融对海平面变化的贡献，陆地水储量变化对海平面的贡献，冰川均衡调整对海平面的贡献以及剩余的非随机信号部分(图 1).由于冰川均衡调整对海平面变化的贡献是线性的(在本文研究十年尺度上，认为冰川均衡调整是线性变化)，Rietbroek et al.(2016)给出的结果为-0.09 mm·a^-1.我们直接将冰川均衡调整的线性速率从总体海平面直接扣除，不作深入讨论(线性变化的瞬时速率是一个固定值).需要强调的是，剩余的非随机信号部分(即图 1中的Others)，是由于反演策略强制海平面平衡方程闭合所引起.

图 1

Fig. 1

图 1 总体海平面及各个贡献成分的变化.蓝色线表示海平面变化的时间序列(季节性信号已经移除，数据作了垂直调整)，红色线表示由EEMD得到的非线性趋势变化. Fig. 1 The total sea level and different contributions. The blue line is the time series of sea level and different contributions (Seasonal signals have been removed. Vertical shift has been added for clarity), and red line indicates the non-linear trend derived by EEMD.

Rietbroek et al.(2016)的反演使用了GRACE数据(用于确定陆地水文和冰川的贡献)，因此提供的时间序列中某些月份的观测缺失.本文使用三次样条函数内插得到完整的时间序列，并利用最小二乘拟合移除季节性信号(周年+半周年).图 1给出了总体海平面和各个贡献成分的时间变化，可以看出，除陆地水文信号对海平面变化为负贡献外，其他各个成分均造成海平面的上升，而且在冰川消融方面，格陵兰的贡献要明显大于南极和山地冰川.值得注意的是，Rietbroek et al.(2016)反演得到的比容线性趋势为1.38 mm·a^-1(表 1)，要高于近年的研究，比如Chen et al.(2013)和Yi et al.(2015)的线性趋势约为0.8~0.9 mm·a^-1.

表 1 Table 1

表 1 全球平均海平面的瞬时速率(单位：mm·a-1) Table 1 The instantaneous rate of global mean sea level (unit: mm·a-1) Total Steric Greenland Antarctica Glaciers Hydrology Others MaxTime of Max 3.832014-06 1.602002-04 0.872012-11 0.342014-06 0.442002-04 0.032002-04 0.842014-06 MinTime of Min 2.542007-02 1.002014-06 0.512002-05 0.132002-04 0.352011-02 -0.252003-10 0.122008-11 2002-04 2.72 1.60 0.51 0.13 0.44 -0.24 0.16 2014-06 3.83 1.00 0.85 0.34 0.39 0.03 0.84 LSE* 2.74 1.38 0.73 0.26 0.38 -0.29 0.22 *LSE indicates the least-squares estimate of linear trend.

表 1 全球平均海平面的瞬时速率(单位：mm·a-1) Table 1 The instantaneous rate of global mean sea level (unit: mm·a-1)

以往的研究通常是利用最小二乘拟合获取某时间段内海平面变化的平均速率，这样的结果存在两方面的缺陷：(1)只能获取平均速率，而海平面是一种非稳态变化，其速率随着时间演变；(2)最小二乘估计受到海平面年际变化的影响(Cazenave et al., 2014).实际上，平均速率表征着海平面在某时间段内的线性变化，而实际上海平面变化具有非线性特征.相比于线性趋势，我们定义非线性趋势为一种具有缓慢变化的精炼趋势(Watson, 2016)，它不再受海平面年际变化的影响，即不再呈现高频的波动.利用非线性趋势，然后再进一步计算海平面变化的瞬时速率.

我们使用EEMD(Wu and Huang, 2009; Chen et al., 2017)来分离得到海平面变化的非线性趋势.EEMD是一种自适应数据处理方法，它不依赖于任何假设和先验信息，特别适用于非稳态变化的时间序列.EEMD的原理是将时间序列分解为一系列本征模态函数，每个本征模态函数表示某一频段的振荡变化.时间序列的高频变化和低频变化被分离到不同的本征模态函数中.本征模态函数的一个显著特点就是其振荡的振幅并非固定不变，而是会随着时间变化，即具有非稳态性.利用不同本征模态函数的算术和可以重新构建时间序列，进一步研究感兴趣的频段变化.通常情况下，最后一个本征模态函数呈现出单调变化，或者只有一个极值点.EEMD的详细计算步骤可参见Wu and Huang (2009).

使用EEMD分解总体海平面和比容海平面总共得到7个本征模态函数，图 2给出了后4个模态的变化，前面3个模态主要代表海平变化的高频振荡.可以看出第7个模态已经非常接近线性变化，代表海平面上升的总体趋势；第6个模态的波动幅度最小，有两个极值点，表示海平面变化的低频振荡；第5和第4模态振荡频率和振幅已经明显要高于第6模态，特别是第4模态在某些时间段内的振幅已经接近4 mm.在本文，我们定义第7和第6模态之和作为非线性趋势，用于计算海平面变化的瞬时速率.

图 2

Fig. 2

图 2 海平面变化的EEMD分解.使用EEMD共得到7个模态，(a)和(b)为总的海平面分解，(c)和(d)为比容海平面的分解；(a)和(c)为第7个模态；(b)和(d)为第4到第6模态. Fig. 2 The decomposition of sea level by EEMD. 7 modes are derived by EEMD, (a) and (b) are derived from total sea level, (c) and (d) are derived from steric sea level; (a) and (c) is the 7th mode, (b) and (d) are the 4th, 5th and 6th mode.

线性速率的本质是线性趋势的一阶导数，因此可以把非线性趋势的一阶导数作为瞬时速率.为得到在数学上连续的一阶导数，使用样条函数平滑(Watson, 2016)拟合非线性趋势，利用拟合曲线计算一阶导数，作为瞬时速率.由于样条函数是三次函数，它能保证一阶导数的在数学上连续性.在拟合时，为减小非线性趋势的信息损失，我们使拟合平滑曲线与非线性变化高度一致，即拟合曲线对非线性趋势的方差解释超过99.99%，即可以认为拟合曲线与非线性趋势重合.我们再次强调，使用样条平滑拟合非线性趋势是为获取在数学上连续的一阶导数，即瞬时速率.

显然，非线性趋势在瞬时速率的计算中起着关键作用.除EEMD外，也可以使用奇异普分析和卡尔曼滤波来分离或者计算非线性趋势，其中卡尔曼滤波也可以直接估计瞬时速率变化(Davis et al., 2012; Didova et al., 2016).奇异普分析也可以用于处理非稳态的时间序列，其效果与EEMD类似，也是将时间序列分解成为不同频率的振荡变化，其的详细计算过程可参见Shen et al.(2018).虽然奇异普分析也是一种自适应数据处理方法，但其结果依赖于一个重要参数：时间窗口，时间序列的矩阵排列取决于时间窗口的大小.通常情况下，时间窗口为总时间长度的1/4到1/2之间.我们发现，在使用奇异普分析分离海平面变化时间序列时，时间窗口大小对分离结果影响较大.卡尔曼滤波不仅可以计算非线性趋势变化，也可以直接计算瞬时速率及瞬时振幅(Davis et al., 2012).作为一种参数估计方法，卡尔曼滤波的计算结果取决于待估参数之间的协方差阵.此外，初始状态对结果也有显著影响(Didova et al., 2016).要获取最优的待估参数之间的统计信息，需要一些复杂的算法，其实现具有一定难度，且计算过程比较耗时.因此本文使用EEMD分离得到非线性趋势变化，与卡尔曼滤波相比，其计算过程要简单快速；与奇异普分析相比，无需考虑时间窗口的设置.

图 3给出了全球平均海平面及其各个贡献成分的瞬时速率.全球平均海平面的瞬时速率在2002—2007年经历了缓慢的下降，从2.72 mm·a^-1下降至2.54 mm·a^-1，在2007年2月达到最小(表 1).而在2009年，特别是2010年后，瞬时速率稳步上升，从2010年1月的2.71 mm·a^-1上升至2014年6月的3.83 mm·a^-1.与总体海平面相反，比容海平面瞬时速率在2002—2014年呈现持续下降的变化(瞬时速率在2009—2010年间有一个缓慢的抬升)，从2002年1月的1.6 mm·a^-1下降至2014年6月的1.0 mm·a^-1.这说明比容海平面对全球平均海平面上升的贡献在下降.

图 3

Fig. 3

图 3 海平面变化的瞬时速率(蓝色线).红色线表示95%的置信水平 Fig. 3 The instantaneous rate of sea level rise (blue line), the red line indicates the 95% confidence level

最近研究表明，格陵兰岛的冰川呈现出加速消融的模式(Velicogna et al., 2014; Wu and Heflin, 2015).瞬时速率的结果也证实了这一点，格陵兰岛冰川消融对海平面上升的贡献从2002年4月的0.51 mm·a^-1持续上升至2012年6月的0.87 mm·a^-1(图 3).于2012年11月达到最大后，瞬时速率有所下降，截止到2014年6月为0.84 mm·a^-1.南极的冰川消融同样对海平面上升有重要贡献.本文结果表明南极冰川消融导致的海平面上升的瞬时速率在持续上升，从2002年4月的0.13 mm·a^-1上升至2014年6月的0.34 mm·a^-1(图 3).速率的持续上升，表明南极的冰川消融在加速(Velicogna et al., 2014; Wu and Heflin, 2015).GRACE数据显示，南极冰川消融主要位于西南极，特别是阿蒙森海湾区域(Wu and Heflin, 2015).需要指出的是，GRACE在南极的结果仍存有较大差异，不同结果之间最大有100 Gt/yr的差异(Mu et al., 2017).

近年来，山地冰川的消融及其对海平面上升的贡献得到广泛关注(Jacob et al., 2012; Shepherd et al., 2012; Gardner et al., 2013; Schrama et al., 2014)，比如高亚洲冰川(Yi and Sun, 2014;朱传东等，2015)以及阿拉斯加冰川(Harig and Simons, 2016)等.陆地冰川消融导致的海平面瞬时速率在2002—2010年期间呈下降趋势，从0.43 mm·a^-1下降至0.35mm·a^-1；在2010年后，陆地冰川的消融又开始缓慢加速，其对海平面上升的瞬时速率上升至2014年的0.39 mm·a^-1(图 3).自从GRACE重力卫星发射成功，确定陆地水储量变化的研究成为可能(Wahr et al., 1998; Tapley et al., 2004;杨元德等，2009；Feng et al., 2013；冯伟等2017)，陆地水在调节全球水循环中具有重要作用.最近，GRACE揭示陆地水储量在缓和海平面的上升(Reager et al., 2016)，换言之，全球陆地水储量在2002—2014年呈现出增加的趋势，并且主要由气候变化所驱动.瞬时速率表明，陆地水储量对海平面上升的缓和作用在下降，从2002年4月的-0.24 mm·a^-1上升到了2012年6月的-0.02 mm·a^-1，2013年5月后，瞬时速率变为正值，即陆地水也开始促进海平面上升.

由于Rietbroek et al.(2016)的反演策略是强制海平面平衡方程闭合，因此得到了一组剩余的非随机信号.之所以称其为非随机信号时因为该信号呈现出明显的季节性和趋势性变化(图 1中的剩余信号Others已经移除了季节性变化).该非随机信号的线性趋势为0.22 mm·a^-1 (表 1)，但是我们发现其瞬时速率在经历了2002—2010年的缓慢变化后(图 3)，在2010年之后持续上升，从2010年0.16 mm·a^-1上升至2014年的0.82 mm·a^-1，剩余信号的瞬时速率与线性速率呈现出较大差异.

从数值看，瞬时速率与线性速率之间呈现出较大的差异(图 3和表 1).全球平均海平面的瞬时速率在2014年已经达到3.8 mm·a^-1，超过平均速率1 mm·a^-1.瞬时速率与线性速率之间的差异，与年际变化有很大关系，比如，卫星测高显示全球平均海平面在2003—2011年的线性速率约为2.4 mm·a^-1，但经过年际改正后，这一速率变为3.3 mm·a^-1 (Cazenave et al., 2014).使用EEMD分离得到非线性趋势变化，实际上就是消除了年际变化(特别是高频部分)的影响.本文瞬时速率的变化与滑动窗口拟合的结果有类似之处.基于5年滑动窗口拟合的结果，Yi et al.(2016)发现全球平均海平面的速率在2010年后快速上升(即超过了1993—2014年的平均速率)，截止到2012年的时候已经接近4 mm·a^-1.

Rietbroek et al.(2016)使用的比容数据来自模型，并非直接的观测，其线性速率与Argo的结果有一定差异.我们分析了Argo比容海平面(数据由日本海洋地球科学技术中心提供，处理过程参见陈威等(2017))的瞬时速率，发现Argo的瞬时速率也呈现出下降的变化，从2002年的2.5 mm·a^-1下降至2014年的0.5 mm·a^-1左右.我们指出，尽管Argo是基于直接观测的结果，但其在时间和空间的采样仍有一定的局限性，比如Argo只能探测0~2000 m的变化，空间覆盖范围在南北纬65°之间(钟玉龙等，2016)，这些可能是导致Argo与模型之间差异的部分原因.由于本文讨论的重点是海平面变化的瞬时速率，对不同数据之间的差异原因不作深入讨论.

剩余信号与总体海平面的瞬时速率变化十分类似(图 3)，特别是2010年后表现为持续上升的趋势.尽管在理论上，卫星测高观测的总体海平面等于海水质量变化(或者全部的冰川消融和陆地水储量变化)与海水比容之和，但截至目前为止，在全球平均尺度上，海平面平衡方程尚不能被完全闭合(Chen et al., 2013).可能的原因有，比如比容数据和GRACE均有观测误差和数据处理的误差，并且卫星测高目前还不能覆盖北冰洋等.考虑到Rietbroek et al.(2016)对比容海平面的估计要高于最近的结果，比如Chen et al.(2013)和Yi et al.(2015)，并且Rietbroek et al.(2016)也指出，其反演结果可能低估了冰川消融.因此我们可以推测，剩余信号在2010—2014年的变化主要由海水质量引起.

我们把格陵兰岛、南极地区、山地冰川和陆地水储量的贡献归为海水质量变化，这样我们可以从海水质量和海水比容的瞬时速率计算其贡献比率(图 4).从图 4中可以看出，海水质量的贡献在2009年超过了比容的贡献.比容的贡献在2002—2014年期间呈现出下降趋势，由57%下降到30%左右，如果我们把2010年后的剩余信号完全归结为海水质量变化，那么海水质量的贡献在2014年就达到了70%，这与Chen et al.(2017)的结果相一致.

图 4

Fig. 4

图 4 海平面变化贡献成分的比率 Fig. 4 The ratio of different contributions to sea level rise

本文提出了一种新的海平面瞬时速率计算的方法.瞬时速率主要取决于海平面的非线性趋势变化.使用EEMD分离海平面时间序列得到具有不同振荡变化的本征模态函数，我们将最后两个模态之和作为非线性趋势.使用三次样条平滑拟合非线性趋势得到瞬时速率变化.

线性速率只能给出海平面上升或下降的总体趋势，而瞬时速率提供了海平面上升随时间演化的特征，提高了速率在时间变化上的分辨率.瞬时速率增进了我们对海平面变化的认识，比如在2010年后，全球平均海平面的瞬时速率在持续上升，截止到2014年为3.8 mm·a^-1，表明在此期间全球平均海平面在加速上升.同时，分析海平面上升各个贡献成分的瞬时速率，使我们认识到比容海平面的贡献在逐渐下降，而海水质量的贡献在增加.线性速率只给出海平面各个成分的平均贡献，而瞬时速率则揭示了海平面的动态变化，即各个成分的贡献随着时间在不断变化.

需要指出的是，瞬时速率只能提供海平面变化的规律，而这些规律变化的机制仍需要被解释，比如剩余信号的归因，其究竟是由质量变化引起，还是由比容变化导致，依赖于更高精度的观测.

本文分析的是全球平均海平面的变化，在区域尺度上，海平面变化与全球平均会有较大差异，将瞬时速率应用至区域海平面变化，特别是空间模态的变化，将有助于增进对海平面变化的认识.本文研究的时间段只有12年，但瞬时速率可以应用于更长的时间段，比如长时间序列的验潮站和重构的20世纪全球平均海平面变化.此外，本文瞬时速率的计算方法也可以用于其他地球物理观测的时间序列，比如地表温度变化等，瞬时速率有助于提高我们对地球物理变量在全球变暖背景下随时间演化特征的理解.