厄尔尼诺是发生在热带太平洋的一种年际异常强信号, 它是大尺度海-气相互作用的产物^[1~3]; 而赤道东太平洋海表温度的异常发展, 也会对大气环流产生影响, 进而牵制包括中国在内的全球许多国家和地区的气候变化, 造成旱涝异常和高低温灾害^[4]; 因此, 提前预知厄尔尼诺的发生发展对于区域气候预测、防灾减灾都具有重要意义。然而当前由于对厄尔尼诺物理机制认知的局限以及观测资料的不足, 厄尔尼诺的预测仍处于试验阶段^[5]: 20世纪80年代以来, 世界各国陆续开展了厄尔尼诺的预报研究, 并在过去的几十年发展出近20种相应的预报模式; 我国于20世纪80年代中期建立了厄尔尼诺监测小组, 20世纪90年代建立了厄尔尼诺监测和预测业务系统, 到今天已发展出了两个厄尔尼诺预报系统:中国科学院大气物理研究所的集合预报版本(IAP Leefs_CDA)和中国科学院海洋研究所的确定性预报版本(IOCAS ICM), 均提供实时的预报结果。目前国际上的近20种厄尔尼诺预报模式, 大部分可以提前6个月给出较为可信的预报, 但由于模式的系统误差, 使得各模式的预报结果间存在较大差距, 同时数值模式对于极端事件的预报效果较差, 难以满足预报需求^[6~8]。近年来, 一些学者开始尝试使用统计方法开展厄尔尼诺预报试验, 并取得了一定成果。韦成宁^[9]通过对历年厄尔尼诺与与太阳黑子、大气环流的关系分析, 提出了厄尔尼诺的概率统计预报方法; 魏松林^[10]提出了一个表征厄尔尼诺/拉尼娜事件的指数——ENI指数, 并在逐月ENI指数谱分析及主周期外推预报的基础上提出了三周期叠加预报方法, 提前6个月成功预报出了1997年的超强厄尔尼诺事件; 冯利华等^[11]将最早在天文学领域提出的“可公度性”概念应用于厄尔尼诺的预报研究中, 成功地预测了1997年及2001年的厄尔尼诺事件。尽管学者们在厄尔尼诺预测方面做了一些尝试, 但这还远远不够, 受预报时效、预报准确度的影响, 厄尔尼诺的预测仍然是气候预报中的一个薄弱环节。因此, 本文尝试利用在极值预测方面有较好表现并能进行多步预测的均生函数模型对厄尔尼诺进行预测试验, 并提出了两种改进方案, 以期为厄尔尼诺的预报提供新的思路与方法。

1 资料与方法 1.1 资料

本文所用的NINO3指数来源于美国国家海洋和大气管理局(National Oceanic and Atmospheric Administration, NOAA)气候预测中心, 时间长度为1970年1月~2017年12月; 同时段月平均海温来源于COBE-SST全球海温数据集, 分辨率1°×1°; 月平均高度场、风场来源于NCEP/NCAR再分析资料, 水平分辨率2.5°×2.5°; 向下太阳辐射通量同样为再分析资料, 由逐日数据计算为月平均值, 资料采用高斯格点分布。在应用过程中, 各类数据都进行了距平或标准化处理, 以滤除季节波动。

1.2 方法

均生函数^[12~14]是本文三种统计预测方法的理论基础, 它是由气候序列按一定时间间隔计算均值而派生出来的。对于长度为n的时间序列:

$x(t) = \{ x(1), x(2)x(1), \cdots , x(n)\} .$

(1)

定义其均值生成函数:

${\bar x_l}(i) = \frac{1}{{{n_l}}}\sum\nolimits_{j = 0}^{{n_l} - 1} {x(i + jl)\;\;\;(i = 1, \cdots , l;1 \le l \le m)} , $

(2)

其中, ${n_l} = {\rm{INT}}(n/l), m = {\rm{INT}}(n/2)$或INT(n/3), INT表示取整数。

