An underwell gas invasion and loss parameter inversion method based on unscented Kalman filter

Abstract

The present application relates to a kind of well gas invasion and loss parameter inversion methods based on unscented Kalman filter, belong to petroleum engineering technical field, comprising: step 1: establish overflow wellbore multiphase flow model;Step 2: establish inversion model;Step 3: solve inversion model;Calculate the inversion result of downhole gas invasion and loss parameter.The application method combines wellbore multiphase flow equation with state space model, uses the nonlinear estimation characteristics of UKF algorithm, carries out dynamic filtering and state updating to wellhead measurable data, realizes the real-time inversion of downhole gas volume fraction, gas invasion flow and loss amount.UKF method can accurately capture the nonlinear characteristics of system without linearization, thereby improving the inversion accuracy and calculation stability under complex conditions.Through example verification, this method has good applicability under single reservoir, multiple reservoir and overflow and leakage coexistence conditions, provides reliable boundary conditions and intelligent control basis for well control system.

Description

Technical Field

[0001] This invention belongs to the field of petroleum engineering technology, specifically relating to a method for inverting downhole gas intrusion and leakage parameters based on unscented Kalman filtering. Background Technology

[0002] In oil drilling, wellbore pressure control is a core element in ensuring drilling safety and efficient operation. The pressure within the wellbore must be consistently maintained between formation pressure and fracture pressure to prevent gas intrusion (gas invasion, overflow) or drilling fluid leakage into the formation. With increasing drilling depth and formation complexity, downhole temperature, pressure, and multiphase flow conditions change more drastically, significantly increasing the difficulty of well control. Once the wellbore pressure becomes unstable, it can easily trigger major safety accidents such as blowouts and lost circulation, posing a serious threat to personal safety and oil and gas resource development. During well control, the uncertainty of downhole gas invasion and leakage information makes it impossible to formulate targeted well control and well control plans, which is the fundamental reason for engineering problems such as "the more you control, the more leakage occurs," simultaneous overflow and leakage, and low success rate of first-time well control.

[0003] Traditional methods primarily rely on parameters such as flow rate difference, mud pool increment, and pressure changes, using threshold judgments (e.g., "mud pool increment exceeding 1 cubic meter is considered a well kick") to determine whether leakage has occurred. However, for leakage information crucial to well control decisions and the success or failure of well operations, such as formation depth and leakage rate, the complex downhole environment makes it difficult to directly measure flow parameters. On-site, inferences are usually made based on human experience, leading to large errors, low well control efficiency, and even ultimately evolving into serious well leakage or blowout accidents. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention proposes a method for inverting downhole gas intrusion and leakage parameters based on unscented Kalman filter (UKF).

[0005] This invention combines the multiphase flow equations of the wellbore with a state-space model, utilizing the nonlinear estimation characteristics of the UKF algorithm to dynamically filter and update the state of measurable data at the wellhead (such as mud pit increments and outlet flow rates), achieving real-time inversion of downhole gas volume fraction, gas intrusion flow rate, and leakage. The UKF method can accurately capture the nonlinear characteristics of the system without linearization, thereby improving the inversion accuracy and computational stability under complex operating conditions. Through numerical examples, this method demonstrates good applicability under single-reservoir, multi-reservoir, and coexisting overflow and leakage conditions, providing reliable boundary conditions and intelligent control basis for well control systems.

[0006] The technical solution of this invention is as follows:

[0007] A method for inverting downhole gas intrusion and leakage parameters based on unscented Kalman filtering includes:

[0008] Step 1: Establish a multiphase flow model for the overflow control wellbore;

[0009] The process involves establishing a basic mathematical model and assumptions through mathematical modeling; describing the general laws of two-phase flow through basic governing equations; supplementing and specifying the basic governing equations through auxiliary equations; establishing the unique solution to the problem through boundary conditions on the initial state and boundary conditions; dividing the continuous wellbore space and continuous time into a finite number of discrete grid points that can be processed by a computer through grid partitioning of the boundary domain; using the finite difference method to transform partial differential equations that cannot be directly solved into algebraic equations that can be solved at discrete grid points through the discretization of the governing equations; and using the numerical solution of the equation system to design an iterative algorithm to actually solve these equations, thus realizing a complete modeling and solution process for the overflow well control wellbore multiphase flow model from physics to mathematics to computation.

[0010] Step 2: Establish the inversion model; set the unmeasurable downhole parameters as state variables, and set the mud pit increment and outlet flow rate as observation variables, and establish the state equation and observation equation;

[0011] Step 3: Solve the inversion model; generate Sigma points by performing an unscented transformation on the state variables using the unscented Kalman filter algorithm, substitute the obtained Sigma points into the observation equation, call the overflow control wellbore multiphase flow model to calculate the observations and perform filtering update, obtain and output the inversion results of downhole gas intrusion and leakage parameters based on the filtered updated state variables, including: bottom hole gas intrusion flow rate, multi-reservoir gas intrusion flow rate, leakage flow rate, and formation pressure.

[0012] According to a preferred embodiment of the present invention, a basic mathematical model and basic assumptions are established through mathematical modeling; including:

[0013] The basic mathematical models include drift flow model, homogeneous flow model, and phase-separated flow model; the basic assumptions are as follows:

[0014] ① The flow of fluid in the annulus is considered as a one-dimensional flow along the wellbore direction;

[0015] ②The gas and liquid phases reach equilibrium instantaneously in the wellbore;

[0016] ③Consider the dissolution and precipitation of CO2 and CH4 gases in the drilling fluid, and ignore the chemical reactions between CO2 and CH4 gases and certain components in the drilling fluid;

[0017] ④ Stable radial heat transfer occurs between the wellbore and the formation.

[0018] According to a preferred embodiment of the present invention, the basic governing equations include a continuity equation and a momentum equation;

[0019] Based on the law of conservation of mass, a continuity equation is established along the flow direction in the wellbore annulus, with the coordinates along the flow direction in the wellbore annulus set as follows: z Extracting micro-segments dz The study found that the cross-sectional area of ​​the annular flow was... A Considering the gas precipitation process during its ascent in the wellbore annulus, the continuity equations for the gas-producing zone, the non-producing zone, and the liquid phase are obtained, as shown in equations (I), (II), and (III), respectively:

[0020] (I);

[0021] (II);

[0022] (III);

[0023] In the formula, The density of free gas at the given temperature and pressure within the annulus, in kg·m³. -3 ; Density of drilling fluid, kg·m -3 ; , The upward return velocities of free gas and drilling fluid, respectively, are in m / s. -1 ; , These are the integrals of free gas and drilling fluid, respectively, and are dimensionless. A Let m be the cross-sectional area of ​​the annulus. 2 , ,in D The inner diameter of the well shaft is in meters (m). d The outer diameter of the drill string for different well sections, in meters; The mass of gas produced per unit time per unit thickness of reservoir, in kg·s -1 ·m -1 ; Rs m represents the solubility of the gas in the drilling fluid. 3 ·m -3 ; The density of the gas under standard conditions is expressed in kg·m³. -3 ; This is the volume coefficient of the drilling fluid in the local area, dimensionless.

[0024] According to the law of conservation of momentum: the rate of change of momentum over time within a unit cell is equal to the sum of all external forces acting on that unit cell, resulting in the momentum equation, as shown in equation (IV):

[0025]

[0026] (IV);

[0027] In the formula, P Pressure, Pa; The angle between the wellbore axis and the direction of gravity is fr, which refers to the pressure gradient caused by friction.

[0028] According to a preferred embodiment of the present invention, the auxiliary equations supplement and specify the basic governing equations; including:

[0029] The auxiliary equations include the two-phase flow regime determination equation, the velocity equation, the fluid property equation, and the annular friction equation;

[0030] The current flow state is determined by the two-phase flow state determination equation, and the corresponding velocity equation and annular friction equation are selected accordingly.

[0031] The convection terms in the continuity and momentum equations are supplemented by introducing the distribution coefficient and slip velocity from the drift flow model in the velocity equation.

[0032] The annular friction equation provides friction models under different flow conditions, enabling the momentum equation to accurately reflect the influence of wellbore roughness and flow state on pressure drop.

[0033] According to a preferred embodiment of the present invention, the boundary conditions establish a unique solution to the problem by constraining the initial state and boundary conditions; including:

[0034] The initial state is defined as the distribution of parameters within the wellbore at the start of the simulation;

[0035] Boundary conditions are defined on the boundaries of the wellbore, and the parameters are constraints that change over time.

[0036] According to a preferred embodiment of the present invention, the definite solution domain meshing divides the continuous wellbore space and continuous time into a finite number of discrete mesh points that can be processed by a computer; including:

[0037] Spatial domain grid partitioning includes:

[0038] Mesh generation is performed using a fixed spatial step size method, where the length of any given mesh is... for:

[0039] (V);

[0040] , It refers to the spatial coordinates of two adjacent grid points in the direction of wellbore depth;

[0041] Then the total number of spatial grids for:

[0042] (VI);

[0043] H refers to the total depth of the well;

[0044] Time-domain gridding includes:

[0045] The time step employs a non-uniform grid partitioning method, based on the velocity of the free gas. and the length of the spatial grid at that location Relationship obtained by time step :

[0046] (VII)

[0047] According to a preferred embodiment of the present invention, the discretization of the governing equations utilizes the finite difference method to transform partial differential equations that cannot be directly solved into algebraic equations solvable at discrete grid points; including:

[0048] The partial differential forms of the continuity equation and the momentum equation are unified as equation (VIII):

[0049] (VIII);

[0050] Where X represents the cumulative term of a physical quantity, Y represents the current term of a physical quantity, and z and t are the spatial coordinate and time coordinate, respectively;

[0051] The four-point difference scheme equation of equation (VIII) is shown in equation (IX):

[0052] (IX);

[0053] Where n and j are the time step index and the spatial grid point index, respectively; and These are the spatial step size and the time step size, respectively. This represents the value of a physical quantity at spatial location j+1 and time point n+1. Let the value of the physical quantity be at spatial location j and time point n+1; and The values ​​of physical quantities at spatial points n, j, and j+1 at the same time point; , Let n+1 be the value of the physical quantity at spatial points j and j+1 at the same time point;

