A tunnel advance detection method based on point spread function correction

Abstract

This invention discloses a tunnel advance detection method based on point spread function correction, comprising: setting up receiver points and shot points, and acquiring observed seismic data; constructing a rectangular grid tunnel geological model; given a velocity model, performing forward modeling to obtain the forward propagation wavefield; using the seismic data to propagate backward in time to obtain the backward propagation wavefield from the receiver points, and obtaining the reverse time migration imaging result through cross-correlation imaging conditions; using the Poynting vector method to obtain the propagation angles of the forward propagation wavefield from the shot point and the backward propagation wavefield from the receiver point; decomposing the wavefield into beam superposition based on the propagation angles, and then obtaining the local illumination matrix; performing coordinate transformation on the local illumination matrix to obtain the angular domain point spread function; and using the point spread function to correct the reverse time migration imaging result to obtain the corrected imaging result m. psf This invention can obtain a more accurate angular domain point spread function and can be used to correct tunnel advance detection imaging.

Description

Technical Field

[0001] This invention belongs to the field of tunnel advanced geological prediction technology, specifically relating to a tunnel advanced detection method based on point spread function correction. Background Technology

[0002] During tunnel construction, complex geological conditions and adverse geological formations can easily trigger disasters such as water inrush, mudslides, and collapses. These can delay construction progress or cause casualties and serious economic losses. To ensure construction safety, it is essential to accurately identify the geological conditions ahead of the tunnel face using advanced detection technologies before excavation. Currently, various geophysical advanced geological detection methods are applied in tunnel construction, enabling rapid and non-destructive geological exploration ahead of the tunnel face. These include seismic wave methods, resistivity methods, transient electromagnetic methods, and ground-penetrating radar methods.

[0003] Among existing tunnel advance detection methods, seismic wave methods have become the mainstream technology due to their non-destructive testing, large-scale forward modeling capabilities, and excellent interface identification. These include tunnel seismic prediction, true reflection tomography, and tunnel seismic monitoring during drilling. Among these, the reverse time migration imaging method based on the two-way wave equation is particularly effective, enabling high-precision imaging of faults, interfaces, and other adverse geological structures ahead of the tunnel face, providing crucial information for construction decisions. However, when using seismic wave reverse time migration imaging for tunnel advance detection, the determination of the location of anomalies in complex geological conditions ahead of the tunnel still requires further study. According to inversion theory, conventional migration imaging results can be viewed as the convolution of the true reflectivity model and a fuzzy operator, also known as the Hessian operator. Theoretically, using the inverse of the Hessian operator can effectively suppress the fuzziness effect in migration imaging results, but the computation and storage of the complete Hessian operator present significant challenges, making the introduction of its approximation essential. The point spread function can be considered a high-precision approximation of the Hessian operator, characterizing the illumination characteristics of the observation system at a specific spatial point. In the field of seismic exploration, the inverse correction migration imaging results using the point spread function operator can effectively eliminate image blur, improve imaging resolution, and balance image amplitude.

[0004] The inventors of this application discovered through research that point spread function correction plays a significant role in seismic exploration, but it has not yet been applied to tunnel advance detection. Current tunnel advance detection methods do not fully utilize the illumination characteristics and angular domain information of the observation system, resulting in insufficient imaging accuracy of anomalies ahead of the tunnel. To obtain more accurate tunnel advance detection results, point spread function correction can be introduced into tunnel advance detection. It should be noted that directly using traditional Fourier transform to calculate the point spread function requires substantial computational resources; therefore, a more efficient angular domain method is used to calculate the point spread function. Summary of the Invention

[0005] The purpose of this invention is to provide a tunnel advance detection method based on point spread function correction. This method uses the Poynting vector method to calculate the wave field propagation angle, constructs a local illumination matrix, and then converts it into an angle-domain point spread function, which is used to correct the tunnel reverse-time migration imaging results. The angle-domain point spread function contains information about the observation system inside the tunnel and information between the tilt angle of the anomalous body and the incident angle and scattering angle of the wave field near the anomalous body, thus obtaining more accurate tunnel migration imaging results.

