用于湍流研究的三维粒子跟踪测速技术：以射流流动为例

摘要

本文介绍了一种基于高速相机和四视角分光器的三维粒子跟踪测速系统（3D-PTV）。该技术应用于雷诺数 Re ≈ 7,000 时，圆形管道下游约十个管径距离范围内的射流流动研究。

摘要

3D-PTV 是一种定量流动测量技术，旨在通过立体成像序列记录，追踪一组粒子在三维空间中的拉格朗日轨迹。本文描述并讨论了由高速相机与四视角分光器构成的 3D-PTV 系统拓扑结构的基本组件、特征、约束条件及优化建议。该技术应用于雷诺数 Re ≈ 7,000 的圆形射流的中间流场区域（5 <x/d <25）。在射流出口下游约十个直径处及距射流中心轴不同径向距离的位置，估算了拉格朗日流动特征和欧拉框架下的湍流参数。拉格朗日参数包括选定粒子的轨迹、速度和加速度，以及流线曲率，这些参数由弗勒内-塞雷方程（Frenet-Serret equation）求得。将射流中心轴周围、位于射流出口下游十个直径处横截面上的三维速度场和湍流场的估算结果与文献数据进行了对比，并在距射流中心不同径向位置获得了大尺度流向速度运动的功率谱。

引言

湍流射流在工程应用中普遍存在。对这类流动进行详细表征，在从大规模环境排放系统到电子微尺度器件的广泛实际问题中至关重要。由于射流流动对众多重要应用具有显著影响，相关研究已十分深入^1-4。多种实验技术已被用于在不同雷诺数和边界条件下表征射流流动，包括热线风速仪^4-8、激光多普勒测速仪（LDV）^4,9-12以及粒子图像测速技术（PIV）^12-16。近年来，已有少数研究采用三维粒子追踪测速（3D-PTV）来研究射流流动中湍流与非湍流的交界面^17,18。3D-PTV 是一种能够从不同视角描述复杂湍流场的特殊技术，它利用多视角立体成像，在拉格朗日参考系下重建流场体积内的粒子运动轨迹。该技术最初由 Chang¹⁹ 提出，随后由 Racca 和 Dewey²⁰ 进一步发展。此后，3D-PTV 的算法和实验装置已得到多项改进^21-24。基于这些进展和前期工作，该系统已成功应用于研究多种流体现象，例如在 4 m × 2 m × 2 m 域内的大尺度流体运动²⁵、室内气流场²⁶、脉动流²⁷以及主动脉血流²⁸。

三维粒子追踪测速（3D-PTV）测量的工作原理包括数据采集系统的搭建、记录/预处理、标定、三维匹配、时间序列追踪以及后处理。精确的标定能够实现对粒子位置的精确定位。通过三个以上视角图像中检测到的粒子对应关系，可基于对极几何重建出粒子的三维位置。连续图像帧之间的关联形成时间序列追踪，从而确定粒子的运动轨迹 s(t)。优化3D-PTV系统对于最大限度提高多粒子可追踪性至关重要。

优化的第一步是获取合适的数据采集系统，包括高速相机、照明光源以及示踪粒子的特性。相机分辨率与 interrogation volume 的尺寸共同决定了像素大小，因此也决定了所需示踪粒子的尺寸，其尺寸应大于单个像素。通过以亮度的平方作为权重对粒子像素的位置取平均值，可实现对检测到的粒子质心的亚像素精度估计²¹。相机的帧率与雷诺数及粒子追踪的关联能力密切相关。更高的帧率有助于解析更快的流动或更多数量的粒子，因为当图像间粒子的平均位移超过粒子间的平均间距时，追踪将变得更加困难。

