A cross-scale surface wave full waveform inversion method based on seismic background noise

Abstract

The invention discloses a cross-scale surface wave full waveform inversion method based on seismic background noise, comprising: S1, blindly sampling seismic background noise data; S2, pre-processing seismic background noise data; S3, extracting Green's function; S4, calculating and obtaining virtual shot gather records; S5, performing surface wave forward simulation to obtain surface wave synthetic shot gather records; S6, constructing a global cross-correlation target function; S7, deriving a gradient formula, and calculating the gradient of the target function with respect to shear wave velocity; S8: calculating the step length, and iteratively updating the current velocity model; S9, judging whether the updated velocity model meets the convergence condition, and if so, executing S10; S10, outputting the shear wave velocity result of the sub-wave inversion of the current frequency band; if it is not the last frequency band data, executing S11; if so, executing S12; S11, increasing the sub-wave frequency band width, setting the shear wave velocity inversion result of the previous frequency band as the initial model, and executing S4-10; S12, outputting the shear wave velocity field.

Description

Technical Field

[0001] The invention relates to the field of geophysical exploration, and in particular to a cross-scale surface wave full waveform inversion method based on seismic background noise. Background Art

[0002] In geological surveys for engineering projects such as highways, rail transit, and urban underground space development, the shallow surface shear wave velocity is a very accurate representation of geological structures and abnormal targets, and is a very important and stable physical parameter. The seismic surface wave method based on background noise records the seismic noise in traditional cognition and converts it into "effective" information. It can adaptively solve interference from human activities, electromagnetics, and transportation vehicles, and does not require artificial seismic sources. It has the advantages of repeatability, green and non-destructiveness, and good observation timeliness. It is currently the most popular technical means to obtain the shear wave velocity of underground space media.

[0003] However, the traditional seismic background noise surface wave method exploration technology simply uses the dispersion characteristics of surface waves to invert the shear wave velocity, and the velocity curve obtained by inversion is the average velocity or equivalent velocity within the entire arrangement or array, which has huge defects in accuracy. Full waveform inversion is based on the wave equation theory. It aims to use all the effective information contained in the wave field, such as wave field travel time, amplitude, dispersion, phase, etc., to directly extract the underground medium parameters from the waveform record excited by the artificial source. It is a geophysical detection technology with extremely high accuracy. Therefore, it is very necessary to combine the advantages of both and seek a seismic detection technology with green and non-destructive data acquisition, strong anti-interference ability, and high imaging resolution. There are currently no reports on this aspect. Summary of the invention

[0004] The purpose of the present invention is to provide a cross-scale surface wave full waveform inversion method based on seismic background noise, which has green and lossless data acquisition, strong anti-interference ability and high imaging resolution.

[0005] To this end, the present invention adopts the following technical solutions:

[0006] A cross-scale surface wave full waveform inversion method based on seismic background noise comprises the following steps:

[0007] S1, arrange node-type seismographs at equal intervals according to the pre-designed linear observation system, and start blindly collecting seismic background noise data for a certain length of time;

[0008] S2, preprocessing the seismic background noise data to enhance and maintain the same phase between multiple channel signals;

[0009] S3, according to the positional relationship between the virtual shot point s and the detection point r in the linear observation system, the Green function G(s, r, t) is extracted using the interference technology:

[0010]

[0011] Where G(s,r,t) is the empirical Green's function between the virtual shot point s and the receiver point r, and τ is the related time shift;

[0012] S4, given a seismic wavelet w(t), use this given seismic wavelet w(t) and the Green function G(s,r,t) obtained in S3 to perform convolution calculation and obtain the virtual shot gather record u obs (s,r,t):

[0013] u obs (s,r,t)=w(t)*G(s,r,t) (4);

[0014] S5, based on the initial shear wave velocity field model and the seismic wavelet given in S4, the surface wave forward simulation is performed using the finite difference forward simulation method to obtain the surface wave synthetic shot gather record u syn (s,r,t);

[0015] S6, according to the virtual shot gather record u obs (s, r, t) and the surface wave synthetic shot gather record u syn (s,r,t), construct the global cross-correlation objective function:

[0016]

[0017] in, represents the sum of the global cross-correlations between all shot points and receiver points, and m represents the shear wave velocity model;

[0018] S7, derive the corresponding gradient formula and calculate the gradient of the objective function with respect to the shear wave velocity:

[0019]

[0020] S8, calculating the step size and iteratively updating the current velocity model;

