CN103604947A - 一种时间分辨率自适应调整的流场状态测量方法 - Google Patents

Abstract

本发明公开了一种时间分辨率自适应调整的流场状态测量方法，该方法包括如下步骤：首先，根据第k+1时刻图像与第k时刻图像之间的时间间隔Δt_k，采集第k+1时刻图像，并对第k时刻和第k+1时刻的图像进行降噪处理；然后，据k时刻图像和第k+1时刻的图像，采用粒子图像测速方法，计算第k+1时刻的流场状态X(k+1)，包括：位移场x(k+1)、速度场x′(k+1)、加速度场a(k+1)，并输出显示；最后，更新第k+2时刻与第k+1时刻图像采集时间之间的时间差Δt_k+1，根据Δt_k+1采集第k+2时刻的图像。本发明提供的流场状态测量方法，突破现有技术的单项计算模式，根据流场动态变化，做时间分辨率自适应调整，鲁棒性好，测量精度高。

Description

一种时间分辨率自适应调整的流场状态测量方法

技术领域

本发明涉及图像处理，特别涉及一种时间分辨率自适应调整的流场状态测量方法。

背景技术

流场状态测量，通常使用粒子图像测速(Particle Image Velocimetry简称PIV)方法，是一种用多次摄像记录流场中粒子的位置，分析摄得图像，从而测出流动速度的方法。测量原理是在流体中添加可以能够与之共同移动的示踪粒子，使用激光片光源将之照亮，使得看不见的流体可视化为离散的亮点，通过相机记录图像并按照一定的规则识别出在同一幅图像或两幅图像上对应的粒子或粒子集合对，得到粒子或粒子集合的运动速度，并以此作为粒子或粒子集合所在位置的流场的速度。如今，PIV作为一种全场非接触无扰动测量方法广泛应用于实验流体力学、生物医学、航空航天、工业制造等诸多领域。

传统PIV相关计算方法就是从独立存在的两幅图像通过一定的判别规则，如互相关计算等，首先建立离散的查询窗口(假设窗口内所有示踪粒子以相同的速度做刚性运动)，通过计算粒子集合的概率最佳配准位移表征该窗口内所有粒子的速度矢量，最终得到流场中各区域的流速矢量，对速度梯度变化剧烈的非定常流场，无法实现高精度测量。

目前，国内专利大部分集中在针对某一特定流场的PIV测量装置和测量流程，图像处理方面，中国专利CN200910109430公开了一种粒子图像测速处理方法，其中公开了通过构建水平集函数和最小化能量函数的方法，实现PIV图像矢量估计的处理方法。然而，该方法仍然存在以下缺陷：首先，该方法虽然在平滑约束项和基本约束项中了做了优化改进，该方法仍是一个单向性的计算过程，对于不同流场，其精度受计算过程参数的影响较大，鲁棒性不够；其次，对于高速流场，由于速度变化剧烈，该方法无法保证算法做出自适应调整，计算精度会降低；最后，该方法计算复杂度也较大，这就决定了该方法在一些时变流场中无法实时跟踪其动态特性，只能离线分析其结果。

发明内容

本发明针对现有的流场状态测量方法，不能适应变化剧烈的流场，同时测量精度不高的问题，本发明提供了一种时间分辨率自适应调整的流场状态测量方法，其目的在与根据流场状态自适应的调整图像采集的时间间隔，从而适应剧烈变化的流场，得出更为精确的测量结果。

按照本发明，提供了一种时间分辨率自适应调整的流场状态测量方法，对于每个图像采集时刻，按照以下步骤对当前流场状态进行测量：

（1）根据第k+1时刻图像与第k时刻图像之间的时间间隔Δt_k，采集第k+1时刻图像，并对第k时刻和第k+1时刻的图像进行降噪处理；

（2）根据k时刻图像和第k+1时刻的图像，采用粒子图像测速方法，计算第k+1时刻的流场状态X(k+1)，包括：位移场x(k+1)、速度场x′(k+1)、加速度场a(k+1)，并输出显示；

（3）更新第k+2时刻与第k+1时刻图像采集时间之间的时间差Δt_k+1第k+2时刻与第k+1时刻图像采集时间之间的时间差Δt_k+1，按照如下公式确定：

