CN104034058B - 基于gpu的塔式太阳能热电系统镜场的成像方法 - Google Patents

Abstract

本发明公开了一种基于GPU的塔式太阳能热电系统镜场的成像方法，实施步骤如下：（1）确定时间、定日镜及镜场参数、吸热器参数；（2）生成镜场，计算太阳位置、主入射光线单位向量、各定日镜的主反射光线单位向量、单位法向量与余弦效率、大气透射效率；（3）利用CUDA计算平台在GPU上多线程地进行每面定日镜的成像计算；（4）对所有定日镜的成像计算进行汇总，得到镜场成像。本发明中，大量的计算和判断借助了CUDA计算平台的高速运算能力，实现GPU高性能并行运算，整个程序镜场成像能够在保证结果精确的前提下，花费较少的时间，有利于吸热器上能量分析和镜场调度，具有一定的实际意义。

Description

基于GPU的塔式太阳能热电系统镜场的成像方法

技术领域

本发明涉及塔式太阳能热电系统领域，特别是涉及一种基于GPU的塔式太阳能热电系统镜场的成像方法。

背景技术

塔式太阳能热电系统是利用很多个独立跟踪太阳的定日镜，将太阳光聚焦到一个固定在接收塔顶部的吸热器上，加热流经吸热器内部的传热介质，再使用高温介质的热能带动汽轮机、发电机来发电。它是所有大规模太阳能发电技术中成本最低的一种，有着广泛的应用前景。

塔式太阳能热电系统中的聚光、集热子系统，包括镜场和安装在接收塔上的吸热器等，一直是研究的关注点。其中，吸热器连接着镜场和储热系统两端，是光热转换的关键部件。吸热器上的成像源自各定日镜反射的太阳光，研究吸热器上的成像有着重要的意义，是提高太阳利用率、降低镜场投资成本的基础，更可以为镜场调度提供依据。在成像的研究中，不单单是关注于吸热器上的光斑形状、大小等因素，更重要的是获取并分析其能流密度分布。由于成像与镜场中各定日镜反射的光线能量、定日镜与吸热器的位置关系等因素有关，其中涉及到各定日镜的余弦效率、阴影遮挡效率、大气透射效率、溢出效率，而阴影遮挡又受入射光线、各定日镜之间的相互位置等影响，所以成像是个复杂的问题。特别是在考虑到太阳光是非平行光时，研究成像的难度又有所增加。

现有的吸热器上成像的研究中，有的技术通过在定日镜表面生成若干锥形光，在成像过程中叠加各锥形光，并利用卷积计算方法来实现。该技术计算量大，时间代价高；而且通常不考虑到定日镜间阴影遮挡，与现实情况有所差别。也有的技术在生成入射光线开始，到判断定日镜间的阴影遮挡，最后到光线反射至吸热器的整个过程中，采用光线追迹法来进行成像。该方法虽然可以直观清晰地判断每一根光线的轨迹，但是全程跟踪光线，这也无法做到精度与速度的兼顾。

发明内容

本发明提供了一种基于GPU的塔式太阳能热电系统镜场的成像方法，在保证结果较精确的情况下，使得计算效率更高。

本发明采用的技术方案如下：

(1)确定时间、定日镜及镜场参数、吸热器参数。

输入仿真的时间点、定日镜的尺寸和中心离地高度、吸热器的尺寸和离地高度，选择镜场的排布方式和定日镜间隔，为后续的计算做好数据准备。

(2)生成镜场，计算太阳位置、主入射光线单位向量、各定日镜的主反射光线单位向量、单位法向量与余弦效率、大气透射效率。

根据定日镜及镜场参数，生成特定的镜场，得到各定日镜中心坐标；由仿真的时间点，求取太阳的高度角和方位角，进而计算主入射光线单位向量，再结合各定日镜的中心坐标和吸热器参数，计算出各定日镜的主反射光线单位向量、单位法向量以及余弦效率、大气透射效率。余弦效率等于主入射光线单位向量跟定日镜表面单位法向量之间夹角的余弦值；大气透射效率是指光线被定日镜反射至吸热器的过程中，经传播后的反射光线强度与反射初的光线强度的比值。

