跑高分7号立体像对做平原DSM最烦的不是匹配参数调不对而是匹配跑完、初始DSM出来之后软件给了一堆“置信度”“质量”图层文档里就写一句“数值越大越可靠”然后就没了。我一开始也这么干看一眼关掉直接进滤波。结果滤波越滤越假——水面被拟合成一个斜坡农田被捋成波浪板建筑阴影里冒出一堆十几米高的“山包”。后来把置信度图层重新打开做了分割用掩膜把不可信区域单独拎出来整个流程才顺畅了。这篇是系列第三篇前两篇分别处理了GF-7立体像对的区域网平差准备以及从立体像对到初始DSM的匹配提取。如果你刚看到这里建议先回翻前文的思路因为“可信度分割”这一步吃的东西就是前面匹配阶段吐出来的副产品——置信度栅格。这篇要解决的核心问题是既然匹配结果不是每个像素都可靠怎么把“可靠”和“不可靠”用一种可复核、可复现的方式分开让后续滤波、插值、编辑和质检都有据可依。1. 可信度分割到底在解决什么问题平原测图里的“隐藏粗差”很多人有一个误解觉得匹配失败会直接导致DSM出现明显的洞或巨大的异常值肉眼扫一眼就能看出来。在山区或许是这样但在平原地区这个前提不成立。1.1 先回顾一下前两步从立体像对到DSM的链路GF-7搭载的是双线阵立体相机前视和后视影像构成一个立体像对。要得到DSM基本链路是影像预处理、有理函数模型RFM区域网平差、核线重采样、密集立体匹配、前方交会、DSM栅格化。在前两篇中我们已经把原始影像处理到了“能匹配”的状态并通过双目密集匹配得到了一版初始DSM。这一步的输入是核线影像对输出是一张高程栅格外加一张隐含的、和DSM同分辨率的置信度图。这版DSM里有真实地貌也混入了大量匹配噪声和粗差绝对不能直接拿去当成果。很多商业软件默认不把置信度图层作为必选输出单独展示它在工程内部以“匹配代价”“质量图”或“权重文件”的形式存在。如果你用的软件没有显式导出要回到匹配参数里翻一翻把它导出来。这一步值得做因为它是后续所有质量控制的基础。1.2 为什么平原地区比山地更需要可信度分割平原地区地表起伏小高差常常只有几米到十几米。摄影测量匹配本身的定位精度在亚像元级如果高层影像的基高比不够大几米的真实地形起伏在视差上就只体现为零点几个像素。这时只要匹配代价出现一点偏差DSM就会把几厘米的视差误差放大成几米的高程跳动。更麻烦的是平原地区的地物类型高度重复。大片冬闲田、收割后的麦茬地、温室大棚连栋排列这些地表在影像上呈现周期性的弱纹理。立体匹配算法遇到重复纹理时会在多个候选视差位置得到极其接近的匹配代价即使算法最终选了一个这一位置的事后置信度也一定不会高。但如果你只盯着DSM看这些区域的表现往往不是“明显错误”而是“说不出来的别扭”。比如田块内部出现规律的条纹状起伏或者整片水田被抹成一个略微倾斜的平面与周边的田埂路面高程完美衔接。这种错误是系统性的滤波去不掉只有结合可信度信息才能识别。1.3 可信度信息的三种来源代价、检核与几何条件做可信度分割之前必须搞清楚手里的置信度图是什么含义。不同来源的置信度物理意义完全不同阈值方向和适用场景也不一样。第一种是匹配代价置信度。密集匹配算法在逐个像素寻找同名点时会构建一个代价体积cost volume记录每个候选视差下的匹配代价。代价越小代表匹配越吻合但由于噪声和重复纹理的存在代价最小的位置不一定是对的。有些算法会输出“最优代价”本身或者归一化互相关NCC峰值这类置信度通常越大越可靠。第二种是峰值形态置信度。算法记录最优代价峰值的高度和次峰的比值或者代价曲线谷底的尖锐程度。主峰比次峰高得越多说明匹配越“坚定”这类置信度也是越大越可靠。峰值比接近1说明有二义性这是重复纹理最典型的症状。第三种是几何检核置信度。比如左右一致性检查L-R check分别从左图往右图和从右图往左图匹配比较两遍得到的视差是否一致残差越小可靠性越高再比如多视匹配中的交会射线命中数、姿态稳定度等。这类置信度更适合筛掉遮挡、运动物体和边缘错误。提示在写分割脚本前先弄清楚软件输出的置信度是“越大越可靠”还是“越小越可靠”。搞反了方向后面所有阈值都是白调的。我见过有人拿着匹配代价图当作“可信度”数值越大越觉得可靠正好反了分割出来的掩膜全反。2. 把可信度图读明白格式对齐与物理含义不同软件置信度图叫法不一样像素深度和存储方式也五花八门。这一节先把数据层面的工作做扎实否则后面阈值分割再漂亮也是空中楼阁。2.1 不同软件输出的置信度通道长什么样以PCI Geomatica为例OrthoEngine的DSM提取模块在生成DSM的同时会输出一个“置信度”文件一般每个像元对应一个浮点值代表匹配质量。NASA Ames Stereo PipelineASP的立体匹配结果中会输出“mapping track”文件表示每个DSM像元被几个匹配投票命中命中数多说明该点可靠。国内一些像素工厂、摄影测量集群软件则会输出“质量标记”文件用枚举值标出0、1、2、3这类质量等级。不管叫法如何这一步最该做的不是急着写代码而是先在图像查看软里把置信度图层加载进来用灰度拉伸和伪彩色翻一翻心里有个直观印象。我常用的做法是把置信度栅格和初始DSM用“卷帘”方式对比找十几个典型的错误区域看它们在置信度图上是否呈现低值。2.2 分辨率与像元对齐最容易忽略的第一步置信度栅格和DSM栅格不一定天然对齐。很多软件在导出时会把置信度降采样到DSM分辨率但有些会保持原始影像分辨率甚至因为采用不同的裁剪范围导致两者行列号不一一对应。直接拿两张错位的栅格做分割边界的错位效应会被后续的形态学处理放大。保险的做法是在GDAL里先确认两个栅格的原点坐标、像元尺寸、行数和列数是否一致。不一致时以DSM为准对置信度做最邻近重采样。不要用双线性或三次卷积做重采样因为置信度本身是分类性质的指标插值会产生中间值干扰后续阈值判断。gdalinfo initial_dsm.tif gdalinfo confidence.tif # 分辨率/原点和行列号一致则跳过不一致时执行 gdalwarp -tr 2.0 2.0 -r near -te xmin ymin xmax ymax \ confidence.tif confidence_aligned.tif2.3 统计与可视化在分割之前先“用眼睛看”对齐之后先做像素统计。用直方图看置信度分布通常有三种形态接近正态分布、长尾分布、双峰或多峰分布。正态分布说明匹配整体稳定阈值可以选在尾部长尾分布说明可靠区占大多数需要重点处理长尾的低置信度部分双峰分布最理想谷底就是最佳分割点。看完直方图再叠加影像做个目视排查。可以把置信度低于某个估计阈值的区域用透明色叠加到正射影像上快速看这些区域落在哪里。如果发现它们正好对应水体、阴影、高反射屋面说明分割方向没问题如果大量出现在道路和田块上则要检核几何校正精度和匹配参数设置。这一步不需要复杂工具QGIS的栅格渲染改一下透明度就行或者直接用Python的rasterio读进来转成RGB叠加PNG。3. 阈值怎么定直方图、分位数与地类分块的组合策略阈值是可信度分割的核心参数也是最容易拍脑袋定错的地方。我见过不少人直接填一个经验值比如“NCC小于0.3就当作不可靠”但换一片影像、换一季地表覆盖这个值就不灵了。3.1 不要一上来就二分先看直方图像形二值掩膜方便但如果阈值没把握直接二分会把“不确定区域”强行归入任意一侧后续很难补救。我习惯的做法是先把置信度分成三块准可靠区、准不可靠区、中间缓冲区。具体操作先用累积分布函数CDF找出置信度10%和60%分位数作为临时下阈值和上阈值。低于10%分位数的一定是不可靠区高于60%分位数的可靠进入后续处理中间区域先保留等出现更多上下文信息再判断。这个临时阈值不是最终结果只是用来把分割推到可操作状态。import numpy as np # conf为读入的置信度二维数组valid为有效像元掩膜 p10 np.percentile(conf[valid], 10) p60 np.percentile(conf[valid], 60) # 方向以“越大越可靠”为例 mask_low valid (conf p10) mask_mid valid (conf p10) (conf p60) mask_high valid (conf p60)3.2 经验阈值与自适应阈值的取舍有些匹配软件输出的置信度有明确的物理含义比如相关系数那么用经验阈值是合理的。NCC低于0.5基本可以判定为匹配不可靠但在不同季节的影像上这个值也会漂移。所以哪怕有经验阈值也建议用分位数做一次交叉验证。如果项目范围内有少量已知高程控制点可以先统计控制点所在位置的置信度看它们落在哪个区间。控制点对应的位置如果置信度都很高说明经验阈值可以适当放宽如果控制点位置置信度有高有低说明阈值设置或匹配参数存在问题。3.3 平原场景中的三类典型低置信度地物平原地区低置信度区域通常不多但一旦出现就是成片出现。常见的有三类地物类型置信度失效原因DSM表现分割建议水体无纹理、镜面反射、风成波纹伪纹理高程无规律跳动、整片洼陷应分割为不可靠区并在后续插值中彻底排除冬闲农田/裸地弱纹理、垄线重复周期性微小起伏、条带状异常结合NDVI等地物指数辅助判断不宜全归低置信度建筑阴影及高反射屋面遮挡、亮度饱和、匹配歧义尖峰、凹陷、边缘鼓起分割为不可靠区配合边缘检测扩展缓冲晴好天气下干涸水塘和裸露耕地最容易出问题因为看上去纹理很弱但并不是完全没纹理算法会强行给出一组“看起来平滑”但实际高程偏离的视差。这类错误在DSM上表现出高度一致性仅靠高程统计完全无法发现必须依靠置信度图。4. 分割执行从置信度图到二值掩膜的完整操作流前面工作做完代码层面其实很简单但有一些坑要避开。4.1 基于GDAL和numpy的最小可用脚本下面给一个最小脚本基于rasterio读取置信度图生成二值低置信度掩膜并做基础形态学清理。如果环境里没有rasterio用GDAL命令行gdal_calc.py也能完成类似操作。import numpy as np import rasterio from scipy import ndimage conf_path confidence_aligned.tif mask_path low_confidence_mask.tif threshold 0.50 # 以“越大越可靠”为例按实际统计调整 with rasterio.open(conf_path) as src: conf src.read(1).astype(np.float32) profile src.profile nodata_value src.nodata valid conf ! nodata_value mask np.zeros_like(conf, dtypenp.uint8) mask[valid (conf threshold)] 1 # 形态学清理先开运算去掉孤立像元再闭运算填补缺口 mask ndimage.binary_opening(mask, iterations2).astype(np.uint8) mask ndimage.binary_closing(mask, iterations3).astype(np.uint8) profile.update(dtyperasterio.uint8, count1, nodata0) with rasterio.open(mask_path, w, **profile) as dst: dst.write(mask, 1)这里有几个细节值得展开。阈值到底是0.50还是别的取决于你的置信度定义脚本中的0.50要替换为上一节统计得到的值。形态学操作中开运算的iterations要从小往大试平原地区孤立噪声通常只有一两个像元2次足够闭运算的iterations去填补小洞但不要太大否则会把低置信度区域边缘扩展到正确地形区。4.2 形态学清理去掉椒盐、补齐碎洞二值掩膜直接使用会碰到两个问题一是可靠区内部离散分布着一些零星低置信度像元形似椒盐噪声二是低置信度区域内部有少量孤立的高置信度小洞看起来像被虫蛀过。开运算先腐蚀再膨胀能去掉小噪声块闭运算先膨胀再腐蚀能补掉小洞。关键参数是迭代次数也就是结构元素大小。平原地区的错误匹配往往成片存在零星的低置信度像元基本不值得保留开运算可以放心做。同样低置信度区域内部如果有少量高值像元它们通常是匹配碰巧成功的孤立点对后续插值意义不大用闭运算补掉更干净。但需要小心处理大面积水体边缘。水体的掩膜边界呈锯齿状如果开闭运算次数过多边界会整体收缩一圈导致后续插值时仍有一部分水面像元混入。一个更稳的做法是先做一次3×3中值滤波再做一次开运算边界收缩量更容易控制。4.3 边界缓冲区与数据编辑接口形态学处理之后的掩膜边界往往是硬边界后续如果要做内插或人工编辑硬边界会带来明显的接缝效应。建议在最终掩膜基础上对低置信度区域边界向外扩展3到5个像元作为缓冲区将边界过渡和后续的滤波衔接问题留给插值阶段处理。缓冲区可以用形态学膨胀实现buffer_mask ndimage.binary_dilation(mask, iterations3).astype(np.uint8)有了带缓冲区的掩膜后续的滤波器就可以只在缓冲区外部采用全权重在缓冲区内采用衰减权重避免高程突变突兀。这个接口设计虽然多写几行代码但在真正进入成果编辑阶段时省下的工作量远超预期。5. 分割结果的质量评估误分率和落地检验掩膜生成了不等于工作完成了。一个没有经过检验的分割掩膜在后续使用中一旦出错错误会被放大到最终成果里而且比原始匹配错误更难排查。所以分割完必须做一轮系统性评估。5.1 与影像叠合目视检查要看到什么程度第一关是目视检查。把低置信度掩膜叠加在正射影像上不要只看几个局部而是全图浏览。重点看三个方面掩膜是否覆盖了所有明显水体掩膜是否避开了大片连续田块和道路掩膜边缘是否与地物边界基本吻合而不是横切过明显的结构线。如果某片水域没有被掩膜覆盖说明阈值定得太严如果大片道路和居民地被掩膜覆盖则阈值又太松。这里没有标准答案因为不同项目对最终DEM的需求不同但目视结果必须与地物直觉一致。5.2 定量指标面积占比、边缘密度与误差相关性除了目视还可以用几个定量指标辅助判断。最常用的是低置信度面积占比。平原地区通常低置信度面积占总面积5%到30%是合理区间。如果低于5%大概率阈值太严大量潜在错误区域漏网如果高于40%说明匹配参数本身有问题或者数据时段选择不当比如大面积云影、冬季裸地这时候先不要急着调阈值回到匹配阶段优化参数更重要。另一个指标是掩膜的边缘密度。低置信度区域的边缘越破碎说明掩膜受噪声影响越严重边缘越平滑越符合地物边界特征。边缘密度可以通过边缘像元数除以总像元数来估算也可以简单观察位图形态。最后有条件的话用少量检查点做定量验证。在一个已知地形的小区块用阈值分割后的低置信度掩膜把检查点分成两组分别统计两组的高程残差。如果低置信度组的高程残差明显大于可靠组说明掩膜方向正确如果两组残差没有显著差异说明低置信度掩膜划分得没有区分度需要调整策略。5.3 一条容易翻车的链路把置信度分割当滤波用置信度分割得到的低置信度掩膜用途是引导后续处理而不是直接对其他数据做修改。最容易翻车的操作是拿掩膜直接对DSM做“清洗”——把低置信度区域的高程全部置为NoData或者简单抹成周边均值。这样做表面上看DSM变干净了实际上抹掉了真实的微地形信息。平原地区高程可靠性很大程度上依赖少量弱纹理区域的匹配这些区域虽然置信度偏低但并非全部错误简单抹除会把几十厘米的真实地形起伏一并抹掉。合理做法是把掩膜作为后续插值或滤波的输入条件用可靠像元约束不可靠区域的重建而不是一删了之。6. 后续衔接掩膜如何参与滤波、插值与质检可信度分割不是终点它是整个DSM/DEM生产链条里的一个中间产物。这里把掩膜往后续流程的接口方式梳理清楚。在滤波阶段加权中值滤波可以通过掩膜为每个像素动态调整窗口内的样本权重。低置信度区域内的像元进入统计时被排除或降权可靠像元的权重保持不变窗口内如果可靠像元太少则扩大搜索半径。直接照搬现成软件中的中值滤波参数是不可取的因为软件并不知道你这个掩膜阈值是怎么来的。在插值阶段掩膜定义了“哪里需要重建”。在低置信度区域内做插值时如果只依赖邻域可靠像元容易产生穿越地物边界的“拉链”式伪地形。配合带缓冲区的掩膜可以把插值范围限制在合理区域内并在缓冲区内部做距离权重衰减减少接缝效应。在质检阶段掩膜的价值更大。测绘项目交付时需要说明哪些区域是实测得来、哪些区域是内插重建得来。低置信度掩膜天然就是一张“内插区域分布图”可以直接作为质检和成果说明的附件。对于平原地区大面积水域掩膜甚至可以直接转化为DXF矢量面指导外业补测或说明测绘范围。我个人实操中的体会是GF-7在平原项目里置信度分割对最终DEM精度的贡献经常被严重低估。很多团队以为把匹配参数调好、滤波做细就万事大吉却在检查点比对时发现误差总下不去。一旦把置信度掩膜纳入处理闭环很多“说不清道不明”的误差来源立刻水落石出——就是这个河边滩地的那片低置信度匹配区在作祟。最后再分享一个小技巧如果时间允许建议把低置信度掩膜和基于高程梯度的粗差检测联合使用。先在高程域做粗差检测把明显跳变区域提取出来再用置信度掩膜缩小范围。两个独立证据链相互印证误分割率会明显降低。这一招在我处理春季融冰期影像时特别有用水面和泥滩交替分布单靠任何一侧都极易出偏。