众所周知, 边坡系统是一个复杂的非线性动力系统, 其演化具有从混沌到有序的特征(黄润秋和许强, 1996; Zhou, 2000; 於崇文, 2002).边坡演化混沌特征分析的首要问题是混沌特征量指标的计算与判别.根据混沌动力学理论, 系统的关联维数、Lyapunov指数、K熵是其混沌特征的主要依据(刘秉正和彭建华, 2004).国内外的学者经过探索给出了混沌判别的基本方法(Qin et al., 2001; Zhou et al., 2001; 吕金虎等, 2002; 刘秉正和彭建华, 2004) : (1) 如果边坡系统吸引子的关联维数为一个分数, 那么就可以判断边坡含有混沌成分; (2) 最大的Lyapunov指数大于零, 就可以肯定混沌的存在; (3) 边坡系统的Kolmogorov熵收敛且为一正数也能表明其混沌特征.但是, 也有学者研究发现, 不能以上述方法中的一种或几种作为诊断或判定边坡是否具有混沌特性的依据, 如对Duffing方程研究发现(黄胜伟, 2002).若Lyapunov指数大于零而关联维数为零, 或关联维数大于零而Lyapunov指数小于零, 系统不是混沌而是周期运动.故更为系统深入地研究混沌判别方法是必须的.对于边坡系统而言, 其系统监测数据少, 含有噪声, 因而, 引入数据量要求较低、抗噪能力较好的混沌特征量计算方法是有益的.

在边坡演化方程未知的情况下, 现今一般采用的Takens (1981)时间延迟重构相空间方法, 因其隐含了测量数据具有无限精度的假设, 对边坡混沌研究并不适用; 并且任意选择延迟时间和独立估计延迟时间与嵌入维数的方法可能造成特征量计算效率低或者相点稀疏, 因而, 适应边坡系统特点的相空间重构方法需要探索.

本文针对边坡监测数据特点, 同时估计延迟时间τ和嵌入维数m进行边坡相空间重构, 采用较新方法计算边坡关联维数、最大Lyapunov指数、Kolmogorov熵, 并引入近似熵ApEn及系统复杂度指标, 从多个角度揭示边坡演化的混沌特性, 并通过实例分析建立边坡系统的多元混沌判据.

1.   边坡演化的非线性时间序列相空间重构

针对边坡系统监测数据特点, 采用Kim et al. (1999)提出的能够同时估计τ和m、需求数据较少、抗噪能力较好的C-C方法进行边坡相空间重构, 以位移监测数据为对象, 步骤如下:

(1) 设δ₁, δ₂, …, δ_i, …, δ_N为监测所得的边坡位移数据序列, 按时间延迟重构方法, 将此时间序列重构成m维相空间如下:

式(1) 中: k为延迟时间指数(k=1, 2, …), 假如测量边坡位移的时间间隔为τ_s (即测量δ_i与δ_i+1之间的时间间隔), 那么延迟时间τ=kτ_s; X_i为m维相空间里的一个相点; {X_i}则为m维相空间里的一个相型.相点的连线即表征了边坡的演变轨迹.

(2) 定义嵌入时间序列的关联积分:

式(2) 中: 当α < 0, Θ (α) =0;当α≥0, Θ (α) =1;N为数据集的大小; k是时间延迟指数; M=N-(m-1).k为m维空间的相点数; ‖*‖为Euclidean距离.

(3) 引入函数:

为了研究非线性依赖关系并且消除与时间有关的伪关联性, 将位移时间序列{δ_i} (i=1, …, N) 分成k个不相交的子序列, 这里k是延迟时间指数, 则:

当N→∞时, 有:

(4) 固定m和k, 选择n个代表性的r_j, 定义:

(5) 定义下列平均量:

(7)~(9) 式中: a为数值实验m上限, b为r_j个数, 一般, r_j=0.5*j*σ, σ为数据集的标准差.

(6) 以S(k)第一个零点, 或ΔS (k) 第一个极小值中寻找数据相互独立的局部最优时间, 此最优时间对应延迟时间τ=kτ_s, τ_s为单位时间; 寻找S_cor (k) 最小值点, 给出延迟时间窗口τ_w=kτ_s, 从而确定延迟时间τ和嵌入维数m.

2.   边坡演化混沌特征判别指标与计算方法

2.1   边坡关联维数D₂

以上为C-C法重构边坡位移时间序列相空间, 关联积分如式(2), 因为固定N、k, 式(2) 简写为C(m, r), 关联维数的定义为:

通常的方法是给定一些m的值, 画出lnC (m, r) -lnr曲线, 通过考察曲线中直线段的斜率确定关联维数.但这种方法对于边坡这种有限数据集, 计算出的D₂可能比真实值小.本文采取排除时间相关性点对的办法(Theiler, 1986; Provenzale et al., 1994) 对式(2) 进行修正, 即假定相点在经过t_min=n_minτ_s后不再产生时间相关性.于是, 式(2) 修正为:

