极限分析法求解含软弱夹层边坡稳定性

徐鹏; 尚念林; 鲍晶晶; 李婷

doi:10.3969/j.issn.0258-2724.20200156

极限分析法求解含软弱夹层边坡稳定性

doi: 10.3969/j.issn.0258-2724.20200156

徐鹏^1,,
尚念林²,
鲍晶晶³,
李婷^{1, 4, ,}

1.
石家庄铁道大学省部共建交通工程结构力学行为与系统安全国家重点实验室，河北石家庄 050043
2.
山东宝泰置业有限公司，山东济南 271100
3.
西南交通大学土木工程学院，四川成都 610031
4.
石家庄铁道大学土木工程学院，河北石家庄 050043

基金项目: 国家自然科学基金(52108331)；河北省自然科学基金(E2021210010)；河北省高等学校科学技术研究项目(ZD2021096，QN2021127)；河北省引进留学人员资助项目（C20210305）

详细信息

作者简介:
徐鹏（1988—），男，讲师，研究方向为道路与铁道工程，E-mail：bk20090201@my.swjtu.edu.cn

通讯作者:
李婷（1992—），女，讲师，研究方向为地基变形控制与路基工程等，E-mail：liting9225@126.com

中图分类号: TU41
计量
- 文章访问数: 386
- HTML全文浏览量: 181
- PDF下载量: 40
- 被引次数: 25
出版历程
- 收稿日期: 2020-04-02
- 修回日期: 2020-09-04
- 刊出日期: 2020-09-16

Stability Analysis of Slopes with Weak Layers Using Limit Analysis Method

XU Peng^1
,,
SHANG Nianlin²,
BAO Jingjing³,
LI Ting^{1, 4
, ,}

1.
State Key Laboratory of Mechanical Behavior and System Safety of Traffic Engineering Structures, Shijiazhuang Tiedao University, Shijiazhuang 050043, China
2.
Shandong Baotai Real Estate Co., Ltd., Jinan 271100, China
3.
School of Civil engineering, Southwest Jiaotong University, Chengdu 610031, China
4.
School of Civil engineering, Shijiazhuang Tiedao University, Shijiazhuang 050043, China

摘要

摘要:
软弱夹层对边坡的稳定性影响显著，目前设计中通常采用极限平衡法计算边坡的稳定性，其在求解中需要建立多个平衡方程. 为了分析含软弱夹层边坡的稳定性，首先，采用极限分析法建立了计算模型；其次，通过极限平衡法验证了求解的准确性；最后，分析了荷载、夹层形状、夹层强度等对稳定性的影响. 研究结果表明：边坡安全系数随着外荷载强度的增大而减小，其中，当加速度放大系数由1.0增大为1.6时，安全系数由1.20降为0.89；当外荷载频率越大时，边坡越易提前产生破坏；软弱夹层形状对边坡安全系数影响显著，特别是当其靠近坡顶与坡面时；安全系数随着软弱夹层摩擦角与黏聚力的减小而近似线性降低，其中，当黏聚力由9 kPa降为5 kPa时，安全系数降低约30%.
- 软弱夹层 /
- 边坡 /
- 极限平衡法 /
- 极限分析法 /
- 安全系数
Abstract:
Weak layers have a significant effect on the stability of slopes. The stability of slopes is usually calculated by the limit equilibrium method in current designs, in which multiple equilibrium equations need to be established and solved. Compared with the limit equilibrium method, the limit analysis method is more rigorous, and requires only one energy equation to be solved. In order to analyze the stability of slopes with weak interlayers, a stability calculation model was established by the limit analysis method, and then the accuracy of the solution was verified by the limit equilibrium method. Finally, the effects of load, weak layer shape, and weak layer strength on stability were analyzed. The results show that the slope safety factor decreases with an increase in the load intensity. When the acceleration amplification factor increases from 1.0 to 1.6 the safety factor decreases from 1.20 to 0.89. With a higher frequency of the external load, the slope is easier to be damaged in advance. Besides, the shape of the weak layer has a significant effect on the slope safety factor, especially when it is close to the top and surface of the slope. The safety factor decreases approximately linearly with the decrease of the friction angle and cohesion of weak layers. When the cohesion strength decreases from 9 kPa to 5 kPa, the safety factor decreases by about 30%.
- weak layer /
- slope /
- limit equilibrium method /
- limit analysis method /
- safety factor

HTML全文

直线拟合不仅是曲线拟合研究的热点，并且在工程实践中被应用广泛^[1-2]. 直线拟合是由n个测量点 $({x_i},{y_i}),\;i = 1,2, \cdots ,n$ ，基于最小二乘准则找到一条最佳拟合直线. x、y分别为自变量、因变量， $a$ 、 $b$ 分别为直线的斜率和截距，拟合直线的方程为y = ax + b.

普通最小二乘认为变量x或y无误差去估计参数 $a$ 和 $b$ ，该方法简单明了，但选择自变量和因变量不同，拟合的直线不同. 正交最小二乘认为x和y具有相同的精度，几何意义是测量点到拟合直线的垂直距离平方和最小，通常被认为是最佳拟合^[3]，在铁路整正工程中被广泛采用^[4-5].

为顾及自变量x存在误差，通常采用EIV （errors-in-variables）模型进行整体最小二乘估计直线参数^[6-7]. 整体最小二乘通过建立变量的随机模型可以实现普通及正交最小二乘拟合直线，同时还可以实现加权整体最小二乘拟合直线^[8]. 但此时还不能考虑自变量x和因变量y间的误差相关性. 当已知观测误差的随机特性时，参数的最优估计应符合这种随机特性^[9]. Amiri-Simkooei等^[10]进一步提出了考虑观测值误差充分相关性的加权整体最小二乘法拟合直线.

