Sound wave full-waveform inversion method and system based on optimized effective boundary storage

Abstract

The invention discloses a sound wave full waveform inversion method and system based on optimization effective boundary storage, and the method comprises the steps: constructing an initial geological velocity model, and obtaining seismic wave parameter information; performing forward modeling according to the initial geological velocity model and the seismic wave parameter information to obtain seismic simulation data and forward seismic source wave field data; acquiring actual seismological observation data of the geologic model to be measured, and constructing an objective function according to the actual seismological observation data and the seismic simulation data; re-sampling the wave field data of the target area according to the forward seismic source wave field data; restoring the resampled wave field data of the target area to an original data state by using an interpolation method, and then reconstructing a forward transmission seismic source wave field; zero-delay cross-correlation is carried out according to the reconstructed forward transmission seismic source wave field and the reconstructed backward transmission residual wave field to obtain a gradient; and carrying out iterative updating on the target function by adopting a gradient optimization algorithm until a final geological velocity model is output.

Description

Technical Field

[0001] The invention relates to the technical field of seismic physical exploration, and in particular to an acoustic full waveform inversion method and system based on optimized effective boundary storage. Background Art

[0002] Since the full waveform inversion theory was proposed by Laily (1983) and Tarantola (1984) in the early 1980s, it has been widely studied in many fields such as oil and gas and mineral resource exploration and development. The full waveform inversion method continuously updates the underground velocity model through a series of iterative processes, aiming to minimize the difference between seismic simulation data and actual seismic observation data. Through repeated iterations, the simulated seismic wave field gradually approaches the actual seismic observation wave field, thereby achieving high-precision reconstruction of the real physical parameters of the underground (such as wave velocity, density, etc.). Compared with traditional seismic imaging methods, full waveform inversion can more finely characterize underground structural characteristics. However, the method has huge computational and storage requirements, and still faces certain limitations and challenges in practical applications. Summary of the invention

[0003] The present invention provides an acoustic full waveform inversion method based on optimized effective boundary storage, which can effectively reduce the calculation and storage requirements in the full waveform inversion process, thereby improving the inversion efficiency.

[0004] The purpose of the embodiment of this specification is to provide an acoustic full waveform inversion method based on optimized effective boundary storage, the method comprising:

[0005] Construct an initial geological velocity model and obtain seismic wave parameter information;

[0006] Perform forward modeling based on the initial geological velocity model and the seismic wave parameter information to obtain seismic simulation data and forward source wave field data;

[0007] Acquire actual seismic observation data of the geological model to be tested, and construct an objective function based on the actual seismic observation data and the seismic simulation data;

[0008] Resampling the target area wavefield data according to the forward source wavefield data; restoring the resampled target area wavefield data to the original data state by using an interpolation method, and then reconstructing the forward source wavefield;

[0009] The gradient is obtained by performing zero-delay cross-correlation based on the reconstructed forward source wave field and the reverse residual wave field; the objective function is iteratively updated using a gradient optimization algorithm until a final geological velocity model is output; the final geological velocity model is the inversion result.

[0010] Preferably, the method for performing the forward simulation includes: performing a first-order stress-acoustic wave equation according to the initial geological velocity model and the seismic wave parameter information, and the expression includes:

[0011]

[0012] Among them, P represents stress; Vx and Vz represent the particle velocity components in the x and z directions respectively; v represents the longitudinal wave velocity; s represents the earthquake source; ρ represents the density; t represents the time; x and z represent the grid x and z coordinate values.

[0013] Preferably, the objective function constructed according to the actual earthquake observation data and the earthquake simulation data is a residual two-norm objective function:

[0014]

[0015] Where D represents the observation target area; u syn and u obs Represent earthquake simulation data and actual earthquake observation data respectively; δ(xx r ) represents the unit pulse function; T max represents the maximum time earthquake record; m represents the velocity model parameter; x r Indicates the location of the detector.

[0016] Preferably, when resampling the target area wavefield data according to the forward source wavefield, the sampling interval formula is obtained according to the Nyquist sampling theorem, and the expression includes:

[0017]

[0018] Where T is the sampling interval obtained according to the Nyquist sampling theorem, f max The maximum frequency of the current sampling.

[0019] The interpolation method is used to restore the resampled wave field data to the original data state, and the steps include:

[0020]

