Seismic wave field numerical simulation method based on wavelet domain Pre-GSOR iteration

Abstract

The invention relates to a seismic wave field numerical simulation method based on wavelet domain Pre-GSOR iteration, and the method comprises the steps: introducing an attenuation factor into a background wave number, rewriting a Helmholz equation, obtaining a complex wave number L-S equation according to a representation theorem, and rewriting the complex wave number L-S equation into a normalized equivalent L-S equation; a discrete equation set corresponding to the normalized equivalent complex wave number L-S equation discretized through the # imgabs0 # method is converted into a wavelet domain, and after threshold processing, a new sparse linear equation set is generated; a wavelet domain sparse linear equation set corresponding to a normalized equivalent complex wave number L-S equation is solved by using a precondition generalized successive hyper-relaxation iteration method, a non-diagonal sparse precondition matrix is generated based on an ARMS decomposition method of block ILU decomposition, and the convergence of iteration series under a strong scattering medium is effectively improved. According to the method, the problem that the convergence speed is slow or stagnant when the generalized successive super-relaxation iteration method is used for solving the discrete linear equation set of the L-S equation under the strong scattering medium is solved is solved.

Description

Technical Field

[0001] The invention belongs to the technical field of numerical simulation of seismic wave fields, and specifically relates to a method for numerical simulation of seismic wave fields by using an equivalent complex wave number Lippmann–Schwinger (LS) integral equation, and in particular to a fast numerical iteration method for solving a sparse discrete equation group of an equivalent complex wave number LS equation in a wavelet domain by using a preconditioned generalized successive over-relaxation (Pre-GSOR) iteration method. Background Art

[0002] With the continuous innovation of oil and gas exploration technology, the global oil and gas exploration and development field has been continuously extended and expanded from shallow to deep and ultra-deep layers, from conventional to unconventional, and from shallow water to deep water and ultra-deep water. This type of resource is buried deeper, has stronger heterogeneity and concealment, and the exploration targets are more complex. The observation data received on the ground not only contains a variety of seismic wave field types, but also has developed scattered waves and weak primary reflections. At this time, the migration imaging and inversion methods based on the theory of primary reflection waves are often difficult to achieve accurate imaging and inversion of such complex structures. Only by using full-wavefield numerical simulation and full-wavefield imaging can the seismic wave propagation laws and imaging problems in complex structural areas be effectively characterized.

[0003] When dealing with seismic scattering problems, the perturbation theory can be used to transform the Helmholtz equation describing the scattering problem into the Lippmann–Schwinger (LS) integral equation. The LS equation can be expressed as The LS equation is solved by numerical methods such as the Galerkin method, the moment method and the iterative method. The Born scattering series can be obtained by solving the LS equation using the iterative method, but the Born series converges only when the scatterer is small or the scattering is weak. For strong scattering media, the Born series either converges very slowly or diverges, and the convergence of the Born series deteriorates with increasing frequency. In the prior art, research on improving the convergence properties of the Born scattering series is disclosed. However, when using the LS equation to solve the numerical simulation problem of the seismic wave field, there is no discrete equation group for the LS equation in the prior art, and the mathematical structure and properties of the discrete coefficient matrix and its influence on the convergence properties of the iterative series are analyzed from a numerical perspective. In addition, when the iterative method is used to solve the LS equation in a medium with strong velocity disturbance, the diagonal preconditioning operator in the prior art often has limited acceleration on the iterative method, while the generation of the non-diagonal preconditioning operator requires a large amount of calculation. Therefore, in the case of strong seismic scattering media, the improvement of the convergence of high-frequency numerical iterative methods is still a difficult point in current research. Summary of the invention

[0004] The purpose of the present invention is to provide a fast numerical calculation method for numerical simulation of seismic wave fields based on a wavelet domain preconditioned generalized successive super-relaxation iterative method to solve the equivalent complex wave number LS equation, so as to solve the problem of iterative series divergence and convergence speed stagnation when the numerical iterative method is used to solve the LS equation in a strong scattering medium.

