层析粒子图像测速技术研究进展

李晓辉, 王宏伟, 黄 湛, 赵俊波

中国航天空气动力技术研究院, 北京 100074

0 引 言

流场速度测量不仅直接关系到诸如飞机、汽车、高铁等流线型流体动力机械的外形设计和优化，还制约着流体力学的发展。对于流场速度的测量，国内外很早就开始了研究，出现了多种技术方法：如激光多普勒测速(Laser Doppler Velocimetry, LDV)、热线风速仪等单点测量技术；平面激光诱导荧光技术(Planar Laser Induced Fluorescence)、激光散斑测速技术(Laser Speckle Velocimetry, LSV)、粒子图像测速技术(Particle Image Velocimetry，PIV)等非接触二维面测量技术。尤其是PIV技术，不仅能定量地对流场进行时间分辨率上的流场演化测量，还能定性地显示流场结构，对理解流动机理及促进型号研制起到了巨大的推动作用。

但无论是目前倍受关注的湍流、涡流等流体力学基本问题，还是飞行器、舰船、导弹等工程型号研制[1]，都面临着复杂三维非定常流动问题。尤其以隐身、高机动性为特征的新一代飞行器面临的复杂三维流场[2-5]，以微型飞行器和机器鱼为背景的生物和仿生流体力学[6-7]，以跨、超声速混合、燃烧为背景的非定常激波膨胀波引起的可压缩剪切复杂流动[8-10]，以微机电系统为背景的微尺度复杂流动等[11-12]，无不具有强非定常性、强三维空间性等复杂流动现象，而这些问题是二维测量技术无法解决的，因此发展新的三维测量技术势在必行。

目前关于攻克3D3C(Three-dimensional three component)测量技术途径的主流仍然是粒子图像测速技术的延伸和发展，主要包括散焦PIV[13-14]、扫描PIV[15]、全息PIV(Holographic PIV,HPIV)[16]及层析PIV(Tomographic PIV, Tomo-PIV)[17-18]等。其中散焦PIV实质上是一种粒子跟踪测速技术，空间分辨率较低；扫描PIV是基于标准的二分量PIV，在相继的深度位置上进行扫描，再将空间上有微小偏移的平面速度场结合起来分析，这种方法具有较高的空间分辨率，但是扫描时间与被观测流场的特征尺度相比小得多，因此并不适用于高速流场；HPIV是一种将数字全息技术和PIV技术相结合产生的三维流场速度测量技术，通过CCD记录空间场的全息干涉条纹，再利用相应的反演公式进行空间流场信息的再现，其光学条件十分苛刻，距离实际应用还需开展更多的研究；Tomo-PIV技术结合了医学CT技术和PIV技术，能够实现空间流场的全场定量测量，对于复杂流场结构的研究能够提供很大的便利，引起了众多研究者的关注。

本文在全面调研的基础上对Tomo-PIV技术进行了综述，介绍了相应的工作原理和技术特点，重点讨论了对重构精度有较大影响的技术因素，总结了国内外关于该技术的典型应用，展示了Tomo-PIV在非定常三维复杂流场测试方面的优势。

1 技术简介

1.1 基本原理

Tomo-PIV是将PIV技术和医学上的CT重构技术结合而产生的一种三维粒子图像测速技术，能够实现空间流场的全场定量测量，其具体原理如图1所示。在测量区域内释放一定浓度的示踪粒子，利用高能量的脉冲激光经体光源照射测量区域内的示踪粒子，同时使用多个(一般为3～6个)相机记录粒子运动图像，然后根据照片中像素灰度的分布，利用倍增代数重构算法(Multiplicative Algebraic Reconstruction Technique, MART)对三维空间粒子分布进行迭代重构，对重构的2个相邻粒子图像进行三维互相关计算，得到粒子的位移信息，再根据曝光时间计算出速度场。

图1 Tomo-PIV工作原理Fig. 1 Principles of Tomo-PIV

1.2 技术特点

Elsinga等[17-18]在加利福尼亚举办的第六届国际PIV专题研讨会上第一次提出Tomo-PIV技术，阐述了其工作原理和相应算法，其中最为关键的为三维空间标定及三维粒子重构。

