基于局部耦合马尔科夫链模型的钻孔优化方法

曾正强; 蔡永昌; 吴江斌

doi:10.11779/CJGE20230927

基于局部耦合马尔科夫链模型的钻孔优化方法 English Version

曾正强^{1, 2,},
蔡永昌^{1, 2},
吴江斌³

1.
同济大学土木工程防灾国家重点实验室, 上海 200092
2.
同济大学土木工程学院, 上海 200092
3.
华东建筑设计研究院有限公司, 上海 200002

基金项目:

国家自然科学基金高铁联合基金重点支持项目 U1934212

详细信息

作者简介:
曾正强(1994—)，男，四川内江人，博士研究生，主要从事岩土工程等领域科研工作。E-mail：zzqiang@tongji.edu.cn

中图分类号: TU42
计量
- 文章访问数: 224
- HTML全文浏览量: 21
- PDF下载量: 60
出版历程
- 收稿日期: 2023-09-19
- 网络出版日期: 2024-04-18
- 刊出日期: 2024-11-30

Borehole optimization method utilizing local coupled Markov chain model

1.
State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University, Shanghai 200092, China
2.
College of Civil Engineering, Tongji University, Shanghai 200092, China
3.
East China Architectural Design & Research Institute Co., Ltd., Shanghai 200002, China

摘要

摘要: 如何合理利用勘查资料进行准确的地质建模和逐级优化钻孔方案，对于降低勘查成本和提高地层数据收集效率具有十分重要的意义。由于岩土勘测初期钻孔数据通常十分稀少，而工程区域的地质结构可能随位置发生显著变化，这大大增加了钻孔设计的难度。CMC模型可以处理稀疏地质勘查数据，并且具有参数简单和理论简洁的优点，因此被广泛用在地质建模和地质不确定分析中。然而，已有的CMC模型无法用于复杂变化地层的钻孔设计，因为它们依赖全局的序列方向和Walther常数，可能导致对地质结构的错误估计。提出了一种基于局部耦合马尔科夫链模型（Local coupled Markov chain model，简称LCMC模型）进行额外钻孔优化的方法。该方法利用地质剖面钻孔数据的片段化处理、局部随机建模和多片段叠加的方式建立适应复杂变化地层的地质模型；使用基于LCMC模型生成的信息熵图定量评价地质单元的不确定性；基于柱平均信息熵曲线对额外钻孔的最佳位置进行逐步预测。研究结果表明，相比于传统方法，提出的方法可以实现更加合理高效的钻孔优化，增加倾角和方向发生变化的复杂地层的建模准确性，以及降低地质剖面模拟的不确定性。提出的方法为复杂变化地质区域的钻孔设计提供了一种有益的参考。
- 局部化 /
- CMC模型 /
- Walther常数 /
- 钻孔优化 /
- 复杂地层
Abstract: The rational use of exploration data for accurate geological modeling and the step-by-step optimization of drilling plans are of great importance for reducing the exploration costs and improving the efficiency of stratigraphic data collection. Due to the typically sparse borehole data at the early stage of geotechnical surveying and the potential for significant variations in the geological structure across different locations within the project area, the difficulty of borehole design is greatly increased. The CMC model can handle sparse geological exploration data and is widely used in geological modeling and geological uncertainty analysis due to its simple parameters and straightforward theory. However, the existing CMC models cannot be used for the borehole design in complexly varying strata because they rely on the global sequence direction and the Walther's constant, which may lead to incorrect estimates of geological structures. A method for optimizing the additional boreholes is proposed based on the local coupled Markov chain model (LCMC model). This method establishes a geological model that adapts to complex and varying strata by fragmenting borehole data from geological profiles, employing local stochastic modeling, and overlaying multiple fragments; it uses an entropy map generated by the LCMC model to quantitatively evaluate the uncertainty of geological units; and it progressively predicts the optimal locations for the additional boreholes based on the column-averaged entropy curve. The research results show that, compared to the traditional methods, the proposed approach can achieve more reasonable and efficient borehole optimization, increase the modeling accuracy of complex strata with changing dips and directions, and reduce the uncertainty in geological profile simulations. The proposed method offers a the beneficial reference for borehole design in geologically complex and variable regions.
- localization /
- CMC model /
- Walther's constant /
- borehole optimization /
- intricate stratigraphy

HTML全文

0. 引言

边坡是人类活动中极为重要的自然地质环境之一。如何保证边坡稳定成为了岩土工程一个重要的研究领域^[1]。在进行边坡稳定分析时，最主要的判断依据就是稳定安全系数。根据安全系数的大小来评判边坡是否稳定，这就涉及安全系数的定义问题。由于各种分析方法所依据的安全系数定义的不同^[2]，对于同一个计算对象，其计算出来的安全系数值一般也有所不同。但现有的一些安全系数的定义式中，很少考虑到滑动的矢量特征，且对边坡真实状态引入了过多的人为假定。为了克服这方面的局限性，葛修润院士从滑动的基本概念出发，指出滑动本身就是一个矢量变化的过程，并从矢量的角度入手，可求出滑裂面的下滑力与抗滑力，从而进一步求出边坡的安全系数。这一求解方法，称之为矢量和分析方法（vector sum method, VSM）^[3-4]。矢量和分析方法是基于力是矢量这一基本概念，安全系数定义为坡体在受各种荷载（包括自重）作用下，潜在滑面上提供的抵抗坡体滑动的合力矢与作用在潜在滑动面上贡献滑动力的合力矢在滑体的整体滑动趋势方向上投影的比值。该方法具有明确的物理意义，考虑了滑动方向，能较好地反映边坡的整体稳定性。

