基于全局位置迭代PCSS算法的光场PTV气泡跟踪测速方法

王粤，孙凯，刘艳，陈龙，朱效宇，许传龙

（1 东南大学能源与环境学院，江苏 南京 210096；2 江苏中科能源动力研究中心，江苏 连云港 222000；3 中国科学院工程热物理研究所，江苏 连云港 222000）

气液两相流广泛存在于能源动力、航空航天、冶金化工等领域，气泡作为气液两相流中的分散相，其参数是表征气液两相流动的重要参数之一[1-2]。因此，准确可靠的气泡特性参数测量对于气液两相流动机理研究和相关系统的设计具有重要意义[3-4]。由于实际过程中气液两相流动多发生在三维空间内，为了获得完整的气液两相流动结构，迫切需要发展先进的气泡特性参数三维测量技术[5]。

为实现气泡形貌以及三维位置测量，传统方法中多采用双目或多目相机的布置方式，从多个视角拍摄气泡投影图像，重建气泡的质心位置和半径[6-7]。该方法测量结果精度较高、视场较大，但由于测量系统包含多台相机，系统的结构较为复杂，标定难度高，无法应用于光学视窗受限的测试场景。近年来，随着计算成像技术的发展，基于光场成像的粒子跟踪测速（PTV）技术在单视角条件下实现了气泡位置的三维重建以及气泡运动轨迹的跟踪，为受限空间条件下气泡三维测量问题提供了一种新的解决方案[8-9]。光场成像气泡特性参数PTV测量技术主要包括光场图像的采集、气泡空间位置的重建以及气泡运动PTV 跟踪三个主要技术环节。其中，气泡运动PTV 跟踪是在气泡位置识别的基础上，对连续两帧或多帧图像中的气泡进行匹配，从而获得气泡运动轨迹和速度。因此，为获得准确的气泡运动速度场，需要发展先进的PTV 匹配算法，以确保多帧图像间气泡的正确匹配。

现有的PTV 气泡匹配方法中，根据匹配过程中搜索半径的大小，可以分为局部和全局优化搜索两类。其中，局部优化搜索算法主要包含最近邻搜索、匹配概率算法以及极坐标系统相似（PCSS）算法等。最近邻搜索是应用最广泛的匹配算法，其原理是以第一帧图像中目标气泡的位置为起点，在第二帧图像中寻找与该点距离最近的气泡，认为此气泡与第一帧图像中的目标气泡相匹配。该算法计算过程简单，但仅适用于气泡浓度低、位移小的场景中[10]。匹配概率算法是在最近邻搜索的基础上，考虑了第一帧图像中的气泡离开光平面无法在第二帧图像中成像的情况，因此匹配准确率较最近邻法有所提高，但该方法对不同含气率的气液两相流场进行气泡匹配时需要人为设置先验参数，较为繁琐[11]。PCSS 算法的原理是先在第一和第二帧图像中分别以目标气泡和候选气泡为圆心作查询区域，选取查询区域内相似系数最大的候选气泡作为匹配气泡。该方法原理简单，匹配准确性对算法参数选取不敏感，但该方法对于气泡在测量平面外的位移难以准确判别，因此同样仅适用于两帧间气泡位移量小的场景中[12]。全局优化搜索算法中，主要包括线性松弛方法和特征向量方法等。松弛匹配方法在高度随机的气泡运动场景表现良好，但是需要保证较小的帧间气泡运动位移。基于特征向量的气泡匹配算法能够解决气泡位移较大问题，但是该方法需要较高的气泡浓度来构建运动特征。