[0021] Where t is time, P(x,t) represents the wavefield information at a specific spatial position x and time t through the interpolation method. P(x,t0) is the resampled wavefield data at time t0, which is the left endpoint of the interpolation process; P(x,t1) is the resampled wavefield data at time t1, which is the right endpoint of the interpolation process; x(x,z) is the spatial position in the x and z directions.

[0022] Preferably, the back propagation residual wave field is calculated according to the difference between the actual seismic observation data and the seismic simulation data, the steps comprising:

[0023] According to the Lagrange principle, a first-order acoustic wave equation constraint is imposed on the residual two-norm objective function:

[0024]

[0025]

[0026] Add the constraint term to the residual two-norm objective function and change the objective function to:

[0027]

[0028] Wherein, λ(x,t) represents the accompanying wave field vector, also known as the Lagrange multiplier; E(m)=0 is the above-mentioned first-order acoustic wave equation;

[0029] According to the above formula, the adjoint equation is derived, and the reverse time propagation equation of the residual wave field with the source as the source is obtained. The formula includes:

[0030]

[0031] Among them, S λp Represents the companion source, which is the difference between the earthquake simulation record and the actual earthquake observation record. Vx represents the Lagrange multiplier in the x direction, λ Vz represents the Lagrange multiplier in the z direction, λ P Denotes the companion wave field.

[0032] Preferably, performing cross-correlation on the back-propagated residual wavefield and the reconstructed source wavefield to obtain the objective function gradient includes:

[0033] The gradient calculation formula of velocity v is:

[0034]

[0035] in, represents the gradient of H2 function with respect to v; represents the gradient of the H2 function with respect to the bulk modulus.

[0036] Preferably, the initial geological velocity model is updated using a gradient optimization algorithm to obtain a final velocity model, including:

[0037] Determining the updating direction of the geological velocity model according to the obtained gradient;

[0038] According to the updating direction of the geological velocity model, the initial geological velocity model is updated to obtain a final velocity model.

[0039] The present invention also provides an acoustic full waveform inversion system based on optimized effective boundary storage, the system is used to implement the above method, including: a construction module, a first calculation module, a second calculation module, a third calculation module and an update module;

[0040] The construction module is used to construct an initial geological velocity model and obtain seismic wave parameter information;

[0041] The first calculation module is used to perform forward simulation according to the initial geological velocity model and the seismic wave parameter information to obtain seismic simulation data and forward source wave field data;

[0042] The second calculation module is used to obtain actual seismic observation data of the geological model to be measured, and to construct an objective function according to the actual seismic observation data and the seismic simulation data;

[0043] The third calculation module is used to resample the target area wavefield data according to the forward source wavefield data; restore the resampled target area wavefield data to the original data state by using an interpolation method, and then reconstruct the forward source wavefield;

[0044] The updating module is used to obtain the gradient by performing zero-delay cross-correlation based on the reconstructed forward source wave field and the reverse residual wave field; the objective function is iteratively updated using a gradient-based optimization algorithm until a final geological velocity model is output; the final geological velocity model is the inversion result.

[0045] Compared with the prior art, the present invention has the following beneficial effects:

[0046] Based on the time-domain acoustic full waveform inversion, a resampling algorithm based on the Nyquist sampling theorem is applied to the effective boundary storage, aiming to significantly reduce the storage capacity requirements of wavefield data. The resampled wavefield data is restored using linear interpolation to complete the source wavefield reconstruction and integrate it into the time-domain full waveform inversion process. This solution can significantly reduce the pressure on data storage and thus accelerate the calculation speed of full waveform inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] In order to more clearly illustrate the technical solution of the present invention, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0048] Figure 1 It is a flow chart of an embodiment of an acoustic full waveform inversion method based on optimized effective boundary storage provided in this specification;

[0049] Figure 2 It is a flowchart of an acoustic full waveform inversion method based on optimized effective boundary storage in a specific scenario example provided in this specification;

[0050] Figure 3 a is a schematic diagram of comparison between original data and interpolated data of wave field data in an embodiment of the present invention when the sampling interval is 2.5 ms;

[0051] Figure 3 b is a schematic diagram of comparison between original data and interpolated data of wave field data in an embodiment of the present invention when the sampling interval is 5 ms;

[0052] Figure 3 c is a schematic diagram of comparison between original data and interpolated data of wave field data in an embodiment of the present invention when the sampling interval is 25 ms;

[0053] Figure 3 d is a schematic diagram showing a comparison of relative errors between original data and interpolated data of wave field data in an embodiment of the present invention when the sampling interval is 2.5 ms;

