Drilling trajectory three-dimensional visualization simulation system

Abstract

The application discloses a drilling trajectory three-dimensional visualization simulation system, comprising the following modules: a multi-source data fusion and probabilistic module, which fuses multi-source data to generate a geological feature probability atlas; a local dynamic geological reconstruction module, which constructs an improved NeRF network and introduces a local orthogonal coordinate system decoupling mechanism to reconstruct a local three-dimensional geological model; a trajectory mechanics correction module, which calculates a lateral deflection force to obtain a corrected trajectory; an uncertainty branch deduction module, which deduces a three-dimensional trajectory tree based on a partially observable Markov decision process; a space risk calibration and penetrating visualization rendering module, which generates a risk boundary and highlights a visualization picture; and an interactive reverse control module, which reversely deduces drilling parameter instruction sets from interactive instructions. The application improves the drilling risk identification accuracy and regulation and control efficiency, and realizes intelligent perception and forward control of drilling risks.

Description

Technical Field

[0001] This invention relates to the fields of geological engineering and three-dimensional visualization, and in particular to a three-dimensional visualization simulation system for drilling trajectories. Background Technology

[0002] Drilling is a crucial means of acquiring underground resources and geological information, and the precise control of the drilling trajectory and the avoidance of geological risks directly affect the safety and efficiency of drilling operations. With the development of measurement-while-drilling (MWD) technology, seismic exploration technology, and rock mechanics testing technology, the geological and engineering data acquired during drilling are becoming increasingly abundant, providing important data support for real-time monitoring and adjustment of the drilling trajectory. Currently, the industry typically integrates multiple data sources such as MWD data and seismic wave data to construct a three-dimensional underground geological model, and designs and displays the drilling trajectory based on this model to assist on-site operators in identifying formation characteristics and potential risks.

[0003] However, existing 3D visualization and control technologies for drilling trajectories still have many shortcomings. First, in terms of multi-source data fusion, existing technologies typically only perform simple data overlay or deterministic spatial interpolation, ignoring the differences in signal-to-noise ratio and sampling density between different data sources. This fails to effectively quantify and separate hidden stratigraphic boundary trends and stress differences, resulting in unclear confidence levels in the fusion results and difficulty in supporting the reliability assessment of unknown strata ahead. Second, in geological model construction, existing 3D modeling methods mostly employ static overall reconstruction, failing to achieve local dynamic updates as the drill bit advances. Furthermore, conventional algorithms are mostly based on isotropic assumptions, failing to introduce a local orthogonal coordinate system decoupling mechanism to distinguish the anisotropic characteristics of stratigraphic alignment and trans-strata, leading to low accuracy in representing thin layers and abrupt stratigraphic changes.

[0004] Furthermore, most existing trajectory prediction and risk assessment methods are based on deterministic models, providing only a single predicted trajectory. They fail to utilize probabilistic branching programming to perform multi-step stochastic extrapolation of spatial attitude to address the uncertainties in the geological data ahead, and lack forward-looking prediction of multiple possible routes. Simultaneously, risk assessment is limited to distance calculation and point-based alarms, failing to generate three-dimensional risk boundaries for spatial constraints. Finally, in terms of visualization and human-computer interaction, existing systems often suffer from occlusion between geological bodies and trajectories, lacking penetrating rendering and visual warning mechanisms. Operators rely heavily on manual calculations based on experience for trajectory adjustments, unable to directly derive the required drilling parameters (pressure, rotation speed, and tool face angle) through intuitive dragging in a 3D view. Human-computer collaboration efficiency and control precision urgently need improvement.

[0005] Therefore, how to provide a three-dimensional visualization simulation system for drilling trajectories is a problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0006] One objective of this invention is to propose a three-dimensional visualization simulation system for drilling trajectories. This invention fully integrates key steps such as probabilistic fusion of multi-source data, local dynamic geological reconstruction, trajectory mechanical correction, uncertainty branch extrapolation, spatial risk calibration, penetrating visualization rendering, and interactive inverse control. It constructs a closed-loop control process for drilling while visualization, encompassing geological feature probabilistic map generation, anisotropic attribute decoupling, multi-branch attitude transfer prediction, three-dimensional risk boundary constraints, and inverse derivation of drilling parameters. This enables real-time perception, dynamic prediction, and proactive intervention of drilling trajectories in complex geological environments. This invention possesses advantages such as accurate quantification of geological feature confidence, accurate expression of formation anisotropy, comprehensive trajectory risk extrapolation coverage, and intuitive and efficient human-computer interactive control. It can significantly improve trajectory control accuracy, risk avoidance capabilities, and interactive decision-making efficiency under complex geological conditions, thereby effectively solving problems in existing methods such as distorted isotropic assumptions in geological models, limited prediction of forward uncertainties, and lack of intuitive inverse feedback for trajectory adjustments.

[0007] The drilling trajectory three-dimensional visualization simulation system according to an embodiment of the present invention includes the following modules:

[0008] The multi-source data fusion and probabilistic module is used to preprocess the collected drilling measurement data, seismic wave data and rock mechanics parameters to obtain a probabilistic map of geological features.

[0009] The local dynamic geological reconstruction module is used to determine the influence radius, estimate the morphology and properties of unknown strata ahead, correct the strata morphology, construct an improved NeRF network, introduce a local orthogonal coordinate system decoupling mechanism, establish a local orthogonal coordinate system, calculate anisotropic coding, spatial density and strata property values, perform volume rendering reconstruction and isosurface extraction, and obtain a local three-dimensional geological model.

[0010] The trajectory mechanics correction module is used to calculate the lateral deflection force generated when the drill bit breaks rock based on the current drilling data and the local three-dimensional geological model, and to correct the true spatial attitude of the current drilling trajectory to obtain the current corrected trajectory.

[0011] The uncertainty branching inference module is used to address the uncertainty of the geological data ahead in the geological feature probability map. It adopts a branching programming algorithm based on a partially observable Markov decision process to perform random sampling and forward inference of the drill bit's multi-step spatial attitude transfer, expand and generate multiple travel routes, and obtain a three-dimensional trajectory tree.

[0012] The spatial risk calibration module is used to calculate the spatial distance between the three-dimensional trajectory tree and the surrounding faults, high-pressure layers and adjacent wells in the local three-dimensional geological model, and to generate risk boundaries.

[0013] The penetrating visualization rendering module is used to display local 3D geological models in a semi-transparent form, display 3D trajectory trees and risk boundaries in a solid line and surface form, and highlight the trajectory segments of the 3D trajectory tree that cross the risk boundary to obtain a 3D visualization of geology and trajectory.

[0014] The interactive reverse control module receives the operator's action commands to drag and select the predicted branch in the 3D trajectory tree or modify the risk boundary in the geological-trajectory 3D visualization screen, and converts them into target spatial coordinates and target trajectory posture. It then reverse-derives and calculates the magnitude and direction of the mechanical force required for the drill bit to reach the target trajectory posture, and obtains the drilling parameter command set for drilling pressure, rotation speed and tool face angle.

[0015] Optionally, the multi-source data fusion and probabilistic module specifically includes:

[0016] The original depth and coordinate information carried by the collected drilling measurement data, seismic wave data and rock mechanics parameters are read, and coordinate translation and rotation transformation calculations are performed to align them to the same spatial coordinate system. The attribute values ​​of the drilling measurement data, seismic wave data and rock mechanics parameters under the unified coordinate system are extracted and normalized to complete the multi-source data fusion standardization.

[0017] In the fusion of standardized multi-source data, the difference between adjacent sampling points of seismic wave impedance is calculated along the depth. When it is greater than the preset impedance threshold, the corresponding depth coordinates are marked as the position of the wave impedance interface. The slope of adjacent sampling points of the logging curve of the measurement while drilling data is calculated along the depth. When it is greater than the preset slope threshold, the corresponding depth coordinates are marked as the position of the abrupt change point of the logging curve.

[0018] The location of the acoustic impedance interface at the same spatial position is matched and compared with the location of the abrupt change point of the logging curve. If the positions coincide, it is determined as a boundary. If they deviate, the vertical and horizontal distances between the acoustic impedance interface and the location of the abrupt change point of the logging curve are calculated. The two-dimensional spatial range of the vertical and horizontal distances is delineated. Within the corresponding two-dimensional spatial range, the acoustic impedance interface and the location of the abrupt change point of the logging curve are connected sequentially along the vertical and horizontal directions to form a broken line segment. The broken line segment is used as the deviation trajectory, and the spatial range of the broken line segment is separated and output as the formation boundary trend.

[0019] In the fusion of standardized multi-source data, the difference between the elastic modulus and Poisson's ratio in rock mechanics parameters is extracted, and the wave velocity difference gradient in seismic wave data is extracted. The difference between the elastic modulus and Poisson's ratio and the wave velocity difference gradient are multiplied and superimposed at the corresponding spatial coordinates. Based on the extreme value distribution of the multiplication and superposition calculation results, the direction of the region of abrupt stress distribution change is delineated. The direction of the region of abrupt stress distribution change is separated and the output is the difference in the direction of rock stress.

[0020] The differences in the trend of the separated stratigraphic boundaries and the direction of rock stress are assigned spatial coordinates in a unified coordinate system. The signal-to-noise ratio and sampling density values ​​of each source data at the corresponding location are statistically compared. The confidence score is calculated by weighted summation of the signal-to-noise ratio and sampling density values, and then normalized to obtain the probability weight of data reliability. The spatial coordinates are bound and combined with the probability weight of data reliability to output a probability map of geological features.

[0021] Optionally, the local dynamic geological reconstruction module specifically includes:

[0022] Using the current drill bit position as the origin of spatial coordinates, and the preset distance threshold as the radius of influence, the reliability probability weight of each spatial point within the radius of influence in the geological feature probability map is read. When the weight is lower than the preset weight threshold, the corresponding spatial point is marked as an unknown stratum ahead.

[0023] Read the current drill bit's drilling pressure, rotation speed, drill bit radius, drill bit cross-sectional area, number of drill bit cutting teeth, and gravity of the overlying rock strata. Calculate the volumetric work done by a single tooth breaking the rock and the horizontal stress, and combine them to obtain the morphology and properties of the unknown strata ahead, and continuously correct the strata morphology.

[0024] An improved NeRF network is constructed, isotropic basic feature vectors are calculated, and a local orthogonal coordinate system decoupling mechanism is introduced to establish a local orthogonal coordinate system and calculate the anisotropic encoded features after decoupling.

[0025] Anisotropic coding features are input into the density prediction branch. The result of multiplying the cross-layer direction features by a preset cross-layer weight coefficient and the result of multiplying the in-layer direction features by a preset in-layer weight coefficient are added element-wise. The spatial density value is calculated by using the Softplus activation function. In the attribute prediction branch, the anisotropic coding features and the spatial density value are concatenated in the channel dimension. After mapping through two fully connected layers, the formation attribute value is output.

