CN101738489B - 一种测量散射性流体横向流速的方法 - Google Patents

Abstract

本发明公开了一种测量散射性流体横向流速的方法，包括以下步骤：选择测量用在前流体截面和在后流体截面；用OCT样品臂扫描光束在前后截面进行交替扫描；重建OCT结构图得到前后两个截面位置处各自的结构图像序列；提取两结构图像序列中对应同一点位置处的散射强度进行互相关运算；确定散射介质流过在前流体截面和在后流体截面的延迟时间；确定在流体截面某一位置处的散射介质的横向流速测量结果；重复上述步骤，确定散射流体在截面各个位置处的横向流速测量结果即测量出流体的横向流速分布。本发明的横向流速测量方法不依赖于对流体夹角进行测量，可以在未知OCT系统样品臂端口发出的扫描光束入射角的情况下测量出流体的横向流速分布。

Description

一种测量散射性流体横向流速的方法

技术领域

本发明涉及流速测量，特别是涉及一种测量散射性流体横向流速的方法。

背景技术

多普勒光学相干层析成像(Doppler optical coherence tomography，缩写为DOCT)技术结合了光学多普勒效应和OCT技术，能以非侵入无损伤方式提供微米级的空间结构分辨率和流体流速信息。一般是从光纤耦合器样品臂端口发出扫描光束，以多普勒角斜入射被测量流体，然后利用流体流速在入射光方向上的分量产生的多普勒效应，通过实测流体中含有散射介质——运动粒子改变了频率的散射光，采集多普勒频移图像数据，由多普勒信号处理器结合多普勒频移图像数据中的入射方向上的分量，确定被测量流体中的散射介质的横向流速。其测量精度直接受到多普勒角、取样定律和频率分辨率等因素的制约。在实际应用中还会由于无法确定多普勒角而导致不能进行测量。至今尚未见有降低直至消除多普勒测速方法中多普勒角影响流速测量的有效方法。

发明内容

本发明所要解决的技术问题是弥补上述现有技术的不足，提出一种基于光学相干层析成像系统的测量散射性流体横向流速的方法。

本发明的技术问题采用以下技术方案予以解决：

这种测量散射性流体横向流速的方法，基于OCT系统测量流体流速。

这种测量散射性流体横向流速的方法的特点是，包括以下步骤：

1)选择测量用在前流体截面和在后流体截面；

2)在设定的采集时间内采用扫描光束交替进行在前流体截面和在后流体截面扫描；

3)重建OCT结构图得到在前流体截面和在后流体截面位置处各自的结构图像序列；

4)将在前流体截面处的结构图像序列中某坐标(i，j)位置处的强度依次提取出，按采集时间先后排成在采集时间内通过在前流体截面坐标(i，j)位置处的不同散射颗粒的在前流体截面散射强度序列a_ij(n)；

又将在后流体截面处的结构图像序列中同为某坐标(i，j)位置处的强度依次提取出，按采集时间先后排成在采集时间内通过在后流体截面坐标(i，j)位置处的不同散射颗粒的在后流体截面散射强度序列b_ij(n)；

所述在前流体截面散射强度序列a_ij(n)和在后流体截面散射强度序列b_ij(n)的长度为正整数N＝t/T，其中t是设定的采集时间，T是扫描光束完成一次交替扫描的周期，两者可依据待测流速预估范围进行调节；

5)对在前流体截面散射强度序列a_ij(n)和在后流体截面散射强度序列b_ij(n)进行互相关运算，确定时域中的互相关函数R_abij(n)；

6)确定散射介质流过在前流体截面和在后流体截面的延迟时间；

7)确定散射介质流过在前流体截面和在后流体截面某坐标(i，j)位置处的平均横向流速测量结果；

8)改变坐标(i，j)，重复步骤4)～7)，确定散射介质在各个位置处的横向流速测量结果，即测量出流体的横向流速分布。

本发明的技术问题采用以下进一步的技术方案予以解决：

所述步骤1)的在前流体截面和在后流体截面之间的距离，即所选流体片段的长度L为几十至几百微米。

所述步骤2)的扫描光束的截面扫描频率至少是百赫兹级。

所述步骤3)是分别对所述在前流体截面采集的N幅截面数据、对所述在后流体截面采集的N幅截面数据进行OCT结构图像重建，分别得到在前流体截面位置处的结构图像序列和在后流体截面位置处的结构图像序列。

