A method, device, and equipment for real-time monitoring and intelligent correction control of borehole trajectory
By using multi-sensor fusion and chaotic control theory, the borehole trajectory deviation vector is decomposed and energy is managed, achieving high-precision adaptive intelligent deviation control of the borehole trajectory. This solves the problems of insufficient accuracy and low energy utilization efficiency of traditional borehole trajectory control methods, and improves the success rate and economic benefits of borehole operations.
Patent Information
- Application Number
- CN202511408820.9
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2025-09-29
- Publication Date
- 2026-01-06
- Estimated Expiration
- 2045-09-29
AI Technical Summary
Traditional borehole trajectory control methods are insufficient in accuracy and adaptability in complex geological environments, have low energy utilization efficiency, and are difficult to achieve fine control, resulting in decreased borehole accuracy and project failure.
By employing multi-sensor fusion measurement technology, the drill bit attitude is generated after compensation through distortion analysis of data from inertial navigation and geomagnetic sensors. The trajectory deviation vector is decomposed into positive and negative deviation components, and a deviation energy pool is constructed. Using chaotic control theory and energy processing, combined with hydraulic cylinder control, precise thrust distribution and deflection torque are generated to achieve intelligent deviation correction control.
It achieves high-precision, adaptive, and intelligent control of the drilling trajectory, improves the accuracy and stability of the deviation control, avoids instability caused by blind control, and enhances the success rate and economic benefits of drilling operations.
Smart Images