[0026] Within the local area of ​​influence radius, samples are uniformly taken along the preset camera ray direction to obtain sampling points. The influence radius is multiplied by 2 to calculate the total ray penetration length. The total ray penetration length is divided by the preset number of sampling points to calculate the ray step size between adjacent sampling points. The spatial density value corresponding to each sampling point is multiplied by the ray step size to obtain the single-point transparency. The single-point transparency is subtracted from 1 to obtain the single-point transmittance. The single-point transmittance of all sampling points upstream of the ray is multiplied together to obtain the cumulative transmittance.

[0027] The cumulative transmittance is multiplied by the single-point transparency and the stratigraphic attribute value to obtain the single-point rendering attribute value. The single-point rendering attribute values ​​of all sampling points on the ray are summed to obtain the ray cumulative rendering value. The cumulative rendering values ​​of each ray are arranged and stitched according to the corresponding spatial sampling grid coordinates to generate a local three-dimensional geological attribute matrix as a local three-dimensional geological body.

[0028] The spatial density values ​​of all sampling points are compared with the preset density threshold. Sampling points with spatial density values ​​greater than or equal to the preset density threshold are retained as boundary points, while sampling points with spatial density values ​​less than the preset density threshold are set to zero. All non-zero boundary points are connected adjacently in space to generate a stratigraphic boundary surface. The local three-dimensional geological body is combined with the stratigraphic boundary surface to obtain a local three-dimensional geological model.

[0029] Optionally, the step of reading the current drill bit's drilling pressure, rotation speed, drill bit radius, drill bit cross-sectional area, number of drill bit cutting teeth, and overlying rock gravity, calculating the single-tooth rock-breaking volumetric work and horizontal geostress, and combining these calculations to obtain the morphology and properties of the unknown strata ahead, and continuously correcting the strata morphology, specifically includes:

[0030] Read the current drill bit's drilling pressure, rotation speed, drill bit radius, drill bit cross-sectional area, number of drill bit cutting teeth, and gravity of the overlying rock strata. Multiply the rotation speed by the drill bit radius to calculate the drill bit cutting linear velocity, divide the drilling pressure by the drill bit cross-sectional area to calculate the drill bit contact pressure, and multiply the cutting linear velocity by the contact pressure and the cutting tooth width to calculate the single-tooth rock-breaking volumetric work.

[0031] Multiply the volumetric work of single-tooth rock breaking by the number of drill bit cutting teeth and add the preset initial cohesion of the rock to calculate the specific work of rock breaking of the unknown strata ahead. Add the specific work of rock breaking to the gravity of the overlying strata to calculate the equivalent density of the strata. Multiply the density by the gravitational acceleration and the depth of the corresponding spatial point to calculate the vertical stress. Multiply the vertical stress by the preset horizontal stress coefficient to calculate the horizontal stress.

[0032] The horizontal geostress is mapped to the spatial surface undulation of the strata, and the equivalent density of the strata and the rock fracture work are mapped to the rock porosity and permeability properties of the strata. The morphology and properties of the unknown strata ahead are obtained. When the probability weight of the data reliability is not lower than the preset weight threshold, the spatial point is marked as the drilled strata behind. The measured strata morphology and property data obtained by drilling measurement are directly read to replace the calculation results and complete the strata morphology correction.

[0033] Optionally, the construction of the improved NeRF network, the calculation of isotropic fundamental feature vectors, and the introduction of a local orthogonal coordinate system decoupling mechanism to establish a local orthogonal coordinate system and calculate the decoupled anisotropic encoded features, specifically include:

[0034] An improved NeRF network is constructed, which concatenates the spatial coordinates and spatial location coordinates within the influence radius as the query input, and maps the query input to an isotropic basic feature vector of a preset dimension through an initial fully connected layer;

[0035] A local orthogonal coordinate system decoupling mechanism is introduced. The isotropic basic feature vector is input into the newly added normal vector prediction branch in the improved NeRF network. After a linear mapping layer, a 3D orientation initial vector is output. The Euclidean norm of the 3D orientation initial vector is calculated as the modulus. Each component of the 3D orientation initial vector is divided by the modulus to obtain the unit normal vector. The unit normal vector is used as the Z-axis of the local orthogonal coordinate system to characterize the stratigraphic cross-layer direction.

[0036] On a plane perpendicular to the unit normal vector, define an arbitrary auxiliary vector that is not parallel to the unit normal vector. Perform a cross product operation between the unit normal vector and the arbitrary auxiliary vector and divide by the magnitude of the cross product result to obtain the first orthogonal unit vector, which serves as the X-axis of the local orthogonal coordinate system, representing the strike of the strata. Perform a cross product operation between the unit normal vector and the first orthogonal unit vector to obtain the second orthogonal unit vector, which serves as the Y-axis of the local orthogonal coordinate system, representing the dip of the strata. The local orthogonal coordinate system is established by combining the first orthogonal unit vector, the second orthogonal unit vector, and the unit normal vector.

[0037] The spatial coordinates within the influence radius are multiplied by the first orthogonal unit vector, the second orthogonal unit vector, and the unit normal vector respectively to obtain the local X coordinate, local Y coordinate, and local Z coordinate in the local orthogonal coordinate system. These coordinates are then added together to form the in-layer coordinates, and the local Z coordinate is used as the trans-layer coordinate.

[0038] The isotropic basic feature vector is concatenated with the in-layer and cross-layer coordinates and input into the newly added anisotropic coding layer in the improved NeRF network. In the in-layer direction, the in-layer coordinates are input into the sine and cosine functions with low-frequency preset parameters for periodic mapping. The mapping results are concatenated and multiplied by the first set of learnable weight matrices to extract the in-layer features.

[0039] In the layer-penetrating direction, the layer-penetrating direction coordinates are respectively input into a sine function and a cosine function with high-frequency preset parameters for periodic mapping. The mapping results are concatenated and multiplied by the second set of learnable weight matrices to extract the layer-penetrating direction features. The in-layer direction features and the layer-penetrating direction features are added element by element to obtain the decoupled anisotropic coding features.

[0040] Optionally, the trajectory mechanics correction module specifically includes:

[0041] Read the drill bit lateral force and drill bit mechanical torque from the current drilling data, extract the formation equivalent density and rock breaking specific work corresponding to the current spatial position of the drill bit from the local three-dimensional geological model, and calculate the rock hardness index by dividing the rock breaking specific work by the formation equivalent density.

[0042] Along the direction of drill bit travel, extract the equivalent density of the formation at a preset distance to the left and a preset distance to the right of the drill bit in the local three-dimensional geological model. Divide the equivalent density of the formation on the left by the equivalent density of the formation on the right to calculate the density ratio of the formation on the left and right. Along the direction of drill bit travel, extract the equivalent density of the formation at a preset distance in front of the drill bit and a preset distance below the drill bit. Divide the equivalent density of the formation in front by the equivalent density of the formation below to calculate the density ratio of the formation in front and below.

[0043] The lateral formation deflection force is calculated by subtracting 1 from the density ratio of the left and right formations and multiplying it by the lateral force of the drill bit. The longitudinal formation deflection force is calculated by subtracting 1 from the density ratio of the front and rear formations and multiplying it by the mechanical torque of the drill bit. The lateral rock-breaking deflection force is calculated by multiplying the rock hardness index by the lateral formation deflection force. The longitudinal rock-breaking deflection force is calculated by multiplying the rock hardness index by the longitudinal formation deflection force.

[0044] The lateral rock-breaking deflection force is calculated by vector addition of the lateral rock-breaking deflection force and the longitudinal rock-breaking deflection force. The reference mechanical vector is constructed by combining the drill bit lateral force and drill bit mechanical torque. The trajectory deflection coefficient is obtained by ratio calculation of the lateral deflection force and the reference mechanical vector. The spatial attitude offset is calculated by multiplying the trajectory deflection coefficient by the three-dimensional direction vector of the current drilling trajectory.

[0045] The corrected three-dimensional direction vector is calculated by subtracting the spatial attitude offset from the three-dimensional direction vector of the current drilling trajectory. The corrected three-dimensional direction vector is then combined with the coordinate position of the current drilling trajectory to obtain the current corrected trajectory.

[0046] Optionally, the uncertainty branch deduction module specifically includes:

[0047] The probability weights of data reliability of each spatial point within the influence radius in front of the current drill bit are read from the geological feature probability map. Spatial points with weights below the preset threshold are extracted as uncertain spatial points. The stratigraphic attribute values ​​of the uncertain spatial points are combined with the probability weights of data reliability to construct a discrete probability distribution table.

[0048] Based on the partially observable Markov decision process, the number of state transition steps is set. In each state transition step, the stratigraphic attribute values ​​of uncertain spatial points are randomly sampled according to the discrete probability distribution table. The stratigraphic attribute values ​​obtained from each random sampling are used as the simulated stratigraphic attributes of the corresponding state transition step.

[0049] Read the drill pressure, rotation speed and tool face angle from the current drilling data, combine the three-dimensional direction vector of the current correction trajectory with the simulated formation properties to calculate the rock hardness index and formation deflection force, divide the formation deflection force by the product of drill pressure and rotation speed to calculate the attitude deflection angle, add the attitude deflection angle to the tool face angle to calculate the resultant force deflection direction, and add the resultant force deflection direction to the three-dimensional direction vector of the current correction trajectory to calculate the next predicted direction vector;

[0050] The coordinate offset is calculated by multiplying the next predicted direction vector by the preset single-step advance length, and the next predicted coordinate position is calculated by adding the current spatial position coordinates. The next predicted direction vector and the next predicted coordinate position are combined to form the single-step spatial attitude transfer result.

[0051] Repeatedly perform random sampling and spatial attitude transfer calculations until the preset number of state transfer steps is reached. Connect the single-step spatial attitude transfer results corresponding to each random sampling to generate a future prediction branch. Summarize all future prediction branches generated by random sampling to expand and generate multiple travel routes. Patch all future prediction branches with the current corrected trajectory at the starting end to obtain a three-dimensional trajectory tree.

[0052] Optionally, the space risk calibration module specifically includes:

[0053] Spatial points whose stratigraphic attribute values ​​reach the preset fault safety threshold are extracted from the local 3D geological model and connected to generate surrounding faults. Spatial points whose stratigraphic equivalent density is greater than the preset high pressure safety threshold are extracted and connected to generate high pressure layers. Pre-stored adjacent well trajectory coordinate data are read to generate adjacent wells.

[0054] The Euclidean distance between each trajectory sampling point in the 3D trajectory tree and the surrounding fault is calculated as the fault distance; the Euclidean distance between each trajectory sampling point in the 3D trajectory tree and the high-pressure layer is calculated as the high-pressure distance; and the Euclidean distance between each trajectory sampling point in the 3D trajectory tree and the adjacent well is calculated as the adjacent well distance.

