CN106683114A - 基于特征光流的流体运动矢量估计方法 - Google Patents

Abstract

本发明公开一种基于特征光流的流体运动矢量估计方法，包括如下步骤：首先对H.264视频码流中的连续两帧图像I₁和I₂分别采用多尺度SUSAN特征点检测算法实现不同分辨率下特征点的提取及亚像素定位，获得特征点集P₁和P₂；然后对P₁中的特征点采用H.264视频压缩域先验的特征光流估计算法提取特征点所在宏块的运动矢量作为设置搜索半径的先验信息，在特征点集P₂中搜索相匹配的特征点，进而估计光流运动矢量，获得稀疏的特征光流场；最后对稀疏的特征光流场采用基于反距离权重插值的流场网格化算法进行网格化处理，得到均匀的二维流速矢量场。本发明能够用于河流水面等天然示踪水流的运动矢量估计，特别适用于高空间分辨率的二维瞬时流场测量。

Description

基于特征光流的流体运动矢量估计方法

技术领域

本发明涉及一种流体运动矢量估计方法，尤其涉及一种基于特征光流的流体运动矢量估计方法，属于流场测量技术领域。

背景技术

流体中存在着大量不具有固定质量并存在相对运动的流体微团，使得流体的运动形态比刚体和固体更为复杂。在过去的二十多年间，光学、电子及计算机技术的发展促进了粒子图像测速(PIV)这种非接触式瞬时全场流速测量技术的实现和应用。它通过粒子图像序列的分析和计算获得局部流体运动位移及速度的大小、方向、特征和分布情况，极大提高了实验室环境下各种复杂流动的测量能力。其中从粒子图像估计流体的运动矢量是PIV技术的核心和难点，方法的选择不仅依赖于硬件系统而且取决于待测流体的特点。

对于利用树叶、树枝等天然漂浮物以及漩滚、表面波等天然水面模式作为水流示踪物的大尺度粒子图像测速(LSPIV)应用而言，河流水面流场的运动矢量估计又面临着特有的挑战。首先，大多数天然水面模式产生于明渠紊流，它们没有明确的粒子特性，通常存在振荡、扩散及形变等非刚性运动，因而缺乏稳定的几何特征和统计特征，使得经典的粒子跟踪测速(PTV)方法难以直接适用。其次，现场条件下的待测流场面积往往覆盖数百至数千平方米，为满足流速测量的目标精度，图像分辨率需要达到百万至千万像素，远大于目前实验室中十万像素的数量级，使得基于直接互相关(DCC)等空域相关匹配的运动矢量估计方法计算量相当可观，难以满足流场、流量连续监测的实时性要求。再则，河流水面的流速分布在断面方向表现出从中泓区域向两岸递减的趋势，流速梯度加之图像透视畸变和水流示踪物分布不均的影响往往导致水面瞬时流场具有复杂的区域特征。采用固定大小的观测窗口不仅会降低测量的空间分辨率，而且会因为对局部流场的适应性较差而降低运动矢量估计的可靠性。

