基于密度擾動的點擴散函數(shù)計算與三維深度域模擬成像

段偉國, 毛偉建, 張慶臣, 石星辰, 張建磊

1 中國科學(xué)院精密測量科學(xué)與技術(shù)創(chuàng)新研究院計算與勘探地球物理研究中心, 大地測量與地球動力學(xué)國家重點實驗室, 武漢 430077 2 中國科學(xué)院大學(xué), 北京 100049 3 中國石油集團公司東方地球物理公司物探技術(shù)研究中心, 河北涿州 072751

0 引言

隨著計算機處理性能的不斷提升,疊前深度偏移技術(shù)已經(jīng)逐步走向成熟并積累了大量的深度域偏移數(shù)據(jù).如何應(yīng)用這些資料直接在深度域開展地震解釋和儲層參數(shù)反演,在實現(xiàn)結(jié)構(gòu)成像的同時獲得巖石物性信息已經(jīng)成為工業(yè)界的迫切需求.然而,深度域地震子波的提取(馮瑋等,2017)和深度域合成地震記錄的困難嚴重限制了深度域反演方法的發(fā)展.這是因為地震反演中模擬偏移成像的傳統(tǒng)方法是在時間域進行一維褶積.這種方法基于水平層狀假設(shè)和對時間偏移數(shù)據(jù)的近似,由于地下介質(zhì)的速度是深度的函數(shù),故深度域“子波”是隨空間變化的,因而難以將時間域合成地震記錄的方法直接應(yīng)用到深度域.

針對上述問題,在傳統(tǒng)穩(wěn)態(tài)褶積模型的基礎(chǔ)上逐漸發(fā)展了非穩(wěn)態(tài)褶積模型(高靜懷等,2009),目前通行的在深度域合成地震數(shù)據(jù)和資料解釋的方法包括：基于深變子波的方法(張雪建等,2000; 何惺華,2004)、基于偽深度轉(zhuǎn)化的方法(林伯香等,2006; 周東勇等,2020)、基于時深轉(zhuǎn)換的方法(周賞等,2017)、基于多屬性變換的方法(張靜等,2010)、基于神經(jīng)網(wǎng)絡(luò)數(shù)據(jù)驅(qū)動的方法(崔鳳林和張向君,2005)、深度域波形反演(Symes and Carazzone,1991; 徐嘉亮等,2015)和基于點擴散函數(shù)的方法(Fletcher et al.,2012).PSF是成像系統(tǒng)對單位強度點散射體的響應(yīng),又稱偏移格林函數(shù)(Schuster and Hu,2000),分辨率函數(shù)(Gelius et al.,2002)和照明響應(yīng)函數(shù)(Xie et al.,2006),是對Hessian 算子的一種近似(Jiang and Zhang,2019).三維空間變化的PSF比一維子波包含了更多的動力學(xué)信息,通過PSF可以直接在深度域合成地震記錄,避免了時深轉(zhuǎn)換的過程,突破了常規(guī)時間域/深度域一維褶積對水平層狀介質(zhì)假設(shè)的限制,能夠更好的處理三維復(fù)雜介質(zhì)的問題.近年來,基于PSF的深度域反演逐漸成為研究的重點和熱點,國內(nèi)外學(xué)者將其應(yīng)用到照明補償(Letki et al.,2015; Fletcher et al.,2016; Dai et al.,2019)、Q因子補償(Cavalca et al.,2015；王雪君等,2019; 姚振岸等,2019)、去鬼波(Caprioli et al.,2014)、同時震源去除串擾噪聲(張攀和毛偉建,2018a)、以及反射系數(shù)和波阻抗反演(He et al.,2019; Leon et al.,2019).

關(guān)于PSF的計算大致可分為三種方法: 1)利用遠場的閉合表達式進行計算(Chen and Schuster,1999); 2)基于射線追蹤的照明分析方式(Lecomte and Gelius,1998); 3)結(jié)合正演和偏移算子進行計算(Toxopeus et al.,2008).在實現(xiàn)PSF計算時,張金陵等(2019)只考慮了一維的情況,可以看作對基于深變子波方法的另一種表述.Fletcher等(2012)基于方法3的思路采用兩次正演的方式計算了二維聲波PSF,該策略以及后來發(fā)展的二維黏聲PSF(姚振岸等,2019)中都有很強的結(jié)構(gòu)信息污染,并且這些文章中都沒有給出三維PSF的計算示例.針對上述問題,本文提出了基于密度擾動的兩步正演法來計算點擴散函數(shù),并對擾動強度對結(jié)構(gòu)信息的污染、成像匹配中的相位校正問題、以及空間褶積的具體實現(xiàn)細節(jié)進行了深入的研究,得到了二維/三維深度域模擬成像結(jié)果,并通過模型實驗驗證了方法的可靠性.

