简介本资源面向光学仿真与激光光束调控领域的科研人员及高校师生提供基于Gerchberg-SaxtonGS算法的高斯光束模式转换完整Matlab实现方案解决高阶拉盖尔-高斯光、一阶空心高斯光及贝塞尔-高斯光三类特殊光场的数值生成与相位反演问题。压缩包共10个文件含3个核心m函数GS.m、frt.m、ifrt.m、1个预存数据mat文件、5张运行结果图jpg及1张空心光束bmp图总大小4.07MB其中主函数main.m驱动全流程调用模块化子函数实现傅里叶变换、迭代优化与强度重构图像文件直观展示各目标光场的空间分布特征。已有138人学习下载配套代码经Matlab 2019b实测可直接运行无需额外配置涵盖光栅衍射、涡旋光束、夫琅禾费衍射等典型光学场景为光束整形、光学镊子及微纳操控等应用提供可复现的算法验证基础。1. 项目概述从高斯光到复杂光场的GS算法实现在光学设计、激光加工和光通信领域我们常常需要特定的光束形态。标准的高斯光束虽然应用广泛但在某些场景下其能量分布过于集中或者相位波前不够特殊无法满足需求。比如在光学镊子中我们需要“空心”的光束来无接触地捕获微粒在材料加工中可能需要一个中心能量为零的环形光斑来避免中心过热而在某些成像和传输应用中贝塞尔光束因其无衍射特性而备受青睐。然而激光器直接输出的往往是基础的高斯光束。如何将这种“标准件”转换成我们想要的“定制件”就是本次项目要解决的核心问题。这个项目就是利用经典的Gerchberg-SaxtonGS算法在Matlab平台上实现从基础高斯光束到高阶高斯光束、一阶空心高斯光束以及贝塞尔高斯光束的数值模拟与转换。简单来说它提供了一个数字化的“光学变换器”。你不需要昂贵的空间光调制器SLM或复杂的衍射光学元件DOE就能在计算机上验证你的光束设计是否可行这对于前期研究和教学演示至关重要。无论你是光学工程的学生还是从事激光应用的工程师通过复现这个项目你都能深入理解GS算法的迭代精髓并掌握在Matlab中构建和操控复杂光场的基本技能。2. 核心原理与GS算法深度解析2.1 认识我们的“原料”与“目标”在开始“烹饪”算法转换之前我们必须清楚“食材”输入光场和“菜品”目标光场的特性。基础高斯光束这是我们的起点也是最常见的激光光束模式TEM00模。其在垂直于传播方向的横截面上的光强分布遵循高斯函数中心最亮向外逐渐衰减。其数学表达式在束腰处相对简单通常作为我们算法的初始复振幅分布同时包含振幅和相位信息初始相位常设为零或平面波前。高阶高斯光束这里通常指拉盖尔-高斯LG光束或厄米-高斯HG光束。它们是波动方程在傍轴近似下的另一组完备解。LG光束在极坐标下描述具有螺旋状的相位波前携带轨道角动量和环形的强度分布。其模式由两个指数径向指数p和角向指数l决定。l≠0的LG光束中心强度为零呈现空心环状。HG光束则在直角坐标下描述光斑呈矩形对称分布。这些高阶模式在光通信、微粒操控等领域有特殊用途。一阶空心高斯光束这可以看作是LG模式的一个特例例如LG_0^1模式或者是一种近似描述。其核心特征是在光束传播轴中心位置光强始终为零形成一个“暗核”被一个亮环包围。这种光束对中心区域的粒子或样品几乎没有光损伤是光学捕获和生命科学研究的理想工具。贝塞尔高斯光束这是一种神奇的光束它虽然不是严格的贝塞尔光束需要无限能量但在有限孔径下其主瓣在一定的传播距离内几乎不发生衍射即具有“无衍射”特性。其横截面光强分布由贝塞尔函数调制的高斯包络构成呈现为中心一个亮斑和周围一系列同心圆环。这种特性在激光钻孔、光学成像和粒子导向等方面极具价值。2.2 GS算法在实空间和傅里叶空间之间“讨价还价”GS算法的核心思想非常直观它是一种迭代傅里叶变换算法用于在已知光场强度分布的两个平面上恢复出其相位信息或者反之。在我们的场景中这两个平面通常是源平面输入平面我们已知或设定振幅如高斯分布和初始猜测的相位。目标平面输出平面我们已知或想要的振幅分布如空心高斯分布但相位未知。算法通过反复在实空间施加振幅约束和傅里叶空间或另一个实空间施加传播约束之间切换让光场信息在两个约束条件下不断迭代最终收敛到一个同时近似满足两个约束条件的解即我们得到了目标平面的完整复振幅信息振幅和相位。其迭代步骤可以精炼为以下四步构成一个循环正向传播从当前迭代的源平面复振幅 ( U_{1}^{k}(x, y) )包含振幅和相位出发通过一个数学变换通常是傅里叶变换FFT模拟夫琅禾费衍射或角谱传播ASP模拟更一般的衍射传播到目标平面得到目标平面的复振幅 ( U_{2}^{k}(u, v) \mathcal{F}{U_{1}^{k}} )。施加目标振幅约束在目标平面我们强行用我们想要的目标振幅分布( A_{target}(u, v) ) 替换掉 ( U_{2}^{k}(u, v) ) 的振幅部分但同时保留其迭代计算得到的相位部分( \phi_{2}^{k}(u, v) )。形成新的复振幅( U_{2}^{k}(u, v) A_{target}(u, v) \cdot \exp(i \cdot \phi_{2}^{k}(u, v)) )。这是算法的关键一步将我们的设计目标“注入”系统。反向传播将施加了约束后的目标平面复振幅 ( U_{2}^{k}(u, v) ) 通过逆变换逆傅里叶变换IFFT或逆角谱传播传回源平面得到更新后的源平面复振幅 ( U_{1}^{k}(x, y) \mathcal{F}^{-1}{U_{2}^{k}} )。施加源振幅约束在源平面我们强行用已知的源振幅分布( A_{source}(x, y) )如高斯振幅替换掉 ( U_{1}^{k}(x, y) ) 的振幅部分同时保留其反向传播回来的相位部分( \phi_{1}^{k}(x, y) )。形成用于下一次迭代的复振幅( U_{1}^{k1}(x, y) A_{source}(x, y) \cdot \exp(i \cdot \phi_{1}^{k}(x, y)) )。循环判断将 ( U_{1}^{k1}(x, y) ) 作为新的起点重复步骤1-4。如此循环往复直到满足停止条件如迭代次数达到预设值或前后两次迭代结果的变化小于某个阈值。注意这个算法最终得到的是源平面所需的相位调制板即每次迭代后源平面的相位 ( \phi_{1}^{k} )。理论上如果我们用这个相位分布去调制初始的高斯光束那么在经过一段距离的衍射由算法中的传播模型定义后我们就能在目标平面得到想要的光强分布。在纯数值模拟中我们通常直接看最后一次正向传播后目标平面的光强是否与目标吻合。为什么GS算法有效你可以把它想象成两个人隔着屏风讨论一幅画的修改。一个人只知道画布的尺寸和想要的最终图案轮廓目标振幅另一个人只知道颜料盘的初始颜色分布源振幅。他们通过传递修改意见复振幅振幅是轮廓相位是笔触细节每次传递都只能按自己知道的规则修正一部分振幅约束保留对方传来的细节相位。经过多次来回他们最终能协同“画”出一幅在两个规则下都尽可能合理的画。相位在这里起到了承载和传递“结构信息”的关键作用算法通过迭代寻找一个最优的相位解来“调和”源和目标之间振幅的不匹配。3. Matlab实现的关键步骤与代码剖析下面我们将结合核心代码片段详细拆解如何在Matlab中实现这一过程。假设我们已经定义了必要的参数如网格大小N、采样间隔dx,dy、波长lambda、传播距离z等。3.1 光场初始化与网格生成任何光学模拟的第一步都是建立计算网格。我们需要在源平面和目标平面创建二维坐标网格。% 定义参数 N 512; % 网格点数通常为2的幂次便于FFT L 0.01; % 物理尺寸10mm x 10mm dx L/N; dy dx; % 采样间隔 x linspace(-L/2, L/2-dx, N); % 生成对称坐标轴 y x; [X, Y] meshgrid(x, y); % 生成二维网格 % 定义波长和传播距离 lambda 632.8e-9; % 氦氖激光波长单位米 z 0.5; % 传播距离单位米 k 2*pi/lambda; % 波数关键点linspace生成的点包含了左端点-L/2但不包含右端点L/2这是为了满足离散傅里叶变换的周期性边界条件避免频谱泄露。使用meshgrid生成网格矩阵X和Y这是后续所有二维函数计算的基础。3.2 构建源光束与目标光束的振幅分布接下来我们需要创建初始高斯光束源振幅和我们想要的目标光束振幅。% 1. 源振幅基础高斯光束 w0 L/10; % 高斯光束束腰半径可根据需要调整 A_source exp(-(X.^2 Y.^2) / w0^2); % 高斯振幅分布 phase_initial zeros(N); % 初始相位设为平面波或随机相位 U_source A_source .* exp(1i * phase_initial); % 初始复振幅 % 2. 目标振幅示例一阶空心高斯光束 (近似为LG01模式) % 首先计算径向坐标r和角向坐标theta r sqrt(X.^2 Y.^2); theta atan2(Y, X); % 拉盖尔-高斯LG_0^1模式的振幅分布 (p0, l1) % 注意这里忽略了归一化常数专注于分布形状 p 0; l 1; % 拉盖尔多项式部分 (对于p0 L_0^|l|(x)1) % 整体振幅分布 r^|l| * exp(-r^2/w0^2) * 角向相位因子exp(i*l*theta)的振幅部分为1 A_target_LG01 (sqrt(2)*r/w0).^abs(l) .* exp(-r.^2 / w0^2); % 这是振幅不含相位 % 对于空心光束我们通常只关心强度图案所以目标振幅就是sqrt(强度) A_target A_target_LG01 / max(max(A_target_LG01)); % 归一化到[0,1]关键点源光束A_source是一个实值矩阵代表光场的振幅。初始相位phase_initial设为零代表平面波前。这是一个常见的起点但有时为了打破对称性、加速收敛也会使用随机相位初始化。目标光束A_target的构建需要根据具体模式选择公式。对于高阶高斯光束LG/HG需调用对应的多项式函数。对于贝塞尔高斯光束则需要结合贝塞尔函数和高斯包络。代码中展示了LG01模式的构建其强度分布中心为零呈现环形。归一化操作(/ max(max(...))) 非常重要。它确保目标振幅的最大值为1避免在迭代过程中出现数值不稳定或能量不匹配的问题。实际物理系统中总能量是守恒的但在GS算法的振幅替换步骤中我们更关注分布形状。3.3 实现GS算法迭代循环这是项目的核心引擎。我们将实现一个完整的GS迭代循环并监控其收敛情况。% GS算法参数 max_iter 200; % 最大迭代次数 error zeros(1, max_iter); % 记录每次迭代的误差 % 初始化迭代变量 U1 U_source; % 从源场开始 for iter 1:max_iter % --- 步骤1: 正向传播 (使用角谱理论适用于任何距离) --- % 计算角谱传递函数 fx (-N/2:N/2-1)/(N*dx); % 频率坐标 [FX, FY] meshgrid(fx, fx); H exp(1i * k * z * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 传递函数 H(isnan(H)) 0; % 处理可能出现的NaN当根号内为负时表示倏逝波可忽略 H fftshift(H); % 将零频移到中心与fft2结果匹配 U1_fft fft2(U1); % 源场频谱 U2 ifft2(U1_fft .* H); % 传播到目标场 % --- 步骤2: 在目标面施加振幅约束 --- amp_U2 abs(U2); phase_U2 angle(U2); % 保留相位替换为目标振幅 U2_prime A_target .* exp(1i * phase_U2); % --- 步骤3: 反向传播回源面 --- U2_prime_fft fft2(U2_prime); U1_prime ifft2(U2_prime_fft .* conj(H)); % 注意反向传播使用传递函数的共轭 % --- 步骤4: 在源面施加振幅约束 --- amp_U1_prime abs(U1_prime); phase_U1_prime angle(U1_prime); % 保留相位替换为源振幅 U1_next A_source .* exp(1i * phase_U1_prime); % --- 计算误差可选用于监控收敛 --- % 常用误差定义为目标面计算振幅与目标振幅的均方根误差 error(iter) sqrt(sum(sum((amp_U2 - A_target).^2)) / (N*N)); % --- 为下一次迭代更新源场 --- U1 U1_next; % --- 可选每50次迭代显示一次进度 --- if mod(iter, 50) 0 fprintf(迭代次数: %d, 当前误差: %.6f\n, iter, error(iter)); end end关键点解析与实操心得传播方法的选择代码中使用了角谱传播Angular Spectrum Method, ASM通过H传递函数实现。与简单的夫琅禾费近似仅一个FFT相比ASM在数学上是严格的适用于近场和远场衍射计算只要采样满足奈奎斯特准则。fftshift的操作是为了让频率坐标与传递函数对齐这是使用ASM时极易出错的地方。反向传播的传递函数反向传播时我们使用了conj(H)即传递函数的共轭。这是因为从目标面回到源面在数学上相当于正向传播的逆过程其传递函数是正向的复共轭。这是实现正确反向传播的关键。误差度量误差函数error用于监控算法收敛。这里使用了目标面计算振幅amp_U2与理想目标振幅A_target之间的均方根误差RMSE。你可以看到误差随着迭代次数增加而下降的趋势。如果误差曲线不再明显下降说明算法已收敛或陷入局部解。陷入停滞与解决方案GS算法可能陷入局部最小值导致误差无法进一步降低。一个经典的改进是输入输出算法Input-Output Algorithm其更新规则为U1_next U1 beta * (A_source.*exp(1i*phase_U1_prime) - U1)其中beta是一个松弛因子通常在0.5到1之间。这相当于在施加源约束时不仅用新相位还混合了上一次迭代的结果能有效跳出局部解加速收敛。在原代码基础上实现这个改进并不复杂但能显著提升性能。3.4 结果可视化与相位提取迭代完成后我们需要分析结果最重要的就是查看最终目标面的光强分布以及计算得到的源面相位板全息图。% 最终结果用最后一次迭代的U1进行最后一次正向传播得到最终的目标场 U1_final_fft fft2(U1); U2_final ifft2(U1_final_fft .* H); I_target_final abs(U2_final).^2; % 最终目标面光强 phase_source_final angle(U1); % 最终源面所需的相位调制分布 % 可视化 figure(Position, [100, 100, 1200, 400]); % 子图1初始高斯光强 subplot(1,3,1); imagesc(x*1e3, y*1e3, abs(U_source).^2); % 单位转换为mm axis image; colorbar; colormap(hot); xlabel(x (mm)); ylabel(y (mm)); title(初始高斯光束强度); % 子图2最终得到的目标光强 subplot(1,3,2); imagesc(x*1e3, y*1e3, I_target_final); axis image; colorbar; colormap(hot); xlabel(x (mm)); ylabel(y (mm)); title(GS算法生成的目标光束强度 (空心高斯)); % 子图3计算得到的源面相位板全息图 subplot(1,3,3); imagesc(x*1e3, y*1e3, mod(phase_source_final, 2*pi)); % 相位包裹在[0, 2π] axis image; colorbar; colormap(jet); xlabel(x (mm)); ylabel(y (mm)); title(源面所需相位分布 (包裹后));关键点相位包裹计算得到的相位phase_source_final取值范围是(-π, π]或更广。为了显示或用于生成全息图通常需要将其“包裹”到[0, 2π)区间使用mod(phase, 2*pi)操作。这个包裹后的相位图理论上可以加载到空间光调制器SLM上用于实际调制入射高斯光。结果评估对比子图1和子图2可以直观判断GS算法的效果。理想情况下子图2应该非常接近我们预设的A_target.^2分布。由于算法是数值逼近且受限于采样、迭代次数等结果边缘可能会有一些噪声或畸变。颜色映射强度图常用hot或gray色图相位图则用jet或hsv这类循环色图能清晰显示相位的周期性变化。4. 扩展到高阶高斯与贝塞尔高斯光束上述流程以空心高斯光束为例。要生成其他光束只需修改A_target的定义。4.1 生成高阶拉盖尔-高斯LG光束% 参数定义 p 1; % 径向指数 l 2; % 角向指数 w0_target w0; % 目标光束束腰可与源不同 % 计算LG模式的振幅分布未归一化 % 首先计算归一化径向坐标 rho sqrt(2)*r/w0_target rho sqrt(2) * r / w0_target; % 广义拉盖尔多项式 L_p^|l| (需要自定义或调用符号工具箱函数) % 这里以p1为例 L_1^|l|(x) (|l|1 - x) L (abs(l)1 - rho.^2); % 对于p1 % LG模式振幅 A_target_LG (rho.^abs(l)) .* L .* exp(-rho.^2/2); % 角向相位因子 exp(i*l*theta) 在GS算法中会被迭代出的相位替代这里只需振幅 A_target A_target_LG / max(max(A_target_LG)); % 归一化注意对于更通用的p需要实现或调用广义拉盖尔多项式的计算函数。Matlab符号数学工具箱有laguerreL函数但用于数值矩阵时需注意向量化操作。4.2 生成贝塞尔高斯光束% 参数定义 k_t 2*pi * 1e4; % 横向波矢决定贝塞尔光束的中心亮斑尺寸需根据需求调整 % 零阶贝塞尔高斯光束振幅分布 (J0) A_target_BG besselj(0, k_t * r) .* exp(-r.^2 / w0_target^2); A_target A_target_BG / max(max(A_target_BG)); % 归一化关键点k_t是一个关键参数k_t * r决定了贝塞尔函数的第一个零点位置从而控制了中心亮斑和环的间距。k_t越大中心斑越小环越密。需要根据模拟的物理尺寸r的范围合理选择避免采样不足导致失真。5. 常见问题、调试技巧与性能优化在实际编写和运行代码时你肯定会遇到各种问题。以下是我从多次实践中总结的“避坑指南”。5.1 算法不收敛或结果杂乱症状迭代误差不下降或最终目标光强图充满噪声与预期形状相差甚远。排查与解决检查传播模型确认正向和反向传播的传递函数H是否正确特别是fftshift和ifftshift的使用顺序。一个快速验证的方法是对一个点源delta函数进行传播看结果是否符合物理直觉如产生球面波。检查振幅约束确保A_source和A_target都是归一化后的非负实数矩阵。目标振幅不能有负值或复数。可视化它们确认分布符合预期。初始相位尝试如果使用全零平面波初始相位不收敛可以尝试使用随机相位初始化phase_initial 2*pi*rand(N);。随机相位能提供更多的自由度帮助算法跳出平庸解尤其对于中心对称的目标如空心光束有时效果更好。引入松弛因子如前所述采用输入输出算法。将源平面的更新规则改为U1_next U1 beta * (A_source.*exp(1i*phase_U1_prime) - U1);beta从0.5开始尝试。增加迭代次数有些复杂模式需要更多迭代才能收敛。将max_iter增加到500或1000。5.2 结果出现“孪生像”或共轭像症状在目标面除了想要的光斑在其对称位置出现一个类似的、较弱的“鬼影”。原因这是GS算法和傅里叶变换性质导致的常见问题。当源面和目标面的振幅约束都是实对称或厄米对称时算法求得的相位解可能不是唯一的会导致共轭像的出现。解决破坏对称性在源振幅或目标振幅中引入轻微的非对称性。例如将高斯光源的束腰中心稍微偏离计算网格中心。使用非对称初始相位使用随机初始相位而非零相位。采用加权GS算法在迭代初期对目标振幅约束进行松弛逐渐加强。5.3 采样与网格设置导致的失真症状生成的光束边缘出现锯齿、断裂或者贝塞尔光束的环不圆、不连续。排查奈奎斯特采样定理确保你的采样间隔dx足够小能够分辨光场中最高的空间频率。对于贝塞尔光束其空间频率与k_t相关。一个经验法则是dx lambda/(2*NA)其中NA是数值孔径的估计。如果出现混叠尝试增大N网格点数或减小L物理尺寸以降低dx。网格大小L确保计算窗口L足够大能够完整包含光束的主要能量部分避免能量在边界处被截断导致衍射效应干扰结果。可以观察初始高斯光束在网格边缘是否已衰减到接近零。传播距离z如果使用角谱传播距离z可以是任意的。但如果使用夫琅禾费近似U2 fft2(U1)则必须满足夫琅禾费条件z (N*dx)^2 / lambda。不满足时结果会错误。5.4 Matlab代码性能优化向量化与预计算所有操作都应使用矩阵运算避免循环。像H这样的传递函数应在循环外计算一次。内存管理对于大网格如N1024或更大复数矩阵会占用较多内存。使用single单精度而非默认的double双精度可以减半内存使用并在大多数情况下保持足够精度U1 single(U_source);。并行计算GS算法本身迭代间有依赖难以并行。但如果你需要为大量不同的目标光束计算相位图可以考虑使用parfor并行循环来处理不同的目标。收敛判断不必总是跑满max_iter次。可以在循环内加入判断当误差error(iter)的变化小于某个容差如1e-6时用break语句提前退出循环节省计算时间。我个人在实现这个项目时最大的体会是GS算法是一个“艺术多于科学”的过程。理论上很简单但调参松弛因子、迭代次数、初始条件对结果质量影响巨大。没有一套放之四海而皆准的参数。最好的方法是先从一个简单目标如将高斯光转换成另一个不同束腰的高斯光开始调试确保整个管道正确无误然后再挑战复杂的空心光或贝塞尔光。可视化每一个中间步骤如每次迭代后的目标面光强对于理解算法的行为非常有帮助。最后记住这个相位解是数值优化的结果它可能不是物理最优解但对于验证概念和指导初步实验已经足够强大。本文还有配套的精品资源点击获取 
长河君 · 2026/9/4 19:39:27 · 阅读 1.2k 基于GS算法与Matlab实现高斯光束到高阶/贝塞尔光束的数值模拟
简介本资源面向光学仿真与激光光束调控领域的科研人员及高校师生提供基于Gerchberg-SaxtonGS算法的高斯光束模式转换完整Matlab实现方案解决高阶拉盖尔-高斯光、一阶空心高斯光及贝塞尔-高斯光三类特殊光场的数值生成与相位反演问题。压缩包共10个文件含3个核心m函数GS.m、frt.m、ifrt.m、1个预存数据mat文件、5张运行结果图jpg及1张空心光束bmp图总大小4.07MB其中主函数main.m驱动全流程调用模块化子函数实现傅里叶变换、迭代优化与强度重构图像文件直观展示各目标光场的空间分布特征。已有138人学习下载配套代码经Matlab 2019b实测可直接运行无需额外配置涵盖光栅衍射、涡旋光束、夫琅禾费衍射等典型光学场景为光束整形、光学镊子及微纳操控等应用提供可复现的算法验证基础。1. 项目概述从高斯光到复杂光场的GS算法实现在光学设计、激光加工和光通信领域我们常常需要特定的光束形态。标准的高斯光束虽然应用广泛但在某些场景下其能量分布过于集中或者相位波前不够特殊无法满足需求。比如在光学镊子中我们需要“空心”的光束来无接触地捕获微粒在材料加工中可能需要一个中心能量为零的环形光斑来避免中心过热而在某些成像和传输应用中贝塞尔光束因其无衍射特性而备受青睐。然而激光器直接输出的往往是基础的高斯光束。如何将这种“标准件”转换成我们想要的“定制件”就是本次项目要解决的核心问题。这个项目就是利用经典的Gerchberg-SaxtonGS算法在Matlab平台上实现从基础高斯光束到高阶高斯光束、一阶空心高斯光束以及贝塞尔高斯光束的数值模拟与转换。简单来说它提供了一个数字化的“光学变换器”。你不需要昂贵的空间光调制器SLM或复杂的衍射光学元件DOE就能在计算机上验证你的光束设计是否可行这对于前期研究和教学演示至关重要。无论你是光学工程的学生还是从事激光应用的工程师通过复现这个项目你都能深入理解GS算法的迭代精髓并掌握在Matlab中构建和操控复杂光场的基本技能。2. 核心原理与GS算法深度解析2.1 认识我们的“原料”与“目标”在开始“烹饪”算法转换之前我们必须清楚“食材”输入光场和“菜品”目标光场的特性。基础高斯光束这是我们的起点也是最常见的激光光束模式TEM00模。其在垂直于传播方向的横截面上的光强分布遵循高斯函数中心最亮向外逐渐衰减。其数学表达式在束腰处相对简单通常作为我们算法的初始复振幅分布同时包含振幅和相位信息初始相位常设为零或平面波前。高阶高斯光束这里通常指拉盖尔-高斯LG光束或厄米-高斯HG光束。它们是波动方程在傍轴近似下的另一组完备解。LG光束在极坐标下描述具有螺旋状的相位波前携带轨道角动量和环形的强度分布。其模式由两个指数径向指数p和角向指数l决定。l≠0的LG光束中心强度为零呈现空心环状。HG光束则在直角坐标下描述光斑呈矩形对称分布。这些高阶模式在光通信、微粒操控等领域有特殊用途。一阶空心高斯光束这可以看作是LG模式的一个特例例如LG_0^1模式或者是一种近似描述。其核心特征是在光束传播轴中心位置光强始终为零形成一个“暗核”被一个亮环包围。这种光束对中心区域的粒子或样品几乎没有光损伤是光学捕获和生命科学研究的理想工具。贝塞尔高斯光束这是一种神奇的光束它虽然不是严格的贝塞尔光束需要无限能量但在有限孔径下其主瓣在一定的传播距离内几乎不发生衍射即具有“无衍射”特性。其横截面光强分布由贝塞尔函数调制的高斯包络构成呈现为中心一个亮斑和周围一系列同心圆环。这种特性在激光钻孔、光学成像和粒子导向等方面极具价值。2.2 GS算法在实空间和傅里叶空间之间“讨价还价”GS算法的核心思想非常直观它是一种迭代傅里叶变换算法用于在已知光场强度分布的两个平面上恢复出其相位信息或者反之。在我们的场景中这两个平面通常是源平面输入平面我们已知或设定振幅如高斯分布和初始猜测的相位。目标平面输出平面我们已知或想要的振幅分布如空心高斯分布但相位未知。算法通过反复在实空间施加振幅约束和傅里叶空间或另一个实空间施加传播约束之间切换让光场信息在两个约束条件下不断迭代最终收敛到一个同时近似满足两个约束条件的解即我们得到了目标平面的完整复振幅信息振幅和相位。其迭代步骤可以精炼为以下四步构成一个循环正向传播从当前迭代的源平面复振幅 ( U_{1}^{k}(x, y) )包含振幅和相位出发通过一个数学变换通常是傅里叶变换FFT模拟夫琅禾费衍射或角谱传播ASP模拟更一般的衍射传播到目标平面得到目标平面的复振幅 ( U_{2}^{k}(u, v) \mathcal{F}{U_{1}^{k}} )。施加目标振幅约束在目标平面我们强行用我们想要的目标振幅分布( A_{target}(u, v) ) 替换掉 ( U_{2}^{k}(u, v) ) 的振幅部分但同时保留其迭代计算得到的相位部分( \phi_{2}^{k}(u, v) )。形成新的复振幅( U_{2}^{k}(u, v) A_{target}(u, v) \cdot \exp(i \cdot \phi_{2}^{k}(u, v)) )。这是算法的关键一步将我们的设计目标“注入”系统。反向传播将施加了约束后的目标平面复振幅 ( U_{2}^{k}(u, v) ) 通过逆变换逆傅里叶变换IFFT或逆角谱传播传回源平面得到更新后的源平面复振幅 ( U_{1}^{k}(x, y) \mathcal{F}^{-1}{U_{2}^{k}} )。施加源振幅约束在源平面我们强行用已知的源振幅分布( A_{source}(x, y) )如高斯振幅替换掉 ( U_{1}^{k}(x, y) ) 的振幅部分同时保留其反向传播回来的相位部分( \phi_{1}^{k}(x, y) )。形成用于下一次迭代的复振幅( U_{1}^{k1}(x, y) A_{source}(x, y) \cdot \exp(i \cdot \phi_{1}^{k}(x, y)) )。循环判断将 ( U_{1}^{k1}(x, y) ) 作为新的起点重复步骤1-4。如此循环往复直到满足停止条件如迭代次数达到预设值或前后两次迭代结果的变化小于某个阈值。注意这个算法最终得到的是源平面所需的相位调制板即每次迭代后源平面的相位 ( \phi_{1}^{k} )。理论上如果我们用这个相位分布去调制初始的高斯光束那么在经过一段距离的衍射由算法中的传播模型定义后我们就能在目标平面得到想要的光强分布。在纯数值模拟中我们通常直接看最后一次正向传播后目标平面的光强是否与目标吻合。为什么GS算法有效你可以把它想象成两个人隔着屏风讨论一幅画的修改。一个人只知道画布的尺寸和想要的最终图案轮廓目标振幅另一个人只知道颜料盘的初始颜色分布源振幅。他们通过传递修改意见复振幅振幅是轮廓相位是笔触细节每次传递都只能按自己知道的规则修正一部分振幅约束保留对方传来的细节相位。经过多次来回他们最终能协同“画”出一幅在两个规则下都尽可能合理的画。相位在这里起到了承载和传递“结构信息”的关键作用算法通过迭代寻找一个最优的相位解来“调和”源和目标之间振幅的不匹配。3. Matlab实现的关键步骤与代码剖析下面我们将结合核心代码片段详细拆解如何在Matlab中实现这一过程。假设我们已经定义了必要的参数如网格大小N、采样间隔dx,dy、波长lambda、传播距离z等。3.1 光场初始化与网格生成任何光学模拟的第一步都是建立计算网格。我们需要在源平面和目标平面创建二维坐标网格。% 定义参数 N 512; % 网格点数通常为2的幂次便于FFT L 0.01; % 物理尺寸10mm x 10mm dx L/N; dy dx; % 采样间隔 x linspace(-L/2, L/2-dx, N); % 生成对称坐标轴 y x; [X, Y] meshgrid(x, y); % 生成二维网格 % 定义波长和传播距离 lambda 632.8e-9; % 氦氖激光波长单位米 z 0.5; % 传播距离单位米 k 2*pi/lambda; % 波数关键点linspace生成的点包含了左端点-L/2但不包含右端点L/2这是为了满足离散傅里叶变换的周期性边界条件避免频谱泄露。使用meshgrid生成网格矩阵X和Y这是后续所有二维函数计算的基础。3.2 构建源光束与目标光束的振幅分布接下来我们需要创建初始高斯光束源振幅和我们想要的目标光束振幅。% 1. 源振幅基础高斯光束 w0 L/10; % 高斯光束束腰半径可根据需要调整 A_source exp(-(X.^2 Y.^2) / w0^2); % 高斯振幅分布 phase_initial zeros(N); % 初始相位设为平面波或随机相位 U_source A_source .* exp(1i * phase_initial); % 初始复振幅 % 2. 目标振幅示例一阶空心高斯光束 (近似为LG01模式) % 首先计算径向坐标r和角向坐标theta r sqrt(X.^2 Y.^2); theta atan2(Y, X); % 拉盖尔-高斯LG_0^1模式的振幅分布 (p0, l1) % 注意这里忽略了归一化常数专注于分布形状 p 0; l 1; % 拉盖尔多项式部分 (对于p0 L_0^|l|(x)1) % 整体振幅分布 r^|l| * exp(-r^2/w0^2) * 角向相位因子exp(i*l*theta)的振幅部分为1 A_target_LG01 (sqrt(2)*r/w0).^abs(l) .* exp(-r.^2 / w0^2); % 这是振幅不含相位 % 对于空心光束我们通常只关心强度图案所以目标振幅就是sqrt(强度) A_target A_target_LG01 / max(max(A_target_LG01)); % 归一化到[0,1]关键点源光束A_source是一个实值矩阵代表光场的振幅。初始相位phase_initial设为零代表平面波前。这是一个常见的起点但有时为了打破对称性、加速收敛也会使用随机相位初始化。目标光束A_target的构建需要根据具体模式选择公式。对于高阶高斯光束LG/HG需调用对应的多项式函数。对于贝塞尔高斯光束则需要结合贝塞尔函数和高斯包络。代码中展示了LG01模式的构建其强度分布中心为零呈现环形。归一化操作(/ max(max(...))) 非常重要。它确保目标振幅的最大值为1避免在迭代过程中出现数值不稳定或能量不匹配的问题。实际物理系统中总能量是守恒的但在GS算法的振幅替换步骤中我们更关注分布形状。3.3 实现GS算法迭代循环这是项目的核心引擎。我们将实现一个完整的GS迭代循环并监控其收敛情况。% GS算法参数 max_iter 200; % 最大迭代次数 error zeros(1, max_iter); % 记录每次迭代的误差 % 初始化迭代变量 U1 U_source; % 从源场开始 for iter 1:max_iter % --- 步骤1: 正向传播 (使用角谱理论适用于任何距离) --- % 计算角谱传递函数 fx (-N/2:N/2-1)/(N*dx); % 频率坐标 [FX, FY] meshgrid(fx, fx); H exp(1i * k * z * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 传递函数 H(isnan(H)) 0; % 处理可能出现的NaN当根号内为负时表示倏逝波可忽略 H fftshift(H); % 将零频移到中心与fft2结果匹配 U1_fft fft2(U1); % 源场频谱 U2 ifft2(U1_fft .* H); % 传播到目标场 % --- 步骤2: 在目标面施加振幅约束 --- amp_U2 abs(U2); phase_U2 angle(U2); % 保留相位替换为目标振幅 U2_prime A_target .* exp(1i * phase_U2); % --- 步骤3: 反向传播回源面 --- U2_prime_fft fft2(U2_prime); U1_prime ifft2(U2_prime_fft .* conj(H)); % 注意反向传播使用传递函数的共轭 % --- 步骤4: 在源面施加振幅约束 --- amp_U1_prime abs(U1_prime); phase_U1_prime angle(U1_prime); % 保留相位替换为源振幅 U1_next A_source .* exp(1i * phase_U1_prime); % --- 计算误差可选用于监控收敛 --- % 常用误差定义为目标面计算振幅与目标振幅的均方根误差 error(iter) sqrt(sum(sum((amp_U2 - A_target).^2)) / (N*N)); % --- 为下一次迭代更新源场 --- U1 U1_next; % --- 可选每50次迭代显示一次进度 --- if mod(iter, 50) 0 fprintf(迭代次数: %d, 当前误差: %.6f\n, iter, error(iter)); end end关键点解析与实操心得传播方法的选择代码中使用了角谱传播Angular Spectrum Method, ASM通过H传递函数实现。与简单的夫琅禾费近似仅一个FFT相比ASM在数学上是严格的适用于近场和远场衍射计算只要采样满足奈奎斯特准则。fftshift的操作是为了让频率坐标与传递函数对齐这是使用ASM时极易出错的地方。反向传播的传递函数反向传播时我们使用了conj(H)即传递函数的共轭。这是因为从目标面回到源面在数学上相当于正向传播的逆过程其传递函数是正向的复共轭。这是实现正确反向传播的关键。误差度量误差函数error用于监控算法收敛。这里使用了目标面计算振幅amp_U2与理想目标振幅A_target之间的均方根误差RMSE。你可以看到误差随着迭代次数增加而下降的趋势。如果误差曲线不再明显下降说明算法已收敛或陷入局部解。陷入停滞与解决方案GS算法可能陷入局部最小值导致误差无法进一步降低。一个经典的改进是输入输出算法Input-Output Algorithm其更新规则为U1_next U1 beta * (A_source.*exp(1i*phase_U1_prime) - U1)其中beta是一个松弛因子通常在0.5到1之间。这相当于在施加源约束时不仅用新相位还混合了上一次迭代的结果能有效跳出局部解加速收敛。在原代码基础上实现这个改进并不复杂但能显著提升性能。3.4 结果可视化与相位提取迭代完成后我们需要分析结果最重要的就是查看最终目标面的光强分布以及计算得到的源面相位板全息图。% 最终结果用最后一次迭代的U1进行最后一次正向传播得到最终的目标场 U1_final_fft fft2(U1); U2_final ifft2(U1_final_fft .* H); I_target_final abs(U2_final).^2; % 最终目标面光强 phase_source_final angle(U1); % 最终源面所需的相位调制分布 % 可视化 figure(Position, [100, 100, 1200, 400]); % 子图1初始高斯光强 subplot(1,3,1); imagesc(x*1e3, y*1e3, abs(U_source).^2); % 单位转换为mm axis image; colorbar; colormap(hot); xlabel(x (mm)); ylabel(y (mm)); title(初始高斯光束强度); % 子图2最终得到的目标光强 subplot(1,3,2); imagesc(x*1e3, y*1e3, I_target_final); axis image; colorbar; colormap(hot); xlabel(x (mm)); ylabel(y (mm)); title(GS算法生成的目标光束强度 (空心高斯)); % 子图3计算得到的源面相位板全息图 subplot(1,3,3); imagesc(x*1e3, y*1e3, mod(phase_source_final, 2*pi)); % 相位包裹在[0, 2π] axis image; colorbar; colormap(jet); xlabel(x (mm)); ylabel(y (mm)); title(源面所需相位分布 (包裹后));关键点相位包裹计算得到的相位phase_source_final取值范围是(-π, π]或更广。为了显示或用于生成全息图通常需要将其“包裹”到[0, 2π)区间使用mod(phase, 2*pi)操作。这个包裹后的相位图理论上可以加载到空间光调制器SLM上用于实际调制入射高斯光。结果评估对比子图1和子图2可以直观判断GS算法的效果。理想情况下子图2应该非常接近我们预设的A_target.^2分布。由于算法是数值逼近且受限于采样、迭代次数等结果边缘可能会有一些噪声或畸变。颜色映射强度图常用hot或gray色图相位图则用jet或hsv这类循环色图能清晰显示相位的周期性变化。4. 扩展到高阶高斯与贝塞尔高斯光束上述流程以空心高斯光束为例。要生成其他光束只需修改A_target的定义。4.1 生成高阶拉盖尔-高斯LG光束% 参数定义 p 1; % 径向指数 l 2; % 角向指数 w0_target w0; % 目标光束束腰可与源不同 % 计算LG模式的振幅分布未归一化 % 首先计算归一化径向坐标 rho sqrt(2)*r/w0_target rho sqrt(2) * r / w0_target; % 广义拉盖尔多项式 L_p^|l| (需要自定义或调用符号工具箱函数) % 这里以p1为例 L_1^|l|(x) (|l|1 - x) L (abs(l)1 - rho.^2); % 对于p1 % LG模式振幅 A_target_LG (rho.^abs(l)) .* L .* exp(-rho.^2/2); % 角向相位因子 exp(i*l*theta) 在GS算法中会被迭代出的相位替代这里只需振幅 A_target A_target_LG / max(max(A_target_LG)); % 归一化注意对于更通用的p需要实现或调用广义拉盖尔多项式的计算函数。Matlab符号数学工具箱有laguerreL函数但用于数值矩阵时需注意向量化操作。4.2 生成贝塞尔高斯光束% 参数定义 k_t 2*pi * 1e4; % 横向波矢决定贝塞尔光束的中心亮斑尺寸需根据需求调整 % 零阶贝塞尔高斯光束振幅分布 (J0) A_target_BG besselj(0, k_t * r) .* exp(-r.^2 / w0_target^2); A_target A_target_BG / max(max(A_target_BG)); % 归一化关键点k_t是一个关键参数k_t * r决定了贝塞尔函数的第一个零点位置从而控制了中心亮斑和环的间距。k_t越大中心斑越小环越密。需要根据模拟的物理尺寸r的范围合理选择避免采样不足导致失真。5. 常见问题、调试技巧与性能优化在实际编写和运行代码时你肯定会遇到各种问题。以下是我从多次实践中总结的“避坑指南”。5.1 算法不收敛或结果杂乱症状迭代误差不下降或最终目标光强图充满噪声与预期形状相差甚远。排查与解决检查传播模型确认正向和反向传播的传递函数H是否正确特别是fftshift和ifftshift的使用顺序。一个快速验证的方法是对一个点源delta函数进行传播看结果是否符合物理直觉如产生球面波。检查振幅约束确保A_source和A_target都是归一化后的非负实数矩阵。目标振幅不能有负值或复数。可视化它们确认分布符合预期。初始相位尝试如果使用全零平面波初始相位不收敛可以尝试使用随机相位初始化phase_initial 2*pi*rand(N);。随机相位能提供更多的自由度帮助算法跳出平庸解尤其对于中心对称的目标如空心光束有时效果更好。引入松弛因子如前所述采用输入输出算法。将源平面的更新规则改为U1_next U1 beta * (A_source.*exp(1i*phase_U1_prime) - U1);beta从0.5开始尝试。增加迭代次数有些复杂模式需要更多迭代才能收敛。将max_iter增加到500或1000。5.2 结果出现“孪生像”或共轭像症状在目标面除了想要的光斑在其对称位置出现一个类似的、较弱的“鬼影”。原因这是GS算法和傅里叶变换性质导致的常见问题。当源面和目标面的振幅约束都是实对称或厄米对称时算法求得的相位解可能不是唯一的会导致共轭像的出现。解决破坏对称性在源振幅或目标振幅中引入轻微的非对称性。例如将高斯光源的束腰中心稍微偏离计算网格中心。使用非对称初始相位使用随机初始相位而非零相位。采用加权GS算法在迭代初期对目标振幅约束进行松弛逐渐加强。5.3 采样与网格设置导致的失真症状生成的光束边缘出现锯齿、断裂或者贝塞尔光束的环不圆、不连续。排查奈奎斯特采样定理确保你的采样间隔dx足够小能够分辨光场中最高的空间频率。对于贝塞尔光束其空间频率与k_t相关。一个经验法则是dx lambda/(2*NA)其中NA是数值孔径的估计。如果出现混叠尝试增大N网格点数或减小L物理尺寸以降低dx。网格大小L确保计算窗口L足够大能够完整包含光束的主要能量部分避免能量在边界处被截断导致衍射效应干扰结果。可以观察初始高斯光束在网格边缘是否已衰减到接近零。传播距离z如果使用角谱传播距离z可以是任意的。但如果使用夫琅禾费近似U2 fft2(U1)则必须满足夫琅禾费条件z (N*dx)^2 / lambda。不满足时结果会错误。5.4 Matlab代码性能优化向量化与预计算所有操作都应使用矩阵运算避免循环。像H这样的传递函数应在循环外计算一次。内存管理对于大网格如N1024或更大复数矩阵会占用较多内存。使用single单精度而非默认的double双精度可以减半内存使用并在大多数情况下保持足够精度U1 single(U_source);。并行计算GS算法本身迭代间有依赖难以并行。但如果你需要为大量不同的目标光束计算相位图可以考虑使用parfor并行循环来处理不同的目标。收敛判断不必总是跑满max_iter次。可以在循环内加入判断当误差error(iter)的变化小于某个容差如1e-6时用break语句提前退出循环节省计算时间。我个人在实现这个项目时最大的体会是GS算法是一个“艺术多于科学”的过程。理论上很简单但调参松弛因子、迭代次数、初始条件对结果质量影响巨大。没有一套放之四海而皆准的参数。最好的方法是先从一个简单目标如将高斯光转换成另一个不同束腰的高斯光开始调试确保整个管道正确无误然后再挑战复杂的空心光或贝塞尔光。可视化每一个中间步骤如每次迭代后的目标面光强对于理解算法的行为非常有帮助。最后记住这个相位解是数值优化的结果它可能不是物理最优解但对于验证概念和指导初步实验已经足够强大。本文还有配套的精品资源点击获取