2.2 径流序列的混沌特性判别原理及方法

2.2.1 功率谱分析

功率谱分析是分析时间序列常用的方法，它通过快速的傅里叶变换研究频率f与自相关函数C（τ）的傅里叶变换之间的关系，根据功率谱图的特点直观地揭示离散数据系列的周期性^{［207，208］}。通常，系统的观测资料在时间历程图可能不规则，但其功率谱图像却可能呈现出某种相似性或具有一定的规则性。因此，通过对时间序列数据进行变换（或逆变换），计算其自相关函数、功率谱，然后在频率域上绘出谱图，通过功率谱图的特点可以了解时间序列的周期性。一般地，不同性质的时间序列的功率谱是不同的，周期性序列的功率谱图具有单峰或几个峰，而非周期混沌运动的功率谱则和噪声序列一样无周期性，其谱图无明显的峰值或峰值连成一片。有些混沌系统，其功率谱超过某一定值后，谱线随频率指数衰减，所以功率谱分析可以定性地反映序列的混沌特征。

一般自相关函数C（τ）为

式中：τ为时间的移动值。

当τ=0时，C（0）表示系统的功率；当τ≠0时，C（τ）的值表示相差为τ的两时刻运动的相互关联或相似的程度。当x（t）的幅值一定时，C（τ）越大，就意味着x（t）与x（t+1）也越相似。而τ越大时，x（t）与x（t+1）差别可能越来越大，最后变成完全无关，C（τ）越来越小直至趋于零。因此，C（τ）可以用于表示时间序列的规律性或系统的确定性程度。为了表示运动的频谱特性，对式（2.1）进行傅里叶变换，即

C（τ）的傅里叶变换S（ω）即表征运动的频谱特性，S（ω）表示功率的（频）谱密度。

功率谱法是一种定性的判定方法，它只能从图形上作直观判断，使用时很难掌握一个明确的标准，难以应用于在线实时的场合。然而，由于功率谱能够把带有噪声看起来很复杂的多周期行为与真正的混沌区别开来，并且计算简单，结果直观，因此，常将它作为混沌特性的初步判别。

2.2.2 BSD统计量

BDS统计^［209］是Brock、Dechert和Scheinkman在关联积分的基础上建立起来的用于检验时间序列是否具有独立同分布的统计量，是一种识别时间序列中确定性混沌特征或者非线性随机特征的较为有效的方法，它能从本质上检验时间序列的非线性依赖特征。具体来说，对时间序列{x₁，x₂，…，x_N}，运用相空间重构方法在m维嵌入空间构造一个时间延迟为τ的嵌入向量X_i=（x_i，x_i+τ，x_i+2τ，…，x_{i+（m-1）τ}），i=1，2，…，N-（m-1）τ，定义其关联函数为

其中

式中：H为Heaviside函数。

Brock、Dechert和Scheinkman在关联积分的基础上提出的BDS统计量定义为

其中，。可以证明在时间序列x={x_i|i=1，2，…，N}是独立同分布 （iid）的原假设H₀ 下，BSD统计量的极限分布为标准正态分布，即

因此，若时间序列为iid序列，则BDS渐进收敛于标准正态分布；反之，则可拒绝原时间序列为iid的假设H₀，从而可判定在该时间序列中可能存在非线性结构。

值得注意的是，为了使用BDS统计量来检验数据中的非线性结构，先要用一个线性自回归模型AR去拟合原时间序列{x₁，x₂，…，x_N}，消去原序列中的线性相关关系，然后再对AR拟和后的残差序列使用BDS统计量检验其是否为独立同分布序列，如果检验结果拒绝残差是独立同分布的，则意味着数据具有内在非线性。

首先用自回归模型提取径流序列中的线性相关成分，模型定阶采用FPE准则，即最小最终预测误差的定阶准则。采用Ljung-BoxQ统计量Q（p）检验AR拟合残差序列的独立性的，记为

式中：为所检验时间序列滞后期为k的样本自相关系数；N为样本容量；p为滞后阶数。

若统计值Q（p）小于检验临界值，则表明不能拒绝残差序列是独立的原假设，从而可以认为残差序列已不包含线性相关成分；否则，表明序列存在相关成分。

实际计算中，用C（1，N，r）取代C（r），K（N，r）取代K（r），得到统计量σ^²（m，N，r），以σ^²（m，N，r）作为σ²（m，N，r）的一致估计。由于BDSL^［210］方法存在计算速度慢、在样本较小时误差大等缺陷，Kanzler Ludwing^［211］对BDSL算法进行改进。本书所用的BDS检验是采用Kanzler Ludwing算法。在进行检验时，考虑到BDS检验为双边检验，因此，在5%的显著水平下，若检验统计量超过1.96或者小于-1.96，则拒绝原假设，认为序列并非独立同分布，该序列存在非线性结构特征。

