Tunnel geological interface dynamic inversion method and system based on seismic while drilling

Abstract

The present application belongs to the field of passive source seismic interference and velocity structure inversion technology, and specifically discloses a tunnel geological interface dynamic inversion method and system based on seismic while drilling, which comprises the following steps: picking up direct surface wave travel time curves based on virtual shot gathers and calculating equivalent surface wave velocities in the observation array; calculating observation coherence curves along the same phase axis corresponding to the equivalent surface wave velocities in the observation array; calculating the normalized error between the observation coherence curves and the simulated coherence curves, and the normalized error between the observation direct surface wave travel time curves and the simulated direct surface wave travel time curves, respectively, to form a single observation multi-target error function; accumulating the multi-target error functions of multiple continuous observations and performing global optimization to obtain the front velocity structure and geological interface position which can simultaneously explain the data of multiple observation positions. The present application can reduce the misjudgment of fault or geological interface distance caused by velocity change and improve the stability and reliability of tunnel advance prediction.

Description

Technical Field

[0001] This invention relates to the fields of passive source seismic interference and velocity structure inversion technology, and particularly to a method and system for dynamic inversion of tunnel geological interfaces based on seismic during drilling. Background Technology

[0002] The statements in this section are merely background information related to the present invention and do not necessarily constitute prior art.

[0003] Seismic inversion is one of the core technical means for advanced geological prediction of tunnels. By studying the propagation law of seismic waves in the rock mass, it infers the velocity structure, interface location, and physical property distribution of the medium in front of the tunnel face, thereby identifying adverse geological bodies such as fault fracture zones, abrupt transitions between soft and hard rocks, water-rich zones, boulders, and underground obstacles.

[0004] However, seismic inversion is inherently a highly uncertain nonlinear problem, and the inversion results suffer from severe ambiguity. Existing tunnel seismic advance prediction methods often approximate the velocity of the medium ahead as constant, or directly perform time-distance conversion using the equivalent velocity from a single observation, further exacerbating the non-uniqueness between velocity and interface distance, leading to incorrect estimation of fault or geological interface locations.

[0005] Continuous passive source seismic methods provide an important data source for velocity-structure detection ahead of tunnels. This method uses seismic interferometry to transform continuous noise without a definite excitation time into a virtual shot gather similar to an active source, and then performs velocity-structure inversion based on the travel time of direct surface waves and the coherence information of reflected waves. However, relying solely on passive source seismic data still cannot avoid the problem of multiple solutions in the inversion. Problems such as the unclear source time of passive sources, unknown source wavelets, and dynamic changes in energy with construction conditions mean that the constraint of a single type of data error function on the velocity model is limited, and different velocity models may still show similar goodness of fit under a single error index. Summary of the Invention

[0006] To address the aforementioned issues, this invention proposes a dynamic inversion method and system for tunnel geological interfaces based on seismic data during drilling. It introduces multi-objective joint error constraints, constructing a single-observation multi-objective error function by calculating the normalized error between the observed coherence curve and the simulated coherence curve, as well as the normalized error between the observed direct surface wave travel time curve and the simulated direct surface wave travel time curve. Furthermore, the error functions from multiple consecutive observations are accumulated and globally optimized to obtain a forward velocity structure and geological interface location that can simultaneously explain data from multiple observation locations.

[0007] In some implementations, the following technical solutions are adopted: A dynamic inversion method for tunnel geological interfaces based on seismic drilling includes: A geophone array is installed on the sidewall of the tunnel section that has been constructed to acquire and preprocess continuous passive source vibration signals generated by tunneling equipment or construction disturbances, and to construct a virtual shot collection. The travel time curves of direct surface waves are picked up based on the virtual shot set, and the equivalent surface wave velocity within the observation array is calculated; the observation coherence curve is calculated along the in-phase axis corresponding to the equivalent surface wave velocity within the observation array. A candidate forward velocity structure model is constructed and forward modeling is performed to obtain the simulated direct surface wave travel time curve and the simulated coherence curve; wherein, the candidate forward velocity structure model includes at least the Rayleigh wave slowness function propagating along the tunnel wall, the S-wave slowness function of the medium ahead of the tunnel, and one or more candidate geological interface locations; The normalized error between the observed coherence curve and the simulated coherence curve, and the normalized error between the observed direct surface wave travel time curve and the simulated direct surface wave travel time curve are calculated separately to form a single observation multi-target error function. By accumulating the multi-target error functions from multiple consecutive observations and performing global optimization, we can obtain the forward velocity structure and geological interface location that can simultaneously explain data from multiple observation locations.

[0008] As a further approach, continuous passive source vibration signals generated by tunneling equipment or construction disturbances are acquired, preprocessed, and used to construct a virtual shot collection, specifically: The continuous passive source vibration signal is divided into multiple short time windows, and each short time window is preprocessed. The detector closest to the tunneling end is selected as the reference channel. In the nth observation, the cross-correlation function between the reference channel and other channels within each short time window is calculated. The cross-correlation results of all short time windows are superimposed to obtain the virtual shot set of the nth observation.

