共炮检距高斯波束偏移的一种快速实现算法

与共炮高斯波束偏移相比,共炮检距高斯波束偏移具有直接抽取炮检距域共成像点道集的优势.过去,共炮检距高斯波束偏移以损失成像精度的代价采用最速下降法来降低积分的维数,从而提高计算效率.但经过最速下降近似简化的偏移公式仍是频率域的,需要在每个频点进行计算.为此,本文提出一种快速实现算法来避免采用最速下降法.本文通过分析一个水平层状速度模型的偏移过程和Marmousi速度模型的成像结果来检验不同插值方法对快速实现算法的成像精度和计算效率的影响,并建议采用二维三次卷积插值方法.同时本文在Marmousi速度模型下验证了快速实现算法相对于最速下降法在成像精度和计算效率上的优势.此外,本文将采用二维三次卷积插值的快速实现算法应用于Sigsbee2A模型并获得了清晰的盐下图像.     

Numerical viscoelastic modelling by the spectral Laguerre method

Seismic wave propagation in a viscoelastic media can be described by a system of integro-differential equations. The solution of such equations requires special methods when using finite-difference techniques in the time domain. In the frequency domain, the integral terms are represented by complex elastic parameters. This paper presents an efficient algorithm for viscoelastic modelling based on the integral Laguerre transform for the approximation of temporal derivatives and for the calculation of convolution integrals. For the calculation of spatial derivatives, it is possible to use various methods: finite-difference and finite-element techniques, spectral and pseudo-spectral methods. We then obtain a system of algebraic equations with a matrix independent of the parameter m , i.e. the degree of the Laguerre polynomials. In this case, only the right-hand side of the system has recurrent dependence on the parameter m , which is an analogue of the temporal frequency in the frequency domain. The obtained system with a large number of right-hand sides can be solved using fast methods, where the matrix is transformed only once, as opposed to the frequency-domain approach, when the matrix is transformed for each temporal frequency.     

Simultaneous perturbation stochastic approximation for tidal models

The Dutch continental shelf model (DCSM) is a shallow sea model of entire continental shelf which is used operationally in the Netherlands to forecast the storm surges in the North Sea. The forecasts are necessary to support the decision of the timely closure of the moveable storm surge barriers to protect the land. In this study, an automated model calibration method, simultaneous perturbation stochastic approximation (SPSA) is implemented for tidal calibration of the DCSM. The method uses objective function evaluations to obtain the gradient approximations. The gradient approximation for the central difference method uses only two objective function evaluation independent of the number of parameters being optimized. The calibration parameter in this study is the model bathymetry. A number of calibration experiments is performed. The effectiveness of the algorithm is evaluated in terms of the accuracy of the final results as well as the computational costs required to produce these results. In doing so, comparison is made with a traditional steepest descent method and also with a newly developed proper orthogonal decomposition-based calibration method. The main findings are: (1) The SPSA method gives comparable results to steepest descent method with little computational cost. (2) The SPSA method with little computational cost can be used to estimate large number of parameters.     

带先验知识的波阻抗反演正则化方法研究

针对波阻抗反演中存在的不适定性问题,本文提出了一种带先验知识的正则化重开始共轭梯度法.该方法的内层循环采用修改的共轭梯度法,并使用重开始技巧;外层循环使用Morozov偏差准则作为停机准则.正则参数的选取采用连续几何选取法.克服了传统共轭梯度法迭代不足或迭代过度的缺点,将迭代步数控制在了合适的范围,使算法能够更快速更准确的收敛.同时考虑了用最速下降法计算先验解和对解施加非均一的规范约束.通过理论模型试算和实际资料处理,并与共轭梯度法进行对比,表明该算法具有精度高、抗病态能力强,运算速度快的优点,具有实用性.     

多频可控源电磁法三维有理函数Krylov子空间模型降阶正演算法研究

