CN111369620A - 一种基于改进piv的水平油水两相流流速测量方法 - Google Patents

Abstract

本发明属于石油工程与测量技术领域，具体涉及一种基于改进PIV的水平油水两相流流速测量方法，1、对连续两帧图像进行处理，将图像信息提取转换为数字信息，确定PIV位移测量过程中查询窗口M；2、固定查询窗口大小，计算两帧图像中查询窗口内中心处粒子的互相关系数

提取最大的互相关系数及对应的位移整像素值(p,q)；3、建立互相关函数R(u,v)，确定粒子的位移亚像素值(u₀,v₀)；4、改变查询窗口大小，对空白位移矢量处进行插值以扩大位移矢量分布；5、重复步骤二～四，直至查询窗口循环结束为止，计算所有位移矢量对应的像素点的速度矢量分布V，最终求得的平均值即为油水两相流在此时的平均流速

Description

一种基于改进PIV的水平油水两相流流速测量方法

技术领域：

本发明属于石油工程与测量技术领域，具体涉及一种基于改进PIV的水平 油水两相流流速测量方法。

背景技术：

目前，我国绝大多数油田开发已进入高含水、低产液阶段。油水两相流的 流型表现为油水分层，且上层含少量油泡的混合状态，此状态易导致油水滑脱。 水平井、小管径以及化学驱等新技术实施后，油水流动状态受这种滑移影响较 大，给测井资料的解释和研究开发中测井仪器的评价带来困难。测井资料的解 释和测井仪器的精度很大程度依赖于流动特性，流型、流场、流速和组分是反 映流动特性的参数，其中流速测量是油田开发动态监测必不可少的一部分，为 油田开采流量计量、提高原油采收率等提供重要的参考依据，故水平井油水两 相流的流速准确测量是目前的急需解决的问题。

PIV算法是常用的一种流体的流速测量方法，不仅可以定量分析整个管道 流动的速度矢量场，而且可以消除传统单点测量的局限性。经典PIV算法是先 独立对水平方向和垂直方向进行位移亚像素拟合再将其结合在一起得到某一点 的精确位移，这一解决方法得到的位移亚像素拟合效果在单个方向是最好的， 但对于某一点而言未必是最好的，因此此方法对位移的计算精度不高。目前油 水两相流PIV图像位移空白矢量处是采用线性插值得到的，故其流体流速的精 度也会降低。针对这两个现象，若将位移亚像素拟合从一维扩大到二维，精确 计算每点的位移变化，同时采用Kriging模型进行位移空白矢量处插值，可提 高位移矢量的精确度，从而降低平均流速测量的误差。

发明内容：

本发明的目的是为了解决采用PIV技术测量水平油水两相流流速过程中， 由于PIV算法过程中空白矢量和分向亚像素拟合造成的测量精度不高问题，提 出一种基于改进PIV的水平油水两相流流速测量方法。

本发明采用的技术方案为：一种基于改进PIV的水平油水两相流流速测量 方法，其特征在于：所述测量方法包括以下步骤：

步骤一：对连续两帧图像进行处理，将图像信息提取转换为数字信息，确 定PIV位移测量过程中查询窗口M；

所述的图像信息提取转换是将两帧图像的有用信息进行裁剪，并进行灰度 处理，转换后的两帧灰度图像函数分别用f(x,y)和g(x,y)表示，x,y分别表示图 像中水平方向和垂直方向上的像素点位置坐标；

查询窗口M的设置可提高位移的分辨率和计算精度，M表示为

式中m₁～m₈均为后续互相关时查询窗口大小，m₁,m₂,m₃,m₄,m₅,m₆,m₇,m₈的查询窗口大小分别为128×128、128×128、64×64、64×64、32×32、32×32、 16×16、16×16；

步骤二：固定查询窗口大小，计算两帧图像中查询窗口内中心处粒子的互 相关系数

提取最大的互相关系数以及对应的位移整像素值(p,q)；

所述两帧图像查询窗口内中心处的相关系数

的计算公式为：

式中矩阵P,Q分别为第一帧图像查询窗口中心与第二帧图像查询窗口中 心在水平方向和垂直方向的位移差；

根据互相关性质可知，

值越大，两帧图像查询窗口中心越相关；因 此，最大的互相关系数对应的位移差即为查询窗口中心的位移；查询窗口中心 位移整像素大小的计算公式为：

