CN103901472B

CN103901472B - Frequency domain forward modeling method and device

Info

Publication number: CN103901472B
Application number: CN201410125866.5A
Authority: CN
Inventors: 唐祥德; 刘洪�
Original assignee: Institute of Geology and Geophysics of CAS
Current assignee: Institute of Geology and Geophysics of CAS
Priority date: 2014-03-31
Filing date: 2014-03-31
Publication date: 2015-03-11
Anticipated expiration: 2034-03-31
Also published as: CN103901472A

Abstract

The invention relates to the technical field of geological exploration, in particular to a frequency domain forward modeling method and device. The methods include establishing a difference formula in 17-point format and constructing a sparse matrix; reading in wavelet parameters and velocity models; frequency cycling to obtain a single-frequency wave field; inverse Fourier transform of all frequency wave fields; and obtaining forward modeling results. Through the method and device of the embodiment of the present invention, the fourth-order high-precision forward modeling overcomes the previous condition that dx must be equal to dz, and can be applied not only to dx=dz, but also to the situation where dx≠dz, and is more suitable for practical applications. It is used in production, and the number of grids required for the minimum wavelength is only 2.4 in this method, which is also better than other similar methods.

Description

A frequency domain forward modeling method and device

技术领域technical field

本发明涉及地质勘探技术领域，具体的讲是一种频率域正演方法及装置。The invention relates to the technical field of geological exploration, in particular to a frequency domain forward modeling method and device.

背景技术Background technique

寻找油气的主要手段有地球物理勘探、石油地质方法等，随着计算机技术的进步，地球物理勘探在油气寻找中越来越受到重视，更多的人力、物力投入其中。作为地球物理勘探的方法之一地震勘探方法在油气田勘探开发中具有极其重要作用，是其他地球物理方法无法比拟的。The main means of searching for oil and gas are geophysical exploration, petroleum geological methods, etc. With the advancement of computer technology, geophysical exploration has received more and more attention in oil and gas search, and more manpower and material resources have been invested in it. As one of the methods of geophysical exploration, seismic exploration method plays an extremely important role in the exploration and development of oil and gas fields, which is unmatched by other geophysical methods.

地震勘探方法每年都有大量新技术新方法涌现，其中全波形反演方法就是近年来的热点问题。这种方法可以利用地震波场的运动学和动力学信息重建地下速度结构，具有揭示复杂地质背景下构造与岩性细节信息的潜力。全波形反演可分为时间域、频率域以及拉普拉斯域。频率域方法以其如下优点备受青睐：首先频率域方法可以同时模拟多炮数据，其次不同于时间域方法，它可以只对某一频率的数据进行反演，第三，这种方法不存在累积误差。A large number of new technologies and methods emerge every year in seismic exploration methods, among which the full waveform inversion method is a hot topic in recent years. This method can use the kinematics and dynamics information of the seismic wave field to reconstruct the subsurface velocity structure, and has the potential to reveal structural and lithological details in complex geological backgrounds. Full waveform inversion can be divided into time domain, frequency domain and Laplace domain. The frequency domain method is favored for its following advantages: firstly, the frequency domain method can simulate multi-shot data at the same time; secondly, unlike the time domain method, it can only invert the data of a certain frequency; thirdly, this method does not exist cumulative error.

频率域正演是频率域全波形反演的基础，只有做好频率域正演，频率域反演方法才能得到有效利用。近几年来，也有很多学者提出各种正演差分方法。Frequency-domain forward modeling is the basis of frequency-domain full waveform inversion. Only when frequency-domain forward modeling is done well, frequency-domain inversion methods can be effectively utilized. In recent years, many scholars have proposed various forward difference methods.

频率域正演的关键是平衡正演精度和正演的计算量，较高的精度必然导致较高的计算量。最早1990年Pratt和Worthington在文献中提出了经典的五点差分格式，但是这种差分格式最小波长内所需要的网格点为13个(在相对误差1％以内，后面所说的最小波长内所需网格点都是基于这个标准)，导致计算量过大，并不能很好的用于频率域正演。为了减少最小波长内所需要的网格点数，Jo等人1996年发表文章提出一种新的九点旋转坐标系下的差分格式，从而有效降低了最小波长内所需要的网格点数，进而减少了计算机内存占有量和计算机运算量。以上两种方法都是基于旋转坐标系，导致了他们都有一个缺陷，就是不能适用于速度模型dx(水平方向的空间采样间隔)不等于dz(深度方向的空间采样间隔)的情况，而在实际勘探生产中，水平采样间隔经常不等于深度方向采样间隔，因此这种弊端必须克服。The key of forward modeling in frequency domain is to balance the forward modeling accuracy and forward modeling calculation amount, and higher precision will inevitably lead to higher calculation amount. As early as 1990, Pratt and Worthington proposed the classic five-point differential scheme in the literature, but the grid points required in the minimum wavelength of this differential scheme are 13 (within a relative error of 1%, the minimum wavelength mentioned later The required grid points are all based on this standard), resulting in an excessive amount of calculation, which cannot be well used for forward modeling in the frequency domain. In order to reduce the number of grid points required in the minimum wavelength, Jo et al. published an article in 1996 and proposed a new differential format under the nine-point rotating coordinate system, which effectively reduced the number of grid points required in the minimum wavelength, thereby reducing The amount of computer memory occupied and the amount of computer operations. The above two methods are based on the rotating coordinate system, which leads to a defect that they cannot be applied to the case where the velocity model dx (spatial sampling interval in the horizontal direction) is not equal to dz (spatial sampling interval in the depth direction), and in In actual exploration and production, the horizontal sampling interval is often not equal to the depth sampling interval, so this disadvantage must be overcome.

发明内容Contents of the invention

本发明的目的在于解决现有技术主要集中在二阶频率域正演，并不能满足高精度正演需求，部分高阶正演又存在适用条件的限制，The purpose of the present invention is to solve the problem that the existing technology mainly focuses on the second-order frequency domain forward modeling, which cannot meet the requirements of high-precision forward modeling, and some high-order forward modeling has limitations in applicable conditions.

本发明实施例提供了一种频率域正演方法，包括，An embodiment of the present invention provides a frequency domain forward modeling method, including:

建立如下的17点格式的差分公式，Establish the difference formula in the following 17-point format,

其中in

$\begin{matrix} {\overset{__}{P P}}_{m m + + 11,, n no} = = {α α}_{22} {P P}_{m m + + 11,, n no} + + \frac{11 - - {α α}_{22}}{22} (({P P}_{m m + + 11,, n no + + 11} + + {P P}_{m m + + 11,, n no + + 11})),, \\ {\overset{__}{P P}}_{m m - - 11,, n no} = = {α α}_{22} {P P}_{m m - - 11,, n no} + + \frac{11 - - {α α}_{11}}{22} (({P P}_{m m - - 11,, n no + + 11} + + {P P}_{m m - - 11,, n no - - 11})),, \\ {\overset{__}{P P}}_{m m + + 22,, n no} = = {α α}_{11} {P P}_{m m + + 22,, n no} + + \frac{11 - - {α α}_{11}}{22} (({P P}_{m m + + 22,, n no + + 22} + + {P P}_{m m + + 22,, n no - - 22})),, \\ {\overset{__}{P P}}_{m m - - 22,, n no} = = {α α}_{11} {P P}_{m m - - 22,, n no} + + \frac{11 - - {α α}_{11}}{22} (({P P}_{m m - - 22,, n no + + 22} + + {P P}_{m m - - 22,, n no - - 22})),, \\ {\overset{__}{P P}}_{m m,, n no} = = {α α}_{33} {P P}_{m m,, n no} + + {α α}_{44} (({P P}_{m m,, n no + + 11} + + {P P}_{m m,, n no - - 11})) + + \frac{11 - - 22 {α α}_{44} - - {α α}_{33}}{22} (({P P}_{m m,, n no + + 22} + + {P P}_{m m,, n no - - 22})),, \end{matrix} - - - - - - ((33))$

和and

