【《弹性波全波形反演理论基础综述》8100字】

弹性波全波形反演理论基础综述目录TOC\o"1-3"\h\u5234弹性波全波形反演理论基础综述 1141001.1弹性波时间域正演理论基础 1167551.1.1弹性波波动方程的差分格式 187551.1.2PML边界条件 6154331.1.3数值稳定性 9177181.2弹性波频率域反演理论基础 10164631.1.1弹性波波场分离 10315081.1.2梯度的计算 11157471.1.3迭代步长的计算 141.1弹性波时间域正演理论基础1.1.1弹性波波动方程的差分格式在进行弹性波正演模拟时，常常需要对弹性波方程进行有限差分计算，有限差分法的基本原理是应用泰勒级数对方程进行展开近似，将弹性波方程中的微分项用离散差分项替换，并进行差分系数的求解。本篇论文给出了较为详细的推导过程，如下述所示。（1）一阶速度-应力弹性波方程通过对弹性波动力学的学习，我们可以通过三大方程整理并推导出弹性波波动方程。其中，几何方程描述了物体位移与物体应变的关系,其表达式为： εxx=本构方程又称为广义胡克定律，是弹性波动理论的一个基本点，即给出了物体应力与应变之间的单值线性关系,具体表达式如下：σ（2-2）不同介质条件下，弹性系数矩阵可以化简为不同的形式，在各向同性介质中，可以对弹性系数矩阵C进行如下形式的简化：C=其中C11=C22=C33根据牛顿第二定律对物体的受力进行分析，则可以得到运动微分方程，表示形式如下: ∂σxx通过几何方程（2-1）、本构方程（2-2）和运动平衡微分方程（2-3）三者之间的关系，在只考虑二维的条件下，一阶速度-应力弹性波方程的表示形式如下： ∂vx在上式（2-4）中，vx,vz代表速度分量，σxx（2）弹性波方程时间上2M阶差分近似在计算弹性波方程时间上2M阶差分近似时，对上式（2-4）进行泰勒展开，以vx为例，将vx(t+v（2-5）同样的方法对vz(t+Δtv（2-6）令t'=t−Δt2，同样利用泰勒展开的思想，带入上式（2-5）στ（2-7）则上式（2-5）、（2-6）以及（2-7）为弹性波方程时间上2M阶差分近似格式。令2M=2时，即可求得弹性波方程时间上2阶差分近似公式，具体表示形式如下： vx(t+（3）弹性波方程空间上2N阶差分近似在计算弹性波方程空间上2N阶差分近似时，对一个给定的函数f（x），如果它的各阶导数都是关于变量x的单值连续函数，则它的泰勒展开可以表示为如下形式[25]：∂f（2-9）对上式（2-9）进行整理，则可以改写为：∂f（2-10）在上式（2-10）中，Δx为网格间距，an、C∂f（2-11）对上式（2-11）进行合并同类项整理可得：∂f（2-12）对等式（2-12）左右两侧通过待定系数法进行整理，可得关于系数的变量方程：C1(N)将上式方程组（2-13）改写成矩阵形式如下：135当N的取值不同时，求解上式方程组（2-20）可得到不同的权系数Cn则弹性波方程空间上2N阶差分近似为：∂f（2-15）（4）弹性波波动方程的差分格式通过上述的推导过程，可以求得弹性波方程的时间上2阶差分近似和空间上2N阶差分近似，在求出弹性波方程的时间上2阶差分近似和空间上2N阶差分近似之后，进而求解弹性波方程的差分格式，设定相应的交错差分网格如图2-1所示，则时间上2阶差分精度、空间上2N阶差分精度的差分格式计算方法如下所述。图2-1O(Δt2+Δ设Ui,jk+1/2、Vi+1/2,j+1/2k+1/2、Ri+1/2,jk、UVRTH（2-16）1.1.2PML边界条件在实际生产生活中，实际地下介质为一个半无限的空间，对地下介质进行数据采集时，实际数据量巨大，通常情况下都以T为单位，普通的计算机内存很难满足需求。因此，在进行地震波模拟时，需要对半无限空间进行人工截断，引入一个边界，但是边界的出现使得地震波在传播到边界上时出现反射现象，形成干扰波。如何引入人工边界，在截断无限空间的同时也能保证尽量减小地震波的反射，削弱不必要的干扰波，是目前需要攻克的难题。1994年，Berenger利用电磁波的传播特点，提出了一种高效的完美匹配层吸收边界条件PML。通过文献阅读与分析，本篇论文选择了效果相对较好的裂化PML边界条件，该方法最早是由Berenger提出，后来，王守东[28]将该思想引入到声波方程中,得出了令人满意的数值模拟结果。PML示意图如下图2-2所示：图2-2PML示意图（据王守东[28]）下面给出在各向同性介质中弹性波一阶速度-应力方程的PML吸收边界条件的推导过程，对上式方程组（2-4）中速度、应力进行变量分离，在二维的条件下，每一个变量都可分为如下的两部分[29]：vx上式（2-17）中上标x代表只与x的空间导数有关，上标z代表只与z的空间导数有关，根据上式（2-17）对弹性波方程组（2-4）进行分裂，同时，引入与x、z方向导数有关的衰减因子dx、d(∂∂t+d在式（2-18）中，dx、dz分别代表与x、z方向导数有关的衰减因子，对于衰减因子的选择十分重要，通常dd其中，xδ区域1：d区域2：d区域3：d设Uxi,jk+1/2、Vxi+1/2,j+1/2k+1/2、RxUUVVRRzTHH（2-19）1.1.3数值稳定性时间域弹性波方程有限差分解法是以一组以差分方程形式近似代替连续方程形式的弹性波方程，即将弹性波方程中的微分项用离散差分项替换，以差分方程组的解代替波动偏微分方程的解，差分格式对最终的计算结果有重要的影响。因此，如何在有限差分数值模拟过程中确保稳定性是一项难题。董良国[27]等人对一阶弹性波波动方程进行了详细的研究，本篇论文引用了董良国给出的稳定性条件，由此我们得出二维稳定性条件：0≤（2-20）其中，在上述方程（2-20）中;LLd=当Δx=Δz=l，最大纵波速度为vmax，则精度为O(Δm=1(2-21)在上式(2-21)中，ρ表示密度，Cn1.2弹性波频率域反演理论基础1.1.1弹性波波场分离由纵波速度、横波速度与拉梅常数的关系λ=ρ(v ∂v为构造等价方程式，引入混合波场新变量U、P波波场新变量Up和S波波场新变量Us，即如下表现形式[32]：U={则可得到混合波场U、P波波场Up和S波波场Us之间的关系： vx=v对一阶速度-应力弹性波方程（2-22）进行分解，即可得到纵横波波场分离后的弹性波方程：ρ（2-24）1.1.2梯度的计算确定搜索方向是弹性全波形反演中至关重要的一步。通常情况之下，确定搜索方向主要有两种方法：梯度引导类方法和牛顿类方法。其中，在梯度引导类方法中，较为常用的方法就是最速下降法和共轭梯度法，相比之下此类方法计算量相对较小、便于实现，但是此类方法只能进行一阶收敛，相对而言收敛精度较低。牛顿类方法中较为常用的方法相对较多，包括牛顿法、拟牛顿法、高斯牛顿法、BFGS法、L-BFGS法等等，此类方法与梯度引导类方法刚好相反，牛顿法计算繁琐、计算量大，但是它的收敛精度较高，具有二阶收敛精度。本篇论文简单阐述其中的两种方法，即牛顿法和共轭梯度法：（1）牛顿法在野外实际地质资料采集过程中，常常会遇到各种各样的问题，例如，地表的噪音会对数据的观测产生干扰以及由于观测系统的不完备使得测量数据缺失等问题，都会影响到反演过程，会使得反演陷入局部极值。因此，如何解决波动方程在反演过程中的所得方程解的不唯一性是目前困惑诸多学者的难题。一般情况下，非线性问题的线性化求解方法分两类，其中一类为全局优化算法，另一类为局部优化算法。全局优化算法的最大优点在于其对目标函数的依赖程度不高，但其缺点也较为明显，就是在反演过程中会出现随机性，计算量也极大，因此，此类方法在实际生产生活中应用率不高。局部优化算法其计算量相对较小，因此，在实际应用中应用范围较为广泛，下面重点介绍牛顿法，即利用牛顿法将非线性问题线性化。在全波形反演方法中，一般定义误差向量为：δdm=dobsmC(m)=1在上述方程（2-25）中，*代表共轭转置。牛顿法认为，目标函数在初始模型附近满足二次型，因此将其进行泰勒级数展开：C（2-26）在上式（2-26）中，m0代表初始模型，δm代表模型扰动量，则真实模型m表示为m= ∂Cm∂m令目标函数的一阶导数为零，即方程（2-27）等于零时，都过整理可解得扰动模型向量δm为： δm=−∂2在上式（2-28）中，令目标函数的梯度为：∇mC=∂Cδm=−H−1将上式（2-29）带入m=mmi+1=（2）共轭梯度法共轭梯度法是在最速下降法的基础之上发展而来的一种非线性的梯度类局部寻优方法，共轭梯度法能够很好地解决最速下降法产生的锯齿效应，同时又能有效地提高收敛精度。共轭梯度法是利用两次梯度方向来构建出一个更新方向，因此，在实际的全波形反演计算中，应用更为广泛。共轭梯度法是迄今为止较为常用且有效的一种方法，在每次迭代的过程中，利用当前位置的最速下降方向生成共轭方向。对于一个给定的函数f(x)，其一般迭代公式为[33]： xk+1=在上式（2-31）中，gk=∇fxPRP公式：βFR公式：βi对牛顿法、最速下降法和共轭梯度法进行分析对比，总结三种方法的优缺点，可以得出结论：牛顿法的收敛精度要求较高，具有二阶收敛精度，但计算量较大，尤其是Hessian矩阵及其逆矩阵的求取往往较为繁琐。共轭梯度法相较与牛顿法，计算效率高、便于实现，并且与最速下降法相比，收敛精度较高。考虑到弹性波全波形反演本身计算量已十分巨大，所以在搜索方向方面，本篇论文选择的是计算高效且便于实现的共轭梯度类方法，将共轭梯度的方向作为搜索方向指引模型更新。同时，考虑到各方面的影响因素，本篇论文在进行梯度处理的过程中，引入衰减补偿因子，以同时满足全波形反演过程中压制浅层噪音和补偿深层能量。进行梯度处理后的速度迭代公式可表示为如下形式：mk+1=在上式（2-32）中，αk表示第k次的迭代步长，gk表示第k次的梯度，β=exp(在上式（2-33）中，f代表频率，z代表深度，C是常数，γ为衰减因子，不同情况下取值不同，一般都在0和1之间进行取