Preconditioned non‐linear conjugate gradient method for frequency domain full‐waveform seismic inversion

We present preconditioned non‐linear conjugate gradient algorithms as alternatives to the Gauss‐Newton method for frequency domain full‐waveform seismic inversion. We designed two preconditioning operators. For the first preconditioner, we introduce the inverse of an approximate sparse Hessian matrix. The approximate Hessian matrix, which is highly sparse, is constructed by judiciously truncating the Gauss‐Newton Hessian matrix based on examining the auto‐correlation and cross‐correlation of the Jacobian matrix. As the second preconditioner, we employ the approximation of the inverse of the Gauss‐Newton Hessian matrix. This preconditioner is constructed by terminating the iteration process of the conjugate gradient least‐squares method, which is used for inverting the Hessian matrix before it converges. In our preconditioned non‐linear conjugate gradient algorithms, the step‐length along the search direction, which is a crucial factor for the convergence, is carefully chosen to maximize the reduction of the cost function after each iteration. The numerical simulation results show that by including a very limited number of non‐zero elements in the approximate Hessian, the first preconditioned non‐linear conjugate gradient algorithm is able to yield comparable inversion results to the Gauss‐Newton method while maintaining the efficiency of the un‐preconditioned non‐linear conjugate gradient method. The only extra cost is the computation of the inverse of the approximate sparse Hessian matrix, which is less expensive than the computation of a forward simulation of one source at one frequency of operation. The second preconditioned non‐linear conjugate gradient algorithm also significantly saves the computational expense in comparison with the Gauss‐Newton method while maintaining the Gauss‐Newton reconstruction quality. However, this second preconditioned non‐linear conjugate gradient algorithm is more expensive than the first one.     

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

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

Frequency‐domain acoustic‐elastic coupled waveform inversion using the Gauss‐Newton conjugate gradient method

We developed a frequency‐domain acoustic‐elastic coupled waveform inversion based on the Gauss‐Newton conjugate gradient method. Despite the use of a high‐performance computer system and a state‐of‐the‐art parallel computation algorithm, it remained computationally prohibitive to calculate the approximate Hessian explicitly for a large‐scale inverse problem. Therefore, we adopted the conjugate gradient least‐squares algorithm, which is frequently used for geophysical inverse problems, to implement the Gauss‐Newton method so that the approximate Hessian is calculated implicitly. Thus, there was no need to store the Hessian matrix. By simultaneously back‐propagating multi‐components consisting of the pressure and displacements, we could efficiently extract information on the subsurface structures. To verify our algorithm, we applied it to synthetic data sets generated from the Marmousi‐2 model and the modified SEG/EAGE salt model. We also extended our algorithm to the ocean‐bottom cable environment and verified it using ocean‐bottom cable data generated from the Marmousi‐2 model. With the assumption of a hard seafloor, we recovered both the P‐wave velocity of complicated subsurface structures as well as the S‐wave velocity. Although the inversion of the S‐wave velocity is not feasible for the high Poisson's ratios used to simulate a soft seafloor, several strategies exist to treat this problem. Our example using multi‐component data showed some promise in mitigating the soft seafloor effect. However, this issue still remains open.     

地震反演成像中的Hessian算子研究

总结了牛顿类地震反演方法中Hessian算子的作用,对其在地震反演成像中的数学物理含义进行了分析.Hessian算子是误差泛函对模型参数的二阶导数,反映了误差泛函对模型变化的二次型特征.分析声波方程下的Hessian算子的格林函数表达形式,发现其表达了整个观测系统和子波频带等因素对地震数据空间到模型空间投影过程的影响.提出了两种分别适用于最小二乘偏移和全波形反演的Hessian算子简化格式.平面波Hessian算子应用于最小二乘偏移能够得到相对保真的成像结果,改善了地震偏移成像的精度.地下偏移距Hessian算子应用于全波形反演能够加快反演迭代的计算效率.最后,对Hessian算子在地震反演成像中的价值进行了讨论和评价.     

