地震偏移反演成像的迭代正则化方法研究

利用伴随算子L^*,直接的偏移方法通常导致一个低分辨率或模糊的地震成像.线性化偏移反演方法需求解一个最小二乘问题.但直接的最小二乘方法的数值不稳定,为目视解译带来困难.本文建立约束正则化数学模型,研究了地震偏移反演成像问题的迭代正则化求解方法.首先对最小二乘问题施加正则化约束,接着利用梯度迭代法求解反演成像问题,特别是提出了共轭梯度方法的混合实现技巧.为了表征该方法的可实际利用性,分别对一维,二维和三维地震模型进行了数值模拟.结果表明该正则偏移反演成像方法是有效的,对于实际的地震成像问题有着良好的应用前景.     

最小二乘偏移方法研究进展综述

最小二乘偏移是一种基于线性化反演理论的真振幅成像方法,其思路是在宏观背景模型的基础上估计出一个最优化的扰动部分对偏移结果进行迭代更新.该方法具有更高的成像精度,是实现地震成像理论由常规地下岩性的几何结构描述向真振幅成像的推进和发展,也是实现高精度储层反演的关键.本文阐述了最小二乘偏移的基本原理,指出了最小二乘偏移与常规偏移的本质区别;介绍了最小二乘偏移的发展历程及研究现状;分析了最小二乘偏移实现过程中的核心问题—Hessian矩阵和反演解的约束条件;探讨了最小二乘偏移存在的问题及今后的的发展趋势,为最小二乘偏移的进一步研究提供参考.     

基于平面波照明的偏移成像补偿

受地下复杂构造和地震数据采集系统的影响,地震波对地下目标的照明出现不均匀,在地震数据的偏移成像中出现成像阴影.根据地震数据最小二乘偏移/反演理论,和把地震波场照明结果作为最小二乘偏移/反演中的Hessian矩阵的近似对偏移成像进行补偿的原理,提出一种应用平面波照明结果对平面波偏移成像结果进行补偿以消除偏移成像阴影的方法.这种基于平面波照明的偏移成像补偿方法相对于局部角度域的照明偏移成像补偿方法具有计算效率上的优势.     

基于照明补偿的单程波最小二乘偏移

最小二乘偏移是一种基于反射地震数据与地下反射率间线性关系而建立起来的地震数据线性反演方法,相比常规偏移成像具有更好的保幅性能.本文提出了一种基于照明补偿的单程波最小二乘偏移方法,首先利用单程波方程的稳定Born近似广义屏波场传播算子构建反射地震数据与地下反射率间的线性算子,然后再应用线性最优化方法求解最小二乘偏移所对应的线性反问题.在迭代求解最优化问题的过程中,以地震波场的地下照明强度作为迭代反演的预条件算子加快迭代的收敛速度.单程波传播过程中考虑了速度分界面产生的透射效应,并用单极震源代替常规偏移中的偶极震源.把本文提出的方法应用于层状理论模型和Marmosi模型地震数据的数值试验中均取得了理想的结果.     

全波形反演在缝洞型储层速度建模中的应用

速度是地震偏移成像准确与否的关键所在.全波形反演综合利用地震波场运动学和动力学信息,能够得到相比传统速度建模方法更高频的成分.全波形反演的理论比较成熟,但实际应用成功的例子相对较少,特别是对于陆上地震资料.塔里木盆地地震地质条件复杂,为了实现缝洞型储层的准确成像,本文开展了针对目标靶区的全波形反演精细速度建场研究.采用一种时间域分层多尺度全波形反演流程:首先通过层析成像建立初始速度模型;其次利用折射波反演浅层速度模型;最后利用反射波反演中深层速度模型.偏移成像结果表明基于全波形反演的速度建模技术能有效改善火成岩下伏构造的成像精度,显示了全波形反演在常规陆上采集资料的应用潜力.     

基于归一化能量谱目标函数的全波形反演方法