式(11) 中: n_min为产生时间相关性的点的最小个数; 其他变量的含义与式(2) 相同.在确定t_min (即确定n_min) 时, 我们采用空间-时间分离图的方法(Grassberger and Procaccia, 1983; 成秋明, 2006) 来确定, 基本思想是: 在时间相关性存在的条件下, 作分离时间Δt和空间距离r为自变量、点对的数目为函数的等值线图来识别.

2.2   边坡最大Lyapunov指数

Wolf et al. (1985)提出的直接基于相轨线的演变估计最大Lyapunov指数算法是现今较通用的方法, 但此法要求较大数据量, 抗噪能力较差, 对边坡系统应用易发生不收敛.近年来, Kautz and Schreiber (1997)提出的方法在其他工程领域的应用取得了较好效果, 本文以此法为基础, 并进行适当改进, 以向量的Euclidean范数作为向量间的距离, 在邻近点集中排除时间相关性强的点, 其主要计算步骤如下:

(1) 按(1) 式重构边坡位移时间序列相空间;

(2) 在X_i的ε邻域内寻找所有的邻近点X_j, 同时排除时间相关性点, 即:

式(12)中: Dist(X_j, X_i; 0) 为向量X_j与X_i之间起始时刻的距离; ‖*‖为Euclidean距离; p为时间序列的平均周期, 可按(11) 式中时间分离图的方法确定.满足上面条件的邻近点构成的集合记为u_j;

(3) 经过t*τ_s后, X_i演变到X_i+1, 这时X_i所有的邻近点X_j演变到X_j+1, 则之间的距离为:

按最大Lyapunov指数的定义, 有下式成立:

两边取对数, 整理可得:

λ₁为与曲线的斜率.为减小误差, 对所有的邻近点以及对所有由X_i点出发的轨道进行平均, 即计算:

式(16) 中: M为相点总数; |u_j|表示集合u_j中的元素个数;

(4) 作曲线, 其斜率就是λ_i.

2.3   Kolmogorov熵

以Kolmogorov熵的定义直接计算边坡系统的Kolmogorov熵是较为困难的, 本文借助2阶Renyi熵K₂对K熵进行估计, 并以关联积分近似代替概率密度函数, 定义为:

由定义: 0≤K₂≤K; 对于随机系统, K₂→∞; 对于混沌系统, K₂≠0.对于边坡位移时间序列, 选取合适的延迟时间τ, 重构m维相空间, 经推导得:

2.4   近似熵

应用K熵或K₂熵研究系统的混沌特性, 较难单独反映规则运动成分的特征, 并且K熵研究中, 要求一定的数据长度N (N一般大于100), 才能达到较好效果, 这对边坡工程是有一定难度的.对这一问题, 本文引入近似熵(ApEn) (Pincus, 1995) 混沌指标, 计算方法如下:

对按(1) 式重构的边坡位移时间序列相空间, 定义矢量X_i与X_j之间的距离d|X_i, X_j|为它们各分量之间最大距离, 即:

类似关联积分, 定义:

再令:

Φ^m (r) 便表示序列{X_i}平均自相关程度.则定义近似熵:

具体计算得:

式(24) 中“ ? ? ”表示对不同i求平均.对于关联维、Lyapunov指数等, 如果N较小, 可能无法收敛或者有较大误差.但ApEn是人为定义的, 对数据量要求较宽松, 通常N=100左右仍可计算; 并且, ApEn具有较好的鲁棒性, 对时间序列中存在的少数异常数据不敏感, 这对于边坡实测数据很有利.对于实际的随机信号(白噪声), ApEn不会像K熵趋于无穷大, 而是趋于有限值, 所以更为合理.

2.5   系统复杂度

系统复杂度是定量化描述系统复杂程度的混沌特征量, 并具有静、动态分布特征, 易于计算, 能够弥补分维、Lyapunov指数、K熵混沌判断指标存在的不足(刘秉正和彭建华, 2004), 笔者将其引入边坡混沌特征研究中.复杂度计算包括统计复杂度和算法复杂度两类概念, 由于无法确定边坡系统的演化方程, 故在边坡系统混沌特征研究中, 算法复杂度具有较大意义.利用符号序列研究系统复杂度, 计算的是L-Z复杂度(Lopez-Ruiz et al., 1995), 记为C_LZ, 具体计算过程如下:

(1) 将要分析的时间序列离散化表示为一数列或字符串(x₁, x₂, …, x_n), 设其平均值为x.对其粗粒化: 用二进制的两个数字(1, 0) 组成的符号序列代替原数列, 即原数列中大于x的数取为1, 小于x的数取为0.

(2) 设字符串: S=S₁S₂…S_r; Q是在S后新加的字符串; SQπ表示在字符串SQ中删去最后一个字符所得到的字符串; ν (SQπ) 表示SQπ中所有不同子字符串的集合.