[0054] From equation (IX), we get:

[0055] For the continuity equation of the gas-phase producing layer, the difference scheme is shown in equation (X):

[0056]

[0057] (X);

[0058]

[0059] For the continuity equation of the non-productive gas phase, the difference scheme is shown in equation (XI):

[0060]

[0061] (XI);

[0062] For the continuity equation of the liquid phase, the difference scheme is shown in equation (XII):

[0063]

[0064] (XII);

[0065] The difference scheme for the momentum equation is shown in equation (XIII):

[0066]

[0067]

[0068] (XIII).

[0069]

[0070]

[0071] According to a preferred embodiment of the present invention, the numerical solution of the system of equations is used to design an iterative algorithm to actually solve these algebraic equations, including:

[0072] When calculating wellbore pressure, an initial pressure value is first assigned to the pressure node to be calculated. This initial pressure value is then substituted into the equation system to obtain the relevant parameters. The flow regime is determined according to the flow pattern discrimination equation. The pressure of the pressure node to be calculated is then calculated using the flow model. The calculated pressure is compared with the initial pressure value. If the absolute deviation between the calculated pressure and the initial pressure value is less than the specified error range, the same method is used to calculate the next pressure node. Otherwise, a new initial pressure value is assigned and the calculation of this pressure node continues until the constraint conditions are met to calculate the next pressure node.

[0073] According to a preferred embodiment of the present invention, the numerical solution of the system of equations is used to design an iterative algorithm to actually solve these algebraic equations, and the specific implementation steps include:

[0074] 1) Input basic parameters, which include the following:

[0075] Wellbore geometric parameters: wellbore depth H, wellbore inner diameter D, drill string outer diameter d;

[0076] Fluid physical properties: Drilling fluid density Drilling fluid plastic viscosity, drilling fluid dynamic shear stress, drilling fluid apparent viscosity, and gas density. Gas solubility Rs, gas viscosity σ, the surface tension of the gas;

[0077] Other parameters: Drilling fluid displacement Inlet temperature, geothermal gradient, and air intrusion flow rate Spatial step size ∆z, time step size ∆t, error threshold ε;

[0078] 2) Given the initial and boundary conditions, t = t 0 When the value is 0, the pressure at each node is calculated starting from the wellhead. And velocity parameters and fluid properties; velocity parameters include apparent gas phase velocity. Apparent velocity in the liquid phase Gas phase true velocity True velocity of liquid phase Fluid properties include drilling fluid density. Drilling fluid plastic viscosity, dynamic shear force, and gas density Gas viscosity ;

[0079] 3) Based on the flow velocity at the bottom node and the length of the divided cell, substitute the flow velocity and cell length into equation (VII) to obtain the time step. ∆ t0 ;

[0080] 4) Assuming the first time node t 1 = t 0 + ∆ t0 The bottom pressure at the well is Calculate the gas solubility, velocity parameters, and fluid properties at the bottom node;

[0081] 5) Assume the pressure at (1, 1) is Calculate the gas solubility, velocity parameters, and fluid properties at this node;

[0082] 6) Assume the gas phase volume fraction at (1,1) is Based on the continuity equation, calculate the apparent flow velocities of the gas and liquid phases at (1,1). , The gas phase volume fraction at this node is calculated based on the velocity equation. ;

[0083] 7) If satisfied If ≤ε, then it is considered If the assumption is reasonable, otherwise it needs to be re-assumed until the condition is met; ε is the error threshold.

[0084] 8) The obtained Substitute the relevant parameters into the momentum equation to calculate the pressure at (1,1). ;

[0085] 9) If satisfied If ≤ε, then it is considered If the assumption is reasonable, then re-assume. Repeat steps 5)-9) until the condition is met;

[0086] 10) The obtained Using the known pressure as the next node, repeat the above calculation until the wellhead. If the absolute difference between the calculated wellhead pressure and atmospheric pressure meets the accuracy requirement, then it is considered... If the assumption is reasonable, then re-assume. Repeat steps 4)-10) until the conditions are met; among them, calculate the absolute value of the difference between the wellhead pressure and the atmospheric pressure. If the result is less than the error threshold ε, the accuracy condition is met.

[0087] By following the steps above, the parameters of all spatial points within the first time point are calculated. These parameters are then used as known parameters for the next time point, and so on, to calculate... t 2 , t 3 … t n All parameters within a given time point.

[0088] According to a preferred embodiment of the present invention, the inversion model takes the overflow process as the analysis object. In the inversion calculation using the unscented Kalman filter algorithm, the air intrusion and leakage flow rates are taken as state variables, and the mud pit increment and outlet flow rate are taken as observations. The corresponding state equation and observation equation are described as equations (XIV) and (XV):

[0089] (XIV);

[0090] (XV);

[0091] In the formula, k express k time; x k This refers to the current inflow or outflow rate, in m. 3 ·s -1 ; This refers to the current monitoring values ​​of the mud pit's incremental volume and outlet flow rate, in m. 3 ·s -1 ; The measurement noise at the current moment; The system noise at the previous moment; f (·) is the state transition function. x k = x k-1 + w k-1 ; h (·) represents the nonlinear observation equation, which is the wellbore multiphase flow model established in step 1.

[0092] According to a preferred embodiment of the present invention, solving the inversion model includes:

[0093] Step 3.1: The unscented Kalman filter algorithm uses unscented transformation to handle the nonlinear propagation problem of mean and covariance, that is, to approximate the probability density distribution of the nonlinear function and use a series of deterministic samples to approximate the posterior probability density of the state.

[0094] Step 3.2: Determine the initial values;

[0095] Assuming initial state variables x 0 And calculate the corresponding statistical values:

[0096] (XVI);

[0097] In the formula, The average of the initial state variables; The initial state quantity covariance; It refers to the expectation in mathematical notation. It is the initial state vector;

[0098] Step 3.3: Calculate the predicted value;

[0099] Substituting the 2n+1 Sigma points into the state transition function, we obtain the 2n+1 Sigma points in... k Predicted value at time :

[0100] (XVII);

[0101] In the formula, Sigma points refer to a set of sample points that are deterministically selected based on the mean and covariance of the current state variables in the unscented transformation step of the unscented Kalman filter algorithm.

[0102] Step 3.4: Calculate the predicted mean and covariance matrix;

[0103] Calculate the system state variables in k Predicted mean at time Covariance Matrix :

[0104] (XVIII);

[0105] In the formula, This represents the predicted value of the state variable at time k based on the observation data at time k-1; k is the time index; n represents the dimension of the vector, i.e., the number of parameters contained in the state variable; i represents the index of the Sigma point. This represents the weight corresponding to the i-th Sigma point; This represents the value of the i-th Sigma point after propagation through the state transition function; Represents the time based on observation data up to time k-1. The predicted state covariance matrix; This represents the covariance weight corresponding to the i-th Sigma point;

[0106] Step 3.5: Perform another UT transformation;

[0107] Based on the obtained predicted mean and covariance matrix, the UT transformation is applied again to calculate... k The Sigma point ζ of the time-time observation (i) and corresponding weights W m (i) and W c (i) This includes: the predicted mean obtained in step 3.4. Covariance Matrix Substituting these values ​​into equation (XVI) as new mean and covariance inputs, the Sigma point ζ for the observation equation at time k is calculated. (i) Weights of Sigma points W m (i) and W c (i) Determined by equation (XVII);

[0108] Step 3.6: Calculate the new predicted value;

[0109] Substituting the new 2n+1 Sigma points into the observation equation, we obtain k Predictive observations for time:

[0110] (XIX);

[0111] Step 3.7: Calculate the predicted mean, variance, and covariance matrix;

[0112] Calculate system observations in k Predicted mean at time ,variance Covariance Matrix :

[0113] (XX);

[0114] In the formula, This represents the predicted mean of the observations at time k based on the observation data at time k-1; This represents the weight corresponding to the i-th Sigma point; This represents the value of the i-th Sigma point after propagation through the nonlinear observation equation; This represents the observation-prediction covariance matrix at time k; This represents the weight corresponding to the i-th Sigma point; This represents the covariance weight corresponding to the i-th Sigma point; Represents the observation noise covariance matrix; Indicates at time The mean of the state prediction;

[0115] Step 3.8: Calculate the Kalman gain matrix;

[0116] (XXI);

[0117] In the formula, This represents the Kalman gain matrix at time k; This represents the observation-prediction covariance matrix at time k; This represents the cross-covariance matrix between the state variables and the observations at time k;

[0118] Step 3.9: Update the system's state variables and covariance:

[0119] (XXII);

[0120] In the formula, That is k Constant air intrusion or leakage flow , .

[0121] Further preferably, the specific implementation process of step 3.1 includes:

[0122] Unscented transform of proportionally corrected symmetric sampling, assuming a nonlinear transform y = f ( x), state vector x Given an n-dimensional random variable and a known state vector. x mean The variance P is used to obtain 2n+1 Sigma points and corresponding weights through an unscented transformation. W ;

[0123] The formula for determining the Sigma point is shown in equation (XXIII):

[0124] (XXIII);

[0125] weight W The formula for obtaining the result is shown in equation (XXIV):

[0126] (XXIV);

[0127] In the formula, (i) Let be the i-th Sigma point; n is the state mean; n is the state dimension; This is a scaling parameter used to reduce the overall prediction error. = (n+κ)-n, κ As a scaling factor, when x When it is a single variable κ =0, when x When there are multiple variables κ =3-n; It refers to the first of the 2n+1 Sigma points generated in the Unscented Transform (UT), i.e., the center point; For matrix The i-th column of the square root matrix obtained by Cholesky decomposition; P is the covariance matrix of the current state; W m (i) , W c (i) These are the weights used for weighting the mean and covariance, respectively. α The scaling factor is 0 ≤ α ≤1, β Higher-order characteristics that reflect historical information of state variables;

[0128] Calculate using equations (XXIII) and (XXIV) k The sigma point of the state quantity at time -1 and the corresponding weights and .

[0129] The beneficial effects of this invention are as follows:

[0130] 1. Analysis of the single reservoir gas intrusion calculation case at the bottom of the well shows that the method can accurately reflect the real change trend of gas intrusion flow rate under both constant wellhead discharge rate and variable wellhead discharge rate conditions. The error is relatively large only in the initial stage of inversion and when the discharge rate is changed, while the overall error after the algorithm is stabilized is within 3%.

[0131] 2. When the key parameters of the reservoir, such as permeability and porosity, are known, this method can directly invert the formation pressure, and the relative error after stabilization is within 0.2%.

[0132] 3. For simultaneous gas intrusion into multiple reservoirs, this method can also track the upward trend of gas intrusion flow in each reservoir. The UKF estimation error for the bottom-hole reservoir (first reservoir) stabilized at less than 2%, but the calculated values ​​were generally smaller than the actual values. The gas intrusion flow inversion accuracy for the second reservoir was relatively low, deviating significantly from the actual values ​​at some time points, but overall remained within 4%.

[0133] 4. For the coexistence of multiple reservoir overflows and leaks, this method can also reflect the process of gas intrusion flow rising and leakage flow rising and then falling well. The overall error after stabilization is within 4%.

[0134] 5. Considering the decrease in the accuracy of the UKF algorithm for inverting gas intrusion and leakage flow in the second reservoir after considering multiple reservoirs, this can be improved by adding additional measured data constraints. Attached Figure Description

[0135] Figure 1 A physical model diagram of the mass conservation equation;

[0136] Figure 2 This is a schematic diagram of a four-point difference grid;

[0137] Figure 3 This is a flowchart of a simulation calculation program;

[0138] Figure 4 Flowchart for inverting and solving downhole gas intrusion / leakage flow rates;

[0139] Figure 5 This is a diagram simulating the wellbore structure and drill string assembly.

[0140] Figure 6 A graph showing the actual value of bottom hole gas intrusion and the estimated value of UKF (under constant wellhead discharge conditions);

[0141] Figure 7 A graph showing the error in the bottom gas intrusion flow inversion results (under constant wellhead discharge conditions);

[0142] Figure 8 The graph shows the results of the observation calculations (constant wellhead discharge conditions);

[0143] Figure 9 A graph showing the actual value of bottom hole gas intrusion and the estimated value of UKF (under varying wellhead discharge conditions);

[0144] Figure 10 A graph showing the error in the bottom gas intrusion flow inversion results (with varying wellhead discharge conditions);

[0145] Figure 11 A graph showing the results of observational calculations (variable wellhead discharge conditions);

[0146] Figure 12 A graph showing the actual formation pressure and the estimated UKF value;

[0147] Figure 13 Error diagram of formation pressure inversion results;

[0148] Figure 14 A graph showing the actual values ​​of gas intrusion flow rates in multiple reservoirs and estimated values ​​of UKF.

[0149] Figure 15 Error diagram of multi-reservoir gas intrusion flow inversion results;

[0150] Figure 16 This is a graph showing the results of the observation calculations;

[0151] Figure 17 A graph showing the actual values ​​of gas intrusion / leakage flow rates and estimated UKF values ​​for multiple reservoir overflows and leaks;

[0152] Figure 18 Error diagram of gas intrusion / leakage flow rate inversion results for multiple reservoir overflows and leaks. Detailed Implementation

[0153] The present invention will be further defined below with reference to the accompanying drawings and embodiments, but is not limited thereto.

[0154] Example 1

[0155] A method for inverting downhole gas intrusion and leakage parameters based on unscented Kalman filtering includes:

[0156] Step 1: Establish a multiphase flow model for the overflow control wellbore;

[0157] The process involves establishing a basic mathematical model and assumptions through mathematical modeling; describing the general laws of two-phase flow through basic governing equations; supplementing and specifying the basic governing equations through auxiliary equations; establishing the unique solution to the problem through boundary conditions on the initial state and boundary conditions; dividing the continuous wellbore space and continuous time into a finite number of discrete grid points that can be processed by a computer through grid partitioning of the boundary domain; using the finite difference method to transform partial differential equations that cannot be directly solved into algebraic equations that can be solved at discrete grid points through the discretization of the governing equations; and using the numerical solution of the equation system to design an iterative algorithm to actually solve these equations, thus realizing a complete modeling and solution process for the overflow well control wellbore multiphase flow model from physics to mathematics to computation.

[0158] Step 2: Establish an inversion model; set downhole unmeasurable parameters such as bottom hole gas inflow, gas inflow in each reservoir, leakage flow or formation pressure as state variables, set mud pit increment and outlet flow as observation variables, and establish state equations and observation equations;

[0159] Step 3: Solve the inversion model; use the Unscented Kalman Filter (UKF) algorithm to perform Unscented Transform (UT) on the state variables to generate Sigma points, substitute the obtained Sigma points into the observation equation h(·), call the overflow control wellbore multiphase flow model to calculate the observations and perform filtering update, obtain and output the inversion results of downhole gas invasion and leakage parameters based on the filtered updated state variables, and present the downhole parameters that cannot be directly measured (time series) in the form of curves, including: bottom hole gas invasion flow rate, multi-reservoir gas invasion flow rate, leakage flow rate, and formation pressure.

[0160] Example 2

[0161] The difference between the downhole gas intrusion and leakage parameter inversion method based on unscented Kalman filtering described in Example 1 and the method described in Example 1 is as follows:

[0162] Establish basic mathematical models and fundamental assumptions through mathematical modeling; including:

[0163] Basic mathematical models include drift flow models, homogeneous flow models, and phase-separated flow models. Homogeneous flow models treat the gas and liquid phases as a uniform mixture of two fluids, assuming no relative velocity (i.e., no slippage effect) and identical flow velocities and average physical properties. They are simple and computationally inexpensive, but neglect interphase interactions and slippage effects, resulting in lower accuracy when there are significant velocity differences between the gas and liquid phases or complex flow patterns. Phase-separated flow models treat the gas and liquid phases as separate flows, each with independent velocity, density, and volume fraction, and consider interphase slippage effects. They are more complex, requiring separate conservation equations for both gas and liquid phases and coupling through interphase forces. Drift flow models are essentially improved phase-separated flow models, introducing drift velocity to characterize the relative motion of the gas with respect to the mixture. This model describes the slippage effect of the gas and liquid phases through drift flux relationships, considering the effects of buoyancy, turbulence, and flow pattern changes. It can more reasonably simulate the slippage effect of the gas and liquid phases, especially in vertical or inclined wellbores where it offers higher accuracy. Compared to other methods that can more reasonably characterize the gas-liquid phase slippage effect, the drift flow model is used as the basis for constructing the multiphase flow model of the wellbore for overflow control. The basic assumptions are as follows:

[0164] ① The flow of fluid in the annulus is considered as a one-dimensional flow along the wellbore direction;

[0165] ②The gas and liquid phases reach equilibrium instantaneously in the wellbore;

[0166] ③Consider the dissolution and precipitation of CO2 and CH4 gases in the drilling fluid, and ignore the chemical reactions between CO2 and CH4 gases and certain components in the drilling fluid;

[0167] ④ Stable radial heat transfer occurs between the wellbore and the formation.

[0168] The fundamental governing equations include the continuity equation and the momentum equation;

[0169] Based on the law of conservation of mass, a continuity equation is established along the flow direction in the wellbore annulus, with the coordinates along the flow direction in the wellbore annulus set as follows: z Extracting micro-segments dz The study found that the cross-sectional area of ​​the annular flow was... A The physical model of the mass conservation equation is as follows: Figure 1 As shown. Considering the gas precipitation process during its ascent in the wellbore annulus, the continuity equations for the gas-producing zone, the gas-non-producing zone, and the liquid phase are obtained, as shown in equations (I), (II), and (III), respectively:

[0170] (I);

[0171] (II);

[0172] (III);

[0173] In the formula, The density of free gas at the temperature and pressure within the annulus, in kg·m³. -3 ; Density of drilling fluid, kg·m -3 ; , The upward return velocities of free gas and drilling fluid, respectively, are in m / s. -1 ; , These are the integrals of free gas and drilling fluid, respectively, and are dimensionless. A Let m be the cross-sectional area of ​​the annulus. 2 , ,in D The inner diameter of the well shaft is in meters (m). d The outer diameter of the drill string for different well sections, in meters; The mass of gas produced per unit time per unit thickness of reservoir, in kg·s -1 ·m -1 ; Rs m represents the solubility of the gas in the drilling fluid. 3 ·m -3 ; The density of the gas under standard conditions is expressed in kg·m³. -3 ; This is the volume coefficient of the drilling fluid in the local area, dimensionless.

[0174] According to the law of conservation of momentum: the rate of change of momentum over time within a unit cell is equal to the sum of all external forces acting on that unit cell, resulting in the momentum equation, as shown in equation (IV):

[0175]

[0176] (IV);

[0177] In the formula, P Pressure, Pa; The angle between the wellbore axis and the direction of gravity is... Figure 1 The label indicates that fr stands for friction, which refers to the pressure gradient caused by friction.

[0178] Auxiliary equations supplement and specify the basic governing equations; including:

[0179] The auxiliary equations include the two-phase flow regime determination equation, the velocity equation, the fluid property equation, and the annular friction equation;

[0180] Two-phase flow regime determination equations: The four most widely used flow regimes in the drilling industry are bubbly flow, slug flow, turbulent flow, and annular flow. This invention applies the flow regime determination equations of Hasan and Kaber as follows.

[0181] (1) Bubble flow:

[0182] ;

[0183] ;

[0184] In the formula, u sg , u sl The apparent flow rates for the gas and liquid phases are respectively, in m / s. -1 ; The limiting ascent velocity of a single bubble, in m·s -1 ; σ The surface tension of the gas and liquid is N·m. -1 .

[0185] (2) Slug flow:

[0186] ;

[0187] (3) Agitated flow;

[0188] Generally, stirred flow is regarded as slug flow, and the flow characteristic parameters are calculated according to the flow law relationship of slug flow.

[0189] ;

[0190] (4) Circular flow:

[0191] ;

[0192] The velocity equation is the drift flow model.

[0193] Fluid property equations include gas property equations and liquid property equations, as detailed below:

[0194] (1) Equations for the physical properties of gases;

[0195] ① Natural gas compressibility factor;