本文采用有理函数Krylov子空间模型降阶算法实现了同时求解多频可控源电磁法三维正演响应的快速计算.首先采用基于Yee氏交错网格的拟态有限体积法实现控制方程的空间离散,将任意频率的电场响应表示为关于频率参数的传递函数.采用有理函数Krylov子空间算法求解该传递函数.针对构建m维有理函数Krylov子空间需要求解m次（几十到上百）关于有理函数极点和离散控制方程系数矩阵的线性方程组的问题,本文提出采用单个重复极点的有理函数Krylov子空间模型降阶算法,结合直接法求解器PARDISO,采用Gram-Schmidt方法,只需要1次系数矩阵分解和m次矩阵回代即可实现有理函数Krylov子空间的构建,极大地减少了计算量.针对最优化有理函数极点选取问题,本文根据传递函数的有理函数Krylov子空间投影算法的误差分析理论,引入关于单个重复极点的收敛率函数,通过求解有理函数的最大收敛率直接给出最优化的单个重复极点公式.最终实现了不同发射频率的可控源电磁法三维正演响应的快速计算.分别计算了典型层状模型多发射频率的CSAMT和海洋CSEM的正演响应,通过与解析解的对比验证了本文算法在多发射频率正演的计算精度和计算效率;并通过一个三维海洋CSEM勘探设计最优化发射频率和接收区域选取的例子进一步说明本文算法的优点.     

Numerical behaviour of the downward continuation of gravity anomalies

The numerical results of downward continuation (DWC) of point and mean gravity anomalies by the Poisson integral using point, single mean, and doubly averaged kernel are examined. Correct evaluation of the integral in its innermost zone is a challenging task. To avoid instabilities, an analytical planar approximation is used in the innermost integration zone. In addition it is shown that the single mean mode has the minimum discretization error. Downward continuation of point and mean anomalies by singly and doubly averaged kernel are the same mean anomalies on the geoid.     

瞬变电磁虚拟波场二阶Born近似成像算法

依据瞬变电磁扩散场与虚拟波动场之间的数学表达式,通过波场反变换得到虚拟波场,进而借用地震勘探上成熟的成像方法定位地质界面的位置以及形态.本文首先将时间域波动方程转换为频率域波动方程,将总波动场分为背景场和散射场、总波速分为背景速度和扰动速度,然后运用Green定理得到了散射场关于地质界面速度扰动量的表达式,并使用Bor...     

Comparison of methods for estimating earth resistivity from airborne electromagnetic measurements

Earth resistivity estimates from frequency domain airborne electromagnetic data can vary over more than two orders of magnitude depending on the half-space estimation method used. Lookup tables are fast methods for estimating half-space resistivities, and can be based on in-phase and quadrature measurements for a specified frequency, or on quadrature and sensor height. Inverse methods are slower, but allow sensor height to be incorporated more directly. Extreme topographic relief can affect estimates from each of the methods, particularly if the portion of the line over the topographic feature is not at a constant height above ground level. Quadrature–sensor height lookup table estimates are generally too low over narrow valleys. The other methods are also affected, but behave less predictably. Over very good conductors, quadrature–sensor height tables can yield resistivity estimates that are too high. In-phase–quadrature tables and inverse methods yield resistivity estimates that are too high when the earth has high magnetic susceptibility, whereas quadrature–sensor height methods are unaffected. However, it is possible to incorporate magnetic susceptibility into the in-phase–quadrature lookup table. In-phase–quadrature lookup tables can give different results according to whether the tables are ordered according to the in-phase component or the quadrature component. The rules for handling negative in-phase measurements are particularly critical. Although resistivity maps produced from the different methods tend to be similar, details can vary considerably, calling into question the ability to make detailed interpretations based on half-space models.     

复杂地表条件下叠前菲涅尔束偏移方法

