熊嫘，孫成禹，蔡瑞乾

(中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580)

0 引言

目前，面波勘探以其低成本、高效率等諸多優(yōu)勢(shì)成為淺地表地震勘探的熱門方法之一。而對(duì)面波基礎(chǔ)理論做充分研究是面波勘探達(dá)到理想效果的前提和條件。1885年，Rayleigh[1]首次在理論上證明Rayleigh面波的存在，自此面波成為地球物理學(xué)家研究的熱點(diǎn)之一。面波具有能量強(qiáng)，傳播遠(yuǎn)，淺層分辨率高，抗干擾能力強(qiáng)等優(yōu)點(diǎn)，且在層狀介質(zhì)中表現(xiàn)出頻散特性。由于面波的頻散曲線能反映地層的速度及厚度，基于頻散方程的數(shù)值模擬算法得以快速發(fā)展。Haskell[2]提出Haskell正演數(shù)值算法，并推導(dǎo)出Rayleigh波頻散方程。Abo-Zena[3]提出Abo-Zena算法，解決了Thomson-Haskell及Knopoff算法出現(xiàn)的高頻數(shù)值頻散問(wèn)題。但基于頻散方程的面波正演方法只適用于波形和走時(shí)信息的計(jì)算，無(wú)法計(jì)算波場(chǎng)各分量的振幅值，更無(wú)法得到全波場(chǎng)。

基于波動(dòng)方程的數(shù)值模擬方法是另一種實(shí)現(xiàn)面波正演的重要手段。其中，有限差分?jǐn)?shù)值模擬法[4-9]能模擬自由表面邊界條件[10-11]下的三維波場(chǎng)及可視化地震波在介質(zhì)中的傳播過(guò)程[12]。Mittet[13]提出基于交錯(cuò)網(wǎng)格有限差分的自由邊界處理方案，但生成的Rayleigh波出現(xiàn)嚴(yán)重的頻散現(xiàn)象。周竹生等[14]使用有限差分實(shí)現(xiàn)二維全波場(chǎng)正演，揭示了面波的形成機(jī)理和傳播規(guī)律。

在波場(chǎng)模擬中，體波波場(chǎng)是由點(diǎn)震源激發(fā)產(chǎn)生的三維球狀波場(chǎng)。但由于三維數(shù)值計(jì)算的工作量大，對(duì)硬件要求高，其計(jì)算效率難以滿足實(shí)際生產(chǎn)的需要。因此，為提高計(jì)算效率，通常在二維條件下研究三維波場(chǎng)的相關(guān)特性。

在二維波動(dòng)方程正演中，Rayleigh面波常以平面縱波與平面垂直極化橫波的干涉波場(chǎng)形式出現(xiàn)，Love波以平面水平極化橫波形式出現(xiàn)，理論研究中所使用的波函數(shù)也是平面波形式。傅承義[15]在平面波的基礎(chǔ)上利用位移函數(shù)和波動(dòng)方程探索面波的頻散關(guān)系及衰減現(xiàn)象。但實(shí)際觀測(cè)到的面波，無(wú)論是其縱波分量還是其橫波分量，都不可能表現(xiàn)為平面波形式。在水平層狀介質(zhì)條件下，面波各分量均表現(xiàn)為柱面波形式，這與其在二維直角坐標(biāo)中呈現(xiàn)的平面波形式明顯不同。因此，在二維直角坐標(biāo)下正演所得的波場(chǎng)存在如下問(wèn)題： ①模擬得到的面波與體波的振幅變化分別為平面波與柱面波特征，而實(shí)際的面波與體波分別為柱面波與球面波，兩者之間不相符； ②模擬得到的面波以平面波形式向外傳播時(shí)，其振幅不具備空間幾何衰減特征，而實(shí)際中的面波是以柱面波形式向外傳播，其振幅存在空間的幾何衰減特征。因此，基于平面波理論的二維直角坐標(biāo)下的波動(dòng)方程對(duì)波場(chǎng)的描述存在不足，正演結(jié)果難以反映真實(shí)波場(chǎng)特征。

