CN103149586A

CN103149586A - Tilted layered viscoelasticity dielectric medium wave field forward modelling method

Info

Publication number: CN103149586A
Application number: CN2013100443842A
Authority: CN
Inventors: 高静怀; 汪超; 王大兴
Original assignee: Xian Jiaotong University
Current assignee: Xian Jiaotong University
Priority date: 2013-02-04
Filing date: 2013-02-04
Publication date: 2013-06-12
Anticipated expiration: 2033-02-04
Also published as: CN103149586B

Abstract

The invention discloses a wave field forward simulation method in an inclined layered viscoelastic medium. The method aims at the geological model as a viscoelastic layered medium with an inclined interface, decomposes the spherical wave excited by a point source into a plane wave, and utilizes the vector wave The equation does not depend on the property of the coordinate system, and the wave propagation is discussed alternately in the fixed coordinate system and the moving coordinate system which changes with the inclined interface. First, the reflection and transmission coefficients of a plane wave on a single inclined interface are deduced from the traditional reflection and transmission coefficients, and the reflection and transmission waves on the interface are obtained; then the reflection and transmission coefficients of plane waves passing through multiple inclined interfaces are obtained recursively. Finally, the wave field excited by the point source is obtained by merging the propagating plane waves. When merging, a fast oscillation integral calculation method and an integral path suitable for viscoelastic media are used.

Description

A Forward Modeling Method for Wavefield in Inclined Layered Viscoelastic Media

技术领域technical field

本发明属于地震勘探技术领域，涉及一种模拟方法，尤其是一种倾斜层状粘弹性介质中波场正演模拟方法，其是一种针对含倾斜界面的层状粘弹性介质中地震波场高效精确的正演模拟技术。The invention belongs to the technical field of seismic exploration, and relates to a simulation method, in particular to a wave field forward simulation method in an inclined layered viscoelastic medium. Accurate forward modeling technology.

背景技术Background technique

波场模拟在油气勘探中占有重要地位，尽管已存在多种类型的波场模拟方法，然而，针对要解决的问题的要求，发展波场模拟理论及方法一直是很活跃的研究领域。另外，在油气勘探等领域，为了解释地震图，或进行数据反演，常常需要仅模拟指定类型的波，譬如仅模拟一次反射波，或仅模拟转换波等。常用的波场模拟方法有多种类型，如基于射线理论的方法，数值解方法，广义反射率法以及其推广而得到的反/透射矩阵法等。这些方法各有其优缺点。Wave field simulation plays an important role in oil and gas exploration. Although there are many types of wave field simulation methods, the development of wave field simulation theory and methods has always been a very active research field according to the requirements of the problems to be solved. In addition, in fields such as oil and gas exploration, in order to interpret seismograms or perform data inversion, it is often necessary to simulate only specified types of waves, such as simulating only one reflection wave, or only simulating converted waves, etc. There are many types of commonly used wave field simulation methods, such as methods based on ray theory, numerical solution methods, generalized reflectivity methods, and reflection/transmission matrix methods obtained by its extension. Each of these methods has its advantages and disadvantages.

射线方法能适用于较复杂的介质模型，可根据事先给定的射线码模拟感兴趣的类型的波，且计算速度快；不足之处是，它只是波动方程的高频近似，用于模拟复杂储层的地震响应具有局限性。采用数值法求解波动方程，是复杂介质中波场模拟的有效途径。常用的数值方法有：有限差分法，有限元法，伪谱法和谱元法等，这类方法可模拟地震波全波场，精度高。但计算量巨大，不能像射线类方法，有选择性的仅模拟我们感性趣的类型的波。The ray method can be applied to more complex medium models, and can simulate waves of interest according to a given ray code, and the calculation speed is fast; the disadvantage is that it is only a high-frequency approximation of the wave equation, which is used to simulate complex The seismic response of reservoirs has limitations. Using numerical method to solve wave equation is an effective way to simulate wave field in complex media. Commonly used numerical methods include: finite difference method, finite element method, pseudospectral method and spectral element method, etc. These methods can simulate the full wave field of seismic waves with high precision. However, the amount of calculation is huge, and it cannot selectively simulate only the types of waves that we are interested in like the ray method.

层状介质中波传播理论及波场模拟技术一直受到人们的重视。Aki和Müller等人发展的反射率法，能快速精确地模拟水平层状粘弹性介质的全波场或指定类型的波。Berkhout提出了波场外推法可用于带倾斜界面的弹性介质波场模拟。Chen等人发展的全局广义反/透射矩阵法和边界元法，能模拟带不规则界面的层状弹性介质中全波场。但这些方法不适用于含倾斜界面的层状粘弹性介质。The wave propagation theory and wave field simulation technology in layered media have been paid attention to by people all the time. The reflectivity method developed by Aki and Müller et al. can quickly and accurately simulate the full wave field or specified types of waves in horizontal layered viscoelastic media. Berkhout proposed that the wave field extrapolation method can be used for wave field simulation of elastic media with inclined interfaces. The global generalized reflection/transmission matrix method and boundary element method developed by Chen et al. can simulate the full wave field in layered elastic media with irregular interfaces. However, these methods are not suitable for layered viscoelastic media with inclined interfaces.

发明内容Contents of the invention

