周兰庭,柳志坤

(河海大学水利水电学院，江苏 南京 210098)

大坝变形监测数据是大坝安全运行的直观体现，其主要受温度、水压、时效三大因素影响[1-2]。由于大坝整体处于一个复杂的动态时空中，由此得到的位移序列也具有较强的非线性和不确定性，增加了资料分析的困难度，本文将引入多重分形理论[3-5]来细致刻画位移序列的内在规律，从多重分形的视角描述大坝位移的变化性态。

分形是自然界的一种常见现象，具有整体与局部的自相似性，可以分为数学分形与统计分形两类，数学分形多由迭代产生，是理想化的严格自相似，并且具有无穷的层次结构，而统计分形则具有弱化的且层次结构有限的自相似性，相比于数学分形而言，其实际研究意义更突出，如海岸线轮廓、污染物浓度等。在大坝变形分形特征领域，赖道平等[6]应用重标度极差分析法(rescaled range analysis，R/S)对大坝变形资料进行了单分形分析，用Hurst指数、分形维数刻画了大坝的运行状况；YIN等[7]在单分形基础上，通过对比双对数散点图提出确定时序统计循环长度的方法，后在大坝监测领域得到推广；田振华等[8]利用R/S方法分析监测资料，通过分形维数将监测数据与大坝损伤联系起来，从分形维数的角度定量刻画大坝的结构损伤程度；赵卿等[9]利用去趋势波动分析法(detrended fluctuation analysis，DFA)对大坝多测点进行分析，从单分形角度刻画了位移序列的变化特征，取得了较好的效果。上述单分形理论往往只在宏观上刻画位移时序的分形特性，但实际大坝位移由于受多种因素的影响导致其实测数据往往具有复杂的多重分形特性，仅靠单分形理论无法完整刻画其内部规律，于是在单分形理论的基础上进一步引入了多重分形理论，该方法在大坝变形性态分析方面得到一定应用。胡江等[10]利用多重分形去趋势波动分析法(multifractal detrended fluctuation analysis，MF-DFA)对大坝位移序列和环境量进行了分析，对大坝变形多重分形特征和整体性态进行了辨识；韦武昌等[11]采用滑动平均优化传统的MF-DFA方法，并在某双曲拱坝的实测资料上进行了应用，对该坝的变形性态进行了分析；杨景文等[12]基于多重分形理论对大坝多个裂缝序列进行了多重分形分析，并探讨了裂缝序列多重分形特征的成因，均取得较好的效果。

目前已有研究主要基于实测数据进行单分形和多重分形分析，从而刻画大坝变形的演化规律和发展趋势，而在位移序列的多元多重分形分析、多重分形非对称性以及算法优化等方面探索较少。故本文以大坝实测位移数据为基础，引入滑动时间窗口对传统多重分形子区间划分方法进行优化，基于MF-DFA和多元多重分形去趋势波动分析法(multivariate multifractal detrended fluctuation analysis，MV-MFDFA)进行单测点和多元测点的联合多重分形分析，并用非对称多重分形去趋势波动分析法(asymmetric multifractal detrended fluctuation analysis，A-MFDFA)进一步分析位移测点不同变化趋势下的分形特征，实现大坝位移从整体波动到局部趋势的多重分形描述和解读，为多重分形与大坝位移性态分析的结合提供更详尽的信息。

1 多重分形计算方法

1.1 多重分形去趋势波动分析法

DFA[13]源于核糖核酸行为机理的研究，其优势在于能够消除非平稳时序的趋势部分，但其对于非平稳时序的分形特征刻画不够完整，后人在DFA的基础上引入多重分形的概念，将DFA改进为MF-DFA[14-15]，经典的MF-DFA计算步骤如下。

(1)

步骤 2给定一时间尺度s，将序列y(t)以s为标准进行等分，总共可划分为m个等长连续且不重叠的子区间，m=int(N/s)。从实际数据来看，N未必能够整除s，即可能出现尾部数据冗余，故此处采用逆序处理法来充分利用数据信息，即在正序划分的基础上，从序列的末端开始重复前述操作，由此得到2m个子区间。

步骤 3对各个子区间进行趋势拟合，并在原基础上减去趋势部分，得到相应残差序列，记为zv(t)。

(2)

步骤 4计算残差序列zv(t)的均方差F2(s,v)

(3)

步骤 5对F2(s,v)数据集取均值，按式(4)计算序列的q阶波动函数Fq(s)

(4)

q的取值与Fq(s)受波动影响的程度有关，原则上可取任意非零实数。特别地，当q=0时，有式(5)的极限形式：

