基于监测时序分解再重构的混凝土拱坝位移预测组合模型

魏博文，罗绍杨，徐富刚,2*，袁冬阳，张婉彤

(1.南昌大学 工程建设学院，江西 南昌 330031；2.江西省水利科学院，江西 南昌 330029；3.河海大学 水文水资源与水利工程科学国家重点实验室，江苏 南京 210098；4.河海大学 水利水电学院，江苏 南京 210098)

服役期混凝土坝长期受静水荷载、温度等循环荷载的作用，且难免遭遇地震等突发灾害及极端气候影响，从而导致筑坝材料的老化与结构性能的退化，进而影响混凝土坝的运行安全与服役寿命[1-2]。变形是最能直观反映混凝土坝服役性态的综合效应量，基于实测变形的健康监控模型是合理解译环境荷载对大坝运行性态影响、评估大坝运行安全与预测未来服役性能的重要途经[3-5]。由于混凝土坝变形机理的复杂性与测值的高度非线性，现有监控理论与方法常难以完全挖掘变形测值的有效信息，造成监控模型的预测能力欠佳。为此，融合多元数据处理手段，深入挖掘测值内蕴有效成分，构建具有高精度的监控模型具有重要意义。

混凝土坝变形主要是由水位和温度变化引起的可逆变形和筑坝材料性能随时间演变引起的不可逆变形组成[6-8]。根据建模方法分类，变形监控模型可分为统计模型、确定性模型和混合模型3类[9]。因统计模型具有函数形式简单、计算高效的优点，在实际工程领域取得了广泛的应用。近年来，伴随人工智能技术的高速发展，诸如人工神经网络（ANN）[10-11]、支持向量机（SVM）[12-14]、极限学习机（ELM）[15-16]、长短期记忆网络（LSTM）[17]等大量的机器学习算法及小波分析（WA）[18]、主成分分析（PCA）[19]、奇异谱分析（SSA）[20]、经验模态分解[21]等诸多数据处理方法被广泛应用于大坝变形性态分析与统计模型中，极大促进了大坝变形监控理论的发展。Shao等[22]将面板数据理论引入大坝变形分析中，从空间和时间两个维度并基于面板数据的随机系数模型建立不同测点之间的潜在关系，解决了传统回归方法的多重共线性问题；de Granrut等[23]通过观察模拟图形结果，选择了人工神经网络（ANN）的最佳参数，以防止模型过度拟合；Zhang等[24]使用改进的粒子群优化算法（PSO），对极限学习机（ELM）的参数进行优化，该模型提高了算法的泛化能力；李明超等[25]在贝叶斯框架下，以加法模型为基础，重构各时序分解项作为模型底层，通过结合参数化检测与直观参数配置，实现交互式建模；Cheng等[26]将支持向量机（SVM）与主成分分析（PCA）相结合，对棉花滩重力坝的响应进行短期预测。

然而，从本质上讲，统计模型为通过构建大坝变形与环境荷载之间数学函数的经验模型，无法有效考虑筑坝材料与结构性能演变对大坝变形的影响。相比于重力坝而言，拱坝这种超静定结构对环境荷载变化尤为敏感，其变形行为更为复杂。限于当前监测技术，大坝变形测值难免存在一定的噪声成分和混沌效应[27]。同时，现有监控模型仅考虑库水位、环境温度及时效变化对大坝变形的影响，然而，大坝实际运行环境复杂，尚存一定的环境影响因素难以有效考虑。Wei等[28-29]研究表明，常规监控模型残差序列常表现出一定的周期性效应，结合多种数据挖掘方法，进一步挖掘残差序列内蕴的有效信息，可有效提升模型的拟合与预测精度。

鉴于单一挖掘手法难以完全有效地提取大坝变形测值的有效信息，为提升监控模型的拟合与预测进度，本文在构建变形预测混合模型的基础上，结合有限元方法计算水压分量进而构建混合模型；同时，考虑到混合模型残差序列的混沌与周期性特征，应用集合经验模式分解（ensemble empirical mode decomposition，EEMD）将其解构为具有不同时域特征的本征模态分量，进一步综合粒子群优化（particle swarm optimization，PSO）的极限学习机（extreme learning machine，ELM）和季节性差分自回归积分滑动平均模型（seasonal autoregressive integrated moving average，SARIMA）模型，对所解构的高频与低频信号分别建模预测，据此，构建了变形组合预测模型。

