基于上下行波分解的分數(shù)階拉普拉斯算子黏滯聲波逆時偏移成像方法

王喜鵬，張慶臣，毛偉建*

（1.中國科學院精密測量科學與技術(shù)創(chuàng)新研究院計算與勘探地球物理研究中心，湖北武漢 430077；2.大地測量與地球動力學國家重點實驗室，湖北武漢 430077； 3.中國科學院大學，北京 100049）

0 引言

近年來，基于雙程波的逆時偏移成像技術(shù)突破了成像角度的限制，對陡傾構(gòu)造成像具有顯著優(yōu)勢［1-2］。但地震波在傳播過程中，由于地下介質(zhì)通常具有黏滯性，導致其能量衰減及相位畸變，尤其當勘探目標區(qū)域存在氣云等強衰減構(gòu)造時，對其做偏移成像會導致分辨率降低和照明不足［3］。為此，學者們進行了大量的研究，模擬地震波在地下介質(zhì)傳播過程中的衰減效應(yīng)。

起初，對黏性逆時偏移的研究是基于近似常Q模型［4］波動方程展開的，通過單個或多個力學元件的串、并聯(lián)或疊加（如標準線性體、Maxwell體等）描述地下介質(zhì)的衰減效應(yīng)［5-8］，但基于該模型推導的波動方程中包含記憶變量，需較多表征參數(shù)，這些參數(shù)還需要計算優(yōu)化選擇，額外增加計算量［9］。此外，在這類波動方程中，振幅衰減與相位頻散相耦合，在逆時偏移成像過程中不能準確地校正相位，會導致深部目標層成像的位置誤差?；诔模型［4］推導的波動方程中含有分數(shù)階時間導數(shù)，對波動方程計算當前時刻波場值時，需依賴該時刻以前所有時刻的波場，因此其計算成本高、效率低且內(nèi)存消耗巨大。

Chen 等［10］采用分數(shù)階拉普拉斯算子代替分數(shù)階時間導數(shù)，將時間導數(shù)轉(zhuǎn)換為空間導數(shù)，運用偽譜法在波數(shù)域求解黏聲波動方程，提高了計算效率。Treeby等［11］基于冪律衰減模型，推導了解耦的分數(shù)階拉普拉斯算子黏滯聲波波動方程，實現(xiàn)了振幅衰減與相位頻散的分離。劉文卿等［12］基于GPU/CPU 協(xié)同系統(tǒng)，通過GPU 實現(xiàn)波場逆時外推，并采用隨機速度邊界提高算法的并行性，解決了逆時偏移中大規(guī)模存儲的I/O 問題。Zhu 等［13］基于常Q頻散關(guān)系，推導了含解耦的分數(shù)階拉普拉斯算子時域黏聲波動方程，在波傳播過程中具有近似常Q的衰減和頻散特征，實現(xiàn)了振幅衰減與相位頻散的分離，在逆時偏移過程中對能量補償、相位校正十分便利，只需在波場外推過程中，反轉(zhuǎn)振幅衰減算子的符號且保持相位頻散算子符號不變，就可較準確地校正由地下介質(zhì)黏滯性引起的振幅衰減與相位頻散效應(yīng)，并基于此方程實現(xiàn)了Q補償逆時偏移（Q-compensated reverse-time migration，Q-RTM）［14］。

然而，采用Q-RTM 在對地震波振幅進行衰減補償時，也會增強數(shù)據(jù)中的高頻噪聲，在波場外推過程中出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象。此前，高頻噪聲通常采用低通Tukey 濾波器壓制［14-15］。但該方法是時不變的，不能適應(yīng)在空間變化的Q值和深度。

