Horizontal hole drilling-while-drilling perception and calibration method and system for tunnel geological exploration

Abstract

The application discloses a horizontal hole drilling-while-sensing and calibration method and system for tunnel geological exploration and belongs to the technical field of tunnel exploration, and specifically comprises the following steps: arranging a vibration sensor at a working face; actively modulating a drill bit thrust according to a complementary binary coding sequence during drilling, recording a thrust modulation time and a vibration signal; pulse compressing the vibration signal, extracting an equivalent seismic source wavelet of each modulation period and giving a virtual seismic source point; performing reverse time migration imaging on the equivalent seismic source wavelet as a seismic source function and the vibration signal as an observed wave field to generate a reflection coefficient distribution; identifying an interface line segment from the reflection coefficient distribution and calculating a spatial migration amount of a lithological interface; generating a directional correction angle and stacking the directional correction angle to a design trajectory azimuth angle and an inclination angle according to the migration amount; and outputting a calibrated drilling trajectory. The application can sense the position of a lithological interface in advance during drilling, realizes forward-looking imaging and active calibration of a trajectory.

Description

Technical Field

[0001] This invention relates to the field of tunnel exploration technology, specifically to a method and system for horizontal borehole sensing and calibration for tunnel geological exploration. Background Technology

[0002] In tunnel construction, horizontal drilling is a crucial technique for obtaining geological information ahead of the tunnel face. By drilling along a designed trajectory into the unexcavated rock mass from the tunnel face, the location and scale of adverse geological features such as lithological interfaces, water-bearing fractured zones, and karst cavities can be determined in advance, providing a basis for adjusting the tunnel excavation plan and optimizing support parameters. During horizontal drilling, the vibration waves generated by the drill bit breaking through the rock propagate backward into the rock mass. These vibration waves carry structural information about the rock mass traversed by the drill bit. By deploying sensors backward to receive and analyze the vibration signals, the geological characteristics of the rock mass surrounding the drilling path can be inverted.

[0003] Existing sensing-while-drilling methods treat drill bit rock-breaking vibrations as passive seismic sources. Vibration sensors are deployed at the tunnel face and on the excavated tunnel walls to continuously record vibration signals during drilling. The velocity structure of the rock mass near the drill bit is then retrieved using the travel time of these vibration signals. This type of method is based on the natural occurrence characteristics of drill bit rock-breaking vibrations. It extracts the initial arrival travel time of the P-wave from the vibration signals, establishes the ray path from the source to the sensor, and obtains the velocity distribution through travel time tomography.

[0004] However, if the rock-breaking vibration of the drill bit is considered a passive source, its intensity is modulated by the inherent fluctuation of the drill bit's rotation period. This inherent fluctuation completely overlaps with the weak reflected signal generated by the geological interface in front of the drill bit in the frequency domain. The two signals lack distinguishable time-frequency characteristics in the received waveform, and the reflected signal component carrying the geological interface location information cannot be separately extracted from the continuously recorded vibration signal. Under these conditions, existing methods can only use the first arrival travel time of the strongest direct P-wave in the vibration signal for tomographic inversion. The ray path of the tomographic inversion is limited to the backward propagation direction from the source to the sensor. The inversion result only reflects the velocity structure distribution of the area already traversed by the drill bit. It does not have the source function and reflected wavefield input conditions required for wavefield extrapolation of the area in front of the drill bit, and cannot generate the reflection coefficient distribution of the geological interface in front of the tunnel face through imaging methods such as reverse time migration. During drilling, the spatial position of adverse geological bodies such as lithological interfaces and water-bearing fractured zones is imperceptible. When the drill bit deviates from the designed trajectory or encounters adverse geological bodies, there is a lack of forward-looking information to support active trajectory calibration. Summary of the Invention

[0005] The purpose of this invention is to provide a method and system for horizontal borehole drilling sensing and calibration for tunnel geological exploration, and to solve the following technical problems: Existing methods suffer from the inherent vibrations caused by drill bit rotation and the reflected signals from geological interfaces overlapping in the frequency domain. This makes it impossible to separate the reflected signal components from the continuous vibration signals, resulting in a lack of imaging capability for the area in front of the drill bit. Furthermore, the methods cannot detect the location of lithological interfaces in advance during drilling and lack forward-looking information to support active trajectory calibration.

[0006] The objective of this invention can be achieved through the following technical solutions: A method for horizontal borehole sensing and calibration for tunnel geological exploration includes the following steps: S1. Before drilling the horizontal hole of the tunnel, multiple vibration sensors are laid out along the tunnel face and the surface of the adjacent excavated tunnel wall, and the multiple vibration sensors are connected to the data acquisition device. S2. Start the drilling rig and control the drill rod propulsion device to apply periodic thrust changes to the drill bit according to the complementary binary coding sequence. Simultaneously record the thrust modulation time sequence, real-time position parameters of the drill bit and drilling time, and continuously record the vibration signals received by multiple vibration sensors. S3. Use the thrust modulation time sequence to perform pulse compression processing on the vibration signal, extract the drill bit equivalent source wavelet corresponding to each thrust modulation cycle, and assign the drill bit equivalent source wavelet to the virtual source point corresponding to the real-time position parameters of the drill bit. S4. Using the drill bit equivalent source wavelet as the source function and the vibration signal as the observation wave field, reverse time migration imaging is performed to generate the reflection coefficient distribution in the area in front of the tunnel face. S5. Identify interface segments from the reflection coefficient distribution, extract the spatial position of the interface segments, and compare them with the interface positions corresponding to the azimuth and dip angles of the designed drilling trajectory to obtain the spatial offset of the lithological interface. S6. Generate the drilling direction correction angle based on the spatial offset of the lithological interface, superimpose the drilling direction correction angle onto the azimuth and dip angle of the designed drilling trajectory, and output the calibrated drilling trajectory coordinate sequence.

[0007] As a further aspect of the present invention: In step S1, the process of deploying multiple vibration sensors along the tunnel face and adjacent excavated tunnel wall surfaces, and connecting the multiple vibration sensors to the data acquisition unit, is as follows: The multiple vibration sensors are distributed acoustic wave sensing optical cables. The distributed acoustic wave sensing optical cables are pasted along the tunnel face to form a spiral shape. The outermost circle of the spiral is attached to the excavation outline of the tunnel face. The distributed acoustic wave sensing optical cables extend axially along the excavated tunnel wall and pass through the arch, the two sides of the arch waist and the two sides of the sidewall. The data acquisition unit is a distributed acoustic wave sensor demodulator. The distributed acoustic wave sensor demodulator emits probe light pulses to the distributed acoustic wave sensor optical cable and receives backscattered Rayleigh light signals. It demodulates the phase difference changes between adjacent spatial resolution units to obtain the vibration signals at each sensor channel position along the distributed acoustic wave sensor optical cable. The spatial coordinates of each sensor channel position of the distributed acoustic wave sensor optical cable are obtained through three-dimensional measurement, and a mapping table is established between the spatial coordinates and the sensor channel labels.

[0008] As a further aspect of the present invention: In S2, the process of controlling the drill pipe propulsion device to apply periodic thrust changes to the drill bit according to a complementary binary coding sequence is as follows: An electronically controlled proportional pressure reducing valve is installed on the hydraulic supply circuit of the drill pipe propulsion device. The electronically controlled proportional pressure reducing valve receives the first and second encoded sequences output by the encoder controller. The sum of the autocorrelation of the first and second encoded sequences is the impulse function. The side lobes of the autocorrelation function of the first and second encoded sequences have the same amplitude and opposite signs as the side lobes of the autocorrelation function of the second encoded sequence. The drill pipe propulsion device outputs thrust modulation according to the first encoding sequence in the first time period and outputs thrust modulation according to the second encoding sequence in the second time period. The high level in the encoding sequence corresponds to the thrust setting base value, and the low level corresponds to the thrust setting base value plus the preset thrust increase value. The thrust modulation period interval is an integer multiple of the width of the encoding sequence chip, and the width of the encoding sequence chip is set to a non-integer multiple of the drill bit rotation period.

