CN103761749A - 一种基于非线性多重网格法的光流场快速估计方法 - Google Patents

Abstract

本发明涉及一种基于非线性多重网格法的光流场快速估计方法，包括：输入图像序列；将图像灰度化，并对图像进行降采样；运用变分法求与光流估计模型等价的离散化Euler-Lagrange方程组；采用非线性多重网格法求解光流；利用OpenMP并行编程模型对能量泛函的求解过程进行加速。本发明利用非线性多重网格法求解非线性光流模型，在每一个多网格循环中消去一个非线性残差，消除了线性多重网格法的求解误差，从而保证了光流的求解精度；本发明的数值求解方法只需迭代20～30次，与现有方法相比，大大提高了光流计算的实时性。同时，本发明利用OpenMP并行编程模型对能量泛函的求解过程进行加速，提高了程序的运行速度。

Description

一种基于非线性多重网格法的光流场快速估计方法

技术领域

本发明属于计算机视觉领域，具体涉及一种利用基于非线性多重网格法的光流场快速估计方法。

背景技术

光流是指空间运动物体在观测成像面上像素运动的瞬时速度，提供了有关场景中物体结构、位置与运动信息，是计算机视觉的重要组成部分。光流是基于像素点定义的，所有光流的集合称为光流场。光流场是运动图像分析技术的重要方法之一，在军事、工业和生活等应用领域具有重要的现实意义。

目前的光流计算方法采用光流模型的能量泛函均由数据项和平滑项组成。数据项是变分光流算法中能量泛函的主要组成部分，主要包含了各种常值守恒假设，例如灰度守恒假设、梯度守恒假设、Hessian矩阵守恒假设、拉普拉斯守恒假设等，这些守恒假设构成的约束条件是光流计算中决定运动模型的关键。平滑项是变分光流算法中能量泛函的又一重要组成部分，主要包含了各种平滑和分段平滑策略，它使变分光流模型取得唯一解。

Horn和Schunck在1981年引入基本光流约束方程及整体平滑约束条件，建立了光流计算的变分模型。该方法数据项采用灰度守恒假设，即：

I₀(x)＝I₁(x+h(x))

对其进行一阶泰勒展开，得：

I_xu(x)+I_yv(x)+I_t＝0

$I_x=\frac{\∂I}{\∂x},I_y=\frac{\∂I}{\∂y},I_t=\frac{\∂I}{\∂t}$

式中，I₀与I₁分别表示前后两帧图像的灰度，x=(x,y)^T表示图像上某个像素点，h(x)=(u(x),v(x))^T表示待求的光流矢量，u(x),v(x)分别代表该点在x，y方向上的光流分量。

光流约束方程中包含两个未知量u(x)和v(x)，因而无法求得方程的唯一解。Horn和Schunck引入全局平滑约束，从而组成光流模型的能量泛函，如下所示：

$E=\∫\∫[{(I_{x}u\left(x\right)+I_{y}v\left(x\right)+I_t)}^2+{\mathrm{\α}}^2({|\▿u\left(x\right)|}^2+{|\▿v\left(x\right)|}^2)]\mathrm{dxdy}$

上述为求解光流场的基本方程，其中α表示影响灰度守恒假设和全局平滑约束两者比例分配的权值。其平滑项为边缘保持的总变差项，该模型为线性模型。

由于光照变化、噪声、遮挡、大位移等因素的影响，采用该模型估计的光流场误差较大，不能满足应用中的要求。为克服这些因素对光流估计的影响，引入了多数据项守恒结合、张量滤波、梯度算子等方法提高了光流场的估计精度，同时使模型变为非线性，增加了算法的复杂度，从而降低了变分光流求解的实时性。

传统的数值求解方法有Gauss–Seidel迭代、Jaccobi迭代和松弛迭代SOR等方法，由于其需要迭代数千次才能求得较理想的结果，因此很难达到实时应用。多重网格法是近三十多年来发展起来的一类新的迭代法，对于求解由偏微分方程离散化得到的大规模线性方程组来说，它是目前最快速最高效的方法。