基于L2范数的常规全波形反演目标函数是一个强非线性泛函,在反演过程中容易陷入局部极小值.本文提出归一化能量谱目标函数来缓解全波形反演过程中的强非线性问题,同时能够有效地缓解噪声和震源子波不准等因素的影响.能量谱目标函数是通过匹配观测数据与模拟数据随频率分布的能量信息来实现最小二乘反演的,其忽略了地震数据波形与相位变化的细节特征,这在反演的过程中能够有效缓解波形匹配错位等问题.数值测试结果表明,基于归一化能量谱目标函数在构建初始速度模型、抗噪性和缓解震源子波依赖等方面都优于归一化全波形反演目标函数.金属矿模型测试结果表明,即使地震数据缺失低频分量,基于归一化能量谱目标函数的全波形反演方法在像金属矿这样的强散射介质反演问题上同样具有一定的优势.     

基于海底电缆数据的声弹耦合介质多参数最小二乘逆时偏移（英文）

在海洋地震勘探中,海底电缆技术采集的多分量地震数据,涉及到多参数的反演与成像问题,本文针对海底电缆多分量地震数据提出了一种弹性波多参数最小二乘逆时偏移方法。该方法的波场延拓算子为混合方程,即在海水介质中采用声波方程进行波场计算,而海底固体介质的波场由纵横波形分离的矢量弹性波方程得到。在海底界面中采用声弹耦合控制方程将两种类型的方程结合起来。通过推导纵横波形分离的弹性波偏移算子、反偏移算子和梯度公式实现基于纵横波形分离的弹性波最小二乘逆时偏移方法。通过模型试算证明了该方法能够得到高质量的纵波速度和横波速度分量的成像剖面,相比与传统弹性波最小二乘逆时偏移方法,多参数耦合造成的成像串扰噪音得到了很好地压制。     

基于Seislet变换的稀疏约束多震源最小二乘逆时偏移

多震源地震采集技术允许一次性激发不同位置处的震源,得到来自多个震源的混合地震数据,该技术采集效率高,能有效降低采集成本.多震源地震数据成像效率高,但在偏移剖面中会引入串扰噪声,影响成像精度.最小二乘偏移常被用于压制多震源地震数据成像中的串扰噪声,但常规的最小二乘偏移并不能很好的消除串扰噪声对成像结果的影响,难以满足成像精度的要求.因此,为了保证反演的稳定性并改善反演结果,根据反射系数在Seislet域的稀疏性,本文引入了Seislet变换作为变换域稀疏约束的变换算子,实现了基于Seislet变换的稀疏约束多震源最小二乘逆时偏移,数值实验表明该方法能有效压制串扰噪声.     

基于高斯束的对角Hessian在保幅成像中的应用

稀疏的接收点采样、较窄的采集孔径以及有限的波场带宽等因素导致偏移成像的振幅往往是欠估计的.基于高斯束叠加积分表征的格林函数,可以显示地计算对角Hessian,并将其逆直接应用于高斯束偏移成像,从而实现模型空间的最小二乘偏移.通过常速度模型的数值试验验证了基于高斯束的对角Hessian显示算法的正确性和有效性,将此算法运用到EAGESEG盐丘模型数据和实际资料的最小二乘偏移可以产生振幅均衡的成像结果,尤其在弱照明和阴影区的振幅补偿效果明显.     

局部倾角约束最小二乘偏移方法研究

随着石油勘探难度的进一步加大,地震数据往往存在采样不规则、地震道缺失等现象,如果不对其进行处理,会对后续的地震成像产生影响,引入成像噪音.针对这一问题,一般是通过地震道插值或数据规则化对叠前数据进行处理,然后采用常规的偏移方法进行成像,本文则是将地震成像看作最小二乘反演问题,在共成像点道集引入平滑算子,在共偏移距/角度道集引入平面波构造算子(PWC)进行约束,通过预条件共轭梯度法使得反偏移后数据与输入数据之间的误差达到最小,最终得到信噪比更高、振幅属性更为可靠的成像结果.理论模型和实际资料处理表明,本文方法不仅可以有效压制数据不规则对成像产生的噪音,而且具有更高的成像精度.     

Reflection multi‐scale envelope inversion