总的来说，现有的局部优化搜索方法的计算速度较快，但由于仅将被测粒子与附近多个气泡进行特征匹配以获取最优匹配气泡，不适用于高气泡浓度和大位移场景，且无法筛除由于气泡重叠效应所产生的幽灵气泡。全局优化搜索方法将整个颗粒流场定义成一个整体，建立与全部粒子匹配相关的目标函数，然后通过多种优化方法获得最终的匹配收敛解，虽然精度较高，但由于需要对整个流场内的气泡进行匹配，搜索维度较高，优化计算过程容易陷入局部极值。因此，有必要将两类方法进行融合，以提高气泡运动参数测量的准确性。为此，本文提出基于全局位置迭代的三维PCSS（GPPIPCSS）方法。在采用PCSS 算法进行气泡三维匹配的基础上，对图像内所有气泡的位置进行迭代更新，使两帧图像中的气泡位置逐渐重合，从而提高两帧图像中气泡匹配以及气泡速度场测量的准确性。通过在鼓泡床装置上开展气泡运动行为光场PTV测量实验，对所提出的GPPI-PCSS方法的准确性进行了评价，并分析了采样间隔和空气流量工况对气泡速度场测量准确性的影响。

1 光场PTV气泡跟踪测速原理

1.1 光场成像及计算重聚焦

光场成像作为一种计算成像技术，通过采集包含光线位置与方向的四维光场信息并对光场信息进行后处理，以获得目标场景的三维空间位置信息[13-14]。光场信息的采集可以由光场相机来实现，其主要光学元件包括主透镜、微透镜阵列（MLA）以及图像传感器CCD。光场相机工作原理示意图如图1所示，由物空间发出的不同方向的光线穿过气泡后，首先经主透镜汇聚于微透镜阵列上。此后，微透镜根据光线的入射方向将其投射到图像探测器的不同像素上。在上述过程中，光线的位置和方向信息分别被对应位置的微透镜和像素所记录。根据所记录的光场信息，利用计算重聚焦技术可以在一定深度范围内自由改变聚焦平面的位置，由此可以获得气泡在不同深度平面上所成的像。

图1 光场成像原理示意图

光场计算重聚焦技术的原理如图2所示，假设相机的聚焦平面上存在两点O1和O2，根据微透镜阵列与聚焦平面的共轭关系，从O1和O2出的光线应该分别汇聚于微透镜阵列上的i1和i2两点。将光线从i1和i2向后继续延长，两条光线汇聚于虚拟点i1'，点i1'所在平面被称为重聚焦像面，而在物空间内与i1'共轭的物点O1'所在的平面则被称重聚焦物平面。利用光场图像中的光线方向信息将不同方向的光线重新投影到新的像面上，即可获得不同空间深度位置对应的重聚焦图像。将所有重聚焦图像上最清晰的部分组合就可以得到全聚焦图像。全聚焦图像中每个气泡投影都清晰成像，通过对全聚焦图像进行特征分割，可以得到图像内所有气泡的轮廓和质心位置。

图2 光场数字重聚焦原理示意图

1.2 气泡空间位置三维重建及PTV测量

气泡三维重建是指根据气泡质心坐标mc和气泡投影面积当量直径de，从图像中还原气泡的真实空间位置的过程[15]。气泡三维重建模型如图3所示，图中D表示物空间内一气泡，其质心坐标为mc(xc,yc,zc)，截面积为A，当量直径为de。气泡在像空间内的聚焦图像为D'，D'的质心坐标为m'c(x'c,y'c,z'c)，投影面积为A'，当量直径为d'e。根据成像放大倍率关系，物空间内气泡的质心位置、当量直径、截面积可由式(1)～式(6)计算获得。

图3 气泡三维重建模型

式中，zc为气泡质心深度，即为气泡投影还原到物空间后各处深度d(x,y)的平均值。放大率M是物空间深度的函数，即M=g(zc)，在像空间则有式(7)～式(10)。

气泡三维重建流程图如图4所示。其中，气泡三维重建所需的图像数据来源于光场相机获取的全聚焦图像和重聚焦图像序列。重聚焦图像序列的深度和放大率数据需要通过毫米尺标定法进行系统标定获取[15]。对全聚焦图像进行处理，得到气泡图像二维质心坐标及尺寸；对重聚焦图像序列进行处理，结合系统标定，确定气泡质心深度和该深度对应的放大率，计算物空间气泡的三维质心坐标以及尺寸参数，进而完成气泡的三维重建。

图4 气泡三维重建流程图

2 基于全局位置迭代的PCSS 气泡跟踪测速算法

2.1 算法原理

