基于多重分形的混凝土重力坝水平位移波动分析

2022-02-12周兰庭柳志坤龚云柱

河海大学学报（自然科学版） 2022年1期

周兰庭，柳志坤，龚云柱

(1.河海大学水利水电学院, 江苏南京 210098; 2.青岛市经济发展研究院, 山东青岛 266000)

时间序列的研究对于揭示事物发展变化规律具有重要的意义，但由于事物的发展常常受多种因素的影响，导致实际时间序列往往是非平稳的，这给其分析研究带来了困难。为了有效挖掘时间序列的内在规律，人们引入分形理论来处理非平稳时间序列，并在环境保护、金融市场、微震监测等[1-3]领域广泛应用。推广到水利领域，混凝土重力坝的水平位移受库水位涨落、温度变化等因素[4]的影响，监测序列表现出高度的非线性和分形特征，分形理论最初在该领域的应用主要集中在大坝变形性态的单分形分析[5-6]，通过单分形来表征位移的长程相关性，描述时间序列波动的宏观特征，但其局限性在于无法对位移过程中某一时刻的局部特征或不同波动层次的复杂精细特征进行刻画，故引入多重分形来完善单分形理论。

在多重分形领域，Su等[7]应用多重分形去趋势波动分析法(multifractal detrended fluctuation analysis，MF-DFA)对某混凝土重力坝现场观测数据进行分析，确定现有大坝时间序列的多重分形标度行为，刻画了大坝的长期行为和结构演变规律；胡江等[8]利用MF-DFA对某混凝土重力坝监测序列的波动特征进行分析，从局部和整体两个波动层次评价了大坝的工作性态及其演变规律；杨景文等[9]基于MF-DFA研究了某混凝土重力坝多个裂缝序列的多重分形特征，并通过序列替换发现裂缝序列的多重分形特征主要由序列自身的长程相关性造成；韦武昌等[10]利用滑动平均法改进MF-DFA，并将其应用于某混凝土坝的变形分析中，探讨了多重分形强度与坝段运行机理的关系，对大坝的工作性态进行了安全评价；周斌等[11]将MF-DFA延伸到某混凝土重力坝的渗漏量监测分析中，通过计算渗漏量的定标指数确定渗漏监测序列的时程相关性，对多重分形在大坝变形分析中的应用有借鉴意义;周兰庭等[12]利用MF-DFA、MV-MFDFA和A-MFDFA解析了大坝变形性态的多重分形特征及其对称性。

前人的研究主要采用传统方法探讨混凝土重力坝监测序列自身的多重分形特征，未考虑环境量对监测序列多重分形波动特征的影响，胡江等[8]虽然在研究中引入了水位、温度等环境量因素，但并未将其与位移的多重分形特征建立联系。本文引入多重分形去趋势相关性分析法(multifractal detrended cross -correlation analysis，MF-DCCA)和多重分形非对称去趋势相关性分析法(multifractal asymmetric detrended cross-correlation analysis，MF-A-DCCA)，在考虑位移序列自身波动特征的同时，兼顾位移与环境量间的相关性，研究位移在环境量影响下整体波动与局部趋势的联合多重分形特征。

1 计算方法

1.1 多重分形去趋势相关性分析法

MF-DCCA[13-14]从多重分形视角考察两个时序的相关性，在金融、环境等众多实证研究中得到广泛应用。将该方法用于混凝土重力坝水平位移与环境量的分析中，算法解析如下。

步骤1记大坝实测位移时间序列为x(1)(t)，某一环境量实测序列为x(2)(t)，序列长度均为N。按式(1)计算相应序列的累计离差和y(i)(t)(i=1，2)：

(1)

步骤2给定时间尺度s，对y(i)(t)进行分割，得到h个长度为s的等长度不重叠连续子区间，其中h=int(N/s)。由于N未必能够整除s，会出现数据冗余问题，故采用逆序处理法回收数据信息，即在正序划分的基础上，从序列的末端开始重复前述操作，总计得到2h个子区间。

步骤3通过多项式拟合去除各子区间的趋势项，按式(2)计算去趋势后序列间的相关系数F2(s,v)：

(2)

式中：v为划分的子区间；pv(i)(k)为第v个子区间的趋势量，i=1，2。

