穩定Q補償梯度的黏滯聲波全波形反演

蔣書琦，周輝*，陳漢明，張明坤，富禹鑫，李紅輝

（1.中國石油大學（北京）油氣資源與工程全國重點實驗室，北京 102249；2.中國石油大學（北京）CNPC 物探重點實驗室，北京 102249；3.中國石油大學（北京）地球物理學院，102249；4.東方地球物理公司物探技術研究中心，河北涿州 072751）

0 引言

準確的地震波速度在地震勘探中起著至關重要的作用，是地震偏移、巖性識別和儲層預測成功的關鍵。作為目前最先進的高分辨率速度建模方法，全波形反演(FWI)已持續得到學術界和產業界的關注。通常在給定的初始速度模型基礎上，FWI以合成地震記錄與觀測記錄差異作為目標函數，計算該目標函數對模型參數的梯度，再結合優化算法迭代更新速度模型參數。FWI 迭代時所用梯度通常由伴隨狀態方法計算得到，即通過正傳波場導數與反傳波場進行互相關獲取[1-2]。然而，當地震波在衰減介質中傳播時，常規的FWI 沒有考慮地震波衰減，無法恢復可靠的模型速度。而梯度隨深度衰減的黏滯聲波FWI 受到深部梯度變小的影響，收斂速度緩慢。本文致力于開發一種梯度補償方法以提高黏滯聲波FWI 的收斂效率和反演結果的精度。

計算效率低一直制約著FWI 的發展。目前加速全波形反演的策略有兩類。一類加速策略是從優化迭代方法的角度出發，利用高階優化算法獲得更為優越的更新方向，如二階優化算法中的高斯—牛頓法和截斷牛頓法[3-4]。但這些方法每次迭代至少需要四次波動方程正演計算更新方向，引入不容小覷的額外計算量。因此高階優化方法務須尋求收斂速度和計算量之間的平衡。此類方法中采用線性搜索技術的LBFGS(Limited-memory Broyden-Fletcher-Goldfarb-Shanno)優化算法是性價比最高的二階優化算法[5-7]。同時，Castellanos 等[8]指出，L-BFGS 方法結合逆 Hessian 近似矩陣的附加信息預處理可以獲得收斂更快的預條件逆Hessian 矩陣，即為P-L-BFGS(Preconditioned LBFGS)算法。

另一類加速策略是通過預處理獲得收斂更快的梯度，如考慮地震波本身的能量特征[9-10]或結合其他信息如地層傾信息[11]進行改進。介質的黏滯效應不僅反映在地震波的能量上[12]，也會影響全波形反演目標函數的梯度。從這一角度出發，通過補償衰減為黏滯聲波FWI 提供更適宜收斂的梯度不失為一種有效的解決策略。該策略的核心在于如何選擇合適的衰減補償方法。Wang 等[13]提出了一種用于衰減補償的自適應穩定逆時偏移方法，其中穩定因子可以用關于補償峰值的經驗公式顯式地表示。與傳統的低通濾波穩定方法相比，自適應穩定策略可以高度提升反向時間外推方法的數值穩定性和抗噪性。Zhao等[14]采取了在完全彈性介質和黏彈性介質中通過外推波場得到的兩個算子實現補償，但該方法受到了特定成像條件的限制。Sun 等[15-16]通過構建穩定的補償算子，在不同條件下進行波場外推，但計算負擔巨大。此外，Chen 等[17]將正則化方法隱式地嵌入波動方程正演模擬中，避免了顯式濾波所增加的額外計算量。

描述黏滯效應的模型大體可以分為基于廣義標準線性體(SLS)[18-21]和基于常Q模型兩類，后者更常用?；诔模型，學者們提出了一系列描述黏滯介質中地震波特性的波動方程[11-12,22-25]。Zhu 等[26]推導了具有解耦變分數階拉普拉斯算子的常Q黏滯聲波方程。該方程可以用來補償振幅衰減和校正相位畸變，非常適用于衰減補償的逆時偏移。Chen 等[27]對該變分數階解耦方程進行了改進，優化了方程的表示方式，獲得了近似常Q固定分數階拉普拉斯算子黏滯波動方程。該方程保留了原有波方程的解耦性，大大降低了方程的數值計算復雜性，但并沒有明顯降低黏滯波場的數值模擬精度。在該方程的基礎上進行黏滯聲波方程的最小二乘逆時偏移可以穩定地補償介質的黏滯性，獲得高分辨率的成像結果[28]。因此，本文基于該常分數階拉普拉斯算子常Q黏滯聲波方程研究Q補償穩定化方案，以獲得穩定波場補償，進而獲得精度高且收斂更快的全波形反演目標函數的梯度，進而提高黏滯聲波方程全波形反演的收斂速度。本文進行速度反演所使用的Q模型是由其他方法獲得的，假設已知而不參與反演。

