肖世鵬，熊高君，袁夢(mèng)雨，毛明秋，王勝藝，韋增濤

(成都理工大學(xué) 地球物理學(xué)院，四川 成都 610000)

0 引言

目前的理論研究和實(shí)際觀測(cè)表明地球介質(zhì)普遍存在非完全彈性特征，介質(zhì)的黏滯效應(yīng)會(huì)顯著影響地震波的動(dòng)力學(xué)特征，包括波性改造、振幅衰減、頻率降低等。這種大地吸收效應(yīng)主要與巖石巖性、含流體性質(zhì)、飽和度以及滲透率等有關(guān)，也是造成高頻成分衰減的主要原因之一[1-2]，針對(duì)不同衰減理論也有著不同的補(bǔ)償方法[3-5]。黏彈性地層對(duì)地震波的吸收效應(yīng)直接導(dǎo)致了地震記錄中深層信息的模糊不清，地震資料分辨率的降低，并間接影響到深層地震數(shù)據(jù)的處理效果[6]。因此通過引入品質(zhì)因子和復(fù)值速度來考慮實(shí)際地層的衰減特性進(jìn)行黏聲波數(shù)值模擬，研究地震波傳播規(guī)律過程并分析衰減特征，具有一定的實(shí)際意義[7]。

自Stocks首次研究黏彈性介質(zhì)及其地震波傳播以來[8]，出現(xiàn)了許多描述黏彈性介質(zhì)的數(shù)學(xué)模型，主要包含Kelvin-Voigt模型、Maxwell模型、標(biāo)準(zhǔn)線性模型等[9]。Auld在1990年指出無論是一維空間還是三維空間的彈性介質(zhì)和黏彈性介質(zhì)，二者之間的本構(gòu)方程和彈性模量都遵循對(duì)應(yīng)規(guī)則。而Kelvin-Voigt模型對(duì)地下介質(zhì)的描述既簡(jiǎn)單又較為符合地下的實(shí)際情況。

黏聲波介質(zhì)相較完全彈性介質(zhì)引入了波的衰減和頻散，導(dǎo)致地震子波發(fā)生畸變[10-11]，因此在數(shù)值模擬方法中受到了較多的限制[12]，射線追蹤法只能保持波的運(yùn)動(dòng)學(xué)特征，不適于模擬衰減和頻散效應(yīng)；波動(dòng)方程法能全面反映地震波的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)特征，模擬精度高[13]。波動(dòng)方程法又分為有限差分法、有限元法、偽譜法和相移法。其中相移法具有簡(jiǎn)單、快速、對(duì)地層傾角無限制好等優(yōu)點(diǎn)[14]，且在頻率—波數(shù)域進(jìn)行更易于引入吸收衰減，但當(dāng)速度有橫向變化時(shí)，就只能得到近似的正演結(jié)果[15]。Stoffa等在相移法的基礎(chǔ)上引入背景速度概念，在頻率—空間域增加誤差校正項(xiàng)，提出分步傅里葉法(SSF)能適應(yīng)更強(qiáng)的速度橫向變化[16]；Ristow等提出了傅里葉有限差分法(FFD)，該方法在SSF的基礎(chǔ)上增加了有限差分補(bǔ)償項(xiàng)[17]；張金海等將FFD法推廣到了黏聲波介質(zhì)的正演領(lǐng)域[15]。

高階傅里葉有限差分法有著結(jié)合相移法的精準(zhǔn)和有限差分法能適應(yīng)速度橫向變化的優(yōu)點(diǎn)。本文用Kelvin-Voigt模型構(gòu)建黏聲波方程，使用梯度下降法對(duì)傅里葉有限差分算子中的有限差分項(xiàng)進(jìn)行了全局參數(shù)優(yōu)化，在不提高方程階次的情況下達(dá)到更高階方程的逼近效果，并用該方法所求近似解對(duì)復(fù)雜模型的零炮檢距地震記錄進(jìn)行了數(shù)值模擬。

1 方法原理

1.1 傅里葉有限差分

Kelvin(開爾文)黏彈性介質(zhì)單元體是由一個(gè)彈性單元與一個(gè)黏性單元體并聯(lián)組成[18]。三維空間時(shí)，介質(zhì)中的每個(gè)體積元就是一個(gè)Kelvin單元體，并且在立體空間與其他單元空間相互聯(lián)結(jié)在一起。如圖1所示。

圖1 開爾文模型Fig.1 Kelvin model

