CN102096100A

Movatterモバイル変換

Info

Publication number: CN102096100A
Application number: CN 201010559213
Authority: CN
Inventors: 王尚旭
Original assignee: China University of Petroleum Beijing; China National Petroleum Corp
Current assignee: China University of Petroleum Beijing; China National Petroleum Corp
Priority date: 2010-11-24
Filing date: 2010-11-24
Publication date: 2011-06-15

Abstract

The embodiment of the invention provides a logging curve and seismographic record full-wave matching method and device, wherein the method comprises the following steps of: acquiring seismographic records; processing the seismographic records by utilizing a reflectivity method to obtain well bypass seismographic records; and matching the well bypass seismographic records with a logging curve. With the embodiment of the invention, multi-angle reflective seismic waves, interlayer repeated waves and interface converting waves can be simultaneously taken into account so that the matching precision of well and ground data is improved.

Description

Translated fromChinese

测井曲线与地震记录全波匹配方法及装置Method and device for full-wave matching of logging curves and seismic records

技术领域technical field

本发明涉及石油地球物理勘探地震模型正演模拟领域，特别涉及一种测井曲线与地震记录全波匹配方法及装置。The invention relates to the field of forward simulation of seismic models for petroleum geophysical exploration, in particular to a method and device for full-wave matching of logging curves and seismic records.

背景技术Background technique

在油气勘探开发过程中，地球物理测井数据常用于约束地震数据进行油气储层预测，地球物理测井数据与地震数据匹配的技术(以下简称井地匹配技术)精度直接影响油气储层预测精度。In the process of oil and gas exploration and development, geophysical logging data is often used to constrain seismic data to predict oil and gas reservoirs. The accuracy of geophysical logging data and seismic data matching technology (hereinafter referred to as well-ground matching technology) directly affects the accuracy of oil and gas reservoir prediction .

现有井地匹配技术属于褶积模型匹配技术，即将声波测井数据转换成声波速度，结合密度测井数据生成反射系数序列，反射系数序列与子波的褶积形成井旁地震记录。The existing well-ground matching technology belongs to the convolution model matching technology, which converts the acoustic logging data into acoustic velocity, combines the density logging data to generate a reflection coefficient sequence, and convolutes the reflection coefficient sequence and wavelet to form a side-hole seismic record.

这种由褶积方法形成的地震记录仅仅是地层一次反射波的结果。主要存在以下问题：(1)没有考虑地震波传播方向并不垂直于地层，事实上，得到的叠后地震记录是由多角度反射地震波叠加而成的；(2)没有考虑层间多次反射波的存在，散射衰减效应被忽略；(3)没有考虑转换波的存在，即地震波在介质分界面上存在纵波与横波能量的相互转化。This seismic record formed by the convolution method is only the result of a reflection wave from the formation. The main problems are as follows: (1) It does not take into account that the direction of seismic wave propagation is not perpendicular to the formation. In fact, the obtained post-stack seismic record is formed by superposition of multi-angle reflected seismic waves; (2) It does not take into account the multiple reflection waves between layers. The existence of the scattering attenuation effect is ignored; (3) The existence of the converted wave is not considered, that is, the mutual transformation of the energy of the longitudinal wave and the shear wave exists in the seismic wave at the interface of the medium.

因此在实现本发明的过程中，发明人发现现有技术的缺陷在于：技术精度不高，即合成地震记录与实际井旁地震记录相似度较低(一般低于90％)。Therefore, in the process of realizing the present invention, the inventor found that the defect of the prior art lies in: the technical precision is not high, that is, the similarity between the synthetic seismic record and the actual side-hole seismic record is low (generally less than 90%).

发明内容Contents of the invention

本发明实施例提供一种测井曲线与地震记录全波匹配方法及装置，目的在于提升井地数据匹配精度。Embodiments of the present invention provide a method and device for full-wave matching of logging curves and seismic records, with the purpose of improving the matching accuracy of well-ground data.

为达到上述目的，本发明实施例提供一种测井曲线与地震记录全波匹配的方法，该方法包括：获取地震记录；利用反射率法对地震记录进行处理，获得井旁道地震记录；将井旁道地震记录与测井曲线进行匹配。In order to achieve the above object, the embodiment of the present invention provides a method for full-wave matching of logging curves and seismic records, the method includes: acquiring seismic records; processing the seismic records by using the reflectivity method to obtain well side channel seismic records; The well bypass seismic records are matched with well logs.

本发明实施例还提供一种测井曲线与地震记录全波匹配的装置，该装置包括：An embodiment of the present invention also provides a device for full-wave matching of logging curves and seismic records, the device comprising:

获取单元，获取地震记录；Get the unit, get the seismic record;

处理单元，利用反射率法对地震记录进行处理，获得井旁道地震记录；The processing unit uses the reflectivity method to process the seismic records to obtain the well bypass seismic records;

匹配单元，将井旁道地震记录与测井曲线进行匹配。The matching unit matches the well bypass seismic records with the logging curves.

本发明实施例的有益效果在于，通过在进行测井数据与地震数据匹配时，同时考虑多角度反射地震波、层间多次波及界面转换波，可以提升井地数据匹配精度。The beneficial effect of the embodiment of the present invention is that, by considering multi-angle reflected seismic waves, interlayer multiple waves and interface converted waves when matching well logging data and seismic data, the matching accuracy of well ground data can be improved.

附图说明Description of drawings

此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，并不构成对本发明的限定。在附图中：The drawings described here are used to provide further understanding of the present invention, constitute a part of the application, and do not limit the present invention. In the attached picture:

图1是本发明实施例的测井曲线与地震记录全波匹配方法的流程图；Fig. 1 is the flow chart of the logging curve and seismic record full-wave matching method of the embodiment of the present invention;

图2是本发明实施例的层间多次波的示意图；Fig. 2 is a schematic diagram of an interlayer multiple wave in an embodiment of the present invention;

图3是本发明实施例的任一区域的反射和透射情况的示意图；Fig. 3 is a schematic diagram of reflection and transmission of any region in an embodiment of the present invention;

图4是本发明实施例的单层区域的含多次波的演示示意图；Fig. 4 is a demonstration schematic diagram of multiple waves in a single-layer region according to an embodiment of the present invention;

图5是本发明实施例的双层区域的传播演示示意图；Fig. 5 is a schematic diagram of a propagation demonstration of a double-layer region according to an embodiment of the present invention;

图6是本发明实施例的地震波场的合成图；Fig. 6 is the synthetic figure of the seismic wave field of the embodiment of the present invention;

图7是本发明实施例的频率汉明窗和慢度汉明窗示意图；Fig. 7 is a schematic diagram of a frequency Hamming window and a slowness Hamming window according to an embodiment of the present invention;

图8是本发明实施例的Ricker子波时间域波形示意图；FIG. 8 is a schematic diagram of a Ricker wavelet time-domain waveform according to an embodiment of the present invention;

图9是本发明实施例的Ricker子波频率域波形示意图；9 is a schematic diagram of a Ricker wavelet frequency domain waveform in an embodiment of the present invention;

图10是本发明实施例的由模拟测井数据提取的反射系数、褶积记录、反射率法正演结果的示意图；Fig. 10 is a schematic diagram of the reflection coefficient, convolution record and reflectivity method forward modeling results extracted from the simulated logging data according to the embodiment of the present invention;

图11是本发明实施例的三层模型一次波、多次波和转换波记录的示意图；Fig. 11 is a schematic diagram of the three-layer model primary wave, multiple wave and converted wave records of the embodiment of the present invention;

图12是本发明实施例的np5-4井速度和密度测井资料及井旁道地震记录的示意图；Fig. 12 is a schematic diagram of well np5-4 velocity and density logging data and side channel seismic records of the embodiment of the present invention;

图13是本发明实施例的对np5-4井测井数据分层及用反射率法模拟记录的示意图；Fig. 13 is a schematic diagram of layering logging data of well np5-4 and simulating records by reflectivity method according to an embodiment of the present invention;

图14是本发明实施例的测井曲线与地震记录全波匹配装置的构成图；Fig. 14 is a composition diagram of a full-wave matching device for logging curves and seismic records according to an embodiment of the present invention;

图15是本发明实施例的处理单元的构成图。Fig. 15 is a configuration diagram of a processing unit according to an embodiment of the present invention.

具体实施方式Detailed ways

为使本发明的目的、技术方案和优点更加清楚明白，下面结合附图对本发明实施例作进一步详细说明。在此，本发明的示意性实施例及其说明用于解释本发明，但并不作为对本发明的限定。In order to make the object, technical solution and advantages of the present invention clearer, the embodiments of the present invention will be further described in detail below in conjunction with the accompanying drawings. Here, the exemplary embodiments and descriptions of the present invention are used to explain the present invention, but not to limit the present invention.

本发明实施例提供一种测井曲线与地震记录全波匹配的方法，如图1所示，所述方法包括：Embodiments of the present invention provide a method for full-wave matching of logging curves and seismic records, as shown in Figure 1, the method includes:

步骤101，获取地震记录；Step 101, obtaining seismic records;

步骤102，利用反射率法对地震记录进行处理，获得井旁道地震记录；Step 102, using the reflectivity method to process the seismic records to obtain the well bypass seismic records;

步骤103，将井旁道地震记录与测井曲线进行匹配。Step 103, matching the well side channel seismic records with the logging curves.

在地震勘探资料中，除了一次反射波、随机干扰以外，还存在着多次波和转换波。层间多次波又称为内部多次波，是指所有下行反射发生在除自由界面以外的其它反射界面的多次波。In seismic exploration data, in addition to the primary reflection wave and random interference, there are also multiple waves and converted waves. Interlayer multiples, also known as internal multiples, refer to the multiples in which all downgoing reflections occur at other reflecting interfaces except the free interface.

图2是层间多次波的示意图，如图2所示，内部多次波的反射界面通常很难精确确定，因此内部多次波很难预测，也较难进行有效地消除。此外内部多次波具有与有效波相似甚至更高的速度，这使得内部多次波的压制更为困难。当假设条件不能满足时，多次波就不但不能被消除，压制措施甚至可能损伤有效波。Figure 2 is a schematic diagram of interlayer multiple waves. As shown in Figure 2, the reflection interface of internal multiple waves is usually difficult to accurately determine, so internal multiple waves are difficult to predict and effectively eliminate. In addition, the internal multiples have similar or even higher velocities than the effective waves, which makes the suppression of internal multiples more difficult. When the assumed conditions cannot be satisfied, the multiple wave cannot be eliminated, and the suppression measures may even damage the effective wave.

