CN112712542A - 一种块匹配与光流法结合的地基云图运动预测方法 - Google Patents

Abstract

本发明公开了一种块匹配和光流法结合的地基云图运动预测方法，针对获取的地基全天空图像，利用云分类信息提取云区域，采用块匹配法计算云区域的主运动矢量，将其结果作为光流法的约束项之一，和亮度恒定项、全局平滑项和亮度梯度项共同构造光流法能量函数，同时在光流法云运动场计算过程中，对每层金字塔加入加权中值滤波，减小运动矢量异常值的影响，实现地基全天空图像的云运动快速准确的预测。本方法针对块匹配预测精度低，光流法计算时间长、光线敏感等问题，通过将块匹配和光流法结合，充分利用块匹配速度快、光流法精度高的优势，提高了全天空图像的云运动预测的速度和精度。

Description

一种块匹配与光流法结合的地基云图运动预测方法

技术领域

本发明属于摄影测量与遥感领域，尤其涉及一种块匹配与光流法结合的地基云图运动预测方法。

背景技术

在短期内到达地面的太阳辐射量可能发生大的波动，云层的遮挡是导致辐射量变化的关键性因素，当云层和太阳的相对位置发生变化时，会造成辐射量的剧烈波动，因此判断短期内云的运动位置和方向是预测太阳辐射量变化的关键。获取云观测数据的地基全天空成像仪安装在某个固定位置，可获取该地区一定范围的全天空图像，其时间分辨率一般为30秒，利用该仪器可实现云识别、云运动场预测和短期太阳辐照度预测。

目前常用的运动目标检测方法包括块匹配法和光流法。

块匹配方法首先将图像划分为若干块，假定每一个块内的所有像素作为一个整体，具有相同的运动模式，围绕某块B₁建立大范围的搜索窗口，在后一幅图像的搜索窗口内，按照一定的搜索原理寻找与块B₁最相似的块B₂，它们中心点在水平和垂直方向上的差值即为块B₁的运动矢量。块匹配的思路简单易懂，对硬件要求低，运算速度快，因其高效的优点而被广泛采用，但算法具有易陷入局部最优的固有缺陷，并且将整块作为一体的假设在一定程度上限制了块匹配算法的精度。

光流被定义为空间运动目标在二维成像平面上的像素运动的瞬时速度，是图像中亮度的表观运动，图像中所有像素点的光流集合组成了光流场，在理想情况下，光流场对应于运动场。光流法不需要预先知道场景的任何信息就能精确计算出目标的运动矢量，但它对光线的变化十分敏感，同时由于采用迭代计算的方式，耗时稍长。

为了克服块匹配精度低、光流法计算时间长且对光线敏感的不足，本发明将块匹配和光流法结合，将块匹配计算出的主运动矢量作为约束项代入光流法，减少光流法的迭代次数，克服光线导致的图像亮度突变，同时在光流计算中加入加权中值滤波，去除云边缘运动矢量异常值，提高计算精度。与单个方法相比，本发明具有云运动预测精度高、速度快、鲁棒性强的优点。

发明内容

本方法主要解决块匹配对云运动预测精度低和光流法预测计算时间长、光线敏感的问题，提出了一种块匹配和光流法结合的地基云图运动预测算法。

本发明核心过程包括以下步骤：

步骤1，读取原始地基全天空图像，首先利用平面结构引导的修复算法对原始地基全天空图像中被遮挡的区域进行图像修复，再以图像中心为原点、设置一定半径保留图像天空区域、裁掉与天空无关区域，最后根据角度的对应关系进一步进行图像畸变校正，得到第k幅和第k+1幅预处理后地基全天空图像I_k,I_k+1,k∈[1,K-1]，K为预处理后地基全天空图像的数量。

步骤2，通过基于多特征的随机森林分类方法将预处理后地基全天空图像的每个像素分为云类别或者晴空类别，从而将预处理后地基全天空图像转换成云或晴空的二值图像，利用逆四叉树法在二值图像I_bw中确定云范围；

步骤3，根据预处理后地基全天空图像中方块坐标的像素亮度值构建均方误差优化目标，以均方误差最小化为优化目标进行块匹配，得到预处理后地基全天空图像的运动矢量；

