A finite element based wavefield simulation and imaging method

Abstract

The application discloses a wave field simulation and imaging method based on finite elements. The method comprises the following steps: establishing a finite element control equation; extrapolating the finite element control equation according to a time integration method to obtain a linear equation group; dividing a three-dimensional space into a plurality of blocks, and representing displacement update values of each element node in each block by q displacement increments; transforming the linear equation group, and substituting the displacement update values of the element node represented by the q displacement increments to obtain a multi-level parallel linear equation group; calling a plurality of slave processors to process data of each block respectively; calling a master processor to receive processing results of each slave processor, and solving the multi-level parallel linear equation group according to the processing results of the plurality of slave processors. The technical scheme of the application can obtain accurate imaging finite element reverse-time migration, and can reduce the requirements of computing resources and storage resources, so that a computer cluster can perform user-acceptable imaging processing, and is suitable for industrial application.

Description

Technical Field

[0001] The embodiments of the present invention relate to the field of seismic imaging technology in oil and gas exploration and development, and particularly to a wave field simulation and imaging method based on finite element method. Background Technology

[0002] Currently, all industrially applied pre-stack reverse-time migration methods for seismic earthquakes are based on finite difference methods with different difference schemes. A small number of finite element-based pre-stack reverse-time migration methods suffer from low efficiency due to their massive storage and computational resource requirements, failing to meet the timeliness requirements of industrial applications and remaining at the experimental research stage for two-dimensional models. Therefore, there is an urgent need to research new, efficient reverse-time migration imaging methods based on finite element methods that can be widely applied in industrial settings. Summary of the Invention

[0003] To address the aforementioned technical problems, at least one embodiment of the present invention provides a wavefield simulation and imaging method based on the finite element method, and a highly scalable pre-stack reverse time migration imaging method based on the finite element method. By using a multi-level parallel scheme with block partitioning of the solution domain, the method effectively utilizes a large cluster system to improve parallel efficiency, extending the high accuracy advantage of finite element wavefield simulation from experimental research to industrial applications, and further improving the accuracy of seismic imaging.

[0004] In some optional embodiments, the method includes the following steps:

[0005] The continuous three-dimensional space is divided into a finite number of interconnected but non-overlapping elements, and spatially discrete finite element control equations are established.

[0006] Extrapolating the finite element control equations using the time integration method yields a system of linear equations that are discrete in both space and time.

[0007] The three-dimensional space is divided into multiple blocks, and the displacement update value of each unit node in each block is represented by q displacement increments;

[0008] The linear equations are transformed according to the principle of minimum potential energy, and the displacement update values ​​of the unit nodes, represented by q displacement increments, are substituted into the linear equations to obtain a multi-level parallel linear equation system:

[0009]

[0010] in, Set(I) is the set of all cell nodes in the I-th block, and Set(J) is the set of all cell nodes in the J-th block. This represents the m-th displacement increment pattern of element node i in the I-th block. For the l-th displacement increment pattern of element node j in the J-th block, k ijFor the corresponding overall stiffness, The coefficient for the l-th displacement increment of the J-th block. s i The pressure exerted by the earthquake source on element node i, Let B be the approximate displacement of node i in element i before displacement update, and B be the total number of blocks.

[0011] Multiple slave processors are invoked to process each data block according to the following formula, where each slave processor processes one data block:

[0012]

[0013] The main control processor receives the processing results from each slave processor and solves the multi-level parallel linear equation system based on the received processing results.

[0014] In some alternative embodiments, the three-dimensional space is divided into multiple blocks by dividing surfaces, wherein the dividing surfaces are placed within a cell.

[0015] In some optional embodiments, the step of dividing the continuous three-dimensional space into a finite number of interconnected but non-overlapping elements and establishing spatially discrete finite element governing equations includes:

[0016] Obtain the equation of sound wave in Cartesian coordinates in the space-time domain:

[0017]

[0018] Where P is the sound pressure, v is the wave velocity, t is the time variable, x, y, z are the spatial variables, and s is the source function;

[0019] By employing partial discretization, the following approximate solution is constructed:

[0020]

[0021] Where, N j (x,y,z) is the interpolation function at element node j, hereinafter abbreviated as N. j d j (t) is the displacement of element node j, hereinafter abbreviated as d. j nd is the total number of unit nodes;