[0054] Figure 3 e is a schematic diagram showing a comparison of relative errors between original data and interpolated data of wave field data in an embodiment of the present invention when the sampling interval is 5 ms;

[0055] Figure 3 f is a schematic diagram showing a comparison of relative errors between original data and interpolated data of wave field data in an embodiment of the present invention when the sampling interval is 25 ms;

[0056] Figure 4 a is a conventional wave field reconstructed image at 250ms in an embodiment of the present invention;

[0057] Figure 4 b is a wave field reconstructed image after resampling at 250ms in the embodiment of the present invention;

[0058] Figure 4 c is a schematic diagram showing a comparison of relative errors between a conventional wave field reconstruction image at 250 ms and a wave field reconstruction after resampling in an embodiment of the present invention;

[0059] Figure 5 a is a conventional wave field reconstructed image at 500ms in the embodiment of the present invention;

[0060] Figure 5 b is a wave field reconstructed image after resampling at 500ms in the embodiment of the present invention;

[0061] Figure 5c is a schematic diagram showing a comparison of relative errors between a conventional wave field reconstruction image at 500ms and a wave field reconstruction after resampling in an embodiment of the present invention;

[0062] Figure 6 a is the initial velocity model diagram based on the two-dimensional Overthrust model;

[0063] Figure 6 b is the real velocity model diagram based on the two-dimensional Overthrust model;

[0064] Figure 6 c is a full waveform inversion result diagram based on the linear interpolation method provided in this application;

[0065] Figure 7 A schematic diagram of the structure of an embodiment of an acoustic full waveform inversion system based on optimized effective boundary storage provided in this specification; DETAILED DESCRIPTION

[0066] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0067] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0068] Embodiment 1

[0069] This specification proposes an acoustic full waveform inversion method based on optimized effective boundary storage. First, an initial geological velocity model is constructed; actual seismic observation data of the geological model to be tested is input; seismic wave parameter information; forward simulation is performed based on the initial geological velocity model and seismic wave parameter information to obtain seismic simulation data and forward source wave field data; an objective function is constructed based on actual seismic observation data and seismic simulation data; the wave field data of the target area is resampled based on the forward source wave field; the resampled wave field data is restored to the original data state using an interpolation method, and then the source wave field is reconstructed; zero-delay cross-correlation is performed based on the reconstructed source wave field and the back-transmitted residual wave field to obtain the gradient. A gradient-based optimization algorithm is used to iteratively update the initial geological velocity model until the final geological velocity model is output;

[0070] like Figure 1 As shown, the embodiment of the present invention provides an acoustic full waveform inversion method based on optimized effective boundary storage. In specific implementation, the method may include the following contents.

[0071] Step 101. Construct an initial geological velocity model and obtain seismic wave parameter information;

[0072] Step 102: Perform forward simulation based on the initial geological velocity model and seismic wave parameter information to obtain seismic simulation data and forward source wave field data;

[0073] Step 103. Acquire actual earthquake observation data of the geological model to be tested, and construct an objective function based on the actual earthquake observation data and earthquake simulation data;

[0074] Step 104. Resample the target area wavefield data according to the forward source wavefield data; restore the resampled target area wavefield data to the original data state using an interpolation method, and then reconstruct the forward source wavefield;

[0075] Step 105. Obtain the gradient by performing zero-delay cross-correlation based on the reconstructed forward source wavefield and the reverse residual wavefield; iteratively update the objective function using a gradient-based optimization algorithm until the final geological velocity model is output; the final geological velocity model is the inversion result.

[0076] Depend on Figure 2 It can be seen from the process shown that the acoustic full waveform inversion method for optimizing effective boundary storage in the embodiment of the present invention directly uses the first-order stress-velocity acoustic wave equation for forward simulation, and applies the resampling algorithm of the Nyquist sampling theorem for effective boundary storage, aiming to significantly reduce the storage capacity requirements of wave field data. The wave field data after resampling is restored using the linear interpolation method to complete the reconstruction of the source wave field, and is integrated into the process of full waveform inversion in the time domain, which can significantly reduce the pressure of data storage and thus accelerate the calculation speed of full waveform inversion.

[0077] 101. Construct an initial geological velocity model and obtain seismic wave parameter information;