1 原理

1.1 基于點擴散函數(shù)的深度域反演原理

定義地球反射率參考模型m和正演算子A,則地震數(shù)據(jù)可以表示為d=Am.假設(shè)觀測數(shù)據(jù)為dobs,通過定義目標函數(shù)Sm=‖Am-dobs‖2可以得到最小二乘意義下的反射率模型為

(1)

這里AT為偏移算子,是正演算子A的伴隨算子,mmig是由觀測數(shù)據(jù)得到的偏移成像,H=ATA是Hessian算子.Hessian算子是對成像系統(tǒng)質(zhì)量的度量,它包括了正演算子和偏移算子的信息,同樣它也可以捕捉振幅隨幾何擴散、傳播損耗和采集系統(tǒng)幾何形狀的變化,完整的Hessian算子代表了一個完備的成像系統(tǒng)(任浩然等,2013).實際上,我們很難直接求解方程(1),因為它需要求解龐大的Hessian矩陣的逆.因此,一種更為可行的方法是對線性方程(2)進行迭代求解得到方程的最小二乘解(周華敏等,2014; 劉玉金和李振春,2015; 張攀和毛偉建,2018b; 胡勇等,2021)：

(2)

迭代的過程如下: 將觀測數(shù)據(jù)進行偏移得到反射率成像,再對反射率成像進行反偏移得到模擬數(shù)據(jù),將模擬數(shù)據(jù)減去觀測數(shù)據(jù)再對殘差進行偏移用來更新反射率模型,這就構(gòu)成了一次迭代.通常,該方法要經(jīng)過15次左右的迭代才能得到一個可接受的解(Nemeth et al.,1999).由于每次迭代需要進行偏移和反偏移操作,這在實際生產(chǎn)中會產(chǎn)生很大的計算量.

PSF 反演方法將上述數(shù)據(jù)域的迭代過程轉(zhuǎn)化到深度域來進行,從而避免了偏移和反偏移的迭代,降低了計算成本.對于一個點散射體產(chǎn)生的地震數(shù)據(jù)dsc=Amsc,其中A是單位散射點的正演算子,msc是單位散射點模型.將偏移算子作用于dsc,可以得到點散射體的偏移成像結(jié)果(也就是PSF)：

fpsf=ATdsc=ATAmsc≈H,

(3)

這里偏移算子AT是正演算子A的伴隨算子,并且我們需要讓正演算子A盡可能地接近地震波的真實傳播過程.單個散射點的響應(yīng)是以該散射點為中心以ATA為算子的體積分布.因此,通過PSF與反射率模型之間進行空間褶積可以得到深度域模擬成像(Toxopeus,2008)：

(4)

1.2 速度和密度擾動對散射場的影響

根據(jù)擾動理論,可以通過設(shè)置速度或密度擾動點來產(chǎn)生散射場,進而求取點擴散函數(shù).詳細推導(dǎo)如下: 無源時壓力場p滿足以下方程(Beylkin and Burridge,1990):

(5)

其中s表示震源位置,x是散射點位置,t表示時間.又因為σ,κ和p可以寫作如下形式:

σ(x)=σ0(x)+σ′(x),

κ(x)=κ0(x)+κ′(x),

p(s,x,t)=p0(s,x,t)+p′(s,x,t),

(6)

這里σ0和κ0為背景比體積(specific volume)和背景壓縮系數(shù)(compressibility),σ′和κ′是它們對應(yīng)的擾動量,p0是背景波場,p′為散射波場.將公式(6)帶入公式(5)并應(yīng)用

(7)

之后將無擾動項和擾動項波場分別放在等號兩邊，可得到

(8)

