段心標(biāo) 王華忠 鄧光校

(①同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院波現(xiàn)象與智能反演成像研究組，上海 200092；②中國(guó)石化石油物探技術(shù)研究院，江蘇南京 211103； ③中國(guó)石化西北油田分公司，新疆烏魯木齊 830011)

0 引言

從地震數(shù)據(jù)中準(zhǔn)確提取并利用復(fù)雜的地下巖性與儲(chǔ)層信息，是油氣勘探的關(guān)鍵環(huán)節(jié)。地震波成像是連接地震觀測(cè)數(shù)據(jù)與儲(chǔ)層解釋的橋梁，常規(guī)偏移成像技術(shù)[1-5]主要目的是獲得能精細(xì)描述地質(zhì)構(gòu)造的地震成果資料，用于部署和落實(shí)井位。

隨著油氣勘探程度的日趨深入和地震成像技術(shù)的不斷進(jìn)步，直接估算地下彈性參數(shù)的反演成像方法受到越來(lái)越多的關(guān)注。早在20世紀(jì)60年代，Backus等[6-7]從數(shù)學(xué)角度研究了地球物理反問(wèn)題，奠定了早期地球物理反演的理論基礎(chǔ)。Tarantola[8-9]基于Bayes估計(jì)理論，建立了一套地震波全波形反演的理論方法體系，推動(dòng)了地震反演成像技術(shù)的快速發(fā)展。然而，全波形反演具有很強(qiáng)的非線性特征，難以由觀測(cè)數(shù)據(jù)直接反演得到較準(zhǔn)確的介質(zhì)彈性參數(shù)，無(wú)法滿足油藏描述的精度需求。實(shí)際中應(yīng)用的依然是經(jīng)典的地震波成像處理和儲(chǔ)層描述流程，而如何獲得更高偏移成像質(zhì)量是其中重要環(huán)節(jié)。

最小二乘偏移[10-18]是常規(guī)偏移的后續(xù)發(fā)展，其核心思想是在線性反演理論下求解模型空間的精確解，相當(dāng)于在常規(guī)偏移構(gòu)造成像基礎(chǔ)上，考慮了振幅補(bǔ)償和波形校正，從而得到更保幅和更高分辨率的成像結(jié)果。

經(jīng)典最小二乘偏移以散射波線性化表達(dá)為基礎(chǔ)，認(rèn)為模型參數(shù)是速度(或慢度)的高波數(shù)擾動(dòng)，其成像實(shí)質(zhì)是估計(jì)介質(zhì)的散射強(qiáng)度。Zhang等[19-20]提出基于逆時(shí)偏移疊加剖面的反偏移方法，并形成相應(yīng)的最小二乘偏移技術(shù)，但并未給出反偏移方法的詳細(xì)理論依據(jù)。陳生昌等[21-22]提出地震反射波數(shù)據(jù)偏移成像方法和相應(yīng)的反射波動(dòng)方程最小二乘偏移方法，并引入“反射率”概念，即最小二乘偏移的實(shí)質(zhì)是估計(jì)反射面上的反射率。本文認(rèn)為：因地下介質(zhì)以層狀為主，獲取高保真反射系數(shù)是地震成像走向儲(chǔ)層刻畫(huà)的關(guān)鍵，估算不同角度的反射系數(shù)是最小二乘偏移的目標(biāo)。然而，考慮到反演不同角度反射系數(shù)的計(jì)算量巨大，且高質(zhì)量偏移疊加成像依然是當(dāng)今實(shí)際應(yīng)用的主要目標(biāo)，因此以最小二乘偏移獲得忽略角度信息的平均反射系數(shù)具有現(xiàn)實(shí)意義。

本文首先簡(jiǎn)介線性最小二乘偏移方法基礎(chǔ)理論，重點(diǎn)討論散射理論線性化表達(dá)方法和反射理論線性化表達(dá)方法，并分別建立散射強(qiáng)度、反射系數(shù)與地震數(shù)據(jù)的線性關(guān)系式； 然后從偏移成像實(shí)際應(yīng)用需求出發(fā)，選擇基于反射理論并忽略角度信息的Kirchhoff正演反偏移表達(dá)方法，實(shí)現(xiàn)與實(shí)際應(yīng)用需求相一致的最小二乘逆時(shí)偏移技術(shù)； 最后，通過(guò)模型數(shù)據(jù)和實(shí)際資料驗(yàn)證本文方法的效果。

1 方法原理

1.1 線性最小二乘偏移方法