[0196] The DPR method was used to calculate the natural gas deviation coefficient. This method was derived by Dranchuk, Purvis, and Robinson by fitting 1500 data points from the Standing-Katz chart using a modified Benedict-Webb-Rubin equation of state. Its expression is shown below:

[0197] ;

[0198] In the formula, Z The natural gas compressibility factor is dimensionless. P pr For the purpose of comparing pressure, the quantity is dimensionless; T pr The temperature is dimensionless for comparison purposes. ρ pr To approximate the dimensionless density, it can be solved using the following formula and Newton's iteration method (assuming...). Z The initial value is 1).

[0199] ;

[0200] In the formula, , , , , , , , , , , , , .

[0201] ② Gas density;

[0202] From the gas law, we get:

[0203] ;

[0204] Under standard conditions:

[0205] ;

[0206] In the formula, P sc Pressure under standard conditions, MPa; m for n Mass of natural gas per mol, in kg; ρ gsc The density of natural gas under standard conditions is expressed in kg·m³. -3 ; Z sc The compressibility factor of natural gas under standard conditions is dimensionless. T sc The temperature under standard conditions is expressed in Kelvin (K).

[0207] Combining the two equations, we get:

[0208] ;

[0209] ③ Gas viscosity

[0210] This invention uses the natural gas viscosity relationship obtained by Carr et al.

[0211] ;

[0212] ;

[0213] In the formula, A 0 -A 15 Viscosity coefficient, dimensionless; μ g The viscosity of natural gas is given in mPa·s. μ ’ g This is the gas viscosity correction value, in mPa·s; T r The temperature is dimensionless for comparison purposes. P r For the purpose of comparing pressure, the quantity is dimensionless; T Temperature, °C.

[0214] ④ Gas solubility;

[0215] If the volumes of base oil, brine, and emulsifier in the drilling fluid are known, the solubility of gas in the drilling fluid can be estimated as the sum of the gas solubility of each liquid phase multiplied by the volume fraction of each liquid phase in the total volume of the drilling fluid, as shown in the following formula.

[0216] ;

[0217] The solubility of a gas in oil, emulsifier, and brine phases can be expressed by the following formula:

[0218] ;

[0219] In the formula, Rs o , Rs br , Rs e The m represents the solubility of the gas in base oil, brine, and emulsifier, respectively. 3 ·m -3 ; φ o , φ br , φ e These represent the volume fractions of base oil, brine, and emulsifier in the drilling fluid, respectively. Rs h , Rs CO2 The m represents the solubility of hydrocarbon gases and CO2 gases in base oil, brine, and emulsifier, respectively.3 ·m -3 ; φ h , φ CO2 These represent the volume fractions of hydrocarbon gases and CO2 gases in the base oil, brine, and emulsifier, respectively.

[0220] O'Bryan et al. provided an empirical equation for the solubility of gases in oils and emulsifiers, as shown below:

[0221] ;

[0222] In the formula, P The pressure at phase equilibrium, expressed in psi; T The temperature at phase equilibrium is °F; the values ​​of a, b, and n are shown in Table 1.

[0223] In addition, the solubility of the gas in the salt water phase was determined by O'Bryan et al. based on the solubility variation law and by curve fitting of experimental data, as shown in equations (XX) and (XXI).

[0224] CO2 gas:

[0225] ;

[0226] ;

[0227] ;

[0228] ;

[0229] ;

[0230] ;

[0231] Hydrocarbon gases:

[0232] ;

[0233] ;

[0234] ;

[0235] ;

[0236] ;

[0237] In the formula, B s This is the salinity correction factor; s The solid content is expressed as %.

[0238] Table 1. Relevant Parameters;

[0239]

[0240] In Table 1, γ g This represents the specific gravity of the gas.

[0241] (2) Liquid property equations

[0242] ① Drilling fluid rheological parameter prediction model:

[0243] ;

[0244] In the formula, P Pressure, Pa; T Temperature, °C; P 0 represents the ground pressure, in Pa; T 0 represents the ground temperature, in °C. f ( P , T ) respectively represent the temperature at which T Pressure is P Performance viscosity under certain conditions ( μ a ) P, T mPa·s, plastic viscosity ( μ p ) P, T mPa·s, dynamic shear force ( τ 0) P, T Pa; f ( P 0, T 0) respectively represent the temperature at which T 0. Pressure is P Performance viscosity under 0 conditions ( μ a ) P0, T0 mPa·s, plastic viscosity ( μ p ) P0, T0 mPa·s, dynamic shear force ( τ 0) P0, T0 Pa; A , B , C , D , E These are parameters related to the properties of drilling fluid, which depend on the formulation of the drilling fluid and can be obtained from rheological experiments of the drilling fluid.

[0245] ② Gas-liquid surface tension:

[0246] ;

[0247] ;

[0248] ;

[0249] In the formula, σ (T) For temperature T Gas-liquid surface tension, N·m -1 ; T Temperature, K; P Pressure, MPa.

[0250] The equation for annular friction is as follows:

[0251] (1) Bubble flow:

[0252] ;

[0253] In the formula, f The coefficient of friction is dimensionless. ρ m For the mixed density, kg·m -3 .

[0254] The coefficient of friction can be determined by equation (XXV):

[0255] ;

[0256] ;

[0257] ;

[0258] In the formula, Δ is the effective roughness of the annulus, in meters (m). Re m The Reynolds number for the mixed fluid is dimensionless. μ m The viscosity is the mixed viscosity, Pa·s.

[0259] The mixture density and mixture viscosity can be calculated using equation (XXVI):

[0260] ;

[0261] ;

[0262] (2) Slug flow:

[0263] ;

[0264] Among them, the coefficient of friction f The calculation method is the same as for bubbly flow, but... Re m as follows:

[0265] ;

[0266] (3) Circular flow:

[0267] ;

[0268] ;

[0269] ;

[0270] In the basic governing equations, the flow is in a uniform or simplified form. However, in actual wellbore, the gas-liquid two-phase flow will exhibit different flow regimes, such as bubbly flow and slug flow. Therefore, the two-phase flow regime judgment equation is used to determine the current flow regime, thereby selecting the corresponding velocity equation and annular friction equation. Since each flow regime corresponds to a different velocity equation and annular friction equation, here we mean first using the two-phase flow regime judgment equation to determine which of the four flow regimes the current flow regime is: bubbly flow, slug flow, stirred flow, or annular flow, and then finding the corresponding velocity equation and annular friction equation.

[0271] In addition, the velocity term in the basic governing equations is abstract. The distribution coefficient and slip velocity introduced from the drift flow model in the velocity equations are used to supplement the convection terms in the continuity and momentum equations. This allows the gas-liquid phase velocities to characterize the slip effect, thus enabling more accurate calculation of volume fraction and pressure gradient. Both the continuity and momentum equations in the basic governing equations contain convection terms centered on gas and liquid phase velocities. To avoid the assumption of homogeneous flow with equal velocities in both gas and liquid phases, this invention introduces a drift flow model, using the gas phase distribution coefficient... C 0 and gas phase slip velocity u gr The actual velocities of the gas and liquid phases are expressed, thus explicitly incorporating the velocity difference between the two phases into the convection term calculation. The fluid property equations concretize the density, viscosity, and other physical property parameters in the fundamental governing equations into measurable physical quantities, enhancing the model's realism; here, density and viscosity are obtained by substituting other measurable physical quantities into the fluid property equations.

[0272] The annular friction equation provides friction models under different flow conditions, enabling the momentum equation to accurately reflect the influence of wellbore roughness and flow state on pressure drop.

[0273] Boundary conditions establish a unique solution to the problem by imposing constraints on the initial state and boundary conditions; these include:

[0274] The initial state is defined as the distribution of parameters within the wellbore at the start of the simulation (t=0). These values ​​are set based on the actual operating conditions of the drilling process, such as normal drilling, pump shutdown and well shut-in, or well control. In the numerical simulation, the initial conditions assign initial values ​​to all spatial grid points as the starting point for iterative calculations.

[0275] (1) Normal drilling conditions;

[0276] During normal drilling, when the gas reservoir is just reached and the gas has not yet entered the wellbore, the flow in the wellbore is considered a single-phase flow. t At time =0, we have:

[0277] ;

[0278] In the formula, P atm Atmospheric pressure, MPa; Q mud For drilling fluid displacement, m 3 ·s -1 ; Density of drilling fluid, kg·m -3 ; u sg , u sl Apparent flow rates in m / s for the gas phase and liquid phase, respectively. -1 ; u g , u l The actual velocities in the gas and liquid phases are shown in m·s, respectively. -1 ; From the wellhead z Total frictional pressure drop at depth of well, MPa.

[0279] (2) Pump shutdown and well closure conditions;

[0280] The initial conditions for shutting down the well during pump shutdown are determined by the overflow situation at the time of well shutdown.

[0281] (3) Well control conditions;

[0282] according to u sl The velocity of each phase is determined by parameters such as =0 and slip velocity, while the volume fraction and wellbore pressure can be determined by the shut-in state.

[0283] Boundary conditions are defined on the boundaries of the wellbore, constraining parameters that change over time. These conditions reflect the physical limitations of actual drilling processes, such as wellhead pressure (atmospheric pressure), bottom hole gas inflow, or formation pressure. Boundary conditions ensure that the model's behavior at spatial boundaries conforms to physical reality.

[0284] (1) Overflow condition;

[0285] At any given moment during the gas invasion process:

[0286] In the formula, H The depth at the bottom of the well is in meters (m). qg The mass of gas produced per unit time per unit thickness of reservoir, in kg·s -1 ·m -1 ; For pressure.

[0287] (2) Pump shutdown and well closure conditions;

[0288] ;

[0289] (3) Well control conditions;

[0290] ;

[0291] In the formula, P e Δ is the formation pressure, Pa; P Additional pressure, Pa, is applied to kill the well.

[0292] The unique solution is established through the combined constraints of the basic governing equations, initial conditions, boundary conditions, and auxiliary equations. The continuity equation and momentum equation describe the fundamental physical laws of gas-liquid two-phase flow. The initial conditions define unique initial values ​​for state parameters such as pressure, volume fraction, and velocity within the wellbore at the start of the calculation. The boundary conditions constrain the physical state of the system at the spatial boundary. The flow regime discrimination equation, drift flow model, and fluid property equations specify the gas-liquid two-phase velocities, slippage effects, and physical parameters. Finally, through finite difference discretization and iterative convergence calculation, a stable solution is obtained under a given error threshold, guaranteeing the uniqueness of the solution in both mathematical and physical terms.