在实际资料的处理和解释中，多次波的类型一般根据对时间、速度的目测和经验来确定。但在有些多次波与有效波没有明显的差异，常常难以确定。层间多次波同相轴通常以不同的正常时差迭加在一次波同相轴之上，和一次有效反射波相干涉，就往往给处理和解释带来很多麻烦。如果浅、中层存在良好的反射界面并产生多次波，就有可能掩盖了中、深层的一次反射波。在剖面上多次波较强时，如果在解释中不能正确地把多次波识别出来，就会造成错误的地质解释。In the processing and interpretation of actual data, the type of multiple waves is generally determined by visual observation and experience of time and velocity. But there is no obvious difference between some multiple waves and effective waves, so it is often difficult to determine. The interlayer multiple event is usually superimposed on the primary event with different normal moveouts, and interferes with an effective reflected wave, which often brings a lot of trouble to processing and interpretation. If there is a good reflection interface in the shallow and middle layers and multiple waves are generated, it may cover the primary reflection waves in the middle and deep layers. When the multiple waves are strong on the section, if the multiple waves cannot be correctly identified in the interpretation, it will cause wrong geological interpretation.

在本实施例中，采用反射率法正演模拟地震记录，然后做试验分析。反射率法是应用最广泛的叠前地震记录计算方法，它能够提供波场的完全解，包括地震波在地层内传播时形成的一次波、多次波和转换波。In this embodiment, the reflectivity method is used to simulate the seismic records forward, and then the test analysis is done. The reflectivity method is the most widely used calculation method for prestack seismic records, which can provide a complete solution of the wave field, including the primary wave, multiple wave and converted wave formed when the seismic wave propagates in the formation.

以下首先对反射率法的原理进行简单的介绍。Firstly, the principle of the reflectance method is briefly introduced as follows.

基本的反射率算法模拟的层状介质，其上边界是半空间或自由表面，下边界为半空间。The layered medium simulated by the basic albedo algorithm has an upper boundary of a half-space or a free surface and a lower boundary of a half-space.

在本实施例中，讨论局限于圆柱坐标系统，x轴(径向)和z轴(纵向)。空间为均匀介质，而且设x为射线参数，z为深度参数。激发和接收点可以在模型的任意处。In this example, the discussion is limited to the cylindrical coordinate system, x-axis (radial) and z-axis (longitudinal). The space is a homogeneous medium, and let x be the ray parameter and z be the depth parameter. The excitation and reception points can be anywhere on the model.

反射率方法和其它的计算合成记录方法一样，是一种对弹性动力学方程The reflectivity method, like other computational synthetic recording methods, is a response to the elastodynamic equation

$- - {ρω ρω}^{22} u u = = &dtri; &dtri; ξ ξ + + f f - - - - - - ((1.1 1.1))$

和关系式and relational

求积分的方法。method of scoring.

这里的ξ和ζ是无限小的压力和张力。ρ是密度，u是弹性模量。而f是体力密度。

是个四阶的弹性张量，在柱空间的均匀介质最通常的形式看作是一个关于垂直轴线对称的横向同性体。

因此有五个独立的参数。Here ξ and ζ are infinitesimal pressure and tension. ρ is the density and u is the modulus of elasticity. And f is physical density.

is a fourth-order elastic tensor, and the most common form of a homogeneous medium in cylindrical space is regarded as a transverse isotropic body symmetric about the vertical axis.

So there are five independent parameters.

一个经常被忽略的方面是式(1.2)根本上是个频率域的关系式，既然

的组成通常是复杂的频率依靠的量。这样假设以时间为变量的一个暂时的傅立叶变换：An often overlooked aspect is that Equation (1.2) is fundamentally a frequency-domain relation, since

The composition of is usually a complex frequency-dependent quantity. This assumes a temporal Fourier transform as a function of time:

F＝∫dte^iωt (1.3)F＝∫dte^iωt (1.3)

既然模型是水平层状的横向各向同性的介质，可以通过傅立叶-汉克尔变换将其扩展成一系列的圆柱状，而使变量变换成半径和方位角。这一过程的详细内容被许多的报纸和书介绍。结论性的等式是：Since the model is a horizontally layered transversely isotropic medium, it can be expanded into a series of cylinders by Fourier-Hankel transform, and the variables transformed into radii and azimuths. The details of this process are described in many newspapers and books. The conclusive equation is:

上等式中，z是深度，ω是频率，是系统矩阵，对于P-SV系统是在半径和垂直的方向，而且是个4*4的矩阵，其元素与地层参数和水平慢度有关；b矢量包含了位移和应力的变换。这里考虑n＝0时的傅立叶-汉克尔逆变换得到：In the above equation, z is the depth, ω is the frequency, is the system matrix, for the P-SV system, it is in the radial and vertical directions, and it is a 4*4 matrix, whose elements are related to formation parameters and horizontal slowness; the b vector contains the transformation of displacement and stress. Considering here the Fourier-Hankel inverse transform when n=0:

$u u ((x x,, ω ω)) = = {&Integral; &Integral;}_{00}^{\infty \infty} {ω ω}^{22} pdp pdp \overset{~ ~}{u u} ((ω ω,, p p)) {J J}_{00} ((ωpx ωpx)) - - - - - - ((1.5 1.5))$

在等式(1.5)中，p是水平慢度，ω是弧度暂时频率，x是炮检距。当积分式(1.5)计算了一系列的频率后，通过傅立叶逆变换：In equation (1.5), p is the horizontal slowness, ω is the radian temporal frequency, and x is the offset. When the integral formula (1.5) calculates a series of frequencies, through the Fourier inverse transform:

$u u ((x x,, t t)) = = \frac{11}{22 π π} {&Integral; &Integral;}_{- - \infty \infty}^{\infty \infty} {dωe dωe}^{- - iωt iωt} u u ((x x,, ω ω)) - - - - - - ((1.6 1.6))$

变换到时域。Transform to the time domain.

这里的J₀是零阶的BESSEL函数，

是频率慢度域的反射透射系数响应。这样反射率法实现的主要工作就包括两部分，一部分是在频率慢度域求解反射透射系数响应

另一部分是积分实现，求解出时间空间域的波场。Here J₀ is the zero-order BESSEL function,

is the reflection-transmission coefficient response in the frequency-slowness domain. In this way, the main work of the reflectivity method includes two parts, one is to solve the reflection and transmission coefficient response in the frequency slowness domain

The other part is integral implementation, which solves the wave field in the time-space domain.

因此，在步骤102实施时，具体包括：在频率慢度域求解反射透射系数响应；根据反射透射系数响应进行慢度积分和频率积分，得到井旁地震记录；Therefore, whenstep 102 is implemented, it specifically includes: solving the reflection and transmission coefficient response in the frequency slowness domain; performing slowness integration and frequency integration according to the reflection and transmission coefficient response to obtain side-hole seismic records;

具体地，在频率慢度域求解反射透射系数响应，包括：分别求取自由界面、固体-固体界面和任意区域的反射透射系数；将反射透射系数合成为总的反射透射系数响应。Specifically, solving the reflection and transmission coefficient response in the frequency slowness domain includes: separately calculating the reflection and transmission coefficient of the free interface, solid-solid interface and any region; synthesizing the reflection and transmission coefficient into the total reflection and transmission coefficient response.

在本实施例中，基于波场以及慢度理论分别求取自由界面、固体-固体界面和一个区域的反射和透射系数，最后合成总的反射透射系数响应。In this embodiment, the reflection and transmission coefficients of a free interface, a solid-solid interface and a region are respectively obtained based on the wave field and slowness theory, and finally the total reflection and transmission coefficient response is synthesized.

对于自由界面：For free interface:

地球表面是一个特殊的分界面，它将无限介质划分为两个半空间，地面以上是空气介质，其密度与地面以下的岩石或海平面以下的海水层相比可以忽略。地球表面可以看成是一个弹性半空间表面，表面以下视为理想弹性介质，这种界面称为自由界面。自由界面上的应力作用为零。The Earth's surface is a special interface that divides an infinite medium into two half-spaces, and above the ground is an air medium whose density is negligible compared to the rock below the ground or the seawater layer below sea level. The surface of the earth can be regarded as an elastic half-space surface, and the surface below the surface is regarded as an ideal elastic medium. This interface is called a free interface. The stress action on the free interface is zero.

在本实施例中，研究一个平面波入射到自由表面时的反射及透射问题。取直角坐标系，x₃＝0为弹性半无限空间的自由表面。这样，就可以利用慢度p建立一个表达式：In this example, the problem of reflection and transmission of a plane wave incident on a free surface is studied. Taking a rectangular coordinate system, x₃ =0 is the free surface of the elastic semi-infinite space. In this way, an expression can be built using the slowness p:

exp(iω[px₁-t]) (2.1.1)exp(iω[px₁ -t]) (2.1.1)

先以简单的SH波为例。考虑上行波H_U入射到界面上产生的SH反射波。这将会产生下行SH波H_D，这里由于界面是水平的，所以没有产生转换波。由于具有相同的慢度，下行波的反射角等于上行波的入射角。Take the simple SH wave as an example first. Consider the SH reflection wave produced by the upgoing wave H_U incident on the interface. This will generate a downgoing SH wave_HD , where no converted wave is generated because the interface is horizontal. Due to the same slowness, the angle of reflection of the downgoing wave is equal to the angle of incidence of the upgoing wave.

利用位移-应力关系可以得到表达式Using the displacement-stress relationship, the expression

$(\begin{matrix} W W \\ T T \end{matrix}) = = {ϵ ϵ}_{H h} (\begin{matrix} 11 & 11 \\ {- - iμ iμ}_{00} {q q}_{β β 00} & {iμ iμ}_{00} {q q}_{β β 00} \end{matrix}) (\begin{matrix} {H h}_{U u} & exp exp [[{- - iωq iωq}_{β β} {x x}_{33}]] \\ {H h}_{D D.} & exp exp [[{+ + iωq iωq}_{β β} {x x}_{33}]] \end{matrix}) - - - - - - ((22 . . 1.2 1.2))$

其中in

ε_H＝1/(2μq_β)^1/2 (2.1.3)ε_H ＝1/(2μq_β )^1/2 (2.1.3)

当处于自由表面时，由于在x₃＝0时应力为零，所以上述表达式就变为When on a free surface, since the stress is zero at x₃ =0, the above expression becomes

$(\begin{matrix} W W ((00)) \\ 00 \end{matrix}) = = {ϵ ϵ}_{H h} (\begin{matrix} 11 & 11 \\ {- - iμ iμ}_{00} {q q}_{β β 00} & {iμ iμ}_{00} {q q}_{β β 00} \end{matrix}) (\begin{matrix} {H h}_{U u} \\ {H h}_{D D.} \end{matrix}) - - - - - - ((2.1.4 2.1.4))$

由应力方程by the stress equation

iμ₀q_β0(H_U-H_D)＝0 (2.1.5)iμ₀ q_β0 (_HU -H_D )＝0 (2.1.5)

得到get

H_U＝H_D (2.1.6)H_U =H_D (2.1.6)

因此反射系数Therefore the reflection coefficient

${R R}_{F f}^{HH HH} = = 11 - - - - - - ((2.1.7 2.1.7))$

