在人类温室气体持续排放等重要外强迫的作用下, 地球气候系统存在长期变化趋势。准确评估地球气候系统中海洋和大气温度、海平面等各要素的长期变化趋势, 是理解地球气候变化和预测未来气候面临的首要关键问题。但是, 目前对地球气候系统各要素长期趋势进行准确估计实际上非常困难, 还存在很大的挑战。一方面, 人类对地球系统的观测资料非常有限。以海洋温度观测为例, 直到20世纪中叶才开始通过抛弃式温度剖面仪(expendable bathy thermograph, XBT)等手段对海洋上层温度进行大面积观测, 20世纪70年代才通过卫星对海表面温度(sea surface temperature, SST)进行大面持续观测, 21世纪初才通过大量布放Argo剖面浮标获取海洋2 000米水深内的温度和盐度数据。这些观测的持续时间都不长, 基本上在几十年左右。早期的观测手段误差很大, 而观测平台之间的转变还引入了系统性误差。历史观测数据中, 在Argo计划实施之前的早期观测数据存在采样空间分布不均和极其稀疏的特点, 而在Argo计划实施之后数据采样状况迅速转好, 这导致在Argo逐渐建立起来的时间段(2001~2003年)内产生虚假的跳跃和系统性误差^[1]。观测资料时间长度和采样率方面的天然缺陷, 使得长期趋势估计中的不确定性非常大^[2]。另一方面, 地球气候系统是一个复杂的耦合系统, 在自然的动力过程(尤其是海洋动力过程)的调控下, 地球气候系统存在强烈的自然振荡。气候系统的自然振荡导致在特定观测期间的趋势估计非常容易受到其影响, 尤其是自然振荡强烈的物理量^[3-4]。

以印度尼西亚贯穿流(简称印尼贯穿流)的多年代际趋势估计为例。印尼贯穿流的多年代际趋势变化通过强烈的热盐输送, 调控太平洋海盆和印度洋海盆的物质能量平衡^[5-9], 对区域和全球气候系统产生重要影响, 大量研究表明印尼贯穿流在厄尔尼诺-南方涛动(El Niño–Southern Oscillation, ENSO)、全球变暖停滞(global warming hiatus)等具有全球性影响的气候过程中扮演了非常关键的角色^[10-13]。准确估计印尼贯穿流多年代际趋势是理解和预测地球气候系统的关键问题之一。

图 1展示了印尼贯穿流流量异常的时间序列, 包括使用EN4^[14]、Ensemble Coupled Data Assimilation(ECDA)^[15]和European Centre for Medium-Range Weather Forecasts Ocean Analysis/Reanalysis System 3(ECMWF-ORA S3)^[16]三种再分析数据反演的印尼贯穿流流量的合成序列和使用格点化Argo数据^[17]反演的印尼贯穿流流量异常时间序列, 数据来自Hu和Sprintall^[7-8]的相关研究。集合反演的印尼贯穿流流量异常时间序列与Argo的结果表现出非常一致的变化, 相互间的整体差异可能是由于二者参考的平均印尼贯穿流流量时间不一致所导致。Hu和Sprintall^[8]的计算结果显示, 1985年之前计算的印尼贯穿流流量存在很大的误差, 因此本文仅选取1985年之后的时间序列进行分析。1985~2010年, 印尼贯穿流流量的粗估线性趋势约为3.8 Sv/20yr, 去除线性趋势后印尼贯穿流流量的变化(标准差)约为4.1 Sv。在该时间序列中年代际变化为主要部分, 印尼贯穿流流量的自然振荡大于粗估趋势20年累计的变化, 自然振荡不同位相的选择对粗估趋势影响很大。可见, 印尼贯穿流的年代际自然振荡对于趋势估计具有非常重要的影响。

由于历史的原因, 观测资料的时间长度缺陷无法弥补, 由此引起的趋势估计不确定性也难以克服。但是, 由于地球气候系统中的自然振荡有许多是周期或准周期振荡, 准确估计气候系统中物理量的长期趋势的一个可行的方法是尽量去除其自然周期/准周期振荡的影响。由于目前海洋和气候的观测资料多数集中在几十到100多年的长度, 去除自然周期/准周期振荡(尤其是年代际振荡)的影响对于研究印尼贯穿流以及地球气候系统中其他重要物理过程的多年代际趋势尤为重要。

