預(yù)條件彈性介質(zhì)最小二乘逆時偏移

劉夢麗，徐興榮，王小衛(wèi)，胡書華，劉金濤

（中國石油勘探開發(fā)研究院西北分院，蘭州 730020）

0 引言

彈性波動理論的發(fā)展以及長期的勘探實踐均已表明，地震波是以矢量波的形式在地球介質(zhì)中傳播的。建立在聲波介質(zhì)假設(shè)上的問題只能主要圍繞縱波速度和密度進行，忽略了橫波，不符合精細勘探的要求［1-3］。為此針對精細巖性油氣藏勘探，發(fā)展彈性介質(zhì)假設(shè)下的成像方法迫在眉睫。逆時偏移是當(dāng)前偏移方法中較為精確的深度偏移方法且發(fā)展迅速，已經(jīng)從疊后發(fā)展到疊前，從二維到三維，從聲波到彈性波，從各向同性到各向異性［4-7］。彈性波逆時偏移主要分為2 個過程，首先構(gòu)建震源波場和檢波波場，利用成像條件得到成像結(jié)果。1986年Chang 等［8］首次完成了各向同性介質(zhì)彈性波疊前逆時偏移；隨后Chang 等［9-10］又深入研究，完成了2D和3D 彈性波逆時偏移；董良國等［11］分析不同地形起伏情況下自由邊界的情況，模擬出地表起伏情況下彈性波復(fù)雜的傳播現(xiàn)象，實現(xiàn)了彈性波起伏地表自由邊界條件的數(shù)值計算；為了深入研究彈性波波場傳播規(guī)律，陳可洋［12］將波印廷矢量應(yīng)用于多波多分量各向異性介質(zhì)彈性波波動方程中，實現(xiàn)了上行波、下行波、左行波和右行波的波場分離。針對彈性波成像的串?dāng)_問題，王維紅等［13］依據(jù)各向同性介質(zhì)一階速度-應(yīng)力方程組，利用Helmholtz 分解提取純縱波和純橫波波場，使用震源歸一化的互相關(guān)成像條件獲得純波成像，避免了縱橫波串?dāng)_問題。張偉等［14］利用高階交錯網(wǎng)格有限差分數(shù)值方法構(gòu)建了矢量的縱橫波波場，將震源歸一化的內(nèi)積成像條件應(yīng)用于分離后的純縱波和橫波矢量場，由此得到的轉(zhuǎn)換波成像避免了傳統(tǒng)彈性波成像方法中出現(xiàn)的極性反轉(zhuǎn)。成像條件是影響彈性逆時偏移的另一個重要因素。杜啟振等［15］基于激發(fā)時間成像條件實現(xiàn)了橫向各向同性介質(zhì)中的多波多分量疊前逆時偏移；張曉語等［16］基于能量守恒定理及能量密度，實現(xiàn)了彈性波能量互相關(guān)成像，能夠壓制背向散射；張智等［17］提出穩(wěn)定的激發(fā)振幅成像條件，在震源波場和檢波器波場的傳播過程中，保存最大能量密度的時刻和相應(yīng)的波場值，歸一化后獲得角度依賴的反射系數(shù)成像剖面，與歸一化成像條件相比，穩(wěn)定激發(fā)振幅成像條件具有更小的計算量。

常規(guī)的逆時偏移使用的是反偏移算子的共軛，而不是它的逆，因此造成成像結(jié)果模糊化，成像精度難以滿足精細油氣勘探的精度要求［18-19］。為了彌補常規(guī)逆時偏移的缺陷，從最小二乘思想的角度重新審視逆時偏移的過程［20-25］。近年來，對最小二乘逆時偏移（LSRTM）的研究逐漸引起了國內(nèi)外學(xué)者的關(guān)注。其中，為了提升效率，減少耗時，Dai 等［26-27］引入了相位編碼；針對高精度成像這一目標(biāo)，黃建平等［28］優(yōu)化了LSRTM 算法；LSRTM 還可應(yīng)用于黏聲介質(zhì)，李振春等［29］完成了黏聲介質(zhì)LSRTM 算法，改善成像質(zhì)量。考慮到LSRTM 的優(yōu)勢，國內(nèi)外學(xué)者將彈性波逆時偏移（ERTM）與最小二乘偏移（LSM）結(jié)合發(fā)展出彈性波最小二乘逆時偏移（EL‐SRTM）來進行復(fù)雜構(gòu)造的保幅成像處理［30-31］。

