基于降维分析的大坝安全特征提取方法及其应用研究

2019-01-23方卫华王润英孙一清

人民珠江 2019年1期

方卫华，王润英，孙一清

(1.水利部南京水利水文自动化研究所，江苏南京 210008 ;2.河海大学水利水电学院，江苏南京 210098)

近年的溃坝失事事件表明，大坝安全风险依然存在。对于运行期的大坝安全问题，目前的安全特征提取方法主要包括基于实测资料的分析方法和数值模拟方法。前者主要包括特征值统计、回归分析和监控指标拟定等，后者主要采用数值模拟法、刚体极限平衡法[1]和强度折减法等获取强度和稳定参数。尽管上述研究已取得很多研究成果，但前者依然受到小样本、多重共线性、非正态分布、测值误差以及分析过程和结果不直观等问题的困扰，而后者由于受计算复杂度加卸载过程等因素的限制，选用的计算模型及参数往往难以完全反映客观实际[2]，特别是难以对材料的本构模型、结构材料分区和多场耦合问题进行准确的考虑。随着大坝监测自动化技术的应用和数据采集时间长度的积累，监测数据多、维数高、计算分析不直观等问题更加突出，而管理现代化对分析的实时性却提出更高的要求。随着多种降维分析方法的不断提出和应用，为解决工程实际问题，特别是信息的传输和可视化问题，对现有的大坝安全高维实测数据或无穷维大坝基于降维的大坝安全特征提取不仅是必要的，而且是可行的。本文在特征正交分解(POD)、等概率粗粒化、主成分/核主成分分析(PCA/KPCA)分析的基础上，将这些降维分析方法与大坝安全特征提取方法结合，提出改进的新方法，并结合典型工程实例检验新方法的有效性。

1 研究现状

1.1 降维分析方法

降维分析是指通过数学、数值或统计方法对高维信息进行低维变换或通过剔出冗余、误差或无用信息，从而获得描述结构整体特征的低维(阶)模型、特征向量或阈值的过程。常用的降维方法主要包括矩阵非负分解、自组织映射、等距映射、局部线性嵌入算法、主成分分析和核主成分分析等，本文涉及的主要降维方法包括特征正交分解(POD)、等概率粗粒化、主成分/核主成分分析(PCA/KPCA)、矩阵范数(如F-范数和谱范数)[3-4]。上述方法介绍的文献较多[5-6]，为节约篇幅，不再赘述。

1.2 大坝安全特征提取方法

特征提取是指基于实测数据或者计算分析，提取表征大坝安全状态的模型、参数向量或阈值的过程。特征提取方法主要包括基于实测数据的回归分析、统计分析、矩阵变换和监控(预警)指标拟定，以及针对大坝结构本身的刚体极限平衡法、数值仿真模拟和强度折减法等。

1.2.1 回归分析

回归分析通过拟合样本来确定自变量和因变量之间的统计关系，通过回归系数向量或方程获得了隐含在高维样本内的大坝安全信息。然而，当前常用的基于最小二乘法的逐步回归分析提取的是一段时间序列的平均信息。实际上，大坝失事是一个典型极端情况下的小概率事件，尾部特征更加重要，因此必须提出考虑最不利工况组合及极端情况下的大坝安全特征提取方法[5-6]。

1.2.2 预警指标拟定

监控指标的拟定实际上也是通过统计分析或数值模拟计算提取实测资料的极值分布、低维统计边界[7]或屈服破坏临界值[8-9]。基于实测资料的监控指标分析目前均利用单一监测变量长时间多维向量，没有考虑到多个变量之间信息的互补性，而基于数值模拟的方法基本都是单调加荷和理想弹塑性模型，与实际情况的多场耦合相差较大。

1.2.3 复杂网络分析

时间序列复杂网络分析是一种直观的特征提取方法，将实测资料转化成复杂网络后，通过最短路径、平均聚集系数、介数等复杂网络有限低维参数提取时间序列隐含的特征[10]。

1.2.4 安全系数获取

刚体极限平衡法和强度折减法都是获取大坝抗滑稳定安全特征提取的常用方法，由于极限平衡法没有考虑材料的特性和应力～应变关系以及边界条件的变化，不能反映应力分布变化及屈服过程，因此逐渐被强度折减法取代。但当前强度折减法主要是针对单一荷载[11]，没考虑多场耦合。实际上大坝都在多场耦合作用下工作，而考虑多场耦合，尤其是全耦合方法计算工作量比较大，因此必须解决数据计算工作量问题。

2 基于降维分析的特征提取方法

2.1 基于POD的分位数回归分析法

为考虑极端水位或气温等自变量对大坝安全的影响，深入挖掘自变量尖峰厚尾等非正态分布以及不同分位点所包含的安全信息，决定选用分位数回归模型对实测数据进行分析，分位数回归具有如下优点：①对模型中的随机扰动项无需做任何分布的假定，整个回归模型就具有很强的稳健性；②本身没有使用一个连接函数来描述因变量的均值和方差的相互关系，因此具有比较好的弹性性质；③由于对所有分位数进行回归，因此对于数据中出现的异常点具有耐抗性；④不同于普通的最小二乘回归，分位数回归对于因变量具有单调变换性；⑤分位数回归估计出来的参数具有在大样本理论下的渐进优良性。