本文使用柱對(duì)稱條件處理線震源所激發(fā)的柱面波，在二維平面內(nèi)實(shí)現(xiàn)球面波的等價(jià)模擬。其中，柱坐標(biāo)由于自然貼合柱面網(wǎng)格及其柱對(duì)稱性質(zhì)，多用于流體運(yùn)動(dòng)研究[16]和儀器制造[17]。相較于地震波場(chǎng)正演中直角坐標(biāo)的廣泛應(yīng)用，柱坐標(biāo)的相關(guān)研究較少。何柏榮等[18]利用二維柱坐標(biāo)研究圓柱狀地質(zhì)體中地震波的繞射問(wèn)題。張碧星等[19]引入柱坐標(biāo)系分析三維水平層狀介質(zhì)的Rayleigh波頻散特征。

本文使用柱對(duì)稱條件將三維問(wèn)題簡(jiǎn)化為二維問(wèn)題，首先推導(dǎo)了柱坐標(biāo)系下的彈性波波動(dòng)方程，利用彈性波的基礎(chǔ)理論，分析正演所得彈性波場(chǎng)，驗(yàn)證使用二維算法實(shí)現(xiàn)三維波場(chǎng)模擬的合理性和可行性。然后分別計(jì)算直角坐標(biāo)系和柱坐標(biāo)系下的面波波場(chǎng)，對(duì)比兩種坐標(biāo)系下面波波場(chǎng)振幅、波形及走時(shí)特征，分析柱對(duì)稱條件下面波的傳播特性和模擬優(yōu)勢(shì)。

1 柱對(duì)稱條件下彈性波波場(chǎng)正演模擬

1.1 柱坐標(biāo)下彈性波波動(dòng)方程

柱坐標(biāo)系rθz中，r、z、θ分別對(duì)應(yīng)水平坐標(biāo)、垂直坐標(biāo)和角度坐標(biāo)。有如下關(guān)系

(1)

(2)

(3)

式中：t為時(shí)間；εm為正應(yīng)變；ηmn為切應(yīng)變(m,n=r、θ、z,m≠n)；σmm為正應(yīng)力；τmn為切應(yīng)力；u、e、w為直角坐標(biāo)系下的位移分量，都是r、θ、z、t的函數(shù)；λ、μ是拉梅常數(shù)；ρ是介質(zhì)密度。

將柱坐標(biāo)下的幾何方程[20](式(1))和本構(gòu)方程(式(2))代入柱坐標(biāo)下的三維運(yùn)動(dòng)平衡方程[21](式(3))中，在柱對(duì)稱條件下各變量與θ無(wú)關(guān)，波動(dòng)方程滿足?/?θ=0，得到

(4)

縱橫波速度與彈性參數(shù)的轉(zhuǎn)換關(guān)系為

(5)

式中vP、vS分別是縱波、橫波速度。

在rOz平面中，震源位于r0處，r-r0表示計(jì)算點(diǎn)與震源之間的水平距離。將式(5)代入式(4)，即得到柱坐標(biāo)下的二維彈性波波動(dòng)方程

(6)

此外，針對(duì)柱坐標(biāo)系下彈性波場(chǎng)在r-r0→0處出現(xiàn)的極點(diǎn)問(wèn)題。本文在差分處理有關(guān)r-r0的方程時(shí)，將常規(guī)差分網(wǎng)格沿徑向偏移半個(gè)網(wǎng)格點(diǎn)[22]以避免r-r0→0處出現(xiàn)數(shù)值奇異的情況。

1.2 ADE-PML邊界

本文選用基于輔助微分方程(ADE)的完全匹配層(PML)方法處理吸收邊界[23]，即為ADE-PML邊界算法。在二維情況下，使用該吸收邊界算法分別處理質(zhì)點(diǎn)的垂直位移分量和水平位移分量。

1.2.1 直角坐標(biāo)形式

以直角坐標(biāo)系中x方向?yàn)槔陬l率域利用變換空間求導(dǎo)算子

(7)

將波動(dòng)方程拉伸至復(fù)數(shù)域。式中sx=dx/jω，是x方向的復(fù)拉伸變量，其中j為虛數(shù)單位，ω是角頻率，dx是人工吸收邊界的衰減參數(shù)。