步骤4，结合步骤2所述的预处理后地基全天空图像的云范围、步骤3所述的块匹配运动矢量得到预处理后地基全天空图像的云范围内的块匹配运动矢量，筛选出频次较高的多个块匹配运动矢量作为主运动矢量；

步骤5，针对步骤2中所述的相邻两幅预处理后地基全天空图像的总云范围内每个像素，分别计算像素点运动矢量与步骤4所述的运动矢量的矢量差，以矢量差的模最小为准则将对应的主运动矢量赋值给该像素，从而重新确定云范围内的运动矢量，得到云的初步运动矢量；

步骤6，根据云的初步运动矢量构造主运动项，连同亮度恒定项、全局平滑项、亮度梯度项和加权中值滤波构建光流法能量函数，对能量函数最小化进行优化求解，得到全天空图像中每个像素点的运动矢量；

步骤7，针对基于微分的光流算法只能计算小位移的光流问题，引入金字塔思想进行由粗到细的多尺度计算，以一定比例将预处理后全天空图像降采样为N个不同分辨率的图像，在金字塔中从上到下按最低分辨率到最高分辨率依次排列，从金字塔顶层开始，对第n层图像进行步骤6的光流法计算，得到运动矢量(u,v)ⁿ，然后作为第n+1层的矢量计算初始值(u₀,v₀)ⁿ⁺¹，直到最底层N，从而将大位移分解为若干个小位移的集合，达到计算大位移光流的目的；

作为优选，步骤2所述预处理后地基全天空图像的每个像素为：

I_k(i,j)

其中，I_k(i,j)表示第k幅预处理后地基全天空图像中第i行第j列的像素，i∈[1,M],j∈[1,N],k∈[1,K-1]，K为预处理后地基全天空图像的数量，M*N为预处理后地基全天空图像的分辨率；

步骤2所述晴空类别定义为C_sky；

步骤2所述云类别定义为C_cloud；

步骤2所述二值图像，具体定义为：

I_k,bw(i,j),type_k(i,j)

其中，I_k,bw(i,j)表示第k幅预处理后地基全天空图像对应的二值图像I_bw中第i行第j列的像素，i∈[1,M],j∈[1,N],k∈[1,K-1]，K为预处理后地基全天空图像的数量，M*N为预处理后地基全天空图像的分辨率，type_k(i,j)表示第k幅预处理后地基全天空图像对应的二值图像中第i行第j列的像素类别，若I_k,bw(i,j)为晴空类别则type_k(i,j)＝C_sky，若I_k,bw(i,j)为云类别则type_k(i,j)＝C_cloud；

步骤2所述利用逆四叉树法在二值图像中确定相邻两幅预处理后地基全天空图像的总云范围，具体为：

对二值图像进行初步处理：

以n*n的像素方阵作为基本单元即U_k,1(i′,j′)；

U_k,1(i′,j′)为第k幅预处理后地基全天空图像对应的二值图像中第i′行第j′列的第1级方阵，即基本单元，方阵的尺寸为n*n，i′∈[1,M/n],j′∈[1,N/n],k∈[1,K-1]，K为预处理后地基全天空图像的数量；

当U_k,1(i′,j′)内属于云类别的像素个数超过单元总像素数n*n的一定比例时，则将U_k,1(i′,j′)内所有像素设置为云类别C_cloud，否则全部设置为晴空类别C_sky；

将基本单元按照四叉树原理自下而上逐步聚拢，得到第k幅全天空图像不同尺寸的云块单元U_k,1,U_k,2,...,U_k,m,...,U_k,L，其中第m级的方阵尺寸为2^m-1n*2^m-1n；

第k幅预处理后地基全天空图像的云范围R_k为第k幅预处理后地基全天空图像对应的二值图像中所有云块单元U_k,1,U_k,2,...,U_k,m,...,U_k,L的集合；

对第k+1幅预处理后地基全天空图像I_k+1重复上述操作得到对应的云范围R_k+1，则第k幅和第k+1幅预处理后地基全天空图像的总云范围为：

R_k,k+1＝R_k∪R_k+1，R_k,k+1为R_k和R_k+1的并集。

作为优选，步骤3所述构建均方误差优化目标为：

其中，

