Data-model coupling driven drilling leakage pattern recognition system and method

Abstract

The present invention relates to a drilling spill pattern recognition system and method driven by data-model coupling, which belongs to the technical field of pattern recognition of drilling conditions; it includes a flow parameter acquisition system and a data intelligent analysis system; the flow parameter acquisition system is used to collect flow parameter data during the drilling process; the data intelligent analysis system is used to process the flow parameter data during the drilling process, perform noise reduction and pattern recognition, analyze the change in wellbore pressure after the accident well overflows, and intuitively display the probability of overflow occurrence. The present invention uses the DTW algorithm for pattern matching, which can quickly identify abnormal conditions before the accident occurs and deteriorates, and effectively shorten the spill diagnosis time.

Description

Technical Field

[0001] The present invention relates to the technical field of pattern recognition of drilling conditions, and in particular to a drilling leakage pattern recognition system and method driven by data-model coupling. Background Art

[0002] Overflow refers to the uncontrolled entry of fluid from the formation into the wellbore driven by pressure differential. Loss refers to the uncontrolled leakage of fluid from the wellbore through the pipe or cracks in the wellbore wall. Overflow will increase the risk of accidents during drilling and cause huge economic and environmental losses. Early diagnosis and evaluation of leakage can effectively ensure the safety and efficiency of drilling and completion operations.

[0003] At present, the commonly used drilling leakage diagnosis technology usually uses the threshold method for judgment, that is, when the outlet flow or mud pool increment exceeds the preset value (such as 1m³), an alarm is issued. Affected by the objective environment and the accuracy of the measurement data, the application of traditional leakage diagnosis methods in deepwater drilling has certain limitations: (1) The threshold-based leakage judgment method is difficult to balance between the accuracy of prediction and the reduction of false alarm rate. In practical applications, a large threshold can reduce the false alarm rate, but it often means that the prediction accuracy will be reduced, while a small threshold is just the opposite. Therefore, under high noise conditions, the prediction accuracy is low, the measurement value accuracy is low, or the prediction effect is poor under complex environmental conditions. (2) Reliance on high-precision computing software and measurement methods. Comparing the calculated data with the measured data values ​​to determine whether overflow has occurred requires high-precision computing software and measurement data, otherwise it will cause a high false alarm rate. However, the complexity of environmental factors leads to uncertainty in the input data and model parameters of multiphase flow calculations, which rely on very accurate automatic calibration. Rely on a large amount of historical leakage data or prior information, such as threshold setting, statistics of leakage events, etc. However, in deepwater drilling, there is usually limited well data and little historical data on spills. Summary of the invention

[0004] In view of the shortcomings of the prior art, the present invention proposes a drilling spill pattern recognition system driven by data-model coupling;

[0005] The present invention also proposes a drilling leakage pattern recognition method driven by data-model coupling.

[0006] The present invention proposes a drilling spill pattern recognition method driven by multi-source data-physical model coupling, establishes a standard pattern for extracting flow parameters from a wellbore multiphase flow model, and uses the DTW algorithm to match the measured data with the standard pattern. This shortens the spill diagnosis time and effectively reduces the false alarm rate, preventing spills from developing into major drilling accidents due to the inability to diagnose in time, and provides protection for deepwater drilling and completion safety.

[0007] Terminology explanation:

[0008] PWD: The full name is Pressure While Drilling, and its Chinese name is the annular pressure measurement tool while drilling.

[0009] Box plot method: It is a method used to detect and clean outliers in a data set. It helps identify the degree of dispersion and symmetry of the data by showing the distribution of the data.

[0010] The technical solution of the present invention is:

[0011] A drilling spill pattern recognition system driven by data-model coupling, including a flow parameter acquisition system and a data intelligent analysis system;

[0012] The flow parameter acquisition system is used to acquire flow parameter data during the drilling process;

[0013] The data intelligent analysis system is used to process flow parameter data during the drilling process, perform noise reduction and pattern recognition, analyze the change in wellbore pressure after overflow of the accident well, and intuitively display the probability of overflow occurrence.

[0014] According to the preferred embodiment of the present invention, the flow parameter acquisition system includes: a downhole PWD tool, a standpipe pressure gauge, a casing pressure gauge, an inlet flow meter, an outlet flow meter, and a mud pool level measuring meter;

[0015] The downhole PWD tool is connected to the logging tool and drill pipe at the top and the drill bit at the bottom. During the drilling process, the downhole PWD tool measures the downhole pressure data in real time.

[0016] The standpipe pressure gauge is installed at the top of the drill pipe to measure the flow pressure change after the drilling fluid is pumped into the standpipe;

[0017] The casing pressure gauge is installed at the front end of the throttle manifold of the wellhead throttle pipeline to measure the fluid pressure change when the drilling fluid is returned from the wellbore to the platform;

[0018] The inlet flow meter is installed between the drilling pump and the drill pipe to measure the flow rate of drilling fluid pumped into the drill pipe;

[0019] The outlet flow meter is connected to the choke manifold and the wellhead spool to monitor the flow rate changes of the fluid returned from the annulus;

[0020] The mud pool level meter is used to measure the changes in the mud pool liquid level during drilling.

[0021] Preferably, according to the present invention, the data intelligent analysis system comprises: a data acquisition pipeline, a computer, a display, and an alarm;

[0022] The data acquisition pipeline connects the computer and the flow parameter acquisition system to transmit the measurement signal of the flow data; the computer is used to analyze the change pattern of the wellbore pressure after the overflow of the current accident well, and to perform noise reduction and pattern recognition on the collected data; the display is connected to the computer to intuitively output the curve of the overflow probability during the drilling process.

[0023] The data-model coupling driven drilling leakage pattern recognition method is implemented by the above-mentioned drilling leakage pattern recognition system, and includes:

[0024] Perform data cleaning, Lagrange interpolation, and normalization on various types of measurement data to generate standard data sets; including:

[0025] Step 1: Use the box plot method to clean outliers;

[0026] Step 2: For the missing values ​​of some eigenvalue information in the sample data, Lagrange interpolation method is used to fill them;

[0027] Step 3: Normalization processing;

[0028] Construct a Kalman filter to reduce the noise of the obtained standard data set;

[0029] Establish a wellbore multiphase flow model to extract the standard mode of flow parameters under abnormal conditions of overflow and leakage as well as normal conditions of drilling, pumping, drilling and drilling;

[0030] The LB_Keogh method based on the constrained search space and lower bound method is combined with the DTW algorithm based on dynamic programming to accelerate the time series similarity search, and the PWD measurement data is pattern matched using the obtained standard patterns of different working conditions.

[0031] According to the DTW distance, the final probability of the fault event is obtained.

[0032] Preferably, according to the present invention, a time series of PWD real-time second data in the drilling process, i.e., multi-type measurement data, is established, and PWD real-time second data is acquired at a certain time interval, and the time series of PWD real-time second data in the drilling event is obtained as follows: ;in, P It is the time series of PWD real-time second data; p i For the i The bottom hole pressure at each time point, m is the length of the time series.

[0033] Preferably, according to the present invention, a box plot method is used to clean outliers; comprising:

[0034] Step 1-1: Calculate the first quartile Q based on the time series of PWD real-time second data in the drilling event1 and the third quartile Q 3 , and the interquartile range IQR = Q 3 -Q 1 ;

[0035] Step 1-2: Define the upper edge Up_Fence and the lower edge Low_Fence:

[0036] ;

[0037] ;

[0038] Where k is the coefficient used to determine outliers;

[0039] Step 1-3: Define the data range:

[0040] According to the upper edge Up_Fence, the lower edge Low_Fence, the first quartile Q 1 and the third quartile Q 3 , the data is divided into three parts:

[0041] The range of the box is: Q 1 ≤p i ≤Q 3 ;

[0042] The box whisker range is: Low_Fence≤pi≤Q 1 and Q 1 ≤pi≤Up_Fence;

[0043] Outlier points: data points that are smaller than Low_Fence and larger than Up_Fence;

[0044] Step 1-4: Delete outlier points.

