一种基于地面激光点云的太阳能潜力评估方法

技术领域

本发明涉及计算机图形学领域，具体涉及一种利用激光点云进行太阳能资源评估的方法。

背景技术

传统的太阳能潜力评估主要是在设备安装过程中采用人为主观估计的方法，然而太阳辐射分布同时取决于时间、气候以及调查环境内物体的空间位置关系，简单的考察和推算难以应对复杂多变的调查环境，如城区、林地、交通路段。对于传统手工测量的方法，其分析效率较低，数据采集和后续处理环节脱离，难以满足需求，因此，迫切需要一种具备时空变换反演能力的精确分析方法。

近年来，LiDAR(Light Detection and Ranging)技术飞速发展，其因高度的场景再现能力而适合于日照遮挡分析。激光扫描仪可以快速，高效，精确地获取调查区域内的三维信息，得到的数据为三维激光点云数据。目前，利用点云数据进行太阳能潜力评估的工作集中于机载点云，且大多处于研究阶段。其处理流程主要是通过原始点云数据生成数字表面模型(DSM)，然后直接在DSM网格上进行阴影计算，这种方法主要用于城市级尺度的太阳能潜力评估。而在微观调查尺度下(例如对于单栋建筑屋顶)，仍需要对数据进行多边形建模来确保分析质量，建筑物目标的还原程度取决于结构细节在建模过程中的保留程度。多边形模型需要大量的手工操作或是先验的几何模型和参数，城市环境的复杂性为该方法带来很大的工作量和挑战，种类繁多的城市部件(树木，电力线，桥梁和围墙等)都需要被置入于模型库中。可以看出，对于点云数据而言，其在太阳能潜力评估中的应用场景还较为单一，主要集中于城区。对于郊区和林区的太阳能资源调查研究工作还较少，植被建模的复杂性也同时制约着该方面研究的发展。

对于机载点云建模类方法。首先，机载点云由于扫描视角的原因难以获取建筑物侧面和凹面结构，导致点云形成层级结构，各层之间存在大量空隙，不利于进行精确的日照分析。机载点云分辨率一般在1米左右，且数据采集及规划需要花费大量时间。其次，已存在的建模方法主要基于表面模型进行阴影计算，容易将表面建模过程中的误差传递到遮挡分析之中。而对建筑物目标的特定建模则更多的依赖于交互式操作，不仅需要先验的几何参数或充足的样本模型库，且对植被，交通设施，山体岩石等难以建模目标的选择性忽略同样容易导致反演时产生误差。

对于直接面向点云的方法，其遮挡分析过程过于简单粗放。直接面向点云的方法都主要依据高程值作为判断依据，忽视物体本身的结构和形状。例如对于植被而言，树冠在高空具有向外延展性，会导致该类方法产生的阴影面积过大，难以适应复杂的场景。

发明内容

本发明的目的在于克服上述现有技术的不足，提供一种基于地面激光点云的太阳能潜力评估方法。

本发明具体采用以下技术方案：

一种基于地面激光点云的太阳能潜力评估方法，包括如下步骤：

S1、对原始点云进行抽稀；

S2、对抽稀后的点云P进行感兴趣区域内的地面点集R的提取，

S3、设置光源偏差控制角；

S4、采用最远点贪心策略(farthest-first traversal)来快速计算基点的位置和个数；

S5、计算基点太阳位置；

S6、采用广义隐藏点移除算法即GHPR(Generalized Hidden Point Removal)算法进行三维点云场景的遮挡分析，从而进行日照模拟计算；

S7、对遮挡分析结果进行二值化阴影绘制；

S8、对点云场景进行太阳辐射计算。

进一步地，S1的具体步骤如下：

输入三维地面激光点云P_raw，首先设置一个抽稀距离L_s，然后采用体素抽稀方法。体素抽稀方法先将点云划分为边长为L_s的正立方体体素，再按体素内点集的中心点作为新的数据点，这些新的点集形成了抽稀后的点云P。

