CN101629966A - 粒子图像测速处理方法 - Google Patents

粒子图像测速处理方法 Download PDF

Info

Publication number
CN101629966A
CN101629966A CN200910109430A CN200910109430A CN101629966A CN 101629966 A CN101629966 A CN 101629966A CN 200910109430 A CN200910109430 A CN 200910109430A CN 200910109430 A CN200910109430 A CN 200910109430A CN 101629966 A CN101629966 A CN 101629966A
Authority
CN
China
Prior art keywords
phi
grid
piv
field
subregion
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN200910109430A
Other languages
English (en)
Other versions
CN101629966B (zh
Inventor
卢宗庆
廖庆敏
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shenzhen Graduate School Tsinghua University
Original Assignee
Shenzhen Graduate School Tsinghua University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shenzhen Graduate School Tsinghua University filed Critical Shenzhen Graduate School Tsinghua University
Priority to CN2009101094306A priority Critical patent/CN101629966B/zh
Publication of CN101629966A publication Critical patent/CN101629966A/zh
Application granted granted Critical
Publication of CN101629966B publication Critical patent/CN101629966B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Image Analysis (AREA)

Abstract

本发明公开了一种粒子图像测速处理方法,包括以下步骤:A.对于N×M像素范围Ω构建水平集函数φ,水平集函数φ为定义在N×M像素范围内的连续函数,用水平集函数φ将区域Ω划分为第一子区域Ω1和第二子区域Ω2,运动场、散度场和旋度场在两个子区域内光滑连续,但在交界处不连续,其中,φ≥0的区域对应第一子区域Ω1,φ<0的区域对应第二子区域Ω2;B.构造区域Ω的能量函数E=E1+λE2,其中E为全局能量泛函,E1为数据约束能量,E2光滑约束能量;C.最小化全局能量泛函,实现运动场、散度场和旋度场的最佳逼近。本发明提供了一种精度高可靠性好的PIV处理方法。

Description

粒子图像测速处理方法
技术领域
本发明涉及图像处理,特别是涉及一种粒子图像测速处理方法。
背景技术
自然界图像的运动形式千差万别,流体运动是一种典型的非刚体运动,流体运动图像计算与分析即粒子图像测速(Particle ImageVelocimitry,简称PIV)是一种新型的非接触式的测量技术。通常PIV是在流体中投入示踪粒子,在激光片光的照射下,在与片光垂直方向用摄像机拍摄随流体运动的粒子图像,再对图像进行处理分析计算,最终得到流场横切面上二维速度场的一种测量方法。作为一种全流场、无接触、无扰动、高精度的流动可视化方法,PIV适用于湍流、非定常流等时变复杂流场的测量。如今,PIV已是一门跨学科的交叉综合技术,其结合了激光技术、视频图像处理技术、计算机技术、流体力学和近代光学技术的最新成果,广泛应用于实验流体力学、生物医学、工业制造等许多领域。PIV将整个运动矢量场作为一个整体连续函数进行估计,其分辨率可以达到图像中的每一个像素,为微观、精确地进行局部图像运动计算和分析提供了可能,同时克服了传统相关方法中的一些固有缺陷。PIV系统在实验流体力学、空气动力学、制造业、医学、工业、飞机制造、水利水电、气象学等领域都有着重要的科学和经济价值。
PIV计算的方法主要分为:基于杨氏干涉条纹法、基于相关(cross-correlation)的方法和粒子跟踪测速(Particle Image Tracing,简称PTV),后两者则是目前PIV计算的主流方法。PIV系统的精度与运动估计算法的性能、实际应用场合和物理器件本身特性有关,抛开后两者就算法而言,传统的算法虽逐渐成熟、实用,但这些方法本身存在一些先天的缺陷。基于相关法的不足归纳如下:
1.由于相关计算时所使用的局部邻域有限,在前后时间内邻域里的粒子出现丢失或出现未知新粒子时,相关计算结果的峰值不一定对应真实的运动位移。
2.根据统计原理,基于相关法得到的运动矢量只是概率上的最大可能位移。
3.当粒子运动的速度很慢时(位移在亚像素级)或有些区域粒子密度过于集中,当采用较大的窗口进行相关搜索时,会出现多峰值的现象。为了克服这个问题常采用尺寸和形状都不固定的局部窗。
与相关方法类似的PTV是通过跟踪单个粒子来完成对运动速度的估计。PTV也存在一些先天的缺陷如:当跟踪过程中粒子消失或别的粒子出现时会带来计算上的不稳定;虽然PTV能得到粒子在一段时间内的运动轨迹(在可靠跟踪的前提下),但PTV对粒子图像相对相关法要求更苛刻,而且PTV不能得到100%密度的运动场,通常需要插值运算。
近年来,基于变分(如光流:物体亮度模式变化的反映,通常指运动矢量场)的PIV计算受到人们越来越多的关注,由于该方法是将整个运动矢量场作为一个整体连续函数进行估计,所以从理论上讲可以克服上面所提到的问题。绝大多数光流计算的分辨率可以达到图像中的每一个像素,这为微观、精确地进行局部图像运动计算和分析提供了可能。同时粒子运动图像大多在时间和空间上都是连续变化的,这与变分方法中图像序列的局部时空可微假设在物理意义上是一致的。光流场计算是通过图像序列获得其运动信息的一种方法,是计算机视觉领域一项重要技术问题。光流不仅携带了被观察物体的运动信息,还携带有被观察物的三维结构、深度、传感器参数、非刚性物体的局部形变,甚至流体运动的矢量结构特征等丰富信息。光流计算在商业领域、科学研究、工业生产中都有广泛的应用,如:视频分割、运动估计、运动目标跟踪、场景的三维结构重建、微观血流运动分析。光流计算的研究真正起始于80年代初Horn和chunck以及Lucas和Kanade等人的奠基性的工作,在其后的30多年时间里出现了许多新方法、新思想。近几年随着偏微分方程(PDE),张量分析,微分几何等数学方法和计算手段在图像分析中的不断渗透,光流计算又有了快速发展,无论在计算精度、可靠性、算法实时性等方面都取得了较大的飞跃。
虽然基于光流技术的PIV计算有着重要的学术和应用价值,但其作为一个新兴的交叉研究领域还存在许多有待解决及改进的地方,归纳如下:
1)对于PIV系统来讲运动矢量计算需要结合实际的流体力学物理意义,如散度和旋度场的数据约束等。另外除矢量场计算外,如何精确计算无散场、无旋场、粘性区域、提取拓扑点(吸收点、发散点)等还存在许多没有解决好的问题。
2)许多变分光流计算需要解决带有边值问题的椭圆偏微分方程,并最终通过求解线性或非线性方程组实现,其需要大量的迭代计算。如何提高光流的计算速度使之达到实时性的要求,一直是一个具有挑战性的课题。
3)对于流体运动图像,除运动矢量本身外矢量的散度和旋度场也是重要的计算与分析对象。运动矢量场的散度和旋度常常被用来计算与分析流体矢量场的结构特征和进行流场可视化(如:分析扩散源、汇聚点、涡旋,以及进行流场的拓扑结构分割等)。通常,会将这些特征单独分析处理再融合判断,而如何基于高维图像的观点对流体运动场进行综合描述和分析仍是一个空白。
发明内容
本发明的主要目的就是针对现有技术的不足,提供一种精度和可靠性高的粒子图像测速处理方法。
本发明还提供一种能进一步有效加快处理过程的粒子图像测速处理方法。
为实现上述目的,本发明采用以下技术方案:
一种粒子图像测速处理方法,包括以下步骤:
A、对于N×M像素范围区域Ω构建水平集函数φ,水平集函数φ为定义在N×M像素范围内的连续函数,用水平集函数φ将区域Ω划分为第一子区域Ω1和第二子区域Ω2,满足运动场、散度场和旋度场在第一子区域Ω1和第二子区域Ω2内光滑连续,但在交界处不连续,其中,φ≥0的区域对应第一子区域Ω1,φ<0的区域对应第二子区域Ω2
其中,第一子区域Ω1对应运动场、散度场和旋度场:[V1,ξ1,η1],第二子区域Ω2区域对应运动场、散度场和旋度场:[V2,ξ2,η2],运动矢量V=[vx,vy],散度 ξ = div ( V ) = ∂ v x x + ∂ v y y , 旋度 η = curl ( V ) = ∂ v x y - ∂ v y x , vx,vy为沿x轴和y轴方向的两个速度分量;
构建指示函数H,定义为:H(x)=1,x≥0 and H(x)=0,x<0;
B、构造区域Ω的能量函数E=E1+λE2,其中E为全局能量泛函,E1为数据约束能量,E2光滑约束能量;
C、最小化全局能量泛函,实现运动场、散度场和旋度场的最佳逼近。
优选地,数据约束能量和光滑约束能量分别采用如下形式:
E 1 = ∫ Ω | I ( X + V 1 , t ) - I ( X , t ) | 2 H ( φ ) + | I ( X + V 2 , t ) - I ( X , t ) | 2 H ( - φ ) dx
+ ∫ Ω | div ( V 1 ) - ξ 1 | 2 H ( φ ) + | div ( V 2 ) - ξ 2 | 2 H ( - φ ) dx
+ ∫ Ω | curl ( V 1 ) - η 1 | 2 H ( φ ) + | curl ( V 2 ) - η 2 | 2 H ( - φ ) dx
E 2 = α ∫ Ω | ▿ V 1 | 2 H ( φ ) + | ▿ V 2 | 2 H ( - φ ) dx
+ α ∫ Ω | ▿ ξ 1 | 2 H ( φ ) + | ▿ ξ 2 | 2 H ( - φ ) dx
+ α ∫ Ω | ▿ η 1 | 2 H ( φ ) + | ▿ η 2 | 2 H ( - φ ) dx
+ β ∫ Ω | ▿ H ( φ ) | 2 dx
优选地,全局能量泛函的最小化通过求解以下线性方程组实现:
AU=F
U=[V1,V2,ξ1,ξ2,η1,η2,φ],A为系数矩阵,F为向量,表示线性方程组的奇次项。
可以通过经典的变分计算方法获得和处理线性方程组AU=F。
优选通过多网格方法处理线性方程组AU=F。
所述处理优选地包括:
设置第一网格,获得方程组AU=F对应在第一网格h下的方程组AhUh=Fh
迭代求解方程组AhUh=Fh,根据所得实际解U
Figure G2009101094306D000410
得到方程误差
Figure G2009101094306D00049
设置尺寸较第一网格更粗的第二网格H,将在第一网格h下对方程组Aheh=Rh的求解尺度映射为在第二网格H下对方程组AHeH=RH的求解;
迭代求解在第二网格H下的高频计算误差eH
对第二网格H下的高频计算误差eH进行尺度映射,得到在第一网格h下的低频计算误差eh;和
用第一网格h下的低频计算误差eh更新实际解: U ‾ new h = e h + U ‾ h .
所述处理还可以包括:
在更新实际解后,再在第一网格h下迭代求解AhUh=Fh,来消除
Figure G2009101094306D000411
的高频计算误差。
优选地,循环多次进行所述处理,其中第二网格的尺度逐次递增。
本发明有益的技术效果是:
与灰度图像类似,矢量场的结构信息表现为局部矢量场的边缘、多矢量交汇、分叉等结构特点,即可以看作运动估计中的多运动和运动边界。本发明粒子图像测速处理方法,对于2维PIV运动矢量场、散度和旋度场获得,采用了基于光流运动方程的计算思想,以整体函数的形式计算运动矢量,并通过构建水平集(level set)函数,将水平集的策略和思想集成于现有的PIV计算框架,解决了运动边界和多运动问题,同时提高了流体运动特征描述的精度,在抑止噪声同时又能保持矢量场、散度和旋度场的局部结构信息,精度高,可靠性好。
进一步地,在图像计算中采用多网格方法进行处理,可以有效加快线性方程组迭代的收敛速度,大大加快图像处理速度,同时不会降低图像处理的可靠性和精度。
附图说明
图1为粒子图像测速处理系统的结构框图;
图2为本发明粒子图像测速处理方法一种实施例的流程图;
图3为用水平集函数φ来划分内部光滑连续而边界不连续的场的示意图;
图4为一种实施例中利用多网格进行加速处理的流程图;
图5a和图5b分别为V型和W型多网格结构的示意图;
图6为一种优选实施例的多网格结构的示意图。
本发明的特征及优点将通过实施例结合附图进行详细说明。
具体实施方式
一个完整的粒子图像测速处理系统其关键部分如图1所示,主要包括数据采集模块、图像预处理模块、运动矢量场(光流场)计算模块和后续数据分析模块,后续数据分析模块包括旋度计算模块和散度计算模块。其中的运动矢量场计算模块是整个PIV系统的核心模块。
一种粒子图像测速处理方法请参考图2。在计算2维运动矢量场、散度和旋度场时,采用了基于光流运动方程的计算思想,以整体函数的形式计算运动矢量。为了能够抑止噪声同时又能保持矢量场、散度和旋度场的局部结构信息,还结合了水平集的思想。
2维运动矢量场是指在大小为N×M像素范围内(如通常熟知的PAL制:768×576或640×480个像素)能够描述每一像素位置(N×M内)下该像素当前时刻的位移矢量,该矢量为一浮点数二维矢量,表示为:V=[vx,vy],其中vx,vy为x轴和y轴方向的两个速度分量。为了统一描述,用大写粗体表示矢量,小写细体表示标量。速度场v实际是坐标位置变量x=[x,y]的函数。流体运动图像计算出了获得瞬时运动矢量场v外,还需计算v的散度场ξ和旋度场η,如(1)(2)式所示。
ξ = div ( V ) = ∂ v x x + ∂ v y y - - - ( 1 )
η = curl ( V ) = ∂ v x y - ∂ v y x - - - ( 2 )
由于实际的流体运动可能会有多矢量流交界的现象,即运动边界和多运动问题,为了解决运动边界和多运动问题,同时为了提高流体运动特征描述的精度,将水平集的策略和思想集成于原有的光流计算框架。
如图3所示,在某一区域Ω内运动矢量、散度和旋度有可能不连续,为此可以将Ω划分为两个子区域Ω1,Ω2,Ω1区域内对应的运动场、散度场和旋度场为:[V1,ξ1,η1];Ω2区域内对应的运动场、散度场和旋度场为:[V2,ξ2,η2]。假设在区域Ω1内[V1,ξ1,η1]光滑连续,在区域Ω2内[V2,ξ2,η2]光滑连续,但在Ω1,Ω2的交界处运动场、散度场和旋度场并不连续。为了保证计算时[V1,ξ1,η1]和[V2,ξ2,η2]的连续性,同时又不破坏Ω1与Ω2之间的结构信息以及便于计算,构建水平集函数φ和指示函数H,φ为定义在图像大小范围内(N×M)的连续函数。H(x)=1,x≥0 and H(x)=0,x<0。
用水平集函数φ来划分Ω1,Ω2,φ≥0的区域属于Ω1,而φ<0的区域属于Ω2。由于采用了水平集的思想,因此在计算[V1,ξ1,η1]和[V2,ξ2,η2]时同时要不断计算更新φ。
为了估计出[V1,ξ1,η1]和[V2,ξ2,η2],构造以下能量函数。
E=E1+λE2    (3)
其中E为全局能量泛函,E1,E2分别为以下式(4)、(5)形式的数据约束能量和光滑约束能量。
E 1 = ∫ Ω | I ( X + V 1 , t ) - I ( X , t ) | 2 H ( φ ) + | I ( X + V 2 , t ) - I ( X , t ) | 2 H ( - φ ) dx
+ ∫ Ω | div ( V 1 ) - ξ 1 | 2 H ( φ ) + | div ( V 2 ) - ξ 2 | 2 H ( - φ ) dx - - - ( 4 )
+ ∫ Ω | curl ( V 1 ) - η 1 | 2 H ( φ ) + | curl ( V 2 ) - η 2 | 2 H ( - φ ) dx
E 2 = α ∫ Ω | ▿ V 1 | 2 H ( φ ) + | ▿ V 2 | 2 H ( - φ ) dx
+ α ∫ Ω | ▿ ξ 1 | 2 H ( φ ) + | ▿ ξ 2 | 2 H ( - φ ) dx (5)
+ α ∫ Ω | ▿ η 1 | 2 H ( φ ) + | ▿ η 2 | 2 H ( - φ ) dx
+ β ∫ Ω | ▿ H ( φ ) | 2 dx
这里I(X,t)代表粒子图像序列。最终估计的光流场可以描述为:
V = V 1 φ ( X ) ≥ 0 V 2 else
通过全局能量泛函的最小化来实现光流场的最佳逼近。
能量泛函(3)的最小化可以通过求解线性方程组(6)实现,。
AU=F    (6)
其中U=[V1,V2,ξ1,ξ2,η1,η2,φ],A为9×9的系数矩阵,F为9×1的向量表示线性方程的奇次项。
线性方程组(6)的获得和计算可以通过经典的变分计算方法实现,这里不详细阐述。
在优选的实施例中,可以将基于多网格加速策略与实际的PIV计算相结合,通过多网格方法来处理线性方程组(6),加速PIV计算。
虽然光流计算无论在计算精度、可靠性等方面都取得了较大的飞跃,但其数值方法和手段并没有重大突破。已在工程数学和有限元数值计算中广泛使用的多网格方法,凭借其优越的加速收敛特性开始得到图像处理领域人员的关注。
多网格计算的一个中心任务是求解各种各样来自实际问题的偏微分方程及线性方程组。由于许多图像分析问题需要解决带有边值问题的椭圆偏微分方程,并最终通过求解线性或非线性方程组实现,如:基于变分方法的图像非线性扩散、主动轮廓计算、图像恢复等,在大尺寸网格上线性方程的高频误差很容易消除,相反即使经过大量迭代低频误差却很难消除,但在小尺度网格上这种低频误差只需经少量迭代就能明显消除。因此多网格策略可通过不同网格的迭代计算来加速收敛过程。对于图像分析而言,网格的定义比较简单,原始图像的每一个像素表示一个网格,即最细的网格,计算都是基于这样的离散网格点上的值运算完成的。
请参考图4,一种优选的实施例的粒子图像测速处理还包括如下步骤:
可以将式(6)写为式(7):
AhUh=Fh    (7)
其中指标h代表了网格尺寸为hx×hy的第一网格,Ah,Uh分别表示在网格尺寸为hx×hy(不同网格尺寸对应不同的图像分辨率)下的方程系数矩阵和解。
假设方程(7)经过n1次迭代后得到的解为
Figure G2009101094306D00082
,并令Uh为理想的真实解,则运动矢量的计算误差为:
Figure G2009101094306D00083
由于计算结果的高频误差经过少量迭代就可以基本消除,因此eh在这里主要包含计算误差的低频能量。如果能知道eh则可以对
Figure G2009101094306D00082
进行补偿达到减少计算误差的目的。虽然直接得到eh并不容易,但可以通过式(9)和式(10)进行间接求解。
通过已知的
Figure G2009101094306D00082
我们可以计算出方程(7)的误差Rh,有式(9)。
Figure G2009101094306D00084
由于Ah是线性的因此有:
Aheh=Rh    (10)
这里的eh主要是在网格h下的低频计算误差。根据前面的分析,知道仅从方程(10)很难得到低频误差eh。多网格计算的巧妙之处就在于eh不是在原先网格h下计算得到的,而是在更粗的网格上计算得到的。
对于较粗网格H的尺寸为Hx×Hy,并有Hx>hx,Hy>hy,在粗网格上可以将(10)写为:
AHeH=RH    (11)
其中AH,RH为Ah,Rh经过尺度映射后得到的新方程系数,映射关系如下:
AH=Ph→H(Ah),RH=Ph→H(Rh)
Ph→H为映射函数,eH反映的是网格H下的高频计算误差,同时又是细网格h下的低频计算误差。由于eH和eh不在一个尺度上,因此还需要将eh映射为eH,如下:
eh=Ph→H(eH)
在网格H下的高频计算误差eH只需通过很少的迭代就能得到,通过映射PH→h(□)最终可以得到eh,有了eh就可以通过(12)来更新
Figure G2009101094306D00082
U ‾ new h = e h + U ‾ h - - - ( 12 )
由于尺度之间的映射也会带来一定的高频误差,因此最后还可以在网格h下再进行n2次迭代来消除Unew h的高频误差。
可将多网格结构与图像计算中的多尺度思想相结合应用于PIV计算。采用多尺度(有可能尺度因子小于2)逐级逼近的策略,可以降低方程模型误差带来的影响,改善计算精度。每次尺度变化都会进行一次基于当前估计结果的补偿更新。将图像尺度的变化看作是不同尺寸网格的变化,本质上多尺度逼近计算可以看作是多网格算法的一种特例。
因此,根据实际计算需要也可以采用更多尺寸的网格。多网格算法在实际中通常分为三类:V型多网格计算、W型多网格计算和完全型多网格计算,图5a和图5b给出了V型和W型多网格计算的示意图。
优选地,采用如图6所示的多网格结构。图6中的圆点代表不同尺度的网格(从下至上网格尺寸依次增大),上升箭头表示从小尺度网格向大尺度网格过渡,向下箭头则相反。每一完成一次V型多网格计算后,便将式(7)中的方程系数A和F进行更新,进行多次循环处理。通过多网格计算可以有效加快线性方程组迭代的收敛速度,同时不会降低计算的可靠性和精度。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明,不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说,在不脱离本发明构思的前提下,还可以做出若干简单推演或替换,都应当视为属于本发明的保护范围。

