Method and device for forward modeling of seismic wave field under undulating surface

Abstract

The invention relates to the technical field of seismic exploration, and particularly discloses a seismic wave field forward modeling method and device under an undulating surface, and the method comprises the steps: obtaining a surface model function of the undulating surface; embedding the surface model function into the difference grid model by using an immersed boundary method to obtain a finite difference grid model; calculating space positions of spanned nodes, pseudo spanned nodes, corresponding projection points and fractional nodes in the finite difference grid model; calculating the wave field value of the global grid node based on a finite difference grid model to obtain the wave field value of each node in the global integer grid nodes; utilizing a parabola interpolation method to determine wave field values at fractional nodes corresponding to the spanning nodes and the pseudo spanning nodes; and on the basis of transmission absorption boundary conditions, according to the wave field values of the fractional nodes corresponding to the spanned node and the pseudo spanned node, updating the wave field values at the spanned node and the pseudo spanned node to obtain a seismic wave forward modeling result. According to the scheme, the accuracy and efficiency of the seismic wave field forward modeling method are improved.

Description

Technical Field

[0001] The present invention relates to the field of seismic exploration technology, and in particular to a forward modeling method and device for seismic wave fields beneath an undulating surface. Background Art

[0002] Wave equation numerical simulation is a simulation method that uses numerical methods to solve wave equations when the structure and characteristics of the underground medium are known, in order to understand the propagation conditions and laws of seismic waves in the underground medium. High-precision and high-efficiency numerical simulation methods affect the accuracy and efficiency of methods such as wave equation migration imaging and full waveform inversion. At present, wave equation numerical simulation has been widely used in industries such as oil, natural gas, and engineering environmental geophysics, among which the finite difference method is the earliest numerical simulation method. The finite difference method uses differential formulas to approximate differential operators and discretizes the original equations, thereby converting continuous problems into discrete problems. After decades of development and research, the finite difference method has become the most widely used mainstream method for wave equation numerical simulation due to its advantages of low computational cost, high computational efficiency, and simple implementation.

[0003] As oil and gas exploration targets become increasingly complex, the industry has placed higher demands on the accuracy of seismic exploration technology. In actual production, the undulating surface is a very important influencing factor in the numerical simulation of wave equations. The traditional finite difference method usually uses rectangular grids to mesh the calculation area in a Cartesian coordinate system. When the interface changes undulatingly, the problem of step approximation will inevitably occur when the rectangular grid is used to approximate the undulating interface, resulting in the appearance of false scattering. In order to suppress this false scattering in the traditional finite difference method, the grid needs to be refined, that is, the large grid is divided into fine grids with extremely small grid spacing. However, the refinement will lead to an increase in the amount of computation for numerical simulation and an increase in computational storage, which seriously affects the computational efficiency.

[0004] Different from the rectangular grid approximation of the undulating surface, the curved grid method uses a curved grid to fit the undulating surface. At the same time, in order to facilitate the calculation using finite differences, the wave equation expressed in the Cartesian coordinate system needs to be mapped accordingly to obtain the wave equation in the curvilinear coordinate system, and then the finite difference method is used to solve it. This idea has been adopted by many scholars and has been studied accordingly. Although the curved grid method can completely fit the smooth undulating surface and can better handle the step approximation problem of the rectangular grid approximation of the inclined interface, this method also has shortcomings. Since the wave equation needs to be mapped from the Cartesian coordinate system to the curvilinear coordinate system, the number of partial derivatives contained in the wave equation after mapping increases, which directly leads to an exponential increase in the amount of calculation and consumes additional calculation costs. In addition, when the surface undulation changes dramatically, it is difficult to obtain stable calculation results using the curved grid method. Later, some scholars conducted research on this problem. Although a finite difference format with good stability after the wave equation mapping was obtained, the calculation accuracy of the difference format was very limited.

[0005] In addition, the immersed boundary method is also an effective means to deal with the problem of undulating surfaces. This method was first applied in fluid dynamics to simulate cardiac dynamics and related blood flow. At present, the immersed boundary method is usually used with Cartesian grids for fluid modeling, and was later introduced into the field of seismic modeling. At present, some scholars have used the immersed boundary method to realize the free boundary conditions of undulations. The immersed boundary method directly embeds the irregular interface into the traditional Cartesian coordinate system. The undulating surface is located on the fractional grid, while the calculation variables in the wave equation are located on the integer grid. The wave equation is controlled in the Cartesian coordinate system. The physical quantities on the undulating surface can be obtained by interpolation and other methods. Compared with the traditional finite difference method and the curved grid method, the immersed boundary method avoids the step approximation problem while having the advantages of high computational efficiency and better adaptability to the undulating surface.

[0006] In actual seismic exploration, the propagation of artificially excited seismic waves in underground media can be regarded as propagation in an infinite medium. However, in the numerical simulation of wave equations, due to the limited computing power and storage capacity of current computers, it is impossible to simulate in an infinite medium, and numerical simulation can only be performed in a limited space, so there will be the problem of artificial truncation of boundaries. Truncated boundaries will produce false reflections, which in turn affect the real reflection waves from the stratum interface. Therefore, absorbing boundary conditions are usually set at the boundaries to absorb and attenuate boundary reflections to ensure the validity of the internal waveform.

[0007] At present, when simulating seismic wave fields under undulating surface conditions, if absorbing boundary conditions are set at the boundary, it is generally to simulate the free boundary conditions of the undulating surface. This boundary condition exists physically, and when seismic waves are incident on the free surface, complex wave field phenomena will be generated. However, when studying some seismic imaging methods, such as reverse time migration and full waveform inversion, most of them regard the reflections related to the free surface as noise, which will reduce the accuracy of migration imaging and inversion.

[0008] To address the above issues, no finite difference scheme that includes absorbing boundary conditions for undulating surfaces has been proposed so far. Summary of the invention

[0009] The embodiments of this specification provide a forward modeling method and device for seismic wave fields under undulating terrain to improve the simulation accuracy and computational efficiency of the seismic wave field simulation method under undulating terrain conditions in the prior art.

[0010] The embodiment of this specification provides a forward modeling method for seismic wave fields under an undulating surface, comprising:

[0011] Obtaining a surface model function of an undulating surface; embedding the surface model function into a differential grid model using an immersed boundary method to obtain a finite difference grid model under the undulating surface; calculating ghost nodes and pseudo-ghost nodes in the finite difference grid model under the undulating surface and spatial positions of projection points and fractional nodes corresponding to the ghost nodes and pseudo-ghost nodes;

[0012] The wave field values ​​of the global grid nodes are calculated based on the finite difference grid model to obtain the wave field values ​​of each node in the global integer grid nodes; the wave field values ​​at the fractional nodes corresponding to the ghost nodes and pseudo ghost nodes are determined by using a parabolic interpolation method;

[0013] Based on the transmission absorption boundary condition, the wave field values ​​at the ghost nodes and pseudo ghost nodes are updated along the surface normal direction according to the wave field values ​​of the fractional nodes corresponding to the ghost nodes and pseudo ghost nodes to obtain the seismic wave forward modeling results.

[0014] The embodiment of this specification also provides a forward modeling device for seismic wave fields under an undulating surface, comprising:

[0015] An acquisition module is used to acquire a surface model function of an undulating surface; embed the surface model function into a differential grid model using an immersed boundary method to obtain a finite difference grid model under the undulating surface; calculate ghost nodes and pseudo-ghost nodes in the finite difference grid model under the undulating surface and spatial positions of projection points and fractional nodes corresponding to the ghost nodes and pseudo-ghost nodes;

[0016] A calculation module is used to calculate the wave field value of the global grid node based on the finite difference grid model to obtain the wave field value of each node in the global integer grid node; and to determine the wave field value at the fractional node corresponding to the ghost node and the pseudo ghost node by using a parabolic interpolation method;

[0017] The updating module is used to update the wave field values ​​at the ghost nodes and pseudo ghost nodes along the surface normal direction according to the wave field values ​​of the fractional nodes corresponding to the ghost nodes and pseudo ghost nodes based on the transmission absorption boundary conditions to obtain seismic wave forward simulation results.

[0018] The present specification also provides a computer device including a processor and a memory for storing processor executable instructions, wherein the processor implements the steps of the forward modeling method for seismic wave fields beneath an undulating surface as described in any of the above embodiments when executing the instructions.

[0019] The embodiments of the present specification also provide a computer-readable storage medium on which computer instructions are stored. When the instructions are executed, the steps of the forward modeling method of the subsurface seismic wave field described in any of the above embodiments are implemented.

