阅读说明：本技术 基于频域的最小二乘共轭梯度迭代的地震正演模拟方法 (Frequency domain-based seismic forward modeling method based on least square conjugate gradient iteration ) 是由 刘文革 周觅路 尹成 牟其松 涂文茂 李崇文 于 2021-08-18 设计创作，主要内容包括：本发明公开了基于频域的最小二乘共轭梯度迭代的地震正演模拟方法,所述方法包括：输入经傅里叶变换至频率域后的时谐波动方程；设置相应的正演模拟参数；根据所述正演模拟参数,通过最小二乘共轭梯度算法对所述时谐波动方程进行迭代求解。本发明的方法相较于基于LU直接求解的正演方法,可在保持同等精度的条件下减少内存消耗、提高计算效率,相较于基于BI-CGSTAB迭代求解的正演方法,可在计算效率接近的情况下获得更宽松的使用条件和更好的计算稳定性。(The invention discloses a frequency domain-based seismic forward modeling method based on least square conjugate gradient iteration, which comprises the following steps: inputting a time harmonic motion equation after Fourier transformation to a frequency domain; setting corresponding forward modeling parameters; and according to the forward modeling parameters, carrying out iterative solution on the time harmonic motion equation by a least square conjugate gradient algorithm. Compared with a forward modeling method based on LU direct solution, the method provided by the invention can reduce the memory consumption and improve the calculation efficiency under the condition of keeping the same precision, and compared with a forward modeling method based on BI-CGSTAB iterative solution, the method provided by the invention can obtain more relaxed use conditions and better calculation stability under the condition that the calculation efficiency is close to that.)

1. The seismic forward modeling method based on least square conjugate gradient iteration of a frequency domain is characterized by comprising the following steps of:

step S1: inputting a time harmonic motion equation after Fourier transformation to a frequency domain;

step S2: setting corresponding forward modeling parameters;

step S3: and according to the forward modeling parameters, carrying out iterative solution on the time harmonic motion equation by a least square conjugate gradient algorithm.

2. The simulation method of claim 1, wherein: it still includes:

and when the calculation error of the target parameter obtained by the iterative solution converges to the range of the preset error, performing forward simulation on the obtained target parameter, and outputting a simulation result.

3. The simulation method of claim 1, wherein: the iterative solution process comprises:

discretizing the time harmonic motion equation by a difference method to obtain a linear system model Ap ═ b (1), wherein A represents a coefficient matrix, p represents a pressure field vector, b represents a discrete form of a seismic source term in the time harmonic motion equation, and b ═ S (omega) delta (x-x)_s)δ(z-z_s) (2) where S (ω) represents a source term in the frequency domain after Fourier transform, δ () represents a Dirac function, x represents an abscissa argument, z represents an ordinate argument, x represents a frequency domain after Fourier transform, and_srepresenting the source abscissa, z_sRepresenting a source ordinate;

carrying out diagonal preprocessing on the linear system model through a preprocessing matrix M with the same sparse structure as the coefficient matrix A to obtain an equivalent linear system model (AM) of the linear system model^-1) v ═ b, v ═ Mp (3), where M denotes the pre-processing matrix and v denotes the intermediate vector;

setting an iteration initial value x₀And a preset error epsilon;

and solving the equivalent linear system model through the least square conjugate gradient iterative algorithm based on the iteration initial value until the iterative calculation error converges to the preset error epsilon, obtaining a result of a parameter to be calculated, namely a target parameter, and performing forward modeling according to the result.

4. The simulation method according to claim 3, wherein: the solving model of the least square conjugate gradient iterative algorithm is as follows:

wherein:

s_k＝α_kd_k (10)；

where k denotes the number of iterations, θ_kThe least squares coefficients for the kth iteration are represented,represents the least squares conjugate gradient parameters at the kth iteration,representing the first conjugate gradient parameter at the k-th iteration,representing the second conjugate gradient parameter at the kth iteration,representing the target parameter x obtained after the kth iteration_kA gradient matrix ofDenotes x_k+1The gradient matrix g_k+1Transpose of (y)_kRepresenting the difference in gradient in two iterations, i.e. y_k＝g_k+1-g_k，Representing the search direction matrix d at the kth iteration_kIs transposed, andi.e. the initial search direction is the initial value x₀Direction of negative gradient of (d), search direction thereafter passing through the parameterAnd the current position is obtained, and,a first search direction matrix representing least squares conjugate gradient parameters at the (k + 1) th iteration,a second search direction matrix, s, representing the least-squares conjugate gradient parameters at the (k + 1) th iteration_kRepresenting the iteration step size alpha_kIn the search direction d_kProjection of (2).

