CN104680564A - 一种灰度增强层析piv重构方法、装置和设备 - Google Patents

Abstract

本发明公开了一种灰度增强层析粒子图像测速(PIV)重构方法，该方法包括：在完成倍增代数重构技术(MART)迭代后，统计空间灰度场中粒子拉长比；依据所述粒子拉长比确定逆扩散强度因子，并依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场；依据统计得到的粒子浓度计算灰度抑制因子，并依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配，之后进行下一次MART迭代。本发明还同时公开了一种实现所述方法的装置和设备。

Description

一种灰度增强层析PIV重构方法、装置和设备

技术领域

本发明涉及基于灰度增强的层析粒子图像测速技术，尤其涉及一种灰度增强层析粒子图像测速(Particle Image Velocimetry，简称PIV)重构方法、装置和设备。

背景技术

PIV是一种现代激光测速技术，主要运用于流场速度测量，通过追踪示踪粒子在流场中的运动来得到速度场。最近兴起的层析粒子图像测速技术(层析PIV)成功地将二维PIV推广到三维流场测量，能够获得瞬时的三维三分量(3D3C)的速度场。该技术通过不同视角下(一般为4个相机)粒子散射成像重构出空间粒子的真实分布，然后采用三维互相关计算相邻曝光粒子间的位移。空间粒子重构是该方法的关键。

层析PIV的相机布置如图1所示，其中相机为‘┼’字型布置。空间测量体E中的示踪粒子被激光照亮以后按照投影关系同时成像在四个不同视角的相机之上。投影成像的粒子灰度一般呈高斯分布，大小与示踪粒子大小和散射强度有关，一般占3×3像素。从测量空间到相机平面的投影关系可以事先标定得到，映射函数决定了粒子之间的相对位置。如果将空间测量体E离散成和像素(pixel)大小相当的体素(voxel)，投影成像可以简化成一系列线性方程组：

WE＝I   (1)

即假设粒子图像是空间粒子灰度沿视线投影积分的结果。其中W为投影权重函数，W_ij表示第j个体素对第i个像素的贡献，该贡献与投影距离有关，距离越小贡献越大。空间灰度重构可以认为是一个逆投影过程，就是已知粒子图像I和权重函数W求空间灰度分布E。由于离散的体素个数远远大于已知的像素个 数，所以该逆问题存在不定解。为了求解该方程组，需要额外添加约束条件。目前，主流的层析PIV都采用基于熵最大的倍增代数重构技术(multiplicative algebraic reconstruction technique，MART)对粒子场进行重构还原。该方法被证明是当前最适合层析PIV重构的，该方法迭代收敛速度快，重构的精度较高，其迭代公式可简单表述为：

E^k+1＝E^k(I/WE^k)^uW   (2) 

其中，k代表迭代次数，I/WE^k代表对投影误差的评估，空间灰度E会根据误差的大小不断逼近最优解。

然而，在实际的使用中发现，重构精度会受到相机个数与视角，粒子浓度和映射函数精度等因素的影响，出现粒子拉长和虚假粒子两个非常棘手的问题。所谓粒子拉长，是指层析PIV在理想情况下重构出来应该是高斯分布的球形粒子，正如平面PIV中的粒子服从高斯分布一样。但由于相机视角以及个数原因，投影灰度并不能完全反应粒子的空间形状，导致了在相机轴线方向(测量体厚度方向)粒子被拉长。

此外，虚假粒子也是层析PIV重构面临的一个难点。图像采集实际上是一个投影积分过程，空间连续的粒子灰度被映射到离散的平面图像之上。MART算法就是希望通过这样的欠采样图像重构出空间粒子灰度分布，这是反投影过程。由于可利用的粒子位置信息有限，导致在所有相机视线的交点处都可能出现粒子，但是这个粒子是否真实存在就不得而知。

通过上面的分析可知，层析PIV方法中的MART算法并没有利用粒子形状、灰度等许多有用的信息，导致了其自身更新速度很慢，甚至不正确。

发明内容

为解决现有存在的技术问题，本发明实施例提供一种灰度增强层析PIV重构方法、装置和设备。

本发明实施例提供了一种灰度增强层析粒子图像测速(PIV)重构方法，该方法包括：

A、在完成倍增代数重构技术MART迭代后，统计空间灰度场中的粒子拉长比；

B、依据所述粒子拉长比确定逆扩散强度因子，并依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场；

C、依据统计得到的粒子浓度计算灰度抑制因子，并依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配；

之后进行下一次MART迭代。

其中，所述统计空间灰度场中粒子拉长比，包括：