基于Bayes估計(jì)理論，地震反演成像可表示為后驗(yàn)概率密度最大化問(wèn)題； 在地震反演通常假設(shè)條件下，進(jìn)而可轉(zhuǎn)化為如下目標(biāo)函數(shù)最小化問(wèn)題

(1)

式中：m為地下介質(zhì)參數(shù)，如速度、散射強(qiáng)度、反射系數(shù)或波阻抗等；dobs為觀測(cè)數(shù)據(jù)；L(·)表示波場(chǎng)傳播算子。該式非線性反演問(wèn)題求解難度大，需對(duì)非線性反演問(wèn)題逐步進(jìn)行線性化，依次反演參數(shù)的低、中、高波數(shù)信息。最小二乘偏移是一種線性化的地震反演成像方法，用于估計(jì)介質(zhì)參數(shù)的高波數(shù)成分，其目標(biāo)函數(shù)可表示為

(2)

式中L為矩陣形式的線性化波場(chǎng)傳播算子。

令目標(biāo)函數(shù)梯度為零，則有如下的法方程

LTLm=LTdobs

(3)

最小二乘偏移成像結(jié)果表示為

(4)

式中：LTL為Hessian矩陣，可用H表示；LTdobs是對(duì)觀測(cè)數(shù)據(jù)進(jìn)行常規(guī)偏移成像。最小二乘偏移是Hessian矩陣的逆H-1作用于常規(guī)偏移成像，能夠消除Hessian矩陣導(dǎo)致的成像模糊化效應(yīng)，提高成像分辨率和保幅性，減少成像假象。

由于“目標(biāo)函數(shù)梯度為零”僅是一種理想假設(shè)，且Hessian矩陣非常龐大，其計(jì)算和存儲(chǔ)都很困難，Hessian矩陣的求逆更難以完成，因而無(wú)法直接求解法方程以實(shí)現(xiàn)最小二乘偏移。為此，通常采用迭代法求解最小二乘偏移問(wèn)題，假設(shè)目標(biāo)函數(shù)在初始值附近滿足二次型，從初始解逐步逼近真實(shí)解。

1.2 地震波場(chǎng)的線性化表達(dá)

1.2.1 散射理論的Born正演線性化表達(dá)

一階Born近似假設(shè)下，散射波場(chǎng)可寫(xiě)成

(5)

(6)

進(jìn)而，散射波場(chǎng)可表示為

G0(x′,xs,ω)dx′

(7)

式(7)定義了一個(gè)Fredholm第一類積分方程，它描述了一個(gè)線性問(wèn)題。其積分核可寫(xiě)成

(8)

在散射假設(shè)下，可認(rèn)為地表接收到的攜帶地下信息的信號(hào)就是散射波場(chǎng)記錄，則可將式(7)寫(xiě)成如下線性算子表達(dá)式

Lm=dobs

(9)

此式中的m設(shè)定為慢度平方擾動(dòng)，則該式表達(dá)了散射強(qiáng)度與地震數(shù)據(jù)之間的線性關(guān)系。

1.2.2 反射理論的Kirchoff正演線性化表達(dá)

反射波理論假設(shè)下，地下界面Ω上入射波場(chǎng)引起的二次場(chǎng)可表示為

(10)

地下界面Ω上波場(chǎng)的法向?qū)?shù)為

(11)

利用WKBJ近似方法[23]，有

(12)

由Green定理，反射面上方地面上任意一點(diǎn)x的反射波場(chǎng)可寫(xiě)成如下積分式

(13)

式中Gr(x,x′,ω)是反射波傳播的Green函數(shù)(圖1)。

圖1 波場(chǎng)傳播示意圖

將式(10)～式(12)代入式(13)，可得

(14)

根據(jù)程函方程關(guān)系，式(14)可進(jìn)一步寫(xiě)成

Gs(x′,xs,ω)Gr(x,x′,ω)dΩ

(15)

式中s(x′)為慢度矢量。

式(15)建立了地震數(shù)據(jù)與不同方向反射系數(shù)之間的聯(lián)系[24-25]。然而，在緩變速度結(jié)構(gòu)中，不同角度的反射波能量較均衡； 在速度劇變結(jié)構(gòu)中，不同角度的反射系數(shù)難以估計(jì)。同時(shí)，考慮到直接基于式(15)進(jìn)行波場(chǎng)模擬的計(jì)算量巨大，因此建立地震數(shù)據(jù)與平均反射系數(shù)之間的關(guān)系就具有現(xiàn)實(shí)意義。若不考慮背景速度對(duì)振幅的影響，式(15)就退化為