1 混凝土拱坝变形混合预测模型

由混凝土坝变形的成因可知，拱坝径向位移主要由水压分量δH、温度分量δT与时效分量δθ3个部分组成。变形水压分量δH是坝体与坝基在库区静水荷载作用下产生的位移，根据其力学机理，水压分量δH可分为3部分：静水压力作用在坝体上引起坝体自身的位移 δ1H，及其作用在地基面上产生的内力，导致坝基变形引起的位移 δ2H，库水重作用使得地基面转动所引起的位移 δ3H，如图1所示。

图1 水压分量组成原理示意图Fig. 1 Principle diagram of water pressure component composition

对于混凝土拱坝，水压分量常采用上游水深H的4次多项式计算：

温度分量δT为坝体混凝土温度变化引起的热变形，根据温度计的布设情况，δT有两种计算形式：

1）当坝体和基岩内布设了足够数量的温度计时，其可表示为：

式（2）和（3）中，b1i、b2i为拟合系数，t为监测日至始测日的累计天数。

时效分量综合反映的是坝体混凝土与基岩徐变、塑性变形以及碱骨料反应和其他可能威胁结构完整性的非线性时变效应，通常δθ由式（4）计算：

式中，c1、c2为拟合系数，θ=t/100。

式中，α0为常数，H0为建模序列初始日上游水深，t0为建模序列初始日至始测日的累计天数，θ0=t0/100。

2 混合模型残差序列的混沌效应预测

鉴于大坝变形机理的复杂性与测值的高度非线性，单一监控模型常难以完全挖掘大坝变形的有效信息，从而限制模型的预测能力。为提升模型的预测精度，本文采用EEMD对混合模型残差序列δr进行多尺度解构，结合解构信号的时域特征，分别利用PSOELM和SARIMA模型，对模型残差有效成分进行挖掘与预测，提出了融合残差有效成分的混凝土拱坝变形预测混合模型，研究总体框架结构如图2所示。

图2 总体框架结构Fig. 2 Overall framework structure

2.1 基于EEMD的混合模型残差多时域信号挖掘

EEMD是一种挖掘非线性、非稳态时间序列信息特征的常用处理方法，其通过在分解过程中多次引入均值为0、标准差为常数的高斯白噪声ni(t)掩盖信号原有的噪声，进而得到更为精准的包络线，有效地解决了EMD等方法存在的模态混叠问题。

式中，cj(t)为EMD分解得到的第j个IMF分量，r(t)为分解得到的余量。

为消除高斯白噪声的多次引入对实际IMF的影响，将上述步骤得到的IMF，取均值得到：

2.2 基于PSO-ELM的模型残差高频信号预测

2.2.1 ELM基本原理

ELM是融合神经网络与支持向量机优势的一种单隐层前馈神经网络算法。该方法无需调整网络的输入权值以及隐含层阈值，仅需设置网络的隐含层节点个数，故ELM具有训练速度快、训练过程简单、调节参数少等优点。ELM原理如图3所示。

图3 ELM原理示意图Fig. 3 Schematic diagram of ELM principle

式中，H+为 隐藏层输出矩阵H的Moore-Penrose广义逆矩阵。

2.2.2 基于PSO优化的ELM模型

鉴于ELM算法中输入层权值和隐含层阈值对模型的训练性能影响较大，为提高模型的预测精度，本文利用PSO寻求上述参数的最优解。在PSO算法[30-31]中，每个粒子的空间坐标对应1组ELM模型中的输入层权值和隐藏层阈值的可能解，通过粒子的空间位置迭代更新计算种群的适应度值，即可得到使得适应度值最优的一组解。粒子的速度与位置更新策略为：