[0293] Boundary domain mesh generation divides the continuous wellbore space and continuous time into a finite number of discrete grid points that a computer can process; including:

[0294] This invention employs the finite difference method to discretize and iteratively solve the above-mentioned wellbore multiphase flow model, using a four-point difference scheme, the specific form of which is as follows: Figure 2 As shown, Figure 2 In the diagram, "●" represents nodes with known parameters, and "▲" represents nodes that are to be found.

[0295] Spatial domain grid partitioning includes:

[0296] Mesh generation is performed using a fixed spatial step size method, where the length of any given mesh is... for:

[0297] (V);

[0298] , It refers to the spatial coordinates of two adjacent grid points in the direction of wellbore depth;

[0299] Then the total number of spatial grids for:

[0300] (VI);

[0301] H refers to the total depth of the well; the total depth divided by the step size and rounded down gives the total number of spatial grids.

[0302] Time-domain gridding includes:

[0303] The time step employs a non-uniform grid partitioning method to track the multiphase flow front, based on the velocity of the free gas. and the length of the spatial grid at that location Relationship obtained by time step :

[0304] (VII)

[0305] Discretization of the governing equations utilizes the finite difference method to transform partially differential equations that cannot be directly solved into algebraic equations solvable at discrete grid points; including:

[0306] The partial differential forms of the continuity equation and the momentum equation are unified as equation (VIII):

[0307] (VIII);

[0308] Where X represents the cumulative term of a physical quantity, Y represents the current term of a physical quantity, and z and t are the spatial coordinate and time coordinate, respectively;

[0309] The four-point difference scheme equation of equation (VIII) is shown in equation (IX):

[0310] (IX);

[0311] Where n and j are the time step index and the spatial grid point index, respectively; and These are the spatial step size and the time step size, respectively. Let be the value of the physical quantity (circulation term) at spatial location j+1 and time point n+1; Let be the value of the physical quantity (circulation term) at spatial location j and time point n+1; and The values ​​of physical quantities (cumulative terms) at spatial points n, j, and j+1 at the same time point; , The values ​​of physical quantities (cumulative terms) at spatial points n+1, j, and j+1 at the same time point;

[0312] From equation (IX), we get:

[0313] For the continuity equation of the gas-phase producing layer, the difference scheme is shown in equation (X):

[0314]

[0315] (X);

[0316]

[0317] For the continuity equation of the non-productive gas phase, the difference scheme is shown in equation (XI):

[0318]

[0319] (XI);

[0320] For the continuity equation of the liquid phase, the difference scheme is shown in equation (XII):

[0321]

[0322] (XII);

[0323] The difference scheme for the momentum equation is shown in equation (XIII):

[0324]

[0325]

[0326] (XIII).

[0327]

[0328]

[0329] Numerical solutions to systems of equations are used to design an iterative algorithm to actually solve these algebraic equations, including:

[0330] Gas intrusion can affect wellbore pressure. In calculating wellbore pressure, this invention first assigns an initial pressure value to the node to be pressured, substitutes this initial pressure value into the equations to obtain relevant parameters, determines the flow regime based on the flow pattern discrimination equation, and uses a flow model to calculate the pressure at the node to be pressured. The calculated pressure is then compared with the initial pressure value. If the absolute deviation between the calculated pressure and the initial pressure value is less than a specified error range, the same method is used to calculate the next pressure node; otherwise, a new initial pressure value is assigned, and the calculation of this pressure node continues until the constraint conditions are met to calculate the next pressure node.

[0331] The numerical solution of the system of equations is used to design an iterative algorithm to actually solve these algebraic equations. The specific implementation steps include:

[0332] 1) Input basic parameters, which include the following:

[0333] Wellbore geometric parameters: wellbore depth H, wellbore inner diameter D, drill string outer diameter d;

[0334] Fluid physical properties: Drilling fluid density Drilling fluid plastic viscosity, drilling fluid dynamic shear stress, drilling fluid apparent viscosity, and gas density. Gas solubility Rs, gas viscosity σ, the surface tension of the gas;

[0335] Other parameters: Drilling fluid displacement Inlet temperature, geothermal gradient, and air intrusion flow rate Spatial step size ∆z, time step size ∆t, error threshold ε;

[0336] 2) Given the initial and boundary conditions, t = t 0 When the value is 0, the pressure at each node is calculated starting from the wellhead. And velocity parameters and fluid properties; velocity parameters include apparent gas phase velocity. Apparent velocity in the liquid phase Gas phase true velocity True velocity of liquid phase Fluid properties include drilling fluid density. Drilling fluid plastic viscosity, dynamic shear force, and gas density Gas viscosity ;

[0337] 3) Based on the flow velocity at the bottom node and the length of the divided cell, substitute the flow velocity and cell length into equation (VII) to obtain the time step. ∆ t0 ;

[0338] 4) Assuming the first time node t 1 = t 0 + ∆ t0 The bottom pressure at the well is Calculate the gas solubility, velocity parameters, and fluid properties at the bottom node; including:

[0339] ①Use the assumed bottom hole pressure as a starting point for iterative calculations, and use the known bottom hole temperature to calculate the fluid properties, including gas solubility, natural gas compressibility factor, gas density, gas viscosity, drilling fluid density, drilling fluid plastic viscosity, dynamic shear force, and gas-liquid surface tension;

[0340] ② Calculate the apparent velocity of the gas phase using the gas intrusion flow rate and the annular cross-sectional area, and calculate the apparent velocity of the liquid phase using the drilling fluid discharge rate and the annular cross-sectional area. Assume a gas phase volume fraction, and then substitute it into the drift flow model to calculate the apparent velocities of the gas and liquid phases. Then update the gas phase volume fraction using the calculated velocities until the gas phase volume fraction calculated twice consecutively is less than the preset threshold.

[0341] ③Use the apparent velocities of the gas and liquid phases obtained in ② and the flow pattern discrimination formula in the auxiliary equation to obtain the actual velocities of the gas and liquid phases;

[0342] ④ Substitute all the parameters calculated in the first three steps into the momentum equation, and compare the calculated result with the initially assumed pressure to see if the assumption is correct. If the absolute value of the difference between the two is less than the threshold, the pressure converges, which means that the equation has been successfully found. t 1 Once the pressure and various parameters at the bottom node are known, calculations for the next spatial node can begin. If the pressure does not converge, a new bottom-hole pressure estimate must be assumed, and the entire loop must be restarted until the pressure converges.

[0343] 5) Assume the pressure at (1,1) (at spatial index 1 and time point 1) is Calculate the gas solubility, velocity parameters, and fluid properties at this node;

[0344] 6) Assume the gas phase volume fraction at (1,1) is Based on the continuity equation, calculate the apparent flow velocities of the gas and liquid phases at (1,1). , The gas phase volume fraction at this node is calculated based on the velocity equation. ;

[0345] (1) Velocity equation – drift flow model;

[0346] ;

[0347] In the formula, u m The apparent flow rate of the mixed phase is in m·s. -1 , u m = u sl + u sg ;C 0 represents the gas phase distribution coefficient, which is dimensionless; u gr Let be the gas phase slip velocity, in m·s -1 .

[0348] (2) Gas phase distribution coefficient ;

[0349] ;

[0350] ;

[0351] ;

[0352] ;

[0353] (3) Gas phase slip velocity:

[0354] ;

[0355] ;

[0356] ;

[0357] 7) If satisfied If ≤ε, then it is considered If the assumption is reasonable, it needs to be re-assumed until the condition is met; ε is the error threshold; and represents the basic parameters to be input, which are set to a value of in this invention. .

[0358] 8) The obtained Substitute the relevant parameters into the momentum equation to calculate the pressure at (1,1). ;

[0359] 9) If satisfied If ≤ε, then it is considered If the assumption is reasonable, then re-assume. Repeat steps 5)-9) until the condition is met;

[0360] 10) The obtained Using the known pressure as the next node, repeat the above calculation until the wellhead. If the absolute difference between the calculated wellhead pressure and atmospheric pressure meets the accuracy requirement, then it is considered... If the assumption is reasonable, then re-assume. Repeat steps 4)-10) until the conditions are met; among them, calculate the absolute value of the difference between the wellhead pressure and the atmospheric pressure. If the result is less than the error threshold ε, the accuracy condition is met.

[0361] By following the steps above, the parameters of all spatial points within the first time point are calculated. These parameters are then used as known parameters for the next time point, and so on, to calculate... t 2 , t 3 … t n All parameters within a given time point.

[0362] This invention uses Visual Basic.NET for programming calculations, and the flowchart of the simulation calculation program is as follows. Figure 3 As shown.

[0363] The inversion model takes the overflow process as the analysis object. In the inversion calculation using the Unscented Kalman Filter (UKF) algorithm, the air intrusion and leakage flow rates are taken as state variables, and the mud pit increment and outlet flow rate are taken as observations. The corresponding state equations and observation equations are described as Equations (XIV) and (XV):

[0364] (XIV);

[0365] (XV);

[0366] In the formula, k express k time; x k This refers to the current inflow or outflow rate, in m. 3 ·s -1 ; This refers to the current monitoring values ​​of the mud pit's incremental volume and outlet flow rate, in m. 3 ·s -1 ; The measurement noise at the current moment; The system noise at the previous moment; f (·) is the state transition function, and the state transition function is: x k = x k-1 + w k-1 Since the inflow and outflow rates are not directly related to time, and the current values ​​are mainly due to changes in pressure difference, it is assumed that... x k = x k-1 + w k-1 ; h (·) represents the nonlinear observation equation, which is the multiphase flow model of the wellbore established in step 1. The mud pit increment and outlet flow rate are obtained by solving the numerical solution method in step 1. This relates to the solution process of the multiphase flow model for the entire wellbore, with the core being the gas intrusion flow rate, which is used as a state variable in UKF. Alternatively, the leakage flow can be used as a source term and substituted into the mathematical model established in step 1, which consists of the basic control equations, auxiliary equations, and boundary conditions. The model can then be solved using a numerical iterative algorithm to finally calculate the predicted value of the corresponding observation.