PTV技术是一种图像测速技术，其主要原理是识别第一帧图像（t时刻）上的每个气泡在第二帧图像（t+Δt）上的位置[16]，通过对两帧图像中的相同气泡进行匹配，以获得气泡在两帧图像拍摄时间间隔Δt内的运动位移Δl并计算气泡的运动速度，如式(11)。

因此，在采用PTV 算法追踪气泡的运动轨迹的过程中，其关键问题在于如何对相邻两帧光场图像中的气泡进行三维位置匹配。针对该问题，简要介绍二维和三维PCSS 算法原理，并在此基础上提出了基于全局位置迭代的PCSS（GPPI-PCSS）气泡跟踪测速算法。

二维PCSS 算法原理示意图如图5 所示。在第一帧图像A 中，气泡i为目标气泡，以其为圆心计算半径R1的邻域气泡，并分别计算气泡i与邻域气泡im的连线长度Rim以及两气泡中心连线与X轴的夹角θim。在第二帧图像B中，气泡i'为帧A中的气泡i在帧B中的匹配气泡，以i'为圆心、R2为半径，筛选虚拟气泡i'在帧B 中的候选气泡。假设气泡j为虚拟气泡i'的候选气泡，再以气泡j为中心、R1为半径，划定气泡j的邻域气泡，分别计算气泡j与邻域气泡jn的连线长度Rin以及两气泡中心连线与X轴的夹角θin。用Sij表示气泡i和气泡j之间的相似系数，Sij计算见式(12)。

图5 二维PCSS算法原理示意图

式中，M表示第一帧中气泡i的邻域气泡个数；N表示第二帧中气泡j的邻域气泡个数；εr和εθ分别表示长度偏差与角度偏差的阈值；H(vr,vθ)是一个阶跃函数，当vr和vθ均大于0时结果为1。

对气泡i在第二帧图像中所有候选气泡j，都可以根据式(12)计算出相似系数Sij，选取其中相似性系数最大的候选气泡j作为气泡i的匹配气泡。对第一帧图像中所有气泡执行上诉匹配过程，即可完成整幅图像的气泡匹配[17-18]。

三维PCSS算法在二维PCSS算法的基础上分别增加气泡连线与Y、Z轴夹角的角度偏差阈值，并计算了气泡i与邻域气泡jm的连线长度rim以及两气泡 中 心 连 线 与X、Y和Z轴 的 夹 角θix、θiy和θiz。三维PCSS 算法原理示意如图6 所示，其相似性系数Sij可以由式(13)计算。

图6 三维PCSS算法原理示意图

三维PCSS 算法原理简单，设置参数少，计算速度快，且不同的实验条件下对参数选取不敏感。然而，三维PCSS方法本质上是一种局部搜索算法，无法解决匹配过程中的幽灵气泡现象，也无法适用于高气泡浓度工况。因此，对三维PCSS 方法进行改进，提出基于全局位置迭代的三维PCSS 测量方法（GPPI-PCSS）。GPPI-PCSS 算法的计算流程如图7所示。首先，获取帧A和帧B图像中气泡的三维位置，并通过PCSS 方法计算帧A 中气泡在帧B中Nump个邻域气泡的相似性系数，继而获得帧A中气泡在帧B中的匹配结果。然后，计算匹配气泡之间的三维速度矩阵，并将速度矩阵插值为三维速度场，计算帧B中所有气泡在三维速度场中的速度变化值。进一步，将帧B中所有气泡的位置进行更新，同时减小Nump的值，缩小检测范围，继续该循环，直至速度矩阵的范数接近于0或者连续3次匹配率大于0.99，则退出循环，当前两帧图像的PTV测量过程结束。最后，对连续多帧图像的检测结果按顺序进行连接，完成多帧图像的PTV测量。

在上述采用GPPI-PCSS 算法进行气泡位置迭代优化的过程中，对邻域气泡数量Nump的更新如式(14)所示。

式中，Ratiom为前一迭代过程的气泡匹配成功数目占总气泡数目的比例。

速度矩阵的范数Normuvw的计算为式(15)。

式中，tempu、tempv和tempw分别为u、v、w方向的速度矩阵；kuvw为速度矩阵的长度。

2.2 幽灵气泡排除

