一种基于变分光流模型的流体速度测量方法

技术领域

本发明属于图像处理相关技术领域，更具体地，涉及一种基于变分光流模型的流体速度测量方法。

背景技术

随着测量技术的发展，流体流动测量技术逐渐取代了流动可视化技术在实验流体动力学研究中的地位。通过测量技术提取流体速度场，可以更深入的分析复杂的流场结构。粒子图像测速技术(Particle Image Velocimetry，PIV)作为一种非接触、全场、瞬时测量技术，已成为流体流动测量的主要技术手段之一，并被广泛应用于航空航天、生物医学和纳米制造等领域。PIV技术测量的原理是将示踪粒子播撒到流体中以跟随流体运动，使用预先设定间隔时间的两束脉冲激光照明流体中某一平面内的示踪粒子，同时利用相机同步捕获两帧粒子图像。最后，通过确定粒子在两帧图像之间的位移和预先预定的间隔时间来计算流体速度场。

从PIV图像中提取粒子场运动速度的测量方法有两类：互相关方法和变分光流方法。互相关方法将粒子图像按照查询窗口的尺寸划分为网格，询问窗口是图像中的子域，通常大小在64×64pixels和16×16pixels之间。通过计算连续两帧图像查询窗口的最大相关系数来确定每个窗口的速度矢量，这个速度矢量代表着查询窗口内所有粒子的平均速度。互相关测量方法得到的速度场分辨率直接依赖于查询窗口的尺寸，通常要比PIV图像的分辨率低一个数量级。同时由于每一个速度矢量都是窗口内的平均速度，在流体速度梯度较大的区域，互相关方法的测量误差就会较大。变分光流方法假设示踪粒子亮度在连续两帧图像中保持不变的(亮度一致性假设)，建立了变分光流模型中的数据项。同时假设流场速度场具有连续平滑的性质(平滑约束假设)，并以此建立了变分光流模型中的平滑项。通过对建立的光流模型进行变分求解，可以得到与图像分辨率一样的速度场，即每像素一个速度矢量。变分光流模型通常是基于计算机视觉理论，并不是从流体力学基本方程中推导出来的。针对复杂流动结构，比如激波、回流区、剪切层、强涡流等，亮度一致性假设和平滑约束假设并不是严格成立的，变分光流方法的测量误差就会较大。同时变分光流方法对PIV实验中出现的图像采集噪声和激光光照不均匀的现象非常敏感，测量鲁棒性低。因此，亟需设计一种适用于高速复杂流体的流体速度测量方法。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种基于变分光流模型的流体速度测量方法，避免了求解过程中对亮度的不合理假设，更加符合实际运动过程，误差更小，同时针对图像采集噪声和激光光照影响具有较好的鲁棒性。

为实现上述目的，按照本发明的一个方面，提供了一种基于变分光流模型的流体速度测量方法，所述方法包括：S1，采用质量守恒方程替代变分光流模型的数据项中的亮度平衡方程，并将所述质量守恒方程中的密度参数替换为亮度参数，而后对该数据项进行噪声修正，同时，将流体的速度分解为多种形式的运动分量，并对分解的每项分别进行平滑约束和噪声修正后求和获得所述变分光流模型的平滑项，以此方式得到优化后的变分光流模型；S2，获取采集的具有示踪粒子的连续两帧粒子图像的亮度参数，将所述亮度参数代入优化后的所述变分光流模型即可求得所述流体的速度。

