熊嫘，孫成禹，蔡瑞乾

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

0 引言

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

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

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

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

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

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

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

1.1 柱坐標下彈性波波動方程

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

(1)

(2)

(3)

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

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

(4)

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

(5)

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

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

(6)

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

1.2 ADE-PML邊界

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

1.2.1 直角坐標形式

以直角坐標系中x方向為例，在頻率域利用變換空間求導(dǎo)算子

(7)

將波動方程拉伸至復(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是計算點與計算邊界之間的距離；cx1、cx2、cx3均為比例系數(shù)； Δx和Δt分別為空間和時間步長。

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

(10)

(11)

1.2.2 柱坐標形式

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

(12)

(13)

1.3 彈性波波正演模擬

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

1.3.1 簡單模型

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

圖1 簡單模型

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

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

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

圖4 水平分量地震記錄

圖5 垂直分量地震記錄

(14)

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

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

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

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

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

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

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

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

1.3.2 Marmousi模型

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

圖9 局部Marmousi模型

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

圖10 兩種坐標系下水平分量1.5 s波場快照

圖11 兩種坐標系下水平分量地震記錄

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

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

2 Rayleigh面波波場正演模擬

2.1 自由表面與面波生成機理

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

(15)

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

2.2 模型試驗與分析

假設(shè)模型頂面為自由表面，均勻上覆地層以下是均勻基底，模型上方為真空。建立坐標軸xOz(或rOz)，為計算方便，假設(shè)x軸(或r軸)沿著頂面向右為正方向，z軸向下為其正方向。模型區(qū)域范圍是7000 m×3500 m，網(wǎng)格采樣步長為2 m×2 m，試驗時間采樣率為0.5 ms，總采樣時間為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可見，兩種坐標系下的波場快照和地震記錄都基本相同。這表明地震波場在柱坐標系和直角坐標系下有相同的波場模擬效果，進一步驗證柱坐標系能夠模擬地震波的傳播。但不同坐標系下波場增益范圍有明顯區(qū)別(參見色標數(shù)值)，柱坐標下的面波波場比直角坐標下小近兩個數(shù)量級，即坐標系變化對面波波場有明顯影響。

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

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

圖16 垂直分量地震記錄

圖17 兩種坐標系下水平分量地震記錄

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

圖18 地表4500 m處質(zhì)點的運動軌跡

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

圖19 垂直分量波形振幅隨炮點埋深變化曲線

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

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

3 結(jié)論

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

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

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

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