简介：这份资源面向光学工程、光电技术及LED封装应用方向的本科与硕士阶段学习者，聚焦LED照明系统中的透镜二次光学设计问题，提供一套可直接运行的Matlab实现方案。包内共7个文件，以m主程序、py脚本、txt结果与说明文档为主，另含少量工程配置文件，压缩包约5KB，体积轻便便于快速获取。代码适配matlab2014a、2019a、2021a等版本，基于光线追迹或非成像光学原理，完成自由曲面或旋转对称透镜的光学参数建模与光强分布优化，输出照度分布图与光效评估数据等典型结果。已有63人学习关注，适合用于课程设计、毕业设计或基础科研验证。下载后无需额外配置即可运行查看透镜设计效果，帮助读者理解配光需求中的建模与仿真流程，掌握从参数设定到结果评估的完整思路，是入门二次光学设计的实用参考。
1. LED透镜二次光学设计：从配光曲线到Matlab可复现的完整链路
拿到一颗LED光源，裸芯的朗伯体配光又圆又散，照到目标面上要么中间过曝、要么边缘发暗，路灯、投光灯、洗墙灯几乎没法直接用。二次光学设计要干的事，就是在LED和被测面之间加一片透镜或反光杯，把光重新分配到你需要的那块区域里。Matlab在这个环节的价值不是画个漂亮三维图，而是把「光源模型→透镜面型→追迹结果」串成一条能反复迭代的数值链路。这篇笔记面向两类人：一类是光学/照明方向的在校生，想搞明白自由曲面透镜的迭代逻辑；另一类是做灯具结构的工程师，手头有Matlab但不知道怎么把配光需求翻译成可算的面型。整条链路我会按「先立模型、再写迭代、最后验证」的顺序拆开，代码能直接跑，参数能直接改。
2. 二次光学设计到底在算什么：光源模型与目标配光的数学描述
二次光学设计的本质是一个逆问题：已知光源的发光分布，已知目标面上的照度或强度分布，反求中间那片透镜的入射面和出射面形状。听起来像解方程，实际上大部分工程做法是「正向迭代逼近」——先猜一个面型，追迹一遍看结果，根据偏差修正面型，再追迹，直到配光曲线落在容差带里。要跑通这条链路，第一步不是写透镜，而是把光源和目标都写成Matlab能算的数学对象。
2.1 朗伯光源的采样模型：为什么不能只用一个点光源
很多新手第一版代码直接拿一个点光源加余弦辐射角去算，结果做出来的透镜装上去边缘总是差一截。原因在于真实LED芯片是有发光面积的，常见1mm×1mm甚至更大，芯片边缘和中心到透镜的距离不同，入射角分布也不同。工程上够用的做法是把芯片离散成N×N个网格点，每个点按朗伯体发射光线，权重按面积分配。
% 朗伯光源离散采样：芯片尺寸1mm x 1mm，网格20x20 chipSize = 1e-3; % 芯片边长，单位米 N = 20; % 每边采样点数 [xg, yg] = meshgrid(linspace(-chipSize/2, chipSize/2, N)); srcPoints = [xg(:), yg(:), zeros(N*N,1)]; % 发光面在z=0平面 dA = (chipSize/(N-1))^2; % 每个采样点代表的面积 weights = dA * ones(N*N,1); % 面积权重，总光通量归一化用 % 朗伯体方向采样：极角theta按cos分布，方位角phi均匀 nRays = 200; % 每个发光点的光线数 theta_max = pi/2; % 半球发射 u = rand(N*N*nRays,1); theta = acos(1 - u*(1-cos(theta_max))); % 逆变换采样，保证cos权重 phi = 2*pi*rand(N*N*nRays,1);这段代码的关键在theta = acos(1 - u*(1-cos(theta_max)))这一行。朗伯体的辐射强度正比于cosθ，直接均匀采样θ会高估大角度光线，逆变换采样才能让每条光线携带相同的能量权重。weights按面积分配，是为了后面统计目标面照度时能正确加权。如果你只做定性验证，N取10、nRays取50就够；要做定量配光，N至少20、nRays至少200，否则目标面会出现明显的采样噪点，看起来像透镜面型有问题，其实是光源采样不够。
2.2 目标配光的两种表达：照度均匀 vs 强度分布
目标面要什么，直接决定透镜面型往哪个方向偏。常见需求分两类：一类是照度均匀，比如投影仪照明、显微镜照明，要求目标面每个点的照度偏差小于某个百分比；另一类是特定强度分布，比如路灯要求蝙蝠翼配光，投光灯要求矩形光斑。Matlab里描述目标分布，我一般用目标面上的网格加权重矩阵。
% 目标面定义：距离光源H米，矩形区域WxL H = 5; % 照明距离，单位米 W = 4; L = 6; % 目标区域宽、长，单位米 Mx = 60; My = 80; % 目标面网格 [xT, yT] = meshgrid(linspace(-W/2, W/2, Mx), linspace(-L/2, L/2, My)); % 目标照度分布：这里以矩形均匀光斑为例 E_target = ones(My, Mx); % 如果是蝙蝠翼配光，可以改成沿x方向的双峰函数 % E_target = exp(-((abs(xT)-W/4).^2)/(2*(W/8)^2)); E_target = E_target / sum(E_target(:)); % 归一化，总能量为1E_target归一化这一步很重要。光源采样出来的总能量也是归一化的，两边能量口径一致，后面算偏差才有意义。如果你做的是强度分布而不是照度，把E_target换成对应角度的强度权重即可，追迹时统计的是光线方向而不是落点位置。这里没有唯一正确的写法，关键是目标分布要能写成矩阵，后面才能和追迹结果做逐点比较。
2.3 自由曲面透镜的参数化：为什么用多项式而不是直接优化每个点
理论上透镜面型可以逐点自由，但那样未知数太多，迭代容易发散。工程上常见做法是把出射面写成多项式或样条，用少量系数控制整体形状。我一般用径向多项式，因为大部分LED透镜是旋转对称或近似对称的，系数少、收敛快。
% 出射面径向多项式：z(r) = sum(c_k * r^(2k)), k=1..K K = 6; % 多项式项数 c = zeros(K,1); % 待优化系数，初值全零即平面 c(1) = -0.02; % 给一个初始曲率，避免完全平面导致全反射 % 面型函数句柄 lensSurface = @(r) sum(c .* (r.^(2*(1:K)'))');c(1)对应r²项，控制近轴曲率；c(2)对应r⁴项，控制边缘偏折；更高阶项用来修局部偏差。初值给一个负曲率是为了让光线在出射面发生会聚而不是发散，如果全零平面，大角度光线会在出射面发生全反射，追迹结果一片黑，新手很容易在这里翻车。K取6是经验值，太少修不动边缘，太多容易过拟合出波浪面型，加工不出来。
3. 用Matlab跑通第一次追迹：从光线发射到目标面照度统计
模型立好之后，下一步是让光线从光源出发，经过透镜入射面折射、出射面折射，最后落到目标面，统计落点分布。这一章把追迹循环写清楚，包括折射计算、全反射判断和照度累加。代码跑通一次，你就能看到当前面型对应的配光长什么样，后面所有迭代都建立在这个循环上。
3.1 入射面与出射面的折射计算：矢量形式的斯涅尔定律
折射计算是追迹的核心，用矢量形式写比角度形式更稳，因为不用反复算三角函数。给定入射方向d、表面法线n、折射率比eta = n1/n2，折射方向t可以直接算出来。
function t = refractVector(d, n, eta) % d: 入射单位矢量, n: 表面单位法线(指向入射侧), eta = n1/n2 d = d / norm(d); n = n / norm(n); cosi = -dot(d, n); if cosi < 0 n = -n; cosi = -cosi; % 保证法线朝向入射侧 end k = 1 - eta^2 * (1 - cosi^2); if k < 0 t = [0 0 0]; % 全反射，返回零矢量标记 return; end t = eta * d + (eta * cosi - sqrt(k)) * n; t = t / norm(t); endcosi是入射角余弦，k < 0对应全反射，实际透镜设计里要尽量避免，因为全反射光线不受出射面控制，会变成杂散光。eta的取值：如果光从空气进入透镜，eta = 1/1.59（常见光学塑料折射率）；从透镜出射到空气，eta = 1.59/1。这个函数被调用两次，一次在入射面，一次在出射面。
3.2 主追迹循环：每条光线的完整路径
主循环把光源采样、两次折射、目标面求交串起来。为了速度，Matlab里尽量向量化，但第一版用for循环更直观，跑通了再优化。
nLens = 1.59; % 透镜折射率 nAir = 1.0; hits = zeros(My, Mx); % 目标面照度累加矩阵 for i = 1:size(srcPoints,1) p0 = srcPoints(i,:); for j = 1:nRays idx = (i-1)*nRays + j; d = [sin(theta(idx))*cos(phi(idx)), ... sin(theta(idx))*sin(phi(idx)), ... cos(theta(idx))]; % 入射面：假设为平面z=0.005，法线[0 0 -1] p1 = p0 + d * (0.005 - p0(3)) / d(3); t1 = refractVector(d, [0 0 -1], nAir/nLens); if all(t1==0), continue; end % 出射面：用多项式面型求交，这里简化用z=0.015平面示意 p2 = p1 + t1 * (0.015 - p1(3)) / t1(3); % 出射面法线由面型梯度决定，示意用平面法线 t2 = refractVector(t1, [0 0 -1], nLens/nAir); if all(t2==0), continue; end % 目标面求交：z = H s = (H - p2(3)) / t2(3); if s <= 0, continue; end pHit = p2 + t2 * s; % 累加到目标面网格 ix = round((pHit(1)+W/2)/W*(Mx-1)) + 1; iy = round((pHit(2)+L/2)/L*(My-1)) + 1; if ix>=1 && ix<=Mx && iy>=1 && iy<=My hits(iy, ix) = hits(iy, ix) + weights(i); end end end这段代码里入射面和出射面都简化成了平面，目的是先把循环跑通。实际设计时，入射面通常是球面或平面，出射面才是自由曲面，求交要用数值方法解z_lens(r) = p_z + t_z * s。s <= 0的判断是防止光线往反方向跑，ix/iy的边界检查是防止落点出界导致索引报错。跑完一遍，hits就是当前面型下的目标面照度分布，可以直接imagesc(hits)看形状。
3.3 配光评价指标：均匀度、光效和偏差怎么算
有了hits，下一步是量化评价。最常用的三个指标：均匀度、光效和相对偏差。
E = hits / sum(hits(:)); % 归一化照度 uniformity = min(E(E>0)) / mean(E(E>0)); % 最小/平均，越接近1越均匀 efficiency = sum(hits(:)) / sum(weights) * size(srcPoints,1); % 简化光效 deviation = norm(E(:) - E_target(:)) / norm(E_target(:)); % 相对偏差uniformity用最小比平均，比最大比最小更稳，因为个别噪点不会把指标拉爆。efficiency这里做了简化，严格算要考虑每条光线的能量权重和全反射损失。deviation是后面迭代的目标函数，越小说明当前配光越接近目标。这三个指标每次迭代都打印出来，你就能看到面型系数往哪个方向调是有效的。
4. 面型迭代与参数调优：让配光曲线收敛到目标带
第一次追迹结果大概率不达标，均匀度可能只有0.3，偏差0.5以上。这一章讲怎么根据偏差修正多项式系数，以及哪些参数最敏感、哪些调了也白调。迭代不是盲目试，每一步都要有依据。
4.1 阻尼最小二乘迭代：系数更新的具体写法
把偏差deviation看成系数c的函数，用数值差分算雅可比，然后做阻尼最小二乘更新。Matlab里不用工具箱也能写。
lambda = 0.01; % 阻尼因子 maxIter = 30; tol = 1e-3; for iter = 1:maxIter [E, dev] = traceAndEvaluate(c); % 封装好的追迹+评价函数 if dev < tol, break; end J = zeros(numel(E), K); for k = 1:K dc = zeros(K,1); dc(k) = 1e-6; [~, dev_k] = traceAndEvaluate(c + dc); J(:,k) = (E - E) / 1e-6; % 实际应存E_k - E，这里示意 end % 阻尼最小二乘更新 delta = (J'*J + lambda*eye(K)) \ (J'*(E_target(:) - E(:))); c = c + delta; fprintf('iter %d, dev=%.4f\n', iter, dev); endlambda是阻尼因子，越大更新越保守，越小收敛越快但容易震荡。我一般从0.01起步，如果dev连续两次不降就乘2，降了就除1.5。1e-6的差分步长对多项式系数够用，太小会被浮点误差淹没，太大线性近似不准。maxIter取30是经验值，大部分情况10到15次就收敛了，超过30还不降说明面型参数化不够或者目标本身不可达。
4.2 三个最敏感的参数：曲率、高阶项和透镜厚度
迭代过程中你会发现，不是所有系数都一样敏感。c(1)曲率影响整体光斑大小，调它光斑缩放最明显；c(2)和c(3)影响边缘光线偏折，调它们均匀度变化最快；c(5)、c(6)高阶项只影响局部，调了容易出波浪，除非前几项已经收敛得差不多，否则不建议动。透镜厚度也是敏感参数，太薄边缘光线入射角大，全反射损失高；太厚材料吸收和成本都上去。常见做法是厚度取口径的0.3到0.5倍，先定厚度再优化面型。
4.3 收敛判据与停止条件：什么时候算调好了
收敛判据不能只看deviation，还要看均匀度和光效。我一般设三个条件同时满足才停：deviation < 0.05、uniformity > 0.7、efficiency > 0.85。如果deviation降了但uniformity反而变差，说明迭代在牺牲均匀性换整体偏差，这时候要检查目标分布是不是本身就不合理。另一个停止条件是系数变化量norm(delta) < 1e-8，说明已经到数值精度极限，再迭代也没意义。
5. 避坑与排查：二次光学设计里最容易翻车的五个地方
这一章是我自己踩过的坑，每个都按「现象→原因→解决」写。新手照着做能省不少时间，熟手可以对照检查自己的代码。
5.1 目标面一片黑或者只有中心一个点
现象：追迹跑完，hits矩阵几乎全是零，只有中心几个像素有值。原因通常是出射面法线方向搞反了，折射计算把光线折回了透镜内部，或者全反射判断把大部分光线标记为零矢量直接跳过。解决：先检查refractVector里法线是不是指向入射侧，再检查eta有没有写反。一个快速验证方法：把出射面暂时改成平面，如果平面能出光斑，说明面型法线有问题；如果平面也黑，说明折射函数本身有bug。
5.2 配光曲线边缘出现明显锯齿或台阶
现象：imagesc(hits)看目标面，光斑边缘不是平滑过渡，而是一圈圈台阶。原因一般是目标面网格太粗，或者光源采样点太少，落点统计时多个光线挤进同一个像素。解决：先把Mx、My加倍，如果锯齿变细但没消失，再把N和nRays加倍。如果还不行，检查ix、iy的取整方式，round在边界处会跳变，可以改成双线性累加，把一条光线的能量按距离分到相邻四个像素。
5.3 迭代不收敛，偏差来回震荡
现象：deviation在0.1附近上下跳，不往下走。原因通常是阻尼因子太小，或者多项式项数太多导致雅可比病态。解决：先把lambda乘10，看震荡是否减弱；如果减弱了但收敛慢，再逐步降回来。另一个办法是减少K，从6降到4，先让低阶项收敛，再逐步加高阶项。我一般用「先粗后细」策略：K=3跑10次，K=5跑10次，K=6跑10次，比一上来就K=6稳得多。
5.4 光效突然掉到0.5以下
现象：均匀度还行，但efficiency很低，大量光线没到目标面。原因通常是透镜边缘太陡，大角度光线在出射面发生全反射。解决：检查出射面边缘斜率，如果dz/dr过大，把c(2)或c(3)的绝对值调小，让边缘平缓一些。另一个原因是透镜口径不够，部分光线根本没打到透镜上，这时候要加大口径或者减小光源到透镜的距离。
5.5 换一个LED型号结果全变了
现象：同一套代码，换一颗芯片尺寸或发光角不同的LED，配光完全不对。原因是你之前的面型是针对特定光源过拟合的，光源模型一变，入射角分布变了，原来的面型就不适用。解决：把光源模型参数化，芯片尺寸、发光角、折射率都做成变量，换LED时重新跑一遍迭代，不要指望一套面型通吃。工程上常见做法是针对每个LED型号单独优化，或者设计时留一定容差，让面型对光源变化不敏感。
6. 从能跑到好用：把追迹速度提上来并做一次交叉验证
第一版代码能跑通，但速度往往让人抓狂。20×20光源点、每点200条光线、30次迭代，for循环跑下来可能要几十分钟。这一章讲两个提速技巧和一个验证方法，让整个链路从「能跑」变成「好用」。
6.1 向量化追迹：把内层循环压成矩阵运算
Matlab的for循环慢，但矩阵运算快。把每个光源点的所有光线打包成矩阵，一次算完折射，速度能提一个数量级。
% 向量化折射：一次处理所有光线 function T = refractVectorBatch(D, N, eta) % D: n x 3 入射方向, N: 1 x 3 法线 D = D ./ vecnorm(D,2,2); N = N / norm(N); cosi = -D * N'; flip = cosi < 0; Nrep = repmat(N, size(D,1), 1); Nrep(flip,:) = -Nrep(flip,:); cosi(flip) = -cosi(flip); k = 1 - eta^2 * (1 - cosi.^2); T = eta * D + (eta * cosi - sqrt(max(k,0))) .* Nrep; T(k<0,:) = 0; % 全反射标记 T = T ./ vecnorm(T,2,2); endvecnorm按行求模，max(k,0)防止负数开方报错，k<0的行直接置零。这个函数一次处理所有光线，比逐条调用快很多。注意Nrep要按flip翻转，保证法线始终朝向入射侧。向量化之后，同样的追迹量从几分钟降到几秒，迭代30次也能接受。
6.2 用能量守恒做一次独立校验
追迹结果对不对，不能只看图好看。一个独立校验方法是能量守恒：所有到达目标面的光线能量加上全反射损失，应该等于光源发出的总能量。
totalEmitted = sum(weights) * nRays; totalHit = sum(hits(:)); totalLost = totalEmitted - totalHit; fprintf(' emitted=%.4f, hit=%.4f, lost=%.4f, loss ratio=%.2f%%\n', ... totalEmitted, totalHit, totalLost, totalLost/totalEmitted*100);如果loss ratio超过15%，说明全反射或口径损失太大，配光再好也不实用。如果loss ratio是负数，说明统计有重复累加，检查hits累加时有没有把同一条光线算两次。这个校验不依赖任何光学软件，纯靠能量账本，能抓出大部分低级错误。
6.3 和商用软件对一次：什么情况下值得做
Matlab追迹再快，也只是近似。如果项目要出正式配光报告，最终还是要用商用光学软件验证。但什么时候值得花时间做交叉验证？我的经验是：如果Matlab算出来的均匀度和商用软件差在5%以内，可以直接用Matlab做迭代优化，最后用商用软件出报告；如果差超过10%，先查光源模型和折射率设置是不是一致，大概率是某个参数口径没对齐。交叉验证不用每次都做，面型定型前做一次就够。
6.4 一个我常犯的错误：忽略透镜材料的色散
最后说一个我踩过的坑。上面所有代码都假设折射率是常数1.59，但实际光学塑料在不同波长下折射率差0.01左右，白光LED追迹时边缘会出现色差，目标面光斑边缘带彩边。解决方法是把折射率写成波长的函数，对RGB三个波长分别追迹再叠加。如果只是做照度均匀性设计，色散影响不大；但如果做成像或对颜色敏感的场景，这一步不能省。我现在的习惯是：第一版用单折射率快速迭代，面型定型后再加色散验证一遍，避免一开始就陷入多波长优化的复杂度里。
希望帮到你。
本文还有配套的精品资源，点击获取
方法：无锁读取)