直线是轨道线形的组成部分，为保证线路平顺性，需要估计直线参数，进行既有轨道线形整正. 采用全站仪采集的数据点，其纵、横坐标不仅精度不同，并且误差具有相关性. 因此，考虑观测值误差相关性的加权整体最小二乘才能处理这种情况，实现直线的最佳拟合. 整体最小二乘法是基于拟合直线得出的，要处理系数矩阵含有误差及与因变量y误差相关等情况，引入了EIV模型及Kronecker积等复杂矩阵运算^[11]. 这在一定程度上造成了理解及工程应用上的困难^[12]. 同时，线路中直线的拟合还受到相邻线元的约束.

因此，本文基于极大似然估计及拉格朗日条件极值原理建立直线拟合模型，引入表征测点位置的附加参数，推导出了顾及约束和观测值误差相关性直线拟合的通用方法. 实验验证了该方法能在任何误差分布情况下顾及约束估计直线参数及其精度.

1. 纵、横坐标的误差相关性

在工程实际中，经常使用全站仪获取点位坐标，测量原理见. 图中：XOY坐标系中，A $({x_A},{y_A})$ 和B $({x_B},{y_B})$ 为控制点，P(x,y)为全站仪测量点， ${\alpha _{AB}}$ 和 ${\alpha _{AP}}$ 分别为AB边和AP边的方位角. 测得独立观测量水平距离s和水平角β，则P点坐标为

图 1 全站仪测量原理

Figure 1. Measuring principle

下载: 全尺寸图片幻灯片

