专利名称：用于粒子图像测速中的多网格处理方法
技术领域：
本发明涉及图像处理，特别是涉及一种用于粒子图像测速中的多网格 处理方法。
背景技术：
自然界图像的运动形式千差万别，流体运动是一种典型的非刚体运
动，流体运动图像计算与分析即粒子图像测速(Particle Image Velocimitry，简称PIV)是一种新型的非接触式的测量技术。通常PIV是 在流体中投入示踪粒子，在激光片光的照射下，在与片光垂直方向用摄像 机拍摄随流体运动的粒子图像，再对图像进行处理分析计算，最终得到流 场横切面上二维速度场的一种测量方法。作为一种全流场、无接触、无扰 动、高精度的流动可视化方法，PIV适用于湍流、非定常流等时变复杂流 场的测量。如今，PIV已是一门跨学科的交叉综合技术，其结合了激光技 术、视频图像处理技术、计算机技术、流体力学和近代光学技术的最新成 果，广泛应用于实验流体力学、生物医学、工业制造等许多领域。PIV将整 个运动矢量场作为一个整体连续函数进行估计，其分辨率可以达到图像中 的每一个像素，为微观、精确地进行局部图像运动计算和分析提供了可能， 同时克服了传统相关方法中的一些固有缺陷。PIV系统在实验流体力学、 空气动力学、制造业、医学、工业、飞机制造、水利水电、气象学等领域 都有着重要的科学和经济价值。
许多图像分析问题需要解决带有边值问题的椭圆偏微分方程，并最终 通过求解线性或非线性方程组实现如基于变分方法的图像非线性扩散、 主动轮廓计算、图像恢复等。虽然光流计算无论在计算精度、可靠性等方 面都取得了较大的飞跃，但其数值方法和手段并没有重大突破。
多网格思想提出于上世纪60年代，70年代后被应用于实际科学计算 中。如今多网格计算已经广泛应用于物理学、力学、电磁场等科学计算，工 程数学和有限元数值计算中广泛使用的多网格方法。网格是偏微分方程数 值解法的基础，网格体系的好坏直接影响计算结果的精度。多网格方法的 研究经历了从结构化到非结构化，从单一网格到混合网格的过程，并不断
3出现新的针对不同情况的网格生成技术。
发明内容
本发明的主要目的就是针对现有技术的不足，提供一种多网格处理方 法，能有效加快粒子图像测速的图像处理过程。 为实现上述目的，本发明采用以下技术方案
一种用于粒子图像测速中的多网格处理方法,所述粒子图像测速包括
求解线性方程组AU-F以实现对运动场的最佳逼近的步骤，U为待求解， A为系数矩阵，其至少包含运动矢量场参数，F为向量，所述方法包括以
下步骤
A、 设置第一网格，获得方程组AU-F对应在第一网格h下的方程组 AhUh = Fh;
B、 迭代求解方程组AhUh-Fh，根据所得实际解tJh得到方程误差 Rh =Fh-AhtJh，
C、 设置尺寸较第一网格更粗的第二网格H，将在第一网格h下对方程 组AV =Rh的求解尺度映射为在第二网格H下对方程组AV1 =RH的求解；
D、 迭代求解在第二网格H下的高频计算误差eH;
E、 对第二网格H下的高频计算误差eH进行尺度映射，得到在第一网 格h下的低频计算误差e、
F、 用第一网格h下的低频计算误差eh更新实际解tI|;e,v=eh+Ch。 优选地，还包括
在更新实际解后，再在第一网格h下迭代求解AhUh = Fh ，来消除t^ew的 高频计算误差。
优选地，循环多次进行所述处理，其中第二网格的尺度逐次递增。 本发明有益的技术效果是
图像分析问题需要解决带有边值问题的椭圆偏微分方程，并最终通过 求解线性或非线性方程组实现，如基于变分方法的图像非线性扩散、主 动轮廓计算、图像恢复等，在大尺寸网格上线性方程的高频误差很容易消 除，相反即使经过大量迭代低频误差却很难消除，但在小尺度网格上这种 低频误差只需经少量迭代就能明显消除。因此，在粒子图像测速的图像处 理过程中，采用多网格方法，通过不同网格的迭代进行相关线性方程组的 求解处理，可以有效加快线性方程组迭代的收敛过程，大大加快图像处理 速度，同时不会降低处理的可靠性和精度。采用网格多尺度(有可能尺度因子小于2)逐级逼近的方法，可以降低 方程模型误差带来的影响，改善计算精度。每次尺度变化都会进行一次基 于当前估计结果的补偿更新。将图像尺度的变化看作是不同尺寸网格的变 化，本质上多尺度逼近计算可以看作是多网格算法的一种特例。
图1为粒子图像测速处理系统的结构框图2为用水平集函数来划分内部光滑连续而边界不连续的场的示意
图3为一种实施例中利用多网格进行加速处理的流程图； 图4a和图4b分别为V型和W型多网格结构的示意图； 图5为一种优选实施例的多网格结构的示意图。 本发明的特征及优点将通过实施例结合附图进行详细说明。
具体实施例方式
首先介绍可应用本发明实施例的多网格处理方法的一种粒子图像测 速方法。
一个完整的粒子图像测速处理系统关键部分如图1所示，主要包括数 据采集模块、图像预处理模块、运动矢量场(光流场)计算模块和后续数据 分析模块，后续数据分析模块包括旋度计算模块和散度计算模块。其中的 运动矢量场计算模块是整个PIV系统的核心模块。
在计算2维运动矢量场、散度和旋度场时，采用了基于光流运动方程 的计算思想，以整体函数的形式计算运动矢量。为了能够抑止噪声同时又 能保持矢量场、散度和旋度场的局部结构信息，还结合了水平集的思想。
2维运动矢量场是指在大小为NXM像素范围内(如通常熟知的PAL制 768 X 576或640 X 480个像素)能够描述每一像素位置(NXM内)下该像素当
前时刻的位移矢量，该矢量为一浮点数二维矢量，表示为v=|;Vx,Vy;|，
其中、，vy为x轴和y轴方向的两个速度分量。为了统一描述，用大写粗体
表示矢量，小写细体表示标量。速度场V实际是坐标位置变量X气X,y]的函数。
流体运动图像计算出了获得瞬时运动矢量场v外，还需计算v的散度场《和 旋度场n，如(1)(2)式所示。
<formula>formula see original document page 5</formula><formula>formula see original document page 6</formula>
由于实际的流体运动可能会有多矢量流交界的现象，即运动边界和多 运动问题，为了解决运动边界和多运动问题，同时为了提高流体运动特征 描述的精度，将水平集的策略和思想集成于原有的光流计算框架。
如图2所示，在某一区域Q内运动矢量、散度和旋度有可能不连续， 为此可以将Q划分为两个子区域^,Q,， Q,区域内对应的运动场、散度场和 旋度场为[Hn,]; Q,区域内对应的运动场、散度场和旋度场为[V2,m2]。 假设在区域q内[Hn,光滑连续，在区域Q,内[V^,T^光滑连续，但在Q^2 的交界处运动场、散度场和旋度场并不连续。为了保证计算时[V^,t!,和 [V",，ri,]的连续性，同时又不破坏^与Q,之间的结构信息以及便于计算，构 建水平集函数(()和指示函数H ， c))为定义在图像大小范围内(NXM)的连续函 数。H(x)-l,x20and H(x)=0,x<0。
用水平集函数小来划分^,Q:， po的区域属于Q,，而"0的区域属于^。 由于采用了水平集的思想，因此在计算[Hr^和[V",、]时同时要不断计
算更新(t)o
为了估计出[Hn,]和[v"^]，构造以下能量函数。
E二E,+XE2 (3) 其中E为全局能量泛函，EpE2分别为以下式(4)、 (5)形式的数据约束
能量和光滑约束能量。
<formula>formula see original document page 6</formula>
这里I(X,t)代表粒子图像序列。最终估计的光流场可以描述为<formula>formula see original document page 7</formula>通过全局能量泛函的最小化来实现光流场的最佳逼近。
能量泛函(3)的最小化可以通过求解线性方程组(6)实现，。
AU = F (6) 其中U-[VpV2,^，^，TipTi2,小]，A为9X9的系数矩阵，F为9X1的向
量表示线性方程的奇次项。
线性方程组(6)的获得和计算可以通过经典的变分计算方法实现，这 里不详细阐述。
对于上述的PIV方法，可以采用多网格加速策略与之相结合，具体而 言，即通过多网格方法来处理线性方程组(6)，加速PIV计算。
多网格计算的一个中心任务是求解各种各样来自实际问题的偏微分 方程及线性方程组。由于许多图像分析问题需要解决带有边值问题的椭圆 偏微分方程，并最终通过求解线性或非线性方程组实现，如基于变分方 法的图像非线性扩散、主动轮廓计算、图像恢复等，在大尺寸网格上线性 方程的高频误差很容易消除，相反即使经过大量迭代低频误差却很难消除， 但在小尺度网格上这种低频误差只需经少量迭代就能明显消除。因此多网 格策略就是通过不同网格的迭代计算来加速收敛过程。对于图像分析而言， 网格的定义比较简单，原始图像的每一个像素表示一个网格，即最细的网 格，计算都是基于这样的离散网格点上的值运算完成的。
请参考图3， 一种实施例的多网格处理方法具体包括如下步骤 在网格尺寸为h,xhy的第一网格h下，可以将式(6) AU-F写为式(7):
AhUh=Fh (7) Ah，Uh分别表示在网格尺寸为h^hy (不同网格尺寸对应不同的图像
分辨率)下的方程系数矩阵和解。
假设方程(7)经过n,次迭代后得到的解为tJh ，并令u-为理想的真实解， 则运动矢量的计算误差为
eh=Uh-!Jh (8)
由于计算结果的高频误差经过少量迭代就可以基本消除，因此e-在这 里主要包含计算误差的低频能量。如果能知道eh则可以对tjh进行补偿达到 减少计算误差的目的。虽然直接得到eh并不容易，但可以通过式(9)和式(10) 进行间接求解。通过已知的Ch可以计算出方程(7)的误差R、有式(9)。
Rh=Fh-AhtIh
(9)
由于Ah是线性的因此有
AV=Rh (10)这里的eh主要是在网格h下的低频计算误差。根据前面的分析，知道仅从方程(10)很难得到低频误差eh 。多网格计算的巧妙之处就在于eh不是在原先网格h下计算得到的，而是在更粗的网格上计算得到的。
对于较粗网格H的尺寸为H,xHy，并有H^^、,H,hy，在粗网格上可
以将(10)写为
AHeH=RH (11)其中AH，RH为V，Rh经过尺度映射后得到的新方程系数，映射关系如
下
AH 二 p"h〖Ab) ， RH = p"h〖R"
Ph—h为映射函数，eH反映的是网格H下的高频计算误差，同时又是细网格h下的低频计算误差。由于eH和eh不在一个尺度上，因此还需要将eh映射为eH，如下
在网格H下的高频计算误差eH只需通过很少的迭代就能得到，通过映
射(□)最终可以得到eh ，有了 eh就可以通过(12)来更新Ob 。
tJL=eh+tJh (12)
由于尺度之间的映射也会带来一定的高频误差，因此最后还可以在网格h下再进行&次迭代来消除tJL的高频误差。
可将多网格结构与图像计算中的多尺度思想相结合应用于PIV计算。采用多尺度(有可能尺度因子小于2)逐级逼近的策略，可以降低方程模型误差带来的影响，改善计算精度。每次尺度变化都会进行一次基于当前估计结果的补偿更新。将图像尺度的变化看作是不同尺寸网格的变化，本质上多尺度逼近计算可以看作是多网格算法的一种特例。
因此，根据实际计算需要也可以采用更多尺寸的网格。多网格算法在实际中通常分为三类V型多网格计算、W型多网格计算和完全型多网格计算，图4a和图4b给出了 V型和W型多网格计算的示意图。
优选地，采用如图5所示的多网格结构。图5中的圆点代表不同尺度的网格(从下至上网格尺寸依次增大)，上升箭头表示从小尺度网格向大尺
度网格过渡，向下箭头则相反。每一完成一次v型多网格计算后，便将式
(7)中的方程系数A和F进行更新，进行多次循环处理。通过多网格计算可以有效加快线性方程组迭代的收敛速度，同时不会降低计算的可靠性和精度。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。
权利要求
1.一种用于粒子图像测速中的多网格处理方法，所述粒子图像测速包括求解线性方程组AU＝F以实现对运动场的最佳逼近的步骤，U为待求解，A为系数矩阵，其至少包含运动矢量场参数，F为向量，其特征在于，所述方法包括以下步骤A、设置第一网格，获得方程组AU＝F对应在第一网格h下的方程组AhUh＝FhB、迭代求解方程组AhUh＝Fh，根据所得实际解<overscore>U</overscore>h得到方程误差Rh＝Fh-Ah<overscore>U</overscore>h，C、设置尺寸较第一网格更粗的第二网格H，将在第一网格h下对方程组Aheh＝Rh的求解尺度映射为在第二网格H下对方程组AHeH＝RH的求解；D、迭代求解在第二网格H下的高频计算误差eH；E、对第二网格H下的高频计算误差eH进行尺度映射，得到在第一网格h下的低频计算误差eh；F、用第一网格h下的低频计算误差eh更新实际解<maths id="math0001" num="0001" ><math><![CDATA[ <mrow><msubsup> <mover><mi>U</mi><mo>‾</mo> </mover> <mi>new</mi> <mi>h</mi></msubsup><mo>=</mo><msup> <mi>e</mi> <mi>h</mi></msup><mo>+</mo><msup> <mover><mi>U</mi><mo>‾</mo> </mover> <mi>h</mi></msup><mo>.</mo> </mrow>]]></math> id="icf0001" file="A2009101094290002C1.tif" wi="27" he="4" top= "134" left = "145" img-content="drawing" img-format="tif" orientation="portrait" inline="yes"/></maths>
2. 如权利要求l所述的多网格处理方法，其特征在于，还包括 在更新实际解后，再在第一网格h下迭代求解VlJh-Fh，来消除C^的高频计算误差。
3. 如权利要求1或2所述的粒子图像测速处理方法，其特征在于,循 环多次进行所述处理，其中第二网格的尺度逐次递增。
全文摘要
本发明公开了一种用于粒子图像测速中的多网格处理方法，包括以下步骤设置第一网格，获得待求方程组对应在第一网格下(对应高分辨图像)的方程组；迭代求解方程组，根据所得实际解得到方程误差；设置尺寸较第一网格更粗的第二网格，将在第一网格下对方程组的求解尺度映射为在第二网格下对方程组的求解；迭代求解在第二网格下的高频计算误差；对第二网格下的高频计算误差进行尺度映射，得到在第一网格下的低频计算误差；用第一网格下的低频计算误差更新实际解。本发明的方法能大大加快粒子图像测速的处理过程，同时不会降低处理的可靠性和精度。
文档编号G01P5/20GK101629965SQ200910109429
公开日2010年1月20日 申请日期2009年8月18日 优先权日2009年8月18日
发明者卢宗庆, 廖庆敏 申请人:清华大学深圳研究生院