1. 地形湿度指数到底在算什么地形湿度指数英文全称 Topographic Wetness Index圈子里一般直接叫 TWI。它描述的是在某个地形位置上水流自然汇聚的趋势有多强。你可以把它理解成一张“地形告诉水该往哪儿待着”的地图指数越高说明这个位置越容易积水、土壤越容易饱和指数越低说明水待不住容易排走。这个指标最早来自水文模型领域后来被土壤学、生态学、景观规划、精准农业大量借用。比如做湿地识别TWI 高值区往往和实际湿地范围高度重合做土壤水分空间分布预测TWI 是性价比极高的地形辅助变量做农业排水规划低洼易涝地块用 TWI 一筛就出来了。它最大的好处是只需要一份 DEM 就能算不需要实测土壤水分数据属于典型的“低成本、高信息量”地形衍生指标。TWI 的经典公式是TWI ln( a / tanβ )其中 a 是单位等高线长度上的汇水面积也就是比汇水面积单位通常是 m²/mβ 是局部坡度tanβ 就是坡度正切值。公式本身不复杂但真正动手算的时候坑几乎全在“a 怎么求”和“β 怎么算才合理”这两步上。ArcGIS 里没有一键出 TWI 的工具必须自己搭流程这也是为什么很多人第一次算出来的结果要么全是异常值要么图面看起来完全不对。这篇文章我按实际项目里的操作顺序把从 DEM 准备到 TWI 出图的完整链路拆开讲包括每一步为什么这么做、参数怎么定、报错怎么排。适合已经会用 ArcGIS 基本操作、但没系统做过水文地形分析的人也适合做过一遍但结果总觉得不对劲、想回头查漏补缺的人。2. 整体流程设计与方案选型2.1 为什么 TWI 必须分步算而不能一步到位ArcGIS 的栅格计算器确实能写表达式但 TWI 没法用一个表达式直接算出来原因是它依赖两个中间结果汇水面积和坡度。而汇水面积又不是 DEM 直接能给的必须先做水文分析把洼地填平、算流向、算流量累积才能得到每个栅格的上游汇水栅格数再换算成实际面积。所以整个流程天然分成三大块DEM 预处理、水文分析求汇水面积、坡度计算与栅格运算。这三块之间有严格的先后依赖顺序错了结果一定不对。我见过有人先算坡度再填洼结果坡度是基于原始 DEM 算的和填洼后的流向对不上最后 TWI 图出现大量条带状伪影。标准流程我建议固定成下面这条链路不要随意调换DEM 填洼Fill流向计算Flow Direction流量累积Flow Accumulation坡度计算Slope比汇水面积换算栅格计算器求 TWI异常值处理与出图2.2 汇水面积两种算法怎么选比汇水面积 a 的求法有两种主流路线选哪种直接决定结果量级和空间形态。第一种是单流向法也就是 D8 算法ArcGIS 的 Flow Direction 默认就是它。每个栅格的水只流向八个邻居中坡度最陡的那一个。优点是计算快、逻辑简单、和 ArcGIS 水文工具链完全兼容缺点是水流方向被强制离散成八个方向在平缓地带容易出现平行条带汇水面积分布不够自然。第二种是多流向法代表是 MFD、D-infinity 这类算法。水按坡度比例分配给多个下游栅格汇水面积分布更平滑理论上更接近真实坡面流。但 ArcGIS 原生工具箱里没有现成的多流向流量累积工具要么用 ArcHydro 扩展要么用 SAGA、WhiteboxTools 这类外部工具算完再导回来。我的建议是如果你只是做区域尺度的相对湿度格局分析D8 完全够用别折腾多流向如果你做的是小流域精细水文模拟坡面汇流路径对结果影响很大那值得上多流向。下面实操部分我以 D8 为主线讲因为它是 ArcGIS 里最稳、最不容易报错的路线。2.3 坡度算法与单位陷阱坡度计算本身简单Slope 工具一跑就出来但单位陷阱特别多。TWI 公式里的 tanβ 是坡度的正切值而 ArcGIS 的 Slope 工具默认输出的是度数不是弧度也不是正切值。所以你不能把 Slope 的输出直接代进公式必须先转成正切。转换方式有两种一是在 Slope 工具里把输出单位选成 PERCENT_RISE得到百分比坡度再除以 100 得到 tanβ二是输出度数后用栅格计算器写 Tan(坡度 × 3.1415926 / 180)。两种都行我个人习惯第二种因为度数中间结果更直观方便检查坡度分布是否合理。还有一个容易忽略的点坡度为 0 的平坦区域。tanβ 等于 0 时a / tanβ 会变成无穷大或直接报错。实际数据里填洼后的 DEM 一定存在大片平坦区所以必须对坡度做下限截断比如把小于 0.001 的 tanβ 统一设成 0.001。这一步不做TWI 图里会出现大量 NoData 或者天文数字。3. DEM 预处理与水文分析实操3.1 DEM 获取与投影选择DEM 来源很多常见的有公开的全球 DEM 数据集、无人机航测生成的 DSM 再转 DEM、以及测绘部门提供的本地 DEM。分辨率方面做流域尺度 TWI10 米到 30 米够用做地块尺度最好用 5 米甚至更高精度。分辨率越粗微地形被抹平TWI 的空间细节越少。投影这一步必须重视。TWI 公式里的汇水面积是实际面积单位是平方米所以 DEM 必须是投影坐标系不能是地理坐标系。如果你拿到的是经纬度坐标的 DEM先用 Project Raster 转成适合当地的投影坐标系比如 UTM 或高斯克吕格。用地理坐标系直接算流量累积出来的像元面积是度数换算的量级完全错。提示投影带号选择要覆盖研究区中心经度跨带的话建议先裁剪再分带投影不要整幅硬转。3.2 填洼处理与参数设置填洼用 Hydrology 工具集里的 Fill。输入就是投影后的 DEM输出是填洼 DEM。这里有个关键参数叫 Z 限制默认是空。它的作用是控制填洼深度上限如果两个洼地之间的高差超过这个值就不填。一般小区域分析直接留空让工具把所有闭合洼地填平。填洼之后一定要做一步检查用填洼 DEM 减去原始 DEM得到填洼深度栅格。正常情况大部分区域差值为 0只有洼地位置有正值。如果整幅图大面积出现明显差值说明 DEM 本身有系统性问题比如拼接缝、噪声这时候要先做滤波或人工修正不能硬填。我踩过的一个坑某次用无人机 DSM 直接当 DEM 用建筑物和树冠导致大量假洼地填洼后整个城区被填成一片平地TWI 完全失去意义。后来先用滤波把建筑和植被去掉再做填洼才正常。所以 DEM 一定要是“裸地高程”不是表面高程。3.3 流向与流量累积流向计算用 Flow Direction输入填洼 DEM输出流向栅格。默认算法就是 D8不用改。输出栅格的取值是 1、2、4、8、16、32、64、128分别代表八个方向这是 ArcGIS 的内部编码不用管它具体含义后续工具会自动识别。流量累积用 Flow Accumulation输入流向栅格输出累积栅格。这里有个选项叫输出数据类型有 INTEGER 和 FLOAT 两种。做 TWI 建议选 FLOAT因为后续要参与浮点运算用整型会在换算面积时丢精度。流量累积的结果是每个栅格上游汇入的栅格数量包括它自己。注意流量累积出来的值在河道里会非常大在坡面上很小这是正常的。不要看到数值跨度大就以为出错。3.4 比汇水面积换算流量累积得到的是栅格个数要换算成比汇水面积 a公式是a 流量累积值 × 栅格面积 / 等高线长度对于规则栅格等高线长度近似等于栅格边长。所以简化成a 流量累积值 × 栅格边长举个例子如果 DEM 分辨率是 10 米栅格面积是 100 m²边长是 10 m那么 a 累积值 × 10单位是 m²/m。这一步用栅格计算器写a FlowAcc * 10如果你用的是 5 米 DEM就把 10 换成 5。这个换算看着简单但很多人直接拿累积值当 a 用导致 TWI 整体偏小量级完全不对。4. 坡度计算与 TWI 栅格运算4.1 坡度输出与正切转换坡度用 Slope 工具输入填洼 DEM输出单位选 DEGREE。得到坡度栅格后用栅格计算器转正切tan_beta Tan(slope_deg * 3.1415926 / 180)这里 3.1415926 是圆周率近似值够用。转换完检查一下取值范围正常应该在 0 到 1 之间陡崖位置接近 1平地接近 0。4.2 坡度下限截断前面说过tanβ 为 0 会导致除零。所以下一步做截断tan_beta_safe Con(tan_beta 0.001, 0.001, tan_beta)Con 是条件函数意思是当 tanβ 小于 0.001 时取 0.001否则取原值。这个 0.001 对应的坡度大约是 0.057 度非常平缓不会影响正常坡面的计算但能避免平坦区报错。4.3 TWI 最终计算现在两个输入都齐了用栅格计算器写最终公式TWI Ln(a / tan_beta_safe)Ln 是自然对数。算完之后先看统计值正常 TWI 范围大概在 3 到 15 之间均值 7 到 9 比较常见。如果出现负值或者超过 20 的极端值回去检查 a 和 tanβ 的换算。4.4 异常值处理与重分类实际数据里总会有少量极端值比如河道出口或者数据边缘。处理方式有两种一是按分位数截断比如把小于 1% 分位和大于 99% 分位的值替换成边界值二是直接对 TWI 做重分类分成低、中、高湿度区用于后续制图。重分类用 Reclassify按自然断点或者手动阈值分 5 到 7 类。我一般用自然断点因为它能根据数据分布自动找分界图面层次比较自然。5. 常见报错与排查速查5.1 填洼报错与内存问题填洼最常见的报错是内存不足尤其是大区域高分辨率 DEM。解决办法是先裁剪研究区只保留需要的范围或者分块填洼再合并但合并处要做接边处理。另一个报错是 DEM 存在 NoData 导致填洼中断先用 IsNull 检查并填补。5.2 流量累积全为 0 或异常流量累积全 0 通常是流向栅格有问题检查流向是否成功生成、是否有大片 NoData。如果累积值在平坦区出现异常大的斑块说明填洼不彻底回去重新填洼并加大 Z 限制。5.3 TWI 出现 NoData 或无穷值NoData 基本是除零导致检查 tanβ 截断是否生效。无穷值也是同样原因确认 Con 函数写对且参与运算的栅格范围一致。如果两个栅格范围不一致栅格计算器会按交集处理边缘出现 NoData先用 Clip 或 Resample 对齐。5.4 结果图面出现条带伪影条带伪影是 D8 算法的固有缺陷在平缓区域尤其明显。缓解办法一是提高 DEM 分辨率让微地形更真实二是对填洼后的 DEM 加一点点随机噪声再算流向打破平行流三是改用多流向算法。如果只是做相对格局分析条带不影响整体趋势可以接受。问题现象可能原因排查与解决填洼报内存不足DEM 范围过大或分辨率过高裁剪研究区分块处理流量累积全 0流向栅格异常或全 NoData检查流向生成填补 NoDataTWI 大量 NoData除零tanβ 未截断用 Con 设置 tanβ 下限TWI 量级偏小汇水面积未换算成实际面积累积值乘栅格边长图面条带伪影D8 单流向固有缺陷加噪声、提分辨率或换多流向边缘出现 NoData参与运算栅格范围不一致统一裁剪范围或重采样对齐6. 实操心得与参数经验6.1 分辨率与结果尺度匹配TWI 对 DEM 分辨率非常敏感。同一区域30 米 DEM 算出来的 TWI 平滑5 米 DEM 算出来的细节丰富但噪声也多。选分辨率要看研究目的做流域湿度格局30 米够做田块排水5 米起步。不要盲目追求高分辨率噪声放大后反而干扰判断。6.2 填洼深度检查不能省每次填洼后我都习惯做一次差值检查这一步花不了几分钟但能提前发现 DEM 质量问题。差值图里如果出现大面积非零说明 DEM 有系统偏差继续往下算就是浪费时间。6.3 坡度截断阈值怎么定0.001 是我常用的默认值但不是唯一选择。如果你的研究区特别平坦比如平原农田可以适当放大到 0.005 甚至 0.01避免平坦区 TWI 被过度放大。阈值越大平坦区 TWI 越保守具体取值建议做几次敏感性测试看结果图是否稳定。6.4 出图配色与分级TWI 出图建议用蓝到红的渐变色低值蓝、高值红符合湿度直觉。分级用自然断点5 到 7 类。如果要做对比分析多期 TWI 图必须用同一套分级阈值否则颜色不可比。6.5 批量处理思路如果研究区多、需要批量算 TWI可以把整个流程写成 ArcGIS 模型构建器或者 Python 脚本。核心是把填洼、流向、累积、坡度、栅格计算串成一条链输入输出用参数控制。模型构建器适合不写代码的人Python 适合需要循环批量的场景。我一般用 Python 加 arcpy一个脚本跑几十个流域没问题。最后分享一个我常用的检查顺序算完 TWI 先看直方图再看均值和分位数最后叠加正射影像看高值区是否落在实际低洼处。三步都过了这份 TWI 基本就能用了。