Treeby［16］提出一種頻域時變Tukey 濾波器以穩(wěn)定衰減補償過程，但采用低通濾波法會損失地震數(shù)據(jù)中部分有用高頻信號。Wang［17］通過引入穩(wěn)定因子，提出一種反Q偏移方法穩(wěn)定補償過程，并認為某深度處的地震波吸收衰減效應(yīng)應(yīng)是從地表到當前深度吸收衰減的累積； Wang 等［18］在此基礎(chǔ)上提出一種新的Q補償逆時偏移自適應(yīng)穩(wěn)定方法，并通過推導分數(shù)階拉普拉斯算子黏滯聲波波動方程及其補償方程和衰減方程的波數(shù)域格林函數(shù)，從理論上解釋了波場外推在高頻端不穩(wěn)定的原因。與傳統(tǒng)基于低通濾波的方法相比，該方法具有更強的時間不一致性和Q依賴性，提高了衰減補償穩(wěn)定性和抗噪性。陳漢明等［19］應(yīng)用最小二乘反演理論，推導了分數(shù)階拉普拉斯算子黏滯聲波方程對應(yīng)的Born正演模擬算子和伴隨方程，利用反演思路補償?shù)卣鸩ǖ奈账p，解決了補償不穩(wěn)定問題。Chen等［20］提出一種隱式穩(wěn)定策略，將時變?yōu)V波器隱式地納入振幅增強波場外推步驟，并在頻譜中加入一個穩(wěn)定因子，避免了Q-RTM 中任何顯式的濾波操作，且該方法不增加額外計算量。

對于逆時偏移成像產(chǎn)生的低頻噪聲，目前主要采用如角度衰減因子成像條件［21-22］、對成像結(jié)果進行濾波（最小二乘濾波器、拉普拉斯濾波等）［23-24］、波場方向分解［25-26］等方法壓制。杜啟振等［27］總結(jié)了各類成像方法，認為采用波場分解的成像方法在復雜構(gòu)造區(qū)能更好地消除低頻噪聲。Liu 等［25-26］提出上下行波場分解的方法，構(gòu)造新的偏移成像條件，實現(xiàn)不同向的震源波場和檢波點波場的互相關(guān)成像，從而消除低頻噪聲。Fei等［28］指出速度模型中存在強速度梯度或反射界面時，震源下行波場和檢波點上行波場在逆時偏移成像過程中可能產(chǎn)生假象，成像時需略去這兩項。

此外，常規(guī)波場分解是在f-k域?qū)崿F(xiàn)，這需將震源波場和檢波點波場進行存儲，然后再進行時空域的傅里葉變換，因此會產(chǎn)生巨大的計算量和I/O 量。針對該問題，胡江濤等［29］提出一種基于聲波波動方程的解析時間波場外推方法，并將其運用于VSP（Vertical Seismic Profile）數(shù)據(jù)的逆時偏移成像。該方法在時間外推過程中實現(xiàn)了波傳播方向的顯式分解，避免了巨大的I/O 量和計算量，實現(xiàn)了低頻噪聲與有效信號的分離，并能有效地消除圖像中由復雜介質(zhì)引起的假像。Zhou 等［30］采用一步法波場外推實現(xiàn)了Q補償逆時偏移波場外推，并基于該方法波場為解析信號的性質(zhì)實現(xiàn)了波場的分離，運用Q補償逆時偏移和分解的波場成像極大程度提高了偏移成像的質(zhì)量，但該方法采用的是傳統(tǒng)低通濾波器穩(wěn)定化衰減補償過程，且一步法外推的計算速度較慢。

目前，Q-RTM 主要采用成像后濾波的方式壓制低頻噪聲，但該方法會改變成像剖面中同相軸振幅大小和波數(shù)成分?；诓▓龇纸獾某上穹椒芨行У貕褐齐p程波逆時偏移成像中的低頻噪聲和速度模型中存在較強速度梯度時產(chǎn)生的假象。

綜上所述，本文采用基于常Q模型解耦的分數(shù)階拉普拉斯算子黏聲波動方程對震源波場和檢波點波場進行外推，避免由于時間歷史記憶變量而需存儲每個時刻的波場，實現(xiàn)衰減算子和頻散算子的分離，可通過反轉(zhuǎn)衰減算子符號和保持頻散算子符號不變達到補償能量和校正頻散的目的，并基于該方程實現(xiàn)Q補償逆時偏移算法。

采用自適應(yīng)穩(wěn)定因子方案穩(wěn)定化衰減補償過程中數(shù)值不穩(wěn)定問題。對于逆時偏移產(chǎn)生的低頻噪聲和假象，采用解析時間波場外推方法，在時間外推過程中實現(xiàn)波傳播方向的顯式分解，選取分解后的震源下行波場和檢波點上行波場進行互相關(guān)成像，減少由復雜介質(zhì)及強速度梯度引起的偏移成像中的低頻噪聲和假象。

