面向模型空間地震反演成像的點(diǎn)擴(kuò)散函數(shù)快速計(jì)算方法

李瀚野,任浩然,陶柳蓉,江金生

(浙江大學(xué)地球科學(xué)學(xué)院,浙江杭州310027)

傳統(tǒng)偏移成像方法可以較好地處理地震波的運(yùn)動(dòng)學(xué)信息,但大多采用線性正演算子的共軛轉(zhuǎn)置,而非逆算子作為其偏移算子[1-2]。因此得到的結(jié)果不能準(zhǔn)確地反映真實(shí)的反射系數(shù),影響最終的成像質(zhì)量。人們對(duì)由此引起的振幅不均、低頻噪聲、低信噪比等各種問題已進(jìn)行過大量研究[3-5]。最小二乘偏移(least-squares migration,LSM)[6-8]通?；谔荻阮惖乃惴?如最速下降法等),通過不斷匹配觀測(cè)數(shù)據(jù)和模擬數(shù)據(jù)之間的殘差,使殘差最小化。該方法可以隱式消除Hessian效應(yīng)[9],有效壓制噪聲,提升成像剖面的分辨率,進(jìn)而獲得真振幅剖面。最小二乘偏移的計(jì)算成本高,以射線類的方法為例,雖然基于射線的偏移算法具有實(shí)現(xiàn)效率高的優(yōu)勢(shì),但是在最小二乘偏移的框架中,每次迭代均需進(jìn)行一次偏移和反偏移,通常需要10次或更多的迭代次數(shù),這將產(chǎn)生至少20倍于傳統(tǒng)偏移方法的計(jì)算成本[10]。與此同時(shí),不夠準(zhǔn)確的背景速度場(chǎng)、含噪數(shù)據(jù)、未知的震源子波等許多因素將導(dǎo)致數(shù)據(jù)空間最小二乘偏移的收斂速度緩慢,甚至不收斂,實(shí)際情況下,計(jì)算成本將會(huì)更高。因此最小二乘偏移方法的實(shí)現(xiàn)效率成為了一個(gè)重要的研究方向。AOKI等[10]構(gòu)建了一種去模糊化濾波器,可以應(yīng)用于最小二乘偏移中的預(yù)條件策略或正則化策略中,相比傳統(tǒng)算法,該濾波器在保證高精度成像的情況下,可以將計(jì)算效率提高約30%。KUEHL等[11]基于多平方根方程提出了最小二乘裂步偏移方法,可以減小觀測(cè)數(shù)據(jù)值和合成數(shù)據(jù)值之間的誤差。YAO等[12]提出了一種非線性最小二乘逆時(shí)偏移算法,該算法直接從反射率中生成預(yù)測(cè)數(shù)據(jù),無需假設(shè)它們之間存在線性關(guān)系。ZHANG等[13]提出了一種新的最小二乘逆時(shí)偏移方法,通過觀測(cè)數(shù)據(jù)和模擬數(shù)據(jù)之間的互相關(guān)最大化,來計(jì)算最小二乘逆時(shí)偏移的最快收斂速度值。YU等[14]利用偏移解卷積濾波器處理局部偏移圖像,僅對(duì)關(guān)注的成像區(qū)域進(jìn)行處理,避免了全局計(jì)算,該方法相當(dāng)靈活。LECOMTE[15]使用基于射線的方法來估計(jì)Hessian的逆,計(jì)算過程簡(jiǎn)潔。TANG[16]改進(jìn)了原有的顯式Hessian公式,使Hessian算子的計(jì)算更加高效,利用隨機(jī)相位編碼的Hessian算子可以得到振幅更為均衡的鹽下成像。

