利用沃罗诺伊划分的粒子追踪测速算法及其在大速度梯度流场中的应用

吴长松,白旭峰,张洋,周立飞

(1.西安交通大学能源与动力工程学院,710049,西安;2.西安陕鼓动力股份有限公司,710075,西安)

随着图像识别技术的提高和精细测量的需要,粒子追踪测速(PTV)算法在近年来得到了不断的进步[1]。在需要观测粒子运动轨迹的情况下[2],PTV则直接跟踪流场中的粒子运动,从而避免了PIV方法的平均效应[3]。在众多PTV算法中,有一类算法基于粒子丛思想,即以目标粒子与其周围粒子的相对位置关系为判断基准来实现前后帧粒子匹配。粒子丛算法预设参数少,因此运行速度较快,而在构建了合理的粒子丛特征后,此类算法在粒子丢失或图像噪音等问题程度较轻时能够获得可观的准确率。典型的粒子丛类算法有Okamoto等提出的弹簧模型算法[4],以及Ishikawa等提出的速度梯度张量算法[5],这两种算法的问题在于粒子丛的大小需要提前规定。因此,当临近粒子所组成粒子丛的大小随流动发生改变时,算法的准确率骤降,而采用适当的网格划分技术则可以避免此问题。Song等提出的基于德劳内划分的算法(DT-PTV)将图像中离散粒子分布转化成为三角形非结构化网格,通过匹配前后两帧的三角形单元来实现粒子的匹配[6],张洋等通过引入逻辑判断模块实现了DT-PTV运算结果中错误矢量的有效剔除[7],贾攀等将这种方法推广到了三维情况,即利用四面体单元的信息进行粒子匹配[8]。然而,上述算法存在两个问题:首先,由于三角形或四面体结构简单,在粒子浓度较高时会出现形状趋同,从而导致匹配难度增加;其次,无论对于二维还是三维问题,利用德劳内划分进行粒子匹配时,需要对三角形或四面体进行表层匹配,然后再通过三角形或四面体的顶点配对来实现里层的粒子匹配,这不仅造成算法结构的复杂化,而且在粒子分布相对均匀时将增加里层匹配的逻辑甄别难度。本文以沃罗诺伊划分为基础,提出了更为简洁高效的粒子丛类PTV算法,尝试构建以原始粒子为中心的匹配单元,从而仅通过一层匹配完成运算,同时引入人工模拟流场和实际复杂流场来检验算法的准确性和实用性。

1 算法准备

1.1 确定候选粒子

粒子丛类算法利用目标粒子与其周围粒子的相对位置关系构建匹配单元。第1帧中的每个目标粒子需要从第2帧挑选出自己的候选粒子,原则如下

|Xn-Ym|

(1)

式中:Xn、Ym分别为第1、2帧粒子的矢量坐标;Rs为查问域半径,其取值需要大于两帧间粒子最大位移,而后者可以通过实际流场的主流速度进行估算。构建与粒子一一对应的匹配单元,然后将目标粒子的匹配单元与候选粒子的匹配单元逐一比较,最匹配单元所对应的候选粒子则被定义为匹配粒子。

1.2 两种对偶划分

对于二维平面中的离散点,沃罗诺伊划分是将平面区域离散点群落转化为互不重叠的多边形,而德劳内划分是将平面区域离散点连接成为成互不重叠的三角形,如图1所示。在构建沃罗诺伊网格时,常先将离散点构成德劳内三角网格,再将德劳内三角形各边的中垂线段进行连接可得到沃罗诺伊划分,因此沃罗诺伊网格和德劳内网格互为对偶,而原始粒子是沃罗诺伊多边形的中心[9]。在三维空间中,沃罗诺伊划分(多面体)和三维德劳内划分(四面体)同样互为对偶。两种结构在运用时各有好处,可根据实际情况选择合适的结构,德劳内网格曾被用于改进松弛迭代类PTV算法中用来圈选匹配粒子的模块[10]。

(a)二维平面粒子点 (b)沃罗诺伊多边形图1 粒子群的沃罗诺伊划分

2 算法原理

2.1 二维算法

图1中原始粒子用“+”符号表示,经划分之后得到的封闭区域被称为沃罗诺伊多边形,它们紧密排列且绝对没有交叠,并与原始粒子呈一一对应的关系,即每一个沃罗诺伊多边形都拥有一个核。对参与匹配的前后两帧粒子群分别作沃罗诺伊划分,基于两帧流动相关性假设,同一粒子在前后两帧中对应的沃罗诺伊多边形变形有限。利用这一特性,将沃罗诺伊多边形定义为基本匹配单元。

(1)提取多边形即匹配单元的特征。二维算法中的匹配单元如图2所示,以原始粒子(x0,y0)为中心建立笛卡尔坐标系,它所对应的匹配单元由若干个三角形Ⅰ,Ⅱ,…,Ⅴ组成,三角形的极半径r与极角α的关系为

