Seismic data interpolation using generalised velocity‐dependent seislet transform

Data interpolation is an important step for seismic data analysis because many processing tasks, such as multiple attenuation and migration, are based on regularly sampled seismic data. Failed interpolations may introduce artifacts and eventually lead to inaccurate final processing results. In this paper, we generalised seismic data interpolation as a basis pursuit problem and proposed an iteration framework for recovering missing data. The method is based on non‐linear iteration and sparse transform. A modified Bregman iteration is used for solving the constrained minimisation problem based on compressed sensing. The new iterative strategy guarantees fast convergence by using a fixed threshold value. We also propose a generalised velocity‐dependent formulation of the seislet transform as an effective sparse transform, in which the non‐hyperbolic normal moveout equation serves as a bridge between local slope patterns and moveout parametres in the common‐midpoint domain. It can also be reduced to the traditional velocity‐dependent seislet if special heterogeneity parametre is selected. The generalised velocity‐dependent seislet transform predicts prestack reflection data in offset coordinates, which provides a high compression of reflection events. The method was applied to synthetic and field data examples, and the results show that the generalised velocity‐dependent seislet transform can reconstruct missing data with the help of the modified Bregman iteration even for non‐hyperbolic reflections under complex conditions, such as vertical transverse isotropic (VTI) media or aliasing.     

Velocity model updating in prestack Kirchhoff time migration for PS converted waves: Part I – Theory

A velocity model updating approach is developed based on moveout analysis of the diffraction curve of PS converted waves in prestack Kirchhoff time migration. The diffraction curve can be expressed as a product of two factors: one factor depending on the PS converted‐wave velocity only, and the other factor depending on all parameters. The velocity‐dependent factor represents the hyperbolic behaviour of the moveout and the other is a scale factor that represents the non‐hyperbolic behaviour of the moveout. This non‐hyperbolic behaviour of the moveout can be corrected in prestack Kirchhoff time migration to form an inverse normal‐moveout common‐image‐point gather in which only the hyperbolic moveout is retained. This hyperbolic moveout is the moveout that would be obtained in an isotropic equivalent medium. A hyperbolic velocity is then estimated from this gather by applying hyperbolic moveout analysis. Theoretical analysis shows that for any given initial velocity, the estimated hyperbolic velocity converges by an iterative procedure to the optimal velocity if the velocity ratio is optimal or to a value closer to the optimal velocity if the velocity ratio is not optimal. The velocity ratio (V_P/V_S) has little effect on the estimation of the velocity. Applying this technique to a synthetic seismic data set confirms the theoretical findings. This work provides a practical method to obtain the velocity model for prestack Kirchhoff time migration.     

Suppressing random noise to broaden the useful bandwidth of seismic reflection data

The signal-to-noise ratio （SNR） of seismic reflection data in many areas is rather poor and conventional two-dimensional filters designed to suppress noise with different moveout from the signal tend to generate artifacts. We have extended a method of multichannel filtering, based on the hypothesis that signals on adjacent channels are similar, for enhancing the SNR on stacked sections. Using only the mid-range frequencies where the SNR is highest, the event trend is found for overlapping windows on the section and the average signal vector is calculated. Then the data from the full bandwidth section are projected onto the spatially varying unit similarity vectors and the results are merged for the overlapping windows. Application of the method to synthetic data containing steeply dipping events and to a stacked section for a marine 2D line has produced good results. The modifications we have introduced carry a small overhead in computing time but they should enable the method to be used effectively even on sections containing steep dips.     

可展宽地震资料有效频宽的随机噪声衰减方法

很多地区地震资料的信噪比较低,而用于压制与信号具有不同方向的随机噪声的常规二维滤波方法常常产生假信息。基于相邻信号具有相干性这一假设,本文提出了一种叠后衰减随机噪声的多道滤波方法。该方法利用信噪比最高的中频段信息(含有主频的这一频率区间)分时窗计算信号单位矢量,并将该时窗内全频段数据向信号单位矢量方向投影,对各时窗(包括时间方向和空间方向)重叠部分按比例进行加权。我们利用这种方法对含有陡倾角的合成地震数据和海上二维实际地震资料进行了处理,处理效果很好。这种方法较为费时,但不受倾角限制,应用范围广。     