反演理論證明,常規(guī)偏移成像的結(jié)果是真實(shí)反射系數(shù)經(jīng)模糊化算子卷積得到的退化結(jié)果[17],該模糊化算子即Hessian算子,許多學(xué)者進(jìn)行過關(guān)于Hessian算子的討論,XIE等[18]提出了照明矩陣、BOSCHI[19]提出了模型分辨率矩陣、YU等[14]提出了偏移解卷積濾波器、SCHUSTER等[20]提出了偏移格林函數(shù)、WANG等[21]提出了解卷積算子、FICHTNER等[22]提出了核函數(shù),以上概念實(shí)質(zhì)上都是Hessian算子的不同稱謂。在物理學(xué)意義上,Hessian算子包含了觀測(cè)系統(tǒng)和速度模型的所有信息[23],研究Hessian算子對(duì)于發(fā)展地震正、反演工作有著重要的啟示作用。應(yīng)用模型空間的最小二乘偏移方法是通過求取Hessian算子,顯式地壓制Hessian效應(yīng)。這種策略避免了在每一輪迭代中進(jìn)行偏移和反偏移計(jì)算,具有更高的實(shí)現(xiàn)效率,但Hessian算子在線性化條件下是一個(gè)巨型的矩陣,故實(shí)際運(yùn)算中難以計(jì)算和保存。實(shí)際應(yīng)用中需要引入Hessian近似,在高頻漸進(jìn)條件下,Hessian矩陣可以看作對(duì)角矩陣[24-25],WANG等[26]利用對(duì)角Hessian的逆來重建地震反演中的深部結(jié)構(gòu),可以平衡偏移圖像振幅以及消除幾何傳播效應(yīng)。點(diǎn)擴(kuò)散函數(shù)(point spreading function,PSF)代表Hessian矩陣的一行元素,通常將PSF元素投影在模型空間,在二維情況下表現(xiàn)為目標(biāo)點(diǎn)周圍的一個(gè)橢圓型花樣,代表地下空間的一個(gè)單點(diǎn)對(duì)于觀測(cè)系統(tǒng)的脈沖響應(yīng)。通過求取地下空間的PSF,可以得到Hessian算子的近似結(jié)果,王雪君等[27]直接求取基于PSF的逆處理黏彈性介質(zhì)逆時(shí)偏移成像成果,有效均衡了成像振幅,提高了分辨率。XU等[28]利用基于射線的方法生成PSF,該方法在提高分辨率的同時(shí),兼具較高的靈活性。姚振岸等[29]構(gòu)建了基于PSF的去模糊濾波器,將其作為最小二乘逆時(shí)偏移迭代過程中的預(yù)條件化因子,得到了高分辨率的偏移結(jié)果。

為了提升計(jì)算效率,本文在兩個(gè)方面進(jìn)行了方法改進(jìn):①建立散射點(diǎn)模型,通過一次正演+偏移快速計(jì)算所有控制點(diǎn)點(diǎn)擴(kuò)散函數(shù),并利用插值迅速得到全成像空間的點(diǎn)擴(kuò)散函數(shù),進(jìn)而得到近似Hessian算子;②在模型空間展開反演。對(duì)于數(shù)據(jù)空間的反演而言,每一次迭代都需要進(jìn)行偏移和反偏移計(jì)算,且輸入龐大的地震數(shù)據(jù)集。模型空間的反演可轉(zhuǎn)化為圖像的去模糊化問題,因其輸入為小尺寸的偏移圖像,因此計(jì)算效率更高。本文采用模型空間反演的方法,只需比數(shù)據(jù)空間小得多的計(jì)算量就可得到常規(guī)偏移成像去模糊化的反演結(jié)果。

1 方法原理

1.1 模型空間的最小二乘偏移成像理論

在反演理論框架下,利用最小二乘偏移方法將數(shù)據(jù)空間和模型空間之間的投影定義為如下線性問題:

d=Lm

(1)

式中:d為Born近似下僅包含一次反射和散射的模擬數(shù)據(jù);L為線性正演算子;m為地下模型的高頻擾動(dòng)項(xiàng)。傳統(tǒng)的反演方法通過提出目標(biāo)函數(shù)J(m)來實(shí)現(xiàn)模擬數(shù)據(jù)和觀測(cè)數(shù)據(jù)之間殘差的最小化。

J(m)=‖Lm-d2‖

(2)

最小二乘偏移時(shí),m的計(jì)算式如下:

m=(L*L)-1L*d=H-1L*d

(3)

