MATLAB哈特曼标定仿真：从质心提取到Zernike波前重建

简介：面向哈特曼波前传感器标定与模式复原需求的MATLAB仿真代码包，适合自适应光学、波前检测方向的研究生、工程师以及相关课程设计人员。资源以泽尼克多项式为核心，完整覆盖哈特曼标定流程、zernike模式复原过程和DM响应函数广义逆的计算，从光斑质心提取到响应矩阵求解均有对应脚本支撑，方便对照理论自行修改或迁移到实际项目中。压缩包共18个文件，以11个m脚本为主，承担标定、质心计算、FFT处理、圆形区域分析等任务；4个mat数据文件用于保存泽尼克系数等初始化参数，另有fig结果图示与m~备份文件，整体仅2.91MB，下载与离线查阅都很方便。脚本细节还包括平面标定、曲面显示与响应矩阵构建等方面，读者可以从数据准备、算法实现到结果验证串起完整仿真链路，适合作为毕设或课题验证的起点。目前已有864人学习下载，对入门哈特曼标定与变形镜响应分析而言是一份实用参考。

1. 哈特曼标定的核心矛盾：像素怎么换算成波前斜率

哈特曼标定这四个字，看着像光学实验里的一个杂活，实际上它是夏克-哈特曼波前传感器从“能拍图”变成“能测波前”的那道门槛。很多人第一次跑通质心提取就以为完事了，结果把斜率重建出来一看，像散和离焦全混在一起，RMS 大得离谱。问题几乎都出在标定上：你并没有真正知道每个微透镜在 CCD 上的准确位置，也没有把“像素偏移”和“波前斜率”之间的比例系数定准。MATLAB 仿真做这件事最值钱的地方，是能让你在没买微透镜阵列和相机之前，就把整套标定链路跑一遍，而且可以随时注入已知像差当“标准答案”。这篇文章适合要做自适应光学、准备搭 SHWFS 硬件、或者只想把波前测量模块封装进系统里的人；读完之后你能自己搭一套可复现的仿真、知道每个参数怎么调、也知道哪几步最容易翻车。

2. 在 MATLAB 里搭一套可调参的哈特曼仿真链路：从波前到焦斑图

2.1 仿真前定四个参数：波长、微透镜间距、焦距与 CCD 像素尺寸

哈特曼标定仿真避免不了先把传感器参数定下来，因为这些参数直接决定焦斑在 CCD 上占多少个像素、质心提取精度能到多少。常见做法是先定四个物理量：波长 lambda、微透镜间距 pitch、微透镜焦距 f、CCD 像素尺寸 dccd。我一般还会加一个子孔径数量 n_sub，它决定波前采样密度。下表是我常用的一组基础参数，后面所有代码都基于这组值展开：

这几个参数算出来有一个关键中间量：单个微透镜的衍射极限焦斑主瓣尺寸。对方形孔径近似，第一个零点位置在 lambda*f/pitch，约等于 0.071 mm，也就是差不多 8 个 CCD 像素。这个数字直接告诉你能用多大窗口做质心提取，也让后面仿真里的 PSF 模板有了物理尺度。另一个要提前算的是微透镜间距在 CCD 上的像素数，即 pitch/dccd，大约 16.7 个像素。这个值后面做网格索引时要用。如果这个值算出来太小，比如小于 5 个像素，焦斑和焦斑会在 CCD 上重叠，实验里这种布局就很难标定，仿真里也一样。

2.2 用 PSF 模板合成焦斑图：又快又不引入衍射仿真的采样陷阱

生成哈特曼焦斑图有两种主流做法。第一种是逐子孔径做 FFT 衍射仿真：把每个子孔径内的复振幅抽出来，做二维 FFT 得到焦面光强。这种做法物理上更严谨，能让焦斑带上子孔径内高阶像差的细节，但也特别容易踩采样坑。比如前面那组参数，每个子孔径口径 0.15 mm，如果用 64×64 点表示，FFT 后焦面相邻像素的物理间隔是 lambdaf/(Ndelta)，按傅里叶变换的坐标换算大约 0.005 mm 量级，需要把零填充做到 1024 以上才能让焦面采样和 CCD 像素尺寸匹配。每帧图要算上千次大尺寸 FFT，仿真一次要等很久。