5. The simulation method of claim 4, wherein: the solving process based on the solving model comprises the following steps:

the iteration initial value x of the target parameter₀Carry-over (4) to obtain the search direction d of the first iteration₀And an iteration step size alpha₀；

Along the search direction d of the first iteration₀Using iteration step size alpha₀Performing a first iterative computation on the equivalent linear system model to obtain a first iterative result x₁；

The first iteration result x₁Bringing inCalculating the current iteration error epsilon by the following error calculation model₁：

Judging the current iteration error epsilon₁Whether the error is smaller than a preset iteration error epsilon or not; if the sum is less than the preset threshold value, outputting a solution matrix for subsequent forward modeling; if yes, the first iteration result x is used₁Repeating the solving process as a new iteration initial value;

repeating the solving process, and obtaining a new iterative searching direction d each time_kAnd an iteration step size alpha_kCalculating a new iteration result x of a new round using a new search direction and a new iteration step size_kAnd a new iteration error epsilon_kUp to epsilon_kJumping out of iteration when the epsilon is less than epsilon to obtain an approximate solution of a linear system;

and (3) participating in forward modeling by using the linear system approximate solution.

6. The simulation method according to claim 3, wherein: the difference method is a nine-point finite difference method.

7. The simulation method according to claim 3, wherein: the preset error epsilon is 0.01%.

8. The simulation method of claim 1, wherein: the forward modeling parameters comprise iteration initial values x₀And a preset error epsilon, and one or more of the following parameters: under a model coordinate system, the space interval dx in the x direction, the space interval dz in the z direction, the space position coordinates of the seismic source, the main frequency MF of the seismic source wavelet, the forward modeling time T, the number of frequency slices, the frequency interval delta F, the number of PML absorption layers and the absorption coefficient.

9. The simulation method of claim 1, wherein: the speed model is selected from one or more of a homogeneous medium model, a layered medium model, a Marmousi model and an Overthrast model.

10. The simulation method of claim 9, wherein: the Marmousi model is selected from a standard Marmousi-1 model or an improved Marmousi model obtained after resampling is carried out on the standard Marmousi model-1.

Technical Field

The invention relates to the technical field of seismic exploration.

Background

Seismic exploration, a basic geophysical exploration method, is the most effective method for solving the problem of oil and gas exploration. Forward modeling based on the wave equation is fundamental work of seismic exploration and has important effects in the stages of acquisition, processing, explanation and the like of the seismic exploration.

The seismic forward modeling method in the prior art can be divided into a time domain and a frequency domain according to a calculation domain, wherein the time domain algorithm is widely applied due to simple implementation and high precision, but under the requirement of full waveform inversion, a wave equation needs to be repeatedly solved by iteration on the basis of forward modeling, and if the forward modeling is carried out in the time domain, huge calculation amount is brought, and a large amount of memory and time are consumed.

Compared with time domain forward modeling, frequency domain forward modeling has the advantages of no time dispersion, no stability problem, no error accumulation, flexible frequency band selection and the like, a time domain wave equation is called as a time harmonic equation after being subjected to Fourier transform to a frequency domain, the solution of the equation is an important ring of the frequency domain forward modeling, and the currently common time harmonic equation solution method mainly comprises two main types: the direct solving method has high precision, large memory consumption and low calculation efficiency; the iteration method comprises the following steps: the accuracy is relatively low, but the memory consumption is low, the calculation efficiency is high, and the direct solution under the three-dimensional condition is very difficult, so the use advantage of the iterative method is more obvious. In order to overcome the defects of an iterative method, different forms of iterative solution methods are proposed in the prior art, most of the iterative solution methods are based on a Krylov subspace method, such as a CG method and a Bi-CGSTAB method, but the CG method is low in convergence speed and has insignificant advantages; the Bi-CGSTAB method has high convergence speed, but iteration is easy to be not converged, and harsh use conditions are required.

