MF-DFA在大坝安全监测序列分析和整体性态识别中的应用

胡 江，苏怀智，马福恒，李子阳

(1.南京水利科学研究院水文水资源与水利工程科学国家重点实验室， 江苏 南京 210029； 2.河海大学水利水电学院，江苏 南京 210098)

变形和应力是水压力、温度、渗流、材料和地质特性等多因素的综合反映，其中变形观测直观可靠，国内外普遍将其作为大坝监测最主要效应量[1-3]。大坝与基岩有机组合成为一个复杂的动力系统，多因素的共同作用使大坝变形观测时间序列(位移和挠度倾斜等)具有显著的时空非线性特征。依据实测变形时间序列，借助先进的数学模型，分析其中蕴含的复杂物理机理，是大坝安全监控的重要手段。

一般采用统计模型、灰色系统理论和模糊数学模型等分析大坝原型观测时间序列，但由于时间序列的复杂非线性动力特性，很难用参数形式描述其中的数字特征。苏怀智等[4-5]应用分形学分析法挖掘观测时间序列的内在规律性，采用Hurst指数H和分形维数D分析大坝结构的工作性态和发现病害隐患，然而采用单一的Hurst指数从整体上识别整个时间序列的结构特征过于笼统，缺乏对时间序列中的小波动的细节分析。

多重分形消除趋势波动分析 (multifractal detrended fluctuation analysis，MF-DFA) 方法是在单分形的基础上提出的，是定量描述复杂序列非线性演化过程及其形成自相似结构特征的有效手段，最先被广泛应用于自然现象、生命科学和医学等多个领域[6-7]。通过MF-DFA，不仅可以根据广义Hurst指数检测序列的长程相关性和确定标度不变性，还可以根据多维分形谱区别不同时间序列演变规律的差异。近年来，MF-DFA方法的潜力被逐步挖掘，Telesca等[8]将MF-DFA方法应用到大地电磁场的研究中，由不同场地的大地电磁数据序列的多重分形分析，辨识了环境因素对序列波动的影响，由此解释序列的动态变化。本文运用MF-DFA方法研究大坝监测时间序列的长程相关性和多重分形特征，通过多测点的异同分析，探索大坝系统整体监控模型，以此对大坝健康进行诊断。

1 MF-DFA方法

消除趋势波动分析(detrended fluctuation analysis，DFA)方法最初是由Peng等[9]在研究DNA序列机理时提出的一种分析幂律函数关系的方法，该方法可有效滤去各阶趋势成分，适合分析非平稳时间序列的长程相关性。在DFA方法的基础上，Kantelhardt等[10]提出了MF-DFA方法，MF-DFA方法提高了非平稳序列波动奇异性的多重分形结构度量的准确性，很快被被运用于时间序列蕴含的多重分形性质的研究。

假定长度为N的原位监测序列为{x1,x2,…,xN}，MF-DFA方法的计算过程如下：

a. 计算序列的累积离差：

(1)

式中：xave为序列的均值。

b. 把序列Yi分割成长度为s的Ns=int(N/s)个互不重叠的等长区间(图1)，由于长度N通常不是尺度s的整数倍，序列尾部可能剩余部分。为了不丢弃这一序列部分，从序列尾部重复这一分割过程。因此，共分割得到2Ns个区间。

图1 消除趋势波动分析中去趋势程序示意图

c. 采用最小二乘拟合法计算序列2Ns个区间中每个区间的局部趋势。定义时间尺度s对应的消除趋势时间序列Ys,i，与原始时间序列和拟合序列不同，消除趋势时间序列Ys,i表示为相对局部平均的波动:

Ys,i=Yi-yv,i

(2)

式中：yv,i是第v个区间的拟合多项式。图1表示了s分别为100和200时采用二次多项式拟合(即DFA2)的情形。拟合过程还可采用线性、三次式或更高阶多项式，相应地，被记为DFA1、DFA3和高阶DFA。确定2Ns个区间中每个区间v的均方差，当v=1,2,…,Ns时,有

(3)

当v=Ns+1,Ns+2,…,2Ns时

(4)

d. 计算所有区间的平均值得到整个序列的q阶波动函数:

(5)

式中：下标q可以取0之外的任何的实值，当q=0时，由洛必达法则有

(6)

e. 在smin≈6和smax≈N/4的范围内变换尺度s的值，确定消除趋势波动函数Fq(s)和尺度s大小间的标度关系，即

Fq(s)-sh(q)

(7)