所述步骤5)的互相关运算，包括时域算法和频域算法。

所述步骤5)的时域算法依次有以下子步骤：

5-1)对a_ij(n)和b_ij(n)两组数据进行补零，得到两组新数据x(n)和y(n)：

其中：a_ij(n)、b_ij(n)、x(n)、y(n)表示为n的函数，a_ijn、b_ijn、x_n、y_n分别表示a_ij(n)、b_ij(n)、x(n)、y(n)中的第n个元素，其他标注意义类似；

5-2)计算 $h_n=\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{x^{*}}_k{y}_{n+k},$ 其中n＝-(N-1)，-(N-2)，...(N-2)，(N-1)；

其中：h_n是下标可以有负数的中间变量，没有特殊的物理或工程上的含义；

5-3)确定x(n)和y(n)即a_ij(n)和b_ij(n)的互相关函数R_abij(n)即R_xy(n)：

R_xyn＝h_n-(N-1)，R_xy(n)是将h(n)向右平移，使自变量n的取值全部变为非负整数，其中n＝0，1，2……2N-2。

所述步骤5)的频域算法依次有以下子步骤：

5′-1)对a_ij(n)和b_ij(n)两组数据进行补零，得到两组新数据x(n)和y(n)：

$x_n=\left\{\begin{array}{cc}{a}_{\mathrm{ijn}},& 0\&le;n<N\\ 0,& N\&le;n<2N\end{array}\right.$

$y_n=\left\{\begin{array}{cc}{b}_{\mathrm{ijn}},& 0\&le;n<N\\ 0,& N\&le;n<2N\end{array}\right.$

5′-2)进行傅里叶变换分解成幅值分量和频率分量

将序列x(n)进行傅里叶变换为对应的谱函数X(k)＝FFT[x(n)]；同时将序列y(n)进行傅里叶变换为对应的谱函数Y(k)＝FFT[y(n)]；

5′-3)进行共轭相乘确定互谱函数

将谱函数X(k)和谱函数Y(k)进行共轭相乘为互谱函数P_xy(k)＝X(k)×Y^*(k)；

5′-4)进行傅里叶逆变换确定时域中的互相关函数

将所述互谱函数P_xy(k)进行傅里叶逆变换为时域中的互相关函数R_xy(n)＝IFFT[P_xy(k)]，R_xy(n)也即为序列a_ij(n)和序列b_ij(n)的互相关函数R_abij(n)；

当数据长度N较大时，互相关函数的频域算法较时域算法的计算量要小。进行互相关运算的两个信号在互相关函数的峰值位置处相关程度最高。所述互相关函数R_abij(n)的峰值对应的横坐标偏离结果数据横坐标中心点(N-1，0)的距离，即表示进行互相关运算的序列a_ij(n)和序列b_ij(n)之间的时移量。

所述步骤6)是根据所述互相关函数R_abij(n)的峰值位置，并在其前后选取几个数据进行至少三次样条插值，由插值结果中的峰值位置可以得到较精确的散射介质流过所述在前流体截面A至在后流体截面B对应位置处的时延估计值τ_ij。样条插值是为了减小时延估计τ_ij的误差。

所述步骤7)确定的散射介质流过在前流体截面和在后流体截面坐标为第i行第j列位置处的平均横向流速测量结果为V_ij＝L/τ_ij。由于τ_ij对应的前后流体截面距离为L，因此，所得测量结果V_ij为流体在距离L上的平均速度。同时由于所选流体片段的长度L只有几十至几百微米，相对于整个流体的长度可以近似为一个截面，因此，V_ij可以视为该横截面上坐标为第i行第j列位置处的流速。

所述步骤8)改变坐标(i，j)，重复步骤4)～7)，确定整个流体截面各个位置处的横向流速分布。

本发明与现有技术对比的有益效果是：

本发明的横向流速测量方法不必依赖于对流体夹角进行测量，可以在未知OCT系统样品臂端口发出的扫描光束入射角的情况下测量出流体的横向流速分布，即测量与该入射角无关。

附图说明

图1是本发明具体实施方式待测量管状流体选择流体截面的示意图；

图2是某一时刻待测量管状流体内流体截面的OCT结构图即散射强度图；

图3是待测量管状流体内流体截面的横向流速分布图；

图4是在待测量管状流体OCT轴向深度为截面中央处的流速分布图；