幽灵气泡指在高截面含气率工况下，由于遮挡、分割失败等原因只存在于一帧图像内的气泡。幽灵气泡会导致相邻两帧图像的气泡出现误匹配现象，严重影响PTV 的测量精度[19-20]。因此，需要对幽灵气泡进行检测排除。通过GPPI-PCSS 气泡匹配迭代，当邻域气泡数量Nump降低到1 时，帧A和帧B图像中的气泡已经高度重合，剩余的未匹配气泡即被认为是幽灵气泡。根据这一特征，当邻域气泡数量Nump降低至1时，对幽灵气泡进行排除。该方法被耦合在GPPI-PCSS 优化过程中，幽灵气泡检测过程共分为两步，其计算见式(16)。

式中，PA和PB为更新的气泡场；PA和PB为原始气泡集合；εd为匹配气泡的临界距离；û为对应的匹配气泡序号。

首先，对帧A气泡场中的所有气泡在帧B中进行最近邻搜索，计算帧A气泡场中气泡与匹配气泡的距离，若匹配气泡间的距离大于阈值，则认为该气泡是帧A图像中的幽灵气泡，将该气泡剔除。同样的，将帧B图像中所有气泡在经过筛选的帧A气泡场中进行最近邻搜索，计算帧B气泡与匹配气泡距离，若出现气泡的距离大于阈值，则认为该气泡是帧B中的幽灵气泡，将该气泡剔除，最终得到两帧图像对应的匹配气泡。剔除幽灵气泡后的气泡融合图像如图8所示。可以看出，幽灵气泡分为两类：一类是由于气泡运动出相机视野或未进入相机视野，导致仅有一帧图像包含该气泡，这种情况出现于相机图像的上下部分；另一类是由于气泡运动过程中产生的重叠或重叠气泡构成发生变化，导致气泡分割失败，产生了幽灵气泡。当邻域气泡数量Nump降低至1时，经过幽灵气泡的排除过程，帧A中的3个幽灵气泡和帧B中的4个幽灵气泡被成功检测出来。值得注意的是，对于帧B和帧C中幽灵气泡的排除，可以重复以上步骤，其结果并不受帧A的影响。

图8 排除幽灵气泡后的气泡融合图像

3 气液两相流气泡光场PTV 测量性能评价

3.1 实验系统

为了进一步验证基于全局三维极坐标相似算法的光场PTV 气泡跟踪测速方法的可行性，在鼓泡床实验装置上开展了气泡特性参数测量实验研究。鼓泡床实验装置如图9所示，主要包括Raytrix R29相机、100mm 定焦镜头、曝气器、水槽、LED 光源、同步控制器、位移台和数据处理计算机。曝气器由10 个并排布置的0.2mm 孔径的喷嘴组成，喷嘴阵列垂直于相机光轴方向，喷嘴中心截面的相对深度为0。Raytrix 相机连续拍摄模式下最大帧率仅为6 帧/s，因此通过相机的双曝光PTV 模式，对气泡流场进行连续拍摄。实验在暗室条件下进行，在拍摄第一帧图像时，相机和光源同时触发，光源点亮时间和相机曝光时间一致，通过调节第二次光源的触发时间，可以实现不同时间间隔下的气泡流场图像的拍摄。

图9 气泡PTV测量实验装置

3.2 结果与讨论

3.2.1 采样间隔对光场PTV测量性能的影响

为了评估气泡位移对于PTV 匹配结果的影响，在0.25L/min 的空气流量下，分别对采样间隔为0.5ms、1ms、3ms、5ms 和7ms 工况下的气泡图像序列进行PTV匹配实验，每种工况均采集了200张气泡图像。图10 为不同采样间隔下的气泡PTV 匹配结果图像。可以看出，随着采样间隔的增加，二维气泡图像上的对应气泡位移显著增加，对于0.5～3ms工况下的气泡图像，GPPI-PCSS方法的检测精度较高，大部分气泡都被成功匹配。然而，当采样间隔大于3ms 时，气泡错误匹配的概率显著增加，这是由于采样间隔较长时，气泡的高速无序运动会导致部分气泡产生较大的位移或与其他气泡发生重叠，严重影响了气泡匹配的精度。因此，气泡光场PTV实验的最大采样间隔应低于3ms。

