VSP上下行反射波联合成像方法研究

VSP资料上下行波场发育丰富.本文在分析VSP直达波、上行反射波、下行反射波传播路径及其照明范围的基础上,指出了常规VSP波动方程偏移方法缺陷,进而通过修改波场延拓方式,提出了上下行反射波联合成像方法,并在高频近似下分析了该方法的成像原理.该方法不需要进行VSP上下行反射波场分离,能够同时对VSP资料中的一次反射波、自由表面多次波、层间多次波进行成像,比常规成像剖面具有更宽的成像范围和更好的成像效果.该方法能够对下行一次反射波进行成像,从而可以实现常规偏移方法难以处理的高陡倾角构造成像.模拟资料和实际资料处理证明了本文方法的正确性.     

Q-LOG DETERMINATION ON DOWNGOING WAVELETS AND TUBE WAVE ANALYSIS IN VERTICAL SEISMIC PROFILES1

Seismic attenuation introduces modifications in the wavelet shape in vertical seismic profiles. These modifications can be quantified by measuring particular signal attributes such as rise-time, period and shape index. Use of signal attributes leads to estimations of a seismic-attenuation log (Q-log). To obtain accurate signal attributes it is important to minimize noise influence and eliminate local interference between upgoing and downgoing waves at each probe location. When tube waves are present it is necessary to eliminate them before performing separation of upgoing and downgoing events. We used a trace-by-trace Wiener filter to minimize the influence of tube waves. The separation of upgoing and downgoing waves was then performed in the frequency domain using a trace-pair filter. We used three possible methods based on signal attribute measurements to obtain g-log from the extracted downgoing wavefield. The first one uses a minimum phasing filter and the arrival time of the first extremum. The two other methods determine the Q-factor from simple relations between the amplitudes of the first extrema and the pseudo-periods of the down-going wavelet. The relations determined between a signal attribute and traveltime over quality factor were then calibrated using field source signature and constant-Q models computed by Ganley's method. Q-logs thus obtained from real data are discussed and compared with geological information, specifically at reservoir level. Analysis of the tube wave arrivals at the level of the reservoir showed a tube wave attenuation that could not be explained by simple transmission effects. There was also a loss of signal coherence. This could be interpreted as tube wave diffusion in the porous reservoir, followed by dispersion. If this interpretation can be verified, tube wave analysis could lead to further characterization of porous permeable zones.     

Marine seismic wavefield measurement to remove sea-surface multiples

We propose a new method for removing sea-surface multiples from marine seismic reflection data in which, in essence, the reflection response of the earth, referred to a plane just above the sea-floor, is computed as the ratio of the plane-wave components of the upgoing wave and the downgoing wave. Using source measurements of the wavefield made during data acquisition, three problems associated with earlier work are solved: (i) the method accommodates source arrays, rather than point sources; (ii) the incident field is removed without simultaneously removing part of the scattered field; and (iii) the minimum-energy criterion to find a wavelet is eliminated. Pressure measurements are made in a horizontal plane in the water. The source can be a conventional array of airguns, but must have both in-line and cross-line symmetry, and its wavefield must be measured and be repeatable from shot to shot. The problem is formulated for multiple shots in a two-dimensional configuration for each receiver, and for multiple receivers in a two-dimensional configuration for each shot. The scattered field is obtained from the measurements by subtracting the incident field, known from measurements at the source. The scattered field response to a single incident plane wave at a single receiver is obtained by transforming the common-receiver gather to the frequency–wavenumber domain, and a single component of this response is obtained by Fourier transforming over all receiver coordinates. Each scattered field component is separated into an upgoing wave and a downgoing wave using the zero-pressure condition at the water-surface. The upgoing wave may then be expressed as a reflection coefficient multiplied by the incident downgoing wave plus a sum of scattered downgoing plane waves, each multiplied by the corresponding reflection coefficient. Keeping the upgoing scattered wave fixed, and using all possible incident plane waves for a given frequency, yields a set of linear simultaneous equations for the reflection coefficients which are solved for each plane wave and for each frequency. To create the shot records that would have been measured if the sea-surface had been absent, each reflection coefficient is multiplied by complex amplitude and phase factors, for source and receiver terms, before the five-dimensional Fourier transformation back to the space–time domain.     

Characterization of seismic hazard and structural response by energy flux