[0020] In an embodiment of the present specification, a method for forward modeling of seismic wave fields under an undulating surface is provided, which can obtain a surface model function of the undulating surface; embed the surface model function into a differential grid model using an immersed boundary method to obtain a finite difference grid model under the undulating surface; calculate ghost nodes and pseudo-ghost nodes in the finite difference grid model under the undulating surface and the spatial positions of projection points and fractional nodes corresponding to the ghost nodes and pseudo-ghost nodes; calculate the wave field values ​​of global grid nodes based on the finite difference grid model to obtain the wave field values ​​of each node in the global integer grid nodes; determine the wave field values ​​at the fractional nodes corresponding to the ghost nodes and pseudo-ghost nodes using a parabolic interpolation method; based on a transmission absorption boundary condition, update the wave field values ​​at the ghost nodes and pseudo-ghost nodes along the surface normal direction according to the wave field values ​​of the fractional nodes corresponding to the ghost nodes and pseudo-ghost nodes to obtain a seismic wave forward modeling result. In the above scheme, the immersed boundary method is used to embed the terrain directly into the Cartesian coordinate system, so that the relevant physical quantities can be calculated directly in the coordinate system, which is conducive to the implementation of the finite difference method and avoids the problem of step approximation of the traditional finite difference method. By combining the immersed boundary method and the transmission absorption boundary condition, the undulating surface can be effectively processed, which has the advantages of good absorption effect and lower calculation cost while better adapting to the undulating changes of the surface. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] The drawings described herein are used to provide a further understanding of this specification, constitute a part of this specification, and do not constitute a limitation of this specification. In the drawings:

[0022] Figure 1 A schematic diagram of the principle of Liao's transmission boundary (hereinafter referred to as Liao boundary) in an embodiment of this specification is shown;

[0023] Figure 2 It shows u in one embodiment of this specification n (x-sΔx) interpolation diagram;

[0024] Figure 3 It shows u in one embodiment of this specification n (x-2sΔx) interpolation diagram;

[0025] Figure 4 A schematic diagram of Liao boundary implementation in a spatial 2L-order finite difference simulation in an embodiment of this specification is shown;

[0026] Figure 5 A schematic diagram of a two-dimensional immersed absorbing boundary condition in an embodiment of this specification is shown;

[0027] Figure 6 A schematic diagram showing the relationship between ghost nodes and surface functions in an embodiment of this specification;

[0028] Figure 7 A schematic diagram showing a fractional node in an embodiment of the present specification being above the ground surface;

[0029] Figure 8 A schematic diagram of a pseudo ghost node in an embodiment of this specification is shown;

[0030] Fig. 9 A schematic diagram of a three-dimensional immersion absorption boundary condition in an embodiment of this specification is shown;

[0031] Fig.10 A schematic diagram of a ghost node in an embodiment of this specification is shown;

[0032] Fig.11 A schematic diagram of fractional node interpolation in an embodiment of this specification is shown.

[0033] Fig.12 A schematic diagram of a three-dimensional pseudo ghost node in an embodiment of this specification is shown;

[0034] Fig.13 A schematic diagram of a two-dimensional undulating surface model in an embodiment of this specification is shown;

[0035] Fig.14 A schematic diagram showing the absorption effect of a two-dimensional immersed absorption boundary in an embodiment of this specification is shown;

[0036] Fig.15 A schematic diagram of a three-dimensional undulating surface model in an embodiment of this specification is shown;

[0037] Fig.16 A schematic diagram showing the absorption effect of a three-dimensional immersed absorption boundary in an embodiment of this specification is shown;

[0038] Fig.17 A schematic diagram showing the absorption effect of a three-dimensional immersed absorption boundary in an embodiment of this specification is shown;

[0039] Fig.18 A flow chart showing a forward modeling method for a seismic wave field beneath an undulating surface in an embodiment of the present specification is shown;

[0040] Fig.19 A schematic diagram showing a forward modeling device for seismic wave fields beneath an undulating surface in an embodiment of the present specification;

[0041] Fig. 20 A schematic diagram of a computer device in an embodiment of the present specification is shown. DETAILED DESCRIPTION

[0042] The principles and spirit of this specification will be described below with reference to several exemplary embodiments. It should be understood that these embodiments are provided only to enable those skilled in the art to better understand and implement this specification, and are not intended to limit the scope of this specification in any way. On the contrary, these embodiments are provided to make this specification more thorough and complete, and to fully convey the scope of this disclosure to those skilled in the art.

[0043] Those skilled in the art will appreciate that the embodiments of this specification may be implemented as a system, device, method, or computer program product. Therefore, this specification may be implemented in the following forms: complete hardware, complete software (including firmware, resident software, microcode, etc.), or a combination of hardware and software.

[0044] When forward modeling the seismic wave field under the undulating surface, absorbing boundary conditions can be used. According to the principle of suppressing the reflected wave energy, absorbing boundary conditions can be divided into two types. The first type of absorbing boundary conditions forces the seismic waves propagating from the source direction to be absorbed and attenuated at the boundary of the calculation area to suppress the reflected wave field at the boundary. The most representative absorbing boundary conditions are perfect matching layer (PML), split perfect matching layer (SPML) and complex frequency shift convolution perfect matching layer (CFS-CPML). Although these absorbing boundary conditions have relatively good absorption and attenuation effects, the computational cost required for three-dimensional seismic simulation is still high, and the PML-type boundary needs to use the grid point wave field values ​​on both sides of the boundary to absorb and attenuate the wave field at the boundary, which is not suitable for processing undulating surfaces. Therefore, the PML-type absorbing boundary conditions still need to be further improved. The second type is to set an impedance matching absorbing boundary condition at the boundary to simulate the propagation of seismic waves through the boundary but without reflection, which is called a transmission absorbing boundary condition.

[0045] At present, when simulating the seismic wave field under the condition of undulating surface, the free boundary condition of undulation is generally simulated. This boundary condition exists physically. When the seismic wave is incident on the free surface, a complex wave field phenomenon will be generated. However, when studying some seismic imaging methods, such as reverse time migration and full waveform inversion, most of the reflections related to the free surface are regarded as noise, which will reduce the accuracy of migration imaging and inversion. For this reason, a forward simulation method of seismic wave field under undulating surface is proposed in the embodiment of this specification. The method uses the immersed boundary method to realize the absorption boundary condition under the condition of undulating surface, so as to suppress the boundary reflection from the undulating surface, and provide a basis for further research and development of seismic imaging and inversion methods. The transmission absorption boundary condition (Liao boundary) based on the plane wave propagation principle adopted in the embodiment of this specification belongs to the second category, that is, the transmission absorption boundary condition. Compared with the PML type absorption boundary condition, the Liao boundary has the advantages of low computational cost, high computational efficiency and simple and easy implementation, and only needs to use the wave field value of the grid point on one side of the undulating surface, which can better adapt to the undulating surface.

[0046] In the embodiments of this specification, a method for realizing an absorbing boundary in numerical simulation of wave equations on undulating surfaces is proposed by combining a transmission absorbing boundary based on the plane wave principle and an immersed boundary method. The technical solution in the embodiments of this specification has the following advantages: 1) The immersed boundary method is used to embed the terrain directly into the Cartesian coordinate system, so that the relevant physical quantities can be calculated directly in the rectangular coordinate system, which is conducive to the implementation of the finite difference method and avoids the problem of step approximation; 2) While having a good absorption effect, it can adapt to the situation where the surface undulation changes more dramatically, and has low calculation cost and high calculation efficiency.

[0047] Based on this, an embodiment of the present specification provides a forward modeling method for seismic wave fields beneath an undulating surface.

[0048] First, in the embodiments of this specification, the boundary is used to realize the update of the wave field at the boundary. The basic principle of the boundary is as follows: Figure 1 shown.

[0049] like Figure 1 As shown, taking the right boundary as an example (the same applies to the following text), a plane wave with a velocity of c is incident along a direction that makes an angle θ with the x-axis and reaches the boundary point Q (assuming its horizontal coordinate is x) at time t+Δt. According to the plane wave propagation theory, under the condition of isotropic perfectly elastic medium, the wave field at point Q at time t+Δt is equal to the wave field at point P on the wave propagation path at time t, and the wave field at point P is equal to the wave field at point P' on the same wave front, that is:

[0050]

[0051] In the formula, u represents the wave field, Δt is the time interval. It can be known from formula (1) that the wave field value at the boundary point Q at time t+Δt can be predicted by the wave field value at point P' on the x-axis. However, the angle θ in formula (1) is unknown, so it cannot be calculated directly. In order to avoid explicitly obtaining the incident angle θ, the approximate N-order plane wave wave field prediction formula can be obtained by establishing the error transfer formula of the wave field derivative:

[0052]