[0009] As a further aspect of the present invention: In step S3, the process of using the thrust modulation time sequence to perform pulse compression processing on the vibration signal and extract the equivalent seismic source wavelet of the drill bit is as follows: The vibration signal is separated into a first vibration signal segment and a second vibration signal segment according to the encoding sequence. The first vibration signal segment corresponds to a first time period, and the second vibration signal segment corresponds to a second time period. The first vibration signal segment is cross-correlated with the first coded sequence copy to obtain the first correlation signal, and the second vibration signal segment is cross-correlated with the second coded sequence copy to obtain the second correlation signal. The first and second correlation signals are aligned along the time axis and then added point by point to obtain the pulse compression signal after sidelobe suppression. A bandpass filter is applied to the pulse compression signal to filter out the periodic vibration component caused by the drill bit rotation. The maximum energy phase axis is tracked along the thrust modulation time sequence in the bandpass filtered pulse compression signal, and the waveform segment at the time delay corresponding to the maximum energy phase axis is extracted as the drill bit equivalent source wavelet. The drill bit equivalent source wavelet is assigned to the corresponding virtual source point spatial coordinates on the drilling path according to the drill bit real-time position parameters.

[0010] As a further aspect of the present invention: In step S4, the process of performing reverse-time migration imaging processing using the drill bit equivalent source wavelet as the source function and the vibration signal as the observed wavefield is as follows: The initial velocity model is a two-dimensional longitudinal profile velocity model that extends from the tunnel face in the drilling direction. The longitudinal wave propagation velocity values ​​of each velocity grid node in the initial velocity model are extrapolated from the wave velocity test values ​​of the exposed rock mass at the tunnel face along the drilling direction. Forward extrapolation wavefield calculation uses the location of the virtual source point as the source grid node. The equivalent source wavelet of the drill bit is loaded onto the right-hand side of the acoustic equation. The wavefield value of each velocity grid node is calculated along the forward time direction using an interleaved grid finite difference scheme. For reverse extrapolation wavefield calculation, the deployment location of multiple vibration sensors is used as the receiving grid node. The vibration signal is loaded onto the right-hand side of the acoustic equation along the reverse time direction. The wavefield value of each velocity grid node is calculated along the reverse time direction using an interleaved grid finite difference scheme. At each velocity grid node, the wavefield values ​​of the forward extrapolation wavefield and the reverse extrapolation wavefield at the same time are multiplied and summed over all time steps. The summation result is used as the reflection coefficient value of the velocity grid node. The reflection coefficient values ​​of all velocity grid nodes constitute the reflection coefficient distribution.

[0011] As a further aspect of the present invention: In step S4, after the reflection coefficient values ​​of all velocity grid nodes constitute a reflection coefficient distribution, the process further includes using the reflection coefficient distribution to correct the longitudinal wave propagation velocity values ​​of each velocity grid node in the initial velocity model. Velocity grid nodes with reflection coefficient values ​​lower than the set background level are marked as background nodes, and the P-wave propagation velocity values ​​at the background nodes are maintained as the P-wave propagation velocity values ​​of the initial velocity model. Velocity grid nodes with reflection coefficient values ​​exceeding the set background level are marked as interface influence nodes. One-dimensional wavenumber domain inversion is applied to the interface influence nodes along the drilling direction to solve for the P-wave propagation velocity correction of the interface influence nodes. The corrected initial velocity model replaces the original initial velocity model. The subsequent forward extrapolation wavefield calculation and reverse time extrapolation wavefield calculation corresponding to the virtual source point adopt the corrected initial velocity model.

[0012] As a further aspect of the present invention: in step S5, the process of obtaining the spatial offset of the lithological interface is as follows: Extract the gradient magnitude distribution of the reflection coefficient values ​​from the reflection coefficient distribution, mark the positions where the gradient magnitude exceeds the set gradient level as interface candidate points, apply Hough transform to the interface candidate points, search for parameter points in the tilt-intercept parameter space where the cumulative value exceeds the set cumulative level, and transform the searched parameter points back to the spatial domain to obtain interface line segments in the form of straight line segments. The coordinates of the intersection of the interface line segment and the drilling axis are extracted as the spatial position representation value of the interface line segment. The interface position corresponding to the azimuth and dip angle of the designed drilling trajectory is compared with the spatial position representation value of the interface line segment. The distance from the interface position corresponding to the azimuth and dip angle of the designed drilling trajectory to the spatial position representation value of the interface line segment along the drilling axis is taken as the spatial offset of the lithological interface.

[0013] As a further aspect of the present invention: the process of extracting the gradient magnitude distribution of the reflection coefficient values ​​from the reflection coefficient distribution and marking the positions where the gradient magnitude exceeds a set gradient level as interface candidate points is as follows: For each velocity grid node in the reflection coefficient distribution, the horizontal gradient component along the drilling direction and the vertical gradient component perpendicular to the drilling direction are calculated using central difference. The gradient magnitude is obtained by taking the square root of the sum of the squares of the horizontal and vertical gradient components. The gradient magnitude of each velocity grid node is compared with the gradient magnitude of the adjacent velocity grid nodes. When the gradient magnitude of a velocity grid node is greater than that of the adjacent velocity grid nodes, the velocity grid node is retained as a local extremum point. The velocity grid nodes with gradient magnitudes exceeding the set gradient level among the local extremum points are marked as interface candidate points, thus obtaining a set of interface candidate points.

[0014] As a further aspect of the present invention: In step S6, the process of generating a drilling direction correction angle based on the spatial offset of the lithological interface, and then superimposing the drilling direction correction angle onto the azimuth and dip angle of the designed drilling trajectory, is as follows: The spatial offset of the lithological interface is decomposed into a lateral offset component perpendicular to the drilling axis and a longitudinal offset component along the drilling axis. The ratio of the lateral offset component to the distance from the current real-time position parameter of the drill bit to the spatial position characterization value of the lithological interface is taken as the tangent value of the drilling direction correction angle. The drilling direction correction angle is decomposed into an azimuth correction component and an inclination correction component. The azimuth correction component is superimposed on the azimuth of the designed drilling trajectory, and the inclination correction component is superimposed on the inclination of the designed drilling trajectory. When the real-time position parameters of the drill bit indicate that the drill bit has entered the preset distance range in front of the spatial position characterization value of the lithology interface, the drilling direction correction angle decreases linearly with the remaining distance. When the drill bit reaches the spatial position characterization value of the lithology interface, the drilling direction correction angle decreases to zero.

[0015] The present invention also includes a horizontal borehole drilling-while-drilling sensing and calibration system for tunnel geological exploration, used to implement the above-described horizontal borehole drilling-while-drilling sensing and calibration method for tunnel geological exploration, comprising: The sensor deployment module is used to deploy multiple vibration sensors along the tunnel face and the adjacent excavated tunnel wall surface before drilling the horizontal hole of the tunnel, and to connect the multiple vibration sensors to the data acquisition unit. The encoding modulation acquisition module is used to start the drilling rig, control the drill rod propulsion device to apply periodic thrust changes to the drill bit according to the complementary binary encoding sequence, synchronously record the thrust modulation time sequence, the real-time position parameters of the drill bit and the drilling time, and continuously record the vibration signals received by multiple vibration sensors. The pulse compression and extraction module is used to perform pulse compression processing on the vibration signal using the thrust modulation time sequence, extract the drill bit equivalent source wavelet corresponding to each thrust modulation cycle, and assign the drill bit equivalent source wavelet to the virtual source point corresponding to the real-time position parameters of the drill bit. The reverse time migration imaging module is used to perform reverse time migration imaging processing with the drill bit equivalent source wavelet as the source function and the vibration signal as the observation wave field to generate the reflection coefficient distribution of the area in front of the tunnel face. The interface recognition and extraction module is used to identify interface segments from the reflection coefficient distribution, extract the spatial position of the interface segments, and compare them with the interface positions corresponding to the azimuth and dip angles of the designed drilling trajectory to obtain the spatial offset of the lithological interface. The trajectory calibration output module is used to generate a drilling direction correction angle based on the spatial offset of the lithological interface, superimpose the drilling direction correction angle onto the azimuth and dip angle of the designed drilling trajectory, and output the calibrated drilling trajectory coordinate sequence.

