METHOD FOR PREDICTING AND OPTIMIZING RATE OF PENETRATION (ROP) FOR OIL AND GAS DRILLING BASED ON BAYESIAN OPTIMIZATION

Abstract

A method for predicting and optimizing an ROP for oil and gas drilling based on Bayesian optimization includes: acquiring raw drilling data according to a preset sampling period, constructing an initial sample data set based on the raw drilling data, constructing an ROP prediction model based on the initial sample data set, and predicting an ROP at a next sample point through a Gaussian process regression based on the ROP prediction model. The present invention realizes rapid analysis of historical drilling data and accurate prediction of an ROP range at a sample point in a feasible domain. The method can obtain an optimized ROP and an optimized engineering parameter through Bayesian optimization, and obtain an engineering parameter corresponding to an optimal ROP. The method has few restrictions on the drilling engineering parameter and the raw formation parameter, improves prediction accuracy, and avoids the problem of blurred parameter value boundaries.

Claims

1. A method for optimizing a rate of penetration of a well drilling operation using a predicted rate of penetration (ROP) for oil and gas drilling based on Bayesian optimization, comprising the following steps: S1: while drilling a well, acquiring raw drilling data according to a preset sampling period by using a sensor, wherein the raw drilling data comprises feature parameters of a drilling engineering parameter, a raw formation parameter, and a sequence parameter, the drilling engineering parameter is selected from the group consisting of weight of bit (WOB), rotation speed, torque, drilling fluid density, displacement, riser pressure, drilling tool, and combinations thereof, the raw formation parameter is selected from the group consisting of natural gamma ray (GR), acoustic travel time (AC), lithology, and combinations thereof, and the sequence parameter is the well depth; S2: constructing an initial sample data set based on the raw drilling data acquired in step S1; S3: constructing an ROP prediction model based on the initial sample data set constructed in step S2; S4: predicting an ROP at a next sample point through a Gaussian process regression based on the ROP prediction model constructed in step S3; S5: constructing an ROP optimization model based on the ROP prediction model and a prediction result in step S4; wherein step S5 comprises: constructing a sample point set based on the prediction result obtained by the ROP prediction model in step S4 and the initial sample data set constructed in step S1; and constructing a sampling function by selecting a maximum value of the ROP function of each sample point that is calculated by the ROP prediction model, and obtaining the ROP optimization model; and S6: obtaining an optimized ROP based on the ROP optimization model constructed in step S5; step S6 comprises: S61: obtaining an ROP fluctuation range of the next sample point through a Bayesian method based on the ROP of the next sample point obtained by the ROP prediction model and a mean and a variance of each feature parameter in the sampling period; S62: obtaining a corresponding ROP by using the raw geological formation parameter and raw well depth of each sample point in the sample point set in step S5 as fixed parameters and randomly changing the engineering parameter; S63: retaining an ROP within the ROP fluctuation range and a corresponding engineering parameter, and updating the retained ROP and the corresponding engineering parameter to the sample point set; S64: obtaining a value of the sampling function of a current sample point set; S65: determining whether a number of obtaining the value of the sampling function in step S64 reaches a preset threshold; and if yes, proceeding to step S66, otherwise returning to step S61 based on the current sample point set; and S66: taking an updated ROP corresponding to the maximum value of the sampling function as the optimized ROP, and taking an updated engineering parameter corresponding to the maximum value of the sampling function as an optimized engineering parameter, changing the engineering parameter of the well drilling operation to the updated engineering parameter so that well drilling continues with the updated engineering parameter to provide the optimized ROP.

2. The method according to claim 1, wherein step S2 comprises: constructing spatial sequence data being spatially continuous based on well depth acquired in step S1, and obtaining the initial sample data set by taking the engineering parameter and the raw formation parameter features and taking an ROP in the raw drilling data as a data label.