针对上述问题，本发明的目的在于针对含倾斜波阻抗界面的层状粘弹性介质，提出一种倾斜层状粘弹性介质中波场正演模拟方法，该方法可以选择性模拟感兴趣的各种类型的波，如一次反射P波、S波和转换波，且模拟的波场非常精确。In view of the above problems, the object of the present invention is to propose a wavefield forward modeling method in an inclined layered viscoelastic medium for a layered viscoelastic medium containing an inclined wave impedance interface, which can selectively simulate various interested Types of waves, such as once reflected P wave, S wave and converted wave, and the simulated wave field is very accurate.

本发明的目的是通过以下技术方案来解决的：The purpose of the present invention is solved by the following technical solutions:

该种倾斜层状粘弹性介质中波场正演模拟方法，包括以下步骤：The wave field forward modeling method in an inclined layered viscoelastic medium comprises the following steps:

1）将点源激励的球面波分解为平面波；1) Decompose the spherical wave excited by the point source into a plane wave;

2）将平面波在倾斜层状粘弹性介质中传播，遇到界面时进行如下处理：以该界面为

轴，界面与全局坐标系Z轴的交点为原点建立局部坐标系，垂直于界面向下，将全局坐标系中的入射波转换到局部坐标系中；在局部坐标系下求得反射波和透射波；再将局部坐标系下的反射波和透射波转换到全局坐标系下，得到全局坐标系下的反射系数和透射系数，以及反射波和透射波；在每个界面上根据需要选取波的类型，P波入射时能够选取的波包括反射PP波、反射PS波、透射PP波、透射PS波；S波入射时能够选取的波包括反射SP波、反射SS波、透射SP波、透射SS波；从而模拟所选取的波类型，得到模拟的平面波；2) Propagate the plane wave in the inclined layered viscoelastic medium, and when encountering the interface, proceed as follows: take the interface as

Axis, the intersection point of the interface and the Z axis of the global coordinate system is the origin to establish a local coordinate system, Vertically downward to the interface, convert the incident wave in the global coordinate system to the local coordinate system; obtain the reflected wave and transmitted wave in the local coordinate system; then convert the reflected wave and transmitted wave in the local coordinate system to the global coordinate system system, the reflection coefficient and transmission coefficient in the global coordinate system, as well as the reflected wave and the transmitted wave are obtained; on each interface, the type of wave is selected according to the needs, and the wave that can be selected when the P wave is incident includes reflected PP wave, reflected PS wave , transmitted PP wave, transmitted PS wave; the waves that can be selected when the S wave is incident include reflected SP wave, reflected SS wave, transmitted SP wave, and transmitted SS wave; thereby simulating the selected wave type and obtaining a simulated plane wave;

3）由步骤2）模拟得到的平面波根据下式合成点源激发的波场：3) The plane wave simulated by step 2) synthesizes the wave field excited by the point source according to the following formula:

$P P = = \frac{11}{22 π π} {&Integral; &Integral;}_{Γ Γ} {Σ Σ}_{i i = = 11}^{N N - - 11} [[{R R}_{i i} ((p p)) \cdot &Center Dot; \frac{exp exp [[jω jω (({p p}_{11 i i} x x - - {q q}_{11 i i} z z))]]}{- - 22 jq jq}]] dp dp$

式中，Γ为积分路径；N为平面波的个数；R_i(p)为平面波在多个界面上的透射和反射系数综合；j是虚数单位，ω是角频率；p_1i、q_1i分别为第i个界面反射回第1层的平面波的水平和垂直慢度；p、q分别为点源处平面波的水平和垂直慢度，x和z表示位置坐标。In the formula, Γ is the integral path; N is the number of plane waves; R _i (p) is the synthesis of transmission and reflection coefficients of plane waves on multiple interfaces; j is the imaginary unit, ω is the angular frequency; p _1i and q _1i are respectively is the horizontal and vertical slowness of the plane wave reflected from the i-th interface back to the first layer; p and q are the horizontal and vertical slowness of the plane wave at the point source, respectively, and x and z represent the position coordinates.

进一步的，上述步骤2）中，模拟震源为P波时的一次反射P波，设震源和检波器均在第一层，将震源处入射波逐层向下传播的递推步骤为：Further, in the above step 2), the simulated seismic source is a reflected P wave when the seismic source is a P wave. Assuming that the seismic source and the geophone are both on the first floor, the recursive steps for propagating the incident wave at the seismic source downward layer by layer are as follows:

①每遇到一个界面i，即以该界面为轴，界面与全局坐标系Z轴的交点为原点建立局部坐标系

垂直于界面向下；将全局坐标系中的入射波φ_i转换到局部坐标系中得

① Whenever an interface i is encountered, the interface is taken as Axis, the intersection of the interface and the Z axis of the global coordinate system is the origin to establish a local coordinate system

Downward perpendicular to the interface; transform the incident wave φ _i in the global coordinate system into the local coordinate system to get

②在局部坐标系下求得反射系数和透射系数

以及反射波

和透射波 ② Obtain the reflection coefficient in the local coordinate system and transmission coefficient

and reflected waves

and transmitted waves

③再将局部坐标系下的反射波和透射波转换到全局坐标系下，得到全局坐标系下的反射系数和透射系数

以及反射波φ′_i，i和透射波φ_i+1；③ Then convert the reflected wave and transmitted wave in the local coordinate system to the global coordinate system to obtain the reflection coefficient in the global coordinate system and transmission coefficient

And reflected wave φ′ _{i, i} and transmitted wave φ _i+1 ;

④反射波φ′_i，i逐层向上递推到第一层，遇到界面时执行①-③步，但只取透射波；④Reflected wave φ′ _{i, i} is recursively pushed up to the first layer layer by layer, and steps ①-③ are performed when encountering an interface, but only the transmitted wave is taken;

