CN104777329B - 一种用于粒子图像测速三维粒子场重构的线性规划算法 - Google Patents

Abstract

本发明公开了一种用于粒子图像测速三维粒子场重构的线性规划算法，包括：获取粒子图像中单个粒子图像的灰度分布参数以及多相机投影的权重矩阵；依据所述灰度分布参数构建两组基函数，并将三维灰度场转换为所述两组基函数的线性组合；对所述权重矩阵和转换后的三维灰度场所构成的线性方程进行优化求解。本发明能够实现三维灰度场高精度的重构，且能明显降低虚假粒子出现的概率。

Description

一种用于粒子图像测速三维粒子场重构的线性规划算法

技术领域

本发明涉及流体力学速度测量技术领域，尤其涉及一种用于粒子图像测速三维粒子场重构的线性规划算法。

背景技术

在粒子图像测速(PIV)技术中，采用多台相机从不同视角进行三维流场的示踪粒子成像能最终测量得到三维空间体内的完整三维速度场，常见的层析PIV(TPIV)方法就属于此类技术。TPIV技术中的三维灰度场重构是实现TPIV速度场测量最关键的一个数据处理环节，业界对此展开了大量算法研究。三维灰度场重构是利用在多个相机的投影图像(projection)反推得到三维灰度分布。

用数学的方法来描述三维粒子重构就是一个线性映射关系下的线性方程组求解。具体来说是将包含示踪粒子的测量体Ω离散成一个以体素(voxel)为单位的三维灰度场其中体素为物理空间离散出的立方体，(X,Y,Z)为体素中心的物理坐标，其中带“～”标记的符号表示三维测量空间内的变量，下文同。为方便描述和实际计算，三维灰度场通常用一维链式数组E存储(E∈Rⁿ，n为的体素总数)。一维灰度数组E的第s个元素用E_s表示，表示第s个空间体素的灰度值，其中下标1≤s≤n为体素索引。三维灰度场在第k个相机(1≤k≤N_c，N_c为相机个数)的投影图像记为其中x,y代表像素坐标。类似地，投影图像用一维灰度数组I_k存储(n_k为的像素总数)。一维灰度数组I_k的元素I_k,t表示第k个相机上第t个像素的灰度，其中下标1≤t≤n_k为像素索引。设为多个相机投影图像的灰度数组，I∈R^m,因此E和I之间满足

WE＝I (1)

的关系。其中W为多相机投影的权重矩阵，包含相机几何成像、光学畸变和粒子散射权重等多方面信息，且W∈R^m×n。W和I均为已知量。

三维灰度场重构就是解上述线性方程组WE＝I,(E_s≥0)。在实际测量中，由于采样的相机个数有限，相机像素总数m远小于体素总数n。因此，方程组WE＝I是欠定的，没有定解。通常情况下求解欠定方程组可以通过增加约束，即方程个数，或者采用优化方法来实现。增加约束需要采用更多的相机从不同的视角对粒子场进行成像来增加方程个数，但是由于m<<n使得增加约束的方法基本没有可操作性。现有TPIV的解决思路是将整个问题通过优化方法来进行求解，即引进目标函数，并将方程组WE＝I作为线性约束，求目标函数的极值。而目前最通用的倍增代数重构技术(Multiplicative Algebraic ReconstructionTechnique，简称MART)将最大熵原理中的信息熵作为目标函数，即其中ln表示自然对数。这样整个问题就变成了一个极值问题的求解，即：

求f(E)极值

研究表明，MART算法重构出的单个三维粒子灰度分布在空间上成椭球形，这主要是由于层析PIV不同相机成像视角有限，景深方向上空间分辨率低造成的。而且MART算法在三维灰度场重构中还会产生大量虚假粒子(ghost particle)，这会给速度场计算引入较大误差。目前有大量基于MART算法的层析PIV三维灰度场重构技术，它们都是从计算效率、空间分辨率和消除虚假粒子等几个方面对MART算法进行优化，但是效果都不明显。

发明内容

为解决现有存在的技术问题，本发明实施例提供了一种用于粒子图像测速三维粒子场重构的线性规划算法。

本发明提供了一种用于粒子图像测速三维粒子场重构的线性规划算法，该方法包括：

获取粒子图像中单个粒子图像的灰度分布参数以及多相机投影的权重矩阵；依据所述灰度分布参数构建两组基函数，并将三维灰度场转换为所述两组基函数的线性组合；对所述权重矩阵和转换后的三维灰度场所构成的线性方程进行优化求解。

其中，所述获取粒子图像中单个粒子图像的灰度分布参数的步骤包括：