(5)

传统的MF-DFA易在区间的划分中产生伪波动，对后续分析造成干扰，此外，时序长度未必能够整除子区间长度，造成数据的冗余，若将冗余部分忽略，会造成数据信息的丢失，若采用逆序处理法，则会打乱原来数据的顺序，影响信息的获取。鉴于此，考虑大坝变形数据的特殊性，此处采用滑动窗口[16-17]对传统MF-DFA划分子区间的方式进行优化，通过一定长度的窗口在序列上按照一定的步长进行滑动取值，以减少数据的伪波动并充分利用数据信息。设窗口长度为s，序列长度为N，取滑动步长为1，则运行一次得到的子区间数目为N-s+1，并用式(6)代替式(4)(5)。

(6)

通过上述步骤可以得到某一尺度s对应的q阶波动函数Fq(s)，变换s的取值，重复上述步骤可得一系列的s-Fq(s)点值，若该时间序列存在长程相关性，则Fq(s)与s间有式(7)所示的幂律关系：

Fq(s)∝sh(q)

(7)

对上式两边取对数，得到式(8)的形式：

lgFq(s)=h(q)lgs+lgb

(8)

式中：Fq(s)为序列的q阶波动函数；h(q)为相应的广义Hurst指数；b为常系数。作lgFq(s)-lgs双对数散点图，并进行拟合，其斜率即为广义Hurst指数h(q)，若h(q)为一固定常数，表明序列为单分形，无多重分形特征；若h(q)为q的非线性减函数，则位移序列具有多重分形特征，当h(q)<0.5时，位移序列表现为具有反持续性的记忆过程，当h(q)=0.5时，位移序列表现为不相关的随机过程，当h(q)>0.5时，位移序列表现为具有正持续性的记忆过程，另外，当h(q)>1时，位移序列表现为具有强非平稳性的长程正相关过程。

时间序列的分形强度和分形奇异性通常可由多重分形谱f(α)来表征，按式(9)～(11)求解：

τ(q)=qh(q)-1

(9)

α=τ′(q)

(10)

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

(11)

式中：τ(q)为Renyi指数，也叫标度函数，若其为q的非线性上凸函数，则位移序列具有多重分形特征，若其为q的线性函数，则位移序列表现为单分形特征，故该特征也常作为判断序列是否具有多重分形的判据。α为奇异强度，f(α)为多重分形谱，当α-f(α)呈单峰上凸状，形似二次函数时，说明位移序列具有多重分形特征；当α-f(α)聚集到一个点附近时，说明序列时单分形的。此外，常用分形谱宽度Δα和分形强度Δh来定量描述位移序列的多重分形特性，Δα主要反映大坝变形的奇异性和性态演变的空间差异性，Δh主要反映序列的分形强度，并与参数值成正比。

Δh=max[h(q)]-min[h(q)]

(12)

Δα=max[α(q)]-min[α(q)]

(13)

1.2 多元多重分形去趋势波动分析法

1.1节所述MF-DFA主要是针对单个测点的一元时间序列的多重分形分析，对于多个测点组成的多元时间序列集，此处引入MV-MFDFA[18-19]进行分析，设一个多元时序集{xt,i}，其中t取值范围为[1,N]，i表示测点个数，计算步骤如下。

步骤 4除去趋势部分，得到测点多元时间序列的局部方差F2(s,v)：

(14)

式中最大范数按式(15)计算：

‖X-Y‖=‖(x1,x2,...,xn)-(y1,y2,...,yn)‖=

max{|x1-y1|,...,|xn-yn|}

(15)

步骤 5按1.1节中式(6)求解q阶波动函数Fq(s)。

采用1.1节中的方法，可以得到多元序列整体的广义Hurst指数h(q)，当h(q)>0.5时，系统呈现出正相关性，即系统今后的增减与过去的增减正相关；当h(q)<0.5时，系统呈现出负相关性，即系统今后的增减与过去的增减负相关；当h(q)=0.5时，系统具有较强的随机性，相关性不明确。多元序列整体的多重分形谱计算方法同1.1节所述。

1.3 非对称多重分形去趋势波动分析法

传统的多重分形理论只能辨识序列的多重分形特征，但无法捕捉序列在不同趋势上的分形特征，此处引入A-MFDFA[20-21]进一步加强时序的非对称多重分形特性研究。该方法前五步同MF-DFA，拟从步骤6开始介绍。

