王雪君,任浩然,江金生,李瀚野

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

吸收衰減是地球介質(zhì)的一種基本屬性。面對復雜山前帶地震成像中小尺度油氣藏、深層油氣藏以及非常規(guī)油氣藏勘探,高分辨率地震成像越來越成為業(yè)界關(guān)注的焦點。對于高分辨率成像,地層的吸收衰減效應是一個重要的影響因素。有效的地震波吸收與衰減補償是地震資料處理工作的一個重要環(huán)節(jié)。

地層的黏滯性會導致地震波傳播過程中振幅衰減和相位變化,傳統(tǒng)的聲波逆時偏移和最小二乘逆時偏移不能對其進行校正,從而導致地下反射體能量不能得到有效的聚焦和歸位。品質(zhì)因子Q描述了介質(zhì)的吸收衰減特性,在地震勘探頻譜范圍內(nèi)通常使用常Q衰減模型[1],在此模型基礎(chǔ)上發(fā)展了很多黏滯聲波方程。在過去的20年中,許多研究使用了頻率域單程波方程來補償上、下行波的衰減和頻散效應[2-7]。這些方法通常將衰減項放在復相速度虛部,以便實現(xiàn)反向補償。與頻率域黏滯聲波方程相比,時間域黏滯聲波方程具有更高的計算效率。ZHU等[8-9]提出了一種吸收衰減介質(zhì)的時間域黏彈性波動方程,描述了常Q吸收衰減介質(zhì)中的地震波傳播,且衰減和頻散運算符在這個方程中能夠被分離,從而通過反轉(zhuǎn)振幅算符符號而頻散算符符號不變來補償振幅衰減和相位頻散[10]。ZHU等[11]提出了基于時間域黏彈性波動方程及其反向傳播(時間反向)建模的黏彈性逆時偏移,在黏彈性波動方程中解耦P波和S波的衰減特性,從而在波場外推期間補償P波和S波的衰減效應,且設(shè)計了一種黏彈性逆時傳播方法,通過反轉(zhuǎn)P波和S波振幅損失算子的符號,來校正P波和S波的衰減效應?；诟飨蛲责椥怨?ZHU[12]將各向同性公式擴展到一般的各向異性黏彈性波動方程來模擬各向異性衰減介質(zhì),且推導出了時間域位移-應力公式。

基于牛頓法的全波形反演方法和基于線性化的高斯-牛頓法的最小二乘偏移方法都需要求解Hessian核函數(shù)。Hessian核函數(shù)在數(shù)學上是反問題泛函對模型參數(shù)的二階導數(shù)。而在物理學上,Hessian核函數(shù)疊合了地震波場正傳播和算子反傳播的所有信息[13]。Hessian核函數(shù)的敏感性更加能夠體現(xiàn)對地下成像空間每個點的照明響應,研究Hessian核函數(shù)對理解地震波場正、反演過程有著重要意義。在地震反演的研究中,Hessian算子有多種稱謂,如SCHUSTER等[14]提出的偏移格林函數(shù)、BOSCHI[15]提出的模型分辨率矩陣、FICHTNER等[16]提出的點傳播矩陣、XIE等[17]提出的照明矩陣、YU等[18]提出的去模糊化算子、WANG等[19]提出的偏移反褶積算子。這些概念都是Hessian算子的各種變稱。在局部線性化情況下,Hessian核函數(shù)又可以寫為線性的Hessian矩陣,當一個Hessian矩陣作用于最小二乘成像,會使得成像結(jié)果模糊,因此常規(guī)疊前偏移成像實際上是模糊的、振幅畸變的成像結(jié)果。點擴散函數(shù)(Point Spreading Function,PSF)是Hessian矩陣的一行,對應于地下一個成像點,反映了觀測系統(tǒng)對該點的照明能力。通常,PSF體現(xiàn)為目標點周圍的一個橢圓。反演理論證明,常規(guī)偏移成像結(jié)果就是該函數(shù)對理論成像結(jié)果的“模糊化”[20]。CHEN等[21]利用一個去模糊化算子來實現(xiàn)黏聲介質(zhì)的最小二乘逆時偏移。研究帶Q因子的Hessian核函數(shù),找出其在不同模型下的特征規(guī)律,解決Hessian矩陣的求解與求逆問題,是發(fā)展更高精度更高分辨率地震反演成像方法的理論基礎(chǔ),可以為地球內(nèi)部結(jié)構(gòu)研究與油氣勘探提供更有利的成像工具。Hessian核函數(shù)的求解及求逆策略可以通過Hessian點響應矩陣的求逆來開展。因此,研究PSF,發(fā)展帶Q因子的PSF的求逆策略,可以將其逆算子直接作用在成像結(jié)果上從而提高成像分辨率。

