CN107220467B - 蓄水期库岸岩质边坡变形的预测方法 - Google Patents

Abstract

本发明公开一种蓄水期库岸岩质边坡变形的预测方法，包括建立包括各变形分量的库岸边坡变形回归分析预测模型，对边坡变形监测数据进行回归分析；基于预测模型，采用数值仿真分析方法对库岸边坡岩体力学参数和热学参数进行反演计算及分析；基于反演计算得到的库岸边坡岩体力学参数和热学参数，利用数值仿真分析方法计算分析当前水位下岸坡的变形，结合变形监测数据进行评估，根据评估结果调整参数，得到调整后的岩体力学参数和热学参数；根据调整后的岩体力学参数和热学参数，预测下一蓄水位下库岸岩质边坡的变形，得到下一蓄水位下库岸岩质边坡的变形预测结果；并对变形预测结果进行检验。基于本发明的混合预测模型，预测结果更为准确、全面。

Description

蓄水期库岸岩质边坡变形的预测方法

技术领域

本发明涉及一种蓄水期库岸岩质边坡变形的预测方法，属于水利水电工程技术领域。

背景技术

目前，水电工程中很多库岸岩质边坡在蓄水后变形启动，边坡虽未发生滑坡失稳破坏，但长时间处于蠕动变形状态，边坡产生了较大的时效变形，存在发生大型滑坡灾害的风险。如我国的锦屏一级、溪洛渡在蓄水后均表现出谷幅缩窄的规律，且在水位稳定后仍持续变形。我国拉西瓦水电站上游右岸果卜边坡蓄水后持续变形，截至目前累计最大变形已达40米左右，边坡整体保持稳定。此外，蓄水后库岸岩质边坡的变形不仅给边坡的稳定安全性带来不利影响，也会影响大坝的工作性态和长期安全状况，难以对大坝的稳定安全性做出科学评估。因此，建立库岸边坡蓄水期变形预测方法，掌握边坡变形规律，科学预测边坡变形破坏，对边坡的稳定与安全具有重要意义。

大型水利水电工程库岸边坡均埋设有变形监测仪器，积累有边坡长期变形监测数据资料，通过研究监测资料时间序列的变化规律，预测边坡的变形趋势。现有的边坡变形预测方法主要有：数值计算、传统的时间序列分析方法、灰色系统、神经网络方法、支持向量机以及遗传算法等，这些方法普遍存在不足之处，包括，第一，无法分离出由水压、渗流、气温变化等诱发的瞬时变形，以及受渗流、岩体蠕变等因素作用所产生的时效变形，没有明确将变形监测量分解为几个主要影响因素贡献的组合，因此，很难对库岸边坡变形的诱发原因和演化规律做出合理的物理力学解释；第二，基于测值资料和各种线性关系，对各监测点的“点”变形进行定量研究，能够反映边坡局部变形特性，由于监测点数量十分有限，其余部位的变形难以掌握，无法完整地反映边坡的工作特性、结构整体性及各部分之间的相互关系。

发明内容

鉴于上述原因，本发明的目的在于提供一种蓄水期库岸岩质边坡变形的预测方法，综合考虑水压、气温、时效(渗流、岩体蠕变)等影响因素，结合变形监测数据的回归分析和数值仿真建立边坡变形的混合预测模型，对边坡变形进行预测，预测结果较为准确、全面。

为实现上述目的，本发明采用以下技术方案：

本发明的优点是：

本发明的蓄水期库岸岩质边坡变形预测方法，结合变形监测数据建立由蓄水瞬时触发的变形分量、气温变形分量、时效变形分量组成的库岸边坡变形回归分析预测模型，基于该变形回归分析预测模型，采用数值仿真分析方法对库岸边坡岩体力学参数和热学参数进行反演，采用反演参数计算分析当前水位下岸坡的变形，结合变形监测数据评估参数结果，若参数合理，则进行下一蓄水位下库岸岩质边坡的变形预测。相对于传统的边坡变形预测方法，本发明提出的方法能更加科学合理地解释边坡变形的力学机制，能够对整个边坡的变形场进行反馈预测，预测结果较为科学准确。

附图说明

图1是本发明的方法流程图。

图2是本发明基于接触单元的法向-切向建立的局部坐标系示意图。

