CN101285894B

CN101285894B - A Direct Prestack Time Migration Method for Seismic Data Acquired Under Undulating Surfaces

Info

Publication number: CN101285894B
Application number: CN2008101141336A
Authority: CN
Inventors: 张剑锋; 董春晖
Original assignee: Institute of Geology and Geophysics of CAS
Current assignee: Institute of Geology and Geophysics of CAS
Priority date: 2008-05-30
Filing date: 2008-05-30
Publication date: 2011-02-09
Anticipated expiration: 2028-05-30
Also published as: CN101285894A

Abstract

The direct pre-stack time migration method of seismic data collected under the undulating surface is applied to reflection seismic data processing in seismic exploration, and it is a pre-stack migration imaging method for seismic data collected under the undulating surface. This method can be directly applied to the seismic data collected under the undulating surface, where the shot point and the receiver point are at different elevations, which can avoid the error caused by the static correction processing of the site, and can also correct and reduce the velocity zone velocity during the migration process. The method can adaptively determine the time-varying migration aperture, can correctly compensate the geometric diffusion effect of seismic waves in the imaging process, and obtain amplitude-preserving common reflection point gathers. The core of this method is to analyze and obtain the travel time and amplitude of seismic waves determined by the velocity and thickness of the low-velocity zone, the elevation of the shot point and the receiver point, and the stacking velocity based on the floating datum by using the one-way wave theory and the principle of the stable phase point. This method has important application value for the exploration of oil, gas and mineral resources in areas with complex surface.

Description

A Direct Prestack Time Migration Method for Seismic Data Acquired Under Undulating Surfaces

技术领域technical field

本发明属于地震勘探中反射地震资料处理技术领域，涉及地震资料处理过程中的叠前偏移成像技术范畴，是一种针对起伏地表下采集的地震资料的叠前时间偏移方法。The invention belongs to the technical field of reflection seismic data processing in seismic exploration, relates to the technical category of pre-stack migration imaging in the process of seismic data processing, and is a pre-stack time migration method for seismic data collected under undulating surfaces.

背景技术Background technique

地震资料处理流程中，叠前偏移成像是最关键的环节，而叠前时间偏移是叠前偏移成像中的一种重要方法。叠前时间偏移方法可对一类倾角、断层较为复杂但速度横向变化不是很剧烈的地质构造较好成像。与叠前深度偏移方法相比，除具有较高的计算效率外，其主要的优点是只需使用叠加速度；这样可简单地通过速度扫描等方式得到恰当的速度模型，回避了使用叠前深度偏移方法面临的一个主要困难：速度见模。因此，叠前时间偏移方法已成为地震勘探领域广泛应用的关键技术。In the process of seismic data processing, prestack migration imaging is the most critical link, and prestack time migration is an important method in prestack migration imaging. The pre-stack time migration method can better image a kind of geological structures with complex dip angles and faults but less severe lateral velocity changes. Compared with the pre-stack depth migration method, in addition to its higher computational efficiency, its main advantage is that it only needs to use stacked velocities; in this way, an appropriate velocity model can be obtained simply by means of velocity scanning, avoiding the use of pre-stack One of the major difficulties faced by the depth migration method is the velocity mode. Therefore, the prestack time migration method has become a key technology widely used in the field of seismic exploration.

在地表条件比较复杂地区进行地震勘探时，地震资料的采集是在起伏的地表上进行的，各个炮点和检波点不在同一个高程上；而现行的叠前时间偏移方法是用双平方根方程求取地震波的旅行时，即应用金字塔面(pyramid)旅行时进行叠加，这要求炮点、检波点是在同一水平面上。对这一问题的解决有两种处理方法，一种是应用静校正方法，即应用静校正技术将所有炮点和检波点静校正到同一的基准面上，然后应用叠前时间偏移技术；第二种是也是应用静校正方法，但将同属于一个CMP道集的炮点和检波点静校正到该CMP点对应的基准面上，采用“浮动”而不是全区同一的基准面，这就是浮动基准面叠前时间偏移方法。前者简单但静校正量过大，因此误差也较大；目前更多地使用的是后者。尽管采用浮动基准面可通过减小静校正量及相互误差抵消来减少静校正误差，但这一技术仍存在如下几个问题：一是当有高速层(老地层)出露时，垂直入射和出射的静校正应用条件已不成立，继续使用静校正方法将带来较大的误差；二是低速层剥离和速度替换的误差难以完全抵消；三是静校正量完全依赖于已知的近地表速度模型，不能在成像过程中修正低降速带速度。When seismic exploration is carried out in areas with complex surface conditions, the acquisition of seismic data is carried out on the undulating surface, and each shot point and receiver point are not on the same elevation; while the current pre-stack time migration method uses the double square root equation To obtain the travel time of the seismic wave, that is, to superimpose the travel time of the pyramid surface, this requires that the shot point and the receiver point are on the same horizontal plane. There are two methods to solve this problem. One is to apply the static correction method, that is, to apply static correction technology to statically correct all shot points and receiver points to the same datum, and then to apply pre-stack time migration technology; The second is to apply the static correction method, but statically correct the shot points and receiver points belonging to the same CMP gather to the datum plane corresponding to the CMP point, and use "floating" instead of the same datum plane in the whole area. It is the floating datum prestack time migration method. The former is simple but the amount of static correction is too large, so the error is also relatively large; the latter is more commonly used at present. Although the use of floating datum can reduce the static correction error by reducing the amount of static correction and mutual error cancellation, there are still several problems in this technology: First, when a high-speed layer (old formation) is exposed, the vertical incidence and The application conditions of the static correction for exiting are no longer established, and continuing to use the static correction method will bring large errors; second, the errors of low-velocity layer stripping and velocity replacement are difficult to completely offset; third, the amount of static correction is completely dependent on the known near-surface velocity model, it is not possible to correct for low deceleration zone velocities during imaging.

偏移孔径的确定对叠前时间偏移是重要的，较小的孔径可减少偏移计算量，但存在着不能对陡倾角构造正确成像的风险，过大的孔径又带来了偏移噪音和较大的计算量；更好地确定偏移孔径仍然是一个需要研究的问题。The determination of the migration aperture is important for prestack time migration. A smaller aperture can reduce the amount of migration calculations, but there is a risk that the image cannot be constructed correctly for steep dip angles. An excessively large aperture brings migration noise and large computational load; better determination of the offset aperture is still a research problem.

现行的叠前时间偏移是以双平方根方程为基础，其成像条件实际上使用了叠前深度偏移的相关成像条件，因此成像没有补偿地震波传播的几何扩散效应，即不保幅。为获得保幅的共反射点(CRP)道集，以服务于叠前反演等利用叠前信息进行油气和流体检测的技术，应发展保幅的成像方法。The current prestack time migration is based on the double square root equation, and its imaging conditions actually use the relevant imaging conditions of prestack depth migration. Therefore, the imaging does not compensate the geometric diffusion effect of seismic wave propagation, that is, it does not preserve amplitude. In order to obtain amplitude-preserved common reflection point (CRP) gathers to serve pre-stack inversion and other technologies that use pre-stack information for oil, gas and fluid detection, amplitude-preserved imaging methods should be developed.

发明内容Contents of the invention