[0006] The technical solution adopted in this invention is a tunnel advance detection method based on point spread function correction, comprising:

[0007] Step 1: Set up shot points and receiver points, and acquire observed seismic data. ;in Represents a time variable; These indicate the locations of the geophone point and the shot point, respectively.

[0008] Step 2: Construct a rectangular grid tunnel geological model for numerical simulation calculations, and set the number of horizontal grid points in the rectangular grid tunnel geological model. and vertical grid point count And set the spatial sampling interval for forward modeling, including the lateral sampling interval. and longitudinal sampling interval Time sampling interval Maximum sampling time ;

[0009] Step 3: Given a velocity model describing the propagation velocity distribution of waves in the underground medium. The velocity model c is determined based on the lateral sampling interval. and longitudinal sampling interval Discretize the signal onto the grid defined in step 2, and use the source signal to numerically solve the wave equation. Using a time sampling interval dt, perform forward propagation within the time range [0, Nt] to obtain the forward modeling simulation from the shot point. Go to a point in space forward wave field ;

[0010] Step 4: Transfer the observed seismic data collected in Step 1 Backpropagation in time generates a backpropagating wavefield from the receiver. The forward transmission of the Tron field and reverse propagation wave field The reverse time migration imaging results were obtained by calculating the cross-correlation imaging conditions. , as input for subsequent correction;

[0011] Step 5: Use the Poynting vector method to obtain the propagating wave field from the shot point. Communication perspective Backpropagating wave field from the detector point The angle of dissemination ;

[0012] Step 6: Based on the propagation angle of the wave field and The wave field can be decomposed into a superposition of beams, and then the local illumination matrix can be obtained. ;

[0013] Step 7: Configure the local illumination matrix Perform coordinate transformation from the incident / scattering angle of the wave field. Transform to the tilt wavenumber domain ,in, The wavenumber corresponding to the dip angle of the anomalous body. Indicates the tilt angle, the specific correspondence is as follows: Figure 2 As shown; then, the point spread function in the angle domain is obtained. ;

[0014] Step 8: Apply the reverse time-shifted image obtained in Step 4. Transform to the wavenumber domain using the point spread function. Reverse time migration imaging results After correction, the result is converted back to the spatial domain to obtain the final corrected imaging result m. psf .

[0015] The invention is further characterized in that,

[0016] In step 3, wave field The following sound wave equation is satisfied:

[0017] (1)

[0018] In equation (1), For spatial location The velocity of the medium at that location; It is the density of the medium. This indicates the speed at which seismic waves propagate in a medium.

[0019] To simplify the formula, omissions are made. ,use express.

[0020] In step 4, back-transfer seismic data is used. The reverse propagation wavefield from the receiver is obtained, and the reverse time migration imaging result is obtained by calculating the zero-delay cross-correlation between the forward and reverse propagation wavefields at each time step. The definition is as follows:

[0021] (2)

[0022] In formula (2) Indicates the firing point inside the tunnel arrive The forward propagation wave field, Indicates from to the detector point inside the tunnel The reverse propagation wave field.

[0023] To remove low-frequency noise, in equation (2) Perform Laplace filtering.

[0024] In step 5, the propagation angle of the forward wave field from the shot point and the propagation angle of the back propagating wave field from the detector point The calculation is performed using the Poynting vector method, as follows:

[0025] (3)

[0026] In equation (3), and It is the Poynting vector of the propagating wave field. The x and z components; and It is the Poynting vector of the anti-propagating wave field. The x and z components, and The definition is as follows:

[0027] (4)

[0028] In equation (4), This represents the spatial derivative of the wave field.

[0029] In step 6, the local illumination matrix is ​​obtained. The specific methods are as follows:

[0030] The propagation angle of the wave field has been obtained from equations (3) and (4). Given a broadband source, such as a Ricker wavelet, the wave field in equation (2) and This can be represented as beam superposition:

[0031] (5)

[0032] In the formula It describes the location of the shot point inside the tunnel. spread to Green's function; It describes the wave from propagation to the receiving point inside the tunnel Green's function; It's frequency. This is the source spectrum. The mean square amplitude within a short time interval T is calculated using equation (5):

[0033] (6)

[0034] During wave field propagation, given a constant dominant frequency, we can obtain:

[0035] (7)

[0036] Based on equation (7), the local illumination matrix can be obtained:

[0037] (8)

[0038] In equation (8), This represents the local wave number.

[0039] In step 7, the angle domain point spread function The definition is as follows:

[0040] (9)

[0041] In equation (9), k g and k i k represents the scattered and incident wavenumbers, respectively. d =k i +k g It is the wavenumber related to the construction dip angle;

[0042] Perform a coordinate transformation on equation (8) to change the incident / scattering angle. Transform to wavenumber domain

[0043] (10)

[0044] In equation (10), For Jacobian matrices, It is the tilt angle. ;

[0045] First, start with the incident / scattering angle Switch to tilt / reflection angle Their relationship is , Then switch to ,in ;this and The relationship between them is ;get:

[0046] (11)

[0047] Substituting equation (11) into equation (9) yields the angular domain point spread function:

[0048] (12)

[0049] In step 8, the tunnel reverse time migration imaging result is regarded as a convolution between the real tunnel velocity model and the point spread function, which is represented as a multiplication in the wavenumber domain. The correction result in the wavenumber domain is obtained using the inverse of the point spread function.

[0050] (13)

[0051] In equation (14), M mig Reverse time migration imaging results Fourier transform;

[0052] By analyzing M psf Perform an inverse Fourier transform to obtain the corrected imaging result m. psf .

[0053] This invention obtains conventional tunnel reverse-time migration imaging results; uses the Poynting vector method to directly decompose the wavefield from the shot and receiver points into angular domain components in the time domain; derives the point spread function by calculating the local illumination matrix; and uses the obtained point spread function to correct the reverse-time migration imaging results, yielding more accurate results. Compared with commonly used tunnel seismic exploration methods, this method can significantly improve the imaging accuracy of anomalies ahead. The reasons for these advantages are: firstly, the use of point spread function correction can suppress the blurring effect in traditional imaging; secondly, the angular domain point spread function contains information about the dip angle of the anomaly ahead of the excavation face, the incident angle and scattering angle of the wavefield near the anomaly, and information from the observation system inside the tunnel, thus obtaining more accurate results. Attached Figure Description

[0054] Figure 1 This is a flowchart of a tunnel advance detection method based on point spread function correction according to an embodiment of the present invention;

[0055] Figure 2 This is a schematic diagram of the local coordinate system of the present invention;

[0056] Figure 3 This is the actual tunnel velocity model in the model example of this invention;

[0057] Figure 4 This is a schematic diagram of the tunnel observation system in the model example of this invention;

[0058] Figure 5 The imaging results are obtained using the traditional reverse time migration method in the model examples of this invention;

[0059] Figure 6 The imaging results are obtained by correcting the method proposed in this invention using the model examples of this invention; Detailed Implementation

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

[0061] like Figure 1 As shown, this embodiment of the invention provides a tunnel advance detection method based on point spread function correction, which is implemented according to the following steps:

[0062] Step 1: Set up shot points and receiver points, and acquire observed seismic data. ;in Represents a time variable; These indicate the locations of the geophone point and the shot point, respectively.

[0063] Step 2: Construct a rectangular mesh tunnel geological model for numerical simulation calculations. This model is used to discretize the velocity model and serves as the computational grid framework for this imaging algorithm. The number of horizontal grid points in the rectangular mesh tunnel geological model is set. and vertical grid point count And set the spatial sampling interval for forward modeling, including the lateral sampling interval. and longitudinal sampling interval Time sampling interval Maximum sampling time ;

