1 引言

近年来，气候系统的非平稳特征已经得到了越来越多的证实和关注(Trenberth, 1990; Tsonis, 1996; Schmutz et al., 2000; Yang et al., 2000; Slonosky et al., 2001；杨培才等, 2003)，而作用于系统的外强迫随时间的变化是导致其非平稳行为产生的根本原因。

目前从数据出发，直接从非平稳时间序列中提取外强迫信号的方法有以下两种，一种是由Verdes et al.（2001）建立的，在交叉预报误差的基础上的反演非平稳系统外强迫因子变化方法；另一种是Wiskott and Sejnowski（2002）建立在奇异谱分析基础上的慢特征分析法（Slow Feature Analysis, SFA）。文中使用Wiskott的方法来对外强迫信号进行提取。

慢特征分析法起源于神经生物学，它的目的是从一个快变的非平稳时间序列中提取缓变的外强迫信号(Wiskott and Sejnowski, 2002)。目前，SFA方法在物理学、经济学等众多技术科学领域均有一定应用。Wiskott (2003a)利用SFA方法从设计的非平稳系统（基于Tent映射与Logistic映射）中提取外强迫信号，将提取出的外强迫信号与真实外强迫进行比较，二者相关系数可达到0.9。Konen and Koch（2009）通过一些试验发现，SFA方法只能给出一个变化最慢的信号或信号分量的组合。潘昕浓等（2017）利用SFA方法对具有层次结构的非平稳系统进行外强迫信号的提取，取得了不错的效果。可以说，SFA方法在理想的一维离散模型中取得了成功的应用，对于这些理想非平稳系统外强迫信号的研究也为认识气候系统的驱动机制提供了新的思路。

Yang et al. (2016)用慢特征分析法重建了北半球月平均气温的外强迫信号，并结合小波变换技术对外强迫信号的尺度特征与物理背景进行了分析，从中发现了太阳11年周期(the Hale cycle)与大西洋年代际振荡(Atlantic Multidecadal Oscillation, AMO)的谐波分量。Wang et al. (2016)基于Arosa臭氧观测资料，利用SFA方法提取出了外强迫信号，认为太阳活动和北大西洋涛动(North Atlantic Oscillation, NAO)对臭氧变化有显著的影响作用。陈潇潇等(2015)与Wang and Chen (2015)将SFA方法应用到臭氧与气溶胶的预测中，建立了包含外强迫信号的预测模型，提高了预测精度。

然而，不管在理想模型还是实际应用中，目前基于SFA方法的研究都只聚焦于对一维映射模型的外强迫信号提取，并得到了成功应用。而对于二维甚至更高维的非平稳系统，SFA方法的适用性值得进一步地试验及讨论。

本文以常见的二维Henon非线性映射为基础，构造更为复杂的非平稳系统，利用慢特征分析法和小波变换技术，对系统产生的非平稳信号进行外强迫信号的重建，对SFA在二维映射中的应用进行尝试；并且利用SFA方法，以北京月平均气温为实际参考序列进行外强迫信号的重建，结合小波变换技术对其尺度结构与可能的物理机制进行简要分析。

2 慢特征分析法

慢特征分析法(Slow Feature Analysis, SFA)旨在从一个已知的非平稳时间序列中提取出变化最慢的分量。它是一种无监督的算法，可以在一般的有限维度函数空间中找到最优函数的最优集合，并且能依据输入信号的大小和尺寸进行合理有效的模拟(Berkes and Wiskott, 2005)。SFA结合状态空间重构理论(Packard et al., 1980)与嵌入定理(Takens, 1981)的思想对函数空间进行重构，经过标准化处理、非线性扩展、球化处理以及主成分分析方法，将已知的信号投影到其中变化最慢的方向上，这个投影分量即可看成原信号的外强迫因子(Wiskott, 2003b)。有关SFA方法的详细的解释可参见(Berkes and Wiskott, 2005)，下面给出慢特征分析法主要的算法步骤：

首先，给定一个非平稳时间序列${\{ x(t)\} _{t = {t_1}, {t_2}, \ldots, {t_n}}}$（n为序列的长度，t为时间），对非平稳时间序列{x(t)}进行空间重构（设置状态空间的嵌入维数为m，时滞参数τ=1）：

