Seismic surface wave array tomography method, storage medium, and device

Abstract

This invention discloses a seismic surface wave array tomography method, storage medium, and device, relating to the field of surface wave tomography technology. The method includes: calculating the cross-correlation function between all pairs of seismic stations based on continuous waveform data recorded by seismic stations; determining the size of the subarray based on the spacing and number of seismic stations; processing all cross-correlation functions within each subarray using a linear Radon transform to obtain the dispersion curve of each subarray; moving the subarray along both the horizontal and vertical directions with a certain step size, and using a linear Radon transform to obtain the dispersion curves of subarrays at other positions; using the dispersion curves of different subarrays as input data, and constructing the phase velocity values ​​of each grid point in the study area using the seismic array tomography method. The solution of this invention can improve the imaging resolution of seismic array processing technology.

Description

Technical Field

[0001] This invention relates to the field of surface wave tomography, and more particularly to a seismic surface wave array tomography method, storage medium, and device. Background Technology

[0002] Surface wave-based tomography has been widely used to study the structure of the Earth's crust and upper mantle. With the development of dense seismic arrays, array-based seismic data processing methods have become highly effective in extracting high-quality broadband surface wave dispersion curves. In recent years, advancements in seismic array data processing methods have facilitated the extraction of multimodal dispersion curves, providing better constraints for studying the Earth's deep structure.

[0003] Compared to single-station or dual-station methods, seismic array data processing techniques can extract high-quality multimodal broadband surface wave dispersion curves from seismic data. However, because seismic array data processing techniques superimpose surface wave signals from a series of subarrays, they typically only obtain the average phase velocity values ​​of the corresponding subarrays. Then, the phase velocity values ​​from different subarrays are interpolated to construct a surface wave phase velocity map of the study area. The resolution of this phase velocity map is limited by the size of the subarrays and the spacing between subarray movements, resulting in relatively low imaging resolution for tomographic imaging. Summary of the Invention

[0004] The purpose of this invention is to address the problem of low phase velocity structure resolution obtained from seismic array data processing techniques by proposing a seismic surface wave array tomography method, comprising the following steps:

[0005] S1. Obtain the raw continuous waveform data recorded by different seismic stations and perform preprocessing to obtain segmented continuous waveform data.

[0006] S2. Using segmented continuous waveform data, calculate the cross-correlation function between all seismic station pairs;

[0007] S3. Determine the size of the subarray based on the spacing and number of seismic stations; use linear Radon transform to process all cross-correlation functions within each subarray to obtain the dispersion curve of the subarray;

[0008] S4. Move the position of the subarray along the horizontal and vertical directions with a certain step size to obtain the subarrays at other positions. Obtain the dispersion curves of the subarrays at other positions based on the method in S3.

[0009] S5. Based on the dispersion curves of different subarrays, the phase velocity values ​​of each grid point in the study area are constructed using the seismic surface wave array tomography method.

[0010] Furthermore, step S1 specifically includes:

[0011] S11. Obtain continuous waveform records containing ground motion signals from different seismic stations;

[0012] S12. Divide the continuous waveform data recorded by each seismic station into segments according to a certain duration;

[0013] S13. Perform time-domain regularization on the segmented continuous waveform data to obtain the regularized data;

[0014] S14. Perform spectral analysis on the regularized data in the frequency domain to obtain segmented continuous waveform data.

[0015] Furthermore, S2 specifically refers to:

[0016] S21. Calculate the cross-correlation function between all station pairs for each same time period using the background noise cross-correlation. The expression is as follows:

[0017]

[0018] Among them, c AB (t) represents the cross-correlation function between seismic stations A and B at time t, f A (τ) represents the continuous waveform record at time τ recorded by seismic station A, f B 9τ+t) represents the continuous waveform record at time τ+t recorded by seismic station B;

[0019] S22. The cross-correlation function of each seismic station pair is obtained by superimposing the cross-correlation functions of each seismic station pair at all times.

[0020] Furthermore, step S3 specifically includes:

[0021] S31. Based on the spacing and number of seismic stations, design a shape of a specific size with a certain location as the center, and select the stations within the shape to form a sub-array.

[0022] S32. Using linear Radon transform, the cross-correlation function between all station pairs within each subarray is processed to obtain the dispersion curve of each subarray. The principle of linear Radon transform is as follows:

[0023] First, define the inverse Radon transform: d = Lm, where d represents the Fourier transform value of waveform data at a certain frequency, m represents the Radon transform result of d, and L is the Radon transform matrix;

[0024] The inversion problem of the Radon transform is: Among them, W m Represents the weighting matrix of the Radon transform model;

[0025] The optimization problem of the Radon transform is: The solution to the optimization problem is: Where I represents the identity matrix, W d It is a data weighted matrix; yes The transpose matrix; parameter λ balances the data fit difference and model weighting;

[0026] All cross-correlation functions within the subarray are sorted and aligned according to station distances to form a virtual source set. Then, the cross-correlation functions are transformed to the frequency domain through linear Radon transform to obtain L and d, resulting in model m. The model m is then calculated, and the results are transformed from the frequency slowness domain to the frequency velocity domain to obtain the dispersion energy map. Based on the maximum amplitude of the dispersion energy map, the dispersion curve of the subarray is extracted.

[0027] Furthermore, S5 specifically refers to:

[0028] S51. Extract the dispersion curves of all subarrays in a specific period;

[0029] S52. Based on the dispersion curves of all subarrays at a specific period, the phase velocity value of each grid point is inverted using the seismic surface wave array tomography method.

[0030] Furthermore, in S52, the principle of the seismic surface wave array tomography method is as follows:

[0031] Constructing the relationship between the geological structure of the study area and observational data:

[0032]

[0033] Among them, C pred Let G(i) represent the average phase velocity value corresponding to the subarray, and G(i) represent the relationship between the geological structure of the i-th grid and the observation data. i This represents the velocity value of the i-th grid, and n represents the number of inverted grids;

[0034] Using the following formula, with the dispersion curves of all subarrays in a specific period as input data, the phase velocity value of each grid point in the study area can be retrieved.

[0035]

[0036] Where Δm represents the change in phase velocity coefficient of the initial model m0 relative to the current model m, G represents the relationship between the geological structure of the grid and the observation data, and C nn Let C represent the covariance matrix, Δd represent the difference between the observed phase and the theoretical phase calculated based on the current model m, and C represent the covariance matrix. mm Let m represent the prior error covariance matrix of model m, m represent the current phase velocity coefficient model, and m0 represent the initial model of the phase velocity coefficient.

[0037] The present invention also proposes a computer-readable storage medium storing a computer program that, when executed by a processor, implements the above-described surface wave tomography method based on array data.

[0038] The present invention also proposes an electronic device, including 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 including computer-readable instructions, and the processor is configured to invoke the computer-readable instructions to execute the above-described surface wave tomography method based on array data.

[0039] The beneficial effects of the technical solution provided by this invention are:

[0040] This invention proposes a seismic surface wave array tomography method. From continuous waveform data recorded by seismic stations, the cross-correlation function between any two stations is calculated. High-quality dispersion curves for each subarray are obtained using array data processing techniques. A correspondence is established between the subarray dispersion curves and the velocity values ​​below each grid point in the study area. The seismic array surface wave tomography method is then used to construct the velocity values ​​below each grid point in the study area. Compared to seismic array data processing techniques that directly interpolate the phase velocity values ​​of the subarrays to construct a surface wave phase velocity map of the study area, this invention, by constructing the velocity values ​​below each grid point in the study area using seismic array tomography, can significantly improve the resolution of the reconstruction of the Earth's internal structure based on seismic array processing techniques. Attached Figure Description

[0041] Figure 1 This is a flowchart of a seismic surface wave array tomography method according to an embodiment of the present invention;

[0042] Figure 2 This is a schematic diagram of continuous seismic waveform data according to an embodiment of the present invention;

[0043] Figure 3 This is a schematic diagram illustrating the principle of seismic station cross-correlation in an embodiment of the present invention;

[0044] Figure 4 This is a schematic diagram of step S3 in an embodiment of the present invention, wherein... Figure 4 (a) shows a schematic diagram of the seismic array and subarrays in the study area. Figure 4 (b) is Figure 4 A schematic diagram of the cross-correlation function between all station pairs within the subarray in (a) is shown. Figure 4 (c) shows the dispersion energy diagram of the subarray obtained using the linear Radon transform;