图10 不同采样间隔下的气泡PTV匹配结果图像

图11 为不同采样间隔下的气泡速度场测量结果。可以看出，不同工况下，气泡的运动速度范围为15.85～47.89cm/s，而由于气泡的空间分布分散，且部分气泡的运动速度值较小，经过插值显示后，这些气泡周边区域的速度值明显偏小。随着采样间隔由0.5ms 增加至7ms，PTV 测量得到的气泡运动速度的整体呈现出逐渐减小的趋势。这是由于三维气泡运动的轨迹多变，当采样间隔较长时，气泡在采样间隔内可能产生复杂的三维运动，导致气泡的实际运动距离较短，气泡运动速度值小于实际值。因此，为了提高PTV 测量的精度，系统应设置尽可能小的采样间隔。

图11 不同采样间隔下的气泡速度场测量结果

为了评价GPPI-PCSS 方法的性能，分别利用GPPI-PCSS 方法、三维PCSS 方法和松弛方法对不同采样间隔工况下的气泡图像进行PTV 匹配，得到的PTV 测量结果曲线如图12 所示。从图12(a)可以看出，在0.5ms 采样间隔工况下，GPPI-PCSS、PCSS和松弛方法的匹配准确率高达95.52%、89.65%和93.71%，而随着采样间隔的增加，三种方法的匹配准确率均相应降低。这表明在较大的气泡位移工况下，PTV方法会出现较多的误匹配现象，匹配精度显著降低，因此PTV 测量的采样间隔应该尽可能短。同时，在不同采样间隔工况下，GPPIPCSS 方法的匹配准确率远高于其他两种方法。对正确匹配的气泡序列计算平均运动速度，结果如图12(b)所示。可以看出，随着采样间隔由0.5ms增加到9ms，气泡的平均运动速度由40.23cm/s 降低至29.45cm/s。综上所述，相较于PCSS 方法和松弛方法，GPPI-PCSS方法具有更高的匹配准确率，但为了提高PTV 测量的精度，应设置尽可能小的采样间隔。

图12 不同采样间隔下的气泡PTV匹配结果曲线

3.2.2 截面含气率对光场PTV测量性能的影响

气液两相流中，随着注入气体流量的增加，会产生更多气泡，流场内的气泡密度和截面含气率也相应变大，这会影响PTV 方法的检测性能。当空气流量大幅增加时，喷嘴产生的气泡会出现剧烈变形，导致单个气泡被误分割或者重叠气泡无法被准备分割。因此，为了评估截面含气率对PTV 匹配结果的影响，设置采样间隔为1ms，对于0.2mm孔径喷嘴，分别在0.15L/min、0.2L/min、0.25L/min，0.3L/min 和0.35L/min 的空气流量下进行了PTV 匹配实验，每种工况均采集了200组气泡图像。

图13 为不同空气流量下的气泡PTV 匹配结果图像。可以看出，当空气流量小于0.25L/min 时，气泡流场截面含气率较低，图像内重叠气泡的数量较少，绝大部分气泡对都被成功检出。而当空气流量增加到0.3L/min 和0.35L/min 时，气泡密度显著增加，有较多气泡出现重叠现象，虽然对重叠气泡进行了分割，但是仍然存在较多的幽灵气泡，影响了GPPI-PCSS 方法的检测精度。同时发现，由于喷嘴无法保证垂直向上，随着空气流量的增加，气泡的初始速度增加，在气泡表面张力、黏性力和惯性力的作用下，气泡在深度方向上的分布也越发扩散，这表明通过三维气泡光场PTV 测量方法，可以获得准确的气泡运动特性参数。

图13 不同空气流量下的气泡PTV测量结果图像