[0078] In some embodiments, the initial geological velocity model is based on prior knowledge of seismic exploration, assumptions of geological structures, and possible geophysical data to construct a preliminary geological velocity model for subsequent iterative optimization. This model should be able to reflect the approximate velocity distribution characteristics of the underground medium.

[0079] In some embodiments, seismic wave parameters contain rich geological information, among which the arrival time of seismic waves is the key basis for determining the propagation path and propagation speed of seismic waves; and the waveform data contains important parameters such as the amplitude and frequency of seismic waves.

[0080] 102. Perform forward simulation based on the initial geological velocity model and seismic wave parameter information to obtain seismic simulation data and forward source wave field data;

[0081] In some embodiments, forward modeling is performed based on the initial geological velocity model and seismic wave parameter information to obtain seismic simulation data and forward source wave field data, specifically including:

[0082] According to the initial geological velocity model and seismic wave parameter information, the first-order stress-velocity acoustic wave equation in the two-dimensional isotropic homogeneous medium in the time domain is used for forward simulation as follows:

[0083]

[0084] Among them, P represents stress; Vx and Vz represent the particle velocity components in the x and z directions respectively; v represents the longitudinal wave velocity; s represents the earthquake source; ρ represents the density; and t represents the time.

[0085] 103. Obtain actual earthquake observation data of the geological model to be tested, and construct an objective function based on the actual earthquake observation data and earthquake simulation data;

[0086] In some embodiments, the actual seismic observation data of the geological model to be measured is used by rationally arranging multiple detectors in the geological area to be measured. When the shot point vibrates, the seismic wave will propagate in the underground medium, and the wave field received by the detection point is used as the actual seismic observation wave field.

[0087] In some embodiments, the residual second norm is a quantitative indicator to measure the degree of such matching. It is obtained by calculating the differences in time and space between the simulated wave field and the observed wave field, and taking the square root of the sum of the squares of these differences. If the residual second norm is small, it means that the seismic simulation wave field is very close to the actual observed wave field, and the model has a good fitting effect on the actual situation; on the contrary, a large residual second norm indicates that there is a large deviation between the model and the actual observation, and it is necessary to use an optimization algorithm to update the velocity model, so as to obtain a seismic simulation wave field that is closer to the actual observed wave field.

[0088] In some embodiments, the residual two-norm objective function of the actual earthquake observation data and the earthquake simulation data is calculated by the following formula:

[0089]

[0090] Where D represents the observation target area; u syn and u obs Represent earthquake simulation data and actual earthquake observation data respectively; δ(xx r ) represents the unit pulse function; T max represents the maximum time earthquake record; m represents the velocity model parameter; x r Indicates the location of the detector.

[0091] 104. Resample the wavefield data of the target area according to the forward source wavefield data; restore the resampled wavefield data of the target area to the original data state by using the interpolation method, and then reconstruct the forward source wavefield;

[0092] In some embodiments, a reasonable sampling interval formula is obtained according to the Nyquist sampling theorem as follows:

[0093]

[0094] Where T is the sampling interval obtained according to the Nyquist sampling theorem, f max The maximum frequency of the current sampling.

[0095] In some embodiments, the resampled boundary wavefield values ​​are interpolated and restored according to a linear interpolation method.

[0096]

[0097] Where t is time, and P(x,t) represents the wavefield information at a specific spatial position x and time t through the interpolation method. P(x,t0) is the resampled wavefield data at time t0, which is used as the left endpoint of the interpolation process, and P(x,t1) is the resampled wavefield data at time t1, which is used as the right endpoint of the interpolation process to solve the wavefield value at any position between the two. x(x,z) is the spatial position in the x and z directions.

[0098] 105. The gradient is obtained by performing zero-delay cross-correlation based on the reconstructed forward source wave field and the reverse residual wave field; the objective function is iteratively updated using a gradient optimization algorithm until the final geological velocity model is output; the final geological velocity model is the inversion result;

[0099] In some embodiments, the gradient is obtained by performing zero-delay cross-correlation based on the reconstructed source wave field and the back-propagated residual wave field. The initial geological velocity model is iteratively updated using a gradient optimization algorithm until the final geological velocity model is output, specifically including:

[0100] According to the difference between the actual seismic observation data and the seismic simulation data, the back propagation residual wave field is calculated;

[0101] The objective function gradient is obtained by cross-correlating the back-propagated residual wave field and the reconstructed source wave field.

[0102] According to the objective function gradient, the initial geological velocity model is updated using the gradient optimization algorithm to obtain the final velocity model.