[0045] Figure 5 This is a schematic diagram of array data processing technology according to an embodiment of the present invention;

[0046] Figure 6 This is a schematic diagram illustrating the principle of seismic surface wave array tomography in an embodiment of the present invention;

[0047] Figure 7 This is a block diagram of an electronic device according to an exemplary embodiment of the present invention. Detailed Implementation

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

[0049] A flowchart of a seismic surface wave array tomography method according to an embodiment of the present invention is shown below. Figure 1 The specific steps are as follows:

[0050] S1. Obtain continuous waveform data recorded by different seismic stations. (See diagram for continuous waveform data.) Figure 2 The data is then preprocessed to obtain segmented continuous waveform data. Specifically:

[0051] S11. Obtain continuous waveform records from different seismic stations, which contain ground motion signals;

[0052] S12. Divide the continuous waveform data recorded by each seismic station into segments according to a certain duration;

[0053] S13. Perform time-domain regularization on the segmented continuous waveform data to obtain the regularized data;

[0054] S14. Perform spectral analysis on the regularized data in the frequency domain to obtain segmented continuous waveform data.

[0055] S2. Using segmented continuous waveform data, calculate the cross-correlation function between all seismic station pairs.

[0056] S21. Calculate the cross-correlation function between all station pairs for each same time period using the background noise cross-correlation. The expression is as follows:

[0057]

[0058] Among them, c AB (t) represents the cross-correlation function between seismic stations A and B at time t, where t also represents the observation time, f A (τ) represents the continuous waveform record at time τ recorded by seismic station A, f B (τ+t) represents the continuous waveform record at time τ+t recorded by seismic station B.

[0059] S22. The cross-correlation function of each seismic station pair is obtained by superimposing the cross-correlation functions of each seismic station pair at all times.

[0060] A schematic diagram of the cross-correlation principle of seismic stations is shown below. Figure 3 In the figure, A and B represent seismic station A and seismic station B, respectively. The black arrows in the figure represent random vibrations inside the Earth. When these vibrations are recorded simultaneously by stations A and B, they contain velocity information of the medium between stations A and B. Then, by processing the continuous waveform data recorded by stations A and B through cross-correlation technology, the cross-correlation function between stations A and B can be obtained. This cross-correlation function can be used to study the velocity structure between stations A and B.

[0061] S3. Determine the size of the subarray based on the spacing and number of seismic stations; use linear Radon transform to process all cross-correlation functions within each subarray to obtain the dispersion curve of the subarray.

[0062] S31, All the seismic stations in the study area represent an array, such as Figure 4 In (a), some seismic stations are selected according to a certain shape, which is called a subarray. Based on the spacing and number of seismic stations in the study area, a certain location is selected as the center in the study area, a shape of a specific size is designed, and the seismic stations within the shape are selected to form a subarray.

[0063] S32. Using linear Radon transform, the cross-correlation function between all station pairs within each subarray is processed to obtain the dispersion curve of each subarray. The principle of linear Radon transform is as follows:

[0064] First, define the inverse Radon transform: d = Lm, where d represents the Fourier transform value of waveform data at a certain frequency, m represents the Radon transform result of d, and L is the Radon transform matrix;

[0065] The inversion problem of the Radon transform is: Among them, W m Represents the weighting matrix of the Radon transform model;

[0066] The optimization problem of the Radon transform is: The solution to the optimization problem is: Where I represents the identity matrix, W d It is a data weighted matrix; yes The transpose matrix; parameter λ balances the data fit difference and model weighting;

[0067] All cross-correlation functions within the subarray are sorted and aligned according to station distances to form a virtual source set. Then, the cross-correlation functions are transformed to the frequency domain through linear Radon transform to obtain L and d, resulting in model m. The model m is then calculated, and the results are transformed from the frequency slowness domain to the frequency velocity domain to obtain the dispersion energy map. Based on the maximum amplitude of the dispersion energy map, the dispersion curve of the subarray is extracted.

