亚栅格地形缺位:流向分配的三重瓶颈
DEM 以规则栅格表达地表起伏,是流域水文模拟、土壤侵蚀评价、地貌过程分析、洪水风险识别与水土资源管理的基础。其中流向确定与流量分配直接决定汇流累积量、单位汇水面积(SCA)、河网提取与流域边界划分结果。开源 DEM 覆盖广、获取便利,已广泛用于大范围水文分析,但其分辨率常为 30 m 甚至更低,单个栅格难以表征内部亚栅格尺度的局部沟谷、坡脚过渡与微地形起伏,如何把栅格内部地形信息引入流向分配,成为关键问题。
单流向算法的离散方向偏向
D8 依据最大坡降方向把水流全部分配给一个下坡邻格,原理简单、计算高效,但只能沿 8 个离散方向之一流动,容易产生方向偏向与水流过度集中。Fairfield 与 Leymarie 提出的 Rho4/Rho8 引入随机因素以削弱系统偏差,却因结果具有随机性,在需要稳定重复计算的定量分析中存在局限。
分流参数难以适配复杂坡面
Freeman 的 FD8 按坡降关系计算分流权重,采用固定分流指数,在沟谷等汇流区易出现水流过度扩散;Quinn 等使分流指数随汇流面积变化,Holmgren 比较了不同分流指数的影响,秦承志等的 MFD-md 引入随局部地形条件变化的自适应分流指数。但分流参数仍多依赖经验设定或特定规则,对不同坡面形态适应能力有限。
栅格内部地形差异被弱化
现有方法虽采用最大坡降、随机扰动、坡降加权、自适应分流指数或连续流向角等不同机制,流向确定与分流比例计算总体上仍建立在中心栅格及邻域栅格高程关系之上。当分配仅依据中心点与邻域高程判断时,栅格内部地形对水流路径的控制作用被弱化,较少考虑栅格内部不同位置连续下坡路径对周围邻格分配比例的影响。
栅格内部差异难以表达
D8 方向编码体系
凸心坡·凹心坡·鞍形坡
考察离散尺度影响
降尺度—连续追踪—比例统计:五步技术链路
三角面连续追踪法与依据中心栅格和邻域栅格高程差、坡降关系或其推导参数确定分配比例的思路不同,它把原始栅格内部大量降尺度栅格的连续流动去向作为分配比例的统计依据。方法先以原始栅格中心点为高程控制点构建新格网并三角剖分,再由局部平面方程为降尺度栅格赋值,随后沿三角面最陡下降方向连续追踪每个起点的首次流出方向,最终以方向频次占比作为该栅格向 8 邻格的流向分配比例。
中心点新格网构建
将原始 DEM 每个栅格的高程视为该栅格中心点高程,以原始栅格中心点为顶点构建中心点新格网。该新格网相对原始格网整体偏移半个栅格,但所有顶点均来自原始 DEM 已知的栅格中心点,因此每个顶点的平面坐标与高程值都是已知量,无需引入额外的插值假设即可作为降尺度地形表达的控制基础。
三角剖分与对角线选取
把新格网单元剖分为 2 个三角面,使降尺度起点能沿连续地形表面逐段传播。定义两条候选对角线高程差 d₁₃=|z₁-z₃|、d₂₄=|z₂-z₄|:d₁₃>d₂₄ 时选取 P₁P₃,d₂₄>d₁₃ 时选取 P₂P₄;二者相等则比较低端点高程,选更低的一条;4 个顶点高程完全相等时视为平地格网,默认采用左上—右下方向对角线。优先选择高程差较大的对角线,可使三角面更好地反映局部主要高程变化方向。
降尺度栅格高程赋值
每个三角面由 3 个已知高程点唯一确定一个局部平面,因而三角面内具有明确的最陡下降方向。对每一个降尺度栅格中心点,先判断其落入哪个三角形内部,再依据该三角面对应的局部平面关系对其高程赋值;对全部降尺度栅格重复该过程,即可构建整个研究区域的降尺度高程场,形成原始栅格内部的连续地形表达。
三角面连续追踪与首次流出方向判定
以降尺度栅格中心为潜在流动起点,沿所在三角面最陡下降方向延伸路径,路径与三角面边界相交处即为碰撞点。若碰撞点仍位于所属原始栅格内部,则以该点为新起点切换至相邻三角面,按新平面方程重新确定最陡下降方向继续追踪;若碰撞点已位于原始栅格外部,则停止追踪,按碰撞点所在原始栅格相对中心栅格的方位,判定该起点的 8 邻域首次流出方向并编码记录。
方向频次统计与汇流累积、SCA 计算
设原始栅格内含 N 个有效降尺度起点,流向第 i 个方向的起点数为 n_i,则该栅格第 i 个方向分配比例 p_i=n_i/N,i=1,2,…,8。输出仍是原始栅格尺度上流向周围 8 邻格的比例值,而非降尺度小栅格的连续流向角。在统一框架下按 TCA_i=a+Σp_ji·TCA_j 逐格累积,再取 SCA_i=TCA_i/r 得到单位汇水面积,用于与 D8、MFD-md、D-inf 的对比评价。
三类数学曲面与两组起伏参数
人工表面由明确数学函数构建,地形形态与理论水流路径可由函数关系确定,可解析求解任意位置的理论 SCA 作为参考值,避免实际 DEM 误差干扰评价。选取凸心坡(椭球面)、凹心坡(逆椭圆面)与鞍形坡(双曲抛物面)3 类曲面,凸心坡与凹心坡取 a=300 m、b=400 m,平面范围 800 m×600 m;鞍形坡取 a=400 m,平面范围 800 m×800 m。起伏参数设 c=20 m 与 c=70 m 两种,分别代表低起伏与高起伏;每种起伏下均生成 5、10、20 m 分辨率的人工 DEM,评价指标为 RMSE、ME 与 SD。
凤阳县样区与低洼区重叠评价
实际样区位于安徽省凤阳县部分区域,以平缓低丘地貌为主,山地、坡地、沟谷、农田及局部水体共同分布;DEM 为 30 m 分辨率 ASTER GDEM V3,遥感影像用于结果对照。4 种方法流向结果统一转换为 8 邻域分配比例,在相同汇流累积框架下计算汇流累积量,河道提取的汇流面积阈值统一设为 45 000 m²,本文方法最小降尺度单元取 1 m。评价借鉴 TPI 思想识别低洼区,以邻域平均高程与中心栅格高程之差大于 0.5 m 判定,进而统计低洼区覆盖率 R_D、河道—低洼区一致性 P_D 与调和平均综合指标 F_D。
鞍形坡误差骤降与河道集中性平衡
人工 DEM 实验以理论 SCA 为参考,比较 D8、MFD-md、D-inf 与三角面连续追踪法在 3 类曲面、2 种起伏与 3 种分辨率下的 RMSE、ME 与 SD;实际样区则在统一的 45 000 m² 汇流阈值下比较河道提取结果,并用低洼区重叠的三项指标定量刻画河道与低洼汇流区的空间对应关系。
低起伏条件最优
相对 MFD-md
高于 D-inf 与 D8
1 222 200 m²
| 方法 | 5 m 分辨率 RMSE / m | 10 m 分辨率 RMSE / m | 20 m 分辨率 RMSE / m | 20 m 分辨率 SD / m |
|---|---|---|---|---|
| D8 | 616.3 | 414.6 | 285.5 | 280.7 |
| MFD-md | 176.4 | 149.0 | 127.6 | 114.3 |
| D-inf | 324.6 | 229.4 | 166.3 | 160.0 |
| 三角面连续追踪法 | 59.2 | 32.2 | 20.7 | 16.3 |
| 方法 | 凸心坡 RMSE / m(20 m) | 凹心坡 RMSE / m(20 m) | 鞍形坡 RMSE / m(20 m) | 起伏条件 |
|---|---|---|---|---|
| D8 | 58.0 | 351.6 | 285.5 | 低起伏 c=20 m |
| MFD-md | 15.7 | 256.6 | 127.6 | 低起伏 c=20 m |
| D-inf | 7.6 | 228.1 | 166.3 | 低起伏 c=20 m |
| 三角面连续追踪法 | 7.1 | 205.3 | 20.7 | 低起伏 c=20 m |
| MFD-md | 10.4 | 254.6 | 115.0 | 高起伏 c=70 m |
| D-inf | 7.6 | 228.1 | 166.3 | 高起伏 c=70 m |
| 三角面连续追踪法 | 7.1 | 205.3 | 20.7 | 高起伏 c=70 m |
| 方法 | 低洼区面积 / m² | 河道面积 / m² | 重叠面积 / m² | R_D 低洼区覆盖率 | P_D 河道—低洼区一致性 | F_D 综合指标 |
|---|---|---|---|---|---|---|
| D8 | 4 752 000 | 1 102 500 | 805 500 | 0.170 | 0.731 | 0.275 |
| MFD-md | 4 752 000 | 1 441 800 | 979 200 | 0.206 | 0.679 | 0.316 |
| D-inf | 4 752 000 | 1 159 200 | 859 500 | 0.181 | 0.741 | 0.291 |
| 三角面连续追踪法 | 4 752 000 | 1 222 200 | 885 600 | 0.186 | 0.725 | 0.296 |
复合地形中 RMSE 与 SD 同步大幅下降
鞍形坡汇聚、发散与流向转换并存,各方法误差差异最为明显。低起伏条件下本文方法在 5、10、20 m 分辨率下的 RMSE 分别为 59.2、32.2 和 20.7 m,较同条件下表现最好的对比方法 MFD-md(176.4、149.0、127.6 m)分别降低约 66.4%、78.4% 和 83.8%;对应 SD 为 59.1、32.9 和 16.3 m,较 MFD-md 分别降低约 65.3%、76.5% 和 85.7%,说明栅格内部连续追踪能显著改善复合地形的流向分配。
优势集中于较低分辨率与汇聚地形
凸心坡上本文方法在 5 m 分辨率的 RMSE 为 20.3 m,高于 MFD-md 的 8.7 m 与 D-inf 的 11.6 m,但在 10 m 和 20 m 分辨率下降至 6.9 m 和 7.1 m,接近或略优于 D-inf;凹心坡上 10 m 与 20 m 分辨率的 RMSE 为 264.9 m 与 205.3 m,均低于 MFD-md 的 297.8、256.6 m 与 D-inf 的 273.6、228.1 m。高分辨率简单光滑地形中,中心栅格与邻域高程关系已足以表达坡面流向,本文方法的额外优势并不突出。
河道范围与低洼覆盖的综合平衡
实际样区中 MFD-md 的 R_D 最高(0.206),但其河道面积达 1 441 800 m² 且 P_D 最低(0.679),表现出汇流扩散特征;D-inf 的 P_D 最高(0.741)而覆盖能力不足。本文方法河道面积 1 222 200 m²,较 MFD-md 减少约 15.2%,P_D 达 0.725,F_D 为 0.296,高于 D-inf 的 0.291 与 D8 的 0.275,综合指标较二者分别提高约 1.7% 与 7.6%。
本文方法并非依靠扩大河道范围提高低洼区重叠程度,而是在控制河道扩张的同时保持较好的低洼汇流区覆盖能力与河道空间集中性。
把栅格内部地形写进分配比例
针对开源低分辨率 DEM 难以表达栅格内部地形差异、进而影响流向分配、汇流累积与河道提取的问题,本文设计了基于 DEM 降尺度三角面连续追踪的地表多流向算法,并以数学曲面实验与凤阳县实际样区实验,与 D8、MFD-md、D-inf 系统对比验证其有效性。
以内部连续路径统计替代邻域高程关系
方法以原始栅格中心点为高程控制点构建新格网,经三角剖分与局部平面方程获得栅格内部的降尺度地形表达,再以降尺度栅格中心为起点沿三角面最陡下降方向连续追踪,依据各起点首次流出所属原始栅格的 8 邻域方向统计分配比例,使原始栅格尺度的分配比例具有更明确的内部地形依据。
数学曲面验证:复杂复合地形适应性
在凸心坡、凹心坡与鞍形坡 3 类人工地形中整体表现稳定:与 D8 相比有效减弱了 8 邻域方向限制与单一路径累积带来的误差;与 MFD-md 和 D-inf 相比,在较低分辨率条件下误差控制能力更好,鞍形坡的 RMSE 与 SD 均明显低于其他方法,对汇聚、发散与流向转换并存的复合地形适应性突出。
真实地形验证:兼顾集中性与覆盖能力
在 30 m 分辨率 ASTER GDEM 条件下,4 种方法主干河道格局总体一致,本文方法低洼区覆盖率、河道—低洼区一致性与综合指标分别为 0.186、0.725 和 0.296,综合指标高于 D-inf 与 D8,同时河道面积较 MFD-md 减少约 15.2%,在减弱河道过度扩散的同时保持了与低洼汇流区的空间对应关系。
开源低分辨率DEM水文分析的新思路
该方法将原始栅格内部地形差异对水流去向的影响引入流向分配过程,在一定程度上弥补开源低分辨率 DEM 难以表达亚栅格尺度地形差异的不足,可为该类数据条件下的汇流累积、SCA 计算与河道提取提供新的方法思路。其输出的仍是原始栅格尺度上流向 8 邻格的比例值,可直接接入既有汇流累积框架,无需改动下游水文分析流程;输出的多方向差异化分配格局在平缓坡面趋于分散、在沟谷低洼带趋于汇聚、在山脊地带自然分化,与区域地形走势保持一致。
适用效果受分辨率与评价指标制约
方法优势与数据分辨率密切相关:在较低分辨率和鞍形坡等复杂复合地形条件下优势明显,而在高分辨率、简单光滑地形条件下,中心栅格及邻域高程关系已能较好表达坡面流向,传统算法亦可取得较小误差,本文方法的额外优势不一定突出。实际样区中亦非所有单项指标均最优,优势主要体现在河道空间集中性与低洼汇流区覆盖能力之间的综合平衡。后续可进一步考察降尺度单元尺度、三角剖分规则与低洼区判定阈值等参数设置的敏感性,并在更多地貌类型与分辨率的开源 DEM 上开展适用性检验。