特高拱坝运行初期变形预测模型温度因子选取方法

2021-02-02王春红马福恒

长江科学院院报 2021年1期

胡江，王春红，马福恒

（1.南京水利科学研究院水文水资源与水利工程科学国家重点实验室，南京 210029；2.南京瑞迪水利信息科技有限公司，南京 210029）

1 研究背景

我国已建、在建一批特高拱坝。这些特高拱坝技术指标多超出现行规范，且修建在高山峡谷，地质地形条件复杂［1］。工程实践表明，特高拱坝对温度边界和基础受力边界条件敏感，其变形特征与200 m以下的拱坝有着较大的差异，蓄水期和运行初期尤其明显［2-3］。

如当水库蓄水并长时间处于540～560 m高程时，溪洛渡上部谷幅收缩、拱坝整体向上游变形［4］。同时，坝体上游主要温度边界——坝前库水温垂直分层还未形成，水泥持续水化使得坝体内部温度回升，影响坝体变形［5-7］。据统计，在蓄水后5 a内发生失事的大坝占总失事大坝数一半以上［8］。变形是特高拱坝工作性态最直观的反映，开展其运行初期预测模型研究对保障工程安全有重要的意义［1-2］。

传统多元回归模型HST（Hydrostatic Season Time）是应用最广泛的大坝安全性态预测模型，但HST模型是针对运行期提出的，采用周期项温度因子表征坝体的稳定温度场，不适合复杂环境条件和非稳定温度场情形下的模型构建。对此，提出了一些解决方法，如基于实测温度因子的多元回归模型HTT（Hydrostatic Temperature Time）［9］，又如基于实测温度因子主成分分析（Principal Component Analysis，PCA）的多元回归模型HTPCAT（Hydrostatic-PCA based Temperature-Time）［10］。此外，过去20多年里非线性数学模型方法也得以大量应用，主要有神经网络模型、支持向量机（Support Vector Machine，SVM）等，两者均具有较高的拟合和预测精度［11］。

库水温对特高拱坝的应力应变和变形影响显著［6］，现行规范的库水温预测主要依赖适用于100 m级水库的经验公式，且未考虑蓄水和运行初期库水温变化的过程。狭长河道型水库坝前库水温垂向对流混接，分布复杂。运行初期变形预测模型须考虑库水温的时空变化［6-8］。同时，封拱后的坝体温度回升也将影响变形。考虑实测温度变量为特高拱坝运行初期变形监测和预报模型提供了思路［7-8］。因此，对众多监测仪器采集得到的海量数据中进行分析［12］，提取其中的关键信息是构建基于实测温度的预测模型需要解决的关键问题。

多变量时间序列是高维Euclidean空间，主成分分析（PCA）法将其近似转换为低维子空间，并保留了主要形态。另一种多变量压缩转换方法是分层聚类（Hierarchical Clustering，HC）［13］。当变量间距离均采用Euclidean距离，可将PCA作为HC的前处理步骤，降低干扰变量影响；还可采用PCA的映射图、因子图等直观地表示HC结果；可见，融合PCA和HC可更稳健地描述变量空间、更直观地表示聚类结果。

为此，本文融合PCA和HC法，提出基于主成分分层聚类法（Hierarchical Clustering on Principal Component，HCPC）的特高拱坝运行初期环境和坝体温度因子的选取方法，同时，考虑谷幅收缩和周期性的加载卸荷恢复变形，引入指数和周期项组合时效因子，在此基础上，构建HTT、SVM等2类模型，并将结果与HST模型进行对比，验证方法的科学性和合理性。

2 基于HCPC法的实测温度因子分类方法

2.1 基于PCA法的变量时间序列预处理

PCA法的核心思想是使用S个不相关变量即主成分描述数据集XK×I（K个变量、I次观测，S＜K）。这些主成分反映了原始数据集中蕴含信息的本质，其余的主成分则可认为为噪音［14-15］。因此，PCA也是一种去噪方法，可作为聚类分析的预处理步骤。

XK×I的元素为xkij，g是X的每列平均值的向量。与某一变量相关的特征值等于该变量的平方因子值的总和。一个变量对于某一主成分的重要性可通过该变量的平方因子值与这一主成分的特征值的比值即变量对该主成分的贡献率求得，变量k对主成分l的贡献率ck，l表示为［14-15］

式中：fk，l是变量k对主成分l的因子值；λl是第l个主成分的特征值。