For example, in a direct solving method based on LU decomposition proposed in the prior patent document CN103901472A, a difference formula of a 17-point format is established, a sparse matrix is constructed, then wavelet parameters and a velocity model are read in, a single-frequency wave field is obtained through frequency cycle, and finally inverse fourier transform is performed on all the frequency wave fields to obtain a forward result. Through the method and the device provided by the embodiment of the invention, the four-order high-precision forward modeling overcomes the limitation that dx must be equal to dz in the past, can be applied to dx ═ dz and also can be applied to the situation that dx ≠ dz, and the number of grids required by the minimum wavelength only needs 2.4 in the method, which is superior to other similar methods.

Or, as disclosed in the prior art document "Frequency-space domain electronic wave simulation with the BiCGstab, iterative method [ J ]" (Journal of geometics & Engineering,2016,13(1):70-77Du Z, Liu J, Liu W, et al.) the iterative method of using Bi-CGSTAB to solve large sparse linear equations greatly reduces memory requirements and computation time as the spatial differential order increases. Compared with LU decomposition, the calculation efficiency of the iteration of the Bi-CGSTAB is improved by 9 times when fourth-order differentiation is adopted, and the numerical simulation result of the iteration is equivalent to the LU decomposition result, but although the memory consumption and the calculation time are obviously reduced, the use condition of the method is harsh, and the unstable calculation condition is easy to occur.

Disclosure of Invention

The invention aims to provide a frequency domain-based least square conjugate gradient iteration seismic forward modeling method, which can reduce memory consumption and improve calculation efficiency under the condition of keeping equivalent precision compared with an LU direct solution-based forward modeling method, and can obtain more relaxed use conditions and better calculation stability under the condition that the calculation efficiency is close to that of a BI-CGSTAB iterative solution-based forward modeling method.

The technical scheme of the invention is as follows:

a seismic forward modeling method based on least square conjugate gradient iteration of a frequency domain comprises the following steps:

1. inputting a time harmonic motion equation after Fourier transformation to a frequency domain;

2. setting corresponding forward modeling parameters;

3. and according to the forward modeling parameters, carrying out iterative solution on the time harmonic motion equation by a least square conjugate gradient algorithm.

According to some embodiments of the invention, the simulation method further comprises:

and when the calculation error of the target parameter obtained by the iterative solution converges to the range of the preset error, performing forward simulation on the obtained target parameter, and outputting a simulation result.

According to some embodiments of the invention, the iterative solution comprises:

discretizing the time harmonic motion equation by a difference method to obtain a linear system model Ap ═ b (1), wherein A represents a coefficient matrix, p represents a pressure field vector, b represents a discrete form of a seismic source term in the time harmonic motion equation, and b ═ S (omega) delta (x-x)_s)δ(z-z_s) (2), wherein S (ω) represents a source term of the frequency domain after Fourier transform, δ () represents a Dirac function, x represents an abscissa argument, z represents an ordinate argument, x represents a frequency domain after Fourier transform, and_srepresenting the source abscissa, z_sRepresenting a source ordinate;

carrying out diagonal preprocessing on the linear system model through a preprocessing matrix M with the same sparse structure as the coefficient matrix A to obtain an equivalent linear system model (AM) of the linear system model^-1) v ═ b, v ═ Mp (3) where M denotes the pre-processing matrix and v denotes the intermediate vector;

setting an iteration initial value x of a target parameter x to be calculated₀And a preset error epsilon;

and solving the equivalent linear system model through the least square conjugate gradient iterative algorithm based on the iteration initial value until the iterative calculation error converges to the preset error epsilon, obtaining a result of a parameter to be calculated, namely a target parameter, and performing forward modeling according to the result.

According to some embodiments of the invention, the solution model of the least squares conjugate gradient iterative algorithm is as follows:

wherein:

s_k＝α_kd_k(10)；