式中p,q分别为查询窗口中心移动的水平方向和垂直方向的位移整像素 值；

步骤三：根据步骤二中提取最大互相关系数及对应的位移整像素值建立互 相关函数R(u,v)，利用最大互相关系数周围与其相邻四个互相关系数及对应的 位移整像素值确定粒子的位移亚像素值(u₀,v₀)；

所述互相关函数R(u,v)的建立方式如下：

因采集的油水两相流粒子图像属于高斯分布，故计算的互相关函数R(u,v) 亦满足高斯分布，互相关函数R(u,v)的计算公式为

式中A为相关峰峰值，σ_u ²,σ_v ²分别为相关峰在水平方向和垂直方向上的方 差，u,v分别为查询窗口中心移动的水平方向和垂直方向的位移矢量，u₀,v₀分 别为查询窗口中心移动的水平方向和垂直方向的精确位移矢量；

根据步骤二确定最大互相关系数为R(p,q)，对应的粒子在水平方向和垂直 方向上位移整像素值分别为p,q，以及周围与其相邻四个相关系数R(p-1,q), R(p+1,q),R(p,q-1),R(p,q+1)，将其带入互相关函数R(u,v)表达式中：

对上述方程组进行处理，消去参数A,σ_u ²,σ_v ²以获取精确的位移亚像素值：

步骤四：改变查询窗口大小，利用Kriging模型对空白位移矢量处进行插 值以扩大位移矢量分布；

改变查询窗口大小时，查询窗口内部分粒子的位移矢量是空白的；该粒子 位置坐标为(x_c,y_c)，且该处的位移矢量用

表示，图像中已知位移矢量的粒子位 置坐标为(x_n,y_n)(n＝1,2,…,h)，且此处的位移矢量用Z_n(n＝1,2,…,h)表示，h为已知 位移矢量的粒子个数，n为已知位移矢量的粒子序号；

Kriging模型的计算方法为：

式中λ_n是指已知位移矢量Z_n对应的权重系数；

根据Kriging模型的无偏最优性，采用拉格朗日算法估算λ_i：

式中F为一拉格朗日量，k,l均为已知位移矢量的粒子序号，λ_k,λ_l分别是 指已知位移矢量Z_k和Z_l对应的权重系数，μ为拉格朗日乘子，γ_kl为位置坐标为 (x_k,y_k)(k＝1,2,…,h)和(x_l,y_l)(l＝1,2,…,h)的粒子间的半方差函数，γ_kc为位置坐标为 (x_k,y_k)(k＝1,2,…,h)和(x_c,y_c)的粒子间的半方差函数，γ_cc为位置坐标为(x_c,y_c)粒子与 其自身间的半方差函数；因图像分布满足高斯分布，半方差函数γ_kl的计算公式 为

式中C₀为粒子与其自身的的半方差，C为随着粒子位置间距离增大基本达 到恒定的半方差值，α为具有相关性的位移矢量对应的位置距离范围；

对F求偏导得到各个粒子对应的权重系数：

步骤五：重复步骤二～四，直至查询窗口循环结束为止，计算所有位移矢 量对应的像素点的速度矢量分布V，最终求得的平均值即为油水两相流在此时 的平均流速v；

所述各粒子的移动速度计算公式为：

式中矩阵V为图像中粒子的移动速度矢量分布，其中

分别为图像 中粒子(x_wβ,y_wβ)处在水平方向和垂直方向的移动速度，i,j分别为图像中粒子在 水平方向和垂直方向的序号，Δx_ij,Δy_ij分别为图像中粒子(x_ij,y_ij)处在水平方向和 垂直方向的移动距离，Δt为两帧图像采集的时间差，β,w分别为图像中粒子在 水平方向和垂直方向的个数；

最后求得图像中油水两相流的平均流速

本发明的有益效果：采用二维位移亚像素拟合方法克服了独立对水平方向 和垂直方向的亚像素拟合造成的低精度，改进了PIV算法，同时采用Kriging 模型对位移进行插值，使得整体的速度矢量分布更为平滑且精确，提高了水平 油水两相流流速测量精度。

附图说明：

图1是实施例一中改进PIV算法的流程图；

图2是实施例一中采用经典PIV算法获取的速度矢量图；

图3图2中围框部分的局部放大图；