[0009] As a further approach, the travel time curves of direct-reaching surface waves are picked up based on virtual shot collections, and the equivalent surface wave velocity within the observation array is calculated, specifically as follows: By picking up the strongest phase arrival time of the direct surface wave in the virtual shot concentration, an approximate linear relationship between the strongest phase arrival time and the observation slowness is constructed; The observation slowness is obtained by fitting the data using the least squares method, and then the reciprocal is taken to obtain the equivalent surface wave velocity within the observation array.

[0010] As a further approach, the observation coherence curve is calculated along the in-phase axis corresponding to the equivalent surface wave velocity within the observation array, specifically as follows: The observed coherence curves are obtained by weighted superposition or normalized coherent superposition of the moving time window along the slope of the direct surface wave travel time curve in the virtual shot set. The length of the moving time window is determined based on the main frequency or effective frequency band of the passive source signal.

[0011] As a further approach, a candidate forward velocity structure model is constructed and forward modeling is performed to obtain the simulated direct surface wave travel time curve and the simulated coherence curve, specifically: For the candidate forward velocity structure model, calculate the arrival time of the simulated direct surface wave. When the reflected wave arrives ; The simulated direct surface wave arrival time This refers to the simulated direct surface wave travel time curve, specifically: ; When the reflected wave arrives Specifically: ; in, This represents the slowness function of Rayleigh waves propagating along the tunnel wall. ; Indicates the slowness of the S-wave in the medium ahead of the tunnel. ; This represents the position coordinates of the i-th detector along the tunnel axis. Represents the equivalent velocity function of Rayleigh surface waves propagating along the tunnel wall; The S-wave velocity function representing the medium ahead of the tunnel; This indicates the location of the tunnel end, face, or equivalent reflection conversion point corresponding to the nth observation; This indicates the location of the k-th candidate forward reflection interface or velocity anomaly interface; A multichannel seismic record is synthesized using candidate wavelets. The simulated virtual record is as follows: ; Using the same coherence calculation operator as the observed data Quantized into a simulated coherence curve, we obtain: ; in, Indicates the arrival time of simulated direct surface waves. The candidate wavelet response centered on the i-th trace is used to simulate the direct surface wave event. This indicates the time when the simulated reflected wave from the k-th candidate reflecting interface arrives. The candidate wavelet response centered on the i-th channel is used to simulate the reflection event of the i-th channel; The amplitude coefficient representing a direct surface wave event; This represents the amplitude coefficient of the reflection event at the k-th candidate reflecting interface; Indicates the delay of the Mth candidate; This indicates the simulated record corresponding to the nth observation in the candidate time delay. The coherence value is calculated at [location].

[0012] As a further approach, a single-observation multi-target error function is formulated, specifically as follows: ; in, Represents the multi-objective error function. This represents the normalized error between the observed coherence curve and the simulated coherence curve. This represents the normalized error between the observed direct surface wave travel time curve and the simulated direct surface wave travel time curve. This represents the error weighting coefficient of the direct surface wave travel time curve.

[0013] As a further approach, the multi-target error function from multiple consecutive observations is accumulated and globally optimized, specifically as follows: The objective function to be optimized is: ; ; in, The weight for the nth observation; This represents the parameters of the model to be inverted, including piecewise S-wave slowness, Rayleigh wave slowness, interface location, and optional reflection amplitude parameters; These are the optimized model parameters; A global optimization algorithm is used to solve the objective function.

[0014] In other embodiments, the following technical solutions are adopted: A dynamic inversion system for tunnel geological interfaces based on seismic drilling includes: The seismic interferometry module is configured to deploy a geophone array on the sidewall of the tunnel section that has been constructed, acquire and preprocess continuous passive source vibration signals generated by tunneling equipment or construction disturbances, and construct a virtual shot collection. The observation feature extraction module is configured to pick up the travel time curve of the direct surface wave based on the virtual shot set and calculate the equivalent surface wave velocity within the observation array; and calculate the observation coherence curve along the in-phase axis corresponding to the equivalent surface wave velocity within the observation array. The velocity structure forward modeling module is configured to construct candidate forward velocity structure models and forward model them to obtain simulated direct surface wave travel time curves and simulated coherence curves; wherein, the candidate forward velocity structure model includes at least the Rayleigh wave slowness function propagating along the tunnel wall, the S-wave slowness function of the medium ahead of the tunnel, and one or more candidate geological interface locations; The error function calculation module is configured to calculate the normalized error between the observed coherence curve and the simulated coherence curve, and the normalized error between the observed direct surface wave travel time curve and the simulated direct surface wave travel time curve, forming a single observation multi-target error function. The joint inversion optimization module is configured to accumulate the multi-objective error functions of multiple consecutive observations and perform global optimization to obtain the forward velocity structure and geological interface location that can simultaneously explain data from multiple observation locations.

[0015] In other embodiments, the following technical solutions are adopted: A terminal device includes a processor and a memory, the processor being used to implement instructions; the memory being used to store multiple instructions, the instructions being adapted to be loaded and executed by the processor to perform the above-described dynamic inversion method for tunnel geological interfaces based on seismic drilling.

[0016] In other embodiments, the following technical solutions are adopted: A computer-readable storage medium storing a plurality of instructions adapted for loading and execution by a processor of a terminal device of the above-described dynamic inversion method for tunnel geological interfaces based on seismic drilling.

