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

技术领域

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

背景技术

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

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

块匹配方法首先将图像划分为若干块，假定每一个块内的所有像素作为一个整体，具有相同的运动模式，围绕某块B

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

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

发明内容

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

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

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

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

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

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

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

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

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

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

I

其中，I

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

步骤2所述云类别定义为C

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

I

其中，I

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

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

以n*n的像素方阵作为基本单元即U

U

当U

将基本单元按照四叉树原理自下而上逐步聚拢，得到第k幅全天空图像不同尺寸的云块单元U

第k幅预处理后地基全天空图像的云范围R

对第k+1幅预处理后地基全天空图像I

R

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

其中，

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

其中，第k幅全天空图像第i行第j列像素为云类别C

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

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

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

其中，

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

其中，

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

光流法的基本约束方程，假设同一像素点在连续两幅全天空图像I

E

其中，I

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

整个图像I

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

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

步骤6所述主运动项为：

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

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

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

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

其中，N

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

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

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

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

以光流法能量函数最小化为优化目标，利用欧拉-拉格朗日方程、逐次超松弛法等数学手段求解，得到全天空图像I

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

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

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

附图说明

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

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

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

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

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

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

具体实施方式

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

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

实施路线如图1所示：

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

得到的第1幅和第2幅预处理后地基全天空图像I

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

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

I

其中，I

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

步骤2所述云类别定义为C

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

I

其中，I

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

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

以5*5的像素方阵作为基本单元即U

U

当U

将基本单元按照四叉树原理自下而上逐步聚拢，得到第1幅全天空图像不同尺寸的云块单元U

第1幅预处理后地基全天空图像的云范围R

对第2幅预处理后地基全天空图像I

R

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

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

其中，

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

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

其中，第1幅全天空图像第i行第j列像素为云类别C

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

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

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

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

其中，

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

其中，

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

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

光流法的基本约束方程，假设同一像素点在连续两幅全天空图像I

E

其中，I

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

整个图像I

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

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

步骤6所述主运动项为：

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

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

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

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

其中，N

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

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

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

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

以光流法能量函数最小化为优化目标，利用欧拉-拉格朗日方程、逐次超松弛法等数学手段求解，得到全天空图像I

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

以0.98的比例将预处理后全天空图像降采样为91个不同分辨率的图像，在金字塔中从上到下按最低分辨率到最高分辨率依次排列，从金字塔顶层开始，对第n层图像进行步骤6的光流法计算，得到运动矢量(u,v)

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

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

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

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