[0068] refer to Figure 4 ,in Figure 4 (a) shows a schematic diagram of the seismic array and seismic subarray in the study area. The size of the subarray is determined by the number of seismic stations that need to be included in the actual calculation. The larger the subarray size, the more cross-correlation functions it contains, and the larger the area represented by the surface wave dispersion curve, i.e., the lower its resolution. Figure 4 (b) is Figure 4 A schematic diagram of the cross-correlation function between all station pairs within the subarray in (a) is shown. Figure 4 (c) is the dispersion energy diagram of the subarray obtained by linear Radon transform.

[0069] S4. Move the subarray positions along the horizontal and vertical directions with a certain step size to obtain other subarray positions. The step size is determined by the spacing and number of seismic stations. Obtain the dispersion curves of other subarray positions according to the method in S3. A schematic diagram of the array data processing technology in this embodiment of the invention is shown below. Figure 5 As shown. Figure 5 There are four circular subarrays, each containing four seismic stations. One is the initial subarray, and the other three are subarrays located at different positions, obtained by moving the initial subarray horizontally and vertically in specific steps. These four subarrays are formed by moving the subarrays, and information about intermediate anomalies is obtained through seismic array tomography. Figure 5 (Which part is shaded in the middle?)

[0070] S5. Based on the dispersion curves of different subarrays, the phase velocity value (velocity structure beneath the subsurface medium) of each grid point in the study area is inverted using the seismic surface wave array tomography method.

[0071] S51. Extract the dispersion curves of all subarrays at a specific period; for example, extract the dispersion curves of all subarrays at a specific period. Figure 5 The average velocity of the four rings is extracted.

[0072] S52. Based on the dispersion curves of all subarrays at a specific period, the phase velocity value of each grid point is inverted using the seismic surface wave array tomography method. The principle of the seismic surface wave array tomography method is as follows:

[0073] Constructing the relationship between the geological structure of the study area and observational data:

[0074]

[0075] Among them, C pred This represents the average phase velocity value corresponding to the subarray, i.e. Figure 6 The velocity value at the location of the central pentagram; G(i) represents the relationship between the geological structure of the i-th grid and the observation data, C i This represents the velocity value of the i-th grid, where the grid corresponds to... Figure 6 The black dots in the diagram represent the number of inversion meshes. In this embodiment of the invention, n represents... Figure 6 The number of black dots inside the middle circle. Figure 6 This is a schematic diagram illustrating the principle of seismic surface wave array tomography in an embodiment of the present invention. Figure 6 The elliptical shaded area represents the range of influence of the ray path on the velocity structure of the subarray. Using existing array data processing techniques, only the velocity value corresponding to the point in the center pentagram can be obtained. However, using the seismic surface wave array tomography technique of this invention, the velocity value corresponding to each black dot within the circle can be obtained.

[0076] Using the following formula, with the dispersion curves of all subarrays as input data, the phase velocity value of each grid point in the study area can be retrieved.

[0077]

[0078] Where Δm represents the change in phase velocity coefficient of the initial model m0 relative to the current model m, G represents the relationship between the geological structure of the grid and the observation data, and C nn Let C represent the covariance matrix, Δd represent the difference between the observed phase and the theoretical phase calculated based on the current model m, and C represent the covariance matrix. mm Let m represent the prior error covariance matrix of model m, where m represents the current phase velocity coefficient model, and m0 represents the initial model of the phase velocity coefficient.

[0079] In one exemplary embodiment, a computer-readable storage medium is included, which stores a computer program that, when executed by a processor, implements the above-described surface wave tomography method based on array data.

[0080] Please see Figure 7 In one exemplary embodiment, the device further includes an electronic device including at least one processor, at least one memory, and at least one communication bus.

[0081] The memory stores a computer program, which includes computer-readable instructions. The processor calls the computer-readable instructions stored in the memory through the communication bus to execute the above-mentioned surface wave tomography method based on array data.

[0082] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those 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 invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A seismic surface wave array tomography method, characterized in that, Includes the following steps: S1. Obtain the raw continuous waveform data recorded by different seismic stations and perform preprocessing to obtain segmented continuous waveform data. S2. Using segmented continuous waveform data, calculate the cross-correlation function between all seismic station pairs; S3. Determine the size of the subarray based on the spacing and number of seismic stations; use linear Radon transform to process all cross-correlation functions within each subarray to obtain the dispersion curve of the subarray; S4. Move the position of the subarray along the horizontal and vertical directions with a certain step size to obtain the subarrays at other positions. Obtain the dispersion curves of the subarrays at other positions based on the method in S3. S5. Based on the dispersion curves of different subarrays, the phase velocity values ​​of each grid point in the study area are constructed using the seismic surface wave array tomography method.