进一步地，S2的具体步骤如下：

S21、首先，采用基于体素高程生长滤波的方法来提取地面点，基于整个点云的最低高程值设置一个阈值高度t_z，该阈值高度通常为感兴趣区域高程值之上3米高度。如果一个体素的中心高程与最低点的相对高度高于该阈值t_z，则认为其不可能是地面点，此外，若一个低于阈值的体素可以生长到该阈值高度t_z，则认为是墙体或垂直物体的一部分，也为非地面点，从而提取所有地面点集；

S22、接着，对地面点集进行欧式距离聚类，对获得的聚类簇标记序号；

S23、最后，计算每一个类簇的平均密度值，选择密度最高的一个作为感兴趣区域(调查对象)，返回抽稀后的点云，根据各点索引值标记感兴趣区域的数据点r_i，r_i∈R。

进一步地，S3的具体步骤如下：

设置模拟光线偏差平行光的最大允许偏差角θ，同时设置太阳距离d，从而可以确定两个基点之间的最大距离l，有：

l＝d/tan(π/2-θ)

进一步地，S4的具体步骤如下:

S41、首先，最远点策略随机选择一个感兴趣区域点作为初始基点o₁，构建基点序列O＝{o₁}；

S42、接着，计算感兴趣区域点集上与该序列最远的点，即并将其添加到O中，为O＝{o₁,o₂}，直到任意一点r_i，r_i∈R，都满足

S43、获得基点序列后O＝{o_k|k∈[1,n]}，对于感兴趣区域点集中的每个点r_i，求取其在基点集上的最近邻点，同属于一个基点的感兴趣区域点构成一个偏差控制区域G_k。

进一步地，S5的具体步骤如下：

S51、对于每个基点o_k，通过一套太阳位置求解模型，来计算其相应的太阳位置s_k，将每个基点从WGS84坐标额外转为一份经纬度坐标，有[o_Ψ,o_Φ,o_h]＝T([o_x,o_y,o_z])，其中，T为坐标变换，[o_Ψ,o_Φ]为基点的经纬度值，o_h为海拔高度值，然后，根据太阳位置模型获得基点o的太阳位置s，s＝[o_Ψ+f_d(y),o_Φ-f_d(x),o_h+z]^T，f_d为将长度化为度数的函数，在经纬度坐标下可以如下计算:

x＝dsin(90°-α)cos(ψ)

y＝dsin(90°-α)sin(ψ)

z＝dcos(90°-α)

S52、将经纬度和高度坐标转回WGS84坐标，获得太阳在点云场景三维空间中的相对位置，s＝T′(s)，T′为T的逆变换。

进一步地，S6的具体步骤如下：

S61、首先，遍历太阳点集S，对于每个太阳位置s_i，以及调查场景抽稀后的三维点集P，P＝{p_i|i∈[1,m]}，其中m是点集的势，执行径向几何变换F，有:

其中，若p_i＝s，则变换后的点集为而f_k是变换的核函数。

S62、然后对点集做一次凸包计算，并保留凸包内的点为阴影点集W。

进一步地，S7的具体步骤如下：

S71、对于每个太阳位置s_k，计算一次W_k；

S72、然后，取W_k与负责区域G_k的交集为各相应区域内的阴影点。遍历完后，将每个区域内的的阴影结果取并集，从而得到整个感兴趣区域的光照结果；

S73、最后将光照结果二值化，设f_i(r_i)为二值化函数用于表示该点r_i的遮挡状况，对感兴趣区域内的阴影点标记为0，记为f_i(r_i)＝0，而非阴影点标记为1，记为f_i(r_i)＝1。

进一步地，S8的具体步骤如下：

