WO2019233475A1 - 拓展有限差分稳定性条件的波场模拟方法、设备及介质 - Google Patents

拓展有限差分稳定性条件的波场模拟方法、设备及介质 Download PDF

Info

Publication number
WO2019233475A1
WO2019233475A1 PCT/CN2019/090303 CN2019090303W WO2019233475A1 WO 2019233475 A1 WO2019233475 A1 WO 2019233475A1 CN 2019090303 W CN2019090303 W CN 2019090303W WO 2019233475 A1 WO2019233475 A1 WO 2019233475A1
Authority
WO
WIPO (PCT)
Prior art keywords
time
wave
wave field
time step
dispersion
Prior art date
Application number
PCT/CN2019/090303
Other languages
English (en)
French (fr)
Inventor
高英杰
张金海
姚振兴
Original Assignee
中国科学院地质与地球物理研究所
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 中国科学院地质与地球物理研究所 filed Critical 中国科学院地质与地球物理研究所
Priority to CN201980036883.2A priority Critical patent/CN112805598B/zh
Priority to US16/973,230 priority patent/US20210239870A1/en
Publication of WO2019233475A1 publication Critical patent/WO2019233475A1/zh

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/32—Transforming one recording into another or one representation into another
    • 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/282—Application of seismic models, synthetic seismograms
    • 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
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364—Seismic filtering
    • G01V1/368—Inverse filtering

