滑坡是黄土区域最常见的地质灾害之一，因其强的破坏力给人民带来较大的生命威胁和财产损失，据了解山西省吕梁地区属滑坡易发区，龙建辉等[1]总结了吕梁山典型顺层岩质滑坡的滑动模式和破坏特征，除此之外土质滑坡更为常见，因此对黄土滑坡预测的研究有很重要的实际意义。至今，关于滑坡预测的研究，物理模拟以及经验预测模型[2-5]为大多数学者使用，除此之外数值模拟方法[6-8]也得到广泛关注，建立监测系统[9-11]则是更直接的手段。由于这些方法不是复杂就是资金消耗，因此广泛应用较为困难。实际上，滑带土的临界状态线能够直观反映其破坏时的应力状态，为预测滑坡发生提供判据。
土的临界状态[12]是指土体在受连续剪切作用下，其应力水平及比容v保持恒定，而土体依然保持塑性流动的状态。通常在p′-q(有效力-偏应力)和v-ln p两个平面内对土的临界状态线进行研究。1990年WOOD提出在p′-q平面内，偏应力-应变图中斜率为0的点为临界状态点，这些点的连线就是临界状态线(CSL线)[13]。
为了能更清晰的认识CSL线，众多学者对黏性土CSL线的影响因素及形态特征进行了探究。戴福初等[14]对松散土体展开了一系列等压和偏压三轴试验，得到土的CSL线，并将土体分为剪缩和剪胀型。临界状态土力学[12]认为应力路径对土的CSL线无影响，而黄质宏等[15]对红黏土进行固结不排水剪试验后得出不同应力路径下其CSL线不唯一。陈能[16]更是对带根重塑黏土在不同应力路径下的CSL线进行了探究，发现幂函数型比直线型更适用。此外，CSL线与含水率的关系也很密切。卞夏等[17]探究黏性土含水率对CSL线的影响，发现初始含水率越小，临界状态线在p′-q坐标系内越靠近上方。此外，李涛等[18]在已构建的二维形式弹塑性双面模型的基础上，构建了三维形式，其适用于饱和黏性土且土体临界状态也被体现。
上述CSL线都是针对黏性土建立的，部分学者也针对砂性土有所探讨。LI等[19]认为，用指数形式来拟合砂土的CSL线更好。张敏等[20]总结出密砂的CSL线是惟一的，只有当多应力路径中存在包含主应力旋转的路径时其临界状态线才会由线性过渡到非线性。而侯悦琪[21]通过对前人研究结果进行归纳，总结出砂土材料的临界状态线不唯一。刘超[22]探讨了非饱和粉砂在达到临界状态时各相关状态参量需达到的条件，并给出不同坐标平面上的临界状态线及其方程。
目前普遍认为滑坡的发生与抗剪强度参数不断弱化存在直接关系，即土体黏聚力c与体内摩擦角φ都会随着滑带土含水量的增多而降低[23-24]，但是直接通过滑带土强度参数的弱化来预测滑坡发生非常困难。因此提出一种新思路，运用土的临界状态线来预测滑坡发生。
通过对正城滑坡重塑滑带土在10组含水率下的CU试验，绘制各含水率下的CSL线。在推导出CSL线与抗剪强度参数相关关系式的基础上，结合滑带土强度参数c、φ与含水率的变化关系式，可以得到CSL线与含水率的数量关系。所以当得知某一时刻含水率下滑带土的临界状态线，再与该时刻土体应力状态进行对比，能够判断滑带土是否破坏，进而为预测滑坡提供判据。最后，用柳林县贺西滑坡来说明临界状态线在预测滑坡中的可行性，并得出了滑坡发生时，滑带上超过对应CSL线的点数占比。
选用的滑带土取自山西省孝义市正城滑坡。根据野外勘察资料确定，该滑坡主要由降雪融水作用引发，属于典型的由水诱发滑坡[25]。
在勘察阶段，通过探井在处于6 m深处的滑带上取得滑带土试样。经观察可知，其主要由Q2的粉质黏土组成，由于水的入渗作用导致天然水分含量较高，相关物理参数见表1。
表1 滑带土的基本物理参数
Table 1 Basic physical properties of slip zone soil
基本物理参数数值基本物理参数数值天然含水率/%21.21塑限/%22.87湿密度/(g·m-3)1.78液限/%35.54干密度/(g·m-3)1.47塑性指数/%12.67比重2.72粉粒含量/%68.12孔隙比0.682黏粒含量/%30.88
由于滑坡是由融雪入渗引发的，在滑坡发生之际坡体内水分来不及排出，且所取滑带土样原始结构基本被破坏，所以采用TFB-1型三轴剪切试验仪进行重塑滑带土试样的CU(固结不排水剪)试验。对土样处理后制备低于塑限的试样6组(土壤含水率w为10.85%、12.27%、15.28%、17.95%、20.25%、22.56%、)以及高于塑限却低于液限的试样4组(土壤含水率w为23.50%、25.32%、27.55%和30.24%)。之后在100、200、300、400 kPa围压下进行试验。
经处理后，绘制10组含水率下的偏应力-应变曲线，篇幅限制只展示22.56%和30.24%两种含水率如图1所示，图2为围压400 kPa下不同含水率滑带土试样应力-应变(q-ε)关系。
图1 2种含水率下滑带土试样的应力-应变曲线
Fig.1 Stress-strain curves of slip zone soil samples under two different water contents
图2 围压为400 kPa的应力-应变曲线
Fig.2 Stress-strain curve of slip zone soil(σ3=400 kPa)
分析图1得，开始时偏应力迅速攀升直至应变达3×10-2速度才有所减小，但是在持续增加，且在同一含水率下，围压增加，滑带土的破坏应力就会随之增大，破坏呈现典型的应变硬化，并且高含水率比低含水率的应变硬化程度要低。出现应变硬化的原因是在剪切发生时土体原有的孔隙被挤压出现减缩现象，许多颗粒相接触使得摩擦力增大进而提高了剪切强度，发生应变硬化[26]。图2表明在同一围压下，随着含水率的增加，滑带土达到破坏时的偏应力先增大后减小，在含水率为15.28%时所需的破坏偏应力最大，即最难破坏。在低应变范围内，含水率较低的土样应力-应变曲线的斜率较高，即应力增长较快。
对于应变硬化型土体的CSL线的确定，WOOD[13]指出可在应力—应变图中找到15%的轴向应变所对应的偏应力即为临界状态点。结合三轴试验结果，将某一含水率下4个围压的临界状态点连起来生成的线即为CSL线，图3为含水率为22.56%时的CSL线图。可以看出，土的临界状态线在p-q(总应力-偏应力)平面内为一条斜直线。
图3 含水率为22.56%时滑带土的临界状态线
Fig.3 Critical state line of slip zone soil samples(w=22.56%)
为建立土的临界状态线与对应含水率间的数量关系，将抗剪强度参数作为桥梁，首先推导出临界状态线与强度参数间的关系式。Mohr-Coulomb破坏理论[27]给出土体达到破坏状态时各主应力及强度参数的几何关系如图4所示。
a—应力莫尔圆的圆心坐标；b—应力莫尔圆的半径；σ1—最大主应力(kPa)；σ3—最小主应力(kPa)；α—任一截面与最小主应力作用方向的夹角(°)；c—土的黏聚力(kPa)；φ—土的内摩擦角(°)
图4 土体极限平衡条件[27]
Fig.4 Limit equilibrium condition of soil[27]
对图4中各个参数的几何关系进行推导可以得到p和q的关系为
(1)
由(1)式可看出，图4的p-q斜直线的截距和斜率均与c和φ有关。
图5为各含水率对应的10条临界状态线。明显观察出在开始时随含水率的小幅增加临界状态线也在上移，在含水率到达15.28%后，随着含水率的增加CSL线逐渐下移，可见含水率在15.28%附近有使得CSL线位于最高的含水率。整体来看，含水率为15.28%试样的CSL线最高，含水率为30.24%试样的CSL线最低，其余含水率对应的CSL线基本都位于这两条之间。这说明含水率在15.28%之后，高含水率水平更容易达到临界状态，即滑带土含有的水分越多越易被破坏，滑坡发生的可能性也越大。同时，CSL线的斜率随含水率的增大而减小，这说明高含水率的土样其破坏偏应力增长较慢，滑带土不需要太大的应力作用便能达到破坏应变量。总体来说，CSL线主要受含水率的控制。
图5 各含水率下的10条临界状态线
Fig.5 Ten critical state lines under different water contents
对各条临界状态线进行拟合，得到相应的q-p关系式，进而得到正城滑坡滑带土p-q曲线与抗剪强度参数的数量关系。
根据试验得到的滑带土偏应力-应变曲线和Mohr-Coulomb破坏理论，在τ-σ(剪应力-法向应力)坐标系下绘制应力莫尔圆，得到各含水率下应力圆公切线的斜率和截距并归纳出滑带土样的φ及c随含水率的变化规律，相应的曲线如图6所示。
图6 抗剪强度参数随含水率的变化规律
Fig.6 Relationship between shear strength parameters and watercontents
可见，c随着含水率的增加先增大后减小，在15.28%时出现最大值134.72 kPa；φ随含水率的增加整体呈下降趋势，并且在到达15.28%含水率前φ陡然下降，之后基本持平但依然有降低趋势，而后在26%左右又开始显著下降。相比之下，c对含水率更敏感，说明含水率的变化主要影响了土颗粒间的黏聚力，尤其是在塑限附近，两参数值降低更明显。
根据室内土工试验结果，可知滑带土的含水率基本分布在15.28%～27.55%(接近饱和)，由图6可看出，含水率在15.28%～27.55%与抗剪强度参数近乎呈线性关系，因此使用最小二乘法，拟合含水率w与
以及c的关系曲线，拟合后的表达式如下：
(2)
c=293.67-10.08w
(3)
将式(2)、式(3)代入式(1)，可推导出滑带土的临界状态线与含水率的关系：
(4)
图7中，表示出含水率为15.28%和27.55%所对应的临界状态线及其表达式，并用点线箭头示意出其随含水率的变化趋势。可以看出，随着含水率的增加，临界状态线逐渐下移并且斜率也在减小，与前述得出的规律完全一致。
图7 临界状态线随含水率的变化趋势
Fig.7 Trend chart of critical state line with water content
由于地域气候等差异，不同滑带土所表现出来的物理力学性质不尽相同，使得滑带土的抗剪强度对水的敏感性也不同，因此其CSL线随含水率的变化会有各种表现。但是从整体趋势来看，必然是含水率越高，土体越容易达到该含水率下的CSL线，即滑带土越容易破坏。
为说明如何应用临界状态线预测滑坡失稳，以山西省柳林县贺西煤矿选煤厂的贺西滑坡为例，选取新滑坡进行讨论，结合相应的勘察资料和室内实验数据，对滑带土临界状态线在滑坡破坏与预测中的应用进行阐明。
首先，根据贺西滑坡滑带土的室内土工试验数据，可知该滑坡的启滑含水率，进而应用第3节所述思想，推导出该含水率下的CSL线公式。该临界状态线即表示滑坡破坏时的临界状态；然后，对滑动面土体在临滑时所受的实际应力大小进行计算，并与启滑含水率下的CSL线进行比较，可以得到在滑坡发生时滑带上有多少比例的土体达到CSL线，从而为类似滑坡预测提供依据。
滑坡体主要由黄土状粉土及粉质黏土堆积而成，根据探井所取得的土样可知滑带土为粉质黏土与红黏土叠加在一起。由贺西滑坡的室内土工试验可知，新滑坡滑带的启滑含水率为17.7%，对应的黏聚力为18.3 kPa，内摩擦角为24.8°，代入式(1)推导出17.7%含水率下临界状态线表达式为
q=0.98p+38.63
(5)
土体的应力状态主要是指土体的平均主应力p及剪应力q，而这两者都与土体所受的上覆土压力σ1及周围压力σ3有关。因此，根据上覆土体的土性和厚度计算出滑带土上覆土压力及周围压力进而表征出滑带土颗粒在临滑时的应力状态。
为保证计算结果的准确性及代表性，将恢复原始地形的剖面图作为计算剖面，在滑带上每隔5 m的水平间距点取一定数量的土颗粒研究点，确定上覆土体土性并量取厚度。本次选取主剖面3—3′进行计算，共计25个土体研究点，如图8所示。
图8 3—3′计算剖面
Fig.8 3—3′ calculation profile
计算上覆土压力及周围压力的公式如下：
σ1=γH
(6)
σ3=ξγH
(7)
式中：γ为上覆土体重度，kN/m3；ξ为土的侧压力系数，无量纲；H为上覆土体厚度，m。
根据勘察探井土样的室内土工试验结果，滑坡体的天然容重为19.5 kN/m3；根据土的经验数据，坚硬状态的粉质黏土侧压力系数取0.33[28]，上覆土体厚度可在剖面图上量取。然后利用公式p=(σ1+σ2+σ3)/3,q=σ1-σ3计算出每个土颗粒所受的应力大小并在图9中表示出来。
由图9可知，3—3′剖面上超过临界状态线的点有15个，占比为60%。由此可得，当滑带土的监测点的应力状态超过临界状态线的比例达到60%，滑坡即有可能发生。由此，如果预先判断出滑坡可能发生破坏的位置，并在该处布置水分计和土压力盒以监测含水量及土体应力的变化，然后通过对某一含水率下滑带土的临界状态线与该含水率下滑带土的应力状态进行比较，可以为滑坡启动破坏提供依据。
图9 滑带土的临界状态线与应力状态的比较
Fig.9 Comparison of critical state line and stress state of slip zone soil
1)试样破坏呈应变硬化，在3%应变内破坏应力会迅速增大，且黏聚力和偏应力均随含水率先增大后减小，都在15.28%时达到顶峰，c的最大值为134.72 kPa；而随含水率的增大，φ呈下降趋势，而且c对水更敏感，尤其是在塑限附近降低的更为明显。
2)在p-q平面内，滑带土的临界状态线在各含水率下表现为不同的斜直线，意味着滑带土在不同含水情况下能被破坏的极限状态不同，在含水率超过15.28%后，高含水率水平更容易达到临界状态。以抗剪强度参数为桥梁，推导出临界状态线与含水率的数量关系。
3)以贺西滑坡为实例说明临界状态线在滑坡预测中的应用，并得出滑坡发生时，滑带上超过临界状态线的点数占比约60%。因此，若能预先判断出坡体在哪处发生滑动，便可以取土展开剪切试验得到各含水率下的临界状态线，并在该处布置水分计和土压力盒以监测含水量及土体应力，然后通过对实时含水率下滑带土的临界状态线与该含水率下滑带土的应力状态进行比较，得到发生破坏的土体所占比例，从而为滑坡预测提供预判依据。
[1] 龙建辉,任 杰,曾凡桂,等.双软弱夹层岩质滑坡的滑动模式及变形规律[J].煤炭学报,2019,44(10):3031-3040.
LONG Jianhui，REN Jie，ZENG Fangui，et al.Sliding mode and deformation law of double weak interlayer rock landslide[J].Journal of China Coal Society,2019,44(10):3031-3040.
[2] 殷坤龙,晏同珍.滑坡预测及相关模型[J].岩石力学与工程学报,1996,15(1):1-8.
YIN Kunlun，YAN Tongzhen.Landslide prediction and relevant models[J].Chinese Journal of Rock Mechanics and Engineering,1996,15(1):1-8.
[3] 黄润秋,许 强.斜坡失稳时间的协同预测模型[J].山地研究,1997,15(1):7-12.
HUANG Runqiu，XU Qiang.Synergetic time-prediction model of slope instability[J].Journal of Mountain Research,1997,15(1):7-12.
[4] 王建锋.两类经典滑坡发生时间预报模型的理论分析[J].地质力学学报,2004,10(1):40-50.
WANG Jianfeng.Theoretical analysis to two classic landslide prediction models:Saito’s Model and Voight’s Model[J].Journal of Geomechnicas,2004,10(1):40-50.
[5] 李秀珍,孔纪名,王成华.中心逼近式灰色GM(1,1)模型在滑坡变形预测中的应用[J].工程地质学报,2007,15(5):673-676.
LI Xiuzhen，KONG Jiming，WANG Chenghua.Application of center approach Gray GM(1,1) Model to prediction of landslide deformation with a case study[J].Journal of Engineering Geology,2007,15(5):673-676.
[6] 朱文彬,刘宝琛.降雨条件下土体滑坡的有限元数值分析[J].岩石力学与工程学报,2002,21(4):509-512.
ZHU Wenbin，LIU Baochen.Forming and development process of soil landslide during rainfall[J].Chinese Journal of Rock Mechanics and Engineering，2002,21(4):509-512.
[7] KIM D,IM S,SANG H L,et al.Predicting the rainfall-triggered landslides in a forested mountain region using TRIGRS model[J].Journal of Mountain Science (in English)，2010,7(1):83-91.
[8] 孙玉进,宋二祥.大位移滑坡形态的物质点法模拟[J].岩土工程学报,2015,37(7):1218-1225.
SUN Yujin，SONG Erxiang.Simulation of large-displacement landslide by material point method[J].Chinese Journal of Geotechnical Engineering,2015,37(7):1218-1225.
[9] KUMAR K V,LAKHERA R C,MARTHA T R,et al.Analysis of the 2003 Varunawat Landslide,Uttarkashi,India using Earth Observation data[J].Environmental Geology，2008,55(4):789-799.
[10] 龙建辉，朱清华，倪向龙.地下资源开采诱发的一种滑坡模式研究——以采空触发滑坡为例[J].煤炭科学技术,2021,49(8)：181-187.
LONG Jianhui,ZHU Qinghua,NI Xianglong.Study on landslide model induced by underground resources exploitation:taking mining-triggered landslide as an example[J].Coal Science and Technology，2021，49(8)：181-187.
[11] DONG Jie,ZHANG Lu,LI Menghua,et al.Measuring precursory movements of the recent Xinmo landslide in Mao County,China with Sentinel-1 and ALOS-2 PALSAR-2 datasets[J].Landslides，2018,15(1):135-144.
[12] SCHOFIELD M A,WROTH C P.Critical state of soil mechanics[M].London:McGraw-Hill,1968.
[13] WOOD D M.Soil behaviour and critical state soil mechanics[M]. London:Cambridge University Press,1990.
[14] 戴福初,李焯芬,王思敬.松散击实火山岩坡残积土的应力应变特性及其对滑坡的意义[J].岩土工程学报,1999,21(3):1-6.
DAI Fuchu，Lee Chackfan，WANG Sijing.Stress-strain behavior of a loosely-compacted volcanic-derived soiland its implications to landsliding[J].Chinese Journal of Geotechnical Engineering,1999,21(3):1-6.
[15] 黄质宏,朱立军,廖义玲,等.不同应力路径下红粘土的力学特性[J].岩石力学与工程学报,2004,23(15):2599-2603.
HUANG Zhihong，ZHU Lijun，LIAO Yiling,et al.Mechanical properties of red clay under different stress paths[J].Chinese Journal of Rock Mechanics and Engineering,2004,23(15):2599-2603.
[16] 陈 能.应力路径对根—土复合体力学特性影响试验研究[D].长沙:中南林业科技大学,2015.
CHEN Neng.Experimental study on the effect of stress path on mechanical characteristics of root-soil composites[D].Changsha:Central South University of Forestry & Technology,2015.
[17] 卞 夏,洪振舜,蔡正银,等.重塑黏土临界状态线随初始含水率的变化规律[J].岩土工程学报,2013,35(1):164-169.
BIAN Xia,HONG Zhenshun,CAI Zhengyin,et al.Change of critical state lines of reconstituted clays with initial water contents[J].Chinese Journal of Geotechnical Engineering,2013,35(1):164-169.
[18] 李 涛,Meissner H.循环荷载作用下饱和黏性土的弹塑性双面模型[J].土木工程学报,2006,39(1):92-97.
LI Tao,Meissner H.Elastoplastic two-surface model for clays under undrained cyclic loading[J].China Civil Engineering Journal,2006,39(1):92-97.
[19] LI X S，WANG Y.Linear Representation of Steady-state Line for Sand[J].Journal of Geotechnical and Geo-environmental Engineering，1998,124(12):1215-1217.
[20] 张 敏,杨蕴明,李 琦,等.含主应力旋转的应力路径对密砂临界状态的影响[J].岩石力学与工程学报,2013,32(12)：2560-2565.
ZHANG Min，YANG Yunming，LI Qi，et al.Influence of stress paths including principal stress rotation on critical state of dense sand[J].Chinese Journal of Rock Mechanics and Engineering,2013,32(12):2560-2565.
[21] 侯悦琪.砂土本构关系与LS-DYNA二次开发应用[D].上海:上海交通大学,2011.
HOU Yueqi.Sand constitutive modeling and its application in LS-DYNA[D].Shanghai:Shanghai Jiao Tong University,2011.
[22] 刘 超.非饱和粉砂临界状态试验研究[D].北京:北京交通大学,2016.
LIU Chao.Experimental study of unsaturated silt critical state[D].Beijing:Beijing Jiaotong University,2016.
[23] HAN Bo,LIAO Hongjian,LI Hangzhou,et al.Experimental study on shear strength characteristics and stress-strain relationship of loess[J].Key Engineering Materials，2013,535-536:574-577.
[24] Shi X S,Herle I.Undrained shear strength and water content distribution of remolded clay mixtures[J].Geotechnik，2017,40(1):60-63.
[25] 王晓娅.自然条件下黄土滑坡预测方法的试验研究：以正城滑坡为例[D].太原:太原理工大学,2018.
WANG Xiaoya.Experimental research on prediction methods of loess landslides under natural condition：A case of Zhengcheng Landslide[D].Taiyuan:Taiyuan University of Technology,2018.
[26] 王靖涛.论岩土塑性体应变与剪应变的相互作用原理[J].中国工程科学,2006,8(9):24-29.
WANG Jingtao.Demonstration of the principle of interaction between plastic volume strain and shear for rock and soil[J].Engineering Sciences,2006,8(9):24-29.
[27] 李广信，张丙印，于玉贞.土力学(第2版)[M].北京：清华大学出版社，2013.
[28] ISBN 978-7-112-08828-7.工程地质手册(第四版)[S].北京：中国建筑工业出版社，2007.