基于XGB算法及全天空成像儀圖像的超短期DNI預測

胡 中，倪杭飛，李 偉，宓霄凌，白 帆，王伊娜

(浙江中控太陽能技術(shù)有限公司，杭州 310053)

0 引言

在塔式太陽能熱發(fā)電站中，由定日鏡鏡場將太陽輻射反射到吸熱器上，然后吸熱器將接收的太陽輻射能轉(zhuǎn)化成熱能，為后續(xù)的太陽能熱發(fā)電提供能量。在該過程中，云層對定日鏡場上方太陽輻射的遮擋是影響吸熱器接收的太陽直接輻射(DNI)的主要因素，當云層遮擋定日鏡上方的太陽輻射時，定日鏡鏡場反射的太陽輻射減少，吸熱器接收的DNI也相應下降；而當云層離開時，定日鏡鏡場反射的太陽輻射增大，吸熱器接收的DNI也隨之升高。因此，云層短時間的遮擋和離開會使吸熱器接收的DNI發(fā)生隨機波動，造成吸熱器表面的受熱情況發(fā)生劇烈變化。而通過超短期(0～30 min)DNI預測可以提前對定日鏡鏡場進行調(diào)控，從而解決這一問題。

基于上述原因，本文提出了一種基于極值梯度提升(eXtreme Gradient Boosting，XGB)算法及全天空成像儀圖像的超短期DNI預測方法。通過使用基于太陽位置的全天空成像儀標定方法來完成全天空成像儀的姿態(tài)校正和畸變矯正；然后對晴空背景參數(shù)進行擬合并建立晴空庫，為云層的檢測提供基礎；根據(jù)云的運動矢量，利用外推法得到未來云層可能會遮擋太陽的圖像點，并提取圖像特征，用于基于XGB算法的DNI超短期預測模型的訓練，最后利用訓練后的模型完成超短期DNI預測。

1 全天空成像儀

本文選用型號為EKO SRF-02的全天空成像儀，其分辨率為2274×1702，視野大于180°，可實現(xiàn)全天空監(jiān)控，其實物圖及安裝效果圖分別如圖1、圖2所示。將該全天空成像儀設置為每40 s采集2幀圖像，其中的一幀圖像為自動曝光，另一幀圖像為欠曝光。本文后續(xù)均基于欠曝光圖像進行處理。

圖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

全天空成像儀在安裝時需進行水平調(diào)準，其水平傾斜度不能超過0.5°；由于全天空成像儀使用的是魚眼鏡頭，會存在畸變，因此需要進行標定。

本文開發(fā)了一種基于太陽位置的全天空成像儀標定方法。該方法不需要通過標定板對全天空成像儀進行標定，而是通過采集一整天的晴空圖像，根據(jù)圖像上太陽的中心位置，并結(jié)合理論太陽高度角及太陽方位角計算全天空成像儀的姿態(tài)參數(shù)及畸變參數(shù)，從而實現(xiàn)全天空成像儀的姿態(tài)校正與畸變矯正，減少了現(xiàn)場工作量。

2 建立晴空庫與云層的檢測

2.1 建立晴空庫

在對云層進行檢測之前，需要建立晴空庫。最好在夏季采集晴空圖像，因為夏季每天的太陽高度角的最大值接近一年中太陽高度角的最大值；同時，需選擇大氣透射率最高的天氣。

同基于太陽位置的全天空成像儀標定方法一樣，建立晴空庫時也需要一整天的晴空圖像，然后分別對每一幀圖像中的紅(R)、綠(G)、藍(B)通道進行背景擬合，從而得到背景擬合參數(shù)。該背景擬合參數(shù)對應該幀圖像的太陽高度角，即每幀圖像的太陽高度角對應一組擬合參數(shù)，據(jù)此建立晴空庫。在進行云層分割時，可以根據(jù)當前時刻的實際太陽高度角找到對應的擬合參數(shù)(選擇與晴空庫中太陽高度角差值最小時對應的擬合參數(shù))，并根據(jù)該擬合參數(shù)生成晴空背景圖像。

根據(jù)文獻[1]，將晴空下天空亮度L建模為關于PZA(像點/天頂角，即像點對應視線和天頂方向的夾角)和SPA(太陽/像點角，即像點對應視線和太陽方向的夾角)的函數(shù)，具體可表示為：

同時，根據(jù)文獻[1]中基于實測數(shù)據(jù)的統(tǒng)計分析及建立的天空亮度仿真模型，可將晴空背景模型(下文簡稱為“晴空模型”)擬合公式表示為：

式中：Ip(PZA,SPA)為天空亮度分布函數(shù)；φ1(PZA)為天空亮度漸變函數(shù)，簡化記作F1；φ2(SPA)為散射指標函數(shù)，簡化記作F2。

其中：

式中：a1、a2、a3均為天空亮度漸變函數(shù)的擬合參數(shù)。

式中：a4、a5、a6、a7均為散射指標函數(shù)的擬合參數(shù)。