[0005] The objective of the present invention is achieved through the following technical solutions:

[0006] A method for numerical simulation of seismic wave fields based on a wavelet domain preconditioned generalized super-relaxation iteration method comprises the following steps: firstly, a small imaginary component is introduced into the background wave number, and then the LS equation is rewritten to obtain an equivalent complex wave number LS equation; the LS equation is rewritten so that the coefficient matrix after numerical discretization has a band structure, and the wavelet coefficient matrix after threshold processing in the wavelet domain has sparsity; then, an algebraic recursive multi-layer decomposition method based on block ILU decomposition is used in the wavelet domain to generate a sparse approximate inverse as a preconditioned operator to accelerate the generalized successive super-relaxation iteration method; finally, the method is used for numerical simulation calculation of single-frequency seismic wave fields in strong scattering media.

[0007] The present invention is based on the wavelet domain Pre-GSOR iterative seismic wave field numerical simulation method, which specifically includes the following steps:

[0008] a. Read the velocity model, the model size is Nx*Nz, set the source wavelet, source position, detection point position and background velocity;

[0009] b. Uniform division of the integral region Ω V j is the jth partitioning element, N is the number of partitioning elements; take the partitioning constant function space L2(Ω) is the Hilbert space; basis function e j (r) Satisfy:

[0010]

[0011] Where r represents the spatial position on the velocity model;

[0012] c. Based on the interpolation basis function e in step b j (r), construct the interpolation operator:

[0013]

[0014] Among them, r j Indicates V j The center of , v(r) is an arbitrary continuous function.

[0015] d. Construct the equivalent complex wave number LS equation and use the interpolation operator P NThe equivalent complex wave number LS equation adopts weak singular kernel The matrix form of the linear equations is obtained by numerical discretization:

[0016]

[0017] in, is the wave field weighted by velocity disturbance at the model space point in the frequency domain, and u(r i ) are space points r i Velocity disturbance values ​​and wave field values ​​on ; is the incident wave field, s(w) represents the source wavelet, r s represents the earthquake source location, For point r i The value of the complex background Green's function on ; The weight function The complex background wave number Green operator is composed of is the complex scattering potential The diagonal operator formed by

[0018] e. Select the Daubechie wavelet function, use the Mallat pyramid decomposition algorithm, and construct the orthogonal wavelet transform matrix W through the Daubechie wavelet filter to satisfy W T W=WW T =I, apply the orthogonal wavelet transform matrix W to the linear equations formula (12) in step d, and obtain the wavelet domain representation of the equivalent complex wave number LS equation discrete equations:

[0019]

[0020] make Then we can get the linear equations in wavelet domain

[0021]

[0022] According to the compact support and high-order vanishing moment characteristics of the wavelet function, the transformed coefficient matrix Most of the wavelet coefficients are zero or very small; set the threshold η, ignore the values ​​less than η without affecting the calculation accuracy, and the coefficient matrix after threshold processing It has a sparse structure. When the coefficient matrix is ​​sparse enough, the fast matrix-vector multiplication allows the filling of matrix elements to be kept at a low level when solving the approximate inverse of the sparse matrix or generating the inverse operator by factorization.

[0023] f. Solve using the generalized successive super-relaxation iterative algorithm with precondition constraints The iteration format is:

[0024]

[0025] or

[0026]

[0027] The preconditioner γ uses an algebraic recursive multi-layer decomposition method based on block ILU decomposition to decompose the coefficient matrix Decompose to generate; α n is the relaxation factor, which is calculated by minimizing the residual ||r n || to get, take The nth iteration takes α n Make ||r n ||Minimize in Hilbert space:

[0028]

[0029] g. The iterative solution obtained in step f Perform inverse wavelet transform,

[0030]

[0031] Then the desired single-frequency seismic wave field can be obtained;

