Tunnel advanced detection method based on point spread function correction

Abstract

The invention discloses a tunnel advanced detection method based on point spread function correction, and the method comprises the steps: setting a detection point and a shot point, and collecting observation seismic data; constructing a rectangular grid tunnel geologic model; giving a speed model, and performing forward modeling to obtain a forward transmission wave field; seismic data is used for back propagation in time, a back propagation wave field from a detection point is obtained, and a reverse time migration imaging result is solved through a cross-correlation imaging condition; obtaining a forward propagation wave field propagation angle from a shot point and a backward propagation wave field propagation angle from a detection point by using a Poynting vector method; decomposing the wave field into superposition of wave beams according to the sum of the propagation angles of the wave field, and further solving a local illumination matrix; performing coordinate transformation on the local illumination matrix, and solving an angle domain point spread function; and correcting the reverse time migration imaging result by using a point spread function to obtain a corrected imaging result mpsf. According to the invention, a more accurate angle domain point spread function can be obtained and is used for correcting tunnel advanced detection imaging.

Description

Technical Field

[0001] The invention belongs to the technical field of tunnel advance geological prediction, and in particular relates to a tunnel advance detection method based on point spread function correction. Background Art

[0002] During tunnel construction, complex geological conditions and poor structural conditions can easily lead to disasters such as water inrush, mud gushing, and landslides, which can delay construction progress at best and cause casualties and significant economic losses at worst. To ensure construction safety, advanced detection technology must accurately identify the geological conditions ahead of the tunnel face before excavation. Currently, a variety of geophysical advanced geological detection methods are used in tunnel construction, enabling rapid and non-destructive geological exploration ahead of the tunnel face. These include seismic wave detection, resistivity, transient electromagnetic detection, and geological radar.

[0003] Among existing tunnel advance detection methods, seismic wave methods have become the mainstream due to their nondestructive testing, large-scale forward imaging capabilities, and excellent interface recognition. These include tunnel seismic prediction, true reflection tomography, and tunnel-while-drilling seismic. Among these, reverse-time migration (RTM) imaging, based on the bidirectional wave equation, is particularly effective, enabling high-precision imaging of unfavorable geological structures such as faults and interfaces ahead of the tunnel face, providing a crucial basis for construction decision-making. However, when using seismic RTM for tunnel advance detection, determining the location of anomalies in the complex geology ahead of the tunnel remains a challenge. According to inversion theory, conventional migration imaging results can be viewed as the convolution of a true reflectivity model with a fuzzy operator, also known as the Hessian operator. In theory, using the inverse of the Hessian operator can effectively suppress fuzziness in the migration imaging results. However, the computation and storage of the full Hessian operator present significant challenges, making its approximation essential. The point spread function (PSF) can be considered a high-precision approximation of the Hessian operator, characterizing the illumination characteristics of the observation system at a specific point in space. In the field of seismic exploration, the use of the point spread function operator to inversely correct the migration imaging results can effectively eliminate imaging blur, improve imaging resolution, and balance image amplitude.

[0004] The inventors of this application have discovered that while point spread function correction (PSF) is highly effective in seismic exploration, it has not yet been applied to tunnel advance detection. Current tunnel advance detection does not fully utilize the illumination characteristics and angular domain information of the observation system, resulting in inaccurate imaging of anomalies ahead of the tunnel. To achieve more accurate tunnel advance detection results, PSF correction can be incorporated into tunnel advance detection. It should be noted that directly calculating the PSF using traditional Fourier transforms is computationally expensive, so a more efficient angular domain method is used to calculate the PSF. Summary of the Invention

[0005] The purpose of the present invention is to provide a tunnel advance detection method based on point spread function correction. The 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. The angle domain point spread function is used to correct the tunnel reverse time migration imaging results. The angle domain point spread function contains the observation system information in the tunnel and the information between the inclination angle of the anomaly and the incident angle and scattering angle of the wave field near the anomaly. Therefore, more accurate tunnel migration imaging results can be obtained.

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