⑤透射波φ_i+1继续向下递推，遇到下一个界面时回到步骤①。⑤ The transmitted wave φ _i+1 continues to deduce downwards, and returns to step ① when the next interface is encountered.

进一步，在上述步骤2）中，模拟其它类型的波采用与模拟震源为P波时的一次反射P波的相同方法和步骤进行。Further, in the above step 2), other types of waves are simulated using the same method and steps as the one-time reflection P-wave simulation when the seismic source is P-wave.

在上述步骤3）中，式中积分方法采用Frazer和Gettrust提出的GFM积分法。In the above step 3), the integral method in the formula adopts the GFM integral method proposed by Frazer and Gettrust.

本发明具有以下几点有益效果：The present invention has the following beneficial effects:

（1）本发明适用于带倾斜界面的层状粘弹性介质，各界面的倾斜方向和倾斜角度可以不同，只需满足界面不相交的条件。(1) The present invention is applicable to layered viscoelastic media with inclined interfaces, and the inclined direction and angle of each interface can be different, as long as the interface does not intersect.

（2）本发明能正确模拟粘弹性介质中波的衰减特征，相对于有限差分法，具有较高的计算效率，更适合层较厚、且层间参数变化剧烈情况。(2) The present invention can correctly simulate the attenuation characteristics of waves in viscoelastic media. Compared with the finite difference method, it has higher calculation efficiency and is more suitable for thicker layers and drastic changes in interlayer parameters.

（3）本发明非常灵活，可以分别模拟各种波，如一次反射P波、SV波、SH波和PS、SP转换波，也可模拟我们感兴趣的多次波，只需将希望模拟的波类型和相应的层次以编码形式给定，入射波入射角超过临界角时会出现首波。(3) The present invention is very flexible and can simulate various waves, such as primary reflection P wave, SV wave, SH wave and PS, SP conversion wave, and can also simulate multiple waves that we are interested in. The wave types and corresponding levels are given in coded form, and the first wave will appear when the incident angle of the incident wave exceeds the critical angle.

（4）本发明模拟的波场非常接近解析解，可为数值方法解提供一个对比参照。所以本发明在地震资料解释和反演中有重要应用价值和应用前景。(4) The wave field simulated by the present invention is very close to the analytical solution, which can provide a comparison reference for the numerical method solution. Therefore, the present invention has important application value and application prospect in seismic data interpretation and inversion.

附图说明Description of drawings

图1是本发明针对的倾斜层状粘弹性介质模型；Fig. 1 is the inclined layered viscoelastic medium model aimed at by the present invention;

图2是波入射到界面上时的反射和透射示意图；Figure 2 is a schematic diagram of reflection and transmission when the wave is incident on the interface;

图3是本发明中平面波合成点源激励波场时的积分路径示意图。Fig. 3 is a schematic diagram of the integration path when the plane wave synthesis point source excites the wave field in the present invention.

具体实施方式Detailed ways

本发明针对的地质模型为含倾斜界面的粘弹性层状介质，将点源激发的球面波分解为平面波，利用矢量波动方程不依赖于坐标系这一性质，在固定坐标系和随倾斜界面变化的动坐标系下交替地讨论波的传播。首先由传统的反射和透射系数推导出单个倾斜界面上平面波的反射和透射系数，得到该界面上的反射及透射波；然后递推求得平面波经过多个倾斜界面时的反射和透射系数。最后由传播后的各平面波合并得到点源激发波场，合并时采用了快速的振荡积分计算方法和一条适合粘弹性介质的积分路径。The geological model targeted by the present invention is a viscoelastic layered medium with an inclined interface, decomposes the spherical wave excited by the point source into a plane wave, and utilizes the property that the vector wave equation does not depend on the coordinate system. Wave propagation is discussed alternately in the moving coordinate system of . First, the reflection and transmission coefficients of a plane wave on a single inclined interface are deduced from the traditional reflection and transmission coefficients, and the reflection and transmission waves on the interface are obtained; then the reflection and transmission coefficients of plane waves passing through multiple inclined interfaces are obtained recursively. Finally, the wave field excited by the point source is obtained by merging the propagating plane waves. When merging, a fast oscillation integral calculation method and an integral path suitable for viscoelastic media are used.

下面结合附图和实施例对本发明进行详细的描述。The present invention will be described in detail below in conjunction with the accompanying drawings and embodiments.

首先本发明的主要符号定义如表1所示：First the main symbols of the present invention are defined as shown in Table 1:

表1本发明中主要符号约定Main symbol convention in the present invention in table 1

图1为倾斜层状粘弹性介质模型示意图，各层介质的密度为ρ_i(i＝1,2,…,N-1)、品质因子为Q_i、P波速度为α_i、S波速度为β_i，各倾斜界面的倾角为θ_i，各倾斜界面与垂直方向的Z轴交点为z_i。点源在第一次层，检波器可以在任意位置。以水平方向为X轴，垂直向下方向为Z轴，建立全局坐标系XOZ。Fig. 1 is a schematic diagram of an inclined layered viscoelastic medium model. The density of each layer medium is ρ _i (i=1,2,…,N-1), the quality factor is Q _i , the P-wave velocity is α _i , and the S-wave velocity is β _i , the inclination angle of each inclined interface is θ _i , and the intersection point of each inclined interface with the vertical Z-axis is z _i . The point source is in the first layer, and the detector can be in any position. With the horizontal direction as the X axis and the vertical downward direction as the Z axis, a global coordinate system XOZ is established.