[0032] h. Calculate the wave fields at all detection points to obtain the seismic wave fields at all detection points.

[0033] Furthermore, step d specifically comprises the following steps:

[0034] d1. Define the velocity disturbance function v0 is the constant background velocity, v(r) is the medium velocity, and according to the wave field superposition principle, the Helmholz equation satisfied by the scalar wave equation in the frequency domain is rewritten as an equivalent source form:

[0035]

[0036] Where ω is the angular frequency, s(ω)δ(rr s ) represents the point source excited at rs, and u(r) is the wave field to be determined. Adding iεu(r) (ε>0) to both ends of formula (3) yields:

[0037]

[0038] Pick where |O| max =max r |O(r)|;

[0039] d2. Define complex background wave number Then there is the complex wave number background Green's function satisfy:

[0040]

[0041] in, r' is the source point of Green's function. The analytical solution is:

[0042]

[0043] In numerical calculations, according to the disturbance intensity of the model, it is necessary to select an appropriate constant a to ensure Greater than 1.0e-6 within the effective calculation area;

[0044] d3, corresponding to the complex background wave number Complex scattering potential Based on Green's theorem, the solution of formula (4) can be expressed in integral form to obtain the complex wave number LS equation:

[0045]

[0046] in, At this time, the integral kernel in the complex wave number LS equation is In addition, it also contains the complex velocity disturbance term;

[0047] d4. Define the wave field with velocity perturbation as weight The equivalent complex wave number LS equation is obtained:

[0048]

[0049] The above formula is an integral equation about u'(r), By normalizing formula (8) point by point as the normalization function, the normalized equivalent LS equation can be obtained:

[0050]

[0051] d5. Define the integral operator Using the interpolation operator P N , replace the continuous part of the integrand with its interpolation function, then the integral operator K Approximate K N It can be expressed as:

[0052]

[0053] Among them, the weight function Then the normalized equivalent complex wave number of the LS equation is The approximate equations can be expressed as:

[0054]

[0055] The matrix (operator) of the above formula is expressed as:

[0056]

[0057] Where I is the identity matrix, The weight function The complex wave number Green operator is composed of is the complex scattering potential The diagonal operator is constructed; by selecting an appropriate attenuation factor ε, the complex velocity perturbation The imaginary component in is not zero, even at the point where the velocity disturbance O(r) is zero, the complex velocity disturbance The amplitude of is never zero.

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

[0059] The present invention proposes a method for solving the equivalent complex wavenumber LS integral equation for numerical simulation of seismic wave field based on the wavelet domain preconditioned generalized successive super-relaxation iteration method, which introduces the complex wavenumber Green's function with imaginary component into the LS integral equation and constructs the normalized equivalent complex wavenumber LS equation. After numerical discretization, the coefficient matrix has a sparse pattern in the wavelet domain. The ARMS decomposition method is further used to generate sparse preconditions to accelerate the convergence speed of the generalized super-relaxation iterative method and improve the convergence of the iterative method. Specifically, it has the following advantages:

[0060] 1. The present invention rewrites the LS integral equation into a normalized equivalent LS equation so that After numerical discretization, the off-diagonal elements of the coefficient matrix are generated only by the complex background Green's function, which can maintain the properties of the Toeplitz matrix. There is no need to store the entire matrix, which can reduce storage.

[0061] 2. The non-diagonal elements of the coefficient matrix in the discrete linear equations of the normalized equivalent LS equation have the Toeplitz property, so that after wavelet transform, the coefficient matrix has a sparse pattern after threshold processing in the wavelet domain, which can effectively reduce the computational complexity of iterative calculation and preconditioning operator generation;

[0062] 3. When solving the wavelet domain linear equation system of the normalized equivalent LS equation, the present invention adopts an algebraic recursive multi-layer decomposition method based on block ILU decomposition to generate a sparse preconditioning operator, which can effectively reduce the computational cost of generating the preconditioning operator and effectively improve the convergence of the generalized successive super-relaxation iterative method under strong scattering media. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for use in the embodiments are briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other related drawings can be obtained based on these drawings without creative work.