[0022] Substituting Formula 2 into Formula 1, we obtain the remainder R:

[0023]

[0024] According to the weighted residual method, the interpolation function of the approximate solution is used as the weighting function, and the weighted integral of the residual R over the solution region Ω is set to zero, i.e.

[0025] ∫ Ω N i RdΩ=0 (i=1,2,…,nd) Formula 4

[0026] Substituting formula 3 into formula 4, we get:

[0027]

[0028] Integrating by parts the three middle terms on the left side of Equation 5, we get...

[0029]

[0030] Where, n x n y n z Let Γ be the direction cosine of the outer normal to the boundary, and let Γ be the outer boundary of Ω.

[0031] Substituting equations 2 and 6 into equation 5, we obtain the spatially discrete finite element governing equations:

[0032]

[0033] Where, d j For d j The second-order partial derivative.

[0034] In some optional embodiments, the step of extrapolating the finite element control equations using the time integration method to obtain a spatially and temporally discrete system of linear equations includes:

[0035] Rewrite the finite element governing equations in matrix form:

[0036]

[0037] in,

[0038]

[0039] Where N is the interpolation function matrix for each element node, and d is the displacement matrix for each element node. ρ is the second-order partial derivative matrix of the displacement of each element node, and ρ is the density. It is the Hamiltonian operator, Ω e The integration region is represented as a single unit;

[0040] In the calculations along the time axis, the implicit NewMark method is used to perform time integration. Within the time region, the NewMark method uses the following assumptions:

[0041]

[0042] Further, according to formula 10, we can obtain...

[0043]

[0044] in, It is the first-order partial derivative matrix of the displacement of each element node, Δt represents the time increment, the subscripts t and t+Δt represent the values ​​of the time variable, and α and δ are constant coefficients;

[0045] By rearranging Formulas 8 and 9 according to Formulas 10 and 11, we get:

[0046]

[0047] Among them are:

[0048]

[0049] After solving Formula 12 and updating, we get:

[0050]

[0051] in, C is the overall stiffness matrix, where c0, c1, c2, c3, c4, c5, c6, and c7 are constant coefficients, and C is the damping matrix.

[0052] Since random velocity boundary conditions are actually used, the damping matrix C does not appear in the calculation, and Equation 12 is simplified to:

[0053] Boundary conditions are

[0054] Where d0 represents the initial value matrix of the displacements of each node. The initial value matrix represents the second-order partial derivatives of the displacements at each node.

[0055] In some optional embodiments, representing the displacement update value of each unit node in each block using q displacement increments includes:

[0056] The displacement d of node i within the element of the I-th block is expressed by the following formula. i :

[0057]

[0058] in, Let be the coefficient of the m-th displacement increment in the I-th block. This represents the m-th displacement increment pattern of node i in the I-th block.

[0059] At least one embodiment of the present invention also provides a wave field simulation and imaging device based on the finite element method, the device comprising:

[0060] The finite element governing equation establishment unit is used to divide a continuous three-dimensional space into a finite number of interconnected but non-overlapping elements, and to establish spatially discrete finite element governing equations.

[0061] The time-integration extrapolation unit is used to extrapolate the finite element control equations according to the time-integration method to obtain a set of linear equations that are discrete in both space and time.

[0062] The block unit is used to divide the three-dimensional space into multiple blocks, and to represent the displacement update value of each unit node in each block with q displacement increments;

[0063] A multi-level parallel equation-building unit is used to transform the linear equation system according to the principle of minimum potential energy, and substitute the displacement update value of the unit node, which is represented by q displacement increments, into the linear equation system to obtain a multi-level parallel linear equation system.

[0064] The processor invocation unit is used to invoke multiple slave processors to process each block of data.

[0065] The main processor calling unit is used to call the main control processor to receive the processing results fed back by each slave processor, and to solve the multi-level parallel linear equation system based on the received processing results.

[0066] At least one embodiment of the present invention also provides an electronic device, characterized in that it comprises:

[0067] At least one processor; and,

[0068] A memory communicatively connected to the at least one processor; wherein,

[0069] The memory stores instructions that can be executed by the at least one processor to enable the at least one processor to perform the finite element-based wave field simulation and imaging method as described above.

[0070] At least one embodiment of the present invention also provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the finite element-based wave field simulation and imaging method as described above.

