A cross-correlation seismic interferometry HSP method, storage medium, and electronic device

Abstract

The present invention discloses a cross-correlation seismic interference HSP method, storage medium, and electronic equipment, and relates to the field of signal processing technology. The method includes: based on the HSP acoustic wave reflection method of TBM construction, a reference unit is arranged on the TBM cutter head, and two acquisition units are arranged on one side of the surrounding rock on both sides of the tunnel wall; the signals collected by the two acquisition units are subjected to LMS minimum mean square error adaptive delay estimation to obtain the time delay between the two acquisition units; the signals collected by the two acquisition units are subjected to a delay differential method to offset the same-frequency interference noise different from the excavation direction to obtain direct waves and reflected waves separated in the time domain; the separated reflected waves are segmented, and the signals collected by the reference unit are subjected to cross-correlation processing and then vertically superimposed to obtain the distance between the reference unit and the poor geological body in front of the tunnel. The method of the present invention can improve the accuracy of HSP method prediction.

Description

Technical Field

[0001] The present invention relates to the field of signal processing technology, and in particular to a cross-correlation seismic interferometry (HSP) method, a storage medium, and an electronic device. Background Art

[0002] At present, the use of full-face tunnel boring machine (TBM) for excavation construction has the advantages of high efficiency, high cost performance, and environmental protection. The horizontal seismic profiling (HSP) prediction technology, which uses the rock-breaking vibration of the TBM cutter head as the excitation source, realizes data collection while the TBM is excavating, and realizes the prediction of the unfavorable geological body in front of the excavation tunnel, which can ensure the safety of the TBM construction tunnel without affecting the TBM construction efficiency. The HSP prediction method uses the TBM excavation noise wave as the source, receives the reflected waves and other coupled noises generated by the geological anomaly, and performs digital filtering, time domain analysis, spectrum analysis, correlation interference analysis, reflection imaging and other processing on the signal data to obtain the spatial position and range of the unfavorable geological body in front of the face, and complete the geological prediction work.

[0003] The cross-correlation seismic interferometry method is often used to process passive source seismic data. It can enhance seismic reflection signals, suppress interference noise, and thus improve the signal-to-noise ratio of seismic signals. It is suitable for the HSP method and is expected to improve the quality of tunnel seismic data. However, the seismic signal generated by the TBM during the excavation process is a continuous random vibration signal. The sensor in the tunnel receives an irregular signal similar to noise, with a low signal-to-noise ratio and severe frequency aliasing. However, the cross-correlation seismic interferometry method cannot handle the situation where the instrument receives a noise signal with the same frequency as the signal, resulting in fluctuations in the signal-to-noise ratio of the recovered signal, and this method cannot distinguish between the direct wave and the surrounding rock reflection wave in the excavation signal. Summary of the invention

[0004] The purpose of the present invention is to solve the problem that the cross-correlation seismic interferometry method cannot process the noise signal with the same frequency as the signal received by the instrument, resulting in fluctuations in the signal-to-noise ratio of the restored signal. The present invention proposes a cross-correlation seismic interferometry HSP method, which can effectively remove the same-frequency noise and improve the signal-to-noise ratio of the reflected wave of the excavation signal. It is of great significance for real-time and accurate prediction of the geological conditions ahead of TBM excavation and ensuring the safety and efficiency of tunnel TBM construction.

[0005] The method of the present invention comprises the following steps:

[0006] S1. Based on the HSP acoustic wave reflection method of TBM construction, two acquisition units x and d are arranged on one side of the surrounding rock on both sides of the tunnel wall, and the signals collected by the two acquisition units x and d are obtained, wherein the signals collected by the two acquisition units x and d both include: seismic signals mixed with noise;

[0007] S2, performing LMS minimum mean square error adaptive delay estimation on the signals collected by the two collection units x and d to obtain the delay between the signals collected by the two collection units;