1 方法原理

1.1 常Q 模型應(yīng)力—應(yīng)變關(guān)系

在一維線性非彈性介質(zhì)中，應(yīng)力與應(yīng)變的關(guān)系［31］為

式中：ψ1-2γ為弛豫函數(shù)；t為時間；t0=1/ω0為參考時間；γ=arctan(1/Q)/π 為無量綱參數(shù)，且0＜γ＜1/2，Q值越小介質(zhì)的衰減效應(yīng)越強，當Q→∞時，γ→0，方程退化為聲波方程；為體積模量，其中ρ0為密度，c0為參考頻率ω0處的相速度；*為卷積算子；ε為應(yīng)變場。

弛豫函數(shù)定義［9］為

式中：Γ為歐拉伽馬函數(shù)；H(t)為Heaviside 階躍函數(shù)。

1.2 解耦的分數(shù)階拉普拉斯算子波動方程

對于均勻聲波介質(zhì)，一階線性動量守恒方程為

應(yīng)變速度方程為

兩式中：σ為應(yīng)力場；v為速度。

結(jié)合式（3），可得時間分數(shù)階波動方程

式中： ?2為拉普拉斯算子；。由于時間分數(shù)階導數(shù)的存在，計算某時刻波場時需存儲該時刻之前每一時刻的波場值，內(nèi)存需求較大，且數(shù)值計算困難。為此，Chen等［10］提出用分數(shù)階拉普拉斯算子(-?2)α/2代替分數(shù)階時間導數(shù)?α/?tα，分數(shù)階拉普拉斯算子可在空間域計算，避免較大的內(nèi)存需求及數(shù)值計算難題。

將平面波解exp(iωt-ik·r)代入式（5），可得對應(yīng)的頻散方程

式中：k=kR+ikI為復波數(shù)矢量；ω為角頻率；r為空間坐標矢量。

由i2γ=cos(πγ)+i sin(πγ)，復波數(shù)代入k，經(jīng)一系列計算及近似，可得新的頻散方程

分數(shù)階拉普拉斯算子的傅里葉變換為

式中：k為空間域波數(shù)； F表示傅里葉變換。

對式（7）進行時空二維反傅里葉變換，可得

1.3 Q-RTM 原理

逆時偏移有三個主要步驟： ①震源波場外推；②檢波點波場外推； ③采用適當?shù)某上駰l件對外推的震源和檢波點波場進行成像，如互相關(guān)成像條件［32］

式中：S(t，x)為震源波場；R(t，x)為檢波點波場；Tmax為波場外推的最大時間值。

Zhu等［14］基于式（9）實現(xiàn)Q補償逆時偏移，并給出統(tǒng)一的震源及檢波點波場外推的黏聲波動方程

式中：p為壓力場；c0為聲波速度；L和H分別為分數(shù)階拉普拉斯算子； {ηL-?2}p和Hp分別為頻散和振幅衰減算子；L=(-?2)γ+1；H=(-?2)γ+1/2。該方程中β1和β2對應(yīng)項分別為頻散算子和振幅損失算子，在正向外推過程中，β2=1 表示振幅衰減，反向外推過程中，時間t被-t代替，衰減項符號反轉(zhuǎn)β2=-1，表示放大振幅。只需保持頻散算子符號不變（β1=1），就可校正頻散效應(yīng)。若β1和β2都為0，則式（11）退化為聲波方程； 若兩者都為1，該方程為黏聲波方程。而補償過程則通過反轉(zhuǎn)衰減算子的符號，使β1=1且β2=-1，式（11）變換為

本文通過式（12）對震源波場和檢波點波場進行外推，校正地下介質(zhì)對地震波造成的能量衰減和頻散效應(yīng)，并采用規(guī)則網(wǎng)格有限差分法計算時間導數(shù)、偽譜法計算拉普拉斯算子。

1.4 自適應(yīng)穩(wěn)定因子

Wang等［18］提出一種新的Q補償逆時偏移自適應(yīng)穩(wěn)定補償策略，推導過程如下。

Q-RTM 的衰減時間傳播算子可表示為

