Well drilling complex event real-time early warning method based on physical attenuation-particle swarm coupling

Abstract

The invention discloses a drilling complex event real-time early warning method based on physical attenuation-particle swarm coupling, and belongs to the technical field of petroleum drilling safety monitoring. The method comprises the steps that firstly, key drilling parameters such as mud pressure, mechanical drilling speed and torque are collected, and physical attenuation reference models of the mud pressure, the ROP and the torque are established respectively; taking an attenuation coefficient and an auxiliary parameter as particle dimensions, constructing a coupling fitness function, and performing real-time optimization by using an improved particle swarm algorithm with a self-adaptive inertia weight and Gaussian disturbance; and directly mapping hyper-parameters of the early warning model by using the optimized global optimal particles, deducing the event probability online and outputting a multi-stage early warning signal. According to the method, a physical mechanism and intelligent optimization are deeply coupled, second-level, high-precision and self-adaptive early warning of complex events such as well leakage, drill jamming and well kick is achieved, and the missing report rate and the false report rate are remarkably reduced.

Description

Technical Field

[0001] This invention relates to the field of oil and gas drilling engineering and intelligent optimization algorithm technology, and in particular to a real-time early warning method for complex drilling events based on physical decay-particle swarm coupling. Background Technology

[0002] With the large-scale development of deep, ultra-deep, and shale gas wells, complex downhole events (lost circulation, stuck pipe, well kick, etc.) occur frequently. If not detected in time, these events can lead to significant economic losses or even safety accidents. Existing technologies have the following main shortcomings:

[0003] 1. Traditional threshold methods rely on manually setting fixed thresholds, which are difficult to adapt to dynamic changes in formation parameters and operating conditions, resulting in high rates of missed and false alarms.

[0004] 2. Conventional machine learning models only perform statistical fitting on monitoring data, lack physical constraints, and have poor generalization ability;

[0005] 3. Single particle swarm optimization is mostly used for offline training, and it does not organically combine decay mechanism, real-time constraints and online learning, resulting in insufficient applicability in the field.

[0006] Therefore, there is an urgent need for a drilling complex event early warning method that can integrate physical attenuation mechanisms with intelligent optimization algorithms, adapt in real time, and meet engineering constraints. Summary of the Invention

[0007] The purpose of this invention is to provide a real-time early warning method for complex drilling events based on physical decay-particle swarm coupling, thereby solving the aforementioned problems existing in the prior art.

[0008] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0009] A real-time early warning method for complex drilling events based on physical decay-particle swarm optimization includes the following steps:

[0010] Step S1: During the drilling process, D-dimensional drilling parameters are collected in a fixed time sequence to form a time series dataset X = {x(1), x(2), ..., x(T)}, where x(t) ∈ R^D includes at least mud pressure, mechanical drilling speed ROP, and torque T;

[0011] Step S2: Construct three types of physical attenuation benchmark models respectively:

[0012] S21, Mud pressure decay model: P(t) = P0·e^(-λt) + P_min,

[0013] S22, Mechanical drilling rate attenuation model ROP(h)=ROP0·(1-α·h),

[0014] S23, Torque decay model T(t)=T0-β·t+δ(t),

[0015] λ, α, and β are then used as attenuation coefficients to be optimized.

[0016] Step S3: Initialize N particles using the attenuation coefficients λ, α, β and D-3 auxiliary parameters as particle dimensions. The position vector of each particle represents a set of parameter combinations to be optimized.

[0017] Step S4: Establish the coupling fitness function F, which is composed of the following weighted sub-functions:

[0018] F1 = 1 / MSE1, where MSE1 is the mean square error between the predicted and measured values ​​of mud pressure.

[0019] F2 = 1 / MSE2, where MSE2 is the mean square error between the ROP predicted value and the measured value;

[0020] F3 = 1 / MSE3, where MSE3 is the mean square error between the predicted torque value and the measured torque value;

[0021] F4 = Penalty terms for λ, α, and β within the preset engineering safety range;

[0022] Step S5: Utilize the adaptive inertia weight ω t An improved particle swarm optimization algorithm with Gaussian perturbation term ε iteratively updates particle positions to maximize F and obtain the globally optimal particle gbest.

[0023] Step S6: Map gbest to the hyperparameters of the real-time early warning model M, perform online inference of the drilling parameters at the current time t, and output the event probability P_event∈[0,1];

[0024] Step S7: When P_event≥θ, trigger a multi-level early warning signal and provide the physical attenuation coefficient offset Δλ, Δα or Δβ for the corresponding working condition to guide the adjustment of drilling parameters.

[0025] Preferably, step S2 includes:

[0026] The physical feasible regions of λ, α, and β are updated in real time using the real-time mud viscosity μ(t) and well depth h(t) such that λ∈[λ_min(μ),λ_max(μ)], α∈[α_min(h),α_max(h)], and β∈[β_min(h),β_max(h)].

[0027] Preferably, the adaptive inertia weight ω_t in step S5 satisfies:

[0028] ω_t=ω_max-(ω_max-ω_min)·t / T_max, and ω_max=0.9, ω_min=0.4.

[0029] Preferably, the Gaussian perturbation term ε in step S5 is ε ~ N(0,σ) 2 ), and σ=0.1·|gbest|.

[0030] Preferably, the penalty term in step S4 is:

[0031] Penalty=∑_{k∈{λ,α,β}}max(0,(k-k_bound) 2 / σ_k 2 ),

[0032] k_bound is the engineering constraint boundary with the corresponding coefficient, and σ_k is the tolerance coefficient.

[0033] Preferably, the multi-level warning signals in step S7 include:

[0034] Yellow alert: Δλ / λ_normal > 1.2 and P_event ≥ 0.7;

[0035] Orange alert: Δα / α_normal > 1.3 and P_event ≥ 0.8;

[0036] Red alert: Δβ / β_normal < 0.5 and torque surge > 20% and P_event ≥ 0.9.

[0037] Preferably, the real-time early warning model M is a support vector machine or a deep neural network, and the parameters of its kernel function or activation function are directly assigned by gbest.

[0038] Preferably, the sampling frequency of step S1 is ≥10Hz, and the data preprocessing includes 3σ anomaly removal and linear interpolation completion.

[0039] Preferably, the online inference latency of step S6 is <200ms, and it is achieved through parallel GPU acceleration.

[0040] Furthermore, a computer-readable storage medium having a program stored thereon, which, when executed by a processor, implements the steps of the above method.

[0041] The beneficial effects of this invention are as follows: This invention discloses a real-time early warning method for complex drilling events based on physical attenuation-particle swarm optimization, which has the following beneficial effects:

[0042] Physics-data dual-drive: The physical decay model provides prior constraints for particle swarm optimization, making the optimization results interpretable and physically reasonable, and improving the early warning accuracy by more than 30% compared with traditional methods.

[0043] Adaptive noise resistance: Adaptive inertial weights and Gaussian perturbation terms are introduced to effectively suppress the impact of high-dimensional noise on optimization, reducing the false negative rate from 15% to 5% under complex geological conditions.

[0044] High real-time performance: Parallel GPU acceleration and lightweight model mapping reduce single-well data processing latency to <200ms, meeting the on-site second-level early warning requirements.

[0045] Engineering constraint embedding: Engineering boundaries such as mud pressure safety range and ROP limit value are written during the particle initialization stage, avoiding infeasible solutions and improving the success rate of on-site implementation.

[0046] Multi-level early warning linkage: Triggering yellow / orange / red three-level alarms based on the attenuation coefficient offset, and providing targeted adjustment suggestions to achieve closed-loop management of "monitoring-early warning-intervention". Attached Figure Description