第二种做法是用预计算的 PSF 模板来做图像合成，这也是我做标定算法验证时的首选。核心思想是：标定本来就要确定质心偏移和斜率的关系，而质心只依赖焦斑整体位置和光强分布，子孔径内的局部高阶相位对质心的贡献在标定精度需求下可以忽略。于是先计算一个方形孔径的理想焦斑模板，再用插值平移到每个微透镜对应的 CCD 位置上。每帧图生成时间从分钟级降到秒级，而且质心偏移与斜率保持严格的线性，非常适合用来验证标定算法本身。

2.3 生成一帧带噪声焦斑图：代码与参数说明

下面的代码生成了完整的一帧哈特曼焦斑图。我们假设入射波前是一个由离焦和像散组成的波前，先求出每个子孔径中心处的局部斜率，再把模板平移过去，最后加泊松噪声和暗电流噪声。代码较长但每一段都对应前面说的一个环节，建议照着敲一遍再改参数。

%-------------------------------------------------------------------------- % 1. 基础参数 %-------------------------------------------------------------------------- lambda = 0.532e-3; % 波长，单位 mm pitch = 0.15; % 微透镜间距，单位 mm f = 20; % 微透镜焦距，单位 mm dccd = 0.009; % CCD 像素尺寸，单位 mm n_sub = 32; % 每边子孔径数量 % 模板与画幅参数 M = 64; % 模板边长（像素），覆盖约 18 个艾里斑主瓣 padding= 10; % 画幅四周留白（像素） img_sz = round(n_sub * pitch / dccd) + 2 * padding; %-------------------------------------------------------------------------- % 2. 生成方形孔径的 PSF 模板（归一化到总能量为 1） %-------------------------------------------------------------------------- x_f = linspace(-3, 3, M) * lambda * f / pitch; % 焦面坐标，单位 mm [X_f, Y_f] = meshgrid(x_f, x_f); arg_x = pi * pitch * X_f / (lambda * f); % 方形孔径的 sinc 参数 arg_y = pi * pitch * Y_f / (lambda * f); psf = (sin(arg_x + 1e-12) ./ (arg_x + 1e-12)) .* ... (sin(arg_y + 1e-12) ./ (arg_y + 1e-12)); psf = psf .^ 2; % 光强 = |sinc|^2 psf = psf / sum(psf(:)); % 归一化，保证每帧总光强一致 %-------------------------------------------------------------------------- % 3. 构造入射波前：离焦 + 像散，单位为“波长” %-------------------------------------------------------------------------- [Xp, Yp] = meshgrid(linspace(-1, 1, n_sub), linspace(-1, 1, n_sub)); R2 = Xp.^2 + Yp.^2; W = 0.8 * (2 * R2 - 1) + 0.25 * (Xp.^2 - Yp.^2); % 系数单位：波长 % 求每个子孔径中心处的局部斜率，单位：波长/归一化坐标 sx = 2 * 0.8 * 2 .* Xp + 0.25 * 2 .* Xp; % dW/dx sy = 2 * 0.8 * 2 .* Yp - 0.25 * 2 .* Yp; % dW/dy % 斜率换算到焦面位移：位移 = f * 斜率，注意斜率要用物理坐标 % 归一化坐标跨度为 2，对应物理宽度 n_sub*pitch，所以物理斜率要乘以因子 scale = n_sub * pitch / 2; % 归一化坐标→物理坐标 sx_mm = sx / scale; % 单位：mm/mm（无量纲） sy_mm = sy / scale; shift_px_x = f * sx_mm / dccd; % 质心平移，单位像素 shift_px_y = f * sy_mm / dccd; %-------------------------------------------------------------------------- % 4. 逐子孔径把模板放进画幅 %-------------------------------------------------------------------------- img = zeros(img_sz, img_sz); for i = 1:n_sub for j = 1:n_sub % 微透镜中心的 CCD 像素坐标（留白后） cx = padding + (j - 0.5) * pitch / dccd + shift_px_x(i, j); cy = padding + (i - 0.5) * pitch / dccd + shift_px_y(i, j); x0 = round(cx - M/2) + 1; y0 = round(cy - M/2) + 1; x1 = x0 + M - 1; y1 = y0 + M - 1; if x0 < 1 || y0 < 1 || x1 > img_sz || y1 > img_sz continue; % 边缘点丢弃，后面正好可以当缺失点练手 end img(y0:y1, x0:x1) = img(y0:y1, x0:x1) + psf; end end %-------------------------------------------------------------------------- % 5. 加泊松光子噪声 + 高斯暗噪声 %-------------------------------------------------------------------------- peak_photons = 3000; % 光斑峰值光子数 img_photons = img / max(img(:)) * peak_photons; noisy = poissrnd(img_photons); % 光子噪声 noisy = noisy + 5 * randn(size(noisy)); % 暗电流与读出噪声 noisy(noisy < 0) = 0; figure; imagesc(noisy); axis image; colormap(gray); title('Simulated Shack-Hartmann spot field');

