CN111208563A - Non-split complete matching layer absorption boundary method - Google Patents

Non-split complete matching layer absorption boundary method Download PDF

Info

Publication number
CN111208563A
CN111208563A CN202010099497.2A CN202010099497A CN111208563A CN 111208563 A CN111208563 A CN 111208563A CN 202010099497 A CN202010099497 A CN 202010099497A CN 111208563 A CN111208563 A CN 111208563A
Authority
CN
China
Prior art keywords
residual
boundary
equation
matching layer
wave
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.)
Granted
Application number
CN202010099497.2A
Other languages
Chinese (zh)
Other versions
CN111208563B (en
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.)
Jilin University
Original Assignee
Jilin University
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 Jilin University filed Critical Jilin University
Priority to CN202010099497.2A priority Critical patent/CN111208563B/en
Publication of CN111208563A publication Critical patent/CN111208563A/en
Application granted granted Critical
Publication of CN111208563B publication Critical patent/CN111208563B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

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/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

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)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention relates to a non-split complete matching layer absorption boundary method, wherein a lower medium is semi-infinite, and a calculation space is limited, so that an absorption boundary needs to be introduced when seismic wave propagation is simulated. The absorption boundary does not exist in the underground medium really, the actual seismic waves propagate to the absorption boundary without forming reflection, and the artificial boundary can bring about certain false reflection. The residual completely-matched layer provided by the invention can effectively absorb incident waves and prevent seismic waves from being transmitted back to the main area. The method does not change the form of the original equation, is beneficial to programming realization, and can be popularized to more complex medium simulation; the auxiliary equation eliminates the time partial derivative term of the original variable, and is more favorable for being popularized to high-precision forward modeling. Meanwhile, complex frequency shift conversion is successfully realized on the basis of a residual completely-matched layer, so that the complex frequency shift conversion can absorb the grazing incoming waves and avoid the generation of low-frequency singular values; and a complete matching layer structure with double attenuation profiles is adopted, so that the stability of the double-attenuation profile is improved.

Description