在空间灰度场中采用局部峰值的方法识别粒子，并对这些粒子的灰度做统计平均得到所有粒子灰度的平均结果；对该结果进行三维高斯拟合，得到所述空间灰度场各方向的高斯拟合的标准差，依据所述标准差计算所述粒子拉长比。

其中，所述三维高斯拟合按照如下公式进行：

$I(x,y,z)=I_0\exp [-\frac{{(x-x_0)}^2}{(1/8)d_{\mathrm{\τx}}^2}-\frac{{(y-y_0)}^2}{(1/8)d_{\mathrm{\τy}}^2}-\frac{{(z-z_0)}^2}{(1/8)d_{\mathrm{\τz}}^2}];$

其中，所述x,y,z代表重构体的三维空间坐标，z为重构体厚度方向，x₀,y₀,z₀代表粒子中心的物理位置，I代表单个粒子的灰度，I₀代表粒子中心的灰度值，d_τx,d_τy,d_τz为x,y,z方向粒子拉长的初始直径。

其中，所述依据空间灰度场各方向的高斯拟合的标准差计算所述粒子拉长比，依据如下公式进行：

$r_d=\frac{d_{\mathrm{\τz}}}{d_{\mathrm{\τ}}}=\frac{{\mathrm{\σ}}_{\mathrm{\τz}}}{{\mathrm{\σ}}_{\mathrm{\τ}}};$

其中，所述r_d为粒子拉长比，所述d_τ为粒子图像直径，所述σ_τ为粒子高斯拟合后在x,y方向标准差的平均值，σ_τ＝1/8(d_τx+d_τy)，所述σ_τz为粒子高斯拟合后在z方向的标准差，σ_τz＝1/4d_τz。

其中，所述依据所述粒子拉长比确定逆扩散强度因子，依据如下公式进行：

$\mathrm{\δ}=p_1\·{r_d}^{p_2}+p_3;$

其中，所述δ表示逆扩散强度因子，所述r_d为粒子拉长比，所述p₁＝-1.40，p₂＝-1.97，p₃＝1.41。

其中，所述依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场，依据如下公式进行：

${\tilde{E}}^k=E^k-\mathrm{\δ}\·E_{\mathrm{zz}}^k;$

该公式为所述逆扩散方程，所述表示灰度关于z方向的二阶导数，E^k表示所述MART迭代后的灰度，表示在z方向逆扩散后的灰度，δ为所述逆扩散强度因子，k为MART迭代次数。

本发明实施例中，所述依据统计得到的粒子浓度计算灰度抑制因子，包括：

依据公式(a)和统计得到的粒子浓度计算得到A_g和B_g；利用灰度矩阵的最大值归一化所述粒子灰度场，并结合公式(b)计算得到虚假粒子的概率分布，并依据所述虚假粒子的概率分布计算得到所述灰度抑制因子，所述灰度抑制因子α与所述虚假粒子的概率分布g之和为1；其中，

$\begin{array}{c}{A}_g=a_1\·{\mathrm{ppp}}^{a_2}+a_3\\ B_g=b_1\·\mathrm{ppp}+b_2\end{array};---\left(a\right)$

其中，a₁＝4.17×10^-5，a₂＝-3，a₃＝6.42×10^-2，b₁＝1.82×10^-1，b₂＝2.65×10^-3，ppp为统计得到的粒子浓度；

$\begin{array}{c}{g}_{(E_r,\mathrm{ppp})}=A_g\·\exp (-E_r/B_g)\\ E_r=E/\max \left(E\right)\end{array};---\left(b\right)$

其中，E为所述灰度矩阵，E_r是用最大灰度归一化后的相对灰度，所述g为虚假粒子的概率分布。

其中，所述依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配方法，依据如下公式进行：

${\tilde{E}}_{m+1}^k={\mathrm{\α}}_{({\tilde{E}}_m^k,\mathrm{ppp})}\·{\tilde{E}}_m^k;$

其中，所述m为所述步骤C的迭代次数，k为MART迭代次数，当m＝1时，所述为步骤B中得到的空间灰度场。

本发明实施例还提供了一种灰度增强层析粒子图像测速(PIV)重构装置，该装置包括：MART迭代处理模块、粒子拉长抑制模块和虚假粒子抑制模块；其中，

所述MART迭代处理模块，用于对空间灰度场进行MART迭代处理；

所述粒子拉长抑制模块，用于在所述MART迭代处理模块完成MART迭代后，统计空间灰度场中的粒子拉长比；依据所述粒子拉长比确定逆扩散强度因子，并依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场；