本文首先給出常分數階拉普拉斯算子黏滯波動方程，闡明該方程的數值解法；再回顧FWI 理論，分析地層黏滯性對黏滯聲波FWI的影響，說明衰減補償對反演的重要性；然后闡述基于該波動方程的補償方案，獲得含補償項的補償黏滯聲波方程，用于伴隨波場的計算；最后應用數值實例驗證基于該補償方案計算的梯度能夠加快黏滯聲波FWI的收斂。

1 方法原理

1.1 常分數階拉普拉斯算子黏滯聲波方程

Zhu等[26]提出的解耦分數階拉普拉斯算子黏滯聲波方程為

式中：p(x,t)為t時刻、空間x處的聲波波場；?2為拉普拉斯算子；γ、c、τ、η定義為

由式(3)可以看出，α和β分別控制著波場的振幅和相位，且振幅隨時間呈指數衰減。在實際情況中，Q和c隨空間變化，故γ也是空間變化的，振幅衰減和相位畸變的程度也會發生對應的變化。應用式(1)進行波場模擬，需要計算變分數階拉普拉斯算子，這將大大降低計算效率。若采用分數階的平均值轉化為常分數階拉普拉斯方程，又將降低模擬的精度。為此采用 Chen 等[27]的簡化方程進行數值模擬和全波形反演。該簡化方程兼具可操作性和較高數值精度的優點，具體為

1.2 基于伴隨狀態法的目標函數梯度計算

全波形反演通常選用L2范數作為目標函數，其梯度通過伴隨狀態法計算獲得[10,29]。若先將式(5)簡單表示為

式中F表示為黏滯聲波波動方程系統，則全波形反演的梯度可以寫為

伴隨波場padj由式(5)對應的伴隨方程獲得。

根據式(5)，式(7)可以寫為

由式（7）和式（8）可知，梯度受震源波場和伴隨波場的黏滯效應影響。若對波場實施衰減補償，則由此計算的梯度能量也會增強。梯度迭代更新方法中，作為二階優化算法的L-BFGS 算法非常有效，因為它能直接計算梯度與逆 Hessian 的乘積，從而通過無矩陣遞歸算法直接確定模型更新量。加入預條件算子的P-L-BFGS 方法相比L-BFGS 方法能更快收斂，故本文采用P-L-BFGS算法迭代更新模型。

1.3 基于穩定化策略的穩定衰減補償方法

由于補償是衰減的逆過程，則式(3)對應的補償形式應為

由上式可知，在補償過程中地震波場的振幅呈現指數型增長，容易導致數值不穩定。為確保波場的穩定補償，本文提出一種含穩定因子的Q補償方法。引入穩定因子σ的波場補償解析表達式為

當σ取0時，式(10)退化為式(9)，穩定化程度隨著σ的增大而增強。利用三角函數的周期性，式(10)可等價為

將式(11)整理成與震源波場傳播類似的表達形式

根據式（12），可以獲得C1和C2的表達式

C2可進一步簡化為

式中

將式(13)由波數域轉換到空間域，即獲得時空域的波動方程。同時利用

可進一步簡化。由于τc隨空間變化，為簡便計算式(17)，取最大值τcmax替代在空間域變化的τc。將最終簡化的C1和C2代入式(13)，得到最終的含穩定因子的波場傳播方程。同時參照 Chen 等[27]文獻中的第一種思路，可得到最終時空域補償方程。該補償方程除最后一項外，其他項與式(5)相同，即右端○取正號且G(p)=Csp，其中Cs為C在空間域的對應響應。

補償穩定效果取決于穩定因子σ2的選取。當穩定因子取0時，該方程退化為無穩定策略的補償方程。穩定因子越大，對補償過程中波場能量就有越強的約束作用，但如果選取過大，可能會損失波場信息。為適應不同黏滯性的衰減模型，通常選用σ2=rexp(-τcmaxΔt)，其中r為常數，一般令r=10-2。

2 數值模擬算例

2.1 穩定衰減補償策略的有效性測試

