一种基于光场成像的层析PIV测量装置及方法[发明专利]

(12)发明专利申请

(10)申请公布号 CN 106908622 A(43)申请公布日 2017.06.30

(21)申请号 201710153642.9(22)申请日 2017.03.15

(71)申请人 东南大学

地址 211189 江苏省南京市江宁区东南大

学路2号(72)发明人 许传龙　曹丽霞　张彪　宋祥磊　(74)专利代理机构 南京众联专利代理有限公司

32206

代理人 许小莉(51)Int.Cl.

G01P 5/26(2006.01)

权利要求书3页 说明书8页 附图3页

(54)发明名称

一种基于光场成像的层析PIV测量装置及方法

(57)摘要

本发明提供一种基于光场成像的层析PIV测量装置及方法。本发明包括用于产生双脉冲激光并照亮三维流场的光源设备和用于拍摄三维流场的成像设备，所述的光源设备的发光端设置有透镜组，所述的光源设备和所述的成像设备连接时序同步控制器，所述的成像设备连接计算机；所述的透镜组用于将所述的光源设备产生的双脉冲激光发展为一定厚度的体光源和1mm的二维片光源，并照亮流场；所述的成像设备，包括主透镜、CCD探测器以及位于所述的主透镜与CCD探测器之间的微透镜阵列。本发明能快速地计算出权重矩阵，无需标定光场相机的内外参数，装置简单，减小了多相机同步耦合时带来的系统误差。

CN 106908622 ACN 106908622 A

权　利　要　求　书

1/3页

1.一种基于光场成像的层析PIV测量装置，包括用于产生双脉冲激光并照亮三维流场的光源设备和用于拍摄三维流场的成像设备，所述的光源设备的发光端设置有透镜组，所述的光源设备和所述的成像设备连接时序同步控制器，所述的成像设备连接计算机；

所述的透镜组用于将所述的光源设备产生的双脉冲激光发展为一定厚度的体光源和1mm的二维片光源，并照亮流场；

其特征在于：所述的成像设备，包括主透镜、CCD探测器以及位于所述的主透镜与CCD探测器之间的微透镜阵列。

2.根据权利要求1所述的基于光场成像的层析PIV测量装置，其特征在于：所述的微透镜阵列数为N×N，CCD探测器像面关于微透镜的共轭面为主透镜平面，光线经过微透镜在CCD探测器上所成的像称为子图像，每个子图像包含M×M个像素点，其中N和M都大于2。

3.一种用上述基于光场成像的层析PIV测量装置进行测量的方法，其特征在于：该方法包括如下步骤：

步骤一、坐标定义及标定板拍摄：流场及标定板所在的坐标定义为世界坐标O-XYZ，光场原图像中子图像坐标v-om-u为光场原图像中的子图像所在矩阵的行数和列数(v,u)，子图像中像素坐标r-op-s为子图像中像素所在矩阵的行数和列数(r,s)，视差图像中的像素坐标rv-ov-sv为视差图像中像素所在矩阵的行数和列数(rv,sv)，视差图像的图像坐标x-o-y为像素在视差图像中的物理位置(x,y)，原点o位于(rv,sv)＝(1,1)的像素中心位置。实验时，将标定板放置于待测流场所在位置，固定光场相机，使标定板平面平行于光场相机的探测器面，沿着待测流场深度方向等间隔地移动标定板，用光场相机拍摄位于至少三个不同深度位置的标定板，获得光场原图像；

步骤二、视差图像的提取：将步骤一拍摄的标定板光场原图像中的像素进行重新组合，使每个子图像中位于相同位置的像素组成一幅新图像，可以获得M×M幅不同的视差图像，剔除由子图像边缘像素组成的视差图像，从而得到M′×M′幅视差图像，记录标定板角点坐标(X,Y,Z)与对应的每幅视差图像中的像素坐标(rv,sv)，利用下式计算角点的图像坐标(x,y)：

式中，p为微透镜的尺寸；步骤三、三维流场空间坐标与视差图像中的图像坐标关系的计算：对每幅视差图像，利用三阶多项式建立由步骤二得到的标定板角点坐标(X,Y,Z)与图像坐标(x,y)之间的关系：