对于P-SV波，情况比较复杂。慢度相同意味着反射前后没有发生波型的转换，反射角与入射角相同。当发生波型转换时，就需要用到Snell定律，在横向各向同性介质中P波入射角i和S波反射角j具有以下关系For P-SV waves, the situation is more complicated. The same slowness means that there is no mode conversion before and after reflection, and the angle of reflection is the same as the angle of incidence. When the mode conversion occurs, Snell's law needs to be used. In the transversely isotropic medium, the P-wave incident angle i and the S-wave reflection angle j have the following relationship

$\frac{sin sin i i}{{α α}_{00}} = = \frac{sin sin j j}{{β β}_{00}} - - - - - - ((2.1.8 2.1.8))$

其中α₀为P波速度，β₀为S波速度。Where α₀ is the P-wave velocity, and β₀ is the S-wave velocity.

在本实施例中，可以用上行波和下行波得到波场的表达式：In this embodiment, the expression of the wave field can be obtained by using the up-going wave and the down-going wave:

$b b (({x x}_{33})) = = (\begin{matrix} w w (({x x}_{33})) \\ t t (({x x}_{33})) \end{matrix}) = = (\begin{matrix} {m m}_{U u} & {m m}_{D D.} \\ {n no}_{U u} & {n no}_{D D.} \end{matrix}) (\begin{matrix} {v v}_{U u} (({x x}_{33})) \\ {v v}_{D D.} (({x x}_{33})) \end{matrix}) - - - - - - ((2.1.9 2.1.9))$

其中，m_U、n_U与P、S波的位移和应力一致。各向同性介质中的位移和应力矩阵具有下面的形式：Among them, m_U , n_U are consistent with the displacement and stress of P and S waves. The displacement and stress matrices in an isotropic medium have the form:

${m m}_{U u,, D D.} = = (\begin{matrix} \overset{&OverBar; &OverBar;}{+ +} {iq iq}_{α α} {ϵ ϵ}_{α α} & {pϵ pϵ}_{β β} \\ {pϵ pϵ}_{α α} & {\overset{&OverBar; &OverBar;}{+ +} iq iq}_{α α} {ϵ ϵ}_{α α} \end{matrix})$

${n no}_{U u,, D D.} = = (\begin{matrix} ρ ρ (({22 β β}^{22} {p p}^{22} - - 11)) {ϵ ϵ}_{α α} & \overset{&OverBar; &OverBar;}{+ +} {22 iρβ iρβ}^{22} {pq pq}_{β β} {ϵ ϵ}_{β β} \\ \overset{&OverBar; &OverBar;}{+ +} {22 iρβ iρβ}^{22} {pq pq}_{α α} {ϵ ϵ}_{α α} & ρ ρ (({22 β β}^{22} {p p}^{22} - - 11)) {ϵ ϵ}_{β β} \end{matrix}) - - - - - - ((2.1.10 2.1.10))$

为了满足应力t在自由界面为零的条件，上行波场v_U将产生下行波场v_D。条件t(0)＝0要求下行波场v_D满足In order to satisfy the condition that the stress t is zero at the free interface, the upgoing wave field v_U will generate the downgoing wave field v_D . The condition t(0)=0 requires the downgoing wave field v_D to satisfy

$(\begin{matrix} w w ((00)) \\ 00 \end{matrix}) = = (\begin{matrix} {m m}_{U u} & {m m}_{D D.} \\ {n no}_{U u} & {n no}_{D D.} \end{matrix}) (\begin{matrix} {v v}_{U u} \\ {v v}_{D D.} \end{matrix}) - - - - - - ((2.1.11 2.1.11))$

因此therefore

n_Uv_U+n_Dv_D＝0 (2.1.12)n_U v_U + n_D v_D = 0 (2.1.12)

引入反射系数矩阵R_F，它包含不同类型的波的反射系数，即Introduce the reflection coefficient matrix R_F , which contains the reflection coefficients for different types of waves, namely

${R R}_{F f} = = (\begin{matrix} {R R}_{F f}^{PP PP} & {R R}_{F f}^{PS P.S.} \\ {R R}_{F f}^{SP SP} & {R R}_{F f}^{SS SS} \end{matrix}) - - - - - - ((2.1.13 2.1.13))$

这样so

v_D＝R_Fv_U (2.1.14)v_D =R_F v_U (2.1.14)

就可以扩展为can be expanded to

$(\begin{matrix} {P P}_{D D.} \\ {S S}_{D D.} \end{matrix}) = = (\begin{matrix} {R R}_{F f}^{PP PP} & {R R}_{F f}^{PS P.S.} \\ {R R}_{F f}^{SP SP} & {R R}_{F f}^{SS SS} \end{matrix}) (\begin{matrix} {P P}_{U u} \\ {S S}_{U u} \end{matrix}) - - - - - - ((2.1.15 2.1.15))$

从(2.1.12)和(2.1.13)可以定义反射系数矩阵为From (2.1.12) and (2.1.13), the reflection coefficient matrix can be defined as

${R R}_{F f} = = - - {n no}_{D D.}^{- - 11} {n no}_{U u} - - - - - - ((2.1.16 2.1.16))$

由此可以得到From this you can get

$(\begin{matrix} {R R}_{F f}^{PP PP} & {R R}_{F f}^{PS P.S.} \\ {R R}_{F f}^{SP SP} & {R R}_{F f}^{SS SS} \end{matrix}) = = \frac{11}{{44 p p}^{22} {q q}_{α α 00} {q q}_{β β 00} + + {v v}^{22}} (\begin{matrix} {44 p p}^{22} {q q}_{α α 00} {q q}_{β β 00} - - {v v}^{22} & 44 ipv ipv {(({q q}_{α α 00} {q q}_{β β 00}))}^{11 / / 22} \\ 44 ipv ipv {(({q q}_{α α 00} {q q}_{β β 00}))}^{11 / / 22} & 44 {p p}^{22} {q q}_{α α 00} {q q}_{β β 00} - - {v v}^{22} \end{matrix}) - - - - - - ((2.1.17 2.1.17))$

其中in

$v v = = (({22 p p}^{22} - - {β β}_{00}^{- - 22})) - - - - - - ((2.1.18 2.1.18))$

对于固体-固体表面：For solid-solid surfaces:

在不同介质的分界面，解决问题时需要两点原则：第一，位移连续；第二，应力连续。对于水平分界面，还需要保持穿过界面前后的水平相位保持一致。因此相互作用的平面波必须具有相同的水平慢度p。At the interface of different media, two principles are needed to solve the problem: first, the displacement is continuous; second, the stress is continuous. For a horizontal interface, it is also necessary to keep the horizontal phase consistent before and after passing through the interface. Therefore the interacting plane waves must have the same horizontal slowness p.

首先，讨论简单的SH波情况。在界面上部的介质1中，用上行波及下行波表示的位移及应力为First, the simple SH wave case is discussed. In the medium 1 at the upper part of the interface, the displacement and stress represented by up-going wave and down-going wave are

$(\begin{matrix} {u u}_{22} \\ {ω ω}^{- - 11} {τ τ}_{23 twenty three} \end{matrix}) = = {ϵ ϵ}_{H h 11} (\begin{matrix} 11 & 11 \\ {- - iμ iμ}_{11} {q q}_{β β 11} & {iμ iμ}_{11} {q q}_{β β 11} \end{matrix}) (\begin{matrix} {H h}_{U u 11} \\ {H h}_{D D. 11} \end{matrix}) - - - - - - ((2.2.1 2.2.1))$

类似的，在界面下部的介质2中，位移及应力的表达式为Similarly, in the medium 2 at the lower part of the interface, the expressions of displacement and stress are

$(\begin{matrix} {u u}_{22} \\ {ω ω}^{- - 11} {τ τ}_{23 twenty three} \end{matrix}) = = {ϵ ϵ}_{H h 22} (\begin{matrix} 11 & 11 \\ {- - iμ iμ}_{22} {q q}_{β β 22} & {iμ iμ}_{22} {q q}_{β β 22} \end{matrix}) (\begin{matrix} {H h}_{U u 22} \\ {H h}_{D D. 22} \end{matrix}) - - - - - - ((2.2 2.2 . . 22))$

由于界面处位移及应力的连续性，所以有Due to the continuity of displacement and stress at the interface, there is

${ϵ ϵ}_{H h 11} (\begin{matrix} 11 & 11 \\ {- - iμ iμ}_{11} {q q}_{β β 11} & {iμ iμ}_{11} {q q}_{β β 11} \end{matrix}) (\begin{matrix} {H h}_{U u 11} \\ {H h}_{D D. 11} \end{matrix}) = = {ϵ ϵ}_{H h 22} (\begin{matrix} 11 & 11 \\ {- - iμ iμ}_{22} {q q}_{β β 22} & {iμ iμ}_{22} {q q}_{β β 22} \end{matrix}) (\begin{matrix} {H h}_{U u 22} \\ {H h}_{D D. 22} \end{matrix}) - - - - - - ((2.2.3 2.2.3))$

所以界面两侧的上、下行波可以表示为So the up-going and down-going waves on both sides of the interface can be expressed as

$(\begin{matrix} {H h}_{U u 11} \\ {H h}_{D D. 11} \end{matrix}) = = \frac{11}{22 {(({μ μ}_{11} {μ μ}_{22} {q q}_{β β 11} {q q}_{β β 22}))}^{11 / / 22}} (\begin{matrix} {μ μ}_{11} {q q}_{β β 11} + + {μ μ}_{22} {q q}_{β β 22} & {μ μ}_{11} {q q}_{β β 11} - - {μ μ}_{22} {q q}_{β β 22} \\ {μ μ}_{11} {q q}_{β β 11} - - {μ μ}_{22} {q q}_{β β 22} & {μ μ}_{11} {q q}_{β β 11} + + {μ μ}_{22} {q q}_{β β 22} \end{matrix}) (\begin{matrix} {H h}_{U u 22} \\ {H h}_{D D. 22} \end{matrix}) - - - - - - ((2.2.4 2.2.4))$

考虑一下行波从介质1入射到分界面上，它将在介质1中产生上行反射波，在介质2中产生下行透射波，由于没有在介质2中产生上行波，所以H_U2＝0。引入下行波反射系数

及透射系数

则有Considering that a traveling wave is incident on the interface from medium 1, it will generate upgoing reflected waves in medium 1 and downgoing transmitted waves inmedium 2. Since no upgoing waves are generated inmedium 2, H_U2 =0. Introducing the downgoing wave reflection coefficient

and transmission coefficient

then there is

${H h}_{U u 11} = = {R R}_{D D.}^{I I} {H h}_{D D. 11} - - - - - - ((2.2.5 2.2.5))$

${H h}_{D D. 22} = = {T T}_{D D.}^{I I} {H h}_{D D. 11} - - - - - - ((2.2 2.2 . . 66))$

式(2.2.4)可以表示为Formula (2.2.4) can be expressed as