在图像采集过程中，快门速度、光圈和感光度是需要考虑的三个因素。快门速度应足够快，以尽量减少粒子周围的模糊，从而降低粒子质心位置的不确定性。相机光圈应调整至与测量区域的景深相匹配，以降低检测到该区域外粒子的概率。由于相机的最大感光度是固定的，因此随着帧率的增加，用于照亮粒子所需的光照强度也应相应增加。与粒子图像测速（PIV）不同，三维粒子跟踪测速（3D-PTV）并不严格要求复杂的光学设置和高功率激光器，只要光源能从示踪粒子充分散射到相机即可。连续发光的LED或卤素灯是具有成本效益的良好选择，且无需同步装置²¹。

在三维粒子追踪测速（3D-PTV）中，与其他光学流场测量技术类似，示踪粒子的速度被假定为局部瞬时流体速度²⁹然而，这种情况仅适用于直径和惯性均为零的理想示踪粒子；实际的示踪粒子应具有足够大的尺寸，以便能被相机捕捉。有限尺寸粒子的示踪保真度可通过斯托克斯数（Stokes number）来确定。 S_t, 即 粒子的弛豫时间尺度与所关注的湍流结构时间尺度之比。通常情况下， S_t 应明显小于1。对于 S_t  ≤0.1 的流动追踪误差低于 1%³⁰详细讨论可参见 Mei 等.^29-313D-PTV 实验推荐的示踪粒子尺寸取决于光源和相机的灵敏度。当使用卤素灯或LED灯作为照明光源时，通常采用相对较大的粒子（例如 50-200 µm³²，而较小的颗粒（例如 1-50 µm^33,34 可与高功率激光器一起使用（例如 80-100 瓦连续波激光器）。对于特定波长的光具有高反射率的粒子，例如在卤素光下带有银涂层的粒子，可将其标记放大为图像。示踪粒子的 seeding density（播种密度）是实现成功三维PTV测量的另一个重要参数。粒子数量过少会导致轨迹数量不足，而粒子密度过高则会引起对应关系建立和追踪过程中的歧义。对应关系建立中的歧义包括粒子重叠以及在定义的极线（epipolar line）上检测到多个候选点。在追踪过程中，由于粒子平均间距较短，高播种密度会导致追踪歧义的发生。

第二步是优化记录/预处理过程中的设置以提高图像质量。摄影参数设置，例如增益 & 黑电平（G&B）在优化图像质量方面起着重要作用。黑电平定义了图像最暗部分的亮度水平，而增益则放大图像的亮度。G 的微小变化&B 水平会显著影响可追溯性的可能性。事实上，较高的 G&B 可能使图像过亮，并最终损坏相机传感器。为说明这一点，以下阐述 G 的影响&本文还研究了#38;B水平对流场重建的影响。在预处理步骤中，图像通过高通滤波器进行滤波，以增强来自粒子的光散射信号。调整像素尺寸和灰度级，以在 interrogation volume 内最大化粒子检测效果。

优化的第三步是立体成像的精确校准，该过程基于对极几何、相机参数（焦距、主点和畸变系数）以及折射率变化。此步骤对于最小化基准目标点的三维重建误差至关重要。对极几何利用相机与测量区域之间的相对距离以及目标图像的倾斜角度。沿相机视场穿过测量区域时的折射率变化可根据Mass提出的方法予以考虑。 等²¹在本实验中，使用具有规则分布靶点的三维阶梯状结构作为靶标。

在三维粒子追踪测速（3D-PTV）实验中，尽管仅需两幅图像即可确定粒子的三维位置，但通常会使用更多的相机以减少不确定性²¹。一种替代昂贵多高速相机配置的方案是视场分光器，该装置由 Hoyer 等³⁵ 提出用于 3D-PTV，最近由 Gulean 等²⁸ 应用于生物医学领域。视场分光器由一个金字塔形反射镜（以下简称主反射镜）和四个可调节反射镜（以下简称次级反射镜）组成。本研究采用一个四视场分光器和单台相机，模拟四台相机的立体成像效果。该系统用于表征管射流在射流起始位置下游约 14.5–18.5 倍管径处的中间流场，射流直径，d_h = 1 cm，雷诺数 Re ≈ 7,000，分别从拉格朗日和欧拉参考系进行分析。

