A method and system for acoustic full waveform inversion based on optimized effective boundary storage

Abstract

The present invention discloses an acoustic full waveform inversion method and system based on optimized effective boundary storage, the method comprising the following steps: constructing an initial geological velocity model and obtaining seismic wave parameter information; performing forward simulation according to the initial geological velocity model and the seismic wave parameter information to obtain seismic simulation data and forward source wavefield data; obtaining actual seismic observation data of the geological model to be measured, and constructing a target function according to 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; performing zero-delay cross-correlation on the reconstructed forward source wavefield and the backward residual wavefield to obtain a gradient; and iteratively updating the target function by using a gradient-based optimization algorithm until a final geological velocity model is output.

Description

Technical Field

[0001] The present 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 its introduction by Laily (1983) and Tarantola (1984) in the early 1980s, full waveform inversion (FWI) theory has been extensively studied in various fields, including oil and gas and mineral resource exploration and development. FWI uses a series of iterative processes to continuously update the subsurface velocity model, aiming to minimize the discrepancy between seismic simulation data and actual seismic observation data. Through repeated iterations, the simulated seismic wavefield gradually approaches the actual seismic observation wavefield, enabling high-precision reconstruction of true subsurface physical parameters (such as wave velocity and density). Compared to traditional seismic imaging methods, FWI can more precisely characterize subsurface structural features. However, this method has significant 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 during the full waveform inversion process, thereby improving the inversion efficiency.

[0004] The purpose of the embodiments 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 simulation based on the initial geological velocity model and the seismic wave parameter information to obtain seismic simulation data and forward source wavefield data;

[0007] Acquiring actual seismic observation data of the geological model to be measured, and constructing an objective function based on the actual seismic observation data and the seismic simulation data;

[0008] Resampling the target area wavefield data based on the forward source wavefield data; restoring the resampled target area wavefield data to the original data state 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 wavefield and the reverse residual wavefield; 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.

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

[0011]

[0012] Where P is stress, Vx and Vz are the particle velocity components in the x and z directions, respectively, v is the longitudinal wave velocity, s is the earthquake source, ρ is the density, t is the time, and x and z are the grid x and z coordinates.

[0013] Preferably, the objective function constructed based on 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 position of the detector.

[0016] Preferably, when resampling the target area wavefield data according to the forward source wavefield, a 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. The steps include:

[0020]

[0021] Where t is time, and P(x,t) represents the wavefield information at a specific spatial location x and time t, obtained through interpolation. P(x,t0) is the resampled wavefield data at time t0, serving as the left endpoint of the interpolation process; P(x,t1) is the resampled wavefield data at time t1, serving as the right endpoint of the interpolation process; and x(x,z) represents 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, and the steps include:

[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] Where λ(x, t) represents the adjoint wave field vector, also known as the Lagrange multiplier; E(m) = 0 is the first-order acoustic wave equation mentioned above;

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

[0030]

[0031] 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 Indicates the accompanying 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 the 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 an update direction of the geological velocity model based on the obtained gradient;

[0038] The initial geological velocity model is updated according to the updating direction of the geological velocity model 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 based on the initial geological velocity model and the seismic wave parameter information to obtain seismic simulation data and forward source wavefield 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 based on the actual seismic observation data and the seismic simulation data;

[0043] The third calculation module is used to resample the target area wavefield data based on 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;

[0044] The updating module is used to obtain the gradient by performing zero-delay cross-correlation based on the reconstructed forward source wavefield and the reverse residual wavefield; 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 time-domain acoustic full-waveform inversion, a resampling algorithm based on the Nyquist sampling theorem is applied to efficient boundary storage, significantly reducing the storage capacity required for wavefield data. The resampled wavefield data is restored using linear interpolation to complete source wavefield reconstruction and integrate it into the time-domain full-waveform inversion process. This solution significantly reduces data storage pressure and accelerates full-waveform inversion computations. 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 any creative work.

[0048] Figure 1 This 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 This 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 showing a comparison between original data and interpolated data of wavefield data in an embodiment of the present invention when the sampling interval is 2.5 ms;

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

[0052] Figure 3 c is a schematic diagram showing the comparison between the original data and the interpolated data of the wavefield data of 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 wavefield data in an embodiment of the present invention at a sampling interval of 2.5 ms;

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

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

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

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

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

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

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

