WO2020057286A1 - 波场正演模拟方法及装置 - Google Patents

波场正演模拟方法及装置 Download PDF

Info

Publication number
WO2020057286A1
WO2020057286A1 PCT/CN2019/099890 CN2019099890W WO2020057286A1 WO 2020057286 A1 WO2020057286 A1 WO 2020057286A1 CN 2019099890 W CN2019099890 W CN 2019099890W WO 2020057286 A1 WO2020057286 A1 WO 2020057286A1
Authority
WO
WIPO (PCT)
Prior art keywords
finite difference
window function
coefficient
solution
function
Prior art date
Application number
PCT/CN2019/099890
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 中国科学院地质与地球物理研究所
Publication of WO2020057286A1 publication Critical patent/WO2020057286A1/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/30—Analysis
    • G01V1/303—Analysis for determining velocity profiles or travel times
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/622—Velocity, density or impedance
    • G01V2210/6222—Velocity; travel time

Definitions

  • the present invention relates to the technical field of data processing, and in particular, to a method and a device for wave field forward modeling.
  • the window function-based optimization method is very flexible.
  • the window function of the truncated pseudospectral method determines the accuracy of the finite difference operator.
  • various window functions are used to obtain the finite difference coefficient.
  • Another disadvantage of the finite difference method is that the calculation cost, variable time step size, and adaptive variable length space operator are two ways to reduce the calculation time from different aspects. However, these methods cannot fundamentally solve the problem of calculation speed.
  • Variable time step and adaptive variable length space operator are two methods to reduce the calculation time from different aspects.
  • these methods cannot fundamentally solve the problem of calculation speed, but GPU technology will bring considerable acceleration effects.
  • Some people have proposed using single GPU technology and window function optimization methods for elastic wave numerical simulation. But for large models or 3D models, single GPU is no longer applicable due to memory limitations.
  • the first type is to construct the objective function and use the optimization algorithm to solve it, but its disadvantage is that the objective function is too complicated (too many parameters), which limits the optimization Effect;
  • the second type is the window function algorithm, which has the disadvantage that the optimization effect of different window functions is difficult to quantitatively control, so that there is a problem of strong numerical dispersion when performing forward modeling of seismic wavefields.
  • An embodiment of the present invention provides a wave field forward simulation method to solve the technical problems that the objective function is complicated during the optimization process of the finite difference operator in the prior art, and that there is a strong numerical dispersion in the forward simulation of the seismic wave field.
  • the method includes:
  • the window function includes a finite difference coefficient and an adjustment coefficient
  • the wave field forward modeling was performed using the optimized uniform grid finite difference operator.
  • An embodiment of the present invention also provides a wave field forward simulation device to solve the technical problems that the objective function is complicated during the finite difference operator optimization process in the prior art, and that there is a strong numerical dispersion in the forward simulation of the seismic wave field.
  • the device includes:
  • a window function acquisition module for acquiring a window function, wherein the window function includes a finite difference coefficient and an adjustment coefficient;
  • An objective function creation module configured to create an objective function according to the maximum norm principle and the window function
  • a solution module for solving the objective function to obtain a solution of the finite difference coefficient and a solution of the tuning coefficient; and substituting the solution of the finite difference coefficient and the solution of the tuning coefficient into the window function ;
  • An optimization module for obtaining an optimized uniform grid finite difference operator by using the window function
  • the simulation module is used to perform forward modeling of the wave field by using the optimized uniform grid finite difference operator.
  • a new window function is proposed, and the window function includes a finite difference coefficient and a tuning coefficient, that is, the window function adds a tuning coefficient relative to a conventional window function in the prior art;
  • the maximum norm principle and the window function create the objective function, and through the curse of the window function, the coefficients to be optimized in the objective function are reduced; by solving the objective function, the solution of the finite difference coefficient and the solution of the tuning coefficient will be limited.
  • the solution of the difference coefficient and the tuning coefficient can be substituted into the window function to obtain an optimized uniform grid finite difference operator based on the window function.
  • the new window function adds the tuning coefficient, which makes it possible to combine the window Function and objective function to optimize the uniform-grid finite-difference operator.
  • the optimized uniform-grid finite-difference operator is beneficial to reduce the numerical dispersion when performing wave field forward simulation.
  • FIG. 1 is a flowchart of a wave field forward simulation method according to an embodiment of the present invention
  • FIG. 2 is a schematic diagram of a corresponding wave field snapshot of a point source according to an embodiment of the present invention
  • FIG. 3 is a single-channel comparison waveform diagram provided by an embodiment of the present invention.
  • FIG. 4 is a structural block diagram of a computer device according to an embodiment of the present invention.
  • FIG. 5 is a structural block diagram of a wave field forward simulation device according to an embodiment of the present invention.
  • a wave field forward simulation method As shown in FIG. 1, the method includes:
  • Step 101 Obtain a window function, where the window function includes a finite difference coefficient and an adjustment coefficient;
  • Step 102 Create an objective function according to the maximum norm principle and the window function.
  • Step 103 Solve the objective function to obtain a solution of the finite difference coefficient and a solution of the tuning coefficient;
  • Step 104 Substituting the solution of the finite difference coefficient and the solution of the tuning coefficient into the window function
  • Step 105 Use the window function to obtain an optimized uniform grid finite difference operator
  • Step 106 Use the optimized uniform grid finite difference operator to perform wave field forward modeling.
  • the window function includes a finite difference coefficient and an adjustment coefficient, that is, the window function is compared with a conventional window in the prior art.
  • the tuning coefficient is added to the function; further, the objective function is created based on the maximum norm principle and the window function.
  • the buffing of the window function reduces the coefficients to be optimized in the objective function.
  • the solution of the finite difference coefficient is obtained.
  • the tuning coefficient solution is substituted into the window function, an optimized uniform grid finite difference operator can be obtained based on the window function.
  • the tuning coefficient makes it possible to optimize the uniform-grid finite-difference operator by combining the window function and the objective function.
  • the optimized uniform-grid finite-difference operator is beneficial to reduce the numerical dispersion when performing wave field forward simulation. , which helps to obtain better simulation accuracy.
  • the uniform grid finite difference operator can be derived by the sinc function interpolation theory.
  • the sinc function can reconstruct the band-limited signal fn after uniform sampling (Diniz et al., 2012).
  • the sinc function can be expressed as:
  • ⁇ x is the spatial sampling interval
  • fn represents the sampled signal.
  • formula (2) and formula (3) can be expressed as formula (6) and formula (7).
  • f 0 represents a function of the intermediate position
  • f n is a function of the positive direction of the intermediate position
  • f -n is a function of the negative direction of the intermediate position
  • ⁇ represents Riemann's ⁇ function, which can be expressed as formula (8) and formula (9) after truncation.
  • n 1,2,3 ...., N / 2
  • w (n) represents the window function
  • c n the forward direction coefficient
  • c -n the reverse direction coefficient
  • the error function can be expressed as formula (12).
  • the error function can be expressed as formula (13).
  • the original binomial window function gives a traditional uniform-grid finite-difference operator, which is equivalent to the uniform-grid finite-difference operator obtained from Taylor series expansion. It retains the good accuracy of the low wavenumber part, but for High wave number components have no improvement.
  • the inventor of the present application proposes a new window function.
  • the new window function adds an adjustment coefficient to the window function in the prior art.
  • the adjustment coefficient may be multiple parameters or one.
  • Parameters for example, this application takes the tuning coefficient as two parameters as an example.
  • the new window function can be expressed as equation (15), from which we can find that in the new window function, we have added two parameters. They are parameters m and h, which can be called tuning coefficients.
  • an optimized uniform mesh finite difference operator with new parameters can be obtained.
  • the optimized uniform mesh finite difference operator is shown in formula (16):
  • bnew n represents the optimized uniform mesh finite difference operator.
  • the window function in the prior art only includes finite difference coefficients, so that the window function in the prior art cannot be combined with the objective function.
  • the new The window function includes a finite difference coefficient and an adjustment coefficient is also added, so that the new window function can be combined with the objective function.
  • This embodiment is based on the principle of maximum norm (for example, 1 norm or 2 norm) and the above-mentioned new The window function creates an objective function and obtains a solution of the finite difference coefficient and a solution of the tuning coefficient by solving the objective function.
  • the objective function is shown in formula (17):
  • k x represents the wave number range
  • ⁇ x represents the space sampling interval
  • T represents the maximum error allowed during optimization.
  • a simulated annealing algorithm is used to solve the above objective function.
  • the optimized uniform grid finite difference operator is obtained.
  • the optimization algorithm and window function optimization algorithm are used to optimize the uniform grid finite difference operator. This algorithm is the first time in the industry, which makes it possible to optimize the uniform grid finite difference. The operator can obtain the characteristics of the window function algorithm and the optimization algorithm at the same time, and can obtain better simulation accuracy in the practical stage.
  • a 600 ⁇ 600 two-dimensional uniform model with a grid spacing of 5m is established, but the model size actually calculated is 700 ⁇ 700, where each side contains 50 absorption boundaries.
  • the P-wave velocity is 2000 m / s
  • the S-wave velocity is 1400 m / s
  • the density is 1000.
  • the main frequency of the Ricker wavelet is 50Hz, which is located in the center of the speed model.
  • the elastic wave equation in a two-dimensional heterogeneous medium is shown in formula (18):
  • Figure 2 shows a snapshot of the corresponding wave field of the point source.
  • the data of the dotted line in Figure 2 is extracted, and a 48-order finite difference operator is used as the reference solution.
  • the dotted line is the data extracted in Fig. 2 processed using the uniform grid finite difference operator optimized by the present application, and the solid line is the reference solution.
  • (A) to (d) in Fig. 3 correspond to the graph
  • the processed data of the dotted waveforms in (a) to (d) of 2 and (e) to (f) in FIG. 3 show that the data extracted in (a) to (d) of FIG. 2 are uniform after the optimization of this application.
  • the data processed by the optimized uniform grid finite difference operator in this application has less numerical dispersion than the traditional method and the improved binomial window.
  • a computer device includes a memory 402, a processor 404, and a computer program stored in the memory and executable on the processor.
  • the processor executes the computer program.
  • obtain a window function wherein the window function includes finite difference coefficients and tuning coefficients; create an objective function according to the maximum norm principle and the window function; solve the objective function to obtain the finite difference The solution of the coefficient and the tuning coefficient; substituting the solution of the finite difference coefficient and the solution of the tuning coefficient into the window function; using the window function to obtain an optimized uniform grid finite difference calculation
  • the wave field forward modeling is performed using the optimized uniform grid finite difference operator.
  • the processor when the processor executes the computer program, the processor further performs the following operations:
  • the simulated annealing algorithm is used to solve the objective function.
  • the computer equipment may be a computer terminal, a server, or a similar computing device.
  • a computer-readable storage medium is provided.
  • the computer program stored on the computer-readable storage medium is executed, at least the following operations are performed: obtaining a window function, where the window function includes a finite difference coefficient and an integer. Tuning coefficient; creating an objective function according to the maximum norm principle and the window function; solving the objective function to obtain a solution of the finite difference coefficient and a solution of the tuning coefficient; and summing the solution of the finite difference coefficient
  • the solution of the tuning coefficient is substituted into the window function; the optimized uniform grid finite difference operator is obtained by using the window function; and the wave field forward simulation is performed by using the optimized uniform grid finite difference operator.
  • the computer program stored on the computer-readable storage medium when executed, it also performs the following operation: using a simulated annealing algorithm to solve the objective function.
  • the computer-readable storage medium includes permanent and non-permanent, removable and non-removable media, and information storage can be achieved by any method or technology.
  • Information may be computer-readable instructions, data structures, modules of a program, or other data.
  • Examples of computer-readable storage media include, but are not limited to, phase change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only Memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, read-only disc read-only memory (CD-ROM), digital versatile disc (DVD), or other optical storage , Magnetic tape cartridges, magnetic tape storage or other magnetic storage devices or any other non-transmission media can be used to store information that can be accessed by computing devices.
  • computer-readable storage media does not include temporary computer-readable media, such as modulated data signals and carrier waves.
  • a wave field forward simulation device is also provided in the embodiments of the present invention, as described in the following embodiments. Since the principle of the wave field forward simulation device solves the problem similar to the wave field forward simulation method, the implementation of the wave field forward simulation device can refer to the implementation of the wave field forward simulation method.
  • the term "unit” or “module” may be a combination of software and / or hardware that realizes a predetermined function.
  • the devices described in the following embodiments are preferably implemented in software, implementation in hardware, or a combination of software and hardware is also possible and conceived.
  • FIG. 5 is a structural block diagram of a wave field forward simulation device according to an embodiment of the present invention. As shown in FIG. 5, the device includes:
  • a window function obtaining module 501 is configured to obtain a window function, where the window function includes a finite difference coefficient and an adjustment coefficient;
  • a solving module 503 is configured to solve the objective function to obtain a solution of the finite difference coefficient and a solution of the tuning coefficient; substitute the solution of the finite difference coefficient and the solution of the tuning coefficient into the window function;
  • the simulation module 505 is configured to perform forward modeling of the wave field by using the optimized uniform grid finite difference operator.
  • the expression of the window function is:
  • the expression of the optimized uniform grid finite difference operator is:
  • bnew n represents the optimized uniform mesh finite difference operator.
  • the expression of the objective function is:
  • k x represents the wave number range
  • ⁇ x represents the space sampling interval
  • T represents the maximum error allowed during optimization.
  • the solving module is specifically configured to solve the objective function by using a simulated annealing algorithm.
  • a new window function is proposed, and the window function includes a finite difference coefficient and a tuning coefficient, that is, the window function adds a tuning coefficient to a conventional window function in the prior art.
  • the objective function is created, and the coefficients to be optimized in the objective function are reduced by the curse of the window function; the solution of the finite difference coefficient and the solution of the tuning coefficient are obtained by solving the objective function
  • an optimized uniform grid finite difference operator can be obtained based on the window function. Because the new window function adds tuning coefficients, it achieves The window function and objective function can be combined to optimize the uniform grid finite difference operator.
  • the optimized uniform grid finite difference operator is beneficial to reduce the numerical dispersion when performing wave field forward simulation.
  • the embodiments of the present invention may be provided as a method, a system, or a computer program product. Therefore, the present invention may take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Moreover, the present invention may take the form of a computer program product implemented on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.
  • computer-usable storage media including, but not limited to, disk storage, CD-ROM, optical storage, etc.
  • These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to work in a particular manner such that the instructions stored in the computer-readable memory produce a manufactured article including an instruction device, the instructions
  • the device implements the functions specified in one or more flowcharts and / or one or more blocks of the block diagram.
  • These computer program instructions can also be loaded on a computer or other programmable data processing device, so that a series of steps can be performed on the computer or other programmable device to produce a computer-implemented process, which can be executed on the computer or other programmable device.
  • the instructions provide steps for implementing the functions specified in one or more flowcharts and / or one or more blocks of the block diagrams.

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

