A method and device for generating an initial velocity model

Abstract

An embodiment of this specification discloses a method and device for generating an initial velocity model. By obtaining input data d input (i) and input parameter σ; performing smoothing processing on the input data using the following formula to obtain the smoothed input data d output (i); using a preset objective function to perform iterative processing according to the smoothed input data d output (i) to generate a target initial velocity model, thereby realizing the extraction of large-scale features of seismic signals using the dilated convolution method when there is a lack of an accurate initial velocity model and insufficient low-frequency components of seismic signals, and then using this feature to construct long-wavelength information of the target area to provide a relatively accurate initial velocity model for full waveform inversion.

Description

Technical Field

[0001] This specification relates to the field of exploration geophysics, and particularly to a method and device for generating an initial velocity model. Background Art

[0002] With the development of oil geophysical exploration in China, the exploration focus has shifted to deep and complex structural oil and gas reservoirs. Such oil and gas reservoirs have the following characteristics: (1) The reservoir burial depth is large, with a burial depth of 6000 - 8000 m or even more than 10000 m; (2) The structure is complex, such as reverse thrust and stepped faults, etc.; (3) The seismic signal acquisition is difficult, and low-frequency information is missing.

[0003] In theory, the full waveform inversion method can effectively construct the velocity information of this type of target structural area. The full waveform inversion requires a relatively accurate initial velocity information as the start of the inversion to prevent the occurrence of cycle skipping phenomenon, which may cause the inversion to fall into local extrema and lead to errors in the full waveform inversion. In actual signal acquisition, due to the interference of uncertain factors such as acquisition cost and environmental noise, the low-frequency components are missing in the seismic signal, and it is difficult to construct an accurate initial velocity field due to the lack of low-frequency components.

[0004] Based on this, a solution that can generate an accurate initial velocity model is needed. Summary of the Invention

[0005] The purpose of the present invention is to provide a solution that can generate an accurate initial velocity model.

[0006] To solve the above technical problems, the present invention adopts the following technical solutions:

[0007] In a first aspect, a method for generating an initial velocity model is provided, including: obtaining input data d input (i) and input parameters, where the input data includes observed data d obs , simulated records and an initial velocity field, and the input parameters include spatial sampling intervals dx and dz, time sampling interval dt, number of time sampling points n t , dominant frequency f0, number of geophones N r , number of PML boundary layers S, convolution rate r, and control factor σ; smoothing the input data using the following formula to obtain smoothed input data d output (i): where i and j are different times of shot records, and g[i, j] is a filter; using a preset objective function, performing iterative processing according to the smoothed input data d output (i) to generate a target initial velocity model.

[0008] In a second aspect, a computer device is provided, including a memory, a processor, and a computer program stored on the memory and executable on the processor. When the processor executes the program, the method described in the first aspect is implemented.

[0009] At least one of the above technical solutions adopted in the embodiments of this specification can achieve the following beneficial effects: By obtaining input data d input (i) and input parameters, where the input data includes observed data d obs , simulation records, and an initial velocity field, and the input parameters include spatial sampling intervals dx and dz, time sampling interval dt, number of time sampling points n t , main frequency f0, number of geophones N r , number of PML boundary layers S, convolution rate r, and control factor σ; the input data is smoothed using the following formula to obtain smoothed input data d output (i); using a preset objective function, iterative processing is performed based on the smoothed input data d output (i) to generate a target initial velocity model, thereby realizing the extraction of large-scale features of seismic signals using the dilated convolution method when an accurate initial velocity model is missing and the low-frequency components of seismic signals are insufficient, and then using this feature to construct the long-wavelength information of the target area, providing a relatively accurate initial velocity model for full waveform inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0010] Figure 1a is a schematic flowchart of a method for generating an initial velocity model provided by an embodiment of this specification;

[0011] Figure 1b is a flowchart of an embodiment of the present invention;

[0012] Figure 2 is the true velocity field of an embodiment of the present invention;

[0013] Figure 3 is the initial velocity field of an embodiment of the present invention;

[0014] Figure 4 is the original observed record of an embodiment of the present invention;

[0015] Figure 5 is the processed observed record of an embodiment of the present invention;

[0016] Figure 6 is the spectrogram of the observed record before and after processing of an embodiment of the present invention;

[0017] Figure 7 is the result of conventional full waveform inversion;

[0018] Figure 8 The inversion background velocity result of an embodiment of the present invention;

