基于高斯束与高斯波包的Gabor框架散射波模拟方法

在给出真实模型和相应光滑背景模型的情况下,如何计算扰动模型(散射体)产生的散射波场是一个有实际意义的正演问题.在Gabor变换域描述散射体,且入射波场为短时宽带信号时,散射波场可以在频率域用高斯束或时间域用高斯波包描述.相对于波动方程方法,高斯束和高斯波包的计算效率更高;背景模型光滑时,高斯束和高斯波包方法的精度也接近波动方程方法.文中导出了声波假设下应用高斯束和高斯波包计算散射波的方法.测试分析了高斯波包的计算精度.给出了一般散射体的散射波模拟策略.同时针对一个理论模型完成了本文方法计算散射波的实验,实验结果表明高斯波包散射波计算方法是有效可行的.     

高斯波包反射走时速度反演方法

扰动高斯波包理论指出,在Gabor域描述模型的扰动成分,且入射波场为短时宽带信号时,扰动波场可在时间域通过高斯波包算子描述.在此基础上通过拟合反射波的走时,提出一种速度反演方法.反射波走时残差利用地震道局部波形的互相关函数表示,以走时残差的二范数作为目标函数,优化目标函数实现对速度场的反演.基于一阶Born近似,利用扰动高斯波包理论推导出目标函数对速度场的梯度是本文理论部分的核心内容.梯度包括两部分:正传的背景波场与反传的扰动高斯波包之间的互相关,反传的背景波场和正传的扰动高斯波包之间的互相关.梯度表达式中背景波场和扰动波场均利用高斯波包算子模拟.计算梯度的具体算法中,如何模拟扰动波场,以及如何计算反射波的走时残差是两个要点,文中对此做了详细的讨论.数值实验进一步阐述了反演的实现策略,实验结果表明高斯波包反射走时速度反演方法和实现策略有效可行,并得到了理想的反演结果.     

Gabor‐frame‐based Gaussian packet migration

We present a Gaussian packet migration method based on Gabor frame decomposition and asymptotic propagation of Gaussian packets. A Gaussian packet has both Gaussian‐shaped time–frequency localization and space–direction localization. Its evolution can be obtained by ray tracing and dynamic ray tracing. In this paper, we first briefly review the concept of Gaussian packets. After discussing how initial parameters affect the shape of a Gaussian packet, we then propose two Gabor‐frame‐based Gaussian packet decomposition methods that can sparsely and accurately represent seismic data. One method is the dreamlet–Gaussian packet method. Dreamlets are physical wavelets defined on an observation plane and can represent seismic data efficiently in the local time–frequency space–wavenumber domain. After decomposition, dreamlet coefficients can be easily converted to the corresponding Gaussian packet coefficients. The other method is the Gabor‐frame Gaussian beam method. In this method, a local slant stack, which is widely used in Gaussian beam migration, is combined with the Gabor frame decomposition to obtain uniform sampled horizontal slowness for each local frequency. Based on these decomposition methods, we derive a poststack depth migration method through the summation of the backpropagated Gaussian packets and the application of the imaging condition. To demonstrate the Gaussian packet evolution and migration/imaging in complex models, we show several numerical examples. We first use the evolution of a single Gaussian packet in media with different complexities to show the accuracy of Gaussian packet propagation. Then we test the point source responses in smoothed varying velocity models to show the accuracy of Gaussian packet summation. Finally, using poststack synthetic data sets of a four‐layer model and the two‐dimensional SEG/EAGE model, we demonstrate the validity and accuracy of the migration method. Compared with the more accurate but more time‐consuming one‐way wave‐equation‐based migration, such as beamlet migration, the Gaussian packet method proposed in this paper can correctly image the major structures of the complex model, especially in subsalt areas, with much higher efficiency. This shows the application potential of Gaussian packet migration in complicated areas.     

Resolution of prestack depth migration