Tomo-PIV是对三维体空间进行重构，因此试验前需要准确知道相机图像坐标和空间物理坐标(重构体)之间的对应关系，以便通过拍摄的二维图像重构出粒子的三维空间分布，这一过程通过标定来完成。标定函数的精度直接决定粒子空间位置的重构误差，进而影响速度场的计算。标定时标定板沿测量体厚度方向遍历整个测量体，每个相机记录不同景深位置Z处的标定靶图像，通过针孔照相机模型[19]或者多项式模型[20]获取标定映射函数。

在Tomo-PIV流场测试中，相机接收的不是流场本身的信息，而是示踪粒子的散射光强。将测量区域进行空间网格离散，得到离散的三维体素；在激光的照射下，假设体素内示踪粒子散射光强度为E(X,Y,Z)，投射到相机平面形成像素灰度I(x,y)，两者的关系可表示为：

(1)

其中，(X,Y,Z)为体素三维空间坐标，(x,y)为图像像素坐标，Ni是对平面像素(xi,yi)灰度有影响的体素数量，i代表相机的第i个像素，ωij是加权系数，表示第j个体素对第i个像素的强度贡献率。

Elsinga等[17-18]在提出Tomo-PIV技术时使用MART算法得到三维空间的粒子分布，该方法通过给定所有体素统一的初始光强E0(X,Y,Z)，迭代出粒子三维分布：

(2)

2 Tomo-PIV技术研究现状

针对Tomo-PIV技术本身的研究主要集中在如何提高其重构精度。在所有影响粒子重构精度的因素中，最主要的是相机布局、示踪粒子密度、标定映射函数及三维重构算法等，众多学者对其开展了研究。

2.1 相机布局

和平面PIV不同，Tomo-PIV采用多相机拍摄，相机的不同布局、相机之间及相机与测量体之间的角度均会影响拍摄的图像质量。Elsinga等[17]通过数值模拟研究了相机角度对重构精度的影响，认为相邻相机之间的最优角度在30°左右。相机之间角度变小，示踪粒子在深度方向被拉长，使得示踪粒子在深度方向的分辨率和其他2个方向差别较大，粒子的灰度峰值位置难以精确确定，对后续数据处理造成不利影响；相机角度变大，沿着相机视线在测量体内的光程变大，在重构时会形成大量的虚假粒子。

常用的相机布置方式有“十字”交叉式和线性布置2种，如图2所示。Scarano等[21]研究了不同相机孔径角β(最远2个相机之间的夹角)对2种布局方式重构质量的影响，发现“十字”交叉式布局的重构精度比线性布局要高，且最远2个相机之间的夹角最好在40°～120°之间(见图3)。

图2 相机布局方式[21]Fig. 2 The camera layout[21]

图3 相机系统孔径角对重构质量的影响[21]Fig. 3 Reconstruction quality factor versus system aperture angle [21]

2.2 示踪粒子浓度

图4为Michaelis等[22]测试的示踪粒子在不同浓度下的拍摄图像，可见示踪粒子的密度对测量精度影响较大。为了从平面图像中重构出高精度的粒子三维空间信息，需要保证二维粒子图像相互不重叠，这通常通过降低示踪粒子浓度或减小测量体空间厚度的方法实现；但是为了提高测量的空间分辨率，示踪粒子浓度不能太低。目前为了平衡示踪粒子浓度变大造成的重构质量降低和示踪粒子浓度变小造成的空间分辨率降低两者之间的矛盾，通常推荐的示踪粒子浓度为0.05 ppp(particles per pixel, ppp)[17]。

图4 不同示踪粒子浓度的实验图像(左为示踪粒子浓度)[22]Fig. 4 Images of experiments at different particle concentrations (the left is the concentration of particles, ppp)[22]

2.3 标定映射函数

Tomo-PIV采用分段式标定，分别记录若干景深位置Z处的图像，使用三阶多项式拟合物理空间坐标和图像像素坐标之间的映射关系，其具体形式为：

F(x)=a0+a1X+a2Y+a3X2+a4XY+

a5Y2+a6X3+a7X2Y+a8XY2+a9Y3

F(y)=b0+b1X+b2Y+b3X2+b4XY+

b5Y2+b6X3+b7X2Y+b8XY2+b9Y3

(3)

