A waveform inversion method and apparatus based on construction shaping regularization

Abstract

This specification discloses a waveform inversion method and apparatus based on construction shaping regularization. It involves acquiring observation data and an initial velocity model; performing forward seismic wave modeling calculations using the initial model to obtain simulated seismic data; calculating the backpropagation residual wavefield using the difference between the simulated data and the observation data; and simultaneously calculating the tomographic mode waveform inversion gradient G. t and the gradient G of the offset mode waveform inversion m Using the offset mode gradient G m Determine the local structural dip angle, and determine the background velocity model m based on the local structural dip angle. n+1 According to the background velocity model m n+1 Waveform inversion is performed to calculate the final inversion result. This enables stable and effective updating of the low wavenumber components in the gradient when the initial velocity model is poor, reconstructing a good background velocity model. The final inverted velocity model has higher accuracy than traditional waveform inversion results and avoids the inversion from quickly getting trapped in local extrema.

Description

TECHNICAL FIELD

[0001] The present application relates to a waveform inversion method based on structure shaping regularization, and belongs to the field of oil and gas geophysical exploration engineering. BACKGROUND

[0002] With the deepening of exploration and development in China, the accuracy requirement of seismic imaging is getting higher and higher, and an accurate velocity model is the key to high-precision seismic imaging. Due to the progress of computing power, full waveform inversion (FWI) has become a powerful and promising velocity modeling method, which updates iteratively to reduce the residual between observed seismic data and simulated seismic data, and finally realizes the reconstruction of high-precision underground velocity model.

[0003] Full waveform inversion (FWI) is a combination of tomographic mode and migration mode. However, when the initial model is far away from the true model, or lacks large offset and low frequency data, waveform inversion will be dominated by migration mode. This will suppress the update of low wave number components in the inversion process, and easily fall into local minimum, which seriously affects the inversion accuracy. SUMMARY

[0004] The purpose of the present application is to provide a waveform inversion method based on structure shaping regularization and its device, which can solve the problems existing in the prior art, alleviate the local extremum problem of waveform inversion, and effectively improve the velocity modeling accuracy.

[0005] To solve the above technical problems, the present application adopts the following technical scheme:

[0006] In a first aspect, a waveform inversion method based on structure shaping regularization is provided, comprising: obtaining observed data and an initial velocity model; performing seismic wave forward simulation calculation using the initial model to obtain simulated seismic data, and calculating the back-propagation residual wave field using the difference between the simulated seismic data and the observed data; simultaneously calculating the tomographic mode waveform inversion gradient G t and the migration mode waveform inversion gradient G m , wherein,

[0007]

[0008]

[0009] wherein, t max is the total time of the seismic record, H represents Hilbert transform, u represents forward-propagation source wave field, r represents back-propagation residual wave field, x s is the horizontal spatial coordinate of the shot point, z s is the vertical spatial coordinate of the shot point, F 2D and represent two-dimensional Fourier transform and two-dimensional inverse Fourier transform, respectively, T kis the threshold function, m is the velocity model parameter, u is the simulated seismic wave field, t is the calculation time of the forward modeling, x is the abscissa, and z is the ordinate.

[0010] Using the offset mode gradient G m Determine the local structural dip angle, and determine the background velocity model m based on the local structural dip angle. n+1 ;

[0011] According to the background velocity model m n+1 Waveform inversion is performed to calculate the final inversion result.

[0012] In a second aspect, an electronic device is provided, comprising:

[0013] At least one processor; and,

[0014] A memory communicatively connected to the at least one processor; wherein,

[0015] The memory stores instructions that can be executed by the at least one processor to enable the at least one processor to perform the method as described in any of the first aspects.