$(\begin{matrix} {H h}_{U u 11} \\ {H h}_{D D. 11} \end{matrix}) = = (\begin{matrix} {Q Q}_{UU UU} & {Q Q}_{UD UD} \\ {Q Q}_{DU DU} & {Q Q}_{DD DD} \end{matrix}) (\begin{matrix} 00 \\ {H h}_{D D. 22} \end{matrix}) - - - - - - ((2.2.7 2.2.7))$

因此therefore

${R R}_{D D.}^{I I} = = {Q Q}_{UD UD} {(({Q Q}_{DD DD}))}^{- - 11}$

${T T}_{D D.}^{I I} = = {(({Q Q}_{DD DD}))}^{- - 11} - - - - - - ((2.2.8 2.2.8))$

所以SH波的反射及透射系数为So the reflection and transmission coefficients of SH wave are

${R R}_{D D.}^{I I} = = \frac{{μ μ}_{11} {q q}_{β β 11} - - {μ μ}_{22} {q q}_{β β 22}}{{μ μ}_{11} {q q}_{β β 11} + + {μ μ}_{22} {q q}_{β β 22}}$

${T T}_{D D.}^{I I} = = \frac{22 {(({μ μ}_{11} {μ μ}_{22} {q q}_{β β 11} {q q}_{β β 22}))}^{11 / / 22}}{{μ μ}_{11} {q q}_{β β 11} + + {μ μ}_{22} {q q}_{β β 22}} - - - - - - ((2.2.9 2.2.9))$

其次，对于P-SV波。在界面上部的介质1(x₃＝z_I-)中，用上行波及下行波表示的位移及应力为Second, for the P-SV wave. In the medium 1 (x₃ =z_I -) at the upper part of the interface, the displacement and stress represented by the upgoing wave and downgoing wave are

$b b (({z z}_{I I} - -)) = = (\begin{matrix} w w (({z z}_{I I} - -)) \\ t t (({z z}_{I I} - -)) \end{matrix}) = = (\begin{matrix} {m m}_{U u 11} & {m m}_{D D. 11} \\ {n no}_{U u 11} & {n no}_{D D. 11} \end{matrix}) (\begin{matrix} {v v}_{U u 11} (({z z}_{I I} - -)) \\ {v v}_{D D. 11} (({z z}_{I I} - -)) \end{matrix}) = = {D D.}_{11} {v v}_{11} - - - - - - ((2.2.10 2.2.10))$

$b b (({z z}_{I I} + +)) = = (\begin{matrix} w w (({z z}_{I I} + +)) \\ t t (({z z}_{I I} + +)) \end{matrix}) = = (\begin{matrix} {m m}_{U u 22} & {m m}_{D D. 22} \\ {n no}_{U u 22} & {n no}_{D D. 22} \end{matrix}) (\begin{matrix} {v v}_{U u 22} (({z z}_{I I} + +)) \\ {v v}_{D D. 22} (({z z}_{I I} + +)) \end{matrix}) = = {D D.}_{22} {v v}_{22} - - - - - - ((2.2.11 2.2.11))$

D₁v₁＝D₂v₂ (2.2.12)D₁ v₁ = D₂ v₂ (2.2.12)

即Right now

${v v}_{11} = = {D D.}_{11}^{- - 11} {D D.}_{22} {v v}_{22} = = {Qv Q}_{22} - - - - - - ((2.2.13 2.2.13))$

将v₁及v₂扩展成上、下行波的形式。类似的，Q扩展为2×2的矩阵形式Expand v₁ and v₂ into the form of up and down waves. Similarly, Q is expanded into a 2×2 matrix form

$(\begin{matrix} {v v}_{U u 11} \\ {v v}_{D D. 11} \end{matrix}) = = (\begin{matrix} {Q Q}_{UU UU} & {Q Q}_{UD UD} \\ {Q Q}_{DU DU} & {Q Q}_{DD DD} \end{matrix}) (\begin{matrix} {v v}_{U u 22} \\ {v v}_{D D. 22} \end{matrix}) = = i i (\begin{matrix} {- - n no}_{D D. 11}^{T T} & {m m}_{D D. 11}^{T T} \\ {n no}_{U u 11}^{T T} & {- - m m}_{U u 11}^{T T} \end{matrix}) (\begin{matrix} {m m}_{U u 22} & {m m}_{D D. 22} \\ {n no}_{U u 22} & {n no}_{D D. 22} \end{matrix}) (\begin{matrix} {v v}_{U u 22} \\ {v v}_{D D. 22} \end{matrix}) - - - - - - ((34.14 34.14))$

考虑下行波系统，同时包含P波和SV波，由介质1入射到分界面上，产生反射上行P、SV波和透射下行P、SV波，由于在介质2中无上行波，所以v_U2＝0。因此得到Consider the downgoing wave system, including both P waves and SV waves, which are incident on the interface from medium 1 to generate reflected upgoing P, SV waves and transmitted downgoing P, SV waves. Since there is no upgoing wave inmedium 2, v_U2 = 0. thus get

$(\begin{matrix} {v v}_{U u 11} \\ {v v}_{D D. 11} \end{matrix}) = = (\begin{matrix} {Q Q}_{UU UU} & {Q Q}_{UD UD} \\ {Q Q}_{DU DU} & {Q Q}_{DD DD} \end{matrix}) (\begin{matrix} 00 \\ {v v}_{D D. 22} \end{matrix}) - - - - - - ((2.2.15 2.2.15))$

引入下行波反射系数矩阵及透射系数

则有Introducing the downgoing wave reflection coefficient matrix and transmission coefficient

then there is

${v v}_{U u 11} = = {R R}_{D D.}^{I I} {v v}_{D D. 11}$

${v v}_{D D. 22} = = {T T}_{D D.}^{I I} {v v}_{D D. 11} - - - - - - ((2.2 2.2 . . 66))$

对比(2.2.15)和(2.2.16)，得到Comparing (2.2.15) and (2.2.16), we get

${R R}_{D D.}^{I I} = = {Q Q}_{UD UD} {(({Q Q}_{DD DD}))}^{- - 11}$

${T T}_{D D.}^{I I} = = {(({Q Q}_{DD DD}))}^{- - 11} - - - - - - ((2.2 2.2 . . 1717))$

考虑上行波系统，由介质2入射到分界面上，产生反射下行P、SV波和透射上行P、SV波，由于在介质1中无下行波，所以v_D1＝0。因此得到Considering the upgoing wave system,medium 2 is incident on the interface to generate reflected downgoing P, SV waves and transmitted upgoing P, SV waves. Since there is no downgoing wave in medium 1, v_D1 =0. thus get

$(\begin{matrix} {v v}_{U u 11} \\ 00 \end{matrix}) = = (\begin{matrix} {Q Q}_{UU UU} & {Q Q}_{UD UD} \\ {Q Q}_{DU DU} & {Q Q}_{DD DD} \end{matrix}) (\begin{matrix} {v v}_{U u 22} \\ {v v}_{D D. 22} \end{matrix}) - - - - - - ((2.2.18 2.2.18))$

引入上行波反射系数矩阵

及透射系数矩阵

则有Introducing the up-going wave reflection coefficient matrix

and the transmission coefficient matrix

then there is

${v v}_{D D. 22} = = {R R}_{U u}^{I I} {v v}_{U u 22}$

${v v}_{U u 11} = = {T T}_{U u}^{I I} {v v}_{U u 22} - - - - - - ((2.2 2.2 . . 1919))$

对比(2.2.18)和(2.2.19)，得到Comparing (2.2.18) and (2.2.19), we get

${R R}_{U u}^{I I} = = - - {(({Q Q}_{DD DD}))}^{- - 11} {Q Q}_{UD UD}$

${T T}_{U u}^{I I} = = {Q Q}_{UU UU} - - {Q Q}_{UD UD} {(({Q Q}_{DD DD}))}^{- - 11} {Q Q}_{DU DU} - - - - - - ((2.2.20 2.2.20))$

由(2.2.17)和(2.2.20)可以得到Q的另一表达式Another expression of Q can be obtained from (2.2.17) and (2.2.20)

$Q Q = = {D D.}_{11}^{- - 11} {D D.}_{22} = = (\begin{matrix} {T T}_{U u}^{I I} - - {R R}_{D D.}^{I I} {(({T T}_{D D.}^{I I}))}^{- - 11} {R R}_{U u}^{I I} & {R R}_{D D.}^{I I} {(({T T}_{D D.}^{I I}))}^{- - 11} \\ - - {(({T T}_{D D.}^{I I}))}^{- - 11} {R R}_{U u}^{I I} & {(({T T}_{D D.}^{I I}))}^{- - 11} \end{matrix}) - - - - - - ((2.2.21 2.2.21))$

所有的界面系数矩阵都依赖于

所以它会出现在反射及透射系数的表达式中，在各向同性介质中，有All interface coefficient matrices depend on

So it will appear in the expressions of reflection and transmission coefficients. In isotropic media, we have

Δ_st＝detQ_DD＝ε_α1ε_α2ε_β1ε_β2_Δst ＝detQ_DD ＝ε_α1 ε_α2 ε_β1 ε_β2

×{[2p²Δμ(q_α1-q_α2)+(ρ₁q_α2+ρ₂q_α1)]×{[2p² Δμ(q_α1 -q_α2 )+(ρ₁ q_α2 +ρ₂ q_α1 )]

×[2p²Δμ(q_β1-q_β2)+(ρ₁q_β2+ρ₂q_β1)]×[2p² Δμ(q_β1 -q_β2 )+(ρ₁ q_β2 +ρ₂ q_β1 )]

+p²[2Δμ(q_α1q_β2+p²)-Δρ][2Δμ(q_β1q_α2+p²)-Δρ]} (2.2.22)+p² [2Δμ(q_α1 q_β2 +p² )-Δρ][2Δμ(q_β1 q_α2 +p² )-Δρ]} (2.2.22)

其中in

Δμ＝μ₁-μ₂，Δρ＝ρ₁-ρ₂ (2.2.23)Δμ=μ₁ -μ₂ , Δρ=ρ₁ -ρ₂ (2.2.23)

各向同性介质中下行波的透射及反射系数表达式为The expressions of the transmission and reflection coefficients of the downgoing wave in the isotropic medium are

$(\begin{matrix} {T T}_{D D.}^{PP PP} & {T T}_{D D.}^{PS P.S.} \\ {T T}_{D D.}^{SP SP} & {T T}_{D D.}^{SS SS} \end{matrix}) = = \frac{11}{{Δ Δ}_{st st}} (\begin{matrix} {q q}_{4444} & {- - q q}_{3434} \\ {- - q q}_{1313} & {q q}_{3333} \end{matrix}) - - - - - - ((2.2.24 2.2.24))$

