2. 上海交通大学 机械与动力工程学院, 上海 200240
2. School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China
在许多自然界和工程领域的研究中,例如鱼在游动时改变自身形态来减小运动阻力、血液在细小血管中的脉动运动、航空发动机内叶栅和高温气体发生剧烈的耦合振动现象[1, 2, 3]等,精确测量复杂固体壁面附近的三维速度场,尤其是流固耦合实验研究中同步测量三维流场信息和固体结构表面的运动/变形特征,是剖析这些复杂流动现象物理本质的关键。激光粒子图像测速技术(Particle Image Velocimetry,简称PIV)由于具有可视化、全场性和非接触测量等特点，目前被广泛应用于固体近壁及远场的流场实验研究[4]。然而，常规PIV算法由于不能准确区别固体边界和示踪粒子，导致在近壁区域的PIV计算中不能精确划分互相关窗口，从而增大了这些区域的速度场计算误差。
针对这一问题，近年来国内外学者对任意物体二维边界的识别及其边界附近流场的测量进行了大量的研究。Jeon和Sung[5]提出了基于Texton卷积的图像处理算法，以精确获取任意运动变形边界的二维信息；同时根 据边界位置信息，对流场粒子图像进行变形，从而提高了传统矩形PIV互相关窗口算法的计算精度。随后曹洪才等[6]采用基于Radon变换的滑动窗口边界识别方法，准确地获得了任意运动变形边界的位置信息，并应用贴体自适应互相关计算窗口测得瞬态流场信息，提高了任意二维边界附近速度场的测量精度。最近Adhikari和Longmire[7]对PIV实验中物体三维边界的识别进行了尝试，采用Visual hull算法对球形、方形等简单模型进行轮廓识别，但最终能识别的仅是物体的三维局部轮廓，识别精度尚待提高。
为了精确识别任意固体边界的三维形状，准确测量固体边界附近三维速度场，本文基于SURF(Speeded Up Robust Features,加速稳健特征)模式识别算法[8]，对三维流场测量中2个相机从不同视角同时拍摄的物体图像进行特征提取和匹配，重构出三维物体的复杂曲面形态，从而获取任意三维物体边界信息。此外，基于MLOS-SMART(Multiplicative Line Of Sight-Simultaneous Multiplicative Algebraic Reconstruction Technique)三维粒子场重构[9]及三维互相关算法计算获得任意三维物体边界附近的空间速度分布。本文首先采用已知不同曲率的圆柱图像进行系统的验证分析，证明了基于SURF模式识别算法的准确性；然后利用计算流体力学(CFD)计算的近壁圆柱绕流速度场，结合四相机标定矩阵，数字合成4组PIV粒子图像对。将提出的新算法应用于此数字 合成Tomo-PIV图像序列，最终准确地获得了圆柱边界的三维形态以及边界附近的瞬态速度场。通过将计算结果与已知圆柱曲率和CFD计算的三 维速度场对比，验证了基于SURF模式识别和MLOS-SMART三维粒子场重构及三维互相关算法的正确性。 1 复杂固体边界流场的三维PIV测试方法 1.1 算法基本流程
三维曲面结构及其流场同步测量算法的计算流程及相关实验系统布置示意图如图 1所示。实验中4个相机将同时记录流场中三维物体与示踪粒子的原始图像I(i,j)，要实现对物体曲面结构和三维速度场的同步测量，首先需要从原始图像中分离出示踪粒子图像和三维曲面图像，即
式中：P(i,j)表示示踪粒子图像，S(i,j)表示三维曲面图像。
 |
| 图 1 实验系统布置示意图及算法流程图Fig. 1 Schematic of experimental system and flow chart of the algorithms |
SURF模式识别算法要求在三维物体表面增加随机分布的人工标记点，标记点直径一般远大于示踪粒子的直径。根据这一特征，可利用二维高斯光滑滤波器(公式2)实现物体图像和粒子图像的分割。
令二维高斯函数的标准差σ大于示踪粒子的直径，可以使图像模糊，滤掉图像中的示踪粒子，从而得到曲面图像S(i,j)。然后利用公式(1)得到示踪粒子图像P(i,j)，实现图像分离。然后，利用分离后相机2，3中的曲面图像，采用SURF模式识别算法便可重构出三维曲面的几何特征。另一方面，利用分离后相机1、2、3和4中的示踪粒子图像，基于MLOS-SMART三维粒子场重构及三维互相关算法计算三维速度矢量场，获得瞬态三维流场。
1.2 基于SURF模式识别算法的三维曲面重构算法SURF算法是图像处理和计算机视觉领域的一个稳健的图像识别和描述算法，由Bay等于2006年提出[8]。SURF算法基于近似的离散小波变换响应，利用卷积积分图像，对图像进行特征提取和描述。它可以处理同一物体两幅图像之间发生平移、旋转、仿射变换等情况下的特征匹配问题，具有很强的匹配能力和很快的匹配速度。针对同一物体两张不同视角的图像，SURF算法首先基于图像Hessian近似矩阵行列式的最大值特性，来提取图像中具有显著特征的关键点。一般关键点处的图像梯度很大，在匹配过程中具有独特性。然后对各个关键点处的特征用一个以关键点附近的图像梯度直方图分布所代表的特征向量进行特征描述，再利用最近邻匹配准则对各个关键点的特征向量进行匹配。所有的匹配利用对极线限制去除错误匹配后，即可得到两张图像关键点的正确匹配。对匹配后的关键点，利用三角化Triangulation方法便可重构出其在空间中的分布[10]。最后利用Delaunay triangulation三角剖分算法[11]对空间点云进行三维曲面重构，得到三维曲面的空间形态。SURF算法的基本流程如图 2所示。
 |
