简介:基于深度学习的波前重建系统是一套面向深度学习、机器学习课程设计及毕业设计的完整工程资源,聚焦利用CNN与残差U-Net网络对波前畸变数据进行自动特征提取与重建,弥补传统Zernike拟合、相位展开等算法依赖先验知识、抗噪性弱的不足,适合希望将人工智能技术落地到光学信号处理方向的研究者参考。压缩包共46个文件,约1.83MB,涵盖Python与MATLAB两类源码、模型权重/输出文件、PNG结果图、配置文件及说明文档,其中py文件实现网络训练与相位恢复,m文件用于Zernike多项式及夏克-哈特曼波前传感器模拟,png文件展示光场成像与网络结构等可视化结果。已有41人学习下载。资源内SH_simulation文件夹可生成波前畸变模拟数据,解决真实数据难获取的问题;SH_resUNet_demo文件夹则给出残差U-Net的完整demo,配套README和说明文档,便于快速理解项目架构、运行流程与复现实验,是完成课程设计或毕业设计的高质量参考。
1. 基于深度学习的波前重建系统:为什么它值得你完整跑一遍
波前重建在自适应光学、光场成像和眼底视网膜成像里都是核心环节——把夏克-哈特曼传感器测到的波前斜率还原成完整相位。传统做法是 Zernike 拟合或区域法积分,依赖先验且对噪声敏感;而深度学习用大量仿真样本端到端学一个映射,鲁棒性明显更好。这个zip里既有 MATLAB 的波前仿真与传感器响应模拟,又有 PyTorch 的残差 U-Net 训练、测试与相位恢复脚本,正好构成一条「仿真出数据 → 网络重建 → 恢复相位」的完整链路。如果你是做深度学习的期末大作业或本科毕设,这套代码能直接跑通并给你一个可展示的实验闭环。
2. SH_simulation 仿真模块:用 Zernike 多项式造出能训练的数据集
深度学习重建波前的第一道门槛不是网络,而是数据。真实的夏克-哈特曼波前传感器数据很难拿到,尤其是带精确相位标签的,所以这个项目选择用 MATLAB 仿真生成训练集——这也是工业界和学术界通用的做法。SH_simulation 文件夹里那一串 zernfun.m、zernikeN.m、RandomWaveFun.m、HartmanDiff.m,就是把「随机像差 → Zernike 相位 → 传感器响应」这条链路完整实现了一遍。
2.1 为什么选 Zernike 多项式而不是直接随机生成相位
Zernike 多项式是定义在单位圆上的正交多项式系,每个项对应一种光学像差形态:离焦、像散、彗差、球差、三叶草等等。波前畸变可以被分解为这些项的线性叠加,系数就代表各像差成分的权重。用 Zernike 表示波前的最大好处是物理意义清晰,且每一项在圆域内正交,系数之间没有冗余,非常便于控制生成数据的分布。
% zernfun.m 的核心调用:在归一化极坐标下计算单个 Zernike 项 x = linspace(-1, 1, 256); [X, Y] = meshgrid(x); [theta, rho] = cart2pol(X, Y); rho(rho > 1) = 1; % 光瞳外半径做掩膜,避免圆外噪声污染 Z = zernfun(n, m, rho, theta); % n 为径向阶数,m 为角向频率这段代码的关键是cart2pol把直角坐标转成极坐标,zernfun接收归一化半径rho和角度theta,返回与 X/Y 同尺寸的相位面。注意rho(rho > 1) = 1这一步:Zernike 多项式定义在单位圆内,圆外区域如果不截断,边缘会出现无意义的高频数值,后续无论是算梯度还是做网络输入都会引入脏数据。
zernikeN.m则是把多个 Zernike 项叠加成一张完整的波前相位图。常见做法是先生成前 N 项的系数数组,然后循环累加:
% zernikeN.m 的叠加逻辑示意 coeffs = zernikeN(15); % 生成前 15 项 Zernike 系数,长度 15 phase = zeros(size(rho)); for k = 1:size(coeffs, 1) phase = phase + coeffs(k) * zernfun(n_list(k), m_list(k), rho, theta); end这里n_list和m_list对应 Noll 编号顺序下每一项的径向阶数与角向频率,项目里应该有完整的编号表。为什么叠加而不是直接生成随机相位面?因为直接随机生成的相位图往往高频噪声占主导,不符合真实光学系统「低阶像差能量大、高阶像差能量小」的物理规律。用 Zernike 叠加,你就能通过控制系数分布来控制数据集的物理合理性。
2.2 RandomWaveFun.m 的随机系数生成:让数据集贴近真实像差分布
RandomWaveFun.m 的作用是生成一组随机的 Zernike 系数,用来构造一张随机的波前畸变。这里的随机不是纯均匀随机,而是要让低阶像差异常项(离焦、像散)系数幅度大,高阶项幅度小,模拟大气湍流或光学系统像差的能量分布特征。
% RandomWaveFun.m 的系数生成策略(常见做法示意) num_terms = 15; coeffs = zeros(1, num_terms); for k = 1:num_terms % 低阶项系数方差大,高阶项按阶数衰减 coeffs(k) = randn * exp(-(k-1) / 6); end指数衰减因子exp(-(k-1)/6)是这里的关键参数:前几项系数有较大概率取到明显非零的值,第 10 项以后的系数基本被压制到接近零。这个衰减速率直接决定数据集的像差复杂度——衰减越快,数据越集中在低频像差,网络训练更容易收敛,但泛化到剧烈畸变场景时能力会弱;衰减越慢,数据越多样,训练难度也越大。实际调参时可以先快后慢,先在简单分布上验证网络能收敛,再逐步放开。
生成系数后,再用 zernfun 累加得到相位图,这一步在 zernikeWaveFun.m 里完成。zernikeWaveFun.m 与 zernikeN.m 的区别在于:前者是给一组具体系数生成对应波前,后者是随机生成系数并返回系数和波前两个输出。建议你训练前先跑一遍 zernikeWavetest.m,把生成的波前图用 imwriteTest.m 存成图片看看分布是否合理。
2.3 HartmanDiff.m 与传感器模拟:把相位转成梯度场
深度学习模型的输入不是相位图本身,而是夏克-哈特曼传感器实际能测到的量——波前斜率。HartmanDiff.m 就是把连续的 Zernike 相位场离散成子孔径网格,模拟传感器每个微透镜下的光斑偏移。
% HartmanDiff.m 的梯度场计算逻辑示意 % 输入 phase: 256x256 的相位图;输出 grad_x, grad_y: 子孔径网格上的斜率 subap_num = 8; % 子孔径数,8x8 网格 block = size(phase, 1) / subap_num; for i = 1:subap_num for j = 1:subap_num % 取当前子孔径区域 patch = phase((i-1)*block+1 : i*block, (j-1)*block+1 : j*block); % 子孔径内 x/y 方向的平均斜率 [gx, gy] = gradient(patch); grad_x(i, j) = mean(gx(:)); grad_y(i, j) = mean(gy(:)); end end这段代码的逻辑很直观:把相位图切成 8×8 的小块,每块内用gradient求局部斜率再取平均,得到该子孔径对应的波前斜率。物理上,这个值正比于传感器上光斑相对参考位置的偏移量。子孔径数subap_num是个重要参数:设得越大,空间采样越密,重建的细节越多,但每个子孔径内的光斑信噪比下降,噪声变大;设得越小,采样越稀疏,网络输入信息量不足。
HartmanWavefront.m 和 WavefrontGrad.m 在这个流程里分别承担从子孔径斜率恢复波前和从完整相位求梯度的角色。PositionFunc.m 和 NormalCoord.m 负责生成子孔径的位置坐标与归一化坐标,是连接相位图和传感器响应的坐标系基础。
到这里,SH_simulation 就把「随机像差 → 相位 → 梯度」这条数据生成链路打通了。一个训练样本就是一对数据:输入是梯度场(传感器能测到的),标签是原始相位图(要重建的目标)。后面深度学习的全部工作,本质上就是让网络学会这个映射。
3. 从梯度到重建:传统算法的定位与深度学习数据管线的搭建
在进入 U-Net 之前,需要先把传统重建方法交代清楚。原因有两点:一是这个项目里保留了完整的传统重建 MATLAB 代码,你不理解它们就无法解释深度学习相对它们的优势;二是传统重建结果可以作为深度学习的参照基准,用来量化网络重建质量提升了多少。
3.1 WavefrontGrad.m 与 FLHFtest.m:传统区域法重建是怎么算的
波前重建的传统思路是:已知子孔径上的斜率测量值,求整个波前相位。这是典型的反问题,核心在于建立差分矩阵 H,使 H·φ = s,其中 φ 是待求的相位向量,s 是斜率向量。WavefrontGrad.m 负责从相位图求梯度,而 FLHFtest.m、FNRdiff.m、FNRSFFTdiff.m 这些脚本实现的是不同的重建求解器。
区域法的思想是逐点建立相位点与邻域斜率之间的关系。对于南极点(Southwell)模型,相位点与斜率网格交错排列,每个斜率测量值连接两个相邻相位点,于是有 φ(i+1,j) - φ(i,j) = s_x(i,j) × d,其中 d 是子孔径间距。把所有这些方程堆起来,得到稀疏矩阵 H,然后用最小二乘解 φ。
% FLHFtest.m 的最小二乘重建核心(逻辑示意) % H 为稀疏差分矩阵,s 为斜率向量,phi 为待重建相位 lambda = 1e-3; % 正则化系数 phi = (H' * H + lambda * speye(size(H,2))) \ (H' * s);lambda是 Tikhonov 正则化系数,它存在的意义是处理 H^T·H 可能奇异的问题。工程上这个参数很敏感:设太小,噪声会被放大,重建面出现高频抖动;设太大,重建面被过度平滑,细节丢失。传统方法的痛点就在这里——lambda 需要针对不同数据集反复调,而且对噪声鲁棒性有限。FNRSFFTdiff.m 走的是另一条路:用 FFT 在频域求解泊松方程,速度快但边界处理麻烦,重建面边缘经常出现振铃。
这些脚本项目里都有可运行的测试(FlhfDifftest2.m、FNRdiff.m 等),建议先跑一遍,保存传统重建的结果,后续和 U-Net 的输出做对比——这是毕设答辩时很有力的实验数据。
3.2 dataset.py:把 MATLAB 仿真结果封装成 PyTorch 能吃的样本
MATLAB 负责产出数据,Python 负责训练网络。中间桥梁是 dataset.py。项目里的 UobPhase_test.mat 就是一组仿真或实测的相位数据,dataset.py 的核心任务是把 .mat 文件里的梯度场和相位图封装成 PyTorch 的 Dataset 对象。
# dataset.py 核心逻辑:把仿真好的 .mat 数据包装成可迭代的数据集 import scipy.io as sio import torch from torch.utils.data import Dataset class WavefrontDataset(Dataset): def __init__(self, mat_path, norm=True): # MATLAB 保存的梯度场和相位图 data = sio.loadmat(mat_path) self.grad = data['grad'].astype('float32') # 形状 (N, H, W) 或 (N, 2, H, W) self.phase = data['phase'].astype('float32') # 形状 (N, H, W) # 归一化到 [0,1],避免数值尺度差异压过损失 if norm: for i in range(self.grad.shape[0]): self.grad[i] = (self.grad[i] - self.grad[i].min()) / \ (self.grad[i].max() - self.grad[i].min() + 1e-8) self.phase[i] = (self.phase[i] - self.phase[i].min()) / \ (self.phase[i].max() - self.phase[i].min() + 1e-8) def __len__(self): return len(self.grad) def __getitem__(self, idx): # 输入是梯度场,标签是相位图 x = torch.from_numpy(self.grad[idx]) y = torch.from_numpy(self.phase[idx]) return x, ynorm=True这个开关值得注意。如果数据集中样本的相位幅度差异很大,全局归一化会让小相位样本的梯度变得极小,网络难以学到有效特征;逐样本归一化则保证了每个样本的数值范围一致,但也意味着网络学到的是「相对分布」而不是「绝对幅度」。我在实际项目中倾向逐样本归一化,因为波前重建关心的是相位分布形态,绝对值可以在恢复阶段通过额外刻度还原。
需要特别注意的是scipy.io.loadmat的版本限制。MATLAB 保存 .mat 时如果用了 v7.3 格式(HDF5),scipy 会直接报错,解决办法是在 MATLAB 里用save('data.mat', 'grad', 'phase', '-v7')强制存为旧版,或者在 Python 端改用h5py读取。这个坑后面避坑章节还会展开。
到这里,数据管线就通了:MATLAB 仿真出带标签的梯度-相位对,dataset.py 把它们封装成 PyTorch Dataset,下一步就是搭建网络并训练。
4. SH_resUNet_demo:残差 U-Net 的工程实现与训练闭环
SH_resUNet_demo 是这套资源的核心部分。它的目标很明确:输入梯度场,输出重建的相位图。之所以选择残差 U-Net 而不是普通 CNN,是因为波前重建本质上是像素级的回归任务,输入输出是逐像素对齐的,这正好是 U-Net 这类编码-解码结构的强项;加上残差连接,深层网络的梯度传导更稳定,训练收敛更快。
4.1 EasyNet / PaperNet / TestNet:三个网络脚本各干什么
项目里有三个网络相关脚本,初学者常常搞混它们的定位。EasyNet.py 是一个简化的 CNN 基线模型,用于快速验证数据管线是否通畅;PaperNet.py 是完整的残差 U-Net 实现,也就是最终用于波前重建的主力网络;TestNet.py 是推理脚本,加载训练好的权重对测试数据做重建并评估指标。
残差 U-Net 的结构可以拆成三块:编码器、解码器、跳跃连接。编码器逐级下采样提取多尺度特征,解码器逐级上采样恢复空间分辨率,跳跃连接把编码器对应层的特征拼接到解码器,保住了细节信息。残差块在这里的典型实现是:
# EasyNet.py 中的残差块实现(示意) import torch.nn as nn class ResBlock(nn.Module): def __init__(self, in_ch, out_ch): super().__init__() self.conv1 = nn.Conv2d(in_ch, out_ch, 3, padding=1) self.bn1 = nn.BatchNorm2d(out_ch) self.conv2 = nn.Conv2d(out_ch, out_ch, 3, padding=1) self.bn2 = nn.BatchNorm2d(out_ch) self.shortcut = nn.Conv2d(in_ch, out_ch, 1) if in_ch != out_ch else nn.Identity() def forward(self, x): out = self.bn1(self.conv1(x)) out = self.bn2(self.conv2(out)) return nn.ReLU()(out + self.shortcut(x))残差块的思路是让网络学习「输入和输出之间的差异」而不是直接学习映射本身,这样即使网络很深,梯度也能通过恒等映射的捷径顺畅回传。shortcut分支在通道数不一致时用 1×1 卷积对齐,这是残差网络的标准做法。编码器和解码器之间,下采样用步长为 2 的卷积或最大池化,上采样用转置卷积或双线性插值加卷积,两种方案在项目里都常见。转置卷积能学到的上采样参数更多,但容易产生棋盘伪影;双线性插值稳但表达能力弱一些。对于波前重建这种平滑连续的任务,我建议用双线性插值,重建面更干净。
4.2 main.py 训练流程与 train_logs 曲线怎么看
main.py 把完整训练流程串起来:加载数据集、实例化网络、定义损失函数和优化器、循环训练、记录 loss、周期性保存模型。参数设置上,这套代码的默认值基本合理但有几个地方需要根据你的数据规模调整:
| 参数 | 典型值 | 调整依据 |
|---|---|---|
| batch_size | 16 | 显存不足时降到 4~8 |
| learning_rate | 1e-3 | 用 Adam 时可以偏大,SGD 需降到 1e-4 |
| epochs | 100 | 看 val_loss 是否还在明显下降 |
| 损失函数 | MSE | 波前重建是回归任务,MSE 最直接 |
| 优化器 | Adam | 默认 β1=0.9, β2=0.999 |
损失函数用 MSE 是因为相位图是连续值回归,逐像素算平方误差最直观。但 MSE 对低频大尺度误差和高频细节误差一视同仁,这会导致重建结果整体形状接近但细节模糊。后续进阶改造时会提到换成 Zernike 系数回归或混合损失。
训练过程中 main.py 会把每个 epoch 的 loss 记录下来写到 train_logs 文件夹。loss 曲线的判读有几个经验:前 10 个 epoch 如果 loss 剧烈震荡且均值不下降,先检查数据归一化和学习率,而不是急着改网络结构;如果训练 loss 持续下降但验证 loss 在第 30 个 epoch 左右开始回升,就是过拟合了,此时优先加数据量或数据增强,其次考虑加 dropout。
4.3 Phase_recovery.py:从网络输出到可评估的相位
训练完成只是第一步,网络的裸输出还不能直接用——因为数据在 dataset.py 里做了逐样本归一化,输出值都在 [0,1] 区间,需要反归一化还原为真实的相位数值。Phase_recovery.py 干的就是这件事。它还做了一件事:去掉 piston 项,也就是把重建相位的整体均值归零。
# Phase_recovery.py 的相位恢复逻辑(示意) import numpy as np def recover_phase(net_out, orig_min, orig_max): # 反归一化:还原到原始相位数值范围 phase = net_out * (orig_max - orig_min) + orig_min # 去掉 piston:整体均值置零,只保留相对畸变 phase = phase - np.mean(phase) return phaseorig_min和orig_max是 dataset.py 归一化时记录的该样本原始最小值和最大值。这里要强调的是:piston 项(整体平移)在波前斜率测量里本身不可见,因为传感器测的是相对斜率而非绝对相位,所以重建时去掉均值是标准操作。如果不去,整体抬升会让 PV 值虚高,导致误判。
TestNet.py 加载训练好的权重,对 UobPhase_test.mat 里的样本做重建,计算 RMSE 和 PV 值。RMSE 是逐像素均方根误差,反映整体重建精度;PV 是峰谷差,反映相位面最差点的偏差。这两个指标必须同时看——RMSE 很小但 PV 很大,说明误差集中在个别点;PV 小但 RMSE 大,说明整体散布着低频偏差。
5. 避坑排查:从 MATLAB 到 PyTorch 的五处翻车点
这套代码我在类似场景下踩过不少坑,挑五个最容易翻车的记录在这里。每条都是「现象 → 原因 → 解决」的结构,你在复现时如果遇到同样现象可以直接对应处理。
5.1 数据侧的翻车点:mat 版本、尺度与样本重复
坑一:loadmat 直接报错。
现象:dataset.py 里sio.loadmat('UobPhase_test.mat')抛出NotImplementedError: Please use HDF reader for matlab v7.3 files。
原因:MATLAB 新版默认的 save 格式是 v7.3,底层是 HDF5,scipy.io 不支持。
解决:两种方案任选。一是去 MATLAB 里重新保存:save('UobPhase_test.mat', 'phase', 'grad', '-v7'),强制用老版本格式;二是在 Python 端改用h5py读取,但要手动处理 key 的结构差异。我建议用方案一,省事且兼容性好。
坑二:loss 不降,梯度爆炸。
现象:训练第一个 epoch loss 就是 NaN,或者 loss 在某个值附近震荡完全不下行。
原因:相位图数值范围可能是 0~10,而梯度场经过归一化后是 0~1,两者量级差太多,MSE 损失对相位误差的梯度被放大,导致优化不稳定。
解决:统一做标准化而不是归一化,即把每个样本的均值减掉、标准差除了,让数据分布在零附近。相位和梯度都要做,且记录均值和标准差用于恢复阶段还原。
坑三:验证集 loss 低得离谱,换真实数据就崩。
现象:训练集和验证集 loss 都很好看,但部署到新数据上重建质量急剧下降。
原因:RandomWaveFun.m 的随机系数生成如果固定了随机种子,多次调用会生成高度相似甚至相同的样本,导致验证集和训练集分布重叠严重,模型实际是记住了样本而非学到了映射。
解决:划分数据集前先对样本做去重,或者给每个样本计算一个特征哈希值,按哈希值而不是随机比例划分。最简单的方式是用生成的随机种子号作为文件名的一部分,让每个样本有唯一标识,再按标识划分。
5.2 训练侧的翻车点:几何失配与显存瓶颈
坑四:训练集和测试集用的子孔径网格数不一致。
现象:网络在仿真数据上 RMSE 不错,但把 HartmanDiff.m 里的subap_num从 8 改成 16 再生成数据,直接用已训练模型测试,重建结果一片混乱。
原因:输入的梯度场空间分辨率变了,网络输入尺寸和特征分布都不匹配。8×8 子孔径的梯度图与 16×16 的梯度图,信息密度完全不同。
解决:训练时就固定一个子孔径数,或者在训练阶段做数据增强——随机在 6 到 16 之间选择子孔径数生成训练样本,让网络学会适应不同采样密度。后者效果更好,网络泛化能力显著提升。
坑五:显存 OOM。
现象:main.py 训练到一半报CUDA out of memory。
原因:高分辨率相位图(512×512 以上)加上深层的 U-Net,中间特征图数量巨大,batch_size=16 很容易塞爆显存。
解决:优先把 batch_size 降到 4 或 2;其次把输入图裁剪成 patch 训练,例如从 256×256 的图中随机裁 128×128 的小块,相当于隐式的数据增强;再不行就用混合精度训练,PyTorch 里一句torch.cuda.amp.autocast()就能把显存占用砍半,代价是精度轻微下降,对波前重建这种任务影响可以忽略。
6. 把 demo 改成自己的毕设:三个能立刻上手的改造
现有代码能跑通,但离一篇合格的毕业设计或期末大作业还差一些让你自己的贡献点。下面是三个不需要重构整体架构就能落地的改造方向,按性价比从高到低排。
第一个改造:把输出从相位图换成 Zernike 系数回归。现网络的输出是像素级相位图,容易在高频区域产生无物理意义的起伏。改动方式是保持 U-Net 编码器不变,把解码器最后一层替换成全局平均池化加全连接层,输出长度等于 Zernike 项数。损失函数从 MSE 相位图误差换成 Zernike 系数误差。这样重建出的相位天然平滑,且直接得到了像差分解结果——答辩时这个输出比一张灰度图更有说服力。
第二个改造:数据管线换成 h5py 加预缓存。当前 dataset.py 每个 epoch 都在重复读 .mat 文件,训练数据量大时 I/O 会成为瓶颈。先把所有训练数据打包进一个 .h5 文件,然后在 Dataset 初始化时一次性读入内存,训练速度能提升一个档次。这在大作业报告里可以写成一节「数据管线优化」,工作量不大但显得很工程化。
第三个改造:评估指标加上斜率残差和 PV 值。现在 TestNet.py 只算 RMSE,这个指标对低频偏差不敏感。建议在测试脚本里加上两项:一是把重建相位重新送入 WavefrontGrad.m 求梯度,与网络输入的梯度场对比,算斜率残差的 RMS;二是记录重建相位的 PV 值与真实相位 PV 值的差。前者直接对应传感器实测数据的一致性,后者对应光学系统评价指标。这三项指标合在一起,你的实验结果部分会扎实很多。
从那以后我做波前重建实验,都会强制走一遍先仿真出数据分布、再对比传统重建基线、最后看多指标评估的流程。这套代码帮我省掉了大量重复造轮子的时间,也让我在对付 Zernike 项数选择、子孔径密度这类细节时有了可调的依据。希望你也能把它跑通,然后改出自己的版本。希望帮到你。
本文还有配套的精品资源,点击获取