A method and device for determining seismic surface wave velocity

Abstract

The present invention provides a method and device for determining seismic surface wave velocity, which belongs to the field of oil and gas geophysical exploration engineering. The method comprises: obtaining seismic observation data to be processed, extracting fundamental surface wave data from the seismic observation data, performing time-frequency conversion on the fundamental surface wave data to obtain frequency domain surface wave data; determining the effective frequency band according to the frequency corresponding to the maximum value of the spectrum energy in the frequency domain surface wave data, reconstructing the wave field of the frequency domain surface wave data within the effective frequency band, and obtaining the amplitude field and phase spectrum of the frequency domain surface wave data after the wave field reconstruction within the effective frequency band; processing the phase spectrum using the phase unwrapping technology to obtain the travel time field; using the amplitude field, the gradient of the amplitude field and the gradient of the travel time field to obtain the seismic surface wave velocity within the effective frequency band. The present invention can be applied to any complex irregular observation system, has high surface wave velocity calculation accuracy, can directly invert the surface wave data of any multiple earthquake sources, and has high surface wave velocity calculation efficiency.

Description

Technical Field

[0001] The invention relates to a method and a device for determining seismic surface wave velocity, and belongs to the field of oil and gas geophysical exploration engineering. Background Art

[0002] Seismic surface waves are the main geophysical properties of the near surface, which can reflect the characteristics of seismic wave propagation in the surface. Surface wave exploration is widely used in velocity structure detection near the surface and is the main means of engineering geophysical exploration. In addition, in the field of petroleum geophysical exploration, surface waves are generally regarded as strong noise, which requires special removal methods, and almost all surface wave removal methods rely on the velocity of surface waves. Therefore, the acquisition of seismic surface wave velocity is very critical. At present, there are two main commonly used surface wave velocity estimation methods:

[0003] The first method is to directly calculate the apparent velocity through the slope on the seismic profile, and then determine the maximum and minimum values ​​of the apparent velocity of the surface wave, such as the method disclosed in the invention patent document with announcement number CN106569281B, which uses an interactive velocity picking tool to calculate the apparent velocity range of the surface wave, groups the apparent velocity range to obtain multiple apparent velocity cones, processes each apparent velocity cone to obtain the starting position of the current surface wave, and performs subsequent processing based on the starting sample point of the current surface wave to remove the surface wave noise. However, this method can only reflect the overall situation of the work area, and the surface wave velocity obtained is only an approximate range, which is not accurate enough.

[0004] The second method is to extract the dispersion curve of the surface wave through spectrum calculation, and then estimate the surface wave velocity, such as the method disclosed in the invention patent document with announcement number CN102759749B, which obtains the fk spectrum of Rayleigh surface wave through two-dimensional Fourier transform, and converts the fk spectrum into depth-velocity (hv) spectrum according to the relationship between the penetration depth, velocity and frequency, wave number, wavelength, and Poisson's ratio of Rayleigh surface wave, and realizes Rayleigh surface wave velocity analysis according to the hv spectrum. However, this method is highly dependent on the distribution of detection points, and requires the observation system to be distributed as regularly as possible. However, in the actual collection process of seismic observation data, due to the existence of villages, rivers, industrial production, living facilities, etc., the distribution of the observation system is very irregular. In this case, the accuracy of the surface wave velocity obtained by this method is limited. In addition, this method can only process the data of one source at a time, and cannot process the data of multiple sources, and the calculation efficiency is low.

[0005] In summary, when calculating the surface wave velocity using the current surface wave velocity estimation method, there are problems of low calculation result accuracy and low calculation efficiency. Summary of the invention

[0006] The object of the present invention is to provide a method and device for determining seismic surface wave velocity, so as to solve the problem of low calculation accuracy and low calculation efficiency when calculating surface wave velocity using the current surface wave velocity estimation method.

[0007] To achieve the above object, the present invention provides a method for determining seismic surface wave velocity, the method comprising the following steps:

[0008] (1) obtaining seismic observation data to be processed, extracting fundamental surface roll data from the seismic observation data, and performing time-frequency conversion on the fundamental surface roll data to obtain frequency domain surface roll data;