[0061] Figure 5c is a schematic diagram showing a comparison of relative errors between a conventional wavefield reconstruction image at 500ms and a wavefield reconstruction image after resampling according to 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 structural diagram 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 clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. 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 making creative efforts 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] Example 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 measured and seismic wave parameter information are input; forward simulation is performed based on the initial geological velocity model and seismic wave parameter information to obtain seismic simulation data and forward source wavefield data; an objective function is constructed based on the actual seismic observation data and seismic simulation data; the wavefield data of the target area is resampled based on the forward source wavefield; the resampled wavefield data is restored to the original data state using an interpolation method, and then the source wavefield is reconstructed; zero-delay cross-correlation is performed based on the reconstructed source wavefield and the back-transmitted residual wavefield to obtain the gradient. A gradient 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, an 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 wavefield data;

[0073] Step 103: Acquire actual earthquake observation data of the geological model to be measured, 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 based on 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: Perform zero-delay cross-correlation based on the reconstructed forward source wavefield and the reverse residual wavefield to obtain the gradient; use a gradient optimization algorithm to iteratively update the objective function until the final geological velocity model is output; the final geological velocity model is the inversion result.

[0076] Depend on Figure 2 As shown in the illustrated process, the acoustic full-waveform inversion method for optimizing effective boundary storage in this embodiment of the present invention directly employs the first-order stress-velocity acoustic wave equation for forward modeling. A resampling algorithm based on the Nyquist sampling theorem is applied to the effective boundary storage, significantly reducing the storage capacity required for wavefield data. The resampled wavefield data is restored using linear interpolation to complete source wavefield reconstruction and integrated into the time-domain full-waveform inversion process, significantly reducing data storage pressure and accelerating the computational 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 constructed based on prior knowledge of seismic exploration, assumptions about geological structures, and possible geophysical data to form a preliminary geological velocity model for subsequent iterative optimization. This model should be able to reflect the general velocity distribution characteristics of the underground medium.

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

[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 wavefield 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] Where 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 density; and t represents 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 obtained 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 dinorm is a quantitative measure of this matching degree. It is calculated by calculating the temporal and spatial differences between the simulated and observed wavefields, taking the square root of the sum of these differences. A small residual dinorm indicates that the simulated seismic wavefield closely matches the observed wavefield, and the model fits the actual situation well. Conversely, a large residual dinorm indicates a significant deviation between the model and the observed wavefield, necessitating an optimization algorithm to update the velocity model and obtain a simulated seismic wavefield that is closer to the observed wavefield.

