一种灰度增强层析piv重构方法、装置和设备的制造方法
【技术领域】
[0001] 本发明涉及基于灰度增强的层析粒子图像测速技术，尤其涉及一种灰度增强层析 粒子图像测速（Particle Image Velocimetry,简称PIV)重构方法、装置和设备。
【背景技术】
[0002] PIV是一种现代激光测速技术，主要运用于流场速度测量，通过追踪示踪粒子在流 场中的运动来得到速度场。最近兴起的层析粒子图像测速技术（层析PIV)成功地将二维 PIV推广到三维流场测量，能够获得瞬时的三维三分量（3D3C)的速度场。该技术通过不同 视角下（一般为4个相机）粒子散射成像重构出空间粒子的真实分布，然后采用三维互相 关计算相邻曝光粒子间的位移。空间粒子重构是该方法的关键。
[0003] 层析PIV的相机布置如图1所不，其中相机为'+'字型布置。空间测量体E中的 示踪粒子被激光照亮以后按照投影关系同时成像在四个不同视角的相机之上。投影成像的 粒子灰度一般呈高斯分布，大小与示踪粒子大小和散射强度有关，一般占3X3像素。从测 量空间到相机平面的投影关系可以事先标定得到，映射函数决定了粒子之间的相对位置。 如果将空间测量体E离散成和像素（pixel)大小相当的体素（voxel)，投影成像可以简化成 一系列线性方程组：
[0004] WE = I (1)
[0005] 即假设粒子图像是空间粒子灰度沿视线投影积分的结果。其中W为投影权重函 数，Wu表示第j个体素对第i个像素的贡献，该贡献与投影距离有关，距离越小贡献越大。 空间灰度重构可以认为是一个逆投影过程，就是已知粒子图像I和权重函数W求空间灰度 分布E。由于离散的体素个数远远大于已知的像素个数，所以该逆问题存在不定解。为了 求解该方程组，需要额外添加约束条件。目前，主流的层析PIV都采用基于熵最大的倍增代 数重构技术（multiplicative algebraic reconstruction technique，MART)对粒子场进 行重构还原。该方法被证明是当前最适合层析PIV重构的，该方法迭代收敛速度快，重构的 精度较高，其迭代公式可简单表述为：
[0006] Ek+1= Ek(I/ffEk)uW (2)
[0007] 其中，k代表迭代次数，I/WEk代表对投影误差的评估，空间灰度E会根据误差的大 小不断逼近最优解。
[0008] 然而，在实际的使用中发现，重构精度会受到相机个数与视角，粒子浓度和映射函 数精度等因素的影响，出现粒子拉长和虚假粒子两个非常棘手的问题。所谓粒子拉长，是指 层析PIV在理想情况下重构出来应该是高斯分布的球形粒子，正如平面PIV中的粒子服从 高斯分布一样。但由于相机视角以及个数原因，投影灰度并不能完全反应粒子的空间形状， 导致了在相机轴线方向（测量体厚度方向）粒子被拉长。
[0009] 此外，虚假粒子也是层析PIV重构面临的一个难点。图像采集实际上是一个投影 积分过程，空间连续的粒子灰度被映射到离散的平面图像之上。MART算法就是希望通过这 样的欠采样图像重构出空间粒子灰度分布，这是反投影过程。由于可利用的粒子位置信息 有限，导致在所有相机视线的交点处都可能出现粒子，但是这个粒子是否真实存在就不得 而知。
[0010] 通过上面的分析可知，层析PIV方法中的MART算法并没有利用粒子形状、灰度等 许多有用的信息，导致了其自身更新速度很慢，甚至不正确。
【发明内容】
[0011] 为解决现有存在的技术问题，本发明实施例提供一种灰度增强层析PIV重构方 法、装置和设备。
[0012] 本发明实施例提供了一种灰度增强层析粒子图像测速（PIV)重构方法，该方法包 括：
[0013] A、在完成倍增代数重构技术MART迭代后，统计空间灰度场中的粒子拉长比；
[0014] B、依据所述粒子拉长比确定逆扩散强度因子，并依据所述逆扩散强度因子和逆扩 散方程更新所述空间灰度场；
[0015] C、依据统计得到的粒子浓度计算灰度抑制因子，并依据所述灰度抑制因子对所述 已更新的空间灰度场的粒子灰度进行重新分配；
[0016] 之后进行下一次MART迭代。
[0017] 其中，所述统计空间灰度场中粒子拉长比，包括：
[0018] 在空间灰度场中采用局部峰值的方法识别粒子，并对这些粒子的灰度做统计平均 得到所有粒子灰度的平均结果；对该结果进行三维高斯拟合，得到所述空间灰度场各方向 的高斯拟合的标准差，依据所述标准差计算所述粒子拉长比。
[0019] 其中，所述三维高斯拟合按照如下公式进行：
【主权项】
1. 一种灰度增强层析粒子图像测速Piv重构方法，其特征在于，该方法包括： A、 在完成倍增代数重构技术MART迭代后，统计空间灰度场中的粒子拉长比； B、 依据所述粒子拉长比确定逆扩散强度因子，并依据所述逆扩散强度因子和逆扩散方 程更新所述空间灰度场； C、 依据统计得到的粒子浓度计算灰度抑制因子，并依据所述灰度抑制因子对所述已更 新的空间灰度场的粒子灰度进行重新分配； 之后进行下一次MART迭代。
2. 根据权利要求1所述的方法，其特征在于，所述统计空间灰度场中粒子拉长比，包 括： 在空间灰度场中采用局部峰值的方法识别粒子，并对这些粒子的灰度做统计平均得到 所有粒子灰度的平均结果；对该结果进行三维高斯拟合，得到所述空间灰度场各方向的高 斯拟合的标准差，依据所述标准差计算所述粒子拉长比。
3. 根据权利要求2所述的方法，其特征在于，所述三维高斯拟合按照如下公式进行：
其中，所述X，y，z代表重构体的三维空间坐标，z为重构体厚度方向，Xtl, y(l，2(|代表粒子 中心的物理位置，I代表单个粒子的灰度，Itl代表粒子中心的灰度值，d τχ，cUy，cUz为X，y，z 方向粒子拉长的初始直径。
4. 根据权利要求3所述的方法，其特征在于，所述依据空间灰度场各方向的高斯拟合 的标准差计算所述粒子拉长比，依据如下公式进行：
其中，所述rd为粒子拉长比，所述(1,为粒子图像直径，所述〇 τ为粒子高斯拟合后在 x，y方向标准差的平均值，〇 τ = l/8(d TX+cUy)，所述〇 τζ为粒子高斯拟合后在ζ方向的标 准差，σ τζ= l/4d τζ。
5. 根据权利要求4所述的方法，其特征在于，所述依据所述粒子拉长比确定逆扩散强 度因子，依据如下公式进行：
其中，所述S表示逆扩散强度因子，所述rd为粒子拉长比，所述Ρι= -1.40, ρ2 =-L 97, p3= L 41〇
6. 根据权利要求5所述的方法，其特征在于，所述依据所述逆扩散强度因子和逆扩散 方程更新所述空间灰度场，依据如下公式进行：
该公式为所述逆扩散方程，所述（表示灰度关于z方向的二阶导数，Ek表示所述MART 迭代后的灰度，扔表示在z方向逆扩散后的灰度，δ为所述逆扩散强度因子，k为MART迭 代次数。
7. 根据权利要求1所述的方法，其特征在于，所述依据统计得到的粒子浓度计算灰度 抑制因子，包括： 依据公式（a)和统计得到的粒子浓度计算得到AjP Bg;利用灰度矩阵的最大值归一化 所述粒子灰度场，并结合公式（b)计算得到虚假粒子的概率分布，并依据所述虚假粒子的 概率分布计算得到所述灰度抑制因子，所述灰度抑制因子α与所述虚假粒子的概率分布g 之和为1 ;其中，
其中，B1= 4· 17ΧΚΓ5, a2= -3, a3= 6· 42ΧΚΓ2, Id1= I. 82ΧΚΓ1，b2= 2· 65ΧΚΓ3， PPP为统计得到的粒子浓度；
其中，E为所述灰度矩阵，民是用最大灰度归一化后的相对灰度，所述g为虚假粒子的 概率分布。
8. 根据权利要求7所述的方法，其特征在于，所述依据所述灰度抑制因子对所述已更 新的空间灰度场的粒子灰度进行重新分配方法，依据如下公式进行：
其中，所述m为所述步骤C的迭代次数，k为MART迭代次数，当m = 1时，所述^为步 骤B中得到的空间灰度场。
9. 一种灰度增强层析粒子图像测速PIV重构装置，其特征在于，该装置包括：MART迭代 处理模块、粒子拉长抑制模块和虚假粒子抑制模块；其中， 所述MART迭代处理模块，用于对空间灰度场进行MART迭代处理； 所述粒子拉长抑制模块，用于在所述MART迭代处理模块完成MART迭代后，统计空间灰 度场中的粒子拉长比；依据所述粒子拉长比确定逆扩散强度因子，并依据所述逆扩散强度 因子和逆扩散方程更新所述空间灰度场； 所述虚假粒子抑制模块，用于依据统计得到的粒子浓度计算灰度抑制因子，并依据 所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配；之后通知所述 MART迭代处理模块进行下一次MART迭代。
10. -种灰度增强层析粒子图像测速PIV重构设备，其特征在于，该设备包括：权利要 求9所述的装置。
【专利摘要】本发明公开了一种灰度增强层析粒子图像测速(PIV)重构方法，该方法包括：在完成倍增代数重构技术(MART)迭代后，统计空间灰度场中粒子拉长比；依据所述粒子拉长比确定逆扩散强度因子，并依据所述逆扩散强度因子和逆扩散方程更新所述空间灰度场；依据统计得到的粒子浓度计算灰度抑制因子，并依据所述灰度抑制因子对所述已更新的空间灰度场的粒子灰度进行重新分配，之后进行下一次MART迭代。本发明还同时公开了一种实现所述方法的装置和设备。
【IPC分类】G06T11-00
【公开号】CN104680564
【申请号】CN201510110129
【发明人】高琪, 王洪平, 王晋军, 魏润杰
【申请人】北京航空航天大学, 魏润杰
【公开日】2015年6月3日
【申请日】2015年3月12日