$(\begin{matrix} {R R}_{D D.}^{PP PP} & {R R}_{D D.}^{PS P.S.} \\ {R R}_{D D.}^{SP SP} & {R R}_{D D.}^{SS SS} \end{matrix}) = = \frac{11}{{Δ Δ}_{st st}} (\begin{matrix} {q q}_{1313} {q q}_{4444} - - {q q}_{1414} {q q}_{3434} & {q q}_{1414} {q q}_{3333} - - {q q}_{1313} {q q}_{3434} \\ {q q}_{23 twenty three} {q q}_{4444} - - {q q}_{24 twenty four} {q q}_{4343} & {q q}_{24 twenty four} {q q}_{3333} - - {q q}_{23 twenty three} {q q}_{3434} \end{matrix}) - - - - - - ((2.2.25 2.2.25))$

其中in

q₁₃＝ε_α1ε_α2[2p²Δμ(q_α1+q_α2)-(ρ₁q_α2-ρ₂q_α1)]，q₁₃ ＝ε_α1 ε_α2 [2p² Δμ(q_α1 +q_α2 )-(ρ₁ q_α2 -ρ₂ q_α1 )],

q₁₄＝iε_α1ε_β2p[2Δμ(q_α1q_β2-p²)+Δρ]，q₁₄ = iε_α1 ε_β2 p[2Δμ(q_α1 q_β2 -p² )+Δρ],

q₂₃＝iε_β1ε_α2p[2Δμ(q_β1q_α2-p²)+Δρ]，q₂₃ = iε_β1 ε_α2 p[2Δμ(q_β1 q_α2 -p² )+Δρ],

q₂₄＝ε_β1ε_β2[2p²Δμ(q_β1+q_β2)-(ρ₁q_β2-ρ₂q_β1)]，q₂₄ =ε_β1 ε_β2 [2p² Δμ(q_β1 +q_β2 )-(ρ₁ q_β2 -ρ₂ q_β1 )],

q₃₃＝ε_α1ε_α2[2p²Δμ(q_α1-q_α2)+(ρ₁q_α2+ρ₂q_α1)]，q₃₃ =ε_α1 ε_α2 [2p² Δμ(q_α1 -q_α2 )+(ρ₁ q_α2 +ρ₂ q_α1 )],

q₃₄＝iε_α1ε_β2p[2Δμ(q_α1q_β2+p²)-Δρ]，q₃₄ =iε_α1 ε_β2 p[2Δμ(q_α1 q_β2 +p² )-Δρ],

q₄₃＝iε_β1ε_α2p[2Δμ(q_β1q_α2+p²)-Δρ]，q₄₃ = iε_β1 ε_α2 p[2Δμ(q_β1 q_α2 +p² )-Δρ],

q₄₄＝ε_β1ε_β2[2p²Δμ(q_β1-q_β2)+(ρ₁q_β2+ρ₂q_β1)]， (2.2.26)q₄₄ =ε_β1 ε_β2 [2p² Δμ(q_β1 -q_β2 )+(ρ₁ q_β2 +ρ₂ q_β1 )], (2.2.26)

对于任意区域for any region

图3为任一区域的反射和透射情况的示意图。考虑下行波v_D从界面Z_a入射到区域AB时，见图3所示，该区域假设嵌在两个均匀半空间a和b中，那么定义在Z_a的上行反射波用反射矩阵

描述，那么反射回a区域的波场可以表示为：Fig. 3 is a schematic diagram of reflection and transmission of any region. Considering that the downgoing wave v_D is incident on the area AB from the interface Z_a , as shown in Fig. 3, the area is assumed to be embedded in two uniform half-spaces a and b, then the upgoing reflected wave defined in Z_a is defined by the reflection matrix

Description, then the wave field reflected back to region a can be expressed as:

${v v}_{U u} (({z z}_{A A})) = = {R R}_{D D.}^{AB AB} {v v}_{D D.} (({z z}_{A A})) - - - - - - ((2.3.1 2.3.1))$

同样的定义在z_b的下行透射波用透射矩阵

描述，那么透射到b区域的下行波可以表示为：The same definition uses the transmission matrix for the downgoing transmitted wave of z_b

Description, then the downgoing wave transmitted to region b can be expressed as:

${v v}_{D D.} (({z z}_{B B})) = = {T T}_{D D.}^{AB AB} {v v}_{D D.} (({z z}_{A A})) - - - - - - ((2.3.2 2.3.2))$

同理可以用来矩阵RT来表示从界面B入射到AB区域的上行波的反射和透射。In the same way, the matrix RT can be used to represent the reflection and transmission of the upgoing wave incident from the interface B to the region AB.

对于耦合的P-SV波，反射矩阵和透射矩阵可以扩展为：For coupled P-SV waves, the reflection and transmission matrices can be extended as:

${R R}_{D D.} = = {(\begin{matrix} {R R}_{pp pp} & {R R}_{ps ps} \\ {R R}_{sp sp} & {R R}_{ss ss} \end{matrix})}_{D D.},, {T T}_{D D.} = = {(\begin{matrix} {T T}_{pp pp} & {T T}_{ps ps} \\ {T T}_{sp sp} & {T T}_{ss ss} \end{matrix})}_{D D.} - - - - - - ((2.3.3 2.3.3))$

假设A和B之间是个均匀的区域，而且该区域和它嵌入的区域半空间AB是同一均匀介质这一特殊情况。显然在波入射到该区域时将没有反射，即：Suppose there is a uniform region between A and B, and this region and its embedded region half-space AB are the special case of the same homogeneous medium. Clearly there will be no reflection when the wave is incident on this region, ie:

${R R}_{D D.}^{AB AB} = = {R R}_{U u}^{AB AB} = = 00 - - - - - - ((2.3.4 2.3.4))$

那么此时的透射项将有：Then the transmission item at this time will have:

${T T}_{D D.}^{AB AB} = = {E E.}_{D D.}^{AB AB},, {T T}_{U u}^{AB AB} = = {E E.}_{U u}^{AB AB} - - - - - - ((2.3.5 2.3.5))$

这里的相位项

在依靠频率和慢度的P-SV波场是通过垂直慢度q_αq_β表示为2×2的矩阵。Phase term here

The frequency- and slowness-dependent P-SV wavefield is expressed as a 2×2 matrix by the vertical slowness q_α q_β .

将上行波场和下行波场在这一阶段区分开是很有必要的，虽然在各向同性的水平层状介质它们是相等的，但是对于各向异性或球面的层状介质它们是不等的。It is necessary to distinguish the upgoing wavefield from the downgoing wavefield at this stage. Although they are equal in isotropic horizontal layered media, they are not equal in anisotropic or spherical layered media. of.

首先，在单层区域。这里讨论一均匀层嵌入和它是不同介质的半空间AB时的情况。这时反射将在界面产生，而在均匀层的内部将只有一个相移的透射。First, in single-layer areas. Here we discuss the case when a uniform layer is embedded and it is a half-space AB of different media. In this case reflection will occur at the interface and there will be only a phase-shifted transmission inside the homogeneous layer.

考虑下行波从区域a入射到界面Z_a时，显然将产生反射用

表示，同时会有透射用

表示，透射波会发生相移E_D到达界面B，而且再次发生反射量Considering that the downgoing wave is incident on the interface Z a from region_a , it is obvious that there will be reflection

Indicates that there will be transmittance at the same time

Indicates that the transmitted wave will undergo a phase shift E_D to reach the interface B, and the reflected wave will again occur

${R R}_{D D.}^{B B} {E E.}_{D D.} {T T}_{D D.}^{A A} - - - - - - ((2.3.6 2.3.6))$

到均匀层和透射部分to the homogeneous layer and the transmitted part

${T T}_{D D.}^{B B} {E E.}_{D D.} {T T}_{D D.}^{A A} - - - - - - ((2.3.7 2.3.7))$

出界面B。Out of interface B.

伴随着在界面B的内部反射，波将会再次产生一个相移E_D。那么就会接着有量With internal reflection at interface B, the wave will again undergo a phase shift E_D . Then there will be a quantity

${T T}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} {T T}_{D D.}^{A A} - - - - - - ((2.3.8 2.3.8))$

透射出界面ZA和内部反射Transmission out interface ZA and internal reflection

${R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} {T T}_{D D.}^{A A} . . - - - - - - ((2.3.9 2.3.9))$

依次类推该过程将会一直延续下去，也就是所说的多次波。那么最后在界面A接收到的响应为：By analogy, this process will continue forever, which is the so-called multiple waves. Then the final response received on interface A is:

${R R}_{D D.}^{AB AB} = = {R R}_{D D.}^{A A} + + {T T}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} {T T}_{D D.}^{A A} + + {T T}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} {R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} {T T}_{D D.}^{A A} + + \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot;$

$= = {R R}_{D D.}^{A A} + + {T T}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} [[11 + + {R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} + + \cdot &Center Dot; \cdot \cdot \cdot &Center Dot;]] {T T}_{D D.}^{A A} - - - - - - ((2.3.10 2.3.10))$

同时透射到区域B的响应为：The response for simultaneous transmission into region B is:

${T T}_{D D.}^{AB AB} = = {T T}_{D D.}^{B B} {E E.}_{D D.} {T T}_{D D.}^{A A} + + {T T}_{D D.}^{B B} {E E.}_{D D.} {R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} {T T}_{D D.}^{A A} + + {T T}_{D D.}^{B B} {E E.}_{D D.} {R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} {R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} {T T}_{D D.}^{A A} + + \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot;$

$= = {T T}_{D D.}^{B B} {E E.}_{D D.} [[11 + + {R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} + + {R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} {R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} + + \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot;]] {T T}_{D D.}^{A A} - - - - - - ((2.3.11 2.3.11))$

事实上，由于层内的反复传播的反射和透射项可以看成是一个几何级数，也就是：In fact, due to the repeated propagation in the layer, the reflection and transmission terms can be regarded as a geometric progression, that is:

$[[11 + + {R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} + + {R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} {R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot;]] = = [[{11 - - {R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.}]]}^{- - 11} - - - - - - ((2.3.12 2.3.12))$

因此上面的公式就可以变成：So the above formula can become:

${R R}_{D D.}^{AB AB} = = {R R}_{D D.}^{A A} + + {T T}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.} {[[11 - - {R R}_{U u}^{A A} {E E.}_{U u} {R R}_{D D.}^{B B} {E E.}_{D D.}]]}^{- - 11} {T T}_{D D.}^{A A} - - - - - - ((2.3.13 2.3.13))$

从上面可以看到，可以在地层内的层间反射进行充分的计算，也就是对于所有的多次反射波可以得到正演。It can be seen from the above that the interlayer reflection in the formation can be fully calculated, that is, forward modeling can be obtained for all multiple reflection waves.

图4为单层区域的含多次波的演示示意图，见图4所示。另外，对于P波和SV波的相互转换也得到了考虑。Figure 4 is a schematic diagram of the demonstration of a single-layer area with multiple waves, as shown in Figure 4. In addition, the mutual conversion of P wave and SV wave has also been considered.

其次，在双层区域。考虑更为复杂一点的区域，在区域(A，C)中引入界面B，这样就考虑双层的情况。Second, in the double layer area. Considering a more complicated area, interface B is introduced in the area (A, C), so that the double-layer situation is considered.

