基于XGB算法及全天空成像仪图像的超短期DNI预测

胡 中，倪杭飞，李 伟，宓霄凌，白 帆，王伊娜

(浙江中控太阳能技术有限公司，杭州 310053)

0 引言

在塔式太阳能热发电站中，由定日镜镜场将太阳辐射反射到吸热器上，然后吸热器将接收的太阳辐射能转化成热能，为后续的太阳能热发电提供能量。在该过程中，云层对定日镜场上方太阳辐射的遮挡是影响吸热器接收的太阳直接辐射(DNI)的主要因素，当云层遮挡定日镜上方的太阳辐射时，定日镜镜场反射的太阳辐射减少，吸热器接收的DNI也相应下降；而当云层离开时，定日镜镜场反射的太阳辐射增大，吸热器接收的DNI也随之升高。因此，云层短时间的遮挡和离开会使吸热器接收的DNI发生随机波动，造成吸热器表面的受热情况发生剧烈变化。而通过超短期(0～30 min)DNI预测可以提前对定日镜镜场进行调控，从而解决这一问题。

基于上述原因，本文提出了一种基于极值梯度提升(eXtreme Gradient Boosting，XGB)算法及全天空成像仪图像的超短期DNI预测方法。通过使用基于太阳位置的全天空成像仪标定方法来完成全天空成像仪的姿态校正和畸变矫正；然后对晴空背景参数进行拟合并建立晴空库，为云层的检测提供基础；根据云的运动矢量，利用外推法得到未来云层可能会遮挡太阳的图像点，并提取图像特征，用于基于XGB算法的DNI超短期预测模型的训练，最后利用训练后的模型完成超短期DNI预测。

1 全天空成像仪

本文选用型号为EKO SRF-02的全天空成像仪，其分辨率为2274×1702，视野大于180°，可实现全天空监控，其实物图及安装效果图分别如图1、图2所示。将该全天空成像仪设置为每40 s采集2帧图像，其中的一帧图像为自动曝光，另一帧图像为欠曝光。本文后续均基于欠曝光图像进行处理。

图1 型号为EKO SRF-02的全天空成像仪实物图Fig.1 Physical image of all-sky imager model EKO SRF-02

图2 全天空成像仪的安装效果图Fig.2 Installation effect image of all-sky imager

全天空成像仪在安装时需进行水平调准，其水平倾斜度不能超过0.5°；由于全天空成像仪使用的是鱼眼镜头，会存在畸变，因此需要进行标定。

本文开发了一种基于太阳位置的全天空成像仪标定方法。该方法不需要通过标定板对全天空成像仪进行标定，而是通过采集一整天的晴空图像，根据图像上太阳的中心位置，并结合理论太阳高度角及太阳方位角计算全天空成像仪的姿态参数及畸变参数，从而实现全天空成像仪的姿态校正与畸变矫正，减少了现场工作量。

2 建立晴空库与云层的检测

2.1 建立晴空库

在对云层进行检测之前，需要建立晴空库。最好在夏季采集晴空图像，因为夏季每天的太阳高度角的最大值接近一年中太阳高度角的最大值；同时，需选择大气透射率最高的天气。

同基于太阳位置的全天空成像仪标定方法一样，建立晴空库时也需要一整天的晴空图像，然后分别对每一帧图像中的红(R)、绿(G)、蓝(B)通道进行背景拟合，从而得到背景拟合参数。该背景拟合参数对应该帧图像的太阳高度角，即每帧图像的太阳高度角对应一组拟合参数，据此建立晴空库。在进行云层分割时，可以根据当前时刻的实际太阳高度角找到对应的拟合参数(选择与晴空库中太阳高度角差值最小时对应的拟合参数)，并根据该拟合参数生成晴空背景图像。

根据文献[1]，将晴空下天空亮度L建模为关于PZA(像点/天顶角，即像点对应视线和天顶方向的夹角)和SPA(太阳/像点角，即像点对应视线和太阳方向的夹角)的函数，具体可表示为：