[0053] Wherein, α is a transmission coefficient set artificially, in the embodiment of this specification, α=1, N is the transmission order, is the binomial coefficient:

[0054]

[0055] In the embodiments of this specification, the transmission order N=2 is set, which is called the second-order material boundary formula. Assuming that the grid spacing in the x direction is Δx, the second-order material boundary can be expressed as:

[0056] u n+1 (x)≈2u n (x-sΔx)-u n-1 (x-2sΔx) (4)

[0057] In the formula, s = αcΔt / Δx, and the superscripts n+1, n, and n-1 represent three adjacent moments. The finite difference numerical simulation of the wave equation generally uses a regular rectangular grid for partitioning. For the two-dimensional constant density acoustic wave equation shown in formula (5), the time partial derivatives in the wave equation are usually discretized using the second-order central difference format, and the explicit three-step recursive format shown in formula (6) can be obtained. Formula (6) is similar to the three-step iterative format in formula (4), so the use of the second-order Liao boundary will not increase the memory consumption.

[0058]

[0059] Where φ(t) is the source wavelet, x = (x, z), x s is the source location. Using the second-order time and 2L-order space central difference operators to approximate the second-order partial derivatives, the difference equation can be obtained:

[0060]

[0061] In the formula, It represents the wave field after grid discretization of time (t) and two-dimensional space (Cartesian coordinate system xoz), n represents the nth time grid node, i represents the i-th grid node in the x direction, and j represents the j-th grid node in the z direction. h is the grid spacing in the x and z directions, (i s ,j s ) is the grid node where the earthquake source is located, δ is the Dirichlet-Delta function used at the earthquake source location, L is the difference operator length, a l (l=0,1,2,…,L) is the differential coefficient, which can be determined using Taylor approximation or least squares method.

[0062] Generally, the explicit finite difference method is used to solve the wave equation. Due to the limitation of stability conditions, it is known that s < 1. At this time, u in equation (4) n (x-sΔx) and u n-1 (x-2sΔx) is located on the fractional grid node. The wave field values ​​of these two items need to be obtained by interpolation using the wave field values ​​of the integer nodes. n (x-sΔx), the specific interpolation process is as follows Figure 2 shown.

[0063] Figure 2 The black “×” in the figure is the score node P to be found. 0 , whose wave field value uses three adjacent integer nodes A 0 ,A 1 ,A 2 It is obtained by parabolic interpolation as shown in formula (7),

[0064] u n (x-sΔx)=T 1,1 u n (x)+T 1,2 u n (x-Δx)+T 1,3 u n (x-2Δx) (7)

[0065]

[0066] The second term u in formula (4) n-1 The interpolation process of (x-2sΔx) is as follows Figure 3 As shown in the figure, the black "*" is the score node F to be found 0 , the solution of its wave field value requires two parabolic interpolations. The first one uses the three adjacent integer nodes to interpolate to obtain the fractional node P 0 ,P 1 ,P 2 The wave field value at the first step is the same as Figure 2 Consistent; the second use of point P 0 ,P 1 ,P 2 The wave field value at is interpolated to obtain F 0 It should be noted that the relative spatial position relationship between the interpolation node and the point to be interpolated is the same during the first and second interpolation, so the weights of the two interpolations are the same, as shown in formula (8).

[0067] Combining equation (7) and equation (8) and Figure 2 and Figure 3 , the second term u in formula (4) n-1 (x-2sΔs) can be expressed as:

[0068]

[0069] Combining the above formula derivation, we can obtain the wave field update formula of the second-order material boundary:

[0070]

[0071] It can be seen from formula (10) that to calculate the wave field value of the point located at the boundary, it is only necessary to use the wave field values ​​of the four points inside the boundary and located on the x-axis.

[0072] The second-order Liao boundary is suitable for the second-order difference operator. At this time, only one boundary point needs to use the Liao boundary formula (10) for wave field update. When the difference operator length is 2L, there will be L grid nodes in the absorption area, and their wave field values ​​can also be updated using formula (10). In order to further improve the absorption effect of the absorption area, this manual sets a grid number of L outside the calculation area except for the undulating surface. c (L c ≥L) buffer, such as Figure 4 shown.

[0073] Figure 4 Midpoint A 0 are the grid nodes located in the inner region, is the expanded boundary point. Figure 4 The wave field update steps of the mid-boundary points are:

[0074] 1) Use formula (6) to calculate the wave field value of the global grid node. When approaching the boundary, the strategy of reducing the difference order is adopted to avoid the problem that the difference template exceeds the model boundary. The calculated wave field value at this time is recorded as u 1 ;

[0075] 2) The wave field value of the grid node located in the absorption area is updated using formula (10), and the calculated wave field value is recorded as u 2 ;

[0076] 3) For u 1 and u 2 Perform weighted summation, and the calculation result is used as the boundary point The wave field value of:

[0077] u(A q )=(1-ω q ) 2 (A q )+ω q u 1 (A q ) (11)

[0078]

[0079] Next, the implementation of the two-dimensional immersed absorbing boundary condition is explained.

[0080] The immersed boundary method can be used to embed the undulating surface directly into the Cartesian coordinate system. When the finite difference method is used to calculate the wave field value at the grid points close to the surface, the nodes located above the surface will be used according to the difference operator length adopted. These nodes located above the surface and required by the internal grid nodes for finite difference calculation are called ghost nodes, such as Figure 5 shown.

[0081] Figure 5 In the figure, the black “×” point F is the difference center point of the fourth-order central difference in space, the black solid square point is the grid node used for the difference calculation under the surface, the black solid circle point G is the ghost node, the white solid circle point G' is the projection point of the ghost node on the surface, which is defined as the point with the shortest distance from the ghost node to the surface, and the white solid circle point A 1 ,A 2 ,A 3 ,A 4 are four fractional nodes with equal spacing between the ghost node and the normal direction of the surface. When using the fourth-order central difference to calculate the wave field value of point F, not only the wave field values ​​of the square solid grid nodes below the surface are needed, but also the wave field values ​​of the black dots above the surface, namely the ghost nodes. Assuming that the surface can be expressed as z = f(x), then the ghost node (x g ,z g ) can be defined in two-dimensional space as follows:

[0082]

[0083] Where Δz is the grid spacing in the z direction, L x and L z They represent the differential operator lengths in the x-direction and z-direction respectively. Formula (13) corresponds to Figure 6 The three situations in .

