Wave-equation migration velocity analysis. I. Theory

We present a migration velocity analysis (MVA) method based on wavefield extrapolation. Similarly to conventional MVA, our method aims at iteratively improving the quality of the migrated image, as measured by the flatness of angle‐domain common‐image gathers (ADCIGs) over the aperture‐angle axis. However, instead of inverting the depth errors measured in ADCIGs using ray‐based tomography, we invert ‘image perturbations’ using a linearized wave‐equation operator. This operator relates perturbations of the migrated image to perturbations of the migration velocity. We use prestack Stolt residual migration to define the image perturbations that maximize the focusing and flatness of ADCIGs. Our linearized operator relates slowness perturbations to image perturbations, based on a truncation of the Born scattering series to the first‐order term. To avoid divergence of the inversion procedure when the velocity perturbations are too large for Born linearization of the wave equation, we do not invert directly the image perturbations obtained by residual migration, but a linearized version of the image perturbations. The linearized image perturbations are computed by a linearized prestack residual migration operator applied to the background image. We use numerical examples to illustrate how the backprojection of the linearized image perturbations, i.e. the gradient of our objective function, is well behaved, even in cases when backprojection of the original image perturbations would mislead the inversion and take it in the wrong direction. We demonstrate with simple synthetic examples that our method converges even when the initial velocity model is far from correct. In a companion paper, we illustrate the full potential of our method for estimating velocity anomalies under complex salt bodies.     

Directional‐oriented wavefield imaging: a new wave‐based subsurface illumination imaging condition for reverse time migration

The key objective of an imaging algorithm is to produce accurate and high‐resolution images of the subsurface geology. However, significant wavefield distortions occur due to wave propagation through complex structures and irregular acquisition geometries causing uneven wavefield illumination at the target. Therefore, conventional imaging conditions are unable to correctly compensate for variable illumination effects. We propose a generalised wave‐based imaging condition, which incorporates a weighting function based on energy illumination at each subsurface reflection and azimuth angles. Our proposed imaging kernel, named as the directional‐oriented wavefield imaging, compensates for illumination effects produced by possible surface obstructions during acquisition, sparse geometries employed in the field, and complex velocity models. An integral part of the directional‐oriented wavefield imaging condition is a methodology for applying down‐going/up‐going wavefield decomposition to both source and receiver extrapolated wavefields. This type of wavefield decomposition eliminates low‐frequency artefacts and scattering noise caused by the two‐way wave equation and can facilitate the robust estimation for energy fluxes of wavefields required for the seismic illumination analysis. Then, based on the estimation of the respective wavefield propagation vectors and associated directions, we evaluate the illumination energy for each subsurface location as a function of image depth point and subsurface azimuth and reflection angles. Thus, the final directional‐oriented wavefield imaging kernel is a cross‐correlation of the decomposed source and receiver wavefields weighted by the illuminated energy estimated at each depth location. The application of the directional‐oriented wavefield imaging condition can be employed during the generation of both depth‐stacked images and azimuth–reflection angle‐domain common image gathers. Numerical examples using synthetic and real data demonstrate that the new imaging condition can properly image complex wave paths and produce high‐fidelity depth sections.     

Wave‐equation migration velocity analysis with time‐shift imaging

Wave‐equation migration velocity analysis is a technique designed to extract and update velocity information from migrated images. The velocity model is updated through the process of optimizing the coherence of images migrated with the known background velocity model. The capacity for handling multi‐pathing of the technique makes it appropriate in complex subsurface regions characterized by strong velocity variation. Wave‐equation migration velocity analysis operates by establishing a linear relation between a slowness perturbation and a corresponding image perturbation. The linear relationship and the corresponding linearized operator are derived from conventional extrapolation operators and the linearized operator inherits the main properties of frequency‐domain wavefield extrapolation. A key step in the implementation is to design an appropriate procedure for constructing an image perturbation relative to a reference image that represents the difference between the current image and a true, or more correct image of the subsurface geology. The target of the inversion is to minimize such an image perturbation by optimizing the velocity model. Using time‐shift common‐image gathers, one can characterize the imperfections of migrated images by defining the focusing error as the shift of the focus of reflections along the time‐shift axis. The focusing error is then transformed into an image perturbation by focusing analysis under the linear approximation. As the focusing error is caused by the incorrect velocity model, the resulting image perturbation can be considered as a mapping of the velocity model error in the image space. Such an approach for constructing the image perturbation is computationally efficient and simple to implement. The technique also provides a new alternative for using focusing information in wavefield‐based velocity model building. Synthetic examples demonstrate the successful application of our method to a layered model and a subsalt velocity update problem.     