为第k幅全天空图像中第i行第j列像素的运动矢量，i∈[1,M],j∈[1,N],k∈[1,K-1]，K为预处理后地基全天空图像的数量，M*N为预处理后地基全天空图像的分辨率；方块的大小为N′*N′，f_k+1(n+u,n+v)表示第k+1幅图像中方块坐标(n+u,n+v)处的像素亮度值，f_k(n,n)则表示第k幅参考图像搜索窗口中方块坐标为(n,n)处的像素亮度值，u和v分别为横向和纵向的运动距离，u∈[-a,a],v∈[-a,a]，a为块匹配搜索窗口的大小。

作为优选，步骤4所述相邻两幅预处理后地基全天空图像的总云范围内像素的块匹配运动矢量为：

其中，第k幅全天空图像第i行第j列像素为云类别C_cloud，其运动矢量为

R_k,k+1为步骤2所述第k幅和第k+1幅预处理后地基全天空图像的总云范围，M*N为预处理后地基全天空图像的分辨率，K为预处理后地基全天空图像的数量；

步骤4所述筛选出频次较高的多个矢量作为主运动矢量为：

统计第k幅全天空图像中每对

出现的频次

选取频次

最高的多个运动矢量作为该图像云层的主运动矢量Ω：

num为主运动矢量的个数。

作为优选，步骤5所述分别计算像素点运动矢量与步骤4所述的运动矢量的矢量差的模为：

其中，

和

分别为第k幅和第k+1幅地基全天空图像的总云范围内第i行第j列像素的块匹配运动矢量的两个分量，

和

分别为第k幅全天空图像第i″个主运动矢量Ω的两个分量；

步骤5中所述以矢量差的模最小为准则将对应的主运动矢量赋值给该像素为：

其中，

为第k幅全天空图像第i行第j列像素的初步运动矢量，(i,j)∈R_k,k+1,k∈[1,K-1]，R_k,k+1为步骤2所述第k幅和第k+1幅预处理后地基全天空图像的总云范围，K为预处理后地基全天空图像的数量。

作为优选，步骤6所述亮度恒定项为：

光流法的基本约束方程，假设同一像素点在连续两幅全天空图像I_k和I_k+1中的亮度值不发生变化，即灰度值不变，其中k∈[1,K-1]，K为地基全天空图像的数量；

E_intensity(u,v)＝∫∫ψ(|I_k+1(x+dx,y+dy)-I_k(x,y)|²)dxdy

其中，I_k(x,y)为某像素点(x,y)在第k幅图像中的亮度值，经过一定时间，该像素点移动了(dx,dy)，此时的亮度值为I_k+1(x+dx,y+dy)。式中

ε＝0.001为Charbonnier惩罚函数，(u,v)表示计算得到的该像素点(x,y)的运动矢量。

步骤6所述全局平滑项为：

整个图像I_k上光流的变化是平滑的，即u和v随着像素点移动而发生的变化是缓慢的；

步骤6所述亮度梯度项为：

由于太阳周围云的运动可能造成图像晴空背景亮度发生较大变化，使得像素亮度值的匹配出现错误，可通过亮度值梯度的约束来增强亮度恒定项的效果，达到进一步精确地预测出后一幅云图I_k+2的目的；

步骤6所述主运动项为：

将块匹配法获得的云初步运动矢量

作为约束项代入光流法，使得光流计算的迭代次数可适当降低，提高计算效率，同时减弱各云块内部亮度不均一的影响；

步骤6所述加权中值滤波为：

对运动矢量(u,v)进行加权中值滤波能有效改善云边缘处异常值的问题；

首先利用Sobel算子检测云边缘并适当扩大范围，对于非边缘的区域，以[m,m]的窗口大小进行传统中值滤波，而在云边缘区域，综合考虑空间距离、Lab颜色空间差距以及遮挡状态这三个因素，在[M,M]的窗口范围内进行加权中值滤波，且M＞m；

其中，N_x,y是以(x,y)为中心的窗口范围内除中心外的点(x′,y′)的集合，

代表两个像素(x,y)和(x′,y′)之间相似性的权重，具体为下式，I代表Lab空间中的颜色向量，n_c为颜色波段数，σ₁＝7,σ₂＝7：

遮挡状态变量o(x,y)被定义为

式中d(x,y)为单边散度函数，div(x,y)代表光流的散度

步骤6所述的光流法能量函数为：

式中，α,γ,β分别为全局平滑项的权重、亮度梯度项的权重、主运动项的权重；