[0071] At least one embodiment of the present invention also provides a computer program product, including a computer program that, when executed by a processor, implements the steps of the finite element-based wave field simulation and imaging method as described above.

[0072] Compared with the prior art, the wave field simulation and imaging method based on the finite element method provided by the embodiments of the present invention has the following beneficial effects:

[0073] 1) The technical solution of this invention uses the finite element method to solve the wave equation, which has higher accuracy than the conventional finite difference wave field simulation method;

[0074] 2) The technical solution of the present invention adopts a multi-level parallel scheme with block-based solution domain, which achieves a time consumption comparable to that of conventional finite difference reverse time migration, effectively improving the efficiency of the algorithm, thereby realizing the promotion of the advantage of finite element reverse time migration in accurate imaging to industrial applications. Attached Figure Description

[0075] The above and other objects, features and advantages of the present invention will become more apparent from the more detailed description of exemplary embodiments of the invention in conjunction with the accompanying drawings, wherein the same reference numerals generally represent the same components in the exemplary embodiments of the invention.

[0076] Figure 1 A flowchart illustrating the wave field simulation and imaging method based on the finite element method according to Embodiment 1 of the present invention is shown.

[0077] Figure 2 This diagram illustrates the depth domain segmentation of Embodiment 1 of the present invention.

[0078] Figure 3 This diagram illustrates the partitioning of the solution domain using the finite element method according to Embodiment 1 of the present invention.

[0079] Figure 4 A block diagram of the wave field simulation and imaging device based on the finite element method according to Embodiment 2 of the present invention is shown.

[0080] Figure 5 This illustrates the Sigsbee2 velocity model according to Embodiment 3 of the present invention;

[0081] Figure 6 This shows a pre-stack depth migration profile of the one-way wave equation of the Sigsbee2 velocity model according to Embodiment 3 of the present invention.

[0082] Figure 7 This shows a finite difference RTM profile of the Sigsbee2 velocity model according to Embodiment 3 of the present invention.

[0083] Figure 8 This shows a finite element RTM result profile of the Sigsbee2 velocity model according to Embodiment 3 of the present invention;

[0084] Figure 9 This illustrates a speed model of a work area according to Embodiment 3 of the present invention;

[0085] Figure 10 This shows a profile of the pre-stack depth migration result of the single-pass wave equation in a certain work area according to Embodiment 3 of the present invention;

[0086] Figure 11 This shows a finite difference RTM result profile of a certain work area according to Embodiment 3 of the present invention;

[0087] Figure 12 The diagram shows a cross-section of the finite element RTM results for a certain work area according to Embodiment 3 of the present invention. Detailed Implementation

[0088] Preferred embodiments of the invention will now be described in more detail with reference to the accompanying drawings. While preferred embodiments of the invention are shown in the drawings, it should be understood that the invention can be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that the invention will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art.

[0089] Example 1

[0090] Please see Figure 1 . Figure 1 A flowchart illustrating a method for establishing a pre-stack time-migration velocity field according to an embodiment of the present invention is shown. As shown, the method includes the following steps.

[0091] Step 102: Divide the continuous three-dimensional space into a finite number of interconnected but non-overlapping elements, and establish spatially discrete finite element control equations.

[0092] The segmented units can be a series of discrete solution domains such as triangles and tetrahedrons, which can match arbitrarily complex regions. Therefore, the imaging accuracy is very high for both the top interface of the rock mass and complex structures under the rock.

[0093] The equations for sound waves in Cartesian coordinates in the space-time domain are obtained as follows:

[0094]

[0095] Where P is the sound pressure, v is the wave velocity, t is the time variable, x, y, z are the spatial variables, and s is the source function;

[0096] By employing partial discretization, the following approximate solution is constructed:

[0097]

[0098] Where, N j (x,y,z) is the interpolation function at element node j, which is a function of spatial location, hereinafter abbreviated as N. j d j (t) is the displacement of element node j, which is a function of time, and is abbreviated as d below. j nd is the total number of unit nodes;

[0099] Substituting Formula 2 into Formula 1, we obtain the remainder R:

[0100]

[0101] According to the weighted residual method, the interpolation function of the approximate solution is used as the weighting function, and the weighted integral of the residual R over the solution region Ω is set to zero, i.e.