$\begin{matrix} {\overset{~ ~}{P P}}_{m m,, n no + + 11} = = {β β}_{22} {P P}_{m m,, n no + + 11} + + \frac{11 - - {β β}_{22}}{22} (({P P}_{m m + + 11,, n no + + 11} + + {P P}_{m m - - 11,, n no + + 11})),, \\ {\overset{~ ~}{P P}}_{m m,, n no - - 11} = = {β β}_{22} {P P}_{m m,, n no + + 11} + + \frac{11 - - {β β}_{22}}{22} (({P P}_{m m + + 11,, n no - - 11} + + {P P}_{m m - - 11,, n no - - 11})),, \\ {\overset{~ ~}{P P}}_{m m,, n no + + 22} = = {β β}_{11} {P P}_{m m,, n no + + 22} + + \frac{11 - - {β β}_{11}}{22} (({P P}_{m m + + 22,, n no + + 22} + + {P P}_{m m - - 22,, n no + + 22})),, \\ {\overset{~ ~}{P P}}_{m m,, n no - - 22} = = {β β}_{11} {P P}_{m m,, n no - - 22} + + \frac{11 - - {β β}_{11}}{22} (({P P}_{m m + + 22,, n no - - 22} + + {P P}_{m m - - 22,, n no - - 22})),, \\ {\overset{~ ~}{P P}}_{m m,, n no} = = {β β}_{33} {P P}_{m m,, n no} + + {β β}_{44} (({P P}_{m m + + 11,, n no} + + {P P}_{m m - - 11,, n no})) + + \frac{11 - - 22 {β β}_{44} - - {β β}_{33}}{22} (({P P}_{m m + + 22,, n no} + + {P P}_{m m - - 22,, n no})),, \end{matrix} - - - - - - ((44))$

其中P_m,n≈P(mΔz,nΔx)，Δx,Δz代表x方向采样间隔和z方向采样间隔，m和n代表z方向的网格坐标和x方向的网格坐标，α_i,β_i,b,c,d,e,f代表求导时候的加权系数；Among them, P _{m, n} ≈ P(mΔz, nΔx), Δx, Δz represent the sampling interval in the x direction and the sampling interval in the z direction, m and n represent the grid coordinates in the z direction and the grid coordinates in the x direction, α _i , β _i ,b,c,d,e,f represent the weighting coefficients when deriving;

构造稀疏矩阵；Construct a sparse matrix;

读入子波参数和速度模型；Read in wavelet parameters and velocity models;

频率循环得到单频波场；Frequency cycle to obtain a single-frequency wave field;

对所有频率波场求反傅里叶变换；Take the inverse Fourier transform of the wavefield at all frequencies;

得到正演结果。Get the forward result.

本发明实施例还提供了一种频率域正演装置，包括，The embodiment of the present invention also provides a frequency domain forward modeling device, including:

差分模块，用于建立如下的17点格式的差分公式，The difference module is used to establish the difference formula in the following 17-point format,

其中in

和and

稀疏矩阵模块，用于构造稀疏矩阵；Sparse matrix module for constructing sparse matrices;

读入模块，用于读入子波参数和速度模型；Read-in module, used to read in wavelet parameters and velocity models;

单频波场模块，用于频率循环得到单频波场；Single-frequency wave field module, used for frequency cycle to obtain single-frequency wave field;

反傅里叶变换模块，用于对所有频率波场求反傅里叶变换；The inverse Fourier transform module is used to calculate the inverse Fourier transform of all frequency wave fields;

结果输出模块，得到正演结果。The result output module to get the forward modeling result.

通过上述实施例的方法和装置，四阶高精度正演克服了以往对dx必须等于dz的条件限制，不仅可应用于dx＝dz,也可以适用于dx≠dz的情况，更适合在实际生产中运用，而且最小波长所需要的网格数在该方法中仅需要2.4个，也优于其他类似方法。Through the method and device of the above embodiment, the fourth-order high-precision forward modeling overcomes the previous condition that dx must be equal to dz, and can be applied not only to dx=dz, but also to the situation where dx≠dz, and is more suitable for actual production. In this method, the number of grids required for the minimum wavelength is only 2.4, which is also better than other similar methods.

附图说明Description of drawings

结合以下附图阅读对实施例的详细描述，本发明的上述特征和优点，以及额外的特征和优点，将会更加清楚。The above-mentioned features and advantages, as well as additional features and advantages of the present invention, will become more apparent when reading the detailed description of the embodiments in conjunction with the following drawings.

图1所示为本发明实施例一种频率域正演方法的流程图；Fig. 1 shows the flow chart of a kind of frequency domain forward modeling method of the embodiment of the present invention;

图2所示为发明实施例一种频率域正演装置的结构图；Fig. 2 is a structural diagram of a frequency domain forward modeling device according to an embodiment of the invention;

图3所示为本发明实施例一种频率域正演方法的具体流程图；Fig. 3 shows the specific flowchart of a kind of frequency domain forward modeling method of the embodiment of the present invention;

图4a和图4b所示为本发明实施例系数分配图；Figure 4a and Figure 4b show the coefficient distribution diagram of the embodiment of the present invention;

表1为在Δx≥Δz的情况下，对于不同的Δx/Δz对应的优化系数表；Table 1 is a table of optimization coefficients corresponding to different Δx/Δz in the case of Δx≥Δz;

表2所示为在Δx<Δz的情况下，对于不同的Δx/Δz对应的优化系数表；Table 2 shows the table of optimization coefficients corresponding to different Δx/Δz in the case of Δx<Δz;

图5a为速度模型及观测系统的分布图；Figure 5a is the distribution diagram of velocity model and observation system;

图5b为正演的结果示意图；Figure 5b is a schematic diagram of the forward modeling results;

图6a为速度模型及其观测系统的布置图；Figure 6a is the layout of the velocity model and its observation system;

图6b为单频波场示意图；Figure 6b is a schematic diagram of a single-frequency wave field;

图6c为采用本发明实施例得到的正演波场在时间域的示意图；Fig. 6c is a schematic diagram of the forward wave field in the time domain obtained by using the embodiment of the present invention;

图6d所示的效果图。The effect diagram shown in Figure 6d.

具体实施方式Detailed ways

下面的描述可以使任何本领域技术人员利用本发明。具体实施例和应用中所提供的描述信息仅为示例。这里所描述的实施例的各种延伸和组合对于本领域的技术人员是显而易见的，在不脱离本发明的实质和范围的情况下，本发明定义的一般原则可以应用到其他实施例和应用中。因此，本发明不只限于所示的实施例，本发明涵盖与本文所示原理和特征相一致的最大范围。The following description will enable any person skilled in the art to utilize the present invention. The descriptions provided in the specific embodiments and applications are examples only. Various extensions and combinations of the embodiments described herein will be apparent to those skilled in the art, and the general principles defined herein may be applied to other embodiments and applications without departing from the spirit and scope of the invention . Thus, the present invention is not limited to the embodiments shown, but the present invention covers the widest scope consistent with the principles and features shown herein.

如图1所示为本发明实施例一种频率域正演方法的流程图。FIG. 1 is a flowchart of a frequency domain forward modeling method according to an embodiment of the present invention.

包括步骤101，建立如下的17点格式的差分公式，Including step 101, establishing the difference formula of the following 17-point format,

其中in

和and

其中P_m,n≈P(mΔz,nΔx)，Δx,Δz代表x方向采样间隔和z方向采样间隔，m和n代表z方向的网格坐标和x方向的网格坐标，α_i,β_i,b,c,d,e,f代表求导时候的加权系数。Among them, P _{m, n} ≈ P(mΔz, nΔx), Δx, Δz represent the sampling interval in the x direction and the sampling interval in the z direction, m and n represent the grid coordinates in the z direction and the grid coordinates in the x direction, α _i , β _i ,b,c,d,e,f represent the weighting coefficients when deriving.

步骤102，构造稀疏矩阵。Step 102, construct a sparse matrix.

步骤103，读入子波参数和速度模型。Step 103, read in wavelet parameters and velocity model.

步骤104，频率循环得到单频波场。Step 104, frequency cycling to obtain a single-frequency wave field.

步骤105，对所有频率波场求反傅里叶变换。Step 105, inverse Fourier transform of all frequency wave fields.

步骤106，得到正演结果。In step 106, the forward modeling result is obtained.