STACKING OF P-SV SEISMIC REFLECTION DATA USING DIP MOVEOUT1

The moveout of P-SV mode-converted seismic reflection events in a common-midpoint gather is non-hyperbolic. This is true even if the medium has constant P- and SV-wave velocities. Furthermore, reflection-point smear occurs even along horizontal reflectors. These effects reduce the resolution of the zero-offset stack. In such a medium, the generalization of the dip moveout transformation to P-SV data can be calculated analytically. The resulting P-SV dip moveout operators solve the problem of reflection-point smear, and image any reflector regardless of dip or depth. The viability of this technique is demonstrated on synthetic and field data.     

起伏地表煤田地震资料静校正

由于地表起伏和近地表结构变化产生的静校正问题严重影响了煤田地震资料的成像质量.为此,首先利用低速带分片拟合的广义线性反演技术进行折射波静校正,解决长波长静校正问题和部分短波长静校正问题,然后,利用叠加能量最大静校正技术进一步解决剩余静校正问题,最后,利用非地表一致性剩余时差校正技术,解决速度和射线等误差引起的非地表一致性剩余时差问题.实验结果表明,在以串连的方式应用了三种校正方法之后,在共炮点道集上,折射渡同相轴的线性形态得到了恢复;在动校正后的共中心点道集上,煤层反射的双曲线同相轴被拉平;在叠加剖面上,煤层反射的信噪比得到了改善.     

APPLICATION OF THE PARTIAL KARHUNEN-LOÈVE TRANSFORM TO SUPPRESS RANDOM NOISE IN SEISMIC SECTIONS1

The Karhunen-Loève (K–L) transform is an effective technique for suppressing spatially uncorrelated noise, but because of its high computational cost, fast transforms, such as the Fourier transform, have been more favoured. Two techniques that combine to make the K–L transform feasible for seismic data processing are discussed. The first technique filters the data for limited dips. For each dip, linear moveout is applied to the seismic sections so that events with this dip are made flat. By interpolation, we can include dips that are fractions of a sample/trace. After linear moveout, zero-lag K–L filtering is applied followed, by inverse linear moveout; the results from all dips are added to form the final filtered data. The second technique is blocking, in which the seismic section is divided into blocks small enough for each block to be processed using relatively small matrices; the processed blocks are assembled to form the final filtered section. Using a combination of these techniques, seismic sections can be filtered at a reasonable cost using the K-L transform.     

平面声波在粗糙界面上的反射特征研究

基于有关粗糙界面的Rayleigh假设,讨论了平面声波按余弦规律快速变化的小尺度粗糙界面上的反射特征.研究表明：这类界面与位于该位置的一个过渡地层的作用相当.该过渡层的厚度为粗糙界面的起伏幅度,速度和密度为上下两层介质相应量的平均值.研究了埋藏很深的微粗糙界面所引起的地震反射（绕射）波的频散特性和走时的构成,即包含零炮检距反射时间、正常时差和界面粗糙时差三部分内容.该粗糙时差与空间坐标和时间坐标无关,与绕射波的阶次有关,绕射波尾随在反射波之后以某一固定的时差出现.且只有当界面的粗糙波长与地震波的波长相当时,才能观测到这类绕射波.该结论为粗糙界面地震反射资料的处理方法提供了理论依据.     

Seismics and optics: hyperbolae and curvatures

We rederive and generalize hyperbolic moveout formulae for the common-midpoint (CMP) gather and for the zero-offset (ZO) section that can be efficiently used for macro-model-independent reflection imaging in two-dimensional media. The hyperbolic moveout formulae for the common-midpoint gather are obtained from different Taylor series expansions of a particular parametric moveout surface defined in the multicoverage data space. Such a moveout surface involves three kinematic wave-field attributes of two hypothetical waves, which have to be determined by a coherency analysis. By using hyperbolic moveout curves in the CMP gather and in the ZO section one can determine these attributes in two steps. The relationships between the local shapes of the interfaces and the attributes of the hypothetical wave-fields attributes are considered by means of geometrical optics. The determination of these attributes allows to perform a macro-model-independent ZO simulation and a subsequent inversion.     