一种波场正演模拟方法及装置，其中，该方法包括：获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；根据最大范数原理和所述窗函数，创建目标函数；求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；将所述最优解代入所述窗函数；利用所述窗函数得到优化后的均匀网格有限差分算子；利用优化后的均匀网格有限差分算子进行波场正演模拟。该方案优化了目标函数中待优化的系数，首次实现了可以结合窗函数和目标函数的方式来优化均匀网格有限差分算子，该优化后的均匀网格有限差分算子在进行波场正演模拟时有利于减少数值频散。

Description

波场正演模拟方法及装置 技术领域
本发明涉及数据处理技术领域，特别涉及一种波场正演模拟方法及装置。
背景技术
基于有限差分格式的地震波场正演模拟在逆时偏移中得到了广泛的应用，全波反演、逆时偏移都需要进行多次正演和波场的反传，因此，正演模拟的精度和速度是非常重要的考量因素，其对成像和反演具有重要意义。许多学者对波动方程的有限差分格式进行了广泛的研究。为了减少数值频散，已经开发了多种有限差分格式，例如，可变网格、不规则网格、交错网格和旋转交错网格等。
然而，空间导数上的高阶有限差分算子的经典系数通常是由空间导数项的泰勒级数展开确定的。仅改变有限差分格式不能最小化数值误差。如果传统的有限差分系数用于地震波场正演模拟，强的数值频散是不可避免的，特别是对于波动方程中的高波数范围。在最近的二十年中，学者们已经提出了许多优化方法，如牛顿方法、隐格式、基于时空域频散关系的方法、最小二乘法等。但是这类算法使用的目标函数参数过多，需要对与有限差分算子阶数相同数量的参数进行优化，对优化算法的要求提出了更高的要求。与上述方法相比，基于窗函数的优化方法是非常灵活的，我们认为，截断伪谱方法的窗函数决定了有限差分算子的精度，如今各种窗函数被用来得到有限差分系数。有限差分法的另一缺点是计算成本、可变时间步长和自适应变长空间算子是从不同方面减少计算时间的两种方式。然而，这些方法不能从根本上解决计算速度的问题。
上世纪60年代末期，有人提出将二阶显格式的有限差分技术应用于层状介质的弹性波数值模拟，拉开了有限差分深入研究的序幕。进而发展了适应非均匀介质的有限差分格式以进行弹性波数值模拟，还发展了高阶数的有限差分格式进行声波方程求解。有人将时间域有限差分方法应用于粘声波方程，并进行了基于各向异性介质的逆时偏移(RTM)。高阶有限差分算子在空间导数上的传统系数通常由泰勒级数展开确定。采用常规的有限差分系数进行地震波场正演模拟，会出现较强的数值频散，特别是波动方程中的高波数范围。近二十年来，不同学者提出了牛顿法、隐式格式、模拟退火算法、最小二乘法。但是这些优化方法实施起来非常复杂，并且由于目标函数的影响而不能广泛 用于偏移和反演。与上述方法相比，基于窗口函数的优化方法实现起来非常灵活，有限差分法是伪谱法空间卷积序列的空间截断形式。不同学者利用不同窗函数得到有限差分系数。有限差分法的另一个缺点是计算成本，可变时间步长和自适应可变长度空间算子，是从不同方面减少计算时间的两种方法。然而，这些方法不能从根本上解决计算速度问题，但是GPU技术会带来相当大的加速效应，有人提出使用单GPU技术和窗函数优化方法进行弹性波数值模拟。但是对于大型模型或者三维模型，由于内存限制，单GPU不再适用。
综上，目前主流的有限差分算子优化算法有两大类，第一类是构建目标函数，利用最优化算法进行求解，但其缺点是目标函数过于复杂(参数过多)，限制了优化的效果；第二类是窗函数算法，其缺点是不同窗函数的优化效果很难定量控制，使得在进行地震波场正演模拟时存在较强的数值频散的问题。
发明内容
本发明实施例提供了一种波场正演模拟方法，以解决现有技术中有限差分算子优化过程中目标函数复杂、地震波场正演模拟时存在较强的数值频散的技术问题。该方法包括：
获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；
根据最大范数原理和所述窗函数，创建目标函数；
求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；
将所述有限差分系数的解和所述整调系数的解，代入所述窗函数；
利用所述窗函数得到优化后的均匀网格有限差分算子；
利用优化后的均匀网格有限差分算子进行波场正演模拟。
本发明实施例还提供了一种波场正演模拟装置，以解决现有技术中有限差分算子优化过程中目标函数复杂、地震波场正演模拟时存在较强的数值频散的技术问题。该装置包括：
窗函数获取模块，用于获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；
目标函数创建模块，用于根据最大范数原理和所述窗函数，创建目标函数；
求解模块，用于求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；将所述有限差分系数的解和所述整调系数的解，代入所述窗函数；
优化模块，用于利用所述窗函数得到优化后的均匀网格有限差分算子；
模拟模块，用于利用优化后的均匀网格有限差分算子进行波场正演模拟。
在本发明实施例中，提出了一种新的窗函数，该窗函数包括有限差分系数和整调系数，即该窗函数相对于现有技术中传统的窗函数添加了整调系数；进而基于最大范数原理和窗函数，创建目标函数，通过窗函数的加持，使得降低了目标函数中待优化的系数；通过求解目标函数，来得到有限差分系数的解和整调系数的解，将有限差分系数的解和整调系数的解，代入窗函数后，即可基于窗函数得到优化后的均匀网格有限差分算子，由于新的窗函数添加了整调系数，使得实现了可以结合窗函数和目标函数的方式来优化均匀网格有限差分算子，该优化后的均匀网格有限差分算子在进行波场正演模拟时有利于减少数值频散。
附图说明
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。在附图中：
图1是本发明实施例提供的一种波场正演模拟方法的流程图；
图2是本发明实施例提供的一种点源相应波场快照的示意图；
图3是本发明实施例提供的一种单道对比波形图；
图4是本发明实施例提供的一种计算机设备的结构框图；
图5是本发明实施例提供的一种波场正演模拟装置的结构框图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例做进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。
在本发明实施例中，提供了一种波场正演模拟方法，如图1所示，该方法包括：
步骤101：获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；
步骤102：根据最大范数原理和所述窗函数，创建目标函数；
步骤103：求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；
步骤104：将所述有限差分系数的解和所述整调系数的解，代入所述窗函数；
步骤105：利用所述窗函数得到优化后的均匀网格有限差分算子；
步骤106：利用优化后的均匀网格有限差分算子进行波场正演模拟。
由图1所示的流程可知，在本发明实施例中，提出了一种新的窗函数，该窗函数包括有限差分系数和整调系数，即该窗函数相对于现有技术中传统的窗函数添加了整调系数；进而基于最大范数原理和窗函数，创建目标函数，通过窗函数的加持，使得降低了目标函数中待优化的系数；通过求解目标函数，来得到有限差分系数的解和整调系数的解，将有限差分系数的解和整调系数的解，代入窗函数后，即可基于窗函数得到优化后的均匀网格有限差分算子，由于新的窗函数添加了整调系数，使得实现了可以结合窗函数和目标函数的方式来优化均匀网格有限差分算子，该优化后的均匀网格有限差分算子在进行波场正演模拟时有利于减少数值频散，有利于获得更好的模拟精度。
本申请发明人发现，现有技术中通过以下方式获得均匀网格有限差分算子：
均匀网格有限差分算子可以通过sinc函数插值理论推导得到，sinc函数可以重构均匀采样后的带限信号fn(Diniz et al.,2012)，sinc函数可以表示为：
Figure PCTCN2019099890-appb-000001
其中， x表示采样点的位置；Δx为空间采样间隔，
Figure PCTCN2019099890-appb-000002
代表截止波数， fn表示采样信号。利用窗函数截断公式(1)的一阶导数和二阶导数，获取了均匀网格有限差分算子，均匀网格有限差分算子求导后可表示为公式(2)和公式(3)，其中，f(n)表示原始信号。
Figure PCTCN2019099890-appb-000003
Figure PCTCN2019099890-appb-000004
利用窗函数
Figure PCTCN2019099890-appb-000005
截断公式(2)和(3)后，可以分别表示为公式(4)和公式(5)。
Figure PCTCN2019099890-appb-000006
Figure PCTCN2019099890-appb-000007
因n＝0为公式(2)和公式(3)的一个奇异点，根据信号采样理论，公式(2)和公式(3)可以表示为公式(6)和公式(7)。
Figure PCTCN2019099890-appb-000008
Figure PCTCN2019099890-appb-000009
其中，f 0表示中间位置的函数，f n为中间位置的正方向的函数，f -n为中间位置的负方向的函数，
Figure PCTCN2019099890-appb-000010
和
Figure PCTCN2019099890-appb-000011
为f n的系数，
Figure PCTCN2019099890-appb-000012
n＝1,2,3....,∞，
Figure PCTCN2019099890-appb-000013
这里ξ代表黎曼ξ函数，加窗截断后可以表示为公式(8)和公式(9)。
Figure PCTCN2019099890-appb-000014
Figure PCTCN2019099890-appb-000015
这里，
Figure PCTCN2019099890-appb-000016
n＝1,2,3....,N/2，
Figure PCTCN2019099890-appb-000017
w(n)表示窗函数，c n表示正方向系数，c -n表示反方向系数。
将公式(8)和公式(9)变换到频率域后分别为公式(10)和公式(11)
Figure PCTCN2019099890-appb-000018
Figure PCTCN2019099890-appb-000019
对于一阶导数，可以将误差函数表示成公式(12)。对于二阶导数，可以将误差函数表示成公式(13)。
Figure PCTCN2019099890-appb-000020
Figure PCTCN2019099890-appb-000021
可见，通过定义不同的窗函数可以得到均匀网格有限差分算子。还可以基于传统二项窗函数定义一个改进二项式窗函数，常规的二项式窗口函数可以表示为方程式(14)，该二项式窗函数改进过程中认为可以用N+M代替N，其中N是阶数，M是偶数。这种方法的最大缺点是不能减少在有限差分阶数较低的情况下的数值频散。
Figure PCTCN2019099890-appb-000022
原始的二项式窗函数给出了传统的均匀网格有限差分算子，相当于从泰勒级数展开式得到的均匀网格有限差分算子，它保留了低波数部分的良好精度，但是对于高波数分量却没有改进效果。
本申请发明人提出了一种新的窗函数，该新的窗函数相对于现有技术中的窗函数添加了整调系数，具体的，该整调系数可以是多个参数也可以是1个参数，例如，本申请以整调系数为2个参数为例，该新的窗函数可以表示为等式(15)，从中我们可以发现，在新的窗函数中，我们增加了两个参数，分别为参数m和h，可以称之为整调系数。
Figure PCTCN2019099890-appb-000023
其中，
Figure PCTCN2019099890-appb-000024
表示窗函数； m和h表示整调系数； n表示窗函数空间点位；N表示有限差分阶数。
基于新的窗函数可以得到包含新参数的优化后的均匀网格有限差分算子，例如，优化后的均匀网格有限差分算子如公式(16)所示：
Figure PCTCN2019099890-appb-000025
其中，bnew n表示优化后的均匀网格有限差分算子。
本申请发明人发现，由于现有技术中的窗函数中的未知数只包含了有限差分系数，使得现有技术中的窗函数无法与目标函数结合，本申请提出新的窗函数之后，该新的窗函数包含有限差分系数的同时还添加了整调系数，使得新的窗函数可以与目标函数进行 结合，本实施例根据最大范数原理(例如，1范数或2范数)和上述新的窗函数创建目标函数，并通过求解目标函数，来得到有限差分系数的解和整调系数的解。该目标函数如公式(17)所示：
Figure PCTCN2019099890-appb-000026
其中，k x表示波数范围；Δx表示空间采样间隔；T表示优化时允许的最大误差。
具体实施时，为了获得更好的求解精度，在本实施例中，采用模拟退火算法求解上述目标函数，求得的有限差分算子如下表1所示，其中，整调系数m＝32.3，h＝15.6。
表1
10阶 12阶 16阶 20阶
0.8571428571 0.8888888889 0.9090909090 1.257430489308
-0.2678571429 -0.311111111 -0.3409090909 -0.126377801689
  0.1131313131 0.13986013986 0.037148969395
  -0.03535353534 -0.05244755244 -0.013909755048
    0.01678321678 0.005509987784
    -0.00437062937 -0.002130526526
      0.000763273991
      -0.000241112023
      0.0000626766642
      -0.000011537920