(3)利用CUDA计算平台在GPU上多线程地进行每面定日镜的成像计算。

首先需要确定每面定日镜被阴影、遮挡的情况。对于某面被计算定日镜，先判断能对它造成阴影、遮挡的其它定日镜有哪些。这是因为，一个镜场中定日镜的数量中众多，对于某面定日镜，能对它造成阴影、遮挡的其它定日镜数量却不多，所以对该定日镜，如果直接计算其它所有定日镜对它造成的阴影、遮挡情况，会有大量不必要的计算，十分耗时。随后在被计算定日镜上进行随机撒点，根据一定的判断准则判断该定日镜上未被阴影、遮挡的随机点。

本发明中考虑太阳光为非平行光，也就是说，入射光实际为光锥，定日镜的反射光也是光锥。被计算定日镜上未被阴影、遮挡的随机点都接收一个太阳光锥，并且反射一个光锥。每个反射光锥中追迹一定数量的光线，被追迹的每根光线携带一定的能量，相当于是利用各光线能量代表光锥形成的光斑中不同点处的能量，即各光线能量总和代表整个光锥能量。这里可以计算被追迹的每根光线的向量表示和包含的能量。

将吸热器表面网格化，即按照吸热器两组边长进行等长划分。根据各反射光锥中被追迹的每根光线的向量表示、吸热器表面的平面方程及其边界范围限定，可求取每根光线与吸热器表面的交点坐标，并确定交点处于哪个吸热器网格中，其所对应的光线的能量就计算在该网格中心点处。至此可获得吸热器上的能流密度分布，从而完成某面被计算定日镜相关的成像计算。其余定日镜的成像计算也可由此得到。

上述这些过程中，除了一些计算量不大的计算和一些必要的负责串行计算的内容是在CPU中实现的，其余计算基本是在CUDA计算平台上，通过GPU并行实现的。

(4)对所有定日镜的成像计算进行汇总，得到镜场成像。

将GPU中各定日镜的成像计算结果调回到CPU中，进行汇总，可以得到镜场成像后吸热器上的能流密度分布，并画出成像图。

本发明的有益效果是：本发明提出的塔式太阳能热电系统镜场的成像方法中，大量的计算和判断借助了CUDA计算平台的高速运算能力，实现GPU高性能并行运算，整个程序镜场成像能够在保证结果精确的前提下，花费较少的时间，有利于吸热器上能量分析和镜场调度，具有一定的实际意义。而且该方法适用于不同排布方式的镜场和不同类型的吸热器，通用性好。

附图说明

图1是基于GPU的塔式太阳能热电系统镜场的成像方法流程图；

图2是确定阴影判断矩形的俯视图；

图3是确定阴影判断矩形的侧视图；

图4是判断镜面随机点是否被阴影遮挡的示意图；

图5是太阳光锥示意图；

图6是反射光锥形成的圆形光斑示意图；

图7是实施例中的麦田型镜场；

图8是实施例中的吸热器上的成像图。

具体实施方式

如图1所示，一种基于GPU的塔式太阳能热电系统镜场的成像方法，实施步骤如下：

(1)确定时间、定日镜及镜场参数、吸热器参数。

输入仿真的时间点、定日镜的尺寸和中心离地高度、吸热器的尺寸和离地高度，选择镜场的排布方式和定日镜间隔，以便为后续的计算做好数据准备。

(2)生成镜场，计算太阳位置、主入射光线单位向量、各定日镜的主反射光线单位向量、单位法向量与余弦效率、大气透射效率。