式中：Vi为粒子的速度；Xi为粒子的空间位置；Pg为当前迭代的最优解；Pi为当前最优解的个体极值；w为惯性权重系数；d为空间的维度，d=1,2,···,D；i为当前计算的粒子编号，i=1,2,···,N；k为当前迭代次数；c1和c2为加速度因子，一般为非负的常数；r1和r2为分布于[0,1]之间的随机数。

采用PSO-ELM算法参数最优解，通过构建适应度函数，经粒子迭代优化，寻求适应度函数最优值对应的搜索代理的位置向量加以实现。本文采用误差平方和作为适应度函数，即

PSO-ELM的构建步骤如下：

Step1：训练数据归一化处理与模型参数初始化。确定PSO种群规模N，加速度因子c1和c2，最大迭代次数Tmax及ELM模型中待求参数的取值范围。

Step2：随机生成初始粒子空间位置，并根据式（21）计算其初始适应度值。

Step3：根据式（20）更新粒子的空间位置，并计算其适应度。若当前迭代粒子空间位置对应的适应度优于前期粒子适应度，则保留当前粒子空间位置；反之，则保留前期最优粒子空间位置。

Step4：PSO迭代寻优过程中，当达到最大迭代次数时，迭代终止，输出最优参数，否则重复Step3。

Step5：通过PSO优化得到的最优参数，建立PSOELM模型，进而对混合模型残差高频信号进行预测。

2.3 基于SARIMA模型的拱坝变形残差序列的低频信号处理

鉴于EEMD分解得到的混合模型残差的低频信号呈现出一定的季节性变化特征，本文采用具有良好的季节性、随机性时间序列处理能力的SARIMA模型，建模预测。SARIMA模型[32]有效综合了随机季节性模型和ARIMA模型的优点，可对既含有趋势成分（T）又含有季节性成分（S）的时间序列进行有效份分析与预测。SARIMA模型可写作SARIMA（p,d,q）×（p',d',q'）S模型，其中，p和q分别为非季节性自回归（AR）过程和移动平均（MA）过程的阶数，p'和q'分别为季节性自回归（SAR）过程和移动平均（SMA）过程的阶数，S为季节性的周期。ARIMA（p,d,q）模型的结构为：

在应用SARIMA模型时需对模型参数予以确定，采用赤池信息准则（Akaike information criterion，AIC）定阶对其阶数进行定阶，AIC准则的定阶方法为：

模型参数的定阶方法可简要概括为：

Step1：混合模型残差低频序列平稳化处理。观察差分后时序的自相关函数（ACF）和偏相关函数（PACF）是否截尾判定序列的平稳性，直至d次差分变换后时序平稳。

Step2：SARIMA模型识别。对混合模型残差低频序列做季节性检验，对可能的p、d、q和p'、d'、q'参数值进行估计，构成不同的备选SARIMA模型集。

Step3：选取最优模型参数进而建模预测。

3 拱坝位移监测混合模型的程序实施方法

融合残差序列有效成分的混凝土拱坝变形混合预测模型的实现流程如图4所示。

图4 组合模型构建流程图Fig. 4 Flowchart of composite model construction

采用平均绝对误差（mean absolute error，MAE）、均方误差（mean square error，MSE）、平均绝对百分误差（mean absolute percentage error，MAPE）以及拟合优度（coefficient of determination，R2）4个衡量指标，综合评估组合模型的拟合精度，并避免单个指标反映总体监测数据的数量特征单一性。其计算如下：

4 实例分析

4.1 工程背景

某混凝土双曲拱坝最大坝高99.1 m，坝顶宽5.0 m，坝顶高程247.6 m，坝顶长287.94 m。该坝正常蓄水位与校核洪水位高程分别为244.0 m和246.72 m，调节库容与总库容分别为0.62×108m3和1.05×108m3。为监测该坝坝体与坝基的水平变形，该坝布置了3组垂线监测仪器，包括6条正垂线和4条倒垂线，布设方案如图5所示。

图5 某混凝土拱坝测点布置图Fig. 5 Layout of measuring points for concrete arch dam

以位于拱冠梁坝段的正垂测点PLA1为例，选取2015-05-01—2020-02-13实测径向位移加以分析。其中，以2015-05-01—2019-11-16监测数据用以模型训练，剩余数据用于检验模型的预测性能。监测时段内测点的径向位移与上游水位和当地气温的实测数据如图6所示。

