基于元胞自动机模型的库岸边坡失稳模拟方法

2021-02-14黄梓莘

水资源与水工程学报 2021年6期

黄梓莘，陈波

(1.中国电建集团中南勘测设计研究院有限公司，湖南长沙 410014； 2.河海大学水利水电学院，江苏南京 210098)

1 研究背景

库岸边坡自身地质结构复杂，且外部荷载条件多变，其稳定性分析一直是水利工程治理和边坡失稳灾害评估的研究重点[1-2]。考虑到水库在运行过程中库岸边坡时刻受到外部环境的影响，导致其不可避免地呈现动态变化的特征，有学者引入直观可靠的实测资料，通过关注特征点位移变化分析库岸边坡的稳定性[3-5]。安可君等[6]结合监测资料分析重点岸坡的变形性态，总结了岸坡变形破坏的发展演化规律以及其对周边环境的影响。裴向军等[7]结合库水位变化曲线分析了岸坡的变形监测数据，发现了岸坡变形与库水位之间存在的强相关性。王国强等[8]基于InSAR技术获取大范围、高精度的库岸变形监测信息，开展了流域库岸的地质灾害研究。然而，水库蓄水打破了库岸边坡的自然演化过程，导致蓄水后库岸边坡呈现出与蓄水前不同的滑坡发育机制，仅仅凭借工程地质勘查资料或宏观变形监测资料，往往难以准确判断其运行稳定情况和未来发展趋势[9-10]。

针对复杂的库岸边坡失稳分析问题，有学者引入数值模拟方法加以研究。晏鄂川等[11]重点考虑了库水的动力作用因素，通过FLAC3D软件构建了滑坡的数值-力学模型。徐寅等[12]采用离散单元法研究了库岸边坡的动态破坏过程，为边坡风险评价提供了全新思路。张坤勇等[13]采用有限元强度折减法对边坡失稳的渐进破坏过程进行了解释和模拟。郑宇豪等[14]考虑库水位和降雨的影响，利用Geo-studio对滑坡进行失稳破坏概率计算。高道路等[15]结合Linear Cohesion接触模型研究了含水堆积体的运动特性，并采用离散元方法计算堆积体的失稳-滑动-堆积全过程。然而，目前的数值模拟方法大多将库岸边坡设定为均质体，而实际库岸高边坡在岩体风化、剥蚀、沉积等外力作用下往往形成上部为堆积体、下部为基岩体的二元结构边坡[16]，无论单独从岩质边坡还是单独从土质边坡的角度展开分析，均无法全面准确地描述库岸高边坡的运行性状，同时库岸边坡的物质组分具有强烈的随机性和不均匀性，导致实际模拟过程中的岩土力学参数取值也存在较大困难[17]。

因此，针对传统边坡失稳分析方法的岩土力学参数取值困难和动态变化特征考虑不足的弊端，本文提出一种基于元胞自动机模型的库岸边坡失稳模拟方法，考虑库岸边坡不同区域的物理力学性质，在变形实测资料驱动下动态辨识库岸边坡失稳路径，从而模拟库岸边坡发生失稳破坏的演化全过程。

2 研究方法

2.1 元胞自动机模型

元胞自动机是在离散动力系统中构造微观个体间的相互作用规则，以描述局部个体在强烈非线性作用下造成整体系统实现自组织演化的过程[18]，其模型具有建模思路简单、模拟功能强大、约束条件宽松等优点。现将元胞自动机模型的主要构成元素介绍如下：

(1)元胞：元胞通常以规则散落在元胞空间的网格形式表示，是模型中的最基础的组成部分。

(2)元胞状态：元胞状态多以布尔型变量进行表示，即采用0或1分别代表元胞发生破坏和元胞维持完好的两种对立状态。

(3)元胞空间：元胞空间是指元胞在空间中分布的网格集合，以二维元胞空间为例，主要的网格划分形式如图1所示。

图1 元胞空间主要网格划分形式

(4)元胞邻居：元胞邻居直接影响了元胞自动机模型模拟的准确度和复杂性，二维空间局部范围内的元胞邻居模型如图2所示，图中黑色网格均代表主元胞，而灰色网格则代表元胞邻居。

图2 二维空间局部范围内的元胞邻居模型

(5)元胞时间：元胞时间是一个概念抽象的离散尺度，是一系列无量纲的整数值，通常假设当前时刻为t，时间步长Δt为1，则下一时刻表示为t+1。

(6)演化规则：根据当前时刻元胞状态以及邻近范围的元胞邻居状态，共同确定下一时刻元胞状态的函数表达式f，即元胞自动机的演化规则，存在：

S1t+1=f(S1t,SNt+1)

(1)

式中：S1t为主元胞t时刻的状态;SNt+1为元胞邻居t+1时刻的状态。

(7)边界条件：理论上元胞空间是呈无限延伸的状态，但由于目前的计算模拟条件尚无法有效处理无限网格，因此需要对模型设置边界条件。

2.2 边坡失稳模拟的元胞自动机模型