本文是在ZHU等[8-10]和羅文山等[22]的研究基礎(chǔ)上,利用常Q模型吸收衰減與補償統(tǒng)一表達的一階“速度-壓力”黏聲波方程實現(xiàn)黏聲介質(zhì)地震波場的格林函數(shù)計算,對Hessian矩陣分解與求逆方法進行研究,基于構(gòu)建去模糊化算子,開展點擴散函數(shù)的逆矩陣研究,將逆Hessian點響應作用在成像結(jié)果上,從而對原始成像結(jié)果有效地去模糊化,提高了成像振幅均衡性和分辨率。最后用模型數(shù)據(jù)對方法的有效性進行了測試。

1 方法原理

1.1 黏聲介質(zhì)地震波格林函數(shù)

有效的地震波吸收與衰減補償首先要解決的問題是如何有效描述地震波吸收與衰減的過程。地震波場的互易性理論揭示,地震波場傳播的正、反過程可以統(tǒng)一利用格林函數(shù)來表達,因此一個合理的、帶吸收衰減效應的格林函數(shù)表達式是進行反演理論框架下的地震成像的基礎(chǔ)。

1.1.1 黏聲介質(zhì)地震波波動方程

吸收衰減介質(zhì)中的應力-應變關(guān)系為:

(1)

式中:σ為應力;ε為應變;p為壓力;ψ1-2γ(t)為弛豫時間函數(shù)。且:

(2)

式中:M0和t0分別為弛豫體積模量和弛豫時間常數(shù);Γ是歐拉γ函數(shù);γ與Q有關(guān);ρ,c0,Q分別為密度、參考角頻率ω0對應的速度和地震品質(zhì)因子。

基于ZHANG等[5]推導的時間域黏聲波動方程,可以寫出吸收衰減介質(zhì)地震波場衰減與補償?shù)慕y(tǒng)一公式:

(3)

其中,

當α=1時,公式(3)為衰減公式;α=-1時,公式(3)為補償公式[10]。(3)式的一階速度-壓力形式為:

(4)

其中,u=?p/?t,vx=?p/?x,vz=?p/?z。當Q趨向無窮大時,γ→0,η→-1,τ→0,則(4)式退化為一階速度-壓力形式的純聲波方程。

1.1.2 數(shù)值計算方法

采用交錯網(wǎng)格有限差分進行波場模擬,對于(4)式中的高階項采用偽譜法進行計算。其公式為:

(5)

波場傳播到計算區(qū)域的邊界,采用最佳匹配層(PML)吸收邊界條件。由于PML是在非吸收衰減方程中得到的邊界條件,本質(zhì)上也是對波場振幅的衰減。為了避免兩種衰減效應剛性過渡造成的邊界反射,在計算區(qū)與PML之間添加一個過渡區(qū),在過渡區(qū)中兩種吸收衰減線性過渡。

1.2 Hessian核函數(shù)的求解與求逆

1.2.1 點擴散函數(shù)與其逆函數(shù)的求取

在背景介質(zhì)中,記炮點位置為xs,檢波點位置為xr,由炮點傳播到散射點的頻率域格林函數(shù)定義為:

(6)

式中:σ(x)=1/v(x),為介質(zhì)點xs上的慢度,v(x)為介質(zhì)點集合x上的速度。類似的可以定義由檢波點到散射點的格林函數(shù)。在Born近似下,觀測系統(tǒng)的地震記錄可以表示為:

(7)

式中:fs(ω)為炮點xs處激發(fā)的能量在圓頻率ω上的譜;r是一個與反射系數(shù)相關(guān)的量。用矩陣-向量記法,(7)式可寫為:

(8)

式中:L表示一個線性模擬算子,與觀測系統(tǒng)、震源子波和地下介質(zhì)模型參數(shù)有關(guān);m0是地下反射率模型。對應正演模擬算子L,可以獲得其伴隨偏移算子LT,則偏移成像結(jié)果m可表示為:

(9)

(9)式用格林函數(shù)表示為:

(10)

式中:Gt表示G的共軛函數(shù);m(x0)為位于x0處的點散射體;Ha(x,y)為Hessian算子,其表達式為:

(11)

可以看出,Hessian矩陣的元素Ha(x,y)對應兩組成像點,分別為x和y。Hessian矩陣的線性項表達成了格林函數(shù)的函數(shù)。而線性Hessian矩陣的每一個元素對應地下的兩個成像點,即線性Hessian矩陣的元素個數(shù)為成像點個數(shù)的平方,這個數(shù)量再加上矩陣求逆運算,計算量巨大。Hessian核函數(shù)的敏感性能夠更加體現(xiàn)在對地下成像空間每個點的照明。抽取線性Hessian矩陣的一行,即點擴散函數(shù):

(12)

(12)式描述了對應于一個成像點x處的點擴散函數(shù),可以看到單點的PSF彌散到整個成像空間的所有點y,反映了觀測系統(tǒng)對該點的照明能力。正向傳播的兩個格林函數(shù)G(xs,x,ω)和G(x,xr,ω)描述了地震波場從震源傳播到點x處經(jīng)由散射傳播到接收點的過程;而反向傳播(共軛即是反向傳播)的兩個格林函數(shù)Gt(xs,y,ω)和Gt(y,xr,ω)描述了地震波場從接收點反傳到y(tǒng)處經(jīng)由散射反傳到震源點的過程。而整個正傳、反傳的過程恰好反映了地震波在背景介質(zhì)中傳播對于散射點的聚焦程度。

針對單個散射體來說,其點擴散函數(shù)響應即為其偏移成像結(jié)果。GUITTON[23]提出了一種局部去模糊濾波器構(gòu)建方法,直接在空間域執(zhí)行去模糊偏移成像。給定參考模型m0,觀測數(shù)據(jù)uobs,可得到常規(guī)偏移成像結(jié)果m1:

(13)

對該結(jié)果進行反偏移,得到新數(shù)據(jù):

(14)

對新數(shù)據(jù)u1再次進行成像,獲得新的偏移結(jié)果m2:

(15)

其中,兩次成像結(jié)果之間的關(guān)系為:

(16)

也就是說,m2是m1經(jīng)過Hessian算子模糊作用后的成像結(jié)果。

假設(shè)采用一組非平穩(wěn)濾波器來解決以下優(yōu)化問題:

(17)

則B就是(LTL)-1的一個最優(yōu)化估計。將估計得到的B作用到成像結(jié)果m1上,得到去模糊化的成像結(jié)果。

1.2.2 解析表達的PSF

聲學格林函數(shù)可以表示為如下解析表達式:

(18)

如果為黏性介質(zhì),則將完全彈性介質(zhì)中的實速度v0替換成復速度v(ω):

(19)

從而得到黏聲介質(zhì)的格林函數(shù):

(20)

其中,ξ={1+(1/πQ)[ln(ω/ω0)]}[1+(1/4Q2)]。因此,可據(jù)此給出聲學和黏聲情況下的PSF:

(21)

利用與聲介質(zhì)情況相同的去模糊濾波器來求取黏聲介質(zhì)點擴散函數(shù)的逆,然后作用于常規(guī)偏移成像結(jié)果上就能達到去模糊化的效果。

2 模型參數(shù)實驗

2.1 黏聲介質(zhì)復雜模型逆時偏移成像測試

針對簡單模型的逆時偏移我們設(shè)計了一個凹陷速度模型,如圖1所示,模型大小為1500m×2000m,水平、垂直網(wǎng)格間距都是5m。速度范圍為3500～5000m/s。

圖1 凹陷速度模型

首先,我們設(shè)Q值為常數(shù),分別對該模型整體設(shè)置Q=5,10,30,100和無窮大,圖2為不同Q值對應的地震記錄。由圖2可以看出,Q值的變化會造成不同程度的振幅衰減,且Q值越小衰減越嚴重。從圖2還可以看出,由于Q值的影響,地震記錄有向下的時移,這表明Q值對地震波傳播具有相位畸變效應。圖3給出了不同Q值下的同一單道地震記錄。從圖3 可以看出,Q值除了引起振幅的衰減,還使得相位發(fā)生了變化:Q越小相位越滯后。