本发明的目的是：提供一种起伏地表下采集的地震资料的叠前时间偏移方法，能处理起伏地表下采集的地震资料，处理结果真实客观地反应出地下构造的走向、断点、断距大小以及对地震波的反射强度。它即不需要使用场地静校正方法将炮点和检波点校正到同一高程或浮动基准面上，又可在偏移过程中修正近地表低降速带的速度；这一方法避免了结合静校正技术的叠前时间偏移方法由于静校正产生的误差，也避免了静校正处理中获取近地表速度场的困难。The purpose of the present invention is to provide a pre-stack time migration method for seismic data collected under the undulating surface, which can process the seismic data collected under the undulating surface, and the processing results can truly and objectively reflect the trend, breakpoint, and discontinuity of the underground structure. The size of the distance and the reflection intensity of seismic waves. It does not need to use the site static correction method to correct the shot point and receiver point to the same elevation or floating datum, and can correct the velocity of the low-velocity zone near the surface during the migration process; this method avoids the combination of static correction The pre-stack time migration method of the technology also avoids the difficulty of obtaining the near-surface velocity field in the static correction process due to the error caused by the static correction.

本发明采用的技术方案是：起伏地表下采集的地震资料的直接叠前时间偏移方法，具体步骤包括：The technical scheme adopted in the present invention is: a direct pre-stack time migration method for seismic data collected under the undulating surface, and the specific steps include:

(1)用单边或双边采集方式采集起伏地表下的地震资料并读取每个检波点记录的地震信号，对炮点和检波点在不同高程的叠前地震资料进行偏移处理。(1) Use unilateral or bilateral acquisition to collect seismic data under the undulating surface and read the seismic signals recorded at each receiver point, and perform migration processing on the pre-stack seismic data of shot points and receiver points at different elevations.

(2)由已知近地表的速度等信息，确定初始的低降速带速度和低降速带时间厚度，将低降速带底部定义为浮动基准面。(2) Based on the known near-surface velocity and other information, determine the initial velocity and time thickness of the low-velocity zone, and define the bottom of the low-velocity zone as the floating datum.

(3)采用基于单程波方程和稳相点原理的查表方法，求得不同高程上炮点和检波点到成像点的地震波走时和振幅，具体是：依据低降速带速度v₁、叠加速度v_rms和炮点或检波点到成像点的水平距离x及成像点到浮动基准面的旅行时T₂和炮点或检波点到浮动基准面的旅行时T₁的变化范围，构建地震波走时和幅值查表算法中利用的表，即按等间距变化的T₁v₁，T₂v_rms，x和v₁/v_rms，求解a＝v_rmsp_x和幅值A_s并将它们填到以T₁v₁，T₂v_rms，x和v₁/v_rms为四个参数的四维表中，其p_x为解得的射线参数；由各成像点的低降速带厚度、速度、距炮点和检波点的水平距离及基于浮动基准面的叠加速度和炮点及检波点的高程，可查表求得a和A_s，进而得到炮点和检波点到成像点的地震波走时和振幅。(3) Using the table look-up method based on the one-way wave equation and the principle of stable phase points, the travel time and amplitude of seismic waves from the shot point and the receiver point to the imaging point at different elevations are obtained, specifically: according to the velocity v ₁ of the low-velocity zone, superposition Velocity v _rms and the horizontal distance x from the shot point or receiver point to the imaging point, the travel time T ₂ from the imaging point to the floating datum and the travel time T ₁ from the shot point or the receiver point to the floating datum, to construct the seismic wave travel time and the table used in the amplitude look-up table algorithm, that is, T ₁ v ₁ , T ₂ v _rms , x and v ₁ /v _rms that change at equal intervals, solve a=v _rms p _x and the amplitude A _s and convert them Fill in the four-dimensional table with T ₁ v ₁ , T ₂ v _rms , x and v ₁ /v _rms as four parameters, where p _x is the obtained ray parameter; Velocity, horizontal distance from shot point and receiver point, superimposed velocity based on floating datum and elevation of shot point and receiver point, a and A _s can be obtained by looking up the table, and then the seismic wave from the shot point and receiver point to the imaging point can be obtained travel time and amplitude.

(4)依据拟成像构造的最大倾角确定时变的偏移孔径，具体是：依据构造在不同深度(旅行时)上的最大倾角等信息，构建时变偏移孔径对应的起始旅行时表，对每一成像道，可查表求得起始旅行时。(4) Determine the time-varying migration aperture according to the maximum dip angle of the proposed imaging structure, specifically: according to the maximum dip angle of the structure at different depths (travel time), construct the initial travel time table corresponding to the time-varying migration aperture , for each imaging channel, the initial travel time can be obtained by looking up the table.

(5)对地震道循环，由每一个地震道生成保幅的成像体，具体是：对每一成像道，由起始旅行时开始，查表求得炮点和检波点到每个成像点的走时和振幅，应用保幅成像条件得到该点的成像幅值，然后将该幅值累加到该成像点对应的共反射点(CRP)道集中相应的偏移距的偏移结果上。(5) For seismic trace circulation, an imaging volume with preserved amplitude is generated from each seismic trace, specifically: for each imaging trace, starting from the initial travel time, look up the table to obtain the shot point and receiver point to each imaging point The travel time and amplitude of the imaging point are obtained by applying the amplitude-preserving imaging condition to obtain the imaging amplitude of the point, and then the amplitude is added to the migration result of the corresponding offset in the common reflection point (CRP) gather corresponding to the imaging point.

(6)在偏移过程中确定基于浮动基准面的叠加速度和修正初始低降速带的速度。(6) Determine the superposition velocity based on the floating datum and correct the velocity of the initial low-speed zone during the migration process.

(7)用偏移过程中得到的新的叠加速度场和低降速带速度场重做偏移，将所有CRP道集用动校正方法拉平得到保幅的CRP道集，将CRP道集叠加得到基于浮动基准面的偏移叠加剖面；计算浮动基准面与选定的水平基准面的垂直旅行时T_f＝T_c+(h₀-h_c)/v₁，其中h_c是地表高程、T_c是低降速带的时间厚度、v₁是修正后的低降速带速度，h₀是选定的水平基准面的高程；用T_f修正偏移叠加剖面，即按T_f移动叠加剖面的各个地震道，得到最终偏移结果。(7) Redo the migration with the new superimposed velocity field obtained during the migration process and the velocity field of the low-speed zone, flatten all the CRP gathers with the dynamic correction method to obtain the amplitude-preserved CRP gathers, and superimpose the CRP gathers Obtain the offset stacked profile based on the floating datum; calculate the vertical travel time between the floating datum and the selected horizontal datum T _f =T _c +(h ₀ -h _c )/v ₁ , where h _c is the surface elevation, T _c is the time thickness of the low-velocity zone, v ₁ is the corrected velocity of the low-velocity zone, and h ₀ is the elevation of the selected horizontal datum; use T _f to correct the offset stacking section, that is, move the stack according to T _f Each seismic trace of the profile is used to obtain the final migration result.

(8)通过常规显示软件将偏移结果转换为地下反射构造的剖面图像，剖面图像将反应地下构造的走向、断点、断距大小以及对地震波的反射强度。(8) Convert the migration result into a section image of the underground reflection structure through conventional display software. The section image will reflect the direction, breakpoint, and distance of the underground structure, as well as the reflection intensity of the seismic wave.

所述的采用基于单程波方程和稳相点原理的查表方法求得不同高程上炮点和检波点到成像点的地震波走时和振幅是这样实现的：令v₁为成像点处浮动基准面以上的低降速带的均方根速度，v_rms为成像点处基于浮动基准面的叠加速度，T₂为成像点到浮动基准面的垂直旅行时，x是炮点或检波点到成像点的水平距离；若成像点处的地表高程为h_c、该处低降速带的时间厚度为T_c、炮点或检波点的地表高程为h_s，有T₁＝T_c+(h_s-h_c)/v₁；以T₁v₁，T₂v_rms，x和v₁/v_rms为四个参数，建立四维表，表中存两个数值，一个值是由方程：The described table look-up method based on the one-way wave equation and the principle of stable phase points is used to obtain the travel time and amplitude of seismic waves from the shot point and the receiver point to the imaging point at different elevations _. The root mean square velocity of the above low-speed zone, v _rms is the stacking velocity based on the floating datum at the imaging point, T ₂ is the vertical travel time from the imaging point to the floating datum, x is the shot point or receiver point to the imaging point horizontal distance; if the surface elevation at the imaging point is h _c , the time thickness of the low velocity zone is T _c , and the surface elevation of the shot point or receiver point is h _s , then T ₁ =T _c +(h _s -h _c )/v ₁ ; with T ₁ v ₁ , T ₂ v _rms , x and v ₁ /v _rms as four parameters, a four-dimensional table is established, two values are stored in the table, and one value is determined by the equation:

$\frac{{T T}_{11} {v v}_{11}^{22} {p p}_{x x}}{\sqrt{11 - - {v v}_{11}^{22} {p p}_{x x}^{22}}} + + \frac{{T T}_{22} {v v}_{rms rms}^{22} {p p}_{x x}}{\sqrt{11 - - {v v}_{rms rms}^{22} {p p}_{x x}^{22}}} + + x x = = 00$

解得的a＝v_rmsp_x，另一个值是由a和b＝v₁/v_rms决定的幅值A_s：The obtained a=v _rms p _x , another value is the amplitude A _s determined by a and b=v ₁ /v _rms :

${A A}_{s the s} = = \sqrt{\frac{22 π π}{{v v}_{rms rms} c c}},,$ $c c = = \frac{{T T}_{11} {b b}^{22}}{((11 - - {a a}^{22} {b b}^{22})) \sqrt{11 - - {a a}^{22} {b b}^{22}}} + + \frac{{T T}_{22}}{((11 - - {a a}^{22})) \sqrt{11 - - {a a}^{22}}}$

对每一成像点，由T₁、v₁、T₂、v_rms和x计算四个参数T₁v₁，T₂v_rms，x和b＝v₁/v_rms，并按表的间距取整，可在表中查得数值a＝v_rmsp_x和幅值A_s，将a带入下式可得地震波走时t_s：For each imaging point, calculate four parameters T ₁ v ₁ , _T ₂ v _rms , x and b=v ₁ /v _rms from T ₁ , v 1 , T ₂ , v rms and x, and take them according to the _interval in the table Integral, the value a=v _rms p _x and the amplitude A _s can be found in the table, and a can be brought into the following formula to obtain the seismic wave travel time t _s :

${t t}_{s the s} = = {T T}_{11} \sqrt{11 - - {((ab ab))}^{22}} + + {T T}_{22} \sqrt{11 - - {a a}^{22}} + + ax ax / / {v v}_{rms rms} . .$

所述的依据拟成像构造的最大倾角确定时变的偏移孔径是这样实现的：由每个地震道生成一个椭圆成像体，根据拟成像构造在不同旅行时上最大倾角决定各成像道的非零值的起始旅行时，用成像道的起始旅行时定义时变的偏移孔径；起始旅行时T_b由下式确定：The described method of determining the time-varying migration aperture based on the maximum dip angle of the quasi-imaging structure is realized in the following way: an ellipse imaging volume is generated from each seismic trace; The initial travel time of zero value defines the time-varying offset aperture with the initial travel time of the imaging track; the initial travel time T _b is determined by the following formula:

${T T}_{b b} = = \frac{22 h h}{{v v}_{rms rms} ((tan the tan [[22 α α (({T T}_{b b})) - - θ θ]] - - tan the tan θ θ))},,$ $tan the tan θ θ = = \frac{x x - - h h}{{T T}_{b b} {v v}_{rms rms}}$

其中α为不同旅行时上的最大倾角，v_rms是叠加速度，h是半偏移距，x是成像点距中心点的距离。实际应用中采用查表算法，表中存放着对应不同偏移距、叠加速度和成像点距中心点距离的成像道的起始旅行时，由各个成像道的对应参数在表中拾取相应的起始旅行时。Where α is the maximum inclination at different travel times, v _rms is the stacking velocity, h is the half-offset distance, and x is the distance between the imaging point and the center point. In the practical application, the look-up table algorithm is used. The table stores the starting travel times of the imaging tracks corresponding to different offsets, stacking speeds, and distances from the imaging point to the center point. The corresponding parameters of each imaging track are picked up from the table. when starting to travel.

所述的生成保幅的成像体是这样实现的：令已得到炮点到成像点的走时为t_s、幅值为A_s，检波点到成像点的走时为t_r、幅值为A_r，令f(t)为该炮检点的地震记录，则成像点处，保幅的成像幅值为：The generation of the amplitude-preserving imaging volume is realized in the following way: let the travel time from the acquired shot point to the imaging point be t _s , the amplitude is A _s , the travel time from the receiver point to the imaging point is t _r , and the amplitude is A _r , let f(t) be the seismic record of the offset point, then at the imaging point, the imaging amplitude of amplitude preservation is:

$I I (({T T}_{22})) = = \frac{{A A}_{r r}}{{A A}_{s the s}} {f f}_{h h} (({t t}_{s the s} + + {t t}_{r r}))$

式中f_h(t)在二维情况时为f(t)的半导数，即

和

分别代表正反傅立叶变换；三维情况时f_h(t)为f(t)的一阶导数。where f _h (t) is the semi-derivative of f(t) in the two-dimensional case, namely

and

Represent the positive and negative Fourier transform respectively; f _h (t) is the first derivative of f(t) in the three-dimensional case.

所述的在偏移过程中确定基于浮动基准面的叠加速度和修正近地表低降速带的速度是这样实现的：首先对地震资料作初步的速度拾取得到初始的、横向均匀的叠加速度场，基于这一速度场和初始的低降速带速度进行偏移；对偏移生成的共反射点道集，依据横向均匀的叠加速度场做反动校，再做动校正确定基于浮动基准面的叠加速度；选定一高信噪比的同相轴，基于新的叠加速度场和依百分比增加和减少初始低降速带速度，对该同相轴的邻域做局部偏移，选取使同相轴最平的速度，对这一速度做空间平滑处理，即得到修正的低降速带速度场。In the migration process, the determination of the stacking velocity based on the floating datum and the correction of the velocity of the near-surface low-velocity zone are realized in the following way: firstly, preliminary velocity picking is performed on the seismic data to obtain an initial, transversely uniform stacking velocity field , based on this velocity field and the initial velocity of the low-speed zone, the migration is performed; for the common reflection point gather generated by the migration, the reverse correction is performed according to the horizontally uniform superimposed velocity field, and then the dynamic correction is performed to determine the Stacking velocity: select an event with a high SNR, based on the new stacking velocity field and increase and decrease the initial low-speed zone velocity by percentage, make a local offset to the neighborhood of the event, and select the event so that the If the speed is flat, the spatial smoothing process is performed on this speed, that is, the corrected low-speed zone velocity field is obtained.

本发明的起伏地表下采集的地震资料的直接叠前时间偏移方法，能直接在偏移过程中正确考虑起伏地表形状和近地表低降速带对偏移成像的影响。The direct pre-stack time migration method of the seismic data collected under the undulating surface of the present invention can directly and correctly consider the influence of the shape of the undulating surface and the low velocity zone near the surface on migration imaging during the migration process.

本发明的起伏地表下采集的地震资料的直接叠前时间偏移方法，能在成像过程中修正低降速带速度，方法对低降速带时间厚度的选择是稳键的，不要求其底面准确对应实际构造的高速层顶面。The direct pre-stack time migration method of the seismic data collected under the undulating surface of the present invention can correct the velocity of the low-velocity zone during the imaging process, and the method is stable for the selection of the time thickness of the low-velocity zone, and does not require its bottom surface Accurately correspond to the top surface of the high-speed layer of the actual structure.

