黏聲波高階傅里葉有限差分法參數(shù)優(yōu)化成像

肖世鵬，熊高君，袁夢雨，毛明秋，王勝藝，韋增濤

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

0 引言

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

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

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

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

1 方法原理

1.1 傅里葉有限差分

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

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

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

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

(1)

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

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

(3)

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

(4)

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

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

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

(5)

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

(6)

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

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

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

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

1.3 偏移算法

頻率空間域的單程波波場延拓公式為：

(7)

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

(8)

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

延拓完成后，對頻率空間域計算結果進行傅里葉反變換:

(9)

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

(10)

上述延拓過程是對每一個ω進行的。每步延拓之后對每一個ω求和，即得偏移結果。

2 模型試算

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

2.1 點繞射模型

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

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

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

2.2 層狀模型

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

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

2.3 Marmousi模型

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

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

a—Maimousi速度模型;b—2階FFD正演(Q=∞);c—2階優(yōu)化FFD正演(Q=∞);d—4階優(yōu)化FFD正演(Q=∞);e—4階優(yōu)化FFD正演(變Q介質)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模型成像結果對比Fig.5 Marmousi model imaging by different mothods

3 結論

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