Non-split complete matching layer absorption boundary method
Technical Field
The invention belongs to the technical field of seismic exploration, particularly relates to an absorption boundary of a time domain, and particularly relates to a non-split complete matching layer absorption boundary method for realizing wave field transformation by solving a residual error.
Background
In order to meet the requirement of seismic wave forward inversion, the research of boundary conditions draws great attention. Currently, the most widely used boundary condition is the Perfect Matching Layers (PMLs) proposed by berenser in studying maxwell's equations. Which is interpreted by Chew et al and Collino et al as a result of a complex coordinate stretch transformation. When the incident wave angle is large, the attenuation coefficient becomes small, so that the conventional perfect matching layer cannot absorb the large-angle incident wave and generates a singular value at a very low frequency. Complex Frequency Shifted (CFS) techniques further process the Complex coordinate stretch transform (CCS), the term originally 1 is replaced by a scale factor and plays a role in the absorption of large angle incident waves, and a frequency shifting factor is added to the denominator to eliminate low frequency singular values. Roden et al propose convolving the perfect matching layer (C-PML) based on CFS; komatitsch et al further generalize it to elastic media.
In addition, the perfect matching layer is obtained by transforming the whole equation and directly transforming the wave field. For example, the near perfect matching layer npml (near pml). The method is rapidly popularized and tried in sound wave media, elastic media and two-phase media after being proposed. However, the method is weak in absorption capacity at large incidence angles and unstable in simulation of complex media, and meanwhile, the transformation equation of the method has time partial derivatives at both ends of the equation, so that the method is not beneficial to being popularized to simulation with high time, and therefore, the mainstream method in seismic forward inversion is still C-PML. For stability of boundary conditions, Festa et al indicate that absorption between the PML boundary and the surface wave is reduced and instability is created. Meza-Fajardo et al believe that conventional PMLs do not meet strict progressive stability in both isotropy and anisotropy, and introduce multiaxial techniques, i.e., the simultaneous introduction of attenuation factors in multiple orthogonal directions.
In short, current perfect match layers can be classified into split and non-split methods, wavefield transforms and equation transforms. The splitting method occupies a large memory and has low calculation efficiency; compared with equation transformation, wave field transformation has a form without changing a wave equation, and is more favorable for realizing and popularizing boundaries, the existing method for wave field transformation is an approximate complete matching layer, the problem of stability and the problem of difficulty in absorbing waves incident at a large angle exist, and meanwhile, original variables and new variables of the transformation equation are in a time partial derivative form, and the method is not favorable for popularizing in seismic wave simulation with higher precision.
Disclosure of Invention
The present invention aims to overcome the defects of the prior art, and provides a non-split complete matching layer absorption boundary method for implementing wave field transformation by solving a residual error, wherein the boundary condition has a better absorption effect and a simpler auxiliary equation form.
The purpose of the invention is realized by the following technical scheme:
the core of the invention is that in the process of solving a new equation after complex frequency shift or complex stretch transformation, a target format of a new equation is assumed when wave fields are transformed, a residual error item is introduced, then the residual error is solved by reverse calculation according to the target format, and a time partial derivative item of an original variable is eliminated in the process, so that the form of an auxiliary equation is simpler. According to the invention, through reverse derivation, the time partial derivative term of an original variable is further eliminated on the basis of non-splitting and wave field transformation, so that the equation is simpler, the method can be popularized to time high-order simulation, and meanwhile, the complex frequency shift transformation is successfully realized on the method in a mode that a scale factor directly acts on a wave field, so that the method can absorb the grazing wave and avoid the generation of low-frequency singular values; and meanwhile, double attenuation profiles are introduced, so that the stability of the composite material is further improved.
The invention discloses a non-split complete matching layer absorption boundary method, which comprises the following steps:
A. carrying out Fourier transform on the wave equation, converting the wave equation into a frequency domain, converting the wave equation into a time domain after carrying out complex stretch transform conversion, and obtaining a new wave equation after introducing a residual complete matching layer;
B. introducing a residual perfect matching layer for forward modeling:
a. expanding a certain number of layers of the calculation variable matrix to the periphery to serve as a complete matching layer for absorbing seismic waves, wherein a main area is a calculation interval and is also a simulated target area, and an expanded range is a PML area;
b. setting an attenuation factor value α, dividing a residual perfect matching layer into three regions, namely an upper boundary and a lower boundary which take the z direction as a main direction, a left boundary and a right boundary which take the x direction as a main direction and an angle interval;
α(x)＝K[φ(x/L)n+γe(-δL/x)]
Figure BDA0002386466770000031
in the above formula, phi (x/L)nThe function meets the basic requirement of the attenuation function, fine adjustment is carried out on the function by gamma exp (-delta L/x), wherein K, gamma, phi and delta are adjustable coefficients, the value of K in the experiment is given, n is the order, and L is the layer thickness of PML; r is the theoretical boundary reflection coefficient; l is the distance between the PML area calculation point and the internal boundary;
c. assigning values to the frequency shift factor and the scale factor in the PML area, wherein the value of the frequency shift factor in the main area is 0, and the value of the scale factor is 1;
Figure BDA0002386466770000032
d. converting a wave equation under a given medium into a frequency domain according to a formula 1, performing complex frequency shift or complex stretch conversion, and converting the wave equation into a time domain;
e. discretizing the whole new wave equation and a residual solving formula, and during iteration, iteratively solving a residual item, subtracting an original variable, and then iterating a new variable to obtain seismic wave simulation capable of rapidly attenuating seismic wave energy in a PML region; and (4) iterating variable and residual items, and only processing the variable and the residual in the PML area.
Further, step a, the process of introducing the residual perfect matching layer is as follows:
Figure BDA0002386466770000033
in the formula: ω is the circle frequency; λ, μ are Lame constants; ρ is the density; vx、VzIs the frequency domain elastic wave field velocity component; t isxx、Tzz、TxzIs a frequency domain elastic wave field stress component;
takes the form of a complex coordinate stretch transform (CCS)
Figure RE-GDA0002443533020000041
Substituting the formula (3) into the formula (1) and directly acting on the velocity component and the stress component to obtain an equation set after complex coordinate stretching transformation
Figure BDA0002386466770000042
The form of the wave equation after the introduction of the residual perfect match layer is assumed to be that of the subtraction of the wavefield and the residual, in which form the residual epsilon is inversely derived.
Further, the residual error epsilon is obtained as follows:
order to
Figure BDA0002386466770000043
Figure BDA0002386466770000044
Figure BDA0002386466770000045
The remaining residual terms can be solved, and the result is transformed into the time domain through inverse Fourier transform:
Figure BDA0002386466770000051
and is provided with
Figure BDA0002386466770000052
In the formula tauxx、τzx、τzzIs the stress, velocity parameter vx、vzAnd T in the frequency domainxx、Tzx、TzzAnd Vx、VzCorresponding; the corresponding residual error is solved through the equation (9) and is brought into the original equation.
Further, step A, realizing complex frequency shift conversion based on the theory of residual perfect matching layer
Figure BDA0002386466770000053
In the formula ηx(x) Is a frequency shift factor βx(x) Is a scale factor;
further, the residual under the transform is obtained as follows:
Figure BDA0002386466770000054
Figure BDA0002386466770000055
Figure BDA0002386466770000056
the remaining residuals are calculated to be transformed into the time domain:
Figure BDA0002386466770000061
is provided with
Figure BDA0002386466770000062
Further, in step a, the number of layers is 10 to 20.
Compared with the prior art, the invention has the beneficial effects that:
the residual perfect matching layer provided by the invention eliminates the time partial derivative term of the original variable by using a reverse derivation method, is a wave field transformation non-splitting perfect matching layer with a simpler auxiliary equation, simplifies the realization of the whole matching layer, enables the whole matching layer to be popularized to a high-precision forward modeling, improves the absorption capacity of the whole matching layer on the grazing incoming wave, and enhances the stability. The method has the following advantages:
1. the method is a non-splitting method, the storage space occupation is small, the calculation efficiency is high, only one residual error item is subtracted on the basis of an original variable because the form of an original equation is not changed, so that a main code is not changed during programming, the programming realization is facilitated, and meanwhile, the result after the residual error complete matching layer is introduced is equivalent to only adding or subtracting the variable, the process of obtaining the residual error item is independent and only related to the original variable, so that the whole process of adding the residual error complete matching layer is simple, and the method is more easily popularized to seismic wave simulation in other media.
2. The method eliminates the time partial derivative term of the original variable in the reverse solving process, so that the auxiliary equation is simpler, compared with other wave field transformation methods, only one time partial derivative term is required, the time partial derivative term is consistent with the original equation set form, and the method can be popularized to time high-order finite difference simulation.
3. The perfect matching layer cannot absorb large-angle incident waves, because the waves propagate along a direction approximately parallel to the boundary, and in order to further improve the absorption capacity of the residual perfect matching layer, complex frequency shift transformation is applied on the basis, and frequency shift factors and scale factors are introduced. The frequency shift factor can avoid the generation of low-frequency singular values, and the scale factor can bend the incident wave to the normal direction, so that the wave can be transmitted to a deeper position of the boundary, the absorption of the grazing incoming wave is enhanced, and the absorption capacity of the residual completely-matched layer is improved.
4. The method adopts the attenuation factor distribution of double attenuation sections in the boundary, namely, a damping section is introduced in the direction parallel to the boundary (as shown in figure 1, the upper boundary in figure 1 is a single section, and the right boundary shows a double section), seismic waves incident into the boundary are absorbed from two directions, the residual energy in the boundary is further eliminated, and the residual energy in the boundary can cause the instability of the boundary and even form false reflections because the energy is excessively accumulated and is reversely transmitted back to the boundary, so that the method has better stability than the rest complete matching layers.
5. The double attenuation sections can prevent waves from entering a boundary to a certain degree under the condition of large-angle incidence, false reflection can be enhanced, and the reflection is enhanced when the stability factor value is larger. However, due to the introduction of the frequency shift factor, the boundary absorption is strengthened, so that the adverse effect is counteracted.
The method realizes the non-split complete matching layer absorption boundary of the wave field transformation by solving the residual error, has high calculation efficiency, simple realization, strong absorption effect and high stability, is beneficial to popularization, keeps consistent with the wave equation form and can be popularized to the finite difference simulation of the time high order.
Drawings
FIG. 1 is a schematic diagram of a residual perfect matching layer architecture;
FIG. 2 residual perfect match layer wavefield snapshots at different times (up: velocity x component, down: velocity z component);
FIG. 3a wave field snapshot at 0.25s for the velocity x component (different residual perfect matching layers, M-DPML: only using dual profiles, P being a stability factor, MC-DPML: simultaneously using dual profiles and complex frequency shift transform, with a maximum value of frequency shift factor of 3 and a maximum value of scale factor of 2);
FIG. 3b wave field snapshot at 0.3s for the velocity x component (different residual perfect matching layers, M-DPML: only using dual profile, P being stability factor, MC-DPML: simultaneously using dual profile and complex frequency shift transform, maximum value of frequency shift factor is 3, maximum value of scale factor is 2);
FIG. 4 residual perfect match layer wavefield snapshots at different times (first three: velocity x component, last three: velocity z component);
FIG. 5a wave field snapshot at 0.2s for the velocity x component (different residual perfect matching layers, M-DPML: only using dual profiles, P being a stability factor, MC-DPML: simultaneously using dual profiles and complex frequency shift transform, with a maximum value of frequency shift factor taken as 3 and a maximum value of scale factor taken as 2);
FIG. 5b wave field snapshot at 0.5s for the velocity x component (different residual perfect match layers, M-DPML: only using dual profile, P is stability factor, MC-DPML: using both dual profile and complex frequency shift transform, max frequency shift factor taken as 3, max scale factor taken as 2).
Detailed Description
The invention is described in further detail below with reference to the following figures and examples:
the non-split perfect match layer absorption boundary method of the present invention is implemented on the MATLAB2014b platform.
The invention discloses a non-split complete matching layer absorption boundary method, which comprises the following steps:
A. and Fourier transformation is carried out on the wave equation, the wave equation is converted into a frequency domain, and then the wave equation is converted into a time domain after complex stretch transformation, so that a new wave equation with a residual completely matched layer introduced is obtained. The process of introducing the residual perfect matching layer is as follows:
Figure BDA0002386466770000081
in the formula: ω is the circle frequency; λ, μ are Lame constants; ρ is the density; vx、VzIs the frequency domain elastic wave field velocity component; t isxx、Tzz、TxzIs the frequency domain elastic wave field stress component.
Takes the form of a complex coordinate stretch transform (CCS)
Figure RE-GDA0002443533020000082
Figure RE-GDA0002443533020000091
Substituting the formula (3) into the formula (1) and directly acting on the velocity component and the stress component to obtain an equation set after complex coordinate stretching transformation
Figure BDA0002386466770000092
The form of the wave equation after the introduction of the residual perfect match layer is assumed to be that of the subtraction of the wavefield and the residual, in which form the residual epsilon is inversely derived. The residual epsilon is obtained as follows:
order to
Figure BDA0002386466770000093
Figure BDA0002386466770000094
Figure BDA0002386466770000095
The same way can get the rest residual terms, transform the result into time domain through inverse fourier transform:
Figure BDA0002386466770000101
and is provided with
Figure BDA0002386466770000102
In the formula tauxx、τzx、τzzIs the stress, velocity parameter vx、vzAnd T in the frequency domainxx、Tzx、TzzAnd Vx、VzAnd correspondingly. And solving the corresponding residual error through the equation 9, and bringing the residual error into the original equation. Then, the complex frequency shift transformation is realized on the theoretical basis of the residual perfect matching layer.
Figure BDA0002386466770000103
In the formula ηx(x) Is a frequency shift factor βx(x) Is a scale factor. The residual under this transformation is found as follows:
the time partial derivative item which can generate an original variable is directly obtained through the method, so that the residual error equation has two time flat partial derivative items which are not beneficial to high-precision forward modeling. Therefore, the following changes are needed when defining the target form:
Figure BDA0002386466770000104
Figure BDA0002386466770000105
Figure BDA0002386466770000106
the remaining residuals are calculated to be transformed into the time domain:
Figure BDA0002386466770000111
is provided with
Figure BDA0002386466770000112
B. The whole process of introducing the residual perfect matching layer and performing forward modeling comprises the following steps:
a. expanding a certain number of layers of the calculation variable matrix to the periphery to serve as complete matching layers for absorbing seismic waves, wherein the number of the layers is generally set to be 10-20, as shown in fig. 1, a main area is a calculation interval and is also a simulated target area, and the expanded range is a PML area;
b. the method comprises the steps of setting attenuation factor values α, dividing a residual complete matching layer into three regions, namely an upper boundary and a lower boundary which take a z direction as a main direction, a left boundary and a right boundary which take an x direction as a main direction, and an angle section, wherein the upper boundary shows a single section and the right boundary shows a double section to know the difference between the single attenuation section and the double attenuation section as shown in figure 1.
α(x)＝K[φ(x/L)n+γe(-δL/x)]
Figure BDA0002386466770000113
The above formula is divided into two parts, phi (x/L)nThe function is made to meet the basic requirements of the attenuation function, and fine adjustment is carried out on the function by gamma exp (-delta L/x). In the formula, K, gamma, phi and delta are all adjustable coefficients, the value of K in the experiment is given, and n is an order. L is the layer thickness of the PML; r is the theoretical boundary reflection coefficient; l is the distance between the PML area calculation point and the internal boundary;
c. assigning frequency shift factors and scale factors in the PML area, wherein the assignment areas are also in the PML area of FIG. 1, and the frequency shift factors and the scale factors are functions related to x and z, but the frequency shift factors are 0 and the scale factors are 1 in the main area;
Figure BDA0002386466770000121
d. the wave equation under a given medium is converted into a frequency domain according to the formula 1, and then is converted into a time domain after complex frequency shift or complex stretch conversion is carried out. For different wave equations, the form of the obtained new equation after the residual complete matching layer is introduced is unchanged, and the same residual solving formula is also the same;
e. discretizing the whole new wave equation and a residual solving formula, and during iteration, iteratively solving a residual item, subtracting an original variable, and then iterating a new variable to obtain seismic wave simulation capable of rapidly attenuating seismic wave energy in a PML region;
when the variables and the residual items are iterated, only the variables and the residual in the PML area can be processed, and the residual in the main area is always 0, so that the processing is not needed, and the storage space is saved and the calculation efficiency is improved.
Example 1
According to the exploration requirement, constructTwo seismic source profiles are immediately common. Various boundary conditions and attenuation functions are applied to the seismic wave numerical simulation, the selected time step is 0.5ms, the loaded seismic source is a Rake wavelet, and the dominant frequency is 25 HZ. And establishing a model and researching the condition of seismic waves under normal incidence. The model size is 1100 × 1100m (including the boundary thickness), and the grid step size x direction and z direction are both 5 m. The number of PML layers is set to 10, and R is 10-6γ is 0.001, Φ and δ are 1 and 0.029, respectively, n is 2, and the seismic source is loaded at x 550m and z 550 m. The longitudinal wave speed is 3000m/s, the transverse wave speed is 1400m/s, and the density is 2000Kg/m 3.
The forward modeling results are shown in fig. 2-3 b, and wavefield simulation snapshots at different times are given in fig. 2, and the propagation process of seismic waves in the subsurface is recorded. The seismic waves are generated by vibration at a seismic source and are transmitted in a target area, the seismic waves are absorbed when reaching a complete matching layer of residual errors, the actual underground medium does not have a boundary, but the seismic waves are continuously transmitted to the boundary without reflection, and the absorption boundary is adopted to absorb the waves transmitted to the boundary of the target area in order to simulate the actual situation. In fig. 2, the seismic wave is reflected, and for the convenience of research and observation of the performance of different residual perfect matching layers, the threshold is reduced when the wave field snapshot is made, so that the reflected wave is displayed on the snapshot. In fig. 3, the double-section has little effect on the seismic waves under the normal incidence, and because the complex frequency shift transformation introduces a frequency shift factor and a scale factor, the absorption capacity is enhanced, and the false reflections under the same threshold value are not shown.
Example 2
Conventional incidence was simulated in example 1, and the case of near-parallel incidence was now designed. The selected time step is 0.5ms, the loading seismic source is a Rake wavelet, and the dominant frequency is 25 HZ. And establishing a model and researching the condition that the seismic waves are incident at a large angle. The model size is 2600 × 1100m (including the boundary thickness), and the grid step size is 5m in both the x-direction and the z-direction. The number of PML layers is set to 10, and R is 10-6γ is 0.001, Φ and δ are 1 and 0.029, respectively, n is 2, and the seismic source is loaded at x 1050m and z 55 m. The longitudinal wave speed is 3000m/s, the transverse wave speed is 1400m/s, and the density is 2000Kg/m 3.
The simulation results are shown in fig. 4 and 5. When incident near parallel, the residual perfect matching layer creates energy accumulation at the upper boundary, causing spurious reflections and instability. As shown in fig. 5a, when a double attenuation profile is introduced (multi-axis residual perfect matching layer M-DPML), the residual energy in the upper boundary is attenuated, the larger the attenuation factor, the more the attenuation, but the spurious reflections are enhanced. This is because the additionally introduced attenuation profile hinders the entry of the wave to some extent. And the multi-axis complex frequency shift residual complete matching layer (MC-DPML) can further eliminate residual energy in the boundary, thereby avoiding the situation that a larger stability factor is adopted to ensure that the energy in the boundary under a complex condition is not accumulated to cause the instability of the system. The above two examples demonstrate that the residual perfect matching layer can effectively absorb the incident wave, and here the reflected wave is amplified to be visible in the wave field snapshot only for the purpose of studying the effect of further work on this basis. Further work on DPML has also been demonstrated, with the final inventive result being better absorption and stability of MC-DPML.
The method is based on a non-splitting method for directly transforming the wave field to derive a new complete matching layer, the complete matching layer has the characteristic of not changing the form of an original equation set, the method is convenient to popularize into more complex medium simulation, and meanwhile, the method is in a non-splitting form and has good calculation performance. In order to enable a residual completely-matched layer to absorb a grazing incoming wave, improve the overall absorption capacity and avoid the generation of low-frequency singular values, complex frequency shift transformation is further adopted on the basis; in order to ensure that the whole system can stably and effectively operate, the attenuation function distribution of the double attenuation profiles is adopted, and the stability of the residual complete matching layer is improved.