双平方根波动方程偏移速度分析

传统的剩余校正（RMO）偏移速度分析方法基于走时原理,在陡倾角和欠照明地区,因为不能得到充分的角度域信息而失效.本文将展示一种基于波场延拓理论的偏移速度分析方法,即波动方程偏移速度分析（WEMVA）.这种方法先利用成像优化方法获得剩余成像,再利用剩余成像反演剩余速度.此类方法继承了波动方程偏移方法的优点和缺点.波动方程偏移速度分析是一种线性反演方法,它要求对Born近似的展开序列作一阶截断.高阶部分的丢失必然带来巨大的截断误差,因此剩余成像必须也进行线性化,以适应大速度扰动和大延拓步长.因此,在此类算法中,剩余成像的获取和线性化是偏移速度分析的关键.在叠前偏移算子中,因为双平方根算子的数学表达式更为简洁,所以本文基于对波动方程偏移速度分析初步讨论,并通过模型验证其原理.     

时间二阶积分波场的全波形反演

通过对波场的时间二阶积分运算以增强地震数据中的低频成分,提出了一种可有效减小对初始速度模型依赖性的地震数据全波形反演方法—时间二阶积分波场的全波形反演方法.根据散射理论中的散射波场传播方程,推导出时间二阶积分散射波场的传播方程,再利用一阶Born近似对时间二阶积分散射波场传播方程进行线性化.在时间二阶积分散射波场传播方程的基础上,利用散射波场反演地下散射源分布,再利用波场模拟的方法构建地下入射波场,然后根据时间二阶积分散射波场线性传播方程中散射波场与入射波场、速度扰动间的线性关系,应用类似偏移成像的公式得到速度扰动的估计,以此建立时间二阶积分波场的全波形迭代反演方法.最后把时间二阶积分波场的全波形反演结果作为常规全波形反演的初始模型可有效地减小地震波场全波形反演对初始模型的依赖性.应用于Marmousi模型的全频带合成数据和缺失4Hz以下频谱成分的缺低频合成数据验证所提出的全波形反演方法的正确性和有效性,数值试验显示缺失4Hz以下频谱成分数据的反演结果与全频带数据的反演结果没有明显差异.     

Strategies for stable attenuation compensation in reverse‐time migration

Attenuation in seismic wave propagation is a common cause for poor illumination of subsurface structures. Attempts to compensate for amplitude loss in seismic images by amplifying the wavefield may boost high‐frequency components, such as noise, and create undesirable imaging artefacts. In this paper, rather than amplifying the wavefield directly, we develop a stable compensation operator using stable division. The operator relies on a constant‐Q wave equation with decoupled fractional Laplacians and compensates for the full attenuation phenomena by performing wave extrapolation twice. This leads to two new imaging conditions to compensate for attenuation in reverse‐time migration. A time‐dependent imaging condition is derived by applying Q‐compensation in the frequency domain, whereas a time‐independent imaging condition is formed in the image space by calculating image normalisation weights. We demonstrate the feasibility and robustness of the proposed methods using three synthetic examples. We found that the proposed methods are capable of properly compensating for attenuation without amplifying high‐frequency noise in the data.     

Full‐wavefield migration: using surface and internal multiples in imaging

Multiple scattering is usually ignored in migration algorithms, although it is a genuine part of the physical reflection response. When properly included, multiples can add to the illumination of the subsurface, although their crosstalk effects are removed. Therefore, we introduce full‐wavefield migration. It includes all multiples and transmission effects in deriving an image via an inversion approach. Since it tries to minimize the misfit between modeled and observed data, it may be considered a full waveform inversion process. However, full‐wavefield migration involves a forward modelling process that uses the estimated seismic image (i.e., the reflectivities) to generate the modelled full wavefield response, whereas a smooth migration velocity model can be used to describe the propagation effects. This separation of modelling in terms of scattering and propagation is not easily achievable when finite‐difference or finite‐element modelling is used. By this separation, a more linear inversion problem is obtained. Moreover, during the forward modelling, the wavefields are computed separately in the incident and scattered directions, which allows the implementation of various imaging conditions, such as imaging reflectors from below, and avoids low‐frequency image artefacts, such as typically observed during reverse‐time migration. The full wavefield modelling process also has the flexibility to image directly the total data (i.e., primaries and multiples together) or the primaries and the multiples separately. Based on various numerical data examples for the 2D and 3D cases, the advantages of this methodology are demonstrated.     