矢量和分析方法是根据力的矢量特性提出的一种边坡稳定性分析方法，该方法在求到了滑体的合力后，通过假定主滑方向，认为滑体沿主滑动方向平面滑动，目前常用如下4种方法来确定主滑方向：①基于潜在滑面各点真实剪应力合力矢单位方向；②基于极限下滑力的合力矢方向；③基于边坡潜在滑面剪出点指向剪入点的位置矢量方向；④基于最小势能原理。然而，在滑坡发生时，滑动面通常呈现弧形、平面、折线或曲面，这意味着坡体滑动不仅涉及平动问题，还需考虑坡体的转动，因此本文针对坡体转动提出了一种新的边坡稳定性分析方法：力矩和法（torque sum method, TSM）。此方法定义边坡稳定安全系数F为：坡体在受各种荷载（包括自重）作用条件下，潜在滑动面所能提供的抵抗坡体绕某点转动的极限强度所提供的力矩之和与滑动体所受各种荷载（包括自重）所引发坡体绕这点转动的力矩之和的比值。力矩和法可以直接采用显式计算，可以对任意滑动面上的安全系数进行求解计算，具有清晰的物理含义，概念清楚，可以直接根据滑动面的实际应力分布以及滑动面所能提供的抗滑能力进行计算，弥补了矢量和方法不能考虑坡体可能发生转动失稳的不足。本文给出了力矩和法的定义表达式，适用于一般的边坡岩土体力学分析。运用本方法计算了两个典型的边坡算例，并与传统分析方法的计算结果进行了对比分析。

1. 基于力矩和法的安全系数分析方法

1.1 安全系数的定义

对于可能发生坡体转动的情况，参照矢量和方法，边坡稳定安全系数F可以通过计算坡体应力状态获得。力矩和方法的安全系数F定义为潜在滑动面所能提供的绕某转动中心的极限抗转动力矩的总和∑M_r与作用在滑动面上计算应力绕转动中心的滑动力矩的总和∑M_s之比，即

$F=\frac{\sum \boldsymbol{M}_{\text{r}}}{\sum \boldsymbol{M}_{\text{s}}} 。$

(1)

而极限抗转动力矩是由滑动面上的极限应力强度提供。由于力是一个矢量，不但有作用的大小，而且有作用的方向，所以式(1)中求和是考虑了力的矢量性质，遵循矢量求和运算法则合成的一个力矩和。

1.2 基本假定

对于一种稳定性分析方法，若要将其所涉及到的所有因素都考虑进去，问题将会是复杂和繁琐的^[5]。如果通过适当的假设，可以让分析方法更加的有针对性，并把握住问题的主要矛盾，进而降低问题的复杂性。对力矩和法，将建立在如下假定的基础之上：

（1）对某一特定的边坡进行稳定性分析的先决条件，是通过地质勘察，确定边坡的基本物理力学参数和边坡的载荷以及边坡的力学边界条件，而坡体内应力分布状态可由其它应力分析方法（如有限元法、离散元法、条分法等）计算获得。

（2）在二维边坡问题中，确定边坡所构成的计算范围为 $\mathit{Ω}$ ，假设在该计算范围内已知边坡的潜在滑动面L，构成滑动体的区域为S。

（3）在边坡潜在滑动面上，假设岩土体的材料强度特性服从Mohr-Coulomb准则，抗剪强度为

${\tau _{\text{f}}} = c + \sigma \tan \varphi 。$

(2)

（4）对绕某点转动的力矩和法安全系数F(x_p, y_p)定义为：转动点p坐标为(x_p, y_p)，滑动面上所提供抵抗坡体转动的各力沿转动点的力矩代数和∑M_r(x_p, y_p)与滑动面上计算应力沿转动点的力矩代数和∑M_s(x_p, y_p)的比值：

$F({x_{\text{p}}},{y_{\text{p}}}) = \frac{{\sum {{M_{\text{r}}}({x_{\text{p}}},{y_{\text{p}}})} }}{{\sum {{M_{\text{s}}}({x_{\text{p}}},{y_{\text{p}}})} }} 。$

(3)

（5）滑动面是滑动体与母体脱离发生绕某点转动的面，可以是任意形状。对于滑动（转动）的滑体，滑动面是约束边界。做力矩和计算时，滑动体假定为刚体。

1.3 滑动体转动中心的确定

当分析某一具体的潜在滑动体时，也即确定了其对应的潜在滑动面。作为平面刚体，刚体最基本的运动形式为沿某直线的平动和绕某点的转动，如果不考虑其后的运动轨迹，从静止时的运动趋势考虑，可以认为刚体只有平动和转动，而且其启动加速度在理论上是可以计算的。

对于边坡滑体，以往的知识认为滑体是沿滑动面滑动（或转动），笔者认为滑动面是边坡滑体失稳后，产生运动形成的轨迹，与边坡初始失稳时形成的失稳面（或破坏面）是有区别的，如果将滑动面的形式与位移或运动的形式相关联，则无法解释现有的分析滑动面为非圆弧或非直线的方法中，滑动体作为刚体时如何发生。当边坡初始失稳时，无论滑动面形式如何，滑动面只是滑动体的一个边界，其相互作用力可以认为是边界上的力边界条件。

因为刚体的基本位移是平动和转动，分析滑动体时，应该按照其基本的位移特征进行分析。矢量和方法考虑了滑动体的平动位移，认为滑动体将沿某一直线运动，并通过力矢量求和等相应的假设，获得其安全系数。作为对矢量和方法的一个补充，本文将只考虑滑动体的转动。对于同时考虑平动和转动的复杂情况，将在往后的工作中开展。

值得强调的是，这里所指的滑动体的平动和转动，均与滑动面无关，仅与滑动体本身和其受力有关，平动和转动也只是指滑动体在初始失稳时的位移趋势，与滑动体失稳后的位移和滑动面的形成无直接关系。实质上，边坡的失稳是指的滑动体在初始位移趋势上的失稳，并在考虑初始位移趋势上计算其失稳安全系数。

