A method for removing artifacts based on multi-mode surface wave dispersion of array passive source

Abstract

The application discloses a kind of based on array passive source multi-mode surface wave dispersion's artifact removal method, it is related to the field of seismic imaging, including the following steps: the relative distance and azimuth between stations in array are calculated;The noise cross-correlation function between stations is calculated;According to the noise cross-correlation function between stations, causal inter-spectral density matrix is constructed;The relative distance and azimuth between stations and causal inter-spectral density matrix are brought into the expression of modified weighted cross-correlation beamforming analysis or modified correction beamforming analysis, to remove the artifact in beam pattern and dispersion chart.The application effectively reduces the appearance of artifact in dispersion chart, significantly improves the accuracy of multi-mode dispersion curve extraction and mode recognition, and then makes the result of subsequent earth velocity structure inversion more accurate and reliable.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of seismic imaging, and particularly relates to a method for removing artifacts based on multi-mode surface wave dispersion of array passive sources. BACKGROUND

[0002] Since the 21st century, the emergence of ambient noise imaging method is undoubtedly a major breakthrough with milestone significance in the field of seismological research. Unlike the traditional research path based on natural earthquake event information, this method relies on ambient noise data and uses seismic interference method to effectively extract surface wave dispersion information and then realize accurate inversion of the earth's velocity structure. Ambient noise imaging method has unique advantages in the study of the earth's internal structure due to its excellent repeatability, high-resolution imaging effect and convenient and efficient operation process. In recent years, with the rapid development of node-type seismograph technology, a large number of dense and super-dense seismic arrays have been deployed worldwide. The establishment of these arrays greatly enriches the channels of seismic data acquisition and provides more abundant data support for seismological research, which has effectively promoted the in-depth study and wide application of array method.

[0003] In practical application, the deployment of seismic arrays faces a series of constraints. Restricted by the current technical level and actual geographical environment, geological conditions and other factors, it is difficult to realize high-density deployment of seismic arrays. And in the actual operation process, in order to ensure the stability and consistency of data acquisition, the seismic array usually adopts a regular distribution mode. Although this mode is convenient for data acquisition and processing to some extent, it also brings some problems. Due to the similarity of the distance and direction between stations, the array repeatedly samples at fixed distances and directions, which makes it difficult to fully and comprehensively sample the spatial noise field.

[0004] The limitations in the array layout and sampling process will directly affect the results of subsequent data processing and analysis. Taking the weighted cross-correlation beamforming analysis (WCBF) method and the modified beamforming analysis (MCBF) method as examples, a large number of artifacts often appear in the beam pattern and dispersion pattern generated by them. The existence of these artifacts seriously interferes with the accurate extraction of the dispersion curve and the accurate identification of the surface wave mode. The dispersion curve and the surface wave mode recognition are crucial for the inversion of the earth velocity model, which indirectly affects the accuracy and reliability of the inversion results of the earth velocity model. At present, CN 116400406A proposes a passive source multi-mode surface wave dispersion curve extraction method based on array, which can extract multi-mode surface wave dispersion curves with sufficient precision from short-time background noise data recorded by the station array. However, the method proposed in the patent still has many artifacts that interfere with the beam pattern and dispersion pattern. How to overcome these problems, further improve the resolution of the multi-mode dispersion pattern extracted by the station array method based on background noise data, and then improve the accuracy of the earth velocity model inversion has become a key problem that needs to be further studied and solved in the field of seismology. SUMMARY

[0005] In view of the above problems in the prior art, the array-based passive source multi-mode surface wave dispersion artifact removal method provided by the present application solves the interference problem of artifacts and power branches in the dispersion pattern.

[0006] In order to achieve the above-mentioned purpose of the application, the technical scheme adopted by the present application is as follows: an array-based passive source multi-mode surface wave dispersion artifact removal method, comprising the following steps:

[0007] S1, calculating the relative distance and orientation between stations in the array;

[0008] S2, calculating the noise cross-correlation function between stations;

[0009] S3, constructing a causal cross-spectral density matrix according to the noise cross-correlation function between stations;

