CN112731533A - Reverse time migration method and device based on finite element - Google Patents
Reverse time migration method and device based on finite element Download PDFInfo
- Publication number
- CN112731533A CN112731533A CN201910973951.XA CN201910973951A CN112731533A CN 112731533 A CN112731533 A CN 112731533A CN 201910973951 A CN201910973951 A CN 201910973951A CN 112731533 A CN112731533 A CN 112731533A
- Authority
- CN
- China
- Prior art keywords
- displacement
- unit
- node
- equation
- block
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 87
- 238000013508 migration Methods 0.000 title claims abstract description 49
- 230000005012 migration Effects 0.000 title claims abstract description 48
- 230000002441 reversible effect Effects 0.000 title claims abstract description 31
- 238000006073 displacement reaction Methods 0.000 claims abstract description 104
- 230000010354 integration Effects 0.000 claims abstract description 18
- 238000012545 processing Methods 0.000 claims abstract description 11
- 230000008569 process Effects 0.000 claims abstract description 9
- 239000011159 matrix material Substances 0.000 claims description 61
- 230000036961 partial effect Effects 0.000 claims description 31
- 238000004364 calculation method Methods 0.000 claims description 21
- 238000005192 partition Methods 0.000 claims description 18
- 238000013016 damping Methods 0.000 claims description 12
- 238000013213 extrapolation Methods 0.000 claims description 8
- 238000005381 potential energy Methods 0.000 claims description 7
- 239000006185 dispersion Substances 0.000 claims description 6
- 230000001131 transforming effect Effects 0.000 claims description 5
- 238000003384 imaging method Methods 0.000 abstract description 34
- 238000003860 storage Methods 0.000 abstract description 8
- 230000002829 reductive effect Effects 0.000 abstract description 4
- 230000006870 function Effects 0.000 description 22
- 238000010586 diagram Methods 0.000 description 10
- 239000011435 rock Substances 0.000 description 7
- 238000004458 analytical method Methods 0.000 description 4
- 230000008901 benefit Effects 0.000 description 4
- 238000004088 simulation Methods 0.000 description 4
- 230000002776 aggregation Effects 0.000 description 3
- 238000004220 aggregation Methods 0.000 description 3
- 238000004590 computer program Methods 0.000 description 3
- YQGOJNYOYNNSMM-UHFFFAOYSA-N eosin Chemical compound [Na+].OC(=O)C1=CC=CC=C1C1=C2C=C(Br)C(=O)C(Br)=C2OC2=C(Br)C(O)=C(Br)C=C21 YQGOJNYOYNNSMM-UHFFFAOYSA-N 0.000 description 3
- 230000001788 irregular Effects 0.000 description 3
- 238000000638 solvent extraction Methods 0.000 description 3
- 238000012360 testing method Methods 0.000 description 3
- 230000003044 adaptive effect Effects 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 2
- 230000000903 blocking effect Effects 0.000 description 2
- 238000004422 calculation algorithm Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000005284 excitation Effects 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 230000001902 propagating effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000001174 ascending effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 239000000835 fiber Substances 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000009776 industrial production Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 150000003839 salts Chemical class 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/51—Migration
- G01V2210/512—Pre-stack
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The application discloses a finite element-based reverse time migration method and a device. The method comprises the following steps: establishing a finite element control equation; extrapolating a finite element control equation according to a time integration method to obtain a linear equation set; dividing the three-dimensional space into a plurality of blocks, and expressing the displacement update value of each unit node in each block by q displacement increments; converting the linear equation set, and substituting the linear equation set into a displacement update value which uses q displacement increments to represent unit nodes to obtain a multi-level parallel linear equation set; calling a plurality of slave processors to process data of each block respectively; and calling the master control processor to receive the processing result of each slave processor. By applying the method and the device, the finite element reverse time migration of accurate imaging can be obtained, and under the bottleneck that the original computer cluster cannot bear, the requirements of computing resources and storage resources are obviously reduced, so that the computer cluster can carry out imaging processing acceptable to users, and the method and the device are very suitable for industrial application and popularization.
Description
Technical Field
The invention belongs to the field of seismic imaging in oil and gas exploration and development, and particularly relates to a reverse time migration method based on finite elements and a reverse time migration device based on the finite elements.
Background
The existing industrially applied seismic prestack reverse time migration methods are all finite difference methods based on different difference formats. The finite difference reverse time migration imaging method is simple in calculation and high in efficiency, and meets the timeliness requirement during industrial production. However, this finite difference reverse time migration imaging method is not accurate when facing a very irregular earth surface or an underground rock mass.
And a few researchers research the finite element-based prestack reverse time migration method, the boundary adaptability is strong, and accurate imaging can be carried out on extremely irregular earth surfaces or underground rock masses. However, the existing finite element-based prestack reverse time migration method has huge calculation amount and low parallel efficiency, cannot meet the requirement of industrial application on timeliness, and stays at the experimental research stage aiming at a two-dimensional model.
Disclosure of Invention
In view of this, the present application provides a method for efficient reverse time migration imaging based on finite elements, which can be generalized in industrial applications, and a corresponding apparatus.
According to an aspect of the application, a finite element based reverse time migration method is proposed, the method comprising:
dividing a continuous three-dimensional space into a limited number of units which are mutually communicated but do not coincide, and establishing a finite element control equation with discrete space;
extrapolating a finite element control equation according to a time integration method to obtain a linear equation set with discrete space and time;
dividing the three-dimensional space into a plurality of blocks through a dividing surface, and expressing the displacement update value of each unit node in each block by q displacement increments;
and transforming the linear equation set according to the minimum potential energy principle, and substituting the displacement update values which use q displacement increments to represent unit nodes to obtain a multi-stage parallel linear equation set:
wherein,set (I) is the set of all unit nodes in the I-th block, set (J) is the set of all unit nodes in the J-th block,for the mth shift increment pattern of the unit node I in the ith block,for the l-th shift increment pattern, k, of cell node J in the J-th blockijFor the sake of a corresponding overall stiffness,the coefficients for the jth shift increment of the jth block,sifor the seismic source to exert pressure on cell node i,the displacement approximate value of the unit node i before the displacement updating is performed, and B is the total number of the blocks;
calling a plurality of slave processors to process data of each block according to the following formula, wherein each slave processor corresponds to one block:
and calling the master processor to receive the processing result of each slave processor so as to solve the multi-stage parallel linear equation system.
In one possible embodiment, the dividing surfaces are placed in cells.
In one possible embodiment, the dividing the continuous three-dimensional space into a finite number of units which are connected with each other but are not overlapped, and establishing a spatially discrete finite element control equation includes:
obtaining a sound wave equation under a rectangular coordinate system in a space-time domain:
wherein P is sound pressure, v is wave velocity, t is time variable, x, y, z are space variables, and s is seismic source function;
with partial dispersion, the following approximate solution is constructed:
wherein N isj(x, y, z) is an interpolation function at cell node j, abbreviated as N hereinafterj,dj(t) is the displacement of the cell node j, abbreviated as d hereinafterjNd is the total number of unit nodes;
substituting equation 2 into equation 1 yields the margin R:
according to the weighted residue method, an interpolation function of approximate solution is used as a weight function, and the weighted integral of the residue R on the solving area omega is zero, namely
∫ΩNiRd Ω ═ 0(i ═ 1,2, …, nd) formula 4
Substituting formula 3 into formula 4 yields:
dividing and integrating the three terms in the middle of the left end of the equal sign of the formula 5 to obtain
Wherein n isx、ny、nzIs the direction cosine of the boundary outer normal, Γ is the outer boundary of Ω;
substituting the formulas 2 and 6 into formula 5 to obtain a finite element control equation with discrete space:
In one possible embodiment, the extrapolating the finite element control equation according to the time integration method to obtain the linear equation system with the discrete space and time includes:
rewriting the finite element control equation into a matrix form:
wherein,
wherein N is an interpolation function matrix of each unit node, d is a displacement matrix of each unit node,is the displacement second order partial derivative matrix of each unit node, p is the density,is the Hamiltonian, ΩeIndicating that the integration area is a unit;
in the calculation along the time axis, the implicit NewMark method is used to achieve time integration, and the NewMark method uses the following assumptions in the time region:
further, it can be obtained from equation 10
Wherein,is a displacement first-order partial derivative matrix of each unit node, delta t represents time increment, subscript t and t + delta t represent values of time variables, and alpha and delta are constant coefficients;
by collating equations 8 and 9 according to equations 10 and 11, we obtain:
among them are:
after solving equation 12, update results in:
wherein,as a global stiffness matrix, c0、c1、c2、c3、c4、c5、c6、c7All are constant coefficients, and C is a damping matrix;
since the random velocity boundary condition is actually used, the damping matrix C does not appear in the calculation, and equation 12 is simplified as:
Wherein d is0A matrix of initial values representing the displacement of each node,and the initial value matrix represents the second-order partial derivatives of the displacement of each node.
In a possible implementation, the representing the displacement update value of each unit node in each partition by q displacement increments includes:
the displacement d of the node I in the cell in the I-th block is expressed byi:
Wherein,is the coefficient of the mth shift increment in the ith block,is the mth shift increment pattern for node I in the ith partition.
According to another aspect of the present application, a finite element based reverse time migration apparatus is provided, the apparatus comprising:
the finite element control equation establishing unit is used for dividing a continuous three-dimensional space into a limited number of units which are mutually communicated but do not coincide, and establishing a finite element control equation with discrete space;
the time integral extrapolation unit is used for extrapolating a finite element control equation according to a time integral method to obtain a linear equation set with discrete space and time;
the block unit is used for dividing the three-dimensional space into a plurality of blocks through a dividing surface, and expressing the displacement update value of each unit node in each block by q displacement increments;
the multilevel parallel equation establishing unit is used for transforming the linear equation set according to the minimum potential energy principle and substituting the displacement update value which uses q displacement increments to represent unit nodes to obtain a multilevel parallel linear equation set:
wherein,set (I) is the set of all unit nodes in the I-th block, set (J) is the set of all unit nodes in the J-th block,for the mth shift increment pattern of the unit node I in the ith block,for the l-th shift increment pattern, k, of cell node J in the J-th blockijFor the sake of a corresponding overall stiffness,the coefficients for the jth shift increment of the jth block,sifor the seismic source to exert pressure on cell node i,the displacement approximate value of the unit node i before the displacement updating is performed, and B is the total number of the blocks;
a slave processor calling unit, configured to call a plurality of slave processors to process data of each partition according to the following formula, where each slave processor corresponds to one partition:
and the master control processor calling unit is used for calling the master control processor to receive the processing result of each slave processor so as to solve the multi-stage parallel linear equation set.
In one possible embodiment, the dividing surfaces are placed in cells.
In a possible embodiment, the finite element control equation establishing unit is specifically configured to:
obtaining a sound wave equation under a rectangular coordinate system in a space-time domain:
wherein P is sound pressure, v is wave velocity, t is time variable, x, y, z are space variables, and s is seismic source function;
with partial dispersion, the following approximate solution is constructed:
wherein N isj(x, y, z) is an interpolation function at cell node j, abbreviated as N hereinafterj,dj(t) is the displacement of the cell node j, abbreviated as d hereinafterjNd is the total number of unit nodes;
substituting equation 2 into equation 1 yields the margin R:
according to the weighted residue method, an interpolation function of approximate solution is used as a weight function, and the weighted integral of the residue R on the solving area omega is zero, namely
∫ΩNiRd Ω ═ 0(i ═ 1,2, …, nd) formula 4
Substituting formula 3 into formula 4 yields:
dividing and integrating the three terms in the middle of the left end of the equal sign of the formula 5 to obtain
Wherein n isx、ny、nzIs the direction cosine of the boundary outer normal, Γ is the outer boundary of Ω;
substituting the formulas 2 and 6 into formula 5 to obtain a finite element control equation with discrete space:
In a possible implementation, the time-integral extrapolation unit is specifically configured to:
rewriting the finite element control equation into a matrix form:
wherein,
wherein N is an interpolation function matrix of each unit node, d is a displacement matrix of each unit node,is the displacement second order partial derivative matrix of each unit node, p is the density,is the Hamiltonian, ΩeIndicating that the integration area is a unit;
in the calculation along the time axis, time integration is achieved using an implicit NewMark device, which uses the following assumptions in the time domain:
further, it can be obtained from equation 10
Wherein,is a displacement first-order partial derivative matrix of each unit node, delta t represents time increment, subscript t and t + delta t represent values of time variables, and alpha and delta are constant coefficients;
by collating equations 8 and 9 according to equations 10 and 11, we obtain:
among them are:
after solving equation 12, update results in:
wherein,as a global stiffness matrix, c0、c1、c2、c3、c4、c5、c6、c7All are constant coefficients, and C is a damping matrix;
since the random velocity boundary condition is actually used, the damping matrix C does not appear in the calculation, and equation 12 is simplified as:
Wherein d is0A matrix of initial values representing the displacement of each node,and the initial value matrix represents the second-order partial derivatives of the displacement of each node.
In a possible implementation, the representing the displacement update value of each unit node in each partition by q displacement increments includes:
the displacement d of the node I in the cell in the I-th block is expressed byi:
Wherein,is the coefficient of the mth shift increment in the ith block,is the mth shift increment pattern for node I in the ith partition.
According to the technical scheme of the application, the wave equation is solved by adopting a finite element method, and the precision is higher compared with the conventional finite difference wave field simulation method; and a solution domain block multistage parallel scheme is adopted, so that the time consumption equivalent to the time-reversal migration of a conventional finite difference method is achieved, the algorithm efficiency is effectively improved, and the popularization of the advantage that the finite element time-reversal migration can accurately image to industrial application is realized.
Drawings
The foregoing and other objects, features and advantages of the application will be apparent from the following more particular descriptions of exemplary embodiments of the application, as illustrated in the accompanying drawings wherein like reference numbers generally represent like parts throughout the exemplary embodiments of the application.
FIG. 1 illustrates a flow diagram of a finite element-based reverse time migration method according to an embodiment of the present application.
Fig. 2 shows a depth domain partitioning diagram according to an exemplary embodiment of the present application.
Fig. 3 illustrates a finite element solution domain partitioning diagram according to an exemplary embodiment of the present application.
FIG. 4 illustrates a block diagram of a finite element-based reverse time migration apparatus according to one embodiment of the present application.
Figure 5 shows the sigbee 2 velocity model.
Figure 6 shows a one-way wave equation prestack depth migration effort profile for the sigbee 2 velocity model.
Figure 7 shows a finite difference RTM outcome profile of the Sigsbee2 velocity model.
Figure 8 shows a finite element RTM outcome profile of the Sigsbee2 velocity model obtained in accordance with the present application.
Fig. 9 shows a velocity model for a work area.
FIG. 10 shows a one-way wave equation prestack depth migration effort profile for this work zone.
Fig. 11 shows a finite difference RTM outcome profile for the work area.
FIG. 12 shows a finite element RTM outcome profile for the work zone obtained in accordance with the present application.
Detailed Description
Preferred embodiments 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 in the drawings, it should be understood that the present application may be embodied in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.
Please refer to fig. 1. FIG. 1 shows a flow diagram of a method of establishing a pre-stack time-migration velocity field according to one embodiment of the present application. As shown, the method includes the following steps.
Each unit obtained by segmentation can be a series of unit discrete solution domains such as triangles and tetrahedrons and can be matched with any complex region, so that the imaging precision is high for the rock top interface and complex underground structure.
Obtaining a sound wave equation under a rectangular coordinate system in a space-time domain:
wherein P is sound pressure, v is wave velocity, t is time variable, x, y, z are space variables, and s is seismic source function;
with partial dispersion, the following approximate solution is constructed:
wherein N isj(x, y, z) is an interpolation function on cell node j, which is a function of spatial position, abbreviated as N belowj,dj(t) is the displacement of the cell node j as a function of time, hereinafter abbreviated as djNd is the total number of unit nodes;
substituting equation 2 into equation 1 yields the margin R:
according to the weighted residue method, an interpolation function of approximate solution is used as a weight function, and the weighted integral of the residue R on the solving area omega is zero, namely
∫ΩNiRd Ω ═ 0(i ═ 1,2, …, nd) formula 4
Substituting formula 3 into formula 4 yields:
dividing and integrating the three terms in the middle of the left end of the equal sign of the formula 5 to obtain
Wherein n isx、ny、nzIs the direction cosine of the boundary outer normal, Γ is the outer boundary of Ω;
substituting the formulas 2 and 6 into formula 5 to obtain a finite element control equation with discrete space:
And 104, extrapolating a finite element control equation according to a time integration method to obtain a linear equation set with both discrete space and time.
In a possible implementation manner, step 104 specifically includes:
rewriting the finite element control equation into a matrix form:
wherein,
wherein N is an interpolation function matrix of each unit node, d is a displacement matrix of each unit node,is the displacement second order partial derivative matrix of each unit node, p is the density,is the Hamiltonian, ΩeIndicating that the integration area is a unit;
in the calculation along the time axis, the implicit NewMark method is used to achieve time integration, and the NewMark method uses the following assumptions in the time region:
further, it can be obtained from equation 10
Wherein,is a displacement first-order partial derivative matrix of each unit node, delta t represents time increment, subscript t and t + delta t represent values of time variables, and alpha and delta are constant coefficients;
by collating equations 8 and 9 according to equations 10 and 11, we obtain:
among them are:
after solving equation 12, update results in:
wherein,referred to as the global stiffness matrix, c0、c1、c2、c3、c4、c5、c6、c7All are constant coefficients, and C is a damping matrix;
since the random velocity boundary condition is actually used, the damping matrix C does not appear in the calculation, and equation 12 is simplified as:
Wherein d is0A matrix of initial values representing the displacement of the respective nodes,and the initial value matrix represents the second-order partial derivatives of the displacement of each node.
And 106, dividing the three-dimensional space into a plurality of blocks through the dividing surfaces, and expressing the displacement update values of all unit nodes in each block by q displacement increments.
Fig. 2 shows a depth domain partitioning diagram according to an exemplary embodiment of the present application. The three-dimensional space is divided into a plurality of blocks, and a foundation can be laid for subsequent parallel computing.
In one possible embodiment, the splitting planes are placed inside the cell, i.e. the splitting planes do not pass through the nodes of the cell, as shown in fig. 3. The cell through which the partition plane passes can be defined as an interface cell, and adjacent partitions share the interface cell, so that a displacement continuous condition is guaranteed. Upon subsequent solutions, the tiles may communicate with adjacent tiles to obtain solutions for external nodes on the interface unit.
In one possible implementation, the displacement d of the node I in the cell in the I-th partition can be expressed by the following formulai:
Wherein,updating the displacement approximation of the unit node i before the displacement,is the coefficient of the mth shift increment in the ith block,is the mth shift increment pattern for node I in the ith partition.
And (2) setting a total number of nd unit nodes, wherein the displacement freedom of each node is e, and the total freedom is e nd, namely the linear equation set obtained in the step 104 has a total number of e nd unknowns to be solved, after the blocks are divided, the number of the degrees of freedom in each block is (e nd)/B, and B is the total number of the blocks. According to the embodiment, the displacement update values of the unit nodes in each block are represented by q displacement increments, the (e x nd)/B fine degrees of freedom in each block are aggregated into q high-order degrees of freedom, and q can take about 10 orders of magnitude, so that the number of unknown numbers to be solved is greatly reduced, the requirements on computing resources and storage resources can be remarkably reduced, and the response speed is remarkably improved.
The displacement increment mode can comprise two types of a conventional displacement increment mode and an adaptive relaxation displacement increment mode. The conventional displacement increment mode is used for capturing the deformation trend of the whole motion of the blocks, and the smaller the number, the better the number. The conventional incremental displacement modes generally include a translational mode and a uniform deformation and rotation mode. The adaptive relaxed-displacement delta mode is used to capture the non-uniform deformation in the blocks, and is usually placed in a slave processor for each block to perform calculation, and needs to be adjusted according to the state after each iteration.
108, transforming the linear equation set according to the minimum potential energy principle, and substituting the displacement update values which use q displacement increments to represent unit nodes to obtain a multi-level parallel linear equation set:
wherein,set (I) is the set of all unit nodes in the I-th block, set (J) is the set of all unit nodes in the J-th block,for the mth shift increment pattern of the unit node I in the ith block,for the l-th shift increment pattern, k, of cell node J in the J-th blockijFor the sake of a corresponding overall stiffness,the coefficients for the jth shift increment of the jth block,sifor the seismic source to exert pressure on cell node i,and B is the displacement approximate value of the unit node i before the displacement update, and B is the total number of the blocks.
Equation 15 can be written as
Equation 16 is based on the principle of minimum potential energy, and the total potential energy is expressed in the form of
From the minimum potential principle, an iterative format is obtained such that the value of the functional (17) is constantly decreasing, and a solution of equation 16 is obtained.
Therefore, substituting equation 18 into equation 17 and solving the partial derivative of equation 17 according to the functional extremum condition yields:
wherein the meaning of the parameters can be found in the above-mentioned relevant description.
Thereby obtaining the following multi-stage parallel linear equation system:
the general inverse time migration is to send the full solution domain wave field from one point excitation to a processor for wave field calculation, that is: a processor calculates the full spatial wavefield of a seismic source (shot) excitation. The inter-shot parallel mode is feasible for finite difference method inverse time migration, but for finite element inverse time migration, the inter-shot parallel efficiency is low, and the industrial requirements cannot be met.
In this embodiment, the three-dimensional space is divided into a plurality of blocks, the linear equation set obtained in step 104 is converted into a multi-stage parallel linear equation set based on solution domain blocks, and a plurality of processors process data of different solution domains in parallel, so that the parallel computing efficiency is significantly improved, and the processing time is greatly shortened.
And step 112, calling the master control processor to receive the processing results of the slave processors so as to solve the multi-stage parallel linear equation system.
I.e. from the respective slave processorAnd substituting the calculation result into the formula 20 to solve.
In the embodiment, the wave equation is solved by adopting a finite element numerical simulation method, so that the imaging precision which is greatly higher than that of the conventional finite difference numerical simulation method can be obtained; and the solution domain is partitioned to obtain a multi-level parallel processing scheme, and a huge unknown number is condensed into a small number of high-order degrees of freedom, so that the time consumption equivalent to that of the conventional finite difference numerical simulation method is achieved, and the requirements of computing resources and storage resources are obviously reduced under the bottleneck that the original computer cluster cannot bear, so that the computer cluster can perform imaging processing acceptable to users, and the method is very suitable for industrial application.
FIG. 4 illustrates a block diagram of a finite element-based reverse time migration apparatus according to one embodiment of the present application. The apparatus shown in FIG. 4 includes a finite element control equation building unit 402, a time integral extrapolation unit 404, a block unit 406, a multi-stage parallel equation building unit 408, a slave processor invocation unit 410, and a master processor invocation unit 412.
The finite element control equation establishing unit 402 is configured to divide a continuous three-dimensional space into a finite number of units which are connected with each other but do not overlap with each other, and establish a spatially discrete finite element control equation.
The time integral extrapolation unit 404 is configured to extrapolate the finite element control equation according to a time integral method to obtain a linear equation set with discrete space and time.
The block unit 406 is configured to divide the three-dimensional space into a plurality of blocks by partition planes, and represent displacement update values of respective unit nodes in each block by q displacement increments.
The multi-stage parallel equation establishing unit 408 is configured to transform the linear equation set according to the minimum potential principle, and substitute the displacement update value representing a unit node by q displacement increments to obtain a multi-stage parallel linear equation set:
wherein,set (I) is the set of all unit nodes in the I-th block, set (J) is the set of all unit nodes in the J-th block,for the mth shift increment pattern of the unit node I in the ith block,for the l-th shift increment pattern, k, of cell node J in the J-th blockijFor the sake of a corresponding overall stiffness,the coefficients for the jth shift increment of the jth block,sifor the seismic source to exert pressure on cell node i,and B is the displacement approximate value of the unit node i before the displacement update, and B is the total number of the blocks.
The slave processor invoking unit 410 is configured to invoke a plurality of slave processors to process data of each of the blocks according to the following formula, respectively, where each slave processor corresponds to one of the blocks:
the master processor calling unit 412 is configured to call the master processor to receive processing results of the slave processors to solve the multi-stage parallel linear equation set.
In one possible embodiment, the dividing surfaces are placed in cells.
In a possible implementation, the finite element control equation establishing unit 402 is specifically configured to:
obtaining a sound wave equation under a rectangular coordinate system in a space-time domain:
wherein P is sound pressure, v is wave velocity, t is time variable, x, y, z are space variables, and s is seismic source function;
with partial dispersion, the following approximate solution is constructed:
wherein N isj(x, y, z) is an interpolation function at cell node j, abbreviated as N hereinafterj,dj(t) is the displacement of the cell node j, abbreviated as d hereinafterjNd is the total number of unit nodes;
substituting equation 2 into equation 1 yields the margin R:
according to the weighted residue method, an interpolation function of approximate solution is used as a weight function, and the weighted integral of the residue R on the solving area omega is zero, namely
∫ΩNiRd Ω ═ 0(i ═ 1,2, …, nd) formula 4
Substituting formula 3 into formula 4 yields:
dividing and integrating the three terms in the middle of the left end of the equal sign of the formula 5 to obtain
Wherein n isx、ny、nzIs the direction cosine of the boundary outer normal, Γ is the outer boundary of Ω;
substituting the formulas 2 and 6 into formula 5 to obtain a finite element control equation with discrete space:
In a possible implementation, the time integral extrapolation unit 404 is specifically configured to:
rewriting the finite element control equation into a matrix form:
wherein,
wherein N is an interpolation function matrix of each unit node, d is a displacement matrix of each unit node,is the displacement second order partial derivative matrix of each unit node, p is the density,is the Hamiltonian, ΩeIndicating that the integration area is a unit;
in the calculation along the time axis, time integration is achieved using an implicit NewMark device, which uses the following assumptions in the time domain:
further, it can be obtained from equation 10
Wherein,is a displacement first-order partial derivative matrix of each unit node,Δ t represents the time increment, subscripts t and t + Δ t represent the values of the time variable, and α and δ are constant coefficients;
by collating equations 8 and 9 according to equations 10 and 11, we obtain:
among them are:
after solving equation 12, update results in:
wherein,as a global stiffness matrix, c0、c1、c2、c3、c4、c5、c6、c7All are constant coefficients, and C is a damping matrix;
since the random velocity boundary condition is actually used, the damping matrix C does not appear in the calculation, and equation 12 is simplified as:
Wherein d is0A matrix of initial values representing the displacement of each node,and the initial value matrix represents the second-order partial derivatives of the displacement of each node.
In a possible implementation, the representing the displacement update value of each unit node in each partition by q displacement increments includes:
the displacement d of the node I in the cell in the I-th block is expressed byi:
Wherein,is the coefficient of the mth shift increment in the ith block,for the mth shift increment pattern of node I in the ith block
Application example
Figure 5 shows the sigbee 2 velocity model. The velocity model contains sedimentary layers of different durations, which are segmented into fault blocks of different sizes by a plurality of normal faults and reverse faults. In addition, a high-speed rock body with a complex structure is embedded in the model, and the high-speed rock body can cause the problems of insufficient illumination of the underlying stratum and difficult imaging. The model is often used to test the accuracy of new offset imaging algorithms.
The sigbee 2 model shown in fig. 5 is 24384 meters long and 9144 meters deep, with a grid size of 7.62 × 7.62 meters. In order to verify the precision and efficiency performance of the finite element reverse time migration method (HEP-FE-RTM) with high expansibility according to the application, one-way wave equation prestack depth migration (WEM), finite difference reverse time migration (FD-RTM) and finite element reverse time migration (HEP-FE-RTM) imaging results are analyzed in a comparative way.
FIG. 6 is a one-way wave equation prestack depth migration imaging result. It can be seen from the figure that depth migration based on the one-way wave equation has dip limitations, cannot be accurately imaged for the steeper formation dip regions (regions pointed by arrows in figure 6), and cannot be reached for the wavefield in the underlying region shielded by the high velocity salt (regions in the box in figure 6), and therefore cannot be imaged for the subsurface formation region.
Figure 7 shows a finite difference RTM outcome profile of the Sigsbee2 velocity model. Because the method adopts the two-way wave equation to describe wave field propagation, no dip angle limitation exists, and imaging of the complex structural area (the area in the box in the figure 7) under the rock is greatly improved compared with the depth migration of the one-way wave equation. However, since the finite difference method mesh subdivision adopts a rectangular or regular hexahedral isometric mesh, when the problem of severe lateral irregular change of a velocity field is handled, an ideal imaging effect cannot be achieved (fig. 7 shows a finite difference RTM result section of a Sigsbee2 velocity model).
Figure 8 shows a finite element RTM outcome profile of the Sigsbee2 velocity model obtained in accordance with the present application. Because the finite element method is adopted to solve the wave field propagation problem, the finite element can adopt a series of unit discrete solution domains such as triangles, tetrahedrons and the like, and can process any complex region, the imaging precision of the rock top interface (the region pointed by the arrow) and the imaging precision of the complex underground structure are high.
Imaging contrast analysis by three methods showed: the high-expansibility finite element reverse time migration method (HEP-FE-RTM) provided by the invention has higher imaging precision and can capture almost all complex construction details.
The above test was performed on an "eosin" cluster system, with an analysis field of 24384 meters horizontally and 9144 meters vertically for each shot; and the calculation domain after the random speed boundary is expanded is 27432 meters in the transverse direction and 10668 meters in the longitudinal direction. The rectangular grid size is 30.48 × 15.24 meters, a 5 × 5 solution domain blocking scheme is adopted, each block comprises 630000 triangular units, and 26 nodes are called on the 'eosin' cluster system for parallel calculation. The three methods, sigbee 2 model, took time to image as shown in table 1. It can be seen that: whether single shot migration or full data volume migration, reverse time migration is significantly more time consuming than single-pass wave migration. Compared with the limited difference method, the reverse time migration of the single shot consumes 46.5 minutes (one shot and one node, and 26 shots are simultaneously operated), the method of the invention consumes 1.9 minutes (one shot and 26 nodes, and one shot is operated at one time) of the single shot, and the reverse time migration of the method is almost equivalent to the reverse time migration of the conventional limited difference method, thereby realizing that the imaging precision is obviously improved and more complex structure detail characteristics can be captured under the condition that the time consumption is not increased compared with the reverse time migration of the conventional limited difference method.
TABLE 1 time-consuming comparison Table for Sigsbee2 model imaging of three methods (5X 5)
FIG. 9 shows a high-precision velocity model of complex small fault blocks in the east of China, which has a transverse direction of 9000 m, a longitudinal direction of 5000 m, a grid size of 10 × 10 m, and comprises a large high-steepness fault with a controlled concave boundary, a plurality of small faults with different dip angles and different directions, and seven reflecting layers from top to bottom. The forward shot gather data comprises 672 shots, and the analysis domain of each shot is 9000 meters in the transverse direction and 5000 meters in the longitudinal direction; the random speed boundary is respectively expanded by 1500 meters in the horizontal direction and the longitudinal direction (except a horizontal plane), the final offset extrapolation calculation domain is 12000 meters in the horizontal direction and 6500 meters in the longitudinal direction, an 8 multiplied by 8 solution domain blocking scheme is adopted, and 65 nodes are called on the 'eosin light' cluster system for parallel calculation.
The imaging results of the three methods are shown in fig. 10 to 12, respectively. As shown in FIG. 10, the one-way wave equation depth shift is not ideal for imaging in both the high and steep fault affected zone (elliptical zone) of the ascending disk and the complex and small fault block gathering zone (square zone) of the descending disk.
As shown in FIG. 11, the finite difference method reverse time migration adopts a two-way wave equation to describe wave field propagation, has no dip angle limitation, can realize rotating wave imaging, and has improved imaging precision in a high and steep large fault affected zone and a complex small fault aggregation zone, and particularly has the most obvious improvement on the contact relation between a horizon and a large fault in a high and steep large fault affected zone of a lifting disc. However, since finite difference method mesh generation adopts rectangular or regular hexahedral mesh, it cannot achieve ideal imaging effect (square area) when dealing with the complex small block aggregation area of the descending disk.
As shown in fig. 12, the finite element method is used for solving the wave field propagation problem, the solution domain unit is diversified in discrete mode, extremely complex velocity field distribution conditions such as complex small block aggregation areas can be processed, and the imaging result is obviously improved in time migration precision relative to finite difference inversion (block area).
Imaging contrast analysis by three methods showed: the method provided by the application has high imaging precision, can accurately depict the fracture system on the accurate imaging of the complex small fracture block in the east of China, and has obvious advantages.
The model tests three methods imaging time consuming pairs as shown in table 2. It can be seen that: whether single shot migration or full data volume migration, reverse time migration is significantly more time consuming than single-pass wave migration. The method disclosed by the application is almost equivalent in time consumption compared with the finite difference method, and the advantage of accurate imaging is more remarkable.
TABLE 2 time-consuming contrast chart for imaging of three methods, namely Subei basin model
(HEP-FE-RTM 8X 8 block scheme)
The present application may be a system, method and/or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions embodied thereon for causing a processor to implement various aspects of the present application.
The computer readable storage medium may be a tangible device that can hold and store the instructions for use by the instruction execution device. The computer readable storage medium may be, for example, but not limited to, an electronic memory device, a magnetic memory device, an optical memory device, an electromagnetic memory device, a semiconductor memory device, or any suitable combination of the foregoing. More specific examples (a non-exhaustive list) of the computer readable storage medium would include the following: a portable computer diskette, a hard disk, a Random Access Memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), a Static Random Access Memory (SRAM), a portable compact disc read-only memory (CD-ROM), a Digital Versatile Disc (DVD), a memory stick, a floppy disk, a mechanical coding device, such as punch cards or in-groove projection structures having instructions stored thereon, and any suitable combination of the foregoing. Computer-readable storage media as used herein is not to be construed as transitory signals per se, such as radio waves or other freely propagating electromagnetic waves, electromagnetic waves propagating through a waveguide or other transmission medium (e.g., optical pulses through a fiber optic cable), or electrical signals transmitted through electrical wires.
Various aspects of the present application 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 application. It will 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.
Having described embodiments of the present application, the foregoing description is intended to be exemplary, not exhaustive, and not limited to the disclosed embodiments. Many modifications and variations will be apparent to those of ordinary skill in the art without departing from the scope and spirit of the described embodiments. The terms used herein were chosen in order to best explain the principles of the embodiments, the practical application, or technical improvements to the techniques in the marketplace, or to enable others of ordinary skill in the art to understand the embodiments disclosed herein.
Claims (10)
1. A finite element-based reverse time migration method, the method comprising:
dividing a continuous three-dimensional space into a limited number of units which are mutually communicated but do not coincide, and establishing a finite element control equation with discrete space;
extrapolating a finite element control equation according to a time integration method to obtain a linear equation set with discrete space and time;
dividing the three-dimensional space into a plurality of blocks through a dividing surface, and expressing the displacement update value of each unit node in each block by q displacement increments;
and transforming the linear equation set according to the minimum potential energy principle, and substituting the displacement update values which use q displacement increments to represent unit nodes to obtain a multi-stage parallel linear equation set:
wherein,set (I) is the set of all unit nodes in the I-th block, set (J) is the set of all unit nodes in the J-th block,for the mth shift increment pattern of the unit node I in the ith block,for the l-th shift increment pattern, k, of cell node J in the J-th blockijFor the sake of a corresponding overall stiffness,the coefficients for the jth shift increment of the jth block,sifor the seismic source to exert pressure on cell node i,the displacement approximate value of the unit node i before the displacement updating is performed, and B is the total number of the blocks;
calling a plurality of slave processors to process data of each block according to the following formula, wherein each slave processor corresponds to one block:
and calling the master processor to receive the processing result of each slave processor so as to solve the multi-stage parallel linear equation system.
2. The method of claim 1, wherein the splitting planes are placed within cells.
3. The method of claim 1, wherein the dividing the continuous three-dimensional space into a finite number of elements that are interconnected but do not coincide and establishing spatially discrete finite element control equations comprises:
obtaining a sound wave equation under a rectangular coordinate system in a space-time domain:
wherein P is sound pressure, v is wave velocity, t is time variable, x, y, z are space variables, and s is seismic source function;
with partial dispersion, the following approximate solution is constructed:
wherein N isj(x, y, z) is an interpolation function at cell node j, abbreviated as N hereinafterj,dj(t) is the displacement of the cell node j, abbreviated as d hereinafterjNd is the total number of unit nodes;
substituting equation 2 into equation 1 yields the margin R:
according to the weighted residue method, an interpolation function of approximate solution is used as a weight function, and the weighted integral of the residue R on the solving area omega is zero, namely
∫ΩNiRd Ω ═ 0(i ═ 1,2, …, nd) formula 4
Substituting formula 3 into formula 4 yields:
dividing and integrating the three terms in the middle of the left end of the equal sign of the formula 5 to obtain
Wherein n isx、ny、nzIs the direction cosine of the boundary outer normal, Γ is the outer boundary of Ω;
substituting the formulas 2 and 6 into formula 5 to obtain a finite element control equation with discrete space:
4. The method of claim 3, wherein extrapolating the finite element control equation according to the time integration method to obtain a spatially and temporally discrete system of linear equations comprises:
rewriting the finite element control equation into a matrix form:
wherein,
wherein N is an interpolation function matrix of each unit node, d is a displacement matrix of each unit node,is the displacement second order partial derivative matrix of each unit node, p is the density,is the Hamiltonian, ΩeIndicating that the integration area is a unit;
in the calculation along the time axis, the implicit NewMark method is used to achieve time integration, and the NewMark method uses the following assumptions in the time region:
further, it can be obtained from equation 10
Wherein,is a displacement first-order partial derivative matrix of each unit node, delta t represents time increment, subscript t and t + delta t represent values of time variables, and alpha and delta are constant coefficients;
by collating equations 8 and 9 according to equations 10 and 11, we obtain:
among them are:
after solving equation 12, update results in:
wherein,as a global stiffness matrix, c0、c1、c2、c3、c4、c5、c6、c7All are constant coefficients, and C is a damping matrix;
since the random velocity boundary condition is actually used, the damping matrix C does not appear in the calculation, and equation 12 is simplified as:
5. The method of claim 1, wherein said representing the displacement update values of the respective unit nodes in each partition by q displacement increments comprises:
the displacement d of the node I in the cell in the I-th block is expressed byi:
6. A finite element-based reverse time migration apparatus, the apparatus comprising:
the finite element control equation establishing unit is used for dividing a continuous three-dimensional space into a limited number of units which are mutually communicated but do not coincide, and establishing a finite element control equation with discrete space;
the time integral extrapolation unit is used for extrapolating a finite element control equation according to a time integral method to obtain a linear equation set with discrete space and time;
the block unit is used for dividing the three-dimensional space into a plurality of blocks through a dividing surface, and expressing the displacement update value of each unit node in each block by q displacement increments;
the multilevel parallel equation establishing unit is used for transforming the linear equation set according to the minimum potential energy principle and substituting the displacement update value which uses q displacement increments to represent unit nodes to obtain a multilevel parallel linear equation set:
wherein,set (I) is the set of all unit nodes in the I-th block, set (J) is the set of all unit nodes in the J-th block,for the mth shift increment pattern of the unit node I in the ith block,for the l-th shift increment pattern, k, of cell node J in the J-th blockijFor the sake of a corresponding overall stiffness,the coefficients for the jth shift increment of the jth block,sifor the seismic source to exert pressure on cell node i,the displacement approximate value of the unit node i before the displacement updating is performed, and B is the total number of the blocks;
a slave processor calling unit, configured to call a plurality of slave processors to process data of each partition according to the following formula, where each slave processor corresponds to one partition:
and the master control processor calling unit is used for calling the master control processor to receive the processing result of each slave processor so as to solve the multi-stage parallel linear equation set.
7. The apparatus of claim 6, wherein the dividing surface is disposed within a cell.
8. The apparatus according to claim 6, wherein the finite element control equation establishing unit is specifically configured to:
obtaining a sound wave equation under a rectangular coordinate system in a space-time domain:
wherein P is sound pressure, v is wave velocity, t is time variable, x, y, z are space variables, and s is seismic source function;
with partial dispersion, the following approximate solution is constructed:
wherein N isj(x, y, z) is an interpolation function at cell node j, abbreviated as N hereinafterj,dj(t) is the displacement of the cell node j, abbreviated as d hereinafterjNd is the total number of unit nodes;
substituting equation 2 into equation 1 yields the margin R:
according to the weighted residue method, an interpolation function of approximate solution is used as a weight function, and the weighted integral of the residue R on the solving area omega is zero, namely
∫ΩNiRd Ω ═ 0(i ═ 1,2, …, nd) formula 4
Substituting formula 3 into formula 4 yields:
dividing and integrating the three terms in the middle of the left end of the equal sign of the formula 5 to obtain
Wherein n isx、ny、nzIs the direction cosine of the boundary outer normal, Γ is the outer boundary of Ω;
substituting the formulas 2 and 6 into formula 5 to obtain a finite element control equation with discrete space:
9. The apparatus according to claim 8, wherein the time-integrated extrapolation unit is specifically configured to:
rewriting the finite element control equation into a matrix form:
wherein,
wherein N is an interpolation function matrix of each unit node, d is a displacement matrix of each unit node,is the displacement second order partial derivative matrix of each unit node, p is the density,is the Hamiltonian, ΩeIndicating that the integration area is a unit;
in the calculation along the time axis, time integration is achieved using an implicit NewMark device, which uses the following assumptions in the time domain:
further, it can be obtained from equation 10
Wherein,is a displacement first-order partial derivative matrix of each unit node, delta t represents time increment, subscript t and t + delta t represent values of time variables, and alpha and delta are constant coefficients;
by collating equations 8 and 9 according to equations 10 and 11, we obtain:
among them are:
after solving equation 12, update results in:
wherein,as a global stiffness matrix, c0、c1、c2、c3、c4、c5、c6、c7All are constant coefficients, and C is a damping matrix;
since the random velocity boundary condition is actually used, the damping matrix C does not appear in the calculation, and equation 12 is simplified as:
10. The apparatus of claim 6, wherein said representing the displacement update values of the respective unit nodes in each partition by q displacement increments comprises:
the displacement d of the node I in the cell in the I-th block is expressed byi:
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910973951.XA CN112731533A (en) | 2019-10-14 | 2019-10-14 | Reverse time migration method and device based on finite element |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910973951.XA CN112731533A (en) | 2019-10-14 | 2019-10-14 | Reverse time migration method and device based on finite element |
Publications (1)
Publication Number | Publication Date |
---|---|
CN112731533A true CN112731533A (en) | 2021-04-30 |
Family
ID=75588559
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910973951.XA Pending CN112731533A (en) | 2019-10-14 | 2019-10-14 | Reverse time migration method and device based on finite element |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112731533A (en) |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20120051179A1 (en) * | 2010-08-26 | 2012-03-01 | Snu R&Db Foundation | Method and apparatus for time-domain reverse-time migration with source estimation |
US20150127311A1 (en) * | 2013-11-06 | 2015-05-07 | Weidlinger Associates, Inc. | Computer Implemented Apparatus and Method for Finite Element Modeling Using Hybrid Absorbing Element |
CN105717539A (en) * | 2016-01-28 | 2016-06-29 | 中国地质大学(北京) | Multi GPU calculation based reverse time migration imaging method of 3D TTI medium |
-
2019
- 2019-10-14 CN CN201910973951.XA patent/CN112731533A/en active Pending
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20120051179A1 (en) * | 2010-08-26 | 2012-03-01 | Snu R&Db Foundation | Method and apparatus for time-domain reverse-time migration with source estimation |
US20150127311A1 (en) * | 2013-11-06 | 2015-05-07 | Weidlinger Associates, Inc. | Computer Implemented Apparatus and Method for Finite Element Modeling Using Hybrid Absorbing Element |
CN105717539A (en) * | 2016-01-28 | 2016-06-29 | 中国地质大学(北京) | Multi GPU calculation based reverse time migration imaging method of 3D TTI medium |
Non-Patent Citations (5)
Title |
---|
LIMIN LIU ET AL.: "A high expansion implicit finite-element prestack reverse time migration method", 《EXPLORATION GEOPHYSICS》 * |
孙培德等: "《环境系统模型及数值模拟》", 31 December 2005, 北京:中国环境科学出版社 * |
王月英: "基于MPI的三维波动方程有限元法并行正演模拟", 《石油物探》 * |
薛东川等: "波动方程有限元叠前逆时偏移", 《石油地球物理勘探》 * |
韩德超: "有限元地震波正演模拟及反演方法研究", 《中国优秀博硕士学位论文全文数据库(硕士) 基础科学辑》 * |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CA2717889C (en) | Method for building a depositional space corresponding to a geological domain | |
Martin et al. | Gravity inversion using wavelet-based compression on parallel hybrid CPU/GPU systems: application to southwest Ghana | |
CA2821482C (en) | Systems and methods for two-dimensional domain decomposition during parallel reservoir simulation | |
CA2862900C (en) | Multi-level solution of large-scale linear systems in simulation of porous media in giant reservoirs | |
WO2009079088A1 (en) | Modeling subsurface processes on unstructured grid | |
CN114117906A (en) | Multi-scale unsupervised seismic wave velocity inversion method based on observation data self-encoding | |
CA2816815A1 (en) | Systems and methods for generating updates of geological models | |
Yagawa | Node‐by‐node parallel finite elements: a virtually meshless method | |
CN105467443A (en) | A three-dimensional anisotropy elastic wave numerical simulation method and system | |
CN106646645A (en) | Novel gravity forward acceleration method | |
EP2603878A1 (en) | Reservoir upscaling method with preserved transmissibility | |
CN104597488B (en) | Optimum design method of finite difference template of non-equiangular long-grid wave equation | |
Muratov et al. | Grid-characteristic method on unstructured tetrahedral meshes | |
Fujisawa et al. | Parallel computing of high‐speed compressible flows using a node‐based finite‐element method | |
US9928315B2 (en) | Re-ordered interpolation and convolution for faster staggered-grid processing | |
CN108983290B (en) | It is a kind of three-dimensional Method in Transverse Isotropic Medium in travel when determine method and system | |
CN109490948B (en) | Seismic acoustic wave equation vector parallel computing method | |
CN112731533A (en) | Reverse time migration method and device based on finite element | |
CN109188522B (en) | Velocity field construction method and device | |
CN109001804B (en) | Method, device and system for determining effective force based on three-dimensional seismic data | |
CN105974471A (en) | Seismic data multi-GPU fast forward computation method based on asynchronous flow | |
US10254441B2 (en) | Method of modelling a subsurface volume | |
US20090043549A1 (en) | Methods, apparatus, and products for seismic ray tracing | |
CN117150835A (en) | Method for simulating seismic wave value of viscous-acoustic anisotropic medium | |
CN110162804B (en) | Wave field forward modeling optimization method based on CPU acceleration |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
RJ01 | Rejection of invention patent application after publication |
Application publication date: 20210430 |
|
RJ01 | Rejection of invention patent application after publication |