所述虚假粒子抑制模块，用于依据统计得到的粒子浓度计算灰度抑制因子，并依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配；之后通知所述MART迭代处理模块进行下一次MART迭代。

本发明实施例还提供了一种灰度增强层析粒子图像测速(PIV)重构设备，该设备包括：上文所述的装置。

本发明实施例提供的灰度增强层析PIV重构方法、装置和设备，在完成MART迭代后，统计空间灰度场中粒子拉长比；依据所述粒子拉长比确定逆扩散强度因子，并依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场；依据统计得到的粒子浓度计算灰度抑制因子，并依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配，之后进行下一次MART迭代。本发明实施例在MART迭代过程中加入了对粒子拉长和虚假粒子的抑制机制，不仅保证了粒子形状的各向同性，而且增加了真实粒子和虚假粒子的灰度对比度，从而减少对速度场测量的影响，有效改善了MART迭代效果。

附图说明

在附图(其不一定是按比例绘制的)中，相似的附图标记可在不同的视图中描述相似的部件。具有不同字母后缀的相似附图标记可表示相似部件的不同示例。附图以示例而非限制的方式大体示出了本文中所讨论的各个实施例。

图1为层析PIV方法的相机布置示意图；

图2为粒子拉长示意图；

图3为虚假粒子示意图；

图4为本发明实施例所述灰度增强层析PIV重构方法实现流程图；

图5为本发明实施例所述灰度增强层析PIV重构装置的结构示意图；

图6为本发明具体实施例所述灰度增强层析PIV重构流程图；

图7为本发明一应用场景中所述粒子浓度为0.15时的重构结果对比图；(a)为已知的正确粒子场；(b)原始MART重构粒子场；(c)重构后的粒子场。

具体实施方式

从现有技术可知，对于粒子拉长，如图2所示，粒子拉长的程度与相机的视角有直接的关系，视角越小，拉长的程度越大。拉长的粒子增加了粒子定位和速度场计算的不确定度。用β表示相机的总视角，d_τ代表粒子图像直径，根据公式d_τz≈d_τ/tan(β/2)可以方便地估计出粒子拉长的初始直径，其中z为重构体厚度方向，图2为该公式的示意图。d_τz/d_τ的值并非固定不变，随着MART迭代，其值会逐渐减小，但在有限次迭代内始终大于1，在有些实验中甚至达到了2。为了更好的描述粒子拉长，单个粒子的灰度需要进行三维高斯拟合：

$I(x,y,z)=I_0\exp [-\frac{{(x-x_0)}^2}{(1/8)d_{\mathrm{\τx}}^2}-\frac{{(y-y_0)}^2}{(1/8)d_{\mathrm{\τy}}^2}-\frac{{(z-z_0)}^2}{(1/8)d_{\mathrm{\τz}}^2}]---\left(1\right)$

其中，所述x,y,z代表重构体的三维空间坐标，z为重构体厚度方向，x₀,y₀,z₀代表粒子中心的物理位置，I₀代表粒子中心的灰度值。从所述公式(1)中可以看出，粒子直径等于4倍高斯拟合的标准差σ。同时，粒子拉长的比例可以用高斯拟合的标准差来表示，即：

$r_d=\frac{d_{\mathrm{\τz}}}{d_{\mathrm{\τ}}}=\frac{{\mathrm{\σ}}_{\mathrm{\τz}}}{{\mathrm{\σ}}_{\mathrm{\τ}}}---\left(2\right)$

由于x,y方向的粒子不受粒子拉长的影响，其直径几乎与真实直径相同，故可认为d_τ＝(d_τx+d_τy)/2。

对于虚假粒子，图3给出了最简单的只包含两个粒子的双相机系统，虚线代表相机视线。黑色菱形代表真实粒子，反投影到体空间后会在视线的四个交 点处都形成粒子，灰色菱形就是虚假粒子。在MART迭代的初始阶段，虚假粒子和真实粒子拥有相同的灰度，但随着迭代的进行，真实粒子的灰度会逐渐超过虚假粒子。从概率统计的意义上可以认为虚假粒子的灰度低于真实粒子，对于单个粒子这一结论可能并不成立。这一事实在实验中已经被证实过多次，只是还未加以利用。虚假粒子受粒子浓度的影响最为明显，粒子浓度越高，增加了相机视线上出现粒子的可能，不论是虚假粒子的个数还是灰度都会快速增加，极大地影响了速度场的测量。

针对上述内容，本发明实施例可分为两部分实现，一部分是运用逆扩散方程在z方向收缩被拉长的粒子，使其逐渐变成球形；另一部分是通过虚假粒子灰度低这一已知事实来重新分配粒子灰度，使得在概率统计上增加真实粒子与虚假粒子的灰度对比度。因此，