首先建立右手直角坐标系，其原点位于接收塔底部中心，x轴、y轴和z轴的正方向分别指向正东、正北和天顶。根据定日镜及镜场参数，生成特定的镜场，镜场排布方式可以是麦田型或辐射型等，也可以是实际电站中的非规则排布。同时得到各定日镜中心坐标。

太阳位置通常用高度角θ_z以及方位角θ_s表示。高度角是指太阳入射光与其在地面投影所成的锐角，范围为0～90°；方位角是指太阳入射光在地面投影光线的走势，θ_s的取值范围定为-180°到180°，0°为正南，顺时针为负。本发明中选取了由美国国家可再生能源实验室研究人员设计的MATLAB程序，输入时间、地理信息等变量直接获得太阳的高度角和方位角。

在获取太阳高度角和方位角之后计算主入射光线单位向量S(x_S,y_S,z_S)：

镜场中某面定日镜的位置用镜面中心坐标P(x_P,y_P,z_P)表示。假设该定日镜中心的主反射光线指向吸热器中心，则计算主反射光线单位向量R(x_R,y_R,z_R)：

$R=\frac{Q-P}{|Q-P|}=\frac{(-x_P,-y_P,H-z_P)}{\sqrt{x_P^2+y_P^2+{(H-z_P)}^2}}---\left(2\right)$

式中，Q(0,0,H)为吸热器中心点坐标，H为吸热器中心高度。

除了定日镜的位置外，还需要对它的朝向进行描述，这里使用定日镜的单位法向量Nr(x_Nr,y_Nr,z_Nr)来表示：

$Nr=\frac{R-S}{|R-S|}---\left(3\right)$

余弦效率η_cos等于主入射光线单位向量跟定日镜表面单位法向量夹角的余弦值：

η_cos＝S·Nr＝x_Sx_Nr+y_Sy_Nr+z_Sz_Nr(4)

在晴朗的天气中，定日镜的大气透射效率η_att可以表示成定日镜中心到吸热器中心距离d_HR的函数：