從彈性波波動方程出發(fā)，推導(dǎo)彈性波線性正演以及逆時偏移算子表達式；然后從誤差函數(shù)出發(fā)，推導(dǎo)Hessian 矩陣的表達式，以此作為梯度預(yù)條件算子，實現(xiàn)一種預(yù)條件彈性波最小二乘逆時偏移（PELSRTM）算法；最后，通過SEG/EAGE 二維鹽丘模型測試P-ELSRTM 對復(fù)雜模型成像時的有效性及其優(yōu)勢，并對比和分析常規(guī)ELSRTM 及P-ELSRTM對低信噪比炮數(shù)據(jù)的適應(yīng)性，以期證明預(yù)處理對反演過程中的改進效果。

1 彈性波最小二乘逆時偏移

1.1 線性化反偏移和偏移算子

在非均勻各向同性彈性介質(zhì)中，二維彈性波一階速度-應(yīng)力方程［30-31］可表示為

式中：Sxx和Szz分別為水平方向和垂直方向上的震源；ux和uz分別為介質(zhì)中在水平方向和垂直方向上質(zhì)點振動速度分量；σxx，σxz，σzz分別代表應(yīng)力分量；ρ為介質(zhì)質(zhì)量密度，kg/cm3，常數(shù)；λ和μ均為拉梅常數(shù)。拉梅常數(shù)λ和μ、密度ρ、P 波速度vP、S 波速度vS之間的關(guān)系為：

與背景介質(zhì)和擾動介質(zhì)相對應(yīng)的是背景波場P=［ux，uz，σxx，σzz，σxz］和擾動波場δ u=［Δux，Δuz，Δσxx，Δσzz，Δσxz］。總波場同樣滿足波動方程，即

將式（2）與式（1）相減，并且在模型擾動足夠小的假設(shè)下，用背景波場P代替總波場Ptotal，可得到基于Born 近似的關(guān)于擾動波場δ u的波動方程，可以表示為

式中：Δux和Δuz均為質(zhì)點振動速度擾動項；Δσxx，Δσxz，Δσzz均為應(yīng)力擾動項；P，S 波反射系數(shù)模型可用速度擾動δ vP，δ vS和背景速度vP，vS表示為mP(x)=2δ vP/vP，mS(x)=2δ vS/vS，反映了介質(zhì)速度擾動的相對大小。

從式（3）可以看出，要得到質(zhì)點振動速度擾動項和應(yīng)力擾動項需要經(jīng)歷2 個過程：第1 個過程是先通過式（1）得到背景波場P0；第2 個過程是利用式（3）進行計算，通過背景慢度模型和新的虛擬震源項再次激發(fā)，再傳播到地下x位置處形成散射波場。另外，從式（3）的右端表達式可以看出，質(zhì)點振動速度擾動項和應(yīng)力擾動項的產(chǎn)生是由背景介質(zhì)參數(shù)與介質(zhì)速度擾動作為新的震源，而且兩者與介質(zhì)速度擾動項以及應(yīng)力擾動項呈線性關(guān)系。

為了便于表示，用矩陣形式將式（3）簡寫成［32-35］

式中：m為反射系數(shù)模型的矩陣形式；d表示矩陣形式的觀測數(shù)據(jù)，d=（dx，dz，0，0，0）T；L為基于Born近似下的線性反偏移算子，它建立了介質(zhì)參數(shù)與正演數(shù)據(jù)之間的聯(lián)系。