基于图像序列的运动目标跟踪技术是图像处理、模式识别和计算机视觉领域的重要研究课题,在国防、工业和交通等领域有着广泛的应用前景。其中光流法利用图像中灰度模式的时域变化和相关性来确定像素点的二维瞬时速度，进而分析场景中目标的结构和运动。通常可以分为微分光流法和特征光流法。前者基于光流的微分基本方程，通过附加一定的约束条件，求得图像上任一点的运动速度，获得一个稠密的光流场；后者通过寻求序列图像间特征点的对应关系，依据相似性准则进行特征匹配，进而获得一个稀疏的光流场。特征点是图像上具有较好的可区分性、位置不变性和特征稳定性的点。相比LSPIV中基于图像灰度特征(PIV)、粒子形态特征或粒子运动特征(PTV)的方法，基于特征光流的方法具有以下优势：①可用于表征不同类型的示踪物，易于提取和分析。②特征相对稳定，对图像噪声、灰度变化、局部形变以及遮挡等具有较好的适应能力。③位置的度量值对局部灰度的变化比较敏感，可提高定位的精确度。④数目比图像像素点或其他特征要少的多，可以显著减少匹配的计算量。然而特征点在计算机视觉领域中主要被用于人物、车辆等大型目标的识别和分类，对于小目标的运动矢量估计，在检测和定位的精度上尚存在不足。例如，基于梯度的Harris算子提取图像灰度梯度的局部极大值点作为角点，具有旋转不变性和较高的稳定性；但抗噪能力较差，且角点的位置可能偏移或存在伪角点。基于灰度的SUSAN算子能够在较强噪声中提取出目标的角点及边缘，定位精度达到像素级；但不具有尺度不变性。基于多尺度思想的SIFT算子能够提取图像中具有尺度、旋转和仿射不变性的局部特征点；但往往并不对应于目标上的真实角点，由于需要从差分高斯(DoG)金字塔中提取128维的特征向量，运算量也很可观。因此研究针对河流水面这类特殊流体的特征光流运动矢量估计方法具有重要的理论意义和应用价值。

发明内容

发明目的：针对现有技术中存在的不足，本发明提供了一种基于特征光流的流体运动矢量估计方法。

技术方案：本发明提出一种基于特征光流的流体运动矢量估计方法，包括如下步骤：

1)对H.264视频码流中的连续两帧图像I₁和I₂分别采用多尺度SUSAN特征点检测算法实现不同分辨率下特征点的提取及亚像素定位，获得特征点集P₁和P₂，具体步骤如下：

1.1)建立图像的高斯金字塔尺度空间：高斯核被证明是实现尺度变换唯一的变换核和线性核，因此利用高斯变换的多尺度特性构建高斯金字塔，获得不同尺度下的分层图像；数字图像I(x,y)的尺度空间L(x,y,σ_i)可以定义为一个变化尺度的高斯函数G(x,y,σ_i)与其自身的卷积，即：

其中，(x,y)表示像素点的位置，m和n表示高斯模板的大小，σ_i是尺度因子，i＝1,…,level，其值越小表示图像被平滑的越少；在小尺度下定位准确，大尺度下抗噪能力强；

1.2)检测尺度空间中的SUSAN特征点：采用SUSAN算子检测尺度空间中的分层图像，获得不同尺度下的候选角点；SUSAN算子采用一种圆形模板在图像上移动，模板内部每个图像像素点的灰度值都和模板中心像素的灰度值作比较，若模板内某个像素的灰度与核的灰度差值小于一定值，则认为该点与核具有相同或相近的灰度，由满足这一条件的像素组成的区域称为吸收核同值区USAN；位于图像上局部最小的USAN区的核像素就是图像的特征角点；图像中某像素的USAN区大小可由下式表示：

n(r₀)＝∑_rC(r₀,r) (2)

其中，C(r₀)是以r₀为圆心的模板；C(r₀,r)是模板内属于USAN区域的像素的判别函数，定义为：

其中，t是亮度差门限，表示所检测的特征的最小对比度，也是噪声的最大容限，决定了提取特征量的多少；用模板扫描整个图像得到每个像素的USAN区域后，再由下式产生USAN特征图像：

式中，g为几何门限，一般取n_max/2，n_max是n(r)所能达到的最大值，g的大小决定了所得特征角点的尖锐程度；

1.3)筛选不同尺度下的均匀化的候选特征点：由于示踪物的密度低、时空分布不均，会造成真实河流水面图像中特征点的聚簇现象，采用图像分块和邻近点剔除的思想对有效特征点进行筛选；首先采用固定大小的均匀网格对图像进行分块，网格大小参考流速场测量的空间分辨率进行选取；然后将子块中检测出的角点按R值大小进行排序，筛选出R值相对较大的角点，以提高特征点检测的重复性；最后选用一个滑动模板对整幅图像进行处理，若在该模板下存在不止一个角点，则只保留R值最大的角点；至此获得具有局部最大特征强度的均匀化候选特征点集；