本发明的实施例中，在完成MART迭代后，统计空间灰度场中粒子拉长比；依据所述粒子拉长比确定逆扩散强度因子，并依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场；依据统计得到的粒子浓度计算灰度抑制因子，并依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配，之后进行下一次MART迭代。

下面结合附图及具体实施例对本发明作进一步详细说明。

图4为本发明实施例所述灰度增强层析PIV重构方法实现流程图，如图4所示，该方法包括：

步骤401：在完成MART迭代后，统计空间灰度场中粒子拉长比；

步骤402：依据所述粒子拉长比确定逆扩散强度因子，并依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场；

步骤403：依据统计得到的粒子浓度计算灰度抑制因子，并依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配；之后进行下一次MART迭代。

本发明实施例在MART迭代过程中加入了对粒子拉长和虚假粒子的抑制机制，不仅保证了粒子形状的各向同性，而且增加了真实粒子和虚假粒子的灰度 对比度，从而减少对速度场测量的影响，有效改善了MART迭代效果。

本发明实施例中，所述统计空间灰度场中粒子拉长比，包括：

在空间灰度场中采用局部峰值的方法识别粒子，并对这些粒子的灰度做统计平均得到所有粒子灰度的平均结果；对该结果进行三维高斯拟合，得到所述空间灰度场各方向的高斯拟合的标准差，依据所述标准差计算所述粒子拉长比。

本发明实施例中，所述三维高斯拟合按照如下公式进行：

$I(x,y,z)=I_0\exp [-\frac{{(x-x_0)}^2}{(1/8)d_{\mathrm{\τx}}^2}-\frac{{(y-y_0)}^2}{(1/8)d_{\mathrm{\τy}}^2}-\frac{{(z-z_0)}^2}{(1/8)d_{\mathrm{\τz}}^2}]---\left(3\right)$

其中，所述x,y,z代表重构体的三维空间坐标，z为重构体厚度方向，x₀,y₀,z₀代表粒子中心的物理位置，I代表单个粒子的灰度，I₀代表粒子中心的灰度值，d_τx,d_τy,d_τz为x,y,z方向粒子拉长的初始直径。

本发明实施例中，由于x,y方向的粒子不受粒子拉长的影响，其直径几乎与真实直径相同，因此，所述依据空间灰度场各方向的高斯拟合的标准差计算所述粒子拉长比，依据如下公式进行：

$r_d=\frac{d_{\mathrm{\τz}}}{d_{\mathrm{\τ}}}=\frac{{\mathrm{\σ}}_{\mathrm{\τz}}}{{\mathrm{\σ}}_{\mathrm{\τ}}}---\left(4\right)$

其中，所述r_d为粒子拉长比，所述d_τ为粒子图像直径，所述σ_τ为粒子高斯拟合后在x,y方向标准差的平均值，σ_τ＝1/8(d_τx+d_τy)，所述σ_τz为粒子高斯拟合后在z方向的标准差，σ_τz＝1/4d_τz。

这里，由于x,y方向的粒子不受粒子拉长的影响，其直径几乎与真实直径相同，故可认为d_τ＝(d_τx+d_τy)/2。

本发明实施例中，所述依据所述粒子拉长比确定逆扩散强度因子，依据如下公式进行：

$\mathrm{\δ}=p_1\·{r_d}^{p_2}+p_3---\left(5\right)$

其中，所述δ表示逆扩散强度因子，p₁＝-1.40，p₂＝-1.97，p₃＝1.41。由于该公式通过数值模拟扩散过程得到，因此具有普遍适用性。

本发明实施例中，所述依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场，依据如下公式进行：

${\tilde{E}}^k=E^k-\mathrm{\δ}\·E_{\mathrm{zz}}^k---\left(6\right)$

该公式为所述逆扩散方程，所述表示灰度关于z方向的二阶导数，E^k表示所述MART迭代后的灰度，表示在z方向逆扩散后的灰度(由于粒子拉长的问题只涉及z方向)，δ为所述逆扩散强度因子，k为MART迭代次数。

本发明实施例中，所述依据统计得到的粒子浓度计算灰度抑制因子，包括：

依据公式(7)和统计得到的粒子浓度ppp计算得到A_g和B_g，所述A_g和B_g表示虚假粒子分布的参数，利用灰度矩阵的最大值归一化所述粒子灰度场，并结合公式(8)计算得到虚假粒子的概率分布g，所述灰度抑制因子与所述虚假粒子的概率分布之和为1，即：依据α＝1-g计算得到所述灰度抑制因子α，其中，

