Reverse time migration method and device based on finite element

[0129]A preferred embodiment of the present application will be described in more detail below with reference to the accompanying drawings. While the preferred embodiments of the present application are shown, it is understood that the present application can be implemented in various forms and should not be embodied herein. Rather, providing these embodiments is to make the present application more thoroughly and complete, and can integrate the scope of the present application to those skilled in the art.

[0130]Seefigure 1.figure 1 A flowchart showing a method of establishing a pre-stack time offset speed field according to an embodiment of the present application. As shown, the method comprises the following steps.

[0131]Step 102, divide the continuous three-dimensional space into finite units that are connected but not coincident, and establish a finite element control equation for spatial discrete.

[0132]The respective units obtained can be a series of a series of units such as triangles, tetrahedrometers, and other components, which can match any complicated areas, so that the imaging accuracy is high whether the rock top interface is also complex structure, its imaging accuracy is high.

[0133]Get the acoustic wave equation under the atmosphere coordinate system in the spatant domain:

[0134]

[0135]Where: P is the sound pressure, V is the wave velocity, T is the time variable, X, Y, Z is spatial variable, and S is a source function;

[0136]Some discrete is used as the following approximate solution:

[0137]

[0138]Among them, Nj(x, y, z) is the interpolation function on the unit node J, which is a function of the spatial position, and the following is simply nj, Dj(t) is the displacement of the unit node J, which is a function of time, here is short-written as DjND is the total number of unit nodes;

[0139]Equation 2 is substituting into formula 1, resulting in a margin R:

[0140]

[0141]According to the weighted balance method, the rational amount of the interpolation function is used, and the balance of weight R is divided into zero, ie

[0142]BambooΩ NiRDω = 0 (i = 1, 2, ..., nd) formula 4

[0143]The formula 3 is introduced into the formula 4, obtained:

[0144]

[0145]Three-in-one branch points for the upper left end of the formula 5 equal sign

[0146]

[0147]Among them, NxNyNzIn the direction of the outer outer outer normal, γ is the outer boundary of Ω;

[0148]Equation 2 and Formula 6 are introduced into Equation 5 to obtain a limited finger element control equation:

[0149]

[0150]among them, DijSecond-order bias.

[0151]Step 104, extract the finite element control equation to obtain a linear equation group of spatial and time-discrete linear equation based on the time integral method.

[0152]In one possible embodiment, step 104 specifically includes:

[0153]Rewritten the finite element control equation as a matrix form:

[0154]

[0155]among them,

[0156]

[0157]Where n is the interpolation function matrix of each unit node, D is the displacement matrix of each unit node, It is a second-order bias matrix of each unit node, ρ is density, Is the Hamiltonian operator, ΩeIndicates that the integral area is a unit;

[0158]In the calculation along the time axis, use the implicit NewMark method to implement time points, in the time zone, the NewMark method uses the following assumptions:

[0159]

[0160]Further, according to the formula 10

[0161]

[0162]among them, It is the first-order deflection matrix of each unit node, ΔT represents time increment, subscript T and T + ΔT indicate the value of the time variable, α and δ are constant coefficients;

[0163]Equation 8 and formula 9 are finished according to Equation 10 and Equation 11, resulting in:

[0164]

[0165]Including:

[0166]

[0167]Solve the formula 12 updated:

[0168]

[0169]among them, Overall stiffness matrix, C0, C1, C2, C3, C4, C5, C6, C7All are constant factors, C is the damping matrix;

[0170]Due to the actual use of random speed boundary conditions, the damping matrix C is not present in the calculation, and the formula 12 is simplified to:

[0171] Border conditions

[0172]Where D0Represents the initial value matrix of each node displacement, Indicates the initial value matrix of the second-order bias of each node.

[0173]Step 106, dividing the three-dimensional space into a plurality of pieces by dividing the surface, and the displacement update value of each unit node in each block is incremented by Q displacement increments.

[0174]figure 2 A depth domain block diagram showing an exemplary embodiment of the present application is shown. The three-dimensional space is divided into multiple pieces, which can be subsequently calculated for subsequent parallel computing.

[0175]In a possible embodiment, the division surface is placed within the unit, i.e., the segmentation surface does not pass the unit, such asimage 3 Indicated. The unit passing through the divided surface can be defined as an interface unit, adjacent to the block shared interface unit, thereby ensuring the displacement continuous condition. When subsequent solving, the block can communicate with adjacent blocks to obtain the solution of the external node on the interface unit.