其中，(x,y)为二维像素坐标，(X,Y)为三维空间物理坐标，系数(ai,bi)是不同标定位置Z的函数。只要有足够多的标定点，通过最小二乘求解即可得到不同标定位置Z上的标定系数。未标定的Z方向位置通过对Z方向上相邻位置上映射的系数线性插值得到。因此，标定得到的映射函数在X、Y方向上具有三阶精度，在Z方向上具有一阶精度。

Tomo-PIV的标定映射函数误差在0.5个像素以下时才能保证三维重构的精度[17]。然而，由于不精确的标定板、不精确的标定板移动、机械的不稳定性、光学畸变以及其他可能的相关因素，标定误差很难控制。当标定误差很大时，不同相机中的同一粒子在沿相机视角计算时其空间位置往往无法重合，需要进行映射函数的体自标定[23]。

图5 标定残差示意图Fig. 5 Schematic diagram of calibration residuals

可以看出，体自标定是通过真实粒子来修正映射函数，因此需要精确匹配粒子的三维位置，这通常通过三角测量法来完成。如图6所示，对于相机1中的每一个粒子，通过标定函数计算其三维空间位置，则其对应在相机2中的粒子位置在宽为2εr、长为Lz的长方形条带内(其中，εr为给定的误差半径，比预计的最大标定误差大，Lz为相机1视角方向在测量体内的长度在相机2视角方向上的投影长度)；结合相机1和相机2的图像确定粒子的三维位置，然后通过相机3及相机4验证相机2条带中的每一个粒子，完成粒子的匹配。从上述方法可以看出，用于体自标定的示踪粒子浓度要较为稀疏，否则不仅计算时间长，还很难真正确定真实粒子的位置，因此通常要在试验之前记录一组稀疏粒子图像。

图6 粒子匹配示意图[24]Fig. 6 Schematic diagram of particle matching[24]

此外，Schanz等[25]利用光学传递函数(Optical Transfer Functions, OTF)来补偿非一致性拍摄条件(如前向散射和后向散射、窗口导致的光学畸变等)形成的图像误差，提高了三维重构的精度。

2.4 重构技术

三维粒子重构是Tomo-PIV的核心，针对粒子重构的算法优化，众多研究者开展了深入研究。本文对近年来典型的三维重构算法进行了汇总，详见表1。

表1 典型三维重构算法Table 1 Typical three-dimensional reconstruction algorithm

重构算法的关键是尽可能地提高重构精度，同时减小重构的时间成本。根据主要目的不同，重构算法可分为提高重构精度、减少重构时间以及两方面同时作用等3类。

1) 提高重构精度。此类方法主要从减少体素的离散误差、减少几何视角误差以及虚假粒子抑制等3个方面进行。Tomo-PIV采用多相机系统进行空间流场拍摄，相机与测量体之间具有一定的角度，导致测量体厚度方向的空间分辨率较低，因此在重构过程中需要根据实际空间分辨率大小来设置体素的形状和大小，减小网格离散误差，提高计算效率和精度。几何视角误差主要体现在相机布局方面，在2.1节中已详细阐述。如图7所示，采用MART算法进行重构时会出现较多的虚假粒子，即在本不该出现粒子的视线交点处出现类似粒子的灰度分布[37]。对于N相机的Tomo-PIV系统，可用真实粒子和虚假粒子的比值来表征信噪比[38]：

(4)

图7 虚假粒子Fig. 7 Ghost particles

其中，Np为真实粒子的数量，Ng为虚假粒子的数量，ppp为单位像素的粒子数量，Ap为单个示踪粒子面积，lz为重构体厚度方向的尺寸大小。

Ap可认为是不变的，则影响的因素主要是示踪粒子浓度和测量体厚度的大小。研究表明，在ppp小于0.05的情况下重构具有较好的效果，随着粒子浓度的提高，虚假粒子所引起的误差变得越来越显著，甚至导致测量速度不可信，这也是Tomo-PIV测量空间通常为扁平长方体的原因。

2) 减少重构时间。MART算法采用相同的初始光强值进行迭代，不仅影响重构的精度,还需要较长的计算时间和存储空间，目前研究的主要方向是寻找更加合适的迭代初值，从而加快迭代的收敛速度。

3) 兼顾三维重构的精度和时间。目前针对此类问题还未见有文献开展专门的研究，有部分算法实现了两方面的同时提高，但是偏重点不同。

2.5 速度场后处理