圖2 不同Q值對應的地震記錄

圖3 不同Q值的同一單道地震記錄

考慮地質(zhì)構(gòu)造特征,我們重新設(shè)置Q值,使得Q值結(jié)構(gòu)與速度模型結(jié)構(gòu)相對應,模型的Q值設(shè)置如圖4所示。模擬100炮地震數(shù)據(jù),炮點位置從450m開始,炮間距為15m。地表位置布滿接收點,總道數(shù)為400道,道間距5m,記錄長度0.88s,時間采樣率1ms。

圖5給出了采用不同數(shù)據(jù)和成像方法得到的成像結(jié)果。

圖4 凹陷Q模型

圖5a為聲介質(zhì)成像結(jié)果;圖5b為正演的衰減數(shù)據(jù)不加以補償?shù)哪鏁r成像結(jié)果;圖5c為正演的衰減數(shù)據(jù)給予補償?shù)哪鏁r成像結(jié)果。由圖5可以看出,與聲介質(zhì)成像結(jié)果相比,偏移時如果不對衰減數(shù)據(jù)進行Q補償,地震波能量明顯減弱,圖像分辨率變低,衰減層越到下方成像效果越差,圖像層位也不準確,層位略向上移動。當采用本文衰減與補償統(tǒng)一表達的黏聲介質(zhì)傳播方程對帶吸收衰減數(shù)據(jù)進行補償成像后,成像結(jié)果的振幅得到補償,下方層位被有效地恢復出來,圖像層位也準確。該實驗結(jié)果證明了吸收與衰減統(tǒng)一表達的黏聲介質(zhì)波動方程能夠有效地計算黏聲介質(zhì)地震波場的格林函數(shù)。

2.2 黏聲介質(zhì)PSF特性分析

圖6為鹽丘速度模型,模型大小為6000m×750m,水平、垂直網(wǎng)格間距都是5m?？紤]常規(guī)地質(zhì)構(gòu)造中,Q值結(jié)構(gòu)與速度模型結(jié)構(gòu)相對應,我們設(shè)置模型中的Q值分布如圖7所示。模擬325炮數(shù)據(jù),炮間距15m。中間放炮,兩邊接收,每一炮覆蓋總道數(shù)為200道,道間距5m,記錄長度2.5s,時間采樣率0.3ms。用數(shù)值實驗來觀測線性Hessian矩陣的單點響應。

圖5 不同數(shù)據(jù)和成像方法得到的成像結(jié)果a 聲介質(zhì)成像結(jié)果; b 正演的衰減數(shù)據(jù)不加以補償?shù)哪鏁r成像結(jié)果; c 正演的衰減數(shù)據(jù)給予補償?shù)哪鏁r成像結(jié)果

按照模型數(shù)據(jù)的觀測系統(tǒng)和頻帶分布,我們根據(jù)公式(21)計算其點擴散函數(shù)。將模型劃分為231個網(wǎng)格大小為21×21的區(qū)域,目標點位于網(wǎng)格區(qū)域中心。圖8給出了圖6中6個參考點的點擴散函數(shù)圖(1～6號參考點的中心坐標分別為:[436,10],[562,31],[751,31],[457,73],[667,73],[688,94])。從圖8中可以看出,點擴散函數(shù)集中在中心點周圍區(qū)域。在這一區(qū)域之外,單點響應相對區(qū)域內(nèi)為極小值。另外,由于淺層速度結(jié)構(gòu)簡單,同時照明也較為均勻,它們的PSF也是簡單的橢圓形。而在速度變化劇烈的區(qū)域,PSF卻是畸變的,而且具有一定的方向特性。對比圖8中聲介質(zhì)單點PSF與黏聲介質(zhì)PSF可以發(fā)現(xiàn):吸收衰減效應顯著降低了照明振幅強度,帶吸收衰減的PSF明顯比聲介質(zhì)PSF能量低,在深層這種能量損耗更為強烈;吸收衰減效應改變了觀測系統(tǒng)對地下的照明圖樣,即各點在兩種情況下的分辨率相差較大。

在中心點位置將PSF圖樣按照空間位置進行排列,將全觀測系統(tǒng)的PSF圖樣顯示在圖9(聲介質(zhì))和圖10(黏聲介質(zhì))中。

