简介：这份资源面向地质工程、土木工程领域的研究人员与工程师，聚焦利用Python实现边坡稳定性数值分析，帮助读者评估山坡、路堑、堤坝等土体结构在自然或人为因素下的稳定性。压缩包共19个文件，约38KB，以15个csv地质参数数据文件为主，另含2个Python脚本、1个xlsx材料计算表及1个for源码文件，覆盖材料分层、水合物、分解区等建模数据与计算逻辑。资源围绕圆弧滑动面法、简化Bishop法、Janbu法、Morgenstern-Price法等常见评价方法，提供数据输入、地形建模、滑移面搜索与安全系数计算所需的脚本与参数支撑，可帮助读者理解土壤重度、内摩擦角、黏聚力等参数如何参与稳定性判定。目前已有600余人学习下载，适合希望将Python科学计算能力应用于边坡工程实践、提升数值分析效率的初中级学习者参考。
1. 边坡稳定性计算文件到底长什么样：从一份 Python 脚本说起
很多人第一次接触边坡稳定性计算，是从一份 Excel 表格或者某个商业软件导出的报告开始的。但当你需要批量处理几十个剖面、反复调整参数、或者把计算逻辑嵌进自己的监测系统时，图形界面就成了累赘。这时候用 Python 写边坡稳定性计算文件，本质上就是把你手动点鼠标的过程翻译成可复现的代码。它解决的核心问题是：同一套地质参数和几何模型，能不能在不同时间、不同人手里跑出完全一致的安全系数。适合谁？适合手里已经有勘察剖面数据、想摆脱重复劳动、又愿意花一个下午把公式敲进编辑器的岩土工程师和地质专业学生。这一章不急着写代码，先把一份典型的计算文件拆开看，它到底由哪几块组成，每块对应现实中的什么。
一份能跑的边坡稳定性计算文件，通常包含四个部分：几何模型、岩土参数、计算方法、输出结果。几何模型就是坡面线、地层线、地下水位线，用坐标点列表表示；岩土参数是容重、黏聚力、内摩擦角，可能还有孔隙水压力系数；计算方法是瑞典条分法、毕肖普法还是简布法，决定了公式里怎么迭代；输出结果是每个滑弧对应的安全系数和圆心坐标。很多人卡在第一步，是因为把 CAD 里的多段线直接复制过来，结果坐标顺序反了，或者单位混了米和毫米。我一般会先用一个极简的均质边坡例子把整条链路跑通，确认安全系数和手算或商业软件对得上，再往里面加地层和水位。这样做的另一个好处是，后面排查错误时你知道问题出在新增的那部分，而不是一开始就面对一团乱麻。
2. 用 Python 搭一个边坡稳定性计算文件的最小骨架
2.1 几何数据怎么组织才不翻车
几何数据最常见的翻车现场是坡面线坐标点顺序混乱。假设你从勘察报告里抄了一组坡面点，从左到右应该是 x 递增，但有人从坡顶往坡脚抄，结果 x 递减，后面算滑弧与坡面交点时就会得到负的条块宽度。我一般用列表套元组的方式存，每个元组是 (x, y)，并且强制自己按 x 从小到大排。下面这段代码定义了一个简单的均质边坡，坡高 10 米，坡比 1:1.5，坡顶和坡脚各延伸一段水平段，避免滑弧从边界外穿出去。
# 边坡几何：从左到右，x 递增 slope_surface = [ (0.0, 10.0), # 坡顶左侧延伸段 (5.0, 10.0), # 坡顶 (20.0, 0.0), # 坡脚，坡比 1:1.5，水平 15m 下降 10m (30.0, 0.0), # 坡脚右侧延伸段 ] # 地层线：这里简化为均质，只有一条地面线 ground_lines = [slope_surface] # 地下水位线：暂时不考虑，设为 None water_table = None逻辑说明：slope_surface 是坡面轮廓，ground_lines 目前和坡面重合，表示均质边坡。water_table 为 None 表示干坡。参数怎么改？坡高改 y 差值，坡比改水平段长度。注意 x 必须严格递增，否则后面条分法分条时会出现负宽度。如果你有多个地层，ground_lines 里放多条线，每条线也是从左到右的坐标点，并且要保证相邻地层线在 x 方向上有重叠区间，不然插值会报错。
2.2 岩土参数和单位统一
岩土参数翻车最多的是单位。勘察报告里黏聚力 c 常用 kPa，内摩擦角 φ 用度，容重 γ 用 kN/m³。但 Python 的 math.tan 要弧度，所以每次算 tan(φ) 之前必须 math.radians(φ)。我见过有人直接把 30 度带进去，结果 tan(30) 算出来是 -6.4，安全系数直接负到离谱。下面定义参数时，我习惯把单位写在注释里，并且把角度转弧度单独做一步。
import math # 岩土参数：均质土坡 gamma = 18.0 # 容重，kN/m³ c = 20.0 # 黏聚力，kPa phi_deg = 25.0 # 内摩擦角，度 phi_rad = math.radians(phi_deg) # 转弧度 # 孔隙水压力系数：干坡取 0 ru = 0.0逻辑说明：gamma 用于计算条块重量，c 和 phi_rad 用于抗剪强度，ru 用于有效应力法。参数怎么改？如果有多层土，每个条块要根据其底部中点落在哪个地层来取对应的 c 和 phi。常见做法是写一个函数，输入条块底部中点的 y 坐标，返回该点所在地层的参数。这一步别偷懒，否则均质坡算惯了，一遇到多层土就全错。
2.3 瑞典条分法的最小实现
瑞典条分法虽然粗，但它是理解边坡稳定性计算文件最好的起点。核心思想是把滑弧以上的土体切成若干竖直条块，对每个条块算重力和抗滑力，然后对圆心取矩。下面这段代码实现了一个固定圆心的安全系数计算，滑弧用圆心和半径定义，条块宽度固定。
def swedish_fs(cx, cy, radius, surface, gamma, c, phi_rad, ru, n_slices=50): """ 瑞典条分法计算安全系数 cx, cy: 滑弧圆心坐标 radius: 滑弧半径 surface: 坡面线坐标列表 n_slices: 条块数量 """ # 找滑弧与坡面的交点，确定积分范围 x_min = min(p[0] for p in surface) x_max = max(p[0] for p in surface) # 简单处理：在坡面 x 范围内均匀分条 dx = (x_max - x_min) / n_slices numerator = 0.0 # 抗滑力矩 denominator = 0.0 # 滑动力矩 for i in range(n_slices): x_left = x_min + i * dx x_mid = x_left + dx / 2.0 # 坡面在 x_mid 处的 y 值，线性插值 y_surface = interp_surface(surface, x_mid) # 滑弧在 x_mid 处的 y 值，取下半圆 y_arc = cy - math.sqrt(max(radius**2 - (x_mid - cx)**2, 0.0)) # 条块高度 h = y_surface - y_arc if h <= 0: continue # 条块重量 W = gamma * h * dx # 条块底部倾角 sin_alpha = (x_mid - cx) / radius cos_alpha = math.sqrt(1.0 - sin_alpha**2) # 抗滑力：c*l + N*tan(phi)，N = W*cos_alpha # 简化：不考虑条块间作用力 l = dx / cos_alpha # 条块底部弧长 N = W * cos_alpha T = W * sin_alpha numerator += c * l + (N - ru * W * cos_alpha) * math.tan(phi_rad) denominator += T if denominator == 0: return float('inf') return numerator / denominator def interp_surface(surface, x): """在坡面线上线性插值求 y""" for i in range(len(surface) - 1): x1, y1 = surface[i] x2, y2 = surface[i + 1] if x1 <= x <= x2: t = (x - x1) / (x2 - x1) return y1 + t * (y2 - y1) return surface[-1][1]逻辑说明：swedish_fs 接收圆心、半径、坡面线、参数和条块数。先确定坡面 x 范围，均匀分条。对每个条块，用 interp_surface 求坡面 y，用圆方程求滑弧 y，得到条块高度。然后算重量 W、底部倾角对应的 sin 和 cos、底部弧长 l。分子累加抗滑力，分母累加下滑力。最后返回比值。参数怎么改？n_slices 越大结果越稳，但超过 200 后变化很小，我一般用 50 到 100。ru 如果非零，表示孔隙水压力系数，会降低有效正应力。注意这个实现没有考虑条块间作用力，所以算出来的安全系数比毕肖普法低 5% 到 10%，这是瑞典条分法本身的保守性，不是代码 bug。
3. 从固定圆心到自动搜索最危险滑弧
3.1 圆心搜索网格怎么设才不漏掉最小值
固定圆心只能算一个滑弧，实际工程要的是最危险滑弧，也就是安全系数最小的那个。常见做法是在坡脚上方和坡顶后方画一个网格，对每个网格点作为圆心，再对每个圆心试算一组半径，取最小安全系数。网格范围怎么定？我一般用坡脚和坡顶的坐标来框：圆心 x 从坡脚左侧 0.5 倍坡高到坡顶右侧 1.5 倍坡高，圆心 y 从坡脚以下 0.5 倍坡高到坡顶以上 2 倍坡高。半径从圆心到坡面的最短距离开始，逐步增加到能覆盖整个坡体。下面是一个搜索框架。
def search_critical_circle(surface, gamma, c, phi_rad, ru): """网格搜索最危险滑弧""" x_min = min(p[0] for p in surface) x_max = max(p[0] for p in surface) y_min = min(p[1] for p in surface) y_max = max(p[1] for p in surface) H = y_max - y_min # 坡高 best_fs = float('inf') best_circle = None # 圆心网格 cx_list = [x_min - 0.5*H + i*0.2*H for i in range(11)] cy_list = [y_min - 0.5*H + j*0.2*H for j in range(16)] for cx in cx_list: for cy in cy_list: # 半径范围：从圆心到坡面的最短距离到 2 倍坡高 r_min = min(math.hypot(cx - p[0], cy - p[1]) for p in surface) r_max = 2.0 * H for k in range(20): r = r_min + (r_max - r_min) * k / 19.0 fs = swedish_fs(cx, cy, r, surface, gamma, c, phi_rad, ru) if fs < best_fs: best_fs = fs best_circle = (cx, cy, r) return best_fs, best_circle逻辑说明：cx_list 和 cy_list 是圆心网格，步长 0.2H，范围覆盖坡脚外侧和坡顶上方。对每个圆心，半径从到坡面最短距离开始，到 2 倍坡高结束，分 20 步。每次调用 swedish_fs 算安全系数，记录最小值。参数怎么改？步长 0.2H 是精度和速度的折中，如果坡形复杂可以降到 0.1H，但计算量翻四倍。半径上限 2H 对大多数边坡够用，如果滑弧很深可以加到 3H。注意这个搜索是暴力网格，没有优化算法，跑一个边坡大概几秒到几十秒，取决于条块数和网格密度。
3.2 毕肖普法迭代为什么容易不收敛
毕肖普法比瑞典条分法多了一个条块间水平作用力，安全系数出现在等式两边，需要迭代。很多人写迭代时直接 while 循环，结果要么不收敛，要么收敛到错误的值。血泪经验是：迭代初值取瑞典条分法的结果，收敛容差设 0.001，最大迭代次数设 100，并且每次迭代后检查安全系数是否在合理范围（0.5 到 5.0），超出就说明参数或几何有问题。下面是一个简化毕肖普法的迭代片段。
def bishop_fs(cx, cy, radius, surface, gamma, c, phi_rad, ru, n_slices=50, tol=0.001, max_iter=100): """简化毕肖普法，迭代求安全系数""" # 先算瑞典法作为初值 fs = swedish_fs(cx, cy, radius, surface, gamma, c, phi_rad, ru, n_slices) for iteration in range(max_iter): numerator = 0.0 denominator = 0.0 dx = (max(p[0] for p in surface) - min(p[0] for p in surface)) / n_slices for i in range(n_slices): x_mid = min(p[0] for p in surface) + (i + 0.5) * dx y_surface = interp_surface(surface, x_mid) y_arc = cy - math.sqrt(max(radius**2 - (x_mid - cx)**2, 0.0)) h = y_surface - y_arc if h <= 0: continue W = gamma * h * dx sin_alpha = (x_mid - cx) / radius cos_alpha = math.sqrt(1.0 - sin_alpha**2) # 毕肖普法分母项 m_alpha = cos_alpha + sin_alpha * math.tan(phi_rad) / fs if m_alpha <= 0: return float('inf') numerator += (c * dx / cos_alpha + (W - ru * W * cos_alpha) * math.tan(phi_rad) / cos_alpha) / m_alpha denominator += W * sin_alpha fs_new = numerator / denominator if abs(fs_new - fs) < tol: return fs_new fs = fs_new return fs # 达到最大迭代次数，返回当前值逻辑说明：bishop_fs 先调用 swedish_fs 得到初值，然后进入迭代。每次迭代重新计算分子和分母，分子里包含 m_alpha 项，它依赖当前 fs。如果 m_alpha 小于等于零，说明迭代发散，直接返回无穷大。收敛后返回 fs_new。参数怎么改？tol 和 max_iter 一般不用动，如果遇到不收敛，先检查几何数据有没有负高度条块，再检查 phi_rad 是不是忘了转弧度。注意这个实现是简化毕肖普，忽略了条块间竖向作用力，对大多数均质土坡精度足够，但遇到软弱夹层要改用完整毕肖普或摩根斯坦-普赖斯法。
4. 边坡稳定性计算文件避坑与排查
4.1 安全系数算出来是负数或者无穷大
现象：运行脚本后，安全系数显示 -3.2 或者 inf。原因：负数通常是因为分母累加的下滑力为负，也就是 sin_alpha 在大部分条块上为负，说明圆心位置在坡面以下或者半径太小，滑弧方向反了。无穷大是因为分母为零，所有条块的 sin_alpha 都接近零，圆心在坡面正上方且半径极大。解决：检查圆心 y 坐标是否低于坡脚，半径是否小于圆心到坡面的最短距离。我一般会在搜索前加一个断言，确保圆心在坡面以上，半径大于最短距离。
4.2 多层土时安全系数突然跳变
现象：均质坡算出来 1.5，加了第二层土后变成 0.8 或者 3.0。原因：条块底部中点落在不同地层时，c 和 phi 取值函数写错了，比如用条块顶部中点判断地层，或者地层线插值范围没覆盖整个滑弧。解决：写一个 get_soil_params(x, y) 函数，输入条块底部中点的 x 和 y，遍历所有地层线，找到该点所在的地层。注意地层线要按从上到下排序，并且相邻地层线在 x 方向上有重叠。测试时先用手算一个两层土的简单例子验证。
4.3 搜索最危险滑弧时结果每次不一样
现象：同样的参数，两次运行得到的最危险圆心差了几米，安全系数差 0.05。原因：网格搜索的步长和半径分步数固定，但浮点数累加顺序不同会导致微小差异，如果安全系数曲面很平，最小值点就会漂移。解决：固定随机种子（如果用了随机搜索），或者把网格步长和半径分步数写死在配置文件里。更稳的做法是用优化算法比如 scipy.optimize.minimize 替代暴力网格，但要注意优化算法可能陷入局部极小，需要多起点。
4.4 条块数量取多少才够
现象：n_slices 从 50 改成 100，安全系数变了 0.02，再改成 200 又变 0.005。原因：条分法本身是数值积分，条块越多越接近连续解，但超过一定数量后变化小于工程精度。解决：我一般用 50 到 100，并且做一个收敛性测试：把 n_slices 从 20 逐步加到 200，画安全系数随条块数的变化曲线，当相邻两次变化小于 0.001 时，取当前值。这个测试跑一次就行，后面同类型边坡直接用这个条块数。
4.5 地下水位线怎么加才不出错
现象：加了水位线后安全系数反而变大。原因：水位线以下应该用浮容重和孔隙水压力，但代码里只改了容重没改 ru，或者 ru 符号搞反了。解决：水位线以下条块，容重用 gamma_sat - gamma_w，同时 ru 取 0.3 到 0.5 的经验值，或者直接用孔隙水压力 u = gamma_w * h_w，在抗滑力里减去 u * l。注意水位线要和坡面线、地层线一起插值，确保每个条块底部中点能查到对应的水头高度。
5. 把计算文件变成能复用的工具：参数化与批量跑
5.1 用配置文件管理多组参数
当你需要算 20 个剖面时，把参数写死在代码里就是灾难。我一般用一个 JSON 文件存每个剖面的几何和岩土参数，Python 读进来循环调用。下面是一个配置示例和读取代码。
{ "slope_1": { "surface": [[0, 10], [5, 10], [20, 0], [30, 0]], "gamma": 18.0, "c": 20.0, "phi_deg": 25.0, "ru": 0.0 }, "slope_2": { "surface": [[0, 12], [6, 12], [24, 0], [36, 0]], "gamma": 19.0, "c": 15.0, "phi_deg": 22.0, "ru": 0.2 } }import json with open('slopes.json', 'r', encoding='utf-8') as f: config = json.load(f) results = {} for name, params in config.items(): surface = [tuple(p) for p in params['surface']] phi_rad = math.radians(params['phi_deg']) fs, circle = search_critical_circle( surface, params['gamma'], params['c'], phi_rad, params['ru'] ) results[name] = {'fs': fs, 'circle': circle} print(f"{name}: FS={fs:.3f}, circle={circle}")逻辑说明：JSON 里每个键是一个剖面名，值包含 surface、gamma、c、phi_deg、ru。读取后转成元组列表，调用 search_critical_circle。结果存字典并打印。参数怎么改？新增剖面就在 JSON 里加一个键。注意 surface 的坐标点顺序必须从左到右，x 递增。如果剖面很多，可以把结果写回 CSV，方便后续出报告。
5.2 用多进程加速批量计算
20 个剖面串行跑可能要几分钟，用 multiprocessing 可以压到几十秒。下面是一个简单的多进程封装。
from multiprocessing import Pool def compute_one(args): name, params = args surface = [tuple(p) for p in params['surface']] phi_rad = math.radians(params['phi_deg']) fs, circle = search_critical_circle( surface, params['gamma'], params['c'], phi_rad, params['ru'] ) return name, fs, circle if __name__ == '__main__': with open('slopes.json', 'r', encoding='utf-8') as f: config = json.load(f) with Pool(processes=4) as pool: results = pool.map(compute_one, config.items()) for name, fs, circle in results: print(f"{name}: FS={fs:.3f}, circle={circle}")逻辑说明：compute_one 接收 (name, params) 元组，返回结果。Pool 创建 4 个进程，map 并行计算。注意 Windows 下必须把入口放在 ifname== 'main' 里，否则会递归创建进程。参数怎么改？processes 设成 CPU 核心数，一般 4 到 8。如果剖面数量少于进程数，有些进程会空闲，但没关系。
5.3 结果验证：和商业软件对一对
写完自己的计算文件，一定要拿一个简单例子和商业软件或者手算结果对。我一般用均质土坡，坡高 10 米，坡比 1:1.5，c=20 kPa，phi=25 度，gamma=18 kN/m³，干坡。瑞典条分法手算安全系数大概 1.35 到 1.45，毕肖普法大概 1.45 到 1.55。如果你的代码算出来在这个范围外，先检查单位，再检查条块底部倾角的正负号，最后检查滑弧与坡面交点是不是取错了。验证通过后，再逐步加水位、加地层、加荷载。这个顺序能帮你把错误隔离在新增的那一步。
5.4 一个我常用的收敛性检查技巧
搜索最危险滑弧时，网格步长和半径分步数决定了结果精度。我习惯在正式跑之前，先用一个粗网格（步长 0.5H，半径 10 步）快速扫一遍，找到安全系数最小的区域，然后在这个区域周围用细网格（步长 0.1H，半径 30 步）再扫一次。这样比全局细网格快好几倍，而且不容易漏掉真正的最小值。具体做法是：第一遍记录前 5 个最小安全系数对应的圆心和半径，第二遍只在这些圆心周围 ±0.3H 范围内细化。这个技巧在坡形复杂时特别有用，因为全局细网格可能跑一晚上，而两阶段搜索几分钟就出结果。希望帮到你。
本文还有配套的精品资源，点击获取