[0047] Figure 1 This is a schematic diagram of the overall process of the present invention;

[0048] Figure 2 This is a schematic diagram of the iterative process of the improved particle swarm optimization algorithm of this invention. Detailed Implementation

[0049] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0050] Reference Figure 1 and Figure 2 The method for real-time early warning of complex drilling events based on physical decay-particle swarm optimization includes the following steps:

[0051] Step S1: During the drilling process, D-dimensional drilling parameters are collected in a fixed time sequence to form a time series dataset X = {x(1), x(2), ..., x(T)}, where x(t) ∈ R^D includes at least mud pressure, mechanical drilling speed ROP, and torque T;

[0052] 1. Sensor selection and deployment:

[0053] Seven types of industrial-grade sensors were deployed on-site, forming a raw data channel with a dimension of D≥7:

[0054] Mud pressure: Honeywell DS type, range 0–70MPa, error ±0.1%F·S, installed in the mud pump outlet manifold;

[0055] Mechanical drilling rate ROP: calculated from the depth encoder (depth accuracy ±0.5cm / 30m) and the tripping timestamp, outputting m·h⁻¹;

[0056] Torque: Honeywell 65 type, range 0–30kN·m, error ±0.1%F·S, mounted on turntable drive shaft;

[0057] Drilling pressure and rotation speed: Feiboer BBVV5500 model, drilling pressure 0–500kN, rotation speed 0–300rpm, error ±0.1%;

[0058] Mud density and temperature: Coriolis density meter + PT100 temperature probe, error ±0.001 g·cm³ -3 ±0.1℃;

[0059] Drilling fluid pool volume: Ultrasonic level gauge, error ±0.5%.

[0060] 2. Sampling Strategy

[0061] Fixed timing: Samples are continuously taken at a uniform frequency of 10Hz to form the original time series.

[0062] X_raw = {x_raw(1), x_raw(2), ..., x_raw(T)}, where T is the drilling time (seconds).

[0063] Clock synchronization: The PTP protocol is used to synchronize the time of all sensors, with a timestamp accuracy of <1ms, ensuring cross-channel data alignment.

[0064] 3. Channel Mapping and Dimension Definition: The above 7 physical quantities are mapped to a 7-dimensional vector: x_raw(t)=[P_mud(t),ROP(t),T(t),WOB(t),

[0065] If higher dimensions are needed, output flow rate, conductivity, etc. can be added, in which case D>7.

[0066] 4. Data Preprocessing (Technical Disclosure: "Data Cleaning and Normalization" section)

[0067] Outlier removal: Outliers are removed for each dimension of data using the 3σ criterion.

[0068] Missing points interpolation: Linear interpolation fills in ≤2% of missing points;

[0069] Normalization: Max-Min Normalization

[0070] x_norm,i(t) = (x_raw,i(t) – x_min,i) / (x_max,i – x_min,i), i = 1…D. Map each dimension to [0,1] to eliminate dimensional differences and obtain the final input.

[0071] X={x(1),x(2),…,x(T)}, where x(t)∈[0,1]^D.

[0072] 5. Interface with subsequent steps

[0073] After processing by S1, X is directly used as the input data stream for steps S2 (physical decay model) and S3–S7 (particle swarm optimization and early warning), meeting the requirement of "real-time data processing latency <200ms" in the technical disclosure.

[0074] Step S2: Construct three types of physical attenuation benchmark models respectively:

[0075] S21, Mud pressure decay model: P(t) = P0·e^(-λt) + P_min,

[0076] S22, Mechanical drilling rate attenuation model ROP(h)=ROP0·(1-α·h),

[0077] S23, Torque decay model T(t)=T0-β·t+δ(t),

[0078] λ, α, and β are then used as attenuation coefficients to be optimized.

[0079] This step is specifically as follows:

[0080] 1. Overall Model Concept

[0081] Step S2 involves independently modeling these three scenarios as the "physical baseline" for subsequent particle swarm optimization. By setting the decay coefficients λ, α, and β as particle dimensions, a dual-driven approach of "mechanism-data" is achieved.

[0082] 2. Construction and Calibration of Three Types of Physical Attenuation Reference Models

[0083] 2.1 Mud Pressure Attenuation Model: P(t) = P0·e^(-λt) + P_min

[0084] Physical meaning:

[0085] P0: Initial riser pressure (measured value at t=0).

[0086] λ: Attenuation coefficient considering the combined effects of frictional resistance along the drilling path, drilling fluid viscosity μ(t), and wellbore roughness; unit: s. -1 .

[0087] P_min: Minimum circulating pressure loss at the bottom of the well, obtained from the annular pressure loss calculation formula or actual measurement.

[0088] Calibration method:

[0089] During normal drilling, a sliding window (e.g., 120s) is taken, and the measured Pt curve is fitted using least squares to obtain the "normal reference value λ_norm".

[0090] λ_norm is used as the initial center value for particle swarm search, and the search interval is [0.5λ_norm, 2λ_norm] to take into account formation variations.

[0091] 2.2 Mechanical Drilling Rate Decay Model: ROP(h) = ROP0·(1-α·h)

[0092] Physical meaning:

[0093] ROP0: Instantaneous drilling velocity (m·h⁻¹) when the drill bit first contacts the formation.

[0094] α: A comprehensive coefficient reflecting drill bit wear, formation abrasiveness, and bottom hole pressure differential; unit: m. -1 .

[0095] h: Relative well depth increment (m), with a window of 0–200m, which can be slid as needed.

[0096] Calibration method:

[0097] Within the 0–50m interval of the newly drilled bit, ROP is calculated in segments with a 1m step size, and α_norm is obtained by linear regression.

[0098] Considering that high drilling pressure will reduce α, it is corrected in real time: α_allow(h)=α_norm·(1–0.1·ΔWOB / WOB_ref).

[0099] 2.3 Torque decay model: T(t) = T0 - β·t + δ(t)

[0100] Physical meaning:

[0101] T0: Initial turntable torque (kN·m).

[0102] β: Linear decay slope caused by friction between drill string and wellbore / cuttings; unit: kN·m·s -1 .

[0103] δ(t): Normal fluctuation noise term, following the order N(0,σ) 2 ), where σ is the standard deviation of historical torque.

[0104] Calibration method:

[0105] Linear fitting was performed on 300s data from wells with good cleanliness and no abnormal sticking to obtain β_norm.

[0106] The search range for β is set to [0.3β_norm, 1.5β_norm] to accommodate different wellbore cleanliness levels.

[0107] 3. Attenuation coefficient → Particle dimension mapping

[0108] The "Particle Coding" section of the technical disclosure states that "each particle position vector represents a set of parameters to be optimized." Therefore, in step S3:

[0109] The particle dimension d = 3 + (D - 3) is fixed in the first three dimensions as [λ, α, β], and the last (D - 3) dimensions are auxiliary parameters (such as mud density, temperature, etc.).

[0110] During initialization, the nominal values ​​of λ, α, and β, λ_norm, α_norm, and β_norm, are used as particle centers, and a search boundary of ±50% is set based on engineering experience to ensure that the particle swarm converges quickly within the physically feasible region.

[0111] 4. Real-time update mechanism

[0112] Every 30 seconds, the latest window data is used to refit λ_norm, α_norm, and β_norm, which serve as the "reference center" for the next round of particle swarm optimization, thus achieving "real-time updating of model parameters".

[0113] When an anomaly is detected (such as a sudden increase in torque), immediately expand the search range of β to [0.1β_norm, 3β_norm] to capture signs of a stuck drill.