Definitions

  • This application relates to the field of wave field simulation, and in particular, to a wave field simulation method, equipment, and medium for expanding finite difference stability conditions.
  • a sufficiently small time step can obtain high-precision numerical simulation results without time dispersion even in a long-term numerical simulation, but a small time step means more calculation amount and lower calculation efficiency; and
  • the maximum time step of the explicit difference is limited by the Courant-Friedrichs-Lewy (CFL) stability condition.
  • CFL Courant-Friedrichs-Lewy
  • the calculation efficiency is higher, but the calculation accuracy is higher. Low, even in short-term simulation results, there will still be more severe time dispersion. Therefore, in a speed model with a fine structure and a high-speed body, the number of iterations of numerical simulation is increased because a small time step can be used, which reduces the calculation efficiency.
  • Time dispersion has been proven to be independent of space dispersion, and various methods have been proposed to remove time dispersion. This allows us to use larger time steps (close to the upper limit of the stable conditions) for numerical simulation without worrying about time dispersion. Although the problem of time dispersion with large time steps has been solved, as of now, the time step of explicit finite difference in wave field simulation is still limited by the stability conditions and cannot be made too large.
  • Implicit finite difference is unconditionally stable. Regardless of calculation accuracy, any time step can be used without being limited by the CFL stability conditions. Implicit finite difference methods need to solve large matrices, which require more calculations than explicit differences. Such methods have lower calculation efficiency.
  • the numerical simulation method of the spatial filtering extended explicit differential stability condition in the field of electromagnetic wave simulation has not been applied to the field of wave field simulation.
  • the spatial filtering method in electromagnetic waves does not solve the problem of time dispersion when using large time steps.
  • the time dispersion of the wave field simulation using a large time step simulation is very serious, so it is not enough to expand the stability conditions.
  • the problem of time dispersion must also be solved.
  • this application proposes a method and equipment for extending the wave field simulation condition of finite difference stability.
  • medium using the spatial filtering method in the wave number domain to filter out wave fields larger than the threshold Kmax, and to implement extended explicit finite difference stability conditions in the field of wave field simulation, and introduce inverse time dispersion transform to remove large time steps
  • the time dispersion produced by the long simulation so as to ensure that the stability conditions are extended, and then a numerical simulation with a large time step is used.
  • the purpose is to develop a new method and tool for numerical simulation of wave field with high precision and high efficiency.
  • the present application provides a method for simulating a wave field with extended finite difference stability conditions, which is characterized in that the method includes: performing a time step iteration based on a wave field numerical simulation model, and using a spatial filtering method to filter out high wave numbers Unstable wave field components in the region; all time step iterations are over and the wave field data obtained by time iteration is saved; the target wave field data that needs to be output in the wave field data is subjected to inverse time dispersion transformation operation to obtain the time frequency elimination Scattered wave field data; storing the wave field data after the inverse time dispersion transform process.
  • the wave-field-based numerical simulation model performs time step iteration, and uses a spatial filtering method to filter out unstable wave field components distributed in a high wave number region, further including: based on a given time step and time In a discrete format, a positive time dispersion transformation is performed for the first source time function to obtain a second source time function; wherein the second source time function is used as an input to the wave field numerical simulation model instead of the first source time function. Time step iteration.
  • the wave field-based numerical simulation model performs time step iteration and uses a spatial filtering method to filter out unstable wave field components distributed in a high wave number region, including the following steps: based on the wave field numerical simulation model, performing time Step size iteration to calculate the maximum wave number threshold corresponding to a given time step that satisfies the CFL stability condition; after the time step iteration is over, transform the space domain to the wave number domain; use a low-pass filter to filter out larger than the The wave field component of the maximum wave number threshold; the wave number domain after the filtering is transformed back to the spatial domain; and the next time step iteration is started.
  • parameters of the wave field numerical simulation model include a spatial grid size, a maximum speed of the model, a time step, a discrete format of time and space variables.
  • the CFL stability condition corresponding to the maximum wave number threshold is a critical state of the CFL stability condition.
  • the transformation from the space domain to the wave number domain uses a Fourier transform; the transformation from the wave number domain to the space domain uses an inverse Fourier transform.
  • F (k) is a low-pass filter
  • k is a wave number
  • an upper limit of the low-pass filter is a maximum wave number threshold Kmax.
  • the method of the inverse time dispersion transform operation includes the steps of calculating an actual phase shift ⁇ ( ⁇ , ⁇ t) for the effective frequency using a dispersion relationship after time dispersion, where ⁇ is a frequency and ⁇ t is a time step; Transformation: Where u (t) is the destination wavefield data to be output; perform inverse Fourier transform: The wave field data u ′ (t) with time dispersion removed is obtained.
  • the present application also provides a wave field simulation device for expanding the stability of the finite difference, which is characterized in that the device includes a filtering unit, an intermediate wave field data storage unit, an inverse time dispersion transform unit, and a resulting wave field data storage unit.
  • the filtering unit is used to perform a time step iteration based on a wave field numerical simulation model, and a spatial filtering method is used to filter out unstable wave field components distributed in a high wave number region; the intermediate wave field data storage unit is used in all After the time step iteration is over, the wave field data obtained by time iteration is saved; the inverse time dispersion transform unit is configured to perform an inverse time dispersion transform operation on the target wave field data that needs to be output in the wave field simulation data to obtain the elimination time. Dispersive wave field data; the resulting wave field data storage unit is configured to store the wave field data after the inverse time dispersion transformation process.
  • the device further includes a positive time dispersion transform unit, and the positive time dispersion transform unit performs a positive time dispersion transform on the first source time function based on a given time step and a time discrete format to obtain a second A source function; wherein the second source time function replaces the first source time function as an input to the wave field numerical simulation model and performs time step iteration.
  • the filtering unit includes a time step iteration unit, a spatial domain transformation unit, a low-pass filter, and a wave number domain transformation unit, and the time step iteration unit is configured to be based on a wave field numerical simulation model and the second seismic source.
  • the function performs time step iteration and calculates the maximum wave number threshold corresponding to a given time step that satisfies the CFL stability condition;
  • the spatial domain transformation unit is configured to transform the spatial domain to a wave number after the time step iteration ends
  • the low-pass filter is used to filter out a wave field component that is greater than the maximum wave number threshold;
  • the wave number domain transformation unit is used to transform the filtered wave number domain back to the spatial domain to start the next time Step size iteration.
  • the present application also provides an electronic device including a memory, a processor, and a computer program stored on the memory and executable on the processor, characterized in that the processor implements the extended finite difference as described above when the program is executed. Wave field simulation method for stability conditions.
  • the present application also provides a computer-readable storage medium on which a processor program is stored, which is characterized in that the processor program is used to execute the wave field simulation method for expanding the finite difference stability condition described above.
  • the technical solution provided in this application is based on wave field numerical simulation.
  • a spatial filtering method is adopted to filter out unstable wave field components distributed in a high wave number region.
  • the wave field simulation is still stable at the time step; the inverse time dispersion transform is used to remove the time dispersion caused by using a large time step, and the final simulation result is almost the same as the simulation result using a small time step, which guarantees that Efficiency and precision of numerical simulation.
  • FIG. 1 is a schematic flowchart of a method according to an embodiment of the present application.
  • FIG. 2 is a schematic flowchart of a method according to another embodiment of the present application.
  • FIG. 3 is a schematic flowchart of a high wave number domain filtering method according to an embodiment of the present application.
  • FIG. 4 is a schematic diagram of a device composition according to an embodiment of the present application.
  • FIG. 5 is a schematic diagram of a device composition according to another embodiment of the present application.
  • FIG. 6 is a schematic diagram of a composition of a filtering unit according to an embodiment of the present application.
  • FIG. 8 is a single-channel waveform record obtained by using different time step simulations according to an embodiment of the present application.
  • FIG. 9 is a single-channel waveform record obtained by using different time step simulations according to another embodiment of the present application.
  • FIG. 4 is a schematic diagram of a device composition according to an embodiment of the present application, including a filtering unit 1, an intermediate wave field data storage unit 2, an inverse time dispersion transform unit 3, and a resulting wave field data storage unit 4.
  • the filtering unit 1 is used to perform time step iteration based on a wave field numerical simulation model, and the spatial filtering method is used to filter out unstable wave field components distributed in a high wave number region; the intermediate wave field data storage unit 2 is used to perform all time steps After the long iteration is over, the wave field data obtained by time iteration is saved; the inverse time dispersion transform unit 3 is used to perform an inverse time dispersion transform operation on the target wave field data that needs to be output in the wave field data to obtain a wave that eliminates time dispersion. Field data; the resulting wavefield data storage unit 4 is used to store wavefield data after inverse time dispersion transform processing.
  • FIG. 5 is a schematic diagram of a device composition according to another embodiment of the present application, including a positive time dispersion transform unit 5, a filtering unit 1, an intermediate wave field data storage unit 2, an inverse time dispersion transform unit 3, and a resulting wave field data storage unit. 4.
  • the positive time dispersion transform unit 5 performs a positive time dispersion transform on the first source time function based on a given time step and time discrete format to obtain a second source function, where the second source time function replaces the first source time function Perform time step iterations as input to a wave field numerical simulation model.
  • the filtering unit 1 is used to perform time step iteration based on a wave field numerical simulation model, and the spatial filtering method is used to filter out unstable wave field components distributed in a high wave number region; the intermediate wave field data storage unit 2 is used to perform all time steps After the long iteration is over, the wave field data obtained by time iteration is saved; the inverse time dispersion transform unit 3 is used to perform an inverse time dispersion transform operation on the target wave field data that needs to be output in the wave field data to obtain a wave that eliminates time dispersion. Field data; the resulting wavefield data storage unit 4 is used to store wavefield data after inverse time dispersion transform processing.
  • FIG. 6 is a schematic diagram of a composition of a filtering unit according to an embodiment of the present application, including a time step iteration unit 11, a spatial domain transformation unit 12, a low-pass filter 13, and a wavenumber domain transformation unit 14.
  • a time step iteration unit 11 is used to perform a time step iteration based on a wave field numerical simulation model to calculate a maximum wave number threshold corresponding to a given time step that satisfies CFL stability conditions; a spatial domain transformation unit 12 is used to After the step size iteration is completed, the spatial domain is transformed into the wave number domain; the low-pass filter 13 is used to filter out wave field components larger than the maximum wave number threshold; the wave number domain transformation unit 14 is used to transform the filtered wave number domain back to space Domain, starting the next time step iteration.
  • FIG. 1 is a schematic flowchart of a method according to an embodiment of the present application, including the following steps.
  • step S1 based on a wave field numerical simulation model, time step iteration is performed, and a spatial filtering method is used to filter out unstable wave field components distributed in a high wave number region.
  • FIG. 3 is a schematic flowchart of a high wave number domain filtering method of the present application, including the following steps.
  • step S11 based on the wave field numerical simulation model, time step iteration is performed to calculate a maximum wave number threshold corresponding to a given time step and satisfying the CFL stability condition.
  • the parameters of the wave field numerical simulation model include the spatial grid size, the maximum speed of the model, the time step, and the discrete format of time and space variables.
  • the CFL stability condition corresponding to the maximum wave number threshold Kmax is a critical state of the CFL stability condition.
  • the CFL condition was named after Courant, Friedrichs, and Lewy. They first proposed this concept in an article about the finite difference method of partial differential equations in 1928, and it was not used to analyze the stability of the difference scheme. In fact, the finite difference method is used as an analysis tool to prove the existence of the solution of some partial differential equations.
  • the basic idea is to first construct a difference equation to obtain a sequence of approximation solutions. As long as you know that the approximation sequence converges under a given grid system, it is easy to prove that the convergence solution is the solution of the original differential equation.
  • the finite difference method and the finite volume method are more and more used in the numerical simulation of fluid mechanics.
  • the CFL condition also becomes very important. But it is worth noting that the CFL condition is only a necessary condition for stability (convergence), not a sufficient condition.
  • step S12 after the time step iteration ends, the space domain is transformed into the wave number domain.
  • step S13 a low-pass filter is used to filter out wave field components larger than the maximum wave number threshold.
  • the upper limit of the low-pass filter is Kmax.
  • the low-pass filter expression is:
  • F (k) is a low-pass filter
  • k is a wave number
  • an upper limit of the low-pass filter is Kmax.
  • the unstable components mainly come from the high wave number domain that exceeds the corresponding wave number threshold.
  • the wave field distributed in space is transformed into the wave number domain by fast Fourier transform, and the wave field above the corresponding threshold is filtered out in the wave number domain, and then inverse Fourier is used.
  • the transformation transforms the wave field into the spatial domain.
  • the unstable wave field components caused by large time steps will be multiplied and not accumulated, and the phenomenon of numerical simulation instability will not occur.
  • the effective wave field components are mainly distributed in the low wave number region and have not been filtered out, so the idea of spatial filtering will not affect the effective wave field.
  • step S14 the inverse Fourier transform is used to transform the filtered wave number domain back to the spatial domain.
  • step S2 all time step iterations end, and the wave field data obtained by the time iteration is saved.
  • Step S3 Perform inverse time dispersion transform operation on the target wave field data that needs to be output from the wave field data to obtain wave field data that eliminates time dispersion.
  • Inverse time dispersion is a post-processing technique that uses inverse time dispersion transform to efficiently remove time dispersion and improve simulation accuracy.
  • the method of inverse time dispersion transform operation includes the following steps:
  • step S4 the wave field data after the inverse time dispersion transformation process is stored.
  • the technical solution provided in this application is based on wave field numerical simulation.
  • a spatial filtering method is adopted to filter out unstable wave field components distributed in a high wave number region.
  • the wave field simulation is still stable at the time step; the inverse time dispersion transform is used to remove the time dispersion caused by using a large time step, and the final simulation result is almost the same as the simulation result using a small time step, which guarantees that Efficiency and precision of numerical simulation.
  • FIG. 2 is a schematic flowchart of a method according to another embodiment of the present application, including the following steps.
  • Step S0 Based on the given time step and time discrete format, perform a positive time dispersion transformation on the first source time function to obtain a second source function, where the second source time function replaces the first source time function as the wave field value.
  • the input of the simulation model is iterated over time steps.
  • the positive time dispersion transform can be regarded as a modified Fourier transform.
  • S (t) is the first source function
  • step S1 based on a wave field numerical simulation model, time step iteration is performed, and a spatial filtering method is used to filter out unstable wave field components distributed in a high wave number region.
  • step S2 until the end of all time step iterations, the wave field data obtained by the time iteration is saved.
  • Step S3 Perform inverse time dispersion transform operation on the target wave field data that needs to be output from the wave field data to obtain wave field data that eliminates time dispersion.
  • step S4 the wave field data after the inverse time dispersion transformation process is stored.
  • the inverse time dispersion transform since the positive time dispersion transform was previously used for the source time function, the inverse time dispersion transform more accurately repairs the time dispersion in the wave field.
  • steps S1, S2, S3, and S4 are the same as those in the foregoing embodiment, and details are not described herein again.
  • An electronic device includes a memory, a processor, and a computer program stored on the memory and can run on the processor.
  • the processor implements the above-mentioned wave field simulation method for extending the finite difference stability condition when the program is executed.
  • a computer-readable storage medium stores a processor program thereon, and the processor program is used to execute the above-mentioned wave field simulation method for expanding finite difference stability conditions.
  • the computer-readable storage medium may be a read-only memory (ROM), a random access memory (RAM), a compact disc read-only memory (CD-ROM), a magnetic tape, a floppy disk, and an optical data storage device. This is the limit.
  • the device includes a filtering unit 1, an intermediate wave field data storage unit 2, an inverse time dispersion transform unit 3, and a resulting wave field data storage unit 4.
  • the filtering unit 1 is used to perform time step iteration based on a wave field numerical simulation model, and the spatial filtering method is used to filter out unstable wave field components distributed in a high wave number region;
  • the intermediate wave field data storage unit 2 is used to perform all time steps After the long iteration is over, the seismic trace obtained by time iteration is saved;
  • the inverse time dispersion transform unit 3 is configured to perform an inverse time dispersion transform operation on the target seismic trace that needs to be output in the wave field data to obtain a seismic trace with time dispersion eliminated;
  • the resulting wavefield data storage unit 4 is used to store seismic traces after inverse time dispersion transform processing.
  • Seismic traces are waveform records of spatial point wave field values distributed along time. Each time point corresponds to a seismic waveform record, which is called a single seismic trace, and one seismic trace constitutes a time seismic wave field.
  • the filtering unit 1 includes a time step iteration unit 11, a spatial domain transformation unit 12, a low-pass filter 13, and a wavenumber domain transformation unit 14.
  • the time step iteration unit 11 is used to perform the time step iteration based on the wave field numerical simulation model, and calculates the maximum wave number threshold corresponding to the given time step that satisfies the CFL stability condition; the spatial domain transformation unit 12 is used to time step After the iteration, the spatial domain is transformed into the wave number domain; the low-pass filter 13 is used to filter out wave field components larger than the maximum wave number threshold; the wave number domain transformation unit 14 is used to transform the filtered wave number domain back to the spatial domain and start Time step iterations.
  • the method includes the following steps.
  • step S1 based on a wave field numerical simulation model, time step iteration is performed, and a spatial filtering method is used to filter out unstable wave field components distributed in a high wave number region.
  • the first seismic source time function is used as an input of a wave field numerical simulation model to perform a time step iteration.
  • step S11 based on the wave field numerical simulation model, time step iteration is performed to calculate a maximum wave number threshold Kmax corresponding to the CFL stability condition corresponding to a given time step.
  • Two-dimensional acoustic wave equations are used for simulation in a homogeneous medium.
  • Pseudospectral method is used for spatial dispersion.
  • Second-order and fourth-order difference formats are used for time dispersion.
  • the space grid size is 10m.
  • a1 is a snapshot of the pseudo-spectral time second-order difference without spatial filtering. It can be seen that the wave field is unstable.
  • Fig. 7 b1 is the wave field of a1 in Fig. 7 after the two-dimensional Fourier transform. The wave field in the wave number domain uses a low-pass filter. As shown by line b1 in FIG. 7, after filtering the wave field in the high wave number region, b 2 in FIG. 7 can be obtained, and the wave field in b 2 in FIG. 7 is two. The inverse Fourier transform can be used to obtain the wave field in the space domain in a2 in Fig. 7. The original instability has been filtered out. For the same reason, spatial filtering is still effective for the unstable wave field of the 4th order difference format of the pseudospectral method in a3 in FIG. 7. The filtered wave field snapshot is shown in b4 in FIG. 7.
  • step S12 after the time step iteration is over, a Fourier transform is used to transform the space domain to the wave number domain.
  • step S13 a low-pass filter is used to filter out wave field components larger than the maximum wave number threshold Kmax.
  • the upper limit of the low-pass filter is Kmax.
  • the low-pass filter expression is:
  • F (k) is a low-pass filter
  • k is a wave number
  • an upper limit of the low-pass filter is Kmax.
  • step S14 the inverse Fourier transform is used to transform the filtered wave number domain back to the spatial domain.
  • step S2 all time step iterations end, and the wave field data obtained by the time iteration is saved.
  • step S3 the inverse time dispersion transform operation is performed on the target seismic track that needs to be output in the wave field data to obtain a seismic track with time dispersion eliminated.
  • Inverse time dispersion is a post-processing technique. Using inverse time dispersion transform can effectively remove time dispersion and improve simulation accuracy.
  • the method of inverse time dispersion transform operation includes the following steps:
  • FIG. 8 is a single-channel waveform record obtained by using different time step simulations according to an embodiment of the present application.
  • (A) in FIG. 8 is a single-channel waveform record obtained by using different time step simulations and performing spatial filtering.
  • the upper three figures ⁇ t1 in (a) and (b) are pseudo-spectral time 2-order difference formats
  • step S4 the seismic traces after the inverse time dispersion transform processing are stored.
  • the above embodiment illustrates that the use of spatial filtering can expand the CFL stability conditions of explicit differentials, and further the use of inverse time dispersion transformation can ensure the accuracy of numerical simulation with large time steps.
  • the device includes a positive time dispersion transform unit 5, a filtering unit 1, an intermediate wave field data storage unit 2, a reverse time dispersion transform unit 3, and a resulting wave field data storage unit 4.
  • the positive time dispersion transform unit 5 performs a positive time dispersion transform on the first source time function based on a given time step and time discrete format to obtain a second source function, where the second source time function replaces the first source time function Perform time step iterations as input to a wave field numerical simulation model.
  • the filtering unit 1 is used to perform time step iteration based on a wave field numerical simulation model, and the spatial filtering method is used to filter out unstable wave field components distributed in a high wave number region;
  • the intermediate wave field data storage unit 2 is used to perform all time steps After the long iteration is over, the seismic trace obtained by time iteration is saved;
  • the inverse time dispersion transform unit 3 is configured to perform an inverse time dispersion transform operation on the target seismic trace that needs to be output in the wave field data to obtain a seismic trace with time dispersion eliminated;
  • the resulting wavefield data storage unit 4 is used to store seismic traces after inverse time dispersion transform processing.
  • Seismic traces are waveform records of spatial point wave field values distributed along time. Each time point corresponds to a seismic waveform record, which is called a single seismic trace, and one seismic trace constitutes a time seismic wave field.
  • the filtering unit 1 includes a time step iteration unit 11, a spatial domain transformation unit 12, a low-pass filter 13, and a wavenumber domain transformation unit 14.
  • the time step iteration unit 11 is used to perform the time step iteration based on the wave field numerical simulation model, and calculates the maximum wave number threshold corresponding to the given time step that satisfies the CFL stability condition; the spatial domain transformation unit 12 is used to time step After the iteration, the spatial domain is transformed into the wave number domain; the low-pass filter 13 is used to filter out wave field components larger than the maximum wave number threshold; the wave number domain transformation unit 14 is used to transform the filtered wave number domain back to the spatial domain and start Time step iterations.
  • the method includes the following steps.
  • Step S0 Based on the given time step and time discrete format, perform a positive time dispersion transformation on the first source time function to obtain a second source function, where the second source time function replaces the first source time function as the wave field value
  • the input of the simulation model is iterated over time steps.
  • the positive time dispersion transform can be regarded as a modified Fourier transform.
  • S (t) is the first source function
  • step S1 based on a wave field numerical simulation model, time step iteration is performed, and a spatial filtering method is used to filter out unstable wave field components distributed in a high wave number region.
  • the second source time function is used as the input of the wave field numerical simulation model to perform time step iteration.
  • Other specific implementations of step S1 are the same as those in the first embodiment, and will not be described again.
  • step S2 all time step iterations end, and the wave field data obtained by the time iteration is saved.
  • step S3 the inverse time dispersion transform operation is performed on the target seismic track that needs to be output in the wave field data to obtain a seismic track with time dispersion eliminated.
  • FIG. 9 is a single-channel waveform record obtained by using different time step simulations according to another embodiment of the present application.
  • (A) in FIG. 9 is a single-channel waveform record obtained by using different time step simulations and performing spatial filtering.
  • the upper five graphs ⁇ t1 in (a) and (b) are pseudo-spectral time second-order difference formats
  • step S4 the seismic traces after the inverse time dispersion transformation process are stored.
  • steps S2, S3, and S4 is the same as that in Embodiment 1, and details are not described again.
  • the inverse time dispersion transform since the positive time dispersion transform was previously used for the source time function, the inverse time dispersion transform more accurately repairs the time dispersion in the wave field.

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