補償時間傳播算子可表示為

根據(jù)上兩式，振幅衰減算子定義為

振幅補償算子定義為

式（16）由于指數(shù)項的存在，它在時間上是發(fā)散的，因此定義如下穩(wěn)定化的振幅補償算子

式中σ2為穩(wěn)定因子。

由于振幅補償算子已體現(xiàn)在波場外推過程中，則穩(wěn)定因子可表示為

考慮當前時刻的衰減效應(yīng)為零時刻到當前時刻衰減效應(yīng)的累積，則每個時間步長穩(wěn)定因子系數(shù)定義為

第n個時間步穩(wěn)定因子系數(shù)為

1.5 上下行波場分解

由于采用雙程波動方程，震源波場和檢波點波場中同時包含上行波和下行波，運用成像條件使波路徑上不必要的震源波場和檢波點波場（震源上行波場與檢波點上行波場，震源下行波場與檢波點下行波場）互相關(guān)會在成像結(jié)果上產(chǎn)生低頻噪聲。此外，F(xiàn)ei等［28］指出當速度模型精度不高，存在較強速度梯度時，在特殊地質(zhì)構(gòu)造區(qū)，震源上行波場和檢波點下行波場成像時可能會產(chǎn)生假象。

為避免上述情況，本文采用的成像條件為

式中： 下標u（upgoing）表示上行波； d（downgoing）表示下行波，如Sd則表示震源下行波場。

在f-k域通過簡單的計算很容易分離上行波與下行波分量［33］

式中：S(ω，kz)為震源波場S(t，z)的二維傅里葉變換；kz為波數(shù)；Su(ω，kz)和Sd(ω，kz)分別為f-k域分解的上行波和下行波分量。由該式可知，波場方向由頻率與波數(shù)乘積的正負號決定。

Liu 等［26］通過固定波數(shù)符號，將震源和檢波點波場分為波數(shù)大于零和小于零的四個分量，使波數(shù)符號相反的震源波場與檢波點波場相乘，再將相乘后的兩個分量相加，最后使用互相關(guān)成像條件成像（如震源波場波數(shù)符號為正，則與之相乘的檢波點波場波數(shù)符號為負。若頻率為正，則分別為震源上行波場及檢波點下行波場；若頻率為負，則分別為震源下行波場和檢波點上行波場），實現(xiàn)傳播方向相反的震源和檢波點波場的互相關(guān)成像。但該方法中上行波場與下行波場互相耦合，不能將波場進行顯式地分解。

為了實現(xiàn)波場的顯式分解，需盡量使用一個變量（頻率或波數(shù)）的符號來確定波場傳播的方向。由于在時間域波場外推中，時間步長上波數(shù)信息較易獲取，若能獲取一個其頻譜僅為正或負的波場，則波場的方向分解就可通過波數(shù)的正、負號確定。為此，胡江濤等［29］提出一種解析時間波場外推方法。解析波場頻譜只包含正頻率，可僅依靠空間波數(shù)的符號實現(xiàn)波場方向的分解。解析波場包括實部和虛部兩部分，實部由信號本身構(gòu)成，虛部由經(jīng)希爾伯特變換后的原信號構(gòu)成。

震源的時間解析波場可表示為

基于式（22）和解析波場僅包含正頻率的性質(zhì)，只需根據(jù)波數(shù)的符號就可達到對波場進行方向分解的目的。震源的上、下行波場表示為

式中： Re表示實部； F-1表示傅里葉反變換；表示解析時間震源波場關(guān)于z方向的傅里葉變換。同理，檢波點波場的上、下行波表示為

值得注意的是，由于采用解析波場外推分離上下行波場，與傳統(tǒng)逆時偏移成像方法相比多了震源波場外推、檢波點波場外推及時間切片上的波場傅里葉變換，因此該方法計算時間會高于傳統(tǒng)Q補償逆時偏移成像。實驗表明，本文方法的計算時間約為采用拉普拉斯濾波Q-RTM 方法的1.3倍。但與傳統(tǒng)的波場分解相比，該方法規(guī)避了I/O 量過大的弊端，解決了波場存儲問題，且提高了計算效率。

2 數(shù)值算例

2.1 兩層模型

