勒夫波頻散曲線計算編程實現(xiàn)及模型試算案例分析目錄TOC\o"1-3"\h\u25509勒夫波頻散曲線計算編程實現(xiàn)及模型試算案例分析 1177891.1勒夫波頻散曲線正演 121521.1.1Knopoff算法 2174051.1.2廣義反射-透射系數(shù)算法 328641.1.3Haskell-Thomson算法 5294291.2正演程序以及成圖 8168951.2.1勒夫波頻散曲線正演流程圖 8180521.2.2理論模型設(shè)計 956261.2.3試算程序 10150341.2.4試算結(jié)果 13216361.2.5成圖分析 151.1勒夫波頻散曲線正演彈性波矢量波動方程可以用以下公式表示：λ+μ式中：λ、μ為拉梅常數(shù)；ρ為密度；S為位移，將位移用勢函數(shù)進行分解可獲得：?式中：φ、χ、ψ為質(zhì)點的位移勢，它們分別代表了P、SH、SV波型，其通解形式可寫為：式中：a=kP2?k2，因為勒夫波是一種由SH波相互干涉而產(chǎn)生的正方向表面波，現(xiàn)僅考慮SH波的方程。假定勒夫波的介質(zhì)為由n層均勻各向同性的彈性水平層狀介質(zhì)組成（參見圖3-1），且勒夫波的頻散曲線只與地層橫波速度VS、層厚h及其密度ρ圖3-1水平層狀介質(zhì)示意圖1.1.1Knopoff算法考慮頻率為ω、相速度為Vr由每一個界面的位移和應(yīng)力分量連續(xù)性條件、自由邊界性條件(其中的應(yīng)力分量為零)、無限遠處的輻射性條件，以及Laplace分解定理，根據(jù)與快速矢量傳遞算法類似的推導過程，可以獲得以下形式的頻散函數(shù):F其中:當1≤m≤n?2時，TTT式中：Qm當m=n?1時且半空間的上界面為固-固界面時：T通過進一步化簡，可以將勒夫波的頻散函數(shù)寫為：F其中：a1.1.2廣義反射-透射系數(shù)算法圖3-2二維水平層狀介質(zhì)模型考慮如圖3-2所示的層狀模型，其位移Uj(z)Vd式中：k為波數(shù)；ω為頻率；λ、μ為拉梅常數(shù)；ρ為密度。其中：ζξ令fjd對于自由邊界處有：P由連續(xù)條件有：f當i=N+1，z→∞時，f當每層內(nèi)部，式3-8可進一步寫為：f式中，Ei和∧i為已知矩陣；E∧式中γP=1?Vrγ=k其中fj令Cjj=1，2，?，N時，式3-12即瑞雷波的解可以寫為：D類似的推導方式，可以推導得到勒夫波的解為：WW其中：∧根據(jù)連續(xù)性條件有：E整理得：T同理，可以得到：T引入廣義反射-透射系數(shù)，其表達形式如下：C當j=1時：C對半空間而言：R在自由邊界，由自由邊界條件有：0=由此可得：R最終獲得如下計算廣義反射-透射系數(shù)的公式：T對于表層：C將兩式合并得到：(1?由于有解存在，可得到勒夫波的頻散方程：1?1.1.3Haskell-Thomson算法Haskell-Thomson算法[26]是Haskell于1953年在Thomson的數(shù)學理論基礎(chǔ)研究成果基礎(chǔ)上，通過對位于相鄰兩個流體界面之間的無限傳遞輻射矩陣邊界公式、自由流體表面的無限邊界輻射條件、以及位于無限遠處的頻散輻射邊界條件進行推導發(fā)展出來的一種層狀介質(zhì)中關(guān)于平面勒夫波的無限頻散輻射方程，其基本原理如下：對水平層狀均勻介質(zhì)模型的第i層有：&?τ上式中i表示層數(shù)（i=1，2，3，…，n，總層數(shù)為l，l=n+1）；x，z，t分別表示模型空間的水平坐標、垂直坐標和時間；?i和ψi為質(zhì)點的位移勢，ui和wi為質(zhì)點的位移場，τix和τiz表示垂向應(yīng)力場，上述參數(shù)均為x，z，t的函數(shù)；Ui，Wi，Xi，Zi表示位移應(yīng)力，均為z的函數(shù)；Ai和Ai'分別表示P波的上下振幅，Bi和Bi'分別表示SV波的上下振幅；pi當c<Vr當c>Vr勒夫波的位移應(yīng)力矢量和振幅為y上式Mi，PME邊界條件：1.在自由表面處：X10=2.半空間的輻射條件為：當i=0時（上半空間）A0=B0=0我們定義一個改進增幅矢量ala對于c=V&對于矩陣Ei和P對于c=VPi，E對于c=VSi，E這些調(diào)整使得當c→V由(3-33)式可以得出yi=QiEi(?z)ai，Q是一個新的矩陣，Q我們定義層傳遞矩陣Ti(z)=QiEi(z)Qi?1，可以看出Ti也為一個4×4的矩陣，盡管QT如果z1和zy由上述兩個矩陣方程結(jié)合消除ai變換可得yizy因此yi和yi代表當前層的頂部和底部的位移應(yīng)力矢量，它們的關(guān)系為yi=Ty這就說明了第一層頂部位移應(yīng)力向量和最后一層位移應(yīng)力向量之間是由于所有層位移應(yīng)力向量之間的乘積T=T1為了得到頻散方程，我們只差對yl和yi應(yīng)用表面條件和射線條件，在任何情況下，邊界條件滿足矩陣方程D上式中f為頻率，f=ck/2π，m為地質(zhì)模型，U'，V為自由表面條件決定的邊界矩陣，對于勒夫波，U'為2×4的矩陣，U（注意：eij'左乘矩陣選擇的是矩陣Q的第i行和第j行，而上述的理論框架是較為成熟的，在半空間找到標準的勒夫波頻散方程比較簡單。在沒有地層介入，即n=0的情況下，頻散方程可簡化為D=det?下面我們簡要的敘述實現(xiàn)頻散函數(shù)Df，c，X這種算法需要子程序計算選取的相速度值c和波數(shù)k的模型中各層的傳遞矩陣Ti，中間矩陣Xi是一個2×4的矩陣，包含若干傳遞矩陣的乘積。對于層數(shù)較多的地質(zhì)模型，傳遞矩陣的計算是很耗時的，為節(jié)省計算時間，我們可以有效地利用傳遞矩陣的元素的許多對稱特性。在頻散方程中，已知頻率f和模型1.2正演程序以及成圖1.2.1勒夫波頻散曲線正演流程圖將面波的頻散曲線方程看作為一個函數(shù):F(或F式中：fj為頻率點；VRj和V把整個頻散曲線的求取過程看作一個隱函數(shù)的求解過程，而該隱函數(shù)可以通過簡單的二分法獲得準確的結(jié)果。同時，該函數(shù)是一個多解的函數(shù)，即一個頻點可能對應(yīng)著多個零值點，即一個頻率可能有多個相速度值。對于某一個給定的頻率，我們將其最低速度稱為基階模式相速度(或第一模式速度)，比基階相速度高的則依次成為一階高階模式相速度(或第二模式速度)、二階高階模式相速度(或第三模式速度)等。根據(jù)勒夫波的傳播特性，一般采用介質(zhì)中的最小橫波速度的某個百分比(比如88%)作為相速度搜索的起始值，將最大橫波速度作為相速度搜索的終止值；速度的搜索步長不宜過大，過大的搜索步長可能導致漏根，但過小的搜索步長可能導致計算量加大。為了方便解讀求取頻散曲線的過程，用程序流程圖的形式(圖3-3)對求取基階頻散曲線的整個流程進行描述。圖3-3頻散曲線正演程序流程示意圖1.2.2理論模型設(shè)計為了進一步驗證上述算法的有效性，本節(jié)建立速度遞增型、含低速層夾層模型、含高速層夾層模型三層模型，各模型參數(shù)如表3-1至表3-6。其中表3-1、3-2、3-3中模型的地層單層厚度為10m，表3-1、3-2、3-3中的模型地層的單層厚度為2m。表3-1速度遞增型地質(zhì)模型參數(shù)值厚度（m）Vp（m/s）Vs（m/s）密度（g/m3）104002002.00106003002.00108004002.00表3-2含低速夾層型地質(zhì)模型參數(shù)值厚度（m）Vp（m/s）Vs（m/s）密度（g/m3）106003002.00104002002.00108004002.00表3-3含高速夾層型地質(zhì)模型參數(shù)值厚度（m）Vp（m/s）Vs（m/s）密度（g/m3）104002002.001010005002.00106003002.00表3-4薄層速度遞增型地質(zhì)模型參數(shù)值厚度（m）Vp（m/s）Vs（m/s）密度（g/m3）24002002.0026003002.0028004002.00表3-5薄層含低速夾層型地質(zhì)模型參數(shù)值厚度（m）Vp（m/s）Vs（m/s）密度（g/m3）26003002.0024002002.0028004002.00表3-6薄層含高速夾層型地質(zhì)模型參數(shù)值厚度（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#計算勒夫波各階次的相速度頻散曲線(根據(jù)輸入模型計算理論頻散曲線)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)#創(chuàng)建等比數(shù)列(周期)#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試算結(jié)果通過程序試算模型可獲得6個模型的頻散曲線圖，如圖3-4至3-9。勒夫波頻散曲線包含基階模式和高階模式。如圖3-4所示，其中勒夫波基階頻散曲線為圖中最外層曲線，其余為1階至12階頻散曲線?；A高階基階高階圖3-4速度遞增型地質(zhì)模型頻散曲線圖3-5含低速夾層地質(zhì)模型頻散曲線圖3-6含高速夾層地質(zhì)模型頻散曲線圖3-7薄層速度遞