$ \mathit{\boldsymbol{X}}(t) = {\{ {x_1}(t), \ldots, {x_m}(t)\} _{t = {t_1}, {t_2}, \ldots, {t_N}}}, $

(1)

其中重构后序列的长度N=n-m+1。

然后，利用X(t)中的一次项和二次项对函数进行非线性扩展，构造一个k维的函数空间H(t)：

$ \begin{array}{l} \mathit{\boldsymbol{H}}(t) = \{ {x_1}(t), \ldots, {x_m}(t), {x_1}(t){x_1}(t), \ldots, \\ \;\;\;\;\;\;\;\;\;\;\;\;{x_1}(t){x_m}(t), \ldots, {\rm{ }}{x_m}(t){x_m}(t){\} _{t = {t_1}, {t_2}, \ldots, {t_N}}}, \end{array} $

(2)

方便起见，记为

$ \mathit{\boldsymbol{H}}(t) = ({h_1}(t), {h_2}(t), \ldots, {h_k}(t)), $

(3)

其中k=m+m(m+1)/2。对H(t)进行标准化与正交化处理，得到

$ \mathit{\boldsymbol{Z}}(t) = {\{ {z_1}(t), {z_2}(t), \ldots, {z_k}(t)\} _{t = {t_1}, {t_2}, \ldots, {t_N}}}, $

(4)

得到的Z(t)满足ZZ^T=1，且Z=0。此时，Z(t)中的每一个分量都可以由z_j的线性组合表示：

$ y(t) = {a_1}{z_1}(t) + {a_2}{z_2}(t), \ldots, {a_k}{z_k}(t), $

(5)

用${\dot z_j}({t_i}) = {z_j}({t_{i + 1}}) - {z_j}({t_i})$表示一阶导数，构建Z(t)的导函数空间：

$ \mathit{\boldsymbol{\dot Z}}(t) = {\{ {\dot z_1}(t), {\dot z_2}(t), \ldots, {\dot z_k}(t)\} _{t = {t_1}, {t_2}, \ldots, {t_N}}}. $

(6)

最后，对矩阵$\left\langle {\mathit{\boldsymbol{\dot Z}}{{\mathit{\boldsymbol{\dot Z}}}^T}} \right\rangle $进行主成分分析，得到$\left\langle {\mathit{\boldsymbol{\dot Z}}{{\mathit{\boldsymbol{\dot Z}}}^T}} \right\rangle $的特征值λ_j和特征向量，进行标准化处理后得到权重向量w₁, …, w_j。对相应权重系数的输出函数

$ {g_j}(t) = {w_j}\mathit{\boldsymbol{Z}}(t), $

(7)

进行积分，求得输出信号：

$ {y_j}(t) = r{w_j}\mathit{\boldsymbol{Z}}(t) + c, $

(8)

获得对应的信号分量，其中r与c均为常数。当选择最小的特征值λ时，带入对应的最小特征的权重向量，积分后可以得到变化最慢的信号分量，也就是所说的外强迫信号。可以看到，由于公式（8）中的r与c为常数，故SFA方法得到的外强迫与实际外强迫只相差一个平移因子和一个放大因子。

3 模型试验

Henon映射是二维空间中产生混沌的一种迭代映射，它只含有一个非线性项，所以它是高维之中最简单的非线性映射（迟洪钦和吴忠英，1994）。Henon映射的动力学方程为