[0055] The fault distance, high pressure distance, and adjacent well distance are compared with their corresponding safety thresholds. When the fault distance is less than the preset fault safety threshold, the trajectory sampling point is marked as a fault risk point. When the high pressure distance is less than the preset high pressure safety threshold, the trajectory sampling point is marked as a high pressure risk point. When the adjacent well distance is less than the preset anti-collision safety threshold, the trajectory sampling point is marked as an anti-collision risk point.

[0056] All fault risk points, high-pressure risk points, and collision avoidance risk points are merged and collectively referred to as risk space points. The Euclidean distance between adjacent risk space points is calculated. Adjacent risk space points with an Euclidean distance less than a preset connection distance are connected to each other. A spatial mesh surface is constructed based on all the connections, and the surface is filled with solids to generate a three-dimensional risk boundary.

[0057] Optionally, the permeable visualization rendering module specifically includes:

[0058] Read the local 3D geological model, multiply the opacity parameter of all voxels in the local 3D geological model by the preset semi-transparency coefficient to calculate the semi-transparency opacity, replace the original opacity with the semi-transparency opacity for ray casting rendering, and display the local 3D geological model in a semi-transparent form.

[0059] Read the 3D trajectory tree and risk boundary, connect the coordinates of each trajectory sampling point in the 3D trajectory tree with preset line segments to generate solid trajectory lines, and stitch together all surface space points in the risk boundary with triangular facets to generate solid boundary surfaces. Rasterize and color the solid trajectory lines and solid boundary surfaces, and display the 3D trajectory tree and risk boundary in solid line surface form.

[0060] Calculate the shortest Euclidean distance between each trajectory line in the 3D trajectory tree and the risk boundary surface. When the shortest Euclidean distance is equal to zero, mark the corresponding trajectory line as an intersecting trajectory segment and extract the coordinates of the two intersection points of the intersecting trajectory segment entering and exiting the risk boundary.

[0061] Extract the intersecting trajectory segment between the coordinates of two intersection points as the trajectory segment that crosses the risk boundary, replace the original rendering color of the trajectory segment that crosses the risk boundary with the preset highlight color, multiply the luminous intensity of the preset highlight color by the sine time function to calculate the dynamic flicker brightness value, and perform light highlighting processing on the trajectory segment that crosses the risk boundary according to the dynamic flicker brightness value.

[0062] The semi-transparent local 3D geological model, the solid line surface 3D trajectory tree and risk boundary, and the trajectory segment passing through the risk boundary with dynamic flashing brightness value light highlighting are deeply buffered and fused in the same 3D spatial coordinate system to obtain a geological-trajectory 3D visualization.

[0063] Optionally, the interactive reverse control module specifically includes:

[0064] The system receives action commands from operators to drag and select predicted branches or modify risk boundaries in the 3D trajectory tree in the geological-trajectory 3D visualization screen, extracts the screen pixel coordinates corresponding to the action commands, multiplies the screen pixel coordinates by the preset camera inverse projection matrix to calculate the normalized ray direction, multiplies the normalized ray direction by the preset depth value and adds the camera spatial coordinates to calculate the target spatial coordinates.

[0065] Construct a tangent vector at the target space coordinates along the extension direction dragged by the operator, calculate the Euclidean norm of the tangent vector as the modulus, divide the tangent vector by the modulus to calculate the unit direction vector, and combine the unit direction vector with the target space coordinates to obtain the target trajectory attitude.

[0066] The angle between the unit direction vector in the target trajectory attitude and the three-dimensional direction vector of the current corrected trajectory is calculated as the attitude deflection angle. The sine and cosine values ​​of the attitude deflection angle are calculated. The mechanical torque and lateral force of the drill bit are decomposed by combining the sine and cosine values ​​of the attitude deflection angle to obtain the required mechanical force in the lateral and longitudinal directions. The required mechanical force in the lateral and longitudinal directions are vector-added to calculate the magnitude and direction of the mechanical force.

[0067] The drill pressure release factor is calculated by dividing the magnitude of the applied mechanical force by the vector sum of the drill bit lateral force and the drill bit mechanical torque. The target drill pressure is calculated by multiplying the drill pressure in the current drilling data by the drill pressure release factor. The target rotational speed is calculated by multiplying the rotational speed in the current drilling data by the drill pressure release factor.

[0068] The tool face deflection angle is calculated by projecting the mechanical force direction onto a plane perpendicular to the three-dimensional direction vector of the current correction trajectory. The tool face deflection angle is then added to the tool face angle in the current drilling data to calculate the target tool face angle. The drilling parameter command set of drilling pressure, rotational speed and tool face angle is obtained by combining the target drilling pressure, target rotational speed and target tool face angle.

[0069] The beneficial effects of this invention are:

[0070] This invention addresses the issues of heterogeneity and uncertain reliability of multi-source data during drilling through a multi-source data fusion and probabilistic module. It employs a signal-to-noise ratio and sampling density weighted calculation of confidence levels, combined with cross-comparison of acoustic impedance and logging curves, to generate a geological feature probability map with probability weights reflecting the reliability of labeled data. Furthermore, through a local dynamic geological reconstruction module, it addresses the problems of lagging geological model updates and distortion of the isotropic assumption. It utilizes an improved NeRF network and introduces a local orthogonal coordinate system decoupling mechanism to decouple input features into anisotropic encodings along and across bedding planes before volume rendering reconstruction and isosurface extraction, resulting in a locally three-dimensional geological model that grows and updates in real-time during drilling. The trajectory mechanics correction module calculates the lateral deflection force to obtain the current corrected trajectory. Combined with the uncertainty branching deduction module, and addressing geological uncertainties ahead, a branch programming algorithm based on a partially observable Markov decision process is used to randomly sample formation attributes and perform forward deduction of multi-step spatial attitude transitions, generating a 3D trajectory tree containing multiple predicted branches. Further, a spatial risk calibration module generates a 3D risk boundary, and a penetration-based visualization rendering module highlights the trajectory segments crossing the risk boundary. An interactive reverse control module then uses the operator's dragging actions to derive and calculate the drilling parameter command set, including drilling pressure, rotational speed, and tool face angle. Ultimately, this achieves closed-loop intelligent monitoring of geological dynamic reconstruction, uncertainty risk deduction, and parameter reverse control during drilling, effectively improving the accuracy of geological attribute representation, the timeliness of trajectory risk warnings, and the execution efficiency of human-machine interactive control. Attached Figure Description

[0071] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings:

[0072] Figure 1 This is a structural diagram of the three-dimensional visualization simulation system for drilling trajectories proposed in this invention;

[0073] Figure 2 This is a flowchart of the improved NeRF network geological reconstruction and attribute prediction based on the local orthogonal coordinate system decoupling mechanism proposed in this invention;

[0074] Figure 3 This is a flowchart of the three-dimensional trajectory tree deduction and inverse control based on a partially observable Markov decision process proposed in this invention. Detailed Implementation

[0075] The present invention will now be described in further detail with reference to the accompanying drawings. These drawings are simplified schematic diagrams, illustrating only the basic structure of the invention, and therefore only show the components relevant to the invention.

[0076] refer to Figures 1-3 The drilling trajectory 3D visualization simulation system includes the following modules:

[0077] The multi-source data fusion and probabilistic module is used to perform spatial mapping and numerical standardization of the collected drilling measurement data, seismic wave data and rock mechanical parameters in a unified coordinate system. Through multi-dimensional data cross-comparison, it separates the hidden stratigraphic boundary trends and differences in rock stress direction, and obtains a geological feature probability map with spatial location coordinates and marked with the probability weight of data reliability.

[0078] The local dynamic geological reconstruction module is used to determine the influence radius with the current drill bit position as the calculation center. Within the influence radius in front of the drill bit, the morphology and properties of the unknown strata ahead are calculated using mechanical formulas. Behind the drill bit, the strata morphology is corrected using the measured data obtained from drilling. An improved NeRF network is constructed, and a local orthogonal coordinate system decoupling mechanism is introduced. The spatial coordinates within the influence radius are used as network input. The strata normal vectors of the predicted spatial points are used to establish a local orthogonal coordinate system. The input features are decoupled into anisotropic encodings of the bedding direction and the cross-bedding direction, and then the spatial density and strata property values ​​are calculated. Only the local area within the influence radius is subjected to volume rendering reconstruction and isosurface extraction to obtain a local three-dimensional geological model that grows and updates in real time as the drill bit advances.

[0079] The trajectory mechanics correction module is used to calculate the lateral deflection force generated by the drill bit when breaking rocks in rocks of different hardness based on the current drilling data and local three-dimensional geological model, and to correct the true spatial attitude of the current drilling trajectory to obtain the current corrected trajectory.

[0080] The uncertainty branch extrapolation module is used to address the uncertainty of the geological data ahead in the geological feature probability map. It adopts a branch programming algorithm based on a partially observable Markov decision process. During the state transition process, it randomly samples the possible geological attributes ahead, simulates various geological changes, and combines drilling mechanics constraints to forward extrapolate the multi-step spatial attitude transition of the drill bit. It expands and generates multiple travel routes, and combines them with the current corrected trajectory to obtain a three-dimensional trajectory tree that includes the current corrected trajectory and multiple future predicted branches.

[0081] The spatial risk calibration module is used to calculate the spatial distance between the three-dimensional trajectory tree and the surrounding faults, high-pressure layers and adjacent wells in the local three-dimensional geological model, and to connect areas with distances below the safety threshold to generate a three-dimensional risk boundary.

[0082] The penetrating visualization rendering module is used to display local 3D geological models in a semi-transparent form, display 3D trajectory trees and risk boundaries in a solid line and surface form, and highlight the trajectory segments of the 3D trajectory tree that cross the risk boundary, resulting in a geological-trajectory 3D visualization screen with a sense of penetration and visual warnings for risk areas.

[0083] The interactive reverse control module is used to receive the operator's action commands to drag and select the predicted branch in the 3D trajectory tree or modify the risk boundary in the geological-trajectory 3D visualization screen, convert the action commands into target spatial coordinates and target trajectory attitude, reverse derive and calculate the magnitude and direction of the mechanical force required for the drill bit to reach the target trajectory attitude, and obtain the drilling parameter command set of drilling pressure, rotation speed and tool face angle required to achieve the target trajectory attitude.