[0102] ∫ Ω N i RdΩ=0 (i=1,2,…,nd) Formula 4

[0103] Substituting equation 3 into equation 4, we get:

[0104]

[0105] Integrating by parts the three middle terms on the left side of Equation 5, we get...

[0106]

[0107] Where, n x n y n z Let Γ be the direction cosine of the outer normal to the boundary, and let Γ be the outer boundary of Ω.

[0108] Substituting equations 2 and 6 into equation 5, we obtain the spatially discrete finite element governing equations:

[0109]

[0110] Where, d j For d j The second-order partial derivative.

[0111] Step 104: Extrapolate the finite element control equations using the time integration method to obtain a set of linear equations that are discrete in both space and time.

[0112] In one possible implementation, step 104 specifically includes:

[0113] Rewrite the finite element governing equations in matrix form:

[0114]

[0115] in,

[0116]

[0117] Where N is the interpolation function matrix for each element node, and d is the displacement matrix for each element node. ρ is the second-order partial derivative matrix of the displacement of each element node, and ρ is the density. It is the Hamiltonian operator, Ωe The integration region is represented as a single unit;

[0118] In the calculations along the time axis, the implicit NewMark method is used to perform time integration. Within the time region, the NewMark method uses the following assumptions:

[0119]

[0120] Further, according to formula 10, we can obtain...

[0121]

[0122] in, It is the first-order partial derivative matrix of the displacement of each element node, Δt represents the time increment, the subscripts t and t+Δt represent the values ​​of the time variable, and α and δ are constant coefficients;

[0123] By rearranging Formulas 8 and 9 according to Formulas 10 and 11, we get:

[0124]

[0125] Among them are:

[0126]

[0127] After solving Formula 12 and updating, we get:

[0128]

[0129] in, This is called the overall stiffness matrix, where c0, c1, c2, c3, c4, c5, c6, and c7 are constant coefficients, and C is the damping matrix.

[0130] Since random velocity boundary conditions are actually used, the damping matrix C does not appear in the calculation, and Equation 12 is simplified to:

[0131] Boundary conditions are

[0132] Where d0 represents the initial value matrix of the displacements of each node. The initial value matrix represents the second-order partial derivatives of the displacements at each node.

[0133] Step 106: Divide the three-dimensional space into multiple blocks by using a segmentation surface, and represent the displacement update value of each unit node in each block with q displacement increments.

[0134] Figure 2This diagram illustrates a depth domain partitioning scheme according to an exemplary embodiment of the present invention. Dividing the three-dimensional space into multiple blocks lays the foundation for subsequent parallel computation.

[0135] In one possible implementation, the dividing surface is placed within the cell, i.e., the dividing surface does not pass through the nodes of the cell, such as... Figure 3 As shown. The elements traversed by the dividing surface can be defined as interface elements. Adjacent blocks share interface elements, thus ensuring displacement continuity. During subsequent solutions, a block can communicate with adjacent blocks to obtain solutions for external nodes on the interface elements.

[0136] In one possible implementation, the displacement d of node i within the element of the I-th block can be expressed by the following formula. i :

[0137]

[0138] in, This represents the approximate displacement of node i in element i before the displacement update. Let be the coefficient of the m-th displacement increment in the I-th block. This represents the m-th displacement increment pattern of node i in the I-th block.

[0139] Suppose there are nd element nodes, and each node has e degrees of freedom in displacement. Then the total degrees of freedom is e*nd, meaning the linear equation system obtained in step 104 has e*nd unknowns to be solved. After dividing the system into blocks, the number of degrees of freedom in each block is (e*nd) / B, where B is the total number of blocks. According to this embodiment, by representing the displacement update value of each element node in each block with q displacement increments, the (e*nd) / B fine degrees of freedom in the block are condensed into q higher-order degrees of freedom. q can be around 10, thereby greatly reducing the number of unknowns to be solved, significantly reducing the demand for computing and storage resources, and significantly improving the response speed.

[0140] Displacement increment patterns can be categorized into two types: conventional displacement increment patterns and adaptive relaxation displacement increment patterns. Conventional displacement increment patterns are used to capture the overall motion deformation trend of the blocks, and the fewer the number of patterns, the better. Conventional displacement increment patterns typically include translational patterns and uniform deformation and rotational patterns. Adaptive relaxation displacement increment patterns are used to capture non-uniform deformation within blocks. They are usually calculated on the slave processors used by each block and need to be adjusted according to the state after each iteration.

