Earthquake surface wave array tomography method, storage medium and equipment

Abstract

The invention discloses a seismic surface wave array tomography method, a storage medium and equipment, and relates to the technical field of surface wave tomography, and the method comprises the steps: calculating a cross-correlation function between all seismic station pairs according to continuous waveform data recorded by seismic stations; determining the size of the sub-array according to the distance and the number of the seismic stations; processing all the cross-correlation functions in each sub-array by using linear Radon transform to obtain a frequency dispersion curve of each sub-array; respectively moving the positions of the sub-arrays according to a certain step length along the horizontal direction and the vertical direction, and obtaining frequency dispersion curves of the sub-arrays at other positions by using linear Radon transformation; and taking the frequency dispersion curves of different sub-arrays as input data, and constructing a phase velocity value of each grid point in the research area by using a seismic array tomography method. According to the scheme of the invention, the imaging resolution of a seismic array processing technology can be improved.

Description

Technical Field

[0001] The present invention relates to the technical field of surface wave tomography, and in particular to a seismic surface wave array tomography method, storage medium and equipment. Background Art

[0002] Tomography based on surface wave signals has been widely used to study the crust and upper mantle structure. With the development of dense seismic arrays, array-based seismic data processing methods have become very effective in extracting high-quality broadband surface wave dispersion curves. In recent years, advances in seismic array data processing methods have helped to extract multimodal dispersion curves, which will provide better constraints for studying the deep structure of the Earth.

[0003] Compared with the single-station method or the dual-station method, the seismic array data processing technology can extract high-quality multi-modal broadband surface wave dispersion curves from seismic data. However, since the seismic array data processing technology uses all surface wave signals in a series of subarrays for superposition, it can usually only obtain the average phase velocity value of the corresponding subarray, and then interpolate the phase velocity values ​​of different subarrays to construct the surface wave phase velocity map of the study area. The resolution of the phase velocity map is limited by the size of the subarray and the spacing of the subarray movement, resulting in low imaging resolution of tomography. Summary of the invention

[0004] The purpose of the present invention is to solve the problem of low resolution of phase velocity structure obtained by seismic array data processing technology, and to propose a seismic surface wave array tomography method, comprising the following steps:

[0005] S1, obtaining original continuous waveform data recorded by different seismic stations, and preprocessing them to obtain segmented continuous waveform data;

[0006] S2, using the segmented continuous waveform data, calculate the cross-correlation function between all pairs of seismic stations;

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

[0008] S4, moving the position of the subarray along the horizontal direction and the vertical direction according to a certain step length to obtain subarrays at other positions, and obtaining dispersion curves of the subarrays at other positions based on the method of S3;

[0009] S5. Based on the dispersion curves of different sub-arrays, the seismic surface wave array tomography method is used to construct the phase velocity value of each grid point in the study area.

[0010] Furthermore, step S1 is specifically as follows:

[0011] S11, obtaining continuous waveform records containing ground vibration signals recorded by different seismic stations;

[0012] S12, segmenting the continuous waveform data recorded by each seismic station according to a certain time length;

[0013] S13, performing time domain regularization on the segmented continuous waveform data to obtain regularized data;

[0014] S14. Spectral whitening is performed on the regularized data in the frequency domain to obtain segmented continuous waveform data.

[0015] Furthermore, S2 is specifically:

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

[0017]

[0018] Among them, c AB (t) represents the cross-correlation function between the seismic station pair 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. Superimposing the cross-correlation functions of each pair of seismic stations at all times obtains the cross-correlation function of each pair of seismic stations.

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

[0021] S31, designing a shape of a specific size based on the distance between seismic stations and the number of stations with a certain location as the center, and selecting the stations within the shape to form a subarray;

[0022] S32. Use linear Radon transform to process the cross-correlation function between all pairs of stations in each subarray to obtain the dispersion curve of each subarray. The principle of linear Radon transform is as follows:

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

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

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

