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

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

0 引言

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

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

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

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

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

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

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

1 方法原理

1.1 一步法波場延拓

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

由式（13）前兩項可得

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

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

因此，式（16）可簡化為

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

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

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

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

式中θ為極化角，便有

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

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

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

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

1.3 Low-rank分解算法

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

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

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

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

式中：W1、W2為矩陣W分解后的矩陣，且W1表示矩陣W中M個線性無關的列向量組成的最大空間體積，表征空間位置x，W2表示矩陣W中N個線性無關的行向量組成的最大空間體積，表征波數(shù)位置k；A＝｛amn｝為連接矩陣W1和W2的權重系數(shù)。特別地，當模型給定時，通過設定容許誤差和Δt，即可求取上述子矩陣和權系數(shù)及對應子矩陣的秩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）求取最佳近似分離權系數(shù)A：

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

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

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

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

2 模型試算

2.1 均勻模型

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

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

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

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

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

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

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

表2 不同介質類型彈性參數(shù)

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

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

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

2.2 簡單TTI模型

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

圖5 二維黏聲TTI Overthrust模型參數(shù)

首先驗證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近似矩陣接近于真實矩陣；此外，其誤差的數(shù)量級僅為10－7，充分表明Low-rank分解能精確地重構真實矩陣。

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

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

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

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

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

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

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

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

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

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

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

2.3 復雜Hess模型

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

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

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

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

2.4 復雜BP模型

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

圖9 二維Hess黏聲TTI介質模型參數(shù)

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

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

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

圖12 二維BP2007黏聲TTI介質模型參數(shù)

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

3 結論與分析

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

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

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

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