S81、采用Hottel模型进行遮挡分析和阴影绘制的结果分析，首先输入太阳时，计算该日大气层外的太阳辐照度I_o，然后输入经纬度计算该地区直射辐射透过率τ_b，再计算斜面修正因子ε，对于修正因子，先求取太阳入射光学与斜面法向量之间的夹角，后采用构建和分解协方差矩阵获得每个感兴趣区域点r_i的法向量接着，连接太阳和当前遍历点构建入射向量为s-r_i，校正因子ε可以通过入射向量与法向量之间的夹角余旋计算:

其中·为向量点乘，而s为该点相关基点所对应的太阳位置；

S82、紧接着计算感兴趣区域内每一点r_i的太阳瞬时直接辐射I_b(KWh/m²)为：

其中，f_i(r_i)为该点的瞬时遮挡状况，α为当时的太阳高度角，然后，计算散射辐射透过率τ_d，从而计算太阳散射辐射I_d(KWh/m²)。计算如下：

I_d＝I_o>d

S83、最后，计算该点瞬时总太阳辐射I_g(KWh/m²)，有I_g＝I_b+I_d，对于任意时间段的单点r_i的太阳能总辐射可以通过对每一小时的总辐射量累加获得，即：

其中，t_begin和t_end为太阳能模拟的起始和截止时间，并将计算结果可视化显示于用户界面。

采用上述技术方案后，本发明与背景技术相比，具有如下优点：

1、本发明采用了GHPR算法，充分利用了点云的三维信息，对于建筑物或植物的延展性结构具有很好的适应性，因此更能适应环境复杂的调查场景；与采用建模的方法相比，本发明直接基于点云运算，避免了人工交互操作和模型建模中的数据失真，有效促进了资源评估流程的自动化，规范化；与传统评估方法相比，本发明可以为用户提供直观，定量，全局的太阳能资源分布结果，该方法易于保存并可重复模拟；

2、本发明采用全局阴影计算并将计算范围限定在感兴趣区域内，从而加快了计算速度，有利于实施现场太阳能资源评估，避免建模过程中次要物体(如窗户，栏杆等)丢失引起的误差，可以为不同行业的太阳能资源考察者提供一步到位的分析结果；

3、本发明提出了一种点云场景中的光源偏差控制方法，用来约束和控制不同光源间产生的阴影偏差，并利用最远点贪心策略快速获取太阳视点数量和位置的近似解，不仅减少偏差，而且提高了效率。

附图说明

图1为本发明实施例的流程示意图；

图2(a)为本发明实施例原始点云数据；(b)为感兴趣区域点集的提取结果；(c)为实施例偏差控制结果；

图3为本发明实施例太阳位置计算结果；

图4为实施例日照遮挡模拟结果：(a)为本方法模拟的结果；(b)为模拟结果与对应时间的真实情况对比(标准结果)；

图5为实施例总辐射模拟结果：(a)为单季度模拟结果(实施例为秋季)；(b)为全年模拟结果。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

实施例

本实施例的流程示意图可见图1，具体操作过程如下：

1.调查区域点集提取

已知三维地面激光点云P_raw的原始点云数据如图2(a)所示，先对点云进行下采样(抽稀)，即减少点云中数据点的数量。通常情况下，地面激光点云的数据精度对于光照计算过于冗余。因此，需要先设置一个抽稀距离L_s，使得抽稀后的点云间距可以保持在L_s左右且点云结构保持不变。本方法采用具有局部结构保持性的体素抽稀方法，该方法首先将点云划分为边长为L_s的正立方体体素，再按体素内点集的中心点作为新的数据点，这些新的点集形成了抽稀后的点云P。

接着，对抽稀后的点云P进行感兴趣区域R的提取，感兴趣区域提取目的是获取调查区域内的数据点集，本方法主要提取的是区域内的地面点，用于将太阳能资源分布结果绘制于地面点上，从而形成太阳能资源分布图。因此，需要采用地面点滤波排除点云中的非地面目标。对于不同的场景，如建筑物天台或是道路，它们的地面点在形态上仍然具有共性，即地面点区域内梯度变化相对平缓，且在垂直方向上没有剧烈的阶跃，因此本方法应用基于体素高程生长滤波的方法来提取地面点。