[0141] Step 108: Transform the linear equation system according to the principle of minimum potential energy, and substitute the displacement update values ​​of the unit nodes into q displacement increments to obtain a multi-level parallel linear equation system:

[0142]

[0143] in, Set(I) is the set of all cell nodes in the I-th block, and Set(J) is the set of all cell nodes in the J-th block. This represents the m-th displacement increment pattern of element node i in the I-th block. For the l-th displacement increment pattern of element node j in the J-th block, k ij For the corresponding overall stiffness, The coefficient for the l-th displacement increment of the J-th block. s i The pressure exerted by the earthquake source on element node i, Let B be the approximate displacement of node i in element before displacement update, and B be the total number of blocks.

[0144] Formula 15 can be written as

[0145]

[0146] Where k ij Global stiffness matrix The corresponding element in.

[0147] Formula 16 is based on the principle of minimum potential energy, and the total potential energy is expressed as follows:

[0148]

[0149] Based on the principle of minimum potential energy, we obtain an iterative scheme that makes the value of the functional (17) continuously decrease, and then we can obtain the solution of formula 16.

[0150] Therefore, substituting Equation 18 into Equation 17, and taking the partial derivative of Equation 17 according to the functional extremum condition, we obtain:

[0151]

[0152] The meaning of the parameters can be found in the relevant description above.

[0153] This results in the following multi-level parallel linear equation system:

[0154]

[0155] Step 110: Invoke multiple slave processors to process each data block according to the following formula, where each slave processor can process one data block:

[0156]

[0157] Typical reverse time migration (RTM) involves sending the entire solution domain wavefield generated from a single point of origin to a single processor for wavefield calculation; that is, one processor calculates the entire space wavefield generated by a single source (shot). This shot-to-shot parallel approach is feasible for finite difference RTM, but for finite element RTM, the shot-to-shot parallel efficiency is low and cannot meet industrial requirements.

[0158] In this embodiment, the three-dimensional space is divided into multiple blocks. Based on this, the linear equation system obtained in step 104 is transformed into a multi-level parallel linear equation system based on the solution domain blocks. The data of different solution domains are processed in parallel by multiple processors, which significantly improves the parallel computing efficiency and greatly shortens the processing time.

[0159] Step 112: The main control processor is invoked to receive the processing results from each slave processor in order to solve the multi-level parallel linear equation system.

[0160] That is, receiving from each processor The calculation results are then substituted into Formula 20 of the multi-level parallel linear equation system for solution.

[0161] In the above embodiments, by using the finite element numerical simulation method to solve the wave equation, the imaging accuracy can be obtained that is much higher than that of the conventional finite difference numerical simulation method. Furthermore, by dividing the solution domain into blocks to obtain a multi-level parallel processing scheme and condensing the huge number of unknowns into a smaller number of high-order degrees of freedom, the time consumption is comparable to that of the conventional finite difference numerical simulation method. Moreover, under the bottleneck that the original computer cluster could not handle, the requirements for computing and storage resources are significantly reduced, enabling the computer cluster to perform user-acceptable imaging processing, which is very suitable for industrial applications.

[0162] Example 2

[0163] Figure 4 A block diagram of the finite element-based reverse time migration device of this embodiment is shown. Figure 4 As shown, the device includes a finite element control equation establishment unit 402, a time integration extrapolation unit 404, a block unit 406, a multi-level parallel equation establishment unit 408, a slave processor calling unit 410, and a main processor calling unit 412. Wherein:

[0164] The finite element control equation establishment element 402 is used to divide the continuous three-dimensional space into a finite number of interconnected but non-overlapping elements and establish spatially discrete finite element control equations.

[0165] The time integration extrapolation unit 404 is used to extrapolate the finite element control equations according to the time integration method to obtain a set of linear equations that are discrete in both space and time.

[0166] The block unit 406 is used to divide the three-dimensional space into multiple blocks and represent the displacement update value of each unit node in each block with q displacement increments.

[0167] The multi-level parallel equation establishment unit 408 is used to transform the linear equation system according to the minimum potential energy principle, and substitute the displacement update value of the unit node represented by q displacement increments into the linear equation system to obtain the multi-level parallel linear equation system.