图3是本发明一具体实施例中K1监测点实测的变形监测数据统计线图。

图4是本发明一具体实施例中LS05监测点实测的变形监测数据统计线图。

图5是本发明一具体实施例中LS04监测点实测的变形监测数据统计线图。

图6是本发明一具体实施例中TP3-3监测点实测的变形监测数据统计线图。

图7是本发明一具体实施例中K1监测点的变形监测数据的回归分析结果线图。

图8是本发明一具体实施例中LS05监测点的变形监测数据的回归分析结果线图。

图9是本发明一具体实施例中LS04监测点的变形监测数据的回归分析结果线图。

图10是本发明一具体实施例中TP3-3监测点的变形监测数据的回归分析结果线图。

具体实施方式

以下结合附图和实施例对本发明作进一步详细的描述。

本发明公开的蓄水期库岸岩质边坡变形的预测方法，包括步骤：

S1：建立蓄水期考虑变形机制的库岸边坡变形回归分析预测模型，对边坡变形监测数据进行回归分析，分离出各影响因素作用下的变形分量；

在蓄水期，影响库岸岩质边坡变形的主要因素包括水压、气温、时效(渗流、岩体蠕变)等。通过分析各影响因素对库岸边坡变形的影响模式，建立蓄水期考虑变形机制的边坡变形回归分析预测模型：蓄水期库岸边坡监测点的实测变形，由蓄水瞬时触发的变形分量δ_HW、气温变化引起的变形分量δ_TW、时效变形分量δ_TD组成，即：

δ＝δ_HW+δ_TW+δ_TD (1)

以下对各变形分量进行详细说明：

1)蓄水瞬时触发的变形分量δ_HW

水库蓄水初期，作用于库盆表面的静水压力增大，库盆及两岸边坡产生变形。库盆压力使得库盆产生沉降变形，库水位以下河谷谷幅增大，库水位以上河谷谷幅收缩。此外，蓄水初期，水会迅速沿渗透性较强的断层、裂隙等构造面渗入岩体，岩体构造面上的有效应力减小，岩体产生变形。水头越高，边坡变形越大。因此，蓄水瞬时触发的变形分量的δ_HW可表示为水头的四次多项式，即：

其中，b_Hi为水头分量回归系数，H为水头。

2)气温变化引起的变形分量δ_TW

气温变化会对库岸边坡的变形产生影响，一般采取周期函数或多项式表示，即：

其中，t为位移观测日到起始监测日的累计天数；b_T1i、b_T2i均为温度分量回归系数。起始监测日是指监测仪器埋设后开始进行监测的日期(起点)，位移观测日是指开始监测以后获取监测数据的日期。

3)时效变形分量δ_TD

水库蓄水后，库水逐渐向两岸边坡岩体内入渗，库水入渗对岩体产生三方面的作用：一是物理化学作用，即岩体材料遇水后物理力学参数有所弱化；二是力学作用，即渗透力作用于岩体使得边坡的受力状态发生改变；三是温度作用，蓄水后库水逐渐向两岸岩体入渗，库水入渗后与岩体进行热量交换，由于库水温度比岩体温度低，从而导致岸坡岩体温度下降。三种作用会诱发库岸边坡产生变形。由于蓄水期库岸边坡的渗流场是由非恒定向恒定逐渐过渡、形成的过程，因此，渗流所引起的边坡变形具有时间效应。

另外，蠕变特性是岩体的重要力学性质之一。在外界条件(如，荷载、水文地质条件、材料物理力学特性等)改变的情况下，岩体的应力状态发生变化，诱发岩体产生变形，且岩体的应力应变状态、变形规律均随时间发生变化，并逐渐趋于稳定，具有明显的时效性。

因此，由渗流、蠕变等因素引起的库岸边坡时效变形分量可统一采用指数函数形式表示，即：

其中，n为蓄水过程中水位明显抬升和下降的次数；τ_i为位移观测日到第i次水位抬升时的累积天数；α_1i、σ_2i、β_i为第i次水位抬升或下降时的时效分量回归系数。蓄水过程中在水位抬升或下降时会打破岸坡原有的平衡状态，引起岩体应力调整重分布，在应力调整过程中会触发岩体产生时效变形。已有岩体时效变形的回归模型采用δ_TD＝α₁τ+α₂(1-e^-βτ)一项式，无法反映引起岸坡变形内在的物理力学机制。本发明提出的公式(5)与水位抬升或下降次数相关，能够反映蓄水过程中水位抬升或下降时岩体的应力的改变对岸坡时效变形的触发作用。