图5是待测量管状流体截面流速分布的三维效果图。

具体实施方式

下面将结合具体实施方式并对照附图对本发明作进一步说明。

一种如图1～5所示的基于OCT系统的测量散射性流体横向流速的方法，OCT系统所用宽带光源中心波长为830nm、带宽为35nm，样品臂端口发出的扫描光束功率约为500μW，轴向扫描频率为20kHz。

包括以下步骤：

1)选择测量用在前流体截面A和在后流体截面B

在平行于流体流场方向即横向待测量的直径为0.4mm的管状流体内随机选择两个相隔距离，即所选流体片段的长度L＝0.25mm的测量用在前流体截面A和在后流体截面B，如图1所示。

2)在设定的采集时间内采用扫描光束交替进行在前流体截面A和在后流体截面B扫描

采用OCT系统样品臂端口发出的扫描光束对在前流体截面A、在后流体截面B交替进行扫描，交替扫描周期为T＝10ms、扫描持续时间即数据采集过程时间为t＝2s，交替扫描次数N＝t/T＝200，分别采集N＝200幅的在前流体截面A的截面数据和N＝200幅的在后流体截面B的截面数据。

3)重建OCT结构图得到在前流体截面A和在后流体截面B位置处各自的结构图像序列

分别将步骤2)中在前流体截面A采集的N＝200幅截面数据、在后流体截面B采集的N＝200幅截面数据进行OCT结构图像重建，得到在前流体截面A位置处的结构图像序列和在后流体截面B位置处的结构图像序列。

4)将在前流体截面A的结构图像序列中某坐标(i，j)位置处的强度依次提取出，按采集时间先后排成在采集时间内通过在前流体截面A坐标(i，j)位置处的不同散射颗粒的在前流体截面散射强度序列a_ij(n)，其长度为正整数N＝t/T＝200；

又将在后流体截面B的结构图像序列中同为某坐标(i，j)位置处的强度依次提取出，按采集时间先后排成在采集时间内通过在后流体截面B坐标(i，j)位置处的不同散射颗粒的在后流体截面散射强度序列b_ij(n)，其长度为正整数N＝t/T＝200；

例如：在前流体截面A和在后流体截面B的某同一点位置的坐标(i，j)为(7，8)，先在前流体截面A的N幅OCT结构图中，提取出每幅图中坐标(7，8)处的强度，也就是说一幅OCT结构图只提取一个数，N幅OCT结构图共提取N个数，这N个数是对应在前流体截面A同一点坐标(i，j)为(7，8)位置的N次测量，按时间先后排列成散射强度序列a₇₈(n)。

5)对在前流体截面散射强度序列a_ij(n)和在后流体截面散射强度序列b_ij(n)进行互相关运算，确定时域中的互相关函数R_abij(n)；

例如：对散射强度序列a₇₈(n)和b₇₈(n)进行互相关运算。

所述互相关运算是频域算法，依次有以下子步骤：

5-1)对a_ij(n)和b_ij(n)两组数据进行补零，得到两组新数据x(n)和y(n)：

$x_n=\left\{\begin{array}{cc}{a}_{\mathrm{ijn}},& 0\&le;n<N\\ 0,& N\&le;n<2N\end{array}\right.$

$y_n=\left\{\begin{array}{cc}{b}_{\mathrm{ijn}},& 0\&le;n<N\\ 0,& N\&le;n<2N\end{array}\right.$

5-2)进行傅里叶变换分解成幅值分量和频率分量

将序列x(n)进行傅里叶变换为对应的谱函数X(k)＝FFT[x(n)]；同时将序列y(n)进行傅里叶变换为对应的谱函数Y(k)＝FFT[y(n)]；

5-3)进行共轭相乘确定互谱函数

将谱函数X(k)和谱函数Y(k)进行共轭相乘为互谱函数P_xy(k)＝X(k)×Y^*(k)；

5-4)进行傅里叶逆变换确定时域中的互相关函数

将所述互谱函数P_xy(k)进行傅里叶逆变换为时域中的互相关函数R_xy(n)＝IFFT[P_xy(k)]，R_xy(n)也即为序列a_ij(n)和序列b_ij(n)的互相关函数R_abij(n)；