[0084] This invention significantly improves the accuracy of risk identification and trajectory control efficiency during drilling in complex formations. Through multi-source data fusion and probabilistic modeling, it achieves unified spatial mapping of heterogeneous data during drilling, effectively extracting hidden formation boundary trends and quantifying data reliability, thus enhancing the ability to perceive unknown uncertainties underground. A local orthogonal coordinate system decoupling mechanism is introduced into local dynamic geological reconstruction to accurately separate the anisotropic features of bedding planes and cross-bedding planes, enabling high-precision real-time growth of local 3D geological models during drilling. Based on geological probability maps, a partially observable Markov decision process is used to randomly sample and mechanically extrapolate uncertainties ahead, proactively generating a multi-branch 3D trajectory tree, transforming traditional passive risk response into active probabilistic extrapolation. Combining spatial risk calibration and penetrating visualization rendering, it achieves three-dimensional constraints on risk boundaries and intuitive high-brightness warnings of dangerous trajectories. Through interactive reverse control, the operator's dragging intentions are instantly decomposed into a set of instructions for drill pressure, rotation speed, and tool face angle, completely establishing a closed loop between 3D visualization interaction and bottom-level drilling execution. This invention demonstrates strong robustness and practicality in complex geological environments of deep wells, significantly improving the lead time for geological risk warnings during drilling, the smoothness of trajectory control, and the level of intelligence in engineering intervention.

[0085] In this embodiment, the multi-source data fusion and probabilistic approach module specifically includes:

[0086] The original depth and coordinate information carried by the collected drilling measurement data, seismic wave data and rock mechanics parameters are read, and coordinate translation and rotation transformation calculations are performed to align them to the same spatial coordinate system. The attribute values ​​of the drilling measurement data, seismic wave data and rock mechanics parameters under the unified coordinate system are extracted and normalized to complete the multi-source data fusion standardization.

[0087] In the standardized multi-source data fusion, the difference in wave impedance between adjacent sampling points of seismic wave data is calculated along the depth. When it is greater than a preset impedance threshold, the corresponding depth coordinates are marked as the wave impedance interface position. The slope of the logging curve between adjacent sampling points of the drilling measurement data is calculated along the depth. When it is greater than a preset slope threshold, the corresponding depth coordinates are marked as the location of the logging curve abrupt change point. The preset impedance threshold is 1.5 times the average absolute value of the historical wave impedance difference, and the preset slope threshold is 2 times the average absolute value of the historical slope of the logging curve.

[0088] The location of the acoustic impedance interface at the same spatial position is matched and compared with the location of the abrupt change point of the logging curve. If the positions coincide, it is determined as a boundary. If they deviate, the vertical and horizontal distances between the acoustic impedance interface and the location of the abrupt change point of the logging curve are calculated. The two-dimensional spatial range of the vertical and horizontal distances is delineated. Within the corresponding two-dimensional spatial range, the acoustic impedance interface and the location of the abrupt change point of the logging curve are connected sequentially along the vertical and horizontal directions to form a broken line segment. The broken line segment is used as the deviation trajectory, and the spatial range of the broken line segment is separated and output as the formation boundary trend.

[0089] In the fusion of standardized multi-source data, the difference between the elastic modulus and Poisson's ratio in different directions is extracted from the rock mechanics parameters, and the wave velocity difference gradient representing the change of confining pressure is extracted from the seismic wave data. The difference between the elastic modulus and Poisson's ratio and the wave velocity difference gradient are multiplied and superimposed at the corresponding spatial coordinates. Based on the extreme value distribution of the product superposition calculation result, the direction of the region of abrupt stress distribution change is outlined. The direction of the region of abrupt stress distribution change is separated and the output is the difference in the direction of rock stress.

[0090] The differences in the trend of the separated stratigraphic boundaries and the direction of rock stress are assigned spatial coordinates in a unified coordinate system. The signal-to-noise ratio and sampling density values ​​of each source data at the corresponding location are statistically compared. The confidence score is calculated by weighted summation of the signal-to-noise ratio and sampling density values, and then normalized to obtain the probability weight of data reliability. The spatial coordinates are bound and combined with the probability weight of data reliability to output a probability map of geological features.

[0091] In this embodiment, the local dynamic geological reconstruction module specifically includes:

[0092] Using the current drill bit position as the origin of spatial coordinates, and a preset distance threshold as the radius of influence, the reliability probability weight of data for each spatial point within the radius of influence in the geological feature probability map is read. When the data weight is lower than the preset weight threshold, the corresponding spatial point is marked as an unknown stratum ahead. The preset distance threshold is 50 meters, and the preset weight threshold is 0.5.

[0093] Read the current drill bit's drilling pressure, rotation speed, drill bit radius, drill bit cross-sectional area, number of drill bit cutting teeth, and gravity of the overlying rock strata. Calculate the volumetric work done by a single tooth breaking the rock and the horizontal stress, and combine them to obtain the morphology and properties of the unknown strata ahead, and continuously correct the strata morphology.

[0094] An improved NeRF network is constructed, isotropic basic feature vectors are calculated, and a local orthogonal coordinate system decoupling mechanism is introduced to establish a local orthogonal coordinate system and calculate the anisotropic encoded features after decoupling.

[0095] Anisotropic coding features are input into the density prediction branch. The result of multiplying the cross-layer direction features by a preset cross-layer weight coefficient and the result of multiplying the in-layer direction features by a preset in-layer weight coefficient are added element-wise. The spatial density value is calculated by using the Softplus activation function. In the attribute prediction branch, the anisotropic coding features and the spatial density value are concatenated in the channel dimension. After mapping through two fully connected layers of dimension 256, the stratigraphic attribute value is output. The preset cross-layer weight coefficient is 1.5, and the preset in-layer weight coefficient is 0.5.

[0096] Within the local area of ​​the influence radius, sampling is uniformly performed along the preset camera ray direction to obtain sampling points. The influence radius is multiplied by 2 to calculate the total ray penetration length. The total ray penetration length is divided by the preset number of sampling points to calculate the ray step size between adjacent sampling points. The spatial density value corresponding to each sampling point is multiplied by the ray step size to obtain the single-point transparency. The single-point transparency is subtracted from 1 to obtain the single-point transmittance. The single-point transmittance of all sampling points upstream of the ray is multiplied together to obtain the cumulative transmittance. The preset camera ray direction is a cluster of scattering rays parallel to the wellbore extension direction and with the drill bit as the vertex. The preset number of sampling points is 64.

[0097] The cumulative transmittance is multiplied by the single-point transparency and the stratigraphic attribute value to obtain the single-point rendering attribute value. The single-point rendering attribute values ​​of all sampling points on the ray are summed to obtain the ray cumulative rendering value. The cumulative rendering values ​​of each ray are arranged and stitched according to the corresponding spatial sampling grid coordinates to generate a local three-dimensional geological attribute matrix as a local three-dimensional geological body.

[0098] The spatial density values ​​of all sampling points are compared with a preset density threshold. Sampling points with spatial density values ​​greater than or equal to the preset density threshold are retained as boundary points, while sampling points with spatial density values ​​less than the preset density threshold are set to zero. All non-zero boundary points are connected adjacently in space to generate a stratigraphic boundary surface. The local three-dimensional geological body is combined with the stratigraphic boundary surface to obtain a local three-dimensional geological model that grows and updates in real time as the drill bit advances. The preset density threshold is 0.5.

[0099] In this embodiment, the drill bit's current drill bit pressure, rotation speed, drill bit radius, drill bit cross-sectional area, number of cutting teeth, and overlying rock gravity are read. The volumetric work done by a single tooth breaking the rock and the horizontal stress are calculated and combined to obtain the morphology and properties of the unknown strata ahead. Continuous strata morphology correction is then performed, specifically including:

[0100] Read the current drill bit's drilling pressure, rotation speed, drill bit radius, drill bit cross-sectional area, number of drill bit cutting teeth, and gravity of the overlying rock strata. Multiply the rotation speed by the drill bit radius to calculate the drill bit cutting linear velocity, divide the drilling pressure by the drill bit cross-sectional area to calculate the drill bit contact pressure, and multiply the cutting linear velocity by the contact pressure and the cutting tooth width to calculate the single-tooth rock-breaking volumetric work.

[0101] Multiply the volumetric work of single-tooth rock breaking by the number of drill bit cutting teeth and add the preset initial cohesion of the rock to calculate the specific work of rock breaking in the unknown strata ahead. Add the specific work of rock breaking to the gravity of the overlying strata to calculate the equivalent density of the strata. Multiply the density by the gravitational acceleration and the depth of the corresponding spatial point to calculate the vertical stress. Multiply the vertical stress by the preset horizontal stress coefficient to calculate the horizontal stress, where the horizontal stress coefficient is 0.8 to 1.2.

[0102] The horizontal geostress is mapped to the spatial surface undulation of the strata, and the equivalent density of the strata and the rock fracture work are mapped to the rock porosity and permeability properties of the strata. The morphology and properties of the unknown strata ahead are obtained. When the probability weight of the data reliability is not lower than the preset weight threshold, the spatial point is marked as the drilled strata behind. The measured strata morphology and property data obtained by drilling measurement are directly read to replace the calculation results and complete the strata morphology correction.

[0103] In this embodiment, an improved NeRF network is constructed, isotropic fundamental feature vectors are calculated, and a local orthogonal coordinate system decoupling mechanism is introduced to establish a local orthogonal coordinate system. The decoupled anisotropic encoded features are then calculated and output, specifically including:

[0104] An improved NeRF network is constructed, which concatenates the spatial coordinates and spatial location coordinates within the influence radius as the query input. The query input is then mapped to an isotropic basic feature vector of a preset dimension of 256 through an initial fully connected layer.

[0105] A local orthogonal coordinate system decoupling mechanism is introduced. The isotropic basic feature vector is input into the newly added normal vector prediction branch in the improved NeRF network. After a linear mapping layer, a 3D orientation initial vector is output. The Euclidean norm of the 3D orientation initial vector is calculated as the modulus. Each component of the 3D orientation initial vector is divided by the modulus to obtain the unit normal vector. The unit normal vector is used as the Z-axis of the local orthogonal coordinate system to characterize the stratigraphic cross-layer direction.

[0106] On a plane perpendicular to the unit normal vector, define an arbitrary auxiliary vector that is not parallel to the unit normal vector. Perform a cross product operation between the unit normal vector and the arbitrary auxiliary vector and divide by the magnitude of the cross product result to obtain the first orthogonal unit vector, which serves as the X-axis of the local orthogonal coordinate system, representing the strike of the strata. Perform a cross product operation between the unit normal vector and the first orthogonal unit vector to obtain the second orthogonal unit vector, which serves as the Y-axis of the local orthogonal coordinate system, representing the dip of the strata. The local orthogonal coordinate system is established by combining the first orthogonal unit vector, the second orthogonal unit vector, and the unit normal vector.

[0107] The spatial coordinates within the influence radius are multiplied by the first orthogonal unit vector, the second orthogonal unit vector, and the unit normal vector respectively to obtain the local X coordinate, local Y coordinate, and local Z coordinate in the local orthogonal coordinate system. These coordinates are then added together to form the in-layer coordinates, and the local Z coordinate is used as the trans-layer coordinate.

[0108] The isotropic basic feature vector is concatenated with the in-layer and cross-layer coordinates and input into the newly added anisotropic coding layer in the improved NeRF network. In the in-layer direction, the in-layer coordinates are respectively input into the sine and cosine functions with low-frequency preset parameters for periodic mapping. The mapping results are concatenated and multiplied by the first set of learnable weight matrices to extract the gradually changing in-layer features. The low-frequency preset parameters are 2 to the power of 0 to 2 to the power of 3, and the high-frequency preset parameters are 2 to the power of 4 to 2 to the power of 7.