[0045] Preferably, according to the present invention, for the values ​​of missing eigenvalue information in the time series of sample data, i.e., PWD real-time second data, Lagrange interpolation method is used to fill in the missing values; including:

[0046] Step 2-1: The Lagrange polynomial is a cubic polynomial; assume that the time series of PWD changes in drilling events is : ;in,( t i ,p i )for t i The bottom hole pressure at the moment, m is the length of the time series;

[0047] Step 2-2: Set the time point at which missing values ​​need to be filled, i.e. the target point t i ;

[0048] Step 2-3: At the target point Select 4 evenly distributed data points near to form a set S : ;

[0049] Step 2-4: For each interpolation point (t i ,p i ) , define the Lagrangian basis function as L i ( t ):

[0050] ;

[0051] in, t For time t , t i and t j is the first i and j The time of each data point;

[0052] Step 2-5: Construct the interpolation polynomial using the Lagrange basis functions and the 4 known data points :

[0053] ;

[0054] in, L i ( t ) is the Lagrangian basis function.

[0055] According to the preferred embodiment of the present invention, the normalization process comprises:

[0056] Assume that the one-dimensional time series after Lagrange interpolation is ;in, is each element in the sequence; the data at each time point is scaled to between [0,1] as follows:

[0057] ;

[0058] in, min ( P ) is a one-dimensional time series P （3） The minimum value in max (P ) is a one-dimensional time series P （3） The maximum value in ;

[0059] After normalization, we finally get a standard data set: .

[0060] Preferably, according to the present invention, a Kalman filter is constructed to reduce noise on the obtained standard data set; comprising:

[0061] Step 4-1: Build the model

[0062] Assume that the one-dimensional time series vector of PWD real-time second data after Lagrange interpolation is P (2) for: ;

[0063] in, Representative m The bottom hole pressure at a time node; that is:

[0064] ;

[0065] in, t k Indicates k The time of the sampling point, n k Representative in the k The noise at the sampling point is obtained by setting the real value of the bottom hole pressure at a certain moment as the state variable and setting the bottom hole pressure at a certain moment with noise as the observed value. k The actual value of bottom hole pressure at the moment is x k , establish the following state space model:

[0066] Step 4-1-1: Establish the state equation:

[0067] ;

[0068] In the formula, x k It's time k Bottom hole pressure value at 10000mA; A is the state transfer matrix, representing k-1 The bottom hole pressure value at time k Moment j Conversion of bottom hole pressure values; B is the influence matrix of the control input; u k It's time k The control input, w k yes kThe prediction noise of bottom hole pressure at any moment;

[0069] Step 4-1-2: Establish the observation equation:

[0070] ;

[0071] In the formula, y k It's time k Observed value of bottom hole pressure at 1000 m / s; H is the observation matrix, v k yes k The prediction noise of bottom hole pressure at any moment;

[0072] Step 4-2: Kalman filter process

[0073] Step 4-2-1: Initialization

[0074] Initialize the system status x 0 , that is: set the initial state estimate and covariance matrix required by the Kalman filter; take the first numerical point in the normalized PWD real-time second data time series standard data set;

[0075] State transition matrix A :Get the state transfer matrix A is a scalar 1;

[0076] State covariance matrix P 0 : Get the state covariance matrix P 0 is a scalar 1;

[0077] Prediction noise covariance matrix Q : Get the prediction noise covariance matrix Q is a scalar 0.05;

[0078] Covariance matrix of observation noise R : Take the covariance matrix of the observation noise R is the scalar 5;

[0079] Step 4-2-2: Prediction step

[0080] According to the state equation, the state and state covariance of the bottom hole pressure at the next moment are predicted:

[0081] ;

[0082] ;

[0083] In the formula, is the estimated value of the predicted bottom hole pressure PWD, is the predicted bottom hole pressure PWD state covariance matrix, Q is the covariance matrix of the prediction noise, u k It's time k Control input of

[0084] Step 4-2-3: Update steps

[0085] According to the observation equation, the predicted value is obtained by combining the measured value of bottom hole pressure PWD Make corrections to get the optimal estimate:

[0086] ;

[0087] ;

[0088] ;

[0089] In the formula, K k is the Kalman gain, R is the covariance matrix of the observation noise, x k is the optimal state estimate, P k is the optimal state covariance matrix, I A matrix of 1;

[0090] Step 4-2-4: Repeat the prediction and update steps

[0091] according to j The one-dimensional time series of bottom hole pressure PWD is predicted and updated repeatedly in sequence, and the state estimation value and state covariance matrix are continuously updated to obtain the optimal state estimation.

[0092] Preferably, according to the present invention, a wellbore multiphase flow model is established to extract the standard mode of flow parameters under abnormal conditions of overflow and leakage and normal conditions of drilling, pumping, drilling and drilling;

[0093] Step 5-1: Establish the mass conservation equation of gas and liquid phase:

[0094] ;

[0095] ;

[0096] The momentum conservation equation for multiphase flow in the wellbore is established as:

[0097] ;

[0098] Where:A is the cross-sectional area of ​​the annulus, m 2 ; α g is the void fraction of the gas, m 3 / m 3 ; α l is the porosity of drilling fluid, m 3 / m 3 ; p is the fluid pressure, MPa; θ is the well inclination angle, rad; K is the reservoir permeability, m 2 ; p res is the reservoir pressure, Pa; p wf is the fluid pressure at the bottom of the well, Pa; S is the skin factor of the reservoir; t is the overflow time, s; s is the distance, m; γ is the Euler constant, γ=0.5772; φ is the reservoir porosity; μ is the viscosity of the gas in the reservoir, Pa∙s; r w is the wellbore radius, m. ρ g is the gas phase density, kg / m 3 : ρ l is the liquid density, kg / m 3 , v g is the gas phase velocity, m / s; v l is the liquid phase velocity, m / s;

[0099] Step 5-2: Set initial and boundary conditions, simulate abnormal conditions of overflow and leakage, and normal conditions of drilling, pumping, drilling and drilling, and obtain the variation law of flow parameters; in overflow conditions:

[0100] Step 5-2-1: Set initial conditions; during normal drilling, when drilling to the top of the reservoir and no formation fluid flows in:

[0101] ;

[0102] ;

[0103] ;

[0104] ;

[0105] ;

[0106] ;

[0107] ;

[0108] ;

[0109] in, Q m is the displacement of the mud pump, m 3 / s; q c is the discharge volume of cuttings, m 3 / s; S is the well depth, m; E H is the hydrate volume fraction; E g is the gas volume fraction, 0-1 ; E c is the volume fraction of cuttings; V sc is the actual flow velocity of cuttings, m / s; E m is the volume fraction of drilling fluid, 0-1; v m is the superficial velocity of drilling fluid, m / s; v sm is the actual flow rate of drilling fluid, m / s; ρ c is the density of rock cuttings, kg / m 3 ; v c is the apparent velocity of cuttings, m / s.

[0110] The phase volume fraction is determined by the above equations. By solving the above equations, the pressure, velocity, and phase volume fraction parameters of the multiphase flow along the well depth at the initial moment are determined;

[0111] Step 5-2-2: Set boundary conditions; the boundaries are:

[0112] Drilling overflow conditions:

[0113] ;

[0114] in, p is pressure, Pa; p s is the wellhead casing pressure, Pa; q g is the volume flow rate of the invading gas, m 3 / s ; q c is the discharge volume of cuttings, m 3 / s.r H is the hydrate phase transition rate, m 3 / s;

[0115] Set simulation parameters, including drilling fluid density, plastic viscosity, dynamic shear force, reservoir permeability, initial pressure difference, displacement,

[0116] Step 5-3: Change the simulation parameters, conduct a large number of numerical simulations, normalize the simulation results, and obtain the time series of the bottom hole pressure and mud pool increment under abnormal conditions of overflow and leakage and normal conditions of drilling, pumping, drilling and drilling, that is, the basic change pattern.

[0117] Preferably, according to the present invention, the LB_Keogh method based on the constrained search space and the lower bound method is combined with the DTW algorithm based on dynamic programming to accelerate the time series similarity search, and the PWD measurement data is pattern matched using the obtained standard modes of different working conditions; including:

[0118] Step 6-1: Initialize the minimum distance variable Min_distance:

[0119] Min_distance is used to store the minimum DTW distance between the sequence to be matched and each standard pattern sequence. For each standard pattern sequence, the operations from step 6-2 to step 6-3 are performed on it in sequence;

[0120] Step 6-2: Calculate the LB_Keogh lower bound distance LB between a standard pattern sequence Q and the sequence to be matched C; it is known that the lower bound distance LB is always smaller than the DTW distance. If LB>Best_Distance, the current calculation must not be the optimal match; if the lower bound distance is smaller than the currently found minimum distance Min_distance, continue to execute step 6-3, otherwise skip the current standard pattern sequence, select the next standard pattern sequence and repeat step 6-2; Best_Distance is a variable used to track the current optimal distance between the sequence to be matched and the standard pattern sequence;

[0121] Step 6-2-1: Assume that the time series of length m of the first standard variation pattern is: ;

[0122] Input the time series of PWD measurement data with length n to be matched:

[0123] ;

[0124] Step 6-2-2: Initialize the distance matrix and calculate the local distance

[0125] Create an n×m matrix D to store the local distance between each pair of data points between the two sequences;

[0126] Step 6-2-3: Set constraints

[0127] Think of global or local constraints as indices of constrained warp paths D(i,j) such that j−r≤i≤j+r, where r is a term that defines the reach or allowed range of the warp at a given point in the sequence;

[0128] Step 6-2-4: Calculate the upper and lower boundaries and the distance between the lower and upper boundaries

[0129] According to the constraints, the upper boundary of the sequence Q near the adjacent time window r is set to U i , the lower boundary is L i , forming an envelope that surrounds the sequence Q above and below, and the upper and lower boundaries are defined as follows:

[0130] ;

[0131] ;

[0132] Among them, for each time node i, there is L i ≤q i ≤U i ;

[0133] According to the upper and lower boundaries, set the lower boundary distance LB to:

[0134] ;

[0135] Step 6-2-5: Compare the lower bound distance LB with the minimum distance Min_distance;

[0136] If the lower bound distance of the standard pattern sequence Q is less than the currently found minimum distance Min_distance, continue to execute step 6-3, otherwise skip the current standard pattern sequence, select the next standard pattern sequence and repeat step 6-2;

[0137] Step 6-3: Calculate the actual DTW distance between the standard pattern sequence Q and the sequence to be matched C. If the actual DTW distance is less than the currently found minimum distance Min_distance, update Min_distance to the actual DTW distance and record the index of the standard pattern sequence Q as index_Best_match;

[0138] Step 6-3-1: Calculate DTW distance

[0139] Dynamic programming calculates the minimum cumulative distance and calculates the cumulative distance matrix A,in, A [ i , j ] indicates a sequence Q Middle i Points and sequences C Middle j The minimum cumulative distance between points; the specific calculation method is as follows:

[0140] ;

[0141] The final A(n,m]) is the DTW distance D between the two sequences. dtw ;

[0142] Step 6-3-2: Compare DTW distance with minimum distance Min_distance

[0143] If the actual DTW distance is less than the currently found minimum distance Min_distance, then Min_distance is updated to the actual DTW distance, and the index of the standard pattern sequence Q is recorded as index_Best_match.

[0144] Step 6-4: Return the minimum distance Min_distance and the corresponding sequence index index_Best_match. Index_Best_match points to the standard pattern sequence Q that is most similar to the sequence C to be matched.

[0145] Preferably, according to the present invention, the final probability of occurrence of a fault event is obtained according to the DTW distance; including:

[0146] Assume that the similarity probability between the optimal solution of an event and the basic pattern should be equal to the DTW distance D dtw Inversely proportional, the similarity probability vector of the basic mode is obtained as:

[0147] ;

[0148] In the formula, M Indicates the number of basic patterns that may cause changes in wellbore flow data; through real-time calculation and matching analysis, the change pattern of fault event probability over time is obtained;

[0149] is the similarity probability of each basic pattern; this vector represents the similarity probability between each basic pattern and the sequence to be matched; is the DTW distance between the sequence to be matched and the i-th standard pattern; i is the number of the standard pattern, ranging from 1 to M; is the inverse of the DTW distance of the iiith standard pattern, used to indicate the similarity with the standard pattern; is the sum of the inverses of the DTW distances of all standard modes and is a normalization factor used to ensure that all sums are 1, thus forming a probability distribution.

[0150] The beneficial effects of the present invention are:

[0151] 1. Shorten the leak diagnosis time: The leak pattern recognition method driven by multi-source data-physical model coupling proposed in the present invention uses the DTW algorithm for pattern matching. It can quickly identify abnormal conditions before an accident occurs and worsens, effectively shortening the leak diagnosis time.

[0152] 2. Reduce false alarm rate: Through comprehensive analysis of bottom hole pressure, flow rate and mud pool increment data, the present invention can accurately determine leakage events, reduce false alarm rate, avoid invalid well control operations caused by false alarms, and ensure the safety of the drilling process.

[0153] 3. Intelligent monitoring and real-time output: This method realizes intelligent monitoring of the spill process and can output the probability of spill occurrence in real time, helping field engineers to quickly understand the development of the spill and provide accurate decision support for well control operations, thereby improving the efficiency and safety of offshore drilling. BRIEF DESCRIPTION OF THE DRAWINGS

[0154] Figure 1 The structural block diagram of the drilling spill pattern recognition system driven by data-model coupling;

[0155] Figure 2 A flowchart of a drilling spill pattern recognition method driven by data-model coupling;

[0156] Figure 3 This is a schematic diagram of the Kalman filter effect;

[0157] Figure 4 This is a schematic diagram of the simulated well body structure;

[0158] Figure 5 It is a schematic diagram of the PWD variation curve during drilling overflow;

[0159] Figure 6 It is a schematic diagram of matching the PWD measurement data time series with the standard mode;

[0160] Figure 7 It is a matching diagram of PWD data;

[0161] Figure 8 It is a schematic diagram of the measured data under overflow conditions;

[0162] Fig. 9It is a schematic diagram of the PWD measured data curve after preprocessing;

[0163] Fig.10 is a schematic diagram of pattern matching for the overflow standard pattern;

[0164] Among them, 1. Downhole PWD tools, 2. Wellbore, 3. Annulus, 4. Wellhead spool, 5. Standpipe pressure gauge, 6. Inlet flow meter, 7. Drilling pump, 8. Casing pressure gauge, 9. Choke manifold, 10. Outlet flow meter, 11. Mud pool level gauge, 12. Data acquisition pipeline, 13. Computer, 14. Display, 15. Alarm. DETAILED DESCRIPTION

[0165] For a better understanding of the present invention, the following is further described in conjunction with specific examples and the accompanying drawings. The experimental methods used in the examples are conventional methods unless otherwise specified. The equipment used in the examples, etc., can be obtained from commercial sources unless otherwise specified.

[0166] Example 1

[0167] A drilling spill pattern recognition system driven by data-model coupling, including a flow parameter acquisition system and a data intelligent analysis system;

[0168] Flow parameter acquisition system, used to collect flow parameter data during drilling;

[0169] The data intelligent analysis system is used to process the flow parameter data during the drilling process, perform noise reduction and pattern recognition, analyze the change in the wellbore 2 pressure after the accident well overflows, and intuitively display the probability of overflow occurrence through the display 14.

[0170] Example 2

[0171] The data-model coupling driven drilling leakage pattern recognition system according to Example 1 is different in that:

[0172] like Figure 1 As shown, the flow parameter acquisition system includes: a downhole PWD tool 1, a standpipe pressure gauge 5, a casing pressure gauge 8, an inlet flow meter 6, an outlet flow meter 10, and a mud pool level measuring meter 11;

[0173] The downhole PWD tool 1 is connected to the logging tool and the drill pipe at the top and the drill bit at the bottom. During the drilling process, the downhole PWD tool 1 measures the downhole pressure data in real time.