为得到回归系数，均值回归可采用最小平方法，中位数回归可采用最小一乘法(最小绝对偏差法)，对于分位数回归可采用线性规划法等。为进一步减少自变量的个数和克服自变量之间的多重共线性，将POD的小样本逼近、特征提取与分位数回归的优势结合起来，首先对实测数据求取特征正交基，再利用特征正交基进行分位数回归分析，从而获得不同分位点上各综合自变量对因变量的影响程度。

2.2 基于实测资料的复杂网络分析法

该方法可以有效克服离散误差对数据分析结论的影响，实现数据分析的可视化，具体步骤如下。

a) 选取典型和重要的大坝测点及其测值。

b) 运用等概率粗粒化的方法，将测值时间序列转化成符号序列[10]。①设xmax和xmin分别是测值时间序列的最大值和最小值，转化后的符号共有N种，记为s1，s2，…，sN，定义:

Si=Sj,xmin+(j-1)d≤xi

(1)

②记大坝测值时间序列为{x(l),l=1,…N}，通过式(2)计算测值时间序列波动k(l)：

(2)

式中 Δl——时间间隔；

③通过式(3)计算不同波动值可能出现的概率p(k):

(3)

式中Num(x)——测值序列的波动模态x发生的次数；

④把测值时间序列的波动k(l)分为5个区间，定义5个特征字符为式：

(4)

式中t——测值数值快速增加；r——测值数值缓慢增加；e——测值数值不变；d——测值数值缓慢降低；f——测值数值快速降低。从而将测值时间序列转换成符号序列为：

S=(S1,S2,S3,…)，Si∈(t，r，e，d，f)

(5)

c) 将符号序列编码通过excel2pajek软件转换成pajek软件可处理的格式构造复杂网路。

d) 分析复杂网络的特征参数，获取该大坝监测项目的敏感测点及其性态演化规律，判别大坝的工程安全性态。

2.3 基于PCA/KPCA的预警指标提取法

置信区间法、典型监测效应量小概率法、极限状态法和极值理论等都是针对单个测点提出的，没有考虑多个测点之间的信息冗余或互补，更没有考虑土石体在失稳破坏中的应变局部化等效应。本文将多个监测项目和测点数据按时间次序形成数据带，选择一定时间长度的窗口，根据主成分和核主成分算法，在85%贡献率的基础上得到PCA/KPCA主成分个数和对应的矩阵范数2个指标作为预警指标。

核函数常用选取方法有Cross-Validation法、多次实验法和混合核函数方法，通过对比分析，本文选择高斯径向基核函数进行核主成分分析[5]。

2.4 基于POD的多场耦合的强度折减法

由于温度荷载对大坝抗滑稳定安全系数影响不大，暂不考虑温度荷载，考虑饱和-非饱和流固耦合方程：

(6)

式中uw,j——孔隙压力；Sw——饱和度；εv——土骨架体应变；kwij——渗透系数；ρw——水的密度；bj——单位体积力；n——孔隙率。坝体及坝基采用非关联准则的弹塑性本构模型、Mohr-Coulomb屈服准则：

(7)

(8)

基于POD的强度折减法首先利用大坝安全监测资料构造瞬像集合，采用POD方法提取特征正交基，结合Galerkin投影法建立有限元方程[12-13]，再结合强度折减法即可获得大坝抗滑稳定安全系数，从而可以大幅度节约计算时间，提高计算效率。

3 工程实例

本文选择坝型分别为常态混凝土重力拱坝、均质土坝和深厚覆盖层上碾压混凝土重力坝3个典型工程实例检验上述方法的有效性。

3.1 安徽某重力拱坝

3.1.1 工程存在问题

3.1.2 研究内容和结果

3.1.2.1 基于POD的分位数回归分析

a) 拟合精度分析：通过对发现裂缝过程中间时段拱冠梁和左右1/4拱坝顶测点径向分析，基于POD的分位数回归分析方法可以很好地适应小样本性质和非平稳性，拟合精度较高，参见图1—3。

b) 大坝安全特征变化分析：将发现裂缝的时段分成前后2个时段，分别进行回归分析，通过模型参数前后变化分析裂缝对大坝安全性态的整体影响。前一时段典型测点回归方程为：

(9)

后一时段回归方程为：

(10)

通过上述系数矩阵可见，采用相同的基Ψ1、Ψ2，裂缝存在前后2个时段大坝典型测点的径向变形并未受到明显影响，说明大坝整体安全性态稳定，这一结论与大坝仿真计算和专家综合评价结果一致。

3.1.2.2 基于实测资料的复杂网络分析

同样利用8、18、26号坝块坝顶径向实测位移数据，分成前后两段，分别采用复杂网络进行分析。