进行互相关运算的两个信号在互相关函数的峰值位置处相关程度最高，所述互相关函数R_abij(n)的峰值对应的横坐标偏离结果数据横坐标中心点(N-1，0)的距离，即表示进行互相关运算的序列a_ij(n)和序列b_ij(n)之间的时移量。

6)确定散射介质分别流过前流体截面A和后流体截面B的延迟时间

例如：确定散射颗粒流过在前流体截面A坐标(i，j)为(7，8)位置再流过在后流体截面B坐标(i，j)为(7，8)位置所用的延迟时间。

为减小延迟时间τ_ij的误差，根据所述互相关函数R_abij(n)的峰值位置，并在其前后选取几个数据进行至少三次样条插值，由插值结果中的峰值位置可以得到较精确的散射介质流过所述在前流体截面A至在后流体截面B对应位置处的延迟时间τ_ij。

7)确定的散射介质流过在前流体截面A和在后流体截面B坐标为第i行第j列位置处的平均横向流速测量结果为V_ij＝L/τ_ij。

例如：用在前流体截面A和在后流体截面B的距离L除以所用的延迟时间就是坐标(i，j)为(7，8)位置处的平均流速。

8)改变坐标(i，j)，重复步骤4)～7)，确定整个流体截面各个位置处的横向流速分布。由于所选流体片段即在前流体截面A和在后流体截面B的距离L为几十至几百微米量级，相对于整个流体的长度可以近似为一个截面，测得的横向流速分布可以视为这个流体截面的横向流速分布。

图2的某一时刻待测量管状流体内流体截面的OCT结构图即散射强度图：白色表示扫描截面的结构，其中白色颗粒为流体中的散射介质；黑色为背景；图3的待测量管状流体内流体截面的横向流速分布图，图中灰度值表示流速大小，实线圆环为等灰度值线，整个截面的流速分布为由管壁向中心逐渐变大，呈圆对称分布；图4的在待测量管状流体OCT轴向深度为截面中央处的流速分布图表明：待测散射流体的平均流速大约为3～5mm/s，且流速分布趋势基本符合描述管状通道中层流流速分布的泊肃叶定律(Poiseuille′s law)；图5的待测量管状流体截面流速分布的三维效果图更直观地反映了截面流速分布。

本具体实施方式的横向流速V的测量范围为L/(N-1)T～L/T。具体应用测量时改变相隔距离L、交替扫描周期T和交替扫描次数N的值，可以使横向流速V的测量范围与测量要求所匹配。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下做出若干等同替代或明显变型，而且性能或用途相同，都应当视为属于本发明的保护范围。

Claims (10)

1.一种测量散射性流体横向流速的方法，基于OCT系统测量流体流速，其特征在于：

包括以下步骤：

1)选择测量用在前流体截面和在后流体截面；

2)在设定的采集时间内采用扫描光束交替进行在前流体截面和在后流体截面扫描；

3)重建OCT结构图得到在前流体截面和在后流体截面位置处各自的结构图像序列；

4)将在前流体截面处的结构图像序列中某坐标(i，j)位置处的强度依次提取出，按采集时间先后排成在采集时间内通过在前流体截面坐标(i，j)位置处的不同散射颗粒的在前流体截面散射强度序列a_ij(n)；

又将在后流体截面处的结构图像序列中同为某坐标(i，j)位置处的强度依次提取出，按采集时间先后排成在采集时间内通过在后流体截面坐标(i，j)位置处的不同散射颗粒的在后流体截面散射强度序列b_ij(n)；

所述在前流体截面散射强度序列a_ij(n)和在后流体截面散射强度序列b_ij(n)的长度为正整数N＝t/T，其中t是设定的采集时间，T是扫描光束完成一次交替扫描的周期，所述设定的采集时间t以及所述扫描光束完成一次交替扫描的周期T依据待测流速预估范围进行调节；

5)对在前流体截面散射强度序列a_ij(n)和在后流体截面散射强度序列b_ij(n)进行互相关运算，确定时域中的互相关函数R_abij(n)；

6)确定散射介质流过在前流体截面和在后流体截面的延迟时间；

7)确定散射介质流过在前流体截面和在后流体截面某坐标(i，j)位置处的平均横向流速测量结果；

8)改变坐标(i，j)，重复步骤4)～7)，确定散射介质在前流体截面和在后流体截面位置处的横向流速测量结果，即测量出流体的横向流速分布。

2.如权利要求1所述的测量散射性流体横向流速的方法，其特征在于：

