A frequency domain forward modeling method based on a heterometric 21-point stretched grid

Abstract

The application discloses a frequency domain forward simulation method based on a different distance 21-point stretching type grid, relates to the technical field of petroleum geophysical exploration, and inputs a velocity field, a frequency domain seismic source wavelet and observation system parameters; derives a DD21-point difference format; solves coefficient optimization and performs frequency dispersion analysis; applies PML boundary conditions to construct a DD21-point format PML wave equation; constructs an impedance matrix based on the DD21-point format PML wave equation; solves the impedance matrix to calculate wave field values in a frequency-space domain; transforms the wave field in the frequency-space domain to a time-space domain; extracts a wave field of a detection point; and outputs a frequency domain forward simulation result. The method is more suitable for a different distance 21-point difference format of rectangular grid sampling, when a grid spacing ratio is greater than 1.5, only 2.13 grid points per wavelength are needed to control a phase velocity error within 1%, the operation efficiency and precision of the frequency domain forward simulation are improved, and a more efficient algorithm is provided for subsequent inversion and imaging processing.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of petroleum geophysical exploration technology, and particularly relates to a frequency domain forward simulation method based on a 21-point stretched grid with different distances. BACKGROUND

[0002] At present, frequency domain finite difference forward simulation plays an important role in seismic inversion and imaging. Compared with time domain forward, since each frequency is calculated separately without interference, it is easy to perform parallel calculation, has great advantages for multi-source simulation, and does not need to consider the stability problem caused by cumulative error.

[0003] However, large memory occupation and large calculation amount have always been the problems of frequency domain forward simulation, which leads to the problem that it cannot be effectively applied to the inversion and imaging of actual data, and the operation efficiency of frequency domain forward simulation is low. SUMMARY

[0004] In order to solve the technical problems of large memory occupation and large calculation amount of the frequency domain simulation method, the present application discloses a frequency domain forward simulation method based on a 21-point stretched grid with different distances, which improves the operation efficiency of frequency domain forward simulation.

[0005] To achieve the above purpose, the present application adopts the following technical scheme:

[0006] A frequency domain forward simulation method based on a 21-point stretched grid with different distances, the specific steps are:

[0007] (1) input velocity field, frequency domain source wavelet, and observation system parameters;

[0008] (2) derive DD21-point difference format;

[0009] (3) optimize the coefficient and perform dispersion analysis;

[0010] (4) apply PML boundary conditions and construct DD21-point format PML wave equation;

[0011] (5) construct impedance matrix based on DD21-point format PML wave equation;

[0012] (6) solve the impedance matrix and calculate the wave field value in the frequency-space domain;

[0013] (7) transform the wave field in the frequency-space domain to the time-space domain;

[0014] (8) extract the wave field of the detection point and output the frequency domain forward simulation result.

[0015] Further, in step (2), the wave equation of two-dimensional isotropic medium in the frequency domain is:

[0016]

[0017] where P is the wavefield value, x represents the horizontal position, z represents the vertical position, ω is the angular frequency, v is the velocity, and S(ω) is the frequency domain source;

[0018] When Δx≥Δz, the DD21 point difference format is represented as:

[0019]

[0020] where c i , d i , b i (i=0, 1, 2…, 7) are weighting coefficients with fixed values, and the following relationship exists:

[0021]

[0022] Let R=Δx / Δz, a i =c i +R 2 d i (i=0, 1, 2…, 7), and substitute into equation (2), then:

[0023]

[0024] According to equation (3), a i satisfies the following relationship:

[0025] a0+2a1+2a2+4a3+2a4+2a5+4a6+4a7=0 (5);

[0026] When Δx<Δz, the DD21 point difference format is represented as:

[0027]

[0028] Let R=Δz / Δx, a i =R 2 c i +d i (i=0, 1, 2…, 7), and substitute into equation (6), then:

[0029]

[0030] Further, in step (3), a plane wave is introduced to perform frequency dispersion analysis research by using the classical frequency dispersion analysis method:

[0031] When Δx≥Δz, the plane wave equation is substituted into equation (4), and the phase velocity dispersion relationship is obtained as:

[0032]

[0033] where v ph represents phase velocity, G represents the number of grid points per wavelength, A = cos(2πsinθ / G), B = cos(4πsinθ / G), C = cos(6πsinθ / G), D = cos(2πcosθ / RG), and θ represents the propagation angle;