方案

1. 实验室安全

1. 查阅所选照明光源的安全使用指南（例如 激光、工业LED、卤素灯
   注意：本实验使用一组五个250瓦卤素聚光灯作为光源。该光源的基本安全事项和使用建议如下所述。
   1. 避免直接接触卤素灯，因其工作温度较高（约3,000 K色温）。
   2. 仅在采集数据时开启光源，以避免加热待测流体。
   3. 远离光源处放置所有易燃材料，包括各种纸张。

2. 实验设置

1. 选择合适的镜头
   1. 选择像差较小的镜头，以避免校准问题。推荐使用远摄镜头或微距镜头。
   2. 通过估算所需放大倍数 M，确保镜头在物距 O 范围内覆盖所需的视场（FOV）。
      注：放大倍数为相机芯片长度与视场之比，物距可由公式 O = f(1/M + 1) 计算得出，其中 f 为镜头焦距。本实验中，相机芯片长度为 20.34 mm，对应的视场（即主反射镜）为 50 mm，且物距受限于 O ≤ 250 mm（由于安装相机和视图分束器的滑轨长度有限）。因此，放大倍数 M = 20.34/50 = 0.41，在给定物距范围内对应的近似焦距为 f ≤ 72.7 mm。故本实验采用焦距为 60 mm、光圈比为 f/2.8D 的微距镜头。
2. 安装并调节带视图分束器的相机
   1. 通过沿垂直支架滑动主反射镜并用支架固定件将其固定，使主反射镜中心与测量区域中心对齐。注意：此步骤应在安装次级反射镜之前完成。
   2. 安装相机，并使图像中心与主反射镜中心重合。
   3. 通过控制采集软件中的感兴趣区域（ROI）设置，调节相机视图的宽度和高度，使其刚好覆盖主反射镜。此过程可减小图像尺寸和图像噪声。注：本实验中，主反射镜尺寸为 5 × 5 cm²（对应 1,728 × 1,728 像素）。
   4. 设计一个（定制的）三维校准靶标。该靶标应完全包含整个测量体积。确保分束器的每个视图都能捕捉到所有靶标标记，以便实现均匀校准。
      注：在本演示中，靶标使用黑色 Vero 塑料 3D 打印而成，呈阶梯状，尺寸为 35 × 35 × 30 mm³，白色靶点直径为 1 mm，在垂直、流向和展向方向上的间距分别为 2.5 mm、5 mm 和 10 mm。靶标的精确几何结构至关重要，因其将用于构建校准模型，并确定相机相对于水槽的位置。
3. 将校准靶标放置于测量区域内
   1. 将校准靶标置于一个可调高度的平台上，并使其面向相机。
   2. 通过调节靶标平台的高度，使校准靶标中心与测量区域中心对齐。
      注：在本示例中，校准靶标中心与射流喷嘴中心对齐，高度为 20 cm。可使用气泡水平仪来校准靶标水平。
4. 安装并调节四视图分束器的次级反射镜
   1. 将主反射镜置于距离测量区域一定距离处，以确保其完整成像。本演示中该距离为 0.2 m（图 1）。
   2. 将次级反射镜大致安装在合适位置，使来自主反射镜两侧的相机视图大致对准相应的次级反射镜，并通过将其固定在次级反射镜的垂直支架上来固定位置。
   3. 对另外三个反射镜重复此过程。检查所有次级反射镜相对于主反射镜的几何对称性。
   4. 通过调节次级反射镜的镜座进行最终调整，确保四个视图均完整包含整个校准靶标。一种高效检查反射镜位置和角度的方法是使用激光笔可视化每个视图的成像光路。
   5. 通过移动一个反射镜来检查子图像之间是否存在重叠。
      注：如果仅有一个视图发生变化，则重叠区域可忽略不计；否则，重复步骤 2.4.2 至 2.4.5，直至重叠区域最小化。