[0007] Step 1: Set the shot points and detection points to collect the observed seismic data ;in represents the time variable; Represent the positions of the detection point and the shot point respectively;

[0008] Step 2: Construct a rectangular grid tunnel geological model for numerical simulation and set the number of horizontal grid points of the rectangular grid tunnel geological model. and the number of vertical grid points , and set the spatial sampling interval for forward simulation, including the lateral sampling interval and longitudinal sampling interval , time sampling interval , maximum sampling time ;

[0009] Step 3: Given a velocity model describing the velocity distribution of underground medium wave propagation , the velocity model c is calculated based on the horizontal sampling interval and longitudinal sampling interval Discretize to the grid defined in step 2, use the source signal, numerically solve the wave equation, use the time sampling interval dt, and perform forward propagation in the time range [0, Nt] to obtain the forward simulation from the shot point To a point in space The original TRON ;

[0010] Step 4: The observed earthquake data collected in step 1 Back propagation is performed in time to generate the back propagation field from the detection point , will be transmitted to the TRON and back-propagation field The reverse time migration imaging results are obtained by calculating the cross-correlation imaging conditions. , as input for subsequent correction;

[0011] Step 5: Use the Poynting vector method to obtain the forward wave field from the shot point Propagation Angle , the back-propagation field from the detection point The transmission angle ;

[0012] Step 6: According to the propagation angle of the wave field and , the wave field can be decomposed into the superposition of beams, and then the local illumination matrix can be obtained ;

[0013] Step 7: Local lighting matrix Perform coordinate transformation from the incident / scattering angle of the wave field Convert to the tilted wavenumber domain ,in, The corresponding wave number representing the inclination angle of the anomaly body, Represents the tilt angle, the specific corresponding relationship is as follows Figure 2 As shown; then obtain the angle domain point spread function ;

[0014] Step 8: Take the reverse time migration image obtained in step 4 Convert to the wavenumber domain and use the point spread function Reverse time migration imaging results Correction is performed, and then the corrected result is converted back to the spatial domain to obtain the final corrected imaging result m psf .

[0015] The present invention is also characterized in that:

[0016] In step 3, TRON The following acoustic wave equation is satisfied:

[0017] (1)

[0018] In formula (1), For spatial location Medium velocity at is the medium density, It represents the speed at which seismic waves propagate in a medium.

[0019] To simplify the formula, omit ,use express.

[0020] In step 4, using back propagation seismic data Obtain the reverse propagation wave field from the detection point, and obtain the reverse time migration imaging result by calculating the zero-delay cross-correlation between the forward propagation wave field and the reverse propagation wave field at each time step , defined as follows:

[0021] (2)

[0022] In formula (2) Indicates the shooting point from the tunnel arrive The forwarding wave field, Indicates from To the detection point in the tunnel The back propagation wave field.

[0023] In order to remove low-frequency noise, the 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-propagation field from the detection point The calculation is done using the Poynting vector method, as follows:

[0025] (3)

[0026] In formula (3), and is the Poynting vector of the forward wave field The x and z components of and is the Poynting vector of the back propagation wave field The x and z components of and The definition is as follows:

[0027] (4)

[0028] In formula (4), Represents the spatial derivative of the wavefield.

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

[0030] According to equations (3) and (4), the propagation angle of the wave field has been obtained. Given a broadband source, such as Ricker wavelet, the wave field in equation (2) is and It can be expressed as a superposition of beams:

[0031] (5)

[0032] In the formula It describes the wave from the shot point position in the tunnel Spread to Green's function; It describes the wave from Propagated to the receiving point in the tunnel Green's function; is the frequency, is the source spectrum. Calculate the mean square amplitude in a short time interval T using equation (5):

[0033] (6)

[0034] In the process of wave field propagation, given that the main frequency remains unchanged, we can get:

[0035] (7)

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

[0037] (8)

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

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

[0040] (9)