[0016] Compared with existing technologies, the method of the present invention acquires observation data and an initial velocity model; performs forward modeling calculations of seismic waves using the initial model to obtain simulated seismic data; calculates the backpropagation residual wavefield using the difference between the simulated seismic data and the observation data; and simultaneously calculates the tomographic mode waveform inversion gradient G. t and the gradient G of the offset mode waveform inversion m Using the offset mode gradient G m Determine the local structural dip angle, and determine the background velocity model m based on the local structural dip angle. n+1 According to the background velocity model m n+1 Waveform inversion is performed to calculate the final inversion result. This enables stable and effective updating of the low wavenumber components in the gradient, reconstructing a good background velocity model, even when the initial velocity model is poor. The final inverted velocity model has higher accuracy than traditional waveform inversion results, avoiding rapid entrapment in local extrema. Furthermore, this invention proposes a new inversion process that, without increasing computational load, fully utilizes the migration mode gradient for shaping regularization, improving the inversion effect of tomographic mode waveform inversion and providing a good initial velocity model for subsequent conventional waveform inversion. Attached Figure Description

[0017] Figure 1a A flowchart illustrating a waveform inversion method based on constructive shaping regularization provided in the embodiments of this specification;

[0018] Figure 1b This is a flowchart of another embodiment of the present invention;

[0019] Figure 2 This is a true speed model of one embodiment of the present invention;

[0020] Figure 3 This is an initial velocity model according to an embodiment of the present invention;

[0021] Figure 4 A single shot in the observed seismic data of one embodiment of the present invention;

[0022] Figure 5 This is a single shot in simulated seismic data calculated from an initial velocity model, according to an embodiment of the present invention.

[0023] Figure 6 This is a conventional waveform inversion result according to an embodiment of the present invention;

[0024] Figure 7 The tomographic mode waveform inversion result based on construction shaping regularization is shown in one embodiment of the present invention (scale 1).

[0025] Figure 8 The tomographic mode waveform inversion result based on construction shaping regularization is shown in one embodiment of the present invention (scale 2).

[0026] Figure 9 The tomographic pattern waveform inversion result based on construction shaping regularization according to an embodiment of the present invention (scale 3);

[0027] Figure 10 This is the final waveform inversion result of one embodiment of the present invention;

[0028] Figure 11 A simulated seismic record of the final waveform inversion result according to one embodiment of the present invention;

[0029] Figure 12 A single-channel comparison of the longitudinal extraction speed model of a real model, an initial model, the method of the present invention, and a conventional method, as an embodiment of the present invention;

[0030] Figure 13 The normalized objective functional convergence curves of the method of the present invention and the conventional waveform inversion method are shown in one embodiment of the present invention. Detailed Implementation

[0031] To enable those skilled in the art to better understand the technical solutions in this specification, the technical solutions in the embodiments of this specification will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this specification, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of this application.

[0032] In the first aspect, such as Figure 1a As shown, Figure 1a A flowchart illustrating a waveform inversion method based on constructive shaping regularization provided in this specification includes:

[0033] S1, acquire observation data and initial velocity model.

[0034] Input the observation data and the initial velocity model required for inversion. The observation data includes, but is not limited to, the observation system, observation records (i.e., shot data), and preset velocity modeling parameters, including the velocity model's lateral sampling points nx and longitudinal sampling points nz, lateral spatial sampling interval dx, longitudinal spatial sampling interval dz, time sampling interval dt, number of time sampling points nt, dominant frequency f0, number of shots ns in the seismic data, number of receivers per shot ng, shot-to-spot interval ds, and receiver-to-spot interval dg.

[0035] S2, use the initial model to perform forward modeling of seismic waves to obtain simulated seismic data, and use the difference between the simulated seismic data and the observed data to calculate the backpropagation residual wave field.

[0036] Seismic wave forward modeling was performed using the initial model to obtain simulated seismic data. The residual backpropagating wavefield was then calculated using the difference between the simulated and observed data. The seismic wave simulation process employed a two-dimensional constant-density acoustic wave equation, and the calculation method was as follows:

[0037]

[0038] Where m represents the velocity model parameters, u represents the simulated seismic wavefield, t represents the calculation time in the forward modeling, and (x, z) represent the spatial coordinates in the horizontal and vertical directions. The residuals are obtained by subtracting the simulated seismic data from the observed seismic data.