[0114] Through the above steps, step S2 not only provides the precise mathematical form of the three types of decay benchmark models, but also clarifies the real-time calibration, update and particle dimension mapping methods of λ, α and β, providing a reliable physical benchmark for subsequent particle swarm optimization.

[0115] Step S3: Initialize N particles using the attenuation coefficients λ, α, β and D-3 auxiliary parameters as particle dimensions. The position vector of each particle represents a set of parameter combinations to be optimized.

[0116] The meaning of "D-3" in D-3 auxiliary parameters

[0117] D: refers to the total number of dimensions of all drilling parameter vectors collected in step S1 (≥7, 7 dimensions or more are taken on site).

[0118] 3: refers to the three physical attenuation coefficients λ, α, and β, which have been modeled separately in step S2.

[0119] D-3: Therefore, it represents "the remaining (D-3) auxiliary parameters other than λ, α, and β", such as mud density, temperature, plastic viscosity, rotary table rotation speed offset, drilling pressure compensation value, etc., which are used to supplement the search space of the particle swarm algorithm.

[0120] 1. Determining the Dimensions of Particles

[0121] The technical disclosure directly maps "D-dimensional drilling parameters" to particle dimensions. Based on this, step S3 is explicitly divided into:

[0122] The first three dimensions: physical attenuation coefficients λ, α, and β (defined in step S2);

[0123] The later (D–3) dimension: auxiliary parameters—including mud density ρ, temperature θ, plastic viscosity μ, rotary table speed RPM offset, and drilling pressure WOB compensation value, etc., have a significant impact on real-time working conditions but are difficult to describe with a single formula.

[0124] Therefore, the particle dimension d satisfies

[0125] d = 3 + (D – 3) = D (≥ 7).

[0126] 2. Particle structure encoding

[0127] The position vector Xi of each particle i is encoded using continuous real numbers:

[0128]

[0129] Where λi, αi, and βi are the attenuation coefficients to be optimized, and the rest are auxiliary parameters.

[0130] The velocity vector Vi is initialized synchronously, with the same dimension as Xi.

[0131] 3. Initialization Strategy

[0132] 3.1 Embedding of Physically Feasible Region

[0133] To avoid generating infeasible solutions, the following boundaries are written during the particle swarm initialization phase:

[0134] λ∈[λ_min,λ_max]: given by the annular pressure loss formula and real-time calculation of mud viscosity;

[0135] α∈[α_min,α_max]: Determined based on drill bit type and formation abrasiveness curve;

[0136] β∈[β_min,β_max]: obtained by linear fitting ±3σ of historical torque segments without abnormalities;

[0137] Auxiliary parameters: such as ρ∈[ρ_l,ρ_u] is the safety window for mud density, and μ is determined by the real-time upper and lower limits of the on-site rheometer.

[0138] 3.2 Population size N

[0139] Field verification: When the data dimension d ≤ 15 for a single well, N = 40 can achieve convergence within 150 iterations;

[0140] Large-scale well sites (d>20) can use N=60–80, and GPU parallelism can be used to maintain real-time performance.

[0141] 3.3 Initial Position and Velocity

[0142] Location: λ, α, and β are uniformly randomized by "nominal value ±10%"; auxiliary parameters are sampled by Gaussian distribution by "real-time on-site measurement value ±5%" to ensure that the initial particles are close to the actual working conditions.

[0143] Speed: Limited to [–v_max, v_max], where v_max is 20% of the feasible region width of the corresponding dimension to prevent out-of-bounds errors in a single iteration.

[0144] 4. Interface with subsequent steps

[0145] After initialization, all particles Xi and Vi are directly fed into the coupling fitness function calculation in step S4; when the physical nominal value in step S2 is refreshed, only the center positions of λi, αi, and βi are translated to maintain population diversity.

[0146] With the above design, step S3 is completed in the "particle encoding" stage:

[0147] The dimensions are clearly defined;

[0148] One-time embedding of the physical and engineering boundaries;

[0149] The initial population covers both the feasible region and is close to reality, significantly shortening the convergence time and improving the global optimization success rate.

[0150] Step S4: Establish the coupling fitness function F, which is composed of the following weighted sub-functions:

[0151] F1 = 1 / MSE1, where MSE1 is the mean square error between the predicted and measured values ​​of mud pressure.

[0152] F2 = 1 / MSE2, where MSE2 is the mean square error between the ROP predicted value and the measured value;

[0153] F3 = 1 / MSE3, where MSE3 is the mean square error between the predicted torque value and the measured torque value;

[0154] F4 = Penalty terms for λ, α, and β within the preset engineering safety range;

[0155] 1. Overall Approach

[0156] Step S4 combines the two to construct a coupled fitness function F with the dual objectives of "error minimization + engineering feasibility", which ensures that the particle swarm can guarantee prediction accuracy and not violate the field safety / process boundaries when searching for the optimal λ, α, β and auxiliary parameters.

[0157] 2. Sub-function construction and mathematical expression

[0158] 2.1 Pressure Error Subfunction F1

[0159] F1 = 1 / MSE1

[0160] MSE1=(1 / K)Σ_{k=1}^{K}[P_pred(t_k)–P_real(t_k)] 2

[0161] P_pred is calculated by substituting λ from the current particle into P(t) = P0·e^(–λt) + P_min;

[0162] K is the length of the sliding window (taken as 120s × 10Hz = 1200 points on site).

[0163] 2.2 ROP Error Subfunction F2

[0164] F 22 =1 / MSE2

[0165] MSE2=(1 / M)Σ_{m=1}^{M}[ROP_pred(h_m)–ROP_real(h_m)] 2

[0166] ROP_pred is calculated by substituting α from the current particle into ROP(h) = ROP0(1–αh);

[0167] M is the number of samples within the current well depth window (typically 200m × 1 point / m = 200 points).

[0168] 2.3 Torque Error Subfunction F3

[0169] F3 = 1 / MSE3

[0170] MSE3=(1 / L)Σ_{l=1}^{L}[T_pred(t_l)–T_real(t_l)] 2

[0171] T_pred is calculated by substituting β in the current particle into T(t) = T0 – βt;

[0172] L is set to 300s × 10Hz = 3000 points.

[0173] 2.4 Engineering safety penalty item F4

[0174] F4=–γ·Penalty

[0175] Penalty=Σ_{j∈{λ,α,β}}w_j·max(0,(x_j–x_{j,bound}) 2 / σ_j 2 )

[0176] x_j is the current dimension value of the particle;

[0177] x_{j,bound} is the engineering safety boundary (λ_bound, α_bound, β_bound are dynamically given by the field mud properties, drill bit manual, and wellbore cleanliness curve);

[0178] σ_j is the tolerance coefficient (taken as 5% of the boundary width on site);

[0179] w_j is the weighting coefficient, with the default w_λ=w_α=w_β=1, which can be adjusted on-site;

[0180] γ is the penalty intensity factor, with an empirical value of 10. 3 -10 4 This ensures that fitness is significantly reduced immediately if the boundary is crossed.

[0181] 3. Weighting and Normalization

[0182] To accommodate errors of different dimensions, min-max normalization is performed on F1, F2, and F3:

[0183] F_i' = (F_i – F_{i,min}) / (F_{i,max} – F_{i,min}), i = 1, 2, 3, final coupling fitness

[0184] F = w1F1' + w2F2' + w3F3' + F4

[0185] The default settings on site are w1=w2=w3=1, which can be adjusted according to the sensitivity of the working conditions (e.g., increase w1 when well leakage is frequent).

[0186] 4. Calculation process (corresponding to the flowchart of "Fitness Function Design" in the handover document)

[0187] ① Read the current particle position Xi = [λ, α, β, ...];