以光流法能量函数最小化为优化目标，利用欧拉-拉格朗日方程、逐次超松弛法等数学手段求解，得到全天空图像I_k每个像素点的运动矢量(u,v)；

作为优选，步骤7所述一定比例为：r,0＜r＜1；

步骤7所述N个不同分辨率的图像为：

本发明具有如下优点：通过综合利用块匹配运算速度快、光流法计算精度高的优点，同时在光流计算中引入的加权中值滤波，有效去除云边缘运动矢量的异常值，实现全天空图像云运动精确快速预测，该方法还具有鲁棒性，克服光流法对光线敏感的问题，对图像亮度突变的情况也能实现精确预测。

附图说明

图1为使用本发明提出的块匹配和光流法结合的地基云图运动预测的作业流程；

图2为图像的预处理过程，包括图像修复和图像校正。

图3为一组连续拍摄的全天空图像，拍摄时间间隔为1分钟。

图4为基于云分类信息采用块匹配计算的主运动矢量。

图5为根据块匹配结果，采用光流法计算出的运动矢量和预测图像。

图6为光流法中无滤波、中值滤波和加权中值滤波的对比图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

下面结合图1至图6介绍本发明的具体实施方式为一种块匹配与光流法结合的地基云图运动预测方法：

实施路线如图1所示：

步骤1，读取原始地基全天空图像，首先利用平面结构引导的修复算法对原始地基全天空图像中被遮挡的区域进行图像修复，再以图像中心为原点、设置一定半径保留图像天空区域、裁掉与天空无关区域，最后根据角度的对应关系进一步进行图像畸变校正，如图2。

得到的第1幅和第2幅预处理后地基全天空图像I₁,I₂，第3幅图像为预测时刻的真值全天空图像，如图3。

步骤2，通过基于多特征的随机森林分类方法将预处理后地基全天空图像的每个像素分为云类别或者晴空类别，从而将预处理后地基全天空图像转换成云或晴空的二值图像，利用逆四叉树法在二值图像I_bw中确定云范围，如图4(a)；

步骤2所述预处理后地基全天空图像的每个像素为：

I₁(i,j)

其中，I₁(i,j)表示第1幅预处理后地基全天空图像中第i行第j列的像素，i∈[1,250],j∈[1,250]，250*250为预处理后地基全天空图像的分辨率；

步骤2所述晴空类别定义为C_sky；

步骤2所述云类别定义为C_cloud；

步骤2所述二值图像，具体定义为：

I_1,bw(i,j),type₁(i,j)

其中，I_1,bw(i,j)表示第1幅预处理后地基全天空图像对应的二值图像中第i行第j列的像素，i∈[1,250],j∈[1,250]，250*250为预处理后地基全天空图像的分辨率，type₁(i,j)表示第1幅预处理后地基全天空图像对应的二值图像中第i行第j列的像素类别，若I_1,bw(i,j)为晴空类别则type₁(i,j)＝C_sky，若I_1,bw(i,j)为云类别则type₁(i,j)＝C_cloud；

步骤2所述利用逆四叉树法在二值图像中确定相邻两幅预处理后地基全天空图像的总云范围，具体为：

对二值图像进行初步处理：

以5*5的像素方阵作为基本单元即U_1,1(i′,j′)；

U_1,1(i′,j′)为第1幅预处理后地基全天空图像对应的二值图像中第i′行第j′列的第1级方阵，即基本单元，方阵的尺寸为5*5，i′∈[1,50],j′∈[1,50]；

当U_1,1(i′,j′)内属于云类别的像素个数超过单元总像素数5*5的一定比例时，则将U_1,1(i′,j′)内所有像素设置为云类别C_cloud，否则全部设置为晴空类别C_sky；

将基本单元按照四叉树原理自下而上逐步聚拢，得到第1幅全天空图像不同尺寸的云块单元U_1,1,U_1,2,...,U_1,m,...,U_1,L，其中第m级的方阵尺寸为(2^m-1*5)*(2^m-1*5)；

第1幅预处理后地基全天空图像的云范围R₁为第1幅预处理后地基全天空图像对应的二值图像中所有云块单元U_1,1,U_1,2,...,U_1,m,...,U_1,L的集合；

对第2幅预处理后地基全天空图像I₂重复上述操作得到对应的云范围R₂，则第1幅和第2幅预处理后地基全天空图像的总云范围为：