采用時間2 階、空間10 階精度的交錯網(wǎng)格有限差分法對式（1）進行差分離散。在進行波場模擬時，由于人工截斷邊界會導(dǎo)致邊界反射，采用完全匹配層（PML）吸收邊界，來壓制邊界反射提高模擬區(qū)域的精度與信噪比。

本文借助伴隨狀態(tài)法得到伴隨算子及相應(yīng)的伴隨波場控制方程，由伴隨狀態(tài)法可得

式中：w為伴隨波場向量，L*為L的伴隨算子。

則伴隨方程可表示為

式中：δ d為殘差向量，δ d=[δ dxδ dz0 0 0]T。

為了得到與記錄數(shù)據(jù)最佳匹配的偏移結(jié)果，把最小二乘逆時偏移引入反演中，其誤差函數(shù)定義為

根據(jù)伴隨狀態(tài)方法，求解式（7）可得λ，μ的梯度表達式為

vP和vS的梯度可由鏈?zhǔn)焦角蟮茫?/p>

梯度計算是最小二乘成像的核心問題之一，對梯度做預(yù)條件，加速問題的收斂。本文中用于后續(xù)迭代更新的梯度為g=(δ vP，)δ vS。將式（9）代入式（7）可得δ vP和δ vS的梯度表達式，即

1.2 梯度預(yù)條件算子

式（7）所對應(yīng)的解的估計稱最小二乘偏移成像的結(jié)果，真實反射模型的最小二乘解：

傳統(tǒng)的偏移方法是通過反偏移算子的共軛轉(zhuǎn)置實現(xiàn)的，會造成地震信息缺失。從式（9）可看出，最小二乘偏移反演通過Hessian 算子L*L對偏移成像進行修正，難點就在于求解一個精確的Hessian并求逆。Hessian 算子H是誤差函數(shù)在m處對模型的二階導(dǎo)數(shù)，即

從式（12）可以看出，H由2 部分構(gòu)成，其中，為其非線性項且依賴于模型，由于非線性項的計算過于復(fù)雜，計算成本過高，通常使用Hessian 算子的線性項L*L來近似代替全漢森矩陣。通過Hessian 算子的逆矩陣作為梯度預(yù)處理算子，可以得到分辨率更高的成像結(jié)果，并且加快收斂。

式（1）用矩陣形式可表示為

對式（13）求導(dǎo)可得

通過漢森矩陣的線性項來近似代替全漢森矩陣，即

用拉梅常數(shù)表示的預(yù)條件算子表達式為

1.3 預(yù)條件彈性波最小二乘逆時偏移

利用共軛梯度法迭代求解式（7）得到模擬數(shù)據(jù)與觀測數(shù)據(jù)的差別最小的最佳反射系數(shù)模型，以P表示預(yù)條件算子，梯度g=()δ vP,δ vS，具體實現(xiàn)流程如下：

（1）輸入初始模型（初始模型賦0 值，即第1 次迭代的結(jié)果為逆時偏移的偏移剖面），觀測數(shù)據(jù)dobs，梯度誤差閾值δ。

（2）計算第i次迭代共軛梯度d ki和共軛方向修正因子β

（3）計算步長α并對成像結(jié)果進行更新，得到mi+1

（4）判斷是否滿足終止迭代條件，若滿足，則退出迭代并保存成像結(jié)果；不滿足，則回到第（2），（3）步繼續(xù)迭代直到滿足迭代終止條件，最后一次迭代即為最終的成像結(jié)果。

2 模型試算

2.1 鹽丘模型