[0188] ② Substituting into the three types of attenuation formulas of S2, we obtain P_pred, ROP_pred, and T_pred;

[0189] ③ Calculate MSE1, MSE2, and MSE3 based on real-time measurements;

[0190] ④ Check whether λ, α, and β exceed the engineering safety range, and calculate Penalty;

[0191] ⑤ Synthesize F according to step 3 and return it to the main loop of the particle swarm.

[0192] 5. On-site example

[0193] In a certain shale gas well section:

[0194] λ_bound = [0.002, 0.008]s -1 α_bound = [0.0005, 0.002]m -1 , β_bound=[0.01,0.05]kN·m·s -1 ;

[0195] γ = 1 × 10 4 ;

[0196] When λ = 0.009s -1 When the boundary is crossed, Penalty = 1 × 10 4 ×(0.009-0.008) 2 / 0.0001 2 =100, which significantly lowers F, forcing the particle swarm to migrate toward the feasible region.

[0197] Through the above coupling design, step S4 achieves simultaneous optimization of "maximizing prediction accuracy" and "engineering feasibility", laying the foundation for subsequent rapid convergence of particle swarm and reliable early warning.

[0198] Step S5: Utilize the adaptive inertia weight ω t An improved particle swarm optimization algorithm with Gaussian perturbation term ε iteratively updates particle positions to maximize F and obtain the globally optimal particle gbest.

[0199] 1. Algorithm Framework

[0200] Traditional Particle Swarm Optimizers (PSOs) are prone to getting trapped in local optima in high-dimensional, noisy drilling data. Step S5 embeds both "adaptive inertial weights" and "Gaussian perturbations" into the standard PSO, forming an "improved particle swarm optimization-physics coupling" closed-loop optimizer. The purpose is to:

[0201] Early-stage global exploration is necessary to prevent premature convergence.

[0202] Later, local fine-tuning search improves the accuracy of the optimal particle.

[0203] With coupling fitness F as the objective, the global optimal particle gbest is output after maximizing F.

[0204] 2. Improve the speed update formula

[0205] 2.1 Adaptive Inertia Weight ω_t

[0206] ω_t=ω_max–t·(ω_max–ω_min) / T_max

[0207] ω_max = 0.9 (global search phase);

[0208] ω_min = 0.4 (partial development stage);

[0209] t: Current iteration number; T_max: Maximum number of iterations (taken as 150 in this case).

[0210] 2.2 Gaussian perturbation term ε

[0211] ε~N(0,σ 2 ),in

[0212] σ=0.1·|gbest|

[0213] |gbest| is the L2 norm of the current globally optimal particle position vector;

[0214] Noise only affects the cognitive part of the speed update, thus avoiding disruption to group information sharing.

[0215] 2.3 Complete velocity / position iterative formula

[0216] v_i^{t+1}=ω_t v_i^t+c1r1(pbest_i–x_i^t)

[0217] +c2r2(gbest–x_i^t+ε) (1)

[0218] x_i^{t+1}=x_i^t+v_i^{t+1} (2)

[0219] c1 = c2 = 2 (default value in the disclosure document);

[0220] r1,r2~U(0,1);

[0221] Velocity limit: |v_i,d|≤v_max,d, where v_max,d is 20% of the width of the feasible region in the corresponding dimension.

[0222] 3. Boundary handling

[0223] If x_i,d exceeds the safe engineering range [x_min,d,x_max,d], then the "bounce" strategy is executed:

[0224] x_i,d = 2·x_bound – x_i,d, and let v_i,d = –v_i,d, to ensure that the particle always falls within the feasible region.

[0225] 4. Termination Criteria

[0226] Use two-condition OR logic:

[0227] The rate of change of F(gbest) is <0.1% after 10 consecutive iterations;

[0228] The maximum number of iterations T_max = 150 is reached.

[0229] Stop when any condition is met, and output gbest as the optimal parameter combination [λ*, α*, β*, ...].

[0230] 5. Parallel acceleration

[0231] GPU-CUDA evaluates the fitness of all particles in parallel, with a single generation taking less than 1ms;

[0232] It runs asynchronously with the on-site data acquisition thread to ensure an overall latency of <200ms.

[0233] 6. On-site verification examples

[0234] Parameters of a shale gas horizontal well: d = 10, N = 40, T_max = 150.

[0235] The ω_t value linearly decreased from 0.9 to 0.73 in the first 50 generations, indicating high population diversity.

[0236] From generation 51 to 150, ω_t continues to decrease to 0.4, enabling fine-grained local search;

[0237] After introducing ε, the number of times the local optimum was escaped increased from 2 to 7, and the final F was improved by 18%.

[0238] With the above improvements, step S5 can still stably converge to the globally optimal particle gbest in a noisy, high-dimensional, and strongly constrained drilling environment, providing high-precision parameters for the subsequent real-time early warning model.

[0239] Step S6: Map gbest to the hyperparameters of the real-time early warning model M, perform online inference of the drilling parameters at the current time t, and output the event probability P_event∈[0,1];

[0240] Step S6: Online Inference and Model Mapping

[0241] 1. Hyperparameter mapping

[0242] The globally optimal particle gbest=[λ*,α*,β*,…] is directly decoded into the hyperparameters of the real-time early warning model M:

[0243] If M is selected as SVM, gbest provides the penalty factor C and kernel parameter γ;

[0244] If M is selected as a deep neural network, gbest provides the learning rate η, the number of hidden layer nodes, the dropout rate, etc.

[0245] Mapping rules: Linearly map each dimension of gbest to the corresponding reasonable engineering range of hyperparameters, such as C = 10^(λ*·2+1), γ = 10^(α*·1.5-3), etc. (on-site calibration of linear scaling factors).

[0246] 2. Online Inference Process

[0247] a) Real-time input: D-dimensional normalized eigenvector x(t) at the current time t.

[0248] b) Forward computation: Model M outputs the soft probability of the event P_event=M(x(t); gbest)∈[0,1].

[0249] c) Parallel acceleration: GPU-batch size = 1, single inference latency < 20ms; overall link (including data I / O) < 200ms, meeting the technical disclosure document's requirement of "single well data processing latency < 200ms".

[0250] 3. Probability smoothing

[0251] To prevent momentary fluctuations, a 5-point moving average is used:

[0252] P_smooth(t)=0.5·P_event(t)+0.3·P_event(t-1)+0.2·P_event(t-2).

[0253] Step S7: When P_event≥θ, trigger a multi-level early warning signal and provide the physical attenuation coefficient offset Δλ, Δα or Δβ for the corresponding working condition to guide the adjustment of drilling parameters.

[0254] 1. Threshold system (yellow): P_smooth ≥ θ1 = 0.7;

[0255] Orange: P_smooth ≥ θ2 = 0.8;

[0256] Red: P_smooth≥θ3=0.9.

[0257] The thresholds θ1–θ3 can be adjusted in real time via the ground software interface.

[0258] 2. Calculation of Physical Attenuation Coefficient Offset

[0259] Δλ=(λ*–λ_norm) / λ_norm

[0260] Δα=(α*–α_norm) / α_norm

[0261] Δβ=(β*–β_norm) / β_norm

[0262] Where λ_norm, α_norm, and β_norm are the "normal attenuation references" that are updated in real time in step S2.

[0263] 3. Early Warning-Decision Mapping Table (Editable on-site)

[0264]

[0265]

[0266] 4. Closed-loop control link

[0267] Early warning signal → SCADA system → PLC control of drilling rig parameters;

[0268] Simultaneously, Δλ, Δα, and Δβ are written into the WITS record for remote experts to review and analyze.

[0269] If P_smooth falls below the threshold within 60 seconds, the alarm will be automatically deactivated and drilling will resume.

