Efficient solution of large-scale matrix of acoustic wave equations in 3D frequency domain

In 3D frequency domain seismic forward and inversion calculation,the huge amount of calculation and storage is one of the main factors that restrict the processing speed and calculation efficiency.The frequency domain finite-difference forward simulation algorithm based on the acoustic wave equation establishes a large bandwidth complex matrix according to the discretized acoustic wave equation,and then the frequency domain wave field value is obtained by solving the matrix equation.In this study,the predecessor's optimized five-point method is extended to a 3D seven-point finite-difference scheme,and then a perfectly matched layer absorbing boundary condition(PML)is added to establish the corresponding matrix equation.In order to solve the complex matrix,we transform it to the equivalent real number domain to expand the solvable range of the matrix,and establish two objective functions to transform the matrix solving problem into an optimization problem that can be solved using gradient methods,and then use conjugate gradient algorithm to solve the problem.Previous studies have shown that in the conjugate gradient algorithm,the product of the matrix and the vector is the main factor that affects the calculation efficiency.Therefore,this study proposes a method that transform bandwidth matrix and vector product problem into some equivalent vector and vector product algorithm,thereby reducing the amount of calculation and storage.

Decoupling coefficients of dilatational wave for Biot's dynamic equation and its Green's functions in frequency domain

Green＇s functions for Blot＇s dynamic equation in the frequency domain can be a highly useful tool for the investigation of dynamic responses of a saturated porous medium. Its applications are found in soil dynamics, seismology, earthquake engineering, rock mechanics, geophysics, and acoustics. However, the mathematical work for deriving it can be daunting. Green＇s functions have been presented utilizing an analogy between the dynamic thermoelasticity and the dynamic poroelasticity in the frequency domain using the u-p formulation. In this work, a special term ＂decoupling coefficient＂ for the decomposition of the fast and slow dilatational waves is proposed and expressed to present a new methodology for deriving the poroelastodynamic Green＇s functions. The correct- ness of the solution is demonstrated by numerically comparing the current solution with Cheng＇s previous solution. The separation of the two waves in the present methodology allows the more accurate evaluation of Green＇s functions, particularly the solution of the slow dilatational wave. This can be advantageous for the numerical implementation of the boundary element method （BEM） and other applications.

;requency domain wave equation forward modeling using gaussian elimination with static pivoting

频率域波动方程正演是求解一个大型线性稀疏方程组问题,其受到计算效率和内存存储问题的限制。常规的高斯消元法不能满足大型数据的并行计算,本文提出基于静主元消元法(GESP)进行稀疏矩阵LU分解和多炮有限差分正演,该方法不仅提高了稳定性,更有利于单频点内LU分解的分布式并行计算。通过Marmousi模型模拟试验,单频波场和转化到时间域地震剖面的试验表明模拟精度和计算效率得到提高,节约并充分利用内存,为波形反演奠定基础。

Implementation of Time-Scale Transformation Based on Continuous Wavelet Theory

The basic objective of time-scale transformation is to compress or expand the signal in time field while keeping the same spectral properties. This paper presents two methods to derive time-scale transformation formula based on continuous wavelet transform. For an arbitrary given square-integrable function f(t),g(t) = f(t/λ) is derived by continuous wavelet transform and its inverse transform. The result shows that time-scale transformation may be obtained through the modification of the time-scale of wavelet function filter using equivalent substitution. The paper demonstrates the result by theoretic derivations and experimental simulation.

Applying a Step Approach Method in Solving the Multi-Frequency Radiation From a Complex Obstacle

In this paper,a step approach method in the time domain is developed to calculate the radiated waves from an arbitrary obstacle pulsating with multiple frequencies.The computing scheme is based on the Boundary Integral Equation and derived in the time domain;thus,the time-harmonic Neumann boundary condition can be imposed.By the present method,the values of the initial conditions are set to zero,and the approach process is carried forward in a loop from the first time step to the last.At each time step,the radiated pressure on each element is updated.After several loops,the correct radiated pressures can be obtained.A sphere pulsating with a monopole frequency in an infinite acoustic domain is calculated first.This result is compared with the analytical solution,and both of them are in good agreement.Then,a complex-shaped radiator is taken as the studied case.The pulsating frequency of this case is multiple,and the waves propagate in half space.It is shown that the present method can treat multiple-frequency pulsation well,even when the radiator is a complex shape,and a robust convergence can be attained quickly.

An imaging condition for reverse time migration based on wavefield decomposition

Reverse Time Migration(RTM) is a high precision imaging method of seismic wavefield at present,but low-frequency noises severely affect its imaging results.Thus one of most important aspect of RTM is to select the proper noise suppression method.The wavefield characteristics of the Poynting vector are analyzed and the upgoing,downgoing,leftgoing and rightgoing waves are decomposed using the Poynting vector of the acoustic wave equation.The normalized wavefield decomposition cross-correlation imaging condition is used to suppress low-frequency noises in RTM and improve the imaging precision.Numerical experiments using the Mamousi velocity model are performed and the results demonstrate that the upgoing,downgoing,leftgoing and rightgoing waves are well decomposed using the Poynting vector.Compared with the normalized cross-correlation imaging and Laplacian filtering method,the results indicate that the low-frequency noises are well suppressed by using the normalized wavefield decomposition cross-correlation imaging condition.

