Matlab实现Gerchberg-Saxton算法：从强度信息恢复相位与波前重建

简介：本资源是一套面向光学工程、计算成像与自适应光学方向初学者及科研人员的Matlab仿真教学工具，聚焦Gerchberg-Saxton（GS）迭代算法在光学相位恢复与波前重建中的原理实现与数值验证。资源包共30个文件，含10个核心Matlab脚本（如Gerchberg_Saxton_Algorithm.m、FT2Dc.m、a_simulate_DP.m等）、5幅算法中间过程可视化fig图、5张重建结果对比jpg图（含target.jpg、reconstractedimage.jpg、相位.jpg等），以及bin格式的预设振幅/相位数据、说明文档与README指南，整体压缩包大小为53.58MB。已有70人学习下载，适用于高校光学课程实验、毕业设计仿真验证及GS算法收敛性探究。用户可直接运行主流程脚本，调整参数观察频域/空域约束交替迭代过程，结合附赠的理论说明文档理解算法推导，并通过fig与jpg结果直观评估重建精度与波前形貌，具备完整可复现、可拓展、可教学的实践闭环。

1. 项目概述：从“丢失的相位”到“完整的波前”

在光学成像、全息术、自适应光学乃至天文观测等领域，我们常常面临一个核心挑战：探测器（如CCD相机）只能记录光波的强度信息，而至关重要的相位信息在传播过程中丢失了。这就好比我们拿到了一张照片，却不知道照片上每个像素点发出的光波是“同步”还是“异步”的，这直接决定了我们能否重建出物体的三维结构或精确的波前形状。Gerchberg-Saxton（GS）迭代算法，就是解决这一“相位恢复”问题的经典且强大的数学工具。它不需要复杂的干涉装置，仅通过采集到的强度信息，通过迭代计算，就能反推出丢失的相位，进而实现波前重建。

这个项目，就是利用Matlab构建一个完整的GS算法仿真系统。它不仅仅是一个算法实现，更是一个用于教学、研究和工程验证的综合性平台。你可以用它来理解GS算法的核心思想，模拟不同条件下的相位恢复过程（比如加入噪声、使用不同的初始猜测），并直观地观察迭代过程中波前是如何一步步被“拼凑”完整的。对于光学工程、物理、图像处理等领域的学生和研究者而言，亲手搭建并运行这样一个系统，远比阅读公式和论文来得深刻。

2. 核心原理拆解：GS算法是如何“猜”出相位的？

GS算法的精妙之处在于其简洁性和有效性。它的核心思想建立在两个基本约束条件上：在像平面（我们测量到强度的地方）和物平面（我们想恢复相位的地方）之间来回迭代，并强制施加各自的已知条件。

2.1 算法的数学与物理基础

整个过程可以概括为四个步骤的循环，我们以一个典型的傅里叶变换关系为例（即物平面和像平面是傅里叶变换对，这在许多衍射问题中成立）：

1. 初始化：在物平面，我们有一个复振幅场 ( U_{obj} = A_{obj} \cdot \exp(i \phi_{guess}) )。其中，( A_{obj} ) 是已知的或假设的物平面振幅（例如，一个均匀照明的孔径），而初始相位 ( \phi_{guess} ) 通常是一个随机猜测，或者设为零（平面波猜测）。

2. 正向传播（约束1）：将物平面的复振幅场通过一个数学变换（如傅里叶变换FFT）传播到像平面，得到 ( U_{img} = \mathcal{F}{U_{obj}} = A_{img} \cdot \exp(i \phi_{img}) )。

