A method for synthesizing seismic records based on wave equation

Abstract

The present invention relates to earthquake numerical simulation technology, and specifically discloses a method for synthesizing earthquake records based on wave equations. The method comprises: step 1, introducing an extreme value obtaining function to obtain the position of wave crests and troughs; step 2, stratifying the earthquake data to obtain the longitudinal wave velocity, transverse wave velocity and density average value of the corresponding stratification; step 3, according to the formula of the vectorized reflectivity method, obtaining the slowness-frequency domain reflection coefficient; step 4, resampling the obtained slowness-frequency domain reflection coefficient to obtain the angle-frequency domain reflection coefficient; step 5, extracting the frequency domain wavelet; step 6, obtaining the angle-time domain synthetic earthquake record. The present invention uses actual earthquake data as basic data, introduces an improved automatic peak-finding algorithm, realizes the reasonable division of the number of layers of the underground medium, reduces the number of layers while ensuring the wave field characteristics, and combines the earthquake wavelet extracted by using a single-channel earthquake record to stably realize the wave field simulation of the underground medium.

Description

Technical Field

[0001] The present invention relates to the technical field of earthquake numerical simulation, and specifically comprises an earthquake record synthesis method based on wave equation. Background Art

[0002] Seismic numerical simulation is a method and means to simulate the propagation of waves in underground media based on the construction of underground models. Through seismic numerical simulation, the mechanism of seismic wave propagation in the actual underground medium can be constructed, and it helps to understand the real form of the underground. At present, numerical simulation plays a very important role in seismic exploration, seismology and various disciplines. There are many methods of numerical simulation. The method based on wave equation is one of the most practical and efficient methods. The vectorized reflectivity method has become the first choice of wave equation methods due to its precise wave field simulation characteristics and excellent computational efficiency.

[0003] The vectorized reflectivity method is a method that solves the elastic wave equation under the one-dimensional assumption to achieve wave field simulation in layered media. It can provide a complete solution to the wave field, accurately simulate multiple waves and conversion waves, and take into account the effects of energy loss during the propagation of seismic waves. For isotropic media, the vectorized reflectivity method does not involve the problem of decoupling of longitudinal and transverse waves, and both longitudinal and transverse waves can be completely separated. It is one of the most effective methods to achieve full wave field simulation of layered vertical inhomogeneous media. The application of the vectorized reflectivity method plays an extremely critical role in fine wave field simulation and pre-stack seismic inversion.

[0004] When using the vectorized reflectivity method for wave field simulation, it is necessary to manually estimate the number of stratified segments and the thickness of the underground medium for the existing seismic data. However, it is found in actual work that the manual estimation and stratification of seismic data with a large amount of data has a strong subjective influence, and the workload of data stratification for the entire work area is extremely large and almost impossible to achieve. In addition, if the entire work area uses the same stratification of well information, the lateral changes of the underground medium cannot be well characterized. Moreover, the number of layers calculated based on well information is much larger than the number of seismic layers, which greatly increases the calculation load and increases the calculation time several times. Therefore, reasonable stratification of seismic data is the key to improving the calculation speed of the vectorized reflectivity method, and it is also an urgent problem to be solved for accurate and rapid use of the vectorized reflectivity method for numerical simulation.

[0005] There are two main existing vector reflectivity layering methods: the first is to use acoustic logging data for division, but because the logging resolution is higher than seismic, the number of layers is extremely large and the layers are thin, which greatly reduces the calculation speed and cannot well characterize the lateral changes of the underground medium. The second is to artificially divide the strata based on prior knowledge. Its disadvantages are that the subjective factors are extremely strong and the workload is too large for actual data, making it difficult to implement. The above-mentioned layering method undoubtedly reduces the efficiency of the vector reflectivity method in forward modeling, hindering the application of this method in geophysical exploration. Summary of the invention

[0006] In order to reduce the subjectivity of stratification and improve the computational efficiency of forward modeling using the vectorized reflectivity method, the present invention proposes a seismic record synthesis method based on the wave equation. This method uses actual seismic data as basic data and combines it with an accurate peak value acquisition algorithm to achieve layer division of underground media, reduce the number of layers while ensuring the wave field characteristics, and achieve accurate and rapid simulation of underground media.