变量贡献率介于0和1之间。对给定的某主成分，所有变量贡献率总和为1。贡献率越大，变量对主成分的贡献越大。一个重要的启发式的方法是，当某一变量贡献率大于所有变量的均值时，该变量可用来解释这一主成分。同时，具有高贡献率但符号不同的变量代表该主成分的2个维度方向，可用来解释该主成分［14-15］。

平方余弦表示一个主成分对某一变量的重要性，它表示了主成分对变量到原点平方距离的贡献，即

式中：xk，i为第k个变量、第i次观测值；gi为质心的第i个元素，即x的第i列平均值。

当数据居中时，式（3）简化为

在PCA映射图上，正相关变量聚集在一起，负相关变量则分布在某一主成分相反的2个方向；变量和原点间的距离表示了变量的重要性，远离原点的变量能很好地表示主成分。平方余弦值用于估计主成分表示变量的质量，高值表示质量好，这种情况下变量位于相关圆的圆周附近，对于解释主成分就越重要；反之，变量未完全由主成分表示，变量接近圆心、不具代表性。

2.2 基于HCPC法的温度因子选取

PCA的降噪和去冗余使得HCPC可获得比在原始变量上直接聚类更稳健的结果。HC法需要确定距离度量和分簇数量的标准。有许多距离度量（Manhattan、Euclidean等）和分簇数量确定（Ward、Single和Centroid等）方法［16-17］。这里，使用应用较广的基于Huygens定理的Ward法。

根据Huygens定理，总方差可以分解为组间方差和组内方差之和［16］，即

式中：xiqk为第q组的变量i在观测点k的值为各变量在观测点k的平均值为第q组内各变量在观测点k的平均值；Iq为第q组内的观测数；K为变量个数；Q为分组个数。

组内方差表征了组内变量的同质性，Ward准则在聚类时使每个步骤中组内方差增长最小，即组间方差减少最小。确定分组数是聚类分析的核心问题。分层聚类本质是一种嵌套分区，最底层上每个变量均是一个小组的，最顶层上所有变量都归属同一个大组。可根据组内方差的增长情况判断最优的聚类分组结果。当分组数从（Q-1）到Q时的组间方差的增加值ΔQ远大于从Q到（Q＋1）时的增加值Δ（Q＋1）时，最优分组数为Q组。

HCPC结果可表示在树形图或主成分映射图上。2种结果可视化方法相互补充、交叉验证。

3 因子选取和监测模型构建

针对特高拱坝布置了众多库水温、坝体混凝土温度测点的情况，基于PCA、HCPC结果，对运行初期温度因子测点进行聚类，并结合PCA映射图，确定各分区的高平方余弦值的变量为同一簇即同一分区的典型温度测点，将其作为变形预测模型的温度因子。

除温度因子外，运行初期的时效因子影响显著。变形时效包括混凝土徐变、刚固地基上坝体徐变、完整岩体外荷下徐变及岩体裂隙节理水压下时效位移等。其特点是运行初期变化急剧，而后随运行时间的增长逐步稳定。同时，混凝土及岩石在卸荷后徐变位移有一定的恢复［18］。

用一阶衰减微分方程表示运行初期随时间单调急剧增加的时效变形为

式中：c、r均是待定常数；t为时间。

当库水位呈年周期变化时，运行初期加载卸荷周期引起的徐变恢复δe也呈近似年周期变化［18］，其表达式为

式中：m＝1～2，这里取1；An、Bn均是待定常数。

当库水位呈年周期变化时，对比式（7）和传统模型中的周期温度因子，两者型式一致。

2）由于社会领域是造成城乡规划非理性化的主要原因，因此，大数据时代可通过对海量数据的高效利用，逐步实现城乡规划社会领域的计量化，为城乡规划中遇到的社会问题的计量处理提供参考信息。

通过平均绝对误差EMAE分析模型的拟合效果，其表达式为

式中：N是训练或者预测样本集的大小；yi是实测的监测效应量值；F（xi）是预测值。

水压则采用库水位的4次项表示；应用HCPC方法选取典型的库水温、坝体温度测点，作为变形预测模型的温度因子；同时，以式（6）、式（7）组合作为时效因子。将上述水压、温度、时效因子作为变量，构建HTT、SVM模型，以式（8）衡量模型的拟合和预测精度。

4 实例验证

4.1 工程和监测概况

