阅读说明：本技术 基于改进的dpiv垂直井油水两相流流速场测量方法 (Improved DPIV vertical well-based oil-water two-phase flow velocity field measurement method ) 是由 韩连福 丛垚 付长凤 刘超 刘兴斌 姜继玉 于 2019-12-09 设计创作，主要内容包括：本发明涉及的是基于改进的DPIV垂直井油水两相流流速场测量方法,这种基于改进的DPIV垂直井油水两相流流速场测量方法采用迭代最近点(Iterative Closest Point,ICP)代替传统互相关进行图像匹配,同时利用移动最小二乘(Moving least squares,MLS)综合整个流速场信息进行边界位移值补充,改善了DPIV的图像匹配效果,提高了垂直井油水两相流的测量精度。本发明在匹配的过程中同时考虑了平面内油滴的平移和旋转,能够提高两相流灰度图像的匹配精度,解决了在垂直井油水两相流流速场测量中由于图像匹配效果差、位移场边界缺失导致的测量精度降低问题。(The invention relates to an improved DPIV-based vertical well oil-water two-phase flow velocity field measurement method, which adopts Iterative Closest Point (ICP) to replace the traditional cross correlation for image matching, and meanwhile, uses Moving Least Squares (MLS) to synthesize the whole flow velocity field information for boundary displacement value supplement, thereby improving the image matching effect of DPIV and improving the measurement precision of the vertical well oil-water two-phase flow. The invention simultaneously considers the translation and rotation of oil drops in a plane in the matching process, can improve the matching precision of two-phase flow gray images, and solves the problem of measurement precision reduction caused by poor image matching effect and displacement field boundary loss in the measurement of the oil-water two-phase flow velocity field of the vertical well.)

基于改进的DPIV垂直井油水两相流流速场测量方法

技术领域：

本发明涉及的是石油工程与测量技术领域，具体涉及的是基于改进的DPIV垂直井油水两相流流速场测量方法。

背景技术：

流速测量是油田开发动态监测必不可少的一部分，油水两相流流速测量方法有超声法、电容法、电导法、电磁法、热示踪法以及光纤法等。当前时期的油田流速测量更倾向于区域流速分布测量，因为从区域流速分布中更容易观察多相流流动特性，从而指导油田开发。传统的流速测量方法已不能满足这一要求，须采取新的流速测量方法。DPIV是流体力学领域里一种有效的全场测速方法，其无扰，非接触，全场测速的特点能够很好地满足油田多相流测速要求，使其在油田测速中拥有可观的应用前景。DPIV技术测量精度的影响因素总体分为两大类：一类是外部测量环境，如示踪粒子分布、光照分布以及各种硬件系统的影响，另一类是DPIV图像后处理技术的实施。窗口迭代变形(Window deformation iterativemultigrid,WIDIM)算法用于提高查询窗口在流场存在大速度梯度示踪粒子情况下的跟随性，是目前广泛使用的DPIV算法。

当前在125mm垂直井油水两相流测量中应用基于WIDIM的DPIV算法测量油水两相流流场存在以下问题：(1)由于油井管径较大导致查询窗口内油滴间相对运动复读较大，对查询窗口的跟随性有更高的要求；(2)由于DPIV技术应用对象为两相流且采用LED照明、油滴示踪，待处理的两相流灰度图像相较于单相流灰度图像信息更复杂，使互相关灰度计算容易出现错误匹配带来测量误差，这些误差会在迭代过程中出现“放大效应”，严重影响图像变形效果，所以迭代计算的精度应提高；(3)在WIDIM算法的插值过程中存在位移场边界缺失的问题，同样对窗口变形效果产生影响。迭代最近点(Iterative Closest Point,ICP)是一种点云数据配准方法，能够实现大量数据的精确配准，将三维ICP技术应用于二维DPIV图像配准中，只需将DPIV灰度图像投影至三维空间，相较于互相关匹配，ICP配准受噪声影响小，且在配准过程中考虑了油滴的平面旋转运动。对缺失的位移场边界，应综合考虑整个位移场信息进行赋值，因此采用移动最小二乘(Moving least squares,MLS)法对DPIV算法进行进一步改进，在每次插值前先对当前位移场进行MLS曲面拟合，根据曲面拟合函数补充边界缺失值。