[0109] In the layer-penetrating direction, the layer-penetrating direction coordinates are respectively input into the sine and cosine functions with high-frequency preset parameters for periodic mapping. The mapping results are concatenated and multiplied by the second set of learnable weight matrices to extract the drastically changing layer-penetrating direction features. The in-layer direction features and the layer-penetrating direction features are added element by element to obtain the decoupled anisotropic coding features.

[0110] This invention achieves dynamic reconstruction and high-precision attribute prediction of unknown strata during drilling by introducing a local orthogonal coordinate system decoupling mechanism and an improved NeRF network. A 50-meter influence radius is defined centered on the drill bit, and unknown and drilled strata are distinguished based on probability weights. For unknown strata, drilling parameters are used to estimate rock fracturing specific work and horizontal in-situ stress, proactively mapping strata morphology and attributes, and continuously correcting based on measured data. The core of the invention is the construction of a local orthogonal coordinate system, predicting strata normal vectors to establish cross-bedding and bedding directions, and decoupling features through an anisotropic coding layer: low-frequency mapping is used to capture gradual changes in bedding, while high-frequency mapping is used to capture dramatic responses in cross-bedding, and differentiated weighting coefficients are assigned to density prediction. Finally, ray volume rendering and isosurface extraction generate a real-time growing local 3D geological model. This invention effectively overcomes the thin-layer distortion problem caused by the isotropic assumption, accurately characterizing anisotropic features in complex structures and deep strata, and significantly improving the boundary resolution and attribute prediction reliability of geological modeling during drilling.

[0111] The improved NeRF network of this invention is similar to the original NeRF network in that both retain the core volume rendering architecture of the neural radiation field, that is, the spatial coordinates are mapped to density and attribute values ​​through a multilayer perceptron, the transparency and transmittance of a single point are calculated by sampling along the ray, and the final rendering result is generated by cumulative summation based on the volume rendering equation. Both also use volume rendering reconstruction and isosurface extraction to generate a three-dimensional model.

[0112] The difference lies in that this invention breaks away from the limitations of the original NeRF network, which uses isotropic feature encoding and a single static mapping for all spatial directions. It introduces a local orthogonal coordinate system decoupling mechanism to construct a direction-aware anisotropic expression system. Building upon the original model's direct input coordinates, this invention pre-decomposes spatial coordinates into in-layer and trans-layer coordinates by calculating unit normals and orthogonal vectors to construct a local coordinate system. Next, in the encoding layer, it breaks away from the original model's sharing of the same frequency mapping. For in-layer directions, low-frequency parameters are used for periodic mapping to extract smooth features, while for trans-layer directions, high-frequency parameters are used to extract dramatic features, and differentiated weight coefficients are assigned for element-wise fusion. Finally, in the prediction branch, density prediction combines differentiated weights to output spatial density values, and the attribute branch concatenates density features to complete a high-precision mapping.

[0113] Based on the above improvements, the beneficial effects of this invention are that by establishing a local orthogonal coordinate system and decoupling the direction, the improved NeRF network can adaptively distinguish between the gentle extension of the strata and the violent fluctuations across the strata, breaking the limitation of the original NeRF using the isotropic assumption to handle geological boundaries and realizing the dynamic differentiation of feature expression. This design significantly enhances the model's resolution depth for complex structures such as thin reservoirs and micro-faults, and more accurately captures abrupt changes in strata boundaries. The differentiated weights and frequency mapping not only improve the accuracy of attribute prediction but also effectively suppress the boundary ambiguity caused by isotropy, enhancing the robustness and accuracy of the model in deep and complex geological reconstruction.

[0114] In this embodiment, the trajectory mechanics correction module specifically includes:

[0115] Read the drill bit lateral force and drill bit mechanical torque from the current drilling data, extract the formation equivalent density and rock breaking specific work corresponding to the current spatial position of the drill bit from the local three-dimensional geological model, and calculate the rock hardness index by dividing the rock breaking specific work by the formation equivalent density.

[0116] Along the direction of drill bit travel, extract the equivalent density of the formation at a preset distance to the left and a preset distance to the right of the drill bit in the local three-dimensional geological model. Divide the equivalent density of the formation on the left by the equivalent density of the formation on the right to calculate the density ratio of the formation on the left and right. Extract the equivalent density of the formation at a preset distance in front of the drill bit and a preset distance below the drill bit along the direction of drill bit travel. Divide the equivalent density of the formation in front by the equivalent density of the formation below to calculate the density ratio of the formation in front and below. The preset distance is 0.5 meters.

[0117] The lateral formation deflection force is calculated by subtracting 1 from the density ratio of the left and right formations and multiplying it by the lateral force of the drill bit. The longitudinal formation deflection force is calculated by subtracting 1 from the density ratio of the front and rear formations and multiplying it by the mechanical torque of the drill bit. The lateral rock-breaking deflection force is calculated by multiplying the rock hardness index by the lateral formation deflection force. The longitudinal rock-breaking deflection force is calculated by multiplying the rock hardness index by the longitudinal formation deflection force.

[0118] The lateral rock-breaking deflection force is calculated by vector addition of the lateral rock-breaking deflection force and the longitudinal rock-breaking deflection force. The reference mechanical vector is constructed by combining the drill bit lateral force and drill bit mechanical torque. The trajectory deflection coefficient is obtained by ratio calculation of the lateral deflection force and the reference mechanical vector. The spatial attitude offset is calculated by multiplying the trajectory deflection coefficient by the three-dimensional direction vector of the current drilling trajectory.

[0119] The corrected three-dimensional direction vector is calculated by subtracting the spatial attitude offset from the three-dimensional direction vector of the current drilling trajectory. The corrected three-dimensional direction vector is then combined with the coordinate position of the current drilling trajectory to obtain the current corrected trajectory.

[0120] In this embodiment, the uncertainty branch deduction module specifically includes:

[0121] The probability weights of data reliability of each spatial point within the influence radius in front of the current drill bit are read from the geological feature probability map. Spatial points with weights below the preset threshold are extracted as uncertain spatial points. The stratigraphic attribute values ​​of the uncertain spatial points are combined with the probability weights of data reliability to construct a discrete probability distribution table.

[0122] Based on the partially observable Markov decision process, the number of state transition steps is set. In each state transition step, the stratigraphic attribute values ​​of uncertain spatial points are randomly sampled according to the discrete probability distribution table. The stratigraphic attribute values ​​obtained from each random sampling are used as the simulated stratigraphic attributes of the corresponding state transition step.

[0123] Read the drill pressure, rotation speed and tool face angle from the current drilling data, combine the three-dimensional direction vector of the current correction trajectory with the simulated formation properties to calculate the rock hardness index and formation deflection force, divide the formation deflection force by the product of drill pressure and rotation speed to calculate the attitude deflection angle, add the attitude deflection angle to the tool face angle to calculate the resultant force deflection direction, and add the resultant force deflection direction to the three-dimensional direction vector of the current correction trajectory to calculate the next predicted direction vector;

[0124] The coordinate offset is calculated by multiplying the next predicted direction vector by the preset single-step advancement length, and the next predicted coordinate position is calculated by adding the current spatial position coordinates. The next predicted direction vector and the next predicted coordinate position are combined to form the single-step spatial attitude transfer result. The preset single-step advancement length is 1 meter.

[0125] Repeatedly perform random sampling and spatial attitude transfer calculations until a preset number of state transition steps are reached. Connect the single-step spatial attitude transfer results corresponding to each random sampling to generate a future prediction branch. Summarize all future prediction branches generated by random sampling to expand and generate multiple travel routes. Patch all future prediction branches with the current corrected trajectory at the starting end to obtain a three-dimensional trajectory tree containing the current corrected trajectory and multiple future prediction branches. The preset number of state transition steps is 10 steps.

[0126] This implementation introduces a branch inference algorithm based on Partially Observable Markov Decision Process (POMDP) ​​as an innovative technology, which has significant differences and advantages compared to traditional deterministic trajectory prediction and standard Markov Decision Process (MDP) techniques. Traditional deterministic inference methods usually assume that formation properties are known or evolve only according to the most probable value, completely ignoring the uncertainty of the formation ahead while drilling. Once an unknown fault or high-pressure layer is encountered, it is very easy to cause the trajectory to go out of control. Although standard MDP supports decision planning, it requires the system state to be completely observable, which cannot adapt to the objective reality of huge blind spots and probability fluctuations in geological information ahead in deep drilling.

[0127] This invention utilizes the POMDP model to construct a discrete probability distribution table by combining the attribute values ​​and probability weights of uncertain spatial points. During multi-step state transitions, random sampling simulates various geological changes, enabling multi-branch probabilistic deduction driven by uncertainty. This mechanism dynamically calculates the attitude deflection angle and resultant force deflection direction in each deduction step, combining mechanical response with geological probability depth to generate a three-dimensional trajectory tree containing multiple travel paths. This design can comprehensively cover and proactively deduce various spatial attitude evolution paths that the drill bit may encounter under complex disturbances of uncertain geological properties, breaking the limitation of traditional single-line prediction that leads to loss of control upon encountering anomalies. While significantly reducing the engineering risks caused by blind drilling, it significantly enhances the adaptability of the trajectory planning system to unknown geology and the robustness of risk avoidance.

[0128] In this embodiment, the space risk assessment module specifically includes:

[0129] Spatial points whose stratigraphic attribute values ​​reach a preset fault safety threshold are extracted from a local 3D geological model and connected to generate surrounding faults. Spatial points whose stratigraphic equivalent density is greater than a preset high pressure safety threshold are extracted and connected to generate high pressure layers. Pre-stored adjacent well trajectory coordinate data are read to generate adjacent wells. The preset fault safety threshold is 50 meters and the preset high pressure safety threshold is 30 meters.

[0130] The Euclidean distance between each trajectory sampling point in the 3D trajectory tree and the surrounding fault is calculated as the fault distance; the Euclidean distance between each trajectory sampling point in the 3D trajectory tree and the high-pressure layer is calculated as the high-pressure distance; and the Euclidean distance between each trajectory sampling point in the 3D trajectory tree and the adjacent well is calculated as the adjacent well distance.

[0131] The fault distance, high pressure distance, and adjacent well distance are compared with their corresponding safety thresholds. When the fault distance is less than the preset fault safety threshold, the trajectory sampling point is marked as a fault risk point. When the high pressure distance is less than the preset high pressure safety threshold, the trajectory sampling point is marked as a high pressure risk point. When the adjacent well distance is less than the preset anti-collision safety threshold, the trajectory sampling point is marked as an anti-collision risk point. The preset anti-collision safety threshold is 10 meters.

