袁士川 宋先海*② 蔡 偉 胡 瑩 魯 鵬
(①中國(guó)地質(zhì)大學(xué)地球物理與空間信息學(xué)院,湖北武漢 430074; ②地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室,湖北武漢 430074)
·地震模擬·
瑞雷波有限差分?jǐn)?shù)值模擬中不同自由表面邊界條件的對(duì)比
袁士川①宋先海*①②蔡 偉①胡 瑩①魯 鵬①
(①中國(guó)地質(zhì)大學(xué)地球物理與空間信息學(xué)院,湖北武漢 430074; ②地球內(nèi)部多尺度成像湖北省重點(diǎn)實(shí)驗(yàn)室,湖北武漢 430074)
自由表面邊界條件是決定瑞雷波數(shù)值模擬效果的關(guān)鍵因素。本文基于標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格高階有限差分法,在二維各向同性彈性介質(zhì)背景下,針對(duì)應(yīng)力鏡像法(SIM)、改進(jìn)應(yīng)力鏡像法(MSIM)、橫向各向同性介質(zhì)替換法(MS)和聲學(xué)—彈性邊界近似法(AEA)等四種最具代表性的自由表面邊界條件進(jìn)行了數(shù)值模擬,并在均勻半空間模型中從波場(chǎng)快照、波形曲線和頻散曲線三個(gè)角度進(jìn)行了對(duì)比分析。在相同條件下,上述四種方法均能生成符合勘探地球物理規(guī)律的波場(chǎng)快照,各自對(duì)應(yīng)的數(shù)值解與解析解的擬合誤差都隨網(wǎng)格剖分精度的提高而減小,SIM和AEA數(shù)值模擬的穩(wěn)定性和精度都明顯高于MSIM和MS。 基于層狀介質(zhì)模型的進(jìn)一步研究表明: 對(duì)于簡(jiǎn)單模型,SIM和AEA都能得到比MSIM和MS更高精度的數(shù)值模擬結(jié)果; 對(duì)于復(fù)雜模型,AEA的精度最高,是最適合瑞雷波數(shù)值模擬的自由表面邊界條件。
瑞雷波 有限差分模擬 自由表面邊界條件 應(yīng)力鏡像法 聲學(xué)—彈性邊界近似法
在石油勘探、淺層反射及折射波人工地震勘探中,瑞雷波是一種強(qiáng)干擾波; 在天然地震中,瑞雷波是危害性最大的一種地震波。因此,在早期研究中,人們主要是根據(jù)瑞雷波的特點(diǎn),采取諸多方法消除其影響或減小其危害。但當(dāng)人們發(fā)現(xiàn)瑞雷波在層狀介質(zhì)中具有頻散特性后,就逐漸將其作為有效波并充分利用該特性研究各種地球物理問(wèn)題。所謂頻散特性,是指瑞雷波相速度隨頻率變化,而此特性主要與地層中縱、橫波速度、密度及厚度有關(guān),尤其對(duì)橫波速度和厚度敏感,因此有人提出了通過(guò)反演瑞雷波相速度獲取橫波速度和厚度的方法[1]。淺地表最關(guān)鍵的地震參數(shù)是橫波速度,故該方法在巖土分層、場(chǎng)地原位測(cè)試、工程質(zhì)量無(wú)損檢測(cè)和地震安全評(píng)價(jià)等淺地表結(jié)構(gòu)探測(cè)中得到廣泛應(yīng)用,其中最具代表性的是瑞雷波多道分析法(Multichannel Analysis of Surface Waves,MASW)[2]。
瑞雷波多道分析法通過(guò)反演瑞雷波相速度獲得淺地表的剪切波速度結(jié)構(gòu),而真實(shí)的淺地表地質(zhì)環(huán)境非常復(fù)雜,通過(guò)地震波場(chǎng)數(shù)值模擬技術(shù)研究典型的淺地表地質(zhì)模型,對(duì)提高實(shí)測(cè)資料采集質(zhì)量與反演解釋精度等都具有重要作用。數(shù)值模擬方法有很多種類,而有限差分法以其計(jì)算速度快、占用內(nèi)存小、能模擬復(fù)雜介質(zhì)等優(yōu)點(diǎn),在地震波場(chǎng)數(shù)值模擬中已得到廣泛應(yīng)用。例如: 董良國(guó)等[3]研究了一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法,并對(duì)其穩(wěn)定性進(jìn)行了分析; 周竹生等[4]利用有限差分算法實(shí)現(xiàn)了彈性介質(zhì)瑞雷波正演模擬; 熊章強(qiáng)等[5]基于瑞雷波數(shù)值模擬實(shí)例分析了差分階數(shù)和網(wǎng)格剖分精度對(duì)計(jì)算效率的影響; 邵廣周等[6]通過(guò)數(shù)值模擬對(duì)瑞雷波頻散曲線特征及各模式能量分布進(jìn)行了詳細(xì)分析; 李振春等[7]提出了變網(wǎng)格差分算法,在兼顧計(jì)算成本的同時(shí)提高了對(duì)復(fù)雜介質(zhì)的刻畫精度; 杜啟振等[8]基于優(yōu)化差分系數(shù)法實(shí)現(xiàn)了橫向各向同性介質(zhì)地震波場(chǎng)數(shù)值模擬。此外,考慮到實(shí)際地質(zhì)環(huán)境的復(fù)雜性,業(yè)界也對(duì)起伏地表、黏彈性和含孔隙性等也做了數(shù)值模擬研究[9-14]。由于瑞雷波是縱波(P波)和橫波(SV波)在自由表面處相長(zhǎng)干涉而形成的沿自由表面?zhèn)鞑サ奶厥獾卣鸩ǎ虼伺c體波相比,其波場(chǎng)數(shù)值模擬實(shí)現(xiàn)起來(lái)更困難,且對(duì)自由表面邊界條件有更高要求。如何處理自由表面問(wèn)題,將直接影響瑞雷波數(shù)值模擬效果。此前,已對(duì)適用于有限差分的自由表面邊界條件做了卓有成效的研究,這些成果可粗略地分為均勻介質(zhì)法和非均勻介質(zhì)法[15]兩大類。
均勻介質(zhì)自由表面處理方法是指基于波動(dòng)方程、通過(guò)在自由表面上直接設(shè)置應(yīng)力連續(xù)性條件而實(shí)現(xiàn)自由邊界處理的方法。Aki等[16]最早給出垂直自由表面應(yīng)力為零的自由表面邊界條件理論表達(dá)式,但該式在離散波動(dòng)方程求解過(guò)程中實(shí)現(xiàn)困難。為此,Levander[17]提出了應(yīng)用鏡像技術(shù)近似處理彈性自由表面的方法。隨后在對(duì)于自由表面以上虛擬層中速度分量的處理問(wèn)題上出現(xiàn)了分歧,出現(xiàn)了三種不同處理方式[17-19]。Robertsson[18]對(duì)比了這三種方式,指出在虛擬層中速度分量設(shè)置為零的鏡像處理方式最好,也就是現(xiàn)在使用最為廣泛的應(yīng)力鏡像法(Stress Image Method,SIM); 同時(shí),還指出Hestholm等[20]使用直接法[21]處理自由表面是不穩(wěn)定且精度較低的。王秀明等[22]基于常規(guī)應(yīng)力鏡像法,提出一種新的自由表面處理方法,稱之為改進(jìn)應(yīng)力鏡像法(Modified Stress Image Method,MSIM),該方法改變了自由表面在差分網(wǎng)格中的采樣位置,并使自由邊界條件的設(shè)置更加簡(jiǎn)單。
非均勻介質(zhì)自由表面處理方法是指在自由表面上,通過(guò)調(diào)整彈性參數(shù)分布間接地實(shí)現(xiàn)自由邊界處理的方法。這種思想最早由Boore[23]提出,之后經(jīng)過(guò)眾多研究與應(yīng)用,發(fā)展了幾種不同的處理方式。真空法較早被提出,它是通過(guò)在自由表面以上網(wǎng)格中設(shè)置拉梅常數(shù)為零和密度接近零的方法處理自由邊界,因其實(shí)現(xiàn)過(guò)程簡(jiǎn)單,故被廣泛使用[24,25]。Mittet[26]提出另一種被稱為“橫向各向同性介質(zhì)替換法”或“Mittet's Scheme(簡(jiǎn)稱為MS)”的非均勻介質(zhì)自由表面處理方法。與真空法相比,MS對(duì)拉梅常數(shù)和密度進(jìn)行了更細(xì)致的處理,使數(shù)值模擬的精度得到了提高。Xu等[27]考慮了自由表面上下橫向應(yīng)力保持連續(xù)的條件,提出了用聲學(xué)—彈性介質(zhì)邊界近似代替自由表面的處理方式,被稱作“聲學(xué)—彈性邊界近似法(Acoustic-elastic Boundary Approach,AEA)”。
不同自由表面處理方法的處理效果和精度存在差異。將多種方法做對(duì)比研究,對(duì)于瑞雷波和其他地震波場(chǎng)模擬研究都具有重要意義。Graves[25]通過(guò)數(shù)值模擬指出應(yīng)力鏡像法比真空法具有更高精度; Bolhen等[15]對(duì)比了應(yīng)力鏡像法與橫向各向同性介質(zhì)替換法后,認(rèn)為應(yīng)力鏡像法的精度更高; Xu等[27]將應(yīng)力鏡像法、橫向各向同性介質(zhì)替換法、真空法及聲學(xué)—彈性邊界近似法進(jìn)行了對(duì)比,指出聲學(xué)—彈性邊界近似法能達(dá)到與應(yīng)力鏡像法相當(dāng)?shù)木龋婵辗ň茸畹停?王周等[28]對(duì)直接法、應(yīng)力鏡像法、改進(jìn)應(yīng)力鏡像法、橫向各向同性介質(zhì)替換法及聲學(xué)—彈性邊界近似法進(jìn)行了對(duì)比研究,認(rèn)為橫向各向同性介質(zhì)替換法的精度和可靠性最高。可見(jiàn): 一方面現(xiàn)有的對(duì)比研究還欠充分、全面; 另一方面各自結(jié)論不統(tǒng)一,且基本上只從均勻半空間波形曲線角度進(jìn)行對(duì)比。然而,瑞雷波波場(chǎng)數(shù)值模擬不僅具有波形曲線解析解,還存在頻散曲線解析解。波形曲線只利用了單道地震記錄信息; 而頻散曲線則綜合考慮了多道地震記錄信息,消除了瑞雷波近場(chǎng)效應(yīng)和遠(yuǎn)場(chǎng)效應(yīng)影響。“頻散”是瑞雷波最重要的特性,故從頻散曲線角度進(jìn)行對(duì)比很有必要。此外,從波場(chǎng)快照可直觀地分析瑞雷波傳播特性,故波場(chǎng)快照也是可選的對(duì)比“角度”。因此,從多方面綜合考慮,全面對(duì)比、評(píng)估各種自由表面邊界條件的數(shù)值模擬精度和穩(wěn)定性,才可得出更可靠結(jié)論。
鑒于此,本文在上述研究基礎(chǔ)上,利用標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格高階有限差分算法,針對(duì)二維各向同性彈性介質(zhì),在均勻半空間、兩層速度遞增、三層含低速軟夾層及四層含高速硬夾層等四種典型地質(zhì)模型中,從波場(chǎng)快照、波形曲線和頻散曲線三個(gè)“角度”對(duì)(應(yīng)力鏡像法、改進(jìn)應(yīng)力鏡像法、橫向各向同性介質(zhì)替換法、聲學(xué)—彈性邊界近似法等)四種最具代表性的自由表面處理方法,進(jìn)行了定性—定量、綜合的數(shù)值模擬對(duì)比研究,以期為瑞雷波和其他地震波場(chǎng)模擬研究提供科學(xué)和有價(jià)值的參考。
本文采用各向同性介質(zhì)一階速度—應(yīng)力彈性波方程進(jìn)行有限差分?jǐn)?shù)值模擬。對(duì)于二維P-SV波 (xoz平面),其差分形式可表示為[27]
(1)