图4是实施例一中采用改进PIV算法获取的速度矢量图；

图5是图4中围框部分的局部放大图；

图6是实施例一中采用经典PIV算法测量不同流速含水的相对误差分析图；

图7是实施例一中采用改进PIV算法测量不同流速含水的相对误差分析图。

具体实施方式：

实施例一

一种基于改进PIV的水平油水两相流流速测量方法，其特征在于：所述测 量方法包括以下步骤：

步骤一：对连续两帧图像进行处理，将图像信息提取转换为数字信息，确 定PIV位移测量过程中查询窗口M；

所述的图像信息提取转换是将两帧图像的有用信息进行裁剪，并进行灰度 处理，转换后的两帧灰度图像函数分别用f(x,y)和g(x,y)表示，x,y分别表示图 像中水平方向和垂直方向上的像素点位置坐标；

查询窗口M的设置可提高位移的分辨率和计算精度，M表示为

式中m₁～m₈均为后续互相关时查询窗口大小，m₁,m₂,m₃,m₄,m₅,m₆,m₇,m₈的查询窗口大小分别为128×128、128×128、64×64、64×64、32×32、32×32、 16×16、16×16；

步骤二：固定查询窗口大小，计算两帧图像中查询窗口内中心处粒子的互 相关系数

提取最大的互相关系数以及对应的位移整像素值(p,q)；

所述两帧图像查询窗口内中心处的相关系数

的计算公式为：

式中矩阵P,Q分别为第一帧图像查询窗口中心与第二帧图像查询窗口中 心在水平方向和垂直方向的位移差；

根据互相关性质可知，

值越大，两帧图像查询窗口中心越相关；因 此，最大的互相关系数对应的位移差即为查询窗口中心的位移；查询窗口中心 位移整像素大小的计算公式为：

式中p,q分别为查询窗口中心移动的水平方向和垂直方向的位移整像素 值；

步骤三：根据步骤二中提取最大互相关系数及对应的位移整像素值建立互 相关函数R(u,v)，利用最大互相关系数周围与其相邻四个互相关系数及对应的 位移整像素值确定粒子的位移亚像素值(u₀,v₀)；

所述互相关函数R(u,v)的建立方式如下：

因采集的油水两相流粒子图像属于高斯分布，故计算的互相关函数R(u,v) 亦满足高斯分布，互相关函数R(u,v)的计算公式为

式中A为相关峰峰值，σ_u ²,σ_v ²分别为相关峰在水平方向和垂直方向上的方 差，u,v分别为查询窗口中心移动的水平方向和垂直方向的位移矢量，u₀,v₀分 别为查询窗口中心移动的水平方向和垂直方向的精确位移矢量；

根据步骤二确定最大互相关系数为R(p,q)，对应的粒子在水平方向和垂直 方向上位移整像素值分别为p,q，以及周围与其相邻四个相关系数R(p-1,q), R(p+1,q),R(p,q-1),R(p,q+1)，将其带入互相关函数R(u,v)表达式中：

对上述方程组进行处理，消去参数A,σ_u ²,σ_v ²以获取精确的位移亚像素值：

步骤四：改变查询窗口大小，利用Kriging模型对空白位移矢量处进行插 值以扩大位移矢量分布；

改变查询窗口大小时，查询窗口内部分粒子的位移矢量是空白的；该粒子 位置坐标为(x_c,y_c)，且该处的位移矢量用

表示，图像中已知位移矢量的粒子位 置坐标为(x_n,y_n)(n＝1,2,…,h)，且此处的位移矢量用Z_n(n＝1,2,…,h)表示，h为已知 位移矢量的粒子个数，n为已知位移矢量的粒子序号；

Kriging模型的计算方法为：

式中λ_n是指已知位移矢量Z_n对应的权重系数；

根据Kriging模型的无偏最优性，采用拉格朗日算法估算λ_i：

式中F为一拉格朗日量，k,l均为已知位移矢量的粒子序号，λ_k,λ_l分别是 指已知位移矢量Z_k和Z_l对应的权重系数，μ为拉格朗日乘子，γ_kl为位置坐标为 (x_k,y_k)(k＝1,2,…,h)和(x_l,y_l)(l＝1,2,…,h)的粒子间的半方差函数，γ_kc为位置坐标为 (x_k,y_k)(k＝1,2,…,h)和(x_c,y_c)的粒子间的半方差函数，γ_cc为位置坐标为(x_c,y_c)粒子与 其自身间的半方差函数；因图像分布满足高斯分布，半方差函数γ_kl的计算公式 为