频率域正演是全波形反演的基础，是全波形反演中的一个重要环节，而且不可缺少，全波形反演是当今地震勘探的热点问题，它可以利用地震波场的运动学和动力学信息重建地下速度结构，具有揭示复杂地质背景下构造与岩性细节信息的潜力，从而更详细的了解地层的构造与分布，为石油勘探预测油气藏。Frequency domain forward modeling is the basis of full waveform inversion, an important link in full waveform inversion, and it is indispensable. Full waveform inversion is a hot issue in seismic exploration today. It can use the kinematics and dynamics of seismic wave field Information reconstruction of underground velocity structure has the potential to reveal structural and lithological details in complex geological backgrounds, so as to understand the structure and distribution of strata in more detail, and predict oil and gas reservoirs for oil exploration.

上述步骤102-步骤106均为现有技术中频率域正演的常规方法，在本申请文中不再赘述。The above steps 102 to 106 are conventional methods of frequency domain forward modeling in the prior art, and will not be described in detail in this application.

上述步骤101的17点格式差分公式可以针对dx≠dz的情况，还可以针对dx＝dz的情况。The 17-point format difference formula in the above step 101 may be for the case of dx≠dz, and may also be for the case of dx=dz.

在所述步骤101之后还包括，求取所述加权系数的最优值。After the step 101, it also includes calculating the optimal value of the weighting coefficient.

所述求取加权系数的最优值进一步包括，由波动方程(1)出发，把平面波代入公式(1)，整理可得频散方程，其中，P₀为波初始的振幅值、i为虚数单位，The optimal value of described seeking weighting coefficient further comprises, by wave equation (1) Starting, take the plane wave Substituting into formula (1), the dispersion equation can be obtained, where P ₀ is the initial amplitude value of the wave, i is the imaginary number unit,

其中k代表的是波数，where k represents the wave number,

当dx≥dz时，定义r＝Δx/Δz，因此上述公式(4)可以改写成如下形式When dx≥dz, define r=Δx/Δz, so the above formula (4) can be rewritten into the following form

在dx≥dz的时候，最小波长内所需要的网格点G的定义为G＝λ/Δx，其中λ为波长，When dx≥dz, the grid point G required in the minimum wavelength is defined as G=λ/Δx, where λ is the wavelength,

至此，归一化的相速度公式为So far, the normalized phase velocity formula is

$\begin{matrix} \frac{{V V}_{ph pH}}{v v} = = \frac{G G}{22 π π} {[[\frac{E E. + + {r r}^{22} F f}{b b + + 22 c c ((A A + + B B)) + + 44 dAB DAB + + 22 e e ((C C + + D D.)) + + 44 fCD f}]]}^{11 / / 22},, \\ A A = = cos cos ((\frac{22 π π cos cos ((θ θ))}{G G})),, \\ B B = = cos cos ((\frac{22 π π sin sin ((θ θ))}{rG rG})),, \\ C C = = cos cos ((\frac{44 π π cos cos ((θ θ))}{G G})),, \\ D D. = = cos cos ((\frac{44 π π sin sin ((θ θ))}{rG rG})) \\ E E. = = \frac{11}{66} {α α}_{11} ((C C - - CD cd)) - - \frac{88}{33} {α α}_{22} ((A A - - AB AB)) + + \frac{55}{22} {α α}_{33} ((11 - - D D.)) + + 55 {α α}_{44} ((B B - - D D.)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} D D.,, \\ F f = = \frac{11}{66} {β β}_{11} ((D D. - - CD cd)) - - \frac{88}{33} {β β}_{22} ((B B - - AB AB)) + + \frac{55}{22} {β β}_{22} ((11 - - C C)) + + 55 {β β}_{44} ((A A - - C C)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} C C,, \end{matrix}- - - - - - - ((77))$

其中k_x＝kcos(θ),k_z＝ksin(θ)，in k _x = kcos(θ), k _z = ksin(θ),

求得上述最优的α_i,β_i,b,c,d,e,f，使得V_ph与v之间差值最小；Find the above optimal α _i , β _i , b, c, d, e, f, so that the difference between V _ph and v is the smallest;

当dx<dz时，定义r＝Δz/Δx，G的定义为G＝λ/Δz；When dx<dz, define r=Δz/Δx, G is defined as G=λ/Δz;

上述公式(4)可以改写为The above formula (4) can be rewritten as

则归一化的相速度公式为Then the normalized phase velocity formula is

$\begin{matrix} \frac{{V V}_{ph pH}}{v v} = = \frac{G G}{22 π π} {[[\frac{E E. + + {r r}^{22} F f}{b b + + 22 c c ((A A + + B B)) + + 44 dAB DAB + + 22 e e ((C C + + D D.)) + + 44 fCD f}]]}^{11 / / 22},, \\ A A = = cos cos ((\frac{22 π π cos cos ((θ θ))}{rG rG})),, \\ B B = = cos cos ((\frac{22 π π sin sin ((θ θ))}{G G})),, \\ C C = = cos cos ((\frac{44 π π cos cos ((θ θ))}{rG rG})),, \\ D D. = = cos cos ((\frac{44 π π sin sin ((θ θ))}{G G})) \\ E E. = = \frac{11}{66} {α α}_{11} ((C C - - CD cd)) - - \frac{88}{33} {α α}_{22} ((A A - - AB AB)) + + \frac{55}{22} {α α}_{33} ((11 - - D D.)) + + 55 {α α}_{44} ((B B - - D D.)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} D D.,, \\ F f = = \frac{11}{66} {β β}_{11} ((D D. - - CD cd)) - - \frac{88}{33} {β β}_{22} ((B B - - AB AB)) + + \frac{55}{22} {β β}_{22} ((11 - - C C)) + + 55 {β β}_{44} ((A A - - C C)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} C C,, \end{matrix}- - - - - - - ((77))$

求得上述最优的α_i,β_i,b,c,d,e,f，使得V_ph与v之间差值最小。Obtain the above optimal α _i , β _i , b, c, d, e, f to minimize the difference between V _ph and v.

在上述求最优α_i,β_i,b,c,d,e,f的步骤中，可以采用最小二乘法或者还可以使用高斯牛顿法求得最优α_i,β_i,b,c,d,e,f。In the above step of finding the optimal α _i , β _i ,b,c,d,e,f, the least square method or the Gauss-Newton method can be used to obtain the optimal α _i ,β _i ,b,c, d, e, f.

作为本发明的一个实施例，在步骤103读入子波参数和速度模型之后还包括，在所述速度模型的边界加载最佳匹配层PML边界条件，用于吸收边界反射的波场，防止边界反射的波场回折到地表。As an embodiment of the present invention, after the wavelet parameters and the velocity model are read in step 103, it also includes loading the best matching layer PML boundary condition on the boundary of the velocity model, which is used to absorb the wave field reflected by the boundary and prevent the boundary condition from The reflected wavefield bounces back to the surface.

作为本发明的一个实施例，所述步骤102中进一步包括，采用如下方式构造稀疏矩阵：As an embodiment of the present invention, step 102 further includes constructing a sparse matrix in the following manner:

c₁P_m+1,n+c₂P_m+1,n+1+c₂P_m+1,n-1+c₁P_m-1,n+c₂P_m-1,n+1+c₂P_m-1,n-1+c₃P_m+2,n+c₄P_m+2,n+2+ (8)c₄P_m+2,n-2+c₃P_m-2,n+c₄P_m-2,n+2+c₄P_m-2,n-2+c₅P_m,n+c₆P_m,n+1+c₆P_m,n-1+c₇P_m,n+2+c₇P_m,n-2＝S,c ₁ P _m+1,n +c ₂ P _m+1,n+1 +c ₂ P _m+1,n-1 +c ₁ P _m-1,n +c ₂ P _m-1,n+1 +c ₂ P _m-1,n-1 +c ₃ P _m+2,n +c ₄ P _m+2,n+2 + (8)c ₄ P _m+2,n-2 +c ₃ P _{m -2,n} +c ₄ P _m-2,n+2 +c ₄ P _m-2,n-2 +c ₅ P _m,n +c ₆ P _m,n+1 +c ₆ P _{m,n- 1} +c ₇ P _m,n+2 +c ₇ P _m,n-2 ＝S,

公式中的系数如下The coefficients in the formula are as follows

每一个(m，n)对应一个方程，共有m*n个方程，方程右边是震源项，是一个行数为m*n，列数为1列的一个矩阵，左边是一个行数m*n，列数也是m*n的一个矩阵，求解该方程组就得到波场每个位置的单频波场，即P_m,n。Each (m, n) corresponds to an equation, a total of m*n equations, the right side of the equation is the source item, which is a matrix with m*n rows and 1 column, and the left side is a row number m*n , and the number of columns is also a matrix of m*n. Solving this system of equations can obtain the single-frequency wave field at each position of the wave field, that is, P _m,n .