这段代码里要特别留意的参数是scale。因为W是在归一化坐标 [-1,1] 上构造的，而物理坐标的斜率必须用“单位 mm 的波前变化 / 单位 mm 的横向坐标”来表示，所以要做一次坐标缩放。如果忘了这一步，后面算出来 f_eff 会差一个因子，整个标定就废了。peak_photons直接控制信噪比，3000 光子对应比较干净的实验条件，标定算法应该轻松处理；把它降到 300，你会发现质心提取就开始抖了，这正是不常用到但很有价值的“压力测试”。

3. 从焦斑图到有序网格：质心提取与几何标定怎么做

3.1 质心提取：为什么阈值重心法比纯重心法更耐噪

哈特曼传感器每个子孔径的波前斜率，是靠焦斑质心相对于参考位置的偏移量来测量的。所以质心提取精度直接决定了后面所有结果。最简单的重心法是把窗口内所有像素的光强加权求平均，公式是 cx = sum(I*x)/sum(I)。这个公式在背景干净、光斑均匀时没问题，但 CCD 上总有暗电流和读出噪声，而这些噪声的位置随机，会把质心往窗口中心方向拉，产生系统性的低估。

我一般建议用阈值重心法：先减去一个背景阈值，只对超过阈值的像素做加权。这样把噪声基底扣掉，质心估计的偏差会明显减小。还有一种做法是用以粗质心为中心的窗口迭代，先算一次全窗口质心，再以这个质心为中心重新开一个更小的窗口再算一次，能有效去掉相邻光斑的串扰。阈值怎么选取决于噪声水平，一个实用的规则是取窗口四周 8 个像素的平均值再乘 1.5 到 2 倍。

下面的函数可以作为你的通用质心工具，它把“阈值+固定窗口”结合在一起，窗口宽度可以直接传入。对于 16 像素间距的焦斑图，窗口取 12 左右比较安全；如果你把窗口开到 18，就可能把邻居光斑卷进来。

function [cx, cy] = centroid_thr(img, win, thr, px, py) % img : 整帧哈特曼图 % win : 窗口宽度（像素，奇数） % thr : 背景阈值 % px,py: 光斑粗略位置（可以是网格预测值） % % 输出：子像素质心坐标 half = floor(win / 2); x0 = round(px); y0 = round(py); % 取窗口 r1 = max(y0 - half, 1); r2 = min(y0 + half, size(img, 1)); c1 = max(x0 - half, 1); c2 = min(x0 + half, size(img, 2)); sub = img(r1:r2, c1:c2); % 扣阈值 sub = sub - thr; sub(sub < 0) = 0; % 重心 total = sum(sub(:)); if total == 0 cx = px; cy = py; return; end % 行列加权，注意坐标原点在窗口左上角 [yy, xx] = meshgrid(1:size(sub,1), 1:size(sub,2)); cx = c1 - 1 + sum(sum(sub .* xx)) / total; cy = r1 - 1 + sum(sum(sub .* yy)) / total; end

调用这个函数之前，必须先把每个光斑的粗略位置找出来。最省事的做法是用imregionalmax找到局部极大值点，或者对图像做一次阈值二值化后用regionprops得到每个连通域的质心作为初值。粗位置不需要很准，只要落在光斑主瓣内即可，因为后面窗口内的重心会把它修正到子像素精度。这个“粗定位+细质心”的两步法是工程里最稳的组合。

3.2 从无序光斑坐标恢复栅格索引：估计间距、旋转角并匹配