Double parameterized regularization inversion method for migration velocity analysis in transversely isotropic media with a vertical symmetry axis

Simultaneous estimation of velocity gradients and anisotropic parameters from seismic reflection data is one of the main challenges in transversely isotropic media with a vertical symmetry axis migration velocity analysis. In migration velocity analysis, we usually construct the objective function using the l₂ norm along with a linear conjugate gradient scheme to solve the inversion problem. Nevertheless, for seismic data this inversion scheme is not stable and may not converge in finite time. In order to ensure the uniform convergence of parameter inversion and improve the efficiency of migration velocity analysis, this paper develops a double parameterized regularization model and gives the corresponding algorithms. The model is based on the combination of the l₂ norm and the non‐smooth l₁ norm. For solving such an inversion problem, the quasi‐Newton method is utilized to make the iterative process stable, which can ensure the positive definiteness of the Hessian matrix. Numerical simulation indicates that this method allows fast convergence to the true model and simultaneously generates inversion results with a higher accuracy. Therefore, our proposed method is very promising for practical migration velocity analysis in anisotropic media.     

Migration velocity analysis and waveform inversion

Least‐squares inversion of seismic reflection waveform data can reconstruct remarkably detailed models of subsurface structure and take into account essentially any physics of seismic wave propagation that can be modelled. However, the waveform inversion objective has many spurious local minima, hence convergence of descent methods (mandatory because of problem size) to useful Earth models requires accurate initial estimates of long‐scale velocity structure. Migration velocity analysis, on the other hand, is capable of correcting substantially erroneous initial estimates of velocity at long scales. Migration velocity analysis is based on prestack depth migration, which is in turn based on linearized acoustic modelling (Born or single‐scattering approximation). Two major variants of prestack depth migration, using binning of surface data and Claerbout's survey‐sinking concept respectively, are in widespread use. Each type of prestack migration produces an image volume depending on redundant parameters and supplies a condition on the image volume, which expresses consistency between data and velocity model and is hence a basis for velocity analysis. The survey‐sinking (depth‐oriented) approach to prestack migration is less subject to kinematic artefacts than is the binning‐based (surface‐oriented) approach. Because kinematic artefacts strongly violate the consistency or semblance conditions, this observation suggests that velocity analysis based on depth‐oriented prestack migration may be more appropriate in kinematically complex areas. Appropriate choice of objective (differential semblance) turns either form of migration velocity analysis into an optimization problem, for which Newton‐like methods exhibit little tendency to stagnate at nonglobal minima. The extended modelling concept links migration velocity analysis to the apparently unrelated waveform inversion approach to estimation of Earth structure: from this point of view, migration velocity analysis is a solution method for the linearized waveform inversion problem. Extended modelling also provides a basis for a nonlinear generalization of migration velocity analysis. Preliminary numerical evidence suggests a new approach to nonlinear waveform inversion, which may combine the global convergence of velocity analysis with the physical fidelity of model‐based data fitting.     

基于Born敏感核函数的速度、密度双参数全波形反演

速度、密度之间的相互耦合使得密度在多参数全波形反演中较难获得.本文将截断高斯-牛顿法用于声介质速度、密度双参数全波形反演,通过考虑近似Hessian矩阵中反映速度、密度相互作用的非主对角块元素,有效解决了多参数全波形反演中速度、密度之间的耦合问题,在不采用反演策略的情况下,仍能够获得精度较高的速度、密度反演结果.常规的截断牛顿类全波形反演通常利用一阶伴随状态法求取目标函数对模型参数的梯度,利用二阶伴随状态法或有限差分法求解Hessian-向量乘,在每一步内循环迭代过程中需要额外求解两次正演问题,计算量较大.本文基于Born近似,将梯度计算中的核函数-向量乘表示为具有明确物理意义的向量-标量乘的累加运算,同时将Hessian-向量乘转化为两次核函数-向量乘,无需额外求解正演问题,有效降低了计算量.数值实验证明了本文提出的方法的有效性.     