1.4)定位候选特征点的亚像素坐标：由于单一尺度下SUSAN算子检测得到的角点坐标是像素级的，利用同一角点在多尺度下的冗余信息实现坐标精定位；首先把不同尺度下的候选特征点还原到原始图像中的相应位置构成候选特征点簇；然后依次扫描各候选特征点，对其邻域内特征点簇的坐标求算术平均作为精化的亚像素坐标，获得图像的特征点集，即具有亚像素坐标的均匀化候选特征点集。

2)对P₁中的特征点采用H.264视频压缩域先验的特征光流估计算法提取特征点所在宏块的运动矢量作为设置搜索半径的先验信息，在特征点集P₂中搜索相匹配的特征点，进而估计光流运动矢量，获得稀疏的特征光流场，具体步骤如下：

2.1)确定匹配空间：首先在图像I₁中以坐标为(x,y)的特征点P_i为中心建立模板图像f(x,y)；然后对H.264视频码流依次进行熵解码、重排序、反量化和反变换分析，定位到宏块层句法元素并提取出特征点P_i所在宏块的运动矢量差值Δu_i；接下来从已解码的参考序列中找到相应宏块的运动矢量和差值相加后得到了当前位置宏块的真实运动矢量大小最后在图像I₂中以(x,y)为中心、2U_i为半径建立匹配空间；

2.2)相似性度量：首先在上述图像I₂中的匹配空间内依次选取一个坐标为(m,n)的特征点P_k建立一个和f(x,y)大小相同的待匹配图像g(m,n)；然后将二维实矩阵扩充成M×N的复数型矩阵，即每个像素由整型的实部和整型的虚部组成，其中虚部为0；接下来分别进行两次二维快速傅立叶变换得到f(x,y)和g(m,n)的频域形式：

最后利用F(p,q)的共轭F^*(p,q)计算f(x,y)和g(m,n)的频域互相关：

C(p,q)＝F^*(p,q)·G(p,q) (6)

2.3)估计光流运动矢量：重复步骤2.1和2.2得到匹配空间内所有特征点的频域互相关，并选取其中的最大者作为匹配特征点，坐标为(a,b)；根据两个特征点间的距离及图像I₁和I₂间的时间间隔Δt计算特征点P_i对应的光流运动矢量

2.4)重复步骤2.1-2.3，对图像I₁中所有的特征点计算运动矢量，获得一幅稀疏的特征光流场。

3)对稀疏的特征光流场采用基于反距离权重插值的流场网格化算法进行网格化处理，得到均匀的二维流速矢量场，具体步骤如下：

3.1)在图像I₁中划分均匀的矩形测速网格；

3.2)采用以下反距离权重法估计网格节点j上的运动矢量：

其中，d_i为特征点i与网格节点j之间的距离，z为网格内特征点的数量，w为权重系数；

3.3)重复步骤3.2对测速网格中所有节点的运动矢量进行插值可获得网格化的流速矢量场。

有益效果：本发明能够用于河流水面等天然示踪水流的运动矢量估计，特别适用于高空间分辨率的二维瞬时流场测量，具体优点如下：

1、鲁棒性强。采用图像中的尺度不变特征点替代传统意义上的水流示踪物，可以在一定程度上克服示踪物非刚性运动引起的不确定性，而频域相关匹配对噪声有较高的容忍度，可以处理图像之间的旋转和尺度变化，因此方法在野外河流的应用中具有较好的鲁棒性。

2、分辨率高。特征点的定位精度可达亚像素级，相比空域相关匹配法中的矩形窗口，有效改善了相机在倾斜视角下拍摄引起的图像远场空间分辨率降低，使测量结果更具有代表性。根据流场运动特征，采用反距离权重插值法对稀疏的特征光流场进行网格化处理，也能够获得均匀分布的流场数据。