一般地，指数h(q)随q不同而变化，对于非平稳序列，h(2)等于Hurst指数，因此，将h(q)定义为广义Hurst指数。对于单分形序列，h(q)不随q不同而变化。这是因为对于所有分割区间v，均方差F2(v,s)相等，且由式(4)表示的平均过程计算得到的波动函数也相等。只有当大、小波动明显不同时，h(q)才会随q显著变化。当q为正值时，由大的均方差F2(v,s)(即与相应的拟合值偏差大)的区间v决定均值Fq(s)；相反地，当q为负值时，由小的均方差F2(v,s)的区间v决定均值Fq(s)。换句话说，q为正值，h(q)刻画了大波动区间的分形行为；q为负值，h(q)则刻画了小波动区间的分形行为[10]。

图2 大坝垂线的布置方案(单位：m)

f. 对于多重分形序列，MF-DFA的最后结果是一系列的分形指数h(q)，是q的减函数。根据标准多重分形分析，由MF-DFA方法计算得到的h(q)与标度指数τ(q)的解析式[10]为

τ(q)=qh(q)-1

(8)

式(8)表明，MF-DFA中由式(7)定义的广义Hurst指数h(q)与传统的多重分形标度指数τ(q)直接相关。另一种描述多重分形序列的方法是奇异谱f(α)，通过Legendre转换得到f(α)和τ(q)间的关系

α=τ′(q)f(α)=qα-τ(q)

(9)

式中：α为奇异强度，或称为Hölder指数；而奇异谱f(α)则表示了由α描述的序列子集的维数。由式(8)可得

α=h(q)+qh′(q)f(α)=q[α-h(q)]+1

(10)

奇异谱f(α)的形状和延伸长度包含了分析序列的分布特征的显著信息。通常奇异谱f(α)有一个下凹的曲率，分形谱宽度Δα(Δα=αmax-αmin)增加，意味着分布的奇异性增大。

2 应用实例

2.1 工程背景和数据来源

某工程属 I 等枢纽工程，主要由碾压混凝土重力坝、坝顶开敞式溢洪道、泄水底孔、左岸输水建筑物等组成。碾压混凝土重力坝最大坝高 113.0 m，坝顶长308.5 m，坝顶高程为179.0 m，水库正常蓄水位173.0 m，属不完全调节水库。主体工程于2000 年12 月18 日下闸蓄水。

为确保大坝及地下厂房的安全运行，在主要建筑物表面(或内部)布置了变形、渗流、温度和应力应变等较为全面的监测项目。大坝变形监测项目主要有正垂线(PL)、倒垂线(IP)、引张线和视准线，垂直位移监测采用水准测量。大坝正、倒垂线用于监测坝顶及坝体内部的水平顺河向、垂直水流向位移量，布置方案如图2所示。为提高可比性，本文选取悬挂点均位于坝顶179.0 m高程的4号坝段(溢流坝段)、2号坝段(挡水坝段)和5号坝段(挡水坝段)的正垂线PL7、PL5和PL6的自动化监测时间序列进行多重分形分析。为方便分析，顺河向水平位移以向下游为正，向上游为负。可用的时间序列自2002年10月开始，选取2003年1月1日至2008年12月31日间的3个测点的顺河向位移序列进行分析，实测的环境量和位移时间序列如图3所示。

图3 环境量和正垂线顺河向位移实测时间序列

2.2 统计描述与正态检验

从图3可以看出，变形呈一定程度的增加趋势，但并非单调增加，而是有涨有落，且变形波动过程呈一定的周期性。计算PL5、 PL6和PL7测点径向位移序列的标准差、峰度和偏度等基本统计量，结果见表1，可看出，各序列的偏度都不为 0，且峰度都小于0，因此均属于非正态分布。PL7和PL5测点位移序列分布的偏度均大于0，这表示其数据分布形态与正态分布相比为右偏；PL6测点位移序列分布的偏度小于0，其数据分布形态与正态分布相比为左偏。同时，所有测点序列的峰度均小于0，这表示测点序列的整体分布与正态分布相比较为平缓，为宽顶峰。由此可知，测点序列为呈右偏或左偏的非正态分布。

表1 正垂线径向位移时间序列的基本统计量

2.3 多重分形分析结果

根据式(1)～(4)计算波动函数Fq(s)，参照文献[10]，这里波动函数的阶数取值范围为-6≤q≤6；对于时间尺度s，取值需适中，取值太小没有实际意义，过大，则分割的区间数太少，结果将与单分形结果相似，一般地，时间尺度s的范围为5≤s≤N/4，根据大坝变形监测资料的周期性波动趋势，这里选取s={6,12,24,48,96,192,384}。图4为采用MF-DFA1计算得到的q=±6、±4、±2、0时温度和库水位及PL5、PL6和PL7测点径向位移的Fq(s)-s对数图。整体来说，小的时间尺度对应的q阶波动函数的差异较大的时间尺度显著，小尺度能够识别局部时间内的大小波动；相反，大尺度在一定程度上均化了大小波动的幅度。