[0034] It is assumed that the phase velocity residual reaches a minimum value, i.e.:

[0035]

[0036] Equation (8) is brought into equation (9), and the optimization coefficients a i and b i are determined by using the optimization program multistart in MATLAB, the value range of θ is [0, π / 2], the value range of (1 / G) max is [0.35, 0.5], when R = 1, the value range of (1 / G) is [0, 0.35], and when R ≥ 2, the value range of (1 / G) is [0, 0.5];

[0037] When Δx < Δz, the plane wave equation is brought into equation (7), and the phase velocity dispersion relation is obtained as:

[0038]

[0039] where A' = cos(2πcosθ / G), B' = cos(4πcosθ / G), C' = cos(6πcosθ / G), and D' = cos(2πsinθ / RG);

[0040] Due to the symmetry of the method, for the same R, the optimization coefficients take the same value whether Δx < Δz or Δx ≥ Δz.

[0041] Further, in step (4), the acoustic wave equation with PML boundary is:

[0042]

[0043] In the formula, η x and η z respectively represent the attenuation functions in the x and z directions:

[0044]

[0045] In the formula, f represents the main frequency of the source, i.e. describes the vibration frequency of the source, i is the imaginary unit, I x and I zwhere L denotes the length of the point in the left and right and the upper and lower absorbing boundary and the adjacent four model boundaries, respectively pml is the width of the PML boundary, a is a constant to control the degree of attenuation of the boundary condition, and the empirical value is 1.79;

[0046] When Δx≥Δz, the PML wave equation of the DD21 point difference format is:

[0047]

[0048] The acoustic wave equation is converted to the frequency-wavenumber domain:

[0049]

[0050] where, denotes the wavenumber in the horizontal direction, denotes the wavenumber in the vertical direction; at the same time, the plane wave equation is brought into equation (2):

[0051]

[0052] Comparing equation (14) and equation (15), let:

[0053]

[0054] Subsequently, by minimizing the error, that is:

[0055]

[0056] Bring equation (16) into equation (17), and also use the optimization program multistart in MATLAB to determine the optimization coefficient c i , and the value range of θ and 1 / G is consistent with the above determination of a i and b i ; the optimization coefficients c i and d i satisfy a i =c i +R 2 d i (i=0, 1, 2…, 7) and equation (3), so the optimization coefficient d i can be easily solved.

[0057] When Δx<Δz, due to the symmetry of the grid, the optimization coefficients c i and d i take the same value.

[0058] Further, in step (5), considering the PML boundary and the source term, taking Δx≥Δz as an example, the construction process of the impedance matrix is derived:

[0059] After adding the source term, formula (13) is arranged as:

[0060]

[0061] Wherein:

[0062]

[0063] When Δx≥Δz, a row-based reading mode is adopted, that is, a two-dimensional matrix is read into a column vector row by row from left to right; when Δx<Δz, a column-based reading mode is adopted, that is, a two-dimensional matrix is read into a column vector column by column from top to bottom; the linear equation set is expressed in a matrix form:

[0064] AP=-S (20);

[0065] Wherein, A is an impedance matrix with a size of (nx×nz)×(nx×nz), also known as a large sparse matrix, P and S are column vectors of (nx×nz).

[0066] Further, in step (6), the wave field value can be obtained by solving the matrix equation (20); using LU decomposition, the matrix A is decomposed into the product of a lower triangular matrix L and an upper triangular matrix U, and the solution of the matrix equation (20) is:

[0067] P=-U -1 L -1 S (21)。

[0068] The method has the advantages that, compared with the prior art, the method makes two improvements on the difference operator, first, the precision of the difference operator is improved by using a general optimization method of finite difference solution; second, the calculation amount of LU decomposition of the impedance matrix is reduced by limiting the sub-diagonal line of the impedance matrix to the two sides of the main diagonal line using a stretched grid.

[0069] The method alleviates the problems of large memory occupation and large calculation amount in frequency domain forward simulation, improves the operation efficiency and precision of frequency domain forward simulation, and provides a more efficient algorithm for subsequent inversion and imaging processing. BRIEF DESCRIPTION OF DRAWINGS

[0070] Figure 1 The figure is a flowchart of the application;

[0071] Figure 2 The figure is a schematic diagram of DD21 point difference format, the left figure is the case of Δx≥Δz, and the right figure is the case of Δx<Δz;