1）将点源激励的球面波分解为平面波1) Decompose the spherical wave excited by the point source into a plane wave

根据Weyl积分，均匀介质中点源激励的球面波表示为一系列平面波的叠加，二维情况下每个平面波的表达式为：According to the Weyl integral, the spherical wave excited by a point source in a homogeneous medium is expressed as a superposition of a series of plane waves, and the expression of each plane wave in two dimensions is:

φ=exp[jω(px+qz)]φ=exp[jω(px+qz)]

上式中省略了时间因子exp(-jωt)，j是虚数单位，ω是角频率，p是水平慢度，q是垂直慢度，x和z表示位置坐标。In the above formula, the time factor exp(-jωt) is omitted, j is the imaginary unit, ω is the angular frequency, p is the horizontal slowness, q is the vertical slowness, and x and z represent the position coordinates.

轴，界面与全局坐标系Z轴的交点为原点建立局部坐标系，

垂直于界面向下，将全局坐标系中的入射波转换到局部坐标系中；在局部坐标系下求得反射波和透射波；再将局部坐标系下的反射波和透射波转换到全局坐标系下，得到全局坐标系下的反射系数和透射系数，以及反射波和透射波。在每个界面上可根据需要选取波的类型，P波入射时能够选取的波包括反射PP波、反射PS波、透射PP波、透射PS波；S波入射时能够选取的波包括反射SP波、反射SS波、透射SP波、透射SS波；从而仅模拟我们感兴趣的波类型。例如需模拟震源为P波时的一次反射P波，设震源和检波器均在第一层，将震源处入射波逐层向下传播的递推步骤为：2) Propagate the plane wave in the inclined layered viscoelastic medium, and when encountering the interface, proceed as follows: take the interface as

Axis, the intersection point of the interface and the Z axis of the global coordinate system is the origin to establish a local coordinate system,

Vertically downward to the interface, convert the incident wave in the global coordinate system to the local coordinate system; obtain the reflected wave and transmitted wave in the local coordinate system; then convert the reflected wave and transmitted wave in the local coordinate system to the global coordinate system In the system, the reflection coefficient and transmission coefficient in the global coordinate system, as well as the reflected wave and transmitted wave are obtained. On each interface, the type of wave can be selected according to the needs. The waves that can be selected when the P wave is incident include reflected PP wave, reflected PS wave, transmitted PP wave, and transmitted PS wave; the waves that can be selected when the S wave is incident include reflected SP wave , reflected SS waves, transmitted SP waves, transmitted SS waves; thereby simulating only the wave types we are interested in. For example, it is necessary to simulate a reflected P wave when the source is a P wave. Assuming that the source and the receiver are both on the first floor, the recursive steps for propagating the incident wave at the source down layer by layer are as follows:

①每遇到一个界面i，即以该界面为

轴，界面与全局坐标系Z轴的交点为原点建立局部坐标系垂直于界面向下。将全局坐标系中的入射波φ_i转换到局部坐标系中得

① Whenever an interface i is encountered, the interface is taken as

Axis, the intersection of the interface and the Z axis of the global coordinate system is the origin to establish a local coordinate system Vertically down the interface. Convert the incident wave φ _i in the global coordinate system to the local coordinate system to get

②在局部坐标系下求得反射系数

和透射系数

以及反射波

和透射波

② Obtain the reflection coefficient in the local coordinate system

and transmission coefficient

and reflected waves

and transmitted waves

③再将局部坐标系下的反射波和透射波转换到全局坐标系下，得到全局坐标系下的反射系数

和透射系数

以及反射波φ′_i，i和透射波φ_i+1；③ Then convert the reflected wave and transmitted wave in the local coordinate system to the global coordinate system to obtain the reflection coefficient in the global coordinate system

and transmission coefficient

And reflected wave φ′ _{i, i} and transmitted wave φ _i+1 ;

④反射波φ′_i，i逐层向上递推到第一层，遇到界面时执行①-③步，但只取透射波（上行波）；④Reflected wave φ′ _{i, i} is recursively pushed up to the first layer layer by layer, and steps ①-③ are performed when encountering the interface, but only the transmitted wave (upgoing wave) is taken;

3）由平面波合成点源激发的波场；3) The wave field excited by a plane wave synthetic point source;

$P P = = \frac{11}{22 π π} {&Integral; &Integral;}_{Γ Γ} {Σ Σ}_{i i = = 11}^{N N - - 11} [[{R R}_{i i} ((p p)) \cdot \cdot \frac{exp exp [[jω jω (({p p}_{11 i i} x x - - {q q}_{11 i i} z z))]]}{- - 22 jq jq}]] dp dp$

