CN103149592A

CN103149592A - Method for separating variable offset vertical seismic profile (VSP) wave fields

Info

Publication number: CN103149592A
Application number: CN201310070899XA
Authority: CN
Inventors: 胡建平; 戴华林; 苑伯达; 高树成; 陈亚东; 王丽
Original assignee: Henan University of Urban Construction
Current assignee: Henan University of Urban Construction
Priority date: 2013-03-07
Filing date: 2013-03-07
Publication date: 2013-06-12

Abstract

A variable-offset VSP wavefield separation method, comprising: inputting the collected original seismic wavefield data of the variable-offset VSP into a computer; performing ^t2 stretch processing on the original seismic wavefield data to form ^t2 Stretch processing data; carry out discrete positive parabolic Radon transformation to described ^t2 stretching processing data, form τ-q domain uplink and downlink wave division data; Describe τ-q domain uplink and downlink wave division data and carry out discrete inversion Parabolic Radon transformation to form xt ² domain data; reverse t ² stretch processing is performed on the xt ² domain data to complete the separation of the up-going wave and down-going wave of the variable offset VSP.

Description

A VSP Wavefield Separation Method with Variable Offset

技术领域technical field

本发明涉及地震勘探数据处理技术，特别是一种变偏移距VSP波场的上行波和下行波分离方法。The invention relates to a seismic exploration data processing technology, in particular to a method for separating upgoing waves and downgoing waves of a variable-offset VSP wave field.

背景技术Background technique

垂直地震剖面(Vertical Seismic Profiling)是与地面观测的地震剖面相对的一种观测方法。地面观测的地震剖面是在地表附近的一些点上激发，同时沿地面布置检波器进行观测；垂直地震剖面也是在地表附近的一些点上激发地震波，但是它是在沿井孔不同深度布置一些检波器进行观测，这些检波器都是三分量的，接收的是质点振动的矢量。垂直地震剖面实际上是早已广泛使用的地震测井方法的变革和发展。它现在已经远远超出了地震测井原来的范围，而发展成为一套完整的、独立的观测方法。垂直地震剖面由于与地面观测相比有它一些固有的优点，所以得到日益广泛的应用，近二十年来，它已经成为勘探地球物理的重要方法。利用VSP可以深入的了解地震子波传播的某些基本特性，帮助了解反射和透射过程，从而反过来又可改善地表地震资料关于构造、地层和岩性的解释。Vertical Seismic Profiling is an observation method opposite to the seismic profile of ground observation. The seismic section of ground observation is excited at some points near the surface, and at the same time, geophones are arranged along the ground for observation; the vertical seismic section also excites seismic waves at some points near the surface, but it is arranged at different depths along the borehole These detectors are all three-component, and what they receive is the vector of particle vibration. The vertical seismic section is actually the transformation and development of the widely used seismic logging method. It has now gone far beyond the original scope of seismic logging, and has developed into a complete and independent observation method. Due to its inherent advantages compared with surface observations, vertical seismic sections have been increasingly used. In the past two decades, it has become an important method for exploration geophysics. The use of VSP can provide an in-depth understanding of some basic characteristics of seismic wavelet propagation, helping to understand the reflection and transmission processes, which in turn can improve the interpretation of surface seismic data on structure, stratum and lithology.

垂直地震剖面方法根据观测方式可分为零偏移距VSP和非零偏移距VSP，非零偏移距VSP又可以分为变偏移距VSP和常数偏移距VSP，后来发展了斜井VSP观测。The vertical seismic section method can be divided into zero offset VSP and non-zero offset VSP according to the observation method, and non-zero offset VSP can be divided into variable offset VSP and constant offset VSP. VSP observations.

变偏移距VSP方法是垂直地震剖面方法中的一种特殊观测方法，这种方法与常规VSP方法相比有波场复杂、波场分离困难的缺陷，但是由于该方法检波器可以固定在储集层的上方，更有利于描述储层性质及其纵横向变化情况。当有油气显示时，用变偏移距VSP探测漏掉的油气层是很有效的。The variable offset distance VSP method is a special observation method in the vertical seismic section method. Compared with the conventional VSP method, this method has the defects of complex wave field and difficult separation of wave fields. However, due to this method, the geophone can be fixed in the storage It is more conducive to describe the properties of reservoirs and their vertical and horizontal changes. When there are oil and gas shows, it is very effective to use the variable offset distance VSP to detect the missing oil and gas layers.

变偏移距VSP观测方法是：检波器固定在储集层上方的某一点，震源沿地表移动。如图1所示，由于检波器靠近勘探目标，因而更有利于研究储层。一般可选择下述情况的井进行观测。The VSP observation method with variable offset distance is: the geophone is fixed at a certain point above the reservoir, and the seismic source moves along the surface. As shown in Figure 1, since the geophone is close to the exploration target, it is more beneficial to study the reservoir. Generally, wells with the following conditions can be selected for observation.

(1)井已穿过储集层，而常规地震剖面方法又难以查清储层的边界。如图1(a)、(b)。(1) The well has passed through the reservoir, and it is difficult to find out the boundary of the reservoir by conventional seismic section methods. Figure 1 (a), (b).

(2)井没有打到真正的储层，而是在水平方向上有偏离，并需进一步钻探重新确定新井位时，如图1(c)。(2) The well does not hit the real reservoir, but deviates in the horizontal direction, and further drilling is required to re-determine the new well location, as shown in Figure 1(c).

变偏移距VSP观测中，检波器被固定在某一深度，震源距等间隔变化，这样就得到在时间-空间域内的波场记录。In variable offset distance VSP observation, the geophone is fixed at a certain depth, and the source distance is changed at equal intervals, so that the wave field record in the time-space domain is obtained.

其中，直达波方程为：Among them, the direct wave equation is:

h²+x²＝v²t² (1.1.1)h ² +x ² = v ² t ² (1.1.1)

上行波方程为：The upgoing wave equation is:

(2H-h)²+x²＝v²t² (1.1.2)(2H-h) ² +x ² ＝v ² t ² (1.1.2)

式中：h为检波点深度；x为井源距；H为目的层深度；v为速度。In the formula: h is the depth of the receiver point; x is the well-source distance; H is the depth of the target layer; v is the velocity.

分析以上两式可知，直达波方程与上行波方程在时空域内都表现为双曲线。由于Analyzing the above two equations, we can see that both the direct wave equation and the upgoing wave equation behave as hyperbolas in the space-time domain. because

$\underset{H h &RightArrow; &Right Arrow; h h}{lim lim} {((22 H h - - h h))}^{22} = = {h h}^{22}$

因此在检波点所在深度附近，上行波方程与直达波方程非常近似，两条曲线趋于重合。这样给上下行波分离带来了一定的难度。Therefore, near the depth where the receiver is located, the upgoing wave equation is very similar to the direct wave equation, and the two curves tend to coincide. This brings certain difficulties to the separation of uplink and downlink waves.

变偏移距VSP这门勘探技术自从提出来已经有20多年的历史了，其间，有大量的地球物理工作者对之进行了大量的研究工作，使得这一技术在上世纪末及本世纪初取得了长足的进展。勘探地球物理学会(SEG)每年的会议都会设立一个独立的专题——井中地球物理(VSP)进行讨论，但是，至今未能看到一个关于变偏移距VSP的成功应用，这主要是处理技术相对滞后的原因。截至目前，关于变偏移据VSP资料的处理，基本上局限在对于下行波场的利用，波场分离技术和成像技术讨论的也很多，但这些成果都是针对固定震源情况下的非零偏移距VSP资料的，关于多偏移距、多接收点的变偏移距VSP处理技术，尤其是上下行波的分离技术至今在国内外均未看到令人满意的结果。究其原因，就是因为变偏移距VSP资料的处理技术有着极大的挑战性的困难。The variable offset distance VSP exploration technology has been proposed for more than 20 years. During this period, a large number of geophysicists have done a lot of research work on it, making this technology a success at the end of the last century and the beginning of this century. made great progress. The annual meeting of the Society for Exploration Geophysics (SEG) will set up an independent topic - well geophysics (VSP) for discussion, but so far, we have not seen a successful application of variable offset VSP, which is mainly the processing technology The reason for the relative lag. Up to now, the processing of variable offset data according to VSP is basically limited to the use of downgoing wave field, and there are many discussions on wave field separation technology and imaging technology, but these achievements are all for non-zero offset in the case of fixed source. For offset VSP data, no satisfactory results have been found so far at home and abroad on the multi-offset, multi-receiving point variable-offset VSP processing technology, especially the separation technology of uplink and downlink waves. The reason is that the processing technology of variable offset distance VSP data has great challenges.

变偏移距VSP方法是垂直地震剖面法中的一种特殊观测方法，主要用于“透视井的周边”，探测井口附近是否有一些重要变化。当有油气显示时，用变偏移距VSP探测漏掉的油气层是很有效的。在国内，有关这方面研究的成果很少，再者由于这种方法与常规VSP方法相比有波场复杂、波场分离困难的缺陷，从而阻止了它发展的步伐。王成礼、焦湘恒1992提出了用接应点法处理小偏移距道集，用中值滤波处理大偏移距道集的方法。在我们的研究中，应用近几年来发展起来的抛物线拉冬变换方法，首先对理论射线追踪模型变偏移距VSP资料进行上下行波分离，进而推广到对实际资料进行处理。The variable offset distance VSP method is a special observation method in the vertical seismic section method, which is mainly used to "see through the periphery of the well" and detect whether there are some important changes near the wellhead. When there are oil and gas shows, it is very effective to use the variable offset distance VSP to detect the missing oil and gas layers. In China, there are very few research results in this area. Compared with the conventional VSP method, this method has the defects of complex wave field and difficult wave field separation, which hinders its development. Wang Chengli and Jiao Xiangheng proposed in 1992 the method of processing small offset gathers with the contact point method and using median filtering to process large offset gathers. In our research, using the parabolic Radon transform method developed in recent years, firstly the uplink and downlink waves are separated from the variable offset distance VSP data of the theoretical ray tracing model, and then extended to the actual data processing.

任何一种地震数字处理方法都经历了提出问题、实际检验、完善理论、再实践、再提高的一个螺旋式的上升过程，拉冬变换方法的发展亦是如此。拉冬变换是由Radon于1917年提出的，它为图像重构问题提供了一个统一的数学基础，目前已被广泛应用于物理、医学、天文、分子生物、材料科学、核磁共振、无损检测、地球物理等方面。1978年Chapman、Schultz和Claerbout的文章把拉冬变换引入到地球物理学，以后一些学者开始把它运用到各种勘探地震学问题中，例如速度分析、压制多次波、波场分离、偏移、反演和道内插等。1985年以来，许多研究者在拉冬变换问题上进行了卓有成效的研究工作，分别提出了频率域中的双曲型拉冬变换，基于剩余动校正和t2拉伸的抛物型拉冬变换，并且抛物型拉冬变换已经被应用到速度滤波和多次波压缩。1995年，Sacchi和Ulrych提出了高分辨率抛物拉冬变换，从而把拉冬变换又推上了一个新的高度。Any seismic digital processing method has gone through a spiral ascending process of raising questions, practical testing, perfecting theory, practice, and improvement, and the development of Radon transform method is the same. The Radon transform was proposed by Radon in 1917. It provides a unified mathematical basis for the image reconstruction problem. It has been widely used in physics, medicine, astronomy, molecular biology, material science, nuclear magnetic resonance, non-destructive testing, Geophysics etc. In 1978, the article of Chapman, Schultz and Claerbout introduced the Ladon transform into geophysics, and later some scholars began to apply it to various exploration seismology problems, such as velocity analysis, suppression of multiple waves, wave field separation, migration , inversion and interpolation etc. Since 1985, many researchers have carried out fruitful research work on the Radon transform problem, respectively proposed the hyperbolic Radon transform in the frequency domain, the parabolic Radon transform based on residual motion correction and t2 stretching, and The parabolic Radon transform has been applied to velocity filtering and multiple compression. In 1995, Sacchi and Ulrych proposed the high-resolution parabolic Radon transform, which pushed the Radon transform to a new height.

