CN107958464B - 噪声环境下光流场快速稳健估计方法 - Google Patents

Abstract

本发明涉及一种噪声环境下光流场快速稳健估计方法。所提算法基于噪声环境下光流场估计方法（ML法），引入惩罚因子以增强光流计算稳健性，并在光流计算迭代公式中加入动量因子缩短光流计算收敛时间以加快光流场计算，而后基于变分方法极小化光流能量函数求解欧拉‑拉格朗日方程，最后通过迭代方法求得速度场。仿真结果表明，对视频中连续两帧图片加入不同高斯噪声后，与M算法及ML方法相比，所提算法可显著增强光流场计算稳健性，缩短光流计算收敛时间，加快光流场计算。

Description

噪声环境下光流场快速稳健估计方法

技术领域

本发明属于无人机视频监测领域，具体涉及一种噪声环境下光流场快速稳健估计方法。

背景技术

光流的概念是Gibson于1950年提出的。光流是图像亮度模式的表观运动，表达了图像变化，由于它包含目标运动信息，因此可被观察者用来确定目标运动情况。

光流技术在目标对象分割、识别、跟踪、机器人导航以及形状信息恢复、模式识别等计算机视觉领域有着广泛应用。1981年，Horn和Schunck创造性地将二维速度场与灰度场结合，在光流基本约束方程基础上引入全局平滑约束，得到光流计算的基本方法(HS法)。HS法是在两组基本假设基础上建立起来的，其中，基本约束方程基于灰度恒定假设即在连续相邻图像上，对应物体上同一点的像素灰度值相同；全局平滑约束方程则基于光流平滑假设即光流场处处平滑。

近些年，随着光流技术的广泛应用，众多学者致力于稳健性更好，速度更快的光流计算方法研究。继Horn和Schunck的HS光流算法之后，其他学者相继提出许多新的光流计算方法，如：Lucas-Kanade的局部平滑法(LK法)、金字塔LK光流法、局部与全局结合光流法、小波光流法等，使得光流计算得到进一步发展。然而，对于实际场景而言，上述方法所得光流场都存在较大误差，主要缘由在于灰度守恒假设的不合理性。实际上，多数应用场景中，即使光照条件不变，连续相邻图像上物体同一元面对应像素灰度仍有可能是变化的。

基于以上原因，日本学者MUKAWA对光流基本约束方程进行修正，并根据物体运动、投影及光流模型对修正项表达式进行推导得到新的光流场计算方法(M法)。此方法通过对图像序列进行光流场估计，取得了优于传统方法的效果。然而，此修正依然是在理想图像环境中进行的，即没有考虑实际视频图像中噪声对光流计算的影响，因此该算法对噪声的鲁棒性较差。针对此问题，马龙等提出噪声环境下光流场估计方法(ML法)，即在M法基础上对基本方程做了进一步修正并加入噪声约束。仿真实验表明，ML方法所得光流较M法而言，对噪声有较好的鲁棒性，但所得光流稳健性较差且计算复杂度较高。

发明内容

本发明提出一种噪声环境下光流场快速稳健估计方法，以解决噪声影响下光流计算稳健性较差及收敛速度慢的问题。

本发明为解决其技术问题所采用的技术方案为提供一种噪声环境下光流场快速稳健估计方法，该方法包括以下步骤：

步骤1：光流约束方程

假设图像上一像素点(x,y)在t时刻的亮度为I(x,y,t)，在Δt时间后所述像素点亮度变为I(x+Δx,y+Δy,t+Δt)，当Δt趋于无穷小时可认为该点亮度不变，得到等式：

将

改写为I_x,I_y,I_t，并令

则重写上式，得到光流计算的基本等式：

I_xu+I_yv+I_t＝0 (2)

u,v表示速度场矢量的两个分量；

步骤2：构建能量函数方程

①光流约束方程修正

设物体面元由位置P运动到位置P'，对应成像投影在像平面上由p移动到p'，假设物体表面为朗伯面，则p和p'的灰度分别为：