仿真里每个光斑生成时都带着(i,j)索引，但真实标定中拿到的是乱序坐标列表，必须把它们重新排成二维网格，这一步通常被称为网格索引或点阵匹配。常见做法是利用先验的微透镜间距做最近邻匹配：先用所有光斑坐标估计实际的像素间距 d_px 和旋转角 theta，然后把每个点归算到整数栅格上。估计 d_px 最实用的方法是对任意一个光斑找它最近的几个邻居，把距离的众数或中位数作为间距。旋转角则可以从两个相邻点的连线角度来估计。

这里给出一个简单的匹配函数。它能处理轻微旋转，但不能处理很大的畸变；如果实验里有明显桶形畸变，得先用多项式拟合做畸变矫正，那是另一个话题了。

function idx = match_grid(pts, d_px, theta) % pts : Nx2 矩阵，每行一个光斑像素坐标 [x, y] % d_px : 微透镜间距的像素值（估计值） % theta : 栅格旋转角（弧度） % idx : Nx2 整数矩阵，每个点的栅格序号 [col, row] % 坐标旋转到栅格坐标系 pts_rot = [cos(theta) sin(theta); -sin(theta) cos(theta)] * pts'; xr = pts_rot(1,:); yr = pts_rot(2,:); % 归一化到整数栅格，并取整 col = round(xr / d_px); row = round(yr / d_px); idx = [col(:), row(:)]; end

这个函数隐含了一个前提：坐标原点恰好落在某个微透镜中心上。真实情况不会这么理想，所以匹配前要先估计原点偏移。估计方法是把所有点的坐标取模 d_px 之后的分布做统计，取众数作为偏移；或者干脆把第一个点的坐标当作原点，先匹配一次，再用残差修正。匹配完成后，最好画一张plot检查：每个光斑位置标上它的栅格序号，看有没有错排。错排的典型表现是某些行或列的序号发生跳变，这时把 d_px 或 theta 微调一下重跑即可。

3.3 参考帧标定与坐标关系：把像素坐标变成微透镜坐标

哈特曼标定里最容易漏的一步是参考帧。你需要先让一束接近理想的平面波入射，记录所有光斑的质心位置，这组位置就叫参考质心。之后测量未知波前时，每个光斑的质心偏移量是当前质心坐标减去参考质心坐标，而不是减去栅格中心。为什么不能直接减栅格中心？因为微透镜阵列在 CCD 上的投影本身就有旋转和位移，还有少量畸变；栅格中心是拟合出来的理想位置，而参考质心包含了真实的安装状态。用参考质心做差，等于把系统的固定偏移全部扣掉了。

参考帧通常取多帧平均，比如拍 20 帧平面波，把对应子孔径的质心做平均，能把随机噪声压下去。这一步做完，你就拥有了一张和微透镜一一对应的参考质心表ref_cx、ref_cy。它的物理含义是“当入射波前斜率为零时光斑应该在的位置”。后续所有斜率计算都变成：

delta_x_px = current_cx - ref_cx delta_y_px = current_cy - ref_cy

到这里，像素坐标和微透镜索引的映射关系已经建立。下一步要解决的就是那个比例系数：一个像素的偏移对应多少波前斜率。

4. 斜率标定与波前重建：把标定系数 f_eff 落实成算法

4.1 标定系数从哪来：用手动加入的倾斜量反推 f_eff

前面说过，质心偏移量 Δx 与子孔径处波前斜率 sx 的关系是 Δx = f_eff * sx，这里的 f_eff 单位是像素。理论上它等于微透镜焦距 f 除以 CCD 像素尺寸 dccd，但实际装配中 f 可能因为胶层、压圈等原因和设计值有百分之几的偏差，所以要用标定把它测出来。

标定方法其实很简单。在仿真里，往入射波前中加入一个已知的倾斜量，也就是给波前加一个线性项 aX + bY，其中 a 和 b 是已知的波前斜率。这时所有子孔径应当产生一致的质心偏移，偏移量除以已知斜率就是 f_eff。写成代码，对第 2 章生成的那帧图做质心提取后，把平均偏移量和已知系数拿来拟合：