Seismic safety of structures depends on the structure's ability to absorb the seismic energy that is transmitted from ground to structure. One parameter that can be used to characterize seismic energy is the energy flux. Energy flux is defined as the amount of energy transmitted per unit time through a cross-section of a medium, and is equal to kinetic energy multiplied by the propagation velocity of seismic waves. The peak or the integral of energy flux can be used to characterize ground motions. By definition, energy flux automatically accounts for site amplification. Energy flux in a structure can be studied by formulating the problem as a wave propagation problem. For buildings founded on layered soil media and subjected to vertically incident plane shear waves, energy flux equations are derived by modeling the building as an extension of the layered soil medium, and considering each story as another layer. The propagation of energy flux in the layers is described in terms of the upgoing and downgoing energy flux in each layer, and the energy reflection and transmission coefficients at each interface. The formulation results in a pair of simple finite-difference equations for each layer, which can be solved recursively starting from the bedrock. The upgoing and downgoing energy flux in the layers allows calculation of the energy demand and energy dissipation in each layer. The methodology is applicable to linear, as well as nonlinear structures.     

ELASTIC EXTRAPOLATION OF PRIMARY SEISMIC P- AND S-WAVES1

The elastic Kirchhoff-Helmholtz integral expresses the components of the monochromatic displacement vector at any point A in terms of the displacement field and the stress field at any closed surface surrounding A. By introducing Green's functions for P- and S-waves, the elastic Kirchhoff-Helmholtz integral is modified such that it expresses either the P-wave or the S-wave at A in terms of the elastic wavefield at the closed surface. This modified elastic Kirchhoff-Helmholtz integral is transformed into one-way elastic Rayleigh-type integrals for forward extrapolation of downgoing and upgoing P- and S-waves. We also derive one-way elastic Rayleigh-type integrals for inverse extrapolation of downgoing and upgoing P- and S-waves. The one-way elastic extrapolation operators derived in this paper are the basis for a new prestack migration scheme for elastic data.     

Seismic migration and absorbing boundaries with a one-way wave system for heterogeneous media1

A first-order one-way wave system has been created based on characteristic analysis of the acoustic wave system and optimization of the dispersion relation. We demonstrate that this system is equivalent to a third-order scalar partial-differential equation which, for a homogeneous medium, reduces to a form similar to the 45° paraxial wave equation. This system describes accurately waves propagating in a 2D heterogeneous medium at angles up to 75°. The one-way wave system representing downgoing waves is used for a modified reverse time migration method. As a wavefield extrapolator in migration, the downgoing wave system propagates the reflection events backwards to their reflectors without scattering at the discontinuities in the velocity model. Hence, images with amplitudes proportional to reflectivity can be obtained from this migration technique. We present examples of the application of the new migration method to synthetic seismic data where P-P reflections P-SV converted waves are present. Absorbing boundaries, useful in the generation of synthetic seismograms, have been constructed by using the one-way wave system. These boundaries absorb effectively waves impinging over a wide range of angles of incidence.     

基于两步SVD变换的零偏VSP资料上下行波场分离方法

垂直地震剖面(Vertical Seismic Profiling,VSP)资料处理中波场分离是关键问题之一.随着属性提取技术的发展,新的属性参数(例如Q值)提取技术对波场分离的保真性要求越来越高.本文改进了传统奇异值分解(Singular Value Decomposition,SVD)法,给出了一种对波场的动力学特征具有更好的保真性,可以作为Q值提取的预处理步骤的零偏VSP资料上下行波场分离方法.该方法通过两步奇异值分解变换实现:第一步,排齐下行波同相轴,利用SVD变换压制部分下行波能量;第二步,在剩余波场中排齐上行波同相轴,使用SVD变换提取上行波场.在该方法的实现过程中,压制部分下行波能量后的剩余波场中仍然存在较强的下行波干扰,使得上行波同相轴的排齐比较困难.本文给出了一种通过极大化多道数据线性相关程度(Maximize Coherence,MC)排齐同相轴的算法,在一定程度上解决了低信噪比下排齐同相轴的问题.将本文提出的方法用于合成数据和实际资料的处理,并与传统SVD法的处理结果进行对比,结果表明本文提出的波场分离方法具有良好的保真性,得到波场的质量明显优于传统SVD法.通过对本文方法和传统SVD法处理合成数据得到的下行波场提取Q值,然后进行对比可知,本文方法可以有效提高所提取Q值的准确性,适合作为Q值提取的预处理步骤.     