式中C₀为粒子与其自身的的半方差，C为随着粒子位置间距离增大基本达 到恒定的半方差值，α为具有相关性的位移矢量对应的位置距离范围；

对F求偏导得到各个粒子对应的权重系数：

步骤五：重复步骤二～四，直至查询窗口循环结束为止，计算所有位移矢 量对应的像素点的速度矢量分布V，最终求得的平均值即为油水两相流在此时 的平均流速

所述各粒子的移动速度计算公式为：

式中矩阵V为图像中粒子的移动速度矢量分布，其中

分别为图像 中粒子(x_wβ,y_wβ)处在水平方向和垂直方向的移动速度，i,j分别为图像中粒子在 水平方向和垂直方向的序号，Δx_ij,Δy_ij分别为图像中粒子(x_ij,y_ij)处在水平方向和 垂直方向的移动距离，Δt为两帧图像采集的时间差，β,w分别为图像中粒子在 水平方向和垂直方向的个数；

最后求得图像中油水两相流的平均流速

将本方法应用于水平油水两相流中测量油水两相流流速过程中，并与经典 PIV算法进行了对比，两种算法求取的速度矢量图如图2-图5所示，图2和图 3均由matlab软件编程绘制获取得到的，如果使用灰度图，只能直观看出速度 矢量图中的速度方向分布，看不出速度矢量图中的速度大小分布，使用颜色图 可清晰明了地观察到速度矢量图中各处的速度大小及方向分布。图中颜色代表 速度绝对值的大小，颜色由蓝变黄再变红，速度绝对值逐渐变大。此外，速度 矢量图中箭头的起始点和终点分别表示油水两相流的起始点和终点。分析图2- 图5可知，图3中速度矢量相比于图2在方向上和大小上变化较为平缓，使整体速度矢量平滑且稳定。

为降低偶然性，进行多组流速测量实验，采用经典PIV算法和改进PIV算 法最终获取的平均流速对应的相对误差分析分别如图6和图7所示。经计算可 知，采用经典PIV算法求不同流速含水状态下平均流速的平均相对误差分别为 12.45％，采用改进PIV算法求不同流速含水状态下平均流速的平均相对误差分 别为2.68％，远低于经典PIV算法，均在测量允许误差范围内，满足油田测量 要求。本方法采用改进PIV算法测量油水两相流平均流速，获取的速度矢量图 平滑稳定，表明该改进PIV算法计算精度高且稳定，具有很强的实用性，实现 水平油水两相流流速准确测量。

固定查询窗口，将采集到的两帧图像进行位移互相关，获取其最大的位移 整像素值；然后改进PIV算法中的亚像素拟合方法，采用二维位移亚像素拟合 进行高斯拟合确定粒子的精确位移；改变查询窗口，采用Kriging插值将位移 空白矢量处进行插值，改变了经典PIV中插值算法；根据获取的精确位移和采 集图像间的时间差计算粒子的速度从而确定水平油水两相流的平均流速，为油 田生产开发提供有效的支撑。

以上内容是结合具体的实施方式对本发明作的进一步详细说明，不能认定 本发明只局限于上述具体实施。在不脱离本发明整体思路和权利要求所保护的 前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (1)

1.一种基于改进PIV的水平油水两相流流速测量方法，其特征在于：所述测量方法包括以下步骤：

步骤一：对连续两帧图像进行处理，将图像信息提取转换为数字信息，确定PIV位移测量过程中查询窗口M；

所述的图像信息提取转换是将两帧图像的有用信息进行裁剪，并进行灰度处理，转换后的两帧灰度图像函数分别用f(x,y)和g(x,y)表示，x,y分别表示图像中水平方向和垂直方向上的像素点位置坐标；

查询窗口M的设置可提高位移的分辨率和计算精度，M表示为

式中m₁～m₈均为后续互相关时查询窗口大小，m₁,m₂,m₃,m₄,m₅,m₆,m₇,m₈的查询窗口大小分别为128×128、128×128、64×64、64×64、32×32、32×32、16×16、16×16；