where k denotes the number of iterations, θ_kThe least squares coefficients for the kth iteration are represented,represents the least squares conjugate gradient parameters at the kth iteration,representing the first conjugate gradient parameter at the k-th iteration,representing the second conjugate gradient parameter at the kth iteration, g_k＝▽f(x_k) Representing the target parameter x obtained after the kth iteration_kA gradient matrix ofDenotes x_k+1The gradient matrix g_k+1Transpose of (y)_kRepresenting the difference in gradient in two iterations, i.e. y_k＝g_k+1-g_k，Representing the search direction matrix d at the kth iteration_kIs transposed, andi.e. the initial search direction is the initial value x₀Direction of negative gradient of (d), search direction thereafter passing through the parameterAnd the current position is obtained, and,a first search direction matrix representing least squares conjugate gradient parameters at the (k + 1) th iteration,a second search direction matrix, s, representing the least-squares conjugate gradient parameters at the (k + 1) th iteration_kRepresenting the iteration step size alpha_kIn the search direction d_kProjection of (2).

According to some embodiments of the invention, the solution process based on the solution model comprises:

the iteration initial value x of the target parameter₀Carry-over (4) to obtain the search direction d of the first iteration₀And an iteration step size alpha₀；

Along the search direction d of the first iteration₀Using iteration step size alpha₀Performing a first iterative computation on the equivalent linear system model to obtain a first iterative result x₁；

The first iteration result x₁Calculating the current iteration error epsilon by the error calculation model₁：

Judging the current iteration error epsilon₁Whether the error is smaller than a preset iteration error epsilon or not; if the sum is less than the preset threshold value, outputting a solution matrix for subsequent forward modeling; if it is largeThen, the first iteration result x is obtained₁Repeating the solving process as a new iteration initial value;

repeating the solving process, and obtaining a new iterative searching direction d each time_kAnd an iteration step size alpha_kCalculating a new iteration result x of a new round using a new search direction and a new iteration step size_kAnd a new iteration error epsilon_kUp to epsilon^kJumping out of iteration when the epsilon is less than epsilon to obtain an approximate solution of a linear system;

and (3) participating in forward modeling by using the linear system approximate solution.

According to some embodiments of the invention, the difference method is selected from a nine point finite difference method.

According to some embodiments of the invention, the predetermined error e is 0.01%.

According to some embodiments of the invention, the forward simulation result comprises one or more of a single frequency data slice, a synthetic wave field snapshot, and a seismic single shot record.

According to some embodiments of the invention, the forward modeling parameters include an iteration initial value x₀And a preset error epsilon, and one or more of the following parameters: under a model coordinate system, the space interval dx in the x direction, the space interval dz in the z direction, the space position coordinates of the seismic source, the main frequency MF of the seismic source wavelet, the forward modeling time T, the number of frequency slices, the frequency interval delta F, the number of PML absorption layers and the absorption coefficient.

According to some embodiments of the invention, the velocity model is selected from one or more of a homogeneous medium model, a stratified medium model, a Marmousi model, an Overthrust model.

According to some embodiments of the invention, the Marmousi model is selected from a raw Marmousi model or a modified Marmousi model after resampling of the raw Marmousi model.

The invention has the following beneficial effects:

the frequency domain forward modeling method based on least square conjugate gradient iteration can quickly and efficiently complete forward modeling of a frequency domain, generate single-frequency slices, wave field snapshots or seismic records, and provide important help for waveform inversion.

The forward modeling method of the invention has the advantages of loose iteration condition, stable calculation and higher efficiency under the condition of ensuring equivalent precision.

The forward modeling method has reliable seismic response simulation characteristics, and compared with the LU direct method, the forward modeling method has the advantages that the calculation memory can be reduced by 15-20% under the calculation of a single frequency slice, the calculation time can be reduced by 58-88%, the calculation result is equivalent in accuracy, and the forward modeling method has excellent performance under the condition that the iteration error is 0.01%.

The simulation precision of the simulation method can be flexibly set according to an actual speed model, a computing hardware environment and the like, the size of the model, the number of frequency slices and the like can be set according to the memory environment of computing equipment, and the method is high in applicability and convenient and rapid to apply.

Drawings

FIG. 1 is a flow chart of a particular seismic modeling forward modeling method;