[0019] Figure 9 The final velocity inversion result of an embodiment of the present invention;

[0020] Figure 10 The single-channel velocity comparison between an embodiment of the present invention and the conventional full-waveform inversion result. Specific embodiments

[0021] To make the objectives, technical solutions and advantages of the present application clearer, the technical solutions of the present application will be clearly and completely described below in conjunction with specific embodiments of the present application and the corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all of the embodiments. Based on the embodiments in this specification, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the scope of protection of the present application.

[0022] In the first aspect, as Figure 1a shown, Figure 1a is a schematic flowchart of a method for generating an initial velocity model provided by an embodiment of this specification, including:

[0023] S101, obtaining input data d input (i) and input parameters.

[0024] The input data includes observed data d obs , simulated records and an initial velocity field, and the input parameters include spatial sampling intervals dx and dz, time sampling interval dt, number of time sampling points n t , main frequency f0, number of geophones N r , number of PML boundary layers S, convolution rate r, and control factor σ. The initial velocity field is modeled based on known geological information and logging information, providing a general understanding of the regional velocity.

[0025] S103, smoothing the input data using the following formula to obtain the smoothed input data d output (i).

[0026] The specific calculation method is:

[0027]

[0028] where i and j are different times of shot records, and g[i,j] is the filter. i and j are different times of shot records. r is the dilation convolution rate, and its values are 3, 5, 7, ···. The larger r is, the lower the low-frequency features obtained, but at the same time, the larger r will also cause the loss of effective information. Therefore, it is more appropriate to use r = 3.

[0029] In the foregoing formula (1), for d input (i) Perform dilated convolution processing to eliminate high-frequency oscillations in the shot record and highlight the features of the shot record on a large scale. Considering that the seismic signal is a continuous signal, there is a strong correlation between the vibration information at the target moment to be obtained and the vibration information at the surrounding moments of this point. Therefore, when introducing the filter g[i, j], a Gaussian filter is adopted. Among them, n t is the number of sampling points, σ is the variance of the Gaussian kernel function, and different frequency band ranges of seismic signals are obtained by controlling the value of the factor σ. Therefore, it can be known that by controlling the value of σ, seismic signals of different scales can be obtained. The low-frequency characteristic signals can be used to restore the background velocity information, while the high-frequency signals can depict the fine structure of the velocity model.

[0030] S105. Using a preset objective function, according to the smoothed input data d output (i) Perform iterative processing to generate an initial target velocity model.

[0031] The specific form of the objective function can be set by itself based on actual needs. For example, a feasible specific form of the objective function can be the L2 objective function Among them, E(v) represents the residual of the data, where d' obs is the smoothed observed data, d' mod is the smoothed simulated record, and T is the matrix transpose. The simulated record d mod can be obtained by using finite-difference forward simulation.

[0032] The two-dimensional acoustic wave equation can be expressed as:

[0033]

[0034]

[0035] Among them, u is the wave field in the time domain, v is the seismic wave velocity, and f is the source term. The filtering part in the foregoing formula (1) is obtained by using second-order time and 2N-order space finite differences for formula (2). Among them, a i is the difference coefficient, k represents the time value, m and n respectively represent the spatial positions of the wave field, dt is the time sampling interval, and N is the difference order.

[0036] Introduce the processed observed record and simulated record to construct the L2 objective function:

[0037]

[0038] Among them, d' obs is the observed data after dilated convolution processing, and d' mod is the simulated record after dilated convolution processing.

[0039] Furthermore, the velocity model iteration can be performed according to the L2 objective function to obtain a high-precision initial velocity model for full waveform inversion. Using the adjoint state method, the inversion gradient of the objective function based on the L2 norm can be obtained as

[0040] where is the wavefield extrapolation operator, S -1 represents the backward propagation process of the wavefield, R is the operator that restricts the wavefield to the geophone positions, and δd = d obs - d mod is the difference between the observed data and the simulated data, which is the adjoint source for inversion.

[0041] Then, the model is updated using formula (8):

[0042]

[0043] where v n-1 represents the model velocity information of the (n - 1)-th update, v n represents the model velocity information of the n-th update, and α represents the update step size. Then, the model is updated to obtain a better initial velocity model of the target.

[0044] In the specific iteration process, the iteration termination condition can also be set based on the control factor σ, that is, if the current control factor σ is less than the preset value, the iteration terminates. Among them, in each round of iteration, if the accuracy of the generated velocity model does not exceed the preset accuracy, the control factor σ is reduced based on the preset value in the current iteration.