[0084] In this embodiment, the positions of all ghost nodes in space are first determined using equation (13), and then the spatial positions of all ghost nodes corresponding to the projection points on the surface are further determined. Each pair of ghost nodes and projection points can determine the internal normal direction ( Figure 5 The solid arrow in the figure is then used to apply boundary conditions along the internal normal direction. Figure 5 As shown, in the internal normal direction, min(Δx,Δz) is selected as the distance to determine the four fractional nodes A with equal spacing. 1 ,A 2 ,A 3 ,A 4 , the wave field update formula of the ghost node can be obtained from the absorption area wave field update formula (10):

[0085]

[0086] In formula (14), the fractional node A 1 ,A 2 ,A 3 ,A 4 The wave field value of needs to be obtained through parabolic interpolation calculation. Figure 5 As shown, the score node A 1The wave field value of is obtained by two-dimensional parabola interpolation of 16 grid nodes in the 4×4 dotted grid in the figure, and the wave field values ​​of the remaining fractional nodes are obtained in a sampled manner.

[0087] The number of ghost nodes above the surface is related to the difference operator length. The longer the difference operator length, the more ghost nodes there are. It may happen that not all the fractional nodes in the internal normal direction are below the surface, such as Figure 7 As shown, the wave field value of the fractional node will be different at this time.

[0088] Figure 7 The black solid circle G is the ghost node, the white solid circle G' is the projection point of the ghost node on the surface, the black "×" F point is the center point of the spatial fourth-order difference operator, and the black solid square point is the grid node required for difference calculation. 1 ,A 2 ,A 3 ,A 4 ,A 5 The fractional nodes are equally spaced in the internal normal direction. 1 When it is located above the surface, the wave field value is no longer obtained by interpolation calculation, but is obtained by fractional node A according to formula (10): 2 ,A 3 ,A 4 ,A 5 Sure:

[0089]

[0090] The TRON value of the ghost node is still determined by the score node A. 1 ,A 2 ,A 3 ,A 4 Calculated using formula (14), formula (15) shows that when the fractional node is located above the surface, the wave field value is solved in the same way as the wave field value of the ghost node. The wave field value of the fractional node below the surface still uses the interpolation method. It should be noted that when there are m fractional nodes above the surface, m fractional nodes need to be extended at equal intervals in the internal normal direction to satisfy the wave field update formula of formula (14).

[0091] When interpolating the wave field values ​​of the fractional nodes below the surface, due to the difference operator length and interpolation window size, the wave field values ​​of some grid nodes above the surface but not ghost nodes need to be used in the interpolation process. Such grid nodes are called pseudo ghost nodes, such as Figure 8 shown.

[0092] Figure 8The black solid dot PG in the middle represents a pseudo ghost node, the black solid dot PG' is the projection point of the pseudo ghost node on the surface, and the white solid dot B l (l=1,2,,6) are pseudo ghost nodes and fractional nodes with equal spacing in the direction of the internal normal of the surface. The black solid squares are all ghost nodes. The black “*” point A is the fractional node on the internal normal of the ghost node. The white solid squares are the internal grid nodes required for differential calculation. Figure 8 As shown, the spatial fourth-order difference format is adopted, and the interpolation window size is 4×4. When calculating the wave field value of the fractional node A, in addition to the ghost nodes and the grid nodes below the surface in the interpolation window, the wave field value of the pseudo ghost node PG is also required. In the interpolation process, if the wave field of the pseudo ghost node PG is not updated, its wave field value is always zero, which will affect the accuracy of the interpolation result and further affect the stability of the numerical simulation. Therefore, in the present invention, the pseudo ghost node wave field value update method is set to be consistent with the ghost node wave field update, that is, the wave field value of the PG point is calculated by the fractional node B. 1 ,B 2 ,B 3 ,B 4 Sure:

[0093]

[0094] The fraction of nodes B located above the ground 1 and B 2 The wave field values ​​are respectively determined by the fractional node B 2 ,B 3 ,B 4 ,B 5 and B 3 ,B 4 ,B 5 ,B 6 The wave field update principle is the same as Figure 7 Middle A 1 Points are consistent. Fractional node B below the surface 3 ,B 4 ,B 5 ,B 6 The wave field value of is still obtained by interpolation calculation.

[0095] Next, the implementation of the three-dimensional immersion absorbing boundary condition is described.

[0096] The three-dimensional second-order constant density acoustic wave equation can be approximately expressed as:

[0097]

[0098] Where: x = (x, y, z), x s is the source location. Using the second-order time and 2L-order space central difference operators to approximate the second-order partial derivatives, the difference equation can be obtained:

[0099]

[0100] Where: It represents the wave field after grid discretization of time (t) and three-dimensional space (Cartesian coordinate system oxyz), n represents the nth time grid node, i represents the i-th grid node in the x-direction, j represents the j-th grid node in the y-direction, and k represents the k-th grid node in the z-direction. n is the number of time steps, h is the grid spacing in the x, y and z directions, (i s ,j s ,k s ) is the grid node where the earthquake source is located, δ is the Dirichlet-delta function used at the earthquake source location. In order to save computational cost, the explicit three-step recursive format similar to equation (18) is maintained, and the second-order Liao boundary is still used to update the absorption region wave field in three-dimensional space.

[0101] In three-dimensional space, the immersed boundary method is used to embed the undulating surface directly into the spatial grid, such as Fig. 9 shown.

[0102] Fig. 9 The black five-pointed star F represents the center point of the second-order spatial difference, the black solid dot G represents the ghost node, the gray surface S represents the undulating ground, the white solid dots are the grid nodes below the surface required for the second-order spatial central difference, and the black solid dot A is the fractional node in the direction of the ghost node and the internal normal of the surface. Assuming that the terrain can be represented by the function z=f(x,y), then the ghost node (x g ,y g ,z g ) can be expressed as:

[0103]

[0104] Where L x ,L y ,L z represents the difference operator length in the three directions of x, y, and z, respectively, and Δx, Δy, and Δz represent the grid spacing in the three directions. Fig.10 The relative position relationship between the ghost nodes and the ground surface is shown.

[0105] Fig.10 The black solid dot G in the middle represents a ghost node, and the surface S represents the ground surface.

[0106] Fig.10In the above example, the wavefield value of the ghost node is still calculated by the wavefield update formula (10) of the second-order Liao boundary, and the spacing between the fractional nodes is selected as min(Δx, Δy, Δz). If all the fractional nodes A located in the direction of the ghost node and the internal normal of the surface are located below the surface, the number of fractional nodes is 4. If there are m fractional nodes located above the surface, m fractional nodes need to be expanded in the direction of the internal normal. There are 4+m fractional nodes on the internal normal (excluding ghost nodes). The wavefield value of the fractional node located above the surface is calculated by the 4 adjacent fractional nodes on its inner side according to formula (10). The wavefield value of the fractional node located below the surface is calculated by interpolation. The interpolation process is as follows: Fig.11 shown.

[0107] Fig.11 The middle left picture is Fig.11 The gray cube in the center of the cube in the right picture, Fig.11 The wave field value of the fractional node A in the left figure adopts Fig.11 The wave field values ​​of the grid nodes contained in the 4×4×4 three-dimensional cube in the right figure are calculated based on parabolic interpolation. In the interpolation process, the wave field values ​​of pseudo ghost nodes may still be used. The wave field values ​​of pseudo ghost nodes in three-dimensional space are as follows: Fig.12 As shown in Figure 2, the calculation of the wave field value of the pseudo ghost node is still consistent with the two-dimensional case, that is, the pseudo ghost node is regarded as a ghost node, and its wave field is updated in the same way.

[0108] Fig.12 The black solid dot PG is a pseudo ghost node, the red solid dot B is a pseudo ghost node and a fractional node in the direction of the internal normal of the surface, the white solid dot G is a ghost node, the surface S represents the surface, the black five-pointed star point F is the center point of the spatial second-order central difference operator, and the white solid dots are the grid nodes below the surface that need to be used when performing finite difference calculations on F.

[0109] In summary, the specific implementation steps of the forward modeling method for seismic wave fields under undulating surfaces (also known as the two-dimensional and three-dimensional immersed absorbing boundary condition implementation method) in the embodiments of this specification are as follows:

[0110] Step 1: Given the surface model function, the terrain function is embedded into the differential grid model using the immersed boundary method, and the spatial positions of all ghost nodes are determined using equation (13) or (19).

[0111] Step 2: Take the point with the shortest distance from the ghost node to the ground as the projection point, find the corresponding projection points of all ghost nodes on the ground, and determine the internal normal direction of the ground.

[0112] Step 3: along the internal normal direction, with the ghost node as the first point, determine four fractional nodes (excluding the ghost node) at equal intervals with a distance of min(Δx, Δz) or min(Δx, Δy, Δz).

[0113] Step 4: Determine whether there are fractional nodes located above the ground surface. If so, expand the same number of fractional nodes in the normal direction and proceed to step 5. If not, proceed directly to step 5.

[0114] Step 5: Set a 4×4 interpolation window for each fractional node, and the size of the three-dimensional interpolation cube window is 4×4×4.

[0115] Step 6, determine whether there are pseudo ghost nodes in the interpolation window or interpolation cube. If so, repeat steps 2 and the following. If not, proceed directly to step 7.

[0116] Step 7: Use parabolic interpolation to calculate the wave field value of each fractional node, and then use equations (14) and (16) to update the wave field values ​​of ghost nodes and pseudo ghost nodes in the space.

[0117] The method in the embodiments of this specification proposes a method for realizing an absorption boundary in the numerical simulation of wave equations on undulating surfaces by combining the transmission absorption boundary and the immersed boundary method based on the plane wave principle. The above technical scheme has the following advantages: 1) The immersed boundary method is used to embed the terrain directly into the Cartesian coordinate system, so that the relevant physical quantities can be calculated directly in the rectangular coordinate system, which is conducive to the implementation of the finite difference method and avoids the problem of step approximation; 2) While having a good absorption effect, it can adapt to the situation where the surface undulation changes more dramatically, and the calculation cost is lower and the calculation efficiency is higher.

[0118] The above method is described below in conjunction with several specific embodiments. However, it should be noted that the specific embodiments are only for better illustrating the present specification and do not constitute an improper limitation on the present specification.

[0119] Example 1. Two-dimensional uniform velocity model calculation example

[0120] The number of grids in the two-dimensional simulation area is 116×425, the grid spacing is 10m, the source is the Ricker wavelet with a main frequency of 10Hz, which is excited at the grid point (44,212) with a velocity of 2000m / s. The finite difference method with second-order time and fourth-order space accuracy is used, and the undulating surface function is:

[0121]

[0122] Among them, (x,z) represents the coordinates of the spatial position.