[0016] The beneficial effects of this invention are: This invention applies periodic thrust variations to the drill bit using a drill pipe propulsion device according to a complementary binary coding sequence, independent of the drill bit's rotation period. This actively modulates the intensity of the drill bit's rock-breaking vibrations, altering the time-frequency characteristics of the vibration signal and allowing the passive rock-breaking signal to be separated from the geological interface reflection signal during cross-correlation processing. In the complementary binary coding sequence, the sum of the autocorrelation of the first and second coding sequences equals the impulse function, with equal sidelobe amplitudes and opposite signs. After time-division cross-correlation, the sidelobes cancel each other out, resulting in a pulse compression signal. The drill bit equivalent source wavelet extracted from the pulse compression signal constitutes a virtual source sequence with known source functions and excitation times, providing deterministic source wavefield input conditions for reverse-time migration. Reverse-time migration employs source-normalized cross-correlation imaging to extract the reflection coefficient distribution from the background scattering by the heterogeneous rock mass, enabling real-time imaging of the geological interface in front of the drill bit during drilling. The interface line segments are identified from the reflection coefficient distribution and their spatial positions are extracted. The spatial offset of the lithological interface is obtained by comparing it with the designed drilling trajectory. Based on this, the drilling direction correction angle is generated to perform feedforward calibration of the trajectory, and the trajectory adjustment is completed before the drill bit contacts the poor geological body. Attached Figure Description

[0017] The invention will now be further described with reference to the accompanying drawings.

[0018] Figure 1This is a flowchart illustrating the horizontal borehole sensing and calibration method for tunnel geological exploration according to the present invention. Figure 2 This is a schematic diagram of the module of the horizontal hole drilling sensing and calibration system for tunnel geological exploration according to the present invention. Detailed Implementation

[0019] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0020] Please see Figure 1 As shown, this invention is a horizontal borehole sensing and calibration method for tunnel geological exploration, comprising the following steps: S1. Before drilling the horizontal hole of the tunnel, multiple vibration sensors are laid out along the tunnel face and the surface of the adjacent excavated tunnel wall, and the multiple vibration sensors are connected to the data acquisition device. S2. Start the drilling rig and control the drill rod propulsion device to apply periodic thrust changes to the drill bit according to the complementary binary coding sequence. Simultaneously record the thrust modulation time sequence, real-time position parameters of the drill bit and drilling time, and continuously record the vibration signals received by multiple vibration sensors. S3. Use the thrust modulation time sequence to perform pulse compression processing on the vibration signal, extract the drill bit equivalent source wavelet corresponding to each thrust modulation cycle, and assign the drill bit equivalent source wavelet to the virtual source point corresponding to the real-time position parameters of the drill bit. S4. Using the drill bit equivalent source wavelet as the source function and the vibration signal as the observation wave field, reverse time migration imaging is performed to generate the reflection coefficient distribution in the area in front of the tunnel face. S5. Identify interface segments from the reflection coefficient distribution, extract the spatial position of the interface segments, and compare them with the interface positions corresponding to the azimuth and dip angles of the designed drilling trajectory to obtain the spatial offset of the lithological interface. S6. Generate the drilling direction correction angle based on the spatial offset of the lithological interface, superimpose the drilling direction correction angle onto the azimuth and dip angle of the designed drilling trajectory, and output the calibrated drilling trajectory coordinate sequence.

[0021] In a preferred embodiment of the present invention, in step S1, the process of deploying multiple vibration sensors along the tunnel face and the adjacent excavated tunnel wall surface, and connecting the multiple vibration sensors to the data acquisition device, is as follows: The multi-channel vibration sensor utilizes a distributed acoustic wave sensing optical cable. This cable uses single-mode optical fiber as the sensing medium and operates by sensing strain or vibration along the cable's path through the phase change of the backscattered Rayleigh light signal within the fiber. The cable contains one or more single-mode optical fibers, which are externally encased in metal armor and a polyethylene sheath. The cable's outer diameter ranges from 3 mm to 5 mm.

[0022] At the tunnel face, the distributed acoustic sensing optical cable is designed with a spiral path. The spiral extends outwards from the center point of the face, with two spirals in total. The radial spacing of the spirals is 20 to 30 centimeters, and the outermost spiral is attached to the inner edge of the excavation outline of the tunnel face. During attachment, a total station is used to mark the path points of the spiral on the tunnel face surface, with these points spaced 15 to 25 centimeters apart along the spiral arc. Then, epoxy resin adhesive is used to press the optical cable segment by segment onto the tunnel face rock wall. After each segment is attached, metal clips are used for auxiliary fixation at the path points. These metal clips are anchored into the rock wall with expansion bolts to a depth of 5 to 8 centimeters. After the optical cable is attached at the tunnel face, the end of the cable is led along the edge of the tunnel face to the excavated tunnel wall.

[0023] On the excavated tunnel wall, the path of the distributed acoustic sensing optical cable extends axially along the tunnel wall. After emerging from the tunnel face, the cable is first pasted along the centerline of the arch towards the tunnel entrance, with a pasting length of 15 to 25 meters. It then passes sequentially through the left arch waist, left sidewall, right sidewall, and right arch waist, finally returning to the arch to form a continuous sensing loop. On the excavated tunnel wall, the cable is pasted in the same way as at the tunnel face, using epoxy resin adhesive and securing it with metal clips. When the cable passes through the arch, arch waist, and sidewall, path points are marked at 20 to 30 centimeters intervals along the cable's axial direction. Optical cable splice boxes are installed at the start and end points of each sensing segment, with 1 to 2 meters of redundant cable length reserved inside the splice boxes.

[0024] The data acquisition unit uses a distributed acoustic wave sensor demodulator, specifically a coherent optical time-domain reflectometer with a pulse repetition frequency of 1 kHz to 10 kHz as an example. The distributed acoustic wave sensor demodulator internally includes a narrow-linewidth laser, an acousto-optic modulator, an erbium-doped fiber amplifier, and a photodetector. The laser outputs a continuous wavelength 1550 nm laser, and the acousto-optic modulator chops the continuous laser into probe light pulses. The pulse width of the probe light pulses is 10 nanoseconds to 50 nanoseconds, corresponding to a spatial resolution unit length of 1 meter to 5 meters. The spatial resolution unit length determines the minimum resolvable vibration channel spacing along the optical cable, and this spatial resolution unit length is set to be smaller than the spacing between adjacent virtual seismic source points on the drilling axis. The spacing between adjacent virtual seismic source points on the drilling axis is jointly determined by the thrust modulation cycle interval and the drilling rate of the drill pipe propulsion device; this spacing is typically 0.5 meters to 2 meters.

[0025] The distributed acoustic wave sensor demodulator is connected to the optical cable junction box via a single-mode fiber optic patch cord, emitting probe light pulses into the optical cable. As the probe light pulses propagate along the fiber, the inhomogeneity of the internal medium generates backscattered Rayleigh light signals. These backscattered signals propagate in the opposite direction along the fiber to the demodulator, where they are received by a photodetector and converted into electrical signals. The demodulator performs analog-to-digital conversion on the electrical signals and extracts the phase difference between adjacent spatial resolution units. The physical meaning of the phase difference is that when a spatial resolution unit of the optical fiber is subjected to vibration, the length of that fiber segment changes slightly, causing a change in the optical path length of the light signal within that segment, which in turn causes a change in the phase difference between the Rayleigh scattered light signals of adjacent spatial resolution units. The demodulator continuously records the time series of phase difference changes corresponding to each spatial resolution unit, using the pulse repetition frequency as the sampling rate. This time series of phase difference changes represents the vibration signal along each sensing path position of the optical cable.

[0026] The 3D measurement of the fiber optic cable bonding path was conducted after the cable bonding was completed. A total station was used, set up on a stable ground point behind the excavated tunnel wall, and a coordinate system was established through centering and leveling. Surveyors aimed at the metal clip positions along the cable bonding path as feature points and recorded the 3D coordinates of each feature point. The spatial coordinates of the sensor channel positions between adjacent metal clips were calculated using linear interpolation, with an interpolation step size of 10 cm. The spatial coordinates of each sensor channel position were mapped one-to-one with the sensor channel labels of the distributed acoustic wave sensor demodulator, establishing a mapping table between spatial coordinates and sensor channel labels, which was stored in the demodulator's storage unit.