R_1,2＝R₁∪R₂，R_1,2为R₁和R₂的并集。

步骤3，根据预处理后地基全天空图像中方块坐标的像素亮度值构建均方误差优化目标，以均方误差最小化为优化目标进行块匹配，得到预处理后地基全天空图像的运动矢量，如图4(b)；

步骤3所述构建均方误差优化目标为：

其中，

为第1幅全天空图像中第i行第j列像素的运动矢量，i∈[1,250],j∈[1,250]，250*250为预处理后地基全天空图像的分辨率；方块的大小为25*25，f₂(n+u,n+v)表示第2幅图像中方块坐标(n+u,n+v)处的像素亮度值，f₁(n,n)则表示第1幅参考图像搜索窗口中方块坐标为(n,n)处的像素亮度值，u和v分别为横向和纵向的运动距离，u∈[-50,50],v∈[-50,50]，50为块匹配搜索窗口的大小。

步骤4，结合步骤2所述的相邻两幅预处理后地基全天空图像的总云范围、步骤3所述的块匹配运动矢量得到相邻两幅预处理后地基全天空图像的总云范围内的块匹配运动矢量，如图4(c)。筛选出频次较高的多个块匹配运动矢量作为主运动矢量；

步骤4所述相邻两幅预处理后地基全天空图像的总云范围内像素的块匹配运动矢量为：

其中，第1幅全天空图像第i行第j列像素为云类别C_cloud，其运动矢量为

R_1,2为步骤2所述第1幅和第2幅预处理后地基全天空图像的总云范围，250*250为预处理后地基全天空图像的分辨率。

步骤4所述筛选出频次较高的多个矢量作为主运动矢量为：

统计第1幅全天空图像中每对

出现的频次

选取频次

最高的多个运动矢量作为该图像云层的主运动矢量Ω：

num为主运动矢量的个数。

步骤5，针对步骤2中所述的相邻两幅预处理后地基全天空图像的总云范围内每个像素，分别计算像素点运动矢量与步骤4所述的运动矢量的矢量差，以矢量差的模最小为准则将对应的主运动矢量赋值给该像素，从而重新确定云范围内的运动矢量，得到云的初步运动矢量，如图4(d)；

步骤5所述分别计算像素点运动矢量与步骤4所述的运动矢量的矢量差的模为：

其中，

和

分别为第1幅和第2幅地基全天空图像的总云范围内第i行第j列像素的块匹配运动矢量的两个分量，

和

分别为第1幅全天空图像第i″个主运动矢量Ω的两个分量。

步骤5中所述以矢量差的模最小为准则将对应的主运动矢量赋值给该像素为：

其中，

为第1幅全天空图像第i行第j列像素的初步运动矢量，(i,j)∈R_1,2，R_1,2为步骤2所述第1幅和第2幅预处理后地基全天空图像的总云范围。

步骤6，根据云的初步运动矢量构造主运动项，连同亮度恒定项、全局平滑项、亮度梯度项和加权中值滤波构建光流法能量函数，对能量函数最小化进行优化求解，得到全天空图像中每个像素点的运动矢量，如图5，图中还展示了根据计算求得的运动矢量预测出的下一时刻的全天空图像I₃以及其与真实全天空图像的差值；

步骤6所述亮度恒定项为：

光流法的基本约束方程，假设同一像素点在连续两幅全天空图像I₁和I₂中的亮度值不发生变化，即灰度值不变；

E_intensity(u,v)＝∫∫ψ(|I₂(x+dx,y+dy)-I₁(x,y)|²)dxdy

其中，I₁(x,y)为某像素点(x,y)在第1幅图像中的亮度值，经过一定时间，该像素点移动了(dx,dy)，此时的亮度值为I₂(x+dx,y+dy)。式中

ε＝0.001为Charbonnier惩罚函数，(u,v)表示计算得到的该像素点(x,y)的运动矢量。

步骤6所述全局平滑项为：

整个图像I₁上光流的变化是平滑的，即u和v随着像素点移动而发生的变化是缓慢的；

步骤6所述亮度梯度项为：

由于太阳周围云的运动可能造成图像晴空背景亮度发生较大变化，使得像素亮度值的匹配出现错误，可通过亮度值梯度的约束来增强亮度恒定项的效果，达到进一步精确地预测出后一幅云图I_k+2的目的；

