塔式太阳能聚热发电系统镜场天空云运动估计研究

李军，张鹏，聂立，赵跃，刘涵

（1.青海宁北发电有限责任公司唐湖分公司，青海海北 810299；2.西安理工大学自动化与信息工程学院，陕西西安 710048；3.东方电气集团东方锅炉股份有限公司，四川自贡 643001）

太阳能是地球上一切能源的源泉，随着技术的进步人们对能源的需求量越来越大，因此对太阳能的开发利用显得极为紧迫。太阳能发电技术主要分为光伏发电与聚热发电两种[1]。光伏发电由于在生产电池时所用的金属对环境污染大且难于并网、无法持续发电等缺点导致其逐渐被聚热发电系统所替代。太阳能聚热发电系统相对于光伏系统具有可持续发电、污染小、能效高等优点。

太阳能聚热发电系统分为塔式、槽式、碟式三种，其中塔式系统具有聚光比大、使用寿命长等优点目前被广泛应用。塔式太阳能热发电系统主要由定日镜场、吸热器、储热系统、汽轮机和发电机5个部分组成[2]。工作时通过定日镜将太阳光线反射到吸热器上，吸热器内部的热导介质吸收太阳能并产生蒸汽，然后蒸汽推动汽轮机发电，当吸热器内部温度变化过快时会产生较大的热应力[3]。对于镜场而言，云层的遮挡是导致辐射量变化的决定性因素，当云层突然进场或离开时，会使镜场接受到的辐射量突变，并引起吸热器内部的温度骤变，从而使吸热器受到热应力损害。因此，有必要对镜场的云运动进行估计，并据此对定日镜的姿态做出调整，这对于整个系统换热的安全与转换效率具有十分重要的意义。

目前，国内外针对镜场天空云运动预测的研究并不多，且大多数都是集中在全天空背景下，而国外学者对于此项研究则较为深入，例如Ricardo Marquez等人基于TSI全天空成像仪对影响太阳辐射的云团路径和辐射的测量数据进行了预测，实现了未来3～15 min内的太阳直接辐射预报[4-6]。但是基于硬件的全天空设备造价高昂，因此研究一套低成本的云层预测系统就显得十分必要。本文采用阈值分辨法将图像分为太阳与云和背景两组，将图片中的太阳区域去除，然后采用图像的匹配技术对云层进行匹配，通过对两张或两张以上的图片寻找特征点进行匹配，最后采用Kalman滤波算法完成云层运动估计和预测。

1 图像处理算法

对于定日镜镜场天空云层的预测，首先需要得到云层所在的区域，因此需要基于光学理论采用图像处理技术对图像进行分析、加工和处理。对定日镜场云图进行处理的流程图如图1所示。

图1 图像处理流程图Fig.1 Flow chart of the image processing

1.1 图像预处理

对云层的运动方向进行预测，必须要通过对天空图像中云层运动的趋势来判定。本文的算法是对镜场云图固定的等时帧间隔图像予以截屏，通过取其连续5幅图像实现对云层运动进行预测。

对图像的预处理是采用自适应中值滤波算法对图像滤波去噪。因Surf匹配算法实现特征点匹配时，对特征点的提取主要是集中在边缘轮廓，所以预处理滤波的主要目的便是既能消除噪声又能保持边缘细节，而中值滤波对脉冲噪声具有良好的滤除作用，特别是在滤除噪声的同时，能够保护信号的边缘，使之不被模糊。因此，对所有输入图像（5幅）都先转换成灰度图像后进行自适应中值滤波处理。

g（x，y）=median{（f（x-i，y-j））}，（i，j）∈W（1）式中：g（x，y）和f（x，y）分别代表滤波处理后和处理前图像；W为二维模板。自适应中值滤波通过判定区域中值是否为极值调整模板大小，对模板中心点为极值的点滤波处理，然后选定0.95作为阈值将灰度图像转换为二值图像[7]。如图2（a）所示，光感最强烈区域为太阳所在范围[8]，此时所得二值图像除了太阳以外区域已全部变为背景，通过计算太阳形心位置分割后利用形心位置去除太阳所在区域。这样将太阳区域RGB三通道赋值其余区域进行三通道复原，便得到了去除太阳的图像如图2b所示。

图2 太阳区域定位图Fig.2SunZoneLocationImage

1.2 云层图像分块处理