頻率空間域的二維單程波延拓方程為[19]：

P(x,zn+1,ω)=P(x,zn,ω)exp(ikzΔz) ，

(1)

式中：P(x,z,ω)為頻率空間域二維地震波場(chǎng)；ω為頻率；(x,z)為二維剖面的坐標(biāo)；kz為波場(chǎng)延拓算子；Δz為深度上的步長(zhǎng)。在已知黏彈性介質(zhì)中相速度V(ω)和品質(zhì)因子Q的情況下，可以完成黏彈性介質(zhì)的波場(chǎng)延拓和偏移成像[8]。

式中：kelas=ω/V，其中V近似為均勻彈性各向同性介質(zhì)中波的相速度；Q為品質(zhì)因子。將復(fù)波數(shù)k(ω)代入式(1)，進(jìn)行單程波分解，并利用分裂算法把高階近似方程分解為N-1個(gè)(N為高階方程的階數(shù))二階串聯(lián)方程，可極大程度地減少計(jì)算量[20]。將其中的根式進(jìn)行泰勒展開得黏彈性介質(zhì)中向下延拓算子kz。

(3)

考慮式(3)，將kz按照多參量全局優(yōu)化的方法拓展為：

(4)

式(4)為本文的高階多參量?jī)?yōu)化傅里葉有限差分算子，其中(an,bn,cn,dn)為待定系數(shù)。本文采用梯度優(yōu)化法確定。

1.2 系數(shù)優(yōu)化方法

為確定式(4)中待定系數(shù)的值，常規(guī)優(yōu)化方法僅考慮速度的對(duì)比度[21]。本文增加多階次速度對(duì)比度系數(shù)并考慮當(dāng)前頻率及延拓步長(zhǎng)等參量的影響，選取單平方根算子kz與真實(shí)值的誤差為目標(biāo)函數(shù)[22]，用多元非線性梯度下降法對(duì)待定系數(shù)做全局優(yōu)化求解。

(5)

對(duì)an、bn、cn及dn估計(jì)的優(yōu)化過程是使積分最小化的過程，所求殘差應(yīng)為最小值，即

(6)

上式中J為近似解與精確解的誤差總和；Φ為波傳播的最大角度，J與精確解的比稱為相對(duì)誤差[23]。優(yōu)化系數(shù)如表1所列；優(yōu)化結(jié)果如圖2。

表1 系數(shù)優(yōu)化Table 1 Coefficient optimization

圖2 相對(duì)誤差隨傳播角變化曲線(p=c/v=0.5)Fig.2 Relative error changes with propagation angle(p=c/v=0.5)

由圖2可以看出，2階優(yōu)化相比FFD法有著低角度的誤差，但在70°多出一個(gè)過零點(diǎn)；4階優(yōu)化平衡了低角度的誤差。

1.3 偏移算法

頻率空間域的單程波波場(chǎng)延拓公式為：

(7)

其中kz如式(3)所示。利用絕對(duì)穩(wěn)定的Crank-Nicolson差分格式[24]，式(7)可化簡(jiǎn)為：

(8)

上式中a為將所有常數(shù)縮減為的系數(shù)；“*”號(hào)表示的值是可調(diào)節(jié)的，并按旁側(cè)邊界條件來調(diào)節(jié)。矩陣(8)是一個(gè)三對(duì)角復(fù)數(shù)矩陣[25]。

延拓完成后，對(duì)頻率空間域計(jì)算結(jié)果進(jìn)行傅里葉反變換:

(9)

根據(jù)t=0成像得：

(10)

上述延拓過程是對(duì)每一個(gè)ω進(jìn)行的。每步延拓之后對(duì)每一個(gè)ω求和，即得偏移結(jié)果。

2 模型試算

為驗(yàn)證黏彈性介質(zhì)高角度傅里葉有限差分正演模擬的正確性和模擬效果，本文分別進(jìn)行了點(diǎn)繞射模型和Marmousi模型的數(shù)值模擬實(shí)驗(yàn)。

2.1 點(diǎn)繞射模型

圖3是分別用相移法和傅里葉有限差分法正演的零炮檢距地震記錄。點(diǎn)繞射模型由強(qiáng)烈橫向變速的簡(jiǎn)單均勻介質(zhì)組成，速度從左至右分別為4000m/s、6000m/s；縱橫向采樣點(diǎn)數(shù)為100×128；縱橫向采樣間隔均為20m；3個(gè)模擬點(diǎn)震源分別位于(33,64)、(50,64)和(66,64)處，記錄網(wǎng)格長(zhǎng)256，采樣間隔4ms。