(2)
自由表面是固體地球與空氣層之間的物性突變面,瑞雷波正是基于它的存在才得以形成和傳播。Aki等[16]最早給出垂直自由表面應(yīng)力為零時(shí)的自由表面邊界條件,二維水平自由表面(z=0)處的表達(dá)式為
(3)
該式在理論上非常簡(jiǎn)單,而在離散波動(dòng)方程求解過(guò)程中卻難以實(shí)現(xiàn),且空間差分階數(shù)越高,處理越困難。但重要的是,它為后續(xù)各種自由表面處理方法,提供了思想基礎(chǔ)。下面將詳細(xì)介紹用于對(duì)比的四種自由表面邊界條件的基本原理及其設(shè)置方法。
應(yīng)力鏡像法[18]的基本思想是把自由表面看作一面鏡子,在鏡面之上設(shè)置虛擬層,并使上、下兩側(cè)與垂向相關(guān)的應(yīng)力(σzz、σxz)關(guān)于自由表面鏡像對(duì)稱,虛擬層中速度分量設(shè)置為零。在二維交錯(cuò)網(wǎng)格中,SIM認(rèn)為自由表面通過(guò)正應(yīng)力σxx、σzz和速度水平分量vx的采樣位置,即自由表面位于整數(shù)網(wǎng)格線(粗線)上。如圖1所示,i和j分別為橫向和縱向網(wǎng)格剖分點(diǎn),j也表示自由表面所在位置,藍(lán)色正方形表示σxx和σzz,藍(lán)色圓形表示σxz,綠色倒三角形表示vz,綠色正三角形表示vx。