每块云的运动方向及速度各不相同，所以在进行特征点匹配后，要将匹配出的特征点分区以得到每片云层区域内的特征点。对特征点分区便需要对云层分块，通过对图像进行分割并把分割出的云块进行标记，如图3所示。本文采用纹理分割结合Otsu算法[9]（最大类间方差法）的方法将二值图像分割为若干连通区域，对若干连通区域分别进行标定如图4所示。

图3 二值化图像Fig.3 Binaryzation image

图4 区域标定图Fig.4 Region segmentation image

如图4所示，对太阳形心位置所在的连通域设置为背景颜色，这样便得到了云层连通域。在此也尝试过其他方法的分割如聚类算法，但因天空图像中光晕与背景天空和云层分割时容易产生混淆，因此所得效果并不理想，如图5所示。由此可以看出本文算法相对聚类算法的优越性。

由于天空中的云形变不可预知，因此我们只针对第一帧图像进行分割，用第一幅图像中的特征点坐标来判定每一个匹配点所处图像中的区域。将图像分区后对区域进行编号，便可对后续图像匹配生成功的匹配对分区处理，以方便对每个区域内的运动趋势计算。

图5 K-means分割图Fig.5 The K-means image segmentation

2 基于Surf算法的特征点匹配

对天空图像处理过后需要对天空云团进行特征点匹配，利用点位的位置信息确定这些云层的运动趋势。本文采用Surf算法[10]实现匹配。Surf（Speeded up robust feature）算法所检测到的特征点具有较强的旋转不变性且对光照条件的鲁棒性较强，因此不会因天空光线的变化而产生过多误匹配现象。为了找到两幅图像对应的云，需要分别提取每幅图像中的稳定点，这些点之间会有相互对应的匹配点。该算法对特征点的提取具有较强的稳定性、鲁棒性，因此在对云层的匹配中利用此算法对特征点的定位具有较好的效果。该算法主要利用积分图像、近似的Hessian矩阵和Haar小波变换，对连续图像中特征点分3步处理：

1）将第一帧图像分别与后4帧连续图像特征匹配，得到第一帧图像的4组匹配成功的特征点；

2）将第一幅图像的4组特征点中相同特征点保留下来，建立起每个特征点在5副图像中的对应关系；

3）通过对云图的分块信息将每组中的特征点根据云图信息分区，这样便得到了每一区域云图中的匹配成功的特征点个数及其在各帧中的坐标等信息。

2.1 箱式滤波器与极值点定位

首先需要输入5幅待匹配图像，并将其转换到灰度空间。为了匹配出5副图像中的特征点对，首先要找到其特征点，而特征点的出处则是极值点，因此利用积分图像结合箱式滤波器找出每一副图像中的极值点，积分图像的意义是图像中任意一点的灰度值，I（i，j）为原图像左上角到（i，j）点相应的对角线区域灰度值的总和，即：

图6所示为灰色区域内的像素值之和。

图6 积分图像Fig.6 Integral image

A，B，C，D各点的积分值进行加减，所得即是S区域内积分图像像素值之和。当得到积分图像后，利用如图7所示的箱式滤波器通过改变自身大小对图像卷积，便可构造出尺度空间。卷积时，由于是矩形区域与像素对应相乘求和，因此可转换为像素求和然后乘以矩形区域中的系数。所以采用箱式滤波器做卷积，就变成了像素求和（积分图像），因此可大幅提高算法的速率从而保证了算法的实时性。对于给定图6中的已知一点，Dyy（X，σ），Dxx（X，σ）是对其卷积结果，X是位置、σ为尺度，下标为不同方向箱式滤波器与图像在该处的卷积，其滤波后得到式（4）的Hessian矩阵：

图7 9×9箱式滤波器（左为Dyy，右为Dxy）Fig.7 9×9 Box filter(left：Dyy;right：DXX)

计算其Hessian矩阵行列式，由（5）所示，其中w为权值。根据Hessian矩阵的性质，若检测点所得行列式值为正则将其判定为局部极大值点并归纳于兴趣点范围。通过改变盒子滤波器的尺寸对图像进行处理构成尺度空间，便得到了每幅图像在不同尺度空间当中的局部极大值点集。

2.2 尺度空间与精确定位

获得5幅图像的极值点信息后便要从中提取特征点，对每幅图像所生成的尺度空间分成几组（Octaves），分别对每一组分层（scale）。其中每一组当中的每一层都是同一副图像经不同尺寸的箱式滤波器滤波处理后的图像，如图8所示，图中的数字为滤波器尺寸的大小。