发明内容

：

本发明的目的是提供基于改进的DPIV垂直井油水两相流流速场测量方法，这种基于改进的DPIV垂直井油水两相流流速场测量方法用于解决基于WIDIM的DPIV算法在垂直井油水两相流测量中存在的图像匹配效果差、位移场边界缺失产生的流速测量不准确的问题。

本发明解决其技术问题所采用的技术方案是：这种基于改进的DPIV垂直井油水两相流流速场测量方法，包括以下步骤：

步骤一：选取时间间隔为Δt的两帧油水两相流图像，图像大小表示为：长×宽，设图像大小为M pixel×N pixel，M为图像长的值，N为图像宽的值，对其进行图像去噪、图像对比度增强后，确定初始查询窗口大小，初始查询窗口大小表示为：长×宽，设初始查询窗口大小为W pixel×W pixel，W为查询窗口长、宽的值，查询步长为W/2pixel，将两帧油水两相流图像按照相应初始查询窗口大小划分为覆盖率为50％的(2M-W)/W×(2N-W)/W个网格，设定当前查询窗口，令当前查询窗口大小等于初始查询窗口大小；

步骤二：在两帧图像中进行区域选择，选定区域表示为：[X坐标范围起点：X坐标范围终点，Y坐标范围起点：Y坐标范围终点]，设在第一帧图像中选定查询区域Area₁＝[i:i+W-1,j:j+W-1]，在第二帧图像中选定查询区域Area₂＝[i:i+W-1,j:j+W-1]，i，j分别表示图像中的X，Y坐标值，i＝1+W(n-1)，j＝1+W(m-1)，n＝1,2,…,(2M-W)/W，m＝1,2,…,(2N-W)/W，m为X方向查询区域序号，n为Y方向查询区域序号；对查询区域Area₁和查询区域Area₂两个区域进行ICP配准，得到查询区域Area₁平均X方向位移u_(i',j')，平均Y方向位移v_(i',j')，i'＝i+W/2为查询区域Area₁中心X坐标，j'＝j+W/2为查询区域Area₁中心Y坐标；

步骤三：当前查询窗口以步长W/2pixel遍历两帧图像，得到初始X方向位移场U_initial、初始Y方向位移场V_initial，设当前X方向位移场为U，当前Y方向位移场为V；

步骤四：根据当前X方向位移场U、当前Y方向位移场V对当前查询窗口变形，重新进行ICP配准得到二次迭代X方向位移场U_new、二次迭代Y方向位移场V_new；

步骤五：更新当前X方向位移场U、当前Y方向位移场V，位移场按下式进行更新：

U＝U_initial+U_new，V＝V_initial+V_new；

步骤六：采用MLS对当前X方向位移场U、当前Y方向位移场V进行拟合，得到X方向边缘补充位移场U′、Y方向边缘补充位移场V′；X方向边缘补充位移场U′、Y方向边缘补充位移场V'按下述方法计算：

X方向边缘补充位移场U′曲面拟合函数f_u(x,y)、Y方向边缘补充位移场V'曲面拟合函数f_v(x,y)，x为X方向坐标变量，y为Y方向坐标变量，表示为k为多项式序号：

其中曲面拟合函数f_u(x,y)的系数阵列α_u(x,y)＝[α_u1(x,y),α_u2(x,y),…,α_uk(x,y)]，α_uk(x,y)为曲面拟合函数f_u(x,y)的第k个系数，曲面拟合函数f_v(x,y)的系数阵列α_v(x,y)＝[α_v1(x,y),α_v2(x,y),…,α_vk(x,y)]，α_vk(x,y)为曲面拟合函数f_v(x,y)的第k个系数，变量阵列μ(x,y)＝[μ₁(x,y),μ₂(x,y),...,μ_k(x,y)]＝[1,x,y,x²,xy,y²]，μ_k(x,y)为两个曲面拟合函数第k个变量，T代表矩阵转置符号；