式中，已知(X,Y,Z)和(x,y)利用MATLAB拟合工具箱计算出a1,a2,…,a20和b1,b2,…,b20，对于每一幅视差图像需要拍摄至少三幅不同深度位置的标定板图像，每一幅视差图像在不同深度位置计算的a1,a2,…,a20和b1,b2,…,b20的值不同，然后利用线性插值计算出任意深度位置的a1,a2,…,a20和b1,b2,…,b20的值，从而确定三维流场空间坐标与每幅视差图像中的图像坐标关系；

2

CN 106908622 A

权　利　要　求　书

2/3页

步骤四、权重矩阵的计算：将流场沿X轴Y轴Z轴划分成m,n,l个体素，将每个体素的中心坐标(Xc,Yc,Zc)分别代入步骤三中的式(3)、式(4)得到相应的图像坐标，再将图像坐标代入步骤二中的式(2)得到每个体素对应的像素坐标，利用高斯函数计算每个体素对每幅视差图像中的像素的贡献W(权重矩阵)，保存权重矩阵；

步骤五、流场的光场原图像的获取与预处理：将标定板移出待测流场，打开双脉冲激光器，使激光经过透镜组产生具有一定厚度的体光源，当激光入射到含有示踪粒子的待测流动对象时，示踪粒子会产生散射光，由光场相机连续两次拍摄三维流场，获得光场原图像，根据步骤二的方法获得一系列三维流场的视差图像的灰度值Iv；

步骤六、层析重建计算三维粒子场：根据步骤四计算得到的权重矩阵W及步骤五得到的视差图像的灰度值Iv，利用MART算法反演三维粒子灰度分布E：

式中，c为体素的序号，c＝1,2,3,…,m×n×l；Xc,Yc,Zc为体素中心坐标；k为迭代次数；g为像素的序号，g＝1,2,3,…,N×N×M′×M′；μ为松弛因子；

步骤七、三维互相关算法计算流速：对步骤六重建所得的三维粒子灰度分布E沿X轴Y轴Z轴分别以xsize,ysize,zsize的尺寸进行网格划分(查询窗口)，然后利用下式计算三维流场速度矢量分布：

式中，R为互相关系数；xsize,ysize,zsize表示查询窗口分别在X,Y,Z方向上的尺寸；Et，Et+Δt分别为t时刻和t+Δt时刻的三维粒子灰度分布；σσt,Δ+t分别为t时刻和t+Δt时刻查询窗口内的体素灰度值的标准差；m′,n′,l′为查询窗口内的体素分别在X,Y,Z方向上的坐标索引(体素的序号)，式(8)的最大值对应的Δm′,Δn′,Δl′为粒子分别在X,Y,Z方向上的位移，已知Δt可以获得流场的三维速度矢量分布。

4.根据权利要求3所述的基于光场成像的层析PIV测量装置进行测量的方法，其特征在于：步骤二中所述M×M幅不同视差图像的提取方法如公式(1)确定：

式中，I(N,N,r,s)为第N行第N列微透镜子图像中的第r行第s列的像素的灰度值；r为每个子图像中的像素所在的行数，r＝1,2,3,…,M；s为每个子图像中的像素所在的列数，s＝1,2,3,…,M。视差图像里相邻两个像素之间的距离为单个微透镜的尺寸p。

5.根据权利要求3或4所述的基于光场成像的层析PIV测量装置进行测量的方法，其特征在于：步骤四中所述利用高斯函数来计算每个体素对每幅视差图像中的像素的贡献(权重矩阵)，其公式为：

式中，“T”为对矩阵转置；h为每个粒子图像灰度峰值；d,e,f表征粒子图像的形状；Wv

3

CN 106908622 A

权　利　要　求　书

3/3页

(rv,sv)为视差图像中第rv行第sv列像素的灰度值，rv＝1,2,3,…,N，sv＝1,2,3,…,N；rv2，sv2为每个粒子图像灰度峰值所在的像素坐标。