上行波和下行波分离不彻底。图2是一个变偏移距VSP的实际资料，图2中能量最强的同相轴几乎都是一次下行直达波及其多次下行波，因此，对于下行波的利用来说几乎没有什么困难，因为我们利用的是一次下行直达波，这是记录中最先到达的、能量最强的一组相位。对于上行波来说，它被掩没在下行波场之中，在动力学和运动学的范畴内，二者除了到达时间有区别外，波场的视速度、频谱范围都是十分接近的，因此，常规的处理技术，无论是一维滤波还是二维滤波都很难将二者分离出来。这是变偏移距处理技术要解决的一个核心问题。Incomplete separation of upgoing and downgoing waves. Fig. 2 is the actual data of a variable offset VSP. The events with the strongest energy in Fig. 2 are almost all downgoing direct waves and multiple downgoing waves. Therefore, there is almost no difficulty in utilizing the downgoing waves, because We took advantage of a downgoing direct wave, which is the first to arrive and the most energetic set of phases on record. For the upgoing wave, it is buried in the downgoing wave field. In the category of dynamics and kinematics, except for the difference in arrival time, the apparent velocity and spectrum range of the wave field are very close. Therefore, conventional processing techniques, whether it is one-dimensional filtering or two-dimensional filtering, are difficult to separate the two. This is a core problem to be solved by the variable offset distance processing technology.

综上所述，现在社会上急需一种变偏移距VSP上行波和下行波的分离方法。To sum up, there is an urgent need in the society for a method for separating upgoing waves and downgoing waves of variable-offset VSPs.

发明内容Contents of the invention

本发明的目的是提供一种变偏移距VSP波场分离方法，特别是一种离散抛物拉冬变换方法，它能克服现有变偏移距VSP分离的上述缺点，更好的分离变偏移距VSP上行波和下行波。The purpose of the present invention is to provide a variable offset distance VSP wave field separation method, especially a discrete parabolic Radon transform method, which can overcome the above-mentioned shortcomings of the existing variable offset distance VSP separation, and better separate variable deviation Shift VSP up-going and down-going waves.

本发明通过如下方法实现变偏移距VSP波场分离：The present invention realizes the separation of VSP wave field with variable offset by the following method:

一种变偏移距VSP波场分离方法包括：在计算机中输入采集到的变偏移距VSP的原始地震波场数据；对所述的原始地震波场数据进行t²拉伸处理，形成t²拉伸处理数据；对所述的t²拉伸处理数据进行对离散正抛物拉冬变换，形成τ-q域上下行波分数据；对所述τ-q域上下行波分数据进行离散反抛物拉冬变换，形成x-t²域数据；对所述x-t²域数据进行反t²拉伸处理，完成变偏移距VSP上行波和下行波的分离。A variable offset VSP wavefield separation method comprises: inputting the original seismic wavefield data of the variable offset VSP collected in a computer; performing ^t2 stretch processing on the original seismic wavefield data to form a ^t2pulse Stretch processing data; Carry out discrete forward parabolic Radon transformation to described ^t2 stretching processing data, form τ-q domain uplink and downlink wave division data; Describe τ-q domain uplink and downlink wave division data and carry out discrete anti-parabolic Radon transformation to form xt ² domain data; reverse t ² stretching processing is performed on the xt ² domain data to complete the separation of up-going waves and down-going waves of the variable offset VSP.

优选包括t²拉伸处理过程：对原始地震波场数据进行读取，输入相关参数，通过iflag＝1的公式进行判断，当iflag＝1成立，对输入参数的原始地震波场数据进行三次样条插值拉伸，形成x-t²数据文件；当iflag＝1不成立，对输入参数的原始地震波场数据进行三次样条反插值拉伸，形成x-t数据文件。形成t²拉伸后的地震数据。Preferably include ^t2 stretching process: read the original seismic wavefield data, input relevant parameters, judge by the formula of iflag=1, when iflag=1 is established, perform cubic spline interpolation on the original seismic wavefield data of the input parameters Stretch to form an xt ² data file; when iflag=1 is not established, perform cubic spline inverse interpolation stretching on the original seismic wave field data of the input parameters to form an xt data file. Form the seismic data after ^t2 stretching.

进一步优选包括离散正抛物拉冬变换的过程：读取t²拉伸后的地震数据，输入参数，计算q的扫描范围rq，对计算后的扫描范围rq进行判断，rq不满足要求返回到输入参数步骤，rq满足要求进行FFT处理，并对处理后的数据进行一维滤波；对滤波后的数据进行判断，当n＜m时利用M＝L^H(LLH+λ²I)^-1D对数据进行处理，当n≥m时利用M＝(L^HL+λ²I)^-1L^HD对数据进行处理，形成处理后的数据；对所述处理后的数据进行二维褶积，并对二维褶积后的数据进行一维滤波并进行IFFT处理，形成τ-q域数据文件，即形成τ-q域上下行波分离。It is further preferred to include the process of discrete positive parabolic Radon transformation: read the seismic data after ^t2 stretching, input parameters, calculate the scanning range rq of q, and judge the calculated scanning range rq, if rq does not meet the requirements, return to the input In the parameter step, rq meets the requirements for FFT processing, and one-dimensional filtering is performed on the processed data; the filtered data is judged, and when n<m, use M=L ^H (LLH+λ ² I) ^-1 D to The data is processed, and when n≥m, the data is processed by using M=(L ^H L+λ ² I) ^-1 ^L HD to form processed data; two-dimensional convolution is performed on the processed data, One-dimensional filtering and IFFT processing are performed on the data after two-dimensional convolution to form a τ-q domain data file, that is, the separation of uplink and downlink waves in the τ-q domain is formed.

进一步优选包括离散反抛物拉冬变换的过程：读取τ-q域数据文件，读取输入的参数文件，对数据进行FFT处理，对处理后的数据进行第一次频率域滤波，然后对第一次频率域滤波后的数据进行D＝LM处理，对处理后的数据进行第二次频率域滤波，对第二次频率域滤波后的数据进行IFFT处理，形成x-t²域数据文件。Further preferably comprise the process of discrete inverse parabolic Radon transformation: read the τ-q domain data file, read the input parameter file, carry out FFT processing to the data, carry out the frequency domain filter to the processed data for the first time, then to the second D=LM processing is performed on the data after the first frequency domain filtering, the second frequency domain filtering is performed on the processed data, and the IFFT processing is performed on the second frequency domain filtering data to form xt ² domain data files.

进一步优选包括反t²拉伸处理过程：对所述经过离散反抛物拉冬变换的数据进行读取，输入相关参数，通过iflag＝1的公式进行判断，当iflag＝1成立，对输入参数的所述经过离散反抛物拉冬变换的数据进行三次样条插值拉伸，形成x-t²数据文件；当iflag＝1不成立，对输入参数的经过反抛物拉冬变换的数据进行三次样条反插值拉伸，形成x-t数据文件。形成反t²拉伸后的地震数据。It further preferably includes an inverse ^t2 stretching process: read the data that has undergone discrete inverse parabolic Radon transformation, input relevant parameters, and judge by the formula of iflag=1, when iflag=1 is established, the input parameter The data through discrete anti-parabolic Radon transformation is subjected to cubic spline interpolation stretching to form an xt ² data file; when iflag=1 is not established, perform cubic spline inverse interpolation on the data of the input parameters after anti-parabolic Radon transformation Extend to form xt data file. Form the seismic data after inverse ^t2 stretching.

进一步优选包括对抛物拉冬变换进行验证：Further preferably includes verifying the parabolic Radon transformation:

首先定义抛物拉冬正变换为：First define the parabolic Radon forward transformation as:

$m m ((q q,, τ τ)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} d d ((x x,, τ τ + + q q {x x}^{22})) dx dx - - - - - - ((2.2.1 2.2.1))$

反变换为：The inverse transform is:

$d d ((x x,, t t)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} m m ((q q,, τ τ + + q q {x x}^{22})) dq dq - - - - - - ((2.2.2 2.2.2))$

由于抛物拉冬变换中t和τ是线性关系，所以可以在频率域中处理，设D(x，ω)和M(q，ω)分别为d(x，t)和m(q，τ)的傅立叶变换，则有变换对：Since t and τ are linear in the parabolic Radon transformation, it can be processed in the frequency domain. Let D(x, ω) and M(q, ω) be d(x, t) and m(q, τ) respectively The Fourier transform of , then there are transform pairs:

$M m ((q q,, ω ω)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} d d ((x x,, t t)) exp exp ((- - iωt iωt)) exp exp ((iωq iωq {x x}^{22})) dxdt wxya$

即频率域正变换为：That is, the forward transformation in the frequency domain is:

$M m ((q q,, ω ω)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} D D. ((x x,, ω ω)) exp exp ((iωq iωq {x x}^{22})) dx dx - - - - - - ((2.2.3 2.2.3))$

而and

$D D. ((x x,, ω ω)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} m m ((q q,, τ τ)) exp exp ((- - iωτ iωτ)) exp exp ((- - iωq iωq {x x}^{22})) dqdτ dqdτ$

即频率域反变换为：That is, the frequency domain inverse transform is:

$D D. ((x x,, ω ω)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} M m ((q q,, ω ω)) exp exp ((- - iωq iωq {x x}^{22})) dq dq - - - - - - ((2.2.4 2.2.4))$

将(2.2.4)代入(2.2.3)得到Substitute (2.2.4) into (2.2.3) to get

