CN115035241B - 基于局部基函数的多方向三维背景纹影层析重建装置及方法 - Google Patents

Abstract

本发明公开了一种基于局部基函数的多方向三维背景纹影层析重建装置及方法，本发明基于背景纹影技术，通过对比有流场和无流场时的背景图像得到背景图像上每个点在流场影响下的偏移，从而进行流场重建；多方向三维背景纹影层析重建装置包含背景板图案、背景纹影图像采集模块以及测量区域，装置简单、易于装调并且测量范围大，有力地支持了背景纹影技术的数据采集和后续的重建。本发明还公开了基于局部基函数的三维背景纹影层析重建方法，其通过建立基于局部基函数的背景纹影投影矩阵，并在代数重建求解过程中使用全变差正则化对结果进行约束，可以定量求解出三维空间区域中的折射率分布，进而可以得到密度分布和温度等参数。

Description

基于局部基函数的多方向三维背景纹影层析重建装置及方法

技术领域

本发明属于流场测量领域，特别是涉及一种基于局部基函数的多方向三维背景纹影层析重建装置及方法。

背景技术

复杂流场的流动显示和关键参数的定量测量可以呈现流动规律和揭示流动现象，同时也能促进现代航空航天、工业设计和能源与环境等领域的发展。其中，基于光学的测量方法具有非接触性、时空分辨率高等优点，以及被用于观测燃烧场、风洞实验和气体射流等流场。其中2000年提出的背景纹影层析技术(Background-oriented schlierentomography,BOST)不需要精密、昂贵的光学仪器，抗震动能力较强。BOS的装置更为灵活，视场较大，且易于获得多角度的纹影图像，容易实现定量化的测量，在国内外的研究中已有广泛的应用。

与传统纹影技术一样，BOS利用光线经过变折射率场产生的偏转进行测量，见图1。光线的偏转角与光线传播路径上的折射率梯度的积分成正比：

从上式中可以看出，BOST与传统的层析重建不同，若直接重建得到的是折射率梯度的分布。如何获得折射率的分布是BOST的研究重点。在BOS发展的早期，一般将BOST的重建和获得折射率分布分为两步，即首先进行重建获得折射率梯度分布，再通过求解泊松方程或者直接线积分的方法获得折射率的分布。但是，求解泊松方程的方法较为复杂并且会使得结果过于平滑，而直接线积分的方法较为简单却会在积分过程中造成噪声的累计。处理(1)式中的并将(1)式离散化为线性方程组是直接重建获得折射率分布的关键步骤。因此，为了直接重建获得折射率的分布，有学者在(1)式的离散化过程中引入了有限差分矩阵，使用差分来代替求导，将(1)式离散化为：

ε_u＝An＝TD_un. (2)

式中，A为BOST投影矩阵；T为BOST权重矩阵；D_u为差分矩阵。但是，差分并不是准确的求导过程，其精度受所用的差分格式影响。

所以，有必要设计一种新的BOS投影模型的离散算法，可以进行重建区域中每一点的精确求导，方便计算BOST的投影矩阵，最后可以通过求解线性方程组的形式直接得到折射率分布，而无需额外的步骤。

发明内容

本发明的目的是针对在现有技术的不足，提供一种基于局部基函数的多方向三维背景纹影层析重建装置及方法，主要针对BOST投影模型的离散化过程，使得重建可以直接获得三维流场的分布。

实现本发明目的的技术解决方案为：第一方面，本发明提供一种基于局部基函数的多方向三维背景纹影层析重建装置，包括：

背景板图案，由印在平面平板上的图案组成，用于显示流场的变化；

背景纹影图像采集模块，由多个相机组成，用于记录所述背景板图案在流场影响下的变化；

测量区域，位于所述背景板图案和所述背景纹影图像采集模块之间，用于限定被测流场的范围。

在其中一个实施例中，所述的背景纹影图像采集模块的多个相机均匀地围绕所述测量区域排列在一个圆弧上。

在其中一个实施例中，所述的背景板图案为横竖交叉的正/余弦条纹、随机点阵或者其它随机噪声图案。