式中:L*為伴隨算子;H=L*L為Hessian算子。Hessian算子在數(shù)學(xué)上表示為目標(biāo)函數(shù)對(duì)模型的二階導(dǎo)數(shù),在物理上它代表觀測(cè)系統(tǒng)對(duì)模型的分辨能力。mmig=L*d表示偏移成像的過程。(3)式可以改寫為L(zhǎng)*d=Hm,觀察等式兩邊可知,偏移成像是特定觀測(cè)系統(tǒng)下成像過程中對(duì)真實(shí)模型模糊化的結(jié)果,Hessian算子即模糊化算子。

在數(shù)據(jù)空間中進(jìn)行最小二乘偏移需要在每次成像過程中,將成像結(jié)果作為二次源生成模擬數(shù)據(jù),每次迭代均需要進(jìn)行一次正演和偏移成像計(jì)算,產(chǎn)生了巨大的計(jì)算量。模型空間的最小二乘偏移算法可以克服計(jì)算成本上的缺陷,目標(biāo)函數(shù)如下:

J(m)=‖Hm-mmig‖2

(4)

無論是在數(shù)據(jù)空間還是模型空間中開展反演,最小二乘偏移的核心都在于求解Hessian算子及其逆算子?；诓▌?dòng)方程得到的Hessian算子,其逆可作為去卷積算子。使用梯度類方法解決該問題時(shí),該算子可作為步長(zhǎng)。Hessian矩陣中對(duì)角元素的重要性已經(jīng)得到了充分的證明[24-26],忽略干擾項(xiàng)后,Hessian算子可寫成如下格林函數(shù)表達(dá)式:

(5)

式中:G為格林函數(shù);ks為震源;kr為檢波點(diǎn);fs(ω)為能量在圓頻率ω上的譜;k和l為不同成像點(diǎn),(1)式中的正演算子L可表示為:

L(ks,kr,ω)|x=ω2|fs(ω)|G(ks,k,ω)·
G(k,kr,ω)

(6)

1.2 點(diǎn)擴(kuò)散函數(shù)與Hessian算子

Hessian算子是一個(gè)巨型矩陣,由(5)式可知,H(k,l)對(duì)應(yīng)兩組成像點(diǎn)k和l。如圖1所示,Hessian矩陣的每一行和每一列均對(duì)應(yīng)于成像空間中的一個(gè)點(diǎn)。抽取線性Hessian矩陣中的一行數(shù)據(jù)即點(diǎn)擴(kuò)散函數(shù)(PSF),如下:

圖1 Hessian作用示意

(7)

通常在二維情況下,PSF在相應(yīng)點(diǎn)周圍展開為一個(gè)橢圓形的斑點(diǎn)。如圖2所示,在反演框架下,疊前偏移的結(jié)果為真實(shí)反射率被PSF模糊化后的結(jié)果,因此PSF被視作構(gòu)建Hessian算子近似的有力工具。PSF的精確求解需要大量的計(jì)算,這里將(7)式分解為如下兩個(gè)部分:

圖2 PSF模糊和去模糊化過程

P1(ω,ks,kr,k)=ω2fs(ω)G(ks,k,ω)G(k,kr,ω)

(8)

P2(ω,ks,kr,l)=ω2fs(ω)G(ks,l,ω)G(l,kr,ω)

(9)

對(duì)于一組震源ks和檢波點(diǎn)kr,(8)式和(9)式中的P1和P2僅與k點(diǎn)和其周圍的l點(diǎn)有關(guān)。如果給定一個(gè)參考點(diǎn)k,那么只需要計(jì)算P1和其周圍的P2,并通過互相關(guān)得到PSF的圖像。

為了得到(8)式和(9)式中格林函數(shù)的解,我們?cè)谡鹪春蜋z波點(diǎn)分別激發(fā)沿地下介質(zhì)傳播的波場(chǎng),并記錄在時(shí)域或頻域中,然后通過頻域正演得到精確的解。當(dāng)模型空間較大時(shí),計(jì)算和存儲(chǔ)的成本巨大。

由(6)式至(9)式可知,L=L(ks,kr,ω),進(jìn)而可得:

P1=P1(ω,ks,kr,k)=L|k

(10)

P2=P2(ω,ks,kr,l)=L|l

(11)