${M m}^{' '} ((q q,, ω ω)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} M m (({q q}^{' '},, ω ω)) exp exp ((iω iω ((q q - - {q q}^{' '})) {x x}^{22} d d {q q}^{' '} dx dx$

＝M(q，ω)*σ(q，ω) (2.2.5)=M(q, ω)*σ(q, ω) (2.2.5)

褶积因子为：The convolution factor is:

$σ σ ((q q,, ω ω)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} exp exp ((iωq iωq {x x}^{22})) dx dx$

$= = [[11 + + isign issign ((q q))]] \sqrt{\frac{π π}{22 ω ω | | q q | |}} - - - - - - ((2.2.6 2.2.6))$

当ωq＝0时σ→∞，为了得到更精确的抛物拉冬正变换，现对q方向进行反滤波，(2.2.5)对q方向作傅立叶变换后为When ωq=0, σ→∞, in order to obtain a more accurate parabolic Radon forward transformation, inverse filtering is now performed on the q direction, and (2.2.5) after Fourier transform on the q direction is

M′(k_q，ω)＝M(k_q，ω)σ(k_q，ω) (2.2.7a)M′(k _q ,ω)=M(k _q ,ω)σ(k _q ,ω) (2.2.7a)

或or

$M m (({k k}_{q q},, ω ω)) = = \frac{{M m}^{' '} (({k k}_{q q},, ω ω))}{σ σ (({k k}_{q q},, ω ω))}$

$= = \frac{\sqrt{ω ω} {M m}^{' '} (({k k}_{q q},, ω ω))}{{σ σ}^{' '} (({k k}_{q q}))} - - - - - - ((2.2.7 2.2.7 b b))$

其中in

${σ σ}^{' '} (({k k}_{q q})) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} [[11 + + isign issign ((q q))]] \sqrt{\frac{π π}{22 | | q q | |}} exp exp ((- - i i {k k}_{q q} q q)) dq dq$

$= = [[11 + + sign sign (({k k}_{q q}))]] \frac{π π}{\sqrt{| | {k k}_{q q} | |}}$

容易验证q变量对应的傅立叶变换的变量k_q为正，因此It is easy to verify that the variable k _q of the Fourier transform corresponding to the q variable is positive, so

${σ σ}^{' '} (({k k}_{q q})) = = \frac{22 π π}{\sqrt{{k k}_{q q}}}$

则得then have

$M m (({k k}_{q q},, ω ω)) = = \frac{11}{22 π π} \sqrt{ω ω {k k}_{q q}} {M m}^{' '} (({k k}_{q q},, ω ω)) - - - - - - ((2.2.8 2.2.8))$

上式显示在q方向做了反褶积，因此提高了抛物拉冬变换域的分辨率。The above formula shows that deconvolution is done in the q direction, thus improving the resolution of the parabolic Radon transform domain.

进一步优选包括对离散抛物拉冬变换进行验证，并形成离散抛物拉冬变换公式：Further preferably comprising verifying the discrete parabolic Radon transform, and forming a discrete parabolic Radon transform formula:

如果d(x，t)是偏移距-时间剖面，m(q，τ)为我们要寻找的速度-时间剖面，变量x，v，t，τ分别表示偏移距、速度、时间和截距时间。则从偏移距-时间剖面到速度-时间剖面的变换为：If d(x, t) is the offset-time profile, m(q, τ) is the velocity-time profile we are looking for, and the variables x, v, t, τ represent offset, velocity, time and intercept respectively from time. Then the transformation from the offset-time profile to the velocity-time profile is:

$m m ((v v,, τ τ)) = = {Σ Σ}_{i i = = 11}^{n no} d d (({x x}_{i i},, t t = = \sqrt{{τ τ}^{22} + + {x x}_{i i}^{22} / / {v v}^{22}})) - - - - - - ((2.2.9 2.2.9))$

n为偏移距个数，从速度-时间剖面到偏移距-时间剖面的变换为：n is the number of offsets, the transformation from velocity-time profile to offset-time profile is:

$d d ((v v,, t t)) = = {Σ Σ}_{j j = = 11}^{m m} m m (({v v}_{j j},, τ τ = = \sqrt{{t t}^{22} + + {x x}^{22} / / {v v}_{j j}^{22}})) - - - - - - ((2.2.10 2.2.10))$

m为速度采样点数。m is the number of velocity sampling points.

假设q＝1/v²，t′＝t²，τ′＝τ²则(2.2.9)和(2.2.10)变为离散抛物拉冬变换对Assuming q=1/v ² , t′=t ² , τ′=τ ² then (2.2.9) and (2.2.10) become discrete parabolic Radon transformation pairs

$m m (({q q}_{j j},, τ τ)) = = {Σ Σ}_{i i = = 11}^{n no} d d (({x x}_{i i},, t t = = τ τ + + {q q}_{j j} {x x}_{i i}^{22})),, j j = = 11,, m m - - - - - - ((2.2.11 2.2.11))$

$d d (({x x}_{i i},, t t)) = = {Σ Σ}_{j j = = 11}^{m m} m m (({q q}_{j j},, τ τ = = t t - - {q q}_{j j} {x x}_{i i}^{22})),, i i = = 11,, n no - - - - - - ((2.2.12 2.2.12))$

经傅立叶变换后为：After Fourier transform, it is:

$M m (({q q}_{j j},, ω ω)) = = {Σ Σ}_{i i = = 11}^{n no} D D. (({x x}_{i i},, ω ω)) exp exp ((iω iω {q q}_{j j} {x x}_{i i}^{22})),, j j = = 11,, m m - - - - - - ((2.2.13 2.2.13))$

$D D. (({x x}_{i i},, ω ω)) = = {Σ Σ}_{i i = = 11}^{m m} M m (({q q}_{j j},, ω ω)) exp exp ((- - iω iω {q q}_{j j} {x x}_{i i}^{22})),, i i = = 11,, n no - - - - - - ((2.2.14 2.2.14))$

对每一个频率成分(2.2.13)和(2.2.14)式可以表示为如下的矩阵形式为：For each frequency component (2.2.13) and (2.2.14) can be expressed as the following matrix form:

M(q_j)＝L^HD(x_i) (2.2.15)M(q _j )＝L ^H D( _xi ) (2.2.15)

D(x_i)＝LM(q_j) (2.2.16)D(x _i )＝LM(q _j ) (2.2.16)

其中 $L = \exp (- iω q_{j} x_{i}^{2});$ i＝1，n；j＝1，min $L = \exp (- iω q_{j} x_{i}^{2});$ i=1,n;j=1,m

写成矩阵形式为：Written in matrix form as:

$L L = = [\begin{matrix} exp exp ((- - iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{11}^{22})) & \cdot \cdot \cdot &Center Dot; \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{11}^{22})) \\ exp exp ((- - iω iω {q q}_{11} {x x}_{22}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot \cdot \cdot &Center Dot; \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{22}^{22})) \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot \cdot \\ \cdot \cdot & \cdot &Center Dot; & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & \cdot \cdot \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ exp exp ((- - iω iω {q q}_{11} {x x}_{n no}^{22})) & exp exp ((- - iω iω {q q}_{22} {x x}_{n no}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & exp exp ((- - iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}] - - - - - - ((2.2.17 2.2.17))$

$L^{H} = \exp (iω q_{j} x_{i}^{2});$ i＝1，n；j＝1，m $L^{h} = \exp (iω q_{j} x_{i}^{2});$ i=1,n; j=1,m

矩阵形式为：The matrix form is:

${L L}^{H h} = = [\begin{matrix} exp exp ((iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{11} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot \cdot \cdot \cdot & exp exp ((iω iω {q q}_{11} {x x}_{n no}^{22})) \\ exp exp ((iω iω {q q}_{22} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((iω iω {q q}_{22} {x x}_{n no}^{22})) \\ \cdot \cdot & \cdot \cdot & \cdot &Center Dot; \\ \cdot \cdot & \cdot \cdot & \cdot \cdot \cdot \cdot \cdot \cdot & \cdot &Center Dot; \\ \cdot \cdot & \cdot \cdot & \cdot \cdot \\ exp exp ((iω iω {q q}_{m m} {x x}_{11}^{22})) & exp exp ((iω iω {q q}_{m m} {x x}_{22}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}] - - - - - - ((2.2.18 2.2.18))$

由(2.2.15)和(2.2.16)式知，L和L^H是长方形矩阵，因此L和L^H并不是真正的互逆算子，这时我们就要利用广义逆来求解。From (2.2.15) and (2.2.16), we know that L and L ^H are rectangular matrices, so L and L ^H are not real reciprocal operators, and we need to use generalized inverse to solve them.

当n≥m，即方程(2.2.12)为超定方程时，m×m矩阵LHL是可逆的，因此L的广义逆矩阵为：When n≥m, that is, when equation (2.2.12) is an overdetermined equation, the m×m matrix LHL is invertible, so the generalized inverse matrix of L is:

L^-1＝(L^HL)^-1L^H (2.2.19)L ^-1 ＝(L ^H L) ^-1 L ^H (2.2.19)

则离散抛物拉冬正变换变为：Then the discrete parabolic Radon forward transformation becomes:

M＝(L^HL)^-1L^HD (2.2.20)M＝(L ^H L) ^-1 L ^H D (2.2.20)

当n＜m，即方程(2.2.12)为欠定方程时，n×n矩阵LL^H是可逆的，因此L的广义逆矩阵为：When n<m, that is, when the equation (2.2.12) is an underdetermined equation, the n×n matrix LL ^H is reversible, so the generalized inverse matrix of L is:

L^-1＝L^H(LL^H)-1 (2.2.21)L ^-1 ＝L ^H (LL ^H )-1 (2.2.21)

M＝L^H(LL^H)^-1D (2.2.22)M＝L ^H (LL ^H ) ^-1 D (2.2.22)

算子L^HL为：The operator L ^H L is:

${L L}^{H h} L L = = [\begin{matrix} exp exp ((iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{11} {x x}_{22}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((iω iω {q q}_{11} {x x}_{n no}^{22})) \\ exp exp ((iω iω {q q}_{22} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & exp exp ((iω iω {q q}_{22} {x x}_{n no}^{22})) \\ \cdot &Center Dot; & \cdot \cdot & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot \cdot & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot \cdot \\ exp exp ((iω iω {q q}_{m m} {x x}_{11}^{22})) & exp exp ((iω iω {q q}_{m m} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & exp exp ((iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}]$

$\times \times [\begin{matrix} exp exp ((- - iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{11}^{22})) & \cdot \cdot \cdot &Center Dot; \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{11}^{22})) \\ exp exp ((- - iω iω {q q}_{11} {x x}_{22}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot \cdot \cdot &Center Dot; \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{22}^{22})) \\ \cdot &Center Dot; & \cdot \cdot & \cdot \cdot \\ \cdot \cdot & \cdot \cdot & \cdot &Center Dot; \cdot \cdot \cdot \cdot & \cdot \cdot \\ \cdot &Center Dot; & \cdot \cdot & \cdot \cdot \\ exp exp ((- - iω iω {q q}_{11} {x x}_{n no}^{22})) & exp exp ((- - iω iω {q q}_{22} {x x}_{n no}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}]$

$= = [\begin{matrix} n no & {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{11} - - {q q}_{22})) {x x}_{i i}^{22}]] & \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; & {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{11} - - {q q}_{m m})) {x x}_{i i}^{22}]] \\ {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{22} - - {q q}_{11})) {x x}_{i i}^{22}]] & n no & \cdot \cdot \cdot &Center Dot; \cdot \cdot & {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{22} - - {q q}_{m m})) {x x}_{i i}^{22}]] \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{m m} - - {q q}_{11})) {x x}_{i i}^{22}]] & {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{m m} - - {q q}_{22})) {x x}_{i i}^{22}]] & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & n no \end{matrix}] - - - - - - ((2.2.23 2.2.23))$

同理in the same way