[0007] In order to solve the above technical problems, the technical solution adopted by the present invention is: a seismic record synthesis method based on wave equation, comprising:

[0008] Step 1, extracting a single-channel observed seismic record, introducing an extreme value finding function, and finding the peak and trough positions;

[0009] Step 2, stratify the seismic data according to the obtained wave crest and trough positions, substitute the P-wave velocity, S-wave velocity and density in the well logging data, and obtain the average P-wave velocity, S-wave velocity and density of the corresponding stratification;

[0010] Step 3, substituting the obtained P-wave velocity, S-wave velocity and density average value of the corresponding layer into the formula of the vector reflectivity method to obtain the slowness-frequency domain reflection coefficient;

[0011] Step 4, resampling the obtained slowness-frequency domain reflection coefficient to obtain the angle-frequency domain reflection coefficient;

[0012] Step 5, extracting frequency domain wavelets from the single-channel observed seismic record;

[0013] Step 6, multiplying the obtained angle-frequency domain reflection coefficient and the extracted frequency domain wavelet to obtain an angle-frequency domain synthetic seismic record, and then converting it to the angle-time domain to obtain an angle-time domain synthetic seismic record.

[0014] Preferably, the step 1 comprises: extracting a single-channel observed seismic record, substituting it into a peak value obtaining function to calculate the position of the wave crest, and then calculating the maximum value of the waveform, and subtracting it from the original waveform to obtain the position of the wave trough;

[0015] P=f(G(t)) (2),

[0016] G'(t)=G(t)-MAX(G(t)) (3),