$ \left\{ \begin{array}{l} x(t + 1) = 1 - ax{(t)^2} + by(t), \\ y(t + 1) = x(t), \end{array} \right. $

(9)

取参数a=1.4、b=0.3，可令系统处于混沌体制下。

Lyapunov指数是用来判断系统是否存在混沌的判据，若系统存在一个Lyapunov指数大于0，则系统一定处于混沌状态；若系统的所有的Lyapunov指数均小于0，则系统为周期运动。本文也利用此方法判断系统模型是否处于混沌状态。图 1是Henon映射中参数a的部分区间解，当最大Lyapunov指数大于0时对应的参数a的值即可使Henon映射处于混沌状态。

本文利用两组单向耦合的Henon映射建立非平稳模型，一组为较为简单的Henon映射${\{ x(t)\} _{t = 1, 2, \ldots, N}}$：

$ \left\{ \begin{array}{l} {x_1}(t + 1) = 1.4 - a{x_1}{(t)^2} + b{x_2}(t){\rm{, }}\\ {x_2}(t + 1) = {x_1}(t){\rm{, }} \end{array} \right. $

(10)

另一组为经过改进的Henon映射{y(t)}：

$ \left\{ \begin{array}{l} {y_1}(t + 1) = 1.4 - (ks(t){y_1}(t) + (1 - k){y_1}{(t)^2}) + b{y_2}(t),\\ {y_2}(t + 1) = {y_1}(t){\rm{, }} \end{array} \right. $

(11)

其中s(t)为外部控制参数。

在下面的试验中，将基于这两组Henon映射构建非平稳系统模型。在这里要说明的是，虽然Henon映射是一个二维映射，但由于它第二维的数组相对第一维数组独立，仅比第一维数组滞后一项，故在试验中，仅选择第一维数组进行外强迫信号的提取。

3.1 单时变参数Henon映射模型试验

在单时变参数Henon映射模型试验中，我们构造了3个变化规律不同的外强迫信号，分别带入到系统模型中，产生非平稳时间序列。利用SFA方法对序列进行外强迫提取。

（1）设置时变参数${\{ a(t)\} _{t = 1, 2, \ldots, N}}$为

$ a(t) = - \cos (2{\rm{ \mathsf{ π} (}}t/{T_1}{\rm{))exp(}}t/{T_2}). $

(12)

利用{a(t)}来构造一个缓变的外强迫信号，并带入公式（11）中，其中T₁=500，T₂=2500，y₁(1)=y₂(1)=0，b=0.3，k=0.02。

将{a(t)}迭代5000次并复制，获得长度为10000的序列${\{ {s_1}(t)\} _{t = 1, 2, \ldots, 10000}}$（图 2a），并将{s₁(t)}作为时变参数代入公式（13）中迭代10000步，得到${\{ y(t)\} _{t = 1, 2, \ldots, 10000}}$。通过计算，获得{y₁(t)}的一个大于0的Lyapunov指数λ=0.3429，说明系统处在混沌状态。取{y₁(t)}的后2000个数据作为试验的非平稳时间序列（图 3b），利用SFA方法对此序列进行外强迫提取，嵌入维数m=13，时滞系数τ=1。将提取出的外强迫信号，记为{as₁(t)}，并与真实外强迫信号{s₁(t)}比较，发现二者相差一个放大因子，所以对提取出的外强迫信号{as₁(t)}进行简单处理，扩大4倍后再与真实外强迫进行比较，二者的趋势一致，有较好的一致性，相关系数可以达到0.98（图 3c）。

（2）利用Henon映射模型[公式（10）]来构造一个跃变的类方波信号作为外强迫信号，带入到Henon映射模型[公式（11）]中，其中，x₁(1)=x₂(1)=0, y₁(1)=y₂(1)=0, b=0.3, k=0.1。

令{x(t)}式迭代10000次，对{x₁(t)}取9901~10000步，令每100步取相同数值，构造长度为10000的序列${\{ {s_2}(t)\} _{t = 1, 2, \ldots, 10000}}$（图 4a），将{s₂(t)}作为时变参数代入公式（15）中迭代10000步。通过计算{y₁(t)}的Lyapunov指数λ，存在一个λ=0.3794(λ>0)，证明{y₁(t)}此时是一个混沌系统。取后2000步作为试验时间序列（图 5b），对其进行外强迫提取，设置嵌入维数m=13、时滞系数τ=1。结果表明，提取出的外强迫信号{as₂(t)}和真实外强迫{s₂(t)}相差一个平移因子，利用SFA提取的外强迫信号较真实的外强迫在跃变处有一定的迟滞，二者相关系数达到0.98（图 5c）。

（3）用Henon映射模型[公式（10）]来构造外强迫信号，与试验2不同，本次试验中使用的外强迫信号变化缓慢，跃变较小。将外强迫带入到Henon映射模型[公式（11）]中，其中，x₁(1)=x₂(1)=0, y₁(1)=y₂(1)=0, b=0.3, k=0.1。

令公式（16）迭代10000次，取{x₁(t)}的9901~10000步，先进行插值延长至1000步，再令每10步取同一个数值，构造长度为10000的序列${\{ {s_3}(t)\} _{t = 1, 2, \ldots, 10000}}$（图 6a），将{s₃(t)}作为时变参数代入公式（17）中迭代10000步，通过计算，获得{y₁(t)}的一个为正数的Lyapunov指数λ=0.2433，说明{y₁(t)}处于混沌状态。取{y₁(t)}的后2000步作为试验时间序列（图 7b），通过SFA方法对其进行外强迫提取，嵌入维数m=13、时滞系数τ=1。分别对获得的外强迫信号{as₃(t)}和真实外强迫{s₃(t)}进行归一化处理并比较，虽较真实的外强迫有一些噪声，但总体趋势一致，二者的相关系数达到0.99（图 7c）。

3.2 双时变参数Henon映射模型试验

在双时变参数Henon映射模型试验中，由于SFA方法只能给出一个变化最慢的分量或一个变化最慢的分量组合，所以我们将结合Morlet小波变换(Torrence and Compo, 1998)尝试对模型中的两个时变参数进行还原。

小波变换（wavelet transform）是建立在傅里叶变换基础上的一种变换分析方法，它可以对信号进行时间和频率的局部变换，有效地将时频信息提取出来，反映信号在频域与时域上的尺度特征，尤其适用于非平稳信号（胡广书，2004）。

设定双时变参数：

$ a(t) = 0.8 + 0.1\cos (4{\rm{ \mathsf{ π} (}}t/{T_1}{\rm{)}}), $

(13)

$ b(t) = 0.3 + 0.09\cos (11{\rm{ \mathsf{ π} (}}t/{T_1}{\rm{)}}), $

(14)

并带入到Henon映射模型[公式（10）]中，其中，T₁=500, x₁(1)=x₂(1)=1。

将参数{a(t)}、{b(t)}代入公式（20），迭代10000次，获得非平稳序列{x₁(t)}（图 8c）。通过计算，发现{x₁(t)}存在一个为正的Lyapunov指数λ=0.2238，系统处于混沌状态。取后2000步作为试验时间序列（图 9c），使用SFA方法对试验序列进行外强迫提取，取嵌入维数m=13、时滞系数τ=1。对获得的外强迫信号{as(t)}进行小波分析，观察外强迫信号的时间平均功率谱，存在两个通过检验的峰值as₁和as₂（图 10a），通过滤波提取出这两个峰值对应的信号分量{as₁(t)}、{as₂(t)}，将提取出的两个信号进行归一化处理并分别与标准化后的{a(t)}、{b(t)}进行比较，相关系数分别为0.99和0.98（图 10b和10c）。

以上两组试验表明，慢特征分析法可以有效地从以Henon映射为基础构造的非平稳系统中提取出外强迫信号。在单时变参数的非平稳模型中，SFA方法可以直接提取出外强迫，并且与真实外强迫相比只相差一个放大因子和一个平移因子；在双时变参数的非平稳模型中，可结合小波变换技术还原出两个时变外强迫序列。在下文中，将选取真实的非平稳系统序列进行外强迫的提取与分析。

4 北京市气温时间序列的外强迫分析

通过利用慢特征分析法重建北京气温的外强迫信号，并利用小波变换技术对外强迫信号的尺度特征与物理机制进行简要分析。资料来自于NCEI (National Centers for Environmental Information)的全球站点月平均数据(https://gis.ncdc.noaa.gov/maps/ncei[2017-06-01])，时间跨度从1951年1月到2012年12月，共744个月。

利用慢特征分析法对北京月平均气温时间序列进行外强迫信号的提取，参照Yang et al. (2016)的工作，取嵌入维数m=13、时滞系数τ=1，外强迫重建结果如图 11b所示。

利用Morlet小波变换技术对提取出的外强迫信号进行分析，得到了小波变换实部和时间平均功率谱如图 12和图 13。小波变换的实部图揭示了外强迫的尺度结构及其随时间的变化，宏观地反映了外强迫在频域和时域的主要特征，图中可以看出，北京气温的外强迫信号存在着明显的周期变化（图 12），其38~45 a尺度振动周期较为明显；20~22 a尺度周期在1980~2000年间有所减弱，2000年后重新开始加强；12~14 a尺度周期在1970~1995年间存在波动，下面将结合时间平均功率谱对其外强迫信号进一步分析。

图 13为其外强迫的时间平均功率谱，从图中可以看到存在着7个峰值，我们用S_n(n=1, 2, …, 7)来表示这7个特征尺度，用L_n来表示它们的特征尺度的周期，f_n来表示频率。其中L₆、L₇两个周期通过了95%置信度的置信检验。它们对应的周期分别为L₆=21.3 a、L₇=42.7 a。L₆与已知的太阳轨道运动长周期（the 22 years rotation cycle）一致（刘复刚等，2013），而L₄恰好是L₆的2倍关系，考虑是太阳活动对气温的影响。太阳质子耀斑也存在着明显的周期（the Schwabe cycle）变化（周树荣等，1996），与成分S₄的周期L₄=11.6 a对应。通过对7个分量频率的分析，发现了各组分之间存在着一定的谐波关系（如表 1）。假设f₄、f₆这两个频率为基本分量，发现其它分量可以由这两个分量的整数倍或整数系数的线性组合来表示，这表明S₁、S₂、S₃和S₅是S₄与S₆的谐波分量。

5 结论与讨论

本文将SFA方法应用于二维非平稳时间序列外强迫信号的重建分析，结果表明，慢特征分析法可以有效地将外强迫信号从Henon映射构造的复杂非平稳系统中提取出来，并结合小波变换技术对外强迫信号进行分析，本文的主要结论如下：

（1）在单时变参数Henon映射模型的试验中，慢特征分析法可以从非平稳的Henon映射模型中提取出外强迫信号；提取出的外强迫信号与真实外强迫的相关系数可以达到0.98，但它们之间相差一个放大系数和一个平移系数，与Wiskott（2003a）等结果相同。

（2）在双时变参数Henon映射模型的试验中，先利用慢特征分析法从中提取出的外强迫信号，接着通过小波变换技术对其进行分析，得到两个分量，分别与已知驱动进行比较，相关系数均达到了0.98。

（3）作为对真实气候系统的应用，本文对62年的北京月平均气温进行外强迫的重建，并简要分析了外强迫的尺度特征与物理机制。其外强迫尺度特征存在2.7 a、3.2 a、5.8 a、11.6 a、13.8 a、21.3 a、42.7 a周期。其中21.3 a对应太阳轨道运动的周期，11.6 a对应太阳黑子耀斑的周期，42.7 a为太阳轨道运动周期的2倍关系，它们代表了太阳活动对气候系统的影响，其他一些分量的物理机制值得在随后的研究中深入探讨。

在SFA方法中，嵌入维数m与时滞参数τ的设定对相空间重构结果有很大的影响。由于在本文中构建的系统为非平稳系统，所以对于以Takens定理为基础的传统求取嵌入维数m跟时滞参数τ的方法：如G-P算法跟自相关函数法将不再适用，需要对嵌入维数m与时滞参数τ进行大量的尝试进而选取。外强迫对气候影响不可忽视，对外强迫的研究与应用也应该得到更多的关注。将外强迫加入预测模型以提高预测精度的研究已经取得很好的效果；从气候要素中提取外强迫，进而从外强迫中分析寻找影响气候的因果关系还需要更多探究，为进一步了解气候变化，甚至利用外强迫建立区域模式做一些贡献。

将SFA方法引入对非平稳气候系统驱动力的研究，还处于尝试阶段。慢特征分析法的扩展及应用需要更多的试验和研究。SFA方法能否适用于其他高维复杂非平稳系统，还需要一些试验。由于SFA方法给出的是一个外强迫分量的组合，所以在实际应用中，对于这样的外强迫分量组合的分析，还需要更多的探索和依据。