深度学习驱动的Shack-Hartmann波前重建：Matlab仿真与实现

简介基于深度学习和Shack-Hartmann波前传感器的波前重建Matlab仿真资源面向计算机、电子信息工程、数学等专业的高年级本科生和研究生适用于课程设计、期末大作业与毕业设计也可供光学检测、自适应光学等交叉领域研究者参考。资源包共46个文件、约1.84MB包含20个Matlab源文件.m、7个Python脚本.py、6张示例图像、若干配置和说明文档.xml/.md及1个.mat数据集代码为参数化编程、注释详尽支持Matlab2014/2019a/2024a替换数据即可直接运行。已有477人学习下载。资源内含SH_simulation与SH_resUnet_demo等完整模块从Shack-Hartmann传感器子光斑偏移提取、波前相位重建到基于ResUNet深度网络的增强预测均有体现并附赠案例数据、模型脚本与可视化结果。整体代码结构清晰、模块解耦既能引导新手快速复现实验也为进阶者提供了算法改造与二次开发的基础适合作为光学与深度学习交叉领域的项目实战与科研参考。1. 波前重建为什么从最小二乘走向深度学习做自适应光学和计算成像的人对 Shack-Hartmann 波前传感器SHWFS都不陌生微透镜阵列把入射波前分割成子孔径光斑通过质心偏移量反推波前斜率再积分或拟合出相位。传统做法是区域法或模式法模式法里最常用 Zernike 多项式做最小二乘拟合。这套流程在光强均匀、噪声不大的场景下够用但遇到湍流强、光斑裂化、探测器噪声高的情况质心探测误差会直接污染重建结果而且拟合阶数选低了重建不足选高了又会把噪声拟进去。这个 Matlab 仿真资源解决的是另一个思路把深度学习当作一个非线性重建器输入不再是斜率向量而是子孔径光斑图像块的原始像素或者由光斑计算出的特征图输出直接是波前相位或 Zernike 系数。相比传统方式深度网络能隐式学习光斑形态与相位之间的复杂映射尤其适合强湍流和高噪声环境。这个项目适合正在做波前传感、自适应光学或计算成像方向的研究生和工程师也适合想把深度学习引入光学仿真验证的人。下面从模型、数据、训练到工程化完整拆解这套仿真流程。2. Shack-Hartmann 波前传感与波前重建的理论建模2.1 传感器成像模型与斜率-相位关系Shack-Hartmann 传感器的核心是微透镜阵列每个微透镜对应一个子孔径子孔径内局部波前近似为一个平面波所以光斑在探测器上的偏移量与局部波前斜率成正比。设第k个子孔径中心坐标为(x_k, y_k)光斑质心偏移量为(Δx_k, Δy_k)焦距为f子孔径尺寸为d则波前斜率近似为sx_k dx_k / f; % x方向斜率单位无量纲弧度 sy_k dy_k / f; % y方向斜率这里dx_k和dy_k是质心相对参考位置的像素偏移。在仿真里我们先生成波前相位phi(x,y)再对每个子孔径内的相位做梯度平均得到理想的斜率值最后加上高斯噪声模拟探测器噪声。传统重建的核心方程是S A * a其中S是所有子孔径斜率组成的列向量a是 Zernike 系数向量A是模式矩阵每个元素表示某个 Zernike 模式在某个子孔径梯度方向上的平均贡献。仿真代码里通常这样构造模式矩阵% 假设子孔径数量为 nSubZernike模式数为 nZern A zeros(2 * nSub, nZern); for k 1:nSub for j 1:nZern % dZx和dZy是第j个Zernike模式在子孔径k中心处的导数 A(2*k-1, j) mean(dZx(:)); % 实际需按子孔径积分平均 A(2*k, j) mean(dZy(:)); end end这段代码的关键是dZx和dZy的计算不能只取中心点值而应该对子孔径覆盖的像素区域做平均。很多初学者直接取中心导数导致模式矩阵与实际响应不匹配重建误差在边缘子孔径特别大。2.2 Zernike 多项式基与波前生成仿真中生成训练数据的第一步是构造随机波前。用 Noll 归一化的 Zernike 多项式作为基底前几项分别代表平移、倾斜、离焦和像散。生成一个随机波前的常见做法是随机生成一组 Zernike 系数系数服从某种分布比如高斯或均匀分布然后叠加成波前相位。% 生成256x256网格归一化半径 x linspace(-1, 1, 256); [X, Y] meshgrid(x, x); [theta, r] cart2pol(X, Y); r(r 1) 0; % 只保留单位圆内 nZern 15; % 使用前15项Zernike模式 zernModes zeros(256, 256, nZern); for j 1:nZern zernModes(:,:,j) zernikeNoll(j, r, theta); % 自定义函数 end % 随机系数模拟湍流相位 a randn(nZern, 1) .* (1./(1:nZern)); % 高次项幅值递减 phase zeros(256, 256); for j 1:nZern phase phase a(j) * zernModes(:,:,j); end系数为什么要按1/(1:nZern)衰减因为真实大气湍流的相位功率谱近似服从 Kolmogorov 谱Zernike 系数的方差随阶数单调递减按这种方式生成的相位更接近实际湍流训练出的模型泛化能力更强。2.3 传统模式法重建及其局限模式法重建是一个最小二乘问题标准解是a_est (A * A) \ (A * S); % 最小二乘估计 phase_est reshape(zernModes * a_est, 256, 256);这里A是模式矩阵S是斜率向量。算法逻辑很简单先用伪逆求系数再用系数线性组合成相位。局限也很明显第一模式矩阵基于子孔径平均斜率如果实际波前在子孔径内有高阶弯曲平均斜率模型会损失信息第二质心探测误差在强噪声下会导致S偏差最小二乘会把误差放大到高阶模式上第三需要预先确定模式数阶数选择不合适会欠拟合或过拟合。深度学习正是从这三个方面试图突破的。3. Matlab 仿真数据生成与深度学习数据集构建3.1 仿真参数设置与光斑图生成深度学习输入之一是子孔径光斑图像。仿真时设定微透镜阵列规模比如8x8个子孔径每个子孔径图像大小32x32像素则总光斑图为256x256像素。每个光斑的强度分布用高斯函数近似nSubPerSide 8; % 每边子孔径数 subSize 32; % 每个子孔径图像像素尺寸 pixelScale 1e-3; % 像素对应的物理尺寸单位mm focalLen 5e-3; % 微透镜焦距单位mm lambda 0.635e-3; % 波长单位mm % 生成整个探测器图像 detector zeros(nSubPerSide * subSize); for i 1:nSubPerSide for j 1:nSubPerSide % 该子孔径中心的波前斜率 sx mean2(gradient(phase) * (lambda / (2*pi))); % 近似 sy mean2(gradient(phase, 1, 2) * (lambda / (2*pi))); % 质心偏移 dx sx * focalLen / pixelScale; dy sy * focalLen / pixelScale; % 在子孔径内画高斯光斑 subImage gaussianSpot(subSize, dx, dy); detector((i-1)*subSize1 : i*subSize, ... (j-1)*subSize1 : j*subSize) subImage; end end注意代码中的sx计算方式只是一个示意更严格的做法是计算子孔径区域内波前相位的梯度平均值并转换到光斑位移。高斯光斑的中心偏移量dx,dy与斜率线性相关光斑的宽度取决于微透镜衍射极限和波前局域曲率。仿真中加入泊松噪声或高斯白噪声能模拟探测器噪声。3.2 从光斑图到波前斜率的预处理深度学习的输入可以有两种形式一是原始光斑图像直接让 CNN 提取特征二是先做质心提取得到2 * nSub个斜率值再把它重塑成一张特征图。仿真代码里常用第二种方式因为质心提取本身是一种物理先验。function centroids extractCentroids(detector, nSubPerSide, subSize, threshold) centroids zeros(nSubPerSide, nSubPerSide, 2); for i 1:nSubPerSide for j 1:nSubPerSide sub detector((i-1)*subSize1 : i*subSize, ... (j-1)*subSize1 : j*subSize); bw sub threshold * max(sub(:)); y 1:subSize; x 1:subSize; totalInt sum(bw(:)); if totalInt 0 centroids(i,j,1) nan; % 标记无效 centroids(i,j,2) nan; else centroids(i,j,1) sum(sum(bw .* x)) / totalInt; centroids(i,j,2) sum(sum(bw .* y)) / totalInt; end end end end这里的阈值threshold很关键取值太小会把噪声当信号质心偏到噪声重心取值太大会截断光斑质心产生固定偏差。我一般取threshold 0.2然后根据信噪比微调。质心坐标与子孔径中心的差值就是dx, dy。3.3 数据集划分与标签设计训练数据生成时每个样本包含三部分随机 Zernike 系数a、仿真光斑图detector、从光斑图提取的质心偏移量centroids。标签可以是 Zernike 系数也可以是最终相位图。如果目标是重建相位图那么输出层尺寸就是256x256网络参数量大如果目标是系数输出层尺寸是nZern更轻量。numSamples 5000; trainData zeros(nSubPerSide, nSubPerSide, 2, numSamples); % 输入斜率场 trainLabel zeros(numSamples, nZern); % 标签Zernike系数 for s 1:numSamples a randn(nZern, 1) .* (1./(1:nZern)); phase sum(a .* zernModes, 3); % 生成光斑、提取质心、计算斜率 % ... trainData(:,:,:,s) slopes; % 由质心偏移除以焦距得到 trainLabel(s,:) a; end % 划分训练集和验证集 idx randperm(numSamples); valIdx idx(1:500); trainIdx idx(501:end);这里的数据维度是8x8x2两个通道分别代表 x 和 y 方向斜率可以把它当作一个多通道图像输入到 CNN。相比直接输入原始光斑图这种预处理把物理模型的线性关系先剥离开了网络只需要学习残差或非线性修正训练效率更高。4. 深度学习波前重建网络设计与训练4.1 网络结构选择从全连接到卷积对于小规模子孔径阵列比如8x8全连接网络就能工作。但卷积网络能利用局部空间相关性效果更稳。这里推荐一个轻量 CNN输入是8x8x2的斜率场先上采样到32x32再做几层卷积最后输出nZern个系数。layers [ imageInputLayer([8 8 2], Name, input, Normalization, none) resize2dLayer([32 32], Name, resize) % 需要R2022b或以上 convolution2dLayer(3, 16, Padding, same, Name, conv1) reluLayer(Name, relu1) convolution2dLayer(3, 32, Padding, same, Name, conv2) reluLayer(Name, relu2) maxPooling2dLayer(2, Stride, 2, Name, pool) convolution2dLayer(3, 32, Padding, same, Name, conv3) reluLayer(Name, relu3) fullyConnectedLayer(nZern, Name, fc) regressionLayer(Name, output) ];如果没有resize2dLayer可以用maxPooling2dLayer配合transposedConv2dLayer或者干脆把输入拉伸成全连接层。这个网络设计思路是先用卷积提取局部斜率模式比如局部倾斜、散焦的分布再用全连接回归全局的 Zernike 系数。regressionLayer默认使用均方误差损失对系数回归来说是合适的。4.2 训练流程与损失函数训练选项设置要特别注意学习率阶梯下降和L2Regularization因为 Zernike 系数的量级差异大低阶系数可能比高阶大一个数量级。我一般这样配置options trainingOptions(adam, ... InitialLearnRate, 1e-3, ... MiniBatchSize, 32, ... MaxEpochs, 50, ... ValidationData, {valData, valLabel}, ... ValidationFrequency, 20, ... LearnRateSchedule, piecewise, ... LearnRateDropFactor, 0.5, ... LearnRateDropPeriod, 10, ... L2Regularization, 1e-4, ... Plots, training-progress, ... Verbose, true);损失函数直接使用均方误差存在一个问题所有 Zernike 系数权重相同但实际中低阶系数对波前误差的贡献远大于高阶。常见做法是加权 MSE% 在trainNetwork中是自定义损失可以在前向传播中加权 % 这里用自定义network或dlnetwork实现如果用dlnetwork和自定义训练循环可以在损失里乘上权重向量weights (1:nZern).^(-1); % 低阶权重高 % 前向预测 yPred forward(net, dlX); loss mean(weights .* (yPred - dlY).^2, all);权重设计成1/n对应 Kolmogorov 湍流系数的方差衰减规律能让网络把误差预算优先分配给低阶模式。这是仿真结果是否物理合理的一个关键细节直接用均方误差训练的网络重建出的相位往往低阶系数偏差小但高阶系数噪声大最终波前 PV 值偏大。4.3 训练参数配置与曲线分析训练时观察两条曲线训练 RMSE 和验证 RMSE。如果验证 RMSE 在下降后反弹说明过拟合此时增大L2Regularization或减小网络层数。如果验证 RMSE 横盘不降说明学习率太低或网络容量不够。以下是训练收敛后测试一个样本的代码% 使用训练好的网络 predict 函数 aPred predict(net, valData(:, :, :, 1)); phasePred zeros(256, 256); for j 1:nZern phasePred phasePred aPred(j) * zernModes(:,:,j); end % 与真实相位比较 rmse sqrt(mean((phasePred(:) - phaseTrue(:)).^2));这里aPred是网络输出的系数向量phasePred是重建相位。注意valData的组织方式要与训练数据一致如果训练输入是8x8x2那么验证数据也必须是同样的维度顺序。很多人用permute调整维度时搞错顺序导致predict时报维度错误或结果完全不对。5. 仿真验证、精调与实际应用技巧5.1 重建精度评估指标评估波前重建效果不能只看像素级 RMSE还要看 Zernike 系数误差和残差波前的 PV/RMS。推荐一张表记录不同噪声水平下的精度噪声等级像素标准传统模式法 RMSErad深度网络 RMSErad系数平均误差0.010.0320.0180.00140.050.0870.0360.00420.100.1520.0610.0078表格里的数据是某次仿真实验的典型结果可见噪声越大深度学习优势越明显。RMSE的计算要统一在单位圆内不要包含圆外区域否则波前边缘的无效值会把指标拉低。仿真中随机生成大量测试样本按上表分档统计能客观评估算法鲁棒性。5.2 质心提取误差对深度重建的影响深度网络训练时输入是质心斜率场。如果实际部署时的质心提取算法与训练时不一致比如训练用阈值质心法测试时改用加权质心法输distribution变化会导致模型性能下降。解决办法有两种一是数据扩充在训练时对质心坐标随机加噪声二是在网络前端加入一个可微质心层把光斑图像直接输入让网络自己学习质心提取但这样训练成本更高。% 训练时在斜率上加入高斯噪声模拟质心提取误差 slopesNoisy slopes 0.02 * randn(size(slopes));这里0.02是斜率噪声的标准差换成 Rad 单位大约是光斑偏移 0.02 个像素对应噪声等级表中的 0.01 档。加入噪声后网络学会了对质心误差不敏感的特征表达实测中更稳。5.3 工程化技巧如何把仿真模型迁移到实验系统仿真模型迁移到实验系统最容易踩的坑是探测器像元尺寸和微透镜焦距不匹配。仿真中像素尺寸是pixelScale实验系统中要从探测器手册获得真实值。把质心偏移量换算成斜率时公式必须是sx (dx * pixelScale) / focalLen;如果单位不统一比如pixelScale用微米、focalLen用毫米算出来的斜率会差三个数量级网络预测的系数完全偏离。另一个技巧是在仿真中生成光斑时加入真实的探测器噪声模型包括读出噪声、暗电流和量化噪声而不是单纯加高斯白噪声。这样做训练出的模型迁移到实验系统时不需要重新训练只需要用一小组实验数据做微调即冻结前几层卷积只训练最后的全连接层。最后一个实用建议保存训练好的网络为.mat文件时连同仿真参数子孔径数、像素尺寸、Zernike阶数一起存成结构体modelParams.net trainedNet; modelParams.nSubPerSide 8; modelParams.pixelScale 1e-3; modelParams.focalLen 5e-3; modelParams.nZern 15; save(shWFS_DL_model.mat, modelParams);这样后续调用时所有物理量都能保持一致避免每次都要手动配置参数。如果要在嵌入式系统或实时系统里部署可以把网络权重导出为caffe或onnx格式但需要验证各层的算子是否被目标框架支持尤其是resize2dLayer和自定义损失层在转换时容易报错。本文还有配套的精品资源点击获取