Tomo-PIV技术的不断完善成熟，让人们对三维复杂流动的试验研究寄予更大的期许，但是和二维PIV相比，其图像质量及测量空间分辨率均有所不足，因此数据后处理技术是Tomo-PIV不可或缺的一部分，其主要包括速度场及其导出量的处理。

常见的速度场后处理有坏点剔除、中值滤波和高斯平滑等。坏点剔除是通过给定速度分布区间，删除流场中明显偏离平均速度的矢量，对全流场或局部流场进行处理，剔除后的流场可以通过33空间插值来进行修补。中值滤波是基于排序统计理论的一种能有效抑制背景噪声的非线性滤波技术，和高斯平滑类似，它对流场有较强的平滑作用。

此外，由于Tomo-PIV采用体互相关算法计算获得速度场，判读体大小即为其空间分辨率，因此算法对速度场有平滑机制，会在一定程度上低估速度梯度。对于不可压缩速度场，如果不能精确获得速度梯度，会造成散度不为0，而速度梯度的偏差将直接影响旋涡的识别和旋涡强度计算[39]。高琪等[40]提出了一种基于不可压缩连续性方程制约机制的流场修正方法，通过求解全流场修正速度场二范数的极小值来获得最优的修正量，使修正后流场完全满足差分形式的连续性方程，对流场中的高斯分布噪声有一定的消除作用，能够改善测量速度场品质。

速度场导出变量主要为涡量的计算和识别。常见的涡识别方法主要有Q准则法[41](速度梯度张量第二不变量)、Δ方法[42](速度梯度张量的特征方程的根的判别式)、λci方法[43](速度梯度张量的共轭复特征值的虚部)、λ2方法[44](压力的海森矩阵的第二特征值)。在这4种方法中，λci方法严格满足流线的基本形状或流体微团轨迹呈螺旋状这一基本性质，直接以描述螺旋运动强度的参数作为涡识别变量，在数学基础和物理意义上更为明确，比其他3种方法应用更为广泛。

复杂流场往往包含不同尺度、不同形态的流动结构，且流动结构间存在非线性的耦合干扰。从复杂流场中辨识出主要流动结构是认识流动特征规律的前提。对于时间解析的层析PIV来说，本征正交分解(Proper Orthogonal Decomposition, POD)、动力学模态分解(Dynamic Mode Decomposition, DMD)等方法能够实现三维速度场的模态分解，为研究复杂流场的动力学特征、获取低维动力学经验模型提供模态分解工具；此外，变分模态分解(Variational Mode Decomposition, VMD)作为一种自适应模态变分处理方法，能够对瞬时速度场进行有效的模态分析，处理流场的非稳态行为，十分适合用来开展复杂流场的分析。

2.6 三维PIV的创新发展

Tomo-PIV是当前最为成熟、应用最为广泛的三维流场测试技术之一，具有空间分辨率高及测量体积大的优点；但其本质属于一种多相机测试系统，平台搭建及系统操作难度较大。为降低系统的复杂度，国内外学者开展了进一步的研究。

高琪等[45]发展了一种单相机三维体视PIV技术，在相机与被测流场之间加装一个三棱镜特效透镜，光线通过该透镜3个棱面的折射能实现多相机不同视角成像的效果，经过三维粒子重构，进而实现了三维流场的测量，并利用该技术获取了零质量射流涡环三维流动结构的时序结果。但由于3个不同视角的成像通过一个相机来实现，因此存在有效测量区域受限的问题。

随着光场三维成像技术的成熟及光场相机的研制，使得单台相机同时记录光线的强度及空间位置成为可能。施圣贤团队[46-47]及Thurow研究团队[48]分别开发了各自的光场相机硬件系统和光场重构算法，构建了单光场相机PIV三维流动测试系统。施圣贤团队[49]成功将该技术应用于射流、湍流边界层等复杂三维流场的试验研究。相对于Tomo-PIV而言，光场PIV在沿测量体厚度方向上的测量精度较低，但其最大优势在于利用单相机进行测量，系统简单，易于布置，适用于光学空间受限情况下的复杂三维流场测量。

3 Tomo-PIV具体应用

Tomo-PIV技术建立之初，主要集中应用在湍流剪切运动(如尾流、边界层、喷流等)方面，用于研究流场三维涡结构演化规律。近年来，逐渐拓展至实际工程应用，并由低速向高超声速发展。

