基于Low-rank一步法波場延拓的黏聲各向異性介質純qP波正演模擬

顧漢明 張奎濤 劉春成 王建花

（①中國地質大學（武漢）地球物理與空間信息學院，湖北武漢430074；②中海油研究總院有限責任公司，北京100028）

0 引言

由于實際地下介質充滿流體、裂縫及裂隙，因而同時表現出各向異性與黏彈性特征，而這種黏彈各向異性的現象也在實驗室和波場測量中被觀測到［1］，因此衰減和各向異性是地震波場數值模擬中越來越不容忽略的因素。

在進行彈性波數值模擬時，得到的混合波場通常同時含有縱波和橫波，而在進行地震偏移成像，尤其是彈性波逆時偏移成像時，須將縱、橫波解耦，以得到物理意義明確的成像剖面［2］。為此，Alkhalifah［3-4］提出了聲學近似假設下的擬聲波方程，人為設定沿對稱軸方向的橫波速度為零，推導出標量形式的四階微分方程。為了簡化方程，提高計算效率，Du等［5］和Zhou等［6］通過不同的降階方法，推導出不同形式的二階耦合qP波微分方程。Duveneck等［7］從胡克定律和運動方程出發，推導得到另一種形式的擬聲波方程，且其波場具有明確的物理意義。程玖兵等［8-9］推導得到各向異性介質的偽純模式波動方程，使其在運動學上同彈性波動方程等價。

考慮到衰減的影響，Zhang等［10］推導出時域擬微分算子的黏聲各向同性方程，并成功地應用于逆時偏移中；隨后，Suh等［11］將該方法拓展到各向異性介質中，補償了各向異性介質中衰減的影響；Xu等［12］基于二階擬微分方程，實現了黏聲各向異性介質的正演模擬及逆時偏移；李金麗等［13］使用有限差分法進行黏彈VTI介質的正演模擬，并實現了該類介質的雙程波照明分析方法。

上述各方法均存在擬聲波方程的突出缺點，即當模型參數中有變化較大的傾角和方位角時，會出現不穩定現象［14］、殘留橫波假象［15］，或只有在滿足Thomsen參數ε≥δ時正演模擬才是穩定的［16］。

為了解決此類問題，有學者另辟蹊徑。從qP波與qSV波的頻散關系出發，直接解耦以構建純qP波控制方程；要實現qP與qSV波的完全分離，關鍵在于導出簡單、靈活的純qP波控制方程和設計有效的擬微分求解算法。例如，Du等［5］和Zhang等［17］基于弱各向異性近似和平方根近似，推導出一種純qP波解耦方程，其表現為時間—波數域形式，能較準確地描述TTI介質中的運動學特征；Zhan等［18］推導了一種新的解耦方程，將P波與qSV波完全分離，并成功將其應用于逆時偏移中；Xu等［19］提出了一種純qP波橢圓微分方程，通過將微分方程分解成一個標量算子和一個微分算子，并用橢圓分解方法修正標量算子，達到校正振幅誤差的目的；楊鵬等［20］改進了Xu等［19］所提方程，使其振幅更準確、更均衡；胡書華等［21］則將Xu等［19］二階微分方程做降階處理，通過Lebedev交錯網格，得到TTI介質下純qP波的穩定傳播的波場；張慶朝等［22］利用近似相速度公式，導出TTI介質三維qP波頻散方程及波動方程，并使用偽譜法做正演模擬。

時間遞歸積分延拓［23］（Recursive integral timeextrapolation，RITE）是一類具有較高穩定性、可實現大時間步長條件下正演模擬的方法。在不同的RITE方法中，Low-rank波場延拓因其高效率和靈活的近似精度控制而尤為突顯。Fomel等［24］首次使用Low-rank兩步法波場延拓，實現了各向同性介質及各向異性介質的正演模擬，并指出該分解算法能準確描述地震波場的動力學特征；隨后，Song等［25］將Low-rank應用于正交各向異性介質的正演模擬；Fang等［26］將交錯網格有限差分與Low-rank相結合，實現了基于交錯網格的Low-rank有限差分波場正演模擬，進一步提高了計算效率和計算精度；Sun等［27］結合Low-rank分解算法和一步法波場延拓，實現了黏彈性介質的Low-rank一步法逆時偏移。在中國國內，該類方法應用較少。黃金強等［28-29］分別使用Low-rank兩步法及Low-rank有限差分法，進行了TI介質純qP波波場的高效正演模擬；袁雨欣等［30］則使用Low-rank分解及Lowrank有限差分法求解二階解耦的彈性波方程，有效提高了計算精度。

