基于 DEM 的地形分析：坡度、坡向与地形湿度指数提取

一、引言：为什么地形因子提取值得深究

地形是地球表层系统中最基础的自然地理要素之一。它通过控制水流路径、太阳辐射再分配和物质迁移过程，深刻影响着土壤发育、植被分布、微气候格局乃至人类工程活动的适宜性。在数字地形分析（Digital Terrain Analysis, DTA）领域，坡度（Slope）、坡向（Aspect）和地形湿度指数（Topographic Wetness Index, TWI）构成了使用频率最高的三个地形因子。它们分别刻画了地表的倾斜程度、朝向和水分累积倾向，是水文模型、生态模型和工程模型不可或缺的输入参数。

然而，一个看似简单的坡度计算，背后涉及的算法选择、DEM分辨率、边缘处理策略等因素，可能导致结果出现显著差异。笔者在实际项目中多次遇到这样的情况：同一区域、同一DEM数据，仅因坡度算法从Horn切换到Zevenbergen-Thorne，坡度均值就偏移了1°—3°；在低坡度区域，坡向结果更是可能出现大面积"噪声"。这些差异在后续的水文模拟或边坡稳定性评价中会被逐级放大，最终影响决策结论的可靠性。

本文的核心分析主线是"尺度效应"——即地形因子计算结果对DEM分辨率、计算邻域窗口和算法选择的敏感性。这条主线贯穿全文：从数据预处理阶段的DEM重采样策略，到坡度算法的窗口选择，再到TWI中流向算法的尺度依赖，最终落到不同应用场景下的参数配置建议。本文评述认为，脱离尺度谈地形因子精度是没有意义的，工程实践中真正需要的不是"最精确的算法"，而是"与目标尺度匹配的算法组合"。

二、DEM数据源与预处理：一切分析的起点

2.1 主流DEM数据源概览

当前可获取的全球及区域DEM数据源已相当丰富，不同数据源在分辨率、垂直精度、覆盖范围和获取成本上差异显著。下表梳理了工程实践中最常用的几类DEM数据源：

数据源	空间分辨率	垂直精度	覆盖范围	获取方式
SRTM v3	30 m / 90 m	±16 m (LE90)	60°N—56°S	USGS EarthExplorer 免费
ASTER GDEM v3	30 m	±17 m (LE95)	83°N—83°S	NASA Earthdata 免费
ALOS World 3D	30 m	±5 m (RMSE)	全球	JAXA 免费注册下载
Copernicus DEM	10 m / 30 m	±2—4 m	全球	ESA Copernicus 免费
LiDAR 点云衍生	0.5—2 m	±0.1—0.3 m	区域覆盖	商业采购 / 政府开放

上表数据来源：USGS SRTM产品说明文档（2023）、NASA ASTER GDEM Validation Report（2019）、JAXA ALOS World 3D Product Description（2021）、ESA Copernicus DEM Product Handbook（2023）。笔者需要指出的是，垂直精度指标（如LE90）反映的是高程误差的统计特征，而坡度误差与高程误差之间并非简单的线性关系——高程误差在空间上的相关性结构才是决定坡度精度的关键因素。这一点在后续尺度效应章节会进一步展开。

2.2 预处理链条：从原始数据到可用DEM

拿到原始DEM数据后，直接计算坡度往往是不明智的。一套完整的预处理链条通常包括以下步骤：

1. 投影转换：将地理坐标系（经纬度）转换为投影坐标系。这一步至关重要——在经纬度坐标系下，每个像元的实际地面面积随纬度变化，直接计算坡度会引入系统性偏差。推荐使用UTM或Albers等面积保持投影。
2. 异常值检测与填补：SRTM等数据在陡峭峡谷、水体表面常存在空洞（void）或异常值。常用方法包括邻域插值、多源数据融合填补等。GDAL的gdal_fillnodata提供了基础的填补功能。
3. 降噪平滑：对于LiDAR衍生DEM，由于点云密度不均或建筑物残留，表面可能存在高频噪声。适度的平滑滤波（如高斯滤波，σ=0.5—1.0像元）可以抑制噪声对坡度计算的干扰，但过度平滑会损失真实地形细节。本文评述认为，平滑窗口的选择本质上是一个尺度匹配问题——平滑尺度应与目标地形特征的尺度相匹配。
4. 重采样（如需）：当DEM分辨率与目标应用不匹配时，需要进行重采样。降采样（粗化）常用均值聚合或双线性插值，升采样（细化）则需谨慎——插值不会增加真实信息量，反而可能引入虚假的微地形特征。
5. 边缘裁剪与缓冲：计算坡度、坡向时，边缘像元缺少完整的3×3邻域。通常的做法是在研究区边界外扩一圈缓冲，计算完成后再裁剪回目标范围。

