尧图网站建设 
简介针对子孔径拼接与孔径拼接应用需求这份MLE程序包基于最大似然估计面向遥感成像、光学检测与SAR图像处理等场景解决多子孔径数据融合时几何配准和辐射一致性问题。压缩包共8个文件总大小约7.59MB核心为5个.m脚本包括主程序和圆形边界检测、外接圆半径计算等辅助函数另有2个.asv自动备份及1个.mat原始数据方便用户直接运行和修改验证。已有364人学习浏览适合具备图像处理与MATLAB基础的工程师和研究者。借助该包可掌握最大似然法完成子孔径拼接的完整流程从数据预处理、特征检测与匹配、位姿配准到辐射校正和整体融合同时通过附带的mat样例数据复现实验为合成孔径或光学拼接项目提供可扩展的代码参考。整体代码结构简洁模块划分清晰便于二次开发与替换实验数据。1. 干涉仪口径不够才轮到 MLE 子孔径拼接出场做大口径光学检测的同行都遇到过同一个瓶颈干涉仪的有效口径摆在那里被测件直径一旦超出常规全口径测量就做不了。预算不够买更大口径的干涉仪时工程里的常见出路就是把镜面拆成若干个重叠的子孔径分别测量再用算法拼回全口径面形这就是子孔径拼接也叫孔径拼接。可拼接不是简单地把图贴在一起每个子孔径自带平移、倾斜、离焦误差重叠区的相位残差往哪儿分配直接决定了拼接面形准不准。MLE_stitching 这类方案解决的核心问题就是在干涉图有明显噪声和相对对准误差的前提下用最大似然估计MLE同时求解全口径面形和每个子孔径的系统误差系数最终交付一张能用于验收的面形图。这篇笔记给做光学检测、干涉测量的工程师按建模、复现、参数、踩坑、验证的顺序写带你从拿到一堆子孔径数据到最后放心的拼接结果。2. 把每帧子孔径里的系统误差拆出来fact9eq 与最大似然建模2.1 子孔径里除了真实面形还有三类误差任何一卷子孔径数据单帧干涉图都能拆成三部分之和全口径面形在这个局部区域上的真实值、测量仪器叠加的系统误差、以及随机噪声。这里说清楚一个关键点第一项是所有子孔径共享的、真正要拼的对象第二项是每一帧各自独有的系统误差主要来自参考面本身的瑕疵、气流扰动、成像畸变和振动。它们在单帧内通常变化平缓可以用一组低阶基函数的线性组合去近似。在工程里这组基函数如果只取平移和倾斜两项拼出来的面形常会残留明显的低阶带状误差。原因很简单Fizeau 干涉仪的标准参考面并不是理想平面离焦、像散这类像差会以不同权重映射进每个子孔径。所以常见的拼接实现会把误差项扩充到 9 个系数命名里出现 fact9eq 这类标记时多数指的就是每个子孔径用 9 个低阶误差系数参与全局估计。选择 9 项而不是更多理由是现实的。前 9 个低阶像差项已经覆盖了活塞、倾斜、离焦、像散、彗差和三叶草的大部分实际误差参数再多一方面求解矩阵条件数变差另一方面系统误差项会开始吸收真实面形里的低频信息拼出来的结果反而不诚实。这个边界要记牢误差模型的核心目的是吸收测量系统带来的偏差不是把镜面本身的形状也拟合掉。2.2 fact9eq 的 9 项基函数怎么选Zernike 子集和多项式坐标一般平面干涉子孔径拼接里9 项配置的长这样第 1 项是活塞piston第 2、3 项是 x、y 方向倾斜第 4 项是离焦第 5、6 项是两个方向的像散第 7、8 项是两个方向的彗差第 9 项通常给球差或三叶草。选 Zernike 子集还是普通多项式取决于你代码里坐标系的写法。Zernike 在圆域上正交理论根底更干净但在拼接代码里子孔径是矩形像素阵列而且很多实现直接用 x、y 的幂次多项式理由是简单、不易搭配错归一化计算量也小。我建议工程上采用多项式并做归一化把每个子孔径的局部坐标压缩到 [-1, 1]避免 x 的量纲差异让矩阵的列之间数值差距过大。基函数必须在子孔径自身坐标系里定义不能用全口径全局坐标。这一点翻车率极高。一旦用了全局坐标系统误差中的倾斜项和全局面形在空间上的倾斜会混杂在一起系数不再可辨识拼接结果会产生虚假的大尺度倾斜。工程常用 9 项误差配置物理含义是否常被误读1活塞/平移容易和全局面形的均值自由度混淆2、3x/y 方向倾斜最容易吸收对准误差4离焦与子孔径沿光轴位置误差强耦合5、60°/45° 像散参考面应力误差的常见出口7、8x/y 彗差高阶对准残差通常幅值较小9球差或三叶草按残差形态二选一不可同时混入2.3 最大似然估计在拼接里做了什么先估计噪声方差再给像素配权重普通最小二乘拼接把重叠区每个像素的残差一视同仁地平方求和隐含假设是每个测量点的噪声方差一致。实际干涉图完全不是这样子孔径边缘解包裹质量差、有杂散光干扰的像素噪声方差可能比中心区大一个数量级。MLE 的做法是把求解过程变成两阶段交替先用一组初始权重做加权最小二乘得到初步面形然后从残差重新估计每个像素的噪声方差再以 1/σ² 为权重重新求解如此迭代直到方差估计收敛。权重低的地方残差对目标的拉动就弱坏点不会再把整个拼接带偏。这里常有一个误区拿 MLE 当普通加权最小二乘只加一次权就结束。最大似然的意义就在于那个权重不是预先猜出来的而是由残差的统计特性驱动更新的。实现上我会用1.4826 * median(abs(residual - median(residual)))这个稳健尺度去估 σ它比标准差更能扛离群点。这样拼出来的面形不再是被坏点拽着走的平均结果而是最可能产生这组观测数据的解。2.4 为什么拼接方程必须先固定基准自由度亏缺与三条出路把模型写成数学形式就会看到全局面形在每个有效像素上的值都是未知量自由度几乎等于全口径采样点数。如果不加额外约束方程有无限多解所有子孔径同时加同一个活塞全局面形相应整体平移重叠区残差一点不变。这就是拼接问题的基准自由度也叫零空间。要消除它工程里有三条常见出路。第一固定第一个子孔径的相对位置作为参考把它对应的全局像素值设为绝对参考。第二对全局面形的均值、x 和 y 倾斜施加弱正则约束把均值拉到零附近同时避免低阶像差被算法随意搬运。第三引入一个绝对测量值例如用另一台小口径干涉仪测某个子孔径的绝对面形把拼接结果锚定到真实面形上。前两种做法实现成本低绝大多数实验室代码选它们要提醒的是选哪种基准直接决定拼接输出的倾斜和离焦量面形图的绝对形状必须结合测量配置去理解算法保证的只是各子孔径之间相对拼接的一致性。3. 用 Python 跑通一次 MLE 子孔径拼接最小代码与必调参数3.1 输入约定相位图、掩模、位置表动手前先把数据格式统一。我一般准备三种输入每个子孔径的相位图是浮点数组单位可以是 nm 或 waves但全部子孔径必须同一单位每个子孔径的掩模是布尔数组标记有效像素每个子孔径的位置表包括它在全口径网格里的圆心坐标和半径。位置信息有的来自运动平台回零位有的靠干涉仪软件的坐标输出。掩模必须提前做边缘收缩用腐蚀操作把最外沿 3 到 5 个像素剔除这比任何算法参数都管用。边缘像素的解包裹跳变和杂散光干扰最严重留着它们拼接残差会明显变大边缘翘起的概率也大幅上升。位置表还有一个常见坑子孔径图的分辨率必须一致像素间距不统一时要先插值到共同网格否则重叠区的对应关系是错的后续配准再怎么调也是白费。我一般写一个预处理函数把相位图、掩模和位置表打包成统一的 list 结构再进拼接主体。3.2 拼接主体代码稀疏观测矩阵与迭代加权下面这个最小实现展示了 MLE 拼接的核心流程。代码省略了部分工程细节但保留了最关键的稀疏矩阵组装和迭代加权逻辑可以直接在小规模数据上跑通。import numpy as np from scipy.sparse import coo_matrix from scipy.sparse.linalg import lsqr def build_basis(x, y): 在子孔径局部坐标上构造 9 项误差基函数 xs x / (x.max() - x.min() 1e-12) ys y / (y.max() - y.min() 1e-12) # 9 列: 1, x, y, x^2, xy, y^2, x^3, y^3, x^2*y return np.column_stack([ np.ones_like(xs), xs, ys, xs**2, xs*ys, ys**2, xs**3, ys**3, xs**2 * ys ]) def stitch_mle(phases, masks, centers, radius, damp1e-3, max_iter30): phases/masks: 按子孔径顺序排列的数组列表 rows, cols, vals, rhs [], [], [], [] n_seg len(phases) row 0 # 先建立全局网格坐标与索引这里假设所有子孔径在同一网格坐标系下 ny, nx phases[0].shape global_idx -np.ones((ny, nx), dtypenp.int64) valid_count 0 for y in range(ny): for x in range(nx): # 凡是被任一掩模覆盖的像素都进入全局未知量 if any(mask[y, x] for mask in masks): global_idx[y, x] valid_count valid_count 1 for i in range(n_seg): mask masks[i] ys, xs np.where(mask) # 子孔径有效像素坐标 phase_i phases[i][mask] basis_i build_basis(xs, ys) # (有效像素数, 9) # 全局面形未知量放在前面9 项系数放在后面 for idx_pt in range(len(xs)): g global_idx[ys[idx_pt], xs[idx_pt]] rows.append(row) cols.append(g) vals.append(1.0) # 全局面形贡献 for j in range(9): rows.append(row) cols.append(valid_count i * 9 j) vals.append(basis_i[idx_pt, j]) # 误差基函数贡献 rhs.append(phase_i[idx_pt]) row 1 A coo_matrix((vals, (rows, cols))).tocsr() b np.array(rhs) # 第一轮普通最小二乘拿初始残差 x lsqr(A, b, dampdamp, iter_lim200, atol1e-6, btol1e-6)[0] residual b - A x for it in range(max_iter): # 稳健估计噪声标准差权重 1 / sigma^2 sigma 1.4826 * np.median(np.abs(residual - np.median(residual))) weight np.sqrt(1.0 / (sigma**2 1e-12)) # 用行权重重新求解 x lsqr(A, b, dampdamp, sqrt_waweight, iter_lim200, atol1e-6, btol1e-6)[0] new_residual b - A x delta np.abs(np.median(np.abs(new_residual)) - np.median(np.abs(residual))) residual new_residual if delta 1e-8: break # 前 valid_count 个未知量是全域面形后面是各子孔径误差系数 surface np.full((ny, nx), np.nan) surface[global_idx 0] x[:valid_count] coefs x[valid_count:].reshape(n_seg, 9) return surface, coefs, residual这段代码的求解逻辑分两条线观测矩阵把每个像素的相位拆成全局面形值和当前子孔径的 9 项误差系数然后用lsqr解这个稀疏线性系统。lsqr的行权重参数sqrt_wa直接实现了 MLE 里按噪声方差加权这一步不需要手动把整个矩阵重新乘一遍。迭代只在噪声尺度变化超过阈值时继续通常 5 到 10 轮就收敛。代码里的damp参数是给方程加一个小正则防止矩阵病态时解出离谱的大数值。要注意的是上面为了可读性用了 Python 循环去组装稀疏矩阵数据量大的时候会慢。实际工程里我会改成向量化构造coo_matrix的三元组但逻辑完全一致。另外这里假设所有子孔径在同一网格分辨率下若像素间距不一致务必先在预处理里统一。3.3 三个必调参数重叠率、基函数阶次、正则强度重叠率是拼接刚度的直接来源。子孔径之间没有重叠的像素方程就是独立的拼不起来重叠太少又会让重叠区里的噪声被过度放大。我的经验值是重叠宽度占子孔径口径的 20% 到 35%最少不要低于 15%。重叠率越高约束越强矩阵越病态计算时间也越长低于 15% 时拼接结果对位置误差变得异常敏感属于左右为难的区域。基函数阶次决定误差模型吸收能力。先跑一版 9 项看重叠区残差图里是否残留明显的条纹或带状图案。如果残差里依然有低阶像差形态试着换成 14 项或更多如果拼出的全口径面形出现之前单帧里看不出来的波浪形多半是阶次太高把真实面形吃进去了。判断准则只有一个残差的 RMS 应该接近单帧测量的重复性而不是一味降到零。正则强度damp的量级要和相位数据的数值单位挂钩。相位值以波为单位的damp 取 1e-3 到 1e-2 之间以纳米为单位时往往要放到 1e-5 到 1e-3。一个可复用的调试办法先设成 0 跑一次观察解的振荡幅度再把 damp 从 1e-6 起按 10 倍递增直到解的 PV 值不再明显下降。damp 太大也会压制真实面形所以记住它只用来约束零空间不是用来平滑面形的。3.4 输出检查重叠区残差 RMS 和收敛曲线跑完代码别急着看面形图漂亮就收工。第一件事是提取所有子孔径重叠区里的残差统计 RMS。理想情况下重叠区残差 RMS 应该在单帧测量重复性的两倍以内若明显偏大说明基函数配置、位置配准或掩模腐蚀有环节出问题。第二件事是看收敛曲线每轮迭代的残差中位数应该单调下降并趋于平缓如果上下跳动说明有像素在干扰权重更新回到稳健估计算子去排查。第三件事是把全局面形里的 NaN 区域画出来确认边界形状符合预期未被掩模覆盖的区域不要参与任何后续分析。4. 子孔径拼接避坑清单条纹、翘边、不收敛的排查顺序4.1 现象拼接面形上出现周期性的水波纹条纹拼接结果里出现类似等高线的规则条纹尤其沿子孔径排列方向周期性重复多半是误差模型缺项。比如只用 5 项系数没有给球差预留位置参考面的球差就会在每个子孔径里以相同的方向叠加上去在全口径上呈现规律的波浪带。另一个常见原因是权重函数在子孔径边缘给得太足边缘噪声被当成真实信号参与约束。遇到条纹先做两件事打开残差图看条纹是集中在重叠区还是全口径均匀分布把误差模型从 9 项扩到 14 项对比残差 RMS 的变化幅度。如果残差显著下降且有规律的条纹消失就是缺项问题如果残差没变则是权重分配问题把边缘像素的权重压低或者在重叠区改用中间带加权。4.2 现象边缘翘起面形四周明显上扬或下沉边缘翘起是子孔径拼接最典型的表现原因很直接掩模没腐蚀干净最外圈的坏像素进入了重叠区。这些位置相位解包裹最容易跳变残差巨大偏偏又处在重叠约束的边缘算法为了降低这些点的残差会人为抬高或压低整个子孔径的倾斜项最后反映成全口径边缘的翘起。另一个次因是子孔径沿光轴方向位置不一致离焦项把这个位置差吸收掉了也会在边缘留下中频翘曲。解决办法优先检查掩模。对掩模做半径为 3 到 5 像素的腐蚀把边缘无效区剔除再看重叠区像素是哪个子孔径贡献的如果有某帧的边缘像素频繁成为残差离群点直接在该帧上做局部掩模裁剪。做完掩模处理仍然翘去检查离焦项系数是否随子孔径位置呈系统变化若是考虑在误差模型里固定离焦项或引入测距信息限制它的幅值。4.3 现象迭代不收敛每一轮残差上下跳动MLE 迭代的本质是用残差估计权重再用权重改解残差又变权重再改。当数据里混入少量极端离群点时稳健尺度估计会失效导致下一轮权重把坏点放大于是残差中枢一路抖动。出现这种情况先在每个子孔径内做一次 3σ 截断把明显超出中位数邻域的相位点标记为无效再重新进入迭代循环。这一步看似简单却能干掉大半不收敛问题。若截断后仍然抖动去查基函数之间的线性相关性。x、y 的多项式在高阶项里有天然的相关性尤其当有效像素区域不是完整圆形时某些高阶项会和低阶项高度共线矩阵条件数恶化稀疏求解器在数值误差里来回震荡。处理办法有两个一是改用 Zernike 子集并做正交化预处理二是提高 damp 正则值同时观察残差的变化。我通常先用正则压住再用残差确认没有损失真实信号最后才考虑换基函数。4.4 现象结果对初始重叠位置极其敏感位置偏一个像素PV 变一倍这几乎可以锁定是子孔径之间的相对位置配准不准或者是重叠率太低。位置误差直接改写观测矩阵里全局像素与子孔径像素的对应关系但不改基函数所以误差会以倾斜和离焦的形式被吸收面形 PV 值跟着剧烈变化。干涉仪运动平台的定位误差一般在微米级但相机像素尺寸可能只有几微米平台回零位误差变成亚像素到几个像素都很常见。解决思路是在进入 MLE 拼接前先用互相关做一次细配准。取相邻子孔径重叠区域的相位图计算两者的二维互相关找到相关峰偏移量修正位置表。这一步做完再进拼接绝大多数位置敏感问题都会消失。如果修正后仍然敏感检查重叠率是不是低于 15%小于这个值时就回退采集方案增加步进密度。4.5 现象拼接出的全口径面形和整镜直接测量对不上这是最难排查的一类问题。拼接算法的先天弱点是绝对低阶面形特别是离焦和球差与子孔径沿光轴方向的位置误差强耦合。你移动干涉仪或被测镜时微小的轴向位移会带来离焦误差MLE 只能把它吸收进每个子孔径的离焦系数无法区分它到底是镜面真实形状还是位置误差导致的伪像。所以拼接面形的中高频成分可信但绝对离焦、绝对球差的数值未必可靠。没有后悔药只能靠标定。定期用标准球面或标准平面做一次绝对标定把轴向位置误差的尺度标出来或者对同一块镜面用旋转拼接、平移拼接各测一遍比对低阶项差异。如果差异超出你预期的绝对精度老老实实引入辅助测距把每个子孔径的轴向位置记录下来再作为先验约束放进误差模型。这条经验几乎成了我的血泪教训曾有一块镜面拼接结果球差异常排了三天最后是样品架上一颗螺丝的端面倾斜导致的。5. 交付前最后一步用重叠区重拼来验证拼接结果拼接算法跑通了面形图也画出来了还差一步证明结果可信。我一般会做一个重叠区独立重拼测试把参与拼接的原始数据随机抽掉一个子孔径用剩余的全部子孔径重新拼接然后拿抽掉的子孔径真实测量值和重拼结果在该区域对比。两者的差值如果与测量重复性相当说明拼接没有为了迁就某帧坏数据而扭曲全局面形。这个验证方法的成本几乎为零只消耗一点计算时间却能把绝大多数误拼抓出来。验证手段做法通过判据单帧剔除重拼随机剔除 1 个子孔径后重拼与原始拼接对比重叠区残差 RMS 漂移小于 20%重叠率扫描用 20%、30%、40% 三组参数重拼面形 PV 变化不超过 RMS 的三倍旋转一致性整镜旋转 90° 后重拼低阶项旋转对称分量一致绝对标定比较用标准面或辅助测长仪校准离焦/球差低阶项差值与标定不确定度匹配另一种值得做的验证是对称配准检查。把子孔径顺序从正序改成逆序也就是反向进场重新测量一组数据再进同一套流程拼接。顺序反了同一个系统误差在正反两次里的符号会变如果拼接结果和原始结果在低阶项上基本一致说明误差模型有可重复性如果符号跟着测序走说明某些误差项没有被模型吸收干净而是随着测量环境在变。这套操作做下来才能放心把拼接面形写进检测报告。最后说一个习惯每次拼接完我会把收敛曲线、重叠区残差图和权重分布快照一并归档而不只是存一张看起来漂亮的面形图。两个月后别人问起这份数据的拼接质量你翻出这些中间产物就能解释一切。子孔径拼接是个看起来简单、实际很吃细节的活9 项系数能解决大部分常规问题真正决定交付质量的永远是掩模、配准和验证这三件小事。希望帮到你。本文还有配套的精品资源点击获取 干涉测量中的MLE子孔径拼接：原理、Python实现与工程避坑
简介针对子孔径拼接与孔径拼接应用需求这份MLE程序包基于最大似然估计面向遥感成像、光学检测与SAR图像处理等场景解决多子孔径数据融合时几何配准和辐射一致性问题。压缩包共8个文件总大小约7.59MB核心为5个.m脚本包括主程序和圆形边界检测、外接圆半径计算等辅助函数另有2个.asv自动备份及1个.mat原始数据方便用户直接运行和修改验证。已有364人学习浏览适合具备图像处理与MATLAB基础的工程师和研究者。借助该包可掌握最大似然法完成子孔径拼接的完整流程从数据预处理、特征检测与匹配、位姿配准到辐射校正和整体融合同时通过附带的mat样例数据复现实验为合成孔径或光学拼接项目提供可扩展的代码参考。整体代码结构简洁模块划分清晰便于二次开发与替换实验数据。1. 干涉仪口径不够才轮到 MLE 子孔径拼接出场做大口径光学检测的同行都遇到过同一个瓶颈干涉仪的有效口径摆在那里被测件直径一旦超出常规全口径测量就做不了。预算不够买更大口径的干涉仪时工程里的常见出路就是把镜面拆成若干个重叠的子孔径分别测量再用算法拼回全口径面形这就是子孔径拼接也叫孔径拼接。可拼接不是简单地把图贴在一起每个子孔径自带平移、倾斜、离焦误差重叠区的相位残差往哪儿分配直接决定了拼接面形准不准。MLE_stitching 这类方案解决的核心问题就是在干涉图有明显噪声和相对对准误差的前提下用最大似然估计MLE同时求解全口径面形和每个子孔径的系统误差系数最终交付一张能用于验收的面形图。这篇笔记给做光学检测、干涉测量的工程师按建模、复现、参数、踩坑、验证的顺序写带你从拿到一堆子孔径数据到最后放心的拼接结果。2. 把每帧子孔径里的系统误差拆出来fact9eq 与最大似然建模2.1 子孔径里除了真实面形还有三类误差任何一卷子孔径数据单帧干涉图都能拆成三部分之和全口径面形在这个局部区域上的真实值、测量仪器叠加的系统误差、以及随机噪声。这里说清楚一个关键点第一项是所有子孔径共享的、真正要拼的对象第二项是每一帧各自独有的系统误差主要来自参考面本身的瑕疵、气流扰动、成像畸变和振动。它们在单帧内通常变化平缓可以用一组低阶基函数的线性组合去近似。在工程里这组基函数如果只取平移和倾斜两项拼出来的面形常会残留明显的低阶带状误差。原因很简单Fizeau 干涉仪的标准参考面并不是理想平面离焦、像散这类像差会以不同权重映射进每个子孔径。所以常见的拼接实现会把误差项扩充到 9 个系数命名里出现 fact9eq 这类标记时多数指的就是每个子孔径用 9 个低阶误差系数参与全局估计。选择 9 项而不是更多理由是现实的。前 9 个低阶像差项已经覆盖了活塞、倾斜、离焦、像散、彗差和三叶草的大部分实际误差参数再多一方面求解矩阵条件数变差另一方面系统误差项会开始吸收真实面形里的低频信息拼出来的结果反而不诚实。这个边界要记牢误差模型的核心目的是吸收测量系统带来的偏差不是把镜面本身的形状也拟合掉。2.2 fact9eq 的 9 项基函数怎么选Zernike 子集和多项式坐标一般平面干涉子孔径拼接里9 项配置的长这样第 1 项是活塞piston第 2、3 项是 x、y 方向倾斜第 4 项是离焦第 5、6 项是两个方向的像散第 7、8 项是两个方向的彗差第 9 项通常给球差或三叶草。选 Zernike 子集还是普通多项式取决于你代码里坐标系的写法。Zernike 在圆域上正交理论根底更干净但在拼接代码里子孔径是矩形像素阵列而且很多实现直接用 x、y 的幂次多项式理由是简单、不易搭配错归一化计算量也小。我建议工程上采用多项式并做归一化把每个子孔径的局部坐标压缩到 [-1, 1]避免 x 的量纲差异让矩阵的列之间数值差距过大。基函数必须在子孔径自身坐标系里定义不能用全口径全局坐标。这一点翻车率极高。一旦用了全局坐标系统误差中的倾斜项和全局面形在空间上的倾斜会混杂在一起系数不再可辨识拼接结果会产生虚假的大尺度倾斜。工程常用 9 项误差配置物理含义是否常被误读1活塞/平移容易和全局面形的均值自由度混淆2、3x/y 方向倾斜最容易吸收对准误差4离焦与子孔径沿光轴位置误差强耦合5、60°/45° 像散参考面应力误差的常见出口7、8x/y 彗差高阶对准残差通常幅值较小9球差或三叶草按残差形态二选一不可同时混入2.3 最大似然估计在拼接里做了什么先估计噪声方差再给像素配权重普通最小二乘拼接把重叠区每个像素的残差一视同仁地平方求和隐含假设是每个测量点的噪声方差一致。实际干涉图完全不是这样子孔径边缘解包裹质量差、有杂散光干扰的像素噪声方差可能比中心区大一个数量级。MLE 的做法是把求解过程变成两阶段交替先用一组初始权重做加权最小二乘得到初步面形然后从残差重新估计每个像素的噪声方差再以 1/σ² 为权重重新求解如此迭代直到方差估计收敛。权重低的地方残差对目标的拉动就弱坏点不会再把整个拼接带偏。这里常有一个误区拿 MLE 当普通加权最小二乘只加一次权就结束。最大似然的意义就在于那个权重不是预先猜出来的而是由残差的统计特性驱动更新的。实现上我会用1.4826 * median(abs(residual - median(residual)))这个稳健尺度去估 σ它比标准差更能扛离群点。这样拼出来的面形不再是被坏点拽着走的平均结果而是最可能产生这组观测数据的解。2.4 为什么拼接方程必须先固定基准自由度亏缺与三条出路把模型写成数学形式就会看到全局面形在每个有效像素上的值都是未知量自由度几乎等于全口径采样点数。如果不加额外约束方程有无限多解所有子孔径同时加同一个活塞全局面形相应整体平移重叠区残差一点不变。这就是拼接问题的基准自由度也叫零空间。要消除它工程里有三条常见出路。第一固定第一个子孔径的相对位置作为参考把它对应的全局像素值设为绝对参考。第二对全局面形的均值、x 和 y 倾斜施加弱正则约束把均值拉到零附近同时避免低阶像差被算法随意搬运。第三引入一个绝对测量值例如用另一台小口径干涉仪测某个子孔径的绝对面形把拼接结果锚定到真实面形上。前两种做法实现成本低绝大多数实验室代码选它们要提醒的是选哪种基准直接决定拼接输出的倾斜和离焦量面形图的绝对形状必须结合测量配置去理解算法保证的只是各子孔径之间相对拼接的一致性。3. 用 Python 跑通一次 MLE 子孔径拼接最小代码与必调参数3.1 输入约定相位图、掩模、位置表动手前先把数据格式统一。我一般准备三种输入每个子孔径的相位图是浮点数组单位可以是 nm 或 waves但全部子孔径必须同一单位每个子孔径的掩模是布尔数组标记有效像素每个子孔径的位置表包括它在全口径网格里的圆心坐标和半径。位置信息有的来自运动平台回零位有的靠干涉仪软件的坐标输出。掩模必须提前做边缘收缩用腐蚀操作把最外沿 3 到 5 个像素剔除这比任何算法参数都管用。边缘像素的解包裹跳变和杂散光干扰最严重留着它们拼接残差会明显变大边缘翘起的概率也大幅上升。位置表还有一个常见坑子孔径图的分辨率必须一致像素间距不统一时要先插值到共同网格否则重叠区的对应关系是错的后续配准再怎么调也是白费。我一般写一个预处理函数把相位图、掩模和位置表打包成统一的 list 结构再进拼接主体。3.2 拼接主体代码稀疏观测矩阵与迭代加权下面这个最小实现展示了 MLE 拼接的核心流程。代码省略了部分工程细节但保留了最关键的稀疏矩阵组装和迭代加权逻辑可以直接在小规模数据上跑通。import numpy as np from scipy.sparse import coo_matrix from scipy.sparse.linalg import lsqr def build_basis(x, y): 在子孔径局部坐标上构造 9 项误差基函数 xs x / (x.max() - x.min() 1e-12) ys y / (y.max() - y.min() 1e-12) # 9 列: 1, x, y, x^2, xy, y^2, x^3, y^3, x^2*y return np.column_stack([ np.ones_like(xs), xs, ys, xs**2, xs*ys, ys**2, xs**3, ys**3, xs**2 * ys ]) def stitch_mle(phases, masks, centers, radius, damp1e-3, max_iter30): phases/masks: 按子孔径顺序排列的数组列表 rows, cols, vals, rhs [], [], [], [] n_seg len(phases) row 0 # 先建立全局网格坐标与索引这里假设所有子孔径在同一网格坐标系下 ny, nx phases[0].shape global_idx -np.ones((ny, nx), dtypenp.int64) valid_count 0 for y in range(ny): for x in range(nx): # 凡是被任一掩模覆盖的像素都进入全局未知量 if any(mask[y, x] for mask in masks): global_idx[y, x] valid_count valid_count 1 for i in range(n_seg): mask masks[i] ys, xs np.where(mask) # 子孔径有效像素坐标 phase_i phases[i][mask] basis_i build_basis(xs, ys) # (有效像素数, 9) # 全局面形未知量放在前面9 项系数放在后面 for idx_pt in range(len(xs)): g global_idx[ys[idx_pt], xs[idx_pt]] rows.append(row) cols.append(g) vals.append(1.0) # 全局面形贡献 for j in range(9): rows.append(row) cols.append(valid_count i * 9 j) vals.append(basis_i[idx_pt, j]) # 误差基函数贡献 rhs.append(phase_i[idx_pt]) row 1 A coo_matrix((vals, (rows, cols))).tocsr() b np.array(rhs) # 第一轮普通最小二乘拿初始残差 x lsqr(A, b, dampdamp, iter_lim200, atol1e-6, btol1e-6)[0] residual b - A x for it in range(max_iter): # 稳健估计噪声标准差权重 1 / sigma^2 sigma 1.4826 * np.median(np.abs(residual - np.median(residual))) weight np.sqrt(1.0 / (sigma**2 1e-12)) # 用行权重重新求解 x lsqr(A, b, dampdamp, sqrt_waweight, iter_lim200, atol1e-6, btol1e-6)[0] new_residual b - A x delta np.abs(np.median(np.abs(new_residual)) - np.median(np.abs(residual))) residual new_residual if delta 1e-8: break # 前 valid_count 个未知量是全域面形后面是各子孔径误差系数 surface np.full((ny, nx), np.nan) surface[global_idx 0] x[:valid_count] coefs x[valid_count:].reshape(n_seg, 9) return surface, coefs, residual这段代码的求解逻辑分两条线观测矩阵把每个像素的相位拆成全局面形值和当前子孔径的 9 项误差系数然后用lsqr解这个稀疏线性系统。lsqr的行权重参数sqrt_wa直接实现了 MLE 里按噪声方差加权这一步不需要手动把整个矩阵重新乘一遍。迭代只在噪声尺度变化超过阈值时继续通常 5 到 10 轮就收敛。代码里的damp参数是给方程加一个小正则防止矩阵病态时解出离谱的大数值。要注意的是上面为了可读性用了 Python 循环去组装稀疏矩阵数据量大的时候会慢。实际工程里我会改成向量化构造coo_matrix的三元组但逻辑完全一致。另外这里假设所有子孔径在同一网格分辨率下若像素间距不一致务必先在预处理里统一。3.3 三个必调参数重叠率、基函数阶次、正则强度重叠率是拼接刚度的直接来源。子孔径之间没有重叠的像素方程就是独立的拼不起来重叠太少又会让重叠区里的噪声被过度放大。我的经验值是重叠宽度占子孔径口径的 20% 到 35%最少不要低于 15%。重叠率越高约束越强矩阵越病态计算时间也越长低于 15% 时拼接结果对位置误差变得异常敏感属于左右为难的区域。基函数阶次决定误差模型吸收能力。先跑一版 9 项看重叠区残差图里是否残留明显的条纹或带状图案。如果残差里依然有低阶像差形态试着换成 14 项或更多如果拼出的全口径面形出现之前单帧里看不出来的波浪形多半是阶次太高把真实面形吃进去了。判断准则只有一个残差的 RMS 应该接近单帧测量的重复性而不是一味降到零。正则强度damp的量级要和相位数据的数值单位挂钩。相位值以波为单位的damp 取 1e-3 到 1e-2 之间以纳米为单位时往往要放到 1e-5 到 1e-3。一个可复用的调试办法先设成 0 跑一次观察解的振荡幅度再把 damp 从 1e-6 起按 10 倍递增直到解的 PV 值不再明显下降。damp 太大也会压制真实面形所以记住它只用来约束零空间不是用来平滑面形的。3.4 输出检查重叠区残差 RMS 和收敛曲线跑完代码别急着看面形图漂亮就收工。第一件事是提取所有子孔径重叠区里的残差统计 RMS。理想情况下重叠区残差 RMS 应该在单帧测量重复性的两倍以内若明显偏大说明基函数配置、位置配准或掩模腐蚀有环节出问题。第二件事是看收敛曲线每轮迭代的残差中位数应该单调下降并趋于平缓如果上下跳动说明有像素在干扰权重更新回到稳健估计算子去排查。第三件事是把全局面形里的 NaN 区域画出来确认边界形状符合预期未被掩模覆盖的区域不要参与任何后续分析。4. 子孔径拼接避坑清单条纹、翘边、不收敛的排查顺序4.1 现象拼接面形上出现周期性的水波纹条纹拼接结果里出现类似等高线的规则条纹尤其沿子孔径排列方向周期性重复多半是误差模型缺项。比如只用 5 项系数没有给球差预留位置参考面的球差就会在每个子孔径里以相同的方向叠加上去在全口径上呈现规律的波浪带。另一个常见原因是权重函数在子孔径边缘给得太足边缘噪声被当成真实信号参与约束。遇到条纹先做两件事打开残差图看条纹是集中在重叠区还是全口径均匀分布把误差模型从 9 项扩到 14 项对比残差 RMS 的变化幅度。如果残差显著下降且有规律的条纹消失就是缺项问题如果残差没变则是权重分配问题把边缘像素的权重压低或者在重叠区改用中间带加权。4.2 现象边缘翘起面形四周明显上扬或下沉边缘翘起是子孔径拼接最典型的表现原因很直接掩模没腐蚀干净最外圈的坏像素进入了重叠区。这些位置相位解包裹最容易跳变残差巨大偏偏又处在重叠约束的边缘算法为了降低这些点的残差会人为抬高或压低整个子孔径的倾斜项最后反映成全口径边缘的翘起。另一个次因是子孔径沿光轴方向位置不一致离焦项把这个位置差吸收掉了也会在边缘留下中频翘曲。解决办法优先检查掩模。对掩模做半径为 3 到 5 像素的腐蚀把边缘无效区剔除再看重叠区像素是哪个子孔径贡献的如果有某帧的边缘像素频繁成为残差离群点直接在该帧上做局部掩模裁剪。做完掩模处理仍然翘去检查离焦项系数是否随子孔径位置呈系统变化若是考虑在误差模型里固定离焦项或引入测距信息限制它的幅值。4.3 现象迭代不收敛每一轮残差上下跳动MLE 迭代的本质是用残差估计权重再用权重改解残差又变权重再改。当数据里混入少量极端离群点时稳健尺度估计会失效导致下一轮权重把坏点放大于是残差中枢一路抖动。出现这种情况先在每个子孔径内做一次 3σ 截断把明显超出中位数邻域的相位点标记为无效再重新进入迭代循环。这一步看似简单却能干掉大半不收敛问题。若截断后仍然抖动去查基函数之间的线性相关性。x、y 的多项式在高阶项里有天然的相关性尤其当有效像素区域不是完整圆形时某些高阶项会和低阶项高度共线矩阵条件数恶化稀疏求解器在数值误差里来回震荡。处理办法有两个一是改用 Zernike 子集并做正交化预处理二是提高 damp 正则值同时观察残差的变化。我通常先用正则压住再用残差确认没有损失真实信号最后才考虑换基函数。4.4 现象结果对初始重叠位置极其敏感位置偏一个像素PV 变一倍这几乎可以锁定是子孔径之间的相对位置配准不准或者是重叠率太低。位置误差直接改写观测矩阵里全局像素与子孔径像素的对应关系但不改基函数所以误差会以倾斜和离焦的形式被吸收面形 PV 值跟着剧烈变化。干涉仪运动平台的定位误差一般在微米级但相机像素尺寸可能只有几微米平台回零位误差变成亚像素到几个像素都很常见。解决思路是在进入 MLE 拼接前先用互相关做一次细配准。取相邻子孔径重叠区域的相位图计算两者的二维互相关找到相关峰偏移量修正位置表。这一步做完再进拼接绝大多数位置敏感问题都会消失。如果修正后仍然敏感检查重叠率是不是低于 15%小于这个值时就回退采集方案增加步进密度。4.5 现象拼接出的全口径面形和整镜直接测量对不上这是最难排查的一类问题。拼接算法的先天弱点是绝对低阶面形特别是离焦和球差与子孔径沿光轴方向的位置误差强耦合。你移动干涉仪或被测镜时微小的轴向位移会带来离焦误差MLE 只能把它吸收进每个子孔径的离焦系数无法区分它到底是镜面真实形状还是位置误差导致的伪像。所以拼接面形的中高频成分可信但绝对离焦、绝对球差的数值未必可靠。没有后悔药只能靠标定。定期用标准球面或标准平面做一次绝对标定把轴向位置误差的尺度标出来或者对同一块镜面用旋转拼接、平移拼接各测一遍比对低阶项差异。如果差异超出你预期的绝对精度老老实实引入辅助测距把每个子孔径的轴向位置记录下来再作为先验约束放进误差模型。这条经验几乎成了我的血泪教训曾有一块镜面拼接结果球差异常排了三天最后是样品架上一颗螺丝的端面倾斜导致的。5. 交付前最后一步用重叠区重拼来验证拼接结果拼接算法跑通了面形图也画出来了还差一步证明结果可信。我一般会做一个重叠区独立重拼测试把参与拼接的原始数据随机抽掉一个子孔径用剩余的全部子孔径重新拼接然后拿抽掉的子孔径真实测量值和重拼结果在该区域对比。两者的差值如果与测量重复性相当说明拼接没有为了迁就某帧坏数据而扭曲全局面形。这个验证方法的成本几乎为零只消耗一点计算时间却能把绝大多数误拼抓出来。验证手段做法通过判据单帧剔除重拼随机剔除 1 个子孔径后重拼与原始拼接对比重叠区残差 RMS 漂移小于 20%重叠率扫描用 20%、30%、40% 三组参数重拼面形 PV 变化不超过 RMS 的三倍旋转一致性整镜旋转 90° 后重拼低阶项旋转对称分量一致绝对标定比较用标准面或辅助测长仪校准离焦/球差低阶项差值与标定不确定度匹配另一种值得做的验证是对称配准检查。把子孔径顺序从正序改成逆序也就是反向进场重新测量一组数据再进同一套流程拼接。顺序反了同一个系统误差在正反两次里的符号会变如果拼接结果和原始结果在低阶项上基本一致说明误差模型有可重复性如果符号跟着测序走说明某些误差项没有被模型吸收干净而是随着测量环境在变。这套操作做下来才能放心把拼接面形写进检测报告。最后说一个习惯每次拼接完我会把收敛曲线、重叠区残差图和权重分布快照一并归档而不只是存一张看起来漂亮的面形图。两个月后别人问起这份数据的拼接质量你翻出这些中间产物就能解释一切。子孔径拼接是个看起来简单、实际很吃细节的活9 项系数能解决大部分常规问题真正决定交付质量的永远是掩模、配准和验证这三件小事。希望帮到你。本文还有配套的精品资源点击获取