5. 将光源直接面向测量区域放置。调节光源时务必遮盖相机，以防损坏相机传感器。
   1. 检查光源是否在整个测量体积内均匀分布。
   2. 如有必要，可在光源正下方放置一个放大镜以增强光照强度。注：本实验中使用焦距 f₀ = 450 mm 的平凸放大镜以增强照明。

3. 设置优化

1. 打开相机并调整设置以提高图像质量。
   1. 调节镜头放大倍数，直至通过主反射镜观察时，四个次级反射镜视图中的反射像均能同时清晰聚焦。
   2. 通过观察分光器四个视图中标定目标图像的对称性，检查各视图图像是否对称，并确认是否覆盖了 interrogation volume（ interrogation 体积）。
   3. 调节f-number，使相机能够同时清晰捕捉到距离相机最近和最远的标定目标点。
      注意：这可确保相机仅在 interrogation volume 的深度范围内捕捉示踪粒子。本示例中，f-number 设置为 11。
   4. 根据实验需求设置所需的帧率为 550 Hz（具体值取决于实际应用，参见引言部分），并在记录软件中相应地最大化光敏感度。
2. 通过实时相机画面观察分光器各视图中粒子密度的差异，检查主反射镜每个视图的照明情况。
   注意：若使用多个光源照射 interrogation volume，则分光器各视图之间可能出现照明差异。本实验中，由于光源位于上方，上方两个次级反射镜接收到的光线较少。在水槽底部使用平面反射镜有助于减少各视图间的光照差异。
3. 在使用 3D-PTV 光源前，关闭实验室内背景光源。
4. 调节相机的 G&B 水平，以更好地捕捉粒子散射光。录制若干段不同 G&B 水平下的短序列，通过观察粒子轨迹的分布与密度，确定最佳设置。
   注意：本实验中 G&B 水平范围为 0–500，黑电平（B）设为 500，以增强较弱的散射光亮度；增益（G）设为中等值 300，以适度放大图像信号并避免图像过亮。

4. 校准

1. 在添加示踪粒子之前，将标定目标放置于检测区域内，并拍摄几张标定图像。使用较弱的光源（例如 LED手电筒）照亮标定目标。
2. 将标定图像划分为四个独立的子图像，并创建一个文本文件，记录标定目标标记点的参考坐标位置。此处使用OpenPTV软件（http://www.openptv.net）完成此步骤。
   注意：后续处理过程与使用多相机系统的用户完全相同。
3. 将第4.2步中获得的图像和文本文件保存至软件的“Cal”文件夹后，点击“Create calibration”选项卡以启动标定过程。
4. 点击“Edit calibration parameters”选项卡，并选择“Calibration Orientation Parameters”选项卡，以定义放大倍数、旋转角度以及每个分视图中心与标定目标原点之间的距离。
   注意：第一行为标定目标原点到相机传感器在x,y,z方向上的距离；第二行为绕x,y,z轴的旋转角度（单位：弧度）；接下来是一个3×3的数据矩阵，表示旋转矩阵；随后两行为针孔在x和y方向上相对于图像中心的距离（单位：mm）以及焦距；最后一行为水槽玻璃相对于标定目标原点在x,y,z方向上的位置。
5. 点击“Detection”和“Show initial guess”，检查“初始估计”点是否与检测到的目标点匹配。
6. 对全部四个视图重复第4.4步，直至“初始估计”点与各标定图像中的对应点对齐。
7. 点击“Orientation”以重建检测区域的空间取向。
   注意：可通过调整镜头畸变和仿射变换进一步优化标定结果。至此，检测区域已完成标定，可进行数据处理。有关标定过程的更多详细说明，请参见作者的学位论文³⁶。

5. 流速设置/数据采集