[0010] S4, bringing the relative distance and orientation between stations and the causal cross-spectral density matrix into the expression of the modified weighted cross-correlation beamforming analysis or the modified modified beamforming analysis to remove the artifacts in the beam pattern and the dispersion pattern.

[0011] Further, the specific method for calculating the noise cross-correlation function between stations in step S2 is as follows: obtaining noise records of each station pair; dividing the noise records of each station pair into noise records of different time periods with a fixed length of time window and step; calculating the noise record cross-correlation of each time period of the station pair and superimposing to obtain the noise cross-correlation function between stations;

[0012] The noise cross-correlation function between stations is expressed as:

[0013] cor ij (ω) = <d * (x i , ω) * d(x j , ω)>,

[0014] where i represents the i-th station in the array, j represents the j-th station in the array; cor ij (ω) represents the noise cross-correlation function of the station i and the station j in the frequency domain, x i represents the vector position of the station i, x j represents the vector position of the station j; ω represents the angular frequency; d(x i , ω) is the Fourier spectrum of the noise record of the station i, d(x j , ω) is the Fourier spectrum of the noise record of the station j; * is the conjugate, and <> represents the long-time average.

[0015] Further, the calculation of the noise record cross-correlation of the stations for each time period is specifically: the noise record of the stations for each time period is normalized in the time domain by using the sliding average to obtain the first normalized result; the first normalized result is subjected to Fourier transform to obtain the Fourier transformed result; the Fourier transformed result is subjected to sliding average normalization to obtain the second normalized result; and the second normalized result is multiplied in the frequency domain to obtain the noise record cross-correlation of the stations for each time period.

[0016] Further, the specific method of step S3 is: Fourier transforming cor ij (ω) to obtain the noise cross-correlation function of the station i and the station j in the time domain; taking the part of the noise cross-correlation function of the station i and the station j in the time domain with t≥0 to constitute the causal response, and taking the part of the noise cross-correlation function of the station i and the station j in the time domain with t≤0 to constitute the acausal response; constructing a causal cross-spectral density matrix based on the causal response and the acausal response according to the station pair, and the expression is:

[0017]

[0018] wherein, is the causal cross-spectral density matrix, θ ij represents the azimuth between the station i and the station j; cor ij (t) represents the causal response between the station i and the station j, cor ij (-t) represents the acausal response between the station i and the station j; t represents the time variable; N represents the number of stations in the array; and FT{} represents the Fourier transform.

[0019] Further, the modified weighted cross-correlation beamforming analysis in step S4 is a linear combination of two base search results, and the expression is:

[0020] WCBF(k,ω,θ) = (CC + SS) / 2

[0021]

[0022] θ - π / 2 ≤ θ ij ≤ θ + π / 2

[0023] Wherein, WCBF(k,ω,θ) represents the modified weighted cross-correlation beamforming analysis result, N represents the number of stations in the array; r ij represents the distance between station i and station j; θ ij represents the azimuth between station i and station j; k represents the wave number of the surface wave, k = ω / c, c represents the phase velocity of the surface wave; Re is the real part, and Im is the imaginary part; θ represents the search azimuth; N(θ) represents the number of stations satisfying θ - π / 2 ≤ θ ij ≤ θ + π / 2, and π is 180 degrees; CC and SS are intermediate parameters; is the causal cross-spectral density matrix; ω represents the angular frequency cos is the cosine function, and sin is the sine function.

[0024] Further, the modified correction beamforming analysis in step S4 is a linear combination of two base search results, and the expression is:

[0025] MCBF(k,ω) = (XX + YY) / 2

[0026]

[0027] Wherein, MCBF(k,ω) represents the modified correction cross-correlation beamforming analysis result, and YY and XX are intermediate parameters.

[0028] The beneficial effects of the present application are:

[0029] 1. By dividing the noise cross-correlation function in the time domain into causal signals and non-causal signals and constructing a causal cross-spectral density matrix, the relationship between signals can be better analyzed, and a dispersion graph with less artifact influence can be obtained.