The resolution of a general 3-D common-shot elastic prestack depth migration in a heterogeneous anisotropic medium is studied approximately, using the ray theory. It is demonstrated that the migrated section can approximately be represented by the convolution of the reflectivity function with the corresponding local resolution function. Alternatively, it can also be approximately represented by the convolution of the spatial distribution of the weak-contrast displacement reflection-transmission coefficient with the corresponding local resolution function. The derived explicit approximate equations enable us to predict the migration resolution approximately without doing the whole and expensive migration. The equations are applicable to 3-D elastic migrations in 3-D isotropic or anisotropic, heterogeneous velocity models. Both the reflectivity function and the spatial distribution of the weak-contrast displacement reflection-transmission coefficient approximately determine the linear combination of the perturbations of elastic moduli and density to which the migrated section is sensitive. The imaged linear combination of the perturbations of elastic parameters depends on the selection of the polarizations (wave types) of the incident and back-propagated wavefields and on the directions of propagation. The resolution of the linear combination of the perturbations of elastic moduli and density in the migrated section is determined by the above mentioned local resolution functions. The local resolution functions depend on the aperture and on the imaging function. The imaging function is determined by the source time function and by the form of the imaging functional. The local resolution functions are considerably sensitive to heterogeneities. The local resolution functions in elastic media differ from their acoustic counterparts, especially by the existence of converted scattered waves in elastic media.     

特征高斯波包叠前深度偏移方法

高斯波包（Gaussian packet）传播算子可在局部时空域高效地计算局部波包的传播.高斯波包叠前深度偏移的基础是在Gabor变换域描述观测数据,再利用高斯波包传播算子计算炮点波场和检波点波场,两者相关即可得到偏移结果.利用炮道集的局部τ-p特征在Gabor变换域表达观测数据,可以仅关注部分高斯波包框架函数上的数据投影,这样既实现了波场的压缩存储,同时可利用高斯波包传播算子反传框架函数以实现整个炮道集的快速反传.这些综合了观测数据局部τ-p特征的高斯波包函数称为特征高斯波包（characteristic Gaussian packet,CGP）,相应的波场反传称为特征高斯波包反传.理论及数值分析证明了上述特征高斯波包反传方法是有效且快速的.炮点正传波场也利用高斯波包传播算子模拟.利用互相关成像条件可实现特征高斯波包叠前深度偏移（characteristic Gaussian packet pre-stack depth migration,CGPM）.由于高斯波包传播算子描述了局部方向及局部空间的波场,所以CGPM可以自然地提取角度域成像道集（ADCIG）,并易于实现面向目标叠前深度偏移,从而作为偏移引擎为偏移速度分析（MVA）服务.数值实验证明了CGPM和面向目标CGPM的有效性和实用性.     

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.     

Elastic wave‐vector decomposition in heterogeneous anisotropic media

The goal of wave‐mode separation and wave‐vector decomposition is to separate a full elastic wavefield into three wavefields with each corresponding to a different wave mode. This allows elastic reverse‐time migration to handle each wave mode independently. Several of the previously proposed methods to accomplish this task require the knowledge of the polarisation vectors of all three wave modes in a given anisotropic medium. We propose a wave‐vector decomposition method where the wavefield is decomposed in the wavenumber domain via the analytical decomposition operator with improved computational efficiency using low‐rank approximations. The method is applicable for general heterogeneous anisotropic media. To apply the proposed method in low‐symmetry anisotropic media such as orthorhombic, monoclinic, and triclinic, we define the two S modes by sorting them based on their phase velocities (S1 and S2), which are defined everywhere except at the singularities. The singularities can be located using an analytical condition derived from the exact phase‐velocity expressions for S waves. This condition defines a weight function, which can be applied to attenuate the planar artefacts caused by the local discontinuity of polarisation vectors at the singularities. The amplitude information lost because of weighting can be recovered using the technique of local signal–noise orthogonalisation. Numerical examples show that the proposed approach provides an effective decomposition method for all wave modes in heterogeneous, strongly anisotropic media.     

True amplitude imaging by inverse generalized Radon transform based on Gaussian beam decomposition of the acoustic Green's function