将均生函数定义域延拓到整个数轴上, 即按照公式(3)作周期性延拓:

${f_l}(t) = {\bar x_l}\left[ {t - l \times {\rm{INT}}\left( {\frac{{t - 1}}{l}} \right)} \right](t = 1, 2, \cdots , n;l = 1, \cdots , m).$

(3)

由此构造均生函数周期延拓序列。简单的均生函数模型就是在这组延拓序列的基础上, 建立与原序列之间的回归, 构建预报方程。

本文在改进均生函数预报模型的过程中, 还用到了小波分析、集合经验模态分解(Ensemble Empirical Mode Decomposition, EEMD)以及奇异值分解(Singular Value Decomposition, SVD)等方法, 详细介绍可参考文献[15]~文献[19], 这里不再赘述。

1.3 建模步骤 1.3.1 方案1——普通均生函数模型

设有长度为n的气候序列, 如公式(1)所示。

(a) 对原序列进行差分计算, 以便拟合序列中的高频部分:

$ \Delta x\left( t \right) = x\left( {t + 1} \right)-x\left( t \right)(t = 1, 2, \ldots , n-1), $

(4)

$ {\Delta ^2}x\left( t \right) = {\rm{ }}\Delta x\left( {t + 1} \right)-\Delta x\left( t \right)(t = 1, 2, \ldots , n-1). $

(5)

得一阶差分序列x⁽¹⁾(t)和二阶差分序列x⁽²⁾(t)。

(b) 对气候序列及其差分序列按照公式(2)分别构造均生函数${\bar x_l}(t)$、$\bar x_l^{(1)}(t)$和$\bar x_l^{(2)}(t)$, 并对三组均生函数利用公式(3)进行延拓, 依次记为$f_l^{(0)}(t)$, $f_l^{(1)}(t)$和$f_l^{(2)}(t)$。

(c) 为了拟合时间序列中的增减趋势, 进一步建立累加延拓序列:

$ f_l^{(3)}(t) = x(1) + \sum\nolimits_{i = 1}^{t - 1} {f_l^{(1)}(i + 1)\;\;\;(t = 2, \cdots , n;l = 1, 2, \cdots , m)} . $

(6)

这样共获得4m个均生函数延拓序列: $f_l^{(0)}(t)$, $f_l^{(1)}(t)$, $f_l^{(2)}(t)$和$f_l^{(3)}(t)$作为自变量供预报方程筛选。

(d) 建立每一个均生函数延拓序列与预报量的一元回归, 计算双评分准则SCS(couple score criterion)值, 将满足SCS> $X_\alpha ^2$的序列粗选为预报因子。设入选了p个延拓序列。

(e) 对所有粗选出的预报因子, 计算2^p个回归子集, 并依据SCS准则选出其中最优的预报方程。若做q步预报, 则对因子序列做q步外延即可。

均生函数模型自20世纪90年代提出以来, 就在气象、水文、航天等领域得到了广泛应用; 但模型中所用的因变量都是由其自身派生构造出来的, 具有一定的局限性。因此, 本文在均生函数模型的基础上, 还应用了两类改进模型, 尝试提高厄尔尼诺的预报技巧。

1.3.2 方案2——基于EEMD的均生函数模型

均生函数模型通过时间序列中隐藏的各种周期和演变趋势进行构造^[20]; 因而本文在提取气候变量周期的基础上, 进一步提出了基于EEMD的均生函数模型。EEMD是分析非线性、非平稳信号序列的新方法, 可以将复杂的信号分解为有限个本征模函数IMF和一个趋势项RES, 所分解出的各个IMF分量均包含了原信号不同时间尺度的局部特征。基于EEMD的预报模型相当于普通模型的备选预报因子中, 加入了预报量EEMD分解结果的均生函数序列; EEMD可以将气候序列平稳化处理, 各IMF周期单一且显著, 具有更好的可预报性。具体建模步骤如下:

(a) 对气候序列两端进行延拓, 这是针对端点效应的处理。EEMD的分解过程中信号两端会出现逐渐向内部扩散的失真现象。通常的处理方法是分解后删去部分端点处序列, 但为了不失去预报信息, 本文首先利用镜像延拓方法扩展序列, 实践表明, 镜像延拓能大幅削减端点失真所造成的影响^[20]。

(b) 对延拓后的序列进行EEMD分解, 得k个IMF; 分解后删去延拓部分。

(c) 取原序列及k个IMF, 按照公式(4)、公式(5)分别计算它们的一、二阶差分序列, 再根据公式(2)、公式(3)求所有均生函数周期性延拓序列, 最后根据公式(6)建立k+1个累加延拓序列; 最终得4m×(k+1)个均生函数周期性延拓序列作为预报方程的备选因子。

(d) 由于预报因子过多, 采用最优子集回归建立预报方程的计算量过大, 因而采用逐步回归的方法筛选预报因子、建立预报模型。对入选方程的因子做q步外延即可得到q步预报值。

1.3.3 方案3——多变量均生函数模型

受普通均生函数模型变量单一的局限, 本文在其基础上, 综合预报量及其影响因子序列, 构建多变量均生函数模型。对模型加入因子变量有两方面的意义^[21]:一方面作为备选的自变量数量增加, 可以通过筛选挑选出最佳因子进入预报方程; 另一方面, 气象要素的变化是多种因素共同作用的结果, 模型中加入多种相关的气象因子, 具有更好的物理基础。

设预报量为y(t), 并有k个因子序列: x₁(t), x₂(t), …, x_k(t), t=1, 2, …, n。n为样本长度。

(a) 为消除变量间的差异, 对预报量和所有因子变量进行标准化处理:

其中, $\bar y, {\bar x_i}, {s_y}, {s_{{x_i}}}$分别为y(t)、x_i(t)的均值和方差。

(b) 对标准化后的k+1个变量序列$y', {x'_1}, {x'_2}, \cdots , {x'_k}$按照公式(2)~公式(6)分别计算它们的均生函数延拓序列、一阶和二阶差分序列的均生函数延拓序列, 以及累加延拓序列。

(c) 通过逐步回归筛选预报因子并建立预报方程, 对入选方程的因子做q步外延即可得到q步预报值。

2 结果与分析

2014年底, 海表温度异常增暖率先出现于赤道中太平洋, 伴随赤道西风的增强和暖水东传, 2015年初, NINO3区海温达到厄尔尼诺标准, 同年冬季增温2~3℃以上, 为本世纪以来最强的一次厄尔尼诺。对于此次事件, 国际上近20种厄尔尼诺模式提前6个月, 只有一个成功预测出本次事件强度将超过2℃; 直到2015年6月, 大部分模式才预报出年底可能出现厄尔尼诺^[7]。因此本文首先利用三种模型对2015年的超强厄尔尼诺进行预报, 1980~2014年为建模时段, 逐年外推三步, 2015~2017年用于预报检验。在此基础上, 建立滑动的预报模型, 对2001~2017年的NINO3指数进行预测, 并给出整体预报评价; 滑动的样本同样取预报年份之前的35 a, 这是因为:数理统计中通常要求样本序列n不小于30, 但均生函数模型并不作严格限制, 理论上有2个以上数据求平均即可, 但时间序列所能构造生成的均生函数个数与样本量n相关。n越大, 不仅可以为模型提供更多备选因子, 而且可以降低均值序列的随机性; 但n增加到一定程度, 预报技巧的提高将变得缓慢, 还会导致计算量增大。因此, 通过反复试验, 本文将样本量n确定为35 a。

2.1 2015年超强厄尔尼诺的个例预报试验 2.1.1 普通均生函数模型在2015年强厄尔尼诺预报中的应用

