GF7立体像对跑出DSM以后下一步到底该干什么很多人会直接上滤波想一步到位把DEM从DSM里抠出来。前两篇我们把GF7数据准备、区域网平差和立体匹配的过程捋了一遍这篇专门讲一个经常被跳过、但在平原地区DSM/DEM项目里绕不开的环节——可信度分割。所谓可信度分割简单说就是把立体匹配阶段产生的质量信息匹配代价、左右一致性、遮挡标记等整理成一张二值掩膜哪些区域高程可信、可以直接交给后续滤波哪些区域不可信、需要单独标记出来处理。对平原测图来说这一刀切得准不准直接决定后面DEM能不能过精度检核。这篇博文适合正在跑GF7数据、或者打算把国产立体卫星数据真正用进1:10000测图流程的朋友参考我会把原理、参数、代码和踩坑一起倒出来。1. 为什么平原地区制作DEM前必须先做可信度分割1.1 GF7立体测图的主流工作流到底长什么样先把大框架摆出来后面的章节都要在这个框架里说话。GF7的双线阵相机获取的是同轨前后视立体像对生产测绘产品的标准流程大致是L1级影像准备与辐射预处理、区域网平差利用RPC模型加控制点、密集立体匹配生成DSM、对DSM做后处理、滤波得到DEM、最后做DOM和等高线。很多人做到第三步“生成DSM”就认为大功告成直接跳到滤波。这个想法在山地丘陵勉强可以拼一拼但在平原地区几乎是必翻车。因为平原地区的立体匹配结果里混着大量错误高程这些错误和真实地形混在一起滤波算法很难自动分辨。可信度分割就应该夹在DSM生成和滤波中间它的作用是告诉后续处理“哪些像素你有把握哪些像素你少碰。”1.2 平原地区DSM里的噪声和错误到底有多“要命”平原地区看名字平地一片但实际并不是。农田、道路、沟渠、低矮房屋、树木、水塘各种地物交织在一起真正高程变化往往只有几米甚至不到一米。这时候DSM里任何一个0.5米的匹配错误放在山区可能只是地形细节的扰动放在平原就直接淹没了真实地形信号。我印象最深的是有一次处理华北某个灌区DSM生成后肉眼看上去很平滑但切一条剖面出来原本应该是平直田埂的地方全是锯齿。后来一查就是匹配窗口跨越田埂和农田边界时算法在纹理弱区选了错误的视差。平原地区的错误通常不是大范围的空洞而是大片“看似合理但其实是错的”高程这种最坑人因为你不做可信度分割根本不知道哪里该信。1.3 可信度分割在整个流程里的位置滤波前的“划红线”你可以把可信度分割理解成给数据“划红线”可信区域是路面滤波可以放心开不可信区域是雷区要么绕开、要么标记出来做局部处理。这一步不做后面滤波会犯两类错误第一把不可信区域的错误高程当成真实地面滤波后DEM上留下一个莫名其妙的鼓包或凹坑第二为了削掉这些错误把滤波阈值调得过于激进结果可信区域的真实地形细节也被磨平。做了可信度分割你就能让滤波只看它该看的地方精度和细节两头都保得住。2. 可信度数据从哪儿来匹配质量信息来源全解析2.1 立体匹配阶段留下的副产品你都用了吗可信度分割听起来玄乎其实判断依据从你生成DSM那一刻就已经有了只是很多人没注意。主流密集匹配算法在执行时都会产生一些副产品最常见的是代价体积cost volume和最小代价图。以半全局匹配SGM为例算法为每个像素计算“左影像视差为d时的匹配代价”最终取代价最小的视差作为该像素的高程来源。那么最小代价本身就是一个天然的可信度指标代价越小说明这个视差在匹配意义上越确定如果最小代价和第二小代价拉不开差距说明这个地方纹理太弱匹配结果随机性很强可信度就低。在OpenCV的SGBM中你会看到uniquenessRatio这个参数本质上就是控制这个“次小代价和最小代价的差距阈值”。差距不够大直接判为无效。做GF7数据时如果用的是自己的匹配流程一定要把代价体积存下来这个信息比DSM本身值钱得多。2.2 左右一致性检查最基础也最有效的可信度指标在摄影测量里左右一致性检查L-R consistency check是立体匹配后处理的老祖宗也是我目前在实战里最依赖的单一指标。方法是这样的先从左影像出发匹配到右影像得到每个像素的视差d_L再从右影像出发匹配到左影像得到d_R。对左影像某个像素来说如果匹配是正确的d_L应该和它在右影像对应像素处的d_R基本相等误差一般控制在1个像素以内。如果左右两个方向算出来的视差差异超过阈值几乎可以断定这个像素的匹配是错的。这套逻辑在平原地区尤其管用因为平原大量弱纹理区域会表现出“左右不一致”这种区域的错误高程基本都是瞎猜出来的。我在项目里通常的做法是生成两张图一张是DSM另一张是“左右一致性差异图”然后直接把差异大于0.5像素对应约0.3米的高程不确定度的像素丢进不可信区。这招简单粗暴但能干掉大部分飞点。2.3 从DSM自身反算可信度的两个土办法不是每次项目都能拿到匹配算法输出的质量图。比如你手上的数据是别人跑好的DSM栅格只有tif文件加上一个附带的质量波段甚至只有裸高程没有辅助信息。这时候可以用两个土办法反算可信度。办法一是局部高程一致性分析对每个像素取周围一定邻域比如5x5或7x7计算局部高程标准差或局部梯度。如果邻域内高程剧烈跳动大概率是匹配噪声。平原地区真实地形变化平缓局部平滑度应该很高所以这个办法在平原的适用性比山区好很多。办法二是多视冗余检查如果同一区域有多轨数据或多次观测把多个DSM叠在一起逐像素做差。差值大的地方就是匹配不稳定区域。GF7重访周期不算长攒同一区域的两期影像在平原项目里可行。这个办法对找系统性错误比如某个时相整体偏移特别有效但对数据量要求高适合精度要求高的控制测区。3. 可信度分割的实操做法与参数选择3.1 方法一阈值分割加形态学精修最直接的做法就是阈值分割。把匹配代价图、相关系数图、左右一致性图等归一化到0到1区间设置一个阈值超过阈值的算可信低于阈值的算不可信。单用阈值不行分割出来的掩膜会像撒芝麻一样全是孤立点和小洞这时候要上形态学处理。我的标准流程是先做二值化再做一次开运算去掉小噪点再做一次闭运算填掉假空洞最后用一个面积阈值抹掉小于N个像素的碎块。N怎么定看产品分辨率GF7全色分辨率约0.8米如果1平方米左右的碎斑块对最终DEM影响很小就把25个像素以下的图斑全部清掉。这一步很多人忽略但形态学处理对后续滤波的影响很大。如果不清理碎块滤波算法在掩膜边缘会产生大量不必要的过渡带拖慢处理速度还会让边缘高程出现奇怪的内插弧线。3.2 方法二区域增长与种子点策略阈值分割属于“全球一刀切”对某些场景不够聪明。我见过更精细的做法是区域增长先从高可信像素里挑种子点再按邻域相似性向四周生长直到遇到不可信边界才停下来。这个方法的物理意义比较清晰真实地形在空间上是连续的可信区域应该成片出现而不是离散的随机点。区域增长把“像素级可信度”升级成了“区域级可信度”滤掉了那种“单个像素看起来可信、但周围全是错误”的孤立点。具体操作上种子点可以选相关系数高于0.95的像素生长条件是“相邻两个像素视差差小于0.3像素且相关系数大于0.8”。我在有大量低矮建筑群的城乡结合部试过这个方法比纯阈值分割要好因为建筑边缘的半遮挡区域会被自然划分到不可信区而不是留下锯齿状的过渡带。3.3 平原地区阈值怎么标定阈值设置是可信度分割里最玄学的部分但也是有规律可循的。我的经验是不要在办公室里拍脑袋定阈值而是先选一块有典型地物的区域大致建模把DSM和原始影像叠在一起人工判读定一个初始阈值再反复迭代。这里给一组我常用的经验参考值供起步使用具体项目请用自己的检查点验证指标来源归一化方式可信阈值参考备注SGM最小代价代价归一化到0~10.70以上可信代价越低越可信可能需要取倒数左右一致性差异差异像素数小于0.5像素GF7前后视交角固定这个值比较稳相关系数/NCC原始值0~10.75以上可信低纹理平原区要适当放宽到0.6局部高程标准差归一化到影像高程范围的0~1小于0.1具体看地形起伏范围有一点要提醒平原地区相关系数低不一定代表匹配错误因为大片麦田、裸土本身就没什么纹理相关系数天然上不去。这时候不要盲目调高阈值否则会把半个测区都切成不可信后续滤波根本无从下手。正确做法是结合左右一致性来判断只要左右两遍匹配结果一致就算相关系数只有0.6也可以先放进候选可信区。4. 掩膜怎么用进DEM生成流程4.1 生成二值掩膜的Python示例把可信度分割落地的第一步是生成一张与DSM同尺寸、同地理参考的掩膜影像。掩膜取值很简单0代表不可信背景1代表可信前景或者反过来看你后续工具的约定。下面是我常用的Python流程依赖GDAL和SciPy。import numpy as np from osgeo import gdal from scipy import ndimage # 读取DSM和左右一致性差异图 dsm_ds gdal.Open(dsm.tif) lr_ds gdal.Open(lr_consistency.tif) if lr_ds is None: raise SystemExit(缺左右一致性图请先由匹配阶段输出) lr lr_ds.ReadAsArray().astype(np.float32) conf lr_ds.GetRasterBand(1).ReadAsArray().astype(np.float32) # 1. 阈值分割左右差异小于0.5像素才算可信 threshold 0.5 mask (lr threshold).astype(np.uint8) # 2. 开运算去噪点 闭运算填洞 mask ndimage.binary_opening(mask, iterations3).astype(np.uint8) mask ndimage.binary_closing(mask, iterations5).astype(np.uint8) # 3. 抹掉面积过小的碎块 label, num ndimage.label(mask) sizes ndimage.sum(mask, label, range(num 1)) small_label np.where(sizes 25)[0] mask[np.isin(label, small_label)] 0 # 4. 写出掩膜保持与DSM相同的地理参考 driver gdal.GetDriverByName(GTiff) out_ds driver.Create(mask.tif, dsm_ds.RasterXSize, dsm_ds.RasterYSize, 1, gdal.GDT_Byte) out_ds.SetGeoTransform(dsm_ds.GetGeoTransform()) out_ds.SetProjection(dsm_ds.GetProjection()) out_ds.GetRasterBand(1).WriteArray(mask) out_ds.FlushCache() out_ds None # 顺手统计可信区域占比 valid_pixels np.count_nonzero(mask) total_pixels mask.size print(f可信区域占比: {valid_pixels / total_pixels * 100:.2f}%)这段代码不复杂但第3步面积阈值那行值得多说一句我用25像素作为默认值对应GF7全色分辨率下约16平方米的面积。如果你做的是重点测区对微小地物敏感建议把这个值降到9约8平方米如果是大面积概查可以升到100减少碎块干扰。4.2 掩膜进入滤波的两种落地方式掩膜生成以后怎么把它“喂”给滤波算法是另一个容易翻车的地方。不同滤波软件接口不一样但思路主要有两条。第一种方式把掩膜作为输入权重层。很多滤波算法支持权重栅格滤波时会显著降低权重低的点对拟合面的影响。这种方式对渐进加密TIN类算法尤其友好你可以在种子点选取阶段就把不可信区域的点排除掉防止错误点污染地面模型。实测效果是加不加掩膜的TIN加密结果在平原地区能差出20到30厘米的高程误差这个量级在平原已经是致命的。第二种方式先滤波后回填。先用完整DSM滤波得到一版“初步DEM”然后把不可信区域的DEM值遮掉用周围可信区域的趋势面做插值回填。这个方式的优点是能保证最终DEM没有数据空洞缺点是如果不可信区域太大插值相当于“无中生有”精度没有保障。我在实际项目里更推荐组合打法大面积不可信区域水域、建筑群用第一种方式直接排除小面积零散不可信区域用第二种方式补。这样既避免了大范围插值带来的假地形又保证了产品连续完整。4.3 参数组合建议与精度验证手段可信度分割做完、滤波跑完一定要有验证手段不能“看着平滑就交付”。我习惯在测区里均匀布设检查点用实测高程或者高精度参考DOM/DSM做差分统计。平原地区重点关注两件事高程中误差是否达标以及是否出现了局部异常斑块。如果发现滤波后DEM在某个区域出现不应该有的凸包或凹陷优先检查那个区域的掩膜是不是有问题。八成是可信度阈值设得太宽放进了错误匹配点。参数组合我现在的起步配置是左右一致性阈值0.5像素相关系数阈值0.65代价阈值0.7面积碎块阈值25像素形态学开闭运算各3次和5次。这个组合在我手上三个平原项目里都经受住了检验但我不建议你直接照搬因为不同区域、不同季节影像的质量差异很大拿着自己的数据多试几组阈值才是正解。5. 常见问题与排查技巧实录5.1 大范围匹配空洞平原地区大范围匹配空洞一般出现在水体、云影、积雪覆盖区域。水体表面几乎是镜面反射前后视看到的目标差异巨大匹配算法直接投降。云影区域纹理被压扁相关系数低得离谱。这类空洞在可信度分割里会被正确标成不可信但你拿到的空洞面积可能大到影响产品成图。我的处理顺序是先看空洞占测区总面积的比例如果超过10%先别急着补直接调原始影像考虑换一景更新鲜、云影更少的数据重做如果比例不大就把它标记出来交给人工作业员结合实地资料编辑。这里有个小技巧GF7的激光测高点在平原水域区域有时能提供参考高程如果数据包里带了激光测高辅助数据可以在水边线上手动加几个高程约束点再让插值算法去拟合水面效果比纯邻近内插好很多。5.2 建筑物边缘外扩这个是平原地区最容易犯的错。立体匹配算法在建筑物边缘这种深度不连续的地方天然会把建筑物顶面的高程向外“摊”一圈导致建筑轮廓比实际大。可信度分割如果只看匹配代价恰恰会把建筑物边缘的高置信错误当成可信点放进掩膜。我在项目里吃过一次大亏一片厂房区DSM上看轮廓都正常但把高程剖面和实际建筑轮廓一对比整体外扩了大约两个像素。解决办法是在可信度分割之前先用边缘检测算法比如Canny加灰度梯度把影像上的强边缘提取出来然后在掩膜上强制把边缘两侧各拓宽1到2个像素的区域设为不可信。这样虽然损失了一点点建筑物边界地带但换来的是DEM不会被错误高程污染划算。5.3 农田重复纹理导致的系统性错位平原农田最怕的是整齐划一的田垄、播种行。这种重复纹理会让匹配算法产生整行整列的视差偏移而且偏移量正好等于纹理周期。更讨厌的是这种错误匹配的代价可能很低左右一致性也可能正常常规可信度指标根本识别不出来。识别这类问题的土办法是把DSM转换成分辨率较粗的网格比如50米看有没有规律的条带状高程异常或者直接对比DOM人工扫视是否存在“整排房屋跟田埂对不齐”的现象。一旦确认是重复纹理错位就要特别注意不能只靠可信度分割解决需要回到匹配参数层增大匹配窗口、提高视差平滑约束权重从源头减少误匹配。5.4 阈值“一刀切”丢掉地形细节阈值调严了错误是少了但平原上的微小地形特征浅沟、微洼地、老河道遗迹也被一刀切掉了。平原地区这些微地貌往往才是水文分析、耕地平整项目的核心关注对象。我自己的做法是分级处理第一级用极严格的阈值切出一块“绝对可信核心区”作为地形趋势面估计的根基第二级用宽松阈值切出一块“候选可信区”交给滤波算法做参考两级掩膜叠加后再由人工在关键地物附近做局部修正。这套流程虽然麻烦但能把平原地区的微地形从噪声里拯救出来。结尾一个不成熟但很实用的建议做了这么多轮GF7平原数据处理我最想强调的其实是可信度分割不是一道“可做可不做”的工序它是平原地区DSM转DEM成败的分水岭。别嫌它多一步操作、多一份掩膜文件这一小步能帮你省下后面几天的人工编辑时间也能让最终成果在精度检核时少流冷汗。最后分享一个小习惯每次跑完可信度分割我都把掩膜叠加在DOM上出一张预览图存到项目目录里。这样事后一旦发现DEM有问题回头检查是阈值问题、匹配问题还是滤波问题一眼就能定位到具体环节。这个习惯救了我至少两次推荐你也试试。