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

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

0 引言

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

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

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

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

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

1 方法原理

1.1 常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程

Zhu等[26]提出的解耦分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程為

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

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

1.2 基于伴隨狀態(tài)法的目標(biāo)函數(shù)梯度計(jì)算

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

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

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

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

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

1.3 基于穩(wěn)定化策略的穩(wěn)定衰減補(bǔ)償方法

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

由上式可知，在補(bǔ)償過程中地震波場的振幅呈現(xiàn)指數(shù)型增長，容易導(dǎo)致數(shù)值不穩(wěn)定。為確保波場的穩(wěn)定補(bǔ)償，本文提出一種含穩(wěn)定因子的Q補(bǔ)償方法。引入穩(wěn)定因子σ的波場補(bǔ)償解析表達(dá)式為

當(dāng)σ取0時，式(10)退化為式(9)，穩(wěn)定化程度隨著σ的增大而增強(qiáng)。利用三角函數(shù)的周期性，式(10)可等價為

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

根據(jù)式（12），可以獲得C1和C2的表達(dá)式

C2可進(jìn)一步簡化為

式中

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

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

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

2 數(shù)值模擬算例

2.1 穩(wěn)定衰減補(bǔ)償策略的有效性測試

為了驗(yàn)證本文衰減補(bǔ)償方法的有效性，選用具有強(qiáng)衰減區(qū)的二維BP氣囪模型(圖1)進(jìn)行測試。該模型網(wǎng)格數(shù)為398×161，空間采樣間隔為10 m。Q模型由速度模型根據(jù)經(jīng)驗(yàn)公式Q=3.516c2.2生成，其中速度c的單位為km/s。該模型在上部的氣囪區(qū)吸收非常強(qiáng)，當(dāng)?shù)卣鸩ù┻^該區(qū)域時會快速衰減，導(dǎo)致下層的反射信號非常弱，應(yīng)用常規(guī)方法在氣囪區(qū)下方無法獲得可靠的反演結(jié)果。

圖1 BP 氣囪模型

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

選用σ2=10-2exp(-τcmaxΔt)為補(bǔ)償時的穩(wěn)定因子，比較補(bǔ)償前、后的反傳波場及其表征補(bǔ)償效果。圖2 為第10 炮不同時刻補(bǔ)償前、后反傳波場快照對比，可見未經(jīng)衰減補(bǔ)償?shù)牟▓稣穹鄬^弱，而經(jīng)補(bǔ)償后恢復(fù)了振幅信息，在高衰減區(qū)域明顯增強(qiáng)，表明該衰減補(bǔ)償算法能實(shí)現(xiàn)有效的振幅補(bǔ)償，且補(bǔ)償后的反傳波場沒有不穩(wěn)定現(xiàn)象。

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

2.2 基于穩(wěn)定梯度衰減補(bǔ)償?shù)姆囱轀y試

為了驗(yàn)證穩(wěn)定Q補(bǔ)償對黏滯聲波FWI的有效性，利用BP氣囪模型進(jìn)行反演測試。圖3為第一次迭代的梯度。由梯度趨勢可以看出，Q補(bǔ)償后，梯度的整體趨勢大體不變，但模型中淺層尤其是高衰減區(qū)域的梯度能量明顯增強(qiáng)，反映的結(jié)構(gòu)信息也更清晰。

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

采用炮點(diǎn)平方預(yù)處理方式，P-L-BFGS法迭代20、40 和100 次反演結(jié)果如圖4所示?？梢钥闯?，基于梯度Q補(bǔ)償?shù)娜ㄐ畏囱菰诔跏茧A段的迭代就可以恢復(fù)大體的速度模型，而利用未補(bǔ)償?shù)奶荻葻o法在反演初期獲得大體的速度趨勢。雖然隨著迭代次數(shù)的增加，梯度未補(bǔ)償黏滯聲波全波形反演效果也得到改善，但即便迭代至100 次，仍然未能反演出準(zhǔn)確的速度。對比x=1000、1500 和3200 m 的最終反演的速度縱向曲線(圖5)可以看出，基于梯度Q補(bǔ)償?shù)娜ㄐ畏囱莶粌H能獲得更精確的層速度，而且較清晰地揭示了深部的大傾角傾斜地層，驗(yàn)證了本文方法在加快黏滯聲波全波形反演的良好效果。