[0072] Figure 3 The figure is the dispersion curve of DD21 point format and ADM25 point format under different sampling interval ratios, the left figure is the DD21 point format, and the right figure is the ADM25 point format;

[0073] Figure 4 Wavefield snapshots at 800 ms for the homogeneous model calculated using different difference schemes, in order: ADM 9-point difference scheme, ADM 25-point difference scheme and DD 21-point difference scheme;

[0074] Figure 5 Comparison of single-channel records of three numerical methods and the analytical method for the homogeneous model, in order: ADM 9-point scheme, ADM 25-point scheme and DD 21-point scheme;

[0075] Figure 6 The tested double-layer model;

[0076] Figure 7 Wavefield snapshots at 750 ms for the double-layer model calculated using different difference schemes, in order: ADM 9-point difference scheme, ADM 25-point difference scheme and DD 21-point difference scheme;

[0077] Figure 8 Forward records of the double-layer model calculated using different difference schemes, in order: ADM 9-point difference scheme, ADM 25-point difference scheme and DD 21-point difference scheme;

[0078] Figure 9 The intercepted partial Marmousi model;

[0079] Figure 10 Wavefield snapshots at 550 ms for the partial Marmousi model calculated using different difference schemes, in order: ADM 9-point difference scheme, ADM 25-point difference scheme and DD 21-point difference scheme;

[0080] Figure 11 Forward records of the partial Marmousi model calculated using different difference schemes, in order: ADM 9-point difference scheme, ADM 25-point difference scheme and DD 21-point difference scheme;

[0081] Figure 12 A local enlarged view of Figure 11 DETAILED DESCRIPTION

[0082] In order to make the objects, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are some but not all of the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative work fall within the scope of protection of the present application.

[0083] ​The application provides a frequency domain forward simulation method based on a different distance 21-point stretching type grid, and the method comprises the following steps:

[0084] (1) input a velocity field, a frequency domain seismic source wavelet and observation system parameters.

[0085] (2) derive a DD21-point difference format.

[0086] The wave equation of a two-dimensional isotropic medium in the frequency domain is as follows:

[0087]

[0088] In the formula, P represents a wave field value, x represents a horizontal position, z represents a vertical position, omega represents an angular frequency, v represents a velocity, and S (omega) represents a frequency domain seismic source.

[0089] It is assumed that the approximation of the wave field value of adjacent points to the spatial derivative is based on the distance from the adjacent points to the center point, and if the distance from the adjacent points to the center point is equal, the contribution weight values of the adjacent points are the same. The acceleration term omega 2 P / v 2 is also expressed as a weighted average of all points based on the distance, so that the discretization processing of the wave equation is realized.

[0090] When Delta x is greater than or equal to Delta z, the DD21 difference format is expressed as:

[0091]

[0092] In the formula, c i , d i , b i (i=0, 1, 2,..., 7) are all weighted coefficients with fixed values, and the following relationship exists:

[0093]

[0094] Let R=Delta x / Delta z, a i = c i +R 2 d i (i=0, 1, 2,..., 7), and the formula (2) is brought into, so that:

[0095]

[0096] According to the formula (3), a i satisfies the following relationship:

[0097] a0+2a1+2a2+4a3+2a4+2a5+4a6+4a7=0 (5)

[0098] When Δx< Δz, the DD21 difference format is expressed as:

[0099]

[0100] Let R = Δz / Δx, a i = R 2 c i + d i (i = 0, 1, 2…, 7), which is brought into equation (6), and then we have:

[0101]

[0102] (3) Optimize the coefficients and analyze the dispersion.

[0103] Using the classical dispersion analysis method, a plane wave is introduced to analyze the dispersion.

[0104] When Δx≥ Δz, the plane wave equation is brought into equation (4), and the phase velocity dispersion relation is obtained as:

[0105]

[0106] where: v ph represents the phase velocity, G represents the number of grid points per wavelength, A = cos(2πsinθ / G), B = cos(4πsinθ / G), C = cos(6πsinθ / G), D = cos(2πcosθ / RG), θ represents the propagation angle. Next, assume that the phase velocity residual reaches a minimum value, i.e.

[0107]

[0108] Bring equation (8) into equation (9), and use the optimization program multistart in MATLAB to determine the optimized coefficients a i and b i . The value range of θ is [0, π / 2]. Since the dispersion degree of different sampling interval ratios R is different, a reasonable value range of 1 / G needs to be found. Through testing, (1 / G) max increases with the increase of R, but reaches a stable value when R = 2. Finally, the range of (1 / G) max is taken as [0.35, 0.5]. When R = 1, the value range of 1 / G is [0, 0.35], and when R ≥ 2, the value range of 1 / G is [0, 0.5].

