CN112731533A - Reverse time migration method and device based on finite element - Google Patents

Reverse time migration method and device based on finite element Download PDF

Info

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
Application number
CN201910973951.XA
Other languages
Chinese (zh)
Inventor
刘立民
刘定进
姚晓龙
李博
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China Petroleum and Chemical Corp, Sinopec Geophysical Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201910973951.XA priority Critical patent/CN112731533A/en
Publication of CN112731533A publication Critical patent/CN112731533A/en
Pending legal-status Critical Current

Links

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

Reverse time migration method and device based on finite element
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:
Figure BDA0002233010220000011
wherein,
Figure BDA0002233010220000021
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,
Figure BDA0002233010220000022
for the mth shift increment pattern of the unit node I in the ith block,
Figure BDA0002233010220000023
for the l-th shift increment pattern, k, of cell node J in the J-th blockijFor the sake of a corresponding overall stiffness,
Figure BDA0002233010220000024
the coefficients for the jth shift increment of the jth block,
Figure BDA0002233010220000025
sifor the seismic source to exert pressure on cell node i,
Figure BDA0002233010220000026
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:
Figure BDA0002233010220000027
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:
Figure BDA0002233010220000028
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:
Figure BDA0002233010220000029
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:
Figure BDA00022330102200000210
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:
Figure BDA0002233010220000031
dividing and integrating the three terms in the middle of the left end of the equal sign of the formula 5 to obtain
Figure BDA0002233010220000032
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:
Figure BDA0002233010220000033
wherein,
Figure BDA0002233010220000034
is djSecond order partial derivatives of (d).
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:
Figure BDA0002233010220000035
wherein,
Figure BDA0002233010220000036
wherein N is an interpolation function matrix of each unit node, d is a displacement matrix of each unit node,
Figure BDA0002233010220000037
is the displacement second order partial derivative matrix of each unit node, p is the density,
Figure BDA0002233010220000038
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:
Figure BDA0002233010220000041
further, it can be obtained from equation 10
Figure BDA0002233010220000042
Wherein,
Figure BDA0002233010220000043
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:
Figure BDA0002233010220000044
among them are:
Figure BDA0002233010220000045
after solving equation 12, update results in:
Figure BDA0002233010220000046
wherein,
Figure BDA0002233010220000047
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:
Figure BDA0002233010220000048
the boundary condition is
Figure BDA0002233010220000049
Wherein d is0A matrix of initial values representing the displacement of each node,
Figure BDA00022330102200000410
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：
Figure BDA0002233010220000051
Wherein,
Figure BDA0002233010220000052
is the coefficient of the mth shift increment in the ith block,
Figure BDA0002233010220000053
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:
Figure BDA0002233010220000054
wherein,
Figure BDA0002233010220000055
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,
Figure BDA0002233010220000056
for the mth shift increment pattern of the unit node I in the ith block,
Figure BDA0002233010220000057
for the l-th shift increment pattern, k, of cell node J in the J-th blockijFor the sake of a corresponding overall stiffness,
Figure BDA0002233010220000058
the coefficients for the jth shift increment of the jth block,
Figure BDA0002233010220000059
sifor the seismic source to exert pressure on cell node i,
Figure BDA00022330102200000512
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:
Figure BDA00022330102200000511
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:
Figure BDA0002233010220000061
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:
Figure BDA0002233010220000062
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:
Figure BDA0002233010220000063
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:
Figure BDA0002233010220000064
dividing and integrating the three terms in the middle of the left end of the equal sign of the formula 5 to obtain
Figure BDA0002233010220000065
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:
Figure BDA0002233010220000071
wherein,
Figure BDA0002233010220000072
is djSecond order partial derivatives of (d).
In a possible implementation, the time-integral extrapolation unit is specifically configured to:
rewriting the finite element control equation into a matrix form:
Figure BDA0002233010220000073
wherein,
Figure BDA0002233010220000074
wherein N is an interpolation function matrix of each unit node, d is a displacement matrix of each unit node,
Figure BDA0002233010220000075
is the displacement second order partial derivative matrix of each unit node, p is the density,
Figure BDA0002233010220000076
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:
Figure BDA0002233010220000077
further, it can be obtained from equation 10
Figure BDA0002233010220000078
Wherein,
Figure BDA0002233010220000079
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:
Figure BDA00022330102200000710
among them are:
Figure BDA0002233010220000081
after solving equation 12, update results in:
Figure BDA0002233010220000082
wherein,
Figure BDA0002233010220000083
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:
Figure BDA0002233010220000084
the boundary condition is
Figure BDA0002233010220000085
Wherein d is0A matrix of initial values representing the displacement of each node,
Figure BDA0002233010220000086
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：
Figure BDA0002233010220000087
Wherein,
Figure BDA0002233010220000088
is the coefficient of the mth shift increment in the ith block,
Figure BDA0002233010220000089
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.
Step 102, dividing the 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.
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:
Figure BDA0002233010220000091
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:
Figure BDA0002233010220000092
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:
Figure BDA0002233010220000101
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:
Figure BDA0002233010220000102
dividing and integrating the three terms in the middle of the left end of the equal sign of the formula 5 to obtain
Figure BDA0002233010220000103
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:
Figure BDA0002233010220000104
wherein,
Figure BDA0002233010220000107
is djSecond order partial derivatives of (d).
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:
Figure BDA0002233010220000106
wherein,
Figure BDA0002233010220000111
wherein N is an interpolation function matrix of each unit node, d is a displacement matrix of each unit node,
Figure BDA0002233010220000112
is the displacement second order partial derivative matrix of each unit node, p is the density,
Figure BDA0002233010220000113
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:
Figure BDA0002233010220000114
further, it can be obtained from equation 10
Figure BDA0002233010220000115
Wherein,
Figure BDA0002233010220000116
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:
Figure BDA0002233010220000117
among them are:
Figure BDA0002233010220000118
after solving equation 12, update results in:
Figure BDA0002233010220000119
wherein,
Figure BDA00022330102200001110
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:
Figure BDA0002233010220000121
the boundary condition is
Figure BDA0002233010220000122
Wherein d is0A matrix of initial values representing the displacement of the respective nodes,
Figure BDA0002233010220000123
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：
Figure BDA0002233010220000124
Wherein,
Figure BDA0002233010220000125
updating the displacement approximation of the unit node i before the displacement,
Figure BDA0002233010220000128
is the coefficient of the mth shift increment in the ith block,
Figure BDA0002233010220000126
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:
Figure BDA0002233010220000127
wherein,
Figure BDA0002233010220000131
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,
Figure BDA0002233010220000132
for the mth shift increment pattern of the unit node I in the ith block,
Figure BDA0002233010220000133
for the l-th shift increment pattern, k, of cell node J in the J-th blockijFor the sake of a corresponding overall stiffness,
Figure BDA0002233010220000134
the coefficients for the jth shift increment of the jth block,
Figure BDA0002233010220000135
sifor the seismic source to exert pressure on cell node i,
Figure BDA0002233010220000136
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
Figure BDA0002233010220000137
Wherein k isijAs a global stiffness matrix
Figure BDA0002233010220000138
The corresponding element in (1).
Equation 16 is based on the principle of minimum potential energy, and the total potential energy is expressed in the form of
Figure BDA0002233010220000139
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:
Figure BDA00022330102200001310
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:
Figure BDA00022330102200001311
step 110, calling a plurality of slave processors to process the data of each block according to the following formula, wherein each slave processor corresponds to one block:
Figure BDA00022330102200001312
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 processor
Figure BDA0002233010220000141
And 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:
Figure BDA0002233010220000142
wherein,
Figure BDA0002233010220000151
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,
Figure BDA0002233010220000152
for the mth shift increment pattern of the unit node I in the ith block,
Figure BDA0002233010220000153
for the l-th shift increment pattern, k, of cell node J in the J-th blockijFor the sake of a corresponding overall stiffness,
Figure BDA0002233010220000154
the coefficients for the jth shift increment of the jth block,
Figure BDA0002233010220000155
sifor the seismic source to exert pressure on cell node i,
Figure BDA0002233010220000156
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:
Figure BDA0002233010220000157
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:
Figure BDA0002233010220000158
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:
Figure BDA0002233010220000159
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:
Figure BDA00022330102200001510
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:
Figure BDA0002233010220000161
dividing and integrating the three terms in the middle of the left end of the equal sign of the formula 5 to obtain
Figure BDA0002233010220000162
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:
Figure BDA0002233010220000163
wherein,
Figure BDA0002233010220000164
is djSecond order partial derivatives of (d).
In a possible implementation, the time integral extrapolation unit 404 is specifically configured to:
rewriting the finite element control equation into a matrix form:
Figure BDA0002233010220000165
wherein,
Figure BDA0002233010220000166
wherein N is an interpolation function matrix of each unit node, d is a displacement matrix of each unit node,
Figure BDA0002233010220000167
is the displacement second order partial derivative matrix of each unit node, p is the density,
Figure BDA0002233010220000168
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:
Figure BDA0002233010220000171
further, it can be obtained from equation 10
Figure BDA0002233010220000172
Wherein,
Figure BDA0002233010220000173
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:
Figure BDA0002233010220000174
among them are:
Figure BDA0002233010220000175
after solving equation 12, update results in:
Figure BDA0002233010220000176
wherein,
Figure BDA0002233010220000177
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:
Figure BDA0002233010220000178
the boundary condition is
Figure BDA0002233010220000179
Wherein d is0A matrix of initial values representing the displacement of each node,
Figure BDA00022330102200001710
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：
Figure BDA0002233010220000181
Wherein,
Figure BDA0002233010220000182
is the coefficient of the mth shift increment in the ith block,
Figure BDA0002233010220000183
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)
Figure BDA0002233010220000191
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)
Figure BDA0002233010220000201
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:
Figure FDA0002233010210000011
wherein,
Figure FDA0002233010210000012
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,
Figure FDA0002233010210000013
for the mth shift increment pattern of the unit node I in the ith block,
Figure FDA0002233010210000014
for the l-th shift increment pattern, k, of cell node J in the J-th blockijFor the sake of a corresponding overall stiffness,
Figure FDA0002233010210000015
the coefficients for the jth shift increment of the jth block,
Figure FDA0002233010210000016
sifor the seismic source to exert pressure on cell node i,
Figure FDA0002233010210000017
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:
Figure FDA0002233010210000018
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:
Figure FDA0002233010210000021
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:
Figure FDA0002233010210000022
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:
Figure FDA0002233010210000023
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:
Figure FDA0002233010210000024
dividing and integrating the three terms in the middle of the left end of the equal sign of the formula 5 to obtain
Figure FDA0002233010210000025
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:
Figure FDA0002233010210000031
wherein,
Figure FDA0002233010210000032
is djSecond order partial derivatives of (d).
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:
Figure FDA0002233010210000033
wherein,
Figure FDA0002233010210000034
wherein N is an interpolation function matrix of each unit node, d is a displacement matrix of each unit node,
Figure FDA0002233010210000035
is the displacement second order partial derivative matrix of each unit node, p is the density,
Figure FDA0002233010210000036
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:
Figure FDA0002233010210000037
further, it can be obtained from equation 10
Figure FDA0002233010210000038
Wherein,
Figure FDA0002233010210000039
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:
Figure FDA0002233010210000041
among them are:
Figure FDA0002233010210000042
after solving equation 12, update results in:
Figure FDA0002233010210000043
wherein,
Figure FDA0002233010210000044
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:
Figure FDA0002233010210000045
the boundary condition is
Figure FDA0002233010210000046
Wherein d is0A matrix of initial values representing the displacement of each node,
Figure FDA0002233010210000047
and the initial value matrix represents the second-order partial derivatives of the displacement of each node.
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：
Figure FDA0002233010210000048
Wherein,
Figure FDA0002233010210000049
is the coefficient of the mth shift increment in the ith block,
Figure FDA00022330102100000410
is the mth shift increment pattern for node I in the ith partition.
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:
Figure FDA0002233010210000051
wherein,
Figure FDA0002233010210000052
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,
Figure FDA0002233010210000053
for the mth shift increment pattern of the unit node I in the ith block,
Figure FDA0002233010210000054
for the l-th shift increment pattern, k, of cell node J in the J-th blockijFor the sake of a corresponding overall stiffness,
Figure FDA0002233010210000055
the coefficients for the jth shift increment of the jth block,
Figure FDA0002233010210000056
sifor the seismic source to exert pressure on cell node i,
Figure FDA0002233010210000057
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:
Figure FDA0002233010210000058
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:
Figure FDA0002233010210000059
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:
Figure FDA0002233010210000061
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:
Figure FDA0002233010210000062
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:
Figure FDA0002233010210000063
dividing and integrating the three terms in the middle of the left end of the equal sign of the formula 5 to obtain
Figure FDA0002233010210000064
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:
Figure FDA0002233010210000065
wherein,
Figure FDA0002233010210000066
is djSecond order partial derivatives of (d).
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:
Figure FDA0002233010210000071
wherein,
Figure FDA0002233010210000072
wherein N is an interpolation function matrix of each unit node, d is a displacement matrix of each unit node,
Figure FDA0002233010210000073
is the displacement second order partial derivative matrix of each unit node, p is the density,
Figure FDA0002233010210000074
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:
Figure FDA0002233010210000075
further, it can be obtained from equation 10
Figure FDA0002233010210000076
Wherein,
Figure FDA0002233010210000077
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:
Figure FDA0002233010210000078
among them are:
Figure FDA0002233010210000079
after solving equation 12, update results in:
Figure FDA0002233010210000081
wherein,
Figure FDA0002233010210000082
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:
Figure FDA0002233010210000083
the boundary condition is
Figure FDA0002233010210000084
Wherein d is0A matrix of initial values representing the displacement of each node,
Figure FDA0002233010210000085
and the initial value matrix represents the second-order partial derivatives of the displacement of each node.
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：
Figure FDA0002233010210000086
Wherein,
Figure FDA0002233010210000087
is the coefficient of the mth shift increment in the ith block,
Figure FDA0002233010210000088
is the mth shift increment pattern for node I in the ith partition.
CN201910973951.XA 2019-10-14 2019-10-14 Reverse time migration method and device based on finite element Pending CN112731533A (en)

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)

* Cited by examiner, † Cited by third party
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

Patent Citations (3)

* Cited by examiner, † Cited by third party
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)

* Cited by examiner, † Cited by third party
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