一種基于改進(jìn)快速掃描法的多尺度近地表層析方法

蔡杰雄,王靜波

(1.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;2.中國石油化工股份有限公司勘探分公司,四川成都610041)

初至波走時(shí)易識(shí)別且穩(wěn)定,在疊前偏移[1]、微震定位[2]等領(lǐng)域得到了廣泛應(yīng)用,尤其在復(fù)雜地區(qū)的初至波層析近地表建模[3-5]中得到了很好的應(yīng)用。但在復(fù)雜、大數(shù)據(jù)地震勘探條件下,目前的初至波層析速度建模精度難以支撐靜校正、真地表偏移等技術(shù)的實(shí)際應(yīng)用。因此,提高初至波層析速度建模精度和效率迫在眉睫。

對(duì)于復(fù)雜工區(qū)的初至波層析,傳統(tǒng)的彎曲法[6]和打靶法[7]計(jì)算效率低,且存在射線陰影區(qū)。VIDALE[8-9]基于波前面求解程函方程計(jì)算初至波走時(shí),采用全局計(jì)算策略保證了計(jì)算的正確性,提高了計(jì)算效率,利用局部計(jì)算方法提高了計(jì)算精度和穩(wěn)定性。在全局計(jì)算策略中,VIDALE[9]提出了盒式擴(kuò)展法,計(jì)算效率高,但該技術(shù)無法適應(yīng)地下介質(zhì)速度劇烈變化情況;SETHIAN[10]提出了快速推進(jìn)法(Fast Marching Method,FMM),通過模擬波前面并不斷排序?qū)ふ易邥r(shí)最小的點(diǎn)向外計(jì)算,能夠適應(yīng)強(qiáng)變速介質(zhì),但極大地降低了計(jì)算效率;ZHAO[11]提出了快速掃描法(Fast Sweeping Method,FSM),在不同方向進(jìn)行掃描,從而獲得初至波走時(shí)在全空間的分布,該方法可以在任意維度進(jìn)行計(jì)算,在局部計(jì)算方法相同的情況下,有效提升了計(jì)算效率[12]。在局部計(jì)算方法中,VIDALE[8]采用一階差分格式;TRIER等[13]引入逆風(fēng)有限差分法求解程函方程,以提高計(jì)算的穩(wěn)定性;ZHANG等[14]和XIONG等[15]分別將差分格式擴(kuò)展到三階和五階,以提高計(jì)算精度;ASAKAWA等[16]提出了線性插值法,采用線性插值和惠更斯原理計(jì)算局部走時(shí);張賽民等[17]采用拋物型插值提高計(jì)算精度。但上述方法都沒有結(jié)合波前擴(kuò)展方法,引入了大量無用計(jì)算,影響計(jì)算效率。孫章慶等[18]利用“窄帶”技術(shù)減少了無用的走時(shí)和射線信息計(jì)算,但是該技術(shù)需要排序,限制了方法的計(jì)算效率。

在反演方面,線性反演方法的主要思想是采用線性近似建立層析方程組,通過不斷修改和迭代獲得最終的速度模型。但是層析方程組數(shù)據(jù)量龐大且嚴(yán)重病態(tài),難以保證計(jì)算的正確性和穩(wěn)定性。另一種方法是基于局部最優(yōu)化思想,通過對(duì)目標(biāo)函數(shù)的局部線性化近似獲得速度的修改方向與步長,然后通過迭代實(shí)現(xiàn)非線性反演問題的線性化求解。獲得梯度的主要方法有伴隨狀態(tài)法[4,19]和散射積分法[20]。前者無需進(jìn)行射線追蹤,但其預(yù)條件在射線走時(shí)層析成像中不易實(shí)現(xiàn);后者通過顯式地計(jì)算核函數(shù),并與走時(shí)殘差向量相乘實(shí)現(xiàn)梯度的計(jì)算,但與傳統(tǒng)的層析成像方法相同,散射積分法面臨內(nèi)存占用大的問題。李勇德等[21]充分利用矩陣元素的物理意義和Hessian矩陣的特性,在射線層析中引入了改進(jìn)散射積分法[22],有效降低了層析過程中的內(nèi)存占用;劉玉柱等[23]利用不同偏移距作為權(quán)函數(shù)提高淺層反演精度;劉小民等[24]將基于胖射線理論的初至波走時(shí)層析流程應(yīng)用于疊前深度偏移。