[0270] 5. On-site verification

[0271] Three shale gas wells in southern Sichuan:

[0272] A 28-second orange warning indicates the drill is stuck. On-site, follow the instructions to reduce drill pressure and torque to prevent the drill from jamming.

[0273] A yellow warning for well leakage was issued 35 seconds in advance; after reducing the discharge rate by 8%, the leakage rate dropped to 2m³. 3 / h;

[0274] The early warning accuracy rate was 92%, and the false alarm rate was 4.3%, which meets the technical disclosure document's target of "improving early warning accuracy by 30% compared to traditional methods".

[0275] Through steps S6 / S7, the system completes the closed loop of "collection-optimization-early warning-intervention", realizing real-time, accurate and controllable management of complex drilling events.

[0276] Preferably, step S2 includes:

[0277] The physical feasible regions of λ, α, and β are updated in real time using the real-time mud viscosity μ(t) and well depth h(t) such that λ∈[λ_min(μ),λ_max(μ)], α∈[α_min(h),α_max(h)], and β∈[β_min(h),β_max(h)].

[0278] 1. Design Purpose

[0279] The physical meaning of the attenuation coefficient in this invention is closely related to real-time drilling condition data (formation density, mud viscosity, well depth, etc.) and must be dynamically updated rather than fixed. This section presents an online algorithm that refreshes the allowable ranges of λ, α, and β in real time with μ(t) and h(t), ensuring that the particle swarm search always operates within a range that is both engineering-safe and physically feasible.

[0280] 2. Real-time data interface

[0281] μ(t) — Outputted every 30s by an online rheometer (Fann35 / 50 model) at the wellhead, with a range of 5–120 mPa·s and an error of ±1%.

[0282] h(t) — measured in real time by a depth encoder, with an accuracy of ±0.5cm / 30m.

[0283] The two data streams enter the central processing unit via the WITS protocol, triggering the feasible domain recalculation thread (cycle 30s).

[0284] 3. Update of the feasible region of λ (for mud pressure decay)

[0285] 3.1 Theoretical Formula

[0286] λ is proportional to the friction coefficient f, and f and μ(t) satisfy the Darcy-Weisbach type empirical formula:

[0287] λ_theory=k1·μ(t)^0.25 / D_h^1.25

[0288] k1 is an empirical constant (given by the hydraulic calculation table of the well team); D_h is the equivalent diameter of the annulus (taken in steps with well depth).

[0289] 3.2 Boundary Scaling

[0290] λ_min(μ)=0.7·λ_theory

[0291] λ_max(μ) = 1.3·λ_theory

[0292] The coefficients 0.7–1.3 are derived from 200 sets of field calibration data, covering a 95% confidence interval.

[0293] 4. Update of the feasible region of α (for ROP decay)

[0294] 4.1 Segmented Experience Model

[0295] α is influenced by both well depth h and formation lithology. The target well section is divided into 100m segments, and a baseline value α_base(h) is given for each segment.

[0296] 4.2 Dynamic Boundary

[0297] α_min(h)=α_base(h)·(1–ε_α)

[0298] α_max(h)=α_base(h)·(1+ε_α)

[0299] ε_α = 0.25 (drill bit wear limit) + 0.05·(h–h_start) / 1000, which increases with well depth.

[0300] 5. Update the feasible region of β (for torque decay)

[0301] 5.1 Wellbore Cleanliness Correction

[0302] β is negatively correlated with the wellbore cleanliness index (CI). CI is calculated in real time based on cuttings concentration and annular return velocity.

[0303] CI(t) = (Amount of returned cuttings / Theoretical amount of cuttings) × 100%.

[0304] 5.2 Dynamic Boundary

[0305] β_min(h)=β_0·(1+0.2·CI(t))

[0306] β_max(h)=β_0·(1+0.8·CI(t))

[0307] β_0 is the baseline slope of the clean well section (kN·m·s) -1 ), obtained by fitting historical segments without anomalies.

[0308] 6. Update trigger and thread flow

[0309] Every 30 seconds: Read the latest μ(t), h(t), CI(t).

[0310] The feasible region is recalculated and broadcast to the particle swarm algorithm within 30ms.

[0311] If a particle crosses the new boundary, it executes a "bounce" strategy:

[0312] x_new = x_bound – (x_old – x_bound) to maintain search continuity.

[0313] 7. Examples of field parameters

[0314]

[0315] 8. Effect Verification

[0316] In a 4500m section of an ultra-deep well, after introducing dynamic boundaries, the feasible domain volume of the particle swarm was reduced by 48%, and the number of iterations decreased from an average of 118 to 72.

[0317] No out-of-bounds particles were observed, and there were no abnormal drilling stoppage records on site, which meets the requirement of "avoiding the generation of infeasible solutions" in the technical disclosure document.

[0318] Through the aforementioned real-time feasible domain update mechanism, λ, α, and β always fit the current mud properties, well depth, and wellbore cleanliness, providing a safe, compact, and physically interpretable search space for subsequent particle swarm optimization.

[0319] Preferably, the adaptive inertia weight ω_t in step S5 satisfies:

[0320] ω_t=ω_max-(ω_max-ω_min)·t / T_max, and ω_max=0.9, ω_min=0.4.

[0321] 1. Mathematical Expression

[0322] ω_t adopts a linear decreasing strategy:

[0323] ω_t=ω_max–(ω_max–ω_min)·t / T_max (Formula-A)

[0324] 2. Parameter values

[0325] ω_max = 0.9 (corresponding to "initial weight ω_max = 0.9" in the original disclosure document)

[0326] ω_min = 0.4 (corresponding to the original text in the disclosure document, "final weight ω_min = 0.4")

[0327] T_max is the maximum number of iterations, with a default value of 150 in practice (consistent with the implementation example in the manual).

[0328] 3. Physical meaning and stage division

[0329] Initially, t≈0: ω_t≈0.9 → Particles maintain a high velocity, which is beneficial for global exploration and prevents premature convergence.

[0330] Mid-term t≈T_max / 2: ω_t≈0.65 → Balances global and local search capabilities.

[0331] In the final stage, t≈T_max: ω_t≈0.4 → velocity decay, fine particle search, and improved accuracy of the optimal solution.

[0332] 4. Calculation Example

[0333]

[0334] 5. Mapping to on-site codes

[0335] In the GPU-CUDA kernel function, rewrite formula -A as follows:

[0336] __device__float omega(int t,int T_max){

[0337] return 0.9f - 0.5f * t / T_max;

[0338] }

[0339] The calculation is performed once before each iteration and broadcast to all particles, taking less than 1 μs.

[0340] 6. Effect Verification

[0341] Scene of 40 particles / 150 generations of shale gas wells in southern Sichuan:

[0342] After adopting linear ω_t, the number of times the local optimum was escaped increased from 2 to 7.

[0343] The final fitness F increased by 18%, which meets the goal of "balancing global search and local development" in the disclosure document.

[0344] Through the above-mentioned linear adaptive inertia weight design, the particle swarm algorithm achieves a smooth transition from wide-area search in the early stage to precise convergence in the later stage in the early warning scenario of complex drilling events. It not only preserves the global optimal probability, but also significantly improves the quality of the final solution.

[0345] Preferably, the Gaussian perturbation term ε in step S5 is ε ~ N(0,σ) 2 ), and σ=0.1·|gbest|.

[0346] 1. Mathematical definition

[0347] The Gaussian perturbation term ε follows a zero-mean normal distribution:

[0348] ε~N(0,σ 2 )

[0349] Where the standard deviation σ is a dynamic quantity:

[0350] σ = 0.1·|gbest|(Formula - B)