表 1为普通均生函数模型的预报结果与实测值, 可以看出, 预报效果整体表现一般。外推一步, 前三个月预报与实测值都处在正常值范围内, 4月NINO3指数开始达到0.5℃以上, 而预报值直到8月才达到厄尔尼诺水平; 此外, 二者在数值上也有较大差距, 预报的厄尔尼诺强度偏弱。外推两步, 距平符号与实况基本一致, 但预报的厄尔尼诺结束时间偏早, 仅在8~11月有较好的预报效果。外推三步, 各月距平符号几乎都不一致, 已经失去了预报意义。

2.1.2 基于EEMD的均生函数模型在2015年强厄尔尼诺预报中的应用

厄尔尼诺作为强气候信号, 每隔2~7 a出现一次, 但并不遵循严格的周期^[22]。近年来厄尔尼诺的发生似乎有愈加频繁的趋势, 其中2002~2010年接连发生四次厄尔尼诺。图 1给出了近38 a年均NINO3指数的小波分析, 可以看出:厄尔尼诺主要表现出两个变化周期:一是2~7 a的年际振荡, 二是9~14 a的年代际振荡; 其中年际振荡更为明显, 为其第一主周期; 值得注意的是, 20世纪90年代后, 年际振荡随时间减弱, 且中心逐渐上移, 但在2~3 a尺度上, 有新的闭合中心发展; 另外, 1982、1998和2015年的强厄尔尼诺事件, 在年际和年代际振荡上表现为一致的正位相(1998年处在更长的年代际振荡正位相上); 不同周期的振荡同位相叠加, 形成强劲的厄尔尼诺, 这可能在某种程度上反映了强厄尔尼诺的特点。

进一步对NINO3指数进行EEMD分解, 得到的四个IMF对应的周期依次为: 3~6 a、6~14 a、15~26 a及30 a以上的振荡, 而趋势项表明NINO3区海温经历了先降后增的变化, 近年来海温正处于持续上升阶段, 转折点在20世纪90年代中期。由于前两模态与原序列的相关性较好, 且方差贡献占比较大, 是NINO3指数变化的主要贡献者, 因而图 2主要给出了EEMD分解的前两个模态: IMF1跟原海温序列变化趋势极其接近, 2000年以后, IMF1振幅减小, 振荡周期缩短, 这正体现了21世纪前十年频繁出现的厄尔尼诺。IMF2反映了更长周期的演变特征, 同样在2000年以后振荡频率加快。此外, 两个不同尺度振荡同正位相叠加, 是出现强厄尔尼诺的主要特征。

上述分析表明, 厄尔尼诺具有显著的多时间尺度特征。为建立预报模型, 对建模时段各月NINO3指数同样进行EEMD分解, 各月均得四个IMF和一个RES, 其与年均值分解结果相似, 仅前两模态的相关系数较高; 不同的是分月的IMF1主要表现2~5 a振荡, IMF2反映5~8 a振荡, IMF3体现10~22 a年代际振荡; 各月模态周期基本一致, 少数月份存在一定差异。

在EEMD分解基础上, 通过逐步回归建立分月预报模型; 从各月回归方程入选的预报因子可以看出, 预报量及其EEMD分解前两模态的均生函数延拓序列所占的比例较高, 这也例证了反映预报量“优势周期”的因子序列更容易入选进方程。

表 2给出了基于EEMD的均生函数模型对2015~ 2017年NINO3指数的预报值及实测值, 预报效果较普通均生函数模型有了明显改进。外推一步, 仅4~6月预报效果相对较差; 7月, 模型预报出NINO3指数开始达到厄尔尼诺水平, 期间预报值与实况值较为接近。外推两步, 模型预报出的厄尔尼诺结束时间与实况仅差两个月; 厄尔尼诺结束后紧随而来的海温负距平, 在预测值上也有很好的体现。外推三步, 预报同样失去意义。

2.1.3 多变量均生函数模型在2015年强厄尔尼诺预报中的应用