[0030] 2. Using the modified correction cross-correlation beamforming analysis method and the modified weighted cross-correlation beamforming analysis method, i.e. using a linear combination of two base search results, the artifacts in the dispersion graph can be more efficiently eliminated, and the stability and reliability of the power branch signal can be greatly improved. BRIEF DESCRIPTION OF DRAWINGS

[0031] Figure 1 is a flowchart of the method of the present application;

[0032] Figure 2 Known earth velocity model plot for the example described;

[0033] Figure 3 Known multi-modal dispersion plot for the example described;

[0034] Figure 4 Station distribution plot for the example described;

[0035] Figure 5 Large number of point sources plot for the example described;

[0036] Figure 6 Station distance and azimuth plot for the example described;

[0037] Figure 7 Partial station noise record plot for the example described using a large number of point sources simulation;

[0038] Figure 8 Computed noise cross-correlation function for a pair of stations for the example described;

[0039] Figure 9 Causal cross-spectral density matrix plot;

[0040] Figure 10 Beam plot for the example described using modified weighted cross-correlation beamforming at 5 Hz;

[0041] Figure 11 Beam plot for the example described using modified weighted cross-correlation beamforming at 7.5 Hz;

[0042] Figure 12 Beam plot for the example described using prior art at 5 Hz;

[0043] Figure 13 Beam plot for the example described using prior art at 7.5 Hz;

[0044] Figure 14 Azimuth-averaged dispersion plot for the example described using modified corrected cross-correlation beamforming directly;

[0045] Figure 15 Azimuth-averaged dispersion plot for the example described using prior art;

[0046] Figure 16 Azimuth-averaged dispersion plot for the example described using Figure 14 Chosen multi-modal dispersion plot;

[0047] Figure 17 ChinArray subarray plot;

[0048] Figure 18 Noise cross-correlation function obtained by using the present application in Comparative Experiment 1;

[0049] Figure 19 Frequency dispersion diagram obtained by using the present application in Comparative Experiment 1;

[0050] Figure 20 Frequency dispersion diagram obtained by using the prior art in Comparative Experiment 1;

[0051] Figure 21 Schematic diagram of Tongzhou subarray;

[0052] Figure 22 Noise cross-correlation function obtained by using the present application in Comparative Experiment 2;

[0053] Figure 23 Frequency dispersion diagram obtained by using the present application in Comparative Experiment 2;

[0054] Figure 24 Frequency dispersion diagram obtained by using the prior art in Comparative Experiment 2;

[0055] Figure 25 Schematic diagram of Tongshan subarray;

[0056] Figure 26 Noise cross-correlation function obtained by using the present application in Comparative Experiment 3;

[0057] Figure 27 Frequency dispersion diagram obtained by using the present application in Comparative Experiment 3;

[0058] Figure 28 Frequency dispersion diagram obtained by using the prior art in Comparative Experiment 3. DETAILED DESCRIPTION

[0059] The specific embodiments of the present application are described below to facilitate the understanding of the present application for those skilled in the art, but it should be clear that the present application is not limited to the scope of the specific embodiments, and for those skilled in the art, it is obvious that various changes are within the spirit and scope of the present application defined and determined by the appended claims, and all the inventions utilizing the concept of the present application are within the scope of protection.

[0060] As shown in Figure 1 The present application provides a method for removing artifacts based on array passive source multi-mode surface wave dispersion, which comprises the following steps:

[0061] S1, calculate the relative distance and direction between stations in the array;

[0062] In one embodiment of the present application, the specific implementation of the present application is introduced by taking the simulated noise data generated by taking the earth velocity model and multi-mode dispersion curve as known conditions as an example, wherein the earth velocity model is shown in Figure 2 , and the multi-mode dispersion curve is shown in Figure 3 . The simulated data simulates the noise record of the velocity record simulation station by using a large number of point sources, and the station array contains 144 stations, the transverse distribution of which is shown in Figure 4 , and the large number of point sources are shown in Figure 5 . The simulated station array position is in the Cartesian coordinate system, and the relative position of the station pair can be directly calculated, as shown in Figure 6 , wherein the two black points give the relative distance and direction of station 00 and station 99. For the longitude and latitude data of the actual laid station array, the relative distance and direction of the station pair need to be calculated according to the ellipsoidal coordinate system.

