Seismic inversion with generalized Radon transform based on local second-order approximation of scattered field in acoustic media

Sound velocity inversion problem based on scattering theory is formulated in terms of a nonlinear integral equation associated with scattered field. Because of its nonlinearity, in practice, linearization algorisms (Born/single scattering approximation) are widely used to obtain an approximate inversion solution. However, the linearized strategy is not congruent with seismic wave propagation mechanics in strong perturbation (heterogeneous) medium. In order to partially dispense with the weak perturbation assumption of the Born approximation, we present a new approach from the following two steps: firstly, to handle the forward scattering by taking into account the second-order Born approximation, which is related to generalized Radon transform (GRT) about quadratic scattering potential; then to derive a nonlinear quadratic inversion formula by resorting to inverse GRT. In our formulation, there is a significant quadratic term regarding scattering potential, and it can provide an amplitude correction for inversion results beyond standard linear inversion. The numerical experiments demonstrate that the linear single scattering inversion is only good in amplitude for relative velocity perturbation ( \( \delta_{c}/c_{0} \) ) of background media up to 10 %, and its inversion errors are unacceptable for the perturbation beyond 10 %. In contrast, the quadratic inversion can give more accurate amplitude-preserved recovery for the perturbation up to 40 %. Our inversion scheme is able to manage double scattering effects by estimating a transmission factor from an integral over a small area, and therefore, only a small portion of computational time is added to the original linear migration/inversion process.     

地球构造反演问题的新途径

本文讨论地球内部构造反演问题的某些新途径.其内容如下:地球构造反问题与固有值反问题;反散射问题中介质间断性的成象与因果广义Radon变换,包含地球构造反问题的新提法,反散射问题的线性化,古典Radon变换与广义Radon变换,线性化反问题的渐近解,渐近解和偏移格式.     

波动方程CT技术在地震数据处理中的应用

本文从畸变的Born近似的微扰技术出发,给出了利用广义Radon变换和Fourier积分算子的理论反演介质间断性的原理.将声学的广义Radon变换与经典Radon变换进行类比,近似地导出了声学广义Radon变换的反演公式. 本文对于反射地震学的情况,提出了一种拟线性化方法,考虑了成象点的一次散射场,从某种程度上减少了Born近似对弱散射的苛求. 利用同一模型的理论记录和物理实验记录的反演计算结果对提出的方法进行了验证,并讨论了进一步提高成象精度的方法.     

各向异性Radon变换及其在多次波压制中的应用

即使采用分辨率很高的双曲Radon变换,对速度各向异性发育介质及长偏移距情况下的地震数据,其Radon域内能量仍不收敛.为了克服此难题,我们在Radon变换的积分路径中考虑了非双曲走时的影响,通过引入非双曲时差公式中的各向异性非椭圆率η参数,可以准确描述出长偏移距条件下来自同一层位的时距曲线,并推导了由偏移距、慢度、非椭圆率三参数控制的积分曲线正反变换公式,我们称之为各向异性Radon变换.离散化求解时,各向异性Radon变换是时变的,频率域快速算法已不适用,本文采用了最优相似系数加权Gauss-Seidel迭代算法,保持其计算精度的同时也有较高的计算效率.将此方法应用在模型数据以及实际长偏移距海上地震数据的多次波压制处理中,收到了较好的处理效果.     

波动方程反问题与广义RADON变换反演的一种关系

本文把波动方程反问题与广义Radon 变换的反演相联系。在假定弱散射条件下,把波动方程反问题转化成广义Radon 变换的反演问题,即如何从一系列关于目标函数在某类子流形上的积分值,去求出目标函数.这种转化提供了一种研究波动方程反问题的途径。     

Inversion-based deblending using migration operators

In this paper, we compare the denoising- and inversion-based deblending methods using Stolt migration operators. We use Stolt operator as a kernel to efficiently compute apex-shifted hyperbolic Radon transform. Sparsity promoting transforms, such as Radon transform, can focus seismic data into a sparse model to separate signals, remove noise or interpolate missing traces. Therefore, Radon transforms are a suitable tool for either the denoising- or the inversion-based deblending methods. The denoising-based deblending treats blending interferences as random noise by sorting the data into new gathers, such as common receiver gather. In these gathers, blending interferences exhibit random structures due to the randomization of the source firing times. Alternatively, the inversion-based deblending treats blending interferences as a signal, and the transform models this signal by incorporating the blending operator to formulate an inversion problem. We compare both methods using a robust inversion algorithm with sparse regularization. Results of synthetic and field data examples show that the inversion-based deblending can produce more accurate signal separation for highly blended data.     