在利用式(2)對全天空成像儀圖像進行晴空背景擬合的過程中，發(fā)現(xiàn)該公式在部分情況下無法對太陽周邊像素進行精確擬合。

將式(2)兩側(cè)各減去φ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，F(xiàn)1 and Fg_b–F1 of channel B

從圖3中可以看出，F(xiàn)g_b–F1與F1并非呈線性關系。這與根據(jù)式(5)得出的結(jié)論不相符，由此可知，式(5)在本文中并不適用。因此，本文參考式(1)、式(2)并結(jié)合圖3，修改了晴空模型，假設其是φ1(PZA)和φ2(SPA)的和，并認為這2個函數(shù)的值獨立且無關聯(lián)，從而得到了新的擬合公式，即：

新擬合公式中用加法替代了式(2)(舊擬合公式)中的乘法，大幅降低了擬合計算過程的復雜度，使晴空模型得到了簡化。

其中：

式中：b1～b3均為新擬合公式的天空亮度漸變函數(shù)的擬合參數(shù)。

式中：b4～b8均為新擬合公式的散射指標函數(shù)的擬合參數(shù)。

利用新的晴空模型對晴空時的全天空成像儀圖像進行擬合，得到擬合后的晴空背景圖像，如圖4所示，以及對應的R通道的和的擬合結(jié)果，如圖5所示。

圖4 全天空成像儀圖像與擬合得到的晴空背景圖像Fig.4 Image of all-sky imager and processed clear sky background

圖5 利用新擬合公式得到的R通道的F′1、F′2值的擬合結(jié)果Fig.5 Fitting result of F′1 and F′2 values of R channel obtained by new fitting formula

對比了分別采用新、舊擬合公式得到的R、G、B通道的F′2的擬合結(jié)果，具體如圖6～圖8所示。

圖6 通過新、舊擬合公式得到的R通道的F′2值和F2值的擬合結(jié)果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值的擬合結(jié)果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值的擬合結(jié)果Fig.8 Fitting results of F′2 values and F2 values of B channel obtained by new and old fitting formulas

由圖6～圖8可知，針對本文所述型號的全天空成像儀拍攝的圖像，在太陽高度角較小時，利用新擬合公式可以得到更好的擬合結(jié)果。

2.2 云層的檢測

云層的檢測范圍為天頂角80°以內(nèi)的區(qū)域。檢測前需要對該區(qū)域內(nèi)的地面建筑進行裁剪并進行插值修復，這是因為該區(qū)域接近地平面，鏡頭畸變最為嚴重，沿圓弧方向進行插值修復可以得到更好的插值結(jié)果。然后利用圖像變換公式分別計算前景NRBR(即歸一化的RB比值，此處為便于分析將值域映射到[0, 1])圖像和晴空背景NRBR圖像，并計算二者的差值圖像，即ΔNRBR圖像。NRBR的計算式可表示為：

由于晴空背景亮度的不均勻性，本文對ΔNRBR圖像進行了歸一化。

在ΔNRBR圖像上的太陽及其臨近區(qū)域會形成空洞，如圖9a所示，因此本文對該區(qū)域進行了插值修復，插值范圍選擇SPA＜30°的區(qū)域，沿圓弧方向進行線性插值，修復后的圖像如圖9b所示，其可用于準確分割太陽及其臨近區(qū)域的云層。

圖9 未處理與處理后的ΔNRBR圖像Fig.9 Images of unprocessed and processed ΔNRBR

使用固定閾值對歸一化的差值圖像進行二值化，將較厚的云層(下文簡稱為“厚云”)的分割閾值設為0.85，薄云的分割閾值設為0.5，此處的分割閾值為差值圖像灰度映射到[0, 1]后的取值。云的分割結(jié)果如圖10所示。

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

3 云的運動矢量計算

在對DNI進行預測前需要先對云的運動矢量進行計算。云的運動矢量是指云的整體運動方向和速度。本文使用歸一化相關系數(shù)匹配算法對全天空成像儀拍攝的前、后2幀圖像進行匹配(匹配前需先對圖像進行矯正，且只對天頂角80°范圍內(nèi)的像素進行匹配，因為該區(qū)域的像素失真較小)，并計算云的運動矢量，計算結(jié)果如圖11所示。圖中：綠線為相關系數(shù)匹配結(jié)果；紅線為計算得到的云的運動矢量，其中，線段長度代表運動速度，箭頭代表運動方向。

本文使用卡爾曼濾波對不同時刻計算得到的云的運動矢量進行了濾波，使計算結(jié)果更為穩(wěn)定。

圖11 云的運動矢量計算結(jié)果Fig.11 Calculation result of cloud motion vector

由于云的運動矢量是后續(xù)進行云層分類及超短期DNI預測特征提取的基礎，因此，根據(jù)云的運動矢量，利用外推法計算未來30 min內(nèi)(時間間隔取1 min)云層可能會遮擋太陽的30個圖像點。從上述30個圖像點提取圖像特征，進行超短期DNI預測模型的訓練與DNI實時預測。