首先，将点云按水平面固定边长粗划分为一系列方形网格区域，并记录整个点云的高程最低点。再对每个水平网格的点集进行一次三维细分，即将网格内的点集按固定边长划分为立方体素，记录每个体素的中心高程值。基于整个点云的最低高程值z_min设置一个阈值高度t_z，该高度根据感兴趣区域的海拔高度而定，也可以根据不同调查对象的基本高度来设置高程阈值t_z。通常激光扫描仪架设于调查区域之上，所以该高度可以通过扫描仪在场景中的高程值加3米获得。如果一个体素的中心高程与最低点的相对高度高于该阈值t_z，则认为其过高，不可能是地面点。若一个低于阈值的体素可以生长到该阈值高度t_z，则认为是墙体或垂直物体的一部分，也为非地面点。生长规则是其体素的上面九个邻域体素存在非空，并迭代查询非空体素的上九邻域体素是否有非空体素，直到某查询的非空体素高于阈值高度t_z。最后，保留地面体素内的数据点，形成地面点集。

接着，对地面点集进行聚类，即将离散点集聚类成一系列对象。本方法采用欧式距离聚类，对获得的聚类簇标记序号。接着，计算每一个类簇的平均密度值。其理论依据是由地面激光扫描仪的作业方式决定的，即离激光扫描仪近的区域扫描密度非常大，通常可以达到远处目标密度的10倍以上。平均密度值是该类簇每个点在原始点云中一定邻域(0.2米)内点数量的均值。最后，将类簇按密度值排序，选择密度最高的一个作为感兴趣区域(调查对象)。返回抽稀后的点云，根据各点索引值标记感兴趣区域的数据点r_i，r_i∈R。其中，感兴趣区域地面点集提取结果如图2(b)所示。

2.基点计算及偏差控制

本发明采用点光源模拟太阳光的方法，因此需要进行光源偏差控制，即控制平行光和点光源之间的阴影偏差。此外，在进行太阳位置计算前，需要给定基点，即基于某一点计算太阳的位置。对于保存后的感兴趣区域点集R，我们提出的偏差控制方法思想是选择感兴趣区域内的一些点作为基点o，每个基点计算一次太阳位置，由每个基点计算出的太阳位置也将作为基点一定范围内邻域点集的太阳参考位置，从而减小阴影绘制时的偏差。因此，设基点点集为O，o_k∈O，太阳位置点集为S，其中S＝{s_k|k∈[1,n]}，每个s_k相应的感兴趣子区域点集为G_k，满足∪_k∈[1,n]G_k＝R且∩_k∈[1,n]G_k＝φ。

首先设置约束的偏差角大小θ，即模拟光线偏差平行光的最大允许角度，从而可以确定两个基点之间的最大距离l，也将由此为依据对提取的感兴趣区域进行划分。此外，还需要设置太阳距离d，本方法中的太阳距离为在点云场景中的太阳位置s_k离感兴趣区域相应基点o_k的距离，有||o_k-s_k||＝d。以上两者关系如下：

l＝d/tan(π/2-θ)

在通常场景中，远距离物体对遮挡的影响很小，因此θ通常设置2°，d设置600米。下一步，进行基点o的个数和位置的计算。本方法提出采用最远点贪心策略来快速计算基点的位置和个数，该子程序的目标是在有限的时间内，尽可能地以最少的基点数n_min使得整个基点集的作用范围覆盖感兴趣点集，即即确保感兴趣区域内的每个点r_i存在一个基点o_k使其满足条件||r_i-o_k||≤l。对于每个划分区域G_k，都有G_k＝{r_i∈R|||o_k-r_i||≤||o_t-r_i||,k≠t且t∈[1,n],i∈[1,n_R]}。n_R为点集R的势。