如何获得滑动体的初始位移趋势或初始转动中心，可以按照边坡稳定性分析的普遍方法，也即求解边坡最小安全系数的方法来获得。定义其初始位移趋势（或转动中心）是使滑动体具有最小安全系数时对应的转动中心。

仅考虑滑动体转动时，如图 1为边坡稳定性分析时的安全系数求解示意图。中可以通过有限元法或其它方法，计算得到潜在滑动面上任何一点i的正应力 ${\sigma _i}$ 和切应力 ${\tau _i}$ 。设滑动面上点i处岩土体内聚力为 ${c_i}$ ，内摩擦角为 ${\varphi _i}$ ，中 ${\alpha _i}$ 为滑动面在i点上的切线相对于坐标系x轴的正向夹角。

图 1 转动点搜索示意图

Figure 1. Schematic diagram of rotating point searching

下载: 全尺寸图片幻灯片

设点p(x_p, y_p)为转动点，为了寻找所计算滑动体具有最小安全系数的转动点，本文采用了较直观的网格（Grid）搜索法，转动点假设位于设定的一矩形搜索区域，矩形区域均布网格节点，作为试算的转动点。每一节点只需计算一次安全系数，就可通过比较获得区域内具有安全系数最小的节点。如果获得的节点位于矩形区域边界上，则需调整区域。该方法与圆弧滑动面用Grid方法搜索滑动圆心类似（图 2）。

图 2 安全系数求解示意图

Figure 2. Solution of factor of safety

下载: 全尺寸图片幻灯片

1.4 力矩和法的安全系数计算

力矩和法安全系数的计算包括转动力矩和计算与抗转动力矩和计算两方面：

（1）转动力矩和计算

滑动体上任一微段 ${{\Delta }}{l_i}$ 的正应力 ${\sigma _i}$ 和切应力 ${\tau _i}$ 是由滑体自重和作用在滑体上的外载荷所产生的，理论上滑动面（也是边界）上的所有正应力 ${\sigma _i}$ 和切应力 ${\tau _i}$ 与滑动体的自重和外荷载是一个静力平衡系。正应力 ${\sigma _i}$ 和切应力 ${\tau _i}$ 对转动中心求力矩并沿滑动面求和，就是滑体的转动力矩和。同时，岩床反作用在滑体上的力为 ${\sigma '_i}$ 和 ${\tau '_i}$ ，它们数值相等，方向相反。如果能够求到滑动体的重心，转动力矩和的求解可以直接求重心对转动中心的力矩与滑体上的外荷载对转动中心的力矩的总和。这与沿滑动面由 ${\sigma _i}$ 和 ${\tau _i}$ 求力矩和是一致的。

将滑动面上的应力对转动点力矩沿滑动面积分（或求和）即为式（3）中的分母项。

滑动面上任何一点i的切应力 ${\tau _i}$ 的斜率为 ${k_{\tau i}}$ ，所在直线方程为

$y={k}_{\tau i}x+{b}_{i} 。$

(4)

式中： ${k_{\tau i}} = \tan {\alpha _i}$ ， ${b_i} = {y_i} - {k_{\tau i}}{x_i}$ ，(x_i, y_i)为i点坐标。

转动点(x_p, y_p)到 ${\tau _i}$ 所在直线的距离为

${d}_{\tau i}=\frac{\mid {k}_{\tau i}{x}_{\text{p}}+\text{(}-1\text{)}{y}_{\text{p}}+{b}_{i}\mid }{\sqrt{{k}_{\tau i}{}^{2}+{\text{(}-1\text{)}}^{2}}} 。$

(5)

切应力 ${\tau _i}$ 提供转动力矩，力矩和为 $\sum _{i=1}^{n}{\tau }_{i}\text{Δ}{l}_{i}{d}_{\tau i}$ 。

滑动面上任何一点i的正应力 ${\sigma _i}$ 的斜率为 ${k_{\sigma i}}$ ，所在直线方程为

$y={k}_{\sigma i}x+{B}_{i} 。$

(6)

其中， ${k_{\sigma i}} = - \frac{1}{{{k_{\tau i}}}} = - \frac{1}{{\tan {\alpha _i}}}$ ， ${B_i} = {y_i} - {k_{\sigma i}}{x_i}$ ，(x_i, y_i)为i点坐标。

转动点(x_p, y_p)到 ${\sigma _i}$ 所在直线的距离为

${d}_{\sigma i}=\frac{\mid {k}_{\sigma i}{x}_{\text{p}}+(-1){y}_{\text{p}}+{B}_{i}\mid }{\sqrt{{k}_{\sigma i}{}^{2}+{(-1)}^{2}}} 。$

(7)

正应力 ${\sigma _i}$ 提供转动力矩，力矩和为 $\sum\limits_{i=1}^{n}{\sigma }_{i}\text{Δ}{l}_{i}{d}_{\sigma i}$ 。

滑动体的总的转动力矩和为

$\sum\limits_{i=1}^{n}{M}_{\text{s}}({x}_{\text{p}},{y}_{\text{p}}) ={\displaystyle \sum _{i=1}^{n}{\tau }_{i}\text{Δ}{l}_{i}{d}_{\tau i}} +\sum\limits_{i=1}^{n}{\sigma }_{i}\text{Δ}{l}_{i}{d}_{\sigma i} 。$

(8)

由于转动点的位置不同，滑动面的形状可能不规则，某一弧段的 ${\sigma _i}$ 可能提供下滑力矩也可能提供抗滑力矩，体现其计算力矩时可正可负（如果定义沿坡面倾向方向转动为正），但作为力矩和方法，不需要区分此时的应力是否是产生下滑力矩或抗滑力矩，只需将抗滑力矩设为负值，进行力矩和计算即可。在计算切应力 ${\tau _i}$ 时也如此，不需要区分，只需根据实际发生的计算应力进行计算即可。不需要区分任意弧段的 ${\sigma _i}$ 和 ${\tau _i}$ 产生的力矩的正负，更符合“和”方法的本意，而且在力学上，滑动面上的 ${\sigma _i}$ 和 ${\tau _i}$ ，等效于滑动体的重力和作用于滑动体上的全部外力，求和时结果是一致的。