Multifocusing revisited ‐ inhomogeneous media and curved interfaces*

We review the multifocusing method for traveltime moveout approximation of multicoverage seismic data. Multifocusing constructs the moveout based on two notional spherical waves at each source and receiver point, respectively. These two waves are mutually related by a focusing quantity. We clarify the role of this focusing quantity and emphasize that it is a function of the source and receiver location, rather than a fixed parameter for a given multicoverage gather. The focusing function can be designed to make the traveltime moveout exact in certain generic cases that have practical importance in seismic processing and interpretation. The case of a plane dipping reflector (planar multifocusing) has been the subject of all publications so far. We show that the focusing function can be generalized to other surfaces, most importantly to the spherical reflector (spherical multifocusing). At the same time, the generalization implies a simplification of the multifocusing method. The exact traveltime moveout on spherical surfaces is a very versatile and robust formula, which is valid for a wide range of offsets and locations of source and receiver, even on rugged topography. In two‐dimensional surveys, it depends on the same three parameters that are commonly used in planar multifocusing and the common‐reflection surface (CRS) stack method: the radii of curvature of the normal and normal‐incidence‐point waves and the emergence angle. In three dimensions the exact traveltime moveout on spherical surfaces depends on only one additional parameter, the inclination of the plane containing the source, receiver and reflection point. Comparison of the planar and spherical multifocusing with the CRS moveout expression for a range of reflectors with increasing curvature shows that the planar multifocusing can be remarkably accurate but the CRS becomes increasingly inaccurate. This can be attributed to the fact that the CRS formula is based on a Taylor expansion, whereas the multifocusing formulae are double‐square root formulae. As a result, planar and spherical multifocusing are better suited to model the moveout of diffracted waves.     

EFFICIENT MULTICHANNEL FILTERING OF SEISMIC DATA1

Multichannel filters are used to eliminate coherent noise from surface seismic data, for wavefield separation from VSP stacks, and for signal enhancement. Their success generally depends on the choice of the filter parameters and the domain of application. Multichannel filters can be applied to shots (monitors), common-receiver traces, CDP traces and stacked sections. Cascaded applications in these domains are currently performed in the seismic industry for better noise suppression and for signal enhancement. One-step shot-domain filtering is adequate for some applications. However, in practice, cascaded applications in shot-and common-receiver domains usually give better results when the S/N ratio is low. Multichannel filtering after stacking (especially after repeated applications in shot and/or receiver domains) may create undesirable results such as artificial continuations, or smearing and smoothing of small features such as small throw faults and fine stratigraphic details. Consequently, multichannel filtering after stacking must be undertaken with the utmost care and occasionally only as a last resort. Multichannel filters with fan-shaped responses (linear moveout filters) should be applied after NMO correction. These are the filters commonly used in the seismic industry where they have such names as velocity filters, moveout filters, f-k filters and coherency filters. Filtering before NMO correction may result in break-up and flattening especially of those shallow reflection events with relatively higher curvatures and diffractions. NMO correction is needed prior to wavefield separation from VSP stacks for the same practical reasons outlined above whenever source-receiver offsets are involved. Creation of artificial lineup and smearing at the outputs of multichannel filters is presently the common practical concern. Optimum multichannel filters with well-defined pass, reject and transition bands overcome the latter problems when applied before stacking and after NMO correction. The trace dimension of these filters must be kept small to avoid such lineups and the smoothing of small structures. Good results can be obtained with only five traces, but seven traces seems to be a better compromise both in surface and well seismic applications. The so-called f-k filtering and τ-p domain filtering are no exceptions to the above practical considerations. Residual static computations after multichannel filtering also need special consideration. Since multichannel filtering improves spatial continuity, residual static algorithms using local correlation, i.e. nonsurface-consistent algorithms, may be impractical especially after multichannel filtering.     