[0027] In another preferred embodiment of the present invention, in S2, the process of controlling the drill pipe propulsion device to apply periodic thrust changes to the drill bit according to a complementary binary coding sequence is as follows: A complementary binary coded sequence consists of a first coded sequence and a second coded sequence. Both the first and second coded sequences are binary sequences composed of elements 1 and -1, and have equal lengths. Complementarity is defined as follows: when the autocorrelation function of the first coded sequence is added point-by-point to the autocorrelation function of the second coded sequence, the result has zero delays except at the zero-delay position (which is twice the sequence length), meaning the sum of the autocorrelation functions equals the impulse function. Simultaneously, the sidelobes of the autocorrelation function of the first coded sequence and the sidelobes of the autocorrelation function of the second coded sequence have equal amplitudes and opposite signs. For example, with a sequence length of 4, the first coded sequence takes the values ​​1,1,1,--1, and the second coded sequence takes the values ​​1,1,--1,1. The autocorrelation result of the first coded sequence is 4,--1,0,1, and the autocorrelation result of the second coded sequence is 4,1,0,--1. Adding them together yields 8,0,0,0. The sequence length is selected based on the drilling rate and imaging resolution requirements, and can be 16, 32, or 64.

[0028] An electrically controlled proportional pressure reducing valve is installed on the hydraulic supply circuit of the drill pipe propulsion device. The drill pipe propulsion device is a hydraulic cylinder-driven propulsion mechanism. The propulsion force of the hydraulic cylinder is supplied by high-pressure hydraulic oil provided by the hydraulic pump station, which is then regulated by the proportional pressure reducing valve. The electrically controlled proportional pressure reducing valve is a valve that controls the output hydraulic pressure through an input voltage signal, with a linear relationship between the input voltage and the output pressure. The electrically controlled proportional pressure reducing valve receives a level sequence output by an encoding controller, which is a programmable logic controller or an embedded microprocessor, internally storing a binary code table containing the first and second encoding sequences.

[0029] During drilling, the encoder controller outputs a first and a second encoded sequence in time intervals. In the first time interval, the encoder controller cyclically outputs the first encoded sequence; in the second time interval, it cyclically outputs the second encoded sequence. The duration of the first and second time intervals is equal, determined based on the horizontal hole design depth and sequence length, and can be selected from 30 to 120 seconds. In the encoder controller's output level sequence, a high level corresponds to the thrust setting base value, and a low level corresponds to the thrust setting base value plus a preset thrust increment. The thrust setting base value is the thrust value of the drill pipe propulsion device under normal drilling conditions, and the preset thrust increment is 10% to 20% of the thrust setting base value. For example, if the thrust setting base value is 50 kN and the preset thrust increment is 7.5 kN, then the input voltage of the electronically controlled proportional pressure reducing valve corresponds to 50 kN when the output is high, and 57.5 kN when the output is low.

[0030] The duration of adjacent elements in the coded sequence is defined as the chip width of the coded sequence. The thrust modulation period interval is an integer multiple of the chip width of the coded sequence multiplied by the sequence length. The chip width of the coded sequence is set to a non-integer multiple of the drill bit rotation period; for example, if the drill bit rotation period is 0.3 seconds, the chip width of the coded sequence is set to 0.4 seconds or 0.7 seconds. This setting ensures that the frequency components of the thrust modulation and the periodic vibration components caused by the drill bit's rotation do not overlap in the frequency domain, providing a frequency domain distinction basis for subsequent signal separation.

[0031] While recording vibration signals, the data acquisition unit simultaneously records the thrust modulation time sequence of the drill pipe propulsion device, the real-time position parameters of the drill bit, and the drilling time. Each recording point in the thrust modulation time sequence includes a timestamp and an coded sequence identifier. The coded sequence identifier is a binary flag bit; a flag bit of 0 indicates that the current thrust modulation time belongs to the first time period, and a flag bit of 1 indicates that the current thrust modulation time belongs to the second time period. The real-time position parameters of the drill bit are output by the displacement sensor built into the drill pipe propulsion device. The displacement sensor is a wire encoder or a magnetostrictive displacement sensor with a measurement accuracy of 1 mm. The drilling time is provided by the internal clock of the data acquisition unit, and the sampling interval is synchronized with the pulse repetition period of the distributed acoustic wave sensor demodulator.

[0032] In another preferred embodiment of the present invention, the process of using the thrust modulation time sequence to perform pulse compression processing on the vibration signal and extract the equivalent seismic source wavelet of the drill bit in step S3 is as follows: The input data for pulse compression processing consists of two parts: multi-channel vibration signals recorded by the distributed acoustic wave sensor demodulator, and thrust modulation time sequence and coded sequence identifiers recorded by the drill pipe propulsion device. Pulse compression processing is performed on a data processing computer, which is connected to the distributed acoustic wave sensor demodulator and data acquisition unit via a network cable to receive the data stream in real time.

[0033] First, the vibration signal is separated into a first vibration signal segment and a second vibration signal segment based on the encoded sequence identifier. The first vibration signal segment corresponds to all vibration signals recorded within a first time period, and the second vibration signal segment corresponds to all vibration signals recorded within a second time period. The separation process is as follows: the encoded sequence identifier for each thrust modulation moment in the thrust modulation time sequence is read. When the encoded sequence identifier indicates the first time period, the vibration signal within the corresponding time window is extracted and stored in the first vibration signal segment buffer; when the encoded sequence identifier indicates the second time period, the vibration signal within the corresponding time window is extracted and stored in the second vibration signal segment buffer. The length of the time window is equal to the duration of the thrust modulation period.

[0034] Subsequently, the first vibration signal segment is cross-correlated with a copy of the first coded sequence. The copy of the first coded sequence is a complete binary sequence of the first coded sequence stored in the memory of the data processing computer, where elements 1 are kept as 1 and elements -1 are kept as -1. The cross-correlation operation is performed on each sampling point in the first vibration signal segment: starting from the current sampling point, a signal segment of equal length to the copy of the first coded sequence is extracted; the values ​​of each sampling point in the signal segment are multiplied by the corresponding element in the copy of the first coded sequence and summed; the sum is used as the value of the first correlated signal at the current sampling point. The cross-correlation calculation for all first vibration signal segments is completed by sliding the sample points one by one, thus obtaining the first correlated signal.

[0035] Similarly, the second vibration signal segment is cross-correlated with the copy of the second coded sequence to obtain the second correlation signal. The copy of the second coded sequence is a complete binary sequence of the second coded sequence stored in the memory of the data processing computer.

[0036] The first and second related signals are aligned along the time axis and then added point by point. The time axis alignment is based on the fact that the origin of the time axis for both the first and second related signals is the start time of their respective vibration signal segments. Since the first and second vibration signal segments have a sequential order on the original time axis, the start time of the second vibration signal segment is later than the start time of the first vibration signal segment by a first time interval. During the alignment operation, the time axis of the second related signal is shifted forward by the entire first time interval, so that the first and second related signals correspond to the same absolute time on the time axis.

[0037] After point-by-point addition, the side lobes of the autocorrelation of the first coding sequence and the side lobes of the autocorrelation of the second coding sequence cancel each other out because they have equal amplitudes and opposite signs. The pulse amplitude at the zero-delay position is enhanced by superposition because the two have the same sign, thus obtaining the pulse compression signal after sidelobe suppression.

[0038] A bandpass filter is applied to the pulse compression signal. The bandpass filter is implemented using an infinite impulse response digital filter. The lower cutoff frequency of the filter passband is set to 0.5 times the thrust modulation frequency, and the upper cutoff frequency is set to 8 to 12 times the drill bit rotation frequency. The thrust modulation frequency is the reciprocal of the code sequence chip width multiplied by the sequence length. The drill bit rotation frequency is the drill bit speed divided by 60, where the speed is in revolutions per minute (rpm). For example, if the thrust modulation frequency is 0.5 Hz and the drill bit rotation frequency is 3.3 Hz, then the passband of the bandpass filter is set to 0.25 Hz to 33 Hz. The signal components filtered out by the bandpass filter include hydraulic system pressure fluctuations below 0.25 Hz, electronic noise above 33 Hz, and periodic vibration components near the drill bit rotation frequency of 3.3 Hz.

[0039] In the bandpass-filtered pulse compression signal, the maximum energy in-phase axis of the pulse compression amplitude is traced sequentially along the thrust modulation time sequence. The tracing process is as follows: on the time axis of the pulse compression signal, with each thrust modulation moment as the center, a search window is extended forward and backward by one times the width of the coded sequence chip. Within the search window, the local maximum value of the pulse compression amplitude is searched, and the time delay corresponding to the local maximum value is taken as the position of the maximum energy in-phase axis. The physical meaning of the maximum energy in-phase axis is the arrival time of the equivalent source wavelet formed by the vibration wave generated by the drill bit breaking the rock during that thrust modulation period after pulse compression.