[0168] The processor calling unit 410 is used to call multiple slave processors to process each block of data respectively.

[0169] The main processor calling unit 412 is used to call the main control processor to receive the processing results fed back by each slave processor, and solve the multi-level parallel linear equation system based on the received processing results.

[0170] Example 3

[0171] The technical solution of the present invention and its beneficial effects will be further illustrated below with a specific example.

[0172] Figure 5 The Sigsbee2 velocity model is shown. This model incorporates sedimentary layers from different periods, which are fragmented into fault blocks of varying sizes by numerous normal and reverse faults. Additionally, a complex high-velocity rock mass is embedded in the model. This high-velocity rock mass can cause insufficient illumination of the underlying strata, leading to imaging difficulties. Therefore, this model is often used to test the accuracy of new migration imaging algorithms.

[0173] like Figure 5 The Sigsbee2 model shown is 24,384 meters long and 9,144 meters deep, with a mesh size of 7.62 × 7.62 meters. To verify the accuracy and efficiency performance of the highly scalable finite element reverse time migration method (HEP-FE-RTM) according to this invention, the imaging results of one-way wave equation pre-stack depth migration (WEM), finite difference reverse time migration (FD-RTM), and finite element reverse time migration (HEP-FE-RTM) were compared and analyzed.

[0174] Figure 6 This is the pre-stack depth migration imaging result based on the one-way wave equation. As can be seen from the figure, depth migration based on the one-way wave equation is subject to dip angle limitations, especially in areas with steeper strata dip angles (…). Figure 6 The area indicated by the middle arrow cannot be accurately imaged, and the underlying area shielded by high-speed salt bodies ( Figure 6 The wave field cannot reach the area within the Chinese box, therefore the underlying structural region cannot be imaged.

[0175] Figure 7The finite-difference RTM profile of the Sigsbee2 velocity model is shown. Because this method uses a two-way wave equation to describe wave field propagation and has no dip angle limitation, it is suitable for complex tectonic zones beneath rocks (…). Figure 7 Imaging of the area within the Chinese box shows a significant improvement over depth migration based on the one-way wave equation. However, because the finite difference method uses regular meshes such as rectangles or hexahedrons, it often fails to achieve ideal imaging results when dealing with problems involving drastic and irregular lateral changes in the velocity field. Figure 7 (This shows a finite difference RTM profile of the Sigsbee2 velocity model).

[0176] Figure 8 The diagram shows the finite element RTM results of the Sigsbee2 velocity model obtained according to the present invention. Since the wave field propagation problem is solved by the finite element method, and the finite element method can use a series of elements such as triangles and tetrahedrons to discretize the solution domain, it can handle arbitrarily complex regions. Therefore, the imaging accuracy is very high for both the top interface of the rock mass (the area pointed to by the arrow) and the imaging of complex structures under the rock.

[0177] The imaging comparison analysis of the three methods shows that the high scalability finite element reverse time migration method (HEP-FE-RTM) proposed in this invention has high imaging accuracy and can capture almost all complex structural details.

[0178] The above tests were conducted on the "Sugon" cluster system. The analysis domain for each shot was 24,384 meters horizontally and 9,144 meters vertically; while the computation domain after random velocity boundary expansion was 27,432 meters horizontally and 10,668 meters vertically. The rectangular mesh size was 30.48 × 15.24 meters, and a 5 × 5 solution domain partitioning scheme was adopted, with each block containing 630,000 triangular elements. Parallel computation was performed on 26 nodes of the "Sugon" cluster system. The comparison of imaging time for the Sigsbee2 model using the three methods is shown in Table 1. It can be seen that, regardless of whether it is single-shot migration or full-volume migration, the reverse-time migration takes significantly longer than the one-way wave migration. Compared to the finite difference method's reverse time migration, which takes 46.5 minutes per shot (one shot per node, 26 shots processed simultaneously), the method according to the present invention takes only 1.9 minutes per shot (26 nodes per shot, one shot run at a time), which is almost equivalent to the conventional finite difference method's reverse time migration time. Thus, it achieves a significant improvement in imaging accuracy without increasing the time consumption compared to the conventional finite difference method's reverse time migration, and can capture more complex structural details.

[0179] Table 1. Comparison of imaging time for the Sigsbee2 model using three methods (5×5)

[0180]

[0181] Figure 9This paper presents a high-precision velocity model of a complex small fault block in eastern China, with a horizontal domain of 9000 meters and a vertical domain of 5000 meters, and a grid size of 10×10 meters. It includes a large, steep fault with a concave boundary, numerous small faults with different dip angles and orientations, and seven reflection layers from top to bottom. The forward modeling shot dataset comprises 672 shots, each with an analysis domain of 9000 meters horizontally and 5000 meters vertically. The random velocity boundary is extended by 1500 meters horizontally and vertically (excluding the horizontal plane). The final offset extrapolation computation domain is 12000 meters horizontally and 6500 meters vertically. An 8×8 solution domain partitioning scheme is used, and 65 nodes are invoked for parallel computation on the "Shuguang" cluster system.

[0182] The imaging results of the three methods are as follows: Figures 10 to 12 As shown. Figure 10 As shown, the imaging effect of single-path wave equation depth migration is not ideal in the steep fault-affected area of ​​the uplift disk (elliptical region) or in the complex small fault-block cluster area of ​​the downlift disk (square region).

[0183] like Figure 11 As shown, the finite difference method (FDM) reverse time migration, due to its use of two-way wave equations to describe wavefield propagation and the absence of dip limitations, can achieve rotating wave imaging. This improves imaging accuracy in areas influenced by steep, large faults and complex small-block clusters, particularly in the uplifted disk where the improvement in stratigraphic and fault contact relationships is most significant. However, because the FDM uses regular grids such as rectangles or hexahedrons, it cannot achieve ideal imaging results (boxed areas) when processing complex small-block clusters in the downlifted disk.

[0184] like Figure 12 As shown, the method of the present invention uses the finite element method to solve the wave field propagation problem. The solution domain element discretization method is diverse, which can handle extremely complex velocity field distributions such as complex small block clusters. Its imaging results are significantly more accurate than the finite difference reverse time migration (box area).

[0185] The comparative analysis of the imaging methods shows that the method proposed in this invention has high imaging accuracy and can accurately characterize the fracture system in the precise imaging of complex small fault blocks in eastern China, demonstrating a significant advantage.

[0186] The comparison of imaging time for the three methods tested in this model is shown in Table 2. It can be seen that, regardless of whether it is single-shot migration or full-volume migration, reverse-time migration significantly increases the time compared to single-path wave migration. The method disclosed in this invention has an almost equivalent time consumption to the finite-difference reverse-time migration method, while its advantage in accurate imaging is more significant.

[0187] Table 2 Comparison of imaging time for three methods in the Subei Basin model (HEP-FE-RTM uses an 8×8 block scheme)

[0188]

[0189]

[0190] Example 4

[0191] Another embodiment of the present invention relates to an electronic device, comprising: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to perform the finite element-based wave field simulation and imaging methods of the above embodiments.

[0192] The memory and processor are connected via a bus, which can include any number of interconnecting buses and bridges, connecting various circuits of one or more processors and memories. The bus can also connect various other circuits, such as peripheral devices, voltage regulators, and power management circuits, which are well known in the art and will not be described further herein. The bus interface provides an interface between the bus and the transceiver. The transceiver can be a single element or multiple elements, such as multiple receivers and transmitters, providing a unit for communicating with various other devices over a transmission medium. Data processed by the processor is transmitted over the wireless medium via an antenna, which further receives data and transmits it to the processor.

[0193] The processor manages the bus and general processing, and also provides various functions, including timing, peripheral interfaces, voltage regulation, power management, and other control functions. Memory is used to store data used by the processor during operation.

[0194] Example 5

[0195] Another embodiment of the present invention relates to a computer-readable storage medium storing a computer program. When executed by a processor, the computer program implements the finite element-based wave field simulation and imaging method of the above embodiments.

[0196] That is, those skilled in the art will understand that all or part of the steps in the methods of the above embodiments can be implemented by a program instructing related hardware. This program is stored in a storage medium and includes several instructions to cause a device (which may be a microcontroller, chip, etc.) or processor to execute all or part of the steps of the methods described in the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0197] Example 6