4

CN 106908622 A

说　明　书

一种基于光场成像的层析PIV测量装置及方法

1/8页

技术领域：[0001]本发明涉及一种基于光场成像的层析PIV测量装置及方法，属于流体测量与可视化技术领域。

背景技术：[0002]流速是描述流场特性的重要参数之一。在科研领域和工业领域都涉及到流场流速的测量，如锅炉燃烧过程中燃料传输速度的控制对于燃烧效率、安全稳定的燃烧具有重要的意义；微流体的传热效率与流场流速有关等。实现非接触、瞬态、三维全场流速的测量是实验流体力学中的热点之一。[0003]随着计算机技术、传感器技术、光电技术的快速发展，PIV技术已经广泛应用于流场流速的测量。该技术利用CCD在短时间内连续拍摄两张流场图像，利用图像处理技术来获取速度矢量分布，该过程无需接触流场，易获得流场中多点流速。目前现有的PIV技术主要归结为三类：2D-2C PIV技术、2D-3C PIV技术(准三维PIV技术)、3D-3C PIV技术。2D-2C PIV技术是流速测量中最基础的方法之一，该方法利用约为1mm厚度的激光片光源照明流场，根据互相关算法计算两幅图像的二维速度。2D-3C PIV技术利用单个或多个激光片光源照明流场。首先根据二维PIV技术计算流场切面的二维速度分量，运用三维流动的连续性方程或其它技术计算出切面法线方向的第三个速度分量。该方法是基于2维PIV技术的方法，又称为准三维PIV技术，但该方法在法线方向上的速度分辨率较低，难以实现体的测量。不同于以上两种方法，3D-3C PIV是利用具有一定厚度的激光体光源照明某个容积内的流场，利用数学方法重构出示踪粒子的三维位置坐标和强度信息(三维粒子场)，利用三维互相关技术计算速度矢量分布。该方法能同时获得流场的三个速度分量。3D-3C PIV技术能获得真正意义上的三维速度矢量分布。层析PIV是3D-3C PIV中的一种，由于该方法具有较高的空间分辨率，适用范围较广，已受广大学者青睐。但层析PIV技术往往需要3～6台相机拍摄流场，大大增加了层析PIV技术的成本、复杂程度及重建计算量。每台CCD需要严格对准流场，目前的层析PIV技术的精度非常依赖于CCD的布局及标定的准确性。

发明内容[0004]本发明的目的是针对基于传统相机的层析PIV技术存在的缺点，提供一种基于光场成像的层析PIV测量装置及方法，利用光场成像能同时记录光场的位置和方向信息，使得单台光场相机代替多台传统相机用于层析PIV测量成为可能，且能克服传统层析PIV的缺点。[0005]上述的目的通过以下技术方案实现：[0006]一种基于光场成像的层析PIV测量装置，包括用于产生双脉冲激光并照亮三维流场的光源设备和用于拍摄三维流场的成像设备，所述的光源设备的发光端设置有透镜组，所述的光源设备和所述的成像设备连接时序同步控制器，所述的成像设备连接计算机；[0007]所述的透镜组用于将所述的光源设备产生的双脉冲激光发展为一定厚度的体光

5

CN 106908622 A

说　明　书

2/8页