[0123] Please refer to Fig.13, showing a two-dimensional undulating surface model. Fig.13 (a) is the established sinusoidal undulating surface model; Fig.13 (b) in Fig.13 A partially enlarged surface model in (a). Fig.13 In (b), the black dots “o” represent ghost nodes, the black dots “x” represent pseudo ghost nodes, the black “*” represents the fractional nodes along the internal normal direction of the ghost nodes, the black “+” represents the fractional nodes along the internal normal direction of the pseudo ghost nodes, the solid dots on the surface curve represent the projection points of the ghost nodes and pseudo ghost nodes on the surface, and the black dotted arrows represent the internal normals of the ghost nodes and pseudo ghost nodes and the surface.

[0124] In the simulation area, the left, right and bottom boundaries are set as second-order material boundaries as absorbing boundaries. Fig.14 , which shows a schematic diagram of the absorption effect of a two-dimensional immersed absorption boundary in an embodiment of the present specification. Fig.14 (a) is a snapshot of the wave field at the 600th time step in the uniform velocity model with no immersed absorbing boundary conditions applied to the undulating surface. The arrows in the figure indicate the reflected waves from the upper surface. Fig.14 (b) is a snapshot of the wave field at the 600th time step in the uniform velocity model, where the immersed absorbing boundary condition in this embodiment is applied on the undulating surface. The absorption effect of the immersed absorbing boundary condition in this embodiment on the undulating surface is shown in FIG. Fig.14 As shown in (b), no obvious reflection wave is observed, which is different from the case where no immersion absorbing boundary condition is applied at the same time. Fig.14 Compared with (a) in FIG. 1 , the reflected wave from the ground surface is well absorbed. This example proves that the method in the embodiment of this specification has a good absorption effect when the ground surface is undulating.

[0125] Example 2. Calculation example of two-dimensional layered velocity model

[0126] The effectiveness and feasibility of the present invention are further demonstrated using a two-layer medium. The number of grids in the simulation area is still 116×425, the grid spacing is 10m, the source is a Ricker wavelet with a main frequency of 10Hz excited at the grid point (44,212), and the finite difference method with second-order time and fourth-order space accuracy is used. The velocity of the upper layer of the two-layered medium is 2000m / s, and the velocity of the lower layer is 3000m / s. The undulating surface function is:

[0127]

[0128] Among them, (x,z) represents the coordinates of the spatial position.

[0129] Please continue to refer to Fig.14 , Fig.14(c) is the layered velocity model used in the examples of this specification. Fig.14 (d) is a wavelength snapshot at the 600th time step when the undulating surface in the layered velocity model does not use the immersed absorbing boundary condition in this embodiment. Fig.14 The arrows in (d) indicate the reflected waves from the upper surface and the interlayer reflected waves from the velocity interface, respectively. Fig.14 (e) is a wave field snapshot at the 600th time step when the immersed absorbing boundary condition in this embodiment is used on the undulating surface in the layered velocity model.

[0130] The velocity model is Fig.14 As shown in (c), the left, right and bottom boundaries of the simulation area are all set with second-order materials as absorbing boundaries. Fig.14 (d) is the case where no immersion absorption boundary condition is applied on the surface. There are obvious surface reflection waves and interlayer reflection waves, which is the same as the case where the immersion absorption boundary condition is applied. Fig.14 Compared with (e), the energy of the reflected wave is well suppressed. This model example further proves the feasibility and effectiveness of this embodiment.

[0131] Example 3. Calculation example of three-dimensional uniform velocity model

[0132] The number of grids in the three-dimensional uniform model is 58×59×59, the grid spacing is 10m, the source is the Ricker wavelet with a main frequency of 15Hz excited at the grid points (17, 29, 29), the time sampling interval is 1ms, the model velocity is 2000m / s, and the finite difference method of the fourth order in space and the second order in time is used. The three-dimensional surface function with relatively gentle fluctuations is:

[0133]

[0134] Among them, (x, y, z) represents the coordinates of the spatial position.

[0135] The left, right and bottom boundaries of the simulation area are all set with second-order materials as absorbing boundaries.

[0136] Please refer to Fig.15 , showing a schematic diagram of a three-dimensional undulating surface model in an embodiment of the present specification. Fig.15 (a) is the established three-dimensional surface model; Fig.15 (b) in Fig.15 A partially enlarged surface model diagram in (a). Fig.15 In (b), the black dot “o” represents a ghost node, the black dot “x” represents a pseudo ghost node, the black “*” represents the projection point of the ghost node on the surface, and the black “+” represents the projection point of the pseudo ghost node on the surface.

[0137] Please refer to Fig.16, which shows a schematic diagram of the absorption effect of a three-dimensional immersed absorption boundary in an embodiment of this specification. Fig.16 (a) is a snapshot of the wave field at the 246th time step in the uniform model with no immersed absorbing boundary conditions applied to the three-dimensional undulating surface. Fig.16 The arrow in (a) indicates the reflected wave from the ground surface. The two-dimensional slice of the plane yoz at the same time is as follows: Fig.16 as shown in (c). Fig.16 (b) is a snapshot of the wave field at the 246th time step in the uniform model after the immersion absorbing boundary condition is applied to the surface. The two-dimensional slice on the plane yoz at the same time is as follows: Fig.16 As shown in (d). Fig.16 It can be seen from (a) and (b) as well as (c) and (d) that, in the three-dimensional uniform model, the method in the embodiments of this specification has a good absorption effect on the reflected waves from the ground surface.

[0138] As can be seen from the figure, Fig.16 (a) is a uniform velocity model. No absorbing boundary is applied to the three-dimensional undulating surface. It can be seen that there are obvious reflection waves from the surface. The two-dimensional wave field slice on the yoz plane at the same time is as follows: Fig.16 as shown in (c). Fig.16 (b) shows a snapshot of the wave field after the immersion absorbing boundary condition is applied. The two-dimensional wave field slice on the yoz plane at the same time is shown in Fig.16 As shown in (d), compared with Fig.16 (a) and Fig.16 (b) and Fig.16 (c) and Fig.16 (d), the reflected wave energy is well suppressed, proving that the immersed absorption boundary condition in the embodiment of this specification also has a good absorption effect in the case of three-dimensional undulating surfaces.

[0139] Example 4. Calculation example of three-dimensional layered velocity model

[0140] The number of grids in the three-dimensional model is 58×59×59, the grid spacing is 10m, the source is the Ricker wavelet with a main frequency of 15Hz excited at the grid points (17, 29, 29), the time sampling interval is 1ms, the upper layer velocity of the layered velocity model is 2000m / s, and the lower layer velocity is 3000m / s. The finite difference method of the fourth order in space and the second order in time is used. The three-dimensional undulating surface function is:

[0141]

[0142] Among them, (x, y, z) represents the coordinates of the spatial position.

[0143] The left, right and bottom boundaries of the simulation area are all set with second-order materials as absorbing boundaries.

[0144] Please refer to Fig.17 , which shows a schematic diagram of the absorption effect of a three-dimensional immersed absorption boundary in an embodiment of this specification. Please continue to refer to Fig.17 , Fig.17 (a) is the layered velocity model used in the examples of this specification, in which the upper layer velocity is 2000 m / s and the lower layer velocity is 3000 m / s. Fig.17 (b) is a snapshot of the wave field at the 280th time step in the layered velocity model when no immersion absorbing boundary condition is applied to the surface. Fig.17 In (b), the black arrows indicate the interlayer reflection wave and the surface reflection wave, respectively. The two-dimensional slices on the plane yoz at the same time are as follows: Fig.17 as shown in (d). Fig.17 (c) is a wavelength snapshot at the 280th time step when the surface is immersed in the layered medium and the absorbing boundary condition is applied. The two-dimensional slice of the plane yoz at the same time is as follows: Fig.17 As shown in (e).

[0145] The layered velocity model is Fig.17 As shown in (a). Fig.17 As shown in (b), in the layered velocity model, no absorbing boundary is imposed on the three-dimensional undulating surface. It can be seen that there are obvious reflection waves from the surface and interlayer reflection waves. The two-dimensional wave field slices on the yoz plane at the same time are as follows: Fig.17 as shown in (d). Fig.17 (c) shows a snapshot of the wave field after the immersion absorbing boundary is applied. The two-dimensional wave field slice on the yoz plane at the same time is shown as Fig.17 As shown in (e).

[0146] Compare Fig.17 It can be seen from (b) and (c) as well as (d) and (e) that after using the immersed absorbing boundary conditions in the embodiments of this specification, both the reflected waves in the layered medium and the interlayer reflected waves are well suppressed, further proving the effectiveness and feasibility of the immersed absorbing boundary conditions in the embodiments of this specification.