[0198] Another embodiment of the present invention relates to a computer program product, including a computer program whose instructions, when executed by a processor, implement the steps of the finite element-based wave field simulation and imaging method of the above embodiments.

[0199] Various aspects of the present invention are described herein with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It should be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer-readable program instructions.

[0200] The various embodiments of the present invention have been described above. These descriptions are exemplary and not exhaustive, nor are they limited to the disclosed embodiments. Many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments. The terminology used herein is chosen to best explain the principles, practical application, or improvement of the technology in the market, or to enable others skilled in the art to understand the embodiments disclosed herein.

Claims

1. A wave field simulation and imaging method based on the finite element method, characterized in that, Includes the following steps: The continuous three-dimensional space is divided into a finite number of interconnected but non-overlapping elements, and spatially discrete finite element control equations are established. Extrapolating the finite element control equations using the time integration method yields a system of linear equations that are discrete in both space and time. The three-dimensional space is divided into multiple blocks, and the displacement update value of each unit node in each block is represented by q displacement increments; The linear equations are transformed according to the principle of minimum potential energy, and the displacement update values ​​of the unit nodes, which are represented by q displacement increments, are substituted into the linear equations to obtain a multi-level parallel linear equation system. Multiple slave processors are invoked to process each block of data separately; The main control processor receives the processing results from each slave processor and solves the multi-level parallel linear equation system based on the received processing results.