Non-monotone spectral projected gradient method applied to full waveform inversion

The seismic inversion problem is a highly non‐linear problem that can be reduced to the minimization of the least‐squares criterion between the observed and the modelled data. It has been solved using different classical optimization strategies that require a monotone descent of the objective function. We propose solving the full‐waveform inversion problem using the non‐monotone spectral projected gradient method: a low‐cost and low‐storage optimization technique that maintains the velocity values in a feasible convex region by frequently projecting them on this convex set. The new methodology uses the gradient direction with a particular spectral step length that allows the objective function to increase at some iterations, guarantees convergence to a stationary point starting from any initial iterate, and greatly speeds up the convergence of gradient methods. We combine the new optimization scheme as a solver of the full‐waveform inversion with a multiscale approach and apply it to a modified version of the Marmousi data set. The results of this application show that the proposed method performs better than the classical gradient method by reducing the number of function evaluations and the residual values.     

Linearized seismic waveform inversion using a multiple re-weighted least-squares method with QR preconditioning

Linearized inversion methods such as Gauss‐Newton and multiple re‐weighted least‐squares are iterative processes in which an update in the current model is computed as a function of data misfit and the gradient of data with respect to model parameters. The main advantage of those methods is their ability to refine the model parameters although they have a high computational cost for seismic inversion. In the Gauss‐Newton method a system of equations, corresponding to the sensitivity matrix, is solved in the least‐squares sense at each iteration, while in the multiple re‐weighted least‐squares method many systems are solved using the same sensitivity matrix. The sensitivity matrix arising from these methods is usually not sparse, thus limiting the use of standard preconditioners in the solution of the linearized systems. For reduction of the computational cost of the linearized inversion methods, we propose the use of preconditioners based on a partial orthogonalization of the columns of the sensitivity matrix. The new approach collapses a band of co‐diagonals of the normal equations matrix into the main diagonal, being equivalent to computing the least‐squares solution starting from a partial solution of the linear system. The preconditioning is driven by a bandwidth L which can be interpreted as the distance for which the correlation between model parameters is relevant. To illustrate the benefit of the proposed approach to the reduction of the computational cost of the inversion we apply the multiple re‐weighted least‐squares method to the 2D acoustic seismic waveform inversion problem. We verify the reduction in the number of iterations in the conjugate'gradient algorithm as the bandwidth of the preconditioners increases. This effect reduces the total computational cost of inversion as well.     

地震反射波反演二阶优化方法及其应用

常规地震数据大都缺乏低频成分与长炮检距信号,经典的全波形反演方法不易获取中、深层弹性参数模型的长波长分量,波动方程反射波形反演作为替代方法近来受到极大关注.然而,现有的反射波形反演方法几乎都采用梯度类的一阶优化算法,收敛性和精度都有待提高.本文在二阶优化理论框架下,推导弹性参数背景与扰动模型的反射波敏感核、泛函梯度以及海森算子,揭示海森矩阵对泛函梯度的去模糊化作用和改善反演的工作机制.推覆体模型合成数据实验表明,相比于常用的共轭梯度法,利用近似海森矩阵的高斯-牛顿法明显提升了反射波形反演的收敛性与宽谱建模能力.在东海实例中,本文方法超越常用的反射走时层析技术,通过改善中、深层偏移速度建模,支撑逆时偏移高分辨率刻画长江坳陷内部复杂的断裂系统,改善了深部基底的成像质量.

     

Image‐domain wavefield tomography with extended common‐image‐point gathers