[0147] Fig.18A flow chart of a forward modeling method for seismic wave fields beneath an undulating surface in an embodiment of the present specification is shown. Although the present specification provides method operation steps or device structures as shown in the following embodiments or drawings, more or fewer operation steps or module units may be included in the method or device based on routine or no creative labor. In the steps or structures where there is no necessary causal relationship logically, the execution order of these steps or the module structure of the device is not limited to the execution order or module structure described in the embodiments of the present specification and shown in the drawings. When the method or module structure is applied to an actual device or terminal product, it can be connected according to the method or module structure shown in the embodiments or drawings for sequential execution or parallel execution (for example, a parallel processor or a multi-threaded processing environment, or even a distributed processing environment).

[0148] Specifically, Fig.18 As shown, a forward modeling method for seismic wave fields under an undulating surface provided in one embodiment of this specification may include the following steps:

[0149] Step S181, obtain the surface model function of the undulating surface; use the immersed boundary method to embed the surface model function into the differential grid model to obtain a finite difference grid model under the undulating surface; calculate the ghost nodes and pseudo-ghost nodes in the finite difference grid model under the undulating surface and the spatial positions of the projection points and fractional nodes corresponding to the ghost nodes and pseudo-ghost nodes.

[0150] Step S182, calculating the wavefield values ​​of the global grid nodes based on the finite difference grid model to obtain the wavefield values ​​of each node in the global integer grid nodes; and determining the wavefield values ​​at the fractional nodes corresponding to the ghost nodes and pseudo ghost nodes by using a parabolic interpolation method.

[0151] Step S183, based on the transmission absorption boundary condition, the wave field values ​​at the ghost nodes and pseudo ghost nodes are updated along the surface normal direction according to the wave field values ​​of the fractional nodes corresponding to the ghost nodes and pseudo ghost nodes to obtain seismic wave forward modeling results.

[0152] The forward modeling method of the seismic wave field under the undulating surface in the embodiments of this specification can be executed by a computer device.

[0153] The surface model function under the undulating surface can be obtained. The surface model function is the function of the undulating surface in the rectangular coordinate system. The surface model function can be embedded into the differential grid model using the immersed boundary method to obtain a finite difference grid model under the undulating surface.

[0154] The immersed boundary method directly embeds the irregular interface into the traditional Cartesian coordinate system. The undulating surface is located on a fractional grid, while the calculation variables in the wave equation are located on an integer grid. The wave equation is discretized in the Cartesian coordinate system, and the physical quantities on the undulating surface can be obtained by interpolation and other methods. Compared with the traditional finite difference method and the curved grid method, the immersed boundary method avoids the step approximation problem while having the advantages of high computational efficiency and better adaptability to the undulating surface. The wave field values ​​of the global grid nodes can be calculated based on the finite difference grid model to obtain the wave field values ​​of each node in the global grid node.

[0155] When the finite difference method is used to calculate the wave field values ​​at the grid points close to the surface, the nodes located above the surface will be used according to the length of the difference operator used. These nodes located above the surface and required for the finite difference calculation of the internal grid nodes are called ghost nodes.

[0156] In one embodiment, assuming that the ground surface can be represented as z=f(x), then the ghost node (x g ,z g ) can be defined in two-dimensional space as follows:

[0157]

[0158] Where Δz is the grid spacing in the z direction, L x and L z represent the difference operator lengths in the x and z directions respectively.

[0159] In some embodiments of the present specification, calculating the spatial positions of ghost nodes in the finite difference grid model under the undulating surface and the projection points and fractional nodes corresponding to the ghost nodes may include: calculating the spatial positions of the ghost nodes in the finite difference grid model after being embedded in the surface; determining the spatial positions of the projection points corresponding to the ghost nodes on the undulating surface based on the spatial positions of the ghost nodes and the surface model function; determining the normal direction inside the surface based on the spatial positions of the ghost nodes and the corresponding projection points; and determining the spatial positions of multiple fractional nodes corresponding to the ghost nodes inward along the normal direction based on a preset interval distance.