[0008] S3, based on the time delay between the collected signals of the two acquisition units, the signals collected by the two acquisition units are used to offset the same-frequency interference noise that is different from the excavation direction using a time delay differential method, so as to obtain the seismic signal after removing the same-frequency noise from the collected signals of the two acquisition units;

[0009] S4. Segment the seismic signals after removing the same-frequency noise from the signals collected by the two acquisition units in S3, cross-correlate the two segmented seismic signals after removing the same-frequency noise to obtain the reflected wave frequency domain signal, vertically superimpose the segmented signals, and calculate the propagation time of the reflected wave to obtain the distance between the reference unit and the unfavorable geological body in front of the tunnel.

[0010] Furthermore, the seismic signal includes a direct wave signal and a reflected wave signal.

[0011] Furthermore, step S2 is specifically as follows:

[0012] S21, establish an adaptive filter, the two acquisition units x and d acquire signals x(n) and d(n), x(n) is used as the input signal of the adaptive filter, d(n) is used as the target signal, and an error signal is obtained:

[0013] e(n)=d(n)-y(n)

[0014]

[0015] Wherein, e(n) is the error signal, n represents the number of signal points, y(n) is the output signal of the adaptive filter with x(n) as the input signal, ω m (n) is the weight coefficient of the m-th order filter, and M is the order of the adaptive filter;

[0016] S22, the error signal e(n) is expressed as:

[0017]

[0018] ω m (n+1)=ω m (n)+μe(n)x(nm)

[0019] Among them, μ is the convergence step size;

[0020] The mean square error performance function is expressed as:

[0021] V(W)=E[e 2 (n)] = E[(d(n)-W T (n)X(n)) 2 ]

[0022] Where V(W) is the mean square error performance function, E[] represents the expectation, W(n)=[ω1(n),ω2(n),...,ω m (n)] T ,X(n)=[x(n-1),x(n-2),...,x(nm)] T ;

[0023] P = E[d(n)X(n)] is the cross-correlation function of the desired output and input vectors, R x =E[X(n)X T (n)] is the autocorrelation function of the input signal, and the mean square error is the quadratic form of the weight vector. When the filter converges, the mean square error V(W) tends to zero, and the weight coefficient reaches its maximum value. At this time, ω a (n) = R x -1 P,a∈[1,2,...,m];

[0024] S23, perform sinc interpolation on the weight coefficient to obtain the time delay:

[0025]

[0026]

[0027] Among them, τ is the time delay and n is the number of signal points.

[0028] Furthermore, step S3 is specifically as follows:

[0029] S31. Establish a time difference model. x(n) and d(n) can be expressed in the time domain as follows:

[0030] x(t)=s(t)+n(t)

[0031] d(t)=s(t-Δt)+n(t+Δt)

[0032] Among them, x(t) and d(t) are the time domain representations of the signals collected by the two acquisition units, s(t) represents the time domain representation of the seismic signal, n(t) represents the time domain representation of the noise signal, and Δt represents the time difference between the signals collected by the two acquisition units;

[0033] S32, subtract the signal x(t) from the signal d(t) after time delay processing to obtain the delayed differential data x0(t) of the signals collected by the two collection units:

[0034] x0(t)=x(t)-d(t-τ)

[0035] x0(t-τ)=s(t-Δt+τ)+n(t+Δt-τ)

[0036] Let τ = Δt, x0(t) = s(t) - s(t-2Δt), and solve the difference equation to obtain the denoised seismic signals s(t) and s(t-Δt). s(t) and s(t-Δt) are the seismic signals after removing the same-frequency noise from the signals x(n) and d(n).

[0037] Furthermore, the cross-correlation processing of the two segmented seismic signals after removing the same-frequency noise is specifically as follows:

[0038]