[0040] The waveform segment corresponding to the time delay at the maximum energy in-phase axis is extracted as the drill bit equivalent source wavelet for the corresponding thrust modulation period. The extraction duration of the waveform segment is twice the width of the coded sequence chip, and it is symmetrically extracted to both sides of the maximum energy in-phase axis position. The extracted drill bit equivalent source wavelet is a short-time waveform containing tens to hundreds of sampling points. The number of sampling points is determined by the demodulator's pulse repetition frequency and the extraction duration.

[0041] Based on the real-time drill bit position parameters, the equivalent source wavelets extracted from the drill bit at different thrust modulation cycles are assigned to the corresponding virtual source point spatial coordinates along the drilling path. The real-time drill bit position parameters provide the three-dimensional coordinates of the drill bit in the tunnel coordinate system, and the drilling path is a spatial curve extending along the designed drilling trajectory. The virtual source point spatial coordinates represent the spatial position of the drill bit at the start of each thrust modulation cycle. The equivalent source wavelets and virtual source point spatial coordinates are stored in a one-to-one correspondence to form a virtual source point sequence, providing deterministic input with known source locations and known source time functions for reverse-time migration imaging.

[0042] In another preferred embodiment of the present invention, the process of performing reverse time migration imaging processing in step S4, using the drill bit equivalent source wavelet as the source function and the vibration signal as the observation wavefield, is as follows: The initial velocity model is a two-dimensional longitudinal profile velocity model extending from the tunnel face towards the drilling direction. The two-dimensional longitudinal profile is a vertical section along the drilling axis, with the horizontal direction corresponding to the drilling axis and the vertical direction perpendicular to the drilling axis. The horizontal extension length of the velocity model is 1.2 times the designed depth of the horizontal borehole, and the vertical extension height is twice the tunnel excavation height. The velocity model is discretized into multiple velocity grid nodes, with a horizontal grid spacing of 0.5 meters and a vertical grid spacing of 0.5 meters. Each velocity grid node stores a longitudinal wave propagation velocity value, representing the propagation speed of the seismic wave in the rock mass at that node location, in meters per second.

[0043] The P-wave propagation velocity values ​​for each velocity grid node in the initial velocity model were extrapolated from the wave velocity test values ​​of the exposed rock mass at the tunnel face along the drilling direction. These wave velocity values ​​were obtained by cross-hole wave velocity testing using ultrasonic probes deployed on the tunnel face surface, yielding P-wave propagation velocities ranging from 2500 m / s to 4500 m / s. The extrapolation method was as follows: starting from the tunnel face, the P-wave propagation velocity value increased or decreased by a gradient of 0.5% to 2% for every 1 meter extended along the drilling direction. The sign of the gradient was determined based on the statistical trend of wave velocity variation with depth at the corresponding depth segment in the regional geological data. In the extrapolated initial velocity model, the velocity grid nodes at the tunnel face were assigned the measured wave velocity values, the farthest velocity grid node in front of the tunnel face was assigned the extrapolated endpoint wave velocity value, and the wave velocity values ​​of intermediate velocity grid nodes were determined by linear interpolation.

[0044] The forward extrapolation wavefield calculation uses the location of the virtual source point as the source grid node. The spatial coordinates of the virtual source point are provided by the virtual source point sequence output in step S3. The row and column indices of the source grid nodes are determined by nearest neighbor matching between the spatial coordinates of the virtual source point and the coordinates of the velocity grid nodes. The drill bit equivalent source wavelet is loaded as a source time function into the right-hand side of the acoustic equation. The acoustic equation in a two-dimensional isotropic medium is in the form that the second-order partial derivative of the wave field with respect to time is equal to the square of the longitudinal wave propagation velocity multiplied by the sum of the second-order partial derivatives of the wave field with respect to the spatial coordinates. The discrete sampled values ​​of the drill bit equivalent source wavelet are added sequentially to the right-hand side of the acoustic equation corresponding to the source grid node location according to the time step.

[0045] Wavefield extrapolation employs a second-order temporal and higher-order spatial staggered-grid finite-difference scheme. The staggered grid means that wavefield values ​​and medium parameters are defined at different spatial locations within the grid; wavefield values ​​are defined at the centers of grid nodes, and P-wave propagation velocities are defined at grid nodes and their boundaries. Temporal discretization uses a second-order central difference scheme, with the time step determined by the stability condition based on the maximum P-wave propagation velocity in the velocity model and the grid spacing, ranging from 0.05 ms to 0.2 ms. Spatial discretization uses an eighth-order central difference scheme, with the spatial partial derivative approximated by a weighted sum of the wavefield values ​​from the surrounding eight grid nodes. Forward extrapolation starts at time zero, calculating the wavefield value for each velocity grid node step-by-step along the positive time direction. The total number of time steps is determined by the length of the vibration signal recording and the time step size.

[0046] The reverse-time extrapolation wavefield calculation uses the locations of the multi-channel vibration sensors as receiving grid nodes. The spatial coordinates of the receiving grid nodes are provided by a mapping table, and the row and column indices of the receiving grid nodes are determined by nearest neighbor matching with the velocity model grid nodes. The vibration signals recorded by the multi-channel vibration sensors are loaded into the right-hand side of the acoustic equation in the reverse time direction. The vibration signal loading method is as follows: at each time step of the reverse-time extrapolation, the amplitude value of the vibration signal at the corresponding moment is read, and this amplitude value is added to the right-hand side of the acoustic equation corresponding to the location of the receiving grid node. The reverse-time extrapolation starts from the end of the vibration signal recording and calculates the wavefield value of each velocity grid node in the reverse time direction step by step, using the same staggered grid finite difference scheme as the forward extrapolation.

[0047] Within each time step of the wavefield extrapolation, source-normalized cross-correlation imaging conditions are applied to both the forward and reverse-time extrapolated wavefields. The application process is as follows: at each velocity grid node, the wavefield value of the forward extrapolated wavefield at the current time step is multiplied by the wavefield value of the reverse-time extrapolated wavefield at the current time step to obtain the cross-correlation value for that time step. Simultaneously, the sum of squares of the wavefield values ​​of the forward extrapolated wavefield at each time step is calculated, and the reciprocal of the sum is used as a normalization factor. This normalization factor is multiplied by the cross-correlation value for each time step to obtain the normalized cross-correlation value. The normalized cross-correlation values ​​for all time steps are summed, and the sum is used as the reflection coefficient value for the velocity grid node. The reflection coefficient value is a dimensionless value representing the magnitude of the change in medium wave impedance at the location of that velocity grid node. The reflection coefficient values ​​of all velocity grid nodes constitute the reflection coefficient distribution.

[0048] In another preferred embodiment of the present invention, after the reflection coefficient values ​​of all velocity grid nodes constitute a reflection coefficient distribution in step S4, the process further includes using the reflection coefficient distribution to correct the longitudinal wave propagation velocity values ​​of each velocity grid node in the initial velocity model: Velocity grid nodes with reflection coefficient values ​​below a set background level are marked as background nodes. The background level is set to the 80th percentile of all reflection coefficient values ​​in the reflection coefficient distribution, from smallest to largest. For example, if the 80th percentile of the reflection coefficient value dataset is 0.02, then velocity grid nodes with reflection coefficient values ​​below 0.02 are marked as background nodes. Background nodes correspond to regions in the rock mass where the wave impedance changes insignificantly, and the P-wave propagation velocity values ​​at background nodes remain unchanged from the initial velocity model.

[0049] Velocity grid nodes with reflection coefficient values ​​exceeding a set background level are marked as interface influence nodes. Interface influence nodes correspond to areas with significant changes in acoustic impedance near lithological interfaces or geological anomaly boundaries. A one-dimensional wavenumber domain inversion is applied to the interface influence nodes along the drilling direction. The one-dimensional wavenumber domain inversion process is as follows: Using each interface influence node as the center, the reflection coefficient values ​​of five velocity grid nodes before and after it are extracted along the drilling direction, forming a one-dimensional reflection coefficient sequence of length 11. A Fourier transform is then performed on this one-dimensional reflection coefficient sequence along the depth direction (the drilling direction). The Fourier transform converts the reflection coefficient sequence from the spatial domain to the wavenumber domain, obtaining the reflection coefficient wavenumber spectrum. The reflection coefficient wavenumber spectrum represents the amplitude distribution of different spatial frequency components.