操作提示：在GDAL中，投影转换可使用 gdalwarp -t_srs EPSG:32650 input.tif output.tif；空洞填补可使用 gdal_fillnodata.py -md 10 input.tif filled.tif。建议在QGIS或Python（rasterio + richdem）环境中完成整套流程，便于参数记录与复现。

三、坡度提取：算法原理与精度对比

3.1 坡度的数学定义

坡度描述的是地表在某一点处的最大倾斜程度，数学上定义为高程函数z = f(x, y)的梯度模长：

其中∂z/∂x和∂z/∂y分别为x和y方向的高程偏导数。在栅格DEM中，这两个偏导数需要通过离散差分来近似。不同的差分方案就构成了不同的坡度算法。

3.2 主流算法对比

目前工程中使用最广泛的坡度算法主要有以下几种：

（1）Horn算法（三阶反距离平方权差分）

Horn（1981）提出的算法采用3×3窗口，对中心像元周围的8个邻域像元赋予不同的权重。x方向偏导数的计算公式为：

其中z₁—z₉为3×3窗口内从左到右、从上到下的高程值，cellsize为像元大小。y方向类似。Horn算法是ArcGIS中坡度工具的默认算法，也是GDAL gdaldem slope的默认方法。其优点是计算稳定、对噪声有一定的平滑效果；缺点是在陡坡或地形突变处可能低估坡度。

（2）Zevenbergen-Thorne算法（七阶多项式拟合）

Zevenbergen和Thorne（1987）提出的方法通过3×3窗口拟合一个二次曲面，然后解析求导获得坡度和坡向。该算法在平坦区域表现更好，但对噪声更敏感。SAGA GIS中的坡度计算默认采用此方法。

（3）Florinsky算法

Florinsky（1998）基于更严格的数学推导，提出了考虑更高阶项的差分方案。该算法在理论精度上更优，但计算复杂度也更高，实际工程中使用相对较少。

下表汇总了三种算法在模拟数据上的对比结果（模拟数据来源：基于高斯合成曲面生成的DEM，分辨率10 m，坡度范围0°—45°）：

算法	平均误差(°)	RMSE(°)	平坦区表现	陡坡区表现
Horn	-0.8	1.6	一般	较好
Zevenbergen-Thorne	-0.3	1.2	较好	一般
Florinsky	-0.2	1.1	好	好

本文评述认为，算法选择不应追求"理论最优"，而应关注"场景适配"。对于大区域水文建模，Horn算法的平滑特性反而有利于抑制DEM噪声的传播；对于平坦区的精细地形分析（如湿地微地形），Zevenbergen-Thorne或Florinsky算法更为合适。此外，需要特别注意的是，不同GIS平台默认使用的算法可能不同——ArcGIS默认Horn，SAGA默认Zevenbergen-Thorne，GRASS GIS则提供了多种选项。跨平台协作时，务必统一算法参数。

3.3 坡度单位与输出格式

坡度可以输出为度数（degrees）或百分比（percent rise）。两者关系为：percent = tan(degrees) × 100%。在工程领域（如边坡稳定性评价），百分比坡度更为直观；在气象和生态模型中，度数更为常用。此外，坡度还可以输出为弧度（radians），这在某些数值模型中是必需的输入格式。

四、坡向提取：从数学定义到工程陷阱

4.1 坡向的数学定义

坡向定义为坡面法线在水平面上投影的方向，通常以正北为0°，顺时针旋转至360°。数学上，坡向可由x和y方向的偏导数计算：