[0109] When Δx< Δz, the DD21 difference format is expressed as:

[0110]

[0111] where A' = cos(2πcosθ / G), B' = cos(4πcosθ / G), C' = cos(6πcosθ / G), D' = cos(2πsinθ / RG). Due to the symmetry of the method, the optimized coefficients take the same values for the same R, regardless of whether Δx < Δz or Δx ≥ Δz.

[0112] (4) Apply PML boundary conditions to construct the DD21-point PML wave equation.

[0113] The PML boundary condition was first introduced by Berenger in electromagnetics. The acoustic wave equation with PML boundary is:

[0114] where η x and η z are the attenuation functions in x and z directions, respectively:

[0115] η x = 1 - 2πaf(i / ω)(I x / L pml )

[0116] η z = 1 - 2πaf(i / ω)(I z / L pml ) (12)

[0117] where f represents the main frequency of the source, i.e., the vibration frequency of the source, I x and I z represent the lengths of the points in the left and right and upper and lower absorbing boundaries, respectively, and L pml is the width of the PML boundary, and a is a constant to control the degree of attenuation of the boundary condition, which is 1.79 according to previous experience.

[0118] Taking Δx ≥ Δz as an example, the DD21-point difference format PML wave equation is:

[0119]

[0120] To create a PML boundary, additional constraints are needed to solve c i and d i . Convert the acoustic wave equation to the frequency-wavenumber domain:

[0121]

[0122] where k is the horizontal wavenumber, and k kzdenotes the wave number in the vertical direction. At the same time, the plane wave equation is brought into equation (2):

[0123]

[0124] Comparing equation (14) with equation (15), let:

[0125]

[0126] Since the optimization coefficients c i and d i satisfy a i =c i +R 2 d i (i=0, 1, 2…, 7) and equation (3), the 16 coefficients can be simplified to 7 independent coefficients. By minimizing the error:

[0127]

[0128] Similarly, the optimization coefficients c i , θ and 1 / G are determined by using the optimization program multistart in MATLAB, and the value ranges of the optimization coefficients are consistent with those of a i and b i . For the case of Δx<Δz, due to the symmetry of the grid, the optimization coefficients c i and d i take the same value.

[0129] (5) Construct the impedance matrix based on the DD21 point format PML wave equation.

[0130] Considering the PML boundary and the source term, taking Δx≥Δz as an example, the construction process of the impedance matrix is derived. After adding the source term, equation (13) is rearranged as:

[0131]

[0132] wherein:

[0133]

[0134] Equation (18) is a linear equation constructed with P m,n as the center point. For each point in the model, a corresponding linear equation can be constructed with the center point, and thus nx×nz linear equations are obtained. By solving the linear equation system composed of the nx×nz linear equations, the wave field values of each grid point in the model can be calculated. In this process, the wave field needs to be rearranged from an nx×nz two-dimensional matrix to an nx×nz column vector.

[0135] It is worth noting that when Δx≥Δz, the reading mode is mainly in the form of row, that is, the two-dimensional matrix is read into a column vector row by row from left to right; and when Δx<Δz, the reading mode is mainly in the form of column, that is, the two-dimensional matrix is read into a column vector column by column from top to bottom. The linear equation set is expressed in the form of a matrix:

[0136] AP=-S (20)

[0137] Wherein, A is an impedance matrix with a size of (nx×nz)×(nx×nz), also known as a large sparse matrix, and P and S are column vectors of (nx×nz).

[0138] (6) Solve the impedance matrix to calculate the wave field value in the frequency-space domain.

[0139] Solving the matrix equation (20) can obtain the wave field value, and usually using LU decomposition, the matrix A is decomposed into the product of a lower triangular matrix L and an upper triangular matrix U, and the solution of the matrix equation (20) is:

[0140] P=-U -1 L -1 S(21)

[0141] The LU decomposition is the most time-consuming part in the frequency domain forward simulation. And formula (19) shows that the matrix A is independent of the source, so for multi-shot simulation, only one LU decomposition of the matrix A is needed, which greatly improves the calculation efficiency of multi-shot simulation. In addition, the bandwidth of the non-zero elements in the matrix A determines the amount of calculation of the LU decomposition, and the larger the bandwidth, the greater the amount of calculation. Under the condition that the number of sampling points is the same, the bandwidth of the DD21 point and the average derivative 9 point difference format is not much different, which is about half of the bandwidth of the average derivative 25 point difference format, thereby reducing the calculation cost and improving the calculation efficiency.