高斯束偏移虽然克服了Kirchhoff偏移不能处理多波至和单程波动方程偏移不能对陡倾构造准确成像的问题,但在复杂地表条件下其偏移精度取决于所选择的初始束宽度,即当初始宽度较小时,近地表成像精度较高,但此时中深层成像质量较差;反之当初始宽度较大时,中深层成像质量提高,但近地表成像精度降低.针对高斯束偏移中深层和浅层成像精度的矛盾,本文发展了一种适用于陆地复杂地表条件的叠前菲涅尔束偏移方法.基于惠更斯-菲涅尔原理,本文首先给出了菲涅尔束的概念及其表征的格林函数,并采用有效邻域波场近似理论和反褶积成像条件,导出了复杂地表条件下叠前保幅深度偏移公式.最后,针对常规旁轴射线追踪中的数值噪音,给出了一种压制策略.同高斯束偏移相比,本文方法不仅解决了中深层和浅层成像精度的矛盾,而且提高了复杂地表条件下平面波的分解精度,使得偏移结果更加准确可靠.典型的模型算例验证了本文方法的有效性和稳健性.     

随机结构动力可靠度估计算法的对比分析

对比分析了随机结构动力可靠度计算的三种估计算法.渐进展开法是基于Laplace算法对概率积分进行渐进估计的,此法通过计算最大被积分式值对应点,并将其代入概率积分的渐进估计表达式求解失效概率.由于概率积分的主要贡献来自于最大被积分式值对应点的周围,因此本文的重要抽样法假定重要抽样函数的最大似然值等于最大被积分式值对应点值.极值分布-泰勒展开法首先通过结构随机参数的极值分布函数给出失效概率的表达式,随后利用泰勒展开法对失效概率进行估计,其中采用中心差分法对极值分布函数的梯度进行估算.最后应用三种算法和Monte Carlo法对受高斯白噪声激励作用的单自由度随机结构进行了计算,结果表明三种方法不但运算简便,而且对比Monte Carlo法计算效率有显著提高.     

垂直地震力作用下弹性半空间上扁壳基础的解析解

用布希涅斯克定义的弹性半空间内的垂直位移包括两项积分,除了积分号前面系数的差别之外,第一项积分是单层位势而第二项积分为双层位势。若扁壳基础是正高斯曲率的几何曲面,则壳底与半空间表面间的挤压强度就是半空间表面作用的分布垂直荷载。当越过边界时,双层势位的函数值和单层势位的法向导数值发生跳跃。利用这些性质,本文得出布希涅斯克积分的反演公式,从而避开要求解偏微分—积分方程组的巨大数学困难而易于得出解析解。以椭园抛物面扁壳为例说明本文方法的应用。     

基于压缩感知的Curvelet域联合迭代地震数据重建

由于野外采集环境的限制,常常无法采集得到完整规则的野外地震数据,为了后续地震处理、解释工作的顺利进行,地震数据重建工作被广泛的研究.自压缩感知理论的提出,相继出现了基于该理论的多种迭代阈值方法,如CRSI方法（Curvelet Recovery by Sparsity-promoting Inversion method）、Bregman迭代阈值算法（the linearized Bregman method）等.CSRI方法利用地震波形在Curvelet的稀疏特性,通过一种基于最速下降的迭代算法在Curvelet变换域恢复出高信噪比地震数据,该迭代算法稳定,收敛,但其收敛速度慢.Bregman迭代阈值法与CRSI最大区别在于每次迭代时把上一次恢复结果中的阈值前所有能量都保留到本次恢复结果中,从而加快了收敛速度,但随着迭代的进行重构数据中噪声干扰越来越严重,导致最终恢复出的数据信噪比低.综合两种经典方法的优缺点,本文构造了一种新的联合迭代算法框架,在每次迭代中将CRSI和Bregman的恢复量加权并同时加回本次迭代结果中,从而加快了迭代初期的收敛速度,又避免了迭代后期噪声干扰的影响.合成数据和实际数据试算结果表明,我们提出的新方法不仅迭代快速收敛稳定,且能得到高信噪比的重建结果.     

基于Sentinel-1A数据的2017年伊拉克哈莱卜杰MW7.3地震同震形变场分析及断层滑动分布反演