3、实时性强。频域相关匹配将空域卷积转换为频域乘积，可显著降低运算量；利用H.264视频编解码的特点，从压缩域码流中提取宏块的运动矢量作为确定匹配空间的先验信息，大大缩小了搜索空间，提高了算法的实时性。

附图说明

图1是本发明涉及的运动矢量估计方法的工作流程图；

图2是本发明涉及的多尺度SUSAN特征点检测算法的示意图；

图3是本发明涉及的基于反距离权重插值的流场网格化算法的示意图。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

本发明的运动矢量估计方法的工作流程如图1所示。首先采用一种多尺度SUSAN特征点检测算法检测图像中的特征点，算法基于尺度空间理论，在传统SUSAN算子的基础上引入连续变化的尺度参数，获得多尺度下的尺度空间表示序列，进而实现不同分辨率下特征点的提取及亚像素定位，提高了特征点检测的鲁棒性和精确性；然后利用视频编解码的特点采用一种H.264视频压缩域先验的特征光流估计算法，从压缩域码流中提取宏块的运动矢量作为设置搜索半径的先验信息以提高频域相关匹配的效率；由于特征点分布不均，特征光流法输出的流场结果为非结构化的离散数据，为了便于观察流体运动的规律，最后采用基于反距离权重插值的流场网格化算法对稀疏的特征光流场进行网格化处理，得到均匀的二维流速矢量场。

本发明的多尺度SUSAN特征点检测算法如图2所示。主要步骤描述如下：

1)建立图像的高斯金字塔尺度空间：高斯核被证明是实现尺度变换唯一的变换核和线性核，因此利用高斯变换的多尺度特性构建高斯金字塔，获得不同尺度下的分层图像；数字图像I(x,y)的尺度空间L(x,y,σ_i)可以定义为一个变化尺度的高斯函数G(x,y,σ_i)与其自身的卷积，即：

其中，(x,y)表示像素点的位置，m和n表示高斯模板的大小，σ_i是尺度因子，i＝1,2,3,…，其值越小表示图像被平滑的越少；在小尺度下定位准确，大尺度下抗噪能力强；

2)检测尺度空间中的SUSAN特征点：采用SUSAN算子检测尺度空间中的分层图像，获得不同尺度下的候选角点；SUSAN算子采用一种圆形模板在图像上移动，模板内部每个图像像素点的灰度值都和模板中心像素的灰度值作比较，若模板内某个像素的灰度与核的灰度差值小于一定值，则认为该点与核具有相同或相近的灰度，由满足这一条件的像素组成的区域称为吸收核同值区USAN。位于图像上局部最小的USAN区的核像素就是图像的特征角点；图像中某像素的USAN区大小可由下式表示：

n(r₀)＝∑_rC(r₀,r) (2)

式中C(r₀)是以r₀为圆心的模板；C(r₀,r)是模板内属于USAN区域的像素的判别函数，定义为：

式中，t是亮度差门限，表示所检测的特征的最小对比度，也是噪声的最大容限，决定了提取特征量的多少；用模板扫描整个图像得到每个像素的USAN区域后，再由下式产生USAN特征图像：

式中，g为几何门限，一般取n_max/2，n_max是n(r)所能达到的最大值。g的大小决定了所得特征角点的尖锐程度；

3)筛选均匀化的候选特征点：由于示踪物的密度低、时空分布不均，会造成真实河流水面图像中特征点的聚簇现象，采用图像分块和邻近点剔除的思想对有效特征点进行筛选；首先采用固定大小的均匀网格对图像进行分块，网格大小参考流速场测量的空间分辨率进行选取；然后将子块中检测出的角点按R值大小进行排序，筛选出R值相对较大的角点，以提高特征点检测的重复性；最后选用一个滑动模板对整幅图像进行处理，若在该模板下存在不止一个角点，则只保留R值最大的角点；至此获得具有局部最大特征强度的均匀化候选特征点集。

