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

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

建立固相，及液相与气相，将孔隙率与弹性模量、内聚率、固有渗透系数相耦合；

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

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

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

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

采用强度折减法来求边坡安全系数，折减系数Ft的定义为当前实际强度参数值与折减后强度之比，相关参数的关系：式中：C和分别为岩土材料的内聚力和摩擦角；

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

步骤S5：得到安全系数。

2.根据权利要求1所述的考虑孔隙率影响的黄土边坡的多场耦合分析方法，其特征在于：步骤S31：引入孔隙率随机场表征黄土边坡非均质特性；

空间相关的非均质孔隙度场的生成基于Davis的研究，并遵循地质统计学理论的半变差函数关系：其中，C0为块金值；C0+C1为台基值；a为相关联的范围；|h|为两互相影响元素的距离；

进一步转化为空间问题，在二阶平稳性条件下获得空间恒定的均值和方差，其协方差函数和半变差函数之间的关系为C(h)＝C(0)-γ(h)，得到：生成一个孔隙率平均值为μ，协方差为C0+C1的孔隙率随机场，其均值和协方差公式为：E(φ)＝μ

Cov(φ)＝C0+C1

向量h做出以下改变，来有效表征空间各向异性：其中,ax、ay和az分别是a对应在轴x、y和z的相应范围；xani、yani、Zani是原始坐标；α、β、γ分别是原坐标轴x、y和z旋转至现坐标轴的所转动的角度。

3.根据权利要求1所述的考虑孔隙率影响的黄土边坡的多场耦合分析方法，其特征在于：所述步骤S2建立与孔隙率相关的本构模型，包括如下：

1)固相；

针对非均质黄土边坡进行力学分析，所考虑的总应变率分为两部分：弹性(dεe/dt)和塑性(dεvp/dt)；

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

对于塑性部分，Fm是屈服面，G是塑性势,Γ是流动性，m是应力功率；

弹性本构是采用具有杨氏模量和泊松比的标准线弹性模型，且弹性模量E随孔隙率变化：式中：E0为参考弹性模量，φ为孔隙率，φ0E为参考孔隙率；

塑性本构采用的是改进的Mohr-Coulomb模型：式中：是摩擦角，ψ是剪胀角，ω是塑性势参数，p′、q、θ分别是有效平均应力、偏应力和罗德角；

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

c＝(ac+bcs)g(φ)；

g(φ)＝(f+|f|)/2；

f(φ)＝1-(φ/φ0c)n；

式中：c为内聚力，ac、bc为内聚力参数，φ0c为参考孔隙率；

2)液相与气相；

用Van-Genuchten模型作为持水曲线来描述黄土的吸力变化，其具体表达式为：式中：Se为实际含水率；Srl为残余含水率，即土体中含水率的最小值；Sls为最大含水率，等于土体饱和时的含水率；S1为当前含水率；Pg为孔隙气压，P1为孔隙水压，二者的差值就是基质吸力；λ为持水曲线的形状控制参数，只与对应的土体相关；P则是土体的进气值，表达式如下：P0(φ)＝P0Exp{a(φ-φ0φ0S)}；

λ(φ)＝λExp{b(φ-φ0S)}；

式中：P0为T0温度下对应的进气值；σ为表面张力；φ为孔隙率，φ0S为参考孔隙率；a、b为相关参数；

达西定律被用来表征液相与气相，具体表达式为：其中：k是固有渗透率；krα是相对渗透率；μα是黏度；是两点之间的孔隙水压梯度；ρα为液体的密度；g是重力加速度矢量；

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

k＝k0Exp{bk(φ-φ0k)}；

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

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

krg＝1-krl；

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

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