3.1 尾流流动

置于流体中的任何物体均会在其下游形成尾迹。在尾迹中，平均速度低于自由流速度，这种速度亏损的程度直接与作用于物体上的阻力有关。不仅如此，尾迹的存在将在很大程度上影响整个上游流场，同时也将影响物体表面的压力分布，从而影响物体受到的升力。因此，分析尾迹及它们与自由流的相互作用对完整地处理物体受到的流体作用力问题具有十分重要的作用。

Tomo-PIV的第一次应用就是测量圆柱尾迹流动[17，50]。图8为Elsinga等[17]测量得到的圆柱后方卡门涡街的涡结构等值面云图，显示了该技术在不稳定分离流方面的应用潜力；Hain等[51]利用高分辨率相机的Tomo-PIV系统研究了竖直放置圆柱上表面的三维流场，获得了详细的剪切层三维结构；Ghaemi等[52]使用时间分辨率Tomo-PIV在开口风洞中测量了NACA0012翼型尾缘三维流场，得到了近尾缘区域发卡涡及发卡涡对的分布情况，提出了用于描述尾缘区域三维不稳定结构的理论模型；许相辉等[53]在低速风洞中对圆柱尾流场进行了试验测量，成功获取了圆柱后方典型的三维卡门涡结构；高琪等[54]利用自主研发的层析PIV技术实现了合成射流的三维测量，观测到双涡流动结构追逐并融合的现象，复现了三维涡环结构的时空演化过程；此外，Zhu等[55]利用6相机的Tomo-PIV系统对高宽比为2的短圆柱尾流场涡结构及动态演化特性进行了研究(如图9所示)，发现有限高圆柱的尾流场被弓形涡所主导，且第一次发现了M形状的弓形涡，并提出了有限高圆柱平均尾流场的概念拓扑模型。

图8 瞬时涡结构等值面(d为圆柱直径)[17]Fig. 8 Iso-surface of instantaneous vortex structure [17]

图9 三维流场瞬时截面图(d为圆柱直径)[55]Fig. 9 Iso-surface of instantaneous field[55]

3.2 湍流边界层

在自然界和实际工程应用中，最经常发生的流动状态是湍流，流动随时间和空间都呈现出不规则的脉动。实现湍流边界层的连续测量不仅需要较高的时空分辨率，还需要同时测量其三维速度分量。

Elsinga等[56]对低速风洞下壁面的边界层开展了Tomo-PIV测量，在测量区域前1 m处安装拌线使气流强制转捩，并采用Q准则计算显示了瞬时三维涡结构(见图10)。王晋军、高琪[57-58]团队采用Tomo-PIV对水洞中平板湍流边界层进行了测量，通过λci准则进行涡识别，统计了展向涡沿法向的变化规律，并给出了在流向-法向平面内高低速区域和横幅展向涡空间位置的关系。姜楠[59-61]团队利用Tomo-PIV在水洞中对平板湍流边界层进行了测量，发现了典型的四极子、六极子式结构，给出了以扫掠事件为中心的壁湍流相干结构局部动力学模型，并观测到了沿流向分布的由发卡涡构成的低速流体；此外，他们还利用Tomo-PIV技术对沟槽被动流动控制方法的减阻效果进行了研究，分析了边界层中流体的运动特点。施圣贤团队[62]利用单相机光场PIV技术对一个自相似的逆压湍流边界层进行了测量，得到了远、近壁面各600组瞬态三维流场，并与相同工况下的2D-PIV进行了对比，发现在近壁面测量结果吻合较好，总体误差小于0.5%，在远壁面则有翘尾现象。

图10 瞬时三维结构[56]Fig. 10 Instantaneous three-dimensional structure[56]

如图11所示，Humble[63]测量了Ma=2.1条件下激波/边界层干扰的三维瞬时结构，能在三维空间内显示出高低速流动区域的相互作用，并根据测量结果给出了超声速激波边界层干扰拓扑结构的概念模型； Ye等[64]在代尔夫特理工大学高超声速风洞中测量了来流马赫数Ma=6.5的过斜劈转捩流动过程，探讨了旋转涡对的上洗及下洗运动对动量输运的影响。