基于耦合反射/透射系数单程波传播算子的地震波模拟研究

单程波近似实际上是一种多次前向散射和单次后向散射近似.利用单程波近似来描述波传播可以极大地节省地震数值模拟的计算时间和内存,实现地震波长距离传播模拟和三维地震模拟快速计算.本文基于单程波近似和波动积分方程的分离变量逼近,从广义Lippmann-Schwinger波动积分方程推导出耦合反射/透射系数的单程波传播算子.该算子由两部分构成：分离变量Fourier单程波传播算子和薄板间的反射/透射系数表达.前者将常规的Fourier分裂步单程波传播算子(SSF)推广适应横向强速度变化介质和大角度传播波场.后者是利用垂直波数来表示反射/透射系数,自然耦合到波场传播的计算过程中,其为地质界面倾角的隐式表达,精确描述振幅随入射角的变化,能适应任意复杂的模型.通过两个数值算例和一个实际地质模型的计算,本文将该方法和边界元法进行了比较,结果表明：在算例给出的介质横向速度变化情况下,本文提出的方法在相位和振幅方面与全波数值方法基本吻合.     

TOWARDS OPTIMUM ONE-WAY WAVE PROPAGATION1

Numerical wavefield extrapolation represents the backbone of any algorithm for depth migration pre- or post-stack. For such depth imaging techniques to yield reliable and interpretable results, the underlying wavefield extrapolation algorithm must propagate the waves through inhomogeneous media with a minimum of numerically induced distortion, over a range of frequencies and angles of propagation. A review of finite-difference (FD) approximations to the acoustic one-way wave equation in the space-frequency domain is presented. A straightforward generalization of the conventional FD formulation leads to an algorithm where the wavefield is continued downwards with space-variant symmetric convolutional operators. The operators can be precomputed and made accessible in tables such that the ratio between the temporal frequency and the local velocity is used to determine the correct operator at each grid point during the downward continuation. Convolutional operators are designed to fit the desired dispersion relation over a range of frequencies and angles of propagation such that the resulting numerical distortion is minimized. The optimization is constrained to ensure that evanescent energy and waves propagating at angles higher than the maximum design angle are attenuated in each extrapolation step. The resulting operators may be viewed as optimally truncated and bandlimited spatial versions of the familiar phase shift operator. They are unconditionally stable and can be applied explicitly. This results in a simple wave propagation algorithm, eminently suited for implementation on pipelined computers and on large parallel computing systems.     

频率-波数域递推计算海底多分量地震记录中的反射系数

在频率-波数域中采用解析法,解出多层条件下海底实测的多分量地震数据分解成上行和下行P波和S波的算法,导出海底各层地震反射系数随入射角变化(简称RVA)的递推计算公式,为海底多波多分量AVO弹性参数的反演及流体因子预测提供基础数据.合成数据的计算结果表明,本文给出的算法能较可靠地从海底多波多分量记录中提取RVA信息.     

基于局部斜率属性的VSP波场分离研究

基于垂直地震剖面(VSP)资料中上、下行波视速度的差异,利用地震剖面同相轴局部斜率属性参数,提出了一种分离上、下行波场的新方法.首先利用Fourier变换初步分离上、下行波场,然后利用平面波分解滤波器(Plane Wave Destruction (PWD) Filter)技术估计初始分离波场的同相轴局部斜率属性参数,在此基础上对VSP原始资料波场分离.该方法是一种时间域最小平方优化分离波场的方法,不存在其他滤波方法阈值滤波器边界的影响,减少了因镶边问题带来的假象.模拟和实际资料处理结果表明,该方法与传统方法相比,分离出的上、下行波噪声假象少,振幅保持好,更好地消除了上、下行波的相互影响.     

用传播矩阵法研究层状正交各向异性地层中多分量感应测井响应

本文采用传播矩阵技术研究并建立了层状正交各向异性地层中多分量感应测井响应的有效算法.首先通过Fourier变换将频率空间域中的Maxwell方程组求解问题转化为频率波数域中关于电磁场水平分量常微分方程组的定解问题.利用该方程组系数矩阵的本征值和归一化本征向量将电磁场分解成上行波和下行波模式的组合,推导出均匀正交各向异性介质中由任意方向磁偶极子产生的电磁波模式解析表达式;在此基础上,利用叠加原理和边界条件研究了电磁波在层状正交各向异性地层中的反射和透射,给出各个界面上的广义反射系数和不同地层中电磁波振幅的递推公式,进而得到电磁波模式的解析解.为了有效确定频率空间域中的电磁场,采用二维Patterson自适应求积算法结合有限连分式展开技术计算傅氏逆变换.最后通过数值模拟结果证明了该算法的有效性,考察了不同各向异性系数、不同井眼倾角以及仪器长度和工作频率变化等情况下的多分量感应测井响应特征.     