解得有限差分系数的解和整调系数的解之后，将有限差分系数的解和整调系数的解代入上述公式(15)中的窗函数，再将窗函数代入公式(16)，即可得到优化后的均匀网格有限差分算子。在该优化均匀网格有限差分算子的过程，同时利用了最优化算法和窗函数优化算法对均匀网格有限差分算子进行优化，该算法在业内尚属首次，使得优化均匀网格有限差分算子可以同时获得窗函数算法和最优化算法的特点，在实用阶段可以获得更好的模拟精度。
以下结合示例说明上述波场正演模拟方法的优点。
例如，建立一个600×600的、网格间距为5m的二维均匀模型，但实际计算的模型尺寸为700×700，其中，每边包含50个吸收边界。纵波速度为2000m/s，横波速度为1400m/s，密度为1000。Ricker子波的主频率为50Hz，位于速度模型的中心。二维非均匀介质中的弹性波方程如公式(18)所示：
Figure PCTCN2019099890-appb-000027
图2显示了点源相应波场的快照，抽取图2中虚线部位的数据，使用48阶有限差分算子作为参考解。在图3中，虚线为图2抽取的数据使用本申请优化后的均匀网格有限差分算子处理的数据，实线为参考解，图3中的(a)至(d)分别对应了图2的(a)至(d)中虚线波形的处理数据，图3中的(e)至(f)显示了图2的(a)至(d)中抽取的数据采用本申请优化后的均匀网格有限差分算子处理的数据与参考解的差值。显然，本申请中优化后的均匀网格有限差分算子处理的数据比传统方法和改进二项式窗有着更少的数值频散。当我们将这些方法与参考解比较时，可以注意到，优化后的均匀网格有限差分算子处理的数据更接近参考解，这也证明了上述误差分析的结论。
在本实施例中提供了一种计算机设备，如图4所示，包括存储器402、处理器404及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时执行以下操作：获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；根据最大范数原理和所述窗函数，创建目标函数；求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；将所述有限差分系数的解和所述整调系数的解，代入所述窗函数；利用所述窗函数得到优化后的均匀网格有限差分算子；利用优化后的均匀网格有限差分算子进行波场正演模拟。
在一个实施例中，所述处理器执行所述计算机程序时还执行以下操作：
采用模拟退火算法求解所述目标函数。
具体的，该计算机设备可以是计算机终端、服务器或者类似的运算装置。
在本实施例中提供了一种计算机可读存储介质，所述计算机可读存储介质存储的计算机程序被执行时至少执行以下操作：获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；根据最大范数原理和所述窗函数，创建目标函数；求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；将所述有限差分系数的解和所述整调系数的解，代入所述窗函数；利用所述窗函数得到优化后的均匀网格有限差分算子；利用优化后的均匀网格有限差分算子进行波场正演模拟。
在一个实施例中，所述计算机可读存储介质存储的计算机程序被执行时还执行以下操作：采用模拟退火算法求解所述目标函数。
具体的，计算机可读存储介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机可读存储介质的例子包括，但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带，磁带磁磁盘存储或其他磁性存储设备或任何其他非传输介质，可用于存储可以被计算设备访问的信息。按照本文中的界定，计算机可读存储介质不包括暂存电脑可读媒体(transitory media)，如调制的数据信号和载波。
基于同一发明构思，本发明实施例中还提供了一种波场正演模拟装置，如下面的实施例所述。由于波场正演模拟装置解决问题的原理与波场正演模拟方法相似，因此波场正演模拟装置的实施可以参见波场正演模拟方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。
图5是本发明实施例的波场正演模拟装置的一种结构框图，如图5所示，该装置包括：
窗函数获取模块501，用于获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；
目标函数创建模块502，用于根据最大范数原理和所述窗函数，创建目标函数；
求解模块503，用于求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；将所述有限差分系数的解和所述整调系数的解，代入所述窗函数；
优化模块504，用于利用所述窗函数得到优化后的均匀网格有限差分算子；
模拟模块505，用于利用优化后的均匀网格有限差分算子进行波场正演模拟。
在一个实施例中，所述窗函数的表达式为：
Figure PCTCN2019099890-appb-000028
其中，
Figure PCTCN2019099890-appb-000029
表示窗函数； m和h表示整调系数； n表示窗函数空间点位；N表示有限差分阶数。
在一个实施例中，优化后的均匀网格有限差分算子的表达式为：
Figure PCTCN2019099890-appb-000030
其中，bnew n表示优化后的均匀网格有限差分算子。
在一个实施例中，所述目标函数的表达式为：
Figure PCTCN2019099890-appb-000031
其中，k x表示波数范围；Δx表示空间采样间隔；T表示优化时允许的最大误差。
在一个实施例中，所述求解模块具体用于采用模拟退火算法求解所述目标函数。
本发明实施例实现了如下技术效果：提出了一种新的窗函数，该窗函数包括有限差分系数和整调系数，即该窗函数相对于现有技术中传统的窗函数添加了整调系数；进而基于最大范数原理和窗函数，创建目标函数，通过窗函数的加持，使得降低了目标函数中待优化的系数；通过求解目标函数，来得到有限差分系数的解和整调系数的解，将有限差分系数的解和整调系数的解，代入窗函数后，即可基于窗函数得到优化后的均匀网格有限差分算子，由于新的窗函数添加了整调系数，使得实现了可以结合窗函数和目标函数的方式来优化均匀网格有限差分算子，该优化后的均匀网格有限差分算子在进行波场正演模拟时有利于减少数值频散。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产 生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (12)

  1. 一种波场正演模拟方法，其特征在于，包括：
    获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；
    根据最大范数原理和所述窗函数，创建目标函数；
    求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；
    将所述有限差分系数的解和所述整调系数的解，代入所述窗函数；
    利用所述窗函数得到优化后的均匀网格有限差分算子；
    利用优化后的均匀网格有限差分算子进行波场正演模拟。
  2. 如权利要求1所述的波场正演模拟方法，其特征在于，所述窗函数的表达式为：
    Figure PCTCN2019099890-appb-100001
    其中，
    Figure PCTCN2019099890-appb-100002
    表示窗函数；m和h表示整调系数；n表示窗函数空间点位；N表示有限差分阶数。
  3. 如权利要求2所述的波场正演模拟方法，其特征在于，优化后的均匀网格有限差分算子的表达式为：
    Figure PCTCN2019099890-appb-100003
    其中，bnew n表示优化后的均匀网格有限差分算子。
  4. 如权利要求3所述的波场正演模拟方法，其特征在于，所述目标函数的表达式为：
    Figure PCTCN2019099890-appb-100004
    其中，k x表示波数范围；Δx表示空间采样间隔；T表示优化时允许的最大误差。
  5. 如权利要求1至4中任一项所述的波场正演模拟方法，其特征在于，求解所述目标函数，包括：
    采用模拟退火算法求解所述目标函数。
  6. 一种波场正演模拟装置，其特征在于，包括：
    窗函数获取模块，用于获取窗函数，其中，所述窗函数包括有限差分系数和整调系数；
    目标函数创建模块，用于根据最大范数原理和所述窗函数，创建目标函数；
    求解模块，用于求解所述目标函数，得到所述有限差分系数的解和所述整调系数的解；将所述有限差分系数的解和所述整调系数的解，代入所述窗函数；
    优化模块，用于利用所述窗函数得到优化后的均匀网格有限差分算子；
    模拟模块，用于利用优化后的均匀网格有限差分算子进行波场正演模拟。
  7. 如权利要求6所述的波场正演模拟装置，其特征在于，所述窗函数的表达式为：
    Figure PCTCN2019099890-appb-100005
    其中，
    Figure PCTCN2019099890-appb-100006
    表示窗函数；m和h表示整调系数；n表示窗函数空间点位；N表示有限差分阶数。
  8. 如权利要求7所述的波场正演模拟装置，其特征在于，优化后的均匀网格有限差分算子的表达式为：
    Figure PCTCN2019099890-appb-100007
    其中，bnew n表示优化后的均匀网格有限差分算子。
  9. 如权利要求8所述的波场正演模拟装置，其特征在于，所述目标函数的表达式为：
    Figure PCTCN2019099890-appb-100008
    其中，k x表示波数范围；Δx表示空间采样间隔；T表示优化时允许的最大误差。
  10. 如权利要求6至9中任一项所述的波场正演模拟装置，其特征在于，所述求解模块具体用于采用模拟退火算法求解所述目标函数。
  11. 一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至5任一项所述的波场正演模拟方法。
  12. 一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有执行权利要求1至5任一项所述的波场正演模拟方法的计算机程序。