步骤6所述主运动项为：

将块匹配法获得的云初步运动矢量

作为约束项代入光流法，使得光流计算的迭代次数可适当降低，提高计算效率，同时减弱各云块内部亮度不均一的影响；

步骤6所述加权中值滤波为：

对运动矢量(u,v)进行加权中值滤波能有效改善云边缘处异常值的问题；

首先利用Sobel算子检测云边缘并适当扩大范围，对于非边缘的区域，以[5,5]的窗口大小进行传统中值滤波，而在云边缘区域，综合考虑空间距离、Lab颜色空间差距以及遮挡状态这三个因素，在[15,15]的窗口范围内进行加权中值滤波；

其中，N_x,y是以(x,y)为中心的窗口范围内除中心外的点(x′,y′)的集合，

代表两个像素(x,y)和(x′,y′)之间相似性的权重，具体为下式，I代表Lab空间中的颜色向量，n_c为颜色波段数，σ₁＝7,σ₂＝7：

遮挡状态变量o(x,y)被定义为

式中d(x,y)为单边散度函数，div(x,y)代表光流的散度

步骤6所述的光流法能量函数为：

式中，α,γ,β分别为全局平滑项的权重、亮度梯度项的权重、主运动项的权重；

以光流法能量函数最小化为优化目标，利用欧拉-拉格朗日方程、逐次超松弛法等数学手段求解，得到全天空图像I₁每个像素点的运动矢量(u,v)；

步骤7，针对基于微分的光流算法只能计算小位移的光流问题，引入金字塔思想进行由粗到细的多尺度计算，如图6分别展示了光流法计算中不使用滤波、使用传统中值滤波和使用加权中值滤波的效果。

以0.98的比例将预处理后全天空图像降采样为91个不同分辨率的图像，在金字塔中从上到下按最低分辨率到最高分辨率依次排列，从金字塔顶层开始，对第n层图像进行步骤6的光流法计算，得到运动矢量(u,v)ⁿ，然后作为第n+1层的矢量计算初始值(u₀,v₀)ⁿ⁺¹，直到最底层第91层，从而将大位移分解为若干个小位移的集合，达到计算大位移光流的目的；

步骤7所述一定比例为：0.98；

步骤7所述91个不同分辨率的图像为：

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (7)

1.一种块匹配和光流法结合的地基云图运动预测方法，其特征在于，包括以下步骤：

步骤1，读取原始地基全天空图像，首先利用平面结构引导的修复算法对原始地基全天空图像中被遮挡的区域进行图像修复，再以图像中心为原点、设置一定半径保留图像天空区域、裁掉与天空无关区域，最后根据角度的对应关系进一步进行图像畸变校正，得到第k幅和第k+1幅预处理后地基全天空图像I_k,I_k+1,k∈[1,K-1]，K为预处理后地基全天空图像的数量；

步骤2，通过基于多特征的随机森林分类方法将预处理后地基全天空图像的每个像素分为云类别或者晴空类别，从而将预处理后地基全天空图像转换成云或晴空的二值图像，利用逆四叉树法在二值图像I_bw中确定云范围；

步骤3，根据预处理后地基全天空图像中方块坐标的像素亮度值构建均方误差优化目标，以均方误差最小化为优化目标进行块匹配，得到预处理后地基全天空图像的运动矢量；

步骤4，结合步骤2所述的预处理后地基全天空图像的云范围、步骤3所述的块匹配运动矢量得到预处理后地基全天空图像的云范围内的块匹配运动矢量，筛选出频次较高的多个块匹配运动矢量作为主运动矢量；

步骤5，针对步骤2中所述的相邻两幅预处理后地基全天空图像的总云范围内每个像素，分别计算像素点运动矢量与步骤4所述的运动矢量的矢量差，以矢量差的模最小为准则将对应的主运动矢量赋值给该像素，从而重新确定云范围内的运动矢量，得到云的初步运动矢量；

步骤6，根据云的初步运动矢量构造主运动项，连同亮度恒定项、全局平滑项、亮度梯度项和加权中值滤波构建光流法能量函数，对能量函数最小化进行优化求解，得到全天空图像中每个像素点的运动矢量；