[0009] (2) determining an effective frequency band according to a frequency corresponding to a maximum value of spectral energy in the frequency-domain surface roll data, reconstructing a wave field of the frequency-domain surface roll data within the effective frequency band, and obtaining an amplitude field and a phase spectrum of the frequency-domain surface roll data after the wave field reconstruction within the effective frequency band;

[0010] (3) using phase unwrapping technology to process the phase spectrum to obtain a travel time field;

[0011] (4) Using the amplitude field, the gradient of the amplitude field and the gradient of the travel time field, the seismic surface wave velocity within the effective frequency band is obtained.

[0012] The present invention also provides a device for determining seismic surface wave velocity, which includes a processor and a memory, wherein the processor executes a computer program stored in the memory to implement the above-mentioned method for determining seismic surface wave velocity.

[0013] The beneficial effects of the present invention are as follows: (1) firstly, starting from the wave field reconstruction of the frequency domain surface wave data, the problem of irregular observation system is solved, and it is ensured that a high surface wave velocity calculation accuracy can be obtained when the observation system is irregularly distributed. The method has strong adaptability, so that the present invention is not only applicable to regular observation system, but also applicable to any complex irregular observation system, which is closer to practical application, can solve actual geophysical problems, and can provide a more accurate velocity structure for subsequent suppression or utilization; (2) a stable travel time field is obtained by phase unwrapping technology, and then the seismic surface wave velocity within the effective frequency band is obtained by using the amplitude field, the gradient of the amplitude field and the gradient of the travel time field. Since there is no point source assumption when calculating the seismic surface wave velocity in the present invention, the present invention can process the data of line source and surface source, and the line source and surface source are equivalent to the combination of multiple point sources. That is to say, the present invention can directly invert the surface wave data of any multiple sources, avoiding the conventional tomography that can only invert the data of one source at a time, thereby greatly improving the efficiency of surface wave velocity calculation. The present invention can provide a fast and accurate near-surface velocity structure.

[0014] Furthermore, in the above method and device, the seismic surface wave velocity within the effective frequency band is obtained using a wave field gradient formula, and the wave field gradient formula is as follows:

[0015]

[0016] In the formula, is the seismic surface wave velocity, It's a walking time field. is the gradient of the field during travel, is the amplitude field, is the square of the gradient of the amplitude field, ω is the frequency, are the coordinates of the detection point.

[0017] Furthermore, in the above method and device, the frequency domain surface wave data after wave field reconstruction is obtained by the following steps: First, by formula Get sparse transform coefficients Then, using the formula Obtain the frequency domain surface wave data after the wave field reconstruction Where m is the sparse transformation coefficient of the frequency domain surface wave data, Ψ is the sampling matrix of the observation system, where there are detection points, it is set to 1, and where there are no detection points, it is set to 0, Φ is the sparse transformation matrix, is the frequency domain surface wave data, and λ is the regularization parameter.

[0018] Furthermore, in the above method and device, the effective frequency band is [fmain-2 / nt / dt, fmain+2 / nt / dt], fmain is the frequency corresponding to the maximum value of the spectral energy in the frequency domain surface roll data, nt is the number of detection points, and dt is the sampling interval.

[0019] Furthermore, in the above method and device, the fundamental-order surface roll data is extracted by using a gun domain Radon transform suppression method or a cross-domain surface roll estimation method. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 is a flow chart of a method for determining seismic surface wave velocity in an embodiment of the method of the present invention;

[0021] Figure 2 A schematic diagram of the arrangement of the observation system in the embodiment of the method of the present invention;

[0022] Figure 3 A schematic diagram of frequency domain surface wave data in an embodiment of the method of the present invention;

[0023] Figure 4 A schematic diagram of frequency domain surface roll data under a regular grid restored in an embodiment of the method of the present invention;

[0024] Figure 5 is a schematic diagram of a phase spectrum in an embodiment of the method of the present invention;

[0025] Figure 6 Schematic diagram of the travel time field in the embodiment of the method of the present invention;

[0026] Figure 7-1 is a schematic diagram of the surface wave velocity field calculated using the method of this embodiment;