[0176]In one possible embodiment, the displacement D of the unit internal node i in the first block can be used in the formula.i:

[0177]

[0178]among them, The displacement approximation of the front unit node i for the displacement update The coefficient of incremental increment in the first block, For the mth displacement increment mode of the node i in the first block.

[0179]Set a common ND cell node, the degree of freedom of each node is e, then the total degree of freedom is E * ND, that is, the linear equation group obtained in step 104 has an E * ND unknown pendant, divided into pieces, each The number of degrees of freedom in the block is (E * ND) / b, B is the total number of pieces. According to the present embodiment, by means of Q displacement increment of each unit node in each of the blocks, the block (E * ND) / B fine free degree is coated into Q high order. Freedom, Q can take 10 left and right levels, which greatly reduces the number of unknown designs, which can significantly reduce the needs of computing resources and the need for storage resources, and significantly increase the response speed.

[0180]The displacement increment mode can contain two types of conventional displacement increment modes and adaptive slack displacement increment mode. Conventional displacement incremental mode is used to capture the overall motion deformation trend of the block, the better the number. Conventional displacement incremental modes typically include a flat mode and uniform deformation and rotation mode. The adaptive slack displacement increment mode is used to capture uneven deformation within the block, which is usually calculated from the processor used in each block, and needs to be adjusted according to each iteration.

[0181]Step 108, the linear equation group is converted according to the minimum energy principle, and the q displacement increment is used to represent the displacement update value of the cell node to obtain a multi-stage parallel linear equation group:

[0182]

[0183]among them, SET (i) is a collection of all unit nodes in the i-th block, and set (j) is a collection of all unit nodes in the j tweh. For the first M-bit shift mode of the unit node i in the first block, For the first displacement increment mode of the unit node J in the J2, kIJ For the corresponding overall stiffness, The coefficient of pulling incremental incremental increments for the J2 block. siThe pressure applied to the cell node i is applied to the source of the source. For the displacement update, the displacement approximation of the front cell node i, B is the total number of pieces.

[0184]Equation 15 can write

[0185]

[0186]Where KIJ For the overall stiffness matrix The corresponding elements in the middle.

[0187]The formula 16 is established according to the principle of the minimum energy, the expression form of the total potential energy is

[0188]

[0189]Depending on the principle of the minimum energy, it is obtained to obtain an iterative format such that the value of the functionality (17) is reduced, and the solution of the formula 16 can be obtained.

[0190]Therefore, the formula 18 is substituted into the formula 17, and the formula 17 is deserved according to the functional value condition, resulting in:

[0191]

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

[0193]Thereby, the following multi-stage parallel linear equations:

[0194]

[0195]Step 110: Calling a plurality of processes from the processor to handle data of each block according to the following formula, each corresponding to a block from the processor:

[0196]

[0197]The usual counterclockwise offset is sent to a processor to perform a wave field calculation from a single excitation full-demand domain wave field, namely: a processor calculates a full space wave field excited by a source (gun). Such artillery parallelism is feasible for finite differences in counterclockwise, but for finite element counterclockwise offset, the guns are less efficient, and industrial needs cannot be met.

[0198]In this embodiment, the three-dimensional space is divided into a plurality of pieces, and the linear equation group obtained in step 104 is converted to a multi-stage parallel linear equation group based on the solution domain block, and is parallel to multiple processors. Treat data different from demand domains, significantly improved parallel computational efficiency, greatly shortened processing time.

[0199]Step 112: Call the master processor to receive each process of processing from the processor to solve the multi-stage parallel linear equations.

[0200]Receive from each individual from the processor The calculation result is to solve the formula 20.

[0201]In the above embodiment, the imaging accuracy that is greatly higher than the conventional apical differential numerical value simulation method is obtained by using a finite element numerical simulation method, and the multi-stage parallel processing scheme is obtained by the quaternary blocking scheme, and it will be large. The unknown aggregate is a small number of high-order degrees, which has reached the time consumption with the conventional differential numerical value simulation method, and under the bottleneck of the original computer cluster, significantly reduces the computational resource and storage resource requirements, so that the computer cluster User acceptable imaging processing can be performed, which is very suitable for industrial applications.