采用统计平均的方法，从实验获得的粒子图像中提取单个粒子图像的灰度分布参数，所述灰度分布参数包括：粒子图像直径d_τ、灰度二维高斯分布标准差σ和灰度峰值I_m；其中，所述粒子图像直径d_τ取值为：d_τ＝3σ，若d_τ不是整数，则取d_τ的整数部分。

其中，所述两组基函数包括：模板基函数和修正基函数相应的，

所述将三维灰度场转换为所述两组基函数的线性组合，为：

将所述三维灰度场转换为所述模板基函数和修正基函数的线性组合：

其中，所述为模板基函数的第s个基函数，为修正基函数的第s个基函数，系数a_1,s和a_2,s分别为模板基函数和修正基函数的第s个投影系数，且a_1,s∈a₁，a_2,s∈a₂，均为非负数。

其中，所述模板基函数采用三维高斯分布的粒子模板或点传播函数的粒子模板。

其中，所述模板基函数采用三维高斯分布的粒子模板时，所述模板基函数为：

其中，所述C₁是模板基函数的最大灰度，σ是高斯分布的标准差，p＝(i,j,k)代表任意体素点坐标，p_s＝(i_s,j_s,k_s)代表索引为s的体素点坐标，e为自然对数函数的底数，D是粒子灰度分布的直径。

其中，所述模板基函数采用三维高斯分布的粒子模板时，所述修正基函数为一组正交的基函数，即：

其中，所述C₂为常数，C₂取值为C₁＝C₂或取任意值。

其中，所述权重矩阵和转换后的三维灰度场所构成的线性方程为：W(M₁a₁+M₂a₂)＝I。

其中，所述对该线性方程进行优化求解过程中包括：设置目标函数为修正基函数的系数和，即：则该优化求解方法为：

求极小值，满足约束

其中，所述求解方法为：采用线性规划方法进行求解，得到最优投影系数a₁和a₂，获得灰度场E及其三维矩阵形式

本发明实施例提供的用于粒子图像测速三维粒子场重构的线性规划算法，获取粒子图像中单个粒子图像的灰度分布参数以及多相机投影的权重矩阵；依据所述灰度分布参数构建两组基函数，并将三维灰度场转换为所述两组基函数的线性组合；对所述权重矩阵和转换后的三维灰度场所构成的线性方程进行优化求解。本发明实施例所述三维粒子场重构的线性规划算法完全脱离现有的MART算法，能够实现三维灰度场高精度的重构，且能明显降低虚假粒子出现的概率。从实验数据可以看出，本发明实施例的算法能更好的还原粒子分布，即重构粒子场和真实粒子场具有更高的相关系数Q，而且本发明实施例的方法能处理粒子浓度相对更高的工况。

附图说明

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

图1为本发明实施例所述的用于粒子图像测速三维粒子场重构的线性规划算法实现流程示意图；

图2为本发明实施例所述模板基函数和修正基函数的原理示意图；

图3为本发明实施例所述算法与现有MART算法的比较结果示意图。

具体实施方式

本发明的实施例中：获取粒子图像中单个粒子图像的灰度分布参数以及多相机投影的权重矩阵；依据所述灰度分布参数构建两组基函数，并将三维灰度场转换为所述两组基函数的线性组合；对所述权重矩阵和转换后的三维灰度场所构成的线性方程进行优化求解。

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

图1为本发明实施例所述的用于粒子图像测速三维粒子场重构的线性规划算法实现流程示意图，包括：

步骤101：获取粒子图像中单个粒子图像的灰度分布参数以及多相机投影的权重矩阵W；

具体的，通过统计平均的方法，从实验获得的粒子图像中提取单个粒子图像的灰度分布参数，所述参数主要包括：粒子图像直径d_τ、灰度二维高斯分布标准差σ和灰度峰值I_m。所述粒子图像直径可取d_τ＝3σ(d_τ也可取其它值)，若d_τ不是整数，则取整数部分。同时，根据多相机成像进行体PIV标准的三维体标定过程获得投影的权重矩阵W。

步骤102：依据所述灰度分布参数构建两组基函数，并将三维灰度场转换为所述两组基函数的线性组合；

这里，假设体PIV中三维重构粒子的直径和灰度分布上都与粒子图像灰度分布相似，则可以用粒子图像的灰度分布参数来构建模板基函数。根据步骤101获取的粒子图像直径d_τ、标准差σ和灰度峰值I_m来构建模板基函数和修正基函数其中D＝d_τ，C₁＝C₂＝I_m(见后续描述)，并将所述三维灰度场转换为所述模板基函数和修正基函数的线性组合。同时，得到两组基函数的矩阵形式M₁和M₂。

具体的，转换为所述两组基函数的线性组合：其中，所述为根据三维粒子灰度分布的先验知识而确定的基函数，称为模板基函数。为模板基函数的第s个基函数。为修正基函数，而 为修正基函数的第s个基函数。系数a_1,s和a_2,s分别为模板基函数和修正基函数的第s个投影系数，且a_1,s∈a₁，a_2,s∈a₂，均为非负数。