其中，ρ_d为物体表面漫反射系数，I_q为入射光强度，设为常数；N和N'分别为物体面元在位置P和P'的单位法向量；L为入射光的单位方向向量；θ为N和L的夹角；θ'为N'和L的夹角；

光流基本方程的修正项为：

其中：

ΔN表示曲面法线运动变化；

从而将基本方程修正为：

cI+I_xu+I_yv+I_t＝0 (5)

定义光流约束因子为：

e≈(cI+I_xu+I_yv+I_t)² (6)

其中，c为参量，连同u,v采用全部平滑约束；全部平滑约束用梯度模的平方和来测度；

定义光流场的全局平滑约束因子为：

e₁＝||▽u||²+||▽v||²＝(u_x ²+v_x ²+u_y ²+v_y ²) (7)

定义参量c的全局平滑约束因子为:

e₂＝||▽c||²＝(c_x ²+c_y ²) (8)

设I'(x,y,z)为I(x,y,z)受加性噪声n(x,y,z)污染后图像点(x,y)的实际灰度值，则：

I(x,y,z)＝I(x,y,z)-n(x,y,z) (9)

假设式(9)中的噪声n(x,y,z)是均值为0，方差为δ²的高斯白噪声，则：

I(x,y,z)((x,y)∈Ω)未受噪声污染，故应满足全局平滑约束，定义灰度全局平滑约束因子为：

e₃＝||▽I||²＝(I_x ²+I_y ²) (12)

②引入惩罚项

引入二次项β²(u²+v²+c²)作为惩罚项：

e₄＝β²(u²+v²+c²) (13)

噪声环境下光流场估计问题可转化为式(11)约束下极小化：

其中，α₁,α₂,α₃为调和4类约束因子的权参数；

步骤3：求解光流场

①光流求解迭代公式

根据拉格朗日乘数法，式(11)约束下极小化式(14)即是极小化：

令：