$\left\{ {\begin{array}{*{20}{l}} {x = {x_A} + s\cos ({\alpha _{AB}} + \beta ) = {x_A} + s\cos \;{\alpha _{AP}}{\text{，}}}\\ {y = {y_A} + s\sin ({\alpha _{AB}} + \beta ) = {y_A} + s\sin \;{\alpha _{AP}}.} \end{array}} \right.$

(1)

由式（1）和误差传播定律可计算点P纵、横坐标的方差-协方差阵为

${\boldsymbol{c}} = \left[ {\begin{array}{*{20}{c}} {\sigma _x^2}&{{\sigma _{xy}}} \\ {{\sigma _{yx}}}&{\sigma _y^2} \end{array}} \right] = {\boldsymbol{K}}\left[ {\begin{array}{*{20}{c}} {\sigma _s^2}&0 \\ 0&{\sigma _\beta ^2} \end{array}} \right]{{\boldsymbol{K}}^{\rm T}},$

(2)

式中： $\sigma _x^{\rm{2}}$ 和 $\sigma _y^{\rm{2}}$ 分别为x、y的方差； ${\sigma _{xy}}$ 、 ${\sigma _{yx}}$ 分别为x和y、y和x间的协方差； ${\sigma _s}$ 、 ${\sigma _\beta }$ （单位：″ ）分别为全站仪测距和测角的先验精度；矩阵 ${\boldsymbol{K}} = \left[ {\begin{array}{*{20}{c}} {\cos \;{\alpha _{AP}}}&{ - \dfrac{s}{\rho }\sin \;{\alpha _{AP}}} \\ {\sin \;{\alpha _{AP}}}&{\dfrac{s}{\rho }\cos \;{\alpha _{AP}}} \end{array}} \right]$ .

将式（2）展开得

$\sigma ^2_x = \sigma _s^2{\cos ^2}{\alpha _{AP}} + \frac{{\sigma _\beta ^2{s^2}}}{{{\rho ^2}}}{\sin ^2}{\alpha _{AP}}{\text{，}}$

(3)

$\sigma _y^2 = \sigma _s^2{\sin ^2}{\alpha _{AP}} + \frac{{\sigma _\beta ^2{s^2}}}{{{\rho ^2}}}{\cos ^2}{\alpha _{AP}}{\text{，}}$

(4)

${\sigma _{xy}} = {\sigma _{yx}} = 0.5\left( {\sigma _s^2 - \frac{{\sigma _\beta ^2{s^2}}}{{{\rho ^2}}}} \right)\sin (2{\alpha _{AP}}),$

(5)

式（3）～（5）中： $\rho {\rm{ \; = 206\;265''}}$ .

点P坐标观测值l = [x y] 的期望在拟合直线上，对应投影位置为 $\boldsymbol{\mu }(\boldsymbol{\hat \varTheta }) = [\hat x\;\;\hat a\hat x + \hat b]$ ，其中：参数估值 $\boldsymbol{\hat \varTheta }= {[\hat a\; \hat b\; \hat x]^{\rm{T}}}$ ， $\hat a$ 和 $\hat b$ 分别为a和b的估值， $\hat x$ 为点P位置x的估值. 坐标改正向量 $\boldsymbol{r}(\boldsymbol{\varTheta }) = {[{r_x}\;{r_y}]^{\rm{T}}} = {\left( { \boldsymbol{l} - \boldsymbol{\mu }(\boldsymbol{\varTheta })} \right)^{\rm{T}}}$ ，其中： $\boldsymbol{\varTheta }$ 为参数向量，r_x和r_y分别为x和y的坐标改正数. 点P坐标观测值正态分布密度函数为

$f(\boldsymbol{l};\boldsymbol{\varTheta }) = {\left| {2{\text{π}} \boldsymbol{c}} \right|^{ - 1/2}}\exp \left( { - 0.5{\boldsymbol{r}^{\rm{T}}}(\boldsymbol{\varTheta }){\boldsymbol{c}^{ - 1}}\boldsymbol{r}(\boldsymbol{\varTheta })} \right).$

(6)

最小二乘原理为 ${\rm{min }}\{{p_x}r_x^2 + {p_y}r_y^2\}$ ，其中 ${p}_{x}$ 和 ${p}_{y}$ 分别为 ${r_x}$ 和 ${r_y}$ 的权重. 当坐标分量同精度时，最小二乘原理为 ${\rm{min }}\{r_x^2 + r_y^2\}$ . 此时，中测量点P的坐标改正向量 ${[{r_{{x_1}}}\;{r_{{y_1}}}]^{\rm{T}}}$ 垂直拟合直线，投影点为P₁，即正交最小二乘. 当纵坐标x无误差时，其权 ${p}_{x}$ 为无穷大，则 ${r_x}$ 必然为0. 此时，测量点P的坐标改正向量 ${[0\;{r_{{y_2}}}]^{\rm{T}}}$ 垂直X轴，在拟合直线上的投影点为P₂，即普通最小二乘.

图 2 拟合原理

Figure 2. Fitting principle

下载: 全尺寸图片幻灯片

对于全站仪测得的采样点，坐标分量不仅精度不同（ $\sigma _x^2 \ne \sigma _y^2$ ），并且误差具有相关性（ ${\sigma _{xy}} \ne 0$ ）. 此时，普通最小二乘或正交最小二乘就不再适合.

2. 约束下顾及误差相关性的直线拟合

2.1 直线拟合模型

直线上n个观测点坐标组成观测向量l，其联合密度函数为

$L(\boldsymbol{l};\boldsymbol{\varTheta }) = \prod\limits_{i = 1}^n {f({\boldsymbol{l}_i};\boldsymbol{\varTheta })} {\text{，}}$

(7)

式中：l_i = [x_i y_i].

$\boldsymbol{\varTheta }$ 的极大似然估计量为

$\boldsymbol{\hat \varTheta }{\rm{ = }}\arg \max \;L(\boldsymbol{l};\boldsymbol{\varTheta }).$

(8)

由式（6）～（8）可推出 $\boldsymbol{\varTheta }$ 的极大似然估计等同于

$\min \left\{ {\sum\limits_{i = 1}^n {\boldsymbol{r}_i^{\rm{T}}(\boldsymbol{\varTheta }){\boldsymbol{c}_i}^{ - 1}{\boldsymbol{r}_i}(\boldsymbol{\varTheta })} } \right\}{\text{，}}$

(9)

式中： ${\boldsymbol{r}_i}(\boldsymbol{\varTheta }) = {[{r_{{x_i}}}\;{r_{{y_i}}}]^{\rm{T}}}$ ；c_i为测点i的方差-协方差矩阵.

令 ${\boldsymbol{r}^{\rm{T}}}\left( \boldsymbol{\varTheta }\right) = \left[ {\boldsymbol{r}_1^{\rm{T}}(\boldsymbol{\varTheta })\;\boldsymbol{r}_2^{\rm{T}}(\boldsymbol{\varTheta })\; \cdots\;\boldsymbol{r}_n^{\rm{T}}(\boldsymbol{\varTheta }]} \right)$ 及n个观测点的方差-协方差矩阵 $\boldsymbol{C} = {\rm{diag}} ({\boldsymbol{c}_1},{\boldsymbol{c}_2}, \cdots ,{\boldsymbol{c}_n})$ . 则 $\boldsymbol{C}$ 、协因数矩阵Q和权阵P即为随机模型，满足

$\boldsymbol{C}{\rm{ = }}\sigma _0^2\boldsymbol{Q} = \sigma _0^2{\boldsymbol{P}^{ - 1}}{\text{，}}$

(10)

式中： $\sigma _0^2$ 为先验单位权方差.

式（9）可写为

${\rm{min}}\left\{ {{\boldsymbol{r}^{\rm{T}}}\left( \boldsymbol{\varTheta }\right)\boldsymbol{Pr}\left( \boldsymbol{\varTheta }\right)} \right\}.$

(11)

即观测值正态分布时，极大似然估计和最小二乘估计等价. 易知，当P为单位矩阵I时，式（11）退化为 ${\rm{min}}\left\{ {{\boldsymbol{r}^{\rm{T}}}\left( \boldsymbol{\varTheta } \right)\boldsymbol{r}\left( \boldsymbol{\varTheta } \right)} \right\} = {\rm{min}}\left\{ {\displaystyle\sum\limits_{i = 1}^n {(r_{{x_i}}^2 + r_{{y_i}}^2} )} \right\}$ ，即正交最小二乘；当x权无穷大时，式（11）退化为 ${\rm{min}}\displaystyle\sum\limits_{i = 1}^n {r_{{y_i}}^2}$ ，即普通最小二乘.

第i个观测点的直线参数方程 ${\boldsymbol{\mu }_i}(\boldsymbol{\varTheta })$ 可表示为

$\left\{ {\begin{array}{l} {{{\hat x}_i} = {{\tilde x}_i} + {\rm{\delta }}{x_i},} \\ {{\hat y}_i} = (\tilde a + {\rm{\delta }}a) ({{\tilde x}_i} + {\rm{\delta }}{x_i}) + (\tilde b + {\rm{\delta }}b)\;{\rm{ =}}\\ \quad{{\tilde y}_i} + {{\tilde x}_i}{\rm{\delta }}a + \tilde a{\rm{\delta }}{x_i} + {\rm{\delta }}b + {\rm{\delta }}a{\rm{\delta }}b, \\ \end{array}} \right.$

(12)

式中： $(\tilde x,\tilde y)$ 和 $(\hat x,\hat y)$ 分别为点的近似坐标和估值坐标； $\tilde a$ 和 $\tilde b$ 分别为斜率和截距的近似值.

舍去高次微小量 ${\rm{\delta }}a{\rm{\delta }}b$ 后，式（12）线性化为

$\left\{ {\begin{array}{l} {{{\hat x}_i} = {{\tilde x}_i} + {\rm{\delta }}{x_i}{\text{，}}}\\ {{{\hat y}_i} = {{\tilde y}_i} + {{\tilde x}_i}{\rm{\delta }}a + \tilde a{\rm{\delta }}{x_i} + {\rm{\delta }}b.} \end{array}} \right.$

(13)

因此，线性化的直线方程可简化成矩阵形式如式（14）.

$\boldsymbol{\mu }({\boldsymbol{\hat \varTheta }}) = \boldsymbol{\mu }({\boldsymbol{\tilde \varTheta }}) + {\boldsymbol{A}}{\rm{\delta }}{\boldsymbol{\varTheta }}{\text{，}}$

(14)

式中： ${\boldsymbol{\tilde \varTheta }}$ 为 ${\boldsymbol{\varTheta}}$ 近似值； ${{\boldsymbol{A}}_{2n\times(n + 2)}} = \left[ {\begin{array}{*{20}{c}} {\bf{0}}&{\boldsymbol{I}} \\ {\boldsymbol{M}}&{\boldsymbol{N}} \end{array}} \right]$ ，其中 ${{\boldsymbol{M}}_{n \times 2}} = {\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{{\tilde x}_1}}&{{{\tilde x}_{\rm{2}}}}& \cdots &{{{\tilde x}_n}} \end{array}} \\ {\begin{array}{*{20}{c}} 1&1& \cdots &1 \end{array}} \end{array}} \right]^{\rm{T}}}$ ， ${{\boldsymbol{N}}_{n \times n}} = {\rm{diag}}(\tilde a)$ .