（2）抗转动力矩和计算

岩土体的材料强度服从Mohr-Coulomb准则，滑动体沿滑动面发生剪切，其抗转动力矩则由滑动面上的切向抗剪强度和法向反力提供。对于滑动面上的任一微段 ${\Delta }{l_i}$ ，岩土体的抗剪强度 ${\tau _{{\text{f}}i}}$ 由莫尔库仑强度理论进行计算，方向与位移趋势方向相反，法向反力 ${\sigma '_i}$ = ${\sigma _i}$ ，方向相反。

${\tau _{{\text{f}}i}} = {c_i} + {\sigma _i}\tan {\varphi _i}。$

(9)

滑动面上切向抗剪强度和法向反力对转动点的力矩积分（或求和）构成了抗转动力矩和，即为式（3）中的分子项。

${\tau _{{\text{f}}i}}$ 所在直线与 ${\tau _i}$ 一致，故转动点(x_p, y_p)到 ${\tau _{{\text{f}}i}}$ 所在直线的距离也为 ${d_{\tau i}}$ ，提供抗转动力矩之和为 $\sum _{i=1}^{n}{\tau }_{\text{f}i}\text{Δ}{l}_{i}{d}_{\tau i}$ 。

转动点(x_p, y_p)到 ${\sigma '_i}$ 所在直线的距离也为 ${d_{\sigma i}}$ ，提供抗转动力矩之和为 $\sum\limits_{i=1}^{n}{{\sigma }^{\prime }}_{i}\text{Δ}{l}_{i}{d}_{\sigma i}$ 。

抗转动力矩和即为

$\sum\limits_{i=1}^n M_{\mathrm{r}}\left(x_{\mathrm{p}}, y_{\mathrm{p}}\right)=\sum\limits_{i=1}^n \tau_{\mathrm{f}i} \Delta l_i d_{\tau i}+\sum\limits_{i=1}^n \sigma_i^{\prime} \Delta l_i d_{\sigma i}。$

(10)

（3）安全系数的计算

由式（8），（10），并代入式（3），得到力矩和法安全系数的表达式（1）或（3）的计算式。

$F({x_{\text{p}}},{y_{\text{p}}}) = \frac{{\sum ({c_i} + {\sigma _i}\tan {\varphi _i}){\Delta }{l_i}{d_{\tau i}} + \sum {{\sigma '}_i}{\Delta }{l_i}{d_{\sigma i}}}}{{\sum {\tau _i}{\Delta }{l_i}{d_{\tau i}} + \sum {\sigma _i}{\Delta }{l_i}{d_{\sigma i}}}}\text{，}$

(11)

写成积分形式为

$F({x_{\text{p}}},{y_{\text{p}}}) = \frac{{\int {\left[ {({c_i} + {\sigma _i}\tan {\varphi _i}){d_{\tau i}} + {{\sigma '}_i}{d_{\sigma i}}} \right]{\text{d}}l} }}{{\int {({\tau _i}{d_{\tau i}} + {\sigma _i}{d_{\sigma i}}){\text{d}}l} }}。$

(12)

1.5 力矩和方法求解步骤

力矩和方法的基本计算步骤如下：

（1）根据边坡地质勘查资料，确定边坡的潜在滑动面和对应的滑动体。

（2）建立考虑所分析边坡和滑动面的力学（应力）分析模型（如有限元、离散元、条分法等可以获得滑动面应力状态的方法均可），根据边坡的载荷、边界条件和岩土体物理力学参数，得到潜在滑动面的计算应力分布状态。

（3）将潜在滑动面分为若干段，从分析模型中获得每段编号、位置坐标和几何参数，从计算所得应力分布状态获得每滑面段的应力。

（4）根据滑动面位置，拟定试算的转动点位置区域，依次以区域内节点为转动点，根据式（8），（10）计算滑动面上的转动力矩和与抗转动力矩和，再按式（11）或（12）计算安全系数。

（5）对所试算区域的所有转动节点，通过计算找到最小安全系数，如果最小安全系数对应节点位于区域内部，可以认为此时的最小安全系数为所分析滑动面的基于力矩和法分析的安全系数，及其对应的转动节点。若最小安全系数对应节点位于区域边界，则回到第（4）步骤，调整试算的转动点位置区域。

2. 基于条分法模型的力矩和方法求解

力矩和法具有物理含义明确、概念清晰的特点，它可以直接采用显式计算，在已经确定了滑动面位置前提下，只要得到边坡在荷载和自重下的计算应力场，就可以通过力矩和法计算边坡的稳定性安全系数。为了获得滑动面上的应力，力矩和法是需要结合其它分析方法的。在边坡稳定性分析方法中，条分法是目前应用最广泛的分析方法之一。为了获得滑动面的应力，本文将采用条分法计算模型，通过条分法的计算获得滑动面上的应力，从而进行力矩和方法的计算。基于条分法的滑动面应力求解，其一便于方法之间的对比特别是与广泛应用的条分法的结果之间的对比分析；其二是更易于实现，避免了有限元计算时对模型和程序的较多依赖。

条分法中，Morgenstern-Price方法（M-P法）是一种同时满足力和力矩平衡条件的通用条分法^[6]，它是边坡极限平衡条分法^[7-9]中最具普遍性和假设条件的一种分析方法。下面将介绍M-P法与力矩和法相结合，求解安全系数。

如图 3所示，将滑体划分成n个竖直条块，取任意条块i，其中条块i的底面长度是 ${l_i}$ ，底边与X轴的夹角是 ${\alpha _i}$ ，转动点坐标是(x_p, y_p)，条块i底面左下脚坐标是(x_i1, y_i1)，条块i底面右下脚坐标是(x_i, y_i)。W_i是条块i的重力，E_i1，E_i分别是条块i左右两侧的水平条间力，X_i1，X_i分别是条块i左右两侧的条间剪力，h_i-1为作用在条块i左侧水平条间力E_i-1离基底的距离，h_i为作用在条块i右侧水平条间力E_i离基底的距离，H是坡高，β是坡角，P_i是底部的法向力, T_i是底部的切向剪力, u_i是底部的孔隙水压， ${\varphi '_i}$ 是滑动表面的有效内摩擦角，有效黏聚力为 ${c'_i}$ ，绕点(x_p, y_p)转动的安全系数F(x_p, y_p)。

图 3 条分法安全系数求解示意图

Figure 3. Schematic diagram of solving factor of safety by slices method

下载: 全尺寸图片幻灯片

M-P法假定条间力函数关系为 ${X_i} = \lambda f(x){E_i}$ 。式中 $\lambda$ 为任意常数， $f(x)$ 可取半正弦函数、常数等。 $f(x) =$ ${\text{sin}}\left[ {(x - a)/(b - a)} \right]{\text{π }}$ （正弦曲线），其中a，b分别为滑动面两端的x坐标。如果 $f(x)$ =1，则条间力( ${X_i}$ 和 ${E_i}$ 的合力)方向平行，即为斯宾塞法对条间力的假设。当力与力矩的平衡被严格满足时，在一个合理的范围内，对条间力的不同假设对其安全系数的影响不大^[10]。

根据莫尔-库仑准则及安全系数F的定义，建立条块的水平和垂直方向力平衡、滑体整体水平向力平衡以及将作用在条块上的力对土条条块底中点取矩建立力矩平衡。为了便于计算力矩和安全系数，得到 ${P_i}$ 和 ${T_i}$ 的表达式，有如下公式：

${E}_{i}={E}_{i-1}+ {P}_{i}\mathrm{sin}{\alpha }_{i}-{T}_{i}\mathrm{cos}{\alpha }_{i} \text{，}$

(13)

${P}_{i}=\frac{{W}_{i}-({X}_{i}-{X}_{i-1})-\frac{1}{F}({{c}^{\prime }}_{i}{l}_{i}-{u}_{i}{l}_{i}\mathrm{tan}{{\varphi }^{\prime }}_{i})}{\mathrm{cos}{\alpha }_{i}+\frac{1}{F}(\mathrm{sin}{\alpha }_{i}\mathrm{tan}{{\varphi }^{\prime }}_{i})} \text{，}$

(14)

${T}_{i}=\frac{{{c}^{\prime }}_{i}{l}_{i}+({P}_{i}-{u}_{i}{l}_{i})\mathrm{tan}{{\varphi}^{\prime }}_{i}}{F} \text{，}$

(15)

$\sum {({P_i}\sin {\alpha _i} - {T_i}\cos {\alpha _i})} = 0 \text{，}$

(16)

${E_{i - 1}}\left( {{h_{i - 1}} + {y_{i - 1}} - \frac{{{y_{i - 1}} + {y_i}}}{2}} \right) - {E_i}({h_i} - \frac{{{y_{i - 1}} + {y_i}}}{2} + {y_i}) -\\ ~~~~~~~~~~ {X}_{i-1}\left(\frac{{x}_{i-1}+{x}_{i}}{2}-{x}_{i-1}\right)-{X}_{i}({x}_{i}-\frac{{x}_{i-1}+{x}_{i}}{2})\text{=}0\text{ }。$

(17)

岩土体的抗剪强度 ${\tau _{{\text{f}}i}}$ 提供抗滑力，其值由莫尔库仑强度准则进行计算，滑动面上条块i的 ${\tau _{{\text{f}}i}}$ 的斜率为 ${k_{\tau {\text{f}}i}}$ ，所在直线方程为

$y={k}_{\tau \text{f}i}x+{b}_{i} \text{，}$

(18)

式中： ${k_{\tau {\text{f}}i}} = \frac{{{y_i} - {y_{i - 1}}}}{{{x_i} - {x_{i - 1}}}}$ ； ${b_i} = {y_i} - \frac{{{y_i} - {y_{i - 1}}}}{{{x_i} - {x_{i - 1}}}}{x_i}$ 。

转动点(x_p, y_p)到 ${\tau _{{\text{f}}i}}$ 所在直线的距离为

${d}_{\tau \text{f}i}=\frac{\mid {k}_{\tau \text{f}i}{x}_{\text{p}}+(-1){y}_{\text{p}}+{b}_{i}\mid }{\sqrt{{k}_{\tau \text{f}i}{}^{2}+{(-1)}^{2}}} 。$

(19)

抗剪强度 ${\tau _{{\text{f}}i}}$ 提供抗转动力矩，之和为

$\sum {\tau _{{\text{f}}i}}{l_i}{d_{\tau {\text{f}}i}} = \sum [{c'_i}{l_i} + ({P_i} - {u_i}{l_i})\tan {\varphi '_i}]{d_{\tau {\text{f}}i}} 。$

(20)

式中： ${l_i}$ 为条块i的基底长度。

滑动面上条块i的基底反力 ${P}_{i}$ 的斜率为 ${k_{Pi}}$ ，所在直线方程为

$y={k}_{Pi}x+{B}_{i} 。$

(21)

式中： ${k_{Pi}} = - \frac{1}{{{k_{\tau {\text{f}}i}}}} = - \frac{{{x_i} - {x_{i - 1}}}}{{{y_i} - {y_{i - 1}}}}$ ； ${B_i} = \frac{{{y_i} + {y_{{\text{i}} - 1}}}}{2} - {k_{Pi}} \cdot$ $\frac{{{x_i} + {x_{i - 1}}}}{2}$ 。

转动点(x_p, y_p)到 ${P}_{i}$ 所在直线的距离为

${d}_{Pi}=\frac{\mid {k}_{Pi}{x}_{\text{p}}+(-1){y}_{\text{p}}+{B}_{i}\mid }{\sqrt{{k}_{Pi}{}^{2}+{(-1)}^{2}}} \text{，}$

(22)

${P}_{i}$ 提供的抗转动力矩和为 $\sum {{P_i}{d_{Pi}}}$ 。

抗转动力矩和即为

$\sum {M}_{\text{r}}\text{(}{x}_{\text{p}},{y}_{\text{p}}\text{)}=\sum {\tau }_{\text{f}i}{l}_{i}{d}_{\tau \text{f}i}+\sum {P}_{i}{d}_{Pi}。$

(23)

${W_i}$ 提供转动力矩，之和为 $\sum {{W_i}{d_{Wi}}}$ 。其中 ${d_{Wi}} = {x_{\text{p}}} -$ ${x_i}$ 。

力矩和安全系数表达式如下：

$F{\text{(}}{x_{\text{p}}},{y_{\text{p}}}{\text{)}} = \frac{{\sum {M_{\text{r}}}{\text{(}}{x_{\text{p}}},{y_{\text{p}}}{\text{)}}}}{{\sum {M_{\text{s}}}{\text{(}}{x_{\text{p}}},{y_{\text{p}}}{\text{)}}}} = \frac{{\sum {\tau _{{\text{f}}i}}{l_i}{d_{\tau fi}} + \sum {P_i}{d_{Pi}}}}{{\sum {W_i}{d_{Wi}}}}。$

(24)

3. 算例

3.1 圆弧滑动面边坡算例

选用ACADS考题1中EX1(c)，此边坡为非均质边坡，由3层土组成，其计算模型见图 4，材料参数见表 1。

图 4 ACADS考题EX1(c)计算模型

Figure 4. Model of problem of EX1(c) of ACADS

下载: 全尺寸图片幻灯片

表 1 ACADS考题EX1(c)材料参数

Table 1. Material parameters of EX1(c) of ACADS

土层	c/kPa	φ/(°)	γ/(kN·m^-3)
1	0	38	19.5
2	5.3	23	19.5
3	7.2	20	19.5

下载: 导出CSV

| 显示表格

本文暂时不对边坡稳定力矩和法的临界滑动面的搜索问题进行讨论，图 4中的滑动面为边坡稳定性分析软件slide求得的临界滑动面，本算例用此滑动面计算边坡稳定的力矩和法安全系数，其中考题EX1(c)的滑动面同时穿过3种土层。

表 2为条分法计算的条块数据。表 3为多种方法计算结果。计算均采用同一个滑动面，从表中可以看出，本文方法计算得到的安全系数最小，其值为1.043，且对于圆弧滑动面，当转动点为圆心时，本文方法计算结果为1.378，与裁判答案基本一致，误差小于2%。

表 2 条方法计算过程

Table 2. Calculation process by slices method

条块i	基底长l_i/ m	条重W_i/ (kN·m^-3)	抗滑力 ${\tau _{{\text{f}}i}}{l_i}$ /(kN·m^-1)	基底反力P_i/ (kN·m^-1)
1	5.23	110.11	75.25	96.31
2	3.97	232.52	143.43	183.58
3	3.30	261.32	110.98	220.29
4	3.35	284.60	119.98	263.44
5	3.19	243.01	117.80	260.58
6	3.13	167.48	95.96	201.65
7	3.17	59.87	50.78	76.88

下载: 导出CSV

| 显示表格

表 3 安全系数计算结果

Table 3. Calculation results of safety factor

计算方法	圆心坐标/m	转动点坐标/m	安全系数
费伦纽斯法	(35.77, 23.54)		1.263
毕肖普法	(35.77, 23.54)		1.372
简布法	(35.77, 23.54)		1.296
M-P法	(35.77, 23.54)		1.379
矢量和法			1.375
力矩和法		(33.5, 4.6)	1.043

下载: 导出CSV

| 显示表格

3.2 非圆弧滑面边坡算例

非圆弧滑面边坡算例选用ACADS提供的考题3中EX3(b)，EX3(b)是一个含软弱夹层的非均质边坡考核题目，边坡计算模型见图 5，岩土参数见表 4。

图 5 ACADS考题EX3(b)计算模型

Figure 5. Model of problem of EX3(b) of ACADS

下载: 全尺寸图片幻灯片

表 4 ACADS考题EX3(b)材料参数

Table 4. Material parameters of EX3(b) of ACADS

土层	c/kPa	φ/(°)	γ/(kN·m^-3)
1	28.5	20	18.84
2	0	10	18.84

下载: 导出CSV

| 显示表格

对于具有软弱夹层的边坡，其滑面必定有一段会沿软弱夹层滑动，此时若依然采用圆弧滑动面，得到的结果往往是不正确的，故这里采用指定的折线滑动面进行计算边坡稳定的力矩和法安全系数。图 5给出剖面几何特征和指定滑动面位置，将滑动体划分为7个条块，滑动面同时穿过2种土层，底部一段沿软弱夹层滑动，滑动面控制点坐标见表 5。

表 5 滑动面控制点坐标

Table 5. Coordinates of control points for sliding surface

点号	X/m	Y/m
1	10.69	20.00
2	20.50	7.00
3	40.00	6.50
4	42.15	7.75

下载: 导出CSV

| 显示表格

表 6为条分法计算的条块数据。表 7为多种方法计算结果。

表 6 条分法计算过程

Table 6. Calculation process by slices method

条块i	基底长l_i/m	条重W_i/ (kN·m^-3)	抗滑力 ${\tau _{{\text{f}}i}}{l_i}$ / (kN·m^-1)	基底反力P_i/ (kN·m^-1)
1	8.14	300.33	299.76	185.96
2	8.14	825.64	491.09	711.65
3	6.50	1158.27	207.78	1178.36
4	6.50	780.68	144.07	817.07
5	6.50	403.10	75.24	426.69
6	0.99	20.82	5.24	29.71
7	1.49	9.21	56.75	39.08

下载: 导出CSV

| 显示表格

表 7 安全系数计算结果

Table 7. Calculated results of factor of safety

计算方法		转动点坐标	安全系数
简布法			1.195
陆军工程师团法			1.364
M-P法			1.286
矢量和法			1.221
力矩和法	条间力函数为常数	(40, 6.5)	1.051
力矩和法	条间力函数为半正弦函数	(40, 6.5)	1.049

下载: 导出CSV

| 显示表格

从表 7可以看出，在非圆弧滑面的情况下，本文方法得到的安全系数较小，而不同的条间力假设对力矩和安全系数的影响并不明显，而对于含软弱夹层的非均质边坡，滑面偏离圆弧较大，本文方法的优越性是能够处理各种非圆弧滑面。

4. 结论

（1）提出了分析边坡稳定性的新方法即力矩和方法。该方法对一确定的滑坡体和滑动面，假定滑坡体绕某点转动，求出滑坡体在自重和外荷载作用下的转动力矩和，同时求出滑动面上基岩法向反力和剪切强度对滑动体的抗转动力矩和，定义抗转动力矩和与转动力矩和的比值为对应的安全系数。本文中，转动点由最小安全系数决定，并通过试算搜索得到。

（2）力矩和方法是对矢量和方法的一个补充。滑动体视为刚体时，它的位移应该包括刚体的平动和转动两部分（只有外力合力通过重心时没有外力引起的转动），矢量和方法认为滑动体是沿某一直线方向平动，不考虑其转动，而力矩和方法考虑了转动，两个方法互为补充。当然，若能同时考虑滑体的平动和转动，将会更合理，也是本文后续将延续的工作。

（3）本文在求解滑动面应力时，采用了更简便的条分法计算（M-P法）。并结合M-P法模型，推导了力矩和方法求解公式。

（4）算例1结果显示，对于所分析的滑动面，由于所选择的转动点不同，力矩和方法所得到的安全系数较其它几种方法更小。这主要是因为最终所得最小安全系数对应的转动点不是所分析滑动面圆弧的圆心点。滑动面圆弧的圆心点只是力矩和方法所分析转动点的一个，其对应的安全系数不一定是滑动体绕某点发生转动时的最小安全系数。

这一结果也提示一点，具有圆弧滑动面的滑坡体，其初始失稳转动时所对应的转动点（或圆心）不一定是圆弧滑动面对应的圆心。进一步可认为，任意形状的滑动面，对于滑动体，滑动面只是其约束边界条件，与滑动体初始失稳或破坏的位移趋势无直接关系。这与人们习惯的认知是不同的，值得后续进一步的研究。

算例1中，当转动点取为圆心时，本文方法求得的安全系数与裁判答案基本一致，说明力矩和方法计算安全系数是可行的。

（5）力矩和方法可以分析任意滑动面。滑动体所对应的滑动面，只是滑动体的边界条件，其形状不会直接影响滑动体作为刚体自身的位移趋势，而是通过边界提供反力影响滑动体的后续位移轨迹。算例2给出了本文方法分析非圆弧滑动面的情况，并可以搜索到比原结果较小的安全系数，说明假定滑动体绕某点转动的位移趋势是值得对边坡进行稳定性分析时考虑。

图 1 某地质剖面的钻孔分布图

Figure 1. Distribution of boreholes in a geological profile

下载: 全尺寸图片幻灯片

图 2 二维耦合马尔科夫链模型模拟地质剖面

Figure 2. Two-dimensional coupled Markov chain model for simulating geological cross-sections

下载: 全尺寸图片幻灯片

图 3 基于常规CMC获得的地层剖面图

Figure 3. Geological cross-sections obtained from conventional CMC model

下载: 全尺寸图片幻灯片

图 4 LCMC模型算法的实现流程图

Figure 4. Flow chart of algorithm of LCMC model

下载: 全尺寸图片幻灯片

图 5 CMC模型分解成片段

Figure 5. Decomposition of CMC model into segments

下载: 全尺寸图片幻灯片

图 6 CMC模型片段的观测场景

Figure 6. Observed scenarios of CMC model segments

下载: 全尺寸图片幻灯片

图 7 基于LCMC模型获得的地质剖面

Figure 7. Geological cross-sections obtained from LCMC model

下载: 全尺寸图片幻灯片

图 8 基于LCMC的地质不确定性结果（信息熵）

Figure 8. Results of geological uncertainty based on LCMC (I_e)

下载: 全尺寸图片幻灯片

图 9 基于柱平均信息熵的钻孔优化和地质剖面更新

Figure 9. Optimization of boreholes and updating of geological cross-sections based on column-wise average information entropy

下载: 全尺寸图片幻灯片

图 10 钻孔数量对平均信息熵的影响以及钻孔优化效率比较

Figure 10. Influences of number of boreholes on average information entropy and comparisons of borehole optimization efficiency

下载: 全尺寸图片幻灯片

表 1 基于钻孔数据获得的转移次数矩阵

Table 1 Transition count matrices obtained from borehole data

(a)垂直转移次数矩阵（VTCM）	状态	砂土	黏土	粉土	粉质黏土
	砂土	245	0	0	12
	黏土	7	455	4	7
	粉土	8	9	399	0
	粉质黏土	0	4	18	392
(b)向左转移次数矩阵（LTCM）	砂土	245	34	64	69
	黏土	32	577	84	33
	粉土	60	58	399	107
	粉质黏土	102	16	87	392
(c)向右转移次数矩阵（RTCM）	砂土	245	32	60	102
	黏土	34	577	58	16
	粉土	64	84	399	87
	粉质黏土	69	33	107	392

下载: 导出CSV

参考文献(20)

[1]	SCHECK D E, CHOU D R. Optimum locations for exploratory drill holes[J]. International Journal of Mining Engineering, 1983, 1(4): 343-355. doi: 10.1007/BF00881550
[2]	HOSSEIN MORSHEDY A, MEMARIAN H. A novel algorithm for designing the layout of additional boreholes[J]. Ore Geology Reviews, 2015, 67: 34-42. doi: 10.1016/j.oregeorev.2014.11.012
[3]	SHURYGIN A M. The probability of finding deposits and some optimal search grids[J]. Journal of the International Association for Mathematical Geology, 1976, 8(3): 323-330. doi: 10.1007/BF01029276
[4]	宣腾. 基于克里金法的地质勘探位置分析[D]. 哈尔滨: 哈尔滨工业大学, 2016. XUAN Teng. Location Analysis of The Site Investigation Using Kriging[D]. Harbin: Harbin Institute of Technology, 2016. (in Chinese)
[5]	LIU Y H, JOURNEL A. Improving sequential simulation with a structured path guided by information content[J]. Mathematical Geology, 2004, 36(8): 945-964. doi: 10.1023/B:MATG.0000048800.72104.de
[6]	DREW L J. Pattern drilling exploration: optimum pattern types and hole spacings when searching for elliptical shaped targets[J]. Journal of the International Association for Mathematical Geology, 1979, 11(2): 223-254. doi: 10.1007/BF01028966
[7]	SHI C, WANG Y. Nonparametric and data-driven interpolation of subsurface soil stratigraphy from limited data using multiple point statistics[J]. Canadian Geotechnical Journal, 2021, 58(2): 261-280. doi: 10.1139/cgj-2019-0843
[8]	SHI C, WANG Y. Data-driven sequential development of geological cross-sections along tunnel trajectory[J]. Acta Geotechnica, 2023, 18(4): 1739-1754. doi: 10.1007/s11440-022-01707-1
[9]	邓志平, 李典庆, 祁小辉, 等. 基于广义耦合马尔可夫链的地层变异性模拟方法[J]. 岩土工程学报, 2018, 40(11): 2041-2050. doi: 10.11779/CJGE201811010 DENG Zhiping, LI Dianqing, QI Xiaohui, et al. Simulation of geological uncertainty using modified generalized coupled Markov chain[J]. Chinese Journal of Geotechnical Engineering, 2018, 40(11): 2041-2050. (in Chinese) doi: 10.11779/CJGE201811010
[10]	ZHANG J Z, HUANG H W, ZHANG D M, et al. Quantitative evaluation of geological uncertainty and its influence on tunnel structural performance using improved coupled Markov chain[J]. Acta Geotechnica, 2021, 16(11): 3709-3724. doi: 10.1007/s11440-021-01287-6
[11]	SHI C, WANG Y. Machine learning of three-dimensional subsurface geological model for a reclamation site in Hong Kong[J]. Bulletin of Engineering Geology and the Environment, 2022, 81(12): 504. doi: 10.1007/s10064-022-03009-y
[12]	李希元, 崔健, 胡望水, 等. 基于多源地球物理数据的机器学习方法在地质体分类中的应用: 以黑龙江多宝山矿集区为例[J]. 地球物理学报, 2022, 65(9): 3634-3649. LI Xiyuan, CUI Jian, HU Wangshui, et al. Application of machine learning method based on multi-source geophysical data to geological body classification—a case study of Duobaoshan ore concentration area (Heilongjiang, China)[J]. Chinese Journal of Geophysics, 2022, 65(9): 3634-3649. (in Chinese)
[13]	KRUMBEIN W C. Statistical models in sedimentology1[J]. Sedimentology, 1968, 10(1): 7–23. doi: 10.1111/j.1365-3091.1968.tb01908.x
[14]	ELFEKI A M M, DEKKING F M. Modelling subsurface heterogeneity by coupled Markov chains: directional dependency, walther's law and entropy[J]. Geotechnical & Geological Engineering, 2005, 23(6): 721-756.
[15]	CHESSA A G. A Markov chain model for subsurface characterization: theory and applications[J]. Mathematical Geology, 2006, 38(4): 503-505. doi: 10.1007/s11004-006-9037-9
[16]	LI J H, CAI Y M, LI X Y, et al. Simulating realistic geological stratigraphy using direction-dependent coupled Markov chain model[J]. Computers and Geotechnics, 2019, 115: 103147. doi: 10.1016/j.compgeo.2019.103147
[17]	ZENG Z Q, CAI Y C, WU J B. Local coupled Markov chain model for simulating varied stratigraphy[J]. Computers and Geotechnics, 2024, 167: 106064. doi: 10.1016/j.compgeo.2023.106064
[18]	程利力, 陈健, 陈睿, 等. 基于二维马尔可夫链的武汉长江公铁隧道地层识别[J]. 土木工程与管理学报, 2021, 38(1): 169-174, 182. CHENG Lili, CHEN Jian, CHEN Rui, et al. Stratigraphic identification of Wuhan Yangtze River railway tunnel based on two-dimensional Markov chain[J]. Journal of Civil Engineering and Management, 2021, 38(1): 169-174, 182. (in Chinese)
[19]	ZHANG J Z, LIU Z Q, ZHANG D M, et al. Improved coupled Markov chain method for simulating geological uncertainty[J]. Engineering Geology, 2022, 298: 106539. doi: 10.1016/j.enggeo.2022.106539
[20]	丁文峰, 韩昊宇, 王协康, 等. 基于信息熵理论的黄土高原小流域雨量站网优化[J]. 工程科学与技术, 2022, 54(2): 113-121. DING Wenfeng, HAN Haoyu, WANG Xiekang, et al. Rainfall networks optimization in small catchment of Loess Plateau based on information entropy[J]. Advanced Engineering Sciences, 2022, 54(2): 113-121. (in Chinese)