基于吸引子特征量的非线性判定法在大坝位移监测数据时间序列中的应用

Application of Nonlinear Decision Method Based on Attractor Characteristic Quantity in Time Series of Dam Displacement Monitoring Data

摘要: 由于大坝位移监测数据存在一定的噪声，为了解决高噪声水平下非线性判定的准确性问题，引入基于吸引子特征量的非线性判定方法。本文对比分析了Lorenz方程和Henon映射2种非线性混沌时间序列在不同噪声干扰情况下，时间反演不可逆量及吸引子特征量分别作为特征统计量时替代数据法的检验性能，发现以Hurst指数作为非线性特征量时的替代数据法显示出较强的鲁棒性。将其应用于实际工程中的大坝三个测点的位移监测数据时间序列，结果表明大坝监测数据位移时间序列呈现出非线性变化的特点，为分析大坝结构形态的变化，及时掌握大坝的运行状态奠定基础。

Abstract: Because of the noise of dam displacement monitoring data, a non-linear judgment method based on attractor characteristic is introduced to solve the accuracy of non-linear judgment under high noise level. In this paper, the test performance of two kinds of nonlinear chaotic time series, Lorenz equation and Henon map, is compared and analyzed when the time inversion irreversibility and attractor characteristic are used as feature statistics respectively under different noise dis-turbances. It is considered that the surrogate data method with Hurst exponent as the non-linear characteristic shows strong robustness. It is applied to the time series of displacement monitoring data of three dam observation points in practical engineering. The results show that the dis-placement time series of dam monitoring data presents the characteristics of non-linear change. This lays the foundation for analyzing the change of dam structure and grasping the operation state of the dam in time.

1. 引言

水利行业是关系到国家发展和人民生活的基础性行业，大坝监测数据的分析 [1] 对大坝系统的安全运行，预防溃坝事件的发生有重要的意义。在设计的情况下，大坝通常根据线性弹性模型来考虑，但由于变形坝体中使用的建筑材料的不同，大坝不是理论上的弹性体。与此同时，运行中的大坝是一个复杂的动力系统，在这个系统中，坝体、库水和坝基间的相互作用以及外部环境(如温度、水压和地震)都将影响大坝 [2] ，导致大坝的实际性态与预期的设计不一致，多种因素的影响使得坝体的变形具有很高的非线性特性 [3] 。

随着科学技术的进步，自动化监测仪器在水利工程安全监测中的应用愈加广泛，但自动化监测仪器对周围环境变化的敏感性，致使监测数据往往含有较大噪声，这对位移监测数据时间序列的非线性特性判定造成不便。因此，准确判定高噪声水平下大坝位移监测数据的非线性特性，才能通过分析大坝结构形态的变化来掌握大坝的运行状态，避免事故的产生。

目前，时间序列的非线性特征判定方法主要分为定性和定量两个方面 [4] 。定性判定方法主要根据时间序列在时域或者频域内所呈现出的某些特性或其主要特征粗略分析 [5] ，如功率谱法、Poincare截面法、相图法等。但其在判定中需要根据图形进行人为的判定，常常混入了很多的主观因素，不利于实际应用。

定量判定方法中，其中一类方法是通过计算关联维数、Lyapunov指数、Hurst指数、分形维数及形貌系数等特征量直接识别时间序列数据中的混沌动力学特性；一类是使用Theiler [6] 等人提出的以替代数据(surrogate data)作为检验时间序列中非线性成分的方法 [7] [8] ，通过检验数据中的非线性成分，间接地判断其非线性特性。此方法一经提出，就引起了各领域在非线性判定时的应用 [9] [10] [11] 。考虑到此方法在高度噪声及数据长度的影响下极易受到影响，所以，基于吸引子特征量 [12] 的替代数据法对大坝位移监测数据时间序列进行非线性特征判定是一种更符合实测数据的办法。

2. 吸引子特征量