[0050] There is a linear relationship between the reflection coefficient wavenumber spectrum and the P-wave propagation velocity perturbation. This relationship is as follows: the P-wave propagation velocity perturbation and the reflection coefficient value are correlated in the wavenumber domain through a product factor, which is equal to half of the reference P-wave propagation velocity value. The reference P-wave propagation velocity value is taken as the P-wave propagation velocity value of the interface-affected node in the initial velocity model. Using this linear relationship, the reflection coefficient wavenumber spectrum is divided by the product factor to obtain the P-wave propagation velocity perturbation spectrum in the wavenumber domain. Then, an inverse Fourier transform is performed on the P-wave propagation velocity perturbation spectrum to obtain the P-wave propagation velocity correction at the location of the interface-affected node in the spatial domain. The P-wave propagation velocity correction is a positive or negative number. A positive number indicates that the P-wave propagation velocity value in the original initial velocity model is too low and needs to be adjusted upwards, while a negative number indicates that the P-wave propagation velocity value in the original initial velocity model is too high and needs to be adjusted downwards.

[0051] The correction amount for the P-wave propagation velocity at each interface-affected node location is added to the corresponding P-wave propagation velocity value in the initial velocity model to obtain the corrected P-wave propagation velocity value. After the P-wave propagation velocity values ​​of all interface-affected nodes are corrected, the corrected initial velocity model replaces the original initial velocity model. Subsequent forward extrapolation wavefield calculations and reverse time extrapolation wavefield calculations corresponding to the virtual source point use the corrected initial velocity model as the new initial velocity model, proceeding to the next round of reverse time migration imaging processing and velocity model correction.

[0052] In another preferred embodiment of the present invention, the process of obtaining the spatial offset of the lithological interface in step S5 is as follows: First, the gradient amplitude of the reflection coefficient distribution is extracted. The reflection coefficient distribution is a two-dimensional data array. The horizontal direction corresponds to the drilling axis, and the vertical direction corresponds to the vertical direction perpendicular to the drilling axis. Each element in the array is the reflection coefficient value of the corresponding velocity grid node. The physical meaning of the gradient amplitude is the degree of drastic change of the reflection coefficient value in space. The larger the gradient amplitude, the more rapid the change of the reflection coefficient value along the spatial direction, corresponding to the interface location where the rock mass wave impedance changes significantly.

[0053] For each velocity grid node in the reflection coefficient distribution, the horizontal gradient component along the drilling direction and the vertical gradient component perpendicular to the drilling direction are calculated using central difference. The horizontal gradient component is calculated as follows: the reflection coefficient value of the adjacent velocity grid node to the right of the current velocity grid node is subtracted from the reflection coefficient value of the adjacent velocity grid node to the left, and the difference is divided by twice the horizontal grid spacing. The vertical gradient component is calculated as follows: the reflection coefficient value of the adjacent velocity grid node above the current velocity grid node is subtracted from the reflection coefficient value of the adjacent velocity grid node below the current velocity grid node, and the difference is divided by twice the vertical grid spacing. The gradient magnitude of the current velocity grid node is obtained by square-squaring the horizontal and vertical gradient components, summing them, and then taking the square root.

[0054] After calculating the gradient magnitudes of all velocity grid nodes, a gradient magnitude distribution of the same size as the reflection coefficient distribution is obtained. Local extrema are then selected from the gradient magnitude distribution. The selection process involves comparing the gradient magnitude of each velocity grid node with the gradient magnitudes of its four neighboring velocity grid nodes (top, bottom, left, and right). When the gradient magnitude of the central velocity grid node is greater than the gradient magnitudes of its four neighboring nodes, that central velocity grid node is retained as a local extremum. The gradient level is set to the 90th percentile of the gradient magnitudes of all local extrema in the gradient magnitude distribution. For example, if the 90th percentile is 0.15, then local extrema with gradient magnitudes exceeding 0.15 are marked as interface candidate points, thus obtaining a set of interface candidate points.

[0055] The candidate points for the interface are distributed as scattered points in the spatial domain. The Hough transform maps these candidate points from the spatial domain to the tilt-intercept parameter space. The Hough transform is a technique for identifying line segments from scattered points. Its principle is as follows: In the spatial domain, countless lines can be drawn through any point, each line uniquely determined by two parameters: tilt angle and intercept. The tilt angle is defined as the angle between the line and the horizontal direction, and the intercept is defined as the coordinates of the intersection point of the line and the vertical axis. Substituting the coordinates of a point in the spatial domain into the equation of the line yields a functional relationship between the tilt angle and the intercept, which is represented as a curve in the tilt-intercept parameter space. Multiple points on the same line in the spatial domain will have their corresponding curves in the tilt-intercept parameter space converge at the same parameter point. The tilt angle and intercept values ​​at this parameter point are the tilt angle and intercept of the corresponding line in the spatial domain.

[0056] In practice, the tilt angle parameter is set to a range of -90 degrees to +90 degrees, with a tilt angle sampling interval of 0.5 degrees. The intercept parameter range is determined based on the spatial dimensions of the reflection coefficient distribution, and the intercept sampling interval is 0.5 times the vertical grid spacing. For each interface candidate point in the set of interface candidate points, its spatial coordinates are substituted into the linear equation. All tilt angle sampling values ​​are traversed, and the corresponding intercept value is calculated. One vote is accumulated at the corresponding parameter point position in the tilt angle-intercept parameter space. After voting for all interface candidate points is completed, parameter points with accumulated values ​​exceeding a set accumulation level are searched in the tilt angle-intercept parameter space. The accumulated level is set to one-tenth of the total number of interface candidate points. For example, if there are 200 interface candidate points, the accumulated level is set to 20, and parameter points with accumulated values ​​exceeding 20 are searched. The tilt angle and intercept values ​​of the searched parameter points are inversely transformed back to the spatial domain. The inverse transformation process involves drawing a straight line in the spatial domain based on the tilt angle and intercept values ​​to obtain interface line segments in the form of straight line segments.

[0057] The coordinates of the intersection point between the interface line segment and the drilling axis are extracted as the spatial position representation value of the interface line segment. The drilling axis is the projection line of the designed drilling trajectory in a two-dimensional longitudinal profile. The spatial position representation value of the interface line segment is a two-dimensional coordinate, which includes the horizontal distance along the drilling axis and the vertical offset from the drilling axis.

[0058] The interface positions corresponding to the azimuth and dip angles of the designed drilling trajectory are compared with the spatial positional representations of the interface line segments. The interface positions corresponding to the azimuth and dip angles of the designed drilling trajectory are pre-determined based on geological exploration data and geophysical interpretation results from the tunnel exploration phase. These interface positions include the distance along the drilling axis and two-dimensional coordinates in a plane perpendicular to the drilling axis. During the comparison, the distance along the drilling axis from the interface position corresponding to the azimuth and dip angles of the designed drilling trajectory to the spatial positional representation of the interface line segment is taken as the spatial offset of the lithological interface.

[0059] In another preferred embodiment of the present invention, in step S6, the process of generating a drilling direction correction angle based on the spatial offset of the lithological interface and superimposing the drilling direction correction angle onto the azimuth and dip angle of the designed drilling trajectory is as follows: The spatial offset of the lithological interface is decomposed into a lateral offset component perpendicular to the drilling axis and a longitudinal offset component along the drilling axis. The direction perpendicular to the drilling axis refers to two orthogonal directions simultaneously perpendicular to the drilling axis in three-dimensional space. The lateral offset component is the vector sum of the projections of the spatial offset of the lithological interface onto these two orthogonal directions. The longitudinal offset component along the drilling axis is the projected length of the spatial offset of the lithological interface along the drilling axis direction.

[0060] The tangent of the drilling direction correction angle is calculated as the ratio of the lateral offset component to the distance from the current real-time drill bit position parameter to the spatial position representation value of the lithology interface. The current real-time drill bit position parameter is provided in real time by the displacement sensor built into the drill pipe propulsion device. The distance to the spatial position representation value of the lithology interface is the straight-line distance from the current drill bit position along the drilling axis to the spatial position representation value of the interface line segment. The drilling direction correction angle is an angle defined in a plane perpendicular to the drilling axis, representing the angle by which the drill bit's current travel direction needs to deflect towards the actual position of the lithology interface.