% 手动加入倾斜：波前 = a*X + b*Y，a、b 为已知斜率 a = 0.2; % 单位：波长/归一化坐标 b = 0.1; scale = n_sub * pitch / 2; % 和第 2 章保持一致 sx_known = a / scale; % 物理斜率，单位 mm/mm sy_known = b / scale; % 假设你已经提取了所有光斑质心 cur_cx, cur_cy 和参考质心 ref_cx, ref_cy delta_x_px = median(cur_cx - ref_cx); % 用中位数抗异常点 delta_y_px = median(cur_cy - ref_cy); f_eff_x = delta_x_px / sx_known; % 单位像素 f_eff_y = delta_y_px / sy_known; f_eff = (f_eff_x + f_eff_y) / 2;

用中位数而不是平均值，是因为边缘子孔径可能已经跑到画幅外，质心变成异常值，平均值会被带跑。这样标出来的 f_eff 包含了整个系统从微透镜到 CCD 的等效比例，后续重建直接拿它做斜率换算就行。注意这个数值和你用 f/dccd 直接算的理论值 2222 像素会有几个百分点的差异，因为仿真里的采样和模板离散化会引入一点偏移。这正是“标定”要解决的实际问题。

4.2 用 Zernike 模式拟合重建波前：前六项设计矩阵的构造

拿到每个子孔径的斜率之后，波前重建有两条路：区域法和模式法。区域法是把相邻子孔径的斜率做最小二乘积分，得到离散相位点；模式法是把波前表示成 Zernike 多项式的线性组合，用斜率数据直接拟合系数。模式法更适合需要关心像差分类的场景，而且实现起来更直观：你要做的是构造一个设计矩阵 A，第 i 行是各个 Zernike 项在对应子孔径位置处的 x 方向导数和 y 方向导数。

Zernike 多项式有无穷多项，工程上常见做法是先取前 6 项或前 10 项。前 6 项对应活塞、倾斜、离焦和像散，已经覆盖了大部分装调误差和基础像差。以下代码构造前 6 项在指定子孔径坐标处的函数值和导数值，然后组装设计矩阵并做最小二乘拟合。

% 子孔径坐标（归一化到单位圆，第 2 章的 Xp, Yp 已经是归一化的 [-1,1]） x = Xp(:); y = Yp(:); N = length(x); % 前 6 项 Zernike（Noll 顺序） % Z1 = 1 % Z2 = 2x % Z3 = 2y % Z4 = sqrt(3) * (2r^2 - 1) % Z5 = 2*sqrt(6) * xy % Z6 = sqrt(6) * (x^2 - y^2) % % 注意：活塞项 Z1 的斜率恒为 0，无法从斜率数据恢复，直接去掉 % 构造斜率设计矩阵，每行 = [dZj/dx, dZj/dy] % 对每一项求导，结果按列排列 A = zeros(2 * N, 5); % 5 项，两个方向 b = zeros(2 * N, 1); % 观测值：斜率 % Z2 = 2x A(1:N, 1) = 2; % dZ2/dx A(N+1:end, 1) = 0; % dZ2/dy % Z3 = 2y A(1:N, 2) = 0; A(N+1:end, 2) = 2; % Z4 = sqrt(3) * (2x^2 + 2y^2 - 1) A(1:N, 3) = sqrt(3) * 4 .* x; A(N+1:end, 3) = sqrt(3) * 4 .* y; % Z5 = 2*sqrt(6) * xy A(1:N, 4) = 2 * sqrt(6) .* y; A(N+1:end, 4) = 2 * sqrt(6) .* x; % Z6 = sqrt(6) * (x^2 - y^2) A(1:N, 5) = sqrt(6) * 2 .* x; A(N+1:end, 5) = sqrt(6) * (-2) .* y; % 观测向量：前面算出的物理斜率（单位已折算） sx_measured = (cur_cx - ref_cx) / f_eff; % 像素偏移 / 像素每斜率 sy_measured = (cur_cy - ref_cy) / f_eff; b(1:N) = sx_measured(:); b(N+1:end) = sy_measured(:); % 最小二乘求解 Zernike 系数 z_fit = A \ b;

最小二乘在这里的行列式条件数一般不大，直接用反斜杠求解即可。有一点要提醒：A的行数远大于列数，这是一个超定问题，反斜杠默认给出最小二乘解。如果你想给边缘子孔径更低的权重，可以通过给A和b乘一个对角权重矩阵来实现，这在后面谈边缘丢点时很有用。

4.3 重建精度评估与残差分析：RMS 和 PV 计算