3. The method according to claim 1, wherein step S3 comprises: S31: calculating a mean vector of each sample point based on the initial sample data set constructed in step S2; S32: constructing a covariance matrix of each sample point through a Gaussian kernel function based on the initial sample data set constructed in step S2: $.Math.\left(x_{1:n},x_{1:n}\right)=\left(\begin{array}{ccc}\mathrm{cov}\left(x_1,x_1\right)& .Math.& \mathrm{cov}\left(x_1,x_n\right)\\ .Math.& & .Math.\\ \mathrm{cov}\left(x_n,x_1\right)& .Math.& \mathrm{cov}\left(x_n,x_n\right)\end{array}\right)=\left(\begin{array}{ccc}k\left(x_1,x_1\right)& .Math.& k\left(x_1,x_n\right)\\ .Math.& & .Math.\\ k\left(x_n,x_1\right)& .Math.& k\left(x_n,x_n\right)\end{array}\right)$ wherein, Σ(x.sub.1:n,x.sub.1:n) denotes the covariance matrix of each sample point in the initial sample data set; k(x.sub.n,x.sub.1) denotes a Gaussian kernel function of a sample point pair (x.sub.n,x.sub.1) in the initial sample data set; and cov(x.sub.n,x.sub.1) denotes a covariance of the sample point pair (x.sub.n,x.sub.1) in the initial sample data set; and S33: constructing an ROP function through a Gaussian function based on the mean vector calculated in step S31 and the covariance matrix constructed in step S32, and obtaining the ROP prediction model.

4. The method according to claim 3, wherein the Gaussian kernel function in step S32 is expressed as: $k\left(x_1,x_2\right)={\alpha }_0\exp \left(-\frac{1}{2{\sigma }^2}{.Math.x_1-x_2.Math.}^2\right)$ wherein, k(x.sub.1,x.sub.2) denotes a Gaussian kernel function of a sample point pair (x.sub.1,x.sub.2) in the initial sample data set; α.sub.0 and σ are kernel function parameters; and ∥.Math.∥.sup.2 are 2-norms.

5. The method according to claim 3, wherein the ROP function constructed in step S33 is expressed as:
ƒ(x.sub.1:n)˜gp(μ(x.sub.1:n),Σ(x.sub.1:n,x.sub.1:n)) wherein, ƒ(x.sub.1:n) denotes an ROP function in a continuous domain; μ(x.sub.1:n) denotes the mean vector of each sample point in the initial sample data set; Σ(x.sub.1:n,x.sub.1:n) denotes the covariance matrix of each sample point in the initial sample data set; gp(.Math.) denotes the Gaussian function; and ˜ denotes mapping.

Description

BRIEF DESCRIPTION OF DRAWINGS

[0059] FIG. 1 is a flowchart of a method for predicting an ROP for oil and gas drilling based on Bayesian optimization according to the present invention;

[0060] FIG. 2 is a flowchart of step S3 of the prediction method according to the present invention;

[0061] FIG. 3 is a flowchart of a method for optimizing an ROP for oil and gas drilling based on Bayesian optimization according to the present invention;

[0062] FIG. 4 shows fitting of a sampling function according to the present invention; and

[0063] FIG. 5 is a flowchart of step S6 of the prediction method according to the present invention.

DETAILED DESCRIPTION OF THE EMBODIMENTS

[0064] The specific implementations of the present invention are described below to facilitate those skilled in the art to understand the present invention, but it should be clear that the present invention is not limited to the scope of the specific implementations. Various obvious changes made by those of ordinary skill in the art within the spirit and scope of the present invention defined by the appended claims should fall within the protection scope of the present invention.

[0065] As shown in FIG. 1, an embodiment of the present invention provides a method for predicting an ROP for oil and gas drilling based on Bayesian optimization, which includes the following sub steps:

[0066] S1: Acquire raw drilling data according to a preset sampling period.

[0067] In a practical application, all oil wells in an oilfield area are sampled in a sampling period defined by an interval of five meters. An engineering parameter of a drilling tool is measured by using a sensor, and a formation parameter at a time of measurement is recorded. The drilling engineering parameter and the raw formation parameter had a total of K dimensions.

[0068] S2: Construct an initial sample data set based on the raw drilling data acquired in step S1.

[0069] In this embodiment, Step S2 includes:

[0070] Construct spatial sequence data being spatially continuous based on well depth data included in the raw drilling data acquired in step S1, and obtain the initial sample data set by taking an engineering parameter included in the raw drilling data and a raw formation parameter as sample features and taking an ROP included in the raw drilling data as a data label.