为了研究自然周期振荡对趋势估计的影响, 本文从基于一个叠加周期振荡和趋势变化的理想数学模型^[18], 探讨趋势估计中去除周期振荡影响的普适性方法, 并将其具体应用在印尼贯穿流多年代际趋势变化的研究中。

1 滑动趋势法的原理和效果

假设一个由自然振荡和长期趋势线性叠加的逐月时间序列ψ(t), 定义其数学模型为:

$ \psi(t) = \sin (f t)+a t, $

(1)

其中系数设置为: $f = \frac{1}{60} \pi$, α = 0.002。则ψ(t)的自然振荡周期为$T_{0} = \frac{2 \pi}{f} = 120$月, 理论的趋势为0.002每月。在实际计算中, 趋势估算为ψ(t)的一阶导数, 即

$ \frac{\mathrm{d} \psi(t)}{\mathrm{d} t} = f \cos (f t)+\alpha, $

(2)

其物理意义为时间序列的瞬时趋势。很显然, 这个瞬时趋势与理论真实趋势之间的差为fcos(ft), 而这项是由于ψ(t)中的周期性振荡引起的。在研究长期趋势时, 需要把fcos(ft)的影响去除掉。对公式(2)进行积分, 可得

$ \int_{0}^{T} \frac{\mathrm{d} \psi(t)}{\mathrm{d} t} \mathrm{d} t = f \int_{0}^{T} \cos (f t) \mathrm{d} t+\alpha T, $

(3)

从而有

$ \int_{0}^{T} \frac{\mathrm{d} \psi(t)}{\mathrm{d} t} \mathrm{d} t = \sin (f T)+\alpha T, $

(4)

当T为自然振荡周期T₀的n倍(n为整数)时, sin(fT) = sin(2nπ) = 0, 则

$ \int_{0}^{T} \frac{\mathrm{d} \psi(t)}{\mathrm{d} t} \mathrm{d} t = \alpha T, $

(5)

因此, 可以看出当样本长度为自然振荡周期的整数倍时, 自然周期性振荡对趋势估计的影响最小。

图 2中蓝色实线给出的是ψ(t)的3个周期长度, 亦即360个月的时间序列, 蓝色虚线给出的是直接对ψ(t)|_{t = 0~360}求线性趋势得到的结果。彩色圆点连接的彩色线为不同整周期时间段的线性趋势。由图 2可见, 各趋势之间差别很大, 趋势值随样本截取起始位相的变化而变化。之所以出现整周期样本趋势间的差异, 是因为各整周期样本的截取位相不同, 而在公式(3)~(5)的积分中, 暗含了一个重要的假设是起始位相点恰好是在振荡的原点, 即ψ(0) = 0。在实际研究过程中, 由于样本量稀缺, 难以保证样本起始点从零位相开始。

为了尽量避免由于起始位相点引起的趋势估计误差, 可以设定一个滑动的窗口, 在窗口内求线性趋势, 并将滑动平均的趋势值作为估算的趋势值。即:

$ D_{\psi(t)} = \frac{\sum\limits_{\tau = 1}^{\tau = N} D_{\psi\left(\tau, \tau+n T_{0}-1\right)}}{N}, $

(6)

式中, ${{D}_{\psi (t)}}$为估算的ψ(t)的趋势, ${{D}_{\psi \left(\tau, \tau +n{{T}_{0}} \right)}}$为滑动窗口内的趋势, 即为整数倍周期以内、时间段$\tau \sim\left(\tau+n T_{0}-1\right)$内的趋势, n为窗口内周期数, τ为样本时间点, N为去除整数倍周期后剩余的样本量时间长度。为便于区分, 我们将公式(6)估算到的趋势称为“滑动趋势”(running trend), 而直接估算且不考虑自然振荡的趋势称为“粗估趋势”(rough trend)。本文所指的粗估趋势, 是传统的线性趋势估计, 使用最小二乘拟合方法所得, 但不考虑自然振荡的影响。

依据该方法, 我们现在计算ψ(t)的滑动趋势, 并与粗估趋势和理论趋势做比较。图 3展示的是计算结果。由于周期振荡的影响, 粗估趋势方法得到的结果非常小, 远远小于真实的趋势值, 即远小于理论趋势值, 这表明粗估趋势是不正确的。而反观滑动趋势法的计算结果, 其趋势线几乎与理论真实趋势线重叠, 表明滑动趋势法能有效去除自然振荡的影响, 估算出接近真实的趋势。