Finite difference time domain analysis of two-dimensional surface acoustic wave piezoelectric phononic crystals at radio frequency

The finite-difference time-domain （FDTD） method is proposed for analyzing the surface acoustic wave （SAW） propagation in two-dimensional （2D） piezoelectric phononic crystals （PCs） at radio frequency （RF）, and also experiments are established to demonstrate its analysis result of the PCs＇ band gaps. The FDTD method takes the piezoelectric effect of PCs into account, in which periodic boundary conditions are used to decrease memory/time consumption and the perfectly matched layer boundary conditions are adopted as the SAW absorbers to attenuate artificial reflections. Two SAW delay lines are established with/without piezoelectric PCs located between interdigital transducers. By removing several echoes with window gating function in time domain, delay lines transmission function is achieved. The PCs＇ transmission functions and band gaps are obtained by comparing them in these two delay lines. When Aluminum/128°YX-LiNbO3 is adopted as scatter and substrate material, the PCs＇ band gap is calculated by this FDTD method and COMSOL respectively. Results show that computational results of FDTD agree well with experimental results and are better than that of COMSOL.

瞬态横向振动桩的Winkler参数

应用Fourier变换,推导了瞬态横向振动桩桩周土阻力在频率域中的Baranov公式,并得到瞬态横向振动桩Winkler参数k和c的表达式。该式因计算复杂不便于桩瞬态横向振动的计算,笔者应用Hankel函数渐近展开式得到k和c的简化式。对桩的瞬态横向振动的计算表明,简化式的精度满足要求,且大大提高了计算速度。

均匀流中声学速度预测的时域解析公式

采用理论解析和数值计算方法研究均匀流中声学速度预测的时域解析公式。对基于可渗透面原理的声学速度预测,理论分析了对流矢量波动方程的右端源项,结果显示部分对声学速度传播有贡献的可渗透面源被忽略。基于对流矢量波动方程时域积分公式V1A-M,进一步考虑丢失可渗透面源的贡献,发展了解析公式CV1A-M。对均匀平均流中静止与旋转的单极子和偶极子,利用公式V1A-M和CV1A-M数值预测了点源辐射的声学速度。公式CV1A-M数值解与理论值吻合一致,公式V1A-M因可渗透面载荷源部分丢失,数值解与理论值差异明显。均匀流中偶极子点源声辐射不会诱发涡波扰动,线化欧拉方程能够描述声学速度与声学压力间关系。

起伏地表条件下基于无网格有限差分方法的频率域声波波动方程数值模拟

随着地震勘探对象日趋复杂,能精确描绘起伏地表、复杂地质构造的数值模拟方法显得十分重要。本文将基于无网格有限差分方法发展至起伏地表条件下频率域数值模拟中:无网格节点理论上能够适应任意地表和地质构造;有限差分方法具有计算效率高、内存占用小、模拟精度较高等优点。因此基于无网格的有限差分方法既保证了高效率又适用于不规则地表。此外,本文将完全匹配层(PML)吸收边界条件引入基于无网格有限差分频率域数值模拟中以压制边界反射,并对比了简化PML、经典PML以及复频移PML方法压制边界反射的效果;随后将起伏地表条件下压制边界反射效果较好的复频移PML方法应用到基于无网格与规则网格数值模拟的对比中,其中波场快照、地震记录以及地震子波波形对比结果都验证了基于无网格有限差分方法数值模拟的有效性。最后,将本文数值模拟方法应用到起伏地表的Marmousi模型与实际资料中,进一步验证了该方法的有效性与应用性。

气泡软化效应下的非线性声场

为研究气泡的存在对介质中非线性声场的影响,本文利用波动方程和Rayleigh-Plesset方程建立的有限差分和有限体积模型,分别求解了在单频波和双频波激励下含气泡水介质中任意时刻的声压分布。计算表明:在非线性条件下,沿声波传播方向长度为半波长的介质内存在共振频移,即介质发生了软化;含气泡水介质的非线性共振频移与振源的振幅有关,在一定范围内,频移随着振幅的增大而增大;并证明了介质中的软化是由气泡平均体积增大而引起的。本研究为提高参量阵声发射和接收效率等提供了一种思路。

基于小波变换的木材声发射驻波信号特征研究