图6 测点PLA1径向位移与上游水位、当地气温的变化过程Fig. 6 Change process of measured displacement of measuring point PLA1 with upstream water level and local temperature

4.2 拱坝数值仿真模型与混合模型

利用ANSYS平台，建立该拱坝有限元仿真模型。模型中，地基模拟范围为坝踵向上游、坝址向下游以及建基面以下各延伸1.5倍坝高，网格单元41 030个，节点46 261个，拱坝有限元模型及坝体的网格划分如图7所示。模型的约束条件为底部边界、左右岸边界、上下游边界分别施加固定约束、横向约束和径向约束。

图7 拱坝有限元模型及坝体网格划分Fig. 7 Finite element model of the arch dam and grid division of dam body

结合有限元数值仿真模型，采用训练段水位测值对其进行加载，即可得到不同静水荷载作用下水压分量，以正常蓄水位（244 m）工况为例，有限元计算得到水压分量 δ1H(t)、 δ2H(t)与 δ3H(t)的位移云图如图8所示。以建模序列初始日水压分量作为初值，得到训练段不同静水荷载作用下的相对水压分量时间序列，通过回归分析，可得水压分量的表达式：

图8 拱坝正常蓄水位下3部分水压位移云图Fig. 8 Three parts of water pressure deformation cloud map of arch dam under normal storage level

基于水压分量数值模拟结果，结合式（6）可得PLA1测点水平位移的混合模型，模型回归系数见表1。建立PLA1测点水平变形统计模型，验证所建混合模型的有效性。测点PLA1径向位移拟合对比及建模残差序列如图9所示。

表1 回归模型系数Tab. 1 Coefficient of regression model

图9 测点PLA1径向位移逐步回归拟合对比与残差序列过程Fig. 9 Comparison of stepwise regression fitting of measuring point PLA1 radial displacement and residual sequence process line

由图9可知：混合模型和统计模型得到的拟合和预测变形与实测变形的变化规律一致，相比于统计模型建模结果，混合模型的拟合与预测变形更贴近实测值，且误差相对较小，表明所建混合模型有较好的预测性能；然而，混合模型与统计模型的残差序列均表现出一定的周期性特征，表明模型残差中仍有一定的有效成分尚待挖掘。

4.3 基于PSO-ELM和SARIMA残差分频修正的拱坝变形监控组合模型

采用EEMD对模型残差序列进行分解，合理挖掘混合模型残差中蕴含的有效成分。采用EEMD，添加白噪声的标准偏差设置为0.2，周期数设置为100，混合模型残差分解结果如图10所示。

图10 测点PLA1回归模型残差序列的EEMD分解结果Fig. 10 EEMD decomposition results of the residual sequence

由式（11）求得高低频信号的临界值为js=4，故IMF1～IMF5为高频信号，IMF6～IMF7以及余量r(t)为低频信号。对于分解得到的高频IMF信号，采用PSO-ELM模型对其逐一建模预报后叠加。其中，PSOELM预测模型中最大迭代次数Tmax、种群粒子数N、隐层节点个数L是影响预测模型的训练速度及精度的重要因素，为实现模型预测的最佳性能，确定预测模型的最优参数：隐层节点个数为30，种群粒子数为20，最大迭代次数为200，惯性权重w为0.9，加速常数c1和c2均设为1.8。同时，为验证对高频IMF信号建模的有效性，对IMF1～IMF5序列叠加得到的IMF（H）构建PSO-ELM预测模型。各高频预测模型的适应度曲线与预测结果如图11所示。

图11 高频预测模型的适应度曲线及预测结果Fig. 11 Fitness curves and prediction results of high frequency prediction model

由图（11）可知，对分解得到的高频IMF逐一建模预测，并将每个分量的最大误差与MAPE计算结果相叠加，得到最大误差和MAPE分别为0.089和14%，而未进行分频修正的高频分量预测结果的最大误差和MAPR分别为0.167和19%。说明分解得到的高频信号逐一建模预测更能有效挖掘序列内蕰的有效信息。