本文基于上述研究成果，引入一步法波場延拓方法，推導了黏聲介質方程在空間—波數域的表達形式，并結合空間—波數域各向異性介質延拓算子，構建了一種適用于黏聲各向異性介質的空間—波數域純qP波波場延拓算子，該算子消除了偽橫波干擾，不受模型參數限制，且地震波場能穩定傳播；為兼顧計算精度與計算效率，引入Low-rank分解算法進行求解，實現了基于Low-rank一步法波場延拓的黏聲各向異性介質純qP波正演模擬。通過簡單和復雜模型測試，驗證了本文方法的正確性及對復雜介質的適用性，為基于Q補償的各向異性介質逆時偏移提供了理論依據。

1 方法原理

1.1 一步法波場延拓

根據一步法波場延拓理論［31］，下一時刻的波場延拓可由空間—波數混合域算子近似表示［32］，即

式中：t為時間坐標；x為空間坐標；k為波數坐標；Φ為角頻率函數；Δt為時間步長；（k，t）為p（x，t）的傅里葉變換，即

將式（1）中的Δt替換成－Δt，可得

結合式（1）與式（3），顯然有

在上述推導中，式（1）稱為一步法波場延拓，式（4）稱為兩步法波場延拓。通常，一步法延拓比兩步法延拓更穩定，但一步法實現較復雜，涉及復數運算；兩步法延拓僅為實數運算，其實現過程則相對簡捷。

對Φ做如下高頻近似［33］，即

式中v為地震波場中的相速度。在各向同性介質中，v僅為空間坐標的函數。而在各向異性介質中，地震波的傳播速度隨傳播方向變化，故此時相速度不僅為空間坐標的函數，還是波數矢量k的函數。

進一步地，忽略式（5）中的高階項，有

將式（6）代入式（1）或式（4），對應可得

1.2 黏聲各向異性介質純qP波波場延拓算子

從式（1）可知，角頻率函數Φ控制著波的傳播，因此可通過構建不同的Φ得到不同介質中的波的傳播。因此，為了構建本文的黏聲各向異性介質純qP波波場延拓算子Φ，引入如下VTI介質純qP波聲學近似下的頻散方程［4］，即

式中：vv為沿對稱軸方向的qP波相速度；為各向同性面內的qP波相速度；為橢圓各向異性參數，其中ε、δ為Thomsen參數。該式可準確地描述qP波的運動學特征，且不受模型參數的限制，不存在偽橫波干擾。

結合式（6）與式（9），可得

該式即為由角頻率函數Φ控制的VTI介質純qP波波場延拓算子。

另一方面，考察如下二維黏聲方程［34］

式中：vx、vz分別為x、z方向的速度場；p為壓力場；e為記憶變量；ρ為密度；vP為縱波速度；τσ、τε分別為應力和應變的松弛時間，其表達式為

式中：Q為品質因子；ω為角頻率。

將式（11）變換到頻率—波數域，即有

由式（13）前兩項可得

將式（14）代入式（13）后兩項，并消去，得到

當Q不是特別小的情況下，τ?1，即有

因此，式（16）可簡化為

即黏聲介質的角頻率函數在空間—波數域表達為

該式即為由Φ控制的黏聲介質波場延拓算子。可見，在空間—波數混合域，由非黏彈性轉變成黏彈性，只需乘以即可，即黏彈性介質在空間—波數混合域的表達為復數形式。

結合式（10）與式（19），得到由Φ控制的黏聲VTI介質純qP波波場延拓算子

更進一步地，對于黏聲TTI介質純qP波波場延拓算子，可做如下坐標變換

式中θ為極化角，便有

因此，通過式（1）和式（20）或式（22），可實現聲各向異性介質純qP波波場的穩定延拓。

為了便于后文模型測試對比，此處給出常規二維黏聲TTI介質擬聲波方程［12］，即