图14 为不同空气流量下的气泡速度场测量结果。可以看出，随着空气流量由0.15L/min 增加到0.35L/min，气泡的最大运动速度由38.34cm/s 增加至49.87cm/s，气泡的整体运动速度大小也有较大提升。这是由于随着空气流量的增加，喷嘴的尺寸保持不变，气泡的初始速度变大，气泡流场的速度值也相应增加。同时发现，气泡流场的速度分布并不均匀，部分区域会出现速度较大或较小的情况，这主要是由气泡运动特性所导致的。相对于流场测量区域，气泡数量相对偏少，由于喷嘴歪斜导致的气泡分布不均以及气泡在脱离喷嘴后的复杂运动，不同气泡的运动速度差异较大。而气泡速度场是由所有气泡的运动速度插值而成，气泡速度差异导致了气泡速度场截面上的速度分布相差较大。因此，在设计气泡发生系统时，需要同时考虑喷嘴孔径、喷嘴数量以及喷嘴布置形式。

图14 不同空气流量下的气泡速度场测量结果

图15 为不同空气流量下的截面含气率和气泡PTV 测量结果曲线。可以看出，当空气流量由0.15L/min 增加到0.35L/min 时，气泡流场的截面含气率由8.35%增加到17.38%，GPPI-PCSS、PCSS和松弛方法的匹配准确率均呈下降趋势，但GPPIPCSS 方法的匹配准确率均高于PCSS 和松弛方法。因此，在低截面含气率工况下，GPPI-PCSS、PCSS和松弛方法都能保持较高的测量精度，但是在高截面含气率工况下，PCSS 方法和松弛方法的匹配准确率远低于GPPI-PCSS 方法。GPPI-PCSS 方法通过对气泡空间位置进行全局迭代更新，使两帧图像内的气泡接近重合，从而大大提升了PTV 测量的准确性，同时也解决了传统PCSS 方法无法处理幽灵气泡和高气泡密度的情况。

图15 不同空气流量下的截面含气率和气泡PTV测量结果曲线

对不同工况下正确匹配的气泡序列计算平均运动速度，结果如图15(b)所示。可以看出，随着空气流量由0.15L/min 增加到0.35L/min，气泡平均运动速度由34.15cm/s 增大至46.45cm/s。综上所述，虽然GPPI-PCSS 方法在高含气率下仍能够保持较高的测量精度，但是由于重叠气泡的复杂构成，不可避免地会出现分割错误的情况，所产生的幽灵气泡会极大地影响PTV 测量的准确率。因此，在设计曝气器时，要设置合理的气体喷射流量，减少复杂重叠气泡的生成。

4 结论

本文提出了一种基于全局位置迭代极坐标系统相似技术（GPPI-PCSS）的三维气泡匹配算法，通过匹配气泡的三维运动速度场，对单帧图像内所有气泡的位置进行迭代更新，使两帧图像中的气泡位置逐渐重合，从而增加两帧图像中匹配气泡的数量，提高传统PCSS 算法气泡匹配的准确性。构建了光场PTV 系统，在气液鼓泡床装置上开展了气泡流场光场PTV 实验测量研究，对所提出方法的准确性进行了评价，结果如下。

（1）在不同采样间隔和截面含气率工况下，GPPI-PCSS 方法的平均匹配准确率为92.65%，高于PCSS方法的81.38%和松弛方法的84.12%，且气泡运动速度的变化趋势与理论相一致。证明了基于GPPI-PCSS方法的光场PTV技术可以实现气泡运动参数的准确测量。

（2）对于采样间隔在0.5～3ms工况下的气泡图像，GPPI-PCSS方法的检测精度较高，大部分气泡都成功匹配。当图像采样间隔过大时，由于气泡的无序运动，气泡匹配准确率和气泡运动速度都相应降低。随着空气流量的增加，重叠气泡的比例增多，会出现较多由于分割失败而产生的幽灵气泡。因此，通过设置合理的喷嘴孔径、喷嘴数量以及喷嘴布置形式，可以降低重叠气泡出现的概率，有效提升气泡匹配的精度。

（3）当空气流量小于0.25L/min 时，气泡流场截面含气率较低，图像内重叠气泡的数量较少，绝大部分气泡对都被成功匹配。当截面含气率过高时，由于重叠气泡的复杂构成，不可避免地会出现分割错误的情况，所产生的幽灵气泡会极大地影响PTV测量的准确率。因此，在设计曝气器时，要设置合理的气体喷射流量，减少复杂重叠气泡的生成。