圖6 SEG/EAGE鹽丘模型

圖7 鹽丘Q模型

從圖9和圖10可以看出,吸收衰減效應的影響體現(xiàn)在全空間中,其顯著降低了鹽丘下方照明振幅。也可以看到,線性Hessian矩陣在不同區(qū)域展示的單點響應不同。淺層速度結(jié)構(gòu)簡單,同時照明也較為均勻,它們的PSF是簡單的橢圓形;而在速度變化劇烈的區(qū)域PSF發(fā)生了畸變,且有一定的方向特性,圖樣更為分散;在鹽丘周圍,吸收衰減效應顯著降低了成像分辨能力。

2.3 基于點擴散函數(shù)的反演測試

圖11和圖12分別給出了聲介質(zhì)和黏聲介質(zhì)的逆PSF成像結(jié)果。圖13給出了常規(guī)聲介質(zhì)偏移成像結(jié)果。雖然沒有Q的影響,但在高速體下方,常規(guī)偏移成像并不能得到顯著的層位信息。圖14給出了逆PSF作用在常規(guī)聲介質(zhì)偏移成像結(jié)果后的成像結(jié)果?？梢钥闯?逆PSF作用在常規(guī)聲介質(zhì)偏移成像結(jié)果上,能夠達到去模糊的作用,特別是在鹽丘下方,能夠有效恢復層位信息。

圖8 圖6中6個參考點的PSF分布(左邊為聲介質(zhì)的PSF分布,右邊為黏聲介質(zhì)的PSF分布)a 參考點1; b 參考點2; c 參考點3; d 參考點4; e 參考點5; f 參考點6

圖9 全觀測系統(tǒng)點擴散函數(shù)分布(聲介質(zhì))

圖10 全觀測系統(tǒng)點擴散函數(shù)分布(黏聲介質(zhì))

圖15給出了采用黏聲數(shù)據(jù)并加以補償?shù)某Ｒ?guī)偏移成像結(jié)果。由圖15可以看出,衰減層及其下部反射體的成像振幅逐漸減弱,成像分辨率逐漸降低。圖16給出了逆PSF作用在圖15 成像結(jié)果后的成像結(jié)果。對比圖15和圖16可以看出,逆PSF顯著增強了深層成像的能量,成像剖面振幅更加均衡,吸收衰減效應能夠通過逆PSF進行有效補償。

圖11 聲介質(zhì)逆PSF成像結(jié)果

圖12 黏聲介質(zhì)逆PSF成像結(jié)果

圖13 直接偏移成像結(jié)果(聲波數(shù)據(jù)+聲波成像)

圖14 逆點擴散函數(shù)作用后的成像結(jié)果(聲波數(shù)據(jù)+聲波成像+聲波PSF)

圖15 直接偏移成像結(jié)果(黏聲數(shù)據(jù)+黏聲成像)

圖16 逆點擴散函數(shù)作用后成像結(jié)果(黏聲數(shù)據(jù)+黏聲成像+黏聲PSF)

3 結(jié)論

本文從衰減與補償統(tǒng)一表達的黏聲介質(zhì)傳播方程出發(fā),利用交錯網(wǎng)格有限差分方法和最佳匹配層吸收邊界條件進行求解,實現(xiàn)了黏聲介質(zhì)衰減與補償?shù)母窳趾瘮?shù)計算?；诘卣鸱囱莩上窭碚摵蛿?shù)值實驗結(jié)果,對點擴散函數(shù)的敏感性以及Q因子對點擴散函數(shù)數(shù)學、物理特征的影響進行了分析。結(jié)果表明,PSF反映了吸收衰減效應對波場穿透能力和角度照明能力的影響。利用去模糊濾波器對PSF進行求逆,并將逆點擴散函數(shù)直接作用在成像結(jié)果上,從而對原始成像結(jié)果有效地去模糊化,并提高成像振幅的分辨率和均衡性。模型數(shù)據(jù)實驗結(jié)果證明了該方法的有效性,在黏聲介質(zhì)中的逆時偏移成像可以達到與聲介質(zhì)成像相當?shù)木取；谶@一系列研究成果,可以進一步開展迭代的最小二乘偏移成像和全波形反演方法研究。

致謝:感謝同濟大學海洋與地球科學學院波現(xiàn)象與反演成像研究組(WPI)對本研究工作的幫助。