每一组位移监测数据生成的时间序列，都蕴含着系统原始的基本信息。那么，对每一组位移监测数据时间序列进行相空间重构，得到其相空间的拓扑结构及原系统的吸引子 [13] ，就可以达到以单一变量研究其运动特性的目的。刻画动力系统吸引子某个方面特征的量被称为动力系统吸引子特征量。最常见的吸引子特征量包括：Lyaupunov指数、分形维数、Kolmogorow熵和Hurst指数等。

1) Lyapunov指数

标志着一个系统是否做混沌运动的物理量是Lyapunov指数 [14] 。当此系统中的最大Lyapunov指数大于零时，系统是混沌的；相反当系统的最大Lyapunov指数小于零时，系统是收敛的。Lyapunov指数的定义式为

$\lambda =\lim \frac{1}{n}{{\displaystyle \underset{i=0}{\overset{n-1}{\sum }}\ln \left|\frac{\text{d}F\left(x\right)}{\text{d}x}\right|}}_{x=x_i},\left(i=1,2,\cdots ,n\right)$ (1)

当n极大时， $\frac{\text{d}F\left(x_0\right)}{\text{d}{x}_0}$ 就可以定义为函数 $F\left(x\right)$ 在 $x_0$ 处的一阶导数。实际计算中往往将Lyapunov指数简化为求窗口大小数据拟合的直线的斜率。

2) 分形维数

“分形”的特征就是具有自相似性，是描述混沌运动的几何语言。人们通过定义容量维，Huasdorff维、盒维数、信息维和关联维等分形维数来表征分形，本文主要计算关联维数。

定义m维嵌入空间内的关联积分，关联函数C_XX(ε)在ε→0时，存在一个常数D₂满足以下关系，意味着它刻画了相空间中可能存在的相似性。

$\underset{r\to 0}{\lim }{C}_{XX}\left(\epsilon \right)\propto {\epsilon }^{D_2}$ (2)

则D₂被称为关联维数，此时重构吸引子具有分形特征。

3) Kolmogorow熵

Kolmogorov熵描述的是混沌轨道随时间演化信息的产生率。当Kolmogorov熵大于0，且为常量时，该系统是混沌系统。采用关联积分法计算，

$K_2=-\underset{\tau \to 0}{\lim }\underset{r\to 0}{\lim }\underset{m\to \infty }{\lim }\frac{1}{m\tau }\ln C\left(r,m\right)$ (3)

由于重构相空间之后，τ及m均为常数，则

$K_2=\frac{1}{m\tau }\ln \frac{r^{D_2}}{C\left(r,m\right)},r\to 0$ (4)

画出关系图ln(r)~ln(K₂)，其最佳线性拟合直线在纵轴上的截距即为Kolmogorov熵K₂的稳定估计。

4) Hurst指数

Hurst指数是英国水利专家H. E. Hurst在1951年提出的一种判别时间序列是否对于时间有依赖的参数。任意在时间序列中取一个长度为n的时窗，求取该时窗内的数据均值 ${\bar{y}}_t$ ，计算相应的累积偏差A(t, n)。将同一个时窗n内累积偏差的最大值和最小值之差称为域，记为

$R\left(t,n\right)=\max A\left(t,n\right)-\min A\left(t,n\right)，\left(t_1\le t\le t_n\right)$ (5)

标准差S(t, n)定义为

$S\left(t,n\right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{j}{\sum }}{(y_{t_i}-{\bar{y}}_t)}^2}$ (6)

引入无量纲的比值R/S对R进行重标度，可以得到以下的经验公式

$\left(\frac{R\left(t,n\right)}{S\left(t,n\right)}\right)\propto n^H$ (7)

式中的H即为Hurst指数。实际计算时，对同一个时窗长度n，在不同的初始时刻t的值取均值， ${\log }_{10}\left(E\left(\frac{R\left(t,n\right)}{S\left(t,n\right)}\right)\right)~{\log }_{10}\left(n\right)$ 曲线的斜率即为Hurst指数。

3. 基于吸引子特征量的替代数据法

