技术领域
[0001] 本
发明属于
流体测量领域，更具体地，涉及一种空间分辨率自适应调整的粒子图像测速矢量估计方法。
背景技术
[0002] 粒子图像测速(Particle Image Velocimetry简称PIV)是一种用多次摄像记录流场中粒子的
位置，分析摄得图像，从而测出流动速度的方法。测量原理是在流体中添加能够与之共同移动的示踪粒子，使用激光片
光源将之照亮，使得看不见的流体
可视化为离散的亮点，通过相机记录图像并按照一定的规则识别出在同一幅图像或两幅图像上对应的粒子或粒子集合对，得到粒子或粒子集合的运动速度，并以此作为粒子或粒子集合所在位置的流场的速度。如今，PIV作为一种全场非
接触无扰动测量方法广泛应用于实验
流体力学、
生物医学、航空航天、工业制造等诸多领域。
[0003] 传统PIV相关计算方法是从独立存在的两幅图像通过一定的判别规则。例如，互相关计算，首先建立离散的查询窗口，并假设窗口内所有示踪粒子以相同的速度做刚性运动，通过计算粒子集合的互相关系数，得到概率最佳配准位移，然后用该位移表征该窗口内所有粒子的速度矢量，最终得到流场中各区域的流速矢量。
[0004] 迄今，人工定义查询窗口的方法，当速度矢量场为非定常流动的时候将会失效。然而实事上定常流动只是一个理想化的状态，任何流动都只可能是在一定范围内的近似定常流动，没有绝对的定常流动。因为流动结构特点在被测区域上的空间位置会随着时间变化移动，所以如果采用人工设置的方法，PIV用户仍然需要具有相当丰富的知识和经验来人为选择最优化的数据记录参数和恰当的数据分析处理参数，特别是相关运算的查询窗口尺寸的设置。在设置这一参数时，须在鲁棒性，
精度和有效空间分辨率等“目标”权衡。此外由于所观测流场各局部区域之间的示踪粒子
密度，图像
质量和流动条件具有不可避免的差异和变化，
对流场采用全局统一单一的分析参数设置将永远无法实现对整个流场的各部分均有效的最优化。在高超声速测量时，查询窗口参数设置问题显得尤为重要。
[0005] 文献号为CN101629966B的中国
专利文献，公开了一种粒子图像测速处理方法，其中公开了利用构建粒子匹配
水平集函数和最小化
能量函数的方法，实现PIV图像矢量估计的处理方法。然而，该方法仍然存在以下
缺陷：首先，该方法虽然在平滑约束项和基本约束项中了做了优化改进，该方法仍是一个单向性的计算过程，对于不同流场，其精度受计算过程参数的影响较大，鲁棒性不够；其次，对于高速流场，由于速度变化剧烈，该方法无法
算法做出自适应调整，计算精度会降低。
发明内容
[0006] 针对
现有技术的以上缺陷或改进需求，本发明提供了一种空间分辨率自适应调整的粒子图像测速矢量估计方法，其目的在于通过在线调整查询窗口参数，由此解决现有技术查询窗口参数设置不能随流场变化而变化，造成现有粒子图像测数矢量估计方法鲁棒性不够或精度较低的技术问题。
[0007] 为实现上述目的，按照本发明的一个方面，提供了一种空间分辨率自适应调整的粒子图像测速矢量估计方法，包括以下步骤：
[0008] （1）对第n时刻和第n+1时刻的图像进行降噪处理；
[0009] （2）根据第n时刻的查询窗口参数更新第n+1时刻的查询窗口参数；
[0010] （3）根据第n+1时刻的查询窗口参数，根据第n+1时刻的查询窗口参数，结合第n和n+1时刻图像进行n到n+1时刻间的流场矢量估计，得到全场速度矢量；
[0011] 所述第n时刻的查询窗口参数，当n为1时，查询窗口参数由人工
指定。
[0012] 优选地，所述的粒子图像测速矢量估计方法，其查询窗口参数包括：查询窗口长度、查询窗口宽度、查询窗口
角度和重叠率。
[0013] 优选地，所述的粒子图像测速矢量估计方法，其步骤（2）采用自适应
迭代算法更新查询窗口参数，具体过程为：
[0014] 将第n时刻的查询窗口参数作为初始查询窗口参数，使用互相关方法计算全场速度矢量、速度梯度大小方向和峰值强度，并根据全场速度矢量、速度梯度大小方向和峰值强度计算查询窗口参数估计，如果查询窗口参数估计与初始查询窗口参数的差值小于给定
阈值，则将更新查询窗口参数作为第n+1时刻的查询窗口参数，否则将查询窗口参数估计作为初始查询窗口参数，再次调整，直至查询窗口参数估计与初始查询窗口参数的差值小于给定阈值，将查询窗口参数估计作为第n+1时刻的查询窗口参数。
[0015] 优选地，所述的粒子图像测速矢量估计方法，其计算查询窗口参数估计的方法包括以下步骤：
[0016] （a）根据初始查询窗口参数，采用互相关计算，得到全场速度矢量U、速度梯度Ux、Uy和峰值强度P；
[0017] （b）计算全场速度矢量U的海瑟张量矩阵的
特征向量，获得当前时刻流场动态特性信息：
[0018] 海瑟张量矩阵：
[0019]
[0020] 其中，Uxx、Uxy、Uxy和Uyy为全场速度矢量U的四个二阶偏导数；
[0021] 根据海瑟张量结构矩阵，计算其特征值：λ1、λ2、特征向量 最小
曲率半径rmin，最大
曲率半径rmax；角度θ；λ1、λ2为det|H-λI|=0的两个根， 为(H-λI)X=0对应的两个
基础解系，最小曲率半径rmin、最大曲率半径rmax、角度θ其计算公式分别为：
[0022]
[0023]
[0024]
[0025] （c）估计查询窗口参数：
[0026] 查询窗口长度lx估计为：
[0027]
[0028] 查询窗口宽度估计为：
[0029] ly=lx/ε
[0030] 查询窗口角度估计为：
[0031]
[0032] 其中，长宽比系数ε=lx/ly=k*rmax/rmin，补偿系数 l为上次使用的窗口大小；
[0033] 查询窗口重叠率（overlap）估计，按照Uxx和Uyy越大，窗口重叠率增大的原则调整。
[0034] 优选地，所述的粒子图像测速矢量估计方法，对所述rmax/rmin进行阈值约束，使rmax/rmin不超过一定阈值。
[0035] 优选地，所述的粒子图像测速矢量估计方法，所述rmax/rmin的阈值等于4。
[0036] 优选地，所述的粒子图像测速矢量估计方法，其窗口重叠率，在0.4、0.5、0.6、0.8中选择。
[0037] 优选地，所述的粒子图像测速矢量估计方法，其对第n时刻和第n+1时刻图进行降噪的处理的方法可采用高斯滤波方法、均值滤波方法或中值滤波方法。
[0038] 优选地，所述的粒子图像测速矢量估计方法，其对第n时刻和第n+1时刻图像进行PIV流场矢量估计可采用基于FFT-CC的互相关算法或者梯度光流法。
[0039] 优选地，所述的粒子图像测速矢量估计方法，其对所述全场速度矢量，在局部邻域内使用
滤波器，如高斯权重滤波器，剔除掉一些误矢量，然后利用线性插值函数将剔除的矢量重新插补，得到最终的全场速度矢量U。
[0040] 总体而言，通过本发明所构思的以上技术方案与现有技术相比，由于采用引入了空间分辨率反馈参数，突破了传统的粒子图像测速方法单向计算的计算模式，能自动适应流畅局部流动的特点，自动调整查询窗口参数，无需人工干预地优化局部查询窗口尺寸和形状；同时计算的时间复杂度降低，能实现在线查询窗口参数调整；本发明提供的空间分辨率自适应调整的粒子图像测速矢量估计方法鲁棒性较好，提高了分析处理的空间分辨率，对低速和告诉流畅测量精度都取得了提高。
附图说明
[0041] 图1是空间分辨率自适应调整的粒子图像测速适量估计方法的
流程图；
[0043] 图3是计算流场速度矢量的互相关函数分布图；
具体实施方式
[0044] 为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0045] 本发明提供的一种空间分辨率自适应调整的粒子图像测速矢量估计方法，包括以下步骤：
[0046] （1）图像降噪预处理
[0047] 对第n时刻图像和第n+1时刻图像进行降噪处理，降噪方法可采用高斯滤波方法、均值滤波方法、中值滤波方法等滤波方法。
[0048] （2）根据第n时刻的查询窗口参数更新第n+1时刻的查询窗口参数。
[0049] 本步骤的技术思路为：根据前一组图像中各区域内的粒子速度矢量，得到各区域内粒子
加速度在X和Y向的变化率，若某区域粒子加速度在X向的变化率相对其Y向较大，则当前时刻在该区域内使用长宽比较大的查询窗口，反之，则当前时刻在该区域内使用长宽比较小的查询窗口，通过上述的查询窗口调整，能够保证查询窗口内的粒子密度恒定，从而保证最终矢量估计的准确性。
[0050] 具体来说是通过自适应迭代算法，更新查询窗口参数：首先将第n时刻的查询窗口参数作为初始查询窗口参数，然后根据初始查询窗口参数，使用互相关方法计算全场速度矢量U、速度梯度Ux、Uy和峰值强度P，计算查询窗口参数估计，如果查询窗口参数估计与初始查询窗口参数的差值小于给定阈值，则将查询窗口参数估计作为第n+1时刻的查询窗口参数，否则将查询窗口参数估计作为初始查询窗口参数，进行再次调整，直至查询窗口参数估计与初始查询窗口参数的差值小于给定阈值，则将查询窗口参数估计作为第n+1时刻的查询窗口参数。
[0051] 其中第1时刻的查询窗口参数，人工指定。
[0052] 计算查询窗口参数估计的具体方法如下：
[0053] （a）根据初始查询窗口参数，采用互相关计算，得到全场速度矢量U、速度梯度Ux、Uy和峰值强度P。
[0054] 同时计算峰值强度P和速度梯度Ux，Uy反馈信息。其计算公式分别为：
[0055] U=(u,v)={(mk,nk)|max(φfg(mk,nk))}
[0056] P=max{φfg(m,n)}，
[0057] Ux(i,j)=(U(i+1,j)-U(i-1,j))/2
[0058] Uy(i,j)=(U(i,j+1)-U(i,j-1))/2
[0059] 其中，矢量场U是一个二维向量，u是U的X向分量，v为Y向分量，k为图像被查询窗口分割成的区域数量，φfg(m,n)表示互相关系数，i,j代表图像内
像素点的坐标。互相关系数由下式计算：
[0060]
[0061] （b）计算全场速度矢量U的海瑟(Hessian)矩阵的特征向量，获得当前时刻流场动态特性信息，从而计算更新查询窗口参数。
[0062] 根据初始查询窗口参数，计算得到的全场速度矢量U，构建海瑟张量结构矩阵（Hessian tensor）：
[0063]
[0064] 其中，Uxx、Uxy、Uxy和Uyy为全场速度矢量U的四个二阶偏导数；根据海瑟张量结构矩阵，计算其特征值：λ1、λ2、特征向量 最小曲率半径rmin，最大曲率半径rmax；角度θ。λ1、λ2为det|H-λI|=0的两个根， 为(H-λI)X=0对应的两个基础解系。最小曲率半径rmin、最大曲率半径rmax的计算公式分别为：
[0065]
[0066]
[0067]
[0068] 估计查询窗口的长宽比系数ε：
[0069] ε=lx/ly=k*rmax/rmin
[0070] 其中，k为长宽比补偿系数，优选范围 对rmax/rmin进行阈值约束以防止窗口长宽比异常产生误矢量。
[0071] 长度lx估计为：
[0072]
[0073] 其中，P为峰值强度，l为上初始窗口长度。
[0074] 则窗口的宽度估计为：
[0075] ly=lx/ε；
[0076] 查询窗口角度估计：
[0077]
[0078] 对rmax/rmin进行阈值约束以防止窗口长宽比异常产生误矢量。
[0079] 查询窗口重叠率(overlap)估计，依据Uxx和Uyy确定，一般Uxx和Uyy越大，则窗口重叠率越大，具体应用时，窗口重叠率可在0.4、0.5、0.6、0.8中按照上述原则选取。
[0080] （3）采用第n+1时刻的查询窗口参数，对第n时刻的图像和第n+1进行PIV流场矢量估计。
[0081] 根据最终查询窗口参数，再次采用传统互相关方法或者基于FFT-CC的互相关算法或者使用梯度光流法，得到粒子速度矢量分量u和v。
[0082] 对得到的全场速度矢量，在局部邻域内使用滤波器，如高斯权重滤波器，剔除掉一些误矢量，然后利用线性插值函数将剔除的矢量重新插补，得到最终的全场速度矢量U，并将最终的全场速度矢量U输出显示。
[0083] 以下为实施例：
[0085] 如图2所示，图像序列描述的是一个微型管道内，流体自左往右流动的图像，第1时刻图像A、第2时刻图像B、间隔时间为Δt。预处理主要是对图像去噪，采用高斯滤波对第一时刻图像A和第二时刻图像B进行降噪处理。
[0086] （2）根据第1时刻的查询窗口参数，采用迭代算法，更新第2时刻的查询窗口参数。
[0087] 迭代算法的具体做法为:
[0088] 将第1时刻的查询窗口参数作为初始查询窗口参数，使用互相关方法计算全场速度矢量、速度梯度大小方向和峰值强度，并根据以上参数计算查询窗口参数估计，如果查询窗口参数估计与初始查询窗口参数的差值小于给定阈值，则将查询窗口参数估计作为第2时刻的查询窗口参数，否则将查询窗口参数估计作为初始查询窗口参数，再次调整，直至查询窗口参数估计与初始查询窗口参数的差值小于给定阈值，将查询窗口参数估计作为第n+1时刻的查询窗口参数。
[0089] 计算查询窗口估计的具体步骤为：
[0090] （a）根据初始查询窗口参数，采用互相关计算，得到全场速度矢量U、速度梯度Ux、Uy和峰值强度P。
[0091] 初始化查询窗口参数为：矩形窗口，长宽lx*ly为64*64、查询窗口角度θ为0、重叠率（overlap）为0，采用传统互相关方法进行速度矢量估计，互相关函数采用下式：
[0092]
[0093] 计算通过遍历(m,n)可以得到一个互相关平面，如图3所示，波峰P所在的位置(m,n)即认为是该区域内流场的位移，同时，互相关平面的峰值强度P可以反映该子图像区域的粒子浓度信息(波峰峰值越大，说明
信号强度越大)，这些信息可在下一步迭代计算中用来调整查询窗口的尺寸参数（长宽）和角度参数。
[0094] 整幅图片计算完成后，可得到全场速度矢量U、速度梯度Ux、Uy和峰值强度P。
[0095] U=(u,v)={(mk,nk)|max(φfg(mk,nk))}
[0096] P=max{φfg(m,n)}
[0097] Ux(i,j)=(U(i+1,j)-U(i-1,j))/2
[0098] Uy(i,j)=(U(i,j+1)-U(i,j-1))/2
[0099] 其中，全场速度矢量U是一个二维向量，u是U的x向分量，v为y向分量，k为图像被查询窗口分割成的区域数量，φfg(m,n)表示互相关系数，i,j代表图像内像素点的坐标。
[0100] （b）计算全场速度矢量U的海瑟(Hessian)矩阵的特征向量，获得当前时刻流场动态特性信息，从而计算更新查询窗口参数。
[0101] 利用上述获得的速度矢量场U反馈信息和峰值强度反馈信息对窗口进行重新调整，再进行对A和B的速度矢量场重新进行计算。
[0102] 基本思想是：在梯度大的区域，例如漩涡处，较小的问询域更为合适，而在粒子密度低并有较强的脱离测量平面运动分量的情况下则需要设置更大的问询域。此外用加长的问询域窗口形状，使长轴垂直于局部流动梯度最大的方向可以使相互垂直的两个方向中的一个方向上的空间分辨率增强，同时保证在长轴方向上粒子有足够大的位移像素数目，以保证速度矢量的测量精度。这种处理方法对沿着流体具有强剪切流动的区域更为重要。在这些区域，问询域窗口可以自动地被调整为沿着边界排列。具体如下：
[0103] 根据上一次计算得到的粒子速度矢量场U，构建海瑟张量结构矩阵(Hessian tensor)
[0104]
[0105] 其中，Uxx、Uxy、Uxy和Uyy为全场速度矢量U的四个二阶偏导数；根据海瑟张量结构矩阵，计算其特征值：λ1、λ2、特征向量 最小曲率半径rmin，最大曲率半径rmax；角度θ。λ1、λ2为det|H-λI|=0的两个根， 为(H-λI)X=0对应的两个基础解系，最小曲率半径rmin、最大曲率半径rmax及角度θ其计算公式分别为：
[0106]
[0107]
[0108]
[0109] 计算查询窗口的长宽比系数估计ε=lx/ly=k*rmax/rmin， 本例取0.8，作为优化，对rmax/rmin进行阈值约束以防止窗口长宽比异常产生误矢量，本例取4。
[0110] 长度lx估计为：
[0111]
[0112] 其中，P为峰值强度，l为上次使用的窗口大小。
[0113] 窗口的宽度估计为：
[0114] ly=lx/ε
[0115] 查询窗口角度估计为： 由于本例算得的自适应角度都比较小，所以为了简化计算，最终都赋值为0。
[0116] 查询窗口重叠率(overlap)估计，依据Uxx和Uyy确定，一般Uxx和Uyy越大，则窗口重叠率越大，本例中取为0.2。
[0117] 迭代的具体过程为：
[0118] 窗口参数自适应调整后，再次采用传统互相关方法，得到粒子速度矢量分量u和v，就可以得到本次迭代计算查询窗口的长宽和角度参数，将这些参数代入本次的迭代计算，得到新的全场速度矢量U。然后再次根据本步骤重复迭代，直到得到的两次迭代的窗口大小和方向变化都小于一定阈值时，本例取为0.1，即可停止迭代过程
[0119] （3）采用第2时刻的查询窗口参数，对第1时刻的图像和第2进行PIV流场矢量估计，得到全场速度矢量U。
[0120] 在局部邻域内使用均值滤波器，剔除掉一些误矢量，然后利用线性插值函数将剔除的矢量重新插补并将最终的全场速度矢量，将不同速度的矢量按不同
颜色输出显示。
[0121] 本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何
修改、等同替换和改进等，均应包含在本发明的保护范围之内。