P(l)|k=L|k·L*|l

(12)

(13)

PSF的計(jì)算過程為這組數(shù)據(jù)的成像過程:

(14)

我們將PSF的計(jì)算過程轉(zhuǎn)換為一次正演+偏移的過程,因?yàn)樵撨^程需要遍歷所有的成像空間,故仍存在巨大的計(jì)算量。從圖1可見,PSF在散射點(diǎn)周圍擴(kuò)散成一個(gè)斑點(diǎn)。因此,我們首先可以選擇一定間隔的離散點(diǎn)作為控制點(diǎn),并保證每個(gè)控制點(diǎn)生成的PSF斑點(diǎn)不會(huì)互相重合,然后使用上述過程一次性計(jì)算所有控制點(diǎn)上的PSF,最后在成像空間上對(duì)與控制點(diǎn)相鄰的點(diǎn)進(jìn)行空間插值,快速獲得全成像空間的PSF。插值方法采用反距離加權(quán)(inverse distance weight,IDW)插值:

(15)

式中:Wi表示權(quán)值,反映插值點(diǎn)受周圍控制點(diǎn)影響的程度;hi表示插值點(diǎn)到控制點(diǎn)的距離;p為任意一個(gè)正實(shí)數(shù),p值越大,則更強(qiáng)調(diào)鄰近控制點(diǎn)的作用,反之相反。

1.3 基于PSF的圖像去模糊方法

解決圖像的去模糊問題,關(guān)鍵在于模糊(卷積)算子的估計(jì)和卷積效應(yīng)的消除,利用1.2節(jié)的方法可以高效求取模糊算子。那么通過對(duì)常規(guī)偏移成像結(jié)果進(jìn)行解卷積處理即可得到去模糊圖像,一個(gè)二維圖像的模糊正過程可由下式描述:

g(x,y)=f(x,y)?k(x,y)

(16)

式中:?表示卷積運(yùn)算;f(x,y)表示真實(shí)圖像;k(x,y)表示模糊算子;g(x,y)表示模糊后的退化圖像。(16)式在波數(shù)域有如下表達(dá)式:

G(u,v)=F(u,v)·K(u,v)

(17)

空間域的卷積對(duì)應(yīng)于波數(shù)域的乘積,其中G(u,v)、F(u,v)、K(u,v)分別由g(x,y)、f(x,y)、k(x,y)經(jīng)二維傅里葉變換得到,u和v分別代表橫、縱方向上的頻率,(17)式可改寫為如下形式:

(18)

(19)

基于點(diǎn)擴(kuò)散函數(shù)快速計(jì)算的模型空間反演成像流程如圖3所示。

圖3 基于點(diǎn)擴(kuò)散函數(shù)快速計(jì)算的模型空間反演成像流程

2 數(shù)值實(shí)驗(yàn)

圖4為SEG/EAGE鹽丘速度模型,模型網(wǎng)格點(diǎn)數(shù)為1200×150,兩個(gè)方向的采樣間隔均為10m,速度為1500～4500m/s。合成地震數(shù)據(jù)的過程中炮間隔為100m,共計(jì)80炮,從地表2000m處開始均勻放炮。采集系統(tǒng)位于地表,排列長(zhǎng)度為整個(gè)模型的長(zhǎng)度,道間距為10m,震源為主頻為12.5Hz的雷克子波。計(jì)算時(shí)采用時(shí)間精度二階,空間精度十階的時(shí)間域有限差分正演及完全匹配吸收邊界條件。

我們將Hessian近似的計(jì)算,即PSF的計(jì)算轉(zhuǎn)換為一次“正演+偏移+插值”過程,目的是保證獲得清晰的PSF斑點(diǎn)并且PSF斑點(diǎn)之間不互相重疊。對(duì)初始速度模型進(jìn)行平滑處理后,將控制點(diǎn)均勻分布在平滑后的速度模型上,得到散射點(diǎn)模型(圖5),控制點(diǎn)處速度為原速度的1.5倍。