申请号为201310174158.6的专利提出了一种基于误差分布式多层网格的快速光流场计算方法，虽然提高了光流的求解速度，但是该方法针对的是线性模型的光流估计，并没有解决非线性模型的光流求解。为了提高光流的求解精度，目前的光流算法虽多采用非线性模型，但因求解方法复杂、速度慢，降低了算法的实时性。

发明内容

针对现有非线性光流求解模型中存在的实时性较差等问题，本发明提出一种基于非线性多重网格法的光流场快速估计方法，引入非线性多重网格法和OpenMP编程模型提高光流的求解速度。

下面给出本发明所述光流场计算方法的原理。

建立光流的非线性求解模型如下所示：

$E\left(h\right)={\∫}_{\mathrm{\Ω}}{[I_0\left(x\right)-I_1(x+h\left(x\right))]}^2\mathrm{dx}+{\mathrm{\λ}\∫}_{\mathrm{\Ω}}g\left(\right|\▿I\left|\right)\left({|\▿h\left(x\right)|}^2\right)\mathrm{dx}$

其中，λ为平滑项系数。该模型的数据项以灰度守恒假设为约束条件，由于像素位移过大时，用一阶泰勒展开对模型线性化使求解的光流误差显著增加，所以只能采用非线性的数据项；平滑项采用各向同性非线性扩散、图像驱动，通常情况下运动边缘是图像边缘的子集，采用该平滑项可以保护运动边缘的光流，使其不被过于平滑和模糊，从而可以保护光流边缘；边缘函数g(r)为传导系数，只与图像数据有关，与光流数据无关，其表达式如下：

$g\left(r\right)=\frac{1}{{1+(r/k)}^2}$

其中，k是选定的常数，用于控制边缘函数g(r)的下降速率，本发明取k=1。

针对选取的光流求解模型，使用最速下降法获得非线性偏微分方程组，从而将求解光流的能量泛函最小化问题转换成求解一个非线性偏微分方程组的问题。然后采用FAS多重网格法求解该微分方程组，并获得光流。

一种基于非线性多重网格法的光流场快速估计方法，主要包含以下步骤：

步骤一，输入图像序列。

步骤二，将图像灰度化，并对图像进行降采样。

步骤三，运用变分法求与光流估计模型等价的离散化Euler-Lagrange方程组。

运用变分法，对非线性光流能量泛函求偏导，得到与光流估计模型等价的Euler-Lagrange方程组：