[0174] The standpipe pressure gauge 5 is installed at the top of the drill pipe to measure the flow pressure change after the drilling fluid is pumped into the standpipe;

[0175] The casing pressure gauge 8 is installed at the front end of the throttle manifold 9 of the wellhead throttle pipeline to measure the fluid pressure change of the drilling fluid returned from the wellbore 2 to the platform; the platform is a structure for drilling operations, usually located at sea or on land.

[0176] The inlet flow meter 6 is installed between the drilling pump 7 and the drill pipe to measure the flow rate of drilling fluid pumped into the drill pipe;

[0177] An outlet flow meter 10 is connected to the choke manifold 9 and the wellhead spool 4 and is used to monitor the change in the fluid flow rate returned from the annulus 3;

[0178] The mud pool level meter 11 is used to measure the change of the mud pool liquid level during the drilling process.

[0179] The data intelligent analysis system includes: a data acquisition pipeline 12, a computer 13, a display 14, and an alarm 15;

[0180] The data acquisition pipeline 12 connects the computer 13 with the flow parameter acquisition system to transmit the measurement signal of the flow data; the computer 13 is used to analyze the pressure change pattern of the wellbore 2 after the overflow of the current accident well, and to perform noise reduction and pattern recognition on the collected data; the display 14 is connected to the computer 13 to intuitively output the curve of the overflow probability during the drilling process.

[0181] During the drilling process, the drilling fluid, driven by the mud pump, flows through the mud pool, the drilling pump 7, into the outflow pipeline, the inlet flow meter 6, the rotary control head and the wellhead spool 4, and is injected into the drill pipe; the drilling fluid flows from top to bottom in the drill pipe, returns from the drill bit water hole to the annulus 3, and gradually flows from bottom to top along the wellbore 2 carrying the broken rock cuttings; after the gas returns to the wellhead and separates the solid phase, it flows through the throttling manifold 9, the outlet flow meter 10, the drilling fluid outlet pipeline, and is re-injected into the mud pool;

[0182] If the fluid pressure at the bottom of the well is lower than the formation pressure, the reservoir fluid, driven by the pressure difference, invades the wellbore 2 and returns to the wellhead through the gas-liquid two-phase flow; the invasion of gas causes the hydrostatic pressure of the gas-liquid two-phase in the wellbore 2 to decrease, the bottom hole pressure to decrease, and the riser pressure to decrease; in addition, the gas phase enters the wellbore 2, the outlet flow rate of the gas-liquid phase gradually increases, and the total amount of drilling fluid in the mud pool gradually increases; the specific measured parameters are: bottom hole pressure PWD, mud pool increment, inlet and outlet flow rates, inlet and outlet temperatures and pressures, etc.

[0183] Example 3

[0184] The data-model coupling driven drilling leakage pattern recognition method is implemented by the drilling leakage pattern recognition system described in Example 1 or 2, such as Figure 2 As shown, including:

[0185] Perform data cleaning, Lagrange interpolation, and normalization on various types of measurement data to generate standard data sets; including:

[0186] Step 1: Use the box plot method to clean outliers;

[0187] Step 2: For the missing values ​​of some eigenvalue information in the sample data, Lagrange interpolation method is used to fill them;

[0188] Step 3: Normalization processing;

[0189] Construct a Kalman filter to reduce the noise of the obtained standard data set;

[0190] Establish a wellbore multiphase flow model to extract the standard mode of flow parameters under abnormal conditions of overflow and leakage as well as normal conditions of drilling, pumping, drilling and drilling;

[0191] The LB_Keogh method based on the constrained search space and lower bound method is combined with the DTW algorithm based on dynamic programming to accelerate the time series similarity search, and the PWD measurement data is pattern matched using the obtained standard patterns of different working conditions.

[0192] According to the DTW distance, the final probability of the fault event is obtained.

[0193] Example 4

[0194] The data-model coupling driven drilling leakage pattern recognition method according to Example 3 is different in that:

[0195] Establish the time series of PWD real-time second data during the drilling process, that is, multi-type measurement data. The time series of PWD real-time second data is "multi-type measurement data" with pressure data as the main type. Some PWD tools can measure drilling fluid density, temperature, and flow rate. Obtain PWD real-time second data at a certain time interval (5 seconds), and the time series of PWD real-time second data in the drilling event is: ;in, P It is the time series of PWD real-time second data; p i For the i The bottom hole pressure at each time point, m is the length of the time series.

[0196] Use the box plot method to clean outliers; including:

[0197] Step 1-1: Based on the time series of PWD real-time second data in drilling events , calculate the first quartile Q 1 and the third quartile Q 3 , and the interquartile range IQR = Q3 -Q 1 ;

[0198] Step 1-2: Define the upper edge Up_Fence and the lower edge Low_Fence:

[0199] ;

[0200] ;

[0201] Wherein, k is a coefficient used to determine outliers; the value of k in the present invention is 1.3.

[0202] Step 1-3: Define the data range:

[0203] According to the upper edge Up_Fence, the lower edge Low_Fence, the first quartile Q 1 and the third quartile Q 3 , the data is divided into three parts:

[0204] The range of the box is: Q 1 ≤p i ≤Q 3 ;

[0205] The box whisker range is: Low_Fence≤pi≤Q 1 and Q 1 ≤pi≤Up_Fence;

[0206] Outlier points: data points that are smaller than Low_Fence and larger than Up_Fence;

[0207] Step 1-4: Delete outlier points.

[0208] For the missing values ​​of some eigenvalue information in the time series of sample data, i.e. PWD real-time second data, the Lagrange interpolation method is used to fill in the missing values, including:

[0209] Step 2-1: The Lagrange polynomial is a cubic polynomial; assume that the time series of PWD changes in drilling events is : ;in,( t i ,p i )for t i The bottom hole pressure at the moment, m is the length of the time series;

[0210] Step 2-2: Set the time point at which missing values ​​need to be filled, i.e. the target point t i ;

[0211] Step 2-3: At the target point Select 4 evenly distributed data points near to form a set S : ;

[0212] Step 2-4: For each interpolation point (t i ,p i ) , define the Lagrangian basis function as L i ( t ):

[0213] ;

[0214] in, t For time t , t i and t j is the first i and j The time of each data point;

[0215] Step 2-5: Construct the interpolation polynomial using the Lagrange basis functions and the 4 known data points :

[0216] ;

[0217] in, L i ( t ) is the Lagrangian basis function.

[0218] Normalization processing; including:

[0219] Assume that the one-dimensional time series after Lagrange interpolation is ;in, is each element in the sequence; the data at each time point is scaled to between [0,1] as follows:

[0220] ;

[0221] in, min ( P ) is a one-dimensional time series P （3） The minimum value in max ( P ) is a one-dimensional time series P （3） The maximum value in ;

[0222] After normalization, we finally get a standard data set: .

[0223] Construct a Kalman filter to reduce the noise of the obtained standard data set; including:

[0224] Step 4-1: Build the model

[0225] Assume that the one-dimensional time series vector of PWD real-time second data after Lagrange interpolation is P (2) for: ;

[0226] in, Representative m The bottom hole pressure at a time node; that is:

[0227] ;

[0228] in, t k Indicates k The time of the sampling point, n k Representative in the k The noise at the sampling point is obtained by setting the real value of the bottom hole pressure at a certain moment as the state variable and setting the bottom hole pressure at a certain moment with noise as the observed value. k The actual value of bottom hole pressure at the moment is x k , establish the following state space model:

[0229] The state space model includes the state equation and the observation equation. These two equations together describe the dynamic behavior and observation process of the system and are the basis for the application of the Kalman filter.

[0230] Step 4-1-1: Establish the state equation:

[0231] ;

[0232] In the formula, x k It's time k Bottom hole pressure value at 10000mA; A is the state transfer matrix, representing k-1 The bottom hole pressure value at time k Moment j Conversion of bottom hole pressure values; B is the influence matrix of the control input; u k It's time k Since there is no control input for time when acquiring bottom hole pressure data, the present invention uk The value of is 0. w k yes k The predicted noise of bottom hole pressure at any time; it can be determined through experiments or instrument measurement analysis.