ADE-PML邊界中，通過(guò)下式

(8)

得到迭代衰減波場(chǎng)的系數(shù)

(9)

式中：dx是x方向的衰減系數(shù)；vmax是彈性波的最大傳播速度；L為PML邊界的厚度；R是理論反射系數(shù)(設(shè)為0.001)；l是計(jì)算點(diǎn)與計(jì)算邊界之間的距離；cx1、cx2、cx3均為比例系數(shù)； Δx和Δt分別為空間和時(shí)間步長(zhǎng)。

由式(8)可知，邊界匹配層越厚，則衰減系數(shù)越小，邊界的吸收效果越好。因此，可通過(guò)改變PML的厚度靈活控制邊界的吸收效果[24-29]。

(10)

(11)

1.2.2 柱坐標(biāo)形式

基于Ma等[23]提出的ADE-PML算法，將其推廣至柱坐標(biāo)系。以柱坐標(biāo)下波動(dòng)方程的水平分量r為例，將波動(dòng)方程的水平分量拉伸至復(fù)數(shù)域

(12)

(13)

1.3 彈性波波正演模擬

在直角坐標(biāo)系xyz下，假設(shè)模型沿y軸無(wú)限延伸，二維波場(chǎng)與坐標(biāo)y無(wú)關(guān)，垂直分量為z分量，水平分量為x分量。在柱坐標(biāo)系rθz下，假設(shè)模型關(guān)于炮點(diǎn)所在軸呈柱對(duì)稱分布，二維波場(chǎng)與坐標(biāo)θ無(wú)關(guān)，垂直分量為z分量，水平分量為r分量。

1.3.1 簡(jiǎn)單模型

為驗(yàn)證柱坐標(biāo)系正演模擬的可行性，建立區(qū)域范圍為4000 m×4000 m、網(wǎng)格數(shù)為800×800，網(wǎng)格采樣步長(zhǎng)為5 m×5 m的模型。為保證模擬結(jié)果具有可比性，該剖面位于相同規(guī)格的三維直角坐標(biāo)模型的y=0處。時(shí)間采樣率為0.5 ms，總采樣時(shí)間為0.5 s。介質(zhì)的縱、橫波速度分別為4000、2000 m/s，介質(zhì)密度為2000 kg/m3。采用ADE-PML吸收邊界。在柱坐標(biāo)系和直角坐標(biāo)系中分別施加主頻為20 Hz的雷克子波，并將震源置于模型中心位置(圖1，其中爆破符號(hào)表示震源，紅點(diǎn)對(duì)應(yīng)檢波點(diǎn)。下同)。

圖1 簡(jiǎn)單模型

結(jié)合彈性波在均勻各向同性介質(zhì)中的傳播特征，分別對(duì)比圖2～圖5可知，無(wú)論是波場(chǎng)快照還是炮集記錄，二維柱坐標(biāo)與三維直角坐標(biāo)的正演結(jié)果都一致。而在二維模型中，雖然兩種坐標(biāo)系下彈性波各分量的波場(chǎng)快照和地震記錄結(jié)果都保持一致，但兩種坐標(biāo)系下二維波場(chǎng)增益強(qiáng)度相差兩個(gè)數(shù)量級(jí)(參見(jiàn)色標(biāo)數(shù)值)。而二維柱坐標(biāo)與三維直角坐標(biāo)的增益強(qiáng)度是相同的。由此可知，采取不同對(duì)稱系統(tǒng)的二維坐標(biāo)系對(duì)正演波場(chǎng)的能量有明顯影響。

圖2 水平分量0.5 s波場(chǎng)快照

圖3 垂直分量0.5 s波場(chǎng)快照

圖4 水平分量地震記錄

圖5 垂直分量地震記錄

(14)

為更直觀地觀察波場(chǎng)正演的模擬精度，提取兩種坐標(biāo)系下地震記錄第500道的檢波器信號(hào)進(jìn)行對(duì)比。通過(guò)圖6、圖7可看出：不同試驗(yàn)中相同位置處接收到的波形基本相同，其中三維直角坐標(biāo)和二維柱坐標(biāo)下波的振幅與相位基本一致，二維直角坐標(biāo)下波的相位明顯滯后于其他兩條波形曲線。比較二維波形曲線最大振幅處的數(shù)據(jù)可知：兩條曲線振幅大小相差約兩個(gè)數(shù)量級(jí)，走時(shí)差為0.0075 s。