Claims (6)

1. A non-split perfect match layer absorption boundary method, comprising the steps of:
A. carrying out Fourier transform on the wave equation, converting the wave equation into a frequency domain, converting the wave equation into a time domain after carrying out complex stretch transform conversion, and obtaining a new wave equation after introducing a residual complete matching layer;
B. introducing a residual perfect matching layer for forward modeling:
a. expanding a certain number of layers of the calculation variable matrix to the periphery to serve as a complete matching layer for absorbing seismic waves, wherein a main area is a calculation interval and is also a simulated target area, and an expanded range is a PML area;
b. setting an attenuation factor value α, dividing a residual perfect matching layer into three regions, namely an upper boundary and a lower boundary which take the z direction as a main direction, a left boundary and a right boundary which take the x direction as a main direction and an angle interval;
α(x)＝K[φ(x/L)n+γe(-δL/x)]
Figure FDA0002386466760000011
in the above formula, phi (x/L)nThe function meets the basic requirement of the attenuation function, fine adjustment is carried out on the function by gamma exp (-delta L/x), wherein K, gamma, phi and delta are adjustable coefficients, the value of K in the experiment is given, n is the order, and L is the layer thickness of PML; r is the theoretical boundary reflection coefficient; l is the distance between the PML area calculation point and the internal boundary;
c. assigning values to the frequency shift factor and the scale factor in the PML area, wherein the value of the frequency shift factor in the main area is 0, and the value of the scale factor is 1;
Figure FDA0002386466760000012
d. converting a wave equation under a given medium into a frequency domain according to a formula 1, performing complex frequency shift or complex stretching conversion, and converting the wave equation into a time domain;
e. discretizing the whole new wave equation and a residual solving formula, and during iteration, iteratively solving a residual item, subtracting an original variable, and then iterating a new variable to obtain seismic wave simulation capable of rapidly attenuating seismic wave energy in a PML region; and (4) iterating variable and residual items, and only processing the variable and the residual in the PML area.
2. The non-split perfect match layer absorption boundary method of claim 1, wherein: step A, the process of introducing the residual perfect matching layer is as follows:
Figure RE-FDA0002443533010000021
in the formula: ω is the circle frequency; λ, μ are Lame constants; ρ is the density; vx、VzIs the frequency domain elastic wave field velocity component; t isxx、Tzz、TxzIs a frequency domain elastic wave field stress component;
takes the form of a complex coordinate stretch transform (CCS)
Figure RE-FDA0002443533010000022
Figure RE-FDA0002443533010000023
Substituting the formula (3) into the formula (1) and directly acting on the velocity component and the stress component to obtain an equation set after complex coordinate stretching transformation
Figure RE-FDA0002443533010000024
The form of the wave equation after the introduction of the residual perfect match layer is assumed to be that of the subtraction of the wavefield and the residual, in which form the residual epsilon is inversely derived.
3. A non-split perfect match layer absorption boundary method according to claim 2, characterized by: the residual error epsilon is obtained through the following process:
order to
Figure FDA0002386466760000031
Figure FDA0002386466760000032
Figure FDA0002386466760000033
The remaining residual terms can be solved, and the result is transformed into the time domain through inverse Fourier transform:
Figure FDA0002386466760000034
and is provided with
Figure FDA0002386466760000035
In the formula tauxx、τzx、τzzIs the stress, velocity parameter vx、vzAnd T in the frequency domainxx、Tzx、TzzAnd Vx、VzCorresponding; the corresponding residual error is solved through the equation (9) and is brought into the original equation.
4. The non-split perfect match layer absorption boundary method of claim 1, wherein: step A, realizing complex frequency shift conversion on the theoretical basis of a residual perfect matching layer
Figure FDA0002386466760000036
In the formula ηx(x) Is a frequency shift factor βx(x) Is a scale factor.
5. The non-split perfect match layer absorption boundary method of claim 4, wherein: the residual under this transformation is found as follows:
Figure FDA0002386466760000037
Figure FDA0002386466760000041
Figure FDA0002386466760000042
the remaining residuals are calculated to be transformed into the time domain:
Figure FDA0002386466760000043
is provided with
Figure FDA0002386466760000044
6. The non-split perfect match layer absorption boundary method of claim 1, wherein: step a, the number of layers is 10 to 20.
CN202010099497.2A 2020-02-18 2020-02-18 Non-split complete matching layer absorption boundary method Expired - Fee Related CN111208563B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010099497.2A CN111208563B (en) 2020-02-18 2020-02-18 Non-split complete matching layer absorption boundary method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010099497.2A CN111208563B (en) 2020-02-18 2020-02-18 Non-split complete matching layer absorption boundary method