α_u(x,y)、α_v(x,y)按下式求：

式中已知X方向位移阵列Z_u＝[u_{(W/2+1,W/2+1)},u_{(3W/2+1,3W/2+1)},…,u_(i',j')]，i'＝1+W(n-1/2)，j'＝1+W(m-1/2)，n＝1,2,…,(2M-W)/W，m＝1,2,…,(2N-W)/W，已知Y方向位移阵列Z_v＝[v_{(W/2+1,W/2+1)},v_{(3W/2+1,3W/2+1)},…,v_(i',j')]，i'＝1+W(n-1/2)，j'＝1+W(m-1/2)，n＝1,2,…,(2M-W)/W，m＝1,2,…,(2N-W)/W，参数阵列G＝[μ^T(W/2+1,W/2+1),μ^T(3W/2+1,3W/2+1),…,μ^T(i',j')]，i'＝1+W(n-1/2)，j'＝1+W(m-1/2)，n＝1,2,…,(2M-W)/W，m＝1,2,…,(2N-W)/W，权值对角矩阵为具有紧支集特性的权函数；

将α_u(x,y)、α_v(x,y)代入f_u(x,y)、f_v(x,y)中，得到X方向边缘补充位移场U'和Y方向边缘补充位移场V'拟合曲面方程：

此时X方向边缘补充位移场U'和Y方向边缘补充位移场V'表示为：

此时当前X方向位移场U＝U′，当前Y方向位移场V＝V'；

步骤七：采用双三次均匀B样条插值扩大当前X方向位移场U、当前Y方向位移场V为原来的4倍；

步骤八：将当前查询窗口大小缩小为原来的1/4，得到新尺寸查询窗口，令当前查询窗口大小等于新尺寸查询窗口大小，根据当前X方向位移场U、当前Y方向位移场V对当前查询窗口变形，重新进行ICP配准得到三次迭代X方向位移场U′_new、三次迭代Y方向V′_new；

步骤九：更新当前X方向位移场U、当前Y方向位移场V；

步骤十：对步骤九处理后的当前查询窗口迭代进行步骤四-步骤九，直至缩小到指定查询窗口的大小，确定最终X方向位移场U_final、最终Y方向位移场V_final，根据时间间隔Δt和位移场U_final、V_final得到流速场f；油水两相流速度场按下式计算：

上述方案中步骤二中进行ICP配准的方法：设置W＝128，初始查询窗口尺寸128pixel×128pixel，查询区域Area₁＝[i:i+127,j:j+127]、查询区域Area₂＝[i:i+127,j:j+127]，设查询区域Area₁的数据点集为p＝{p_ii|ii＝1,2,….,16384}，ii为查询区域Area₁的数据点序号，p_ii为数据点集p中第ii个数据点，查询区域Area₂的数据点集为q＝{q_jj|jj＝1,2,….,16384}，jj为查询区域Area₂的数据点序号，q_jj为数据点集q中第jj个数据点，该区域平均X方向位移u_(i',j')，平均Y方向位移v_(i',j')按如下方法计算：

数据点集变换关系式q'为：

q′＝rp+t

匹配目标函数E为：

其中r为旋转矩阵，t为平移向量，q'_jj为数据点集p经数据点集变换关系式q'变换后的数据点集中第jj个数据点；

采用SVD求解旋转矩阵r和平移向量t，对点集p、q进行如下变换，P_ii为p_ii经变换后的数据，Q_jj为q_jj经变换后的数据：

有：

对最优解矩阵H进行奇异值分解，将H分解为左奇异矩阵D，右奇异矩阵L，和奇异值矩阵Λ，T代表矩阵转置符号：

H＝DΛL^T

r和t的计算公式如下：

该区域的平均位移为：

u_(i′,j′)＝t(1)，v_(i′,j′)＝t(2)。

上述方案步骤三中，初始X方向位移场U_initial，初始Y方向位移场V_initial表示为:

其中i'＝1+128(n-1/2)，j'＝1+128(m-1/2)，n＝1,2,…,(2M-128)/128，m＝1,2,…,(2N-128)/128；当前X方向位移场U＝U_initial，当前Y方向位移场V＝V_initial。

上述方案步骤四中当前查询窗口进行变形后，第一帧图像的查询区域为Area₁＝[i:i+W-1,j:j+W-1]，第二帧图像的查询区域为Area₂＝[i+u(i',j'):i+W-1+u(i',j'),j+v(i',j'):j+W-1+v(i',j')]。

上述方案步骤七中双三次均匀B样条插值按下式计算：

式中Θ_u为当前X方向位移场U的双三次均匀B样条曲面片，Θ_v为当前Y方向位移场V的双三次均匀B样条曲面片，X方向坐标变量矩阵λ＝[ξ³ξ²ξ1]，ξ为X坐标变量，Y方向坐标变量矩阵γ＝[ζ³ζ²ζ1]，ζ为Y坐标变量，B_w为B样条基函数，C_u为当前X方向位移场U的双三次均匀B样条曲面片的控制顶点，C_v为当前Y方向位移场V的双三次均匀B样条曲面片的控制顶点，当X坐标变量ξ、Y坐标变量ζ在[0,1]之间遍历时，得到曲面片上任意一点的位移值。

上述方案步骤八中当前查询窗口变形后，第一帧图像的查询区域Area₁＝[i:i+W-1,j:j+W-1]，第二帧图像的查询区域Area₂＝[i+u(i',j'):i+W-1+u(i',j'),j+v(i',j'):j+W-1+v(i',j')]。

上述方案步骤九中位移场按下式更新：

U＝U′+U′_new，V＝V′+V′_new。

本发明具有以下有益效果：

(1)本发明在测量过程中使用LED背光光源代替传统激光片光源，同时采用油滴作为示踪粒子，避免了传统示踪剂在油水两相流中受油滴遮挡影响测量的问题，提高了流场和流速测量精度；

(2)本发明采用ICP代替互相关算法进行图像匹配，在匹配的过程中同时考虑了平面内油滴的平移和旋转，能够提高两相流灰度图像的匹配精度；

(3)本发明采用MLS算法综合整个流场的信息，对缺失的位移场边界值进行补充，提高基于WIDIM的DPIV算法的查询窗口跟随性，提高测量精度。

附图说明：

图1为图像坐标示意图。

图2为查询窗口遍历图像示意图，查询窗口由左至右，从上到下遍历整幅图像。

图3采用改进的DPIV算法测得的垂直井油水两相流流速场。

图4采用改进的DPIV算法测得的垂直井油水两相流流流速，并与经典DPIV测量结果做了对比，(a)为在含水率为90％工况下垂直井油水两相流平均流速测量结果，(b)为在含水率为80％工况下垂直井油水两相流平均流速测量结果。

具体实施方式

下面结合附图对本发明做进一步的说明：

具体实施方式一：这种基于改进的DPIV垂直井油水两相流流速场测量方法，包括以下步骤：

步骤一：选取时间间隔为Δt的两帧油水两相流图像，图像大小表示为：长×宽，设图像大小为M pixel×N pixel，M为图像长的值，N为图像宽的值，对其进行图像去噪、图像对比度增强后，确定初始查询窗口大小，初始查询窗口大小表示为：长×宽，设初始查询窗口大小为Wpixel×Wpixel，W为查询窗口长、宽的值，查询步长为W/2pixel，将两帧油水两相流图像按照相应初始查询窗口大小划分为覆盖率为50％的(2M-W)/W×(2N-W)/W个网格，设定当前查询窗口，令当前查询窗口大小等于初始查询窗口大小；