圖1 應(yīng)力鏡像法自由表面取樣位置
在2L階標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格有限差分條件下,虛擬層寬度為L(zhǎng)個(gè)網(wǎng)格點(diǎn),則應(yīng)力設(shè)置方式為
(4)
式中k(=1,2,…,L)為虛擬層中正參與鏡像處理的網(wǎng)格點(diǎn)距離自由表面的網(wǎng)格點(diǎn)數(shù)。式(4)中正應(yīng)力和切應(yīng)力在邊界上都是奇函數(shù),保證它們?cè)谶吔缟蠟榱恪?/p>
同時(shí),在自由表面上,正應(yīng)力σxx采用下式(據(jù)式(3))進(jìn)行迭代
(5)
利用該式進(jìn)行迭代可使數(shù)值模擬更加穩(wěn)健。
改進(jìn)應(yīng)力鏡像法[22](MSIM)與SIM的不同之處在于,在二維交錯(cuò)網(wǎng)格中,MSIM認(rèn)為自由表面是通過(guò)切應(yīng)力σxz和速度垂直分量vz的采樣位置,即自由表面位于半網(wǎng)格線上(圖2中細(xì)線)。

圖2 改進(jìn)應(yīng)力鏡像法自由表面取樣位置
在2L階標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格有限差分條件下,虛擬層具有L個(gè)網(wǎng)格點(diǎn)寬度,其應(yīng)力設(shè)置方式為
(6)
與SIM相比,MSIM具有兩大優(yōu)點(diǎn)[22]: ①由于正應(yīng)力σxx、σzz不在自由表面上采樣,因此不需單獨(dú)對(duì)σxx做式(5)的處理,直接采用式(6)就能滿足條件; ②當(dāng)邊界形狀構(gòu)成一個(gè)頂角時(shí),SIM需正應(yīng)力σxx、σzz為零,而MSIM只要求切應(yīng)力為零,因此效率有所提高,在物理意義上也是合理的。
橫向各向同性介質(zhì)替換法[26](MS)的基本思想是利用橫向各向同性介質(zhì)近似代替自由表面,不直接對(duì)式(3)的所有應(yīng)力條件進(jìn)行處理,而僅令正應(yīng)力σzz在自由表面處直接為零,切應(yīng)力σxz為零的條件通過(guò)對(duì)自由表面上物性參數(shù)的設(shè)定,在波動(dòng)方程交錯(cuò)網(wǎng)格有限差分迭代求解過(guò)程中實(shí)現(xiàn)。二維情況下MS的處理方式為
(7)
式中:σzz、ρx和λ、μ對(duì)應(yīng)表示自由表面上的正應(yīng)力、與vx相關(guān)的密度和拉梅常數(shù);ρ1和μ1分別表示自由表面以下彈性介質(zhì)的密度和剪切模量。圖3展示各分量和參數(shù)在交錯(cuò)網(wǎng)格中分布情況。

圖3 橫向各向同性介質(zhì)替換法自由表面取樣位置
MS相對(duì)于鏡像法而言,處理方式更簡(jiǎn)單,更易于編程實(shí)現(xiàn),計(jì)算效率也更高。MS不需對(duì)虛擬層中的應(yīng)力做鏡像處理,只需在自由表面上對(duì)正應(yīng)力σzz和部分彈性參數(shù)進(jìn)行處理,就能滿足式(3)自由表面邊界條件的要求。
聲學(xué)—彈性邊界近似法(AEA)[27]考慮了自由表面上、下橫向應(yīng)力保持連續(xù)的條件,其他與MS思路大體相同。在二維情況下AEA的處理方式為
(8)
圖4展示了自由表面在交錯(cuò)網(wǎng)格中的取位位置及各分量和參數(shù)在交錯(cuò)網(wǎng)格中的分布情況。

圖4 聲學(xué)—彈性邊界近似法自由表面取樣位置
不同于MS,AEA主要考慮自由表面上、下橫向應(yīng)力保持連續(xù)的條件,具體在式(8)中表現(xiàn)為2μ=2μ1,即拉梅常數(shù)μ不變。Xu等[27]指出MS違背了部分基本物理原理,而AEA較符合物理規(guī)律,且模擬精度更高。
本文選取標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格高階有限差分算法做瑞雷波波場(chǎng)數(shù)值模擬: 采用時(shí)間2階、空間12階差分格式;力源為自由表面激發(fā)的z方向點(diǎn)力源,震源函數(shù)為主頻20Hz、延遲時(shí)間50ms的高斯一階導(dǎo)數(shù);采集道數(shù)為60,最小炮檢距為20m,道間距為1m;采用完美匹配層(Perfectly Matched Layer,PML)作為吸收邊界條件,吸收寬度為20m,理論反射系數(shù)為0.000001(下同)。利用瑞雷波最小波長(zhǎng)λmin與空間步長(zhǎng)dh之比(points per minimum wavelength,ppw)度量網(wǎng)格剖分精度[15],用Cagniard-De Hoop技術(shù)求解均勻半空間波場(chǎng)解析解[29],以快速矢量傳遞算法正演理論頻散曲線[30],通過(guò)高分辨率線性拉東變換(Linear Radon Transform,LRT)提取頻散能量圖[31]。
本文首先設(shè)計(jì)均勻半空間模型(模型A),其地質(zhì)模型參數(shù)見(jiàn)表1。針對(duì)該模型分別從波場(chǎng)快照、波形曲線及頻散曲線三個(gè)角度對(duì)自由表面邊界條件進(jìn)行定性—定量對(duì)比分析。

