CN106153977A - 一种基于单光场相机的三维流场测试方法 - Google Patents

Abstract

本发明公开了一种基于单光场相机的三维流场测试方法，A1，利用单光场相机获取待测流场中示踪粒子的时序粒子光场图像；A2，预处理时序粒子光场图像，去除背景噪声；A3，采用CPU和/或GPU并行处理方法，对时序粒子光场图像进行重构，获得时序三维粒子图像；A4，对时序三维粒子图像进行互相关计算，获得三维流场速度场分布；A5，对三维流场速度场进行后处理，剔除错误速度矢量，同时根据相邻正确矢量插值替代被剔除的错误矢量。本发明能通过单个光场相机的时序图像采集而获取待测流场的三维速度分布，相比于现有多视角三维流场测试方法，极大的减少了硬件系统配置、简化了硬件系统调节步骤，特别适用于受限空间下的三维流场测量。

Description

一种基于单光场相机的三维流场测试方法

技术领域

本发明属于流场观测技术领域，特别涉及一种基于单光场相机的三维流场测试方法。

背景技术

专利申请CN201210260127.8，公开了“一种单CCD相机三维粒子图像测速方法，单CCD相机三维粒子图像测速技术是基于半凸透镜成像原理，在PIV系统中示踪粒子的散射光线通过两上下并排放置的半凸透镜，在像平面或CCD上成两个实像；两个像的距离与物距以及两半凸透镜的光轴距离成函数关系；示踪粒子在景深内无需改变像距便可得到清晰的像；在PIV系统应用中BB1为第一幅照片内同一示踪粒子所成的像之间的距离，CC1为第二幅照片内同一示踪粒子像的距离；采用两半凸透镜做相机的镜头，两半凸透镜间填充不透光的密封介质；示踪粒子运动的二维速度，采用当今通用的互相关算法得到”。但是该方法实际使用中操作复杂，计算精度也不能满足需要。

本发明涉及的参考文献：

[1]Rafael C.Gonzalez,Richard E.Woods.Digital Image Processing(3rdEdition),Pearson Education,2011

[2]Herman GT,Lent A(1976)Iterative reconstructionalgorithms.Computers in Biology and Medicine,6:273–294

[3]Raffel,M.,Willert,C.E.,wereley,s.,Kompenhans,J.,Particle ImageVelocimetryAPractical Guide,Springer-Verlag Berlin Heidelberg,2007

发明内容

本发明的目的是提供一种基于单光场相机的三维流场测试方法。

本发明的技术方案是，一种基于单光场相机的三维流场测试方法，包含以下步骤：

A1，利用单光场相机获取待测流场中示踪粒子的时序粒子光场图像；

A2，预处理时序粒子光场图像，去除背景噪声；

A3，采用CPU和/或GPU并行处理方法，对时序粒子光场图像进行重构，获得时序三维粒子图像；

A4，对时序三维粒子图像进行互相关计算，获得三维流场速度场分布；

A5，对三维流场速度场进行后处理，剔除错误速度矢量，同时根据相邻正确矢量插值替代被剔除的错误矢量。

步骤A1中，在待测流场中撒布示踪粒子，用高亮度光源提供体式照明，然后利用高分辨或者高速图像传感器与微透镜阵列(MLA)的精密封装体形成的单光场相机拍摄示踪粒子的时序光场图像。

步骤A2中，对所拍摄的时序粒子光场图像进行预处理，其中去除背景噪声的降噪运算采用全局阈值、局部阈值、高斯光滑滤波、滑动最小值滤波算法的一种或多种算法组合。

步骤A3中，所述的重构包括：

根据几何光学计算空间光束经主透镜、微透镜阵列后到达图像传感器的准确位置，对空间每一个体素(voxel)所发出的光束，计算其与相交微透镜单元的重叠面积，以及其与相交CCD/CMOS像素的重叠面积，从而计算获得空间每个体素与对应像素的权重系数W_i，j

$Y_{CCD}=\frac{V_B-(s_y+Y_l)}{s_i}{f}_l+s_y=\frac{(V_B-s_y)f_l}{f_m(1-M)}+s_y-\frac{f_l}{f_m(1-M)}{Y}_l---\left(1\right)$

式中各变量的含义如下：

Y_CCD为光线与图像传感器相交点的坐标；

V_B为光线与主镜头相交点的坐标；

s_y为微透镜阵列中某一微透镜单元相对主光轴的偏移值；

Y_l为光线与某一微透镜单元相交点的坐标；

s_i为主透镜与微透镜阵列之间的距离；

f_l为微透镜焦距；