第二方面，本发明提供一种基于第一方面所述装置的重建方法，重建装置用于记录多方向背景纹影图像，重建方法从多方向纹影图像中计算重建出流场的三维折射率分布，包括基于局部基函数的BOST投影矩阵的计算、BOST的重建求解以及全变差正则化，具体步骤如下：

1)使用背景纹影图像采集模块分别采集无流场和有流场时的背景板图案的图像，无流场时的图像作为参考图像，有流场时的图像作为实验图像；

2)对比参考图像和实验图像，使用互相关算法或者光流算法获得图像偏移量；

3)将图像偏移量转化为光线偏转角；

4)计算BOST投影矩阵，将待测场的三维折射率分布用局部基函数展开；

5)通过求解线性方程组获得各个网格上的局部基函数系数，从而获得折射率的三维分布；

6)ART迭代求解过程中采用全变差正则化对结果进行约束。

第三方面，本发明提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现第二方面所述的方法的步骤。

第四方面，本发明提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现第二方面所述的方法的步骤。

与现有技术相比，本发明的显著优点为：通过利用局部基函数的可微分性，通过其微分的解析形式来求解BOST的投影矩阵，从而将BOST中沿光线路径的偏导积分的偏导项“消融”掉，构建BOST的投影矩阵。通过我们的方法可以直接计算出待测场的折射率，既不需要再去对折射率梯度进行积分，也不需要额外的有限差分矩阵。

附图说明

图1是本发明背景纹影技术原理图。

图2是本发明三维背景纹影图像采集装置图。

图3为本发明基于局部基函数的多方向三维背景纹影层析重建方法流程图。

图4是本发明实施例中提取的背景图像的偏移矢量。

图5中的(a)是本发明实施例重建结果沿着纵向的截面。

图5中的(b)是本发明实施例重建结果沿着横向的截面。

图6是本发明实施例不同时刻重建出来的三维空间中的温度分布。

具体实施方式

如图1所示，本发明提供一种基于局部基函数的多方向三维背景纹影层析重建装置，包括：

背景板图案，由印在平面平板上的图案组成，用于显示流场的变化；所述的背景板图案为横竖交叉的正/余弦条纹、随机点阵或者其它随机噪声图案；

背景纹影图像采集模块，由多个相机组成，多个相机均匀地围绕所述测量区域排列在一个圆弧上，用于记录所述背景板图案在流场影响下的变化；

测量区域，位于所述背景板图案和所述背景纹影图像采集模块之间，用于限定被测流场的范围。

第二方面，本发明提供一种基于局部基函数的多方向三维背景纹影层析重建方法，如图2所示，所述方法从多方向纹影图像中计算重建出流场的三维折射率分布，包括基于局部基函数的BOST投影矩阵的计算、BOST的重建求解以及全变差正则化，具体步骤如下：

1)使用背景纹影图像采集模块分别采集无流场和有流场时的背景板图案的图像，无流场时的图像作为参考图像，有流场时的图像作为实验图像；由于流场中存在折射率的变化，因此光线穿过流场后会发生偏转，那么实验图像就会发生变形。

2)对比参考图像和实验图像，使用互相关算法或者光流算法获得图像偏移量；即图3中的Δx和Δy。

3)将图像偏移量转化为光线偏转角，在旁轴近似条件下有：

式中，l_A为背景板图案与测量区域之间的距离；l_B为测量区域与背景纹影图像采集模块之间的距离；l_C为相机焦距；x，y和z分别为坐标系的三个方向；

4)计算BOST投影矩阵，将待测场的三维折射率分布用局部基函数展开，假设折射率分布为n(x,y,z)，选取的局部基函数为g(x,y,z)，则待测场可以表示为：

式中，m、n和p分别为x,y和z方向上的网格数；d为每个网格的长度；

g(x-id,y-jd,z-kd)表示(i,j,k)网格上的局部基函数；f_ijk为系数。根据BOS理论可得：

式中，ε_u为u方向的偏折角。可以将上式写成矩阵的形式：

ε＝Tf (4)

其中，ε为各采样点的偏移量在三个轴上的分量；f为各局部基函数对应的系数值f_ijk，即为方程组待求解的未知数；n₀为环境折射率；T为投影矩阵，其中元素的具体形式为：

其中，表示编号为(i,j,k)网格u方向上的权重值。若所有图像一共有K个采样点，测量区域被划分为M×N×P个网格，则ε的大小为3K×1；f的大小为MNP×1；T的大小为3K×MNP；