$L L {L L}^{H h} = = [\begin{matrix} exp exp ((- - iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{11}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{11}^{22})) \\ exp exp ((- - iω iω {q q}_{11} {x x}_{22}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{22}^{22})) \\ \cdot \cdot & \cdot \cdot & \cdot \cdot \\ \cdot \cdot & \cdot \cdot & \cdot \cdot \cdot \cdot \cdot \cdot & \cdot \cdot \\ \cdot \cdot & \cdot \cdot & \cdot \cdot \\ exp exp ((- - iω iω {q q}_{11} {x x}_{n no}^{22})) & exp exp ((- - iω iω {q q}_{22} {x x}_{n no}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}]$

$\times \times [\begin{matrix} exp exp ((iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{11} {x x}_{22}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((iω iω {q q}_{11} {x x}_{n no}^{22})) \\ exp exp ((iω iω {q q}_{22} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot \cdot \cdot \cdot & exp exp ((iω iω {q q}_{22} {x x}_{n no}^{22})) \\ \cdot \cdot & \cdot \cdot & \cdot \cdot \\ \cdot \cdot & \cdot \cdot & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & \cdot \cdot \\ \cdot &Center Dot; & \cdot \cdot & \cdot \cdot \\ exp exp ((iω iω {q q}_{m m} {x x}_{11}^{22})) & exp exp ((iω iω {q q}_{m m} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & exp exp ((iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}]$

$= = [\begin{matrix} m m & {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{22}^{22} - - {x x}_{11}^{22})) {q q}_{j j}]] & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{n no}^{22} - - {x x}_{11}^{22})) {q q}_{j j}]] \\ {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{11}^{22} - - {x x}_{22}^{22})) {q q}_{j j}]] & m m & \cdot \cdot \cdot &Center Dot; \cdot \cdot & {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{n no}^{22} - - {x x}_{22}^{22})) {q q}_{j j}]] \\ \cdot \cdot & \cdot \cdot & \cdot \cdot \\ \cdot &Center Dot; & \cdot \cdot & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{11}^{22} - - {x x}_{n no}^{22})) {q q}_{j j}]] & {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{22}^{22} - - {x x}_{n no}^{22})) {q q}_{j j}]] & \cdot \cdot \cdot \cdot \cdot \cdot & m m \end{matrix}] - - - - - - ((2.2.24 2.2.24))$

由于L^HL和LL^H有可能为奇异阵或接近奇异阵，特别是对于低频成分。因此，为了保证数值解的稳定，常要加一阻尼因子，使得解的形式变为：Because L ^H L and LL ^H may be singular or close to singular, especially for low frequency components. Therefore, in order to ensure the stability of the numerical solution, it is often necessary to add a damping factor, so that the form of the solution becomes:

M＝(L^HL+λ²I)^-1L^HD (2.2.25)M＝(L ^H L+λ ² I) ^-1 L ^H D (2.2.25)

M＝L^H(LL^H+λ²I)^-1D (2.2.26)M＝L ^H (LL ^H +λ ² I) ^-1 D (2.2.26)

λ为阻尼因子，取值一般为0.1～1。由此，(2.2.20)和(2.2.22)式的最小二乘问题变为(2.2.25)和(2.2.26)式的阻尼最小二乘问题。λ is the damping factor, and its value is generally 0.1~1. Thus, the least squares problem of (2.2.20) and (2.2.22) becomes the damped least squares problem of (2.2.25) and (2.2.26).

在方程(2.2.20)和(2.2.22)中要分别计算L^HL和LL^H，由(2.2.23)和(2.2.24)可知这两个矩阵具有特殊的结构，都是Hermit矩阵。且L^HL具有Toeplitz结构，因此又称为Toeplitz-Hermit矩阵。对于该矩阵可以用Levison递推算法求解。In equations (2.2.20) and (2.2.22), it is necessary to calculate L ^H L and LL ^H respectively. From (2.2.23) and (2.2.24), we can see that these two matrices have a special structure, both are Hermit matrices . And L ^H L has a Toeplitz structure, so it is also called a Toeplitz-Hermit matrix. The matrix can be solved by Levison recursive algorithm.

进一步优选包括对优选参数进行选择：Further preferably includes selecting the preferred parameters:

由信号分析理论可以知道，对于有限带宽B_f的任意函数x(t)，其采样率为Δt，则应该满足如下关系式It can be known from the signal analysis theory that for any function x(t) with a finite bandwidth B _f and its sampling rate Δt, the following relationship should be satisfied

ΔtB_f≤2π (2.2.27)ΔtB _f ≤ 2π (2.2.27)

对于离散抛物拉冬变换，其带宽为

把它代入(2.2.27)式中，则又对于频率为ω的信号，离散抛物拉冬变换中参数q的采样率应该满足下式For the discrete parabolic Radon transform, the bandwidth is

Substituting it into (2.2.27), then for a signal with frequency ω, the sampling rate of parameter q in the discrete parabolic Radon transform should satisfy the following formula

$Δq Δq \leq \leq \frac{22 π π}{ω ω (({x x}_{max max}^{22} - - {x x}_{min min}^{22}))} - - - - - - ((2.2.28 2.2.28))$

式中x_max和x_min分别为最大炮检距和最小炮检距，由此可知离散抛物拉冬变换中参数q的临界采样率Δq_c为where x _max and x _min are the maximum and minimum offsets respectively, from which it can be seen that the critical sampling rate Δq _c of the parameter q in the discrete parabolic Radon transform is

$Δ Δ {q q}_{c c} = = \frac{22 π π}{ω ω (({x x}_{max max}^{22} - - {x x}_{min min}^{22}))} - - - - - - ((2.2.29 2.2.29))$

设原始数据中信号的最高频率为ω_max，欲使所有的有效频率成分均成立，则将上式的ω替换为有ω_max，则有Assuming that the highest frequency of the signal in the original data is ω _max , and to make all effective frequency components valid, replace ω in the above formula with ω _max , then we have

$Δ Δ {q q}_{c c} = = \frac{22 π π}{{ω ω}_{max max} (({x x}_{max max}^{22} - - {x x}_{min min}^{22}))} - - - - - - ((2.2.30 2.2.30))$

由于在离散抛物拉冬变换过程中出现假频，使得变换的质量降低。这里通过选择适当的q的扫描范围来减小假频产生。对于离散抛物拉冬变换设其真实的同相轴对应的曲率参数为q₀，则抗假频条件为：The quality of the transform is degraded due to aliasing during the discrete parabolic Radon transform. Here by selecting the appropriate scanning range of q to reduce false frequency generation. For the discrete parabolic Radon transform, if the curvature parameter corresponding to its real event is q ₀ , then the anti-aliasing condition is:

$| | q q - - {q q}_{00} | | \leq \leq \frac{π π}{ω ω | | {x x}_{max max} | | Δx Δx} - - - - - - ((2.2.31 2.2.31))$

考虑到所有有效波的成分，则上式可以写成Considering all effective wave components, the above formula can be written as

$Δx Δx \leq \leq \frac{π π}{{ω ω}_{max max} | | {x x}_{max max} | | (({q q}_{max max} - - {q q}_{min min}))} - - - - - - ((2.2.32 2.2.32))$

由此可以得到曲率参数q的扫描范围满足的条件为：From this, it can be obtained that the scanning range of the curvature parameter q satisfies the following conditions:

$(({q q}_{max max} - - {q q}_{min min})) \leq \leq \frac{π π}{{ω ω}_{max max} | | {x x}_{max max} | | Δx Δx} - - - - - - ((2.2.33 2.2.33))$

因此离散抛物拉冬变换的实现过程中，参数的选择不依赖于信号的频率，而只与原始数据的空间采样率、最大偏移距和信号的最高频率有关。Therefore, in the realization of discrete parabolic Radon transform, the selection of parameters does not depend on the frequency of the signal, but only on the spatial sampling rate of the original data, the maximum offset and the highest frequency of the signal.

本发明的有益效果在于，利用离散抛物拉冬变换，高分辨率并快速的实现变偏移距VSP的上行波和下行波的分离，同时由于使用的时离散抛物拉冬变换，克服了现有技术中使用抛物拉冬变换计算量大、对设备要求高、计算速度低的种种缺陷。提高地震分辨率，实现波场的完全分离，还可以提高计算效率，为构造勘探、岩性勘探和油气开发提供更准确快捷的数据支持。同时该方法还具有实现变偏移距VSP资料在地下空间的准确归位，同时既可以实现叠前时间偏移也可以实现叠前深度偏移。The beneficial effect of the present invention is that, using the discrete parabolic Radon transform, the separation of the upgoing wave and the downgoing wave of the VSP with variable offset distance can be realized quickly with high resolution, and at the same time, due to the use of the discrete parabolic Radon transform, it overcomes the existing The use of parabolic Radon transformation in the technology has various defects such as large amount of calculation, high requirements for equipment, and low calculation speed. Improving the seismic resolution and realizing the complete separation of the wave field can also improve the calculation efficiency and provide more accurate and fast data support for structural exploration, lithology exploration and oil and gas development. At the same time, this method also has the ability to realize the accurate homing of variable offset distance VSP data in the underground space, and can realize both pre-stack time migration and pre-stack depth migration.

附图说明Description of drawings

图1是变偏移距VSP观察方法示意图；Figure 1 is a schematic diagram of the VSP observation method with variable offset distance;

图2是某油田原始变偏移距VSP资料；Figure 2 is the original VSP data of an oilfield;

图3是变偏移距VSP波场分离方法流程图；Fig. 3 is a flow chart of the VSP wave field separation method with variable offset distance;

图4是离散正抛物拉冬变换流程图；Fig. 4 is the flow chart of discrete positive parabolic Radon transformation;

图5是离散反抛物拉冬变换流程图；Fig. 5 is the flowchart of discrete inverse parabolic Radon transformation;

图6是t²拉伸处理和反t²拉伸处理流程图；Fig. 6 is a flowchart of ^t2 stretching process and reverse ^t2 stretching process;

图7是有多次下行波及转换SV波的离散抛物拉冬变换波场分离图；Fig. 7 is the discrete parabolic Radon transform wavefield separation diagram with multiple downgoing waves and converted SV waves;

图8是某油田变偏移距VSP实际资料的离散抛物拉冬波场分离图。Fig. 8 is the separation diagram of the discrete parabolic Radon wave field of the variable offset VSP actual data in an oilfield.

具体实施方式Detailed ways

有效的分离上行波和下行波是VSP资料处理的一项基本任务，只有正确的识别、分离和提取上行波和下行波，才能充分利用VSP的有用信息，补充地面勘探的不足，为构造勘探、岩性勘探和油气开发服务。分离VSP记录的上行波和下行波主要依据是二者的视速度不同。Effective separation of up-going waves and down-going waves is a basic task of VSP data processing. Only by correctly identifying, separating and extracting up-going waves and down-going waves can the useful information of VSP be fully utilized to supplement the insufficiency of ground exploration and contribute to structural exploration, Lithology exploration and oil and gas development services. The main basis for separating the upgoing and downgoing waves recorded by the VSP is the difference in their apparent velocities.

离散拉冬变换在地震勘探中早期以线性拉冬变换为主的情况下进行的一种改进，抛物拉冬变化对压制线性干扰波、多次波、分离纵横波都有明显的效果。但因为地震资料的反射波同相轴多为双曲线，它在变换域内变为椭圆，尽管速度的差异使不同的波场有所分离，但并未完全分离。双曲拉冬变换是将地震波同相轴在变换域聚焦，提高了分辨率，能成功的实现波场的完全分离。但是直接进行双曲拉冬变换是时变的，计算量相当大。而离散抛物拉冬变换不仅具有时不变性，同时还可使地震同相轴在变换域内聚焦，这样不仅可以提高地震分辨率，实现波场的完全分离，还可以提高计算效率。Discrete Radon transform is an improvement in the early stage of seismic exploration when linear Radon transform is the main method. Parabolic Radon transform has obvious effects on suppressing linear interference waves, multiple waves, and separating longitudinal and transverse waves. But because the reflected wave event of seismic data is mostly a hyperbola, it becomes an ellipse in the transform domain. Although the difference in velocity makes the different wave fields separate, they are not completely separated. The hyperbolic Radon transform focuses the seismic wave event in the transform domain, improves the resolution, and successfully realizes the complete separation of the wave field. However, the direct hyperbolic Radon transform is time-varying, and the amount of calculation is quite large. The discrete parabolic Radon transform not only has time invariance, but also can focus the seismic event in the transform domain, which can not only improve the seismic resolution, realize the complete separation of the wave field, but also improve the computational efficiency.