同时，根据文献[1]中基于实测数据的统计分析及建立的天空亮度仿真模型，可将晴空背景模型(下文简称为“晴空模型”)拟合公式表示为：

式中：Ip(PZA,SPA)为天空亮度分布函数；φ1(PZA)为天空亮度渐变函数，简化记作F1；φ2(SPA)为散射指标函数，简化记作F2。

其中：

式中：a1、a2、a3均为天空亮度渐变函数的拟合参数。

式中：a4、a5、a6、a7均为散射指标函数的拟合参数。

在利用式(2)对全天空成像仪图像进行晴空背景拟合的过程中，发现该公式在部分情况下无法对太阳周边像素进行精确拟合。

将式(2)两侧各减去φ1(PZA)，可得到：

由式(5)可知，Ip(PZA,SPA)–F1与F1呈线性关系。

B通道的Fg_b、Bg_b、F1及Fg_b–F1生成的曲线如图3所示。其中：Fg_b为前景图像B通道的值；Bg_b为背景图像B通道的值。

图3 B通道的Fg_b、Bg_b、F1及Fg_b–F1生成的曲线Fig.3 Curves generated by Fg_b，Bg_b，F1 and Fg_b–F1 of channel B

从图3中可以看出，Fg_b–F1与F1并非呈线性关系。这与根据式(5)得出的结论不相符，由此可知，式(5)在本文中并不适用。因此，本文参考式(1)、式(2)并结合图3，修改了晴空模型，假设其是φ1(PZA)和φ2(SPA)的和，并认为这2个函数的值独立且无关联，从而得到了新的拟合公式，即：

新拟合公式中用加法替代了式(2)(旧拟合公式)中的乘法，大幅降低了拟合计算过程的复杂度，使晴空模型得到了简化。

其中：

式中：b1～b3均为新拟合公式的天空亮度渐变函数的拟合参数。

式中：b4～b8均为新拟合公式的散射指标函数的拟合参数。

利用新的晴空模型对晴空时的全天空成像仪图像进行拟合，得到拟合后的晴空背景图像，如图4所示，以及对应的R通道的和的拟合结果，如图5所示。

图4 全天空成像仪图像与拟合得到的晴空背景图像Fig.4 Image of all-sky imager and processed clear sky background

图5 利用新拟合公式得到的R通道的F′1、F′2值的拟合结果Fig.5 Fitting result of F′1 and F′2 values of R channel obtained by new fitting formula

对比了分别采用新、旧拟合公式得到的R、G、B通道的F′2的拟合结果，具体如图6～图8所示。

图6 通过新、旧拟合公式得到的R通道的F′2值和F2值的拟合结果Fig.6 Fitting results of F′2 values and F2 values of R channel obtained by new and old fitting formulas

图7 通过新、旧拟合公式得到的G通道的F′2值和F2值的拟合结果Fig.7 Fitting results of F′2 values and F2 values of G channel obtained by new and old fitting formulas

图8 通过新、旧拟合公式得到的B通道的F′2值和F2值的拟合结果Fig.8 Fitting results of F′2 values and F2 values of B channel obtained by new and old fitting formulas

由图6～图8可知，针对本文所述型号的全天空成像仪拍摄的图像，在太阳高度角较小时，利用新拟合公式可以得到更好的拟合结果。

2.2 云层的检测

云层的检测范围为天顶角80°以内的区域。检测前需要对该区域内的地面建筑进行裁剪并进行插值修复，这是因为该区域接近地平面，镜头畸变最为严重，沿圆弧方向进行插值修复可以得到更好的插值结果。然后利用图像变换公式分别计算前景NRBR(即归一化的RB比值，此处为便于分析将值域映射到[0, 1])图像和晴空背景NRBR图像，并计算二者的差值图像，即ΔNRBR图像。NRBR的计算式可表示为：