所述步骤1)的在前流体截面和在后流体截面之间的距离L，即所选流体片段的长度L为几十至几百微米。

3.如权利要求2所述的测量散射性流体横向流速的方法，其特征在于：

所述步骤2)的扫描光束的截面扫描频率至少是百赫兹级。

4.如权利要求3所述的测量散射性流体横向流速的方法，其特征在于：

所述步骤3)是分别对所述在前流体截面采集的N幅截面数据、对所述在后流体截面采集的N幅截面数据进行OCT结构图像重建，分别得到在前流体截面位置处的结构图像序列和在后流体截面位置处的结构图像序列。

5.如权利要求4所述的测量散射性流体横向流速的方法，其特征在于：

所述步骤5)的互相关运算，包括时域算法和频域算法。

6.如权利要求5所述的测量散射性流体横向流速的方法，其特征在于：

所述步骤5)的时域算法依次有以下子步骤：

5-1)对a_ij(n)和b_ij(n)两组数据进行补零，得到两组新数据x(n)和y(n)：

其中：a_ij(n)、b_ij(n)、x(n)、y(n)表示为n的函数，a_ijn、b_ijn、x_n、y_n分别表示a_ij(n)、b_ij(n)、x(n)、y(n)中的第n个元素，其他标注意义类似；

5-2)计算 $h_n=\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{x^{*}}_k{y}_{n+k},$ 其中n＝-(N-1)，-(N-2)，...(N-2)，(N-1)；

其中：h_n是下标可以有负数的中间变量，没有特殊的物理或工程上的含义；

5-3)确定x(n)和y(n)即a_ij(n)和b_ij(n)的互相关函数R_abij(n)即R_xy(n)：

R_xyn＝h_n-(N-1)，R_xy(n)是将h(n)向右平移，使自变量n的取值全部变为非负整数，其中n＝0，1，2……2N-2。

7.如权利要求6所述的测量散射性流体横向流速的方法，其特征在于：

所述步骤5)的频域算法依次有以下子步骤：

5′-1)对a_ij(n)和b_ij(n)两组数据进行补零，得到两组新数据x(n)和y(n)：

$x_n=\left\{\begin{array}{cc}{a}_{\mathrm{ijn}},& 0\&le;n<N\\ 0,& N\&le;n<2N\end{array}\right.$

$y_n=\left\{\begin{array}{cc}{b}_{\mathrm{ijn}},& 0\&le;n<N\\ 0,& N\&le;n<2N\end{array}\right.$

5′-2)进行傅里叶变换分解成幅值分量和频率分量

将序列x(n)进行傅里叶变换为对应的谱函数X(k)＝FFT[x(n)]；同时将序列y(n)进行傅里叶变换为对应的谱函数Y(k)＝FFT[y(n)]；

5′-3)进行共轭相乘确定互谱函数

将谱函数X(k)和谱函数Y(k)进行共轭相乘为互谱函数P_xy(k)＝X(k)×Y^*(k)；

5′-4)进行傅里叶逆变换确定时域中的互相关函数

将所述互谱函数P_xy(k)进行傅里叶逆变换为时域中的互相关函数R_xy(n)＝IFFT[P_xy(k)]，R_xy(n)也即为序列a_ij(n)和序列b_ij(n)的互相关函数R_abij(n)；

所述互相关函数R_abij(n)的峰值对应的横坐标偏离结果数据横坐标中心点(N-1，0)的距离，即表示进行互相关运算的序列a_ij(n)和序列b_ij(n)之间的时移量。

8.如权利要求6或7所述的测量散射性流体横向流速的方法，其特征在于：

所述步骤6)是为根据所述互相关函数R_abij(n)的峰值位置，并在其前后选取几个数据进行至少三次样条插值，由插值结果中的峰值位置可以得到较精确的散射介质流过所述在前流体截面至在后流体截面对应位置处的时延估计值τ_ij。

9.如权利要求8所述的测量散射性流体横向流速的方法，其特征在于：

所述步骤7)确定的散射介质流过在前流体截面和在后流体截面坐标为第i行第j列位置处的平均横向流速测量结果为V_ij＝L/τ_ij。

10.如权利要求9所述的测量散射性流体横向流速的方法，其特征在于：

所述步骤8)改变坐标(i，j)，重复步骤4)～7)，确定整个流体截面各个位置处的横向流速分布。

Images