預條件彈性介質最小二乘逆時偏移

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

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

0 引言

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

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

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

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

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

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

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

與背景介質和擾動介質相對應的是背景波場P=［ux，uz，σxx，σzz，σxz］和擾動波場δ u=［Δux，Δuz，Δσxx，Δσzz，Δσxz］?？偛▓鐾瑯訚M足波動方程，即

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

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

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

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

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

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

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

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

則伴隨方程可表示為

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

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

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

vP和vS的梯度可由鏈式公式求得，即

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

1.2 梯度預條件算子

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

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

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

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

對式（13）求導可得

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

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

1.3 預條件彈性波最小二乘逆時偏移

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

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

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

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

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

2 模型試算

2.1 鹽丘模型

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

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

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

圖2 ERTM 成像結果對比Fig.2 Comparison of ERTM imaging results

圖3 ELSRTM 成像結果對比Fig.3 Comparison of ELSRTM imaging results

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

分別從常規(guī)ELSRTM 和P-ELSRTM 的成像結果中抽取偏移距為4 900 m 的地震單道以對比各方法的保幅性（迭代30 次），結果如圖4 所示。從圖4可看出：①在模型的淺部區(qū)域由于低頻噪音和震源效應的存在，常規(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 成像結果單道振幅對比（第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 的成像結果與真實結果的歸一化殘差隨迭代次數(shù)變化的結果如圖5所示。從圖5 可以看出：P-ELSRTM 的成像反演結果在前60 次迭代中都能穩(wěn)定收斂，最終收斂到22.6%；常規(guī)ELSRTM 隨著迭代次數(shù)的增加只能收斂到52% 便停止收斂。因此，使用梯度預條件算子能夠使復雜模型成像結果的殘差收斂更快且達到更低的值。

在這次測試中，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結果的對比如表1 所列。由于計算漢森算子所花費的計算量相對于偏移運算較小，因此表1 中忽略了漢森算子的計算量，僅對比了偏移運算的計算量。P-ELSRTM 算法與ELSRTM 算法的計算量相近，但前者具有更高的成像質量、更快的收斂速度及更好的穩(wěn)定性。

表1 ELSRTM 成像結果對比Table 1 Comparison of ELSRTM imaging results

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

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

圖7 低信噪比數(shù)據(jù)成像結果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 結論與討論

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

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

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