Characterization of thermodynamic non-equilibrium of plasma flow using coherent anti-Stokes Raman scattering
-
摘要:
流场温度是流体最基本的参数之一，流场温度测量精度直接影响到高超声速飞行器气动力、气动热以及热防护等性能预测准确性。本文基于CARS技术基本原理，开发了面向非平衡流场的CARS光谱计算及振转温度反演算法，并在
1000 K～2300 K静态环境下对其准确性开展了实验验证。搭建了基于微波等离子体的非平衡流场，获取了不同压强、流量以及组分下等离子体流场的振转温度信息。结果显示在实验条件范围内，振动温度随着压强的增加而降低，随N2体积流量和Ar体积分数呈现先增大后降低的趋势；转动温度随着压强的增加而降低，随N2体积流量呈现先增大后降低的趋势，随Ar体积分数的增加而增大；非平衡度随压强的增加而增加但其变化率随压强增加而减小，随N2体积流量的增加呈现先降低后增大的趋势，随Ar体积分数的增加而降低。-
关键词:
- 相干反斯托克斯拉曼散射 /
- 热力学非平衡 /
- 等离子体流场 /
- 振动温度 /
- 转动温度
Abstract:Temperature is considered to be one of the most concerned parameters to quantitatively describe flow characteristics, of which the measurement accuracy directly affects the prediction of aerodynamic, aerothermal and thermal protection performance of hypersonic vehicles. Based on the principles of Coherent Anti-Stokes Raman Scattering (CARS), a CARS spectral computation and vib-rotational temperature inversion program is proposed for characterizing the thermodynamic non-equilibrium properties of the high-temperature gas flow field. And corresponding accuracy from 1000 K to 2300 K is verified in a static environment. A non-equilibrium microwave plasma flow is built and its vibrational temperature and rotational temperature with different pressures, N2 volumetric flow rate, and compositions are obtained by using the developed program. The results show that within the range of experimental conditions, with pressure increasing, the vibrational temperature and rotational temperature decrease, while the thermodynamic non-equilibrium degree increases but corresponding increase rate decreases. With N2 volumetric flow rate increasing, the vibrational temperature and rotational temperature first increase and then decrease, while the thermodynamic non-equilibrium degree exhibits an opposite trend. With Ar volume fraction increasing, the vibrational temperature first increases and then decreases, and the rotational temperature increases, while the thermodynamic non-equilibrium degree decreases.
-
0 引 言
粒子图像测速（Particle Image Velocimetry, PIV）技术[1]从平面二维（two-dimensional two-component, 2D−2C）测量发展到全体积三维（three-dimensional three-component, 3D−3C）测量，为探索流体的运动规律和特性提供了更强大和全面的工具[2-3]。3D−3C PIV技术可以提供瞬态三维速度场测量，能够更全面地揭示流体运动的细节和结构[4]。然而实现3D−3C测量的同时，往往需要增加相机的数量来获得更多的拍摄视角。实验系统的复杂性给受限空间内的三维流场测量带来了挑战。
面对这些挑战，研究者们积极探索各种解决方案。Gao等[5]在单个相机与被测流场之间加装一个三棱面特效透镜，光线通过该透镜三个棱面的折射后将相机感光芯片分为三个区域记录粒子的三个视角图像，实现多相机不同视角成像的效果。该技术极大简化了系统组成及实验操作，但该技术公共成像区域不足相机成像区域的1/3，导致有效测量区域受限。无法进行完整画幅测量。光场PIV[6]采用集成微透镜阵列的单台光场相机对体照明的示踪粒子进行成像获得粒子光场图片。但是微透镜阵列导致单个视角图像分辨率呈数量级下降，实现了完整画幅测量，但无法实现全分辨率测量，需要采用极高分辨率感光芯片，难以进行高频响流场测量。彩虹PIV[7]通过使用连续的色彩变化光平面组合成的体光源照亮示踪粒子，采用彩色相机记录被颜色信息编码的粒子深度信息。在该技术中光源照射方向与相机光轴方向夹角需保证一定的垂直角度，以实现颜色信息与深度信息较好的单一映射关系，但这一要求限制了其在有限光学窗口场景中的应用。
Zhang等[8]提出的三色掩膜PIV技术实现了单台相机对速度场进行3D−3C测量。其成像原理如图1所示，三色掩膜PIV采用三色掩膜通光孔调制彩色相机成像光路，将颜色信息与视角信息结合，实现感光芯片的RGB三通道分别对示踪粒子的三个视角图像进行完整画幅全分辨率记录，可以有效降低对相机分辨率的要求，提升相机采样频率。在此过程中，高质量的粒子三视角图像提取对提高速度场测量精度至关重要。
理想情况下，彩色相机每个像素处应该配备三个图像传感器，用于记录红、绿、蓝三种波段所传递的信息。然而，高昂的制造成本以及复杂的设计限制了三传感器彩色相机的普及和应用范围。目前大多数彩色相机使用单个图像传感器完成色彩信息的记录。单传感器相机通过在图像传感器上覆盖彩色滤波阵列（Color Filter Array，CFA）进行图像信息的采集。以最常见的Bayer CFA为例[9]，如图2所示。CFA阵列由一组光谱阵列排列而成，每个像素点位置只允许一种颜色分量通过。这样，每个像素点的值只能记录一种颜色分量，而缺失了另外两种颜色分量，从而生成原始马赛克图像。对马赛克图像中每个位置缺失的色彩信息进行重建的过程被称为图像的去马赛克。由于三色掩膜PIV技术将色彩信息与视角信息相结合，因此如何从粒子马赛克图像中重建高质量的全彩色图像是至关重要的。
Ramanath等[10]提出最近邻域法通过采用最近的像素点的采样值来插值缺失的像素值。以红色采样点为例，在2 × 2的领域内，复制采样点处红色像素的灰度值到其他三个位置缺失的红色灰度值，完成色彩的恢复。这种方法在图像边缘引入了明显的颜色误差，并不适用于高质量的颜色恢复需求。双线性插值[11]对最近邻域插值进行了改进，对图像3 × 3邻域上进行一个类似于平均的计算。这种插值方法在颜色变化缓慢的光滑区域表现良好，然而当在细节纹理或者颜色变化明显的边缘位置执行计算时，会引入拉链伪影。三次插值则是在一个7 × 7的邻域上计算，与双线性插值类似，采用两个卷积核，分别对缺失的信息进行估计计算。三次插值对伪影效果存在了一定的改善，但是成像质量还是不尽人意。Malvar等[12]提出了高质量线性插值（High Quality Linear Interpolation, HQLI），通过利用不同颜色通道之间的通道间相关性来扩展和改进线性插值。使用5 × 5邻域，对对应颜色通道的邻近像素取平均值，然后将其添加到根据不同颜色通道中的信息计算的校正项中。尽管与线性插值和三次插值相比，计算复杂性并没有增加多少，但该方法被证明优于许多更复杂的非线性方法，并且大大减少了边缘伪影。然而，HQLI算法没有考虑图像的边缘方向，导致在细节纹理和图像的高频区域可能出现一些错误。基于对细节纹理的考虑，Zhang等[13]提出了基于方向线性最小均方误差估计的去马赛克算法。该算法通过计算图像中相邻像素沿水平和垂直方向的颜色差异，使用线性最小均方误差框架将这些估计值进行最佳组合。Pekkucuksen等[14]提出了基于梯度的无阈值插值算法（Gradient-Based Threshold-Free Color Filter Array Interpolation，GBTF）对Zhang等人提出的算法进行了优化。他们将每个邻近像素都包括在给定的局部窗口内，以改进颜色差异估计。此外他们将南北方向和东西方向解耦，并将其分别考虑。这种方法利用图像的梯度信息来决定插值的方式，而不受阈值选择的影响。
基于深度学习的去马赛克方法可以从大量的数据集中提取图像的特征，并通过不断迭代来寻找原始图像和马赛克图像之间的映射关系，从而获得更好的去马赛克效果。受Bayer图案中绿色像素是红蓝像素之和的两倍这一事实的启发，Tan等[15]提出了一种基于CNN的两阶段去马赛克网络。第一阶段重建G通道的中间估计，第二阶段在重建的G通道的引导下重建R和B通道。该方法实现了图像的去马赛克，重建精度高，但网络结构复杂，参数多，运行速度慢，不适合实际应用。Kokkinos[16]提出了一种基于残差网络的联合去马赛克和去噪网络，其新颖的网络结构设计使其即使在较小的训练数据集上也能训练出良好的模型。虽然作者认为他们提出的网络结构的参数数量比其他性能更好的方法要少，但他们的网络中大量的跨层连接仍然使其运行时间不占优势，甚至比其他几个参数数量更多的网络运行时间更长。该方法需要首先对输入图像进行双线性插值，使其也不适合在边缘计算的实际应用。Wang等[17]提出的基于U-Net++神经网络的去马赛克算法通过增加密集连接的层块和使用深度可分离卷积，充分利用图像本身特征之间的相关性进行计算，大大减少了网络的参数数量，并仍取得了优异的性能。
大多数商业相机的去马赛克算法是保密的，去马赛克结果在细节上往往存在模糊、伪影等缺陷[18]。此外，这些算法大多数是基于自然场景环境开发的，并不适用于粒子图像的去马赛克。
鉴于此，本文采用了高质量线性插值算法、基于梯度的无阈值算法以及基于U-Net++神经网络的去马赛克算法，对Bayer格式为GRBG的粒子图像进行了去马赛克算法研究。通过人工合成高斯涡环三维流场分析三种算法对测量误差的影响。并通过对低雷诺数零质量射流的实验测量，探索了三种去马赛克算法在零质量射流实验中的实用性。
1 计算方法与评价
1.1 HQLI高质量线性插值算法
HQLI在双线性插值的基础上结合了不同通道之间的相关性，计算时增加了二阶拉普拉斯梯度校正。具体计算方法以红色位置$ R(i,j) $计算缺失的绿色信息$ \widehat g(i,j) $为例：
$$ \widehat g(i,j) = {\widehat g_B}(i,j) + \alpha {{\text{Δ}} _R}(i,j) $$ (1) 其中$ {{\text{Δ}} _R}(i,j) $为R 在该位置的梯度，由二阶拉普拉斯计算，具体公式为：
$$ \begin{gathered} {{\text{Δ}} _R}(i,j) \triangleq r(i,j) - \frac{1}{4}\sum {r(i + m,j + n)} \\ \left( {m,n} \right) = \left\{ {\left( {0, - 2} \right),\left( {0,2} \right),\left( { - 2,0} \right),\left( {2,0} \right)} \right\} \\ \end{gathered} $$ (2) 在公式（1）中$ {\widehat g_B}(i,j) $为该像素点位置对绿色的双线性插值计算结果，$ {{\text{Δ}} _R}(i,j) $为已知颜色R灰度值的梯度，$ \alpha $为增益因子控制梯度的校正量。
对于在绿色像素位置$ G(i,j) $处来插值缺失的红色分量$ \widehat r(i,j) $的公式和在蓝色像素位置$ B(i,j) $处来插值缺失的红色分量$ \widehat r(i,j) $的公式与公式（1）类似，不过在要使用周围已知信息来进行梯度校正计算。
增益参数的确定是通过Wiener方法来实现的。Wiener方法是一种信号处理技术，通常用于估计线性滤波器的系数，以最小化输出信号与目标信号之间的均方误差。作者使用了Kodak数据集计算了二阶统计量，使用1/2的整数次幂近似最优Wiener系数，最终确定增益参数的取值。
1.2 GBTF基于梯度的无阈值算法
GBTF算法主要通过插值绿色通道，计算色差图。对于红色像素位置，水平和垂直绿色通道的计算如下：
$$ \widetilde G_{i,j}^H = {{({G_{i,j - 1}} + {G_{i,j + 1}})} \mathord{\left/ {\vphantom {{({G_{i,j - 1}} + {G_{i,j + 1}})} 2}} \right. } 2} + {{(2 * {R_{i,j}} - {R_{i,j - 2}} - {R_{i,j + 2}})} \mathord{\left/ {\vphantom {{(2 * {R_{i,j}} - {R_{i,j - 2}} - {R_{i,j + 2}})} 4}} \right. } 4} $$ (3) $$ \widetilde G_{i,j}^V = {{({G_{i - 1,j}} + {G_{i + 1,j}})} \mathord{\left/ {\vphantom {{({G_{i - 1,j}} + {G_{i + 1,j}})} 2}} \right. } 2} + {{(2*{R_{i,j}} - {R_{i - 2,j}} - {R_{i + 2,j}})} \mathord{\left/ {\vphantom {{(2*{R_{i,j}} - {R_{i - 2,j}} - {R_{i + 2,j}})} 4}} \right. } 4} $$ (4) 同样对于绿色像素位置处，垂直方向和水平方向对红色像素和蓝色像素计算方法与公式（3）（4）类似。通过第一轮计算后，便可使用原始和方向估计的像素值求得水平和垂直方向的色差估计，公式如下:
$$ \begin{array}{c}{\tilde{{\text{Δ}} }}_{g,r}^{V}(i,j)=\left\{\begin{array}{c}{\tilde{G}}_{i,j}^{V}-{R}_{i,j},G是插值计算的\\ {G}_{i,j}-{\tilde{R}}_{i,j}^{V},R是插值计算的\end{array}\right.\\ {\tilde{{\text{Δ}} }}_{g,r}^{H}(i,j)=\left\{\begin{array}{c}{\tilde{G}}_{i,j}^{H}-{R}_{i,j},G是插值计算的\\ {G}_{i,j}-{\tilde{R}}_{i,j}^{H},R是插值计算的\end{array}\right.\end{array} $$ (5) 水平和垂直方向的G-B色差计算与公式（5）类似。在色差计算结束后，便可得出水平方向和垂直方向的色差图，如图3所示。
下一步结合方向色差，形成最终的目标像素值估计：
$$ \begin{aligned} {\widetilde {\text{Δ}} _{g,r}}(i,j) &= \left[ {{\omega _N}*f*\widetilde {\text{Δ}} _{g,r}^V(i - 4:i,j) } \right.\\ &\left.{+ \omega _S}*f*\widetilde {\text{Δ}} _{g,r}^V(i:i + 4,j) \right.\\ &\left.+ {\omega _E}*\widetilde {\text{Δ}} _{g,r}^H(i,j - 4:j)*f'\right. \\ &\left.{ + {\omega _W}*\widetilde {\text{Δ}} _{g,r}^H(i,j:j + 4)*f'} \right] * \frac{1}{{{\omega _T}}} \\ & {\omega _T} = {\omega _N} + {\omega _S} + {\omega _E} + {\omega _W} \\ & f = {{\left[ {\begin{array}{*{20}{c}} 1&1&1&1 \end{array}} \right]} \mathord{\left/ {\vphantom {{\left[ {\begin{array}{*{20}{c}} 1&1&1&1 \end{array}} \right]} 5}} \right. } 5} \\ \end{aligned} $$ (6) 每个方向$ ({w_N},{w_S},{w_E},{w_W}) $的权重是通过在局部窗口上添加该方向的色差梯度来计算的。对于5 × 5的窗口来讲：
$$ \begin{array}{l}{ \omega }_{N}= 1/{\left({\displaystyle \sum _{a=i-4}^{i}{\displaystyle \sum _{b=j-2}^{j + 2}{D}_{a,b}^{V}}}\right)}^{2} ,{\omega }_{S}= 1/{\left({\displaystyle \sum _{a=i}^{i + 4}{\displaystyle \sum _{b=j-2}^{j + 2}{ D}_{a,b}^{V}}}\right)}^{ 2}\\ { \omega }_{W}= 1/{\left({\displaystyle \sum _{a=i-2}^{i + 2}{\displaystyle \sum _{b=j-4}^{j}{D}_{a,b}^{H}}}\right)}^{2} ,{\omega }_{E}= 1/{\left({\displaystyle \sum _{a=i-2}^{i + 2}{\displaystyle \sum _{b=j}^{j + 4}{ D}_{a,b}^{H}}}\right)}^{ 2}\end{array} $$ (7) 其中梯度定义为：
$$ \begin{array}{*{20}{c}} {D_{i,j}^V = \left| {\widetilde {\text{Δ}} _{i - 1,j}^V - \widetilde {\text{Δ}} _{i + 1,j}^V} \right|} \\ {D_{i,j}^H = \left| {\widetilde {\text{Δ}} _{i,j - 1}^H - \widetilde {\text{Δ}} _{i,j + 1}^H} \right|} \end{array} $$ (8) 最后，对可用的（红色或蓝色）目标像素加上估计的色差，计算目标绿色的像素值：
$$ \begin{array}{*{20}{c}} {\widetilde G(i,j) = R(i,j) + {{\widetilde {\text{Δ}} }_{g,r}}(i,j)} \\ {\widetilde G(i,j) = B(i,j) + {{\widetilde {\text{Δ}} }_{g,b}}(i,j)} \end{array} $$ (9) 对于计算蓝色像素位置缺失的红色信息和红色像素位置缺失的蓝色信息是参考文献[19]所提出的插值计算方式。对于绿色位置的红色和蓝色像素，在最近的四个邻域上采用双线性插值。
1.3 基于U-Net++神经网络的去马赛克算法
本文采用的基于深度学习的去马赛克算法是Wang等[17]提出的紧凑高质量图像去马赛克神经网络，该神经网络是基于U-Net++结构设计的。网络结构由图像特征提取和图像重构两部分组成，如图4所示。白色框内部分为图像特征提取部分，它由每个级别的高斯平滑（红色箭头）和特征提取模块（绿色节点）组成。灰色框内部分为图像重建部分，它包含通道连接操作（白色圆圈）、重建节点（蓝色节点）和最后的上采样层（蓝色箭头）。
模型训练使用的损失函数（Loss Function）分为两部分。在初始的10次迭代中，使用图像的结构相似性指标（Structural Similarity Index Measurement, SSIM）作为网络的损失函数，以帮助网络快速学习图像的结构信息。经过这十次迭代（epoch）后，切换为使用图像的均方误差（Mean Squared Error，MSE）作为损失函数，以便网络逐渐学习图像的细节部分。由于U-Net++网络包含三层输出，因此在训练时需要考虑将L1、L2、L3输出的损失结合考虑作为最终的损失函数。因此，网络训练的最终损失函数表示为：
$$ L=\left\{\begin{array}{*{20}{c}}L\mathrm{^{SS\text{I}\text{M}_1}}+L^{\mathrm{SS\text{I}\text{M}_2}}+L^{\mathrm{SS\text{I}\text{M}_3}},epoch\leqslant10 \\ L^{\mathrm{M S E_1}}+L^{\mathrm{M S E_2}}+L^{\mathrm{M S E_3}},epoch > 10\end{array}\right. $$ (10) 使用U-Net++神经网络的优势在于当以这种深度监督的方式进行训练时，可以得到不同层次的输出，并且可以从整个模型中选择一个子网络，作为去马赛克计算成本和精度之间的权衡。此外，采用了高斯平滑层代替池化层来扩大输入拼接图像的接受域，并保持图像尺寸不变。此外，在每一层的开始处插入采用深度可分卷积的密集连接层单元，以充分提取模型参数较少的拼接图像的特征。在进入网络之前，输入的拼接图像被排列成RGGB的四通道图案，并根据相应层的网络级别通过高斯核分别模糊次数。通过结合来自不同尺度的特征图来增强模型的表达能力和性能，每一层的特征映射将与前一层和下层的特征映射连接起来，然后逐层通过重构节点，直到最顶层。最后，它们将作为L1、L2和L3来输出。
1.4 算法评价
本文采用彩色峰值信噪比（Color Peak Signal-to-Noise Ratio，CPSNR）和SSIM来评估去马赛克算法对粒子Bayer图像色彩重建的效果。
CPSNR的值越高表示图像的重构质量越好，即与原始图像的差异越小。CPSNR是对每个颜色通道的PSNR进行平均，公式如下：
$$ CPSNR = \frac{1}{3}(RSN{R_R} + RSN{R_G} + RSN{R_B}) $$ (11) 其中，PSNR的下标R、G、B分别代表的是彩色图像中红绿蓝三个颜色通道。单通道的PSNR定义为：
$$ PSNR = 10 \cdot {\log _{10}}(\frac{{M A{X^2}}}{{M SE}}) $$ (12) 其中MSE代表两张图像的均方误差，MAX代表图像中像素的最大值。对于常见的8位图像，MAX一般是255。
SSIM是一种用于衡量两幅图像之间结构相似性的指标，它从亮度、对比度和图像结构三个角度分析两张图像的相似性，能够更好地反映人眼对图像质量的感知。SSIM指数的取值范围为0到1，其中1表示两幅图像完全相同，0表示两幅图像完全不同。SSIM的公式如下：
$$ SSIM(x,y) = \frac{{(2{{\text{μ}} _x}{{\text{μ}}_y} + {C_1})(2{\sigma _{xy}} + {C_2})}}{{({\text{μ}} _x^2 + {\text{μ}} _y^2 + {C_1})(\sigma _x^2 + \sigma _y^2 + {C_2})}} $$ (13) 其中x、y是两个图像的对应块，$ {\mu _x} $、$ {\mu _y} $分别是x和y的平均值，$ {\sigma _x} $、$ {\sigma _y} $分别是x和y的方差，$ {\sigma _{xy}} $是x和y之间的协方差。C1和C2是是用于稳定公式的小常数。
2 仿真分析
2.1 粒子三视角图像分离评价
合成图像使用的三色掩膜相机分辨率为
1000 像素$ \times $1000 像素，像素尺寸为0.004 mm，三个掩膜孔中心外接圆半径为32 mm。利用光线追踪技术生成体校准图像，将三维示踪粒子投影到二维图像平面上，以获得三色掩膜的三个视角映射矩阵。随后，在测量体积内随机生成一组粒子，并通过相机映射矩阵计算每个视角下粒子的投影，以确定粒子中心的图像坐标。在这些坐标上应用高斯分布，生成直径约为3个像素的粒子图像。上述过程在每个视角上重复进行，得到三色掩膜全彩色粒子图像。本节的合成图像粒子浓度为0.05 ppp（particles per pixel，ppp）。在基于U-Net++神经网络的去马赛克算法训练时，选取了5000 张模拟的粒子图像作为网络输入，并进行了300次迭代，最终得到了去马赛克的网络模型。在合成三色掩膜全彩色粒子图像后，生成Bayer格式的马赛克图像，分别使用三种算法分别进行去马赛克处理。采用CPSNR和SSIM表征三种算法色彩重建质量，结果见表1。结果表明，U-Net++网络模型的去马赛克效果从CPSNR和SSIM的表现都优于传统去马赛克算法。图5所示的图像细节也表明基于U-Net++网络的去马赛克结果更好的还原了真实粒子所在的位置和对应像素的灰度值。
表 1 去马赛克质量Table 1 Demosaicing quality算法 CPSNR（dB） SSIM HQLI 17 0.74 GBTF 14 0.58 U-Net++ 31 0.95 此外，对去马赛克粒子图像进行了残差分析，包括对残差图像的均值和标准差的计算与分析。残差均值用于评估去噪处理的系统性偏差，而标准差则用于衡量残差的变异程度，从而进一步验证去马赛克效果的稳定性。表2展示了残差图像均值和标准差的计算结果。结果显示，U-Net++算法在去马赛克处理后的粒子图像中，残差均值和标准差均低于其他两种算法，表明U-Net++去马赛克算法在稳定性方面优于其他两种算法。
表 2 残差图均值与方差Table 2 The mean and variance of the residual map算法 残差均值 残差标准差 HQLI 6.8 16.1 GBTF 9.2 21 U-Net++ 3.3 7.8 为模拟真实实验场景中的粒子图像，本研究引入了不同强度的随机噪声，分别为10%、25%和40%，这些噪声强度以粒子峰值强度的百分比表示。对加噪数据集中的图像进行灰度值的平均计算，将其作为背景进行去噪处理。随后，对去噪后的图像进行了去马赛克处理。图6展示了去马赛克处理结果，计算的SSIM值见表3所示。结果表明，U-Net++算法在抗噪能力上展示出一定的优势。随着噪声强度的增加，其SSIM值虽有所下降，但仍优于HQLI和GBTF算法。
表 3 噪声对算法性能的影响Table 3 The impact of noise on algorithm performance算法 Noise（10%） Noise（25%） Noise（40%） HQLI 0.64 0.61 0.53 GBTF 0.46 0.44 0.38 U-Net + + 0.82 0.73 0.58 2.2 粒子三维重构质量影响
采用MLOS−SMART（Multiplicative Light of Sight-Simultaneous MART） [20]算法对提取的粒子三视角图像进行三维重构。获取三种算法提取的粒子三视角图像的三维重构质量Q值[21]。Q值表示粒子重构分布矩阵与真实分布矩阵的相似程度，如式（14）：
$$ Q = \dfrac{{\sum\limits_{X,Y,Z} {{E_1}(X,Y,Z) \cdot {E_0}(X,Y,Z)} }}{{\sqrt {\sum\limits_{X,Y,Z} {E_1^2(X,Y,Z)} \cdot \sum\limits_{X,Y,Z} {E_0^2(X,Y,Z)} } }} $$ (14) 其中，E0人工生成粒子场的分布矩阵，E1是MLOS−SMART重构后的粒子的分布矩阵。通过HQLI、GBTF、U-Net++网络模型分别对模拟的粒子图像进行色彩重建，并将重建后的粒子图像进行MLOS−SMART三维重构，计算粒子重构质量Q值。
Q值计算结果见图7，人工合成原始全彩色粒子场图像经MLOS−SMART三维重构后得到的重构质量Q值为0.52。经U-Net++去马赛克处理最终计算的Q值为0.45。结果表明U-Net++网络模型提取的三视角图像重构质量优于传统算法，并与全彩色三视角图像重构质量非常接近。
2.3 人工合成高斯涡环三维流场仿真实验研究
仿真采用人工合成高斯涡环速度场仿真实验进行速度场测量精度分析。人工合成高斯涡环流场的涡核为直径1.4 mm的圆，位于中心平面Z = 0 mm处。速度场d的表达式为：
$$ d = \left\| {\begin{array}{*{20}{c}} u \\ v \\ w \end{array}} \right\| = \dfrac{{0.03R}}{l}{e^{ - (\dfrac{R}{l})}} $$ (15) 式中：R是粒子与涡核的距离，l为涡宽度的尺度量（本文l = 0.3 mm）；u、v、w为速度矢量在3个方向的分量，w方向为相机景深方向。最大粒子位移为2.9 voxel。测量区域尺寸为（4$ \times $4$ \times $0.8）mm3。在测量区域内以0.05 ppp的示踪粒子浓度随机均匀播撒
50000 个粒子，并合成三色掩膜粒子图像。粒子三维速度场由多重网格三维互相关计算获得，第一重网格和第二重网格互相关窗口大小分别为64 $ \times $ 64 $ \times $ 64 voxel和32 $ \times $ 32 $ \times $ 32 voxel，重叠因子为50%。计算结果由3 × 3 × 3中值滤波器检测，异常值通过线性插值进行替换。图8左图为仿真真实涡环流场的速度矢量及涡量为0.065的等值面，右图为X = 0 mm处的Y–Z视角速度场和涡量场截面图。
图9（a）为真实图像经计算所得涡环流场的速度矢量及涡量图。图9（b）和（c）展示了不同质量的图像输入测量的人工合成流场结果。基于U-Net++神经网络的去马赛克算法基本可以准确还原涡环流场结构；传统的算法计算结果出现较大的偏差。
![]()
图 9 不同质量图片人工合成流场测量结果图Fig. 9 Measurement results of artificially synthesized flow field for images of different quality在不同质量的插值图像输入下，速度场u、v、w测量误差的累积分布函数（Cumulative Distribution Function, CDF）曲线如图10所示。由不同质量图像计算后，u和v的误差累积分布函数几乎完全重合。但是基于HQLI和GBTF插值算法的图像出现的误差明显偏大，在w方向尤其明显，其中80%矢量误差绝对值分别对应于1.5 voxels、1.8 voxels。基于U-Net++算法结果的误差控制效果明显优于传统插值图像，在w方向80%矢量误差绝对值为0.4 voxels。误差散点图如图11所示。其中u和v的误差以零误差点为中心聚集，正负两边整体对称。由于粒子重构会在景深方向拉长，低质量的图像会在景深方向的重构出现较多错误。但是从CDF曲线和误差散点图来看，基于U-Net++的去马赛克算法结果和真实图像输入的结果十分接近，在对粒子图像的去马赛克表现优异。
![]()
图 10 u、v、w测量误差的累积分布函数Fig. 10 Cumulative distribution function of u、v、w measurement errors3 零质量射流实验分析
上节利用人工合成高斯涡环流场研究了不同去马赛克算法对最终测量结果的影响，结果表明基于U-Net++神经网络的去马赛克算法测量误差表现优于传统算法。考虑到实验真实拍摄存在环境噪声、光照亮度不足等因素，基于实验图像去马赛克具有一定的挑战。本节利用零质量射流实验验证去马赛克算法在三色掩膜PIV应用中的实用性。
实验中采用了三波段白光激光器，其光束功率为3 W，波长分别为450 nm、532 nm和650 nm。将定制的三色掩膜安装在ZEISS Milvus 2/100M镜头前，并与Revealer M120彩色相机连接，搭建了三色掩膜成像系统，如图12所示。相机的分辨率为
1024 像素 $ \times $1280 像素，像元尺寸为10 μm × 10 μm。实验装置如图13所示，水箱由15 mm厚的亚克力板制成，尺寸为
1000 mm长、500 mm宽和500 mm高。活塞的半行程设定为13 mm，驱动频率f为0.25 Hz。实验的斯特劳哈尔数为0.77，雷诺数为32。射流圆形喷嘴直径D = 10 mm。测量区域位于喷嘴左侧0.1D处，测量体积定义为沿x、y、z方向分布的2D $ \times $ 2D $ \times $ 0.4D。采用聚酰胺示踪粒子，直径为20 µm，密度为1.03 g/cm³。粒子的浓度大约为0.04 ppp，实验拍摄相机设定为150 帧/s。记录实验粒子图像后，传统算法选用表现较好的高质量线性插值算法，与基于U-Net++神经网络的去马赛克算法进行比较，主要分析对实验结果的处理效果，观察射流的速度分布状态以及Z方向的涡量。三视角图像提取结果如图14所示，图14（a）中粒子像素大小大部分控制在4像素 × 4像素范围。与图14（b）相比，粒子像素的灰度值控制较好，没有出现太多过曝的像素点。
![]()
图 14 实验粒子图片去马赛克以及三视角分离图Fig. 14 Demosaicing of experimental particle images and separation of three perspectives在MLOS−SMART三维重构过程中，测量区域被离散化为
1000 $ \times $1000 $ \times $ 200个体素，体素分辨率为−0.004 mm $ \times $ 0.004 mm $ \times $ 0.004 mm/M3（其中M为镜头的放大倍率，取−0.5），松弛因子μ设置为2.0，并进行了100次迭代。瞬时速度场通过多网格三维互相关算法计算得到，第一层网格的窗口大小为64 $ \times $ 64 $ \times $ 64体素，第二层网格的窗口大小为32 $ \times $ 32 $ \times $ 32体素，重叠率为50%。重构结果经过3 $ \times $ 3 $ \times $ 3的中值滤波处理，并对异常值进行了线性插值替换。图15所示为两种算法最终计算的瞬态速度场图。在图中箭头表示速度矢量，每个箭头的长度与速度的大小成正比，箭头的颜色代表速度值分布。等值面根据Z方向的涡量进行绘制。类似于人工合成高斯涡环速度场的测量结果，经HQLI算法计算结果的涡量结构不完整。参考Zhong等[22]的零质量射流测量结果，U-Net++神经网络去马赛克得到的结果计算的涡量更为完整，速度分布也基本正确。
对于不可压缩流体，其散度为零（$\nabla $∙u=0），即$ {{\partial u}/{\partial x }+{\partial v}/{\partial y} + {\partial w}/ {\partial z = 0}}$。散度误差通过${\partial u}/{\partial x } $与$-({\partial v}/{\partial y} + {\partial w}/ {\partial z}) $的联合概率密度函数（JPDF）表示，其中数据偏离零散度线（即各分量相等的线）的程度反映了散度的大小。散度通过分量间的Pearson相关系数$ Q_{\mathrm{div}} $进行量化[23]。计算的联合概率密度如图16所示，经U-Net++提取的三视角图像最终计算的$ Q_{\mathrm{div}} $值为0.77，而通过HQLI提取的三视角图像所获得的$ Q\mathrm{_{div}} $值为0.63。结果表明，经U-Net++计算提取的三视角图像更加接近实际的三视角图像。此外Worth等[23]通过层析PIV结果计算得到的$ Q\mathrm{_{div}} $为0.66，Ganapathisubramani等[24]通过Stereo-PIV结果计算得到的$ Q\mathrm{_{div}} $为0.82。三色掩模PIV计算得到的$ Q\mathrm{_{div}} $略大于层析PIV计算得到的$ Q\mathrm{_{div}} $，这可能是由于测量流场的不同三维流动特性所致。
![]()
图 16 ${\partial u}/{\partial x } $与$-({\partial v}/{\partial y} + {\partial w}/ {\partial z}) $的联合概率密度函数图Fig. 16 JPDF between ${\partial u}/{\partial x } $ and $-({\partial v}/{\partial y} + {\partial w}/ {\partial z}) $4 结 论
本文采用HQLI、GBTF、基于U-Net++神经网络方法分别对三色掩膜PIV粒子Bayer图像进行了去马赛克研究，完成了粒子图像三视角提取。采用CPSNR、SSIM评价了三视角图像的提取质量；计算了不同算法提取的三视角图像的粒子重构质量Q值；通过对人工合成高斯涡环三维流场进行仿真实验，分析了三种算法对测量误差的影响；利用零质量射流三维测量实验研究真实实验中粒子图像三视角提取，分析了瞬态速度场的结果。主要结论如下：
1）使用基于U-Net++神经网络的去马赛克算法可以对粒子Bayer图像进行去马赛克处理，能够完成三色掩膜PIV技术中的三视角提取。其中，基于U-Net++神经网络提取的三视角图像的CPSNR和SSIM均优于传统方法；粒子三维重构质量Q也与全彩色图像三维重构质量Q值接近。
2）从人工合成高斯涡环三维流场来看，三视角图像的提取质量对测量误差影响较大。基于U-Net++的去马赛克算法提取的三视角图像与真实三视角图像接近，重建三维流场的涡结构较为完整。测量速度误差最小，其中在z方向，80%的速度误差绝对值为0.4 voxels。
3）通过低雷诺数零质量射流真实实验测量研究，相较于传统算法，基于U-Net++的去马赛克算法提取的三视角图像粒子像素大小和像素灰度值分布更符合粒子成像规律，能够更准确还原粒子三视角图像，测量的流场结果与Zhong等[22]的零质量射流测量结果更加相符。
-
![]()
图 2 CARS光谱计算和拟合过程示意图
Fig. 2 Flow diagram of CARS spectrum calculation and fitting program
![]()
图 3 （a） 非平衡等离子体振转温度测量实验；（b）USED相位匹配方式下各个光束的空间位置示意图；（c） CARS测量系统示意图
Fig. 3 (a) Experiment on the measurement of vibrational-rotational temperature of non-equilibrium plasma; (b) Schematic diagram of the spatial positions of various beams in the USED phase matching method; (c) Schematic diagram of the CARS measurement system
![]()
图 4 （a）实测CARS光谱及拟合结果；（b）宽范围CARS测温结果
Fig. 4 (a) Measured CARS spectrum and corresponding fitted results, (b) measurement results of vibrational and rotational temperature in wide-range
![]()
图 5 不同压强下的（a）实测CARS光谱及拟合结果、（b）振转温度及非平衡度
Fig. 5 (a) The measured CARS spectra and corresponding fitted results, and (b) vibrational-rotational temperatures and thermal non-equilibrium with different pressure
![]()
图 6 不同N2流量下的（a）实测CARS光谱及拟合结果、（b）振转温度及非平衡度
Fig. 6 (a) The measured CARS spectra and corresponding fitted results, and (b) vibrational-rotational temperatures and thermal non-equilibrium with different N2 flow rates
![]()
图 7 不同Ar流量下的（a）实测CARS光谱及拟合结果、（b）振转温度及（c）非平衡度
Fig. 7 (a) The measured CARS spectra and corresponding fitted results, (b) vibrational-rotational temperatures, and (c) thermal non-equilibrium with different Ar flow rates
表 1 微波等离子体参数
Table 1 Parameters of microwave plasma
工况 频率f/MHz 功率P/W 压强p/kPa 流量v/（L·min−1） 实验气体 1 2450 160 1～5 2 N2 2 2450 160 3 0.5～3.0 N2 3 2450 160 3 2 N2、Ar 
下载: 导出CSV
-
[1] 姜宗林, 李进平, 胡宗民, 等. 高超声速飞行复现风洞理论与方法[J]. 力学学报, 2018, 50(6): 1283–1291. DOI: 10.6052/0459-1879-18-238 JIANG Z L, LI J P, HU Z M, et al. Shock tunnel theory and methods for duplicating hypersonic flight conditions[J]. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(6): 1283–1291. doi: 10.6052/0459-1879-18-238
[2] KNISELY C P, ZHONG X L. Impact of vibrational nonequilibrium on the supersonic mode in hypersonic boundary layers[J]. AIAA Journal, 2020, 58(4): 1704–1714. doi: 10.2514/1.j058758
[3] CHEN X L, WANG L, FU S. Secondary instability of the hypersonic high-enthalpy boundary layers with thermal-chemical nonequilibrium effects[J]. Physics of Fluids, 2021, 33(3): 034132. doi: 10.1063/5.0045184
[4] LI W J, HUANG J, ZHANG Z W, et al. A model for thermal protection ablative material with local thermal non-equilibrium and thermal radiation mechanisms[J]. Acta Astronautica, 2021, 183: 101–111. doi: 10.1016/j.actaastro.2021.03.001
[5] PARK C. Hypersonic aerothermodynamics: past, present and future[J]. International Journal of Aeronautical and Space Sciences, 2013, 14(1): 1–10. doi: 10.5139/IJASS.2013.14.1.1
[6] HOLDEN M, WADHAMS T, MACLEAN M, et al. Review of basic research and development programs conducted in the LENS facilities in hypersonic flows[C]//Proceedings of the 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition. 2012. doi: 10.2514/6.2012-469
[7] PARK C. Stagnation-region heating environment of the Galileo probe[J]. Journal of Thermophysics and Heat Transfer, 2009, 23(3): 417–424. doi: 10.2514/1.38712
[8] MATSUYAMA S, OHNISHI N, SASOH A, et al. Numerical simulation of Galileo probe entry flowfield with radiation and ablation[J]. Journal of Thermophysics and Heat Transfer, 2005, 19(1): 28–35. doi: 10.2514/1.10264
[9] GRISCH F, BOUCHARDY P, JOLY V, et al. Coherent anti-stokes Raman scattering measurements and computational modeling of nonequilibrium flow[J]. AIAA Journal, 2000, 38(9): 1669–1675. doi: 10.2514/2.1152
[10] KIM M K, ESSER B, KOCH U, et al. Numerical and experimental study of high enthalpy flows in a hypersonic plasma wind tunnel: L3K[C]//Proc of the Proceedings of the 42nd AIAA Thermophysics Conference. 2011. doi: 10.2514/6.2011-3777
[11] GÜLHAN A, ESSER B, KOCH U, et al. Characterization of high-enthalpy-flow environment for ablation material tests using advanced diagnostics[J]. AIAA Journal, 2018, 56(3): 1072–1084. doi: 10.2514/1.j056312
[12] SAKURAI K, OSADA T, NIINOMI S, et al. High-speed visualization of total radiation and CARS measurement of vibrational/rotational temperatures behind hypervelocity shock waves of 5 km/s[J]. Journal of JSEM, 2010, 10: 85–89. doi: 10.11395/jjsem.10.s85
[13] VENIGALLA H B, NIINOMI S, OGURO M, et al. The spatial distribution of ro/vibrational temperatures for nitrogen molecules behind hypervelocity shock waves by CARS measurement[J]. Journal of JSEM, 2012, 12: 46–50. doi: 10.11395/jjsem.12.s46
[14] DOGARIU A, DOGARIU L E, SMITH M S, et al. Single Shot Temperature Measurements using Coherent Anti-Stokes Raman Scattering in Mach 14 Flow at the Hypervelocity AEDC Tunnel 9[C]//Proceedings of the AIAA Scitech 2019 Forum. 2019. doi: 10.2514/6.2019-1089
[15] DOGARIU A, DOGARIU L E, SMITH M S, et al. Velocity and temperature measurements in Mach 18 nitrogen flow at tunnel 9[C]//Proceedings of the AIAA Scitech 2021 Forum. 2021. doi: 10.2514/6.2021-0020
[16] MONTELLO A, NISHIHARA M, RICH J W, et al. Picosecond USED-CARS for simultaneous rotational/translational and vibrational temperature measurement of nitrogen in a nonequilibrium Mach 5 flow[C]//Proceedings of the 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition. 2012. doi: 10.2514/6.2012-239
[17] 杨文斌, 齐新华, 李猛, 等. 超燃冲压发动机燃烧室出口温度场分布CARS测量[J]. 推进技术, 2022, 43(9): 210190. DOI: 10.13675/j.cnki.tjjs.210190 YANG W B, QI X H, LI M, et al. Outlet temperature measurements of scramjet combustor using Coherent Anti-Stokes Raman scattering[J]. Journal of Propulsion Technology, 2022, 43(9): 210190. doi: 10.13675/j.cnki.tjjs.210190
[18] 杨文斌, 齐新华, 王林森, 等. 基于CARS技术的超燃冲压发动机点火过程温度测量[J]. 气体物理, 2020, 5(2): 8–13. DOI: 10.19527/j.cnki.2096-1642.0800 YANG W B, QI X H, WANG L S, et al. Temperature measurement based on CARS in scramjet ignition process[J]. Physics of Gases, 2020, 5(2): 8–13. doi: 10.19527/j.cnki.2096-1642.0800
[19] HUANG A, XU Z Y, DENG H, et al. High-pressure gas temperature sensing for exit plane of aero-engine combustor using tunable diode laser absorption spectroscopy[J]. Microwave and Optical Technology Letters, 2024, 66(1): 1–9. doi: 10.1002/mop.33897
[20] HU Z Y, LIU J R, YE J F, et al. Laser-based measurements of temperature, species, and velocity in engine combustor[C]//Proc of the 2nd International Symposium on Laser Interaction with Matter. 2013. doi: 10.1117/12.2011241
[21] 张虎. 基于CARS的火焰温度测量技术研究[D]. 哈尔滨: 哈尔滨工业大学, 2009. ZHANG H. Study on flame temperature measurement based on CARS[D]. Harbin: Harbin Institute of Technology, 2009.
[22] DEDIC C E. Hybrid fs/ps coherent anti-Stokes Raman scattering for multiparameter measurements of combustion and nonequilibrium[D]. Iowa State University, 2017.
[23] 赵松原. 模拟退火结合正交分解算法的气动外形最优化设计[D]. 南京, 南京航空航天大学, 2006. ZHAO S Y. Research on optimum aerodynamic design using simulated annealing algorithm and POD[D]. NanJing: Nanjing University of Aeronautics and Astronautics, 2006.
[24] WEIGAND P, LÜCKERATH R, MEIER W. Documentation of flat premixed laminar CH4/air standard flames: temperatures and species concentrations[J/OL]. http//www. dlr. de/VT/Datenarchiv, 2003.
计量
- 文章访问数: 49
- HTML全文浏览量: 11
- PDF下载量: 2