True amplitude migration is one of the most important procedures of seismic data processing. As a rule it is based on the decomposition of the velocity model of the medium into a known macrovelocity component and its sharp local perturbations to be determined. Under this decomposition the wavefield can be considered as the superposition of an incident and reflected/scattered waves. The single scattering approximation introduces the linear integral operator that connects the sharp local perturbations of the macrovelocity model with the multishot/multioffset data formed from reflected/scattered waves. We develop the pseudoinverse of this operator using the Gaussian beam based decomposition of acoustic Green's functions. The computation of this pseudoinverse operator is done pointwise by shooting Gaussian beams from the target area towards the acquisition system. The numerical implementation of the pseudoinverse operator was applied to the synthetic data Sigsbee2A. The results obtained demonstrate the high quality of the true amplitude images computed both in the smooth part of the model and under the salt body.     

Analysis of the resolution function in seismic prestack depth imaging

We consider the problem of estimating subsurface quantities such as velocity or reflectivity from seismic measurements. Because of a limited aperture and band-limited signals, the output from a seismic prestack reconstruction method is a distorted or blurred image. This distortion can be computed using the concept of resolution function, which is a quantity readily accessible in the Fourier space of the model. The key parameter is the scattering wavenumber, which at a particular image point is defined by the incident and scattered ray directions in a given background model. Any location in any background model can be considered. In general, the resolution function will depend on the following four quantities: the background velocity model, the frequency bandwidth, the wavefield type and the acquisition geometry.
We first establish the resolution function for a general scattering model assuming local reaction. We then adapt this result for two well-known scattering models: Born and Kirchhoff. For each of these approximations the corresponding resolution function is derived and discussed. Finally, by employing a simple synthetic data example we demonstrate the ability of the resolution function to predict the image distortions.     

Decomposition of the wave field into optimized Gaussian packets

The decomposition of the wave field into optimized Gaussian packets represents a crucial step of the Gaussian packet prestack depth migration algorithm. The shape of optimized Gaussian packets, in the plane perpendicular to the central ray of the packet, depends not only on the frequency, but also on the coordinate of the intersection of the central ray of a Gaussian packet with the profile, on its arrival time, and on the component of the slowness vector along the profile. We express the amplitude of a Gaussian packet in the form of an integral transform similar to the forward coherentstate transform. Our method is suitable for a smooth distribution of the parameter determining the shape of a packet in the plane perpendicular to its central ray.     

Calibration of cross-line components for sea-bed 4C acquisition systems

Faithful recording of the elastic wavefield at the sea‐bed is required for quantitative applications of 4C seismic. The accuracy of the recorded vectorial wavefield depends on factors that vary from deployment to deployment. This paper focuses on one such factor: the interaction of the acquisition system with the sea‐bed, which is referred to here as coupling. We show, using multi‐azimuth data recorded with a cable‐based sea‐bed acquisition system, whose sensor housing is cylindrically shaped and with the in‐line geophone fixed to the cable, that coupling depends on the propagation direction and wave type (P‐ or S‐waves) of the incident wavefield. We show that coupling is more critical for S‐waves than for P‐waves. Detection of inconsistent coupling using both P‐ and S‐waves is therefore mandatory. A data‐driven processing method to compensate for the frequency‐dependent coupling response of the cross‐line geophone is derived. Its application to field data verifies the effectiveness of the method.     

Regional Seismic Wavefield Computation on a 3-D Heterogeneous Earth Model by Means of Coupled Traveling Wave Synthesis

-- I present a new algorithm for calculating seismic wave propagation through a three-dimensional heterogeneous medium using the framework of mode coupling theory originally developed to perform very low frequency (f < ~0.01т.05 Hz) seismic wavefield computation. It is a Greens function approach for multiple scattering within a defined volume and employs a truncated traveling wave basis set using the locked mode approximation. Interactions between incident and scattered wavefields are prescribed by mode coupling theory and account for the coupling among surface waves, body waves, and evanescent waves. The described algorithm is, in principle, applicable to global and regional wave propagation problems, but I focus on higher frequency (typically f S ~0.25 Hz) applications at regional and local distances where the locked mode approximation is best utilized and which involve wavefields strongly shaped by propagation through a highly heterogeneous crust. Synthetic examples are shown for P-SV-wave propagation through a semi-ellipsoidal basin and SH-wave propagation through a fault zone.     