通过上述的方法，四阶高精度正演克服了以往对dx必须等于dz的条件限制，不仅可应用于dx＝dz,也可以适用于dx≠dz的情况，更适合在实际生产中运用，而且最小波长所需要的网格数在该方法中仅需要2.4个，也优于其他类似方法。Through the above method, the fourth-order high-precision forward modeling overcomes the previous condition that dx must be equal to dz, and can be applied not only to dx=dz, but also to the situation where dx≠dz, and is more suitable for use in actual production. The number of grids required for the minimum wavelength is only 2.4 in this method, which is also better than other similar methods.

如图2所示为发明实施例一种频率域正演装置的结构图。FIG. 2 is a structural diagram of a frequency domain forward modeling device according to an embodiment of the invention.

包括差分模块201，用于建立如下的17点格式的差分公式，Including difference module 201, for establishing the following difference formula in 17-point format,

其中in

和and

稀疏矩阵模块202，用于构造稀疏矩阵。The sparse matrix module 202 is configured to construct a sparse matrix.

读入模块203，用于读入子波参数和速度模型。The read-in module 203 is used to read in wavelet parameters and velocity models.

单频波场模块204，用于频率循环得到单频波场。The single-frequency wave field module 204 is used for frequency cycling to obtain a single-frequency wave field.

反傅里叶变换模块205，用于对所有频率波场求反傅里叶变换。The inverse Fourier transform module 205 is used to calculate the inverse Fourier transform of all frequency wave fields.

结果输出模块206，得到正演结果。The result output module 206, to obtain the forward modeling result.

上述稀疏矩阵模块202-结果输出模块206均可以采用现有技术中频率域正演的常规方法，在本申请文中不再赘述。The aforementioned sparse matrix module 202-result output module 206 can all adopt conventional methods of frequency-domain forward modeling in the prior art, which will not be repeated in this application.

上述差分模块201得到的17点格式差分公式可以针对dx≠dz的情况，还可以针对dx＝dz的情况。The 17-point format difference formula obtained by the above-mentioned difference module 201 may be for the case of dx≠dz, and may also be for the case of dx=dz.

还包括加权系数模块207，连接于所述差分模块201和稀疏矩阵模块202之间，用于求取所述加权系数的最优值。It also includes a weighting coefficient module 207, connected between the difference module 201 and the sparse matrix module 202, and used to calculate the optimal value of the weighting coefficient.

所述加权系数模块207求取加权系数的最优值进一步包括，由波动方程(1)出发，把平面波代入公式(1)，整理可得频散方程，其中，P0为波初始的振幅值、i为虚数单位，The weighting coefficient module 207 to obtain the optimal value of the weighting coefficient further includes, by the wave equation (1) Starting, take the plane wave Substituting into formula (1), the dispersion equation can be obtained, where P0 is the initial amplitude value of the wave, i is the imaginary number unit,

其中k代表的是波数，where k represents the wave number,

其中k_x＝kcos(θ),k_z＝ksin(θ)，in k _x = kcos(θ), k _z = ksin(θ),

上述公式(4)可以改写为The above formula (4) can be rewritten as

则归一化的相速度公式为Then the normalized phase velocity formula is

所述加权系数模块207在上述求最优α_i,β_i,b,c,d,e,f的过程中，可以采用最小二乘法或者还可以使用高斯牛顿法求得最优α_i,β_i,b,c,d,e,f。In the process of obtaining the optimal α _i , β _i , b, c, d, e, f, the weighting coefficient module 207 can use the least square method or the Gauss-Newton method to obtain the optimal α _i , β _i ,b,c,d,e,f.

作为本发明的一个实施例，还包括吸收边界模块208，连接于所述读入模块203和单频波场模块204之间，用于在所述速度模型的边界加载最佳匹配层PML边界条件，用于吸收边界反射的波场，防止边界反射的波场回折到地表。As an embodiment of the present invention, it also includes an absorption boundary module 208, which is connected between the read-in module 203 and the single-frequency wave field module 204, and is used to load the best matching layer PML boundary condition on the boundary of the velocity model , which is used to absorb the wave field reflected by the boundary and prevent the wave field reflected by the boundary from returning to the surface.

作为本发明的一个实施例，所述稀疏矩阵模块202还进一步用于，采用如下方式构造稀疏矩阵：As an embodiment of the present invention, the sparse matrix module 202 is further configured to construct a sparse matrix in the following manner:

公式中的系数如下The coefficients in the formula are as follows

通过上述的装置，四阶高精度正演克服了以往对dx必须等于dz的条件限制，不仅可应用于dx＝dz,也可以适用于dx≠dz的情况，更适合在实际生产中运用，而且最小波长所需要的网格数在该方法中仅需要2.4个，也优于其他类似方法。Through the above-mentioned device, the fourth-order high-precision forward modeling overcomes the previous condition that dx must be equal to dz, and can be applied not only to dx=dz, but also to the situation where dx≠dz, and is more suitable for use in actual production. The number of grids required for the minimum wavelength is only 2.4 in this method, which is also better than other similar methods.

如图3所示为本发明实施例一种频率域正演方法的具体流程图。FIG. 3 is a specific flow chart of a frequency domain forward modeling method according to an embodiment of the present invention.

包括步骤301，建立如下的17点格式的差分公式。Step 301 is included to establish the difference formula in the following 17-point format.

频率域的波动方程可以表示成The wave equation in the frequency domain can be expressed as

P代表波场，代表圆频率，v是波在介质中的传播速度。P stands for wave field, Represents the circular frequency and v is the propagation speed of the wave in the medium.

利用有限差分方法和平均导数的方法，推导出新的17点格式的差分公式，如下：Using the finite difference method and the method of the average derivative, a new difference formula of the 17-point format is derived, as follows:

其中in

和and

其中P_m,n≈P(mΔz,nΔx)，Δx,Δz代表x方向采样间隔和z方向采样间隔，与文中的dx与dz含义一样，m和n代表z方向的网格坐标和x方向的网格坐标。α_i,β_i,b,c,d,e,f代表求导时候的加权系数，具体的系数分配如图4a和图4b所示。Among them, P _{m, n} ≈ P(mΔz, nΔx), Δx, Δz represent the sampling interval in the x direction and the sampling interval in the z direction, which is the same as the meaning of dx and dz in the text, m and n represent the grid coordinates in the z direction and the grid coordinates. α _i , β _i , b, c, d, e, and f represent the weighting coefficients when deriving, and the specific coefficient distribution is shown in Figure 4a and Figure 4b.

这种新的格式可以适用于dx≠dz的情况，增加了该差分公式的适用性。在推导dx≥dz中包含了dx＝dz的情形，而且在表格1中也求取出了dx＝dz情况下的优化参数，有了优化参数就有了dx＝dz情形下的频率域波动方程差分公式，有了这个差分公式就可以做dx＝dz情形下的频率域正演。后面的实施例中就是用本发明做的dx＝dz情况的频率域正演。This new format can be applied to the case of dx≠dz, increasing the applicability of the difference formula. The case of dx=dz is included in the derivation of dx≥dz, and the optimization parameters in the case of dx=dz are also obtained in Table 1. With the optimization parameters, there is the difference of the frequency domain wave equation in the case of dx=dz Formula, with this differential formula, the forward modeling in the frequency domain under the condition of dx=dz can be done. In the following embodiments, the forward modeling in the frequency domain of the case of dx=dz is done by using the present invention.

步骤302，求取优化的加权系数。Step 302, obtain optimized weighting coefficients.

基于这种差分格式，进行频率域声波正演，还需要解决频散的问题，小的频散可以获得较好的正演效果。因此进行频散公式的推导：Based on this differential format, it is necessary to solve the problem of dispersion for sound wave forward modeling in the frequency domain. Small dispersion can obtain better forward modeling effect. Therefore, the derivation of the dispersion formula is carried out:

首先由波动方程公式(1)出发，把平面波代入公式(1)，整理可得频散方程，其中，P₀为波初始的振幅值、i为虚数单位Firstly, starting from the wave equation formula (1), the plane wave Substituting into formula (1), the dispersion equation can be obtained, where P ₀ is the initial amplitude value of the wave, and i is the imaginary number unit