某混凝土双曲拱坝坝顶高程610.00 m，建基面高程324.50 m，共31个坝段。河床坝基地质条件相对复杂，采用了“基础＋大坝”的扩大基础整体结构的方法（图1）。水库的正常蓄水位、死水位和汛限水位分别为600.00、540.00、560.00 m。水库具有年调节能力。

大坝内埋设了7 200支监测仪器以全方位地监测坝体的环境和效应量。在5＃、10＃、15＃、22＃、27＃坝段各布置1组垂线系统。垂线系统在高程610.00、563.25、527.25、470.25、395.25、347.25 m各布置正垂线测点，正垂线对应的坝基布置倒垂线。在典型高程上、下游坝面各布置1支温度计观测环境温度；坝体内部温度则采用位于坝体上游侧、中部和下游侧应变计测点的测温传感器获取。15＃坝段的垂线测点的布置和16＃坝段的温度测点布置如图1所示。

图1 15#坝段垂线测点和16#坝段温度测点示意图Fig.1 Layout of installed pendulums in Monolith 15 and temperature-measuring points in Monolith 16

2013年5月4 日，大坝导流底孔下闸，库水位从高程441.25 m开始逐步抬升。2014年3月6日，大坝完成最后一仓浇筑。2014年9月28日，水位第1次蓄至正常蓄水位600.00 m，蓄水和运行初期库水位过程如图2所示。坝顶垂线测点于2014年7月4日开始观测，15＃坝段垂线测值过程线如图2所示。以典型高程测点为例，16＃坝段的环境温度和混凝土内部温度的时空变化过程如图3所示，其中温度计测的是环境温度，应变计如S3-1测的是混凝土温度。以16＃坝段的库水温作为坝前库水温，2015年、2016年和2017年的典型时间的坝前库水温垂直分层如图4所示。

图2 库水位和15#坝段的垂线径向位移实测过程线Fig.2 Variations of reservoir water level and radial displacements of Monolith 15

图4 16#坝段库水温垂直分层情况示意图Fig.4 Vertical stratifications of reservoir water before Monolith 16

4.2 模型因子选取

由于坝顶测点初始观测日期为2014年7月4日，本次分析的垂线时间序列为2014年7月—2018年6月。在此时段内，坝体经历了4次完整的加载和卸荷周期。其中，库水位、温度为日均值，各1 460个测值；垂线测点的监测频次为2～3 d/次，各621个测值。

采用HCPC法，得到库水温各测点在第一、第二主成分维度上的贡献率及分簇，如图5所示，树状图和映射图的分区结果如图6所示，Ta表示气温。

图5 库水温测点的贡献图和分簇图Fig.5 Contributions and grouping results of reservoir water temperature points

HCPC法将库水温分为3类，且在映射图上与气温存在较大的差异。结合图3至图6可看出，480 m高程以下区域的库水温测点（T1—T17），库水温明显异于气温，但仍存在年变化，该区域水温年变幅小，且与气温存在较大的相位差；每年1—7月份温度持续处于较低的温度，在15℃左右。500.00～560.00 m高程区域的库水温测点（T21—T27），库水温变化基本类似，存在明显的年变化；蓄水后，该区域水温年变幅小于气温，且与气温存在相位差。

由于深孔底高程在490.00 m左右，并且每年均有泄洪，从而加强了上下层水流的热交换，使得该区域温度变化大于底部高程，T19测点既类似于低高程测点，又表现出明显的差异。560.00 m高程以上区域的库水温测点（T29、T31）受上游来水温度和气温双重影响，水温变化较为复杂。由图4可知，汛期垂向混合，温跃层高程降低，导致深水层水温上升；同时，来流密度大于表水层和深水层，进入库区后下潜，贴河床流动，直到遇到水温更低的深层水体，不再下潜，转为中间流，导致汛期垂向水温分布不规则［7］。总之，2014年以后，坝前垂向库水温的表温层、温跃层、滞温层层化结构明显，但汛期较为复杂。

图6 库水温测点的树状和映射图Fig.6 Hierarchical clustering and projection based on grouping results

同样地，采用HCPC法，得到PL15-5垂线段及其以下各温度测点在第一、第二主成分维度上的贡献率及分簇图，如图7所示，树状图和映射图的分区结果如图8所示。