[0233] Step 4-1-2: Establish the observation equation:

[0234] ;

[0235] In the formula, y k It's time k Observed value of bottom hole pressure at 1000 m / s; H is the observation matrix, which describes how to obtain the observation value from the bottom hole pressure data. Since the collected PWD data set is a one-dimensional time series, H The value of is scalar 1. v k yes k The predicted noise of bottom hole pressure at any time; it can be determined through experiments or instrument measurement analysis.

[0236] Step 4-2: Kalman filter process

[0237] Step 4-2-1: Initialization

[0238] Initialize the system status x 0 , that is: set the initial state estimate and covariance matrix required by the Kalman filter; take the first numerical point in the normalized PWD real-time second data time series standard data set;

[0239] State transition matrix A :Get the state transfer matrix A It is a scalar 1; PWD data is one-dimensional time series data and can be considered as a scalar.

[0240] State covariance matrix P 0 : Get the state covariance matrix P 0 is a scalar 1; the initial covariance matrix has little effect on the filtering effect of the data, so a non-zero value can be used.

[0241] Prediction noise covariance matrix Q : Get the prediction noise covariance matrix Q It is a scalar value of 0.05; generally Q The value is inversely proportional to the filtering effect, but Q Too small will cause a certain deviation between the filtering result curve and the initial curve. Therefore, multiple experimental simulations can be used to obtain an appropriate Q value.

[0242] Covariance matrix of observation noise R : Take the covariance matrix of the observation noise R is a scalar 5; analysis found that R If the value of is too large or too small, the filtering effect will be reduced. Therefore, it is necessary to conduct multiple experimental simulations to determine the covariance matrix of the observation noise. R The value of .

[0243] Step 4-2-2: Prediction step

[0244] According to the state equation, the state and state covariance of the bottom hole pressure at the next moment are predicted:

[0245] ;

[0246] ;

[0247] In the formula, is the estimated value of the predicted bottom hole pressure PWD, is the predicted bottom hole pressure PWD state covariance matrix, Q is the covariance matrix of the prediction noise, u k It's time k Since there is no control input system for time when collecting bottom hole pressure PWD, the present invention u k-1 The value of is 0.

[0248] Step 4-2-3: Update steps

[0249] According to the observation equation, the predicted value is obtained by combining the measured value of bottom hole pressure PWD Make corrections to get the optimal estimate:

[0250] ;

[0251] ;

[0252] ;

[0253] In the formula, K k is the Kalman gain, R is the covariance matrix of the observation noise, x k is the optimal state estimate, P k is the optimal state covariance matrix, I The matrix of 1, in the present invention I =1;

[0254] Step 4-2-4: Repeat the prediction and update steps

[0255] according to j The one-dimensional time series of bottom hole pressure PWD is predicted and updated repeatedly in sequence, and the state estimation value and state covariance matrix are continuously updated to obtain the optimal state estimation.

[0256] The Kalman filter will dynamically fuse the predicted value and observed value of the bottom hole pressure PWD data at each moment according to the system model and the measured value, thereby achieving the denoising of the one-dimensional time series data of the bottom hole pressure PWD and the state estimation of the bottom hole pressure PWD at each moment.

[0257] The noise reduction and smoothing effects of Kalman filtering are as follows: Figure 3 shown.

[0258] Establish a wellbore multiphase flow model to extract the standard mode of flow parameters under abnormal conditions of overflow and leakage as well as normal conditions of drilling, pumping, drilling and drilling;

[0259] Step 5-1: The flow parameters of the gas and liquid fluids (such as phase distribution, velocity and pressure) can be obtained by solving the continuity equation and momentum conservation equation. The mass conservation equation of the gas and liquid phase is established as:

[0260] ;

[0261] ;

[0262] The momentum conservation equation for multiphase flow in wellbore 2 is established as:

[0263] ;

[0264] Where: A is the cross-sectional area of ​​annulus 3, m 2 ; α g is the void fraction of the gas, m 3 / m 3 ; α l is the porosity of drilling fluid, m 3 / m 3 ; p is the fluid pressure, MPa; θ is the well inclination angle, rad; K is the reservoir permeability, m 2 ; p res is the reservoir pressure, Pa; p wf is the fluid pressure at the bottom of the well, Pa; S is the skin factor of the reservoir;t is the overflow time, s; s is the distance, m; γ is the Euler constant, γ=0.5772; φ is the reservoir porosity; μ is the viscosity of the gas in the reservoir, Pa∙s; r w is the radius of wellbore 2, m. ρ g is the gas phase density, kg / m 3 : ρ l is the liquid density, kg / m 3 , v g is the gas phase velocity, m / s; v l is the liquid phase velocity, m / s;

[0265] Step 5-2: Set initial and boundary conditions, simulate abnormal conditions of overflow and leakage, and normal conditions of drilling, pumping, drilling and drilling, and obtain the variation law of flow parameters; in overflow conditions:

[0266] Step 5-2-1: Set initial conditions; during normal drilling, when drilling to the top of the reservoir and no formation fluid flows in:

[0267] ;

[0268] ;

[0269] ;

[0270] ;

[0271] ;

[0272] ;

[0273] ;

[0274] ;

[0275] in, Q m is the displacement of the mud pump, m 3 / s; q c is the discharge volume of cuttings, m 3 / s; S is the well depth, m; E H is the hydrate volume fraction; E gis the gas volume fraction, 0-1 ; E c is the volume fraction of cuttings; V sc is the actual flow velocity of cuttings, m / s; E m is the volume fraction of drilling fluid, 0-1; v m is the superficial velocity of drilling fluid, m / s; v sm is the actual flow rate of drilling fluid, m / s; ρ c is the density of rock cuttings, kg / m 3 ; v c is the apparent velocity of cuttings, m / s.

[0276] The phase volume fraction is determined by the above equations. By solving the above equations, the pressure, velocity, and phase volume fraction parameters of the multiphase flow along the well depth at the initial moment are determined;

[0277] Step 5-2-2: Set boundary conditions; The integration of multiphase flow including drill string, drill bit, formation seepage, wellbore annulus, riser and choke line is a nonlinear equation system with multiple equations and multiple unknowns. To obtain a certain solution, boundary conditions are also required. For different engineering stages, the boundary conditions are different. Here the boundaries are:

[0278] Drilling overflow conditions:

[0279] ;

[0280] in, p is pressure, Pa; p s is the wellhead casing pressure, Pa; q g is the volume flow rate of the invading gas, m 3 / s ; q c is the discharge volume of cuttings, m 3 / s. r H is the hydrate phase transition rate, m 3 / s;

[0281] Set simulation parameters, including drilling fluid density, plastic viscosity, dynamic shear force, reservoir permeability, initial pressure difference, displacement,

[0282] The main simulation parameters are set as follows: drilling fluid density, 1.83 g / cm 3 ; plastic viscosity, 49 mPa.s; dynamic shear force, 14 Pa; reservoir permeability, 100 md; initial pressure difference, 1.5 MPa; displacement, 0.015 m3 / s; Well structure such as Figure 4 As shown in the figure, the calculated mud pool increment is Figure 5 shown.

[0283] in, Figure 4 The main parameters of the well depth structure are as follows:

[0284] Riser: top depth 0m, bottom depth 520m.

[0285] Casing: top depth 0m, bottom depth 3801m.

[0286] Naked hole section: diameter 215.9mm.

[0287] The first section of drill pipe: outer diameter 129.7mm, inner diameter 118.6mm; the second section of drill pipe: outer diameter 127mm, inner diameter 101.6mm; the third section of drill pipe: outer diameter 127mm, inner diameter 76mm; the fourth section of drill pipe: outer diameter 127mm, inner diameter 101.6mm; the fifth section of drill collar: outer diameter 171mm, inner diameter 83mm.

[0288] Step 5-3: Change the simulation parameters, conduct a large number of numerical simulations, normalize the simulation results, and obtain the time series of the bottom hole pressure and mud pool increment under abnormal conditions of overflow and leakage and normal conditions of drilling, pumping, drilling and drilling, that is, the basic change pattern.