一种拓展有限差分稳定性条件的波场模拟方法，包括：基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除出现在高波数区域的波场分量（S1）；所有时间步长迭代结束，保存时间迭代所得的波场数据（S2）；针对波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据（S3）；存储逆时间频散变换处理之后的波场数据（S4）。还提供用于执行拓展有限差分稳定性条件的波场模拟方法的电子设备和计算机可读存储介质。

Description

拓展有限差分稳定性条件的波场模拟方法、设备及介质 技术领域
本申请涉及波场模拟领域，具体涉及拓展有限差分稳定性条件的波场模拟方法、设备及介质。
背景技术
显式有限差分由于其形式简单易于实现，是目前波场模拟领域应用最广泛的时间偏导数求解方法。其中模拟精度和计算效率是显式有限差分所关心的最重要的两个问题。而时间步长的选取会同时影响数值模拟的精度和效率。
足够小的时间步长即使在长时程数值模拟中仍然可以得到高精度的数值模拟结果而不出现时间频散，但是小的时间步长意味着更多的计算量，计算效率较低；与之相反，显式差分的最大时间步长受到Courant-Friedrichs-Lewy(CFL)稳定性条件的限制，采用接近于CFL稳定性条件上限的大时间步长时，计算效率较高，但是计算精度较低，即使在短时程的模拟结果中仍然会出现较为严重时间频散。于是，在具有精细结构、高速体的速度模型中，由于只能采用小时间步长而增加了数值模拟的迭代次数，降低了计算效率。
时间频散被证明独立于空间频散，人们提出了多种去除时间频散的方法。这使得我们可以采用更大的时间步长(接近于稳定条件上限)进行数值模拟而不用担心时间频散。虽然大时间步长的时间频散问题得到了解决，截至目前，波场模拟领域显式有限差分的时间步长仍然受到稳定条件的限制而不能取得过大。
隐式有限差分是无条件稳定的，不考虑计算精度的情况下，可以采用任意的时间步长而不受CFL稳定性条件的限制。隐式有限差分法需要求解大型的矩阵，比显示差分需要更多的计算量，这类方法的计算效率较低。
电磁波模拟领域出现了采用空间滤波的思路实现显式有限差分无条件稳定的方法，这是首次针对显式差分进行稳定性条件拓展的数值模拟方法。
电磁波模拟领域的空间滤波拓展显式差分稳定性条件的数值模拟方法还没有应用于波场模拟领域。同时，电磁波中的空间滤波方法并没有解决采用大时间步长所面临的时间频散问题。而波场模拟采用大时间步长模拟的时间频散非常严重，因此仅实现稳定性条件拓展是不够的，还必须解决时间频散的问题。
发明内容
为了解决波场模拟领域的时间步长受CFL稳定性条件限制的问题，并且解决大时间步长模拟时的时间频散的问题，本申请提出拓展有限差分稳定性条件的波场模拟方法、设备及介质，在波数域内采用空间滤波的方法，将大于阈值Kmax的波场滤除，在波场模拟领域实现拓展显式有限差分稳定性条件，并且将逆时间频散变换引入以去除大时间步长模拟所产生的时间频散，从而保证稳定性条件拓展之后采用大时间步长的数值模拟。旨在发展一种高精度、高效率的波场数值模拟新方法、新工具。
本申请提供一种拓展有限差分稳定性条件的波场模拟方法，其特征在于，所述方法包括：基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量；所有时间步长迭代结束，保存时间迭代所得的波场数据；针对所述波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据；存储所述的逆时间频散变换处理之后的波场数据。
进一步地，所述基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量之前，还包括：基于给定的时间步长和时间离散格式，针对第一震源时间函数进行正时间频散变换得到第二震源时间函数；其中，所述第二震源时间函数代替所述第一震源时间函数作为所述波场数值模拟模型的输入进行时间步长迭代。
进一步地，所述基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量，包括如下步骤：基于波场数值模拟模型，进行时间步长迭代，计算出给定时间步长对应的满足CFL稳定性条件的最大波数阈值；所述时间步长迭代结束后，将空间域变换到波数域；采用低通滤波器滤除大于所述最大波数阈值的波场分量；所述滤波之后的波数域变换回所述空间域；开始下次的时间步长迭代。
进一步地，所述波场数值模拟模型的参数包括空间网格大小、模型最大速度、时间步长、时间和空间变量的离散格式。
进一步地，所述最大波数阈值对应的CFL稳定性条件为CFL稳定性条件的临界态。
进一步地，所述空间域变换到所述波数域采用傅立叶变换；所述波数域变换到所述空间域采用傅立叶逆变换。
进一步地，所述低通滤波器表达式为：
Figure PCTCN2019090303-appb-000001
其中，F(k)为低通滤波器，k是波数，所述低通滤波器的上限是最大波数阈值Kmax。
进一步地，所述逆时间频散变换操作的方法包括如下步骤：利用时间离散之后频散关系对有效频率计算实际相移θ(ω,Δt)，其中ω是频率，Δt是时间步长；应用变换：
Figure PCTCN2019090303-appb-000002
其中u(t)是需要输出的目的波场数据；进行逆傅里叶变换：
Figure PCTCN2019090303-appb-000003
得到消除时间频散的波场数据u′(t)。
本申请还提供一种拓展有限差分稳定性条件的波场模拟装置，其特征在于，所述装置包括滤波单元、中间波场数据存储单元、逆时间频散变换单元、结果波场数据存储单元，所述滤波单元用于基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量；所述中间波场数据存储单元用于在所有时间步长迭代结束后，保存时间迭代所得的波场数据；所述逆时间频散变换单元用于针对波场模拟数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据；所述结果波场数据存储单元用于存储所述的逆时间频散变换处理之后的波场数据。
进一步地，所述装置还包括正时间频散变换单元，所述正时间频散变 换单元基于给定的时间步长和时间离散格式，针对第一震源时间函数进行正时间频散变换得到第二震源函数；其中，所述第二震源时间函数代替所述第一震源时间函数作为所述波场数值模拟模型的输入进行时间步长迭代。
进一步地，所述滤波单元包括时间步长迭代单元、空间域变换单元、低通滤波器、波数域变换单元，所述时间步长迭代单元用于基于波场数值模拟模型和所述第二震源函数进行时间步长迭代，计算出给定时间步长对应的满足CFL稳定性条件的最大波数阈值；所述空间域变换单元用于在所述时间步长迭代结束后，将空间域变换到波数域；所述低通滤波器用于滤除大于所述最大波数阈值的波场分量；所述波数域变换单元用于将所述滤波之后的波数域变换回所述空间域，开始下次的时间步长迭代。
本申请还提供一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述程序时实现如上所述拓展有限差分稳定性条件的波场模拟方法。
本申请还提供一种计算机可读存储介质，其上存储有处理器程序,其特征在于，该处理器程序用于执行上述所述拓展有限差分稳定性条件的波场模拟方法。
本申请提供的技术方案，基于波场数值模拟，在采用大时间步长时，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量，使得即使采用超过CFL稳定性条件的大时间步长时波场模拟仍然稳定；应用逆时间频散变换来去除由于采用大时间步长所产生的时间频散，最终的模拟结果与采用小时间步长的模拟结果精度几乎一致，保证了数值模拟的效率及精度。
附图说明
图1是本申请一实施例提供的方法流程示意图；
图2是本申请另一实施例提供的方法流程示意图；
图3是本申请一实施例提供的高波数域滤波方法流程示意图；
图4是本申请一实施例提供的装置组成示意图；
图5是本申请另一实施例提供的装置组成示意图；
图6是本申请一实施例提供的滤波单元组成示意图；
图7是本申请一实施例提供的采用时间步长Δt＝3ms进行模拟，333ms时的波场快照；
图8是本申请一实施例提供的采用不同时间步长模拟得到的单道波形记录；
图9是本申请另一实施例提供的采用不同时间步长模拟得到的单道波形记录。
具体实施方式
为使本申请实施例的目的、技术方案和优点更加清楚，以下将结合附图和实施例，对本申请技术方案的具体实施方式进行更加详细、清楚的说明。然而，以下描述的具体实施方式和实施例仅是说明的目的，而不是对本申请的限制。其只是包含了本申请一部分实施例，而不是全部的实施例，本领域技术人员对于本申请的各种变化获得的其他实施例，都属于本申请保护的范围。
图4是本申请一实施例提供的装置组成示意图，包括滤波单元1、中间波场数据存储单元2、逆时间频散变换单元3、结果波场数据存储单元4。
滤波单元1用于基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量；中间波场数据存储单元2用于在所有时间步长迭代结束后，保存时间迭代所得的波场数据；逆时间频散变换单元3用于针对波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据；结果波场数据存储单元4用于存储逆时间频散变换处理之后的波场数据。
图5是本申请另一实施例提供的装置组成示意图，包括正时间频散变换单元5、滤波单元1、中间波场数据存储单元2、逆时间频散变换单元3、结果波场数据存储单元4。
正时间频散变换单元5基于给定的时间步长和时间离散格式，针对第一震源时间函数进行正时间频散变换得到第二震源函数，其中，第二震源 时间函数代替第一震源时间函数作为波场数值模拟模型的输入进行时间步长迭代。滤波单元1用于基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量；中间波场数据存储单元2用于在所有时间步长迭代结束后，保存时间迭代所得的波场数据；逆时间频散变换单元3用于针对波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据；结果波场数据存储单元4用于存储逆时间频散变换处理之后的波场数据。
图6是本申请一实施例提供的滤波单元组成示意图，包括时间步长迭代单元11、空间域变换单元12、低通滤波器13、波数域变换单元14。
时间步长迭代单元11，用于基于波场数值模拟模型进行时间步长迭代，计算出给定时间步长对应的满足CFL稳定性条件的最大波数阈值；空间域变换单元12，用于在时间步长迭代结束后，将空间域变换到波数域；低通滤波器13，用于滤除大于最大波数阈值的波场分量；波数域变换单元14，用于将滤波之后的波数域变换回空间域，开始下次的时间步长迭代。
图1是本申请一实施例提供的方法流程示意图，包括如下步骤。
步骤S1，基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量。
如图3所示，图3是本申请高波数域滤波方法流程示意图，包括如下步骤。
步骤S11，基于波场数值模拟模型，进行时间步长迭代，计算出给定时间步长对应的满足CFL稳定性条件的最大波数阈值。
波场数值模拟模型的参数包括空间网格大小、模型最大速度、时间步长、时间和空间变量的离散格式。
最大波数阈值Kmax对应的CFL稳定性条件为CFL稳定性条件的临界态。
CFL条件是以Courant，Friedrichs，Lewy三个人的名字命名的，他们最早在1928年一篇关于偏微分方程的有限差分方法的文章中首次提出这个概念的时候，并不是用来分析差分格式的稳定性，而是仅仅以有限差分方法作为分析工具来证明某些偏微分方程的解的存在性的。其基本思想 是先构造差分方程得到一个逼近解的序列，只要知道在给定的网格系统下这个逼近序列收敛，那么就很容易证明这个收敛解就是原微分方程的解。
要使这个逼近序列收敛，必须满足一个条件，那就是著名的CFL条件，记述如下：一个数值方法只有在其依赖的数值域内包含双曲型方程的差分格式收敛的必要条件是差分格式的依赖域包含了微分方程的依赖域。
随着计算机的迅猛发展，有限差分方法和有限体积方法越来越多的应用于流体力学的数值模拟中，CFL条件作为一个格式稳定性和收敛性的判据，也随之显得非常重要了。但值得注意的是，CFL条件仅仅是稳定性(收敛性)的必要条件，而不是充分条件。
步骤S12，时间步长迭代结束后，空间域变换到波数域。
步骤S13，采用低通滤波器滤除大于最大波数阈值的波场分量。
低通滤波器的上限是Kmax。
低通滤波器表达式为：
Figure PCTCN2019090303-appb-000004
其中，F(k)为低通滤波器，k是波数，低通滤波器的上限是Kmax。
实验证明，采用超过CFL稳定性条件限制的时间步长进行数值模拟时，不稳定分量主要来自超过相应波数阈值的高波数域。只要在每个时间步长迭代结束后，通过快速傅里叶变换将空间内分布的波场变换到波数域，在波数域内将高于相应阈值的波场滤除掉，之后采用逆傅里叶变换将波场变换至空间域，由于大时间步长所造成的不稳定波场分量就倍率出而不会累积，数值模拟不稳定的现象就不会产生。与此同时，有效的波场分量主要分布在低波数区域并没有被滤除，因此空间滤波的思路并不会影响有效波场。
步骤S14，采用傅立叶逆变换将滤波之后的波数域变换回空间域。
步骤S2，所有时间步长迭代结束，保存时间迭代所得的波场数据。
步骤S3，针对波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据。
采用大时间步长会产生时间频散，降低数值模拟的精度。逆时间频散是一种后处理的技术采用逆时间频散变换可以有效率去除时间频散，提高模拟精度。
逆时间频散变换操作的方法包括如下步骤：
利用时间离散之后频散关系对有效频率计算实际相移θ(ω,Δt)，其中ω是频率，Δt是时间步长；
应用变换：
Figure PCTCN2019090303-appb-000005
其中u(t)是需要输出的目的波场数据；
进行逆傅里叶变换：
Figure PCTCN2019090303-appb-000006
得到消除时间频散的波场数据u′(t)。
步骤S4，存储逆时间频散变换处理之后的波场数据。
本申请提供的技术方案，基于波场数值模拟，在采用大时间步长时，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量，使得即使采用超过CFL稳定性条件的大时间步长时波场模拟仍然稳定；应用逆时间频散变换来去除由于采用大时间步长所产生的时间频散，最终的模拟结果与采用小时间步长的模拟结果精度几乎一致，保证了数值模拟的效率及精度。
图2是本申请另一实施例提供的方法流程示意图，包括如下步骤。
步骤S0，基于给定的时间步长和时间离散格式，针对第一震源时间函数进行正时间频散变换得到第二震源函数，其中，第二震源时间函数代替第一震源时间函数作为波场数值模拟模型的输入进行时间步长迭代。
正时间频散变换可以看作修正的傅里叶变换。
具体的变换公式为：
Figure PCTCN2019090303-appb-000007
其中S(t)是第一震源函数，
Figure PCTCN2019090303-appb-000008
是第一震源函数在频率域中的量。
Figure PCTCN2019090303-appb-000009
是利用时间离散之后频散关系对有效频率计算的理论相移，其中ω是频率，Δt是时间步长。
继续进行逆傅里叶变换，得到新的第二震源函数，具体的变换公式为：
Figure PCTCN2019090303-appb-000010
步骤S1，基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量。
步骤S2，至所有时间步长迭代结束，保存时间迭代所得的波场数据。
步骤S3，针对波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据。
步骤S4，存储逆时间频散变换处理之后的波场数据。
在本实施例中，由于在先对震源时间函数采用了正时间频散变换，使得逆时间频散变换更准确的修复了波场中的时间频散。
在本实施例中，步骤S1、S2、S3、S4与上述实施例相同，不再赘述。
一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，处理器执行程序时实现上述的拓展有限差分稳定性条件的波场模拟方法。
一种计算机可读存储介质，其上存储有处理器程序,该处理器程序用于执行上述的拓展有限差分稳定性条件的波场模拟方法。
其中，计算机可读存储介质可以是只读存储器(ROM)、随机存取存储器(RAM)、光盘只读存储器(CD-ROM)、磁带、软盘和光数据存储设备等，根据实际情况选择，并不以此为限。
实施例1
在此提供本申请拓展显式差分稳定性条件的波场模拟方法、装置、设备和存储介质在地震波场的具体应用实施例。
如图4所示，装置包括滤波单元1、中间波场数据存储单元2、逆时间频散变换单元3、结果波场数据存储单元4。
滤波单元1用于基于波场数值模拟模型，进行时间步长迭代，采用空 间滤波的方法滤除分布在高波数区域的不稳定波场分量；中间波场数据存储单元2用于在所有时间步长迭代结束后，保存时间迭代所得的地震道；逆时间频散变换单元3用于针对波场数据中需要输出的目的地震道进行逆时间频散变换操作，得到消除时间频散的地震道；结果波场数据存储单元4用于存储逆时间频散变换处理之后的地震道。
地震道，是空间点波场值沿着时间分布的波形记录，每一个时间点对应一道地震波形记录，称为单个地震道，一道一道地震道组成了时间地震波场。
如图6所示，滤波单元1包括时间步长迭代单元11、空间域变换单元12、低通滤波器13、波数域变换单元14。
时间步长迭代单元11用于基于波场数值模拟模型进行时间步长迭代，计算出给定时间步长对应的满足CFL稳定性条件的最大波数阈值；空间域变换单元12用于在时间步长迭代结束后，将空间域变换到波数域；低通滤波器13用于滤除大于最大波数阈值的波场分量；波数域变换单元14用于将滤波之后的波数域变换回空间域，开始下次的时间步长迭代。
如图1、图3所示，所述方法包括如下步骤。
步骤S1，基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量。
在本实施例中，第一震源时间函数作为波场数值模拟模型的输入，进行时间步长迭代。
步骤S11，基于波场数值模拟模型，进行时间步长迭代，计算出给定时间步长对应的满足CFL稳定性条件的最大波数阈值Kmax。
均匀介质中采用二维声波方程进行模拟，空间离散采用伪谱法，时间离散分别采用2阶差分和4阶差分格式。空间网格大小采用10m。伪谱法(时间2阶差分)和伪谱法(时间4阶差分)由CFL稳定性条件所限定的最大时间步长分别为Δt＝1.125ms和Δt＝1.949ms。这里采用时间步长Δt＝3ms进行模拟，两种格式均会出现波场不稳定的现象。
图7是本申请一实施例提供的采用时间步长Δt＝3ms进行模拟，333ms时的波场快照。
图7中a1是伪谱法时间2阶差分未采用空间滤波的快照，可以看到波场出现不稳定现象，图7中b1是图7中a1的波场经过二维傅里叶变换之后在波数域的波场，采用低通滤波器，图7中b1线所示，将高波数区域的波场波场滤除以后，可以得到图7中b2，对图7中b2中波场做二维逆傅里叶变换，可以得到图7中a2中空间域的波场，原来的不稳定量已经被滤除。同理，针对图7中a3中伪谱法时间4阶差分格式的不稳定波场采用空间滤波仍然有效，滤波后的波场快照如图7中b4所示。
步骤S12，时间步长迭代结束后，采用傅立叶变换将空间域变换到波数域。
步骤S13，采用低通滤波器滤除大于最大波数阈值Kmax的波场分量。低通滤波器的上限是Kmax。
低通滤波器表达式为：
Figure PCTCN2019090303-appb-000011
其中，F(k)为低通滤波器，k是波数，低通滤波器的上限是Kmax。
步骤S14，采用傅立叶逆变换将滤波之后的波数域变换回空间域。
步骤S2，所有时间步长迭代结束，保存时间迭代所得的波场数据。
步骤S3，针对波场数据中需要输出的目的地震道进行逆时间频散变换操作，得到消除时间频散的地震道。
采用大时间步长会产生时间频散，降低数值模拟的精度。逆时间频散是一种后处理的技术，采用逆时间频散变换可以有效率去除时间频散，提高模拟精度。
逆时间频散变换操作的方法包括如下步骤：
利用时间离散之后频散关系对有效频率计算实际相移θ(ω,Δt)，其中ω是频率，Δt是时间步长；
应用变换：
Figure PCTCN2019090303-appb-000012
其中u(t)是需要输出的目的地震道；
进行逆傅里叶变换：
Figure PCTCN2019090303-appb-000013
得到消除时间频散的地震道u′(t)。
图8是本申请一实施例提供的采用不同时间步长模拟得到的单道波形记录。
图8中的(a)为采用不同时间步长模拟，并进行空间滤波所得到的单道波形记录。其中(a)和(b)中上面三个图Δt1是伪谱法时间2阶差分格式，下面三个图Δt2是伪谱法时间4阶差分格式。可以看到即使Δt1＝Δt2＝7ms时，仍未出现不稳定现象。但是时间频散会随着时间步长的不断增大而变大，计算精度随着时间步长的增大而降低。图8中的(b)是对图8中的(a)中的波形采用逆时间频散变换之后得到的波形记录，可以看出时间频散得到了有效的压制。从图8中可知，尤其是伪谱法时间4阶差分格式，在Δt2＝7ms时的计算结果精度仍然非常高。
步骤S4，存储所述的逆时间频散变换处理之后的地震道。
以上实施例说明采用空间滤波可以拓展显式差分的CFL稳定性条件，进一步采用逆时间频散变换可以确保采用大时间步长的数值模拟精度。
实施例2
在此提供本申请拓展显式差分稳定性条件的波场模拟方法、装置、设备和存储介质在地震波场的具体应用实施例。
如图5所示，装置包括正时间频散变换单元5、滤波单元1、中间波场数据存储单元2、逆时间频散变换单元3、结果波场数据存储单元4。
正时间频散变换单元5基于给定的时间步长和时间离散格式，针对第一震源时间函数进行正时间频散变换得到第二震源函数，其中，第二震源时间函数代替第一震源时间函数作为波场数值模拟模型的输入进行时间步长迭代。滤波单元1用于基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量；中间波场数据存储单元2用于在所有时间步长迭代结束后，保存时间迭代所得的地震道；逆时间频散变换单元3用于针对波场数据中需要输出的目的地震道 进行逆时间频散变换操作，得到消除时间频散的地震道；结果波场数据存储单元4用于存储逆时间频散变换处理之后的地震道。
地震道，是空间点波场值沿着时间分布的波形记录，每一个时间点对应一道地震波形记录，称为单个地震道，一道一道地震道组成了时间地震波场。
如图6所示，滤波单元1包括时间步长迭代单元11、空间域变换单元12、低通滤波器13、波数域变换单元14。
时间步长迭代单元11用于基于波场数值模拟模型进行时间步长迭代，计算出给定时间步长对应的满足CFL稳定性条件的最大波数阈值；空间域变换单元12用于在时间步长迭代结束后，将空间域变换到波数域；低通滤波器13用于滤除大于最大波数阈值的波场分量；波数域变换单元14用于将滤波之后的波数域变换回空间域，开始下次的时间步长迭代。
如图2、图3所示，所述方法包括如下步骤。
步骤S0，基于给定的时间步长和时间离散格式，针对第一震源时间函数进行正时间频散变换得到第二震源函数，其中，第二震源时间函数代替第一震源时间函数作为波场数值模拟模型的输入进行时间步长迭代。
正时间频散变换可以看作修正的傅里叶变换。
具体的变换公式为：
Figure PCTCN2019090303-appb-000014
其中S(t)是第一震源函数，
Figure PCTCN2019090303-appb-000015
是第一震源函数在频率域中的量。
Figure PCTCN2019090303-appb-000016
是利用时间离散之后频散关系对有效频率计算的理论相移，其中ω是频率，Δt是时间步长。
继续进行逆傅里叶变换，得到新的第二震源函数，具体的变换公式为：
Figure PCTCN2019090303-appb-000017
步骤S1，基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量。
在本实施例中，第二震源时间函数作为波场数值模拟模型的输入，进 行时间步长迭代。步骤S1的其他具体实现与实施例1相同，不再赘述。
步骤S2，所有时间步长迭代结束，保存时间迭代所得的波场数据。
步骤S3，针对波场数据中需要输出的目的地震道进行逆时间频散变换操作，得到消除时间频散的地震道。
图9是本申请另一实施例提供的采用不同时间步长模拟得到的单道波形记录。
图9中的(a)为采用不同时间步长模拟，并进行空间滤波所得到的单道波形记录。其中(a)和(b)中上面五个图Δt1是伪谱法时间2阶差分格式，下面五个图Δt2是伪谱法时间4阶差分格式。可以看到即使Δt1＝Δt2＝7ms时，基本稳定，直到Δt1＝Δt2＝11ms时，才出现不稳定现象。但是时间频散会随着时间步长的不断增大而变大，计算精度随着时间步长的增大而降低。图9中的(b)是对图9中的(a)中的波形采用逆时间频散变换之后得到的波形记录，可以看出时间频散得到了有效的压制。从图9中可知，尤其是伪谱法时间4阶差分格式，在Δt2＝11ms时的计算结果精度仍然非常高。
步骤S4，存储逆时间频散变换处理之后的地震道。
在本实施例中，步骤S2、S3、S4具体实现与实施例1相同，不再赘述。
在本实施例中，由于在先对震源时间函数采用了正时间频散变换，使得逆时间频散变换更准确的修复了波场中的时间频散。
需要说明的是，以上参照附图所描述的各个实施例仅用以说明本申请而非限制本申请的范围，本领域的普通技术人员应当理解，在不脱离本申请的精神和范围的前提下对本申请进行的修改或者等同替换，均应涵盖在本申请的范围之内。此外，除上下文另有所指外，以单数形式出现的词包括复数形式，反之亦然。另外，除非特别说明，那么任何实施例的全部或一部分可结合任何其它实施例的全部或一部分来使用。