综上所述，蓄水期考虑变形机制的库岸边坡变形回归分析预测模型表示为：

利用该变形回归分析预测模型，对边坡变形监测数据进行回归分析，分离出各影响因素作用下的变形分量。

S2：基于变形回归分析预测模型，采用数值仿真分析方法对库岸边坡岩体力学参数和热学参数进行反演；

以公式(6)所示变形回归分析预测模型为基础，结合岩体设计力学参数和热学参数，采用黏弹塑性损伤本构模型对岩体进行模拟，进行岸坡岩体力学和热学参数反演。

1)岩体中结构面的本构模型；

自然界的岩体内部存在大量结构面，如断层、节理、层间错动面、软弱夹层等。这些结构面的存在导致岩体具有各向异性、非均匀性以及强非线性，在荷载作用下可能出现局部脱开、滑动、错位以及开闭等非连续变化现象，岩体的变形主要受结构面控制。因此，真实模拟岩体中结构面的物理力学特性，对岸坡岩体变形的准确模拟至关重要。

本发明采用有限元方法，基于有厚度、带强度的接触单元进行岩体结构面力学特性的模拟。如图2所示，以接触单元的法向-切向建立局部坐标系，在局部坐标系中，接触单元的应力-节点位移关系为：

其中，τ_s为切向剪应力，σ_n为法向正应力；h为接触单元的厚度；Δu为接触单元的剪切相对位移，Δv为接触单元的法向相对位移；{δ}^e为单元节点位移；[N]为形函数矩阵；[K]为接触单元的刚度矩阵，其由式(8)确定：

其中，k_ss为切向刚度；k_nn为法向刚度；k_sn为考虑切向耦合效应的刚度，k_ns为考虑法向耦合效应的刚度，由于目前很难确定耦合刚度项的数值，一般假定k_sn和k_ns均为零。

形函数矩阵[N]可由下式确定：

式中：

如图2所示，l为接触单元的长度，x′为局部坐标。

根据虚功原理，由式(7)可得到以下有限元方程：

其中，[K]^e为接触单元在局部坐标系的刚度矩阵，转换为全局坐标系的刚度矩阵为：

[K]^G＝[T]^-1[K]^e[T] (11)

其中，[K]^G为全局坐标系下的刚度矩阵，[T]为与结构面倾角θ相关的坐标转换矩阵。

当接触单元的法向应力达到其拉伸应力阈值，即σ_n≥f_t0时，接触单元产生拉伸损伤，法向刚度损伤变量D_n表示为：

其中：f_t0为接触单元的初始抗拉强度；f_tu为接触单元的拉伸损伤残余强度。

当接触单元产生拉伸损伤时，认为其切向刚度同时弱化，其损伤量与法向刚度相同，即：

D_n＝D_s (13)

D_s为切向刚度损伤变量

当接触单元的剪应力τ_s达到Mohr-Coulomb强度破坏准则时，即：

F＝|τ_s|-(c₀+f₀σ_n)≥0 (14)

接触单元产生剪切损伤，从而导致切向刚度弱化，相应的切向刚度损伤变量D_s表示为：

其中：τ_su为接触单元的剪切残余强度。当发生剪切损伤时，接触单元处于剪切滑移状态，其法向刚度保持不变。在数值仿真计算过程中，先判断接触单元的法向应力是否大于其抗拉强度，若不满足拉伸损伤破坏条件，再按照摩尔库仑准则判断其剪切应力是否达到其最大抗剪强度。

2)岩块的本构模型

边坡岩体中除了在宏观尺度上存在断层、节理、层间错动面、软弱夹层等结构面外，岩块内部还存在大量微裂隙。因此，本发明采用Drucker-Prager屈服准则模拟岩块的物理力学特性，即：

其中，α、k为常数；I₁为应力张量的第一不变量；J₂为应力偏量的第二不变量。

Drucker-Prager强度准则的破坏曲面为圆锥体，圆锥体的大小可通过α、k两个参数来调整。若圆锥面与莫尔受压子午线外接时，则有：