這里D是由封閉曲面?D圍成的區(qū)域,在這個區(qū)域中存在滿足方程(8)的擾動量σ′和κ′,*t表示時間卷積.為了得到單散射波場的最奇異項,我們考慮格林函數(shù)及其導(dǎo)數(shù)的最奇異部分：

(10)

(11)

(12)

將公式(10)—(13)帶入公式(9)可以得到

(14)

(15)

(16)

(17)

這里θ=θ(x,s,r)是震源s、散射點x和檢波器r所夾的散射張角(如圖1所示).

圖1 震源-散射點-檢波器之間的散射張角示意圖

通過公式(17)和(6),將公式(14)中括號內(nèi)的表達式寫成如下形式：

(18)

將公式(18)帶入(14)最終得到單散射波場的最奇異項表達式：

(19)

其中

(20)

從公式(19)可以看出,散射場和f(x,θ)存在線性關(guān)系,f是積分中惟一依賴擾動參數(shù)σ和κ的項,我們把f稱為散射勢.

散射勢f和角度θ有關(guān),它表明了振幅輻射模式(amplitude radiation patterns).因為比體積是密度ρ的倒數(shù)σ=1/ρ,壓縮系數(shù)是體積模量的倒數(shù)κ=σ/c2.結(jié)合公式(6)可以將公式(20)寫成由密度和速度參數(shù)表示的形式：

(21)

當擾動點處密度ρ等于背景密度ρ0,即ρ=ρ0時,公式(21)可以化簡為僅與速度項有關(guān)的形式

(22)

當擾動點處的速度c等于背景速度c0,即c=c0時,公式(21)可以化簡為僅與密度有關(guān)的形式

(23)

1.3 點擴散函數(shù)的計算

有效地生成準確的PSF是應(yīng)用該方法進行深度域反演的關(guān)鍵.本文采用結(jié)合正演算子A和偏移算子AT的方式計算PSF.根據(jù)公式(3)，有

psf=ATAmsc=ATA(m-m0)=AT(Am-Am0),

(24)

其中m0為地球模型,m為在地球模型的基礎(chǔ)上加入單位強度擾動點的含散射點模型.通過兩次有限差分正演分別產(chǎn)生地震數(shù)據(jù)Am和Am0.有限差分正演是一個多重散射算子,在一階近似的情況下,通過兩次正演數(shù)據(jù)的差A(yù)m-Am0來模擬單位強度散射點的波場Amsc.之后對散射場數(shù)據(jù)Amsc進行深度偏移ATAmsc即可得到PSF.

本文提出的PSF計算方法基于Born近似理論,而常規(guī)的反射理論基于Kirchhoff近似.這兩種假設(shè)之間存在一定的差別: Born近似所描述的是散射場與散射勢之間的關(guān)系,而Kirchhoff近似則是反射場與反射系數(shù)模型的關(guān)系(Jaramillo and Bleistein,1999).基于Born近似假設(shè)條件獲得的成像參數(shù)與基于Kirchhoff近似假設(shè)條件獲得的成像參數(shù)之間的差異等價于一次空間濾波(Beylkin and Burridge,1990).因此,在對散射數(shù)據(jù)進行深度偏移之前,需要先對其進行90°相位校正,以此來保證Born近似條件下PSF與反射系數(shù)褶積得到的模擬成像剖面和實際偏移剖面的相位一致.下面通過Kirchhoff面積分公式改寫為類似Born體積分公式的形式,對二者進行比較來證明上述結(jié)論.

對于Kirchhoff積分公式(Beylkin and Burridge,1990)：

(25)

R(x,θs)=f(x,2θs)/4cos2θs,

(26)

由于θs=θ/2,公式(26)可以改寫成

R(x,θ/2)=f(x,θ)/4cos2(θ/2).

(27)

又根據(jù)(Beylkin and Burridge,1990)

(28)

利用公式(27)和(28)可以將公式(25)表示成

定義γ(x)為積分面的奇異函數(shù)(singular function of surface),它表示以該積分面為支撐集的狄拉克函數(shù)(Bleistein,1987)：

(30)

使得對于給定的函數(shù)g,滿足

(31)

利用公式(31)可以將公式(29)轉(zhuǎn)化為體積分形式：

將公式(19)轉(zhuǎn)化到頻率域可以得到與(32)式對應(yīng)的基于Born 近似的單散射場積分公式：