(2)

最后,此多边形就被转换成特征曲线,结果如图3所示。

图2 二维算法中的匹配单元

图3 沃罗诺伊多边形特征曲线

(3)

式中:Γ为协方差运算;D为方差运算。考虑到多边形跨帧过程中有可能进行旋转,需要比较两条曲线在0°～360°不同相位差下的相似性,本文选择的相位变化步长为2°。根据式(3),相比较的两条曲线长度相同,因此当相位差不为零时,仅取两条曲线相位重叠的部分进行比较,任意一对特征曲线会计算出多个相关系数,选择最大值作为最终相关系数。最后,所有候选者中对应最大相关系数的粒子被认定为匹配粒子。

考虑到候选粒子中真实匹配粒子已经丢失的情况,例如粒子在第2帧中逃逸至图像之外,只有足够相似的匹配才会被认为有效,若某一粒子与其所有候选粒子的相关系数均小于0.5,则认为此粒子在第2帧中没有匹配。由于实际流场变形方式及强度均不相同,0.5实际上对应均匀概率,即两个多边形完全相同与完全不同两种极端情况之间的中间点。对于变形较强的流动,可适当放宽匹配条件即相关性阈值小于0.5,并基于匹配结果的人眼判断返回修正这一阈值,直至达到满意结果,反之亦然。

(3)遍历第1帧所有粒子,完成整场粒子匹配,基于匹配关系和时间间隔可得向前差分格式的瞬时速度场,即实现追踪测速,而二维算法仅有Rs这一个独立预设参数。

2.2 三维算法

(a)空间中的粒子n (b)粒子n的参考四面体图4 三维算法中的匹配单元

对三维空间离散的粒子点,采用三维德劳内划分将其转化为空间四面体群,此时名义上的匹配单元为所有包含同一粒子n的四面体所聚合而成的多面体,或者此多面体可以被划分成若干个以粒子n为共同顶点的四面体,如图4中所示。但是,实际操作中的基本匹配单元则是这些四面体,如此拆分的原因在于多面体的参数化远比二维多边形的要困难。此处引入一种全新的四面体投票机制,所涉及的多面体匹配与三维DT-PTV中的表层匹配有根本不同,这些作为基本匹配单元的四面体被称为粒子的参考四面体。如果一个四面体单元可在空间中做任意旋转且具有一个基点,那么它一共有6个自由度。因此,选取6个线性无关量作为参考四面体的特征量,包括从粒子n指向其余3个顶点的向量的长度,以及代表这3个向量间夹角的点积。将上述特征量汇总为如下两个四面体特征向量

(4)

参考四面体之间的差异dij,计算式为

dij=‖lni-lmj‖1+‖ani-amj‖1

(5)

式中:下标n和m分别为第1和第2帧中的粒子;下标i和j分别为第1帧和第2帧中粒子的参考四面体;‖‖1为向量的1范数。dim为dij的最小值,取得最小值的参考四面体为j′,对应的四面体特征向量为lmj′和amj′。对dim进行无量纲化处理,即

(‖lni‖1+‖ani‖1+‖lmj′‖1+‖amj′‖1)-1

(6)

为了避免多面体局部大变形或者粒子缺失所带来的影响,引入参考四面体投票机制来确定匹配粒子。粒子n的全部参考四面体即为选民,粒子n的候选粒子即为候选人。选民根据如下公式来决定投票与否,获得最多认可的候选人当选,即

(7)

vm/Ni>β

(8)

式中:β∈[0,1]为投票系数。如果式(8)成立,对粒子n的所有候选粒子进行相同操作,否则认为本次投票结果无效,此时vm为0。vm值最大的候选粒子即为粒子n在第2帧中的匹配粒子,若vm的最大值为0,则认为粒子n在第2帧粒子中没有匹配结果。

α代表两个参考四面体间的差异,在算法中决定是否能投票,β代表两个粒子之间相似性的预判,在算法中决定投票最终是否有效。如果α增大,根据式(7)可知,投票的可能性增大,即总票数增加,此时β也相应增大,从而提高了投票的有效性,使算法在逻辑上更加合理,反之亦然。由式(7)和vm的定义可知,α增加一个较小量就可能会导致总票数有较大增长,即β有较大增长,β的增速应大于α的。本文通过人工流场测试证明α∈(0,0.25]和β=2α是较为有效的选择,三维算法仅有Rs和α两个独立的预设参数。

3 人工流场测试

(a)人工流场速度矢量图

(b)两种算法的准确率随查问域半径的变化图5 人工流场及测试结果