为研究木材损伤断裂时的声发射(AE)信号所激发的驻波信号特征与木材固有特性之间的关系,采用薄木条折断的方式产生AE源,在小波变换的基础上分析驻波频率,并计算纵波传播速率,依据弹性波理论计算出木材顺纹弹性模量(MOE)。首先,在2种不同长度的木材试件一端分别加工出8根80 mm×10 mm的薄木条,通过外加冲击力折断木条以产生AE源,通过放置在试件端面的2个传感器采集原始AE信号,采样频率设定为500 kHz。然后,根据驻波特性确定原始信号的驻波阶段,进而对该阶段AE信号进行4层小波分解,依据分解后信号的时频域特征析取驻波信号波形。最后,依据驻波产生原理计算纵波传播速率,并结合弹性波理论计算试件的MOE。结果表明,拉伸试验测得樟子松和榉木试件的MOE分别为9.30 GPa和11.63 GPa,800 mm樟子松和榉木试件通过驻波计算所得MOE分别为9.37 GPa和12.34 GPa,与实测MOE的误差分别为0.75%和5.24%;600 mm的樟子松和榉木试件通过驻波计算所得MOE分别为9.31 GPa和11.81 GPa,与实测MOE的误差分别为0.10%和1.55%。

EQUIVALENCE BETWEEN EXACT INTERNAL CONTROLLABILITY OF THE KIRCHHOFF PLATE-LIKE EQUATION AND THE WAVE EQUATION

When the rotatory inertia is taken into account, vibrations of a linear plate can be described by the Kirchhoff plate equation. Consider this equation with locally distributed control forces and some boundary condition which is the simply supported boundary condition for a rectangular plate. In this paper, the authors establish exact controllability of the system in terms of the equivalence to exact internal controllability of the wave equation, by means of a frequency domain characterization of exact controllability introduced recently in [11].

溶洞模型逆散射成像技术

针对常规叠前深度偏移成像技术难以将溶洞型储集层的多次绕射波完全归位、难以准确成像溶洞分布这一难题,引入了频率域声波逆散射成像方法。该方法直接从炮集地震资料出发,利用散射记录的逆时传播波场与入射波场褶积得到扰动速度的分布,将速度模型迭代更新与构造成像融为一体。溶洞模型数值模拟和物理模拟资料的处理效果表明,即使初始速度模型存在12%的误差,该方法仍可以将溶洞产生的多次绕射波较好地归位,较准确地将溶洞的形状、位置成像。同时,对常规偏移成像方法和逆散射成像方法影响溶洞成像精度的主要因素进行了定量分析,为提高溶洞型碳酸盐岩油气藏勘探精度提供了理论指导。

运动声源多普勒畸变信号的一种时域校正方法

运动声源信号中广泛存在多普勒效应的影响,为了解决其引起的频移、频带扩展问题,提出了一种时域信号校正方法。该方法首先从运动关系和莫尔斯声学理论出发推导得出离散插值时间序列和幅值还原公式,然后对畸变信号进行幅值还原和插值拟合。仿真信号和实验信号处理结果表明多普勒效应带来的频移和频带扩展问题得到了较好的解决,畸变信号得到了有效校正。该方法在时域内对离散信号进行单点操作,与以往的频域处理方法相比具有无需已知声源特征频率、不受频带混叠因素影响的优点,该方法能够处理多频信号,是一种简单有效的多普勒效应消除方法。

超声兰姆波的时频分析

声-超声方法在被测板件厚度较小时,板中传播的超声波模式主要是多模式的兰姆波。由于兰姆波的复杂性,声-超声方法回避了对信号进行进一步深入的分析,而采用所谓“应力波因子”对材料中能量耗散进行相对度量。当多模式与频散同时存在时,信号的识别异常困难,本文采用时-频分析的方法对声-超声信号进行了分离,只需一次测量就可满足分析要求。通过对兰姆波信号的伪Wigner-Ville分布进行时频重排和图像特征提取,能够精确地获得兰姆波的局部群速度频散曲线。

利用单频双程波动方程计算初至走时及其振幅

通过在频率域双程波动方程模拟算法中加入一个复数频率(实部表示频率,虚部表示衰减因子)压制地震波初至走时之后的能量,从而把初至走时及其振幅的计算问题转换为单一频率波场中最大能量走时和振幅的拾取问题,然后利用单一频率域波场的相位项和振幅项分别计算初至走时及其振幅.本文还提出利用参数分析方法求取最优的复数频率,并给出数值计算例子,将本方法的计算结果与有限差分程函方程初至走时和最大能量走时振幅进行比较,结果表明,该方法具有适应于任意复杂介质和多炮多接收点走时和振幅的计算.

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

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

频率-空间域有限差分法叠前深度偏移

为了处理横向强变速介质中的深度成像问题 ,本文提出一种基于共炮道集的优化系数的傍轴近似方程叠前深度偏移算子 ,并在基于反射系数估算的成像条件下 ,可实现叠前深度偏移成像 .该算子具有方程阶数低且能对陡倾角成像的特征 ,并采用有限差分法波场延拓 ,能适应速度场的任意变化 .当在频率 -空间域进行计算时 ,相对于纯粹的时间 -空间域有限差分算法有计算效率高、成像方便的优点 .脉冲响应测试和对Marmousi模型进行的叠前深度偏移结果表明 。