THE SPECTRAL FUNCTION OF A LAYERED SYSTEM AND THE DETERMINATION OF THE WAVEFORMS AT DEPTH*

Consider the mathematical model of a horizontally layered system subject to an initial downgoing source pulse in the upper layer and to the condition that no upgoing waveforms enter the layered system from below the deepest interface. The downgoing waveform (as measured from its first arrival) in each layer is necessarily minimum-phase. The net downgoing energy in any layer, defined as the difference of the energy spectrum of the downgoing wave minus the energy spectrum of the upgoing wave, is itself in the form of an energy spectrum, that is, it is non-negative for all frequencies. The z-transform of the autocorrelation function corresponding to the net downgoing energy spectrum is called the net downgoing spectral function for the layer in question. The net downgoing spectral functions of any two layers A and B are related as follows: the product of the net downgoing spectral function of layer A times the overall transmission coefficient from A to B equals the product of the net downgoing spectral function of layer B times the overall transmission coefficient from B to A. The net downgoing spectral function for the upper layer is called simply the spectral function of the system. In the case of a marine seismogram, the autocorrelation function corresponding to the spectral function can be used to recursively generate prediction error operators of successively increasing lengths, and at the same time the reflection coefficients at successively increasing depths. This recursive method is mathematically equivalent to that used in solving the normal equations in the case of Toeplitz forms. The upgoing wave-form in any given layer multiplied by the direct transmission coefficient from that layer to the surface is equal to the convolution of the corresponding prediction error operator with the surface seismogram. The downgoing waveform in this given layer multiplied by the direct transmission coefficient from that layer to the surface is equal to the convolution of the corresponding hindsight error operator (i.e., the time reverse of the prediction error operator) with the surface seismogram.     

用于零偏移距VSP资料的自适应波形反演方法研究

提出一种利用零偏移距VSP资料初至下行波(即直达波)反演介质品质因子Q及层速度V等参数的方法, 称为自适应时域波形反演法(ATWI). 为了充分地利用有效信息, 该方法根据实际VSP资料的信噪比及直达波与上行波干涉的程度,自适应最大限度地选取未受干扰的初至波片段,并用该片段构造目标函数; 通过恰当地构造数据加权矩阵提高目标函数对Q值变化的敏感性;为克服非线性反演的病态问题,采用近来发展的乘性正则化方法,并通过约束条件限制待求参数的取值范围;文中推导出了雅可比矩阵各元素的解析表达式,从而减小了反问题的计算量.合成数据反演结果表明,与谱比值法和子波包络峰值瞬时频率法相比较,ATWI法受上行波影响相对较小、抗噪性能更强.实际资料算例进一步证明了ATWI方法的有效性.     

VERTICAL SEISMIC PROFILING IN COAL*

The intellection of seismic wave propagation in coal measures demands direct observation of the wavefield progression. Two vertical seismic profiles with high spatial and temporal sampling, were recently recorded in the Sydney Basin coalfields as part of an experimental coal seismic program. Static corrections and interval velocities were obtained by an automated system to determine first kicks and pulse rise times. Upgoing and downgoing waves were separated in the f—k-plane using a novel technique of contour slice filtering. The isolated upgoing waves clearly display reflections from the major coal seams within the stratigraphic sequence. The downgoing wave spectra were subjected to attenuation analysis. The deduced specific quality factor Q for Permian coal measure rocks lies in the range 20–70. Similar estimates were obtained in the time domain from measurements of pulse broadening. Synthetic VSP seismograms, computed using an exact recursive formulation, are an indispensable aid to interpretation. They illustrate the filtering effects of coal seams and sequences, and the effects of the contribution of internal and free-surface multiple reflections in the recorded wavetrains.     

Quadratic normal moveouts of symmetric reflections in elastic media: A quick tutorial