[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 position of the detector.

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

[0092] In some embodiments, a formula for obtaining a reasonable sampling interval according to the Nyquist sampling theorem is 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 location x and time t, obtained through interpolation. P(x,t0) is the resampled wavefield data at time t0, serving as the left endpoint of the interpolation process. P(x,t1) is the similarly resampled wavefield data at time t1, serving as the right endpoint of the interpolation process. Wavefield values ​​at any location between the two are thus determined. x(x,z) represents 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 wavefield and the reverse residual wavefield. 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, a gradient is obtained by performing zero-delay cross-correlation based on the reconstructed source wavefield and the back-propagated residual wavefield. A gradient-based optimization algorithm is used to iteratively update the initial geological velocity model until the final geological velocity model is output, specifically including:

[0100] Calculate the backpropagation residual wave field based on the difference between actual earthquake observation data and earthquake simulation data;

[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 Indicates the accompanying wave field.

[0111] Expanding Formula 7 yields the following form:

[0112]

[0113] To obtain the minimum value of Formula 7, the following three optimization conditions must 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 Indicates the accompanying wave field.

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

[0119]

[0120] The gradient calculation formula about 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 obtain that the H2 function is the gradient of v, which can be 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 computational complexity of wavefield data at all times during the entire time period is enormous. Equations 3 and 4 show that the effective boundary wavefield resampling algorithm, based on the Nyquist sampling theorem, significantly reduces the storage capacity required for wavefield data. Resampling the resampled wavefield data using linear interpolation effectively reduces computational complexity and storage requirements, thereby improving the computational speed of 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 update direction of the geological velocity model based on 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, updating the initial geological velocity model according to the update direction of the geological velocity model specifically includes:

[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 the velocity model parameters.

[0136] In some embodiments, to verify the effectiveness of the interpolation method, a uniform velocity model is designed with a length and depth of 0.804 km, a horizontal and vertical grid spacing of 4 meters, and 22 absorbing boundaries. The wave field propagation in this model is simulated first with a duration of 0.5 seconds and a time interval of 0.5 milliseconds. The wavelet used is a Ricker wavelet with a main frequency of 10 Hz. According to the Nyquist sampling theorem, the maximum time sampling interval can be obtained to be 10.5 milliseconds. Figure 3 a is a schematic diagram showing a comparison between original data and interpolated data of wavefield data in an embodiment of the present invention when the sampling interval is 2.5 ms; Figure 3 b is a schematic diagram showing the comparison between the original data and the interpolated data of the wavefield data of an embodiment of the present invention when the sampling interval is 5 ms; Figure 3 c is a schematic diagram showing the comparison between the original data and the interpolated data of the wavefield data of 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 wavefield data in an embodiment of the present invention at a sampling interval of 2.5 ms; Figure 3 e is a schematic diagram showing a comparison of relative errors between original data and interpolated data of wavefield data in an embodiment of the present invention at a sampling interval of 5 ms; Figure 3 f is a schematic diagram showing 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 are 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 As can be seen from f, a 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 according to an embodiment of the present invention; Figure 4 b is a wave field reconstructed image after resampling at 250ms in an embodiment of the present invention; Figure 4 c is a schematic diagram showing the relative error comparison between the conventional wavefield reconstruction image at 250ms and the wavefield reconstruction after resampling according to the embodiment of the present invention; Figure 4 a to Figure 4As can be seen from Figure c, by comparing the conventional wavefield reconstruction at 250ms with the wavefield reconstruction image after resampling, it is found that the reconstructed images before and after sampling are highly consistent in morphology. To further understand the difference between the two, a relative error calculation was performed, 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 an embodiment of the present invention; Figure 5 b is a wave field reconstructed image after resampling at 500ms in an embodiment of the present invention; Figure 5 c is a schematic diagram showing the relative error comparison between the conventional wavefield reconstruction image at 500ms and the wavefield reconstruction after resampling in the embodiment of the present invention; Figure 5 a to Figure 5 As can be seen from Figure c, by comparing the conventional wavefield reconstruction at 500ms with the wavefield reconstruction image after resampling, it is found that the reconstructed images before and after sampling are highly consistent in morphology. To further understand the difference between the two, a relative error calculation was performed, 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, which 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 6The initial velocity model is a full-waveform inversion result diagram based on linear interpolation provided by this application. The model has a length of 6.528 km, a depth of 1.544 km, and a total simulation time of 2.8 s. The horizontal and vertical grid spacing of the velocity model is 8 m, the absorbing boundary is 22 layers, the time interval is 0.8 ms, and the main frequency settings are 2.5 Hz, 5.5 Hz, 12 Hz, 16 Hz, and 20 Hz, respectively. The maximum time sampling intervals are 43.2 ms, 19.2 ms, 8.8 ms, 6.4 ms, and 4.8 ms, respectively. The effective boundary wave field data after resampling is restored by linear interpolation and wave field reconstruction is performed, and the reconstructed wave field is integrated into the time domain acoustic full-waveform inversion process. This inversion process is iteratively optimized, gradually approximating the initial model, and the final inversion result is in good agreement 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, this application provides an acoustic full-waveform inversion method based on optimized effective boundary storage. This method employs a first-order stress-velocity acoustic wave equation and, based on the Nyquist theorem, derives the range of sampling time steps. The effective boundary wavefield is then resampled in the time domain. Interpolation is then used to restore the resampled wavefield data to its original state. The source wavefield is then reconstructed, and the gradient is derived by performing zero-delay cross-correlation between the reconstructed source wavefield and the back-propagated residual wavefield. This method can effectively reduce data storage and improve data transmission efficiency.

[0141] Example 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 based on the initial geological velocity model and seismic wave parameter information to obtain seismic simulation data and forward source wavefield data;

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

[0146] The third calculation module 704 is used to resample the target area wavefield data based on 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;

[0147] Update module 705 is used to obtain the gradient by performing zero-delay cross-correlation based on the reconstructed forward source wavefield and the reverse residual wavefield; it uses a gradient optimization algorithm to iteratively update the objective function 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 may be specifically configured to perform forward simulation based on the initial geological velocity model and seismic wave parameter information to obtain seismic simulation data and forward source wavefield data.

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

[0150] In some embodiments, the third calculation module 704 can be specifically used to resample the target area wavefield data based on the forward source wavefield; use the interpolation method to restore the resampled wavefield data to the original data state, and then reconstruct the source wavefield, including: finding a reasonable sampling interval according to the Nyquist sampling theorem, resampling the effective boundary wavefield in the time domain, and then interpolating and restoring the resampled boundary wavefield values ​​according to the linear interpolation method.

[0151] In some embodiments, the updating module 705 can be specifically configured to perform zero-delay cross-correlation on the reconstructed source wavefield and the back-propagated residual wavefield to obtain a gradient. The initial geological velocity model is iteratively updated using a gradient-based optimization algorithm until a final geological velocity model is output. The update module 705 includes the following steps: calculating the back-propagated residual wavefield based on the difference between actual seismic observation data and seismic simulation data; performing cross-correlation on the back-propagated residual wavefield and the reconstructed source wavefield to obtain a target function gradient; and updating the initial geological velocity model using a gradient-based optimization algorithm based on the obtained target function gradient to obtain a final velocity model.

[0152] In summary, the acoustic full-waveform inversion method and system based on optimized effective boundary storage in the embodiments of the present invention directly uses the first-order stress-velocity acoustic wave equation for forward modeling. Based on the Nyquist theorem, the sampling time step range is determined. The effective boundary wavefield is resampled in the time domain. Subsequently, the resampled wavefield data is restored to its original state using interpolation. The source wavefield is then reconstructed, and the gradient is derived by performing zero-delay cross-correlation between the reconstructed source wavefield and the back-propagated residual wavefield. This effectively reduces data storage and improves data transmission efficiency.

[0153] The embodiments described above are merely descriptions of preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by persons skilled in the art should fall within the scope of protection defined by the claims of the present invention.

Claims

1. A full-waveform inversion method for 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 simulation based on the initial geological velocity model and the seismic wave parameter information to obtain seismic simulation data and forward source wavefield data; Acquiring actual seismic observation data of the geological model to be measured, and constructing an objective function based on the actual seismic observation data and the seismic simulation data; Resampling the target area wavefield data based on the forward source wavefield data; restoring the resampled target area wavefield data to the original data state using an interpolation method, and then reconstructing the forward source wavefield; 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: 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. The steps include: Where t is time, P(x,t) represents the wavefield information at a specific spatial position x and time t through interpolation; 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 resampled wavefield data 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; The gradient is obtained by performing zero-delay cross-correlation based on the reconstructed forward source wavefield and the reverse residual wavefield; 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.

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

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 position of the detector.

4. The acoustic full waveform inversion method based on optimized effective boundary storage according to claim 3 is characterized in that: Calculating the backpropagation residual wavefield according to the difference between the actual seismic observation data and the seismic simulation data includes: 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: Where λ(x, t) represents the adjoint wave field vector, also known as the Lagrange multiplier; E(m) = 0 is the first-order acoustic wave equation mentioned above; The improved objective function is differentiated to obtain the adjoint equation, which is the reverse time propagation equation with the source as the residual wave field. 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 Indicates the accompanying wave field.

5. The acoustic full waveform inversion method based on optimized effective boundary storage according to claim 4 is characterized in that: The objective function gradient is obtained by performing cross-correlation based on the back propagation residual wavefield and the reconstructed forward propagation source wavefield, including: The gradient calculation formula of velocity v is: 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.

6. The acoustic full waveform inversion method based on optimized effective boundary storage according to claim 4 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 an update direction of the geological velocity model based on the obtained gradient; The initial geological velocity model is updated according to the updating direction of the geological velocity model to obtain a final velocity model.

7. An acoustic full waveform inversion system based on optimized effective boundary storage, the system being used to implement the method according to any one of claims 1 to 6, 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 based on the initial geological velocity model and the seismic wave parameter information to obtain seismic simulation data and forward source wavefield 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 based on the actual seismic observation data and the seismic simulation data; The third calculation module is used to resample the target area wavefield data based on 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; The updating module is used to obtain the gradient by performing zero-delay cross-correlation based on the reconstructed forward source wavefield and the reverse residual wavefield; 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