a—均勻介質(zhì)中的點(diǎn)繞射模型;b—四階有限差分算子優(yōu)化正演;c—相移法正演;d—不同Q值波場(chǎng)變化對(duì)比a—point diffraction model;b—forward results of fourth-order optimized;c—phase shift method forward;d—comparison of different Q values圖3 點(diǎn)繞射模型Fig.3 Point diffraction model

從圖中可以明顯看出：①黏滯性對(duì)振幅、頻帶和波性都有著明顯的改造作用，隨著傳播時(shí)間的增加，子波的主峰振幅和高頻成分都出現(xiàn)明顯的衰減趨勢(shì)，且Q值越小，地下黏滯性特征越明顯，如圖3d。②由于常規(guī)的頻率—波數(shù)域相移法外推算子只計(jì)算水平界面的平均速度，考慮放射點(diǎn)在垂向上的移動(dòng)，故對(duì)于橫向速度有變化的地層不敏感，水平速度界面在零炮檢距模擬記錄上形態(tài)將保持不變，如圖3c，而傅里葉有限差分法兼顧相移法成像精度高、穩(wěn)定性好和計(jì)算速度快之外還能適應(yīng)速度的縱橫向任意變化，并在時(shí)間剖面上正確成像。

2.2 層狀模型

圖4為簡(jiǎn)單層狀模型零炮檢距合成地震記錄。其模型是由背斜構(gòu)造、向斜構(gòu)造及一個(gè)水平界面組成，3層介質(zhì)速度由上至下分別為3000m/s、4000m/s、5000m/s；縱、橫向采樣點(diǎn)數(shù)為100×128；縱橫向采樣間隔皆為20m；時(shí)間采樣間隔為4ms；記錄長(zhǎng)度為1024ms。在記錄圖4b上均可以清晰地看到背斜和向斜構(gòu)造產(chǎn)生的“蝴蝶結(jié)”形地震反射特征。由于背斜頂部凸界面的反射存在發(fā)散現(xiàn)象，如圖4b，分配到單位面積上的波的能量會(huì)減弱。凸度越大，埋藏越深，射線發(fā)散越嚴(yán)重，地震波的振幅也越小[26]。如圖4a所示的速度模型中，凹界面的向斜模型曲率中心剛好在地層下，這時(shí)射線發(fā)生交叉，同相軸出現(xiàn)回轉(zhuǎn)，最凹點(diǎn)處波的旅行路徑比兩邊短，因此在時(shí)間剖面上將出現(xiàn)凸起，如圖4b，向斜兩翼的反射特征會(huì)左右顛倒，形成回轉(zhuǎn)波，在實(shí)際為向斜的部位深處又出現(xiàn)一個(gè)“背斜”。

a—黏彈性介質(zhì)速度模型;b—系數(shù)優(yōu)化后傅里葉有限差分正演結(jié)果;c—傅里葉有限差分偏移結(jié)果;d—優(yōu)化系數(shù)后4階傅里葉有限差分偏移結(jié)果a—velocity model of viscoclastic media;b—parameter optimization forward results;c—migration results;d—migration results of fourthorder optimized圖4 層狀模型Fig.4 Layered model

圖4c為常規(guī)傅里葉有限差分法的偏移結(jié)果，圖中因回轉(zhuǎn)波所造成的“背斜”消失，背斜構(gòu)造及向斜構(gòu)造都得到完全歸位，這表明常規(guī)FFD法能很好地解決橫向速度變化不大及頻散現(xiàn)象的問題。圖4d是優(yōu)化系數(shù)后的傅里葉有限差分法偏移結(jié)果，與圖4c相比，背斜構(gòu)造及向斜構(gòu)造完全歸位證明了該方法的正確性，成像結(jié)果更干凈，佐證了該方法的優(yōu)越性。兩種偏移方法所得結(jié)果基本相同，這和模型的復(fù)雜程度有關(guān)，簡(jiǎn)單的層狀模型橫向傾角變化低，屬于常規(guī)FFD法的傾角包含范圍，故復(fù)雜模型的成像結(jié)果能更好地體現(xiàn)兩種方法的差異性。

2.3 Marmousi模型