圖6 第500道接收的水平分量

圖7 第500道接收的垂直分量

因此，均勻介質(zhì)中的彈性波在二維直角坐標(biāo)下表現(xiàn)為柱面波形式，二維柱坐標(biāo)下彈性波波場(chǎng)表現(xiàn)為球面波。通過(guò)該結(jié)果能夠解釋振幅差是由于二維柱坐標(biāo)下波場(chǎng)以球面形式擴(kuò)散，二維直角坐標(biāo)下波場(chǎng)以柱面形式擴(kuò)散，而在同等震源條件下球面波波前上一點(diǎn)能量比柱面上小。時(shí)差主要是二維直角坐標(biāo)的平面波近似導(dǎo)致的相位差。進(jìn)而驗(yàn)證二維直角坐標(biāo)下正演波場(chǎng)與實(shí)際三維分布不符。

結(jié)合劉君等[30]關(guān)于波動(dòng)方程坐標(biāo)變換會(huì)引入誤差的相關(guān)研究可知，柱坐標(biāo)和直角坐標(biāo)下波動(dòng)方程的數(shù)值計(jì)算結(jié)果會(huì)存在微小誤差，這是由于兩種坐標(biāo)系波場(chǎng)在每個(gè)時(shí)間片計(jì)算產(chǎn)生的誤差來(lái)源不同、累計(jì)程度也不同。本文的正演結(jié)果也驗(yàn)證這一結(jié)論。從圖中可見(jiàn)，二維柱坐標(biāo)與三維直角坐標(biāo)的正演結(jié)果間存在著約0.00125 s的微小走時(shí)誤差，這正是波動(dòng)方程坐標(biāo)變換的離散采樣造成的。

由圖8可知，隨炮檢距增加，彈性波振幅在柱坐標(biāo)下比在二維直角坐標(biāo)下衰減得更快。試驗(yàn)結(jié)果滿足同等震源條件下球面擴(kuò)散比柱面擴(kuò)散速度更快的傳播規(guī)律，進(jìn)一步驗(yàn)證柱坐標(biāo)下正演結(jié)果符合實(shí)際波場(chǎng)傳播中能量的變化規(guī)律。

圖8 相對(duì)振幅隨炮檢距變化曲線

對(duì)比簡(jiǎn)單模型中兩種二維坐標(biāo)系波場(chǎng)模擬的結(jié)果可知，柱坐標(biāo)不僅能夠表征現(xiàn)有二維直角坐標(biāo)下的彈性波場(chǎng)特征，而且能夠表征實(shí)際彈性波傳播的能量、相位特征。通過(guò)波場(chǎng)快照及波形曲線兩個(gè)角度的對(duì)比分析可見(jiàn)，二維柱坐標(biāo)與三維直角坐標(biāo)下的波場(chǎng)特征基本一致，但二維柱坐標(biāo)下的波場(chǎng)正演更節(jié)省計(jì)算成本，正演速度更快。因此，對(duì)簡(jiǎn)單模型下利用柱坐標(biāo)系進(jìn)行波場(chǎng)正演研究具有準(zhǔn)確、高效的優(yōu)點(diǎn)。

1.3.2 Marmousi模型

為進(jìn)一步驗(yàn)證柱坐標(biāo)系應(yīng)用的可行性，采用Marmousi模型。網(wǎng)格尺寸為248×248，網(wǎng)格采樣步長(zhǎng)為5 m×5 m。時(shí)間采樣率為0.5 ms，總采樣時(shí)間為1.5 s。邊界為ADE-PML吸收邊界。在柱坐標(biāo)系和直角坐標(biāo)系中分別施加主頻為20 Hz雷克子波。為適應(yīng)實(shí)際的炮集處理流程，模擬單邊接收，將震源放置在模型左上邊界處(圖9)，紅點(diǎn)表示檢波點(diǎn)。