2. The wave field simulation and imaging method based on the finite element method according to claim 1, characterized in that, The step of dividing the three-dimensional space into multiple blocks includes: dividing the three-dimensional space into multiple blocks by dividing surfaces.

3. The wave field simulation and imaging method based on the finite element method according to claim 2, characterized in that, The dividing surface is placed within the cell.

4. The wave field simulation and imaging method based on the finite element method according to claim 1, characterized in that, The process of representing the displacement update value of each unit node in each block using q displacement increments includes expressing the displacement d of node i within the unit of the I-th block using the following formula. i : in, Let be the coefficient of the m-th displacement increment in the I-th block. This represents the m-th displacement increment pattern of node i in the I-th block.

5. The wave field simulation and imaging method based on the finite element method according to claim 4, characterized in that, The multi-level parallel linear equation system is as follows: in, Set(I) is the set of all cell nodes in the I-th block, and Set(J) is the set of all cell nodes in the J-th block. This represents the m-th displacement increment pattern of element node i in the I-th block. For the l-th displacement increment pattern of element node j in the J-th block, k ij For the corresponding overall stiffness, The coefficient for the l-th displacement increment of the J-th block. s i The pressure exerted by the earthquake source on element node i, Let B be the approximate displacement of node i in the element before the displacement update, and let B be the total number of blocks.

6. The wave field simulation and imaging method based on the finite element method according to claim 5, characterized in that, The data for each block is as follows:

7. The wave field simulation and imaging method based on the finite element method according to claim 1, characterized in that, Each processor processes one block of data.

8. A wave field simulation and imaging device based on the finite element method, characterized in that, include: The finite element governing equation establishment unit is used to divide a continuous three-dimensional space into a finite number of interconnected but non-overlapping elements, and to establish spatially discrete finite element governing equations. The time-integration extrapolation unit is used to extrapolate the finite element control equations according to the time-integration method to obtain a set of linear equations that are discrete in both space and time. The block unit is used to divide the three-dimensional space into multiple blocks, and to represent the displacement update value of each unit node in each block with q displacement increments; A multi-level parallel equation-building unit is used to transform the linear equation system according to the principle of minimum potential energy, and substitute the displacement update value of the unit node, which is represented by q displacement increments, into the linear equation system to obtain a multi-level parallel linear equation system. The processor invocation unit is used to invoke multiple slave processors to process each block of data. The main processor calling unit is used to call the main control processor to receive the processing results fed back by each slave processor, and to solve the multi-level parallel linear equation system based on the received processing results.

9. An electronic device, characterized in that, include: At least one processor; as well as, A memory communicatively connected to the at least one processor; wherein, The memory stores instructions that can be executed by the at least one processor to enable the at least one processor to perform the finite element-based wave field simulation and imaging method as described in any one of claims 1 to 7.

10. A computer-readable storage medium storing a computer program, characterized in that, When the computer program is executed by the processor, it implements the wave field simulation and imaging method based on the finite element method as described in any one of claims 1 to 7.

Images