其中k代表的是波数，其他参数含义已在上文中指出，不再赘述。Among them, k represents the wave number, and the meanings of other parameters have been pointed out above, so I won’t repeat them here.

分两种情况：There are two situations:

(1)当dx≥dz时，定义r＝Δx/Δz，因此上述公式(4)可以改写成如下形式(1) When dx≥dz, define r=Δx/Δz, so the above formula (4) can be rewritten into the following form

在dx≥dz的时候，最小波长内所需要的网格点G的定义为G＝λ/Δx(λ为波长)，When dx≥dz, the grid point G required in the minimum wavelength is defined as G=λ/Δx (λ is the wavelength),

其中k_x＝kcos(θ),k_z＝ksin(θ)in k _x ＝kcos(θ),k _z ＝ksin(θ)

减少频散就是寻找最优系数α_i,β_i,b,c,d,e,f，使得V_ph与v之间差值最小，可以利用诸如最小二乘法求得最优值，也可以通过Min等人在2000年文献中提出的高斯牛顿方法。从而获得公式(2)(3)(4)中不同dx和dz对应的α_i,β_i,b,c,d,e,f，具体参见表1(见后附列表)为在Δx≥Δz的情况下，对于不同的Δx/Δz对应的优化系数表。To reduce dispersion is to find the optimal coefficients α _i , β _i , b, c, d, e, f, so that the difference between V _ph and v is the smallest, and the optimal value can be obtained by using the least square method, or by The Gauss-Newton method proposed by Min et al. in the 2000 paper. Thus, α _i , β _i , b, c, d, e, f corresponding to different dx and dz in formula (2)(3)(4) can be obtained. For details, see Table 1 (see attached list) for Δx≥Δz In the case of different Δx/Δz corresponding optimization coefficient table.

当dx<dz时，r＝Δz/Δx，G(最小波长内所需要的网格点)的定义为G＝λ/Δz(λ为波长)，可以得到dx<dz时候的不同dx和dz比值下的优化系数，如表2(见后附列表)所示为在Δx<Δz的情况下，对于不同的Δx/Δz对应的优化系数表。When dx<dz, r=Δz/Δx, G (the grid point required in the minimum wavelength) is defined as G=λ/Δz (λ is the wavelength), and different ratios of dx and dz can be obtained when dx<dz The optimization coefficient below, as shown in Table 2 (see attached list), is the table of optimization coefficients corresponding to different Δx/Δz under the condition of Δx<Δz.

步骤303，构造稀疏矩阵。Step 303, construct a sparse matrix.

构造稀疏矩阵，该步骤可用公式解释如下To construct a sparse matrix, this step can be explained by the formula as follows

公式中的系数如下The coefficients in the formula are as follows

上述公式(8)实质上代表一个方程组，每一个(m，n)对应一个方程，因此共有m*n个方程。在给定速度模型和观测系统的情况下上述方程系数是确定的，所求的即是P_m,n，方程右边是震源项，是一个行数为m*n，列数为1列的一个矩阵，左边是一个行数m*n，列数也是m*n的一个矩阵，求解该方程组就得到波场每个位置的单频波场，即P_m，n，在此要用到LU分解来解该方程组。The above formula (8) essentially represents a system of equations, and each (m, n) corresponds to an equation, so there are m*n equations in total. In the case of a given velocity model and observation system, the coefficients of the above equation are determined, and what is sought is P _{m,n .} Matrix, on the left is a matrix with m*n rows and m*n columns. Solving the equations will get the single-frequency wave field at each position of the wave field, that is, P _{m, n} . LU is used here factorize to solve the system of equations.

其中，LU分解是将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积，主要应用在数值分析中，用来解线性方程、求反矩阵或计算行列式，具体公式和程序可以参考现有技术，在本申请中不再赘述。Among them, LU decomposition is to decompose a matrix into the product of a lower triangular matrix and an upper triangular matrix. It is mainly used in numerical analysis to solve linear equations, invert matrices or calculate determinants. For specific formulas and procedures, please refer to There are technologies, which will not be described in detail in this application.

步骤304，读入子波参数和速度模型。子波参数是时间域的，通过傅里叶变化使其变换到频率域，得到每个频率对应的幅值，对应上面公式(8)中的S(1列的一个矩阵)。Step 304, read in wavelet parameters and velocity model. The wavelet parameters are in the time domain, and they are transformed into the frequency domain by Fourier transformation to obtain the corresponding amplitude of each frequency, which corresponds to S (a matrix of 1 column) in the above formula (8).

步骤305，加载PML(最佳匹配层)边界条件，即在边界加一层吸收波场的介质，防止边界反射的波场回折到地表。Step 305 , load PML (best matching layer) boundary conditions, that is, add a layer of wavefield-absorbing medium on the boundary to prevent the wavefield reflected by the boundary from returning to the surface.

步骤306，对频率由低到高进行循环，每次循环得到的单频波场都存储在一个矩阵中。In step 306, the frequency is cycled from low to high, and the single-frequency wavefield obtained in each cycle is stored in a matrix.

步骤307，对波场内每个点进行所有频率的反傅里叶变换就换算成时间域的波场。In step 307, the inverse Fourier transform of all frequencies is performed on each point in the wave field to convert it into a wave field in the time domain.

步骤308，通过提取地表的波场记录得到正演记录。In step 308, the forward modeling record is obtained by extracting the surface wave field record.

为了清晰的看到上述方法的得到正演记录中的同相轴位置和频散情况，先用一个简单的层状模型，层状模型速度结构简单，易于分析，该模型的网格间距dx＝dz＝5m，水平方向和深度方向的网格点数均为100(nx＝nz＝100)，对于这样的模型，公式(8)和(9)中系数α_i,β_i,b,c,d,e,f采用表格1中r＝1的参数。震源放在(250m，20m)的位置，子波采用雷克子波，主频为20Hz，检波器位于Z＝20m的位置，时间采样间隔为2ms，接收记录时间是1s。图5a为速度模型及观测系统的分布图，图5b为正演的结果示意图，其中，图5a代表这一块区域地下构造的地震波传播速度，不同颜色代表不同的速度，如右侧图例所示，颜色趋于红色代表速度值高，最高是2000m/s，蓝色代表速度低，最低为1000m/s；白色的倒三角代表检波器的分布，检波器中间红色的点代表爆炸点所在位置，通过图5a和图5b表明同相轴位置正确，而且边界吸收较好，不存在干扰波，由于采用了四阶的差分格式，因此频散很弱，较好的压制了频散。In order to clearly see the event position and dispersion in the forward modeling records obtained by the above method, a simple layered model is used first. The layered model has a simple velocity structure and is easy to analyze. The grid spacing of the model is dx=dz =5m, the number of grid points in the horizontal and depth directions is 100 (nx=nz=100), for such a model, the coefficients α _i , β _i ,b,c,d, e, f use the parameters of r=1 in Table 1. The seismic source is placed at the position of (250m, 20m), the wavelet adopts the Leike wavelet, the main frequency is 20Hz, the geophone is located at the position of Z=20m, the time sampling interval is 2ms, and the receiving and recording time is 1s. Figure 5a is the distribution diagram of the velocity model and observation system, and Figure 5b is the schematic diagram of the forward modeling results, in which Figure 5a represents the seismic wave propagation velocity of the underground structure in this area, and different colors represent different velocities, as shown in the legend on the right, The color tends to red to represent the high speed value, the highest is 2000m/s, blue represents the low speed, the lowest is 1000m/s; the white inverted triangle represents the distribution of the geophone, and the red point in the middle of the geophone represents the location of the explosion point. Figure 5a and Figure 5b show that the position of the event axis is correct, and the boundary absorption is good, and there is no interference wave. Because the fourth-order difference scheme is used, the dispersion is very weak, and the dispersion is suppressed better.