FIG. 2 is a diagram of a specific coefficient matrix structure;

FIG. 3 is a schematic view of a specific layer model;

FIG. 4 is a schematic diagram of a specific Marmous-1 model;

FIG. 5 is a diagram showing the forward modeling result of the layered model at 30Hz data slice;

FIG. 6 is a diagram illustrating the results of the layered model synthesizing a wave field snapshot at 312 ms;

FIG. 7 is a diagram showing forward modeling results of the Marmousi-1 model at 15Hz data slice earthquake;

FIG. 8 is a diagram showing the Marmousi-1 model synthesizing a wave field snapshot result in 625 ms;

FIG. 9 is a diagram showing the Marmousi-1 model seismic single shot record comparison under two methods.

Detailed Description

The present invention is described in detail below with reference to the following embodiments and the attached drawings, but it should be understood that the embodiments and the attached drawings are only used for the illustrative description of the present invention and do not limit the protection scope of the present invention in any way. All reasonable variations and combinations that fall within the spirit of the invention are intended to be within the scope of the invention.

Referring to fig. 1, according to the technical solution of the present invention, a specific seismic forward modeling method based on least square conjugate gradient iteration includes:

s1, inputting a speed model time harmonic motion equation after Fourier transformation to a frequency domain;

s2, setting corresponding forward modeling parameters according to the input speed model;

s3, according to the set parameters, iterative solution is carried out on the time harmonic motion equation through a least square conjugate gradient algorithm;

s4, judging whether the relative error of the current iteration result meets the preset error range or not, and if not, continuing to perform iteration solution;

and S5, outputting a forward modeling result including a single-frequency data sheet, a synthesized wave field snapshot or an earthquake single shot record after the relative error meets a preset error range.

More specifically, the speed model may be any one of simple models such as a homogeneous medium model and a stratified medium model, or any one of standard models such as a Marmousi model and an Overthrust model.

Some more specific forward modeling parameters may include, for example:

spatial separation dx (unit: m) in the x direction;

the spatial separation dz (unit: m) in the z direction;

a seismic source spatial location;

the dominant frequency MF (unit: Hz) of the seismic source wavelet;

forward modeling duration T (unit: s);

the number of frequency slices;

a frequency interval Δ F;

the number of PML absorption layers and the absorption coefficient;

iterative initial value x of unknown quantity to be solved₀And its preset error epsilon.

In particular implementations, the corresponding parameter settings may be slightly adjusted as needed for forward output, such as single frequency slices, wavefield snapshots, or seismic recordings.

In the above process, a specific iterative solution process of the least square conjugate gradient algorithm includes:

s31 discretizing the time-harmonic motion equation by a difference method to obtain a linear system model as follows:

ap ═ b (1), where a denotes the coefficient matrix, p denotes the pressure field vector, and b denotes the dispersion of the seismic source terms in the equation, specifically, b ═ S (ω) δ (x-x)_s)δ(z-z_s) (2) where S (ω) represents a source term in the frequency domain after Fourier transform, δ () is a Dirac function, x represents an abscissa argument, z represents an ordinate argument, x represents a frequency domain after Fourier transform, and_srepresenting the source abscissa, z_sRepresenting a source ordinate;

preferably, the difference method is a nine-point finite difference method in the prior art, and the time-harmonic motion equation is processed according to the difference format, the optimized difference coefficient and the like to obtain the discretization form.

The structure of the discrete time-harmonic motion equation obtained in the preferred mode is a sparse diagonal band matrix, for example, the number of grid points of the speed model in the x direction is N_xThe number of grid points in the z direction is N_zThen the matrix is one (N)_x*N_z)*(N_x*N_z) A matrix of size in which non-zero elements do not exceed 9 × N at most_x*N_zThe majority of the elements are 0 and 0.

S32 preprocesses the linear system model Ap ═ b, and obtains an equivalent linear system model as follows:

(AM^-1)v＝b,v＝Mp (3)

where M denotes a preprocessing matrix having the same sparse structure as a, v denotes an intermediate vector, and the solution p of the equivalent linear system is M^-1v。

Preferably, the preprocessing selects diagonal preprocessing according to the diagonal dominance property of the coefficient matrix to be solved, and can accelerate the solution of convergence speed.