[0160] In this embodiment, the positions of all ghost nodes in space are first determined using equation (13), and then the spatial positions of all ghost nodes corresponding to the projection points on the surface are further determined. Each pair of ghost nodes and projection points can determine the internal normal direction ( Figure 5 The solid arrow in the figure is then used to apply boundary conditions along the internal normal direction. Figure 5 As shown, in the internal normal direction, min(Δx,Δz) is selected as the distance to determine the four fractional nodes A with equal spacing. 1 ,A 2 ,A3 ,A 4 , and get the coordinates of the four fractional nodes.

[0161] In some embodiments of the present specification, determining the wavefield value at the fractional node corresponding to the ghost node by using the parabola interpolation method may include: determining whether there is a fractional node located above the ground surface among the fractional nodes corresponding to the ghost node; when it is determined that there is at least one fractional node located above the ground surface among the fractional nodes corresponding to the ghost node, expanding the same number of fractional nodes toward the inside of the ground surface along the normal direction; obtaining the wavefield value at each fractional node among the multiple fractional nodes located below the ground surface corresponding to the ghost node by using parabola interpolation; and calculating the wavefield value of each fractional node among the multiple fractional nodes located below the ground surface corresponding to the ghost node based on the wavefield value at each fractional node among the multiple fractional nodes located below the ground surface corresponding to the ghost node.

[0162] Specifically, it can be determined whether there is a fractional node located above the ground surface among the fractional nodes corresponding to the ghost node. In the case where all fractional nodes corresponding to the ghost node are located below the ground surface, the wavefield value of each fractional node corresponding to the ghost node can be calculated using the parabolic interpolation method. Then, based on the transmission absorption boundary condition, the wavefield value of the ghost node is calculated according to the wavefield value of each fractional node corresponding to the ghost node.

[0163] In the case where there are fractional nodes located above the ground surface among the fractional nodes corresponding to the ghost node, the same number of fractional nodes can be expanded toward the inside of the ground surface along the normal direction. In one embodiment, in the case where two of the four fractional nodes corresponding to the ghost node are located above the ground surface, the two fractional nodes can be expanded toward the inside of the ground surface. The interval between the two expanded fractional nodes is the same as the interval between the original fractional nodes. The parabolic interpolation method can be used to calculate the wavefield values ​​of multiple (for example, four) fractional nodes located below the ground surface. For the first fractional node located above the ground surface, the wavefield value of the first fractional node can be calculated based on the transmission absorption boundary condition according to the wavefield values ​​of the four fractional nodes located below the ground surface. For the second fractional node located above the ground surface, the wavefield value of the second fractional node can be calculated based on the transmission absorption boundary condition according to the wavefield values ​​of the three fractional nodes located below the ground surface and the first fractional node located above the ground surface.

[0164] In some embodiments of the present specification, updating the wavefield value of the ghost node along the surface normal direction according to the wavefield value of the fractional node of the ghost node based on the transmission absorption boundary condition may include: setting an interpolation window or an interpolation cube for each of the multiple fractional nodes below the surface corresponding to the ghost node; determining whether there is a pseudo ghost node in the interpolation window or the interpolation cube corresponding to each of the multiple fractional nodes; and calculating the wavefield value at the ghost node according to the wavefield value of the fractional node corresponding to the ghost node based on the transmission absorption boundary condition when it is determined that there is no pseudo ghost node in the interpolation window or the interpolation cube corresponding to the multiple fractional nodes.

[0165] Specifically, when interpolating the wavefield values ​​of the fractional nodes below the ground surface, due to the difference operator length and the interpolation window size, some grid node wavefield values ​​that are above the ground surface but are not ghost nodes are needed, and such grid nodes are called pseudo ghost nodes. When it is determined that there are no pseudo ghost nodes in the interpolation windows or interpolation cubes corresponding to the multiple fractional nodes below the ground surface, the wavefield values ​​at the ghost nodes can be calculated based on the transmission absorption boundary conditions according to the wavefield values ​​of the fractional nodes corresponding to the ghost nodes.

[0166] In some embodiments of the present specification, calculating the pseudo ghost nodes in the finite difference grid model under the undulating surface and the spatial positions of the projection points and fractional nodes corresponding to the pseudo ghost nodes may include: when it is determined that there are pseudo ghost nodes in the interpolation windows or interpolation cubes corresponding to the multiple fractional nodes, calculating the spatial positions of the pseudo ghost nodes and the projection points and fractional nodes corresponding to the pseudo ghost nodes; correspondingly, determining the wave field values ​​at the fractional nodes corresponding to the pseudo ghost nodes by using a parabolic interpolation method, including: determining whether there are fractional nodes located above the surface among the fractional nodes corresponding to the pseudo ghost nodes; when it is determined that there is at least one fractional node located above the surface among the fractional nodes corresponding to the pseudo ghost nodes, expanding the same number of fractional nodes toward the inside of the surface along the normal direction; obtaining the wave field values ​​at each fractional node among the multiple fractional nodes located below the surface corresponding to the pseudo ghost nodes by using parabolic interpolation; and calculating the wave field values ​​of each fractional node among the multiple fractional nodes located below the surface corresponding to the pseudo ghost nodes based on the wave field values ​​at each fractional node among the multiple fractional nodes located below the surface corresponding to the pseudo ghost nodes.

[0167] Specifically, when it is determined that there are pseudo ghost nodes in the interpolation windows or interpolation cubes corresponding to the multiple score nodes, the spatial positions of the corresponding projection points of all ghost nodes on the surface are determined, and each pair of ghost nodes and projection points can determine the internal normal direction (such as Figure 8 As shown), in the internal normal direction, select min(Δx,Δz) as the distance to determine the four fractional nodes B with equal spacing distribution 1 ,B2 ,B 3 ,B 4 , and obtain the coordinates of the four fractional nodes. After that, the wavefield value at the pseudo ghost node can be calculated based on the transmission absorption boundary condition according to the wavefield value of the fractional node corresponding to the pseudo ghost node. After the update, the wavefield value of the fractional node in the interpolation window or interpolation cube can be calculated according to the wavefield value of the pseudo ghost node, and then the wavefield value of the ghost node corresponding to the fractional node in the interpolation window or interpolation cube can be calculated.

[0168] In some embodiments of the present specification, based on the transmission absorption boundary condition, updating the wavefield value of the ghost node along the surface normal direction according to the wavefield value of the fractional node of the ghost node includes:

[0169] The wave field value at the ghost node is calculated according to the following formula:

[0170]

[0171]

[0172] T 1,2 =s(2-s);

[0173]

[0174] Among them, u represents the wave field value, the superscripts n+1, n, and n-1 represent three adjacent moments respectively, the subscript G represents the ghost node, and the subscript A 1 , A 2 , A 3 , A 4 Represents the four fractional nodes corresponding to the ghost nodes, s = αcΔt / Δx, α is the transmission coefficient, c is the velocity of the seismic wave, Δt is the time interval between adjacent moments, and Δx is the grid spacing along the x-axis.

[0175] In order to improve the absorption effect of other absorption areas except the undulating surface, a grid number L is set outside the calculation area in the embodiment of this specification. c (L c ≥L) buffer, such as Figure 4In some embodiments of the present specification, the finite difference grid model is provided with a buffer area outside the calculation area except for the undulating surface, and the width of the buffer area is greater than or equal to the differential operator length of the finite difference grid model; accordingly, the method further comprises: calculating the first wave field of each node in the grid nodes of the buffer area by using a differential equation; calculating the second wave field of each node in the grid nodes of the buffer area by using a ghost node wave field update formula (i.e., the formula for calculating the wave field value at the ghost node mentioned above); based on the first wave field and the second wave field, performing weighted summation of the first wave field and the second wave field to determine the final wave field value of each node in the buffer area.

[0176] In some embodiments of the present specification, based on the first wave field and the second wave field, performing weighted summation of the first wave field and the second wave field to determine a final wave field value of each node in the buffer area includes:

[0177] The final wavefield value of each node in the buffer area is calculated according to the following formula:

[0178] u(A q )=(1-ω q ) 2 (A q )+ω q u 1 (A q );

[0179]

[0180] Among them, u is the final wave field value, A q is the qth grid node in the buffer area, L c is the total number of grid nodes in the buffer area, u 1 is the first wave field value, u 2 is the second wave field value.

[0181] Each embodiment in this specification is described in a progressive manner, and the same or similar parts between the embodiments can be referred to each other, and each embodiment focuses on the differences from other embodiments. For details, please refer to the description of the above-mentioned related processing related embodiments, and no further description is given here.

[0182] The above is a description of a specific embodiment of the specification. Other embodiments are within the scope of the appended claims. In some cases, the actions or steps recorded in the claims can be performed in an order different from that in the embodiments and still achieve the desired results. In addition, the processes depicted in the drawings do not necessarily require the specific order or continuous order shown to achieve the desired results. In some embodiments, multitasking and parallel processing are also possible or may be advantageous.

[0183] Based on the same inventive concept, a forward modeling device for seismic wavefields under undulating surface is also provided in the embodiments of this specification, as described in the following embodiments. Since the principle of solving the problem by the forward modeling device for seismic wavefields under undulating surface is similar to that of the forward modeling method for seismic wavefields under undulating surface, the implementation of the forward modeling device for seismic wavefields under undulating surface can refer to the implementation of the forward modeling method for seismic wavefields under undulating surface, and the repeated parts will not be repeated. As used below, the term "unit" or "module" can be a combination of software and / or hardware that implements predetermined functions. Although the device described in the following embodiments is preferably implemented in software, the implementation of hardware, or a combination of software and hardware, is also possible and conceived. Fig.19 is a structural block diagram of a forward modeling device for earthquake wave fields under an undulating surface according to an embodiment of the present specification, such as Fig.19 As shown, it includes: an acquisition module 191, a calculation module 192 and an update module 193, and the structure is described below.

[0184] The acquisition module 191 is used to obtain the surface model function of the undulating surface; use the immersed boundary method to embed the surface model function into the differential grid model to obtain a finite difference grid model under the undulating surface; calculate the ghost nodes and pseudo-ghost nodes in the finite difference grid model under the undulating surface and the spatial positions of the projection points and fractional nodes corresponding to the ghost nodes and pseudo-ghost nodes.

[0185] The calculation module 192 is used to calculate the wave field value of the global grid node based on the finite difference grid model to obtain the wave field value of each node in the global integer grid node; and determine the wave field value at the fractional node corresponding to the ghost node and the pseudo ghost node by using the parabolic interpolation method.

[0186] The updating module 193 is used to update the wave field values ​​at the ghost nodes and pseudo ghost nodes along the surface normal direction according to the wave field values ​​of the fractional nodes corresponding to the ghost nodes and pseudo ghost nodes based on the transmission absorption boundary conditions to obtain seismic wave forward modeling results.

[0187] From the above description, it can be seen that the embodiments of this specification achieve the following technical effects: the immersed boundary method is used to directly embed the terrain into the Cartesian coordinate system, so that the relevant physical quantities can be directly calculated in the rectangular coordinate system, which is conducive to the implementation of the finite difference method, while avoiding the problem of step approximation of traditional finite differences; by combining the immersed boundary method and the transmission absorption boundary condition, the undulating surface can be effectively processed, and while having the advantages of good absorption effect and lower calculation cost, it can also better adapt to the undulating changes of the surface.

[0188] This specification also provides a computer device, which can be found in Fig. 20The computer device structure diagram of the forward modeling method of the seismic wave field under the undulating surface provided by the embodiment of this specification is shown, and the computer device may specifically include an input device 201, a processor 202, and a memory 203. Among them, the memory 203 is used to store processor executable instructions. When the processor 202 executes the instructions, the steps of the forward modeling method of the seismic wave field under the undulating surface described in any of the above embodiments are implemented.

[0189] In this embodiment, the input device may specifically be one of the main devices for information exchange between the user and the computer system. The input device may include a keyboard, a mouse, a camera, a scanner, a light pen, a handwriting input board, a voice input device, etc.; the input device is used to input the original data and the program for processing these numbers into the computer. The input device can also obtain and receive data transmitted from other modules, units, and devices. The processor can be implemented in any appropriate manner. For example, the processor can take the form of a computer-readable medium, a logic gate, a switch, an application-specific integrated circuit (Application Specific Integrated Circuit, ASIC), a programmable logic controller, and an embedded microcontroller, etc., such as a microprocessor or a processor and a computer-readable program code (such as software or firmware) that can be executed by the (micro) processor. The memory may specifically be a memory device used to store information in modern information technology. The memory may include multiple levels. In a digital system, anything that can store binary data can be a memory; in an integrated circuit, a circuit with a storage function without a physical form is also called a memory, such as a RAM, a FIFO, etc.; in a system, a storage device with a physical form is also called a memory, such as a memory stick, a TF card, etc.

[0190] In this embodiment, the functions and effects specifically realized by the computer device can be explained in comparison with other embodiments and will not be described in detail here.

[0191] A computer storage medium based on a forward modeling method for seismic wave fields beneath an undulating surface is also provided in an implementation manner of the present specification. The computer storage medium stores computer program instructions, and when the computer program instructions are executed, the steps of the forward modeling method for seismic wave fields beneath an undulating surface described in any of the above embodiments are implemented.

[0192] In this embodiment, the storage medium includes, but is not limited to, a random access memory (RAM), a read-only memory (ROM), a cache, a hard disk (HDD) or a memory card. The memory may be used to store computer program instructions. The network communication unit may be an interface for network connection communication set in accordance with the standard specified by the communication protocol.

[0193] In this embodiment, the functions and effects specifically implemented by the program instructions stored in the computer storage medium can be explained in comparison with other embodiments and will not be repeated here.

[0194] Obviously, those skilled in the art should understand that the modules or steps of the above-mentioned embodiments of this specification can be implemented by a general computing device, they can be concentrated on a single computing device, or distributed on a network composed of multiple computing devices, and optionally, they can be implemented by a program code executable by a computing device, so that they can be stored in a storage device and executed by the computing device, and in some cases, the steps shown or described can be executed in a different order from that here, or they can be made into individual integrated circuit modules, or multiple modules or steps therein can be made into a single integrated circuit module for implementation. In this way, the embodiments of this specification are not limited to any specific combination of hardware and software.

[0195] It should be understood that the above description is for illustration and not for limitation. Many embodiments and many applications beyond the examples provided will be apparent to those skilled in the art upon reading the above description. Therefore, the scope of this specification should not be determined with reference to the above description, but should be determined with reference to the preceding claims and the full scope of equivalents to which these claims belong.

[0196] The above description is only the preferred embodiment of this specification and is not intended to limit this specification. For those skilled in the art, the embodiments of this specification may have various changes and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of this specification shall be included in the protection scope of this specification.

Claims

1. A forward modeling method for seismic wave fields under an undulating surface, characterized in that: include: Obtaining the surface model function of the undulating surface; The surface model function is embedded into the differential grid model by using the immersed boundary method to obtain a finite difference grid model under the undulating surface; Calculating ghost nodes and pseudo ghost nodes in the finite difference grid model under the undulating surface and spatial positions of projection points and fractional nodes corresponding to the ghost nodes and pseudo ghost nodes; The wave field values ​​of the global grid nodes are calculated based on the finite difference grid model to obtain the wave field values ​​of each node in the global integer grid nodes; the wave field values ​​at the fractional nodes corresponding to the ghost nodes and pseudo ghost nodes are determined by using a parabolic interpolation method; Based on the transmission absorption boundary condition, the wave field values ​​at the ghost nodes and pseudo ghost nodes are updated along the surface normal direction according to the wave field values ​​of the fractional nodes corresponding to the ghost nodes and pseudo ghost nodes to obtain the seismic wave forward modeling results.

2. The forward modeling method for seismic wave fields under an undulating surface according to claim 1, characterized in that: Calculating the ghost nodes in the finite difference grid model under the undulating surface and the spatial positions of the projection points and fractional nodes corresponding to the ghost nodes, including: Calculating the spatial positions of ghost nodes in the finite difference grid model embedded in the surface; Determine the spatial position of the projection point corresponding to the ghost node on the undulating surface according to the spatial position of the ghost node and the surface model function; According to the spatial positions of the ghost nodes and the corresponding projection points, the normal direction inside the surface is determined; The spatial positions of a plurality of fractional nodes corresponding to the ghost node are determined inwardly along the normal direction according to a preset interval distance.

3. The forward modeling method for seismic wave fields under an undulating surface according to claim 2, characterized in that: Determining the wave field value at the fractional node corresponding to the ghost node by using a parabolic interpolation method includes: Determine whether there is a score node located above the ground among the score nodes corresponding to the ghost node; When it is determined that at least one of the fractional nodes corresponding to the ghost node is located above the ground surface, the same number of fractional nodes are expanded toward the inside of the ground surface along the normal direction; Using parabolic interpolation to obtain the wave field value at each of the multiple fractional nodes below the ground surface corresponding to the ghost node; Based on the wavefield value at each fractional node among the plurality of fractional nodes below the ground surface corresponding to the ghost node, the wavefield value of each fractional node among at least one fractional node above the ground surface corresponding to the ghost node is calculated.

4. The forward modeling method for seismic wave fields under an undulating surface according to claim 1, characterized in that: Based on the transmission absorption boundary condition, updating the wavefield value of the ghost node along the surface normal direction according to the wavefield value of the fractional node of the ghost node, including: Setting interpolation windows or interpolation cubes respectively for a plurality of fractional nodes below the ground surface corresponding to the ghost nodes; Determine whether there is a pseudo ghost node in the interpolation window or interpolation cube corresponding to each score node in the multiple score nodes; When it is determined that there is no pseudo ghost node in the interpolation window or interpolation cube corresponding to the plurality of fractional nodes, the wavefield value at the ghost node is calculated according to the wavefield values ​​of the fractional nodes corresponding to the ghost node based on the transmission absorption boundary condition.

5. The forward modeling method for seismic wave fields under an undulating surface according to claim 4, characterized in that: Calculating pseudo ghost nodes in the finite difference grid model under the undulating surface and spatial positions of projection points and fractional nodes corresponding to the pseudo ghost nodes, including: When it is determined that a pseudo ghost node exists in the interpolation window or the interpolation cube corresponding to the plurality of score nodes, calculating the pseudo ghost node and the spatial positions of the projection points and score nodes corresponding to the pseudo ghost node; Accordingly, the parabolic interpolation method is used to determine the wave field value at the fractional node corresponding to the pseudo ghost node, including: Determine whether there is a score node located above the ground among the score nodes corresponding to the pseudo ghost node; When it is determined that at least one of the fractional nodes corresponding to the pseudo ghost node is located above the ground surface, the same number of fractional nodes are expanded toward the inside of the ground surface along the normal direction; Using parabolic interpolation to obtain the wave field value at each of the multiple fractional nodes below the ground surface corresponding to the pseudo ghost node; Based on the wavefield value at each fractional node among the plurality of fractional nodes below the ground surface corresponding to the pseudo ghost node, the wavefield value of each fractional node among at least one fractional node above the ground surface corresponding to the pseudo ghost node is calculated.

6. The forward modeling method for seismic wave fields under an undulating surface according to claim 1, characterized in that: Based on the transmission absorption boundary condition, updating the wavefield value of the ghost node along the surface normal direction according to the wavefield value of the fractional node of the ghost node, including: The wave field value at the ghost node is calculated according to the following formula: T 1,2 =s(2-s); Where u represents the wave field value, the superscripts n+1, n, and n-1 represent three adjacent moments respectively, the subscript G represents the ghost node, the subscripts A1, A2, A3, and A4 represent the four fractional nodes corresponding to the ghost node, and s = αcΔt / Δx, where α is the transmission coefficient, c is the velocity of the seismic wave, Δt is the time interval between adjacent moments, and Δx is the grid spacing along the x-axis.

7. The forward modeling method for seismic wave fields under an undulating surface according to claim 1, characterized in that: The finite difference grid model is provided with a buffer area outside the calculation area except for the undulating surface, and the width of the buffer area is greater than or equal to the difference operator length of the finite difference grid model; Accordingly, the method further includes: Calculating the first wave field of each node in the grid nodes of the buffer area by using a difference equation; Calculating the second wave field of each node in the buffer area grid node by using the ghost node wave field update formula; Based on the first wave field and the second wave field, a final wave field value of each node in the buffer area is determined by weighted summing the first wave field and the second wave field.

8. The forward modeling method for seismic wave fields under an undulating surface according to claim 7, characterized in that: Based on the first wave field and the second wave field, weighted summing the first wave field and the second wave field to determine a final wave field value of each node in the buffer area includes: The final wavefield value of each node in the buffer area is calculated according to the following formula: u(A q )=(1-ω q )u2(A q )+ω q u1(A q ); Among them, u is the final wave field value, A q is the qth grid node in the buffer area, L c is the total number of grid nodes in the buffer area, u1 is the first wavefield value, and u2 is the second wavefield value.

9. A forward modeling device for seismic wave fields under an undulating surface, characterized in that: include: An acquisition module, used for acquiring a surface model function of an undulating surface; The surface model function is embedded into the differential grid model by using the immersed boundary method to obtain a finite difference grid model under the undulating surface; ghost nodes and pseudo ghost nodes in the finite difference grid model under the undulating surface and the spatial positions of projection points and fractional nodes corresponding to the ghost nodes and pseudo ghost nodes are calculated; A calculation module is used to calculate the wave field value of the global grid node based on the finite difference grid model to obtain the wave field value of each node in the global integer grid node; and to determine the wave field value at the fractional node corresponding to the ghost node and the pseudo ghost node by using a parabolic interpolation method; The updating module is used to update the wave field values ​​at the ghost nodes and pseudo ghost nodes along the surface normal direction according to the wave field values ​​of the fractional nodes corresponding to the ghost nodes and pseudo ghost nodes based on the transmission absorption boundary conditions to obtain seismic wave forward simulation results.

10. A computer-readable storage medium having computer instructions stored thereon, characterized in that: When the instructions are executed by a processor, the steps of the method according to any one of claims 1 to 8 are implemented.

Images