[0367] Solving the inversion model; including:

[0368] Step 3.1: The Unscented Kalman Filter (UKF) algorithm does not use Taylor expansion to linearize the nonlinear system. Instead, it uses Unscented Transform (UT) to handle the nonlinear propagation problem of mean and covariance. That is, it approximates the probability density distribution of the nonlinear function and uses a series of deterministic samples to approximate the posterior probability density of the state.

[0369] Step 3.2: Determine the initial values;

[0370] Assuming initial state variables x 0 And calculate the corresponding statistical values:

[0371] (XVI);

[0372] In the formula, The average of the initial state variables; The initial state quantity covariance; It refers to the expectation in mathematical notation. It is the initial state vector;

[0373] Step 3.3: Calculate the predicted value;

[0374] Substituting the 2n+1 Sigma points into the state transition function, we obtain the 2n+1 Sigma points in... k Predicted value at time :

[0375] (XVII);

[0376] In the formula, Sigma points refer to a set of sample points that are deterministically selected based on the mean and covariance of the current state variables in the unscented transformation step of the unscented Kalman filter algorithm.

[0377] Step 3.4: Calculate the predicted mean and covariance matrix;

[0378] Calculate the system state variables in k Predicted mean at time Covariance Matrix :

[0379] (XVIII);

[0380] In the formula, This represents the predicted value of the state variable at time k based on the observation data at time k-1; k is the time index; n represents the dimension of the vector, i.e., the number of parameters contained in the state variable; i represents the index of the Sigma point. This represents the weight corresponding to the i-th Sigma point; This represents the value of the i-th Sigma point after propagation through the state transition function; Represents the time based on observation data up to time k-1. The predicted state covariance matrix; This represents the covariance weight corresponding to the i-th Sigma point;

[0381] Step 3.5: Perform another UT transformation;

[0382] Based on the obtained predicted mean and covariance matrix, the UT transformation is applied again to calculate... k The Sigma point ζ of the time-time observation (i) and corresponding weights W m (i) and W c (i) This includes: the predicted mean obtained in step 3.4. Covariance Matrix Substituting these new mean and covariance inputs into equation (XVI), the Sigma point ζ for the observation equation at time k is calculated. (i) Weights of Sigma points W m (i) and W c (i) Determined by equation (XVII);

[0383] Step 3.6: Calculate the new predicted value;

[0384] Substituting the new 2n+1 Sigma points into the observation equation, we obtain k Predictive observations for time:

[0385] (XIX);

[0386] Step 3.7: Calculate the predicted mean, variance, and covariance matrix;

[0387] Calculate system observations in k Predicted mean at time ,variance Covariance Matrix :

[0388] (XX);

[0389] In the formula, This represents the predicted mean of the observations at time k based on the observation data at time k-1; This represents the weight corresponding to the i-th Sigma point; This represents the value of the i-th Sigma point after propagation through the nonlinear observation equation; This represents the observation-prediction covariance matrix at time k; This represents the weight corresponding to the i-th Sigma point; This represents the covariance weight corresponding to the i-th Sigma point; Represents the observation noise covariance matrix; Indicates at time The mean of the state prediction;

[0390] Step 3.8: Calculate the Kalman gain matrix;

[0391] (XXI);

[0392] In the formula, This represents the Kalman gain matrix at time k; This represents the observation-prediction covariance matrix at time k; This represents the cross-covariance matrix between the state variables and the observations at time k;

[0393] Step 3.9: Update the system's state variables and covariance:

[0394] (XXII);

[0395] In the formula, That is k Constant air intrusion or leakage flow , .

[0396] The above steps constitute one cyclic step in the UKF inversion of downhole gas intrusion or leakage flow. The solution process is as follows: Figure 4 As shown.

[0397] The specific implementation process of step 3.1 includes:

[0398] Unscented transform of proportionally corrected symmetric sampling, assuming a nonlinear transform y = f ( x ), state vector x Given an n-dimensional random variable and a known state vector. x mean The variance P is used to obtain 2n+1 Sigma points and corresponding weights through an unscented transformation. W ;

[0399] The formula for determining the Sigma point is shown in equation (XXIII):

[0400] (XXIII);

[0401] weight W The formula for obtaining the result is shown in equation (XXIV):

[0402] (XXIV);

[0403] In the formula, (i) Let be the i-th Sigma point; n is the state mean; n is the state dimension; This is a scaling parameter used to reduce the overall prediction error. = (n+κ)-n, κ As a scaling factor, when x When it is a single variable κ =0, when x When there are multiple variables κ =3-n; It refers to the first of the 2n+1 Sigma points generated in the Unscented Transform (UT), i.e., the center point; For matrix The i-th column of the square root matrix obtained by Cholesky decomposition; P is the covariance matrix of the current state; W m (i) , W c (i) These are the weights used for weighting the mean and covariance, respectively. α The scaling factor is 0 ≤ α ≤1, the value reflects the range of the Sigma point set; β It reflects the higher-order properties of the state variables' historical information; for a Gaussian distribution, the optimal value is 2.

[0404] Calculate using equations (XXIII) and (XXIV) k The sigma point of the state quantity at time -1 and the corresponding weights and .

[0405] This embodiment uses basic information from a well (designated as PT-1 well) in a block in western my country and combines it with a multiphase flow model of the wellbore to formulate multiple sets of test data, and analyzes the feasibility of the algorithm under the conditions of single reservoir overflow, simultaneous overflow of multiple reservoirs, and coexistence of overflow and leakage.

[0406] Simulation parameter settings:

[0407] The parameters of the PT-1 well used in the calculation are shown in Table 2. Its well structure is as follows: surface casing Φ273.05mm×1500 m + technical casing Φ200.03 mm×6348 m + open hole Φ171.50 mm×892 m.

[0408] The drill string assembly is as follows: Φ171.50 mm roller cone bit × 0.21 m + Φ127.00 mm adapter × 0.55 m + Φ88.90 mm inclined weight drill pipe × 411.74 m + Φ88.90 mm inclined drill pipe × 4699.10 m + Φ136.00 mm adapter × 0.50 m + Φ114.30 mm inclined drill pipe × 2127.90 m. Figure 5 As shown.

[0409] Table 2

[0410]

[0411] The UKF calculation parameters are shown in Table 3.

[0412] Table 3

[0413]

[0414] Bottom gas inflow inversion calculation:

[0415] First, a common bottom-hole negative pressure differential gas intrusion condition is simulated. Since the formation pressure of this simulated well is known, the gas intrusion flow rate can be calculated using the bottom-hole pressure during the overflow process through the underbalanced gas intrusion non-Darcy flow gas well productivity equation. This calculated flow rate is then substituted into the wellbore multiphase flow model as a boundary condition to obtain the trends in mud pit increment and outlet flow rate. Noise is added to this data to simulate actual field conditions. Finally, the proposed data are used as observations, and the bottom-hole gas intrusion flow rate is used as a state variable. The UKF algorithm is used for inversion, and the results are compared with the set values ​​to verify the feasibility of the method.

[0416] Productivity equations for underbalanced gas-infiltrated non-Darcy flow gas wells (both binomial and exponential forms):

[0417] Binomial: ;

[0418] Exponential: ;

[0419] In the formula, P e Formation pressure, MPa; P w is the bottom hole pressure, MPa; A, B, and C are coefficients related to rock and fluid properties; n is the seepage index, 0.5≤n≤1.

[0420] (1) Constant displacement condition

[0421] The inversion results under constant wellhead discharge conditions are as follows: Figure 6 As shown, from Figure 6 As can be seen, the UKF estimate generally reflects the approximately linear increase in air intrusion flow well, but it fluctuates around the true value. Figure 7 The relative error shown indicates that the error is relatively large in the initial stage of gas intrusion. This is because the initial value of the gas intrusion flow rate is set to 0, thus requiring a certain period of adjustment to approximate the true value. After stabilization, the error remains within 3%, indicating that this method can accurately invert the downhole gas intrusion flow rate.

[0422] Figure 8 The UKF calculation results represent the increase in mud pit volume and the outlet flow rate. Figure 8 (a) shows the increase in mud volume in the mud pit. Figure 8 (b) is a graph showing the export flow rate, from... Figure 8 As can be seen, the error between the calculated and measured parameter values ​​is very small, and it accurately reflects the changing trend of the observed values. Furthermore, the fluctuations in the observed data are eliminated, achieving a noise reduction effect.

[0423] (2) Variable displacement conditions

[0424] The wellhead discharge rate is set to be changed every 3.3 minutes, with the changed discharge rates being 25, 20, and 15 L·s, respectively. -1 The inversion results under these conditions are as follows: Figure 9 As shown, the reduced displacement leads to decreased circulating friction, resulting in lower bottom hole pressure and thus an increase in bottom hole gas influx. The results indicate that the UKF estimate can track abrupt changes in gas influx well, but it exhibits a certain lag and requires adjustment over a period of time. This phenomenon is reflected in the relative error plot (…). Figure 10 This means a sudden increase in error, but the method will react quickly to reduce the error value, and the error in the steady phase is also less than 3%.

[0425] Similarly, Figure 11 This study compares the trends in the incremental volume of the mud pit and the effluent flow rate under variable discharge conditions. Figure 11 (a) shows the increase in mud volume in the mud pit under variable discharge conditions. Figure 11 (a) is the outlet flow diagram under variable displacement conditions, from Figure 11 As can be seen, the UKF calculated value can reflect the changes in the measured value very well, with the error being relatively large only at the point of displacement change.

[0426] (3) Given the formation parameters, invert the formation pressure.

[0427] If key parameters such as reservoir permeability and porosity are known in advance, the state variable can be set as formation pressure. A formation gas invasion characterization equation (underbalanced gas invasion non-Darcy flow gas well productivity equation) can be added to the original observation equation (multiphase flow model), thereby directly inverting the formation pressure using the UKF algorithm. The inversion results are as follows: Figure 12 As shown, the UKF estimate quickly approaches the true value after initial adjustment and fluctuates within a small range around it. The relative error after stabilization is as follows: Figure 13 As shown, all values ​​are within 0.2%, indicating that the method has high feasibility and accuracy when sufficient data is available.