为进一步验证算法的有效性，采用地震勘探行业认知度较高的Marmousi速度模型进行正演，该模型dx＝12.5m，dz＝4m，因此dx≠dz，我们采用表格1中r＝dx/dz＝3.125的系数参数。该模型我们取其中一部分正演，网格数为nx＝nz＝301，子波同样采用雷克子波，主频为12Hz，放置在x＝500m，z＝16m的地方，检波器位于z＝16m的地方。图6a为速度模型及其观测系统的布置图，该图中的颜色如上述图5a中所述，图6b为单频波场示意图，其中较红的部分表示振幅值较高，偏黄色部分表示振幅值较弱，越红说明振幅值越高，越黄说明振幅值越弱，图6c为采用本发明实施例得到的正演波场在时间域的示意图，为了验证其正确性，用时间域正演的方法，对该模型进行正演，得到了图6d所示的效果图，两图同相轴位置一致，说明本专利方法是正确的。In order to further verify the effectiveness of the algorithm, the Marmousi velocity model, which is well known in the seismic exploration industry, is used for forward modeling. The model dx=12.5m, dz=4m, so dx≠dz, we use r=dx/ Coefficient parameter of dz=3.125. We take a part of the model for forward modeling, the number of grids is nx=nz=301, the wavelet also adopts the Reker wavelet, the main frequency is 12Hz, it is placed at x=500m, z=16m, and the detector is located at z=16m The place. Figure 6a is the layout of the velocity model and its observation system. The colors in this figure are as described in Figure 5a above, and Figure 6b is a schematic diagram of the single-frequency wave field, where the redder part indicates higher amplitude values, and the yellowish part indicates The amplitude value is weaker, the redder it is, the higher the amplitude value is, and the yellower it is, the weaker the amplitude value is. Figure 6c is a schematic diagram of the forward wave field in the time domain obtained by using the embodiment of the present invention. In the forward modeling method, the forward modeling is carried out on the model, and the effect diagram shown in Fig. 6d is obtained. The positions of the event axes in the two diagrams are consistent, which shows that the method of this patent is correct.

通过上述的方法及装置，四阶高精度正演克服了以往对dx必须等于dz的条件限制，不仅可应用于dx＝dz,也可以适用于dx≠dz的情况，更适合在实际生产中运用，而且最小波长所需要的网格数在该方法中仅需要2.4个，也优于其他类似方法。Through the above method and device, the fourth-order high-precision forward modeling overcomes the previous condition that dx must be equal to dz, and can be applied not only to dx=dz, but also to the situation where dx≠dz, and is more suitable for use in actual production , and the number of grids required for the minimum wavelength is only 2.4 in this method, which is also better than other similar methods.

本发明可以以任何适当的形式实现，包括硬件、软件、固件或它们的任意组合。本发明可以根据情况有选择的部分实现，比如计算机软件执行于一个或多个数据处理器以及数字信号处理器。本文的每个实施例的元素和组件可以在物理上、功能上、逻辑上以任何适当的方式实现。事实上，一个功能可以在独立单元中、在一组单元中、或作为其他功能单元的一部分来实现。因此，该系统和方法既可以在独立单元中实现，也可以在物理上和功能上分布于不同的单元和处理器之间。The invention can be implemented in any suitable form including hardware, software, firmware or any combination of these. The invention can be implemented in select parts, as appropriate, as computer software executing in one or more data processors and digital signal processors. The elements and components of each embodiment herein may be physically, functionally and logically implemented in any suitable way. In fact a function may be implemented in a separate unit, in a group of units or as part of other functional units. Thus, the systems and methods can be implemented in separate units or can be physically and functionally distributed between different units and processors.

在相关领域中的技术人员将会认识到，本发明的实施例有许多可能的修改和组合，虽然形式略有不同，仍采用相同的基本机制和方法。为了解释的目的，前述描述参考了几个特定的实施例。然而，上述的说明性讨论不旨在穷举或限制本文所发明的精确形式。前文所示，许多修改和变化是可能的。所选和所描述的实施例，用以解释本发明的原理及其实际应用，用以使本领域技术人员能够最好地利用本发明和各个实施例的针对特定应用的修改、变形。Those skilled in the relevant art will recognize that there are many possible modifications and combinations of the embodiments of the invention, albeit in slightly different forms, still employing the same basic mechanisms and methods. The foregoing description, for purposes of explanation, has referred to a few specific embodiments. However, the illustrative discussions above are not intended to be exhaustive or to limit the precise forms of the inventions herein. As indicated above, many modifications and variations are possible. The embodiments were chosen and described to explain the principles of the invention and its practical application, and to enable those skilled in the art to best utilize the invention and various modifications and variations of the embodiments for specific applications.

Claims

1. A frequency domain forward modeling method is characterized in that comprising,

Establish the difference formula in the following 17-point format,

in

{\overset{&OverBar; &OverBar;}{P P}}_{m m + + 11,, n no} = = {α α}_{22} {P P}_{m m + + 11,, n no} + + \frac{11 - - {α α}_{22}}{22} (({P P}_{m m + + 11,, n no - - 11})),,

{\overset{&OverBar; &OverBar;}{P P}}_{m m - - 11,, n no} = = {α α}_{22} {P P}_{m m - - 11,, n no} + + \frac{11 - - {α α}_{22}}{22} (({P P}_{m m - - 11,, n no + + 11} + + {P P}_{m m - - 11,, n no - - 11})) - - - - - - ((33))

{\overset{&OverBar; &OverBar;}{P P}}_{m m + + 22,, n no} = = {α α}_{11} {P P}_{m m + + 22,, n no} + + \frac{11 - - {α α}_{11}}{22} (({P P}_{m m + + 22,, n no + + 22} + + {P P}_{m m + + 22,, n no - - 22})),,

{\overset{&OverBar; &OverBar;}{P P}}_{m m - - 22,, n no} = = {α α}_{11} {P P}_{m m - - 22,, n no} + + \frac{11 - - {α α}_{11}}{22} (({P P}_{m m - - 22,, n no + + 22} + + {P P}_{m m - - 22,, n no - - 22})),,

{\overset{&OverBar; &OverBar;}{P P}}_{m m,, n no} = = {α α}_{33} {P P}_{m m,, n no} + + {α α}_{44} (({P P}_{m m,, n no + + 11} + + {P P}_{m m,, n no - - 11})) + + \frac{11 - - {22 α α}_{44} - - {α α}_{33}}{22} (({P P}_{m m,, n no + + 22} + + {P P}_{m m,, n no - - 22})),,

and

{\overset{~ ~}{P P}}_{m m,, n no + + 11} = = {β β}_{22} {P P}_{m m,, n no + + 11} + + \frac{11 - - {β β}_{22}}{22} (({P P}_{m m + + 11,, n no + + 11} + + {P P}_{m m - - 11,, n no + + 11})),,

{\overset{~ ~}{P P}}_{m m,, n no - - 11} = = {β β}_{22} {P P}_{m m,, n no - - 11} + + \frac{11 - - {β β}_{22}}{22} (({P P}_{m m + + 11,, n no - - 11} + + {P P}_{m m - - 11,, n no - - 11})),, - - - - - - ((44))

{\overset{~ ~}{P P}}_{m m,, n no + + 22} = = {β β}_{11} {P P}_{m m,, n no + + 22} + + \frac{11 - - {β β}_{11}}{22} (({P P}_{m m + + 22,, n no + + 22} + + {P P}_{m m - - 22,, n no + + 22})),,

{\overset{~ ~}{P P}}_{m m,, n no - - 22} = = {β β}_{11} {P P}_{m m,, n no - - 22} + + \frac{11 - - {β β}_{11}}{22} (({P P}_{m m + + 22,, n no - - 22} + + {P P}_{m m - - 22,, n no - - 22})),,

{\overset{~ ~}{P P}}_{m m,, n no} = = {β β}_{33} {P P}_{m m,, n no} + + {β β}_{44} (({P P}_{m m + + 11,, n no} + + {P P}_{m m - - 11,, n no})) + + \frac{11 - - {22 β β}_{44} - - {β β}_{33}}{22} (({P P}_{m m + + 22,, n no} + + {P P}_{m m - - 22,, n no})),,

Among them, P _{m, n} ≈ P(mΔz, nΔx), Δx, Δz represent the sampling interval in the x direction and the sampling interval in the z direction, m and n represent the grid coordinates in the z direction and the grid coordinates in the x direction, α _i , β _i ,b,c,d,e,f represent the weighting coefficients when deriving;

Construct a sparse matrix;

Read in wavelet parameters and velocity models;

Frequency cycle to obtain a single-frequency wave field;

Take the inverse Fourier transform of the wavefield at all frequencies;

Get the forward result.

2. A frequency-domain forward modeling method according to claim 1, characterized in that, after establishing the difference formula in the 17-point format, it also includes calculating the optimal value of the weighting coefficient.