| 图 2 基于SURF模式识别的三维曲面重构算法基本流程Fig. 2 SURF based three-dimensional boundary recognition algorithm |
本文中三维瞬态速度场测量基于三维层析流场测试技术(Tomographic PIV，Tomo-PIV)。Tomo-PIV由Elsinga等[12]于2006年提出，该技术实现了三维流场测量中三维三分量速度场的精确测量，提高了PIV技术对复杂三维流场的定量测量能力。典型的Tomo-PIV系统由4个相机和脉冲激光组成，多相机系统根据Scheimpflug准则[13]，均对焦于实验测量区域，实验区域由激光通过柱棱镜产生的体光源照亮。在处理Tomo-PIV粒子图像之前，需首先对相机系统进行标定，以确定相机平面和流场测量空间的映射关系；同时需采用体自标定算法以减小相机标定误差[14]。获取准确的相机标定矩阵后，即可采用SMART或者MLOS-SMART算法[9]对粒子图像进行重构(本文采用更加高效的MLOS-SMART算法)以获取相邻两时刻的三维粒子图像，最后通过三维互相关算法处理重构后的粒子图像体，便可得到瞬态三维速度矢量场。 2 算法验证及应用 2.1 算法验证 2.1.1 实验设置和相机标定
为验证SURF模式识别算法的准确性，本文采用图 3(a)、(b)所示的实验系统分别采集40、60和80mm(见图 3(c))3种圆柱曲面2个不同视角的图像，并将计算获得的结果与圆柱已知曲率进行对比。实验中采用大功率卤素灯作为光源，圆柱不同视角的图像分别由2个互成60°夹角 的MIKROTRON EoSens MC1362(1024pixel×1280pixel)相机采集，实验测量区域的几何尺寸为100mm×100mm×50mm。
 |
| 图 3 相机布置示意图、实验系统图及带有人工标记点的圆柱模型Fig. 3 Sketch of experimental system,two-camera experimental system and circular cylinder with artificial markers |
实验的第一步是相机标定，相机标定的目的是确定每个相机的三维投影变换矩阵，即相机拍摄图像上 的二维平面坐标和空间中三维空间坐标的对应关系。
相机模型采用针孔模型(Pinhole model)，每个相机的投影变换矩阵可以由一个34的矩阵表示[15]，运用直接线性变换法(Direct linear transformation)，投影关系如公式(3)所示，
式中：(u,v)代表相机拍摄图像的二维坐标，λ代表尺度因子，(X,Y,Z)代表空间三维真实坐标，A 代表投影变换矩阵。
实验中使用的相机标定靶板为一个互成90°的不共面V形标定板，如图 4(a)所示，标定板上有等间距排列、已知空间三维坐标的标定控制点。相机拍摄获得标定板的图像后，即可采用阈值法将标记点图像与黑色背景进行区分，同时使用二维高斯函数拟合标记点图像的中心位置，图 4(b)是标定过程中被识别出的标记控制点。标定过程中，每个相机拍摄的标定图像都有超过300个标定控制点被识别出来。标定过程使用最小二乘法拟合相机投影变换矩阵，若最后计算投影误差的均方根误差小于1 pixel，则相机标 定成功。实际实验中相机系统各相机投影误差的均方根误差分别为：0.80、0.78、0.77和0.87pixel，满足标定要求。
 |
