2.3 应用实例

径流系统一般都是确定性、随机性和混沌性的统一体。因此，描述径流系统演化特性，既可以用确定性的方法，也可以用非线性动力学方法，但只有将两者有机地结合起来才能客观、真实地认识和描述径流系统演化特性。以河北省南水北调受水区邯郸市漳河观台水文站月径流量为例，对该站的1963—2000年月径流量系列资料进行分析。考虑到资料获得的难易，选取降雨量作为主要影响因素，构成径流系列的多变量时间序列X₁，X₂，…，X₄₅₆，其中X_n=（x_1，n，x_2，n，x_3，n，x_4，n），n=1，2，…，456，x_1，n为月径流量时间序列，x_2，n、x_3，n、x_4，n分别为匡门口、白土和观台站降雨量时间序列。根据2.2节所述的混沌特性判别方法，对观台站降雨径流系列进行混沌特性分析。

2.3.1 功率谱定性分析

首先，从定性的角度来分析观台站月径流降雨时间序列是否具有混沌特性。由上述实际观测资料根据式（2.1）、式（2.2）计算各时间序列的功率谱，绘制功率谱图。观台站月径流量及其附近各雨量站降雨时间序列的功率谱图如图2.1所示。观察月径流量、降雨时间序列的功率谱图可以看出，4个站点的功率谱有连续的谱线，且均为宽带的连续谱，但在某些频率附近均具有明显的峰值，像具拟周期的图，从形状上总体看接近混沌序列，初步分析该降雨径流时间序列可能具有一定的混沌特征。

当然，由于实测数据有限，仅从功率谱图来判定时序数据是否是混沌序列或长周期序列还不十分精确，因此需用其他方法定量分析计算，以确定其是否具有混沌特征。以下我们再用BDS统计量、E₂（m）、Lyapunov指数、关联维数以及Kolmogorov熵对降雨径流系统进行混沌特性的定量分析。

图2.1 降雨径流时间序列功率谱图

2.3.2 BSD非线性检验

首先，我们对降雨径流时间序列进行BDS检验。在使用BDS进行非线性检查前，需要消除原时间序列自相关的影响，通常对原时间序列拟合AR（p）模型，在寻找到合适的阶数p后，计算AR（p）的残差序列并对该残差序列进行BDS检验。BDS统计量满足渐近的标准正态分布，零假设为该残差序列是iid的。若计算的统计量大于1.96，则在5%的显著性水平下拒绝数据是iid的原假设，如果结果拒绝零假设，则意味着原时间序列在该显著水平下是非线性的。给定显著性水平α=0.05下，按式（2.7）计算原序列的Ljung-BoxQ统计量Q（p），结果见表2.1，括号内为该统计量的p值（均大于0.05）。Ljung-Box检验结果表明，所有检验统计值在给定显著性水平α=0.05下均较为显著，因而，降雨径流序列存在线性相关。对原时间序列用AR（p）模型拟合，再对AR（p）模型的残差序列进行Ljung-BoxQ检验，结果见表2.1。从表中可以看出，所计算的Q（p）值均小于相应的检验临界值，不能拒绝残差序列是独立的原假设（在5%的显著性水平下），从而可以认为AR（p）已经消除原时间序列自相关的影响，残差序列已不包含线性相关成分。应用AR模型分别对降雨径流时间序列进行拟和，残差序列的Ljung-BoxQ统计量在5%显著水平条件下均显示很好克服了线性相关关系，可以对残差序列进行BDS非线性检验。采用Kanzler Ludwing算法计算BDS统计量。在BDS统计量计算过程中，取嵌入维数m=2～7，阈值r分别为拟合模型残差序列的标准差σ的0.25、0.50、0.75、1、1.25、1.5倍，对残差序列进行BDS检验，计算结果见表2.2。

由表2.2的BDS统计结果来看，BDS统计量均大于5%显著性水平的标准正态分布的临界值1.96，均可以在较高的显著水平拒绝AR拟合模型的残差是iid的零假设，从而可以认为貌似随机的降雨径流数据不是完全随机的，其中存在着非线性相关结构。因此，为了正确地描述径流序列的变动，必须使用非线性模型，这就为进一步建立径流量变化的动态模型提供了一定的依据。另一方面也说明径流系统演化过程中所出现的波动并非仅简单地由外界随机因素的扰动所造成，其本身还存在着复杂的非线性相关关系，正是因为存在着这种非线性，才使径流系统的演化表现出混沌特性。

2.3.3 混沌特性的定量识别

2.3.3.1 延迟时间τ和嵌入维数m的选取