3. A kind of frequency domain forward modeling method according to claim 2, it is characterized in that, the optimal value of described seeking weighting coefficient further comprises, by wave equation (1) Starting, take the plane wave Substituting into formula (1), the dispersion equation can be obtained, where P ₀ is the initial amplitude value of the wave, i is the imaginary number unit,

A=cos(k _x Δx),

B＝cos(k _z Δz), (5)

C＝cos(2k _x Δx),

D＝cos(2k _z Δz),

E E. = = \frac{11}{66} {α α}_{11} ((C C - - CD cd)) - - \frac{88}{33} {α α}_{22} ((A A - - AB AB)) + + \frac{55}{22} {α α}_{33} ((11 - - D D.)) + + {55 α α}_{44} ((B B - - D D.)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} D D.,,

F f = = \frac{11}{66} {β β}_{11} ((D D. - - CD cd)) - - \frac{88}{33} {β β}_{22} ((B B - - AB AB)) + + \frac{55}{22} {β β}_{33} ((11 - - C C)) + + {55 β β}_{44} ((A A - - C C)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} C C,,

where k represents the wave number,

When dx≥dz, define r=Δx/Δz, so the above formula (4) is rewritten into the following form

A=cos(k _x Δx),

B=cos(k _z Δz),

C＝cos(2k _x Δx), (6)

D＝cos(2k _z Δz),

E E. = = \frac{11}{66} {α α}_{11} ((C C - - CD cd)) - - \frac{88}{33} {α α}_{22} ((A A - - AB AB)) + + \frac{55}{22} {α α}_{33} ((11 - - D D.)) + + {55 α α}_{44} ((B B - - D D.)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} D D.,,

F f = = \frac{11}{66} {β β}_{11} ((D D. - - CD cd)) - - \frac{88}{33} {β β}_{22} ((B B - - AB AB)) + + \frac{55}{22} {β β}_{33} ((11 - - C C)) + + {55 β β}_{44} ((A A - - C C)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} C C,,

r=Δx/Δz,

When dx≥dz, the grid point G required in the minimum wavelength is defined as G=λ/Δx, where λ is the wavelength,

So far, the normalized phase velocity formula is

\frac{{V V}_{ph pH}}{v v} = = \frac{G G}{22 π π} {[[\frac{E E. + + {r r}^{22} F f}{b b + + 22 c c ((A A + + B B)) + + 44 dAB DAB + + 22 e e ((C C + + D D.)) + + 44 fCD f}]]}^{11 / / 22},,

A A = = cos cos ((\frac{22 π π cos cos ((θ θ))}{G G}))

B B = = cos cos ((\frac{22 π π sin sin ((θ θ))}{rG rG})),, - - - - - - ((77))

C C = = cos cos ((\frac{44 π π cos cos ((θ θ))}{G G})),,

D D. = = cos cos ((\frac{44 π π sin sin ((θ θ))}{rG rG}))

E E. = = \frac{11}{66} {α α}_{11} ((C C - - CD cd)) - - \frac{88}{33} {α α}_{22} ((A A - - AB AB)) + + \frac{55}{22} {α α}_{33} ((11 - - D D.)) + + {55 α α}_{44} ((B B - - D D.)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} D D.,,

F f = = \frac{11}{66} {β β}_{11} ((D D. - - CD cd)) - - \frac{88}{33} {β β}_{22} ((B B - - AB AB)) + + \frac{55}{22} {β β}_{33} ((11 - - C C)) + + {55 β β}_{44} ((A A - - C C)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} C C,,

in

Find the above optimal α _i , β _i , b, c, d, e, f, so that the difference between V _ph and v is the smallest;

When dx<dz, define r=Δz/Δx, G is defined as G=λ/Δz;

Rewrite the above formula (4) as

A=cos(k _x Δx),

B=cos(k _z Δz),

C＝cos(2k _x Δx),

D＝cos(2k _z Δz),

r=Δz/Δx

E E. = = \frac{11}{66} {α α}_{11} ((C C - - CD cd)) - - \frac{88}{33} {α α}_{22} ((A A - - AB AB)) + + \frac{55}{22} {α α}_{33} ((11 - - D D.)) + + {55 α α}_{44} ((B B - - D D.)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} D D.,,

F f = = \frac{11}{66} {β β}_{11} ((D D. - - CD cd)) - - \frac{88}{33} {β β}_{22} ((B B - - AB AB)) + + \frac{55}{22} {β β}_{33} ((11 - - C C)) + + {55 β β}_{44} ((A A - - C C)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} C C,,

Then the normalized phase velocity formula is

\frac{{V V}_{ph pH}}{v v} = = \frac{G G}{22 π π} {[[\frac{E E. + + {r r}^{22} F f}{b b + + 22 c c ((A A + + B B)) + + 44 dAB DAB + + 22 e e ((C C + + D D.)) + + 44 fCD f}]]}^{11 / / 22},,

A A = = cos cos ((\frac{22 π π cos cos ((θ θ))}{rG rG}))

B B = = cos cos ((\frac{22 π π sin sin ((θ θ))}{G G})),,

C C = = cos cos ((\frac{44 π π cos cos ((θ θ))}{rG rG})),,

D D. = = cos cos ((\frac{44 π π sin sin ((θ θ))}{G G}))

E E. = = \frac{11}{66} {α α}_{11} ((C C - - CD cd)) - - \frac{88}{33} {α α}_{22} ((A A - - AB AB)) + + \frac{55}{22} {α α}_{33} ((11 - - D D.)) + + {55 α α}_{44} ((B B - - D D.)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} D D.,,

F f = = \frac{11}{66} {β β}_{11} ((D D. - - CD cd)) - - \frac{88}{33} {β β}_{22} ((B B - - AB AB)) + + \frac{55}{22} {β β}_{33} ((11 - - C C)) + + {55 β β}_{44} ((A A - - C C)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} C C,,

Obtain the above optimal α _i , β _i , b, c, d, e, f to minimize the difference between V _ph and v.

4. A kind of frequency domain forward modeling method according to claim 3, it is characterized in that, in the process of above-mentioned optimum α _i , β _i , b, c, d, e, f, adopt least square method or Use the Gauss-Newton method to find the optimal α _i , β _i ,b,c,d,e,f.

5. A kind of frequency domain forward modeling method according to claim 1, is characterized in that, after reading in wavelet parameter and velocity model, also comprises, loads best matching layer PML boundary condition at the boundary of described velocity model, The wavefield used to absorb boundary reflections.

6. A frequency domain forward modeling device is characterized in that comprising,

The difference module is used to establish the difference formula in the following 17-point format,

in

{\overset{&OverBar; &OverBar;}{P P}}_{m m + + 11,, n no} = = {α α}_{22} {P P}_{m m + + 11,, n no} + + \frac{11 - - {α α}_{22}}{22} (({P P}_{m m + + 11,, n no + + 11} + + {P P}_{m m + + 11,, n no - - 11})),,

{\overset{&OverBar; &OverBar;}{P P}}_{m m - - 11,, n no} = = {α α}_{22} {P P}_{m m - - 11,, n no} + + \frac{11 - - {α α}_{22}}{22} (({P P}_{m m - - 11,, n no + + 11} + + {P P}_{m m - - 11,, n no - - 11})) - - - - - - ((33))

{\overset{&OverBar; &OverBar;}{P P}}_{m m + + 22,, n no} = = {α α}_{11} {P P}_{m m + + 22,, n no} + + \frac{11 - - {α α}_{11}}{22} (({P P}_{m m + + 22,, n no + + 22} + + {P P}_{m m + + 22,, n no - - 22})),,

{\overset{&OverBar; &OverBar;}{P P}}_{m m - - 22,, n no} = = {α α}_{11} {P P}_{m m - - 22,, n no} + + \frac{11 - - {α α}_{11}}{22} (({P P}_{m m - - 22,, n no + + 22} + + {P P}_{m m - - 22,, n no - - 22})),,

{\overset{&OverBar; &OverBar;}{P P}}_{m m,, n no} = = {α α}_{33} {P P}_{m m,, n no} + + {α α}_{44} (({P P}_{m m,, n no + + 11} + + {P P}_{m m,, n no - - 11})) + + \frac{11 - - {22 α α}_{44} - - {α α}_{33}}{22} (({P P}_{m m,, n no + + 22} + + {P P}_{m m,, n no - - 22})),,

