1.一种基于牵引式斜坡变形破坏机理的稳定性分析和预测预警方法，其特征在于包含如下步骤：(1)分析牵引式斜坡的变形机理，将斜坡前缘以临界滑面为基准，定义为破坏后区、临界状态面及稳定状态区，牵引式斜坡高度与变形的关系呈现出S曲线特征；分析可能发生破坏的三种情况，情况I为沿最弱的滑面发生破坏，情况II为斜坡发生新的剪切破坏，情况III为坡体或坡体上的节理或裂隙发生拉破坏；

(2)计算获得滑面各点的临界应变值，并计算出现行滑面临界应力状态点及位移值；

(3)针对步骤(1)的三种情况分别计算稳定系数，进行稳定性分析；具体为利用不平衡拉力法计算稳定系数：(3.1)先对选定坡体的滑面，结合坡面特征点决定滑面特征点，对滑面特征点的切线作垂线，交滑体坡面形成不同特征点；

(3.2)连接各特征点形成直线，形成不平衡拉力法的计算条块；

(3.3)不平衡拉力法的假设同不平衡推力法假设；

(3.4)计算时，按条块编号从下往上进行；

(3.5)计算方法为使斜坡体最上面的一个条块推力为零，从而获得不平衡拉力法的各条块下滑力、阻力、剩余拉力和稳定系数；

上述不平衡拉力法的条块划分是为了决定坡体各特征点对应的最小断面，比较各断面所对应的剩余拉力Pi和滑体相应断面对应的最大抗拉力Pit＝σctSd大小，从而决定各断面是否发生拉破坏，当Pi≥Pit时，滑体发生拉破坏；

式中σct,Sd分别为滑体或节理、裂隙的抗拉强度和对应的断面面积；

(4)利用边坡稳定性计算的滑面边界法获取斜坡坡面位移，进行预测预警。

2.如权利要求1所述的基于牵引式斜坡变形破坏机理的稳定性分析和预测预警方法，其特征在于：所述情况II的剪切破坏由滑面的剪应力与剪应变的全过程曲线决定，情况III的拉破坏由坡体的拉应力与拉应变全过程曲线决定。

3.如权利要求1所述的基于牵引式斜坡变形破坏机理的稳定性分析和预测预警方法，其特征在于所述步骤(2)的具体计算过程如下：步骤(2.1)分析斜坡的基本形态、特征，试验获得滑面的基本物理力学参数G、S、m、ρ、C、φ、a1,a2,a3,ξN和坡体的弹性模量力学参数，计算相对应的位移场和应力场，通过应力场、位移场决定相对应的稳定系数；

步骤(2.2)将步骤(2.1)得到的参数代入公式τ＝Gγ[1+γm/S]ρ，式中τ、γ分别表示材料的剪应力和剪应变，τ、G的单位为MPa或kPa或Pa，S、m和ρ为无单位参数，且-1＜ρ≤0和1+mρ≠0；

peak peak

临界应力空间τ 采用摩尔库伦准则τ ＝C+σntanφ，式中C为凝聚力，σn为法向应力，C和σn的单位为MPa或kPa或Pa，φ为滑面摩擦角；

临界应变空间γpeak采用公式 式中σn为法向应力，单位

为MPa或kPa或Pa；

临界应力空间与临界应变空间的关系为τpeak/γpeak＝G[1-1/(1+mρ)]ρ，临界应变空间符合公式S+(1+mρ)γmpeak＝0；

参数 式中ρ0为法向应力σn为零值的ρ值，ρc为σn等于σnc时的ρ值，为常系数。

4.如权利要求1所述的基于牵引式斜坡变形破坏机理的稳定性分析和预测预警方法，其特征在于所述步骤(3)利用综合下滑力—抗滑力法计算稳定系数：对于情况I或情况II选定的滑体的每个条块，对沿滑面的单元的下滑力Psi和临界摩阻力Ti分别求其矢量和，形成综合下滑力Ps和综合摩阻力T，且分别与水平轴形成的最小夹角为αs和αf，定义稳定性系数为：

式(1)的物理意义为单位下滑力所产生的抗力；

对于情况III，假定计算所得e点为临界状态点，选取可能破坏的滑体为研究对象，对于选定滑体的每个条块，对沿滑面的单元的下滑力Psi和可能破坏面的剩余拉力与临界摩阻力Ti和可能破坏面的最大拉力Pit分别求其矢量和，形成综合下滑力Ps和综合抗滑力T，且分别与水平轴形成的最小夹角为αs和αf，定义稳定性系数为式(1)；

对于情况III，以最小稳定系数作为可能发生拉破坏的稳定系数；在此基础上分析下一步可能的破坏形式。

5.如权利要求1所述的基于牵引式斜坡变形破坏机理的稳定性分析和预测预警方法，其特征在于所述步骤(3)利用牵引式斜坡的主拉力法计算稳定系数：对于情况I或情况II选定的滑体的每个条块，m-1条块为临界状态条块，求1～m-1条块滑面的剩余拉力Psi矢量和，形成综合剩余拉力Pm，且与X轴的夹角为αsm；求m至n每个条块的下滑力Psi与临界摩阻力Ti差值的矢量和，形成综合剩余抗滑力Tm-n，其与水平轴的夹角为αfm-n，定义富余稳定性系数为：

式(2)的物理意义为沿主滑方向单位下滑力的剩余抗力；