[0041] In formula (9), k g and k i denote the scattered and incident wave numbers, respectively, k d =k i +k g is the wave number associated with the tectonic dip;

[0042] Perform coordinate transformation on equation (8) and convert the incident / scattering angle into Convert to wavenumber domain

[0043] (10)

[0044] In formula (10), is the Jacobian matrix, is the tilt angle, ;

[0045] First, the incident / scattering angle Convert to oblique / reflective angle , and the relationship is , ; and then convert to ,in ;this and The relationship between ;get:

[0046] (11)

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

[0048] (12)

[0049] In step 8, the tunnel RTM imaging result is considered as the convolution between the true tunnel velocity model and the point spread function, expressed 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 formula (14), M mig The reverse time migration imaging result The Fourier transform of

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

[0053] The present invention obtains conventional tunnel reverse time migration imaging results; uses the Poynting vector method to directly decompose the wave field from the shot point and the receiver point 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 to obtain more accurate results. Compared with commonly used tunnel seismic exploration methods, this method can significantly improve the imaging accuracy of the anomaly ahead. The reasons for the above advantages of this method are: first, because the point spread function correction is used, the blurring effect in traditional imaging can be suppressed; second, the angle domain point spread function contains information about the inclination angle of the anomaly ahead of the excavation face and the incident angle and scattering angle of the wave field near the anomaly, as well as information about the observation system in the tunnel, which can obtain more accurate results. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

[0056] Figure 3 It is the real velocity model of the tunnel in the model calculation example of the present invention;

[0057] Figure 4 Schematic diagram of the observation system in the tunnel in the model calculation example of the present invention;

[0058] Figure 5 This is the imaging result obtained using the traditional reverse time migration method in the model example of the present invention;

[0059] Figure 6 is the imaging result after correction using the method proposed in the present invention in the model example of the present invention; DETAILED DESCRIPTION

[0060] To make the objectives, technical solutions, and advantages of the embodiments of the present invention more clear, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts shall fall within the scope of protection of the present invention.

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

[0062] Step 1: Set the shot points and detection points to collect the observed seismic data ;in represents the time variable; Represent the positions of the detection point and the shot point respectively;

[0063] Step 2: Construct a rectangular grid tunnel geological model for numerical simulation calculations, which is used to discretize the velocity model and serve as the computational grid framework for this imaging algorithm. Set the number of horizontal grid points of the rectangular grid tunnel geological model. and the number of vertical grid points , and set the spatial sampling interval for forward simulation, 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 underground medium , the velocity model c is calculated based on the horizontal sampling interval and longitudinal sampling interval Discretize it onto the grid defined in step 2, use the source signal, numerically solve the wave equation, use the time sampling interval dt, and perform forward propagation in the time range [0, Nt] to obtain the forward simulation from the shot point To a point in space The original TRON ;

[0065] In step 3, TRON The following acoustic wave equation is satisfied:

[0066] (1)

[0067] In formula (1), For spatial location Medium velocity at is the medium density, It represents the speed at which seismic waves propagate in a medium.

[0068] To simplify the formula, omit ,use express.

[0069] Step 4: The observed earthquake data collected in step 1 Back propagation is performed in time to generate the back propagation field from the detection point , will be transmitted to the TRON and back-propagation field The reverse time migration imaging results are obtained by calculating the cross-correlation imaging conditions. , as input for subsequent correction;

[0070] In step 4, the reverse time migration imaging result 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, which is defined as follows:

[0071] (2)

[0072] In formula (2) represents the earthquake source x in the tunnel i To the positive wave field of x, Indicates the transmission from x to the receiving side x in the tunnel g The back propagation field;

[0073] In order to remove low-frequency noise, the Perform Laplace filtering;

[0074] Step 5: Use the Poynting vector method to obtain the forward wave field from the shot point Propagation Angle , the back-propagation field from the detection point The transmission angle ;

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

[0076] (3)

[0077] In formula (3), and is the Poynting vector of the forward wave field The x and z components of and is the Poynting vector of the back propagation wave field The x and z components of and The definition is as follows:

[0078] (4)

[0079] In formula (4), represents the spatial derivative of the wave field;

[0080] Step 6: According to the propagation angle of the wave field and , the wave field can be decomposed into the superposition of beams, and then the local illumination matrix can be obtained ;

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

[0082] According to equations (3) and (4), the propagation angle of the wave field has been obtained. Given a broadband source, such as Ricker wavelet, the wave field in equation (2) is and It can be expressed as a superposition of beams:

[0083] (5)

[0084] In the formula It describes the wave from the shot point position in the tunnel Spread to Green's function; It describes the wave from Propagated to the receiving point in the tunnel Green's function; is the frequency, is the source spectrum. Calculate the mean square amplitude in a short time interval T using equation (6):

[0085] (6)

[0086] In the process of wave field propagation, given that the main frequency remains unchanged, we can get:

[0087] (7)

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

[0089] (8)

[0090] In formula (8), represents the local wave number;

[0091] Step 7: Local lighting matrix Perform coordinate transformation from incident / scattering angle of wave field Convert to the tilted wavenumber domain ,in, The corresponding wave number representing the inclination angle of the anomaly body, Represents the tilt angle, the specific corresponding relationship is as follows Figure 2 As shown; then obtain the angle domain point spread function ;

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

[0093] (9)

[0094] In formula (9), k g and k i denote the scattered and incident wave numbers, respectively, k d =k i +k g is the wave number associated with the tectonic dip;

[0095] Perform coordinate transformation on equation (8) and convert the incident / scattering angle into Convert to wavenumber domain :

[0096] (10)

[0097] In formula (11), is the Jacobian matrix, is the tilt angle, ;

[0098] First, the incident / scattering angle Convert to oblique / reflective angle , and the relationship is , ; and then convert to ,in ;this and The relationship between ;get:

[0099] (11)

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

[0101] (12)

[0102] Step 8: Take the reverse time migration image obtained in step 4 Convert to the wavenumber domain and use the point spread function Reverse time migration imaging results Correction is performed, and then the corrected result is converted back to the spatial domain to obtain the final corrected imaging result m psf .

[0103] In step 8, the tunnel RTM imaging result is considered as the convolution between the true tunnel velocity model and the point spread function, expressed 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 formula (13), M mig The reverse time migration imaging result The Fourier transform of

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

[0107] Model example:

[0108] The specific implementation process of the present invention is applied to a tunnel model. The model has a total length of 200 meters horizontally and 150 meters vertically. The model contains low-velocity anomalies, different media, etc. The length of the tunnel body is 80 meters and the height is 10 meters. The velocity model is as follows Figure 3 shown.

[0109] The arrangement 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 face, there are 11 sources spaced 1 m apart and 6 receivers spaced 2 m apart.

[0110] Figure 5 Results of conventional reverse time migration imaging of the tunnel model. For the first layer of interface, artifacts can be seen, which will interfere with the judgment of the seismic information in front of the tunnel. For the cave in the second layer of medium, in the area indicated by the two white arrows, it can be seen that the interface imaging behind the cave is not obvious.