图5为双层区域的传播演示示意图，见图5所示。在区域AC，可以定义反射和透射矩阵为：

和

Figure 5 is a schematic diagram of the propagation demonstration in the double-layer area, as shown in Figure 5 . In region AC, the reflection and transmission matrices can be defined as:

and

假设下行波[D^a]由A入射，那么来自AB的反射可以表示为：Assuming that the downgoing wave [D^a ] is incident by A, then the reflection from AB can be expressed as:

${R R}_{D D.}^{AB AB} [[{D D.}^{a a}]] - - - - - - ((2.3.14 2.3.14))$

透射过AB的部分为：The part transmitted through AB is:

${T T}_{D D.}^{AB AB} [[{D D.}^{a a}]] - - - - - - ((2.3.15 2.3.15))$

它将会成为BC区域的入射波。该部分将会一部分It will be the incident wave in the BC region. part of the

${R R}_{D D.}^{BC BC} {T T}_{D D.}^{AB AB} [[{D D.}^{a a}]] - - - - - - ((2.3.16 2.3.16))$

反射回区域AB，而另一部分reflected back to area AB, while another part

${T T}_{D D.}^{BC BC} {T T}_{D D.}^{AB AB} [[{D D.}^{a a}]] - - - - - - ((2.3.17 2.3.17))$

将会透射到区域C。will be transmitted into region C.

反射回区域AB的在界面B将会一部分透射出区域AB到区域AReflected back into region AB at interface B will be partially transmitted out of region AB into region A

${T T}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {T T}_{D D.}^{AB AB} [[{D D.}^{a a}]] - - - - - - ((2.3.18 2.3.18))$

另一部分another part

${R R}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {T T}_{D D.}^{AB AB} [[{D D.}^{a a}]] - - - - - - ((2.3.19 2.3.19))$

会反射回AB。will be reflected back to AB.

依次类推And so on

${R R}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {T T}_{D D.}^{AB AB} [[{D D.}^{a a}]]$

${T T}_{D D.}^{BC BC} {R R}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {T T}_{D D.}^{AB AB} [[{D D.}^{a a}]] - - - - - - ((2.3.20 2.3.20))$

透射出界面C，另一部分transmitted out of interface C, another part

${R R}_{D D.}^{BC BC} {R R}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {T T}_{D D.}^{AB AB} [[{D D.}^{a a}]] - - - - - - ((2.3.21 2.3.21))$

反射回B。这个过程将会不间断下去。Reflect back to B. This process will continue without interruption.

由此会得到界面A的总波场为：From this, the total wave field at interface A can be obtained as:

${R R}_{D D.}^{AC AC} [[{D D.}^{a a}]] = = {R R}_{D D.}^{AB AB} [[{D D.}^{a a}]] + + {T T}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {T T}_{D D.}^{AB AB} [[{D D.}^{a a}]] + + {T T}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {T T}_{D D.}^{BC BC} {R R}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {T T}_{D D.}^{AB AB} [[{D D.}^{a a}]] + + . . . . . . - - - - - - ((2.3.22 2.3.22))$

那么

可以表示为So

It can be expressed as

${R R}_{D D.}^{AC AC} = = {R R}_{D D.}^{AB AB} + + {T T}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {T T}_{D D.}^{AB AB} + + {T T}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {T T}_{D D.}^{BC BC} {R R}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {T T}_{D D.}^{AB AB} + + . . . . . . . . . . . . - - - - - - ((2.3.23 2.3.23))$

它还可以用简化的形式It can also be used in a simplified form

${R R}_{D D.}^{AC AC} = = {R R}_{D D.}^{AB AB} + + {T T}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {((I I - - {R R}_{U u}^{AB AB} {R R}_{D D.}^{BC BC}))}^{- - 11} {T T}_{D D.}^{AB AB} - - - - - - ((2.3.24 2.3.24))$

这里的I是个单位矩阵。而且Here I is an identity matrix. and

${((I I - - {R R}_{U u}^{AB AB} {R R}_{D D.}^{BC BC}))}^{- - 11} = = I I + + {R R}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} + + {R R}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {R R}_{U u}^{AB AB} {R R}_{D D.}^{BC BC} {((I I - - {R R}_{U u}^{AB AB} {R R}_{D D.}^{BC BC}))}^{- - 11} - - - - - - ((2.3.25 2.3.25))$

同样对于透射部分也得到了相似的结果Similar results were also obtained for the transmitted part

${T T}_{D D.}^{AC AC} = = {T T}_{D D.}^{BC BC} {((I I - - {R R}_{U u}^{AB AB} {R R}_{D D.}^{BC BC}))}^{- - 11} {T T}_{D D.}^{AB AB} - - - - - - ((2.3.26 2.3.26))$

在本实施例中，在得到单个地层和两个地层的计算公式以后，对于整个的地层模型就可以进行叠代计算，从而得到整个地层的反射和透射系数。In this embodiment, after obtaining the calculation formulas for a single formation and two formations, iterative calculations can be performed on the entire formation model, so as to obtain the reflection and transmission coefficients of the entire formation.

在本实施例中，可就具体的情况进行波场合成的讨论。即震源位于自由表面之下，震源位于水平线S处，自然可以以此线为界面将整个地层分为上下两部分。In this embodiment, the discussion of wave field synthesis can be carried out based on specific situations. That is, the seismic source is located under the free surface, and the seismic source is located at the horizontal line S. Naturally, the entire stratum can be divided into upper and lower parts with this line as the interface.

图6为地震波场的合成图，如图6所示。同时震源产生的波场也分为上行波场[U^S]和下行波场[D^S]。由此可以用前面提到的公式来表示区域fS和SL的波场。Fig. 6 is a synthetic diagram of the seismic wave field, as shown in Fig. 6 . At the same time, the wave field generated by the seismic source is also divided into the upgoing wave field [U^S ] and the downgoing wave field [D^S ]. The wavefields in the regions fS and SL can thus be represented by the formulas mentioned above.

首先，在自由表面接收的第一个波可以表示为：First, the first wave received at a free surface can be expressed as:

${W W}_{U u}^{fS f} {U u}^{S S} - - - - - - ((2.4.1 2.4.1))$

这里的是个矩阵，表示上行波在自由表面产生的位移。同时上行波也会有部分被区域fS反射回来为：here is a matrix representing the displacement produced by the upgoing wave on the free surface. At the same time, part of the upgoing wave will be reflected back by the area fS as:

${R R}_{U u}^{fS f} {U u}^{S S} - - - - - - ((2.4.2 2.4.2))$

这时总的下行波场将包括D^S和

可表示为：At this time the total downgoing wave field will include D^S and

Can be expressed as:

${D D.}^{S S} + + {R R}_{U u}^{fS f} {U u}^{S S} - - - - - - ((2.4.3 2.4.3))$

它将入射到区域SL，并发生第一次的反射为：It will be incident on area SL, and the first reflection will be:

${R R}_{D D.}^{SL SL} [[{D D.}^{S S} + + {R R}_{U u}^{fS f} {U u}^{S S}]] - - - - - - ((2.4.4 2.4.4))$

同时在表面将接收到该反射的位移贡献为：At the same time, the surface will receive the displacement contribution of the reflection as:

${W W}_{U u}^{fS f} {R R}_{D D.}^{SL SL} [[{D D.}^{S S} + + {R R}_{U u}^{fS f} {U u}^{S S}]] - - - - - - ((2.4 2.4 . . 55))$

然而之后还要考虑好多来自于表面的反射，而且每个在震源平面上下的一系列的反射都要用合成式

来表示。However, there are a lot of reflections from surfaces to consider later, and each series of reflections above and below the source plane must be synthesized

To represent.

这样可以将表面位移项表示为：This allows the surface displacement term to be expressed as:

${w w}_{00} = = {W W}_{U u}^{fS f} {U u}^{S S} + + {W W}_{U u}^{fS f} [[I I + + {R R}_{D D.}^{SL SL} {R R}_{U u}^{fS f} + + {R R}_{D D.}^{SL SL} {R R}_{U u}^{fS f} {R R}_{D D.}^{SL SL} {R R}_{U u}^{fS f} + + \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot;]] {R R}_{D D.}^{SL SL} [[{D D.}^{S S} + + {R R}_{U u}^{fS f} {U u}^{S S}]] - - - - - - ((2.4.6 2.4.6))$

它也可以简化为规则公式：It can also be reduced to a regular formula:

${w w}_{00} = = {W W}_{U u}^{fS f} {U u}^{S S} + + {W W}_{U u}^{fS f} {[[I I - - {R R}_{D D.}^{SL SL} {R R}_{U u}^{fS f}]]}^{- - 11} {R R}_{D D.}^{SL SL} [[{D D.}^{S S} + + {R R}_{U u}^{fS f} {U u}^{S S}]] - - - - - - ((2.4.7 2.4.7))$

对于震源在自由表面的情况，就可以将面S移到自由表面，此时下行波场

将只有D^S。For the case where the seismic source is on the free surface, the surface S can be moved to the free surface, at this time the downgoing wave field

There will be only^DS .

在本实施例中，积分的具体实现可如下所述。首先进行慢度积分，然后进行频率积分。这里来讨论慢度积分的数值计算。在方程：In this embodiment, the specific implementation of integration can be described as follows. Slowness integration is performed first, followed by frequency integration. Here we discuss the numerical computation of the slowness integral. In the equation:

$u u ((x x,, ω ω)) = = {&Integral; &Integral;}_{00}^{\infty \infty} {ω ω}^{22} pdp pdp {J J}_{00} ((ωpx ωpx)) \overset{~ ~}{u u} ((ω ω,, p p)) - - - - - - ((3.1.1 3.1.1))$

利用Bessel函数的远场近似：Far-field approximation using the Bessel function:

${J J}_{00} ((ωpx ωpx)) = = \frac{11}{22} [[{H h}_{00}^{((11))} ((ωpx ωpx)) + + {H h}_{00}^{((22))} ((ωpx ωpx))]] - - - - - - ((3.1.2 3.1.2))$

对于远场体波，第二项的贡献可以忽略不计，因此(3.1.1)就变成For far-field body waves, the contribution of the second term is negligible, so (3.1.1) becomes

$u u ((x x,, ω ω)) = = {&Integral; &Integral;}_{00}^{\infty \infty} \frac{11}{22} {ω ω}^{22} pdp pdp {H h}_{00}^{((11))} ((ωpx ωpx)) \overset{~ ~}{u u} ((ω ω,, p p)) - - - - - - ((3.1 3.1 . . 33))$

该积分是震荡的，而且当ωx很大时，要求许多的步长来准确计算它。因为步长的大小是与的ωx大小成反比的。因此对于较大的ω或x，方程(3.1.3)运算就会很耗时。The integral is oscillatory and requires many steps to compute it accurately when ωx is large. Because the size of the step size is inversely proportional to the size of ωx. Therefore, for larger ω or x, the operation of equation (3.1.3) will be very time-consuming.