步骤4计算两序列的q阶相关系数Fq(s)：

(3)

式中：q为波动阶数。当q>0时，表明大波动的影响占优势；当q<0时，表明小波动的影响占优势；特别地，当q=0时，可得极限形式：

(4)

步骤5若序列间存在多重分形相关关系，则有如下幂律关系成立：

Fq(s)∝sH12(q)

(5)

式中H12(q)为序列间的广义相关Hurst指数。对式(5)两边取对数得到式(6)，对式(6)最小二乘拟合，其斜率值即为H12(q)。

lgFq(s)=H12(q)lgs+lgC

(6)

式中C为常数。

观察拟合得到的H12(q)，若其不随q的变化而变化，说明序列间的相关关系是单分形的；若其随q的变化呈非线性变化，说明序列间存在多重分形相关关系。当H12(q)处于不同的区间范围时，可得到两序列间相关关系的不同结论：若H12(q)>0.5，两个序列间呈现正持续的相关性；H12(q)<0.5，两个序列间呈现反持续的相关性；H12(q)=0.5，则两个序列间的关系是随机、不确定的。

步骤6求解两序列间的相关性多重分形谱：

τ(q)=qH12(q)-1

(7)

α=τ′(q)

(8)

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

(9)

式中：τ(q)为标度函数；α为奇异性指数；f(α)为多重分形谱，其表征序列具有相同奇异性指数子区间的维数值。此外，常用分形谱宽度Δα和分形强度ΔH来定量描述位移序列的多重分形特性。Δα主要反映大坝变形的奇异性和性态演变的空间差异性，ΔH主要反映序列的分形强度，并与参数值成正比:

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

(10)

ΔH=max(H(q))-min(H(q))

(11)

MF-DCCA的计算流程如图1所示。

图1 MF-DCCA和MF-A-DCCA计算流程Fig.1 Calculation steps of MF-DCCA and MF-A-DCCA

1.2 多重分形非对称去趋势相关性分析法

MF-DCCA只能辨识序列间相关性的多重分形特征，但无法进一步刻画当位移序列处于上升或下降趋势时，其与环境量间相关性的多重分形特征，故引入MF-A-DCCA[15-16]来实现这一目的，算法步骤如下。

步骤1同MF-DCCA的步骤1，得到两序列的累计离差和y(i)(t)，i=1，2。

步骤2选取位移序列x(1)(t)、位移累计离差和y(1)(t)、某一环境量实测序列累计离差和y(2)(t)，采取MF-DCCA步骤2的方法，得到2h个子区间。

步骤3利用式(12)(13)对步骤2中的子区间进行趋势拟合，L为线性拟合函数，a为常数部分，b为拟合斜率，参数取值范围：v=1，2,…,2h；i=1，2；k=1，2,…,s。

Lxv(1)(k)=bxv(1)k+axv(1)

(12)

Lyv(i)(k)=bxv(i)k+axv(i)

(13)

步骤4对子区间进行去趋势处理，并计算协方差F2(s，v)：

(14)

步骤5考察式(12)的斜率值bxv(1)，当bxv(1)为正时，表明位移序列在该子区间为上升趋势，记为正趋势，反之，则记为负趋势。按式(15)(16)计算当位移处于不同趋势时，位移与环境量的q阶相关系数：

(15)

(16)

式中：M+为正趋势序列的数目，M-为负趋势序列的数目。当bxv(1)≠0时，二者之和即为子区间数：

(17)

(18)

步骤6若不同趋势下位移与环境量的多重分形相关关系仍存在，则有式(5)所示的幂律关系。同理两边取对数：

(19)

(20)

|ΔHA|越大，表明不同趋势下位移与环境量相关关系的多重分形非对称程度越高，趋势变化对多重分形强度的影响越显著。当|ΔHA|为正时，说明上升趋势下的位移与环境量相关性的持续性更强；当|ΔHA|为负时，说明下降趋势下的位移与环境量相关性的持续性更强。

1.3 多尺度滑动时间窗

传统多重分形分析易在区间划分中产生伪波动和冗余数据，若将该部分忽略，会造成数据信息的丢失；若采用逆序处理法，则会打乱原数据的顺序，影响信息获取。鉴于此，结合大坝位移监测数据的特殊性，采用滑动时间窗[17]对区间划分方式进行优化。