[0103] In some embodiments, the back-propagated residual wavefield and the reconstructed source wavefield are cross-correlated to obtain the objective function gradient.

[0104] In some embodiments, the calculation process of the back-propagated residual wave field and the calculation process of the gradient are described:

[0105] According to the Lagrange principle, a first-order acoustic wave equation constraint is imposed on the objective function:

[0106]

[0107]

[0108] By adding E(m)=0 (Formula 6) as a constraint term to the constraint of the residual two-norm objective function (i.e., Formula 2), the objective function can be changed to:

[0109]

[0110] Where λ(x,t)=(λ Vx ,λ Vz ,λ P ) T represents the adjoint wave field vector, also known as the Lagrange multiplier, λ Vx represents the Lagrange multiplier in the x direction, λ Vz represents the Lagrange multiplier in the z direction, λ P Denotes the companion wave field.

[0111] Expanding formula 7 gives the following form:

[0112]

[0113] To obtain the minimum value of formula 7, the following three optimization conditions need to be met:

[0114]

[0115] make Combined with the boundary conditions and termination conditions, the calculation formula of the accompanying wave field, i.e. the back-propagation residual wave field, can be derived:

[0116]

[0117] Among them, S λp represents the accompanying source, which means the difference between the earthquake simulation record and the actual earthquake observation record; Vx represents the Lagrange multiplier in the x direction, λ Vz represents the Lagrange multiplier in the z direction, λ P Denotes the companion wave field.

[0118] For formula (7)ρv 2 Taking partial derivatives, we can get:

[0119]

[0120] The gradient calculation formula for velocity v is:

[0121]

[0122] in, represents the gradient of the H2 function with respect to v, represents the gradient of the H2 function with respect to the bulk modulus.

[0123] According to formula 15, we can get: The H2 function is the gradient of v, which is obtained by integrating the zero-delay cross-correlation of the back-propagated residual wavefield and the reconstructed source wavefield along time.

[0124] In some embodiments, the gradient is calculated using Formula 15. The full-time cross-correlation calculation of the current technology is from 0 to T max The wave field data at all times in the entire time period requires a huge amount of calculation. From formulas 3 and 4, it can be seen that the resampling algorithm of the effective boundary wave field based on the Nyquist sampling theorem greatly reduces the storage capacity requirements of the wave field data. The resampled wave field data is restored using the linear interpolation method, which can effectively reduce the amount of calculation and storage, thereby improving the calculation speed of the full waveform inversion.

[0125] In some embodiments, the initial geological velocity model is updated using a gradient optimization algorithm to obtain a final velocity model, specifically including:

[0126] Determine the updating direction of the geological velocity model based on the obtained gradient;

[0127] The initial geological velocity model is updated according to the updating direction of the geological velocity model.

[0128] In some embodiments, determining the updating direction of the geological velocity model according to the obtained gradient specifically includes:

[0129] The update direction of the velocity model is calculated as follows:

[0130]

[0131] Where U represents the update direction of the velocity model.

[0132] In some embodiments, the initial geological velocity model is updated according to the update direction of the geological velocity model, specifically including:

[0133] The initial geological velocity model is updated as follows:

[0134] m i+1 =m i +α i U (17)

[0135] Among them, α i represents the update step size, m i Represents velocity model parameters.

[0136] In some embodiments, to verify the effectiveness of the interpolation method, a uniform velocity model is designed, the length and depth of the model are both 0.804km, the horizontal and vertical grid spacing is 4m, and the absorption boundary is 22 layers. First, the wave field propagation in the model is simulated, the duration is 0.5s, the time interval is 0.5ms, and the sub-wave used is the Ricker sub-wave with a main frequency of 10Hz. According to the Nyquist sampling theorem, the maximum time sampling interval can be obtained to be 10.5ms, Figure 3 a is a schematic diagram of comparison between original data and interpolated data of wave field data in an embodiment of the present invention when the sampling interval is 2.5 ms; Figure 3 b is a schematic diagram of comparison between original data and interpolated data of wave field data in an embodiment of the present invention when the sampling interval is 5 ms; Figure 3 c is a schematic diagram of comparison between original data and interpolated data of wave field data in an embodiment of the present invention when the sampling interval is 25 ms; Figure 3 d is a schematic diagram showing a comparison of relative errors between original data and interpolated data of wave field data in an embodiment of the present invention when the sampling interval is 2.5 ms; Figure 3 e is a schematic diagram showing a comparison of relative errors between original data and interpolated data of wave field data in an embodiment of the present invention when the sampling interval is 5 ms; Figure 3 f is a schematic diagram of the relative error comparison between the original data and the interpolated data of the wave field data of the embodiment of the present invention when the sampling interval is 25ms; Figure 3 a to Figure 3 c It can be seen that when the seismic data is resampled using a time step smaller than the sampling interval, the original data and the interpolated data are closely matched. On the contrary, if the sampling interval is larger than the maximum sampling interval allowed by the Nyquist sampling theorem, there will be a significant deviation between the interpolated data and the original data. Figure 3 e to Figure 3 It can be seen from f that the time step less than or equal to the Nyquist sampling interval should be selected as the sampling rate. The sampling interval accuracy obtained within the Nyquist sampling theorem is close.