[0064] Figure 1 A flow chart of the steps of the seismic wave field numerical simulation method based on the wavelet domain preconditioned generalized super-relaxation iteration method of the present invention;

[0065] Figure 2 Velocity model;

[0066] Figure 3 Sparse pattern of wavelet domain coefficient matrix;

[0067] Figure 4 Normalized residual error curve with the number of iterations;

[0068] Figure 5 Single-frequency source wave field at frequency f = 30 Hz, Figure 5 a and 5b are the real and imaginary parts of the source wave field obtained by the direct method corresponding to the real wave number LS equation, Figure 5 c and 5d are the real and imaginary parts of the source wavefield obtained by the Pre-GSOR iteration method of the normalized equivalent LS equation;

[0069] Figure 6 The scattered wave field of the detection point with frequency f = 30Hz, Figure 6 a and Figure 6 b is the real and imaginary parts of the scattered wave field at the detection point at a frequency of 30 Hz. DETAILED DESCRIPTION

[0070] The present invention will be further described below in conjunction with embodiments:

[0071] The present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments. It is to be understood that the specific embodiments described herein are only used to explain the present invention, rather than to limit the present invention. It should also be noted that, for ease of description, only parts related to the present invention, rather than all structures, are shown in the accompanying drawings.

[0072] It should be noted that similar reference numerals and letters represent similar items in the following drawings, so once an item is defined in one drawing, it does not need to be further defined and explained in the subsequent drawings. At the same time, in the description of the present invention, the terms "first", "second", etc. are only used to distinguish the description and cannot be understood as indicating or implying relative importance.

[0073] The present invention proposes a seismic wave field forward modeling method based on the wavelet domain preconditioned generalized successive over-relaxation iterative method. First, the attenuation factor is introduced into the background wave number and the Helmholz equation is rewritten to obtain the complex wave number LS equation according to the representation theorem. In order to make the coefficient matrix of the complex wave number LS equation have a better mathematical structure after numerical discretization, the complex wave number LS equation is further rewritten as a normalized equivalent LS equation. Then, The discrete equations corresponding to the normalized equivalent complex wave number LS equation discretized by the method are transformed into the wavelet domain, and after threshold processing, a new sparse linear equation system is generated. Finally, the preconditioned generalized successive super-relaxation iterative method is used to solve the wavelet domain sparse linear equation system corresponding to the normalized equivalent complex wave number LS equation. The ARMS decomposition method based on block ILU decomposition generates a sparse preconditioned matrix, which can effectively improve the convergence of the iterative series in strong scattering media.

[0074] The present invention is based on the wavelet domain Pre-GSOR iterative seismic wave field numerical simulation method, comprising the following steps:

[0075] a. Read parameters and velocity model (model size Nx*Nz), set source wavelet, source, detection point parameters and background velocity.

[0076] b. Uniform division of integral region (velocity model) (V j is the jth partitioning element, N is the number of partitioning elements); take the partitioning constant function space (L2(Ω) is Hilbert space); basis function e j (r) Satisfy:

[0077]

[0078] Wherein, r represents the spatial position (spatial point) on the velocity model.

[0079] c. Based on the interpolation basis function e in step b j (r), construct an orthogonal constant interpolation operator:

[0080]

[0081] Among them, r j Indicates V j The center of , v(r) is an arbitrary continuous function.

[0082] d. Construct the normalized equivalent complex wave number LS equation and use the interpolation operator P N The equivalent complex wave number LS equation adopts weak singular kernel The matrix form of the linear equations is obtained by numerical discretization:

[0083]