SEG的Marmousi模型存在背斜、尖滅和斷層等復(fù)雜地質(zhì)構(gòu)造和較強(qiáng)的橫向速度變化，是數(shù)值模擬中廣泛使用的標(biāo)準(zhǔn)模型之一。本文將使用該模型驗(yàn)證黏彈性介質(zhì)高階傅里葉有限差分法對(duì)較強(qiáng)橫向速度變化和Q值變化的適應(yīng)能力。

圖5是用不同方法對(duì)Marmousi模型進(jìn)行實(shí)驗(yàn)和對(duì)比。速度模型網(wǎng)格點(diǎn)數(shù)為300×941；網(wǎng)格間距10m；模型最小相速度1000m/s，最大相速度2600m/s，其背景參考為每層最小相速度；所用零炮檢距記錄共941道，每道1022個(gè)時(shí)間采樣點(diǎn)，采樣間隔為4ms。由圖5b可以看出，常規(guī)二階傅里葉有限差分法可以較為準(zhǔn)確地對(duì)地下界面進(jìn)行簡(jiǎn)單成像，如地層分界面及隆起進(jìn)行成像，但對(duì)模型復(fù)雜構(gòu)造及陡傾角斷層存在較大誤差，圖5c為本文二階優(yōu)化法(N=1)，其計(jì)算效率與常規(guī)二階傅里葉有限差分法相當(dāng)，在不提高方程階次的情況下兼顧陡傾角誤差，明顯提高了成像信噪比。圖5d為本文四階優(yōu)化法(N=2)，偏移結(jié)果具有較高的信噪比，各復(fù)雜構(gòu)造邊界清晰，反射界面清晰，陡傾角斷層及下伏隆起界面同相軸連續(xù)且?guī)缀醪淮嬖跉w位誤差(如圖5c和圖5d箭頭所示)，成像精度進(jìn)一步提高。圖5e為本文四階優(yōu)化法在縱向變Q值介質(zhì)中的應(yīng)用效果，Q值由淺至深取20～100，對(duì)比完全彈性介質(zhì)可知，其振幅和相位特征只有在傳播時(shí)間短且距地面距離較近處的部位才較吻合，表明黏滯吸收作用會(huì)對(duì)地震波波性、相位及振幅等產(chǎn)生很大影響，隨著傳播距離的增加，波場(chǎng)高頻成分明顯衰減，波性逐漸變寬、能量逐漸變?nèi)?如圖5d和圖5e紅框所示)，這些規(guī)律直接反映了黏滯性吸收的基本特征及本文方法在黏聲介質(zhì)中的應(yīng)用效果。

a—Maimousi速度模型;b—2階FFD正演(Q=∞);c—2階優(yōu)化FFD正演(Q=∞);d—4階優(yōu)化FFD正演(Q=∞);e—4階優(yōu)化FFD正演(變Q介質(zhì))a—Marmousi velocity model;b—second-order FFD forward(Q=∞);c—second-order optimized FFD forward(Q=∞);d—fourth-order optimized FFD forward(Q=∞);e—fourth-order optimized FFD forward(changing Q value)圖5 不同方法的Marmousi模型成像結(jié)果對(duì)比Fig.5 Marmousi model imaging by different mothods

3 結(jié)論

本文使用的黏聲波高階傅里葉有限差分法是一種混合域的計(jì)算方法，其波場(chǎng)外推算子由相移項(xiàng)、折射項(xiàng)(時(shí)移項(xiàng))和有限差分項(xiàng)組成。用梯度法近似確定的差分項(xiàng)外推算子，可在階數(shù)不變的情況下求解高階非線性方程組，適應(yīng)強(qiáng)空間變速介質(zhì)。數(shù)值實(shí)驗(yàn)的對(duì)比分析表明，優(yōu)化后的黏聲波高階傅里葉有限差分法相比常規(guī)FFD法及低階優(yōu)化的FFD法有著頻散低、模擬結(jié)果反射特征清楚、復(fù)雜地質(zhì)結(jié)構(gòu)更精確成像等優(yōu)點(diǎn)。并通過黏聲模型模擬結(jié)果對(duì)比，充分說明了優(yōu)化后黏聲波高階傅里葉有限差分法成像的優(yōu)越性。本文方法可以處理二維復(fù)雜地質(zhì)構(gòu)造中波傳播問題，一般情況下二階優(yōu)化差分項(xiàng)已達(dá)到成像精度要求，四階有限差分項(xiàng)反而會(huì)影響計(jì)算效率，且在三維模型中差分項(xiàng)存在的雙向分裂會(huì)引起方位各向異性誤差，直接影響成像精度。