基于局部余弦基小波束的观测系统沉降法叠前深度偏移

将局部余弦基小波束波场分解、传播与观测系统沉降法叠前深度偏移相结合,推导了源-检波器观测系统沉降法传播算子.本算法中,先对频率域的共点源和共点检波器道集做局部余弦小波束分解,然后分别沿共小波束源和共小波束检波器在深度方向延拓得到下一层波场.每个深度的波场,都等效于把源和检波器放在该层后所能接收到的地震记录,每点的像值由炮点和检波点重合时的零时刻波场值给出.通过二维SEG/EAGE盐丘模型和Marmousi模型的偏移成像结果验证该方法理论推导的正确性.另外,结果显示该方法继承了小波束域波场延拓在速度扰动较大情况下波传播及成像精度高的优点.     

Increasing illumination and sensitivity of reverse‐time migration with internal multiples

Reverse‐time migration is a two‐way time‐domain finite‐frequency technique that accurately handles the propagation of complex scattered waves and produces a band‐limited representation of the subsurface structure that is conventionally assumed to be linear in the contrasts in model parameters. Because of this underlying linear single‐scattering assumption, most implementations of this method do not satisfy the energy conservation principle and do not optimally use illumination and model sensitivity of multiply scattered waves. Migrating multiply scattered waves requires preserving the non‐linear relation between the image and perturbation of model parameters. I modify the extrapolation of source and receiver wavefields to more accurately handle multiply scattered waves. I extend the concept of the imaging condition in order to map into the subsurface structurally coherent seismic events that correspond to the interaction of both singly and multiply scattered waves. This results in an imaging process referred to here as non‐linear reverse‐time migration. It includes a strategy that analyses separated contributions of singly and multiply scattered waves to a final non‐linear image. The goal is to provide a tool suitable for seismic interpretation and potentially migration velocity analysis that benefits from increased illumination and sensitivity from multiply scattered seismic waves. It is noteworthy that this method can migrate internal multiples, a clear advantage for imaging challenging complex subsurface features, e.g., in salt and basalt environments. The results of synthetic seismic imaging experiments, including a subsalt imaging example, illustrate the technique.     

扩展成像条件下的最小二乘逆时偏移

逆时偏移(RTM)是复杂介质条件下地震成像的重要手段.因受观测系统限制、上覆地层影响以及波场带宽有限等因素的影响,现行的常规RTM所采用的互相关成像条件通常对地下构造进行模糊成像.最小二乘逆时偏移(LSRTM)通过最小化线性Born近似正演数据和采集数据之间的波形差异,采用梯度类反演算法优化反射系数模型,获得的成像结果具有更高的分辨率和更可靠的振幅保真度.然而,基于波形拟合的LSRTM对背景速度模型的依赖性很强.误差太大的速度模型容易产生周波跳跃现象,导致LSRTM难以获得全局最优解.为了克服这一问题,本文基于扩展模型的思想,在线性Born近似下,推导得到RTM扩展成像条件.并基于最小二乘反演理论,提出扩展成像条件下的LSRTM方法.理论模型试算表明,本文方法不仅可以提供分辨率更高、振幅属性更为可靠的成像结果,而且能够在一定程度上消除速度误差对反演成像的影响.     

基于高斯束传播算子的成像域走时层析成像方法

层析反演是速度建模中最重要的方法之一,结合偏移成像在成像域进行走时层析速度反演是当前比较成熟有效且广泛应用的技术.本文从高斯束偏移成像条件出发,在波动方程的一阶Born近似和Rytov近似下,推导了成像域走时扰动与速度扰动的线性关系,建立了成像域走时层析方程及其显式表达的层析核函数.该核函数的本质是有限频层析核函数,利用该核函数替换常规射线层析核函数可以明显提高层析反演精度.该核函数的计算关键是背景波场格林函数的计算,本文利用高斯束传播算子计算格林函数进而得到走时层析核函数,实现方式灵活高效且计算精度较高.基于高斯束传播算子的偏移成像与层析成像相结合进行深度域建模迭代,体现了速度建模与偏移成像一体化的思想.数值计算及实际数据应用证明了基于高斯束传播算子的成像域走时层析方法的有效性.     