and

{\overset{~ ~}{P P}}_{m m,, n no + + 11} = = {β β}_{22} {P P}_{m m,, n no + + 11} + + \frac{11 - - {β β}_{22}}{22} (({P P}_{m m + + 11,, n no + + 11} + + {P P}_{m m - - 11,, n no + + 11})),,

{\overset{~ ~}{P P}}_{m m,, n no - - 11} = = {β β}_{22} {P P}_{m m,, n no - - 11} + + \frac{11 - - {β β}_{22}}{22} (({P P}_{m m + + 11,, n no - - 11} + + {P P}_{m m - - 11,, n no - - 11})),, - - - - - - ((44))

{\overset{~ ~}{P P}}_{m m,, n no + + 22} = = {β β}_{11} {P P}_{m m,, n no + + 22} + + \frac{11 - - {β β}_{11}}{22} (({P P}_{m m + + 22,, n no + + 22} + + {P P}_{m m - - 22,, n no + + 22})),,

{\overset{~ ~}{P P}}_{m m,, n no - - 22} = = {β β}_{11} {P P}_{m m,, n no - - 22} + + \frac{11 - - {β β}_{11}}{22} (({P P}_{m m + + 22,, n no - - 22} + + {P P}_{m m - - 22,, n no - - 22})),,

{\overset{~ ~}{P P}}_{m m,, n no} = = {β β}_{33} {P P}_{m m,, n no} + + {β β}_{44} (({P P}_{m m + + 11,, n no} + + {P P}_{m m - - 11,, n no})) + + \frac{11 - - {22 β β}_{44} - - {β β}_{33}}{22} (({P P}_{m m + + 22,, n no} + + {P P}_{m m - - 22,, n no})),,

Sparse matrix module for constructing sparse matrices;

Read-in module, used to read in wavelet parameters and velocity models;

Single-frequency wave field module, used for frequency cycle to obtain single-frequency wave field;

The inverse Fourier transform module is used to calculate the inverse Fourier transform of all frequency wave fields;

The result output module to get the forward modeling result.

7. A kind of frequency domain forward modeling device according to claim 6, is characterized in that, also comprises weighting coefficient module, is connected between described difference module and sparse matrix module, is used for obtaining the maximum of described weighting coefficient Merit.

8. A kind of frequency domain forward modeling device according to claim 7, it is characterized in that, said weight coefficient module obtains the optimal value of weight coefficient further comprises, by wave equation (1) Starting, take the plane wave Substituting into formula (1), the dispersion equation can be obtained, where P ₀ is the initial amplitude value of the wave, i is the imaginary number unit,

A=cos(k _x Δx),

B＝cos(k _z Δz), (5)

C＝cos(2k _x Δx),

D＝cos(2k _z Δz),

E E. = = \frac{11}{66} {α α}_{11} ((C C - - CD cd)) - - \frac{88}{33} {α α}_{22} ((A A - - AB AB)) + + \frac{55}{22} {α α}_{33} ((11 - - D D.)) + + {55 α α}_{44} ((B B - - D D.)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} D D.,,

F f = = \frac{11}{66} {β β}_{11} ((D D. - - CD cd)) - - \frac{88}{33} {β β}_{22} ((B B - - AB AB)) + + \frac{55}{22} {β β}_{33} ((11 - - C C)) + + {55 β β}_{44} ((A A - - C C)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} C C,,

where k represents the wave number,

A=cos(k _x Δx),

B=cos(k _z Δz),

C＝cos(2k _x Δx), (6)

D＝cos(2k _z Δz),

E E. = = \frac{11}{66} {α α}_{11} ((C C - - CD cd)) - - \frac{88}{33} {α α}_{22} ((A A - - AB AB)) + + \frac{55}{22} {α α}_{33} ((11 - - D D.)) + + {55 α α}_{44} ((B B - - D D.)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} D D.,,

F f = = \frac{11}{66} {β β}_{11} ((D D. - - CD cd)) - - \frac{88}{33} {β β}_{22} ((B B - - AB AB)) + + \frac{55}{22} {β β}_{33} ((11 - - C C)) + + {55 β β}_{44} ((A A - - C C)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} C C,,

r=Δx/Δz,

So far, the normalized phase velocity formula is

\frac{{V V}_{ph pH}}{v v} = = \frac{G G}{22 π π} {[[\frac{E E. + + {r r}^{22} F f}{b b + + 22 c c ((A A + + B B)) + + 44 dAB DAB + + 22 e e ((C C + + D D.)) + + 44 fCD f}]]}^{11 / / 22},,

A A = = cos cos ((\frac{22 π π cos cos ((θ θ))}{G G}))

B B = = cos cos ((\frac{22 π π sin sin ((θ θ))}{rG rG})),, - - - - - - ((77))

C C = = cos cos ((\frac{44 π π cos cos ((θ θ))}{G G})),,

D D. = = cos cos ((\frac{44 π π sin sin ((θ θ))}{rG rG}))

E E. = = \frac{11}{66} {α α}_{11} ((C C - - CD cd)) - - \frac{88}{33} {α α}_{22} ((A A - - AB AB)) + + \frac{55}{22} {α α}_{33} ((11 - - D D.)) + + {55 α α}_{44} ((B B - - D D.)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} D D.,,

F f = = \frac{11}{66} {β β}_{11} ((D D. - - CD cd)) - - \frac{88}{33} {β β}_{22} ((B B - - AB AB)) + + \frac{55}{22} {β β}_{33} ((11 - - C C)) + + {55 β β}_{44} ((A A - - C C)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} C C,,

in

When dx<dz, define r=Δz/Δx, G is defined as G=λ/Δz;

Rewrite the above formula (4) as

A=cos(k _x Δx),

B=cos(k _z Δz),

C＝cos(2k _x Δx),

D＝cos(2k _z Δz),

r=Δz/Δx

E E. = = \frac{11}{66} {α α}_{11} ((C C - - CD cd)) - - \frac{88}{33} {α α}_{22} ((A A - - AB AB)) + + \frac{55}{22} {α α}_{33} ((11 - - D D.)) + + {55 α α}_{44} ((B B - - D D.)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} D D.,,

F f = = \frac{11}{66} {β β}_{11} ((D D. - - CD cd)) - - \frac{88}{33} {β β}_{22} ((B B - - AB AB)) + + \frac{55}{22} {β β}_{33} ((11 - - C C)) + + {55 β β}_{44} ((A A - - C C)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} C C,,

Then the normalized phase velocity formula is

\frac{{V V}_{ph pH}}{v v} = = \frac{G G}{22 π π} {[[\frac{E E. + + {r r}^{22} F f}{b b + + 22 c c ((A A + + B B)) + + 44 dAB DAB + + 22 e e ((C C + + D D.)) + + 44 fCD f}]]}^{11 / / 22},,

A A = = cos cos ((\frac{22 π π cos cos ((θ θ))}{rG rG})),,

C C = = cos cos ((\frac{44 π π cos cos ((θ θ))}{rG rG})),,

D D. = = cos cos ((\frac{44 π π sin sin ((θ θ))}{G G}))

E E. = = \frac{11}{66} {α α}_{11} ((C C - - CD cd)) - - \frac{88}{33} {α α}_{22} ((A A - - AB AB)) + + \frac{55}{22} {α α}_{33} ((11 - - D D.)) + + {55 α α}_{44} ((B B - - D D.)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} D D.,,

F f = = \frac{11}{66} {β β}_{11} ((D D. - - CD cd)) - - \frac{88}{33} {β β}_{22} ((B B - - AB AB)) + + \frac{55}{22} {β β}_{33} ((11 - - C C)) + + {55 β β}_{44} ((A A - - C C)) - - \frac{88}{33} AB AB + + \frac{11}{66} CD cd + + \frac{55}{22} C C,,

9. A kind of frequency domain forward modeling device according to claim 8, it is characterized in that, in the above-mentioned process of seeking optimal α _i , β _i , b, c, d, e, f of said weighting coefficient module, Use the least squares method or the Gauss-Newton method to find the optimal α _i , β _i , b, c, d, e, f.

10. A frequency-domain forward modeling device according to claim 1, further comprising an absorbing boundary module, connected between the read-in module and the single-frequency wavefield module, for The boundary-loaded best-matched-layer PML boundary condition is used to absorb wavefields reflected by the boundary.