(33)

通過對比Kirchhoff積分公式(32)和Born積分公式(33)可以得到

(34)

公式(34)中iw項表明與基于Kirchhoff近似得到的反射場相比,基于Born近似得到的散射場存在90°的相移(Yang and Zhang,2019).因此在利用本文方法得到散射場后,需要對其進行90°的相位校正,再進行偏移操作,否則利用PSF得到的模擬成像剖面在相位上無法與實際偏移成像剖面相匹配.

在1.2節(jié)中,已經(jīng)推導(dǎo)了單個速度或密度擾動點對散射場的影響.然而正演模擬空間中每一個點的單散射波場需要很大的計算量,因此在實際應(yīng)用中不是每次只設(shè)置一個速度或密度散射點然后逐個計算PSF,而是假設(shè)在一定空間范圍內(nèi)PSF變化不大并且是線性變化的.這樣就可以同時設(shè)置多個散射點進行正演,以一種經(jīng)濟有效的方式計算PSF.雖然使用兩步正演法會在散射點之間產(chǎn)生一些噪聲,但這些噪聲屬于二階效應(yīng),偏移后呈現(xiàn)出非相干性,且相對較弱.

在設(shè)置擾動時需要保證所有散射點的散射勢f大小相等.因此,既可以根據(jù)公式(21)同時存在速度和密度擾動,也可以根據(jù)公式(22)或(23)只存在速度或密度一種擾動,這三種方式在理論上都是可行的.考慮到實現(xiàn)的可操作性,一般只需要設(shè)置一種擾動.為了更進一步分析速度擾動和密度擾動對有限差分正演中數(shù)值計算準確性的影響,本研究進行了如下測試.

首先分析速度擾動和密度擾動產(chǎn)生的散射場的特性.如圖2所示，在網(wǎng)格大小為401×401、網(wǎng)格間距為10 m×10 m的均勻模型中心分別設(shè)置速度和密度擾動點.該模型的背景速度為2000 m·s-1,背景密度為2000 kg·m-3,在擾動點處設(shè)置3000 m·s-1的速度擾動或3000 kg·m-3的密度擾動.之后分別進行有限差分正演,炮點位置設(shè)在(0,0)處,震源子波為主頻25 Hz的雷克子波,滿排列接收,進而對比由速度擾動和密度擾動產(chǎn)生的散射場.

圖2 速度/密度擾動示意圖

圖3a為由速度擾動產(chǎn)生的散射場記錄,并抽取了其中3個地震道進行展示.對比這3個地震道并結(jié)合公式(22)可知: 速度擾動產(chǎn)生的散射場是各項同性的,不同角度的散射場振幅大小相等,并且散射場不是零相位的.對圖3a中的地震記錄做90°相位校正可以得到圖3b,可以看到經(jīng)過校正后,散射場變成零相位的.

圖3 (a) 速度擾動產(chǎn)生的散射場; (b) 90°相位校正的速度散射場

同樣地,我們可以得到由密度擾動產(chǎn)生的散射場記錄(見圖4a),對比圖4a中3個地震道并且結(jié)合公式(23)可知: 由密度擾動產(chǎn)生的散射場的振幅是隨散射張角變化的,該角度與波場入射方向有關(guān)，而入射方向受背景速度的影響.并且密度散射場也存在90°的相移,對其做90°相位校正后,可以得到圖4b,此時散射場是零相位的.

圖4 (a) 密度擾動產(chǎn)生的散射場; (b) 90°相位校正的密度散射場

為了測試擾動點對波場傳播的影響,設(shè)定網(wǎng)格大小為1001×1001,網(wǎng)格間距為10 m×10 m的雙層反射模型.分別對三次有限差分正演的數(shù)值實驗進行對比: 第一次正演采用無擾動點的背景模型,如圖5a所示,背景模型上層速度為2000 m·s-1,密度為2000 kg·m-3,下層速度為3000 m·s-1,密度為3000 kg·m-3; 第二次正演在第一次正演背景模型的基礎(chǔ)上設(shè)置間距為200 m×200 m的速度擾動點,擾動點處速度大小是背景速度的兩倍; 第三次正演在第一次正演背景模型的基礎(chǔ)上設(shè)置間距為200 m×200 m的密度擾動點,擾動點處密度大小為背景密度的兩倍(見圖5b).三次正演實驗的采集裝置相同,震源為位于(0,0)處主頻25 Hz的爆炸震源,滿排列接收.通過對比三次正演接收到的地震數(shù)據(jù)來分析擾動點對波場傳播的影響.