图8 尺度空间图Fig.8 Scale space figure

第一组的滤波器尺寸大小由9×9，增大到27×27，每次增大6像素；第二组由15开始每次增加12（每一组的第一个滤波器与上一组第二个滤波器尺寸相同，相邻两组层步长增加值为二倍关系），通过滤波模板生成尺度空间与极值点后。为了确定其中的特征点及其位置信息，以检测点为中心构建了3×3×3的在三维空间内进行非极大值抑制，只有当检测点的Hessian行列式值大于固定阈值T，且是三维空间内的极大值点才将其确定为特征点，其后通过对特征点进行三维线性插值确定其精确位置，如图9所示。

图9 非极大值抑制图Fig.9 Non-maximum inhibition of figure

2.3 构造描述子

获得特征点后需要对其匹配，找到这个点在5幅图像中的位置，从而得到其运动、位置关系。为了实现匹配，需要构建5幅图像中所有特征点的描述信息。因云层在运动过程中可能发生旋转，为了实现特征点的旋转不变性，首先需要确定其主方向信息。将特征点作为中心，圈取一个半径为6倍尺度值大小的圆形区域，对整个邻域内计算Haar小波分别在垂直和水平方向的响应值。根据其与中心点的远近进行加权计算，然后对响应进行统计。将圆形区域内每60°内的所有响应值求和得到一个矢量值。然后将步长设定为5°，每隔5°进行求和计算，直到覆盖整个区域为止，将矢量值最大的方向确定为主方向。确定主方向后需要得到特征点的描述子。在特征点周围取一个与主方向同向的正方形框，框的边长大小为20s（s为特征点所在的尺度值），并将其分为大小相等的16个子块，用尺寸大小为2s的Haar小波模板对子块滤波处理，得到（以特征点主方向建立的坐标轴）水平及垂直方向的响应值dx、dy，后生成如式（6）所示V子块内的4维特征矢量统计和。由上所述每个特征点可以得到包含16×4个特征向量，这样便得到了每幅图像内的特征点。如图10所示，对单幅天空图像进行的特征点进行了提取，由图可以看出其包含了特征点的位置、主方向信息。

图10 特征点提取图Fig.10 Feature point extraction

2.4 特征点匹配

由上所述得到了每幅图像中的特征点及其特征信息，接下来便要对所有特征点进行匹配，寻找每两幅图像中最为相似的点。对每一个特征点求其Hessian的迹，如图11所示。根据其亮度将其分为两种，一种是特征点所在的圆形区域比背景亮度更亮；另一种特征点所在的圆形区域比背景亮度更暗。通过实验得出当Hessian矩阵的迹为正时得到的是亮点区域，而当Hessian矩阵的迹为负值时得到的是暗点区域。倘若待匹配两点的迹为同号，说明具有相同的对比度则为可能的匹配对；倘若迹为异号则说明对比度较低，不再进行后续处理。对同号的进行下一步相似度测量，通过对上节所述生成的64维描述子进行欧式距离计算得到待匹配点的相似程度如式（7）所示。找到一个特征点另一幅图像中距离最小和次最小的两点的欧式距离。倘若两个欧式距离的比值小于设定阈值则将与待匹配点欧式距离最小的点确定为匹配点。

式（7）中，两个待匹配的特征点集分别含有i和j个点，第一个点集中i点的第k个特征描述子为Xik，同理可得Xjk，n表示特征向量的维数。因云层运动缓慢，所以匹配时加入了位移的限定阈值，当两帧图像同一特征点移动超过30像素时认定为误匹配，在此使用Sift算法[11]与Surf算法分别对前后两帧图像（图像大小为800×450）进行匹配得到下表所示结果：

图11 对比亮度图Fig.11 Brightness contrast figure

表1 Sift与Surf算法的比较Tab.1 ComparisonofSiftandSurfalgorithmsinthispaper

由此可以看出，经Surf匹配结果的速率与准确率及有效匹配对数量均高于Sift算法。对处理过的5幅图像采用Surf算法进行特征点匹配，对4组匹配后的坐标重新排序，得到同一特征点在5幅图像中的具体位置。通过对匹配成功的匹配对计算，可以得到每一个特征点在图片中的位移。相邻两幅图像的时间间隔取为5 s，图12为对其中两帧图像的匹配，从匹配中可看出未出现误匹配现象。这样就可得到多副图像中每个特征点所在不同帧数的具体位置，及相邻两帧间特征点的具体位移。