本文首先利用二轨法对欧洲空间局Sentinel-1A雷达卫星影像进行差分干涉处理,获取了覆盖2017年伊拉克哈莱卜杰（Halabjah）M_W7.3地震震区的同震形变场,结果表明:哈莱卜杰地震造成的地表形变影响范围约为60 km×70 km,形变场基本沿扎格罗斯主前缘断层展布;形变场的西南盘呈现隆升趋势,最大视线向形变值为88 cm,东北盘呈现下降趋势,最大视线向形变值为37 cm;隆升形变值远大于沉降值,反映出发震断层以逆冲运动为主的特征。然后基于弹性半平面空间矩形位错模型,分别采用多峰值粒子群算法和最速下降法确定了发震断层的几何参数和滑动分布结果。反演结果显示发震断层以逆冲运动为主,兼少量右旋走滑运动,最大滑动量为3.34 m,释放的地震矩为1.68×1020 N·m （M_W7.4）,与地震学的研究结果一致。      

Computation of ray integrals and ray amplitudes in radially symmetric media

Summary A new approximation of the velocity-depth distribution in radially symmetric media is suggested. This approximation guarantees the continuity of velocity and its first and second derivatives, and does not generate false low-velocity layers. It removes false anomalies from the amplitude-distance curve and considerably increases its stability. The evaluation of ray integrals and ray amplitudes using this velocity-depth approximation does not require the computation of any transcendental function and is, therefore, very fast. Numerical examples are presented.     

Computation of Ray Integrals and Ray Amplitudes in Radially Symmetric Media

A new approximation of the velocity-depth distribution in radially symmetric media is suggested. This approximation guarantees the continuity of velocity and its first and second derivatives, and does not generate false low-velocity layers. It removes false anomalies from the amplitude-distance curve and considerably increases its stability. The evaluation of ray integrals and ray amplitudes using this velocity-depth approximation does not require the computation of any transcendental function and is, therefore, very fast. Numerical examples are presented.     

THE CAGNIARD METHOD IN COMPLEX TIME REVISITED1

The Cagniard-de Hoop method is ideally suited to the analysis of wave propagation problems in stratified media. The method applies to the integral transform representation of the solution in the transform variables (s, p) dual of the time and transverse distance. The objective of the method is to make the p-integral take the form of a forward Laplace transform, so that the cascade of the two integrals can be identified as a forward and inverse transform, thereby making the actual integration unnecessary. That is, the exponent (–sw(p)) is set equal to –sτ, with τ varying from some (real) finite time to infinity. As usually presented, the p-integral is deformed onto a contour on which the exponent is real and decreases to –∞ as p tends to infinity. We have found that it is often easier to introduce a complex variable τ for the exponent and carry out the deformation of contour in the complex τ-domain. In the τ-domain the deformation amounts to ‘closing down’ the contour of integration around the real axis while taking due account of singularities off this axis. Typically, the method is applied to an integral that represents one body wave plus other types of waves. In this approach, the saddle point of w(p) that produces the body wave plays a crucial role because it is always a branch point of the integrand in the τ-domain integral. Furthermore, the paths of steepest ascent from the saddle point are always the tails of the Cagniard path along which w(p) →∞. That is, the image of the pair of steepest ascent paths in the p-domain is a double covering of a segment of the Re τ-axis in the τ-domain. The deformed contour in the p-domain will be the only pair of steepest ascent paths unless the original integrand had other singularities in the p-domain between the imaginary axis and this pair of contours. This motivates the definition of a primary p-domain, i.e. the domain between the imaginary axis and the steepest descent paths, and its image in the τ-domain, the primary τ-domain. In terms of these regions, singularities in the primary p-domain have images in the primary τ-domain and the deformation of contour on to the real axis in the τ-domain must include contributions from these singularities. This approach to the Cagniard-de Hoop method represents a return from de Hoop's modification to Cagniard's original method, but with simplifications that make the original method more tractable and straightforward. This approach is also reminiscent of van der Waerden's approach to the method of steepest descents, which starts exactly the same way. Indeed, after the deformation of contour in the τ-domain, the user has the choice of applying asymptotic analysis to the resulting ‘loop’ integrals (Watson's lemma) or continuing to obtain the exact, although usually implicit, time-domain solution by completing the Cagniard-de Hoop analysis. In developing the method we examine the transformation from a frequency-domain representation of the solution (ω) to a Laplace representation (s). Many users start from the frequency-domain representation of solutions of wave propagation problems. In this case issues arising from the movement of singularities under the transformation from ω to s must be considered. We discuss this extension in the context of the Sommerfeld half-plane problem.     