PARAMETRIZATION OF GRT INVERSION FOR ACOUSTIC AND P–P SCATTERING1

The generalized Radon transform (GRT) inversion contains an explicit relationship between seismic amplitude variations, the reflection angle and the physical parameters which can be used to describe the earth efficiently for inversion purposes. Using this relationship, we have derived parametrizations for acoustic and P–P scattering so that the variations in seismic amplitude with reflection angle for each parameter are sufficiently independent. These parametrizations show that small offset and large offset amplitudes are related to different physical parameters. In the case of acoustic scattering, the small-offset amplitudes are related to impedance variations while large-offset amplitudes are related to velocity variations. A similar result has been established for P–P scattering. The Born approximation (which is used to derive the GRT inversion) does not correctly predict the amplitude due to velocity variations at large offsets, and thus the inversion of velocity is not as satisfactory as the inversion of impedance.     

2.5D modelling, inversion and angle migration in anisotropic elastic media

2.5D modelling approximates 3D wave propagation in the dip‐direction of a 2D geological model. Attention is restricted to raypaths for waves propagating in a plane. In this way, fast inversion or migration can be performed. For velocity analysis, this reduction of the problem is particularly useful. We review 2.5D modelling for Born volume scattering and Born–Helmholtz surface scattering. The amplitudes are corrected for 3D wave propagation, taking into account both in‐plane and out‐of‐plane geometrical spreading. We also derive some new inversion/migration results. An AVA‐compensated migration routine is presented that is simplified compared with earlier results. This formula can be used to create common‐image gathers for use in velocity analysis by studying the residual moveout. We also give a migration formula for the energy‐flux‐normalized plane‐wave reflection coefficient that models large contrast in the medium parameters not treated by the Born and the Born–Helmholtz equation results. All results are derived using the generalized Radon transform (GRT) directly in the natural coordinate system characterized by scattering angle and migration dip. Consequently, no Jacobians are needed in their calculation. Inversion and migration in an orthorhombic medium or a transversely isotropic (TI) medium with tilted symmetry axis are the lowest symmetries for practical purposes (symmetry axis is in the plane). We give an analysis, using derived methods, of the parameters for these two types of media used in velocity analysis, inversion and migration. The kinematics of the two media involve the same parameters, hence there is no distinction when carrying out velocity analysis. The in‐plane scattering coefficient, used in the inversion and migration, also depends on the same parameters for both media. The out‐of‐plane geometrical spreading, necessary for amplitude‐preserving computations, for the TI medium is dependent on the same parameters that govern in‐plane kinematics. For orthorhombic media, information on additional parameters is required that is not needed for in‐plane kinematics and the scattering coefficients. Resolution analysis of the scattering coefficient suggests that direct inversion by GRT yields unreliable parameter estimates. A more practical approach to inversion is amplitude‐preserving migration followed by AVA analysis. SYMBOLS AND NOTATION A list of symbols and notation is given in Appendix D .     

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.     

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

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

Time-domain hyperbolic Radon transform for separation of P-P and P-SV wavefields

A time-domain hyperbolic Radon transform based method for separating multicomponent seismic data into P-P and P-SV wavefields is presented. This wavefield separation method isolates P-P and P-SV wavefields in the Radon panel due to their differences in slowness, and an inverse transform of only part of the data leads to separated wavefields. A problem of hyperbolic Radon transform is that it works in the time domain entailing the inversion of large operators which is prohibitively time-consuming. By applying the conjugate gradient algorithm during the inversion of hyperbolic Radon transform, the computational cost can be kept reasonably low for practical application. Synthetic data examples prove that P-P and P-SV wavefield separation by hyperbolic Radon transform produces more accurate separated wavefields compared with separation by high-resolution parabolic Radon transform, and the feasibility of the proposed separation scheme is also verified by a real field data example.     

各向异性介质Low-rank有限差分法纯qP波叠前平面波最小二乘逆时偏移