图12 图像匹配图Fig.12 Adjacent image matching

3 基于Kalman滤波器的云运动估计

为了预测云块区域的运动速度，需要得到云层区域之前的运动趋势。在特征点的匹配中已经得到了过去连续5帧内各区域中的每个特征点的坐标信息，对前后两帧每个区域内所有点依次计算坐标位移，因对区域内所有点取均值后计算坐标位移具有同等作用，所以先对区域内的点位坐标取均值作为整体区域的坐标。因Surf方法在寻找兴趣点时利用了线性插值，因此所匹配到的特征点会有最大1像素的误差，当云团运动过慢时会产生较大的误差值。为了消除误差，在预测时选用了Kalman滤波器，利用Kalman滤波算法对云层的运动做出最优估计[12]。

3.1 Kalman滤波器

卡尔曼（Kalman）滤波是基于状态空间方法的递推滤波算法，卡尔曼滤波器主要由状态空间模型和观测模型所组成。状态模型对系统的运动规律进行控制，观测模型为系统提供修正信息，通过不断的迭代修改使系统得到最优的运动估计信息。在对云层速度的预测时设动态系统的状态空间和观测模型为：

因为云层运动缓慢，因此短时内将每块云层运动看做匀速运动，其中X（t）为系统在t时刻的状态，将其设置为二维状态包括坐标及速度信息，此处将第一幅图像中各区域坐标（x，y方向分开计算）及速度作为状态变量，初始速度为第二帧与第一帧图像坐标之差，初始坐标为第一幅图像中的坐标；Y（t）为状态的观测值，即连续几帧图像对应区域的坐标值；V（t）为观测噪声，由于像素误差为1，故将方差R也设置为1；Ф为状态转移矩阵；H为观测矩阵。

3.2 Kalman滤波器实现云层运动的最优估计

将图像匹配分割后，选定一个云块。这里对输入的5幅图像进行卡尔曼滤波的计算，分别以第1幅，第2幅，第3幅图像为作为初始图像，3幅连续图像作为图像序列进行卡尔曼滤波，因此可得到3组运动速度值。

考虑到3幅图像为一组，速度误差较大，为减少误差，对卡尔曼滤波后所得的坐标和速度值，再次利用卡尔曼滤波器进行滤波处理，然后对上述所得三组速度值进行一维滤波，滤波后得到3组速度的最优预测值。这样便可得到所取5幅图像的下一时刻每个云块x，y方向的运动速度的最优预测值，这时云层的运动方向也可计算求得，这样就实现了天空图像中云层的短时运动估计。

4 实验及结果分析

为验证算法，实验时采用自制的Flash视频。以天空图像及太阳区域作为背景，以云层作为运动元件通过在Flash中加入监控程序得到云层在每一帧的坐标位置，以方便后续的验证工作。利用Flash程序使云层的运动保持在x方向上每帧3.15像素、y方向上每帧负4.25像素的匀速运动。对Flash中连续截屏得到连续帧，图像其中一帧如图13所示，得到每帧图像并得到云层在连续帧中的位置集合。对所截图集进行特征点匹配，匹配后对特征点的位置信息进行滤波处理得到运动速度，再对速度滤波处理后得到对云层的最优预测速度。图14所示为匹配程序对两帧截图的匹配效果。这样便可通过云图实际移动速度与匹配滤波后得到的预测速度作比较，得到匹配的误差率。实验中Flash图像大小为550×381，且采用单片云作为运动元件，所以图像中的特征点信息并不如实际天空图像中一样丰富，匹配后得到了75个匹配特征点。

图13 Flash截图Fig.13 Flash screenshot

图14 特征匹配图Fig.14 Feature matching figure

表2 实际值与使用两次卡尔曼滤波的预测值和观测值的比较Tab.2 The actual value compared with the value with Kalman filter prediction used twice