[0137] In some embodiments, Figure 4 a is a conventional wave field reconstructed image at 250ms in an embodiment of the present invention; Figure 4 b is a wave field reconstructed image after resampling at 250ms in the embodiment of the present invention; Figure 4 c is a schematic diagram showing the relative error comparison between the conventional wave field reconstruction image at 250ms and the wave field reconstruction after resampling in the embodiment of the present invention; Figure 4 a to Figure 4c It can be seen that by comparing the conventional wave field reconstruction at 250ms with the wave field reconstruction image after resampling, it is found that the reconstructed images before and after sampling are highly consistent in morphology. In order to further see the difference between the two, the relative error calculation was performed on the two, and it was found that the difference between the two is extremely small.

[0138] In some embodiments, Figure 5 a is a conventional wave field reconstructed image at 500ms in the embodiment of the present invention; Figure 5 b is a wave field reconstructed image after resampling at 500ms in the embodiment of the present invention; Figure 5 c is a schematic diagram showing a comparison of relative errors between a conventional wave field reconstruction image at 500ms and a wave field reconstruction after resampling in an embodiment of the present invention; Figure 5 a to Figure 5 c It can be seen that by comparing the conventional wave field reconstruction at 500ms and the wave field reconstruction image after resampling, it is found that the reconstructed images before and after the sampling process are highly consistent in morphology. In order to further see the difference between the two, the relative error calculation was performed on the two, and it was found that the difference between the two is extremely small.

[0139] In a specific scenario example, Figure 6 a is the initial velocity model diagram based on the two-dimensional Overthrust model; Figure 6 b is a true velocity model diagram based on the two-dimensional Overthrust model. The true velocity model is used to obtain earthquake observation records; Figure 6 c is a full waveform inversion result diagram based on the linear interpolation method provided in this application; Figure 6 c means use Figure 6a The initial velocity model adopts the result diagram of the full waveform inversion based on linear interpolation provided by this application. The length of the model is 6.528km, the depth is 1.544km, and the total simulation time is 2.8s. The horizontal and vertical grid spacings of the velocity model are both 8m, the absorption boundary is 22 layers, the time interval is 0.8ms, the main frequency settings are 2.5Hz, 5.5Hz, 12Hz, 16Hz and 20Hz, and the maximum time sampling intervals are 43.2ms, 19.2ms, 8.8ms, 6.4ms and 4.8ms, respectively. The effective boundary wave field data after resampling is restored by linear interpolation and the wave field is reconstructed, and the reconstructed wave field is integrated into the time domain acoustic full waveform inversion process. The inversion process is gradually approached from the initial model through iterative optimization, and the final inversion result is consistent with the true velocity model. The method of this application determines the memory saving based on the sampling interval size obtained by the Nyquist sampling theorem. The larger the sampling interval, the more memory is saved. This result confirms that the interpolation sparsification method can not only effectively implement full waveform inversion, but also has high inversion accuracy, significantly reduces data storage requirements, accelerates the calculation process, and thus improves inversion efficiency.

[0140] Based on the above embodiments, the present application provides an acoustic full waveform inversion method based on optimized effective boundary storage, which adopts the first-order stress-velocity acoustic wave equation, obtains the value range of the sampling time step according to the Nyquist theorem, resamples the effective boundary wave field in the time domain, and then uses the interpolation method to restore the resampled wave field data to the original data state, and then reconstructs the source wave field, and performs zero-delay cross-correlation between the reconstructed source wave field and the back-transmitted residual wave field to obtain the gradient. It can effectively reduce the storage amount of data and speed up the transmission efficiency of data.

