1. 从DSM到DEM建筑物点过滤到底在解决什么问题拿到一幅遥感图像生成的DSM数字表面模型你会发现它跟DEM数字高程模型最大的区别就一个词地表附着物。DSM记录的是地表最上层的反射面树冠、桥梁、建筑、甚至高压线塔全都算在里面DEM只关心裸地高程。做洪水淹没分析、坡度坡向计算、地形断面提取你拿DSM去算结果全被建筑物带偏。所以建筑物点过滤这件事本质上是把DSM里属于人工建筑的像元识别出来、剔除掉再通过插值填补回裸地高程最终得到一张干净的DEM。我最早接触这个需求是在做城市内涝模拟的时候。当时手里有一批0.5米分辨率的DSM直接拿去做汇水分析结果小区里的水全“汇”到了楼顶上出来的淹没图跟实际完全对不上。后来花了整整两周时间做建筑物点过滤才把数据修正过来。这个教训让我意识到DSM建筑物点过滤不是可选项而是所有基于DSM做地形分析的前置必做步骤。适合看这篇内容的人大概分三类一是做遥感数据处理的技术人员手里有DSM但不知道怎么清洗二是做GIS分析的朋友被建筑物高程干扰搞得头疼三是做遥感图像语义分割的算法工程师想把自己的分割结果用到DSM后处理上。不管你是哪一类下面这些实操细节和踩坑经验应该都能帮到你。2. 建筑物点过滤的整体思路与方案选型2.1 为什么不能一刀切建筑物点的特征分析建筑物点在DSM上有什么特征我总结下来主要是三条高程突变、区域连通、边缘规则。高程突变是指建筑物边缘处DSM值会有一个明显的跳变从地面高程直接跳到楼顶高程区域连通是指同一栋楼的楼顶高程基本一致形成一块连通的“高地”边缘规则是指建筑物的轮廓通常比较规整不像树冠那样边缘模糊。但问题在于不是所有高地都是建筑物。山丘、土堆、立交桥、大型广告牌这些在DSM上也会表现为局部高地。所以过滤算法不能只看高程阈值必须结合形状、面积、纹理等多维度特征。我试过单纯用高程阈值做过滤结果把城市里的小土坡全删了后来不得不加回人工检查环节效率极低。2.2 三条主流技术路线对比目前做建筑物点过滤主流方案可以归为三类方案类型核心原理优点缺点适用场景形态学滤波通过开运算、闭运算等形态学操作分离地物实现简单计算快参数敏感对复杂建筑群效果差地形起伏小的平坦城区渐进式形态学逐步增大结构元素尺寸模拟地面点筛选对地形适应性强计算量大大区域处理慢山地城市、地形复杂区域语义分割辅助用U-Net等网络分割建筑物再映射到DSM精度高可区分建筑与树木需要标注数据模型训练成本高有标注样本的城市区域我个人的建议是如果只是做小区域快速处理形态学滤波够用如果是生产级项目语义分割辅助方案更靠谱。特别是现在基于U-Net的遥感图像语义分割已经比较成熟拿分割结果做建筑物掩膜再叠加到DSM上做过滤精度比纯几何方法高出一大截。2.3 方案选型的核心考量因素选方案的时候我一般会问自己四个问题数据分辨率是多少0.5米和2米分辨率的DSM建筑物特征差异巨大。高分辨率下建筑边缘清晰可以用形态学低分辨率下建筑跟地面混在一起必须上语义分割。地形复杂度如何平原城市和山地城市完全是两个难度级别。山地城市的建筑物往往依山而建高程变化本身就大单纯靠高程阈值根本分不开。有没有标注数据如果有建筑物轮廓的矢量数据或者标注好的遥感图像语义分割方案可以直接上如果没有要么先标注一批要么退回到形态学方案。精度要求多高做科研发论文精度要求高值得花时间做语义分割做工程预研快速出结果更重要形态学滤波先跑一版看看效果。3. 核心细节解析从DSM到建筑物掩膜的完整链路3.1 数据预处理别急着过滤先把DSM洗干净很多人拿到DSM就直接开始做过滤结果发现异常值一大堆。DSM预处理这一步绝对不能省。我一般会做三件事第一去除无效值。DSM里经常有-9999或者NaN的像元这些是数据缺失区域。处理方式很简单用周围有效像元的均值填充或者直接用OpenCV的inpaint函数做修复。第二平滑滤波。原始DSM往往有噪声特别是基于图像匹配生成的DSM在纹理弱的区域会有明显的“椒盐噪声”。我一般用高斯滤波或者中值滤波做一次平滑窗口大小取3×3或5×5。注意不要用太大的窗口否则建筑物边缘会被模糊掉。第三坡度计算。坡度是后续过滤的重要辅助特征。建筑物边缘的坡度通常很大而平地坡度接近0。用GDAL的gdaldem slope命令可以快速生成坡度图gdaldem slope input_dsm.tif slope.tif -p -s 1.0这里的-p表示输出为百分比坡度-s 1.0是垂直单位与水平单位的比例尺。如果你的DSM高程单位是米水平坐标也是米这个参数就设1.0。3.2 基于形态学滤波的建筑物点初筛形态学滤波的核心思想是用一个“结构元素”在DSM上滑动通过开运算先腐蚀后膨胀去掉比结构元素小的凸起。建筑物通常比结构元素大所以会被保留树木、车辆这些小凸起会被去掉。具体操作步骤确定结构元素尺寸。这个尺寸应该略大于最大树木冠幅但小于最小建筑物面积。城市区域一般取15×15到25×25像元。我通常先用几个不同尺寸试跑看哪个尺寸能把树木去掉但保留建筑物。执行开运算。用OpenCV的morphologyEx函数import cv2 import numpy as np dsm cv2.imread(dsm.tif, cv2.IMREAD_UNCHANGED) kernel cv2.getStructuringElement(cv2.MORPH_RECT, (21, 21)) opened cv2.morphologyEx(dsm, cv2.MORPH_OPEN, kernel)计算差值。原始DSM减去开运算结果得到的就是“凸起物”的高程diff dsm - opened building_mask diff 2.0 # 高程差大于2米的认为是建筑物这里的2.0米是个经验阈值。我试过1.5米和3.0米1.5米会把一些矮墙也框进来3.0米会漏掉一些低矮建筑。2.0米在城市区域比较稳妥但具体项目还是要根据当地建筑高度分布调整。3.3 基于U-Net语义分割的建筑物精细提取形态学滤波有个硬伤它分不清建筑物和土堆。土堆在DSM上也是凸起形态学滤波会把它当成建筑物删掉。这时候就需要语义分割出场了。基于U-Net的遥感图像语义分割与地物面积估算系统核心思路是用标注好的遥感图像训练一个U-Net网络让它学会区分建筑物、道路、植被、水体等地物。然后把分割结果映射到DSM上只过滤建筑物类别的像元。U-Net的结构这里不展开讲网上资料很多。重点说一下训练数据的准备标注工具我用的是LabelMe标注建筑物轮廓导出JSON再转成掩膜图。样本数量至少500张512×512的切片覆盖不同建筑类型高层、多层、棚户区、工业厂房。数据增强随机旋转、翻转、亮度调整把样本量扩到2000张以上。类别平衡建筑物像元通常只占10%-20%需要做类别加权否则网络会偏向预测背景。训练好之后推理阶段输出的是建筑物概率图取阈值0.5得到二值掩膜。然后把这个掩膜叠加到DSM上building_mask (prob_map 0.5).astype(np.uint8) dsm_filtered dsm.copy() dsm_filtered[building_mask 1] np.nan # 建筑物像元置为无效3.4 建筑物点剔除后的DEM插值重建建筑物像元被置为NaN之后需要插值填补回裸地高程。常用的插值方法有三种反距离加权IDW简单快速适合平坦区域。GDAL的gdal_fillnodata命令可以直接用gdal_fillnodata.py -md 100 -si 3 input.tif output.tif薄板样条TPS适合地形起伏区域但计算量大。普通克里金精度最高但需要拟合变异函数操作复杂。我一般先用gdal_fillnodata快速填一版看看效果。如果填出来的DEM在建筑物边缘有明显的不自然过渡再换TPS或者克里金。注意插值窗口不要太大否则会把远处的地形特征带进来导致建筑物位置出现虚假的坡面。4. 实操过程一套完整的建筑物点过滤流程4.1 环境准备与工具选型我平时用的工具链是这样的GDAL处理栅格数据的基础库命令行工具和Python绑定都好用。OpenCV做形态学运算和图像滤波。PyTorch segmentation_models_pytorch训练U-Net做语义分割。QGIS可视化检查中间结果和最终成果。安装命令conda create -n dsm_filter python3.9 conda activate dsm_filter conda install -c conda-forge gdal opencv pytorch torchvision -y pip install segmentation-models-pytorch4.2 完整处理流程与参数设置下面是我在实际项目中跑通的一套流程以某城市0.5米分辨率DSM为例第一步DSM预处理from osgeo import gdal import numpy as np import cv2 # 读取DSM ds gdal.Open(city_dsm.tif) dsm ds.ReadAsArray().astype(np.float32) nodata ds.GetRasterBand(1).GetNoDataValue() # 无效值处理 dsm[dsm nodata] np.nan mask np.isnan(dsm) dsm_filled dsm.copy() dsm_filled[mask] np.nanmean(dsm) # 中值滤波去噪 dsm_smooth cv2.medianBlur(dsm_filled.astype(np.float32), 3)第二步形态学滤波初筛kernel_size 21 kernel cv2.getStructuringElement(cv2.MORPH_RECT, (kernel_size, kernel_size)) dsm_opened cv2.morphologyEx(dsm_smooth, cv2.MORPH_OPEN, kernel) diff dsm_smooth - dsm_opened building_mask_morph (diff 2.0).astype(np.uint8)第三步U-Net语义分割精细提取import torch from segmentation_models_pytorch import Unet model Unet(encoder_nameresnet34, encoder_weightsimagenet, in_channels3, classes2) model.load_state_dict(torch.load(best_model.pth)) model.eval() # 假设已经准备好了遥感图像切片 # 推理得到建筑物概率图 prob_map model.predict(image_tensor) building_mask_unet (prob_map 0.5).astype(np.uint8)第四步掩膜融合与建筑物点剔除# 形态学结果和U-Net结果取交集减少误判 building_mask_final building_mask_morph building_mask_unet # 剔除建筑物像元 dsm_filtered dsm_smooth.copy() dsm_filtered[building_mask_final 1] np.nan第五步DEM插值重建# 用gdal_fillnodata插值 driver gdal.GetDriverByName(GTiff) out_ds driver.Create(dem_filled.tif, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32) out_ds.SetGeoTransform(ds.GetGeoTransform()) out_ds.SetProjection(ds.GetProjection()) out_band out_ds.GetRasterBand(1) out_band.WriteArray(dsm_filtered) out_band.SetNoDataValue(np.nan) out_ds None # 调用gdal_fillnodata import subprocess subprocess.run([gdal_fillnodata.py, -md, 100, -si, 3, dem_filled.tif, dem_final.tif])4.3 结果验证与精度评估过滤完之后怎么知道效果好不好我一般做三个检查目视检查在QGIS里把原始DSM、建筑物掩膜、最终DEM叠在一起看。重点看建筑物边缘有没有残留的“台阶”以及插值区域有没有明显的异常坡面。高程剖面检查在建筑物密集区域拉一条剖面线对比过滤前后的高程变化。理想情况下过滤后的剖面应该是一条平滑的地面线没有楼顶的凸起。定量指标如果有参考DEM比如从LiDAR点云生成的DEM可以计算RMSErmse np.sqrt(np.nanmean((dem_final - dem_ref) ** 2))我做过的一个项目形态学滤波单独用的RMSE是1.2米加上U-Net语义分割后降到0.6米效果提升很明显。5. 常见问题与排查技巧实录5.1 建筑物边缘残留“台阶”怎么办这是最常见的问题。建筑物被剔除后边缘处插值出来的高程跟周围地面不连续形成一圈“台阶”。原因通常是建筑物掩膜比实际建筑物大了一圈把边缘的地面像元也误删了。解决办法对建筑物掩膜做一次腐蚀操作把边缘收缩1-2个像元kernel_erode cv2.getStructuringElement(cv2.MORPH_RECT, (3, 3)) building_mask_eroded cv2.erode(building_mask_final, kernel_erode, iterations1)这样虽然会漏掉一点建筑物边缘但能保证地面像元不被误删插值出来的DEM更自然。5.2 树木被误判为建筑物怎么处理树木在DSM上也是凸起形态学滤波很容易把大树冠当成建筑物。U-Net语义分割可以区分树木和建筑物但前提是训练数据里要有足够的树木样本。如果不想训练模型可以用形状特征做后处理建筑物的形状通常比较规整长宽比接近1:1到1:3之间而树冠的形状更接近圆形或不规则。计算每个连通区域的紧凑度面积/周长²建筑物紧凑度一般在0.5以上树冠在0.3以下。from skimage.measure import label, regionprops labeled label(building_mask_final) for region in regionprops(labeled): compactness region.area / (region.perimeter ** 2) if compactness 0.4: building_mask_final[labeled region.label] 05.3 大区域处理内存不够怎么优化0.5米分辨率的DSM一个城市可能有好几GB。直接读进内存做形态学运算很容易爆内存。我的做法是分块处理把DSM切成512×512的瓦片瓦片之间留32像元的重叠区。每个瓦片单独做形态学滤波和语义分割。拼接时只取每个瓦片的中心区域重叠区用来消除边界效应。tile_size 512 overlap 32 for i in range(0, height, tile_size - overlap): for j in range(0, width, tile_size - overlap): tile dsm[i:itile_size, j:jtile_size] # 处理tile # 拼接时去掉overlap区域5.4 常见问题速查表问题现象可能原因排查方法解决方案建筑物边缘有台阶掩膜过大误删地面像元叠加掩膜和DSM目视检查对掩膜做腐蚀操作树木被误删形态学滤波分不清树木和建筑检查掩膜中的圆形区域加形状特征过滤或U-Net分割插值后DEM有异常坡面插值窗口太大检查插值区域的高程梯度减小插值窗口或用TPS插值大区域处理内存溢出一次性读入整个DSM监控内存使用分块处理瓦片间留重叠建筑物漏删高程阈值太高统计建筑物高度分布降低阈值或改用语义分割5.5 几个我踩过的坑坑一直接用原始DSM做形态学滤波。原始DSM噪声大形态学滤波会把噪声当成小凸起去掉结果建筑物边缘变得锯齿状。一定要先做平滑滤波。坑二U-Net训练数据只用一个区域的样本。我试过用A城市的样本训练直接用到B城市精度掉了一半。不同城市的建筑风格差异很大训练数据一定要覆盖多种建筑类型。坑三插值后不做检查直接出图。有一次插值参数设错了建筑物位置插出来一个“坑”比周围地面低了3米。后来养成了习惯每次插值完都要拉几条剖面线检查。坑四忽略坐标系统。DSM的坐标系如果是地理坐标系经纬度坡度的单位会不对。一定要先投影到投影坐标系比如UTM再做坡度和形态学运算。6. 进阶技巧把语义分割和传统方法结合起来6.1 为什么纯语义分割也不够U-Net语义分割虽然精度高但有个问题它只看图像纹理不看高程。有些建筑物在遥感图像上跟地面颜色接近分割网络会漏掉有些地面上的深色斑块会被误判成建筑物。所以纯语义分割的掩膜直接用到DSM上会有不少误判。我的做法是把语义分割结果和形态学滤波结果做融合语义分割结果作为主掩膜提供建筑物的大致范围。形态学滤波结果作为辅助提供高程突变信息。两者取交集既保留了语义分割的精度又利用了高程信息做验证。6.2 融合策略的具体实现# 语义分割掩膜 mask_unet (prob_map 0.5).astype(np.uint8) # 形态学掩膜 mask_morph (diff 2.0).astype(np.uint8) # 融合语义分割为主形态学为辅 mask_final mask_unet mask_morph # 对于语义分割置信度高的区域放宽形态学约束 high_conf prob_map 0.8 mask_final[high_conf] mask_unet[high_conf]这样处理之后建筑物漏删和误删的情况都少了很多。我实测下来融合后的F1分数比单独用U-Net提高了约8个百分点。6.3 后处理让DEM更自然过滤完建筑物点之后DEM上还是会留下一些不自然的痕迹。我一般再做两步后处理第一步高斯平滑。对插值后的DEM做一次轻度高斯滤波窗口取3×3sigma取0.5。这样能消除插值带来的微小震荡让DEM更平滑。第二步坡度约束。计算DEM的坡度如果某个像元的坡度超过35度自然地面很少有这么陡的坡说明这里可能是插值异常用周围像元的均值替换。slope calculate_slope(dem_final) steep_mask slope 35 dem_final[steep_mask] np.nan dem_final fill_nodata(dem_final)这两步做完DEM基本就能直接用于后续分析了。7. 写在最后的一些个人体会做DSM建筑物点过滤这些年最大的感受是没有一种方法能通吃所有场景。形态学滤波在平原城市快准狠但到了山地城市就歇菜U-Net语义分割精度高但需要大量标注数据而且换一个城市就要重新训练。实际项目中我基本都是组合拳先用形态学快速跑一版看效果如果精度不够再上语义分割最后用形状特征和高程约束做后处理。另外参数没有绝对的最优值。同样是2米的高程阈值在高层住宅区可能偏低在别墅区可能偏高。我一般会先统计一下研究区域内建筑物高度的分布直方图取一个能覆盖80%建筑物的阈值剩下的20%靠语义分割和后处理来补。最后分享一个小技巧如果你手头有建筑物的矢量轮廓数据比如从OpenStreetMap下载的可以直接把矢量转成栅格掩膜叠加到DSM上做过滤。这比任何算法都准而且速度快。没有矢量数据的时候再考虑用形态学或语义分割。这个思路我在好几个项目中都用过效果很稳。