Figure CN120867709B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of intelligent technology in borehole engineering, and in particular to a method, device, and equipment for real-time monitoring and intelligent correction control of borehole trajectory. Background Technology
[0002] Drilling trajectory control is a key technology in modern drilling engineering, directly affecting drilling quality and construction efficiency. In complex geological environments, the drill bit is prone to deviating from the predetermined trajectory, leading to decreased drilling accuracy and even project failure. Current trajectory correction technologies mainly rely on manual experience and basic control methods, which have significant shortcomings in terms of accuracy and automation.
[0003] Traditional borehole trajectory monitoring typically uses a single sensor, which is easily affected by the complex downhole environment and results in significant measurement errors. Conventional deviation control methods often employ simple feedback mechanisms, lacking in-depth analysis of the complex drilling dynamics and thus offering limited control effectiveness. Existing energy management methods are relatively crude, unable to achieve precise energy allocation and utilization. Furthermore, traditional technologies lack sufficient understanding of the nonlinear characteristics and fluctuations in the drilling process, making it difficult to handle control challenges under complex conditions. These limitations become more pronounced as borehole depth increases or complex formations are encountered, severely impacting the success rate and economic efficiency of drilling operations. Therefore, there is an urgent need to develop more advanced intelligent borehole trajectory control technologies. Summary of the Invention
[0004] This invention provides a method, device, and equipment for real-time monitoring and intelligent correction control of borehole trajectory, aiming to solve the technical problems of insufficient accuracy, poor adaptability, and low energy utilization efficiency of traditional borehole trajectory control methods. Through multi-sensor fusion measurement, trajectory deviation energy processing, application of chaotic control theory, and intelligent correction control, it achieves high-precision, adaptive, and intelligent control of borehole trajectory.
[0005] The first aspect of this invention proposes a method for real-time monitoring and intelligent correction control of borehole trajectory, comprising the following steps:
[0006] Inertial navigation data and geomagnetic sensor data are collected during the drilling process. Distortion analysis is performed on the geomagnetic sensor data to obtain magnetic field disturbance characteristics. Waveform compensation and fusion are then performed using the magnetic field disturbance characteristics to generate the compensated drill bit attitude.
[0007] Based on the compensated drill bit attitude, the trajectory deviation vector is extracted, and the trajectory deviation vector is decomposed into positive deviation components and negative deviation components. Potential energy is accumulated in the negative deviation components to generate a deviation energy pool.
[0008] The deviation energy pool is charged and discharged to identify stable and unstable boundaries. A chaotic control band is determined between the stable and unstable boundaries, and the optimal response range is extracted through the chaotic control band.
[0009] Based on the optimal response range and the positive deviation component, an active correction vector is determined, the deviation energy pool is controlled to release energy to obtain the energy release amount, and the active correction vector and the energy release amount are vector-superimposed to generate a composite correction command.
[0010] The composite correction command controls the hydraulic cylinder to generate thrust distribution, and the thrust distribution acts on the drill bit to generate deflection torque. The deflection torque is used to change the drilling direction to obtain the actual correction amount, and the actual correction amount is subjected to oscillation feature extraction to form system oscillation parameters.
[0011] A stability index is generated by mapping the system oscillation parameters with the chaotic control zone. A fine correction force is formed by superimposing the stability index with the thrust distribution. A chaotic edge control field is constructed based on the fine correction force and the negative deviation component.
[0012] The conversion efficiency between the deviation energy pool and the actual correction amount is obtained, and a real-time correction control amount is generated based on the chaotic edge control field and the conversion efficiency.
[0013] A second aspect of this invention provides a real-time monitoring and intelligent correction control device for borehole trajectory, comprising:
[0014] The data fusion module is used to collect inertial navigation data and geomagnetic sensor data during the drilling process, perform distortion analysis on the geomagnetic sensor data to obtain magnetic field disturbance characteristics, and use the magnetic field disturbance characteristics to perform waveform compensation fusion to generate the compensated drill bit attitude.
[0015] The energy conversion module is used to extract the trajectory deviation vector based on the compensated drill bit attitude, decompose the trajectory deviation vector into positive deviation components and negative deviation components, and accumulate the potential energy of the negative deviation components to generate a deviation energy pool.
[0016] The chaos analysis module is used to monitor the charging and discharging of the deviation energy pool to identify stable and unstable boundaries, determine a chaotic control band between the stable and unstable boundaries, and extract the optimal response range through the chaotic control band.
[0017] The instruction generation module is used to determine an active correction vector based on the optimal response range and the positive deviation component, control the release of the deviation energy pool to obtain the energy release amount, and vector superimpose the active correction vector and the energy release amount to generate a composite correction instruction.
[0018] The execution control module is used to control the hydraulic cylinder to generate thrust distribution through the composite correction command, generate deflection torque based on the thrust distribution acting on the drill bit, change the drilling direction using the deflection torque to obtain the actual correction amount, and extract the oscillation features of the actual correction amount to form system oscillation parameters.
[0019] The stability adjustment module is used to generate a stability index by mapping the system oscillation parameters with the chaotic control band, to form a fine correction force by superimposing the stability index with the thrust distribution, and to construct a chaotic edge control field based on the fine correction force and the negative deviation component.
[0020] An optimized control module is used to obtain the conversion efficiency between the deviation energy pool and the actual correction amount, and to generate a real-time correction control amount based on the chaotic edge control field and the conversion efficiency.
[0021] A third aspect of the present invention provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of a real-time monitoring and intelligent correction control method for borehole trajectory disclosed in the first aspect.
[0022] The beneficial effects of this invention are reflected in the following aspects: First, a complete drill bit attitude measurement and trajectory deviation analysis system is established. High-precision attitude data is acquired through multi-sensor fusion, and the trajectory deviation vector is innovatively decomposed into positive and negative deviation components. The negative deviation component generates a deviation energy pool through potential energy accumulation, realizing energy management of the correction process. Second, a chaotic control zone theoretical framework and optimal response interval extraction technology are constructed. By monitoring the charging and discharging of the deviation energy pool to form an energy fluctuation curve, the stable boundary and instability boundary of the system are accurately identified. Based on this, the chaotic control zone is determined and the optimal response interval is extracted, enabling the correction control to be carried out within the optimal working range, avoiding instability problems caused by blind control. Finally, a complete control link from active correction vector generation to real-time control execution is realized. Precise thrust distribution and deflection torque are generated through hydraulic cylinder control. A chaotic edge control field is constructed by combining system oscillation parameter analysis, and real-time correction control quantity is generated by obtaining energy conversion efficiency, forming a closed-loop intelligent control loop. This allows the entire correction process to be dynamically adjusted according to the actual drilling state, improving the accuracy and stability of the correction control.
[0023] It should be understood that the above general description and the following detailed description are exemplary and explanatory only, and do not limit this application. Attached Figure Description
[0024] The accompanying drawings illustrate specific examples of the technical solutions described in this invention and, together with the detailed embodiments, form part of the specification, serving to explain the technical solutions, principles, and effects of this invention.
[0025] Unless otherwise specified or defined, the same reference numerals in different figures represent the same or similar technical features, and different reference numerals may be used to represent the same or similar technical features.
[0026] Figure 1 This is a flowchart illustrating a real-time monitoring and intelligent correction control method for borehole trajectory according to the present invention.
[0027] Figure 2 This is a structural block diagram of a drilling trajectory real-time monitoring and intelligent correction control device according to the present invention.
[0028] Figure 3 This is a schematic diagram of the structure of a computer device according to the present invention. Detailed Implementation
[0029] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of the embodiments. Based on the embodiments of this application, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of this application.
[0030] It should be noted that all directional indicators (such as up, down, left, right, front, back, etc.) in the embodiments of this application are only used to explain the relative positional relationship and movement of each component in a certain specific posture (as shown in the figure). If the specific posture changes, the directional indicator will also change accordingly.
[0031] It should also be noted that when a component is described as "fixed to" or "set on" another component, it can be directly on the other component or there may be an intervening component present. When a component is described as "connected to" another component, it can be directly connected to the other component or there may be an intervening component present.
[0032] Furthermore, the use of terms such as "first" and "second" in this application is for descriptive purposes only and should not be construed as indicating or implying their relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined as "first" or "second" may explicitly or implicitly include at least one of that feature. Additionally, the technical solutions of the various embodiments can be combined with each other, but only on the basis of being achievable by those skilled in the art. When the combination of technical solutions is contradictory or impossible to implement, such a combination of technical solutions should be considered non-existent and not within the scope of protection claimed in this application.
[0033] The technical solutions of the embodiments of this application will be described below.
[0034] like Figure 1 As shown, this embodiment of the invention provides a method for real-time monitoring and intelligent correction control of borehole trajectory, including the following steps S110-S170:
[0035] Step S110: Collect inertial navigation data and geomagnetic sensor data during the drilling process, perform distortion analysis on the geomagnetic sensor data to obtain magnetic field disturbance characteristics, and use the magnetic field disturbance characteristics to perform waveform compensation and fusion to generate the compensated drill bit attitude.
[0036] Specifically, multi-dimensional data acquisition during the drilling process is achieved through a high-precision inertial measurement unit (IMU) and a three-axis geomagnetic sensor array. Inertial navigation data acquisition utilizes a combination of MEMS gyroscopes and accelerometers. The gyroscopes have a measurement range of ±2000° / s and a resolution of 0.01° / s, with a sampling frequency set to 200Hz to capture rapid attitude changes of the drill bit. The accelerometers have a measurement range of ±16g and a sensitivity of 4096 LSB / g, simultaneously acquiring three-axis linear acceleration. Data preprocessing uses a Kalman filter to eliminate random noise. The filter's state equation includes four state variables: position, velocity, attitude angle, and angular velocity. The observation equation incorporates GPS positioning information for position correction. Geomagnetic sensor data acquisition employs fluxgate sensors with a measurement accuracy of ±0.5nT and an operating frequency of 1kHz. Sensor nodes are deployed at the drill bit, drill rod, and drill platform to form a three-point magnetic field monitoring network. The data synchronization mechanism uses hardware timestamps to ensure that the time alignment error between inertial data and magnetic field data is less than 1ms. A coordinate system was established with the borehole starting point as the origin, the vertically downward direction as the positive Z-axis, the north as the X-axis, and the east as the Y-axis—a right-handed coordinate system. A temperature compensation algorithm was used to correct sensor drift in real time based on the temperature gradient at the borehole depth (average 3°C / 100m). Multi-sensor collaborative acquisition and precise time synchronization resulted in inertial navigation data streams and geomagnetic sensor data streams throughout the drilling process.
[0037] Distortion analysis was performed on geomagnetic sensor data to obtain magnetic field disturbance characteristics. The distortion analysis employed a frequency domain decomposition method, decomposing the raw magnetic field data into four frequency bands using Fast Fourier Transform (FFT): DC component, low-frequency component (<1Hz), mid-frequency component (1-10Hz), and high-frequency component (>10Hz). The DC component reflects the static characteristics of the Earth's magnetic field; the low-frequency component mainly originates from drill bit rotation and geological structural changes; the mid-frequency component corresponds to drill bit vibration and hydraulic system pulsation; and the high-frequency component includes electrical equipment interference and sensor noise. Anomaly detection algorithms, based on statistical methods, calculated the mean μ and standard deviation σ for each frequency band, marking data points exceeding the μ±3σ range as outliers. Interference source classification used a machine learning clustering algorithm, categorizing anomalous magnetic field patterns into three main types: steel casing interference, motor magnetic field, and formation magnetic minerals. The magnetic field gradient tensor was calculated by differencing adjacent sensor data, with a gradient threshold set to 500 nT / m; regions exceeding this threshold were defined as strong disturbance zones. Time-series correlation analysis revealed a significant correlation between magnetic field disturbances and drilling parameters (rotation speed, drilling speed, and pump pressure). Parameter combinations with a correlation coefficient r > 0.7 were identified as the main sources of disturbance.
[0038] In some embodiments, the step of using the magnetic field disturbance features to perform waveform compensation and fusion to generate the compensated drill bit attitude includes: constructing a disturbance spatiotemporal map based on the magnetic field disturbance features; extracting the disturbance propagation path from the disturbance spatiotemporal map; constructing a compensation waveform in reverse along the disturbance propagation path; and fusing the compensation waveform with the inertial navigation data to generate the compensated drill bit attitude.
[0039] A spatiotemporal map of disturbances was constructed based on magnetic field disturbance characteristics. A three-dimensional tensor representation T(x, y, t) was used, where x and y are spatial coordinates, t is the time coordinate, and the tensor element values are the magnetic field disturbance intensity at the corresponding spatiotemporal point. The spatial dimension was based on the geometric configuration of a three-point magnetic field monitoring network, using the sensor positions of the drill bit, drill pipe, and drill platform as reference points, and triangulation to determine the spatial distribution of the disturbance source. The temporal dimension used the marked anomalous moments in the disturbance characteristics as key nodes, and B-spline interpolation was used to construct a continuous temporal evolution curve. Disturbance intensity quantification employed normalization, mapping the three identified disturbance categories (steel casing interference, motor magnetic field, and formation magnetic minerals) uniformly to the [0, 1] interval. The map resolution was set to 10cm × 10cm spatially and 0.1s temporally to ensure the capture of detailed disturbance features. A data fusion algorithm weighted averaged the disturbance information from multiple sensors, with weight coefficients determined based on the distance from the sensor to the disturbance source; strong disturbance areas with gradient tensor values greater than 500 nT / m received higher weights. The smoothing of the graph uses a Gaussian convolution kernel with a convolution radius of 2cm to eliminate the influence of local noise.
[0040] Perturbation propagation paths are extracted from the perturbation spatiotemporal map. A gradient pursuit algorithm is used to calculate the gradient vector at each spatiotemporal point in the map. ,in and These are the partial derivatives of the disturbance intensity T with respect to the spatial coordinates x and y, respectively. The gradient is the partial derivative of the perturbation intensity with respect to time t, and the gradient direction indicates the direction of perturbation energy propagation. Path origin identification is performed by finding local maxima in the graph, i.e., points that satisfy... Furthermore, points with negative definite Hessian matrices correspond to the locations of perturbation sources. Path tracing starts from the starting point and iteratively calculates along the gradient descent direction with a step size of 0.5 cm until the gradient magnitude is less than the threshold of 0.01 or the map boundary is reached. Path branching uses the watershed algorithm; when a starting point corresponds to multiple paths, the main path is selected based on path length and energy decay rate. Propagation speed is calculated based on the distance and time difference between adjacent spatiotemporal points, v = Δs / Δt, where Δs is the spatial distance and Δt is the time interval. Path clustering analysis groups similar propagation paths into three categories: fast decay path groups, medium decay path groups, and slow decay path groups. Each path node records its spatial coordinates (x, y, z) and timestamp t.
[0041] A compensation waveform is constructed in reverse along the disturbance propagation path. Based on the principle of reverse propagation of disturbance energy, an anti-phase compensation signal is calculated at each node coordinate (x, y, z) along the path according to the propagation velocity v and the attenuation exponent α. The compensation amplitude is determined using the least squares method by minimizing the residual disturbance energy ∑(T_original - T_compensated). ² The optimal compensation coefficients are obtained. Phase compensation employs Hilbert transform, shifting the original disturbance signal by 90° and then taking the negative value as the compensation component, ensuring that the compensation signal is completely out of phase with the disturbance signal. Frequency domain compensation uses filters designed according to the propagation modes obtained from path clustering analysis: a high-pass filter with a cutoff frequency of 5Hz is used for fast attenuation path groups to target high-frequency pulse interference; a band-pass filter with a passband of 1-10Hz is used for medium attenuation path groups to target mid-frequency vibration interference; and a low-pass filter with a cutoff frequency of 1Hz is used for slow attenuation path groups to target low-frequency drift interference. The spatial interpolation algorithm interpolates based on the three-dimensional coordinates of the path nodes to generate a continuous spatial compensation field. Time synchronization processing ensures that the compensation waveform is precisely aligned with the original disturbance in time, with the synchronization error controlled within 0.1ms.
[0042] The compensated waveform and inertial navigation data are fused to generate the compensated drill bit attitude. First, the compensated waveform is applied to the original geomagnetic sensor data at the corresponding location according to the spatial coordinates (x, y, z) of the path nodes. Vector superposition eliminates the disturbance at that location, obtaining the purified magnetic field vector. The magnetic field attitude is calculated using a quaternion method, matching the purified magnetic field vector with the Earth's magnetic field reference model to calculate the magnetic field attitude angles (magnetic declination and magnetic inclination). Inertial attitude extraction obtains the noise-filtered and GPS-corrected attitude angle estimates (pitch, roll, and yaw) from the Kalman filter's state output. The attitude fusion algorithm uses a weighted average method, dynamically adjusting the fusion ratio of the magnetic field attitude and inertial attitude based on the compensation effect. The fusion process involves multiplying the purified magnetic field attitude angle by its corresponding weight, multiplying it by the attitude angle obtained from inertial navigation by its corresponding weight, and adding the two to obtain the final fused attitude angle. When the compensation effect is good, the magnetic field attitude dominates (weight 0.7), and the inertial attitude serves as an auxiliary factor (weight 0.3); when the compensation effect is poor, the weight of the magnetic field attitude is reduced to 0.3, and the weight of the inertial attitude is increased to 0.7. The attitude angle output includes pitch angle (angle between the drill bit and the horizontal plane), azimuth angle (direction of the drill bit in the horizontal plane), and tool face angle (rotation angle of the drill bit around its axis), with an angle resolution of 0.1°.
[0043] Step S120: Based on the compensated drill bit attitude, extract the trajectory deviation vector, decompose the trajectory deviation vector into positive deviation components and negative deviation components, and accumulate potential energy of the negative deviation components to generate a deviation energy pool.
[0044] Specifically, the trajectory deviation vector is extracted using the compensated drill bit attitude. Through the three-dimensional spatial vector difference method, the drill bit attitude angle θ_actual(t)=[α(t), β(t), γ(t)] at the current time t is compared component-by-component with the target attitude angle θ_target(t)=[α_0(t), β_0(t), γ_0(t)] of the corresponding designed trajectory, where α is the pitch angle, β is the azimuth angle, and γ is the tool face angle. The deviation vector Δθ(t)=θ_actual(t)-θ_target(t) is then calculated. The magnitude of the deviation vector is calculated as |Δθ|=√[(α-α_0)] ² +(β-β_0) ² +(γ-γ_0) ²The deviation is represented by the overall degree of deviation. Time series analysis constructs historical records of the deviation vector: Δθ(tn), Δθ(t-n+1), ..., Δθ(t), where n is the length of the historical window, set to 50 sampling points. Deviation trend analysis calculates the deviation rate and acceleration through numerical difference to identify the dynamic characteristics of deviation changes. Coordinate system transformation converts the attitude angle deviation of the coordinate system (Z-axis vertically downward, X-axis northward, Y-axis eastward) into the position deviation in the borehole trajectory space. Real-time monitoring of the drilling trajectory and comparison with the target form a precise sequence of trajectory deviation vectors.
[0045] The trajectory deviation vector is decomposed into positive and negative deviation components. Using a geometric projection method, the projection of the deviation vector Δθ onto the tangent direction of the design trajectory is defined as the positive deviation component, and the projection onto the normal direction is defined as the negative deviation component. The local tangent vector of the design trajectory is calculated using the time derivative of the target attitude angle and normalized to a unit tangent vector. The positive deviation component is obtained by multiplying the deviation vector by the unit tangent vector (dot product), representing the advance or lag deviation along the design trajectory direction. The negative deviation component is obtained through vector subtraction, i.e., subtracting the positive deviation component from the deviation vector, representing the lateral deviation from the design trajectory. Component amplitude analysis calculates the magnitude of each component; the magnitude of the positive deviation component reflects the degree of trajectory progress deviation, and the magnitude of the negative deviation component reflects the degree of trajectory direction deviation. Time-domain feature extraction analyzes the mean, variance, and autocorrelation function of each component to identify the periodic and random characteristics of the deviation. Component correlation analysis calculates the cross-correlation function between the positive and negative deviation components, revealing the coupling relationship between the two types of deviations.
[0046] In some embodiments, the step of accumulating potential energy to generate a deviation energy pool from the negative deviation component includes: mapping the negative deviation component to a potential energy gradient field; finding an energy convergence point in the potential energy gradient field; performing energy concentration processing based on the potential energy gradient field and the energy convergence point to form concentrated potential energy; and accumulating the concentrated potential energy to obtain a deviation energy pool.
[0047] The negative deviation component is mapped to a potential energy gradient field. Using the Lagrange potential energy model, the negative deviation component is used as a generalized coordinate to construct the potential energy function U = ½k(Δθ_neg). ² Where k is the system stiffness coefficient, U is the potential energy value, and Δθ_neg is the negative deviation component. The gradient field is calculated using partial derivatives. The gradient vector indicates the direction of the fastest potential energy growth. The field strength distribution is described using the equipotential surface method; different equipotential surfaces correspond to different potential energy levels, forming a hierarchical energy distribution structure. Spatial grid partitioning discretizes the value space of the negative deviation component, forming a three-dimensional grid structure. Each grid point records the corresponding potential energy value and gradient vector. Topological analysis of the field identifies critical points, saddle points, and extreme points in the gradient field. The spatial characteristics of the potential energy distribution are characterized by changes in gradient intensity and direction, revealing the energy concentration pattern of the negative deviation component. The energy characteristics and spatial distribution of the negative deviation component are transformed into a structured potential energy gradient field.
[0048] The energy convergence point is located within the potential energy gradient field. A gradient flow tracing algorithm is employed, starting from various grid points in the gradient field and iteratively calculating along the descent direction of the potential energy gradient. Here, x represents the position coordinates of the negative deviation component in the 3D attitude space, comprising three components: pitch angle negative deviation Δα_neg, roll angle negative deviation Δβ_neg, and yaw angle negative deviation Δγ_neg. The energy convergence point is determined by a gradient magnitude less than a threshold ε and a position change less than δ within five consecutive iteration steps, where ε = 0.001 N·m / rad is the gradient magnitude threshold, and δ = 0.0001° is the position change threshold. Multi-startpoint tracing involves uniformly selecting 1000 starting points from the grid for parallel streamline calculations, and statistically analyzing the distribution of the final convergence position. Convergence region clustering uses the DBSCAN algorithm, grouping convergence points less than 0.005° apart into the same convergence region, identifying primary and secondary convergence regions. Stability analysis of the convergence region is conducted through perturbation testing, applying small-amplitude random perturbations near the convergence point and observing the system's convergence behavior. Basin analysis divides the potential energy gradient field into different potential energy basins, each basin corresponding to a convergence region, and the basin boundary is determined by the watershed line.
[0049] Concentrated potential energy is formed by focusing energy based on the potential energy gradient field and energy convergence points. Using the virtual potential well method, deep potential wells are established at each major convergence point to concentrate the dispersed potential energy from the surrounding area. The potential well model adopts a harmonic oscillator form, V_well(r) = ½k_well·r. ² Where r is the radial distance from the convergence point, and k_well is the potential well stiffness coefficient, set to k_well = 500 N·m / rad. ²Energy transfer paths are established along gradient streamlines, with the reverse direction of the streamlines indicating the energy transfer direction. Energy allocation across multiple convergence points is weighted according to the convergence region classification and the potential energy value of the potential field at the convergence point: the weight coefficient for primary convergence regions is set to 0.7, for secondary convergence regions to 0.3, and among convergence points of the same type, allocation is proportional based on their numerical values in the potential energy function U. Potential energy superposition employs the principle of linear superposition; concentrated potential energy is obtained by weighted superposition of the potential well functions of each convergence point, with the weight coefficients determined according to the convergence region type. The concentration effect is evaluated by calculating the variance of the potential energy distribution before and after concentration; the degree of variance reduction reflects the concentration effect.
[0050] A bias energy pool is obtained by accumulating concentrated potential energy. The distribution of concentrated potential energy is quantified using numerical integration, and the stored energy value is obtained through three-dimensional volume integral calculation. The integration region represents the spatial range where potential energy is significant, and a potential energy threshold is used for truncation; regions below the threshold are ignored. A high-precision Gauss-Legendre numerical integration formula is employed to ensure computational accuracy. Finer meshes are used in regions with high concentrated potential energy density, and the mesh spacing is dynamically adjusted based on the local potential energy gradient. Error control is achieved through adaptive mesh refinement; the mesh is automatically refined when the potential energy difference between adjacent mesh points exceeds a set threshold. The calculated stored energy value is normalized to obtain a standardized stored energy coefficient. A circular buffer structure is used to construct the bias energy pool, with a buffer capacity of 1000 stored energy values, storing historical stored energy data in chronological order. The accumulation strategy employs a sliding window method with a window length of 100 time steps; the accumulated energy is the sum of the stored energy values over the past 100 time steps, enabling dynamic updates to the energy pool. Energy pool status monitoring includes three key indicators: total energy, energy change rate, and energy distribution characteristics. Total energy reflects the currently available corrective energy reserves, while the energy change rate indicates the rate of energy accumulation or consumption. Energy pool capacity management employs a priority-based elimination strategy, prioritizing the elimination of historical data with the lowest energy storage value when the buffer is full. Upper and lower thresholds are set for the energy pool; exceeding these thresholds triggers corresponding energy regulation strategies.
[0051] Step S130: Monitor the charging and discharging of the deviation energy pool to identify stable and unstable boundaries, determine the chaotic control zone between the stable and unstable boundaries, and extract the optimal response interval through the chaotic control zone.
[0052] Specifically, the deviation energy pool is monitored for charging and discharging to identify stable and unstable boundaries. A bidirectional flow metering method is used for continuous monitoring of the deviation energy pool. Charging is identified by detecting a cumulative energy growth rate greater than zero, and discharging is identified by detecting a cumulative energy growth rate less than zero. The charging and discharging flow rates are calculated using a moving average filter: charging flow rate Q_in = Σ(ΔE_positive) / Δt, and discharging flow rate Q_out = Σ(ΔE_negative) / Δt, where ΔE_positive represents positive energy change, ΔE_negative represents negative energy change, and Δt is the time interval. Energy fluctuation amplitude is reflected by the standard deviation σ_E of the energy fluctuation, indicating the stability of the energy pool. Frequency analysis uses Fast Fourier Transform to identify the dominant frequency components. A charge / discharge cycle counter records the number of complete cycles, constructing an energy change time series, and continuously monitoring this data to plot the energy pool's fluctuation curve. Boundaries are identified based on the fluctuation curve. The stable boundary is defined using a statistical control method as μ + 2σ, where μ is the mean fluctuation amplitude and σ is the standard deviation of the fluctuation amplitude. Fluctuations within this boundary are considered normal operating conditions. The instability boundary is determined using extreme value statistics theory and fitted with a generalized extreme value distribution model, set as the threshold corresponding to a high confidence level. Boundary detection employs a sliding window algorithm to identify states near the boundary. A charge-discharge imbalance exceeding a set threshold serves as an instability warning signal. Frequency domain stability is evaluated using power spectral density distribution; predominance of low-frequency components indicates stability, while enhanced high-frequency components indicate a tendency towards instability.
[0053] A chaotic control zone is constructed between the identified stable and unstable boundaries to achieve precise dynamic regulation. Using fractal geometry, the interval between the stable and unstable boundaries is divided into three sub-zones according to the golden ratio of 0.618: the inner control zone, the core control zone, and the outer control zone. The inner control zone is adjacent to the stable boundary, with a range of [B_stable, B_stable + 0.618 × (B_unstable - B_stable)], where B_unstable is the unstable boundary value and B_stable is the stable boundary value. The core control zone is located in the middle and undertakes the main control function. The outer control zone is close to the unstable boundary and requires enhanced control measures. The width of the control zone is determined based on the boundary spacing, with a base width W_base = (B_unstable - B_stable) / 3, ensuring that the widths of the three sub-zones are equal. A control strategy mapping establishes the correspondence between different sub-zones and control parameters: linear control is used in the inner control zone, nonlinear control in the core control zone, and enhanced control in the outer control zone. The dynamic boundary of the control zone is adjusted by monitoring changes in the boundary position in real time, and the control zone division is automatically updated when the boundary shifts significantly. A boundary buffer is set at the intersection of each sub-zone, with a width of 5% of the sub-zone width, to avoid abrupt changes during control switching. Control priority increases from the inner control area to the outer control area, with the outer control area having the highest control priority.
[0054] In some embodiments, the step of extracting the optimal response interval through the chaotic control band includes: identifying energy-rich regions from the chaotic control band; obtaining response sensitivity coefficients based on the energy-rich regions; performing interval optimization processing based on the response sensitivity coefficients to form candidate response intervals; and selecting the interval with the highest stability from the candidate response intervals as the optimal response interval.
[0055] Enriched regions with high energy distribution density were identified from the constructed chaotic control zone. A density clustering algorithm was used to divide the energy distribution within the control zone into a two-dimensional grid based on spatial location and time, with the grid resolution set to 1% of the control zone width. Density calculation was performed using kernel density estimation. The enrichment region was defined as a density value exceeding 1.5 times the average density; regions with a density greater than 1.5 times the average density were identified as energy-rich regions. Spatial connectivity analysis employed the 8-neighborhood connectivity rule, merging adjacent high-density grids into continuous enriched regions. Enriched regions were classified into three levels based on their maximum density value: Level 1 enriched regions were those with a maximum density exceeding 2 times the average density; Level 2 enriched regions were those with a maximum density between 1.5 and 2 times the average density; and Level 3 enriched regions were those with a maximum density between 1.2 and 1.5 times the average density. The stability of enriched regions was quantified by density changes within a time window; small changes indicated stable enriched regions, while large changes indicated unstable enriched regions. The enrichment intensity calculation reflected the total energy content of the enriched region. Boundary clarity analysis evaluates the boundary by calculating the gradient of the enriched region's edge; a larger gradient indicates a clearer boundary.
[0056] For example, obtaining the response sensitivity coefficient based on the energy-rich region includes: decomposing the energy-rich region into a core-rich region and a peripheral-rich region; setting a reference sensitivity based on the core-rich region; dynamically adjusting the reference sensitivity using the peripheral-rich region to form an adjusted sensitivity; and generating a response sensitivity coefficient based on the adjusted sensitivity.
[0057] The energy-rich region is decomposed into a core-rich region and a peripheral-rich region with different response characteristics. Isodensity lines are used for layering, with the maximum density value of the enriched region as the benchmark. Regions with a density greater than 80% of the maximum density value are defined as core-rich regions, and regions with a density between 50% and 80% of the maximum density value are defined as peripheral-rich regions. Geometric feature analysis calculates the area, perimeter, and shape factor of each region. Core-rich regions typically exhibit compact ellipses or circles, while peripheral-rich regions exhibit ring-shaped or irregular shapes. Density distribution fitting uses a two-dimensional Gaussian model to describe the spatial distribution characteristics of the enriched regions. The core region boundary is determined by density gradient calculation, with the isodensity line corresponding to the maximum gradient value serving as the core region boundary. The effective thickness of the peripheral region is determined by radial density decay analysis, with the decay length constant defined as the distance at which the density drops to one-third of the natural logarithm of the peak value. Connectivity checks ensure that the decomposed regions maintain spatial continuity; disconnected regions are repaired through morphological connection operations. Spatial weight allocation is determined based on the area ratio of the core and peripheral regions.
[0058] A baseline sensitivity is set based on the core enriched region. The baseline sensitivity is set using a density-sensitivity linear relationship model, establishing a quantitative mapping relationship between the density characteristics of the core region and the response sensitivity: S_baseline = k × ρ_max + b, where S_baseline is the baseline sensitivity, k is a scaling factor set to 0.5, b is a bias parameter set to 0.1, and ρ_max is the peak density of the core enriched region. Geometric correction considers the ratio of the major and minor axes of the core region. When the ratio is greater than 2, it indicates that the core region is elongated, and the sensitivity needs directional correction. The correction coefficient is adjusted according to the ellipticity. The spatial position effect is described by a distance attenuation function; the sensitivity decreases with distance from the core center, and the attenuation follows an exponential decay characteristic. The frequency characteristics of the baseline sensitivity adopt a single-pole model, with the corner frequency positively correlated with the core region density. High-density regions have a wider frequency response bandwidth.
[0059] The adjustment sensitivity is formed by dynamically adjusting the reference sensitivity using the peripheral enrichment region. A feedback control principle is employed, monitoring the density change in the peripheral enrichment region as the feedback signal; the adjustment amplitude is proportional to the density change. Radial adjustment weights are allocated based on the density gradient at different radial locations in the peripheral region, with the peripheral region closer to the core region receiving a larger adjustment weight. The time constant describes the dynamic response characteristics of the adjustment process. Adjustment saturation limiting prevents over-adjustment, with the adjustment amplitude limited to ±20% of the reference sensitivity. Frequency-dependent adjustment considers the differences in adjustment effects across different frequency components, using different adjustment gains for high and low frequencies. Adaptive threshold setting is dynamically adjusted based on the statistical characteristics of the peripheral region's density; adjustment is only triggered when the density change exceeds the threshold. Adjustment stability is monitored by evaluating the sensitivity coefficient of variation before and after adjustment; a decrease in the coefficient of variation indicates effective adjustment.
[0060] Based on the established adjustment sensitivity, the final response sensitivity coefficient is generated through nonlinear correction and normalization. Single-input processing is employed, directly conditioning and optimizing the signal based on the adjustment sensitivity. Nonlinear correction addresses sensitivity saturation under large signal conditions, using a hyperbolic tangent model to prevent over-response. Frequency characteristic optimization considers the frequency response characteristics of the adjustment sensitivity, using digital filters for frequency compensation to ensure consistent response across the entire frequency band. Temperature drift compensation uses a temperature coefficient correction, based on the difference between the current and reference temperatures. Noise suppression is achieved through low-pass filtering, with the filter cutoff frequency set to twice the signal bandwidth, effectively suppressing high-frequency noise. Linearity correction is achieved through multi-point calibration, setting multiple calibration points within the operating range of the adjustment sensitivity and obtaining correction parameters through linear fitting. Dynamic range compression compresses the wide dynamic range adjustment sensitivity to the standard output range, with an adjustable compression ratio ranging from 1:1 to 10:1. Normalized output normalizes the processed sensitivity to the 0-1 range.
[0061] Based on the obtained response sensitivity coefficients, interval optimization is performed to form candidate response intervals. A multi-objective optimization method is adopted, with optimization objectives including maximizing the peak sensitivity coefficient, maximizing the sensitivity coefficient gradient, and optimizing the stability of the sensitivity coefficient. Multi-objective optimization is solved using the non-dominated sorting genetic algorithm NSGA-II, with a population size of 100 and 200 generations. Constraints include: a lower limit for the sensitivity coefficient greater than 0.1, monotonicity requirement for the sensitivity coefficient within the interval, and a minimum interval length limit. The boundaries of candidate intervals are determined through cluster analysis. The response sensitivity coefficient values are clustered according to their numerical values using K-means, with 5 clusters, and each cluster boundary corresponds to a candidate interval. Interval length optimization is based on the local extreme value distribution of the response sensitivity coefficients, defining interval boundaries near extreme points to ensure that each interval contains at least one peak sensitivity coefficient. Overlapping intervals are automatically merged into a single interval when the difference in the maximum response sensitivity coefficient of adjacent candidate intervals is less than 10%. The interval scoring mechanism takes into account the weighted combination of three objective functions, with weight coefficients of 0.5, 0.3, and 0.2, respectively.
[0062] The optimal response interval with the highest stability was selected from the generated set of candidate response intervals. A multi-indicator system was adopted, including short-term stability, long-term stability, and anti-interference capability. Short-term stability was calculated through analysis of variance; the smaller the variance, the better the short-term stability. Long-term stability was evaluated through trend analysis; intervals with low drift rates exhibited good long-term stability. Anti-interference capability was quantified by signal-to-noise ratio (SNR); intervals with high SNR demonstrated strong anti-interference capability. The comprehensive stability index weighted the three aspects. Operability assessment considered the control difficulty and response delay of the interval location; intervals located in the core control area had better operability. Interval ranking was performed using a weighted combination of stability index and operability, and a final score was calculated. The optimal interval was selected using a threshold discrimination method, choosing the interval with the highest score and a stability index exceeding a preset threshold as the optimal response interval. Based on multi-dimensional stability analysis and comprehensive evaluation, the optimal response interval was finally determined, effectively preventing oscillations in the correction control.
[0063] Step S140: Determine the active correction vector based on the optimal response range and positive deviation component, control the release of the deviation energy pool to obtain the energy release amount, and vector superimpose the active correction vector and the energy release amount to generate a composite correction command.
[0064] Specifically, the active correction vector is determined based on the optimal response range and the positive deviation component. A proportional-integral-derivative (PID) control algorithm is employed, using the positive deviation component as the control error input. The correction vector is defined as F_correction = K_p × Δθ_positive + K_i × ∫Δθ_positivedt + K_d × d(Δθ_positive) / dt, where K_p is the proportional gain coefficient, K_i is the integral gain coefficient, K_d is the derivative gain coefficient, and Δθ_positive represents the positive deviation component. The gain coefficient is adaptively adjusted based on the position of the optimal response range: when the response range is within the core control region, K_p = 2.5, K_i = 0.8, and K_d = 0.3; when within the inner control region, the gain is reduced by 20%; and when within the outer control region, the gain is increased by 30%. The vector direction is determined using gradient descent, with the correction direction opposite to the gradient direction of the positive deviation component, ensuring that the correction force points towards the target trajectory. The vector amplitude is limited using a saturation function, restricting the correction vector amplitude to within the maximum correction force F_max (500N) to prevent over-correction from causing system instability. The response interval weighting factor is dynamically adjusted based on the interval's stability index; higher stability results in a larger weight. Vector component decomposition decomposes the three-dimensional correction vector into radial component F_r, tangential component F_t, and axial component F_a, corresponding to the radial, tangential, and axial hydraulic cylinder control channels, respectively. Time window smoothing employs a 5-point moving average filter to eliminate the impact of high-frequency noise on correction accuracy.
[0065] The energy storage release is obtained by controlling the release of the deviation energy pool. Based on the energy-control mapping principle, the hydraulic system control parameters are modulated by the energy storage state. The release amount, Q_release, represents the equivalent control energy extracted from the energy pool, calculated as Q_release = η × E_equivalent, where η is the energy utilization efficiency coefficient set to 0.75, and E_equivalent is the equivalent control energy corresponding to the correction vector. The formula for calculating the equivalent control energy is E_equivalent = |F_correction|. ²×Δt / (2×k_equiv), where |F_correction| is the magnitude of the correction vector, Δt is the control time interval, and k_equiv is the system equivalent stiffness coefficient determined according to the drill pipe parameters. Release priority is based on the timestamp and importance ranking of the energy storage values in the energy pool, prioritizing the release of historical data with longer storage times and smaller energy storage values, while ensuring that the cumulative released energy storage value reaches the required equivalent control energy. The release rate control adopts an exponential decay model: dQ / dt=(Q_target-Q_current) / τ, where Q_target is the target release amount, Q_current is the current release amount, and τ is a time constant set to 2 seconds. Energy pool status monitoring tracks changes in total energy E_total, available energy E_available, and reserved energy E_reserved in real time. When available energy falls below 20% of the total energy, an energy protection mode is triggered, limiting the extraction of equivalent control energy. Release path optimization uses an optimal allocation algorithm to determine the proportion of energy extracted from different energy storage units, minimizing control performance loss. Release precision control is achieved through closed-loop feedback, with the deviation between the actual extracted equivalent energy and the target demand being less than 5%. The batch release mechanism employs parallel extraction when a large amount of control energy is required, simultaneously extracting equivalent energy from multiple energy storage units according to weighted proportions to ensure timely control response.
[0066] In some embodiments, the step of vector superimposing the active correction vector with the energy storage release amount to generate a composite correction command includes: constructing an energy release timing sequence based on the energy storage release amount; embedding the active correction vector into the energy release timing sequence to form a time-varying vector sequence; performing envelope modulation on the time-varying vector sequence to output a modulation vector; and converting the modulation vector into a composite correction command.
[0067] Based on the acquired energy release volume, an energy release timing sequence is constructed to coordinate corrective actions. Based on the piecewise linear principle, the release process is divided into three stages: a start-up stage, a stabilization stage, and a termination stage. The start-up stage lasts 0.5 seconds, with the release rate linearly increasing from 0 to the target rate; the stabilization stage maintains a constant release rate, and its duration is determined by the equivalent control energy Q_release; the termination stage lasts 0.3 seconds, with the release rate linearly decaying to 0. Timing precision control uses a high-resolution clock with a time resolution of 1ms to ensure accurate release timing. A synchronization triggering mechanism is implemented through hardware interrupts, automatically triggering the next stage when the energy release reaches a preset threshold. Timing compensation considers system latency and response time, with an advance factor set to 1.2 times the system response time. Multi-channel timing supports simultaneous control of multiple release channels, with each channel timing independently but maintaining overall synchronization. A timing fault tolerance mechanism automatically switches to a backup timing sequence when an abnormality occurs in a certain stage. Periodic analysis of the release timing identifies recurring patterns for predicting and optimizing future release strategies.
[0068] The determined active correction vector is embedded into the constructed energy release time sequence to form a time-varying vector sequence. Using time interpolation, at each time point t_i (the i-th time point in the release time sequence), the correction vector value is equal to the product of the active correction vector and the time weighting function. The weighting function is designed with a Gaussian distribution, reaching its maximum value at the center of the release time sequence and gradually decaying towards both ends, ensuring synchronization between the correction force and the energy release process. Vector amplitude modulation is dynamically adjusted according to the real-time state of energy release; the vector amplitude increases when the release rate is high and decreases when the release rate is low. Component synchronization ensures that all vector components remain consistent in time, avoiding deviations in the correction effect caused by component asynchrony. Time-varying characteristic analysis identifies the rate and trend of vector change by calculating the time derivative of the vector sequence. Sequence smoothing uses a Butterworth low-pass filter with a cutoff frequency set to 10Hz to eliminate high-frequency jitter. Periodicity detection of the vector sequence is performed through autocorrelation function analysis to identify the existence of periodic patterns.
[0069] The constructed time-varying vector sequence is envelope-modulated to output a modulated vector. Based on digital signal processing methods, the amplitude envelope features of the vector sequence are extracted using filtering and envelope detection techniques. Envelope signal extraction is achieved by calculating the amplitude of the time-varying vector sequence; the envelope signal reflects the intensity distribution of the vector sequence. Frequency analysis identifies the main frequency components of the vector sequence using Fast Fourier Transform (FFT) and optimizes modulation parameters to enhance the useful signal. Modulation depth control employs an automatic gain control algorithm to maintain a stable dynamic range of the output signal and avoid over-saturation or under-modulation. Frequency compensation considers the system's frequency response characteristics and performs differentiated processing on different frequency components. Phase correction ensures that the modulated vector remains synchronized with the original timing sequence, with phase deviation controlled within ±5 degrees. The modulated vector output is standardized, with the amplitude normalized to the standard output range for easy subsequent digitization. Signal quality assessment is performed through signal-to-noise ratio (SNR) calculation to ensure that the modulated signal quality meets control requirements.
[0070] The obtained modulation vector is converted into a composite correction command. Digital signal processing methods are used to discretize the continuous modulation vector into digital control codes. The encoding format uses binary two's complement representation, with positive values corresponding to the extension direction and negative values corresponding to the retraction direction. Multi-channel command generation generates independent control commands for the radial, tangential, and axial hydraulic cylinders, with each channel corresponding to one component of the modulation vector. Command timing synchronization uses a unified clock reference to ensure time consistency of multi-channel commands, with synchronization errors controlled within 0.1ms. Safety limit checks perform boundary checks on the generated commands, automatically limiting commands exceeding the safety range to the allowable interval. Command smoothing uses a first-order inertial element to avoid hydraulic shocks caused by sudden command changes; the smoothing time constant is set to 50ms. Command priority management prioritizes commands in case of conflict, executing them according to their urgency. The communication protocol uses the CAN bus standard, supporting real-time command transmission and status feedback.
[0071] Step S150: The hydraulic cylinder is controlled by the composite correction command to generate thrust distribution. Based on the thrust distribution, a deflection torque is generated on the drill bit. The deflection torque is used to change the drilling direction to obtain the actual correction amount. The oscillation characteristics of the actual correction amount are extracted to form the system oscillation parameters.
[0072] Specifically, the hydraulic cylinders generate thrust distribution through composite correction commands. A multi-channel parallel drive method is employed, controlling the extension and retraction of the radial, tangential, and axial hydraulic cylinders separately according to the control parameters in the composite correction commands. The hydraulic system response is achieved through electro-hydraulic proportional valves, with the valve opening linearly corresponding to the command amplitude to ensure precise control. Thrust calculation is based on Pascal's law F=P×A, where F is the hydraulic cylinder thrust, P is the working pressure, and A is the effective piston area. Radial, tangential, and axial hydraulic cylinders have different piston areas to meet control requirements in different directions. Multi-cylinder thrust coordination is achieved through synchronous control, keeping the thrust output error of each hydraulic cylinder within acceptable limits to avoid side effects caused by thrust imbalance. Dynamic response optimization, through feedforward compensation technology, significantly shortens the system response time and improves control real-time performance. The thrust monitoring system uses high-precision force sensors to measure the actual output thrust of each hydraulic cylinder in real time, ensuring control accuracy. The thrust vectors in the three directions are combined to form a spatial thrust distribution containing three components: radial thrust F_r, tangential thrust F_t, and axial thrust F_a. This thrust distribution directly acts on the drill bit to achieve deviation control.
[0073] The generated thrust distribution acts on the drill bit to produce a deflection torque. Employing the principle of combined thrust, the thrust distribution generated in the initial stage is applied to the geometric center of the drill bit, producing a deflection torque M = r_center × F_distribution, where r_center is the position vector of the drill bit's geometric center, and F_distribution is the thrust distribution vector. The drill bit's geometric parameters include a diameter of 216 mm, a length of 300 mm, and a radial distance of 180 mm from the hydraulic cylinder to the drill bit's center. The radial component F_r of the thrust distribution primarily generates the pitching torque, the tangential component F_t generates the yaw torque, and the axial component F_a generates the rolling torque. These three components work together to achieve deflection control of the drill bit in three rotational degrees of freedom. The torque amplification effect is achieved through the lever principle; the radial arrangement of the hydraulic cylinder forms an effective torque arm, and the amplification factor is determined based on the ratio of the arm length to the drill bit radius. The drill string stiffness influence correction considers the impact of drill pipe elastic deformation on torque transmission. The stiffness coefficient K_drill is calculated based on the drill pipe material's elastic modulus and geometry, with a typical range of 1000-5000 N·m / rad. Formation resistance coupling analysis examines the reaction torque generated by the drill bit's contact with the formation. The formation resistance coefficient is determined based on lithology: 0.3 for soft rock, 0.5 for medium-hard rock, and 0.8 for hard rock. Effective thrust distribution transmission is assessed through drill string connection stiffness, with a transmission efficiency of 85%-95%, depending on the tightness of the drill string joints. Transient torque response analysis examines the dynamic process from thrust distribution to the establishment of deflection torque; the response time is influenced by hydraulic system characteristics and the spatial configuration of the thrust distribution.
[0074] The actual correction amount is obtained by changing the drilling direction using the generated deflection torque. The angular momentum theorem is used for analysis, where angular acceleration α = M / J, and M is the deflection torque and J is the drill bit's moment of inertia. The change in direction angle is obtained through integration, with the integration time window set to the duration of the correction action. The drilling trajectory correction amount is calculated geometrically based on the angle change caused by the deflection torque; the radial and axial offsets together constitute the spatial correction vector. The actual correction amount is verified using a borehole trajectory measurement system, including gyroscope inclination measurement, magnetometer azimuth measurement, and accelerometer gravity tool face measurement, with a measurement accuracy of ±0.1°. The correction effect is evaluated by comparing the trajectory parameters before and after correction; the correction efficiency is defined as the ratio of the actual correction amount to the target correction amount, ideally approaching 100%. Multi-factor coupling analysis considers the influence of drilling parameters on the correction effect, including rotational speed, drilling pressure, and displacement, establishing a multivariate regression model to predict the correction effect. Time delay compensation is applied to address signal transmission delays caused by the measurement system and drill string length. The correction measure verification ensures data reliability through repeated measurements and statistical analysis, with measurement repeatability better than 0.05°.
[0075] The system oscillation parameters are formed by extracting oscillation characteristics from the actual correction amount. A multi-scale decomposition of the correction amount's time series is performed using a joint time-frequency domain analysis method. Time-domain feature analysis calculates the statistical characteristic parameters of the correction amount: the mean reflects the system's static bias, the standard deviation characterizes the oscillation amplitude, the skewness describes the symmetry of the oscillation, and the kurtosis reflects the sharpness of the oscillation. Frequency-domain features are obtained through power spectral density analysis, identifying the dominant oscillation frequency and spectral bandwidth. The oscillation period is identified using the autocorrelation function; the first non-zero maximum corresponds to the dominant oscillation period. Damping characteristics are evaluated using the logarithmic decline rate; the damping ratio reflects the attenuation characteristics of the system oscillation. Nonlinear feature analysis uses the Lyapunov exponent to evaluate the system's chaotic characteristics. Wavelet transform is used to extract the time-frequency localization features of the oscillation, reflecting the oscillation intensity distribution at different time scales. Oscillation mode identification uses a mode decomposition algorithm to decompose complex oscillations into multiple single-mode components, each corresponding to a specific physical mechanism. Comprehensive analysis forms a complete set of system oscillation parameters, including the dominant oscillation frequency, damping ratio, oscillation amplitude, mode order, and nonlinear coefficients, comprehensively characterizing the system's dynamic characteristics.
[0076] Step S160: Based on the system oscillation parameters and the chaotic control zone, position mapping is performed to generate a stability index. The stability index is superimposed with the thrust distribution to form a fine correction force. The chaotic edge control field is constructed based on the fine correction force and the negative deviation component.
[0077] In some embodiments, the step of generating a stability index based on the position mapping between the system oscillation parameters and the chaotic control zone includes: decomposing the system oscillation parameters into periodic oscillation components and random oscillation components; performing trajectory tracking on the periodic oscillation components within the chaotic control zone to form a periodic trajectory; performing probability distribution analysis on the random oscillation components to determine the diffusion boundary; and generating a stability index based on the relative positional relationship between the periodic trajectory and the diffusion boundary.
[0078] The system oscillation parameters are decomposed into periodic oscillation components and random oscillation components. Periodic oscillation component identification uses spectral analysis to determine the main periodic components and extract stable oscillation modes. Random oscillation component extraction employs residual signal analysis; the remaining portion after removing the periodic components constitutes the random component. Component energy distribution typically involves periodic components accounting for 60%-80% of the total energy, and random components accounting for 20%-40%. Characteristics of periodic components include period length, amplitude range, and phase characteristics. Characteristics of random components are obtained through statistical analysis, including statistical features such as mean, variance, and distribution shape. Decomposition integrity is ensured using signal reconstruction methods to guarantee that the sum of the energies of each component after decomposition equals the energy of the original signal. The time window selection is determined based on the main oscillation period, with the window length set to 5-10 times the main period.
[0079] The decomposed periodic oscillation components are tracked within the chaotic control zone to form periodic orbits. Orbit tracking maps the amplitude characteristics of the periodic oscillation components to a one-dimensional interval of the chaotic control zone, determining their positional distribution within the inner, core, and outer control regions based on the oscillation amplitude magnitude. A linear interpolation method is used to normalize the numerical range of the oscillation amplitude and map it to the positional coordinates within the control zone. Orbit construction connects the mapping points at consecutive time points to form the motion trajectory within the control zone. Periodic orbit identification is based on the closure of the trajectory, with the starting and ending points coinciding within a complete cycle. Orbit shape classification is based on geometric characteristics, categorized into elliptical, spiral, and composite types, each reflecting different oscillation characteristics. The orbit boundary is determined by calculating the envelope of the trajectory points, forming the spatial boundary of the periodic orbit. Orbit stability is evaluated through comparison of orbits over multiple consecutive cycles; stable orbits exhibit good reproducibility. Orbit parameter extraction includes orbit area, center position, and main geometric dimensions.
[0080] The diffusion boundary is determined by probability distribution analysis of the decomposed random oscillating components. Statistical methods are used to estimate the distribution characteristics of the random oscillating components within the chaotic control zone. Distribution fitting is performed using kernel density estimation to establish the probability density function of the random components. Distribution parameters include the center location, diffusion range, and distribution shape, reflecting the statistical characteristics of the random oscillations. The diffusion boundary is defined as the boundary of the region containing 95% probability mass, determined by probability contour lines. The boundary shape is typically an ellipse or an irregular closed curve, depending on the correlation structure of the random components. Multidimensional diffusion considers the diffusion differences of the random components in different directions, forming an anisotropic diffusion boundary. Boundary calculation uses the isoprobability density line method, connecting points with the same probability density value to form the boundary curve. Diffusion intensity is characterized by the size of the boundary area; a larger area indicates stronger randomness. Boundary dynamics are monitored over time using sliding window analysis.
[0081] Stability indices are generated based on the relative positional relationship between the periodic orbit and the diffusion boundary. Positional relationship analysis calculates the shortest distance (d_stable) from the periodic orbit to the diffusion boundary and the shortest distance (d_unstable) to the instability boundary of the chaotic control zone through geometric distance calculation. The stability index is defined as S_stability = d_stable / (d_stable + d_unstable), with a value ranging from 0 to 1; values closer to 1 indicate greater system stability. Relative position assessment considers the degree of offset between the orbit center and the diffusion center; the offset reflects the system's bias. Overlap analysis calculates the ratio of the intersection area between the periodic orbit and the diffusion boundary to the orbital area; the overlap ratio affects stability assessment. Spatial distribution coordination is assessed through the geometric matching degree between the orbit and the boundary; good coordination indicates harmonious system operation. A comprehensive stability evaluation is formed by considering multiple parameters, including periodicity intensity, randomness level, and spatial distribution characteristics. Index dynamics are analyzed through trend changes over continuous time windows; the trend direction predicts the development of system stability. The stability threshold is set according to system safety requirements; control intervention is triggered when the value falls below the threshold.
[0082] A fine-grained correction force is formed by superimposing stability indices and thrust distribution. A tiered control strategy is adopted based on the numerical range of the stability indices: small-amplitude fine-tuning in high-stability regions, medium-amplitude adjustment in medium-stability regions, and large-amplitude enhanced adjustment in low-stability regions. The gradient adjustment amount is calculated using a proportional control algorithm: ΔP_gradient = K_stability × (S_target - S_stability), where K_stability is the stability gain coefficient, S_target is the target stability, and S_stability is the current stability index. The gradient adjustment amount is allocated to the radial, tangential, and axial hydraulic cylinder channels, with the allocation weights determined according to the current correction requirements. A filtering algorithm is used for pressure gradient smoothing to avoid pressure abrupt changes impacting the system. A vector synthesis method is used to convert the pressure increment generated by gradient adjustment into a thrust increment, which is then multiplied by the piston area to obtain the thrust increment. The fine-grained correction force is calculated as F_fine = F_distribution + ΔF, where F_distribution is the original thrust distribution, ΔF is the adjusted thrust increment, and vector superposition maintains the independence of each component. The superposition weights are dynamically adjusted according to the correction requirements. During emergency corrections, the adjustment weight is increased, while during smooth corrections, the original thrust weight is maintained. Force vector normalization ensures that the fine correction force remains within the system's load-bearing capacity.
[0083] In some embodiments, constructing a chaotic edge control field based on the fine correction force and the negative deviation component includes: obtaining the center of action of the fine correction force; generating a field strength gradient distribution based on the center of action; guiding the negative deviation component to redistribute energy along the field strength gradient distribution to form a redistribution field; and performing edge sharpening processing on the redistribution field to form a chaotic edge control field.
[0084] The process involves several steps: First, the center of action of the fine-grained corrective force is determined. Considering the magnitude and location of each component of the force, the equivalent point of action of the resultant force is determined using a torque balance method. Geometric center correction considers the actual geometric layout of the drill bit and hydraulic cylinder and is adjusted based on the system's symmetry. Dynamic center tracking monitors the trajectory of the center of action over time, and center stability is assessed using positional variance. Multi-scale center analysis calculates the center of action at different time scales; short-term scales reflect instantaneous characteristics, while long-term scales reflect trend characteristics. Center offset is calculated relative to the theoretical geometric center, measuring the distance and direction of the offset; the magnitude of the offset reflects the degree of system imbalance. A center stability domain defines the allowable range of variation for the center of action; the radius of the stability domain is set to 20% of the drill bit radius, triggering a corrective action when the range is exceeded. The center position record saves the historical trajectory of the center of action at a frequency of 1 Hz for system behavior analysis and fault prediction. The spatial distribution characteristics of the center of action are obtained through statistical analysis, including the mean, standard deviation, and distribution range of the center position.
[0085] A field strength gradient distribution is generated based on a defined center of action. A control strength distribution C(r) = k × |F_fine| × exp(-r / λ) is established with the center of action as the origin, where k is the control strength coefficient set to 0.1, |F_fine| is the modulus of the fine correction force, λ is the characteristic attenuation length set to the drill bit radius, and r is the distance from the center of action. The gradient distribution calculates the spatial derivative of the control strength, with the gradient direction pointing towards the direction of the fastest increase in control strength, used to determine the priority direction of correction control. The field strength function is selected according to control requirements; linear functions are suitable for uniform control, while exponential functions are suitable for localized intensification control. Iso-field strength lines are distributed by numerical calculation of the field strength, forming concentric circles or ellipses of contour lines, with the spacing between contour lines reflecting the gradient strength. The gradient strength distribution is calculated numerically, and the spatial distribution of the gradient magnitude reflects the drastic change in field strength. Field strength boundaries are determined by setting a threshold; regions with field strength below the threshold are considered the boundaries of the control field. A multi-level gradient structure divides the control field into a core region, a transition region, and an edge region, with different control strategies employed in different regions. Gradient continuity is checked through numerical differentiation to ensure the smoothness of the field strength distribution and avoid control instability caused by abrupt gradient changes.
[0086] A redistribution field is formed by guiding the negative deviation component along the constructed field strength gradient distribution to redistribute energy. The negative deviation component Δθ_negative is used as a control reference quantity, and its weights are redistributed under the guidance of the control strength gradient. A weighted allocation algorithm is used in the redistribution process, determining the influence weight of the negative deviation component based on the local control strength; the weight is proportional to the control strength. Control weight flow analysis calculates the trend of weight allocation changes, with the weight adjustment direction shifting from low control strength regions to high control strength regions along the gradient direction. Redistribution efficiency is evaluated using the principle of weight sum conservation, ensuring that the total weight remains unchanged before and after redistribution. Energy accumulation zone identification uses density analysis to identify high-density regions after energy redistribution, with the accumulation threshold set to 1.5 times the average density. Steady-state analysis of the redistribution field involves long-term simulation to observe the characteristics of energy distribution reaching equilibrium, with an equilibrium time constant of 5-10 seconds. Redistribution rate control balances redistribution speed and stability by adjusting the intensity of the field strength gradient. Field boundary effects consider the influence of boundary conditions on energy redistribution, employing reflective boundaries to avoid energy loss. The diffusion mechanism analysis reveals the diffusion characteristics of energy during the redistribution process, showing that the degree of diffusion is related to the local gradient strength.
[0087] Edge sharpening is applied to the redistribution field to form a chaotic edge control field. Chaotic control leverages the extreme sensitivity of chaotic phenomena to initial conditions; when the borehole trajectory approaches a critical instability state, applying a small control action can produce a significant correction effect, similar to how a slight push at the equilibrium critical point can change the direction of motion. Edge sharpening employs numerical differentiation methods to calculate the second derivative of the control intensity of the redistribution field, identifying regions with drastic gradient changes as control boundaries. The sharpening algorithm adjusts the sharpening coefficient to enhance boundary contrast while maintaining the continuity of the field intensity distribution. Multi-scale processing performs edge detection and enhancement at different spatial scales, ensuring that both local details are captured and the overall structure is preserved. The sharpening process focuses on strengthening the transition region of control intensity, making the originally blurred control boundaries clear and sharp, forming a distinct strong control region, weak control region, and transition zone. Nonlinear transformations enhance chaotic characteristics, keeping the control field smooth in stable regions and exhibiting abrupt changes in boundary regions. The resulting chaotic edge control field has a clear hierarchical structure and sensitive boundary response characteristics. It can promptly identify and trigger corresponding control strategies when the borehole trajectory approaches the critical point, and prevent the borehole from deviating from the target path through precise control of the edge region.
[0088] Step S170: Obtain the conversion efficiency between the deviation energy pool and the actual correction amount, and generate a real-time correction control amount based on the chaotic edge control field and the conversion efficiency.
[0089] Specifically, the conversion efficiency between the deviation energy pool and the actual correction amount is obtained. An energy audit method is used to statistically analyze the total energy, used energy, and remaining energy storage in the deviation energy pool to assess energy utilization. By tracking the charging and discharging history of the energy pool in real time, a correlation between energy consumption and correction effect is established. Energy consumption is mapped to the actual correction amount, and the conversion coefficient K = Δθ / E is calculated, where Δθ is the actual correction angle and E is the corresponding energy consumed; this coefficient reflects the correction effect produced per unit of energy. The conversion coefficient is normalized, using the historical best case as a reference benchmark. The mapping function is selected by comparing multiple models to find the best fit, and the parameters are calibrated using the least squares method. Considering the differences in conversion characteristics under different energy storage levels, a piecewise approach is used to describe the nonlinear relationship. The reliability range of the coefficients is determined through statistical analysis. Based on the correlation analysis between energy consumption patterns and correction effects, dynamically updated conversion efficiency parameters are obtained.
[0090] Real-time correction control quantities are generated based on chaotic edge control fields and conversion efficiency. An adaptive modulation method is employed, using conversion efficiency as a modulation factor and the chaotic edge control field as a control template; these two are combined to generate real-time control commands. The spatial distribution of the control quantity is determined by the control field, and the control intensity is adjusted by conversion efficiency to ensure efficient energy utilization. Spatial interpolation transforms discrete control field data into a continuous distribution, ensuring control smoothness. A time synchronization mechanism ensures that conversion efficiency updates and control field evolution are coordinated. The control quantity is limited within a safe range to prevent over-control from causing system instability. The total control quantity is decomposed into three channels: radial, tangential, and axial, and dynamically allocated according to correction requirements. The control quantity undergoes smoothing processing to eliminate high-frequency disturbances. The real-time monitoring system acquires drill bit status information through multi-sensor fusion and employs filtering algorithms to improve measurement accuracy. The generated control quantity drives the hydraulic actuator, achieving precise control through an electro-hydraulic conversion system. The correction effect is evaluated in real-time by comparing trajectory deviations. Through the synergistic effect of conversion efficiency and the control field, real-time monitoring and intelligent correction control of the borehole trajectory are achieved.
[0091] To implement the above-described method embodiments, a real-time borehole trajectory monitoring and intelligent deviation control method is proposed to achieve the corresponding functions and technical effects. See also... Figure 2 , Figure 2 This diagram illustrates a structural block diagram of a real-time borehole trajectory monitoring and intelligent deviation correction control device 200 according to an embodiment of this application. For ease of explanation, only the parts relevant to this embodiment are shown. The real-time borehole trajectory monitoring and intelligent deviation correction control device 200 provided in this embodiment includes:
[0092] The data fusion module 201 is used to collect inertial navigation data and geomagnetic sensor data during the drilling process, perform distortion analysis on the geomagnetic sensor data to obtain magnetic field disturbance characteristics, and use the magnetic field disturbance characteristics to perform waveform compensation fusion to generate the compensated drill bit attitude.
[0093] Energy conversion module 202 is used to extract trajectory deviation vector based on the compensated drill bit attitude, decompose the trajectory deviation vector into positive deviation component and negative deviation component, and accumulate potential energy of the negative deviation component to generate deviation energy pool.
[0094] Chaos analysis module 203 is used to monitor the charging and discharging of the deviation energy pool to identify stable and unstable boundaries, determine a chaotic control band between the stable and unstable boundaries, and extract the optimal response interval through the chaotic control band.
[0095] The instruction generation module 204 is used to determine an active correction vector based on the optimal response range and the positive deviation component, control the release of the deviation energy pool to obtain the energy release amount, and vector superimpose the active correction vector and the energy release amount to generate a composite correction instruction.
[0096] The execution control module 205 is used to control the hydraulic cylinder to generate thrust distribution through the composite correction command, generate deflection torque based on the thrust distribution acting on the drill bit, change the drilling direction using the deflection torque to obtain the actual correction amount, and extract the oscillation features of the actual correction amount to form system oscillation parameters.
[0097] The stabilization adjustment module 206 is used to generate a stability index by mapping the system oscillation parameters with the chaotic control band, form a fine correction force by superimposing the stability index with the thrust distribution, and construct a chaotic edge control field based on the fine correction force and the negative deviation component.
[0098] The optimization control module 207 is used to obtain the conversion efficiency between the deviation energy pool and the actual correction amount, and to generate a real-time correction control amount based on the chaotic edge control field and the conversion efficiency.
[0099] The aforementioned real-time borehole trajectory monitoring and intelligent deviation control device 200 can implement a real-time borehole trajectory monitoring and intelligent deviation control method according to the above-described method embodiments. The options in the above method embodiments are also applicable to this embodiment, and will not be detailed here. The remaining content of this application embodiment can be referred to the content of the above method embodiments, and will not be repeated in this embodiment.
[0100] like Figure 3As shown, the third embodiment of the present invention also provides a computer device, including a memory 301, a processor 302, and a computer program stored in the memory 301 and executable on the processor 302. When the processor 302 executes the program, it implements the steps of the drilling trajectory real-time monitoring and intelligent correction control method described in the first embodiment of the present invention.
[0101] The above embodiments are not an exhaustive list based on the present invention, and there may be many other embodiments not listed. Any substitutions and improvements made without departing from the concept of the present invention are within the protection scope of the present invention.
Claims
1. A method for real-time monitoring and intelligent deviation control of a drilling trajectory, characterized in that, The method comprises the following steps: Collecting inertial navigation data and geomagnetic sensor data during drilling, performing distortion analysis on the geomagnetic sensor data to obtain magnetic field disturbance characteristics, and using the magnetic field disturbance characteristics to perform waveform compensation fusion to generate a compensated drill bit attitude; Extracting a trajectory deviation vector based on the compensated drill bit attitude, decomposing the trajectory deviation vector into a positive deviation component and a negative deviation component, and performing potential energy accumulation on the negative deviation component to generate a deviation energy pool; Monitoring the deviation energy pool to identify a stable boundary and an unstable boundary, determining a chaos control band between the stable boundary and the unstable boundary, and extracting an optimal response interval through the chaos control band; Determining an active correction vector based on the optimal response interval and the positive deviation component, controlling the release of the deviation energy pool to obtain an energy release amount, and performing vector superposition on the active correction vector and the energy release amount to generate a composite correction instruction; Controlling a hydraulic cylinder to generate a thrust distribution through the composite correction instruction, generating a deflection torque based on the thrust distribution acting on the drill bit, changing the drilling direction using the deflection torque to obtain an actual correction amount, and extracting oscillation characteristics of the actual correction amount to form a system oscillation parameter; Generating a stability index based on the position mapping of the system oscillation parameter and the chaos control band, superimposing the stability index and the thrust distribution to form a fine correction force, and constructing a chaos edge control field based on the fine correction force and the negative deviation component; Obtaining the conversion efficiency of the deviation energy pool and the actual correction amount, and generating a real-time correction control amount based on the chaos edge control field and the conversion efficiency.
2. The method of claim 1, wherein, The method comprises the following steps: Constructing a disturbance space-time graph based on the magnetic field disturbance characteristics; Extracting a disturbance propagation path from the disturbance space-time graph; Reversely constructing a compensation waveform along the disturbance propagation path; Fusing the compensation waveform with the inertial navigation data to generate a compensated drill bit attitude.
3. The method of claim 1, wherein, The method comprises the following steps: Mapping the negative deviation component to a potential energy gradient field; Finding an energy convergence point in the potential energy gradient field; Performing energy concentration processing based on the potential energy gradient field and the energy convergence point to form concentrated potential energy; Cumulatively processing the concentrated potential energy to obtain a deviation energy pool.
4. The method of claim 1, wherein, The method comprises the following steps: Identifying an energy-rich area from the chaos control band; Obtaining a response sensitivity coefficient based on the energy-rich area; Performing interval optimization processing based on the response sensitivity coefficient to form a candidate response interval; Selecting the interval with the highest stability from the candidate response interval as the optimal response interval.
5. The method of claim 1, wherein, The method comprises the following steps: Constructing an energy release time sequence based on the energy release amount; Embedding the active correction vector into the energy release time sequence to form a time-varying vector sequence; Performing envelope modulation on the time-varying vector sequence to output a modulation vector; Converting the modulation vector into a composite correction instruction.
6. The method of claim 1, wherein, The stability index is generated based on position mapping of the system oscillation parameter and the chaos control band, and the stability index includes: The system oscillation parameter is decomposed into a periodic oscillation component and a random oscillation component; The periodic oscillation component is orbitally tracked in the chaos control band to form a periodic orbit; The random oscillation component is subjected to probability distribution analysis to determine a diffusion boundary; The stability index is generated based on the relative position relationship between the periodic orbit and the diffusion boundary.
7. The method of claim 1, wherein, The chaos edge control field is constructed based on the fine correction force and the negative deviation component, and the chaos edge control field includes: An action center of the fine correction force is obtained; A field strength gradient distribution is generated based on the action center; The negative deviation component is guided along the field strength gradient distribution to perform energy redistribution to form a redistribution field; The redistribution field is subjected to edge sharpening processing to form a chaos edge control field.
8. The method of claim 4, wherein, The response sensitivity coefficient is obtained based on the energy enrichment area, and the response sensitivity coefficient includes: The energy enrichment area is decomposed into a core enrichment area and a peripheral enrichment area; A reference sensitivity is set based on the core enrichment area; The reference sensitivity is dynamically adjusted using the peripheral enrichment area to form an adjusted sensitivity; The response sensitivity coefficient is generated according to the adjusted sensitivity.
9. A device for real-time monitoring and intelligent deviation control of drilling trajectory, characterized in that, It includes: A data fusion module is used to collect inertial navigation data and geomagnetic sensor data during drilling, to obtain magnetic field disturbance characteristics by analyzing the distortion of the geomagnetic sensor data, and to generate a compensated drill bit attitude by waveform compensation fusion using the magnetic field disturbance characteristics; An energy conversion module is used to extract a trajectory deviation vector based on the compensated drill bit attitude, to decompose the trajectory deviation vector into a positive deviation component and a negative deviation component, and to generate a deviation energy pool by accumulating potential energy of the negative deviation component; A chaos analysis module is used to monitor the deviation energy pool to identify stable boundaries and unstable boundaries, to determine a chaos control band between the stable boundaries and the unstable boundaries, and to extract an optimal response interval through the chaos control band; An instruction generation module is used to determine an active correction vector according to the optimal response interval and the positive deviation component, to control the release of the deviation energy pool to obtain a stored energy release amount, and to generate a composite correction instruction by vector superposition of the active correction vector and the stored energy release amount; An execution control module is used to control a hydraulic cylinder to generate a thrust distribution through the composite correction instruction, to generate a deflection torque based on the thrust distribution acting on a drill bit, to change a drilling direction using the deflection torque to obtain an actual correction amount, and to extract a system oscillation parameter from the actual correction amount; A stability adjustment module is used to generate a stability index based on position mapping of the system oscillation parameter and the chaos control band, to form a fine correction force by superposition of the stability index and the thrust distribution, and to construct a chaos edge control field based on the fine correction force and the negative deviation component; An optimization control module is used to obtain a conversion efficiency of the deviation energy pool and the actual correction amount, and to generate a real-time correction control amount based on the chaos edge control field and the conversion efficiency.
10. A computer device, comprising: A computer program product comprising a computer readable medium, the computer readable medium having stored thereon the computer program of claim 9. A computer program comprising program code adapted to perform the method of any one of claims 1 to 8 when the program is executed on a computer. A computer program comprising program code adapted to perform the method of any one of claims 1 to 8 when the program is executed on a computer. A computer program comprising program code adapted to perform the method of any one of claims 1 to 8 when the program is executed on a computer. A computer
Citation Information
Patent Citations
Method for controlling well track by rotary steering tool
CN104453713A
Drilling machine positioning and adjusting device
CN118088038A