[0064] Step 3: Given a velocity model describing the velocity distribution of wave propagation in the subsurface medium. The velocity model c is determined based on the lateral sampling interval. and longitudinal sampling interval Discretize the signal onto the grid defined in step 2, and use the source signal to numerically solve the wave equation. Using a time sampling interval dt, perform forward propagation within the time range [0, Nt] to obtain the forward modeling simulation from the shot point. Go to a point in space forward wave field ;

[0065] In step 3, wave field The following sound wave equation is satisfied:

[0066] (1)

[0067] In equation (1), For spatial location The velocity of the medium at that location; It is the density of the medium. This indicates the speed at which seismic waves propagate in a medium.

[0068] To simplify the formula, omissions are made. ,use express.

[0069] Step 4: Transfer the observed seismic data collected in Step 1 Backpropagation in time generates a backpropagating wavefield from the receiver. The forward transmission of the Tron field and reverse propagation wave field The reverse time migration imaging results were obtained by calculating the cross-correlation imaging conditions. , as input for subsequent correction;

[0070] In step 4, the reverse time migration imaging results It can be obtained by calculating the zero-delay cross-correlation between the source-side wavefield and the receiver-side wavefield at each time step, as defined below:

[0071] (2)

[0072] In formula (2) Indicates the source x inside the tunnel i The propagating wave field to x, Indicates the distance from x to the receiving side x inside the tunnel. g The reverse propagation wave field;

[0073] To remove low-frequency noise, in equation (2) Perform Laplace filtering;

[0074] Step 5: Use the Poynting vector method to obtain the propagating wave field from the shot point. Communication perspective Backpropagating wave field from the detector point The angle of dissemination ;

[0075] In step 5, the propagation angle of the forward wave field from the shot point and the propagation angle of the back propagating wave field from the detector point The calculation is performed using the Poynting vector method, as follows:

[0076] (3)

[0077] In equation (3), and It is the Poynting vector of the propagating wave field. The x and z components; and It is the Poynting vector of the anti-propagating wave field. The x and z components, and The definition is as follows:

[0078] (4)

[0079] In equation (4), Represents the spatial derivative of the wave field;

[0080] Step 6: Based on the propagation angle of the wave field and The wave field can be decomposed into a superposition of beams, and then the local illumination matrix can be obtained. ;

[0081] In step 6, the local illumination matrix is ​​obtained. The specific methods are as follows:

[0082] The propagation angle of the wave field has been obtained from equations (3) and (4). Given a broadband source, such as a Ricker wavelet, the wave field in equation (2) and This can be represented as beam superposition:

[0083] (5)

[0084] In the formula It describes the location of the shot point inside the tunnel. spread to Green's function; It describes the wave from propagation to the receiving point inside the tunnel Green's function; It's frequency. This is the source spectrum. The mean square amplitude within a short time interval T is calculated using equation (6):

[0085] (6)

[0086] During wave field propagation, given a constant dominant frequency, we can obtain:

[0087] (7)

[0088] Based on equation (7), the local illumination matrix can be obtained:

[0089] (8)

[0090] In equation (8), Indicates the local wavenumber;

[0091] Step 7: Configure the local illumination matrix Perform coordinate transformation from the incident and scattering angles of the wave field. Transform to the tilt wavenumber domain ,in, The wavenumber corresponding to the dip angle of the anomalous body. Indicates the tilt angle, the specific correspondence is as follows: Figure 2 As shown; then, the point spread function in the angle domain is obtained. ;

[0092] In step 7, the angle domain point spread function The definition is as follows:

[0093] (9)

[0094] In equation (9), k g and k i k represents the scattered and incident wavenumbers, respectively. d =k i +k g It is the wavenumber related to the construction dip angle;