拟声波最小二乘逆时偏移是一种极具潜力的地震波成像工具,但该方法遭受各向异性拟声波近似的限制,TTI介质正演模拟不稳定、反偏移记录中遭受伪横波二次扰动及数值频散假象,另外拟声波最小二乘逆时偏移还面临计算效率低、收敛速度慢、对速度等模型参数依赖性高等问题.为了克服各向异性拟声波最小二乘逆时偏移的缺陷,在反演框架下,本文借助Low-rank有限差分算法首次提出并实现了TTI介质纯qP波线性正演模拟及纯qP波最小二乘逆时偏移;为了进一步提升反演成像效率,同时改善反演成像方法对模型参数误差的依赖性及对地震数据噪声的适应性,通过引入叠前平面波优化策略,发展了TTI介质纯qP波叠前平面波最小二乘逆时偏移成像方法.在编程实现方法的基础上,通过开展模型成像测试,展示了本方法的优势和潜力：一方面加快了反演成像效率,另一方面也提升了方法的抗噪性,同时还降低了方法对模型参数的依赖性.     

井孔高分辨率Radon变换算子研究

本文针对井间和3D VSP波场的线性特征,研究井孔地震波场线性高分辨率Radon变换算子,用于井孔地震波场分析与纵横波分离.在Radon变换原理分析基础上,采用基于柯西分布的高分辨率线性Radon变换对井孔数据进行Radon变换,其间通过对离散倾角叠加算子求取的研究,及对影响Radon能量收敛的重要参数阻尼因子算法的改进,使数据在Radon域以能量团的形式呈现,得到很好的收敛效果,基本解决了Radon域数据的一定程度的拖尾现象,消除了各能量团之间的平滑效应,采用柯西分布来规则化数据,提高了Radon域的分辨率,Radon域能量也收敛到一个点上,有利于上下行波或纵横波波场分离.最后通过反演结果和模型试算验证了该方法的可行性和稳定性.     

Azimuthal Seismic Amplitude Variation with Offset and Azimuth Inversion in Weakly Anisotropic Media with Orthorhombic Symmetry

Seismic amplitude variation with offset and azimuth (AVOaz) inversion is well known as a popular and pragmatic tool utilized to estimate fracture parameters. A single set of vertical fractures aligned along a preferred horizontal direction embedded in a horizontally layered medium can be considered as an effective long-wavelength orthorhombic medium. Estimation of Thomsen’s weak-anisotropy (WA) parameters and fracture weaknesses plays an important role in characterizing the orthorhombic anisotropy in a weakly anisotropic medium. Our goal is to demonstrate an orthorhombic anisotropic AVOaz inversion approach to describe the orthorhombic anisotropy utilizing the observable wide-azimuth seismic reflection data in a fractured reservoir with the assumption of orthorhombic symmetry. Combining Thomsen’s WA theory and linear-slip model, we first derive a perturbation in stiffness matrix of a weakly anisotropic medium with orthorhombic symmetry under the assumption of small WA parameters and fracture weaknesses. Using the perturbation matrix and scattering function, we then derive an expression for linearized PP-wave reflection coefficient in terms of P- and S-wave moduli, density, Thomsen’s WA parameters, and fracture weaknesses in such an orthorhombic medium, which avoids the complicated nonlinear relationship between the orthorhombic anisotropy and azimuthal seismic reflection data. Incorporating azimuthal seismic data and Bayesian inversion theory, the maximum a posteriori solutions of Thomsen’s WA parameters and fracture weaknesses in a weakly anisotropic medium with orthorhombic symmetry are reasonably estimated with the constraints of Cauchy a priori probability distribution and smooth initial models of model parameters to enhance the inversion resolution and the nonlinear iteratively reweighted least squares strategy. The synthetic examples containing a moderate noise demonstrate the feasibility of the derived orthorhombic anisotropic AVOaz inversion method, and the real data illustrate the inversion stabilities of orthorhombic anisotropy in a fractured reservoir.     

n—D Radon变换与波动方程偏移

本文把Ｒａｄｏｎ变换公式推广到任意ｎ维的情况。同时结合ｎ维Ｒａｄｏｎ变换和摄动理论提出了一种既能用于地面资料又能用于ＶＳＰ资料的偏移方法。     

平面Radon变换的反演公式

