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

权利要求书

1.基于物理学的改善粒子图像测速稳健性光流方法，其特征在于：该方法包括如下步骤：

步骤1：基于物理学的光流求解

①基于物理学的光流约束方程

∂g∂t+▿·(gu)=f--(1)]]>

其中，u＝(u,v)代表图像平面上光流速度，g表示归一化图像强度，f中包含忽略一阶近似的扩散和边界条件，为拉普拉斯算子，D₁是扩散系数，c代表粒子散射或吸收系数，表示边界项，表示零磁通条件，ψ代表纯量密度，通常情况下，当ψ为常量时，项可被忽略，表示投影在平面(x₁,x₂)上的梯度算子，Γ₁,Γ₂的取值与x₁,x₂相关，λ₁表示拉格朗日乘子；

②光流平滑约束方程

极小化平滑约束项E_s表示为：

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)]]>

其中，β₂是一个小的正常数；

⑥基于物理学的稳健光流求解

根据变分原理将式(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)得uⁿ⁺¹和vⁿ⁺¹：

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均值计算公式如下：

u‾i,j,k=16(u‾i,j-1,k+u‾i,-1j,k+u‾i,j+1,k+u‾i+1,j,k)+112(u‾i-1,j-1,k+u‾i-1,j+1,k+u‾i+1,j+1,k+u‾i,j-1,k)v‾i,j,k=16(v‾i,j-1,k+v‾i,-1j,k+v‾i,j+1,k+v‾i+1,j,k)+112(v‾i-1,j-1,k+v‾i-1,j+1,k+v‾i+1,j+1,k+v‾i,j-1,k)]]>

使用sobel算子对滤波后的图像求x方向、y方向的导数；利用两帧中对应像素相减的方式求t方向的导数，所用的sobel模板如下所示：

10-120-2101121000-1-2-1]]>

步骤3：迭代求解

根据式(14)，基于Gauss-Seidel方法迭代求解，若相邻两次迭代光流之差小于给定阈值，或迭代次数超过给定值，则终止迭代。