[0142] (7) Transform the wave field in the frequency-space domain to the time-space domain.

[0143] (8) Extract the wave field of the detection point, and output the frequency domain forward simulation result.

[0144] Application example

[0145] The frequency domain forward simulation method based on the DD21 point stretched grid of the application can be applied to uniform model, double-layer model and part of Marmousi model, and ideal calculation effect is obtained. Figure 1 and Figure 2 The flowchart of the frequency domain forward simulation method based on the DD21 point stretched grid is shown, and the DD21 point difference grid schematic diagram is shown.

[0146] Figure 3The phase velocity dispersion curves of the DD21-point and ADM25-point differential schemes for Δx ≥ Δz under different sampling interval ratios are shown. Under the constraint of a phase velocity error of less than 1%, the number of grid points required per wavelength for the DD21-point scheme varies with the grid spacing ratio. R The increase and decrease, in R It reaches a stable value around 2. R When R = 1, 3.17 grid points are needed per wavelength; when R = 1.5, 2.13 grid points are needed per wavelength. R≥2 At that time, 2.04 grid points are needed per wavelength. However, the number of grid points required per wavelength for ADM25 is essentially independent of the grid spacing ratio. R The frequency response varies, requiring 2.78 grid points per wavelength. Compared to the ADM25 grid, the DD21 grid format provides better dispersion suppression when R ≥ 1.5.

[0147] The uniform model has a size of 3600×3600m and a velocity of 2000m / s. The source is a Ricker wavelet located at the center of the model, with a dominant frequency of 20Hz. The time sampling interval is 1ms, and the sampling time is 2000ms. The grid spacing is set to Δx = 15m and Δz = 10m, and the grid size is 241×361. Figure 4 800ms wavefield snapshots calculated using the ADM9-point difference scheme, ADM25-point difference scheme, and DD21-point difference scheme are presented respectively. The ADM9-point result exhibits the most severe numerical dispersion, the ADM25-point result also shows significant dispersion, while the DD21-point result is almost unaffected by numerical dispersion. To further verify the accuracy of the method, a receiver point was set 450m to the right of the source. The analytical and numerical solutions at this receiver point were compared. Figure 5 The comparison between numerical and analytical solutions calculated using three difference schemes is shown, with the DD21-point calculation results showing almost perfect agreement with the analytical solution. This more intuitively illustrates the accuracy of the method.

[0148] The two-layer model measures 3000×3000m. The upper 1000m layer has a speed of 2000m / s, and the lower 2000m layer has a speed of 3000m / s. Figure 6 As shown. The Reichsson source wavelet is located at x = 1500 m and z = 100 m, with a dominant frequency of 20 Hz. The time sampling interval is 1 ms, and the sampling time is 1500 ms. The grid spacing is set to Δx = 15 m and Δz = 10 m, and the grid size is 201 × 301. Figure 7 A 750ms wavefield snapshot obtained through three difference schemes is shown. Figure 8The shot records calculated by the three difference formats are shown. As can be seen from the figures, the ADM9 point and the ADM25 point all have different degrees of dispersion phenomenon; and the DD21 point difference format shows good adaptability to the layered model and has higher accuracy than the other two difference formats.

[0149] The intercepted part of the Marmousi model has a size of 1845*1065m, the grid spacing is set to be Δx=7.5m and Δz=5m, the grid size is 247*214, as shown in Figure 9 To improve the resolution of the seismic record, the Ricker source wavelet main frequency is taken as 30Hz, the position is set at x=922.5m and z=50m. The time sampling interval is 1ms, and the sampling time is 2000ms. The 550ms wave field snapshot calculated by the three difference formats ( Figure 10 ) can clearly show that the DD21 point has higher accuracy than the ADM9 point and the ADM25 point. Figure 11 and Figure 12 The shot records calculated by the three difference formats and the local enlarged view of the shot record are also shown, and it can be seen that the DD21 point difference format also has good applicability and forward accuracy to the complex model. In addition, compared with the ADM25 point, the DD21 point uses 21-point grid difference, reduces the storage requirement and the calculation amount of matrix decomposition, and further realizes the high-efficiency and high-precision frequency domain forward simulation.