步骤 6通过线性拟合估计原序列x(t)各个子区间的斜率bv，当bv>0时，记为正趋势，当bv<0时，记为负趋势，并按式(16)～(17)计算正、负趋势序列的q阶波动函数。

其中

分别表示正、负趋势序列的数目，当bv≠0时，二者之和为子区间数目。

若各类序列具有长程相关性，则有形如上文式(7)的幂律关系，同理两边取对数，可得式(18)的形式：

(18)

通过拟合得到相应的h(q)、h+(q)、h-(q)，若正、负趋势序列的广义Hurst指数图重合，即h+(q)=h-(q)，说明位移序列在多重分形特征上具有对称性，反之，则具有非对称性，即在位移序列处于不同趋势时，其表现出的多重分形特征也不尽相同。非对称多重分形分析正、负趋势序列的多重分形谱计算方法同1.1节中所述。

2 工程实例

某水电站枢纽布置如图1所示，选取图中河床坝段⑦号、⑧号、⑨号对应的EX5、EX6、EX7测点自2008年1月1日至2014年12月31日的坝顶水平位移监测数据，对各测点部分缺失数据进行3次样条插值，保证数据每天1个测值，去噪后的各测点实测过程线如图2所示，并规定向下游位移为正、向上游位移为负。

图1 某水电站枢纽布置图

图2 位移实测过程线

2.1 数据特征描述

在研究大坝位移的多重分形特征之前，应先了解样本位移数据的基本统计特征，此处对EX5～EX7测点数据进行描述性统计如表1所示。

表1 EX5～EX8测点位移描述性统计

根据K-S准则，考虑样本均值、标准差、偏度和峰度四个评价指标，各测点的位移均值均小于0，表现为向上游位移；从标准差来看，EX7标准差最大，EX5和EX6标准差接近，说明EX7位移序列偏离均值的离散程度比其他各测点大，而其他2测点的偏离程度和波动程度均较小，这一点也可从图2中直观印证；各测点的偏度均大于0，说明数据分布均为右偏态，即均值右侧的离散度比左侧强，表现为一定程度的右侧长尾，同时各测点峰度均显著大于0，说明数据高峰比正态分布的高峰更陡峭，处于尖峰状态，综合J-B统计量的结果，在峰度和偏度的共同作用下，各测点的位移序列均不服从正态分布；此外，对各测点进行长记忆性分析，所得长记忆性指标均大于0.5，说明各测点位移序列均具有较好的长程相关性，并非随机游走序列。综上所述，各测点位移序列均不符合正态分布特征，并呈现出Levy分形分布所特有的尖峰胖尾特点，表明位移序列具有明显的分形特征，可用分形理论的相关方法进行研究。

2.2 多重分形分析

采用滑动时间窗口优化MFDFA和MV-MFDFA，分别对测点和由测点构成的多元测点序列进行多重分形分析，此处波动阶数q取值范围为[-10,10]，尺度s的取值范围为[10,100]，滑动窗口步长为1，分析成果如图3～5所示。

图3 广义Hurst指数

图3和图4分别为测点整体及单个测点的广义Hurst指数和标度函数τ(q)变化情况，由图3可知，当q在[-10,10]之间变化时，无论是测点整体还是各测点分量序列，其广义Hurst指数均非常数，而是随着q的变化呈非线性递减趋势，说明从测点整体到内部各测点位移序列均具有明显的多重分形特性，仅通过单分形理论进行描述是不充分的；在不同的波动阶数q下，整体分析的广义Hurst指数曲线集中在各单测点的下部波动，表现出较弱的多重分形性，但无论整体或分量的h(q)值仍明显大于0.5，说明变形监测序列从整体到局部分量均具有较好的记忆性和长程相关性，兼具非平稳性和随机性；当q<0时，h(q)主要受小波动方差的影响，当q>0时，h(q)主要受大波动方差的影响，由图3可知，当q<0时，各测点h(q)值分布较离散，正持续性强弱差异明显，而当q>0时，各测点h(q)值一致性较好，具有相近的正持续性，此外，当q>0时测点间h(q)值的差异比q<0时小，说明在水位等大波动因素的影响下，各测点递减的趋势基本相同，整体变化趋势也相同，测点所在区域保持较稳定的长程相关性，坝体状态良好。在h(q)的基础上，通过式(9)计算Renyi指数即标度函数τ(q)，如图4所示，各测点标度函数一致性较好，中部呈上凸状，满足τ(0)=-1，且整体为非线性关系，这一点进一步印证了各位移测点具有多重分形特征。

图4 标度函数