[0289] By using the multiphase flow model established by the present invention, the change patterns of various flow parameters under different drilling conditions can be analyzed and the standard pattern can be extracted.

[0290] The LB_Keogh method based on the constrained search space and lower bound method is combined with the DTW algorithm based on dynamic programming to accelerate the time series similarity search, and the PWD measurement data is pattern matched using the obtained standard modes of different working conditions; including:

[0291] The present invention adopts the LB_Keogh algorithm, and uses the lower bound method and the constrained search space method to optimize DTW, quickly identifying and excluding those candidate sequences that are far away from the query sequence, avoiding DTW calculations, and thus reducing the complexity of the calculations.

[0292] Step 6-1: Initialize the minimum distance variable Min_distance:

[0293] Min_distance is used to store the minimum DTW distance between the sequence to be matched and each standard pattern sequence. For each standard pattern sequence, the operations from step 6-2 to step 6-3 are performed on it in sequence;

[0294] Step 6-2: Calculate the LB_Keogh lower bound distance LB between a standard pattern sequence Q and the sequence to be matched C; it is known that the lower bound distance LB is always less than the DTW distance. If LB>Best_Distance, the current calculation must not be the optimal match; if the lower bound distance is less than the currently found minimum distance Min_distance, continue to step 6-3, otherwise skip the current standard pattern sequence, select the next standard pattern sequence and repeat step 6-2; Best_Distance is a variable used to track the current optimal distance between the sequence to be matched and the standard pattern sequence; this distance is usually the smallest dynamic time warping (DTW) distance found so far. When matching sequences, Best_Distance helps determine whether to continue calculating more complex DTW distances or terminate the calculation early.

[0295] Step 6-2-1: Assume that the time series of length m of the first standard variation pattern is: ;

[0296] Input the time series of PWD measurement data with length n to be matched:

[0297] ;

[0298] Step 6-2-2: Initialize the distance matrix and calculate the local distance

[0299] Create an n×m matrix D to store the local distance between each pair of data points between the two sequences;

[0300] The present invention uses the Euclidean distance to calculate the local distance D (q i , c j ). D(q i ,c j ) = (q i – c j ) 2 。 The final square root function is omitted in the Euclidean distance definition, which speeds up the actual calculation without affecting the final similarity measure.

[0301] Step 6-2-3: Set constraints

[0302] A global constraint or a local constraint is thought of as an index into a constrained warp path D(i,j) such that j−r≤i≤j+r, where r is a term defining the reach or allowed range of the warp for a given point in the sequence; that is, for each time point j in the time series, the range of allowed indices of the time point i corresponding to it in the time warping matrix is ​​[j - r, j + r].

[0303] Step 6-2-4: Calculate the upper and lower boundaries and the distance between the lower and upper boundaries

[0304] According to the constraints, the upper boundary of the sequence Q near the adjacent time window r is set to U i , the lower boundary is L i , forming an envelope that surrounds the sequence Q above and below, and the upper and lower boundaries are defined as follows:

[0305] ;

[0306] ;

[0307] Among them, for each time node i, there is L i ≤q i ≤U i ;

[0308] According to the upper and lower boundaries, set the lower boundary distance LB to:

[0309] ;

[0310] Step 6-2-5: Compare the lower bound distance LB with the minimum distance Min_distance;

[0311] If the lower bound distance of the standard pattern sequence Q is less than the currently found minimum distance Min_distance, continue to execute step 6-3, otherwise skip the current standard pattern sequence, select the next standard pattern sequence and repeat step 6-2;

[0312] Step 6-3: Calculate the actual DTW distance between the standard pattern sequence Q and the sequence to be matched C. If the actual DTW distance is less than the currently found minimum distance Min_distance, update Min_distance to the actual DTW distance and record the index of the standard pattern sequence Q as index_Best_match;

[0313] Step 6-3-1: Calculate DTW distance

[0314] Dynamic programming calculates the minimum cumulative distance and calculates the cumulative distance matrix A ,in, A [ i , j ] indicates a sequence Q Middle i Points and sequences C Middle j The minimum cumulative distance between points; the specific calculation method is as follows:

[0315] ;

[0316] The final A(n,m]) is the DTW distance D between the two sequences. dtw ;

[0317] Step 6-3-2: Compare DTW distance with minimum distance Min_distance

[0318] If the actual DTW distance is less than the currently found minimum distance Min_distance, then Min_distance is updated to the actual DTW distance, and the index of the standard pattern sequence Q is recorded as index_Best_match.

[0319] Step 6-4: Return the minimum distance Min_distance and the corresponding sequence index index_Best_match. Index_Best_match points to the standard pattern sequence Q that is most similar to the sequence C to be matched.

[0320] The time complexity of calculating the lower bound distance is O(n), which is much smaller than the time complexity of calculating the DTW distance, O(n 2 ), so before calculating the DTW distance, the lower bound distance can be used to quickly exclude obviously dissimilar sequences, thereby reducing the number of sequences that need to be calculated for DTW and improving the efficiency of similarity search. The classification method of LB-Keogh distance can reduce the time complexity of the original DTW metric from O(n 2 ) is reduced to O(n).

[0321] According to the DTW distance, the final probability of the fault event is obtained, including:

[0322] Assume that the similarity probability between the optimal solution of an event and the basic pattern should be equal to the DTW distance D dtw Inversely proportional, the similarity probability vector of the basic mode is obtained as:

[0323] ;

[0324] In the formula, M Indicates the number of basic modes that may cause changes in the flow data of the wellbore 2; obtains the change law of the probability of a fault event over time through real-time calculation and matching analysis;

[0325] is the similarity probability of each basic pattern; this vector represents the similarity probability between each basic pattern and the sequence to be matched; is the DTW (Dynamic Time Warping) distance between the sequence to be matched and the i-th standard pattern; i is the number of the standard pattern, ranging from 1 to M (indicating that there are M standard patterns); It is the reciprocal of the DTW distance of the iiith standard pattern, which is used to indicate the similarity with the standard pattern; the smaller the distance, the higher the similarity. is the sum of the inverses of the DTW distances of all standard modes and is a normalization factor used to ensure that all sums are 1, thus forming a probability distribution.

[0326] The PWD measured data of SHI153 well are used for analysis. The PWD measured data of SHI153 well are obtained, and the measured data segment curve is selected as follows: Figure 6 shown.

[0327] According to the methods of step 1 to step 3, the PWD measured data is cleaned, Lagrange interpolated, and normalized to obtain the standard data set.

[0328] According to the method in step 4, the standard data set is subjected to Kalman filtering for denoising. The denoised PWD measured data curve is as follows: Figure 7 shown.

[0329] Establishing a multiphase flow model can extract standard modes for various working conditions. The standard mode diagrams for various working conditions are as follows: Figure 8 shown.

[0330] The LB_Keogh and DTW algorithms were established to match the overflow standard mode. The DTW distance with the normal drilling standard mode was 11.13623899867073; the DTW distance with the pump-on standard mode was 4.08359353178429; the DTW distance with the drilling standard mode was 8.65990000651971; the DTW distance with the drilling standard mode was 3.42227282731567; the DTW distance with the overflow standard mode was 2.63613005267348; and the DTW distance with the leakage standard mode was 8.58086544885356. (Keep 14 decimal places) After automatic sorting by the algorithm, it can be seen that the DTW distance between the selected PWD measured data and the overflow standard mode is the smallest, that is, the selected PWD measured data is identified as the overflow mode. Fig. 9 It is a schematic diagram of the PWD measured data curve after preprocessing; Fig.10 Schematic diagram of pattern matching for the overflow standard mode.

[0331] The DTW distance calculation shows that the probability of a fault event is 91.32%.

[0332] Compared with the traditional leakage diagnosis method based on threshold judgment, the leakage pattern recognition method driven by multi-source data-physical model coupling established in the present invention can quickly identify abnormal working conditions and divide the probability of occurrence by integrating bottom hole pressure, flow rate, and mud pool incremental data before the accident occurs and deteriorates, which can effectively shorten the leakage diagnosis time and reduce the leakage false alarm rate, and realize intelligent monitoring of the leakage process. The method is simple to operate, the measurement parameters are easy to obtain, and the real-time output of the probability of leakage can be realized, which provides a reference for field engineers to quickly understand the development process of leakage and accurately complete subsequent well control operations, effectively improving the efficiency and safety of offshore drilling.

