CN112859162B - Three-dimensional acoustic wave equation space domain forward modeling method and device based on pseudo-spectrum method - Google Patents
Three-dimensional acoustic wave equation space domain forward modeling method and device based on pseudo-spectrum method Download PDFInfo
- Publication number
- CN112859162B CN112859162B CN201911187548.0A CN201911187548A CN112859162B CN 112859162 B CN112859162 B CN 112859162B CN 201911187548 A CN201911187548 A CN 201911187548A CN 112859162 B CN112859162 B CN 112859162B
- Authority
- CN
- China
- Prior art keywords
- order derivative
- derivative operator
- operator
- spatial
- dimensional
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 48
- 238000001228 spectrum Methods 0.000 title claims abstract description 18
- 238000004364 calculation method Methods 0.000 claims abstract description 60
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 31
- 238000005070 sampling Methods 0.000 claims abstract description 19
- 230000009466 transformation Effects 0.000 claims abstract description 17
- 239000011159 matrix material Substances 0.000 claims description 24
- 238000004590 computer program Methods 0.000 claims description 13
- 238000005457 optimization Methods 0.000 claims description 7
- 230000008569 process Effects 0.000 claims description 7
- 238000003860 storage Methods 0.000 claims description 6
- 238000005429 filling process Methods 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 13
- 238000006243 chemical reaction Methods 0.000 description 8
- 238000004458 analytical method Methods 0.000 description 6
- 230000006870 function Effects 0.000 description 6
- 238000009795 derivation Methods 0.000 description 5
- 238000004611 spectroscopical analysis Methods 0.000 description 4
- 230000006872 improvement Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 125000004122 cyclic group Chemical group 0.000 description 1
- 238000003874 inverse correlation nuclear magnetic resonance spectroscopy Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
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. analysis, for interpretation, for correction
- G01V1/30—Analysis
- G01V1/301—Analysis for determining seismic cross-sections or geostructures
- G01V1/302—Analysis for determining seismic cross-sections or geostructures in 3D data cubes
-
- 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. analysis, for interpretation, for correction
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/624—Reservoir parameters
Abstract
The application discloses a three-dimensional acoustic wave equation space domain forward modeling method and device based on a pseudo-spectrum method, wherein the method comprises the following steps: aiming at a specified coordinate direction of a three-dimensional space, calculating a second-order derivative operator of the specified coordinate direction according to the sampling point number and an m-order derivative operator formula of the specified coordinate direction; the m-order derivative operator formula is obtained by deduction according to a Fourier positive transformation algorithm and an inverse Fourier transformation algorithm, and m is more than or equal to 1; and determining the spatial partial derivative of the three-dimensional wave equation according to the second order derivative operator and the spatial wave field, and taking the spatial partial derivative as a spatial domain forward result of the three-dimensional wave equation. The application can reduce the calculation amount of the computer, is more suitable for GPU parallel calculation and improves the calculation efficiency.
Description
Technical Field
The application relates to the technical field of petroleum exploration, in particular to a three-dimensional acoustic wave equation space domain forward modeling method and device based on a pseudo-spectrum method.
Background
This section is intended to provide a background or context to the embodiments of the invention that are recited in the claims. The description herein is not admitted to be prior art by inclusion in this section.
Pseudo-spectroscopy is an important method for forward modeling of earthquakes. In pseudo-spectroscopy, the time domain finds the second partial derivative of the wavefield with respect to time by finite difference, and the space domain finds the partial derivative of the wavefield with respect to space by: converting the spatial domain to the wavenumber domain by a fast algorithm (Fast Fourier Transformation, FFT) of the discrete fourier transform; then multiplying the space domain shifted to the wavenumber domain by the square of the negative wavenumber; the second partial derivative of the wavefield with respect to space is then found by inverse fast fourier transform (INVERSE FAST Fourier Transform, IFFT).
The three-dimensional acoustic wave equation is:
Where u (x, y, z, t) is the pressure wavefield and v (x, y, z) is the velocity field in space.
When the pseudo-spectrum method is used for carrying out earthquake forward modeling on the space domain of the three-dimensional acoustic wave equation, the space partial derivatives in the directions of the x-axis, the y-axis and the z-axis of the three-dimensional space coordinate are all realized in 3 steps, and taking the direction of the x-axis as an example, the calculation steps are as follows: the space domain is subjected to FFT, the FFT result is multiplied by-k x 2, and the FFT result after being multiplied by-k x 2 is subjected to IFFT. Each step needs to perform a large amount of operations during the operation, which consumes a large amount of time for a computer, and simultaneously, a plurality of steps execute parallel computing processes which are difficult to adapt to a graphics processor (Graphics Processing Unit, GPU) one by one, so that the operation efficiency is reduced.
Disclosure of Invention
The embodiment of the application provides a three-dimensional acoustic wave equation space domain forward modeling method and device based on a pseudo-spectrum method, which are used for reducing the calculated amount of a computer, are more suitable for GPU parallel calculation and improve the operation efficiency, and comprise the following steps:
aiming at a specified coordinate direction of a three-dimensional space, calculating a second-order derivative operator of the specified coordinate direction according to the sampling point number and an m-order derivative operator formula of the specified coordinate direction; the m-order derivative operator formula is obtained by deduction according to a Fourier positive transformation algorithm and an inverse Fourier transformation algorithm, and m is more than or equal to 1; and determining the spatial partial derivative of the three-dimensional wave equation according to the second order derivative operator and the spatial wave field, and taking the spatial partial derivative as a spatial domain forward result of the three-dimensional wave equation.
The embodiment of the application also provides a three-dimensional acoustic wave equation space domain forward modeling device based on a pseudo-spectrum method, which is used for reducing the calculation amount of a computer, is more suitable for GPU parallel calculation and improves the calculation efficiency, and comprises the following steps:
The calculation module is used for calculating a second order derivative operator of the specified coordinate direction according to the sampling point number and the m order derivative operator formula of the specified coordinate direction aiming at the specified coordinate direction of the three-dimensional space; the m-order derivative operator formula is obtained by deduction according to a Fourier positive transformation algorithm and an inverse Fourier transformation algorithm, and m is more than or equal to 1; and the determining module is used for determining the spatial partial derivative of the three-dimensional wave equation according to the second-order derivative operator and the spatial wave field obtained by calculation of the calculating module and taking the spatial partial derivative as a forward result of the spatial domain of the three-dimensional wave equation.
In the embodiment of the application, an m-order derivative operator formula is obtained through derivation by a Fourier forward conversion algorithm and a Fourier inverse conversion algorithm, a second-order derivative operator in a specified coordinate direction can be directly determined through the m-order derivative operator formula, and then the spatial partial derivative of the three-dimensional wave equation can be determined according to the second-order derivative operator and the spatial wave field. Compared with the pseudo-spectrum method in the prior art, FFT calculation and IFFT calculation are not needed, and meanwhile, the calculation steps are reduced, so that the calculation amount of a computer is reduced, the calculation efficiency is improved, and the calculation of a single step is more suitable for GPU parallel calculation.
Drawings
In order to more clearly illustrate the embodiments of the application or the technical solutions in the prior art, the drawings that are required in the embodiments or the description of the prior art will be briefly described, it being obvious that the drawings in the following description are only some embodiments of the application, and that other drawings may be obtained according to these drawings without inventive effort for a person skilled in the art. In the drawings:
FIG. 1 is a flow chart of a three-dimensional acoustic wave equation space domain forward modeling method based on a pseudo-spectrum method in an embodiment of the application;
FIG. 2 is a flow chart of another three-dimensional acoustic wave equation space domain forward modeling method based on pseudo-spectroscopy in an embodiment of the application;
fig. 3 (a) to 3 (e) are schematic diagrams of intermediate point derivative operator curves with n=128 to 2048 in the embodiment of the present application;
Fig. 4 (a) to fig. 4 (f) are schematic diagrams of intermediate point derivative operator curves after extracting about 50 points of the intermediate point in the embodiment of the present application;
FIGS. 5 (a) -5 (c) are schematic diagrams of single shot analysis of operator cutoff lengths 64, 128, and 256, respectively, in an embodiment of the present application;
FIG. 5 (d) is a schematic diagram of a single shot analysis of an optimized truncator in an embodiment of the present application;
FIG. 6 is a schematic diagram of a velocity model according to an embodiment of the present application;
FIG. 7 is a schematic diagram of a three-dimensional forward single shot display in an embodiment of the present application;
fig. 8 is a schematic structural diagram of a three-dimensional acoustic wave equation space domain forward modeling device based on a pseudo-spectrum method in an embodiment of the application.
Detailed Description
For the purpose of making the objects, technical solutions and advantages of the embodiments of the present application more apparent, the embodiments of the present application will be described in further detail with reference to the accompanying drawings. The exemplary embodiments of the present application and their descriptions herein are for the purpose of explaining the present application, but are not to be construed as limiting the application.
The embodiment of the application provides a three-dimensional acoustic wave equation space domain forward modeling method based on a pseudo-spectrum method, which comprises the following steps of 101 and 102 as shown in fig. 1:
Step 101, calculating a second order derivative operator of the specified coordinate direction according to the sampling point number and the m order derivative operator formula of the specified coordinate direction aiming at the specified coordinate direction of the three-dimensional space.
Wherein, the m-order derivative operator formula is obtained by deduction according to a Fourier positive transformation algorithm and an inverse Fourier transformation algorithm, and m is more than or equal to 1.
The m-order derivative operator matrix in the x-axis direction is obtained through the following derivation process:
The formula of the fourier positive transform algorithm is as follows:
the formula of the inverse fourier transform algorithm is shown below:
wherein X (k) represents a fourier positive transform function with k as a variable; x (n) represents an inverse fourier transform function with n as a variable; n represents the number of sampling points;
The above equation of fourier positive transform is expressed in matrix form as:
the formula of the inverse fourier transform algorithm is expressed in the form of a matrix as:
Order the
Then
The formula for the fourier positive transform can be expressed as:
X＝Wx
the formula for the inverse fourier transform can be expressed as:
x＝W-1X
from the binding properties of the matrix (AB) c=a (BC), and using the fourier transform space derivative formula, it is possible to obtain:
First derivative of x: x' =w -1((ikx) Wx), first order derivative operator: w -1(ikx) W;
Second derivative of x: x "=w -1((ikx)2 Wx), second order derivative operator: w -1(ikx)2 W;
Third derivative of x: x' "=w -1((ikx)3 Wx), third order derivative operator: w -1(ikx)3 W;
Fourth derivative of x: x "" = W -1((ikx)4 Wx), fourth order derivative operator: w -1(ikx)4 W;
…
The m-order derivative operator in the x-axis direction is W -1(ikx)m W.
Wherein i represents an imaginary number, and N represents a sampling point number; k x is the wave number in the x-axis direction, which is the diagonal matrix of order N:
For a determined direction, such as the x-axis direction, the number of sampling points N is known, so that W, W -1 and k x can be determined according to the m-order derivative operator formula W -1(ikx)m W, thereby determining the x-axis direction second derivative operator.
Note that, the second order derivative operator W -1(ikx)2 W is a cyclic matrix, and N uniquely determines W, W -1 and k x, thereby uniquely determining one derivative operator from N.
The m-order derivative operator matrix in the y-axis direction is W -1(iky)n W; the m-order derivative operator matrix in the z-axis direction is W -1(ikz)n W.
Wherein, the matrix of wavenumbers k y in the y-axis direction is expressed as:
The matrix of wavenumbers k z in the z-axis direction is expressed as:
and 102, determining the spatial partial derivative of the three-dimensional wave equation according to the second-order derivative operator and the spatial wave field, and taking the spatial partial derivative as a forward result of the spatial domain of the three-dimensional wave equation.
And multiplying the second order derivative operator by the spatial wave field to obtain the spatial partial derivative of the three-dimensional wave equation.
In the embodiment of the application, an m-order derivative operator formula is obtained through derivation by a Fourier forward conversion algorithm and a Fourier inverse conversion algorithm, a second-order derivative operator in a specified coordinate direction can be directly determined through the m-order derivative operator formula, and then the spatial partial derivative of the three-dimensional wave equation can be determined according to the second-order derivative operator and the spatial wave field. Compared with the pseudo-spectrum method in the prior art, FFT calculation and IFFT calculation are not needed, and meanwhile, the calculation steps are reduced, so that the calculation amount of a computer is reduced, the calculation efficiency is improved, and the calculation of a single step is more suitable for GPU parallel calculation.
In the embodiment of the application, considering that the truncated of the derivative operator can generate a virtual source due to the periodicity of Fourier transform, in order to solve the problem of the virtual source, in the embodiment of the application, the derivative operators with different truncated lengths are weighted and averaged. As shown in fig. 2, the following steps 201 to 203 may be performed:
Step 201, calculating a second order derivative operator of the specified coordinate direction according to the sampling point number and the m order derivative operator formula of the specified coordinate direction aiming at the specified coordinate direction of the three-dimensional space.
Wherein, the m-order derivative operator formula is obtained by deduction according to a Fourier positive transformation algorithm and an inverse Fourier transformation algorithm, and m is more than or equal to 1.
The calculation method of the second order derivative operator in this step is the same as that in step 101.
And 202, performing truncation and zero filling processing on the second order derivative operator in the appointed coordinate direction to obtain initial second order derivative operators with different truncation lengths.
In the embodiment of the application, the derivative operator is analyzed by using a middle point derivative operator curve, and the middle point derivative operator curve of N=128 to 2048 is shown in fig. 3 (a) to 3 (e). From the derivative operator curve, the intermediate point derivative operator curve is oscillated and converged around the intermediate point.
As shown in fig. 4 (a) to 4 (f), the relationship between the intermediate point derivative operator curve and the spatial sampling point N is small, and the relationship between the intermediate point derivative operator curve and the intermediate point is mainly limited. Therefore, only the space derivative operator needs to be truncated by taking the space sampling point as the center, so that the operand is reduced, and the efficiency is improved.
In the embodiment of the application, for the second order derivative operator expressed in a matrix form, taking each row of peak values as a center, taking the designated number as a cut-off length, sequentially selecting the designated number of points to reserve left and right of the center, filling zero into other points except the reserved points, and calculating to obtain an initial second order derivative operator according to all reserved points; adding 1to the appointed number as a new cut-off length, and calculating again to obtain an initial second-order derivative operator; and repeating the process of adding 1to the specified number as a new truncated length to calculate the initial second order derivative operator until the value of the truncated length is equal to the number threshold value, and calculating to obtain the initial second order derivative operators with different truncated lengths.
Where the specified number is typically a positive power of 2, say 16, 32, 64 or 128, etc. Illustratively, the center of a row of the matrix is (A, A), the designated number is 16, the designated number of points sequentially selected on the left side of the center are (A, A-1), (A, A-2), (A, A-3), …, (A, A-16), the designated number of points sequentially selected on the right side of the center are (A, A+1), (A, A+2), (A, A+3), …, (A, A+16), then the initial second order derivative operator is calculated according to the above selected (A, A-16), …, (A, A), …, (A, A+16), and the other points such as (A, A-17), (A, A-18), (A+17) and the like are marked as 0. Where (a, a), (a, a-1), etc. are used to represent the location of a point, not to represent the actual numerical value of the point.
The number threshold may be set by the user, but considering that the larger the selected cut length is, the larger the calculated amount of the cut zero-filling process is, and the range of improvement of the calculation accuracy is limited, in order to balance the calculated amount and the calculation accuracy, the designated number used in the initial calculation is generally set to be 16, the number threshold is set to be 32, and the calculated amount of the cut zero-filling process is smaller and the accuracy can also meet the requirement in the range of the cut length.
And 203, carrying out weighted average processing on the initial second-order derivative operators with different cut-off lengths to obtain a second-order optimization derivative operator.
In the embodiment of the application, the same weight of the initial second-order derivative operators with different cut-off lengths can be set.
The points adopted in the calculation of each initial derivative algorithm comprise the appointed number of points adopted in the calculation of each initial calculation of the left and right of the center, other points except the appointed number are 0, for example, the appointed number adopted in the calculation of the initial calculation is 16, and the number threshold is 32, so that the center point and the 16 points of each calculation are used, the use frequency of the other points except the 16 points of the left and right is far smaller than that of the 16 points of the left and right, the center point and the points of each appointed number of the left and right are reinforced in the calculation of the initial derivative algorithm, the other points are weakened, the initial second derivative operators with different cut lengths are averaged, and the virtual source problem caused by the cut-off of the derivative operators is suppressed.
And 204, determining the spatial partial derivative of the three-dimensional wave equation according to the second-order optimization derivative operator and the spatial wave field, and taking the spatial partial derivative as a spatial domain forward result of the three-dimensional wave equation.
In the embodiment of the application, an m-order derivative operator formula is obtained through derivation by a Fourier forward conversion algorithm and a Fourier inverse conversion algorithm, a second-order derivative operator in a specified coordinate direction can be directly determined through the m-order derivative operator formula, and then the spatial partial derivative of the three-dimensional wave equation can be determined according to the second-order derivative operator and the spatial wave field. Compared with the pseudo-spectrum method in the prior art, FFT calculation and IFFT calculation are not needed, and meanwhile, the calculation steps are reduced, so that the calculation amount of a computer is reduced, the calculation efficiency is improved, and the calculation of a single step is more suitable for GPU parallel calculation.
In the embodiment of the application, the virtual seismic source problem is optimally suppressed by cutting off the two-dimensional acoustic wave forward operation verification operator. The model used was a four-layer horizontally uniform layered media model. The calculation area was 5120m×5120m. The first layer is 1000m thick and has a speed of 2000m/s; the thickness of the second layer is 1000m, and the speed is 3000m/s; the thickness of the third layer is 1000m, and the speed is 4000m/s; the fourth layer is 2120m thick at a speed of 5000m/s. Spatial step Δx=Δy=Δz=10m, time sampling 4000ms, time step 0.5ms. The source is located at (2560 m,0 m), the source main frequency is 30HZ. The single shot analysis of operator cut-off length 64 is shown in fig. 5 (a); the single shot analysis of operator cut-off length 128 is shown in fig. 5 (b); the single shot analysis of operator cut length 256 is shown in fig. 5 (c); the single shot analysis of the optimized truncation operator is shown in fig. 5 (d).
The three-dimensional acoustic wave equation forward model is a four-layer horizontal uniform layered medium model. The calculation area is 1280m×1280m. The first layer has a thickness of 200m and a speed of 2000m/s; the thickness of the second layer is 200m, and the speed is 3000m/s; the thickness of the third layer is 200m, and the speed is 4000m/s; the fourth layer is 680m thick at a speed of 5000m/s. Spatial step Δx=Δy=Δz=10m, time sampling 1000ms, time step 0.5ms. The source is located at (640 m,320m,0 m) and the source is at a main frequency of 30HZ. Fig. 6 is a velocity model, and fig. 7 is a three-dimensional forward single shot display.
In order to verify the above-mentioned beneficial effects of reducing the calculation amount, in the embodiment of the present application, taking N sampling points as an example, the calculation amounts of the pseudo-spectrum method in the prior art and the three-dimensional acoustic wave equation space domain forward modeling method (i.e. the improved algorithm) based on the pseudo-spectrum method in the embodiment of the present application are compared, and the comparison results are shown in the following table one and table two:
List one
Watch II
As can be seen from the first and second tables, compared with the prior art, the three-dimensional acoustic wave equation space domain forward modeling method based on the pseudo-spectrum method provided by the application has the advantages that the operation amount (comprising the times of addition operation and the times of multiplication operation) is greatly reduced, and correspondingly, the operation efficiency is also greatly improved.
The embodiment of the application also provides a three-dimensional acoustic wave equation space domain forward modeling device based on pseudo-spectroscopy, as shown in fig. 8, the device 800 comprises a calculation module 801 and a determination module 802.
The calculation module 801 is configured to calculate, for a specified coordinate direction of the three-dimensional space, a second-order derivative operator of the specified coordinate direction according to a sampling point number and an m-order derivative operator formula of the specified coordinate direction; the m-order derivative operator formula is obtained by deduction according to a Fourier positive transformation algorithm and an inverse Fourier transformation algorithm, and m is more than or equal to 1.
The determining module 802 is configured to determine a spatial partial derivative of the three-dimensional wave equation according to the second-order derivative operator and the spatial wave field obtained by the calculating module 801, and use the spatial partial derivative as a forward result of the spatial domain of the three-dimensional wave equation.
In one implementation of an embodiment of the present application, the apparatus 800 further includes:
the truncation processing module 803 is configured to perform truncation zero-filling processing on the second order derivative operator in the specified coordinate direction calculated by the calculation module 801, so as to obtain initial second order derivative operators with different truncation lengths.
The truncation processing module 803 is further configured to perform weighted average processing on the initial second order derivative operators with different truncation lengths, so as to obtain a second order optimized derivative operator.
The determining module 802 is further configured to determine a spatial partial derivative of the three-dimensional wave equation according to the second-order optimization derivative operator and the spatial wave field calculated by the truncation processing module 803, as a spatial domain forward result of the three-dimensional wave equation.
In one implementation of the embodiment of the present application, the truncation processing module 803 is configured to:
sequentially selecting a specified number of points to reserve by taking each row of peak values as a center and taking the specified number as a cut-off length, filling zero into other points except the reserved points, and calculating according to all reserved points to obtain an initial second order derivative operator; adding 1 to the appointed number as a new cut-off length, and calculating again to obtain an initial second-order derivative operator; and repeating the process of adding 1 to the specified number as a new truncated length to calculate the initial second order derivative operator until the value of the truncated length is equal to the number threshold value, and calculating to obtain the initial second order derivative operators with different truncated lengths.
In one implementation manner of the embodiment of the application, an m-order derivative operator formula of the x-axis direction of the three-dimensional space is W -1(ikx)m W; the m-order derivative operator matrix in the y-axis direction is W -1(iky)m W; the m-order derivative operator matrix in the z-axis direction is W -1(ikz)m W;
wherein k x、ky、kz represents wave numbers in the x-axis direction, the y-axis direction, and the z-axis direction, respectively; i represents an imaginary number; w represents the fourier transform matrix and, W -1 denotes the fourier transform inverse matrix,N represents the number of sampling points.
In the embodiment of the application, an m-order derivative operator formula is obtained through derivation by a Fourier forward conversion algorithm and a Fourier inverse conversion algorithm, a second-order derivative operator in a specified coordinate direction can be directly determined through the m-order derivative operator formula, and then the spatial partial derivative of the three-dimensional wave equation can be determined according to the second-order derivative operator and the spatial wave field. Compared with the pseudo-spectrum method in the prior art, FFT calculation and IFFT calculation are not needed, and meanwhile, the calculation steps are reduced, so that the calculation amount of a computer is reduced, the calculation efficiency is improved, and the calculation of a single step is more suitable for GPU parallel calculation.
The embodiment of the application also provides a computer device, which comprises a memory, a processor and a computer program stored on the memory and capable of running on the processor, wherein the processor executes the computer program to realize the method from step 101 to step 102.
The embodiment of the present application further provides a computer readable storage medium storing a computer program for executing the method of any one of steps 101 to 102.
It will be appreciated by those skilled in the art that embodiments of the present application may be provided as a method, system, or computer program product. Accordingly, the present application may take the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects. Furthermore, the present application may take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, and the like) having computer-usable program code embodied therein.
The present application is described 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 flow and/or block of the flowchart illustrations and/or block diagrams, and combinations of flows and/or blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.
These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instruction means which implement the function specified in the flowchart flow or flows and/or block diagram block or blocks.
These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.
The foregoing description of the embodiments has been provided for the purpose of illustrating the general principles of the application, and is not meant to limit the scope of the application, but to limit the application to the particular embodiments, and any modifications, equivalents, improvements, etc. that fall within the spirit and principles of the application are intended to be included within the scope of the application.
Claims (8)
1. A three-dimensional acoustic wave equation space domain forward modeling method based on a pseudo-spectrum method, which is characterized by comprising the following steps:
Aiming at a specified coordinate direction of a three-dimensional space, calculating a second-order derivative operator of the specified coordinate direction according to the sampling point number and an m-order derivative operator formula of the specified coordinate direction; the m-order derivative operator formula is obtained by deduction according to a Fourier positive transformation algorithm and an inverse Fourier transformation algorithm, and m is more than or equal to 1;
determining a spatial partial derivative of the three-dimensional wave equation according to the second order derivative operator and the spatial wave field, and taking the spatial partial derivative as a spatial domain forward result of the three-dimensional wave equation;
The formula of an m-order derivative operator in the x-axis direction of the three-dimensional space is W -1(ikx)m W; the m-order derivative operator matrix in the y-axis direction is W -1(iky)m W; the m-order derivative operator matrix in the z-axis direction is W -1(ikz)m W;
wherein k x、ky、kz represents wave numbers in the x-axis direction, the y-axis direction, and the z-axis direction, respectively; i represents an imaginary number; w represents the fourier transform matrix and, W -1 denotes the fourier transform inverse matrix,N represents the number of sampling points, W (N-1)*(N-1) J represents an imaginary number.
2. The method of claim 1, wherein after computing a second order derivative operator for the specified coordinate direction from the number of samples and the m-order derivative operator formula for the specified coordinate direction, the method further comprises:
performing truncation and zero filling processing on the second order derivative operator in the appointed coordinate direction to obtain initial second order derivative operators with different truncation lengths;
Carrying out weighted average treatment on the initial second order derivative operators with different cut-off lengths to obtain a second order optimization derivative operator;
And determining the spatial partial derivative of the three-dimensional wave equation according to the second-order optimization derivative operator and the spatial wave field, and taking the spatial partial derivative as a spatial domain forward result of the three-dimensional wave equation.
3. The method according to claim 2, wherein the performing the truncated zero-filling process on the second order derivative operator of the specified coordinate direction to obtain the initial second order derivative operators with different truncated lengths includes:
Sequentially selecting a specified number of points to reserve by taking each row of peak values as a center and taking the specified number as a cut-off length, filling zero into other points except the reserved points, and calculating according to all reserved points to obtain an initial second order derivative operator;
adding 1 to the appointed number as a new cut-off length, and calculating again to obtain an initial second-order derivative operator;
And repeating the process of adding 1 to the specified number as a new truncated length to calculate the initial second order derivative operator until the value of the truncated length is equal to the number threshold value, and calculating to obtain the initial second order derivative operators with different truncated lengths.
4. A pseudo-spectral method-based three-dimensional acoustic wave equation space domain forward modeling apparatus, the apparatus comprising:
the calculation module is used for calculating a second order derivative operator of the specified coordinate direction according to the sampling point number and the m order derivative operator formula of the specified coordinate direction aiming at the specified coordinate direction of the three-dimensional space; the m-order derivative operator formula is obtained by deduction according to a Fourier positive transformation algorithm and an inverse Fourier transformation algorithm, and m is more than or equal to 1;
The determining module is used for determining the spatial partial derivative of the three-dimensional wave equation according to the second-order derivative operator and the spatial wave field obtained by calculation of the calculating module and taking the spatial partial derivative as a forward result of the spatial domain of the three-dimensional wave equation;
The formula of an m-order derivative operator in the x-axis direction of the three-dimensional space is W -1(ikx)m W; the m-order derivative operator matrix in the y-axis direction is W -1(iky)m W; the m-order derivative operator matrix in the z-axis direction is W -1(ikz)m W;
wherein k x、ky、kz represents wave numbers in the x-axis direction, the y-axis direction, and the z-axis direction, respectively; i represents an imaginary number; w represents the fourier transform matrix and, W -1 denotes the fourier transform inverse matrix,N represents the number of sampling points, W (N-1)*(N-1) J represents an imaginary number.
5. The apparatus of claim 4, wherein the apparatus further comprises:
The truncation processing module is used for carrying out truncation zero-filling processing on the second-order derivative operator in the appointed coordinate direction calculated by the calculation module to obtain initial second-order derivative operators with different truncation lengths;
The truncation processing module is further used for carrying out weighted average processing on the initial second order derivative operators with different truncation lengths to obtain a second order optimization derivative operator;
the determining module is further used for determining the spatial partial derivative of the three-dimensional wave equation according to the second-order optimization derivative operator and the spatial wave field obtained by calculation of the truncation processing module, and the spatial partial derivative is used as a spatial domain forward result of the three-dimensional wave equation.
6. The apparatus of claim 5, wherein the truncation processing module is configured to:
Sequentially selecting a specified number of points to reserve by taking each row of peak values as a center and taking the specified number as a cut-off length, filling zero into other points except the reserved points, and calculating according to all reserved points to obtain an initial second order derivative operator;
adding 1 to the appointed number as a new cut-off length, and calculating again to obtain an initial second-order derivative operator;
And repeating the process of adding 1 to the specified number as a new truncated length to calculate the initial second order derivative operator until the value of the truncated length is equal to the number threshold value, and calculating to obtain the initial second order derivative operators with different truncated lengths.
7. A computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, characterized in that the processor implements the method of any of claims 1 to 3 when executing the computer program.
8. A computer readable storage medium, characterized in that the computer readable storage medium stores a computer program for executing the method of any one of claims 1 to 3.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911187548.0A CN112859162B (en) | 2019-11-28 | 2019-11-28 | Three-dimensional acoustic wave equation space domain forward modeling method and device based on pseudo-spectrum method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911187548.0A CN112859162B (en) | 2019-11-28 | 2019-11-28 | Three-dimensional acoustic wave equation space domain forward modeling method and device based on pseudo-spectrum method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112859162A CN112859162A (en) | 2021-05-28 |
| CN112859162B true CN112859162B (en) | 2024-04-30 |
Family
ID=75985142
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911187548.0A Active CN112859162B (en) | 2019-11-28 | 2019-11-28 | Three-dimensional acoustic wave equation space domain forward modeling method and device based on pseudo-spectrum method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112859162B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113656750B (en) * | 2021-10-20 | 2021-12-17 | 中南大学 | Magnetic induction intensity calculation method of strong magnetic medium based on space wave number domain |
| CN114509755A (en) * | 2022-01-17 | 2022-05-17 | 深圳力维智联技术有限公司 | Coal monitoring method, equipment and storage medium based on reverse time migration imaging algorithm |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103699798A (en) * | 2013-12-25 | 2014-04-02 | 中国石油天然气股份有限公司 | Method for realizing seismic wave field numerical simulation |
| CN106526677A (en) * | 2016-10-26 | 2017-03-22 | 中海石油(中国)有限公司 | Marine self-adaptive ghost reflection-suppressing broadband reverse time migration imaging method |
| CN106772585A (en) * | 2017-01-26 | 2017-05-31 | 中国科学院地质与地球物理研究所 | Analysis method and device is intended in a kind of optimization based on elastic wave decoupling equation |
| CN106842320A (en) * | 2017-01-19 | 2017-06-13 | 北京大学 | The parallel 3-D seismics wave field generation methods of GPU and system |
| CN107179549A (en) * | 2017-07-11 | 2017-09-19 | 中海石油(中国)有限公司 | A kind of acoustic wave equation in time domain Explicit finite difference seismic response analogy method |
| CN107561585A (en) * | 2017-09-19 | 2018-01-09 | 北京大学 | A kind of multinuclear multi-node parallel 3-D seismics wave field generation method and system |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9285491B2 (en) * | 2010-05-12 | 2016-03-15 | Shell Oil Company | Seismic P-wave modelling in an inhomogeneous transversely isotropic medium with a tilted symmetry axis |
| US8700372B2 (en) * | 2011-03-10 | 2014-04-15 | Schlumberger Technology Corporation | Method for 3-D gravity forward modeling and inversion in the wavenumber domain |
-
2019
- 2019-11-28 CN CN201911187548.0A patent/CN112859162B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103699798A (en) * | 2013-12-25 | 2014-04-02 | 中国石油天然气股份有限公司 | Method for realizing seismic wave field numerical simulation |
| CN106526677A (en) * | 2016-10-26 | 2017-03-22 | 中海石油(中国)有限公司 | Marine self-adaptive ghost reflection-suppressing broadband reverse time migration imaging method |
| CN106842320A (en) * | 2017-01-19 | 2017-06-13 | 北京大学 | The parallel 3-D seismics wave field generation methods of GPU and system |
| CN106772585A (en) * | 2017-01-26 | 2017-05-31 | 中国科学院地质与地球物理研究所 | Analysis method and device is intended in a kind of optimization based on elastic wave decoupling equation |
| CN107179549A (en) * | 2017-07-11 | 2017-09-19 | 中海石油(中国)有限公司 | A kind of acoustic wave equation in time domain Explicit finite difference seismic response analogy method |
| CN107561585A (en) * | 2017-09-19 | 2018-01-09 | 北京大学 | A kind of multinuclear multi-node parallel 3-D seismics wave field generation method and system |
Non-Patent Citations (6)
| Title |
|---|
| Parallel 3-D pseudospectral simulation of seismic wave propagation;Furumura T 等;Geophysics;第63卷(第01期);第279-288页 * |
| 伪谱法地震波正演模拟的多线程并行计算;谢桂生 等;地球物理学进展;第20卷(第01期);第17-23页 * |
| 伪谱法运算效率的提高及井间波场的正演模拟;田仁飞 等;物探与化探;第32卷(第04期);第430-433、441页 * |
| 基于PML边界的时间高阶伪谱法弹性波场模拟;李青阳 等;地球物理学进展;第33卷(第01期);第228-235页 * |
| 基于伪谱法的VSP逆时偏移及其应用研究;孙文博 孙赞东;地球物理学报;第53卷(第09期);第2196-2203页 * |
| 基于伪谱法的三维声波方程正演算法优化;闫海洋 等;第四届油气地球物理学术年会论文集;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112859162A (en) | 2021-05-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112859162B (en) | Three-dimensional acoustic wave equation space domain forward modeling method and device based on pseudo-spectrum method | |
| RU2012154894A (en) | METHOD AND SYSTEM FOR CREATING CONTROL POINTS DURING MODELING | |
| CN106168942B (en) | A kind of fluctuation types dynamic data reconstructing method based on singular boundary method | |
| CN111596347B (en) | Method and device for rapidly obtaining surface layer longitudinal and transverse wave speed ratio | |
| CN109506135A (en) | Pipe leakage independent positioning method and device | |
| Ohtani et al. | Fast computation of quasi-dynamic earthquake cycle simulation with hierarchical matrices | |
| US9564115B2 (en) | Producing sounds in a virtual world and related sound card | |
| CN107561586B (en) | A kind of method and apparatus of bubble compacting | |
| CN108762719A (en) | A kind of parallel broad sense inner product reconfigurable controller | |
| CN105068119A (en) | Attenuation method and apparatus for surface wave in low-frequency earthquake data | |
| Petrov et al. | Higher-order grid-characteristic schemes for the acoustic system | |
| CN107220212B (en) | Boundary element calculation method for two-dimensional non-compact boundary acoustic scattering | |
| CN113866827B (en) | Interpretation velocity modeling seismic imaging method, system, medium and equipment | |
| CN112419493B (en) | Shale reservoir three-dimensional attribute model building method and device | |
| Norwood | Transient response of an elastic plate to loads with finite characteristic dimensions | |
| Wang et al. | Real-Time CUDA Based Collision detection and Physics based collision response simulation | |
| Lin et al. | A time-domain coupled scaled boundary isogeometric approach for earthquake response analysis of dam-reservoir-foundation system | |
| CN103440228B (en) | A kind of method for accelerating FFT to calculate based on the multiply-add instruction of fusion | |
| RU2611892C1 (en) | Method of three-dimensional simulation of specified hydrogeological feature implemented in computer system | |
| Zhang et al. | Multiple Prediction through Inversion-GPU Acceleration | |
| CN110516341B (en) | Noise reduction method for additional damping of gearbox based on modal strain energy | |
| CN117724798A (en) | Parameter determination method, device and equipment of three-dimensional finite difference forward modeling method | |
| CN113935119A (en) | Sound source calculation method, calculation device, and computer-readable storage medium | |
| Wu et al. | A holistic approach for prediction of chatter stability using piecewise Hermite interpolation | |
| Degtyarev et al. | Virtual testbed: ship motion simulation for personal workstations |
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 | ||
| GR01 | Patent grant |