图4 波动函数和时间尺度的对数关系

图5为广义Hurst指数h(q)随q的变化情况，h(q)-q曲线能反映序列是否存在多重分形行为。当h(q)不随q变化时，表明序列是单分形；反之，序列为多重分形。由图5可以看出，由环境量序列和各测点位移序列计算得到的广义指数h(q)均随q增加而单调减少，这意味着环境量序列和各测点位移序列都存在较显著的多重分形特征。当q=2时，计算得到的温度和库水位序列h(2)的值分别为1.31和1.46，各测点位移序列h(2)的值分别为1.38、1.43和1.5。

图5说明位移序列存在明显的分形特征，因此有必要进行多重分形谱分析。由式(8)和式(9)分别计算得到多重分形标度指数τ(q)和分形谱宽度Δα、奇异谱高度Δf(α)，得到的τ(q)-q和f(α)-α曲线分别如图6和图7所示。其中，温度和库水位序列的分形谱宽度Δα分别为0.612和1.227，奇异谱高度Δf(α)分别为0.715和1.122；PL5、PL6和PL7测点序列的分形谱宽度Δα分别为0.614、0.728和0.679，奇异谱高度Δf(α)分别为0.72、0.79和0.69。

图5 广义Hurst指数和q的关系

图6 位移日测值序列质量指数谱曲线

图7 位移日测值序列多重分形奇异谱曲线

2.4 结果讨论和大坝的整体性态识别

从图4(a)(b)可以看出，温度序列的Fq(s)-s曲线在s=96处扭曲，对应时间尺度大约为一个季度；水位序列在s=192处扭曲，对应时间尺度大约为半年周期。从图4(c)(d)(e)可以看出，PL5测点序列的Fq(s)-s曲线都在s=96处扭曲；PL6测点序列有所不同，在s=24处接近月周期附近扭曲；PL7测点曲线扭曲发生在s=384接近年周期处。同时，由图3(b)可以看出，各测点位移序列大小波动存在一定的差异，Fq(s)-s曲线较好地反映了测点序列局部时间内和较大尺度范围内的日测值的波动情况。

对测点序列进行统计分析，其中水压因子取3次方，温度取半年周期，时效取对数形式[1]，统计模型分析得到2007年典型测点年变幅的结果(表2)，计算得到PL5、PL6和PL7测点序列的温度分量约占总位移年变幅的61.7%、63.5%和63.0%，温度变化对坝体上下游方向水平位移变化影响最大，其次为库水位变化；时效分量相对较小，表明坝体变形趋于稳定。

表2 统计模型分析结果 mm

测点序列的h(2)值大于1且接近于1.5，这说明测点序列的波动相关性和布朗运动有一定相似性，换句话说，监测序列存在较大的不确定性。同时，各序列h(2)值接近，由文献[4]可知，整体上大坝结构处于弹性状态。由h(q)-q曲线q>0和q<0两段可以区分各测点的大小波动的差异。图5(a)表明温度序列的大波动标度无显著变化，库水位序列的大波动标度相对较大；温度和库水位序列尤其是库水位序列的小波动标度单调递减较显著。图5(b)表明PL7测点序列波动较显著，其次为PL6测点系列，PL5测点序列波动较小；各测点序列的大波动标度单调递减较明显，这是由于蓄水和运行初期，环境量和坝体自身引起的，各测点递减的趋势基本相同，表明整体变化趋势相同；PL6和PL7测点序列的小波动标度单调递减趋势也基本相同；PL5测点序列h(q)序列小波动的标度较其他两测点小，表明PL5测点序列的波动较小。

f(α)-α曲线可以反映位移序列在分形结构上不均匀分布的性质。由图7(a)可见温度和库水位序列尤其是库水位序列的奇异谱呈明显的右钩状，小波动的影响在库水位序列中占显著优势。由图7(b)可见各测点序列尤其是PL5和PL7的奇异谱基本对称；而PL6测点序列奇异谱则呈不明显的右钩状，即对于PL6测点序列，小波动的影响在日测值序列中稍占优势。如前所述，分形谱宽度Δα的值可以反映大坝性态演变复杂度的空间差异，各测点分形谱宽度Δα相差不大，相对而言，PL6测点序列的奇异谱谱宽Δα最大，PL5次之，PL7最小。从图2可以看出，PL6测点高程较低，水位波动相对坝段坝高在3个测点中最大；PL7测点位于河床，水位波动相对坝段坝高在3个测点中最小。初步推断，各坝段不同的运行机理是导致坝顶变形测点序列多重分形特征差异的根本原因。

