基于時(shí)空域自適應(yīng)高階有限差分的聲波疊前逆時(shí)偏移

基于時(shí)空域自適應(yīng)高階有限差分的聲波疊前逆時(shí)偏移

1疊前逆時(shí)偏移思想20世紀(jì)80年代的逆偏移技術(shù)開始了，該技術(shù)可以恢復(fù)和返回復(fù)雜路徑上的信號(hào)，以便能夠準(zhǔn)確地繪制復(fù)雜的地下結(jié)構(gòu)。因此，作為目前偏移方法中最精確的成像方法之一，反偏移是一種基于波動(dòng)理論的深度偏移方法。常用的波動(dòng)方程包括有限差分法、有限差分法、偽光譜法等。逆時(shí)位移思想的第一個(gè)起源是利用有限差分法求解波動(dòng)方程。疊前逆時(shí)偏移實(shí)現(xiàn)過程中的主要步驟是波場(chǎng)的正向延拓和反向延拓、成像條件以及成像后的低頻噪聲的消除,其中,波場(chǎng)正向延拓和反向延拓的本質(zhì)就是求解波動(dòng)方程.求解波動(dòng)方程是為了描述地震波場(chǎng)的傳播特征,理論上利用全波場(chǎng)的信息對(duì)地下構(gòu)造進(jìn)行準(zhǔn)確成像,所以精確、快速地求解波動(dòng)方程至關(guān)重要.Karazincir和Gerrard然而,常規(guī)的有限差分法求解波動(dòng)方程時(shí),其求空間導(dǎo)數(shù)時(shí)的差分系數(shù)一般都是采用傳統(tǒng)的高階有限差分系數(shù)2方法原理2.1數(shù)值模擬結(jié)果分析逆時(shí)偏移的核心問題之一是求解波動(dòng)方程.在縱波勘探中,一般采用聲波方程近似描述地震波在地下介質(zhì)中的傳播,二維聲波方程為其中,p為標(biāo)量波場(chǎng),v為速度,t為時(shí)間,x、z為空間坐標(biāo).一般情況下,計(jì)算時(shí)間二階導(dǎo)數(shù)時(shí),如果采用高階精度差分形式,則計(jì)算量會(huì)很大,并且很容易造成不穩(wěn)定其中,τ為時(shí)間采樣間隔,p方程(1)中空間二階導(dǎo)數(shù)的2M階精度差分形式可以寫為其中,a將式(2)和式(3)代入聲波方程(1)中,可以得到:逆時(shí)偏移的正向延拓就是差分方程正演模擬的過程,由式(4)可以推導(dǎo)出正向延拓的表達(dá)式:其中,f逆時(shí)偏移波場(chǎng)的反向延拓過程就是已知后一時(shí)刻的波場(chǎng)值來求前一時(shí)刻的波場(chǎng)值,具體計(jì)算過程與正向延拓類似,其計(jì)算表達(dá)式為L(zhǎng)iu和Sen其中,令f(θ)=cos為了計(jì)算和分析基于時(shí)空域有限差分的計(jì)算精度,我們分析差分的頻散,由式(7)定義頻散速度v用頻散速度與真實(shí)速度的比值來衡量二維聲波方程的頻散現(xiàn)象,則可以得到二維聲波方程的頻散分析表達(dá)式由于δ(θ)=δ(θ+π/2),所以δ(θ)是一個(gè)以π/2為周期的函數(shù).同時(shí)考慮到δ(θ)=δ(π/2-θ),所以在計(jì)算δ(θ)時(shí),θ的變化是從0到π/4.從式(11)中可以看出,如果δ=1,則模擬時(shí)不存在數(shù)值頻散,如果δ越是偏離1,則存在越大的數(shù)值頻散.由式(7)變形可得到:從式(12)可以看出,由于由于a其中M設(shè)計(jì)一個(gè)具體模型進(jìn)行精度分析:介質(zhì)的速度為3000m/s,模型網(wǎng)格大小為10m×10m,時(shí)間步長(zhǎng)為0.001s,空間導(dǎo)數(shù)精度階數(shù)2M分別為4、8、16、32.圖1(a,b)是波在傳播方向θ=π/8上,分別采用傳統(tǒng)有限差分和時(shí)空域有限差分?jǐn)?shù)值模擬時(shí)的頻散分析曲線圖.從圖1中可以看出:在相同參數(shù)情況下,時(shí)空域有限差分的頻散程度比傳統(tǒng)有限差分的頻散程度要小,即時(shí)空域有限差分的模擬精度要高于傳統(tǒng)有限差分如果空間導(dǎo)數(shù)的精度階數(shù)2M固定為16,而波的傳播角度θ分別取0、π/16、2π/16、3π/16、4π/16,其按順序?qū)?yīng)的標(biāo)識(shí)為angle1、angle2、angle3、angle4、angle5.圖2是有限差分?jǐn)?shù)值模擬在不同方向上頻散分析曲線圖.從圖2(a,b)的頻散分析曲線圖對(duì)比也可以看出:在精度階數(shù)相同時(shí),在波的各個(gè)傳播方向上,時(shí)空域有限差分的模擬精度比傳統(tǒng)有限差分高為比較傳統(tǒng)有限差分和基于時(shí)空域有限差分的穩(wěn)定性,由式(12)可令穩(wěn)定因子s為圖3是聲波方程有限差分在不同的精度階數(shù)(M分別為2、4、8、16)時(shí)隨數(shù)值模擬參數(shù)變化的穩(wěn)定性分析曲線圖.從圖3中可以看到,當(dāng)精度階數(shù)越高時(shí),越容易造成不穩(wěn)定,且時(shí)空域有限差分的穩(wěn)定因子與精度階數(shù)2M和數(shù)值模擬參數(shù)r有關(guān),而傳統(tǒng)有限差分的穩(wěn)定因子s只與精度階數(shù)2M有關(guān),同時(shí)對(duì)比兩類差分的穩(wěn)定性分析圖,可以看到,在保持求解方程穩(wěn)定的前提下(s≥r),時(shí)空域的有限差分對(duì)應(yīng)的r可取更大更廣的范圍.所以圖3中表明,時(shí)空域的有限差分的穩(wěn)性比傳統(tǒng)有限差分的穩(wěn)定性更好2.2自適應(yīng)差分算子長(zhǎng)度有限差分策略在給定一個(gè)地質(zhì)模型后,采用高階有限差分進(jìn)行數(shù)值模擬時(shí),通常都是采用固定長(zhǎng)度的空間差分算子來計(jì)算空間導(dǎo)數(shù).差分算子長(zhǎng)度的選擇必須同時(shí)滿足穩(wěn)定性、精度、計(jì)算量等的要求,而這些需要滿足的要求又與介質(zhì)模型的速度有關(guān).對(duì)于給定的速度模型,如果整體采用一個(gè)固定長(zhǎng)度的空間差分算子來進(jìn)行數(shù)值模擬,則這個(gè)固定的空間差分算子必須同時(shí)滿足模型中所有速度各自應(yīng)該滿足的條件.因此,固定的空間差分算子長(zhǎng)度一般會(huì)要求比較長(zhǎng),導(dǎo)致求解波動(dòng)方程耗費(fèi)大量計(jì)算時(shí)間.Liu和Sen求取自適應(yīng)差分算子長(zhǎng)度時(shí),先定義有限差分?jǐn)?shù)值模擬的誤差計(jì)算公式給定最大頻率f設(shè)計(jì)一個(gè)速度范圍為1500~5500m/s,網(wǎng)格大小為10m×10m,時(shí)間步長(zhǎng)為0.001s,其最大頻率為50Hz,而計(jì)算誤差范圍分別控制在102.3成像結(jié)果(s)成像條件也是逆時(shí)偏移處理過程中的一個(gè)重要步驟,成像條件的好壞直接影響著最終成像剖面的質(zhì)量.目前疊前逆時(shí)偏移中常用的成像條件主要有三種類型:激發(fā)時(shí)間成像條件其中,I(x,z)是點(diǎn)(x,z)的成像結(jié)果;s(x,z,t)表示時(shí)間域的正向延拓波場(chǎng);r(x,z,ω)表示時(shí)間域的反向延拓波場(chǎng);ε為無窮小實(shí)數(shù),其作用是為了保證計(jì)算的穩(wěn)定性.由于時(shí)間域計(jì)算相互關(guān)的速度相對(duì)比較慢,為了提高其計(jì)算效率,本文在實(shí)現(xiàn)過程中將波場(chǎng)轉(zhuǎn)換到頻率域計(jì)算互相關(guān),式(16)可以轉(zhuǎn)換為其中,ω表示角頻率,S(x,z,ω)表示頻率域的正向延拓波場(chǎng),R(x,z,ω)表示頻率域的反向延拓波場(chǎng),ue494(x,z,ω)表示S(x,z,ω)的共軛.2.4拉普拉斯算子去噪疊前逆時(shí)偏移的過程中,不可避免地會(huì)產(chǎn)生低頻噪聲,所以要得到最終的清晰成像結(jié)果,還要對(duì)偏移的初始剖面進(jìn)行去噪處理.拉普拉斯算子去噪方法簡(jiǎn)單且易操作,一般將其寫成二階差分形式,對(duì)逆時(shí)偏移成像結(jié)果進(jìn)行濾波,對(duì)低頻噪聲去除效果比較明顯拉普拉斯算子表達(dá)式為其在波數(shù)域表示為其中,k應(yīng)用余弦定理,由波數(shù)域矢量計(jì)算方法可以得到:式(20)中θ是入射角,k同時(shí),從式(21)也可以看出,在拉普拉斯算子去噪的同時(shí)也會(huì)損害成像的有用信息,如分子中的頻率項(xiàng)應(yīng)該在成像前的數(shù)據(jù)上進(jìn)行補(bǔ)償.本文采用的是時(shí)間域積分方法對(duì)成像前數(shù)據(jù)進(jìn)行補(bǔ)償,利用的關(guān)系表達(dá)式為其中,g(t)是時(shí)間域原始信號(hào).3計(jì)算值的示例3.1介質(zhì)模型數(shù)值模擬為了進(jìn)一步驗(yàn)證時(shí)空域有限差分相對(duì)傳統(tǒng)有限差分有更好的數(shù)值模擬效果,設(shè)計(jì)一個(gè)大小為201×201(網(wǎng)格數(shù))、網(wǎng)格大小為10m×10m、速度為3000m/s的均勻介質(zhì)模型.數(shù)值模擬時(shí),震源采用一個(gè)周期的正弦子波,且其頻率為50Hz,時(shí)間采樣間隔為0.001s,接收的記錄長(zhǎng)度為2s.震源位于(1000m,1000m)處,分別采用傳統(tǒng)有限差分和時(shí)空域有限差分,以時(shí)間上二階精度、空間上二十階精度進(jìn)行數(shù)值模擬,在數(shù)值模擬邊界處理時(shí),采用完全匹配層(PML)吸收邊界條件圖5(a,b)分別是采用傳統(tǒng)有限差分和時(shí)空域有限差分模擬得到的0.3s時(shí)刻的波場(chǎng)快照.從快照波前面的形狀來看,傳統(tǒng)有限差數(shù)值分模擬得到的波場(chǎng)快照的波形發(fā)生了畸變,有較強(qiáng)頻散,而時(shí)空域有限差分模擬得到的波場(chǎng)快照的波形保持得比較好,頻散相對(duì)較弱.圖6(a,b)分別是采用傳統(tǒng)有限差分和時(shí)空域有限差分模擬得到的單炮記錄,從(a)、(b)圖對(duì)比也可以看出,時(shí)空域有限差分模擬得到記錄的頻散更小,波形保持得更好.因此介質(zhì)模型中的數(shù)值模擬結(jié)果表明,時(shí)空域有限差分?jǐn)?shù)值模擬精度要高于傳統(tǒng)有限差分的數(shù)值模擬精度3.2數(shù)值模擬對(duì)比斷層模型如圖7所示.模型大小為301×201(網(wǎng)格數(shù)),網(wǎng)格大小為10m×10m,斷層上、下部分的速度分別為1600m/s和2300m/s.接收記錄長(zhǎng)度為3s,時(shí)間步長(zhǎng)為0.001s,震源采用一個(gè)周期的25Hz的正弦子波,一共58炮,炮點(diǎn)位置在模型90m的深度面從距離最左邊100m處開始向右移動(dòng),炮間距為50m,每一炮都是全地表等距(10m)接收.在數(shù)值模擬和逆時(shí)偏移中,都以時(shí)間的二階精度、空間的十二階精度有限差分求解方程.圖8(a,b)是震源在(1500m,90m)時(shí)分別采用傳統(tǒng)有限差分和時(shí)空域有限差分模擬得到的1.5s時(shí)刻的波場(chǎng)快照.從圖8中對(duì)比傳統(tǒng)有限差分和時(shí)空域有限差分?jǐn)?shù)值模擬的結(jié)果表明,時(shí)空域有限差分模擬得到的波場(chǎng)的頻散更小,精度更高.圖9和圖10分別是采用傳統(tǒng)有限差分和時(shí)空域有限差分對(duì)斷層速度模型進(jìn)行逆時(shí)偏移而得到的剖面.對(duì)比去噪前后的偏移剖面,可以看到,高通濾波和拉普拉斯算子去噪都能有效地去除低頻噪聲,但是高通濾波在去除噪聲的同時(shí),也把斷層的垂直部分濾掉了,即對(duì)陡傾角信號(hào)有損害;而拉普拉斯算子去噪能在去除噪聲的同時(shí)保留好斷層的垂直部分,所以拉普拉斯算子去噪在逆時(shí)偏移的成像噪聲消除中,優(yōu)于傳統(tǒng)高通濾波方法.從圖9和圖10中也能看到,逆時(shí)偏移使水平地層和垂直斷層都能較好地成像,且與模型相比,界面成像位置準(zhǔn)確,尤其是時(shí)空域有限差分逆時(shí)偏移得到的最終剖面,斷層的垂直部分成像更清楚、更準(zhǔn)確.3.3自適應(yīng)差分算子長(zhǎng)度有限差分的計(jì)算效率圖11是Sigsbee2A偏移速度模型.該模型的特點(diǎn)是擁有較為復(fù)雜的鹽丘形狀以及頂部構(gòu)造,因此經(jīng)常被用來驗(yàn)證逆時(shí)偏移的效果以及比較噪聲去除效果.模型大小為3201×1201(網(wǎng)格數(shù)),網(wǎng)格大小為7.62m×7.62m,速度范圍為1437~4511m/s.接收記錄的時(shí)間長(zhǎng)度為12s,時(shí)間步長(zhǎng)為0.0005s,震源采用一個(gè)周期的正弦子波,其頻率為30Hz.一共200炮,炮點(diǎn)位置是從距離模型最左邊1524m處開始向右移動(dòng),炮間距為106.68m,每一炮都是全地表等距(7.62m)接收,震源和接收點(diǎn)位于同一深度面15.24m.將求解波動(dòng)方程的誤差控制在10數(shù)值模擬和逆時(shí)偏移的過程中,基于時(shí)空域有限差分方法,分別采用固定差分算子長(zhǎng)度和自適應(yīng)差分算子長(zhǎng)度求解波動(dòng)方程.將自適應(yīng)差分算子長(zhǎng)度應(yīng)用到偏移的正向延拓和反向延拓計(jì)算時(shí),從理論上分析,可以看到自適應(yīng)差分算子長(zhǎng)度平均上要短一些,能有效改進(jìn)偏移效率.我們?cè)谔幚砥脚_(tái)為ThinkPadR400筆記本電腦(Core2雙核,2.0GHz)上做了一個(gè)計(jì)算效率測(cè)試:統(tǒng)計(jì)5個(gè)炮點(diǎn)正向延拓的耗時(shí),如表1所示,從表中統(tǒng)計(jì)可以看出,自適應(yīng)差分算子長(zhǎng)度有限差分的計(jì)算效率相對(duì)固定差分算子長(zhǎng)度有限差分可以節(jié)約用時(shí)33%左右;統(tǒng)計(jì)了5個(gè)炮點(diǎn)各自單炮逆時(shí)偏移成像整個(gè)過程的耗時(shí),如表2所示,從表中也可以看到,在逆時(shí)偏移的過程中,自適應(yīng)差分算子長(zhǎng)度有限差分策略能提高效率23%左右.圖13(a,b)是分別采用兩類算子長(zhǎng)度正演模擬得到的波場(chǎng)快照(4s時(shí)刻).從圖13可以看出,采用兩類算子長(zhǎng)度正演模擬的結(jié)果基本一致.因此采用自適應(yīng)有限差分法,是一種在能保證精度要求的前提下,提高計(jì)算效率的策略圖14a是Sigsbee2A速度模型采用時(shí)空域自適應(yīng)有限差分方法逆時(shí)偏移的結(jié)果.雖然該速度模型對(duì)速度場(chǎng)做了平滑,但在鹽丘的邊界處還是存在著強(qiáng)烈的速度間斷,會(huì)對(duì)逆時(shí)偏移帶來較強(qiáng)烈的低頻噪聲,可以看到很強(qiáng)的低頻噪聲集中在鹽丘邊界反射處,如圖14a所示,因此難以直接分析逆時(shí)偏移的效果.圖14b是采用高通濾波對(duì)逆時(shí)偏移結(jié)果進(jìn)行去噪后的結(jié)果,從該結(jié)果中可以看到,經(jīng)過高通濾波,低頻噪聲得到了有效的壓制,逆時(shí)偏移的效果也顯現(xiàn)出來了,但還是沒能很干凈地去除噪聲,損害陡傾角信息.圖14c是應(yīng)用拉普拉斯算子去噪后的逆時(shí)偏移剖面,可以看到低頻噪聲得到了很好的消除,與模型對(duì)比,偏移成像結(jié)果比較好,成像位置準(zhǔn)確,特別是在鹽丘邊界與鹽丘底部成像清晰.由于受觀測(cè)系統(tǒng)和地下的地質(zhì)構(gòu)造的影響,可能會(huì)導(dǎo)致在巖丘底部有些區(qū)域照明度不夠,而要克服此問題,則在偏移過程中需要進(jìn)行照明度補(bǔ)償3.4邊界條件、自適應(yīng)差分和自適應(yīng)差分逆時(shí)偏移的應(yīng)用效果圖15是Marmousi速度模型,對(duì)其進(jìn)行重新采樣后,模型大小為921×300(網(wǎng)格數(shù)),網(wǎng)格大小為10m×10m.該模型的速度范圍1500.0~5500.0m/s,總層數(shù)超過100層,有大量的反射界面,并含有一些陡傾角、逆沖斷層、角度不整合面、地層隆起以及強(qiáng)烈的橫向、縱向的速度變化.震源采用一個(gè)周期的正弦子波,其頻率為30Hz,一共240炮,炮點(diǎn)位置在從距離最左邊的2500m處開始向右移動(dòng),炮間距為20m,右邊放炮,左邊接收,每一炮都是120道等距(20m)接收.接收記錄長(zhǎng)度為3s,時(shí)間步長(zhǎng)為0.001s.將求解波動(dòng)方程的誤差控制在10分別采用固定差分算子長(zhǎng)度和自適應(yīng)差分算子長(zhǎng)度的時(shí)空域有限差分法求解波動(dòng)方程和進(jìn)行逆時(shí)偏移成像.表3是統(tǒng)計(jì)兩類差分算子長(zhǎng)度進(jìn)行逆時(shí)偏移的計(jì)算效率,可以看到,自適應(yīng)差分算子長(zhǎng)度有限差分逆時(shí)偏移能節(jié)時(shí)24%左右,有效提高了計(jì)算效率.圖17a是基于時(shí)空域自適應(yīng)高階有限差分逆時(shí)偏移后的原始結(jié)果,由于Marmousi模型存在大量的反射界面,所以在逆時(shí)偏移之后,其整個(gè)剖面上存在較強(qiáng)的低頻噪聲.圖17b是應(yīng)用拉普拉斯算子去噪方法對(duì)圖17a進(jìn)行處理后的成像剖面,低頻噪聲得到了比較好的壓制.從圖17b中也可看出,逆時(shí)偏移對(duì)模型的主要構(gòu)造都進(jìn)行了較好的成像,特別是在中淺部的大傾角構(gòu)造、斷裂及深部的角度不整合地層、目的層(地層隆起)均得到了比較好的成像.4數(shù)值模擬結(jié)果分析本文采用基于時(shí)空域有限差分的差分系數(shù)求解聲波方程,并分析了其數(shù)值頻散和穩(wěn)