Claims (8)

1.一种粒子图像测速处理方法,其特征在于包括以下步骤:
A、对于N×M像素范围区域Ω,构建水平集函数φ,水平集函数φ为定义在N×M像素范围内的连续函数,用水平集函数φ将区域Ω划分为第一子区域Ω1和第二子区域Ω2,满足运动场、散度场和旋度场在第一子区域Ω1和第二子区域Ω2内光滑连续,但在交界处不连续,其中,φ≥0的区域对应第一子区域Ω1,φ<0的区域对应第二子区域Ω2
其中,第一子区域Ω1对应运动场、散度场和旋度场为[V1,ξ1,η1],第二子区域Ω2区域对应运动场、散度场和旋度场:[V2,ξ2,η2],运动矢量V=[vx,vy],散度 ξ = div ( V ) = ∂ v x x + ∂ v y y , 旋度 η = curl ( V ) = ∂ v x y - ∂ v y x , vx,vy为沿x轴和y轴方向的两个速度分量;
构建指示函数H,定义为:H(x)=1,x≥0 and H(x)=0,x<0;
B、构造区域Ω的能量函数E=E1+λE2,其中E为全局能量泛函,E1为数据约束能量,E2光滑约束能量;
C、最小化全局能量泛函,实现运动场、散度场和旋度场的最佳逼近。
2.如权利要求1所述的粒子图像测速处理方法,其特征在于,数据约束能量和光滑约束能量分别采用如下形式:
E 1 = ∫ Ω | I ( X + V 1 , t ) - I ( X , t ) | 2 H ( φ ) + | I ( X + V 2 , t ) - I ( X , t ) | 2 H ( - φ ) dx
+ ∫ Ω | div ( V 1 ) - ξ 1 | 2 H ( φ ) + | div ( V 2 ) - ξ 2 | 2 H ( - φ ) dx
+ ∫ Ω | curl ( V 1 ) - η 1 | 2 H ( φ ) + | curl ( V 2 ) - η 2 | 2 H ( - φ ) dx
E 2 = α ∫ Ω | ▿ V 1 | 2 H ( φ ) + | ▿ V 2 | 2 H ( - φ ) dx
+ α ∫ Ω | ▿ ξ 1 | 2 H ( φ ) + | ▿ ξ 2 | 2 H ( - φ ) dx
+ α ∫ Ω | ▿ η 1 | 2 H ( φ ) + | ▿ η 2 | 2 H ( - φ ) dx
+ β ∫ Ω | ▿ H ( φ ) | 2 dx .
3.如权利要求1所述的粒子图像测速处理方法,其特征在于,全局能量泛函的最小化通过求解以下线性方程组实现:
AU=F
U=[V1,V2,ξ1,ξ2,η1,η2,φ],A为系数矩阵,F为向量,表示线性方程组的奇次项。
4.如权利要求3所述的粒子图像测速处理方法,其特征在于,通过经典的变分计算方法获得和处理线性方程组AU=F。
5.如权利要求3所述的粒子图像测速处理方法,其特征在于,通过多网格方法处理线性方程组AU=F。
6.如权利要求5所述的粒子图像测速处理方法,其特征在于,所述处理包括:
设置分辨率较高的第一网格,获得方程组AU=F对应在第一网格h下的方程组AhUh=Fh
迭代求解方程组AhUh=Fh,根据所得实际解
Figure A2009101094300003C3
得到方程误差
Figure A2009101094300003C2
设置尺寸较第一网格更粗的第二网格H,将在第一网格h下对方程组Aheh=Rh的求解尺度映射为在第二网格H下对方程组AHeH=RH的求解;
迭代求解在第二网格H下的高频计算误差eH
对第二网格H下的高频计算误差eH进行尺度映射,得到在第一网格h下的低频计算误差eh;和
用第一网格h下的低频计算误差eh更新实际解: U ‾ new h = e h + U ‾ h .
7.如权利要求6所述的粒子图像测速处理方法,其特征在于,所述处理还包括:
在更新实际解后,再在第一网格h下迭代求解AhUh=Fh,来消除
Figure A2009101094300003C4
的高频计算误差。
8.如权利要求6或7所述的粒子图像测速处理方法,其特征在于,循环多次进行所述处理,其中第二网格的尺度逐次递增。
CN2009101094306A 2009-08-18 2009-08-18 粒子图像测速处理方法 Expired - Fee Related CN101629966B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN2009101094306A CN101629966B (zh) 2009-08-18 2009-08-18 粒子图像测速处理方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN2009101094306A CN101629966B (zh) 2009-08-18 2009-08-18 粒子图像测速处理方法