如圖6(a—c)所示,分別為三次正演實驗接收到的地震記錄,其中(b)圖和(c)圖中能量較弱的部分為散射場信號.抽取檢波點位置為4000 m、5000 m、6000 m處,采樣時間為9.2～9.7 s的地震道進行對比(見圖7).其中黑色線為無擾動模型的地震記錄,紅色線為速度擾動模型的地震記錄,綠色線為密度擾動模型的地震記錄.可以看出不論是密度擾動還是速度擾動,都會對經(jīng)過該擾動點的波場的振幅產(chǎn)生一定的影響.在本模型測試中,密度擾動點對整個波場的走時幾乎沒有任何影響,速度擾動點使得接收到的波場約有9 ms的走時誤差,且由于該誤差的存在使得采用波場相減計算散射場時,對后續(xù)得到的PSF產(chǎn)生更多的結(jié)構(gòu)污染.同時,考慮到速度模型對有限差分正演是至關(guān)重要的,對速度值的改變可能會使模型不再滿足正演的穩(wěn)定性條件,必須要調(diào)整原模型網(wǎng)格間距或者采樣步長重新進行計算.雖然密度擾動產(chǎn)生的散射場振幅和散射張角有關(guān),對于疊前PSF的計算,可以通過對不同的角道集除以1+cosθ來消除散射張角對振幅的影響.綜上,由于速度擾動對經(jīng)過該擾動點的波場影響更大,以及速度模型受限于有限差分正演的穩(wěn)定性,因此本研究認為采用密度擾動要優(yōu)于速度擾動.

結(jié)合上述分析,可以給出計算PSF的具體方式: 將模型劃分為多個相同大小的窗口,一般來講窗口長度為2～3倍的波長,每個窗口的中心放置一個密度散射點.每個散射點的密度擾動為δρ=ρ-ρ0.因為散射場是散射勢的近似線性函數(shù),對單個散射點來說ρ可以是比ρ0大的任意值,但對于多散射點模型來說,理論上我們只需要保證每個散射點δρ/ρ是相同的,以此確保由不同散射點產(chǎn)生的PSF擁有相同的散射勢.之后對經(jīng)過90°相位校正后的散射場采用和常規(guī)偏移成像相同的偏移算子進行偏移就可以得到這些散射點的PSF.然后,對擾動點處的PSF進行空間插值就可以得到任意位置的PSF.

當模型結(jié)構(gòu)復(fù)雜時,通過兩次正演波場相減得到的散射場中可能仍含有很強的結(jié)構(gòu)信息,進而導(dǎo)致偏移后生成的PSF被結(jié)構(gòu)信息嚴重污染.Fletcher和Cavalca(2018)選擇將被結(jié)構(gòu)信息污染的PSF進行剔除,之后用被剔除點附近的PSF作為替代再進行下一步的處理,防止被污染的PSF影響處理效果; 姚振岸等(2019)則沒有進一步考慮被污染的PSF對結(jié)果的影響.本文通過測試發(fā)現(xiàn),造成PSF被結(jié)構(gòu)信息污染的很重要一個原因在于擾動強度過大,當減小擾動強度時,結(jié)構(gòu)信息的污染效應(yīng)可以得到很大程度的緩解.這是因為散射點產(chǎn)生的波場在遇到模型界面時會發(fā)生二次反射,當擾動強度過大時接收到的高階能量變強,導(dǎo)致PSF被結(jié)構(gòu)信息污染,通過減小擾動強度可以降低二次反射的能量,進而減弱結(jié)構(gòu)信息對PSF的污染.但是這并不意味著擾動強度越小越好,當擾動強度過小導(dǎo)致一階散射場的能量很弱時,又會使散射波場容易受到微小噪聲的干擾,降低PSF的信噪比.因此擾動強度的設(shè)置取決于模型的復(fù)雜程度,很難確定一個對所有模型都適用的最優(yōu)擾動值.當模型結(jié)構(gòu)復(fù)雜而導(dǎo)致PSF中含有很強的結(jié)構(gòu)信息污染時,可以通過適當?shù)臏p小擾動強度進行改善.經(jīng)驗表明當δρ=0.1ρ0時可以很好地避免結(jié)構(gòu)信息的污染,同時又不會受到散射體之間的噪聲干擾(詳見模型測試中圖14和圖15).