Waveform inversion is a velocity‐model‐building technique based on full waveforms as the input and seismic wavefields as the information carrier. Conventional waveform inversion is implemented in the data domain. However, similar techniques referred to as image‐domain wavefield tomography can be formulated in the image domain and use a seismic image as the input and seismic wavefields as the information carrier. The objective function for the image‐domain approach is designed to optimize the coherency of reflections in extended common‐image gathers. The function applies a penalty operator to the gathers, thus highlighting image inaccuracies arising from the velocity model error. Minimizing the objective function optimizes the model and improves the image quality. The gradient of the objective function is computed using the adjoint state method in a way similar to that in the analogous data‐domain implementation. We propose an image‐domain velocity‐model building method using extended common‐image‐point space‐ and time‐lag gathers constructed sparsely at reflections in the image. The gathers are effective in reconstructing the velocity model in complex geologic environments and can be used as an economical replacement for conventional common‐image gathers in wave‐equation tomography. A test on the Marmousi model illustrates successful updating of the velocity model using common‐image‐point gathers and resulting improved image quality.     

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.     

频率多尺度全波形速度反演

以二维声波方程为模型,在时间域深入研究了全波形速度反演.全波形反演要解一个非线性的最小二乘问题,是一个极小化模拟数据与已知数据之间残量的过程.针对全波形反演易陷入局部极值的困难,本文提出了基于不同尺度的频率数据的"逐级反演"策略,即先基于低频尺度的波场信息进行反演,得出一个合理的初始模型,然后再利用其他不同尺度频率的波场进行反演,并且用前一尺度的迭代反演结果作为下一尺度反演的初始模型,这样逐级进行反演.文中详细阐述和推导了理论方法及公式,包括有限差分正演模拟、速度模型修正、梯度计算和算法描述,并以Marmousi复杂构造模型为例,进行了MPI并行全波形反演数值计算,得到了较好的反演结果,验证了方法的有效性和稳健性.     

Two‐dimensional waveform inversion of multi‐component data in acoustic‐elastic coupled media

In order to account for the effects of elastic wave propagation in marine seismic data, we develop a waveform inversion algorithm for acoustic‐elastic media based on a frequency‐domain finite‐element modelling technique. In our algorithm we minimize residuals using the conjugate gradient method, which back‐propagates the errors using reverse time migration without directly computing the partial derivative wavefields. Unlike a purely acoustic or purely elastic inversion algorithm, the Green's function matrix for our acoustic‐elastic algorithm is asymmetric. We are nonetheless able to achieve computational efficiency using modern numerical methods. Numerical examples show that our coupled inversion algorithm produces better velocity models than a purely acoustic inversion algorithm in a wide variety of cases, including both single‐ and multi‐component data and low‐cut filtered data. We also show that our algorithm performs at least equally well on real field data gathered in the Korean continental shelf.     

Non‐linear prestack seismic inversion with global optimization using an edge‐preserving smoothing filter

Estimating elastic parameters from prestack seismic data remains a subject of interest for the exploration and development of hydrocarbon reservoirs. In geophysical inverse problems, data and models are in general non‐linearly related. Linearized inversion methods often have the disadvantage of strong dependence on the initial model. When the initial model is far from the global minimum, inversion iteration is likely to converge to the local minimum. This problem can be avoided by using global optimization methods. In this paper, we implemented and tested a prestack seismic inversion scheme based on a quantum‐behaved particle swarm optimization (QPSO) algorithm aided by an edge‐preserving smoothing ( EPS) operator. We applied the algorithm to estimate elastic parameters from prestack seismic data. Its performance on both synthetic data and real seismic data indicates that QPSO optimization with the EPS operator yields an accurate solution.     

Exploring some issues in acoustic full waveform inversion