表1 均勻半空間模型A參數(shù)
波場(chǎng)快照包含豐富波場(chǎng)信息,包括波的種類及其傳播特性。通過(guò)波場(chǎng)快照,不僅可直觀地再現(xiàn)瑞雷波在不同時(shí)刻的傳播特性,而且可分析整個(gè)地震波場(chǎng)在不同時(shí)刻縱波、橫波、轉(zhuǎn)換波等與瑞雷波的相互關(guān)系、刻畫各種波的能量分配方式,進(jìn)而從地震波場(chǎng)傳播角度評(píng)價(jià)四種自由表面邊界條件用于瑞雷波波場(chǎng)數(shù)值模擬的可行性和優(yōu)越性。設(shè)置模擬區(qū)域?yàn)?00m×50m,時(shí)間步長(zhǎng)dt=0.1ms,空間步長(zhǎng)dh=0.25m,分別采用四種自由表面邊界條件對(duì)模型A進(jìn)行模擬,得到時(shí)間t=0.18s的波場(chǎng)快照。圖5分別為SIM、MSIM、MS及AEA對(duì)應(yīng)的vz和vx分量的波場(chǎng)快照。
由圖5可知: 瑞雷波(RW)在自由表面生成且沿自由表面?zhèn)鞑ィ且环N面波; 橫波(S)以震源為中心向外傳播,且充滿整個(gè)空間,是一種體波; 縱波速度遠(yuǎn)高于橫波速度,在該時(shí)刻縱波已被PML完全吸收; 瑞雷波速度稍慢于橫波,但其能量在整個(gè)波場(chǎng)中占據(jù)主導(dǎo)地位。通過(guò)對(duì)比可知,四種方法均能生成符合勘探地球物理規(guī)律的波場(chǎng)快照,均適合作為瑞雷波波場(chǎng)數(shù)值模擬的自由表面邊界條件。
瑞雷波的模擬精度與網(wǎng)格剖分精度密切相關(guān),為對(duì)比不同網(wǎng)格剖分精度下四種自由表面邊界條件的數(shù)值模擬效果,設(shè)計(jì)了瑞雷波最小波長(zhǎng)對(duì)應(yīng)網(wǎng)格點(diǎn)數(shù)為10、20和40ppw的三組網(wǎng)格剖分參數(shù)(表2)。通過(guò)對(duì)模型A做模擬,得到對(duì)應(yīng)炮集記錄。圖6是以AEA和20ppw網(wǎng)格剖分精度為例,模擬得到的均勻半空間vz和vx分量炮集記錄。分別抽取炮集記錄上炮檢距為49m(第30道)處的單道地震記錄,并與解析解進(jìn)行對(duì)比,如圖7所示。

表2 均勻半空間模型網(wǎng)格剖分參數(shù)
從圖6可知,因采用標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格高階有限差分?jǐn)?shù)值模擬技術(shù)和PML吸收邊界條件,大大削弱了數(shù)值頻散和人工邊界反射等數(shù)值模擬假象。炮集記錄上既有能量占主導(dǎo)地位的瑞雷波,也有能量較弱的體波: 瑞雷波清晰可見(jiàn),同相軸為一條傾斜直線,在所有地震道上均具有很強(qiáng)能量; 但體波僅在近道上可見(jiàn)明顯波形,且其能量隨著炮檢距的增大而迅速減小。
為準(zhǔn)確表征有限差分?jǐn)?shù)值解與解析解的相對(duì)誤差,應(yīng)用L2范數(shù)誤差量化不同情況下的模擬誤差[15]

(9)
式中:f(lΔt)代表數(shù)值解;q(lΔt)代表解析解;MΔt和NΔt分別為瑞雷波的起跳時(shí)間和結(jié)束時(shí)間。顯然,E越小,數(shù)值模擬精度則越高。利用式(9)可計(jì)算圖7中瑞雷波波形曲線數(shù)值解與解析解的L2范數(shù)誤差,計(jì)算結(jié)果見(jiàn)圖8的誤差變化曲線所示。

圖5 均勻半空間模型t=0.18s時(shí)刻質(zhì)點(diǎn)速度場(chǎng)快照
圖6 均勻半空間模型的炮集記錄
(a)AEA,vz分量; (b)AEA,vx分量

圖7 四種自由邊界條件在三種網(wǎng)格剖分精度下瑞雷波波形曲線數(shù)值解與解析解的對(duì)比
由圖7可看出: SIM的波形擬合效果最好,其不同網(wǎng)格剖分精度對(duì)應(yīng)的數(shù)值解與解析解都基本重合(圖7a和圖7b); AEA的模擬精度僅次于SIM,其瑞雷波只存在輕微相移(圖7g和圖7h); MSIM(圖7c和圖7d)和MS(圖7e和圖7f)的數(shù)值解不僅具有明顯相移,也伴隨明顯的振幅改變,且網(wǎng)格剖分精度越低,該現(xiàn)象越明顯。由圖8可知,四種方法對(duì)應(yīng)的L2范數(shù)相對(duì)誤差都隨網(wǎng)格剖分精度的提高而減小,SIM和AEA的相對(duì)誤差變化較平緩,MSIM和MS的相對(duì)誤差變化較大,說(shuō)明SIM和AEA的穩(wěn)定性更好。在粗網(wǎng)格剖分(10ppw)時(shí),MS和MSIM的擬合相對(duì)誤差極大,MS對(duì)應(yīng)的vz和vx分量的相對(duì)誤差分別高達(dá)63.90%和110.41%,MSIM對(duì)應(yīng)的相對(duì)誤差分別高達(dá)48.45%和83.42%,而AEA和SIM的擬合相對(duì)誤差卻較小,AEA對(duì)應(yīng)的相對(duì)誤差分別為9.94%和33.45%,SIM對(duì)應(yīng)的相對(duì)誤差最低,分別為1.28%和9.85%; 在20ppw網(wǎng)格剖分情況下,MSIM和MS的擬合相對(duì)誤差雖有明顯減小,但與SIM和AEA的模擬精度仍有差距;當(dāng)網(wǎng)格剖分精度為40ppw時(shí),四種方法的誤差較接近,但相對(duì)誤差由低到高的順序(SIM↗AEA↗MSIM↗MS)保持不變(圖8)。 通過(guò)定性—定量分析都可看出,在三種網(wǎng)格剖分精度下,SIM和AEA的數(shù)值模擬精度和穩(wěn)定性都明顯高于MSIM和MS,即使當(dāng)網(wǎng)格剖分精度很低時(shí),SIM和AEA仍能保持較高的模擬精度,而MSIM和MS卻不能。