[0084] in, is the wave field weighted by velocity disturbance at the model space point in the frequency domain, and u(r i ) are space points r i The velocity disturbance value and wave field value on it. is the incident wave field, s(w) represents the source wavelet, r s represents the earthquake source location, For point r i The value of the complex background Green's function on . The weight function The complex background wave number Green operator is composed of is the complex scattering potential The diagonal operator.

[0085] Specifically, define the velocity perturbation function (v0 is the constant background velocity, v(r) is the medium velocity). According to the wave field superposition principle, the Helmholz equation satisfied by the scalar wave equation in the frequency domain is rewritten in the form of an equivalent source:

[0086]

[0087] Where ω is the angular frequency, s(ω)δ(rr s ) indicates excitation at r s is the point source at , and u(r) is the wave field to be determined.

[0088] Adding iεu(r)(ε>0) to both ends of formula (3) yields:

[0089]

[0090] Pick Where |O max =max r |O(r)|.

[0091] Define complex background wavenumber The complex wave number background Green's function satisfy:

[0092]

[0093] in, r' is the source point of Green's function. The analytical solution is:

[0094]

[0095] In numerical calculations, according to the disturbance intensity of the model, it is necessary to select an appropriate constant a to ensure Greater than 1.0e-6 within the effective calculation area.

[0096] Corresponding to the complex background wave number Complex scattering potential Based on Green's theorem, the solution of formula (4) can be expressed in integral form to obtain the complex wave number LS equation:

[0097]

[0098] in, At this time, the integral kernel in the complex wave number LS equation above is divided by In addition, it also contains the complex velocity disturbance term.

[0099] Define the wave field weighted by velocity perturbations The equivalent complex wave number LS equation is obtained:

[0100]

[0101] The above equation is an integral equation about u'(r). By normalizing formula (8) point by point as the normalization function, the normalized equivalent LS equation can be obtained:

[0102]

[0103] Defining the integral operator Using the interpolation operator P N , replace the continuous part of the integrand with its interpolation function, then the integral operator K Approximate K N It can be expressed as:

[0104]

[0105] Among them, the weight function Then the normalized equivalent complex wave number of the LS equation is The approximate equations can be expressed as:

[0106]

[0107] The matrix (operator) of the above formula is expressed as:

[0108]

[0109] Where I is the identity matrix, The weight function The complex wave number Green operator is composed of is the complex scattering potential By selecting an appropriate attenuation factor ε, the complex velocity perturbation The imaginary component in is not zero, even at the point where the velocity disturbance O(r) is zero, the complex velocity disturbance The amplitude of is never zero, which can avoid the singularity problem during normalization.

[0110] e. Select the Daubechie wavelet function, use the Mallat pyramid decomposition algorithm, and use the Daubechie wavelet filter to perform wavelet transform on the equations to construct an orthogonal wavelet transform matrix W that satisfies W. T W=WW T =I. Apply the orthogonal wavelet transform matrix W to the linear equation system (12) in step d to obtain the wavelet domain representation of the equivalent complex wave number LS equation discrete equation system:

[0111]

[0112] make Then we can get the linear equations in wavelet domain

[0113]

[0114] According to the compact support and high-order vanishing moment characteristics of the wavelet function, the transformed coefficient matrix Most of the wavelet coefficients are zero or very small. Set the threshold η and ignore the values ​​less than η without affecting the calculation accuracy. After threshold processing, the coefficient matrix Sparse structure. When the coefficient matrix is ​​sparse enough, fast matrix-vector multiplication allows the padding of matrix elements to be kept low when solving the approximate inverse of the sparse matrix or generating the inverse operator by factorization.

[0115] f. Use the generalized successive super-relaxation iterative algorithm with preconditioning (Pre-GSOR) to solve the linear equations in the wavelet domain The iteration format is:

[0116]

[0117] or

[0118]

[0119] The preconditioner γ uses an algebraic recursive multilevel solver (ARMS) decomposition method based on block ILU decomposition to decompose the coefficient matrix Decompose to generate. n is the relaxation factor, which is calculated by minimizing the residual ||r n|| to get, take The nth iteration takes α n Make ||r n ||Minimize in Hilbert space:

[0120]

[0121] g. For the result obtained in step f Perform inverse wavelet transform,

[0122]

[0123] Then the desired single-frequency seismic wave field can be obtained;

[0124] h. Substitute the source wave field into the real wave number LS equation, calculate the scattered wave field at all detection points, and then you can get the seismic wave field at all detection points.

[0125] In order to better illustrate the implementation effect of the above specific algorithm process, a set of specific examples are given below:

[0126] Example 1

[0127] like Figure 1 The present invention is based on the wavelet domain Pre-GSOR iteration seismic wave field numerical simulation method, comprising the following steps:

[0128] a. Read in the velocity model, such as Figure 2 The model grid size is 129×65, and the spatial sampling interval is 10m. The background velocity is 1500m / s, the source wavelet is the Ricker wavelet with a main frequency of 15Hz, and the source position is r s =(640m,0m), the detection point is placed on the plane of z=0m, the channel spacing is 10m, and there are 129 channels in total.

[0129] b. Uniformly divide the integral region Ω and construct a constant interpolation basis function.

[0130] c. Construct a constant interpolation operator based on the basis function of step b.

[0131] d. Construct the normalized equivalent complex wave number LS equation. In order not to affect the calculation accuracy of the numerical results, the attenuation factor ε is taken as 0.03|O max Using the interpolation operator P N The normalized equivalent complex wave number LS equation adopts weak singular kernel The method is used for numerical discretization to obtain the linear equations corresponding to the normalized equivalent complex wave number LS equation.

[0132] e. Select the Daubechie wavelet function and construct an orthogonal wavelet transform matrix W through the Daubechie wavelet filter to satisfy WT W=WW T =I. Apply the orthogonal wavelet transform matrix W to the linear equations in step d to obtain the wavelet domain discrete equations corresponding to the normalized equivalent LS equation. Select the Daubechie12 wavelet, Figure 3 a and Figure 3 b shows the sparse pattern of the coefficient matrix in the wavelet domain after numerical discretization of the real wave number LS equation (attenuation factor ε = 0) and the normalized equivalent complex wave number LS equation at frequency f = 30 Hz, where the truncation threshold of the wavelet coefficient is set to η = Wf max *1.0e-3(Wf max Indicates the maximum value of the wavelet coefficient). For complex media, even after threshold processing, the density of the wavelet coefficient matrix of the real wave number LS equation is still 0.4936, and the number of non-zero elements is still close to 50% of the total number of matrix elements. However, under the same threshold processing principle, the density of the wavelet domain coefficient matrix of the normalized equivalent complex wave number LS equation is only 0.0887, and only about 9% of the matrix elements are non-zero values.

[0133] f. Solve the wavelet domain linear equations in step e using the Pre-GSOR iterative algorithm. Figure 4 The convergence speed of the generalized successive super-relaxation iterative series is compared when no preconditioner, multi-layer filling ILU decomposition (ILUK) preconditioner, and ARMS preconditioner are used. For the SEG / EAGE model, the iterative series does not converge when there is no preconditioner. It takes 6 iterations to reach the convergence error when using the ILUK (1) preconditioner, and 64 iterations to converge to the error range when using the ARMS preconditioner. Table 1 shows the time required to generate the above two preconditioners, the iterative calculation time, and the filling ratio of non-zero values. Although the number of iterations of the ILUK (1) preconditioner is less, the iterative calculation time required is about 1 / 5 of that of the ARMS preconditioner, but the time to generate the preconditioner is about 23 times that of the ARMS preconditioner, and the non-zero elements filled are equivalent to more than 3 times the non-zero elements of the original coefficient matrix. Considering the computational efficiency of the two preconditioner iterative calculations, the ARMS preconditioner has obvious advantages in complex strong scattering media.

[0134] Table 1