源和1mm的二维片光源，并照亮流场；[0008]所述的成像设备，包括主透镜、CCD探测器以及位于所述的主透镜与CCD探测器之间的微透镜阵列。[0009]所述的基于光场成像的层析PIV测量装置，所述的微透镜阵列数为N×N，CCD探测器像面关于微透镜的共轭面为主透镜平面，光线经过微透镜在CCD探测器上所成的像称为子图像，每个子图像包含M×M个像素点，其中N和M都大于2。[0010]用上述基于光场成像的层析PIV测量装置进行测量的方法，该方法包括如下步骤：[0011]步骤一、坐标定义及标定板拍摄：流场及标定板所在的坐标定义为世界坐标O-XYZ，光场原图像中子图像坐标v-om-u为光场原图像中的子图像所在矩阵的行数和列数(v,u)，子图像中像素坐标r-op-s为子图像中像素所在矩阵的行数和列数(r,s)，视差图像中的像素坐标rv-ov-sv为视差图像中像素所在矩阵的行数和列数(rv,sv)，视差图像的图像坐标x-o-y为像素在视差图像中的物理位置(x,y)，原点o位于(rv,sv)＝(1,1)的像素中心位置。实验时，将标定板放置于待测流场所在位置，固定光场相机，使标定板平面平行于光场相机的探测器面，沿着待测流场深度方向等间隔地移动标定板，用光场相机拍摄位于至少三个不同深度位置的标定板，获得光场原图像；[0012]步骤二、视差图像的提取：将步骤一拍摄的标定板光场原图像中的像素进行重新组合，使每个子图像中位于相同位置的像素组成一幅新图像，可以获得M×M幅不同的视差图像，剔除由子图像边缘像素组成的视差图像，从而得到M′×M′幅视差图像。记录标定板角点坐标(X,Y,Z)与对应的每幅视差图像中的像素坐标(rv,sv)，利用下式计算角点的图像坐标(x,y)：

[0013]

式中，p为微透镜的尺寸。

[0015]步骤三、三维流场空间坐标与视差图像中的图像坐标关系的计算：对每幅视差图像，利用三阶多项式建立由步骤二得到的标定板角点坐标(X,Y,Z)与图像坐标(x,y)之间的关系：

[0016][0017]

[0014]

式中，已知(X,Y,Z)和(x,y)利用MATLAB拟合工具箱计算出a1,a2,…,a20和b1,

b2,…,b20。对于每一幅视差图像需要拍摄至少三幅不同深度位置的标定板图像，每一幅视差图像在不同深度位置计算的a1,a2,…,a20和b1,b2,…,b20的值不同，然后利用线性插值计算出任意深度位置的a1,a2,…,a20和b1,b2,…,b20的值，从而确定三维流场空间坐标与每幅视差图像中的图像坐标关系；[0019]步骤四、权重矩阵的计算：将流场沿X轴Y轴Z轴划分成m,n,l个体素，将每个体素的中心坐标(Xc,Yc,Zc)分别代入步骤三中的式(3)、式(4)得到相应的图像坐标，再将图像坐标代入步骤二中的式(2)得到每个体素对应的像素坐标，利用高斯函数计算每个体素对每幅视差图像中的像素的贡献W(权重矩阵)，保存权重矩阵；

6

[0018]

CN 106908622 A[0020]

说　明　书

3/8页

步骤五、流场的光场原图像的获取与预处理：将标定板移出待测流场，打开双脉冲

激光器，使激光经过透镜组产生具有一定厚度的体光源，当激光入射到含有示踪粒子的待测流动对象时，示踪粒子会产生散射光，由光场相机连续两次拍摄三维流场，获得光场原图像，根据步骤二的方法获得一系列三维流场的视差图像的灰度值Iv；[0021]步骤六、层析重建计算三维粒子场：根据步骤四计算得到的权重矩阵W及步骤五得到的视差图像的灰度值Iv，利用MART算法反演三维粒子灰度分布E：

[0022]

式中，c为体素的序号，c＝1,2,3,…,m×n×l；Xc,Yc,Zc为体素中心坐标；k为迭代次数；g为像素的序号，g＝1,2,3,…,N×N×M′×M′；μ为松弛因子。[0024]步骤七、三维互相关算法计算流速：对步骤六重建所得的三维粒子灰度分布E沿X轴Y轴Z轴分别以xsize,ysize,zsize的尺寸进行网格划分(查询窗口)，然后利用下式计算三维流场速度矢量分布：

[0025]

[0023]

式中，R为互相关系数；xsize,ysize,zsize表示查询窗口分别在X,Y,Z方向上的尺寸；