本文借助球面平均法和Ｈｉｌｂｅｒｔ空间算子理论,给出了平面上Ｒａｄｏｎ变换的反演公式,此公式是具体的、构造性的,便于数值计算,或者进一步探讨重建图象的性质。     

地震逆散射波场和算子的谱分解

本文对地震逆散射的研究，旨在于为抑制层间多次波和地震波场多重散射对一次反射干扰效应提供理论依据.这对薄互层地层滤波的高频恢复、保幅弹性反演、衍射地震勘探及海洋地震勘探中的干扰消除皆具重要意义.本文基于上下行波分解及弹性波互易定理，导出横向变速介质条件下线性预测算子的表达式和反射数据的广义谱分解方程. 文中先由上覆地层广义反射透射矩阵的元素定义线性预测算子，并将其表示成一系列单程波算子的线性组合,之后将横向变速介质条件下线性预测方程表达为反射数据与线性预测算子及其逆的乘积. 对该方程的求解可获得上覆地层的线性预测算子，从而可借以求出相应的反射透射算子. 本文先将水平层状介质条件下垂直入射的一维线性预测方程推广到斜入射的情况，以此为参照，导出横向非均匀介质条件下反射数据的地震逆散射广义谱分解方程.文中也揭示了单程波地震逆散射算子、反射透射算子的性态.本文还针对水平层状介质条件，给出斜入射的数值结果.     

3D modelling of VLF radio wave propagation in terrestrial waveguide allowing for localized large-scale ionosphere perturbation

The problem of radio wave propagation allowing for 3D localized lower ionosphere irregularity appears in accordance with the necessity of the theoretical interpretation of VLF remote sensing data. The various processes in the Earth's crust and in space (earthquakes, magnetic storms, sporadic E-layers, lightning induced electron precipitations, rocket launches, artificial ionosphere heating, nuclear explosions, etc.) may cause different power and size ionospheric disturbances. This paper presents a further development of the numerical–analytical method for 3D problem solving. We consider a vector problem of VLF vertical electric dipole field in a plane Earth-ionosphere waveguide with a localized anisotropic ionosphere irregularity. The possibility of lowering (elevating) of the local region of the upper waveguide wall is taken into account. The field components on the boundary surfaces obey the Leontovich impedance conditions. The problem is reduced to a system of 2D integral equations taking into account the depolarization of the field scattered by the irregularity. Using asymptotic (kr⪢1) integration along the direction perpendicular to the propagation path, we transform this system to a system of 1D integral equations. The system is solved in the diagonal approximation, combining direct inversion of the Volterra integral operator and the subsequent iterations. The proposed method is useful for study of both small-scale and large-scale irregularities. We obtained estimates of the TE field components that originate entirely from field scattering by a 3D irregularity.     

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.     

Seismic velocity analysis in the scattering-angle/azimuth domain

Migration velocity analysis is carried out by analysing the residual moveout and amplitude variations in common image point gathers (CIGs) parametrized by scattering angle and azimuth. The misfit criterion in the analysis is of the differential-semblance type. By using angles to parametrize the imaging we are able to handle and exploit data with multiple arrivals, although artefacts may occur in the CIGs and need to be suppressed. The CIGs are generated by angle migration, an approach based on the generalized Radon transform (GRT) inversion, and they provide multiple images of reflectors in the subsurface for a range of scattering angles and azimuths. Within the differential semblance applied to these CIGs, we compensate for amplitude versus angle (AVA) effects. Thus, using a correct background velocity model, the CIGs should have no residual moveout nor amplitude variation with angles, and the differential semblance should vanish. If the velocity model is incorrect, however, the events in the CIGs will appear at different depths for different angles and the amplitude along the events will be non-uniform. A standard, gradient-based optimization scheme is employed to develop a velocity updating procedure. The model update is formed by backprojecting the differential semblance misfits through ray perturbation kernels, within a GRT inverse. The GRT inverse acts on the data, subject to a shift in accordance with ray perturbation theory. The performance of our algorithm is demonstrated with two synthetic data examples using isotropic elastic models. The first one allows velocity variation with depth only. In the second one, we reconstruct a low-velocity lens in the model that gives rise to multipathing. The velocity model parametrization is based upon the eigentensor decomposition of the stiffness tensor and makes use of B-splines.