非序列光线追迹程序照度分布计算的随机误差分析_郑晓东
第37卷第10期
2008年10月光 子 学 报
Vo l . 37N o . 10O ctober 2008
非序列光线追迹程序照度分布计算的随机误差分析
郑晓东1, 汪扬春1, 秦文红2
(1浙江大学现代光学仪器国家重点实验室, 杭州310027)
(2杭州锐力光学有限公司, 杭州310012)
摘 要:以点光源在平面上产生的照度分布为基准, 研究了非序列光线追迹程序计算照度分布时的随机误差. 结果表明, 随机误差和计算面上网格的划分方式无关, 和各网格内光线条数的平方根成反比. 给出了网格内光线条数与随机误差之间的关系公式. 根据该公式, 为了使结果的标准差优于2%,则单个网格内的光线条数需要大于1000条. 根据具体计算目标, 合理选择网格划分方法及追
迹光线条数才能获得合理的计算结果.
关键词:非序列; 光线追迹; 照度分布; 蒙特卡洛; 随机误差
中图分类号:O439 文献标识码:A 文章编号:1004-4213(2008) 10-1970-5
0 引言
随着航天、天文、精细加工等高科技领域的不断发展, 对光学系统的性能提出了前所未有的高要求. 对应于这些要求, 计算光、机结构整体特性的非序列光线追迹方法越来越广泛地应用于照明系统设计、杂散光控制、照度均匀性计算、鬼像分析和非成像系统的一般设计等领域, 这些都是用序列光线追迹法所无法完成的.
非序列光线追迹的具体算法可参见文献[1]. 但很少有对非序列光线追迹的计算准确度进行研究的文章. 利用非序列追迹进行设计的文章也没有对所使用程序的计算准确度进行评估
[3-6]
[2-3]
产生大量的随机数, 通过统计其中部分随机数的分布特性获得复杂系统的自然特性, 解决大量用解析
法无法求解的复杂问题. MC 法可以仿真更加复杂的系统, 仿真的结果也更接近实际情况, 但M C 法本身给出的是问题的近似解. 近似程度常常取决于如何合理地定义所要处理的问题.
M C 法的基础是概率论和数理统计中的两个重要数学定理:大数定律和中心极限定理. Khintchine 大数定律为:设X 1, X 2, …, X n 为随机变量序列, 若期望μ=E (X j ) (j =1, 2, …) 有限, 则对任意ε>0有
n ※∞
. 本文以应用非
lim P
n
∑X i -
常广泛的LightToo ls 软件为例, 通过理论公式可以计算的例子对软件的随机误差进行评价. 以此来分析和验证非序列光线追迹程序的计算准确度和使用中需要注意的问题, 并对正确使用非序列仿真程序并得到合理的结果提出建议.
即:当n ※∞时, 算术平均值收敛到期望(或统计平均) μ. 要进一步研究收敛的程度, 做出误差估计, 则要用到中心极限定理:设X 1, X 2, …, X n 为独立同
2分布的随机变量序列, 均值为μ, 方差为σ, 则
1 非序列光线追迹的计算原理
目前通用的三个照明设计软件(ASA P ,
Lig htTools 和TracePro ) 均基于蒙特卡洛(M onte Ca rlo , MC ) 算法. 当所要求解的问题是某种事件出现的概率时, 可以通过某种“试验”的方法, 得到这种事件出现的频率, 或者这个随机变量的平均值, 并用它们作为问题的解. 这是MC 方法的基本思想. 它以概率模型为基础, 按照这个模型所描绘的过程, 用模拟试验的结果作为问题的近似解. M C 法首先
1n Z n 渐近于正态分布N μ, .
利用MC 法的光学设计软件按照某种规律向系统中随机投射大量光线, 但不指定光线和系统内每个物体、表面相交的顺序及相互位置关系. 每条光线在它和物体的相交处可以被吸收、反射、折射、衍射或散射. 光线在系统内的任意方向和空间传播, 软件一直追踪每条光线经过各种材料和表面后所携带的辐射通量发生的变化. 辐照度分布的计算是把接收面划分成矩形网格单元, 把各网格单元内所有光线的辐射功率相加除以网格面积. 真实系统产生的辐照度分布可以被看成是连续的, 而非序列光学追
T el :0571-87951681 Email :xiao do ng zheng @zju . edu . cn
收稿日期:2007-01-14
迹法所得到的照度分布是离散的. 理论上, 根据采
[1]
[7]
频率的两倍即可完全再现照度分布的精确值. 但
实际计算中, 目标面上的照度分布通常是未知的, 是
要计算的对象, 所以合理划分网格本身是个复杂的问题. 本文主要研究由M C 算法引起的随机误差, 量化误差的问题将另文讨论.
在很多情况下, MC 法需要大量的时间和计算机资源来完成复杂系统的计算. 为了对多种可能的情况进行比较以得出最佳方案, 必须尝试计算各种方案. 研究如何估计计算结果的准确度, 在保证所需准确度的前提下, 减少计算时间具有非常现实的意义.
本课题所研究的是光线追迹条数与计算结果的随机误差之间的关系, 通过该研究可以使使用者能够合理选择光线追迹的条数以得到所期望的计算准确度, 同时最大限度地减少计算时间, 提高计算效率.
图1 点光源在平面上产生的照度分布Fig . 1 I lluminance dist ribution on a plane produced
by a point so urce
2. 2 照度分布随机误差的实验验证
LightTo ols 模拟光源时, 从光源发出的每条光线所携带的初始光通量是相同的, 在发光强度强的方向, 发射的光线空间密度高, 发光强度低的方向则光线密度低. 在计算照度分布时, LightTo ols 将接收面划分为若干个小网格, 追迹从光源发出经过整个光学系统中所有部件的折射、反射等变化之后入射到网格单元内的所有光线, 将网格内所有光线所携带的光通量相加除以网格单元面积作为该网格的平均照度.
模拟时分两种情况:1) 固定接收面上划分的网格数, 改变光源发出的光线条数, 即追迹光线的条数; 计算最大随机误差随追迹光线条数变化的情况; 2) 光源发出的光线条数一定, 改变接收面上划分网格的数目, 计算最大随机误差和网格划分数目之间的关系. 计算结果的随机误差用所有网格中的最大误差来表示.
由于本文只考虑随机误差, 确定网格划分数量之后, 用M A TLAB 软件的数值积分功能对网格内各微小面元的照度值对面积进行积分然后除以网格面积作为每个网格照度的理论平均值. 即
E av S mesh
(2)
2 照度分布计算准确度的评价
2. 1 验证照度分布计算准确度的方法
为了验证程序的计算准确度, 可以利用实物实
验的方法, 搭建实验装置, 通过实际测得的照度分布来验证计算准确度. 但利用实际装置的实验不可避免地将引入各种测量误差. 本文要验证的是程序计算值和理论值之间的误差. 故选择了一种在理论上非常容易精确计算的情况:在空间放置一个点光源和一个正方形平面接收面, 根据平方反比定律, 点光源在距离其有限距离的平面上产生的照度是可以精确计算的[8]. 用Light To ols 对同样情况进行模拟, 然后和理论计算的结果相比较, 可以在一定程度上验证软件的计算准确度.
以下计算中所使用的参量为:点光源发光强度1000cd , 平面大小100×100mm 2, 光源与平面中心的连线垂直于平面, 光源距为离平面的垂直距离100mm . 根据平方反比定律y ) 的光照度可表示为
E co s θ
r
[8]
, 平面上任一点(x ,
(1)
式中E 为点光源在距离为r 的面元上产生的照度, 单位lx ; I 为光源在照射方向的发光强度, 单位cd ; r
为点光源到计算面元的距离, 单位m ; θ为被照面元法线与入射光线的夹角.
图1为根据式(1) 计算得到的结果. 中心照度为100000lx ; 四个角上光照度最低, 只有54433. 1lx
.
用式(2) 计算得到的平均照度已经是将连续的
照度分布进行量化之后得到的一系列里离散点, 点的数量和网格划分的数量相对应. 显然网格划分的数量越少, 量化所产生的误差就越大. 极端的情况是将整个计算面作为一个网格, 这时的量化误差通常会很大, 计算结果就是整个面上的平均照度. 对所计算的点光源照射平面的情况, 最大的量化误差发生在被照射面的四个角上.
M C 法所计算的是网格内含有随机误差的平均
照度, 随机误差的计算公式可表示为
simu av
×100(3)
E av
式中E av 为利用式(2) 所得到的网格内照度分布的
Error =
理论平均值, E simu 为Lig htToo ls 软件模拟得到的各网格内的照度计算值. 最大随机误差是按式(3) 计算每个网格内照度的随机误差, 得出整个平面上所有网格的随机误差, 然后选择其中绝对值最大的数作为整个计算的最大随机误差.
模拟时先固定网格数, 光线条数从100k 、400k 、700k 到1M , 得到各种情况下的最大误差; 然后固定光线条数, 改变接收平面上的网格数分别为7×7、11×11、15×15、19×19, 得出各种情况下的最大随机误差, 见表1.
表1 各种情况下模拟计算的最大随机误差误差
网%7×72. 821. 501. 091. 08
11×1115×155. 633. 74
2. 362. 02
7. 824. 84
4. 072. 84
19×1911. 446. 56
6. 144. 46
为了获得随机误差与网格划分及光线追迹条数之间更清晰的关系, 本文以总光线条数除以总网格
数得到的网格的平均光线条数为横坐标, 所有网格中的最大计算误差为纵坐标画图, 表1的计算结果表示为图3.
100k
400k 700k 1M
图3 网格中的平均光线数与最大随机误差的关系Fig . 3 T he maximum rando m erro r a s a function of
ave rage r ay -number s in each cell
从图中可以看出, 无论如何划分网格, 随机误差
与网格内的平均光线数有非常直接的对应关系. 图中的环形标记为计算结果, 实线是对结果的最小二乘法拟合, 结果为
0. 47
E max (%) =6. 0/(N /1000)
(4)
从表1的计算结果可以看出, 当光线条数一定时, 最大误差随网格划分数的增加而增加. 其原因是当追迹光线的条数一定时, 随着网格划分数量的增加, 每个网格内平均光线条数不断减少. 图2是表1
结果的图形化表示, 其纵坐标为根据式(3
) 计算各网格的误差后得到的所有网格中误差的最大值, 横坐标为计算面上每个边所分割的数目, 整个平面上所分割的网格数为两个边分割数目的乘积. 因为本文所计算的是正方形, 两个边所分割的数量相同, 总的网格数为单边分割数目的平方. 从图2可以看出, 最大误差和单边分割的数目有近似线性的关系, 随机误差随每个边网格划分数目的增加近似线性增加, 这意味着误差应该与网格中光线条数的平方根成反比. 因为总的网格数为两个边网格分割数的乘积, 而每个网格中光线平均数与总网格数成反比.
式中N 为每个网格中的平均光线数. 从式(4) 可以
看出, 计算结果的最大随机误差近似与网格内光线数平均值的平方根成反比.
本文希望通过该研究得到一个经验性的光线条数与随机误差之间的关系, 但由于M C 计算方法本身的随机性, 最大误差也是随机分布的.
图3中数据的分布有比较强的离散性, 对估算准确度有一定影响. 其原因如图1, 计算对象的中心和边缘处照度不同, 不同位置网格内光线条数也不相同. 中心处照度高, 光线条数多; 边缘处照度低, 光线条数少. 而本文所使用的是整个平面内网格的平均光线数. 为了消除该误差, 选择了计算面上离中心的距离一定, 理论平均照度相同的八个对称点, 重复进行多次计算. 这时理论上每个网格内光线的条数均应相同. 可以确定, 这样获得的计算结果的误差完全是由M C 方法的统计性所引起. 将计算面划分为13×13个网格, 在计算面上选取距离中心对称的8个网格分别计算了总光线条数为100, 200…1000k 的10种情况.
和最大误差相比较, 标准差与光线条数的关系更具有普遍意义. 为了获得随机误差的标准差, 对每种光线条数, 重复计算21次, 获得21×8=168个数
图2 最大随机误差和每个边的网格数的关系
Fig . 2 T he re
latio nship be tween maximum r andom er ro r
number in o ne side of square
据. 然后对这168个结果进行正态分布拟合得出每种情况的标准差. 根据十种、每种168个数据所得到
以减少随机误差, 但这以牺牲细节分辨能力, 即增加系统误差为代价的. 因为程序所计算的是网格内照
度的平均值, 如果网格过大, 则在网格内照度的变化被平均化, 无法反映实际分布, 所以如何合理选择网格是个非常值得研究的问题.
同样研究了存在折、反射表面的情况, 结果表明这些表面的存在完全不影响随机误差. 以上研究结果广泛适用于各种折、反射光学系统. 但本文研究中
图4 网格内光线数平方根的倒数与标准差的关系
F ig . 4 T he standard err or v s the reciprocal of square r oo t
of the ray numbe rs in a ce ll
没有考虑任何光散射模型本身的误差. 式(4) 的拟合公式是根据整个计算平面上的平均光线条数得出的, 而最大误差基本出现在照度最低的4个角上. 4个角上实际光线的条数小于整个计算面上的平均光线条数, 故式(4) 过多估算了获得一定随机误差所需要的光线条数. 式(5) 更精确地反映了随机误差与网格内光线条数的关系, 是推荐使用的结果. 因为随机误差符合正态分布, 故标准差为2%意味着结果的随机误差小于2%的概率为68. 3%,小于6%的概率为99. 7%.
参考文献
[1] AN Lian -sheng , LI Guo -dong . Compu ter simulation and
analy sis for the inten sity dis trib ution of illumination sy stems [J ]. Opt Tech , 1998, 24(6) :45-47, 73.
安连生, 李国栋. 照明光学系统照度分布的计算机模拟分析[J ]. 光学技术, 1998, 24(6) :45-47, 73.
[2] CASS ARLY W J , H AYFORD M J . Illumnation optimization :
the revolution h as begun [C ]. S P IE , 2002, 4832:258. [3] SH EN M o , LI Hai -feng , LU W ei , et al . T he method of
reflective fly eye lens design for LED illuminating projection sys tem [J ]. Acta Photonica S inica , 2006, 35(1) :93-95. 沈默, 李海峰, 陆巍, 等. 用于LE D 照明的反射型复眼设计方法[J ]. 光子学报, 2006, 35(1) :93-95.
[4] ZHENG Zhen -ron g . Th e analysis of relative illumination for
projection len s based on etendu e [J ]. Acta Photon ica S inica , 2005, 34(1) :55-58.
郑臻荣. 基于etendue 量的液晶投影物镜相对照度分析[J ]. 光子学报, 2005, 34(1) :55-58.
[5] ZHANG Zen g -bao , SH ENG Yi -qiang , W ENG Zhi -cheng , et a l .
Off -axis illumination sys tem design for projection display [J ]. I n fr ared and Laser E ng ineer in g , 2004, 33(3) :316-319. 张增宝, 盛益强, 翁志成, 等. 用于投影显示的离轴照明系统设计[J ]. 外与激光工程, 2004, 13(3) :316-319.
[6] JIN Xin -hu a , ZH ENG Xiao -dong , W ANG Zheng . S tudy on
illu mination u nifo rmity of the end su rface ′s radiation of polygon rod in tegratin g lens [J ]. Acta P hotonica S inica , 2006, 35(12) :1964-1968.
金新华, 郑晓东, 王政. 多边形棒状透镜出端辐射均匀性的研究[J ]. 光子学报, 2006, 35(12) :1964-1968.
[7] OPPENH EIM A V , WI LLS KY A S . Signals &sys tems [M ].
New Jersey :Prentice H all , 1997:Chapter 7.
[8] COATON J R , M A RSDEN M M . Lamps and lighting [M ]. 4th
ed . New York :A rchitectural Press , 2004:Chapter 1.
准差, 横坐标是网格中的光线条数平方根的倒数. 可
以看出, 计算结果的标准差与网格内光线条数的平方根成反比关系. 计算结果的拟合公式为
0. 247E rms (%) =71. 3/-(5)
式中N 为网格内的光线条数. 该拟合公式与原始数据的相关度高于99%.根据式(5) , 当要求某网格处的统计准确度为标准差小于2%时, 网格内的光线数应不少于1000条.
3 结论
利用LightTo ols 等软件进行设计的目的是实际应用. 设计人员付出大量的时间和精力使其仿真结果最大限度地满足设计目标和要求. 设计人员相信软件的计算结果是正确的, 一旦实际制作和测试的结果与计算结果不相符合, 人们常常会认为仅仅靠软件的仿真是靠不住的. 实际上软件的结果是否值得信赖, 在很大程度上取决于使用软件的设计师对软件本身的理解和熟悉的程度.
本文根据非序列光线追击程序的计算原理分析了计算照度分布时所获得的计算结果的随机误差. 对于光刻等对照明均匀性要求较高的领域一般要求照明的不均匀性小于±5%.为了获得真实可信的设计结果, 必须对程序的计算准确度有比较深入的了解, 才能获得真正有意义的结果. 根据以上研究可知, 无论计算时网格如何划分, 随机误差均与网格内实际光线条数的平方根成反比, 当光线条数足够多时, 本文所得出的规律是计算结果中随机误差的标247, N 为网格内准差为E rms (%) =71. 3/N -0. 光线的条数. 即为了获得小于2%的标准差, 网格内光线条数要大于N =[71. 3/(2+0. 247) ]≈1007
条. 由于网格内光线条数和该网格处的照度值成正比, 计算时要特别注意目标面上照度较低的部位, 该处通常随机误差较大. 由于随机误差和网格内的光, , 2
Random Error Analysis of Illumination Distribution Calculated by
Non -sequential Ray Tracing Programs
ZH ENG Xiao -dong , WANG Yang -chun , QIN Wen -hong
(2Realopti x Inc , Hangz hou 310012, China )
Received da te :2007-01-14
1
1
2
(1State K ey L aboratory o f Modern Optical I nstrumentation , Z hejiang University , Hangzhou 310027, China )
A bstract :The com putation random erro r of illuminance distribution w as studied in case o f a plane
illuminated by a point source . The data show that the random error is inversely propor tional to the square roo t of ray hits in one pix el . Every single pixel area has to be hit by 1000ray s at least to achieve a result with standard deviation better than 2%.It can be co ncluded that a prope r selection o f grid density and ray number are impo rtant fo r necessary accuracy in different cases .
Key words :Non -sequential ; Ray tracing ; Illumination distribution ; Mo nte Carlo ; Random
error
ZHENG Xiao -Dong was bo rn in 1962. H e received his Ph . D . degree in electrical and electronics engineering fro m Niigata U niversity , Japan 1997. Now he is w o rking at
Depar tment of Optical Engineering as an associate professo r at Zhejiang University . His majo r research inte rests focus o n no n -im aging o ptics , solid illumination and hig h -speed o pto -mo dule packaging .