The least‐squares error measures the difference between observed and modelled seismic data. Because it suffers from local minima, a good initial velocity model is required to avoid convergence to the wrong model when using a gradient‐based minimization method. If a data set mainly contains reflection events, it is difficult to update the velocity model with the least‐squares error because the minimization method easily ends up in the nearest local minimum without ever reaching the global minimum. Several authors observed that the model could be updated by diving waves, requiring a wide‐angle or large‐offset data set. This full waveform tomography is limited to a maximum depth. Here, we use a linear velocity model to obtain estimates for the maximum depth. In addition, we investigate how frequencies should be selected if the seismic data are modelled in the frequency domain. In the presence of noise, the condition to avoid local minima requires more frequencies than needed for sufficient spectral coverage. We also considered acoustic inversion of a synthetic marine data set created by an elastic time‐domain finite‐difference code. This allowed us to validate the estimates made for the linear velocity model. The acoustic approximation leads to a number of problems when using long‐offset data. Nevertheless, we obtained reasonable results. The use of a variable density in the acoustic inversion helped to match the data at the expense of accuracy in the inversion result for the density.     

基于GPU并行的时间域全波形优化共轭梯度法快速GPR双参数反演

探地雷达（GPR）时间域全波形反演计算量巨大,内存要求高,在微机上计算难度大.本文中作者基于GPU并行加速的维度提升反演策略,采用优化的共轭梯度法,避免了Hessian矩阵的计算,在普通微机上实现了时间域全波形二维GPR双参数（介电常数和电导率）快速反演.论文首先推导了二维TM波的时域有限差分法（FDTD）的交错网格离散差分格式及波场更新策略.然后,基于Lagrange乘数法,将约束问题转化为无约束最小问题,构建了共轭梯度法反演目标函数,采用Fletcher-Reeves公式与非精确线搜索Wolfe准则,确保了梯度方向修正因子及迭代步长选取的合理性.而GPU并行计算及维度提升反演策略的应用,数倍地提升了反演速度.最后,开展了3个模型的合成数据的反演实验,分别从观测方式、梯度优化及天线频率等方面,分析了这些因素对雷达全波形反演的影响,说明双参数的反演较单一的介电常数反演,能提供更丰富的信息约束,有效提高模型重建的精度.

     

Research Note: Full‐waveform inversion of the unwrapped phase of a model

Reflections in seismic data induce serious non‐linearity in the objective function of full‐ waveform inversion. Thus, without a good initial velocity model that can produce reflections within a half cycle of the frequency used in the inversion, convergence to a solution becomes difficult. As a result, we tend to invert for refracted events and damp reflections in data. Reflection induced non‐linearity stems from cycle skipping between the imprint of the true model in observed data and the predicted model in synthesized data. Inverting for the phase of the model allows us to address this problem by avoiding the source of non‐linearity, the phase wrapping phenomena. Most of the information related to the location (or depths) of interfaces is embedded in the phase component of a model, mainly influenced by the background model, while the velocity‐contrast information (responsible for the reflection energy) is mainly embedded in the amplitude component. In combination with unwrapping the phase of data, which mitigates the non‐linearity introduced by the source function, I develop a framework to invert for the unwrapped phase of a model, represented by the instantaneous depth, using the unwrapped phase of the data. The resulting gradient function provides a mechanism to non‐linearly update the velocity model by applying mainly phase shifts to the model. In using the instantaneous depth as a model parameter, we keep track of the model properties unfazed by the wrapping phenomena.     

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.     

Migration velocity analysis using pre‐stack wave fields

Using both image and data domains to perform velocity inversion can help us resolve the long and short wavelength components of the velocity model, usually in that order. This translates to integrating migration velocity analysis into full waveform inversion. The migration velocity analysis part of the inversion often requires computing extended images, which is expensive when using conventional methods. As a result, we use pre‐stack wavefield (the double‐square‐root formulation) extrapolation, which includes the extended information (subsurface offsets) naturally, to make the process far more efficient and stable. The combination of the forward and adjoint pre‐stack wavefields provides us with update options that can be easily conditioned to improve convergence. We specifically use a modified differential semblance operator to split the extended image into a residual part for classic differential semblance operator updates and the image (Born) modelling part, which provides reflections for higher resolution information. In our implementation, we invert for the velocity and the image simultaneously through a dual objective function. Applications to synthetic examples demonstrate the features of the approach.