圖8 四種自由邊界條件對(duì)應(yīng)的L2范數(shù)相對(duì)誤差隨網(wǎng)格剖分精度的變化情況
瑞雷波不僅具有波形曲線解析解,還存在頻散曲線解析解。波形曲線對(duì)比只利用了單道地震記錄信息,所選擇的目標(biāo)地震道可能受瑞雷波近場(chǎng)效應(yīng)和遠(yuǎn)場(chǎng)效應(yīng)的影響,因此僅從該角度進(jìn)行不同自由表面邊界條件對(duì)比研究,得出的結(jié)論可能會(huì)有片面性。而頻散曲線是瑞雷波最重要的特性,頻散曲線綜合考慮了多道地震記錄信息,消除了瑞雷波近場(chǎng)效應(yīng)和遠(yuǎn)場(chǎng)效應(yīng)的影響,故有必要從頻散曲線角度進(jìn)行對(duì)比。二維瑞雷波數(shù)值模擬可采集到vz和vx分量的炮集記錄,兩個(gè)分量都包含模型信息,充分利用兩個(gè)分量的信息,能更好地反映模型的真實(shí)情況。vz和vx炮集記錄通過(guò)高分辨率LRT可分別得到vz和vx分量頻散能量數(shù)據(jù),再通過(guò)加權(quán)平均得到二分量頻散能量數(shù)據(jù)。本文均采用該類型數(shù)據(jù)進(jìn)行對(duì)比研究,權(quán)值均取為0.5。這樣綜合考慮了自由表面邊界條件對(duì)兩個(gè)分量模擬效果的影響。
為此,分別對(duì)圖7中對(duì)應(yīng)的炮集記錄進(jìn)行高分辨率LRT計(jì)算,得到相應(yīng)的頻散能量數(shù)據(jù)。圖9a、圖9c、圖9e和圖9g是以20ppw網(wǎng)格剖分精度,分別采用SIM、MSIM、MS、AEA得到的頻散能量圖。頻散能量圖上能量極大值點(diǎn)組成的曲線即為瑞雷波頻散曲線數(shù)值解,分別從此四頻散能量圖中提取頻散曲線數(shù)值解,并與解析解進(jìn)行對(duì)比,得到圖9b、圖9d、圖9f和圖9h,而數(shù)值解與解析解的L2范數(shù)相對(duì)誤差被繪制成誤差變化曲線(圖10)。
瑞雷波在均勻半空間不發(fā)生頻散,其相速度為一定值,頻散曲線解析解為一條水平直線。模擬得到的數(shù)值解只是近似解,并不能與解析解完全一致,其擬合誤差會(huì)隨著模擬精度的提高而減小。由圖9a、圖9c、圖9e和圖9g可直觀看出: 在同一網(wǎng)格剖分精度(20ppw)下,AEA模擬效果最好,解析解恰好從其頻散能量極值區(qū)域中心處穿過(guò)(圖9g); 其次是SIM,它對(duì)應(yīng)的頻散能量在高頻處更收斂,但有輕微上揚(yáng)(圖9a); 而MSIM和MS在高頻段存在不同程度的下跌(圖9c和圖9e)。
由圖9b、圖9d、圖9f和圖9h可知: 三種網(wǎng)格剖分精度下,AEA對(duì)應(yīng)的頻散曲線數(shù)值解與解析解偏差都很小(圖9h); 在粗網(wǎng)格剖分(10ppw)時(shí),SIM(圖9b)、MSIM(圖9d)及MS(圖9f)對(duì)應(yīng)的數(shù)值解與解析解都存在不同程度的較大偏差,對(duì)應(yīng)相對(duì)誤差分別為8.1684E-05、3.3763E-04和1.5812E-04,而AEA的誤差卻僅有2.2193E-05(圖10)。由圖10可看出: 隨著網(wǎng)格剖分精度的提高,四種方法對(duì)應(yīng)的誤差都在減小; AEA的誤差變化近似一條水平直線,SIM只是在10ppw時(shí)誤差比AEA高,在20ppw和40ppw時(shí),SIM和AEA數(shù)值模擬精度相當(dāng); 從整體上看,MSIM和MS的誤差都高于SIM和AEA; 在40ppw時(shí),四種方法的誤差均很小,但SIM (1.5246E-05)和AEA(1.9670E-05)的數(shù)值模擬精度仍高于MSIM(2.3724E-05)和MS(2.9539E-05)。通過(guò)定性—定量分析易知,在三種網(wǎng)格剖分精度下,SIM和AEA的模擬精度和穩(wěn)定性都明顯高于MSIM和MS,即使在粗網(wǎng)格剖分(10ppw)時(shí),SIM和AEA仍能保持較高的模擬精度,而MSIM和MS卻不盡人意。

圖9 四種自由邊界條件對(duì)應(yīng)的頻散能量圖和頻散曲線數(shù)值解與解析解的對(duì)比圖