图11 激波边界层概念模型[63]Fig. 11 Conceptual model of shock/boundary interaction[63]

3.3 工程应用

随着Tomo-PIV的深入发展，不仅在湍流等基础研究方面展现了巨大优势，在实际的工程应用中也发挥了重要作用。

风洞流场速度的均匀分布是飞行器精细化研制和空气动力研究的基础保障。李晓辉等[65]利用Tomo-PIV技术对亚跨声速风洞流场速度均匀性进行了校测，来流马赫数0.6，测量马赫数均方根偏差小于0.005，并通过测量超临界翼型OAT15a的尾缘速度场比较了小肋减阻的控制效果；Stolt等[66]利用Tomo-PIV研究了低雷诺数下NACA0015翼型酒窝状粗糙前缘对流动分离及失速迎角的影响；Avallone利用时间解析的Tomo-PIV测量了[67]NACA0018翼型有/无锯齿状后缘的三维流场，发现锯齿状后缘在根部产生了沿流向的涡对，且在展向呈周期性分布，并结合声学测量探讨了锯齿的不同尺度对翼型噪声的影响。

推进器、风力涡轮机、直升机旋翼等转子系统的尾流不稳定性机理研究与设备性能、振动、噪声和结构问题直接相关，在工程应用中扮演了重要角色。意大利海洋研究中心的Felli[68]在空化水槽中利用Tomo-PIV开展了推进器尾流近场的研究，突出显示了端部旋涡附近二次细丝的盘绕机制；德国达姆施塔特工业大学Weinkauff等[69]使用8台相机重构得到了时间平均的火焰成像，研究了火焰与湍动的相互作用，增进了对燃烧过程的理解，有助于更加高效清洁的燃烧装置的研发；Peterson等[70]通过向内燃机缸内添加油滴颗粒，成功测量了火花诱导内燃机内的三维流场，对现代内燃机引擎的研发设计提供了重要的数据支撑。

图12 Tomo-PIV相关工程应用Fig. 12 Some engineering applications of Tomo-PIV

4 展 望

1) 算法的进一步优化。Tomo-PIV三维重构和三维互相关的计算耗费时间长，算法的优化从Tomo-PIV提出之初就一直是研究的热点。如何在不影响计算精度的基础上大幅度减少计算时间、提高计算效率，仍是未来研究的重点。

2) 复杂外形流动。目前Tomo-PIV仅仅在低速至高超声速风洞中实现了初步的应用，模型比较简单，大多为平板或翼型。针对诸如多段翼起降构型的缝道流动、多相流、微流动以及具有复杂外形模型的流场测量还需要更多的研究和尝试。

3) 高可靠性的流场/压力场/声场一体化研究。三维时间解析的速度场对于不稳定空气动力学和声学的研究具有重要意义。对于不可压缩流动而言，时间解析3D-3C测量的引入形成了控制方程连续性[71]，动量方程的所有变量均可以测量得到，仅有压力梯度张量未知：

(5)

采用合理的数值分析方法和相应的边界条件，即可进行压力场的三维重构。

Tomo-PIV测量得到的压力场可以用于预测固面的声场噪声。Violato等[72]对其进行了初步的研究，并进行了转捩射流的测量。目前国内外针对二维PIV技术重构压力场的方法进行了一定的研究，但是三维压力场及声场的重构方法还处于探索阶段。建立基于Tomo-PIV技术的压力场和声场重构技术，对开展流场/压力场/声场一体化研究具有十分重大的工程价值。

4) 高超声速非定常流场测量。由于Tomo-PIV采用多相机布局，相机和测量区域具有一定的角度，高超声速非定常流动的流场密度变化大，导致光线发生折射，影响示踪粒子的成像，拍摄图像和初始标定映射函数存在一定的误差，对测量体构成影响。因此，针对高超声速非定常三维流场的测量，还需要针对密度场导致的光场畸变开展进一步的探索研究。

5) 基于深度学习的PIV/Tomo-PIV研究。目前PIV/Tomo-PIV分析方法采用的是均匀移动线性假设，不仅对噪声敏感，还会出现测量异常值。将人工智能中的深度学习技术引入到示踪粒子空间分布的重构和追踪以及速度场信息提取中，提高处理的速度和精度，扩展PIV技术的应用场景，是PIV/Tomo-PIV最新的研究方向。