圖4 BP 氣囪模型不補(bǔ)償（左）與補(bǔ)償（右）不同迭代次數(shù)的反演結(jié)果對比

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

2.3 不同補(bǔ)償穩(wěn)定因子選取測試

穩(wěn)定因子σ2=rexp(-τcmaxΔt)中r分別取不同值時，第10炮0.4s 時刻的反傳波場如圖6所示。對比圖2c與圖6 可以看出，選用的穩(wěn)定因子都可以使補(bǔ)償?shù)姆磦鞑▓霰３址€(wěn)定，同時補(bǔ)償后的反傳波場均較未補(bǔ)償?shù)姆磦鞑▓瞿芰坑忻黠@增強(qiáng)，大穩(wěn)定因子比小穩(wěn)定因子補(bǔ)償后的波場能量弱。

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

圖7展示了由圖6所示的反傳波場計(jì)算的第一次迭代的梯度。對比圖7 與圖3 可見：與小穩(wěn)定因子相比，大穩(wěn)定因子補(bǔ)償?shù)奶荻葦?shù)值較小，但依然大于未補(bǔ)償波場計(jì)算的梯度。

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

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

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

應(yīng)用平均相對誤差定量評價反演結(jié)果的精度

圖9為不同r迭代17次和66次的反演結(jié)果誤差對比。從圖9a 可以看到，當(dāng)r=10-1時，補(bǔ)償穩(wěn)定因子過大，抑制了波場振幅，獲得的反演結(jié)果偏差大。相對于傳統(tǒng)未補(bǔ)償?shù)姆囱萁Y(jié)果，使用較小穩(wěn)定因子的梯度補(bǔ)償方法的反演結(jié)果的誤差更小。由圖9b可以看出，當(dāng)補(bǔ)償穩(wěn)定因子過小時，如r=10-32，即便最終的反演誤差低于傳統(tǒng)未補(bǔ)償?shù)恼`差，但由于未能保證補(bǔ)償過程的穩(wěn)定，出現(xiàn)明顯的數(shù)值不穩(wěn)定現(xiàn)象，影響反演結(jié)果。對不同穩(wěn)定因子的反演結(jié)果進(jìn)行誤差分析可以看出，隨著穩(wěn)定因子的增大，反演誤差呈現(xiàn)出先降后升的現(xiàn)象。在全波形反演過程中進(jìn)行梯度補(bǔ)償時，應(yīng)選用不抑制補(bǔ)償、同時保證數(shù)值計(jì)算穩(wěn)定的補(bǔ)償穩(wěn)定因子。

圖9 反演結(jié)果誤差與穩(wěn)定因子中r 的關(guān)系曲線

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

圖10 BP 鹽丘模型

圖11 BP 鹽丘模型主頻為3 Hz 子波正演數(shù)據(jù)的反演結(jié)果

圖12 BP 鹽丘模型主頻為5 Hz 子波正演數(shù)據(jù)的反演結(jié)果

3 結(jié)論

本文基于常分?jǐn)?shù)階拉普拉斯算子黏滯聲波方程提出了一種基于Q補(bǔ)償?shù)奶荻阮A(yù)處理方法，能加速黏滯聲波全波形反演，改善反演效果。該梯度通過Q補(bǔ)償后的反傳波場與正傳波場互相關(guān)獲得。針對特定常分?jǐn)?shù)階拉普拉斯算子黏滯波動方程的反傳波場Q補(bǔ)償，本文提出了一種有效的穩(wěn)定補(bǔ)償策略，可以確保計(jì)算過程的穩(wěn)定，且提升反演效果。得到如下結(jié)論。

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

(2)BP 氣囪模型和鹽丘模型全波形反演算例表明，實(shí)施所提出的穩(wěn)定的Q補(bǔ)償梯度預(yù)處理方法能為黏滯聲波全波形反演提供更有利于收斂的目標(biāo)函數(shù)梯度，相比于未補(bǔ)償?shù)酿暡ㄈㄐ畏囱莸奶荻?，補(bǔ)償?shù)奶荻饶芰空w有所提升，不僅能提高反演的收斂速度，還能提供更精確的反演結(jié)果，尤其是在高衰減區(qū)域。

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