从滑动趋势法的定义可以看出, 滑动趋势本身是一组线性趋势的组合。每一个线性趋势本身的置信水平是可以基于t-检验等方法进行评估的。对滑动趋势来说, 其对应的置信水平其实是一组置信水平的集合, 假设这组线性趋势各自独立平等, 则滑动趋势的置信水平可以近似为这组线性趋势置信水平的平均值。本文着重刻画趋势估计中的物理过程, 因此对该问题不再进行更多讨论。

2 样本量和周期估计的影响

在滑动趋势方法的使用中, 有两个因素对趋势估计的准确度有非常重要的影响, 即样本量的大小和样本自然振荡周期估计的准确性。

我们用理论模型ψ(t)进行敏感性测试。分别选取不同样本长度的ψ(t)时间序列进行研究, 同时选用不同周期T₀, 分别计算相应的滑动趋势。如图 4所示, 当样本量从3个周期长度增加到4个和5个周期长度时, 滑动趋势的整体误差大幅度减少, 粗估趋势误差均相对较大。由此可见, 样本量越大, 得到的趋势估计越能接近真实值, 但这里所谓的样本量大是指相对于自然振荡周期而言。在一个给定的时间序列中, 其自然存在的高频变化因为样本量较大, 整体来说对于趋势估计的影响较小(图 4c)。对于趋势估计影响最大的信号, 就是周期接近于样本总长度的自然振荡(或者说周期不远小于样本总长度), 比如周期为样本总长度约三分之一的自然振荡(图 4a)。

另一个影响滑动趋势估计准确性的因素是自然振荡周期估计的准确性。图 4横轴为滑动趋势计算过程中使用的不同周期。不论样本量多大, 当给定的样本周期等于理论真实周期T₀时, 估算的滑动趋势等于真实值。但是, 当给定的样本周期偏离自然振荡的真实周期时, 滑动趋势法给出的趋势也会逐渐偏移真实的趋势值。

由此可见, 在使用滑动趋势法时, 一方面应尽量选择足够长的样本以保证整体的准确性, 重点需要关注周期接近总样本长度(如一半或三分之一长度)的自然振荡, 在固定时间段内高频信号的样本数相对较多, 对于趋势估计的影响较小。另一方面应尽量准确估算样本的周期, 以减少估算的误差。当前研究中, 对于样本周期的估计主要采用谱分析的方法, 本文采用的是功率谱分析。

3 在印尼贯穿流多年代际趋势研究中的应用

滑动趋势法具有相当的普适性, 可以应用在很多领域。本文仅将其运用于作者非常关注的印尼贯穿流流量的多年代际趋势变化研究中。如前文所述, 人类对海洋(尤其是深层海洋)的观测时间不过数十年而已, 但是现有研究已表明, 海洋气候系统存在显著的年代际变化。对于时间长度为几十年的海洋观测样本来说, 高频信号(如季节内到季节乃至年际信号)对趋势估计的影响较小(图 4c), 但是年代际信号却能在趋势估计中引入巨大的误差, 因此不得不考虑采用滑动趋势法来去除自然振荡的影响。

对印尼贯穿流的长期连续观测非常缺乏, 仅在望加锡海峡有长达约14年的潜标观测数据^[19]。但由于受到太平洋年代振荡(Pacific decadal oscillation, PDO)等影响^[8], 印尼贯穿流自身存在强烈的年代际变化, 14a的样本长度不过是自然振荡的一个周期左右, 由此估算得到的趋势可想而知很大一部分是自然振荡信号, 而非真实的趋势。为此, 我们采用多种再分析温盐数据反演的印尼贯穿流流量的集合平均值作为印尼贯穿流时间序列^[7-8]。首先对印尼贯穿流流量时间序列进行功率谱分析(图 5a), 可以看到该时间序列在周期为12a处存在非常显著的谱峰, 可以认为印尼贯穿流的主要年代际振荡周期为12a。于是, 进而将T₀设为144个月进行滑动趋势计算。结果发现, 1985~2010年, 印尼贯穿流的趋势变化约为7.3 Sv/20yr, 而粗估趋势约为3.8 Sv/20yr, 因此实际的趋势变化远大于粗估趋势值, 是后者的大约2倍。