The influence of sea water velocity variation on seismic traveltimes, ray paths, and amplitude

The main factors affecting seismic exploration is the propagation velocity of seismic waves in the medium. In the past, during marine seismic data processing, the propagation velocity of sea water was generally taken as a constant 1500 m/s. However, for deep water exploration, the sound velocity varies with the season, time, location, water depth, ocean currents, and etc.. It also results in a layered velocity distribution, so there is a difference of seismic traveltime, ray paths, and amplitude, which affect the migration imaging results if sea water propagation velocity is still taken as constant for the propagation wavefield. In this paper, we will start from an empirical equation of seismic wave velocity in seawater with changes of temperature, salinity, and depth, consider the variation of their values, build a seawater velocity model, and quantitatively analyze the impact of seawater velocity variation on seismic traveltime, ray paths, and amplitude in the seawater velocity model.     

地震多次波成像

地震资料含有各种类型多次波,而传统成像方法仅利用地震一次反射波成像,在地震成像前需将多次波去除.然而,多次波携带了丰富的地下结构信息,多次波偏移能够提供除反射波外的额外地下照明.修改传统逆时偏移方法,用包含一次反射波和多次波的原始记录代替震源子波,将SRME方法预测的表面多次波代替一次反射波作为输入数据,可将表面多次波成像.多次波成像的挑战和困难在于大量串扰噪声的产生,针对表面多次波成像中的成像噪声问题,将最小二乘逆时偏移方法与多次波分阶思想结合起来,发展可控阶数的表面多次波反演成像方法,有望初步实现高精度的表面多次波成像.在消除原始记录中的表面多次波后,通过逆散射级数方法预测得到层间多次波,将层间多次波作为逆时偏移方法的输入数据可将其准确归位到地下反射位置.数值实验表明,多次波成像能够有效地为地下提供额外照明,而可控阶表面多次波最小二乘逆时偏移成像方法几乎完全避免成像噪声.     

Reverse time migration of ground-penetrating radar with full wavefield decomposition based on the Hilbert transform

The conventional reverse time migration of ground-penetrating radar data is implemented with the two-way wave equation. The cross-correlation result contains low-frequency noise and false images caused by improper wave paths. To eliminate low-frequency noise and improve the quality of the migration image, we propose to separate the left-up-going, left-down-going, right-up-going and right-down-going wavefield components in the forward- and backward-propagated wavefields based on the Hilbert transform. By applying the reverse time migration of ground-penetrating radar data with full wavefield decomposition based on the Hilbert transform, we obtain the reverse time migration images of different wavefield components and combine correct imaging conditions to generate complete migration images. The proposed method is tested on the synthetic ground-penetrating radar data of a tilt-interface model and a complex model. The migration results show that the imaging condition of different wavefield components can highlight the desired structures. We further discuss the reasons for incomplete images by reverse time migration with partial wavefields. Compared with the conventional reverse time migration methods for ground-penetrating radar data, low-frequency noise can be eliminated in images generated by the reverse time migration method with full wavefield decomposition based on the Hilbert transform.     

Wavefield tomography based on local image correlations

The estimation of a velocity model from seismic data is a crucial step for obtaining a high‐quality image of the subsurface. Velocity estimation is usually formulated as an optimization problem where an objective function measures the mismatch between synthetic and recorded wavefields and its gradient is used to update the model. The objective function can be defined in the data‐space (as in full‐waveform inversion) or in the image space (as in migration velocity analysis). In general, the latter leads to smooth objective functions, which are monomodal in a wider basin about the global minimum compared to the objective functions defined in the data‐space. Nonetheless, migration velocity analysis requires construction of common‐image gathers at fixed spatial locations and subsampling of the image in order to assess the consistency between the trial velocity model and the observed data. We present an objective function that extracts the velocity error information directly in the image domain without analysing the information in common‐image gathers. In order to include the full complexity of the wavefield in the velocity estimation algorithm, we consider a two‐way (as opposed to one‐way) wave operator, we do not linearize the imaging operator with respect to the model parameters (as in linearized wave‐equation migration velocity analysis) and compute the gradient of the objective function using the adjoint‐state method. We illustrate our methodology with a few synthetic examples and test it on a real 2D marine streamer data set.     

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

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