${\mathrm{\η}}_{att}=\left\{\begin{array}{cc}0.99321-0.0001176\×d_{HR}+1.97\×{10}^{-8}\×d_{HR}^2& (d_{HR}\≤1000m)\\ e^{-0.0001106\×d_{HR}}& (d_{HR}\≥1000m)\end{array}\right.---\left(5\right)$

(3)利用CUDA计算平台在GPU上多线程地进行每面定日镜的成像计算。

由于镜场中各个定日镜之间的相互关系，成像时不得不考虑定日镜间的阴影、遮挡问题。

一个镜场中定日镜的数量比较多，少则几百面，多则上万面。对于某一面定日镜，能对它造成阴影、遮挡的其它定日镜数量其实不多，因此对某一面定日镜，如果直接计算其它所有定日镜对它造成的阴影、遮挡，会有大量不必要的计算，十分耗时。所以，应该先确定能对被计算的定日镜造成阴影、遮挡的有哪些定日镜，这样会提高效率。

本发明采用一种借助阴影判断矩形和遮挡判断矩形分别来确定能对某面被计算定日镜造成阴影和遮挡的其它定日镜的方法。如果其它定日镜中心落在阴影判断矩形和遮挡判断矩形之内，则它会对被计算镜子分别造成阴影和遮挡。这里以阴影判断矩形为例来说明。

以定日镜A作为被计算定日镜。在镜场俯视图上建立一个矩形K，使镜面A的中心位于矩形边L1的中点位置，如图2所示，并使另外一对边L2与太阳入射光在x-y平面的分量平行。图中的等价圆表示定日镜转动时，线段端点形成的轨迹。矩形两组边长的确定方式如下。假设定日镜B的镜面边缘的投影刚好和A边缘相切，此时过镜面B中心做到L1的垂线，交点到A中心点距离为l，当A,B均和阳光在x-y平面的分量垂直，即阳光分别与A,B的等价圆相切时，l取最大值lm，lm为定日镜的长，如果l大于lm则B任何时刻都不会遮挡镜面A，因此令L1＝2lm。

之后确定L2的长度。图3为图2的侧视图，假设通过镜面B上边缘的太阳入射光刚好和镜面A的下边缘相切，为求L₂的最大值我们假设镜面A、B和太阳入射光垂直，即A、B的等价圆刚好和太阳入射光在不同侧相切，此时有(0.5w_m)/(0.5L₂)＝sinθ_z，w_m为定日镜的宽。

因此可以确定阴影判断矩形的两组边长L₁、L₂：

$\left\{\begin{array}{c}{L}_1=2l_m\\ L_2=w_m/{\mathrm{sin\θ}}_z\end{array}\right.---\left(6\right)$

将入射光换成反射光，可以类似地确定遮挡判断矩形的两组边长L₃、L₄：

$\left\{\begin{array}{c}{L}_3=2l_m\\ L_4=w_m/{\mathrm{sin\θ}}_r\end{array}\right.---\left(7\right)$

其中，θ_r为定日镜反射光线的高度角。

为了确定某面定日镜被阴影、遮挡的具体情况，需要在定日镜表面进行随机撒点。从理论上说，撒点数量越多，精度就越高；如果撒点数过少，就会跟实际情况偏差大。接下来介绍随机点是否被阴影、遮挡的判断方法。

图4中的投影四边形是指能对被计算定日镜造成阴影、遮挡的某面定日镜沿着入射光或反射光在被计算定日镜平面上的投影，P_R是被计算定日镜上的某一随机点，其x坐标为，y坐标为。过随机点P_R作平行于x轴的直线和投影四边形两组对边或其延长线相交于四个点，对应x坐标为，如果满足：

$\left\{\begin{array}{c}(x_{p_1}-x_{p_R})(x_{p_3}-x_{p_R})<0\\ (x_{p_2}-x_{p_R})(x_{p_4}-x_{p_R})<0\end{array}\right.---\left(8\right)$

则点P_R在投影四边形内，也就是说P_R被阴影、遮挡了，否则P_R不会被阴影、遮挡。对于过随机点P_R作平行于x轴的直线跟投影四边形的其中一组对边平行的情况，则需过随机点P_R再作与y轴平行的直线，此时，过随机点P_R作平行于x轴的直线与投影四边形的一组对边或其延长线相交于两个点，对应x坐标为 ，过随机点P_R作平行于y轴的直线与投影四边形的另一组对边或其延长线相交于两个点，对应y坐标为，如果满足：

$\left\{\begin{array}{c}(x_{p_5}-x_{p_R})(x_{p_6}-x_{p_R})<0\\ (y_{p_7}-y_{p_R})(y_{p_8}-y_{p_R})<0\end{array}\right.---\left(9\right)$

则点P_R在投影四边形内，也就是说P_R被阴影、遮挡了，否则P_R不会被阴影、遮挡。

本发明中考虑太阳光为非平行光。在进行成像仿真时，吸热器与定日镜之间的距离相对较远，由于光线的发散性，距离越远的定日镜在吸热器上的光斑会偏大，距离近的所形成的光斑影像偏小。因此，若采用平行光会造成较大的误差，非平行光下吸热器上的成像仿真更具实际意义。太阳上各点向任意方向都有光线射出，也就是说太阳光照射到地球上某一点实际上为光锥，如图5所示，太阳光锥的锥角为9.3mrad。

地球上一点观察太阳的像为圆形，称为SolarDisk，该像上能量分布并不均匀，由内而外逐渐减小。其能流密度分布可由如下公式描述：

$f\left(\mathrm{\&alpha;}\right)=\left\{\begin{array}{cc}{S}_0\{1-\mathrm{\&lambda;}{\left(\frac{\mathrm{\&alpha;}}{{\mathrm{\&alpha;}}_s}\right)}^4\},& \mathrm{\&alpha;}\&le;{\mathrm{\&alpha;}}_s\\ 0,& \mathrm{\&alpha;}>{\mathrm{\&alpha;}}_s\end{array}\right.---\left(10\right)$

f(α)表示的是SolarDisk上一点的能流密度，其中λ＝0.5138，α_s是光锥锥角的一半，为4.6mrad，α是SolarDisk上一点到观察点的连线同SolarDisk中心点和观察点连线的夹角，当α大于α_s时，f(α)等于0，即SolarDisk外部是黑体。S₀同f(α)单位一样，也为W/m²，其数值取决于观察点到太阳的距离。

镜面上未被阴影遮挡的每个随机点都接收一个太阳光锥，并且反射一个光锥。每个反射光锥中追迹一定数量的光线，被追迹的每根光线携带一定的能量，相当于是利用各光线能量代表光锥形成的光斑中不同点处的能量，即各光线能量总和代表整个光锥能量。这里需要计算每根光线的能量。

假设一个反射光锥照射在某个垂直于光锥中心线的平面上形成圆形光斑，如图6所示。该光斑的能流密度分布服从式(10)。图中的点表示被追迹的光线跟平面的交点，每个交点又可代表一个面积微元，它的能量就是光线携带的能量，于是计算光线能量的问题就转换成计算各面积微元的能量。

每个面积微元的大小ds：

$ds=\frac{{\mathrm{\&pi;R}}_s^2}{n}---\left(11\right)$

式中，R_s为光斑的半径，n为的一个反射光锥中追迹的光线数。

在该情况下，可以确定式(10)中的S₀：

$S_0=\frac{G_{bn}}{{\mathrm{\&pi;R}}_s^2(1-\frac{\mathrm{\&lambda;}}{3})}---\left(12\right)$

上式中G_bn为反射光锥的能量：

$G_{bn}=\frac{DNI\&times;S_h\&times;{\mathrm{\&eta;}}_{cos}\&times;\mathrm{\&rho;}}{N}---\left(13\right)$

其中，N是镜面撒点数，DNI为太阳光垂直入射强度，S_h为定日镜表面积，ρ为镜面反射率。

由式(10～13)可以得到单个面积微元的能量，也即反射光锥中单根光线的能量dE：

$dE=S_0\{1-\mathrm{\&lambda;}{\left(\frac{\mathrm{\&alpha;}}{{\mathrm{\&alpha;}}_s}\right)}^4\}ds=\frac{DNI\&times;S_h\&times;{\mathrm{\&eta;}}_{\cos }\&times;\mathrm{\&rho;}}{N\&times;n(1-\frac{\mathrm{\&lambda;}}{3})}\{1-\mathrm{\&lambda;}{\left(\frac{\mathrm{\&alpha;}}{{\mathrm{\&alpha;}}_s}\right)}^4\}---\left(14\right)$

接着将吸热器表面网格化，即按照吸热器两组边长进行等长划分。根据反射光锥中被追迹的每根光线的向量表示、吸热器表面的平面方程及其边界范围限定，可求取每根光线与吸热器表面的交点坐标。若交点不在吸热器的边界范围内则舍弃该交点。

对于有效范围内的交点，它处在某个吸热器网格内，其所对应的光线的能量就计算在该网格中心点处。网格中心点的所有能量加和之后，除以网格的面积，就是网格中心点处的能流密度大小。这一过程用公式可表示成：

$f_{i,j}=\frac{m\underset{k=1}{\overset{K}{\&Sigma;}}{\mathrm{dE}}_k}{S_r}---\left(15\right)$

式中，f_i,j为第i行、第j列网格中心点处的能流密度，S_r为吸热器表面面积，m为网格数，表示处在该网格内的所有交点对应的光线能量总和，k为处在该网格内的各交点的编号，K为处在该网格内的交点总数。

至此可获得吸热器上的能流密度分布，从而完成某面被计算定日镜相关的成像计算。其余定日镜的成像计算也可由此得到。

上述这些过程中，除了一些计算量不大的计算和一些必要的负责串行计算的内容是在CPU中实现的，其余计算基本是在CUDA计算平台上，通过GPU并行实现的。

(4)对所有定日镜的成像计算进行汇总，得到镜场成像。

将GPU中各定日镜的成像计算结果调回到CPU中，进行汇总，可以得到全镜场成像后吸热器上的能流密度分布，并画出成像图。

本发明实施例应用于一个包含2500面定日镜的麦田型镜场，如图7所示(实心点代表定日镜，圆圈代表装有吸热器的接收塔的所在位置)。使用的计算机操作系统为Windows7，使用的GPU设备为支持CUDA计算的NVIDIAGeforceGTX570显卡，在MATLAB(R2010b)中引入CUDA计算平台，实现基于GPU的塔式太阳能热电系统镜场的成像的并行计算。每面定日镜上随机撒点1000个，每个随机点反射的光锥中追迹100根光线，即每面定日镜反射了107根光线。实施例中，塔式太阳能热电系统镜场的成像耗时1.5秒，得到吸热器各网格中心点处的能流密度后，画出的成像图如图8所示，该图为能流密度等高线图，能流密度由内而外逐渐减小,平均能流密度为179.8W/m²。

Claims (1)

1.一种基于GPU的塔式太阳能热电系统镜场的成像方法，其特征在于它的步骤如下：

1)输入仿真的时间点、定日镜参数、镜场参数和吸热器参数，定日镜参数包括定日镜的尺寸和中心离地高度，镜场参数包括镜场的排布方式和定日镜间隔，吸热器参数包括吸热器的尺寸、离地高度和位置；

2)根据定日镜参数、镜场参数，生成特定的镜场，得到每面定日镜中心坐标，由仿真的时间点，求取太阳的高度角和方位角，进而计算主入射光线单位向量，再结合每面定日镜的中心坐标和吸热器参数，计算出每面定日镜的主反射光线单位向量、单位法向量以及余弦效率、大气透射效率；

所述的步骤2)具体为：

建立右手直角坐标系，其原点位于接收塔底部中心，x轴、y轴和z轴的正方向分别指向正东、正北和天顶，根据定日镜的尺寸和中心离地高度、镜场的排布方式和定日镜间隔，生成特定的镜场，得到每面定日镜中心坐标；

太阳位置用高度角θ_z以及方位角θ_s表示，高度角是指太阳入射光与其在地面投影所成的锐角，范围为0～90°；方位角是指太阳入射光在地面投影光线的走势，θ_s的取值范围定为-180°到180°，0°为正南，顺时针为负；

在获取太阳高度角和方位角之后计算主入射光线单位向量S(x_S,y_S,z_S)：

镜场中某面定日镜的位置用其镜面中心坐标P(x_P,y_P,z_P)表示，假设该定日镜中心的主反射光线指向吸热器中心，则主反射光线单位向量R(x_R,y_R,z_R)为：

$R=\frac{Q-P}{|Q-P|}=\frac{(-x_P,-y_P,H-z_P)}{\sqrt{x_P^2+y_P^2+(H-z_P}{)}^2}---\left(2\right)$

式中，Q(0,0,H)为吸热器中心点坐标，H为吸热器中心高度；

除了定日镜的位置外，还需要对它的朝向进行描述，使用定日镜的单位法向量Nr(x_Nr,y_Nr,z_Nr)来表示：

$Nr=\frac{R-S}{|R-S|}---\left(3\right)$

余弦效率η_cos等于主入射光线单位向量跟镜面单位法向量夹角的余弦值：

η_cos＝S·Nr＝x_Sx_Nr+y_Sy_Nr+z_Sz_Nr(4)

大气透射效率是指光线被定日镜反射至吸热器的过程中，经传播后的反射光线强度与反射初的光线强度的比值，在晴朗的天气中，大气透射效率η_att表示成定日镜中心到吸热器中心距离d_HR的函数：

${\mathrm{\&eta;}}_{att}=\left\{\begin{array}{cc}0.99321-0.0001176\&times;d_{HR}+1.97\&times;{10}^{-8}\&times;d_{HR}^2& \left(d_{HR}\&le;1000m\right)\\ e^{-0.0001106\&times;d_{HR}}& \left(d_{HR}>1000m\right)\end{array}\right.---\left(5\right);$

3)利用CUDA计算平台在GPU上多线程地进行每面定日镜的成像计算；

所述的步骤(3)为：

在CUDA计算平台上利用GPU实现并行计算,利用阴影判断矩形和遮挡判断矩形分别来确定能对某面被计算定日镜造成阴影和遮挡的其它定日镜有哪些，如果其它定日镜中心落在阴影判断矩形和遮挡判断矩形之内，则它会对被计算镜子分别造成阴影和遮挡，阴影判断矩形的两组边长L1、L2和遮挡判断矩形的两组边长L3、L4分别由式(6)、(7)得到：

$\left\{\begin{array}{c}{L}_1=2l_m\\ L_2=w_m/{\mathrm{sin\&theta;}}_z\end{array}\right.---\left(6\right)$

$\left\{\begin{array}{c}{L}_3=2l_m\\ L_4=w_m/{\mathrm{sin\&theta;}}_r\end{array}\right.---\left(7\right)$

其中，l_m为定日镜的长，w_m为定日镜的宽，θ_z为太阳高度角，θ_r为定日镜反射光线的高度角；

为了确定被计算定日镜被阴影、遮挡的具体情况，需要在定日镜表面进行随机撒点，判断随机点是否被阴影、遮挡的方法如下：把能对被计算定日镜造成阴影和遮挡的其它定日镜分别沿着入射光和反射光在被计算定日镜平面上进行投影，形成投影四边形；P_R是被计算定日镜上的某一随机点，其x坐标为y坐标为过随机点P_R作平行于x轴的直线和投影四边形两组对边或其延长线相交于四个点，对应x坐标为如果满足：

$\left\{\begin{array}{c}(x_{p_1}-x_{p_R})(x_{p_3}-x_{p_R})<0\\ (x_{p_2}-x_{p_R})(x_{p_4}-x_{p_R})<0\end{array}\right.---\left(8\right)$

则点P_R在投影四边形内，也就是说P_R被阴影、遮挡了，否则P_R不会被阴影、遮挡；对于过随机点P_R作平行于x轴的直线跟投影四边形的其中一组对边平行的情况，则需过随机点P_R再作与y轴平行的直线，此时，过随机点P_R作平行于x轴的直线与投影四边形的一组对边或其延长线相交于两个点，对应x坐标为 过随机点P_R作平行于y轴的直线与投影四边形的另一组对边或其延长线相交于两个点，对应y坐标为如果满足：

$\left\{\begin{array}{c}(x_{p_5}-x_{p_R})(x_{p_6}-x_{p_R})<0\\ (y_{p_7}-y_{p_R})(y_{p_8}-y_{p_R})<0\end{array}\right.---\left(9\right)$

则点P_R在投影四边形内，也就是说P_R被阴影、遮挡了，否则P_R不会被阴影、遮挡；

在太阳光为非平行光的条件下，入射光实际为光锥，定日镜的反射光也是光锥，锥角为9.3mrad，地球上一点观察太阳的像为圆形，称为SolarDisk，该像上能量分布并不均匀，由内而外逐渐减小，其能流密度分布由如下公式描述：

$f\left(\mathrm{\&alpha;}\right)=\left\{\begin{array}{cc}{S}_0\{1-\mathrm{\&lambda;}{\left(\frac{\mathrm{\&alpha;}}{{\mathrm{\&alpha;}}_s}\right)}^4\},& \mathrm{\&alpha;}\&le;{\mathrm{\&alpha;}}_s\\ 0,& \mathrm{\&alpha;}>{\mathrm{\&alpha;}}_s\end{array}\right.---\left(10\right)$

其中，f(α)表示的是SolarDisk上一点的能流密度，单位为W/m²，λ＝0.5138，α_s是光锥锥角的一半，为4.6mrad，α是SolarDisk上一点到观察点的连线同SolarDisk中心点和观察点连线的夹角，当α大于α_s时，f(α)等于0，即SolarDisk外部是黑体；S₀同f(α)单位一样，也为W/m²，其数值取决于观察点到太阳的距离；

镜面上未被阴影、遮挡的每个随机点都接收一个太阳光锥，并且反射一个光锥，每个反射光锥中追迹一定数量的光线，被追迹的每根光线携带一定的能量，相当于是利用各光线能量代表光锥形成的光斑中不同点处的能量，即各光线能量总和代表整个光锥能量；

假设一个反射光锥照射在某个垂直于光锥中心线的平面上形成圆形光斑，该光斑的能流密度分布服从式(10)；光锥中被追迹的各光线与平面的每个交点又可各代表一个面积微元，它的能量就是光线携带的能量，计算光线能量的问题就转换成计算各面积微元的能量；

每个面积微元的大小ds：

$ds=\frac{{\mathrm{\&pi;R}}_s^2}{n}---\left(11\right)$

式中，R_s为圆形光斑的半径，n为的一个反射光锥中追迹的光线数；

在该情况下，可以确定式(10)中的S₀：

$S_0=\frac{G_{bn}}{{\mathrm{\&pi;R}}_s^2(1-\frac{\mathrm{\&lambda;}}{3})}---\left(12\right)$ 上式中G_bn为反射光锥的能量：

$G_{bn}=\frac{DNI\&times;S_h\&times;{\mathrm{\&eta;}}_{cos}\&times;\mathrm{\&rho;}}{N}---\left(13\right)$

其中，N是一面镜面上撒的随机点数，DNI为太阳光垂直入射强度，S_h为定日镜表面积，ρ为镜面反射率；

由式(10～13)得到单个面积微元的能量，也即反射光锥中单根光线的能量dE：

$dE=S_0\{1-\mathrm{\&lambda;}{\left(\frac{\mathrm{\&alpha;}}{{\mathrm{\&alpha;}}_s}\right)}^4\}ds=\frac{DNI\&times;S_h\&times;{\mathrm{\&eta;}}_{cos}\&times;\mathrm{\&rho;}}{N\&times;n(1-\frac{\mathrm{\&lambda;}}{3})}\{1-\mathrm{\&lambda;}{\left(\frac{\mathrm{\&alpha;}}{{\mathrm{\&alpha;}}_s}\right)}^4\}---\left(14\right)$

接着将吸热器表面网格化，即按照吸热器两组边长进行等长划分，根据反射光锥中被追迹的每根光线的向量表示、吸热器表面的平面方程及其边界范围限定，求取每根光线与吸热器表面的交点坐标，若交点不在吸热器的边界范围内则舍弃该交点，对于有效范围内的交点，它处在某个吸热器网格内，其所对应的光线的能量就计算在该网格中心点处，网格中心点的所有能量加和之后，除以网格的面积，就是网格中心点处的能流密度大小，这一过程用公式可表示成：

$f_{i,j}=\frac{m\underset{k=1}{\overset{K}{\&Sigma;}}{\mathrm{dE}}_k}{S_r}---\left(15\right)$

式中，f_i,j为第i行、第j列网格中心点处的能流密度，S_r为吸热器表面面积，m为网格数，表示处在该网格内的所有交点对应的光线能量总和，k为处在该网格内的各交点的编号，K为处在该网格内的交点总数；

至此可获得吸热器上的能流密度分布，从而完成某面被计算定日镜相关的成像计算，其余定日镜的成像计算也可由此得到；

4)将GPU中每面定日镜的成像计算结果调回到CPU中，进行汇总，得到镜场成像后吸热器上的能流密度分布，并画出成像图。