步骤二：固定查询窗口大小，计算两帧图像中查询窗口内中心处粒子的互相关系数

提取最大的互相关系数以及对应的位移整像素值(p,q)；

所述两帧图像查询窗口内中心处的相关系数

的计算公式为：

式中矩阵P,Q分别为第一帧图像查询窗口中心与第二帧图像查询窗口中心在水平方向和垂直方向的位移差；

根据互相关性质可知，

值越大，两帧图像查询窗口中心越相关；因此，最大的互相关系数对应的位移差即为查询窗口中心的位移；查询窗口中心位移整像素大小的计算公式为：

式中p,q分别为查询窗口中心移动的水平方向和垂直方向的位移整像素值；

步骤三：根据步骤二中提取最大互相关系数及对应的位移整像素值建立互相关函数R(u,v)，利用最大互相关系数周围与其相邻四个互相关系数及对应的位移整像素值确定粒子的位移亚像素值(u₀,v₀)；

所述互相关函数R(u,v)的建立方式如下：

因采集的油水两相流粒子图像属于高斯分布，故计算的互相关函数R(u,v)亦满足高斯分布，互相关函数R(u,v)的计算公式为

式中A为相关峰峰值，σ_u ²,σ_v ²分别为相关峰在水平方向和垂直方向上的方差，u,v分别为查询窗口中心移动的水平方向和垂直方向的位移矢量，u₀,v₀分别为查询窗口中心移动的水平方向和垂直方向的精确位移矢量；

根据步骤二确定最大互相关系数为R(p,q)，对应的粒子在水平方向和垂直方向上位移整像素值分别为p,q，以及周围与其相邻四个相关系数R(p-1,q),R(p+1,q),R(p,q-1),R(p,q+1)，将其带入互相关函数R(u,v)表达式中：

对上述方程组进行处理，消去参数A,σ_u ²,σ_v ²以获取精确的位移亚像素值：

步骤四：改变查询窗口大小，利用Kriging模型对空白位移矢量处进行插值以扩大位移矢量分布；

改变查询窗口大小时，查询窗口内部分粒子的位移矢量是空白的；该粒子位置坐标为(x_c,y_c)，且该处的位移矢量用

表示，图像中已知位移矢量的粒子位置坐标为(x_n,y_n)(n＝1,2,…,h)，且此处的位移矢量用Z_n(n＝1,2,…,h)表示，h为已知位移矢量的粒子个数，n为已知位移矢量的粒子序号；

Kriging模型的计算方法为：

式中λ_n是指已知位移矢量Z_n对应的权重系数；

根据Kriging模型的无偏最优性，采用拉格朗日算法估算λ_i：

式中F为一拉格朗日量，k,l均为已知位移矢量的粒子序号，λ_k,λ_l分别是指已知位移矢量Z_k和Z_l对应的权重系数，μ为拉格朗日乘子，γ_kl为位置坐标为(x_k,y_k)(k＝1,2,…,h)和(x_l,y_l)(l＝1,2,…,h)的粒子间的半方差函数，γ_kc为位置坐标为(x_k,y_k)(k＝1,2,…,h)和(x_c,y_c)的粒子间的半方差函数，γ_cc为位置坐标为(x_c,y_c)粒子与其自身间的半方差函数；因图像分布满足高斯分布，半方差函数γ_kl的计算公式为

式中C₀为粒子与其自身的的半方差，C为随着粒子位置间距离增大基本达到恒定的半方差值，α为具有相关性的位移矢量对应的位置距离范围；

对F求偏导得到各个粒子对应的权重系数：

步骤五：重复步骤二～四，直至查询窗口循环结束为止，计算所有位移矢量对应的像素点的速度矢量分布V，最终求得的平均值即为油水两相流在此时的平均流速

所述各粒子的移动速度计算公式为：

式中矩阵V为图像中粒子的移动速度矢量分布，其中

分别为图像中粒子(x_wβ,y_wβ)处在水平方向和垂直方向的移动速度，i,j分别为图像中粒子在水平方向和垂直方向的序号，Δx_ij,Δy_ij分别为图像中粒子(x_ij,y_ij)处在水平方向和垂直方向的移动距离，Δt为两帧图像采集的时间差，β,w分别为图像中粒子在水平方向和垂直方向的个数；

最后求得图像中油水两相流的平均流速

Images