DEPTH IMAGING OF OFFSET VERTICAL SEISMIC PROFILE DATA1

Depth migration consists of two different steps: wavefield extrapolation and imaging. The wave propagation is firmly founded on a mathematical frame-work, and is simulated by solving different types of wave equations, dependent on the physical model under investigation. In contrast, the imaging part of migration is usually based on ad hoc‘principles’, rather than on a physical model with an associated mathematical expression. The imaging is usually performed using the U/D concept of Claerbout (1971), which states that reflectors exist at points in the subsurface where the first arrival of the downgoing wave is time-coincident with the upgoing wave. Inversion can, as with migration, be divided into the two steps of wavefield extrapolation and imaging. In contrast to the imaging principle in migration, imaging in inversion follows from the mathematical formulation of the problem. The image with respect to the bulk modulus (or velocity) perturbations is proportional to the correlation between the time derivatives of a forward-propagated field and a backward-propagated residual field (Lailly 1984; Tarantola 1984). We assume a physical model in which the wave propagation is governed by the 2D acoustic wave equation. The wave equation is solved numerically using an efficient finite-difference scheme, making simulations in realistically sized models feasible. The two imaging concepts of migration and inversion are tested and compared in depth imaging from a synthetic offset vertical seismic profile section. In order to test the velocity sensitivity of the algorithms, two erroneous input velocity models are tested. We find that the algorithm founded on inverse theory is less sensitive to velocity errors than depth migration using the more ad hoc U/D imaging principle.     

基于解析时间波场外推与波场分解的逆时偏移方法研究

双程波方程逆时深度偏移是复杂介质高精度成像的有效技术,但其结果中通常包含成像方法引起的噪音和假象,一般的滤波方法会破坏成像剖面上的振幅,其中的假象也会给后续地质解释带来困扰.将波场进行方向分解然后实现入射波与反射波的相关成像能够有效地消除这类成像噪音,并提高逆时偏移成像质量.波传播方向的分解通常在频率波数域实现,它会占用大量的存储和计算资源,不便于在沿时间外推的逆时深度偏移中应用.本文提出解析时间波场外推方法,可以在时间外推的每个时间片上实现波传播方向的显式分解,逆时深度偏移中利用分解后的炮检波场进行对应的相关运算,实现成像噪音和成像信号的分离.在模型和实际数据上的测试表明,相比于常规互相关逆时偏移成像结果,本文方法能够有效地消除低频成像噪音和特殊地质构造导致的成像假象.     

退化的Fourier偏移算子及其在复杂断块成像中的应用

波动方程宽角抛物逼近得到的通常是非常系数的单程波传播算子,其系数是速度横向变化的函数,因此需要利用有限差分(FD)进行数值实施. 通过对Lippmann Schwinger单程波动积分方程的退化核逼近,本文研究了一类宽角退化算子的偏移成像. 这种退化偏移算子只用快速Fourier变换进行波场延拓,将常规的Fourier分裂步地震偏移方法(SSF)推广适应强速度横向变化介质和大角度传播波场. 退化的Fourier偏移算子通过在两个分裂步项之间作波数域线性插值来实现波场延拓,每延拓一层需要比常规的SSF地震偏移方法多一次快速Fourier变换(FFT). 通过SEG/EAGE盐丘模型和实际地震资料的应用表明,退化Fourier偏移算子能很好地对盐下的陡倾角断层和实际地震剖面上的复杂小断块和大断裂地质构造成像.     

Angle-domain imaging of relative impedance perturbation

Impedance is a physical parameter that plays an important role in seismic data processing and interpretation. A relative impedance perturbation (the ratio of the impedance perturbation and the impedance for the background models) imaging method in depth domain based on the reflection wave equation is proposed. Under the small perturbation assumption, primary wave and high-frequency approximation condition, a linear propagation equation of the primary reflection waves based on the relative impedance perturbation was first derived. On this basis, we further derived the imaging formula of the relative impedance perturbation using a linear inversion theory. Then, the source–receiver bidirectional illumination compensation was used to improve the image quality of the subsurface structures. The image result obtained by this method can be used to estimate the relative impedance perturbation. In the angle domain, the extracted near-angle-domain image gather with amplitude compensation can estimate the relative impedance perturbation, and the far-angle image gather provides the estimation of the relative velocity perturbation (the ratio of the velocity perturbation and the background velocity). Finally, several numerical tests demonstrate the effectiveness of the method.