[0026] All the cross-correlation functions within the subarray are sorted and aligned according to the station distance to form a virtual source set. Then, the cross-correlation function is transformed into the frequency domain by linear Radon transform to obtain L and d, and the model m is obtained and calculated. The result is transformed from the frequency slowness domain to the frequency velocity domain to obtain the dispersion energy diagram. Based on the maximum amplitude of the dispersion energy diagram, the dispersion curve of the subarray is extracted.

[0027] Furthermore, S5 is specifically:

[0028] S51, extracting dispersion curves of all sub-arrays at a specific period;

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

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

[0031] Construct the relationship between the geological structure of the study area and the observed data:

[0032]

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

[0034] The following formula is used to invert the phase velocity value of each grid point in the study area using the dispersion curves of all subarrays at a specific period as input data;

[0035]

[0036] Among them, Δm represents the initial model m 0 The change in the phase velocity coefficient relative to the current model m, G represents the relationship between the geological structure of the grid and the observed data, C nn represents the covariance matrix, Δd represents the difference between the observed phase and the theoretical phase calculated based on the current model m, and C mm represents the prior error covariance matrix of model m, m represents the current phase velocity coefficient model, m 0 Represents the initial model of the phase velocity coefficients.

[0037] The present invention further provides a computer-readable storage medium, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the above-mentioned surface wave tomography method based on array data is implemented.

[0038] The present invention also proposes an electronic device, comprising a processor and a memory, wherein the processor and the memory are interconnected, wherein the memory is used to store a computer program, the computer program comprises computer-readable instructions, and the processor is configured to call the computer-readable instructions to execute the above-mentioned array data-based surface wave tomography method.

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

[0040] The present invention proposes a seismic surface wave array tomography method, which calculates the cross-correlation function between any two stations from the continuous waveform data recorded by the seismic station, obtains the high-quality dispersion curve of each sub-array by using the array data processing technology, establishes the corresponding relationship between the sub-array dispersion curve and the velocity value under each grid point in the study area, and constructs the velocity value under each grid point in the study area by using the seismic array surface wave tomography method. Compared with the seismic array data processing technology that directly interpolates the phase velocity value of the sub-array to construct the surface wave phase velocity map of the study area, the technology of the present invention constructs the velocity value under each grid point in the study area by the seismic array tomography method, which can greatly improve the resolution of the internal structure of the earth constructed based on the seismic array processing technology. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

[0043] Figure 3 is a schematic diagram of the cross-correlation principle of seismic stations according to an embodiment of the present invention;

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

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

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

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

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

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

[0050] S1. Obtain continuous waveform data recorded by different seismic stations. For the schematic diagram of continuous waveform data, refer to Figure 2 , and preprocess to obtain segmented continuous waveform data. Specifically:

[0051] S11, obtaining continuous waveform records recorded by different seismic stations, the records containing ground vibration signals;

[0052] S12, segmenting the continuous waveform data recorded by each seismic station according to a certain time length;

[0053] S13, performing time domain regularization on the segmented continuous waveform data to obtain regularized data;

[0054] S14. Spectral whitening is performed on the regularized data in the frequency domain to obtain segmented continuous waveform data.

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

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

[0057]

[0058] Among them, c AB (t) represents the cross-correlation function between the seismic station pair A and B at time t, where t also represents the observation time, and 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. Superimposing the cross-correlation functions of each pair of seismic stations at all times obtains the cross-correlation function of each pair of seismic stations.

[0060] Schematic diagram of the cross-correlation principle of seismic stations Figure 3 , A and B in the figure represent seismic station A and seismic station B. The black arrows in the figure represent random vibrations inside the earth. When these vibrations are recorded by stations A and B at the same time, they contain the velocity information of the medium between stations A and B. Then, through the cross-correlation technology, the continuous waveform data recorded by stations A and B are processed to obtain the cross-correlation function between stations A and B. 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 according to the spacing and number of seismic stations; use linear Radon transform to process all cross-correlation functions in each subarray to obtain the dispersion curve of the subarray.