以往研究报道中常见对物理量线性趋势开展空间模态的研究^[20], 但是, 由于粗估趋势存在很大的误差, 以往研究中常见的线性趋势空间分布也存在很大误差, 粗估趋势的所谓空间分布实际上是周期振荡混合了实际趋势的结果, 在用于研究物理量长期变化的空间模态时可能会产生严重的虚假信号。图 6a是用滑动趋势法计算的海表面高度异常(sea level anomaly, SLA)在1993~2010年的趋势, 由于样本长度为18a, 功率谱分析表明主要振荡信号为周期为4a的年际信号, 因此滑动窗口设为4a。热带印太海洋表现出显著的增高趋势, 但是西太平洋SLA的增加趋势明显大于东印度洋的趋势, 这导致西太平洋向印度洋的压力梯度增大, 印尼贯穿流增强。

Susanto和Song^[21]指出, 位于西太平洋162°E, 11°N处SLA和位于印度洋80°E, 0°N处SLA的差可以作为印尼贯穿流流量的指示因子。如图 6a所示, 西太平洋162°E, 11°N处的趋势约为9 cm/10yr, 而东印度洋80°E, 0°N的趋势约为2 cm/10yr, 这说明西太平洋向印度洋的压力梯度差有增大趋势, 约为7 cm/10yr, 这对应着印尼贯穿流的增强趋势, 与现有的观测结果吻合。

但需要说明的是, 滑动趋势和粗估趋势之差形成的误差空间模态值得关注。图 6b给出了AVISO SLA在1993~2010年间的滑动趋势减去粗估趋势的差。可以很清晰地看到, 二者的差存在显著的空间模态, 在以印尼海域为中心的热带印太海洋为正值, 即粗估趋势中高估了该区域SLA的升高趋势, 但是在西南印度洋和热带中太平洋区域却表现为负值, 即粗估趋势中低估了这些海域SLA的升高速度。由于这种空间模态的存在, 导致粗估趋势中西太平洋和东印度洋之间SLA差偏大, 进而对印尼贯穿流增长趋势的估计也偏大。

这种差异的空间模态的形成可能是粗估趋势中包含强烈的年际变化信号导致的。在1993~2010年间, 发生了6次El Niño事件, 包括1次超强El Niño事件(1997~1998年), 3次中等强度El Niño事件(1994~1995年、2002~2003年、2009~2010年), 和2次弱El Niño事件(2004~2005年和2006~2007年)。在此期间发生了7次La Niña事件, 包括3次强La Niña事件(1998~1999年、1999~2000年、2007~2008年), 1次中等强度La Niña事件(1995~1996年), 和3次弱La Niña事件(2000~2002年、2005~2006年和2008~2009年)。可见在1993~2010年发生了更多的La Niña事件, 而且1993~2010年平均的Niño3.4指数为-0.04℃。这些数据表明在此期间主要的年际变化表现为La Niña信号, 粗估估计中受到了La Niña的影响, 相对于滑动趋势表现出La Niña模态, 印尼贯穿流表现出相对偏强的特点, 这与图 6b的特点一致(即海表面高度在西太平洋为正异常而中东太平洋为负异常的特点), 同时也说明去除自然振荡对于趋势估计存在重要意义。

4 结论

如何在物理量的趋势估计中去除掉自然周期性振荡的影响是气候变化研究中的关键科学问题之一。本文针对该问题, 基于理想化的数学模型和对印尼贯穿流为代表的实际研究工作, 提出了滑动趋势法来去除或减少自然周期振荡对趋势的影响, 从方法的基本原理和效果进行了阐述, 通过理论模型对该方法进行了检验, 表明滑动趋势法可以得到接近真实趋势的估计值, 值得推广使用。

本文利用滑动趋势法估算了印尼贯穿流1985~ 2010年的多年代际趋势, 发现该方法得到的趋势比粗估趋势大了近一倍, 这对理解印尼贯穿流的真实变化趋势具有非常重要的意义。但需要说明的是, 该方法具有普适性, 并不局限于印尼贯穿流研究, 也不仅仅在地球气候变化研究中使用, 而是可以应用在任何变量的趋势研究中。