圖9 局部Marmousi模型

分析對(duì)比圖10、圖11，可知在兩種坐標(biāo)系下復(fù)雜模型的1.5 s全波場(chǎng)快照和地震記錄除波場(chǎng)增益外差別不大。由圖12a可見(jiàn)，全波場(chǎng)中直達(dá)波能量過(guò)強(qiáng)無(wú)法清晰得到反射波波形，且柱坐標(biāo)下反射波振幅小于直角坐標(biāo)下的振幅。為進(jìn)一步研究復(fù)雜模型下的反射波，分析去除直達(dá)波后的相對(duì)波形記錄(圖12b)。在復(fù)雜模型中柱坐標(biāo)下波的相位仍滯后于直角坐標(biāo)，且隨著時(shí)間增加，柱坐標(biāo)下波形振幅衰減比直角坐標(biāo)快。復(fù)雜模型下的反射波場(chǎng)與簡(jiǎn)單模型下的正演結(jié)果一致。

圖10 兩種坐標(biāo)系下水平分量1.5 s波場(chǎng)快照

圖11 兩種坐標(biāo)系下水平分量地震記錄

圖12 第30道反射波水平分量

至此，通過(guò)均勻介質(zhì)模型和Marmousi模型的正演波場(chǎng)，驗(yàn)證柱坐標(biāo)系對(duì)兩種模型下均能保證正演效果，實(shí)現(xiàn)波場(chǎng)正演，并體現(xiàn)三維彈性波波場(chǎng)的能量擴(kuò)散特征及相位特征。

2 Rayleigh面波波場(chǎng)正演模擬

2.1 自由表面與面波生成機(jī)理

地表空氣與地球接觸面是一個(gè)強(qiáng)阻抗界面。由于空氣阻抗與固體地球阻抗相比非常小，故可將地表看作是自由表面[31-33]，利用矩形網(wǎng)格可直接描述水平自由表面。本文使用隱式差分法處理邊界條件，規(guī)避差分邊界溢出問(wèn)題?；趶椥圆ú▌?dòng)方程(式(6))，在自由邊界處利用下式生成面波[5，12，34]

(15)

式中：i為自由表面任一橫向網(wǎng)格點(diǎn)位置；γ=λ/(λ+2μ)。

2.2 模型試驗(yàn)與分析

假設(shè)模型頂面為自由表面，均勻上覆地層以下是均勻基底，模型上方為真空。建立坐標(biāo)軸xOz(或rOz)，為計(jì)算方便，假設(shè)x軸(或r軸)沿著頂面向右為正方向，z軸向下為其正方向。模型區(qū)域范圍是7000 m×3500 m，網(wǎng)格采樣步長(zhǎng)為2 m×2 m，試驗(yàn)時(shí)間采樣率為0.5 ms，總采樣時(shí)間為1.5 s。模型中上覆地層深10 m，上覆地層縱波波速為2000 m/s，橫波波速為1100 m/s； 下伏地層縱波波速為2200 m/s，橫波波速為1200 m/s，介質(zhì)密度為2000 kg/m3。震源使用主頻25 Hz的雷克子波，置于上覆地層底面中心處(圖13)所示。

圖13 面波試算模型

由圖14～圖17可見(jiàn)，兩種坐標(biāo)系下的波場(chǎng)快照和地震記錄都基本相同。這表明地震波場(chǎng)在柱坐標(biāo)系和直角坐標(biāo)系下有相同的波場(chǎng)模擬效果，進(jìn)一步驗(yàn)證柱坐標(biāo)系能夠模擬地震波的傳播。但不同坐標(biāo)系下波場(chǎng)增益范圍有明顯區(qū)別(參見(jiàn)色標(biāo)數(shù)值)，柱坐標(biāo)下的面波波場(chǎng)比直角坐標(biāo)下小近兩個(gè)數(shù)量級(jí)，即坐標(biāo)系變化對(duì)面波波場(chǎng)有明顯影響。

圖14 水平分量1.5 s波場(chǎng)快照

圖15 垂直分量1.5 s波場(chǎng)快照

圖16 垂直分量地震記錄

圖17 兩種坐標(biāo)系下水平分量地震記錄