[0062] S31. All seismic stations in the study area represent an array, e.g. Figure 4 In (a), some seismic stations are selected according to a certain shape, which is called a subarray. According to the spacing and number of seismic stations in the study area, a certain position 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. Use linear Radon transform to process the cross-correlation function between all pairs of stations in each subarray to obtain the dispersion curve of each subarray. The principle of linear Radon transform is as follows:

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

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

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

[0067] All the cross-correlation functions within the subarray are sorted and aligned according to the station distance to form a virtual source set. Then, the cross-correlation function is transformed into the frequency domain by linear Radon transform to obtain L and d, and the model m is obtained and calculated. The result is transformed from the frequency slowness domain to the frequency velocity domain to obtain the dispersion energy diagram. Based on the maximum amplitude of the dispersion energy diagram, the dispersion curve of the subarray is extracted.

[0068] refer to Figure 4 ,in Figure 4 (a) is a schematic diagram of the seismic array and seismic subarray in the study area. The size of the subarray is determined by the 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 obtained surface wave dispersion curve, that is, the lower its resolution. Figure 4 (b) is Figure 4 Schematic diagram of the cross-correlation function between all pairs of stations in the subarray (a). Figure 4 (c) is the dispersion energy diagram of the subarray obtained using linear Radon transform.

[0069] S4, moving the subarray position along the horizontal direction and the vertical direction according to a certain step length to obtain subarrays at other positions, the moving step length is determined by the spacing between seismic stations and the number of stations, and obtaining dispersion curves of subarrays at other positions according to the method of S3. Figure 5 shown. Figure 5 There are 4 circular subarrays in the system, each of which contains 4 seismic stations, one of which is the initial subarray, and the other 3 are obtained by moving the initial subarray horizontally and vertically according to a certain step length. By moving the subarray positions to form 4 subarrays, the information of the intermediate anomaly can be obtained through seismic array tomography ( Figure 5 which part of the middle shadow).

[0070] S5. Based on the dispersion curves of different subarrays, the seismic surface wave array tomography method is used to invert the phase velocity value (velocity structure below the underground medium) of each grid point in the study area.

[0071] S51, extract the dispersion curves of all sub-arrays at a specific period; for example, Figure 5 The average speed of the four rings is extracted.

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

[0073] Construct the relationship between the geological structure of the study area and the observed data:

[0074]

[0075] Among them, C pred represents the average phase velocity value corresponding to the subarray, that is, Figure 6 The velocity value at the position of the five-pointed star; G(i) represents the relationship between the geological structure of the i-th grid and the observed data, C i Represents the velocity value of the i-th grid, where the grid corresponds to Figure 6 The black dots in the figure represent the number of inversion grids. In the embodiment of the present invention, n represents Figure 6 The number of black dots in the middle circle, Figure 6 Schematic diagram of the principle of seismic surface wave array tomography according to an embodiment of the present invention, Figure 6 The elliptical shaded part in the figure represents the influence range of the ray path in the subarray on the subarray velocity structure. If the existing array data processing technology is used, only the velocity value of the point corresponding to the middle five-pointed star can be obtained, but after using the seismic surface wave array tomography technology of the present invention, the velocity value corresponding to each black dot in the circle can be obtained.

[0076] The following formula is used to invert the phase velocity value of each grid point in the study area using the dispersion curves of all subarrays as input data;

[0077]

[0078] Among them, Δm represents the initial model m 0 The change in the phase velocity coefficient relative to the current model m, G represents the relationship between the geological structure of the grid and the observed data, C nn represents the covariance matrix, Δd represents the difference between the observed phase and the theoretical phase calculated based on the current model m, and C mm represents the prior error covariance matrix of model m, m represents the current phase velocity coefficient model, m 0 Represents the initial model of the phase velocity coefficients.