步骤二：在两帧图像中进行区域选择，选定区域表示为：[X坐标范围起点：X坐标范围终点，Y坐标范围起点：Y坐标范围终点]，设在第一帧图像中选定查询区域Area₁＝[i:i+W-1,j:j+W-1]，在第二帧图像中选定查询区域Area₂＝[i:i+W-1,j:j+W-1]，i，j分别表示图像中的X，Y坐标值，i＝1+W(n-1)，j＝1+W(m-1)，n＝1,2,…,(2M-W)/W，m＝1,2,…,(2N-W)/W，m为X方向查询区域序号，n为Y方向查询区域序号；对查询区域Area₁和查询区域Area₂两个区域进行ICP配准，得到查询区域Area₁平均X方向位移u_(i',j')，平均Y方向位移v_(i',j')，i'＝i+W/2为查询区域Area₁中心X坐标，j'＝j+W/2为查询区域Area₁中心Y坐标；

步骤三：当前查询窗口以步长W/2pixel遍历两帧图像，得到初始X方向位移场U_initial、初始Y方向位移场V_initial，设当前X方向位移场为U，当前Y方向位移场为V；

步骤四：根据当前X方向位移场U、当前Y方向位移场V对当前查询窗口变形，窗口变形是改变查询区域，重新进行ICP配准得到二次迭代X方向位移场U_new、二次迭代Y方向位移场V_new；

步骤五：更新当前X方向位移场U、当前Y方向位移场V，位移场按下式进行更新：

U＝U_initial+U_new，V＝V_initial+V_new；

步骤六：采用MLS对当前X方向位移场U、当前Y方向位移场V进行拟合，得到X方向边缘补充位移场U′、Y方向边缘补充位移场V′；X方向边缘补充位移场U′、Y方向边缘补充位移场V'按下述方法计算：

X方向边缘补充位移场U′曲面拟合函数f_u(x,y)、Y方向边缘补充位移场V'曲面拟合函数f_v(x,y)，x为X方向坐标变量，y为Y方向坐标变量，表示为k为多项式序号：

其中曲面拟合函数f_u(x,y)的系数阵列α_u(x,y)＝[α_u1(x,y),α_u2(x,y),…,α_uk(x,y)]，α_uk(x,y)为曲面拟合函数f_u(x,y)的第k个系数，曲面拟合函数f_v(x,y)的系数阵列α_v(x,y)＝[α_v1(x,y),α_v2(x,y),…,α_vk(x,y)]，α_vk(x,y)为曲面拟合函数f_v(x,y)的第k个系数，变量阵列μ(x,y)＝[μ₁(x,y),μ₂(x,y),...,μ_k(x,y)]＝[1,x,y,x²,xy,y²]，μ_k(x,y)为两个曲面拟合函数第k个变量，T代表矩阵转置符号；

α_u(x,y)、α_v(x,y)按下式求：

式中已知X方向位移阵列Z_u＝[u_{(W/2+1,W/2+1),}u_{(3W/2+1,3W/2+1)},…,u_(i',j')]，i'＝1+W(n-1/2)，j'＝1+W(m-1/2)，n＝1,2,…,(2M-W)/W，m＝1,2,…,(2N-W)/W，已知Y方向位移阵列Z_v＝[v_{(W/2+1,W/2+1),}v_{(3W/2+1,3W/2+1)},…,v_(i',j')]，i'＝1+W(n-1/2)，j'＝1+W(m-1/2)，n＝1,2,…,(2M-W)/W，m＝1,2,…,(2N-W)/W，参数阵列G＝[μ^T(_W/2+1,W/2+1),μ^T(3W/2+1,3W/2+1),…,μ^T(i',j')]，i'＝1+W(n-1/2)，j'＝1+W(m-1/2)，n＝1,2,…,(2M-W)/W，m＝1,2,…,(2N-W)/W，权值对角矩阵为具有紧支集特性的权函数；

将α_u(x,y)、α_v(x,y)代入f_u(x,y)、f_v(x,y)中，得到X方向边缘补充位移场U'和Y方向边缘补充位移场V'拟合曲面方程：

此时X方向边缘补充位移场U'和Y方向边缘补充位移场V'表示为：

此时当前X方向位移场U＝U′，当前Y方向位移场V＝V'；

步骤七：采用双三次均匀B样条插值扩大当前X方向位移场U、当前Y方向位移场V为原来的4倍；