本发明的起伏地表下采集的地震资料的直接叠前时间偏移方法，能在成像过程中正确补偿地震波的几何扩散效应，得到保幅的共反射点道集。The direct pre-stack time migration method of the seismic data collected under the undulating surface of the present invention can correctly compensate the geometric diffusion effect of the seismic wave in the imaging process, and obtain the amplitude-preserving common reflection point gather.

本发明的起伏地表下采集的地震资料的直接叠前时间偏移方法，能依据拟成像构造在不同深度上的最大倾角，自适应地确定时变的偏移孔径。The direct pre-stack time migration method of the seismic data collected under the undulating surface of the present invention can adaptively determine the time-varying migration aperture according to the maximum dip angle of the pseudo-imaging structure at different depths.

本发明的起伏地表下采集的地震资料的直接叠前时间偏移方法，可灵活应用于各类地震道集，适宜于并行实现。The direct pre-stack time migration method of seismic data collected under the undulating surface of the present invention can be flexibly applied to various seismic gathers and is suitable for parallel implementation.

本发明的具体实现原理如下：Concrete realization principle of the present invention is as follows:

本发明的核心思想是从叠前深度偏移方法出发分析和研究叠前时间偏移，将叠前时间偏移看作是叠前深度偏移方法在简化速度模型下的一个近似(用旅行时替代深度)。The core idea of the present invention is to analyze and study the pre-stack time migration from the pre-stack depth migration method, and regard the pre-stack time migration as an approximation of the pre-stack depth migration method under the simplified velocity model (using travel time alternative depth).

首先基于单程波理论和稳相点原理，得到计算跑点或检波点至成像点的地震波走时和幅值的解析表达式，这一公式是以成像点处基于浮动基准面的叠加速度和浮动基准面顶部的降速带均方根速度为基础的(以下讨论将以二维问题为例)。Firstly, based on the theory of one-way wave and the principle of stable phase point, the analytical expressions for calculating the travel time and amplitude of seismic waves from the running point or receiver point to the imaging point are obtained. This formula is based on the stacking velocity and floating datum at the imaging point Based on the root mean square velocity of the deceleration band at the top of the surface (the following discussion will take the two-dimensional problem as an example).

在波数-频率域，单个炮点或检波点的地震波场的深度延拓可表示为：In the wavenumber-frequency domain, the depth continuation of the seismic wavefield of a single shot point or receiver point can be expressed as:

$P P (({k k}_{x x},, ω ω,, {Σ Σ}_{i i = = 11}^{n no} {t t}_{i i})) = = F f ((ω ω)) exp exp ((- - j j {Σ Σ}_{i i = = 11}^{n no} {t t}_{i i} \sqrt{{ω ω}^{22} - - {α α}_{i i}^{22} {k k}_{x x}^{22}})) - - - - - - ((11))$

其中F(ω)是炮点的震源或检波点的地震记录的付氏变换，t_i和α_i分别是各层介质的时间厚度和层间速度。式(1)中的相移量可近似为Where F(ω) is the Fourier transform of the source or receiver seismic records of the shot point, and t _i and α _i are the time thickness and interlayer velocity of each layer of media, respectively. The phase shift in formula (1) can be approximated as

${Σ Σ}_{i i = = 11}^{n no} {t t}_{i i} \sqrt{{ω ω}^{22} - - {α α}_{i i}^{22} {k k}_{x x}^{22}} \approx \approx {T T}_{11} \sqrt{{ω ω}^{22} - - {v v}_{11}^{22} {k k}_{x x}^{22}} + + {T T}_{22} \sqrt{{ω ω}^{22} - - {v v}_{rms rms}^{22} {k k}_{x x}^{22}} - - - - - - ((22))$

若设第m层界面是浮动基准面，则 $T_{1} = Σ_{i = 1}^{m} t_{i},$ $v_{1}^{2} = \frac{Σ_{i = 1}^{m} α_{i}^{2} t_{i}}{T_{1}},$ $T_{2} = Σ_{i = m + 1}^{n} t_{i},$ $v_{rms}^{2} = \frac{Σ_{i = m + 1}^{n} α_{i}^{2} t_{i}}{T_{2}} .$ 式(2)表明，浮动基准面并不必须对应实际构造的高速层顶面；但若这一基准面若准确分隔强速度对比界面(即分隔低降速带和高速层顶面)，将有效提高式(2)的近似精度。将(2)式的近似相移量代入(1)式的深度延拓公式并作空间付氏反变换，得空间-频率域波场为If the interface of layer m is assumed to be a floating reference plane, then $T_{1} = Σ_{i = 1}^{m} t_{i},$ $v_{1}^{2} = \frac{Σ_{i = 1}^{m} α_{i}^{2} t_{i}}{T},$ $T_{2} = Σ_{i = m + 1}^{no} t_{i},$ $v_{rms}^{} = \frac{Σ_{i = m + 1}^{no} α_{i}^{2} t_{i}}{T} .$ Equation (2) shows that the floating datum level does not necessarily correspond to the top surface of the actual high-velocity layer; Improve the approximate accuracy of formula (2). Substituting the approximate phase shift in Eq. (2) into the depth extension formula in Eq. (1) and performing inverse Fourier transform in space, the wave field in the space-frequency domain is obtained as

$P P ((x x,, ω ω,, {T T}_{11} + + {T T}_{22})) = = \frac{11}{22 π π} &Integral; &Integral; F f ((ω ω)) exp exp [[- - jω jω (({T T}_{11} \sqrt{11 - - {v v}_{11}^{22} {(({k k}_{x x} / / ω ω))}^{22}} + + {T T}_{22} \sqrt{11 - - {v v}_{rms rms}^{22} {(({k k}_{x x} / / ω ω))}^{22}} - - \frac{x x {k k}_{x x}}{ω ω}))]] d d {k k}_{x x} - - - - - - ((33))$

令p_x＝k_x/ω，式(3)可转换为关于p_x的积分，即Let p _x ＝k _x /ω, formula (3) can be transformed into the integral about p _x , that is

$P P ((x x,, ω ω,, {T T}_{11} + + {T T}_{22})) = = \frac{ω ω}{22 π π} &Integral; &Integral; F f ((ω ω)) exp exp [[- - jω jω (({T T}_{11} \sqrt{11 - - {v v}_{11}^{22} {p p}_{x x}^{22}} + + {T T}_{22} \sqrt{11 - - {v v}_{rms rms}^{22} {p p}_{x x}^{22}} - - {p p}_{x x} x x))]] d d {p p}_{x x} - - - - - - ((44))$

式中p_x是射线参数，与频率无关，因此对(4)式可应用稳相点原理求得渐近解为In the formula, p _x is the ray parameter, which has nothing to do with frequency, so the asymptotic solution can be obtained by applying the principle of stable phase point to formula (4) as

$P P ((x x,, ω ω,, {T T}_{11} + + {T T}_{22})) = = F f ((ω ω)) \sqrt{ω ω} exp exp ((- - j j \frac{π π}{44})) \sqrt{\frac{11}{22 π π | | {φ φ}^{' '' '} (({p p}_{x x}^{00})) | |}} exp exp ((- - jωφ jωφ (({p p}_{x x}^{00})))) - - - - - - ((55))$

式中 $φ (p_{x}) = T_{1} \sqrt{1 - v_{1}^{2} p_{x}^{2}} + T_{2} \sqrt{1 - v_{rms}^{2} p_{x}^{2}} - p_{x} x, φ^{''} (p_{x})$ 是其二阶导数，p_x ⁰是φ(p_x)的一阶导数的零点，φ(p_x ⁰)和 $\sqrt{\frac{1}{2 π | φ^{''} (p_{x}^{0}) |}}$ 即是地震波的走时和幅值；而p_x ⁰可由如下方程求得：In the formula $φ (p_{x}) = T_{1} \sqrt{1 - v_{1}^{2} p_{x}^{2}} + T_{2} \sqrt{1 - v_{rms}^{2} p_{x}^{2}} - p_{x} x, φ^{''} (p_{x})$ is its second derivative, p _x ⁰ is the zero point of the first derivative of φ(p _x ), φ(p _x ⁰ ) and $\sqrt{\frac{1}{2 π | φ^{''} (p_{x}^{0}) |}}$ That is, the travel time and amplitude of the seismic wave; and p _x ⁰ can be obtained by the following equation:

$\frac{{T T}_{11} {v v}_{11}^{22} {p p}_{x x}}{\sqrt{11 - - {v v}_{11}^{22} {p p}_{x x}^{22}}} + + \frac{{T T}_{22} {v v}_{rms rms}^{22} {p p}_{x x}}{\sqrt{11 - - {v v}_{rms rms}^{22} {p p}_{x x}^{22}}} + + x x = = 00 - - - - - - ((66))$

式(5)和(6)即是炮点或检波点至成像点的地震波走时和幅值的解析计算式，式中v₁为该浮动基准面以上的低降速带的均方根速度，v_rms为基于浮动基准面的叠加速度，T₂是成像点距浮动基准面的旅行时，T₁是由成像点处的地表高程h_c、低降速带的时间厚度T_c、炮点(或检波点)的地表高程h_s(或h_r)确定的，有T₁＝T_c+(h_s-h_c)/v₁。Equations (5) and (6) are the analytical calculation formulas of seismic wave traveltime and amplitude from the shot point or receiver point to the imaging point, where _v1 is the root mean square velocity of the low velocity zone above the floating datum, v _rms is the stacking velocity based on the floating datum, T ₂ is the travel time between the imaging point and the floating datum, T ₁ is the surface elevation h _c at the imaging point, the time thickness T _c of the low velocity zone, and the shot point ( Or receiver point) surface elevation h _s (or h _r ) determined, there is T ₁ =T _c +(h _s -h _c )/v ₁ .

本发明基于直射线假设并忽略地表起伏来确定偏移孔径，以下讨论将以二维问题为例。本发明采用输入道成像方式，每个地震道将生成以炮点和检波点为一对焦点的椭圆成像体，时变偏移孔径是用椭圆成像体中各成像道的非零值的起始旅行时来决定的。若要对倾角为α的构造成像，这一成像椭圆在这一深度的最大张开角度就应大于等于α(半椭圆将对应着90度)；而从炮点发出的地震射线将入射到这一倾角为α的构造并恰好反射到检波点，因此依据下面的几何关系确定偏移孔径对应的起始旅行时T_h：The present invention determines the offset aperture based on the straight ray assumption and ignores the surface fluctuation. The following discussion will take the two-dimensional problem as an example. The present invention adopts the input channel imaging method, and each seismic channel will generate an elliptical imaging volume with the shot point and the receiver point as a pair of focal points, and the time-varying offset aperture is the initial value of each imaging channel in the elliptical imaging volume. It is decided when traveling. If the structure with dip angle α is to be imaged, the maximum opening angle of this imaging ellipse at this depth should be greater than or equal to α (the semi-ellipse will correspond to 90 degrees); and the seismic rays emitted from the shot point will be incident on this A structure with an inclination angle of α just reflects to the receiver point, so the initial travel time T _h corresponding to the migration aperture is determined according to the following geometric relationship:

T_bv_rmstan[2α(T_b)-θ]-T_bv_rmstanθ＝2h (7)T _b v _rms tan[2α(T _b )-θ]-T _b v _rms tanθ＝2h (7)

T_bv_rmstanθ＝x-hT _b v _rms tanθ=xh

式中α为不同深度(旅行时)上成像构造的最大倾角，v_rms是成像点的叠加速度，θ是炮点出射的地震波的入射角，而2α-θ是界面反射波的出射角，2h是偏移距，x是成像点距中心点的距离。由(7)式求得对应于不同x的T_b，成像时每一成像道由T_b开始成像，即实现时变的偏移孔径。where α is the maximum dip angle of the imaging structure at different depths (travel time), v _rms is the stacking velocity of the imaging point, θ is the incident angle of the seismic wave emitted from the shot point, and 2α-θ is the exit angle of the interface reflection wave, 2h is the offset, and x is the distance from the imaging point to the center point. T _b corresponding to different x is obtained from formula (7), and each imaging track starts to be imaged from T _b during imaging, that is, a time-varying offset aperture is realized.

对每个地震道，既然已基于相移法和稳相点原理求得了炮点和检波点到各成像点的幅值和走时，利用波动方程深度偏移的反褶积成像条件获得真幅值成像。设炮点和检波点至成像点的走时和幅值分别为t_s、A_s和t_r、A_r，若设震源是一时间脉冲，有真幅值成像结果(二维情况)：For each seismic trace, since the amplitude and travel time from the shot point and receiver point to each imaging point have been obtained based on the phase shift method and the principle of stable phase point, the true amplitude value can be obtained by using the deconvolution imaging condition of wave equation depth migration imaging. Let the travel time and amplitude from the shot point and the receiver point to the imaging point be t _s , A _s and t _r , A _r respectively, if the seismic source is a time pulse, there is a true amplitude imaging result (two-dimensional case):

$I I ((x x,, T T)) = = \frac{{A A}_{r r}}{{A A}_{s the s}} &Integral; &Integral; F f ((ω ω)) \sqrt{ω ω} exp exp ((- - j j \frac{π π}{44})) exp exp ((- - jω jω (({t t}_{s the s} + + {t t}_{r r})))) dω dω - - - - - - ((88))$

积分中的前三项可看作是该道地震记录的半导数的付氏变换，记其在时间域为f_h(t)，则(8)式可简化为 $I (x, T) = \frac{A_{r}}{A_{s}} f_{h} (t_{s} + t_{r}),$ 即在做求半导数后的地震记录中拾取t_s+t_r时刻的值，然后乘上由幅值确定的权系数。由于实际地震记录是离散的，上述拾取可通过四点插值实现。The first three items in the integral can be regarded as the Fourier transform of the semi-derivative of the seismic record, which is recorded as f _h (t) in the time domain, then the formula (8) can be simplified as $I (x, T) = \frac{A_{r}}{A_{the s}} f_{h} (t_{the s} + t_{r}),$ That is, pick the value at time t _s + t _r from the seismic record after calculating the semi-derivative, and then multiply it by the weight coefficient determined by the amplitude. Since the actual seismic records are discrete, the above picking can be realized by four-point interpolation.

本发明的有益效果：该方法可直接应用于起伏地表下采集的、炮点和检波点在不同高程的地震资料，生成叠前时间偏移道集和图像，这不仅省却了场地静校正的处理环节，也避免了静校正技术带来的误差。该方法能在偏移过程中修正降低速带速度，且允许降低速带底面不准确对应实际构造的高速层顶面，这避免了使用静校正等处理方法中准确获取近地表速度模型的困难。该方法能生成保幅的共反射点(CRP)道集，更好地服务于叠前反演等油气和流体检测技术。显示的剖面图像反应出地下构造的走向、断点、断距大小以及对地震波的反射强度。该方法对地表复杂地区的油气、矿产资源勘探有重要应用价值。Beneficial effects of the present invention: the method can be directly applied to the seismic data collected under the undulating surface, where the shot point and the receiver point are at different elevations, to generate pre-stack time migration gathers and images, which not only saves the processing of site static correction link, and also avoid the error caused by the static correction technology. This method can correct the velocity of the deceleration zone during the migration process, and allows the bottom surface of the deceleration zone to inaccurately correspond to the top surface of the actual high-velocity layer, which avoids the difficulty of accurately obtaining the near-surface velocity model in processing methods such as static correction. This method can generate amplitude-preserved common reflection point (CRP) gathers, which can better serve oil and gas and fluid detection technologies such as pre-stack inversion. The displayed section image reflects the direction, breakpoint, and distance of the underground structure, as well as the reflection intensity of seismic waves. This method has important application value for the exploration of oil, gas and mineral resources in areas with complex surface.