(3) 规定: ①对于原字符串S, 如果后加的字符串Q包含在S中(Q与S中某一子字符串S相同), 则Q称为是S的复制, SQ被认为是与S属于同一段; ②如果Q不包含在S中, 则Q认为是插入.每遇插入, 便认为分段数增加1, 分段与分段之间用分隔号“.”分隔开; ③字符串的未归一化的复杂度按分段数计, 即复制时复杂度不变, 而每次插入使复杂度加1;④每一字符串的第一字符总是算作新字符而且单独占一分段.

(4) 按如上规定, 从字符串S的第一个字符开始计算整个字符串的分段数, 也就是它的未归一化的复杂度数c (N).

(5) c (N) 的归一化:

(6) 计算平均复杂度: 将时间序列分为N/L个长度为L的窗口, 对窗口内的数据按上述方法计算复杂度C_LZ (L).对这N/L个窗口分别计算其复杂度, 即整个时间序列平均复杂度.

3.   基于时间序列的边坡系统多元混沌判别

选取某高速公路工程4个重点监控的高边坡实际位移监测资料, 研究边坡的混沌特征, 建立边坡系统的多元混沌判据.

3.1   边坡概况

广东省某高速公路工程4个重点边坡分别为K97+640、K96+170、K17+320及K87+450, 边坡所在的地貌单元均属剥蚀丘陵地貌区, 其中K97+640和K96+170右侧边坡由第四系覆盖层和下伏混合花岗岩组成, 地下水发育; K17+320左侧边坡地层为第四系覆盖层, 基岩为混合花岗岩, 地下水不发育; K87+450边坡表层为第四系坡积亚粘土, 下伏全-弱风化花岗岩, 地下水不发育.4个边坡于2004年9月-10月开始监测, 根据监测报告, K96+170边坡与K87+450边坡分别于2005年6月中下旬发生小面积滑坡, K97+640边坡位移监测亦明显增大, 而同期K17+320边坡位移无明显变化.

3.2   相空间重构

对4个边坡位移监测资料按C-C法进行相空间重构, 取得基本相同结果.图 1为K97+640右侧边坡位移监测结果.

图  1  x, y, z方向分量位移曲线

Fig.  1.  Displacement curves in x, y, z directions

下载: 全尺寸图片 幻灯片

如图 2, ΔS较S首先出现极值, 此时k≈7;当k=15时, S_∞出现最小值, 则延迟时间窗口为15;由7 (m-1) τ_s=15τ_s, 取最小整数, 得m=3.则以延迟时间τ=7τ_s、m=3对地表整体位移时间序列, 重构相空间.

图  2  相空间重构分析曲线

Fig.  2.  Cuurves show the reconstruction of phase space

下载: 全尺寸图片 幻灯片

3.3   混沌特征量的求取与混沌特征判别

分别按照本文混沌特征量的求取方法计算关联维数、最大Lyapunov指数、K熵、近似熵和L-Z复杂度, 计算结果见表 1.

表  1  不同边坡混沌特征量计算值

Table  Supplementary Table   Calculation data show the chaotic features of different slopes

下载: 导出CSV 

| 显示表格

非线性动力学研究表明(刘秉正和彭建华, 2004) : 非线性系统内部特征量的取值情况可以用来表征系统所处的运动状态(表 2).因此, 我们通过综合考虑关联维数、最大Lyapunov指数、K熵、近似熵和L-Z复杂度的取值情况来进行边坡系统的混沌判别.参照表 2可知, 除K17+320边坡外, 计算所得的边坡特征量值全部在期盼标准范围内.表明多数边坡演化具有混沌特点, 并且依据监测资料显示, 除K17+320左侧边坡外, 其他边坡均处于临近破坏状态, 说明在临滑失稳阶段边坡系统混沌特征显著.

表  2  不同状态运动的特征量取值比较

Table  Supplementary Table   Comparison of different characteristic values from different motion types of slopes

下载: 导出CSV 

| 显示表格

4.   结论

(1) 研究表明: 对于边坡系统, 最优相空间延迟时间τ介于7~15τ_s之间, 而嵌入维数约为3;边坡系统关联维数为非整数, 最大Lyapunov指数、熵值均大于零, 复杂度位于(0, 1) 区间偏小值, 说明了影响边坡演变的因素中除了有通常认为的确定性成分和随机成分外, 还存在混沌成分, 边坡演化具有混沌特性.

(2) 近似熵ApEn以及系统复杂度混沌特征量的引入, 是边坡混沌特征判断的新方法, 与关联维数, 最大Lyapunov指数和K熵一起, 能够建立起较为完善的边坡多元混沌判断理论, 并且为其他非线性系统混沌研究提供借鉴.

(3) 基于监测的位移时间序列数据进行边坡的混沌特征分析与判别是一项较有意义的工作, 可为边坡演化规律研究提供重要的基础和研究新途径.