一、VSP波场分离方法：1. VSP wave field separation method:

变偏移距VSP波场分离方法包括：在计算机中输入采集到的变偏移距VSP的原始地震波场数据；对所述的原始地震波场数据进行t²拉伸处理，形成t²拉伸处理数据；对所述的t²拉伸处理数据进行对离散正抛物拉冬变换，形成τ-q域上下行波分数据；对所述τ-q域上下行波分数据进行离散反抛物拉冬变换，形成x-t²域数据；对所述x-t²域数据进行反t²拉伸处理，完成变偏移距VSP上行波和下行波的分离。VSP波场分离方法如图3所示，反映出变偏移距VSP数据通过离散正抛物线拉冬变换与离散反抛物拉冬变换以及t²拉伸和t²反拉伸程序，进行变偏移距VSP的上、下行波分离；The variable offset VSP wavefield separation method includes: inputting the original seismic wavefield data of the variable offset VSP collected in the computer; performing ^t2 stretch processing on the original seismic wavefield data to form a ^t2 stretch process Data; carry out discrete forward parabolic Radon transformation to described ^t2 stretch processing data, form τ-q domain uplink and downlink wave division data; Described τ-q domain uplink and downlink wave division data carry out discrete inverse parabolic Radon transformation to form xt ² domain data; reverse t ² stretch processing is performed on the xt ² domain data to complete the separation of upgoing waves and downgoing waves of variable offset VSP. The VSP wave field separation method is shown in Figure 3, which reflects that the VSP data with variable offsets are subjected to variable offsets through discrete forward parabolic Radon transform and discrete inverse parabolic Radon transform, as well as ^t2 stretching and ^t2 reverse stretching procedures. Separation of up-going and down-going waves from VSP;

t²拉伸处理过程：对原始地震波场数据进行读取，输入相关参数，通过iflag＝1的公式进行判断，当iflag＝1成立，对输入参数的原始地震波场数据进行三次样条插值拉伸，形成x-t²数据文件；当iflag＝1不成立，对输入参数的原始地震波场数据进行三次样条反插值拉伸，形成x-t数据文件。形成t²拉伸后的地震数据。其中t²拉伸程序如图6所示，主要功能为对原始地震记录进行t²拉伸，处理过程中运用了三次样条进行插值；t ² Stretching process: read the original seismic wavefield data, input relevant parameters, judge through the formula of iflag=1, when iflag=1 is established, perform cubic spline interpolation stretching on the original seismic wavefield data of the input parameters , to form an xt ² data file; when iflag=1 is not established, perform cubic spline inverse interpolation stretching on the original seismic wave field data of the input parameters, and form an xt data file. Form the seismic data after ^t2 stretching. Among them, the ^t2 stretching program is shown in Figure 6, the main function is to carry out ^t2 stretching on the original seismic records, and the cubic spline is used for interpolation in the process of processing;

离散正抛物拉冬变换的过程：读取t²拉伸后的地震数据，输入参数，计算q的扫描范围rq，对计算后的扫描范围rq进行判断，rq不满足要求返回到输入参数步骤，rq满足要求进行FFT处理，并对处理后的数据进行一维滤波；对滤波后的数据进行判断，当n＜m时利用M＝L^H(LL^H+λ²1)^-1D对数据进行处理，当n≥m时利用M＝(L^HL+λ²I)^-1L^HD对数据进行处理，形成处理后的数据；对所述处理后的数据进行二维褶积，并对二维褶积后的数据进行一维滤波并进行IFFT处理，形成τ-q域数据文件，即形成τ-q域上下行波分离。其中离散正抛物拉冬变换过程如图4所示，主要功能是完成正变换，并且形成τ-q域的数据文件；The process of discrete forward parabolic Radon transformation: read the seismic data stretched by ^t2 , input parameters, calculate the scanning range rq of q, and judge the calculated scanning range rq, if rq does not meet the requirements, return to the input parameter step, rq meets the requirements to perform FFT processing, and perform one-dimensional filtering on the processed data; judge the filtered data, and use M=L ^H (LL ^H +λ ² 1) ^-1 D to perform data processing when n<m Processing, when n≥m, use M=(L ^H L+λ ² I) ^-1 ^L HD to process the data to form processed data; perform two-dimensional convolution on the processed data, and The data after two-dimensional convolution is subjected to one-dimensional filtering and IFFT processing to form a τ-q domain data file, that is, to form the separation of uplink and downlink waves in the τ-q domain. The discrete forward parabolic Radon transformation process is shown in Figure 4, the main function is to complete the forward transformation and form a data file in the τ-q domain;

离散反抛物拉冬变换的过程：读取τ-q域数据文件，读取输入的参数文件，对数据进行FFT处理，对处理后的数据进行第一次频率域滤波，然后对第一次频率域滤波后的数据进行D＝LM处理，对处理后的数据进行第二次频率域滤波，对第二次频率域滤波后的数据进行IFFT处理，形成x-t²域数据文件。其中离散反抛物拉冬变换如图5所示，主要功能是利用τ-q域数据文件重建x-t²域数据文件；The process of discrete inverse parabolic Radon transformation: read the τ-q domain data file, read the input parameter file, perform FFT processing on the data, perform the first frequency domain filtering on the processed data, and then filter the first frequency D=LM processing is performed on the data after domain filtering, a second frequency domain filtering is performed on the processed data, and IFFT processing is performed on the data after the second frequency domain filtering to form an xt ² domain data file. Among them, the discrete anti-parabolic Radon transform is shown in Figure 5, and its main function is to use the τ-q domain data file to reconstruct the xt ² domain data file;

反t²拉伸处理过程：对经过离散反抛物拉冬变换的数据进行读取，输入相关参数，通过iflag＝1的公式进行判断，当iflag＝1成立，对输入参数的经过离散反抛物拉冬变换的数据进行三次样条插值拉伸，形成x-t²数据文件；当iflag＝1不成立，对输入参数的经过反抛物拉冬变换的数据进行三次样条反插值拉伸，形成x-t数据文件。形成反t²拉伸后的地震数据。其中t²反拉伸程序如图6所示，主要功能为对经过离线反抛物拉冬变换所形成的x-t²域的数据文件进行反t²拉伸，处理过程中运用了三次样条进行插值。Inverse ^t2 stretch processing process: read the data that has undergone discrete inverse parabolic pull, input relevant parameters, and judge through the formula of iflag=1. When iflag=1 is established, the discrete inverse parabolic pull Perform cubic spline interpolation stretching on the winter transformed data to form an xt ² data file; when iflag=1 is not established, perform cubic spline inverse interpolation stretching on the input parameters after inverse parabolic Radon transformation to form an xt data file. Form the seismic data after inverse ^t2 stretching. Among them, the t ² inverse stretching program is shown in Figure 6. Its main function is to inverse t ² stretch the data files in the xt ² domain formed by off-line inverse parabolic Radon transformation. During the processing, cubic splines are used for interpolation .

二、常规抛物拉冬变换：2. Conventional parabolic Radon transformation:

在前人提出的理论及方法的基础上，对变偏移距VSP的波场特征、抛物拉冬变换的基本方法理论和算法进行了分析。研究中对大量的变偏移距VSP模型进行了试算，证实了算法的合理性和计算结果的正确性。通过对实际资料进行处理，证明该算法具有实用性和有效性。On the basis of the theories and methods proposed by the predecessors, the wave field characteristics of the variable offset VSP, the basic method theory and algorithm of the parabolic Radon transformation are analyzed. In the study, a large number of VSP models with variable offset distances were tested, which proved the rationality of the algorithm and the correctness of the calculation results. The practicality and effectiveness of the algorithm are proved by processing the actual data.

反变换为：The inverse transform is:

而and

即频率域反变换为：That is, the frequency domain inverse transform is:

将(2.2.4)代入(2.2.3)得到Substitute (2.2.4) into (2.2.3) to get

＝M(q，ω)*σ(q，ω) (2.2.5)= M(q, ω)*σ(q, ω) (2.2.5)

褶积因子为：The convolution factor is:

或or

其中in

${σ σ}^{' '} (({k k}_{q q})) = = \frac{22 π π}{\sqrt{{k k}_{q q}}}$

则得then have

三、本申请创造性应用的离散抛物拉冬变换：3. The discrete parabolic Radon transform of the creative application of this application:

如果d(x，t)是偏移距-时间剖面，m(q，τ)为我们要寻找的速度-时间剖面，变量x，v，t，τ分别表示偏移距、速度、时间和截距时间。则从偏移距-时间剖面到速度-时间剖面的变换为：If d(x, t) is the offset-time profile, m(q, τ) is the velocity-time profile we are looking for, and the variables x, v, t, τ represent offset, velocity, time and intercept respectively from time. Then the transformation from offset-time profile to velocity-time profile is:

m为速度采样点数。m is the number of velocity sampling points.

经傅立叶变换后为：After Fourier transform, it is:

M(q_j)＝L^HD(x_i) (2.2.15)M(q _j )＝L ^H D( _xi ) (2.2.15)

D(x_i)＝LM(q_j) (2.2.16)D(x _i )＝LM(q _j ) (2.2.16)

其中 $L = \exp (- iω q_{j} x_{i}^{2});$ i＝1，n；j＝1，min $L = \exp (- iω q_{j} x_{i}^{2});$ i=1,n; j=1,m

写成矩阵形式为：Written in matrix form as:

$L L = = [\begin{matrix} exp exp ((- - iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{11}^{22})) & \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; & exp exp ((- - iω iω {q q}_{m m} {x x}_{11}^{22})) \\ exp exp ((- - iω iω {q q}_{11} {x x}_{22}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot \cdot \cdot &Center Dot; \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{22}^{22})) \\ \cdot \cdot & \cdot \cdot & \cdot &Center Dot; \\ \cdot \cdot & \cdot \cdot & \cdot \cdot \cdot &Center Dot; \cdot \cdot & \cdot \cdot \\ \cdot &Center Dot; & \cdot \cdot & \cdot &Center Dot; \\ exp exp ((- - iω iω {q q}_{11} {x x}_{n no}^{22})) & exp exp ((- - iω iω {q q}_{22} {x x}_{n no}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & exp exp ((- - iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}] - - - - - - ((2.2.17 2.2.17))$

$L^{H} = \exp (iω q_{j} x_{i}^{2});$ i＝1，n；j＝1，m $L^{h} = \exp (iω q_{j} x_{i}^{2});$ i=1,n;j=1,m

矩阵形式为：The matrix form is:

${L L}^{H h} = = [\begin{matrix} exp exp ((iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{11} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & exp exp ((iω iω {q q}_{11} {x x}_{n no}^{22})) \\ exp exp ((iω iω {q q}_{22} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & exp exp ((iω iω {q q}_{22} {x x}_{n no}^{22})) \\ \cdot \cdot & \cdot \cdot & \cdot \cdot \\ \cdot \cdot & \cdot \cdot & \cdot \cdot \cdot \cdot \cdot \cdot & \cdot \cdot \\ \cdot \cdot & \cdot \cdot & \cdot \cdot \\ exp exp ((iω iω {q q}_{m m} {x x}_{11}^{22})) & exp exp ((iω iω {q q}_{m m} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot \cdot \cdot \cdot & exp exp ((iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}] - - - - - - ((2.2.18 2.2.18))$

当n≥m，即方程(2.2.12)为超定方程时，m×m矩阵L^HL是可逆的，因此L的广义逆矩阵为：When n≥m, that is, when the equation (2.2.12) is an overdetermined equation, the m×m matrix L ^H L is reversible, so the generalized inverse matrix of L is:

L^-1＝(L^HL)^-1L^H (2.2.19)L ^-1 ＝(L ^H L) ^-1 L ^H (2.2.19)

M＝(L^HL)^-1L^HD (2.2.20)M＝(L ^H L) ^-1 L ^H D (2.2.20)

L^-1＝L^H(LL^H)^-1 (2.2.21)L ^-1 = L ^H (LL ^H ) ^-1 (2.2.21)

M＝L^H(LL^H)^-1D (2.2.22)M＝L ^H (LL ^H ) ^-1 D (2.2.22)

算子L^HL为：The operator L ^H L is:

${L L}^{H h} L L = = [\begin{matrix} exp exp ((iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{11} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & exp exp ((iω iω {q q}_{11} {x x}_{n no}^{22})) \\ exp exp ((iω iω {q q}_{22} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot \cdot \cdot \cdot & exp exp ((iω iω {q q}_{22} {x x}_{n no}^{22})) \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot \cdot \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot \cdot & \cdot &Center Dot; & \cdot \cdot \\ exp exp ((iω iω {q q}_{m m} {x x}_{11}^{22})) & exp exp ((iω iω {q q}_{m m} {x x}_{22}^{22})) & \cdot \cdot \cdot &Center Dot; \cdot &Center Dot; & exp exp ((iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}]$

$\times \times [\begin{matrix} exp exp ((- - iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{11}^{22})) & \cdot \cdot \cdot &Center Dot; \cdot &Center Dot; & exp exp ((- - iω iω {q q}_{m m} {x x}_{11}^{22})) \\ exp exp ((- - iω iω {q q}_{11} {x x}_{22}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & exp exp ((- - iω iω {q q}_{m m} {x x}_{22}^{22})) \\ \cdot \cdot & \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot \cdot & \cdot \cdot \cdot &Center Dot; \cdot &Center Dot; & \cdot \cdot \\ \cdot &Center Dot; & \cdot \cdot & \cdot \cdot \\ exp exp ((- - iω iω {q q}_{11} {x x}_{n no}^{22})) & exp exp ((- - iω iω {q q}_{22} {x x}_{n no}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}]$

$= = [\begin{matrix} n no & {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{11} - - {q q}_{22})) {x x}_{i i}^{22}]] & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{11} - - {q q}_{m m})) {x x}_{i i}^{22}]] \\ {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{22} - - {q q}_{11})) {x x}_{i i}^{22}]] & n no & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{22} - - {q q}_{m m})) {x x}_{i i}^{22}]] \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & \cdot \cdot \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{m m} - - {q q}_{11})) {x x}_{i i}^{22}]] & {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{m m} - - {q q}_{22})) {x x}_{i i}^{22}]] & \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; & n no \end{matrix}] - - - - - - ((2.2.23 2.2.23))$

同理in the same way

$L L {L L}^{H h} = = [\begin{matrix} exp exp ((- - iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{11}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & exp exp ((- - iω iω {q q}_{m m} {x x}_{11}^{22})) \\ exp exp ((- - iω iω {q q}_{11} {x x}_{22}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot \cdot \cdot &Center Dot; \cdot &Center Dot; & exp exp ((- - iω iω {q q}_{m m} {x x}_{22}^{22})) \\ \cdot \cdot & \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & \cdot &Center Dot; \\ \cdot \cdot & \cdot &Center Dot; & \cdot &Center Dot; \\ exp exp ((- - iω iω {q q}_{11} {x x}_{n no}^{22})) & exp exp ((- - iω iω {q q}_{22} {x x}_{n no}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}]$

$\times \times [\begin{matrix} exp exp ((iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{11} {x x}_{22}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((iω iω {q q}_{11} {x x}_{n no}^{22})) \\ exp exp ((iω iω {q q}_{22} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((iω iω {q q}_{22} {x x}_{n no}^{22})) \\ \cdot \cdot & \cdot &Center Dot; & \cdot \cdot \\ \cdot &Center Dot; & \cdot \cdot & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot \cdot & \cdot \cdot \\ exp exp ((iω iω {q q}_{m m} {x x}_{11}^{22})) & exp exp ((iω iω {q q}_{m m} {x x}_{22}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}]$

$= = [\begin{matrix} m m & {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{22}^{22} - - {x x}_{11}^{22})) {q q}_{j j}]] & \cdot &Center Dot; \cdot \cdot \cdot \cdot & {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{n no}^{22} - - {x x}_{11}^{22})) {q q}_{j j}]] \\ {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{11}^{22} - - {x x}_{22}^{22})) {q q}_{j j}]] & m m & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{n no}^{22} - - {x x}_{22}^{22})) {q q}_{j j}]] \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot \cdot \\ \cdot \cdot & \cdot \cdot & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & \cdot \cdot \\ \cdot &Center Dot; & \cdot \cdot & \cdot &Center Dot; \\ {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{11}^{22} - - {x x}_{n no}^{22})) {q q}_{j j}]] & {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{22}^{22} - - {x x}_{n no}^{22})) {q q}_{j j}]] & \cdot \cdot \cdot &Center Dot; \cdot \cdot & m m \end{matrix}] - - - - - - ((2.2.24 2.2.24))$

M＝(L^HL+λ²I)^-1L^HD (2.2.25)M＝(L ^H L+λ ² I) ^-1 L ^H D (2.2.25)

M＝L^H(LL^H+λ²I)^-1D (2.2.26)M＝L ^H (LL ^H +λ ² I) ^-1 D (2.2.26)

四、有关参数选择4. Relevant parameter selection

在地震数据处理中，参数的选择是很重要的一个环节，好的参数的选择对方法的应用和效果都是很重要的，离散抛物拉冬变换也不例外。In seismic data processing, the selection of parameters is a very important link. A good selection of parameters is very important to the application and effect of the method, and the discrete parabolic Radon transform is no exception.

ΔtB_f≤2π (2.2.27)ΔtB _f ≤ 2π (2.2.27)

对于离散抛物拉冬变换，其带宽为

五、具体计算实例5. Specific calculation examples

变偏移距VSP方法中波场复杂，在地震记录中不仅含有直达波、一次下行波、多次下行波，而且含有转换波。尤其对于靠近检波点地层其反射波与直达波逐渐重合，给波场的完全分离带来了很大的难度。因此，我们首先选择模型的简单理论记录进行试算，完成简单的上下行波分离，与原始模型记录进行对比，以验证程序的准确性及可行性。然后，对复杂波场进行分离，用它进行反变换所得到的数据与原始模型记录进行对比，以进一步检验程序的正确性和实用性。The wave field in the variable offset VSP method is complex, and the seismic records not only contain direct waves, primary downgoing waves, multiple downgoing waves, but also converted waves. Especially for the formation near the receiver point, the reflected wave and the direct wave gradually overlap, which brings great difficulty to the complete separation of the wave field. Therefore, we first select the simple theoretical records of the model for trial calculation, complete the simple separation of uplink and downgoing waves, and compare them with the original model records to verify the accuracy and feasibility of the program. Then, the complex wave field is separated, and the data obtained by inverse transformation are compared with the original model records to further test the correctness and practicability of the program.

(一)、简单理论计算(1) Simple theoretical calculation

在一般情况下，变偏移距VSP数据中都会含有一次下行波和上行波，因此首先我们对几个只包含有一次下行波和上行波，利用不同采样参数的原始记录剖面进行上下行波分离，以检验所实现方法的正确性。In general, variable offset VSP data will contain primary downgoing waves and upgoing waves, so first we use the original record profiles with different sampling parameters to separate the upgoing and downgoing waves. , to verify the correctness of the implemented method.

选用的模型采用不同的检波点深度、炮间距、偏移距，并在进行射线追踪合成记录时分别加入了一次上下行波，二次下行波和转换波成分，通过试算表明，对各种模型的记录所进行的上下行波分离结果与理论结果能够比较精确的吻合，并且适用于存在速度差异的波场分离。The selected model adopts different depths of receiver points, shot distances, and offsets, and adds the first up-down wave, the second down-going wave, and the converted wave components when performing ray tracing synthesis records. The up-going wave separation results recorded by the model are in good agreement with the theoretical results, and are suitable for wave field separation with velocity differences.

(二)、复杂理论计算(2) Complicated theoretical calculations

然而我们知道，在变偏移距VSP方法所记录的数据中并不仅仅只含有单一的上下行波，或者仅仅含有纵波和转换波，而是含有多次上行纵波、多次下行波、转换波等。因此我们应该关心含有复杂波场记录的波场分离。However, we know that the data recorded by the variable offset VSP method does not only contain a single up-and-down wave, or only a longitudinal wave and a converted wave, but contains multiple up-going longitudinal waves, multiple down-going waves, and converted waves. wait. We should therefore be concerned with wavefield separation for records containing complex wavefields.

下列模型在射线追踪合成记录时同时加入了二次下行波和转换SV波。图7为有多次下行波及转换SV波的离线抛物拉冬变换波场分离。图7(a)为射线追踪合成记录，其中包括一次下行波、二次下行波、两组上行纵波和两组上行转换SV波。对其进行阵变换后得到了τ-q域结果如图7(b)，从图7中可以完全分辨出这几组波，经过切除处理得到了上行纵波如图7(c)和上行转换SV波记录如图7(d)。通过以上各个模型的试算说明该程序的计算是正确的，所得到的结果是可靠的，能达到预期的波场分离效果。The following models incorporate both secondary downgoing and converted SV waves in ray-traced synthetic recordings. Figure 7 shows the separation of off-line parabolic Radon transformed waves with multiple downgoing waves and transformed SV waves. Figure 7(a) is the ray tracing synthetic record, which includes primary downgoing waves, secondary downgoing waves, two groups of upgoing longitudinal waves and two groups of upgoing converted SV waves. After the matrix transformation, the result in the τ-q domain is obtained as shown in Figure 7(b). From Figure 7, these groups of waves can be completely distinguished. The wave recording is shown in Fig. 7(d). The trial calculation of the above models shows that the calculation of the program is correct, the obtained results are reliable, and the expected wave field separation effect can be achieved.

(三)、实际资料计算(3) Actual data calculation

对于任何一种地震数字处理方法，只有其在实际资料处理过程中具有明显的效果才能说明该方法的实用性，才可以得到更进一步应用和推广。基于此，在完成各种理论模型的试算后，我们利用该方法对国内某油田的变偏移距VSP资料进行了试处理如图8。处理过程中只选用了Z分量，实际资料见图8(a)。该剖面在采集过程中，测线长度为3780m；炮间距为60m；采样率为1ms；采样长度为4s。我们选用了2s的记录进行了处理。在原始记录中，能够看到的有效波全部是下行直达波及其多次波，上行波被极强能量的下行波场所淹没几乎无法观察到。要在如此强的下行波背景中提取出上行波场，毫无疑问有一定的难度，尤其对于用常规的波场分离方法来说，将是更大的难度。For any seismic digital processing method, only when it has obvious effect in the actual data processing process can it prove the practicability of the method, and it can be further applied and promoted. Based on this, after completing the trial calculation of various theoretical models, we used this method to conduct a trial processing of the VSP data of a domestic oilfield, as shown in Figure 8. Only the Z component was selected during the processing, and the actual data is shown in Figure 8(a). During the acquisition process of this profile, the survey line length is 3780m; the shot spacing is 60m; the sampling rate is 1ms; the sampling length is 4s. We selected 2s records for processing. In the original records, the effective waves that can be seen are all downgoing direct waves and their multiples, and the upgoing waves are submerged by the extremely powerful downgoing wave field, which is almost impossible to observe. It is undoubtedly difficult to extract the upgoing wave field from such a strong downgoing wave background, especially for conventional wave field separation methods.