标定做没做准，必须靠残差说话。所谓残差，是把拟合出的 Zernike 系数代回设计矩阵，算出“预测斜率”，再和“实测斜率”相减。理想情况下残差应该只是随机噪声，而不是有结构性的趋势。如果有趋势，说明你少标了旋转角、窗口不对、或者有子孔径点阵错位。

斜率残差可以直接用向量算：

residual = A * z_fit - b; res_slope = reshape(residual(1:N), n_sub, n_sub); % x 方向残差 % 残差的标准差 sigma_slope = std(residual);

不过，斜率残差不够直观，大家更关心的是波前重建误差。方法是在网格上把真值波前和重建波前对比。重建波前由 z_fit 和 Zernike 基函数合成，真值波前就是仿真一开始构造的 W。两者相减得到误差面，再算 PV（峰谷值）和 RMS（均方根）：

% 重建波前 W_rec = zeros(size(Xp)); W_rec = W_rec + z_fit(2) * 2 .* Xp; % Z2 W_rec = W_rec + z_fit(3) * 2 .* Yp; % Z3 W_rec = W_rec + z_fit(4) * sqrt(3) .* (2 * (Xp.^2 + Yp.^2) - 1); W_rec = W_rec + z_fit(5) * 2 * sqrt(6) .* Xp .* Yp; W_rec = W_rec + z_fit(6) * sqrt(6) .* (Xp.^2 - Yp.^2); err = W_rec - W; % 单位：波长 % 去掉活塞项的影响 err = err - mean(err(:)); wavefront_rms = sqrt(mean(err(:).^2)); wavefront_pv = max(err(:)) - min(err(:)); fprintf('RMS = %.3f wave, PV = %.3f wave\n', wavefront_rms, wavefront_pv);

一个标定流程如果做得好，在 3000 光子噪声水平下，RMS 应该小于 0.05 个波长。如果这个值超过 0.1 个波长，别急着调重建算法，先回头检查：参考帧是不是用多帧平均了？f_eff 是不是用中位数拟合的？网格索引有没有错排？大多数情况下问题都出在这些“看起来早就处理过”的地方。

5. 避坑：哈特曼标定仿真里最常见的 4 个“玄学”

5.1 焦斑图仿真看着清晰，质心却系统性偏向像素中心

表现：图像上光斑明明很锐利，但把大量随机位移的仿真帧拿去统计质心，发现质心总是往光斑所在像素的整数中心靠拢，偏移量被低估。

原因：用 PSF 模板平移时如果直接用round做了像素取整，再做质心提取时亚像素信息已经丢了。另一个常见来源是质心窗口内像素太少，光斑只覆盖 3×3 或 4×4 个像素，离散化误差会把真正的连续重心朝像素中心拉。

解决：生成焦斑图时用亚像素平移，也就是在x0 = round(cx - M/2)这一步保留小数部分的偏移，用插值或频域移位实现亚像素平移。如果只是想验证标定算法，可以用双线性插值把模板整体平移；如果追求高精度，直接用imtranslate并指定双三次插值。质心窗口最好保证光斑主瓣直径占到窗口宽度的 60% 以上，让窗口内至少有 5×5 个有效像素。

5.2 微透镜旋转角没标，重建波前多出一倍像散

表现：重建出来的像散系数明显偏大，而且无论怎么改变真实像散，拟合结果里总有一对固定大小的像散项。

原因：微透镜阵列和 CCD 像素行之间有一个小角度旋转。如果你忽略旋转角，直接按理想栅格去做索引和坐标映射，每个子孔径的实际坐标和理想坐标之间会产生一个随位置线性变化的偏差。这个偏差的梯度恰好长得像像散，会被 Zernike 拟合误解出来。

解决：匹配栅格前用第 3.2 节的match_grid先估计旋转角再做索引。旋转角的估计办法很简单：取相邻两个光斑的连线夹角，和理想栅格的角度做差。在仿真里你可以在生成焦斑图前给整个微透镜阵列的坐标加一个 0.5 度左右的旋转，然后跑一遍完整流程，观察像散残差的变化，这是验证标定流程是否稳健的最快方法。

5.3 质心窗口没开对：截断与串扰二选一

表现：质心提取结果忽大忽小，波前重建残差在个别子孔径上突然飙升。