由于晴空背景亮度的不均匀性，本文对ΔNRBR图像进行了归一化。

在ΔNRBR图像上的太阳及其临近区域会形成空洞，如图9a所示，因此本文对该区域进行了插值修复，插值范围选择SPA＜30°的区域，沿圆弧方向进行线性插值，修复后的图像如图9b所示，其可用于准确分割太阳及其临近区域的云层。

图9 未处理与处理后的ΔNRBR图像Fig.9 Images of unprocessed and processed ΔNRBR

使用固定阈值对归一化的差值图像进行二值化，将较厚的云层(下文简称为“厚云”)的分割阈值设为0.85，薄云的分割阈值设为0.5，此处的分割阈值为差值图像灰度映射到[0, 1]后的取值。云的分割结果如图10所示。

图10 云分割前、后的图像Fig.10 Images of before and after cloud segmentation

3 云的运动矢量计算

在对DNI进行预测前需要先对云的运动矢量进行计算。云的运动矢量是指云的整体运动方向和速度。本文使用归一化相关系数匹配算法对全天空成像仪拍摄的前、后2帧图像进行匹配(匹配前需先对图像进行矫正，且只对天顶角80°范围内的像素进行匹配，因为该区域的像素失真较小)，并计算云的运动矢量，计算结果如图11所示。图中：绿线为相关系数匹配结果；红线为计算得到的云的运动矢量，其中，线段长度代表运动速度，箭头代表运动方向。

本文使用卡尔曼滤波对不同时刻计算得到的云的运动矢量进行了滤波，使计算结果更为稳定。

图11 云的运动矢量计算结果Fig.11 Calculation result of cloud motion vector

由于云的运动矢量是后续进行云层分类及超短期DNI预测特征提取的基础，因此，根据云的运动矢量，利用外推法计算未来30 min内(时间间隔取1 min)云层可能会遮挡太阳的30个图像点。从上述30个图像点提取图像特征，进行超短期DNI预测模型的训练与DNI实时预测。

本文的假设是：在短时间(0～30 min)内云的形态保持不变，且云的运动形式为整体平移。但事实上，云的变化方式多种多样，且目前暂无可以对云的变化进行精确预测的方法，这也是影响DNI预测结果精度的主要原因之一。

4 云层的分类

云层的光学厚度是决定其对太阳辐射衰减产生影响的主要因素。本文利用XGB算法对云层的光学厚度进行了分类，具体为：晴空、气溶胶、薄云、可能为厚云、厚云、蓝色的厚云(下文简称为“蓝厚云”)。由于气溶胶是影响DNI预测结果精度的主要因素之一，因此将其单独作为一类进行分类。此外，部分云用肉眼很难判断是薄云或厚云，因此添加了“可能为厚云” 这一类型。

在提取图像特征之前对太阳及其邻近区域进行了插值处理，这样有助于得到较好的DNI预测结果。

在保证数据准确性的基础上，本文采用半人工方式提取了一定数量的云层分类特征样本，用于云层分类模型的训练，并对云层分类模型进行了验证。不同类型云层的分类结果占比如表1所示。

表1 不同类型云层的分类结果占比Table 1 Proportion of classification results of different types of cloud layers

有时会存在某些云层分类错误的情况，比如：可能为厚云被归类为厚云，蓝厚云被归类为厚云等，但这是符合预期的。在这种情况下，云层分类的预期范围准确率可达90%以上，如表2所示。

表2 云层分类的预期范围准确率Table 2 Accuracy of expected range of cloud classification

将通过云层分类训练后得到的云层分类模型用于DNI预测时的云层分类，将分类后得到的云层类型作为超短期DNI预测模型的特征之一。

5 超短期DNI预测

5.1 基于XGB算法的超短期DNI预测模型训练