图5为整体分析及各测点分量的多重分形谱，各多重分形谱图像呈单峰凸分布，形似二次函数曲线，位移多重分形的局部标度不是常数，刻画了不同时刻局部变化的多样性，奇异强度α主要集中于图像的两侧，反映了位移序列的分形结构分布不均匀，α分布的不均匀也印证了测点序列的多重分形特性；无论整体或测点分量的多重分形谱基本对称，整体协同性较好，发展状态稳定，其中EX6的分形谱呈现不明显的右钩状，说明对该点序列而言，小波动的影响略占优势；f(α)表征位移序列具有相同奇异指数α的子区间的分形维，与位移序列的分布特征和分形强度相关。由式(12)(13)计算多元序列及其分量的分形特征统计量如表2所示，从表2可以看到多元序列及其分量的分形强度Δh和多重分形谱宽度Δα均为最小，这是由于各测点位于同一高程的相邻坝段，测点变形之间并非相互独立的，而是存在一定的关联性，这种关联性削弱了多元测点整体的多重分形性，使得整体分析的分形特征统计量不会表现的比分量更强；整体来看，三个测点的Δα值接近，变形具有相似的奇异性和空间差异性，各测点水平位移的概率分布范围相对接近，但分形强度Δh值略有差异，表现为EX5、EX6的分形强度均强于EX7，结合实际情况初步推断，三个测点均位于河床中部同一高程的相邻坝段，变形均受水位波动等因素的影响，从而表现出一定的奇异性和多重分形特征，但由于EX5、EX6所在坝段的运行机理与EX7所在坝段不同，使得测点间的多重分形特征存在一定的差异。

图5 多重分形谱

表2 分形特征统计量

2.3 非对称多重分形分析

在上述研究的基础上，进一步引入多重分形的非对称分析法，即A-MFDFA，并同样采用滑动窗口进行优化，该方法考虑到位移序列在上升或下降趋势时的分形特征，故能够更全面刻画位移的分形特性和非对称性，正、负趋势序列的计算方法见前文所述，运行程序计算正、负趋势序列及原序列的广义Hurst指数、Renyi函数和多重分形谱，如图6所示。

图6 非对称多重分形分析

图6(a)(c)(e)为各测点不同趋势下的广义Hurst指数图，图6(b)(d)(f)为各测点不同趋势下的多重分形谱。从图6中可以看出，各测点位移原序列、正趋势、负趋势下的广义Hurst指数和多重分形谱均存在不同程度的偏差，图像并非重合的，说明位移序列在多重分形特征上具有明显的非对称性；位移序列的整体波动是正、负趋势波动综合作用的结果，正趋势和负趋势的h(q)均随q的增加而减小，同样表现出明显的多重分形特征，且h(q)均大于0.5，说明位移序列无论是处于正趋势还是负趋势时其发展变化仍具有良好的记忆性和长程相关性；综合各测点正趋势和负趋势的变化情况来看，q>0时h(q)的离散程度普遍高于q<0时的离散程度，说明水位等大波动因素对多重分形非对称性的影响最显著。

此处仍用多重分形普宽度Δα和分形强度Δh两个参数来定量评价位移序列的非对称多重分形特性，计算结果如表3所示。由表3可知，位移序列在各趋势上的分形特征均不相同，有明显的非对称特点，其中，EX5和EX6各自正、负趋势的分形特征接近，正负趋势各有强弱但分形强度差异较小，二者位移发展走势相近，相比之下由于所在坝段运行机理的差异，EX7测点正趋势的分形强度要明显大于负趋势的分形强度，原序列分形强度与正趋势分形强度接近，说明EX7位移序列在上升时的分形强度大于下降时的分形强度，即在上升时倾向于波动上升，而在下降时倾向于稳步下降。

表3 非对称分形特征统计量

3 结 论

本文将多重分形理论应用到某混凝土重力坝的变形性态分析中，并利用滑动时间窗改进传统多重分形子区间的划分方式，联合MF-DFA和MV-MFDFA对单测点和多测点序列进行多重分形特征辨识，并利用A-MFDFA进一步刻画序列局部处于上升或下降状态时的多重分形特征，分析结果表明，某大坝变形实测数据具有明显的分形序列特征，在水位等大波动因素的影响下，单测点和多测点均保持稳定的多重分形特征和长程相关性，坝体变形整体性态较好，处于稳定的安全状态。研究成果从变形的整体波动和局部趋势两个层面解析了大坝变形的多重分形特征及其非对称性，从一个新的视角看待大坝变形的波动特征，实现对大坝变形性态的准确把握，有利于大坝运行的管理和维护。