[0027] Figure 7-2 is a schematic diagram of the error field of the surface wave velocity calculated using the method of this embodiment;

[0028] Figure 8-1 This is a schematic diagram of the surface wave velocity field calculated using conventional algorithms;

[0029] Figure 8-2 Schematic diagram of the error field of surface wave velocity calculated using conventional algorithms;

[0030] Fig. 9 It is a schematic diagram of the structure of the seismic surface wave velocity determination device in the device embodiment of the present invention. DETAILED DESCRIPTION

[0031] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments.

[0032] Method Example:

[0033] like Figure 1 As shown, the method for determining seismic surface wave velocity in this embodiment includes the following steps:

[0034] (1) Obtaining earthquake observation data to be processed;

[0035] In this embodiment, the earthquake observation data is Figure 2 The observation system shown collects, Figure 2 The points in the figure represent detection points. It can be seen that the arrangement of detection points in this observation system is very irregular, so this observation system belongs to an irregular observation system.

[0036] (2) Extracting fundamental surface wave data from seismic observation data is the coordinate of the detection point, t is the round-trip travel time;

[0037] Among them, conventional surface wave extraction methods are used to extract fundamental surface wave data, such as the Radon transform suppression method in the gun domain or the cross-domain surface wave estimation method. Both methods are existing technologies. For details, please refer to the following literature: Su Yue, Yang Zhentao, etc. published a paper entitled "Research on the Application of Active Source Rayleigh Wave Multi-mode Extraction Method Based on Vector Wavenumber Transform in Near-surface Stratum Structure Detection" in the Journal of Peking University (Natural Science Edition) Volume 56, Issue 3 in 2020; Xue Chao published a paper entitled "Research and Application of Cross Arrangement Denoising Method" in Inner Mongolia Petrochemical Industry Volume 45, Issue 7 in 2019.

[0038] (3) Fundamental surface wave data Perform time-frequency conversion (such as Fourier transform) to obtain frequency domain surface wave data ω is the frequency;

[0039] Frequency domain surface wave data The wave field diagram of the real part is as follows Figure 3 As shown, Figure 3 The color bar on the right represents the amplitude energy. Figure 3 The blank area in the middle indicates that there is no detector at that point, that is, Figure 2 The irregular observing system shown does not receive data at this location.

[0040] (4) Based on frequency domain surface wave data The frequency fmain corresponding to the maximum value of the spectrum energy in determines the effective frequency band, which is [fmain-2 / nt / dt, fmain+2 / nt / dt], where nt is the number of sampling points (i.e. the number of detection points) and dt is the sampling interval;

[0041] (5) Frequency domain surface wave data within the effective frequency band Perform wave field reconstruction to obtain frequency domain surface wave data after wave field reconstruction (hereinafter referred to as frequency domain surface wave data under regular grid); the process of wave field reconstruction is as follows:

[0042] First, use SPGL1 (i.e., spectral projection gradient method) or CVX (i.e., convex optimization toolbox method) to solve the formula Get sparse transform coefficients Where m is the sparse transformation coefficient of the frequency domain surface wave data; Ψ is the sampling matrix of the observation system, which is set to 1 where there is a detection point and 0 where there is no detection point; Φ is the sparse transformation matrix, which can be specifically the sparse transformation matrix in DCT discrete cosine transform or curvelet transform; λ is the regularization parameter, which can be obtained using the L curve method; ‖‖2 represents the L2 norm, and ‖‖1 represents the L1 norm;

[0043] Then, using the formula Get the frequency domain surface wave data under regular grid

[0044] Frequency domain surface wave data under regular grid in this embodiment See Figure 4 , Figure 4 The color bar on the right represents the amplitude energy.

[0045] (6) Obtaining frequency domain surface wave data on a regular grid within the effective frequency band The amplitude field and phase spectrum

[0046] Among them, the amplitude field Phase spectrum express The imaginary part of express The real part of .

[0047] The phase spectrum in this embodiment like Figure 5 As shown, the horizontal axis represents east (Km), the vertical axis represents north (Km), and the color depth in the figure represents the size of the phase value, that is, the size of the take-off time.