[0071] In a practical application, all oil wells in an oilfield area are sampled in a sampling period defined by an interval of five meters. An engineering parameter of a drilling tool is measured by using a sensor, and raw drilling data X={x.sub.n}.sub.1.sup.N is acquired, N being a dimension of a sample feature, that is, a number of sample features. The drilling engineering parameter and the raw formation parameter in the drilling data are taken as the sample features L={L.sub.k}.sub.1.sup.K, K being the dimension of the sample features in the drilling data and the raw formation data. The ROP is taken as a data label Y={y.sub.n}.sub.1.sup.N, and the parameters are divided into the drilling engineering parameter, the raw formation parameter and a sequence parameter. The drilling engineering parameter affecting the ROP includes WOB, rotation speed, torque, drilling fluid density, flow rate, etc. The raw formation parameter affecting the ROP includes natural gamma ray GR, acoustic travel time difference AC, etc. The sequence parameter includes well depth. All drilling data included in the drilling engineering parameter, the raw formation parameter and the sequence parameter are input. With the ROP label, an initial sample data set D is constructed, as shown in Table 1.

TABLE-US-00001 TABLE 1 Data sample parameters in the initial sample data set Type of input parameter Feature name Drilling engineering parameter WOB Rotation speed Torque Drilling fluid density Displacement Riser pressure Drilling tool . . . Raw formation parameter Natural gamma ray GR Acoustic travel time difference AC Lithology . . . Sequence parameter Well depth S3: Construct an ROP prediction model based on the initial sample data set constructed in step S2. As shown in FIG. 2, Step S3 includes: S31: Calculating a mean vector of each sample point based on the initial sample data set constructed in step S2.

[0072] In a practical application, a mean function μ(x) is used to calculate a mean vector of each sample point x in the initial sample data set, that is, μ(x.sub.1:n)=μ(x.sub.1) . . . μ(x.sub.n)).sup.T, where the mean vector of each sample point is expressed as: μ(x.sub.1:n), ( ).sup.T denoting a transposition.

[0073] S32: Construct a covariance matrix of each sample point through a Gaussian kernel function based on the initial sample data set constructed in step S2:

[00004]$.Math.\left(x_{1:n},x_{1:n}\right)=\left(\begin{array}{ccc}\mathrm{cov}\left(x_1,x_1\right)& .Math.& \mathrm{cov}\left(x_1,x_n\right)\\ .Math.& & .Math.\\ \mathrm{cov}\left(x_n,x_1\right)& .Math.& \mathrm{cov}\left(x_n,x_n\right)\end{array}\right)=\left(\begin{array}{ccc}k\left(x_1,x_1\right)& .Math.& k\left(x_1,x_n\right)\\ .Math.& & .Math.\\ k\left(x_n,x_1\right)& .Math.& k\left(x_n,x_n\right)\end{array}\right)$

[0074] where, Σ(x.sub.1:n,x.sub.1:n) denotes the covariance matrix of each sample point in the initial sample data set; k(x.sub.n,x.sub.1) denotes a Gaussian kernel function of a sample point pair (x.sub.n,x.sub.1) in the initial sample data set; and cov(x.sub.n,x.sub.1) denotes a covariance of the sample point pair (x.sub.n,x.sub.1) in the initial sample data set.

[0075] In this embodiment, the Gaussian kernel function in step S32 is expressed as:

[00005]$k\left(x_1,x_2\right)={\alpha }_0\exp \left(-\frac{1}{2{\sigma }^2}{.Math.x_1-x_2.Math.}^2\right)$

[0076] where, k(x.sub.1,x.sub.2) denotes a Gaussian kernel function of a sample point pair (x.sub.1,x.sub.2) in the initial sample data set; α.sub.0 and σ are kernel function parameters; and ∥.Math.∥.sup.2 are 2-norms.

[0077] In the present invention, the parameters α.sub.0 and σ of the kernel function both may take a value of 1.

[0078] S33: Construct an ROP function through a Gaussian function based on the mean vector calculated in step S31 and the covariance matrix constructed in step S32, and obtain the ROP prediction model.

[0079] In this embodiment, the ROP function constructed in step S33 is expressed as:

ƒ(x.sub.1:n)˜gp(μ(x.sub.1:n),Σ(x.sub.1:n,x.sub.1:n))