1. 根据帧率和最大流速估算每个相机视场中捕获的粒子最大数量。在本示例中，参考流速为 U ≈ 0.4 m/s，帧率为 550 Hz， interrogation volume 约为 4 × 4 × 4 cm³，每帧图像中约包含 1,000 个粒子。
2. 使用步骤 3 中获得的优化设置开启相机。
3. 加入示踪粒子，并等待数个平均停留时间，使流动达到稳态。如需可补充添加粒子，但应避免过高的粒子密度（在步骤 5.1 中已估算），以免引起歧义。
   注：本示例中，使用约 1.6 g 的银涂层空心陶瓷微球（粒径 100 µm，密度 1.1 g/cm³）作为流体介质（2 × 0.4 × 0.4 m³）的示踪粒子。
4. 记录所需数量的流动图像。
   注：本实验中，使用采集软件以 550 Hz 的帧率采集了 9,000 张图像。若相机和/或视场分光器发生移动（即使轻微移动也可能严重影响结果），需重复步骤 2.4 至 5.3。

6. 数据处理（通过 OpenPTV 软件）

1. 将步骤 5.4 中获得的原始图像分割为四个独立的子图像。
2. 在“开始”选项卡下点击“初始化/重启”，以从四个视角加载初始图像。
3. 右键单击“Run”目录，然后点击“主参数”，以控制相机数量、折射率、粒子识别、序列图像数量、观测体积以及对应关系的判定标准。
   1. 在“常规”选项卡下定义实验所用的相机（视角）数量。本实验中，将相机数量设为 4。
   2. 在“折射率”选项卡下定义沿相机视角方向的折射率。
   3. 在“粒子识别”选项卡下设定像素检测的最小值和最大值，以及灰度值阈值，以优化四个视角中的粒子检测数量。像素检测的最小值、最大值及灰度阈值决定了用于粒子检测的像素尺寸和亮度水平，可消除噪声及失焦粒子。
   4. 在“序列处理参数”选项卡下定义需处理的图像数量。
   5. 在“观测体积”选项卡下定义观测体积。
   6. 在“对应关系标准”选项卡下定义对应关系的相关性，包括用于立体匹配的总带宽参数（mm）。
4. 在“预处理”选项卡下点击“高通滤波”。该操作可增强四个视角中粒子的光散射信号。
5. 点击“粒子检测”，以在亚像素级别确定四个视角中被检测粒子的质心。重复步骤 6.2 和 6.3，直至检测到的粒子数量接近步骤 5.1 中计算出的预期粒子数量。
6. 点击“对应关系”，以在各视角中建立立体对应关系。
   注意：为重建检测粒子的三维位置，至少需从三个视角确定对应关系。
7. 点击“3D 位置”，基于标定结果获取检测粒子的三维坐标。
8. 点击“无显示序列处理”，对所有图像序列重复执行步骤 6.4 至 6.7。
   注意：此过程将为每组图像生成一个“rt_is”文件，其中包含帧内检测粒子的汇总信息，采用制表符分隔的文本文件格式。
9. 右键单击“Run”目录，然后点击“跟踪参数”，以定义用于搜索跟踪候选粒子的球体半径参数（例如 dvxmin 和 dvymin，单位为 mm/帧）。
10. 点击“无显示跟踪”，为步骤 6.7 中重建的粒子定义粒子识别编号（ID）。
    注意：该过程采用四帧预测-校正方案²⁴，关联相邻帧序列以实现粒子跟踪。此步骤将为每组图像生成一个 ptv_is 文件，其中包含帧内检测粒子的跟踪信息；前两列分别显示粒子在前一帧和后一帧中的 ID。
11. 点击“显示轨迹”，以在每个相机视角中可视化粒子运动轨迹。

7. 后处理（可选）

注意：后处理的范围和类型取决于个人需求，因此可进行定制。此处以点基计算为例进行简要说明。