步骤八：将当前查询窗口尺寸缩小为原来的1/4，得到新尺寸查询窗口，令当前查询窗口大小等于新尺寸查询窗口大小，根据当前X方向位移场U、当前Y方向位移场V对当前查询窗口变形，重新进行ICP配准得到三次迭代X方向位移场U′_new、三次迭代Y方向V′_new；

步骤九：更新当前X方向位移场U、当前Y方向位移场V；

步骤十：对当前查询窗口迭代进行步骤四-步骤九，直至缩小到指定查询窗口的大小，确定最终X方向位移场U_final、最终Y方向位移场V_final，根据时间间隔Δt和位移场U_final、V_final得到流速场f；油水两相流速度场按下式计算：

具体实施方式二：本实施方式为对实施方式一的进一步说明，所述步骤二中，W＝128，初始查询窗口尺寸128pixel×128pixel，查询区域分别为Area₁＝[i:i+127,j:j+127]、Area₂＝[i:i+127,j:j+127]，设查询区域Area₁的数据点集为p＝{p_ii|ii＝1,2,….,16384}，ii为查询区域Area₁的数据点序号，p_ii为数据点集p中第ii个数据点，查询区域Area₂的数据点集为q＝{q_jj|jj＝1,2,….,16384}，jj为查询区域Area₂的数据点序号，q_jj为数据点集q中第jj个数据点，该区域平均X方向位移u_(i',j')，平均Y方向位移v_(i',j')按如下方法计算：

数据点集变换关系式q'为：

q′＝rp+t

匹配目标函数E为：

其中r为旋转矩阵，t为平移向量。

采用SVD求解旋转矩阵r和平移向量t，对点集p、q进行如下变换，P_ii为p_ii经变换后的数据，Q_jj为q_jj经变换后的数据：

有：

对最优解矩阵H进行奇异值分解，将H分解为左奇异矩阵D，右奇异矩阵L，和奇异值矩阵Λ：

H＝DΛL^T

r和t的计算公式如下：

该区域的平均位移为：

u_(i′,j′)＝t(1)，v_(i′,j′)＝t(2)。

具体实施方式三：本实施方式为对实施方式二的进一步说明，所述步骤三中，初始X方向位移场U_initial，初始Y方向位移场V_initial表示为:

其中i'＝1+128(n-1/2)，j'＝1+128(m-1/2)，n＝1,2,…,(2M-128)/128，m＝1,2,…,(2N-128)/128。

当前X方向位移场U＝U_initial，当前Y方向位移场V＝V_initial；

具体实施方式四：本实施方式为对实施方式三的进一步说明，所述步骤四中经过当前查询窗口变形后，第一帧图像的查询区域为Area₁＝[i:i+W-1,j:j+W-1]，第二帧图像的查询区域为Area₂＝[i+u(i',j'):i+W-1+u(i',j'),j+v(i',j'):j+W-1+v(i',j')]。

具体实施方式五：本实施方式为对实施方式六的进一步说明，所述步骤七中双三次均匀B样条插值按下式计算：

式中Θ_u为当前X方向位移场U的双三次均匀B样条曲面片，Θ_v为当前Y方向位移场V的双三次均匀B样条曲面片，X方向坐标变量矩阵λ＝[ξ³ξ²ξ1]，ξ为X坐标变量，Y方向坐标变量矩阵γ＝[ζ³ζ²ζ1]，ζ为Y坐标变量，B_w为B样条基函数，C_u为当前X方向位移场U的双三次均匀B样条曲面片的控制顶点，C_v为当前Y方向位移场V的双三次均匀B样条曲面片的控制顶点，当X坐标变量ξ、Y坐标变量ζ在[0,1]之间遍历时，得到曲面片上任意一点的位移值。

具体实施方式六：本实施方式为对实施方式七的进一步说明，所述步骤八中，经当前查询窗口变形后第一帧图像的查询区域Area₁＝[i:i+W-1,j:j+W-1]，而第二帧图像的查询区域Area₂＝[i+u(i',j'):i+W-1+u(i',j'),j+v(i',j'):j+W-1+v(i',j')]。