通过运用抛物拉冬变换对该变偏移距VSP地震资料进行处理后得到图8(b)所示的结果。在图8中我们发现，规则的、成组出现的波列行波均匀地展布在剖面上，这些波列都具有双曲线的特性，且曲率比较小(视速度比较大)。这就是我们所希望得到的变偏移距VSP上行波资料。原始资料中的下行波除在图边缘地带能看到下行直达波的影子外，其它的几乎看不到了，这表明我们很成功地在强背景中提取出了极弱能量的有效信号。再次体现了抛物拉冬变换方法在变偏移距VSP资料的处理中起到的重要作用。The results shown in Fig. 8(b) are obtained after processing the VSP seismic data with parabolic Radon transformation. In Fig. 8, we found that the traveling waves of regular and grouped wave trains are evenly distributed on the section. These wave trains all have the characteristics of hyperbola, and the curvature is relatively small (the apparent velocity is relatively large). This is the VSP upgoing wave data we hope to obtain. The downgoing wave in the original data can hardly be seen except the shadow of the downgoing direct wave at the edge of the figure, which shows that we have successfully extracted the effective signal with extremely weak energy in the strong background. Again, the parabolic Radon transformation method plays an important role in the processing of variable offset distance VSP data.

尽管上述通过举例说明，已经描述了本发明最佳的具体实施方式，本发明的保护范围并不仅限于上述说明，本领域一般技术人员可以理解的是，在不背离本发明所教导的实质和精髓的前提下，任何修改和变化都落入本发明的保护范围中。Although the above has described the best specific implementation of the present invention by way of example, the protection scope of the present invention is not limited to the above description, and those skilled in the art can understand that without departing from the essence and essence taught by the present invention Under the premise, any modifications and changes fall within the protection scope of the present invention.

Claims

1. A variable offset VSP wavefield separation method is characterized in that comprising: the original seismic wavefield data of the variable offset VSP collected is carried out ^t2 Stretch processing, forms ^t2 Stretch processing data; The above ^t2 stretch processing data is carried out to discrete forward parabolic Radon transformation, forming the uplink and downlink wave division data in the τ-q domain; performing discrete inverse parabolic Radon transformation on the uplink and downlink wave division data in the τ-q domain, forming xt ^2- domain data; the inverse t ² stretching process is performed on the xt ^2- domain data to complete the separation of up-going waves and down-going waves of the variable-offset VSP.

2. variable offset distance VSP wavefield separation method according to claim 1, is characterized in that: ^t2 stretching process is to read original seismic wavefield data, input relevant parameters, carry out by the formula of iflag=1 Judgment, when iflag=1 is established, perform cubic spline interpolation and stretching on the original seismic wave field data of the input parameters to form xt ² data files; when iflag=1 is not established, perform cubic spline inverse interpolation on the original seismic wave field data of the input parameters Stretch to form an xt data file, that is, form the seismic data after ^t2 stretching.

3. variable offset distance VSP wave field separation method according to claim 1, is characterized in that: the process of discrete positive parabolic Radon transformation is to read the seismic data after ^t2 stretching, input parameters, calculate the scanning of q Range rq, judge the calculated scanning range rq, if rq does not meet the requirements, return to the input parameter step, if rq meets the requirements, perform FFT processing, and perform one-dimensional filtering on the processed data; judge the filtered data, when When n<m, use M=L ^H (LL ^H +λ ² 1) ^-1 D to process the data; when n≥m, use M=(L ^H L+λ ² I) ^-1 L ^HD to process the data processing to form processed data; two-dimensional convolution is performed on the processed data, and one-dimensional filtering is performed on the two-dimensional convolution data and IFFT processing is performed to form a τ-q domain data file, that is, to form τ Separation of uplink and downlink waves in the -q domain.

4. variable offset distance VSP wave field separation method according to claim 1, is characterized in that: the process of discrete inverse parabolic Radon transformation is to read τ-q domain data file, read the parameter file of input, to data Perform FFT processing, perform frequency domain filtering for the first time on the processed data, then perform D=LM processing on the data after the first frequency domain filtering, perform second frequency domain filtering on the processed data, and perform second frequency domain filtering on the second The filtered data in the secondary frequency domain is processed by IFFT to form xt ² domain data files.

5. variable offset distance VSP wave field separation method according to claim 1, it is characterized in that: comprise reverse t ² Stretching process is to read the data described through discrete anti-parabolic Radon transformation, input correlation Parameters, judged by the formula of iflag=1, when iflag=1 is established, perform cubic spline interpolation and stretching on the data of the input parameters that have undergone discrete inverse parabolic Radon transformation, and form xt ² data files; when iflag=1 If it is not established, carry out cubic spline inverse interpolation stretching on the data that has undergone inverse parabolic Radon transformation of the input parameters, and form an xt data file, that is, form seismic data after inverse ^t2 stretching.

6. According to the described variable offset distance VSP wave field separation method according to one of claims 1, 3 or 4, it is characterized in that: it also includes verifying the parabolic Radon transformation:

First define the parabolic Radon forward transformation as:

m m ((q q,, τ τ)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} d d ((x x,, τ τ + + q q {x x}^{22})) dx dx - - - - - - ((2.2.1 2.2.1))

The inverse transform is:

d d ((x x,, t t)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} m m ((q q,, τ τ + + q q {x x}^{22})) dq dq - - - - - - ((2.2.2 2.2.2))

Since t and τ are linear in the parabolic Radon transformation, it can be processed in the frequency domain. Let D(x, ω) and M(q, ω) be d(x, t) and m(q, τ) respectively The Fourier transform of , then there are transform pairs:

M m ((q q,, ω ω)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} d d ((x x,, t t)) exp exp ((- - iωt iωt)) exp exp ((iωq iωq {x x}^{22})) dxdt wxya

That is, the forward transformation in the frequency domain is:

M m ((q q,, ω ω)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} D D. ((x x,, ω ω)) exp exp ((iωq iωq {x x}^{22})) dx dx - - - - - - ((2.2.3 2.2.3))

and

D D. ((x x,, ω ω)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} m m ((q q,, τ τ)) exp exp ((- - iωτ iωτ)) exp exp ((- - iωq iωq {x x}^{22})) dqdτ dqdτ

That is, the frequency domain inverse transform is:

D D. ((x x,, ω ω)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} M m ((q q,, ω ω)) exp exp ((- - iωq iωq {x x}^{22})) dq dq - - - - - - ((2.2.4 2.2.4))

Substitute (2.2.4) into (2.2.3) to get