优选地，步骤S2具体包括：S21，对采集的连续两帧粒子图像进行降采样获得对应层中分辨率大小一致的第一帧粒子图像金字塔和第二帧粒子图像金字塔；S22，根据所述第一帧粒子图像和第二帧粒子图像在该层的亮度获取亮度梯度，将所述亮度梯度代入所述变分光流模型获得速度场(u

优选地，步骤S1中所述对该数据项进行噪声修正具体为采用双边滤波器进行修正；对分解的每项进行噪声修正具体为采用双边滤波器或基于图像亮度信息的函数进行修正。

优选地，步骤S1中修正后的数据项E

其中，Ω为整个二维图形空间，x和y分别为像素坐标，I为输入的粒子图像，u(x)＝(u(x)，v(x))

优选地，所述将流体的速度分解为多种形式的运动分量具体为采用亥姆霍兹速度分解定理将流体的速度分解为将所述流体的运动分为平移运动分量u、线变形运动分量χ、剪切变形运动分量θ和旋转运动分量ζ。

优选地，所述步骤S1中所述并对分解的每项进行平滑约束的具体表达式为：

其中，

优选地，修正后的所述平滑项E

其中，BF

优选地，步骤S1还包括对优化后的所述变分光流模型引入中间变量而后对其进行变分求解获得速度场增量方程和中间变量方程。

优选地，步骤S22中求解速度场增量时具体步骤为：S221，根据所述速度场增量方程获取速度场增量；S222，将所述速度场增量代入所述中间变量方程获得中间变量；S223，将步骤S222获得的中间变量代入步骤S221中所述速度场增量方程重新计算速度场增量和中间变量，以此方式不断重复步骤S221和S222，直至迭代次数达到预设值则输出的速度场增量为最终的速度场增量。

优选地，所述速度场增量方程表达式为：

其中，

所述中间变量方程的表达式为：

其中，

总体而言，通过本发明所构思的以上技术方案与现有技术相比，本发明提供的一种基于变分光流模型的流体速度测量方法至少具有如下有益效果：

1.采用质量守恒方程替换变分光流模型的数据项中的亮度平衡方程，由于浓度与亮度成正比，因此将质量守恒方程中的密度参数可以替换为亮度参数，进而可以避免在求解过程中使用并不是严格成立的亮度一致性假设，降低了求解误差。

2.数据项基于流体力学基本方程推出而不是采用基于计算机视觉理论的亮度平衡方程推出，可以适用于复杂流动结构，比如激波、回流区、强涡流的流动情况，极大的扩展了应用范围。

3.通过将速度方程分解为多种形式的运动分量，而不仅仅是相互垂直的两个速度方向，使得速度分解的多种分量均进行了平滑处理，得到的速度场误差更小。

4.对数据项和平衡项均进行了噪声修正，对图像采集噪声和激光光照隐形具有较好的鲁棒性。

5.将每帧图像均处理为图像金字塔，使得可以逐层迭代计算每层中的速度场，直至最终得到与原图像分辨率一致的图像，整个过程没有进行速度场的简化和假设，进而获得速度场可以准确表征真实的速度场。

6.在对变分光流模型进行求解时，可以采用引用中间变量的形式，将变分光流模型分解为速度场增量方程和中间变量方程，极大的简化了运算。

7.采用速度场增量方程求取速度场增量时，不断的对中间变量进行迭代修正，并将修正后的中间变量代入到速度场增量方程中进行速度场增量方程的迭代计算，进而可以变分光流模型的最小化求解，得到最优的速度场增量。

8.采用亥姆霍兹速度分解定理将流体的速度进行分解，并对分解的每项进行平滑约束可以保护流程中的非均匀流动结构，可以适用于湍流等复杂流场和激波的高精度高分辨率测量。

附图说明

图1示意性示出了本实施例的基于变分光流模型的流体速度测量方法的流程图；

图2示意性示出了本实施例的单层图像金字塔内的图像变形操作的流程图；

图3A示意性示出了本实施例的第一帧粒子图像；

图3B示意性示出了本实施例的第二帧粒子图像；

图4示意性示出了本实施例的第一帧粒子图像和第二帧粒子图像构建的图像金字塔。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

请参阅图1，本发明提供了一种基于变分光流模型的流体速度测量方法，所述方法包括步骤S1～S2。

S1，采用质量守恒方程替代变分光流模型的数据项中的亮度平衡方程，并将所述质量守恒方程中的密度参数替换为亮度参数，而后对该数据项进行噪声修正，同时，将流体的速度分解为多种形式的运动分量，并对分解的每项分别进行平滑约束和噪声修正后求和获得所述变分光流模型的平滑项，以此方式得到优化后的变分光流模型。

现有技术中变分光流模型包括数据项和平滑项，其表达式为：

其中，

本实施例首先对该变分光流模型进行优化，包括对数据项的优化和对平滑项的优化。

基于流体力学中的质量守恒方程对数据项进行优化，由于PIV技术中密度与亮度成正比，因此可以将质量守恒方程中的密度参数全部替换为亮度参数。具体过程如下。

S11，求解数据项

在二维空间平面中，针对某一位置x＝(x，y)

在PIV测量中，示踪粒子对激光进行散射，散射光在相机中成像粒子图像。粒子图像的亮度是与示踪粒子的浓度呈正比的，即图像的亮度I正比于流体密度ρ。在此基础上添加双边滤波器BF

其中，Ω为整个二维图形空间，x和y分别为像素坐标，I为输入的粒子图像，u(x)＝(u(x)，v(x))

S12，求解平滑项

在二维平面中，针对某一位置x＝(x，y)

u(x+δx)＝u(x)+χ·δx+θ·δx+ζ·δx

其中，χ为线性变形张量，θ为剪切变形张量，ζ为旋转张量，

其中，

在此基础上添加基于图像亮度梯度信息的单调递减函数

其中，BF

对上述数据项平滑项求和获得变分光流模型E’，如下：

其中λ为平滑参数。可以对以上公式进行直接迭代求解，也可以如下所示引入中间变量

通过E

对E

其中，

对于

对

同样地对

通过对

S2，获取采集的具有示踪粒子的连续两帧粒子图像的亮度参数，将所述亮度参数代入优化后的所述变分光流模型即可求得所述流体的速度。

可以采用PIV技术获取连续两帧粒子图像的亮度参数，将该亮度参数代入上述优化后的变分光流模型E’或者公式(1)～(6)进行最小化求解即可获得流体的速度。

步骤S2优选采用多尺度迭代的方式对速度场进行求解，包括以下步骤S21～S23具体为粒子图像的多尺度迭代操作，包括图像金字塔处理和图像变形操作，具体如图2所示，具体如下。

S21，对采集的连续两帧粒子图像进行降采样获得对应层中分辨率大小一致的第一帧粒子图像金字塔和第二帧粒子图像金字塔。

对PIV采集的连续两帧粒子图像进行降采样获得对应层中分辨率大小一致的第一帧粒子图像金字塔和第二帧粒子图像金字塔，从顶层到底层分别为第1层，第2层……第K层。

通过对PIV采集的连续两帧粒子图像(如图3A和图3B所示)进行一系列的低通滤波和下采样处理，构建粒子图像的图像金字塔，金字塔中的图像分辨率由下向上依次降低，金字塔最底层(第K层)为原始粒子图像。并且第一帧图像金字塔和第二帧图像金字塔的每一层相互对应，如图4所示。

S22，根据所述第一帧粒子图像和第二帧粒子图像在该层的亮度获取亮度梯度，将所述亮度梯度代入所述变分光流模型获得速度场(u

S22a，从最顶层的第1层金字塔进行计算，首先初始化速度场u和v，初始化迭代参数n、m和k，其中，n为步骤S1中求解速度增量时的迭代次数，m为每层金字塔层中图像变形操作的迭代次数，k为金字塔中的层数。其中，u和v初始化为零矩阵，其规模大小与图像金字塔顶层的分辨率一致，迭代参数n、m和k初始化为1。

S22b，通过迭代求解(δu，δv)和

S221，根据所述速度场增量方程获取速度场增量；

对公式(1)和(2)中的偏微分方程进行离散化和迭代求解，获得速度场(u，ν)，其中该速度场(u，v)即为速度场增量(δu，δv)。

S222，将所述速度场增量代入所述中间变量方程获得中间变量；

根据步骤S221求解的速度场增量(δu，δv)对公式(3)～(6)进行求解获得中间变量

S223，将步骤S222获得的中间变量代入步骤S221中所述速度场增量方程重新计算速度场增量和中间变量，以此方式不断重复步骤S221和S222，直至迭代次数达到预设值则输出的速度场增量为最终的速度场增量。

将该中间变量

S22c，将步骤S22b计算的速度场增量(δu，δv)添加至速度场(u，v)，得到更新后的速度场，即u＝u+δu，v＝v+δv。将更新后的速度场对当前金字塔层中的第二帧粒子图像进行图像变形操作，图像变形的插值方法使用双三次插值，迭代参数m＝m+1，并返回步骤S22b，以此方式不断的迭代，直至迭代次数m的值达到预设次数M，则此时的速度场为本层的最终速度场，并将该速度场作为下一层图像的初始速度场。

S23，对所述下一层图像重复执行步骤S22直至第K层图像执行结束所得速度场即为流体速度场。

以上一层计算的速度场为初始速度场重新执行步骤S22中的操作，以此方式不断迭代，直至第K层图像执行结束，所得速度场即为流体速度场。

综上所述，本发明所提出的方法，建立了基于多种形式速度分解的光流模型平滑项，为相应的平滑项设计自适应的平滑项参数，保护流场中的非均匀流动结构；同时建立了基于流体质量守恒的光流模型数据项，使用双边滤波器作为数据项参数消除粒子图像噪声和光照不均对测量的影响。该方法综合考虑了粒子图像信息和流体物理理论，不仅适用于一般低速不可压缩湍流，也适用于高超声速可压缩流场的精确测量，同时针对粒子图像的噪声和光照影响具有较好的鲁棒性。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。