限于篇幅，仅以PL15-5垂线段为例说明。结合图3、图4与图7、图8可知，395.25 m高程以下坝体温度整体趋于稳定，但各个部位有略有差异，HCPC法将该区域的环境和坝体温度测点分为6类。具体地，水化热引起的温度回升现象主要表现在坝体中部，测点主要为S6-4—S6-7、S3-2；坝体中部的S6-6测点又略有差异，仍呈较小幅度的温升现象；S6-9、S6-10则更为显著；坝体上游部分受库水温影响，已基本稳定，测点主要为T1、T3和S3-1；坝体下游部分呈较小幅度的温升，但受年周期的气温影响明显，测点主要为T2、T4和T6。从图1（b）可看出，扩大基础采用大坝混凝土，使得S6-4类似于坝体内部温度。

图7 PL15-5区间及其以下高程温度测点的贡献图和分簇图Fig.7 Contributions and grouping results of temperature-measuring points for PL15-5

图8 PL15-5垂线及其以下高程温度测点的树状和映射图Fig.8 Hierarchical clustering and projection based on grouping results of temperature measuring points for PL15-5

通过HCPC对其余垂线段温度测点分析可发现，PL15-4垂线段以上坝体仍受水化热影响，普遍呈温升趋势；坝体中部以单调温升为主；坝体上下游还受库水温、气温的影响，既表现出年周期变化，也呈温升趋势，但受环境温度影响的程度随高程的增加逐渐降低，如570 m高程左右，距上下游面4.2 m处的测点温度同坝体中部相同，近似呈单调温升状态；418 m高程以上都处于该阶段。

以16＃坝段的环境和坝体温度温度测点作为相邻15＃坝段的温度因子，采用HCPC法选取每一垂线段的典型温度计测点，结果列于表1。

表1 模型构建考虑的实测温度因子Table 1 Measured temperature points considered in model construction

最终考虑的所有环境、荷载和边界变量列于表2。一般地，选取样本总数的10%～20%作为预测集［19］，考虑到运行初期的环境、荷载条件的复杂性，这里选取15%作为预测集。

表2 模型构建考虑的因子Table 2 Variables considered in model construction

4.3 模型构建

采用HST、HTT、SVM这3类模型构建运行初期15＃坝段垂线径向位移的预测模型。HTT与SVM模型采用的变量相同，为表2中列举的所有因子；HST模型的时效为式（6）型式，其温度周期项同式（7），可见其与HTT模型的区别在于是否包含了实测的环境和温度变量。SVM采用了优化算法选择最佳的惩罚因子C、不敏感系数gamma分别为1 000和0.000 1。HST和HTT模型采用最小二乘法拟合。得到各模型的效果列于表3，将各变量标准化，HTT模型各变量拟合系数如表4。

表3 各模型效果对比Table 3 Comparison of performance parameters among constructed models

表4 HTT模型各变量的拟合系数Table 4 Fitted coefficients of variables of HTT models

以PL15-1测点为例，其拟合和预测效果如图9。分别统计库水位、库水温和气温等环境温度、坝体内部温度、时效等各影响因素对监测效应量的相对影响大小，结果如图10所示。

图9 PL15-1的3类模型的训练和预测效果Fig.9 Training and prediction performances of three constructed models for PL15-1

4.4 结果分析与讨论

从表3和图9可看出，HTT、SVM模型均具有较高的拟合和预测精度。HST模型精度较差，具有较大的残差和EMAE值，以PL15-1为例，测点的最大年变幅为2017年的34.18 mm，训练集和预测集的MAE值分别为5.96、6.62 mm，远超了规范要求的拱坝径向位移观测的中误差±2.0 mm的要求。SVM模型的预测精度较HTT模型要高，但HTT模型可解释各个变量的影响大小。利用kappa函数计算自变量矩阵的条件数，以PL15-1测点为例，计算得到的条件数为206.8，表明变量间存在中等相似共线性。

图10 各变量对PL15-1测点径向位移的相对影响Fig.10 Relative importance of variables to radial displacements of PL15-1

图10可知，在假设各变量相互独立的条件下，坝体温度（如S6-4、S6-9和S5-5测点）对15＃坝段的坝顶径向位移影响最显著，坝体温度的缓慢回升使得坝体产生了向上游变形［3］，为31.77%，这也是HST模型精度不足的重要原因；库水位波动引起的水压变化、时效及环境温度分别占24.63%、23.25%和20.35%。由此可见，特高拱坝运行初期变形受到加载卸荷、库水温垂直分层、坝体内部水化热温升、时效等多因素影响，造成了变形的非线性特征。