{M m}^{' '} ((q q,, ω ω)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} M m (({q q}^{' '},, ω ω)) exp exp ((iω iω ((q q - - {q q}^{' '})) {x x}^{22} d d {q q}^{' '} dx dx

= M(q, ω)*σ(q, ω) (2.2.5)

The convolution factor is:

σ σ ((q q,, ω ω)) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} exp exp ((iωq iωq {x x}^{22})) dx dx

= = [[11 + + isign issign ((q q))]] \sqrt{\frac{π π}{22 ω ω | | q q | |}} - - - - - - ((2.2.6 2.2.6))

When ωq=0, σ→∞, in order to obtain a more accurate parabolic Radon forward transformation, inverse filtering is now performed on the q direction, and (2.2.5) after Fourier transform on the q direction is

M′(k _q ,ω)=M(k _q ,ω)σ(k _q ,ω) (2.2.7a)

or

M m (({k k}_{q q},, ω ω)) = = \frac{{M m}^{' '} (({k k}_{q q},, ω ω))}{σ σ (({k k}_{q q},, ω ω))}

= = \frac{\sqrt{ω ω} {M m}^{' '} (({k k}_{q q},, ω ω))}{{σ σ}^{' '} (({k k}_{q q}))} - - - - - - ((2.2.7 2.2.7 b b))

in

{σ σ}^{' '} (({k k}_{q q})) = = {&Integral; &Integral;}_{- - \infty \infty}^{+ + \infty \infty} [[11 + + isign issign ((q q))]] \sqrt{\frac{π π}{22 | | q q | |}} exp exp ((- - i i {k k}_{q q} q q)) dq dq

= = [[11 + + sign sign (({k k}_{q q}))]] \frac{π π}{\sqrt{| | {k k}_{q q} | |}}

It is easy to verify that the variable k _q of the Fourier transform corresponding to the q variable is positive, so

{σ σ}^{' '} (({k k}_{q q})) = = \frac{22 π π}{\sqrt{{k k}_{q q}}}

then have

M m (({k k}_{q q},, ω ω)) = = \frac{11}{22 π π} \sqrt{ω ω {k k}_{q q}} {M m}^{' '} (({k k}_{q q},, ω ω)) - - - - - - ((2.2.8 2.2.8))

The above formula shows that deconvolution is done in the q direction, thus improving the resolution of the parabolic Radon transform domain.

7. According to the described variable offset VSP wave field separation method according to one of claims 1, 3 or 4, it is characterized in that: it also includes discrete parabolic Radon transformation for verification, and forms a discrete parabolic Radon transformation formula:

If d(x, t) is the offset-time profile, m(q, τ) is the velocity-time profile we are looking for, and the variables x, v, t, τ represent offset, velocity, time and intercept respectively distance time, the transformation from the offset-time profile to the velocity-time profile is:

m m ((v v,, τ τ)) = = {Σ Σ}_{i i = = 11}^{n no} d d (({x x}_{i i},, t t = = \sqrt{{τ τ}^{22} + + {x x}_{i i}^{22} / / {v v}^{22}})) - - - - - - ((2.2.9 2.2.9))

n is the number of offsets, the transformation from velocity-time profile to offset-time profile is:

d d ((v v,, t t)) = = {Σ Σ}_{j j = = 11}^{m m} m m (({v v}_{j j},, τ τ = = \sqrt{{t t}^{22} + + {x x}^{22} / / {v v}_{j j}^{22}})) - - - - - - ((2.2.10 2.2.10))

m is the number of velocity sampling points;

Assuming q=1/v ² , t′=t ² , τ′=τ ² then (2.2.9) and (2.2.10) become discrete parabolic Radon transformation pairs

m m (({q q}_{j j},, τ τ)) = = {Σ Σ}_{i i = = 11}^{n no} d d (({x x}_{i i},, t t = = τ τ + + {q q}_{j j} {x x}_{i i}^{22})),, j j = = 11,, m m - - - - - - ((2.2.11 2.2.11))

d d (({x x}_{i i},, t t)) = = {Σ Σ}_{j j = = 11}^{m m} m m (({q q}_{j j},, τ τ = = t t - - {q q}_{j j} {x x}_{i i}^{22})),, i i = = 11,, n no - - - - - - ((2.2.12 2.2.12))

After Fourier transform, it is:

M m (({q q}_{j j},, ω ω)) = = {Σ Σ}_{i i = = 11}^{n no} D D. (({x x}_{i i},, ω ω)) exp exp ((iω iω {q q}_{j j} {x x}_{i i}^{22})),, j j = = 11,, m m - - - - - - ((2.2.13 2.2.13))

D D. (({x x}_{i i},, ω ω)) = = {Σ Σ}_{i i = = 11}^{m m} M m (({q q}_{j j},, ω ω)) exp exp ((- - iω iω {q q}_{j j} {x x}_{i i}^{22})),, i i = = 11,, n no - - - - - - ((2.2.14 2.2.14))

For each frequency component (2.2.13) and (2.2.14) can be expressed as the following matrix form:

M(q _j )＝L ^H D( _xi ) (2.2.15)

D(x _i )＝LM(q _j ) (2.2.16)

in

L = \exp (- iω q_{j} x_{i}^{2});

i=1,n; j=1,m

Written in matrix form as:

L L = = [\begin{matrix} exp exp ((- - iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{11}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{11}^{22})) \\ exp exp ((- - iω iω {q q}_{11} {x x}_{22}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot \cdot \cdot \cdot \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{22}^{22})) \\ \cdot &Center Dot; & \cdot \cdot & \cdot &Center Dot; \\ \cdot \cdot & \cdot &Center Dot; & \cdot \cdot \cdot &Center Dot; \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot \cdot & \cdot &Center Dot; & \cdot &Center Dot; \\ exp exp ((- - iω iω {q q}_{11} {x x}_{n no}^{22})) & exp exp ((- - iω iω {q q}_{22} {x x}_{n no}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & exp exp ((- - iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}] - - - - - - ((2.2.17 2.2.17))

L^{h} = \exp (iω q_{j} x_{i}^{2});

i=1,n; j=1,m

The matrix form is:

{L L}^{H h} = = [\begin{matrix} exp exp ((iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{11} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; & exp exp ((iω iω {q q}_{11} {x x}_{n no}^{22})) \\ exp exp ((iω iω {q q}_{22} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & exp exp ((iω iω {q q}_{22} {x x}_{n no}^{22})) \\ \cdot \cdot & \cdot &Center Dot; & \cdot \cdot \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ exp exp ((iω iω {q q}_{m m} {x x}_{11}^{22})) & exp exp ((iω iω {q q}_{m m} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; & exp exp ((iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}] - - - - - - ((2.2.18 2.2.18))

From (2.2.15) and (2.2.16), we know that L and L ^H are rectangular matrices, so L and L ^H are not real reciprocal operators, and we need to use generalized inverse to solve;

When n≥m, that is, when the equation (2.2.12) is an overdetermined equation, the m×m matrix L ^H L is reversible, so the generalized inverse matrix of L is:

L ^-1 ＝(L ^H L) ^-1 L ^H (2.2.19)

Then the discrete parabolic Radon forward transformation becomes:

M＝(L ^H L) ^-1 L ^H D (2.2.20)

When n<m, that is, when the equation (2.2.12) is an underdetermined equation, the n×n matrix LL ^H is reversible, so the generalized inverse matrix of L is:

L ^-1 = L ^H (LL ^H ) ^-1 (2.2.21)

Then the discrete parabolic Radon forward transformation becomes:

M＝L ^H (LL ^H ) ^-1 D (2.2.22)

The operator L ^H L is:

{L L}^{H h} L L = = [\begin{matrix} exp exp ((iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{11} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot \cdot \cdot \cdot & exp exp ((iω iω {q q}_{11} {x x}_{n no}^{22})) \\ exp exp ((iω iω {q q}_{22} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; & exp exp ((iω iω {q q}_{22} {x x}_{n no}^{22})) \\ \cdot \cdot & \cdot \cdot & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot \cdot \cdot \cdot \cdot \cdot & \cdot &Center Dot; \\ \cdot \cdot & \cdot \cdot & \cdot &Center Dot; \\ exp exp ((iω iω {q q}_{m m} {x x}_{11}^{22})) & exp exp ((iω iω {q q}_{m m} {x x}_{22}^{22})) & \cdot \cdot \cdot &Center Dot; \cdot \cdot & exp exp ((iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}]

\times \times [\begin{matrix} exp exp ((- - iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{11}^{22})) & \cdot &Center Dot; \cdot \cdot \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{11}^{22})) \\ exp exp ((- - iω iω {q q}_{11} {x x}_{22}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot \cdot \cdot &Center Dot; \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{22}^{22})) \\ \cdot \cdot & \cdot \cdot & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot \cdot & \cdot \cdot \cdot \cdot \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot \cdot & \cdot &Center Dot; & \cdot &Center Dot; \\ exp exp ((- - iω iω {q q}_{11} {x x}_{n no}^{22})) & exp exp ((- - iω iω {q q}_{22} {x x}_{n no}^{22})) & \cdot \cdot \cdot \cdot \cdot &Center Dot; & exp exp ((- - iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}]

= = [\begin{matrix} n no & {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{11} - - {q q}_{22})) {x x}_{i i}^{22}]] & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{11} - - {q q}_{m m})) {x x}_{i i}^{22}]] \\ {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{22} - - {q q}_{11})) {x x}_{i i}^{22}]] & n no & \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; & {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{22} - - {q q}_{m m})) {x x}_{i i}^{22}]] \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot \cdot & \cdot &Center Dot; & \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot \cdot & \cdot &Center Dot; & \cdot &Center Dot; \\ {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{m m} - - {q q}_{11})) {x x}_{i i}^{22}]] & {Σ Σ}_{i i = = 11}^{n no} exp exp [[iω iω (({q q}_{m m} - - {q q}_{22})) {x x}_{i i}^{22}]] & \cdot \cdot \cdot \cdot \cdot \cdot & n no \end{matrix}] - - - - - - ((2.2.23 2.2.23))

in the same way

L L {L L}^{H h} = = [\begin{matrix} exp exp ((- - iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{11}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & exp exp ((- - iω iω {q q}_{m m} {x x}_{11}^{22})) \\ exp exp ((- - iω iω {q q}_{11} {x x}_{22}^{22})) & exp exp ((- - ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot \cdot \cdot &Center Dot; \cdot &Center Dot; & exp exp ((- - iω iω {q q}_{m m} {x x}_{22}^{22})) \\ \cdot \cdot & \cdot &Center Dot; & \cdot \cdot \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot \cdot \cdot &Center Dot; \cdot \cdot & \cdot \cdot \\ \cdot &Center Dot; & \cdot \cdot & \cdot &Center Dot; \\ exp exp ((- - iω iω {q q}_{11} {x x}_{n no}^{22})) & exp exp ((- - iω iω {q q}_{22} {x x}_{n no}^{22})) & \cdot \cdot \cdot &Center Dot; \cdot &Center Dot; & exp exp ((- - iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}]

\times \times [\begin{matrix} exp exp ((iω iω {q q}_{11} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{11} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & exp exp ((iω iω {q q}_{11} {x x}_{n no}^{22})) \\ exp exp ((iω iω {q q}_{22} {x x}_{11}^{22})) & exp exp ((ω ω {q q}_{22} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot \cdot \cdot \cdot & exp exp ((iω iω {q q}_{22} {x x}_{n no}^{22})) \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & \cdot &Center Dot; \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ exp exp ((iω iω {q q}_{m m} {x x}_{11}^{22})) & exp exp ((iω iω {q q}_{m m} {x x}_{22}^{22})) & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & exp exp ((iω iω {q q}_{m m} {x x}_{n no}^{22})) \end{matrix}]

= = [\begin{matrix} m m & {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{22}^{22} - - {x x}_{11}^{22})) {q q}_{j j}]] & \cdot &Center Dot; \cdot \cdot \cdot &Center Dot; & {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{n no}^{22} - - {x x}_{11}^{22})) {q q}_{j j}]] \\ {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{11}^{22} - - {x x}_{22}^{22})) {q q}_{j j}]] & m m & \cdot &Center Dot; \cdot \cdot \cdot \cdot & {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{n no}^{22} - - {x x}_{22}^{22})) {q q}_{j j}]] \\ \cdot \cdot & \cdot \cdot & \cdot \cdot \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \cdot &Center Dot; \cdot \cdot & \cdot \cdot \\ \cdot &Center Dot; & \cdot &Center Dot; & \cdot &Center Dot; \\ {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{11}^{22} - - {x x}_{n no}^{22})) {q q}_{j j}]] & {Σ Σ}_{j j = = 11}^{m m} exp exp [[iω iω (({x x}_{22}^{22} - - {x x}_{n no}^{22})) {q q}_{j j}]] & \cdot &Center Dot; \cdot &Center Dot; \cdot &Center Dot; & m m \end{matrix}] - - - - - - ((2.2.24 2.2.24))

Since L ^H L and LL ^H may be singular or close to singular, especially for low-frequency components, in order to ensure the stability of the numerical solution, a damping factor is often added, so that the form of the solution becomes:

M＝(L ^H L+λ ² I) ^-1 L ^H D (2.2.25)

M＝L ^H (LL ^H +λ ² I) ^-1 D (2.2.26)

λ is the damping factor, and its value is generally 0.1~1. Therefore, the least squares problem of (2.2.20) and (2.2.22) becomes the minimum damping of (2.2.25) and (2.2.26). The quadratic problem.

8. according to one of claim 1,3 or 4 described variable offset distance VSP wavefield separation method, it is characterized in that: also comprise selecting preferred parameter:

It can be known from the signal analysis theory that for any function x(t) with a finite bandwidth B _f and its sampling rate Δt, the following relationship should be satisfied

AtB _f ≤ 2π (2.2.27)

For the discrete parabolic Radon transform, the bandwidth is

Δq Δq \leq \leq \frac{22 π π}{ω ω (({x x}_{max max}^{22} - - {x x}_{min min}^{22}))} - - - - - - ((2.2.28 2.2.28))

where x _max and x _min are the maximum and minimum offsets respectively, from which it can be seen that the critical sampling rate Δq _c of the parameter q in the discrete parabolic Radon transform is

Δ Δ {q q}_{c c} = = \frac{22 π π}{ω ω (({x x}_{max max}^{22} - - {x x}_{min min}^{22}))} - - - - - - ((2.2.29 2.2.29))

Assuming that the highest frequency of the signal in the original data is ω _max , and to make all effective frequency components valid, replace ω in the above formula with ω _max , then we have

Δ Δ {q q}_{c c} = = \frac{22 π π}{{ω ω}_{max max} (({x x}_{max max}^{22} - - {x x}_{min min}^{22}))} - - - - - - ((2.2.30 2.2.30))

Due to the aliasing frequency in the process of discrete parabolic Radon transformation, the quality of the transformation is reduced. Here, the aliasing frequency is reduced by selecting an appropriate scanning range of q. For the discrete parabolic Radon transformation, the curvature corresponding to the real event axis is set The parameter is q ₀ , then the anti-aliasing condition is:

| | q q - - {q q}_{00} | | \leq \leq \frac{π π}{ω ω | | {x x}_{max max} | | Δx Δx} - - - - - - ((2.2.31 2.2.31))

Considering all effective wave components, the above formula can be written as

Δx Δx \leq \leq \frac{π π}{{ω ω}_{max max} | | {x x}_{max max} | | (({q q}_{max max} - - {q q}_{min min}))} - - - - - - ((2.2.32 2.2.32))

From this, it can be obtained that the scanning range of the curvature parameter q satisfies the following conditions:

(({q q}_{max max} - - {q q}_{min min})) \leq \leq \frac{π π}{{ω ω}_{max max} | | {x x}_{max max} | | Δx Δx} - - - - - - ((2.2.33 2.2.33))

Therefore, in the realization of discrete parabolic Radon transform, the selection of parameters does not depend on the frequency of the signal, but only on the spatial sampling rate of the original data, the maximum offset and the highest frequency of the signal.