[0095] Perform a coordinate transformation on equation (8) to change the incident / scattering angle. Transform to wavenumber domain :

[0096] (10)

[0097] In equation (11), For Jacobian matrices, It is the tilt angle. ;

[0098] First, start with the incident / scattering angle Switch to tilt / reflection angle Their relationship is , Then switch to ,in ;this and The relationship between them is ;get:

[0099] (11)

[0100] Substituting equation (11) into equation (9) yields the angular domain point spread function:

[0101] (12)

[0102] Step 8: Apply the reverse time-shifted image obtained in Step 4. Transform to the wavenumber domain using the point spread function. Reverse time migration imaging results After correction, the result is converted back to the spatial domain to obtain the final corrected imaging result m. psf .

[0103] In step 8, the tunnel reverse time migration imaging result is regarded as a convolution between the real tunnel velocity model and the point spread function, which is represented as a multiplication in the wavenumber domain. The correction result in the wavenumber domain is obtained using the inverse of the point spread function.

[0104] (13)

[0105] In equation (13), M mig Reverse time migration imaging results Fourier transform;

[0106] By analyzing M psf Perform an inverse Fourier transform to obtain the corrected imaging result m. psf .

[0107] Model example:

[0108] The specific implementation process of this invention is applied to a tunnel model. The model is 200 meters wide laterally and 150 meters long longitudinally. The model includes low-velocity anomalies and different media. The tunnel body is 80 m long and 10 m high. The velocity model is as follows... Figure 3 As shown.

[0109] The layout of the observation system is as follows Figure 4 As shown, on each side of the tunnel, there are 35 sources spaced 2 m apart and 7 receivers spaced 10 m apart. On the tunnel surface, there are 11 sources spaced 1 m apart and 6 receivers spaced 2 m apart.

[0110] Figure 5 The results of conventional reverse time migration imaging of the tunnel model show that artifacts exist in the first interface layer, which interfere with the judgment of seismic information in front of the tunnel. In the second medium layer, the interface behind the cave is not clearly visible in the areas indicated by the two white arrows.

[0111] Figure 6The image is obtained after point spread function correction in steps 1 to 8 of the method proposed in this invention. For the first interface, artifacts are shown to be suppressed; for the cave, the imaging of the front and rear interfaces is more balanced and the location is clearer. It can be seen that the method proposed in this invention can effectively image anomalies in front of the tunnel and improve image resolution.

[0112] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. A tunnel advance detection method based on point spread function correction, characterized in that, Includes the following steps: Step 1: Set up receiver points and shot points, and acquire observed seismic data. ;in Represents a time variable; These indicate the locations of the geophone point and the shot point, respectively. Step 2: Construct a rectangular mesh tunnel geological model for numerical simulation calculations. Set the number of horizontal grid points in the rectangular mesh tunnel geological model. and vertical grid point count And set the spatial sampling interval for forward modeling, including the lateral sampling interval. and longitudinal sampling interval Time sampling interval Maximum sampling time ; Step 3: Given a velocity model describing the velocity distribution of wave propagation in the subsurface medium. , velocity model According to the horizontal sampling interval and longitudinal sampling interval Discretize the data onto the rectangular grid tunnel geological model constructed in step 2. Using the seismic source signal, solve the wave equation numerically. Employ a time sampling interval dt and perform forward propagation within the time range [0, Nt] to obtain the forward modeling results from the shot point. Go to a point in space forward wave field ; Step 4: Transfer the observed seismic data collected in Step 1 Backpropagation is performed over a time span to generate a backpropagating wavefield from the receiver. The forward transmission of the Tron field and reverse propagation wave field The reverse time migration imaging results were obtained by calculating the cross-correlation imaging conditions. ; Step 5: Use the Poynting vector method to obtain the propagating wave field from the shot point. Communication perspective and the backpropagating wave field from the detector point The angle of dissemination ; Step 6: Based on the propagation angle of the wave field and The wave field is decomposed into a superposition of beams, and then the local illumination matrix is ​​obtained. ; Step 7: Configure the local illumination matrix Perform coordinate transformation from the incident / scattering angle of the wave field. Transform to the tilt wavenumber domain ,in, The wavenumber corresponding to the dip angle of the anomalous body. To represent the tilt angle, we then obtain the angular domain point spread function. ; Step 8: Apply the reverse time-shifted image obtained in Step 4. Transform to the wavenumber domain using the point spread function. For reverse time-shifted imaging After correction, the result is converted back to the spatial domain to obtain the final corrected imaging result m. psf .