${\mathrm{\&Delta;t}}_{k+1}=\left\{\begin{array}{cc}\frac{\sqrt{x^{\&prime;}{\left(k\right)}^2+2a\left(k\right)}-x^{\&prime;}\left(k\right)}{a\left(k\right)}& (x(k+2)\&GreaterEqual;a)\\ {\mathrm{\&Delta;t}}_k& (b\&le;x(k+2)<a)\\ \frac{\sqrt{x^{\&prime;}{\left(k\right)}^2+1.2a\left(k\right)}-x^{\&prime;}\left(k\right)}{a\left(k\right)}& (x(k+2)<b)\end{array}\right.$

其中，x′(k)为第k时刻的速度场，a(k)为第k时刻的加速度场，x(k+2)为采用卡尔曼一步预测算法根据第k+1时刻对k+2时刻的流场状态预测值

的位移矢量场，a、b为阈值；

其中第0时刻图像与第1时刻图像之间的时间间隔Δt₁，第0时刻的流场状态X(0)由人工指定。

优选地，所述的流场状态测量方法，其步骤（1）的降噪处理，可采用高斯滤波方法、均值滤波方法或中值滤波方法。

优选地，所述的流场状态测量方法，其步骤（2）采用的粒子图像测速方法为互相关算法或光流法。

优选地，所述的流场状态测量方法，其步骤（2）采用的粒子图像测速方法，采用卡尔曼一步预测法，根据第k时刻对k+1时刻的流场状态预测值

和第k+1时刻的流场状态测定值Z(k+1)，估计第k+1时刻流场状态最佳估计作为第k+1时刻流场状态X(k+1)。

优选地，所述的流场状态测量方法，其步骤（3）中阈值a处于1至2.5之间，b处于0.4至0.7之间。

总体而言，按照本发明的自适应测量方法相对于现有技术，主要可获得以下方面的技术效果：

1.突破传统PIV单向的计算模式，引入时间分辨率自适应调整，最终实现对流场动态变化的时间分辨率自适应调整，鲁棒性较好，同时对低速和高速流场，计算精度都很高(<0.01pixel)；

2.PIV矢量估计使用基于卡尔曼一步预测的理论，对粒子的下一时刻动态信息进行预测，使PIV矢量估计计算复杂度降低低、能够有效自适应流场的动态特性复杂变化，高频的测量结果输出能够实现对时变流场的高精度测量。

附图说明

图1是本发明的粒子图像测速时间分辨率自适应测量方法原理图；

图2是实施例所用的例子图像示意图；

图3是本发明的粒子图像测速时间分辨率自适应测量方法的参数计算流程图；

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明提供的一种空间分辨率自适应调整的粒子图像测速矢量估计方法，如图1所示，包括以下步骤：

（1）图像序列预处理

对第k时刻图像和第k+1时刻图像进行降噪处理，降噪方法可采用高斯滤波方法、均值滤波方法、中值滤波方法等滤波方法。

（2）估计当前流场状态

采用粒子图像测速方法估计当前流场状态，如PIV互相关算法和光流法，参见MarkusRaffel、Christian E.Willert等出版的《Particle ImageVelocimetry A Practical Guide》一书的第5章。计算第k+1时刻流场状态X(k+1)，包括：位移场x(k+1)、速度场x′(k+1)、加速度场a(k+1)。

本发明另外提供一种更为精确的粒子图像测速方法，用于精确估计流场状态：

第0时刻图像与第1时刻图像之间的时间间隔Δt₁，第0时刻的流场状态X(0)由人工指定，采用卡尔曼一步预测法，根据第k时刻对k+1时刻的流场状态预测值和第k+1时刻的流场状态测定值Z(k+1)，对第k+1时刻流场状态作出最佳估计

并将第k+1时刻流场状态最佳估计

作为第k+1时刻流场状态X(k+1)，即第k+1时刻流场状态X(k+1)可表示为：

$X(k+1)=\hat{X}(k+1|k+1)$

其中，

为第k+1时刻流场状态最佳估计，根据第k时刻对k+1时刻的流场状态预测值

和第k+1时刻的流场状态测定值Z(k+1)计算，其公式为：

$\hat{X}(k+1|k+1)=\hat{X}(k+1|k)+K(k+1)[Z(k+1)-H(k+1)\hat{X}(k+1|k)]$

其中：为第k时刻对k+1时刻的流场状态预测值；K(k+1)为第k时刻到第k+1时刻的卡尔曼增益；Z(k+1)为第k+1时刻的流场状态测定值，根据第k时刻图像和第k+1时刻图像，利用PIV互相关或者光流法计算；H(k+1)为测量系统的模型参数。

其计算公式分别如下：

第k时刻对k+1时刻的流场状态预测值

$\hat{X}(k+1|k)=\mathrm{\&phi;}(k+1,k)\hat{X}\left(k\right|k)+U\left(k\right)$

第k时刻到第k+1时刻的卡尔曼增益K(k+1)：

K(k+1)=P(k+1|k)H^T(k+1)[H(k+1)P(k+1|k)H^T(k+1)+R(k+1)]^-1

测量流场的模型参数H(k+1)：

H(k+1)=(Z(k+1)-V(k+1))(X(k+1)^-1)

其中：φ(k+1,k)为系统参数， $\mathrm{\&phi;}(k+1,k)=\left[\begin{array}{cc}1& {\mathrm{\&Delta;t}}_k\\ 0& 1\end{array}\right],$ Δt_k为第k时刻图像与第k+1时刻图像采集间隔；

为第k时刻流场状态最佳估计；U(k)为第k时刻对系统的控制量， $U\left(k\right)={\left[\begin{array}{cc}\frac{1}{2}a\left(k\right)\mathrm{\&Delta;}{t_k}^2& a\left(k\right){\mathrm{\&Delta;t}}_k\end{array}\right]}^T,$ a(k)为第k时刻流场状态中的加速度；P(k+1|k)为

对应的协方差矩阵；H^T(k+1)为H(k+1)的转置矩阵；R(k+1)为第k时刻到第k+1时刻测量过程的噪声V(k+1)的协方差，假设V(k+1)为高斯白噪声。

对应的协方差矩阵的P(k+1|k)获取方法为：

P(k+1|k)=φ(k+1,k)P(k|k)φ^T(k+1,k)

P(k|k)=[I-K(k)H(k)]P(k|k-1)

其中，φ(k+1,k)为系统参数；φ^T(k+1,k)为φ(k+1,k)的转置矩阵；P(k|k)为

对应的协方差矩阵，I为单位矩阵；K(k)第k-1时刻到第k时刻的卡尔曼增益；H(k)为第k时刻测量系统的模型参数；

对应的协方差矩阵。

（3）更新第k+2时刻与第k+1时刻图像采集时间之间的时间差Δt_k+1。

第0时刻图像与第1时刻图像之间的时间间隔Δt₁，第0时刻的流场状态X(0)由人工指定，根据卡尔曼一步预测算法预测的第k+1时刻对k+2时刻的流场状态预测值x(k+2)和第k时刻的速度场x′(k)，更新第k+2时刻与第k+1时刻图像采集时间之间的时间差Δt_k+1，具体方法如下：

${\mathrm{\&Delta;t}}_{k+1}=\left\{\begin{array}{cc}\frac{\sqrt{x^{\&prime;}{\left(k\right)}^2+2a\left(k\right)}-x^{\&prime;}\left(k\right)}{a\left(k\right)}& (x(k+2)\&GreaterEqual;a)\\ {\mathrm{\&Delta;t}}_k& (b\&le;x(k+2)<a)\\ \frac{\sqrt{x^{\&prime;}{\left(k\right)}^2+1.2a\left(k\right)}-x^{\&prime;}\left(k\right)}{a\left(k\right)}& (x(k+2)<b)\end{array}\right.$

其中，x′(k)为第k时刻的速度场，x(k+2)为第k+1时刻对k+2时刻的流场状态预测值

的位移矢量场，a(k)为第k时刻的加速度场，a、b为阈值，a处于1至2.5之间，b处于0.4至0.7之间。

以下为实施例：

实施例1

如附图2所示，图像序列描述的是一个微型管道内，流体自左往右流动的图像。假定k时刻到k+1时刻图像采集间隔为Δt_k，图像大小为512*512pixel，由于高速采集装置采图时间较短，在图像序列中相邻的若干帧里，目标可以看作匀速直线运动或定加速直线运动，目标粒子运动模型可用下列差分方程来描述(以直角坐标系统中x轴方向为例，y轴方向与之类似)。粒子模型可描述为：

$x(k+1)=x\left(k\right)+{\mathrm{\&Delta;t}}_k{x}^{\&prime;}\left(k\right)+\frac{1}{2}{a}_x\left(k\right){{\mathrm{\&Delta;t}}_k}^2$

x′(k+1)=x′(k)+Δt_ka_x(k)

其中，x(k)、x′(k)分别为k时刻PIV图像中示踪粒子x坐标的目标位置和速度，a_x(k)为加速度，认为状态变量，加速度的随机扰动大多数情况下并不是白噪声,是平稳随机序列，服从正态分布，均值为0，由于对象的时变动态特性，所以假定某一时刻的加速度与另一时刻的加速度不相关图像序列预处理。

(1)图像序列预处理

对第k时刻图像和第k+1时刻图像进行降噪处理，降噪方法可采用高斯滤波方法。

(2)估计当前流场速度矢量

整个参数计算流程图如附图3所示。

根据第k时刻对k+1时刻的流场状态预测值

和第k+1时刻的流场状态测定值Z(k+1)，对沿着匀速直线附加随机加速度的PIV粒子目标第k+1时刻的流场状态计算最佳估计

并将第k+1时刻流场状态最佳估计

作为第k+1时刻流场状态即第k+1时刻流场状态X(k+1)可表示为：

$X(k+1)=\hat{X}(k+1|k+1)$

其中，

为第k+1时刻流场状态最佳估计，根据第k时刻对k+1时刻的流场状态预测值

和第k+1时刻的流场状态测定值Z(k+1)计算，其公式为：

$\hat{X}(k+1|k+1)=\hat{X}(k+1|k)+K(k+1)[Z(k+1)-H(k+1)\hat{X}(k+1|k)]$

其中：为第k时刻对k+1时刻的流场状态预测值；K(k+1)为第k时刻到第k+1时刻的卡尔曼增益；Z(k+1)为第k+1时刻的流场状态测定值，根据第k时刻图像和第k+1时刻图像，利用PIV互相关或者光流法计算；H(k+1)为测量系统的模型参数。

其计算公式分别如下：

第k时刻对k+1时刻的流场状态预测值

$\hat{X}(k+1,k)=\mathrm{\&phi;}(k+1,k)\hat{X}\left(k\right|k)+U\left(k\right)$

第k时刻到第k+1时刻的卡尔曼增益K(k+1)：

K(k+1)=P(k+1|k)H^T(k+1)[H(k+1)P(k+1|k)H^T(k+1)+R(k+1)]^-1

测量流场的模型参数H(k+1)：

H(k+1)=(Z(k+1)-V(k+1))(X(k+1)^-1)

其中：φ(k+1,k)为系统参数，本实例中， $\mathrm{\&phi;}(k+1,k)=\left[\begin{array}{cc}1& {\mathrm{\&Delta;t}}_k\\ 0& 1\end{array}\right],$ Δt_k为第k时刻图像与第k+1时刻图像采集间隔；为第k时刻流场状态最佳估计；U(k)为第k时刻对系统的控制量，本实例中， $U\left(k\right)={\left[\begin{array}{cc}\frac{1}{2}a\left(k\right)\mathrm{\&Delta;}{t_k}^2& a\left(k\right){\mathrm{\&Delta;t}}_k\end{array}\right]}^T,$ a(k)为第k时刻流场状态中的加速度；P(k+1|k)为

对应的协方差矩阵；H^T(k+1)为H(k+1)的转置矩阵；本实例中， $H(k+1)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right],$ R(k+1)为第k时刻到第k+1时刻测量过程的噪声V(k+1)的协方差，本实例中，R(k+1)=0；假设V(k+1)为高斯白噪声。

对应的协方差矩阵的P(k+1|k)获取方法为：

P(k+1|k)=φ(k+1,k)P(k|k)φ^T(k+1,k)

P(k|k)=[I-K(k)H(k)]P(k|k-1)

其中，φ(k+1,k)为系统参数；φ^T(k+1,k)为φ(k+1,k)的转置矩阵；P(k|k)为

对应的协方差矩阵，I为单位矩阵；K(k)第k-1时刻到第k时刻的卡尔曼增益；H(k)为第k时刻测量系统的模型参数；

对应的协方差矩阵。

其中计算的时候涉及到初始化值，取一个像素点(i,j)为例，本实例中初值赋值如下：

$X_{i,j}\left(0\right|0)=\left[\begin{array}{c}{x}_{i,j}\left(0\right)\\ x_{i,j}^{\&prime;}\left(0\right)\end{array}\right]=\left[\begin{array}{c}1\\ \frac{1}{2}\end{array}\right]$

$P_{i,j}\left(0\right)=\mathrm{var}\left[\begin{array}{c}{x}_{i,j}\left(0\right)\\ x_{i,j}^{\&prime;}\left(0\right)\end{array}\right]=\left[\begin{array}{cc}0.5& 0\\ 0& 0.045\end{array}\right]$

系统测量值Z(k+1)使用PIV快速矢量测量算法，本例使用互相关算法，互相关系数φ_fg(m,n)由下式计算：

${\mathrm{\&phi;}}_{\mathrm{fg}}(m,n)=\frac{\stackrel{M}{\mathrm{\&Sigma;}}\stackrel{N}{\mathrm{\&Sigma;}}[f(k,l)-\stackrel{\&OverBar;}{f}][g(k+m,l+n)-\stackrel{\&OverBar;}{g}]}{\sqrt{{[\stackrel{M}{\mathrm{\&Sigma;}}\stackrel{N}{\mathrm{\&Sigma;}}f(k,l)-\stackrel{\&OverBar;}{f}]}^2}{[\stackrel{M}{\mathrm{\&Sigma;}}\stackrel{N}{\mathrm{\&Sigma;}}g(k,l)-\stackrel{\&OverBar;}{g}]}^2}$

(3)更新第k+2时刻与第k+1时刻图像采集时间之间的时间差Δt_k+1。

第0时刻图像与第1时刻图像之间的时间间隔Δt₁，第0时刻的流场状态X(0)由人工指定，根据卡尔曼一步预测算法预测的第k+1时刻对k+2时刻的流场状态预测值x(k+2)和第k时刻的速度场x′(k)，更新第k+2时刻与第k+1时刻图像采集时间之间的时间差Δt_k+1，具体方法如下：

${\mathrm{\&Delta;t}}_{k+1}=\left\{\begin{array}{cc}\frac{\sqrt{x^{\&prime;}{\left(k\right)}^2+2a\left(k\right)}-x^{\&prime;}\left(k\right)}{a\left(k\right)}& (x(k+2)\&GreaterEqual;a)\\ {\mathrm{\&Delta;t}}_k& (b\&le;x(k+2)<a)\\ \frac{\sqrt{x^{\&prime;}{\left(k\right)}^2+1.2a\left(k\right)}-x^{\&prime;}\left(k\right)}{a\left(k\right)}& (x(k+2)<b)\end{array}\right.$

其中，x′(k)为第k时刻的速度场，a(k)为第k时刻的加速度场，x(k+2)为第k+1时刻对k+2时刻的流场状态预测值

的位移矢量场，a、b为阈值，a为1，b为0.4。

实施例2

如附图2所示，图像序列描述的是一个微型管道内，流体自左往右流动的图像。假定k时刻到k+1时刻图像采集间隔为Δt_k，图像大小为512*512pixel，由于高速采集装置采图时间较短，在图像序列中相邻的若干帧里，目标可以看作匀速直线运动或定加速直线运动，目标粒子运动模型可用下列差分方程来描述(以直角坐标系统中x轴方向为例，y轴方向与之类似)。粒子模型可描述为：

$x(k+1)=x\left(k\right)+{\mathrm{\&Delta;t}}_k{x}^{\&prime;}\left(k\right)+\frac{1}{2}{a}_x\left(k\right){{\mathrm{\&Delta;t}}_k}^2$

x′(k+1)=x′(k)+Δt_ka_x(k)

其中，x(k)、x′(k)分别为k时刻PIV图像中示踪粒子x坐标的目标位置和速度，a_x(k)为加速度，认为状态变量，加速度的随机扰动大多数情况下并不是白噪声,是平稳随机序列，服从正态分布，均值为0，由于对象的时变动态特性，所以假定某一时刻的加速度与另一时刻的加速度不相关图像序列预处理。

(1)图像序列预处理

对第k时刻图像和第k+1时刻图像进行降噪处理，降噪方法可采用高斯滤波方法。

(2)估计当前流场速度矢量

整个参数计算流程图如附图3所示。

根据第k时刻对k+1时刻的流场状态预测值

和第k+1时刻的流场状态测定值Z(k+1)，对沿着匀速直线附加随机加速度的PIV粒子目标第k+1时刻的流场状态计算最佳估计

并将第k+1时刻流场状态最佳估计

作为第k+1时刻流场状态

即第k+1时刻流场状态X(k+1)可表示为：

$X(k+1)=\hat{X}(k+1|k+1)$

其中，

为第k+1时刻流场状态最佳估计，根据第k时刻对k+1时刻的流场状态预测值

和第k+1时刻的流场状态测定值Z(k+1)计算，其公式为：

$\hat{X}(k+1|k+1)=\hat{X}(k+1|k)+K(k+1)[Z(k+1)-H(k+1)\hat{X}(k+1|k)]$

其中：

为第k时刻对k+1时刻的流场状态预测值；K(k+1)为第k时刻到第k+1时刻的卡尔曼增益；Z(k+1)为第k+1时刻的流场状态测定值，根据第k时刻图像和第k+1时刻图像，利用PIV互相关或者光流法计算；H(k+1)为测量系统的模型参数。

其计算公式分别如下：

第k时刻对k+1时刻的流场状态预测值

$\hat{X}(k+1|k)=\mathrm{\&phi;}(k+1,k)\hat{X}\left(k\right|k)+U\left(k\right)$

第k时刻到第k+1时刻的卡尔曼增益K(k+1)：

K(k+1)=P(k+1|k)H^T(k+1)[H(k+1)P(k+1|k)H^T(k+1)+R(k+1)]^-1

测量流场的模型参数H(k+1)：

H(k+1)=(Z(k+1)-V(k+1))(X(k+1)^-1)

其中：φ(k+1,k)为系统参数，本实例中， $\mathrm{\&phi;}(k+1,k)=\left[\begin{array}{cc}1& {\mathrm{\&Delta;t}}_k\\ 0& 1\end{array}\right],$ Δt_k为第k时刻图像与第k+1时刻图像采集间隔；

为第k时刻流场状态最佳估计；U(k)为第k时刻对系统的控制量，本实例中， $U\left(k\right)={\left[\begin{array}{cc}\frac{1}{2}a\left(k\right)\mathrm{\&Delta;}{t_k}^2& a\left(k\right){\mathrm{\&Delta;t}}_k\end{array}\right]}^T,$ a(k)为第k时刻流场状态中的加速度；P(k+1|k)为

对应的协方差矩阵；H^T(k+1)为H(k+1)的转置矩阵；本实例中， $H(k+1)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right],$ R(k+1)为第k时刻到第k+1时刻测量过程的噪声V(k+1)的协方差，本实例中，R(k+1)=0；假设V(k+1)为高斯白噪声。

对应的协方差矩阵的P(k+1|k)获取方法为：

P(k+1|k)=φ(k+1,k)P(k|k)φ^T(k+1,k)

P(k|k)=[I-K(k)H(k)]P(k|k-1)

其中，φ(k+1,k)为系统参数；φ^T(k+1,k)为φ(k+1,k)的转置矩阵；P(k|k)为

对应的协方差矩阵，I为单位矩阵；K(k)第k-1时刻到第k时刻的卡尔曼增益；H(k)为第k时刻测量系统的模型参数；

对应的协方差矩阵。

其中计算的时候涉及到初始化值，取一个像素点(i,j)为例，本实例中初值赋值如下：

$X_{i,j}\left(0\right|0)=\left[\begin{array}{c}{x}_{i,j}\left(0\right)\\ x_{i,j}^{\&prime;}\left(0\right)\end{array}\right]=\left[\begin{array}{c}1\\ \frac{1}{2}\end{array}\right]$

$P_{i,j}\left(0\right)=\mathrm{var}\left[\begin{array}{c}{x}_{i,j}\left(0\right)\\ x_{i,j}^{\&prime;}\left(0\right)\end{array}\right]=\left[\begin{array}{cc}0.5& 0\\ 0& 0.045\end{array}\right]$

系统测量值Z(k+1)使用PIV快速矢量测量算法，本例使用互相关算法，互相关系数φ_fg(m,n)由下式计算：

${\mathrm{\&phi;}}_{\mathrm{fg}}(m,n)=\frac{\stackrel{M}{\mathrm{\&Sigma;}}\stackrel{N}{\mathrm{\&Sigma;}}[f(k,l)-\stackrel{\&OverBar;}{f}][g(k+m,l+n)-\stackrel{\&OverBar;}{g}]}{\sqrt{{[\stackrel{M}{\mathrm{\&Sigma;}}\stackrel{N}{\mathrm{\&Sigma;}}f(k,l)-\stackrel{\&OverBar;}{f}]}^2}{[\stackrel{M}{\mathrm{\&Sigma;}}\stackrel{N}{\mathrm{\&Sigma;}}g(k,l)-\stackrel{\&OverBar;}{g}]}^2}$

(3)更新第k+2时刻与第k+1时刻图像采集时间之间的时间差Δt_k+1。

第0时刻图像与第1时刻图像之间的时间间隔Δt₁，第0时刻的流场状态X(0)由人工指定，根据卡尔曼一步预测算法预测的第k+1时刻对k+2时刻的流场状态预测值x(k+2)和第k时刻的速度场x′(k)，更新第k+2时刻与第k+1时刻图像采集时间之间的时间差Δt_k+1，具体方法如下：

${\mathrm{\&Delta;t}}_{k+1}=\left\{\begin{array}{cc}\frac{\sqrt{x^{\&prime;}{\left(k\right)}^2+2a\left(k\right)}-x^{\&prime;}\left(k\right)}{a\left(k\right)}& (x(k+2)\&GreaterEqual;a)\\ {\mathrm{\&Delta;t}}_k& (b\&le;x(k+2)<a)\\ \frac{\sqrt{x^{\&prime;}{\left(k\right)}^2+1.2a\left(k\right)}-x^{\&prime;}\left(k\right)}{a\left(k\right)}& (x(k+2)<b)\end{array}\right.$

其中，x′(k)为第k时刻的速度场，a(k)为第k时刻的加速度场，x(k+2)为第k+1时刻对k+2时刻的流场状态预测值

的位移矢量场，a、b为阈值，a为1.5，b为0.6。

实施例3

如附图2所示，图像序列描述的是一个微型管道内，流体自左往右流动的图像。假定k时刻到k+1时刻图像采集间隔为Δt_k，图像大小为512*512pixel，由于高速采集装置采图时间较短，在图像序列中相邻的若干帧里，目标可以看作匀速直线运动或定加速直线运动，目标粒子运动模型可用下列差分方程来描述(以直角坐标系统中x轴方向为例，y轴方向与之类似)。粒子模型可描述为：

$x(k+1)=x\left(k\right)+{\mathrm{\&Delta;t}}_k{x}^{\&prime;}\left(k\right)+\frac{1}{2}{a}_x\left(k\right){{\mathrm{\&Delta;t}}_k}^2$

x′(k+1)=x′(k)+Δt_ka_x(k)

其中，x(k)、x′(k)分别为k时刻PIV图像中示踪粒子x坐标的目标位置和速度，a_x(k)为加速度，认为状态变量，加速度的随机扰动大多数情况下并不是白噪声,是平稳随机序列，服从正态分布，均值为0，由于对象的时变动态特性，所以假定某一时刻的加速度与另一时刻的加速度不相关图像序列预处理。

(1)图像序列预处理

对第k时刻图像和第k+1时刻图像进行降噪处理，降噪方法可采用高斯滤波方法。

(2)估计当前流场速度矢量

整个参数计算流程图如附图3所示。

根据第k时刻对k+1时刻的流场状态预测值

和第k+1时刻的流场状态测定值Z(k+1)，对沿着匀速直线附加随机加速度的PIV粒子目标第k+1时刻的流场状态计算最佳估计

并将第k+1时刻流场状态最佳估计

作为第k+1时刻流场状态

即第k+1时刻流场状态X(k+1)可表示为：

$X(k+1)=\hat{X}(k+1|k+1)$

其中，

为第k+1时刻流场状态最佳估计，根据第k时刻对k+1时刻的流场状态预测值和第k+1时刻的流场状态测定值Z(k+1)计算，其公式为：

$\hat{X}(k+1|k+1)=\hat{X}(k+1|k)+K(k+1)[Z(k+1)-H(k+1)\hat{X}(k+1|k)]$

其中：

为第k时刻对k+1时刻的流场状态预测值；K(k+1)为第k时刻到第k+1时刻的卡尔曼增益；Z(k+1)为第k+1时刻的流场状态测定值，根据第k时刻图像和第k+1时刻图像，利用PIV互相关或者光流法计算；H(k+1)为测量系统的模型参数。

其计算公式分别如下：

第k时刻对k+1时刻的流场状态预测值

$\hat{X}(k+1|k)=\mathrm{\&phi;}(k+1,k)\hat{X}\left(k\right|k)+U\left(k\right)$

第k时刻到第k+1时刻的卡尔曼增益K(k+1)：

K(k+1)=P(k+1|k)H^T(k+1)[H(k+1)P(k+1|k)H^T(k+1)+R(k+1)]^-1

测量流场的模型参数H(k+1)：

H(k+1)=(Z(k+1)-V(k+1))(X(k+1)^-1)

其中：φ(k+1,k)为系统参数，本实例中， $\mathrm{\&phi;}(k+1,k)=\left[\begin{array}{cc}1& {\mathrm{\&Delta;t}}_k\\ 0& 1\end{array}\right],$ Δt_k为第k时刻图像与第k+1时刻图像采集间隔；

为第k时刻流场状态最佳估计；U(k)为第k时刻对系统的控制量，本实例中， $U\left(k\right)={\left[\begin{array}{cc}\frac{1}{2}a\left(k\right)\mathrm{\&Delta;}{t_k}^2& a\left(k\right){\mathrm{\&Delta;t}}_k\end{array}\right]}^T,$ a(k)为第k时刻流场状态中的加速度；P(k+1|k)为

对应的协方差矩阵；H^T(k+1)为H(k+1)的转置矩阵；本实例中， $H(k+1)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right],$ R(k+1)为第k时刻到第k+1时刻测量过程的噪声V(k+1)的协方差，本实例中，R(k+1)=0；假设V(k+1)为高斯白噪声。

对应的协方差矩阵的P(k+1|k)获取方法为：

P(k+1|k)=φ(k+1,k)P(k|k)φ^T(k+1,k)

P(k|k)=[I-K(k)H(k)]P(k|k-1)

其中，φ(k+1,k)为系统参数；φ^T(k+1,k)为φ(k+1,k)的转置矩阵；P(k|k)为

对应的协方差矩阵，I为单位矩阵；K(k)第k-1时刻到第k时刻的卡尔曼增益；H(k)为第k时刻测量系统的模型参数；

对应的协方差矩阵。

其中计算的时候涉及到初始化值，取一个像素点(i,j)为例，本实例中初值赋值如下：

$X_{i,j}\left(0\right|0)=\left[\begin{array}{c}{x}_{i,j}\left(0\right)\\ x_{i,j}^{\&prime;}\left(0\right)\end{array}\right]=\left[\begin{array}{c}1\\ \frac{1}{2}\end{array}\right]$

$P_{i,j}\left(0\right)=\mathrm{var}\left[\begin{array}{c}{x}_{i,j}\left(0\right)\\ x_{i,j}^{\&prime;}\left(0\right)\end{array}\right]=\left[\begin{array}{cc}0.5& 0\\ 0& 0.045\end{array}\right]$

系统测量值Z(k+1)使用PIV快速矢量测量算法，本例使用互相关算法，互相关系数φ_fg(m,n)由下式计算：

${\mathrm{\&phi;}}_{\mathrm{fg}}(m,n)=\frac{\stackrel{M}{\mathrm{\&Sigma;}}\stackrel{N}{\mathrm{\&Sigma;}}[f(k,l)-\stackrel{\&OverBar;}{f}][g(k+m,l+n)-\stackrel{\&OverBar;}{g}]}{\sqrt{{[\stackrel{M}{\mathrm{\&Sigma;}}\stackrel{N}{\mathrm{\&Sigma;}}f(k,l)-\stackrel{\&OverBar;}{f}]}^2}{[\stackrel{M}{\mathrm{\&Sigma;}}\stackrel{N}{\mathrm{\&Sigma;}}g(k,l)-\stackrel{\&OverBar;}{g}]}^2}$

(3)更新第k+2时刻与第k+1时刻图像采集时间之间的时间差Δt_k+1。

第0时刻图像与第1时刻图像之间的时间间隔Δt₁，第0时刻的流场状态X(0)由人工指定，根据卡尔曼一步预测算法预测的第k+1时刻对k+2时刻的流场状态预测值x(k+2)和第k时刻的速度场x′(k)，更新第k+2时刻与第k+1时刻图像采集时间之间的时间差Δt_k+1，具体方法如下：

${\mathrm{\&Delta;t}}_{k+1}=\left\{\begin{array}{cc}\frac{\sqrt{x^{\&prime;}{\left(k\right)}^2+2a\left(k\right)}-x^{\&prime;}\left(k\right)}{a\left(k\right)}& (x(k+2)\&GreaterEqual;a)\\ {\mathrm{\&Delta;t}}_k& (b\&le;x(k+2)<a)\\ \frac{\sqrt{x^{\&prime;}{\left(k\right)}^2+1.2a\left(k\right)}-x^{\&prime;}\left(k\right)}{a\left(k\right)}& (x(k+2)<b)\end{array}\right.$

其中，x′(k)为第k时刻的速度场，a(k)为第k时刻的加速度场，x(k+2)为第k+1时刻对k+2时刻的流场状态预测值的位移矢量场，a、b为阈值，a为2.5，b为0.7。

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

Claims (5)

1.一种时间分辨率自适应调整的流场状态测量方法，其特征在于，对于每个图像采集时刻，按照以下步骤对当前流场状态进行测量：

（1）根据第k+1时刻图像与第k时刻图像之间的时间间隔Δt_k，采集第k+1时刻图像，并对第k时刻和第k+1时刻的图像进行降噪处理；

（2）根据k时刻图像和第k+1时刻的图像，采用粒子图像测速方法，计算第k+1时刻的流场状态X(k+1)，包括：位移场x(k+1)、速度场x′(k+1)、加速度场a(k+1)；

（3）更新第k+2时刻与第k+1时刻图像采集时间之间的时间间隔Δt_k+1：

${\mathrm{\&Delta;t}}_{k+1}=\left\{\begin{array}{cc}\frac{\sqrt{x^{\&prime;}{\left(k\right)}^2+2a\left(k\right)}-x^{\&prime;}\left(k\right)}{a\left(k\right)}& (x(k+2)\&GreaterEqual;a)\\ {\mathrm{\&Delta;t}}_k& (b\&le;x(k+2)<a)\\ \frac{\sqrt{x^{\&prime;}{\left(k\right)}^2+1.2a\left(k\right)}-x^{\&prime;}\left(k\right)}{a\left(k\right)}& (x(k+2)<b)\end{array}\right.$

其中，x′(k)为第k时刻的速度场，a(k)为第k时刻的加速度场，x(k+2)为采用卡尔曼一步预测算法根据第k+1时刻对k+2时刻的流场状态预测值

的位移矢量场，a、b为阈值。

2.如权利要求1所述的流场状态测量方法，其特征在于，所述步骤（1）的降噪处理，可采用高斯滤波方法、均值滤波方法或中值滤波方法。

3.如权利要求1或2所述的流场状态测量方法，其特征在于，所述步骤（2）采用的粒子图像测速方法为互相关算法或光流法。

4.如权利要求1或2所述的流场状态测量方法，其特征在于，所述步骤（2）采用的粒子图像测速方法的具体实现方式为：采用卡尔曼一步预测法，根据第k时刻对k+1时刻的流场状态预测值

和第k+1时刻的流场状态测定值Z(k+1)，估计第k+1时刻流场状态最佳估计

作为第k+1时刻流场状态X(k+1)。

5.如权利要求1所述的流场状态测量方法，其特征在于，所述步骤（3）中阈值a处于1至2.5之间，b处于0.4至0.7之间。

Images