| 图 4 V型标定板及标定过程中被识别的标记点Fig. 4 V-shape calibration target and detected markers in the calibration images |
采用SURF重构算法对3种不同曲率圆柱曲面图像处理的结果如表 1所示，重构后的三维曲面形态和原始曲面基本一致，曲率大小的误差均低于10%(见表 1)，证明计算结果正确。
| 编号 | 实际圆柱面的 曲率/mm-1 | 三维重构曲面的 曲率/mm-1 | 误差 |
| 圆柱1(40mm) | 0.025 | 0.02368 | 5.28% |
| 圆柱2(60mm) | 0.01667 | 0.01675 | 4.79% |
| 圆柱3(80mm) | 0.0125 | 0.01142 | 8.64% |
在3种圆柱曲面分析计算中，误差最小的是直径 为60mm的圆柱，其图像处理过程如图 5所示：(a)为SURF模式识别算法提取的关键点特征分布，特 征描述以关键点为中心不同半径的圆来代表，(b)代 表的是2张不同视角图像关键点的特征匹配，(c)对匹配后的特征点进行三角化重构后获得的点云空间分布。图 6所示为利用Delaunay triangulation三角剖分算法重构计算出的三维曲面，曲面已经过光滑处理，图中颜色代表曲面的空间坐标大小。
 |
| 图 5 直径60mm圆柱的图像处理过程：SURF算法提取关键点、关键点的特征匹配及重构点的空间分布Fig. 5 Image processing of circular cylinder (d=60mm): detected keypoints by SURF，feature matching of keypoints and distribution of reconstructed points in space |
 |
| 图 6 三维曲面重构得到的三维曲面Fig. 6 Three dimensional surface by surface reconstruction |
根据CFD计算获得的圆柱绕流速度场数字合成示踪粒子图像，同时采用图 4(b)中相机标定过程的布置方式实际拍摄圆柱曲面图像，将这2组图像叠加作为原始图像来检验基于SURF模式识别算法和三维速度场测量方法的计算精度。
首先，按照图 4(b)所示的四相机系统标定布置方式，对圆柱靠近尾流区域的曲面进行图像采集，得到4个视角的圆柱尾缘曲面图像，其中左右相机位于同一水平线，呈60°夹角，上下相机在同一竖直线，呈30°夹角；然后，利用四相机标定矩阵人工合成粒子图像。示踪粒子图像数字合成的基本步骤简述如下：首先利用CFD计算获得近壁圆柱绕流的三维速度场分布，所采用的计算几何模型如图 7所示。CFD计算参数设定为：自由来流速度U=1m/s，圆柱的直径D=66mm，雷诺数Re=UD/υ=6540，圆柱侧面离壁面距离为1D，计算采用SST模型，边界条件设置为：进口速度条件，出口压力条件。计算获得CFD速度场后，在图 7所示的三维PIV计算区域内随机生成t0时刻的粒子三维空间坐标(x0,y0,z0)，根据CFD三维速度场插值获得每个粒子的瞬态速度。在给定的 时间间隔内(Δt=0.0005s)，根据每个粒子的瞬态速度及其在t0时刻的位置，计算获得粒子在t1=t0+Δt时刻的三维空间坐标(x1,y1,z1)。最后根据多相机系统的4个标定矩阵，分别将粒子在t0和t1时刻的空间位置投射至4个相机的CCD图像平面，计算获得粒子图像在4个相机平面的像素坐标。根据粒子图像中像素灰度值按高斯分布的规律及粒子图像直径为3 pixel的设定，即可生成4个相机在2个时刻所拍摄到的粒子图像，并且设定示踪粒子图像中粒子浓度为0.05particle/pixel。为简化计算，数字合成图像没有考虑相机噪声的影响。最后，将数字合成PIV粒子图像和近壁圆柱尾缘曲面图像叠加，形成最终的数字合成图像，如图 8所示。
 |
| 图 7 近壁圆柱绕流几何模型示意图Fig. 7 Geometric model of near-wall flow of a cylinder wake |
 |
| 图 8 图像叠加过程：圆柱图像、示踪粒子图像及合成图像Fig. 8 Image addition process: surface image,particle image and synthetic image after image addition |
数字合成的4组粒子图像分别通过公式(1)和(2)进行图像分离，得到三维曲面图像和示踪粒子图像。然后根据三维边界识别算法从三维曲面图像中重构出流场中的圆柱尾缘曲面；另一方面，采用基于MLOS-SMART粒子场重构和三维互相关算法的三维PIV算法对示踪粒子图像进行处理，三维互相关 算法计算时，互相关判读窗口大小为64pixel× 64pixel×64pixel，窗口重叠为0，最终计算获得具有三维物体边界的三维速度矢量场，如图 9(a)所示。将计算所得流场与CFD流场比较发现，三维瞬态速度场分布与CFD计算结果基本一致(见图 9(b))，图 9(b)中所示为三维PIV计算区域圆柱展向中心截面处沿流向方向速度及垂直流向方向速度的分布，所计算的各空间点上速度的相对误差在5%之内，证明三维流场处理算法正确。
 |