$\begin{array}{c}{A}_g=a_1\·{\mathrm{ppp}}^{a_2}+a_3\\ B_g=b_1\·\mathrm{ppp}+b_2\end{array}---\left(7\right)$

其中，a₁＝4.17×10^-5，a₂＝-3，a₃＝6.42×10^-2，b₁＝1.82×10^-1，b₂＝2.65×10^-3，ppp为粒子浓度；

g_(Er,ppp)＝A_g·exp(-E_r/B_g)    

E_r＝E/max(E)   (8)

其中，E为所述灰度矩阵，E_r是用最大灰度归一化后的相对灰度，A_g和B_g的取值与粒子浓度ppp有关，如公式(7)所示。

本发明实施例中，所述依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配方法依据如下公式(9)进行：

${\tilde{E}}_{m+1}^k={\mathrm{\α}}_{({\tilde{E}}_m^k,\mathrm{ppp})}\·{\tilde{E}}_m^k---\left(9\right)$

其中，所述m为所述步骤403(不包括之后进行下一次MART迭代)的迭代次数，当m＝1时，所述为步骤402中得到的空间灰度场

本发明实施例还提供了一种灰度增强层析PIV重构装置，如图5所示，该装置包括：MART迭代处理模块501、粒子拉长抑制模块502和虚假粒子抑制 模块503；其中，

所述MART迭代处理模块501，用于对空间灰度场进行MART迭代处理；

所述粒子拉长抑制模块502，用于在所述MART迭代处理模块501完成MART迭代后，统计空间灰度场中的粒子拉长比；依据所述粒子拉长比确定逆扩散强度因子，并依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场；

所述虚假粒子抑制模块503，用于依据统计得到的粒子浓度计算灰度抑制因子，并依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配；之后通知所述MART迭代处理模块501进行下一次MART迭代。

本发明实施例在MART迭代过程中加入了对粒子拉长和虚假粒子的抑制机制，不仅保证了粒子形状的各向同性，而且增加了真实粒子和虚假粒子的灰度对比度，从而减少对速度场测量的影响，有效改善了MART迭代效果。

本发明实施例中，所述粒子拉长抑制模块502统计空间灰度场中粒子拉长比，包括：

在空间灰度场中采用局部峰值的方法识别粒子，并对这些粒子的灰度做统计平均得到所有粒子灰度的平均结果；对该结果进行三维高斯拟合，得到所述空间灰度场各方向的高斯拟合的标准差，依据所述标准差计算所述粒子拉长比。

本发明实施例中，所述粒子拉长抑制模块502进行三维高斯拟合按照如下公式进行：

$I(x,y,z)=I_0\exp [-\frac{{(x-x_0)}^2}{(1/8)d_{\mathrm{\τx}}^2}-\frac{{(y-y_0)}^2}{(1/8)d_{\mathrm{\τy}}^2}-\frac{{(z-z_0)}^2}{(1/8)d_{\mathrm{\τz}}^2}]---\left(5\right)$

其中，所述x,y,z代表重构体的三维空间坐标，z为重构体厚度方向，x₀,y₀,z₀代表粒子中心的物理位置，I代表单个粒子的灰度，I₀代表粒子中心的灰度值，d_τx,d_τy,d_τz为x,y,z方向粒子拉长的初始直径。

本发明实施例中，由于x,y方向的粒子不受粒子拉长的影响，其直径几乎与真实直径相同，因此，所述粒子拉长抑制模块502依据空间灰度场各方向的高斯拟合的标准差计算所述粒子拉长比，依据如下公式进行：

$r_d=\frac{d_{\mathrm{\τz}}}{d_{\mathrm{\τ}}}=\frac{{\mathrm{\σ}}_{\mathrm{\τz}}}{{\mathrm{\σ}}_{\mathrm{\τ}}}---\left(6\right)$

其中，所述r_d为粒子拉长比，所述d_τ为粒子图像直径，所述σ_τ为粒子高斯拟合后在x,y方向标准差的平均值，σ_τ＝1/8(d_τx+d_τy)，所述σ_τz为粒子高斯拟合后在z方向的标准差，σ_τz＝1/4d_τz。

这里，由于x,y方向的粒子不受粒子拉长的影响，其直径几乎与真实直径相同，故可认为d_τ＝(d_τx+d_τy)/2。

本发明实施例中，所述粒子拉长抑制模块502依据所述粒子拉长比确定逆扩散强度因子，依据如下公式进行：

$\mathrm{\δ}=p_1\·{r_d}^{p_2}+p_3---\left(5\right)$