[0039] Among them, C(s x ,s d ,t) is the signal s x With signal s d The cross-correlation coefficient in the time domain, s x is the segmented seismic signal collected by acquisition unit x after removing the same-frequency noise, s d is the segmented seismic signal collected by acquisition unit d after removing the same-frequency noise, G(s x ,s d ,t) is the Green function, which represents the displacement response of the unit force in the time domain t generated by the acquisition units x and d, G(s x ,s d ,-t) is the Green's function, which represents the displacement response in the time domain -t produced by the unit force on the acquisition unit x and d;

[0040] Furthermore, C(s x ,s d ,t) is expressed in the frequency domain as:

[0041]

[0042] In the formula, C(s x ,s d ,ω) represents C(s x ,s d ,t) in the frequency domain, u1(s x ,ω) represents the seismic frequency domain signal received by the acquisition unit x; u2(s d ,ω) represents the seismic frequency domain signal received by the acquisition unit d; G(sx ,s d ,ω) represents G(s x ,s d ,t) is expressed in the frequency domain; * indicates the complex number is conjugated, || indicates the complex number is real part; P(ω) indicates the expression of the reflected wave p(t) in the frequency domain;

[0043]

[0044] Wherein, T is the collection time of the two collection units;

[0045] Vertically superimpose each signal of P(ω) to obtain the signal:

[0046]

[0047]

[0048] Where N is the total number of segments; floor is the number of rounded down when the data does not meet the segment length; T N Indicates the segment length; P n (ω,n) is the frequency domain signal of the reflected wave obtained by superimposing n segments P(ω). The propagation time of the reflected wave is calculated to obtain the distance between the reference unit and the unfavorable geological body in front of the tunnel.

[0049] The present invention also proposes a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the above-mentioned LMS-TDOA-based cross-correlation seismic interferometry HSP method is implemented.

[0050] The present invention also proposes an electronic device, comprising a processor and a memory, wherein the processor and the memory are interconnected, wherein the memory is used to store a computer program, the computer program comprises computer-readable instructions, and the processor is configured to call the computer-readable instructions to execute the above-mentioned LMS-TDOA-based cross-correlation seismic interferometry HSP method.

[0051] The beneficial effects brought by the technical solution provided by the present invention are:

[0052] The present invention uses the time delay estimated by the adaptive time delay estimation method for the time delay differential method, removes the same frequency noise received by the acquisition unit, and separates the TBM excavation direct wave and the formation reflection wave signal in the HSP seismic wave detection data. Combined with the signal segmentation superposition method, the shortcomings of the cross-correlation seismic interferometry method used for the HSP method are improved, and the formation reflection signal submerged by the TBM tunnel background noise is extracted, thereby improving the accuracy of the HSP method prediction. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 is a flow chart of a cross-correlation seismic interferometry HSP method based on LMS-TDOA according to an embodiment of the present invention;

[0054] Figure 2 It is a principle diagram of the HSP acoustic wave reflection method based on TBM construction in an embodiment of the present invention;

[0055] Figure 3 This is a simulation result diagram of delay estimation performed in an embodiment of the present invention when the sampling rate is set to 500 Hz, the FIR filter order is set to 32, and the collection unit spacing is 200 m;

[0056] Figure 4 The comparison between the excavation data calculated by the time-domain delay difference algorithm and the original data in the embodiment of the present invention is shown in the figure;

[0057] Figure 5 is a curve diagram showing the change of the signal-to-noise ratio of data with the number of superpositions when the embodiment of the present invention processes a tunneling signal with a signal-to-noise ratio of 0 dB;

[0058] Figure 6 It is a block diagram of an electronic device in an exemplary embodiment of the present invention. DETAILED DESCRIPTION

[0059] To make the objectives, technical solutions and advantages of the present invention more clear, the embodiments of the present invention will be further described below with reference to the accompanying drawings.

[0060] The flowchart of the cross-correlation seismic interferometry HSP method based on LMS-TDOA in the embodiment of the present invention is as follows: Figure 1 , specifically including the following steps:

[0061] S1. Based on the HSP acoustic wave reflection method of TBM construction, two acquisition units x and d are arranged on one side of the surrounding rock on both sides of the tunnel wall, and the signals collected by the two acquisition units x and d are obtained, wherein the signals collected by the two acquisition units x and d both include: seismic signals mixed with noise, and the seismic signals include direct wave signals and reflected wave signals.

[0062] The schematic diagram of the HSP acoustic reflection method based on TBM construction in the embodiment of the present invention is shown in FIG. Figure 2 The strong vibration generated by the cutterhead cutting the rock during the TBM excavation process is used as a passive excitation source. The spatial array seismic wave sensor arranged behind the tunnel and the reference sensor near the cutterhead simultaneously receive seismic wave signals and other coupled noises. Among them, two sensors at both ends of one side of the surrounding rock on both sides of the tunnel are selected as two acquisition units, and the reference sensor is used as the reference unit. The acquisition unit r1 on the right is closer to the reference unit than the acquisition unit r2 on the right.

[0063] The vibration signal generated by TBM cutting rocks is a direct wave when it propagates backward to the right-end collection unit, and a reflected wave when it propagates forward to the reflection surface and then backward to the right-end collection unit.

[0064] S2. Perform LMS minimum mean square error adaptive delay estimation on the signals collected by the two collection units x and d to obtain the delay between the signals collected by the two collection units.

[0065] The signals collected by the two acquisition units include seismic wave signals and noise signals. Specifically:

[0066] S21, establish an adaptive filter, discretize the signals collected by the two acquisition units into x(n) and d(n), where x(n) is the signal collected by the left acquisition unit, d(n) is the signal collected by the right acquisition unit, d(n) is the target signal, d(n) is delayed for p sampling periods and input into the desired input end of the adaptive filter, x(n) is used as the input of the adaptive filter, and the error signal is obtained:

[0067] e(n)=d(n)-y(n)

[0068]

[0069] Where e(n) is the error signal, n represents the number of signal points, y(n) is the output signal of the adaptive filter with x(n) as the input signal, ω m (n) is the weight coefficient of the m-th order filter of the adaptive filter, and M is the order of the adaptive filter;

[0070] S22, initialize the weight coefficient ω of the delay estimator m (0), the error signal e(n) is expressed as:

[0071]

[0072] ω m (n+1)=ω m (n)+μe(n)x(nm)

[0073] Among them, μ is the convergence step size;

[0074] The mean square error performance function is expressed as:

[0075] V(W)=E[e 2 (n)] = E[(d(n)-W T (n)X(n)) 2 ]

[0076] Where V(W) is the mean square error performance function, E[] represents the expectation, W(n)=[ω1(n),ω2(n),...,ω m (n)] T ,X(n)=[x(n-1),x(n-2),...,x(nm)] T ;

[0077] According to the LMS criterion, when ω a (n) = R x -1 When P, the mean square error reaches the minimum value, the filter converges, and the weight coefficient vector ω is read. m (n) change curve.

[0078] P = E[d(n)X(n)] is the cross-correlation function of the desired output and input vectors, R x =E[X(n)X T (n)] is the autocorrelation function of the input signal, and the mean square error is the quadratic form of the weight vector. When the filter converges, the mean square error V(W) tends to zero, and the weight coefficient reaches its maximum value. At this time, ω a (n) = R x -1 P,a∈[1,2,...,m].

[0079] Theoretically, the delay is the a value corresponding to the maximum value of the weight coefficient minus the delay p, but the weight coefficient curve is composed of discrete points. The sinc interpolation of the weight coefficient is to obtain a more accurate a value in the discrete points.

[0080] S23, the above weight coefficient vector ω m The variation curve of (n) is the variation curve of the discrete signal. In order to obtain a more accurate maximum value of the weight coefficient vector, the discrete weight coefficient vector signal is resampled using the sinc interpolation algorithm to obtain more discrete points.