為了測試本文方法的有效性，首先采用水平層狀模型對解析波場分離上下行波的有效性進行測試。圖1a 所示速度模型中，第一、第二層層速度分別為2000、3000 m/s，模型網(wǎng)格尺寸為300×300，x和z方向的網(wǎng)格間距為dx=dz=10 m，得到該層狀模型對應(yīng)的Q模型（圖1b）。圖2a為1.0 s時刻未分離的波場快照，用解析波場分解方法得到該時刻經(jīng)波場分離后的下行波場（圖2b）和上行波場（圖2c）。可見本文方法能有效分離上下行波場。

圖1 兩水平層速度模型（a）及其Q 模型（b）

圖2 1.0 s 時刻波場快照

2.2 抽稀的Sigsbee 模型

為驗證采用的新成像方法與Q補償逆時偏移對振幅補償和相位校正的有效性及對復雜模型的適應(yīng)性，采用抽稀后的Sigsbee 模型做數(shù)值模擬測試。速度模型（圖3a）及對應(yīng)的Q模型（圖3b，由經(jīng)驗公式得出）的大小為1067×401，網(wǎng)格間距為dx=dz=10 m，時間采樣間隔為1 ms，采樣點數(shù)為4001，震源采用主頻為30 Hz的Ricker子波，總共104 炮，均勻分布在深度50 m 處，炮間距為100 m。

圖3 Sigsbee 速度模型（a）及其Q 模型（b）

圖4a 為2.0 s 時刻未分離的波場快照，用解析波場分解方法分別得到該時刻經(jīng)波場分離后的下行波場（圖4b）和上行波場（圖4c）。

圖4 2.0 s 時刻波場快照

圖5a 為未經(jīng)波場分離的傳統(tǒng)聲波逆時偏移成像結(jié)果，分別得到波場分離后的聲波逆時偏移成像結(jié)果（圖5b，參考）、波場分離后的未補償衰減數(shù)據(jù)逆時偏移成像結(jié)果（圖5c）及Q補償逆時偏移成像結(jié)果（圖5d）。圖5b右側(cè)子圖為選取水平距離8000 m 處作為參考的單道信息（黑色實線），圖5c、圖5d 右側(cè)子圖則分別為水平距離8000 m 處（圖5c虛線位置）相應(yīng)的單道信息和參考單道信息（黑色實線）。圖6為從上述子圖中截取的0.6～1.6 km 層段單道信息。其對應(yīng)的圖5b～圖5d的振幅譜如圖7a～圖7c所示。

圖5 偏移成像結(jié)果

圖6 截取的0.6～1.6 km 層段單道信息

圖7 不同逆時偏移的振幅譜

與圖5b（聲波）相比，圖5c（衰減）偏移成像層位及鹽丘下邊界能量減弱；通過對應(yīng)的單道信息子圖與振幅譜對比，與作為參考的聲波數(shù)據(jù)相比，衰減數(shù)據(jù)的振幅減小，相位發(fā)生畸變，頻帶變窄。而圖5d（Q補償逆時偏移）中層位及鹽丘下邊界更清晰； 通過與參考圖像（5b）截取的子圖及振幅譜對比，可見其損失的振幅與產(chǎn)生的相移都得到一定程度校正，與參考圖像吻合度較高，頻帶得到拓寬。從圖5a與圖5b 可知，采用傳統(tǒng)的逆時偏移成像方法在成像剖面上存在嚴重的噪聲及假象，降低成像精度，而采用上、下行波場分離成像的偏移結(jié)果則避免了這種影響，提高了分辨率。

2.3 抽稀的Pluto 模型

為進一步驗證Q補償逆時偏移對深層能量補償?shù)男Ч芭c傳統(tǒng)的拉普拉斯濾波方法對偏移成像噪聲壓制的效果，采用抽稀后的Pluto 模型進行數(shù)值模擬測試。其速度模型（圖8a）及對應(yīng)的Q模型（圖8b，由經(jīng)驗公式得出）的大小為1160×201，網(wǎng)格間距為dx=dz=10 m，時間采樣間隔為0.5 ms，采樣點數(shù)為4002，震源采用主頻為30 Hz的Ricker子波，總共114炮，炮間距為100 m，均勻分布在深度50 m 處。