对于情况III，假定计算所得e点为临界状态点，选取可能破坏的滑体为研究对象，m-1条块为临界状态条块，求1～m-1条块滑面的剩余拉力Psi矢量和，形成综合剩余拉力Pm，且与X轴的夹角为αsm；求m条块至bf断面对应的每个条块的下滑力Psi与临界摩阻力Ti差值，以及bf断面最大抗拉力与此时具有的拉力差值矢量和，形成综合剩余抗滑力Tm-n，其与水平轴的m-n夹角为αf ，定义富余稳定性系数为式(2)。

6.如权利要求1所述的基于牵引式斜坡变形破坏机理的稳定性分析和预测预警方法，其特征在于所述步骤(3)利用牵引式斜坡变形稳定系数法计算稳定系数：定义牵引式边坡现状稳定系数，X轴方向的现状稳定系数为边坡后缘关键块发生破坏时沿滑面在X轴方向的位移矢量和Sc-t除以现状边坡前缘至滑面现状临界状态单元或条块在X轴方向的位移矢量和Sp-t，采用相同的方法定义Y、Z轴方向的现状稳定系数；综合稳定系数定义为边坡后缘关键块发生破坏时沿滑面在X、Y、Z轴方向的位移矢量和除以现状边坡前缘至滑面现状临界状态单元或条块在X、Y、Z轴方向的位移矢量和；存在三个不同方向的稳定系数，X、Y、Z轴三方向的稳定系数分别为Fs-x＝Sc-tx/Sp-tx、Fs-y＝Sc-ty/Sp-ty、Fs-z＝Sc-tz/Sp-tz。

7.如权利要求1所述的基于牵引式斜坡变形破坏机理的稳定性分析和预测预警方法，其特征在于步骤(4)中所述边坡稳定性计算的滑面边界法，对于牵引式斜坡，取滑体作为计算对象，潜在滑动面上的力与位移以边界条件加以处理，滑面以梁单元或其它单元加以处理，基本假设为梁单元沿横向没有变形或变形遵循不同本构关系，坡体采用现行的各种单元；具体计算步骤如下：第一步，选取坡体作为有限元计算对象，除滑面的正应力与剪应力外，坡体表面边界条件计算方法与常规一致，以滑面底边梁单元为第1单元，沿着滑面底边直至滑面顶面的最后一个梁单元依次排序定义为第N单元；

第二步，在滑面垂直方向赋应变为零，并在此基础上进行斜坡的第一次计算，从而获得滑面梁单元底面第一次计算的法向应力、切向应力和应变，如式(3，4)，

式中Tn,si，Sn,si分别表示滑面第i次计算所得的法向应力、切向应力和法向应变、切向应变矩阵，σnm,i，τsm,i，εnm,i和εsm,i分别表示第m单元第i次计算所得的法向应力、切向应力、法向应变和切向应变；n表示法向，s表示切向，i表示计算次数；

比较每一单元第一次计算所得的切向应力与计算所得的法向应力相对应的临界摩阻力的大小及方向，当切向应力与临界摩阻力方向相反时，且切向应力绝对值大于临界摩阻力时，相应单元切向应力作为第二次计算时的边界条件等于临界摩阻力加上该单元计算所得切向力，作为该单元下次计算的初始切向应力边界条件，假定1～k单元切向应力与摩阻力方向相反；对于某单元计算所得切向应力方向与摩阻力方向一致或相反，但绝对值小于临界摩阻力，下次计算时不作为初始条件，假定为k+1～N单元，如式(5)，

第三步，进行第二次计算，相对应1～k及1～N个单元力的初始应力及应变边界条件分别如下式，

按式(6)和式(7)，进行第二次计算所得的滑面单元的法向应力及切向应力矩阵见式(8)，

首先对1～N单元进行判定，按式(8)的法向应力以下式计算对应的临界摩阻力：i,2 i,2

τcrit ＝Ci+σn tanφi i∈(1,N)   (9)式中，τcriti,j为第i单元第j次计算所得临界摩阻力，σni,j为第i单元第j次计算所得法向应力；

对于已经施加了临界摩阻力的滑面对应单元，比较法向应力大小，当此时计算的法向应力之差大于某定值D1时，重新对该单元切向应力按式(11)赋值，|σnL,2-σnL,1|≥D1，L∈(1,k)   (10)ΔτsL,2＝τsL,1+τcritL,1+(τcritL,2-τcritL,1)/2，L∈(1,k)   (11)

对于未赋临界摩阻力的单元，当切向应力与临界摩阻力方向相反时，且绝对值大于临界摩阻力时，对这样的单元施加临界摩阻力如下式(12)，如第i单元第2次计算，施加的切向应力((τcriti,2+τsi,2)，假定i∈(k,m))；

Δτsi,2＝τcriti,2+τsi,2,i∈(k,m)   (12)按照上述计算，在进行第三次计算时，1～N个滑面应力及应变边界条件分别为：

第四步：重复第三步，计算直至多次，对于已经施加了临界摩阻力的滑面对应单元，比较法向应力大小，当此时与上一次计算的法向应力之差最大值小于某定值D1时，则可以结束计算，|σnL,MM+1-σnL,MM|max＜D1,L∈(1,m)   (15)式中，MM为计算次数，m为施加临界摩阻力单元数；

此时记下法向应力σni,MM+1和切向应力τsi,MM+1，其中i∈(1,N)，以及i,MM+1

切向应变εs ，其中i∈(1,N)。