$\left\{\begin{array}{c}{I}_x\·(I_1-I_0)+\mathrm{\λg}\left(\right|\▿I\left|\right)|\▿u|-\frac{\mathrm{\λ}}{2}\mathrm{div}(g\left(\right|\▿I\left|\right)\▿u)=0\\ I_y\·(I_1-I_0)+\mathrm{\λg}\left(\right|\▿I\left|\right)|\▿v|-\frac{\mathrm{\λ}}{2}\mathrm{div}\left(g\right|\▿I|\▿v)=0\end{array}\right.$

式中，I为图像的灰度值，I₀与I₁分别表示前后两帧图像的灰度，u,v分别代表该点在x，y方向上的光流分量，λ为平滑项系数，g(·)为边缘函数。

对方程组进行离散化，设定需要满足的精度要求并用数值方法进行逼近求解。Euler-Lagrange方程组的离散形式如下：

$\left\{\begin{array}{c}{I}_{\mathrm{xi},j}{u}_{i,j}\·(I_{1i,j}-I_{0i,j})+{\mathrm{\λg}}_{i,j}{m}_1-\frac{\mathrm{\λ}}{2}\·(g_{i+1,j}{u}_{i+1,j}+g_{i-1,j}{u}_{i-1,j}+g_{i,j+1}{u}_{i,j+1}+g_{i,j-1}{u}_{i,j-1}-{4g}_{i,j}{u}_{i,j})=0\\ I_{\mathrm{yi},j}{v}_{i,j}\·(I_{1i,j}-I_{0i,j})+{\mathrm{\λg}}_{i,j}{m}_2-\frac{\mathrm{\λ}}{2}\·(g_{i+1,j}{v}_{i+1,j}+g_{i-1,j}{v}_{i-1,j}+g_{i,j+1}{v}_{i,j+1}+g_{i,j-1}{v}_{i,j-1}-{4g}_{i,j}{v}_{i,j})=0\end{array}\right.$

式中，i，j表示像素点在图像中的位置。

步骤四，运用非线性多重网格法求解上述的Euler-Lagrange方程组。

多重网格是加速光流迭代收敛的有效方法，现有的光流加速计算中，大部分都是利用线性多重网格方法针对线性光流模型加速计算。由于线性模型求解的光流精度不够，而对线性模型改进后的非线性模型具有更好的鲁棒性，因此不能运用线性多重网格方法对非线性光流计算模型的求解进行加速。解决方法是采用非线性多重网格方法来计算非线性变分光流模型，从而在保持非线性模型所求解光流精度的情况下，提高算法迭代收敛的效率。采用FAS多重网格法的具体方法如下：

（1）定义非线性多重网格法求解过程中的参数。

非线性多重网格的方程为：

A_h(x_h)＝f_h

式中，A_h是一个非线性算子，x_h是光流矢量，f_h是与光流矢量无关的非线性矢量，上式对应的离散形式为：

$\left\{\begin{array}{c}[I_{\mathrm{xi},j}\·(I_{1i,j}-I_{0i,j})+{2\mathrm{\λg}}_{i,j}]u_{i,j}+{\mathrm{\λg}}_{i,j}{m}_1=\frac{\mathrm{\λ}}{2}\·(g_{i+1,j}{u}_{i+1,j}+g_{i-1,j}{u}_{i-1,j}+g_{i,j+1}{u}_{i,j+1}+g_{i,j-1}{u}_{i,j-1})\\ [I_{\mathrm{yi},j}\·(I_{1i,j}-I_{0i,j})+{2\mathrm{\λg}}_{i,j}]v_{i,j}+{\mathrm{\λg}}_{i,j}{m}_2=\frac{\mathrm{\λ}}{2}\·(g_{i+1,j}{v}_{i+1,j}+g_{i-1,j}{v}_{i-1,j}+g_{i,j+1}{v}_{i,j+1}+g_{i,j-1}{v}_{i,j-1})\end{array}\right.$

设

是在细网格上经过非线性光滑子处理得到的结果，非线性残量方程定义为：

$A_h{x}_h-A_h{\stackrel{\‾}{x}}_h=A_h({\stackrel{\‾}{x}}_h+e_h)-A_h{\stackrel{\‾}{x}}_h=f_h-A_h{\stackrel{\‾}{x}}_h=r_h$

式中，e_h表示误差，r_h表示残差。

为了在细网格上校正光流，需要计算误差。但是在细网格上不能直接计算误差e_h，需将残差方程转移到粗网格上，粗网格校正就是将粗网格上的残差插值回细网格得到e_h，从而得到光流的一个更新的近似值。粗网格上的残差e_H为：

$e_H=x_H-{\stackrel{\‾}{x}}_H$

式中，x_H为粗网格上的光流矢量，为x_H在粗网格上经非线性光滑子处理后的结果。

（2）进行前光滑，消除高频分量。

以x_h为初值，对细网格方程迭代m次得到

${\stackrel{\‾}{x}}_h={\mathrm{smoother}}^m(A_h,x_h,f_h)$

大多数光滑迭代都有非线性形式，如Gauss-Seidel迭代、Jaccobi迭代、松弛迭代SOR、半隐式不动点迭代等迭代方法。这里选择Gauss-Seidel迭代。

（3）进行粗网格校正，将余量限制在粗网格上，消除低频分量。

残差为：

$r_h={A_{h}x}_h-A_h{\stackrel{\‾}{x}}_h$

将残差限制到粗网格上：

$r_H=A_H\left(I_h^H{x}_h\right)-I_h^H{A}_H\left(x_h\right)$

式中，

为限制算子。

将细网格上的余量r_h限制为粗网格上的f_H，网格限制示意图如图2所示。限制算子可以采用直接映射或者加权平均方法，这里选择直接映射法，如图3所示。图2和图3中圆形网格点和菱形网格点均代表图像上的像素点，经过限制算子处理后的网格点用圆形网格点表示。直接映射法的数学表达示为：

$v_{i,j}^H=v_{2i,2j}^h$

式中，i和j表示网格点的位置。

从而在粗网格上求得近似解：

${\stackrel{\‾}{x}}_H=I_h^H{\stackrel{\‾}{x}}_h$

粗网格上线性方程组右边的表达式为：

$f_H=r_H+A_H{\stackrel{\‾}{x}}_H$

（4）将步骤（3）求得的误差量延拓到细网格上。

如果当前网格尺度为最细网格，采用直接迭代法或快速迭代法求解方程：A_He_H＝f_H；否则，以为初始解，执行FAS多重网格迭代法求解e_H。校正误差为：

${\stackrel{\‾}{e}}_H=e_H-{\stackrel{\‾}{x}}_H$

（5）进行细网格校正。

对粗网格上传回的误差进行插值修正：

${\stackrel{\‾}{e}}_h=I_H^h{\stackrel{\‾}{e}}_H$

式中，

为延拓算子，选用双线性插值算子。

模板的中心对准一个细网格点，它的值就是该点周围粗网格点的加权平均值。所有的权值可以从模板中对应的位置得到，定义为：

$v_{2i,2j}^h=({9v}_{i,j}^H{+3v}_{i+1,j}^H+{3v}_{i,j+1}^H+v_{i+1,j+1}^H)/16$

$v_{2i+1,2j}^h=({9v}_{i+1,j}^H{+3v}_{i,j}^H+{3v}_{i+1,j+1}^H+v_{i,j+1}^H)/16$

$v_{2i,2j+1}^h=({9v}_{i,j+1}^H+{3v}_{i,j}^H+{3v}_{i+1,j+1}^H+v_{i+1,j}^H)/16$

$v_{2i+\mathrm{1,2}j+1}^h=({9v}_{i+1,j+1}^H+{3v}_{i+1,j}^H+{3v}_{i,j+1}^H+v_{i,j}^H)/16$

用模板可以表示为：

计算新的近似解，得到细网格校正后的解：

${\hat{x}}_h={\stackrel{\‾}{x}}_h+{\stackrel{\‾}{e}}_h$

（6）进行后光滑。

细网格校正之后，由于加入了粗网格传回的误差，以

为初值，对细网格方程迭代n次得到x_h，消除误差中的高频分量：

$x_h={\mathrm{smoother}}^n(A_h,{\hat{x}}_h,f_h)$

步骤五，利用OpenMP并行编程模型对能量泛函的求解过程进行加速。

OpenMP能够为编写多线程应用程序提供一种简单的方法，针对单主机上多核/多CPU并行计算而设计的工具，效率很高、内存开销小、编程语句简洁直观，因此编程容易、编译器实现也容易（现在最新版的C、C++、Fortran编译器基本上都内置OpenMP支持），无须程序员进行复杂的线程创建、同步、负载平衡和销毁工作。OpenMP是一个外部的编程模型，而不是自动编程模型，它能够使程序员完全控制并行化。当程序开始执行时，只存在一个主线程，程序执行为串行模式，直至遇到并行域才开始执行并行运算。

在每个并行域中，分配多个任务到不同的线程中，采用标准的OpenMP任务分配编译指导语句，借助OpenMP库函数omp_get_num_threads()和pragma omp parallel进行任务划分，借助编译制导语句parallel for循环访问一个索引范围，并在每次迭代时以并行方式执行循环内的函数，从而使编译器可以自动将程序进行并行化。

与现有技术相比，本发明具有以下优点：

1.利用非线性多重网格法求解非线性光流模型，在每一个多网格循环中消去一个非线性残差，消除了线性多重网格法的求解误差，从而保证了光流的求解精度；

2.与现有的原始对偶、Gauss-Seidel迭代等数值求解方法相比，现有方法需要进行百次或千次迭代，而本发明只需二十至三十次迭代，提高了光流计算的实时性；

3.利用OpenMP并行编程模型，在现有的算法基础上，提高了程序的运行速度。

附图说明

图1为基于非线性多重网格法的光流场快速估计方法流程图；

图2为网格限制示意图；

图3为直接映射法示意图；

图4为Middlebury标准库中两组图像序列及其彩色光流图。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步阐释。

本发明的实施在英特尔Core i3-2100双核CPU，主频3.10GHz，内存2G的PC机上进行实验，操作系统为Windows XP，软件开发环境为Visual Studio2008软件平台。

本实施方式的具体过程为，首先对输入的彩色图像灰度化并采用高斯滤波去噪，预处理后的图像运用到后面的光流估计中；其次确定非线性光流估计模型，采用变分法得到模型的离散形式，并建立3～4层不同粗细程度的网格图像，使用非线性多重网格法求解光流；然后利用OpenMP并行编程模型加速该算法以提高算法实时性；最后利用光流误差评价函数计算该算法的误差。图1是本发明的方法流程，具体包括以下几个步骤：

步骤一，输入连续两帧图像数据。

步骤二，图像灰度化，并将图像降采样。

步骤三，运用变分法，求出与光流估计模型等价的离散化Euler-Lagrange方程组。

步骤四，采用非线性多重网格法求解光流。

步骤五，利用OpenMP并行编程模型对能量泛函的求解过程进行加速。

下面给出本发明的一个应用实例。

为验证本发明提出的基于非线性多重网格法的光流场快速估计方法，从标准光流数据库Middlebury中选取了2组序列图，如图4所示，a1，b1和a2，b2分别为Grove3和Hydrangea两组图像序列的连续两帧，所画出的彩色光流图如c1、c2所示。

为了与现有技术进行比较，分别采用传统光流求解方法（方法一）、在本发明提出的方法基础上去除OpenMP并行编程模型部分（方法二）和本发明提出的基于非线性多重网格法的光流场快速估计方法（方法三）进行实验，实验仍然采用图4中的a1、b1和a2、b2两组图像序列。

采用平均角误差(average angular error,AAE)和平均端点误差(average endpoint error,EPE)两种光流精度评价指标，从而比较两种不同方法所估计的光流精度。

AAE反映了计算的光流矢量场整体偏离标准光流矢量场的程度，其计算公式如下：

$\mathrm{AAE}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\φ}}_{\mathrm{aae}}\left(i\right)$

其中：

${\mathrm{\φ}}_{\mathrm{aae}}\left(i\right)=\arccos \left[\frac{u_i^s{u}_i^n+v_i^s{v}_i^n+k^2}{\sqrt{{\left(u_i^s\right)}^2+{\left(v_i^s\right)}^2+k^2\sqrt{{\left(u_i^n\right)}^2+{\left(v_i^n\right)}^2+k^2}}}\right]$

N表示一帧图像的像素个数，

表示第i个像素的标准光流矢量，

表示计算得到的第i个像素的光流矢量，k表示相隔的帧数。

EPE用以衡量计算的光流场的矢量长度与标准光流场的矢量长度之间的误差，其计算公式如下：

$\mathrm{EPE}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\φ}}_{\mathrm{epe}}\left(i\right)$