其中，所述δ表示逆扩散强度因子，p₁＝-1.40，p₂＝-1.97，p₃＝1.41。由于该公式通过数值模拟扩散过程得到，因此具有普遍适用性。

本发明实施例中，所述粒子拉长抑制模块502依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场，依据如下公式进行：

${\tilde{E}}^k=E^k-\mathrm{\δ}\·E_{\mathrm{zz}}^k---\left(6\right)$

该公式为所述逆扩散方程，所述表示灰度关于z方向的二阶导数，E^k表示所述MART迭代后的灰度，表示在z方向逆扩散后的灰度，δ为所述逆扩散强度因子，k为MART迭代次数。

本发明实施例中，所述虚假粒子抑制模块503依据统计得到的粒子浓度计算灰度抑制因子，包括：

依据公式(7)和统计得到的粒子浓度ppp计算得到A_g和B_g，利用灰度矩阵的最大值归一化所述粒子灰度场，并结合公式(8)计算得到虚假粒子的概率分布g，所述灰度抑制因子与所述虚假粒子的概率分布之和为1，即：依据α＝1-g计算得到所述灰度抑制因子α，其中，

$\begin{array}{c}{A}_g=a_1\·{\mathrm{ppp}}^{a_2}+a_3\\ B_g=b_1\·\mathrm{ppp}+b_2\end{array}---\left(7\right)$

其中，a₁＝4.17×10^-5，a₂＝-3，a₃＝6.42×10^-2，b₁＝1.82×10^-1，b₂＝2.65×10^-3，ppp为粒子浓度；

$\begin{array}{c}{g}_{(E_r,\mathrm{ppp})}=A_g\·\exp (-E_r/B_g)\\ E_r=E/\max \left(E\right)\end{array}---\left(8\right)$

其中，E为所述灰度矩阵，E_r是用最大灰度归一化后的相对灰度，A_g和B_g的取值与粒子浓度ppp有关，如公式(7)所示。

本发明实施例中，所述虚假粒子抑制模块503依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配方法依据如下公式(9)进行：

${\tilde{E}}^k=E^k-\mathrm{\δ}\·E_{\mathrm{zz}}^k---\left(9\right)$

其中，所述m为所述虚假粒子抑制模块503所执行操作的迭代次数，当m＝1时，所述为粒子拉长抑制模块502中计算得到的空间灰度场

本发明实施例还提供了一种灰度增强层析PIV重构设备，该设备包括：上文所述的装置。

综上所述，对于粒子拉长，本发明实施例从图像处理的角度出发，将粒子拉长看成灰度沿z向的一维高斯平滑，而这种高斯平滑已经被证明可以用扩散方程表示。其中δ是一个始终大于0的常数，代表扩散强度。E_zz是灰度在z方向的二阶导数。当δ等于0时表示灰度没有发生任何变化，当δ大于0时表示灰度在z方向有类似温度的从高灰度区域传播到低灰度区域的效果。相反，其逆扩散方程就可以实现灰度的反向传播，使低灰度向高灰度聚集。本发明实施例在MART的迭代中加入逆扩散方程就可以很好的抑制粒子灰度沿z向的传播，使粒子趋近球形。考虑到程序的稳定性，δ在本发明实施例中的取值范围为0到1。

对于虚假粒子，本发明实施例以虚假粒子灰度低于真实粒子灰度为基础，对不同灰度乘以不同的灰度抑制因子α，即E＝α_(E,ppp)·E，从而达到从概率统计的角度抑制虚假粒子的生长。α与重构体的灰度E和粒子浓度ppp(particle per pixel)都有关系。粒子浓度ppp代表单位像素上粒子的个数，一般通过统计的 方法得到。灰度越低，成为虚假粒子的可能性越大，对其的抑制作用就该越强，也就是说α应该越小。同时，α对灰度的作用范围应该和虚假粒子出现的范围一致，尽量减少对真实粒子(高灰度粒子)的干扰。因此，抑制因子α应该与虚假粒子概率g呈负相关，因此α＝1-g。本发明实施例对不同粒子浓度下的虚假粒子概率g分布做了大量统计模拟，发现其服从公式(8)中所述的指数分布：具体模拟方法为：人工生成四幅不同视角的粒子图像，用原始MART算法迭代五次重构该图像得到粒子灰度，统计其中的虚假粒子分布，并用指数函数拟合得到如公式(8)所示的经验公式。进一步生成不同粒子浓度ppp的粒子图像，找出(A_g,B_g)和ppp的依赖关系，如公式(7)所示。