式中：ph、pv分別為壓力場水平分量和垂直分量；Hx、Hz為坐標變換算子，其表達式為

式（23）即為常規二維黏聲TTI介質擬聲波波動方程。因該式令對稱軸方向的橫波速度為零，會產生菱形偽橫波干擾，其本質原因在于對qP波相速度公式中的根式項做了平方處理，進而引入了增根，并在傾角變化劇烈的強各向異性區域極易出現數值模擬不穩定現象，從而限制了其廣泛應用。

1.3 Low-rank分解算法

由式（1）可知，在時間方向每延拓一步，需做一次傅里葉變換和Nxz次傅里葉反變換（Nxz＝NxNz，Nx、Nz分別為模型x、z方向上的網格點數），其計算復雜度為，這對于大規模復雜模型地震波波場延拓來說，計算成本極其昂貴。為了降低計算成本，本文引入Low-rank分解算法［35］，只需幾次（通常為2或3次）傅里葉反變換，計算復雜度為mNxzO（lgNxz）（m為很小的實數），從而顯著降低了計算成本。

Low-rank基本思想是對空間—波數域的傳播矩陣W進行分解，即

式中W（x，k）刻畫了波的所有傳播特征，包括了動力學與運動學特征，且式（25a）為兩步法延拓波場矩陣，只涉及實數運算，式（25b）為一步法波場延拓矩陣，還涉及到復數運算。

當Δt固定時，有如下分解

式中：W1、W2為矩陣W分解后的矩陣，且W1表示矩陣W中M個線性無關的列向量組成的最大空間體積，表征空間位置x，W2表示矩陣W中N個線性無關的行向量組成的最大空間體積，表征波數位置k；A＝｛amn｝為連接矩陣W1和W2的權重系數。特別地，當模型給定時，通過設定容許誤差和Δt，即可求取上述子矩陣和權系數及對應子矩陣的秩M和N，且在整個波場延拓過程中，上式分解過程只需做一次分解即可，適用于并行計算。

Low-rank分解具體實施步驟如下。

（1）求取子矩陣W1：

1）從全矩陣W中隨機選取βR行（通常β取3或4，R為全矩陣W的秩）組成一個新矩陣S1；

2）對S1進行旋轉QR分解，即

式中Π1為包含矩陣S1各列置換信息的置換矩陣；

3）取Π1中前R個置換序號組成序列Λ1＝｛k1，k2，…，kR｝；

4）取全矩陣W中對應序列Λ1的列向量，組成的矩陣W1即為所求

（2）求取子矩陣W2：

1）從全矩陣W中隨機選取βR列組成一個新矩陣S2，其轉置為；

式中Π2為包含矩陣各列置換信息的置換矩陣。

3）取Π2中前R個置換序號組成序列Λ2＝｛x1，x2，…，xR｝；

4）取全矩陣W中對應序列Λ2的行向量，組成的矩陣W2即為所求

（3）求取最佳近似分離權系數A：

選取子矩陣W1和子矩陣W2相互交叉的部分求偽逆，即

（4）結合以上步驟就可得矩陣的Low-rank分解

經過上述步驟Low-rank分解后，可將式（26）分別代入式（8）和式（7）中，就有

式（32）為基于Low-rank分解的兩步法波場延拓公式，式（33）為基于Low-rank分解的一步法波場延拓公式。

2 模型試算

2.1 均勻模型

為了驗證本文構建的純qP波算子不受模型參數的限制，且無偽橫波干擾，設計了尺寸為2000m×2000m的均勻介質模型。其空間網格單元為5m×5m，時間步長為1ms，總時間采樣點數為1000；縱波震源位于模型中心，采用主頻為25 Hz的Ricker子波震源；邊界條件為海綿邊界條件；接收排列范圍是0～2000m，排列深度為750m，道間距為10m，共201道接收。均勻黏聲各向異性純qP波介質彈性參數見表1。

表1 均勻黏聲各向異性介質純qP波彈性參數

