MATLAB高斯光束到平顶光束整形:GS算法与SLM相位分布计算

发布时间:2026/10/1 11:59:26
MATLAB高斯光束到平顶光束整形:GS算法与SLM相位分布计算 我一直觉得光束整形是光学实验里“看起来简单、做起来全是细节”的典型课题。这篇东西聊的是用MATLAB实现高斯光束到平顶光束的转变核心手段是GS算法和直接计算SLM相位分布这两条路。简单说就是激光器出来的光斑强度是中间亮、边缘暗的高斯分布但很多场景需要的是强度均匀的平顶光斑比如激光加工、荧光显微、光镊、全息显示。要完成这个转变最灵活的办法是用空间光调制器SLM加载特定的相位图而MATLAB就是用来算这张相位图的主力工具。这篇文章适合正在做光路实验的研究生、搞激光加工的工程师以及刚接触SLM不太清楚算法流程的新手。我会把GS算法的迭代逻辑、直接计算相位分布的几何映射思路、MATLAB的代码骨架、还有实际调试中容易踩的坑全部摊开讲尽量让你看完之后能在自己的平台上跑出第一版结果而不是只停留在概念层面。1. 先把问题拆清楚高斯光束凭什么要变成平顶1.1 高斯光束的“不均匀”到底意味着什么大多数激光器输出的基模光束横截面的光强分布都符合高斯函数中心最亮往外指数衰减。对某些应用来说这种分布是优点比如高斯光束聚焦后能得到很小的光斑适合精密加工和单分子激发。但对另外一些应用高斯分布反而是麻烦。举个例子激光焊接和退火工艺往往希望加热区域能量均匀否则材料中心已经过烧、边缘还没达到阈值。再比如荧光显微成像里如果激发光强度不均匀图像亮度就会有明显的中心亮斑后期校正很头疼。医学光动力治疗也是同理照射区域剂量不均会影响疗效。所以就有了“光束整形”这个方向想办法把高斯分布重新分配成接近矩形的平顶分布。传统做法有双折射透镜组、微透镜阵列、衍射光学元件但这些方案都是固定的结构做完了就很难改。SLM最大的价值在于可编程同一套光路换一张相位图就能切换成不同的输出形状甚至可以做动态扫描。这也是为什么近几年SLM几乎成了光学实验室的标配。1.2 SLM在这里扮演什么角色SLM本身不发光它是一个纯相位调制器件。液晶像素上加载不同的灰度值对应的相位延迟就不同。激光束经过SLM反射后波前的形状被改变再经过透镜在焦平面或者远场形成目标强度分布。关键在于我们肉眼看不到相位只能看到强度。所谓“计算SLM相位分布”本质上是求一个问题在输入面振幅已知高斯的情况下相位应该长什么样才能让光传播到输出面后振幅变成我们想要的平顶形状这个问题在光学里叫相位恢复问题。它不是直接测量出来的而是靠算法反推的。GS算法和几何映射法就是两种最常用的反推思路。在MATLAB里整个过程可以用几十行代码完成难的不是语法而是你选哪一种思路、以及参数和边界条件怎么定。1.3 两条技术路线的基本差异先说结论方便你心里有个谱GS算法是通过迭代在“输入面-输出面”之间来回做傅里叶变换反复施加振幅约束最终让相位收敛到一个能产生目标光强图案的解。直接计算法是基于几何光学的能量映射关系用一个解析或者半解析的公式直接把相位分布算出来不需要迭代。GS算法灵活任意形状的目标都能试但收敛后的相位图往往存在高频噪声和振铃。直接计算法速度快、相位光滑、对SLM像素量化更友好但只适合处理有明确对称性的整形目标比如高斯转平顶这种圆形分布。我在实际项目中两种方法都用过后面会详细讲各自的代码和适配场景。2. GS算法MATLAB里的核心迭代到底在迭代什么2.1 从傅里叶光学理解GS算法的物理直觉GS算法最早是Gerchberg和Saxton在1972年提出来的最初是为了解决电子显微学里的相位问题。后来被光学全息领域大量使用用来计算纯相位全息图。它的物理基础是夫琅禾费衍射。一束光经过SLM调制后再经过一个焦距为f的透镜在透镜后焦面上的复振幅分布近似等于SLM平面复振幅分布的傅里叶变换。用公式写就是E_out(u, v) F { E_in(x, y) * exp(iφ(x, y)) }这个关系意味着我们在SLM平面上改动相位输出面上的光强振幅的模平方就会跟着变。问题是这个变换是双向耦合的没法直接解析求逆所以GS算法采用“循环逼近”的策略。每次迭代做四件事输入面施加振幅约束光场幅度固定为高斯分布相位取当前迭代值。做一次正向傅里叶变换得到输出面的复振幅。输出面施加振幅约束保留计算出的相位但把幅度强行替换成目标平顶分布。做一次逆向傅里叶变换得到新的输入面相位然后回到第一步。听起来有点绕其实通俗说就是先猜一个相位看结果差多少然后把不匹配的那部分振幅“纠正”掉保留相位信息再反推反复纠正直到结果不再明显变化。2.2 一段可以直接跑的MATLAB代码骨架很多第一次接触GS算法的同学卡住的不是原理而是FFT的细节。下面这段代码是经过调通的版本基于二维FFT模拟SLM到焦平面的衍射过程。% GS.m - 高斯光束转平顶光束的基本实现 clear; close all; % 基本参数 lambda 532e-9; % 波长 532nm f 0.2; % 透镜焦距 200mm N 512; % 计算网格 pixel 8e-6; % SLM像素尺寸 8um % 坐标网格以SLM平面为准 x linspace(-N/2*pixel, N/2*pixel, N); [X, Y] meshgrid(x, x); r sqrt(X.^2 Y.^2); % 输入光场高斯振幅1/e半径 w0 2e-3; % 高斯光束半径 2mm A_gauss exp(-(r/w0).^2); % 振幅分布 % 目标光场平顶超高斯近似阶数越高边缘越锐利 R_target 3e-3; % 目标平顶半径 3mm n_order 8; A_target exp(-(r/R_target).^(2*n_order)); % 也可以用硬边界A_target double(r R_target); % 初始相位随机相位或者透镜相位 phase rand(N, N) * 2 * pi; % 随机相位作为起点 % GS迭代主循环 iter 100; q zeros(N, N); % 记录收敛误差 for k 1:iter % 输入面高斯振幅 当前相位 field_in A_gauss .* exp(1i * phase); % 正向傅里叶变换相当于透镜焦平面 spec fftshift(fft2(fftshift(field_in))); % 输出面施加强度约束保留相位替换振幅为目标平顶 spec_corrected A_target .* exp(1i * angle(spec)); % 逆向傅里叶变换回到SLM平面 field_back ifftshift(ifft2(ifftshift(spec_corrected))); % 提取新的相位 phase angle(field_back); % 计算误差输出面实际强度与目标强度的差值 diff abs(spec).^2 - A_target.^2; q(k) sqrt(mean(diff(:).^2)); end % 输出相位图wrap到0~2pi便于加载SLM phase_out mod(phase, 2*pi); figure; imagesc(x*1e3, x*1e3, phase_out); axis image; colormap gray; title(GS算法计算的SLM相位分布); xlabel(x (mm)); ylabel(y (mm));这段代码跑完之后你会得到一张黑白渐变的相位图。把这张图按灰度映射加载到SLM上光路中加上焦距匹配的透镜焦平面就能看到平顶光斑。2.3 三个最容易出错的细节第一fftshift不要乱删除。在没有做光轴离轴处理的情况下正确顺序是fftshift(fft2(fftshift(A)))。有些版本写成fft2(ifftshift(A))也可以但你要保证正向和逆向的shift操作是配套的。如果方向反了相位图会出现奇怪的高频条纹而且怎么调都调不好。第二振幅半径定义要统一。高斯光束的标准写法里光强分布是exp(-2r²/w²)振幅分布是exp(-r²/w²)。很多资料把w0定义混用导致仿真里高斯光斑比实际小或者大。建议统一采用振幅的1/e半径这样代入光路反推SLM上光束直径时更直接。第三迭代次数不是越多越好。我的经验是50到100次基本足够。超过200次以后相位图会开始出现明显的斑点噪声因为算法在强行匹配目标的每一个像素细节而这种细节在真实光路里根本达不到。如果你发现收敛曲线后期反而上升了说明迭代过拟合了可以把迭代次数降回去或者对目标分布加一点平滑。3. 直接计算SLM相位分布几何映射法的实现思路3.1 能量守恒是一切推导的起点第二种方法不走迭代而是直接从能量重新分配的角度去算。核心思路是输入面每个径向位置r1处的能量通过相位偏折被引导到输出面某个径向位置r2处。只要保证映射前后能量守恒输出强度就可以按需要分布。高斯光束的光强分布是I1(r1) (2P / πw0²) * exp(-2r1² / w0²)这里P是总光功率。目标平顶光强是I2(r2) P / (πR²)当r2 ≤ R其中R是平顶半径。能量守恒的意思是从r1处往外到r1dr1这个环带里的能量全部落进r2处对应的环带里。写成积分式∫0^r1 I1(s) * 2πs ds ∫0^r2 I2(s) * 2πs ds左边积出来是P * (1 - exp(-2r1²/w0²))右边是P * (r2² / R²)。两边一对比就能得到映射关系r2 R * sqrt(1 - exp(-2r1² / w0²))这个式子的物理图像很清晰r1等于0的光线正好落在输出面中心r1趋近无穷大的光线落在输出面边缘R处。高斯中心的能量密度高所以稍微往外移动一点就能填满目标中心附近的大面积区域高斯边缘能量少需要更大的偏折角度才能把剩余能量铺开。相位偏折的角度和相位梯度之间的关系是dφ/dr1 k * r2 / f。这里的k 2π/λ是波数f是透镜焦距。也就是说要知道相位φ(r1)不是直接套公式而是需要对r2进行径向积分。3.2 MATLAB实现几何映射法的代码由于映射积分没有解析解直接用数值积分处理。下面的代码演示了完整流程% GeometricalMapping.m - 几何映射法计算平顶整形相位 clear; close all; % 参数设置 lambda 532e-9; f 0.2; N 512; pixel 8e-6; w0 2e-3; % 高斯光束半径 R_target 3e-3; % 平顶目标半径 % 一维径向坐标 r1 linspace(0, N/2*pixel, N/2); dr r1(2) - r1(1); % 高斯能量累积 integral 1 - exp(-2 * r1.^2 / w0^2); integral max(integral, 0); % 映射关系r1 - r2 r2 R_target * sqrt(integral); % 相位梯度 k 2 * pi / lambda; dphi_dr k * r2 / f; % 数值积分得到径向相位分布 phi_r cumtrapz(r1, dphi_dr); % 扩展到二维 x linspace(-N/2*pixel, N/2*pixel, N); [X, Y] meshgrid(x, x); r_2d sqrt(X.^2 Y.^2); phi_2d interp1(r1, phi_r, r_2d(:), linear, phi_r(end)); phi_2d reshape(phi_2d, N, N); % 相位wrap到0~2pi phase_out mod(phi_2d, 2*pi); % 模拟输出光斑 input_field exp(-(r_2d/w0).^2) .* exp(1i * phi_2d); output_field fftshift(fft2(fftshift(input_field))); intensity abs(output_field).^2; % 绘制相位和模拟输出 figure; subplot(1,2,1); imagesc(x*1e3, x*1e3, phase_out); axis image; colormap gray; title(几何映射法相位分布); subplot(1,2,2); imagesc(x*1e3, x*1e3, intensity); axis image; title(模拟输出光强);注意几个实现细节。第一interp1里我指定的末尾值是phi_r(end)这样r1超出范围时相位不再变化避免边缘出现不合理的跃变。第二中心处r2在r10时等于零但相位梯度的导数在中心附近偏大量化到SLM灰度后中心区域会有比较密集的条纹这是正常的实验中这些细节会被像素尺寸平滑掉。3.3 两种方法到底该选哪个这里放一张我做过的对比经验表纯属个人实践感觉不是绝对的结论对比维度GS算法几何映射法计算时间迭代100次约几十秒到几分钟一次积分秒出结果目标任意性高任意图案都支持低适合圆形对称整形相位光滑度有噪声和振铃光滑连续对SLM友好输出边缘锐度可调容易过冲边缘过渡平滑实验鲁棒性对系统误差敏感相对更稳定代码复杂度中等低我自己的筛选原则是如果只是做高斯转平顶优先用几何映射法因为相位光滑、加载到液晶SLM上量化误差小而且不需要微调迭代参数。如果需要生成任意形状比如矩形平顶、环形光斑、字母形状的衍射图案GS算法是更通用的选择。很多项目最终会两条路都走先用几何映射法快速验证实验系统再用GS算法优化特殊目标。4. 从仿真相位到真实光路参数设计的完整检查清单4.1 仿真参数的标定逻辑很多人在MATLAB里调出一张漂亮的相位图一到实验就完全不是那么回事大部分原因是参数标定不一致。最核心的一个参数是仿真网格尺寸与SLM物理尺寸的对应关系。假设SLM像素尺寸是8μm有效区域是512×512那你仿真时网格就应该设置为512×512每个网格对应一个SLM像素。坐标向量x的范围应当是从-2.048mm到2.048mm而不是随便选一个“看着舒服”的归一化坐标。还有一个非常容易忽略的点FFT之后输出面的像素间隔不等于输入面的像素间隔而是满足δx_out λ * f / (N * δx_in)代入上面的参数λ532nmf200mmN512δx_in8μm算出来δx_out大约是26μm。也就是说输出面上的图案物理尺寸会被放大或者缩小取决于焦距和像素尺寸的组合。你在看模拟结果的平顶半径时不要直接拿像素数当毫米数要乘以这个缩放系数。这一步错掉加载到SLM上的相位图对应的平顶尺寸就会偏得离谱。4.2 相位图加载到SLM时的灰度映射SLM的控制软件一般接受8位或10位灰度图灰度值0到255对应相位0到2π。问题在于不同厂家SLM的相位响应曲线并不是理想的线性关系而且灰度变大对应的相位可能增加也可能减小。实操上建议按这几个步骤走先查SLM的出厂相位标定曲线或者自己搭一个迈克尔逊干涉仪标定。在MATLAB里生成相位图后用mat2gray(phase_out, [0, 2*pi])转成0到1的归一化灰度再乘255得到8位灰度。如果发现实验光斑是“反相”的或者中心出现暗斑而边缘亮大概率是灰度方向反了把灰度图取反再试一次。这一步千万别偷懒直接用屏幕截图去加载相位图必须按原始矩阵导出成bmp或png且不能经过任何压缩和缩放。微信传图、QQ截图这类操作基本都会破坏灰度细节加载上去相位就完全不是那么回事了。4.3 光学链路中最容易出问题的三个位置SLM前端的偏振态非常关键。绝大多数液晶SLM只对特定偏振方向的光有相位调制作用所以SLM前面一定要加偏振片或半波片把激光偏振调到SLM要求的轴方向。我第一次做的时候没仔细看手册结果SLM反射率倒是很高但相位调制几乎为零折腾了一天才发现是偏振方向差了45度。SLM到透镜的距离也要注意。脚本里的傅里叶变换默认的是“SLM平面到透镜后焦面”的关系实际上透镜放在SLM之后、与SLM距离越近越好。如果两者间距过大实际光路会引入额外的菲涅耳衍射传播相位图的整形效果就会变差。稳妥的做法是采用4f系统把SLM放在第一个透镜的前焦面目标面放在第二个透镜的后焦面这样傅里叶关系最干净。还有个容易被忽略的是透镜口径。计算相位图时很多人在仿真里把网格边界设置为振幅为零的区域觉得没有影响。但实际透镜如果口径不够大边缘的高频相位成分会被切掉输出平顶的边缘就会出现明显的强度抖动。选择透镜时通光口径至少要是仿真区域尺寸的1.2倍以上才安全。5. 实际调试中我踩过的那些坑5.1 输出面中心总有一个刺眼的亮斑这个问题十个用SLM的人九个会遇到。原因是SLM的相位调制不是理想的“纯相位”液晶分子对光还有残余的振幅调制加上像素间有间隙部分光没有被调制直接沿光轴方向透射或反射在透镜焦平面形成一个零级亮斑。处理思路有两个方向。一是尽量压低零级把SLM放在偏振干涉模式下通过调整相位灰度范围避免残余调制二是光学上把零级光滤掉给相位图叠加一个离轴的闪耀光栅相位让目标光斑偏到一边然后用挡板或光阑把中心的零级挡住。离轴相位可以用mod(k_x * X k_y * Y, 2*pi)叠加到已有相位图上。代价是SLM的有效分辨率会损失一点但实验效果好很多。5.2 平顶的边缘总有一圈一圈的振铃振铃的根源是目标平顶的边界太锐利。仿真里的“硬边界”平顶在频域上带宽很宽而SLM的像素有限相当于低通滤波自然会出现吉布斯现象也就是边缘过冲和振荡。最简单的解决办法是把目标分布改得平滑一点用超高斯代替硬边界。我一般在GS算法中使用阶数8左右的超高斯边缘的过渡带大约占平顶半径的10%到20%。如果你对边缘锐度要求特别高可以试试增加网格分辨率或者改用混合算法比如把GS迭代的结果作为初始值再跑几轮模拟退火或梯度优化来压制振铃。5.3 仿真里平顶很均匀实验里却是毛毛糙糙的仿真用的是理想单色平面波但真实激光器输出不是完美高斯可能有像散、有灰尘衍射、有杂散光。SLM的相位量化误差和像素填充率也会给结果增加噪声。这种问题没有一步到位的解药我的习惯是先分步排查不用SLM直接在焦平面看激光光斑确认原始光束质量。加载一个纯透镜相位也就是焦点相位图看焦点是否锐利对称检查SLM本身有没有坏区或者响应不均匀。加载优化好的整形相位观察整体包络和理论模拟的差异。如果原始光束质量不够通常要在SLM前面加太空滤波器和准直系统把光束先整干净。5.4 相位图上微小差异对输出影响巨大GS算法算出的相位图迭代50次和迭代300次之间肉眼看起来差异可能不大但输出强度分布会有明显变化。几何映射法对r2映射积分时累积误差带来的相位偏差同样会影响边缘过渡带的宽度。我的建议是把仿真当成“相对值”而不是“最终值”。相位图的绝对形状并不重要重要的是算法得到的相位梯度是否平滑、输出光斑是否满足应用需求。到了实验阶段真正的微调手段其实是灰度映射范围的选择很多情况下稍微改变相位调制深度比如把相位范围压缩到1.8π而不是2π平顶均匀度会有意想不到的提升这一点需要在实际光路上耐心扫参数。结尾一点个人体会我最初做这个项目时一上来就迷信GS算法总觉得迭代算法“高级、万能”结果在实验里被零级亮斑和振铃折磨了很久。后来老老实实把几何映射法也用起来才发现不同算法不是替代关系而是互补关系几何映射快速稳定GS算法处理特殊目标两者的仿真结果还能互相验证。如果你正在做类似的整形项目我强烈建议先用几何映射法跑通整条光路再根据需求决定要不要上GS算法优化。最后提醒一句所有仿真结果都要换算到真实物理尺寸再上光路验证这一步省不掉也最容易出错。