[0132] All fault risk points, high-pressure risk points, and collision avoidance risk points are merged and collectively referred to as risk space points. The Euclidean distance between adjacent risk space points is calculated. Adjacent risk space points with an Euclidean distance less than a preset connection distance are connected to each other. A spatial grid surface is constructed based on all the connections, and the surface is filled with solids to generate a three-dimensional risk boundary. The preset connection distance is 5 meters.

[0133] In this embodiment, the through-view visualization rendering module specifically includes:

[0134] The local 3D geological model is read, and the opacity parameters of all voxels in the local 3D geological model are multiplied by a preset semi-transparency coefficient to calculate the semi-transparency opacity. The semi-transparency opacity is used to replace the original opacity for ray casting rendering, and the local 3D geological model is displayed in a semi-transparent form. The preset semi-transparency coefficient is 0.2.

[0135] Read the 3D trajectory tree and risk boundary, connect the coordinates of each trajectory sampling point in the 3D trajectory tree with preset line segments to generate solid trajectory lines, and stitch together all surface space points in the risk boundary with triangular facets to generate solid boundary surfaces. Rasterize and color the solid trajectory lines and solid boundary surfaces, and display the 3D trajectory tree and risk boundary in solid line surface form.

[0136] Calculate the shortest Euclidean distance between each trajectory line in the 3D trajectory tree and the risk boundary surface. When the shortest Euclidean distance is equal to zero, mark the corresponding trajectory line as an intersecting trajectory segment and extract the coordinates of the two intersection points of the intersecting trajectory segment entering and exiting the risk boundary.

[0137] The intersecting trajectory segment between the coordinates of two intersection points is extracted as the trajectory segment that crosses the risk boundary. The original rendering color of the trajectory segment that crosses the risk boundary is replaced with a preset highlight color. The luminous intensity of the preset highlight color is multiplied by a sine time function to calculate the dynamic flicker brightness value. The trajectory segment that crosses the risk boundary is then highlighted according to the dynamic flicker brightness value. The preset highlight color is red.

[0138] The semi-transparent local 3D geological model, the solid line surface 3D trajectory tree and risk boundary, and the trajectory segment passing through the risk boundary with dynamic flashing brightness value light highlighting are deeply buffered and blended in the same 3D spatial coordinate system to produce a geological-trajectory 3D visualization screen with a sense of penetration and visual warning of risk areas.

[0139] In this embodiment, the interactive reverse control module specifically includes:

[0140] The system receives action commands from operators to drag and select predicted branches or modify risk boundaries in the 3D trajectory tree in the geological-trajectory 3D visualization screen, extracts the screen pixel coordinates corresponding to the action commands, multiplies the screen pixel coordinates by a preset camera inverse projection matrix to calculate the normalized ray direction, multiplies the normalized ray direction by a preset depth value and adds the camera spatial coordinates to calculate the target spatial coordinates, where the preset depth value is 50 meters.

[0141] Construct a tangent vector at the target space coordinates along the extension direction dragged by the operator, calculate the Euclidean norm of the tangent vector as the modulus, divide the tangent vector by the modulus to calculate the unit direction vector, and combine the unit direction vector with the target space coordinates to obtain the target trajectory attitude.

[0142] The angle between the unit direction vector in the target trajectory attitude and the three-dimensional direction vector of the current corrected trajectory is calculated as the attitude deflection angle. The sine and cosine values ​​of the attitude deflection angle are calculated. The mechanical torque and lateral force of the drill bit are decomposed by combining the sine and cosine values ​​of the attitude deflection angle to obtain the required mechanical force in the lateral and longitudinal directions. The required mechanical force in the lateral and longitudinal directions are vector-added to calculate the magnitude and direction of the mechanical force.

[0143] The drill pressure release factor is calculated by dividing the magnitude of the applied mechanical force by the vector sum of the drill bit lateral force and the drill bit mechanical torque. The target drill pressure is calculated by multiplying the drill pressure in the current drilling data by the drill pressure release factor. The target rotational speed is calculated by multiplying the rotational speed in the current drilling data by the drill pressure release factor.

[0144] The tool face deflection angle is calculated by projecting the mechanical force direction onto a plane perpendicular to the three-dimensional direction vector of the current correction trajectory. The tool face deflection angle is then added to the tool face angle in the current drilling data to calculate the target tool face angle. The target drill pressure, target rotation speed, and target tool face angle are combined to obtain the drilling parameter instruction set of drill pressure, rotation speed, and tool face angle required to achieve the target trajectory attitude.

[0145] Example 1: To verify the feasibility of this invention in the field of geological risk identification and trajectory dynamic control during drilling in complex formations, this invention was deployed in a deep shale gas horizontal well project in the Southwest Oil and Gas Field, managed by a certain limited company. The well was designed to reach a depth of 6580 meters, with a horizontal section length of 2200 meters, targeting the deep shale gas reservoir of the Longmaxi Formation. The work area is located in the complex tectonic zone of the eastern Sichuan Basin, characterized by dense underground fault development, dramatic micro-structural changes, and the interaction of multiple high-pressure gas layers with weak interlayers. During actual drilling, traditional measurement-while-drilling (MWD) and geological steering systems faced significant challenges. Due to the limited resolution of seismic data and the transmission delay and blind spots in MWD signals, on-site engineers found it difficult to accurately and promptly perceive uncertain geological risks ahead of the drill bit. This often led to serious leakage when encountering unknown faults, or overflows when encountering high-pressure layers, and even a significant risk of wellbore collisions during horizontal section extension due to untimely updates of anti-collision data from adjacent wells. Traditional drilling early warning and trajectory control modes rely excessively on the static geological experience and post-hoc comprehensive analysis of field experts, and do not make sufficient use of multi-source heterogeneous data. The isotropic assumption of the geological model deviates from the anisotropic characteristics of the actual strata, and the trajectory adjustment lacks forward-looking probability extrapolation and intuitive human-machine collaboration mechanisms, which can easily lead to a significant extension of the drilling cycle and an extremely high proportion of non-productive time.

[0146] In the actual deployment of this project, the system of this invention comprehensively integrates real-time drilling parameter data from the integrated logging tool, gamma and electromagnetic resistivity data while drilling, 3D seismic data, and logging and trajectory data from adjacent wells, processing and fusing more than 850GB of multi-source heterogeneous geological and engineering data daily. Firstly, addressing the strong anisotropy of deep shale formations, the system abandons the traditional static overall modeling approach and adopts an improved NeRF network to establish a local dynamic coordinate system with a radius of 50 meters around the drill bit in real time. By introducing a local orthogonal coordinate system decoupling mechanism, the cross-bedding direction and bedding direction of the formation are physically decoupled and their anisotropic characteristics encoded. This allows the bedding direction to capture low-frequency, gradual changes, while the cross-bedding direction captures high-frequency, dramatic responses. This achieves high-precision local 3D geological model reconstruction that grows in real time as the drill bit advances, completely solving the problem of thin-layer and micro-fault identification distortion caused by the isotropic assumption in traditional methods. Building upon this foundation, and addressing the uncertainties in the geological attributes ahead, the system employs a branch programming algorithm based on a partially observable Markov decision process. It randomly samples the unknown space according to the probability weights of data reliability, generating a 3D trajectory tree containing multiple possible routes through forward inference. Simultaneously, the system automatically extracts surrounding faults, high-pressure zones, and adjacent well risk sources, generating a 3D risk boundary to spatially constrain the trajectory tree. This risk overlap area is then visually presented in a 3D visualization using penetrating rendering and dynamic red highlighting. When on-site operators select a safe target branch or avoid risk boundaries by dragging in the 3D view, the interactive reverse control module instantly decomposes the target attitude, automatically deriving the target drilling pressure, target rotation speed, and target tool face angle drilling parameter instruction set. This completely changes the outdated mode of manually calculating and adjusting parameters based on experience.

[0147] During a four-month field application, the method of this invention demonstrated significant advantages in terms of timeliness of geological risk early warning during drilling, accuracy of trajectory control, and reliability of information disclosure in complex formations. Table 1 below shows the core performance comparison data between the method of this invention and the traditional static geological steering method under typical drilling risk scenarios during the actual drilling of this shale gas horizontal well:

[0148] Table 1. Comparison of the overall performance of the present invention and traditional methods.

[0149] Based on the comparative data shown in Table 1, it can be seen that the drilling three-dimensional visualization simulation system based on the improved NeRF and partially observable Markov decision process proposed in this invention has significant performance advantages over the traditional static guidance method in terms of complex geological risk identification and dynamic trajectory control. In particular, it has achieved comprehensive improvement in key indicators such as early warning accuracy, response lead time, trajectory control quality, and false alarm and missed alarm control.

[0150] In terms of early warning accuracy and response timeliness, this invention maintains an accuracy rate of over 92% in all four typical drilling risk scenarios, far exceeding the average of approximately 58% for traditional methods. Furthermore, this invention achieves a significant improvement in early warning response distance through the POMDP branch extrapolation mechanism. For example, in the "adjacent well collision avoidance early warning" scenario, traditional methods rely on post-hoc logging and static collision avoidance scanning, resulting in signal delays and spatial blind spots, with an average early response distance of only 22.4 meters. In contrast, this invention, through uncertain random sampling and multi-step attitude transfer extrapolation, achieves an average early response distance of up to 65.8 meters, providing drillers with nearly three times the decision-making and adjustment space, and greatly enhancing real-time risk avoidance capabilities in complex well networks.

[0151] Regarding trajectory control and reservoir encounter quality, this invention maintains a trajectory smoothness compliance rate and a high-quality reservoir encounter rate of over 90% and 84% respectively across all scenarios, while traditional methods typically hover around 65% and 60%. Traditional methods often employ isotropic assumptions and empirical adjustments, which can easily lead to excessive dogleg severity or frequent breakouts in thin reservoirs. This invention utilizes a local orthogonal coordinate system decoupling mechanism to accurately characterize formation anisotropy, combined with interactive inverse control to directly output the target drilling parameter command set, resulting in more refined and smoother trajectory adjustments and effectively ensuring long-distance horizontal penetration within high-quality reservoirs.

[0152] The invention also demonstrates significant advantages in controlling false alarm and false negative rates. The average false alarm rate is controlled to within 5%, and the false negative rate is around 3%. Compared with the traditional method, which has a false alarm rate of nearly 20% and a false negative rate of about 16%, this invention significantly reduces redundant early warnings and blind spot omissions, and lowers the costs of ineffective drilling shutdowns and complex accident handling.

[0153] Overall, this invention achieves efficient, accurate, and visualized risk identification and trajectory control in key drilling scenarios such as "unknown fault identification," "early avoidance of high-pressure layers," "collision warning with adjacent wells," and "traversing lost circulation zones." Through multi-source probabilistic fusion, anisotropic dynamic reconstruction, and uncertainty risk extrapolation, it demonstrates broad practical value and promising prospects for deployment and promotion. The improvements in spatial early warning distance and reservoir encounter quality are particularly significant, filling the technical shortcomings of traditional methods in the forward-looking extrapolation of dynamic risks in deep wells.