[0428] Inversion calculation of gas influx and leakage flow rates under complex conditions with multiple reservoirs:

[0429] The above analysis shows that the UKF algorithm performs well in parameter inversion after gas invasion in a single reservoir at the bottom of the well. Therefore, we propose the existence of two reservoirs downhole and simulate the simultaneous occurrence of gas invasion and leakage to explore the applicability of this algorithm under complex conditions.

[0430] (1) Dual reservoir simultaneous gas intrusion condition

[0431] First, the depth of the second reservoir was set to 6540 m and the reservoir pressure to 90 MPa. The corresponding observation parameters were calculated using a multiphase flow model, and noise was added to simulate the field conditions.

[0432] In the algorithm's solution process, the gas influx rates of the two reservoirs are used as state variables, which increases the number of inversion parameters and thus increases the difficulty. Furthermore, since the observations are all wellhead parameters, two different combinations of gas influx rates may yield the same wellhead results. Therefore, this invention adds bottomhole pressure as an additional observation parameter (this method is applicable to situations where downhole PWD equipment is available). The inversion results are as follows: Figure 14 As shown, Figure 14 (a) is a schematic diagram showing the actual value and UKF estimate of the first reservoir. Figure 14(b) is a schematic diagram of the true value and UKF estimate of the second reservoir. For the bottom reservoir (first reservoir), the UKF estimate can track the true trend of change well, but it is generally too small. The inversion accuracy of the second reservoir is relatively low. It shows an upward trend in the same way as the true value, but the increase is significantly higher than the true value in the early stage of gas invasion. The algorithm then made adjustments to reduce the estimate, but it still deviates significantly at some time points. Figure 15 (a) shows the error diagram of the gas intrusion flow inversion results for the first reservoir. Figure 15 (b) shows the error diagram of the second reservoir gas intrusion flow inversion results, combined with... Figure 15 Error analysis also shows that after the calculation stabilizes, the gas intrusion flow inversion error of the first reservoir is less than 2%, while the second reservoir has several instances of relatively large errors, but overall it is still within 4%.

[0433] Figure 16 The results of the observation calculations are presented. Figure 16 Figure (a) shows the calculation results of the mud pit increment. Figure 16 Figure (b) shows the calculated export flow rate. Figure 16 Figure (c) shows the calculated bottom hole pressure results, where the mud pit increment and outlet flow rate increase faster compared to the single reservoir scenario. The UKF algorithm can also describe this process well. Furthermore, the calculated bottom hole pressure results also match the trend of the measured values.

[0434] (2) Overflow and leakage coexisting conditions

[0435] The basic parameters of the second reservoir were modified to simulate leakage conditions. Its depth was also set to 6540m, but the leakage pressure was set to 80 MPa. The state variables for UKF calculation were set as the gas influx flow rate of the first reservoir and the leakage flow rate of the second reservoir, while the observations remained the mud pit increment, outlet flow rate, and bottom hole pressure.

[0436] UKF inversion results are as follows Figure 17 As shown, Figure 17 (a) shows the actual value of air intrusion and the estimated value by UKF. Figure 17 Figure (b) shows the actual and UKF estimated values ​​of leakage. For gas-infiltrated reservoirs, the UKF estimated value can also track the growth trend of the actual value, and is generally lower than the actual value. For leaking reservoirs, the gas intrusion at the bottom of the well causes an increase in the fluid velocity in the annulus, which increases the annular friction and leads to a slow increase in wellbore pressure in the pure liquid phase region. This corresponds to the slow increase in leakage flow rate from 0 to 4 minutes in the figure. During this period, the UKF estimated value also increases, but at a relatively large rate. After 4 minutes, due to the gradual migration and expansion of the gas, the wellbore pressure at the leaking layer decreases, resulting in a reduction in leakage flow rate. The UKF algorithm can also simulate this process relatively accurately. The overall error after stabilization is as follows: Figure 18 All of them are below 4%. Figure 18 (a) shows the error diagram of the gas intrusion inversion results. Figure 18 (b) is the error graph of the leakage inversion result.

[0437] In summary, the inversion accuracy of the UKF algorithm decreases somewhat when considering multiple reservoirs, particularly in the inversion results of the second reservoir. However, the overall error remains below 4%, indicating that the method can handle relatively complex situations. Further improvements in the algorithm's accuracy and stability can be achieved by adding additional measured data constraints.

Claims

1. A method for inverting downhole gas intrusion and leakage parameters based on unscented Kalman filtering, characterized in that, include: Step 1: Establish a multiphase flow model for the overflow control wellbore; The process involves establishing a basic mathematical model and assumptions through mathematical modeling; describing the general laws of two-phase flow through basic governing equations; supplementing and specifying the basic governing equations through auxiliary equations; establishing the unique solution to the problem through boundary conditions on the initial state and boundary conditions; dividing the continuous wellbore space and continuous time into a finite number of discrete grid points that can be processed by a computer through grid partitioning of the boundary domain; using the finite difference method to transform partial differential equations that cannot be directly solved into algebraic equations that can be solved at discrete grid points through the discretization of the governing equations; and using the numerical solution of the equation system to design an iterative algorithm to actually solve these equations, thus realizing a complete modeling and solution process for the overflow well control wellbore multiphase flow model from physics to mathematics to computation. Step 2: Establish the inversion model; set the unmeasurable downhole parameters as state variables, and set the mud pit increment and outlet flow rate as observation variables, and establish the state equation and observation equation; Step 3: Solve the inversion model; Generate Sigma points by performing an unscented transformation on the state variables using the unscented Kalman filter algorithm, substitute the obtained Sigma points into the observation equation, call the overflow control wellbore multiphase flow model to calculate the observations and perform filtering updates, obtain and output the inversion results of downhole gas intrusion and leakage parameters based on the filtered and updated state variables, including: bottom hole gas intrusion flow rate, multi-reservoir gas intrusion flow rate, leakage flow rate, and formation pressure; The inversion model takes the overflow process as the analysis object. In the inversion calculation using the unscented Kalman filter algorithm, the air intrusion and leakage flow rates are taken as state variables, and the mud pit increment and outlet flow rate are taken as observations. The corresponding state equation and observation equation are described as Equation (XIV) and Equation (XV): (XIV); (XV); In the formula, k express k time; x k This refers to the current inflow or outflow rate, in m. 3 ·s -1 ; This refers to the current monitoring values ​​of the mud pit's incremental volume and outlet flow rate, in m. 3 ·s -1 ; The measurement noise at the current moment; The system noise at the previous moment; f (·) is the state transition function. x k = x k-1 + w k-1 ; h (·) represents the nonlinear observation equation, i.e., the wellbore multiphase flow model established in step 1; Solving the inversion model; including: Step 3.1: The unscented Kalman filter algorithm uses unscented transformation to handle the nonlinear propagation problem of mean and covariance, that is, to approximate the probability density distribution of the nonlinear function and use a series of deterministic samples to approximate the posterior probability density of the state. Step 3.2: Determine the initial values; Assuming initial state variables x 0 And calculate the corresponding statistical values: (XVI); In the formula, The average of the initial state variables; The initial state quantity covariance; It refers to the expectation in mathematical notation. It is the initial state vector; Step 3.3: Calculate the predicted value; Substituting the 2n+1 Sigma points into the state transition function, we obtain the 2n+1 Sigma points in... k Predicted value at time : (XVII); In the formula, Sigma points refer to a set of sample points that are deterministically selected based on the mean and covariance of the current state variables in the unscented transformation step of the unscented Kalman filter algorithm. Step 3.4: Calculate the predicted mean and covariance matrix; Calculate the system state variables in k Predicted mean at time Covariance Matrix : (XVIII); In the formula, This represents the predicted value of the state variable at time k based on the observation data at time k-1; k is the time index; n represents the dimension of the vector, i.e., the number of parameters contained in the state variable; i represents the index of the Sigma point. This represents the weight corresponding to the i-th Sigma point; This represents the value of the i-th Sigma point after propagation through the state transition function; Represents the time based on observation data up to time k-1. The predicted state covariance matrix; This represents the covariance weight corresponding to the i-th Sigma point; Step 3.5: Perform another UT transformation; Based on the obtained predicted mean and covariance matrix, the UT transformation is applied again to calculate... k The Sigma point ζ of the time-time observation (i) and corresponding weights W m (i) and W c (i) This includes: the predicted mean obtained in step 3.

4. Covariance Matrix Substituting these new mean and covariance inputs into equation (XVI), the Sigma point ζ for the observation equation at time k is calculated. (i) Weights of Sigma points W m (i) and W c (i) Determined by equation (XVII); Step 3.6: Calculate the new predicted value; Substituting the new 2n+1 Sigma points into the observation equation, we obtain k Predictive observations for time: (XIX); Step 3.7: Calculate the predicted mean, variance, and covariance matrix; Calculate system observations in k Predicted mean at time ,variance Covariance Matrix : (XX); In the formula, This represents the predicted mean of the observations at time k based on the observation data at time k-1; This represents the weight corresponding to the i-th Sigma point; This represents the value of the i-th Sigma point after propagation through the nonlinear observation equation; This represents the observation-prediction covariance matrix at time k; This represents the weight corresponding to the i-th Sigma point; This represents the covariance weight corresponding to the i-th Sigma point; Represents the observation noise covariance matrix; Indicates at time The mean of the state prediction; Step 3.8: Calculate the Kalman gain matrix; (XXI); In the formula, This represents the Kalman gain matrix at time k; This represents the observation-prediction covariance matrix at time k; This represents the cross-covariance matrix between the state variables and the observations at time k; Step 3.9: Update the system's state variables and covariance: (XXII); In the formula, That is k Constant air intrusion or leakage flow , .