利用高斯求积和连分式展开计算电磁张量格林函数积分

针对于利用快速汉克尔变换求解电磁张量格林函数中含有贝塞尔函数积分项时，积分核函数不总是随贝塞尔函数（Jv(λr)）参数λ的增大而趋于收敛，从而产生错误结果的问题，我们利用一种直接数值积分——结合连分式展开的高斯求积方法进行计算，通过对不同模型的试算表明这种方法总能够保证电磁张量格林函数积分的正确计算.     

水平层状介质中基于DTA的三维电磁波逆散射快速模拟算法

为提高水平层状介质中三维电磁波散射和逆散射数值模拟的效率,在对角张量近似(DTA)的基础上根据不同回代方式得到了求解积分方程的DTA1和DTA2两种近似. 这两种近似可以作为计算积分方程稳定型双共轭梯度快速Fourier变换(BCGS-FFT)算法的初始猜测值和预条件因子,从而形成效率更高的混合DTA-BCGS算法. 散射实例说明了DTA2的高精度和混合DTA-BCGS算法尤其是混合DTA2-BCGS算法的高效率. 由于DTA2近似程度更高,将DTA2与变型Born迭代反演方法(DBIM)相结合形成了一种对三维异常体进行重构的快速电磁波逆散射技术. 文中的逆散射实例说明所开发的逆散射技术对重构水平层状介质中的任意三维异常体是非常有效的.     

Determination of porosity and saturation using seismic waveform inversion

Characterization of a reservoir model requires determination of its petrophysical parameters, such as porosity and saturation. We propose a new method to determine these parameters directly from seismic data. The method consists of the computation and inversion of seismic waveforms. A high frequency method is presented to model wave propagation through an attenuative and dispersive poroelastic medium. The high frequency approximation makes it possible to efficiently compute sensitivity functions. This enables the inversion of seismic waveforms for porosity and saturation. The waveform inversion algorithm is applied to two laboratory crosswell datasets of a water saturated sand. The starting models were obtained using travel time tomography. The first dataset is inverted for porosity. The misfit reduction for this dataset is approximately 50%. The second dataset was obtained after injection of a nonaqueous-phase liquid (NAPL), possibly with some air, which made the medium more heterogeneous. This dataset was inverted for NAPL and air saturation using the porosity model obtained from the first inversion. The misfit reduction of the second experiment was 70%. Regions of high NAPL and high air saturation were found at the same location. These areas correlate well with the position of one of the injection points as well as regions of higher NAPL concentrations found after excavation of the sand. It is therefore possible to directly invert waveforms for pore fluid saturation by taking into account the attenuation and dispersion caused by the poroelasticity.     

Stationary approximation of the magnetotelluric field in a halfspace with a vertical boundary and three-dimensional perturbing body

Summary Based on the generalized potential of a double layer, integral formulae have been derived for calculating the stationary approximation of the magnetotelluric field in a halfspace divided by a vertical boundary into two quarterspaces, one of which contains a three-dimensional perturbing body. The appropriate boundary integral equation and other surface integrals have been computed for a perturbing body in the shape of the three-dimensional prism located at the vertical boundary, or in contact with the said boundary. The exciting electrical field is assumed to be homogeneous and perpendicular to the vertical boundary. Isoline graphs of the electrical and magnetic fields on the surface of the halfspace have been plotted and their anomalies are discussed.