[0079] In an exemplary embodiment, a computer-readable storage medium is included, wherein the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, the above-mentioned surface wave tomography method based on array data is implemented.

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

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

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

Claims

1. Based on a seismic surface wave array tomography method, characterized in that: The following steps are involved: S1, obtaining original continuous waveform data recorded by different seismic stations, and preprocessing them to obtain segmented continuous waveform data; S2, using the segmented continuous waveform data, calculate the cross-correlation function between all pairs of seismic stations; S3. Determine the size of the subarray according to the spacing and number of seismic stations; use linear Radon transform to process all cross-correlation functions in each subarray to obtain the dispersion curve of the subarray; S4, moving the position of the subarray along the horizontal direction and the vertical direction according to a certain step length to obtain subarrays at other positions, and obtaining dispersion curves of the subarrays at other positions based on the method of S3; S5. Based on the dispersion curves of different sub-arrays, the phase velocity value of each grid point in the study area is constructed using the seismic surface wave array tomography method.

2. A seismic surface wave array tomography method according to claim 1, characterized in that: Step S1 is specifically as follows: S11, obtaining continuous waveform records containing ground vibration signals recorded by different seismic stations; S12, segmenting the continuous waveform data recorded by each seismic station according to a certain time length; S13, performing time domain regularization on the segmented continuous waveform data to obtain regularized data; S14. Spectral whitening is performed on the regularized data in the frequency domain to obtain segmented continuous waveform data.

3. A seismic surface wave array tomography method according to claim 1, characterized in that: S2 is specifically: S21. Calculate the cross-correlation function of each same time period between all pairs of stations using the background noise cross-correlation. The expression is as follows: Among them, c AB (t) represents the cross-correlation function between the seismic station pair 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. Superimposing the cross-correlation functions of each pair of seismic stations at all times obtains the cross-correlation function of each pair of seismic stations.

4. The seismic surface wave array tomography method according to claim 1, characterized in that: Step S3 is specifically as follows: S31, designing a shape of a specific size based on the distance between seismic stations and the number of stations with a certain location as the center, and selecting the stations within the shape to form a subarray; S32. Use linear Radon transform to process the cross-correlation function between all pairs of stations in each subarray to obtain the dispersion curve of each subarray. The principle of linear Radon transform is as follows: First, define the Radon inverse transform: d = Lm, where d represents the Fourier transform value of the waveform data of a certain frequency, m represents the Radon transform result of d, and L is the Radon transform matrix; The inversion problem of Radon transform is: Among them, W m represents the weighting matrix of the Radon transform model; The optimization problem of Radon transform is: The solution to the optimization problem is: Where I represents the identity matrix, W d is the data weighting matrix; yes The transposed matrix of ; the parameter λ balances the data fit difference and model weight; All the cross-correlation functions within the subarray are sorted and aligned according to the station distance to form a virtual source set. Then, the cross-correlation function is transformed into the frequency domain by linear Radon transform to obtain L and d, and the model m is obtained and calculated. The result is transformed from the frequency slowness domain to the frequency velocity domain to obtain the dispersion energy diagram. Based on the maximum amplitude of the dispersion energy diagram, the dispersion curve of the subarray is extracted.

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

6. A 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: Construct the relationship between the geological structure of the study area and the observed data: Among them, C pred represents the average phase velocity value corresponding to the subarray, G(i) represents the relationship between the geological structure of the i-th grid and the observed data, and C i represents the velocity value of the i-th grid, and n represents the number of inversion grids; The following formula is used to invert the phase velocity value of each grid point in the study area using the dispersion curves of all subarrays at a specific period as input data; Where Δm represents the change in the 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 observed data, and C nn represents the covariance matrix, Δd represents the difference between the observed phase and the theoretical phase calculated based on the current model m, and C mm represents the prior error covariance matrix of model m, m represents the current phase velocity coefficient model, and m0 represents 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, the method according to any one of claims 1 to 6 is implemented.

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

Images