[0141] Embodiment 2

[0142] Based on the above-mentioned acoustic full waveform inversion method based on optimized effective boundary storage, the present invention also provides an embodiment of an acoustic full waveform inversion system based on optimized effective boundary storage, such as Figure 7 As shown in the figure, the acoustic full waveform inversion system based on optimized effective boundary storage specifically includes the following modules:

[0143] Construction module 701: used to construct an initial geological velocity model and obtain seismic wave parameter information;

[0144] The first calculation module 702 is used to perform forward simulation according to the initial geological velocity model and seismic wave parameter information to obtain seismic simulation data and forward source wave field data;

[0145] The second calculation module 703 is used to obtain the actual seismic observation data of the geological model to be measured, and to construct the objective function according to the actual seismic observation data and the seismic simulation data;

[0146] The third calculation module 704 is used to resample the target area wavefield data according to the forward source wavefield data; restore the resampled target area wavefield data to the original data state by using the interpolation method, and then reconstruct the forward source wavefield;

[0147] Update module 705: used to obtain the gradient by performing zero-delay cross-correlation based on the reconstructed forward source wave field and the reverse residual wave field; iteratively update the objective function using a gradient optimization algorithm until the final geological velocity model is output; the final geological velocity model is the inversion result.

[0148] In some embodiments, the first calculation module 702 can be specifically used to perform forward simulation according to the initial geological velocity model and seismic wave parameter information to obtain seismic simulation data and forward source wave field data.

[0149] In some embodiments, the second calculation module 703 may be specifically used to construct a residual two-norm objective function according to actual earthquake observation data and earthquake simulation data;

[0150] In some embodiments, the third calculation module 704 can be specifically used to resample the wave field data of the target area according to the forward source wave field; use the interpolation method to restore the resampled wave field data to the original data state, and then reconstruct the source wave field, including: according to the Nyquist sampling theorem, a reasonable sampling interval is calculated, the effective boundary wave field is resampled in the time domain, and the resampled boundary wave field value is interpolated and restored according to the linear interpolation method.

[0151] In some embodiments, the updating module 705 can be specifically used to obtain a gradient based on the zero-delay cross-correlation of the reconstructed source wave field and the back-propagation residual wave field. The initial geological velocity model is iteratively updated using a gradient-based optimization algorithm until the final geological velocity model is output, including: calculating the back-propagation residual wave field based on the difference between the actual seismic observation data and the seismic simulation data; performing cross-correlation based on the back-propagation residual wave field and the reconstructed source wave field to obtain the objective function gradient; and updating the initial geological velocity model using a gradient-based optimization algorithm based on the obtained objective function gradient to obtain the final velocity model.

[0152] In summary, the acoustic full waveform inversion method and system based on the optimized effective boundary storage in the embodiment of the present invention directly adopts the first-order stress-velocity acoustic wave equation for forward simulation, and obtains the value range of the sampling time step according to the Nyquist theorem, resamples the effective boundary wave field in the time domain, and then uses the interpolation method to restore the resampled wave field data to the original data state, and then reconstructs the source wave field, and performs zero-delay cross-correlation between the reconstructed source wave field and the back-transmitted residual wave field to obtain the gradient. It can effectively reduce the storage amount of data and speed up the transmission efficiency of data.

[0153] The embodiments described above are only descriptions of the preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope determined by the claims of the present invention.

Claims

1. A method for full-waveform inversion of acoustic waves based on optimized effective boundary storage, characterized in that the steps include: Construct an initial geological velocity model and obtain seismic wave parameter information; Perform forward modeling based on the initial geological velocity model and the seismic wave parameter information to obtain seismic simulation data and forward source wave field data; Acquire actual seismic observation data of the geological model to be tested, and construct an objective function based on the actual seismic observation data and the seismic simulation data; Resampling the target area wavefield data according to the forward source wavefield data; restoring the resampled target area wavefield data to the original data state by using an interpolation method, and then reconstructing the forward source wavefield; The gradient is obtained by performing zero-delay cross-correlation based on the reconstructed forward source wave field and the reverse residual wave field; the objective function is iteratively updated using a gradient optimization algorithm until a final geological velocity model is output; the final geological velocity model is the inversion result.

