1.一种考虑孔隙率影响的黄土边坡的多场耦合分析方法，包括如下的步骤：步骤S1：有限元软件几何建模；

步骤S2：建立与孔隙率相关的本构模型；

建立固相，及液相与气相，将孔隙率与弹性模量、粘聚力、固有渗透系数相耦合；包括如下：

1）固相；

针对非均质黄土边坡进行力学分析，所考虑的总应变率 分为两部分：弹性（）和塑性（ ）；

；

对于弹性部分， 是弹性刚度矩阵， 是弹性模量， 是泊松比；

对于塑性部分，  是屈服面，  是塑性势, 是流动性， 是应力功率， 是有效应力；

弹性本构是采用具有杨氏模量和泊松比的标准线弹性模型，且弹性模量 随孔隙率变化：；

式中： 为参考弹性模量， 为孔隙率， 为弹性模量参考孔隙率；

塑性本构采用的是改进的Mohr‑Coulomb模型： ；

式中： 是摩擦角， 是剪胀角， 是塑性势参数， 、 、 分别是有效平均应力、偏应力和罗德角；

同时，考虑孔隙率对粘聚力的影响：；

；

；

式中：为粘聚力， 、 为粘聚力参数， 为粘聚力参考孔隙率；

2）液相与气相；

用Van‑Genuchten模型作为持水曲线来描述黄土的吸力变化，其具体表达式为：；

式中： 为实际含水率； 为残余含水率，即土体中含水率的最小值； 为最大含水率，等于土体饱和时的含水率； 为当前含水率； 为孔隙气压， 为孔隙水压，二者的差值就是基质吸力； 为持水曲线的形状控制参数，只与对应的土体相关； 则是土体的进气值，表达式如下：；

；

；

式中： 为 温度下对应的进气值； 为表面张力； 为孔隙率， 为参考孔隙率； 、为相关参数；

达西定律被用来表征液相通量 与气相通量 ，具体表达式为：；

其中： 是固有渗透率； 是相对渗透率； 是黏度； 是两点之间的孔隙水压梯度； 为液体的密度； 是重力加速度矢量； 时，所有参数表示液相参数，当  时，所有参数表示气相参数；

采用固有渗透率的指数定律： ；

其中：  为渗透率参考孔隙率， 是对应的固有渗透率；

液相的相对渗透率 和气相的相对渗透率 定义如下：；

；

所使用的本构模型考虑了孔隙率的影响，将孔隙率与弹性模量、粘聚力、固有渗透系数相耦合；

步骤S3：引入孔隙率随机场和引入大气边界条件；

步骤S31：引入孔隙率随机场表征黄土边坡非均质特性；

空间相关的非均质孔隙度场的生成，并遵循地质统计学理论的半变差函数关系：；

其中， 为块金值； 为基台值； 为相关联的范围； 为两互相影响元素的距离；

进一步转化为空间问题，在二阶平稳性条件下获得空间恒定的均值和方差，其协方差函数和半变差函数之间的关系为 ，得到：；

生成一个孔隙率平均值为 ，协方差为 的孔隙率随机场，其均值和协方差公式为：；

；

向量 做出以下改变，来有效表征空间各向异性：；

；

； ； ；

其中, 、 和 分别是 对应在轴x、y和z的相应范围； 、 、 是向量   的坐标； 、 、 是原始坐标； 、 、 分别是原坐标轴x、y和z旋转至现坐标轴的所转动的角度；

步骤S32：引入大气边界条件模拟复杂的土壤‑大气相互作用；

步骤S4：改进的强度折减法求边坡安全系数；

采用强度折减法来求边坡安全系数，折减系数 的定义为当前实际强度参数值（ 、）与折减后强度（ 、 ）之比，相关参数的关系：；

式中：和 分别为岩土材料的粘聚力和内摩擦角；

将塑性剪切带贯穿作为边坡失稳的判据，计算黄土边坡的安全系数；

步骤S5：得到安全系数。

2.一种计算机系统，包括存储器、处理器及存储在存储器上的计算机程序，其特征在于：所述处理器执行所述计算机程序以实现权利要求1所述多场耦合分析方法的步骤。