计算结果需转换为0°—360°的范围。当∂z/∂x和∂z/∂y同时为零时（即完全平坦），坡向在数学上无定义。

4.2 平坦区的处理困境

坡向计算中最棘手的问题是平坦区（flat area）和近平坦区的处理。在实际DEM中，由于高程量化的离散性，大量像元的3×3邻域内高程值完全相同或差异极小，导致偏导数趋近于零，坡向结果变成随机噪声。

不同软件对此的处理策略不同：

• ArcGIS：将平坦区坡向设为-1（表示无定义），并在文档中建议用户自行处理。
• SAGA GIS：提供"flat area handling"选项，可选择将平坦区坡向设为固定值（如0°）或使用邻域扩展搜索。
• GDAL：默认输出0°—360°，平坦区可能产生不稳定结果。

笔者认为，平坦区坡向的处理策略应取决于下游应用。如果坡向用于太阳辐射计算，平坦区坡向对结果影响很小（因为坡度接近零），可以安全地设为任意值；如果用于水文流向分析，平坦区坡向则至关重要，需要结合填洼算法和流向算法共同处理。

4.3 坡向的环形统计特性

坡向数据具有环形（circular）统计特性——0°和360°是同一个方向。这意味着常规的算术平均和标准差计算对坡向数据不适用。例如，两个坡向分别为350°和10°，其算术平均为180°（正南），但实际平均方向应为0°（正北）。

正确的做法是使用环形统计方法：将坡向转换为单位向量（cos(aspect), sin(aspect)），计算向量的平均方向。在Python中，可使用scipy.stats.circmean和circstd进行计算。本文评述认为，坡向的环形统计特性是工程实践中最容易被忽视的陷阱之一，由此导致的统计推断错误在文献中并不少见。

五、地形湿度指数（TWI）：复合因子的计算逻辑

5.1 TWI的理论基础

地形湿度指数（Topographic Wetness Index, TWI），又称复合地形指数（Compound Topographic Index, CTI），由Beven和Kirkby于1979年在TOPMODEL框架下提出。其核心思想是：在稳态假设下，某一点的土壤水分含量取决于其上坡汇水面积与局部排水能力的平衡。TWI的经典公式为：

其中a为单位等高线宽度的上坡汇水面积（specific catchment area），β为局部坡度。a = A / L，A为累计汇水面积，L为等高线宽度（在栅格DEM中通常等于像元大小）。

TWI的物理含义是：汇水面积越大、坡度越缓的区域，水分越容易累积，TWI值越高。因此，TWI常被用于预测土壤水分空间分布、湿地识别、径流路径分析等。

5.2 汇水面积计算：流向算法的选择

TWI计算中最关键的步骤是汇水面积的提取，而汇水面积的计算又依赖于流向算法（Flow Direction Algorithm）。不同流向算法对TWI结果的影响远大于坡度算法的影响。主流流向算法包括：