具体实施方式七：本实施方式为对实施方式八的进一步说明，所述步骤九中位移场按下式更新：

U＝U′+U′_new，V＝V′+V′_new

具体实施方式八：下面结合图1、图2、图3、图4说明本实施方式：

一、以油滴作为示踪粒子，经LED背光光源照亮待测区域，采用高速摄像机连续拍摄多幅垂直井油水两相流图片；

二、选取时间间隔为Δt的两帧油水两相流图像，图像大小表示为：长×宽，如图1所示。设图像大小为M pixel×N pixel，M为图像长的值，N为图像宽的值，对其进行图像去噪、图像对比度增强后，确定初始查询窗口大小，查询窗口大小表示为：长×宽，设查询窗口大小为Wpixel×Wpixel，W为查询窗口长、宽的值，查询步长为W/2pixel，将两帧油水两相流图像按照相应查询窗口大小划分为覆盖率为50％的(2M-W)/W×(2N-W)/W个网格，设定当前查询窗口，令当前查询窗口大小等于初始查询窗口大小；

三、在两帧图像中进行区域选择，选定的区域表示为：[X坐标范围起点：X坐标范围终点，Y坐标范围起点：Y坐标范围终点]，设在第一帧图像中选定区域Area₁＝[i:i+127,j:j+127]，在第二帧图像中选定区域Area₂＝[i:i+127,j:j+127]，i，j分别表示图像中的X，Y坐标值，i＝1+128(n-1)，j＝1+128(m-1)，n＝1,2,…,(2M-128)/128，m＝1,2,…,(2N-128)/128，m为X方向查询区域序号，n为Y方向查询区域序号。对两个区域进行ICP配准，得到该区域平均X方向位移u_(i',j')，平均Y方向位移v_(i',j')，i'＝i+W/2为该区域中心X坐标，j'＝j+W/2为该区域中心Y坐标；

四、当前查询窗口以步长W/2pixel遍历两帧图像，移动方向如图2所示，得到初始X方向位移场U_initial、初始Y方向位移场V_initial，设当前X方向位移场为U，当前Y方向位移场为V。

五、根据当前X方向位移场U、当前Y方向位移场V对进行当前查询窗口变形，重新进行ICP配准得到二次迭代X方向位移场U_new、二次迭代Y方向位移场V_new。

六、更新当前X方向位移场U、当前Y方向位移场V。

七、采用MLS对X方向位移场U、Y方向位移场V进行拟合，得到X方向边缘补充位移场U'、Y方向边缘补充位移场V'。

八、采用双三次均匀B样条插值扩大当前位移场为原来的4倍。

九、将当前查询窗口大小缩小为原来的1/4(X方向、Y方向坐标范围分别缩小为原来的1/2)，得到新尺寸的查询窗口，令当前查询窗口大小等于新尺寸查询窗口大小，根据当前X方向位移场U、Y方向位移场V对进行当前查询窗口变形，重新进行ICP配准得到三次迭代X方向位移场U'_new、三次迭代Y方向V'_new。

十、更新当前X方向位移场U、当前Y方向位移场V。

十一、对当前查询窗口迭代进行五-十(对当前查询窗口重复步骤五-步骤十，又得到不同尺寸的查询窗口，再次重复步骤五-步骤十，如此重复)，直至查询窗口缩小到指定大小32pixel×32pixel，确定最终X方向位移场U_final、最终Y方向位移场V_final，根据时间间隔Δt和位移场U_final、V_final得到流速场f。

本发明解决了基于窗口迭代变形(Window deformation iterative multigrid,WIDIM)的DPIV算法在垂直井油水两相流流速场测量中由于图像匹配效果差、位移场边界缺失导致的测量精度降低问题。采用迭代最近点(Iterative Closest Point,ICP)代替传统互相关进行图像匹配，同时利用移动最小二乘(Moving least squares,MLS)综合整个流速场信息进行边界位移值补充，改善了DPIV的图像匹配效果，提高了垂直井油水两相流的测量精度。