[0017] Compared with the prior art, the beneficial effects of the present invention are: This invention unifies the travel time curve error of direct surface waves and the reflection coherence curve error into a multi-objective error function, enabling the velocity structure and reflection interface to be constrained simultaneously. By combining the error function through multiple observations, the instability of the inverse problem of a single observation is reduced, allowing the velocity model and interface position to gradually converge. This avoids simply treating the velocity structure ahead of the tunnel as an invariant constant, reducing misjudgments of fault or geological interface distances caused by velocity changes, and improving the stability and reliability of tunnel advance prediction.

[0018] Other features and advantages of the invention will be set forth in part in the description which follows, and in part will be obvious from the description, or may be learned by practice of the invention. Attached Figure Description

[0019] Figure 1 This is a flowchart of the dynamic inversion method for tunnel geological interfaces based on seismic drilling in an embodiment of the present invention. Detailed Implementation

[0020] It should be noted that the following detailed description is illustrative and intended to provide further explanation of the invention. Unless otherwise specified, all technical and scientific terms used in this invention have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0021] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the scope of exemplary embodiments according to the invention. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, devices, components, and / or combinations thereof.

[0022] Example 1 Potential adverse geological risks exist ahead of the tunnel face, including lithological changes, weak interlayers, fault fracture zones, water-rich structures, or abrupt changes in rock mass integrity. To obtain lithological information ahead of the tunnel face without additional construction time, this embodiment utilizes drilling disturbances during advanced geological drilling as a passive source, uses direct surface waves to constrain the equivalent propagation velocity near the tunnel wall, uses reflection coherence curves to constrain the position of the leading interface and the slow S-wave structure, and reduces velocity-distance non-uniqueness through multiple observations and joint optimization, thereby improving the stability and reliability of tunnel advance prediction.

[0023] As a specific implementation, in one or more embodiments, a dynamic inversion method for tunnel geological interfaces based on seismic data during drilling is disclosed, combined with... Figure 1 Specifically, it includes: S101: Install a geophone array on the sidewall of the tunnel section that has been constructed to acquire and preprocess continuous passive source vibration signals generated by tunneling equipment or construction disturbances, and construct a virtual shot collection.

[0024] In this embodiment, 12–48 geophones are axially arranged on the sidewall or surface of the tunnel segments of the constructed section. This embodiment selects 24 geophones, with a spacing of 0.5–1.0 m. The geophones can be horizontal component or three-component sensors, preferably perpendicular to the tunnel wall direction to enhance Rayleigh wave energy reception. During each observation, 10–30 minutes of continuous tunneling vibration are recorded, with the sampling rate set according to the target frequency band.

[0025] Record the position of the tunneling end during observation. Coordinates of each detector The sampling rate, sensor orientation, observation period, and construction status are all considered. If the tunneling end moves forward, the geophone array is repeatedly deployed or moves forward with the already constructed section to form multiple continuous observations.

[0026] Here, a one-dimensional coordinate system is established along the tunnel axis. The direction of tunneling is Positive direction. No. The location of the tunnel end at the time of the next observation was: Side wall layout A detector, its coordinates are... Select a detector located near the tunnel end. Alternatively, the detector with the highest signal-to-noise ratio can be used as the reference channel. Continuous passive source recording is denoted as... ,in, .

[0027] The data obtained above is preprocessed, including: dividing the continuous passive source records into fixed time windows and removing windows with obvious saturation, impact anomalies, or construction intervals; performing mean removal, detrending, bandpass filtering, spectral whitening, and 1-bit normalization on the retained short windows; selecting the geophone closest to the tunneling end as the reference channel and cross-correlating the other channels with the reference channel; and superimposing all the cross-correlation results of the short windows to obtain the virtual shot set for each observation.

[0028] Specifically, the continuous passive source records are divided into Several short-time windows, each with a length of [missing information]. .

[0029] For each short time window, mean removal, detrending, bandpass filtering, spectral whitening, 1-bit normalization, and outlier window removal are performed to obtain the preprocessed records. The above preprocessing procedures are all relatively conventional technical methods, so they will not be described in detail.

[0030] In the In the second observation, the first The cross-correlation function between the Dao and the reference Dao is: (1) in, For cross-correlation delay, i.e. the first The amount of time delay of the channel signal relative to the reference channel signal.

[0031] By superimposing all short-time-window cross-correlation results, random noise is weakened, and coherent energy associated with stable propagation paths is enhanced. The resulting... Constituting the first Virtual shot collection for the second observation.

[0032] In the frequency domain, cross-correlation can also be expressed as: (2) in, This represents the frequency domain response of the reference channel. Indicates the first The frequency domain response of the channel; This indicates complex conjugation. This formula represents the reference path and the first... The frequency domain responses of the channels are multiplied and conjugate, thereby eliminating the influence of unknown source triggering time and some common propagation paths.

[0033] Through the above processing, construction noise without a clear excitation time is converted into an interpretable seismic record with the reference trace as a virtual source, providing input for subsequent direct surface wave moveout curve picking and reflection coherence curve calculation.