厄尔尼诺是复杂的海-气相互作用的结果, 为了了解厄尔尼诺与其他因子的相关关系, 本文运用SVD分析技术, 主要从大气环流和太阳辐射的角度, 为厄尔尼诺的预测提供参考。将赤道太平洋海温设为左场, 而右场依次取1 000 hPa高度场、500 hPa高度场、850 hPa纬向风场、850 hPa经向风场、500 hPa垂直速度场以及到达地面的向下太阳辐射通量, 第一模态的方差贡献至少在70%以上(表 3), 远大于其他模态; 同时, 各因子第一模态的左异类相关图表示的海温分布形式相似, 大值中心均位于赤道中、东太平洋; 同时左场时间系数与NINO3指数的相关非常高, 至少在0.94以上, 并全部都通过了99%的显著性检验, 反映了厄尔尼诺的海温分布特征。因此, 在以下分析中仅讨论每种分解结果的第一特征向量。

第一模态的右异类相关图, 即1 000 hPa高度场(图 3a)上, 与厄尔尼诺相关的分布形式表现为明显的东西反相特征, 正负交界处大致在180°附近, 正值中心位于西太平洋的印尼、菲律宾等地, 负值中心位于赤道东太平洋; 这种以日界线为轴的东西高度场符号相反的特征与南方涛动类似。500 hPa高度场(图 3b), 赤道低纬度地区为正值带; 两侧的印度洋及亚洲大陆同为正值区; 但在太平洋上, 正值带两侧为沿纬线分布的负距平带, 这种从北到南的“–, +, –”高度场分布形式, 反映了赤道太平洋西风偏强的特征, 有利于厄尔尼诺的形成。850 hPa纬向风场(图 3c), 太平洋大部分地区盛行西风, 大值中心在赤道中太平洋附近, 有利于西太平洋的暖水向东积聚, 而赤道印度洋仍然为东风距平。850 hPa经向风场(图 3d), 太平洋上, 赤道北侧为负距平, 盛行北风, 而赤道南侧的太平洋中、西部为南风; 经向风应力距平是从两侧向赤道辐合的, 有利于表层暖水在赤道太平洋积聚, 并随西风移向东太平洋。500 hPa垂直速度场(图 3e), 伴随海表温度正值中心的东移, 上升区也移向赤道中、东太平洋, 与厄尔尼诺有很好的对应关系。最后来看到达地表的向下太阳辐射通量(图 3f), 太阳辐射是地气系统能量的总来源, 是气候形成的决定性因子; 众所周知, 由于地气系统吸收太阳辐射的区域不均性, 加之地球自转作用, 产生了三圈环流, 低纬度盛行偏东信风, 当赤道中、东太平洋入射的太阳辐射减弱, 将导致信风减弱, 甚至有可能出现风向翻转, 有利于厄尔尼诺的形成。

再来分析SVD第一模态左、右场时间系数间的时滞相关, 除850 hPa纬向风场为超前海温一月相关系数最大外, 其余各物理因子与海温耦合的相关系数均为同时相关最强。因此, 结合SVD分解的空间型与时间系数, 取各月与厄尔尼诺相关的预报因子(表 4)。

根据多变量均生函数的建模步骤, 建立了各月NINO3指数的预报模型, 并对2015~2017年进行了预报试验, 预报结果如下表:外推一步的预报, 与方案二预报结果类似, 预报与实测值非常接近, 但预报的厄尔尼诺开始时间较晚; 外推两步, 预报与实测值正负演变趋势基本一致, 对厄尔尼诺结束时间以及2016年下半年弱的拉尼娜都有较好的预报; 外推三步, 预报模型失去意义。

2.2 近17年厄尔尼诺的滑动预报试验

为了进一步探究均生函数及其两种改进模型在厄尔尼诺预报中的适用性, 本文在2015年强厄尔尼诺个例预报的基础上, 滑动选取35 a预报样本(例如对2006年外推两步的预报, 1970~2004年为建模时段; 2010年外推一步的预报, 则1975~2009年取为样本), 对2001~2017年进行预测, 给出整体预报评价。