F(x,y,I,u,v,c,I_x,I_y,u_x,u_y,v_x,v_y,c_x,c_y)＝e+α₁e₁+α₂e₂+α₃e₃+e₄+λ(I-I')-λδ²(16)

由变分法，式(15)取极小值的必要条件为：

将式(16)带入式(17)，化简后可得：

其中，s_xx+s_yy＝▽²s(s＝I,u,v,c)为s分量的拉普拉斯算子，可采用下式估算：

和s_i,j,k分别为t_k时刻(x_i,y_i)处s分量的取值及其邻域平均值，K＝1；

定义为：

由式(18)式中(d)可得：

(I_xu+I_yv+I_t+cI)＝α₂(c_xx+c_yy)/I (21)

化简得：

其中：

式(23)两端乘以(I-I')并在图像平面Ω上积分得：

将式(11)带入式(24)可得:

由式(22)和(25)进一步可得：

②添加动量因子后的光流求解迭代公式

在光流迭代公式中添加动量因子以加快光流场计算，即：

其中：

在对图像求x方向、y方向和t方向的导数时，采用的公式如下：

步骤4：迭代求解

根据式(28)，用Gauss-Seidel方法进行迭代，并计算两次迭代光流值的误差，如果小于0.005，或者迭代次数超过200次，则此时所得即为所求光流值。

针对噪声影响下光流计算稳健性较差及收敛速度慢的问题，提出了一种噪声环境下光流场快速稳健估计方法，在ML算法基础上，引入惩罚因子以增强光流计算稳健性，在光流计算迭代公式中加入动量因子缩短光流计算收敛时间以加快光流场计算。

附图说明

图1为本发明实现的流程图。

图2为运动引起的面元成像灰度变化图。

图3为不同算法得到的光流图。

具体实施方式

以下结合附图和具体实施方式对本发明作进一步说明。

如图1所示为本发明实现的流程图。该方法包括以下步骤：

步骤1：光流约束方程

假设图像上一点(x,y)在t时刻的亮度为I(x,y,t)，在Δt时间后该像素点亮度变为I(x+Δx,y+Δy,t+Δt)，当Δt趋于无穷小时可认为该点亮度不变，得到等式：

将

改写为I_x,I_y,I_t，并令

则重写上式，得到光流计算的基本等式：

I_xu+I_yv+I_t＝0 (2)

u,v表示速度场矢量的两个分量。计算光流的过程就是求解u,v的过程，两个参量只有一个等式，因此为欠估计问题，为了确定光流速度，必须加入其他约束进行求解。

步骤2：构建能量函数方程

①光流约束方程修正

如图2所示，物体面元由位置P运动到位置P'，对应成像投影在像平面上由p移动到p'，假设物体表面为朗伯面，则p和p'的灰度分别为：

其中，ρ_d为物体表面漫反射系数，I_q为入射光强度，设为常数；N和N'分别为物体面元在位置P和P'的单位法向量；L为入射光的单位方向向量；θ为N和L的夹角；θ'为N'和L的夹角；由于面元在位置P和P'时的表面法线朝向不同，因此I_p≠I_p'。通常应用中，即使光照条件不变，连续相邻图像上物体同一元面对应像素灰度仍有可能是变化的。

MUKAWA根据物体运动、投影及光照模型证明灰度恒定假设在绝大多数情况下都不成立，并推导的基本方程的修正项为：

其中：

ΔN表示曲面法线运动变化；

从而将基本方程修正为：

cI+I_xu+I_yv+I_t＝0 (5)

基本方程约束就是要使光流场偏离基本等式的误差最小，定义光流约束因子为：

e≈(cI+I_xu+I_yv+I_t)² (6)

其中，c为参量，连同u,v采用全部平滑约束；全部平滑约束用梯度模的平方和来测度；

定义光流场的全局平滑约束因子为：

e₁＝||▽u||²+||▽v||²＝(u_x ²+v_x ²+u_y ²+v_y ²) (7)

定义参量c的全局平滑约束因子为:

e₂＝||▽c||²＝(c_x ²+c_y ²) (8)

实际应用中，视频图像不可避免会受到噪声污染，因此，准确估计光流场需要考虑噪声因素。设I'(x,y,z)为I(x,y,z)受加性噪声n(x,y,z)污染后图像点(x,y)的实际灰度值，则：

I(x,y,z)＝I(x,y,z)-n(x,y,z) (9)

高斯噪声在通信系统以及信号处理中最为常见。根据中心极限定理，在自然界中，一些现象受到许多相互独立随机因素的影响，如果每个因素所产生的影响都很微小，那么总的影响可以看作是服从正态分布。

基于此，假设式(9)中的噪声n(x,y,z)是均值为0，方差为δ²的高斯白噪声，则：

I(x,y,z)((x,y)∈Ω)未受噪声污染，故应满足全局平滑约束，定义灰度全局平滑约束因子为：

e₃＝||▽I||²＝(I_x ²+I_y ²) (12)

②引入惩罚项

为了增强变分模型稳健性，求得更精确的光流场，引入二次项β²(u²+v²+c²)作为惩罚项：

e₄＝β²(u²+v²+c²) (13)

求解出的光流场使得约束因子在上式中e,e₁,e₂,e₃,e₄约束噪声条件下达到最小。

综上，噪声环境下光流场估计问题可转化为式(11)约束下极小化：

其中，α₁,α₂,α₃为调和4类约束因子的权参数；

步骤3：求解光流场

①光流求解迭代公式

根据拉格朗日乘数法，式(11)约束下极小化式(14)即是极小化：

令：

F(x,y,I,u,v,c,I_x,I_y,u_x,u_y,v_x,v_y,c_x,c_y)＝e+α₁e₁+α₂e₂+α₃e₃+e₄+λ(I-I')-λδ²(16)

由变分法，式(15)取极小值的必要条件为：

将式(16)带入式(17)，化简后可得：

其中，s_xx+s_yy＝▽²s(s＝I,u,v,c)为s分量的拉普拉斯算子，可采用下式估算：

和s_i,j,k分别为t_k时刻(x_i,y_i)处s分量的取值及其邻域平均值，K＝1；

定义为：

由式(18)式中(d)可得：

(I_xu+I_yv+I_t+cI)＝α₂(c_xx+c_yy)/I (21)

化简得：

其中：

式(23)两端乘以(I-I')并在图像平面Ω上积分得：

将式(11)带入式(24)可得:

由式(22)和(25)进一步可得：

②添加动量因子后的光流求解迭代公式

由于光流场的求解是一个不断迭代的耗时过程，要得到稳定的方程解，一般都需要上百次迭代，因此为了缩短计算的收敛时间，同时得到十分稳定的光流场，故在光流迭代公式中添加动量因子以加快光流场计算，即：

显然要使(28)中uⁿ⁺¹，vⁿ⁺¹收敛，动量因子μ就必须满足|μ|<1的条件，这样才能保证

即只有在|μ|<1的条件下，随着迭代次数增加，才能保证算法最终收敛。其中：

至此，总结噪声环境下光流场快速稳健估计方法的迭代计算方法如下：

(1)读取视频中连续两帧图像，并将所需参数初始化；

(2)图像求x方向、y方向和t方向的导数，如公式(31)、(32)、(33)所示；

(3)对于I,u,v,c的均值，采用公式(20)进行计算。

(4)根据式(28)，并用Gauss-Seidel方法进行迭代，并计算两次迭代光流值的误差0.005，如果小于给定误差阈值，或者迭代次数超过200次，则此时所得即为所求光流值。

本发明的效果可通过以下仿真进一步说明：

仿真条件：

通过与HS光流算法、M法、ML法对比，本发明将在以下两个方面验证噪声环境下所提光流场快速稳健估计方法的有效性：稳健性、以及高效性。

仿真参数设置如下：α₁＝1，α₂＝0.5，α₃＝0.5，β＝0.02，初值u⁰、v⁰、c⁰均取0，I⁰取I'，迭代次数n＝200；当|uⁿ⁺¹-uⁿ|<β₁且|vⁿ⁺¹-vⁿ|<β₂时迭代终止，β₁，β₂为非负的阈值参数，β₁＝0.005，β₂＝0.005。

仿真内容：

图3(a)，(b)为视频中连续两帧图像，图3(c)，(d)，(e)依次是M法，ML法及本发明所提算法对加入均值为0，方差为5的高斯白噪声的视频图像计算所得光流场效果图。由于噪声原因，M法所得光流比较模糊，边缘光流扩散严重，ML法所得光流较M法而言，光流扩散明显减少，但所得光流仍较模糊，本发明所提算法较以上两种算法扩散明显减少，模糊度有所降低，因而所得光流场稳健性更好。

为了定量评估所提算法得到的光流场的稳健性，我们可采用如下评价指标：平均角误差

平均角度标准差σ_ψ和绝对平均误差

其中，角误差计算公式为：

上式中的

分别为实际光流场和估计光流场的时域扩展，假设v_c、v_e分别为时域间隔δt的两帧图像间的实际光流场和估计光流场，则：

令v_c＝(u_c,v_c)^T，v_e＝(u_e,v_e)^T，则：

平均角度误差：

平均角度标准差为：

绝对平均误差：

基于上述公式，可得HS、ML、M以及所提算法的各项指标如表1所示：

表1光流场评价指标

表1中，5,10,15代表噪声方差，平均角误差用E表示，平均角度标准差用F表示，绝对平均误差用H表示。由表1可以看出，对不同噪声水平图像进行光流估计，本发明方法的误差均为最小，因此能较大程度增强噪声环境下光流算法稳健性。

为了满足应用需求，算法需要具有较好的实时性，因此下面将对算法运行时间进行评估。考虑不同噪声环境，上述算法求得光流场所需时间如表2所示，其中5,10,15代表噪声方差。

表2光流场计算消耗时间

由表2可看出，HS、M、ML以及所提算法求解光流场所需时间随噪声方差的增加而增加，且与HS及ML相比，M算法由于在光源及目标表面进行了光几何约束，因而可以实现较快估计，进而运算时间较小，而本发明所提算法由于在迭代公式中加入了动量因子，可以加速收敛，因而运算时间叫M算法更小。

本发明从提高光流场计算稳健性及高效方面进行了详细的描述和研究。仿真结果表明，对视频中连续两帧图片加入不同高斯噪声后，与M算法及ML方法相比，所提算法可显著增强光流场计算稳健性，缩短光流计算收敛时间，加快光流场计算。

Claims (1)

1.噪声环境下光流场快速稳健估计方法，其特征在于，该方法包括以下步骤：

步骤1：光流约束方程

假设图像上一像素点(x,y)在t时刻的亮度为I(x,y,t)，在Δt时间后所述像素点亮度变为I(x+Δx,y+Δy,t+Δt)，当Δt趋于无穷小时认为该点亮度不变，得到等式：

将

改写为I_x,I_y,I_t，并令

则重写上式，得到光流计算的基本等式：

I_xu+I_yv+I_t＝0 (2)

u,v表示速度场矢量的两个分量；

步骤2：构建能量函数方程

①光流约束方程修正

设物体面元由位置P运动到位置P'，对应成像投影在像平面上由p移动到p'，假设物体表面为朗伯面，则p和p'的灰度分别为：

其中，ρ_d为物体表面漫反射系数，I_q为入射光强度，设为常数；N和N'分别为物体面元在位置P和P'的单位法向量；L为入射光的单位方向向量；θ为N和L的夹角；θ'为N'和L的夹角；

光流基本方程的修正项为：

其中：

ΔN表示曲面法线运动变化；

从而将基本方程修正为：

cI+I_xu+I_yv+I_t＝0 (5)

定义光流约束因子为：

e≈(cI+I_xu+I_yv+I_t)² (6)

其中，c为参量，连同u,v采用全部平滑约束；全部平滑约束用梯度模的平方和来测度；

定义光流场的全局平滑约束因子为：

定义参量c的全局平滑约束因子为:

设I'(x,y,z)为I(x,y,z)受加性噪声n(x,y,z)污染后图像点(x,y)的实际灰度值，则：

I(x,y,z)＝I(x,y,z)-n(x,y,z) (9)

假设式(9)中的噪声n(x,y,z)是均值为0，方差为δ²的高斯白噪声，则：

I(x,y,z)((x,y)∈Ω)未受噪声污染，故应满足全局平滑约束，定义灰度全局平滑约束因子为：

②引入惩罚项

引入二次项β²(u²+v²+c²)作为惩罚项：

e₄＝β²(u²+v²+c²) (13)

噪声环境下光流场估计问题转化为式(11)约束下极小化：

其中，α₁,α₂,α₃为调和4类约束因子的权参数；

步骤3：求解光流场

①光流求解迭代公式

根据拉格朗日乘数法，式(11)约束下极小化式(14)即是极小化：

令：

F(x,y,I,u,v,c,I_x,I_y,u_x,u_y,v_x,v_y,c_x,c_y)＝e+α₁e₁+α₂e₂+α₃e₃+e₄+λ(I-I')-λδ² (16)

由变分法，式(15)取极小值的必要条件为：

将式(16)带入式(17)，化简后得：

其中，

为s分量的拉普拉斯算子，采用下式估算：

和s_i,j,k分别为t_k时刻(x_i,y_i)处s分量的取值及其邻域平均值，K＝1；

定义为：

由式(18)式中(d)可得：

(I_xu+I_yv+I_t+cI)＝α₂(c_xx+c_yy)/I (21)

化简得：

其中：

式(23)两端乘以(I-I')并在图像平面Ω上积分得：

将式(11)带入式(24)可得:

由式(22)和(25)进一步可得：

②添加动量因子后的光流求解迭代公式

在光流迭代公式中添加动量因子以加快光流场计算，即：

其中：

其中，在对图像求x方向、y方向和t方向的导数时，采用的公式如下：

步骤4：迭代求解

根据式(28)，用Gauss-Seidel方法进行迭代，并计算两次迭代光流值的误差，如果小于0.005，或者迭代次数超过200次，则此时所得即为所求光流值。

Images