利用特征点匹配结合卡尔曼滤波的算法对云的运动进行预测，实际情况下太阳作为背景并没有运动，因此可将运动速度在0左右的区域归纳为太阳区域，从而对其余区域的运动进行记录。表2为云层实际速度与使用Surf算法结合卡尔曼滤波器计算的速度值对比，表2中第一组数据是对表3所得的4组数据滤波处理后所得。而表3为未使用卡尔曼滤波器、仅使用Surf算法通过帧间位移计算得到的预测值，误差率最高可达29%。表2和表3结果表明，经过两次卡尔曼滤波后对运动速度预测的精度有了大幅提高。经大量数据测算比对后得出本方法可以将速度预测误差控制在7.5%以内，将误差率减少了21.5%。

表3 实际值与未加入卡尔曼滤波器测量值比较Tab.3 The actual value compared with the value with Kalman filter not used

表4 速度放大4倍后实际值与预测值比较Tab.4 The actual value compared with the prediction value after the speed is amplified by four times

表4是表2中速度放大4倍后得到的实际值与预测值的比较，此时的误差已可控制在5%以内。经过对比可知，当云层运动速度更快时，可以获得更好的预测效果。由于本方法对天空图像进行了分块处理，所以当天空中有多片云层时可得到每一片云层的运动估计信息，针对理想天空环境下可实现误差精度在7.5%以内的镜场天空云层运动估计，可为塔式太阳能聚热发电定日镜场提供较为可靠的前馈控制信息。

5 结语

本文提出了一种基于数字图像处理技术的云层运动估计方法，采用自适应阈值分割的方法对云层进行了分区处理，通过特征点匹配的方法得到了云层的运动信息，然后结合两次卡尔曼滤波器方法得到了在理想天空条件下对镜场天空中每一块云层的运动进行了监测预测，最后通过对云层运动的监测和预测验证了方法的有效性，该方法可为塔式太阳能聚热发电定日镜场的协调控制策略提供可靠的前馈控制信息。

[1]刘静静，杨帆.太阳能热发电系统的研究开发现状[J].电力与能源，2012，33（6）：573-576.LIU Jingjing，YANG Fan.Present situation of research and development of solar thermal power generation system[J].Power and Energy，2012，33（6）：573-576（in Chinese）.

[2]廖葵，龙新峰.塔式太阳能热力发电技术进展 [J].广东电力，2007，20（4）：6-11.LIAO Kui，LONG Xinfeng.Technology progress in solar power tower[J].Guangdong Electric Power，2007，20（4）：6-1（in Chinese）.

[3]王建楠.塔式太阳能热发电站熔融盐吸热器故障诊断专家系统研究[D].北京：中国科学院研究生院，2010：10-13.

[4]MARQUEZ R，COIMA C F M.Intra-hour DNI forecasting based on cloud tracking image analysis[J].Solar Energy，2013，91：95-99.

[5]漆随平，刘涛.地基云观测技术及装备研究进展[J].山东科学，2014，6（27）：1-8.QI Suiping，LIU Tao.Ground cloud observation technology and its instrument research advances[J]，Shandong Science，2014，6（27）：1-8（in Chinese）.

[6]M BLANCO-MURIEL，D.C.ALARCON-PADILLA.Computing the solar vector[J].Solar Energy，2000，70（5）：431-441.

[7]霍娟，吕达仁.全天空云识别阈值法的数值模拟初步研究[J].自然科学进展，2006，4（16）：480-484.HUO J，LÜ D.The preliminary research on numerical simulation about the whole-sky cloud threshold recognition[J].Progress Innatural Science，2006，4（16）：480-484（in Chinese）.

[8]陶小红，胡以华.基于CCD图像的太阳定位技术[J].能源研究与信，2005，21（1）：50-53.TAO Xiaohong，HU Yihua.A technique for sun location based on CCD imaging[J].Enerty Reserach and Information，2005，21（1）：50-53.

[9]王祥科，郑志强.Otsu多阈值快速分割算法及其在彩色图像中的应用[J].计算机应用，2006，6（26）；14-16.WANG Xiangke，ZHENG Zhiqiang，Otsu multi threshold fast segmentation algorithm and its application in color images[J].Computer Applications，2006，6（26）；14-16.

[10]BAY H，ESS A，et al.Speeded up robust features[J].Computer Vision and Imaging，2008，110（3）：346-359.

[11]DAVIS G lOWE.Distinctive image features from scale-invariant key points[J].International Journalof Computer Vision，2004，60（2）：91-110.

[12]KALMAN R E.A new approach to linear filtering and prediction problems[J].Transactions of the ASME-Journal of Basic Engineering，1960，82（Series D）：35-45.