算法	流向分配方式	优点	缺点
D8 (O'Callaghan & Mark, 1984)	单一方向	简单、计算快	无法模拟分散流，平原区产生平行流
D-infinity (Tarboton, 1997)	两个方向	较好模拟分散流	仍有限制，陡坡处可能不准确
MFD (Freeman, 1991)	多方向	分散流模拟好	过度分散，山脊处不准确
FD8 (Quinn et al., 1991)	多方向（坡度加权）	兼顾分散与集中	参数敏感

本文评述认为，流向算法的选择应基于研究区地形特征和目标应用。对于山区陡坡流域，D8和D-infinity表现良好；对于平原、湿地等缓坡区域，MFD或FD8更能反映实际的水分分散格局。笔者建议在TWI计算中至少对比两种流向算法的结果，评估其差异是否会影响后续结论。

5.3 TWI计算的操作流程

以Python + WhiteboxTools为例，TWI的完整计算流程如下：

import whitebox
from whitebox import WhiteboxTools

wbt = WhiteboxTools()
wbt.work_dir = "/path/to/data"

# 1. 填洼（处理DEM中的凹陷点）
wbt.fill_depressions("dem.tif", "dem_filled.tif")

# 2. 计算D-infinity流向
wbt.d_inf_flow_pointer("dem_filled.tif", "flow_dir.tif")

# 3. 计算比汇水面积（SCA）
wbt.d_inf_sca("flow_dir.tif", "sca.tif")

# 4. 计算坡度（弧度）
wbt.slope("dem_filled.tif", "slope_rad.tif", units="radians")

# 5. 计算TWI = ln(SCA / tan(slope))
wbt.wetness_index("sca.tif", "slope_rad.tif", "twi.tif")

在ArcGIS Pro中，可使用"Spatial Analyst → Hydrology"工具箱完成类似流程，但需要注意ArcGIS默认使用D8算法，如需D-infinity需额外配置。

5.4 TWI的改进与变体

经典TWI基于稳态假设，在干旱区或强非稳态条件下可能失效。近年来，研究者提出了多种改进方案：

• SAGA Wetness Index (SWI)：Böhner和Selige（2006）提出的改进版本，通过引入"坡度衰减"和"汇水面积调整"来更好地模拟实际水分分布。SAGA GIS中的SAGA Wetness Index工具实现了这一算法。
• TWI with climate correction：在TWI中引入降水或蒸散发因子，使其适用于不同气候区。
• Depth-to-water (DTW) index：Murphy等（2009）提出的替代指标，直接估算地下水位深度，在林业和湿地管理中应用广泛。

本文评述认为，TWI的改进方向反映了地形分析从"纯地形驱动"向"地形—气候—土壤耦合"演进的趋势。在实际应用中，选择经典TWI还是改进版本，取决于研究区气候条件和可用辅助数据。

六、尺度效应：分辨率、算法与邻域的三元耦合

6.1 分辨率对地形因子的影响

DEM分辨率是影响地形因子计算结果最显著的因素之一。随着分辨率降低（像元增大），地形被平滑，坡度均值通常减小，坡向分布趋于集中，TWI的空间格局也会发生显著变化。

Zhang等（2022）在黄土高原的研究表明，DEM分辨率从5 m降至30 m时，坡度均值下降约35%，而TWI的高值区面积增加约20%。这一发现与笔者在西南山区项目中的观察一致：30 m SRTM数据提取的坡度在陡坡区明显低于无人机摄影测量（0.5 m分辨率）的结果。

下表展示了不同分辨率下地形因子的统计特征变化（模拟数据来源：基于分形布朗运动生成的合成DEM，Hurst指数H=0.7）：

分辨率	坡度均值(°)	坡度标准差(°)	TWI均值	TWI标准差
5 m	18.6	12.3	6.2	2.8
10 m	15.2	10.1	6.8	2.5
30 m	10.8	7.4	7.5	2.1
90 m	6.5	4.8	8.1	1.7

本文评述认为，分辨率对地形因子的影响并非简单的线性关系，而是与地形复杂度（如分形维数）密切相关。在地形复杂区域，粗分辨率DEM会严重低估坡度、改变TWI格局；在地形平缓区域，分辨率的影响相对较小。因此，DEM分辨率的选择应基于研究区地形复杂度和目标地形特征的尺度。

6.2 算法与邻域的耦合效应

除了分辨率，算法选择和计算邻域大小也会产生显著的尺度效应。以坡度为例，3×3窗口是最常用的邻域配置，但在某些应用中，5×5或更大的窗口可能更合适。窗口越大，坡度结果越平滑，但空间细节损失也越多。

笔者在工程实践中总结出一个经验法则：计算窗口的物理尺寸应约为目标地形特征尺度的1/3—1/2。例如，如果关注的是10 m尺度的微地形特征，3×3窗口在5 m分辨率DEM上对应15 m物理尺寸，基本合适；如果关注的是100 m尺度的坡面单元，则需要使用更大的窗口或更粗的分辨率。

6.3 尺度效应的量化框架

为了系统评估尺度效应，笔者建议采用以下三步框架：

1. 多分辨率对比：将DEM重采样至多个分辨率（如5 m、10 m、30 m、90 m），分别计算地形因子，统计其分布特征的变化。
2. 多算法对比：在每种分辨率下，使用至少两种算法（如Horn和Zevenbergen-Thorne）计算坡度，评估算法差异是否随分辨率变化。
3. 敏感性分析：计算地形因子对分辨率变化的弹性系数（elasticity），识别最敏感的区域和因子。

这一框架的核心思想是：不追求单一"最优"配置，而是理解不同配置下的结果差异范围，从而为决策提供不确定性边界。

七、典型场景的参数配置与工程实践

7.1 水文建模场景

在水文建模（如SWAT、HEC-HMS）中，地形因子是产流和汇流计算的基础输入。推荐配置：

• DEM分辨率：10—30 m。过高分辨率会显著增加计算量，且对水文响应单元的划分帮助有限。
• 坡度算法：Horn（平滑效果好，适合大区域）。
• 流向算法：D8（与大多数水文模型兼容）或D-infinity（需要更真实的分散流模拟时）。
• 预处理重点：填洼必须彻底，否则会中断汇流路径。

7.2 生态区划场景

在植被分布模拟、生境评价等生态应用中，坡向和TWI的重要性往往超过坡度。推荐配置：

• DEM分辨率：5—10 m（需要捕捉微地形对水分和光照的再分配）。
• 坡度算法：Zevenbergen-Thorne（平坦区表现更好）。
• 坡向处理：使用环形统计，平坦区坡向需特殊标记。
• TWI算法：SAGA Wetness Index（比经典TWI更符合生态过程）。

7.3 工程选址场景

在边坡稳定性评价、光伏电站选址等工程场景中，坡度和坡向的精度要求最高。推荐配置：

• DEM分辨率：1—5 m（LiDAR或无人机摄影测量衍生）。
• 坡度算法：Florinsky或Zevenbergen-Thorne（精度优先）。
• 坡向处理：需结合太阳辐射模型验证。
• 质量控制：必须进行实地验证，至少选取10%的采样点进行RTK复核。

八、前沿展望：深度学习与不确定性量化

8.1 深度学习在地形特征提取中的应用

近年来，深度学习技术开始渗透到数字地形分析领域。与传统差分算法不同，深度学习方法可以从DEM中自动学习多尺度地形特征。例如，卷积神经网络（CNN）已被用于从LiDAR点云中直接提取地形特征（如冲沟、陡坎），而生成对抗网络（GAN）则被用于DEM超分辨率重建。

本文评述认为，深度学习在地形分析中的优势在于其多尺度特征学习能力和对噪声的鲁棒性，但当前面临两大挑战：一是训练数据获取困难（高精度地形标签稀缺），二是模型可解释性不足（难以像差分算法那样给出明确的物理含义）。未来，物理约束的深度学习（Physics-Informed Deep Learning）可能是突破方向。

8.2 不确定性量化与传播

地形因子提取中的不确定性来源多样：DEM高程误差、算法近似误差、尺度效应等。这些不确定性会通过后续模型逐级传播，最终影响决策结论。近年来，蒙特卡洛模拟、贝叶斯方法和模糊集理论被引入地形分析的不确定性量化中。

笔者建议，在关键工程应用中，应至少进行以下不确定性分析：

• 基于DEM误差模型的蒙特卡洛模拟，评估坡度、坡向的置信区间。
• 多分辨率、多算法组合的敏感性分析，识别结果的不确定性范围。
• 关键位置的实地验证，校准模型偏差。

九、结语

坡度、坡向和TWI的提取看似是GIS中的基础操作，但其背后涉及的算法选择、尺度匹配和不确定性传播问题，远比表面复杂。本文以"尺度效应"为主线，系统梳理了从DEM预处理到地形因子计算的全链路技术细节，并针对水文、生态和工程三类场景给出了具体的参数配置建议。

本文评述认为，地形因子提取的核心挑战不在于"能不能算"，而在于"算得对不对、合不合适"。在实际项目中，笔者建议遵循"明确目标尺度 → 选择匹配分辨率 → 对比多种算法 → 量化不确定性"的工作流，而非盲目追求高分辨率或复杂算法。随着深度学习与不确定性量化技术的成熟，地形分析正从"静态计算"走向"动态自适应"，这将是未来五到十年最值得关注的方向。