采用SEG/EAGE 提供的鹽丘模型來測試常規(guī)ELSRTM 和P-ELSRTM 對模型的成像效果，縱橫波真實速度場和反射系數(shù)模型如圖1 所示。設(shè)定縱橫波速度比為1.73，密度為常數(shù)。模型橫向網(wǎng)格點數(shù)為645，縱向網(wǎng)格點數(shù)為150，采樣間隔為10 m。利用該模型主要考察P-ELSRTM 成像算法是否能夠使鹽下構(gòu)造準(zhǔn)確成像以及對振幅的補償效果：由于鹽丘的屏蔽作用，導(dǎo)致鹽下構(gòu)造的振幅較小，另外，鹽丘側(cè)翼發(fā)育的高陡小構(gòu)造會引起嚴(yán)重的構(gòu)造假象。觀測系統(tǒng)采用全接收方式，每炮的網(wǎng)格間隔數(shù)為10 m，共60 炮。震源類型為爆炸震源，震源函數(shù)采用主頻為25 Hz 的雷克子波。

圖1 SEG/EAGE 二維鹽丘模型Fig.1 SEG-EAGE 2-D salt dome model

ERTM 成像結(jié)果如圖2 所示，逆時偏移成像中使用全波動方程模擬波場，可以實現(xiàn)對陡傾角構(gòu)造的成像，但是在ERTM 剖面中含有較強的低波數(shù)噪聲和震源效應(yīng)［圖2（a）和圖2（b）］。應(yīng)用拉普拉斯算子去噪之后，相比未作濾波的ERTM 結(jié)果得到了改善［圖2（c）和圖2（d）］，模型淺部的低頻噪音得到了壓制，但成像能量不均衡。使用預(yù)條件算子和濾波之后，相比濾波后的ERTM 成像剖面質(zhì)量得到進一步提高［圖2（e）和圖2（f）］，鹽下構(gòu)造清晰，且模型淺部的低頻噪音和震源效應(yīng)也得到了更好地壓制。

圖2 ERTM 成像結(jié)果對比Fig.2 Comparison of ERTM imaging results

圖3 ELSRTM 成像結(jié)果對比Fig.3 Comparison of ELSRTM imaging results

為了討論預(yù)條件ELSRTM 的有效性及優(yōu)點，分別對常規(guī)ELSRTM 和預(yù)條件P-ELSRTM 成像結(jié)果進行對比（圖3）：①未作濾波的常規(guī)ELSRTM 在迭代次數(shù)為5 和30 次的成像結(jié)果中都殘留較嚴(yán)重的低波數(shù)噪音，影響對鹽丘之上沉積層的識別［圖3（a）—（d）］，而應(yīng)用高通濾波預(yù)條件算子后，成像剖面中的低頻噪音得到了壓制，剖面振幅更加均衡［圖3（e）—（h）］，證明了濾波算子的有效性；②濾波后的常規(guī)ELSRTM成像結(jié)果中鹽體邊界清晰，鹽下成像相對于ERTM 結(jié)果更為準(zhǔn)確，但仍存在深部成像能量不足的現(xiàn)象，并且對頂部低頻噪聲的壓制不徹底［圖3（e）—（h）］；③應(yīng)用照明補償預(yù)條件算子和高通濾波后的結(jié)果［圖3（i）—（l）］表明，采用預(yù)條件成像方法實現(xiàn)了對復(fù)雜模型深部小斷層較為清晰的成像，鹽下斷層及基底均能正確成像，并且成像剖面信噪比更高、能量更加均衡。

分別從常規(guī)ELSRTM 和P-ELSRTM 的成像結(jié)果中抽取偏移距為4 900 m 的地震單道以對比各方法的保幅性（迭代30 次），結(jié)果如圖4 所示。從圖4可看出：①在模型的淺部區(qū)域由于低頻噪音和震源效應(yīng)的存在，常規(guī)ELSRTM 的振幅曲線會有明顯跳動，而P-ELSRTM 的振幅曲線抖動幅度較小；②在模型的500 m，830 m，1 430 m，4 500 m（反射界面）處，常規(guī)ELSRTM 的振幅與真實反射系數(shù)之間有較大的差異，而P-ELSRTM 的振幅與真實反射系數(shù)更加接近。因此，在淺—中深層P-ELSRTM 都比常規(guī)ELSRTM 具有更好的保幅性。

圖4 ELSRTM 與P-ELSRTM 成像結(jié)果單道振幅對比（第1 次及第30 次迭代）Fig.4 Amplitude of single trace of ELSRTM and P-ELSRTM imaging results