[0021] S9, judging whether the iteratively updated velocity model satisfies the convergence condition; if so, executing S10; if not, repeating S5-9 using the updated shear wave velocity model;

[0022] S10, output the shear wave velocity result of the sub-wave inversion of the current frequency band; if it is not the last frequency band data, execute S11; if it is the last frequency band data, execute S12;

[0023] S11, increase the sub-waveband width, set the shear wave velocity inversion result of the previous band as the initial model, and execute S4-10;

[0024] S12, the program is executed and the final shear wave velocity field is output.

[0025] Preferably, the certain time in S1 is ≥30 minutes; and the preprocessing in S1 includes bandpass filtering and One-Bit normalization.

[0026] The One-Bit normalization is:

[0027]

[0028]

[0029] Where s(t) represents the background noise signal, index1 and index2 represent the index positions where the background noise amplitude is greater than 0 and less than 0, respectively.

[0030] In one embodiment of the present invention, the effective frequency band range of the bandpass filter is 5-30 Hz; the first stage wavelet frequency band is 5-10 Hz, the second stage wavelet frequency band is 5-20 Hz, and the third stage wavelet frequency band is 5-30 Hz.

[0031] In another embodiment of the present invention, the effective frequency band of the bandpass filter is in the range of 5-25 Hz; the effective frequency band in S4 is divided into two stages: the first stage sub-wave band is 5-15 Hz, and the second stage sub-wave band is 5-25 Hz.

[0032] The cross-scale surface wave full waveform inversion method based on seismic background noise of the present invention is mainly used to obtain the shear wave velocity in the engineering survey area with high precision, so as to locate abnormal geological structures such as soil-rock boundaries, active faults, karst, and boulders, optimize the geological interpretation effect, and ensure the safety of design and construction.

[0033] Compared with the prior art, the present invention has the following beneficial effects:

[0034] (1) High detection accuracy: The full waveform inversion method of the present invention has the advantages of repeatability, green and non-destructiveness, and good observation timeliness in terms of acquisition; in terms of processing, the full waveform inversion method makes full use of the full wave field information such as surface wave travel time, amplitude, phase dispersion, etc., and has obvious advantages in inversion accuracy and precision compared with traditional dispersion curve inversion;

[0035] (2) High inversion stability: The cross-scale surface wave full waveform inversion adopted by the present invention is a step-by-step updating process from large scale (low-frequency band data) to small scale (high-frequency band data), which effectively improves the stability of waveform inversion and avoids the problem of the iterative process falling into local extreme values. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1It is a flow chart of the cross-scale surface wave full waveform inversion method based on background noise of the present invention;

[0037] Figure 2 The figure shows 21 seismic background noises collected simultaneously in a certain arrangement in Example 1;

[0038] Figure 3 for Figure 2 The empirical Green's function of the background noise in the middle is obtained with the first trace as the shot point;

[0039] Figure 4 and Figure 5 They are the sub-waves of different frequency bands and Figure 3 Virtual shot gather after convolution of Green's function;

[0040] Figure 6 and Figure 7 These are the cross-scale shear wave velocity inversion results of the first and second stages in implementation 1, respectively. DETAILED DESCRIPTION

[0041] The technical solution of the present invention is described in detail below in conjunction with the accompanying drawings and embodiments.

[0042] See also Figure 1 The cross-scale surface wave full waveform inversion method based on seismic background noise of the present invention comprises the following steps:

[0043] S1, according to the pre-designed linear observation system, multiple node-type seismographs are arranged at equal intervals, and the seismic background noise data for a certain length of time are blindly collected. The "certain time" is generally at least 30 minutes. In order to ensure the signal-to-noise ratio, the collection time can be appropriately extended; the observation system should pre-design the corresponding virtual shot point s-detection point r position relationship to achieve the best detection effect on the underground target body.

[0044] The smaller the spacing between node seismometers, the higher the detection accuracy. However, considering production efficiency and cost, when designing the observation system, the spacing between seismometers should be determined based on the size of the detection target body known in advance to ensure that the spacing between detectors can meet the detection target requirements.

[0045] S2, preprocessing the original seismic background noise data, including: (1) bandpass filtering, which is determined according to the engineering detection depth and the background noise frequency band range, generally 5-30 Hz is the effective frequency band range; (2) One-bit normalization, the main purpose of which is to enhance and maintain the phase coherence between multi-channel signals:

[0046]

[0047]

[0048] Where s(t) represents the background noise signal, index1 and index2 represent the index positions where the background noise amplitude is greater than 0 and less than 0, respectively.

