基于密度擾動(dòng)的點(diǎn)擴(kuò)散函數(shù)計(jì)算與三維深度域模擬成像

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

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

0 引言

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

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

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

1 原理

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

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

(1)

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

(2)

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

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

fpsf=ATdsc=ATAmsc≈H,

(3)

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

(4)

1.2 速度和密度擾動(dòng)對散射場的影響

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

(5)

其中s表示震源位置,x是散射點(diǎn)位置,t表示時(shí)間.又因?yàn)棣?κ和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)的擾動(dòng)量,p0是背景波場,p′為散射波場.將公式(6)帶入公式(5)并應(yīng)用

(7)

之后將無擾動(dòng)項(xiàng)和擾動(dòng)項(xiàng)波場分別放在等號(hào)兩邊，可得到

(8)

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

(10)

(11)

(12)

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

(14)

(15)

(16)

(17)

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

圖1 震源-散射點(diǎn)-檢波器之間的散射張角示意圖

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

(18)

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

(19)

其中

(20)

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

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

(21)

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

(22)

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

(23)

1.3 點(diǎn)擴(kuò)散函數(shù)的計(jì)算

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

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

(24)

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

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

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

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

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

圖2 速度/密度擾動(dòng)示意圖

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

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

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

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

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

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

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

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

1.4 深度域模擬成像

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

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

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

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

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

(35)

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

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

2 模型測試

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

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

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

圖10 三維Salt模型的一部分

圖11 Salt模型的點(diǎn)擴(kuò)散函數(shù)

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

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

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

圖13 Seam速度模型

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

圖15 Seam 模型的PSFs

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

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

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

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

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

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

3 結(jié)論

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