A FREQUENCY DOMAIN MAPPING APPROXIMATION OF MOVEOUT FILTERS WITH APPLICATIONS*

Wavenumber aliasing is the main limitation of conventional optimum least-squares linear moveout filters: it prevents adequate reject domain weighting for efficient coherent noise rejection. A general frequency domain multichannel filter design technique based on a one-to-one mapping method between two-dimensional (2D) space and one-dimensional (1D) space is presented. The 2D desired response is mapped to the 1D frequency axis after a suitable sorting of the coefficients. A min-max or Tchebycheff approximation to the desired response is obtained in the 1D frequency domain and mapped back to the 2D frequency domain. The algorithm is suitable for multiband 2D filter design. No aliasing damage is inherent in the linear moveout filters designed using this technique because the approximation is done in the frequency-wavenumber (f, k)-domain. Linear moveout filters designed by using the present coefficient mapping technique achieve better pass domain approximations than the corresponding conventional least-squares filters. Compatible reject domain approximations can be obtained from suitable mappings of the origin coefficient of the desired (f k)-response to the 1D frequency axis. The (fk)-responses of linear moveout filters designed by using the new technique show equi-ripple behavior. Synthetic and real data applications show that the present technique is superior to the optimum least-squares filters and straight stacking in recovering and enhancing the signal events with relatively high residual statics. Their outputs also show higher resolution than those of the optimum least-squares filters.     

Goos-Hänchen效应对自由表面全反射SV波校正时差的影响

利用SV地震波（偏振化方向在入射面内的横波）在自由表面反射系数的附加相角导出了SV波Goos-Hänchen效应所引起的横向偏移和横向偏移渡越时间,给出了Goos-Hänchen效应正常时差动校正量,讨论了Goos-Hänchen效应对反射SV波正常时差的影响,绘出了横向偏移、Goos-Hänchen效应正常时差及Goos-Hänchen效应正常时差校正量曲线.曲线表明：反射地震波存在负横向偏移效应,在大多数入射角度范围,横向偏移（横向偏移渡越时间）与波长（周期）为同一个数量级.对掠入射波或入射角在临界角附近的入射波,Goos-Hänchen效应对正常时差有较大的测量误差,对反射P、SV波的传播走时产生了不可忽略的影响,因此在实际的地震资料处理中应进行横向偏移效应误差校正.     

A new parameterization for generalized moveout approximation,based on three rays

A simple and accurate traveltime approximation is important in many applications in seismic data processing, inversion and modelling stages. Generalized moveout approximation is an explicit equation that approximates reflection traveltimes in general two-dimensional models. Definition of its five parameters can be done from properties of finite offset rays, for general models, or by explicit calculation from model properties, for specific models. Two versions of classical finite-offset parameterization for this approximation use traveltime and traveltime derivatives of two rays to define five parameters, which makes them asymmetrical. Using a third ray, we propose a balance between the number of rays and the order of traveltime derivatives. Our tests using different models also show the higher accuracy of the proposed method. For acoustic transversely isotropic media with a vertical symmetry axis, we calculate a new moveout approximation in the generalized moveout approximation functional form, which is explicitly defined by three independent parameters of zero-offset two-way time, normal moveout velocity and anellipticity parameter. Our test shows that the maximum error of the proposed transversely isotropic moveout approximation is about 1/6 to 1/8 of that of the moveout approximation that had been reported as the most accurate approximation in these media. The higher accuracy is the result of a novel parameterization that do not add any computational complexity. We show a simple example of its application on synthetic seismic data.     

全反射SH地震波的Goos-Hänchen效应动校正时差

利用SH地震波（偏振化方向垂直入射面的横波）在地层界面反射系数的附加相角导出了SH波Goos-Hänchen效应所引起的横向偏移和横向偏移渡越时间,给出了Goos-Hänchen效应正常时差公式,讨论了Goos-Hänchen效应对反射SH波正常时差的影响,绘出了横向偏移、横向偏移渡越时间、Goos-Hänchen效应正常时差及Goos-Hänchen效应正常时差校正量曲线.数值算例表明：对掠入射波或入射角在临界角附近的入射波,SH反射波的横向偏移、横向偏移渡越时间非常大,Goos-Hänchen效应对正常时差会产生较大的测量误差,在其他角度的入射波,横向偏移（横向偏移渡越时间）与波长（周期）为同一个数量级.横向偏移效应对SH反射波的传播走时影响是不可忽略的,因此在实际的地震资料处理中应进行横向偏移效应误差校正.     