最远点策略随机选择一个感兴趣区域点作为初始基点o₁，然后构建基点序列O＝{o₁}。然后计算感兴趣区域点集上与该序列最远的点，即并将其添加到O中，为O＝{o₁,o₂}。直到任意一点r_i，r_i∈R，都满足获得基点序列后O＝{o_k|k∈[1,n]}，对于感兴趣区域点集中的每个点r_i，求取其在基点集上的最近邻点，同属于一个基点的感兴趣区域点构成一个偏差控制区域G_k。图2(c)为本实施例的偏差控制结果，其中参数设置为d＝400，θ＝2°，小黑点表示基点，分别对应基点划分的三个负责区域。

3.太阳在点云场景中的位置计算

对于每个基点o_k，我们通过一套太阳位置求解模型，来计算其相应的太阳位置s_k。该模型的输入为该基点的经纬度和当地太阳时。因此，第一个步为获取基点的经纬度坐标，本方法首先将每个基点从WGS84坐标额外转为一份经纬度坐标，有[o_Ψ,o_Φ,o_h]＝T([o_x,o_y,o_z])。其中，T为坐标变换，[o_Ψ,o_Φ]为基点的经纬度值，o_h为海拔高度值。接着，设置太阳时间，根据日期计算赤纬角：

δ＝23.45°sin(360(n_d-80)/365)

其中，n_d为查询日期所在全年中的天数。再根据设置的具体时间(时分秒)计算太阳时角：

ω＝15°(12-T)

其中T为查询的具体时间，需要注意为当地太阳时间。从而获得太阳方位角ψ和高度角α，有：

sinα＝sinδsino_Φ+cosδcoso_Φcosω

然后，基于太阳方位角ψ，太阳高度角α和我们之前设置的太阳距离d，则基点o的太阳位置s，s＝[o_Ψ+f_d(y),o_Φ-f_d(x),o_h+z]^T，f_d为将长度化为度数的函数，在经纬度坐标下可以如下计算:

x＝d sin(90°-α)cos(ψ)

y＝d sin(90°-α)sin(ψ)

z＝d cos(90°-α)

最后，我们将经纬度和高度坐标转化WGS84坐标，获得太阳在点云场景三维空间中的相对位置，s＝T′(s)，T′为T的逆变换。对于每个基点o_k，我们计算一个s_k，构成了太阳点集，S＝{s_k|k∈[1,n]}，在空间上形成了太阳点阵。图3为本实施例太阳位置计算结果，表示本实施例模拟2017年冬至日(12月22日)点云场景上空的太阳视运动轨迹(小球颜色从明到暗为从6点到18点的太阳位置)。

4.遮挡分析和阴影绘制

我们应用GHPR算法绘制日照阴影，从而进行日照模拟计算。首先，遍历太阳点集S，对于每个太阳位置s_i，以及调查场景抽稀后的三维点集P，P＝{p_i|i∈[1,m]}，其中m是点集的势，执行径向几何变换F，有:

其中，若p_i＝s，则变换后的点集为而f_k是变换的核函数:

f_k(||p_i-s||)＝2λmax_pi∈P||p_i-s||-||p_i-s||

其中λ是径长放大参数，有λ≥1。然后对点集做一次凸包计算，并保留凸包内的点为阴影点集W。

在遍历的过程中，对于每个太阳位置s_k，我们计算一次W_k。然后，我们取W_k与负责区域G_k的交集为各相应区域内的阴影点。遍历完后，我们将每个区域内的的阴影结果取并集，从而得到整个感兴趣区域的光照结果。

最后一步是将光照结果二值化，即将阴影点标记为0，记为f_i(r_i)＝0，而非阴影点标记为1，记为f_i(r_i)＝1，因为标记的阴影只存在于感兴趣区域，感兴趣区域外可不做计算，一律置零，从而减少交集运算时间。本实施例日照遮挡模拟结果如图4所示，其中图4(a)为本方法模拟的结果，图4(b)为模拟结果与对应时间的真实情况对比(标准结果)。