由以上分析和讨论可以看出，整体上，大坝结构变形处于正常状态，坝体处于弹性状态；对比来看，岸坡5号坝段的PL6测点序列表现出较为复杂的性态，河床4号坝段的PL7测点序列较为稳定。综上，统计模型方法能够得到监测效应量的主要影响因子；MF-DFA方法可以识别监测效应量序列不同时间尺度上的分形特征和大坝整体性态时空差异和演变规律。

3 结 论

a. 各垂线测值序列分别在月、季度和年周期时间尺度上，波动函数存在一个转折点，反映了各测点测值时间序列在不同时间尺度上的波动特征，说明各坝段运行机理具有一定的差异。环境量的波动函数表明温度和水位存在较为显著的季度和年周期波动，同时广义Hurst指数说明水位大小波动较为显著；通过统计模型分析表明，温度分量和水位分量占变形量年变幅的大部分。综合统计模型分析和MF-DFA方法分析结果，整体上，温度和库水位的大波动决定了坝体变形幅度，水位频繁涨落控制各坝段尤其是岸坡坝段变形的小波动。

b. 广义Hurst指数h(q)均随q的增加而减小，大坝变形测值序列呈显著的多重分形特征。由Hurst指数大于1且接近于1.5说明运行初期大坝整体上处于正常线弹性状态。测点序列波动具有较大的不确定性，这与大坝运行初期受到较为复杂的环境、水文和运行管理等多重因素影响有关。

c. 综合环境量和测点序列的奇异谱形状和宽度，可以刻画大坝不同坝段和局部运行性态的动力学特征，其物理意义明确，为此通过大坝多个测点的奇异谱宽度分析可反映大坝长期性态演变的时空差异。进一步挖掘造成不同测点复杂程度差异的可能原因，如坝段的工作条件和材料的衰变等，可以指导大坝的运行管理和维护。

参考文献：

[1] 顾冲时，吴中如.大坝与坝基安全监控理论和方法及其应用[M].南京:河海大学出版社，2006.

[2] 远近，郑东健，李波，等.模糊可拓模型在坝体变形评价中的应用[J].水利水电科技进展，2009，29(1)：23-25.(YUAN Jin,ZHENG Dongjian,LI Bo,et al.The application of fuzzy extension model in dam deformation evaluation[J].Advances in Science and Technology of Water Resources,2009,29(1):23-25.(in Chinese)

[3] 胡江，苏怀智.再生缝对拱坝变形的影响分析与预测[J].水电能源科学，2008，26(5)：96-100.(HU Jiang,SU Huaizhi.Influence of regenerated crack to arch dam deformation and displacement forecast[J].Water Resources and Power,2008,26(5):96-100.(in Chinese)

[4] SU Huaizhi,HU Jiang,WU Zhongru.A study of safety evaluation and early-warning method for dam global behavior[J].Structural health Monitoring:An International Journal,2012,11(3):269-279.

[5] 赖道平，吴中如，周红.分形学在大坝安全监测资料分析中的应用[J].水利学报，2004(1)：100-104.(LAI Daoping,WU Zhongru,ZHOU Hong.Application of fractal theory to analyze dam safety monitoring data[J].Journal of Hydraulic Engineering,2004(1):100-104.(in Chinese)

[6] TELESCA L,PIERINI J O,SCIAN B.Investigating the temporal variation of the scaling behavior in rainfall data measured in central Argentina by means of the detrended fluctuation analysis[J].Physica A,2012,391(4):1553-1562.

[7] 陈莹，许有鹏，尹义星，等.长江干流日径流序列的多重分形特征[J].地理研究，2008，27(4)：819-828.(CHEN Ying,XU Youpeng,YIN Yixing,et al.Multifractal characteristics of daily discharge series in the Yangtze River[J].Geographical Research,2008,27(4):819-828.(in Chinese)

[8] TELESCA L,LOVALLO M.,HSU H.L,et al.Analysis of site effects in magnetotelluric data by using the multifractaldetrended fluctuation analysis[J].Journal of Asian Earth Sciences,2012,54-55:72-77.

[9] PENG C K,BULDYREV S V,HAVLIN S S,et al.Mosaic organization of DNA nucleotides[J].Physical Review E,1994,49(2):1685-1689.

[10] KANTELHARDT J W,ZSCHIEGNER S A,KOSCIELNY-BUNDE E,et al.Multifractal detrended fluctuation analysis of nonstationary time series[J].Physica A,2002,316(1):87-114.