由式（11）、（14），目标函数可写为

$\min \left\{ {{{\left[ {\boldsymbol{A}{\rm{\delta }}\boldsymbol{\varTheta } - \boldsymbol{r}\left( {\boldsymbol{\tilde \varTheta }} \right)} \right]}^{\rm{T}}}\boldsymbol{P}\left[ {\boldsymbol{A}{\rm{\delta }}\boldsymbol{\varTheta } - \boldsymbol{r}\left( {\boldsymbol{\tilde \varTheta }} \right)} \right]} \right\}.$

(15)

线路中直线的拟合还要满足与相邻线元相切等约束条件，如已知相邻线元为圆曲线，其圆心位置（x₀，y₀）和半径R已知，则直线拟合受到的约束为

$\frac{{{y_0} - a{x_0} - b}}{{\sqrt {1 + {a^2}} }} - R = 0.$

(16)

线性化后为

$\begin{split} &{\frac{{ - {x_0}(1 + {{\tilde a}^2}) - \left( {{y_0} - \tilde a{x_0} - \tilde b} \right)\tilde a}}{{1{\rm{ + }}{{\tilde a}^2}}}{\rm{\delta }}a - }\\ &\quad{{\rm{\delta }}b + {y_0} - \tilde a{x_0} - \tilde b - R(1 + {{\tilde a}^2}) = 0}. \end{split}$

(17)

约束条件线性化后，矩阵表示为

$\boldsymbol{B}{\rm{\delta }}\boldsymbol{\varTheta } + \boldsymbol{w} = \bf{0}{\text{，}}$

(18)

式中：B为系数矩阵，其秩为约束条件的个数；w为闭合差向量.

构造拉格朗日极值函数为

$\begin{split} &\min \left\{ {{\left[ {\boldsymbol{A}{\rm{\delta }}\boldsymbol{\varTheta } - \boldsymbol{r}\left( {\boldsymbol{\tilde \varTheta }} \right)} \right]}^{\rm{T}}}\boldsymbol{P}\left[ {\boldsymbol{A}{\rm{\delta }}\boldsymbol{\varTheta } - \boldsymbol{r}\left( {\boldsymbol{\tilde \varTheta }} \right)} \right] +\right.\\ &\quad\bigg.{\boldsymbol{\gamma }^{\rm{T}}}\left( {\boldsymbol{B}{\rm{\delta }}\boldsymbol{\varTheta } + \boldsymbol{w}} \right) \bigg\}, \end{split}$

(19)

式中： $\boldsymbol{\gamma }$ 为拉格朗日乘向量.

为得到极小值，对式（19）求一阶导数并令其为0，可得

${\boldsymbol{A}^{\rm{T}}}\boldsymbol{PA}{\rm{\delta }}\boldsymbol{\varTheta } - {\boldsymbol{A}^{\rm{T}}}\boldsymbol{Pr}\left( {\boldsymbol{\tilde \varTheta }} \right) + {\boldsymbol{B}^{\rm{T}}}\boldsymbol{\gamma } = \bf{0}\bf{.}$

(20)

令 $\boldsymbol{N} = {\boldsymbol{A}^{\rm{T}}}\boldsymbol{PA}$ 及 $\boldsymbol{M} = {\boldsymbol{A}^{\rm{T}}}\boldsymbol{Pr}\left( {\boldsymbol{\tilde \varTheta }} \right)$ ，式（20）写为

$\boldsymbol{N}{\rm{\delta }}\boldsymbol{\varTheta } - \boldsymbol{M} + {\boldsymbol{B}^{\rm{T}}}\boldsymbol{\gamma } = \bf{0}\bf{.}$

(21)

式（21）乘以 $\boldsymbol{B}{\boldsymbol{N}^{ - 1}}$ ，再减去式（18）可得

$\boldsymbol{B}{\boldsymbol{N}^{ - 1}}{\boldsymbol{B}^{\rm{T}}}\boldsymbol{\gamma } - \boldsymbol{B}{\boldsymbol{N}^{ - 1}}\boldsymbol{M} - \boldsymbol{w} = \bf{0}\bf{.}$

(22)

令 ${\boldsymbol{N}_{\rm{c}}} = \boldsymbol{B}{\boldsymbol{N}^{ - 1}}{\boldsymbol{B}^{\rm{T}}}$ ，由式（22）可求得

$\boldsymbol{\gamma } = \boldsymbol{N}_{\rm{c}}^{ - 1}\left( {\boldsymbol{B}{\boldsymbol{N}^{ - 1}}\boldsymbol{M}{\rm{ + }}\boldsymbol{w}} \right).$

(23)

将式（23）代入式（21）得

$\begin{split} &{\rm{\delta }}\boldsymbol{\varTheta } = {\boldsymbol{N}^{ - 1}}\left( {\boldsymbol{I} - {\boldsymbol{B}^{\rm{T}}}\boldsymbol{N}_{\rm{c}}^{ - 1}\boldsymbol{B}{\boldsymbol{N}^{ - 1}}} \right)\boldsymbol{M} \; - \\ &\quad{\boldsymbol{N}^{ - 1}}{\boldsymbol{B}^{\rm{T}}}\boldsymbol{N}_{\rm{c}}^{ - 1}\boldsymbol{w}. \end{split}$

