1. 引言
随着遥感测量技术与空间信息获取手段的不断发展，基于多源数据融合的地表液态水储存量调查已成为水资源精细化管理的重要技术途径。在调查实施过程中，数字高程模型(DEM)作为描述地表及水下地形特征的基础空间数据，是开展水体容积计算、水位变化分析及统计评价的重要依据，其数据质量直接关系到调查成果的可靠性与应用价值。
在地表液态水调查中，通常需要综合利用水下地形测量获取的点云数据与机载激光雷达等遥感手段获取的陆地区域DEM数据，构建水陆一体DEM。然而，由于数据来源多样、获取方式差异明显，加之测量环境复杂，在DEM生产与融合过程中不可避免地会引入噪声点及异常高程值[1]。已有研究指出，多源遥感数据在融合过程中易受数据精度差异和空间一致性影响，从而对DEM数据质量产生不利影响，这类粗差若未得到有效识别和处理，不仅会降低DEM的整体精度，还可能在后续水体储量计算中被放大，从而影响最终调查成果的科学性。
针对DEM数据质量控制问题，国内外学者已在滤波处理、异常检测及误差分析等方面开展了大量研究。其中，基于空间滤波的粗差检测方法因实现简便、适用性强，被广泛应用于DEM数据预处理与质量检查环节[2]。然而，现有方法多侧重于单一滤波手段，难以兼顾复杂地形条件下的精度保持与异常识别需求，且在实际业务应用中仍以人工判读或半自动处理为主，效率与稳定性有待进一步提升。
基于此，本文结合辽宁省地表液态水储存量调查的实际需求，在分析水陆一体DEM数据特征的基础上，研究了一种基于多尺度滤波分析的DEM粗差检测方法，并依托Python编程环境开发了相应的自动化处理工具。通过将多种滤波方法与残差分析相结合，实现DEM粗差的自动识别与定位，为地表液态水调查数据质量控制提供一种可操作性较强的技术方案。
2. 研究方法
2.1. DEM粗差特征分析
DEM粗差通常表现为明显偏离周围地形变化规律的异常高程值，其形成原因主要包括测量误差、传感器噪声、数据拼接不一致以及人为处理失误等[2]。在水陆一体DEM中，由于水陆交界区域地形变化复杂，粗差问题尤为突出。因此，有必要采用适应性较强的检测方法对异常高程进行识别。
2.2. 多尺度滤波检测原理
滤波技术是DEM数据处理中的常用方法，其基本思想是通过对邻域像元进行空间平滑，削弱随机噪声影响，同时突出异常信息。本文在综合分析不同滤波方法特性的基础上，引入均值滤波、中值滤波、高斯滤波和双边滤波等多种方式，对DEM数据进行多尺度处理，以提高粗差检测的完整性和可靠性。
均值滤波：以邻域像元高程均值替代中心像元，去除整体的、平缓的随机噪声，运算高效，适用于平坦区整体降噪，但易弱化地形边缘[3]。
中值滤波：基于邻域排序统计，用邻域内所有像元的中间值替代中心像元，抑制孤立尖峰/凹陷，保边缘性能优于均值滤波，是平原/丘陵区常用最佳滤波[4]-[6]。
高斯滤波：根据高斯函数对邻域像元进行加权平均，距离越远权重越小。采用高斯加权核实现渐进平滑，去噪自然，适合丘陵区地形过渡区处理[7]-[9]。
双边滤波：同时考虑空间邻近度与高程相似度进行加权平均，在保留边缘的同时平滑噪声，严格保护陡坡、山脊、沟谷，效果最接近真实地形，是山地/高山的最佳保边滤波[10]。
形态学滤波：通过开/闭运算剔除伪地形，对LiDAR-DEM中建筑物、植被干扰抑制效果显著，可作为多地形组合滤波环节[11] [12]。
2.3. 粗差检测算法设计
综合滤波处理与残差分析，本文构建的DEM粗差检测流程主要包括：数据导入、多尺度滤波、残差计算、阈值判别以及结果输出等步骤，实现了DEM粗差检测的自动化处理。
3. 工具开发实现
基于上述方法，本文采用Python 3.7作为主要开发语言，结合arcpy与OpenCV等类库，实现DEM粗差检测功能的集成。工具在设计过程中充分考虑批量处理需求，引入图形化界面以提升操作便捷性，可完成DEM数据从导入到结果输出的全流程处理。本文设计的粗差检测核心算法包括以下几个步骤：
1) 滤波处理：首先通过多种滤波方法对原始DEM数据进行平滑处理。常见的滤波方法包括均值滤波、中值滤波、高斯滤波、双边滤波和形态学滤波。不同的滤波方法对噪声的去除效果不同，因此选择合适的滤波方法是确保后续粗差检测准确性的基础。
2) 残差计算：通过计算滤波后的DEM数据与原始DEM数据之间的残差(差异)，可以有效地突出高程的异常值。残差的计算是通过对比原始数据和滤波后的数据来完成的，即：
Residual = OriginalDEM − FilteredDEM
残差的大小直接反映了DEM数据中高程值的异常程度，较大的残差通常代表较大的粗差。
3) 阈值判定：通过设定一个阈值，将计算出的残差与阈值进行比较，识别出超过阈值的部分作为粗差。阈值的选择直接影响到粗差检测的准确性。本文实验中按照地形类别设置了不同的阈值(平原0.8米、丘陵1.0米、山地2.4米，高山地3.0米)，以适应不同地形条件下的粗差检测需求。
4) 结果输出：最终，通过对超出阈值的残差进行处理，识别出粗差位置并输出相应的结果。这些粗差可以以栅格和矢量(点)格式进行保存，便于后续分析和修正。
5) 主要代码如下：
def process_one(self, ras, out_filter, out_diff, out_shp, th, win, ftype):
dem = Raster(ras)
desc = arcpy.Describe(ras)
base = desc.baseName
ext = "." + desc.extension if desc.extension else ""
# 选择滤波方法进行处理
if ftype == "mean":
fs = FocalStatistics(dem, NbrRectangle(win, win, "CELL"), "MEAN", "DATA")
elif ftype == "median":
fs = FocalStatistics(dem, NbrRectangle(win, win, "CELL"), "MEDIAN", "DATA")
else:
arr = arcpy.RasterToNumPyArray(dem, nodata_to_value=None)
if ftype == "gaussian":
arr_f = cv2.GaussianBlur(arr, (win, win), 0)
elif ftype == "bilateral":
arr_f = cv2.bilateralFilter(arr.astype("float32"), win, 50, 50)
origin = arcpy.Point(dem.extent.XMin, dem.extent.YMin)
fs = arcpy.NumPyArrayToRaster(arr_f, origin, dem.meanCellWidth, dem.meanCellHeight, dem.noDataValue)
# 计算残差
diff = dem - fs
absd = Abs(diff) # 取残差的绝对值
diff_masked = SetNull(absd < th, diff) # 按阈值去除小于阈值的差异
# 输出残差栅格
out_d = os.path.join(out_diff, "D_" + base + ext)
diff_masked.save(out_d)
# 输出粗差点SHP
shp = os.path.join(out_shp, base + "_outlier.shp")
arcpy.RasterToPoint_conversion(diff_masked, shp, "Value")
4. 应用实例
4.1. 研究区及实验数据
实验选取辽宁省地表液态水储存量调查中的水陆一体DEM生产过程数据，数据陆地部分采用机载激光雷达数据获取的DEM成果，成果为1:2000比例尺，空间分辨率为5 m；水下DEM没有做要求。
数据成果采集要求：
平原：高程中误差0.40米，最大误差0.80米(粗差 ≥ 0.8米) [13]；
丘陵地：高程中误差0.50米，最大误差1.00米(粗差 ≥ 1.0米) [13]；
山地：高程中误差1.20米，最大误差2.40米(粗差 ≥ 2.4米) [13]；
高山地：高程中误差1.50米，最大误差3.00米(粗差 ≥ 3.0米) [13]。
研究区为辽宁省西部，范围包括朝阳市、葫芦岛市、锦州市和盘锦市，面积约42,556平方千米，地形类别包括平原、丘陵、山地和高山地。
4.2. 滤波参数设置
结合实验数据的地形特征，设置不同滤波方式及窗口参数，采用多尺度策略对DEM进行平滑处理。具体参数如下：
1) 平原DEM
2) 丘陵DEM
3) 山地DEM
4) 高山地DEM
4.3. 检测结果分析
4.3.1. 程序粗差筛查情况
通过对研究区不同地形类别设置相应滤波、阈值和窗口，程序筛查平原202个粗差、丘陵738个粗差、山地4610个粗差和高山地6155个粗差，程序筛查合计11,705个粗差。
平原高程差统计见图1，丘陵高程差统计见图2，山地高程差统计见图3，高山地高程差统计见图4。
Figure 1. Statistical chart of elevation differences in plain areas
图1. 平原高程差统计图
Figure 2. Statistical chart of elevation differences in hilly areas
图2. 丘陵高程差统计图
(a) 2.4米 ≤ 粗差 < 20米统计 (b) 20米 ≤ 粗差统计
Figure 3. Statistical chart of mountain elevation differences
图3. 山地高程差统计图
(a) 3米 ≤ 粗差 < 20米统计 (b) 20米 ≤ 粗差 < 60米统计
(c) 60米 ≤ 粗差统计
Figure 4. Statistical chart of elevation differences in alpine regions
图4. 高山地高程差统计图
4.3.2. 人工排查
排查方法主要是利用同期高分辨率卫星影像，对比滤波前后DEM和滤波前垂直剖面，对全部粗差像元进行人工核对，修正程序筛查的粗差数据，排除统计误报粗差点。
通过人工排查，共筛查出研究区误判粗差点559个，多分布在地貌破碎的采矿场区域；粗差点大部分为树冠点分层错误，主要分布在山地和高山地；小部分粗差点为独立建构筑物和噪声点分层错误。
4.3.3. 实验结论
通过不同阈值条件下的实验对比分析，结果表明本文方法在不同粗差幅度条件下均具有较高识别精度，整体检测效率较传统方法提升约40%以上。其中：
1) 平原：阈值0.8米下识别粗差点202个，均为有效粗差；
2) 丘陵：阈值1.0米下识别粗差点738个(其中粗差709个，误判29个)；
3) 山地：阈值2.4米下识别粗差点4610个(其中粗差4394个，误判216个)；
4) 高山地：阈值3.0米下识别粗差点6155个(其中粗差5841个，误判314个)。
与传统检测方法相比，本文方法在保证识别精度的同时，整体处理效率提升约42.5%，粗差检出率超过95%。
5. 结束语
本文围绕地表液态水储存量调查中水陆一体DEM数据质量控制问题，研究实现了基于多尺度滤波的DEM粗差检测方法并开发了自动化处理工具。
1) 本文以高程差为核心指标，以排查超规范最大误差像元为滤波核心目标，构建平原、丘陵、山地、高山地四大地形的DEM粗差滤波参数体系，成果完全满足规范精度要求；
2) 最佳滤波方法随地形显著变化：平原以中值滤波最优，丘陵以高斯滤波最优，山地与高山地以双边滤波为最佳，可针对性抑制粗差；
3) 滤波窗口随地形复杂度逐级增大：从平原3 × 3~5 × 5到高山地9 × 9~13 × 13，适配地形表达与滤波强度；
实验结果表明，该方法在精度与效率方面均具有较好表现，可有效控制DEM噪声，同时保留真实地形结构，适用于DEM生产、质量检查与成果交付，可为测绘、水利、国土等领域高精度应用提供技术支撑。