圖5 P-ELSRTM 與ELSRTM 誤差隨迭代次數(shù)的變化圖（迭代60 次）Fig.5 Residual convergence curves of P-ELSRTM and ELSRTM

常規(guī)ELSRTM 和P-ELSRTM 的成像結(jié)果與真實結(jié)果的歸一化殘差隨迭代次數(shù)變化的結(jié)果如圖5所示。從圖5 可以看出：P-ELSRTM 的成像反演結(jié)果在前60 次迭代中都能穩(wěn)定收斂，最終收斂到22.6%；常規(guī)ELSRTM 隨著迭代次數(shù)的增加只能收斂到52% 便停止收斂。因此，使用梯度預(yù)條件算子能夠使復(fù)雜模型成像結(jié)果的殘差收斂更快且達到更低的值。

在這次測試中，CPU 型號為Intel（R）Xeon（R）CPU E5-2650 v2@2.60 GHz，P-ELSRTM和常規(guī)EL‐SRTM 均串行運行，迭代一次所花費的時間分別為4 667 s 和4 642 s。基于不同最優(yōu)化算法的LSRTM結(jié)果的對比如表1 所列。由于計算漢森算子所花費的計算量相對于偏移運算較小，因此表1 中忽略了漢森算子的計算量，僅對比了偏移運算的計算量。P-ELSRTM 算法與ELSRTM 算法的計算量相近，但前者具有更高的成像質(zhì)量、更快的收斂速度及更好的穩(wěn)定性。

表1 ELSRTM 成像結(jié)果對比Table 1 Comparison of ELSRTM imaging results

野外地震資料采集的數(shù)據(jù)中不可避免地含有環(huán)境以及人為因素造成的隨機噪聲。為了測試PELSRTM 對低信噪比數(shù)據(jù)的適應(yīng)性，對炮記錄添加低信噪比的隨機噪聲，使其X分量和Z分量信噪比為0.5 dB（圖6）。從ELSRTM 和P-ELSRTM 的成像結(jié)果（圖7）可看出，對于低信噪比數(shù)據(jù)，ELSRTM和P-ELSRTM 結(jié)果中都含有類似于隨機噪聲的偏移假象，但相比常規(guī)ELSRTM，P-ELSRTM 能夠更好地壓制低頻噪音。從殘差收斂曲線可以看出，P-ELSRTM 的收斂速度較常規(guī)ELSRTM 快，并且迭代60 次后其殘差較小（圖8）。因此，P-ELSRTM 在處理含噪數(shù)據(jù)方面有較大優(yōu)勢，抗噪性更好。

圖6 低信噪比單炮數(shù)據(jù)Fig.6 Low SNR shot data

圖7 低信噪比數(shù)據(jù)成像結(jié)果Fig.7 Imaging results of low SNR data

圖8 含噪音數(shù)據(jù)的P-ELSRTM 與ELSRTM 誤差隨迭代次數(shù)的變化圖（迭代60 次）Fig.8 Residual convergence curves of P-ELSRTM and ELSRTM with noise

3 結(jié)論與討論

（1）P-ELSRTM 和ELSRTM 能夠有效壓制淺層低頻噪音，消除震源效應(yīng)，有效改善常規(guī)ERTM的成像質(zhì)量，具有更高的成像分辨率。

（2）P-ELSRTM 在中深部具有更好的保幅性和分辨率，有利于深部及鹽下構(gòu)造的識別；P-ELSRTM比常規(guī)ELSRTM 的殘差曲線更快地收斂。在對含噪數(shù)據(jù)的成像測試中，P-ELSRTM 比常規(guī)ELSRTM可以收斂到一個更低的殘差水平。

（3）本文使用一階位移波動方程得到的縱橫波是耦合在一起的，因此通過迭代無法完全消除縱橫波串?dāng)_噪音。針對該問題，進行縱波和橫波波場矢量波場分離可進一步提高成像精度。