[0154] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A three-dimensional visualization simulation system for drilling trajectories, characterized in that, Includes the following modules: The multi-source data fusion and probabilistic module is used to preprocess the collected drilling measurement data, seismic wave data and rock mechanics parameters to obtain a probabilistic map of geological features. The local dynamic geological reconstruction module is used to determine the influence radius, estimate the morphology and properties of unknown strata ahead, correct the strata morphology, construct an improved NeRF network, introduce a local orthogonal coordinate system decoupling mechanism, establish a local orthogonal coordinate system, calculate anisotropic coding, spatial density and strata property values, perform volume rendering reconstruction and isosurface extraction, and obtain a local three-dimensional geological model. The trajectory mechanics correction module is used to calculate the lateral deflection force generated when the drill bit breaks rock based on the current drilling data and the local three-dimensional geological model, and to correct the true spatial attitude of the current drilling trajectory to obtain the current corrected trajectory. The uncertainty branching inference module is used to address the uncertainty of the geological data ahead in the geological feature probability map. It adopts a branching programming algorithm based on a partially observable Markov decision process to perform random sampling and forward inference of the drill bit's multi-step spatial attitude transfer, expand and generate multiple travel routes, and obtain a three-dimensional trajectory tree. The spatial risk calibration module is used to calculate the spatial distance between the three-dimensional trajectory tree and the surrounding faults, high-pressure layers and adjacent wells in the local three-dimensional geological model, and to generate risk boundaries. The penetrating visualization rendering module is used to display local 3D geological models in a semi-transparent form, display 3D trajectory trees and risk boundaries in a solid line and surface form, and highlight the trajectory segments of the 3D trajectory tree that cross the risk boundary to obtain a 3D visualization of geology and trajectory. The interactive reverse control module receives the operator's action commands to drag and select the predicted branch in the 3D trajectory tree or modify the risk boundary in the geological-trajectory 3D visualization screen, and converts them into target spatial coordinates and target trajectory posture. It then reverse-derives and calculates the magnitude and direction of the mechanical force required for the drill bit to reach the target trajectory posture, and obtains the drilling parameter command set for drilling pressure, rotation speed and tool face angle.

2. The three-dimensional visualization simulation system for drilling trajectories according to claim 1, characterized in that, The multi-source data fusion and probabilistic module specifically includes: The original depth and coordinate information carried by the collected drilling measurement data, seismic wave data and rock mechanics parameters are read, and coordinate translation and rotation transformation calculations are performed to align them to the same spatial coordinate system. The attribute values ​​of the drilling measurement data, seismic wave data and rock mechanics parameters under the unified coordinate system are extracted and normalized to complete the multi-source data fusion standardization. In the fusion of standardized multi-source data, the difference between adjacent sampling points of seismic wave impedance is calculated along the depth. When it is greater than the preset impedance threshold, the corresponding depth coordinates are marked as the position of the wave impedance interface. The slope of adjacent sampling points of the logging curve of the measurement while drilling data is calculated along the depth. When it is greater than the preset slope threshold, the corresponding depth coordinates are marked as the position of the abrupt change point of the logging curve. The location of the acoustic impedance interface at the same spatial position is matched and compared with the location of the abrupt change point of the logging curve. If the positions coincide, it is determined as a boundary. If they deviate, the vertical and horizontal distances between the acoustic impedance interface and the location of the abrupt change point of the logging curve are calculated. The two-dimensional spatial range of the vertical and horizontal distances is delineated. Within the corresponding two-dimensional spatial range, the acoustic impedance interface and the location of the abrupt change point of the logging curve are connected sequentially along the vertical and horizontal directions to form a broken line segment. The broken line segment is used as the deviation trajectory, and the spatial range of the broken line segment is separated and output as the formation boundary trend. In the fusion of standardized multi-source data, the difference between the elastic modulus and Poisson's ratio in rock mechanics parameters is extracted, and the wave velocity difference gradient in seismic wave data is extracted. The difference between the elastic modulus and Poisson's ratio and the wave velocity difference gradient are multiplied and superimposed at the corresponding spatial coordinates. Based on the extreme value distribution of the multiplication and superposition calculation results, the direction of the region of abrupt stress distribution change is delineated. The direction of the region of abrupt stress distribution change is separated and the output is the difference in the direction of rock stress. The differences in the trend of the separated stratigraphic boundaries and the direction of rock stress are assigned spatial coordinates in a unified coordinate system. The signal-to-noise ratio and sampling density values ​​of each source data at the corresponding location are statistically compared. The confidence score is calculated by weighted summation of the signal-to-noise ratio and sampling density values, and then normalized to obtain the probability weight of data reliability. The spatial coordinates are bound and combined with the probability weight of data reliability to output a probability map of geological features.

3. The three-dimensional visualization simulation system for drilling trajectories according to claim 1, characterized in that, The local dynamic geological reconstruction module specifically includes: Using the current drill bit position as the origin of spatial coordinates, and the preset distance threshold as the radius of influence, the reliability probability weight of each spatial point within the radius of influence in the geological feature probability map is read. When the weight is lower than the preset weight threshold, the corresponding spatial point is marked as an unknown stratum ahead. Read the current drill bit's drilling pressure, rotation speed, drill bit radius, drill bit cross-sectional area, number of drill bit cutting teeth, and gravity of the overlying rock strata. Calculate the volumetric work done by a single tooth breaking the rock and the horizontal stress, and combine them to obtain the morphology and properties of the unknown strata ahead, and continuously correct the strata morphology. An improved NeRF network is constructed, isotropic basic feature vectors are calculated, and a local orthogonal coordinate system decoupling mechanism is introduced to establish a local orthogonal coordinate system and calculate the anisotropic encoded features after decoupling. Anisotropic coding features are input into the density prediction branch. The result of multiplying the cross-layer direction features by a preset cross-layer weight coefficient and the result of multiplying the in-layer direction features by a preset in-layer weight coefficient are added element-wise. The spatial density value is calculated by using the Softplus activation function. In the attribute prediction branch, the anisotropic coding features and the spatial density value are concatenated in the channel dimension. After mapping through two fully connected layers, the formation attribute value is output. Within the local area of ​​influence radius, samples are uniformly taken along the preset camera ray direction to obtain sampling points. The influence radius is multiplied by 2 to calculate the total ray penetration length. The total ray penetration length is divided by the preset number of sampling points to calculate the ray step size between adjacent sampling points. The spatial density value corresponding to each sampling point is multiplied by the ray step size to obtain the single-point transparency. The single-point transparency is subtracted from 1 to obtain the single-point transmittance. The single-point transmittance of all sampling points upstream of the ray is multiplied together to obtain the cumulative transmittance. The cumulative transmittance is multiplied by the single-point transparency and the stratigraphic attribute value to obtain the single-point rendering attribute value. The single-point rendering attribute values ​​of all sampling points on the ray are summed to obtain the ray cumulative rendering value. The cumulative rendering values ​​of each ray are arranged and stitched according to the corresponding spatial sampling grid coordinates to generate a local three-dimensional geological attribute matrix as a local three-dimensional geological body. The spatial density values ​​of all sampling points are compared with the preset density threshold. Sampling points with spatial density values ​​greater than or equal to the preset density threshold are retained as boundary points, while sampling points with spatial density values ​​less than the preset density threshold are set to zero. All non-zero boundary points are connected adjacently in space to generate a stratigraphic boundary surface. The local three-dimensional geological body is combined with the stratigraphic boundary surface to obtain a local three-dimensional geological model.

4. The three-dimensional visualization simulation system for drilling trajectories according to claim 3, characterized in that, The process involves reading the current drill bit's drilling pressure, rotation speed, drill bit radius, drill bit cross-sectional area, number of cutting teeth, and overlying rock gravity; calculating the volumetric work done by a single tooth breaking the rock and the horizontal stress; combining these calculations to obtain the morphology and properties of the unknown strata ahead; and continuously correcting the strata morphology, specifically including: Read the current drill bit's drilling pressure, rotation speed, drill bit radius, drill bit cross-sectional area, number of drill bit cutting teeth, and gravity of the overlying rock strata. Multiply the rotation speed by the drill bit radius to calculate the drill bit cutting linear velocity, divide the drilling pressure by the drill bit cross-sectional area to calculate the drill bit contact pressure, and multiply the cutting linear velocity by the contact pressure and the cutting tooth width to calculate the single-tooth rock-breaking volumetric work. Multiply the volumetric work of single-tooth rock breaking by the number of drill bit cutting teeth and add the preset initial cohesion of the rock to calculate the specific work of rock breaking of the unknown strata ahead. Add the specific work of rock breaking to the gravity of the overlying strata to calculate the equivalent density of the strata. Multiply the density by the gravitational acceleration and the depth of the corresponding spatial point to calculate the vertical stress. Multiply the vertical stress by the preset horizontal stress coefficient to calculate the horizontal stress. The horizontal geostress is mapped to the spatial surface undulation of the strata, and the equivalent density of the strata and the rock fracture work are mapped to the rock porosity and permeability properties of the strata. The morphology and properties of the unknown strata ahead are obtained. When the probability weight of the data reliability is not lower than the preset weight threshold, the spatial point is marked as the drilled strata behind. The measured strata morphology and property data obtained by drilling measurement are directly read to replace the calculation results and complete the strata morphology correction.

5. The three-dimensional visualization simulation system for drilling trajectories according to claim 3, characterized in that, The improved NeRF network is constructed, isotropic fundamental feature vectors are calculated, and a local orthogonal coordinate system decoupling mechanism is introduced to establish a local orthogonal coordinate system. The decoupled anisotropic encoded features are then calculated and output, specifically including: An improved NeRF network is constructed, which concatenates the spatial coordinates and spatial location coordinates within the influence radius as the query input, and maps the query input to an isotropic basic feature vector of a preset dimension through an initial fully connected layer; A local orthogonal coordinate system decoupling mechanism is introduced. The isotropic basic feature vector is input into the newly added normal vector prediction branch in the improved NeRF network. After a linear mapping layer, a 3D orientation initial vector is output. The Euclidean norm of the 3D orientation initial vector is calculated as the modulus. Each component of the 3D orientation initial vector is divided by the modulus to obtain the unit normal vector. The unit normal vector is used as the Z-axis of the local orthogonal coordinate system to characterize the stratigraphic cross-layer direction. On a plane perpendicular to the unit normal vector, define an arbitrary auxiliary vector that is not parallel to the unit normal vector. Perform a cross product operation between the unit normal vector and the arbitrary auxiliary vector and divide by the magnitude of the cross product result to obtain the first orthogonal unit vector, which serves as the X-axis of the local orthogonal coordinate system, representing the strike of the strata. Perform a cross product operation between the unit normal vector and the first orthogonal unit vector to obtain the second orthogonal unit vector, which serves as the Y-axis of the local orthogonal coordinate system, representing the dip of the strata. The local orthogonal coordinate system is established by combining the first orthogonal unit vector, the second orthogonal unit vector, and the unit normal vector. The spatial coordinates within the influence radius are multiplied by the first orthogonal unit vector, the second orthogonal unit vector, and the unit normal vector respectively to obtain the local X coordinate, local Y coordinate, and local Z coordinate in the local orthogonal coordinate system. These coordinates are then added together to form the in-layer coordinates, and the local Z coordinate is used as the trans-layer coordinate. The isotropic basic feature vector is concatenated with the in-layer and cross-layer coordinates and input into the newly added anisotropic coding layer in the improved NeRF network. In the in-layer direction, the in-layer coordinates are input into the sine and cosine functions with low-frequency preset parameters for periodic mapping. The mapping results are concatenated and multiplied by the first set of learnable weight matrices to extract the in-layer features. In the layer-penetrating direction, the layer-penetrating direction coordinates are respectively input into a sine function and a cosine function with high-frequency preset parameters for periodic mapping. The mapping results are concatenated and multiplied by the second set of learnable weight matrices to extract the layer-penetrating direction features. The in-layer direction features and the layer-penetrating direction features are added element by element to obtain the decoupled anisotropic coding features.