[0081] The basic idea of ​​sinc interpolation is to use the shape of the sinc function to approximate the value of the unknown point. For each interpolation point, calculate its corresponding sinc function and use it as a weight. Generally, the farther the point is from the interpolation point, the smaller the weight.

[0082] Perform sinc interpolation on the weight coefficients to obtain the time delay:

[0083]

[0084]

[0085] Among them, τ is the time delay and n is the number of signal points.

[0086] Specifically, the sampling rate is set to 500 Hz, the FIR filter order is set to 32, and the collection unit spacing is 200m for delay estimation. The simulation results are as follows: Figure 3 As shown in the figure, it can be seen that the delay estimation value is 0.088s, which is close to the theoretical delay of 0.064s (the iteration time corresponding to the maximum value of the normalized weight coefficient is 0.192s minus the total FIR iteration time of 0.128s). Table 1 shows the error relationship between the adaptive delay estimation value and the simulation value under different spacings.

[0087] Table 1

[0088]

[0089] S3. Based on the time delay between the signals collected by the two acquisition units, the time delay differential method is used to offset the co-frequency interference noise that is different from the excavation direction, and the seismic signal after the co-frequency noise is removed from the signals collected by the two acquisition units is obtained. The time domain delay differential algorithm delays the seismic signal to obtain a time domain delayed signal. This time domain delayed signal has a certain difference in phase with the reference signal, so destructive interference or additive interference will occur when it is merged, thereby reducing or eliminating the co-frequency noise in the excavation signal and achieving noise suppression.

[0090] Specifically:

[0091] S31. Establish a time differential model, taking the signal x(n) as the signal reference point, x(n) and d(n) can be expressed in the time domain as follows:

[0092] x(t)=s(t)+n(t)

[0093] d(t)=s(t-Δt)+n(t+Δt)

[0094] in,

[0095] Among them, x(t) and d(t) are the time domain representations of the signals collected by the two acquisition units, s(t) represents the time difference between the signals collected by the two acquisition units, Δl is the distance between two acquisition units, c is the wave speed;

[0096] S32, subtract the signal x(t) from the signal d(t) after time delay processing to obtain the delayed differential data x0(t) of the signals collected by the two collection units:

[0097] x0(t)=x(t)-d(t-τ)

[0098] x0(t-τ)=s(t-Δt+τ)+n(t+Δt-τ)

[0099] Let τ = Δt, x0(t) = s(t) - s(t-2Δt), and solve the difference equation to obtain the denoised seismic signals s(t) and s(t-Δt). s(t) and s(t-Δt) are the seismic signals after removing the same-frequency noise from the signals x(n) and d(n).

[0100] There is a path difference between the direct wave and the reflected wave received at the same time. If the distance between the two sensors is much smaller than the path of the seismic wave, the attenuation of the seismic wave signal between the two sensors can be ignored, and the seismic signal is expressed as:

[0101] s(t)=p(t)+q(t-Δt)

[0102] Among them, p(t) is the reflected wave, and q(t-Δt) is the direct wave;

[0103] The time-delay difference equation of the reflected wave is:

[0104] s(t)-s(t+Δt)=p(t)+q(t-Δt)-p(t+Δt)-q(t)

[0105] =p(t)-p(t+Δt)+q(t-Δt)-q(t)

[0106] Δt is the time difference between the direct wave and the reflected wave. s(t) and s(t+Δt) are aligned in time. The phase of the direct wave in the seismic signal remains roughly the same, while the phase of the reflected wave is different from the phase in the direct wave due to the delay. The direct wave usually arrives before the reflected wave, and in simple cases, the path of the direct wave is shorter and the shape is simpler. The reflected wave undergoes multiple reflections and refractions before reaching the receiver, which may cause its shape and phase to be different from the direct wave. The purpose of time delaying and differential processing of the signal is actually to try to use these differences between the direct wave and the reflected wave to reduce or eliminate the influence of the reflected wave. The reflected wave is almost eliminated in the differential equation, which mainly contains the direct wave signal. Solving the differential equation can obtain the direct wave signal q(t-Δt); then substitute q(t-Δt) into the de-noised seismic signal s(t) to obtain the separated reflected wave.