| 图 9 具有三维曲面边界的圆柱绕流三维速度云图及三维PIV算法计算结果和CFD结果的对比 Fig. 9 Three-dimensional velocity contour of circular cylinder wake with surface boundary and result comparison of 3D-PIV and CFD calculation |
采用二维高斯滤波器模糊PIV粒子图像的方法，分离了流场中固体边界和示踪粒子图像。针对不同视角所拍摄的固体边界图像，提出了一种基于SURF模式识别算法，来精确测量任意三维物体边界的三维信息；针对多视角粒子图像，本文采用基于MLOS-SMART粒子场重构及三维互相关的三维Tomo-PIV算法，计算获得三维瞬态流场。文章首先采用基于SURF模式识别算法对已知3种不同曲率的圆柱曲面进行系统的计算分析，发现所获得的三维重构曲面与已知曲面形态基本一致，验证了所提出算法的合理性。然后利用由CFD计算获得的近壁圆柱绕流流场数字生成四相机PIV图像序列，证明了本文所发展的SURF模式识别算法能很好地结合三维Tomo-PIV算法，同步获取固体边界的三维曲面形态及其固体边界附近的流场分布。基于SURF模式识别算法的任意三维边界识别算法和三维Tomo-PIV算法，为复杂固体边界流场测试以及流固耦合现象的实验研究提供了一种有效的实验测量技术。
| [1] | Dowell E H, Hall K C. Modeling of fluid-structure interaction[J]. Annual Review of Fluid Mechanics, 2001, 33(1): 450-480. |
| [2] | Thomas J P, Dowell E H, Hall K C. Nonlinear inviscid aerodynamic effects on transonic divergence, flutter, and limit-cycle oscillation[J]. AIAA journal, 2002, 40(4): 638-646. |
| [3] | 陈焕龙, 李得英, 俞建阳, 等. 缝隙扩压叶栅近壁流场与流动损失实验[J]. 推进技术, 2013, 34(3): 332-338. Chen Huanlong, Li Deying, Yu Jianyang, et al. Experimental investigation of near-wall flow and flow loss in compressor cascades with slots injection[J]. Journal of Propulsion Technology, 2013, 34(3): 332-338. |
| [4] | Adrian R J. Particle-imaging techniques for experimental mechanics[J]. Annual Review of Fluid Mechanics, 1991, 23(1): 261-304. |
| [5] | Jeon Y J, Sung H J. PIV measurement of flow around an arbitrarily moving body[J]. Experiments in Fluids, 2011, 50(4): 787-798. |
| [6] | 曹洪才, 施圣贤, 刘应征. 任意运动变形边界流场测量PIV算法的研究: 薄膜涡激振动[J]. 实验流体力学, 2013, 27(6): 58-63. Cao Hongcai, Shi Shengxian, Liu Yingzheng. Improved PIV algorithm for flow around arbitrarily moving and deforming boundary: vortex-induced vibration of flexible membrane[J]. Journal of Experiments in Fluid Mechanics, 2013, 27(6): 58-63. |
| [7] | Adhikari D, Longmire E K. Visual hull method for tomographic PIV measurement of flow around moving objects[J]. Experiments in Fluids, 2012, 53(4): 943-964. |
| [8] | Bay H, Tuytelaars T, Van Gool L. SURF: Speeded up robust features[C]. Computer Vision–ECCV 2006, Springer. 2006: 404-417. |
| [9] | Atkinson C, Soria J. An efficient simultaneous reconstruction technique for tomographic particle image velocimetry[J]. Experiments in Fluids, 2009, 47(4-5): 553-568. |
| [10] | Hartley R I, Sturm P. Triangulation[J]. Computer Vision and Image Understanding, 1997, 68(2): 146-157. |
| [11] | Lee D T, Schachter B J. Two algorithms for constructing a Delaunay triangulation[J]. International Journal of Computer & Information Sciences, 1980, 9(3): 219-242. |
| [12] | Elsinga G, Scarano F, Wieneke B. Tomographic particle image velocimetry[J]. Experiments in Fluids, 2006, 41(6): 933-947. |
| [13] | Kahler C J, Kompenhans J. Fundamentals of multiple plane stereo particle image velocimetry[J]. Experiments in Fluids, 2000, 29(1): S070-S077. |
| [14] | Wieneke B. Volume self-calibration for 3D particle image velocimetry[J]. Experiments in Fluids, 2008, 45(4): 549-556. |
| [15] | Hartley R, Zisserman A. Multiple view geometry in computer vision[M]. Cambridge University Press, 2003: 153-164. |