[0202]Figure 4 A block diagram of a counterclockwise offset apparatus based on one embodiment of the present application is shown. Such asFigure 4 The apparatus shown includes a finite element control equation establishing unit 402, a time integrating external push unit 404, a blocking unit 406, a multi-stage parallel equation establishing unit 408, from the processor modulation unit 410, and the main processor modulation unit 412.

[0203]The finite element control equation establishing unit 402 is used to divide a continuous three-dimensional space into limited units that communicate but do not coincide, and establish a finite element control equation for spatial discrete.

[0204]The time integral external push unit 404 is used to extract the finite element control equation based on the time integral method to obtain a linear equation group of spatial and time-discrete.

[0205]The block unit 406 is configured to divide the three-dimensional space into a plurality of blocks by dividing the surface, and representing the displacement update value of each unit node in each block by Q disable increment.

[0206]Multi-stage parallel equation establishing unit 408 is used to transform the linear equation in accordance with the minimum energy principle, and substitutively use Q displacement increment to represent the displacement update value of the cell node, obtain multi-stage parallel linear equations:

[0207]

[0208]among them, SET (i) is a collection of all unit nodes in the i-th block, and set (j) is a collection of all unit nodes in the j tweh. For the first M-bit shift mode of the unit node i in the first block, For the first displacement increment mode of the unit node J in the J2, kIJ For the corresponding overall stiffness, The coefficient of pulling incremental incremental increments for the J2 block. siThe pressure applied to the cell node i is applied to the source of the source. For the displacement update, the displacement approximation of the front cell node i, B is the total number of pieces.

[0209]The processor calling unit 410 is used to call a plurality of data from the processor to process each block according to the following formula, each corresponding to a block from the processor:

[0210]

[0211]The master processor calling unit 412 is configured to call the master processor to receive each processor of the processor to solve the multi-stage parallel linear equation.

[0212]In one possible embodiment, the division surface is placed within the unit.

[0213]In one possible embodiment, the finite element control equation establishment unit 402 is specifically used:

[0214]Get the acoustic wave equation under the atmosphere coordinate system in the spatant domain:

[0215]

[0216]Where: P is the sound pressure, V is the wave velocity, T is the time variable, X, Y, Z is spatial variable, and S is a source function;

[0217]Some discrete is used as the following approximate solution:

[0218]

[0219]Among them, Nj(x, y, z) are interpolation functions on the unit node J, here is shortly writtenj, Dj(t) is the displacement of the unit node J, here is short-writtenjND is the total number of unit nodes;

[0220]Equation 2 is substituting into formula 1, resulting in a margin R:

[0221]

[0222]According to the weighted balance method, the rational amount of the interpolation function is used, and the balance of weight R is divided into zero, ie

[0223]BambooΩ NiRDω = 0 (i = 1, 2, ..., nd) formula 4

[0224]The formula 3 is introduced into the formula 4, obtained:

[0225]

[0226]Three-in-one branch points for the upper left end of the formula 5 equal sign

[0227]

[0228]Among them, NxNyNzIn the direction of the outer outer outer normal, γ is the outer boundary of Ω;

[0229]Equation 2 and Formula 6 are introduced into Equation 5 to obtain a limited finger element control equation:

[0230]

[0231]among them, DijSecond-order bias.

[0232]In one possible embodiment, the time integrating external push unit 404 is specifically used:

[0233]Rewritten the finite element control equation as a matrix form:

[0234]

[0235]among them,

[0236]

[0237]Where n is the interpolation function matrix of each unit node, D is the displacement matrix of each unit node, It is a second-order bias matrix of each unit node, ρ is density, Is the Hamiltonian operator, ΩeIndicates that the integral area is a unit;

[0238]In the calculation along the time axis, the Time integral is implemented using the implicit NewMark device, and the following assumptions are used in the time zone:

[0239]

[0240]Further, according to the formula 10

[0241]

[0242]among them, It is the first-order deflection matrix of each unit node, ΔT represents time increment, subscript T and T + ΔT indicate the value of the time variable, α and δ are constant coefficients;

[0243]Equation 8 and formula 9 are finished according to Equation 10 and Equation 11, resulting in:

[0244]

[0245]Including:

[0246]

[0247]Solve the formula 12 updated:

[0248]

[0249]among them, For the overall stiffness matrix, C0, C1, C2, C3, C4, C5, C6, C7All are constant factors, C is the damping matrix;