Publications (2)

Publication Number Publication Date
CN101629966A true CN101629966A (zh) 2010-01-20
CN101629966B CN101629966B (zh) 2011-09-28

Family

ID=41575144

Family Applications (1)

Application Number Title Priority Date Filing Date
CN2009101094306A Expired - Fee Related CN101629966B (zh) 2009-08-18 2009-08-18 粒子图像测速处理方法

Country Status (1)

Country Link
CN (1) CN101629966B (zh)

Cited By (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102831616A (zh) * 2012-07-10 2012-12-19 华东师范大学 一种视频流体运动矢量计算方法
CN104680564A (zh) * 2015-03-12 2015-06-03 北京航空航天大学 一种灰度增强层析piv重构方法、装置和设备
CN104777329A (zh) * 2014-01-13 2015-07-15 北京航空航天大学 一种用于粒子图像测速三维粒子场重构的线性规划算法
CN104915928A (zh) * 2014-03-12 2015-09-16 北京航空航天大学 一种基于本征正交分解的速度场坏矢量识别和修正方法
CN105139423A (zh) * 2015-08-18 2015-12-09 上海交通大学 基于旋度和散度的运动特征提取方法及系统
CN105807093A (zh) * 2016-04-18 2016-07-27 北京航空航天大学 一种基于粒子图像测速技术的加速度测量方法及装置
CN106097396A (zh) * 2016-06-29 2016-11-09 华东师范大学 一种基于流体运动矢量场的并行分析方法
CN103618904B (zh) * 2013-11-20 2017-02-22 华为技术有限公司 基于像素的运动估计方法及装置
CN110223328A (zh) * 2019-05-14 2019-09-10 焦作大学 一种基于物理学的改善粒子图像测速稳健性光流方法
CN111369620A (zh) * 2019-12-23 2020-07-03 东北石油大学 一种基于改进piv的水平油水两相流流速测量方法
CN111398625A (zh) * 2020-03-19 2020-07-10 西安理工大学 一种物理模型试验中的测速方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN2606895Y (zh) * 2003-03-19 2004-03-17 申功炘 数字式粒子图像测速系统
CN1588092A (zh) * 2004-08-18 2005-03-02 浙江大学 微观流场粒子图像显微测速系统
JP2008199477A (ja) * 2007-02-15 2008-08-28 Matsushita Electric Ind Co Ltd 撮像装置
CN201037868Y (zh) * 2007-03-13 2008-03-19 魏润杰 一种高速粒子图像测速装置

Cited By (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102831616B (zh) * 2012-07-10 2014-12-10 华东师范大学 一种视频流体运动矢量计算方法
CN102831616A (zh) * 2012-07-10 2012-12-19 华东师范大学 一种视频流体运动矢量计算方法
CN103618904B (zh) * 2013-11-20 2017-02-22 华为技术有限公司 基于像素的运动估计方法及装置
CN104777329B (zh) * 2014-01-13 2018-06-05 北京航空航天大学 一种用于粒子图像测速三维粒子场重构的线性规划算法
CN104777329A (zh) * 2014-01-13 2015-07-15 北京航空航天大学 一种用于粒子图像测速三维粒子场重构的线性规划算法
CN104915928A (zh) * 2014-03-12 2015-09-16 北京航空航天大学 一种基于本征正交分解的速度场坏矢量识别和修正方法
CN104915928B (zh) * 2014-03-12 2018-02-06 北京航空航天大学 一种基于本征正交分解的速度场坏矢量识别和修正方法
CN104680564A (zh) * 2015-03-12 2015-06-03 北京航空航天大学 一种灰度增强层析piv重构方法、装置和设备
CN104680564B (zh) * 2015-03-12 2017-11-10 北京航空航天大学 一种灰度增强层析piv重构方法、装置和设备
CN105139423A (zh) * 2015-08-18 2015-12-09 上海交通大学 基于旋度和散度的运动特征提取方法及系统
CN105139423B (zh) * 2015-08-18 2018-04-20 上海交通大学 基于旋度和散度的运动特征提取方法及系统
CN105807093A (zh) * 2016-04-18 2016-07-27 北京航空航天大学 一种基于粒子图像测速技术的加速度测量方法及装置
CN106097396A (zh) * 2016-06-29 2016-11-09 华东师范大学 一种基于流体运动矢量场的并行分析方法
CN110223328A (zh) * 2019-05-14 2019-09-10 焦作大学 一种基于物理学的改善粒子图像测速稳健性光流方法
CN111369620A (zh) * 2019-12-23 2020-07-03 东北石油大学 一种基于改进piv的水平油水两相流流速测量方法
CN111398625A (zh) * 2020-03-19 2020-07-10 西安理工大学 一种物理模型试验中的测速方法
CN111398625B (zh) * 2020-03-19 2022-04-12 西安理工大学 一种物理模型试验中的测速方法

Also Published As

Publication number Publication date
CN101629966B (zh) 2011-09-28

Similar Documents

Publication Publication Date Title
CN101629966B (zh) 粒子图像测速处理方法
CN108010103A (zh) 复杂河道地形快速精细生成方法
Weinkauf et al. Stable feature flow fields
Ferstl et al. Interactive separating streak surfaces
CN109308308A (zh) 基于三维动态可视化排水管网模拟和结果分析方法及装置
CN106886980A (zh) 一种基于三维激光雷达目标识别的点云密度增强的方法
CN103700117A (zh) 一种基于tv-l1变分模型的鲁棒光流场估计方法
CN101629965B (zh) 用于粒子图像测速中的多网格处理方法
Suvorov Wave-optical effects in the microlensing of continuous gravitational waves by star clusters
CN103514600A (zh) 一种基于稀疏表示的红外目标快速鲁棒跟踪方法
CN104091352A (zh) 基于结构相似度的视觉跟踪方法
CN112834432A (zh) 一种基于遥感技术与运动学原理的滑坡厚度反演方法
Wang et al. Remote sensing image semantic segmentation method based on small target and edge feature enhancement
Ma et al. A comprehensive deep learning geometric shape optimization framework with field prediction surrogate and reinforcement learning
Chen et al. Motion estimation for complex fluid flows using helmholtz decomposition
Liao et al. Physics-based optical flow estimation under varying illumination conditions
Li et al. Research on PIV image matching algorithm of turbulent velocity field based on circular projection
CN110555826B (zh) 一种基于局部离群因子的三维点云特征提取方法
CN102323965A (zh) 面向非定常三维流场涡结构的智能分析方法
Zhang et al. Applicability of deep learning optical flow estimation for PIV methods
CN106846481A (zh) 一种地质剖面图的生成方法
Pan et al. Embedded U-Net: combines multiple feature fusion encode and subpixel reconstruction for microcracks salient object detection
CN111967151A (zh) 一种通过添加k和ε耦合源项来校正标准k-ε模型的方法
Ying et al. Methods of target motion estimation for AUV target tracking
CN103413341B (zh) 一种基于Gerstner模型的波浪折射模拟方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20110928

Termination date: 20150818

EXPY Termination of patent right or utility model