[0034] S102: Based on the virtual shot collection, pick up the direct surface wave travel time curve and calculate the equivalent surface wave velocity in the observation array; calculate the observation coherence curve along the in-phase axis corresponding to the equivalent surface wave velocity in the observation array.

[0035] In this embodiment, when the strongest phase of the direct-reaching surface wave is picked up by the virtual gun concentration, it is denoted as... For local arrays processed with narrow bandwidth, the travel time curve of the direct surface wave can be approximated as a linear relationship: (3) in, To observe the slowness or the slope of the in-phase axis, The intercept is obtained through least squares fitting. : (4) in, This indicates that in the nth observation, relative to the reference detector... The position of the i-th detector Due to the slowness of observation The resulting linear time increment; This indicates that during the least squares fitting process, candidate slowness p and relative position are used. The calculated predicted travel time increment; 'a' represents the intercept of the linear travel time curve or the reference zero time. This represents the position coordinates of the i-th detector along the array direction; This indicates the position coordinates of the reference detector.

[0036] No. The equivalent surface wave velocity within the observation array for this observation is: (5) This velocity is not the layered velocity obtained from the complete surface wave dispersion inversion, but rather the local equivalent propagation velocity within a specified effective frequency band and detector array range, used to determine the in-phase axis slope of subsequent coherent scans.

[0037] Since the reflected wave, after returning to the tunnel wall, has the same or similar propagation speed along the array as the direct surface wave, the reflection event can be scanned along the slope determined by equation (3). For candidate time delays... Construct a moving in-phase axis: (6) Candidate delay It is not obtained directly from a single channel, but rather by scanning within the allowed reflection time range of the observation record at time sampling intervals or with a given step size. For example: first, based on the arrival time curve of the direct surface wave, the total recording duration, the effective signal-to-noise ratio time window, and the estimated forward detection distance, the possible time delay range of the reflection event [τ] is determined. min ,τ max Then, candidate time delay sequences are generated according to the sampling interval Δt or several times the sampling interval Δτ. Where τ min The time after the effective window of the direct surface wave ends or the minimum resolvable reflection delay can be taken; τ max It can be determined by the recording length, the maximum estimated detection distance, and the minimum reasonable S-wave velocity; Δτ can be the sampling interval Δt, half a main cycle, or the scanning step size set by the engineering.

[0038] A half-window length is set near the phase axis. or Sampling points were used, and a Hamming window was employed. Weighted.

[0039] The normalized coherence curve in an engineering implementation can be written as: (7) in, For sampling points, The sampling interval is... To prevent extremely small positive numbers with a denominator of zero. This represents the amplitude value recorded by the virtual shot gather of the i-th channel in the n-th observation at a sampling point offset by q points near the candidate phase axis.

[0040] The physical meaning of equation (7) is: if multiple records can be superimposed along the same slope at the candidate time delay, the molecular superposition energy is high and the coherence value is close to 1; if it is just random noise or incoherent events, the coherence value is low. The numerator of equation (7) represents the coherence energy after multi-channel weighted superposition along the candidate in-phase axis, that is, first superimpose each record at... Take values ​​from nearby areas, then apply window weights. The weighted summation is performed across channels, and the summation result is squared. If a real reflection event exists at the candidate time delay, the multi-channel signals can be superimposed in phase, and the numerator will increase significantly. The denominator of equation (7) represents the total energy or normalized energy benchmark of each channel within the same window. That is, instead of superimposing across channels in phase, the energy of each channel is squared and summed, then multiplied by the number of channels Nr and a stability term ε is added. The denominator is used to offset the differences in amplitude, number of channels, and window energy, so that the coherence value mainly reflects the degree of in-phase consistency, rather than simply reflecting the magnitude of the signal energy.

[0041] For all candidate delays Calculate formula (7) to obtain the first... The coherence curve of the observation: (8) The strongest phase of the direct Rayleigh wave is identified in the virtual gun focus, and picked up using equation (3). The slope of the in-phase axis is obtained by fitting using equation (4). and equivalent surface wave velocity .

[0042] Then, according to equation (6), a moving in-phase axis parallel to the direct wave is constructed, and the candidate time delay is... The normalized coherence of equation (7) is calculated point by point to obtain the observed coherence curve. .

[0043] If the data exhibits significant dispersion, it can be divided into multiple narrow frequency bands, and the above processing can be performed separately. The error can be calculated or weighted fusion can be applied to the results of different frequency bands.

[0044] S103: Construct candidate forward velocity structure models and perform forward modeling to obtain simulated direct surface wave travel time curves and simulated coherence curves; wherein, the candidate forward velocity structure model includes at least the Rayleigh wave slowness function propagating along the tunnel wall, the S-wave slowness function of the medium ahead of the tunnel, and one or more candidate geological interface locations.

[0045] In this embodiment, the forward velocity structure is not fixed as a constant during implementation, but is represented as a slowness function to be inverted.