Claims

1. A drilling spill pattern recognition method driven by data-model coupling, characterized in that: The drilling leakage pattern recognition system is realized through the data-model coupling driven drilling leakage pattern recognition system, which includes a flow parameter acquisition system and a data intelligent analysis system; The flow parameter acquisition system is used to acquire flow parameter data during the drilling process; The data intelligent analysis system is used to process the flow parameter data during the drilling process, perform noise reduction and pattern recognition, analyze the change in wellbore pressure after the accident well overflows, and intuitively display the probability of overflow occurrence; The flow parameter acquisition system includes: downhole PWD tools, standpipe pressure gauge, casing pressure gauge, inlet flow meter, outlet flow meter, mud pool level meter; The downhole PWD tool is connected to the logging tool and drill pipe at the top and the drill bit at the bottom. During the drilling process, the downhole PWD tool measures the downhole pressure data in real time. The standpipe pressure gauge is installed at the top of the drill pipe to measure the flow pressure change after the drilling fluid is pumped into the standpipe; The casing pressure gauge is installed at the front end of the throttle manifold of the wellhead throttle pipeline to measure the fluid pressure change when the drilling fluid is returned from the wellbore to the platform; The inlet flow meter is installed between the drilling pump and the drill pipe to measure the flow rate of drilling fluid pumped into the drill pipe; The outlet flow meter is connected to the choke manifold and the wellhead spool to monitor the flow rate changes of the fluid returned from the annulus; Mud pool level meter, used to measure the change of mud pool liquid level during drilling; The data intelligent analysis system includes: data acquisition pipeline, computer, display, alarm instrument; The data acquisition pipeline connects the computer and the flow parameter acquisition system to transmit the measurement signal of the flow data; the computer is used to analyze the change pattern of the wellbore pressure after the overflow of the current accident well, and perform noise reduction and pattern recognition on the collected data; the display is connected to the computer to intuitively output the curve of the overflow probability during the drilling process; including: Perform data cleaning, Lagrange interpolation, and normalization on various types of measurement data to generate standard data sets; including: Step 1: Use the box plot method to clean outliers; Step 2: For the missing values ​​of some eigenvalue information in the sample data, Lagrange interpolation method is used to fill them; Step 3: Normalization processing; Construct a Kalman filter to reduce the noise of the obtained standard data set; Establish a wellbore multiphase flow model to extract the standard mode of flow parameters under abnormal conditions of overflow and leakage as well as normal conditions of drilling, pumping, drilling and drilling; The LB_Keogh method based on the constrained search space and lower bound method is combined with the DTW algorithm based on dynamic programming to accelerate the time series similarity search, and the PWD measurement data is pattern matched using the obtained standard patterns of different working conditions. According to the DTW distance, the final probability of the fault event is obtained; Establish a wellbore multiphase flow model to extract the standard mode of flow parameters under abnormal conditions of overflow and leakage as well as normal conditions of drilling, pumping, drilling and drilling; Step 5-1: Establish the mass conservation equation of gas and liquid phase: The momentum conservation equation for multiphase flow in the wellbore is established as: Where: A2 is the cross-sectional area of ​​the annulus, m 2 ; α g is the void fraction of the gas, m 3 / m 3 ; αl is the porosity of drilling fluid, m 3 / m 3 ; p is fluid pressure, MPa; θ is well inclination, rad; K is reservoir permeability, m 2 ;p res is the reservoir pressure, Pa; p wf is the fluid pressure at the bottom of the well, Pa; S is the skin coefficient of the reservoir; t is the overflow time, s; s is the distance, m; γ is the Euler constant, γ = 0.5772; φ is the reservoir porosity; μ is the viscosity of the gas in the reservoir, Pa·s; r w is the wellbore radius, m; ρ g is the gas phase density, kg / m 3 :ρ l is the liquid density, kg / m 3 , v g is the gas phase velocity, m / s; v l is the liquid phase velocity, m / s; Step 5-2: Set initial and boundary conditions, simulate abnormal conditions of overflow and leakage, and normal conditions of drilling, pumping, drilling and drilling, and obtain the variation law of flow parameters; in overflow conditions: Step 5-2-1: Set initial conditions; during normal drilling, when drilling to the top of the reservoir and no formation fluid flows in: E H (S,0)=E g (S,0)=0 E=1-E mc p(S,0)=p(S) Among them, Q m is the displacement of the mud pump, m 3 / s;q c is the discharge volume of cuttings, m 3 / s; S is the well depth, m; E H is the hydrate volume fraction; E g is the gas volume fraction, 0-1; E c is the volume fraction of rock cuttings; V sc is the actual flow velocity of the cuttings, m / s; E m is the volume fraction of drilling fluid, 0-1; v m is the superficial velocity of drilling fluid, m / s; v sm is the actual flow rate of drilling fluid, m / s; ρ c is the density of rock cuttings, kg / m 3 ;v c is the superficial velocity of the cuttings, m / s; Step 5-2-2: Set boundary conditions; the boundaries are: Drilling overflow conditions: Where p is pressure, Pa; p s is the wellhead casing pressure, Pa; q g is the volume flow rate of the invading gas, m 3 / s;q c is the discharge volume of cuttings, m 3 / s;r H is the hydrate phase transition rate, m 3 / s; Set simulation parameters, including drilling fluid density, plastic viscosity, dynamic shear force, reservoir permeability, initial pressure difference, displacement, Step 5-3: Change the simulation parameters, conduct a large number of numerical simulations, normalize the simulation results, and obtain the time series of the bottom hole pressure and mud pool increment under abnormal conditions of overflow and leakage and normal conditions of drilling, pumping, drilling and drilling, that is, the basic change mode; The LB_Keogh method based on the constrained search space and lower bound method is combined with the DTW algorithm based on dynamic programming to accelerate the time series similarity search, and the PWD measurement data is pattern matched using the obtained standard modes of different working conditions; including: Step 6-1: Initialize the minimum distance variable Min_distance: Min_distance is used to store the minimum DTW distance between the sequence to be matched and each standard pattern sequence. For each standard pattern sequence, the operations from step 6-2 to step 6-3 are performed on it in sequence; Step 6-2: Calculate the LB_Keogh lower bound distance LB between a standard pattern sequence Q and the sequence to be matched C; it is known that the lower bound distance LB is always smaller than the DTW distance. If LB>Best_Distance, the current calculation must not be the optimal match; if the lower bound distance is smaller than the currently found minimum distance Min_distance, continue to execute step 6-3, otherwise skip the current standard pattern sequence, select the next standard pattern sequence and repeat step 6-2; Best_Distance is a variable used to track the current optimal distance between the sequence to be matched and the standard pattern sequence; Step 6-2-1: Assume that the time series of length m of the first standard variation pattern is: Q = [q1, q2, ..., q m ]; Input the time series of PWD measurement data with length n to be matched: C=[c1,c2,...,c n ]; Step 6-2-2: Initialize the distance matrix and calculate the local distance Create an n×m matrix D to store the local distance between each pair of data points between the two sequences; Step 6-2-3: Set constraints Think of a global constraint or a local constraint as an index into a constrained warp path D(i,j) such that jr≤i≤j+r, where r is a term that defines the reach or allowed range of the warp at a given point in the sequence; Step 6-2-4: Calculate the upper and lower boundaries and the distance between the lower and upper boundaries According to the constraints, the upper boundary of the sequence Q near the adjacent time window r1 is set to U i , the lower boundary is L i , forming an envelope that surrounds the sequence Q above and below, and the upper and lower boundaries are defined as follows: U i =max(q i-r1 :q i+r1 ) L i =min(q i-r1 :q i+r1 ) Among them, for each time node i, there is L i ≤q i ≤U i ; According to the upper and lower boundaries, set the lower boundary distance LB to: Step 6-2-5: Compare the lower bound distance LB with the minimum distance Min_distance; If the lower bound distance of the standard pattern sequence Q is less than the currently found minimum distance Min_distance, continue to execute step 6-3, otherwise skip the current standard pattern sequence, select the next standard pattern sequence and repeat step 6-2; Step 6-3: Calculate the actual DTW distance between the standard pattern sequence Q and the sequence to be matched C. If the actual DTW distance is less than the currently found minimum distance Min_distance, update Min_distance to the actual DTW distance and record the index of the standard pattern sequence Q as index_Best_match; Step 6-3-1: Calculate DTW distance Dynamic programming is used to calculate the minimum cumulative distance and the cumulative distance matrix A, where A[i,j] represents the minimum cumulative distance between the i-th point in sequence Q and the j-th point in sequence C. The specific calculation method is as follows: The final A(n,m) is the DTW distance D between the two sequences. dtw ; Step 6-3-2: Compare DTW distance with minimum distance Min_distance If the actual DTW distance is less than the currently found minimum distance Min_distance, update Min_distance to the actual DTW distance, and record the index of the standard pattern sequence Q as index_Best_match; Step 6-4: Return the found minimum distance Min_distance and the corresponding sequence index index_Best_match; index_Best_match points to the standard pattern sequence Q that is most similar to the sequence C to be matched.