其中，R_i(p)为平面波在多个界面上的透射和反射系数综合，如模拟一次P-P反射波时 $R_{i} (p) = Π_{n = 1}^{i - 1} T_{pp, n}^{u} \cdot R_{pp, i}^{d} \cdot Π_{m = 1}^{i - 1} T_{pp, m}^{d};$ p和q为震源处入射平面波的水平和垂直慢度；p_1i和q_1i为平面波经第i个界面反射后到达检波器所在层时的水平和垂直慢度，它们是入射波水平慢度p的函数。Γ为积分路径，本发明选择如图2所示的积分路径，1/α_i和1/β_i是P波和S波的慢度，路径在第一象限中的拐点B的实部应稍大于max(1/β_i)，对单个接收器，θ=tan^-1(xz_r-z_s|)时积分收敛最快（x、z_r、z_s分别是偏移距和接收器与源的z轴坐标），如是多个接收器同时接收，可取某个折中值。线段OB的倾角应使得OB段非常靠近第一象限内最小极点的下方。积分方法采用Frazer和Gettrust提出GFM(Generalization ofFilon’s Method)积分法。Among them, R _i (p) is the synthesis of transmission and reflection coefficients of plane waves on multiple interfaces, such as when simulating a PP reflection wave $R_{i} (p) = Π_{no = 1}^{i - 1} T_{pp, no}^{u} &Center Dot; R_{pp, i}^{d} &Center Dot; Π_{m = 1}^{i - 1} T_{pp, m}^{d};$ p and q are the horizontal and vertical slowness of the incident plane wave at the source; p _1i and q _1i are the horizontal and vertical slownesses of the plane wave when it reaches the layer where the detector is located after being reflected by the i-th interface, and they are the horizontal slowness of the incident wave p The function. Γ is an integral path, the present invention selects the integral path as shown in Figure 2, 1/α _i and 1/β _i are the slowness of P wave and S wave, and the real part of the inflection point B of path in the first quadrant should be slightly is greater than max(1/β _i ), for a single receiver, when θ=tan ^-1 (xz _r -z _s |) the integral converges fastest (x, z _r , z _s are the offset distance and the receiver and source z-axis coordinates), if multiple receivers receive at the same time, a certain compromise value can be taken. The inclination of line segment OB should be such that segment OB is very close below the smallest pole in the first quadrant. The integral method adopts the GFM (Generalization of Filon's Method) integral method proposed by Frazer and Gettrust.

步骤2）中入射波在界面上产生的反射波和透射波的具体计算方法为：The specific calculation method of the reflected wave and transmitted wave generated by the incident wave on the interface in step 2) is:

如图2，设有一水平慢度为p₁的P波φ₁从原点出发在介质I中向下传播，振幅为1，则该平面波可写为：As shown in Figure 2, suppose a P-wave φ ₁ with horizontal slowness p ₁ propagates downward in the medium I from the origin, with an amplitude of 1, then the plane wave can be written as:

φ₁=exp[jω(p₁x+q₁z)],φ ₁ =exp[jω(p ₁ x+q ₁ z)],

将入射P波由全局坐标系转换到局部坐标系，得到局部坐标系下的入射P波为：Transform the incident P wave from the global coordinate system to the local coordinate system to obtain the local coordinate system The incident P wave under is:

${\overset{^^}{φ φ}}_{11} = = {A A}_{11} exp exp [[jω jω (({\overset{^^}{p p}}_{11} \overset{^^}{x x} + + {\overset{^^}{q q}}_{11} \overset{^^}{z z}))]],,$

其中：in:

A₁=exp[jω(p₁x₀+q₁z₀)], $[\begin{matrix} {\hat{p}}_{1} \\ {\hat{q}}_{1} \end{matrix}] = C [\begin{matrix} p_{1} \\ q 1 \end{matrix}],$ $C = [\begin{matrix} \cos θ & \sin θ \\ - \sin θ & \cos θ \end{matrix}] .$ A ₁ =exp[jω(p ₁ x ₀ +q ₁ z ₀ )], $[\begin{matrix} {\hat{p}}_{1} \\ {\hat{q}}_{1} \end{matrix}] = C [\begin{matrix} p_{1} \\ q 1 \end{matrix}],$ $C = [\begin{matrix} \cos θ & \sin θ \\ - \sin θ & \cos θ \end{matrix}] .$

在局部坐标系

下，P波

从介质I向下入射到界面，则在界面处会产生反射回介质I的上行P波

和SV波

以及透射到介质II中的下行P波

和SV波

各反射波和透射波的表达式分别为：in the local coordinate system

Down, P wave

From the medium I down to the interface, there will be an upward P wave reflected back to the medium I at the interface

and SV waves

and the downgoing P-wave transmitted into medium II

and SV wave

The expressions of each reflected wave and transmitted wave are respectively:

${\overset{^^}{φ φ}}_{11}^{' '} = = {\overset{^^}{R R}}_{pp pp}^{d d} \cdot &Center Dot; {A A}_{11} exp exp [[jω jω (({\overset{^^}{p p}}_{11} \overset{^^}{x x} - - {\overset{^^}{q q}}_{11 p p} \overset{^^}{z z}))]],,$

${\overset{^^}{ψ ψ}}_{11}^{' '} = = {\overset{^^}{R R}}_{ps ps}^{d d} \cdot &Center Dot; {A A}_{11} exp exp [[jω jω (({\overset{^^}{p p}}_{11} \overset{^^}{x x} - - {\overset{^^}{q q}}_{11 s the s} \overset{^^}{z z}))]],,$

${\overset{^^}{φ φ}}_{22} = = {\overset{^^}{T T}}_{pp pp}^{d d} \cdot \cdot {A A}_{11} exp exp [[jω jω (({\overset{^^}{p p}}_{11} \overset{^^}{x x} + + {\overset{^^}{q q}}_{22 p p} \overset{^^}{z z}))]],,$

${\overset{^^}{ψ ψ}}_{22} = = {\overset{^^}{T T}}_{ps ps}^{d d} \cdot &Center Dot; {A A}_{11} exp exp [[jω jω (({\overset{^^}{p p}}_{11} \overset{^^}{x x} + + {\overset{^^}{q q}}_{22 s the s} \overset{^^}{z z}))]],,$

其中

是反射和透射系数。局部坐标系下各波水平慢度相等，垂直慢度由下式给出：in