Et，Et+Δt分别为t时刻和t+Δt时刻的三维粒子灰度分布；σσt,Δ+t分别为t时刻和t+Δt时刻查询窗口内的体素灰度值的标准差；m′,n′,l′为查询窗口内的体素分别在X,Y,Z方向上的坐标索引(体素的序号)。式(8)的最大值对应的Δm′,Δn′,Δl′为粒子分别在X,Y,Z方向上的位移，已知Δt可以获得流场的三维速度矢量分布。[0027]所述的基于光场成像的层析PIV测量装置进行测量的方法，步骤二中所述M×M幅不同视差图像的提取方法如公式(1)确定：

[0026]

[0028]

式中，I(N,N,r,s)为第N行第N列微透镜子图像中的第r行第s列的像素的灰度值；r

为每个子图像中的像素所在的行数，r＝1,2,3,…,M；s为每个子图像中的像素所在的列数，s＝1,2,3,…,M。视差图像里相邻两个像素之间的距离为单个微透镜的尺寸p。[0030]所述的基于光场成像的层析PIV测量装置进行测量的方法，步骤四中所述利用高斯函数来计算每个体素对每幅视差图像中的像素的贡献(权重矩阵)，其公式为：

[0031]

[0029]

式中，“T”为对矩阵转置；h为每个粒子图像灰度峰值；d,e,f表征粒子图像的形状；Wv(rv,sv)为视差图像中第rv行第sv列像素的灰度值，rv＝1,2,3,…,N，sv＝1,2,3,…,N；rv2，sv2为每个粒子图像灰度峰值所在的像素坐标。[0033]有益效果：[0034]1.本发明采用单光场相机拍摄流场中示踪粒子发出的散射光。相比于传统相机，光场相机在主透镜和CCD探测器之间加了一微透镜阵列。流场本身不会发光，所以需要往流

7

[0032]

CN 106908622 A

说　明　书

4/8页

场中加入示踪粒子，示踪粒子被入射光照射时会产生散射光。散射光经过镜头投影到微透镜阵列上，再经过微透镜阵列二次成像投影到CCD探测器上，形成一系列子图像，每个子图像记录了光场的位置信息，每个子图像中位于相同位置的像素记录了光场的某一视角信息，本发明成像装置能同时实现光场方向信息、位置信息及散射光强度的采集。[0035]2.相比于传统层析PIV装置，由于光场相机能同时记录光场的位置和方向信息，采用单光场相机代替多台传统相机拍摄流场，从而实现单相机多方向的采集，系统更紧凑，标定更方便，能减小多相机同步耦合时带来的系统误差。[0036]3.在测量方法上，光场相机中微透镜阵列数为N×N，每个子图像包含M×M个像素点，N和M都大于2。光场相机可以被看成以微透镜阵列为镜头的许多微型相机对物体在不同视角成像的装置。由于相同方向的光线被微透镜下相同位置的像素接收，本发明将每个子图像中相同位置的像素组成一幅新图像，从而将光场图像转换为一系列不同的视差图像，不同的视差图像代表对示踪粒子不同方向的成像，能提取出M×M幅不同的视差图像。通过标定进一步建立三维流场空间坐标与每幅视差图像中像素的图像坐标之间的映射关系求解权重矩阵，建立层析重建模型，利用MART算法反演三维粒子灰度分布，利用三维互相关实现3D-3C速度场测量。[0037]4.对于光场相机的权重矩阵计算，本发明对标定板图像中每个子图像相同位置的像素进行重新组合提取出一系列视差图像，利用三阶多项式建立三维流场空间与这些视差图像中的像素之间的对应关系，然后利用高斯函数来计算三维体素对各个视差图像中的像素的贡献(权重矩阵)，该方法计算速度快。[0038]5.本发明通过建立标定板角点坐标与视差图像中像素的图像坐标之间的映射关系来确定三维流场的每个体素坐标与像素之间的关系，该方法不需要标定光场相机的内外参数。