若圆锥面与莫尔受拉子午线相吻合，则有：

其中，

为岩块与结构面的内摩擦角，c为岩块与结构面的黏聚力。

3)岩体的蠕变模型

岩体是具有弹性、塑性和粘滞性的黏弹塑性体，即具有蠕变特性。岩体的蠕变模型可表示为：

其中，ε^c(t，τ)为蠕变应变，Δσ₀为在t₀时刻施加的应力增量，t为计算时间，τ为加载时间，A₁、A₂为蠕变度参数，k₁为蠕变速率参数。

4)反演分析方法

采用上述数值计算模型，结合岩体设计力学和热学参数，反演得到岩体当前的力学和热学参数。

以分离后蓄水瞬时触发的岸坡变形分量为基础，对岸坡岩体的变形模量和强度参数(抗拉强度、黏聚力、内摩擦角)进行反演，得到岸坡短期内的水头分量的力学参数。

以分离后气温周期性作用下变形分量为基础，对岸坡岩体的热学参数进行反演。

以分离后岸坡时效变形分量为基础，对岸坡岩体的时效变形参数进行反演，即对式(19)、(20)中的A₁、A₂、k₁等蠕变参数进行反演。

S3：基于反演计算得到的库岸边坡岩体力学参数和热学参数，以水库开始蓄水为计算起点，考虑实际蓄水过程和边界条件，利用数值仿真分析方法，计算得到当前水位下岸坡的变形场；

将变形计算结果与变形监测数据资料进行对比，分析二者之间的共同点和差异点，评价数值仿真计算结果的合理性和准确性，如果差异较大，需要调整岩体力学参数和热学参数，并复核边界条件重新进行计算。

S4：采用当前水位下调整后的岩体力学参数和热学参数，根据规划的未来短时间内水位变化情况，利用数值仿真分析方法计算得到下一蓄水位下库岸岩质边坡的变形场，预测边坡下一蓄水位时的工作性态，预评价下一蓄水位下岸坡的稳定安全性；

S5：对下一蓄水位下库岸岩质边坡的变形预测结果进行检验。

将下一蓄水位下库岸岩质边坡的变形预测结果与实测的新的水位下变形监测数据进行对比研究，分析预测结果的合理性和准确性。若预测结果与监测结果吻合良好，说明根据调整后的力学参数和热学参数取值合理，可直接进行下一阶段岸坡变形的仿真预测分析；若预测结果与监测结果相差较大，重新执行步骤S1至S3，反演计算得到新的力学参数和热学参数，重新进行预测。

以下结合一具体实施例对本发明的方法及效果进行详细说明。

国内西部某水电站上游库岸边坡，岸坡总体为花岗岩岩质，岸坡高700m，属于超高边坡。由于蓄水后岸坡3#山梁变形最为典型，且变形量较大，因此选取3#山梁典型剖面进行分析。该岸坡岩体由表及里分为四种结构类型，依次为：散体结构，碎裂结构，块裂结构和原岩。岸坡花岗岩体中断裂构造较发育，规模较大，性状较差，有四组结构面对岸坡形成过程及后期岸坡变形破坏起主要控制作用。前期勘测工作时发现岸坡上部存在一错落体，其底部控制面为规模较大的缓倾角断层HF104，断层上盘裂隙较密集，断层带内有明显蠕滑迹象，较为破碎，下盘岩体稍完整；错落体后缘控制面以顶部平台后缘断层LF1为界。

早期勘探结果表明，岸坡整体处于稳定状态，2009年3月初水库蓄水后，在库岸地质巡视过程中发现，边坡顶部平台有新的倾倒变形，岸坡前沿局部发生崩塌破坏。变形监测点布设四个：K1监测点、LS05监测点、LS04监测点及TP3-3监测点，前三个监测点自蓄水开始就有观测数据，TP3-3监测点安装时间较晚，2011年3月开始有监测数据。各监测点实测变形过程线如图3～6所示。监测数据资料表明：1)变形沿岸坡由下至上逐渐增大，岸坡顶部前缘变形最大，截止2015年1月，岸坡前缘K1监测点的最大综合变形为39.5m；2)边坡的变形与库水位有很强的相关性，蓄水是蓄水期岸坡产生新的倾倒变形的诱发因素；3)从观测到的变形过程线的趋势来看，早期库水位抬升时，岸坡测点变形速率在1～7天之后开始增大，且持续发展，当水位稳定后，变形速率又有所减小，目前变形并未完全收敛，还以1.8mm/d左右的变形速率变形。