2. The method for drilling spill pattern recognition driven by data-model coupling according to claim 1, characterized in that: Establish the time series of PWD real-time second data in the drilling process, that is, multi-type measurement data, obtain PWD real-time second data at a certain time interval, and obtain the time series of PWD real-time second data in the drilling event: P = <p1,p2,...,p i ,...,p m >; Wherein, P is the time series of PWD real-time second data; pi is the bottom hole pressure at the i-th time node, and m is the length of the time series; Use the box plot method to clean outliers; including: Step 1-1: Based on the time series P of the PWD real-time second data in the drilling event <p1,p2,...,p i ,...,p m >, calculate the first quartile Q1 and the third quartile Q3, as well as the interquartile range IQR = Q3-Q1; Step 1-2: Define the upper edge Up_Fence and the lower edge Low_Fence: Up_Fence=Q3+k·IQR Low_Fence=Q1-k·IQR Where k is the coefficient used to determine outliers; Step 1-3: Define the data range: According to the upper edge Up_Fence, the lower edge Low_Fence, the first quartile Q1 and the third quartile Q3, the data is divided into three parts: The box range is: Q1≤pi≤Q3; The box whisker range is: Low_Fence≤pi≤Q1 and Q1≤pi≤Up_Fence; Outlier points: data points that are smaller than Low_Fence and larger than Up_Fence; Step 1-4: Delete outlier points.

3. The method for drilling leakage pattern recognition driven by data-model coupling according to claim 1, characterized in that: For the missing values ​​of some eigenvalue information in the time series of sample data, i.e. PWD real-time second data, the Lagrange interpolation method is used to fill in the missing values, including: Step 2-1: The Lagrange polynomial is a cubic polynomial; assuming that the time series P of the PWD change in the drilling event is (1) For: P (1) ={(t1,p1),(t2,p2),...,(t i ,p i ),...,(t m ,p m )}; where (t i ,p i ) is t i The bottom hole pressure at time, m is the length of the time series; Step 2-2: Set the time point at which missing values ​​need to be filled, i.e. the target point, as t i ; Step 2-3: At the target point t i Select 4 uniformly distributed data points near to form a set S: S = {(t0, p0), (t1, p1), (t2, p2), (t3, p3)}; Step 2-4: For each interpolation point (t i ,p i ), define the Lagrangian basis function as L i (t): Where t is the time t, t i and t j is the time of the i-th and j-th data points in the set S; Step 2-5: Construct the interpolation polynomial P(t) using the Lagrange basis function and the four known data points: Among them, L i (t) is the Lagrangian basis function; Normalization processing; including: Assume that the one-dimensional time series after Lagrange interpolation is P (3) =[p (3) 1,p (3) 2,...,p (3) i ,...,p (3) m ]; where p (3) i is each element in the sequence; the data at each time point is scaled to between [0,1] as follows: Where min(P) is the one-dimensional time series P (3) The minimum value in , max(P) is the one-dimensional time series P (3) The maximum value in ; After normalization, we finally get a standard data set: P (4) =[p (4) 1,p (4) 2,...,p (4) m ].

4. The method for drilling spill pattern recognition driven by data-model coupling according to claim 1, characterized in that: Construct a Kalman filter to reduce the noise of the obtained standard data set; including: Step 4-1: Build the model Assume that the one-dimensional time series vector P of the PWD real-time second data after Lagrange interpolation is (2) For: P (2) =[p (2) 1,p (2) 2,...,p (2) m ]; where p (2) m Represents the bottom hole pressure at the mth time node; that is: p k =f(t k )+n k Among them, t k represents the time of the kth sampling point, n k Represents the noise at the kth sampling point. By setting the true value of the bottom hole pressure at a certain moment as the state variable and the bottom hole pressure at a certain moment with noise as the observed value, the following state space model is established: Step 4-1-1: Establish the state equation: x k =Ax k-1 +Bu k +w k In the formula, x k is the bottom hole pressure value at time k; A is the state transfer matrix, which represents the conversion of the bottom hole pressure value at time k-1 to the bottom hole pressure value at time k; B is the influence matrix of the control input; u k is the control input at time k, w k is the prediction noise of bottom hole pressure at time k; Step 4-1-2: Establish the observation equation: y k =Hx k +v k In the formula, y k is the observed value of the bottom hole pressure at time k; H is the observation matrix, v k is the prediction noise of bottom hole pressure at time k; Step 4-2: Kalman filter process Step 4-2-1: Initialization Initialize the system state x0, that is, set the initial state estimate and covariance matrix required by the Kalman filter; take the first numerical point in the normalized PWD real-time second data time series standard data set; State transfer matrix A1: Take the state transfer matrix A1 as scalar 1; State covariance matrix P0: Take the state covariance matrix P0 as scalar 1; Prediction noise covariance matrix Q: Take the prediction noise covariance matrix Q as a scalar 0.05; Covariance matrix R of observation noise: Take the covariance matrix R of observation noise as scalar 5; Step 4-2-2: Prediction step According to the state equation, the state and state covariance of the bottom hole pressure at the next moment are predicted: In the formula, is the estimated value of the predicted bottom hole pressure PWD, is the predicted bottom hole pressure PWD state covariance matrix, Q is the predicted noise covariance matrix, u k is the control input at time k; Step 4-2-3: Update steps According to the observation equation, the predicted value is obtained by combining the measured value of bottom hole pressure PWD Make corrections to get the optimal estimate: In the formula, K k is the Kalman gain, R is the covariance matrix of the observation noise, x k is the optimal state estimate, P k is the optimal state covariance matrix, a matrix with I being 1; Step 4-2-4: Repeat the prediction and update steps The prediction and update are repeated in the order of the one-dimensional time series of the bottom hole pressure PWD, and the state estimation value and the state covariance matrix are continuously updated to obtain the optimal state estimation.

5. The method for drilling leakage pattern recognition driven by data-model coupling according to any one of claims 1 to 4, characterized in that: According to the DTW distance, the final probability of the fault event is obtained, including: Assume that the similarity probability between the optimal solution of an event and the basic pattern should be equal to the DTW distance D dtw Inversely proportional, the similarity probability vector of the basic mode is obtained as: Where M represents the number of basic modes that may cause changes in wellbore flow data. Through real-time calculation and matching analysis, the change law of the probability of fault events over time is obtained. P i is the similarity probability of each basic pattern; this vector represents the similarity probability between each basic pattern and the sequence to be matched; D i dtw is the DTW distance between the sequence to be matched and the i-th standard pattern; i is the number of the standard pattern, ranging from 1 to M; is the reciprocal of the DTW distance of the iiith standard pattern, used to indicate the similarity with the standard pattern; is the sum of the inverses of the DTW distances of all standard modes, and is a normalization factor used to ensure that all sums are 1, thus forming a probability distribution.

Images