BDS统计量检验的原假定是数据生成过程为独立同分布过程。若拒绝该原假定表明序列不是独立的，然而我们却不能进一步根据该统计量判断这种非独立的具体形式是非线性相关还是混沌，它只能够对序列的非线性依赖提供间接证据。因此，确定序列不独立的具体形式是进一步的任务。

2.2.3 Cao方法

Cao方法^［212］是一种直观而简洁的判定时序数据非线性混沌特性的方法，从本质上来说，它是对虚假邻域算法的改进。我们知道，在虚假邻域算法中，进行虚假近邻点判断时用到的两个阈值，一般根据经验给出，具有较大的随意性。从理论上来说，相空间中每个不同的点，判断它的近邻点是否是虚假的时，需要不同的阈值，算法只是根据经验给出一个统一的阈值，这与实际情况有一定程度的偏差。为了解决这一问题，Liangyue Cao提出了改进的方法，该方法定义了两个参数E₁（m）和E₂（m），其中，E₁（m）用于确定最小的嵌入维数，而E₂（m）用于判别一个时间序列是否具有分维特性的混沌序列。

对时间序列{x₁，x₂，…，x_N}和延迟时间τ，运用相空间重构方法我们可以构造一个m维状态空间向量X_i（m）=（x_i，x_i+τ，x_i+2τ，…，x_{i+（m-1）τ}），i=1，2，…，N-（m-1）τ。与虚假邻域算法相类似我们定义如下统计量：

式中：为某种意义下的距离量度，本书采用无穷范数，见式 （2.9）；n（i，m）为在上述距离定义下与m维状态空间第i个向量距离最近的点的标号，1≤n（i，m）≤N-mτ。

用E（m）表示嵌入维数为m情况下时间序列各点与其邻近点之间距离的统计平均，即

为了观察E（m）随嵌入维数m的变化情况，定义如下统计量

从式（2.11）中可以看出，如果时间序列具有混沌特性，则当嵌入维数m达到最小嵌入维数时，E₁（m）将保持不变，此时的嵌入维数m就是我们所寻找的最小嵌入维数。

为了区分出随机噪声序列，我们定义了另外一个统计量E₂（m）。以E^*（m）表示从嵌入维数m变化到m+1的时距离增加的统计平均，即

式中各项意义如前。

对于随机噪声而言，各个时刻的观测值是相互独立的，故E^*（m）的大小和嵌入维数m没有关系，即E^*（m）在各维情况下大致相等；而对于其他混沌信号，随着m的增加，假邻近点数目的不断减少会导致E^*（m）随之减小，直到增大到最小嵌入维数m时，E^*（m）的大小将不再随着m变化而变化，因此我们可以根据E^*（m）变化情况判断待判序列是否为混沌信号。定义统计量E₂（m）为E^*（m）在相邻维上的相对值，即

对于随机时间序列，原则上来讲，E₁（m）不会随着嵌入维数m的增加而达到一个稳定值，但是，在实际计算中，当m足够大时，很难判断E₁（m）有没有增加，事实上，由于实际观测样本有限，随机序列也可能在m达到一定值后停止变化。这个时候我们可以通过E₂（m）值来辨别随机噪声。对于任意的嵌入维数m，如果对应的E₂（m）值都等于1（或在1附近），则表明该时间序列是随机的；而对于非线性混沌序列，E₂（m）和嵌入维数m之间存在某种关系，不可能对任何嵌入维数m，E₂（m）都等于1，一般对确定性混沌序列E₂（m）会渐渐趋向1。因而，可以通过嵌入维数m与函数E₂（m）的曲线图，直观地判定时间序列是否具有非线性混沌特性。

2.2.4 Lyapunov指数

Lyapunov指数是描述混沌运动对初值条件的敏感性吸引子特征量，它从整体上定量描述轨道的平均分离或者收缩的快慢，表征了系统的混沌性质^［207］。对于一般的n维动力系统，定义Lyapunov指数如下：设F是Rⁿ→Rⁿ上的n维映射，若系统的初始条件用一个无穷小的n维的球表示，随着时间的演变系统变为椭球。将n维椭球的n个主轴按其长度顺序排列，λ₁≥λ₂≥…≥λ_n，那么第i个Lyapunov指数根据第i个主轴的长度P_i（n）的增加速率定义为

n个Lyapunov指数表示了系统在相空间的n个方向的收缩或者扩张的性质。最大的Lypunov指数决定了轨道发散覆盖整个吸引子的快慢，最小的Lyapunov指数决定了轨道收缩的快慢，而所有Lyapunov指数之和大体上表征了轨道总的平均发散快慢。

若λ_i<0则意味着相邻点最终要靠投合并成一点，这对应于稳定的不动点和周期运动；若λ_i>0则意味着相邻点最终要分离，这对应于轨道的局部不稳定，如果轨道还有整体的稳定因素（如整体有界、耗散、存在捕捉区域等），则在此作用下反复折叠并形成混沌吸引子。因此，区分不同类型的吸引子的重要标志就在于混沌吸引子具有正的Lyapunov指数。故λ_i>0可作为系统混纯行为的一个判据。若系统最大Lyapunov指数λ₁>0，则系统一定是混沌的。所以，通常将最大Lyapunov指数λ₁为正作为判断混沌性质的重要条件。