[0135] Pre condition Pre-Time(s) Iter-Time(s) Non-zero ratio Iter-number Rsd-Norm ILUK(1) 2149.48 4.18 3.24 6 1.80E-08 ARMS 93.74 19.53 0.855 64 6.30E-08

[0136] g. Solve the linear equations in the wavelet domain in step f to obtain Perform inverse wavelet transform to obtain the desired source wave field.

[0137] Figure 5a and 5b are the real and imaginary parts of the source wavefield obtained by the direct method (Singular Value Decomposition, SVD) corresponding to the real wave number LS equation. Figure 5 c and Figure 5 d are the real and imaginary parts of the source wavefield obtained by the Pre-GSOR iteration method of the normalized equivalent LS equation. It can be seen that the new method can obtain numerical results consistent with the direct method of the real wave number LS equation.

[0138] h. Substitute the single-frequency source scattered wave field calculated in step g into the real wave number LS equation and remove the direct wave to obtain the single-frequency scattered wave field at all detection points. Figure 6 a and Figure 6 b is the real and imaginary parts of the scattered wave field at the detection point at a frequency of 30 Hz. It can be seen that the numerical solution of the normalized equivalent complex wavenumber LS equation in the wavelet domain Pre-GSOR iterative method matches the numerical solution obtained by the SVD direct method of the real wavenumber LS equation highly.

[0139] The seismic wavefield forward modeling method based on wavelet domain preconditioned generalized successive super-relaxation iteration method proposed in the present invention can be further applied to seismic migration imaging and full waveform inversion to calculate the migration Green's function and forward and inversion Green's function of the full wavefield.

[0140] Note that the above are only preferred embodiments of the present invention and the technical principles used. Those skilled in the art will understand that the present invention is not limited to the specific embodiments described herein, and that various obvious changes, readjustments and substitutions can be made by those skilled in the art without departing from the scope of protection of the present invention. Therefore, although the present invention has been described in more detail through the above embodiments, the present invention is not limited to the above embodiments, and may include more other equivalent embodiments without departing from the concept of the present invention, and the scope of the present invention is determined by the scope of the appended claims.

Claims

1. A numerical simulation method of seismic wave field based on wavelet domain Pre-GSOR iteration, characterized in that: The following steps are involved: Firstly, a small imaginary component is introduced into the background wavenumber, and then the LS equation is rewritten to obtain the equivalent complex wavenumber LS equation; the rewriting of the LS equation makes the coefficient matrix have a band structure after numerical discretization, and its wavelet coefficient matrix has sparseness after threshold processing in the wavelet domain; then, the algebraic recursive multi-layer decomposition method based on block ILU decomposition is used in the wavelet domain to generate a sparse approximate inverse as a preconditioning operator to accelerate the generalized successive super-relaxation iterative method; finally, this method is used for the numerical simulation calculation of single-frequency seismic wave fields in strong scattering media.

