CN103353952A - 一种基于地基云图的光伏功率预测方法 - Google Patents

Abstract

本发明提供一种基于地基云图的光伏功率预测方法，包括以下步骤：步骤1：云团运动预估；步骤2：对地面辐照度进行预测；步骤3：建立综合光电转换模型；步骤4：进行光伏功率预测。本发明综合运用多种数字图像处理技术，突出问题需求，有很强的针对性和实用性，与传统的统计方法结合为光伏功率精细化预测提供了有效的技术支撑和方法参考。且本发明综合多因素多步建模，方法的描述已涵盖光伏预测领域所有关键环节，有较强的普适性和鲁棒性。本发明为了在使用中结合使用的具体场景和应用对象，有些环节预留了调节参数，操作方便并保证预测效果。

Description

一种基于地基云图的光伏功率预测方法

技术领域

本发明属于新能源光伏发电技术领域，具体涉及一种基于地基云图的光伏功率预测方法。

背景技术

近年来，大规模光伏电站接入电网，由于光伏发电输出功率具有随机性和波动性，对电网安全稳定和经济运行造成影响。对光伏电站的输出功率进行准确预测，为电力调度提供重要的决策支持，统筹安排常规电源和光伏发电的协调配合，不仅降低电力系统运行成本，而且有利于充分利用光伏资源，获得更大的经济效益和社会效益。

地面辐照度作为光伏电站输出功率的主要影响因素之一，它的不确定性直接导致输出功率的随机性和波动性。而云作为影响太阳地面辐照量的主要气象要素，其生消和移动变化是地面辐照度变化不确定性的主要原因。目前，光伏功率预测方法主要有统计分析或机器学习方法，它们基于历史气象要素数据和光伏电站输出功率数据进行预测，但均未涉及云图数据，对云的遮挡导致的功率突变很难学习导致预测效果欠佳。国外有学者根据气象卫星云图预估云层移动情况进行太阳辐照度预测，但其时空分辨率相对较低,为了满足区域内光伏电站超短期0-4小时内分钟级功率预测的要求，需要局地高时空分辨率的云图像数据。

近年来，随着电荷耦合器件CCD等硬件技术的发展和数字图像处理技术的不断完善，很多地基遥感测云仪器研制成功，比如全天空成像仪(total sky imager，TSI)，该仪器能实现对全天空云的自动化观测和实时采集，最小时间分辨率可达30秒，能很好地实现光伏电站上空云图的实时监测和采集。

发明内容

为了克服上述现有技术的不足，本发明提供一种基于地基云图的光伏功率预测方法，利用数字图像处理技术对地基云图进行处理和分析，对未来时刻云团运动和太阳遮挡状况进行预测，相当于将云这一主要气象要素从众多影响辐照度的随机因素中单独分离出来进行建模。同时结合衡量太阳通过大气层时衰减情况的辐射清晰度指数实现未来时刻地面辐照度的预测，最后通过光电转换模型实现光伏功率超短期1小时内分钟级光伏功率精细化预测。

为了实现上述发明目的，本发明采取如下技术方案：

提供一种基于地基云图的光伏功率预测方法，所述方法包括以下步骤：

步骤1：云团运动预估；

步骤2：对地面辐照度进行预测；

步骤3：建立综合光电转换模型；

步骤4：进行光伏功率预测。

所述步骤1包括以下步骤：

步骤1-1：确定图像上太阳位置；

所述步骤1‐1中，用天顶角和方位角来确定图像上太阳位置；

假设观测地的纬度为

，赤纬角为δ，时角为ω，天顶角和方位角分别用θ_z和γ_s表示为：

其中：

ω＝(t-12)×15°  （3）

其中，dayth表示一年中的日期序号，t为时间；

天顶角为和方位角分别为θ_z和γ_s时，距离图像中心点的距离为R，于是即可确定某一时刻图像上太阳位置；距离图像中心点的距离R表示为：

$R=-17.59{\mathrm{\θ}}_z^2+\mathrm{118.1.35}{\mathrm{\θ}}_z---\left(4\right).$

步骤1-2：对地基云图原始图像进行预处理；

所述步骤1-2包括以下步骤：

步骤1-2-1：对地基云图原始图像进行图像复原；

首先，通过遮光带覆盖像素点范围和镜头支臂的像素点范围确定遮光带和镜头支臂的掩模图像；根据已确定的图像上太阳位置，标定出复原图像对应时刻太阳在图像中的太阳位置点，然后由中心点到太阳位置点至到圆周边沿点确定一条直线，然后根据图像形态学膨胀运算扩展至遮光带宽度，从而确定遮光带覆盖像素点范围；同时由中心点垂直向上至到圆周边沿点确定另一条直线，然后根据图像形态学膨胀运算把该直线扩展至宽度为镜头支臂宽度，从而确定镜头支臂的像素点范围。

然后在地基云图原始图像上，与掩模图像对于遮光带覆盖范围的像素点和镜头支臂范围的像素点的周边像素点进行线性插值，最终生成复原图像I；

步骤1-2-2：对所述复原图像I进行中值滤波，得到中值滤波图像；

M(i,j)＝Med{I(i+r,j+s)}  （5）