[0063] S2, calculating the noise cross-correlation function between stations: obtaining the noise record of each station pair; dividing the noise record of each station pair into noise records of different time periods by using a fixed length time window and step; calculating the noise record cross-correlation of each time period of the station pair and stacking to obtain the noise cross-correlation function between stations.

[0064] Part of the simulated station noise record by using a large number of point sources is shown in Figure 7 , and the noise record in the part of the station noise record is divided into different time periods by using a time window of 20s and a step of 10s. For each station pair, the noise record cross-correlation in each time period is calculated: the noise record of each time period of the station pair is normalized in the time domain by using sliding average to obtain the first normalized result; the first normalized result is subjected to Fourier transform to obtain the Fourier transformed result; the Fourier transformed result is subjected to sliding average normalization to obtain the second normalized result; the second normalized result is multiplied in the frequency domain to obtain the noise record cross-correlation of each time period of the station pair. The noise record cross-correlation of each time period is stacked to obtain the noise cross-correlation function of the station pair. Figure 8 The noise cross-correlation function of part of the station pairs calculated is shown in , wherein the black solid line is the noise cross-correlation function between station 00 and station 99, and according to the distance between the station pairs, the separation of different mode wave packets can be seen, and the trend is indicated by a black dotted line. For the actual data, the corresponding parameters need to be set according to the pre-extracted dispersion curve.

[0065] The noise cross-correlation function between stations is expressed as:

[0066] cor ij (ω)=<d * (x i ,ω)d(xj ,ω)>

[0067] Where i represents the i-th station in the array, j represents the j-th station in the array; cor ij (ω) represents the noise cross-correlation function between stations i and j in the frequency domain, x i represents the vector position of station i, x j represents the vector position of station j; ω represents the angular frequency; d(x i ,ω) is the Fourier spectrum of the noise record of station i, d(x j ,ω) is the Fourier spectrum of the noise record of station j, * means taking the conjugate, and <> means taking the long-term average.

[0068] S3, constructing the causal cross-spectral density matrix based on the noise cross-correlation function between stations;

[0069] Cor ij (ω) is Fourier transformed to obtain the noise cross-correlation function of station i and station j in the time domain; the part of the noise cross-correlation function of station i and station j in the time domain with t≥0 is taken to constitute the causal response, and the part of the noise cross-correlation function of station i and station j in the time domain with t≤0 is taken to constitute the non-causal response; based on the causal response and non-causal response, the causal cross-spectral density matrix is ​​constructed according to the station pair, as follows: Figure 9 As shown, its expression is:

[0070]

[0071] in, is the causal mutual spectral density matrix, θ ij represents the direction between station i and station j; cor ij (t) represents the causal response between station i and station j, cor ij (-t) represents the non-causal response between station i and station j; t represents the time variable; N represents the number of stations in the array; FT{} represents the Fourier transform.

[0072] S4. The relative distances and orientations between stations and the causal cross-spectral density matrix are introduced into the modified weighted cross-correlation bunching analysis or the modified corrected bunching analysis expression to remove artifacts in the beam pattern and dispersion diagram.

[0073] The modified weighted cross-correlation clustering analysis in step S4 is a linear combination of the two basis search results, and its expression is:

[0074] WCBF(k,ω,θ)=(CC+SS) / 2

[0075]

[0076] θ - π / 2 ≤ θ ij ≤ θ + π / 2

[0077] where WCBF(k, ω, θ) denotes the modified weighted cross-correlation beamforming result, N denotes the number of stations in the array; r ij denotes the distance between station i and station j; θ ij denotes the azimuth between station i and station j; k denotes the wave number of the surface wave, k = ω / c, c denotes the phase velocity of the surface wave; Re denotes the real part, and Im denotes the imaginary part; θ denotes the search azimuth; N(θ) denotes the number of stations satisfying θ - π / 2 ≤ θ ij ≤ θ + π / 2, and π denotes 180 degrees; CC and SS are intermediate parameters. is the causal cross-spectral density matrix; cos denotes the cosine function, and sin denotes the sine function.