2. The method for numerical simulation of seismic wave fields based on wavelet domain Pre-GSOR iteration according to claim 1, characterized in that: The specific steps include: a. Read the velocity model, the model size is Nx*Nz, set the source wavelet, source position, detection point position and background velocity; b. Uniform division of the integral region Ω V j is the jth partitioning element, N is the number of partitioning elements; take the partitioning constant function space L2(Ω) is the Hilbert space; basis function e j (r) Satisfy: Where r represents the spatial position on the velocity model; c. Based on the interpolation basis function e in step b j (r), construct the interpolation operator: Among them, r j Indicates V j The center of , v(r) is an arbitrary continuous function. d. Construct the equivalent complex wave number LS equation and use the interpolation operator P N The equivalent complex wave number LS equation adopts weak singular kernel The matrix form of the linear equations is obtained by numerical discretization: in, (i=1,2,…,N) is the wave field with velocity disturbance as weight at the model space point in the frequency domain, and u(r i ) are space points r i Velocity disturbance values ​​and wave field values ​​on ; is the incident wave field, s(w) represents the source wavelet, r s represents the earthquake source location, For point r i The value of the complex background Green's function on ; The weight function The complex background wave number Green operator is composed of is the complex scattering potential The diagonal operator formed by e. Select the Daubechie wavelet function, use the Mallat pyramid decomposition algorithm, and construct the orthogonal wavelet transform matrix W through the Daubechie wavelet filter to satisfy W T W=WW T =I, apply the orthogonal wavelet transform matrix W to the linear equations formula (12) in step d, and obtain the wavelet domain representation of the equivalent complex wave number LS equation discrete equations: make Then we can get the linear equations in wavelet domain According to the compact support and high-order vanishing moment characteristics of the wavelet function, the transformed coefficient matrix Most of the wavelet coefficients are zero or very small; set the threshold η, ignore the values ​​less than η without affecting the calculation accuracy, and the coefficient matrix after threshold processing It has a sparse structure. When the coefficient matrix is ​​sparse enough, the fast matrix-vector multiplication allows the filling of matrix elements to be kept at a low level when solving the approximate inverse of the sparse matrix or generating the inverse operator by factorization. f. Solve using the generalized successive super-relaxation iterative algorithm with precondition constraints The iteration format is: or The preconditioner γ uses an algebraic recursive multi-layer decomposition method based on block ILU decomposition to decompose the coefficient matrix Decompose to generate; α n is the relaxation factor, which is calculated by minimizing the residual ||r n || to get, take The nth iteration takes α n Make ||r n ||Minimize in Hilbert space: g. The iterative solution obtained in step f Perform inverse wavelet transform, Then the desired single-frequency seismic wave field can be obtained; h. Calculate the wave fields at all detection points to obtain the seismic wave fields at all detection points.

3. The method for numerical simulation of seismic wave fields based on wavelet domain Pre-GSOR iteration according to claim 2, characterized in that: Step d specifically comprises the following steps: d1. Define the velocity disturbance function v0 is the constant background velocity, v(r) is the medium velocity, and according to the wave field superposition principle, the Helmholz equation satisfied by the scalar wave equation in the frequency domain is rewritten as an equivalent source form: Where ω is the angular frequency, s(ω)δ(rr s ) indicates excitation at r s The point source at , u(r) is the wave field to be determined. Add iεu(r) (ε>0) to both ends of formula (3), and we get: Pick (0≤a<1.0), where |O| max =max r |O(r)|; d2. Define complex background wave number Then there is the complex wave number background Green's function satisfy: in, r' is the source point of Green's function. The analytical solution is: In numerical calculations, according to the disturbance intensity of the model, it is necessary to select an appropriate constant a to ensure Greater than 1.0e-6 within the effective calculation area; d3, corresponding to the complex background wave number Complex scattering potential Based on Green's theorem, the solution of formula (4) can be expressed in integral form to obtain the complex wave number LS equation: in, At this time, the integral kernel in the complex wave number LS equation is removed In addition, it also contains the complex velocity disturbance term; d4. Define the wave field with velocity perturbation as weight The equivalent complex wave number LS equation is obtained: The above formula is an integral equation about u'(r), By normalizing formula (8) point by point as the normalization function, the normalized equivalent LS equation can be obtained: d5. Define the integral operator Using the interpolation operator P N , replace the continuous part of the integrand with its interpolation function, then the integral operator K Approximate K N It can be expressed as: Among them, the weight function Then the normalized equivalent complex wave number of the LS equation is The approximate equations can be expressed as: The matrix (operator) of the above formula is expressed as: Where I is the identity matrix, The weight function The complex wave number Green operator is composed of is the complex scattering potential The diagonal operator is constructed; by selecting an appropriate attenuation factor ε, the complex velocity perturbation The imaginary component in is not zero, even at the point where the velocity disturbance O(r) is zero, the complex velocity disturbance The amplitude of is never zero.

Images