are the reflection and transmission coefficients. The horizontal slowness of each wave in the local coordinate system is equal, and the vertical slowness is given by the following formula:

${\overset{^^}{q q}}_{ip ip} = = \sqrt{{α α}_{i i}^{- - 22} - - {\overset{^^}{p p}}_{11}^{22}} . .$ ${\overset{^^}{q q}}_{is is} = = \sqrt{{β β}_{i i}^{- - 22} - - {\overset{^^}{p p}}_{11}^{22}} . . ((i i = = 1,2 1,2))$

相对于局部坐标系，界面是水平的，这种情况的反射和透射系数计算方法可参照相关文献（Muller，1985）。将所有反射波和透射波由局部坐标系转换到全局坐标系，整理变换后各种波的表达式则得到全局坐标系下各反射波和透射波为：Relative to the local coordinate system, the interface is horizontal, the reflection and transmission coefficients for this case The calculation method can refer to relevant literature (Muller, 1985). Transform all reflected waves and transmitted waves from the local coordinate system to the global coordinate system, and sort out and transform the expressions of various waves to obtain the reflected waves and transmitted waves in the global coordinate system as:

${φ φ}_{11}^{' '} = = {R R}_{pp pp}^{d d} \cdot \cdot exp exp [[jω jω (({p p}_{11 p p} x x - - {q q}_{11 p p} z z))]],,$

${ψ ψ}_{11}^{' '} = = {R R}_{ps ps}^{d d} \cdot \cdot exp exp [[jω jω (({p p}_{11 s the s} x x - - {q q}_{11 s the s} z z))]],,$

${φ φ}_{22} = = {T T}_{pp pp}^{d d} \cdot \cdot exp exp [[jω jω (({p p}_{22 p p} x x + + {q q}_{22 p p} z z))]],,$

${ψ ψ}_{22} = = {T T}_{ps ps}^{d d} \cdot \cdot exp exp [[jω jω (({p p}_{22 s the s} x x + + {q q}_{22 s the s} z z))]],,$

式中：In the formula:

${R R}_{pp pp}^{d d} = = {\overset{^^}{R R}}_{pp pp}^{d d} \cdot &Center Dot; {A A}_{11} exp exp [[- - jω jω (({p p}_{11 p p} {x x}_{00} - - {q q}_{11 p p} {z z}_{00}))]],,$

${R R}_{ps ps}^{d d} = = {\overset{^^}{R R}}_{ps ps}^{d d} \cdot &Center Dot; {A A}_{11} exp exp [[- - jω jω (({p p}_{11 s the s} {x x}_{00} - - {q q}_{11 s the s} {z z}_{00}))]],,$

${T T}_{pp pp}^{d d} = = {\overset{^^}{T T}}_{pp pp}^{d d} \cdot \cdot {A A}_{11} exp exp [[- - jω jω (({p p}_{22 p p} {x x}_{00} + + {q q}_{22 p p} {z z}_{00}))]],,$

${T T}_{ps ps}^{d d} = = {\overset{^^}{T T}}_{ps ps}^{d d} \cdot \cdot {A A}_{11} exp exp [[- - jω jω (({p p}_{22 s the s} {x x}_{00} + + {q q}_{22 s the s} {z z}_{00}))]],,$

$[\begin{matrix} {p p}_{11 p p} \\ {q q}_{11 p p} \end{matrix}] = = C C [\begin{matrix} {\overset{^^}{p p}}_{11} \\ {\overset{^^}{q q}}_{11 p p} \end{matrix}] . .$ $[\begin{matrix} {p p}_{11 s the s} \\ {q q}_{11 s the s} \end{matrix}] = = C C [\begin{matrix} {\overset{^^}{p p}}_{11} \\ {\overset{^^}{q q}}_{11 s the s} \end{matrix}] . .$ $[\begin{matrix} {p p}_{22 p p} \\ {q q}_{22 p p} \end{matrix}] = = {C C}^{- - 11} [\begin{matrix} {\overset{^^}{p p}}_{22} \\ {\overset{^^}{q q}}_{22 p p} \end{matrix}] . .$ $[\begin{matrix} {p p}_{22 s the s} \\ {q q}_{22 s the s} \end{matrix}] = = {C C}^{- - 11} [\begin{matrix} {\overset{^^}{p p}}_{22} \\ {\overset{^^}{q q}}_{22 s the s} \end{matrix}] . .$

其中

即分别为全局坐标系XOZ下粘弹性倾斜界面的P-P反射系数、P-S反射系数、P-P透射系数和P-S透射系数。in

That is, PP reflection coefficient, PS reflection coefficient, PP transmission coefficient and PS transmission coefficient of the viscoelastic inclined interface in the global coordinate system XOZ, respectively.

对于波由下方向上入射到界面的情况，设全局坐标系下有一P波φ′₂在介质II中向上传播：For the case where the wave is incident upward on the interface from below, it is assumed that there is a P wave φ′ ₂ propagating upward in the medium II in the global coordinate system:

φ′₂=exp[jω(p₂x-q₂z)].φ′ ₂ =exp[jω(p ₂ xq ₂ z)].