圖1上部為黏聲TTI介質常規擬聲波方程分別在介質類型Ⅲ（圖1a）、Ⅳ（圖1b）、Ⅴ（圖1c）時400ms時刻波場快照，圖1下部為黏聲TTI介質純qP波方程分別在介質類型Ⅲ（圖1d）、Ⅳ（圖1e）、Ⅴ（圖1f）時400ms時刻波場快照。從圖中可見常規擬聲波方程（即直接將橫波速度設為0）在ε＞δ時雖能穩定傳播，但會出現菱形假象，嚴重干擾了有效波；在ε＜δ時會出現不穩定現象；只有在ε＝δ時才能穩定傳播且無假象。這樣就極大地限制了其應用。而本文構建純qP波方程（圖1下部），無論是ε＜δ、ε＞δ還是ε＝δ時，都能穩定地傳播且不受橫波干擾，表明本文方法突破了該限制，具有更強適用性。

圖2為純qP波方程分別在介質類型Ⅰ、Ⅱ、Ⅲ時400ms時刻波場快照。對比對應為黏聲VTI（圖2a）、HTI（圖2b）和TTI（圖2c）三種介質發現，通過引入坐標旋轉，使得TTI介質具有更強適用性，更能反映實際地下真實介質的情況。

采用表2所示介質參數做正演模擬，研判本文構建的純qP波算子是否同時具有各向異性特征及黏滯性特征，并分別取四種介質的四分之一波場拼成一張波場快照以突顯對比效果（圖3）。從圖中對比發現，橫向對比可體現出各向異性特征，即聲波波前為一個圓，在VTI介質中波前為橢圓；縱向對比可體現出黏滯性特征，即波前走時差異及能量差異。為了更直觀地分析，采用表1介質類型Ⅲ中的參數，并設定不同的Q值（無窮大、200、50和20），得到不同Q值時的地震記錄。

圖1 均勻黏聲各向異性介質擬聲波方程（上）和純qP波方程（下）400ms時刻波場快照

圖2 純qP波方程在均勻黏聲各向異性介質中400ms時刻波場快照

表2 不同介質類型彈性參數

從不同Q值下的單道波形記錄及頻譜（圖4）對比可知，由于地下介質的黏滯性，使得地震波的振幅得到明顯的衰減，品質因子Q越低，其振幅衰減的越多，波形記錄的相位延遲越嚴重，波形變寬，則畸變也越嚴重。從頻譜圖看到，品質因子Q越低，其能量越低，子波主頻也越低，符合真實實際地下介質中地震波由于地層的吸收衰減能量逐漸減弱，高頻衰減快，主頻降低的特征，說明了本文構建的純qP波算子的正確性，即能同時體現出各向異性特征及黏滯性特征。

圖3 不同介質類型在400ms時刻的波場快照對比

圖4 不同Q值下的單道波形記錄（a）及頻譜（b）對比

2.2 簡單TTI模型

選用Overthrust模型測試Low-rank分解算法的穩定性。該模型構造較簡單，有利于觀測波的傳播規律，特別是運動學特性；此外，模型中部區域各向異性特征明顯，且傾角變化劇烈，是測試分析算法穩定性的標準模型（圖5）。該模型尺寸為7000m×2000m，空間網格單元為10m×10m，時間步長為2ms，總時間采樣點數為1250，記錄長度為2.5s，縱波震源位置為（3500m，10m），震源是主頻為30Hz的Ricker子波，選取海綿邊界條件，接收排列范圍為0～7000m，排列深度為0，道間距為20m，共351道接收。模型參數如圖5所示，其中品質因子由經驗公式［36］Q＝14v2.2得到。

圖5 二維黏聲TTI Overthrust模型參數

首先驗證Low-rank分解算法重構矩陣的有效性及正確性。選取圖5模型中心處（vP＝3000m/s，θ＝50°，ε＝0.15，δ＝0.08）的真實矩陣，即得到該真實矩陣的實部（圖6a）、虛部（圖6d）和模（圖6g），Low-rank近似矩陣的實部（圖6b）、虛部（圖6e）和模（圖6h），以及真實矩陣與Low-rank近似矩陣的實部（圖6c）、虛部（圖6f）和模（圖6i）的誤差。對比可知Low-rank近似矩陣接近于真實矩陣；此外，其誤差的數量級僅為10－7，充分表明Low-rank分解能精確地重構真實矩陣。

圖6 Overthrust模型中心處的傳播矩陣W（x，k）＝eiΦ（x，k）Δt