[0351] |gbest| represents the L2 norm (Euclidean length) of the current globally optimal particle, which is updated in real time and recalculated in each iteration.

[0352] 2. Introduction of Position and Function

[0353] The social cognitive component embedded in the particle swarm velocity update formula:

[0354] v_i^{t+1}=ω_t v_i^t+c1r1(pbest_i–x_i^t)

[0355] +c2r2(gbest–x_i^t+ε) (Formula-C)

[0356] It only applies to the gbest direction to avoid disrupting individual cognitive (pbest) information.

[0357] 3. The physical meaning of dynamic standard deviation σ

[0358] The larger |gbest| is, the larger σ is, which enhances the disturbance and can prevent the area from falling into a flat region;

[0359] The smaller the |gbest| value, the smaller the σ value, which reduces the perturbation and facilitates fine convergence.

[0360] This "adaptive noise intensity" is proportional to the drilling data noise level, which aligns with the disclosure document's objective of "reducing the impact of data noise".

[0361] 4. Calculation process

[0362] Step 1: Before each iteration, |gbest| is calculated by the host machine.

[0363] Step 2: Calculate σ using formula -B;

[0364] Step 3: GPU-CUDA parallel generation of ε:

[0365] __device__float epsilon(curandState*state,float gbest_norm){

[0366] return 0.1f*gbest_norm*curand_normal(state); / / N(0,1)→N(0,σ 2 )}

[0367] Step 4 adds ε to the social term of the formula, taking <1μs / particle.

[0368] 5. Boundary and Stability

[0369] ε only affects the direction of velocity and does not directly modify position;

[0370] The speed limit remains |v_i,d|≤v_max,d to prevent out-of-bounds errors due to large disturbances.

[0371] 6. On-site verification

[0372] Shale gas well in southern Sichuan (40 particles, 150 generations)

[0373] Without ε: It escapes the local optimum twice, and the final F value is 7.82;

[0374] Adding ε: escaped local optima 7 times, and the final F value was 9.23, an improvement of 18%;

[0375] The false negative rate decreased from 6.8% to 4.3%, and the false positive rate decreased from 5.1% to 3.0%.

[0376] 7. Correspondence with the disclosure document

[0377] Original text on "Perturbation-based speed updates":

[0378] "Where c1 = c2 = 2, r1, r2 ~ U(0, 1), ε ~ N(0, 0. 1)"

[0379] In this embodiment, "0.1" is upgraded to "0.1·|gbest|", which makes the perturbation intensity adaptive with the solution space scale, further enhancing the noise resistance and global escape capability.

[0380] Through the aforementioned Gaussian perturbation design, the particle swarm optimization algorithm maintains the diversity of global search when facing high-dimensional and noisy drilling data, and automatically reduces perturbation when approaching the optimal solution, thus achieving robust and efficient parameter optimization.

[0381] Preferably, the penalty term in step S4 is:

[0382] Penalty=∑_{k∈{λ,α,β}}max(0,(k-k_bound) 2 / σ_k 2 ),

[0383] k_bound is the engineering constraint boundary with the corresponding coefficient, and σ_k is the tolerance coefficient.

[0384] 1. Design Purpose

[0385] This technique emphasizes that during particle swarm optimization, it is essential to "embed parametric physical constraints to avoid generating infeasible solutions." This penalty term explicitly writes the engineering safety boundaries of λ, α, and β into the fitness function, immediately penalizing any out-of-bounds particle, thus ensuring that the optimization result always remains within the field-executable range.

[0386] 2. Mathematical Expressions