1. 通过拉格朗日坐标系获取数据（使用 Matlab）。
   1. 从 ptv_is 文件中提取每个粒子的三维位置及其对应的 ID。该步骤可实现图像序列中检测到的粒子之间的关联，从而重建粒子轨迹。
   2. 根据给定的帧率，计算每条轨迹中粒子的速度和加速度。在本示例中，通过对位置信号应用移动三次样条低通滤波来计算粒子的速度和加速度^34,37。
   3. 构建一个结构体数组，其字段包含每条轨迹的三维位置、三维速度、三维加速度、时间戳以及轨迹 ID。在此数据格式中，struct 数组的长度表示轨迹的数量。
2. 通过欧拉坐标系获取数据（使用 Matlab）。
   1. 利用每个粒子的时间戳，将结构体数组（步骤 7.1.3）转换为时间序列形式。这将生成与步骤 7.1.3 中类似的结构体数组结构，但此时结构体数组的长度表示帧数，在本实验中为 9,000。
   2. 将每个时间帧的时序结构体数组插值到三维网格上，以获得欧拉坐标系下的瞬时速度场。在本示例中，使用 Matlab 中的 griddata 函数进行插值。

结果

装置的照片和示意图见图 1 和图 2。校准靶标、视场分光器上反射的基准标记以及三维校准重建结果如图 3 所示。在流向 x、展向 y 和深度 z 方向上，识别出的校准靶标的均方根误差（RMS）分别为 7.3 µm、5.7 µm 和 141.7 µm。在z 坐标方向上相对较高的 RMS 误差是由于该方向上的靶标点数量较其他方向少，且四个视角相对于 z 轴的夹角较小，与 x 和 y 坐标相比分辨率较低所致。在任意时刻，每个四视角图像中检测到的粒子数量约为 10³ 量级。在所有被检测到的粒子中，由于仅有位于交叠区域内的粒子才能被成功重建，因此成功的三维重建粒子数量减少至大约一半。视频 1 展示了来自四视角分光器的射流流动的高速视频样本。

讨论

3D-PTV 在揭示多种湍流的复杂物理机制方面具有巨大潜力，例如近地面大气中的大尺度湍流运动²⁵，室内空气分布²⁶或主动脉拓扑结构中的脉动流²⁸ 以及其他众多因素。然而，充分理解其优势与局限性并积累相关经验，对于最大化发挥其潜力至关重要。为确定最佳设置（包括帧率、照明光源、G）而进行的试错性初步测试和 exhaustive iterations（详尽迭代）是必不可少的。&B 水平和图像滤波方法与重建一组（的拉格朗日路径的能力直接相关例如，示踪）颗粒。需要注意的是，关键的实验步骤（如本实验所示）是调整G&B 水平以及视场的照明情况（包括卤素点光源、放大透镜和水槽底部反射镜的组合）。

这些调整有助于将光散射优化至四个视图的检测范围内。在确定了高保真测量的实验参数后，应进行全面的调整和故障排查，以根据帧率、相机分辨率和检测区域体积的大小，计算出可准确追踪的粒子轨迹的最大数量。尽管提高帧率可以增加捕获的粒子数量，但需要注意的是，3D-PTV 中可追踪的粒子数量远低于PIV。3...

材料

本文使用的材料清单

姓名	公司	目录编号	评论
Mikrotron 4CXP 4 通道 CXP-6 CoaXPress	ImageOps	CAMMC4082	高速相机
Active Silicon FireBird CoaX 帧采集卡	ImageOps	FBD-4XCXP6	帧采集卡
100 μm 镀银中空陶瓷微球	Potters Industries LLC	AG-SL150-30-TRD	示踪粒子
StreamPix6	Upstate Technical Equipment CO.,INC	MISNOR-STP-6-S-CL	相机应用软件
四视角分光器	Photrack AG		定制部件，若使用单台相机进行三维粒子追踪测速（3D-PTV）则为必需
250 瓦卤素聚光灯	General Electrics	23719	光源
OpenPTV（软件）	OpenPTV（http://www.openptv.net）		开源粒子追踪软件（注：可为希望免于繁琐安装过程或缺乏足够计算资源的用户提供服务形式的使用支持）