4)定位候选特征点的亚像素坐标：由于单一尺度下SUSAN算子检测得到的角点坐标是像素级的，利用同一角点在多尺度下的冗余信息实现坐标精定位；首先把不同尺度下的候选特征点还原到原始图像中的相应位置构成候选特征点簇；然后依次扫描各候选特征点，对其邻域内特征点簇的坐标求算术平均，作为精化的亚像素坐标。

本发明的H.264视频压缩域先验的特征光流估计算法主要步骤描述如下：

1)确定匹配空间：首先在图像I₁中以坐标为(x,y)的特征点P_i为中心建立模板图像f(x,y)；然后对H.264视频码流依次进行熵解码、重排序、反量化和反变换分析，定位到宏块层句法元素并提取出特征点P_i所在宏块的运动矢量差值Δu_i；接下来从已解码的参考序列中找到相应宏块的运动矢量和差值相加后得到了当前位置宏块的真实运动矢量大小最后在图像I₂中以(x,y)为中心、2U_i为半径建立匹配空间；

2)相似性度量：首先在上述图像I₂中的匹配空间内依次选取一个坐标为(m,n)的特征点P_k建立一个和f(x,y)大小相同的待匹配图像g(m,n)；然后将二维实矩阵扩充成M×N的复数型矩阵，即每个像素由整型的实部和整型的虚部组成，其中虚部为0；接下来分别进行两次二维快速傅立叶变换得到f(x,y)和g(m,n)的频域形式：

最后利用F(p,q)的共轭F^*(p,q)计算f(x,y)和g(m,n)的频域互相关：

C(p,q)＝F^*(p,q)·G(p,q) (6)

3)估计光流运动矢量：重复步骤1和2得到匹配空间内所有特征点的频域互相关，并选取其中的最大者作为匹配特征点，坐标为(a,b)；根据两个特征点间的距离及图像I₁和I₂间的时间间隔Δt计算特征点P_i对应的光流运动矢量

本发明的基于反距离权重插值的流场网格化算法如图3所示。反距离权重插值基于一种假设，即彼此距离较近的事物比彼此距离较远的事物更相似，并且假定每个测量点都有一种局部影响，这种影响会随着距离的增大而减小。这与流体运动连续性的定律相吻合。主要步骤描述如下：

1)采用特征光流估计算法对图像I₁中所有的离散特征点(x_i,y_i)计算运动矢量获得一幅稀疏的特征光流场；

2)在流场图像中划分均匀的矩形测速网格；

3)采用以下反距离权重法估计网格节点j上的运动矢量：

其中，d_i为特征点i与网格节点j之间的距离；z为网格内特征点的数量；w为权重系数，w越大表明距离相关性越明显；

4)重复步骤3对测速网格中所有节点的运动矢量进行插值可获得网格化的流速矢量场。

Claims (4)

1.一种基于特征光流的流体运动矢量估计方法，其特征在于，包括如下步骤：

1)对H.264视频码流中的连续两帧图像I₁和I₂分别采用多尺度SUSAN特征点检测算法实现不同分辨率下特征点的提取及亚像素定位，获得特征点集P₁和P₂；

2)对P₁中的特征点采用H.264视频压缩域先验的特征光流估计算法提取特征点所在宏块的运动矢量作为设置搜索半径的先验信息，在特征点集P₂中搜索相匹配的特征点，进而估计光流运动矢量，获得稀疏的特征光流场；

3)对稀疏的特征光流场采用基于反距离权重插值的流场网格化算法进行网格化处理，得到均匀的二维流速矢量场。

2.根据权利要求1所述的基于特征光流的流体运动矢量估计方法，其特征在于：所述步骤1采用多尺度SUSAN特征点检测算法提取特征点集的具体步骤如下：

1.1)建立图像的高斯金字塔尺度空间：利用高斯变换的多尺度特性构建高斯金字塔，获得不同尺度下的分层图像；

1.2)检测尺度空间中的SUSAN特征点：采用SUSAN算子检测尺度空间中的分层图像，获得不同尺度下的候选角点；