在利用SARIMA模型建模预报之前，需对IMF6、IMF7和余量r(t)进行平稳性检验与处理。经分析，3个低频分量经1阶差分后均转化为平稳序列。由于IMF6、IMF7经过1阶差分后仍具有明显的周期性，故选取周期S=12予以差分。差分后，结合AIC定阶方法确定相关参数取值，建立备选SARIMA模型，备选模型预测优度监测见表2。通过残差ACF、PACF图及各模型的AIC值、R2、MAPE值，比较选取最优模型。

表2 备选模型预测优度监测Tab. 2 Alternative model optimality detection

由表2可知：对于IMF6序列，当模型为SARIMA(2,1,0)(1,1,1)12时，MAPE最小，且R2最大，表明该模型拟合程度最优：绘制并观察该模型残差ACF、PACF图，如图12（a）、（b）所示。可以看出两个图形数值趋向于0，表明该模型残差序列不具有相关性，据此可知IMF6的最优模型为SARIMA(2,1,1)(1,1,1)12。同理，可分别确定IMF7和r(t)序列的最优模型分别为SARIMA(1,1,1)(1,1,0)12、ARIMA(0,1,1)(0,0,0)。将3组序列得到的预测值累加，得到低频信号的预测结果，如图12（c）所示。

图12 低频序列SARIMA模型预测结果Fig. 12 Prediction results of low frequency by SARIMA

将残差高频与低频项预测值叠加，得到混合模型预测值，即得到融合残差有效成分的混凝土拱坝位移预测混合模型。采用统计回归模型、混合模型与不考虑残差成分的EEMD-PSO-ELM模型进行预测，检验该组合模型的有效性，各模型的拟合与预报结果如图13所示。

图13 测点PLA1建模结果及残差对比Fig. 13 Displacement model prediction results at monitoring point PLA1

以位于拱冠梁坝段的正垂测点PLA2为参照点，采用本文方法预测其相应时段的水平位移验证所建模型的有效性与适用性，测点PLA2建模结果及残差对比如图14所示。同时，为直观地对3种模型的预测精度进行比较，分别量化计算4种模型预测结果的统计指标，统计指标见表3。

图14 测点PLA2建模结果及残差对比Fig. 14 Displacement model prediction results at monitoring point PLA2

表3 统计指标比较Tab. 3 Comparison of statistical indicators

由图13与14可知：采用4种模型构建的PLA1与PLA2测点水平位移监控模型，拟合与预测结果均与实测变形的变化规律一致，表明本文所建模型的有效性。与单一预报模型相比，本文所建组合混合模型的拟合与预测结果更接近实测位移，且拟合与预报误差更小，表明该模型较单一预报模型的预测精度更优。由此，验证了所建组合预报模型的合理性。由表3中的统计指标可知，相比回归模型，PLA1与PLA2测点水平位移混合模型R2相较于统计模型约提高了3%和2%，且MAE、MSE、MAPE也明显减小。由此表明，混合模型更能有效融合材料性能变化对拱坝变形的影响效应。此外，本文所建组合模型各项统计指标均优于其余3种单一模型。由此可知，融合多元数据挖掘手段，可有效地辨识原型监测信号中的时频非线性特征，并提取混合模型残差序列所蕴含的有效成分，进一步提升模型的预测能力。

5 结 论

1）考虑筑坝材料性能演变对变形预测模型的影响，结合有限元方法，构建混凝土拱坝变形预测混合模型；同时，针对混凝土坝变形混合模型残差序列呈现出的周期性特征，结合EEMD、ELM与SARIMA模型，有效挖掘了残差序列中所蕴含的有效信息成分。较单一模型而言，该模型可以捕捉到变形监测信号中绝大部分数据特征，有效弥补了常规模型信息挖掘手法单一的缺陷。

2）相比单一监控模型，考虑残差分频修正的组合预报模型有良好的拟合与预测精度，证明所提方法的有效性与合理性。本研究提出组合模型，可为水工施工安全监控模型中位移、渗流等效应的预测分析提供一种新方法。