其中，M(i,j)表示中值滤波图像在(i,j)点的灰度值，I(i+r,j+s)表示复原图像I在(i+r,j+s)点的灰度值；滤波窗口为5×5领域时，r,s∈{-2，-1，0，1，2}；滤波窗口为3×3领域时，r,s∈{-1，0，1}；Med{I(i+r,j+s)表示随着r,s取不同值而得到的25个序列或9个序列的中值；

步骤1-2-3：对中值滤波图像进行锐化；

对中值滤波图像进行锐化后的锐化图像假定为g(x,y)，则有：

$g(x,y)=M(x,y)-{\▿}^{2}M(x,y)---\left(6\right)$

其中，M(x,y)表示中值滤波图像，且有

${\▿}^{2}M(x,y)=\frac{{\∂}^{2}M(x,y)}{\∂x^2}+\frac{{\∂}^{2}M(x,y)}{\∂y^2}---\left(7\right).$

步骤1-3：对预处理后的锐化图像g(x,y)进行二值化处理，并提取云团质心；

所述步骤1-3包括以下步骤：

步骤1-3-1：对预处理最终得到的锐化图像g(x,y)进行二值化处理；

由下式得到(i,j)点的二值图像Binary：

$\mathrm{Binary}(i,j)=\left\{\begin{array}{cc}1;& R(i,j)/B(i,j)>\mathrm{\γ}\\ 0;& R(i,j)/B(i,j)\≤\mathrm{\γ}\end{array}\right.---\left(8\right)$

其中，Binary(i,j)表示二值图像Binary在(i,j)点的像素值，R(i,j)，B(i,j)分别表示锐化图像g(x,y)红通道分量和蓝通道分量在(i,j)点的灰度值；γ表示锐化图像g(x,y)红蓝通道分量比的阈值，由如下公式计算：

γ＝max_level(RBR_hist)+TH  （9）

其中，RBR_hist表示红蓝分量比图像RBR的灰度直方图，max_level表示灰度直方图中频度最大的灰度值，TH表示阈值微调量；

步骤1-3-2：进行云团检测和标注；

二值图像Binary中，若像素值为1且坐标为(i,j)的像素点为P，则搜索P的8连接点，8连接点的坐标集为：{(i-1,j-1)，(i-1,j)，(i-1,j+1)，(i,j-1)，(i,j+1)，(i+1,j-1)，(i+1,j)，(i+1,j+1)}；

对P点的8连接点的坐标集对应的像素点中选择像素值为1的点，并分别对其再次展开8连接点搜索，如此反复，并且每个像素点不能重复被搜索，直至所有连接点的8连接点均为已搜索的点或者像素值为0的点，则结束搜索并标注云团；

云团的标注按照检索的顺序，第一次检索到的团块标记为1，第二次检索到的团块标记为2，以此类推；

步骤1-3-3：云团质心提取；

设Al₁，Al₂为顶点的矩形框经过云团A1最小外接矩形的左上角和右下角像素点，若Al₁点坐标为(Alx₁,Aly₁)，Al₂点坐标为(Alx₂,Aly₂)，云团A1的质心坐标为(Alx,Aly)，则有：

$\left\{\begin{array}{c}\mathrm{Alx}=\frac{\mathrm{Al}{x}_1+\mathrm{Al}{x}_2}{2}\\ \mathrm{Aly}=\frac{\mathrm{Al}{y}_1+\mathrm{Al}{y}_2}{2}\end{array}\right.---\left(10\right)$

按照上述云团A1质心提取的方法，依次提取云团A1,A2,…和Am共m个云团的质心。

步骤1-4：云团运动预估：

所述步骤1-4包括以下步骤：

步骤1-4-1：定义云团运动方向和速度；

建立以球面图像中心为坐标原点的极坐标系，球面图像中心的坐标为(x_c,y_c)，像素点P的原始坐标为(x,y)，其中x≠x_c；则像素点P在极坐标下的极径ρ和极角θ分别表示为：

$\left\{\begin{array}{c}\mathrm{\ρ}=\sqrt{{(x-x_c)}^2+{(y-y_c)}^2}\\ \mathrm{\θ}=\arctan (-\frac{y-y_c}{x-x_c})\end{array}\right.---\left(11\right)$

(x_t,y_t)和(x_t+1,y_t+1)分别为t时刻和t+1时刻某云团的质心坐标；则t时刻该云团移动速度v表示为：

$v=\sqrt{{(x_{t+1}-x_t)}^2+{(y_{t+1}-y_t)}^2}/\mathrm{\Δt}---\left(12\right)$

其中Δt为t时刻和t+1时刻之间的采样时间间隔；

云团移动方向定义如下：

（1）如果x_t+1＞x_t，且y_t+1≤y_t，则有

（2）如果x_t+1＝x_t，且y_t+1≤y_t，则有θ=π/2；

（3）如果x_t+1＜x_t，且y_t+1≤y_t，则有

（4）如果x_t+1＜x_t，且y_t+1≥y_t，则有

（5）如果x_t+1＝x_t，且y_t+1≥y_t，则有θ=3π/2；

（6）如果xt+1＞xt，且yt+1≥yt，则有

步骤1-4-2：进行云团运动预估；

假设二值图像Binary有n张，每张二值图像Binary中云团个数有m个，第j个云团根据第i张二值图像Binary中的质心位置和第i+1张二值图像Binary中的质心位置所得到的速度和角度分别为v_ji和θ_ji，其中：j＝1,2,...,m；i＝1,2,...,n-1；

假设云团在相对较短时间内的运动为刚体运动，位于全天空成像仪视野范围内的云团近似匀速直线运动，通过求取v_ji和θ_ji平均值预估未来第j个云团的移动速度v_j和移动方向θ_j，具体表示为：

$\left\{\begin{array}{c}{v}_j=\frac{\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}{v}_{\mathrm{ji}}}{n-1}\\ {\mathrm{\θ}}_j=\frac{\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{\θ}}_{\mathrm{ji}}}{n-1}\end{array}\right.---\left(13\right).$

所述步骤2包括以下步骤：

步骤2-1：计算学习阶段太阳被云团遮挡时刻对应的辐射清晰度指数；

t时刻天顶辐射H₀(t)和大气质量AM(t)分别表示为：

H₀(t)＝I_sc(1+0.033cos(2π·dayth/365))cosθ_Z  （14）

$\mathrm{AM}\left(t\right)=\sqrt{1229+{\left(614\cos {\mathrm{\θ}}_Z\right)}^2-614\cos {\mathrm{\θ}}_Z}{---\left(15\right)}_{}$

其中，I_sc＝1367±5W/m²；

于是：t时刻地面理论辐射值I_theory(t)表示为：

I_theory(t)＝H₀(t)/AM^C(t)  （16）

在晴空状况下，令t时刻地面理论辐射值I_theory(t)等于地面实测辐射值I_obs(t)根据式（16）计算出C，进而确定t时刻地面理论辐射值I_theory(t)；

学习阶段太阳被云团遮挡时刻对应的辐射清晰度指数用来衡量太阳通过大气层时云团对辐照度的衰减情况，其表示为：

$\mathrm{KT}\left(t\right)=\frac{I_{\mathrm{obs}}\left(t\right)}{I_{\mathrm{theory}\left(t\right)}}---\left(17\right)$

其中，KT(t)为学习阶段太阳被云团遮挡时刻对应的辐射清晰度指数；

步骤2-2：进行太阳遮挡预测；

通过判断二值图像Binary上太阳位置处像素点是否为二值图像Binary中云团的像素点，进而判断此刻太阳是否被云团遮挡，具体表示为：

其中sign(t)表示t时刻太阳的遮挡信息，sign(t)=0表示太阳未被云团遮挡，sign(t)=1表示太阳被云团遮挡；(x_S(t),y_S(t))表示t时刻太阳在图像中的位置，ex的取值为{-2，-1，0，1，2}；

步骤2-3：对地面辐照度进行预测；

假设预测阶段的辐射清晰度指数为pKT，其由学习阶段太阳被云团遮挡时刻对应的所有辐射清晰度指数的加权平均确定，表示为：

$\mathrm{pKT}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\λ}}_i\mathrm{KT}\left(t_i\right)---\left(19\right)$

其中，N为学习阶段太阳被云团遮挡的时刻；i表示太阳被云团遮挡时刻时的排序，i＝1,2,...,N-1；λ_i为学习样本中第i个辐射清晰度指数加权平均的权重系数，且满足λ_i＜λ_i+1；KT(t_i)为学习阶段太阳被云团遮挡时刻对应的第i个辐射清晰度指数；

于是，预测阶段t时刻的地面理论辐照度表示为：

$\mathrm{pI}\left(t\right)=\left\{\begin{array}{cc}\mathrm{pKT}\·I_{\mathrm{theory}}\left(t\right)& \mathrm{sign}\left(t\right)=1\\ I_{\mathrm{theory}\left(t\right)}& \mathrm{sign}\left(t\right)=0\end{array}\right.---\left(20\right)$

其中，pI(t)表示预测阶段t时刻的地面理论辐照度。

所述步骤3包括以下步骤：

步骤3-1：建立基础光电转换模型；

通过所读取的历史功率数据和辐射数据绘制散点图，若P₀(t)、I(t)分别为t时刻的实测有功功率和实测地面辐照度，并根据散点图确定曲线类型和曲线拟合方程式中的阶数q，所述曲线拟合方程式为：

P₀(t)＝a_qI^q(t)+a_q-1I^q-1(t)+…+a₀  （21）

其中a_q,a_q-1,…和a₀均为曲线拟合方程式的系数；

根据上述曲线拟合方程式建立目标函数：

$\underset{a_1\·\·\·a_q}{\min }\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(a_q{I}^q\left(t\right)+a_{q-1}{I}^{q-1}\left(t\right)+\·\·\·+a_0-P_0\left(t\right))}^2---\left(22\right)$

通过最小二乘法确定a_q,a_q-1,…和a₀，从而确定t时刻的实测有功功率P₀(t)，表示为：

P₀(t)＝f(I(t))＝a_qI^q(t)+a_q-1I^q-1(t)+…+a₀  （23）

f(I(t))为t时刻的实测有功功率P₀(t)关于t时刻的实测地面辐照度I(t)的函数；

步骤3-2：结合基础光电转换模型建立综合光电转换模型；

设t时刻的组建温度分别为T(t)，由温度导致的t时刻的转换效率η(t)表示为：

η(t)＝1+0.0035(T(t)-T(t-1))  （24）

由此可得预测阶段t时刻的实测有功功率，具体表示为：

P(t)＝η(t)f(I(t))  （25）

其中，P(t)表示预测阶段t时刻的实测有功功率。

所述步骤4中，通过建立光伏功率预测模型进行光伏功率预测；

所述光伏功率预测模型为：

$P\left(t\right)=\left\{\begin{array}{cc}\mathrm{\η}\left(t\right)f(\mathrm{pKT}\·I_{\mathrm{theory}}\left(t\right))& \mathrm{sign}\left(t\right)=1\\ \mathrm{\η}\left(t\right)f\left(I_{\mathrm{theory}}\left(t\right)\right)& \mathrm{sign}\left(t\right)=0\end{array}\right.---\left(26\right).$

与现有技术相比，本发明的有益效果在于：

（1）本发明综合运用多种数字图像处理技术，突出问题需求，有很强的针对性，与传统的统计方法结合为光伏功率精细化预测提供了有效的技术支撑和方法参考。

（2）本发明综合多因素多步建模，方法的描述已涵盖光伏预测领域所有关键环节，有较强的普适性和鲁棒性。

（2）本发明采用必要且合理的假设，实现对问题简化的同时又能最大限度地准确描述问题。

（3）本发明为了在使用中结合使用的具体场景和应用对象，有些环节预留了调节参数，操作方便并保证预测效果。

（4）流程化的处理，详尽的算法描述和形式化表达，为工程实现提供了方便，具有很好的工程实用性。

附图说明

图1是本发明实施例中某一时刻图像上太阳的位置示意图；

图2是图像形态学膨胀运算原理示意图；

图3是本发明实施例中掩模图像示意图；

图4是本发明实施例中对地基云图原始图像进行图像复原示意图；

图5是本发明实施例中8连接点的连接方式示意图；

图6是本发明实施例中云团质心提取示意图；

图7是本发明实施例中建立以球面图像中心为坐标原点的极坐标系中云团运动预估示意图。

具体实施方式

下面结合附图对本发明作进一步详细说明。

提供一种基于地基云图的光伏功率预测方法，所述方法包括以下步骤：

步骤1：云团运动预估；

步骤2：对地面辐照度进行预测；

步骤3：建立综合光电转换模型；

步骤4：进行光伏功率预测。

所述步骤1包括以下步骤：

步骤1-1：确定图像上太阳位置；

所述步骤1‐1中，用天顶角和方位角来确定图像上太阳位置；

假设观测地的纬度为

赤纬角为δ，时角为ω，天顶角和方位角分别用θ_z和γ_s表示为：

其中：

ω＝(t-12)×15°  （3）

其中，dayth表示一年中的日期序号，t为时间；

天顶角为和方位角分别为θ_z和γ_s时，距离图像中心点的距离为R，于是即可确定某一时刻图像上太阳位置；由于天顶角在与图像半径有相关性，根据式（4）可将距离图像中心点的距离R表示为：

$R=-17.59{\mathrm{\θ}}_z^2+\mathrm{118.1.35}{\mathrm{\θ}}_z---\left(4\right).$

如图1，Center点为图像中心点，这样就可以标定出某一时刻图像上太阳的位置。

步骤1-2：对地基云图原始图像进行预处理；

所述步骤1-2包括以下步骤：

步骤1-2-1：对地基云图原始图像进行图像复原；

首先，通过遮光带覆盖像素点范围和镜头支臂的像素点范围确定遮光带和镜头支臂的掩模图像；根据已确定的图像上太阳位置，标定出复原图像对应时刻太阳在图像中的太阳位置点，然后由中心点到太阳位置点至到圆周边沿点确定一条直线，然后根据图像形态学膨胀运算扩展至遮光带宽度，从而确定遮光带覆盖像素点范围；同时由中心点垂直向上至到圆周边沿点确定另一条直线，然后根据图像形态学膨胀运算把该直线扩展至宽度为镜头支臂宽度，从而确定镜头支臂的像素点范围。

如图2，图像形态学膨胀运算相当于用B的中心在A的边界走一圈，B的边界扫过的区域加上A本身的区域就是A膨胀B的结果。

然后在地基云图原始图像上，与掩模图像对于遮光带覆盖范围的像素点和镜头支臂范围的像素点的周边像素点进行线性插值，最终生成复原图像I；

步骤1-2-2：对所述复原图像I进行中值滤波，得到中值滤波图像；

M(i,j)＝Med{I(i+r,j+s)}  （5）

其中，M(i,j)表示中值滤波图像在(i,j)点的灰度值，I(i+r,j+s)表示复原图像I在(i+r,j+s)点的灰度值；滤波窗口为5×5领域时，r,s∈{-2，-1，0，1，2}；滤波窗口为3×3领域时，r,s∈{-1，0，1}；Med{I(i+r,j+s)表示随着r,s取不同值而得到的25个序列或9个序列的中值；

步骤1-2-3：对中值滤波图像进行锐化；

对中值滤波图像进行锐化后的锐化图像假定为g(x,y)，则有：

$g(x,y)=M(x,y)-{\▿}^{2}M(x,y)---\left(6\right)$

其中，M(x,y)表示中值滤波图像，且有

${\▿}^{2}M(x,y)=\frac{{\∂}^{2}M(x,y)}{\∂x^2}+\frac{{\∂}^{2}M(x,y)}{\∂y^2}---\left(7\right).$

步骤1-3：对预处理后的锐化图像g(x,y)进行二值化处理，并提取云团质心；

所述步骤1-3包括以下步骤：

步骤1-3-1：对预处理最终得到的锐化图像g(x,y)进行二值化处理；

由下式得到(i,j)点的二值图像Binary：

$\mathrm{Binary}(i,j)=\left\{\begin{array}{cc}1;& R(i,j)/B(i,j)>\mathrm{\γ}\\ 0;& R(i,j)/B(i,j)\≤\mathrm{\γ}\end{array}\right.---\left(8\right)$

其中，Binary(i,j)表示二值图像Binary在(i,j)点的像素值，R(i,j)，B(i,j)分别表示锐化图像g(x,y)红通道分量和蓝通道分量在(i,j)点的灰度值；γ表示锐化图像g(x,y)红蓝通道分量比的阈值，由如下公式计算：

γ＝max_level(RBR_hist)+TH  （9）

其中，RBR_hist表示红蓝分量比图像RBR的灰度直方图，max_level表示灰度直方图中频度最大的灰度值，TH表示阈值微调量；

步骤1-3-2：进行云团检测和标注；

二值图像Binary中，若像素值为1且坐标为(i,j)的像素点为P，则搜索P的8连接点，其中8连接点连接方式如图5所示；

8连接点的坐标集为：{(i-1,j-1)，(i-1,j)，(i-1,j+1)，(i,j-1)，(i,j+1)，(i+1,j-1)，(i+1,j)，(i+1,j+1)}；

对P点的8连接点的坐标集对应的像素点中选择像素值为1的点，并分别对其再次展开8连接点搜索，如此反复，并且每个像素点不能重复被搜索，直至所有连接点的8连接点均为已搜索的点或者像素值为0的点，则结束搜索并标注云团；

云团的标注按照检索的顺序，第一次检索到的团块标记为1，第二次检索到的团块标记为2，以此类推；

步骤1-3-3：云团质心提取；

如图6，设Al₁，Al₂为顶点的矩形框经过云团A1最小外接矩形的左上角和右下角像素点，若Al₁点坐标为(Alx₁,Aly₁)，Al₂点坐标为(Alx₂,Aly₂)，云团A1的质心坐标为(Alx,Aly)，则有：

$\left\{\begin{array}{c}\mathrm{Alx}=\frac{\mathrm{Al}{x}_1+\mathrm{Al}{x}_2}{2}\\ \mathrm{Aly}=\frac{\mathrm{Al}{y}_1+\mathrm{Al}{y}_2}{2}\end{array}\right.---\left(10\right)$

按照上述云团A1质心提取的方法，依次提取云团A1,A2,…和Am共m个云团的质心。

步骤1-4：云团运动预估：

所述步骤1-4包括以下步骤：

步骤1-4-1：定义云团运动方向和速度；

如图7，建立以球面图像中心为坐标原点的极坐标系，球面图像中心的坐标为(x_c,y_c)，像素点P的原始坐标为(x,y)，其中x≠x_c；则像素点P在极坐标下的极径ρ和极角θ分别表示为：

$\left\{\begin{array}{c}\mathrm{\ρ}=\sqrt{{(x-x_c)}^2+{(y-y_c)}^2}\\ \mathrm{\θ}=\arctan (-\frac{y-y_c}{x-x_c})\end{array}\right.---\left(10\right)$

(x_t,y_t)和(x_t+1,y_t+1)分别为t时刻和t+1时刻某云团的质心坐标；则t时刻该云团移动速度v表示为：

$v=\sqrt{{(x_{t+1}-x_t)}^2+{(y_{t+1}-y_t)}^2}/\mathrm{\Δt}---\left(12\right)$

其中Δt为t时刻和t+1时刻之间的采样时间间隔；

云团移动方向定义如下：

（1）如果x_t+1＞x_t，且y_t+1≤y_t，则有

（2）如果x_t+1＝x_t，且y_t+1≤y_t，则有θ=π/2；

（3）如果x_t+1＜x_t，且y_t+1≤y_t，则有

（4）如果x_t+1＜x_t，且y_t+1≥y_t，则有

（5）如果x_t+1＝x_t，且y_t+1≥y_t，则有θ=3π/2；

（6）如果x_t+1＞x_t，且y_t+1≥y_t，则有

步骤1-4-2：进行云团运动预估；

假设二值图像Binary有n张，每张二值图像Binary中云团个数有m个，第j个云团根据第i张二值图像Binary中的质心位置和第i+1张二值图像Binary中的质心位置所得到的速度和角度分别为v_ji和θ_ji，其中：j＝1,2,...,m；i＝1,2,...,n-1；

假设云团在相对较短时间内的运动为刚体运动，位于全天空成像仪视野范围内的云团近似匀速直线运动，通过求取v_ji和θ_ji平均值预估未来第j个云团的移动速度v_j和移动方向θ_j，具体表示为：

$\left\{\begin{array}{c}{v}_j=\frac{\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}{v}_{\mathrm{ji}}}{n-1}\\ {\mathrm{\θ}}_j=\frac{\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{\θ}}_{\mathrm{ji}}}{n-1}\end{array}\right.---\left(13\right).$

所述步骤2包括以下步骤：

步骤2-1：计算学习阶段太阳被云团遮挡时刻对应的辐射清晰度指数；

t时刻天顶辐射H₀(t)和大气质量AM(t)分别表示为：

H₀(t)＝I_sc(1+0.033cos(2π·dayth/365))cosθ_Z  （14）

$\mathrm{AM}\left(t\right)=\sqrt{1229+{\left(614\cos {\mathrm{\θ}}_Z\right)}^2-614\cos {\mathrm{\θ}}_Z}{---\left(15\right)}_{}$

其中，I_sc＝1367±5W/m²；

于是：t时刻地面理论辐射值I_theory(t)表示为：

I_theory(t)＝H₀(t)/AM^C(t)  （16）

在晴空状况下，令t时刻地面理论辐射值I_theory(t)等于地面实测辐射值I_obs(t)根据式（16）计算出C，进而确定t时刻地面理论辐射值I_theory(t)；

学习阶段太阳被云团遮挡时刻对应的辐射清晰度指数用来衡量太阳通过大气层时云团对辐照度的衰减情况，其表示为：

$\mathrm{KT}\left(t\right)=\frac{I_{\mathrm{obs}}\left(t\right)}{I_{\mathrm{theory}\left(t\right)}}---\left(17\right)$

其中，KT(t)为学习阶段太阳被云团遮挡时刻对应的辐射清晰度指数；

步骤2-2：进行太阳遮挡预测；

通过判断二值图像Binary上太阳位置处像素点是否为二值图像Binary中云团的像素点，进而判断此刻太阳是否被云团遮挡，具体表示为：

其中sign(t)表示t时刻太阳的遮挡信息，sign(t)=0表示太阳未被云团遮挡，sign(t)=1表示太阳被云团遮挡；(x_S(t),y_S(t))表示t时刻太阳在图像中的位置，ex的取值为{-2，-1，0，1，2}；

步骤2-3：对地面辐照度进行预测；

假设预测阶段的辐射清晰度指数为pKT，其由学习阶段太阳被云团遮挡时刻对应的所有辐射清晰度指数的加权平均确定，表示为：

$\mathrm{pKT}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\λ}}_i\mathrm{KT}\left(t_i\right)---\left(19\right)$

其中，N为学习阶段太阳被云团遮挡的时刻；i表示太阳被云团遮挡时刻时的排序，i＝1,2,...,N-1；λ_i为学习样本中第i个辐射清晰度指数加权平均的权重系数，且满足λ_i＜λ_i+1；KT(t_i)为学习阶段太阳被云团遮挡时刻对应的第i个辐射清晰度指数；

于是，预测阶段t时刻的地面理论辐照度表示为：

$\mathrm{pI}\left(t\right)=\left\{\begin{array}{cc}\mathrm{pKT}\·I_{\mathrm{theory}}\left(t\right)& \mathrm{sign}\left(t\right)=1\\ I_{\mathrm{theory}\left(t\right)}& \mathrm{sign}\left(t\right)=0\end{array}\right.---\left(20\right)$

其中，pI(t)表示预测阶段t时刻的地面理论辐照度。

所述步骤3包括以下步骤：

步骤3-1：建立基础光电转换模型；

通过所读取的历史功率数据和辐射数据绘制散点图，若P₀(t)、I(t)分别为t时刻的实测有功功率和实测地面辐照度，并根据散点图确定曲线类型和曲线拟合方程式中的阶数q，所述曲线拟合方程式为：

P₀(t)＝a_qI^q(t)+a_q-1I^q-1(t)+…+a₀  （21）

其中a_q,a_q-1,…和a₀均为曲线拟合方程式的系数；

根据上述曲线拟合方程式建立目标函数：

$\underset{a_1\·\·\·a_q}{\min }\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(a_q{I}^q\left(t\right)+a_{q-1}{I}^{q-1}\left(t\right)+\·\·\·+a_0-P_0\left(t\right))}^2---\left(22\right)$

通过最小二乘法确定a_q,a_q-1,…和a₀，从而确定t时刻的实测有功功率P₀(t)，表示为：

P₀(t)＝f(I(t))＝a_qI^q(t)+a_q-1I^q-1(t)+…+a₀  （23）

f(I(t))为t时刻的实测有功功率P₀(t)关于t时刻的实测地面辐照度I(t)的函数；

步骤3-2：结合基础光电转换模型建立综合光电转换模型；

设t时刻的组建温度分别为T(t)，由温度导致的t时刻的转换效率η(t)表示为：

η(t)＝1+0.0035(T(t)-T(t-1))  （24）

由此可得预测阶段t时刻的实测有功功率，具体表示为：

P(t)＝η(t)f(I(t))  （25）

其中，P(t)表示预测阶段t时刻的实测有功功率。

所述步骤4中，通过建立光伏功率预测模型进行光伏功率预测；

所述光伏功率预测模型为：

$P\left(t\right)=\left\{\begin{array}{cc}\mathrm{\η}\left(t\right)f(\mathrm{pKT}\·I_{\mathrm{theory}}\left(t\right))& \mathrm{sign}\left(t\right)=1\\ \mathrm{\η}\left(t\right)f\left(I_{\mathrm{theory}}\left(t\right)\right)& \mathrm{sign}\left(t\right)=0\end{array}\right.---\left(26\right).$

最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求范围当中。

Claims (10)

1.一种基于地基云图的光伏功率预测方法，其特征在于：所述方法包括以下步骤：

步骤1：云团运动预估；

步骤2：对地面辐照度进行预测；

步骤3：建立综合光电转换模型；

步骤4：进行光伏功率预测。

2.根据权利要求1所述的基于地基云图的光伏功率预测方法，其特征在于：所述步骤1包括以下步骤：

步骤1-1：确定图像上太阳位置；

步骤1-2：对地基云图原始图像进行预处理；

步骤1-3：对预处理后的锐化图像g(x,y)进行二值化处理，并提取云团质心；

步骤1-4：云团运动预估。

3.根据权利要求2所述的基于地基云图的光伏功率预测方法，其特征在于：所述步骤1-1中，用天顶角和方位角来确定图像上太阳位置；

假设观测地的纬度为，赤纬角为δ，时角为ω，天顶角和方位角分别用θ_z和γ_s表示为：

其中：

ω＝(t-12)×15°  （3）

其中，dayth表示一年中的日期序号，t为时间；

天顶角为和方位角分别为θ_z和γ_s时，距离图像中心点的距离为R，于是即可确定某一时刻图像上太阳位置；距离图像中心点的距离R表示为：

$R=-17.59{\mathrm{\θ}}_z^2+\mathrm{118.1.35}{\mathrm{\θ}}_z---\left(4\right).$

4.根据权利要求2所述的基于地基云图的光伏功率预测方法，其特征在于：所述步骤1-2包括以下步骤：

步骤1-2-1：对地基云图原始图像进行图像复原；

首先，通过遮光带覆盖像素点范围和镜头支臂的像素点范围确定遮光带和镜头支臂的掩模图像；然后在地基云图原始图像上，与掩模图像对于遮光带覆盖范围的像素点和镜头支臂范围的像素点的周边像素点进行线性插值，最终生成复原图像I；

步骤1-2-2：对所述复原图像I进行中值滤波，得到中值滤波图像；

M(i,j)＝Med{I(i+r,j+s)}  （5）

其中，M(i,j)表示中值滤波图像在(i,j)点的灰度值，I(i+r,j+s)表示复原图像I在(i+r,j+s)点的灰度值；滤波窗口为5×5领域时，r,s∈{-2，-1，0，1，2}；滤波窗口为3×3领域时，r,s∈{-1，0，1}；Med{I(i+r,j+s)表示随着r,s取不同值而得到的25个序列或9个序列的中值；

步骤1-2-3：对中值滤波图像进行锐化；

对中值滤波图像进行锐化后的锐化图像假定为g(x,y)，则有：

$g(x,y)=M(x,y)-{\▿}^{2}M(x,y)---\left(6\right)$

其中，M(x,y)表示中值滤波图像，且有

${\▿}^{2}M(x,y)=\frac{{\∂}^{2}M(x,y)}{\∂x^2}+\frac{{\∂}^{2}M(x,y)}{\∂y^2}---\left(7\right).$

5.根据权利要求4所述的基于地基云图的光伏功率预测方法，其特征在于：根据已确定的图像上太阳位置，标定出复原图像对应时刻太阳在图像中的太阳位置点，然后由中心点到太阳位置点至到圆周边沿点确定一条直线，然后根据图像形态学膨胀运算扩展至遮光带宽度，从而确定遮光带覆盖像素点范围；

同时由中心点垂直向上至到圆周边沿点确定另一条直线，然后根据图像形态学膨胀运算把该直线扩展至宽度为镜头支臂宽度，从而确定镜头支臂的像素点范围。

6.根据权利要求2所述的基于地基云图的光伏功率预测方法，其特征在于：所述步骤1-3包括以下步骤：

步骤1-3-1：对预处理最终得到的锐化图像g(x,y)进行二值化处理；

由下式得到(i,j)点的二值图像Binary：

$\mathrm{Binary}(i,j)=\left\{\begin{array}{cc}1;& R(i,j)/B(i,j)>\mathrm{\γ}\\ 0;& R(i,j)/B(i,j)\≤\mathrm{\γ}\end{array}\right.---\left(8\right)$

其中，Binary(i,j)表示二值图像Binary在(i,j)点的像素值，R(i,j)，B(i,j)分别表示锐化图像g(x,y)红通道分量和蓝通道分量在(i,j)点的灰度值；γ表示锐化图像g(x,y)红蓝通道分量比的阈值，由如下公式计算：

γ＝max_level(RBR_hist)+TH  （9）

其中，RBR_hist表示红蓝分量比图像RBR的灰度直方图，max_level表示灰度直方图中频度最大的灰度值，TH表示阈值微调量；

步骤1-3-2：进行云团检测和标注；

二值图像Binary中，若像素值为1且坐标为(i,j)的像素点为P，则搜索P的8连接点，8连接点的坐标集为：{(i-1,j-1)，(i-1,j)，(i-1,j+1)，(i,j-1)，(i,j+1)，(i+1,j-1)，(i+1,j)，(i+1,j+1)}；

对P点的8连接点的坐标集对应的像素点中选择像素值为1的点，并分别对其再次展开8连接点搜索，如此反复，并且每个像素点不能重复被搜索，直至所有连接点的8连接点均为已搜索的点或者像素值为0的点，则结束搜索并标注云团；

云团的标注按照检索的顺序，第一次检索到的团块标记为1，第二次检索到的团块标记为2，以此类推；

步骤1-3-3：云团质心提取；

设Al₁，Al₂为顶点的矩形框经过云团A1最小外接矩形的左上角和右下角像素点，若Al₁点坐标为(Alx₁,Aly₁)，Al₂点坐标为(Alx₂,Aly₂)，云团A1的质心坐标为(Alx,Aly)，则有：

$\left\{\begin{array}{c}\mathrm{Alx}=\frac{\mathrm{Al}{x}_1+\mathrm{Al}{x}_2}{2}\\ \mathrm{Aly}=\frac{\mathrm{Al}{y}_1+\mathrm{Al}{y}_2}{2}\end{array}\right.---\left(10\right)$

按照上述云团A1质心提取的方法，依次提取云团A1,A2,…和Am共m个云团的质心。

7.根据权利要求2所述的基于地基云图的光伏功率预测方法，其特征在于：所述步骤1-4包括以下步骤：

步骤1-4-1：定义云团运动方向和速度；

建立以球面图像中心为坐标原点的极坐标系，球面图像中心的坐标为(x_c,y_c)，像素点P的原始坐标为(x,y)，其中x≠x_c；则像素点P在极坐标下的极径ρ和极角θ分别表示为：

$\left\{\begin{array}{c}\mathrm{\ρ}=\sqrt{{(x-x_c)}^2+{(y-y_c)}^2}\\ \mathrm{\θ}=\arctan (-\frac{y-y_c}{x-x_c})\end{array}\right.---\left(11\right)$

(x_t,y_t)和(x_t+1,y_t+1)分别为t时刻和t+1时刻某云团的质心坐标；则t时刻该云团移动速度v表示为：

$v=\sqrt{{(x_{t+1}-x_t)}^2+{(y_{t+1}-y_t)}^2}/\mathrm{\Δt}---\left(12\right)$

其中Δt为t时刻和t+1时刻之间的采样时间间隔；

云团移动方向定义如下：

（1）如果x_t+1＞x_t，且y_t+1≤y_t，则有

（2）如果x_t+1＝x_t，且y_t+1≤y_t，则有θ=π/2；

（3）如果x_t+1＜x_t，且y_t+1≤y_t，则有

（4）如果x_t+1＜x_t，且y_t+1≥y_t，则有

（5）如果x_t+1＝x_t，且y_t+1≥y_t，则有θ=3π/2；

（6）如果x_t+1＞x_t，且y_t+1≥y_t，则有

步骤1-4-2：进行云团运动预估；

假设二值图像Binary有n张，每张二值图像Binary中云团个数有m个，第j个云团根据第i张二值图像Binary中的质心位置和第i+1张二值图像Binary中的质心位置所得到的速度和角度分别为v_ji和θ_ji，其中：j＝1,2,...,m；i＝1,2,...,n-1；

假设云团在相对较短时间内的运动为刚体运动，位于全天空成像仪视野范围内的云团近似匀速直线运动，通过求取v_ji和θ_ji平均值预估未来第j个云团的移动速度v_j和移动方向θ_j，具体表示为：

$\left\{\begin{array}{c}{v}_j=\frac{\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}{v}_{\mathrm{ji}}}{n-1}\\ {\mathrm{\θ}}_j=\frac{\underset{i=1}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{\θ}}_{\mathrm{ji}}}{n-1}\end{array}\right.---\left(13\right).$

8.根据权利要求1所述的基于地基云图的光伏功率预测方法，其特征在于：所述步骤2包括以下步骤：

步骤2-1：计算学习阶段太阳被云团遮挡时刻对应的辐射清晰度指数；

t时刻天顶辐射H₀(t)和大气质量AM(t)分别表示为：

H₀(t)＝I_sc(1+0.033cos(2π·dayth/365))cosθ_Z  （14）

$\mathrm{AM}\left(t\right)=\sqrt{1229+{\left(614\cos {\mathrm{\θ}}_Z\right)}^2-614\cos {\mathrm{\θ}}_Z}{---\left(15\right)}_{}$

其中，I_sc＝1367±5W/m²；

于是：t时刻地面理论辐射值I_theory(t)表示为：

I_theory(t)＝H₀(t)/AM^C(t)  （16）

在晴空状况下，令t时刻地面理论辐射值I_theory(t)等于地面实测辐射值I_obs(t)根据式（16）计算出C，进而确定t时刻地面理论辐射值I_theory(t)；

学习阶段太阳被云团遮挡时刻对应的辐射清晰度指数用来衡量太阳通过大气层时云团对辐照度的衰减情况，其表示为：

$\mathrm{KT}\left(t\right)=\frac{I_{\mathrm{obs}}\left(t\right)}{I_{\mathrm{theory}\left(t\right)}}---\left(17\right)$

其中，KT(t)为学习阶段太阳被云团遮挡时刻对应的辐射清晰度指数；

步骤2-2：进行太阳遮挡预测；

通过判断二值图像Binary上太阳位置处像素点是否为二值图像Binary中云团的像素点，进而判断此刻太阳是否被云团遮挡，具体表示为：

其中sign(t)表示t时刻太阳的遮挡信息，sign(t)=0表示太阳未被云团遮挡，sign(t)=1表示太阳被云团遮挡；(x_S(t),y_S(t))表示t时刻太阳在图像中的位置，ex的取值为{-2，-1，0，1，2}；

步骤2-3：对地面辐照度进行预测；

假设预测阶段的辐射清晰度指数为pKT，其由学习阶段太阳被云团遮挡时刻对应的所有辐射清晰度指数的加权平均确定，表示为：

$\mathrm{pKT}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\λ}}_i\mathrm{KT}\left(t_i\right)---\left(19\right)$

其中，N为学习阶段太阳被云团遮挡的时刻；i表示太阳被云团遮挡时刻时的排序，i＝1,2,...,N-1；λ_i为学习样本中第i个辐射清晰度指数加权平均的权重系数，且满足λ_i＜λ_i+1；KT(t_i)为学习阶段太阳被云团遮挡时刻对应的第i个辐射清晰度指数；

于是，预测阶段t时刻的地面理论辐照度表示为：

$\mathrm{pI}\left(t\right)=\left\{\begin{array}{cc}\mathrm{pKT}\·I_{\mathrm{theory}}\left(t\right)& \mathrm{sign}\left(t\right)=1\\ I_{\mathrm{theory}\left(t\right)}& \mathrm{sign}\left(t\right)=0\end{array}\right.---\left(20\right)$

其中，pI(t)表示预测阶段t时刻的地面理论辐照度。

9.根据权利要求1所述的基于地基云图的光伏功率预测方法，其特征在于：所述步骤3包括以下步骤：

步骤3-1：建立基础光电转换模型；

通过所读取的历史功率数据和辐射数据绘制散点图，若P₀(t)、I(t)分别为t时刻的实测有功功率和实测地面辐照度，并根据散点图确定曲线类型和曲线拟合方程式中的阶数q，所述曲线拟合方程式为：

P₀(t)＝a_qI^q(t)+a_q-1I^q-1(t)+…+a₀  （21）

其中a_q,a_q-1,…和a₀均为曲线拟合方程式的系数；

根据上述曲线拟合方程式建立目标函数：

$\underset{a_1\·\·\·a_q}{\min }\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(a_q{I}^q\left(t\right)+a_{q-1}{I}^{q-1}\left(t\right)+\·\·\·+a_0-P_0\left(t\right))}^2---\left(22\right)$

通过最小二乘法确定a_q,a_q-1,…和a₀，从而确定t时刻的实测有功功率P₀(t)，表示为：

P₀(t)＝f(I(t))＝a_qI^q(t)+a_q-1I^q-1(t)+…+a₀  （23）

f(I(t))为t时刻的实测有功功率P₀(t)关于t时刻的实测地面辐照度I(t)的函数；

步骤3-2：结合基础光电转换模型建立综合光电转换模型；

设t时刻的组建温度分别为T(t)，由温度导致的t时刻的转换效率η(t)表示为：

η(t)＝1+0.0035(T(t)-T(t-1))  （24）

由此可得预测阶段t时刻的实测有功功率，具体表示为：

P(t)＝η(t)f(I(t))  （25）

其中，P(t)表示预测阶段t时刻的实测有功功率。

10.根据权利要求1所述的基于地基云图的光伏功率预测方法，其特征在于：所述步骤4中，通过建立光伏功率预测模型进行光伏功率预测；

所述光伏功率预测模型为：

$P\left(t\right)=\left\{\begin{array}{cc}\mathrm{\η}\left(t\right)f(\mathrm{pKT}\·I_{\mathrm{theory}}\left(t\right))& \mathrm{sign}\left(t\right)=1\\ \mathrm{\η}\left(t\right)f\left(I_{\mathrm{theory}}\left(t\right)\right)& \mathrm{sign}\left(t\right)=0\end{array}\right.---\left(26\right).$

Images