通过图4可知，裂缝存在前后2个时段大坝典型测点的复杂网络未出现显著变化，说明大坝整体安全性态稳定，这一结论与上文分析方法得到的结论一致[14]。

3.2 某均质试验土坝

3.2.1 研究目的

为研究溃坝过程中大坝安全性态的变化，建造了一座水库并进行了溃坝试验。该试验水库最大库容10万m3，大坝为黏性均质土坝，坝高9.7 m，坝顶长120 m、宽3 m。坝体分层碾压填筑。试验过程包括大坝填筑(仪器埋设)、首次蓄水、水位稳定、第2次蓄水和溃决等过程，共历时46 d[14]。

3.2.2 研究内容和结果

3.2.2.1 基于POD的分位数回归分析

假设实测数据分析前，并不知晓溃坝开始时间，为了检验本方法正确性，将得到的变形数据分成前80组进行基于POD的回归分析。前80组位移进行模拟(图5)。

由图5可知，模型的拟合值与实测值拟合较好，表明所建模型能够较好地反映前80组位移监测值的变化。为了检验模型的预测能力，预测了后15 d的位移值，见图6。

由图6可知，前29组预测值与实测值相差较小，说明所建立的模型精度较高，能够较好的预测后期数据。而前29组数据之后预测值与实测值相差较大，表明此时大坝开始处于溃坝前期，濒临溃坝。第59组数据的模型预测值突然跳动很大，此时已经溃坝。

将上述第59组数据开始将数据分成两段，得到位移与POD正交基的回归方程分别为：

y=-6.841Ψ1-9.151Ψ2-168.16

(11)

y=-3.286Ψ1-5.473Ψ2-288.52

(12)

式(11)与式(12)存在明显差异，说明通过降维分析获得的方程系数反映了土坝安全性态的变化。

3.2.2.2 基于实测资料的复杂网络分析

根据上文时间结点，将土坝渗流实测资料分成溃坝前后两段，对应的复杂网络见图7、8。

通过图7、8不难发现，两者之间存在明显差别，也就是说基于复杂网络的降维分析方法能发现的溃坝时间与实际一致。

3.2.2.3 溃坝预警指标拟定

土坝溃坝实验过程见文献[15]，通过在溃坝过程中埋设安装在坝体上的63台仪器进行主成分分析。将63个测点每连续两天(每天48次)的测值构成一个窗口大小为96×63的数据窗口，其中63为测点数，对应列向量个数；96为测次，对应行向量个数。从2012年10月16日—2012年11月17日共32个窗口。

应用PCA对监测数据进行主成分分析，结果见图9，其中横坐标为时间窗口，纵坐标为主成分累计贡献率(下同)。通过图9可知，10月16—31日满足一定累计贡献率的主成分个数一直增加，10月31日至11月14日逐渐下降，在11月17日又突然增加。同时第一主成分贡献率呈现先缓慢减少、再缓慢增加又突然减少的变化过程。可见主成分个数反映了大坝安全性态变化，随着大坝安全性态的变化，主成分个数先增加、后减少，而在溃坝时突然增加。根据分形理论，在溃坝前主成分数量减少对应了大坝整体降维和应变局部化过程，而17日主成分个数突然增加则反映了大坝结构的离散化和增维过程[15]。

一定累计贡献率下，96×63窗口核主成分个数随时间变化过程线见图10。图10与图9结果基本一致，10月16—30日核主成分个数一直增加，11月3—14日逐渐下降，在11月17日又突然增加。同时第一主成分贡献率呈现先缓慢减少、再增加又突然减少的趋势。即随着大坝安全性态的变化，主成分个数先增加、后减少，而在溃坝时突然增加。

由图9、10可知，总体上10月17—31日主成分个数一直增加，10月31—11日逐渐下降，在11月17日又突然增加。同时主成分范数也呈现先缓慢减少、再增加又突然减少的过程，上述过程和溃坝过程吻合。在统计主成分个数的同时还计算相应的矩阵。通过图11可知，尽管采用不同的范数计算，但计算结果基本一致，即范数的变化对具体范数计算方法不敏感，因此主成分范数可以作为溃坝分级预警指标。

3.3 山西某碾压混凝土重力坝

3.3.1 工程存在的问题

大坝为深厚覆盖层上碾压混凝土坝，由于在大坝主体工程竣工时坝体坝肩固结和帷幕灌浆未完成，导致廊道内漏水严重。2014年开始进行以坝肩帷幕灌浆、坝基固结灌浆和坝后连续墙施工等应急除险加固措施。为了解除险加固后大坝的安全性态，采用数值方法进行分析研究。

3.3.2 研究内容和结果

大坝埋设安装了大量的位移、渗流、温度和应变传感器，为基于POD的多场耦合方法提供了条件。分别用除险加固前后的多场实测数据建立瞬像集合，采用基于特征正交分解方法提取正交基，并将正交基作为试函数，利用有限元法获得大坝除险加固后的变形、渗流分布和抗滑稳定安全系数。表1给出了不同计算工况下采用POD的多场耦合方法与直接有限元法的计算机机时比统计。