(24)

当无约束时，B为空矩阵，则式（24）变为

${\rm{\delta }}\boldsymbol{\varTheta } = {\boldsymbol{N}^{ - 1}}\boldsymbol{M}.$

(25)

约束条件下参数的估值为

${\boldsymbol{\hat \varTheta }} = {\boldsymbol{\tilde \varTheta }} + {\rm{\delta }}{\boldsymbol{\varTheta }}.$

(26)

2.2 参数精度

求得参数的最佳估值，需进一步评定其精度. 依据协因数传播定律，由式（24）得到参数估值 ${\boldsymbol{\hat \varTheta }}$ 的协因数阵为

${\boldsymbol{Q}_{\boldsymbol{\hat \varTheta }}} = {\boldsymbol{N}^{ - 1}}(\boldsymbol{I} - {\boldsymbol{B}^{\rm{T}}}\boldsymbol{N}_{\rm{c}}^{ - 1}\boldsymbol{B}{\boldsymbol{N}^{ - 1}}).$

(27)

同理，当无约束时，式（27）变为

${\boldsymbol{Q}_{\boldsymbol{\hat \varTheta }}} = {\boldsymbol{N}^{ - 1}}{\rm{ = }}{({\boldsymbol{A}^{\rm{T}}}\boldsymbol{PA})^{ - 1}}.$

(28)

因此，附加约束后，参数的协因数将减小.

由式（10）可得 ${\boldsymbol{\hat \varTheta }}$ 的方差-协方差如式（29）所示.

${\boldsymbol{C}_{\boldsymbol{\hat \varTheta }}} = \hat \sigma _0^2{\boldsymbol{Q}_{\boldsymbol{\hat \varTheta }}}{\text{，}}$

(29)

$\hat \sigma _0^2 = \frac{{\boldsymbol{r}{{(\boldsymbol{\hat \varTheta })}^{\rm{T}}}\boldsymbol{Pr}(\boldsymbol{\hat \varTheta })}}{{2n - t + m}},$

(30)

式中： ${\hat \sigma _0}$ 为后验单位权中误差，反映了观测值与模型之间的拟合程度；t为 ${\boldsymbol{\hat \varTheta }}$ 中元素个数；m为约束个数.

2.3 迭代寻优过程

由于拟合模型舍去了高次项，参数初值为近似值. 因此，需要使用高斯-牛顿迭代寻优算法进行迭代计算，直至 ${\rm{\delta }}{\boldsymbol{\varTheta }}$ 向量趋于0. 迭代终止条件设为 $\left\| {{\rm{\delta }}{\boldsymbol{\varTheta }}} \right\| \leqslant \varepsilon$ ， $\varepsilon$ 为阈值，取为10⁻⁶. 算法流程见图3.

图 3 高斯-牛顿迭代算法

Figure 3. Gauss-Newton iteration algorithm

下载: 全尺寸图片幻灯片

3. 试验算例

选取国内某既有专用线（设计速度30 km/h）直线段的复测数据，每20 m测一个点，共20个测点，坐标数据见表1. 使用的全站仪的测角精度为2″，测距精度为（2 + 2 × 10⁻⁶D） mm，D为距离，mm. 控制点A的坐标为（1000.012，1500.023），后视方位角为45°00′05″，由式（3）～（5）可计算出各测点纵、横坐标的方差及协方差，进而组成方差-协方差矩阵C. 取先验单位权中误差 ${\sigma _0}{\rm{ = }}1$ ，由式（10）可计算出实际测点的随机模型P₁，其各项元素列于表1. 表中：p_x、p_y和p_xy分别为式（2）中c⁻¹对应对角线元素和非对角线元素. 结果表明观测点的坐标分量不仅精度不同，并且误差具有相关性.

表 1 实地观测点坐标及采用的3种随机模型

Table 1. Coordinate pairs of field surveying data and three stochastic models for fitting

点号	x/m	y/m	C 中非零元素			P₁中非零元素			P₂中非零元素	P₃中非零元素
点号	x/m	y/m	$\sigma _x^2$ /mm²	$\sigma _y^2$ /mm²	${\sigma _{xy}}$ /mm²	p_x	p_y	p_xy	${p_x}/{p_y}$	${p_x}/{p_y}$
1	688.639	1398.869	7.3371	9.7881	−0.8902	0.1378	0.1033	0.0125	1	10000
2	701.467	1383.525	7.3289	9.3014	−0.9080	0.1381	0.1088	0.0135	1	10000
3	714.294	1368.180	7.3340	8.8888	−0.9115	0.1381	0.1140	0.0142	1	10000
4	727.121	1352.835	7.3547	8.5485	−0.9080	0.1378	0.1185	0.0146	1	10000
5	739.953	1337.495	7.3948	8.2767	−0.9044	0.1371	0.1225	0.0150	1	10000
6	752.783	1322.152	7.4602	8.0684	−0.9071	0.1359	0.1257	0.0153	1	10000
7	765.609	1306.806	7.5584	7.9166	−0.9210	0.1342	0.1281	0.0156	1	10000
8	778.434	1291.460	7.6977	7.8126	−0.9487	0.1319	0.1299	0.0160	1	10000
9	791.262	1276.115	7.8867	7.7477	−0.9909	0.1289	0.1312	0.0165	1	10000
10	804.088	1260.770	8.1340	7.7132	−1.0462	0.1251	0.1320	0.0170	1	10000
11	816.915	1245.425	8.4471	7.7012	−1.1111	0.1207	0.1324	0.0174	1	10000
12	829.740	1230.078	8.8323	7.7054	−1.1805	0.1156	0.1325	0.0177	1	10000
13	842.564	1214.731	9.2939	7.7206	−1.2485	0.1100	0.1324	0.0178	1	10000
14	855.389	1199.384	9.8346	7.7437	−1.3087	0.1040	0.1321	0.0176	1	10000
15	868.213	1184.036	10.4556	7.7729	−1.3546	0.0979	0.1316	0.0171	1	10000
16	881.043	1168.694	11.1562	7.8077	−1.3803	0.0916	0.1309	0.0162	1	10000
17	893.875	1153.353	11.9355	7.8489	−1.3806	0.0855	0.1301	0.0150	1	10000
18	906.703	1138.009	12.7917	7.8981	−1.3511	0.0796	0.1289	0.0136	1	10000
19	919.537	1122.670	13.7218	7.9569	−1.2880	0.0740	0.1276	0.0120	1	10000
20	932.371	1107.331	14.7236	8.0277	−1.1886	0.0687	0.1261	0.0102	1	10000