2. The seismic surface wave array tomography method according to claim 1, characterized in that, Step S1 is as follows: S11. Obtain continuous waveform records containing ground motion signals from different seismic stations; S12. Divide the continuous waveform data recorded by each seismic station into segments according to a certain duration; S13. Perform time-domain regularization on the segmented continuous waveform data to obtain the regularized data; S14. Perform spectral analysis on the regularized data in the frequency domain to obtain segmented continuous waveform data.

3. The seismic surface wave array tomography method according to claim 1, characterized in that, S2 specifically refers to: S21. Calculate the cross-correlation function between all station pairs for each same time period using the background noise cross-correlation. The expression is as follows: Among them, c AB (t) represents the cross-correlation function between seismic stations A and B at time t, f A (τ) represents the continuous waveform record at time τ recorded by seismic station A, f B (τ+t) represents the continuous waveform record at time τ+t recorded by seismic station B; S22. The cross-correlation function of each seismic station pair is obtained by superimposing the cross-correlation functions of each seismic station pair at all times.

4. The seismic surface wave array tomography method according to claim 1, characterized in that, Step S3 is as follows: S31. Based on the spacing and number of seismic stations, design a shape of a specific size with a certain location as the center, and select the stations within the shape to form a sub-array. S32. Using linear Radon transform, the cross-correlation function between all station pairs within each subarray is processed to obtain the dispersion curve of each subarray. The principle of linear Radon transform is as follows: First, define the inverse Radon transform: d = Lm, where d represents the Fourier transform value of waveform data at a certain frequency, m represents the Radon transform result of d, and L is the Radon transform matrix; The inversion problem of the Radon transform is: Among them, W m Represents the weighting matrix of the Radon transform model; The optimization problem of the Radon transform is: The solution to the optimization problem is: Where I represents the identity matrix, W d It is a data weighted matrix; yes The transpose matrix; parameter λ balances the data fit difference and model weighting; All cross-correlation functions within the subarray are sorted and aligned according to station distances to form a virtual source set. Then, the cross-correlation functions are transformed to the frequency domain through linear Radon transform to obtain L and d, resulting in model m. The model m is then calculated, and the results are transformed from the frequency slowness domain to the frequency velocity domain to obtain the dispersion energy map. Based on the maximum amplitude of the dispersion energy map, the dispersion curve of the subarray is extracted.

5. The seismic surface wave array tomography method according to claim 1, characterized in that, S5 specifically refers to: S51. Extract the dispersion curves of all subarrays in a specific period; S52. Based on the dispersion curves of all subarrays at a specific period, the phase velocity value of each grid point is inverted using the seismic surface wave array tomography method.

6. The seismic surface wave array tomography method according to claim 5, characterized in that, In S52, the principle of the seismic surface wave array tomography method is as follows: Constructing the relationship between the geological structure of the study area and observational data: Among them, C pred Let G(i) represent the average phase velocity value corresponding to the subarray, and G(i) represent the relationship between the geological structure of the i-th grid and the observation data. i This represents the velocity value of the i-th grid, and n represents the number of inverted grids; Using the following formula, with the dispersion curves of all subarrays in a specific period as input data, the phase velocity value of each grid point in the study area can be retrieved. Where Δm represents the change in phase velocity coefficient of the initial model m0 relative to the current model m, G represents the relationship between the geological structure of the grid and the observation data, and C nn Let C represent the covariance matrix, Δd represent the difference between the observed phase and the theoretical phase calculated based on the current model m, and C represent the covariance matrix. mm Let m represent the prior error covariance matrix of model m, m represent the current phase velocity coefficient model, and m0 represent the initial model of the phase velocity coefficient.

7. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, it implements the method as described in any one of claims 1-6.

8. An electronic device, characterized in that, The device includes a processor and a memory, the processor being interconnected with the memory, wherein the memory is used to store a computer program, the computer program including computer-readable instructions, and the processor is configured to invoke the computer-readable instructions to perform the method as described in any one of claims 1-6.

Images