Publications (2)

Publication Number Publication Date
CN111208563A true CN111208563A (en) 2020-05-29
CN111208563B CN111208563B (en) 2021-08-06

Family

ID=70782136

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010099497.2A Expired - Fee Related CN111208563B (en) 2020-02-18 2020-02-18 Non-split complete matching layer absorption boundary method

Country Status (1)

Country Link
CN (1) CN111208563B (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112285788A (en) * 2020-10-29 2021-01-29 吉林大学 CPML (continuous phase markup language) absorption boundary condition loading method based on electromagnetic wave equation
CN112904417A (en) * 2021-01-21 2021-06-04 中国石油大学(华东) Finite difference simulation method for seismic wave propagation of prepressing solid medium
CN113917526A (en) * 2020-07-10 2022-01-11 中国石油化工股份有限公司 Forward modeling method based on non-split complete matching layer absorption boundary

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102880590A (en) * 2012-09-25 2013-01-16 复旦大学 Method for constructing non-split complete matching layer of second-order fluctuation equation
CN104749628A (en) * 2015-03-30 2015-07-01 西安交通大学 Absorbing boundary reflection method based on dispersal viscosity wave equation
CN106353797A (en) * 2015-07-17 2017-01-25 中国石油化工股份有限公司 High-precision earthquake forward modeling method
EP3256886B1 (en) * 2015-02-13 2020-01-29 Exxonmobil Upstream Research Company Efficient and stable absorbing boundary condition in finite-difference calculations

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102880590A (en) * 2012-09-25 2013-01-16 复旦大学 Method for constructing non-split complete matching layer of second-order fluctuation equation
EP3256886B1 (en) * 2015-02-13 2020-01-29 Exxonmobil Upstream Research Company Efficient and stable absorbing boundary condition in finite-difference calculations
CN104749628A (en) * 2015-03-30 2015-07-01 西安交通大学 Absorbing boundary reflection method based on dispersal viscosity wave equation
CN106353797A (en) * 2015-07-17 2017-01-25 中国石油化工股份有限公司 High-precision earthquake forward modeling method

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
RENÉ MATZEN: "An efficient finite element time-domain formulation for the elastic second-order wave equation: A non-split complex frequency shifted convolutional PML", 《INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING》 *
罗玉钦等: "多轴复频移近似完全匹配层弹性波模拟", 《石油地球物理勘探》 *
罗玉钦等: "近似完全匹配层边界条件吸收效果分析及衰减函数的改进", 《石油地球物理勘探》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113917526A (en) * 2020-07-10 2022-01-11 中国石油化工股份有限公司 Forward modeling method based on non-split complete matching layer absorption boundary
CN112285788A (en) * 2020-10-29 2021-01-29 吉林大学 CPML (continuous phase markup language) absorption boundary condition loading method based on electromagnetic wave equation
CN112285788B (en) * 2020-10-29 2021-09-14 吉林大学 CPML (continuous phase markup language) absorption boundary condition loading method based on electromagnetic wave equation
CN112904417A (en) * 2021-01-21 2021-06-04 中国石油大学(华东) Finite difference simulation method for seismic wave propagation of prepressing solid medium

Also Published As

Publication number Publication date
CN111208563B (en) 2021-08-06

Similar Documents

Publication Publication Date Title
CN111208563B (en) Non-split complete matching layer absorption boundary method
Alkhalifah Scanning anisotropy parameters in complex media
Carcione Theory and modeling of constant-Q P-and S-waves using fractional time derivatives
Yan et al. Improving the efficiency of elastic wave-mode separation for heterogeneous tilted transverse isotropic media
CN110058307B (en) Full waveform inversion method based on fast quasi-Newton method
CN109633752B (en) Offshore towing cable data self-adaptive ghost wave compression method based on three-dimensional fast Radon transformation
CN111580163B (en) Full waveform inversion method and system based on non-monotonic search technology
Ma et al. Nonsplit complex-frequency shifted perfectly matched layer combined with symplectic methods for solving second-order seismic wave equations—Part 1: Method
Qin et al. The implementation of an improved NPML absorbing boundary condition in elastic wave modeling
Ren et al. Elastic full-waveform inversion using the second-generation wavelet and an adaptive-operator-length scheme
Xue et al. Reverse-time migration using multidirectional wavefield decomposition method
Liu Removing the stability limit of the time-space domain explicit finite-difference schemes for acoustic modeling with stability condition-based spatial operators
WO2018067014A1 (en) Seismic modeling
CN116520418A (en) Efficient extraction method for elastic wave angle domain common imaging point gather
CN110658558A (en) Front-of-stack depth reverse time migration imaging method and system for absorption attenuation medium
CN113552633B (en) Elastic wave frequency dispersion pressing method for optimizing difference coefficient and longitudinal and transverse wave separation FCT
Drijkoningen et al. Seismic velocity structure of oceanic crust by inversion using genetic algorithms
Wang* et al. An unsplit convolutional perfectly matched layer for visco-acoustic wave equation with fractional time derivatives
CN111208568B (en) Time domain multi-scale full waveform inversion method and system
CN114942472A (en) Offset imaging method and equipment based on uplink ray tracing strategy
Yang et al. An improved convolution perfectly matched layer for elastic second-order wave equation
Ying et al. Denoise investigation on prestack reverse time migration based on GPU/CPU collaborative parallel accelerating computation
CN112649848A (en) Method and apparatus for solving seismic wave impedance using wave equation
CN113917526B (en) Forward modeling method based on non-split complete matching layer absorption boundary
Bregman et al. A noniterative procedure for inverting plane‐wave reflection data at several angles of incidence using the Riccati equation

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
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210806

CF01 Termination of patent right due to non-payment of annual fee