然后探討基于Low-rank分解算法的波場延拓的穩定性。應用一步法波場延拓得到不同子波主頻、不同時間步長的600ms時刻波場快照（圖7）。對比子波主頻為30Hz、時間步長為2ms（圖7a）與子波主頻為60 Hz、相同時間步長（圖7b）的波場可知，提高地震子波主頻使子波波長變短，空間子波的分辨率明顯提高，其地震波場特征更清晰，也說明了子波主頻與波長呈反比關系；若使用有限差分正演，一個波長內約有兩個采樣點，必然會嚴重干擾有效波，并形成顯著的頻散，而Low-rank分解算法能在不減小網格單元尺寸的情況下適應高頻正演且不損失精度，證明Low-rank分解算法是高精度波場延拓算法。同樣地，將子波主頻為60Hz、時間步長為4ms（圖7c）的波場與圖7b對比，可知在相同高頻子波主頻下，Low-rank算法能通過提高時間步長提升計算效率且不損失計算精度。據穩定性條件，有限差分法時間步長不大于1.5ms，說明Low-rank分解算法比有限差分算法能更有效地通過提高地震子波主頻和時間步長改善地震記錄分辨率和計算效率，且不損失計算精度，表明Low-rank分解算法能同時兼顧計算精度和計算效率。

圖7 不同參數下的Low-rank一步法600ms時刻波場快照

一般而言，Low-rank分解算法的計算效率與秩R有很大關系。R值大小決定正反傅里葉變換的次數，即R越大，正、反傅里葉變換次數越多，計算量越大，導致計算效率降低。

表3為不同時間步長和不同容許誤差下的秩R。從表中可知，當容許誤差一定時，隨時間步長增大，R逐漸增加。值得注意的是，雖然隨著時間步長的增大，R也增大，其計算效率會降低，但增大時間步長所提高的計算效率足以彌補R增加導致的計算效率降低，即增大時間步長還是會在總體上提高計算效率；當時間步長一定時，R隨著容許誤差的降低而增大。特別地，時間步長的選取一般通過試驗得到，即在滿足精度及穩定性要求的前提下，試驗選取最大時間步長。通常，在設定容許誤差并滿足計算精度要求的前提下，以R不超過12時的時間步長為宜。

表3 不同時間步長和容許誤差下的秩R

再次設計尺寸為2000m×2000m的兩層均勻介質模型，空間網格單元為5m×5m，縱波震源置于（1000m，900m）處，震源選取主頻為30Hz的Ricker子波。為了便于對比，采用常密度聲學介質，上、下層厚度均為1000m，上、下層縱波速度分別為1700、2300m/s。

圖8為有限差分法、Low-rank兩步法及Lowrank一步法在不同時間步長下得到的400ms時刻波場快照。對比可知：當時間步長為1ms時，三種方法的波場都能穩定地傳播；當步長為2ms時，有限差分法出現不穩定，難以進行波場延拓；當步長增至5ms時，Low-rank兩步法出現頻散，達到能正演的極限，而此時Low-rank一步法所得波場仍能穩定地傳播，并保持很高計算精度；當步長高達15ms時，Low-rank一步法出現頻散，計算精度開始下降，而其他兩種方法根本無法進行波場延拓。該對比結果表明：Low-rank一步法的穩定性最高，其次是Low-rank兩步法，有限差分法則是最低。再次證實Low-rank算法是一種能兼顧計算精度和計算效率的較新穎算法。

若針對上述雙層均勻介質模型進行計算效率分析，會因模型網格點數較少，單炮計算時間很短，而不便于比較，故取10炮CPU耗時進行計算效率分析、對比，其中單炮記錄時長為2s。

從上述三種正演模擬法在不同時間步長時的CPU耗時（表4）可見：當時間步長為1ms時，有限差分法計算效率最高，其次為Low-rank兩步法，而Low-rank一步法因要做復數運算，使得其計算效率最低。另一方面，由于Low-rank兩步法和一步法的穩定性比有限差分法高較多，因此可使用較大時間步長。結合圖8及表4可知，在同等計算精度條件下，Low-rank兩步法和一步法因可使用較大時間步長，使其計算效率高于有限差分法，且隨著時間步長的增加，無論是Low-rank兩步法還是一步法，其計算效率均降低，即說明在滿足計算精度前提下，增加時間步長，能提高計算效率。