圖10 四種自由邊界條件對(duì)應(yīng)的L2范數(shù)相對(duì)誤差隨網(wǎng)格剖分精度的變化
通過(guò)均勻半空間模型中定性—定量對(duì)比可知,在瑞雷波波場(chǎng)數(shù)值模擬中,SIM和AEA的模擬精度和穩(wěn)定性都明顯高于MSIM和MS。為驗(yàn)證該結(jié)論,并進(jìn)一步對(duì)比SIM與AEA的優(yōu)劣,又設(shè)計(jì)了三種典型的層狀介質(zhì)模型對(duì)四種自由表面邊界條件的模擬效果進(jìn)行對(duì)比。
設(shè)計(jì)兩層速度遞增模型(模型B),其地質(zhì)模型參數(shù)如表3所示。設(shè)置時(shí)間步長(zhǎng)dt=0.05ms,空間步長(zhǎng)dh=0.25m,采集時(shí)長(zhǎng)為0.9s,分別采用四種自由表面邊界條件進(jìn)行數(shù)值模擬,可得到對(duì)應(yīng)的炮集記錄。圖11是以AEA為例,模擬得到的vz和vx分量炮集記錄。炮集記錄通過(guò)高分辨率LRT可得到對(duì)應(yīng)的頻散能量圖(圖12)。
由圖11可看出,因采用標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格高階有限差分?jǐn)?shù)值模擬技術(shù)和PML吸收邊界條件,大大削弱了數(shù)值頻散和人工邊界反射等數(shù)值模擬假象。圖中瑞雷波清晰可見(jiàn),且呈發(fā)散的掃帚狀。此外,還有直達(dá)波、折射波、反射波等體波,但瑞雷波具有更高的信噪比,其能量在整個(gè)地震記錄中占主導(dǎo)地位。在圖12頻散能量圖上,既有能量占主導(dǎo)地位的基階模式,也有多條高階模式發(fā)育,高階模式隨著階數(shù)的增加,能量逐漸減弱; 對(duì)于數(shù)值解與解析解的擬合程度,基階模式擬合得最好,高階模式隨著階數(shù)的增加,擬合度逐漸減小。四種自由表面邊界條件的模擬結(jié)果中,AEA(圖12d)的模擬效果最好,高階模式的能量更為連續(xù),且擬合度都較其他三種方法高,其次是SIM(圖12a)和MS(圖12c),MSIM(圖12b)的擬合情況相對(duì)最差。該模擬是在20ppw網(wǎng)格剖分精度情況下進(jìn)行的,如果進(jìn)一步提高網(wǎng)格剖分精度,四種方法的模擬效果會(huì)更好。

表3 兩層速度遞增模型B的參數(shù)

圖11 兩層速度遞增模型的炮集記錄

圖12 四種自由表面邊界條件對(duì)應(yīng)的兩層速度遞增模型的頻散能量圖
設(shè)計(jì)三層含低速軟夾層模型(模型C),其地質(zhì)模型參數(shù)如表4所示。設(shè)置時(shí)間步長(zhǎng)dt=0.05ms,空間步長(zhǎng)dh=0.25m,采集時(shí)長(zhǎng)為0.9s,通過(guò)數(shù)值模擬可得到對(duì)應(yīng)的炮集記錄。圖13是以AEA為例模擬得到的vz和vx分量炮集記錄。炮集記錄通過(guò)高分辨率LRT可得到對(duì)應(yīng)的頻散能量圖(圖14)。
在圖14頻散能量圖上既有基階模式,也有高階模式,基階模式與高階模式能量都很強(qiáng),且均有截止頻率,如基階模式截止頻率約為20Hz,第一高階模式截止頻率約為28Hz。由于低速夾層的存在,高階模式存在嚴(yán)重的“模式接吻”現(xiàn)象[32](圖14a紅圈圈出部分)。在5~20Hz頻段內(nèi),四種方法對(duì)應(yīng)的基階模式與解析解均擬合得很好; 在20~40Hz頻段內(nèi),四種方法對(duì)應(yīng)的第一和第二高階模式與解析解都存在一定程度的偏差; 在40~90Hz頻段內(nèi),AEA的數(shù)值解與解析解擬合情況最好,解析解正好從頻散能量圖極大值區(qū)域的中心穿過(guò),其次是SIM,但是其能量圖分辨率和連續(xù)性較差,MSIM和MS的擬合情況相當(dāng),均與解析解有一定程度的偏差。

表4 三層含低速軟夾層模型C的參數(shù)

圖13 三層含低速軟夾層模型的炮集記錄

圖14 四種自由表面邊界條件對(duì)應(yīng)的三層含低速軟夾層模型的頻散能量圖
設(shè)計(jì)四層含高速硬夾層模型(模型D),其地質(zhì)模型參數(shù)如表5所示。設(shè)置時(shí)間步長(zhǎng)dt=0.04ms,空間步長(zhǎng)dh=0.125m, 采集時(shí)長(zhǎng)為0.8s, 經(jīng)數(shù)值模擬可得到對(duì)應(yīng)的炮集記錄。圖15是以AEA為例模擬得到的vz和vx分量炮集記錄。炮集記錄通過(guò)高分辨率LRT可得到對(duì)應(yīng)的頻散能量圖(圖16)。
該模型既是高泊松比又含高速硬夾層,數(shù)值模擬難度有所增加,故采用40ppw網(wǎng)格剖分精度進(jìn)行模擬。由圖16可知,四種方法對(duì)應(yīng)的頻散能量圖上,基階模式在所有頻段內(nèi)能量都很強(qiáng),第一高階模式僅在10~20Hz頻段內(nèi)能量較強(qiáng),其他高階模式能量都很弱,甚至未發(fā)育。在該精度條件下,四種方法對(duì)應(yīng)的數(shù)值模擬效果相當(dāng),僅在10~30Hz頻段內(nèi)可發(fā)現(xiàn)較小差異,即此頻段內(nèi)AEA的模擬精度最高,其基階和第一高階的頻散能量最大值比另三種方法都更接近于解析解(圖16d)。

表5 四層含高速硬夾層模型D的參數(shù)

圖15 四層含高速硬夾層模型的炮集記錄