利用散射點(diǎn)模型進(jìn)行正演,并對(duì)得到的地震數(shù)據(jù)進(jìn)行逆時(shí)偏移(RTM)成像,所得結(jié)果為控制點(diǎn)上的脈沖響應(yīng),即PSF,圖6展示了圖5散射點(diǎn)模型生成的PSF圖像。對(duì)于控制點(diǎn)之外的點(diǎn),我們使用IDW插值生成其PSF,圖7為冪值參數(shù)p在不同取值下的全成像點(diǎn)PSF圖像(橫向網(wǎng)格點(diǎn)數(shù):1200×41、縱向網(wǎng)格點(diǎn)數(shù):150×21)。當(dāng)p值為1時(shí),所得的插值圖像存在明顯的網(wǎng)格狀間斷,隨著取值的增大,圖像更加平滑。在后續(xù)測(cè)試中我們將繼續(xù)討論不同p值對(duì)結(jié)果產(chǎn)生的影響。

為了驗(yàn)證計(jì)算得到的PSF是否可以起到近似Hessian算子造成的模糊化效果,首先從圖4的速度模型中計(jì)算反射率,得到如圖8所示的鹽丘反射率模型。我們將反射率模型與計(jì)算得到的PSF進(jìn)行卷積,得到如圖9a所示的圖像(圖像中空白為卷積運(yùn)算無法覆蓋區(qū)域,顯示時(shí)置零),此時(shí)p取5,可以看到圖9a的分辨率低于圖7,PSF的模糊化效應(yīng)平滑了邊界,起到了模糊化的作用。觀察圖9b至圖9e不同p值下的卷積局部放大結(jié)果,可發(fā)現(xiàn)當(dāng)插值冪次選取不當(dāng),會(huì)在模糊化(卷積)的正過程中引入網(wǎng)格狀間斷,而去模糊化問題的實(shí)現(xiàn)有賴于PSF的精度,即正過程效果越好,則去模糊化(解卷積)的效果越好,這些間斷會(huì)在后續(xù)反演過程中引入不必要的噪聲,圖9d 的正過程效果最為理想,故本實(shí)驗(yàn)選取p=5作為插值參數(shù)。

圖4 SEG/EAGE鹽丘速度模型

圖5 圖4速度模型對(duì)應(yīng)的散射點(diǎn)模型

圖6 圖5散射點(diǎn)模型生成的PSF圖像

圖7 不同p值下的全成像點(diǎn)PSF圖像

圖8 鹽丘反射率模型

在對(duì)RTM成像結(jié)果去模糊化前,可以先利用卷積模型驗(yàn)證波數(shù)域?yàn)V波方法的有效性,將圖9a作為圖像去模糊化問題的輸入,采用波數(shù)域?yàn)V波算法進(jìn)行解卷積運(yùn)算,得到如圖10所示的去模糊化結(jié)果,該模型與圖8中的原始反射率模型相似,證明了當(dāng)PSF模糊效應(yīng)的正反過程匹配時(shí),波數(shù)域?yàn)V波算法有效。

圖9 鹽丘反射率模型卷積結(jié)果及不同p值下的卷積結(jié)果局部放大顯示

圖10 鹽丘反射率模型的去模糊化結(jié)果

采用與計(jì)算PSF相同的觀測(cè)系統(tǒng)和震源子波,根據(jù)圖4的速度模型生成地震數(shù)據(jù),并進(jìn)行RTM成像,得到如圖11所示的結(jié)果。將波數(shù)域?yàn)V波方法應(yīng)用于圖11的鹽丘模型RTM成像結(jié)果,得到的該模型RTM去模糊化結(jié)果如圖12所示,可以看出,處理后的結(jié)果中平層及鹽丘頂部分辨率明顯提升,層位邊界得到銳化。淺地表震源效應(yīng)引起的線性噪聲在濾波后得到了較好的壓制,同時(shí)深部層位的振幅能量一定程度上被抬高,一些層位的連續(xù)性得到了提升。觀察淺藍(lán)線框可以發(fā)現(xiàn),濾波補(bǔ)齊了弱照明區(qū)域的層位(受觀測(cè)孔徑的影響,補(bǔ)出的層位向上傾斜),但會(huì)在鹽丘內(nèi)部引入噪聲。圖13為鹽丘模型波數(shù)域反演前、后的波數(shù)譜,不難發(fā)現(xiàn),反演后的結(jié)果在波數(shù)域有著更寬的展布。