[0039] f r (x r ,z r ,t)=u(x r ,z r ,t)-u obs (x r ,z r ,t), (2)

[0040] The residual f r (x r ,z r ,t) is used as the source of the back propagation wave field. The residual back propagation wave field is calculated using formula (1) for subsequent gradient calculation.

[0041] S3, simultaneously calculate the tomographic mode waveform inversion gradient G t and the gradient G of the offset mode waveform inversion m .

[0042] The gradient of waveform inversion can be decomposed into tomographic and offset components. Therefore, the gradient decomposition can be expressed as:

[0043] G FWI =G t +G m (3)

[0044] Among them G FWI It is a traditional waveform inversion gradient, G t It is a chromatographic component used to update the background velocity, G m The offset component is used to update the velocity perturbation.

[0045] The gradient field of conventional waveform inversion is derived from the adjoint state method as the source propagation wave field u(x,z,t; x s ,z s ) and residual inverse propagated wavefield r(x,z,t; x s ,z s The cross-correlation of ) is represented as,

[0046]

[0047] Where t max For the total time of earthquake records, (x s ,z s ) represents the spatial coordinates of the gun point in the horizontal and vertical directions.

[0048] By separating the ascending and descending wave fields using the Hilbert transform and cross-correlated wave fields with the same propagation direction, the decomposed tomographic mode gradient can be derived. To mitigate the leakage of the offset component in the tomographic mode gradient and suppress this high wavenumber noise, this method applies a multi-scale wavenumber domain two-dimensional filter to obtain an optimized tomographic mode gradient, which is then used to update the background velocity field.

[0049]

