技术领域
[0001] 本
发明涉及
太阳能发电模拟技术领域，特别是涉及一种峰值稳定的蒙特卡洛太阳辐射能密度仿真方法。
背景技术
[0002] 随着近年来化石
能源的资源短缺和造成的环境问题日益严重，清洁
可再生能源引起了学术界和工业圈的广泛关注。太阳能就是一种常见的可再生资源中。作为地球的主要来源提供者的太阳，太阳能有着天然的安全、清洁、易获取、可再生且储量丰富特征。若存在一种技术能充分地将太阳能转换为人类所需的其他
能量，那么理论上太阳能可以满足全世界对能源的需求。尽管太阳能有着巨大的利用价值，现阶段太阳能对全球能源供应的贡献仍然微不足道。
[0003] 太阳能光热发电技术是目前最常见的太阳能应用方式之一。而现代太阳能光热发电系统中
定日镜镜场反射投影效率极大程度的决定了太阳能的利用效率。因此为了提升太阳能光热发电系统的发电效率，需要对定日镜镜场的设计布局和聚焦策略进行优化。其中对接收器表面的辐射能密度分布进行仿真是解决这些问题的
基础。
[0004] 目前主要采用两种辐射能密度仿真方法：
光线跟踪和解析方法。
光线跟踪方法对大量光线的路径进行模拟，仿照光路在定日镜镜场中与其他物体的相交、折射、反射等关系，从而模拟推导光线与接收器的交点和能量值。在
计算机图形学应用中，
光线跟踪算法可以绘制出真实感极强的图像，因此在定日镜投影反射应用中也可以获得与真实分布接近的结果。但一方面，由于
光线跟踪算法需要对大量光线路径进行计算，因此计算效率较低；另一方面，由于光线跟踪算法的离散性，其辐射能密度仿真的结果存在不确定性噪音。
[0005] 解析方法的本质是通过一个可以解析计算的表达式，直接计算接收器表面任何一点的辐射能密度值。由于解析方法计算效率较高，在大型镜场的仿真中有很大优势。但由于任何一个解析模型都是对真实光斑的一个逼近。因此，解析方法不但
精度理论低于光线跟踪方法，而且由于为了保证模拟的精度，解析方法需要提前确定针对每个定日镜的解析函数的参数，而这些参数值通常由光线跟踪的结果得出。所以，高效且精确的光线跟踪算法是解析方法相关研究的基础，在辐射能密度仿真领域研究中有着不可替代的地位。
[0006] 目前光线跟踪算法仿真结果的不确定性噪音有两方面：
[0007] 一、即使针对单个定日镜的仿真，光线跟踪反射投影在接收器表面的辐射
能量密度峰值(简称峰值，指接收器表面在单位时间内、单位面积上接收的最大辐射能能量)总是不稳定的。峰值的准确性是评估辐射能密度仿真结果的重要因素。一方面，峰值的错误估计会影响接收器的使用寿命；另一方面，峰值也是许多解析计算模型使用的拟合参数。而由于光线跟踪是一种离散的计算方法，仿真结果的峰值十分不稳定，有些情况下可以达到10％的误差。不稳定的峰值必然会对镜场设计结果造成很大影响。
[0008] 二、光线跟踪反射投影在接收器表面的辐射能量分布与真值分布(ground truth，简称GT)的均方根误差(Root Mean Square Error，简称RMSE)较大。这一点由光线跟踪的每次运行生成不同的随机数造成的。RMSE较大同样会降低仿真结果的可信度，从而对诸如镜场布局设计、场地参数优化、聚焦优化等研究中提供的相关参数与真实情况偏差。
发明内容
[0009] 为发明的目的在于提供一种峰值稳定的蒙特卡洛太阳辐射能密度仿真方法，降低光线跟踪算法仿真结果中RMSE，使得峰值与真值更接近且更稳定。
[0010] 一种峰值稳定的蒙特卡洛太阳辐射能密度仿真方法，包括以下步骤：
[0011] (1)根据太阳入射光模型，对入射光方向进行
采样；
[0012] (2)对入射光与定日镜表面的交点进行采样；
[0013] (3)对入射光进行遮挡测试以获得接收器表面的能量密度分布初始仿真结果；
[0014] (4)将接收器表面的能量密度分布初始仿真结果进行滤波计算，获得蒙特卡洛太阳辐射能密度仿真结果。
[0015] 作为优选，所述太阳入射光模型为Buie太阳光分布模型、PillBox太阳光分布模型或高斯太阳光分布模型。
[0016] 作为优选，步骤(2)中采样使用均匀采样的随机数生成器。使用随机数生成器的进行均匀采样的方法实现简单，计算效率高。
[0017] 作为优选，步骤(3)中遮挡测试方法包括以下步骤：
[0018] (a)如果入射光被场景中其他物体遮挡，停止计算；
[0019] (b)否则，根据入射光方向采样的微表面法向方向，计算其反射光方向；
[0020] (c)对反射光进行遮挡测试；
[0021] (d)如果反射光被场景中其他物体遮挡，停止计算；
[0022] (e)否则，计算反射光与接收器相交的
像素位置，将该入射光能量密度累加到该像素。
[0023] 本
申请是一种仿真方法，所以此处计算反射光与接收器相交的像素位置中像素位置是使用数组进行存储的虚拟位置，将入射光能量密度累加到该像素就是对相应的数组进行处理。
[0024] 作为优选，所述滤波计算为高斯滤波计算、均值滤波计算、截尾均值滤波计算、中位数滤波计算或双边滤波计算。更优选的。所述滤波计算为截尾均值滤波计算。
[0025] 截尾均值滤波计算的方法为：对于接收器表面的能量密度分布初始仿真结果中的每个元素，截尾均值滤波计算时将毗邻的k×k个元素中根据去除比例p删除相应数量的最大值和最小值，再计算剩余元素的平均值作为当前元素处理后的值。其中参数的选择对结果具有很大影响，优选k值为9～13。优选P值为0.01～0.05。
[0026] 本发明峰值稳定的蒙特卡洛太阳辐射能密度仿真方法通过对光线跟踪算法仿真计算的结果加上一步后处理：滤波计算，在保证接收器表面接受的总能量不变的前提下，使得峰值与真值更接近且更稳定，同时降低了RMSE。稳定的峰值和与GT更接近的分布为后续的相关算法研究和仿真提供了可靠保障。
附图说明
[0027] 图1为验证提出的模型准确性设计单面平面定日镜与平面接收器的位置分布的示意图。
[0028] 图2为对每个像素进行截尾均值滤波计算的
流程图。
[0029] 图3为本申请为验证提出方法有效性，对单面平面定日镜投影在平面接收器表面生成的辐射能密度分布进行处理前后的对比图，其中，图3(a)展示的是经典蒙特卡洛算法和GT的辐射能密度分布仿真轮廓图对比图；图3(b)展示的是图3(a)中经典蒙特卡洛算法结果加上本申请处理后，与GT的辐射能密度分布仿真轮廓图对比图。
[0030] 图4为图3(a)和(b)在z坐标值等于0平面的截面图结果对比图。
[0031] 图5为100次光线跟踪结果在进行截尾均值滤波操作前后的统计结果图，其中图5(a)为峰值统计图；图5(b)为RMSE统计图。
具体实施方式
[0033] 以现有的双向蒙特卡洛光线跟踪算法对接收器表面辐射能密度仿真算法为例，其计算步骤如下：
[0034] 1.根据太阳入射模型，对入射光方向进行采样。光线方向由高度
角(设为θ)和方位角(设为 )确定。例如，使用Buie太阳光分布模型时，服从[0,2π]区间内的均匀分布，而θ服从Buie分布模型，分别如公式1和公式2所示：
[0035]
[0036]
[0037] 其中κ和γ的计算方法见公式3，
[0038]
[0039] 其中χ代表控制Buie太阳光分布模型在光线大于太阳张角(4.65mrad)时分布情况的CSR参数。
[0040] 2.对入射光与定日镜表面的交点进行采样。定日镜表面随机生成一个点，其在空间中的位置作为入射光线与定日镜的交点；
[0041] 3.对入射光进行遮挡测试：
[0042] 3.1.如果入射光线被场景中其他物体遮挡，停止计算。
[0043] 3.2.否则，计算其反射光方向。微表面法向扰动方向也由参数高度角(设为α)和方位角(设为β)确定，α和β的分布见公式4和公式5：
[0044] α～N(0,σ2) 公式4
[0045] β′～U(0,2π),β＝arctan(β′) 公式5
[0046] 其中σ是标准差，它的值取决于定日镜表面的平坦程度。
[0047] 3.3.对反射光线进行遮挡测试：
[0048] 3.3.1.如果反射光线被被场景中其他物体遮挡，停止计算。
[0049] 3.3.2.否则，计算反射光线与接收器相交的像素位置，将其能量密度累加到该像素，能量密度计算公式见公式6：
[0050]
[0051] 其中ID是太阳光能量密度，单位是W/m2；Shsub和Srsub分别是定日镜和像素的面积，单位是m2；φ是入射光方向与微表面法向之间的夹角；ρ是反射率，在[0,1]区间内；Nc是仿真中与定日镜相交的光线数目；ηaα是大气衰减系数，其值由大气衰减模型确定，如公式7：
[0052]
[0053] 其中d是发射点到接收器表面的距离，单位是m。
[0054] 4.重复以上1-3步骤，直至生成的光线数目达到设置的Nc。即对每一条入射光均进行上述处理，直到所有的入射光均处理完成。
[0055] 由于原有的蒙特卡洛光线跟踪算法的计算是离散的，且过程中涉及到随机数采样，因此每次计算的峰值结果会有较大浮动，且RMSE值较大。其中RMSE的计算方法见公式8：
[0056]
[0057] 本申请需要解决的技术问题是如何对光线跟踪辐射能密度仿真结果进行后处理，并在保证接收器表面总能量不变的前提下，使得峰值与GT更接近且更稳定。同时，处理后的结果与GT分布之间的RMSE减小。图1为本申请为验证提出的模型准确性设计单面平面定日镜与平面接收器的位置分布的示意图，其具体实验参数如表1所示。其中方位角以北为0°，顺
时针为正。
[0058] 表1单面定日镜反射投影的蒙特卡洛光线跟踪仿真结果实验参数
[0059]
[0060] 有五种广泛使用的滤波方法：高斯滤波，均值滤波，截尾均值滤波，中位数滤波和双边滤波。表2所示的是以上五种滤波操作处理蒙特卡洛仿真结果的RMSE和峰值统计值。其中，未进行滤波操作是指以现有的双向蒙特卡洛光线跟踪算法对接收器表面辐射能密度仿真算法得到的结果。统计结果表明截尾均值滤波优于其他四种方法。因此，本申请中着重对截尾均值滤波进行介绍。
[0061] 表2五种滤波操作处理MCRT仿真结果的RMSE和峰值统计
[0062]
[0063] 对于蒙特卡洛光线跟踪仿真得到的辐射能密度分布，其每个像素都进行截尾均值滤波计算。对于每个像素点，截尾均值滤波将毗邻的k×k个像素中最大和最小的p元素删除(p∈[0,0.5))，其余值进行平均计算，计算的结果作为当前像素新的值。
[0064] 截尾均值滤波中使用两个相同大小的最大堆和最小堆来分别保存最小元素和最大元素集合。最小堆和最大堆为完全二叉树，其中最小堆中每个父
节点的值小于或等于其两个孩
子节点；最大堆与之相反，每个
父节点的值大于或等于其两个孩子节点。值得注意的是，最大堆中保存的是最小的元素集合，而最小堆则是保存最大的元素集合。对光线跟踪计算出的辐射能密度分布中每个像素像素的截尾均值计算步骤如下：
[0065] 1.设置sum的值为0。
[0066] 2.初始化两个大小为 的最小堆和最大堆；
[0067] 3.循环读取该像素毗邻k×k个像素，读取的值设为currentValue：
[0068] 3.1.如果最大堆和最小堆未满，将currentValue插入最小堆和最大堆中；
[0069] 3.2.否则：
[0070] 3.2.1.如果currentValue大于最小堆堆顶元素，删除最小堆堆顶元素，将其累加到sum，并将currentValue插入堆中；
[0071] 3.2.2.如果currentValue小于最大堆堆顶元素，删除最大堆堆顶元素，将其累加到sum，并将currentValue插入堆中；
[0072] 3.2.3.如同时不满足步骤3.2.1和3.2.2，则将currentValue累加到sum。
[0073] 4.前像素新的值等于sum除以保留元素的数目。
[0074] 上述过程的流程图见图2，对光线跟踪计算出的辐射能密度分布中每个像素都使用该截尾均值方法进行计算。特别的，该算法十分适合进行并行计算。在GPU上实现该算法时，可以将算法中的堆存储在共享内存，从而
加速计算效率。GPU中的时间复杂度仅仅为O(k2·log(k·p))。
[0075] 表3罗列了截尾均值滤波在不同k和p取值时的RMSE和峰值统计。可以发现在p在区间0.01～0.05和k值在区间9～13中效果较好，RMSE的值较低，且峰值与GT较为接近。图3(a)是对应的蒙特卡洛光线跟踪仿真结果与GT分布的轮廓对比图，其中GT峰值为380.74W/m2。2
在未加入截尾均值滤波计算时RMSE的值为10.80W/m。图3(b)则是在图3(a)计算的结果上加入截尾均值滤波计算后获得的新分布，其中k和p分别为11和2％，此时RMSE的值为0.83W/m2。图4则为图3在z坐标值等于0平面的截面图结果。不难发现，在未经过截尾均值滤波操作处理的结果充满噪音且边界毛糙，这些噪音就是导致峰值不稳定和RMSE值较大的原因。处理之后，仿真结果分布更加接近GT的分布。
[0076] 表3截尾均值滤波在不同k和p取值时的RMSE和峰值统计
[0077]
[0078] 图5(a)和图5(b)所示的分别是100次光线跟踪结果的峰值和RMSE的值在后处理前后的统计图。可以发现后处理后(增加截尾均值滤波计算，k和p分别为11和2％)的峰值更加稳定，RMSE的值也明显降低。
[0079] 在截尾均值滤波中，由于删除元素的不确定性，总能量会发生变化，但能量的变化在0.5％以内。因此，在严格要求总能量不变的应用场景中，可以将p设为0，将截尾滤波退
化成均值滤波的方式进行计算，从而保持总能量不变。