根据之前计算得到的云的运动矢量，利用外推法获得未来0～30 min内(间隔1 min)云层可能会遮挡太阳的30个图像点，然后对这30个图像点进行图像特征提取，并利用云层分类模型完成云层分类。以通过云层分类提取到的图像特征和云层类型分类结果作为超短期DNI预测的特征。特征提取过程完全由软件自动完成。

全天空成像仪旁边安装有DNI测量仪器，用于采集DNI数据。受云层和其他因素的影响，DNI会产生衰减，实测DNI的衰减率DDNI可表示为：

式中：DNIt为理论的DNI值；DNIm为实测的DNI值。

建立超短期DNI预测模型，并将DDNI作为其输出值。该模型的训练数据选用2018～2019年期间收集的所有多云天气中以天为单位筛选出的云的形态较为稳定时的数据，然后根据云层类型分类结果和实测DNI的衰减率对筛选出的数据进行进一步过滤，过滤掉云层类型和实测DNI的衰减率明显不符的数据，使通过训练得到的超短期DNI预测模型更为稳定。最后利用XGB算法完成超短期DNI预测模型的训练。

不同云层类型时DDNI的取值范围如表3所示。

表3 不同云层类型时DDNI的取值范围Table 3 Value range of DDNI for different cloud layers types

5.2 基于XGB算法的超短期DNI预测

超短期DNI预测时同样根据云的运动矢量，采用外推法获得未来0～30 min内(间隔1 min)云层可能会遮挡太阳的30个图像点，对应30个时刻，并提取这30个时刻的特征；然后利用超短期DNI预测模型对这30个时刻的DNI进行预测。

相较于实测DNI，利用超短期DNI预测模型得到不同天气(晴天、多云、阴天)下提前10 min时的预测DNI曲线，并将该预测曲线与理论DNI曲线和实测DNI曲线进行对比，具体如图12所示。

从图12中可以看出，不同天气情况下，超短期DNI预测模型得到的提前10 min的预测DNI与实测DNI曲线较为吻合。

针对2018～2019年中的152天且每天太阳高度角大于5°的时段，将基于XGB算法的超短期DNI预测模型得到的所有1～30 min 的DNI预测结果进行统计分析，具体如表4所示。

图12 不同天气下超短期DNI预测模型提前10 min的预测DNI曲线与理论及实测DNI曲线的对比Fig.12 Comparison between predicted DNI curve of ultrashort-term DNI prediction model 10 minutes ahead of time and theoretical and actual DNI curves under different weather conditions

表4 不同提前预测时下DNI预测结果的MAE与RMSETable 4 MAE and RMSE of DNI prediction results under different advance prediction time

(续表)

从表4可以看出，不同提前预测时间的总体DNI预测值的平均绝对误差(MAE)为116.0 W/m2，均方根误差(RMSE)为203.0 W/m2。在未使用可持续模型的情况下，从这2个值来看，预测结果的精度较高。

将实测的DNI数据与超短期DNI预测模型的预测结果进行对比，当以“某个时刻的DNI值是否大于或小于400 W/m2”作为判断标准时，超短期DNI预测模型的预测结果的正确率均大于85%，符合塔式太阳能热发电站运行过程中对DNI预测的指标要求。而通过进一步分析发现，影响DNI预测结果精度的主要因素为气溶胶及云的不规则变化。

6 结论

本文未使用可持续模型，而是采用基于XGB算法及全天空成像仪图像实现了超短期(0～30 min)的DNI预测。通过使用基于太阳位置的全天空成像仪标定方法，完成了全天空成像仪的姿态校正和畸变矫正；准确拟合了晴空背景参数并建立晴空库，为云层的检测提供了基础；通过对太阳及其临近区域进行插值修复，使1 min内的DNI预测成为可能。根据云的运动矢量，利用外推法得到未来云层可能会遮挡太阳的图像点，并提取图像特征，用于基于XGB算法的超短期DNI预测模型的训练，最终以训练后的该模型完成超短期DNI预测，并得到了较好的结果，符合塔式太阳能热发电站运行过程中对DNI预测的指标要求。