f_m为主透镜焦距；

p_p为像素单元的尺寸；

M为放大系数，

利用密集光线逆追踪方法，定位空间重构体素(Voxel)所对应的像素集合，以每个像素为单元，根据式(1)来计算该像素所采集空间光线的位置，即反向追踪光线，以建立起像素与空间某一点光源所发出光线的对应关系，

利用流场示踪粒子的稀疏性，其对应所需重构的体素区域也为三维稀疏矩阵，将空间某一体素所对应的像素值相乘，如果乘积非零或者大于某一阈值，则表示该体素处可能存在一个示踪粒子，可以对其进行后续的重构计算；反之如果乘积为零或者小于某一阈值，则表示该体素处没有示踪粒子，后续重构计算中可以忽略该体素；

根据乘法代数重构算法，利用计算得到的权重系数W_i，j、以及所筛选出的非零体素，根据所记录的粒子光场像素矩阵(I(x_i，y_i))，迭代计算稀疏重构矩阵的体素值(E(X_j，Y_j，Z_j)^k+1)，

$E{(X_j,Y_j,Z_j)}^{k+1}=E{(X_j,Y_j,Z_j)}^k{\left(\frac{I(x_i,y_i)}{\underset{j\∈N_i}{\Σ}{w}_{i,j}E{(X_j,Y_j,Z_j)}^k}\right)}^{{\mathrm{\μw}}_{i,j}}---\left(2\right)$

式中：

E(X_j,Y_j,Z_j)为空间体素(X_j,Y_j,Z_j)的数值；上标k表示第k次迭代计算所得到该体素的数值；

I(x_i，y_i)为像素(x_i，y_i)处的数值，该数值由光场相机所拍摄的示踪粒子光场图像获得；

w_i,j为体素j与像素i所对应的权重系数；

μ为MART计算参数。

所述步骤A4中，将相邻时刻的示踪粒子三维体素矩阵按照式(3)进行三维互相关计算，获取流场三维速度场分布，

$\mathrm{\Φ}(m,n,l)=\underset{i=1}{\overset{M}{\Σ}}\underset{j=1}{\overset{N}{\Σ}}\underset{k=1}{\overset{L}{\Σ}}{E}_1(i,j,k)\·E_2(i+m,j+n,k+l)---\left(3\right).$

步骤A5中，采用全局阈值计算式(4)、局部中值滤波式(5)或者局部平均值滤波式(6)对获得的三维流场U(i，j，k)，挑选并剔除出其中的错误矢量，并采用线性插值或者三次样条插值的方法，根据正确速度矢量插值获得替代矢量，

$\left\{\begin{array}{c}U(i,j,k)>\stackrel{\&OverBar;}{U(i,j,k)}+T*STD\left(U\right(i,j,k\left)\right)\\ U(i,j,k)<\stackrel{\&OverBar;}{U(i,j,k)}+T*STD\left(U\right(i,j,k\left)\right)\end{array}\right.---\left(4\right)$

$\left\{\begin{array}{c}U(i,j,k)>\stackrel{\&OverBar;}{U(i-1:i+1,j-1:j+1,k-1:k+1)}+T*STD\left(U\right(i-1:i+1,j-1:j+1,k-1:k+1\left)\right)\\ U(i,j,k)<\stackrel{\&OverBar;}{U(i-1:i+1,j-1:j+1,k-1:k+1)}+T*STD\left(U\right(i-1:i+1,j-1:j+1,k-1:k+1\left)\right)\end{array}\right.---\left(5\right)$

其中，是三维速度场的平均值；

STD(U(i，j，k))是三维速度场的标准方差；

T为滤波阈值；

是局部(3×3×3)区域内的三维速度场平均值；

STD(U(i-1：i+1，j-1；j+1，k-1：k+1)是局部(3×3×3)区域内的三维速度场标准方差；

U(i-1：i+1，j-1：j+1，k-1：k+1)是局部(3×3×3)区域内的三维速度场中值。

一种三维流场测试系统，包括由图像传感器与微透镜阵列(MLA)的精密封装体构成的单光场相机，在拍摄示踪粒子的时序光场图像时，在示踪粒子图像的反面具有高亮度体式照明，该体式照明是由包括光源、透镜组合和镜头的光学系统提供的。

本发明能通过单个光场相机的时序图像采集而获取待测流场的三维速度分布，相比于现有多视角三维流场测试方法，极大的减少了硬件系统配置、简化了硬件系统调节步骤，特别适用于受限空间下的三维流场测量。

附图说明

图1本发明的单光场相机三维流场测试系统原理示意图。

图2本发明实施例中光线追踪示意图。

图3本发明实施例中权重系数计算示意图。

图4本发明实施例中密集光线逆追踪方法原理示意图(Dense Ray Tracing，DRT)。

图5是本发明实施例中测试方法的流程图。

具体实施方式

本发明的一种基于单个光场相机的流场三维速度场测试方法，通过单个视角测量获取流场的三维速度场分布。包含以下步骤：

1)在待测流场中撒布示踪粒子(水中一般使用20～30微米的空心玻璃珠；空气中一般使用1微米左右的液滴或者几百纳米的二氧化钛颗粒)，按照图1的方式，用高亮度光源(激光、LED或者高亮投影仪)提供体式照明；然后利用高分辨或者高速图像传感器(CCD或者CMOS)与微透镜阵列(MLA)的精密封装体(即光场相机，以下简称单光场相机)拍摄示踪粒子的时序光场图片。

2)对所拍摄的时序粒子光场图片进行前处理，降噪运算采用全局阈值、局部阈值、高斯光滑滤波、滑动最小值滤波算法的一种或多种算法组合[1]。

3)根据几何光学计算空间光束经主透镜、微透镜阵列后到达图像传感器的准确位置(式1，图2)。按图3所示，对空间每一个体素(voxel)所发出的光束，计算其与相交微透镜单元的重叠面积，以及其与相交CCD/CMOS像素的重叠面积，从而计算获得空间每个体素与对应像素的权重系数W_i，j