为了克服这一困难，可以采用Filon方法。该方法要求的步长大小与(ωx)^-1/2成比例，这使得对于一个给定的步长可以用较大的ωx。为了用广义Filon方法，将(3.1.1)式变形为：To overcome this difficulty, the Filon method can be used. The method requires a step size proportional to (ωx)^-1/2 , which allows larger ωx for a given step size. In order to use the generalized Filon method, formula (3.1.1) is transformed into:

${&Integral; &Integral;}_{Γ Γ} \frac{{ω ω}^{22}}{22} {pdpH pdpH}_{00}^{((11))} ((ωpx ωpx)) \overset{~ ~}{u u} ((ω ω,, p p)) - - - - - - ((3.1.4 3.1.4))$

Γ在这里是积分区间。

表示在1形式和0阶时的汉克尔函数。Γ here is the integration interval.

Represents the Hankel function in form 1 andorder 0.

${H h}_{00} ((ωpx ωpx)) = = ((\frac{22}{πωpx πωpx})) exp exp ((iωpx iωpx - - \frac{iπ iπ}{44})) - - - - - - ((3.1.5 3.1.5))$

当ωpx为0时

为奇异的，所以变换时要注意。When ωpx is 0

is singular, so care should be taken when transforming.

现在变换(3.1.4)式为：Now transform (3.1.4) into:

${&Integral; &Integral;}_{Γ Γ} f f ((p p)) {e e}^{sg sg ((p p))} dp dp - - - - - - ((3.1.6 3.1.6))$

这里here

$f f ((p p)) = = \frac{{ω ω}^{22}}{22} {pH pH}_{00}^{((11))} ((ωpx ωpx)) {e e}^{- - iωpx iωpx} \overset{~ ~}{u u} ((ω ω,, p p)),,$

s＝iωx，g(p)＝ps=iωx, g(p)=p

应用复化梯形积分公式于等式(3.1.6)，在a到b的范围得到：Applying the complex trapezoidal integral formula to Equation (3.1.6), in the range from a to b, we get:

${&Integral; &Integral;}_{a a}^{b b} f f ((p p)) {e e}^{sg sg ((p p))} dp dp = = \frac{11}{22} [[f f ((a a)) {e e}^{sg sg ((a a))} + + f f ((b b)) {e e}^{sg sg ((b b))}]] δp δp - - - - - - ((3.1.7 3.1.7))$

复化梯形积分公式对于任意给定的积分式可以采用离散化，将其积分区间按照一定的步长等分，然后对每一个步长内进行梯形近似来求面积，最后再进行求和的一种近似方法。The complex trapezoidal integral formula can be discretized for any given integral formula, and its integral interval is equally divided according to a certain step size, and then the area is calculated by trapezoidal approximation within each step size, and finally a summation is performed an approximation method.

当Im{s(g(b)-g(a))}≥π/4，该等式是不成立的，因为它是假定f(p)e^sg(p)在区间(a，b)用线性函数很好的近似的。如果假定f(p)和g(p)在区间(a，b)可以被不同的线形函数近似，那么就可以得到另外一种形式：When Im{s(g(b)-g(a))}≥π/4, the equation is not valid, because it is assumed that f(p)e^sg(p) is linear in the interval (a, b) function is a good approximation. If it is assumed that f(p) and g(p) can be approximated by different linear functions in the interval (a, b), then another form can be obtained:

${&Integral; &Integral;}_{a a}^{b b} f f ((p p)) {e e}^{Sg S g ((p p))} dp dp = = \{\begin{matrix} \frac{δp δp}{Sδg Sδg} [[δ δ (({fe fe}^{Sg S g})) - - \frac{δ δ ((f f)) δ δ (({e e}^{Sg S g}))}{Sδg Sδg}]] ((11)) \\ \frac{δp δp}{22} [[f f ((a a)) {e e}^{Sg S g ((a a))} + + f f ((b b)) {e e}^{Sg S g ((b b))}]] ((22)) \end{matrix}$

(1)δ(g)≠0，(2)δ(g)＝0 (3.1.8)(1) δ(g)≠0, (2) δ(g)=0 (3.1.8)

这里的δ(p)表示p₂-p₁，对任意给的函数h(p)，δ(h)表示h(b)-h(a)。那么δ(g)就表示g(p₂)-g(p₁)。当δ(g)→0时等式(3.1.8)的右边趋近于等式(3.1.7)的右边。Here δ(p) represents p₂ -p₁ , and for any given function h(p), δ(h) represents h(b)-h(a). Then δ(g) means g(p₂ )-g(p₁ ). When δ(g)→0, the right side of equation (3.1.8) tends to the right side of equation (3.1.7).

既然为了生成(3.1.8)等式用一阶多项式来近似f(p)，就将(3.1.8)等式的积分作为一阶Filon的方法。用一阶Filon方法来代替复化梯形积分公式的方法极大的提高了处理的质量。然而对相同的精度(3.1.8)允许更大的步长。步长越大，反射率函数要求算的越少。两种方法都给出了计算机上矢量化得很好的编码。Since f(p) is approximated by a first-order polynomial in order to generate the (3.1.8) equation, the integral of the (3.1.8) equation is used as the method of the first-order Filon. The method of replacing the complex trapezoidal integral formula with the first-order Filon method greatly improves the processing quality. However larger step sizes are allowed for the same precision (3.1.8). The larger the step size, the reflectivity function Ask for less. Both methods give well-vectorized encodings on the computer.

在本实施例中，反射率函数在计算的时候，必须通过傅立叶逆变换In this example, the reflectance function When calculating, it is necessary to pass the inverse Fourier transform

$u u ((x x,, t t)) = = \frac{11}{22 π π} {&Integral; &Integral;}_{- - \infty \infty}^{\infty \infty} {dωe dωe}^{- - iωt iωt} u u ((x x,, ω ω)) - - - - - - ((3.2.1.1 3.2.1.1))$

变换到时间-空间域(x-t)空间。然而在实际的情况下，要受带宽的限制。截断积分等式为：Transform to time-space domain (x-t) space. However, in practical situations, it is limited by the bandwidth. The truncated integral equation is:

$\frac{11}{22 π π} {&Integral; &Integral;}_{- - Ω Ω}^{Ω Ω} {e e}^{- - iωt iωt} {&Integral; &Integral;}_{00}^{p p} {ω ω}^{22} pdp pdp \overset{~ ~}{u u} ((ω ω,, p p)) {J J}_{00} ((ωpx ωpx)) - - - - - - ((3.2.1.2 3.2.1.2))$

截断的影响是要将函数与一个核函数进行褶积，该核函数使信号发生畸变，并且产生了截断相位。结果是脉冲响应也产生了畸变。The effect of truncation is to convolve the function with a kernel function that distorts the signal and produces a truncated phase. The result is that the impulse response is also distorted.

为改进这一效果，必须用汉明窗使信号在边界处逐渐变陡，在频率域将信号与汉明窗相乘。To improve this effect, the signal must be steepened at the boundaries with a Hamming window, and the signal is multiplied by the Hamming window in the frequency domain.

图7为频率汉明窗和慢度汉明窗示意图。如图7所示，(a)为频率汉明窗，(b)为慢度汉明窗。窗函数为：Fig. 7 is a schematic diagram of a frequency Hamming window and a slowness Hamming window. As shown in Figure 7, (a) is the frequency Hamming window, and (b) is the slowness Hamming window. The window function is:

${W W}_{H h} ((\frac{ω ω}{Ω Ω})) = = \frac{11}{22} [[11 + + cos cos ((\frac{ωπ ωπ}{Ω Ω}))]],, | | ω ω | | < < Ω Ω$

${W W}_{H h} ((\frac{H h}{Ω Ω})) = = 00,, | | ω ω | | > > Ω Ω - - - - - - ((3.2.1.3 3.2.1.3))$

对慢度积分，可以用双边窗，范围分别从p₁到p₂和从p₃到p₄。p₁，p₄别是慢度的最大和最小值。p₂和p₃是p₁和p₄之间的慢度值。p₂通常要比最大速度对应的慢度值要大，而p₃要比最低速度对应的慢度值要小。对应这样一个窗的数学表达式为：For slowness integration, bilateral windows can be used, ranging from p₁ to p₂ and from p₃ to p₄ . p₁ , p₄ are the maximum and minimum values of slowness. p₂ and p₃ are slowness values between p₁ and p₄ . p₂ is usually larger than the slowness value corresponding to the maximum speed, and p₃ is smaller than the slowness value corresponding to the lowest speed. The mathematical expression corresponding to such a window is:

${P P}_{H h} ((p p)) = = \frac{11}{22} ((11 + + cos cos π π \frac{{p p}_{22} - - p p}{{p p}_{22} - - {p p}_{11}})),, {p p}_{11} \leq \leq p p \leq \leq {p p}_{22}$

P_H(p)＝1，p₂＜p＜p₃P_H (p) = 1, p₂ < p < p₃

${P P}_{H h} ((p p)) = = \frac{11}{22} ((11 + + cos cos π π \frac{p p - - {p p}_{33}}{{p p}_{44} - - {p p}_{33}})),, {p p}_{33} \leq \leq p p \leq \leq {p p}_{44} - - - - - - ((3.2.1.4 3.2.1.4))$

在慢度域，将该汉宁窗与反射率函数相乘就可以消除截断效应。加窗后的记录在较小的慢度和较大的慢度处的同相轴振幅变得较弱，且消除了假频，从而中间慢度处的同相轴更加清晰。In the slowness domain, the truncation effect can be eliminated by multiplying the Hanning window with the reflectance function. The windowed recordings have weaker event amplitudes at smaller slownesses and larger slownesses, and aliasing has been eliminated, resulting in clearer events at intermediate slownesses.

下面对反射率法的实现加以试算，模拟采用的震源子波是Ricker子波。Ricker子波是波场数值模拟中最常用的一种零相位子波，在波场数值模拟中起到了不可替代的作用。The following is a trial calculation of the realization of the reflectivity method. The source wavelet used in the simulation is the Ricker wavelet. Ricker wavelet is the most commonly used zero-phase wavelet in wave field numerical simulation, and it plays an irreplaceable role in wave field numerical simulation.

Ricker子波的时间域表达式为The time domain expression of the Ricker wavelet is

$s the s ((t t)) = = ((11 - - {22 π π}^{22} {f f}_{M m}^{22} {t t}^{22})) {e e}^{- - {π π}^{22} {f f}_{M m}^{22} {t t}^{22}} - - - - - - ((5.1.1 5.1.1))$

对上式进行傅立叶变换，得到Ricker子波的频率域表达式Perform Fourier transform on the above formula to get the frequency domain expression of the Ricker wavelet

$S S ((f f)) = = \frac{22}{\sqrt{π π}} \frac{{f f}^{22}}{{f f}_{M m}^{33}} {e e}^{- - {f f}^{22} / / {f f}_{M m}^{22}} - - - - - - ((5.1.2 5.1.2))$

其中，t为时间，f为频率，f_M为主频。Among them, t is time, f is frequency, and f_M is the main frequency.

图8为Ricker子波时间域波形示意图；图9为Ricker子波频率域波形示意图。FIG. 8 is a schematic diagram of Ricker wavelet time domain waveform; FIG. 9 is a schematic diagram of Ricker wavelet frequency domain waveform.

图10为由模拟测井数据提取的反射系数、褶积记录、反射率法正演结果的示意图。由图10的结果对比可知，褶积模拟的结果只是正确的模拟了一次反射结果，而反射率法模拟的结果包含了转换波、一次多次波以及层间多次波，从简单模型的实验结果可以看出反射率实现比较准确，方法效果不错，可以用于实际数据模拟的井地匹配。Fig. 10 is a schematic diagram of the reflection coefficient extracted from the simulated logging data, the convolution record, and the forward modeling results of the reflectivity method. From the comparison of the results in Figure 10, it can be seen that the result of the convolution simulation is only the correct simulation of the primary reflection result, while the simulation result of the reflectivity method includes converted waves, primary multiples and interlayer multiples. From the simple model experiment It can be seen from the results that the reflectivity is relatively accurate and the method works well, which can be used for well-ground matching of actual data simulation.

下面就一个简单的三层的水平层状介质模型，进行地震记录合成。其具体的厚度，密度和纵横波速度参数如表1所示。The following is a simple three-layer horizontal layered medium model for synthesizing seismic records. Its specific thickness, density and P-S wave velocity parameters are shown in Table 1.

表1Table 1

层 layer 纵波速度(m/s)P-wave velocity (m/s) 横波速度(m/s)Shear wave velocity (m/s) 密度(m/cm³)Density(m/cm³ ) 厚度(m)Thickness (m) 1 1 24002400 12001200 2.02.0 600600 2 2 28002800 14001400 2.12.1 600600 33 33003300 17001700 2.22.2 800800 44 38003800 21002100 2.32.3

在计算中采用了主频为35HZ的雷克子波来模拟炸药震源，记录的采样率为2ms，时间长度是3.5秒，道数为60道，道距为30m。模型专门设计的三个纵波反射系数分别出现在0.5秒，0.9秒和1.4秒。In the calculation, the Reker wavelet with the main frequency of 35HZ is used to simulate the source of explosives. The sampling rate of the record is 2ms, the time length is 3.5 seconds, the number of traces is 60, and the trace distance is 30m. The three longitudinal wave reflection coefficients specially designed by the model appear at 0.5 seconds, 0.9 seconds and 1.4 seconds respectively.

图11为三层模型一次波、多次波和转换波记录的示意图；从图11可以看出，这三个一次反射波平P1，P2，P3都在相应的时刻准确出现。三层的转换波SV1，SV2，SV3的出现时间应该在0.75秒，1.4秒和2.1秒，也都准确无误。Figure 11 is a schematic diagram of the primary wave, multiple wave and converted wave records of the three-layer model; it can be seen from Figure 11 that the three primary reflected wave levels P1, P2, and P3 all appear accurately at the corresponding time. The appearance time of the three layers of conversion waves SV1, SV2, and SV3 should be 0.75 seconds, 1.4 seconds and 2.1 seconds, and they are all accurate.

在本实施例中，还对反射率法正演的井地匹配的实际效果进行分析，模型为大港油田一条测井曲线得到的数据。In this embodiment, the actual effect of the well-ground matching in the forward modeling of the reflectance method is also analyzed, and the model is the data obtained from a logging curve in Dagang Oilfield.

图12为np5-4井速度和密度测井资料及井旁道地震记录的示意图，见图12所示，左边分别是纵波速度和密度曲线，深度300米到5300米，横波速度由经验公式求取。该模型最大的特点是中间1800ms处速度密度变化剧烈，反射系数很强。震源还是用主频为35赫兹的Ricker子波。对数据进行反射率法井地匹配模拟。Figure 12 is a schematic diagram of well np5-4 velocity and density logging data and well side channel seismic records, as shown in Figure 12, the left side is the compressional wave velocity and density curve, the depth is 300 meters to 5300 meters, the shear wave velocity is obtained by empirical formula Pick. The biggest feature of this model is that the velocity density changes sharply at the middle 1800ms, and the reflection coefficient is very strong. The source still uses the Ricker wavelet with a main frequency of 35 Hz. Well-ground matching simulation is carried out on the data.

图13为对np5-4井测井数据分层及用反射率法模拟记录的示意图，由图13模拟的结果对比图11中的由褶积方法获得的的井旁道记录可以看出，由反射率法模拟的记录与测井数据匹配的效果更好。Fig. 13 is a schematic diagram of layering logging data of well np5-4 and simulating records by reflectivity method. It can be seen from the simulated results in Fig. 13 compared with the well side channel records obtained by the convolution method in Fig. 11 that by The records simulated by the reflectance method are better matched with the logging data.

由上述实施例可知，通过在进行测井数据与地震数据匹配时，同时考虑多角度反射地震波、层间多次波及界面转换波，可以提升井地数据匹配精度。It can be seen from the above embodiments that, when matching well logging data and seismic data, considering multi-angle reflected seismic waves, interlayer multiple waves and interface converted waves simultaneously, the matching accuracy of well ground data can be improved.

本发明实施例还提供一种测井曲线与地震记录全波匹配的装置，如图14所示，所述装置包括：获取单元1401、处理单元1402和匹配单元1403；其中，The embodiment of the present invention also provides a device for full-wave matching of logging curves and seismic records, as shown in Figure 14, the device includes: anacquisition unit 1401, aprocessing unit 1402 and amatching unit 1403; wherein,

获取单元1401获取地震记录；Theacquisition unit 1401 acquires seismic records;

处理单元1402利用反射率法对地震记录进行处理，获得井旁道地震记录；Theprocessing unit 1402 uses the reflectivity method to process the seismic records to obtain the well bypass seismic records;

匹配单元1403将井旁道地震记录与测井曲线进行匹配。Thematching unit 1403 matches the side channel seismic records with the logging curves.

进一步地，如图15所示，处理单元1402具体可包括：求解单元1501和积分单元1502；其中，Further, as shown in FIG. 15 , theprocessing unit 1402 may specifically include: a solvingunit 1501 and an integratingunit 1502; wherein,

求解单元1501在频率慢度域求解反射透射系数响应；Thesolving unit 1501 solves the reflection and transmission coefficient response in the frequency slowness domain;

积分单元1502根据反射透射系数响应进行慢度积分和频率积分，得到井旁地震记录。Theintegration unit 1502 performs slowness integration and frequency integration according to the reflection and transmission coefficient responses to obtain side-hole seismic records.

进一步地，求解单元1501具体用于：分别求取自由界面、固体-固体界面和任意区域的反射透射系数；将反射透射系数合成为总的反射透射系数响应。Further, thesolving unit 1501 is specifically configured to: respectively calculate the reflection and transmission coefficients of the free interface, the solid-solid interface and any region; synthesize the reflection and transmission coefficients into a total reflection and transmission coefficient response.

进一步地，积分单元1502具体用于：在频率域将信号与汉明窗相乘；在慢度域采用汉明窗与反射率函数相乘来消除截断效应。Further, the integratingunit 1502 is specifically configured to: multiply the signal by the Hamming window in the frequency domain; and multiply the Hamming window by the reflectivity function in the slowness domain to eliminate the truncation effect.

本领域普通技术人员还可以进一步意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本发明的范围。Those of ordinary skill in the art can further appreciate that the units and algorithm steps of the examples described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, computer software, or a combination of the two. In order to clearly illustrate the hardware and Interchangeability of software. In the above description, the components and steps of each example have been generally described according to their functions. Whether these functions are executed by hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art may use different methods to implement the described functions for each specific application, but such implementation should not be regarded as exceeding the scope of the present invention.

结合本文中所公开的实施例描述的方法或算法的步骤可以用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储器(RAM)、内存、只读存储器(ROM)、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM、或技术领域内所公知的任意其它形式的存储介质中。The steps of the methods or algorithms described in connection with the embodiments disclosed herein may be implemented by hardware, software modules executed by a processor, or a combination of both. Software modules can be placed in random access memory (RAM), internal memory, read-only memory (ROM), electrically programmable ROM, electrically erasable programmable ROM, registers, hard disk, removable disk, CD-ROM, or any other Any other known storage medium.

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。The specific embodiments described above have further described the purpose, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above descriptions are only specific embodiments of the present invention and are not intended to limit the scope of the present invention. Protection scope, within the spirit and principles of the present invention, any modification, equivalent replacement, improvement, etc., shall be included in the protection scope of the present invention.

Claims

1. A method of full wave matching a log to a seismic record, the method comprising:

acquiring a seismic record;

processing the seismic record by using a reflectivity method to obtain a well side channel seismic record;

and matching the well side channel seismic records with the logging curves.

2. The method of claim 1, wherein the processing the seismic records using a reflectivity method to obtain well side-channel seismic records comprises:

solving the reflection transmission coefficient response in a frequency slowness domain;

and performing slowness integration and frequency integration according to the reflection and transmission coefficient response to obtain the well side seismic record.

3. The method of claim 2, wherein solving for the reflection-transmission coefficient response in the frequency slowness domain comprises:

respectively obtaining the reflection and transmission coefficients of a free interface, a solid-solid interface and any area;

the reflection-transmission coefficients are combined into a total reflection-transmission coefficient response.

4. The method of claim 2, wherein performing slowness integration and frequency integration based on the reflection-transmission coefficient response comprises:

multiplying the signal with a hamming window in the frequency domain;

and in the slowness domain, a Hamming window is multiplied by the reflectivity function to eliminate the truncation effect.

5. An apparatus for full wave matching of a log to a seismic record, the apparatus comprising:

an acquisition unit that acquires a seismic record;

the processing unit is used for processing the seismic record by utilizing a reflectivity method to obtain a well side channel seismic record;

and the matching unit is used for matching the well side channel seismic record with the logging curve.

6. The apparatus according to claim 5, wherein the processing unit specifically comprises:

the solving unit is used for solving the response of the reflection and transmission coefficients in a frequency slowness domain;

and the integration unit is used for performing slowness integration and frequency integration according to the reflection and transmission coefficient response to obtain the well side seismic record.

7. The apparatus according to claim 6, wherein the solving unit is specifically configured to: respectively obtaining the reflection and transmission coefficients of a free interface, a solid-solid interface and any area; the reflection-transmission coefficients are combined into a total reflection-transmission coefficient response.

8. The apparatus of claim 6, wherein the integration unit is specifically configured to: multiplying the signal with a hamming window in the frequency domain; and in the slowness domain, a Hamming window is multiplied by the reflectivity function to eliminate the truncation effect.