原因：窗口开小了会截断光斑主瓣，质心计算丢掉一部分有效能量，导致偏移量被低估；窗口开大了则把相邻子孔径的光斑边缘卷进来，引入“额外光强”把质心拉向邻居方向。两者看起来都是质心不准，但误差规律完全不同。

解决：先用imagesc单帧图量一下相邻光斑中心的像素间距，再量光斑主瓣直径。窗口宽度取主瓣直径的 1.5 到 2 倍，同时必须小于光斑间距的 0.7 倍。以第 2 章参数为例，间距约 16.7 像素，主瓣直径约 8 像素，窗口取 12 是安全的。调整完窗口后，重跑标定流程，看残差是否明显下降。

5.4 边缘子孔径丢几个点，重建残差立刻爆表

表现：画幅边缘的光斑因为超出 CCD 范围被丢弃，拟合出的波前在边缘区域出现剧烈振荡，RMS 数值一下子涨上去。

原因：Zernike 多项式在单位圆边缘的导数值往往比中心大很多，如果边缘子孔径缺失，这些高权重位置没有观测数据约束，最小二乘解为了平衡残差会把边缘的系数推得很大，产生边缘振铃。

解决：第一个方案是拟合时给边缘缺失的子孔径剔除，直接去掉那些没有有效质心的行。第二个方案是采用加权最小二乘，给靠近边缘的子孔径一个较小的权重，避免边缘主导拟合结果。第三个方案更暴力：重新设计画幅尺寸，把 CCD 尺寸放大到比微透镜阵列投影大 10% 以上，保证边缘子孔径不会轻易跑丢。如果你要做的实验里边缘光斑注定会出画幅，那仿真里也刻意丢弃边缘点，用加权拟合来验证重建算法的鲁棒性，这样后面拿到实测数据才不会手忙脚乱。

6. 闭环验证：注入随机像差，确认你的标定流程真的能用

6.1 用随机 Zernike 系数自检标定系统

整套标定链路搭好之后，别急着拿去处理实测数据，先做一次闭环验证：随机生成一组 Zernike 系数当作真值，用仿真链路生成焦斑图，再跑一遍质心提取、网格索引、斜率标定和波前重建，最后把反推的系数和真值比较。这个流程相当于给标定系统做一次“全链路体检”，能一次性发现前面任何一步的问题。

% 随机生成真值系数，离焦和像散为主，倾斜量保持较小 z_true = [0.6 * randn(), 0.1 * randn(), 0.3 * randn(), ... 0.2 * randn(), 0.15 * randn()]; % 对应 Z2~Z6 % 用 z_true 构造波前，得到子孔径斜率 W_truth = z_true(1) * 2 .* Xp + z_true(2) * 2 .* Yp + ... z_true(3) * sqrt(3) .* (2 * R2 - 1) + ... z_true(4) * 2 * sqrt(6) .* Xp .* Yp + ... z_true(5) * sqrt(6) .* (Xp.^2 - Yp.^2); sx_truth = 梯度计算; sy_truth = 梯度计算; % 生成焦斑图（代码与第 2 章一致） img = generate_sh_image(sx_truth, sy_truth, psf, n_sub, ...); % 然后做完整标定流程，得到 z_fit z_rec = run_full_calibration(img, ref_cx, ref_cy, f_eff); % 对比 err_z = z_rec - z_true; fprintf('Zernike coefficient RMS error: %.4f wave\n', ... sqrt(mean(err_z.^2)));

判断标准：如果系数误差都在 0.05 个波长以内，标定流程是健康的；如果某个方向误差特别大，看是不是那个方向的斜率灵敏度接近零，比如 Z6 在 x 轴附近时 y 方向导数为零，此时该方向观测不足，拟合会退化。真值生成时不妨多跑 20 组随机系数，统计误差分布，这样能得到比单次更可靠的结论。仿真里的闭环验证能给你一个“基线”：拿到实测数据后，如果重建残差远差于这个基线，就能确定问题出在硬件装配或者图像采集，而不是算法链路本身。

我自己的习惯是每次改完代码，都先跑这一遍闭环，误差超标就回头查参数，超标几次之后慢慢就形成了条件反射：残差大先查参考帧，再查旋转角，最后查窗口。这套流程给不了你什么花哨的技巧，但能让你在真正面对一片模糊焦斑图时，第一时间知道该动哪里。希望帮到你。

本文还有配套的精品资源，点击获取