2. The acoustic full waveform inversion method based on optimized effective boundary storage according to claim 1 is characterized in that: The method for performing the forward modeling includes: performing a first-order stress-acoustic wave equation according to the initial geological velocity model and the seismic wave parameter information, and the expression includes: Among them, P represents stress; Vx and Vz represent the particle velocity components in the x and z directions respectively; v represents the longitudinal wave velocity; s represents the earthquake source; ρ represents the density; t represents the time; x and z represent the grid x and z coordinate values.

3. The acoustic full waveform inversion method based on optimized effective boundary storage according to claim 1 is characterized in that: According to the actual earthquake observation data and the earthquake simulation data, the objective function constructed is the residual two-norm objective function: Where D represents the observation target area; u syn and u obs Represent earthquake simulation data and actual earthquake observation data respectively; δ(xx r ) represents the unit pulse function; T max represents the maximum time earthquake record; m represents the velocity model parameter; x r Indicates the location of the detector.

4. The acoustic full waveform inversion method based on optimized effective boundary storage according to claim 1 is characterized in that: When resampling the wavefield data of the target area according to the forward source wavefield, the sampling interval formula is obtained according to the Nyquist sampling theorem, and the expression includes: Where T is the sampling interval obtained according to the Nyquist sampling theorem, f max is the maximum frequency of the current sampling; The interpolation method is used to restore the resampled wave field data to the original data state, and the steps include: Where t is time, P(x,t) represents the wavefield information at a specific spatial position x and time t through the interpolation method; P(x,t0) is the resampled wavefield data at time t0, which serves as the left endpoint of the interpolation process; P(x,t1) is the wavefield data that has also been resampled at time t1, which serves as the right endpoint of the interpolation process; x(x,z) is the spatial position in the x and z directions.

5. The acoustic full waveform inversion method based on optimized effective boundary storage according to claim 3 is characterized in that: Calculating the back-propagation residual wave field according to the difference between the actual seismic observation data and the seismic simulation data, the steps include: According to the Lagrange principle, a first-order acoustic wave equation constraint is imposed on the residual two-norm objective function: Add the constraint term to the residual two-norm objective function and change the objective function to: Wherein, λ(x,t) represents the accompanying wave field vector, also known as the Lagrange multiplier; E(m)=0 is the above-mentioned first-order acoustic wave equation; The improved objective function is differentiated to obtain the adjoint equation, which is the reverse time propagation equation of the residual wave field with the source as the source. The formula includes: Among them, S λp represents the accompanying source, which is the difference between the earthquake simulation record and the actual earthquake observation record; Vx represents the Lagrange multiplier in the x direction, λ Vz represents the Lagrange multiplier in the z direction, λ P Denotes the companion wave field.

6. The acoustic full waveform inversion method based on optimized effective boundary storage according to claim 5 is characterized in that: The objective function gradient is obtained by cross-correlating the back-propagated residual wave field and the reconstructed source wave field, including: The gradient calculation formula of velocity v is: in, represents the gradient of H2 function with respect to v; represents the gradient of the H2 function with respect to the bulk modulus.

7. The acoustic full waveform inversion method based on optimized effective boundary storage according to claim 5 is characterized in that: According to the obtained gradient, the initial geological velocity model is updated using a gradient optimization algorithm to obtain a final velocity model, including: Determining the updating direction of the geological velocity model according to the obtained gradient; According to the updating direction of the geological velocity model, the initial geological velocity model is updated to obtain a final velocity model.

8. An acoustic full waveform inversion system based on optimized effective boundary storage, the system is used to implement the method according to any one of claims 1 to 7, characterized in that: include: A construction module, a first calculation module, a second calculation module, a third calculation module and an update module; The construction module is used to construct an initial geological velocity model and obtain seismic wave parameter information; The first calculation module is used to perform forward simulation according to the initial geological velocity model and the seismic wave parameter information to obtain seismic simulation data and forward source wave field data; The second calculation module is used to obtain actual seismic observation data of the geological model to be measured, and to construct an objective function according to the actual seismic observation data and the seismic simulation data; The third calculation module is used to resample the target area wavefield data according to the forward source wavefield data; restore the resampled target area wavefield data to the original data state by using an interpolation method, and then reconstruct the forward source wavefield; The updating module is used to obtain the gradient by performing zero-delay cross-correlation based on the reconstructed forward source wave field and the reverse residual wave field; the objective function is iteratively updated using a gradient-based optimization algorithm until a final geological velocity model is output; the final geological velocity model is the inversion result.

Images