[0050] Where H represents the Hilbert transform, u represents the forward propagation source wave field, r represents the backward propagation residual wave field, (x s ,z sF represents the spatial coordinates of the shot point in the horizontal and vertical directions. 2D and T represents the two-dimensional forward Fourier transform and the two-dimensional inverse Fourier transform, respectively. k This is a threshold function.

[0051] By cross-correlating wavefields with opposite propagation directions, the decomposed migration mode gradient G can be derived. m :

[0052]

[0053] Since the two gradient fields after decomposition differ only in the sign of the second term, both tomographic and migration gradients can be obtained simultaneously with almost no increase in computation.

[0054] S4, using the offset mode gradient G m Determine the local structural dip angle, and determine the background velocity model m based on the local structural dip angle. n+1 .

[0055] For example, the offset mode gradient G can be used. m And based on the plane wave deconstruction filter (PWD) method, it performs local construction tilt angle estimation.

[0056]

[0057] σ is the local structural dip angle, and W is the local window for calculating the dip angle.

[0058] Then, the velocity model is updated based on the conjugate gradient algorithm, and a constructive regularization is performed. After iteration, the background velocity model is obtained. The new iterative formula is written as follows:

[0059] m n+1 =S -1 T s S(m n +α n G t,n (10)

[0060] Where n is the number of iterations, m n For the velocity model of the nth iteration, α n This is the step size for the nth iteration. The gradient used in the iteration is the tomographic gradient G. t,n This emphasizes the update amount of low-wavelength array components in the velocity model, avoiding a velocity model update dominated by high-wavelength array components, which could cause the inversion to quickly fall into local extrema. S and S -1 These are the Seislet forward and inverse transforms, respectively, T sIt is a threshold function. The construction of a local tilt angle σ is introduced into the Seislet transform; this introduction is reflected in the prediction operator P and the update operator U in the Seislet transform. The prediction and update operators defined by the PWD-based Seislet transform are as follows:

[0061]

[0062]

[0063] Where e is an even sequence, r = oP(e) is the residual between the predicted values ​​of the odd sequence o and the even sequence e, and i is the sequence number. and This represents a prediction operator that drifts adjacent seismic traces based on the local structural dip angle σ. + and - correspond to drifting on the left and right sides along the local structural dip angle, respectively.

[0064] S5, based on the background velocity model m n+1 Waveform inversion is performed to calculate the final inversion result.

[0065] The background velocity model m output by iteration n+1 Input a regular waveform for inversion and calculate the final inversion result.

[0066] The method of this invention decomposes the gradient into two parts: (1) an optimized tomographic mode gradient, which enhances the update of long-wavelength components in the early stage of inversion, avoids the leakage of offset components, and thus alleviates the problem of inversion getting trapped in local extrema; (2) an offset mode gradient, which can be used to obtain structural dip information, add structural shaping regularization, enhance the stability of inversion, and improve the reconstruction accuracy of the background velocity model. Based on the inversion process proposed in this invention, a more accurate background velocity model can be obtained as the initial model for conventional waveform inversion. Even when the initial model differs greatly from the real model, the method of this invention can still obtain high-precision waveform inversion results.

[0067] By acquiring observational data and an initial velocity model, forward modeling of seismic waves is performed using the initial model to obtain simulated seismic data. The difference between the simulated seismic data and the observational data is used to calculate the backpropagation residual wavefield. Simultaneously, the tomographic mode waveform inversion gradient G is calculated. t and the gradient G of the offset mode waveform inversion m Using the offset mode gradient G m Determine the local structural dip angle, and determine the background velocity model m based on the local structural dip angle. n+1 According to the background velocity model m n+1Waveform inversion is performed to calculate the final inversion result. This enables stable and effective updating of the low wavenumber components in the gradient, reconstructing a good background velocity model, even when the initial velocity model is poor. The final inverted velocity model has higher accuracy than traditional waveform inversion results, avoiding rapid entrapment in local extrema. Furthermore, this invention proposes a new inversion process that, without increasing computational load, fully utilizes the migration mode gradient for shaping regularization, improving the inversion effect of tomographic mode waveform inversion and providing a good initial velocity model for subsequent conventional waveform inversion.

[0068] To illustrate the method of this invention more specifically, the Marmousi model is used as an example. The specific process is as follows: Figure 1b As shown. The real velocity model used in this embodiment is as follows. Figure 2 As shown, the initial velocity model is as follows: Figure 3 As shown. The velocity model has a grid size of 384×136 and a grid spacing of 12 meters.

[0069] The observation system is distributed as follows: 128 shots are evenly distributed across the surface at 36-meter intervals. Each shot receives all 384 geophones, with a geophone interval of 12 meters. The time sampling interval is 1 millisecond, with 2000 time sampling points and a dominant frequency of 10 Hz. The observation records are calculated based on a true velocity model, and the seismic records of a single shot are as follows: Figure 4 As shown. In contrast, simulated earthquake records based on the initial model are as follows: Figure 5 As shown. Then, the residuals between the observed and simulated records are calculated, and the residual backpropagated wavefield is obtained by reverse time propagation.

[0070] Figure 6 The results are from a conventional waveform inversion performed directly on the initial model. After 200 iterations, a significant high-speed anomaly is observed on the left side of the model. This indicates that the conventional inversion method is trapped in local extrema, resulting in inversion results that deviate significantly from the true model and exhibit low inversion accuracy.

[0071] This invention calculates the gradients of the tomographic and migration modes optimized in the wavenumber domain, and employs a constructed shaping regularized iterative formula to extract low wavenumber components and update the background velocity model in the early stages of inversion. This embodiment utilizes a multi-scale inversion strategy, and the multi-scale inversion results are as follows: Figure 7 , Figure 8 , Figure 9 As shown, each of the three inversion scales was iterated 20 times.

[0072] After multi-scale inversion, the tomographic waveform inversion results based on construction shaping regularization are obtained, such as... Figure 9 As shown. Using this model as the initial model for conventional waveform inversion, the final inversion result is as follows. Figure 10As shown in the figure, the velocity modeling accuracy of this invention is significantly better than that of conventional waveform inversion methods, based on the velocity model inversion results. For the low-velocity region on the left side of the model, shallow steep faults, and deep anticline structures, this invention can achieve relatively accurate reconstruction.

[0073] Figure 11 For simulated seismic data based on the inversion results of the method of this invention, and Figure 4 The comparison shows that the simulated seismic data in the final inversion results are similar to the observed seismic data, and the waveforms in the single-shot records have consistent amplitude and phase.

[0074] Figure 12 The velocity curve extracted from the Distance = 1200 meters in the inversion results of the real velocity model, the initial velocity model, the conventional waveform inversion method, and the inversion results of the present invention is similar to that of the real model and has higher inversion accuracy than that of the conventional waveform inversion results.

[0075] Figure 13 The figures show the convergence curves of the normalized objective functional for the conventional waveform inversion method and the method of this invention. While the conventional waveform inversion method converges stably and rapidly in the first 30 iterations, it eventually converges to a large value due to getting trapped in local extrema. The objective functional of this invention can stably converge to a small value close to zero, indicating that the final inversion result is very close to the true model.

[0076] In a second aspect, embodiments of this specification also provide an electronic device, the device comprising:

[0077] At least one processor; and,

[0078] A memory communicatively connected to the at least one processor; wherein,

[0079] The memory stores instructions that can be executed by the at least one processor to enable the at least one processor to perform the method as described in the first aspect.

[0080] The foregoing has described specific embodiments of this specification. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recited in the claims may be performed in a different order than that shown in the embodiments and may still achieve the desired result. Furthermore, the processes depicted in the drawings do not necessarily require the specific or sequential order shown to achieve the desired result. In some embodiments, multitasking and parallel processing are possible or may be advantageous.

[0081] The above description is merely one or more embodiments of this specification and is not intended to limit this specification. Various modifications and variations can be made to the one or more embodiments of this specification by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principle of one or more embodiments of this specification should be included within the scope of the claims of this specification.

Claims

1. A waveform inversion method based on structural shape regularization, comprising: obtaining observed data and an initial velocity model; performing seismic wave forward simulation calculation using the initial model to obtain simulated seismic data, and calculating inverse transmission residual wave field using the difference between the simulated seismic data and the observed data; Simultaneously compute tomographic model waveform inversion gradient G t and migration model waveform inversion gradient G m wherein, where t max is the total time of the seismic record, H represents the Hilbert transform, u represents the forward propagating source wavefield, r represents the backward propagating residual wavefield, x s is the horizontal spatial coordinate of the shot point, z s is the vertical spatial coordinate of the shot point, F 2D and represent the two-dimensional Fourier forward transform and the two-dimensional Fourier inverse transform, respectively, T k is the threshold function, m is the velocity model parameter, u is the simulated seismic wavefield, t is the calculation time of the forward simulation, x is the horizontal coordinate, and z is the vertical coordinate. using the offset pattern gradient G m determining a local structural dip angle, determining a background velocity model m from the local structural dip angle n+1 ; According to the background velocity model m n+1 Waveform inversion is performed to obtain a final inversion result.

2. The method of claim 1, wherein, using the offset pattern gradient G m determining the local structural dip angle comprises: performing local structural dip estimation in the following manner, where σ is the local configuration angle and W is the local window over which the angle is computed.

3. The method of claim 2, wherein, determining a background velocity model m from the local configuration dip angle n+1 comprising: The background velocity model is iterated n times in the following way, m n+1 = S -1 T s S(m n + a n G t,n ), where n is the iteration number, m n is the velocity model of the nth iteration, a n is the step size of the nth iteration, S and S -1 are the forward and inverse Seislet transforms, respectively, T s is a threshold function, where the prediction operator P and the update operator U in the Seislet transform are as follows: where e is an even sequence, r = o - P(e) is the residual of the odd sequence o and the prediction value of the even sequence e, i is the sequence number, and denotes the prediction operator that shifts the adjacent seismic traces according to the local structural dip angle σ, + and - correspond to the shift along the left and right sides of the local structural dip angle, respectively. 4.An electronic device, comprising: at least one processor; and a memory connected to the at least one processor in communication; wherein the memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to perform the method of any one of claims 1 to 3.

Images