下载: 导出CSV

| 显示表格

为验证直线的拟合方法的通用性，选用了对应正交最小二乘和普通最小二乘的两种随机模型对比计算，随机模型P₂为单位矩阵，表示观测点的x坐标和y坐标具有相同的精度，即正交最小二乘拟合；模型P₃表示x坐标的精度远高于y坐标，即普通最小二乘拟合.

由设计资料知该直线与圆曲线直接相接，圆心C的坐标为（1693.970，1443.300），半径R = 800 m. 为测试高斯-牛顿算法的效率和稳健性，选择的初始直线远离最优位置，初始直线参数列于表2中的第0次迭代，闭合差w为圆心到直线的距离减去半径，初始直线与圆曲线相离479.525 m. 算法经过6次迭代收敛，w变为0，说明拟合直线满足约束，与圆曲线相切于直圆点（Z），图4（a）展示了该拟合过程，验证了方法的正确性.

表 2 顾及约束和相关误差的直线拟合过程

Table 2. Process of straight-line fitting with both a constraint and correlated noise

迭代数/次	$\hat a$	$\hat b$ /m	${\hat \sigma _0}$ /mm	${\sigma _a}$	${\sigma _b}$ /mm	w/mm
0	−0.681654264	2210.153480	59093.318880	6.14777 × 10⁻²	93044.287650	−479524.603400
1	−0.935165106	2013.498778	5844.096808	6.07990 × 10⁻³	9201.714082	−59428.755230
2	−1.146887043	2183.687057	448.743795	9.66425 × 10⁻⁴	1148.244207	−9789.990571
3	−1.195130784	2221.779805	10.362946	3.19907 × 10⁻⁵	35.137584	−394.045005
4	−1.197201392	2223.403509	16.492003	5.52543 × 10⁻⁵	59.714483	−0.691388
5	−1.197204886	2223.406205	16.525002	5.55595 × 10⁻⁵	60.003229	−0.000002
6	−1.197204886	2223.406205	16.525002	5.55598 × 10⁻⁵	60.003515	0

下载: 导出CSV

| 显示表格

图 4 约束及无约束的直线拟合

Figure 4. Straight-line fitting with a constraint and without constraint

下载: 全尺寸图片幻灯片

3种随机模型的非零元素见表1. 约束下3种随机模型的拟合直线参数和精度列于，算法均经过6次迭代收敛，耗时均在0.6 s以内. 基于P₁得到的后验单位权中误差 ${\hat \sigma _0}$ = 16.5 mm，截距的中误差 ${\sigma _b}$ = 60.0 mm，均为最小，相对于P₂和P₃得到的 ${\sigma _b}$ = 66.1 mm，精度提升了9.2%. 3种随机模型下的直圆点坐标（x_Z，y_Z）列于表3，对应位置如图4（b），分别表示为Z₁、Z₂和Z₃，其中Z₂和Z₃点重合.

表 3 约束下3种随机模型拟合直线的参数估值及其精度

Table 3. Parameter estimation of fitting line and their precisions of three stochastic models with constraints

随机模型	$\hat a$	$\hat b$ /m	${\hat \sigma _0}$ /mm	${\sigma _a}$	${\sigma _b}$ /mm	x_Z/m	y_Z/m	迭代数/次	耗时/s
P₁	−1.19720489	2223.4062	16.525	5.55598 × 10⁻⁵	60.0	1079.9809	930.4478	6	0.494
P₂	−1.19722236	2223.4251	49.031	6.12012 × 10⁻⁵	66.1	1079.9772	930.4522	6	0.503
P₃	−1.19722233	2223.4250	76.478	6.12012 × 10⁻⁵	66.1	1079.9772	930.4522	6	0.496

下载: 导出CSV

| 显示表格

令B为空矩阵，用同样的参数初值进行无约束的直线拟合，3种随机模型的拟合直线参数和精度列于表4，算法均经过3次迭代收敛，耗时均在0.5 s以内，表明无约束拟合直线具有更高的效率. 且无约束时的单位权中误差及参数精度均高于约束下的直线拟合，说明在无约束情况下，测点与拟合直线的贴合度更好. 但是，无约束拟合获得的直线不相切于已知圆曲线，如图4（b）所示，相离距离d = 0.18 m，不符合线路连续性的要求. 另外，在无约束拟合时，基于P₁得到的后验单位权中误差 ${\hat \sigma _0}$ = 2.3 mm，截距的中误差 ${\sigma _b}$ = 25.0 mm，相对P₂和P₃精度提升了2.7%. 说明不论是否有约束，考虑坐标相关误差时，获得的直线参数精度均为最高.

表 4 无约束下3种随机模型拟合直线的参数估值及其精度

Table 4. Parameter estimation of line fitting and their precisions of three stochastic models without constraint