[0045] For example, the following steps can be used for iteration:

[0046] S41. The initial control factor σ = 64 (the initial value can be set by yourself according to actual needs). After processing the observed records and the simulated records, inversion is performed. It is judged whether the accuracy of the currently obtained initial velocity model meets the requirements. If it meets, the control factor σ is reduced by half of the control factor of the current stage. The inversion result of the current stage is used as the initial model for inversion with the reduced control factor σ, and inversion is performed. If not, the inversion of the current step continues.

[0047] S42. When the control factor σ is reduced multiple times, it is judged whether the control factor σ is less than 2. If so, the inversion of the current stage is stopped and the inversion result is output. If not, step S42 is repeated.

[0048] After obtaining the target initial velocity model, the obtained high-precision target initial velocity model can also be used for full waveform inversion to finally obtain a high-precision full waveform inversion result. The inversion result obtained in S42 is used as the input for full waveform inversion to perform full waveform inversion.

[0049] By obtaining input data d input (i) and input parameters, where the input data includes observed data d obs , simulated records, and an initial velocity field, and the input parameters include spatial sampling intervals dx and dz, a time sampling interval dt, the number of time sampling points n t , dominant frequency f0, the number of geophones N r , the number of PML boundary layers S, convolution rate r, and control factor σ; the following formula is used to smooth the input data to obtain the smoothed input data d output (i); using a preset objective function, iterative processing is performed based on the smoothed input data d output (i) to generate a target initial velocity model, so as to realize extracting the large-scale features of seismic signals using the dilated convolution method when there is no accurate initial velocity model and the low-frequency components of seismic signals are insufficient, and then using this feature to construct the long-wavelength information of the target area to provide a relatively accurate initial velocity model for full waveform inversion.

[0050] To more specifically illustrate the method of the present invention, the method of the present invention is illustrated by taking the Marmousi model as an example. The specific process is as Figure 1b shown. The true model (as Figure 2 shown) has 596 horizontal grid points and 180 vertical grid points, and the spatial interval of grid points is 8 meters. The model contains 3 sets of complex structures such as faults. The input observed records (as Figure 4 shown), the initial velocity field (as Figure 3 shown), and the initial velocity changes linearly from shallow to deep, which is only the general change trend of the velocity and does not contain structural information.

[0051] Set the relevant parameters for the inversion of the method of the present invention, and set the initial control factor σ to 64. The distribution of the observation system is as follows: at a depth of 24 meters from the surface, 60 shots are evenly distributed horizontally, with an interval of 80 meters between each shot. Each shot is received by 596 geophone points, and the interval between geophone points is 8 meters. The inversion parameters are as follows: use a 20 Hz Ricker wavelet as the source, the time sampling interval is 0.5 ms, the number of time sampling points is 4500, and the recording duration is 2.25 s. After each stage of inversion, determine whether the current inversion result meets the accuracy requirements. If the judgment condition is met, reduce the value of the control factor σ, and use the inversion result of the current stage as the initial model for the inversion of the next control factor σ. If the judgment condition is not met, continue the inversion. When the inversion result meets the accuracy requirements of the initial velocity model of full waveform inversion, stop the inversion, and use this inversion result as the initial velocity model of conventional full waveform inversion. Finally, obtain a high-precision full waveform inversion result.

[0052] Figure 4 It is a raw single-shot observation record, where the energy decreases with the increase of depth, and the reflected wave waveform is relatively complex.

[0053] Figure 5 It is the single-shot record processed by the method of the present invention, and this method extracts the large-scale features in the original observation record.

[0054] Figure 6 It is the spectrogram of the observation record before and after processing. It can be seen from the spectrogram that as the value of the control factor σ gradually increases, the frequency of the shot record gradually decreases. The large-scale features of the signal correspond to the low-frequency information of the signal. For seismic wave velocity inversion, the low-frequency components can be used to restore the background velocity information of the geological model. Prevent full waveform inversion from falling into local extrema due to inaccurate background velocity, thereby reducing inversion errors and improving inversion accuracy.

[0055] Figure 7 It is for Figure 3 As the conventional full waveform inversion result with the initial model, it can be seen from the figure. Under the condition of low initial model accuracy, errors occur in the inversion, velocity anomalies appear in the shallow layer, and the velocity information in the deep layer cannot be inverted.