[0017] L=f(G'(t)) (4),

[0018] In equations (2), (3), and (4), f represents the extreme value finding function, G(t) is the single-channel observed seismic record, P is the position of the crest, L is the position of the trough, and MAX represents the maximum value finding function.

[0019] Preferably, in step 2, the method for obtaining the average values ​​of the longitudinal wave velocity, the shear wave velocity and the density of the corresponding layer is:

[0020]

[0021] In formula (5), represents the average parameter value of the nth layer, N represents the position of the nth extreme point, M represents the position of the n+1th extreme point, and r i It represents any parameter value within this range, and the parameters include longitudinal wave velocity, shear wave velocity and density.

[0022] Preferably, in step 3, the method for obtaining the slowness-frequency domain reflection coefficient is:

[0023]

[0024] In formula (6), R pp (p,ω) represents the slowness-frequency domain reflection coefficient, p is the horizontal slowness, ω is the angular frequency, υ0 is a six-element vector, υ0(4) and υ0(1) represent the 4th and 1st elements of υ0 respectively;

[0025] υ0 is calculated as:

[0026] υ0=Q0Q1···Q N-1 υ N , (7),

[0027] υ n =Q n υ n+1 , (8),

[0028] υ N =[1 0 0 0 0 0] T . (10),

[0029]

[0030] In equations (6), (7), (8), (10), and (11), υ0 represents the total reflection response of the top layer of the medium, n represents any layer of the medium, N represents the total number of layers of the medium, and υ n+1 represents the response of the n+1th underground layer, υ N represents the initial response of the lowest layer of the medium, Q n represents the transfer matrix, and represents the influence of the nth layer of medium on the amplitude of seismic waves, The superscript + indicates the lower interface of the layer interface. The superscript - indicates the upper interface of the layer interface. and Each contains 16 independent components; the matrix E n Represents the change in phase of a seismic wave when it passes through the nth layer of medium.

[0031] Preferably, step 4 comprises:

[0032] Introducing p n = sinθ / α n , where p n represents the horizontal slowness of the nth layer, θ represents the incident angle, α n represents the longitudinal wave velocity of the nth layer;

[0033] Resample the slowness-frequency domain reflection coefficient to obtain the angle-frequency domain reflection coefficient:

[0034] R pp (θ,ω)=resample(R pp (p,ω)) (12),

[0035] In formula (12), resample represents the resampling process, R pp (θ, ω) represents the angle-frequency domain reflection coefficient.

[0036] Preferably, in step 5, the method of extracting the frequency domain wavelet from the single-channel observed seismic record is:

[0037] S(ω)=waveletextract(G(t)) (13),

[0038] In formula (13), S(ω) represents the frequency domain wavelet, G(t) represents the single-channel observed seismic record, and waveletextract represents the wavelet extraction function;

[0039] The wavelet extraction function is as follows:

[0040] Z = log(fft(G(t))) (14),

[0041] wl=ifft(Z)*W (15),

[0042] S(ω)=exp(fft(wl)); (16),

[0043] In formulas (14), (15), and (16), log represents logarithm, fft represents Fourier transform, ifft represents inverse Fourier transform, exp represents exponential, and W represents the time window range.

[0044] Preferably, step 6 comprises: multiplying the obtained angle-frequency domain reflection coefficient and the wavelet in the frequency domain to obtain an angle-frequency domain synthetic seismic record, and performing an inverse Fourier transform to the angle-time domain to obtain an angle-time domain synthetic seismic record:

[0045]

[0046] In formula (17), d(θ,t) is the angle-time domain synthetic seismic record, S(ω) represents the frequency domain wavelet, and R pp (θ, ω) is the angle-frequency domain reflection coefficient, and ω is the angular frequency.

[0047] The method of the present invention uses actual seismic data as basic data and introduces an improved automatic peak-finding algorithm to achieve a reasonable division of the number of layers of the underground medium, reducing the number of layers while ensuring the wave field characteristics. Combined with the seismic wavelets extracted from single-channel seismic records, the underground medium wave field simulation can be stably realized, thereby improving the calculation efficiency and accuracy of the forward modeling. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 The present invention provides a flow chart of an embodiment of a method for synthesizing seismic records based on a wave equation.

[0049] Figure 2 This is a schematic diagram of the P-wave velocity, S-wave velocity and density actually measured in a certain well.

[0050] Figure 3 Schematic diagram of the peak and trough positions obtained for a single-channel seismic record.

[0051] Figure 4 This is a schematic diagram of the parameter stratification results of a certain stratum based on the locations of peaks and troughs.

[0052] Figure 5 It is a schematic diagram of the wavelet (time domain) extracted from a single-channel seismic record.

[0053] Figure 6a This is a schematic diagram of synthesizing seismic records using traditional methods.

[0054] Figure 6b It is a schematic diagram of synthesizing seismic records by the method of the present invention.

[0055] Figure 6c It is a residual diagram of the seismic records synthesized by the traditional method and the seismic records synthesized by the method of the present invention. DETAILED DESCRIPTION

[0056] In order to facilitate a further understanding of the method of the present invention and the effects achieved, the present invention is described in detail below with reference to the accompanying drawings. However, the accompanying drawings are only provided for reference and illustration and are not intended to limit the present invention.

[0057] The present invention constructs a seismic record synthesis method based on wave equation. Since the peaks and troughs of a single-channel waveform of a seismic record represent the stratum interface experienced by the seismic wave during propagation, the extreme points of the waveform can be used as the basis for dividing the strata. Under the horizontal layered assumption, the strata between the two interfaces have the same seismic wave propagation velocity and density. By introducing and improving the peak search algorithm, an accurate waveform extreme value search result can be obtained. According to the position of the extreme point, the corresponding stratum thickness is obtained, and the longitudinal wave velocity, shear wave velocity and density data are layered, and the reflection coefficient is obtained by substituting the vectorized reflectivity method, and multiplied with the obtained seismic wavelet in the frequency domain to realize the forward simulation of the underground medium.

[0058]

[0059] In formula (1), R pp (p, ω) represents the slowness-frequency domain reflection coefficient, S(ω) is the frequency domain wavelet, d(p, t) represents the slowness-time domain synthetic seismic record, e is a natural constant, ω is the angular frequency, and i is an imaginary unit. From formula (1), it can be seen that the accuracy of the synthetic seismic record is closely related to the reflection coefficient and the wavelet.

[0060] Furthermore, according to Snell's law, slowness P is equal to sinθ / α, where θ is the angle and α is the seismic wave velocity. The slowness-time domain synthetic seismic record can be converted into an angle-time domain synthetic seismic record. Only by converting the synthetic seismic record to the angle domain can it be compared and calculated with the actual seismic record.

[0061] like Figure 1 As shown, the process of an embodiment of a seismic record synthesis method based on a wave equation provided by the present invention includes:

[0062] Step 1, extracting a single-channel observed seismic record, introducing an extreme value finding function, and finding the peak and trough positions;

[0063] Step 2, stratify the seismic data according to the obtained wave crest and trough positions, substitute the P-wave velocity, S-wave velocity and density in the well logging data, and obtain the average P-wave velocity, S-wave velocity and density of the corresponding stratification;

[0064] Step 3, substituting the obtained P-wave velocity, S-wave velocity and density average value of the corresponding layer into the formula of the vector reflectivity method to obtain the slowness-frequency domain reflection coefficient;

[0065] Step 4, resampling the obtained slowness-frequency domain reflection coefficient to obtain the angle-frequency domain reflection coefficient;

[0066] Step 5, extracting frequency domain wavelets from the single-channel observed seismic record;

[0067] Step 6, multiplying the obtained angle-frequency domain reflection coefficient and the extracted frequency domain wavelet to obtain an angle-frequency domain synthetic seismic record, and then converting it to the angle-time domain to obtain an angle-time domain synthetic seismic record.

[0068] Step 1, extracting a single-channel observed seismic record, introducing an extreme value obtaining function, and obtaining the position of the wave crest and the wave trough, specifically comprising: extracting a single-channel observed seismic record, substituting the peak value obtaining function into the peak value obtaining function to calculate the position of the wave crest, and then calculating the maximum value of the waveform, and subtracting it from the original waveform to obtain the position of the wave trough;

[0069] P=f(G(t)) (2),

[0070] G'(t)=G(t)-MAX(G(t)) (3),

[0071] L=f(G'(t)) (4),

[0072] In formulas (2), (3), and (4), f represents the extreme value obtaining function, G(t) is the single-channel observed seismic record (i.e., the actual seismic record), P is the position of the crest, L is the position of the trough, and MAX represents the maximum value obtaining function.

[0073] Step 2: According to the desired wave crest and trough positions and their corresponding depths, the seismic data is layered, and the velocity and density values ​​in the well logging data are substituted to calculate the parameter values ​​(P-wave velocity, S-wave velocity and density) of each layer. The specific calculation process can be expressed as follows:

[0074]

[0075] In formula (5), represents the average parameter value of the nth layer, N represents the position of the nth extreme point, M represents the position of the n+1th extreme point, and r i Indicates any parameter value within this range, including longitudinal wave velocity, shear wave velocity and density;

[0076] Well logging data is obtained from actual well logging and is considered accurate data in geophysics. The relationship between the measured time and depth is the standard for verifying the accuracy of synthetic seismic records.

[0077] Step 3, according to the formula of the vector reflectivity method, the slowness-frequency domain reflection coefficient is obtained, which specifically includes:

[0078]

[0079] In formula (6), R pp (p, ω) represents the slowness-frequency domain reflection coefficient, p is the horizontal slowness, ω is the angular frequency, and the reflection coefficient is composed of the ratio of two elements of a six-element vector υ0. υ0(4) and υ0(1) represent the fourth and first elements of the six-element vector υ0, respectively. The calculation method of the six-element vector is as follows:

[0080] υ0=Q0Q1···Q N-1 υ N (7)

[0081] υ n =Q n υ n+1 (8)

[0082] Q n =E n F n (9)

[0083] υ N =[1 0 0 0 0 0] T (10)

[0084] In the above formulas, all subscripts represent the layer number, n represents any layer of the medium, and N represents the total number of layers of the medium; Q n represents the transfer matrix of the nth layer, the six-element vector υ n represents the response of the nth underground layer, υ n+1 represents the response of the n+1th underground layer, while υ0 represents the total reflection response of the top layer of the medium, υ N represents the initial response of the bottom layer of the medium, and the superscript T represents the transpose of the vector. n It expresses the change of the phase of the seismic wave when it passes through the nth layer of medium, F n The effect of the nth layer of medium on the amplitude of seismic waves is described;

[0085] In order to facilitate subsequent calculations, Q n Written in the following form

[0086]

[0087] and It represents the influence of the nth layer of medium on the amplitude of seismic waves. The superscript + represents the lower interface of the layer, and - represents the upper interface of the layer. Each of them contains 16 independent components. The matrix E n It expresses the change in phase of a seismic wave when it passes through the nth layer of medium.

[0088] Step 4: After obtaining the slowness-frequency domain reflection coefficient, the relationship between horizontal slowness and incident angle p is introduced. n = sinθ / α n , here p n represents the horizontal slowness of the nth layer, θ represents the incident angle, α n represents the longitudinal wave velocity of the nth layer, for R pp (p,ω) is resampled to obtain the reflection coefficient in the angle-frequency domain.

[0089] R pp (θ,ω)=resample(R pp (p,ω)) (12),

[0090] In formula (12), resample represents the resampling process, R pp (θ, ω) represents the angle-frequency domain reflection coefficient.

[0091] Step 5: Extract frequency domain wavelets from single-channel seismic records:

[0092] S(ω)=waveletextract(G(t)) (13),

[0093] In formula (13), S(ω) represents the frequency domain wavelet, G(t) represents the single-channel seismic record, and waveletextract represents the wavelet extraction function. The specific process is as follows:

[0094] Z = log(fft(G(t))) (14),

[0095] wl=ifft(Z)*W (15),

[0096] S(ω)=exp(fft(wl)); (16),

[0097] In equations (14), (15), and (16), log represents logarithm, fft represents Fourier transform, ifft represents inverse Fourier transform, exp represents exponential, and W represents the time window range.

[0098] Step 6, multiply the obtained angle-frequency domain reflection coefficient and the wavelet in the frequency domain to obtain the angle-frequency domain synthetic seismic record, and perform inverse Fourier transform to the angle-time domain to obtain the angle-time domain synthetic seismic record:

[0099]

[0100] In formula (17), d(θ,t) is the angle-time domain synthetic seismic record, S(ω) represents the frequency domain wavelet, and R pp (θ, ω) is the reflection coefficient in the angle-frequency domain obtained by equation (12), e is a natural constant, ω is the angular frequency, i is an imaginary unit, and π is pi. The entire forward modeling process is now completed.

[0101] A block of seismic data in a work area and well data in a corresponding area are selected, and the method of the present invention is used to perform seismic data stratification and complete forward modeling. Figure 2 The P-wave velocity, S-wave velocity and density in the well logging curve are displayed. Figure 3 It represents the result of calculating the peaks and troughs of a single-channel seismic record. The hollow points represent the peaks and the solid points represent the troughs. It can be found that this method accurately calculates the position of each peak and trough, and can clearly distinguish the position of the stratum interface. According to the corresponding positions of the peaks and troughs, the corresponding data of the well logging are substituted to obtain the stratified results of velocity and density, such as Figure 4 As shown in the figure, this process can adjust the fineness of the stratification based on the prior knowledge of geological data. This figure is only divided according to the results of the actual seismic record stratification without any adjustment. Figure 3 The single-channel seismic record shown can be used to extract seismic wavelets, and the results are as follows Figure 5 To verify the accuracy of this method, the synthetic records of traditional logging layers were compared with this method, and the results were obtained. Figure 6a , Figure 6b , Figure 6c By comparing the residuals, we can find that the accuracy of synthesizing a single gather (within 30°) by this method is almost the same as that of the traditional method. However, by recording the operation time of the traditional method and this method (as shown in Table 1), this method has increased the speed by about 45 times.

[0102] Table 1 Comparison of operation time between traditional method and this method

[0103] name Traditional methods This method Operation time(S) 4.5716 0.1002

[0104] It can be seen from the above charts that the method of the present invention not only has extremely high accuracy, but also has a great improvement in speed, which proves the rationality and correctness of the method of the present invention.

[0105] The above embodiments are only preferred embodiments for fully illustrating the present invention, and the protection scope of the present invention is not limited thereto. Any equivalent substitution or change made by a person skilled in the art based on the present invention is within the protection scope of the present invention.

Claims

1. A method for synthesizing seismic records based on wave equation, characterized in that: include: Step 1, extracting a single-channel observed seismic record, introducing an extreme value finding function, and finding the peak and trough positions; Step 2, stratify the seismic data according to the obtained wave crest and trough positions, substitute the P-wave velocity, S-wave velocity and density in the well logging data, and obtain the average P-wave velocity, S-wave velocity and density of the corresponding stratification; Step 3, substituting the obtained P-wave velocity, S-wave velocity and density average value of the corresponding layer into the formula of the vector reflectivity method to obtain the slowness-frequency domain reflection coefficient; Step 4, resampling the obtained slowness-frequency domain reflection coefficient to obtain the angle-frequency domain reflection coefficient; Step 5, extracting frequency domain wavelets from the single-channel observed seismic record; Step 6, multiplying the obtained angle-frequency domain reflection coefficient and the extracted frequency domain wavelet to obtain an angle-frequency domain synthetic seismic record, and then converting it to the angle-time domain to obtain an angle-time domain synthetic seismic record; The step 1 comprises: extracting a single-channel observed seismic record, substituting it into a peak value obtaining function to calculate the position of the wave crest, and then calculating the maximum value of the waveform, and subtracting it from the original waveform to obtain the position of the wave trough; P = f(G(t))(2), G'(t) = G(t) - MAX(G(t))(3), L=f(G'(t)) (4), In formulas (2), (3), and (4), f represents the extreme value finding function, G(t) is the single-channel observed seismic record, P is the position of the crest, L is the position of the trough, and MAX represents the maximum value finding function; In step 2, the method for obtaining the average values ​​of the longitudinal wave velocity, the shear wave velocity and the density of the corresponding layer is: In formula (5), represents the average parameter value of the nth layer, N represents the position of the nth extreme point, M represents the position of the n+1th extreme point, and r i represents any parameter value within this range, and the parameters include longitudinal wave velocity, shear wave velocity and density; In step 5, the method of extracting frequency domain wavelets from single-channel observed seismic records is as follows: S(ω)=waveletextract(G(t)) (13), In formula (13), S(ω) represents the frequency domain wavelet, G(t) represents the single-channel observed seismic record, and waveletextract represents the wavelet extraction function; The wavelet extraction function is as follows: Z = log(fft(G(t))) (14), wl=ifft(Z)*W (15), S(ω)=exp(fft(wl)); (16), In formulas (14), (15), and (16), log represents logarithm, fft represents Fourier transform, ifft represents inverse Fourier transform, exp represents exponential, and W represents the time window range.

2. A method for synthesizing seismic records based on wave equation as claimed in claim 1, characterized in that: In step 3, the method for obtaining the slowness-frequency domain reflection coefficient is: In formula (6), R pp (p,ω) represents the slowness-frequency domain reflection coefficient, p is the horizontal slowness, ω is the angular frequency, υ0 is a six-element vector, υ0(4) and υ0(1) represent the 4th and 1st elements of υ0 respectively; υ0 is calculated as: υ0=Q0Q1···,Q N-1 u N (7), u n =Q n u n+1 , (8), u N =[1 0 0 0 0 0] T . (10), In equations (6), (7), (8), (10), and (11), υ0 represents the total reflection response of the top layer of the medium, n represents any layer of the medium, N represents the total number of layers of the medium, and υ n+1 represents the response of the n+1th underground layer, υ N represents the initial response of the lowest layer of the medium, Q n represents the transfer matrix, and represents the influence of the nth layer of medium on the amplitude of seismic waves, The superscript + indicates the lower interface of the layer interface. The superscript - indicates the upper interface of the layer interface. and Each contains 16 independent components; the matrix E n Represents the change in phase of a seismic wave when it passes through the nth layer of medium.

3. A method for synthesizing seismic records based on wave equation as claimed in claim 1, characterized in that: The step 4 comprises: Introducing p n = sinθ / α n , where p n represents the horizontal slowness of the nth layer, θ represents the incident angle, α n represents the longitudinal wave velocity of the nth layer; Resample the slowness-frequency domain reflection coefficient to obtain the angle-frequency domain reflection coefficient: R pp (θ,ω)=resample(R pp (p,ω)) (12), In formula (12), resample represents the resampling process, R pp (θ, ω) represents the angle-frequency domain reflection coefficient.

4. A method for synthesizing seismic records based on wave equation as claimed in claim 1, characterized in that: The step 6 comprises: multiplying the obtained angle-frequency domain reflection coefficient and the wavelet in the frequency domain to obtain an angle-frequency domain synthetic seismic record, and performing an inverse Fourier transform to the angle-time domain to obtain an angle-time domain synthetic seismic record: In formula (17), d(θ,t) is the angle-time domain synthetic seismic record, S(ω) represents the frequency domain wavelet, and R pp (θ, ω) is the angle-frequency domain reflection coefficient, and ω is the angular frequency.

Images