3. 像平面约束：这是我们掌握的唯一实测数据。用实际测量到的像平面强度 ( I_{measured} ) 的平方根，替换掉计算得到的像平面振幅 ( A_{img} )，同时保留计算得到的相位 ( \phi_{img} )。即，更新像平面复振幅为：( U'{img} = \sqrt{I{measured}} \cdot \exp(i \phi_{img}) )。这一步强制计算出的波前符合我们观测到的强度分布。

4. 反向传播（约束2）：将施加了约束后的像平面复振幅 ( U'{img} ) 通过逆变换（如逆傅里叶变换IFFT）传播回物平面，得到 ( U'{obj} = \mathcal{F}^{-1}{U'{img}} = A'{obj} \cdot \exp(i \phi'_{obj}) )。

5. 物平面约束：在物平面，我们通常知道光场的支持域（即光只存在于某个形状的孔径内）。因此，我们将 ( U'{obj} ) 的振幅 ( A'{obj} ) 替换为已知的物平面振幅 ( A_{obj} )（或在其支持域内保持不变，域外强制为零），同时保留新计算出的相位 ( \phi'{obj} )。即，更新物平面复振幅为：( U^{(new)}{obj} = A_{obj} \cdot \exp(i \phi'_{obj}) )。

至此，完成一次迭代。将更新后的 ( U^{(new)}_{obj} ) 作为下一次迭代的起点，重复步骤2-5。

注意：这里的“物平面”和“像平面”是广义的，取决于具体的光学系统模型。它们也可以是菲涅耳衍射下的两个平面，甚至是多次散射介质前后的两个平面。核心是找到一个可逆的数学变换来描述光场在两个平面间的传播。

2.2 为什么迭代会收敛？

算法的收敛性并非严格数学证明，但可以从物理上直观理解。每次迭代，我们都用“已知的事实”（像平面的实测强度、物平面的支持域或振幅）去修正计算出的波前。虽然初始猜测的相位可能是完全错误的，但经过多次这样的“事实修正”，计算出的波前会越来越同时满足两个平面的约束，其相位也就越来越接近真实情况。我们可以定义一个误差函数，如像平面计算强度与实测强度的均方根误差（RMSE），来监控收敛过程。通常，随着迭代进行，该误差会逐渐减小并趋于稳定。

3. 系统设计与Matlab实现框架

一个健壮的GS仿真系统不能只是一个算法循环，它需要包含数据生成、算法核心、可视化监控和性能评估等多个模块。

3.1 系统模块划分

我们的Matlab仿真系统主要包含以下模块：

1. 仿真数据生成模块：负责创建“地面真值”。即，先定义一个理想的物平面相位分布（如球面波、泽尼克像差、随机相位板等），结合已知的物平面振幅（如圆形孔径），通过正向传播（如FFT）计算出“理论上”的像平面复振幅，取其强度作为“实测数据” ( I_{measured} )。这样我们就拥有了一个完全可控的、已知答案的测试案例。

2. GS算法核心迭代模块：这是系统的心脏。它接收“实测强度”和“物平面振幅约束”作为输入，进行前述的迭代循环。该模块需要高度可配置，包括：

   • 最大迭代次数
   • 收敛阈值（误差低于某值则停止）
   • 初始相位猜测方式（随机相位、零相位等）
   • 传播模型选择（傅里叶变换、菲涅耳衍射等）

3. 可视化与监控模块：在迭代过程中实时或分阶段显示关键信息，对于理解算法至关重要。需要显示：

   • 每次迭代后恢复的物平面相位分布。
   • 每次迭代后计算出的像平面强度与“实测”强度的对比。
   • 误差曲线（如RMSE vs. 迭代次数）。

4. 性能分析与评估模块：在迭代结束后，定量评估恢复效果。常用指标包括：

   • 相位恢复误差：计算恢复相位与真实相位之间的差值（需处理相位包裹问题）。
   • 相关系数：计算恢复相位与真实相位的二维相关系数。
   • 收敛速度分析。

3.2 Matlab实现的关键技巧与代码结构

下面给出一个高度简化的核心框架，并附上关键技巧说明。

%% 主函数框架 function [recovered_phase, error_history] = GS_PhaseRetrieval(simulated_intensity, object_amplitude, max_iter, threshold) % simulated_intensity: 模拟的像平面强度（实测数据） % object_amplitude: 物平面的振幅约束（如孔径函数） % max_iter: 最大迭代次数 % threshold: 收敛阈值 [M, N] = size(simulated_intensity); target_amplitude = sqrt(simulated_intensity); % 像平面振幅约束 % 1. 初始化：随机相位猜测 initial_phase = 2 * pi * rand(M, N); U_object = object_amplitude .* exp(1i * initial_phase); error_history = zeros(max_iter, 1); for iter = 1:max_iter % 2. 正向传播 (FFT) U_image = fft2(U_object); current_phase_image = angle(U_image); % 3. 像平面约束：替换振幅，保留相位 U_image_constrained = target_amplitude .* exp(1i * current_phase_image); % 4. 反向传播 (IFFT) U_object_back = ifft2(U_image_constrained); current_phase_object = angle(U_object_back); % 5. 物平面约束：替换振幅，保留相位 U_object = object_amplitude .* exp(1i * current_phase_object); % 计算当前误差（像平面强度误差） computed_intensity = abs(fft2(U_object)).^2; error = sqrt(sum(sum((computed_intensity - simulated_intensity).^2))) / (M*N); error_history(iter) = error; % 可视化监控（每50次迭代显示一次） if mod(iter, 50) == 0 figure(1); subplot(2,2,1); imagesc(angle(U_object)); title(['Recovered Phase, Iter: ', num2str(iter)]); axis image; colorbar; subplot(2,2,2); plot(error_history(1:iter)); title('Error Convergence'); xlabel('Iteration'); ylabel('RMSE'); drawnow; end % 检查收敛 if error < threshold fprintf('Converged at iteration %d with error %.4e\n', iter, error); break; end end recovered_phase = angle(U_object); end

实操心得1：FFT的缩放与移位在光学仿真中，FFT/IFFT默认的输出是“数学中心”在矩阵的(1,1)角点。为了更直观地显示以零频为中心的光学衍射图样，我们经常需要使用fftshift和ifftshift函数。关键区别：fftshift用于将零频分量移动到数组中心，适用于显示；而在进行连续的FFT和IFFT运算时，为了保持数学正确性，通常需要在正向变换前使用ifftshift，在反向变换后使用fftshift。一个常见的传播模型是：U_f = fftshift(fft2(ifftshift(U_o)))。在我们的核心循环中，为了代码简洁和速度，有时可以省略shift，但必须清楚最终显示时需要做相应调整。

实操心得2：初始猜测的影响GS算法可能陷入局部极小值。使用随机相位作为初始猜测通常比零相位（平面波猜测）更容易收敛到全局解，尤其是在物平面相位变化剧烈的情况下。可以尝试多次使用不同的随机种子运行算法，选择恢复效果最好的一次。

4. 仿真案例深度实操：从简单到复杂

我们通过三个逐步深入的案例，来演示系统的强大功能和分析能力。

4.1 案例一：恢复简单倾斜相位

这是最基础的验证。假设物平面相位是一个简单的倾斜平面（相位随x方向线性增加）。

%% 案例1：生成倾斜相位 M = 512; N = 512; x = linspace(-1, 1, N); y = linspace(-1, 1, M); [X, Y] = meshgrid(x, y); % 定义物平面振幅（圆形孔径） aperture_radius = 0.4; object_amplitude = sqrt((X.^2 + Y.^2) < aperture_radius^2); % 定义真实的物平面相位（倾斜波前） true_phase = 10 * pi * X; % 相位随X线性变化 % 生成“地面真值”像平面强度 U_obj_true = object_amplitude .* exp(1i * true_phase); U_img_true = fftshift(fft2(ifftshift(U_obj_true))); simulated_intensity = abs(U_img_true).^2; % 加入少量噪声模拟真实情况 simulated_intensity = simulated_intensity + 0.01 * max(simulated_intensity(:)) * randn(M, N); simulated_intensity = max(simulated_intensity, 0); % 强度非负 % 调用GS算法 [recovered_phase, err_hist] = GS_PhaseRetrieval(simulated_intensity, object_amplitude, 200, 1e-6); % 评估：计算相位差（需要解包裹） phase_diff = recovered_phase - true_phase; phase_diff_unwrapped = unwrap(phase_diff, [], 2); % 沿x方向解包裹 figure; subplot(1,3,1); imagesc(true_phase); title('True Phase'); axis image; colorbar; subplot(1,3,2); imagesc(recovered_phase); title('Recovered Phase'); axis image; colorbar; subplot(1,3,3); imagesc(phase_diff_unwrapped); title('Unwrapped Phase Error'); axis image; colorbar; fprintf('Phase RMSE: %.4f radians\n', sqrt(mean(phase_diff_unwrapped(object_amplitude>0).^2)));

结果分析：对于这种简单的线性相位，GS算法通常能快速（几十次迭代内）且高精度地恢复。误差主要来源于数值计算精度和添加的噪声。通过观察误差曲线，可以看到典型的指数衰减式收敛。

4.2 案例二：恢复包含泽尼克像差的复杂波前

泽尼克多项式是描述光学像差的标准工具。我们模拟一个包含慧差（Z3）和球差（Z4）的复杂波前。

%% 案例2：泽尼克像差恢复 % 使用zernike函数生成泽尼克多项式（需自行实现或使用工具箱） % 假设我们有函数 Z = zernike(n, m, rho, theta) 返回归一化的泽尼克模式 [theta, rho] = cart2pol(X, Y); rho = rho / aperture_radius; % 归一化到孔径内 rho(rho > 1) = 0; % 定义真实相位 = a3*Z3 + a4*Z4 Z3 = zernike(1, 1, rho, theta); % 垂直慧差 (n=1, m=±1) Z4 = zernike(2, 0, rho, theta); % 初级球差 (n=2, m=0) true_phase_complex = 2*pi * (0.5*Z3 + 0.3*Z4); true_phase_complex(rho == 0) = 0; % 孔径外置零 % 生成仿真数据（同上） U_obj_true = object_amplitude .* exp(1i * true_phase_complex); % ... 后续传播、加噪声、调用GS算法与案例1类似 % 评估：将恢复的相位投影到泽尼克基上，分析像差系数 mask = object_amplitude > 0.5; phase_to_fit = recovered_phase(mask); % 使用最小二乘法拟合泽尼克系数 % A = [Z3(mask), Z4(mask), ...]; % coeffs = A \ phase_to_fit;

注意事项：复杂像差的收敛挑战对于复杂像差，尤其是高阶像差，GS算法可能对初始猜测更敏感。随机相位初始猜测可能无法收敛到正确解，或者收敛速度很慢。此时，可以采用“逐步增加复杂度”的策略：先用低分辨率或平滑的相位作为初始猜测运行少量迭代，再将结果作为高分辨率仿真的初始值。另一种方法是使用“输入-输出”算法等GS的变种，它们通过引入一个松弛参数来加速收敛或避免停滞。

4.3 案例三：在噪声与部分遮挡下的鲁棒性测试

真实实验数据充满噪声，且探测器可能存在坏点或部分遮挡。我们测试系统的鲁棒性。

%% 案例3：加入噪声和遮挡 % 使用案例1的简单倾斜相位 % 1. 加入更强的高斯噪声 noise_level = 0.1; % 噪声强度为峰值强度的10% simulated_intensity_noisy = simulated_intensity + noise_level * max(simulated_intensity(:)) * randn(M, N); simulated_intensity_noisy = max(simulated_intensity_noisy, 0); % 2. 加入中心遮挡（模拟探测器中心有灰尘或遮挡物） occlusion_radius = 20; % 遮挡像素半径 center = [M/2, N/2]; [Y_idx, X_idx] = meshgrid(1:N, 1:M); occlusion_mask = sqrt((X_idx - center(1)).^2 + (Y_idx - center(2)).^2) > occlusion_radius; simulated_intensity_occluded = simulated_intensity_noisy .* occlusion_mask; % 分别用干净数据、噪声数据、遮挡数据运行GS算法 % ...

结果对比与技巧：

• 噪声影响：中等程度的加性噪声通常不会阻止GS算法收敛，但会降低恢复相位的精度，并在最终结果中引入高频“毛刺”。可以在迭代循环中引入简单的低通滤波（例如，在每次迭代后对恢复的相位进行轻微的高斯滤波），但需谨慎，以免过滤掉真实的相位细节。
• 遮挡影响：像平面数据的部分缺失是更严峻的挑战。GS算法在像平面的约束是“用实测强度替换计算振幅”，在遮挡区域，实测强度为零，这相当于施加了一个非常强的错误约束，可能导致算法发散或恢复出完全错误的相位。解决策略：在像平面约束步骤，只对未被遮挡的区域进行振幅替换，遮挡区域保持计算出的振幅和相位不变。这需要修改核心算法中的约束步骤。

5. 高级话题与系统扩展

一个基础的GS系统搭建完成后，我们可以从多个方向对其进行扩展和深化研究。

5.1 GS算法的变种与性能对比

经典GS算法虽然稳定，但收敛速度有时较慢，且可能停滞。以下是两个重要的变种：

1. 输入-输出算法：由Fienup提出，是GS的高效改进版。它在物平面约束步骤引入了一个反馈机制：U_obj_new = U_obj + β * [Constraint(U_obj_back) - U_obj]其中，Constraint()代表物平面约束操作，β是一个介于0.5到1之间的反馈系数。这个简单的修改能显著加速收敛，尤其对于复杂物体。

2. 混合输入-输出算法：在输入-输出算法基础上，根据每次迭代的结果动态选择使用“输出”约束（GS型）还是“输入”约束（更激进），进一步提升了逃离局部极小值的能力。

在Matlab中实现这些变种并与经典GS对比，绘制收敛曲线，是深入理解迭代优化算法的绝佳练习。

5.2 结合实际应用场景：非干涉型表面形貌测量

GS算法的一个直接应用是非干涉光学轮廓测量。假设我们用一束准直光照射一个粗糙表面，散射光在远场形成散斑图样（像平面强度）。通过GS算法，可以从单幅散斑强度图中恢复出物体表面的高度起伏引起的相位变化，进而重建表面形貌。在仿真中，我们需要将物体表面高度h(x,y)映射为相位φ(x,y) = (4π/λ) * h(x,y)（对于正入射反射），然后按照前述流程仿真。这需要将菲涅耳衍射或角谱理论作为传播模型，而不仅仅是傅里叶变换。

5.3 系统性能的定量评估体系

建立一个完整的评估体系对于研究至关重要，应包括：

• 收敛性指标：像平面强度误差（RMSE）随迭代次数的变化。
• 准确性指标：
  • 恢复相位与真实相位的均方根误差（RMSE）。
  • 斯特列尔比：恢复波前与理想波前聚焦光斑峰值强度之比。
  • 泽尼克系数拟合残差。
• 鲁棒性指标：在不同信噪比（SNR）的噪声下，上述准确性指标的变化曲线。
• 效率指标：达到指定误差阈值所需的迭代次数和计算时间。

将这些指标封装成函数，可以对不同算法、不同参数进行自动化批量测试和比较。

6. 常见问题、调试技巧与避坑指南

在实际编写和运行GS仿真时，你会遇到各种问题。以下是一些典型问题及解决方案。

一个关键的调试技巧：分阶段验证不要一开始就用复杂数据和完整算法。构建一个“金字塔”式的调试流程：

1. 单元测试：先单独测试你的正向传播函数。给定一个已知相位，计算像平面强度，再手动计算（或用小规模FFT验证）是否正确。
2. 单步迭代：关闭循环，手动执行一次GS迭代。检查每个步骤（FFT、振幅替换、IFFT）的输入输出数据维度、类型和量级是否合理。
3. 简单案例：用倾斜相位这种有解析解或一目了然的案例进行测试。确保算法能正确恢复。
4. 逐步增加复杂度：在简单案例成功的基础上，再加入噪声、复杂像差、非对称孔径等元素。

最后，相位恢复是一个充满挑战和乐趣的领域。GS算法是打开这扇大门的钥匙。通过这个Matlab仿真系统，你获得的不只是一个工具，更是一套分析、解决波前重建问题的思维方法。我个人的体会是，多动手调整参数、多观察中间结果的变化、多思考每个约束背后的物理意义，比单纯追求算法的最终结果更重要。当你看到随机猜测的相位，经过几十次迭代后，神奇地收敛到预设的复杂波形时，那种对数学和物理之美的体会，是任何教科书都无法给予的。

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