S33 setting iteration initial value x₀And an iterative preset error epsilon;

wherein the initiation isValue x₀Can be set as a zero vector, and can also be set as other values; the preset error epsilon can be flexibly selected according to the actual situation, and the inventor finds that the preset error epsilon has a better comprehensive effect when the preset error epsilon is set to be 0.01%, namely the iteration frequency, the calculation time and the calculation precision are balanced, if the requirement on the precision of a calculation result is higher, the set relative error can be properly reduced, but the iteration frequency and the calculation time can be increased.

S34 iteratively calculating a solution p-M of the equivalent linear system by a Least-Squares Conjugate Gradient iteration (LSCG) method^-1v；

Wherein, the conjugate gradient parameter of LSCGThe following were used:

where k denotes the number of iterations, θ_kRepresenting a least squares parameter (which varies accordingly with the number of iterations k),the Conjugate Gradient parameters in the CG (Conjugate Gradient) method proposed by hestees and Stiefel in 1952 are shown as follows:

wherein, g_k＝▽f(x^k) Denotes x^kA gradient matrix ofDenotes x^k+1The gradient matrix g_k+1Transpose of (x)^k+1Representing the value of x, y, for the (k + 1) th iteration_kRepresents the parameter g_kDifference of (a) y_k＝g_k+1-g_k,Representing the search direction matrix d_kIs transposed, andi.e. the initial search direction is the initial value x₀Direction of negative gradient of (d), search direction thereafter passing through the parameterSearching iteration step length and determining iteration direction from the current position, and further approximating the solution of a linear equation;

the conjugate gradient parameters proposed by Dai and Yuan in 1999 are expressed as follows:

θ_kthe following least squares problem is solved:

wherein the content of the first and second substances,

s_k＝α_kd_k(10)

wherein the content of the first and second substances,representing the search direction of the LSCG conjugate gradient parameter,representing the search direction, s, of the extended CG conjugate gradient parameter_kRepresenting the step size alpha_kIn the search direction d_kProjection of_kThe iteration step size is indicated.

The calculation solving process may further include:

s341 iterates the initial value x₀Carry-over (4) to obtain the search direction d of the first iteration₀And an iteration step size alpha₀；

S342 along the search direction d of the first iteration₀Using iteration step size alpha₀Performing iterative calculation of the first linear system to obtain the result x of the first iteration₁；

S343 the result x of the first iteration₁Bringing inCalculating to obtain the current iteration error epsilon₁；

S344 judging the current iteration error epsilon₁Whether the error is smaller than a preset iteration error epsilon or not; if the sum of the sums of the sum; if so, repeating S341-S343 to obtain a new iterative search direction d_kAnd an iteration step size alpha_kCalculating a new iteration result x of a new round using a new search direction and a new iteration step size_kCalculating a new iteration error epsilon_kUp to epsilon_kJumping out of iteration when the epsilon is less than epsilon to obtain an approximate solution of a linear system;

and S345, performing subsequent forward modeling on the obtained linear system approximate solution.

According to the above embodiments, the present invention further provides some simulation examples as follows.

The LU direct method used in the following embodiments is a method of solving by decomposing an object matrix into a product of a unit lower triangular matrix and an upper triangular matrix, or a product of them and a permutation matrix.

Example 1

SelectingGrid size N_x＝N_zThe discretized coefficient matrix structure of the uniform medium model of 100 is shown in fig. 2, and it can be seen that the coefficient matrix is highly sparse and has a size of 10000 × 10000, and the specific structure of the 9-point difference can be seen from the partially enlarged part in the figure.

Example 2

Selecting a two-dimensional layered model as shown in fig. 3, wherein the model mesh size is as follows: n is a radical of_x＝N_z300, and 5m, setting parameters including:

forward time is 1s, PML absorption layer number is 10, seismic source coordinates (150 ), seismic source main frequency is 30Hz, preset error is 0.01%, and initial vector x₀Is zero.