Claims (20)

  1. 一种拓展有限差分稳定性条件的波场模拟方法，其特征在于，所述方法包括：
    基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量；
    所有时间步长迭代结束，保存时间迭代所得的波场数据；
    针对所述波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据；
    存储所述的逆时间频散变换处理之后的波场数据。
  2. 根据权利要求1所述的方法，其特征在于，所述基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量之前，还包括：
    基于给定的时间步长和时间离散格式，针对第一震源时间函数进行正时间频散变换得到第二震源时间函数；其中，
    所述第二震源时间函数代替所述第一震源时间函数作为所述波场数值模拟模型的输入进行时间步长迭代。
  3. 根据权利要求1所述的方法，其特征在于，所述基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量，包括：
    基于波场数值模拟模型，进行时间步长迭代，计算出给定时间步长对应的满足CFL稳定性条件的最大波数阈值；
    所述时间步长迭代结束后，将空间域变换到波数域；
    采用低通滤波器滤除大于所述最大波数阈值的波场分量；
    所述滤波之后的波数域变换回所述空间域，开始下次的时间步长迭代。
  4. 根据权利要求3所述的方法，其特征在于，所述波场数值模拟模型的参数包括空间网格大小、模型最大速度、时间步长、时间和空间变量的离散格式；所述最大波数阈值对应的CFL稳定性条件为CFL稳定性条件的临界态。
  5. 根据权利要求3所述的方法，其特征在于，所述空间域变换到所述波数域采用傅立叶变换；所述波数域变换到所述空间域采用傅立叶逆变换。
  6. 根据权利要求3所述的方法，其特征在于，所述低通滤波器表达式为：
    Figure PCTCN2019090303-appb-100001
    其中，F(k)为低通滤波器，k是波数，所述低通滤波器的上限是最大波数阈值Kmax。
  7. 根据权利要求1所述的方法，其特征在于，所述逆时间频散变换操作的方法包括如下步骤：
    利用时间离散之后频散关系对有效频率计算实际相移θ(ω,Δt)，其中ω是频率，Δt是时间步长；
    应用变换：
    Figure PCTCN2019090303-appb-100002
    其中u(t)是需要输出的目的波场数据；
    进行逆傅里叶变换：
    Figure PCTCN2019090303-appb-100003
    得到消除时间频散的波场数据u′(t)。
  8. 一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行以下方法：
    基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量；
    所有时间步长迭代结束，保存时间迭代所得的波场数据；
    针对所述波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据；
    存储所述的逆时间频散变换处理之后的波场数据。
  9. 根据权利要求8所述的设备，其特征在于，所述基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量之前，还包括：
    基于给定的时间步长和时间离散格式，针对第一震源时间函数进行正时间频散变换得到第二震源时间函数；其中，
    所述第二震源时间函数代替所述第一震源时间函数作为所述波场数值模拟模型的输入进行时间步长迭代。
  10. 根据权利要求8所述的设备，其特征在于，所述基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量，包括：
    基于波场数值模拟模型，进行时间步长迭代，计算出给定时间步长对应的满足CFL稳定性条件的最大波数阈值；
    所述时间步长迭代结束后，将空间域变换到波数域；
    采用低通滤波器滤除大于所述最大波数阈值的波场分量；
    所述滤波之后的波数域变换回所述空间域，开始下次的时间步长迭代。
  11. 根据权利要求10所述的设备，其特征在于，所述波场数值模拟模型的参数包括空间网格大小、模型最大速度、时间步长、时间和空间变量的离散格式；所述最大波数阈值对应的CFL稳定性条件为CFL稳定性条件的临界态。
  12. 根据权利要求10所述的设备，其特征在于，所述空间域变换到所述波数域采用傅立叶变换；所述波数域变换到所述空间域采用傅立叶逆变换。
  13. 根据权利要求10所述的设备，其特征在于，所述低通滤波器表达式为：
    Figure PCTCN2019090303-appb-100004
    其中，F(k)为低通滤波器，k是波数，所述低通滤波器的上限是最大波数阈值Kmax。
  14. 根据权利要求8所述的设备，其特征在于，所述逆时间频散变换操作的方法包括如下步骤：
    利用时间离散之后频散关系对有效频率计算实际相移θ(ω,Δt)，其中ω是频率，Δt是时间步长；
    应用变换：
    Figure PCTCN2019090303-appb-100005
    其中u(t)是需要输出的目的波场数据；
    进行逆傅里叶变换：
    Figure PCTCN2019090303-appb-100006
    得到消除时间频散的波场数据u′(t)。
  15. 一种计算机可读存储介质，其上存储有处理器程序,其特征在于，该处理器执行以下方法：
    基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量；
    所有时间步长迭代结束，保存时间迭代所得的波场数据；
    针对所述波场数据中需要输出的目的波场数据进行逆时间频散变换操作，得到消除时间频散的波场数据；
    存储所述的逆时间频散变换处理之后的波场数据。
  16. 根据权利要求15所述的介质，其特征在于，所述基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量之前，还包括：
    基于给定的时间步长和时间离散格式，针对第一震源时间函数进行正时间频散变换得到第二震源时间函数；其中，
    所述第二震源时间函数代替所述第一震源时间函数作为所述波场数值模拟模型的输入进行时间步长迭代。
  17. 根据权利要求15所述的介质，其特征在于，所述基于波场数值模拟模型，进行时间步长迭代，采用空间滤波的方法滤除分布在高波数区域的不稳定波场分量，包括：
    基于波场数值模拟模型，进行时间步长迭代，计算出给定时间步长对应的满足CFL稳定性条件的最大波数阈值；
    所述时间步长迭代结束后，将空间域变换到波数域；
    采用低通滤波器滤除大于所述最大波数阈值的波场分量；
    所述滤波之后的波数域变换回所述空间域，开始下次的时间步长迭代。
  18. 根据权利要求17所述的介质，其特征在于，所述波场数值模拟模型的参数包括空间网格大小、模型最大速度、时间步长、时间和空间变量的离散格式；所述最大波数阈值对应的CFL稳定性条件为CFL稳定性条件的临界态；所述空间域变换到所述波数域采用傅立叶变换；所述波数域变换到所述空间域采用傅立叶逆变换。
  19. 根据权利要求17所述的介质，其特征在于，所述低通滤波器表达式为：
    Figure PCTCN2019090303-appb-100007
    其中，F(k)为低通滤波器，k是波数，所述低通滤波器的上限是最大波数阈值Kmax。
  20. 根据权利要求15所述的介质，其特征在于，所述逆时间频散变换操作的方法包括如下步骤：
    利用时间离散之后频散关系对有效频率计算实际相移θ(ω,Δt)，其中ω是频率，Δt是时间步长；
    应用变换：
    Figure PCTCN2019090303-appb-100008
    其中u(t)是需要输出的目的波场数据；
    进行逆傅里叶变换： 得到消除时间频散的波场数据u′(t)。