為提高初至波層析的計(jì)算精度和效率,在正演方面,在傳統(tǒng)快速掃描算法的基礎(chǔ)上,提出基于改進(jìn)快速掃描法的初至波走時(shí)計(jì)算方法;在反演方面,在改進(jìn)散射積分算法的基礎(chǔ)上,實(shí)現(xiàn)了多尺度反演策略。最后將提出的近地表建模方法應(yīng)用于丁山工區(qū)實(shí)際資料處理,驗(yàn)證了方法的有效性。

1 改進(jìn)的走時(shí)計(jì)算方法

1.1 基于雙線性插值的局部計(jì)算方法

如圖1所示,當(dāng)A,B,C,D點(diǎn)的走時(shí)tA,tB,tC,tD已知時(shí),需要計(jì)算點(diǎn)F的走時(shí)。為了簡單起見,僅討論網(wǎng)格間距(h)相等的情況。假設(shè)走時(shí)符合雙線性假設(shè)(如圖1所示,平面波在x,y兩個(gè)方向上均為線性變化),則E1,E2,E點(diǎn)的走時(shí)tE1,tE2,tE可以表示為:

(1)

此時(shí)點(diǎn)F的走時(shí)可以表示為:

(2)

其中,s表示慢度。將公式(1)代入公式(2)可得:

(3)

根據(jù)費(fèi)馬原理,點(diǎn)F的走時(shí)是公式(3)的最小值,對(duì)公式(3)分別求x,y的偏導(dǎo)數(shù)并令其偏導(dǎo)數(shù)等于0可得:

(4)

方程(4)是一個(gè)走時(shí)計(jì)算的超越方程,沒有解析解。張東等[25]采用網(wǎng)格界面剖分法求得固定節(jié)點(diǎn)上的近似解;劉鋒等[26]通過最速下降法求得數(shù)值解;孫章慶等[18]利用平面波假設(shè)構(gòu)造更簡單的公式求得解析解。本文借鑒孫章慶等提出的方法進(jìn)行求解。

圖1 雙線性插值法局部圖解[18]

平面波先經(jīng)過C點(diǎn),然后經(jīng)過D,B點(diǎn)到達(dá)A點(diǎn),由平面波假設(shè)可得這4個(gè)點(diǎn)的走時(shí)方程:

tA-tD=tC-tB

(5)

則公式(4)化簡為:

(6)

分析其單調(diào)性并整理得:

(7)

但是這種方法僅適用于平面波的情況,為了使其精度更高,在不符合平面波假設(shè)的區(qū)域可以利用一階迎風(fēng)差分格式進(jìn)行計(jì)算。

為了考慮射線所有可能的方向,需要計(jì)算點(diǎn)F相鄰的所有面元,但計(jì)算效率很低,可以采用迎風(fēng)策略,選擇射線最可能的方向進(jìn)行計(jì)算。具體方法如下:

1) 假設(shè)點(diǎn)F的空間位置為(i,j,k);

2) 選擇在z方向上靠近點(diǎn)F的兩點(diǎn)中走時(shí)較小的為點(diǎn)A(如圖1所示),其走時(shí)tA=min(ti,j,k+1,ti,j,k-1);

3) 選擇在x方向上靠近點(diǎn)A的兩點(diǎn)中走時(shí)較小的為點(diǎn)B(如圖1所示),其走時(shí)tB=min(ti+1,j,Ak,ti-1,j,Ak),其中,下標(biāo)Ak表示點(diǎn)A的z坐標(biāo);