為了驗證本文衰減補償方法的有效性，選用具有強衰減區的二維BP氣囪模型(圖1)進行測試。該模型網格數為398×161，空間采樣間隔為10 m。Q模型由速度模型根據經驗公式Q=3.516c2.2生成，其中速度c的單位為km/s。該模型在上部的氣囪區吸收非常強，當地震波穿過該區域時會快速衰減，導致下層的反射信號非常弱，應用常規方法在氣囪區下方無法獲得可靠的反演結果。

圖1 BP 氣囪模型

初始速度模型(圖1c)的頂層60 m 設為準確速度，其下速度由準確模型光滑而得，只保留了速度變化趨勢。以200 m 間隔在模型頂部布設20炮，首炮位于模型左邊界右側90 m 處。接收方式為雙邊可變觀測系統，檢波器的位置隨震源位置移動而移動，檢波點間距為10 m。設置最大炮檢距為1800 m，以期望接收到一定角度范圍內的深層反射信息。震源子波選用主頻為10 Hz 的 Ricker 子波，參考角頻率為200 rad/s。為了充分記錄波場信息及滿足穩定性條件，設置記錄的時間樣點數為2501，采樣間隔設置為0.8 ms。

選用σ2=10-2exp(-τcmaxΔt)為補償時的穩定因子，比較補償前、后的反傳波場及其表征補償效果。圖2 為第10 炮不同時刻補償前、后反傳波場快照對比，可見未經衰減補償的波場振幅相對較弱，而經補償后恢復了振幅信息，在高衰減區域明顯增強，表明該衰減補償算法能實現有效的振幅補償，且補償后的反傳波場沒有不穩定現象。

圖2 BP 氣囪模型第10 炮不補償（左）與補償（右）不同時刻反傳波場的切片對比

2.2 基于穩定梯度衰減補償的反演測試

為了驗證穩定Q補償對黏滯聲波FWI的有效性，利用BP氣囪模型進行反演測試。圖3為第一次迭代的梯度。由梯度趨勢可以看出，Q補償后，梯度的整體趨勢大體不變，但模型中淺層尤其是高衰減區域的梯度能量明顯增強，反映的結構信息也更清晰。

圖3 BP 氣囪模型不補償（a）與補償（b）黏滯聲波FWI 第一次迭代的梯度對比

采用炮點平方預處理方式，P-L-BFGS法迭代20、40 和100 次反演結果如圖4所示?？梢钥闯?，基于梯度Q補償的全波形反演在初始階段的迭代就可以恢復大體的速度模型，而利用未補償的梯度無法在反演初期獲得大體的速度趨勢。雖然隨著迭代次數的增加，梯度未補償黏滯聲波全波形反演效果也得到改善，但即便迭代至100 次，仍然未能反演出準確的速度。對比x=1000、1500 和3200 m 的最終反演的速度縱向曲線(圖5)可以看出，基于梯度Q補償的全波形反演不僅能獲得更精確的層速度，而且較清晰地揭示了深部的大傾角傾斜地層，驗證了本文方法在加快黏滯聲波全波形反演的良好效果。

圖4 BP 氣囪模型不補償（左）與補償（右）不同迭代次數的反演結果對比

圖5 BP 氣囪模型不同位置的黏滯聲波FWI 第100 次反演速度曲線對比

2.3 不同補償穩定因子選取測試

穩定因子σ2=rexp(-τcmaxΔt)中r分別取不同值時，第10炮0.4s 時刻的反傳波場如圖6所示。對比圖2c與圖6 可以看出，選用的穩定因子都可以使補償的反傳波場保持穩定，同時補償后的反傳波場均較未補償的反傳波場能量有明顯增強，大穩定因子比小穩定因子補償后的波場能量弱。

圖6 r 取不同值時第10 炮0.4 s 時刻的反傳波場對比

圖7展示了由圖6所示的反傳波場計算的第一次迭代的梯度。對比圖7 與圖3 可見：與小穩定因子相比，大穩定因子補償的梯度數值較小，但依然大于未補償波場計算的梯度。

圖7 r 取不同值時BP 氣囪模型的黏滯聲波FWI 第一次迭代的梯度

反演的最大迭代次數設為100，如果目標函數停止下降也可提前終止迭代。r取不同值的反演結果如圖8所示。對比最終反演結果(圖4c與圖8)可以看出，由于衰減效應，梯度未補償的反演結果存在明顯的錯誤。相對而言，Q補償的反演結果更優越?；诓煌姆€定因子進行梯度Q補償的反演結果不盡相同。當補償穩定因子過大時，如r=10-1，補償過程被抑制，獲得的反演結果（圖8a）反而出現偏差，只迭代17次就陷入局部極小，無法繼續迭代。當補償穩定因子過小時，如r=10-32，反傳過程出現不穩定現象，終止于第66次迭代，獲得極不準確的反演結果(圖8d)。當補償穩定因子較小但處于合理范圍時，既能保證反演順利進行同時又保證補償過程穩定，如r=10-2、r=10-4及r=10-16，迭代100 次的反演效果（圖4c 右、圖8b和圖8c）均比未補償時的結果更準確。