[0046] The candidate forward velocity structure model is as follows: Let the slowness of the Rayleigh wave propagating along the tunnel wall be... The slowness of the S-wave in the medium ahead of the tunnel is ,but: , (9) For the This observation simulates the arrival time of direct surface waves. The result is obtained by integrating the slowness of the Rayleigh wave along the array: (10) Discrete implementation is as follows: (11) If the first The candidate reflective interfaces are located at Then the reflected wave simulates time for: (12) In equation (12), the first term represents the propagation time of the reflected energy within the array after it is converted into a surface wave propagating along the tunnel wall; the second term represents the round-trip time of the S-wave propagating from the tunneling end to the front interface and back to the vicinity of the tunneling end. If we consider... If the error is considered a constant, then the same Multiple can be matched and Therefore, the distance needs to be calculated by inverting the velocity structure rather than by converting a fixed velocity.

[0047] in, This represents the slowness function of Rayleigh waves propagating along the tunnel wall. ; Indicates the slowness of the S-wave in the medium ahead of the tunnel. ; This represents the position coordinates of the i-th detector along the tunnel axis. Represents the equivalent velocity function of Rayleigh surface waves propagating along the tunnel wall; The S-wave velocity function representing the medium ahead of the tunnel; This indicates the location of the tunnel end, face, or equivalent reflection conversion point corresponding to the nth observation; This indicates the location of the k-th candidate forward reflection interface or velocity anomaly interface; ΔT represents the length of the j-th spatial grid segment when the direct surface wave arrives in the discrete simulation; ΔT represents the additional time difference or reflection delay of the reflected wave arrival time relative to the arrival time of the direct wave; L represents the distance of the forward reflecting interface or velocity anomaly relative to the tunnel end.

[0048] To generate a simulated record, either a Ricker wavelet (or Ricker source wavelet) or a Gaussian wavelet (or Gaussian wavelet) can be used, which is consistent with the measured frequency band. To synthesize a multichannel seismic record, the first... The simulated virtual record is as follows: (13) in, The amplitude coefficient represents the amplitude of the simulated direct surface wave event and is used to control the energy of the direct wave in the simulated record. This represents the reflection event amplitude coefficient of the k-th candidate reflective interface, which can be used to characterize the reflection intensity of that interface; Indicates the arrival time of simulated direct surface waves. The wavelet response centered on the wavelet; Indicates the arrival time of the k-th candidate reflected wave The wavelet response centered on the wavelet.

[0049] right Using the same coherence calculation operator (7) as the observed data, the simulated coherence curve can be obtained: (14) in, These represent the normalized coherence values ​​calculated from the simulated record corresponding to the nth observation at the 1st, 2nd, ..., Mth candidate time delays, respectively; τ1 to τ M Indicates the candidate time-delay scan sequence; This represents the simulated coherence curve composed of the simulated coherence values ​​at all candidate time delays. This represents the delay of the Mth candidate.

[0050] The velocity structure to be inverted is parameterized as an S-wave slowness function in piecewise constant, piecewise linear, or B-spline form. and set one or more candidate interface locations. .

[0051] For any candidate model, calculate the arrival times of the simulated direct surface wave and the reflected wave according to equations (10) and (12), synthesize the simulated record according to equation (13), and then obtain the simulated coherence curve using the coherence operator consistent with the observation. .

[0052] S104: Calculate the normalized error between the observed coherence curve and the simulated coherence curve respectively. And the normalized error between the observed direct surface wave travel time curve and the simulated direct surface wave travel time curve. This forms a single-observation multi-target error function.

[0053] In this embodiment, introducing multi-objective joint error constraints is an effective way to reduce the ambiguity of inversion. A single error function usually only measures the degree of agreement between observed and simulated data from one aspect, such as only considering the travel time of reflected waves or only considering waveform residuals, and its constraint dimension on the velocity model is limited. If we can simultaneously utilize multiple observation features with different physical meanings and sensitive directions, such as the travel time of direct surface waves and the coherence energy of reflected waves, to construct independent error terms and form a joint objective function, then different error terms can impose constraints on different aspects of the velocity model, which can significantly compress the feasible solution space and reduce the non-uniqueness of coupling between velocity and distance.

[0054] To avoid the dominance of the error function by differences in energy levels between observations, the maximum value of the observed and simulated coherence curves is first normalized: (15) (16) in, This represents the observed coherence curve. This represents the simulated coherence curve. To prevent extremely small positive numbers with a denominator of zero.

[0055] No. The coherence error term for each observation is defined as the normalized... Norm: (17) The discrete form is: (18) If p=2, then it is the normalized root mean square difference; if p=2, then it is the difference. This makes it more robust to abnormal spikes.

[0056] coherence error The smaller the value, the closer the reflection peak delay, peak shape, and energy distribution simulated by the candidate velocity model are to the measured data.

[0057] Direct surface wave travel time curve error This is used to evaluate the interpretability of candidate models for the travel time curves of direct surface waves. In this embodiment, instead of comparing a single average velocity, the arrival time curves of direct waves for each trace are compared.

[0058] Specifically, in the In this observation, the travel time vector was observed. and simulated travel vector They are respectively: (19) (20) To eliminate the difference at reference zero time, the two travel time curves can be first subjected to relative time shift processing: (twenty one) (twenty two) in, This indicates the arrival time of the direct surface wave in the i-th channel of the n-th observation relative to the reference channel. Relative travel time; This indicates the arrival time of the simulated direct surface wave at the i-th channel relative to the reference channel under the candidate model corresponding to the n-th observation. Relative travel time; This indicates the position coordinates of the reference detector or the starting detector of the array.