[0061] The drilling direction correction angle is weighted and allocated based on the drilling rate of the drill pipe propulsion device and the tangent of the drilling direction correction angle. The drilling rate is the instantaneous propulsion speed of the drill bit along the drilling axis, obtained by dividing the displacement difference output by the displacement sensor by the time interval. The weighting process is as follows: the drilling direction correction angle is multiplied by a weighting coefficient, which is the drilling rate divided by the sum of the drilling rate and the tangent of the drilling direction correction angle. The weighted drilling direction correction angle is then decomposed into an azimuth correction component and a dip correction component. The azimuth correction component is the projection angle of the drilling direction correction angle in the horizontal plane, and the dip correction component is the projection angle of the drilling direction correction angle in the vertical plane.

[0062] The azimuth correction component is superimposed on the azimuth of the designed drilling trajectory, by directly adding the azimuth correction component to the current azimuth value of the designed drilling trajectory. Similarly, the dip correction component is superimposed on the dip of the designed drilling trajectory, by directly adding the dip correction component to the current dip value of the designed drilling trajectory. The superimposed azimuth and dip values ​​constitute the output parameters of the calibrated drilling trajectory coordinate sequence.

[0063] When the real-time drill bit position parameters indicate that the drill bit has entered a preset distance range ahead of the spatial position representation value of the lithology interface, the weighting coefficient of the drilling direction correction angle decreases linearly with the remaining distance. The preset distance range is set to 2 meters ahead of the spatial position representation value of the lithology interface. The remaining distance is the distance from the current drill bit position along the drilling axis to the spatial position representation value of the lithology interface. The linear decrease is achieved by multiplying the weighting coefficient by the remaining distance and dividing by the preset distance range. When the drill bit reaches the spatial position representation value of the lithology interface, the remaining distance is zero, the weighting coefficient decreases to zero, the drilling direction correction angle is reset to zero, and the drill bit resumes drilling according to the original azimuth and dip angles of the designed drilling trajectory.

[0064] Please see Figure 2 As shown, the present invention also includes a horizontal hole drilling-while-drilling sensing and calibration system for tunnel geological exploration, used to implement the above-described horizontal hole drilling-while-drilling sensing and calibration method for tunnel geological exploration, including: The sensor deployment module is used to deploy multiple vibration sensors along the tunnel face and the adjacent excavated tunnel wall surface before drilling the horizontal hole of the tunnel, and to connect the multiple vibration sensors to the data acquisition unit. The encoding modulation acquisition module is used to start the drilling rig, control the drill rod propulsion device to apply periodic thrust changes to the drill bit according to the complementary binary encoding sequence, synchronously record the thrust modulation time sequence, the real-time position parameters of the drill bit and the drilling time, and continuously record the vibration signals received by multiple vibration sensors. The pulse compression and extraction module is used to perform pulse compression processing on the vibration signal using the thrust modulation time sequence, extract the drill bit equivalent source wavelet corresponding to each thrust modulation cycle, and assign the drill bit equivalent source wavelet to the virtual source point corresponding to the real-time position parameters of the drill bit. The reverse time migration imaging module is used to perform reverse time migration imaging processing with the drill bit equivalent source wavelet as the source function and the vibration signal as the observation wave field to generate the reflection coefficient distribution of the area in front of the tunnel face. The interface recognition and extraction module is used to identify interface segments from the reflection coefficient distribution, extract the spatial position of the interface segments, and compare them with the interface positions corresponding to the azimuth and dip angles of the designed drilling trajectory to obtain the spatial offset of the lithological interface. The trajectory calibration output module is used to generate a drilling direction correction angle based on the spatial offset of the lithological interface, superimpose the drilling direction correction angle onto the azimuth and dip angle of the designed drilling trajectory, and output the calibrated drilling trajectory coordinate sequence.

[0065] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A method for horizontal borehole sensing and calibration while drilling for tunnel geological exploration, characterized in that, Includes the following steps: S1. Before drilling the horizontal hole of the tunnel, multiple vibration sensors are laid out along the tunnel face and the surface of the adjacent excavated tunnel wall, and the multiple vibration sensors are connected to the data acquisition device. S2. Start the drilling rig and control the drill rod propulsion device to apply periodic thrust changes to the drill bit according to the complementary binary coding sequence. Simultaneously record the thrust modulation time sequence, real-time position parameters of the drill bit and drilling time, and continuously record the vibration signals received by multiple vibration sensors. S3. Use the thrust modulation time sequence to perform pulse compression processing on the vibration signal, extract the drill bit equivalent source wavelet corresponding to each thrust modulation cycle, and assign the drill bit equivalent source wavelet to the virtual source point corresponding to the real-time position parameters of the drill bit. S4. Using the drill bit equivalent source wavelet as the source function and the vibration signal as the observation wave field, reverse time migration imaging is performed to generate the reflection coefficient distribution in the area in front of the tunnel face. S5. Identify interface segments from the reflection coefficient distribution, extract the spatial position of the interface segments, and compare them with the interface positions corresponding to the azimuth and dip angles of the designed drilling trajectory to obtain the spatial offset of the lithological interface. S6. Generate the drilling direction correction angle based on the spatial offset of the lithological interface, superimpose the drilling direction correction angle onto the azimuth and dip angle of the designed drilling trajectory, and output the calibrated drilling trajectory coordinate sequence.

2. The horizontal borehole drilling sensing and calibration method for tunnel geological exploration according to claim 1, characterized in that, In step S1, the process of deploying multiple vibration sensors along the tunnel face and the adjacent excavated tunnel wall surface, and connecting the multiple vibration sensors to the data acquisition unit, is as follows: The multiple vibration sensors are distributed acoustic wave sensing optical cables. The distributed acoustic wave sensing optical cables are pasted along the tunnel face to form a spiral shape. The outermost circle of the spiral is attached to the excavation outline of the tunnel face. The distributed acoustic wave sensing optical cables extend axially along the excavated tunnel wall and pass through the arch, the two sides of the arch waist and the two sides of the sidewall. The data acquisition unit is a distributed acoustic wave sensor demodulator. The distributed acoustic wave sensor demodulator emits probe light pulses to the distributed acoustic wave sensor optical cable and receives backscattered Rayleigh light signals. It demodulates the phase difference changes between adjacent spatial resolution units to obtain the vibration signals at each sensor channel position along the distributed acoustic wave sensor optical cable. The spatial coordinates of each sensor channel position of the distributed acoustic wave sensor optical cable are obtained through three-dimensional measurement, and a mapping table is established between the spatial coordinates and the sensor channel labels.

3. The horizontal borehole drilling sensing and calibration method for tunnel geological exploration according to claim 1, characterized in that, In S2, the process of controlling the drill pipe propulsion device to apply periodic thrust changes to the drill bit according to the complementary binary coding sequence is as follows: An electronically controlled proportional pressure reducing valve is installed on the hydraulic supply circuit of the drill pipe propulsion device. The electronically controlled proportional pressure reducing valve receives the first and second encoded sequences output by the encoder controller. The sum of the autocorrelation of the first and second encoded sequences is the impulse function. The side lobes of the autocorrelation function of the first and second encoded sequences have the same amplitude and opposite signs as the side lobes of the autocorrelation function of the second encoded sequence. The drill pipe propulsion device outputs thrust modulation according to the first encoding sequence in the first time period and outputs thrust modulation according to the second encoding sequence in the second time period. The high level in the encoding sequence corresponds to the thrust setting base value, and the low level corresponds to the thrust setting base value plus the preset thrust increase value. The thrust modulation period interval is an integer multiple of the width of the encoding sequence chip, and the width of the encoding sequence chip is set to a non-integer multiple of the drill bit rotation period.