附图说明[0039]图1为本发明所述的基于光场成像的层析PIV测量装置结构示意图；[0040]图中的标记含义：1、双脉冲激光器，2、透镜组；3、时序同步控制器；4、光场相机；5、计算机；6、待测流场；7、示踪粒子；8、传输线。[0041]图2为光场相机结构及成像原理图；图中的标记含义：9、虚拟物面；10、光线；11、主透镜；12、微透镜阵列；13、CCD探测器；14、像素；15、子图像。[0042]图3(a)为光场原图像坐标定义示意图；图中的标记含义：16、光场原图像；横坐标omu轴为子图像所在矩阵的列数，纵坐标omv轴为子图像所在矩阵的行数；图3(b)为光场原图像中子图像像素坐标定义示意图；图中的标记含义：17、子图像边缘像素；横坐标ops轴为子图像中的像素所在矩阵的列数，纵坐标opr轴为子图像中的像素所在矩阵的行数；图3(c)为视差图像坐标定义示意图；图中的标记含义：18、视差图像，两条横坐标ovsv轴,oy轴分别为视差图像中像素所在矩阵的列数及像素在横向的物理位置，两条纵坐标ovrv轴,ox轴分别为视差图像中像素所在矩阵的行数及像素在纵向的物理位置。[0043]图4(a)为三维体标定示意图；图中的标记含义：19、标定板；图4(b)为标定板的示意图；图中的标记含义：20、标定板角点。[0044]图5为体素和查询窗口划分示意图；图中的标记含义：21、体素；22查询窗口。

8

CN 106908622 A

说　明　书

5/8页

具体实施方式[0045]下面结合具体实施方式，进一步阐明本发明，应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。[0046]本发明的基于光场成像的层析PIV装置工作原理为：由于流场本身不会发光，所以在待测流场中加入跟随性、散射性较好的示踪粒子，将具有一定厚度的激光入射到流场中，使示踪粒子产生散射光。将光场相机设置在双帧工作模式下，利用时序同步控制器调节激光器的双脉冲信号的时间间隔，使两个激光脉冲分别落在光场相机的双曝光时间上，从而使光场相机连续两次拍摄到三维流场图像。然后利用计算机对采集到的图像进行视差图像提取、层析重建、三维互相关计算，从而得到3D-3C速度场分布。[0047]本发明装置能够实现单相机的3D-3C速度场分布测量，且无需标定光场相机内外参数。测量装置如图1所示，主要包括双脉冲激光器1、透镜组2、时序同步控制器3、光场相机4、计算机5、待测流场6、示踪粒子7、传输线8。双脉冲激光器1是整个测量装置的光源(波长为532nm，光束直径D～7mm)。透镜组2的作用是将双脉冲激光发展为具有一定体积(如100mm×100mm×100mm)的激光光束或厚度为1mm的二维片激光光束(如1mm×100mm×100mm)。时序同步控制器3产生多个信号通过传输线8分别传输给双脉冲激光器1和光场相机4。光场相机4与双脉冲激光器1成正交放置。计算机5的作用是图像存储、视差图像提取、层析重建、三维互相关计算。待测流场6放置于脉冲激光器1和光场相机4正交的交点处，使激光能垂直入射到待测流场6中。示踪粒子7被均匀地撒入待测流场中。[0048]不同于传统的层析PIV装置，本发明采用单光场相机4代替多台传统相机来拍摄流场。光场相机结构如图2所示，光场相机4在主透镜11和CCD探测器13之间加了一微透镜阵列12，使经过镜头的光线10投影到微透镜阵列12上，光线经过微透镜阵列12二次投影到CCD探测器13上，产生一幅光场原图像。光场原图像里包含一系列子图像15，每个子图像15记录了光场的位置信息，子图像15内每个像素点14记录了光场的方向信息，能同时实现光场方向和位置信息的采集。这一特性使单光场相机4重建流场中的示踪粒子7的三维位置和强度信息成为可能。[0049]基于光场成像的层析PIV的测量方法，其步骤包括如下：[0050]步骤一、坐标定义及标定板拍摄：流场6及标定板19所在的坐标定义为世界坐标O-XYZ。光场原图像16中子图像坐标v-om-u为光场原图像16中的子图像15所在矩阵的行数和列数(v,u)，子图像15中像素坐标r-op-s为子图像15中像素14所在矩阵的行数和列数(r,s)，视差图像18中的像素坐标rv-ov-sv为视差图像18中像素14所在矩阵的行数和列数(rv,sv)，视差图像18的图像坐标x-o-y为像素14在视差图像18中的物理位置(x,y)，原点o位于(rv,sv)＝(1,1)的像素14中心位置，如图3所示。实验时，将标定板19放置于待测流场6所在位置，固定光场相机4，使标定板19平面平行于光场相机4的探测器面13，沿着待测流场深度方向等间隔地移动标定板19，用光场相机4拍摄位于至少三个不同深度位置的标定板，获得光场原图像16，该过程如图4所示；[0051]步骤二、视差图像的提取：将步骤一拍摄的标定板光场原图像16中的像素14按式(1)进行重新组合：

