熊嫘，孫成禹，蔡瑞乾

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

0 引言

目前，面波勘探以其低成本、高效率等諸多優(yōu)勢成為淺地表地震勘探的熱門方法之一。而對面波基礎(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ù)值頻散問題。但基于頻散方程的面波正演方法只適用于波形和走時(shí)信息的計(jì)算，無法計(jì)算波場各分量的振幅值，更無法得到全波場。

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

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

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

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

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

1 柱對稱條件下彈性波波場正演模擬

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

柱坐標(biāo)系rθz中，r、z、θ分別對應(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))中，在柱對稱條件下各變量與θ無關(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)

此外，針對柱坐標(biāo)系下彈性波場在r-r0→0處出現(xiàn)的極點(diǎ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)槔?，在頻率域利用變換空間求導(dǎo)算子

(7)

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

ADE-PML邊界中，通過下式

(8)

得到迭代衰減波場的系數(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í)間步長。

由式(8)可知，邊界匹配層越厚，則衰減系數(shù)越小，邊界的吸收效果越好。因此，可通過改變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軸無限延伸，二維波場與坐標(biāo)y無關(guān)，垂直分量為z分量，水平分量為x分量。在柱坐標(biāo)系rθz下，假設(shè)模型關(guān)于炮點(diǎn)所在軸呈柱對稱分布，二維波場與坐標(biāo)θ無關(guān)，垂直分量為z分量，水平分量為r分量。

1.3.1 簡單模型

為驗(yàn)證柱坐標(biāo)系正演模擬的可行性，建立區(qū)域范圍為4000 m×4000 m、網(wǎng)格數(shù)為800×800，網(wǎ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)對應(yīng)檢波點(diǎn)。下同)。

圖1 簡單模型

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

圖2 水平分量0.5 s波場快照

圖3 垂直分量0.5 s波場快照

圖4 水平分量地震記錄

圖5 垂直分量地震記錄

(14)

為更直觀地觀察波場正演的模擬精度，提取兩種坐標(biāo)系下地震記錄第500道的檢波器信號(hào)進(jìn)行對比。通過圖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)下彈性波波場表現(xiàn)為球面波。通過該結(jié)果能夠解釋振幅差是由于二維柱坐標(biāo)下波場以球面形式擴(kuò)散，二維直角坐標(biāo)下波場以柱面形式擴(kuò)散，而在同等震源條件下球面波波前上一點(diǎn)能量比柱面上小。時(shí)差主要是二維直角坐標(biāo)的平面波近似導(dǎo)致的相位差。進(jìn)而驗(yàn)證二維直角坐標(biāo)下正演波場與實(shí)際三維分布不符。

結(jié)合劉君等[30]關(guān)于波動(dòng)方程坐標(biāo)變換會(huì)引入誤差的相關(guān)研究可知，柱坐標(biāo)和直角坐標(biāo)下波動(dòng)方程的數(shù)值計(jì)算結(jié)果會(huì)存在微小誤差，這是由于兩種坐標(biāo)系波場在每個(gè)時(shí)間片計(jì)算產(chǎn)生的誤差來源不同、累計(jì)程度也不同。本文的正演結(jié)果也驗(yàn)證這一結(jié)論。從圖中可見，二維柱坐標(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í)際波場傳播中能量的變化規(guī)律。

圖8 相對振幅隨炮檢距變化曲線

對比簡單模型中兩種二維坐標(biāo)系波場模擬的結(jié)果可知，柱坐標(biāo)不僅能夠表征現(xiàn)有二維直角坐標(biāo)下的彈性波場特征，而且能夠表征實(shí)際彈性波傳播的能量、相位特征。通過波場快照及波形曲線兩個(gè)角度的對比分析可見，二維柱坐標(biāo)與三維直角坐標(biāo)下的波場特征基本一致，但二維柱坐標(biāo)下的波場正演更節(jié)省計(jì)算成本，正演速度更快。因此，對簡單模型下利用柱坐標(biāo)系進(jìn)行波場正演研究具有準(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)格采樣步長為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模型

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

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

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

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

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

2 Rayleigh面波波場正演模擬

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

地表空氣與地球接觸面是一個(gè)強(qiáng)阻抗界面。由于空氣阻抗與固體地球阻抗相比非常小，故可將地表看作是自由表面[31-33]，利用矩形網(wǎng)格可直接描述水平自由表面。本文使用隱式差分法處理邊界條件，規(guī)避差分邊界溢出問題?；趶椥圆ú▌?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)格采樣步長為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可見，兩種坐標(biāo)系下的波場快照和地震記錄都基本相同。這表明地震波場在柱坐標(biāo)系和直角坐標(biāo)系下有相同的波場模擬效果，進(jìn)一步驗(yàn)證柱坐標(biāo)系能夠模擬地震波的傳播。但不同坐標(biāo)系下波場增益范圍有明顯區(qū)別(參見色標(biāo)數(shù)值)，柱坐標(biāo)下的面波波場比直角坐標(biāo)下小近兩個(gè)數(shù)量級(jí)，即坐標(biāo)系變化對面波波場有明顯影響。

圖14 水平分量1.5 s波場快照

圖15 垂直分量1.5 s波場快照

圖16 垂直分量地震記錄

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

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

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

此外，采用不同深度激發(fā)、同一點(diǎn)接收的方式研究模擬激發(fā)深度對面波能量的關(guān)系，得到面波振幅與激發(fā)深度之間的關(guān)系(圖19)。從圖中可見地表激發(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)能量越來越小，平面波沒有幾何擴(kuò)散效應(yīng)，其波前上一點(diǎn)能量不會(huì)發(fā)生改變。結(jié)合鄧瑞[35]利用勢函數(shù)求解波動(dòng)方程，其在均勻介質(zhì)中的計(jì)算結(jié)果顯示面波在水平方向具有衰減現(xiàn)象。上述結(jié)果說明柱坐標(biāo)的模擬結(jié)果體現(xiàn)三維空間中面波波場能量的幾何擴(kuò)散特征，符合實(shí)際情況。

3 結(jié)論

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

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

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

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