图2 滑动时间窗Fig.2 Sliding time window

滑动时间窗通常指长度一定、起始时间和结束时间变化的窗口。给定窗口初始长度为s，两端均以相同步长r推进，即一端加入r个数据，另一端依次分离出r个数据，运行原理如图2所示。滑动窗口在保持数据总量不变的同时，能够实时更新窗口数据，运行一次得到子区间的数目为(N-s)/r+1，通过改变s的取值，使滑动窗的长度随s的改变而改变，建立多尺度滑动时间窗模型。取r=1，并用改进后的式(21)～(23)代替式(3)(4)、式(15)(16):

(21)

(22)

(23)

2 实例应用

图3 测点布置Fig.3 Layout of measuring points

某水电站主要由混凝土重力坝、坝后溢流式厂房、埋设于坝内的输水系统、泄水底孔及过木筏等组成。大坝坝顶高程280 m，最大坝高78 m，坝顶全长253 m，水平位移引张线测点布设于坝顶处，测点布置见图3。选取河床坝段7～10号对应的引张线测点EX5～EX8自动化监测序列为样本，样本时间跨度自2008年1月1日至2014年12月31日，规定向下游位移为正、向上游位移为负，实测位移过程线和环境量如图4所示。

图4 实测位移过程线和环境量Fig.4 Hydrographs of measured displacement and environmental quantities

表1 主成分贡献率

由图4(a)可知，各测点监测数据同升同降，同方向、同高程、不同位置各测点位移间存在关联性。传统方法基于权重思想[18-19]建立各测点的权重关系，但计算量大且权重确定方法有待验证。而运用PCA[20-21]进行数据处理时，可保证在有效信息丢失最少的情况下，提取最能表征大坝多测点变形整体关系的主成分综合效应量，该量是各测点效应量的线性组合，可将其视为单测点效应量建模分析，降低了分析的复杂度。表1为PCA的主成分贡献率。

根据文献[20-21]，PCA主成分的阈值为85%，由表1可知，第一主成分贡献率1已达到PCA的阈值，能够在反映有效信息的同时将数据最简化，故选择第一主成分表征多测点位移信息。根据SPSS分析给出的评分系数，带入式(24)求解多测点位移主成分。

F=m1Y1+m2Y2+m3Y3+m4Y4

(24)

式中：F为主成分；Y1，Y2,Y3,Y4分别为位移；m1，m2,m3,m4分别为评分系数。

2.1 水平位移与环境量的多重分形相关性分析

将MF-DCCA应用到水平位移与环境量的多重分形相关性分析中，在matlab环境下运行MF-DCCA程序，分别计算该混凝土重力坝水平位移与环境量的波动函数Fq(s)、广义相关Hurst指数H12(q)、标度函数τ(q)和多重分形谱f(α)，如图5、图6所示。

图5 位移与水位的多重分形相关性分析Fig.5 Multifractal correlation analysis of displacement and water level

图6 位移与温度的多重分形相关性分析Fig.6 Multifractal correlation analysis of displacement and temperature

a.结合图5(a)和图6(a)，在q一定时，波动函数Fq(s)与s近似呈正相关关系。整体来看，q越小，Fq(s)受波动影响程度越大，位移与环境量相关性的波动随着s的增大而逐渐有序，并趋于均化。

b.结合图5(b)和图6(b)，位移与环境量的广义相关Hurst指数H12(q) 并非恒定不变，均随着q的增加而减小，说明位移的多重分形复杂性与水位、温度的影响均有关。位移与环境量的H12(q)均大于0.5，说明位移与环境量间的相关性具有明显的正持续性特征，总体表现出良好的长程相关性，即位移序列数值的上升或下降伴随着影响因子的上升或下降。当q较大或较小时，位移与温度间的H12(q)变化趋势相较于与水位间的H12(q)变化趋势更剧烈，说明位移与温度间的多重分形相关性对波动的敏感性更强。但从定量的角度来看，位移与水位间H12(q)明显大于位移与温度间的H12(q)，说明位移与水位间的正持续性强于位移与温度间的正持续性，水位变化对位移的影响更显著。