[0107] Assume the sensor sampling rate is f s , then the above formula can be expressed as:

[0108] p(n)-p(n+m)=s(n)-s(n+m)

[0109] m=Δtf s

[0110] The separated reflection wave seismic record data p(n) can be obtained by solving the differential equation.

[0111] The tunneling data calculated by the time-domain delay difference algorithm is compared with the original data. Figure 4 As shown, the correlation coefficient between the solved data and the excavation data without noise is 0.9992, which proves that the time domain delay difference algorithm can effectively suppress the noise interference from different sources in the excavation data.

[0112] S4. Segment the seismic signals after removing the same-frequency noise from the signals collected by the two acquisition units in S3, cross-correlate the two segmented seismic signals after removing the same-frequency noise to obtain the reflected wave frequency domain signal, vertically superimpose the segmented signals, and calculate the propagation time of the reflected wave to obtain the distance between the reference unit and the unfavorable geological body in front of the tunnel.

[0113] The improvement of the signal-to-noise ratio after cross-correlation interference is related to the acquisition time and bandwidth. According to Landrum's analysis of the vibroseis scanning signal, the signal-to-noise ratio formula after cross-correlation interference can be expressed as:

[0114]

[0115] Where ω represents the bandwidth; T represents the acquisition time.

[0116] Since the vibration of TBM excavation is uncontrollable and the bandwidth cannot be changed, the signal-to-noise ratio can only be improved by extending the acquisition time. However, the recording time T of HSP seismic wave detection data is limited by the performance of the acquisition instrument, the on-site excavation conditions and the data processing method, and unlimited data acquisition cannot be achieved. In practical applications, the acquisition time is increased by cross-correlating the data collected multiple times and then superimposing them in sections, and the signals collected multiple times are considered to be cut into the same time length according to the actual detection needs.

[0117] When processing a tunneling signal with a signal-to-noise ratio of 0 dB, the curve of the signal-to-noise ratio of the data changes with the number of stacking times is as follows: Figure 5 shown.

[0118] The reflected waves separated in S3 are segmented, cross-correlated with the signals collected by the reference unit, and then vertically superimposed to obtain the distance between the reference unit and the unfavorable geological body in front of the tunnel. Specifically:

[0119] By cross-correlating the signals collected from the two observation points, it is equivalent to a seismic wave field record with one point as the excitation source and the other point as the virtual source. This method can construct a seismic response Green's function composed of two different observation points. The cross-correlation processing of the two segmented seismic signals after removing the same-frequency noise is specifically as follows:

[0120] The Green function relationship is constructed to represent the cross-correlation function between two acquisition units:

[0121]

[0122] Among them, C(s x ,s d ,t) is the signal s x With signal s d The cross-correlation coefficient in the time domain, s x is the segmented seismic signal collected by acquisition unit x after removing the same-frequency noise, s d is the segmented seismic signal collected by acquisition unit d after removing the same-frequency noise, G(s x ,s d ,t) is the Green function, which represents the displacement response of the unit force in the time domain t generated by the acquisition units x and d, G(s x ,s d ,-t) is the Green's function, which represents the displacement response of the unit force in the time domain -t produced by the acquisition unit x and d.

[0123] Assuming that the detected noise signal sources are uncorrelated, C(s x ,s d ,t) is expressed in the frequency domain as:

[0124]