其中，所述模板基函数可以采用不同的粒子灰度分布模板，如三维高斯分布的粒子模板、点传播函数的粒子模板等。下文采用三维高斯分布的粒子灰度模板来说明本发明实施例。通过三维高斯分布的粒子灰度模板来构建模板基函数，其每一个基函数是只含有单个粒子灰度模板的三维灰度场，并且粒子灰度模板中心在第s个体素上，则所述模板基函数即为：

其中C₁是模板基函数的最大灰度，σ是高斯分布的标准差，p＝(i,j,k)代表任意体素点坐标，p_s＝(i_s,j_s,k_s)代表索引为s的体素点坐标，e为自然对数函数的底数，D是粒子灰度分布的直径。

在实际应用中，所述C₁、σ和D由二维投影图像通过统计获得，为已知量。则所述修正基函数为一组正交的基函数，即：

其中，C₂为常数，可以为任意值，通常情况下可以取C₁＝C₂。

由于基函数和与三维灰度场是维度相同的三维数组，同样可用一维链式数组和存储其组成的基函数矩阵分别为M₁和M₂。

步骤103：对所述权重矩阵W和三维灰度场转换所得的方程E构成的线性方程WE＝I进行优化求解；

具体的，所述转换所得的方程为E＝M₁a₁+M₂a₂，因为模板基函数M₁和修正基函数M₂为已知基函数集合，因此优化问题转变成求解投影系数a₁和a₂。此时，未知数的个数变为2n个，而方程的个数不变，依然是m个，因此线性约束

W(M₁a₁+M₂a₂)＝I (5)

变得更加欠定。但是由于M₁为粒子灰度模板，相当于额外添加了约束，因此优化问题采用了更多的先验知识来逼近真实情况。

基于上面的方案，本发明实施例还设置目标函数为修正基函数的系数和，即：在求解过程中求目标函数的极小值，这意味着在极小化修正基的同时，保证极大化模板基函数和真实灰度分布的相关性。以一维离散的灰度分布为例，本发明实施例中所述模板基函数和修正基函数的叠加来拟合真实粒子灰度分布的示意图见图2。图2中真实灰度分布曲线离散到(-1，0，+1)三个网格节点上的灰度等于模板基函数和修正基函数相加的和。图2中所述修正基灰度即为修正基函数的灰度，模板基灰度即为模板基函数的灰度。

最终，优化问题转化为

求极小值，满足约束

其中，除投影系数a₁和a₂外其余都是已知。计算获得投影系数a₁和a₂后根据方程E＝M₁a₁+M₂a₂完成三维粒子灰度场的重构。这里，由于优化问题的目标函数和约束条件W(M₁a₁+M₂a₂)＝I均为线性关系，因此本优化问题可以采用通用的线性规划方法进行求解，如内点法等求解最优投影系数a₁和a₂，获得灰度场E及其三维矩阵形式完成三维粒子重构。

下面结合一具体实施例对本发明实施例所述方法进行论述。

在TPIV算法研究中，大都采用二维人工粒子场进行重构算法的数值实验评估。此时，粒子成像由人工模拟生成，因此粒子图像直径d_τ、标准差σ和灰度峰值I_m都是已知值，而非实际TPIV实验中通过实验数据统计平均分析获得。

本实施例采用较常见的二维人工粒子场来进行算例说明。该二维人工粒子场的算例模拟了四个相机成像，总视角为60°，相邻两个相机的视角为20°，每个相机的像素总数为1020。二维计算域的尺寸为200*1000体素，根据不同的粒子浓度(ppp)生成对应数量的人工粒子随机分布在计算域内，同时生成该粒子场在四个相机上的一维粒子图像投影。本算例测试了ppp＝0.05，0.1，0.15，0.2，0.25，0.3，0.35共7个工况，具体计算步骤如下：

步骤一：在二维空间生成随机分布的粒子灰度场(虽然该算例以二维粒子场为例，但等效于三维问题，因此完全采用上文所述的三维问题的符号系统)。本算例中人工粒子大小为3×3个体素，即d_τ＝3，粒子灰度高斯分布的标准差为σ＝1，峰值灰度最大为I_m＝200。粒子个数为N_p＝1000×ppp，粒子中心(x_p,y_p)随机生成，具有亚像素精度。生成灰度场的投影图及其对应的一维链式数组E和I。根据相机的几何关系获得权重矩阵W，使得I＝WE。

步骤二：依据步骤102所述的方法构造两个基函数矩阵M₁和M₂。由于二维数值实验的粒子是人工生成的，因此粒子大小d_τ，灰度峰值I_m和高斯分布的σ都是已知的，因此模板基函数的粒子直径D为3个体素，高斯分布的σ为1，C₁为200，而修正基函数的C₂也为200。