[0080] where, ƒ(x.sub.1:n) denotes an ROP function in a continuous domain, which defines the ROP prediction model; μ(x.sub.1:n) denotes the mean vector of each sample point in the initial sample data set, which is a vector composed of the mean of each feature parameter in each sampling period; Σ(x.sub.1:n,x.sub.1:n) denotes the covariance matrix of each sample point in the initial sample data set, which is a matrix composed of the covariance of each feature parameter in each sampling period; gp(.Math.) denotes the Gaussian function; and ˜ denotes mapping.

[0081] S4: Predict an ROP at a next sample point through a Gaussian process regression based on the ROP prediction model constructed in step S3.

[0082] In a practical application, the corresponding function value ƒ(x.sub.1:n) of the known sampling group x.sub.1:n in the initial sample data set is calculated by the ROP function in the ROP prediction model. According to the corresponding function value ƒ(x.sub.1:n) of the known sampling group x.sub.1:n in the initial sample data set, a Gaussian process regression is performed through the ROP prediction model to predict a new set of sample values x.sub.n+1 of the drilling parameter and corresponding ROP values. In other words, a posterior probability distribution p(ƒ(x.sub.n+1)|x.sub.n+1, D) is calculated. According to the properties of the Gaussian distribution, both the marginal distribution and the conditional distribution of the Gaussian distribution are Gaussian distributions. Therefore, the posterior probability distribution satisfies a multivariate Gaussian distribution, thereby obtaining the ROP at the next sample point.

[0083] As shown in FIG. 3, the present invention further provides a method for optimizing an ROP, which is based on the method for predicting an ROP for oil and gas drilling based on Bayesian optimization, and includes steps S5 and S6:

[0084] S5: Construct an ROP optimization model based on the ROP prediction model and a prediction result in step S4.

[0085] In this embodiment, step S5 includes:

[0086] Construct a sample point set based on the prediction result obtained by the ROP prediction model in step S4 and the initial sample data set constructed in step S1; and construct a sampling function by selecting a maximum value of the ROP function of each sample point that is calculated by the ROP prediction model, and obtain the ROP optimization model.

[0087] In a practical application, the initial sample data set D is taken as a sample point set D.sub.m. According to the data calculated by the ROP prediction model, that is, the mean and variance at each sample point in the feasible domain and the maximum value corresponding to the ROP function at each sample point, a sampling function is constructed:

[00006]$u_{\mathrm{EI}}\left(x\right)=\left(\mu \left(x\right)-f_n^{*}\right)-\left(\mu \left(x\right)-f_n^{*}\right)\Phi \left(\frac{\left(\mu \left(x\right)-f_n^{*}\right)}{\sigma \left(x\right)}\right)+\sigma \left(x\right)\phi \left(\frac{\left(\mu \left(x\right)-f_n^{*}\right)}{\sigma \left(x\right)}\right)$

[0088] where, u.sub.EI(x) denotes the sampling function; μ(x) denotes a mean of the sampling function at searched n sample points; σ(x) denotes a variance of the sampling function at the searched n sample points; ƒ.sub.n* denotes a maximum value of the ROP function at the searched n sample points, which is expressed as: ƒ.sub.n*=max(ƒ(x.sub.1), . . . , ƒ(x.sub.n)), φ(.Math.) denotes a probability density function of a standard normal distribution (SND); and Φ(.Math.) denotes a distribution function of the SND.

[0089] In a practical application, the sampling function is an expected improvement (EI) function, which considers the mean and variance of the function value at each point, and considers the optimal function value that has been found in previous iterations. Through the mean, variance and sampling function, the sampling times are significantly reduced, and the next sample point is screened, as shown in FIG. 4.

[0090] S6: Obtain an optimized ROP based on the ROP optimization model constructed in step S5.

[0091] As shown in FIG. 5, in this embodiment, Step S6 includes:

[0092] S61: Obtain an ROP fluctuation range of the next sample point through Bayesian method based on the ROP of the next sample point obtained by the ROP prediction model and a mean and a variance of each feature parameter in the sampling period.

[0093] S62: Obtain a corresponding ROP by using a raw geological parameter and a raw well depth of each sample point in the sample point set in step S5 as fixed parameters and randomly changing the engineering parameter.

