3. DEM质量检查与空洞修复

在构造线性体提取流程中，DEM 不只是"底图"，它直接决定坡度、坡向、山脊线与谷线等地形派生量的可靠性。SRTM、ASTER GDEM、ALOS World 3D 等产品在陡崖、水体、云覆盖区、雷达阴影区常出现无值像元（NoData），在城区与林区则可能出现"尖峰""凹坑"等异常。若不做质量检查就直接提取线性体，这些缺陷会被边缘检测算子放大成虚假的线性构造。因此本节把质量检查与修复作为一条独立工序来处理：先量化缺陷，再选择修复策略，最后用指标验证修复是否"改对了地方"。

一、数据空洞的识别

空洞识别本质上是"有效值掩膜"的构建。对单波段浮点 DEM，先统计全图有效像元比例，再生成二值掩膜：有效为 1，NoData 或超出合理高程范围为 0。合理高程范围应结合研究区设定，例如青藏高原东部可设 -100～8000 m，超出即视为异常。随后对掩膜做连通域分析，得到每个空洞斑块的面积、周长、位置与外接矩形。这一步很关键：面积小于 3×3 像元的孤立空洞可用邻域插值直接填补，而面积达数百像元、呈条带状分布的空洞往往对应雷达阴影或水体，需要结合辅助数据判断是否应当填补——水体本就不该有地形高程，强行插值反而制造假地形。

实际检查中建议同时输出三类统计：空洞像元总数与占比、空洞斑块数量及面积分布直方图、空洞在坡度分级上的分布。若空洞高度集中在坡度大于 35° 的区域，说明数据源在陡坡失相干，修复时应优先采用能保持坡面连续性的方法，而不是简单均值填充。

二、异常值检测

异常值分两类。第一类是"孤立尖峰/凹坑"，表现为某像元与其 8 邻域均值差异极大。可用 3×3 窗口计算局部中值 M 与局部标准差 σ，若 |z − M| > 3σ 且该像元与邻域高差超过阈值（如 50 m），标记为可疑。第二类是"系统性条带或台阶"，常见于拼接产品，表现为沿某一方向的高程突变。可用方向梯度检测：计算 x、y 方向一阶差分，若某一行或列的差分均值显著偏离全图均值（如超过 5 倍中误差），则判定为接缝异常。

# 示例：基于局部中值的异常检测（参数为教学假设值）
win = 3
M = median_filter(dem, size=win)
S = std_filter(dem, size=win)
resid = dem - M
mask_spike = (abs(resid) > 3 * S) & (abs(resid) > 50)
print("可疑尖峰像元数:", mask_spike.sum())

检测到异常后不要急于全部替换。建议先叠加正射影像或等高线核对：若尖峰位于建筑密集区，可能是真实地物（如高层建筑在 DSM 中的表现），此时应保留并在后续线性体提取中通过尺度滤波抑制；若位于自然山坡且与邻域明显不连续，则判定为噪声并纳入修复。

三、插值修复方法的选择

修复方法需与空洞形态匹配。小面积、各向同性空洞优先用反距离加权（IDW）或局部多项式；条带状、沿坡面延伸的空洞宜用各向异性插值或坡面线性外推，以保持坡向连续；大面积空洞若周边地形复杂，可用薄板样条（TPS）或基于辅助 DEM 的融合填补。下面这张图给出了一条从缺陷识别到修复验证的处理链，以及不同缺陷类型对应的策略分支。

以 IDW 为例，其实现要点是：搜索半径取空洞等效半径的 2～3 倍，幂指数 p 取 2，参与计算的已知点不少于 8 个；若邻域点不足，则扩大半径直至满足。TPS 则需控制平滑参数 λ，λ 过小会过拟合噪声，过大则地形被过度平滑。教学示例中可先用 10% 的已知像元人为挖空，再用不同方法回填，比较回填值与真值的差异，从而为真实空洞选择参数——这是在没有真值情况下最实用的参数标定思路。

四、质量评估指标

修复是否合格，不能只看"空洞被填满了"。建议从三个层面评估。第一是数值精度：在有真值或人工挖空样本时，计算 RMSE、平均绝对误差 MAE 与偏差 Bias，RMSE 反映整体误差量级，Bias 反映系统性偏高或偏低。第二是地形一致性：比较修复前后坡度、坡向的分布，若修复区坡度直方图与周边未修复区显著不同，说明插值引入了不自然坡面。第三是应用级验证：用修复前后 DEM 分别提取线性体，统计线性体密度、走向玫瑰图与长度分布的变化；若修复后线性体数量突增且方向杂乱，往往意味着插值产生了伪边缘。

完成本节工序后，DEM 应达到"无异常尖峰、空洞填补合理、坡度连续、线性体提取结果稳定"的状态，才可进入后续的增强与线性体提取环节。