考虑到其他几个特征量的计算均是建立在相空间重构基础上的，因此，我们首先要先确定相空间重构的2个关键参数：参数延迟时间τ和嵌入维数m。这里，我们采用改进的C-C法同时求取降雨径流时间序列的τ和m，取达到局部到最小作为最佳时延τ，的第一次达到局部峰值作为最佳嵌入窗τ_w，计算结果如图2.2所示 （具体计算原理见3.2.2小节）。

从图2.2上看，观台水文站径流序列的曲线在t=3时达第一次过零点，同时达到局部最小，故取最佳时延 τ=4，而相应的最佳嵌入窗选取在的第一次达到局部峰值处为t_w=12，结合前人的研究成果，即认为径流过程是一个低维混沌过程，选择嵌入维数为5。三个雨量站降雨量序列的曲线基本上都在t=2时第一次过零点，同时达到局部最小，故取最佳时延τ=2，而相应的最佳嵌入窗选取在的第一次达到局部峰值处为t_w=12，这与观台水文站流量序列的时间窗口相一致。根据时延窗口与嵌入维数的关系，三个降雨序列的嵌入维数均选择为7。

在确定了延迟时间的基础上，我们首先应用Cao方法分析降雨时间序列是否具有非线性混沌特性。选取改进C-C法计算的径流最佳延迟时间τ=4、降雨序列的最佳延迟时间τ=2，根据2.2节所述方法计算E₁（m）与E₂（m），得到月降雨径流量时序数据的嵌入维数m与函数E₁（m）、E₂（m）曲线图，如图2.3所示。由图2.3可以看到，E₂（m）、E₁（m）均随着嵌入维数m的增大而增大，并最终趋于1，因而排除了月降雨径流序列为随机序列的可能性。特别的，从E₁（m）与嵌入维数m的关系图中我们还可以得到稳定的嵌入维数，结果与用C-C法求得的嵌入维数基本一致。

2.3.3.2 关联维数计算

接下来我们采用G-P法计算关联维数。取由改进C-C算得的延迟时间，对于径流序列取τ=4，降雨系列取τ=2，固定时间延迟，对嵌入维数从1～15进行重构相空间，按照式（2.23）分别计算出给定一系列标度r时的关联积分函数的值C（r，m），绘制lnC（r，m）-lnr关系曲线。如果系统具有混沌特性，则随m的增大而D₂（m）逐渐趋于一个稳定值，此时的m即为最佳嵌入维数，它所对应的D₂（m）为饱和关联维数。然后选取不同m值及其对应曲线上无标度区间内的所有点，通过最小二乘法计算关联维值和Kolmogorov熵的稳定估计。降雨径流时间序列lnC（r，m）-lnr关系曲线如图2.4所示。

图2.4给出降雨径流时间序列序列对应嵌入维数m=3，4，5，…，15时，lnC（r，m）-lnr之间的变化关系。从图2.4可以看出，C（r，m）与r之间的标度关系在lnC（r，m）-lnr曲线图中表现为一段直线，而且随着嵌入维数的增大，其图形的形状和走向也基本一致。其中，观台水文站径流序列当嵌入维数m=5以后，线性部分逐渐趋于平行，这说明关联维已趋于饱和，此时，饱和关联维数趋近于2.66；而匡门口、白土以及观台雨量站的降雨序列在嵌入维数m=7以后，线性部分逐渐趋于平行，关联维趋于饱和，饱和关联维数分别为4.36、4.16和4.19。

图2.3 降雨径流时间序列E₁（m）-m与E₂（m）-m关系曲线

2.3.3.3 Kolmogorov熵计算

通过拟合降雨径流时间序列lnC（r，m）-lnr关系曲线上无标度区间内的所有点，得到Kolmogorov熵分别为0.22、0.42，0.53和0.53。由此可说明，该降雨径流量的变化存在着低维的混沌吸引子，月径流量变化的时间序列存在无标度区间，吸引子关联维也在这一区间内趋于收敛和稳定，说明了降雨径流系列的变化并不是毫无规律的，而是周期性变化和非周期性变化的叠加，具有复杂的内在动力学机制。

2.3.3.4 Lyapunov指数计算

最后利用2.2节中给出的方法计算出重构相空间中吸引子的最大Lyapunov指数λ₁，结果见表2.3。由表2.3可以看出，计算出来的最大Lyapunov指数λ₁>0，再一次印证了径流系统具有混沌特性。

上述分析方法应用结果表明，观台站径流变化存在着内在动力学机制，是由周期性和非周期性影响因子共同作用的结果，具有明显的混沌特性。径流系统演化过程中所出现的波动并非仅简单地由外界随机因素的扰动所造成，其本身还存在着复杂的非线性相关关系，正是因为存在着这种非线性，才使径流系统的演化表现出混沌特性。因此，为了正确地描述径流序列的变动，必须使用非线性模型，这就为进一步建立径流量变化的动态模型提供了一定的依据。

图2.4 降雨径流时间序列lnC（r，m）-lnr关系曲线