替代数据检验法是一种鲁棒性强且易操作的非线性数据检验方法。其基本思想是：首先假设待检验的时间序列是一个特定的线性时间序列，以此作为零假设H₀；再根据零假设使用某种生成算法产生一组保留待检验时间序列特定统计特性的替代数据；然后分别计算原始数据和替代数据的某一特征统计量；最后使用Sigma统计检验方法 [15] ，根据原始数据和替代数据检验统计量的显著性差异程度，在一定置信度内接受或拒绝零假设H₀。用替代数据检验法进行时间序列非线性检验的关键之处在于替代数据生成方法的选择与特征统计量的选取。

3.1. 替代数据的生成

目前基于替代数据法的非线性检验中，不同文献均采用不同的替代数据生成方法和检验统计量。本文采用迭代幅度匹配的傅立叶变换(Iterated Amplitude Adjusted Fourier Transform, IAAFT)算法产生替代数据。

IAAFT算法是Schreiber和Schmitz在提出的改进的迭代生成算法，此算法解决了幅值匹配的傅里叶变换(Amplitude Adjusted Fourier Transform, FT)生成替代数据时的功率谱白化问题。IAATA算法能很好的匹配原始数据的傅立叶频谱和概率密度分布，并且被广泛应用于数据的非线性判定 [16] 。

3.2. 基于吸引子特征量的特征统计量

实际工程的位移监测数据，不可避免地混入了噪声，采用常用的特征统计量可能会影响其检验性能。故本文使用4种吸引子特征量作为替代数据法中的特征统计量，与时间反演不可逆量形成对比，检验不同的统计量在Sigma统计检验方法下的检验性能。该方法假设满足零假设的时间序列的特征统计量服从正态分布， $T_o$ 是原始数据统计量的值，替代数据的统计量的均值 $T_s$ 和方差 ${\sigma }_s^2$ ，计算Sigma检验量 [17] 。

$Sigma=\frac{\left|T_0-T_s\right|}{{\sigma }_s^2}$ (8)

在实践中，将0.05显著性水平作为标准。数理统计中正态分布表中，检验统计量依95%概率拒绝零假设时的 $Sigma=1.96$ ，即当 $Sigma>1.96$ 时拒绝零假设，原始序列是非线性的； $Sigma<1.96$ 时接受零假设，原始序列是线性的。

4. 数值模拟

4.1. 数据来源

混沌时间序列是一种典型的非线性时间序列，选用Lorenz方程和Henon映射在不同噪声干扰情况下的非线性混沌时间序列作为分析对象，分析将时间反演不可逆量及吸引子特征量分别作为特征统计量时的替代数据法的检验性能。

1) Lorenz方程

Lorenz吸引子是1963年美国气象学家洛伦兹在研究天气预报问题时得到的第1个表现奇异吸引子的动力学系统，本文对第二个状态变量 $\stackrel{\dot{}}{y}$ 进行研究：

$\{\begin{array}{c}\eta \stackrel{\dot{}}{x}=s\left(y-x\right)\\ \eta \stackrel{\dot{}}{y}=rx-y-xz\\ \eta \stackrel{\dot{}}{z}=xy-bz\end{array}$ (9)

其中，设定系统参数s = 10，r = 28，b = 8/3，η =1。

2) Henon映射

Henon映射是天文学家Henon从研究球状云团以及从洛伦兹吸引子中得到启发而得出的 [5] ，其离散型为：