1.4 深度域模擬成像

如果將偏移成像看作一個模糊的反射率模型(Yu et al.,2006),模糊核就是隨空間變化的PSF.所以當我們有了PSF之后,深度域的模擬成像就可以用公式(4)所示的空間褶積來表示.假定在一定的區(qū)域內(nèi)PSF是隨空間均勻變化的,通過對密度擾動點處的PSF進行空間插值便可以得到任意位置的點擴散函數(shù).這里我們采取了線性插值的方法,并在插值過程中進行了跨邊界處理(Fletcher and Cavalca,2018).

線性插值有一階、二階、三階,對應(yīng)為單線性插值、雙線性插值和三線性插值.三者皆為線性插值,不同之處在于單線性插值可以對2點之間的任一點的插值,雙線性插值可以對4點形成的正方形中間中任一點進行插值,三線性插值可以對8點形成的立方體內(nèi)任一點進行插值.

如圖8(a—c)所示,單線性插值中已知x1,x2兩個點的值,對于兩點之間的任一點x,有比例因子a=(x-x1)/(x2-x1),則x=a(x2-x)+(1-a)(x-x1).在雙線性插值中,對于4點之間任意點(x,y),可以求出x方向和y方向的比例因子ax,ay然后根據(jù)ax在x方向進行單線插值得到如圖8b中(x,y1)和(x,y2)兩個點處的值,再根據(jù)ay在y方向進行單線性插值得到插值點的值.在三線性插值中,對于8點之間的任意點(x,y,z),可以求出x、y、z方向的比例因子ax、ay、az,然后在任意兩個相對平面上進行雙線性插值,之后再對另一個方向進行單線性插值便可以得到點(x,y,z)處的值.

為了驗證插值方法的準確性,本文進行了PSF插值測試，結(jié)果如圖9所示.用立方體空間上8個角處的PSF對其中心位置的PSF進行插值.圖9a為立方體中用來插值的8個點處的PSF,圖9b為待插值點處的真實PSF,圖9c為插值得到的PSF,將插值結(jié)果和真實值進行比較可以看出,空間線性插值較好地還原了真實PSF的大小和形態(tài).

反射率模型與PSF進行空間褶積可以看作將反射率模型中每一個點的反射率用對應(yīng)位置處的PSF進行模糊的過程.將每個位置的反射率與對應(yīng)的PSF相乘,并將重疊部分進行相加求和便可以得到深度域模擬成像.在利用疊后PSF進行模擬成像時,可以采用公式(35)計算反射率模型：

(35)

在實現(xiàn)過程中，作者采用了動態(tài)插值的方式,即在對PSF進行插值的同時進行褶積運算,這樣就避免了存儲每個插值點的PSF,有效節(jié)省了內(nèi)存空間.并且這種動態(tài)處理的方式便于對局部目標區(qū)域進行模擬成像,為實時成像提供了可能.

進行矩陣的褶積運算時,往往要考慮褶積核是否需要翻轉(zhuǎn)的問題.在數(shù)學(xué)中嚴格的褶積定義要求褶積核圍繞中心進行180°翻轉(zhuǎn)之后再進行褶積,這種褶積方式廣泛應(yīng)用于信號處理,以及求兩個隨機變量的分布等問題中.但是在其他領(lǐng)域,特別是卷積神經(jīng)網(wǎng)絡(luò)中,不需要對褶積核進行翻轉(zhuǎn)操作.因為這種褶積是為了提取圖像的特征,其本質(zhì)上是一種空間濾波,只借鑒了前面加權(quán)求和的思想.上述嚴格意義上兩種不同的運算往往被不加區(qū)別地稱為卷積或褶積(convolution),因此在涉及(尤其是使用開源程序)褶積或反褶積的過程中,褶積核是否需要翻轉(zhuǎn)可能影響結(jié)果的好壞.點擴散函數(shù)為單位強度散射點的成像響應(yīng),空間褶積是將整個反射率模型用點擴散函數(shù)進行加權(quán)求和,所以在這個過程中不需要將點擴散函數(shù)進行翻轉(zhuǎn).