步骤三：用matlab自带的函数linprog解上述线性规划方程(6)，解得a₁和a₂的最优解。

步骤四：通过投影系数a₁和a₂得到重构灰度场E及其二维形式

图3给出了本发明实施例所述算法与现有MART算法的比较结果示意图(横坐标为被测粒子场浓度(particle per pixel，ppp)，纵坐标为被测粒子场和重构粒子场的相关系数Q)。从图3中可以看出，本发明实施例的算法能更好的还原粒子分布，即重构粒子场和真实粒子场具有更高的相关系数Q。且本发明的算法能处理粒子浓度更高的工况，即PPP大于0.3。

可见，本发明所述三维粒子场重构的线性规划算法完全脱离现有的MART算法，能够实现三维灰度场高精度的重构，且能明显降低虚假粒子出现的概率。本发明基于体PIV中示踪粒子几何和光学成像特性的先验知识，即所有三维粒子灰度分布符合三维高斯分布，首次提出了将三维灰度场进行模态分解并采用线性规划的方法进行求解。

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

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

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

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

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

Claims (9)

1.一种用于粒子图像测速三维粒子场重构的线性规划算法，其特征在于，该方法包括：

获取粒子图像中单个粒子图像的灰度分布参数以及多相机投影的权重矩阵；依据所述灰度分布参数构建两组基函数，并将三维灰度场转换为所述两组基函数的线性组合；对所述权重矩阵和转换后的三维灰度场所构成的线性方程进行优化求解。

2.根据权利要求1所述的用于粒子图像测速三维粒子场重构的线性规划算法，其特征在于，所述获取粒子图像中单个粒子图像的灰度分布参数的步骤包括：

采用统计平均的方法，从实验获得的粒子图像中提取单个粒子图像的灰度分布参数，所述灰度分布参数包括：粒子图像直径d_τ、灰度二维高斯分布标准差σ和灰度峰值I_m；其中，所述粒子图像直径d_τ取值为：d_τ＝3σ，若d_τ不是整数，则取d_τ的整数部分。

3.根据权利要求1或2所述的用于粒子图像测速三维粒子场重构的线性规划算法，其特征在于，所述两组基函数包括：模板基函数和修正基函数相应的，

所述将三维灰度场转换为所述两组基函数的线性组合，为：

将所述三维灰度场转换为所述模板基函数和修正基函数的线性组合：

其中，所述为模板基函数的第s个基函数，为修正基函数的第s个基函数，系数a_1,s和a_2,s分别为模板基函数和修正基函数的第s个投影系数，且a_1,s∈a₁，a_2,s∈a₂，均为非负数，n为的体素总数，s为体素索引，a₁和a₂分别为模板基函数和修正基函数的投影系数。

4.根据权利要求3所述的用于粒子图像测速三维粒子场重构的线性规划算法，其特征在于，所述模板基函数采用三维高斯分布的粒子模板或点传播函数的粒子模板。

5.根据权利要求4所述的用于粒子图像测速三维粒子场重构的线性规划算法，其特征在于，所述模板基函数采用三维高斯分布的粒子模板时，所述模板基函数为：

其中，所述C₁是模板基函数的最大灰度，σ是高斯分布的标准差，p＝(i,j,k)代表任意体素点坐标，p_s＝(i_s,j_s,k_s)代表索引为s的体素点坐标，e为自然对数函数的底数，D是粒子灰度分布的直径。

6.根据权利要求5所述的用于粒子图像测速三维粒子场重构的线性规划算法，其特征在于，所述模板基函数采用三维高斯分布的粒子模板时，所述修正基函数为一组正交的基函数，即：

其中，所述C₂为常数，C₂取值为C₁＝C₂或取任意值。

7.根据权利要求3所述的用于粒子图像测速三维粒子场重构的线性规划算法，其特征在于，所述权重矩阵和转换后的三维灰度场所构成的线性方程为：W(M₁a₁+M₂a₂)＝I；

其中，W为权重矩阵，M₁为所述模板基函数的矩阵形式，M₂为所述修正基函数的矩阵形式，I为多个相机投影图像的灰度数组。

8.根据权利要求7所述的用于粒子图像测速三维粒子场重构的线性规划算法，其特征在于，所述对该线性方程进行优化求解过程中包括：设置目标函数为修正基函数的系数和，即：则该优化求解方法为：

求极小值，满足约束

9.根据权利要求8所述的用于粒子图像测速三维粒子场重构的线性规划算法，其特征在于，所述求解方法为：采用线性规划方法进行求解，得到最优投影系数a₁和a₂，获得灰度场E及其三维矩阵形式