圖16 四種自由表面邊界條件對(duì)應(yīng)的四層含高速硬夾層模型的頻散能量圖
本文基于標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格高階有限差分算法,分別在均勻半空間、兩層速度遞增、三層含低速軟夾層和四層含高速硬夾層模型中,對(duì)瑞雷波波場(chǎng)數(shù)值模擬中最具代表性的四種自由表面邊界條件:應(yīng)力鏡像法(SIM)、改進(jìn)應(yīng)力鏡像法(MSIM)、橫向各向同性介質(zhì)替換法(MS)、聲學(xué)—彈性邊界近似法(AEA),進(jìn)行了數(shù)值模擬對(duì)比研究。結(jié)果表明:
(1)在均勻半空間模型中,四種方法均能生成符合勘探地球物理規(guī)律的波場(chǎng)快照;四種方法對(duì)應(yīng)的瑞雷波波形曲線和頻散曲線的數(shù)值解與解析解的擬合誤差都隨網(wǎng)格剖分精度的提高而減小。在相同網(wǎng)格剖分情況下,SIM和AEA的模擬精度和穩(wěn)定性都明顯高于MSIM和MS。即使在粗網(wǎng)格剖分時(shí),SIM和AEA仍有很高模擬精度,而MSIM和MS對(duì)應(yīng)的相對(duì)誤差卻較大。
(2)層狀介質(zhì)模型的模擬結(jié)果,驗(yàn)證了均勻半空間模型中的結(jié)論,并進(jìn)一步對(duì)比了SIM和AEA的模擬精度。結(jié)果表明,在復(fù)雜層狀介質(zhì)模型中,AEA的模擬精度一定程度上優(yōu)于SIM。
綜上所述,在瑞雷波有限差分?jǐn)?shù)值模擬中,對(duì)于簡(jiǎn)單模型,SIM和AEA都能得到比MSIM和MS更高精度的模擬結(jié)果; 對(duì)于復(fù)雜地質(zhì)模型,AEA的模擬效果最好、精度最高,是最適合瑞雷波數(shù)值模擬的自由表面邊界條件。
本文是在各向同性彈性介質(zhì)條件下展開研究,雖然所設(shè)計(jì)模型考慮了實(shí)際淺地表高泊松比的特點(diǎn),但實(shí)際淺地表環(huán)境更復(fù)雜,既具有各向異性,又具有黏彈性。那么,對(duì)于各向異性黏彈性介質(zhì)瑞雷波數(shù)值模擬,其自由表面邊界條件該如何設(shè)置?這有待后續(xù)深入研究、探討。
感謝中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院汪利民老師分享計(jì)算彈性均勻半空間波場(chǎng)解析解的程序。
[1] Xia J H,Miller R D,Park C B.Estimation of near-surface shear-wave velocity by inversion of Rayleigh waves.Geophysics,1999,64(3):691-700.
[2] Park C B,Miller R D,Xia J H.Multichannel analysis of surface waves (MASW).Geophysics,1999,64(3):800-808.
[3] 董良國(guó),馬在田,曹景忠.一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法穩(wěn)定性研究.地球物理學(xué)報(bào),2000,43(6):856-864.
Dong Liangguo,Ma Zaitian,Cao Jingzhong.A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation.Chinese Journal of Geophysics,2000,43(6):856-864.
[4] 周竹生,劉喜亮,熊孝雨.彈性介質(zhì)中瑞雷面波有限差分法正演模擬.地球物理學(xué)報(bào),2007,50(2):567-573.
Zhou Zhusheng,Liu Xiliang,Xiong Xiaoyu.Finite difference modelling of Rayleigh surface wave in elastic media.Chinese Journal of Geophysics,2007,50(2):567-573.
[5] 熊章強(qiáng),張大洲,秦臻等.瑞雷波數(shù)值模擬中的邊界條件及模擬實(shí)例分析.中南大學(xué)學(xué)報(bào)(自然科學(xué)版),2008,39(4):824-830.
Xiong Zhangqiang,Zhang Dazhou,Qin Zheng et al.Boundary conditions and case analysis of numerical modeling of Rayleigh wave.Journal of Central South University (Science and Technology Edition),2008,39(4):824-830.
[6] 邵廣周,李慶春,吳華.基于波場(chǎng)數(shù)值模擬的瑞利波頻散曲線特征及各模式能量分布.石油地球物理勘探,2015,50(2):306-315.
Shao Guangzhou,Li Qingchun,Wu Hua.Dispersion curves and mode energy distribution of Rayleigh wave based on wavefield numerical simulation.OGP,2015,50(2):306-315.
[7] 李振春,張慧,張華.一階彈性波方程的變網(wǎng)格高階有限差分?jǐn)?shù)值模擬.石油地球物理勘探,2008,43(6):711-716.
Li Zhenchun,Zhang Hui,Zhang Hua.Variable-grid high-order finite-difference numeric simulation of first-order elastic wave equation.OGP,2008,43(6):711-716.
[8] 杜啟振,白清云,李賓.橫向各向同性介質(zhì)優(yōu)化差分系數(shù)法地震波場(chǎng)數(shù)值模擬.石油地球物理勘探,2010,45(2):170-176.
Du Qizhen,Bai Qingyun,Li Bin.Optimized difference coefficient method seismic wave field numerical simulation in vertical transverse isotropic medium.OGP,2010,45(2):170-176.
[9] Wang L M,Luo Y H,Xu Y X.Numerical investigation of Rayleigh-wave propagation on topography surface.Journal of Applied Geophysics,2012,86(11):88-97.
[10] 高靜懷,何洋洋,馬逸塵.黏彈性與彈性介質(zhì)中Rayleigh面波特性對(duì)比研究.地球物理學(xué)報(bào),2012,55(1):207-218.
Gao Jinghuai,He Yangyang,Ma Yichen.Comparison of the Rayleigh wave in elastic and viscoelastic media.Chinese Journal of Geophysics,2012,55(1):207-218.
[11] 楊宇,黃建平,雷建設(shè)等.Lebedev網(wǎng)格黏彈性介質(zhì)起伏地表正演模擬.石油地球物理勘探,2016,51(4):698-706.
Yang Yu,Huang Jianping,Lei Jianshe et al.Numerical simulation of Lebedev grid for viscoelastic media with irregular free-surface.OGP,2016,51(4):698-706.
[12] 張煜,徐義賢,夏江海等.含流體孔隙介質(zhì)中面波的傳播特性及應(yīng)用.地球物理學(xué)報(bào),2015,58(8):2759-2778.
Zhang Yu,Xu Yixian,Xia Jianghai et al.Characteristics and application of surface wave propagation in fluid-filled porous media.Chinese Journal of Geophy-sics,2015,58(8):2759-2778.
[13] Yuan S Y,Wang S X,Sun W J et al.Perfectly matched layer on curvilinear grid for the second-order seismic acoustic wave equation.Exploration Geophysics,2014,45(2):94-104.
[14] 裴正林.任意起伏地表彈性波方程交錯(cuò)網(wǎng)格高階有限差分法數(shù)值模擬.石油地球物理勘探,2004,39(6):629-634.
Pei Zhenglin.Staggered-grid high-order finite-difference numerical simulation of elastic wave equation for any irregular surface.OGP,2004,39(6):629-634.
[15] Bohlen T,Saenger E H.Accuracy of heterogeneous staggered-grid finite-difference modeling of Rayleigh waves.Geophysics,2006,71(4):T109-T115.
[16] Aki K,Richards P G.Quantitative Seismology:Theory and Methods.W H Freeman,New York,1980.
[17] Levander A R.Fourth-order finite-difference P-SV seismograms.Geophysics,1988,53(11):1425-1436.
[18] Robertsson J O A.A numerical free-surface condition for elastic/viscoelastic finite-difference modeling in the presence of topography.Geophysics,1996,61(6):1921-1934.
[19] Crase E.High-order (space and time) finite-difference modeling of the elastic wave equation.SEG Technical Program Expanded Abstracts,1990,9:987-991.
[20] Hestholm S,Ruud B O.3-D finite-difference elastic wave modeling including surface topography.Geophysical Prospecting,1994,42(5):371-390.
[21] Kosloff D,Kessler D,F(xiàn)ilho A Q et al.Solution of the equations of dynamic elasticity by a Chebychev spectral method.Geophysics,1990,55(6):734-748.
[22] 王秀明,張海瀾.用于具有不規(guī)則起伏自由表面的介質(zhì)中彈性波模擬的有限差分算法.中國(guó)科學(xué)G輯,2004,34(5):481-493.
Wang Xiuming,Zhang Hailan.Modeling the seismic wave in the media with irregular free interface by the finite-difference method.Science in China (Series G),2004,34(5):481-493.
[23] Boore D.Finite difference methods for seismic wave propagation in heterogeneous materials∥Methods in Computational Physics:Advances in Research and Applications,Academic Press,London,1972,11:1-37.
[24] Zahradnik J,Priolo E.Heterogeneous formulations of elastodynamic equations and finite-difference schemes.Geophysical Journal International,1995,120(3):663-676.
[25] Graves R W.Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences.Bulletin of Seismological Society of America,1996,86(4):1091-1106.
[26] Mittet R.Free-surface boundary conditions for elastic staggered-grid modeling schemes.Geophysics,2002,67(5):1616-1623.
[27] Xu Y,Xia J,Miller R D.Numerical investigation of implementation of air-earth boundary by acoustic-elastic boundary approach.Geophysics,2007,72(5):SM147-SM153.
[28] 王周,李朝暉,龍桂華等.求解彈性波有限差分法中自由邊界處理方法的對(duì)比.工程力學(xué),2012,29(4):77-83.
Wang Zhou,Li Chaohui,Long Guihua et al.Comparison among implementations of free-surface boundary in elastic wave simulation using the finite-difference method.Engineering Mechanics,2012,29(4):77-83.
[29] Berg P,If F,Nielsen P et al.Analytical reference solutions,Advanced seismic modeling//Modeling the Earth for Oil Exploration (Helbig K Ed).Pergamon Press,New York,1994,421-427.
[30] 凡友華,劉家琦,肖柏勛.計(jì)算瑞利波頻散曲線的快速矢量傳遞算法.湖南大學(xué)學(xué)報(bào)(自然科學(xué)版),2002,29(5):25-30.
Fan Youhua,Liu Jiaqi,Xiao Boxun.Fast vector-transfer algorithm for computation of Rayleigh wave dispersion curves.Journal of Hunan University (Natural Sciences Edition),2002,29(5):25-30.
[31] Luo Y H,Xia J H,Miller R D et al.Rayleigh-wave dispersive energy imaging using a high-resolution linear Radon transform.Pure and Applied Geophysics,2008,165(5):903-922.
[32] 夏江海,高玲利,潘雨迪等.高頻面波方法的若干新進(jìn)展.地球物理學(xué)報(bào),2015,58(8):2591-2605.
Xia Jianghai,Gao Lingli,Pan Yudi et al.New findings in high-frequency surface wave method.Chinese Journal of Geophysics,2015,58(8):2591-2605.
*湖北省武漢市洪山區(qū)魯磨路388號(hào)中國(guó)地質(zhì)大學(xué)地球物理與空間信息學(xué)院,430074。 Email:xhsong@cug.edu.cn
本文于2017年1月2日收到,最終修改稿于2017年10月4日收到。
本項(xiàng)研究受國(guó)家自然科學(xué)基金項(xiàng)目(41574114、41174113)和中國(guó)地質(zhì)大學(xué)(武漢)教學(xué)實(shí)驗(yàn)室開放基金項(xiàng)目(SKJ2016098)聯(lián)合資助。
1000-7210(2017)06-1156-14
袁士川,宋先海,蔡偉,胡瑩,魯鵬.瑞雷波有限差分?jǐn)?shù)值模擬中不同自由表面邊界條件的對(duì)比.石油地球物理勘探,2017,52(6):1156-1169.
P631
A
10.13810/j.cnki.issn.1000-7210.2017.06.005
(本文編輯:朱漢東)

袁士川 碩士研究生,1990年生;2015年本科畢業(yè)于長(zhǎng)江大學(xué)勘查技術(shù)與工程專業(yè),獲學(xué)士學(xué)位;目前在中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院攻讀地球探測(cè)與信息技術(shù)專業(yè)碩士學(xué)位,主要從事復(fù)雜介質(zhì)瑞雷波正演模擬及特性分析等方面的學(xué)習(xí)和研究。