[0125] In the formula, C(s x ,s d ,ω) represents C(s x ,s d ,t) in the frequency domain, u1(s x ,ω) represents the seismic frequency domain signal received by the acquisition unit x; u2(s d ,ω) represents the seismic frequency domain signal received by the acquisition unit d; G(s x ,s d ,ω) represents G(s x ,s d ,t) is expressed in the frequency domain; * indicates the complex number is conjugated, || indicates the complex number is real part; P(ω) represents the expression of the reflected wave p(t) in the frequency domain.

[0126]

[0127] Wherein, T is the collection time of the two collection units;

[0128] Vertically superimpose each signal of P(ω) to obtain the signal:

[0129]

[0130]

[0131] Where N is the total number of segments; floor is the number of rounded down when the data does not meet the segment length; T N Indicates the segment length; P n (ω,n) is the frequency domain signal of the reflected wave obtained by superimposing n segments P(ω). The propagation time of the reflected wave is calculated to obtain the distance between the reference unit and the unfavorable geological body in front of the tunnel.

[0132] After delaying the tunneling data, the reflected wave restored by solving the differential equation was compared with the simulated reflected wave. The correlation coefficient between the restored reflected wave and the simulation result was 0.9929, and the signal-to-noise ratio of the solved reflected wave was 41.47dB. This proves that the time-domain delay differential algorithm has a good effect in distinguishing reflected waves from direct waves from tunneling data.

[0133] In an exemplary embodiment, a computer-readable storage medium is included, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the above-mentioned cross-correlation seismic interferometry HSP method is implemented.

[0134] See also Figure 6 In an exemplary embodiment, an electronic device is also included, including at least one processor, at least one memory, and at least one communication bus.

[0135] Wherein, a computer program is stored in the memory, and the computer program includes computer-readable instructions. The processor calls the computer-readable instructions stored in the memory through the communication bus to execute the above-mentioned cross-correlation seismic interferometry HSP method.

[0136] The above description of the disclosed embodiments enables one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A cross-correlation seismic interferometry HSP method, characterized in that: The following steps are involved: S1. Based on the HSP acoustic wave reflection method of TBM construction, two acquisition units x and d are arranged on one side of the surrounding rock on both sides of the tunnel wall, and the signals collected by the two acquisition units x and d are obtained, wherein the signals collected by the two acquisition units x and d both include: seismic signals mixed with noise; S2, performing LMS minimum mean square error adaptive delay estimation on the signals collected by the two collection units x and d to obtain the delay between the signals collected by the two collection units; Specifically: S21, establish an adaptive filter, the two acquisition units x and d acquire signals x(n) and d(n), x(n) is used as the input signal of the adaptive filter, d(n) is used as the target signal, and an error signal is obtained: e(n)=d(n)-y(n) Wherein, e(n) is the error signal, n represents the number of signal points, y(n) is the output signal of the adaptive filter with x(n) as the input signal, ω m (n) is the weight coefficient of the m-th order filter, and M is the order of the adaptive filter; S22, the error signal e(n) is expressed as: ω m (n+1)=ω m (n)+μe(n)x(n-m) Among them, μ is the convergence step size; The mean square error performance function is expressed as: V(W)=E[e 2 (n)]=E[(d(n)-W T (n)X(n)) 2 ] Where V(W) is the mean square error performance function, E[] represents the expectation, W(n)=[ω1(n),ω2(n),...,ω m (n)] T ,X(n)=[x(n-1),x(n-2),...,x(nm)] T ; P = E[d(n)X(n)] is the cross-correlation function of the desired output and input vectors, R x =E[X(n)X T (n)] is the autocorrelation function of the input signal, and the mean square error is the quadratic form of the weight vector. When the filter converges, the mean square error V(W) tends to zero, and the weight coefficient reaches its maximum value. At this time, ω a (n) = R x -1 P,a∈[1,2,...,m]; S23, perform sinc interpolation on the weight coefficient to obtain the time delay: Among them, τ is the time delay between the two acquisition units collecting signals, and n represents the number of signal points; S3, based on the time delay between the collected signals of the two acquisition units, the signals collected by the two acquisition units are used to offset the same-frequency interference noise that is different from the excavation direction using a time delay differential method, so as to obtain the seismic signal after removing the same-frequency noise from the collected signals of the two acquisition units; Specifically: S31. Establish a time difference model. x(n) and d(n) can be expressed in the time domain as follows: x(t)=s(t)+n(t) d(t)=s(t-Δt)+n(t+Δt) Among them, x(t) and d(t) are the time domain representations of the signals collected by the two acquisition units, s(t) represents the time domain representation of the seismic signal, n(t) represents the time domain representation of the noise signal, and Δt represents the time difference between the signals collected by the two acquisition units; S32, subtract the signal x(t) from the signal d(t) after time delay processing to obtain the delayed differential data x0(t) of the signals collected by the two collection units: x0(t)=x(t)-d(t-τ) x0(t-τ)=s(t-Δt+τ)+n(t+Δt-τ) Let τ = Δt, x0(t) = s(t) - s(t-2Δt), and solve the difference equation to get the de-noised seismic signals s(t) and s(t-Δt). s(t) and s(t-Δt) are the seismic signals after removing the same-frequency noise from the signals x(n) and d(n); S4. Segment the seismic signals after removing the same-frequency noise from the signals collected by the two acquisition units in S3, cross-correlate the two segmented seismic signals after removing the same-frequency noise to obtain the reflected wave frequency domain signal, vertically superimpose the segmented signals, and calculate the propagation time of the reflected wave to obtain the distance between the reference unit and the unfavorable geological body in front of the tunnel.