c.进一步分析图5(c)和图6(c)，τ(q)呈非线性递增趋势，且图5(c)的上凸程度和非线性程度均强于图6(c)，说明位移与环境量间的相关性具有多重分形特征，且位移与水位的这一特征要强于位移与温度的组合，进一步印证了前述分析。

d.结合图5(d)和图6(d)，位移与环境量的相关性多重分形谱图像呈单峰凸分布，刻画了不同时刻局部变化的多样性。α主要集中于图像的两侧，α分布的不均匀也印证了位移与环境量间的多重分形相关性。相关性多重分形谱基本对称，整体协同性较好，发展状态稳定，其中位移与温度的相关性多重分形谱呈现不明显的右钩状，说明在该情况下小波动的影响略占优势。进一步计算位移与环境量间的分形特征统计量，即位移-水位：ΔH=1.255 2，Δα=1.596 4；位移-温度：ΔH=0.476 8，Δα=0.713 8，位移与环境量间的ΔH和Δα均不为0。由于Δα越大，相关性的分形奇异性越强，ΔH越大，相关性的分形强度越强，故可知此时水平位移与水位间的多重分形相关性显著强于位移与温度间的多重分形相关性，水位变化对水平位移的发展影响更显著。加之水位的涨落变化在周期性上不如温度的变化稳定，属于高频非线性影响量，故其与位移间的多重分形相关性往往更复杂，波动性更强。

2.2 水平位移与环境量的非对称多重分形相关性分析

基于MF-A-DCCA进一步探究当位移序列处于正趋势或负趋势时，其与环境量间相关性的多重分形特征，分析成果见图7、图8，非对称分形特征统计量见表2。

图7 位移与水位的非对称多重分形相关性分析Fig.7 Correlation analysis of asymmetric multifractal of displacement and water level

图8 位移与温度的非对称多重分形相关性分析Fig.8 Correlation analysis of asymmetric multifractal of displacement and temperature

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

由图7、图8和表2可知，随着q的变化，当位移分别处于正趋势和负趋势时，其与环境量的广义相关Hurst指数曲线互不重叠，说明此时位移与环境量间相关性的多重分形特征具有明显的非对称性；无论是位移原序列、正趋势还是负趋势，其与环境量间的H12(q)并非定值，而是随着q的变化而变化，且均大于0.5，说明位移与环境量的相关性无论在位移序列上升还是下降时都具备多重分形特性，同时兼具良好的长程相关性。从表2的结果可知，位移与环境量的正趋势分形强度均大于负趋势分形强度，说明当位移序列处于上升趋势时，其与环境量间的正持续性比处于下降趋势时的正持续性强；进一步对比表3同一序列类型下位移-水位、位移-温度的计算结果可知，无论原序列、正趋势还是负趋势，前者的分形奇异性和分形强度都明显高于后者，说明即使在不同的变化趋势下，水位变化对位移多重分形特征的影响仍占优势，与前文结论对应。

图9 位移与环境量相关性的非对称度量Fig.9 Asymmetric metric of displacement and environmental quantities

3 结论

a.运用多重分形理论刻画了混凝土重力坝水平位移与环境量相关性的多重分形特征，并研究了其非对称性，利用广义相关Hurst指数、相关性多重分形谱等分析了多重分形相关性的强度，剖析了在环境量的影响下，大坝变形性态的动力学特征和演变规律。

b.引入多尺度滑动时间窗模型消除由子区间划分而产生的伪波动，相较于文献[9]中未优化的情形，优化后所绘图像平滑性较好，趋势性明显，减小了伪波动的影响，提高了多重分形分析的信息利用率。

c.实例分析表明，该混凝土重力坝水平位移与环境量的相关性具有较强的多重分形特征，且这种特征在位移处于不同趋势时具有非对称性，位移变化从整体波动到局部趋势均表现出良好的记忆性和长程相关性。水平位移的多重分形特征受环境量的影响，其中水位的频繁涨落起主导作用，在环境量的影响下，该混凝土重力坝顶水平位移性态良好，分析结果符合实际工程情况，表明本文的方法在混凝土重力坝水平位移的分析中具有较好的可行性和适应性。