对于库岸边坡而言，在空间尺度上，基本单元之间存在相互作用的关系，基本单元的状态变化会受到邻近单元状态的影响；在时间尺度上，基本单元的状态变化则会受到前一时刻状态的影响，这种时空局部性与元胞自动机模型的基本原理相契合，故本文引入元胞自动机模型模拟库岸边坡发生失稳破坏的演化过程。现结合图3所示的边坡失稳元胞自动机模型流程图，将方法主要步骤介绍如下：

图3 边坡失稳元胞自动机模型流程图

(1)构建元胞空间。选定边坡典型断面作为分析对象，将三维边坡简化为二维平面，然后根据边坡的形态尺寸，将边坡离散成规则的正方形网格，设定包括邻居型式、时间步长、边界条件、能量传递系数等模型参数，构建用于计算模拟的元胞空间。其中，元胞i的坐标为(xi,yi)，i=1,2,…,M，M为元胞空间内的元胞总数。考虑到边坡能量传递方向具有随机性，采用图2(b)所示的摩尔型邻居型式，以元胞i为中心，定义局部范围内的元胞为元胞邻居j(j=1,2,…,8)。

(2)赋予初始能量。边坡在形成过程中存储有一定的弹性应变能，因此在边坡失稳的元胞自动机模型中需要考虑元胞初始能量E(xi,yi, 0)。考虑到边坡监测信息是边坡运行过程中最直观的表现形式，其变化规律可以反映边坡自身的结构性态及外部环境的荷载特性和变化规律，此处结合边坡工程的实时监测资料对各个边坡元胞进行初始能量赋值，其典型力-变形-能量曲线如图4所示[19]。

图4 边坡元胞典型力-变形-能量关系曲线

将位移测点映射到边坡典型断面，选定同一时间断面的监测资料对不同位移测点位进行赋值。在根据边坡断面轮廓尺寸确定的边界条件约束下，利用克里格法对输入的测点监测信息进行插值处理，形成监测信息的插值网格δi(xi,yi)，i= 1, 2, …,M，分别计算各个元胞应变能作为边坡失稳模拟的初始能量，其计算式为：

(2)

式中：Ui为元胞i的应变能;δi为元胞i的变形值;σi为元胞i的应力值。

对图4中曲线三角形OAB进行简化考虑，近似计算为U(xi,yi, 0)=0.5·Eiδi2，其中Ei为元胞i的杨氏模量，由此形成元胞空间的初始状态分布。

(3)设置破坏阈值。由于边坡的岩土材料具有各向异性，同时堆积体区和岩体区的材料参数也存在显著差异，因此在考虑边坡物理力学性质的情况下，此处采用两参数的Weibull概率分布对边坡的岩土材料参数进行分区模拟[20]，设置各个边坡元胞的破坏阈值：

(3)

式中：x为边坡岩土的材料参数；μi为元胞i的弹模特征值；mi为元胞i的材料均质度；λ为破坏系数。

(4)确定失稳起始位置。搜索监测信息插值网格确定边坡位移测值的极大值，将该极值点作为模拟边坡失稳破坏的起始位置，并将起始位置所对应的元胞状态调整为0，显示为黑色网格，代表元胞发生破坏；其余元胞状态仍设置为1，显示为白色网格，代表元胞维持完好。

(5)设定演化规则。在元胞空间中，所有单位元胞按照一定的演化规则进行状态转移，这也是边坡失稳路径的发展条件。图5为边坡元胞空间演化规则示意图。

结合图5所示对元胞空间的演化规则具体介绍如下：

① 能量积累。当元胞存储能量E(xi,yi,t)小于破坏阈值Ecr(xi,yi)时，元胞维持完好，元胞状态为1，此时在外部荷载作用下，元胞将持续积累能量，则元胞i在时刻t的能量E(xi,yi,t)表达式为：

E(xi,yi,t)=E(xi,yi,0)+e(t)

(4)

式中：E(xi,yi,t)为t时刻元胞i的能量，kJ；E(xi,yi,0)为元胞i的初始能量，kJ；e(t)为外界环境的输入能量，kJ，此处采用e(t)=0.05t·Rnd模拟能量输入方式，其中Rnd为范围在[0, 1]的随机数，表征边坡元胞能量积累的非均匀性质。

②能量传递。当元胞存储能量E(xi,yi,t)大于破坏阈值Ecr(xi,yi)时，元胞发生破坏，状态由1调整成0。根据能量守恒规则，破坏元胞i所存储的能量中有一部分被耗散，另一部分按能量传递概率向元胞邻居j传递，此时元胞邻居j在下一时刻的存储能量演化为：

E(xi,yi,t)=Ecr(xi,yi)→E(xi,yi,t)=0

(5)

ΔE(xi,yi)=(1-α)E(xi,yit,t)

(6)

E(xj,yj,t)→E(xj,yj,t+1)=E(xj,yj,t)+

w(i,j)ΔE(xi,yi)

(7)