PCT/CN2019/099890 2018-09-20 2019-08-09 波场正演模拟方法及装置 WO2020057286A1 (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201811117484.2A CN109490954B (zh) 2018-09-20 2018-09-20 波场正演模拟方法及装置
CN201811117484.2 2018-09-20

Publications (1)

Publication Number Publication Date
WO2020057286A1 true WO2020057286A1 (zh) 2020-03-26

Family

ID=65690555

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2019/099890 WO2020057286A1 (zh) 2018-09-20 2019-08-09 波场正演模拟方法及装置

Country Status (2)

Country Link
CN (1) CN109490954B (zh)
WO (1) WO2020057286A1 (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109490954B (zh) * 2018-09-20 2019-12-20 中国科学院地质与地球物理研究所 波场正演模拟方法及装置
CN116578825A (zh) * 2022-12-28 2023-08-11 上海勘测设计研究院有限公司 气象预测误差修正方法、装置、介质及电子设备

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103499835A (zh) * 2013-10-13 2014-01-08 中国石油集团西北地质研究所 初至波波形反演近地表速度模型方法
US20150362610A1 (en) * 2014-06-16 2015-12-17 Petroleum Geo-Services Inc. Method of suppressing spectral artefacts of wavefield decomposition caused by imperfect extrapolation
CN106033125A (zh) * 2016-06-29 2016-10-19 中国石油化工股份有限公司 压制叠前大角度道集干涉的提频方法
CN106842306A (zh) * 2017-04-18 2017-06-13 中国科学院地质与地球物理研究所 一种全局优化的交错网格有限差分正演模拟方法和装置
CN109490954A (zh) * 2018-09-20 2019-03-19 中国科学院地质与地球物理研究所 波场正演模拟方法及装置

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103853930B (zh) * 2014-03-19 2017-01-18 中国科学院地质与地球物理研究所 地震矢量波场数值模拟方法和装置
CN106353797A (zh) * 2015-07-17 2017-01-25 中国石油化工股份有限公司 一种高精度地震正演模拟方法
CN108196303B (zh) * 2017-12-29 2019-10-01 中国石油天然气集团公司 弹性波波场分离方法、装置、存储介质及设备

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103499835A (zh) * 2013-10-13 2014-01-08 中国石油集团西北地质研究所 初至波波形反演近地表速度模型方法
US20150362610A1 (en) * 2014-06-16 2015-12-17 Petroleum Geo-Services Inc. Method of suppressing spectral artefacts of wavefield decomposition caused by imperfect extrapolation
CN106033125A (zh) * 2016-06-29 2016-10-19 中国石油化工股份有限公司 压制叠前大角度道集干涉的提频方法
CN106842306A (zh) * 2017-04-18 2017-06-13 中国科学院地质与地球物理研究所 一种全局优化的交错网格有限差分正演模拟方法和装置
CN109490954A (zh) * 2018-09-20 2019-03-19 中国科学院地质与地球物理研究所 波场正演模拟方法及装置

Also Published As

Publication number Publication date
CN109490954A (zh) 2019-03-19
CN109490954B (zh) 2019-12-20

Similar Documents

Publication Publication Date Title
Etgen et al. The pseudo-analytical method: Application of pseudo-Laplacians to acoustic and acoustic anisotropic wave propagation
AU2012231651B2 (en) System and method for seismic data modeling and migration
CN104977607B (zh) 利用变步长网格声波波场模拟的时间域全波形反演方法
WO2020057286A1 (zh) 波场正演模拟方法及装置
CN110471113B (zh) 基于非稳态地震资料的反演动校正方法、装置及存储介质
CN110109177B (zh) 基于旋转时空双变网格有限差分法的地震波正演模拟方法
CN109635513B (zh) 三维TTI介质qP波波场模拟方法及装置
CN109164488B (zh) 一种梯形网格有限差分地震波场模拟方法
CN108828659B (zh) 基于傅里叶有限差分低秩分解的地震波场延拓方法及装置
CN107561588B (zh) 一种地震数据噪声压制方法及装置
WO2020063131A1 (zh) 拟微分算子储存方法及装置
US11686870B2 (en) Interpretive-guided velocity modeling seismic imaging method and system, medium and device
CN115270579A (zh) 二阶声波方程有限差分数值模拟参数选取方法
US20150142831A1 (en) Computerized method and a computer program rpoduct for determining a resulting data set representative of a geological region of interest
Kim et al. Least-squares reverse time migration using analytic-signal-based wavefield decomposition
CN112906267A (zh) 频率域粘弹性波正演模拟方法、装置、设备及存储介质
CN104977608B (zh) 利用固定网格声波波场模拟的时间域全波形反演方法
Shao‐Lin et al. Symplectic RKN schemes for seismic scalar wave simulations
CN115184986B (zh) 不依赖震源的全局包络互相关全波形反演方法
Pan et al. A fast minimum weighted norm interpolation
CN117077480A (zh) 一种分数阶粘声波动方程数值模拟方法、系统及设备
Nickerson Developments in minimum entropy deconvolution
CN112649874A (zh) 基于低秩分解的粘声逆时偏移方法及系统
CN112379422A (zh) 垂变网格地震波场外推方法及装置
CN117607951A (zh) 一种适用于三维隧道模型的多尺度全波形反演方法及系统

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: 19862741

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: 19862741

Country of ref document: EP

Kind code of ref document: A1