HORIZONTAL STACKING AND MULTICHANNEL FILTERING APPLIED TO COMMON DEPTH POINT SEISMIC DATA*

The common depth point method of shooting in oil exploration provides a series of seismic traces which yield information about the substrata layers at one location. After normal moveout and static corrections have been applied, the traces are combined by horizontal stacking, or linear multichannel filtering, into a single record in which the primary reflections have been enhanced relative to the multiple reflections and random noise. The criterion used in optimum horizontal stacking is to maximize the signal to noise power ratio, where signal refers to the primary reflection sequence and noise includes the multiple reflections. It is shown when this criterion is equivalent to minimizing the mean square difference between the desired signal (primary reflection sequence) and the weighted horizontally stacked traces. If the seismic traces are combined by multichannel linear filtering, the primary reflection sequence will have undergone some phase and frequency distortion on the resulting record. The signal to noise power ratio then becomes less meaningful a criterion for designing the optimum linear multichannel filter, and the mean square criterion is adopted. In general, however, since more a priori information about the seismic traces is required to design the optimum linear multichannel filter than required for the optimum set of weights of the horizontal stacking process, the former will be an improvement over the latter. It becomes evident that optimum horizontal stacking is a restricted form of linear multichannel filtering.     

HORIZONTAL STACKING AND MULTICHANNEL FILTERING APPLIED TO COMMON DEPTH POINT SEISMIC DATA*

The common depth point method of shooting in oil exploration provides a series of seismic traces which yield information about the substrata layers at one location. After normal moveout and static corrections have been applied, the traces are combined by horizontal stacking, or linear multichannel filtering, into a single record in which the primary reflections have been enhanced relative to the multiple reflections and random noise. The criterion used in optimum horizontal stacking is to maximize the signal to noise power ratio, where signal refers to the primary reflection sequence and noise includes the multiple reflections. It is shown when this criterion is equivalent to minimizing the mean square difference between the desired signal (primary reflection sequence) and the weighted horizontally stacked traces. If the seismic traces are combined by multichannel linear filtering, the primary reflection sequence will have undergone some phase and frequency distortion on the resulting record. The signal to noise power ratio then becomes less meaningful a criterion for designing the optimum linear multichannel filter, and the mean square criterion is adopted. In general, however, since more a priori information about the seismic traces is required to design the optimum linear multichannel filter than required for the optimum set of weights of the horizontal stacking process, the former will be an improvement over the latter. It becomes evident that optimum horizontal stacking is a restricted form of linear multichannel filtering.     

Moveout approximation for horizontal transversely isotropic and vertical transversely isotropic layered medium. Part I: 1D ray propagation‡

Anisotropy in subsurface geological models is primarily caused by two factors: sedimentation in shale/sand layers and fractures. The sedimentation factor is mainly modelled by vertical transverse isotropy (VTI), whereas the fractures are modelled by a horizontal transversely isotropic medium (HTI). In this paper we study hyperbolic and non‐hyperbolic normal reflection moveout for a package of HTI/VTI layers, considering arbitrary azimuthal orientation of the symmetry axis at each HTI layer. We consider a local 1D medium, whose properties change vertically, with flat interfaces between the layers. In this case, the horizontal slowness is preserved; thus, the azimuth of the phase velocity is the same for all layers of the package. In general, however, the azimuth of the ray velocity differs from the azimuth of the phase velocity. The ray azimuth depends on the layer properties and may be different for each layer. In this case, the use of the Dix equation requires projection of the moveout velocity of each layer on the phase plane. We derive an accurate equation for hyperbolic and high‐order terms of the normal moveout, relating the traveltime to the surface offset, or alternatively, to the subsurface reflection angle. We relate the azimuth of the surface offset to its magnitude (or to the reflection angle), considering short and long offsets. We compare the derived approximations with analytical ray tracing.     