[0049] s(t) is a time series that fluctuates between positive and negative. Formula (1) is used to find the index of positive and negative values ​​in the signal s(t), that is, at what position s(t) is positive or negative (it can be understood as finding the positive and negative values ​​in it); Formula (2) is used to assign 1 to the positive values ​​found and -1 to the negative values.

[0050] S3, according to the positional relationship between the virtual shot point s and the detector point r in the observation system, the empirical Green's function is extracted using the interference technology:

[0051]

[0052] Where G(s,r,t) is the empirical Green's function between the virtual shot point s and the receiver point r, and τ is the related time shift.

[0053] S4, given a seismic wavelet w(t), use this given seismic wavelet w(t) and the Green function G(s,r,t) obtained in S3 to perform convolution calculation and obtain the virtual shot gather record u obs (s,r,t):

[0054] u obs (s,r,t)=w(t)*G(s,r,t) (4)

[0055] In the formula, the frequency band range of the given seismic wavelet w(t) gradually widens with the advancement of the cross-scale inversion process. Taking the frequency band range of 5-30 Hz described in S2 as an example, during the implementation process, it can be divided into the first stage wavelet frequency band of 5-10 Hz, the second stage wavelet frequency band of 5-20 Hz, and the third stage wavelet frequency band of 5-30 Hz.

[0056] S5, based on the initial shear wave velocity field model and the seismic wavelet given in S4, the surface wave forward simulation is performed using the finite difference forward simulation method to obtain the surface wave synthetic shot gather record u syn (s,r,t);

[0057] S6, constructing an objective function based on the virtual shot gather record obtained in S4 and the surface wave synthetic shot gather record obtained in S5. In order to reduce the influence of waveform amplitude on the objective function, the present invention constructs a global cross-correlation objective function:

[0058]

[0059] in, represents the sum of the global cross-correlations between all shot points and receiver points, and m represents the shear wave velocity model;

[0060] S7, derive the corresponding gradient formula and calculate the gradient of the objective function with respect to the shear wave velocity:

[0061]

[0062] S8, calculating the step size and iteratively updating the current velocity model;

[0063] S9, judging whether the iteratively updated velocity model satisfies the convergence condition; if so, executing S10; if not, repeating S5-9 using the updated shear wave velocity model;

[0064] S10, output the shear wave velocity model of the current frequency band wavelet inversion; if it is not the last frequency band data, execute S11; if it is the last frequency band data, execute S12;

[0065] S11, increase the sub-waveband width, set the shear wave velocity inversion result of the previous stage as the initial model, and execute S4-10;

[0066] S12, the program is executed and the final shear wave velocity field is output.

[0067] Example 1

[0068] In order to verify the application effect of the method and system of the present invention, an application test was carried out in a certain engineering survey area. The regional engineering overview is: through the preliminary survey, it is known that there are small underground cavities in the area, and the cavities are shallowly buried. For this reason, a survey line with a total length of 34m was laid, and the interval between adjacent seismometers was 1m to detect the depth and position of the cavities in the area, so as to verify the correctness of the method of the present invention.

[0069] S1, using node-type digital seismograph to observe seismic background noise data, the continuous acquisition time is 30 minutes. The frequency band range of this test is 5-25Hz, Figure 2 Shown is 21 seismic background noises collected simultaneously in a certain arrangement.

[0070] S2, preprocessing the original seismic background noise obtained in S1, including: (1) bandpass filtering, the frequency band range in this test is set to 5-25 Hz; (2) one-bit normalization to enhance the phase consistency between multi-channel signals.

[0071] S3, based on the virtual shot point-receiver point position relationship, the Green function is extracted using interferometry technology. Figure 3 for Figure 2 The background noise in the middle is the empirical Green's function obtained with the first seismic background noise as the shot point.

[0072] S4, from S1, we know that the frequency band range of this test is 5-25Hz. The implementation process can be divided into two stages. The first stage sub-wave frequency band is 5-15Hz, and the second stage sub-wave frequency band is 5-25Hz. Figure 4 and Figure 5 are the sub-waves and Green's functions of different frequency bands ( Figure 3 ) Virtual shot gather after convolution.

[0073] S5, based on the initial shear wave velocity field model and the first stage wavelet (frequency band 5-15 Hz), the finite difference forward modeling method is used to perform surface wave forward simulation to obtain the surface wave synthetic shot gather record.

[0074] S6, constructing a global cross-correlation objective function based on the virtual shot gather records and the surface wave synthetic shot gather records obtained in S5.

