简介面向地理信息系统与遥感领域学习者的六种坡度计算方法程序包为需要分析数字高程模型地形数据的开发者提供可直接运行的算法实现与配套实验资料。包内覆盖简单差分、二阶差分、三阶反距离平方权差分、三阶反距离权差分、三阶不带权差分及边框差分六类方案并给出相应的源代码与说明文档。压缩包体积约1.53MB以RAR文件形式发布适合配合实验指导书逐步理解每个算法的差分原理和适用场景。已有2038人在CSDN下载学习尤其适合高校地理信息系统相关课程实践、地形分析入门以及要处理栅格数据边界的开发人员。通过阅读实验五文档和栅格数据示例读者能快速掌握不同坡度算法的计算差异并根据实际数据特点选择或改进方法。1. 六种坡度计算方法同一段地形为什么能算出六个数同一段坡有人说它是 5.7°有人说它是 10%还有人写成 1:10——三个数描述的是同一个物理对象。让坡度计算程序变复杂的从来不是除法而是“高差从哪取”离散两点、规则 DEM、带噪点云各自的差分路径和噪声模型完全不同工程上于是沉淀出六种常用算法。下面按数据形态把六种方法分成三组比例/百分比/角度三种工程定义Horn 与 Zevenbergen-Thorne 两种栅格差分以及最小二乘平面拟合。每组都给出可直接跑的 Python 代码、参数边界和校验套路适合写地形分析工具、处理高程数据或做道路纵断面设计的人参考。2. 比例、百分比与角度坡度计算代码三种定义与一个函数跑通2.1 三种坡度定义的共同分母水平距离不是斜距三种基础定义共享同一个分子——高差 rise分母也都取水平投影距离 run。用斜距去算得到的是 sin(θ) 而不是 tan(θ)30° 坡上两者相差约 15%落到土方量或设计速度上就是系统性偏差。所以坡度计算程序的第一步都是先确认输入的高程差和水平距离处于同一投影、同一单位工程上推荐统一用米。坡度定义表达式6% 对应的值比例坡度1 : (run / rise)1 : 16.67百分比坡度(rise / run) × 100%6%角度坡度atan(rise / run)≈ 3.43°接口设计上建议三种定义同时返回因为下游消费方完全不一样道路与铁路纵断面用百分比或千分率‰边坡稳定分析用角度公式里直接出现 tan φ竣工图则习惯写比例。只返回一种调用方就会自己写换算换算代码一多单位错误就进来了。2.2 一个函数同时输出三种坡度的 python 代码import math def slope_defs(rise, run, wantall): 坡度三种定义一次算齐。 rise/run 必须是水平投影下的同单位数值推荐用米。 if run 0: raise ValueError(run 必须大于 0) percent rise / run * 100.0 degree math.degrees(math.atan2(rise, run)) ratio None if rise 0 else run / rise out {percent: percent, degree: degree, ratio: ratio} return out if want all else out[want] print(slope_defs(6.0, 100.0)) # {percent: 6.0, degree: 3.433629385640829, ratio: 16.666666666666668}实现说明atan2(rise, run)比atan(rise / run)更稳rise 为 0 时直接给出 0° 而不是除零rise 恰好为 0 时 ratio 退化为None工程上表达为“平坡”。返回 dict 而不是元组是让调用方按名字取字段避免两个浮点数的顺序在传参时被记混。2.3 单位陷阱与近似公式的适用范围最常见的坑是经纬度 DEM 的 x/y 单位是度、z 单位是米直接进这个函数算出来的“角度”没有任何物理意义。遇到这种数据先投影到 UTM 等米制坐标系或者在经度方向乘 cos(纬度) 修正水平距离再调用本函数。提示经纬度 DEM 必须投影成米制坐标系后再算坡度否则 x/y 与 z 单位不一致输出角度没有物理意义。另一个坑是百分比与角度的近似换算。小角度下 percent ≈ degree × 100 / 57.3很多人把它当精确公式用实际上 10% 坡度时误差约 0.5°到 30% 时误差已经接近 2°超过 30% 的坡千万别用线性近似老老实实走atan2。3. DEM 栅格坡度算法Horn 三阶差分与四点二阶差分的 Python 代码3.1 3×3 窗口与两条差分路径栅格 DEM 上一个像元的坡度不属于它自己而属于它所在的 3×3 邻域。窗口小于 3×3两个点只能定一条线定不了一个面窗口大于 3×3真实地形起伏会被平滑掉。以中心像元 z5 为原点邻域编号如下z1z2z3z4z5z6z7z8z9Horn 法与 Zevenbergen-Thorne 法下文简称 ZT 法最终都走同一条公式slope arctan(sqrt((dz/dx)² (dz/dy)²))区别只在 dz/dx、dz/dy 的估计方式。算法流程图里这段通常长这样输入 DEM → 差分 → arctan → 坡度图真正需要拍板的只有差分那一步的权重取法。3.2 Horn 法ArcGIS 默认的三阶加权差分Horn 法的核心是 [1,2,1] 加权后再差分dz/dx (z3 2·z6 z9 − z1 − 2·z4 − z7) / (8·d) dz/dy (z7 2·z8 z9 − z1 − 2·z2 − z3) / (8·d)这里的 d 是像元边长。1/(8d) 的来源每行先做 [1,2,1] 加权和权重合计 4再取两侧行的差对应的距离是 2d合起来除以 8d。y 方向差分的符号只影响坡向不影响坡度幅值要对齐 ArcGIS 的坡向输出把 dzdy 取反即可。import numpy as np from scipy.ndimage import convolve def horn_slope(dem, cellsize, degreesTrue, modereflect): Horn 三阶加权差分坡度ArcGIS 默认算法。 dem: 二维 ndarray 高程; cellsize: 像元边长, 与高程同单位。 dem np.asarray(dem, dtypefloat) kx np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]]) / (8.0 * cellsize) ky np.array([[-1, -2, -1], [ 0, 0, 0], [ 1, 2, 1]]) / (8.0 * cellsize) dzdx convolve(dem, kx, modemode) dzdy convolve(dem, ky, modemode) slope np.arctan(np.hypot(dzdx, dzdy)) return np.degrees(slope) if degrees else slope参数说明modereflect用镜像外推补齐边界让输出与输入形状一致改成modenearest会让边界坡度整体偏平缓改成modeconstant则会在边界制造假坡。如果后续分析只认内部像元建议把边界像元标记为无效而不是换 mode 硬凑。np.hypot先放大再开方能避免高程值很大时先平方后溢出整型 DEM 会被np.asarray(dtypefloat)自动提升但调用方最好自己转换避免在传参链路上丢失精度。3.3 ZT 法只用四个直接邻居的快速实现ZT 法假设局部地表可以用低阶二次多项式逼近在这个假设下中心像元的一阶偏导数恰好退化成最简单的四点差分dz/dx (z6 − z4) / (2d) dz/dy (z8 − z2) / (2d)def zt_slope(dem, cellsize, degreesTrue): Zevenbergen-Thorne 四点二阶差分坡度。 只使用东西南北四个直接邻居输出比输入小一圈的内域。 dem np.asarray(dem, dtypefloat) dzdx (dem[1:-1, 2:] - dem[1:-1, :-2]) / (2.0 * cellsize) dzdy (dem[2:, 1:-1] - dem[:-2, 1:-1]) / (2.0 * cellsize) slope np.arctan(np.hypot(dzdx, dzdy)) out np.full_like(dem, np.nan, dtypefloat) out[1:-1, 1:-1] np.degrees(slope) if degrees else slope return out切片写法的要点dem[1:-1, 2:]取的是右列、dem[1:-1, :-2]取的是左列顺序反了坡向会跟着翻转南北差分同理。输出用 NaN 填充边界这正是 ZT 与 Horn 的一个关键差异——Horn 给出边界外推值ZT 直接声明“边界没有坡度”。两种栅格算法的行为差异如下表。在理想平面上两者内部像元完全一致只要高程带一点噪声LiDAR 或摄影测量的典型情况ZT 的单像元噪声会直接进入结果Horn 因为做了 8 格加权等效于一个轻度低通滤波。实测中 30° 以上的陡坡区两者结果可以差 2° 到 4°。比较项Horn 法ZT 法参与差分像元9 格加权4 格抗噪性较好差噪声敏感边界输出外推值内缩一圈NaN计算方式2 次卷积4 次切片减法4. 点云坡度计算代码最小二乘平面拟合与邻域半径选择4.1 为什么点云和不规则邻域绕不开平面拟合激光雷达、摄影测量和野外实测得到的是离散 (x, y, z) 点没有规则栅格。先把点插值成 DEM 再算坡度等于把插值误差带进结果——点密度不均匀的区域插值本身就能给地形“造”出 2° 的坡度。更稳的常见做法是以目标点为中心取一个邻域在邻域内用最小二乘拟合平面 z a·x b·y c平面法向量的仰角就是局部坡度。道路纵断面也常用这个思路断面测量点高差跳动剧烈时最小二乘比任意两点法稳定得多。4.2 最小二乘坡度计算代码从法方程到坡向对包含 n 个点的邻域列出线性方程组 A·[a, b, c]ᵀ Z用np.linalg.lstsq求解。坡度只与 a、b 有关slope arctan(sqrt(a² b²))。import numpy as np def lsq_slope(xyz, wantboth): 最小二乘平面拟合坡度。 xyz: (N, 3) 数组列依次为 x, y, z建议先转成米制。 want: slope 返回坡度角(度); grad 返回平面系数(a,b,c); both 返回 (坡度, 坡向)。 pts np.asarray(xyz, dtypefloat) if pts.ndim ! 2 or pts.shape[1] ! 3 or pts.shape[0] 3: raise ValueError(需要 (N,3) 点集且至少 3 个点) A np.column_stack([pts[:, 0], pts[:, 1], np.ones(pts.shape[0])]) (a, b, c), *_ np.linalg.lstsq(A, pts[:, 2], rcondNone) if want grad: return a, b, c slope np.degrees(np.arctan(np.hypot(a, b))) if want slope: return slope # 坡向最大下坡方向的方位角正北顺时针 aspect (90.0 - np.degrees(np.arctan2(-b, -a))) % 360.0 return slope, aspect参数与实现说明rcondNone是 NumPy 处理病态矩阵的标准写法3 个点是理论下限邻域点数低于 10 时坡度方差很大工程上建议至少 5 到 10 个点再算。坡向换算结果一律% 360归一到 0~360°否则负角度会让下游色带映射出现断层。如果下游只需要“比率”而不是角度直接返回np.hypot(a, b)即可。4.3 邻域半径与距离加权平面拟合的隐藏参数平面拟合真正的调参项不是求解器而是邻域半径。半径太小噪声占主导半径太大弯曲地形被拟平成平均坡。对点云取平均点距的 2~3 倍是常见起点对道路纵断面取竖曲线长度的 1/4 到 1/2拟合出来的才是“设计纵坡”而不是局部毛刺。另一种增强是把距离纳入权重远点对平面的贡献随距离衰减def weighted_lsq(xyz, sigma): 以 xyz[0] 为锚点的高斯加权平面拟合 sigma 控制衰减半径单位与 x/y 一致。 pts np.asarray(xyz, dtypefloat) dx pts[:, 0] - pts[0, 0] dy pts[:, 1] - pts[0, 1] w np.exp(-(dx**2 dy**2) / (2.0 * sigma**2)) A np.column_stack([pts[:, 0], pts[:, 1], np.ones(pts.shape[0])]) W np.diag(w) coef np.linalg.solve(A.T W A, A.T W pts[:, 2]) return coef加权前后对照σ 取平均点距量级时结果接近普通最小二乘σ 缩小到 0.5 倍点距以下远点几乎不参与拟合结果会趋向 ZT 的行为——坡度由最近邻决定噪声也随之被放大。权重矩阵用np.diag(w)构造数据量大时改成广播乘法避免显式建矩阵。场景邻域半径建议说明机载 LiDAR 点云平均点距 × 2~3兼顾信噪比与地形保真地面实测断面竖曲线长度 1/4~1/2拟合的是“设计纵坡”混合密度点云半径 高斯加权用 σ 控制衰减避免远点主导5. 坡度计算程序的选型边界与理想斜面校验技巧5.1 六种方法的选型速查数据形态首选方法备选关键参数离散两点角度/百分比比例水平距必须投影修正规则 DEMHorn 法ZT 法cellsize 与高程同单位点云/断面最小二乘加权最小二乘邻域半径、σ常见的工作流是组合使用先用 Horn 出全图坡度再在局部特征点陡坎、坡脚用最小二乘复核两者差异超过 2° 时先怀疑数据噪声而不是算法。5.2 用解析可导的曲面校验坡度计算代码校验坡度程序最稳妥的基线是用一个数学上能求出精确坡度的曲面。平面太平凡用高斯山体f(x, y) A·exp(−(x² y²)/σ²)它的梯度有闭式解任何差分权重的改动都会被精确暴露出来。x np.linspace(-10, 10, 201) X, Y np.meshgrid(x, x) A, sigma 12.0, 7.0 dem A * np.exp(-(X**2 Y**2) / sigma**2) fx -2.0 * X / sigma**2 * dem fy -2.0 * Y / sigma**2 * dem slope_true np.degrees(np.arctan(np.hypot(fx, fy))) s_horn horn_slope(dem, 1.0) s_zt zt_slope(dem, 1.0) m np.isfinite(s_zt) # ZT 边界是 NaN顺势当成内域掩膜 print(np.abs(s_horn[m] - slope_true[m]).mean()) print(np.abs(s_zt[m] - slope_true[m]).mean())在 201×201 网格、1 米像元的设置下Horn 与 ZT 在内部像元的平均误差应该在 0.1° 量级如果均值误差超过 0.5°先查 cellsize 单位再查卷积核符号最后查坐标系是否是米制。这个验证可以写成 demo 程序放进 CI任何一次差分权重改动都会在回归测试里立刻报红。实际项目里我一般还会加一道单位用例把同一组经纬度数据在投影前后各算一遍确认坡度变化落在分辨率对应的量级内——单位修正永远比算法选型更先被测试。本文还有配套的精品资源点击获取