以二维算法为例进行人工流场测试。通过叠加流动解析解生成复杂人工流场。首先在第1帧随机布置粒子,然后赋予一定的流动规则,从而在第2帧形成新的粒子位置,并使得前后两帧联合描述一个瞬时流动。流动规则分多步实施:首先形成具有一定速度梯度的剪切流,然后设置若干旋涡和偶极子,最后增大流场速度和粒子浓度,并添加相对误差为5%的位置扰动,以及随机抹除第2帧数量比例为5%的粒子。流场大小为256×256像素,旋涡区大小为128×128像素,偶极子的距离为80像素,两帧之间最大位移为10像素,由于人工流场真实匹配结果已知,任一种PTV的匹配结果与其比较之后即可得到计算准确率,人工流场及测试结果如图5所示。随着Rs的增大,新二维算法的准确率更加稳定,这是因为所有的三角形在拓扑结构上是相同的,Rs值的增大将包纳更多形状相似、但实际上毫无关系的候选匹配单元,这将对DT-PTV的表层匹配产生较大干扰,而多边形的个体差异远大于三角形之间的差异,所以能够更好抵抗Rs的变化所带来的影响。因此,采用一个足够大的固定Rs值可使算法满足一个图像序列中所有瞬时或者一个瞬时中所有粒子的匹配需要,这两种情况分别对应流动的非定常性和非均匀性。

4 实验测试

4.1 二维算法测试

采用脉冲流检验算法,流动通过模拟腹主动脉瘤中血液流动实现,流动循环系统如图6所示。

图6 人工腹主动脉瘤中脉冲流动循环系统

系统主要包括位于测试段的腹主动脉瘤模型、一个离心泵以及两个电磁阀。其中离心泵用来给流体提供动力,两个电磁阀用来模拟周期约为1 s的人体血流脉冲波形,电磁流量计用来测试主流速度。腹主动脉瘤模型由透明聚甲基丙烯酸甲酯制成,人造血液为45%蒸馏水加55%甘油混合而成,并添加了一定浓度的荧光示踪粒子。通过高速激光器和狭缝装置在测试段的流向中心纵剖面上形成片光,然后采用高速相机以最高100 Hz的频率在侧面记录粒子图像。实验共记录20个脉冲周期的图像,选择其中一个周期测试二维算法,每个周期包含20对图像,每对图像对应一个瞬时流场。记人造血液的流动方向为x,铅垂方向为z。用动态阈值算法[11]识别图像中的示踪粒子并将其转化为单像素点。模型总长L=80 mm,选择固定查问域半径为L/5,这已超过主流速度所对应的两帧间最大位移,选取2个时刻进行分析,结果与文献[12]的分析结果一致,二维算法对人工脉冲流动的重构如图7所示。新算法可以在选取较大Rs值的条件下应对流动序列的非定常性以及瞬时流场的非均匀性。

(a)模拟脉冲流动波形

(b)时刻1的流场重构

(c)时刻2的流场重构图7 二维算法对人工脉冲流动的重构

4.2 三维算法测试

本实验在循环水洞中进行。在水中播撒直径为10 μm的荧光粒子作为示踪粒子,在测试段上游布置电磁流量计,测得主流速度约为0.16 m/s。尾迹流实验段示意图如图8所示,图中透明测试段长约200 cm,截面尺寸为30×30 cm2。V3V系统[13]安装在测试段的顶部,通过3个镜头重构示踪粒子的瞬时三维坐标,拍摄频率为15 Hz。V3V的测量区域为一长方体,其厚度为7 cm,截面积为20 cm×20 cm,距离水洞底部10 cm。拍摄区域上游10 cm处放置一个高为30 cm、截面积为6 cm×6 cm的长方体障碍物以生成尾迹流。

图8 尾迹流实验段示意图

(a)流场全局图

(b)流场俯视图

将V3V系统计算得到的粒子三维坐标导入新发展的三维算法中,预设参数α=β/2=0.05,综合主流速度和采集频率信息设置查问域半径为1.2 cm,通过三维算法重构得到的瞬时流场如图9所示,图中x方向为主流方向。由图9可知:算法成功重构了速度梯度较大的区域,即势流-尾迹流结合区域以及尾迹内部的旋涡结构;图9c所示的少量错误匹配可通过基于逻辑判断的后处理方法移除[7,14-15]。目前,V3V系统的后处理模块采用基于松弛迭代的粒子追踪算法进行粒子匹配,其预设参数较多,用户友好性一般[11],当前的三维流场测试结果表明,预设参数少、结构简洁的新三维算法亦可用来与V3V等三维成像硬件系统搭配使用。

5 结 论

本文提出了一种基于沃罗诺伊划分的粒子丛类PTV算法,并针对匹配单元的构建原则提出了相应的二维及三维算法。二维算法充分利用了多边形形状的特异性,是一种结构简单、稳定性与准确率兼得的算法;三维算法则采用沃罗诺伊的对偶划分来构建最小匹配单元,从而规避了三维多面体自由度过多、难以参数化的问题。人工脉冲流场和尾迹流流场的测试结果表明,新发展的算法能很好处理速度梯度较大的复杂流场,具有良好的实用价值。