则在界面处会产生反射回介质II的下行P波φ₂和SV波Ψ₂，以及透射到介质I中的上行P波φ′₁和SV波Ψ′₁。按照P波从界面上方入射时的处理流程：先将入射波φ′₂转换到局部坐标系下，在局部坐标系里求得各反射波和透射波的传统反射和透射系数；再将局部坐标系下各反射波和透射波转换到全局坐标系，得到全局坐标系下反射波和透射波的表达式分别为：Then at the interface, there will be downgoing P wave φ ₂ and SV wave Ψ ₂ reflected back to medium II, and upgoing P wave φ′ ₁ and SV wave Ψ′ ₁ transmitted into medium I. According to the processing flow when the P wave is incident from above the interface: first transform the incident wave φ′ ₂ into the local coordinate system, and obtain the traditional reflection and transmission coefficients of each reflected wave and transmitted wave in the local coordinate system; then convert the local coordinate The reflected waves and transmitted waves in the global coordinate system are converted to the global coordinate system, and the expressions of the reflected waves and transmitted waves in the global coordinate system are respectively:

${φ φ}_{22} = = {R R}_{pp pp}^{u u} \cdot \cdot exp exp [[jω jω (({p p}_{22 p p} x x + + {q q}_{22 p p} z z))]],,$

${ψ ψ}_{22} = = {R R}_{ps ps}^{u u} \cdot \cdot exp exp [[jω jω (({p p}_{22 s the s} x x + + {q q}_{22 s the s} z z))]],,$

${φ φ}_{11}^{' '} = = {T T}_{pp pp}^{u u} \cdot \cdot exp exp [[jω jω (({p p}_{11 p p} x x - - {q q}_{11 p p} z z))]],,$

${ψ ψ}_{11}^{' '} = = {T T}_{ps ps}^{u u} \cdot &Center Dot; exp exp [[jω jω (({p p}_{11 s the s} x x - - {q q}_{11 s the s} z z))]],,$

其中：in:

${R R}_{pp pp}^{u u} = = {\overset{^^}{R R}}_{pp pp}^{u u} \cdot \cdot {A A}^{' '} exp exp [[- - jω jω (({p p}_{22 p p} {x x}_{00} + + {q q}_{22 p p} {z z}_{00}))]],,$

${R R}_{ps ps}^{u u} = = {\overset{^^}{R R}}_{ps ps}^{u u} \cdot &Center Dot; {A A}^{' '} exp exp [[- - jω jω (({p p}_{22 s the s} {x x}_{00} + + {q q}_{22 s the s} {z z}_{00}))]],,$

${T T}_{pp pp}^{u u} = = {\overset{^^}{T T}}_{pp pp}^{u u} \cdot &Center Dot; {A A}^{' '} exp exp [[- - jω jω (({p p}_{11 p p} {x x}_{00} - - {q q}_{11 p p} {z z}_{00}))]],,$

${T T}_{ps ps}^{u u} = = {\overset{^^}{T T}}_{ps ps}^{u u} \cdot \cdot {A A}^{' '} exp exp [[- - jω jω (({p p}_{11 s the s} {x x}_{00} - - {q q}_{11 s the s} {z z}_{00}))]],,$

A′=exp[jω(p₂x₀-q₂z₀)]A′=exp[jω(p ₂ x ₀ -q ₂ z ₀ )]

$[\begin{matrix} {p p}_{22 p p} \\ {q q}_{22 p p} \end{matrix}] = = {C C}^{- - 11} [\begin{matrix} {\overset{^^}{p p}}_{22} \\ {\overset{^^}{q q}}_{22 p p} \end{matrix}],,$ $[\begin{matrix} {p p}_{22 s the s} \\ {q q}_{22 s the s} \end{matrix}] = = {C C}^{- - 11} [\begin{matrix} {\overset{^^}{p p}}_{22} \\ {\overset{^^}{q q}}_{22 s the s} \end{matrix}],,$ $[\begin{matrix} {p p}_{11 p p} \\ {q q}_{11 p p} \end{matrix}] = = C C [\begin{matrix} {\overset{^^}{p p}}_{11} \\ {\overset{^^}{q q}}_{11 p p} \end{matrix}],,$ $[\begin{matrix} {p p}_{11 s the s} \\ {q q}_{11 s the s} \end{matrix}] = = C C [\begin{matrix} {\overset{^^}{p p}}_{11} \\ {\overset{^^}{q q}}_{11 s the s} \end{matrix}],,$

${\hat{q}}_{ip} = \sqrt{α_{i}^{- 2} - {\hat{p}}_{1}^{2}},$ ${\hat{q}}_{is} = \sqrt{β_{i}^{- 2} - {\hat{p}}_{1}^{2}},$ (i=1，2) ${\hat{p}}_{1} = {\hat{p}}_{2},$ $[\begin{matrix} {\hat{p}}_{2} \\ {\hat{q}}_{2} \end{matrix}] = C^{- 1} [\begin{matrix} p_{2} \\ q_{2} \end{matrix}] .$ ${\hat{q}}_{ip} = \sqrt{α_{i}^{- 2} - {\hat{p}}_{1}^{2}},$ ${\hat{q}}_{is} = \sqrt{β_{i}^{- 2} - {\hat{p}}_{1}^{2}},$ (i=1, 2) ${\hat{p}}_{1} = {\hat{p}}_{2},$ $[\begin{matrix} {\hat{p}}_{2} \\ {\hat{q}}_{2} \end{matrix}] = C^{- 1} [\begin{matrix} p_{2} \\ q_{2} \end{matrix}] .$

其中

即分别为全局坐标系下倾斜界面的P-P反射系数、P-S反射系数、P-P透射系数和P-S透射系数。in

That is, PP reflection coefficient, PS reflection coefficient, PP transmission coefficient and PS transmission coefficient of the inclined interface in the global coordinate system, respectively.