Imaging Offsets in the Moho: Synthetic Tests using Gaussian Beams with Teleseismic Waves

We carry out a sequence of numerical tests to understand conditions under which rapid changes in crustal thickness can be reliably imaged by teleseismic body waves. Using the finite-difference method over a 2-D grid, we compute synthetic seismograms resulting from a planar P-wavefield incident below the grid. We then image the Moho using a migration scheme based on the Gaussian beam representation of the wavefield. The use of Gaussian beams for the downward propagation of the wavefield is particularly advantageous in certain geologically critical cases such as overthrusting of continental lithosphere, resulting in the juxtaposition of high-velocity mantle material over crustal rocks. In contrast to ray-based methods, Gaussian beam migration requires no special treatment to handle such heterogeneities. Our results suggest that with adequate station spacing and signal-to-noise ratios, offsets of the Moho, on the order of 10 km in height, can be reliably imaged beneath thickened crust at depths of about 50 km. Furthermore, even sharp corners and edges are faithfully imaged when precise values of seismic wave speeds are available. Our tests also demonstrate that flexibility in choices of different types of seismic phases is important, because any single phase has trade-offs in issues such as spatial resolution, array aperture, and amplitude of signals.     

Scattering of plane harmonic SH,SV, P and rayleigh waves by non-axisymmetric three-dimensional canyons: A wave function expansion approach

Scattering of elastic waves by three-dimensional canyons embedded within an elastic half-space is investigated by using a wave function expansion technique. The geometry of the canyon is assumed to be non-axisymmetric. The canyon is subjected to incident plane Rayleigh waves and oblique incident SH, SV and P waves. The unknown scattered wavefield is expressed in terms of spherical wave functions which satisfy the equations of motion and radiation conditions at infinity, but they do not satisfy stress-free boundary conditions at the half-space surface. The boundary conditions are imposed locally in the least-squares sense at several points on the surface of the canyon and the half-space. Through a comparative study the validity and limitations of two-dimensional approximations (antiplane strain and plane strain models) have been examined. It is shown that scattering of waves by three-dimensional canyons may cause substantial change in the surface displacement patterns in comparison to the two-dimensional models. These results emphasize the need for three-dimensional modelling of realistic problems of interest in strong ground motion seismology and earthquake engineering.     

各向异性半空间中浅埋孔洞对地表 反平面运动的影响

本文采用复变函数方法研究了稳态水平剪切波(SH波)在各向异性弹性半空间中任意形状孔洞上的散射及其对地面运动的影响.在变换空间中构造出自动满足各向异性半空间水平表面上应力为零的散射波函数,并根据孔洞表面应力为零的边界条件,采用最小二乘法求解散射波函数的系数.用介质的各向异性性质来模拟地质条件,给出了SH波作用下含圆孔、椭圆孔和方孔的各向异性半空间表面位移幅值的数值结果,并分析了介质特性、孔洞的形状、埋深、入射波波数及入射角度等因素对地面运动的影响规律.数值结果表明:介质的各向异性对含有孔洞的半空间表面的地表位移具有显著的影响;沿一定角度的入射波在某一频段内所引起的地表位移幅值比各向同性介质的可能要大,且随着孔洞埋深的增加,地表位移的幅值逐渐减小.      

中层耗散大气中Rossby波包的非线性相互作用理论

采用弱非线性近似得出中层耗散大气连续谱Ｒｏｓｓｂｙ波包的非线性时空演化方程，讨论了Ｒｏｓｓｂｙ波包的三波相互作用问题．数值计算表明，耗散和非线性的共同效应决定了Ｒｏｓｓｂｙ波包的演变．当一个Ｒｏｓｓｂｙ波包通过大气传播时，它的振幅若超过某个阈值，空间尺度分别比它大和比它小的两个次级Ｒｏｓｓｂｙ波包的振幅会随时间增长．特别当这两个次级波包同时随时空变化时，仅当主波的振幅超过一个更大的阈值，且其群速度介于两次级波包的群速度之间时，两次级波包的振幅才会随时空同时增长，即出现绝对不稳定现象，耗散和３个波包的频率失配都会增大不稳定的阈值．     