9

CN 106908622 A

说　明　书

6/8页

[0052]

式中，I(N,N,r,s)为第N行第N列微透镜子图像中的第r行第s列的像素的灰度值；r

为每个子图像中的像素所在的行数，r＝1,2,3,…,M；s为每个子图像中的像素所在的列数，s＝1,2,3,…,M。视差图像里相邻两个像素之间的距离为单个微透镜的尺寸p。[0054]式(1)使每个子图像15中位于相同位置的像素组成一幅新图像，可以获得M×M幅不同的视差图像18V，剔除由子图像边缘像素17组成的视差图像，从而获得M′×M′幅视差图像。视差图像里的像素坐标(rv,sv)与图像坐标(x,y)的关系按公式(2)计算：

[0055]

[0053]

记录标定板角点20坐标(X,Y,Z)与对应的每幅视差图像18中的图像坐标(x,y)。[0057]步骤三、三维流场空间坐标与视差图像中的图像坐标关系的计算：对每幅视差图像18，利用三阶多项式建立由步骤二得到的标定板角点20坐标(X,Y,Z)与图像坐标(x,y)之间的关系：

[0058][0059]

[0056]

式中，已知(X,Y,Z)和(x,y)利用MATLAB拟合工具箱计算出a1,a2,…,a20和b1,

b2,…,b20。每幅视差图像18在不同深度位置计算的a1,a2,…,a20和b1,b2,…,b20的值不同，已知至少三个不同深度位置的a1,a2,…,a20和b1,b2,…,b20的值，利用线性插值计算出任意深度位置的a1,a2,…,a20和b1,b2,…,b20的值。从而确定三维流场空间坐标与每幅视差图像18中的图像坐标的函数关系。[0061]步骤四、权重矩阵的计算：[0062](a)高斯函数系数的标定：移除标定板19，调整透镜组2使激光厚度为1mm(称为片光源)，将片光源入射到含有示踪粒子7的待测流场区域，且片光源平面平行于光场相机4的探测器面13。当激光入射到含有示踪粒子7的待测区域时，示踪粒子会产生散射光，由光场相机4拍摄至少三幅不同深度位置的示踪粒子的光场原图像16，且使每个平面都有约为50个互不重叠的粒子被光场相机记录，下面以虚拟物面这个深度图像为例。根据步骤二中的式(1)获得位于虚拟物面处粒子的M′×M′幅视差图像18。对每一幅视差图像18进行如下处理：识别这50个粒子图像的灰度峰值p及峰值所在的像素坐标14(rv2,sv2)，即利用高斯函数分别对50个粒子进行拟合，高斯函数公式为：

[0063]

[0060]

式中，“T”为对矩阵转置；h为每个粒子图像灰度峰值；d,e,f表征粒子图像的形状；Wv(rv,sv)为视差图像中第rv行第sv列像素的灰度值，rv＝1,2,3,…,N，sv＝1,2,3,…,N；rv2，sv2为每个粒子图像灰度峰值所在的像素坐标。

[0064]

10

CN 106908622 A[0065]

说　明　书

7/8页

已知h,rv,sv,rv2,sv2,Wv(rv,sv)的值，可以拟合出50个粒子对应的系数d,e,f，然后