目前，计算混沌时间序列的最大Lyapunov指数最常用的方法是Wolf法^［213］，但这一方法适用于时间序列无噪声、切空间中小向量的演变高度非线性，对于含动态噪声的数据序列或是数据系列较少的情况得到的结果并不理想。因此，本书采用Rosenstein等^［214］提出一种根据小数据量来计算最大Lyapunov指数。与Wolf法作对比，该算法很容易应用，无需计算向量间的夹角，充分利用了所有的相点，而且对重构相空间的参数变化以及相应的噪声具有鲁棒性。具体步骤如下：

（1）对时间序列{x_i，i=1，2，…，n}进行快速傅立叶变换（FFT），计算出平均周期T。

（2）找相空间中每个点X_j的最近邻点X_j^，并限制短暂分离，即

（3）对相空间中每个点X_j，计算出该邻点对的第i个离散时间步的距离为

（4）对每个i，求出所有j的lnd_j（i）平均值y（i），即

式中：p为非零d_j（i）的数目。

用最小二乘法作出回归直线，该直线的斜率就是最大Lyapunov指数。

2.2.5 关联维数

关联维数D₂是分形的维数中最具代表性的一种，可对吸引子的几何特征及基于吸引子的轨道随时间演化情况进行数量上的描述，是序列混沌特征的重要不变量。关联维数的取值是系统混沌与否的一个表征。一般地，确定系统的关联维数是一个1～3之间的整数；随机系统的维数为无穷大；而混沌系统的维数为一个正的分数，所以称之为分维。

从单变量时间序列（x₁，x₂，…，x_N）出发，构造一批m维矢量支起一个嵌入空间。只要嵌入维数m足够高（m≥2D+1，D是吸引子维数），就可以在相差拓扑变换的意义下恢复原来的动力系统。1983年，Grassberger和Procaccia^［49］提出了从时间序列计算吸引子关联维数的G-P，该算法具有重要的理论意义和实践意义。对于从m维嵌入空间重构的混沌动力系统，构造m维矢量的办法常采用时间延迟法，即按间隔τ从时间序列（x₁，x₂，…，x_N）中取数作为分量：

式中：M=N-（m-1）τ。构造好矢量Y_i之后，定义它们之间的距离：

凡是距离小于给定数r的矢量，称为有关联的矢量，若一共构造了M个矢量，数一下有多少对关联矢量。它在一切可能的M²种配对中所占比例称为关联积分C（r，m）：

其中H为Heaviside单位函数：

当r→0时，C（r，m）与r存在以下关系：

式中：D₂为关联维，随着m的增加，极限值D₂也增加并稳定于某一有限值。

在实际计算中一般是通过相空间重构和G-P算法，得到反映嵌入维数与关联维数的曲线图；从小到大改变m的值，当m大于等于一个整数时，观察D是否保持一个常数，即双对数坐标lnC（r，m）-lnr中的直线段，由此可以判定系统具有混沌系统的分维特征。分维数除了标志着混沌过程结构的自相似规律外，并不能完全揭示出产生相应结构的动力学特征，要证明系统的混沌性质还必须借助其他方法。

2.2.6 Kolmogorov熵

Kolmogorov熵又称测度熵，是动力系统相空间轨道分裂数目渐近增长率的一种量度，它实际表征系统信息的损失程度，是精确化的信息熵，常用来度量系统运行的混沌或无序的程度^{［207，208］}。和Lyapunov指数、关联维数一样，Kolmogorov熵也不是混沌系统特有的，它的取值可以表征系统的混沌特征，但不是系统混沌性的充分条件。对于守恒系统，Kolmogorov熵为0；对于混沌系统，Kolmogorov熵为0到正无穷之间的数；对于随机系统，Kolmogorov熵趋于无穷大。

事实上，Kolmogorov熵是q阶广义熵（Renyi熵）的一个特殊情形。Kolmogorov熵定义为

q阶Renyi熵的定义为

Grassberger和Procaccia提出的K₂熵与关联积分存在如下关系：

当固定时间延迟τ，嵌入相空间维数分别为m和m+Δm时，有

一般取Δm=1，从理论上讲，当m→∞时，K₂→K。实际上，当m达到一定值时，K₂更趋于稳定，可将此时相对稳定的K₂作为K的估计值

因此，根据前述探讨，在嵌入维数按等间隔不断增加的情况下，对双对数坐标（ε）-lnε中的点，在无标度区间内，作等斜率线性回归，可同时得到关联维和K的稳定估计。使用K值可判断系统运动的性质，K=0，表明系统是规则运动的；K等于正无穷，表明系统做随机运动；0＜K＜+∞，表明系统做混沌运动，K值越大，那么信息损失速度越大，系统的混沌程度越大，系统也就越复杂。