步骤7，针对基于微分的光流算法只能计算小位移的光流问题，引入金字塔思想进行由粗到细的多尺度计算，以一定比例将预处理后全天空图像降采样为N个不同分辨率的图像，在金字塔中从上到下按最低分辨率到最高分辨率依次排列，从金字塔顶层开始，对第n层图像进行步骤6的光流法计算，得到运动矢量(u,v)ⁿ，然后作为第n+1层的矢量计算初始值(u₀,v₀)ⁿ⁺¹，直到最底层N，从而将大位移分解为若干个小位移的集合，达到计算大位移光流的目的。

2.根据权利要求1所述的块匹配和光流法结合的地基云图运动预测方法，其特征在于：

步骤2所述预处理后地基全天空图像的每个像素为：

I_k(i,j)

其中，I_k(i,j)表示第k幅预处理后地基全天空图像中第i行第j列的像素，i∈[1,M],j∈[1,N],k∈[1,K-1]，K为预处理后地基全天空图像的数量，M*N为预处理后地基全天空图像的分辨率；

步骤2所述晴空类别定义为C_sky；

步骤2所述云类别定义为C_cloud；

步骤2所述二值图像，具体定义为：

I_k,bw(i,j),type_k(i,j)

其中，I_k,bw(i,j)表示第k幅预处理后地基全天空图像对应的二值图像I_bw中第i行第j列的像素，i∈[1,M],j∈[1,N],k∈[1,K-1]，K为预处理后地基全天空图像的数量，M*N为预处理后地基全天空图像的分辨率，type_k(i,j)表示第k幅预处理后地基全天空图像对应的二值图像中第i行第j列的像素类别，若I_k,bw(i,j)为晴空类别则type_k(i,j)＝C_sky，若I_k,bw(i,j)为云类别则type_k(i,j)＝C_cloud；

步骤2所述利用逆四叉树法在二值图像中确定相邻两幅预处理后地基全天空图像的总云范围，具体为：

对二值图像进行初步处理：

以n*n的像素方阵作为基本单元即U_k,1(i′,j′)；

U_k,1(i′,j′)为第k幅预处理后地基全天空图像对应的二值图像中第i′行第j′列的第1级方阵，即基本单元，方阵的尺寸为n*n，i′∈[1,M/n],j′∈[1,N/n],k∈[1,K-1]，K为预处理后地基全天空图像的数量；

当U_k,1(i′,j′)内属于云类别的像素个数超过单元总像素数n*n的一定比例时，则将U_k,1(i′,j′)内所有像素设置为云类别C_cloud，否则全部设置为晴空类别C_sky；

将基本单元按照四叉树原理自下而上逐步聚拢，得到第k幅全天空图像不同尺寸的云块单元U_k,1,U_k,2,...,U_k,m,...,U_k,L，其中第m级的方阵尺寸为2^m-1n*2^m-1n；

第k幅预处理后地基全天空图像的云范围R_k为第k幅预处理后地基全天空图像对应的二值图像中所有云块单元U_k,1,U_k,2,...,U_k,m,...,U_k,L的集合；

对第k+1幅预处理后地基全天空图像I_k+1重复上述操作得到对应的云范围R_k+1，则第k幅和第k+1幅预处理后地基全天空图像的总云范围为：

R_k,k+1＝R_k∪R_k+1，R_k,k+1为R_k和R_k+1的并集。

3.根据权利要求1所述的块匹配和光流法结合的地基云图运动预测方法，其特征在于：

步骤3所述构建均方误差优化目标为：

其中，

为第k幅全天空图像中第i行第j列像素的运动矢量，i∈[1,M],j∈[1,N],k∈[1,K-1]，K为预处理后地基全天空图像的数量，M*N为预处理后地基全天空图像的分辨率；方块的大小为N′*N′，f_k+1(n+u,n+v)表示第k+1幅图像中方块坐标(n+u,n+v)处的像素亮度值，f_k(n,n)则表示第k幅参考图像搜索窗口中方块坐标为(n,n)处的像素亮度值，u和v分别为横向和纵向的运动距离，u∈[-a,a],v∈[-a,a]，a为块匹配搜索窗口的大小。

4.根据权利要求1所述的块匹配和光流法结合的地基云图运动预测方法，其特征在于：

步骤4所述相邻两幅预处理后地基全天空图像的总云范围内像素的块匹配运动矢量为：