PCT/CN2019/090303 2018-06-08 2019-06-06 拓展有限差分稳定性条件的波场模拟方法、设备及介质 WO2019233475A1 (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN201980036883.2A CN112805598B (zh) 2018-06-08 2019-06-06 拓展有限差分稳定性条件的波场模拟方法、设备及介质
US16/973,230 US20210239870A1 (en) 2018-06-08 2019-06-06 Wave-field simulation method for extending finite-difference stability conditions, and apparatus and medium for implementing same

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201810585618.7A CN110579796A (zh) 2018-06-08 2018-06-08 拓展有限差分稳定性条件的波场模拟方法、装置及设备
CN201810585618.7 2018-06-08

Publications (1)

Publication Number Publication Date
WO2019233475A1 true WO2019233475A1 (zh) 2019-12-12

Family

ID=68770801

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2019/090303 WO2019233475A1 (zh) 2018-06-08 2019-06-06 拓展有限差分稳定性条件的波场模拟方法、设备及介质

Country Status (3)

Country Link
US (1) US20210239870A1 (zh)
CN (2) CN110579796A (zh)
WO (1) WO2019233475A1 (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115034103A (zh) * 2022-05-11 2022-09-09 中国石油大学(华东) 一种基于优化差分系数的正演模拟方法及其设备
CN115081267A (zh) * 2022-05-18 2022-09-20 内蒙古农业大学 基于优化的时空变步长有限差分地震波数值模拟方法
CN116719088A (zh) * 2023-05-30 2023-09-08 长安大学 一种航空瞬变电磁数据噪声压制方法

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111460593B (zh) * 2020-04-24 2023-03-31 安徽大学 一种空间域电磁分量确定方法及系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090204375A1 (en) * 2006-09-12 2009-08-13 Kikuji Hirose Numerical simulation apparatus for time dependent schrodinger equation
US20170336522A1 (en) * 2014-12-31 2017-11-23 Landmark Graphics Corporation Seismic Elastic Wave Simulation For Tilted Transversely Isotropic Media Using Adaptive Lebedev Staggered Grid
CN107942375A (zh) * 2017-11-17 2018-04-20 河海大学 一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10185046B2 (en) * 2014-06-09 2019-01-22 Exxonmobil Upstream Research Company Method for temporal dispersion correction for seismic simulation, RTM and FWI
CN106547023B (zh) * 2017-01-16 2017-11-28 青岛海洋地质研究所 一种精度高、计算稳定的复杂介质地震波场延拓方法
CN107976710B (zh) * 2017-11-17 2019-05-28 河海大学 一种基于声波方程的线性优化隐式时空域有限差分数值模拟方法
GB2572797B (en) * 2018-04-11 2020-08-12 Equinor Energy As Methods and systems for finite-difference wave equation modelling

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090204375A1 (en) * 2006-09-12 2009-08-13 Kikuji Hirose Numerical simulation apparatus for time dependent schrodinger equation
US20170336522A1 (en) * 2014-12-31 2017-11-23 Landmark Graphics Corporation Seismic Elastic Wave Simulation For Tilted Transversely Isotropic Media Using Adaptive Lebedev Staggered Grid
CN107942375A (zh) * 2017-11-17 2018-04-20 河海大学 一种基于声波方程的非线性优化隐式时空域有限差分数值模拟方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
GAO, YINGJIE ET AL.: "Long-term Numerical Simulation Based on Inverse Time Frequency Dispersion Transformation and Third-Order Symplectic Integral", ANNUAL MEETING OF CHINESE GEOSCIENCE UNION 2016, 15 October 2016 (2016-10-15), pages 1190 - 1192 *
HU , SHUHUA ET AL.: "Method for Simulating Pure qP Wave Field in Stable Complex TTI Medium", OIL GEOPHYSICAL PROSPECTING, vol. 53, no. 02, 15 April 2018 (2018-04-15), pages 280 - 286 *
ZHOU, XIYAN ET AL.: "3D Elastic Reverse Time Migration Based on P-and S-Wave Decoupling", CHINESE JOURNAL OF GEOPHYSICS, vol. 61, no. 03, 15 March 2018 (2018-03-15), pages 1038 - 1050 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115034103A (zh) * 2022-05-11 2022-09-09 中国石油大学(华东) 一种基于优化差分系数的正演模拟方法及其设备
CN115081267A (zh) * 2022-05-18 2022-09-20 内蒙古农业大学 基于优化的时空变步长有限差分地震波数值模拟方法
CN115081267B (zh) * 2022-05-18 2023-08-29 内蒙古农业大学 基于优化的时空变步长有限差分地震波数值模拟方法
CN116719088A (zh) * 2023-05-30 2023-09-08 长安大学 一种航空瞬变电磁数据噪声压制方法
CN116719088B (zh) * 2023-05-30 2024-05-14 长安大学 一种航空瞬变电磁数据噪声压制方法

Also Published As

Publication number Publication date
CN110579796A (zh) 2019-12-17
CN112805598A (zh) 2021-05-14
US20210239870A1 (en) 2021-08-05
CN112805598B (zh) 2023-10-31

Similar Documents

Publication Publication Date Title
WO2019233475A1 (zh) 拓展有限差分稳定性条件的波场模拟方法、设备及介质
Schötzau et al. A robust a-posteriori error estimator for discontinuous Galerkin methods for convection–diffusion equations
Sarris Extending the stability limit of the FDTD method with spatial filtering
WO2020038007A1 (zh) 拓展显式差分稳定性条件的波场模拟方法、装置及设备
Mohammadi et al. A Galerkin-reproducing kernel method: Application to the 2D nonlinear coupled Burgers' equations
San et al. A posteriori analysis of low-pass spatial filters for approximate deconvolution large eddy simulations of homogeneous incompressible flows
CN109490955B (zh) 一种基于规则网格的声波波动方程正演模拟方法及装置
Fillion-Gourdeau et al. A split-step numerical method for the time-dependent Dirac equation in 3-D axisymmetric geometry
Benner et al. Numerical linear algebra methods for linear differential-algebraic equations
Liu et al. Global stability for solutions to the exponential PDE describing epitaxial growth
Duru et al. On energy stable discontinuous Galerkin spectral element approximations of the perfectly matched layer for the wave equation
Rani et al. Numerical inversion of Laplace transform based on Bernstein operational matrix
Li et al. Engineering analysis error estimation when removing finite-sized features in nonlinear elliptic problems
Bernstein Optimal prediction of Burgers’s equation
Lyzell et al. Handling certain structure information in subspace identification
Zhang et al. Splitting algorithms for the high-order compact finite-difference schemes in wave-equation modeling
Li et al. Second‐order defeaturing error estimation for multiple boundary features
Aimi et al. BEM‐FEM coupling for the 1D Klein–Gordon equation
Bognar et al. Spectral method for time dependent Navier-Stokes equation
Salehi et al. The use of a Legendre pseudospectral viscosity technique to solve a class of nonlinear dynamic Hamilton–Jacobi equations
Lambers Krylov subspace spectral methods for the time-dependent Schrödinger equation with non-smooth potentials
Achleitner et al. Analysis and numerics of traveling waves for asymmetric fractional reaction-diffusion equations
Khan et al. An optimally scaled polynomial-Fourier-series method for the numerical solution of the Duffing oscillator
Duru et al. Stable and conservative time propagators for second order hyperbolic systems
Asiri et al. Cross-Convolution Approach for Delay Estimation in Fractional-Order Time-Delay Systems

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 19815870

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 19815870

Country of ref document: EP

Kind code of ref document: A1