附图说明Description of drawings

图1是采集地震资料并进行叠前时间偏移的具有起伏地表的断陷盆地质构造图，构造表层没有明显的低速带，在近地表存在一层最大厚度20m左右的风化层，速度为1400m/s；其它各部分的速度由图中数字所示(单位为m/s)。我们将基于这一地质构造，由声波方程全离散数值方法进行正演模拟，生成模拟炮记录，以此验证和解释起伏地表下采集的地震资料的直接叠前时间偏移方法。Figure 1 is a geological structure diagram of a faulted basin with undulating surface obtained by seismic data and pre-stack time migration. There is no obvious low-velocity zone on the surface of the structure. There is a weathered layer with a maximum thickness of about 20m near the surface and a velocity of 1400m. /s; the speed of other parts is shown by the figures in the figure (the unit is m/s). Based on this geological structure, we will perform forward modeling by the full discrete numerical method of the acoustic wave equation to generate simulated shot records, so as to verify and interpret the direct prestack time migration method of seismic data collected under the undulating surface.

图2是由正演生成的151个炮记录中抽取的零偏移距剖面，浅层面波、直达波及近地表低速带的反射波已被剔除。图中第一个反射同相轴对应图1中水平界面的反射。白线是选定的浮动基准面。Figure 2 is a zero-offset profile extracted from 151 shot records generated by forward modeling. Shallow layer waves, direct arrival waves, and reflected waves from the near-surface low-velocity zone have been eliminated. The first reflection event in the figure corresponds to the reflection of the horizontal interface in Figure 1. The white line is the selected floating datum.

图3给出了在初始低降速带速度和横向均匀的初始叠加速度场下CDP＝7500m点的CRP道集；从图中可看出，同相轴存在剩余动校，也存在表面和层间多次波对应的同相轴。Figure 3 shows the CRP gather at point CDP=7500m under the initial low-velocity zone velocity and the transversely uniform initial superimposed velocity field; it can be seen from the figure that there are residual dynamic corrections in the event, as well as surface and interlayer The event corresponding to the multiple.

图4给出了在新的叠加速度场下对图3的CDP点作部分偏移的速度扫描结果；图中横轴是选取的速度与初始速度的比值；从图中可确定，该CDP点的最佳低降速带速度应为初始低降速带速度的95％。Fig. 4 shows the velocity scanning results of the partial offset of the CDP point in Fig. 3 under the new superimposed velocity field; the horizontal axis in the figure is the ratio of the selected velocity to the initial velocity; it can be determined from the figure that the CDP point The optimal low-speed belt speed should be 95% of the initial low-speed belt speed.

图5是修正后的低降速带速度与初始的均匀低降速带速度的比较。Figure 5 is a comparison of the corrected low deceleration belt speed to the original uniform low deceleration belt speed.

图6给出了与图3相同的CDP＝7500m点，在更新后的叠加速度和降速带速度下的CRP道集，由图中可见同相轴更加平直，聚焦也更好(较小的剩余动校量可在叠加流程中修正)。Fig. 6 shows the same CDP=7500m point as Fig. 3, and the CRP gather under the updated stacking velocity and deceleration zone velocity. It can be seen from the figure that the event axis is straighter and the focus is better (smaller The remaining dynamic calibration can be corrected in the superposition process).

图7是以浮动基准面为零时刻的偏移叠加剖面。Fig. 7 is the offset superposition section at the time when the floating datum is zero.

图8是以264m处水平面为基准面的修正的偏移叠加剖面。在图8中图1水平界面对应的偏移同相轴也变得平直，这表明修正后的低降速带速度是正确的。图1与图8对比可知，三条规模较大的断层，两个小构造均被很好成像，断面也被很好地成像。Figure 8 is the corrected offset superimposed section with the horizontal plane at 264m as the reference plane. In Fig. 8, the offset event corresponding to the horizontal interface in Fig. 1 also becomes flat, which indicates that the corrected velocity of the low deceleration zone is correct. Comparing Figure 1 with Figure 8, it can be seen that the three large-scale faults and the two small structures are well imaged, and the cross-sections are also well imaged.

图9是贵州水城地区实际采集的地震资料的偏移结果的局部比较。上图是采用常规静校正方法的偏移结果，下图是本发明的直接叠前时间偏移方法的偏移结果，图9中直线段是断层解释的结果，由断层解释可知，本发明方法很好地刻画了地下构造的走向和断点及断距。Figure 9 is a partial comparison of the migration results of the seismic data actually collected in the Shuicheng area of Guizhou. The upper figure is the migration result of the conventional static correction method, and the lower figure is the migration result of the direct pre-stack time migration method of the present invention. The straight line segment in Figure 9 is the result of fault interpretation. From the fault interpretation, it can be seen that the method of the present invention The trend, breakpoint and distance of the underground structure are well described.

具体实施方式Detailed ways

起伏地表下采集的地震资料的直接叠前时间偏移方法，针对图1的具有起伏地表的断陷盆地地质构造，用数值方法合成了151炮的中间放炮两边采集的地震记录，以这组炮点和检波点在不同高程的叠前地震资料为例，具体为以下步骤：The direct pre-stack time migration method of seismic data collected under the undulating surface, aiming at the geological structure of the faulted basin with undulating surface in Fig. Taking the pre-stack seismic data of the points and receiver points at different elevations as an example, the specific steps are as follows:

(1)针对双边采集方式得到的叠前地震资料，读取每个炮记录的每个检波点的地震记录，对炮点和检波点在不同高程的叠前地震资料直接进行偏移处理；对全部地震记录求半导数；采用并行处理，按偏移距范围将地震记录分配到各个计算机节点，图2是零偏移距的地震记录。(1) For the pre-stack seismic data obtained by bilateral acquisition, read the seismic records of each receiver point recorded by each shot, and directly perform migration processing on the pre-stack seismic data of shot points and receiver points at different elevations; Calculate the semi-derivative of all seismic records; use parallel processing to distribute seismic records to each computer node according to the offset range. Figure 2 shows the zero offset seismic records.

(2)确定初始的低降速带速度和低降速带时间厚度，将低降速带底部定义为浮动基准面。图2中的白线是我们定义的低降速带底部，也是浮动基准面，而白线在图2上的时间就是低降速带时间厚度。用横向均匀的速度(2313m/s)作为初始的低降速带速度。(2) Determine the initial velocity of the low deceleration zone and the time thickness of the low deceleration zone, and define the bottom of the low deceleration zone as the floating datum. The white line in Figure 2 is the bottom of the low-speed zone we defined, which is also the floating reference level, and the time of the white line in Figure 2 is the time thickness of the low-speed zone. Use the horizontal uniform speed (2313m/s) as the initial low speed belt speed.

(3)依据低降速带速度v₁、叠加速度v_rms和炮点或检波点到成像点的水平距离x及成像点到浮动基准面的旅行时T₂和炮点或检波点到浮动基准面的旅行时T₁的变化范围，构建地震波走时和幅值查表算法中利用的表，即按等间距变化的T₁v₁，T₂v_rms，x和v₁/v_rms的值的组合，求解方程：(3) According to the velocity v ₁ of the low-velocity zone, the stacking velocity v _rms , the horizontal distance x from the shot point or receiver point to the imaging point, and the travel time T ₂ from the imaging point to the floating reference plane and the distance from the shot point or receiver point to the floating reference plane The variation range of travel time T ₁ of the surface, construct the table used in the table look-up algorithm of seismic wave travel time and amplitude, that is, the values of T ₁ v ₁ , T ₂ v _rms , x and v ₁ /v _rms that change at equal intervals Combine, solving the equations:

得到a＝v_rmsp_x，将a和b＝v₁/v_rms代入下式得幅值A_s：Obtain a=v _rms p _x , substituting a and b=v ₁ /v _rms into the following formula to obtain the amplitude A _s :

然后将a和A_s填到以T₁v₁，T₂v_rms，x和v₁/v_rms为四个参数的四维表中。Then fill a and A _s into a four-dimensional table with T ₁ v ₁ , T ₂ v _rms , x and v ₁ /v _rms as four parameters.

(4)依据构造在不同深度(旅行时)上的最大倾角等信息构建时变偏移孔径对应的起始旅行时表。令α(T)为不同旅行时上的最大倾角，v_rms是叠加速度，h是半偏移距，x是成像点距中心点的距离，按等间距变化的v_rms、h和x的值的组合，求解方程：(4) Construct the initial travel time table corresponding to the time-varying migration aperture based on information such as the maximum inclination at different depths (travel time). Let α(T) be the maximum inclination angle at different travel times, v _rms is the stacking velocity, h is the half-offset distance, x is the distance from the imaging point to the center point, and the values of v _rms , h and x vary at equal intervals A combination of , solving the equation:

将解得的T_h填到以v_rms、h和x为三个参数的三维表中。成像时由各个成像道的对应参数在表中拾取相应的起始旅行时T_b。Fill the obtained T _h into the three-dimensional table with v _rms , h and x as three parameters. During imaging, the corresponding initial travel time T _b is picked up from the table by the corresponding parameters of each imaging track.

(5)对地震道循环，具体步骤包括：1)由地震道的偏移距，成像点的位置和叠加速度查表确定各道的起始旅行时T_b，若T_b大于最大成像深度或地震记录的长度，表明成像点超出了成像范围，不需进行成像。2)对每一成像点，由成像点顶部的低降速带速度v₁和时间厚度T_c、成像点处的地表高程h_c、炮点(或检波点)的地表高程h_s(或h_r)确定T₁＝T_c+(h_s-h_c)/v₁，令T₂＝T_b+nΔT(ΔT是选定的偏移剖面的时间深度采样)，根据成像点处的降速带速度v₁和基于浮动基准面的叠加速度v_rms，以及成像点距炮点和检波点的水平距离x，由查表算法得到a＝v_rmsp_x和幅值A_s，令b＝v₁/v_rms，将a和b带入下式可得地震波走时t_s：(5) For the seismic trace cycle, the specific steps include: 1) by the offset distance of the seismic trace, the position of the imaging point and the stacking velocity look-up table to determine the initial travel time T _b of each trace, if T _b is greater than the maximum imaging depth or The length of the seismic record, indicating that the imaging point is beyond the imaging range and need not be imaged. 2) For each imaging point, the velocity v ₁ and the time thickness T _c of the low-speed zone at the top of the imaging point, the surface elevation h _c at the imaging point, and the surface elevation h _s (or h _r ) Determine T ₁ =T _c +(h _s -h _c )/v ₁ , let T ₂ =T _b +nΔT (ΔT is the time-depth sampling of the selected migration profile), according to the deceleration at the imaging point With the velocity v ₁ and the stacking velocity v _rms based on the floating datum, and the horizontal distance x between the imaging point and the shot point and the receiver point, a=v _rms p _x and the amplitude A _s are obtained by the look-up table algorithm, and b=v ₁ /v _rms , put a and b into the following formula to get the seismic wave travel time t _s :

${t t}_{s the s} = = {T T}_{11} \sqrt{11 - - {((ab ab))}^{22}} + + {T T}_{22} \sqrt{11 - - {a a}^{22}} + + ax ax / / {v v}_{rms rms}$

这样可得到炮点和检波点至成像点的走时和幅值t_s、A_s和t_r、A_r(当t_s+t_r大于地震记录的长度或T₂大于等于最大成像深度时停止n的增加)，应用保幅的成像条件得成像幅值I(T₂)：In this way, the travel time and _amplitude t _s , A _s and t _r , A _r from the shot point and the receiver point to the imaging point can _be obtained ₍ stop n increase), the imaging amplitude I(T ₂ ) is obtained by applying the amplitude-preserving imaging condition:

式中f_h(t)为该道地震记录的半导数。3)依据查表得到的v_rmsp_x，计算炮点和检波点发出的地震射线与该成像点处浮动基准面的交点，据此求得在这一基准面上的偏移距。4)根据偏移距值，将成像幅值I(T₂)累加到该成像点对应的CRP道集中相应的偏移距的结果上。where f _h (t) is the semi-derivative of the seismic record. 3) According to the v _rms p _x obtained by looking up the table, calculate the intersection point of the seismic rays emitted by the shot point and the receiver point and the floating datum at the imaging point, and obtain the offset on this datum accordingly. 4) According to the offset value, the imaging amplitude I(T ₂ ) is added to the corresponding offset result in the CRP gather corresponding to the imaging point.

(6)在偏移过程中确定基于浮动基准面的叠加速度和修正初始低降速带速度。具体步骤包括：1)沿粗网格选定CMP道集，通过NMO技术和横向平均，确定一个横向均匀的叠加速度场，基于这一速度场和初始的低降速带速度进行偏移，图3是偏移生成的一个CRP道集。2)对全部CRP道集循环，依据横向均匀的初始叠加速度场做反动校，再通过动校正确定基于浮动基准面的叠加速度场。3)沿较细网格选定成像点，由对应的CRP道集中选定一高信噪比的同相轴，基于新的叠加速度场和依百分比增加和减少初始低降速带速度，重新对同相轴的邻域做局部偏移成像，图4是依百分比增加和减少初始低降速带速度时的一个CRP道集局部的部分偏移的结果。4)选取使同相轴最平的速度为修正后的低降速带速度，如图4中能确定百分比为95％，对这一速度做空间平滑处理，得到修正的低降速带速度场。参阅图5。(6) Determine the stacking velocity based on the floating reference plane and correct the initial low-speed zone velocity during the migration process. The specific steps include: 1) Select the CMP gathers along the coarse grid, determine a horizontally uniform superimposed velocity field through NMO technology and horizontal averaging, and perform migration based on this velocity field and the initial velocity of the low velocity zone, as shown in Fig. 3 is a CRP gather generated by migration. 2) For all the CRP gather cycles, the reverse correction is made according to the horizontally uniform initial superimposed velocity field, and then the superimposed velocity field based on the floating datum is determined through dynamic correction. 3) Select imaging points along the finer grid, select an event with high SNR from the corresponding CRP gather, based on the new superimposed velocity field and increase and decrease the velocity of the initial low-velocity zone according to the percentage, re-image The neighborhood of the event is locally migrated and imaged. Figure 4 shows the results of a local partial migration of a CRP gather with percentage increases and decreases in initial low-slow zone velocity. 4) Select the velocity that makes the event the flattest as the corrected low-velocity zone velocity, as shown in Figure 4, the percentage can be determined as 95%, and spatial smoothing is performed on this velocity to obtain the corrected low-velocity zone velocity field. See Figure 5.