其中，第k幅全天空图像第i行第j列像素为云类别C_cloud，其运动矢量为

R_k,k+1为步骤2所述第k幅和第k+1幅预处理后地基全天空图像的总云范围，M*N为预处理后地基全天空图像的分辨率，K为预处理后地基全天空图像的数量；

步骤4所述筛选出频次较高的多个矢量作为主运动矢量为：

统计第k幅全天空图像中每对

出现的频次

选取频次

最高的多个运动矢量作为该图像云层的主运动矢量Ω：

num为主运动矢量的个数。

5.根据权利要求1所述的块匹配和光流法结合的地基云图运动预测方法，其特征在于：

步骤5所述分别计算像素点运动矢量与步骤4所述的运动矢量的矢量差的模为：

其中，

和

分别为第k幅和第k+1幅地基全天空图像的总云范围内第i行第j列像素的块匹配运动矢量的两个分量，

和

分别为第k幅全天空图像第i″个主运动矢量Ω的两个分量；

步骤5中所述以矢量差的模最小为准则将对应的主运动矢量赋值给该像素为：

其中，

为第k幅全天空图像第i行第j列像素的初步运动矢量，(i,j)∈R_k,k+1,k∈[1,K-1]，R_k,k+1为步骤2所述第k幅和第k+1幅预处理后地基全天空图像的总云范围，K为预处理后地基全天空图像的数量。

6.根据权利要求1所述的块匹配和光流法结合的地基云图运动预测方法，其特征在于：

步骤6所述亮度恒定项为：

光流法的基本约束方程，假设同一像素点在连续两幅全天空图像I_k和I_k+1中的亮度值不发生变化，即灰度值不变，其中k∈[1,K-1]，K为地基全天空图像的数量；

E_intensity(u,v)＝∫∫ψ(|I_k+1(x+dx,y+dy)-I_k(x,y)|²)dxdy

其中，I_k(x,y)为某像素点(x,y)在第k幅图像中的亮度值，经过一定时间，该像素点移动了(dx,dy)，此时的亮度值为I_k+1(x+dx,y+dy)；式中

ε＝0.001为Charbonnier惩罚函数，(u,v)表示计算得到的该像素点(x,y)的运动矢量；

步骤6所述全局平滑项为：

整个图像I_k上光流的变化是平滑的，即u和v随着像素点移动而发生的变化是缓慢的；

步骤6所述亮度梯度项为：

由于太阳周围云的运动可能造成图像晴空背景亮度发生较大变化，使得像素亮度值的匹配出现错误，可通过亮度值梯度的约束来增强亮度恒定项的效果，达到进一步精确地预测出后一幅云图I_k+2的目的；

步骤6所述主运动项为：

将块匹配法获得的云初步运动矢量

作为约束项代入光流法，使得光流计算的迭代次数可适当降低，提高计算效率，同时减弱各云块内部亮度不均一的影响；

步骤6所述加权中值滤波为：

对运动矢量(u,v)进行加权中值滤波能有效改善云边缘处异常值的问题；

首先利用Sobel算子检测云边缘并适当扩大范围，对于非边缘的区域，以[m,m]的窗口大小进行传统中值滤波，而在云边缘区域，综合考虑空间距离、Lab颜色空间差距以及遮挡状态这三个因素，在[M,M]的窗口范围内进行加权中值滤波，且M＞m；

其中，N_x,y是以(x,y)为中心的窗口范围内除中心外的点(x′,y′)的集合，

代表两个像素(x,y)和(x′,y′)之间相似性的权重，具体为下式，I代表Lab空间中的颜色向量，n_c为颜色波段数，σ₁＝7,σ₂＝7：

遮挡状态变量o(x,y)被定义为

式中d(x,y)为单边散度函数，div(x,y)代表光流的散度

步骤6所述的光流法能量函数为：

式中，α,γ,β分别为全局平滑项的权重、亮度梯度项的权重、主运动项的权重；

以光流法能量函数最小化为优化目标，利用欧拉-拉格朗日方程、逐次超松弛法等数学手段求解，得到全天空图像I_k每个像素点的运动矢量(u,v)。

7.根据权利要求1所述的块匹配和光流法结合的地基云图运动预测方法，其特征在于：

步骤7所述一定比例为：r,0＜r＜1；

步骤7所述N个不同分辨率的图像为：

Images