简介:这份资源面向光学遥感、天文观测及图像处理方向的研究生与工程师,提供大气湍流退化图像复原的硕士论文及配套MATLAB仿真代码,帮助读者理解湍流对成像质量的影响机理并掌握复原算法的工程实现。压缩包共25个文件,约20.72MB,包含1篇PDF论文、3个m脚本、20张bmp实验图像及1个asv备份文件;其中论文系统阐述瑞利散射、折射指数不均匀性与像差效应,脚本则覆盖图像采集、噪声添加、相位恢复、去模糊与重构等关键环节,bmp图像用于直观对比不同复原策略的效果。目前已有2699人学习下载,读者可借助代码修改参数、复现卡尔曼滤波自适应光学、光强波动统计模型与傅立叶变换域反卷积等方法的仿真结果,进而对比算法优缺点与适用场景,为相关课题研究或工程开发积累可复用的实现思路与调试经验。
1. 从一张糊掉的远距离照片说起:大气湍流退化图像复原到底在解决什么
拍过远距离目标的人都有体会:明明镜头对焦准确,画面却像隔着一层晃动的水。这不是镜头脏了,也不是对焦失败,而是光在到达传感器之前穿过了折射率随机起伏的大气。温度梯度让空气密度不均,密度不均让折射率不均,波前被随机相位扰动,成像系统的点扩散函数(PSF)就不再是一个理想冲激,而是一个随时间和空间变化的模糊核。这类退化在长焦监控、天文观测、遥感成像里非常普遍。
大气湍流退化图像复原要做的,就是从观测到的退化图像反推清晰图像。它和普通去模糊的区别在于:模糊核不是固定的,而是由湍流统计特性决定的;噪声往往同时存在;单帧图像信息不足,所以工程上常走多帧配准加融合,或者单帧盲去卷积两条路。这篇内容面向做图像处理、遥感、光电工程的从业者,也适合正在用 MATLAB 做课程设计或论文复现的人。下面从退化模型讲起,一路落到可运行的 MATLAB 代码、参数设置和排错。
2. 大气湍流退化模型与 MATLAB 建模:从 PSF 到退化图像
2.1 湍流退化的数学描述与常见 PSF 选择
大气湍流对成像的影响,经典描述来自 Kolmogorov 湍流理论。折射率结构常数 Cn² 决定湍流强度,相干直径 r0 决定横向相干尺度,比值 D/r0(D 为孔径直径)决定退化程度。在长曝光条件下,湍流引起的平均 PSF 常近似为高斯函数;在短曝光条件下,则更接近由相位结构函数导出的复杂核。工程复现里,最常用的两种建模方式是高斯 PSF 和基于 Zernike 多项式的相位屏法。
高斯 PSF 简单、参数少,适合快速验证复原算法;相位屏法更贴近物理,能生成随帧变化的退化,适合做多帧复原的仿真数据。选哪种取决于你的目标:如果论文重点是复原算法本身,用高斯 PSF 造数据就够了;如果重点是湍流建模,就得用相位屏。
| 建模方式 | 核心参数 | 适用场景 | 计算量 |
|---|---|---|---|
| 高斯 PSF | 核尺寸、sigma | 算法快速验证 | 低 |
| 相位屏法 | r0、D、Zernike 阶数 | 物理仿真、多帧 | 中高 |
| 实测标定 PSF | 标定点扩散函数 | 有先验的复原 | 取决于标定 |
2.2 用 MATLAB 生成高斯型湍流退化图像
下面这段代码生成一个高斯 PSF,对清晰图像做卷积,再加高斯噪声,得到退化图像。它是后续所有复原实验的输入。
% 读取清晰图像并转灰度、归一化 img = imread('lena.png'); if size(img,3) == 3 img = rgb2gray(img); end img = im2double(img); % 生成高斯 PSF,ksize 为核尺寸,sigma 控制模糊强度 ksize = 21; sigma = 3.5; [X, Y] = meshgrid(-(ksize-1)/2:(ksize-1)/2, -(ksize-1)/2:(ksize-1)/2); psf = exp(-(X.^2 + Y.^2) / (2*sigma^2)); psf = psf / sum(psf(:)); % 能量归一化,保证亮度不漂移 % 卷积退化,'symmetric' 缓解边界振铃 blurred = imfilter(img, psf, 'conv', 'symmetric'); % 加高斯噪声,noise_var 为噪声方差 noise_var = 1e-4; blurred_noisy = imnoise(blurred, 'gaussian', 0, noise_var); figure; subplot(1,3,1); imshow(img); title('清晰图像'); subplot(1,3,2); imshow(blurred); title('湍流模糊'); subplot(1,3,3); imshow(blurred_noisy); title('模糊加噪声');逻辑说明:先归一化避免量纲问题;PSF 求和归一化是关键,否则卷积后整体亮度会随 sigma 变化。imfilter用conv而非默认相关,符合卷积退化模型;symmetric边界比零填充更少引入暗边。噪声方差 1e-4 对应约 20dB 信噪比,是湍流图像里比较典型的量级。
参数说明:ksize一般取 6 倍 sigma 以上,太小会截断 PSF 导致复原出现振铃;sigma越大模糊越重,D/r0 越大时对应 sigma 越大。做消融实验时固定 ksize,只调 sigma,才能公平比较算法。
2.3 相位屏法生成随帧变化的退化
多帧复原需要每帧不同的 PSF。相位屏法通过生成随机相位并传播到成像面得到 PSF,核心是相位结构函数要匹配湍流统计。
% 简化相位屏:用低阶 Zernike 近似,生成 N 帧不同 PSF N = 10; % 帧数 D = 0.2; % 孔径直径 r0 = 0.05; % 相干直径,D/r0=4,中等湍流 psfStack = zeros(ksize, ksize, N); for k = 1:N % 用随机相位经傅里叶变换近似相位屏 phase = randn(256) * (D/r0)^(5/6); pupil = zeros(256); [xx, yy] = meshgrid(-128:127, -128:127); pupil(xx.^2 + yy.^2 <= 128^2) = 1; % 圆形孔径 field = pupil .* exp(1i * phase); psfFull = abs(fftshift(fft2(field))).^2; % 裁剪并归一化到 ksize c = 128; psfStack(:,:,k) = psfFull(c-10:c+10, c-10:c+10); psfStack(:,:,k) = psfStack(:,:,k) / sum(sum(psfStack(:,:,k))); end逻辑说明:(D/r0)^(5/6)来自 Kolmogorov 相位结构函数的标度关系,控制相位起伏幅度。圆形孔径用掩膜实现,避免方形孔径带来的衍射伪影。每帧独立生成相位,得到统计独立但同分布的 PSF 序列,符合短曝光多帧的假设。
参数说明:r0越小湍流越强,PSF 越弥散;D/r0是决定退化程度的核心无量纲数,做实验时按 2、4、8 分档。相位屏分辨率 256 是精度和速度的折中,要更高精度可加到 512,但 FFT 耗时约增四倍。
3. 单帧复原:维纳滤波、Richardson-Lucy 与盲去卷积的 MATLAB 实现
3.1 维纳滤波:已知 PSF 时的最快基线
维纳滤波在频域做最小均方误差估计,公式为 F_hat = conj(H)/(|H|² + K) · G,K 是噪信比。它假设 PSF 已知,是评估其他算法时的基线。
% 频域维纳滤波,K 为噪信比 G = fft2(blurred_noisy); H = fft2(psf, size(blurred_noisy,1), size(blurred_noisy,2)); K = noise_var / var(img(:)); % 噪信比估计 F_hat = conj(H) ./ (abs(H).^2 + K) .* G; restored_wiener = real(ifft2(F_hat)); restored_wiener = min(max(restored_wiener, 0), 1); % 截断到有效范围逻辑说明:fft2(psf, M, N)把 PSF 补零到图像尺寸,保证频域相乘对应空域循环卷积。K 用噪声方差比图像方差估计,比手工试值更稳。最后截断到 [0,1] 抑制负值和过冲。
参数说明:K 偏大会过度平滑,偏小会放大噪声。实际中可对 K 做几档扫描,用 PSNR 选最优。维纳滤波在 PSF 准确、噪声不大时效果好,PSF 有误差时性能下降明显。
3.2 Richardson-Lucy 迭代复原与迭代次数控制
RL 算法基于泊松噪声假设,迭代式更新估计,适合光子计数受限的成像。
% Richardson-Lucy 去卷积 J = blurred_noisy; est = J; % 初始估计 numIter = 30; % 迭代次数 psfFlip = rot90(psf, 2); % PSF 翻转用于相关运算 for i = 1:numIter convEst = imfilter(est, psf, 'conv', 'symmetric'); ratio = J ./ max(convEst, 1e-8); % 防止除零 est = est .* imfilter(ratio, psfFlip, 'conv', 'symmetric'); est = min(max(est, 0), 1); end restored_rl = est;逻辑说明:max(convEst, 1e-8)避免暗区除零产生 NaN,这是 RL 最常见的崩溃点。PSF 翻转对应相关运算,是 RL 更新公式的要求。每轮截断保证非负且有界。
参数说明:迭代次数是核心。次数太少残留模糊,太多会放大噪声并产生斑点。经验上 20 到 50 之间,用 PSNR 曲线找拐点。噪声大时迭代次数要调低,或先做去噪再 RL。
3.3 盲去卷积:PSF 未知时的交替估计
实际场景里 PSF 往往未知。盲去卷积交替更新图像和 PSF,常用交替最小化或基于先验的迭代。
% 简化盲去卷积:交替更新 PSF 和图像 estImg = blurred_noisy; estPsf = ones(ksize) / ksize^2; % PSF 初始化为均匀核 outerIter = 15; for o = 1:outerIter % 固定 PSF,用 RL 更新图像 for i = 1:5 convEst = imfilter(estImg, estPsf, 'conv', 'symmetric'); ratio = blurred_noisy ./ max(convEst, 1e-8); estImg = estImg .* imfilter(ratio, rot90(estPsf,2), 'conv', 'symmetric'); estImg = min(max(estImg, 0), 1); end % 固定图像,更新 PSF convImg = imfilter(estImg, estPsf, 'conv', 'symmetric'); ratioP = blurred_noisy ./ max(convImg, 1e-8); estPsf = estPsf .* imfilter(rot90(estImg,2), ratioP, 'conv', 'symmetric'); estPsf = estPsf / sum(estPsf(:)); % 保持能量 end restored_blind = estImg;逻辑说明:外层循环交替,内层用少量 RL 迭代更新图像,再更新 PSF。PSF 每轮归一化防止能量漂移。这种朴素盲去卷积对初始化和噪声敏感,实际论文里会加稀疏先验或 TV 正则。
参数说明:外层 15 次、内层 5 次是速度与效果的折中。PSF 初始核尺寸要和真实核接近,否则容易收敛到错误解。噪声大时建议先维纳粗复原再盲去卷积。
4. 多帧配准融合复原:把湍流随机性变成优势
4.1 多帧退化数据的生成与配准
单帧信息不足,多帧的关键是各帧 PSF 不同,融合后能互补。先要配准,消除帧间整体位移。
% 生成多帧退化并做互相关配准 frames = zeros(size(img,1), size(img,2), N); for k = 1:N frames(:,:,k) = imfilter(img, psfStack(:,:,k), 'conv', 'symmetric'); frames(:,:,k) = imnoise(frames(:,:,k), 'gaussian', 0, noise_var); end % 以第一帧为参考,互相关求位移 ref = frames(:,:,1); aligned = frames; for k = 2:N c = normxcorr2(ref, frames(:,:,k)); [~, imax] = max(abs(c(:))); [ypeak, xpeak] = ind2sub(size(c), imax); shift = [ypeak - size(ref,1), xpeak - size(ref,2)]; aligned(:,:,k) = circshift(frames(:,:,k), -shift); end逻辑说明:normxcorr2对亮度变化鲁棒,比直接互相关稳。峰值位置减参考尺寸得到位移量,circshift反向补偿。循环位移会引入边缘环绕,必要时改用imtranslate加边界填充。
参数说明:配准精度直接影响融合效果,亚像素位移要用imregcorr或相位相关。帧数 N 越大融合增益越高,但配准误差累积也越明显,一般 10 到 30 帧。
4.2 多帧融合与去噪的联合处理
配准后做融合,简单平均能降噪但也会平均掉高频;加权融合按帧质量给权。
% 按帧清晰度加权融合 weights = zeros(1, N); for k = 1:N weights(k) = var(aligned(:,:,k), 0, 'all'); % 方差作为清晰度代理 end weights = weights / sum(weights); fused = zeros(size(ref)); for k = 1:N fused = fused + weights(k) * aligned(:,:,k); end % 融合后再做一次 RL 提升细节 restored_multi = deconvlucy(fused, psf, 15);逻辑说明:方差大的帧通常含更多高频,给更高权重。融合后噪声被平均降低,再 RL 能恢复更多细节。deconvlucy是 MATLAB 自带函数,内部已处理边界和迭代。
参数说明:权重用方差是简化做法,更严谨可用梯度能量或频域高频占比。RL 迭代次数在融合后可比单帧少,因为噪声已降。
| 方法 | 需要 PSF | 抗噪 | 细节恢复 | 计算量 |
|---|---|---|---|---|
| 维纳滤波 | 是 | 中 | 中 | 低 |
| Richardson-Lucy | 是 | 中低 | 高 | 中 |
| 盲去卷积 | 否 | 低 | 中高 | 高 |
| 多帧融合 | 部分 | 高 | 高 | 高 |
5. 复原质量评估与 MATLAB 调参排错技巧
5.1 用 PSNR、SSIM 和频域指标量化复原效果
主观看着好不够,论文和工程都要量化。PSNR 和 SSIM 是标配,湍流复原里还常看频域能量恢复。
% 计算 PSNR 和 SSIM psnrVal = psnr(restored_rl, img); ssimVal = ssim(restored_rl, img); % 频域高频能量占比,反映细节恢复 F_restored = fftshift(abs(fft2(restored_rl))); F_orig = fftshift(abs(fft2(img))); [M, Nn] = size(img); radius = M/8; [XX, YY] = meshgrid(-Nn/2:Nn/2-1, -M/2:M/2-1); maskHigh = (XX.^2 + YY.^2) > radius^2; highRatio = sum(F_restored(maskHigh)) / sum(F_orig(maskHigh));逻辑说明:PSNR 对亮度偏移敏感,SSIM 更贴近人眼。高频能量占比能看出算法是否真的恢复了细节,而不是只做了平滑。radius取 M/8 是经验分界,可按图像内容调整。
参数说明:PSNR 低于 25dB 通常说明复原不足或过冲;SSIM 低于 0.7 要检查配准和 PSF 估计。高频占比接近 1 说明细节恢复充分,超过 1 可能是噪声放大。
5.2 常见报错与调参排错清单
MATLAB 做图像复原时,几个坑反复出现。下面按现象给排查方向。
- 复原图出现网格状振铃:PSF 核尺寸太小或边界处理不当,加大 ksize,改用
symmetric或replicate。 - 结果全黑或全白:归一化丢失或除零,检查 PSF 求和是否为 1,
max(...,1e-8)是否加上。 - RL 迭代后噪声爆炸:迭代次数过多,降到 20 以内,或先做非局部均值去噪。
- 盲去卷积不收敛:PSF 初始核与真实差距大,先用维纳粗估核尺寸。
- PSNR 不升反降:可能过冲,检查是否截断到 [0,1],K 值是否偏小。
- 多帧融合后更糊:配准错误,用
normxcorr2可视化峰值确认位移方向。
提示:调参时固定随机种子,
rng(0),否则每次生成的噪声和相位不同,PSNR 波动会让你误判参数好坏。
5.3 从仿真到实测的迁移注意点
仿真里 PSF 已知,实测里未知,迁移时先做 PSF 标定,用点光源或已知靶标估计核。实测噪声往往不是高斯而是泊松加读出噪声,RL 比维纳更合适。另外实测有帧间非刚性形变,配准要上非刚性方法。把仿真调好的参数直接搬到实测,通常要重新扫一遍迭代次数和 K 值。
本文还有配套的精品资源,点击获取