[0387] Penalty=Σ_{k∈{λ,α,β}}w_k·max(0,(k–k_bound) 2 / σ_k 2 (Formula-D)

[0388] Variable description:

[0389] k — the value of the current particle in dimension k (k is λ, α or β);

[0390] k_bound — the engineering constraint boundary of this dimension, calculated in real time on site;

[0391] σ_k — Tolerance coefficient, reflecting the allowable fluctuation width at the boundary;

[0392] w_k — Penalty weight (default w_λ=w_α=w_β=1, which can be adjusted in real time in the ground software).

[0393] 3. Determination of Boundaries and Tolerance Coefficients

[0394]

[0395] 4. Real-time calculation process

[0396] Step 1: Read μ(t), h(t), CI(t) every 30 seconds and recalculate λ_bound, α_bound, and β_bound.

[0397] Step 2: In each iteration of the particle swarm, substitute the k value of each particle into the formula -D;

[0398] Step 3: If k is within the feasible region, the max term is 0, and there is no penalty; if it exceeds the boundary, the penalty value is proportional to the square of the distance beyond the boundary.

[0399] 5. Adjustment of punishment intensity

[0400] The penalty weight w_k can be adjusted in real time through the ground software interface:

[0401] During the leak-sensitive phase, w_λ↑2× immediately removes particles from the high-pressure danger zone.

[0402] During the normal drilling phase, w_k = 1, maintaining the degree of freedom in the search.

[0403] 6. On-site verification

[0404] Examples of shale gas wells in southern Sichuan:

[0405]

[0406] The results show that:

[0407] The penalty reduces the proportion of out-of-bounds particles to 0%;

[0408] Fitness decreased only slightly (<8%), and the search remained efficient within the feasible domain;

[0409] There were no records of abnormal drilling stoppages due to parameters exceeding the limits, which meets the requirement of "high engineering practicality".

[0410] 7. Code snippet (GPU-CUDA kernel function)

[0411]

[0412] The process is incremented dimension by dimension during invocation, with a total time consumption of <1μs / particle.

[0413] Through the above penalty term design, particle swarm optimization can ensure prediction accuracy and always meet engineering safety and process feasibility constraints in the early warning scenario of complex drilling events.

[0414] Preferably, the multi-level warning signals in step S7 include:

[0415] Yellow alert: Δλ / λ_normal > 1.2 and P_event ≥ 0.7;

[0416] Orange alert: Δα / α_normal > 1.3 and P_event ≥ 0.8;

[0417] Red alert: Δβ / β_normal < 0.5 and torque surge > 20% and P_event ≥ 0.9.

[0418] 1. Design Purpose

[0419] This technical requirement mandates that "when P ≥ θ, a three-level early warning should be sent via an audible and visual alarm system and a mobile app." This embodiment maps three typical complex events (well leakage, sudden drop in drill bit wear / ROP, and signs of stuck drill bit) to real-time offsets Δλ, Δα, and Δβ of λ, α, and β, respectively, forming a physically interpretable and manageable "yellow-orange-red" three-level early warning logic.

[0420] 2. Definition and Triggering Conditions of Warning Levels

[0421]

[0422] illustrate:

[0423] λ_normal, α_normal, and β_normal are refreshed every 30 seconds in step S2;

[0424] Δλ, Δα, and Δβ are calculated in real time by step S6;

[0425] A sudden torque surge is defined as a torque increase of more than 20% of the baseline value within a 5-second sliding window.

[0426] 3. Basis for threshold setting

[0427] Yellow: Δλ>1.2 corresponds to a pressure decay rate >1.2 times the normal rate, with an empirical lower limit of 1.15, leaving a 5% safety margin;

[0428] Orange: Δα>1.3 corresponds to ROP attenuation > 30% of normal, which is consistent with the drill bit wear test curve;

[0429] Red: Δβ<0.5 corresponds to the torque "gradually decreases then suddenly increases" window period, which is 20-40 seconds earlier according to on-site statistics.

[0430] 4. Real-time calculation and judgment process

[0431] ① Read the latest λ*, α*, β* and real-time torque values ​​every 1 second;

[0432] ② Calculate Δλ, Δα, Δβ and the 5s slope of the torque;

[0433] ③ Determine the level according to the conditions in the table above;

[0434] ④ If multiple conditions are met simultaneously, the highest level shall prevail to avoid duplicate alarms.

[0435] 5. On-site output and intervention actions

[0436]

[0437] 6. Adjustable threshold interface

[0438] Ground monitoring software provides sliders:

[0439] θ_yellow∈[0.6,0.8], default 0.7;

[0440] θ_orange∈[0.7,0.9], default 0.8;

[0441] θ_red∈[0.8,1.0], default 0.9;

[0442] Simultaneously, the coefficients of Δλ, Δα, and Δβ can be modified to achieve "one well, one policy".

[0443] 7. On-site verification

[0444] Three shale gas wells in southern Sichuan:

[0445] Yellow alerts were triggered 12 times, well leakage was confirmed on site 11 times, and there was 1 false alarm;

[0446] Six orange alerts were triggered, indicating excessive drill bit wear was detected 25–35 seconds in advance.

[0447] The red alert was triggered twice, with the drill getting stuck 28 seconds and 33 seconds earlier each time. The drill was successfully pulled out on site without any jamming incidents.

[0448] Through the aforementioned three-level early warning system based on the dual conditions of "physical offset + probability threshold", a closed-loop management of "monitoring-early warning-intervention" was achieved on-site, significantly reducing the economic losses and safety risks caused by complex events.

[0449] Preferably, the real-time early warning model M is a support vector machine or a deep neural network, and the parameters of its kernel function or activation function are directly assigned by gbest.

[0450] 1. Model Selection

[0451] M adopts the following two interchangeable structures, which can be switched on-site with a single button via software switch:

[0452] Support Vector Machine (SVM) with RBF kernel – suitable for small-scale well drilling scenarios with fewer than 50,000 samples and a dimension of ≤20;

[0453] Deep Neural Networks (DNN) (3 fully connected layers + ReLU) – suitable for scenarios with large amounts of data (>100,000 records) or high online self-learning requirements.

[0454] 2. Mapping relationship between hyperparameters and gbest

[0455] After linear / logarithmic scaling, gbest = [λ*, α*, β*, par1, par2, ...] is directly mapped to the core hyperparameters of M. The mapping table is as follows:

[0456]

[0457] All mappings are monotonic functions, ensuring that the direction of gbest change is consistent with the trend of model performance.

[0458] If a gbest element exceeds [0,1], it is first clipped to the boundary and then mapped to prevent extreme values.

[0459] 3. Real-time assignment process

[0460] ① After the particle swarm converges, the mapping calculation is completed on the host side within 1ms;

[0461] ② Write hyperparameters into the SVM / DNN instance via memory pointers;

[0462] ③ For DNN, if the number of hidden nodes changes, the pre-compiled CUDA dynamic parallel kernel function is used to reallocate GPU memory in real time, taking less than 3ms.

[0463] ④ Once mapping is complete, the online inference thread is entered immediately, with an overall switching delay of <5ms, which does not affect the 200ms real-time link.

[0464] 4. On-site verification

[0465] Actual measurements of shale gas wells in southern Sichuan:

[0466] gbest=[0.63,0.42,0.58,0.31]→DNN(hidden1=31, hidden2=14, dropout=0.158, η=0.0038)

[0467] The validation set accuracy was 92.4%, a 12% improvement over the fixed hyperparameter baseline.

[0468] If SVM is used instead, the same gbest mapping yields C=42.7, γ=0.84, and an accuracy of 90.1%, satisfying the requirement of "model interchangeability".

[0469] 5. Adaptive retraining

[0470] For each well completed, the latest fault sample is added to the training set;

[0471] Sample size < 50,000: Trigger SVM retraining, taking 2–3 seconds;

[0472] Sample size ≥ 50,000: Trigger incremental DNN training (Adam, 3 epochs), taking 8–10 seconds;

[0473] After retraining, gbest re-participates in a new round of optimization, achieving a closed loop of "well-to-well self-learning".

[0474] The "gbest→hyperparameter" direct assignment mechanism eliminates the need for manual parameter tuning on-site, maintaining real-time performance while ensuring that model performance adapts and improves with changing working conditions.

[0475] Preferably, the sampling frequency of step S1 is ≥10Hz, and the data preprocessing includes 3σ anomaly removal and linear interpolation completion.

[0476] 1. Sampling frequency setting

[0477] Fixed master clock: 10Hz (100ms cycle) to trigger all sensors uniformly;

[0478] Clock source: PTP timing module of downhole drilling system, drift <1ms / 24h;

[0479] Redundancy: If a channel loses one sampling point, the system marks it as "missing" instead of filling it in on average, ensuring that subsequent 3σ determinations are true and effective.

[0480] 2. Data Channel and Raw Frame Format

[0481] One frame vector is generated every 100ms.

[0482] x_raw(t)=[P_mud,ROP,T,WOB,RPM,ρ_mud,V_pit] T ,

[0483] It has 7 dimensions, each with 32-bit float, a frame size of 28 bytes, and a network bandwidth requirement of 280B / s, which is far lower than the 9600bps margin of the WITS link.

[0484] 3. 3σ anomaly removal

[0485] 3.1 Sliding Window

[0486] Window length 60s (600 points), step size 1 point;

[0487] Calculate the mean μ_i and standard deviation σ_i independently for each dimension.

[0488] 3.2 Elimination Rules

[0489] If |x_raw,i(t)–μ_i|>3σ_i, mark it as an anomaly and set it to NaN;

[0490] Rejection rate control: Single channel rejection rate <0.3%, total rejection rate <2%.

[0491] 3.3 Real-time implementation

[0492] Using GPU shared memory to compute 600 points μ_i and σ_i in parallel, the time taken is <0.5ms;

[0493] An exception is immediately written to the ring-buffer after being marked for use in the next interpolation step.

[0494] 4. Linear interpolation completion

[0495] 4.1 Interpolation Object

[0496] Interpolation is performed only on points marked as NaN with 3σ to avoid over-smoothing.

[0497] 4.2 Interpolation Method

[0498] If there are ≤5 consecutive missing points: use linear interpolation with the nearest valid points before and after;

[0499] If more than 5 consecutive points are missing: retain NaN, trigger "sensor offline" alarm, and do not force completion.

[0500] 4.3 Interpolation Formula

[0501] x_interp(t)=x_left+(x_right–x_left)·(t–t_left) / (t_right–t_left).

[0502] 5. Delay and Accuracy

[0503] 3σ calculation + elimination delay: <1ms per 600 points;

[0504] Linear interpolation delay: <0.05ms per missing point;

[0505] The overall latency is <2ms, which meets the requirement that the total link latency is <200ms.

[0506] 6. On-site verification

[0507] Actual measurements at three wells:

[0508] Raw 10Hz data 1.2×10 8 point;

[0509] 3σ anomalies were eliminated at a rate of 0.21%.

[0510] Linear interpolation completion is 0.18%;

[0511] After processing, the data integrity is 99.61%, with no additional phase distortion, meeting the requirements of subsequent physical decay models and particle swarm optimization.

[0512] By combining 10Hz high-frequency sampling with 3σ anomaly removal and linear interpolation completion, the system effectively removes sensor spikes and transient interruptions while ensuring real-time performance, providing clean and continuous input data for high-precision early warning.

[0513] Preferably, the online inference latency of step S6 is <200ms, and it is achieved through parallel GPU acceleration.

[0514] 1. Delayed Target

[0515] This technology explicitly requires that "single-well data processing latency < 200ms". This section explains how to compress the entire process (data acquisition → preprocessing → model inference → probability output) to < 200ms and implement it in parallel using GPU.

[0516] 2. Link decomposition and time budgeting

[0517]

[0518] 3. GPU parallel architecture

[0519] Hardware: NVIDIA Jetson Xavier NX (6-core ARM + 384 CUDA) or PCIe RTX3060 Laptop;

[0520] Software: CUDA 11.7 + cuDNN 8.5 + TensorRT 8.4;

[0521] Model format: ONNX → TensorRT engine (FP16, maximum batch size 1).

[0522] 4. Key Parallel Implementation

[0523] 4.1 Zero-copy data-model

[0524] Sensor DMA writes the raw frame to the GPU global memory;

[0525] The preprocessing kernel performs 3σ, interpolation, and normalization within the same memory block, avoiding... copy.

[0526] 4.2 Inferring the kernel

[0527] Using the TensorRT auto-tuned engine, batch=1, input 7×1 vector;

[0528] The output P_event∈[0,1] is written directly back to the video memory for the next smoothing step.

[0529] 4.3 Probabilistic Smoothing Kernel

[0530] 5. Point moving average:

[0531] shared float buf[5];

[0532] __shfl_sync enables register swapping within a warp, which takes 1ms to complete.

[0533] 4.4 Early warning judgment kernel

[0534] Atomic operations update the global flag; the CPU polls the flag, triggering an audio-visual / PLC signal within 0.5ms.

[0535] 5. End-to-end measurement

[0536] Three shale gas wells in southern Sichuan have been operating continuously for 72 hours: the average end-to-end delay is 187ms (maximum 199ms);

[0537] GPU utilization was 38%, power consumption was 15W, and temperature rise was <55℃; there was no packet loss and no timeout.

[0538] 6. Expansion and Fault Tolerance

[0539] If the GPU fails, it will automatically slow down to the CPU path (latency 400–500ms) and report a "performance degradation" alarm; it supports multi-well parallelism: the same GPU can serve 2–3 wells at the same time, and the latency is still <200ms / well.

[0540] Through the aforementioned GPU parallel acceleration link, this system compresses the entire process of "data acquisition → early warning output" to within 200ms, meeting the real-time requirements of second-level early warning at the drilling site.

[0541] Furthermore, a computer-readable storage medium having a program stored thereon, which, when executed by a processor, implements the steps of the above method.

[0542] By adopting the above-disclosed technical solution of this invention, the following beneficial effects are obtained:

[0543] 1. Economic benefits: On average, well leakage or stuck drill bit is detected 25-40 seconds earlier per well, increasing the success rate of on-site intervention by 50%; According to statistics of shale gas wells in southern Sichuan, each avoided stuck drill bit accident can save 1.8-2.2 million yuan (including retrieval, downtime, and mud loss); The entire well team reduces non-production time by an average of 12% per year, directly improving the profitability of a single well.

[0544] 2. Safety benefits: The false alarm rate has been reduced from 15% to 4%, and the false alarm rate has been reduced from 12% to 3%; there are zero records of fatal accidents caused by stuck drill pipe on site, meeting Sinopec's "zero fatalities in complex events" red line requirement.

[0545] 3. Intelligent benefits: No manual parameter adjustment is required on site: gbest → automatic hyperparameter mapping, the model adapts to working conditions; supports OTA upgrades and multi-well parallel operation, one algorithm covers all scenarios of deep wells, shale gas and horizontal wells.

[0546] 4. Real-time performance benefits: GPU parallel link end-to-end latency <200ms; it can still maintain <250ms in a 4G public network environment, meeting the "second-level early warning" industrial standard.

[0547] In summary, this invention has significant advantages in three aspects: real-time performance, accuracy, and engineering feasibility, providing a replicable technical paradigm for safe and efficient drilling in deep wells, ultra-deep wells, and complex formations.

[0548] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A real-time early warning method for complex drilling events based on physical decay-particle swarm optimization, characterized in that, Includes the following steps: Step S1: During the drilling process, D-dimensional drilling parameters are collected in a fixed time sequence to form a time series dataset X = {x(1), x(2), ..., x(T)}, where x(t) ∈ R^D includes at least mud pressure, mechanical drilling rate ROP, and torque T; Step S2: Construct three types of physical attenuation benchmark models respectively: S21, Mud pressure decay model: P(t) = P0·e^(-λt) + P_min, S22, Mechanical drilling rate attenuation model ROP(h)=ROP0·(1-α·h), S23, Torque decay model T(t)=T0-β·t+δ(t), λ, α, and β are then used as attenuation coefficients to be optimized. Step S3: Initialize N particles using the attenuation coefficients λ, α, β and D-3 auxiliary parameters as particle dimensions. The position vector of each particle represents a set of parameter combinations to be optimized. Step S4: Establish the coupling fitness function F, which is composed of the following weighted sub-functions: F1 = 1 / MSE1, where MSE1 is the mean square error between the predicted and measured values ​​of mud pressure. F2 = 1 / MSE2, where MSE2 is the mean square error between the ROP predicted value and the measured value; F3 = 1 / MSE3, where MSE3 is the mean square error between the predicted torque value and the measured torque value; F4 = Penalty terms for λ, α, and β within the preset engineering safety range; Step S5: Utilize the adaptive inertia weight ω t An improved particle swarm optimization algorithm with Gaussian perturbation term ε iteratively updates particle positions to maximize F and obtain the globally optimal particle gbest. Step S6: Map gbest to the hyperparameters of the real-time early warning model M, perform online inference of the drilling parameters at the current time t, and output the event probability P_event∈[0,1]; Step S7: When P_event≥θ, trigger a multi-level early warning signal and provide the physical attenuation coefficient offset Δλ, Δα or Δβ for the corresponding working condition to guide the adjustment of drilling parameters.

2. The method according to claim 1, wherein step S2 comprises: The physical feasible regions of λ, α, and β are updated in real time using the real-time mud viscosity μ(t) and well depth h(t) such that λ∈[λ_min(μ),λ_max(μ)], α∈[α_min(h),α_max(h)], and β∈[β_min(h),β_max(h)].

3. The method according to claim 2, wherein the adaptive inertia weight ω_t in step S5 satisfies: ω_t=ω_max-(ω_max-ω_min)·t / T_max, and ω_max=0.9, ω_min=0.

4.

4. The method according to claim 3, wherein the Gaussian perturbation term ε ~ N(0,σ) in step S5 2 ), and σ=0.1·|gbest|.

5. The method according to claim 4, wherein the penalty term in step S4 is: Penalty=∑_{k∈{λ,α,β}}max(0,(k-k_bound) 2 / s_k 2 ), k_bound is the engineering constraint boundary with the corresponding coefficient, and σ_k is the tolerance coefficient.

6. The method according to claim 5, wherein the multi-level early warning signal in step S7 includes: Yellow alert: Δλ / λ_normal > 1.2 and P_event ≥ 0.7; Orange alert: Δα / α_normal > 1.3 and P_event ≥ 0.8; Red alert: Δβ / β_normal < 0.5 and torque surge > 20% and P_event ≥ 0.

9.

7. The method according to claim 6, wherein the real-time early warning model M is a support vector machine or a deep neural network, and the parameters of its kernel function or activation function are directly assigned by gbest.

8. The method according to claim 7, wherein the sampling frequency of step S1 is ≥10Hz, and the data preprocessing includes 3σ anomaly removal and linear interpolation completion.

9. The method according to any one of claims 1-8, wherein the online inference latency of step S6 is <200ms and is accelerated by parallel GPU.

10. A computer-readable storage medium having a program stored thereon, which, when executed by a processor, implements the steps of the method according to any one of claims 1-8.

Images