Sufficient low‐frequency information is essential for full‐waveform inversion to get the global optimal solution. Multi‐scale envelope inversion was proposed using a new Fréchet derivative to invert the long‐wavelength component of the model by directly using the low‐frequency components contained in an envelope of seismic data. Although the new method can recover the main structure of the model, the inversion quality of the model bottom still needs to be improved. Reflection waveform inversion reduces the dependence of inversion on low‐frequency and long‐offset data by using travel‐time information in reflected waves. However, when the underground medium contains strong contrast or the initial model is far away from the true model, it is hard to get reliable reference reflectors for the generation of reflected waves. Here, we propose a combination inversion algorithm, i.e., reflection multi‐scale envelope inversion, to overcome the limitations of multi‐scale envelope inversion and reflection waveform inversion. First, wavefield decomposition was introduced into the multi‐scale envelope inversion to improve the inversion quality of the long‐wavelength components of the model. Then, after the initial model had been established to be accurate enough, migration and de‐migration were introduced to achieve multi‐scale reflection waveform inversion. The numerical results of the salt‐layer model and the SEG/EAGE salt model verified the validity of the proposed approach and its potential.     

Elastic full waveform inversion with probabilistic petrophysical clustering

Full waveform inversion aims to use all information provided by seismic data to deliver high-resolution models of subsurface parameters. However, multiparameter full waveform inversion suffers from an inherent trade-off between parameters and from ill-posedness due to the highly non-linear nature of full waveform inversion. Also, the models recovered using elastic full waveform inversion are subject to local minima if the initial models are far from the optimal solution. In addition, an objective function purely based on the misfit between recorded and modelled data may honour the seismic data, but disregard the geological context. Hence, the inverted models may be geologically inconsistent, and not represent feasible lithological units. We propose that all the aforementioned difficulties can be alleviated by explicitly incorporating petrophysical information into the inversion through a penalty function based on multiple probability density functions, where each probability density function represents a different lithology with distinct properties. We treat lithological units as clusters and use unsupervised K-means clustering to separate the petrophysical information into different units of distinct lithologies that are not easily distinguishable. Through several synthetic examples, we demonstrate that the proposed framework leads full waveform inversion to elastic models that are superior to models obtained either without incorporating petrophysical information, or with a probabilistic penalty function based on a single probability density function.     

Full waveform inversion and the truncated Newton method: quantitative imaging of complex subsurface structures

Full waveform inversion is a powerful tool for quantitative seismic imaging from wide‐azimuth seismic data. The method is based on the minimization of the misfit between observed and simulated data. This amounts to the solution of a large‐scale nonlinear minimization problem. The inverse Hessian operator plays a crucial role in this reconstruction process. Accounting accurately for the effect of this operator within the minimization scheme should correct for illumination deficits, restore the amplitude of the subsurface parameters, and help to remove artefacts generated by energetic multiple reflections. Conventional minimization methods (nonlinear conjugate gradient, quasi‐Newton methods) only roughly approximate the effect of this operator. In this study, we are interested in the truncated Newton minimization method. These methods are based on the computation of the model update through a matrix‐free conjugate gradient solution of the Newton linear system. We present a feasible implementation of this method for the full waveform inversion problem, based on a second‐order adjoint state formulation for the computation of Hessian‐vector products. We compare this method with conventional methods within the context of 2D acoustic frequency full waveform inversion for the reconstruction of P‐wave velocity models. Two test cases are investigated. The first is the synthetic BP 2004 model, representative of the Gulf of Mexico geology with high velocity contrasts associated with the presence of salt structures. The second is a 2D real data‐set from the Valhall oil field in North sea. Although, from a computational cost point of view, the truncated Newton method appears to be more expensive than conventional optimization algorithms, the results emphasize its increased robustness. A better reconstruction of the P‐wave velocity model is provided when energetic multiple reflections make it difficult to interpret the seismic data. A better trade‐off between regularization and resolution is obtained when noise contamination of the data requires one to regularize the solution of the inverse problem.     

Wave equation least square imaging using the local angular Hessian for amplitude correction

