CN110887976B - Improved DPIV vertical well-based oil-water two-phase flow velocity field measurement method - Google Patents

Improved DPIV vertical well-based oil-water two-phase flow velocity field measurement method Download PDF

Info

Publication number
CN110887976B
CN110887976B CN201911253732.0A CN201911253732A CN110887976B CN 110887976 B CN110887976 B CN 110887976B CN 201911253732 A CN201911253732 A CN 201911253732A CN 110887976 B CN110887976 B CN 110887976B
Authority
CN
China
Prior art keywords
displacement field
current
direction displacement
query
size
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201911253732.0A
Other languages
Chinese (zh)
Other versions
CN110887976A (en
Inventor
韩连福
丛垚
付长凤
刘超
刘兴斌
姜继玉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Dongying Ruigang Pipeline Engineering Co ltd
Original Assignee
Northeast Petroleum University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northeast Petroleum University filed Critical Northeast Petroleum University
Priority to CN201911253732.0A priority Critical patent/CN110887976B/en
Publication of CN110887976A publication Critical patent/CN110887976A/en
Application granted granted Critical
Publication of CN110887976B publication Critical patent/CN110887976B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01P—MEASURING LINEAR OR ANGULAR SPEED, ACCELERATION, DECELERATION, OR SHOCK; INDICATING PRESENCE, ABSENCE, OR DIRECTION, OF MOVEMENT
    • G01P5/00—Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft
    • G01P5/001—Full-field flow measurement, e.g. determining flow velocity and direction in a whole region at the same time, flow visualisation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01P—MEASURING LINEAR OR ANGULAR SPEED, ACCELERATION, DECELERATION, OR SHOCK; INDICATING PRESENCE, ABSENCE, OR DIRECTION, OF MOVEMENT
    • G01P5/00—Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft
    • G01P5/18—Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft by measuring the time taken to traverse a fixed distance
    • G01P5/22—Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft by measuring the time taken to traverse a fixed distance using auto-correlation or cross-correlation detection means

Abstract

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.

Description

Improved DPIV vertical well-based oil-water two-phase flow velocity field measurement method
The technical field is as follows:
the invention relates to the technical field of petroleum engineering and measurement, in particular to an improved DPIV vertical well oil-water two-phase flow velocity field measurement method.
Background art:
the flow velocity measurement is an essential part for dynamic monitoring of oil field development, and the oil-water two-phase flow velocity measurement methods include an ultrasonic method, a capacitance method, an electric conductivity method, an electromagnetic method, a heat tracing method, an optical fiber method and the like. Oilfield flow rate measurements of the current time period are more favored over zonal flow rate profile measurements because multiphase flow characteristics are more easily observed from zonal flow rate profiles, thereby guiding oilfield development. The traditional flow rate measurement method cannot meet the requirement, and a new flow rate measurement method is adopted. DPIV is an effective full-field speed measurement method in the field of hydrodynamics, and the characteristics of no disturbance, no contact and full-field speed measurement can well meet the requirements of multiphase flow speed measurement of oil fields, so that the DPIV has considerable application prospect in oil field speed measurement. The influence factors of the DPIV technology measurement accuracy are generally divided into two categories: one is the effect of external measurement environments such as trace particle distribution, illumination distribution, and various hardware systems, and the other is the implementation of DPIV image post-processing techniques. The Window Iterative Deformation (WIDIM) algorithm is used for improving the followability of the query Window under the condition that large-speed gradient tracer particles exist in a flow field, and is a DPIV algorithm widely used at present.
The oil-water two-phase flow field measurement by applying the WIDIM-based DPIV algorithm in the oil-water two-phase flow measurement of the 125mm vertical well has the following problems: (1) due to the fact that the diameter of the oil well is large, relative motion between oil drops in the query window is large in repeated reading, and the follow-up performance of the query window is high; (2) because the DPIV technology is applied to two-phase flow and adopts LED illumination and oil drop tracing, the gray level image of the two-phase flow to be processed is more complex than the gray level image information of the single-phase flow, so that the cross-correlation gray level calculation is easy to generate error matching to bring measurement errors, the errors can generate an amplification effect in the iteration process to seriously influence the image deformation effect, and the precision of the iteration calculation is improved; (3) the problem of displacement field boundary loss exists in the interpolation process of the WIDIM algorithm, and the window deformation effect is also influenced. Iterative Closest Point (ICP) is a Point cloud data registration method, which can realize accurate registration of a large amount of data, and applies a three-dimensional ICP technique to two-dimensional DPIV image registration, only a DPIV gray image needs to be projected to a three-dimensional space, and compared with cross-correlation matching, ICP registration is less affected by noise, and planar rotational motion of oil droplets is considered in the registration process. And assigning values to the missing displacement field boundary by comprehensively considering the whole displacement field information, so that a DPIV algorithm is further improved by adopting a Moving Least Square (MLS) method, MLS (Moving least square) surface fitting is carried out on the current displacement field before each interpolation, and the missing boundary value is supplemented according to a surface fitting function.
The invention content is as follows:
the invention aims to provide an improved DPIV vertical well oil-water two-phase flow velocity field measuring method which is used for solving the problems of poor image matching effect and inaccurate flow velocity measurement caused by the loss of a displacement field boundary in the vertical well oil-water two-phase flow measurement of a DPIV algorithm based on WIDIM.
The technical scheme adopted by the invention for solving the technical problems is as follows: the improved DPIV vertical well-based oil-water two-phase flow velocity field measurement method comprises the following steps:
the method comprises the following steps: selecting two frames of oil-water two-phase flow images with the time interval delta t, wherein the image size is expressed as: length is multiplied by width, the image size is set to be M pixel multiplied by N pixel, M is the value of the image length, N is the value of the image width, after image denoising and image contrast enhancement are carried out on the image width, the size of an initial query window is determined, and the size of the initial query window is expressed as: setting the length multiplied by the width, setting the size of an initial query window as W pixel multiplied by W pixel, setting W as the value of the length and the width of the query window, setting the query step length as W/2pixel, dividing two frames of oil-water two-phase flow images into (2M-W)/Wx (2N-W)/W grids with the coverage rate of 50% according to the size of the corresponding initial query window, setting the current query window, and making the size of the current query window equal to the size of the initial query window;
step two: region selection is performed in two images, and the selected region is expressed as: [ X coordinate Range Start: end of X-coordinate range, start of Y-coordinate range: end of range of Y coordinates]Setting the selected query Area in the first frame image1＝[i:i+W-1,j:j+W-1]Selecting a query Area in the second frame image2＝[i:i+W-1,j:j+W-1]I and j respectively represent X and Y coordinate values in the image, i is 1+ W (N-1), j is 1+ W (M-1), N is 1,2, …, (2M-W)/W, M is 1,2, …, (2N-W)/W, M is an X-direction query region number, and N is a Y-direction query region number; for the Area of query1And query Area2Performing ICP registration on the two areas to obtain an Area to be inquired1Average X-direction displacement u(i',j')Average Y-direction displacement v(i',j')I ═ i + W/2 is the query Area1The central X coordinate, j' ═ j + W/2 is the Area of inquiry1A center Y coordinate;
step three: the current query window is in step size W/2pTraversing two frames of images by ixel to obtain an initial X-direction displacement field UinitialInitial Y-direction displacement field VinitialSetting a current X-direction displacement field as U and a current Y-direction displacement field as V;
step four: deforming the current query window according to the current X-direction displacement field U and the current Y-direction displacement field V, and performing ICP registration again to obtain a secondary iteration X-direction displacement field UnewSecond iteration Y direction displacement field Vnew；
Step five: updating the current X-direction displacement field U and the current Y-direction displacement field V, and updating the displacement fields according to the following formula:
U＝Uinitial+Unew，V＝Vinitial+Vnew；
step six: fitting the current X-direction displacement field U and the current Y-direction displacement field V by adopting an MLS (Multi-level modeling System) to obtain an X-direction edge supplementary displacement field U 'and a Y-direction edge supplementary displacement field V'; the X-direction edge supplementary displacement field U 'and the Y-direction edge supplementary displacement field V' are calculated according to the following method:
x-direction edge supplementary displacement field U' curved surface fitting function fu(x, Y) and Y-direction edge complementary displacement field V' surface fitting function fv(X, Y), wherein X is an X-direction coordinate variable, Y is a Y-direction coordinate variable, and k is a polynomial serial number:
Figure BDA0002309732600000031
wherein the surface fitting function fuCoefficient array alpha of (x, y)u(x,y)＝[αu1(x,y),αu2(x,y),…,αuk(x,y)]，αuk(x, y) is a surface fitting function fu(x, y) kth coefficient, surface fitting function fvCoefficient array alpha of (x, y)v(x,y)＝[αv1(x,y),αv2(x,y),…,αvk(x,y)]，αvk(x, y) is a surface fitting function fv(x, y) k-th coefficient, and variable array μ (x, y) ═ μ1(x,y),μ2(x,y),...,μk(x,y)]＝[1,x,y,x2,xy,y2]，μk(x, y) is the kth variable of two surface fitting functions, and T represents a matrix transposition symbol;
αu(x,y)、αv(x, y) is calculated as:
Figure BDA0002309732600000041
wherein the known X-direction displacement array Zu＝[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, and it is known that the Y-direction displacement array Z is a linear displacement arrayv＝[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, and the parameter array G ═ μ MT(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, weight diagonal matrix
Figure BDA0002309732600000042
Is a weight function with tightly-supported characteristics;
will be alphau(x,y)、αvSubstitution of (x, y) into fu(x,y)、fvIn the step (X, Y), obtaining a fitting surface equation of an X-direction edge supplementary displacement field U 'and a Y-direction edge supplementary displacement field V':
Figure BDA0002309732600000043
the X-direction edge supplemental displacement field U 'and the Y-direction edge supplemental displacement field V' are now expressed as:
Figure BDA0002309732600000044
Figure BDA0002309732600000045
at this time, the current X-direction displacement field U is equal to U ', and the current Y-direction displacement field V is equal to V';
step seven: expanding the current X-direction displacement field U and the current Y-direction displacement field V by adopting bicubic uniform B spline interpolation to be 4 times of the original displacement field V;
step eight: reducing the size of the current query window to 1/4 of the original size to obtain a new-size query window, enabling the size of the current query window to be equal to the size of the new-size query window, deforming the current query window according to the current X-direction displacement field U and the current Y-direction displacement field V, and performing ICP (inductively coupled plasma) registration again to obtain a three-time iteration X-direction displacement field U'newAnd V 'in Y direction in three iterations'new；
Step nine: updating the current X-direction displacement field U and the current Y-direction displacement field V;
step ten: iterating the current query window processed in the step nine to perform the step four to the step nine until the current query window is reduced to the size of the specified query window, and determining the final X-direction displacement field UfinalFinal Y-direction displacement field VfinalAccording to the time interval Δ t and the displacement field Ufinal、VfinalObtaining a flow velocity field f; the oil-water two-phase flow velocity field is calculated according to the following formula:
Figure BDA0002309732600000051
the ICP registration method in step two of the above scheme: set W-128, initial query window size 128 pixels by 128 pixels, query Area1＝[i:i+127,j:j+127]Query Area2＝[i:i+127,j:j+127]Setting a query Area1Is p ═ pii1,2, …,16384, ii is the query Area1Data point number of piiQuerying the Area for the ii data point in the set of data points p2Is q ═ qjjJj ═ 1,2, …,16384}, where jj is the query Area2Data point number of (1), qjjFor the jj data in the data point set qPoint, average X-direction displacement u of the region(i',j')Average Y-direction displacement v(i',j')The calculation is as follows:
the data point set transformation relation q' is:
q′＝rp+t
the matching objective function E is:
Figure BDA0002309732600000052
where r is a rotation matrix, t is a translation vector, q'jjThe jj data point in the data point set after the data point set p is transformed by the data point set transformation relational expression q';
solving a rotation matrix r and a translational vector t by adopting SVD, and transforming point sets P and q as follows, PiiIs piiTransformed data, QjjIs qjjTransformed data:
Figure BDA0002309732600000061
comprises the following steps:
Figure BDA0002309732600000062
performing singular value decomposition on the optimal solution matrix H, and decomposing H into a left singular matrix D, a right singular matrix L and a singular value matrix Lambda, T representing a matrix transposition symbol:
H＝DΛLT
the calculation of r and t is as follows:
Figure BDA0002309732600000063
the average displacement of this region is:
u(i′,j′)＝t(1)，v(i′,j′)＝t(2)。
in the third step of the above-mentioned scheme,initial X-direction displacement field UinitialInitial Y-direction displacement field VinitialExpressed as:
Figure BDA0002309732600000071
wherein i ═ 1+128(N-1/2), j ═ 1+128(M-1/2), N ═ 1,2, …, (2M-128)/128, M ═ 1,2, …, (2N-128)/128; the current X-direction displacement field U is equal to UinitialThe current Y-direction displacement field V is equal to Vinitial。
After the current query window is deformed in the fourth step of the scheme, the query Area of the first frame image is Area1＝[i:i+W-1,j:j+W-1]The query Area of the second frame image is Area2＝[i+u(i',j'):i+W-1+u(i',j'),j+v(i',j'):j+W-1+v(i',j')]。
The bicubic uniform B-spline interpolation in the seventh step of the scheme is calculated according to the following formula:
Figure BDA0002309732600000072
Figure BDA0002309732600000081
in the formula (c) (-)uIs a bicubic uniform B-spline surface patch of the current X-direction displacement field UvIs a bicubic uniform B-spline surface sheet of the current Y-direction displacement field V, and the X-direction coordinate variable matrix lambda is ([ xi ])3ξ2ξ1]Xi is X coordinate variable, Y direction coordinate variable matrix gamma ═ zeta3ζ2ζ1]Zeta is a Y coordinate variable, BwAs a B-spline basis function, CuIs the control peak, C, of the bicubic uniform B-spline surface of the current X-direction displacement field UvIs the control vertex of a bicubic uniform B-spline curved surface sheet of a current Y-direction displacement field V, and when an X coordinate variable xi and a Y coordinate variable zeta are in [0,1 ]]And obtaining the displacement value of any point on the curved surface sheet during the process of the motion.
The current inquiry window in the step eight of the schemeAfter the mouth is deformed, the Area of the first frame image is inquired1＝[i:i+W-1,j:j+W-1]Query Area of second frame image2＝[i+u(i',j'):i+W-1+u(i',j'),j+v(i',j'):j+W-1+v(i',j')]。
In the above scheme step nine, the displacement field is updated according to the following formula:
U＝U′+U′new，V＝V′+V′new。
the invention has the following beneficial effects:
(1) according to the invention, the LED backlight light source is used for replacing the traditional laser sheet light source in the measuring process, and oil drops are used as tracer particles, so that the problem that the traditional tracer is shielded by the oil drops to influence the measurement in the oil-water two-phase flow is solved, and the measurement precision of the flow field and the flow rate is improved;
(2) the invention adopts ICP to replace cross-correlation algorithm to carry out image matching, simultaneously considers the translation and rotation of oil drops in a plane in the matching process, and can improve the matching precision of two-phase flow gray level images;
(3) the invention adopts MLS algorithm to synthesize the information of the whole flow field, supplements the missing displacement field boundary value, improves the follow-up property of the query window of the DPIV algorithm based on WIDIM, and improves the measurement precision.
Description of the drawings:
fig. 1 is a schematic diagram of image coordinates.
Fig. 2 is a schematic diagram of an image traversed by a query window, wherein the query window traverses the whole image from top to bottom from left to right.
FIG. 3 is a vertical well oil-water two-phase flow velocity field measured using a modified DPIV algorithm.
FIG. 4 shows the flow velocity of the two-phase flow of oil and water in the vertical well measured by the improved DPIV algorithm and compared with the classical DPIV measurement results, wherein (a) is the measurement result of the average flow velocity of the two-phase flow of oil and water in the vertical well under the condition that the water content is 90%, and (b) is the measurement result of the average flow velocity of the two-phase flow of oil and water in the vertical well under the condition that the water content is 80%.
Detailed Description
The invention is further described below with reference to the accompanying drawings:
the first embodiment is as follows: the improved DPIV vertical well-based oil-water two-phase flow velocity field measurement method comprises the following steps:
the method comprises the following steps: selecting two frames of oil-water two-phase flow images with the time interval delta t, wherein the image size is expressed as: length is multiplied by width, the image size is set to be M pixel multiplied by N pixel, M is the value of the image length, N is the value of the image width, after image denoising and image contrast enhancement are carried out on the image width, the size of an initial query window is determined, and the size of the initial query window is expressed as: setting the length multiplied by the width, setting the size of an initial query window as Wpixel multiplied by Wpixel, setting W as the value of the length and the width of the query window, setting the query step length as W/2pixel, dividing two frames of oil-water two-phase flow images into (2M-W)/Wx (2N-W)/W grids with the coverage rate of 50% according to the size of the corresponding initial query window, setting the current query window, and making the size of the current query window equal to the size of the initial query window;
step two: region selection is performed in two images, and the selected region is expressed as: [ X coordinate Range Start: end of X-coordinate range, start of Y-coordinate range: end of range of Y coordinates]Setting the selected query Area in the first frame image1＝[i:i+W-1,j:j+W-1]Selecting a query Area in the second frame image2＝[i:i+W-1,j:j+W-1]I and j respectively represent X and Y coordinate values in the image, i is 1+ W (N-1), j is 1+ W (M-1), N is 1,2, …, (2M-W)/W, M is 1,2, …, (2N-W)/W, M is an X-direction query region number, and N is a Y-direction query region number; for the Area of query1And query Area2Performing ICP registration on the two areas to obtain an Area to be inquired1Average X-direction displacement u(i',j')Average Y-direction displacement v(i',j')I ═ i + W/2 is the query Area1The central X coordinate, j' ═ j + W/2 is the Area of inquiry1A center Y coordinate;
step three: traversing two frames of images by the current query window by step length W/2pixel to obtain an initial X-direction displacement field UinitialInitial Y-direction displacement field VinitialSetting a current X-direction displacement field as U and a current Y-direction displacement field as V;
step four: according to the current X-direction displacement field U and the current Y-direction displacement field V, the current query window is deformed, and the window deformation is changedInquiring the region, and performing ICP registration again to obtain a secondary iteration X-direction displacement field UnewSecond iteration Y direction displacement field Vnew；
Step five: updating the current X-direction displacement field U and the current Y-direction displacement field V, and updating the displacement fields according to the following formula:
U＝Uinitial+Unew，V＝Vinitial+Vnew；
step six: fitting the current X-direction displacement field U and the current Y-direction displacement field V by adopting an MLS (Multi-level modeling System) to obtain an X-direction edge supplementary displacement field U 'and a Y-direction edge supplementary displacement field V'; the X-direction edge supplementary displacement field U 'and the Y-direction edge supplementary displacement field V' are calculated according to the following method:
x-direction edge supplementary displacement field U' curved surface fitting function fu(x, Y) and Y-direction edge complementary displacement field V' surface fitting function fv(X, Y), wherein X is an X-direction coordinate variable, Y is a Y-direction coordinate variable, and k is a polynomial serial number:
Figure BDA0002309732600000101
wherein the surface fitting function fuCoefficient array alpha of (x, y)u(x,y)＝[αu1(x,y),αu2(x,y),…,αuk(x,y)]，αuk(x, y) is a surface fitting function fu(x, y) kth coefficient, surface fitting function fvCoefficient array alpha of (x, y)v(x,y)＝[αv1(x,y),αv2(x,y),…,αvk(x,y)]，αvk(x, y) is a surface fitting function fv(x, y) k-th coefficient, and variable array μ (x, y) ═ μ1(x,y),μ2(x,y),...,μk(x,y)]＝[1,x,y,x2,xy,y2]，μk(x, y) is the kth variable of two surface fitting functions, and T represents a matrix transposition symbol;
αu(x,y)、αv(x, y) is calculated as:
Figure BDA0002309732600000111
wherein the known X-direction displacement array Zu＝[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, and it is known that the Y-direction displacement array Z is a linear displacement arrayv＝[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, and the parameter array G ═ μ MT(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, weight diagonal matrix
Figure BDA0002309732600000112
Is a weight function with tightly-supported characteristics;
will be alphau(x,y)、αvSubstitution of (x, y) into fu(x,y)、fvIn the step (X, Y), obtaining a fitting surface equation of an X-direction edge supplementary displacement field U 'and a Y-direction edge supplementary displacement field V':
Figure BDA0002309732600000113
the X-direction edge supplemental displacement field U 'and the Y-direction edge supplemental displacement field V' are now expressed as:
Figure BDA0002309732600000114
Figure BDA0002309732600000121
at this time, the current X-direction displacement field U is equal to U ', and the current Y-direction displacement field V is equal to V';
step seven: expanding the current X-direction displacement field U and the current Y-direction displacement field V by adopting bicubic uniform B spline interpolation to be 4 times of the original displacement field V;
step eight: reducing the size of the current query window to 1/4 of the original size to obtain a new-size query window, enabling the size of the current query window to be equal to that of the new-size query window, deforming the current query window according to the current X-direction displacement field U and the current Y-direction displacement field V, and performing ICP (inductively coupled plasma) registration again to obtain a three-time iteration X-direction displacement field U'newAnd V 'in Y direction in three iterations'new；
Step nine: updating the current X-direction displacement field U and the current Y-direction displacement field V;
step ten: iterating the current query window to carry out the steps four to nine until the current query window is reduced to the size of the specified query window, and determining the final X-direction displacement field UfinalFinal Y-direction displacement field VfinalAccording to the time interval Δ t and the displacement field Ufinal、VfinalObtaining a flow velocity field f; the oil-water two-phase flow velocity field is calculated according to the following formula:
Figure BDA0002309732600000122
the second embodiment is as follows: in the second step, W is 128, the initial query window size is 128 pixels × 128 pixels, and the query areas are Area respectively1＝[i:i+127,j:j+127]、Area2＝[i:i+127,j:j+127]Setting a query Area1Is p ═ pii1,2, …,16384, ii is the query Area1Data point number of piiQuerying the Area for the ii data point in the set of data points p2Is q ═ qjjJj ═ 1,2, …,16384}, where jj is the query Area2Data point number of (1), qjjFor the jj data point in the data point set q, the average X-direction displacement u of the region(i',j')Average Y-direction displacement v(i',j')The calculation is as follows:
the data point set transformation relation q' is:
q′＝rp+t
the matching objective function E is:
Figure BDA0002309732600000131
where r is the rotation matrix and t is the translation vector.
Solving a rotation matrix r and a translational vector t by adopting SVD, and transforming point sets P and q as follows, PiiIs piiTransformed data, QjjIs qjjTransformed data:
Figure BDA0002309732600000132
comprises the following steps:
Figure BDA0002309732600000133
performing singular value decomposition on the optimal solution matrix H, and decomposing H into a left singular matrix D, a right singular matrix L and a singular value matrix Lambda:
H＝DΛLT
the calculation of r and t is as follows:
Figure BDA0002309732600000134
the average displacement of this region is:
u(i′,j′)＝t(1)，v(i′,j′)＝t(2)。
the third concrete implementation mode: in the third step, the initial X-direction displacement field U is used to further explain the second embodimentinitialInitial Y-direction displacement field VinitialExpressed as:
Figure BDA0002309732600000141
where i ═ 1+128(N-1/2), j ═ 1+128(M-1/2), N ═ 1,2, …, (2M-128)/128, M ═ 1,2, …, (2N-128)/128.
The current X-direction displacement field U is equal to UinitialThe current Y-direction displacement field V is equal to Vinitial；
The fourth concrete implementation mode: in the fourth step, after the current query window is deformed, the query Area of the first frame image is Area1＝[i:i+W-1,j:j+W-1]The query Area of the second frame image is Area2＝[i+u(i',j'):i+W-1+u(i',j'),j+v(i',j'):j+W-1+v(i',j')]。
The fifth concrete implementation mode: in this embodiment, to further explain the sixth embodiment, the bicubic uniform B-spline interpolation in the seventh step is calculated according to the following formula:
Figure BDA0002309732600000142
Figure BDA0002309732600000151
in the formula (c) (-)uIs a bicubic uniform B-spline surface patch of the current X-direction displacement field UvIs a bicubic uniform B-spline surface sheet of the current Y-direction displacement field V, and the X-direction coordinate variable matrix lambda is ([ xi ])3ξ2ξ1]Xi is X coordinate variable, Y direction coordinate variable matrix gamma ═ zeta3ζ2ζ1]Zeta is a Y coordinate variable, BwAs a B-spline basis function, CuIs the control peak, C, of the bicubic uniform B-spline surface of the current X-direction displacement field UvIs the control vertex of a bicubic uniform B-spline curved surface sheet of a current Y-direction displacement field V, and when an X coordinate variable xi and a Y coordinate variable zeta are in [0,1 ]]And obtaining the displacement value of any point on the curved surface sheet during the process of the motion.
The sixth specific implementation mode: in the eighth step, the query Area of the first frame image after the deformation of the current query window is obtained1＝[i:i+W-1,j:j+W-1]And the Area of the query of the second frame image2＝[i+u(i',j'):i+W-1+u(i',j'),j+v(i',j'):j+W-1+v(i',j')]。
The seventh embodiment: in this embodiment, the eighth embodiment is further described, and the displacement field in the ninth step is updated according to the following formula:
U＝U′+U′new，V＝V′+V′new
the specific implementation mode is eight: the present embodiment will be described below with reference to fig. 1,2, 3, and 4:
firstly, oil drops are used as tracer particles, an LED backlight light source illuminates an area to be measured, and a high-speed camera is adopted to continuously shoot a plurality of vertical well oil-water two-phase flow pictures;
selecting two frames of oil-water two-phase flow images with the time interval delta t, wherein the image size is represented as: length x width as shown in fig. 1. Setting the image size as M pixel multiplied by N pixel, wherein M is the value of the image length, N is the value of the image width, and after image denoising and image contrast enhancement are carried out on the image size, determining the size of an initial query window, wherein the size of the query window is represented as: the method comprises the steps of length-width setting, wherein the size of a query window is Wpixel-Wpixel, W is the value of the length and the width of the query window, the query step length is W/2pixel, dividing two frames of oil-water two-phase flow images into (2M-W)/Wx (2N-W)/W grids with the coverage rate of 50% according to the size of the corresponding query window, setting a current query window, and enabling the size of the current query window to be equal to the size of an initial query window;
and thirdly, selecting areas in the two images, wherein the selected areas are expressed as follows: [ X coordinate Range Start: end of X-coordinate range, start of Y-coordinate range: end of range of Y coordinates]Setting in a selected Area in the first frame image1＝[i:i+127,j:j+127]Selecting Area in the second frame image2＝[i:i+127,j:j+127]And i and j respectively represent X and Y coordinate values in the image, i is 1+128(N-1), j is 1+128(M-1), N is 1,2, …, (2M-128)/128, M is 1,2, …, (2N-128)/128, M is an X-direction query region number, and N is a Y-direction query region number. Carrying out ICP registration on the two regions to obtain the average X-direction displacement u of the region(i',j')Average Y-direction displacement v(i',j')I ═ i + W/2 is the center X of the zoneThe mark j ═ j + W/2 is the Y coordinate of the center of the region;
fourthly, the current query window traverses two frames of images by step length W/2pixel, the moving direction is as shown in figure 2, and an initial X-direction displacement field U is obtainedinitialInitial Y-direction displacement field VinitialLet the current X-direction displacement field be U and the current Y-direction displacement field be V.
Fifthly, according to the current X-direction displacement field U and the current Y-direction displacement field V, the current query window is deformed, and ICP registration is carried out again to obtain a secondary iteration X-direction displacement field UnewSecond iteration Y direction displacement field Vnew。
And sixthly, updating the current X-direction displacement field U and the current Y-direction displacement field V.
And seventhly, fitting the displacement field V in the direction of the X-direction displacement field U, Y by adopting an MLS (Multi-level quadratic System) to obtain an edge supplementary displacement field U 'in the X direction and an edge supplementary displacement field V' in the Y direction.
And eighthly, expanding the current displacement field to be 4 times of the original displacement field by adopting bicubic uniform B spline interpolation.
Ninthly, reducing the size of the current query window to 1/4 (the coordinate ranges in the X direction and the Y direction are respectively reduced to 1/2) to obtain a query window with a new size, enabling the size of the current query window to be equal to the size of the query window with the new size, deforming the current query window according to the current X-direction displacement field U, Y-direction displacement field V, and performing ICP registration again to obtain a three-time iteration X-direction displacement field U'newAnd V 'in Y direction in three iterations'new。
And tenthly, updating the current X-direction displacement field U and the current Y-direction displacement field V.
Eleven, iterating the current query window for five-ten times (repeating the steps five-ten for the current query window, obtaining query windows with different sizes again, repeating the steps five-ten again, and repeating the steps in the same way) until the query window is reduced to a specified size of 32 pixels multiplied by 32 pixels, and determining a final X-direction displacement field UfinalFinal Y-direction displacement field VfinalAccording to the time interval Δ t and the displacement field Ufinal、VfinalA flow velocity field f is obtained.
The invention solves the problem that the DPIV algorithm based on Window Iterative Deformation (WIDIM) reduces the measurement precision in the measurement of the oil-water two-phase flow velocity field of the vertical well due to poor image matching effect and displacement field boundary loss. And an Iterative Closest Point (ICP) is adopted to replace the traditional cross-correlation for image matching, and meanwhile, Moving Least Squares (MLS) is utilized to synthesize the information of the whole flow velocity field for boundary displacement value supplement, so that the image matching effect of the DPIV is improved, and the measurement precision of the oil-water two-phase flow of the vertical well is improved.