$\{\begin{array}{l}{x}_{n+1}=b{y}_n+1-a{x}_n^2\\ y_{n+1}=x_n\end{array}$ (10)

其中，设定系统参数a = 1.4，b = 0.3，初始值为(0.1，0.2)。

4.2. 实验结果分析

利用替代数据法对Lorenz方程和Henon映射在不同噪声干扰情况下的非线性混沌时间序列，生成100组替代数据，进行Sigma检验，结果如表1所示。

Table 1. The test results of nonlinear chaotic time series under different noise levels

表1. 不同噪声水平下混沌时间序列非线性检验结果

注：^*表示检验统计量小于1.96，即未能正确检验的非线性时间序列。

由表1可知，随着噪声水平的提高，各种特征量的非线性检验结果都有一定的波动；当噪声水平较高时，时间反演不可逆量已经不能准确判定混沌时间序列的非线性，而除分形维数外的其余三种吸引子特征量对噪声显示出较高的鲁棒性；从检验性能来看，Hurst指数作为特征统计量时的Sigma统计量最大，其图形检验结果见图1，说明原始时间序列Hurst指数与替代数据的统计特性差异最大。

(a) Lorenz方程 (b) Henon映射

Figure 1. Nonlinear testing of chaotic time series with hurst index (The linear representation of the original sequence of the Hurst index)

图1. 用Hurst指数对混沌时间序列进行非线性检验(直线表示原始序列的Hurst指数)

所以，对于噪声干扰强烈的大坝位移监测数据时间序列而言，选择Hurst指数作为特征统计量时检验性能最好。

5. 工程应用实例

5.1. 工程概况

某大坝位于福建省永定县，碾压混凝土重力坝最大坝高113.0 m，坝顶全长308.5 m，坝顶高程179.0 m。该水库正常蓄水位17 3m，校核洪水位177.80 m。自1997年以来，该工程的监测自动化系统已经建成使用 [18] 。为了监测坝顶及坝体内部的水平径向、切向位移量，该监测系统在2#、4#、5#坝段布置了7条正垂线，5条倒垂线。在4#坝段依次布置有三条正垂线PL7、PL4和PL1，其监测点分别为PP4、PP3和PP2。具体监测测点分布见图2。

Figure 2. Layout of deformation monitoring for a dam

图2. 某大坝变形监测布置图

本文以该大坝4^#坝段三个测点自2003年1月1日至2008年12月31日的自动化水平位移数据为例，应用基于吸引子特征量的非线性判定法。位移监测数据时间序列如图3所示。

Figure 3. Time series of displacement monitoring data of PP2, PP3 and PP4 measuring points

图3. PP2、PP3及PP4测点位移监测数据时间序列

5.2. 计算结果

1) 对三个测点的位移监测数据时间序列在0.05显著性水平下，零假设为原始数据由线性相关随机过程经静态、非线性变换产生，利用IAAFT算法产生100组替代数据，其Sigma检验量计算结果如表2所示。

Table 2. The nonlinear test results on time series of displacement monitoring data

表2. 基于位移监测数据时间序列的非线性检验结果

由表2可知，在0.05显著性水平下，利用IAAFT算法产生的替代数据时间序列，其Sigma检验量均大于1.96，拒绝监测时间序列是线性时间序列这一零假设，说明吸引子特征量可以作为特征统计量来判定大坝位移监测数据时间序列的非线性。

2) 采用Hurst指数的替代数据法检验结果如图4所示。

(a) PP2测点 (b) PP3测点 (c) PP4测点

Figure 4. Nonlinear testing of chaotic time series with hurst index (The linear representation of the original sequence of the Hurst index)

图4. 基于Hurst指数的替代数据法检验结果(直线表示原始序列的Hurst指数)

根据图4可知，当Hurst指数作为非线性检验量，替代数据和原始数据存在显著差异，因此大坝监测数据时间序列为非线性序列。

6. 结论

通过实验分析可知，随着噪声水平的提高，各种特征量的非线性检验结果都有一定的波动；当噪声水平较高时，时间反演不可逆量已经不能准确判定混沌时间序列的非线性，而除分形维数外的其余三种吸引子特征量对噪声显示出较高的鲁棒性；从检验性能来看，Hurst指数作为特征统计量时的Sigma统计量最大，说明原始时间序列Hurst指数与替代数据的统计特性差异最大。因此对高噪声水平下的时间序列进行非线性检验时，选用Hurst指数作为特征统计量最佳。

在对大坝监测数据水平位移时间序列进行非线性判定时，使用IAAFT算法生成替代数据，将Hurst指数作为特征统计量，按照Sigma检验原理进行检验。结果表明大坝监测数据位移时间序列呈现出非线性变化的特点，为分析大坝结构形态的变化，及时掌握大坝的运行状态奠定基础。