2. The method for inverting downhole gas intrusion and leakage parameters based on unscented Kalman filtering according to claim 1, characterized in that, Establish basic mathematical models and fundamental assumptions through mathematical modeling; including: The basic mathematical models include drift flow model, homogeneous flow model, and phase-separated flow model; the basic assumptions are as follows: ① The flow of fluid in the annulus is considered as a one-dimensional flow along the wellbore direction; ②The gas and liquid phases reach equilibrium instantaneously in the wellbore; ③Consider the dissolution and precipitation of CO2 and CH4 gases in the drilling fluid, and ignore the chemical reactions between CO2 and CH4 gases and certain components in the drilling fluid; ④ Stable radial heat transfer occurs between the wellbore and the formation.

3. The method for inverting downhole gas intrusion and leakage parameters based on unscented Kalman filtering according to claim 1, characterized in that, The fundamental governing equations include the continuity equation and the momentum equation; Based on the law of conservation of mass, a continuity equation is established along the flow direction in the wellbore annulus, with the coordinates along the flow direction in the wellbore annulus set as follows: z Extracting micro-segments dz The study found that the cross-sectional area of ​​the annular flow was... A Considering the gas precipitation process during its ascent in the wellbore annulus, the continuity equations for the gas-producing zone, the non-producing zone, and the liquid phase are obtained, as shown in equations (I), (II), and (III), respectively: （I）； (II); (III); In the formula, The density of free gas at the temperature and pressure within the annulus, in kg·m³. -3 ; Density of drilling fluid, kg·m -3 ; , The upward return velocities of free gas and drilling fluid, respectively, are in m / s. -1 ; , These are the integrals of free gas and drilling fluid, respectively, and are dimensionless. A Let m be the cross-sectional area of ​​the annulus. 2 , ,in D The inner diameter of the well shaft is in meters (m). d The outer diameter of the drill string for different well sections, in meters; The mass of gas produced per unit time per unit thickness of reservoir, in kg·s -1 ·m -1 ; Rs m represents the solubility of the gas in the drilling fluid. 3 ·m -3 ; The density of the gas under standard conditions is expressed in kg·m³. -3 ; The volume coefficient of the drilling fluid in the local area is dimensionless. According to the law of conservation of momentum: the rate of change of momentum over time within a unit cell is equal to the sum of all external forces acting on that unit cell, resulting in the momentum equation, as shown in equation (IV): ； (IV); In the formula, P Pressure, Pa; The angle between the wellbore axis and the direction of gravity is fr, which refers to the pressure gradient caused by friction.

4. The method for inverting downhole gas intrusion and leakage parameters based on unscented Kalman filtering according to claim 1, characterized in that, Auxiliary equations supplement and specify the basic governing equations; including: The auxiliary equations include the two-phase flow regime determination equation, the velocity equation, the fluid property equation, and the annular friction equation; The current flow state is determined by the two-phase flow state determination equation, and the corresponding velocity equation and annular friction equation are selected accordingly. The convection terms in the continuity and momentum equations are supplemented by introducing the distribution coefficient and slip velocity from the drift flow model in the velocity equation. The annular friction equation provides friction models under different flow conditions, enabling the momentum equation to accurately reflect the influence of wellbore roughness and flow state on pressure drop.

5. The method for inverting downhole gas intrusion and leakage parameters based on unscented Kalman filtering according to claim 1, characterized in that, Boundary conditions establish a unique solution to the problem by imposing constraints on the initial state and boundary conditions; these include: The initial state is defined as the distribution of parameters within the wellbore at the start of the simulation; Boundary conditions are defined on the boundaries of the wellbore, and the parameters are constraints that change over time. Boundary domain mesh generation divides the continuous wellbore space and continuous time into a finite number of discrete grid points that a computer can process; including: Spatial domain grid partitioning includes: Mesh generation is performed using a fixed spatial step size method, where the length of any given mesh is... for: (V)； , It refers to the spatial coordinates of two adjacent grid points in the direction of wellbore depth; Then the total number of spatial grids for: (WE); H refers to the total depth of the well; Time-domain gridding includes: The time step employs a non-uniform grid partitioning method, based on the velocity of the free gas. and the length of the spatial grid at that location Relationship obtained by time step : (VII)。 6. The method for inverting downhole gas intrusion and leakage parameters based on unscented Kalman filtering according to claim 3, characterized in that, Discretization of the governing equations utilizes the finite difference method to transform partially differential equations that cannot be directly solved into algebraic equations solvable at discrete grid points; including: The partial differential forms of the continuity equation and the momentum equation are unified as equation (VIII): (VIII); Where X represents the cumulative term of a physical quantity, Y represents the current term of a physical quantity, and z and t are the spatial coordinate and time coordinate, respectively; The four-point difference scheme equation of equation (VIII) is shown in equation (IX): (IX); Where n and j are the time step index and the spatial grid point index, respectively; and These are the spatial step size and the time step size, respectively. This represents the value of a physical quantity at spatial location j+1 and time point n+1. Let the value of the physical quantity be at spatial location j and time point n+1; and The values ​​of physical quantities at spatial points n, j, and j+1 at the same time point; , Let n+1 be the value of the physical quantity at spatial points j and j+1 at the same time point; From equation (IX), we get: For the continuity equation of the gas-phase producing layer, the difference scheme is shown in equation (X): ； （X）； ； For the continuity equation of the non-productive gas phase, the difference scheme is shown in equation (XI): ； (XI); For the continuity equation of the liquid phase, the difference scheme is shown in equation (XII): ； (XII); The difference scheme for the momentum equation is shown in equation (XIII): ； ； (XIII); ； 。 7. The method for inverting downhole gas intrusion and leakage parameters based on unscented Kalman filtering according to claim 1, characterized in that, Numerical solutions to systems of equations are used to design an iterative algorithm to actually solve these algebraic equations, including: When calculating wellbore pressure, an initial pressure value is first assigned to the pressure node to be calculated. This initial pressure value is then substituted into the equation system to obtain the relevant parameters. The flow regime is determined according to the flow pattern discrimination equation. The pressure of the pressure node to be calculated is then calculated using the flow model. The calculated pressure is compared with the initial pressure value. If the absolute deviation between the calculated pressure and the initial pressure value is less than the specified error range, the same method is used to calculate the next pressure node. Otherwise, a new initial pressure value is assigned and the calculation of this pressure node continues until the constraint conditions are met to calculate the next pressure node. The numerical solution of the system of equations is used to design an iterative algorithm to actually solve these algebraic equations. The specific implementation steps include: 1) Input basic parameters, which include the following: Wellbore geometric parameters: wellbore depth H, wellbore inner diameter D, drill string outer diameter d; Fluid physical properties: Drilling fluid density Drilling fluid plastic viscosity, drilling fluid dynamic shear stress, drilling fluid apparent viscosity, and gas density. Gas solubility Rs, gas viscosity σ, the surface tension of the gas; Other parameters: Drilling fluid displacement Inlet temperature, geothermal gradient, and air intrusion flow rate Spatial step size ∆z, time step size ∆t, error threshold ε; 2) Given the initial and boundary conditions, t = t 0 When the value is 0, the pressure at each node is calculated starting from the wellhead. And velocity parameters and fluid properties; velocity parameters include apparent gas phase velocity. Apparent velocity in the liquid phase Gas phase true velocity True velocity of liquid phase Fluid properties include drilling fluid density. Drilling fluid plastic viscosity, dynamic shear force, and gas density Gas viscosity ; 3) Based on the flow velocity at the bottom node and the length of the divided cell, substitute the flow velocity and cell length into equation (VII) to obtain the time step. ∆ t0 ; 4) Assuming the first time node t 1 = t 0 + ∆ t0 The bottom pressure at the well is Calculate the gas solubility, velocity parameters, and fluid properties at the bottom node; 5) Assume the pressure at (1, 1) is Calculate the gas solubility, velocity parameters, and fluid properties at this node; 6) Assume the gas phase volume fraction at (1,1) is Based on the continuity equation, calculate the apparent flow velocities of the gas and liquid phases at (1,1). , The gas phase volume fraction at this node is calculated based on the velocity equation. ; 7) If satisfied If ≤ε, then it is considered If the assumption is reasonable, otherwise it needs to be re-assumed until the condition is met; ε is the error threshold. 8) The obtained Substitute the relevant parameters into the momentum equation to calculate the pressure at (1,1). ; 9) If satisfied If ≤ε, then it is considered If the assumption is reasonable, then re-assume. Repeat steps 5)-9) until the condition is met; 10) The obtained Using the known pressure as the next node, repeat the above calculation until the wellhead. If the absolute difference between the calculated wellhead pressure and atmospheric pressure meets the accuracy requirement, then it is considered... If the assumption is reasonable, then re-assume. Repeat steps 4)-10) until the conditions are met; among them, calculate the absolute value of the difference between the wellhead pressure and the atmospheric pressure. If the result is less than the error threshold ε, the accuracy condition is met. By following the steps above, the parameters of all spatial points within the first time point are calculated. These parameters are then used as known parameters for the next time point, and so on, to calculate... t 2 , t 3 … t n All parameters within a given time point.

8. The method for inverting downhole gas intrusion and leakage parameters based on unscented Kalman filtering according to claim 1, characterized in that, The specific implementation process of step 3.1 includes: Unscented transform of proportionally corrected symmetric sampling, assuming a nonlinear transform y = f ( x ), state vector x Given an n-dimensional random variable and a known state vector. x mean The variance P is used to obtain 2n+1 Sigma points and corresponding weights through an unscented transformation. W ; The formula for determining the Sigma point is shown in equation (XXIII): (XXIII); weight W The formula for obtaining the result is shown in equation (XXIV): (XXIV); In the formula, (i) Let be the i-th Sigma point; n is the state mean; n is the state dimension; This is a scaling parameter used to reduce the overall prediction error. = (n+κ)-n, κ As a scaling factor, when x When it is a single variable κ =0, when x When there are multiple variables κ =3-n; It refers to the first of the 2n+1 Sigma points generated in the Unscented Transform (UT), i.e., the center point; For matrix The i-th column of the square root matrix obtained by Cholesky decomposition; P is the covariance matrix of the current state; W m (i) , W c (i) These are the weights used for weighting the mean and covariance, respectively. α The scaling factor is 0 ≤ α ≤1, β Higher-order characteristics that reflect historical information of state variables; Calculate using equations (XXIII) and (XXIV) k The sigma point of the state quantity at time -1 and the corresponding weights and .

Images