1.3)筛选不同尺度下的均匀化的候选特征点：首先采用固定大小的均匀网格对图像进行分块，网格大小参考流速场测量的空间分辨率进行选取；然后将子块中检测出的角点按R值大小进行排序，筛选出R值相对较大的角点，以提高特征点检测的重复性；最后选用一个滑动模板对整幅图像进行处理，若在该模板下存在不止一个角点，则只保留R值最大的角点，至此获得具有局部最大特征强度的均匀化候选特征点集；

1.4)定位候选特征点的亚像素坐标：首先把不同尺度下的候选特征点还原到原始图像中的相应位置构成候选特征点簇；然后依次扫描各候选特征点，对其邻域内特征点簇的坐标求算术平均作为精化的亚像素坐标，获得图像的特征点集。

3.根据权利要求1所述的基于特征光流的流体运动矢量估计方法，其特征在于：所述步骤2采用H.264视频压缩域先验的特征光流估计算法获取稀疏的特征光流场的具体步骤如下：

2.1)确定匹配空间：首先在图像I₁中以坐标为(x,y)的特征点P_i为中心建立模板图像f(x,y)；然后对H.264视频码流依次进行熵解码、重排序、反量化和反变换分析，定位到宏块层句法元素并提取出特征点P_i所在宏块的运动矢量差值Δu_i；接下来从已解码的参考序列中找到相应宏块的运动矢量和差值相加后得到了当前位置宏块的真实运动矢量大小最后在图像I₂中以(x,y)为中心、2U_i为半径建立匹配空间；

2.2)相似性度量：首先在上述图像I₂中的匹配空间内依次选取一个坐标为(m,n)的特征点P_k建立一个和f(x,y)大小相同的待匹配图像g(m,n)；然后将二维实矩阵扩充成M×N的复数型矩阵，即每个像素由整型的实部和整型的虚部组成，其中虚部为0；接下来分别进行两次二维快速傅立叶变换得到f(x,y)和g(m,n)的频域形式：

$\left\{\begin{array}{c}F(p,q)=FFT2\left(f\left(x,y\right)\right)=\frac{1}{MN}\underset{x=0}{\overset{M-1}{\Σ}}\underset{y=0}{\overset{N-1}{\Σ}}f(x,y)e^{-j2\mathrm{\π}(\frac{xp}{M}+\frac{yq}{N})}\\ G(p,q)=FFT2\left(g\left(m,n\right)\right)=\frac{1}{MN}\underset{m=0}{\overset{M-1}{\Σ}}\underset{n=0}{\overset{N-1}{\Σ}}g(m,n)e^{-j2\mathrm{\π}(\frac{mp}{M}+\frac{nq}{N})}\end{array}\right.$

最后利用F(p,q)的共轭F^*(p,q)计算f(x,y)和g(m,n)的频域互相关：

C(p,q)＝F^*(p,q)·G(p,q)

2.3)估计光流运动矢量：重复步骤2.1和2.2得到匹配空间内所有特征点的频域互相关，并选取其中的最大者作为匹配特征点，坐标为(a,b)；根据两个特征点间的距离及图像I₁和I₂间的时间间隔Δt计算特征点P_i对应的光流运动矢量

2.4)重复步骤2.1-2.3，对图像I₁中所有的特征点计算运动矢量，获得一幅稀疏的特征光流场。

4.根据权利要求1所述的基于特征光流的流体运动矢量估计方法，其特征在于：所述步骤3获取均匀的二维流速矢量场的具体步骤如下：

3.1)在图像I₁中划分均匀的矩形测速网格；

3.2)采用以下反距离权重法估计网格节点j上的运动矢量：

其中，d_i为特征点i与网格节点j之间的距离，z为网格内特征点的数量，w为权重系数；

3.3)重复步骤3.2对测速网格中所有节点的运动矢量进行插值可获得网格化的流速矢量场。