[0048] (7) Use phase unwrapping technology to analyze the phase spectrum Processing to obtain travel time field The specific process is as follows:

[0049] First, Then from Starting from the second point, calculate the difference between the current point and the previous point. If the difference is greater than pi / 2, then Subtract pi from the current point and all subsequent points. If the difference is less than -pi / 2, then Add pi to the current point and all subsequent points, and repeat the above steps until the last point. Here, pi represents the radian value of π.

[0050] The travel time field in this embodiment like Figure 6 As shown, the horizontal axis represents east (km), the vertical axis represents north (km), and the light and dark colors in the figure represent the time, that is, the length of travel time.

[0051] (8) The wave field gradient formula is used to obtain the seismic surface wave velocity within the effective frequency band.

[0052] The wave field gradient formula is as follows:

[0053]

[0054] In the formula, is the seismic surface wave velocity, It's a walking time field. is the gradient of the field during travel, is the amplitude field, is the square of the gradient of the amplitude field, ω is the frequency, are the coordinates of the detection point.

[0055] In this embodiment, the seismic surface wave velocity field and the surface wave velocity error field obtained by using the wave field gradient formula are respectively as follows: Figure 7-1 and Figure 7-2 As shown in the figure, the seismic surface wave velocity field and the error field of surface wave velocity calculated by conventional algorithms are shown in Figure 8-1 and Figure 8-2 As shown, Figure 7-1 and Figure 8-1 The color bar on the right side of the figure represents the surface wave velocity value. Figure 7-2 and Figure 8-2 The color bar on the right side of the center represents the surface wave velocity error value, and the darker the color, the larger the error.

[0056] By comparison Figure 7-2 and Figure 8-2 It can be seen that Figure 8-2 The darker the color, the greater the error, which indicates that the surface wave velocity obtained by the method of this embodiment is more accurate and has higher precision.

[0057] The derivation process of the wave field gradient formula of this embodiment is as follows:

[0058] The wave equation formula is:

[0059]

[0060] Through Fourier transform, the frequency domain form of the wave equation can be obtained:

[0061]

[0062] Will Written in plane wave form, we can get U(r,w)=A(r)e iw(t-τ(r)) , substitute into the above formula, and assume is smooth, and can be obtained

[0063] In the above formula, is the surface wave velocity, and They are the time domain and frequency domain forms of surface waves respectively.

[0064] The method for determining the seismic surface wave velocity of the present embodiment, on the one hand, starts with the reconstruction of the surface wave field in the frequency domain, solves the problem of irregular observation systems, ensures that a high surface wave velocity calculation accuracy can be obtained when the observation system is irregularly distributed, and has strong adaptability, so that the method of the present embodiment can be applied not only to regular observation systems, but also to any complex irregular observation systems, is closer to practical applications, can solve actual geophysical problems, and can provide a more accurate velocity structure for subsequent suppression or utilization; on the other hand, a stable travel time field is obtained by phase unwrapping technology, and then the seismic surface wave velocity within the effective frequency band is obtained by using the wave field gradient formula. Since there is no point source assumption when calculating the seismic surface wave velocity in the method of the present embodiment, the method of the present embodiment can process the data of line sources and surface sources, and the line source and surface source are equivalent to the combination of multiple point sources, that is, the method of the present embodiment can directly invert the surface wave data of any multiple sources, avoiding the conventional tomography that can only invert the data of one source at a time, thereby greatly improving the efficiency of surface wave velocity calculation. Therefore, the method for determining the seismic surface wave velocity of the present embodiment can provide a fast and accurate near-surface velocity structure.

[0065] In this embodiment, the wave field gradient formula is used to obtain the seismic surface wave velocity within the effective frequency band. This is only a specific implementation method for calculating the seismic surface wave velocity and does not constitute a limitation on the protection scope of the present invention. In practical applications, as long as the seismic surface wave velocity within the effective frequency band is calculated using the amplitude field, the gradient of the amplitude field and the gradient of the travel time field of the frequency domain surface wave data within the effective frequency band, it is within the protection scope of the present invention. For example, a new wave field gradient formula is obtained by simply deforming the wave field gradient formula of this embodiment (such as square root, etc.), and then the seismic surface wave velocity is obtained using the new wave field gradient formula, which is also within the protection scope of the present invention.