The design of reflection traveltime approximations for optimal stacking and inversion has always been a subject of much interest in seismic processing. A most prominent role is played by quadratic normal moveouts, namely reflection traveltimes around zero-offset computed as second-order Taylor expansions in midpoint and offset coordinates. Quadratic normal moveouts are best employed to model symmetric reflections, for which the ray code in the downgoing direction coincides with the ray code in the upgoing direction in reverse order. Besides pure (non-converted) primaries, many multiply reflected and converted waves give rise to symmetric reflections. We show that the quadratic normal moveout of a symmetric reflection admits a natural decomposition into a midpoint term and an offset term. These, in turn, can be be formulated as the traveltimes of the one-way normal (N) and normal-incidence-point (NIP) waves, respectively. With the help of this decomposition, which is valid for propagation in isotropic and anisotropic elastic media, we are able to derive, in a simple and didactic way, a unified expression for the quadratic normal moveout of a symmetric reflection in its most general form in 3D. The obtained expression allows for a direct interpretation of its various terms and fully encompasses the effects of velocity gradients and Earth surface topography.     

全波震相分析的应用

全波震相分析法是以弹性波传播理论为基础,以现代数字技术为手段,对地震波进行全面综合分析的理论和方法.它既包括对单一分量震波观测记录中的各种震相的识别、确认和分析,又包括对多分量震波记录的合成、分析和图示.全波震相分析法涉及震源、射线路径、初动、走时、振幅、波形、时域和空域的瞬态谱、质点振动矢量等诸多方面,因而能更全面地揭示地震波场与地下介质之间的关系.本文以实例展示在不同领域的研究中全波震相分析的应用情况.全波震相分析法为实现多种波型的联合应用奠定了理论基础,并提供了新的方法.     

The One-Way Wave Equation: A Full-Waveform Tool for Modeling Seismic Body Wave Phenomena

The study of seismic body waves is an integral aspect in global, exploration and engineering scale seismology, where the forward modeling of waves is an essential component in seismic interpretation. Forward modeling represents the kernel of both migration and inversion algorithms as the Green’s function for wavefield propagation and is also an important diagnostic tool that provides insight into the physics of wave propagation and a means of testing hypotheses inferred from observational data. This paper introduces the one-way wave equation method for modeling seismic wave phenomena and specifically focuses on the so-called operator-root one-way wave equations. To provide some motivation for this approach, this review first summarizes the various approaches in deriving one-way approximations and subsequently discusses several alternative matrix narrow-angle and wide-angle formulations. To demonstrate the key strengths of the one-way approach, results from waveform simulation for global scale shear-wave splitting modeling, reservoir-scale frequency-dependent shear-wave splitting modeling and acoustic waveform modeling in random heterogeneous media are shown. These results highlight the main feature of the one-way wave equation approach in terms of its ability to model gradual vector (for the elastic case) and scalar (for the acoustic case) waveform evolution along the underlying wavefront. Although not strictly an exact solution, the one-way wave equation shows significant advantages (e.g., computational efficiency) for a range of transmitted wave three-dimensional global, exploration and engineering scale applications.     

Joint viscoacoustic waveform inversion with the separated upgoing and downgoing wavefields in zero-offset vertical seismic profile data

The attenuation of seismic waves propagating in reservoirs can be obtained accurately from the data analysis of vertical seismic profile in terms of the quality-factor Q. The common methods usually use the downgoing wavefields in vertical seismic profile data. However, the downgoing wavefields consist of more than 90% energy of the spectrum of the vertical seismic profile data, making it difficult to estimate the viscoacoustic parameters accurately. Thus, a joint viscoacoustic waveform inversion of velocity and quality-factor is proposed based on the multi-objective functions and analysis of the difference between the results inverted from the separated upgoing and downgoing wavefields. A simple separating step is accomplished by the reflectivity method to obtain the individual wavefields in vertical seismic profile data, and then a joint inversion is carried out to make full use of the information of the individual wavefields and improve the convergence of viscoacoustic full-waveform inversion. The sensitivity analysis of the different wavefields to the velocity and quality-factor shows that the upgoing and downgoing wavefields contribute differently to the viscoacoustic parameters. A numerical example validates our method can improve the accuracy of viscoacoustic parameters compared with the direct inversion using full wavefield and the separate inversion using upgoing or downgoing wavefield. The application on real field data indicates our method can recover a reliable viscoacoustic model, which helps reservoir appraisal.