[0250]Due to the actual use of random speed boundary conditions, the damping matrix C is not present in the calculation, and the formula 12 is simplified to:

[0251] Border conditions

[0252]Where D0Indicates the initial value matrix of the displacement of each node, Indicates the initial value matrix of the second-order bias of each node.

[0253]In one possible embodiment, the displacement update value of each unit node in each block is represented by Q-Disperating increments, including:

[0254]Use the following formula to indicate the displacement D of the unit in the unit in the first blocki:

[0255]

[0256]among them, The coefficient of incremental increment in the first block, The first M-displacement increment mode of node i in the first block

[0257]Application example

[0258]Figure 5 A SigsBee2 speed model is shown. This speed model contains a deposited layer of different periods, which are divided into different sizes of broken blocks by numerous positive tomography and inverse laminates. In addition, in the model, a complicated high-speed rock mass is embedded, and this high-speed rock body causes insufficient lighting, imaging difficult problems. Therefore, the model is often used to test the accuracy of the new offset imaging algorithm.

[0259]Such asFigure 5 The SigsBee2 model shown is 24,384 meters long, 9144 meters deep, and the mesh size is 7.62 × 7.62 meters. In order to verify the precision and efficiency performance of the finite element counterclockwise offset method (Hep-Fe-RTM) having high expansion according to the present application, the comparison analysis of the depth offset (WEM), finite differences in the single-way wavefront Offset (FD-RTM) and finite element counterclockwise offset (Hep-Fe-RTM) imaging results.

[0260]Figure 6 Deep offset imaging results for single-way wavefront programs. As can be seen from the figure, there is an inclination limit based on the depth offset of the one-way wave equation, for a steep formation inclination zone (Figure 6 The area pointed to by the arrow) cannot be accurately imaging, and the cover area of ​​the high-speed salt shield (Figure 6 The area in the middle box) The wave field cannot be reached, so the rock construction area cannot be imaged.

[0261]Figure 7 A finite differential RTM results section of the SigsBee2 speed model are shown. Since the method uses a two-way wave equation description wave field propagation, there is no inclination limit, the complex construction area of ​​the rock (Figure 7 The area in the middle box is imaging larger than the one-way wave equation depth offset. However, since the finite difference method mesh is used to use a rectangular or positive six-sided body equation, it is often unable to achieve the ideal imaging effect when processing speed field transverse irregular variations.Figure 7 A finite differential RTM outcome profile of the SigsBee2 speed model is shown).

[0262]Figure 8 The finite element RTM results profile are shown in accordance with the SIGSBee2 speed model of the present application. Since the finite element method is used to resolve the wave field propagation, the finite element can be used with a series of units such as triangles, tetrahedrous, and other components such as the dispersion domain, and can handle any complex area, so regardless of the rock top interface (region of arrows), Or is the complex structure of the rocks.

[0263]The imaging comparison analysis of the three methods indicated that the highly extended finite element counterclockwise offset method proposed by the present invention has a higher imaging accuracy, which can capture almost all complex construct details.

[0264]The above test is completed on the "dawn" cluster system, and each cannon is 24,384 meters, longitudinal to 9144 meters; the random speed boundary has been extended to a horizontal direction of 27432 meters, and portrait 10668 meters. The rectangular mesh is 30.48 × 15.24 meters, and the 5 × 5 solution domain block scheme is used, and 630000 triangular units are included in each piece, and 26 nodes are called parallel calculations on the "Twilight" cluster system. Three methods SigsBee2 model imaging is shown in Table 1. It can be seen that whether it is a single shot shift, or a full data shift, the counterclockwise offset is greatly increased relative to one-way wave offset consumption. Compared with the finite difference method, it is 46.5 minutes (one cannon, while the operation 26 cannon), according to the method of the present application, according to the method of the present application, the single cannon (26 nodes, one run can be run once) When the conventional finite differential method is counter-off, it is almost comparable to the time consumption, thereby achieving the imaging accuracy of the imaging accuracy is significantly improved relative to the conventional finite differences, and the imaging accuracy is significantly improved, and more complex construct details can be captured. feature.

[0265]Table 1 Three Methods SigsBee2 Model Imaging Match Table (5 × 5)

[0266]

