CN109630092B - Data-based multi-model soft measurement method for pumping well pump efficiency - Google Patents
Data-based multi-model soft measurement method for pumping well pump efficiency Download PDFInfo
- Publication number
- CN109630092B CN109630092B CN201811350625.5A CN201811350625A CN109630092B CN 109630092 B CN109630092 B CN 109630092B CN 201811350625 A CN201811350625 A CN 201811350625A CN 109630092 B CN109630092 B CN 109630092B
- Authority
- CN
- China
- Prior art keywords
- data
- value
- ite
- follows
- output
- 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
Links
Images
Classifications
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH DRILLING; MINING
- E21B—EARTH DRILLING, e.g. DEEP DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B47/00—Survey of boreholes or wells
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH DRILLING; MINING
- E21B—EARTH DRILLING, e.g. DEEP DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B47/00—Survey of boreholes or wells
- E21B47/008—Monitoring of down-hole pump systems, e.g. for the detection of "pumped-off" conditions
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH DRILLING; MINING
- E21B—EARTH DRILLING, e.g. DEEP DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B47/00—Survey of boreholes or wells
- E21B47/06—Measuring temperature or pressure
Abstract
A multi-model soft measurement method for pumping well pump efficiency based on data comprises the following steps: building a pump efficiency soft measuring system, and respectively acquiring data of an upper stroke average load, a lower stroke average load, wellhead pressure, stroke and stroke frequency of the pumping unit, wellhead liquid outlet quantity, an upper stroke average current of the motor and a lower stroke average current of the motor; training and establishing a pumping well pump efficiency multi-model soft measurement model according to historical production data; and substituting the newly acquired data into the trained multi-model to obtain pump efficiency data, thereby realizing the soft measurement of the pump efficiency of the pumping well. The invention can realize the on-line continuous measurement of production data, has low measurement cost and strong real-time performance, can avoid the problems of data loss, recording errors and the like caused by manual recording, can avoid the problem that the decision process lags behind the actual working condition of the pumping well caused by sampling and testing the oil pumped out of the wellhead, has high prediction accuracy and can effectively process the problem of the influence of noise contained in the collected data sample on the prediction accuracy.
Description
Technical Field
The invention belongs to the technical field of oil extraction of oil wells, and particularly relates to a data-based multi-model soft measurement method for pumping well pump efficiency.
Background
The oil well pump is an important production part of the pumping well in the underground, and the pump efficiency of the oil well pump is an important parameter for evaluating the production capacity of the pumping well. The pump efficiency can reflect the working efficiency of the oil pumping equipment, the reasonable degree of pumping parameters and the management level. In the production process of the pumping well, the factors influencing the pump efficiency mainly comprise: the elastic expansion of the sucker rod and the oil pipe column, the gas influence, the insufficient liquid supply, the loose valve or the fixed valve, and the like. For the measurement of pumping efficiency of an oil pumping well, an oil liquid sampling analysis method, a mechanism modeling calculation method and the like exist at present. The oil sampling analysis method is to sample and test the oil pumped out of the wellhead, and calculate through a corresponding analysis method, but has the problem that the decision process lags behind the actual working condition of the pumping well; the mechanism modeling calculation method establishes a mechanism model for pump efficiency calculation according to a production process, but has the problems of excessive underground parameters which cannot be directly measured, lag in sampling time, strong uncertainty, larger calculation error and the like.
Disclosure of Invention
The invention aims to provide a data-based pumping well pump efficiency multi-model soft measurement method, which realizes remote acquisition of pumping well operation data and automatic calculation of pump efficiency.
The technical solution of the invention is as follows: a multi-model soft measurement method for pumping well pumping efficiency based on data comprises the following steps:
step 1: the method comprises the following steps of building a pump efficiency soft measurement system which consists of an indicator diagram wireless acquisition module, a wellhead pressure wireless acquisition module, a wellhead liquid output wireless acquisition module, a pumping unit motor power parameter wireless acquisition module, wireless data remote transmission equipment and a pumping well pump efficiency calculation server, and respectively acquiring upper stroke average load, lower stroke average load, wellhead pressure, pumping unit stroke and stroke frequency, wellhead liquid output, motor upper stroke average current and motor lower stroke average current data;
step 2: training and establishing multi-model soft measurement model for pumping well pump efficiency according to historical production data
Step 2.1: preprocessing collected historical production data of the pumping well, filling missing data, and taking an average value of a previous value and a next value of the missing data points as a substitute value;
step 2.2: taking the average load of the upper stroke, the average load of the lower stroke, the pressure of a well head, the stroke and the stroke frequency of a pumping unit, the liquid outflow quantity of oil flowing out of the well head, the average current of the upper stroke of a motor and the average current of the lower stroke of the motor in historical data as model input variables, taking the pump efficiency in the historical data as a model output variable, training and establishing a soft measurement model, and firstly performing dimensionality reduction on the data by a Principal Component Analysis (PCA) method, wherein the method comprises the following steps:
step 2.2.1: suppose the dataset is D = { (x) 1 ,y 1 ),(x 2 ,y 2 ),…,(x i ,y i ),…,(x m ,y m ) Where m is the number of samples in the data set, x i Input variables (n-dimensional), y, for the ith set of data i An output variable (1 dimension) for the ith group of data; the normalization process is performed on the m groups of data in the data set D as follows:
wherein i =1,2, …, m; j =1,2, …, n;
the average value of j dimension of m groups of data; x j ＝[x 1j ,x 2j ,…,x mj ] T A vector composed of the jth dimension of the m groups of data,
step 2.2.2: the covariance matrix of the normalized dataset D is found as follows:
of these, cov (X) j ,X j ) The calculation formula of (c) is as follows:
step 2.2.3: the eigenvalues of the covariance matrix and the corresponding eigenvectors are calculated from the following equation,
λ·q＝covD·q (4)
wherein, lambda and q are respectively a characteristic value and a characteristic vector;
step 2.2.4: the principal component cumulative contribution rate is calculated by the following formula,
wherein G is the accumulated contribution rate of the main component; p is the number of the characteristic values;
when the current cumulative contribution rate of the K principal components is greater than 80%, the K principal components can be considered to contain most information of the original data, and the feature vectors corresponding to the K principal components are used as feature vectors of an original data set, so that the data set is mapped to K dimensions from n dimensions;
step 2.2.5: a new data set is obtained from the K principal components as follows:
Z＝q T ·X (6)
wherein q = (q) 1 ,q 2 ,…,q K ) The first K feature vectors are obtained; x is the original input variable data set, X = [) 1 ,X 2 ,…,X n ](ii) a Z is a new input variable data set, Z = [ Z = [ [ Z ] 1 ,Z 2 ,…,Z K ]；
Step 2.2.6: mapping the original dataset D to a new dataset Z = { (Z) 1 ,y 1 ),(z 2 ,y 2 ),…,(z i ,y i ),…,(z m ,y m ) In which z is i Input variables (K dimension) for the ith group of data;
step 2.3: a density peak value clustering method based on a variant group fruit fly optimization algorithm for optimizing truncation distance is adopted to divide a new data set Z into a plurality of data subsets, and the method comprises the following steps:
step 2.3.1: initializing parameters and cutting off the distance d c As parameters to be optimized, the initial fruit fly population size Sizepop and the maximum iteration number Maxgen are given according to d c Value range of [ d ] cmin ,d cmax ]Wherein: d cmin And d cmax Randomly generating Sizepop fruit fly population positions [ X _ axis, Y _ axis ] as the minimum value and the maximum value of the truncation distance respectively]Calculating scale N, population change weight omega and step change weight xi by using the information entropy;
step 2.3.2: each fruit fly moves in random directions by smell, defined as follows:
wherein, [ X _ new _ axis [ ] i ,Y_new_axis i ]New location for the ith individual drosophila, i =1,2, …, sizepop; randomValue is the search distance;
step 2.3.3: calculating the distance Dist between the position of each individual fruit fly and the origin i Then, the taste concentration determination value S of the new position is calculated i Will S i As the truncated distance value to be optimized, the following is defined:
step 2.3.4: judging the taste concentration value S i Substituting into fitness function to obtain taste concentration Smell of each individual fruit fly position i The definition is as follows:
Smell i ＝f(S i )(9)
where f (-) is a fitness function defined as follows:
wherein, the first and the second end of the pipe are connected with each other,
ρ g the local density of data points g in the sample data set Z is represented as follows:
wherein d is gg’ Representing the Euclidean distance between any other g' point and the g point in the data set;
δ g the distance between all data points with local density larger than g point in the sample data set Z and the data point with the minimum distance from the g point is defined as follows:
step 2.3.5: the drosophila with the best taste concentration (designated as the best individual) in the drosophila population was found as follows:
[bestSmell,bestindex]＝min(Smell i ) (13)
step 2.3.6: the best taste concentration value bestsmlll is recorded and retained while the location of the best individual is recorded as follows:
Smellbest＝bestSmell (14)
step 2.3.7: after iterating Num _ ite (Num _ ite < Maxgen) times, the value of the entropy is calculated from the optimal taste concentration values of the last Num _ ite iteration as follows:
wherein EH (·) represents information entropy; cur _ ite represents the current cur _ ite iteration;
step 2.3.8: comparing the information entropy value obtained by the last Num _ ite iteration with the information entropy value obtained by the last Num _ ite iteration, and if the entropy value becomes larger, updating the population number as follows:
Sizepop(cur_ite+1)＝Sizepop(cur_ite)-w*Sizepop(cur_ite) (18)
wherein, sizepop (cur _ ite) and Sizepop (cur _ ite + 1) respectively represent cur _ ite and cur _ ite +1 iterations;
if the entropy becomes smaller, the population number and the new position of the individual drosophila are updated as follows:
Sizepop(cur_ite+1)＝Sizepop(cur_ite)+w*Sizepop(cur_ite) (19)
when the entropy value is constant for γ iterations, the update population number and the new location of the drosophila individual are as follows:
Sizepop(cur_ite+1)＝Sizepop(cur_ite)-w*Sizepop(cur_ite) (21)
step 2.3.9: repeating the steps 2.3.2-2.3.8 until the maximum iteration number is reached;
step 2.3.10: outputting the optimal truncation distance d c Obtaining the local density ρ of each data point according to the formula 11-formula 12 i And relative distance delta i Selecting rho i Value sum delta i The data points with the same label in the first 5% of the values arranged from large to small are taken as the clustering center c i (i＝1,2,...,r)；
Step 2.3.11: for the points except the clustering center in the data set, respectively calculating the Euclidean distance between the points and each clustering center, and dividing each point into the cluster class where the clustering center closest to the point is located, then, the data set Z is divided into r data subsets, which are expressed as: c 1 ,C 2 ,…,C r ；
Step 2.4: establishing a sub-model for each sample subset, which comprises the following steps:
step 2.4.1: for each data subset C i (i =1,2.., r), any given one input x j (j＝1,2,…,m i ) The relationship between input and output is generated by:
y ij ＝F i (x j )+ε i (23)
wherein, F i (x j ) A function representing the ith subset of data; epsilon i Is a mean of 0 and a variance of σ 2 Gaussian noise of (2);
step 2.4.2: each data subset C i (i =1,2.., r), for a new input
The corresponding probability prediction output is
Assuming that their relationship also satisfies the gaussian distribution, the output mean function and covariance function of the new input can be obtained as:
wherein, the first and the second end of the pipe are connected with each other,
a covariance matrix which is defined symmetrically positively; x is the number of j And y ij I =1,2, …, r, j =1,2, …, m, respectively, of the training set i ；k i Training a covariance matrix between data for the ith data subset;
in order to predict the input variables of a point,
in order to predict the covariance of the point itself,
covariance of the predicted point and the training data input variables;
and
the average of the predicted output and the predicted output, respectively, from which the training set output y can be derived ij And the predicted value output
The joint prior distribution between is:
wherein x is i Input variables for all training samples of the ith data subset; k (x) i ,x i )＝k(x iα ,x iβ )(α,β＝1,2,…,m i ) Is m i *m i Order symmetric positive definite covariance matrix for measuring x iα And x iβ Correlation between, k (x) iα ,x iβ ) Using the radial basis function as the covariance function, the following is defined:
wherein upsilon is 0 A metric representing local correlation; upsilon is 1 A variance representing a gaussian-distributed-compliant noise; theta iκ A dimension parameter representing an ith data subset; psi αβ Is Kronecker operator, if α = β, ψ αβ =1, otherwise ψ αβ ＝0；
Hyper-parameters of covariance function
The maximum likelihood estimation is adopted to obtain the maximum likelihood estimation, and the log likelihood function of the hyper-parameter is as follows:
wherein det (-) represents the determinant of the matrix;
training an output vector of data for the ith data subset;
can obtain new input
The posterior distribution of the output of (a) is:
then, can be made of
And
to obtain
Output value of
Step 2.4.3: and (3) performing Gaussian process regression modeling on the clustered r data subsets respectively according to formulas (23) to (30) to obtain r submodels, and recording as: model 1 ，Model 2 ，…，Model r ；
Step 2.5: the method comprises the following steps of fusing output values of a plurality of submodels to obtain a final output value, wherein the steps are as follows:
step 2.5.1: the weight for each sub-model is calculated as follows:
all weights satisfy the following relationship:
step 2.5.2: and carrying out weighted average on the output values of all the submodels to obtain a final output value as follows:
wherein y is the final output value; y is i Is the output value of the ith sub-model;
and step 3: and substituting the newly acquired upper stroke average load, lower stroke average load, wellhead pressure, stroke and stroke frequency of the pumping unit, the liquid outlet quantity of oil flowing out of the wellhead, the upper stroke average current of the motor and the lower stroke average current of the motor into the trained multi-model to obtain pump efficiency data, and realizing the soft measurement of the pump efficiency of the pumping well.
Further, the wellhead pressure includes oil pressure and casing pressure.
The invention has the beneficial effects that:
(1) The established pumping well pump efficiency soft measurement system adopts a wireless remote data acquisition mode to realize the acquisition of various production data of a pumping well, can realize the online continuous measurement of the production data, has low measurement cost and strong real-time performance, does not need technicians to operate relevant instruments on site to acquire the data, can reduce the labor cost, reduce the risk of site construction, can continuously and automatically store the data, and can avoid the problems of data loss, record errors and the like caused by manual recording.
(2) The adopted data-based soft measurement method can avoid the problem that the decision process caused by sampling and testing the oil pumped out of the wellhead lags behind the actual working condition of the pumping well, and can also effectively solve the problems of excessive underground non-directly measurable parameters, delayed sampling time, strong uncertainty, larger calculation error and the like in the process of establishing a pump efficiency calculation mechanism model by a production process. The continuous on-line measurement of the pumping well pump efficiency can be realized.
(3) The pumping well pump efficiency multi-model soft measurement method and system based on data, which are established by the invention, have the advantages of simple principle and small calculation complexity. The established multi-model soft measurement structure has high prediction accuracy, and can effectively solve the problem of influence of noise contained in the acquired data sample on the prediction accuracy.
Drawings
FIG. 1 is a schematic diagram of a pumping well remote wireless data acquisition system for pumping well pump efficiency calculation according to an embodiment of the present invention;
FIG. 2 is a 831 group historical on-trip average load data plot of an embodiment of the present invention;
FIG. 3 is a 831 set of historical downstroke average load data plots for an embodiment of the present invention;
FIG. 4 is a graph of 831 historical pumping unit stroke data for an embodiment of the present invention;
FIG. 5 is a graph of 831 sets of historical pumping unit stroke frequency data for an embodiment of the present invention;
FIG. 6 is a graph of 831 group historical wellhead oil pressure data for an embodiment of the present invention;
FIG. 7 is a 831 set historical wellhead casing pressure data plot of an embodiment of the present invention;
FIG. 8 is a graph of oil throughput data for a 831 group historical wellhead outflow of an embodiment of the present invention;
FIG. 9 is a graph of historical motor upstroke average current data for 831 groups according to an embodiment of the present invention;
FIG. 10 is a graph of 831 sets of historical motor downstroke average current data for an embodiment of the present invention;
FIG. 11 is a graph of 831 sets of historical pump efficiency data for an embodiment of the present invention;
FIG. 12 is a graph showing the results of principal component analysis according to the embodiment of the present invention;
fig. 13 is a clustering result diagram according to the embodiment of the present invention.
Detailed Description
The following detailed description of embodiments of the invention refers to the accompanying drawings.
Step 1: the adopted soft pump efficiency measuring system comprises an indicator diagram wireless acquisition module, a wellhead pressure wireless acquisition module, a wellhead liquid output wireless acquisition module, a pumping unit motor power parameter wireless acquisition module, wireless data remote transmission equipment and a pumping well pump efficiency calculation server. The indicator diagram wireless acquisition module is responsible for acquiring the average load of the upper stroke, the average load of the lower stroke, the stroke of the pumping unit and stroke frequency data; the wellhead pressure wireless acquisition module is responsible for acquiring wellhead oil pressure and casing pressure data; the wireless wellhead liquid outlet quantity acquisition module is used for acquiring the liquid outlet quantity of oil flowing out of a wellhead; the wireless acquisition module of the electric power parameters of the motor of the oil pumping unit is responsible for acquiring the average current of the upper stroke and the average current of the lower stroke of the motor; the wireless data remote transmission equipment is responsible for remotely transmitting various data acquired on an oil field production site to a pumping well pump efficiency calculation server located in an operation area management center; the pumping well pump efficiency calculation server is responsible for receiving and storing various production data acquired on an oil field production field and carrying out online calculation on the pumping well pump efficiency.
As shown in fig. 1, a wireless indicator diagram acquisition module, a wireless wellhead pressure acquisition module, a wireless wellhead liquid output acquisition module and a wireless pumping unit motor power parameter acquisition module are respectively installed at a wellhead of a pumping well, and an upper stroke average load, a lower stroke average load, a pumping unit stroke frequency, a wellhead oil pressure, a wellhead casing pressure, a liquid output amount of oil flowing out of the wellhead, an upper stroke average current of a motor and a lower stroke average current of the motor are respectively acquired through a wireless network 1 (which can be Zigbee, bluetooth or Wifi); the wireless data remote transmission equipment 1 packages various collected data and transmits the data to the wireless data remote transmission equipment 2 positioned in a management center of an oil field operation area through a wireless network 2 (which can be a wireless data transmission radio station, a GRPS (general packet radio service) or a wireless AP (access point); the wireless data remote transmission equipment 2 unpacks the received field production data and sends the unpacked field production data to the pumping well pump efficiency calculation server; and the pumping well pump efficiency calculation server receives and stores various production data acquired on a production field and performs on-line calculation on the pumping well pump efficiency.
Step 2: and training and establishing a multi-model soft measurement model of the pumping well pump efficiency according to historical production data.
Step 2.1: and preprocessing the collected historical production data of the pumping well, filling missing data, and taking the average value of the previous value and the next value of the missing data point as a substitute value.
Step 2.2: taking the average load of the upper stroke, the average load of the lower stroke, the pressure (including oil pressure and casing pressure) of a well head, the stroke and stroke frequency of a pumping unit, the liquid outlet amount of oil flowing out of the well head, the average current of the upper stroke of a motor and the average current of the lower stroke of the motor in historical data as model input variables, taking the pump efficiency in the historical data as a model output variable, training and establishing a soft measurement model, and firstly performing dimensionality reduction processing on the data by a Principal Component Analysis (PCA) method, wherein the method comprises the following steps:
step 2.2.1: establishing a data set D = { (x) by using 831 groups of historical data 1 ,y 1 ),(x 2 ,y 2 ),…,(x i ,y i ),…,(x m ,y m ) Where m is the number of samples in the dataset, m =831; x is the number of i Input variables (9D), y for the ith set of data i Is the output variable (1 dimension) of the ith group of data. Wherein, each dimension variable of 831 group data is shown in fig. 2 to fig. 11.
The data of 831 groups in the data set D are normalized as follows:
wherein i =1,2, …,831; j =1,2, …,9;
the average value of 831 groups of data in the j dimension; x j ＝[x 1j ,x 2j ,…,x mj ] T A vector composed of the jth dimension of 831 sets of data,
step 2.2.2: the covariance matrix of the normalized dataset D is found as follows:
of these, cov (X) j ,X j ) The calculation formula of (a) is as follows:
step 2.2.3: the eigenvalues of the covariance matrix and the corresponding eigenvectors are calculated from the following equation,
λ·q＝covD·q (4)
wherein λ and q are eigenvalues and eigenvectors, respectively.
Step 2.2.4: the principal component cumulative contribution rate is calculated by the following formula,
wherein G is the accumulated contribution rate of the main component; p is the number of eigenvalues.
When the current cumulative contribution rate of the K principal components is greater than 80%, the K principal components can be considered to contain most information of the original data, and the feature vectors corresponding to the K principal components are used as feature vectors of the original data set, so that the data set is mapped to the K dimension from the 11 dimension.
Step 2.2.5: a new data set is obtained from the K principal components as follows:
Z＝q T ·X (6)
wherein q = (q) 1 ,q 2 ,…,q K ) The first K eigenvectors; x is the original input variable data set, X = [) 1 ,X 2 ,…,X 9 ](ii) a Z is a new input variable data set, Z = [ Z = [) 1 ,Z 2 ,…,Z K ]。
Step 2.2.6: mapping the original dataset D to a new dataset Z = { (Z) 1 ,y 1 ),(z 2 ,y 2 ),…,(z i ,y i ),…,(z 831 ,y 831 ) In which z is i Is the input variable (K dimension) of the ith group of data.
From fig. 12, the cumulative contribution rate of the first 5 principal components is already greater than 80%, so the first 5 principal components can be considered to contain most of the information of the original data, and then the data set is mapped from 9 dimensions to 5 dimensions.
Step 2.3: a density peak value clustering method based on a variant group fruit fly optimization algorithm for optimizing truncation distance is adopted to divide a new data set Z into a plurality of data subsets, and the method comprises the following steps:
step 2.3.1: initializing parameters and cutting off the distance d c As parameters to be optimized, the initial fruit fly population size Sizepop and the maximum iteration number Maxgen are given according to d c Value range of [ d ] cmin ,d cmax ](wherein: d) cmin And d cmax Minimum and maximum of the truncation distance, respectively)) randomly generating Sizepop population positions [ X _ axis, Y _ axis]The information entropy calculates scale N, population change weight omega and step change weight xi.
Step 2.3.2: each fruit fly moves in random directions by olfaction, defined as follows:
wherein, [ X _ new _ axis [ ] i ,Y_new_axis i ]New location for ith individual drosophila, i =1,2, …, sizepop; randomValue is the search distance.
Step 2.3.3: calculating the distance Dist between the position of each individual fruit fly and the origin i Then, the taste concentration determination value S of the new position is calculated i (will S) i As a truncated distance value to be optimized), the following is defined:
step 2.3.4: judging taste concentration value S i Substituting into fitness function to obtain taste concentration Smell of each individual fruit fly position i Defined as follows:
Smell i ＝f(S i ) (9)
where f (-) is a fitness function defined as follows:
wherein the content of the first and second substances,
ρ g the local density of data points g in the sample data set Z is represented as follows:
wherein d is gg’ Representing the euclidean distance between any other g' point in the data set and the g point.
δ g The distance between all the data points with the local density larger than the g point in the sample data set Z and the data point with the minimum distance from the g point is defined as follows:
step 2.3.5: the drosophila with the best taste concentration (designated as the best individual) in the drosophila population was found as follows:
[bestSmell,bestindex]＝min(Smell i ) (13)
step 2.3.6: the best taste concentration value bestsmlll is recorded and retained while the location of the best individual is recorded as follows:
Smellbest＝bestSmell (14)
step 2.3.7: after iterating Num _ ite (Num _ ite < Maxgen) times, the value of the entropy is calculated from the optimal taste concentration values of the last Num _ ite iteration as follows:
wherein EH (·) represents information entropy; cur _ ite represents the current cur _ ite iteration.
Step 2.3.8: comparing the information entropy value obtained by the last Num _ ite iteration with the information entropy value obtained by the last Num _ ite iteration, and if the entropy value becomes larger, updating the population number as follows:
Sizepop(cur_ite+1)＝Sizepop(cur_ite)-w*Sizepop(cur_ite) (18)
where Sizepop (cur _ ite) and Sizepop (cur _ ite + 1) represent the cur _ ite and cur _ ite +1 iterations, respectively.
If the entropy becomes smaller, the population number and the new location of the Drosophila individual are updated as follows:
Sizepop(cur_ite+1)＝Sizepop(cur_ite)+w*Sizepop(cur_ite) (19)
when the entropy value is constant for γ iterations, the update population number and the new location of the drosophila individual are as follows:
Sizepop(cur_ite+1)＝Sizepop(cur_ite)-w*Sizepop(cur_ite) (21)
step 2.3.9: repeat steps 2.3.2-2.3.8 until the maximum number of iterations is reached.
Step 2.3.10: outputting the optimal truncation distance d c Obtaining the local density rho of each data point according to the formula (11) to the formula (12) i And relative distance delta i Choosing rho i Value sum delta i The data points with the same label in the first 5% of the values arranged from large to small are taken as the clustering center c i (i＝1,2,...,r)。
Step 2.3.11: for the points except the clustering center in the data set, respectively calculating the Euclidean distance between the points and each clustering center, and dividing each point into the cluster class where the clustering center closest to the point is located, then, the data set Z is divided into r data subsets, which are expressed as: c 1 ,C 2 ,…,C r 。
From fig. 13, data set Z is divided into 4 data subsets: c 1 、C 2 、C 3 And C 4 And respectively contains 351, 142, 181 and 157 sets of data.
Step 2.4: establishing a sub-model for each sample subset, which comprises the following steps:
step 2.4.1: for each data subset C i (i =1,2,3,4), any given one input x j (j＝1,2,…,m i ) The relationship between input and output is generated by:
y ij ＝F i (x j )+ε i (23)
wherein, F i (x j ) A function representing the ith subset of data; epsilon i Is a mean of 0 and a variance of σ 2 Gaussian noise.
Step 2.4.2: each data subset C i (i =1,2,3,4), for a new input
The corresponding probability prediction output is
Assuming that their relationship also satisfies the gaussian distribution, the output mean function and covariance function of the new input can be obtained as:
wherein the content of the first and second substances,
a covariance matrix which is defined symmetrically positively; x is a radical of a fluorine atom j And y ij I =1,2, …, r, j =1,2, …, m, respectively, the input variables and output variables of the training set i ；k i Training a covariance matrix between data for the ith data subset;
in order to predict the input variables of a point,
to predictThe covariance of the point itself is determined,
covariance of the predicted point and the training data input variables;
and
the average of the predicted output and the predicted output, respectively, from which the training set output y can be derived ij And the predicted value output
The joint prior distribution between is:
wherein x is i Input variables for all training samples of the ith data subset; k (x) i ,x i )＝k(x iα ,x iβ )(α,β＝1,2,…,m i ) Is m i *m i An order-symmetric positive definite covariance matrix for measuring x iα And x iβ Correlation between, k (x) iα ,x iβ ) Using the radial basis function as the covariance function, the following is defined:
wherein upsilon is 0 A metric representing local correlation; upsilon is 1 A variance representing a gaussian-distributed-compliant noise; theta iκ A dimension parameter representing an ith data subset; psi αβ Is Kronecker operator, if α = β, ψ αβ =1, otherwise ψ αβ ＝0。
Hyper-parameters of covariance function
The maximum likelihood estimation is adopted to obtain the maximum likelihood estimation, and the log likelihood function of the hyper-parameter is as follows:
wherein det (-) represents the determinant of the matrix;
an output vector of data is trained for the ith data subset.
Can obtain new input
The posterior distribution of the output of (a) is:
then, can be made of
And
to obtain
Output value of
Step 2.4.3: and (3) performing Gaussian process regression modeling on the 4 clustered data subsets according to formulas (23) to (30) to obtain 4 sub-models, and recording as: model 1 、Model 2 、Model 3 And a Model 4 。
Step 2.5: and fusing the output values of the 4 submodels to obtain a final output value. The method comprises the following steps:
step 2.5.1: the weight for each sub-model is calculated as follows:
all weights satisfy the following relationship:
step 2.5.2: and carrying out weighted average on the output values of all the submodels to obtain a final output value as follows:
wherein y is the final output value; y is i Is the output value of the ith sub-model.
And 3, step 3: and (2) substituting the newly acquired upper stroke average load, lower stroke average load, well head pressure (including oil pressure and casing pressure), stroke and stroke frequency of the pumping unit, the liquid outlet quantity of oil flowing out of the well head, the upper stroke average current of the motor and the lower stroke average current of the motor into the multi-model trained offline according to the historical production data in the step 2 to obtain pump efficiency data, and realizing the soft measurement of the pump efficiency of the pumping well.
Production data of six pumping wells in certain oil field operation area in China are adopted for verification, and relevant data are shown in a table 1.
TABLE 1 production data for six-well pumping wells
| Oil well 1 | Oil well 2 | Oil well 3 | Oil well 4 | Oil well 5 | Oil well 6 | |
| Upstroke mean load/KN | 75.0 | 76.2 | 82.4 | 83.5 | 75.3 | 82.7 |
| Down stroke mean load/KN | 46.0 | 47.1 | 43.4 | 48.0 | 46.6 | 44.6 |
| Oil pressure/Mpa | 0.3 | 0.2 | 0.2 | 0.2 | 0.25 | 0.2 |
| Sleeve pressure/Mpa | 1.9 | 2.5 | 0.3 | 0.6 | 1.2 | 2.2 |
| Stroke/ |
3 | 3 | 3 | 3 | 3 | 3 |
| Number of strokes/ |
4 | 6 | 3 | 4 | 6 | 4 |
| Oil liquid outlet amount/t | 6.3 | 5.1 | 17 | 25 | 18 | 5.4 |
| Upstroke average current/ |
45 | 42 | 56 | 42 | 53 | 52 |
| Mean of down strokeCurrent/mA | 27 | 38 | 34 | 38 | 27 | 30 |
The production data in table 1 are respectively substituted into the multi-model which is trained in step 2 in an off-line mode according to historical production data, and the obtained soft measurement results of the oil-gas-oil ratio are shown in table 2.
TABLE 2 Soft measurement results of pump efficiency for six-port pumping well
The above description is only exemplary of the present invention and is not intended to limit the present invention, and various modifications and changes may be made by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.
Claims (1)
1. A multi-model soft measurement method for pumping well pump efficiency based on data is characterized by comprising the following steps:
step 1: the method comprises the following steps of building a pump efficiency soft measurement system consisting of an indicator diagram wireless acquisition module, a wellhead pressure wireless acquisition module, a wellhead liquid output wireless acquisition module, a pumping unit motor power parameter wireless acquisition module, wireless data remote transmission equipment and a pumping well pump efficiency calculation server, and respectively acquiring upper stroke average load, lower stroke average load, wellhead pressure, pumping unit stroke and stroke frequency, wellhead liquid output, motor upper stroke average current and motor lower stroke average current data, wherein the wellhead pressure comprises oil pressure and casing pressure;
step 2: training and establishing multi-model soft measurement model of pumping well pump efficiency according to historical production data
Step 2.1: preprocessing collected historical production data of the pumping well, filling missing data, and taking an average value of a previous value and a next value of the missing data points as a substitute value;
step 2.2: taking the average load of the upper stroke, the average load of the lower stroke, the pressure of a well head, the stroke and the stroke frequency of a pumping unit, the liquid outflow quantity of oil flowing out of the well head, the average current of the upper stroke of a motor and the average current of the lower stroke of the motor in historical data as model input variables, taking the pump efficiency in the historical data as a model output variable, training and establishing a soft measurement model, and firstly performing dimensionality reduction on the data by a Principal Component Analysis (PCA) method, wherein the method comprises the following steps:
step 2.2.1: suppose the dataset is D = { (x) 1 ,y 1 ),(x 2 ,y 2 ),…,(x i ,y i ),…,(x m ,y m ) Where m is the number of samples in the dataset, x i For input variables of the ith group of data, i.e. n-dimensional input variables, y i Is an output variable of the ith group of data, namely a 1-dimensional output variable; the normalization process is performed on the m groups of data in the data set D as follows:
wherein i =1,2, …, m; j =1,2, …, n;
the average value of j dimension of m groups of data; x j ＝[x 1j ,x 2j ,…,x mj ] T A vector composed of the jth dimension of the m groups of data,
step 2.2.2: the covariance matrix of the normalized dataset D is found as follows:
of these, cov (X) j ,X j ) The calculation formula of (a) is as follows:
step 2.2.3: the eigenvalues of the covariance matrix and the corresponding eigenvectors are calculated from the following equation,
λ·q＝covD·q (4)
wherein, lambda and q are respectively a characteristic value and a characteristic vector;
step 2.2.4: the principal component cumulative contribution rate is calculated by the following formula,
wherein G is the accumulated contribution rate of the main component; p is the number of the characteristic values;
when the current accumulated contribution rate of the K principal components is greater than 80%, the K principal components are considered to contain most information of original data, the feature vectors corresponding to the K principal components are used as feature vectors of an original data set, and then the data set is mapped to K dimensions from n dimensions;
step 2.2.5: a new data set is obtained from the K principal components as follows:
Z＝q T ·X (6)
wherein q = (q) 1 ,q 2 ,…,q K ) The first K eigenvectors; x is the original input variable data set, X = [) 1 ,X 2 ,…,X n ](ii) a Z is a new input variable data set, Z = [ Z = [) 1 ,Z 2 ,…,Z K ]；
Step 2.2.6: mapping the original dataset D to a new dataset Z = { (Z) 1 ,y 1 ),(z 2 ,y 2 ),…,(z i ,y i ),…,(z m ,y m ) In which z is i The input variables of the ith group of data are K-dimensional input variables;
step 2.3: a density peak value clustering method based on a variety group fruit fly optimization algorithm for optimizing truncation distance is adopted to divide a new data set Z into a plurality of data subsets, and the method comprises the following steps:
step 2.3.1: initializing parameters and cutting off the distance d c As parameters to be optimized, the initial fruit fly population size Sizepop and the maximum iteration number Maxgen are given according to d c Value range of [ d ] cmin ,d cmax ]Wherein: d cmin And d cmax Randomly generating Sizepop fruit fly population positions [ X _ axis, Y _ axis ] for the minimum value and the maximum value of the truncation distance respectively]Calculating scale N, population change weight omega and step change weight xi by using the information entropy;
step 2.3.2: each fruit fly moves in random directions by olfaction, defined as follows:
wherein, [ X _ new _ axis [ ] i ,Y_new_axis i ]New location for ith individual drosophila, i =1,2, …, sizepop; randomValue is the search distance;
step 2.3.3: calculating the distance Dist between the position of each individual drosophila and the origin i Then, the taste concentration determination value S of the new position is calculated i Will S i As the truncated distance value to be optimized, the following is defined:
step 2.3.4: judging the taste concentration value S i Substituting fitness functionTo determine the taste concentration Smell of each individual location of Drosophila i The definition is as follows:
Smell i ＝f(S i ) (9)
where f (-) is a fitness function defined as follows:
wherein the content of the first and second substances,
ρ g represents the local density of data points g in the sample data set Z, defined as follows:
wherein d is gg’ Representing the Euclidean distance between any other g' point and the g point in the data set;
δ g the distance between all data points with local density larger than g point in the sample data set Z and the data point with the minimum distance from the g point is defined as follows:
step 2.3.5: the drosophila with the best taste concentration in the drosophila population was found and scored as the best individual as follows:
[bestSmell,bestindex]＝min(Smell i ) (13)
step 2.3.6: the best taste concentration value bestsmll is recorded and retained, while the location of the best individual is recorded, as follows:
Smellbest＝bestSmell (14)
step 2.3.7: after iterating Num _ ite (Num _ ite < Maxgen) times, the value of the entropy is calculated from the optimal taste concentration values of the last Num _ ite iteration as follows:
wherein EH (·) represents information entropy; cur _ ite represents the current cur _ ite iteration;
step 2.3.8: comparing the information entropy value obtained by the last Num _ ite iteration with the information entropy value obtained by the last Num _ ite iteration, and if the entropy value becomes larger, updating the population number as follows:
Sizepop(cur_ite+1)＝Sizepop(cur_ite)-w*Sizepop(cur_ite) (18)
where Sizepop (cur _ ite) and Sizepop (cur _ ite + 1) represent the cur _ ite and cur _ ite +1 iterations, respectively;
if the entropy becomes smaller, the population number and the new position of the individual drosophila are updated as follows:
Sizepop(cur_ite+1)＝Sizepop(cur_ite)+w*Sizepop(cur_ite) (19)
when the entropy value is constant for γ iterations, the update population number and the new location of the drosophila individual are as follows:
Sizepop(cur_ite+1)＝Sizepop(cur_ite)-w*Sizepop(cur_ite) (21)
step 2.3.9: repeating the steps 2.3.2-2.3.8 until the maximum iteration number is reached;
step 2.3.10: outputting the optimal truncation distance d c Obtaining the local density ρ of each data point according to the formula 11-formula 12 i And relative distance delta i Selecting rho i Value sum delta i The data points with the same label in the first 5% of the values arranged from large to small are taken as the clustering center c i (i＝1,2,...,r)；
Step 2.3.11: for the points except the clustering center in the data set, respectively calculating the Euclidean distance between the points and each clustering center, and dividing each point into the cluster class where the clustering center closest to the point is located, then, the data set Z is divided into r data subsets, which are expressed as: c 1 ,C 2 ,…,C r ；
Step 2.4: establishing a sub-model for each sample subset, comprising the following steps:
step 2.4.1: for each data subset C i (i =1,2.., r), any given one input x j (j＝1,2,…,m i ) The relationship between input and output is generated by:
y ij ＝F i (x j )+ε i (23)
wherein, F i (x j ) A function representing the ith subset of data; epsilon i Is a mean of 0 and a variance of σ 2 Gaussian noise of (2);
step 2.4.2: each data subset C i (i =1,2.., r), for a new input
The corresponding probability prediction output is
Assuming that their relationship also satisfies the gaussian distribution, the output mean function and covariance function of the new input can be obtained as:
wherein the content of the first and second substances,
a covariance matrix which is defined symmetrically positively; x is the number of j And y ij I =1,2, …, r, j =1,2, …, m, respectively, of the training set i ；k i Training a covariance matrix between data for the ith data subset;
in order to predict the input variables of a point,
in order to predict the covariance of the point itself,
covariance of the predicted point and the training data input variables;
and
the average of the predicted output and the predicted output, respectively, from which the training set output y can be derived ij And the predicted value output
The joint prior distribution between is:
wherein x is i Input variables for all training samples of the ith data subset; k (x) i ,x i )＝k(x iα ,x iβ )(α,β＝1,2,…,m i ) Is m i *m i An order-symmetric positive definite covariance matrix for measuring x iα And x iβ Correlation between, k (x) iα ,x iβ ) Using the radial basis function as the covariance function, the following is defined:
wherein upsilon is 0 A metric representing local correlation; upsilon is 1 A variance representing a gaussian-distributed-compliant noise; theta iκ A dimension parameter representing an ith data subset; psi αβ Is Kronecker operator, if α = β, ψ αβ =1, otherwise ψ αβ ＝0；
Hyper-parameters of covariance function
The maximum likelihood estimation is adopted to obtain the maximum likelihood estimation, and the log likelihood function of the hyper-parameter is as follows:
wherein det (-) represents the determinant of the matrix;
is the ith dataAn output vector of the subset training data;
can obtain new input
The posterior distribution of the output of (a) is:
then, can be made of
And
to obtain
Output value of
Step 2.4.3: and (3) performing Gaussian process regression modeling on the clustered r data subsets respectively according to formulas (23) to (30) to obtain r submodels, and recording as: model 1 ，Model 2 ，…，Model r ；
Step 2.5: the method comprises the following steps of fusing output values of a plurality of submodels to obtain a final output value, wherein the steps are as follows:
step 2.5.1: the weight for each sub-model is calculated as follows:
all weights satisfy the following relationship:
step 2.5.2: and carrying out weighted average on the output values of all the submodels to obtain a final output value as follows:
wherein y is the final output value; y is i Is the output value of the ith sub-model;
and step 3: and substituting the newly acquired upper stroke average load, lower stroke average load, wellhead pressure, stroke and stroke frequency of the pumping unit, the liquid outlet quantity of oil flowing out of the wellhead, the upper stroke average current of the motor and the lower stroke average current of the motor into the trained multi-model to obtain pump efficiency data, and realizing the soft measurement of the pump efficiency of the pumping well.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811350625.5A CN109630092B (en) | 2018-11-14 | 2018-11-14 | Data-based multi-model soft measurement method for pumping well pump efficiency |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811350625.5A CN109630092B (en) | 2018-11-14 | 2018-11-14 | Data-based multi-model soft measurement method for pumping well pump efficiency |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109630092A CN109630092A (en) | 2019-04-16 |
| CN109630092B true CN109630092B (en) | 2023-02-10 |
Family
ID=66067947
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811350625.5A Active CN109630092B (en) | 2018-11-14 | 2018-11-14 | Data-based multi-model soft measurement method for pumping well pump efficiency |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109630092B (en) |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103061749B (en) * | 2013-01-25 | 2015-07-15 | 东北大学 | Soft measuring method of dynamic liquid level under sucker rod oil well |
| CN105631554B (en) * | 2016-02-22 | 2019-11-26 | 渤海大学 | A kind of oil well oil liquid moisture content multi-model prediction technique based on time series |
| CN107909202B (en) * | 2017-11-03 | 2021-12-21 | 渤海大学 | Time sequence-based oil well liquid production integrated prediction method |
| CN107944607B (en) * | 2017-11-03 | 2022-01-18 | 渤海大学 | Time sequence-based pumping well shut-down time integrated prediction method |
| CN108805215B (en) * | 2018-06-19 | 2021-06-11 | 东北大学 | Dynamic liquid level soft measurement method for sucker-rod pump pumping well based on improved drosophila algorithm |
-
2018
- 2018-11-14 CN CN201811350625.5A patent/CN109630092B/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN109630092A (en) | 2019-04-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112101480B (en) | Multivariate clustering and fused time sequence combined prediction method | |
| CN107305601B (en) | Method for analyzing efficiency factors of oil pumping well system | |
| CN106677763B (en) | Dynamic integrated modeling-based oil well working fluid level prediction method | |
| CN110472689B (en) | Sucker-rod pump pumping well moving liquid level soft measurement method based on integrated Gaussian process regression | |
| CN106980910B (en) | Medium-and-long-term power load measuring and calculating system and method | |
| CN115270965A (en) | Power distribution network line fault prediction method and device | |
| CN107037306A (en) | Transformer fault dynamic early-warning method based on HMM | |
| CN107463993A (en) | Medium-and Long-Term Runoff Forecasting method based on mutual information core principle component analysis Elman networks | |
| CN104481496A (en) | Fault diagnosis method of sucker-rod pump well | |
| CN111401599A (en) | Water level prediction method based on similarity search and L STM neural network | |
| CN111985610A (en) | System and method for predicting pumping efficiency of oil pumping well based on time sequence data | |
| CN112364560B (en) | Intelligent prediction method for working hours of mine rock drilling equipment | |
| CN111695288B (en) | Transformer fault diagnosis method based on Apriori-BP algorithm | |
| Mishra et al. | Rainfall-runoff modeling using clustering and regression analysis for the river brahmaputra basin | |
| CN113780420A (en) | Method for predicting concentration of dissolved gas in transformer oil based on GRU-GCN | |
| CN112529053A (en) | Short-term prediction method and system for time sequence data in server | |
| CN110260914B (en) | Engineering safety monitoring system region division method based on measuring point space-time characteristics | |
| CN115127192A (en) | Semi-supervised water chilling unit fault diagnosis method and system based on graph neural network | |
| CN109630092B (en) | Data-based multi-model soft measurement method for pumping well pump efficiency | |
| CN114429238A (en) | Wind turbine generator fault early warning method based on space-time feature extraction | |
| CN113205125A (en) | XGboost-based extra-high voltage converter valve operation state evaluation method | |
| CN116303626A (en) | Well cementation pump pressure prediction method based on feature optimization and online learning | |
| CN109611077A (en) | A kind of pumpingh well oil liquid gas-oil ratio online soft sensor method based on multi-model | |
| CN112329813A (en) | Energy consumption prediction feature extraction method and system | |
| CN109000924B (en) | Method for monitoring state of ball screw pair based on K mean value |
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 |