According to the forward modeling method of the present invention, after the relative error satisfies the preset error range, the output 30Hz frequency slice is as shown in fig. 5, and the wave field snapshot at the time of 312ms is as shown in fig. 6, so that it can be seen that: FIG. 5 is a clear single-frequency data sheet, in which the wave number of the upper low-speed layer is greater than that of the lower high-speed layer, and the wave number relationship between the upper low-speed layer and the lower high-speed layer is also clearly shown in accordance with the velocity interface of the theoretical relationship of the velocity model, and the simulation result is correct; in the wave field snapshot of fig. 6, the waveform characteristics are clear, the reflected wave transmitted waves can be seen from the upper part and the lower part of the interface respectively, and the simulation result is correct.

Example 3

A Marmousi model (Marmousi-1 model) as shown in fig. 4 is selected, and the model mesh size is: n is a radical of_x＝737，N_z751, the model mesh size is N, since samples of 1/4 are resampled in both the horizontal and vertical directions, respectively_x＝184，N_z187, 12.5m for the spatial interval dx and 10m for dz, and setting parameters including:

the forward time is 2s, the number of PML absorption layers is 50, the coordinates of a seismic source (50, 50) and the main frequency of the seismic source is 15Hz, a receiving point is positioned above the shot point, z is 49, the preset error is 0.01 percent, and the initial vector x₀Is zero.

According to the forward modeling method of the present invention, after the relative error satisfies the preset error range, the output 15Hz frequency slice is as shown in fig. 7, and the wave field snapshot at 625ms is as shown in fig. 8, so that it can be seen that: the high-speed zone at the right side of the single-frequency data sheet model in FIG. 7 has a certain display on the data sheet, and because the extracted main frequency is lower, the characteristics are not clear, and the simulation result is correct; in the wave field snapshot of FIG. 8, the first strong velocity interface near the seismic source is characterized clearly, the projected reflected wave is clear, and the simulation result is correct.

The output seismic single shot record with the PML processing area removed is shown in fig. 9, in which graph (a) is the record obtained by the LSCG iterative method of the present invention, and graph (b) is the record obtained by the LU direct method as a comparison, and it can be seen that the waveform characteristics of both are substantially consistent, demonstrating that the present invention can obtain the same accuracy as the direct method.

Example 4

The model Marmousi-1 as in example 3, i.e. the first model (original model mesh size: N)_x＝737，N_z751, the grid size is obtained by resampling 1/4 samples in both the horizontal and vertical directions, respectively, to N_x＝184，N_zA Marmousi-1 model of 187);

the Marmousi-2 model, the second model, is as follows:

the original model mesh size is: n is a radical of_x＝737，N_z751, the grid size is obtained by resampling 1/2 samples in both the horizontal and vertical directions, respectively, to N_x＝368，N_zMarmousi-2 model 375;

and the mesh size is N_x＝737、N_z751, the original Marmousi model, the third model.

The space intervals are set to be dx-12.5 m and dz-10 m;

the simulation process of each model is characterized in that hardware parameters are as follows: intel (R) Xeon (R) Silver 4110CPU @2.1GHz (8-core 16 thread), which is processed by a single thread under RAM 32 GB; the parameters were set as in example 3, except that the forward time was extended in equal proportion to the model size, and the LU direct method was used for comparison.

According to the forward modeling method of the present invention, after the relative error satisfies the predetermined error range, the memory occupied by the calculation and the calculation time of each model are as follows:

marmousi model 30Hz frequency slice calculation data table

It can be seen that, as the size of the model increases, the memories occupied by the LU method and the LSCG iteration method both increase linearly, but the memory occupied by the LSCG is reduced by 15-20% compared with the LU, and the larger the model is, the more the memory is reduced; the model size is increased by 4 times, the LU method calculation time is increased by nearly 15 times, while the LSCG iteration method calculation time is related to the iteration number, and the increased calculation time is not much due to good convergence, compared with the LU method, the calculation time is reduced by more than 70% in the second and third models, and even reduced by 87.98% in the third model.

The above examples are merely preferred embodiments of the present invention, and the scope of the present invention is not limited to the above examples. All technical schemes belonging to the idea of the invention belong to the protection scope of the invention. It should be noted that modifications and embellishments within the scope of the invention may be made by those skilled in the art without departing from the principle of the invention, and such modifications and embellishments should also be considered as within the scope of the invention.