[0267]Figure 9 It is shown that the high-precision speed model of China's eastern complicated small block is 9,000 meters, the longitudinal direction is 5,000 meters, the grid size is 10 × 10 meters, including the high steeper fault of controlled concave boundary, many different tilt ports, and from Up to seven reflective layers. It is a total of 672 guns. Each cannon is a horizontal direction of 9,000 meters, and the vertical direction is 5,000 meters; the random speed boundary extends 1500 meters respectively in the lateral and longitudinal (except horizontal plane), the final offset extrapolation calculation 12,000 meters, longitudinal direction 6500 meters, 8 × 8 solution domain block scheme, 65 nodes are called parallel calculations on "Dawning" cluster system.

[0268]The imaging results of the three methods are respectivelyFigure 10 to 12Indicated. Such asFigure 10 As shown, the one-way wave equation depth offset is influencing the region (elliptical region), or the imaging effect is not satisfactory in the high-and-steepler fault, or the imaging effect is not satisfactory in the low-end discharging area (elliptical region).

[0269]Such asFigure 11 As shown, the finite difference method reversal of the counterclockwise due to the propagation of the two-way wave equation, there is no inclination limit, and the rotary wave imaging can be realized, and the imaging accuracy of the high and steep changing block aggregation area. The improvement in the improvement in the high-steeblex, especially in the upper and steep fairing of the rising plate is most obvious. However, since the finite difference method mesh is used to use a rectangular or positive hexahedral and other rule grids, it cannot achieve the ideal imaging effect (block area) when processing the lower disc complicated small block aggregation area.

[0270]Such asFigure 12 As shown, according to the method of the present application, the finite element method is used to solve the wave field propagation problem, and the discrete method of decrease in the domain unit can handle the distribution of extremely complex speed field distribution such as complex small block aggregation regions, and its imaging results are distributed relative to finite difference distribution. The counter-off offset accuracy is significantly improved (box area).

[0271]The imaging comparison analysis of three methods indicates that the method proposed in this application has a high imaging accuracy, and it can accurately engrave the fracture system in the exact imaging of the complex small block of China, which can accurately engrave the fracture system.

[0272]This model tests three ways to imaging, such as Table 2. It can be seen that whether it is a single shot shift, or a full data shift, the counterclockwise offset is greatly increased relative to one-way wave offset consumption. The method disclosed in this application is almost comparable to the finite difference of the finite difference method, and the advantage of its precise imaging is more significant.

[0273]Table 2 Three Methods Northern Basin model imaging time comparison table

[0274](Hep-Fe-RTM adopts 8 × 8 block scheme)

[0275]

[0276]This application can be a system, method, and / or computer program product. The computer program product can include a computer readable storage medium that carries a computer readable program instruction for implementing the various aspects of the present application.

[0277]The computer readable storage medium can be a tangible device that can hold and store instructions used by the instruction execution device. Computer readable storage media can be, for example, - but is not limited to, - electric storage devices, magnetic storage devices, optical storage devices, electromagnetic storage devices, semiconductor storage devices, or any suitable combination of the above. More specific examples of computer readable storage media (non-exhaustive list) include: Portable Computer Disc, Hard Disk, Random Access Memory (RAM), read-only memory (ROM), removable programmable read-only memory (EPROM Or flash memory), static random access memory (SRAM), portable compression disk read only memory (CD-ROM), digital multi-function disk (DVD), memory stick, floppy disk, mechanical encoding device, such as on it stores instructions The convex structure in the hole or groove, and any suitable combination of the above. The computer readable storage medium used herein is not interpreted as an instantaneous signal itself, such as radio waves or other free propagated electromagnetic waves, electromagnetic waves propagated by waveguides or other transport media (e.g., through the optical pulse of the fiber optic cable), or pass the wire The transmitted electrical signal.

[0278]The various aspects of the present application will be described herein with reference to the flowcharts and / or block diagrams of the method, device (system) and computer program products, in accordance with the method of the present application. It should be understood that each block of the flowchart and / or block diagram and the combination of all boxes in the flowchart and / or block diagram can be implemented by a computer readable program instruction.

[0279]Embodiments of the present application have been described above, and the above description is exemplary, non-exhaustive, and is also not limited to the disclosed embodiments. Many modifications and changes in the art will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the various embodiments described. The selection of the terms used herein is intended to be, in the practice of the principles, practical applications of the various embodiments, or techniques for techniques in the market, or other of the art will appreciate the embodiments disclosed herein.