VTI介质角度域转换波高斯束偏移成像方法研究

转换波偏移可以利用纵横波波场信息,得到高分辨率的成像结果,从而为油藏描述提供高质量的地震资料.目前的研究主要是利用纵波波场信息进行偏移成像,然而,传统的纵波方法在复杂探区成像时具有一定的局限性.为此,本文在各向异性介质声波射线追踪算法的基础上,推导出各向异性介质转换波射线追踪方程,发展了一种转换波射线追踪算法;并将研究的追踪算法应用到偏移成像中,提出了一种各向异性VTI介质角度域转换波高斯束偏移成像方法.通过各向异性VTI介质断块模型和复杂构造模型试算,说明了本文方法的正确性和有效性.模型试算的结果表明,在考虑地下各向异性时,本文研究的方法具有更好的成像效果,提取的角道集结果可以为偏移速度分析提供依据.     

正交各向异性介质中qP波入射的二阶近似反射系数与透射系数

地震波在各向异性介质中以一个准P波(qP)和两个准S波(qS1和qS2)的形式传播.研究三种波的相速度、群速度以及偏振方向等传播性质能够为各向异性介质中的正反演问题提供有效支撑.具有比横向各向同性(TI)介质更一般对称性的正交各向异性介质通常需要9个独立参数对其进行描述,这使得对传播特征的计算更为复杂.当两个准S波速度相近时具有耦合性,从而令慢度的计算产生奇异性.因此,奇异点(慢度面的鞍点和交叉点)附近的反射与透射(R/T)系数的求解不稳定,会导致波场振幅不准确.本文首次通过结合耦合S波射线理论和基于迭代的各向异性相速度与偏振矢量的高阶近似解,得到了适用于正交各向异性介质以qP波入射所产生的二阶R/T系数的计算方法.与基于一阶近似的结果相比,基于二阶近似的方法提高了qP波R/T系数的精度,能得到一阶耦合近似无法表达的准确的qP-qS转换波的R/T系数解,且方法适用于较强的各向异性介质.     

基于delta波包叠加的时间域深度偏移

delta波包可由高斯波束经傅里叶逆变换得到,是高斯波束在时空域的表达.它最早出现在合成理论地震图的研究中,本文将其应用于偏移领域.通过delta波包叠加表达时间域格林函数,可将高斯波束偏移由频率域转换到时间域,再结合Rayleigh积分和激励时间成像条件,本文给出了基于delta波包叠加的深度偏移算法.该偏移算法可在时间域直接计算,但因包含褶积运算,成像时将耗费大量的计算时间.针对这一问题,本文提出了将褶积简化为乘积的近似公式.近似后的偏移算法,不仅保留了高斯波束偏移的优点,而且计算效率得到显著提升.文中通过两个数值算例验证了上述结论.     

Earthquake site effect modeling in the Granada basin using a 3-D indirect boundary element method

In this work, we develop the indirect boundary element method (IBEM) to simulate the seismic site response in a realistic, large-scale 3-D sedimentary basin. Most previous applications of boundary element method have used full-space Green’s functions for wave propagation between element points. We use half-space Green’s functions, which include the seismic wavefield interactions at the free surface and require only the boundary elements of the basin interface. In this way, the size of the matrix equation for solution in the IBEM can be reduced to approximately a quarter of that using full-space Green’s functions. The site response modeling of the Granada basin in southern Spain using the IBEM shows that the basin-induced scattering waves were identified as propagating back and forth inside the basin. The scattered waves also generate surface waves that are weakly propagated outside of the basin. The wave propagation inside and outside of the basin shows different patterns. We observe that the scattered wave is locally amplified, and its propagation direction deviates from that of the incident waves propagation direction. Therefore, the computed seismic response in the basin could provide us with good estimates of the seismic motion.