随机模型	$\hat a$	$\hat b$ /m	${\hat \sigma _0}$ /mm	${\sigma _a}$	${\sigma _b}$ /mm	迭代数/次	耗时/s
P₁	−1.196269322	2222.672000	2.293655	3.10233 × 10⁻⁵	25.008852	3	0.461
P₂	−1.196259074	2222.663879	6.710739	3.16316 × 10⁻⁵	25.743943	3	0.414
P₃	−1.196259065	2222.663872	0.462488	3.16316 × 10⁻⁵	25.743942	3	0.438

下载: 导出CSV

| 显示表格

4. 结　论

在铁路维护中，全站仪获取的轨道坐标点具有误差相关性，普通最小二乘或正交最小二乘拟合直线不能考虑观测值之间的这种误差相关性. 提出的直线拟合通用模型，可以考虑坐标分量间的误差相关性实现直线拟合. 通过指定随机模型，可以实现普通最小二乘、正交最小二乘或加权整体最小二乘直线拟合，并揭示了其对应的几何意义.

通常的直线拟合方法未考虑约束条件，线路重构中直线的拟合受到相邻线元的约束. 提出的直线拟合通用模型，可以在约束下同时顾及误差相关性实现直线拟合. 顾及坐标相关误差时，可以提升估计参数的精度：约束及无约束下参数估计精度分别提高了9.2%和2.7%.

采用的高斯-牛顿算法运行效率高，能够快速得到直线参数的最佳估值及其精度. 在约束及无约束情况下分别仅6次及3次迭代就搜索出最优直线.

致谢：中南大学土木工程国家级实验教学示范中心项目（201905406）.

图 1 边坡模型

Figure 1. Slope model

下载: 全尺寸图片幻灯片

图 2 速度矢量图

Figure 2. Diagram of velocity vectors

下载: 全尺寸图片幻灯片

图 3 力学模型

Figure 3. Mechanical model

下载: 全尺寸图片幻灯片

图 4 f_a对F的影响

Figure 4. Effect of f_a on F

下载: 全尺寸图片幻灯片

图 5 k_p对F的影响

Figure 5. Effect of k_p on F

下载: 全尺寸图片幻灯片

图 6 f对F的影响

Figure 6. Effect of f on F

下载: 全尺寸图片幻灯片

图 7 软弱夹层形状对F的影响

Figure 7. Effect of slip zone shapes on F

下载: 全尺寸图片幻灯片

图 8 c和F的关系曲线

Figure 8. Ralation curves of c and F

下载: 全尺寸图片幻灯片

表 1 计算结果

Table 1. Calculation results

工况	k_P	d₁ /m	d₂ /m	h₁ /m	h₂ /m	γ /（kN•m⁻³）	c /kPa	φ /（^o）	安全系数
工况	k_P	d₁ /m	d₂ /m	h₁ /m	h₂ /m	γ /（kN•m⁻³）	c /kPa	φ /（^o）	本文方法	摩根斯坦-普莱斯法
1	0	11.59	4.95	3.11	4.95	20	6	30	2.02	2.00
2	0.2	11.59	4.95	3.11	4.95	20	6	30	1.35	1.32
3	0.4	11.59	4.95	3.11	4.95	20	6	30	0.99	0.97
4	0.2	11.59	4.95	3.11	4.95	20	6	20	1.05	1.03
5	0.2	9.78	3.54	2.08	3.54	20	6	20	1.22	1.19
注：表中本文方法计算结果为时间 t 内的最小值.

下载: 导出CSV

表 2 计算结果

Table 2. Calculation results

k_P	d₁ /m	d₂ /m	h₁ /m	h₂ /m	γ /（kN•m⁻³）	c /kPa	φ /（^o）	安全系数
k_P	d₁ /m	d₂ /m	h₁ /m	h₂ /m	γ /（kN•m⁻³）	c /kPa	φ /（^o）	本文方法	二楔块法^[16]	有限元
0	53.56	90.18	1.03	30.00	10.2	0	17	1.36	1.35	1.39

下载: 导出CSV

参考文献(17)