圖11 鹽丘模型RTM成像結(jié)果

圖12 鹽丘模型RTM成像去模糊化結(jié)果

采用圖14的Marmousi速度模型進(jìn)一步驗(yàn)證本文方法的有效性,Marmousi速度模型的RTM成像結(jié)果如圖15所示。采用本文方法計(jì)算得到的PSF圖像如圖16所示,應(yīng)用波數(shù)域?yàn)V波算法處理圖15的RTM剖面,得到如圖17所示的去模糊化結(jié)果,圖18為Marmousi模型波數(shù)域反演前、后的波數(shù)譜。對(duì)比圖17和圖15可以發(fā)現(xiàn),本文方法可以銳化邊界,提高層位分辨率,同時(shí)提高深部振幅能量,相較于RTM成像結(jié)果,采用本文方法得到的反演結(jié)果具有更加均衡的振幅。

圖14 Marmousi速度模型

圖15 Marmousi速度模型RTM成像結(jié)果

圖16 Marmousi散射點(diǎn)模型生成的PSF圖像

圖17 Marmousi模型RTM成像去模糊化結(jié)果

圖18 Marmousi模型波數(shù)域反演前(a)、后(b)的波數(shù)譜

3 討論

1) 正問題是反問題的基礎(chǔ),圖像去模糊化工作有賴于PSF的精度,本文通過散射點(diǎn)模型一次“正演+偏移+插值”的過程獲取PSF。在速度變化劇烈的區(qū)域,PSF的結(jié)果會(huì)有較大變動(dòng)。因此,控制點(diǎn)布設(shè)越密,插值結(jié)果的準(zhǔn)確性越高。控制點(diǎn)的疏密不會(huì)影響計(jì)算量,但是布設(shè)過密時(shí)會(huì)導(dǎo)致PSF樣式互相重疊,影響精度。故在求取PSF時(shí)通常需要1～2次試驗(yàn)(正演+偏移)來獲得最佳的散射點(diǎn)布設(shè)間隔。

2) 對(duì)于插值方法中的冪值參數(shù)選取,以正過程(反射率模型同PSF卷積)結(jié)果不會(huì)引入額外噪聲為原則。

3) 模型空間的反演基于偏移成像結(jié)果展開,初始偏移成像結(jié)果的精度制約了最終的反演效果,這個(gè)過程與數(shù)據(jù)空間一樣需要相對(duì)精確的偏移速度。

4) 本文通過RTM的結(jié)果證明了方法的有效性,該方法的核心在于PSF的求取和圖像去模糊,因此該策略同樣適用于其它偏移算子方法進(jìn)行求解。

4 結(jié)論

本文分析了PSF和Hessian算子的數(shù)學(xué)物理性質(zhì),提出了一種對(duì)散射點(diǎn)模型進(jìn)行一次“正演+偏移+插值”獲取點(diǎn)擴(kuò)散函數(shù)的快速算法?？刂泣c(diǎn)的PSF僅需要一次正演和偏移,而后可以通過IDW插值快速獲得全成像空間的PSF。該方法具有較高的計(jì)算效率和靈活性。與反射率模型的卷積結(jié)果對(duì)比發(fā)現(xiàn),該方法生成的PSF可以很好的起到模糊化作用。在得到模糊算子后,將反演目標(biāo)轉(zhuǎn)化為模型空間下的圖像去模糊化問題,通過波數(shù)域反演算法一次運(yùn)算就可以得到反演結(jié)果。反演過程僅需進(jìn)行兩組二維傅里葉正反變換和相應(yīng)的矩陣運(yùn)算,具有較高的實(shí)現(xiàn)效率。模型測(cè)試結(jié)果表明,本文方法可以提高常規(guī)偏移成像的分辨率,均衡弱照明區(qū)域振幅。

致謝:感謝中國(guó)石油化工股份有限公司石油物探技術(shù)研究院的資助。感謝同濟(jì)大學(xué)波現(xiàn)象與智能反演成像研究組(WPI)徐鵬博士的有益討論。