5.太阳能辐射模型计算

由于地面激光扫描点云本身在前期数据采集阶段，可以通过GPS等定位设备获取大地坐标WGS84，因此可以直接执行太阳能辐射估计。太阳能辐射评估需要建立在太阳能辐射模型和遮挡分析结果的基础上。对于太阳能辐射模型，本文主要使用Hottel模型。输入太阳时，可以通过以下公式获取任意一天大气层外的太阳辐照度：

I_o＝I_s[1+0.033cos(360°n/365)]

其中，I_s为太阳常数，即大气层外与辐射方向垂直时所接收的太阳辐照度，有I_s＝1353W/m²。则点云场景中感兴趣区域点集中某点r_i的太阳瞬时直接辐射I_b(KWh/m²)为：

其中，τ_b为直射辐射透过率，ε为斜面修正因子，f_i(r_i)为该点的瞬时遮挡状况。此外，还存在太阳散射辐射I_d(KWh/m²)。其关系式如下：

I_d＝I_o>d

其中τ_d为散射辐射透过率。直射辐射和散射辐射透过率的计算如下：

τ_b＝a₀+a₁e^-k/cosα

τ_d＝0.2710-0.293τ_b

其中，对于τ_b的参数可以通过如下计算：

a₀＝(0.4327-0.00821(o_h-6)²)r₀

a₁＝(0.5055+0.00595(o_h-6.5)²)r₁

k＝(0.2711+0.01858(o_h-2.5)²)r_k

其中，[r₀,r₁,r_k]为所在地的气候修正，该三个参数可以通过建立经纬度气候查询表自动获得。例如，以厦门亚热带地区为例，其三个参数分别为r₀＝0.95，r₁＝0.98，r_k＝1.02。

接着，计算太阳直接辐射的瞬时倾斜修正因子，该步骤首先需要考虑入射斜面的斜率，在太阳能资源评估中考虑倾斜因素才能使结果更接近真实的情况。因此，我们求取太阳入射光学与斜面法向量之间的夹角，从而进行倾斜因素修正。首先，是获得法向量，这一步可以在提取感兴趣区域后就进行。我们遍历每个感兴趣区域的点，每遍历到一个感兴趣区域点r_i，就用kd-tree搜索其邻域，从而构建邻域点集R'，R'＝{r_i'|i∈[1,K_r]}。在这个点集的基础上，先计算邻域点集的重心再构建协方差矩阵M。

然后对协方差矩阵进行矩阵分解，获取协方差矩阵的特征向量和特征值矩阵。提取最小特征根对应的单位向量为法向量该法向量就是局部曲面的法向量。

接着，连接太阳和当前遍历点构建入射向量为s-r_i。校正因子ε可以通过入射向量与法向量之间的夹角余旋计算，关系式计算如下:

其中·为向量点乘，而s为该点相关基点所对应的太阳位置。

接着，计算该点瞬时总太阳辐射I_g(KWh/m²)，有I_g＝I_b+I_d，则对于任意时间段的单点太阳能总辐射可以通过对每一小时的总辐射量累加获得，即：

其中，t_begin和t_end为太阳能模拟的起始和截止时间，按需求可以是数小时，数天甚至数年。

最后，将计算结果可视化显示于用户界面，通过不同色彩渲染来显示各个点之间的辐射值差异，从而为用户提供可视化、定量化、全局化的太阳能潜力评估结果。图5为本实施例总辐射模拟结果，其中图5(a)为单季度模拟结果(秋季)，图5(b)为全年模拟结果，数据点颜色从浅到深表示太阳能资源分布结果中从最小数值到最大数值的辐射值，大图为俯视图视角，右上角黑框中的小图为侧视图视角，右侧的色域栏显示了辐射值(KWh/m²)的变化范围。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。