[1]	刘汉香,许强,周飞,等. 含软弱夹层斜坡地震动力响应特性的振动台试验研究[J]. 岩石力学与工程学报,2015,34(5): 994-1005. LIU Hanxiang, XU Qiang, ZHOU Fei, et al. Shaking table test for seismic responses of slopes with a weak interlayer[J]. Chinese Journal of Rock Mechanics and Engineering, 2015, 34(5): 994-1005.
[2]	闫孔明,刘飞成,朱崇浩,等. 地震作用下含倾斜软弱夹层斜坡场地的动力响应特性研究[J]. 岩石力学与工程学报,2017,36(11): 2686-2698. YAN Kongming, LIU Feicheng, ZHU Chonghao, et al. Dynamic responses of slopes with intercalated soft layers under seismic excitations[J]. Chinese Journal of Rock Mechanics and Engineering, 2017, 36(11): 2686-2698.
[3]	唐云波,刘炎,相晨琳. 汶川地区含软弱层岩质边坡地震响应研究[J]. 路基工程,2018(6): 114-118. TANG Yunbo, LIU Yan, XIANG Chenlin. Research on seismic response of rock slope with weak layer in Wenchuan region[J]. Subgrade Engineering, 2018(6): 114-118.
[4]	刘晋南, 蒋鑫, 邱延峻, 等. 水平地震荷载作用下斜坡路基动力稳定性分析. 中外公路, 2010, 30（5）: 42-45. LIU Jinnan, JIANG Xin, QIU Yanjun. Analysis of dynamic stability of subgrade slope under horizontal earthquake load [J]. Journal of China & Foreign Highway, 2010, 30(5): 42-45.
[5]	CHEN W F. Limit analysis and soil plasticity[M]. Amsterdam: Elsevier Scientific Company, 1975.
[6]	XU P, HATAMI K, JIANG G. Seismic rotational stability analysis of reinforced soil retaining walls[J]. Computers and Geotechnics, 2020, 118(4): 103297.1-103297.10.
[7]	MICHALOWSKI R L. Limit analysis and stability charts for 3D slope failures[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2010, 136(4): 583-593. doi: 10.1061/(ASCE)GT.1943-5606.0000251
[8]	何思明,张晓曦,吴永. 基于上限定理的边坡潜在破裂面确定方法与稳定性判识研究[J]. 岩土力学,2012,33(1): 162-166. doi: 10.3969/j.issn.1000-7598.2012.01.025 HE Siming, ZHANG Xiaoxi, WU Yong. Study of stability discriminant of slope and position determination of potential sliding surface based on upper bound theorem[J]. Rock and Soil Mechanics, 2012, 33(1): 162-166. doi: 10.3969/j.issn.1000-7598.2012.01.025
[9]	UTILI S. Investigation by limit analysis on the stability of slopes with cracks[J]. Géotechnique, 2013, 63(2): 140-154.
[10]	XU P, HATAMI K, JIANG G. Seismic sliding stability analysis of reinforced soil retaining walls. Geosynthetics International, 2019, 26(5), 485-496.
[11]	CHOUDHURY D, NIMBALKAR S. Seismic passive resistance by pseudo-dynamic method[J]. Geotechnique, 2006, 55(7): 517-520.
[12]	BASHA B M, BABU G L S. Seismic reliability assessment of external stability of reinforced soil walls using pseudo-dynamic method[J]. Geosynthetics International, 2009, 16(3): 197-215. doi: 10.1680/gein.2009.16.3.197
[13]	GHOSH P G. Seismic active earth pressure behind a nonvertical retaining wall using pseudo-dynamic analysis[J]. Canadian Geotechnical Journal, 2008, 45(1): 117-123. doi: 10.1139/T07-071
[14]	钱家欢. 土工原理与计算[M]. 北京: 中国水利水电出版社, 1996.
[15]	中华人民共和国住房和城乡建设部. 建筑边坡工程技术规范: GB 50330—2013[S]. 北京: 中国建筑工业出版社, 2014.
[16]	QIAN X D, KOERNER R M, GRAY D H. Translational failure analysis of landfills[J]. Journal of Geotechnical and Geoenvironmental Engineering, 2003, 129(6): 506-519. doi: 10.1061/(ASCE)1090-0241(2003)129:6(506)
[17]	周炜,李海波,刘亚群,等. 地震作用下顺层岩质边坡锚固特性的拟动力分析[J]. 岩石力学与工程学报,2016,35(增刊2): 3570-3576. ZHOU Wei, LI Haibo, LIU Yaqun, et al. Pseudo-dynamic analysis of anchored characteristics of layered rock slopes subjected to seismic loads[J]. Chinese Journal of Rock Mechanics and Engineering, 2016, 35(S2): 3570-3576.

施引文献

期刊类型引用(8)

1.	吴功勇，聂兴信，李宗利，侯展娜，江松，阮顺领，赵林海. 岩-土复合边坡研究现状、问题及展望. 有色金属(矿山部分). 2025(01): 1-12+37 . 百度学术
2.	郑晋溪. 顺层软弱夹层影响下的边坡稳定性分析与加固研究. 福建建设科技. 2025(02): 40-43 . 百度学术
3.	张标，瞿凡，蒋毅. 饱水软化效应下砂质泥岩边坡三维稳定性上限分析. 土木工程学报. 2025(04): 108-123 . 百度学术
4.	袁家好，鲁祖德，陈从新，孙朝燚，范凯，刘才华. 含双软弱夹层岩质边坡稳定性上限解. 岩石力学与工程学报. 2024(02): 412-423 . 百度学术
5.	叶鲁青，王雁杰. 帮坡角变化对露天煤矿帮坡稳定状态和最终境界的影响. 矿业研究与开发. 2024(07): 57-63 . 百度学术
6.	陈伟明，郁金平. 含泥质粉砂夹层的土质边坡稳定性分析. 交通科学与工程. 2024(04): 57-64+92 . 百度学术
7.	解则斌，苏子龙. 基于GeoStudio的海螺峪隧道边坡稳定性影响因素的分析. 四川建筑. 2024(04): 191-193 . 百度学术
8.	严彬华，王万平，刘瑞辉，曹校勇. 风积沙地层隧道洞口边坡稳定性分析. 公路. 2022(05): 267-271 . 百度学术

其他类型引用(17)

附加材料(0)

访问统计

点击查看大图

图(8) / 表(2)

计量

文章访问数: 386
HTML全文浏览量: 181
PDF下载量: 40
被引次数: 25

1. 纵、横坐标的误差相关性
2. 约束下顾及误差相关性的直线拟合
2.1 直线拟合模型
2.2 参数精度
2.3 迭代寻优过程
3. 试验算例
4. 结　论

极限分析法求解含软弱夹层边坡稳定性

doi: 10.3969/j.issn.0258-2724.20200156

作者简介: 徐鹏（1988—），男，讲师，研究方向为道路与铁道工程，E-mail：bk20090201@my.swjtu.edu.cn

通讯作者: 李婷（1992—），女，讲师，研究方向为地基变形控制与路基工程等，E-mail：liting9225@126.com

计量

出版历程

Stability Analysis of Slopes with Weak Layers Using Limit Analysis Method

1. 纵、横坐标的误差相关性

2. 约束下顾及误差相关性的直线拟合

2.1 直线拟合模型

2.2 参数精度

2.3 迭代寻优过程

3. 试验算例

4. 结 论

期刊类型引用(8)

其他类型引用(17)

计量

出版历程

目录

1. 纵、横坐标的误差相关性

2. 约束下顾及误差相关性的直线拟合

2.1 直线拟合模型

2.2 参数精度

2.3 迭代寻优过程

3. 试验算例

4. 结 论

作者简介:
徐鹏（1988—），男，讲师，研究方向为道路与铁道工程，E-mail：bk20090201@my.swjtu.edu.cn

通讯作者:
李婷（1992—），女，讲师，研究方向为地基变形控制与路基工程等，E-mail：liting9225@126.com

4. 结　论

4. 结　论