根据公式(6)对各监测点的变形监测数据进行回归分析，由于该岸坡变形与库水位有很强的相关性，蓄水期的变形主要是由水位抬升触发的，气温对岸坡变形影响很小，因此回归分析中只考虑水位变化和时效变形分量。

如图7～10所示，根据各监测点实测的变形监测数据的回归分析结果，表明：1)采用本发明提出的回归分析方法得到的各监测点变形与实测的变形监测数据一致，且能够有效分离出各影响因素作用下的变形分量；2)通过对各监测点的变形分量进行分离，观测期岸坡的变形中由蓄水扰动直接引起的变形量很小，大多数为时效变形。

对岸坡3#山梁典型剖面进行有限元离散，基于各监测点的变形监测数据回归分析结果，结合岩体设计力学参数，进行力学参数反演，得到岸坡岩体及结构面的力学参数，如表1～表5所示。

表1不同分区岩体力学参数

表2不同分区结构面力学参数

表3断层力学参数

表4不同分区岩体蠕变参数

表5不同分区结构面蠕变参数

基于反演分析得到的岩体力学参数，考虑实际蓄水过程和边界条件，对岸坡蓄水期变形全过程进行数值仿真分析，对岸坡不同蓄水阶段的变形进行预测。根据绘制的变形云图，计算得到的2011年12月20日(水位2430m)时岸坡的最大综合变形约为32.16m，出现在岸坡顶部平台前缘，实测的综合变形监测数据为33.84m；计算得到的2012年4月20日(水位2440m)时岸坡的最大综合变形约为34.61m，出现在岸坡顶部平台前缘，实测的综合变形监测数据为35.19m；计算得到的2014年12月20日(水位2448m)时岸坡的最大综合变形约为40.94m，出现在岸坡顶部平台前缘，实测的综合变形监测数据为39.48m；通过仿真结果数据与实测的变形监测数据的对比可得，仿真计算得到的岸坡变形与实测值吻合较好。

根据仿真计算得到的岸坡变形形态结果，可得，岸坡中上部表层岩块有明显的倾倒变形，尤其是顶部平台前缘倾倒变形最为剧烈；岸坡顶部平台后缘有明显的下沉变形；由于缓倾角断层HF104的蠕滑作用，在缓倾角断层HF104出露附近有岩块挤出，甚至可能发生滚落。数值仿真得到的岸坡的整体变形形态与实际情况较为接近。

以上所述是本发明的较佳实施例及其所运用的技术原理，对于本领域的技术人员来说，在不背离本发明的精神和范围的情况下，任何基于本发明技术方案基础上的等效变换、简单替换等显而易见的改变，均属于本发明保护范围之内。

Claims (3)

1.蓄水期库岸岩质边坡变形的预测方法，其特征在于，包括步骤：

S1：建立包括各变形分量的库岸边坡变形回归分析预测模型，对边坡变形监测数据进行回归分析，该模型为：

δ＝δ_HW+δ_TW+δ_TD (1)

其中，δ_HW为蓄水瞬时触发的变形分量，δ_TW为气温变形分量，δ_TD为时效变形分量；δ为实测变形；

1)蓄水瞬时触发的变形分量δ_HW

蓄水瞬时触发的变形分量的δ_HW表示为水头的四次多项式，即：

其中，b_Hi为水头分量回归系数，H为水头；

2)气温变化引起的气温变形分量δ_TW

采取周期函数或多项式表示，即：

其中，t₁为位移观测日到起始监测日的累计天数；b_T1i、b_T2i均为温度分量回归系数；起始监测日是指监测仪器埋设后开始进行监测的日期，位移观测日是指开始监测以后获取监测数据的日期；

3)时效变形分量δ_TD

由渗流、蠕变因素引起的库岸边坡时效变形分量统一采用指数函数形式表示，即：

其中，n为蓄水过程中水位明显抬升和下降的次数；τ_i为位移观测日到第i次水位抬升时的累积天数；α_1i、α_2i、β_i为第i次水位抬升或下降时的时效变形分量回归系数；