圖8 Pluto 速度模型（a）及其Q 模型（b）

圖9a 為1.0 s 時刻未分離的波場快照，用上述方法分別得到該時刻經(jīng)波場分離后的下行波場（圖9b）和上行波場（圖9c）波場快照。

圖9 1.0s 時刻波場快照

圖10a為未經(jīng)波場分離的傳統(tǒng)聲波逆時偏移成像結(jié)果，分別得到聲波逆時偏移經(jīng)拉普拉斯濾波后成像結(jié)果（圖10c）、Q補償逆時偏移經(jīng)拉普拉斯濾波后成像結(jié)果（圖10e）、波場分離后的聲波逆時偏移成像結(jié)果（圖10b，參考）、波場分離后的衰減數(shù)據(jù)逆時偏移成像結(jié)果（圖10d）及波場分離后的Q補償逆時偏移成像結(jié)果（圖10f）。圖10b右側(cè)子圖為水平距離6000 m 處作為參考的單道信息（黑色實線），圖10d、圖10f 右側(cè)子圖則分別為水平距離6000 m 處相應(yīng)的單道信息（紅色虛線）和參考單道信息（黑色實線）。圖11a～圖11c對應(yīng)為10b、圖10d、圖10f 中水平距離5～8 km、深度1.3～1.8 km 范圍的振幅譜。

圖10 不同的逆時偏移成像結(jié)果

圖11 圖10d 中紅框?qū)?yīng)位置振幅譜

與圖10b（參考）相比，圖10d（衰減未補償）因地層的衰減效應(yīng)，其能量減弱，偏移成像層位不清晰，尤其是鹽下邊界以下； 圖10f 由于采用Q補償逆時偏移，其能量得到補償。對比10b、圖10d、圖10f 的單道信息與圖11a～圖11c 的振幅譜，可見緣于地層的非完全彈性效應(yīng)，對成像結(jié)果造成振幅衰減及相位畸變，而采用Q補償逆時偏移（圖11c）能在一定程度上減緩這種效應(yīng)的影響，拓寬其振幅譜，補償衰減的能量，并對畸變的相位進行校正，提高了成像質(zhì)量。由圖10a可知，采用傳統(tǒng)的互相關(guān)偏移方法，成像結(jié)果中淺層的有效信息完全湮沒在噪聲中，嚴重干擾成像精度；而采用拉普拉斯濾波后的成像結(jié)果（圖10c、圖10e）中低頻噪聲受到一定程度壓制，但對淺層信息仍有較明顯的干擾。

綜上所述，本文采用的Q補償逆時偏移方法能有效補償由地下介質(zhì)非彈性特性造成的振幅損失和相移，采用該成像方法所得偏移成像結(jié)果（圖10d～圖10f）中，有效地壓制由傳統(tǒng)成像條件所引起的低頻噪聲及假象。

值得注意的是，該方法是在進行兩次震源波場和檢波點波場外推的情況下，還要在時間切片上對波場進行傅里葉變換，所以計算時間稍長于采用拉普拉斯濾波的Q-RTM，測試結(jié)果（表1）顯示約為其1.3 倍。但該方法不用對波場進行存儲，避免了巨大的I/O量。

表1 不同方案下的計算時間

3 結(jié)論

本文提出一種上、下行波分解的解耦分數(shù)階拉普拉斯算子黏滯聲波逆時偏移成像方法：采用解耦的分數(shù)階拉普拉斯算子黏滯聲波方程對震源波場和檢波點波場進行外推，能顯著減緩由于地下介質(zhì)的黏滯性所引起的振幅衰減和頻散效應(yīng)，提高了偏移成像的照明能量；采用分數(shù)階拉普拉斯算子代替分數(shù)階時間導數(shù)，不需在求解波動方程計算當前時刻波場值時依賴該時刻以前所有時段的波場。采用解析時間波場能顯式地分解上下行波場。新的成像條件能避免同向傳播的震源波場和檢波點波場在應(yīng)用互相關(guān)成像條件時所產(chǎn)生的低頻噪聲，消除由于偏移速度模型存在較強速度梯度時由震源上行波場和檢波點下行波場互相關(guān)產(chǎn)生的假象，從而提高偏移成像的分辨率和質(zhì)量。