对于气候预测来说, “趋势”预报正确尤为重要, 业务中往往将预报对象划分若干等级, 预报与实况级别一致即为预报正确。因而本文采用如下评分标准:预报与实况距平符号相同、或符号不同但二者绝对值均小于0.5℃, 则为预报正确, 其余情况为预报错误。表 6、表 7分别给出了三种模型外推一步、外推两步的相对误差及距平准确率, 可以看出:外推一步, 三种方案17 a平均的相对误差最小为1.19%, 最大4.26%, 各月间存在一定差异, 春季预报效果较好, 冬季较差; 整体来看, 两种改进方案不仅多年平均的相对误差有所减小, 而且预报准确率高于普通模型约5个百分点, 预报技巧有一定提升。外推两步, 模型的预报效果略有降低, 但三种方案平均的预报准确率均在65%左右, 仍有一定参考价值。

对于NINO3指数的预报, 能否正确预报出两类海温异常事件是检验模型优劣的关键。因此, 进一步分析三类模型对2001年以来5次厄尔尼诺和3次拉尼娜的预报情况(表 8、表 9), 表中事件的实际起止日期来源于国家气候中心; 预报值仅给出连续的结果, 中间允许1个月的间断, 非连续值即便达到厄尔尼诺/拉尼娜标准也不计入统计。由表中可以看出:尽管对个别事件的预报效果较差(预报效果欠佳的个例往往是三个模型都有较差的预报结果), 预报的异常事件持续时间明显短于实际持续时间, 甚至出现完全预报失败的情况; 但对于大部分事件, 模型都能给出具有一定参考意义的预报。整体而言两种改进模型的预报结果, 都明显好于普通均生函数模型。此外, 分析发现, 三种方案对中部型海温异常事件的预报效果均欠佳, 尤其是持续时间短、强度弱的个例。例如2004年的中部型厄尔尼诺和2011年的中部型拉尼娜, 三类模型或预报失败或预报的持续时间仅2~3个月, 与实况差别较大; 而2015年的东部型厄尔尼诺及2010年的东部型拉尼娜则有与实况非常接近的预报结果, 其中2015年的厄尔尼诺, 虽然气候中心定义其开始于2014年10月, 但此次事件增温首先出现于赤道中太平洋, NINO3区海温持续大于0.5℃则开始于2015年春季, 预报与实况是接近的。总地来说, 均生函数及其改进模型对厄尔尼诺和拉尼娜的预报是可行的, 两种改进模型也能有效提高模型的预报准确度。

3 结论

本文利用三种基于均生函数的预报方法, 首先对2015年的超强厄尔尼诺进行了预报试验, 为了进一步检验模型的适用性, 还建立了滑动的预报模型对2001~2017年的NINO3区海温指数进行预报, 主要得到以下结论:

(1) 均生函数模型对于极端值的拟合和预测效果较好, 试验表明利用均生函数进行厄尔尼诺/拉尼娜的气候预测是可行的, 模型具有一定实用价值; 基于EEMD的均生函数预测方法和多变量均生函数预测方法相对于普通模型, 预报技巧有了一定的提高, 尤其是外推一步以内, 是改进一般均生函数模型的有效手段。

(2) 研究中发现, 三类模型对于中部型厄尔尼诺/拉尼娜事件的预报效果较差, 尤其是对于持续时间短、强度较弱的个例。这主要是因为文章中仅考虑了NINO3区的海温变化。在今后的分析中, 区分厄尔尼诺类型, 并分别建立预报模型, 是改进厄尔尼诺预报的有效手段。

(3) 一般国际上广泛采用的数值模式, 通常只能提前6个月对厄尔尼诺给出较为正确的预报, 而应用统计预报方法, 外推两年以内的预报具有一定的参考意义, 可以适当延长厄尔尼诺的预报时效; 均生函数模型具有进行多步预测的优点, 可以制作滚动的预报并应用于业务中, 为厄尔尼诺预报提供新的思路与参考。