求平均值可以确定公式(5)中的系数d,e,f。其它深度位置的系数d,e,f按该方法计算，不同深度位置计算所得的d,e,f的值不同，可以利用线性插值计算任意深度位置的d,e,f的值。从而确定了位于任意深度位置每幅视差图像18对应的高斯函数系数。[0066](b)数据处理：将流场沿X轴Y轴Z轴划分成m,n,l个体素21(如m＝n＝1001，l＝101，每个体素的边长为0.1mm)，将每个体素21的中心坐标(Xc,Yc,Zc)(Xc,Yc,Zc分别表示第c个体素的中心在世界坐标系上的坐标值，c＝1,2,3,…,m×n×l)的值代入步骤三的式(3)、式(4)中的变量X,Y,Z计算每个体素21的中心坐标对应的图像坐标，然后利用式(2)将图像坐标转化为像素坐标，从而确定了每个体素21对应的视差图像18中的像素坐标。将得到的像素坐标代入到公式(5)中的变量rv2,sv2，从而得到每个体素21对应每幅视差图像18中的所有像素14的灰度值Wv，并转换成一系列一维列向量V′(r,s,c)＝[Wv(1,1),Wv(1,2),Wv(1,3),...,Wv(N,N)]T，保存形式如下：

[0067]

式中，W为空间体素21对各幅视差图像18中的所有像素14的贡献；V′(M′,M′,m×n×l)为第m×n×l个体素21对由每个子图像中像素坐标为(M′,M′)的像素所组成的视差图像18的贡献。[0069]步骤五、流场的光场原图像的获取与预处理：调整透镜组2，使激光经过透镜组2产生具有一定厚度的体光源入射到待测流场中。由光场相机4连续两次拍摄三维流场，获得光场原图像16，根据步骤二中的式(1)获得M′×M′幅三维流场的视差图像18V(r,s)。[0070]步骤六、层析重建计算三维粒子场：[0071](a)数据处理：将步骤五得到的三维流场的每幅视差图像18的灰度值V(r,s)转化为一维列向量的形式V(r,s)＝[I(1,1,r,s),I(1,2,r,s),I(1,3,r,s),…,I(N,N,r,s)]T，然后将M′×M′幅不同视差图像灰度值列向量组成一个新的列向量Iv(xg,yg)＝[V(1,1),V(1,2),V(1,3),…,V(M′,M′)]T。[0072](b)层析重建：初始化三维粒子灰度分布E1＝[1,2,3,…,m×n×l]′；已知权重矩阵W和视差图像18灰度值Iv的值，利用MART算法反演三维粒子灰度分布E：

[0073]

[0068]

式中，c为体素的序号，c＝1,2,3,…,m×n×l；Xc,Yc,Zc为体素中心坐标；k为迭代次数；g为像素的序号，g＝1,2,3,…,N×N×M′×M′；μ为松弛因子。[0075]步骤七、三维互相关算法计算流速：对步骤六重建所得的三维粒子灰度分布E沿X轴Y轴Z轴分别以xsize,ysize,zsize(如xsize＝ysize＝zsize＝16个体素)的尺寸进行网格划分(查询窗口)22。然后利用下式计算三维流场速度矢量分布：

[0076][0077]

[0074]

式中，R为互相关系数；xsize,ysize,zsize表示查询窗口分别在X,Y,Z方向上的尺寸；

11

CN 106908622 A

说　明　书

8/8页

Et，Et+Δt分别为t时刻和t+Δt时刻的三维粒子灰度分布；σσt,Δ+t分别为t时刻和t+Δt时刻查询窗口内的体素灰度值的标准差；m′,n′,l′为查询窗口内的体素分别在X,Y,Z方向上的坐标索引(体素的序号)。式(8)的最大值对应的Δm′,Δn′,Δl′为粒子分别在X,Y,Z方向上的位移，已知Δt可以获得流场的三维速度矢量分布。

12

CN 106908622 A

说　明　书　附　图

1/3页

图1

图2

13

CN 106908622 A

说　明　书　附　图

2/3页

图3

14

CN 106908622 A

说　明　书　附　图

3/3页

图4

图5

15