4) 選擇在y方向上靠近點(diǎn)A的兩點(diǎn)中走時(shí)較小的為點(diǎn)D(圖1),其走時(shí)tD=min(ti,j+1,Ak,ti,j-1,Ak);

5) 由點(diǎn)D和點(diǎn)B共同確定點(diǎn)C,其走時(shí)tC=tBi,Dj,Ak,其中,Bi表示點(diǎn)B的x坐標(biāo),Dj表示點(diǎn)D的y坐標(biāo);

6) 按照步驟2)到步驟5)的方法可以確定點(diǎn)F的走時(shí)值,比較該走時(shí)值與計(jì)算前的走時(shí)值,選取較小的值作為該點(diǎn)計(jì)算結(jié)果。

1.2 數(shù)值模型計(jì)算

為比較改進(jìn)快速掃描法與傳統(tǒng)一階差分形式快速掃描法的優(yōu)劣,本文分別將兩種方法應(yīng)用于簡單均勻模型、復(fù)雜起伏地表模型。

1.2.1 簡單均勻模型

模型大小為1000m×1000m×1000m,炮點(diǎn)位于(500m,500m,500m),網(wǎng)格間距為10m×10m×10m,采用1000m/s均勻速度。圖2給出了在不同截面位置本文方法與傳統(tǒng)方法的誤差分布。由圖2可見,本文方法零誤差的范圍和角度都好于傳統(tǒng)方法。傳統(tǒng)方法僅在水平和豎直方向誤差較小,但本文方法在水平方向、45°方向和135°方向誤差都較小。在計(jì)算效率不受影響的前提下,本文方法的最大相對(duì)誤差降低為傳統(tǒng)方法的48%,同時(shí)本文方法計(jì)算時(shí)間為0.82s,傳統(tǒng)方法為0.83s,計(jì)算效率幾乎沒有改變。

1.2.2 復(fù)雜起伏地表模型

圖3給出了2.5維起伏地表Marmousi模型,模型大小為7450m×1000m×3000m,網(wǎng)格間距為25m×10m×10m,炮點(diǎn)位于(2500m,500m,0)。共5條接收線,線間距200m。每條測線檢波器x坐標(biāo)最小為1250m,最大為6250m,道間距25m。圖4a為采用本文方法得到的y=500m處的走時(shí)等值線,在起伏地表情況下,本文方法計(jì)算結(jié)果沒有出現(xiàn)奇異值,且正常產(chǎn)生了回轉(zhuǎn)波,說明本文方法適應(yīng)起伏地表復(fù)雜速度模型。

為了提高計(jì)算效率,僅計(jì)算炮檢距范圍內(nèi)地下2000m深度內(nèi)的走時(shí),結(jié)果如圖4b所示。本文方法計(jì)算全部走時(shí)場需4.21s,計(jì)算2000m深度范圍走時(shí)場需3.36s,傳統(tǒng)方法計(jì)算全部走時(shí)場需6.11s,本文方法計(jì)算時(shí)間為傳統(tǒng)方法的55%。需要注意的是,計(jì)算范圍的選擇應(yīng)該根據(jù)實(shí)際情況確定,否則會(huì)影響回轉(zhuǎn)波和折射波的產(chǎn)生。

圖2 不同截面位置本文方法(左)與傳統(tǒng)方法(右)的相對(duì)誤差分布a z=0; b x=500m; c x=y

圖3 2.5維起伏地表Marmousi模型

圖4 采用本文方法得到的y=500m處的全局(a)與部分(b)速度模型與走時(shí)等值線

2 多尺度改進(jìn)散射積分法走時(shí)層析

2.1 改進(jìn)散射積分法預(yù)條件最速下降方向計(jì)算

初至波走時(shí)層析通過計(jì)算走時(shí)與拾取走時(shí)的差計(jì)算出模型修正量,不斷迭代逼近真實(shí)模型。假設(shè)某個(gè)模型空間大小為m,數(shù)據(jù)空間大小為n,根據(jù)走時(shí)與慢度的關(guān)系可以建立走時(shí)矩陣:

Ks=t

(8)

其中,K為走時(shí)層析中的敏感核函數(shù)(Fréchet導(dǎo)數(shù)矩陣),其規(guī)模為n×m,s為空間慢度向量,t為初至波走時(shí)向量。對(duì)公式(8)進(jìn)行線性近似可以獲得慢度的修正量與走時(shí)差之間的關(guān)系:

KΔs=Δt

(9)

建立預(yù)條件最速下降方向迭代公式:

sl+1=sl+αlpl

(10)

式中:αl表示第l輪迭代的步長;pl表示第l輪迭代的梯度;sl表示第l輪空間內(nèi)各點(diǎn)的慢度。

(11)

p=(KTK)-1KTΔt

(12)

其中,Ha表示Hessian矩陣。改進(jìn)的散射積分法[21]將(12)式中的KTΔt和KTK表示為多個(gè)向量-標(biāo)量乘積的累加運(yùn)算,由于(KTK)-1主對(duì)角線占優(yōu),這里用其對(duì)角線H0代替。為了避免由于對(duì)角元素中存在很小的數(shù)而導(dǎo)致的不穩(wěn)定性,需要引入正則化項(xiàng),所以最終的梯度為:

p=(H0+λI)-1KTΔt

(13)

其中,I表示單位矩陣,λ表示一個(gè)較小的數(shù)。將(13)式代入(10)式,即可實(shí)現(xiàn)預(yù)條件最速下降法初至波走時(shí)層析。

2.2 多尺度層析策略

對(duì)于地下介質(zhì)速度隨深度線性增加的初始模型,地下深度z處的速度可以表示為v=v0+βz,射線可以達(dá)到的最大深度h可以表示為:

(14)

其中,v0表示地表速度,β表示速度梯度,h表示深度,x表示最大炮檢距?？梢园l(fā)現(xiàn),當(dāng)使用同一梯度模型作為初始模型進(jìn)行初至層析建模時(shí),隨著偏移距的增大,射線可以達(dá)到的深度也隨之增大;對(duì)于同樣偏移距,隨著速度梯度的增大,射線達(dá)到的深度也隨之增大。

鑒于射線必須經(jīng)過近地表淺層,在更新過程中會(huì)模糊淺層細(xì)節(jié),所以在深層更新結(jié)束后逐步減小最大炮檢距,刻畫淺層速度模型,實(shí)現(xiàn)反演深度的多尺度策略。為了建立更高精度的速度模型,采用先深后淺的反演策略,先用全炮檢距數(shù)據(jù)進(jìn)行層析,更新近地表深層信息。由于高頻假設(shè),射線層析集中在一條沒有寬度的“線”上,層析矩陣大量零元素會(huì)導(dǎo)致反演不穩(wěn)定,為此,采用較大的網(wǎng)格減少反演的不穩(wěn)定性,并隨著迭代次數(shù)的增加,逐步減小網(wǎng)格大小來提高精度,實(shí)現(xiàn)反演精度的多尺度策略。這種反演深度與反演網(wǎng)格同時(shí)變化的多尺度建模策略能有效提高速度建模精度與穩(wěn)定性。

2.3 理論數(shù)據(jù)測試

采用SEG起伏地表SEAM模型數(shù)據(jù)對(duì)本文方法進(jìn)行了測試,圖5為y坐標(biāo)6500m處的真實(shí)速度模型。該數(shù)據(jù)x,y,z方向網(wǎng)格采樣點(diǎn)分別為1447,1250,800個(gè),網(wǎng)格大小均為10m。起伏地表上部填充的速度為340m/s。觀測數(shù)據(jù)是高精度彈性波模擬數(shù)據(jù)的初至拾取結(jié)果。

