简介：本资源是一套面向计算机、电子信息工程及数学等专业本科生与研究生的波前重建教学仿真工具，聚焦Shack-Hartmann传感器原理与深度学习融合建模，解决光学波前相位实时重建中的噪声鲁棒性与计算效率问题，适用于课程设计、期末大作业及毕业设计等实践场景。压缩包共46个文件，含20个核心Matlab（.m）脚本实现参数化波前仿真与神经网络训练流程，7个Python（.py）文件支持数据预处理与模型验证，6张PNG图像用于结果可视化，另有MATLAB数据文件（.mat）、说明文档（README.md、.txt）及配置文件，整体仅1.84MB，轻量易部署。已有478人学习下载，代码注释详尽、参数接口清晰、案例数据即开即用，特别适合初学者快速理解波前传感物理建模、SH光斑位移反演、ResUNet等深度网络结构设计及Matlab+Python协同仿真工作流。
1. 这不是传统Zernike拟合：用ResUNet在Matlab里跑通Shack-Hartmann波前重建，3分钟出相位图，新手改两行参数就能复现
你手头有一组Shack-Hartmann传感器拍下来的光斑图，想反推原始波前相位——但传统Zernike多项式拟合在高阶像差、低信噪比或非均匀照明下容易发散；而商用波前重建软件要么黑匣子不透明，要么授权贵得离谱。这个Matlab仿真包，就是为这种“卡在实验最后一公里”的场景准备的：它把ResUNet（带跳跃连接的U型卷积网络）完整移植进Matlab环境，不依赖Python后端，不调用外部深度学习框架，所有训练、推理、可视化全在.m文件里闭环完成。附赠的SH_resUNet_demo.m脚本开箱即用，加载img/下的50张模拟光斑图，3分钟内输出重建相位图、残差热力图、PSF对比图三件套。适用对象非常明确——电子信息工程专业做光学课程设计的大三学生、数学系要交数值仿真大作业的研究生、某高校光学实验室刚接手波前校正课题的助教。它不讲抽象的深度学习理论，只解决一个具体问题：给定一组光斑位移数据，怎么用Matlab原生工具链，稳定、可调参、可验证地重建出物理意义明确的波前相位分布？所有代码用matlab2014a及以上版本实测通过，2024a用户无需降级兼容。
2. 从光斑图到相位图：ResUNet重建流程拆解与Matlab实现细节
2.1 Shack-Hartmann传感器建模：为什么必须自己生成仿真数据？
真实Shack-Hartmann传感器输出的是CCD图像上的光斑阵列，每个子孔径对应一个光斑中心坐标。但直接拿实测图像训练网络极不可靠：光斑形变、背景噪声、CCD非线性响应、微透镜阵列畸变都会污染位移估计。因此，该仿真包的第一步是可控建模：在SH_simulation/目录下，generate_SH_pattern.m脚本基于几何光学原理，将输入波前相位（如Zernike组合）映射为理想光斑位移场。
% generate_SH_pattern.m 关键段落 lambda = 632.8e-9; % 波长，单位米 f_lenslet = 5e-3; % 微透镜焦距，单位米 pitch = 0.4e-3; % 微透镜间距，单位米 N_sub = 16; % 子孔径数（16x16网格） % 核心公式：光斑位移 dx = (lambda * f_lenslet / (2*pi)) * dPhi/dx % 实际实现中，用有限差分近似梯度，并叠加高斯白噪声（SNR=25dB可调） dx = -lambda * f_lenslet / (2*pi) * gradient(Phi, 1) + noise_x; dy = -lambda * f_lenslet / (2*pi) * gradient(Phi, 2) + noise_y;提示：这段代码里的
gradient(Phi, 1)计算相位Φ在x方向的偏导，本质是波前斜率，再乘以光学放大系数得到像素位移。noise_x和noise_y是可控的加性噪声，SNR值在config.m中统一配置。这步建模决定了后续网络学的是“物理可解释的映射”，而非黑盒拟合——如果你跳过这步直接喂实测图，网络大概率学到的是传感器硬件缺陷，不是波前本身。
2.2 ResUNet网络结构：Matlab原生实现的轻量级U-Net变体
网络定义在HWDSS-master/SH_resUNet/目录，核心是resunet_architecture.m。它不是简单套用U-Net，而是针对波前重建任务做了三点关键裁剪：
- 输入通道压缩：传统U-Net输入RGB三通道，这里输入是单通道光斑位移图（dx或dy），但网络设计为双通道输入（[dx, dy]），强制网络联合学习两个正交方向的耦合关系；
- 跳跃连接注入物理约束：在每个下采样层后，不直接拼接特征图，而是将对应尺度的Zernike基函数系数预测值（由小分支网络生成）加权融合进跳跃特征，让深层特征自带低阶像差先验；
- 输出层物理归一化：最终输出不是任意矩阵，而是经过
phase_normalization_layer.m处理的相位图，确保均值为0、PV值（峰谷值）落在合理光学范围（±2π rad以内）。
% resunet_architecture.m 中跳跃连接融合逻辑（简化版） function fused_feat = fuse_skip_with_zernike(low_feat, zernike_coeff) % low_feat: 下采样后的特征图，size=[H,W,C] % zernike_coeff: Zernike系数向量，size=[1, N_zernike] zernike_map = zernike_to_phase_map(zernike_coeff, size(low_feat(1:2))); % 将系数转为同尺寸相位图，再经1x1卷积映射到C通道 zernike_proj = conv2d(zernike_map, weights_zernike); fused_feat = low_feat + 0.3 * zernike_proj; % 0.3为可调融合权重 end参数说明：
zernike_to_phase_map调用zernike_poly.m生成标准Zernike多项式网格；weights_zernike是可学习参数，在train_resunet.m中随主网络一同优化。这个设计让网络即使在训练数据不足时，也能退化为“Zernike拟合+残差修正”模式，避免完全失效。
2.3 训练与推理全流程：从config.m到SH_resUNet_demo.m
整个流程由config.m统一管理超参数，这是新手最该盯住的文件。它不像Python项目那样分散在多个yaml里，所有关键开关集中在此：
| 参数名 | 默认值 | 作用说明 | 新手建议 |
|---|---|---|---|
TRAIN_EPOCHS | 150 | 总训练轮数 | 首次运行可设为30快速验证流程 |
BATCH_SIZE | 8 | 每批样本数 | 显存不足时调至4，避免OOM |
LEARNING_RATE | 1e-4 | 初始学习率 | 若loss不降，尝试1e-5 |
NOISE_SNR | 25 | 仿真噪声信噪比(dB) | 实测数据质量差，可降至20 |
ZERNIKE_ORDER | 5 | Zernike展开最高阶数 | 影响zernike_to_phase_map精度 |
执行入口是SH_resUNet_demo.m，它按顺序调用：
load_data.m：读取img/下预生成的.mat文件（含dx,dy,true_phase三字段）；preprocess_data.m：对位移图做零均值化、归一化到[-1,1]，对相位图做PV归一化；train_resunet.m：构建网络、设置优化器（Adam）、启动训练；test_resunet.m：用训练好模型推理，输出recon_phase.mat。
% SH_resUNet_demo.m 中关键调用链 cfg = config(); % 加载全部参数 [data_train, data_val, data_test] = load_data(cfg); % 数据加载 data_train = preprocess_data(data_train, cfg); % 预处理 net = train_resunet(data_train, data_val, cfg); % 训练 results = test_resunet(net, data_test, cfg); % 推理 save('recon_phase.mat', 'results'); % 保存结果逻辑说明：
preprocess_data.m中的归一化不是随意的。位移图归一化到[-1,1]是为了适配tanh激活函数的输出范围；相位图PV归一化（true_phase = (true_phase - mean(true_phase(:))) / max(abs(true_phase(:)))）则保证网络输出能直接反推物理相位值。若跳过此步，网络会因量纲混乱而无法收敛。
3. 数据准备与格式规范：如何用自己的光斑图替换案例数据
3.1 案例数据结构解析：img/目录下的.mat文件到底装了什么？
打开img/SH_data_001.mat，你会看到三个变量：
dx: 大小为[N_sub, N_sub]的double矩阵，每个元素是第(i,j)个子孔径的x方向位移（单位：像素）；dy: 同尺寸y方向位移矩阵；true_phase: 大小为[256, 256]的double矩阵，是生成该位移数据所用的真实波前相位（单位：弧度）。
注意：
dx和dy是整数坐标位移，不是亚像素精度。这意味着传感器CCD分辨率已隐含在pitch和f_lenslet参数中。若你的实测数据是亚像素坐标（如[12.34, 56.78]），需先用round()取整，否则网络输入维度错乱。
3.2 替换自有数据的四步操作法
假设你有一组实测光斑图，存为my_data/spot_001.tif,my_data/spot_002.tif...，需转换为包要求的.mat格式：
步骤1：批量提取光斑中心用img/extract_spots.m脚本（已提供）处理TIFF序列：
% extract_spots.m 示例 folder = 'my_data/'; tif_files = dir(fullfile(folder, '*.tif')); for i = 1:length(tif_files) img = imread(fullfile(folder, tif_files(i).name)); [dx_i, dy_i] = find_spot_centers(img, N_sub); % 内置函数，返回16x16位移矩阵 % 保存为.mat save(fullfile('my_img/', ['spot_' num2str(i, '%03d') '.mat']), 'dx_i', 'dy_i'); end关键点：
find_spot_centers.m使用高斯拟合定位光斑，比质心法抗噪更强。若你的光斑严重重叠，需手动调整其内部gaussian_fit_sigma参数（默认1.5像素）。
步骤2：生成伪真值相位（仅用于验证）实测数据没有true_phase，但评估重建质量需要参考。此时启用config.m中的GENERATE_PSEUDO_TRUE_PHASE = true，系统会用Zernike拟合你的dx_i/dy_i，生成一个物理自洽的相位作为伪真值：
% 在train_resunet.m中，当无true_phase时触发 if ~isfield(data, 'true_phase') pseudo_phase = zernike_fit_from_shifts(dx_i, dy_i, cfg); data.true_phase = pseudo_phase; end步骤3：组织数据目录将生成的.mat文件放入新文件夹my_img/，并在config.m中修改：
cfg.DATA_PATH = 'my_img/'; % 原为'img/' cfg.TRAIN_RATIO = 0.7; % 训练集占比，实测数据少时可设为0.9步骤4：验证数据加载运行test_data_loading.m（包内提供），它会加载第一个文件并显示dx,dy,true_phase的热力图，确认尺寸和数值范围正确：
% test_data_loading.m 输出示例 disp(['dx size: ', num2str(size(dx))]); % 应显示 16 16 disp(['dx range: [', num2str(min(dx(:))), ', ', num2str(max(dx(:))), ']']); % 合理范围[-5,5]像素避坑提醒：若
dx矩阵出现全零或全NaN，大概率是find_spot_centers.m未能检测到光斑——检查TIFF是否为8位灰度图（非RGB），或调整extract_spots.m中的threshold参数（默认0.3*max_intensity）。
4. 避坑指南：五个让新手当场翻车的Matlab深度学习陷阱
4.1 现象：train_resunet.m报错“未定义函数或变量 'dlnetwork'”
原因：dlnetwork是Matlab R2018a引入的深度学习网络类，而包声明支持matlab2014a。实际测试发现，2014a/2016a用户必须使用旧版网络定义方式（network对象），但包内代码默认走dlnetwork路径。
解决：打开config.m，将cfg.USE_DLN_NETWORK = true改为false，并确保HWDSS-master/legacy_network/目录在Matlab路径中。旧版网络不支持自动微分，需手动实现backward函数，但包已提供legacy_resunet_backward.m。
4.2 现象：训练loss曲线剧烈震荡，100轮后仍高于0.15
原因：dx/dy数据未归一化。原始位移值可能达±20像素，而网络权重初始化默认适配[-1,1]输入，导致梯度爆炸。
解决：在preprocess_data.m中，确认以下代码未被注释：
dx = (dx - mean(dx(:))) / std(dx(:)); % 必须做Z-score归一化 dy = (dy - mean(dy(:))) / std(dy(:));血泪经验：曾有A同学跳过此步，折腾三天以为网络结构有问题，最后发现只是归一化漏了一行。
4.3 现象：SH_resUNet_demo.m运行到test_resunet时内存溢出（Out of Memory）
原因：BATCH_SIZE设得过大，且N_sub=16时，单个batch需加载8张16x16位移图+对应相位图，显存占用超限。
解决：在config.m中将BATCH_SIZE从8改为2，并在test_resunet.m开头添加：
% 强制使用CPU推理，避免GPU显存争抢 net = dlnetwork(net.Learnables, 'ExecutionEnvironment', 'cpu');4.4 现象：重建相位图全是平滑渐变，丢失高频细节（如彗差、球差）
原因：ZERNIKE_ORDER参数过低（如设为3），导致网络只能拟合低阶像差，高频信息被当作噪声丢弃。
解决：将config.m中ZERNIKE_ORDER从默认5提升至7或8，并同步增大zernike_to_phase_map的网格分辨率（修改zernike_poly.m中N_grid = 256为512）。注意：阶数每+1，训练时间约+40%，需权衡。
4.5 现象：recon_phase.mat中results.recon_phase尺寸为[1,1]，不是预期[256,256]
原因：test_resunet.m中resize_image函数调用错误。包内utils/resize_image.m是Matlab原生imresize的封装，但若用户Matlab版本低于2017b，imresize不支持'bicubic'插值。
解决：打开utils/resize_image.m，将第12行：
out = imresize(in, scale, 'bicubic');改为：
out = imresize(in, scale, 'bilinear'); % 兼容所有版本5. 重建结果验证：三重交叉验证法确保物理可信度
5.1 相位残差分析：不只是看MSE，要看Zernike分解谱
单纯计算mean((true_phase - recon_phase).^2)（MSE）会掩盖系统性偏差。更可靠的做法是：将重建相位recon_phase和真值true_phase分别做Zernike分解，对比各阶系数误差。包内validate_zernike_fit.m提供一键分析：
% validate_zernike_fit.m 核心逻辑 coeff_true = zernike_decompose(true_phase, 7); % 分解到7阶 coeff_recon = zernike_decompose(recon_phase, 7); error_coeff = coeff_true - coeff_recon; % 绘制误差柱状图，重点关注低阶（活塞、倾斜、离焦）和高阶（球差、三叶草） bar([1:length(error_coeff)], abs(error_coeff)); xlabel('Zernike Term Index'); ylabel('Coefficient Error (rad)'); title('Zernike Coefficient Reconstruction Error');参数说明：
zernike_decompose使用最小二乘法拟合，coeff_true(1)是活塞项（应接近0，因相位均值已归零），coeff_true(2:3)是x/y倾斜，coeff_true(4)是离焦——这些低阶项误差应<0.05 rad，否则重建失效。
5.2 PSF重构验证：用重建相位反推点扩散函数
波前重建的终极检验是能否复现光学系统的点扩散函数（PSF）。包内compute_psf_from_phase.m实现菲涅尔衍射计算：
% compute_psf_from_phase.m 关键步骤 lambda = 632.8e-9; k = 2*pi/lambda; % 对重建相位做FFT，乘以传播因子，再IFFT psf_recon = ifft2(fft2(exp(1i * recon_phase)) .* propagation_kernel); psf_recon = abs(psf_recon).^2; % 强度分布 % 与真值PSF（由true_phase计算）对比 psf_true = compute_psf_from_phase(true_phase, lambda); ssim_value = ssim(psf_recon, psf_true); % 结构相似性指数逻辑说明：
propagation_kernel是菲涅尔衍射的传递函数，由config.m中PROPAGATION_DISTANCE（传播距离）和PIXEL_SIZE（CCD像素尺寸）决定。若ssim_value < 0.85，说明重建相位的高频信息严重失真，需检查ZERNIKE_ORDER或增加训练轮数。
5.3 实验一致性验证：用同一组dx/dy输入，对比不同算法输出
包内compare_algorithms.m提供三算法并行验证：
- Zernike拟合：传统方法，调用
zernike_fit_from_shifts.m - SVD拟合：基于奇异值分解的线性方法，调用
svd_wavefront_recon.m - ResUNet：本包深度学习方法
% compare_algorithms.m 输出表格 algorithms = {'Zernike', 'SVD', 'ResUNet'}; mse_list = zeros(1,3); for i = 1:3 phase_i = eval([algorithms{i}, '_recon(dx, dy, cfg)']); mse_list(i) = mean((true_phase - phase_i).^2); end T = table(algorithms', mse_list', 'VariableNames', {'Algorithm', 'MSE'}); disp(T); % 示例输出： % Algorithm MSE % _________ ________ % 'Zernike' 0.1245 % 'SVD' 0.0892 % 'ResUNet' 0.0327 ← 深度学习优势在此体现技巧：从那以后我每次交付波前重建结果，都强制跑一遍
compare_algorithms.m，把三行MSE值写进报告附录。导师一眼就能看出深度学习是否真带来了提升，而不是在炫技。如果ResUNet的MSE不比Zernike低30%以上，我就知道该回去调参了——不是模型不行，是参数没喂对。希望帮到你。
本文还有配套的精品资源，点击获取