上述对粒子拉长和虚假粒子的抑制都是在MART的基础上完成的，额外增加的这两步抑制操作和MART本身所需的时间相比非常少，能够在极短的时间内完成对粒子灰度的再次分配，这两步合并到一起统称为灰度增强层析PIV技术。在实际测试中发现，本发明实施例能有效的提高重构精度。

下面结合具体实施例对本发明进行详细描述。

本实施例灰度增强层析PIV重构技术按图6所示流程进行，具体步骤如下：

步骤601：按照MART的公式(2)对流场进行重构得到E^k，并对其进行3×3×3的高斯平滑，平滑强度为0.65；

其中，k为迭代次数，当k为最后一次迭代时，不再进行任何处理。这一步骤完全按照现有MART的算法不进行任何修正。

步骤602：在灰度场E_k中统计粒子拉长比r_d；

具体统计方法为：在E_k中用局部峰值方法识别粒子并对这些粒子的灰度做统计平均，得到所有粒子灰度的平均结果，该结果能够反映粒子拉长的整体情况，对该结果按照公式(3)进行三维高斯拟合，最后将各方向的标准差带入公式(4)得到粒子的拉长比r_d。

步骤603：依据粒子拉长比r_d，根据公式：计算逆扩散强度因子δ，其中p₁＝-1.40，p₂＝-1.97，p₃＝1.41；

该公式通过数值模拟扩散过程得到，具有普适性。具体模拟的方法：人为生成大量标准差为1的3×3二维高斯粒子，用不同δ的扩散方程模拟粒子拉长的过程，最后统计粒子拉长比r_d和δ的关系，用上面的公式进行拟合。

步骤604：将逆扩散强度因子带入逆扩散方程中更新灰度场；

表示灰度关于z向的二阶导数，表示在z向逆扩散后的结果，和E^k相比，在z方向的粒子直径更小了。

步骤605：依据统计得到的粒子浓度计算灰度抑制因子，并依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配；

为了提高对虚假粒子的抑制效果，步骤605可迭代多次，设迭代次数为m，本实施例中m≤20。此步骤的具体实施流程为：

步骤6051：计算灰度抑制因子α；

所述α与每个体素的灰度和整体粒子浓度有关。首先将统计得到的粒子浓度ppp带入公式(7)得到A_g和B_g，然后用灰度矩阵的最大值归一化粒子灰度场 并带入公式(8)得到虚假粒子的概率分布g，最后带入公式α＝1-g即可。

步骤6052：依据公式更新灰度场；

当m＝1时，为步骤604得到的结果如果没有达到迭代次数m则继续步骤6051，如果m已经达到则进入步骤606。

步骤606：完成步骤605后进入下一循环，用替代E^k，重复步骤601。

下面结合具体应用场景对本发明进行描述。

本实例按图1所示配置人工模拟粒子图像投影的过程。重构体大小为35×35×7mm³，每个体素(voxel)的大小为0.05mm，所以整个重构体包含700×700×140个体素。空间四个相机的像素空间为700×700pixel²，对称布置在重构体的上下左右，每台相机的视角都为30度，映射函数通过数学计算得到。空间粒子的直径为3个体素，随机均匀地分布在测量体内，每个粒子的最大灰度都为255。得到真实空间粒子场后，按公式(1)投影生成平面粒子图像，该图像用于空间重构。粒子图像用I1，I2，I3，I4表示，真实空间粒子场用E_real表示，重 构粒子场用E表示，重构次数用k表示。

步骤一：将映射函数和粒子图像I1，I2，I3，I4带入MART重构公式(2)中按步骤进行无修正的粒子场重构，得到当前的重构粒子场E_k。

步骤二：在E_k中用局部3×3×3的区域最大值的方法识别粒子，将所有的粒子灰度相加取平均得到统计平均的结果，该结果可以很好的反映粒子整体被拉长的情况。对该结果按照公式(3)进行三维高斯拟合，最后将各方向的标准差带入公式(4)得到粒子的拉长比r_d。

步骤三：将粒子拉长比r_d带入公式：得出逆扩散强度因子δ。其中p₁＝-1.40，p₂＝-1.97，p₃＝1.41。

步骤四：将步骤三得到的δ带入逆扩散方程中更新灰度场，收缩粒子在z向的灰度。

步骤五：根据步骤605计算灰度抑制因子α并更新灰度场，本实例中该步骤始终迭代更新10次，即m＝10得到灰度增强后的粒子场

步骤六：完成步骤五后，用替代E^k，重复步骤一，直到进行了10次循环。