(16)

式中：D表示散射體；R′表示平均反射系數(shù)。

在反射假設(shè)下，認(rèn)為地表接收到的攜帶地下信息的信號(hào)就是反射波場(chǎng)記錄。式(16)同樣可寫(xiě)成線性算子(式(9))表達(dá)形式，此時(shí)m表示忽略地下角度信息的平均反射系數(shù)，其積分核為

(17)

1.3 最小二乘逆時(shí)偏移實(shí)現(xiàn)方法

逆時(shí)偏移技術(shù)現(xiàn)已在實(shí)際處理中得到廣泛應(yīng)用，其成像基礎(chǔ)是Claerbout[26]給出的成像方法：波場(chǎng)外推+成像條件。成像條件為：成像點(diǎn)處，入射波的到達(dá)時(shí)等于反射波的出發(fā)時(shí)。顯然，常規(guī)逆時(shí)偏移是估計(jì)地下介質(zhì)的反射系數(shù)。由于地下介質(zhì)以層狀為主，以定位反射界面的位置或產(chǎn)生角度反射系數(shù)道集為目標(biāo)的反射系數(shù)成像是當(dāng)前地震成像的主流需求。

如前所述，由于角度反射系數(shù)反演的計(jì)算量巨大，因而實(shí)際應(yīng)用中，最小二乘偏移的目的是在常規(guī)逆時(shí)偏移基礎(chǔ)上進(jìn)一步提高偏移疊加成像質(zhì)量，使成像結(jié)果分辨率更高，成像假象更少。

在反演成像方法中，反偏移采用Kirchhoff正演思路，具體迭代反演步驟如下。

(1)以設(shè)定的偏移速度場(chǎng)對(duì)觀測(cè)數(shù)據(jù)dobs作逆時(shí)偏移成像，作為最小二乘偏移的初始值m0。

(2)對(duì)第i次迭代后偏移成像值mi(i=0,1,2,…)做反偏移，正演模擬預(yù)測(cè)數(shù)據(jù)dcal，反偏移公式為

(18)

(3)計(jì)算觀測(cè)數(shù)據(jù)dobs與預(yù)測(cè)數(shù)據(jù)dcal的殘差波場(chǎng)δd，當(dāng)殘差滿足設(shè)定條件時(shí)迭代停止。

(4)對(duì)δd進(jìn)行逆時(shí)偏移成像，得到目標(biāo)泛函的梯度gi。

(5)計(jì)算共軛梯度方向

(19)

(6)計(jì)算迭代步長(zhǎng)

(20)

(7)更新偏移成像值

mk+1=mk-μk+1hk

(21)

(8)轉(zhuǎn)到步驟(2)，展開(kāi)新一輪迭代。

由于迭代算法計(jì)算量巨大，三維最小二乘逆時(shí)偏移中采用了GPU加速計(jì)算技術(shù)，在此不再贅述。

2 數(shù)值模型試驗(yàn)及實(shí)際應(yīng)用

首先利用二維Sigsbee2a模型(圖2a)做數(shù)值試驗(yàn)，驗(yàn)證本文方法效果； 進(jìn)而將三維最小二乘逆時(shí)偏移技術(shù)應(yīng)用于實(shí)際地震資料成像處理。

圖2 Sigsbee2a模型及其成像結(jié)果

2.1 Sigsbee2a模型數(shù)據(jù)

為了便于計(jì)算，把Sigsbee2a經(jīng)典模型速度值單位轉(zhuǎn)為m/s，并重新定義速度網(wǎng)格，新網(wǎng)格的尺寸為3200×1200，間隔為5m×5m。對(duì)該模型正演模擬120炮數(shù)據(jù)，炮間距為100m，檢波點(diǎn)間隔為20m，炮檢距范圍是-2000～2000m，記錄長(zhǎng)度為7s。

利用本文所提最小二乘逆時(shí)偏移方法進(jìn)行反演迭代，得到第20次迭代結(jié)果(圖2c)，并與常規(guī)逆時(shí)偏移剖面(圖2b)進(jìn)行對(duì)比； 圖3是目標(biāo)泛函迭代誤差分析曲線，可見(jiàn)反演逐步收斂。從整體剖面及鹽丘左下局部放大圖(圖4)的對(duì)比可見(jiàn)，相對(duì)于常規(guī)逆時(shí)偏移方法，本文方法的成像效果(圖4b)顯著改善，分辨率變高，深部成像更清晰。