6. The three-dimensional visualization simulation system for drilling trajectories according to claim 1, characterized in that, The trajectory mechanics correction module specifically includes: Read the drill bit lateral force and drill bit mechanical torque from the current drilling data, extract the formation equivalent density and rock breaking specific work corresponding to the current spatial position of the drill bit from the local three-dimensional geological model, and calculate the rock hardness index by dividing the rock breaking specific work by the formation equivalent density. Along the direction of drill bit travel, extract the equivalent density of the formation at a preset distance to the left and a preset distance to the right of the drill bit in the local three-dimensional geological model. Divide the equivalent density of the formation on the left by the equivalent density of the formation on the right to calculate the density ratio of the formation on the left and right. Along the direction of drill bit travel, extract the equivalent density of the formation at a preset distance in front of the drill bit and a preset distance below the drill bit. Divide the equivalent density of the formation in front by the equivalent density of the formation below to calculate the density ratio of the formation in front and below. The lateral formation deflection force is calculated by subtracting 1 from the density ratio of the left and right formations and multiplying it by the lateral force of the drill bit. The longitudinal formation deflection force is calculated by subtracting 1 from the density ratio of the front and rear formations and multiplying it by the mechanical torque of the drill bit. The lateral rock-breaking deflection force is calculated by multiplying the rock hardness index by the lateral formation deflection force. The longitudinal rock-breaking deflection force is calculated by multiplying the rock hardness index by the longitudinal formation deflection force. The lateral rock-breaking deflection force is calculated by vector addition of the lateral rock-breaking deflection force and the longitudinal rock-breaking deflection force. The reference mechanical vector is constructed by combining the drill bit lateral force and drill bit mechanical torque. The trajectory deflection coefficient is obtained by ratio calculation of the lateral deflection force and the reference mechanical vector. The spatial attitude offset is calculated by multiplying the trajectory deflection coefficient by the three-dimensional direction vector of the current drilling trajectory. The corrected three-dimensional direction vector is calculated by subtracting the spatial attitude offset from the three-dimensional direction vector of the current drilling trajectory. The corrected three-dimensional direction vector is then combined with the coordinate position of the current drilling trajectory to obtain the current corrected trajectory.

7. The three-dimensional visualization simulation system for drilling trajectories according to claim 1, characterized in that, The uncertainty branch deduction module specifically includes: The probability weights of data reliability of each spatial point within the influence radius in front of the current drill bit are read from the geological feature probability map. Spatial points with weights below the preset threshold are extracted as uncertain spatial points. The stratigraphic attribute values ​​of the uncertain spatial points are combined with the probability weights of data reliability to construct a discrete probability distribution table. Based on the partially observable Markov decision process, the number of state transition steps is set. In each state transition step, the stratigraphic attribute values ​​of uncertain spatial points are randomly sampled according to the discrete probability distribution table. The stratigraphic attribute values ​​obtained from each random sampling are used as the simulated stratigraphic attributes of the corresponding state transition step. Read the drill pressure, rotation speed and tool face angle from the current drilling data, combine the three-dimensional direction vector of the current correction trajectory with the simulated formation properties to calculate the rock hardness index and formation deflection force, divide the formation deflection force by the product of drill pressure and rotation speed to calculate the attitude deflection angle, add the attitude deflection angle to the tool face angle to calculate the resultant force deflection direction, and add the resultant force deflection direction to the three-dimensional direction vector of the current correction trajectory to calculate the next predicted direction vector; The coordinate offset is calculated by multiplying the next predicted direction vector by the preset single-step advance length, and the next predicted coordinate position is calculated by adding the current spatial position coordinates. The next predicted direction vector and the next predicted coordinate position are combined to form the single-step spatial attitude transfer result. Repeatedly perform random sampling and spatial attitude transfer calculations until the preset number of state transfer steps is reached. Connect the single-step spatial attitude transfer results corresponding to each random sampling to generate a future prediction branch. Summarize all future prediction branches generated by random sampling to expand and generate multiple travel routes. Patch all future prediction branches with the current corrected trajectory at the starting end to obtain a three-dimensional trajectory tree.

8. The three-dimensional visualization simulation system for drilling trajectories according to claim 1, characterized in that, The space risk assessment module specifically includes: Spatial points whose stratigraphic attribute values ​​reach the preset fault safety threshold are extracted from the local 3D geological model and connected to generate surrounding faults. Spatial points whose stratigraphic equivalent density is greater than the preset high pressure safety threshold are extracted and connected to generate high pressure layers. Pre-stored adjacent well trajectory coordinate data are read to generate adjacent wells. The Euclidean distance between each trajectory sampling point in the 3D trajectory tree and the surrounding fault is calculated as the fault distance; the Euclidean distance between each trajectory sampling point in the 3D trajectory tree and the high-pressure layer is calculated as the high-pressure distance; and the Euclidean distance between each trajectory sampling point in the 3D trajectory tree and the adjacent well is calculated as the adjacent well distance. The fault distance, high pressure distance, and adjacent well distance are compared with their corresponding safety thresholds. When the fault distance is less than the preset fault safety threshold, the trajectory sampling point is marked as a fault risk point. When the high pressure distance is less than the preset high pressure safety threshold, the trajectory sampling point is marked as a high pressure risk point. When the adjacent well distance is less than the preset anti-collision safety threshold, the trajectory sampling point is marked as an anti-collision risk point. All fault risk points, high-pressure risk points, and collision avoidance risk points are merged and collectively referred to as risk space points. The Euclidean distance between adjacent risk space points is calculated. Adjacent risk space points with an Euclidean distance less than a preset connection distance are connected to each other. A spatial mesh surface is constructed based on all the connections, and the surface is filled with solids to generate a three-dimensional risk boundary.

9. The three-dimensional visualization simulation system for drilling trajectories according to claim 1, characterized in that, The penetration-based visualization rendering module specifically includes: Read the local 3D geological model, multiply the opacity parameter of all voxels in the local 3D geological model by the preset semi-transparency coefficient to calculate the semi-transparency opacity, replace the original opacity with the semi-transparency opacity for ray casting rendering, and display the local 3D geological model in a semi-transparent form. Read the 3D trajectory tree and risk boundary, connect the coordinates of each trajectory sampling point in the 3D trajectory tree with preset line segments to generate solid trajectory lines, and stitch together all surface space points in the risk boundary with triangular facets to generate solid boundary surfaces. Rasterize and color the solid trajectory lines and solid boundary surfaces, and display the 3D trajectory tree and risk boundary in solid line surface form. Calculate the shortest Euclidean distance between each trajectory line in the 3D trajectory tree and the risk boundary surface. When the shortest Euclidean distance is equal to zero, mark the corresponding trajectory line as an intersecting trajectory segment and extract the coordinates of the two intersection points of the intersecting trajectory segment entering and exiting the risk boundary. Extract the intersecting trajectory segment between the coordinates of two intersection points as the trajectory segment that crosses the risk boundary, replace the original rendering color of the trajectory segment that crosses the risk boundary with the preset highlight color, multiply the luminous intensity of the preset highlight color by the sine time function to calculate the dynamic flicker brightness value, and perform light highlighting processing on the trajectory segment that crosses the risk boundary according to the dynamic flicker brightness value. The semi-transparent local 3D geological model, the solid line surface 3D trajectory tree and risk boundary, and the trajectory segment passing through the risk boundary with dynamic flashing brightness value light highlighting are deeply buffered and fused in the same 3D spatial coordinate system to obtain a geological-trajectory 3D visualization.

10. The three-dimensional visualization simulation system for drilling trajectories according to claim 1, characterized in that, The interactive reverse control module specifically includes: The system receives action commands from operators to drag and select predicted branches or modify risk boundaries in the 3D trajectory tree in the geological-trajectory 3D visualization screen, extracts the screen pixel coordinates corresponding to the action commands, multiplies the screen pixel coordinates by the preset camera inverse projection matrix to calculate the normalized ray direction, multiplies the normalized ray direction by the preset depth value and adds the camera spatial coordinates to calculate the target spatial coordinates. Construct a tangent vector at the target space coordinates along the extension direction dragged by the operator, calculate the Euclidean norm of the tangent vector as the modulus, divide the tangent vector by the modulus to calculate the unit direction vector, and combine the unit direction vector with the target space coordinates to obtain the target trajectory attitude. The angle between the unit direction vector in the target trajectory attitude and the three-dimensional direction vector of the current corrected trajectory is calculated as the attitude deflection angle. The sine and cosine values ​​of the attitude deflection angle are calculated. The mechanical torque and lateral force of the drill bit are decomposed by combining the sine and cosine values ​​of the attitude deflection angle to obtain the required mechanical force in the lateral and longitudinal directions. The required mechanical force in the lateral and longitudinal directions are vector-added to calculate the magnitude and direction of the mechanical force. The drill pressure release factor is calculated by dividing the magnitude of the applied mechanical force by the vector sum of the drill bit lateral force and the drill bit mechanical torque. The target drill pressure is calculated by multiplying the drill pressure in the current drilling data by the drill pressure release factor. The target rotational speed is calculated by multiplying the rotational speed in the current drilling data by the drill pressure release factor. The tool face deflection angle is calculated by projecting the mechanical force direction onto a plane perpendicular to the three-dimensional direction vector of the current correction trajectory. The tool face deflection angle is then added to the tool face angle in the current drilling data to calculate the target tool face angle. The drilling parameter command set of drilling pressure, rotational speed and tool face angle is obtained by combining the target drilling pressure, target rotational speed and target tool face angle.

Images