$Y_{CCD}=\frac{V_B-(s_y+Y_l)}{s_i}{f}_l+s_y=\frac{(V_B-s_y)f_l}{f_m(1-M)}+s_y-\frac{f_l}{f_m(1-M)}{Y}_l---\left(1\right)$

式中各变量的含义如下(参考图2)：

Y_CCD为光线与图像传感器相交点的坐标；

V_B为光线与主镜头相交点的坐标；

s_y为微透镜阵列中某一微透镜单元相对主光轴的偏移值；

Y_l为光线与某一微透镜单元相交点的坐标；

s_i为主透镜与微透镜阵列之间的距离；

f_l为微透镜焦距；

f_m为主透镜焦距；

p_p为像素单元的尺寸；

M为放大系数。

4)利用密集光线逆追踪方法(Dense RayTracing，DRT)，定位空间重构体素(Voxel)所对应的像素集合。密集光线逆追踪方法采用图2所示的方式，以每个像素为单元，根据式1来计算该像素所采集空间光线的位置。即反向追踪光线，以建立起像素与空间某一点光源所发出光线的对应关系。因为流场示踪粒子的稀疏性，其对应所需重构的体素区域也为三维稀疏矩阵。利用这一特征，将空间某一体素(如图4中红色方框标识的E_j体素)所对应的像素值(如图4中红色线条标识的像素集合)相乘。如果乘积非零或者大于某一阈值，则表示该体素处可能存在一个示踪粒子，可以对其进行后续的重构计算；反之如果乘积为零或者小于某一阈值，则表示该体素处没有示踪粒子，后续重构计算中可以忽略该体素。采用此DRT方法，可以过滤掉稀疏重构矩阵中不含示踪粒子的体素，提高重构计算效率。

5)根据乘法代数重构算法(multiplicativealgebraicreconstructiontechnique,MART，式2，[2])，利用第3步所计算的权重系数W_i，j、以及第4步所筛选出的非零体素，根据所记录的粒子光场像素矩阵(I(x_i，y_i))，迭代计算稀疏重构矩阵的体素值(E(X_j，Y_j，Z_j)^k+1)。一般给定E(X_j，Y_j，Z_j)⁰的初始估计值为1，迭代计算5-10次。

$E{(X_j,Y_j,Z_j)}^{k+1}=E{(X_j,Y_j,Z_j)}^k{\left(\frac{I(x_i,y_i)}{\underset{j\&Element;N_i}{\&Sigma;}{w}_{i,j}E{(X_j,Y_j,Z_j)}^k}\right)}^{{\mathrm{\&mu;w}}_{i,j}}---\left(2\right)$

式中：

E(X_j,Y_j,Z_j)为空间体素(X_j,Y_j,Z_j)的数值；上标k表示第k次迭代计算所得到该体素的数值；

I(x_i，y_i)为像素(x_i，y_i)处的数值，该数值由光场相机所拍摄的示踪粒子光场图像获得；

w_i,j为体素j与像素i所对应的权重系数，其值根据第3)步所描述的方法计算；

μ为MART计算参数[2]。