Local angular Hessian can be used to improve wave equation least square migration images. By decomposing the original Hessian operator into the local wavenumber domain or the local angle domain, the least square migration image is obtained as the solution of a linearized least‐squares inversion in the frequency and local angle domains. The local angular Hessian contains information about the acquisition geometry and the propagation effects based on the given velocity model. The inversion scheme based on the local angular Hessian avoids huge computation on the exact inverse Hessian matrix. To reduce the instability in the inversion, damping factors are introduced into the deconvolution filter in the local wavenumber domain and the local angle domain. The algorithms are tested using the SEG/EAGE salt2D model and the Sigsbee2A model. Results show improved image quality and amplitudes.     

基于截断牛顿法的VTI介质声波多参数全波形反演

不同类别参数间的相互耦合使多参数地震全波形反演的非线性程度显著增加,地震波速度与各向异性参数取值数量级的巨大差异也会使反演问题的性态变差.合理使用Hessian逆算子可以减弱这两类问题对反演的影响,提高多参数反演的精度,而截断牛顿法是一种可以比较准确地估计Hessian逆算子的优化方法.本文采用截断牛顿法在时间域进行了VTI介质的声波双参数同时反演的研究.不同模型的反演试验表明,在VTI介质声波双参数同时反演中,截断牛顿法比有限内存BFGS(Limited-memory Broyden-Fletcher-Goldfarb-Shanno,L-BFGS)法能更准确地估计Hessian逆算子,进而较好地平衡两类不同参数的同时更新,得到了比较精确的反演结果.     

Velocity model building from seismic reflection data by full‐waveform inversion

Full‐waveform inversion is re‐emerging as a powerful data‐fitting procedure for quantitative seismic imaging of the subsurface from wide‐azimuth seismic data. This method is suitable to build high‐resolution velocity models provided that the targeted area is sampled by both diving waves and reflected waves. However, the conventional formulation of full‐waveform inversion prevents the reconstruction of the small wavenumber components of the velocity model when the subsurface is sampled by reflected waves only. This typically occurs as the depth becomes significant with respect to the length of the receiver array. This study first aims to highlight the limits of the conventional form of full‐waveform inversion when applied to seismic reflection data, through a simple canonical example of seismic imaging and to propose a new inversion workflow that overcomes these limitations. The governing idea is to decompose the subsurface model as a background part, which we seek to update and a singular part that corresponds to some prior knowledge of the reflectivity. Forcing this scale uncoupling in the full‐waveform inversion formalism brings out the transmitted wavepaths that connect the sources and receivers to the reflectors in the sensitivity kernel of the full‐waveform inversion, which is otherwise dominated by the migration impulse responses formed by the correlation of the downgoing direct wavefields coming from the shot and receiver positions. This transmission regime makes full‐waveform inversion amenable to the update of the long‐to‐intermediate wavelengths of the background model from the wide scattering‐angle information. However, we show that this prior knowledge of the reflectivity does not prevent the use of a suitable misfit measurement based on cross‐correlation, to avoid cycle‐skipping issues as well as a suitable inversion domain as the pseudo‐depth domain that allows us to preserve the invariant property of the zero‐offset time. This latter feature is useful to avoid updating the reflectivity information at each non‐linear iteration of the full‐waveform inversion, hence considerably reducing the computational cost of the entire workflow. Prior information of the reflectivity in the full‐waveform inversion formalism, a robust misfit function that prevents cycle‐skipping issues and a suitable inversion domain that preserves the seismic invariant are the three key ingredients that should ensure well‐posedness and computational efficiency of full‐waveform inversion algorithms for seismic reflection data.     

Full waveform inversion of SH‐ and Love‐wave data in near‐surface prospecting