[0059] This approach can eliminate the differences in reference zero time between different observation batches or forward modeling, so that the surface wave travel time error term mainly compares the differences in the shape of the travel time curve and the equivalent propagation velocity within the array.

[0060] Then the first Surface wave travel time error term of the second observation Defined as normalization Norm: (twenty three) Its discrete form is: (twenty four) When data quality is poor or only fast computation is required, the average velocity within the array can also be used as an approximation: (25) in, Indicates the first Equivalent surface wave velocity within the observation array for this observation It can be obtained by linear fitting of the simulated direct surface wave arrival time curve. Recommended value to .

[0061] Equations (23) and (24) preserve the complete travel time curve information and are applicable to situations where the velocity varies within the array range; Equation (25) is simpler to calculate and is applicable to situations where the direct wave is approximately linear.

[0062] No. The single-observation multi-objective error function for this observation is:

[0063] Where λ is the weight of the surface wave travel time error term. Since equations (18) and (24) have been normalized, It has a clear meaning of relative weight. If the direct wave pickup is stable and the bandwidth is narrow, a larger value can be chosen. If the direct wave is significantly affected by segment waves, mechanical noise, or dispersion, a smaller value can be used. ; The recommended value is between 0.1 and 3, or determined by cross-validation. In engineering implementation, this can be determined through validating the dataset, cross-validation, or consistency of multiple observations. .

[0064] S105: The multi-target error functions from multiple consecutive observations are accumulated and globally optimized to obtain the forward velocity structure and geological interface location that can simultaneously explain data from multiple observation locations.

[0065] Tunnel excavation is a dynamic process. As the tunnel face moves forward, the relative position of the geophone array and the geological interface ahead continuously changes, and the same geological interface is repeatedly observed at different excavation locations. This continuous series of observations provides natural redundancy and multi-perspective constraints for inversion. An incorrect velocity model might fit the data characteristics of a single observation, but it is difficult to simultaneously fit the results of multiple observations of the same interface at different excavation distances. Therefore, jointly accumulating and optimizing the error functions of multiple consecutive observations can further tighten the constraints in the temporal dimension, allowing the velocity structure and interface position to gradually converge to the true solution through mutual verification of multiple observations.

[0066] Specifically, assuming there are a total of Multiple observations. Since each observation corresponds to a different tunneling end position and a different array position, the same forward interface will be repeatedly observed with different time delays. To reduce the velocity-distance non-uniqueness of a single observation, all observation errors are jointly optimized: (27) Right now: (28) in, For the first The weight of each observation can be determined based on signal-to-noise ratio, reflectance peak sharpness, pickup confidence, or observation coverage. If the quality of each observation is similar, a weight can be selected. .

[0067] The objective function to be optimized is: (29) in, This represents the parameters of the model to be inverted, including piecewise S-wave slowness, Rayleigh wave slowness, interface location, and optional reflection amplitude parameters. Global optimization methods such as particle swarm optimization, genetic algorithms, simulated annealing, and grid search can be used to solve the problem and avoid local minima.

[0068] The key role of the joint error of multiple observations in this embodiment is that: an incorrect velocity structure may fit the fit and moveout of a certain observation, but it is difficult to fit multiple observations at different tunneling positions at the same time; by forcing the model to interpret all observations at the same time through equation (28), the misjudgment of interface position can be significantly reduced.

[0069] As a further implementation, model error can be used to evaluate the inversion accuracy in synthetic data or known verification sections. Let the true S-wave slowness be... The inverted S-wave slowness is The last interface position that can be constrained by reflection information is Then the model error is: (30) For on-site engineering, since the actual speed and structure are unknown, the following evaluation indicators can be used: (1) Coherence fitting error: Evaluate whether the reflection peak delay, peak shape and peak distribution match; (2) Travel time fitting error: Evaluate whether the travel time curves of the direct surface wave match; (3) Consistency of multiple observations: Whether the positions obtained by inversion at different observation locations of the same interface converge; (4) Candidate model dispersion: Select several candidate models with the smallest objective function and calculate the standard deviation of interface position and velocity; (5) Construction verification indicators: Compare the predicted interface location with the excavation exposure, borehole data, face sketch, spoil change or sudden change in tunneling parameters; (6) Advance prediction distance error: evaluate the spatial deviation between the predicted interface location and the actual revealed location.

[0070] After optimization, the optimal S-wave velocity structure is output. Slow-speed structure Candidate interface or fault location And the range of uncertainty at each interface. The location of abrupt velocity changes can serve as a candidate explanation for the location of soft and hard rock interfaces, fault fracture zone boundaries, or other geological anomalies.

[0071] Furthermore, the minimum error model can be used as the optimal model, and the objective function can be ranked first. The model is used as a set of feasible models, and its interface location distribution represents the inversion uncertainty.

[0072] Because this method updates the forward velocity structure and interface position simultaneously during the inversion process, instead of directly converting the reflection delay into distance using a fixed velocity, it can avoid the problem of predicting the fault location too close or too far under conditions of significant changes in velocity structure.