(7)用新的叠加速度场和低降速带速度场重做偏移，图6是生成的新的CRP道集，将所有CRP道集用动校正方法拉平即得到保幅的共反射点(CRP)道集；将CRP道集叠加，即得到基于浮动基准面的偏移成像剖面，如图7所示；由地表高程h_c、低降速带的时间厚度T_c和修正后的低降速带速度v₁，计算浮动基准面与选定的高程在h₀的水平基准面的垂直旅行时T_f＝T_c+(h₀-h_c)/v₁，对图7的各个地震道按T_f作时间移动，获得基于同一水平基准面的偏移结果。(7) Redo the migration with the new superimposed velocity field and the velocity field of the low-velocity zone. Figure 6 shows the generated new CRP gathers. Flatten all the CRP gathers with the dynamic correction method to obtain the amplitude-preserved common reflection points (CRP ₎ gathers; the CRP gathers are superimposed to obtain the migration imaging section based on the floating datum, as shown in Fig. 7 _; Velocity v ₁ of the deceleration zone, calculate the vertical travel time T _f =T _c +(h ₀ -h _c )/v ₁ of the floating datum and the selected elevation at the horizontal datum of h ₀ , for each earthquake in Fig. 7 The trace moves in time according to T _f to obtain migration results based on the same horizontal datum.

(8)通过显示软件将偏移结果转换为地下反射构造的剖面图像。参阅图8。(8) Convert the migration result into a section image of the underground reflection structure through the display software. See Figure 8.

Claims

1. The direct pre-stack time migration method of the seismic data collected under the undulating surface, is characterized in that the following steps are adopted: A) collect the seismic data under the undulating surface and read each receiver point with unilateral or bilateral acquisition mode The recorded seismic signals are migrated to the pre-stack seismic data of the shot point and the receiver point at different elevations; B) the initial low-velocity zone velocity and the time thickness of the low-velocity zone are determined from the known near-surface information, and the The bottom of the low-velocity zone is defined as the floating datum; C) Using the table look-up method based on the one-way wave equation and the principle of stable phase points to obtain the travel time and amplitude of seismic waves from the shot point and the receiver point to the imaging point at different elevations; D) According to the simulated The maximum dip angle of the imaging structure determines the time-varying migration aperture; E) For the seismic trace cycle, an amplitude-preserving imaging volume is generated from each seismic trace, and the amplitude of the imaging volume is added to the common reflection point gather according to the offset distance ; F) Determining the superimposed velocity based on the floating datum and correcting the velocity of the near-surface low-speed zone during the migration process; G) Redoing with the final superimposed velocity field and the low-velocity zone velocity field obtained during the migration process Offset, use the dynamic correction method to eliminate the remaining dynamic correction amount, obtain the intermediate offset result, calculate the vertical travel time between the floating datum and the selected horizontal datum, use this travel time to correct the intermediate offset result, and obtain the final offset Result; H) The migration result is converted into a section image of the underground reflection structure through conventional display software, and the section image reflects the direction, breakpoint, and distance of the underground structure, as well as the reflection intensity of the seismic wave.

2. the direct pre-stack time migration method of the seismic data collected under a kind of undulating surface according to claim 1, is characterized in that: described adopting the look-up table method based on one-way wave equation and stable phase point principle obtains The travel time and amplitude of seismic waves from the shot point and the receiver point to the imaging point at different elevations are realized as follows: let v ₁ be the root mean square velocity of the low-velocity zone above the floating datum at the imaging point, and v _rms be the imaging point at the imaging point based on The stacking speed of the floating datum, T ₂ is the vertical travel time from the imaging point to the floating datum, x is the horizontal distance from the shot point or receiver point to the imaging point; T ₁ is the travel time from the shot point or receiver point to the floating datum ; If the surface elevation at the imaging point is h _c , the time thickness of the low velocity zone is T _c , and the surface elevation of the shot point or receiver point is h _s , then T ₁ =T _c +(h _s -h _c )/v ₁ ; with T ₁ v ₁ , T ₂ v _rms , x and v ₁ /v _rms as four parameters, a four-dimensional table is established, two values are stored in the table, and one value is determined by the equation:

\frac{{T T}_{11} {v v}_{11}^{22} {p p}_{x x}}{\sqrt{11 - - {v v}_{11}^{22} {p p}_{x x}^{22}}} + + \frac{{T T}_{22} {v v}_{rms rms}^{22} {p p}_{x x}}{\sqrt{11 - - {v v}_{rms rms}^{22} {p p}_{x x}^{22}}} + + x x = = 00

The obtained a=v _rms p _x , p _x is the obtained ray parameter, and the other value is the amplitude A _s determined by a and b=v ₁ /v _rms :

{A A}_{s the s} = = \sqrt{\frac{22 π π}{{v v}_{rms rms} c c}},,

c c = = \frac{{T T}_{11} {b b}^{22}}{((11 - - {a a}^{22} {b b}^{22})) \sqrt{11 - - {a a}^{22} {b b}^{22}}} + + \frac{{T T}_{22}}{((11 - - {a a}^{22})) \sqrt{11 - - {a a}^{22}}}

For each imaging point, calculate four parameters T ₁ v ₁ , _T ₂ v _rms , x and b=v ₁ /v _rms from T ₁ , v 1 , T ₂ , v rms and x, and take them according to the _interval in the table Integral, the value a=v _rms p _x and the amplitude A _s can be found in the table, and a can be brought into the following formula to obtain the seismic wave travel time t _s :

{t t}_{s the s} = = {T T}_{11} \sqrt{11 - - {((ab ab))}^{22}} + + {T T}_{22} \sqrt{11 - - {a a}^{22}} + + ax ax / / {v v}_{rms rms} . .

3. the direct pre-stack time migration method of seismic data collected under a kind of undulating surface according to claim 1, is characterized in that: the described migration aperture that determines time-varying according to the maximum dip angle of the imaginary imaging structure is like this Realized: Generate an elliptical imaging volume from each seismic trace, determine the non-zero initial travel time of each imaging trace according to the maximum dip angle of the quasi-imaging structure at different travel times, and use the initial travel time of the imaging trace to define time-varying The offset aperture of ; _Tb at the initial travel time is determined by:

{T T}_{b b} = = \frac{22 h h}{{v v}_{rms rms} ((tan the tan [[22 α α (({T T}_{b b})) - - θ θ]] - - tan the tan θ θ))},,

tan the tan θ θ = = \frac{x x - - h h}{{T T}_{b b} {v v}_{rms rms}}

Among them, α is the maximum inclination angle at different travel times, v _rms is the stacking velocity, h is the half-offset distance, and x is the distance between the imaging point and the center point. The initial travel time of the imaging track, the stacking speed, and the distance from the imaging point to the center point, the corresponding starting travel time is picked up from the table by the corresponding parameters of each imaging track.

4. the direct pre-stack time migration method of seismic data collected under a kind of undulating surface according to claim 1, is characterized in that: the imaging volume of described generating amplitude preservation is realized in this way: let the obtained shot point The travel time to the imaging point is t _s , the amplitude is A _s , the travel time from the receiver point to the imaging point is t _r , and the amplitude is A _r , let f(t) be the seismic record of the shot detection point, then at the imaging point, The imaging amplitude of amplitude preservation is:

I I (({T T}_{22})) = = \frac{{A A}_{r r}}{{A A}_{s the s}} {f f}_{h h} (({t t}_{s the s} + + {t t}_{r r}))

where f _h (t) is the semi-derivative of f(t) in the two-dimensional case, namely

and

5. The method for direct pre-stack time migration of seismic data collected under a undulating surface according to claim 1, characterized in that: in the migration process, determining the stacking velocity based on the floating datum and correcting the near The velocity of the low-velocity zone on the surface is realized in the following way: firstly, preliminary velocity picking is performed on the seismic data to obtain an initial, transversely uniform stacked velocity field, and migration is performed based on this velocity field and the initial velocity of the low-velocity zone; The common-reflection point gather generated by migration is reversely corrected according to the horizontally uniform stacking velocity field, and then the stacking velocity based on the floating datum is determined by dynamic correction; an event with a high signal-to-noise ratio is selected, and based on the new stacking velocity The field sum increases and decreases the velocity of the initial low-velocity zone according to the percentage, makes a local offset to the neighborhood of the event, selects the velocity that makes the event the flattest, and performs spatial smoothing on this velocity, that is, the corrected low-velocity zone is obtained Speed belt velocity field.