We develop a two‐dimensional full waveform inversion approach for the simultaneous determination of S‐wave velocity and density models from SH ‐ and Love‐wave data. We illustrate the advantages of the SH/Love full waveform inversion with a simple synthetic example and demonstrate the method's applicability to a near‐surface dataset, recorded in the village ?achtice in Northwestern Slovakia. Goal of the survey was to map remains of historical building foundations in a highly heterogeneous subsurface. The seismic survey comprises two parallel SH‐profiles with maximum offsets of 24 m and covers a frequency range from 5 Hz to 80 Hz with high signal‐to‐noise ratio well suited for full waveform inversion. Using the Wiechert–Herglotz method, we determined a one‐dimensional gradient velocity model as a starting model for full waveform inversion. The two‐dimensional waveform inversion approach uses the global correlation norm as objective function in combination with a sequential inversion of low‐pass filtered field data. This mitigates the non‐linearity of the multi‐parameter inverse problem. Test computations show that the influence of visco‐elastic effects on the waveform inversion result is rather small. Further tests using a mono‐parameter shear modulus inversion reveal that the inversion of the density model has no significant impact on the final data fit. The final full waveform inversion S‐wave velocity and density models show a prominent low‐velocity weathering layer. Below this layer, the subsurface is highly heterogeneous. Minimum anomaly sizes correspond to approximately half of the dominant Love‐wavelength. The results demonstrate the ability of two‐dimensional SH waveform inversion to image shallow small‐scale soil structure. However, they do not show any evidence of foundation walls.     

Full waveform inversion using oriented time‐domain imaging method for vertical transverse isotropic media

Full waveform inversion for reflection events is limited by its linearised update requirements given by a process equivalent to migration. Unless the background velocity model is reasonably accurate, the resulting gradient can have an inaccurate update direction leading the inversion to converge what we refer to as local minima of the objective function. In our approach, we consider mild lateral variation in the model and, thus, use a gradient given by the oriented time‐domain imaging method. Specifically, we apply the oriented time‐domain imaging on the data residual to obtain the geometrical features of the velocity perturbation. After updating the model in the time domain, we convert the perturbation from the time domain to depth using the average velocity. Considering density is constant, we can expand the conventional 1D impedance inversion method to two‐dimensional or three‐dimensional velocity inversion within the process of full waveform inversion. This method is not only capable of inverting for velocity, but it is also capable of retrieving anisotropic parameters relying on linearised representations of the reflection response. To eliminate the crosstalk artifacts between different parameters, we utilise what we consider being an optimal parametrisation for this step. To do so, we extend the prestack time‐domain migration image in incident angle dimension to incorporate angular dependence needed by the multiparameter inversion. For simple models, this approach provides an efficient and stable way to do full waveform inversion or modified seismic inversion and makes the anisotropic inversion more practicable. The proposed method still needs kinematically accurate initial models since it only recovers the high‐wavenumber part as conventional full waveform inversion method does. Results on synthetic data of isotropic and anisotropic cases illustrate the benefits and limitations of this method.     

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

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

Assumptions and goals for least squares migration

Least squares migration uses the assumption that, if we have an operator that can create data from a reflectivity function, the optimal image will predict the actual recorded data with minimum square error. For this assumption to be true, it is also required that: (a) the prediction operator must be error-free, (b) model elements not seen by the operator should be constrained by other means and (c) data weakly predicted by the operator should make limited contribution to the solution. Under these conditions, least squares migration has the advantage over simple migration of being able to remove interference between different model components. Least squares migration does that by de-convolving or inverting the so-called Hessian operator. The Hessian is the cascade of forward modelling and migration; for each image point, it computes the effects of interference from other image points (point-spread function) given the actual recording geometry and the subsurface velocity model. Because the Hessian contains illumination information (along its diagonal), and information about the model cross-correlation produced by non-orthogonality of basis functions, its inversion produces illumination compensation and increases resolution. In addition, sampling deficiencies in the recording geometry map to the Hessian (both diagonal and non-diagonal elements), so least squares migration has the potential to remove sampling artefacts as well. These (illumination compensation, resolution and mitigating recording deficiencies) are the three main goals of least squares migration, although the first one can be achieved by cheaper techniques. To invert the Hessian, least squares migration relies on the residual errors during iterations. Iterative algorithms, like conjugate gradient and others, use the residuals to calculate the direction and amplitudes (gradient and step size) of the necessary corrections to the reflectivity function or model. Failure of conditions (a), (b) or (c) leads the inversion to calculate incorrect model updates, which translate to noise in the final image. In this paper, we will discuss these conditions for Kirchhoff migration and reverse time migration.