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

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

0 引言

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

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

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

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

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

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

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

1 方法原理

1.1 一步法波場(chǎng)延拓

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

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

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

結(jié)合式（1）與式（3），顯然有

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

對(duì)Φ做如下高頻近似［33］，即

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

進(jìn)一步地，忽略式（5）中的高階項(xiàng)，有

將式（6）代入式（1）或式（4），對(duì)應(yīng)可得

1.2 黏聲各向異性介質(zhì)純qP波波場(chǎng)延拓算子

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

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

結(jié)合式（6）與式（9），可得

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

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

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

式中：Q為品質(zhì)因子；ω為角頻率。

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

由式（13）前兩項(xiàng)可得

將式（14）代入式（13）后兩項(xiàng)，并消去，得到

當(dāng)Q不是特別小的情況下，τ?1，即有

因此，式（16）可簡化為

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

該式即為由Φ控制的黏聲介質(zhì)波場(chǎng)延拓算子?？梢姡诳臻g—波數(shù)混合域，由非黏彈性轉(zhuǎn)變成黏彈性，只需乘以即可，即黏彈性介質(zhì)在空間—波數(shù)混合域的表達(dá)為復(fù)數(shù)形式。

結(jié)合式（10）與式（19），得到由Φ控制的黏聲VTI介質(zhì)純qP波波場(chǎng)延拓算子

更進(jìn)一步地，對(duì)于黏聲TTI介質(zhì)純qP波波場(chǎng)延拓算子，可做如下坐標(biāo)變換

式中θ為極化角，便有

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

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

式中：ph、pv分別為壓力場(chǎng)水平分量和垂直分量；Hx、Hz為坐標(biāo)變換算子，其表達(dá)式為

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

1.3 Low-rank分解算法

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

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

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

當(dāng)Δt固定時(shí)，有如下分解

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

Low-rank分解具體實(shí)施步驟如下。

（1）求取子矩陣W1：

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

2）對(duì)S1進(jìn)行旋轉(zhuǎn)QR分解，即

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

3）取Π1中前R個(gè)置換序號(hào)組成序列Λ1＝｛k1，k2，…，kR｝；

4）取全矩陣W中對(duì)應(yīng)序列Λ1的列向量，組成的矩陣W1即為所求

（2）求取子矩陣W2：

1）從全矩陣W中隨機(jī)選取βR列組成一個(gè)新矩陣S2，其轉(zhuǎn)置為；

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

3）取Π2中前R個(gè)置換序號(hào)組成序列Λ2＝｛x1，x2，…，xR｝；

4）取全矩陣W中對(duì)應(yīng)序列Λ2的行向量，組成的矩陣W2即為所求

（3）求取最佳近似分離權(quán)系數(shù)A：

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

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

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

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

2 模型試算

2.1 均勻模型

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

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

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

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

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

圖1 均勻黏聲各向異性介質(zhì)擬聲波方程（上）和純qP波方程（下）400ms時(shí)刻波場(chǎng)快照

圖2 純qP波方程在均勻黏聲各向異性介質(zhì)中400ms時(shí)刻波場(chǎng)快照

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

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

圖3 不同介質(zhì)類型在400ms時(shí)刻的波場(chǎng)快照對(duì)比

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

2.2 簡單TTI模型

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

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

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

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

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

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

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

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

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

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

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

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

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

表4 不同方法在不同時(shí)間步長時(shí)的CPU耗時(shí)

2.3 復(fù)雜Hess模型

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

圖8 三種方法在不同時(shí)間步長下得到的400ms時(shí)刻波場(chǎng)快照

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

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

2.4 復(fù)雜BP模型

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

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

圖10 二維Hess黏聲TTI介質(zhì)800ms時(shí)刻波場(chǎng)快照（a）及地震記錄（b）

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

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

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

圖13 二維BP2007黏聲TTI介質(zhì)4s時(shí)刻波場(chǎng)快照（a）及地震記錄（b）

3 結(jié)論與分析

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

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

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

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