圖8 r 取不同值時BP 氣囪模型的黏滯聲波FWI 結果

應用平均相對誤差定量評價反演結果的精度

圖9為不同r迭代17次和66次的反演結果誤差對比。從圖9a 可以看到，當r=10-1時，補償穩定因子過大，抑制了波場振幅，獲得的反演結果偏差大。相對于傳統未補償的反演結果，使用較小穩定因子的梯度補償方法的反演結果的誤差更小。由圖9b可以看出，當補償穩定因子過小時，如r=10-32，即便最終的反演誤差低于傳統未補償的誤差，但由于未能保證補償過程的穩定，出現明顯的數值不穩定現象，影響反演結果。對不同穩定因子的反演結果進行誤差分析可以看出，隨著穩定因子的增大，反演誤差呈現出先降后升的現象。在全波形反演過程中進行梯度補償時，應選用不抑制補償、同時保證數值計算穩定的補償穩定因子。

圖9 反演結果誤差與穩定因子中r 的關系曲線

以圖9 的誤差曲線推算出最優的補償穩定因子為r=10-2，故進一步測試其對二維BP 鹽丘模型(圖10)的適用性。該模型網格數為625×120，空間采樣間隔為10 m，其Q值同樣由速度模型根據經驗公式生成。模型的空間采樣間隔為10 m。初始模型(圖10c)保留了頂部60 m以上的準確速度，其余由準確速度模型平滑獲得。以100 m間隔在模型表面布設60炮，首個炮點距模型左邊界170 m；接收方式為雙邊可變觀測系統，為接收大入射角的反射信息，設置最大炮檢距為1200 m，道間距為10 m。選用主頻為3、5 Hz的Ricker子波進行正、反演，參考角頻率為200 rad/s，記錄的樣點數為2001，時間間隔為0.8 ms。主頻為3 Hz子波的正演數據的反演結果作為主頻為5 Hz子波的正演數據反演的初始模型，分別反演15和10次。對比反演結果(圖11、圖12)，梯度Q補償可以有效加速反演收斂。即便選用較低頻的數據，經過梯度補償后，FWI結果的鹽丘速度精度有所提高，斷層附近層間微弱的變化得到凸顯，對于鹽丘下方高衰減區域的表征也更加清晰，證實了本文方法對不同模型的適應性。

圖10 BP 鹽丘模型

圖11 BP 鹽丘模型主頻為3 Hz 子波正演數據的反演結果

圖12 BP 鹽丘模型主頻為5 Hz 子波正演數據的反演結果

3 結論

本文基于常分數階拉普拉斯算子黏滯聲波方程提出了一種基于Q補償的梯度預處理方法，能加速黏滯聲波全波形反演，改善反演效果。該梯度通過Q補償后的反傳波場與正傳波場互相關獲得。針對特定常分數階拉普拉斯算子黏滯波動方程的反傳波場Q補償，本文提出了一種有效的穩定補償策略，可以確保計算過程的穩定，且提升反演效果。得到如下結論。

(1)基于常分數階拉普拉斯算子黏滯聲波方程，通過在時間—波數域的解析解中引入穩定因子，以確保Q補償反傳波場計算的穩定，保證了梯度計算過程中的數值穩定。數值算例表明，未經穩定因子作用的直接補償會使反傳波場不穩定，所提方法能有效解決補償過程中因高頻成分指數放大而導致的數值不穩定問題。

(2)BP 氣囪模型和鹽丘模型全波形反演算例表明，實施所提出的穩定的Q補償梯度預處理方法能為黏滯聲波全波形反演提供更有利于收斂的目標函數梯度，相比于未補償的黏滯聲波全波形反演的梯度，補償的梯度能量整體有所提升，不僅能提高反演的收斂速度，還能提供更精確的反演結果，尤其是在高衰減區域。

(3)在本文提出的Q補償全波形反演中，穩定因子需要選在一定范圍內才有提升效果。在滿足波場穩定延拓的前提下，應該依據所反演區域的黏滯性(τcmax)盡量選擇合適的穩定因子以確保相位信息準確。本文試算結果顯示，穩定因子可選用σ2=10-2exp(-τcmaxΔt)。