[0078] The modified modified beamforming is a linear combination of two base search results, and the expression is as follows:

[0079] MCBF(k, ω) = (XX + YY) / 2

[0080]

[0081] where MCBF(k, ω) denotes the modified modified beamforming result, and YY and XX are intermediate parameters.

[0082] In the present embodiment, Figure 10 is the beam pattern obtained by using the modified weighted cross-correlation beamforming (WCBF) at 5 Hz, and the power represents the intensity of the incident surface wave in the X velocity and the Y velocity. Figure 11 is the beam pattern obtained by using the modified weighted cross-correlation beamforming (WCBF) at 7.5 Hz. Figure 12 is the beam pattern obtained by using the prior art at 5 Hz. Figure 13 is the beam pattern obtained by using the prior art at 7.5 Hz. Compared with the prior art, the power in the beam pattern obtained by using the present application changes more continuously along the azimuth, has less impact of artifacts, and has clearer mode branches. By superimposing the beam patterns along the azimuth for all frequencies and arranging them according to the frequency, an azimuth-averaged dispersion map without the impact of artifacts can be obtained.

[0083] Figure 14 is the azimuth-averaged dispersion map obtained by directly using the modified modified cross-correlation beamforming (MCBF), Figure 15 is the azimuth-averaged dispersion map obtained by using the prior art. As shown in the results in the figure, compared with the prior art, the dispersion map obtained by the method of the present application has less impact of artifacts, more continuous mode branches, and higher resolution.

[0084] After obtaining the dispersion map, the frequency-velocity values can be selected from the dispersion map along the power extreme branch to form a multi-mode dispersion curve, Figure 16 In order to utilize Figure 14 the selected multi-mode dispersion curve. The multi-mode dispersion curve is obtained and used for subsequent inversion of the earth velocity structure.

[0085] In order to verify the beneficial effects of the present application, the following comparative experiments are carried out for different subarrays.

[0086] Comparative Experiment 1: The dispersion maps corresponding to the ChinArray subarray are obtained by using the method proposed in the present application and the prior art respectively. As shown in Figure 17 , the ChinArray subarray is composed of 104 broadband seismographs, the station spacing is about 70 km, and the noise recording length is about 30 days. The time window length used for calculating the noise cross-correlation function is 1000s, and the step length is 500s, as shown in Figure 18 , the noise cross-correlation function can be observed to have obvious surface waveforms arranged according to the station pair distance. The dispersion map obtained by using the present application is as shown in Figure 19 , and the dispersion map obtained by using the prior art is as shown in Figure 20 .

[0087] Comparative Experiment 2: The dispersion maps corresponding to the Tongzhou subarray are obtained by using the method proposed in the present application and the prior art respectively. As shown in Figure 21 , the Tongzhou subarray is composed of 56 short-period seismographs, the station spacing is about 1 km, and the noise recording length is about 45 days. The time window length used for calculating the noise cross-correlation function is 300s, and the step length is 150s, as shown in Figure 22 , the noise cross-correlation function can be observed to have two obviously separated modes of surface waveforms arranged according to the station pair distance. The dispersion map obtained by using the present application is as shown in Figure 23 , and the dispersion map obtained by using the prior art is as shown in Figure 24 .

[0088] Comparative Experiment 3: The dispersion maps corresponding to the Tangshan subarray are obtained by using the method proposed in the present application and the prior art respectively. As shown in Figure 25 , the Tangshan subarray is composed of 30 node-type seismographs, the station spacing is about 30 m, and the noise recording length is about 8 days. The time window length used for calculating the noise cross-correlation function is 50s, and the step length is 25s, as shown in Figure 26 , the noise cross-correlation function can be observed to have surface waveforms arranged according to the station pair distance, but the wave packets are not separated. The dispersion map obtained by using the present application is as shown in Figure 27 , and the dispersion map obtained by using the prior art is as shown in Figure 28 .

