做航测项目最怕什么?外业辛辛苦苦测了一堆控制点,内业空三一跑,平面精度超限,整体精度被拉低,你还得像个挨个查——到底是哪个控制点坏了事。这个排查过程我用Agisoft Metashape手工做过太多次,先看Ground Control面板,逐个人工核对每个点的误差,再回到照片上找刺点的位置,效率低不说,碰到几十个控制点的大测区,眼睛都快看花。后来我干脆用Metashape的Python API写了一套控制点粗差探测脚本,专门干这件事:把控制点重投影误差全部抽出来,统计筛选,直接告诉我哪个点有问题、问题出在哪几张照片上。今天就把这套思路和源码完整分享出来,给正在被粗差折磨的同行省点时间。
这篇文章适合两类人:一类是经常用Metashape做无人机航测、地形图生产的从业者,另一类是刚接触Metashape Python API想搞明白控制点对象模型怎么操作的人。文中代码不需要你懂多深的Python,能在Metashape自带控制台跑起来就行,我会把每一段逻辑干什么、为什么要这么写讲清楚。
1. 控制点粗差问题为什么值得专门写代码解决
1.1 粗差是怎么产生的
粗差这个词在摄影测量里是"错误误差"的代名词,它和平常说的偶然误差完全是两码事。偶然误差是测量过程中不可避免的微小波动,比如GPS定位的cm级抖动、刺点时的亚像素偏差,它们的分布是有规律的,平差算法能消化掉。粗差就不同了——它是带方向性的系统性错误,一个粗差控制点就能把空三结果拉偏一大截。
我在实际项目里总结过,控制点粗差主要来自三个环节。第一个是外业采集环节,RTK测量时忘记等待固定解、杆高输入错误、GPS信号受遮挡引起的跳变,都会让坐标高程偏离真实值。第二个是内业刺点环节,影像分辨率高,同名点在照片上很小,一不小心就把点刺到了地物的另一侧,和相邻影像的实际位置差了十几甚至几十个像素。第三个是影像本身的问题,比如无人机拍摄的运动模糊、大阴影区域的弱纹理、地物高程突变造成的投影差,这些都容易让人误判特征点的位置。
粗差值如果顺利通过平差,问题非常隐蔽。因为Metashape空三时,控制点作为约束条件参与整体解算,单个粗差点会被分配到整个测区的相机位姿误差里,看起来好像每个点的误差都有一点超标,但没有任何一个点表现出"异常"。初学者遇到这种情况,往往会怀疑整个控制网布设有问题,实际上只需要剔除一两个点,整体精度立刻恢复正常。
1.2 Metashape自带的误差面板够用吗
Metashape在Ground Control面板里提供了每个控制点的三维误差显示,估计算法也会给出每个控制点的平方根误差。但说实话,这个面板在实际项目中只能作为参考,不能当作粗差探测的主要手段。
原因有几个。第一,面板显示的是控制点在X、Y、Z三个方向上的误差合成值,反映的是控制点整体坐标差异,但并没有告诉你这个误差到底是坐标测量错了还是影像刺点错了。第二,当控制点本身参与空三优化的过程中,部分粗差会被"吸收"到相机位姿里,导致面板显示的误差偏小,看起来一切正常。第三,当控制点多的时候,人工肉眼看几十上百行数字,很难形成统计意识,往往只能抓住最大的那一两个点,漏掉次大的。
另外,Metashape在优化之后把控制点当作固定值处理,它的误差统计并不会主动告诉你"某个点需要删除"。如果你反复优化、反复看误差,每次结果都差不多,但精度就是上不去,你的第一反应应该是有粗差,而不是去动平差参数。
1.3 用脚本处理的核心优势
写脚本做粗差探测,核心不是自动化,而是把"误差"拆到足够细的粒度去看。Metashape面板只给你一个汇总表格,脚本可以把每个控制点在每一张照片上的重投影误差全部抽出来,单独分析。
重投影误差的概念很简单:控制点有一个三维坐标,根据相机内外方位元素,这个三维坐标理论上应该投影到影像的某个像素位置,而你在影像上实际刺点的位置和理论位置的差,就是重投影误差,单位是像素。如果某个控制点在所有照片上的重投影误差普遍很大,说明是空间坐标有问题。如果只在某几张照片上大,说明那几张照片的刺点位置标错了。这个精细度是面板做不到的。
脚本还会做统计判别:计算所有控制点平均误差的均值和标准差,那些落在3倍标准差之外的点,大概率是粗差。这样比单纯看极值更可靠,因为粗差往往隐藏在"次大值"里,而统计方法能从全局分布中识别出离群的个体。
2. 粗差探测的原理与算法设计
2.1 重投影误差:最可靠的粗差判据
在写任何代码之前,得先把"重投影误差"这个核心指标吃透。Metashape在做空三平差时,本质上是在求解一组相机位姿和三维点坐标,使得所有观测值的重投影误差总和最小。重投影误差越小,说明模型和影像观测越匹配。
对控制点来说,重投影误差的计算路径是:控制点三维坐标 → 由控制测量得到,经过坐标变换到chunk局部坐标系 → 用相机内参、畸变参数和外方位元素投影 → 得到理论像素坐标 → 与影像上实际刺点坐标求差。这个差值包含了两部分信息:控制点坐标本身的误差,以及刺点操作的误差。恰好,这两部分正是粗差的全部来源。
有一个细节比较关键:Metashape里camera.project()函数的输入参数是chunk局部坐标系下的三维点,而控制点的position属性一般存储在chunk的坐标系下(可能是地理坐标、本地投影坐标)。所以必须先做坐标变换,不能直接用point.position塞进camera.project()里。很多自己写脚本的人踩的第一个坑就是这个,后面源码部分我会给完整的转换逻辑。
2.2 判定规则:统计阈值与绝对阈值双保险
粗差判定的经典方法是3σ法则:假设正常控制点的平均重投影误差服从近似正态分布,那么超过平均值3倍标准差的值就可以判为异常。但这个法则有一个前提——数据中的异常值不能太多。如果测区里有一半控制点都有问题,均值和方差本身就被污染了,再用3σ法则可能一个都识别不出来,或者把所有点都识别成异常。
所以我采用了双阈值策略。第一是统计阈值,也就是mean+3*std;第二是绝对阈值,直接指定一个像素上限,比如3像素。最终判定规则是:只要超过任何一个阈值,就标记为候选粗差。这样做的好处是,即使整体误差水平普遍偏高(比如影像质量差、相机标定不完美),也能通过绝对阈值兜底,防止把所有点都当成可疑点。
绝对阈值的设置没有统一的万能值,它取决于测区的地面分辨率和项目精度要求。举个例子,1:500地形图航测,飞行高度设计的地面分辨率是3.5cm/px,如果要求平面精度5cm,那么3.5cm的1个像素误差就已经占了很大比例,阈值定在2~3像素比较合理。如果只是做正射影像拼接,30cm的分辨率下,一个像素的误差也就30cm的平面位移,阈值可以放宽到5像素甚至更多。这部分我会在常见问题里再展开。
2.3 区分点级粗差和影像标记粗差
脚本输出的不该只是一个"有问题"的结论,更好的设计是区分两种粗差类型。第一种是点级粗差:某控制点的所有影像重投影误差整体偏大,误差值相差不大且都远高于其他点。这种情况基本可以确定是控制点三维坐标测量错误,比如坐标记错、杆高错误,需要回到外业重新测量。
第二种是影像级粗差:某个控制点在大部分照片上误差都很小,但在某两张照片上误差突然显著增大,其他照片正常。这种情况说明影像刺点位置错了,可能是刺到了相似纹理上,也可能是粗刺。修复方法是回到对应的照片,在影像上重新精确定位控制点,而不是去动三维坐标。
两种粗差的处理方式不同,源码里必须把每个点在不同影像上的误差明细输出出来,只给一个均值报告会掩盖掉影像级粗差的信息。我在脚本里专门写了detail_errors()函数,传入控制点标签后,能把该点在所有照片上的误差按从大到小排列,方便定位问题出在哪张照片。
3. 准备环境与理解API对象模型
3.1 版本兼容和运行环境
Metashape的Python API在不同版本之间有过调整,最典型的是控制点组的属性名。Metashape 2.x系列使用chunk.control_points作为控制点组的入口,而早期1.x系列还叫chunk.ground_control_points。为了让脚本跨版本可用,代码里我对两种情况做了兼容处理,抛出AttributeError时自动切换备选接口。
运行环境方面,Metashape自带Python解释器,不需要单独安装外部Python环境,直接打开Metashape,在Tools菜单里找到Python Console,把代码粘贴进去运行就行。但要注意,Metashape的Headless模式(命令行批处理)和图形界面模式下,控制台访问方式略有差异。如果用的是命令行批处理,需要在脚本开头多写一段加载项目文件的逻辑;图形界面下直接用Metashape.app.document就能拿到当前打开的项目。我下面给的代码默认支持图形界面运行,也比较容易改造成批处理模式。
关于Metashape的授权,建议使用正式授权版本。脚本本身不需要额外库,只用到了Metashape模块和Python标准库math、csv,没有任何外部依赖,这对我来说很重要。在客户现场或者没有网络的环境下,装第三方库很麻烦,能不依赖就不依赖。
3.2 控制点相关API对象模型
写代码前,把Metashape里控制点相关的对象层级理清楚。顶部是Metashape.Document,对应项目文件,里面包含一个或多个Chunk(数据块)。每个Chunk包含影像、传感器、相机位姿、密集点云等数据,控制点组挂在Chunk下。
控制点组是ControlPoints对象,它的入口属性在不同版本叫points或者items,返回一个点列表。每一个点是一个Point对象,里面有label(点号)、position(三维坐标)、enabled(启用状态)、selected(选中状态)、projections(投影列表)等属性。其中projections是重投影误差计算的关键——它包含该控制点在所有可见照片上的刺点信息。
每个投影对象是Projection,关键属性有两个:camera(所属相机对象)和coord(影像上的像素坐标)。camera对象有project()方法,可以把三维点投影到二维像素坐标;还有sensor.width和sensor.height两个属性记录影像宽度和高度。这些API组合起来,正好能支撑我们的粗差探测算法。
# Metashape对象层级速查 doc = Metashape.app.document # 文档对象 chunk = doc.chunk # 当前chunk cp = chunk.control_points # 控制点组 point = cp.points[0] or cp.items[0] # 第一个控制点 projection = point.projections[0] # 该点的第一个投影 camera = projection.camera # 投影对应的照片 coord = projection.coord # 照片上的像素坐标4. 源码实现与逐段讲解
4.1 主流程设计
整个脚本的主流程分四步。第一步是获取当前chunk的控制点组,并兼容旧版本API。第二步是遍历每个控制点,对每个点遍历它所有投影,计算重投影误差并汇总。第三步是调用统计判别函数,用双阈值策略找出候选粗差。第四步是输出报告,把每个点的平均误差、最大误差、标准差、可疑点名单全部打印出来,支持导出CSV存档。
代码中我把核心的误差计算函数和统计判别函数拆开,是为了方便在Metashape控制台交互调用。第一次运行脚本输出汇总报告后,如果你对某个可疑点想细看是哪个照片的问题,直接在控制台调用detail_errors(chunk, "点号")就能得到逐影像误差明细。
4.2 坐标转换细节
坐标转换是整套代码里最容易出问题的地方。Metashape控制点的position属性存储的坐标系取决于chunk.crs的设置。如果chunk的CRS是UTM等投影坐标系,position就是(E, N, H);如果CRS是WGS-84经纬度坐标,position会是纬度、经度、高程的组合。而camera.project()需要的输入是chunk局部坐标系的三维点,这个坐标系由chunk.transform.matrix定义。两者通常不一致,所以必须转换。
转换路径是:先用chunk.crs.unproject(position)把地理坐标换成地心直角坐标(ECEF,单位米),再用chunk.transform.matrix.inv().mulp(world)把地心坐标转到chunk局部坐标。这里有一个细节是chunk.transform.matrix可能为None,比如项目还没有执行过对齐操作,这种情况我就直接返回世界坐标,让后续投影计算自己去判断。另外,chunk.crs也可能为None,处理方式类似。
还有一点要提醒,camera.project()返回的理论像素坐标可能与projection.coord的坐标系原点定义存在差异。Metashape内部会做图像坐标系归一化,理论上两者的参考点是一致的,我在多个版本的实测中都能直接相减取模。如果遇到异常情况,可以通过camera.sensor.width和camera.sensor.height换算来校验。
4.3 误差统计与判定
误差统计的实现思路是:对每个控制点,先算所有投影误差的平均值、最大值、标准差。然后在整体层面,再对所有控制点的平均误差求一次均值和标准差。这里有一个统计学上的小陷阱——当某个点的误差特别大时,整体均值和标准差都会被拉高,导致原本正常的点显得符合分布,原本超差的点反而被隐藏。所以我在统计判别函数identify_outliers()里保留了两个维度的判断,并且把统计阈值和绝对阈值分开输出,方便你自己判断当前数据到底属于"整体正常、个别异常"还是"整体混乱、需要重新检核"。
关于标准差计算,我用了总体标准差公式,除以n而不是n-1。因为控制点数量通常只有几个到几十个,属于样本统计而不是全量调查,但为了简单起见,总体方差的计算方式也能满足粗差判别需求。追求严谨的可以改成样本标准差。
4.4 完整源码:可直接在Metashape控制台运行
# -*- coding: utf-8 -*- """ Agisoft Metashape 控制点粗差探测脚本 适用版本:Metashape 1.x / 2.x 运行位置:Metashape Python Console 或通过 Tools > Run Script 执行 """ import math import csv import Metashape def get_control_points(chunk): """获取控制点组的兼容包装,支持不同版本API""" if chunk is None: return None try: # Metashape 2.x return chunk.control_points except AttributeError: try: # Metashape 1.x return chunk.ground_control_points except AttributeError: print("当前Metashape版本不支持控制点访问") return None def to_chunk_local_coord(chunk, point): """把控制点的坐标转换为chunk局部坐标,供camera.project()使用""" if point.position is None: return None try: if chunk.crs is not None: world = chunk.crs.unproject(point.position) else: world = point.position # 注意:只有transform.matrix存在时才做逆变换 if chunk.transform is not None and chunk.transform.matrix is not None: local = chunk.transform.matrix.inv().mulp(world) return local return world except Exception as e: print("坐标转换失败: {}".format(e)) return None def calc_point_errors(chunk, point): """计算单个控制点在所有影像上的重投影误差""" local_pos = to_chunk_local_coord(chunk, point) if local_pos is None: return None img_errors = [] for proj in point.projections: cam = proj.camera # 跳过未启用的相机和没有外方位元素的相机 if cam is None or not cam.enabled: continue if cam.reference is None or cam.reference.location is None: continue theory_coord = cam.project(local_pos) if theory_coord is None: continue observed_coord = proj.coord err = (observed_coord - theory_coord).norm() img_errors.append((cam.label, err, observed_coord, theory_coord)) if len(img_errors) == 0: return None avg_err = sum(e for _, e, _, _ in img_errors) / len(img_errors) max_err = max(e for _, e, _, _ in img_errors) std_err = math.sqrt( sum((e - avg_err) ** 2 for _, e, _, _ in img_errors) / len(img_errors) ) return { "label": point.label, "avg": avg_err, "max": max_err, "std": std_err, "n": len(img_errors), "detail": img_errors, } def detect_all_errors(chunk): """遍历所有控制点,返回误差统计列表""" cp = get_control_points(chunk) if cp is None: return [] try: pts = cp.points except AttributeError: pts = cp.items if len(pts) == 0: return [] results = [] for point in pts: if point.position is None: print("跳过控制点 {}(没有坐标)".format(point.label)) continue item = calc_point_errors(chunk, point) if item is not None: results.append(item) return results def identify_outliers(results, sigma_threshold=3.0, abs_threshold=3.0): """用双阈值识别粗差点,返回可疑点列表和整体统计""" if len(results) == 0: return [], 0.0, 0.0 avg_all = sum(p["avg"] for p in results) / len(results) std_all = math.sqrt( sum((p["avg"] - avg_all) ** 2 for p in results) / len(results) ) outliers = [] for p in results: stat_flag = p["avg"] > avg_all + sigma_threshold * std_all abs_flag = p["avg"] > abs_threshold if stat_flag or abs_flag: outliers.append({ "label": p["label"], "avg": p["avg"], "std": p["std"], "max": p["max"], "stat_flag": stat_flag, "abs_flag": abs_flag, }) return outliers, avg_all, std_all def print_report(results, outliers, avg_all, std_all): """打印控制点粗差检测报告""" print("=" * 90) print("控制点重投影误差统计报告") print("=" * 90) header = "{:<24} {:>10} {:>10} {:>10} {:>6}".format( "点号", "平均误差(px)", "最大误差(px)", "标准差(px)", "影像数") print(header) print("-" * 90) for r in sorted(results, key=lambda x: x["avg"], reverse=True): print("{:<24} {:10.3f} {:10.3f} {:10.4f} {:6d}".format( r["label"], r["avg"], r["max"], r["std"], r["n"])) print("-" * 90) print("整体平均误差: {:.3f} px,整体标准差: {:.4f} px".format(avg_all, std_all)) if len(outliers) == 0: print("未发现可疑粗差控制点") else: print("可疑粗差控制点:") for o in outliers: reason = [] if o["stat_flag"]: reason.append("超过3sigma统计阈值") if o["abs_flag"]: reason.append("超过绝对阈值{:.1f}px".format(3.0)) print(" [{}] 平均误差 {:.3f}px,最大误差 {:.3f}px,原因:{}".format( o["label"], o["avg"], o["max"], ", ".join(reason))) print("=" * 90) def export_errors_to_csv(chunk, path): """把每个控制点在每张照片上的误差导出到CSV文件""" results = detect_all_errors(chunk) if len(results) == 0: print("没有可导出的误差数据") return with open(path, "w", newline="", encoding="utf-8-sig") as f: writer = csv.writer(f) writer.writerow(["点号", "影像名", "重投影误差(px)"]) for r in sorted(results, key=lambda x: x["label"]): for cam_label, err, _, _ in r["detail"]: writer.writerow([r["label"], cam_label, "{:.3f}".format(err)]) print("误差明细已导出到: {}".format(path)) def detail_errors(chunk, point_label): """查看指定控制点在每张照片上的投影误差(按误差从大到小排列)""" results = detect_all_errors(chunk) for r in results: if r["label"] == point_label: print("控制点 {} 在各影像上的投影误差:".format(point_label)) for cam_label, err, _, _ in sorted(r["detail"], key=lambda x: x[1], reverse=True): print(" {:<40} {:8.3f} px".format(cam_label[:40], err)) return print("未找到控制点 {}".format(point_label)) def main(): doc = Metashape.app.document if doc is None: print("没有打开的文档") return chunk = doc.chunk if chunk is None: print("当前chunk为空") return print("正在处理Chunk: {}".format(chunk.label)) results = detect_all_errors(chunk) if len(results) == 0: print("没有有效控制点(可能没有投影或相机未标定)") return outliers, avg_all, std_all = identify_outliers(results) print_report(results, outliers, avg_all, std_all) print("") print("提示:使用 detail_errors(chunk, '控制点号') 查看单点逐影像误差") print("提示:使用 export_errors_to_csv(chunk, 'C:/tmp/gcp_errors.csv') 导出明细") print("") # 在Metashape控制台运行时,自动执行main() if __name__ == "__main__": main()4.5 快速改造成批处理接口
上面这段代码在图形界面下直接运行就行。如果你希望把脚本集成到自己的批处理流程里,比如每天晚上自动跑一遍刚对齐完的测区,可以把main()函数里的前端控制台访问改成sandbox方式:用Metashape.app.param传入项目文件路径,脚本加载文档后对指定chunk做检测,最后把报告写入日志文件。基本改动只有两三行,核心的误差计算逻辑不用动。
批处理还有一个额外好处:可以不依赖图形界面,在一台配置较高的服务器上把所有测区的控制点粗差检测排队跑完。我在项目里就是这么用的,空三完成之后自动触发脚本,把可疑点列表作为告警信息推送出来,再去人工复核,效率提升很大。
5. 实测案例与输出解读
5.1 干净数据的输出
先看一个正常测区的效果。这个测区有12个控制点,影像分辨率约3.2cm/px,飞行高度150m,控制点全部由RTK实测。运行脚本后,输出报告的典型情况如下:
整体平均误差一般在0.3到0.8像素之间,标准差在0.1到0.3像素之间,最大的控制点平均误差也很少超过1.0像素。所有点的stat_flag和abs_flag都是False,报告输出"未发现可疑粗差控制点"。这种数据可以直接进入下一阶段的空三优化,不需要额外干预。
干净数据本身也能反映一些问题。比如你发现某个测区整体平均误差虽然不大,但标准差偏大,说明某些影像上的刺点精度不稳定,可能是有几张照片的对焦或清晰度不够。虽然不用剔除控制点,但可以提醒内业人员回头检查这些影像,减少潜在隐患。
5.2 含粗差数据的输出
同一个测区,如果其中一个控制点的平面坐标被错误地偏移了3米,运行脚本后,报告会明确显示该点的平均误差显著高于其他点。比如其他点平均0.6像素,该点平均达到5.7像素,而且该点会同时触发统计阈值和绝对阈值两个标记。
更有意思的情况是影像级粗差。有个控制点整体平均误差是1.9像素,看着好像还能接受,但通过detail_errors()逐影像查看,会发现该点在20张照片里的误差都稳定在0.4像素左右,唯独在某两张照片上误差达到8.2像素和7.6像素。这正是典型的刺点偏移,手动回到那两张照片就能很快修复,修复后重新优化,整体精度就能恢复。
5.3 如何根据报告处理粗差
拿到报告后,处理流程分三步。第一步,先用detail_errors()确认粗差类型——点级问题还是影像级问题。第二步,针对不同情况操作。如果是影像刺点错误,直接在Metashape的Ground Control面板中找到该控制点,双击打开照片列表,定位到误差偏大的照片,重新刺点。如果是坐标问题,建议直接在面板中禁用该控制点(取消enabled),而不是删除,因为删掉后再想恢复就需要重新输入坐标。第三步,剔除或修复可疑点后,重新执行一次"Align Chunks"或"Optimize Cameras",让空三重新收敛,再跑一遍脚本验证。
有一个重要提醒:粗差探测脚本最好在每次空三优化之后都跑一遍,而不是只在最终阶段跑。空三迭代过程中,如果粗差一直存在,每次优化都会被这个错误的约束拖拽,导致最终结果整体偏斜。早发现早处理,能省掉后续大量返工。
6. 踩坑记录与经验总结
6.1 常见报错排查
写脚本的过程中我踩过不少坑,挑几个最常见的分享出来,省得大家重复踩。
第一个坑是camera.project()返回None。这个函数不是对任何三维点都能成功投影,如果点位于相机后方或超出了有效投影范围,它会返回None。代码里我做了容错处理,遇到None直接跳过,不会中断主流程。如果发现某个控制点的有效投影数明显少于照片数,要留意是不是有照片的位姿精度特别差,导致点被投影到了影像范围之外。
第二个坑是控制点没有投影数据。用户在Ground Control面板里只输入了控制点坐标但没有刺点,或者刺点后没有执行对齐,point.projections就是空的。脚本会跳过这类点并打印提示,不用慌张。
第三个坑是chunk.crs.unproject()和chunk.transform.matrix的坐标系匹配问题。如果chunk的坐标参考在后期被修改过,比如从WGS-84换成了当地投影坐标系,之前存储的控制点位置可能无法被unproject()正确处理,函数会抛出异常。我在代码中对异常做了捕获并打印错误信息,方便定位。遇到这种情况,建议先在Metashape里重新设置chunk的CRS,或者用Transform Coordinates功能把控制点坐标更新到新坐标系下。
第四个坑是误用point.position直接投影。这个是新手最容易犯的错误。我第一次写这段代码的时候也直接用了camera.project(point.position),结果出来的误差值大得离谱,后来才意识到坐标系完全不对。建议在代码里加入按坐标转换后的投影坐标和实际刺点坐标的距离打印,快速比对验证。
下面用表格汇总一下常见的几个问题和处理方式:
| 现象 | 可能原因 | 排查方法 |
|---|---|---|
| 某个点所有照片误差都大 | 控制点坐标错误 | 检查外业记录、重新测量 |
| 某个点在个别照片误差大 | 影像刺点错误 | 用detail_errors()定位照片后重刺 |
| 所有点误差都偏大 | 相机位姿或云台参数问题 | 先检查空三质量,再判断控制点 |
camera.project()返回None | 投影点在相机后方或影像范围外 | 检查该照片位姿精度,考虑剔除 |
point.projections为空 | 控制点没有刺点 | 在Ground Control面板中完成刺点 |
| 坐标转换异常 | chunk的CRS被修改 | 重设CRS后再运行脚本 |
6.2 关于阈值的进一步思考
我写脚本的时候,把sigma_threshold和abs_threshold默认值分别设置为3.0和3.0,但这只是一个起点值,实际项目里几乎都需要调。
如果测区控制点数量特别少,比如只有6个控制点,3σ统计阈值的意义就很有限,因为样本量太小,标准差估计不稳定。这种情况下,优先相信绝对阈值。如果影像质量非常高,控制点分布均匀,误差通常在0.3像素左右,这时绝对阈值可以收紧到1.5像素,防止一个误差1.8像素的点混进去。
如果是复杂地形测区,比如山区带陡坎、高密度植被,刺点误差本身就比平地大,绝对阈值可以放宽一些,但要结合检查点(checkpoints)来验证。我通常的做法是:先按阈值跑一遍,剔除粗差,然后随机抽3到5个检查点评估平面和高程中误差,再根据结果微调阈值。这个循环比单纯相信统计公式可靠得多。
6.3 项目级自动化扩展方向
这套脚本除了当独立工具用,还可以扩展成项目级质检流程的一部分。比如在空三处理完后自动调用粗差探测脚本,把报告输出成CSV格式存档,同时在控制面板中自动禁用被标记为粗差的控制点。更进一步,还可以把控制点误差的变化趋势记录到数据库里,长期积累后用来评估外业测量水平、内业刺点质量。
我自己的经验是,粗差探测不能只做一次,而是应该在空三优化的每个迭代周期都跑一遍。尤其当测区范围大、控制点数量多的时候,第一轮粗差剔除后,重新优化可能又暴露出新的粗差点,因为之前粗差点的约束被去除后,整体平差结果发生了变化,其他点的误差分布也随之改变。这种情况下,第一次运行脚本可能只能识别出最明显的两个点,剔除后重新优化,再跑一次,又识别出两个点。不用担心这是脚本出错,这是粗差逐层暴露的正常过程。
7. 实测中的一句真心话
脚本写完之后,我在几个真实项目里做了验证,效果最明显的是那个有12个控制点的测区。外业布点时有一组坐标被RTK的杆高参数干扰了,偏差约2.8米,内业怎么优化都控制不住平面精度。用这套脚本跑了不到一分钟,直接把那个点钉在了报告首位。人工排查时,透过detail_errors()很快定位到问题点是坐标问题,回外业重测后重新优化,平面精度从差了7厘米压回到2.3厘米,整个项目顺利提交。
最后分享一个小技巧:如果你要在项目现场快速用这个脚本,建议把代码保存成.py文件,放到Metashape安装目录或者工作目录下,每次打开Metashape后直接用Tools > Run Script加载,不用复制粘贴。如果想要更快,还可以把它注册成一个Metashape自定义菜单项,鼠标点一下就弹出报告,操作体验会舒服很多。