[0066] Device Example:

[0067] The seismic surface wave velocity determination device of this embodiment is as follows: Fig. 9 As shown, the device includes a processor and a memory, wherein a computer program that can be run on the processor is stored in the memory, and the processor implements the method in the above method embodiment when executing the computer program.

[0068] That is, the method in the above method embodiment should be understood as a flow of a method for determining seismic surface wave velocity that can be implemented by computer program instructions. These computer program instructions can be provided to a processor so that the processor executes these instructions to generate functions specified in the above method flow.

[0069] The processor referred to in this embodiment refers to a processing device such as a microprocessor MCU or a programmable logic device FPGA.

[0070] The memory referred to in this embodiment includes a physical device for storing information, which is usually to digitize the information and then store it in a medium using electrical, magnetic or optical methods. For example: various memories that use electrical energy to store information, such as RAM, ROM, etc.; various memories that use magnetic energy to store information, such as hard disks, floppy disks, tapes, magnetic core memories, bubble memories, and USB flash drives; various memories that use optical methods to store information, such as CDs or DVDs. Of course, there are other types of memories, such as quantum memories, graphene memories, and so on.

[0071] The device is composed of the above-mentioned memory, processor and computer program, and is implemented by the processor executing corresponding program instructions in the computer. The processor can be equipped with various operating systems, such as Windows operating system, Linux system, Android, iOS system, etc.

Claims

1. A method for determining seismic surface wave velocity, characterized in that: The method comprises the following steps: (1) obtaining seismic observation data to be processed, extracting fundamental surface roll data from the seismic observation data, and performing time-frequency conversion on the fundamental surface roll data to obtain frequency domain surface roll data; (2) determining an effective frequency band according to a frequency corresponding to a maximum value of spectral energy in the frequency-domain surface roll data, reconstructing a wave field of the frequency-domain surface roll data within the effective frequency band, and obtaining an amplitude field and a phase spectrum of the frequency-domain surface roll data after the wave field reconstruction within the effective frequency band; The effective frequency band is [fmain-2 / nt / dt, fmain+2 / nt / dt], where fmain is the frequency corresponding to the maximum value of the spectrum energy in the frequency domain surface wave data, nt is the number of detection points, and dt is the sampling interval; (3) using phase unwrapping technology to process the phase spectrum to obtain a travel time field; (4) Using the amplitude field, the gradient of the amplitude field and the gradient of the travel time field, the seismic surface wave velocity within the effective frequency band is obtained.

2. The method for determining seismic surface wave velocity according to claim 1, characterized in that: The seismic surface wave velocity within the effective frequency band is obtained using the wave field gradient formula, which is as follows: In the formula, is the seismic surface wave velocity, It's a walking time field. is the gradient of the field during travel, is the amplitude field, is the square of the gradient of the amplitude field, ω is the frequency, are the coordinates of the detection point.

3. The method for determining seismic surface wave velocity according to claim 1, characterized in that: The frequency domain surface wave data after wave field reconstruction is obtained by the following steps: First, by formula Get sparse transform coefficients Then, using the formula Obtain the frequency domain surface wave data after the wave field reconstruction Where m is the sparse transformation coefficient of the frequency domain surface wave data before wave field reconstruction, Ψ is the sampling matrix of the observation system, where there are detection points, it is set to 1, and where there are no detection points, it is set to 0, Φ is the sparse transformation matrix, is the frequency domain surface wave data before wave field reconstruction, λ is the regularization parameter, ω is the frequency, are the coordinates of the detection point.

4. The method for determining seismic surface wave velocity according to claim 1 or 2, characterized in that: The fundamental-order surface roll data are extracted using a gun domain Radon transform suppression method or a cross-domain surface roll estimation method.

5. A device for determining seismic surface wave velocity, characterized in that: The device comprises a processor and a memory, wherein the processor executes a computer program stored in the memory to implement the method for determining seismic surface wave velocity according to any one of claims 1 to 4.

Images