本实施在层析MART迭代的基础上额外添加了逆扩散方程抑制粒子拉长和虚假粒子的灰度抑制。在相同的迭代次数下，本方法能有效的改善MART迭代的效果。图7是当粒子浓度为0.15时重构结果对比图，图7(a)为已知正确粒子场；图7(b)为原始MART重构粒子场；图7(c)为本实例重构粒子场。从中可以看出本实例不仅消除了粒子拉长效应，同时明显的削弱了虚假粒子，提高了重构的精度。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用硬件实施例、软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器和光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品 的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述，仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。

Claims (10)

1.一种灰度增强层析粒子图像测速PIV重构方法，其特征在于，该方法包括：

A、在完成倍增代数重构技术MART迭代后，统计空间灰度场中的粒子拉长比；

B、依据所述粒子拉长比确定逆扩散强度因子，并依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场；

C、依据统计得到的粒子浓度计算灰度抑制因子，并依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配；

之后进行下一次MART迭代。

2.根据权利要求1所述的方法，其特征在于，所述统计空间灰度场中粒子拉长比，包括：

在空间灰度场中采用局部峰值的方法识别粒子，并对这些粒子的灰度做统计平均得到所有粒子灰度的平均结果；对该结果进行三维高斯拟合，得到所述空间灰度场各方向的高斯拟合的标准差，依据所述标准差计算所述粒子拉长比。

3.根据权利要求2所述的方法，其特征在于，所述三维高斯拟合按照如下公式进行：

其中，所述x,y,z代表重构体的三维空间坐标，z为重构体厚度方向，x₀,y₀,z₀代表粒子中心的物理位置，I代表单个粒子的灰度，I₀代表粒子中心的灰度值，d_τx,d_τy,d_τz为x,y,z方向粒子拉长的初始直径。

4.根据权利要求3所述的方法，其特征在于，所述依据空间灰度场各方向的高斯拟合的标准差计算所述粒子拉长比，依据如下公式进行：

其中，所述r_d为粒子拉长比，所述d_τ为粒子图像直径，所述σ_τ为粒子高斯 拟合后在x,y方向标准差的平均值，σ_τ＝1/8(d_τx+d_τy)，所述σ_τz为粒子高斯拟合后在z方向的标准差，σ_τz＝1/4d_τz。

5.根据权利要求4所述的方法，其特征在于，所述依据所述粒子拉长比确定逆扩散强度因子，依据如下公式进行：

其中，所述δ表示逆扩散强度因子，所述r_d为粒子拉长比，所述p₁＝-1.40，p₂＝-1.97，p₃＝1.41。

6.根据权利要求5所述的方法，其特征在于，所述依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场，依据如下公式进行：

该公式为所述逆扩散方程，所述表示灰度关于z方向的二阶导数，E^k表示所述MART迭代后的灰度，表示在z方向逆扩散后的灰度，δ为所述逆扩散强度因子，k为MART迭代次数。

7.根据权利要求1所述的方法，其特征在于，所述依据统计得到的粒子浓度计算灰度抑制因子，包括：

依据公式(a)和统计得到的粒子浓度计算得到A_g和B_g；利用灰度矩阵的最大值归一化所述粒子灰度场，并结合公式(b)计算得到虚假粒子的概率分布，并依据所述虚假粒子的概率分布计算得到所述灰度抑制因子，所述灰度抑制因子α与所述虚假粒子的概率分布g之和为1；其中，

其中，a₁＝4.17×10^-5，a₂＝-3，a₃＝6.42×10^-2，b₁＝1.82×10^-1，b₂＝2.65×10^-3，ppp为统计得到的粒子浓度；

其中，E为所述灰度矩阵，E_r是用最大灰度归一化后的相对灰度，所述g为 虚假粒子的概率分布。

8.根据权利要求7所述的方法，其特征在于，所述依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配方法，依据如下公式进行：

其中，所述m为所述步骤C的迭代次数，k为MART迭代次数，当m＝1时，所述为步骤B中得到的空间灰度场。

9.一种灰度增强层析粒子图像测速PIV重构装置，其特征在于，该装置包括：MART迭代处理模块、粒子拉长抑制模块和虚假粒子抑制模块；其中，

所述MART迭代处理模块，用于对空间灰度场进行MART迭代处理；

所述粒子拉长抑制模块，用于在所述MART迭代处理模块完成MART迭代后，统计空间灰度场中的粒子拉长比；依据所述粒子拉长比确定逆扩散强度因子，并依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场；

所述虚假粒子抑制模块，用于依据统计得到的粒子浓度计算灰度抑制因子，并依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配；之后通知所述MART迭代处理模块进行下一次MART迭代。

10.一种灰度增强层析粒子图像测速PIV重构设备，其特征在于，该设备包括：权利要求9所述的装置。