[0073] Example 2 In one or more embodiments, a dynamic inversion system for tunnel geological interfaces based on seismic drilling is disclosed, specifically including: The seismic interferometry module is configured to deploy a geophone array on the sidewall of the tunnel section that has been constructed, acquire and preprocess continuous passive source vibration signals generated by tunneling equipment or construction disturbances, and construct a virtual shot collection. The observation feature extraction module is configured to pick up the travel time curve of the direct surface wave based on the virtual shot set and calculate the equivalent surface wave velocity within the observation array; and calculate the observation coherence curve along the in-phase axis corresponding to the equivalent surface wave velocity within the observation array. The velocity structure forward modeling module is configured to construct candidate forward velocity structure models and forward model them to obtain simulated direct surface wave travel time curves and simulated coherence curves; wherein, the candidate forward velocity structure model includes at least the Rayleigh wave slowness function propagating along the tunnel wall, the S-wave slowness function of the medium ahead of the tunnel, and one or more candidate geological interface locations; The error function calculation module is configured to calculate the normalized error between the observed coherence curve and the simulated coherence curve, and the normalized error between the observed direct surface wave travel time curve and the simulated direct surface wave travel time curve, forming a single observation multi-target error function. The joint inversion optimization module is configured to accumulate the multi-objective error functions of multiple consecutive observations and perform global optimization to obtain the forward velocity structure and geological interface location that can simultaneously explain data from multiple observation locations.

[0074] It should be noted that the specific implementation methods of the above modules are exactly the same as those in Example 1, and will not be described in detail again.

[0075] Example 3 In one or more embodiments, a terminal device is disclosed, comprising a processor and a memory, wherein the processor is used to implement instructions; and the memory is used to store multiple instructions adapted to be loaded by the processor and executed by the processor to perform the dynamic inversion method for tunnel geological interfaces based on seismic drilling as described in Embodiment 1.

[0076] It should be understood that in this embodiment, the processor can be a central processing unit (CPU), or it can be other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. A general-purpose processor can be a microprocessor or any other conventional processor.

[0077] Memory may include read-only memory and random access memory, and provides instructions and data to the processor. A portion of the memory may also include non-volatile random access memory. For example, memory may also store information about the device type.

[0078] In the implementation process, each step of the above method can be completed by the integrated logic circuits in the processor hardware or by software instructions.

[0079] Example 4 In one or more embodiments, a computer-readable storage medium is disclosed, wherein a plurality of instructions are stored, the instructions being adapted to be loaded by a processor of a terminal device and executed by the method for dynamic inversion of tunnel geological interfaces based on seismic drilling as described in Embodiment 1.

[0080] While the specific embodiments of the present invention have been described above in conjunction with the accompanying drawings, this is not intended to limit the scope of protection of the present invention. Those skilled in the art should understand that various modifications or variations that can be made by those skilled in the art without creative effort based on the technical solutions of the present invention are still within the scope of protection of the present invention.

Claims

1. A dynamic inversion method for tunnel geological interfaces based on seismic drilling, characterized in that, include: A geophone array is installed on the sidewall of the tunnel section that has been constructed to acquire and preprocess continuous passive source vibration signals generated by tunneling equipment or construction disturbances, and to construct a virtual shot collection. The travel time curves of direct surface waves are picked up based on the virtual shot set, and the equivalent surface wave velocity within the observation array is calculated; the observation coherence curve is calculated along the in-phase axis corresponding to the equivalent surface wave velocity within the observation array. A candidate forward velocity structure model is constructed and forward modeling is performed to obtain the simulated direct surface wave travel time curve and the simulated coherence curve; wherein, the candidate forward velocity structure model includes at least the Rayleigh wave slowness function propagating along the tunnel wall, the S-wave slowness function of the medium ahead of the tunnel, and one or more candidate geological interface locations; The normalized error between the observed coherence curve and the simulated coherence curve, and the normalized error between the observed direct surface wave travel time curve and the simulated direct surface wave travel time curve are calculated separately to form a single observation multi-target error function. By accumulating the multi-target error functions from multiple consecutive observations and performing global optimization, we can obtain the forward velocity structure and geological interface location that can simultaneously explain data from multiple observation locations.

2. The method for dynamic inversion of tunnel geological interfaces based on seismic drilling as described in claim 1, characterized in that, Acquire and preprocess continuous passive source vibration signals generated by tunneling equipment or construction disturbances to construct a virtual shot set, specifically: The continuous passive source vibration signal is divided into multiple short time windows, and each short time window is preprocessed. The detector closest to the tunneling end is selected as the reference channel. In the nth observation, the cross-correlation function between the reference channel and other channels within each short time window is calculated. The cross-correlation results of all short time windows are superimposed to obtain the virtual shot set of the nth observation.

3. The method for dynamic inversion of tunnel geological interfaces based on seismic drilling as described in claim 1, characterized in that, Based on the virtual shot gather, the travel time curves of direct surface waves are picked up, and the equivalent surface wave velocity within the observation array is calculated. Specifically: By picking up the strongest phase arrival time of the direct surface wave in the virtual shot concentration, an approximate linear relationship between the strongest phase arrival time and the observation slowness is constructed; The observation slowness is obtained by fitting the data using the least squares method, and then the reciprocal is taken to obtain the equivalent surface wave velocity within the observation array.