SEISMOGRAM SYNTHESIS AND RECOMPRESSION OF DISPERSIVE IN-SEAM SEISMIC MULTIMODE DATA USING A NORMAL-MODE SUPERPOSITION APPROACH

In spite of a geometrical rotation into radial and transverse parts, two- or three-component in-seam seismic data used for underground fault detection often suffer from the problem of overmoding ‘noise’. Special recompression filters are required to remove this multimode dispersion so that conventional reflection seismic data processing methods, e.g. CMP stacking techniques, can be applied afterwards. A normal-mode superposition approach is used to design such multimode recompression filters. Based on the determination of the Green's function in the far-field, the normal-mode superposition approach is usually used for the computation of synthetic single- and multi-mode (transmission) seismograms for vertically layered media. From the filter theory's point of view these Green's functions can be considered as dispersion filters which are convolved with a source wavelet to produce the synthetic seismograms. Thus, the design of multimode recompression filters can be reduced to a determination of the inverse of the Green's function. Two methods are introduced to derive these inverse filters. The first operates in the frequency domain and is based on the amplitude and phase spectrum of the Green's function. The second starts with the Green's function in the time domain and calculates two-sided recursive filters. To test the performance of the normal-mode superposition approach for in-seam seismic problems, it is first compared and applied to synthetic finite-difference seismograms of the Love-type which include a complete solution of the wave equation. It becomes obvious that in the case of one and two superposing normal modes, the synthetic Love seam-wave seismograms based on the normal-mode superposition approach agree exactly with the finite-difference data if the travel distance exceeds two dominant wavelengths. Similarly, the application of the one- and two-mode recompression filters to the finite-difference data results in an almost perfect reconstruction of the source wavelet already two dominant wavelengths away from the source. Subsequently, based on the dispersion analysis of an in-seam seismic transmission survey, the normal-mode superposition approach is used both to compute one- and multi-mode synthetic seismograms and to apply one- and multimode recompression filters to the field data. The comparison of the one- and two-mode synthetic seismograms with the in-seam seismic transmission data reveals that arrival times, duration and shape of the wavegroups and their relative excitation strengths could well be modelled by the normal-mode superposition approach. The one-mode recompressions of the transmission seismograms result in non-dispersive wavelets whose temporal resolution and signal-to-noise ratio could clearly be improved. The simultaneous two-mode recompressions of the underground transmission data show that, probably due to band-limitation, the dispersion characteristics of the single modes could not be evaluated sufficiently accurately from the field data in the high-frequency range. Additional techniques which overcome the problem of band-limitation by modelling all of the enclosed single-mode dispersion characteristics up to the Nyquist frequency will be mandatory for future multimode applications.     

APPROXIMATIONS TO SHEAR-WAVE VELOCITY AND MOVEOUT EQUATIONS IN ANISOTROPIC MEDIA1

Backus and Crampin derived analytical equations for estimating approximate phase-velocity variations in symmetry planes in weakly anisotropic media, where the coefficients of the equations are linear combinations of the elastic constants. We examine the application of similar equations to group-velocity variations in off-symmetry planes, where the coefficients of the equations are derived numerically. We estimate the accuracy of these equations over a range of anisotropic materials with transverse isotropy with both vertical and horizontal symmetry axes, and with combinations of transverse isotropy yielding orthorhombic symmetry. These modified equations are good approximations for up to 17% shear-wave anisotropy for propagations in symmetry planes for all waves in all symmetry systems examined, but are valid only for lower shear-wave anisotropy (up to 11%) in off-symmetry planes. We also obtain analytical moveout equations for the reflection of qP-, qSH-, and qSV- waves from a single interface for off-symmetry planes in anisotropic symmetry. The moveout equation consists of two terms: a hyperbolic moveout and a residual moveout, where the residual moveout is proportional to the degree of anisotropy and the spread length of the acquisition geometry. Numerical moveout curves are computed for a range of anisotropic materials to verify the analytical moveout equations.