式中：E(xj,yj,t+1)为t+1时刻元胞邻居j的能量积累，kJ；E(xj,yj,t)为t时刻元胞邻居j的能量积累，kJ；α为能量耗散系数，此处根据元胞邻居的破坏个数将一定能量作为耗散能量处理； ΔE(xi,yi)为破坏元胞i所传递的能量，kJ；w(i,j)为元胞i与元胞邻居j之间的能量传递概率，同时表示元胞邻居吸收能量的能力，故取值参考各元胞邻居的破坏阈值占所有元胞邻居阈值之和的比例，即：

(8)

③发生破坏。在发生能量传递后，若每个元胞的存储能量均满足E(xi,yi,t)≤Ecr(xi,yi)，则元胞维持完好状态，并根据公式(4)再次获取能量；若某个元胞的存储能量满足E(xi,yi,t)>Ecr(xi,yi)，则元胞发生破坏，根据公式(5)～(7)发生能量传递。

(6)输出模拟结果。重复第(4)步的演化过程，分别记录元胞空间在不同时刻的破坏情况，即可模拟边坡发生失稳破坏的演化过程，以达到动态辨识边坡失稳路径的目的，当元胞空间不再发生能量传递，或元胞破坏数目最终恒定时，即代表演化过程停止，判定模型计算结束，并输出最终的边坡失稳区域分布。同时，为了直观展现边坡元胞空间的破坏情况，此处引入边坡失稳破坏指标I(t)=N/M(N为t时刻发生失稳破坏的元胞个数)，从而可定量计算不同时刻破坏元胞面积占整个边坡元胞空间面积的比例。

图5 边坡元胞空间演化规则示意图

3 工程实例模拟与分析

3.1 工程概况

本文所选取的研究对象为某水库左岸边坡一古老滑坡堆积区，该边坡上覆岩体在河谷不断下切过程中，受重力作用曾多次出现顺层滑动的现象，且目前仍然存在失稳的可能性，情况严重时甚至可能危及整个上部边坡堆积体的稳定性。该水库岸边坡滑坡体(I区)的位置及工程地质剖面图分别如图6(a)、6(b)所示。

3.2 元胞自动机模型参数设置

采用元胞自动机模型对所选库岸边坡进行失稳路径辨识和失稳破坏模拟，首先对边坡失稳元胞自动机模型的各项参数进行设置：

(1)元胞空间：在图6(b)所示的边坡工程地址剖面图的基础上，根据边坡形态尺寸将边坡典型断面离散成单位长度为1 m的正方形网格，所构建的元胞空间为900×460的元胞列阵，如图7所示。

(2)边界条件：选取断面轮廓作为定值型边界，在失稳路径辨识图中以黑色网格表示。

(3)元胞状态：各个边坡元胞被赋予两种基本状态{0,1}，其中0表示元胞发生破坏，网格以黑色网格表示；1表示元胞维持完好状态，网格以白色网格表示。

(4)邻居型式：采用标准摩尔型定义元胞的邻居型式，并设定邻居半径为1，即定义与主元胞邻近的8个元胞为元胞邻居。

(5)元胞破坏阈值：工程实例边坡在河谷不断下切的过程中暴露出假整合面，上覆岩体在重力作用下发生多次顺层滑动，因此在考虑边坡地质构造和演化历史的情况下，参考《地基与基础》[21]分别选取堆积体和岩体的岩土材料参数，并采用两参数的Weibull概率分布计算各边坡元胞的破坏阈值，如图8所示。

图6 工程实例的库岸边坡滑坡体位置及地质剖面图

图7 边坡典型断面的元胞空间示意图图8 边坡元胞破坏阈值分布图

(6)元胞初始能量：将滑坡体(Ⅰ区)位移测点投影到图6(b)所示的边坡典型断面上，基于各测点2020年全年的位移监测数据，利用克里格法对输入监测信息进行插值处理，初步形成图9所示的监测信息插值网格，以此为基础分别计算各个元胞的应变能作为边坡失稳模拟的初始能量，进一步形成元胞空间的初始状态分布。

图9 克里格插值结果及边坡失稳模拟起始位置

(7)失稳起始位置：为模拟边坡发生失稳破坏的变形演化过程，此处通过搜索克里格插值结果中的极大值，假设边坡失稳元胞自动机模型的起始位置，选取的边坡失稳模拟起始位置如图9所示。

汇总边坡失稳的元胞自动机模型参数如表1所示，由此构建边坡失稳元胞自动机模型的元胞空间。

表1 边坡失稳的元胞自动机模型参数

3.3 边坡失稳破坏模拟结果

在设定模型参数后，构建边坡失稳元胞自动机模型的元胞空间，按照设定的演化规律进行元胞空间的能量传递，并随时间步长增加记录不同时刻的元胞破坏状态。图10展示了t在10、50、55、60、70和90时步时，采用边坡失稳元胞自动机模型动态辨识边坡失稳路径的模拟结果，图11展示了边坡失稳破坏指标I(t)随时步长增加的判别结果。