4. The method for dynamic inversion of tunnel geological interfaces based on seismic drilling as described in claim 1, characterized in that, The observation coherence curve is calculated along the in-phase axis corresponding to the equivalent surface wave velocity within the observation array, specifically as follows: The observed coherence curves are obtained by weighted superposition or normalized coherent superposition of the moving time window along the slope of the direct surface wave travel time curve in the virtual shot set. The length of the moving time window is determined based on the main frequency or effective frequency band of the passive source signal.

5. The method for dynamic inversion of tunnel geological interfaces based on seismic drilling as described in claim 1, characterized in that, A candidate forward velocity structure model was constructed and forward modeling was performed to obtain the simulated direct surface wave travel time curve and the simulated coherence curve, specifically: For the candidate forward velocity structure model, calculate the arrival time of the simulated direct surface wave. When the reflected wave arrives ; The simulated direct surface wave arrival time This refers to the simulated direct surface wave travel time curve, specifically: ； When the reflected wave arrives Specifically: ； in, This represents the slowness function of Rayleigh waves propagating along the tunnel wall. ; Indicates the slowness of the S-wave in the medium ahead of the tunnel. ; This represents the position coordinates of the i-th detector along the tunnel axis. Represents the equivalent velocity function of Rayleigh surface waves propagating along the tunnel wall; The S-wave velocity function representing the medium ahead of the tunnel; This indicates the location of the tunnel end, face, or equivalent reflection conversion point corresponding to the nth observation; This indicates the location of the k-th candidate forward reflection interface or velocity anomaly interface; A multichannel seismic record is synthesized using candidate wavelets. The simulated virtual record is as follows: ； Using the same coherence calculation operator as the observed data Quantized into a simulated coherence curve, we obtain: ； in, Indicates the arrival time of simulated direct surface waves. The candidate wavelet response centered on the i-th trace is used to simulate the direct surface wave event. This indicates the time when the simulated reflected wave from the k-th candidate reflecting interface arrives. The candidate wavelet response centered on the i-th channel is used to simulate the reflection event of the i-th channel; The amplitude coefficient representing a direct surface wave event; This represents the amplitude coefficient of the reflection event at the k-th candidate reflecting interface; Indicates the delay of the Mth candidate; This indicates the simulated record corresponding to the nth observation in the candidate time delay. The coherence value is calculated at [location].

6. The method for dynamic inversion of tunnel geological interfaces based on seismic drilling as described in claim 1, characterized in that, The single-observation multi-target error function is formed as follows: ； in, Represents the multi-objective error function. This represents the normalized error between the observed coherence curve and the simulated coherence curve. This represents the normalized error between the observed direct surface wave travel time curve and the simulated direct surface wave travel time curve. This represents the error weighting coefficient of the direct surface wave travel time curve.

7. The method for dynamic inversion of tunnel geological interfaces based on seismic drilling as described in claim 1, characterized in that, The multi-target error function from multiple consecutive observations is accumulated and then globally optimized, specifically as follows: The objective function to be optimized is: ； ； in, The weight for the nth observation; This represents the parameters of the model to be inverted, including piecewise S-wave slowness, Rayleigh wave slowness, interface location, and optional reflection amplitude parameters; These are the optimized model parameters; A global optimization algorithm is used to solve the objective function.

8. A dynamic inversion system for tunnel geological interfaces based on seismic drilling, characterized in that, include: The seismic interferometry module is configured to deploy a geophone array on the sidewall of the tunnel section that has been constructed, acquire and preprocess continuous passive source vibration signals generated by tunneling equipment or construction disturbances, and construct a virtual shot collection. The observation feature extraction module is configured to pick up the travel time curve of the direct surface wave based on the virtual shot set and calculate the equivalent surface wave velocity within the observation array; and calculate the observation coherence curve along the in-phase axis corresponding to the equivalent surface wave velocity within the observation array. The velocity structure forward modeling module is configured to construct candidate forward velocity structure models and forward model them to obtain simulated direct surface wave travel time curves and simulated coherence curves; wherein, the candidate forward velocity structure model includes at least the Rayleigh wave slowness function propagating along the tunnel wall, the S-wave slowness function of the medium ahead of the tunnel, and one or more candidate geological interface locations; The error function calculation module is configured to calculate the normalized error between the observed coherence curve and the simulated coherence curve, and the normalized error between the observed direct surface wave travel time curve and the simulated direct surface wave travel time curve, forming a single observation multi-target error function. The joint inversion optimization module is configured to accumulate the multi-objective error functions of multiple consecutive observations and perform global optimization to obtain the forward velocity structure and geological interface location that can simultaneously explain data from multiple observation locations.

9. A terminal device comprising a processor and a memory, the processor for implementing instructions; the memory for storing multiple instructions, characterized in that, The instructions are adapted to be loaded by a processor and executed by the dynamic inversion method for tunnel geological interfaces based on seismic drilling as described in any one of claims 1-7.

10. A computer-readable storage medium storing a plurality of instructions, characterized in that, The instructions are adapted to be loaded and executed by the processor of the terminal device using the dynamic inversion method for tunnel geological interfaces based on seismic drilling as described in any one of claims 1-7.

Images