圖6為y坐標(biāo)6500m處全偏移距層析反演的速度模型?？梢钥闯?通過改進(jìn)散射積分法初至波走時(shí)層析可以獲得平滑的背景速度。為了提高反演分辨率,采用多尺度反演策略,即隨著迭代次數(shù)的增加,輸入的最大偏移距分別為4000,2000,1000和500m,網(wǎng)格大小分別為100,50,20,10m。每一個(gè)反演參數(shù)內(nèi)部迭代20次。抽取y坐標(biāo)為6500m處的剖面進(jìn)行分析。圖7為反演的速度模型,可以看出,多尺度反演策略增加了近地表的高頻成分。圖8對(duì)比了不同位置縱向抽取的速度函數(shù),可以看出,采用分偏移距反演策略之后,淺部的反演精度得到了有效提高。

圖5 y坐標(biāo)為6500m處的真實(shí)速度模型

圖6 y坐標(biāo)為6500m處全偏移距反演的速度模型

圖7 y坐標(biāo)為6500m處分偏移距反演的速度模型

圖8 不同位置縱向抽取的速度函數(shù)a x坐標(biāo)為5000m,y坐標(biāo)為6500m; b x坐標(biāo)為7000m,y坐標(biāo)為6500m; c x坐標(biāo)為9000m,y坐標(biāo)為6500m

3 實(shí)際應(yīng)用

丁山工區(qū)為復(fù)雜山地頁巖氣勘探區(qū),地表起伏劇烈,高速巖體出露地表,導(dǎo)致地表橫向速度變化快。前期處理采用常規(guī)初至波層析方法,得到的最終單炮初至走時(shí)擬合平均誤差達(dá)30ms(圖9),嚴(yán)重影響后續(xù)偏移成像精度,導(dǎo)致較大的井震誤差,達(dá)不到頁巖氣勘探水平井鉆探對(duì)成像精度的要求。為此,在該工區(qū)采用本文方法進(jìn)行了近地表速度建模。輸入的最大偏移距分別為全炮檢距,4000,2000和1000m,網(wǎng)格大小分別為400,200,100和100m,每一個(gè)反演參數(shù)迭代25次。該單炮的初至波走時(shí)擬合平均誤差為11.2ms,降低到以前誤差值的37.33%,說明本文方法建模精度更高。圖10為利用傳統(tǒng)方法和本文方法得到的INLINE 2400線速度模型,可以看出,與傳統(tǒng)建模方法相比,本文建模方法獲得的速度模型近地表信息更加豐富,速度橫向變化更加劇烈。圖11為采用傳統(tǒng)方法和本文方法建模的偏移結(jié)果,可以看出,利用本文方法建立的近地表速度模型,偏移結(jié)果更加清晰,層位更加連續(xù)(圖11中黑圈所示)。

圖9 某單炮傳統(tǒng)方法建模(藍(lán)線)與本文方法建模(紅線)走時(shí)誤差對(duì)比

圖10 INLINE 2400線采用傳統(tǒng)方法(a)與本文方法(b)建模得到的速度模型對(duì)比

圖11 INLINE 2400線采用傳統(tǒng)方法(a)與本文方法(b)建模的偏移結(jié)果對(duì)比

4 結(jié)論

本文從正演和反演兩個(gè)方面對(duì)初至波走時(shí)射線層析進(jìn)行了改進(jìn)。正演方面引入雙線性插值技術(shù),有效提高了快速掃描法的計(jì)算精度,同時(shí)保持了該方法計(jì)算效率高的優(yōu)勢;反演方面采用多尺度策略,有效提高了改進(jìn)散射積分法層析反演精度。理論模型數(shù)據(jù)測試結(jié)果表明,本文的正、反演方法可以有效提高近地表建模精度。丁山工區(qū)近地表速度建模精度的提高有效改善了成像質(zhì)量,說明了本文速度建模思路的有效性和可行性。

本文方法基于射線理論,計(jì)算精度不易再進(jìn)一步提升,下一步應(yīng)逐步開發(fā)依托波動(dòng)理論的近地表建模方法,充分挖掘初至波振幅、相位等信息,進(jìn)一步提高近地表建模精度。