[0150] Of course, the above description is not a limitation on the present application, and the present application is not limited to the above examples. Changes, modifications, additions or substitutions made by those skilled in the art within the essential scope of the present application should also be within the protection scope of the present application.

Claims

1. A frequency domain forward modeling method based on a heterogeneous 21-point stretched grid, characterized in that, The specific steps are as follows: (1) input velocity field, frequency domain source wavelet, and observation system parameters; (2) derive DD21 point difference format; (3) optimize the coefficient and perform dispersion analysis; (4) apply PML boundary condition to construct DD21 point format PML wave equation; (5) construct impedance matrix based on DD21 point format PML wave equation; (6) solve the impedance matrix to calculate the wave field value in the frequency-space domain; (7) transform the wave field in the frequency-space domain to the time-space domain; (8) extract the wave field of the detection point and output the frequency domain forward simulation result; In step (2), the wave equation of two-dimensional isotropic medium in the frequency domain is: （1）； In the formula, P is the wave field value, x represents the horizontal position, z represents the vertical position, ω is the angular frequency, v is the velocity, and S(ω) is the frequency domain source; When the DD21 point-difference format is represented as: （2）； wherein , , are weighting factors that are constant values, and the following relationship exists: （3）； Let , , into equation (2), then we have: （4）； According to formula (3), satisfies the following relationship: （5）； When the DD21 point-difference format is represented as: （6）； Let , , into equation (6), then we have: （7）。 2. The method of claim 1, wherein, In step (3), a classical dispersion analysis method is used to introduce a plane wave for dispersion analysis research: When The plane wave equation is brought into equation (4), and the phase velocity dispersion relation is obtained as follows: （8）； wherein, denotes phase velocity, denotes the number of grid points per wavelength, , , , , denotes the propagation angle; It is assumed that the phase velocity residual reaches the minimum value, that is: （9）； Substituting equation (8) into equation (9), the optimization coefficients are determined using the multistart optimization program in MATLAB. and , The range of values ​​is , The range is taken as ,when hour, The range of values ​​is ,when hour, The range of values ​​for is ; When The plane wave equation is brought into equation (7), and the phase velocity dispersion relation is obtained as follows: （10）； wherein , , , ; Due to the symmetry of the method, the optimization coefficients take the same value for the same or , for the same .

3. The method of claim 2, wherein, In step (4), the acoustic wave equation with PML boundary is: （11）； wherein with respectively represent and attenuation functions in the direction （12）； where, denotes the source dominant frequency, i.e., describes the vibration frequency of the source, is the imaginary unit, and denotes the length of the point in the left and right and up and down absorbing boundary and the adjacent four model boundaries, respectively, is the width of the PML boundary, is a constant to control the degree of attenuation of the boundary condition, and the empirical value is 1.79; The DD21 point-difference format PML wave equation is: （13）； Convert the acoustic wave equation to the frequency-wavenumber domain: （14）； wherein, denotes the wave number in the horizontal direction, denotes the wave number in the vertical direction; and the plane wave equation is brought into equation (2): （15）； By comparing formula (14) and formula (15), let: （16）； Subsequently, by minimizing the error, that is: （17）； Substituting equation (16) into equation (17), the optimization coefficients are determined using the MATLAB optimization program multistart. ,and and The range of values ​​is the same as that obtained above. and Maintain consistency; optimize coefficients and satisfy By combining equation (3), the optimization coefficients can be obtained. ; When the optimization coefficients and take the same value due to the symmetry of the grid.

4. The method of claim 3, wherein, In step (5), the PML boundary and the source term are taken into account to derive the construction of the impedance matrix: when After adding the source term, formula (13) is arranged as: （18）； Wherein: （19）； In case, the reading-in mode is mainly by row, i.e. reading the two-dimensional matrix into column vector row by row from left to right; while in case, the reading-in mode is mainly by column, i.e. reading the two-dimensional matrix into column vector column by column from top to bottom. Express the linear equation set in the form of a matrix: （20）； wherein, is an impedance matrix of size also known as a large sparse matrix, and is a column vector of size .

5. The method of claim 4, wherein, In step (6), the wavefield values are obtained by solving the matrix equation (20) using Decomposing, the matrix is decomposed into a product of a lower triangular matrix and an upper triangular matrix The solution of the matrix equation (20) is then given by （21）。

Images