2. A cross-correlation seismic interferometry HSP method according to claim 1, characterized in that: The seismic signal includes a direct wave signal and a reflected wave signal.

3. A cross-correlation seismic interferometry HSP method according to claim 1, characterized in that: The cross-correlation processing of the two segmented seismic signals after removing the same-frequency noise is specifically as follows: Among them, C(s x ,s d ,t) is the signal s x With signal s d The cross-correlation coefficient in the time domain, s x is the segmented seismic signal collected by acquisition unit x after removing the same-frequency noise, s d is the segmented seismic signal collected by acquisition unit d after removing the same-frequency noise, G(s x ,s d ,t) is the Green function, which represents the displacement response of the unit force in the time domain t generated by the acquisition units x and d, G(s x ,s d ,-t) is the Green's function, which represents the displacement response of the unit force in the time domain -t produced by the acquisition unit x and d.

4. A cross-correlation seismic interferometry HSP method according to claim 3, characterized in that: C(s x ,s d ,t) is expressed in the frequency domain as: In the formula, C(s x ,s d ,ω) represents C(s x ,s d ,t) in the frequency domain, u1(s x ,ω) represents the seismic frequency domain signal received by the acquisition unit x; u2(s d ,ω) represents the seismic frequency domain signal received by the acquisition unit d; G(s x ,s d ,ω) represents G(s x ,s d ,t) is expressed in the frequency domain; * indicates the complex number is conjugated, || indicates the complex number is real part; P(ω) indicates the expression of the reflected wave p(t) in the frequency domain; Wherein, T is the collection time of the two collection units; Vertically superimpose each signal of P(ω) to obtain the signal: Where N is the total number of segments; floor is the number of rounded down when the data does not meet the segment length; T N Indicates the segment length; P n (ω,n) is the frequency domain signal of the reflected wave obtained by superimposing n segments P(ω). The propagation time of the reflected wave is calculated to obtain the distance between the reference unit and the unfavorable geological body in front of the tunnel.

5. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 4 are implemented.

6. An electronic device, characterized in that: The method comprises a processor and a memory, wherein the processor and the memory are interconnected, wherein the memory is used to store a computer program, the computer program comprises computer-readable instructions, and the processor is configured to call the computer-readable instructions to execute the method according to any one of claims 1 to 4.

Images