2 模型測試

為了驗證本文提出的基于密度擾動的兩步正演法計算點擴散函數(shù)的準確性,以及利用點擴散函數(shù)進行深度域模擬成像的效果,進行以下模型測試.

首先以三維SEG/EAGE Salt模型的一部分為例進行測試,如圖10a所示為速度模型,模型網(wǎng)格大小為200×200×150,網(wǎng)格間距均為20 m.正演數(shù)據(jù)統(tǒng)一采用有限差分方法獲得,以(x,y)=(400 m,400 m)處為第一個炮點和檢波點的位置,在地表上均勻布置81個炮點,每炮由6561個檢波點接收,炮點空間間距為400 m×400 m,檢波點間距為40 m×40 m.采用爆炸震源激發(fā),震源子波為主頻25 Hz的雷克子波,采樣間隔為4.8 ms,采樣時間為3.6 s.

第一次有限差分正演采用速度模型(見圖10a)和常密度(1000 kg·m-3)模型進行正演得到無散射點的常規(guī)地震數(shù)據(jù),第二次有限差分正演仍采用相同的速度模型和密度模型,只是在原來常密度模型基礎(chǔ)上加上一些散射點(見圖10b,散射點間距為500 m×500 m×500 m),采用相同的觀測系統(tǒng)進行正演.將兩次正演的波場記錄相減,作為散射點波場,之后對散射點波場進行90°相位校正,并對校正之后的數(shù)據(jù)進行高斯束偏移成像,便可以得到密度擾動點處散射點的PSF.如圖所示,圖11a為Salt模型的三維點擴散函數(shù),圖11b為從圖11a中截取的一個二維PSF剖面.由圖11可以看出: PSF的形狀和角度照明有關(guān),在圖像空間中接近圓形(所有角度)的PSF表示大范圍的傾角照明,扁的PSF表示小范圍的傾角照明; PSF的振幅反映了照明的強弱和偏移孔徑的大小; 此外PSF的大小是隨空間變化的,可以看做隨著深度變化的三維子波.以上分析表明了PSF可以很好地捕獲照明、帶寬和采集系統(tǒng)幾何形態(tài)的影響,并且有效避免了對變化子波進行估計的問題.

圖10 三維Salt模型的一部分

圖11 Salt模型的點擴散函數(shù)

有了擾動點處的PSFs,便可以通過空間線性插值的方法得到每個網(wǎng)格點的PSF,之后將反射率模型與插值后的PSF進行空間褶積處理就可以得到模擬成像結(jié)果.如圖所示,圖12a為實際的高斯束偏移結(jié)果,其淺部噪聲是因為受到了直達波的影響,圖12b為通過PSF得到的深度域模擬成像(僅對數(shù)值的量級做了歸一化處理),可以看到模擬成像剖面和實際的高斯束深度偏移剖面在相位和振幅上可以很好地進行匹配.

圖12 Salt模型模擬成像和偏移成像的比較

為了檢驗本文所提出方法的適應(yīng)性,又對更為復(fù)雜的Seam模型進行了測試.如圖13所示為二維Seam速度模型,模型大小為3501×1501個網(wǎng)格點,網(wǎng)格間距設(shè)為10 m以第一個網(wǎng)格點為初始炮點和檢波點位置,均勻放置了175炮,炮間距為200 m,共3501個檢波點,檢波點間距為10 m.震源為爆炸震源,震源子波為主頻為25 Hz的雷克子波,采樣間隔為4 ms,采樣時間為12 s.采用同樣的方式我們可以得到Seam模型的PSF和深度域模擬成像.如圖14所示為當密度擾動強度過大時(擾動點間隔為600 m×600 m,密度是背景密度的10倍)得到的PSF,可以看出此時PSF受到了很強的結(jié)構(gòu)信息污染,嚴重影響了后續(xù)的插值和空間褶積的準確性.通過調(diào)整擾動強度的大小,我們發(fā)現(xiàn)當擾動量減小到δρ=0.1ρ0時可以很好的壓制模型的結(jié)構(gòu)信息,得到更加清晰的PSF(見圖15).