本文的假設是：在短時間(0～30 min)內(nèi)云的形態(tài)保持不變，且云的運動形式為整體平移。但事實上，云的變化方式多種多樣，且目前暫無可以對云的變化進行精確預測的方法，這也是影響DNI預測結(jié)果精度的主要原因之一。

4 云層的分類

云層的光學厚度是決定其對太陽輻射衰減產(chǎn)生影響的主要因素。本文利用XGB算法對云層的光學厚度進行了分類，具體為：晴空、氣溶膠、薄云、可能為厚云、厚云、藍色的厚云(下文簡稱為“藍厚云”)。由于氣溶膠是影響DNI預測結(jié)果精度的主要因素之一，因此將其單獨作為一類進行分類。此外，部分云用肉眼很難判斷是薄云或厚云，因此添加了“可能為厚云” 這一類型。

在提取圖像特征之前對太陽及其鄰近區(qū)域進行了插值處理，這樣有助于得到較好的DNI預測結(jié)果。

在保證數(shù)據(jù)準確性的基礎上，本文采用半人工方式提取了一定數(shù)量的云層分類特征樣本，用于云層分類模型的訓練，并對云層分類模型進行了驗證。不同類型云層的分類結(jié)果占比如表1所示。

表1 不同類型云層的分類結(jié)果占比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預測模型訓練

根據(jù)之前計算得到的云的運動矢量，利用外推法獲得未來0～30 min內(nèi)(間隔1 min)云層可能會遮擋太陽的30個圖像點，然后對這30個圖像點進行圖像特征提取，并利用云層分類模型完成云層分類。以通過云層分類提取到的圖像特征和云層類型分類結(jié)果作為超短期DNI預測的特征。特征提取過程完全由軟件自動完成。

全天空成像儀旁邊安裝有DNI測量儀器，用于采集DNI數(shù)據(jù)。受云層和其他因素的影響，DNI會產(chǎn)生衰減，實測DNI的衰減率DDNI可表示為：

式中：DNIt為理論的DNI值；DNIm為實測的DNI值。

建立超短期DNI預測模型，并將DDNI作為其輸出值。該模型的訓練數(shù)據(jù)選用2018～2019年期間收集的所有多云天氣中以天為單位篩選出的云的形態(tài)較為穩(wěn)定時的數(shù)據(jù)，然后根據(jù)云層類型分類結(jié)果和實測DNI的衰減率對篩選出的數(shù)據(jù)進行進一步過濾，過濾掉云層類型和實測DNI的衰減率明顯不符的數(shù)據(jù)，使通過訓練得到的超短期DNI預測模型更為穩(wěn)定。最后利用XGB算法完成超短期DNI預測模型的訓練。

不同云層類型時DDNI的取值范圍如表3所示。

表3 不同云層類型時DDNI的取值范圍Table 3 Value range of DDNI for different cloud layers types

5.2 基于XGB算法的超短期DNI預測

超短期DNI預測時同樣根據(jù)云的運動矢量，采用外推法獲得未來0～30 min內(nèi)(間隔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預測結(jié)果進行統(tǒng)計分析，具體如表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預測結(jié)果的MAE與RMSETable 4 MAE and RMSE of DNI prediction results under different advance prediction time

(續(xù)表)

從表4可以看出，不同提前預測時間的總體DNI預測值的平均絕對誤差(MAE)為116.0 W/m2，均方根誤差(RMSE)為203.0 W/m2。在未使用可持續(xù)模型的情況下，從這2個值來看，預測結(jié)果的精度較高。

將實測的DNI數(shù)據(jù)與超短期DNI預測模型的預測結(jié)果進行對比，當以“某個時刻的DNI值是否大于或小于400 W/m2”作為判斷標準時，超短期DNI預測模型的預測結(jié)果的正確率均大于85%，符合塔式太陽能熱發(fā)電站運行過程中對DNI預測的指標要求。而通過進一步分析發(fā)現(xiàn)，影響DNI預測結(jié)果精度的主要因素為氣溶膠及云的不規(guī)則變化。

6 結(jié)論

本文未使用可持續(xù)模型，而是采用基于XGB算法及全天空成像儀圖像實現(xiàn)了超短期(0～30 min)的DNI預測。通過使用基于太陽位置的全天空成像儀標定方法，完成了全天空成像儀的姿態(tài)校正和畸變矯正；準確擬合了晴空背景參數(shù)并建立晴空庫，為云層的檢測提供了基礎；通過對太陽及其臨近區(qū)域進行插值修復，使1 min內(nèi)的DNI預測成為可能。根據(jù)云的運動矢量，利用外推法得到未來云層可能會遮擋太陽的圖像點，并提取圖像特征，用于基于XGB算法的超短期DNI預測模型的訓練，最終以訓練后的該模型完成超短期DNI預測，并得到了較好的結(jié)果，符合塔式太陽能熱發(fā)電站運行過程中對DNI預測的指標要求。