同理，对于SV波从倾斜界面上方（或下方）入射的情况，其求解与P波入射时相同。Similarly, for the case where the SV wave is incident from above (or below) the inclined interface, the solution is the same as that of the P wave incident.

上述步骤3）中GFM积分方法为：The GFM integration method in the above step 3) is:

为了利用GFM积分法，先将步骤3）中积分式改写成如下形式(这里只给出求和号里的一项)：In order to use the GFM integral method, first rewrite the integral formula in step 3) into the following form (only one item in the summation is given here):

P=∫_Γf(p)exp[sg(p)]dpP= _∫Γ f(p)exp[sg(p)]dp

其中：in:

$f (p) = \frac{{- R}_{i} (p)}{2 π \cdot 2 jq},$ s=jω·max(x,z), $g (p) = \frac{p_{1 i} x - q_{1 i} z}{\max (x, z)} .$ $f (p) = \frac{{- R}_{i} (p)}{2 π &Center Dot; 2 jq},$ s=jω·max(x,z), $g (p) = \frac{p_{1 i} x - q_{1 i} z}{\max (x, z)} .$

GFM梯形法则为：The GFM trapezoidal rule is:

${&Integral; &Integral;}_{a a}^{b b} f f ((p p)) exp exp [[sg sg ((p p))]] dp dp = = \{\begin{matrix} \frac{δp δp}{sδg sδg} [[δ δ (({fe fe}^{sg sg})) - - \frac{δ δ ((f f)) δ δ (({e e}^{sg sg}))}{sδ sδ ((g g))}]] & δg δ g &NotEqual; &NotEqual; 00 \\ \frac{δp δp}{22} [[f f ((a a)) {e e}^{sg sg ((a a))} + + f f ((b b)) {e e}^{sg sg ((b b))}]] & δg δ g = = 00 \end{matrix} . .$

这里b-a等于积分步长，对任意函数h，δh=h(b)-h(a)。Here b-a is equal to the integral step size, for any function h, δh=h(b)-h(a).

Claims

1. A wave field forward modeling method in an inclined layered viscoelastic medium, is characterized in that, comprises the following steps:

1) Decompose the spherical wave excited by the point source into a plane wave;

2) Propagate the plane wave in the inclined layered viscoelastic medium, and when encountering the interface, proceed as follows: take the interface as

Vertically downward to the interface, convert the incident wave in the global coordinate system to the local coordinate system; obtain the reflected wave and transmitted wave in the local coordinate system; then convert the reflected wave and transmitted wave in the local coordinate system to the global coordinate system system, the reflection coefficient and transmission coefficient in the global coordinate system, as well as the reflected wave and the transmitted wave are obtained; on each interface, the type of wave is selected according to the needs, and the wave that can be selected when the P wave is incident includes reflected PP wave, reflected PS wave , transmitted PP wave, transmitted PS wave; the waves that can be selected when the S wave is incident include reflected SP wave, reflected SS wave, transmitted SP wave, and transmitted SS wave; thereby simulating the selected wave type and obtaining a simulated plane wave;

3) The plane wave simulated by step 2) synthesizes the wave field excited by the point source according to the following formula:

P P = = \frac{11}{22 π π} {&Integral; &Integral;}_{Γ Γ} {Σ Σ}_{i i = = 11}^{N N - - 11} [[{R R}_{i i} ((p p)) \cdot \cdot \frac{exp exp [[jω jω (({p p}_{11 i i} x x - - {q q}_{11 i i} z z))]]}{- - 22 jq jq}]] dp dp

In the formula, Γ is the integral path; N is the number of plane waves; R _i (p) is the synthesis of transmission and reflection coefficients of plane waves on multiple interfaces; j is the imaginary unit, ω is the angular frequency; p _1i and q _1i are respectively is the horizontal and vertical slowness of the plane wave reflected from the i-th interface back to the first layer; p and q are the horizontal and vertical slowness of the plane wave at the point source, respectively, and x and z represent the position coordinates.

2. The wavefield forward modeling method in inclined layered viscoelastic media according to claim 1, characterized in that, in step 2), the simulated seismic source is a P wave when the primary reflected P wave, and the seismic source and the geophone are both In the first layer, the recursive steps of propagating the incident wave at the source down layer by layer are:

①When encountering an interface i, take this interface as

Axis, the intersection of the interface and the Z axis of the global coordinate system is the origin to establish a local coordinate system Downward perpendicular to the interface; transform the incident wave φ _i in the global coordinate system into the local coordinate system to get

② Obtain the reflection coefficient in the local coordinate system

and transmission coefficient

and reflected waves

and transmitted waves

③ Then convert the reflected wave and transmitted wave in the local coordinate system to the global coordinate system to obtain the reflection coefficient in the global coordinate system

and transmission coefficient And reflected wave φ′ _{i, i} and transmitted wave φ _i+1 ;

④Reflected wave φ′ _{i, i} is recursively pushed up to the first layer layer by layer, and steps ①-③ are performed when encountering an interface, but only the transmitted wave is taken;

⑤ The transmitted wave φ _i+1 continues to deduce downwards, and returns to step ① when the next interface is encountered.

3. The wavefield forward modeling method in inclined layered viscoelastic media according to claim 1, characterized in that, in step 2), simulating other types of waves adopts the primary reflection P wave when the simulated source is P wave The same method and steps are carried out.

4. The wavefield forward modeling method in inclined layered viscoelastic media according to claim 1, characterized in that, in step 3), the integration method in the formula adopts the GFM integration method proposed by Frazer and Gettrust.