5)通过求解(4)式的线性方程组则可以获得各个网格上的局部基函数系数f_ijk，从而获得折射率的三维分布。由于(4)式一般是病态的线性方程组，因此，求解过程可以采用代数重建法(Algebra reconstruction technique,ART)进行求解，ART迭代过程遵循二范数最小原则，其迭代过程为：

其中f_i ^(k,j)的上标表示第k次迭代的第j条射线，下标i为第i个网格，λ为松弛因子；ε_j第j条射线的偏移量；t_ij为第j条射线，第i个网格的权重值。

6)由于ART得到的结果噪声比较多，因此在ART迭代求解过程中采用了全变差(Total variation,TV)正则化对结果进行约束。流场折射率全变差的离散形式可以写为：

其中，为折射率在i方向上的偏导数；在每一次ART迭代后，都对其结果进行TV正则化约束，即最小化TV(n)：

其中，β为松弛因子；n^*和n分别为正则化后和正则化前的折射率；n'为进行正则化的变量。

下面结合附图和实施例描述对本发明的装置及方法进行详细说明。

首先，需要搭建三维BOS图像采集装置，如图1所示，测量区域位于背景板图案和图像采集模块之间，待测流场位于测量区域中，图像采集模块由12个CCD相机组成，12个CCD相机均匀地分布在大约160°的圆弧上，并且12个CCD相机都聚焦在背景板图案上。

使用多相机标定方法对12个CCD相机进行标定，目的是为了在同一世界坐标系下求得所有相机参数，包括：内部参数、旋转矩阵和平移向量。世界坐标系中的点、相机坐标系中的点和图像坐标系中的点的齐次坐标分别可以表示为M_w＝(x_w,y_w,z_w,1)^T，M_c＝(x_c,y_c,z_c,1)^T和m＝(u,v,1)^T。从世界坐标点到图像点的成像过程为：

sm＝PM_w,

式中，R为3×3旋转矩阵；t为3×1平移向量；K为相机内部参数；P为相机投影矩阵。得到相机参数后，就可以计算出图像中一点(u,v)对应的光线方程为：

因此，得到光线方程后，就可以根据(5)式和图像中选取的采样点计算出BOST投影矩阵。

对比参考图像和实验图像，使用互相关算法或者光流算法可以获得图像中采样点的偏移量。由于光线的偏转角是三维的，即一条光线有三个偏转角分量(ε_x,ε_y,ε_z)，但是图像采用点的偏移量是二维的，因此需要将偏移量分解到世界坐标系的三个方向上，分解过程为：

Δ_x＝r₁₁Δx+r₂₁Δy

Δ_y＝r₁₂Δx+r₂₂Δy

Δ_z＝r₁₃Δx+r₂₃Δy

并且在旁轴近似的条件下，偏移量与偏转角的关系为：

得到BOST投影矩阵和光线偏转角后就可以使用ART算法进行重建。在重建过程中，为了抑制噪声在每次ART迭代后都进行了TV正则化。TV的最小化过程具体采用的是分裂布雷格曼算法。

本实施例中背景板图案与测量区域的距离为l_A＝750mm，测量区域与图像采集模块的距离为l_B＝650mm，相机焦距为l_C＝f＝650mm。使用该装置对蜡烛火焰上方的热气流进行了测量。测量区域位于蜡烛火焰的上方，测量区域大小为120×120×192mm，并且被划分为200×200×320个网格。使用的背景板图案为横竖交叉的正/余弦条纹。

本实施例中，在每个相机图像的中间选取了185×126个采样点，采样点的间隔为5个像素，采样点的偏移量如图4所示。选取完采样点后，可以根据前述方法计算出采样点对应光线的方程。从而可以沿着光线路径按照前述方法计算出BOST投影矩阵。在本实施例中，投影矩阵的大小为839160×12800000，是一个比较庞大的矩阵，因此为了减少内存消耗，本实施例中采用了稀疏矩阵三元组的存储方式。

用ART迭代和TV正则化相结合的方式可以计算出折射率分布，本实施例中迭代200次可以收敛。求出折射率后，可以通过Gladstone–Dale公式和理想气体状态方程由折射率计算出温度分布：