2. The tunnel advance detection method based on point spread function correction according to claim 1, characterized in that, In step 3, the forward propagation wave field The following sound wave equation is satisfied: (1)； In equation (1), Spatial location The velocity of the medium at that location; It is the density of the medium. This indicates the speed at which seismic waves propagate in a medium; in the formula express .

3. The tunnel advance detection method based on point spread function correction according to claim 2, characterized in that, In step 4, the reverse time migration imaging results are obtained by calculating the zero-delay cross-correlation between the forward and reverse propagation wavefields at each time step. The definition is as follows: (2)； In formula (2) Indicates the firing point inside the tunnel arrive The forward propagation wave field, Indicates from to the detector point inside the tunnel The reverse propagation wave field; To remove low-frequency noise, in equation (2) Perform Laplace filtering.

4. The tunnel advance detection method based on point spread function correction according to claim 3, characterized in that, In step 5, the Poynting vector method is used to obtain the propagating wave field from the shot point. Communication perspective and the backpropagating wave field from the detector point The angle of dissemination Specifically, it includes: (3)； In equation (3), and It is the Poynting vector of the propagating wave field. The x and z components; and It is the Poynting vector of the anti-propagating wave field. The x and z components, and The definition is as follows: (4)； In equation (4), This represents the spatial derivative of the wave field.

5. The tunnel advance detection method based on point spread function correction according to claim 4, characterized in that, In step 6, the local illumination matrix is ​​obtained. The specific methods are as follows: The propagation angle of the wave field has been obtained from equations (4) and (5). Given a broadband source, the wave field in equation (2) and Represented as beam superposition: (5)； In the formula It describes the location of the shot point inside the tunnel from the wave. spread to Green's function; It describes the wave from propagation to the receiving point inside the tunnel Green's function; It's frequency. It is the source spectrum; calculate the mean square amplitude within a short time interval T using equation (5): (6)； During wave field propagation, given a constant dominant frequency, we obtain: (7)； Based on equation (7), the local illumination matrix is ​​obtained: (8)； In equation (8), This represents the local wave number.

6. The tunnel advance detection method based on point spread function correction according to claim 5, characterized in that, In step 7, the angle domain point spread function The definition is as follows: (9)； In equation (9), k g and k i k represents the scattered and incident wavenumbers, respectively. d =k i +k g It is the wavenumber related to the construction dip angle; Perform a coordinate transformation on equation (8) to change the incident / scattering angle. Transform to wavenumber domain : (10)； In equation (10), For Jacobian matrices, It is the tilt angle. ; First, start with the incident / scattering angle Switch to tilt / reflection angle Their relationship is , Then switch to ,in ; this and The relationship between them is ;get: (11)； Substituting equation (11) into equation (9) yields the angular domain point spread function: (12)。 7. The tunnel advance detection method based on point spread function correction according to claim 6, characterized in that, In step 8, the tunnel reverse time migration imaging result is regarded as a convolution between the real tunnel velocity model and the point spread function, which is represented as a multiplication in the wavenumber domain. The correction result in the wavenumber domain is obtained using the inverse of the point spread function. (13)； In equation (13), M mig Reverse time migration imaging results Fourier transform; By analyzing M psf Perform an inverse Fourier transform to obtain the corrected imaging result m. psf .

Images