4. The horizontal borehole drilling sensing and calibration method for tunnel geological exploration according to claim 1, characterized in that, In step S3, the process of using the thrust modulation time sequence to perform pulse compression processing on the vibration signal and extract the equivalent seismic source wavelet of the drill bit is as follows: The vibration signal is separated into a first vibration signal segment and a second vibration signal segment according to the encoding sequence. The first vibration signal segment corresponds to a first time period, and the second vibration signal segment corresponds to a second time period. The first vibration signal segment is cross-correlated with the first coded sequence copy to obtain the first correlation signal, and the second vibration signal segment is cross-correlated with the second coded sequence copy to obtain the second correlation signal. The first and second correlation signals are aligned along the time axis and then added point by point to obtain the pulse compression signal after sidelobe suppression. A bandpass filter is applied to the pulse compression signal to filter out the periodic vibration component caused by the drill bit rotation. The maximum energy phase axis is tracked along the thrust modulation time sequence in the bandpass filtered pulse compression signal, and the waveform segment at the time delay corresponding to the maximum energy phase axis is extracted as the drill bit equivalent source wavelet. The drill bit equivalent source wavelet is assigned to the corresponding virtual source point spatial coordinates on the drilling path according to the drill bit real-time position parameters.

5. The horizontal borehole drilling sensing and calibration method for tunnel geological exploration according to claim 1, characterized in that, In step S4, the process of performing reverse time migration imaging, using the drill bit equivalent source wavelet as the source function and the vibration signal as the observed wavefield, is as follows: The initial velocity model is a two-dimensional longitudinal profile velocity model that extends from the tunnel face in the drilling direction. The longitudinal wave propagation velocity values ​​of each velocity grid node in the initial velocity model are extrapolated from the wave velocity test values ​​of the exposed rock mass at the tunnel face along the drilling direction. Forward extrapolation wavefield calculation uses the location of the virtual source point as the source grid node. The equivalent source wavelet of the drill bit is loaded onto the right-hand side of the acoustic equation. The wavefield value of each velocity grid node is calculated along the forward time direction using an interleaved grid finite difference scheme. For reverse extrapolation wavefield calculation, the deployment location of multiple vibration sensors is used as the receiving grid node. The vibration signal is loaded onto the right-hand side of the acoustic equation along the reverse time direction. The wavefield value of each velocity grid node is calculated along the reverse time direction using an interleaved grid finite difference scheme. At each velocity grid node, the wavefield values ​​of the forward extrapolation wavefield and the reverse extrapolation wavefield at the same time are multiplied and summed over all time steps. The summation result is used as the reflection coefficient value of the velocity grid node. The reflection coefficient values ​​of all velocity grid nodes constitute the reflection coefficient distribution.

6. The horizontal borehole drilling sensing and calibration method for tunnel geological exploration according to claim 1, characterized in that, In step S4, after the reflection coefficient values ​​of all velocity grid nodes form a reflection coefficient distribution, the process also includes using the reflection coefficient distribution to correct the longitudinal wave propagation velocity values ​​of each velocity grid node in the initial velocity model: Velocity grid nodes with reflection coefficient values ​​lower than the set background level are marked as background nodes, and the P-wave propagation velocity values ​​at the background nodes are maintained as the P-wave propagation velocity values ​​of the initial velocity model. Velocity grid nodes with reflection coefficient values ​​exceeding the set background level are marked as interface influence nodes. One-dimensional wavenumber domain inversion is applied to the interface influence nodes along the drilling direction to solve for the P-wave propagation velocity correction of the interface influence nodes. The corrected initial velocity model replaces the original initial velocity model. The subsequent forward extrapolation wavefield calculation and reverse time extrapolation wavefield calculation corresponding to the virtual source point adopt the corrected initial velocity model.

7. The horizontal borehole drilling sensing and calibration method for tunnel geological exploration according to claim 1, characterized in that, In S5, the process of obtaining the spatial offset of the lithological interface is as follows: Extract the gradient magnitude distribution of the reflection coefficient values ​​from the reflection coefficient distribution, mark the positions where the gradient magnitude exceeds the set gradient level as interface candidate points, apply Hough transform to the interface candidate points, search for parameter points in the tilt-intercept parameter space where the cumulative value exceeds the set cumulative level, and transform the searched parameter points back to the spatial domain to obtain interface line segments in the form of straight line segments. The coordinates of the intersection of the interface line segment and the drilling axis are extracted as the spatial position representation value of the interface line segment. The interface position corresponding to the azimuth and dip angle of the designed drilling trajectory is compared with the spatial position representation value of the interface line segment. The distance from the interface position corresponding to the azimuth and dip angle of the designed drilling trajectory to the spatial position representation value of the interface line segment along the drilling axis is taken as the spatial offset of the lithological interface.

8. The horizontal borehole drilling sensing and calibration method for tunnel geological exploration according to claim 7, characterized in that, The process of extracting the gradient magnitude distribution of reflection coefficient values ​​from the reflection coefficient distribution and marking locations where the gradient magnitude exceeds a set gradient level as interface candidate points is as follows: For each velocity grid node in the reflection coefficient distribution, the horizontal gradient component along the drilling direction and the vertical gradient component perpendicular to the drilling direction are calculated using central difference. The gradient magnitude is obtained by taking the square root of the sum of the squares of the horizontal and vertical gradient components. The gradient magnitude of each velocity grid node is compared with the gradient magnitude of the adjacent velocity grid nodes. When the gradient magnitude of a velocity grid node is greater than that of the adjacent velocity grid nodes, the velocity grid node is retained as a local extremum point. The velocity grid nodes with gradient magnitudes exceeding the set gradient level among the local extremum points are marked as interface candidate points, thus obtaining a set of interface candidate points.

9. The horizontal borehole drilling sensing and calibration method for tunnel geological exploration according to claim 1, characterized in that, In step S6, the process of generating a drilling direction correction angle based on the spatial offset of the lithological interface and then superimposing the drilling direction correction angle onto the azimuth and dip angle of the designed drilling trajectory is as follows: The spatial offset of the lithological interface is decomposed into a lateral offset component perpendicular to the drilling axis and a longitudinal offset component along the drilling axis. The ratio of the lateral offset component to the distance from the current real-time position parameter of the drill bit to the spatial position characterization value of the lithological interface is taken as the tangent value of the drilling direction correction angle. The drilling direction correction angle is decomposed into an azimuth correction component and an inclination correction component. The azimuth correction component is superimposed on the azimuth of the designed drilling trajectory, and the inclination correction component is superimposed on the inclination of the designed drilling trajectory. When the real-time position parameters of the drill bit indicate that the drill bit has entered the preset distance range in front of the spatial position characterization value of the lithology interface, the drilling direction correction angle decreases linearly with the remaining distance. When the drill bit reaches the spatial position characterization value of the lithology interface, the drilling direction correction angle decreases to zero.

10. A horizontal hole drilling-while-drilling sensing and calibration system for tunnel geological exploration, used to implement the horizontal hole drilling-while-drilling sensing and calibration method for tunnel geological exploration as described in any one of claims 1-9, characterized in that, include: The sensor deployment module is used to deploy multiple vibration sensors along the tunnel face and the adjacent excavated tunnel wall surface before drilling the horizontal hole of the tunnel, and to connect the multiple vibration sensors to the data acquisition unit. The encoding modulation acquisition module is used to start the drilling rig, control the drill rod propulsion device to apply periodic thrust changes to the drill bit according to the complementary binary encoding sequence, synchronously record the thrust modulation time sequence, the real-time position parameters of the drill bit and the drilling time, and continuously record the vibration signals received by multiple vibration sensors. The pulse compression and extraction module is used to perform pulse compression processing on the vibration signal using the thrust modulation time sequence, extract the drill bit equivalent source wavelet corresponding to each thrust modulation cycle, and assign the drill bit equivalent source wavelet to the virtual source point corresponding to the real-time position parameters of the drill bit. The reverse time migration imaging module is used to perform reverse time migration imaging processing with the drill bit equivalent source wavelet as the source function and the vibration signal as the observation wave field to generate the reflection coefficient distribution of the area in front of the tunnel face. The interface recognition and extraction module is used to identify interface segments from the reflection coefficient distribution, extract the spatial position of the interface segments, and compare them with the interface positions corresponding to the azimuth and dip angles of the designed drilling trajectory to obtain the spatial offset of the lithological interface. The trajectory calibration output module is used to generate a drilling direction correction angle based on the spatial offset of the lithological interface, superimpose the drilling direction correction angle onto the azimuth and dip angle of the designed drilling trajectory, and output the calibrated drilling trajectory coordinate sequence.

Images