[0111] Figure 6This is the image obtained after point spread function correction (PSF) in steps 1 through 8 of the proposed method. For the first interface, artifacts are suppressed. For the cave, the images of the front and rear interfaces are more balanced and their positions are clearer. This demonstrates that the proposed method can effectively image the anomaly ahead 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, ordinary technicians in the field should understand that the specific implementation methods of the present invention can still be modified or replaced by equivalents. Any modification or equivalent replacement that does not depart from the spirit and scope of the present invention should be covered by 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: The steps include: Step 1: Set the detection points and shot points to collect the observed seismic data ;in represents the time variable; Represent the positions of the detection point and the shot point respectively; Step 2: Construct a rectangular grid tunnel geological model for numerical simulation and set the number of horizontal grid points of the rectangular grid tunnel geological model. and the number of vertical grid points , and set the spatial sampling interval for forward simulation, 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 underground medium , the velocity model According to the horizontal sampling interval and longitudinal sampling interval Discretize it onto the grid of the rectangular grid tunnel geological model constructed in step 2, use the source signal, numerically solve the wave equation, use the time sampling interval dt, and perform forward propagation in the time range [0, Nt] to obtain the forward simulation from the shot point To a point in space The original TRON ; Step 4: The observed earthquake data collected in step 1 Perform back propagation in time to generate the back propagation field from the detection point , will be transmitted to the TRON and back-propagation field The reverse time migration imaging results are obtained by calculating the cross-correlation imaging conditions. ; Step 5: Use the Poynting vector method to obtain the forward wave field from the shot point Propagation Angle and the back-propagation field from the detection point The transmission angle ; Step 6: According to the propagation angle of the wave field and , decompose the wave field into the superposition of beams, and then obtain the local illumination matrix ; Step 7: Local lighting matrix Perform coordinate transformation from the incident / scattering angle of the wave field Convert to the tilted wavenumber domain ,in, The corresponding wave number representing the inclination angle of the anomaly body, Represents the tilt angle, and then obtains the angle domain point spread function ; Step 8: Take the reverse time migration image obtained in step 4 Convert to the wavenumber domain and use the point spread function Reverse time migration imaging Correction is performed, and then the corrected 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 wave field The following acoustic wave equation is satisfied: (1)； In formula (1), For spatial location Medium velocity at is the medium density, It represents the speed of seismic waves propagating in the 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 result is obtained by calculating the zero-delay cross-correlation of the forward wave field and the reverse wave field at each time step. , defined as follows: (2)； In formula (2) Indicates the shooting point from the tunnel arrive The forwarding wave field, Indicates from To the detection point in the tunnel The back propagation field; In order to remove low-frequency noise, the Perform Laplace filtering.

4. The method for tunnel advance detection 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 forward wave field from the shot point Propagation Angle and the back-propagation field from the detection point The transmission angle , specifically including: (3)； In formula (3), and is the Poynting vector of the forward wave field The x and z components of and is the Poynting vector of the back propagation wave field The x and z components of and The definition is as follows: (4)； In formula (4), Represents the spatial derivative of the wavefield.

5. The method for tunnel advance detection based on point spread function correction according to claim 4, characterized in that: In step 6, the local illumination matrix is ​​obtained The specific method is as follows: According to equations (4) and (5), the propagation angle of the wave field has been obtained. Given a broadband source, the wave field in equation (2) is and Expressed as a superposition of beams: (5)； In the formula It describes the wave from the shot point position in the tunnel Spread to Green's function; It describes the wave from Propagated to the receiving point in the tunnel Green's function; is the frequency, is the source spectrum; calculate the mean square amplitude in a short time interval T for equation (5): (6)； In the process of wave field propagation, given that the main frequency remains unchanged, we get: (7)； Based on formula (7), the local illumination matrix is ​​obtained: (8)； In formula (8), 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 is defined as follows: (9)； In formula (9), k g and k i denote the scattered and incident wave numbers, respectively, k d =k i +k g is the wave number associated with the tectonic dip; Perform coordinate transformation on equation (8) and convert the incident / scattering angle into Convert to wavenumber domain : (10)； In formula (10), is the Jacobian matrix, is the tilt angle, ; First, the incident / scattering angle Convert to oblique / reflective angle , and the relationship is , ; and then convert to ,in ; this and The relationship between ;get: (11)； Substituting formula (11) into formula (9) yields the angle domain point spread function: (12)。 7. The method for tunnel advance detection based on point spread function correction according to claim 6, characterized in that: In step 8, the tunnel RTM imaging result is considered as the convolution between the true tunnel velocity model and the point spread function, expressed 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 formula (13), M mig The reverse time migration imaging result The Fourier transform of Through M psf Perform inverse Fourier transform to obtain the corrected imaging result m psf .

Images