[0075] S7, derive the corresponding gradient formula and calculate the gradient of the objective function with respect to the shear wave velocity.

[0076] S8, calculate the step size and iteratively update the current velocity model.

[0077] S9, judging whether the iteratively updated velocity model satisfies the convergence condition; if so, executing S10; if not, repeating S5-9 using the updated shear wave velocity model.

[0078] S10, output the current shear wave velocity result; if it is not the last frequency band data, execute S11; if it is the last frequency band data, execute S12.

[0079] S11, increase the sub-wave bandwidth and enter the next stage of inversion. At this time, the shear wave velocity inversion result of the previous stage is the initial model, and S4-10 is further executed.

[0080] S12, the program is executed and the final shear wave velocity field is output. Figure 6 and Figure 7 These are the cross-scale shear wave velocity inversion results of the first and second stages of this embodiment, respectively. It can be seen that as the frequency band becomes wider, the accuracy of the inversion result becomes higher.

Claims

1. A cross-scale surface wave full waveform inversion method based on seismic background noise, characterized in that: The following steps are involved: S1, arrange node-type seismographs at equal intervals according to the pre-designed linear observation system, and start blindly collecting seismic background noise data for a certain length of time; S2, preprocessing the seismic background noise data to enhance and maintain the same phase between multiple channel signals; S3, according to the positional relationship between the virtual shot point s and the detection point r in the linear observation system, the Green function G(s, r, t) is extracted using the interference technology: Where G(s,r,t) is the empirical Green's function between the virtual shot point s and the receiver point r, and τ is the related time shift; S4, given a seismic wavelet w(t), use this given seismic wavelet w(t) and the Green function G(s,r,t) obtained in S3 to perform convolution calculation and obtain the virtual shot gather record u obs (s,r,t): u obs (s,r,t)=w(t)*G(s,r,t) (4); S5, based on the initial shear wave velocity field model and the seismic wavelet given in S4, the surface wave forward simulation is performed using the finite difference forward simulation method to obtain the surface wave synthetic shot gather record u syn (s,r,t); S6, according to the virtual shot gather record u obs (s, r, t) and the surface wave synthetic shot gather record u syn (s,r,t), construct the global cross-correlation objective function: in, represents the sum of the global cross-correlations between all shot points and receiver points, and m represents the shear wave velocity model; S7, derive the corresponding gradient formula and calculate the gradient of the objective function with respect to the shear wave velocity: S8, calculating the step size and iteratively updating the current velocity model; S9, judging whether the iteratively updated velocity model satisfies the convergence condition; if so, executing S10; if not, repeating S5-9 using the updated shear wave velocity model; S10, output the shear wave velocity result of the sub-wave inversion of the current frequency band; if it is not the last frequency band data, execute S11; if it is the last frequency band data, execute S12; S11, increase the sub-waveband width, set the shear wave velocity inversion result of the previous band as the initial model, and execute S4-10; S12, the program is executed and the final shear wave velocity field is output.

2. The cross-scale surface wave full waveform inversion method based on seismic background noise according to claim 1 is characterized in that: The certain time mentioned in S1 is ≥30 minutes.

3. The cross-scale surface wave full waveform inversion method based on seismic background noise according to claim 1 is characterized in that: The preprocessing described in S1 includes bandpass filtering and One-Bit normalization.

4. The cross-scale surface wave full waveform inversion method based on seismic background noise according to claim 3 is characterized in that: The One-Bit normalization is: Where s(t) represents the background noise signal, index1 and index2 represent the index positions where the background noise amplitude is greater than 0 and less than 0, respectively.

5. The cross-scale surface wave full waveform inversion method based on seismic background noise according to claim 3 is characterized in that: The effective frequency band of the bandpass filter is 5-30 Hz.

6. The cross-scale surface wave full waveform inversion method based on seismic background noise according to claim 5, characterized in that: The effective frequency band in S4 is divided into three stages: the first stage sub-waveband is 5-10Hz, the second stage sub-waveband is 5-20Hz, and the third stage sub-waveband is 5-30Hz.

7. The cross-scale surface wave full waveform inversion method based on seismic background noise according to claim 3 is characterized in that: The effective frequency band range of the bandpass filter is 5-25 Hz.

8. The cross-scale surface wave full waveform inversion method based on seismic background noise according to claim 7, characterized in that: The effective frequency band in S4 is divided into two stages: the first stage sub-wave band is 5-15Hz, and the second stage sub-wave band is 5-25Hz.

Images