[发明专利]基于物理学的改善粒子图像测速稳健性光流方法

摘要

本发明针对可视化流动图像边缘扩散及噪声和异常点影响使得光流计算稳健性较差的问题，提出一种基于物理学的稳健光流计算方法以改善光流计算稳健性。所提算法在基于物理学光流方法中引入各向异性滤波器以增强边缘光流稳健性，并增加惩罚因子以减少噪声及异常点对光流计算的影响，而后基于变分方法极小化光流能量函数以求解欧拉‑拉格朗日方程，最后通过迭代方法求得速度场。仿真结果表明，与传统的Lucas‑Kanade、Horn‑Schunck、金字塔Lucas‑Kanade以及基于物理的光流计算方法相比，所提算法可显著减少边缘和角落区域光流扩散，改善针对噪声及异常点的稳健性，从而得到具有较好稳健性的速度场。

主权项

基于物理学的改善粒子图像测速稳健性光流方法，其特征在于：该方法包括如下步骤：步骤1：基于物理学的光流求解①基于物理学的光流约束方程∂g∂t+▿·(gu)=f--(1)]]>其中，u＝(u,v)代表图像平面上光流速度，g表示归一化图像强度，f中包含忽略一阶近似的扩散和边界条件，为拉普拉斯算子，D1是扩散系数，c代表粒子散射或吸收系数，表示边界项，表示零磁通条件，ψ代表纯量密度，通常情况下，当ψ为常量时，项可被忽略，表示投影在平面(x1,x2)上的梯度算子，Γ1,Γ2的取值与x1,x2相关，λ1表示拉格朗日乘子；②光流平滑约束方程极小化平滑约束项Es表示为：Es=[(∂u∂x)2+(∂u∂y)2+(∂v∂x)2+(∂v∂y)2]2dxdy---(2)]]>③基于物理学的光流求解结合基于物理学的光流方程(1)及光流平滑约束方程(2)，可得如下光流能量函数：E(u,v)=∫Ω[∂g∂t+▿.(gu)-f]2dx1dx2+λ∫Ω(||▿u||2+||▿v||2)dx1dx2---(3)]]>其中，λ是拉格朗日乘子，Ω代表图片中一个小区域；步骤2：基于物理学的稳健光流求解④各向异性滤波器为求得稳健性更高的光流，本发明采用各向异性扩散的光流计算方法，各向异性扩散的光流计算方法的扩散系数方程为：D(||▿I||)=e-α||▿I||β1---(4)]]>其中，常量K控制边缘的敏感度，通常由实验确定其数值；基于式(4)，将光流能量函数式(3)进一步改写为:E(u,v)=∫Ω[∂g∂t+▿.(gu)-f]2dx1dx2+λ∫Ω(||D▿u||2+||D▿v||2)dx1dx2---(5)]]>⑤惩罚项在式(5)中加入二次项作为惩罚项，加入惩罚项后，式(5)可改写为：E(u,v)=∫Ω[∂g∂t+▿.(gu)-f]2dx1dx2+β22(u2+v2)dx1dx2+λ∫Ω(||D▿u||2+||D▿v||2)dx1dx2---(6)]]>其中，β2是一个小的正常数；⑥基于物理学的稳健光流求解根据变分原理将式(6)最小化，可得如下欧拉‑拉格朗日(Euler‑Lagrange)方程：Lu-∂Lux∂x-∂Luy∂y=0Lv-∂Lvx∂x-∂Lvy∂y=0---(7)]]>其中：L=L(u,v,ux,uy,vx,vy)=[∂g∂t+▿·(gu)-f]2+β22(u2+v2)+λ(||D▿u||2+||D▿v||2)---(8)]]>Lu=2·[gt+▿·(gu)-f]*∂▿·(gu)∂u-2β22u∂Lux∂x+∂Luy∂y=2λD2▿2u---(9)]]>将式(8)及(9)带入式(7)可化简为：[gt+▿·(gu)-f]*∂▿·(gu)∂u-β22u-λD2▿2u=0[gt+▿·(gu)-f]*∂▿·(gu)∂v-β22v-λD2▿2v=0---(10)]]>其中，▽·(gu)＝g▽·(u)+u·▽g，可用某点与其周围速度平均值之差近似表示，即：▿2u=un+1-u‾n▿2v=vn+1-v‾n---(11)]]>结合式(11)，式(10)可重新表示为：(gt+ugx+vgy+gux+gvy-f)·gx-β2u-λD2(un+1-u‾n)=0(gt+ugx+vgy+gux+gvy-f)·gy-β2v-λD2(vn+1-v‾n)=0---(12)]]>对式(12)整理可得：(λD2+gx2-β22)un+1+vn+1.gxgy+gtgx+ggxux+ggxvy-f·gx-λD2u‾n=0(λD2+gy2-β22)vn+1+un+1.gxgy+gtgy+ggyux+ggyvy-f·gy-λD2v‾n=0---(13)]]>求解式(13)得un+1和vn+1：un+1=(λD2-β22)[gx(gt+g·u‾xn+g·v‾yn-f)-λD2u‾n]-λD2(u‾ngyn-v‾ngxgy)(λD2-β22)[(λD2-β22)2+gy2+gx2]vn+1=(λD2-β22)[gy(gt+g·u‾xn+g·v‾yn-f)-λD2v‾n]-λD2(v‾ngx2-u‾ngxgy)(λD2-β22)[(λD2-β22)2+gy2+gx2]---(14)]]>其中，u,v均值计算公式如下：使用sobel算子对滤波后的图像求x方向、y方向的导数；利用两帧中对应像素相减的方式求t方向的导数，所用的sobel模板如下所示：10-120-2101121000-1-2-1]]>步骤3：迭代求解根据式(14)，基于Gauss‑Seidel方法迭代求解，若相邻两次迭代光流之差小于给定阈值，或迭代次数超过给定值，则终止迭代。