6)对所记录的示踪粒子时序光场图像按照步骤3～5进行三维重构，将相邻时刻的示踪粒子三维体素矩阵按照式3进行三维互相关计算[3]，获取流场三维速度场分布。

$\mathrm{\&Phi;}(m,n,l)=\underset{i=1}{\overset{M}{\&Sigma;}}\underset{j=1}{\overset{N}{\&Sigma;}}\underset{k=1}{\overset{L}{\&Sigma;}}{E}_1(i,j,k)\&CenterDot;E_2(i+m,j+n,k+l)---\left(3\right)$

7)采用全局阈值(式4)、局部中值滤波(式5)或者局部平均值滤波(式6)[3]对第6步处理获得的三维流场U(i，j，k)，挑选并剔除出其中的错误矢量。并采用线性插值或者三次样条插值的方法，根据正确速度矢量插值获得替代矢量。

$\left\{\begin{array}{c}U(i,j,k)>\stackrel{\&OverBar;}{U(i,j,k)}+T*STD\left(U\right(i,j,k\left)\right)\\ U(i,j,k)<\stackrel{\&OverBar;}{U(i,j,k)}+T*STD\left(U\right(i,j,k\left)\right)\end{array}\right.---\left(4\right)$

$\left\{\begin{array}{c}U(i,j,k)>\stackrel{\&OverBar;}{U(i-1:i+1,j-1:j+1,k-1:k+1)}+T*STD\left(U\right(i-1:i+1,j-1:j+1,k-1:k+1\left)\right)\\ U(i,j,k)<\stackrel{\&OverBar;}{U(i-1:i+1,j-1:j+1,k-1:k+1)}+T*STD\left(U\right(i-1:i+1,j-1:j+1,k-1:k+1\left)\right)\end{array}\right.---\left(5\right)$

上式中是三维速度场的平均值；

STD(U(i，j，k))是三维速度场的标准方差；

T为滤波阈值；

是局部(3×3×3)区域内的三维速度场平均值；

STD(U(i-1：i+1，j-1；j+1，k-1：k+1))是局部(3×3×3)区域内的三维速度场标准方差；

U(i-1：i+1，j-1：j+1，k-1：k+1)是局部(3×3×3)区域内的三维速度场中值。

Claims (7)

1.一种基于单光场相机的三维流场测试方法，包含以下步骤：

A1，利用单光场相机获取待测流场中示踪粒子的时序粒子光场图像；

A2，预处理时序粒子光场图像，去除背景噪声；

A3，采用CPU和/或GPU并行处理方法，对时序粒子光场图像进行重构，获得时序三维粒子图像；

A4，对时序三维粒子图像进行互相关计算，获得三维流场速度场分布；

A5，对三维流场速度场进行后处理，剔除错误速度矢量，同时根据相邻正确矢量插值替代被剔除的错误矢量。

2.如权利要求1所述的基于单光场相机的三维流场测试方法，其特征在于，

步骤A1中，在待测流场中撒布示踪粒子，用高亮度光源提供体式照明，然后利用高分辨或者高速图像传感器与微透镜阵列(MLA)的精密封装体形成的单光场相机拍摄示踪粒子的时序光场图像。

3.如权利要求1所述的基于单光场相机的三维流场测试方法，其特征在于，

步骤A2中，对所拍摄的时序粒子光场图像进行预处理，其中去除背景噪声的降噪运算采用全局阈值、局部阈值、高斯光滑滤波、滑动最小值滤波算法的一种或多种算法组合。

4.如权利要求1所述的基于单光场相机的三维流场测试方法，其特征在于，

步骤A3中，所述的重构包括：

根据几何光学计算空间光束经主透镜、微透镜阵列后到达图像传感器的准确位置，对空间每一个体素(voxel)所发出的光束，计算其与相交微透镜单元的重叠面积，以及其与相交CCD/CMOS像素的重叠面积，从而计算获得空间每个体素与对应像素的权重系数W_i，j

$Y_{CCD}=\frac{V_B-(s_y+Y_l)}{s_i}{f}_l+s_y=\frac{(V_B-s_y)f_l}{f_m(1-M)}+s_y-\frac{f_l}{f_m(1-M)}{Y}_l---\left(1\right)$

式中各变量的含义如下：

Y_CCD为光线与图像传感器相交点的坐标；

V_B为光线与主镜头相交点的坐标；

s_y为微透镜阵列中某一微透镜单元相对主光轴的偏移值；

Y_l为光线与某一微透镜单元相交点的坐标；

s_i为主透镜与微透镜阵列之间的距离；

f_l为微透镜焦距；

f_m为主透镜焦距；

p_p为像素单元的尺寸；

M为放大系数，

利用密集光线逆追踪方法，定位空间重构体素(Voxel)所对应的像素集合，以每个像素为单元，根据式(1)来计算该像素所采集空间光线的位置，即反向追踪光线，以建立起像素与空间某一点光源所发出光线的对应关系，

利用流场示踪粒子的稀疏性，其对应所需重构的体素区域也为三维稀疏矩阵，将空间某一体素所对应的像素值相乘，如果乘积非零或者大于某一阈值，则表示该体素处可能存在一个示踪粒子，可以对其进行后续的重构计算；反之如果乘积为零或者小于某一阈值，则表示该体素处没有示踪粒子，后续重构计算中可以忽略该体素；

根据乘法代数重构算法，利用计算得到的权重系数W_i，j、以及所筛选出的非零体素，根据所记录的粒子光场像素矩阵(I(x_i，y_i))，迭代计算稀疏重构矩阵的体素值(E(X_j，Y_j，Z_j)^k+1)，

$E{(X_j,Y_j,Z_j)}^{k+1}=E{(X_j,Y_j,Z_j)}^k{\left(\frac{I(x_i,y_i)}{\underset{j\&Element;N_i}{\&Sigma;}{w}_{i,j}E{(X_j,Y_j,Z_j)}^k}\right)}^{{\mathrm{\&mu;w}}_{i,j}}---\left(2\right)$

式中：

E(X_j,Y_j,Z_j)为空间体素(X_j,Y_j,Z_j)的数值；上标k表示第k次迭代计算所得到该体素的数值；

I(x_i，y_i)为像素(x_i，y_i)处的数值，该数值由光场相机所拍摄的示踪粒子光场图像获得；

w_i,j为体素j与像素i所对应的权重系数；

μ为MART计算参数。

5.如权利要求4所述的基于单光场相机的三维流场测试方法，其特征在于，

所述步骤A4中，将相邻时刻的示踪粒子三维体素矩阵按照式(3)进行三维互相关计算，获取流场三维速度场分布，

$\mathrm{\&Phi;}(m,n,l)=\underset{i=1}{\overset{M}{\&Sigma;}}\underset{j=1}{\overset{N}{\&Sigma;}}\underset{k=1}{\overset{L}{\&Sigma;}}{E}_1(i,j,k)\&CenterDot;E_2(i+m,j+n,k+l)---\left(3\right).$

6.如权利要求5所述的基于单光场相机的三维流场测试方法，其特征在于，

步骤A5中，

采用全局阈值计算式(4)、局部中值滤波式(5)或者局部平均值滤波式(6)对获得的三维流场U(i，j，k)，挑选并剔除出其中的错误矢量，并采用线性插值或者三次样条插值的方法，根据正确速度矢量插值获得替代矢量，

$\left\{\begin{array}{c}U(i,j,k)>\stackrel{\&OverBar;}{U(i,j,k)}+T*STD\left(U\right(i,j,k\left)\right)\\ U(i,j,k)<\stackrel{\&OverBar;}{U(i,j,k)}+T*STD\left(U\right(i,j,k\left)\right)\end{array}\right.---\left(4\right)$

$\left\{\begin{array}{c}U(i,j,k)>\stackrel{\&OverBar;}{U(i-1:i+1,j-1:j+1,k-1:k+1)}+T*STD\left(U\right(i-1:i+1,j-1:j+1,k-1:k+1\left)\right)\\ U(i,j,k)<\stackrel{\&OverBar;}{U(i-1:i+1,j-1:j+1,k-1:k+1)}+T*STD\left(U\right(i-1:i+1,j-1:j+1,k-1:k+1\left)\right)\end{array}\right.---\left(5\right)$

其中，是三维速度场的平均值；

STD(U(i，j，k))是三维速度场的标准方差；

T为滤波阈值；

是局部(3×3×3)区域内的三维速度场平均值；

STD(U(i-1：i+1，j-1：j+1，k-1：k+1))是局部(3×3×3)区域内的三维速度场标准方差；

U(i-1：i+1，j-1：j+1，k-1：k+1)是局部(3×3×3)区域内的三维速度场中值。

7.一种三维流场测试系统，其特征在于，包括由图像传感器与微透镜阵列(MLA)的精密封装体构成的单光场相机，在拍摄示踪粒子的时序光场图像时，在示踪粒子图像的反面具有高亮度体式照明，该体式照明是由包括光源、透镜组合和镜头的光学系统提供的。