式中，n为折射率；ρ为密度；K_GD为Gladstone-Dale常数，对空气通常取2.26×10^{- 4}m³/kg；P为压强；M为气体摩尔质量；R为大气常数。

计算出的温度分布如图5和图6所示，图5显示了一个重建结果不同位置截面的温度分布，图6显示了不同时刻重建出来的三维空间中的温度分布。从图中可以看出，本实施例采用的基于局部基函数的多方向三维背景纹影层析重建方法可以重建出复杂的流场结构。

Claims (6)

1.一种基于局部基函数的多方向三维背景纹影层析重建方法，其特征在于，该方法基于多方向三维背景纹影层析重建装置实现，装置包括：背景板图案，由印在平面平板上的图案组成，用于显示流场的变化；背景纹影图像采集模块，由多个相机组成，用于记录所述背景板图案在流场影响下的变化；测量区域，位于所述背景板图案和所述背景纹影图像采集模块之间，用于限定被测流场的范围；

该方法从多方向纹影图像中计算重建出流场的三维折射率分布，包括基于局部基函数的BOST投影矩阵的计算、BOST的重建求解以及全变差正则化，具体步骤如下：

1)使用背景纹影图像采集模块分别采集无流场和有流场时的背景板图案的图像，无流场时的图像作为参考图像，有流场时的图像作为实验图像；

2)对比参考图像和实验图像，使用互相关算法或者光流算法获得图像偏移量；

3)将图像偏移量转化为光线偏转角，在旁轴近似条件下有：

式中，l_A为背景板图案与测量区域之间的距离；l_B为测量区域与背景纹影图像采集模块之间的距离；l_C为相机焦距；ε_u为u方向的偏折角；x，y和z分别为坐标系的三个方向；

4)计算BOST投影矩阵，将待测场的三维折射率分布用局部基函数展开；假设折射率分布为n(x,y,z)，选取的局部基函数为g(x,y,z)，则待测场表示为：

式中，m、n和p分别为x,y和z方向上的网格数；d为每个网格的长度；

g(x-id,y-jd,z-kd)表示(i,j,k)网格上的局部基函数；f_ijk为系数；根据BOS理论得：

将上式写成矩阵的形式：

ε＝Tf (4)

其中，ε为各采样点的偏移量在三个轴上的分量；f为各局部基函数对应的系数值f_ijk，即为方程组待求解的未知数；n₀为环境折射率；T为投影矩阵，其中元素的具体形式为：

其中，表示编号为(i,j,k)网格u方向上的权重值；若所有图像一共有K个采样点，测量区域被划分为M×N×P个网格，则ε的大小为3K×1；f的大小为MNP×1；T的大小为3K×MNP；

5)通过求解(4)式的线性方程组则获得各个网格上的局部基函数系数f_ijk，从而获得折射率的三维分布；求解过程采用代数重建法进行求解，ART迭代过程遵循二范数最小原则，其迭代过程为：

其中f_i ^(k,j)的上标表示第k次迭代的第j条射线，下标i为第i个网格；λ为松弛因子；ε_j第j条射线的偏移量；t_ij为第j条射线，第i个网格的权重值；

6)ART迭代求解过程中采用全变差正则化对结果进行约束。

2.根据权利要求1所述的方法，其特征在于，所述的背景纹影图像采集模块的多个相机均匀地围绕所述测量区域排列在一个圆弧上。

3.根据权利要求2所述的方法，其特征在于，所述的背景板图案为横竖交叉的正/余弦条纹、随机点阵或者其它随机噪声图案。

4.根据权利要求1所述的方法，其特征在于，ART迭代求解过程中采用全变差正则化对结果进行约束；流场折射率全变差的离散形式写为：

其中，为折射率在i方向上的偏导数；在每一次ART迭代后，均对其结果进行TV正则化约束，即最小化TV(n)：

其中，β为松弛因子；n^*和n分别为正则化后和正则化前的折射率；n'为进行正则化的变量。

5.一种计算机设备，包括存储器、处理器及存储在存储器上并在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现权利要求1-4中任一所述的方法的步骤。

6.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，该程序被处理器执行时实现权利要求1-4中任一所述的方法的步骤。