其中：

${\mathrm{\φ}}_{\mathrm{epe}}\left(i\right)=\sqrt{{(u_i^s-u_i^n)}^2+{(u_i^s-v_i^n)}^2}$

表1给出了两种方法的AAE、EPE及计算所用时间。由表1可知，本发明所述方法与传统光流求解方法的AAE、AEPE差别很小，但本发明所述方法的计算时间却明显小于传统光流求解方法。因此，本发明在保证光流精度的同时，提高了光流的求解速度，具有良好的实时性。

表1本发明与现有技术计算误差和速度的对比

Claims (3)

1.一种基于非线性多重网格法的光流场快速估计方法，其特征在于，引入非线性多重网格法和OpenMP编程模型提高光流的求解速度；所述方法包括以下步骤：

步骤一，输入图像序列；

步骤二，将图像灰度化，并对图像进行降采样；

步骤三，运用变分法求与光流估计模型等价的离散化Euler-Lagrange方程组；

运用变分法，对非线性光流能量泛函求偏导，得到与光流估计模型等价的Euler-Lagrange方程组：

$\left\{\begin{array}{c}{I}_x\·(I_1-I_0)+\mathrm{\λg}\left(\right|\▿I\left|\right)|\▿u|-\frac{\mathrm{\λ}}{2}\mathrm{div}(g\left(\right|\▿I\left|\right)\▿u)=0\\ I_y\·(I_1-I_0)+\mathrm{\λg}\left(\right|\▿I\left|\right)|\▿v|-\frac{\mathrm{\λ}}{2}\mathrm{div}\left(g\right|\▿I|\▿v)=0\end{array}\right.$

式中，I为图像的灰度值，I₀与I₁分别表示前后两帧图像的灰度，u,v分别代表该点在x，y方向上的光流分量，λ为平滑项系数，g(·)为边缘函数；

对方程组进行离散化，设定需要满足的精度要求并用数值方法进行逼近求解；Euler-Lagrange方程组的离散形式如下：

$\left\{\begin{array}{c}{I}_{\mathrm{xi},j}{u}_{i,j}\·(I_{1i,j}-I_{0i,j})+{\mathrm{\λg}}_{i,j}{m}_1-\frac{\mathrm{\λ}}{2}\·(g_{i+1,j}{u}_{i+1,j}+g_{i-1,j}{u}_{i-1,j}+g_{i,j+1}{u}_{i,j+1}+g_{i,j-1}{u}_{i,j-1}-{4g}_{i,j}{u}_{i,j})=0\\ I_{\mathrm{yi},j}{v}_{i,j}\·(I_{1i,j}-I_{0i,j})+{\mathrm{\λg}}_{i,j}{m}_2-\frac{\mathrm{\λ}}{2}\·(g_{i+1,j}{v}_{i+1,j}+g_{i-1,j}{v}_{i-1,j}+g_{i,j+1}{v}_{i,j+1}+g_{i,j-1}{v}_{i,j-1}-{4g}_{i,j}{v}_{i,j})=0\end{array}\right.$

式中，i，j表示像素点在图像中的位置；

步骤四，运用非线性多重网格法求解所述步骤3得到的Euler-Lagrange方程组；

采用非线性多重网格方法计算非线性变分光流模型，在保持非线性模型所求解光流精度的情况下，提高算法迭代收敛的效率；

步骤五，利用OpenMP并行编程模型对能量泛函的求解过程进行加速。

2.根据权利要求1所述的一种基于非线性多重网格法的光流场快速估计方法，其特征在于，所述步骤四采用FAS多重网格法求解光流，具体方法如下：

（1）定义非线性多重网格法求解过程中的参数；

非线性多重网格的方程为：

A_h(x_h)＝f_h

式中，A_h是一个非线性算子，x_h是光流矢量，f_h是与光流矢量无关的非线性矢量，上式对应的离散形式为：

$\left\{\begin{array}{c}[I_{\mathrm{xi},j}\·(I_{1i,j}-I_{0i,j})+{2\mathrm{\λg}}_{i,j}]u_{i,j}+{\mathrm{\λg}}_{i,j}{m}_1=\frac{\mathrm{\λ}}{2}\·(g_{i+1,j}{u}_{i+1,j}+g_{i-1,j}{u}_{i-1,j}+g_{i,j+1}{u}_{i,j+1}+g_{i,j-1}{u}_{i,j-1})\\ [I_{\mathrm{yi},j}\·(I_{1i,j}-I_{0i,j})+{2\mathrm{\λg}}_{i,j}]v_{i,j}+{\mathrm{\λg}}_{i,j}{m}_2=\frac{\mathrm{\λ}}{2}\·(g_{i+1,j}{v}_{i+1,j}+g_{i-1,j}{v}_{i-1,j}+g_{i,j+1}{v}_{i,j+1}+g_{i,j-1}{v}_{i,j-1})\end{array}\right.$

设是在细网格上经过非线性光滑子处理得到的结果，非线性残量方程定义为：

$A_h{x}_h-A_h{\stackrel{\‾}{x}}_h=A_h({\stackrel{\‾}{x}}_h+e_h)-A_h{\stackrel{\‾}{x}}_h=f_h-A_h{\stackrel{\‾}{x}}_h=r_h$

式中，e_h表示误差，r_h表示残差；

为了在细网格上校正光流，需要计算误差；但是在细网格上不能直接计算误差e_h，需将残差方程转移到粗网格上，粗网格校正就是将粗网格上的残差插值回细网格得到e_h，从而得到光流的一个更新的近似值；粗网格上的残差e_H为：

$e_H=x_H-{\stackrel{\‾}{x}}_H$

式中，x_H为粗网格上的光流矢量，

为x_H在粗网格上经非线性光滑子处理后的结果；

（2）进行前光滑，消除高频分量；

以x_h为初值，对细网格方程进行m次Gauss-Seidel迭代得到xh：

${\stackrel{\‾}{x}}_h={\mathrm{smoother}}^m(A_h,x_h,f_h)$

（3）进行粗网格校正，将余量限制在粗网格上，消除低频分量；

残差为：

$r_h=f_h-A_h{\stackrel{\‾}{x}}_h$

将残差限制到粗网格上：

$r_H=A_H\left(I_h^H{x}_h\right)-I_h^H{A}_H\left(x_h\right)$

式中，

为限制算子；

将细网格上的余量r_h限制为粗网格上的f_H，限制算子采用直接映射或者加权平均方法；直接映射法的数学表达示为：

$v_{i,j}^H=v_{2i,2j}^h$

式中，i和j表示网格点的位置；

从而在粗网格上求得近似解：

${\stackrel{\‾}{x}}_H=I_h^H{\stackrel{\‾}{x}}_h$

粗网格上线性方程组右边的表达式为：

$f_H=r_H+A_H{\stackrel{\‾}{x}}_H$

（4）将步骤（3）求得的误差量延拓到细网格上；

如果当前网格尺度为最细网格，采用直接迭代法或快速迭代法求解方程：A_He_H＝f_H；否则，以

为初始解，执行FAS多重网格迭代法求解e_H；校正误差为：

${\stackrel{\‾}{e}}_H=e_H-{\stackrel{\‾}{x}}_H$

（5）进行细网格校正；

对粗网格上传回的误差进行插值修正：

${\stackrel{\‾}{e}}_h=I_H^h{\stackrel{\‾}{e}}_H$

式中，为延拓算子，选用双线性插值算子；

模板的中心对准一个细网格点，它的值就是该点周围粗网格点的加权平均值；所有的权值可以从模板中对应的位置得到，定义为：

$v_{2i,2j}^h=({9v}_{i,j}^H{+3v}_{i+1,j}^H+{3v}_{i,j+1}^H+v_{i+1,j+1}^H)/16$

$v_{2i+1,2j}^h=({9v}_{i+1,j}^H{+3v}_{i,j}^H+{3v}_{i+1,j+1}^H+v_{i,j+1}^H)/16$

$v_{2i,2j+1}^h=({9v}_{i,j+1}^H+{3v}_{i,j}^H+{3v}_{i+1,j+1}^H+v_{i+1,j}^H)/16$

$v_{2i+\mathrm{1,2}j+1}^h=({9v}_{i+1,j+1}^H+{3v}_{i+1,j}^H+{3v}_{i,j+1}^H+v_{i,j}^H)/16$

用模板可以表示为：

计算新的近似解，得到细网格校正后的解：

${\hat{x}}_h={\stackrel{\‾}{x}}_h+{\stackrel{\‾}{e}}_h$

（6）进行后光滑；

细网格校正之后，由于加入了粗网格传回的误差，以

为初值，对细网格方程迭代n次得到x_h，消除误差中的高频分量：

$x_h={\mathrm{smoother}}^n(A_h,{\hat{x}}_h,f_h).$

3.根据权利要求1所述的一种基于非线性多重网格法的光流场快速估计方法，其特征在于，所述步骤五利用OpenMP并行编程模型对能量泛函的求解过程进行加速的方法如下：

在每个并行域中，分配多个任务到不同的线程中，采用标准的OpenMP任务分配编译指导语句，借助OpenMP库函数omp_get_num_threads()和pragma omp parallel进行任务划分，借助编译制导语句parallel for循环访问一个索引范围，并在每次迭代时以并行方式执行循环内的函数，从而使编译器可以自动将程序进行并行化。

Images