抽取4500 m處質(zhì)點(diǎn)在不同時(shí)刻的位移(圖18)，可見(jiàn)質(zhì)點(diǎn)沿逆時(shí)針?lè)较蜻\(yùn)動(dòng)，其方向與波的傳播方向相反，其軌跡形狀近似為主軸垂直的橢圓。該結(jié)果表明地表質(zhì)點(diǎn)運(yùn)動(dòng)軌跡與理論一致，模擬所得波場(chǎng)為面波。

圖18 地表4500 m處質(zhì)點(diǎn)的運(yùn)動(dòng)軌跡

此外，采用不同深度激發(fā)、同一點(diǎn)接收的方式研究模擬激發(fā)深度對(duì)面波能量的關(guān)系，得到面波振幅與激發(fā)深度之間的關(guān)系(圖19)。從圖中可見(jiàn)地表激發(fā)時(shí)會(huì)產(chǎn)生較強(qiáng)的面波，隨著炮點(diǎn)埋深的增加，產(chǎn)生的面波振幅迅速減小。從圖16中提取面波的垂直振幅，得到振幅隨炮檢距的變化關(guān)系如圖20所示。

圖19 垂直分量波形振幅隨炮點(diǎn)埋深變化曲線

圖20 垂直分量波形振幅隨炮檢距變化曲線

可以看出，隨炮檢距的增加，柱坐標(biāo)系下模擬得到的面波振幅發(fā)生衰減，而直角坐標(biāo)系下模擬得到的面波振幅基本不變。由于二維柱坐標(biāo)下面波以柱面形式傳播，二維直角坐標(biāo)下面波以平面形式傳播，而隨著波前的擴(kuò)散，柱面波波前上一點(diǎn)能量越來(lái)越小，平面波沒(méi)有幾何擴(kuò)散效應(yīng)，其波前上一點(diǎn)能量不會(huì)發(fā)生改變。結(jié)合鄧瑞[35]利用勢(shì)函數(shù)求解波動(dòng)方程，其在均勻介質(zhì)中的計(jì)算結(jié)果顯示面波在水平方向具有衰減現(xiàn)象。上述結(jié)果說(shuō)明柱坐標(biāo)的模擬結(jié)果體現(xiàn)三維空間中面波波場(chǎng)能量的幾何擴(kuò)散特征，符合實(shí)際情況。

3 結(jié)論

本文研究實(shí)現(xiàn)柱對(duì)稱條件下彈性波傳播及面波正演，并對(duì)比二維柱坐標(biāo)系與二維直角坐標(biāo)系下的地震波場(chǎng)。模型試驗(yàn)和理論分析結(jié)果表明：

(1)在均勻模型和Marmousi模型中，兩種坐標(biāo)系下模擬結(jié)果的波形基本一致，即利用柱坐標(biāo)實(shí)現(xiàn)二維波場(chǎng)模擬具有可行性；

(2)在均勻模型中，二維直角坐標(biāo)下波形曲線比二維柱坐標(biāo)的相位超前45°、振幅小兩個(gè)數(shù)量級(jí)。結(jié)合理論分析可驗(yàn)證二維直角坐標(biāo)下的模擬波場(chǎng)在能量和相位特征上與實(shí)際不符。而利用柱坐標(biāo)系可等效模擬三維空間中地震波傳播的振幅幾何衰變，在能量和相位的變化上具備基本的三維特征，可提高正演效率，降低三維模擬的計(jì)算成本。

本文提出利用柱坐標(biāo)實(shí)現(xiàn)波場(chǎng)正演的方法能夠推動(dòng)對(duì)實(shí)際波場(chǎng)的相關(guān)研究，有可能獲得新的認(rèn)識(shí)，或新的實(shí)際應(yīng)用方案，但其相對(duì)于常規(guī)二維直角坐標(biāo)的波場(chǎng)在計(jì)算上更復(fù)雜。因此，在后續(xù)的研究中可考慮將柱坐標(biāo)方法拓展到層狀介質(zhì)正反演中，進(jìn)一步發(fā)揮柱坐標(biāo)的優(yōu)勢(shì)，提高計(jì)算效率。