【《勒夫波频散曲线计算编程实现及模型试算案例分析》7300字】

勒夫波频散曲线计算编程实现及模型试算案例分析目录TOC\o"1-3"\h\u25509勒夫波频散曲线计算编程实现及模型试算案例分析 1177891.1勒夫波频散曲线正演 121521.1.1Knopoff算法 2174051.1.2广义反射-透射系数算法 328641.1.3Haskell-Thomson算法 5294291.2正演程序以及成图 8168951.2.1勒夫波频散曲线正演流程图 8180521.2.2理论模型设计 956261.2.3试算程序 10150341.2.4试算结果 13216361.2.5成图分析 151.1勒夫波频散曲线正演弹性波矢量波动方程可以用以下公式表示：λ+μ式中：λ、μ为拉梅常数；ρ为密度；S为位移，将位移用势函数进行分解可获得：∇式中：φ、χ、ψ为质点的位移势，它们分别代表了P、SH、SV波型，其通解形式可写为：式中：a=kP2−k2，因为勒夫波是一种由SH波相互干涉而产生的正方向表面波，现仅考虑SH波的方程。假定勒夫波的介质为由n层均匀各向同性的弹性水平层状介质组成（参见图3-1），且勒夫波的频散曲线只与地层横波速度VS、层厚h及其密度ρ图3-1水平层状介质示意图1.1.1Knopoff算法考虑频率为ω、相速度为Vr由每一个界面的位移和应力分量连续性条件、自由边界性条件(其中的应力分量为零)、无限远处的辐射性条件，以及Laplace分解定理，根据与快速矢量传递算法类似的推导过程，可以获得以下形式的频散函数:F其中:当1≤m≤n−2时，TTT式中：Qm当m=n−1时且半空间的上界面为固-固界面时：T通过进一步化简，可以将勒夫波的频散函数写为：F其中：a1.1.2广义反射-透射系数算法图3-2二维水平层状介质模型考虑如图3-2所示的层状模型，其位移Uj(z)Vd式中：k为波数；ω为频率；λ、μ为拉梅常数；ρ为密度。其中：ζξ令fjd对于自由边界处有：P由连续条件有：f当i=N+1，z→∞时，f当每层内部，式3-8可进一步写为：f式中，Ei和∧i为已知矩阵；E∧式中γP=1−Vrγ=k其中fj令Cjj=1，2，⋯，N时，式3-12即瑞雷波的解可以写为：D类似的推导方式，可以推导得到勒夫波的解为：WW其中：∧根据连续性条件有：E整理得：T同理，可以得到：T引入广义反射-透射系数，其表达形式如下：C当j=1时：C对半空间而言：R在自由边界，由自由边界条件有：0=由此可得：R最终获得如下计算广义反射-透射系数的公式：T对于表层：C将两式合并得到：(1−由于有解存在，可得到勒夫波的频散方程：1−1.1.3Haskell-Thomson算法Haskell-Thomson算法[26]是Haskell于1953年在Thomson的数学理论基础研究成果基础上，通过对位于相邻两个流体界面之间的无限传递辐射矩阵边界公式、自由流体表面的无限边界辐射条件、以及位于无限远处的频散辐射边界条件进行推导发展出来的一种层状介质中关于平面勒夫波的无限频散辐射方程，其基本原理如下：对水平层状均匀介质模型的第i层有：&ϕτ上式中i表示层数（i=1，2，3，…，n，总层数为l，l=n+1）；x，z，t分别表示模型空间的水平坐标、垂直坐标和时间；ϕi和ψi为质点的位移势，ui和wi为质点的位移场，τix和τiz表示垂向应力场，上述参数均为x，z，t的函数；Ui，Wi，Xi，Zi表示位移应力，均为z的函数；Ai和Ai'分别表示P波的上下振幅，Bi和Bi'分别表示SV波的上下振幅；pi当c<Vr当c>Vr勒夫波的位移应力矢量和振幅为y上式Mi，PME边界条件：1.在自由表面处：X10=2.半空间的辐射条件为：当i=0时（上半空间）A0=B0=0我们定义一个改进增幅矢量ala对于c=V&对于矩阵Ei和P对于c=VPi，E对于c=VSi，E这些调整使得当c→V由(3-33)式可以得出yi=QiEi(−z)ai，Q是一个新的矩阵，Q我们定义层传递矩阵Ti(z)=QiEi(z)Qi−1，可以看出Ti也为一个4×4的矩阵，尽管QT如果z1和zy由上述两个矩阵方程结合消除ai变换可得yizy因此yi和yi代表当前层的顶部和底部的位移应力矢量，它们的关系为yi=Ty这就说明了第一层顶部位移应力向量和最后一层位移应力向量之间是由于所有层位移应力向量之间的乘积T=T1为了得到频散方程，我们只差对yl和yi应用表面条件和射线条件，在任何情况下，边界条件满足矩阵方程D上式中f为频率，f=ck/2π，m为地质模型，U'，V为自由表面条件决定的边界矩阵，对于勒夫波，U'为2×4的矩阵，U（注意：eij'左乘矩阵选择的是矩阵Q的第i行和第j行，而上述的理论框架是较为成熟的，在半空间找到标准的勒夫波频散方程比较简单。在没有地层介入，即n=0的情况下，频散方程可简化为D=det⁡下面我们简要的叙述实现频散函数Df，c，X这种算法需要子程序计算选取的相速度值c和波数k的模型中各层的传递矩阵Ti，中间矩阵Xi是一个2×4的矩阵，包含若干传递矩阵的乘积。对于层数较多的地质模型，传递矩阵的计算是很耗时的，为节省计算时间，我们可以有效地利用传递矩阵的元素的许多对称特性。在频散方程中，已知频率f和模型1.2正演程序以及成图1.2.1勒夫波频散曲线正演流程图将面波的频散曲线方程看作为一个函数:F(或F式中：fj为频率点；VRj和V把整个频散曲线的求取过程看作一个隐函数的求解过程，而该隐函数可以通过简单的二分法获得准确的结果。同时，该函数是一个多解的函数，即一个频点可能对应着多个零值点，即一个频率可能有多个相速度值。对于某一个给定的频率，我们将其最低速度称为基阶模式相速度(或第一模式速度)，比基阶相速度高的则依次成为一阶高阶模式相速度(或第二模式速度)、二阶高阶模式相速度(或第三模式速度)等。根据勒夫波的传播特性，一般采用介质中的最小横波速度的某个百分比(比如88%)作为相速度搜索的起始值，将最大横波速度作为相速度搜索的终止值；速度的搜索步长不宜过大，过大的搜索步长可能导致漏根，但过小的搜索步长可能导致计算量加大。为了方便解读求取频散曲线的过程，用程序流程图的形式(图3-3)对求取基阶频散曲线的整个流程进行描述。图3-3频散曲线正演程序流程示意图1.2.2理论模型设计为了进一步验证上述算法的有效性，本节建立速度递增型、含低速层夹层模型、含高速层夹层模型三层模型，各模型参数如表3-1至表3-6。其中表3-1、3-2、3-3中模型的地层单层厚度为10m，表3-1、3-2、3-3中的模型地层的单层厚度为2m。表3-1速度递增型地质模型参数值厚度（m）Vp（m/s）Vs（m/s）密度（g/m3）104002002.00106003002.00108004002.00表3-2含低速夹层型地质模型参数值厚度（m）Vp（m/s）Vs（m/s）密度（g/m3）106003002.00104002002.00108004002.00表3-3含高速夹层型地质模型参数值厚度（m）Vp（m/s）Vs（m/s）密度（g/m3）104002002.001010005002.00106003002.00表3-4薄层速度递增型地质模型参数值厚度（m）Vp（m/s）Vs（m/s）密度（g/m3）24002002.0026003002.0028004002.00表3-5薄层含低速夹层型地质模型参数值厚度（m）Vp（m/s）Vs（m/s）密度（g/m3）26003002.0024002002.0028004002.00表3-6薄层含高速夹层型地质模型参数值厚度（m）Vp（m/s）Vs（m/s）密度（g/m3）24002002.00210005002.0026003002.001.2.3试算程序此次编程采用Python语言，下载并安装专门用于面波正演程序包disba，disba是一个计算效率高的Python库，用于正演模拟面波频散。并且使用Thomson-Haskell方法计算勒夫波相/群速度频散曲线。以下为频散曲线正演程序：#!/usr/bin/envpython3#-*-coding:utf-8-*-#usage:ThefollowingexamplecomputestheRayleigh-andLove-wave#phasevelocitydispersioncurvesforthe3firstmodes#计算勒夫波各阶次的相速度频散曲线(根据输入模型计算理论频散曲线)importnumpyasnpfromdisbaimportPhaseDispersion#计算相速度频散#fromdisbaimportGroupDispersion#计算群速度频散importmatplotlibfrommatplotlibimportpyplotasplt#Velocitymodel#thickness，Vp，Vs，density#km，km/s，km/s，g/cm3#velocity_model=numpy.array##0.readvelocitymodel读取速度模型velocity_model=np.loadtxt('simple_mod.txt'，skiprows=1)#print('velocity_model='，velocity_model)#Periodsmustbesortedstartingwithlowperiods#t=numpy.logspace(0.0，1.0，100)t=np.logspace(0.0，1.0，100)#创建等比数列(周期)#print('t='，t)#Computethe3firstRayleigh-andLove-wavemodaldispersioncurves计R/L频散曲线#Fundamentalmodecorrespondstomode0(基阶mode=0)pd=PhaseDispersion(*velocity_model.T)#print(pd)#pd=GroupDispersion(*velocity_model.T)cpr=[pd(t，mode=i，wave="rayleigh")foriinrange(3)]cpl=[pd(t，mode=i，wave="love")foriinrange(3)]##开始画图fig，(ax1，ax2)=plt.subplots(1，2，figsize=(12，6)，sharey=False)##1.RayleighDispersionax1.set_xscale("log")ax1.set_xlim(1，1000)ax1.set_xlabel("Period[s]"，fontsize=15)ax1.set_ylim(1.2，4.8)ax1.grid(axis="y")ax1.set_ylabel("PhaseVelocity[km/s]"，fontsize=15)ax1.plot(cpr[0][0]，cpr[0][1]，color="blue"，linewidth=1，label="Fundamental")ax1.plot(cpr[1][0]，cpr[1][1]，color="orange"，linewidth=1，label="Mode1")ax1.plot(cpr[2][0]，cpr[2][1]，color="green"，linewidth=1，label="Mode2")ax1.legend(loc="lowerright"，fontsize=10)ax1.set_title("Rayleigh-wave"，fontsize=20)print("cpr="，cpr)print("mode0cpr="，cpr[0])print("mode1cpr="，cpr[1])print("mode2cpr="，cpr[2]) print("mode2cprperiod="，cpr[2][0]) #periodprint("mode2cprvelocity="，cpr[2][1]) #velocityprint("mode2cprmode="，cpr[2][2]) #modeprint("mode2cprwave="，cpr[2][3]) #waveprint("mode2cprtype="，cpr[2].type) #type#pdreturnsanamedtuple(period，velocity，mode，wave，type)##2.LoveDispersion#print("cpl="，cpl)ax2.set_xscale("log")ax2.set_xlim(1，1000)ax2.set_xlabel("Period[s]"，fontsize=15)ax2.set_ylim(1.2，4.8)ax2.grid(axis="y")ax2.set_ylabel("LoveVelocity[km/s]"，fontsize=15)ax2.plot(cpl[0].period，cpl[0].velocity，color="blue"，linewidth=1，label="Fundamental")ax2.plot(cpl[1].period，cpl[1].velocity，color="orange"，linewidth=1，label="Mode1")ax2.plot(cpl[2].period，cpl[2].velocity，color="green"，linewidth=1，label="Mode2")ax2.legend(loc="lowerright"，fontsize=10)ax2.set_title("Love-wave"，fontsize=20)##1.保存图片#plt.subplot_tool()#plt.show()#plt.tight_layout()plt.savefig('图片.png'，dpi=300)1.2.4试算结果通过程序试算模型可获得6个模型的频散曲线图，如图3-4至3-9。勒夫波频散曲线包含基阶模式和高阶模式。如图3-4所示，其中勒夫波基阶频散曲线为图中最外层曲线，其余为1阶至12阶频散曲线。基阶高阶基阶高阶图3-4速度递增型地质模型频散曲线图3-5含低速夹层地质模型频散曲线图3-6含高速夹层地质模型频散曲线图3-7薄层速度递