Claims (3)

1. An improved DPIV vertical well-based oil-water two-phase flow velocity field measurement method is characterized by comprising the following steps:
the method comprises the following steps: selecting two frames of oil-water two-phase flow images with the time interval delta t, wherein the image size is expressed as: length is multiplied by width, the image size is Mpixel multiplied by N pixel, M is the value of the image length, N is the value of the image width, after image denoising and image contrast enhancement are carried out on the image width, the size of an initial query window is determined, and the size of the initial query window is expressed as: length is multiplied by width, the size of an initial query window is set to be Wpixel multiplied by Wpixel, W is a value of the length and the width of the query window, the query step length is W/2pixel, two frames of oil-water two-phase flow images are divided into a plurality of grids with fixed coverage rate according to the size of the corresponding initial query window, a current query window is set, the size of the current query window is equal to the size of the initial query window, and the size of the window coverage rate and the number of the grids are as follows:
the window coverage rate is 50%, and the two frames of oil-water two-phase flow images are divided into (2M-W)/Wx (2N-W)/W grids according to the window coverage rate;
step two: region selection is performed in two images, and the selected region is expressed as: [ X coordinate Range Start: end of X-coordinate range, start of Y-coordinate range: end of range of Y coordinates]Setting the selected query Area in the first frame image1＝[i:i+W-1,j:j+W-1]Selecting a query Area in the second frame image2＝[i:i+W-1,j:j+W-1]I and j respectively represent X and Y coordinate values in the image, i is 1+ W (N-1), j is 1+ W (M-1), N is 1,2, …, (2M-W)/W, M is 1,2, …, (2N-W)/W, and M is an X-direction query regionThe sequence number, n is the sequence number of the query area in the Y direction;
for the Area of query1And query Area2Performing ICP registration on the two areas to obtain an Area to be inquired1Average X-direction displacement u (i ', j'), average Y-direction displacement v (i ', j'), i + W/2 is the query Area1The central X coordinate, j' ═ j + W/2 is the Area of inquiry1A center Y coordinate;
step three: traversing two frames of images by the current query window by step length W/2pixel to obtain an initial X-direction displacement field UinitialInitial Y-direction displacement field VinitialSetting a current X-direction displacement field as U and a current Y-direction displacement field as V;
step four: deforming the current query window according to the current X-direction displacement field U and the current Y-direction displacement field V, and performing ICP registration again to obtain a secondary iteration X-direction displacement field UnewSecond iteration Y direction displacement field Vnew；
Step five: updating the current X-direction displacement field U and the current Y-direction displacement field V, and updating the displacement fields according to the following formula:
U＝Uinitial+Unew，V＝Vinitial+Vnew；
step six: fitting the current X-direction displacement field U and the current Y-direction displacement field V by adopting an MLS (Multi-level modeling System) to obtain an X-direction edge supplementary displacement field U 'and a Y-direction edge supplementary displacement field V'; the X-direction edge supplementary displacement field U 'and the Y-direction edge supplementary displacement field V' are calculated according to the following method:
x-direction edge supplementary displacement field U' curved surface fitting function fu(x, Y) and Y-direction edge complementary displacement field V' surface fitting function fv(X, Y), X is X direction coordinate variable, Y is Y direction coordinate variable, k is polynomial serial number:
Figure FDA0003047480950000021
wherein the surface fitting function fuCoefficient array alpha of (x, y)u(x,y)＝[αu1(x,y),αu2(x,y),…,αuk(x,y)]，αuk(x, y) is a surface fitting function fu(x, y) kth coefficient, surface fitting function fvCoefficient array alpha of (x, y)v(x,y)＝[αv1(x,y),αv2(x,y),…,αvk(x,y)]，αvk(x, y) is a surface fitting function fv(x, y) k-th coefficient, and variable array μ (x, y) ═ μ1(x,y),μ2(x,y),...,μk(x,y)]＝[1,x,y,x2,xy,y2]，μk(x, y) is the kth variable of two surface fitting functions, and T represents a matrix transposition symbol;
αu(x,y)、αv(x, y) is calculated as:
Figure FDA0003047480950000031
wherein the known X-direction displacement array Zu＝[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, and it is known that the Y-direction displacement array Z is a linear displacement arrayv＝[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, and the parameter array G ═ μ MT(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, weight diagonal matrix
Figure FDA0003047480950000032
Figure FDA0003047480950000033
Is a weight function with tightly-supported characteristics;
will be alphau(x,y)、αvSubstitution of (x, y) into fu(x,y)、fvIn (X, Y), fitting is carried out to obtain an X-direction edge supplementary displacement field U' and a Y-direction edge supplementary displacement field VThe curved surface equation:
Figure FDA0003047480950000034
the X-direction edge supplemental displacement field U 'and the Y-direction edge supplemental displacement field V' are now expressed as:
Figure FDA0003047480950000041
Figure FDA0003047480950000042
at this time, the current X-direction displacement field U is equal to U ', and the current Y-direction displacement field V is equal to V';
step seven: expanding the current X-direction displacement field U by adopting bicubic uniform B spline interpolation, enabling the current Y-direction displacement field V to be 4 times of the original Y-direction displacement field V, reducing the size of a current query window to be 1/4 of the original size to obtain a new-size query window, enabling the size of the current query window to be equal to the size of the new-size query window, deforming the current query window according to the current X-direction displacement field U and the current Y-direction displacement field V, and performing ICP registration again to obtain a triple-iteration X-direction displacement field U'newAnd three times of iteration Y-direction displacement field V'new；
Step eight: updating the current X-direction displacement field U and the current Y-direction displacement field V;
step nine: iterating the current query window in the fourth step to the eighth step until the current query window is reduced to the size of the specified query window, and determining the final X-direction displacement field UfinalFinal Y-direction displacement field VfinalAccording to the time interval Δ t and the displacement field Ufinal、VfinalObtaining a flow velocity field f; the oil-water two-phase flow velocity field is calculated according to the following formula:
Figure FDA0003047480950000043
2. the improved DPIV vertical well oil-water two-phase flow velocity field measurement method as claimed in claim 1, wherein the method comprises the following steps: in the third step, the displacement field U in the X direction is initiatedinitialInitial Y-direction displacement field VinitialExpressed as:
Figure FDA0003047480950000051
the size of the current query window is 128 × 128; the current X-direction displacement field U is equal to UinitialThe current Y-direction displacement field V is equal to Vinitial。
3. The improved DPIV vertical well oil-water two-phase flow velocity field measurement method as claimed in claim 1, wherein the method comprises the following steps: after the current query window is deformed in the fourth step, the query Area of the first frame image is Area1＝[i:i+W-1,j:j+W-1]The query Area of the second frame image is Area2＝[i+u(i′,j′):i+W-1+u(i′,j),j+v(i′,j′):j+W-1+v(i′,j′)]。
CN201911253732.0A 2019-12-09 2019-12-09 Improved DPIV vertical well-based oil-water two-phase flow velocity field measurement method Active CN110887976B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911253732.0A CN110887976B (en) 2019-12-09 2019-12-09 Improved DPIV vertical well-based oil-water two-phase flow velocity field measurement method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911253732.0A CN110887976B (en) 2019-12-09 2019-12-09 Improved DPIV vertical well-based oil-water two-phase flow velocity field measurement method

Publications (2)

Publication Number Publication Date
CN110887976A CN110887976A (en) 2020-03-17
CN110887976B true CN110887976B (en) 2021-06-01

Family

ID=69751199

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911253732.0A Active CN110887976B (en) 2019-12-09 2019-12-09 Improved DPIV vertical well-based oil-water two-phase flow velocity field measurement method

Country Status (1)

Country Link
CN (1) CN110887976B (en)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5792962A (en) * 1994-07-05 1998-08-11 Institut Francais Du Petrole Device and method for measuring velocity profiles in a multiphase fluid
CN2606895Y (en) * 2003-03-19 2004-03-17 申功炘 Digital particle image velocity measurement system
CN104599268A (en) * 2015-01-06 2015-05-06 广州医科大学附属肿瘤医院 Local area accurate deformation registration algorithm combining point registration
CN108280850A (en) * 2018-03-05 2018-07-13 北京中科嘉宁科技有限公司 A kind of fast image registration method optimized based on B-spline and Levenberg-Marquardt

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5792962A (en) * 1994-07-05 1998-08-11 Institut Francais Du Petrole Device and method for measuring velocity profiles in a multiphase fluid
CN2606895Y (en) * 2003-03-19 2004-03-17 申功炘 Digital particle image velocity measurement system
CN104599268A (en) * 2015-01-06 2015-05-06 广州医科大学附属肿瘤医院 Local area accurate deformation registration algorithm combining point registration
CN108280850A (en) * 2018-03-05 2018-07-13 北京中科嘉宁科技有限公司 A kind of fast image registration method optimized based on B-spline and Levenberg-Marquardt

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Flow Velocity Field Measurement of Vertical Upward Oil–Water Two-Phase Immiscible Flow Using the Improved DPIV Algorithm Based on ICP and MLS;Lianfu Han 等;《APPLIED SCIENCES》;20190811;第9卷(第16期);第1-22页 *

Also Published As

Publication number Publication date
CN110887976A (en) 2020-03-17

Similar Documents

Publication Publication Date Title
Cai et al. Dense motion estimation of particle images via a convolutional neural network
CN108416840B (en) Three-dimensional scene dense reconstruction method based on monocular camera
CN107204009B (en) Three-dimensional point cloud registration method based on affine transformation model CPD algorithm
Anjyo et al. Scattered data interpolation for computer graphics
CN106228544A (en) A kind of significance detection method propagated based on rarefaction representation and label
CN112258658B (en) Augmented reality visualization method based on depth camera and application
CN104091350B (en) A kind of object tracking methods of utilization motion blur information
CN101504770B (en) Structural light strip center extraction method
CN110926339A (en) Real-time three-dimensional measurement method based on one-time projection structured light parallel stripe pattern
CN103826032A (en) Depth map post-processing method
CN113390605B (en) Full-field measurement method for wing deformation of wind tunnel test airplane
CN110688440B (en) Map fusion method suitable for less sub-map overlapping parts
CN110887976B (en) Improved DPIV vertical well-based oil-water two-phase flow velocity field measurement method
Yu et al. Deep dual recurrence optical flow learning for time-resolved particle image velocimetry
CN102982523A (en) Multisource and multi-focus color image fusion method
CN112446179A (en) Fluid velocity measuring method based on variable split optical flow model
CN103411562B (en) A kind of structured light strip center extraction method based on dynamic programming and average drifting
CN111724428A (en) Depth map sampling and reconstructing method based on-map signal model
CN110363713A (en) High spectrum image noise-reduction method based on recurrence sample scaling and bilinearity Factorization
CN106600629B (en) A kind of light stream estimation method towards Large Scale Motion
CN112508007B (en) Space target 6D attitude estimation method based on image segmentation Mask and neural rendering
CN111693729B (en) Particle image velocity measurement method and device based on global optimization
CN112927169B (en) Remote sensing image denoising method based on wavelet transformation and improved weighted kernel norm minimization
CN114298950A (en) Infrared and visible light image fusion method based on improved GoDec algorithm
CN111640084A (en) High-speed pixel matching method based on LK optical flow

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right

Effective date of registration: 20221103

Address after: No. 19, Gangcheng Road, Dongying Port Economic Development Zone, Dongying City, Shandong Province 257237

Patentee after: Dongying Ruigang Pipeline Engineering Co.,Ltd.

Address before: 163319 No. 99 Xuefu Street, Daqing Hi-tech Development Zone, Heilongjiang Province

Patentee before: NORTHEAST PETROLEUM University

TR01 Transfer of patent right