圖13 Seam速度模型

圖14 被結(jié)構(gòu)信息污染的PSFs

圖15 Seam 模型的PSFs

為了更加直觀地對比高斯束深度偏移結(jié)果(見圖16)和由PSF得到的模擬成像結(jié)果(見圖17)的匹配程度,將兩幅圖進行拼接對比(見圖18): 其中兩個箭頭中間的部分為深度域模擬成像結(jié)果,箭頭兩邊的部分為實際高斯束偏移結(jié)果,可以看到在拼接處二者幾乎完全重合在一起,進一步驗證了利用PSF直接在深度域進行模擬成像的可行性.在鹽丘邊界,速度模型存在高對比度間斷,使得偏移成像本身存在低頻噪聲.同時,速度不連續(xù)引起子波拉伸和照明變化使得PSF在該區(qū)域變化劇烈,簡單的線性插值難以準確的表達該位置的PSF.盡管可以通過設(shè)置掩膜算子將鹽丘兩邊分離開,避免跨邊界插值來提高模擬成像的準確度,但是仍會產(chǎn)生一定的誤差使得模擬成像和偏移成像在鹽丘邊界位置匹配欠佳.

圖16 Seam模型的高斯束偏移結(jié)果

圖17 Seam模型的PSF模擬成像結(jié)果

圖18 Seam模型模擬成像結(jié)果與實際偏移結(jié)果的對比

盡管高斯束偏移結(jié)果和模擬成像結(jié)果存在一定的振幅差異,這是由于偏移算法不是保幅的,并且在模擬成像時采用垂直反射率(公式(35))與PSF空間褶積無法完全反映實際的偏移過程.針對鹽丘底部照明不足的問題,可以將PSF作為反褶積算子對偏移結(jié)果進行照明補償.此外,可以通過疊前角度域PSF對真實反射系數(shù)進行迭代求解,提高反演結(jié)果的準確性.

綜上所述,通過PSF與反射率模型進行空間褶積可以直接在深度域得到模擬成像,所得到的模擬成像結(jié)果可以用于建立深度成像域反演的目標函數(shù).與數(shù)據(jù)域最小二乘反演相比,基于PSF的成像域反演方法通過PSF模擬成像和實際偏移成像的殘差來建立目標泛函,不需要進行反復(fù)的偏移和反偏移迭代,有助于大幅提高反演的計算效率.同時,因為PSF的計算和深度偏移擁有相同的采集系統(tǒng)幾何形狀、偏移算子和參數(shù)模型,PSF可捕獲成像系統(tǒng)的弱項或不足,因此可以通過PSF對由成像系統(tǒng)造成的模糊進行去除,并且降低了對偏移算法保幅性的要求(理論上可以采用任何一種深度偏移方法來進行處理).

3 結(jié)論

本文從擾動理論出發(fā),通過分析速度擾動和密度擾動對散射場的影響,在前人工作的基礎(chǔ)上提出了基于密度擾動的兩步正演法來計算點擴散函數(shù),避免了速度擾動對運動學(xué)信息的破壞.進一步研究了利用點擴散函數(shù)進行深度域模擬成像時存在的問題: 提出可以通過調(diào)整擾動強度的方式減少結(jié)構(gòu)信息對PSF的污染,提高了復(fù)雜模型中PSF計算的準確性; 指出散射場存在的相位畸變,通過進行90°相位校正使得模擬成像與偏移成像的相位得到很好的匹配; 實現(xiàn)了PSF的動態(tài)插值,明確了空間褶積過程中不需要對褶積核(PSF)進行翻轉(zhuǎn).模型實驗表明: PSF可以很好的捕獲照明、帶寬和采集系統(tǒng)幾何形態(tài)的影響,通過PSF得到的模擬成像在振幅和相位上與實際偏移成像匹配良好.這些都進一步證明了利用PSF直接在深度域進行反演的可行性.當然,對于PSF深度域反演方法仍然有許多問題尚待解決.本文只是針對PSF的計算給出了一個相當可行的方法,并對利用PSF進行深度域模擬成像過程中存在的問題進行了深入的探討.當前的研究主要是基于聲波方程,關(guān)于彈性/各項異性PSF計算和復(fù)雜模型的PSF插值算法仍需要進行更加深入的研究.