表4 不同方法在不同時間步長時的CPU耗時

2.3 復雜Hess模型

采用二維Hess黏聲TTI介質模型（圖9）測試Low-rank一步法純qP波波場延拓對復雜介質的適應性。該模型尺寸為3000m×2500m，空間網格單元為5m×5m，時間步長為5ms，總時間采樣點為500，縱波震源置于（2000m，10m）處，仍采用前述震源子波和邊界條件，接收排列范圍是0～3000m，排列深度為0，接收點間距為10m，共301道接收。模型參數見圖9，品質因子同樣可由經驗公式［36］Q＝14v2.2得到，極化角θ為45°。

圖8 三種方法在不同時間步長下得到的400ms時刻波場快照

從黏聲TTI介質純qP波800ms時刻波場快照（圖10a）可見，波場中只含有純P波，無橫波干擾，且波場特征復雜；同時，地震記錄（圖10b）上表現出淺部能量強、深部能量弱的特征，與實際采集的地震記錄特征相同。這為地下波場特征研究及該類介質的逆時偏移提供了更多理論依據，也表明本文構建的純qP波方程能適應復雜介質和復雜構造，能夠穩定地對其進行波場延拓。

對比不同時間步長時所得地震記錄及其差值（圖11），可見10ms步長地震記錄（圖11a）的淺部反射同相軸不連續，出現輕微頻散，說明5ms步長地震記錄的精度高于10ms步長，即隨著時間步長的增加，計算精度降低。因此，在可接受計算精度前提下，采用較大時間步長，可顯著提高計算效率。

2.4 復雜BP模型

為檢驗所提方法對傾角劇變復雜介質的穩定性，截取一段傾角變化范圍是－55°～55°的二維BP2007黏聲TTI介質模型（圖12）進行測試。該模型規模為17km×11.25km，空間網格單元為6.25m×6.25m，時間步長為5ms，總時間采樣點為1600，記錄長度為8s，縱波震源位于（68.5km，0），采用主頻為20Hz的Ricker子波震源，仍為上述邊界條件，接收排列范圍是0～17km，排列深度為0，接收點間距為25m，共681道接收。模型參數見圖12，品質因子仍沿用上述方式得到。

圖9 二維Hess黏聲TTI介質模型參數

圖10 二維Hess黏聲TTI介質800ms時刻波場快照（a）及地震記錄（b）

圖11 時間步長為5ms（a）和10ms（b）時的地震記錄及其差值（c）

圖13為采用二維BP2007模型得到的黏聲TTI介質純qP波4s時刻波場快照及其地震記錄。從波場快照（圖13a）可見，地震波場能穩定地傳播，說明所提方法能適應傾角劇烈變化的情形，具有較高穩定性；同時，在地震記錄（圖13b）中，地震波能量表現出淺部強、深部弱的特征，體現了地下介質的黏滯性特征，即說明本文方法能在兼顧計算精度的同時，使用較大時間步長，以提高計算效率，并能夠適用于復雜介質、復雜構造，包括傾角劇變情形，是一種穩定高效的正演模擬方法。

圖12 二維BP2007黏聲TTI介質模型參數

圖13 二維BP2007黏聲TTI介質4s時刻波場快照（a）及地震記錄（b）

3 結論與分析

本文在前人研究基礎上，引入一步法波場延拓方法，構建了一種適用于黏聲各向異性介質的空間—波數域純qP波波場延拓算子，該算子消除了偽橫波干擾，不受模型參數限制，且地震波場能穩定傳播；為兼顧計算精度與計算效率，引入Low-rank分解算法進行求解，實現了基于Low-rank一步法波場延拓的黏聲各向異性介質純qP波正演模擬。通過多種模型的測試，得出以下認識和結論：

（1）通過本文構建的純qP波算子計算得到的波場，能同時表現出各向異性特征和黏滯性特征，更符合實際地下介質情況，也證明所構建純qP波算子的正確性；

（2）本文方法克服了擬聲波方程的局限性，消除了偽橫波干擾，不受模型參數限制且地震波場能穩定地傳播，具有較為廣泛的應用前景；

（3）本文方法可適應較大時間步長，無數值頻散，能同時兼顧計算效率和計算精度，是一種穩定、高效的正演模擬方法。