[0089] Through the comparison and analysis of the present application and the prior art, it can be seen that the present application can significantly and efficiently eliminate the artifacts in the dispersion map, greatly reduce the result deviation and the possibility of mode error discrimination caused by the existence of artifacts. At the same time, the present application can also effectively reduce the disturbance phenomenon on the power branch, greatly improve the stability and reliability of the power branch signal. Through actual verification and data comparison, the present application has significant effect in improving the power branch resolution, can provide more accurate and reliable data support for subsequent data analysis and processing, and thus fundamentally improves the performance and application value of the whole method.

Claims

1. A method for removing artifacts based on array passive source multi-mode surface wave dispersion, characterized in that, The method comprises the following steps: S1, calculating the relative distance and azimuth between stations in the array; S2, calculating the noise cross-correlation function between stations; S3, constructing the causal cross-spectral density matrix according to the noise cross-correlation function between stations; S4, inputting the relative distance and azimuth between stations and the causal cross-spectral density matrix into the expression of the modified weighted cross-correlation beamforming analysis or the modified correction beamforming analysis to remove the artifacts in the beam pattern and dispersion pattern; The specific method of step S3 is: i and stations j Noise cross-correlation function Perform Fourier transform to obtain the station in the time domain i and stations j The noise cross-correlation function of the station in the time domain i and stations j The noise cross-correlation function t The part ≥0 constitutes the causal response, taking the station in the time domain i and stations j The noise cross-correlation function t The part ≤0 constitutes the non-causal response. Based on the causal response and the non-causal response, the causal cross-spectral density matrix is ​​constructed according to the station pair, and its expression is: wherein is the causal cross-spectral density matrix, denotes a station i and a station j between; denotes a causal response between a station i and a station j , denotes an acausal response between a station i and a station j ; t denotes a time variable; N denotes the number of stations in the array; denotes a Fourier transform; denotes an angular frequency; The modified weighted cross-correlation beamforming analysis in step S4 is a linear combination of two base search results, and the expression is as follows: wherein, represents the modified weighted cross-correlation beam analysis result, N represents the number of stations in the array; represents the distance between the stations i and j ; represents the azimuth between the stations i and j ; represents the wave number of the surface wave, , c represents the phase velocity of the surface wave; Re is the real part, and Im is the imaginary part; represents the searched azimuth; represents the number of stations satisfying , π is an angle of 180 degrees; CC and SS are intermediate parameters; is the causal cross-spectral density matrix; represents the angular frequency; cos is the cosine function, and sin is the sine function.

2. The method of claim 1, wherein, The specific method for calculating the noise cross-correlation function between stations in step S2 is as follows: obtaining the noise record of each station pair; dividing the noise record of each station pair into noise records of different time periods with a fixed length of time window and step; calculating the noise record cross-correlation of each time period of the station pair and superimposing to obtain the noise cross-correlation function between stations; The noise cross-correlation function between stations is expressed as: where i denotes the j-th station of the array, i denotes the j-th station of the array, j denotes the j-th station of the array, j denotes the j-th station of the array; denotes the noise cross-correlation function of the stations i and j in the frequency domain, denotes the vector position of the station i , denotes the vector position of the station j ; denotes the angular frequency; is the Fourier spectrum of the noise record of the station i , is the Fourier spectrum of the noise record of the station j , * is taken conjugate, denotes the long-time average.

3. The method of claim 2, wherein, The specific method for calculating the noise record cross-correlation of each time period of the station pair is as follows: normalizing the noise record of each time period of the station pair in the time domain by using the sliding average to obtain the first normalized result; performing Fourier transform on the first normalized result to obtain the Fourier transformed result; performing sliding average normalization on the Fourier transformed result to obtain the second normalized result; conjugate multiplying the second normalized result in the frequency domain to obtain the noise record cross-correlation of each time period of the station pair.

4. The method of claim 1, wherein, The modified correction beamforming analysis in step S4 is a linear combination of two base search results, and the expression is as follows: wherein, represents the modified corrected bunch analysis result, both YY and XX are intermediate parameters.

Images