[0056] Figure 8 It is for the method of the present invention with Figure 3 As the initial model inversion result, compared with the linear velocity field (as shown in Figure 3 ), this result restores the changing trend of the background velocity information of the true model (as shown in Figure 2 ), and has a certain degree of restoration of the structural information in the model.

[0057] Figure 9 It is the final velocity inversion result of an embodiment 2 of the present invention; compared with Figure 7For the conventional full waveform inversion results, the method of the present invention can effectively recover the shallow-deep velocity information in the model. Faults at the near surface and in the middle and deep layers are clearer, and the main structural information of the Marmousi model is restored.

[0058] Figure 10 For the single-channel velocity information at Distance = 3000 meters of the true velocity model, the initial velocity model, the conventional full-wave inversion result, and the inversion result of the method of the present invention. Compared with the conventional full waveform inversion result, the inversion result of the method of the present invention is more similar to the inversion result of the true velocity model and can better reflect the change of the true model velocity.

[0059] In a second aspect, correspondingly, the embodiments of the present application further provide a computer device, including a memory, a processor, and a computer program stored on the memory and executable on the processor. Wherein, when the processor executes the program, it implements the foregoing imaging method based on wave field decomposition.

[0060] Each embodiment in this specification is described in a progressive manner. The same or similar parts between the embodiments can be referred to each other. Each embodiment focuses on the differences from other embodiments. In particular, for the embodiments of devices, equipment, and media, since they are basically similar to the method embodiments, the description is relatively simple. The relevant parts can be referred to the partial description of the method embodiments, and will not be elaborated here one by one.

[0061] Each embodiment in this specification is described in a progressive manner. The same or similar parts between the embodiments can be referred to each other. Each embodiment focuses on the differences from other embodiments. In particular, for the embodiments of devices, equipment, and media, since they are basically similar to the method embodiments, the description is relatively simple. The relevant parts can be referred to the partial description of the method embodiments, and will not be elaborated here one by one.

[0062] The above describes specific embodiments of this specification. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps or modules recited in the claims can be executed in a different order than in the embodiments and still achieve the desired result. Additionally, the processes depicted in the figures do not necessarily require the specific order or sequential order shown to achieve the desired result. In certain embodiments, multitasking and parallel processing are also possible or may be advantageous.

Claims

1. A method for generating an initial velocity model, comprising: Obtain the input data d input (i) and input parameters, where the input data includes the observed data d obs , simulation records, and the initial velocity field, and the input parameters include the spatial sampling intervals dx and dz, the time sampling interval dt, the number of time sampling points n t , the dominant frequency f0, the number of geophones N r , the number of PML boundary layers S, the convolution rate r, and the control factor σ; The input data is smoothed using the following formula to obtain the smoothed input data d output (i): Among them, i and j record different moments, and g[i, j] is the filter; Using a preset objective function, based on the smoothed input data d output (i) Perform iterative processing to generate an initial target velocity model; Among them, the preset objective function includes: L2 objective function where d' obs is the smoothed observed data, d' mod is the smoothed simulation record, and T is the matrix transpose.

2. The method according to claim 1, wherein Based on the smoothed input data d output (i) Perform iterative processing to generate a target initial velocity model, including: Determine the inversion gradient of the L2 objective function where is the wavefield continuation operator, S -1 represents the backpropagation process of the wavefield, R is the operator that restricts the wavefield to the geophone positions, and δd = d obs - d mod , which is the difference between the observed data d obs and the simulated data d mod . Iterative update is performed according to the inversion gradient. When the iterative termination condition is reached, a target initial velocity model is generated: v n = v n-1 - α▽ v E(v), where v n-1 represents the model velocity information updated at the (n - 1)-th time, v n represents the model velocity information updated at the n-th time, and α represents the update step size.

3. The method according to claim 2, wherein, The iteration termination condition includes: Whether the current control factor σ is less than a preset value. Among them, in each round of iteration, if the accuracy of the generated velocity model does not exceed the preset accuracy, the control factor σ is reduced based on the preset value in the current iteration.

4. The method according to claim 1, further comprising: Performing full waveform inversion according to the target initial velocity model to generate a full waveform inversion imaging result.

5. A computer device, comprising a memory, a processor, and a computer program stored on the memory and executable on the processor, wherein, When the processor executes the program, the method according to any one of claims 1 to 4 is implemented.

Images