[0094] In a practical application, during iteration, the raw geological parameter and the well depth are taken as fixed constants, and only the engineering parameter is taken as a variable. The sampling function is calculated in the fluctuation range, and the sampling function value at each sample point is taken as the function value of this point. Taking x=[3890, a.sub.2, a.sub.3, a.sub.4, a.sub.5, a.sub.6, a.sub.7, 1, a.sub.9, 48.39, 43.843] as an example, the well depth is 3890, the lithology category is 1, that is, sandstone, the acoustic travel time difference AC is 48.39, and the natural gamma GR is 43.843. Under the above fixed well depth and geological parameters, a.sub.2, a.sub.3, a.sub.4, a.sub.5, a.sub.6, a.sub.7, a.sub.9, denote engineering parameter variables, respectively. The sampling function u.sub.EI(x) is calculated by randomly changing the engineering parameter.

[0095] S63: Retain an ROP within the ROP fluctuation range and a corresponding engineering parameter, and update the retained ROP and the corresponding engineering parameter to the sample point set.

[0096] In a practical application, the maximum value of the sampling function obtained randomly within the fluctuation range is compared. The next sample point x.sub.n+1 is determined according to the maximum value of the sampling function, x.sub.n+1=argmax.sub.xu.sub.EI(x). argmax.sub.x(.Math.) is an extreme value function, which is used to obtain the value of the sample point corresponding to the maximum value of the sampling function, and obtain the next sample point x.sub.n+1. The retained ROP and the corresponding engineering parameter are updated to the sample point set, which is expressed as: D.sub.m={D.sub.m,x.sub.n+1}.

[0097] S64: Obtain a value of the sampling function of a current sample point set.

[0098] S65: Determine whether a number of obtaining the value of the sampling function in step S64 reaches a preset threshold; and if yes, proceed to step S66, otherwise return to step S61 based on the current sample point set.

[0099] S66: Take an updated ROP corresponding to the maximum value of the sampling function as the optimized ROP, and take an updated engineering parameter corresponding to the maximum value of the sampling function as an optimized engineering parameter.

[0100] In a practical application, the updated ROP of the sample point corresponding to the maximum value of the sampling function is obtained as the optimized ROP, so as to obtain the optimized engineering parameter combination under the fixed well depth and fixed raw geological parameter.

[0101] The present invention is described with reference to the flowcharts and/or block diagrams of the method, the device (system), and the computer program product according to the embodiments of the present invention. It should be understood that computer program instructions may be used to implement each process and/or each block in the flowcharts and/or the block diagrams and a combination of a process and/or a block in the flowcharts and/or the block diagrams. These computer program instructions may be provided for a general-purpose computer, a dedicated computer, an embedded processor, or a processor of any other programmable data processing device to generate a machine, so that the instructions executed by a computer or a processor of any other programmable data processing device generate an apparatus for implementing a specific function in one or more processes in the flowcharts and/or in one or more blocks in the block diagrams.

[0102] These computer program instructions may also be stored in a computer readable memory that can instruct the computer or any other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory generate an artifact that includes an instruction apparatus. The instruction apparatus implements a specific function in one or more processes in the flowcharts and/or in one or more blocks in the block diagrams.

[0103] These computer program instructions may also be loaded onto a computer or another programmable data processing device, so that a series of operations and steps are performed on the computer or another programmable device, thereby generating computer-implemented processing. Therefore, the instructions executed on the computer or another programmable device provide steps for implementing a specific function in one or more processes in the flowcharts and/or in one or more blocks in the block diagrams.

[0104] In this specification, specific embodiments are used to describe the principle and implementations of the present invention, and the description of the embodiments is only intended to help understand the method and core idea of the present invention. Meanwhile, a person of ordinary skill in the art may, based on the idea of the present invention, make modifications with respect to the specific implementations and the application scope. Therefore, the content of this specification shall not be construed as a limitation to the present invention.

[0105] Those of ordinary skill in the art will understand that the embodiments described herein are intended to help readers understand the principles of the present invention, and it should be understood that the protection scope of the present invention is not limited to such special statements and embodiments. Those of ordinary skill in the art may make other various specific modifications and combinations according to the technical teachings disclosed in the present invention without departing from the essence of the present invention, and such modifications and combinations still fall within the protection scope of the present invention.