S2：基于变形回归分析预测模型，采用数值仿真分析方法对库岸边坡岩体力学参数和热学参数进行反演计算及分析；

1)岩体中结构面的本构模型；

采用有限元方法，基于有厚度、带强度的接触单元进行岩体结构面力学特性的模拟；以接触单元的法向-切向建立局部坐标系，在局部坐标系中，接触单元的应力-节点位移关系为：

其中，{σ}为单元应力，τ_s为切向剪应力，σ_n为法向正应力；h为接触单元的厚度；Δu为接触单元的剪切相对位移，Δv为接触单元的法向相对位移；{δ}^e为单元节点位移；[N]为形函数矩阵；[K]为接触单元的刚度矩阵，其由式(7)确定：

其中，k_ss为切向刚度；k_nn为法向刚度；k_sn为考虑切向耦合效应的刚度，k_ns为考虑法向耦合效应的刚度，k_sn和k_ns均为零；

形函数矩阵[N]可由下式确定：

式中：

l为接触单元的长度，x′为局部坐标；

根据虚功原理，由式(6)可得到以下有限元方程：

其中，[K]^e为接触单元在局部坐标系的刚度矩阵，转换为全局坐标系的刚度矩阵为：

[K]^G＝[T]^-1[K]^e[T] (10)

其中，[K]^G为全局坐标系下的刚度矩阵，[T]为与结构面倾角θ相关的坐标转换矩阵；

当接触单元的法向正应力达到其初始抗拉强度，即σ_n≥f_t0时，接触单元产生拉伸损伤，法向刚度损伤变量D_n表示为：

其中：f_t0为接触单元的初始抗拉强度；f_tu为接触单元的拉伸损伤残余强度；

当接触单元产生拉伸损伤时，认为其切向刚度同时弱化，其损伤量与法向刚度相同，即：

D_n＝D_s (12)

D_s为切向刚度损伤变量

当接触单元的切向剪应力τ_s达到Mohr-Coulomb强度破坏准则时，即：

F＝|τ_s|-(c₀+f₀σ_n)≥0 (13)

接触单元产生剪切损伤，从而导致切向刚度弱化，相应的切向刚度损伤变量D_s表示为：

其中：τ_su为接触单元的剪切残余强度；

2)岩块的本构模型

采用Drucker-Prager屈服准则模拟岩块的物理力学特性，即：

其中，α、k为常数；I₁为应力张量的第一不变量；J₂为应力偏量的第二不变量；

3)岩体的蠕变模型

岩体的蠕变模型表示为：

其中，ε^c(t,τ)为蠕变应变，Δσ₀为在t₀时刻施加的应力增量，t为计算时间，τ为加载时间，A₁、A₂为蠕变度参数，k₁为蠕变速率参数；

S3：基于反演计算得到的库岸边坡岩体力学参数和热学参数，利用数值仿真分析方法计算分析当前水位下岸坡的变形，将计算得到的变形与监测结果进行对比，分析参数的合理性和准确性；如果计算值与监测值吻合良好，说明参数取值合理，可以进行下一水位下岸坡变形的预测分析，如果计算得到的变形与监测值吻合较差，需重新调整参数，得到岩体的力学参数和热学参数；

S4：根据调整后的岩体力学参数和热学参数，预测下一蓄水位下库岸岩质边坡的变形，得到下一蓄水位下库岸岩质边坡的变形预测结果；

S5：对下一蓄水位下库岸岩质边坡的变形预测结果进行检验。

2.根据权利要求1所述的蓄水期库岸岩质边坡变形的预测方法，其特征在于，所述步骤S2中，采用黏弹塑性损伤本构模型对岩体进行模拟，进行岸坡岩体力学参数和热学参数的反演计算。

3.根据权利要求1所述的蓄水期库岸岩质边坡变形的预测方法，其特征在于，所述步骤S2中，对边坡岩体力学参数和热学参数进行反演计算及分析的方法为：

对岸坡监测点的变形进行分离，以分离后蓄水瞬时触发的岸坡变形分量为基础，结合岩体设计力学参数，对岸坡岩体的变形模量和强度参数进行反演；

以分离后气温周期性作用下变形分量为基础，对岸坡岩体的热学参数进行反演；

以分离后岸坡时效变形分量为基础，对岸坡岩体的时效变形参数进行反演。

Images