简介一套基于MATLAB编写的光线追迹计算程序面向太阳能热发电、光伏与光学系统设计人员可模拟光线传播路径分析反射、折射等行为辅助优化集热器或光伏面板布局。压缩包共含7个文件其中5个m脚本为核心计算模块覆盖主流程、反射/折射计算、光线与表面交点求解等功能另附README说明文档与许可证便于理解代码结构和合规使用。整包仅6KB轻巧易用目录结构清晰适合具备一定光学或MATLAB基础的读者直接运行与二次开发。目前已有566人学习下载。通过该程序使用者可直观观察光线追迹过程自由调整折射率、反射率等参数并理解Chebyshev多项式在光线路径近似中的加速作用从而掌握几何光学算法要点为太阳能热发电或光伏系统的效率评估提供高效实用工具。1. 光线追迹程序不是黑匣子先搞清楚它在算什么光线追迹程序的核心不是“画图”而是用海量光线样本去模拟光在场景中的反射、折射和能量分配。我做光学系统像质评估时最早就用现成渲染器出图但它只能给我漂亮的图片给不了焦平面上精确的坐标分布和能量统计。只有自己写追迹程序才能看到每一条光线到底从哪里来、往哪里去、带了多重的一份能量。这篇笔记记录我从零实现一个最小光线追迹程序的过程包含射线生成、几何求交、反射与折射、能量累积和解析验证的完整路径。如果你正在做光学仿真、辐射传热模拟或者需要对渲染器做底层参数控制下文会按我的实现顺序讲清楚。新手可以直接照抄第3章的代码熟手建议优先看第4章的参数标定和第5章的5个坑那几处都是我实际调程序时踩过的。2. 从物理到代码光线追迹的数学模型与程序选型2.1 光线的数学表示参数方程、方向归一化与有效区间几何光学里的光线可以用一条射线表示起点 O方向 D写成参数方程 P(t)OtD其中 t 表示沿 D 方向传播的距离。写代码前必须定死两个约定D 必须归一化t 必须大于等于 0。归一化之后t 才真正对应光线走过的路径长度平面求交公式里的分母 N·D、折射公式里的角度余弦才有正确的物理含义t 为负表示交点在起点后方追迹时应该直接丢弃。我习惯把方向归一化和射线有效性写进数据类让所有外部传入的方向都自动规范化。这样不仅减少调用方的负担也避免在求交、折射这些地方反复检查向量长度。from dataclasses import dataclass import numpy as np dataclass class Ray: origin: np.ndarray direction: np.ndarray def __post_init__(self): self.origin np.asarray(self.origin, dtypefloat) self.direction np.asarray(self.direction, dtypefloat) norm np.linalg.norm(self.direction) if norm 1e-12: raise ValueError(光线方向不能是零向量) self.direction / norm逻辑说明__post_init__是 dataclass 初始化后自动执行的钩子我在里面完成类型转换和归一化。参数说明origin是三维浮点坐标direction传入任意非零向量构造后变成单位向量。这样做的好处是后续所有函数拿到 Ray 时都默认方向长度为 1求交公式可以少做很多次np.linalg.norm。2.2 三种必写的求交与方向变换公式追迹循环里调用最频繁的操作就是求交。先写平面和球面这两个最基础的几何体再加一个光学仿真常用的旋转抛物面。平面求交平面过点 P0法线 N交点满足 N·(P-P0)0代入射线方程得到 t N·(P0-O)/(N·D)。当 N·D 的绝对值接近 0 时光线与平面平行视为不相交。球面求交把射线方程代入球面方程 |OtD-C|²R²展开后是关于 t 的一元二次方程。判别式小于 0 无交点否则取两个根中较小的非负根。抛物面求交旋转抛物面方程写为 z(x²y²)/(4f)代入射线方程后同样得到一个二次方程。这类面在平行光聚焦、望远镜副镜这类场景里能解析求交完全不需要网格化。反射和折射是方向变换的核心。反射方向公式为 R D - 2(D·N)NN 和 D 都要求是单位向量。折射方向用斯涅尔定律的矢量形式设入射介质折射率 n1、出射介质折射率 n2N 指向入射介质一侧cosθi-D·Nηn1/n2则出射方向 T ηD (ηcosθi - √(1-η²(1-cosθi²)))N。根号内部小于 0 表示发生全反射此时应返回空。def intersect_plane(ray: Ray, p0, normal): denom np.dot(normal, ray.direction) if abs(denom) 1e-12: return None t np.dot(normal, p0 - ray.origin) / denom return t if t 0 else None def intersect_sphere(ray: Ray, center, radius): oc ray.origin - center a 1.0 # 方向已经归一化 b 2.0 * np.dot(oc, ray.direction) c np.dot(oc, oc) - radius * radius disc b * b - 4.0 * a * c if disc 0: return None t (-b - np.sqrt(disc)) / (2.0 * a) return t if t 0 else None def reflect_dir(d, n): return d - 2.0 * np.dot(d, n) * n def refract_dir(d, n, n1, n2): cosi -np.dot(d, n) eta n1 / n2 sin2t eta * eta * (1.0 - cosi * cosi) if sin2t 1.0: return None return eta * d (eta * cosi - np.sqrt(1.0 - sin2t)) * n逻辑说明球面求交里因为 Ray 方向已归一化二次项系数 a 直接等于 1省一次乘法。平面求交里denom过小说明光线几乎平行于平面返回 None 避免除零。折射方向里sin2t 1.0对应全反射调用处必须把它当作镜面反射处理而不是继续追迹折射光线。参数说明n1、n2是两种介质的折射率n必须指向入射介质一侧如果你的场景法线统一朝外从空气进入玻璃前要先将法线取反。2.3 自己写还是用现成库选型边界这个程序要不要自己写取决于你手里拿到的到底是“面型数据”还是“物理假设”。如果只是要一套能跑复杂 CAD 模型的非序列追迹直接上商业光学设计软件或专门的光线追迹库更省事它们对镜头库、镀膜数据、像差报表的支持不是手写几十行能用完的。但如果要往结果里嵌入自定义散射模型比如测量过的粗糙度分布、镀膜的波长响应或者要输出每一条光线的完整路径用于自研算法商业软件往往不够透明这时候必须自己写核心循环。我的一般做法是先手写一个最小循环跑通一个解析解可验证的场景例如第 6 章的抛物面聚焦。验证通过后再决定是否迁移到现成库。直接上手现成库的年轻人最容易遇到的问题是不知道某个像差值为什么算出来是负的也不知道网格设置怎么影响收敛。手写一次很多“玄学”就变成了算术。3. 用 Python 写一个最小可跑的光线追迹程序3.1 场景数据结构Ray、几何体与 Scene有了 Ray 之后下一步是定义几何体对象。我定义一个 Geometry 基类要求所有子类必须实现intersect和normal两个方法前者返回光线与表面的最近交点距离后者返回交点处的单位法线。Scene 负责把所有几何体放在一起并按顺序找最近交点。class Geometry: def intersect(self, ray: Ray): raise NotImplementedError def normal(self, point): raise NotImplementedError class Sphere(Geometry): def __init__(self, center, radius): self.center np.asarray(center, dtypefloat) self.radius radius def intersect(self, ray: Ray): oc ray.origin - self.center b 2.0 * np.dot(oc, ray.direction) c np.dot(oc, oc) - self.radius * self.radius disc b * b - 4.0 * c if disc 0: return None t (-b - np.sqrt(disc)) / 2.0 return t if t 1e-8 else None def normal(self, point): n point - self.center return n / np.linalg.norm(n) class Paraboloid(Geometry): def __init__(self, f): self.f float(f) def intersect(self, ray: Ray): x0, y0, z0 ray.origin dx, dy, dz ray.direction a (dx * dx dy * dy) / (4.0 * self.f) b (x0 * dx y0 * dy) / (2.0 * self.f) - dz c (x0 * x0 y0 * y0) / (4.0 * self.f) - z0 if abs(a) 1e-12: if abs(b) 1e-12: return None t -c / b return t if t 1e-8 else None disc b * b - 4.0 * a * c if disc 0: return None t (-b - np.sqrt(disc)) / (2.0 * a) return t if t 1e-8 else None def normal(self, point): x, y, _ point n np.array([x / (2.0 * self.f), y / (2.0 * self.f), -1.0]) return n / np.linalg.norm(n)逻辑说明Sphere 的求交直接复用 2.2 的数学公式Paraboloid 求解时a0对应光线方向与抛物面轴线平行的退化情况需要退化成线性方程处理否则过顶点的中心光线会被漏掉。这个边界我第一版写漏过结果光斑中心总有一个洞。参数说明f是抛物面焦距顶点在原点焦点在 (0,0,f)。法线方向取 ±(x/(2f), y/(2f), -1) 的朝向这里选指向镜面反射一侧的符号。class Scene: def __init__(self, objects): self.objects objects def closest_hit(self, ray: Ray): nearest_t None nearest_obj None for obj in self.objects: t obj.intersect(ray) if t is not None: if nearest_t is None or t nearest_t: nearest_t t nearest_obj obj return nearest_t, nearest_objScene 的代码不需要额外说明逻辑就是遍历所有几何体保留最小的t。这里有一个隐藏约定所有几何体共享同一个坐标系后续如果引入旋转平移最好在求交前先把光线变换到局部坐标不要在求交公式里混着世界坐标算。3.2 核心追迹循环发射光线、求交、反射、记录追迹循环的主体是生成一条光线找到最近交点计算交面法线用反射公式改变方向然后把反射光线与接收面的交点记录下来。这里我选的验证场景是平行光入射抛物面凹面镜反射后汇聚到焦点。程序直接在焦平面 zf 上求交点不继续递归这样最简单。def reflect_ray_at_scene(ray: Ray, scene: Scene): t, obj scene.closest_hit(ray) if t is None: return None point ray.origin t * ray.direction normal obj.normal(point) reflected reflect_dir(ray.direction, normal) return point, reflected def run_focal_plane_spot(f10.0, n_rays20000, aperture8.0): scene Scene([Paraboloid(f)]) rng np.random.default_rng(0) hits [] for _ in range(n_rays): x rng.uniform(-aperture, aperture) y rng.uniform(-aperture, aperture) ray Ray(np.array([x, y, 50.0]), np.array([0.0, 0.0, -1.0])) result reflect_ray_at_scene(ray, scene) if result is None: continue point, reflected result if abs(reflected[2]) 1e-9: continue s (f - point[2]) / reflected[2] if s 0: continue hit point s * reflected hits.append(hit) return np.asarray(hits)逻辑说明光线从高度 z50 处竖直向下发出方向为 (0,0,-1)在抛物面上交点后反射。s表示反射光线从交点继续传播到焦平面 zf 所需的参数距离若s0说明反射光线背向焦平面不该记录。参数说明n_rays是发射样本数aperture是入射光在 x-y 平面上的半径f是抛物面焦距。这套流程可以直接复用到球面镜、平面镜甚至多反射结构只改变 Scene 里的几何体列表。3.3 焦平面光斑可视化与数值验证入口把记录的焦点附近交点投影到 x-y 平面画散点图就是最常见的光斑图工程上叫 spot diagram。程序跑完只需要几行代码就能验证结果import matplotlib.pyplot as plt hits run_focal_plane_spot(f10.0, n_rays20000, aperture8.0) fig, ax plt.subplots(figsize(6, 6)) ax.scatter(hits[:, 0], hits[:, 1], s0.1, alpha0.5) ax.set_xlim(-0.01, 0.01) ax.set_ylim(-0.01, 0.01) ax.set_xlabel(x / mm) ax.set_ylabel(y / mm) ax.set_title(Focal plane spot diagram) plt.show()逻辑说明理想抛物面对平行光精确聚焦光斑应该缩成一个点。这里把坐标范围限制在 ±0.01是为了观察数值误差造成的弥散。如果程序实现正确散点会集中在中心附近如果反射公式符号写反光斑会在毫米尺度散开一眼就能看出来。参数说明s0.1让每个点很小alpha0.5控制透明度光线重叠多时颜色加深便于观察泊松噪声。4. 影响结果正确性的关键参数光线数量、波长与表面散射4.1 光线数量与蒙特卡洛噪声光线追迹本质是蒙特卡洛采样落在某个探测器区域内的光线数是泊松过程相对误差大约正比于 1/√kk 是落在该区域内的样本数。所以不是无条件把 N 拉大就有用还要看每个小格子里分到的样本量。加入粗糙反射后这个效应会直接体现在焦面光斑的 RMS 半径上。def run_rough_spot(f, n_rays, aperture, sigma): scene Scene([Paraboloid(f)]) rng np.random.default_rng(1) hits [] for _ in range(n_rays): x rng.uniform(-aperture, aperture) y rng.uniform(-aperture, aperture) ray Ray(np.array([x, y, 50.0]), np.array([0.0, 0.0, -1.0])) t, obj scene.closest_hit(ray) if t is None: continue p ray.origin t * ray.direction n obj.normal(p) n n sigma * rng.standard_normal(3) n n / np.linalg.norm(n) r reflect_dir(ray.direction, n) if abs(r[2]) 1e-9: continue s (f - p[2]) / r[2] if s 0: continue hit p s * r hits.append(hit) return np.asarray(hits) for n in [500, 2000, 8000, 32000]: h run_rough_spot(10.0, n, 8.0, sigma0.002) r np.sqrt(h[:, 0] ** 2 h[:, 1] ** 2) print(fN{n:6d}, mean{r.mean():.5e}, std{r.std():.5e})现象随着 N 增大均方根半径的均值基本不再漂移而标准差按约 1/√N 收缩。这说明光斑的物理宽度由粗糙度决定统计涨落由样本量决定。调参时如果发现均值随 N 一直变优先检查采样是否偏置而不是继续堆光线数。4.2 折射率、吸收系数与菲涅尔反射率不要写死常数很多入门程序把反射率写成一个固定值比如 0.9这在单波长、单一入射角下勉强能看但换一个入射角结果就失真。反射率在透明介质表面必须用菲涅尔公式计算它同时依赖入射角、两种介质的折射率以及偏振方向。渲染里通常取 s 和 p 偏振的平均值。def fresnel_reflectance(cosi, n1, n2): sin2t (n1 / n2) ** 2 * (1.0 - cosi * cosi) if sin2t 1.0: return 1.0 # 全反射 cost np.sqrt(1.0 - sin2t) rs (n1 * cosi - n2 * cost) / (n1 * cosi n2 * cost) rp (n1 * cost - n2 * cosi) / (n1 * cost n2 * cosi) return 0.5 * (rs * rs rp * rp)逻辑说明cosi是入射光线方向与法线夹角的余弦正值表示光线从法线同侧进入。参数说明n1是入射介质折射率n2是出射介质折射率。全反射直接返回 1.0此时菲涅尔反射率就是 1。实际项目中折射率要按波长查表可见光波段通常用柯西公式拟合不要用一个常数跑宽谱仿真。4.3 表面粗糙度模型把平面法线变成随机微元法线真实表面不是理想镜面。粗糙度会让反射方向在一个锥角内散开工程上最简单的做法是给法线加一个零均值高斯扰动扰动的标准差 σ 就是均方根粗糙度角单位是弧度。这个模型是微面元模型的最简形态适合快速评估光斑弥散和耦合效率。def perturb_normal(n, sigma, rng): m n sigma * rng.standard_normal(3) return m / np.linalg.norm(m)逻辑说明sigma越大法线方向抖动越剧烈反射光锥越宽。参数说明对于光学抛光面σ 通常在 0.0001~0.001 弧度普通机械加工金属面可能到 0.01 弧度以上。要注意的是加扰动后必须重新归一化否则法线长度变化会让菲涅尔公式的余弦计算失真。我经常用下面这张表作为参数标定起点实际值应根据自己的实验数据或加工公差调整参数建议起点主要影响常见误用光线数量 N10^4~10^6统计涨落、运行时间只加数量不改善采样分布递归深度4~8间接反射次数、性能设成上百导致指数级膨胀粗糙度 σ0.001 rad散射波瓣宽度和表面抛光指标不匹配折射率按波长查表反射率、透射光偏折用单一常数跑宽谱5. 光线追迹程序避坑指南5 个常见问题与排查方法5.1 自相交浮点误差导致的斑点与鬼影现象反射光线刚离开表面下一轮求交又在原始位置打到自己画面上出现周期性黑点、亮点或者闪烁。原因交点坐标由浮点运算得来到 1e-16 量级误差反射后的起点严格来说落在表面内侧或外侧。下一轮求交时光线起点和表面的距离小于误差系统t会算出一个极小正数于是光线在原点附近反复反射。解决在求交结果里过滤掉过近的t并在生成新光线时沿反射方向偏移一个微小量。EPSILON 1e-6 offset_point point reflected * EPSILON next_ray Ray(offset_point, reflected)注意EPSILON不能设太大否则对薄透镜、邻近曲面会造成真实偏移一般取当前场景尺寸的 1e-6 到 1e-4 量级。5.2 全反射分支写漏光线“穿模”或亮度突变现象光线从玻璃射向空气时在某些角度看会出现折射光线射进错误方向或者表面某部分突然消失。原因refract_dir里根号内为负表示发生了全反射。很多代码只写了折射公式忘了检查sin2t 1.0于是向上取整、取 NaN 或复数值继续算。解决把全反射当成反射处理。t_dir refract_dir(d, n, n1, n2) if t_dir is None: reflected reflect_dir(d, n) else: # 使用 t_dir 继续追迹 pass这个分支在透镜边缘经常触发因为边缘入射角大。第一版程序我在光线追到透镜边缘时发现画面内圈亮外圈暗查半天就是这里漏了分支。5.3 能量不守恒追迹次数越多结果越亮现象把递归深度从 2 改成 6光斑整体亮度明显升高甚至超过入射光能量。原因每次反射或折射都给能量乘一个系数但这个系数可以是反射率、透射率各个独立配置如果两者的和大于 1能量就被“无中生有”。更隐蔽的问题是反射率用了固定值没有考虑入射角导致能量计算和方向计算用的是两套物理。解决保证每个表面满足 RTA1A 是吸收率透射能量用 1-R-A 推导而不是单独写一个可调的透射率。对无吸收介质表面直接调用fresnel_reflectance算反射率透射率为1.0 - R。5.4 法线方向与坐标系约定不一致现象同一个场景改一个几何体顺序或者旋转一个子部件部分表面反射突变折射方向完全错误。原因反射和折射公式里对法线方向的符号约定是敏感的。reflect_dir不关心法线指向哪一侧只要单位长度即可但refract_dir要求法线必须指向入射介质一侧。如果场景几何体统一把法线指向物体外部光线从外部打向表面时恰好相反必须取反。解决在程序里固定一条规则入射侧的约定写进注释并在 debug 模式里对关键光线的法线、方向做断言。cosi -np.dot(d, n) assert cosi 0, 法线方向约定被破坏cosi 应为正这个坑在涉及多层光学系统时尤其频繁。某次我在做一个三层膜结构程序第二层反射方向始终反着排查到最后发现第一层返回的法线在局部坐标系下没做变换。5.5 性能瓶颈纯 Python 循环追迹 100 万条光线的救法现象N 从 1 万涨到 100 万运行时间从几秒涨到几小时内存被列表对象挤爆。原因每一条光线都在 Python 层循环里逐个求交每次求交产生多个临时 NumPy 数组列表动态增长还会造成大量内存碎片。解决先把 Numpy 的向量化用在发射端把一批光线的方向、起点放进固定数组一次算完与一个几何体的交点进一步可以直接用 Numba 把整个求交函数编译成机器码。from numba import njit njit def intersect_paraboloid_batch(origins, directions, f, t_min1e-8): n origins.shape[0] ts np.full(n, -1.0) for i in range(n): x0, y0, z0 origins[i] dx, dy, dz directions[i] a (dx * dx dy * dy) / (4.0 * f) b (x0 * dx y0 * dy) / (2.0 * f) - dz c (x0 * x0 y0 * y0) / (4.0 * f) - z0 if abs(a) 1e-12: if abs(b) 1e-12: continue t -c / b else: disc b * b - 4.0 * a * c if disc 0: continue t (-b - np.sqrt(disc)) / (2.0 * a) if t t_min: ts[i] t return ts逻辑说明这个函数把所有光线的起点和方向打包成二维数组一次调用循环处理全部样本。参数说明t_min是自相交过滤阈值在批处理版本里表现为明显优势因为它避免了对每个交点重新建 Ray 对象。我的经验是 Numba 版本比纯 Python 快 50~200 倍实际提速取决于循环里临时数组的规模。6. 验证你的光线追迹程序解析解对照、能量守恒与采样优化6.1 用抛物面焦点做解析解对照凡是写过追迹程序的人都会有一个共识程序看着能跑不代表算得对。最好的验证方式不是看渲染图漂不漂亮而是找一个有解析解的场景把数值结果直接和理论值比较。抛物面就是这样的场景平行于轴线的平行光经理想抛物面反射后严格经过焦点。f 10.0 hits run_focal_plane_spot(ff, n_rays10000, aperture8.0) focus np.array([0.0, 0.0, f]) err np.linalg.norm(hits - focus, axis1) print(max err:, err.max()) print(rms err:, np.sqrt(np.mean(err ** 2)))如果程序实现正确rms err应该在 1e-7 量级主要来自双精度浮点累积误差。如果跑出来是毫米级优先检查反射公式里的法线方向和求交时的t选择。6.2 能量守恒曲线与收敛性检查第二个必须检查的是能量守恒。在没有吸收的纯镜面场景里一条光线经过任意次反射总能量应该保持不变。加入菲涅尔反射后能量按反射率逐级衰减但总能量曲线应该单调递减且不会反弹。def trace_energy(ray: Ray, scene: Scene, n1, n2, depth0): if depth 5: return 0.0 t, obj scene.closest_hit(ray) if t is None: return 0.0 p ray.origin t * ray.direction n obj.normal(p) cosi -np.dot(ray.direction, n) reflectance fresnel_reflectance(cosi, n1, n2) reflected reflect_dir(ray.direction, n) next_ray Ray(p reflected * 1e-6, reflected) return reflectance * trace_energy(next_ray, scene, n1, n2, depth 1)跑完打印每一深度下的能量如果第 3 层能量反而比第 2 层高说明反射率、吸收系数或法线方向有一处写错了。这类问题在双面表面模型中很常见因为光线从两个方向进出时菲涅尔公式的cosi符号会变化。6.3 进阶圆盘均匀采样与重要性采样当入射光束口径是圆时很多人直接用uniform(-aperture, aperture)分别采样 x 和 y这会让光线在正方形内均匀分布而不是在圆盘内均匀分布。结果就是中心区域单位面积的光线密度高于边缘光斑出现异常的边缘暗、中间亮。正确做法是在圆盘上按面积均匀采样半径要用平方根变换。r aperture * np.sqrt(rng.random(n_rays)) theta 2.0 * np.pi * rng.random(n_rays) x r * np.cos(theta) y r * np.sin(theta)逻辑说明直接对半径均匀采圆盘会导致中心密度过高因为同样一条半径对应的小面积元宽度随半径增大。用平方根变换后单位面积的光线密度才均匀。参数说明n_rays是总样本数aperture是光束半径theta在 [0, 2π) 均匀采样。这个技巧属于很重要的一类“低成本优化”不增加光线总数只改变采样策略高频噪声就能明显下降。我后来在耦合效率计算里把光源发射也改成了按辐射强度分布采样同样样本数下结果的标准差降了一半。光线追迹程序的调优很多时候不是拼机器性能而是拼你对采样分布的理解。那次改完圆盘采样后再看之前直接用正方形采样的光斑图我才意识到之前不少数据波动其实是采样偏置造成的假象。希望这套从建模、实现到验证的流程能帮你少走我走过的弯路。本文还有配套的精品资源点击获取