圖3 迭代誤差曲線

圖4 Sigsbee2a模型的常規(guī)逆時(shí)偏移(a)和最小二乘逆時(shí)偏移(b)成像的局部放大

2.2 實(shí)際資料

本次試驗(yàn)所用實(shí)際數(shù)據(jù)取自中國(guó)西北M油田，測(cè)試處理的目的是提高儲(chǔ)層的成像精度，力求對(duì)奧陶系縫洞體及小斷裂精確成像。所選數(shù)據(jù)滿覆蓋面積為401km2，共48000炮，成像面元是30m×30m。利用100GPU節(jié)點(diǎn)(每節(jié)點(diǎn)4塊K10卡)3次迭代，共用時(shí)59天。

圖5和圖6是兩條不同測(cè)線的常規(guī)逆時(shí)偏移與最小二乘逆時(shí)偏移成像剖面對(duì)比，可見(jiàn)最小二乘逆時(shí)偏移剖面中“串珠”成像收斂更好，對(duì)深層大斷裂、層間小斷層的刻畫(huà)更清晰，T78(奧陶系蓬萊壩組頂界面)以下成像能量更強(qiáng)，成像分辨率得到提高。對(duì)比奧陶系目的層的成像頻譜(圖7)，可見(jiàn)最小二乘逆時(shí)偏移成像拓寬了頻譜(圖7b，頻帶展寬5Hz)； T74(一間房組頂界面)下40～120ms均方根振幅，也證實(shí)本文方法成像結(jié)果(圖8b)對(duì)“串珠”刻畫(huà)更清晰。

以上結(jié)果均展示了最小二乘逆時(shí)偏移成像在提高成像分辨率、振幅均衡性及成像精度等方面的效果。

圖5 L1測(cè)線的常規(guī)逆時(shí)偏移(a)和最小二乘逆時(shí)偏移(b)成像剖面

圖6 L2測(cè)線的常規(guī)逆時(shí)偏移(a)和最小二乘逆時(shí)偏移(b)成像剖面

圖7 常規(guī)逆時(shí)偏移(藍(lán)線)和最小二乘逆

3 結(jié)論與討論

最小二乘偏移是一種線性化的地震反演成像方法，線性近似可以是“背景速度+速度擾動(dòng)”模式，也可以是“背景速度+反射系數(shù)”模式。由于地下介質(zhì)以層狀為主，現(xiàn)今業(yè)界常規(guī)偏移成像的主要目標(biāo)是反射系數(shù)。本文分析了兩種線性化模式的差異，選擇基于忽略角度信息的反射系數(shù)的反偏移正問(wèn)題表達(dá)方法，提出一種估計(jì)平均反射系數(shù)的最小二乘逆時(shí)偏移方法。理論模型測(cè)試結(jié)果表明本方法能提高成像分辨率，深部成像更清晰。將研發(fā)的三維最小二乘逆時(shí)偏移技術(shù)應(yīng)用于中國(guó)西部探區(qū)資料處理，所得“串珠”成像比常規(guī)逆時(shí)偏移結(jié)果收斂更好，深層大斷裂、層間小斷層的刻畫(huà)更清晰，成像分辨率獲得提高，充分體現(xiàn)了最小二乘逆時(shí)偏移成像在提高成像分辨率和振幅均衡性及提高成像精度方面的優(yōu)勢(shì)。

圖8 常規(guī)逆時(shí)偏移(a)與最小二乘逆時(shí)偏移(b)的均方根振幅平面屬性對(duì)比

然而，目前最小二乘逆時(shí)偏移實(shí)際應(yīng)用效果仍不甚理想，主要原因包括反偏移得到的無(wú)噪聲一次波模擬數(shù)據(jù)與實(shí)際觀測(cè)記錄差異較大、實(shí)際地震子波未知、現(xiàn)有建模精度達(dá)不到最小二乘逆時(shí)偏移要求及三維技術(shù)計(jì)算量大等。為了降低最小二乘逆時(shí)偏移對(duì)速度模型和地震子波的依賴，近年也發(fā)展了動(dòng)態(tài)拉伸時(shí)差校正方法和空變子波估計(jì)方法，一定程度上減弱了反演的病態(tài)性。但考慮到反問(wèn)題的復(fù)雜性，這些應(yīng)用瓶頸問(wèn)題并未得到有效解決，它們會(huì)持續(xù)伴隨最小二乘逆時(shí)偏移的發(fā)展過(guò)程，有待逐步克服和解決。