疊前衰減補償時間偏移及GPU實現(xiàn)

吳吉忠 左 虎

(中國石油冀東油田公司，河北唐山 063004)

0 引言

地下介質(zhì)的黏滯特性會導(dǎo)致地震波能量耗散、相位畸變、地震剖面分辨率降低?；谄焚|(zhì)因子Q與頻率無關(guān)或隨頻率變化弱的假設(shè)[1-2]，一般采用反Q濾波[3-5]或黏聲波方程疊前深度偏移[6-7]兩種方法進行衰減補償。對疊后剖面應(yīng)用反Q濾波時，沒有考慮地震波真實的傳播路徑，補償精度不能得到保證；疊前深度偏移在波場延拓中利用黏聲波動方程模擬地震波的傳播，精度高但層Q模型求取十分困難。

關(guān)于Q值的求取，前人已做了大量研究[8-12]。常用譜比法和頻移法基于VSP數(shù)據(jù)求取地層的Q值，但對于地面地震而言，受薄互層調(diào)諧效應(yīng)及資料信噪比影響，Q值求取精度仍需進一步提高。在利用地震數(shù)據(jù)與井?dāng)?shù)據(jù)聯(lián)合求取Q上，業(yè)界也有不少嘗試，但在兩種數(shù)據(jù)生成的Q值標(biāo)定融合方面沒有形成統(tǒng)一的認識與標(biāo)準。

對于橫向速度變化不是很劇烈的復(fù)雜構(gòu)造，疊前時間偏移是一種有效的成像方法。成像過程中偏移孔徑直接影響疊前時間偏移的信噪比，偏移孔徑過小，成像信噪比較高，但大傾角構(gòu)造面臨成像不足的風(fēng)險；過大的偏移孔徑可以保證陡傾構(gòu)造成像但是增加了偏移噪聲和計算量。

通常偏移孔徑是根據(jù)地層傾角估算的，但在成像前獲取準確的地層傾角十分困難。基于地層傾角與偏移孔徑的關(guān)系，很多學(xué)者利用地層傾角約束偏移孔徑，進而達到改善成像質(zhì)量的目的。一般采用傾角掃描的方式獲取地下地層傾角，如果掃描時窗內(nèi)有多個同相軸，且每個同相軸對應(yīng)的傾角不同，則掃描結(jié)果會出現(xiàn)不穩(wěn)定的現(xiàn)象。此外傾角掃描對數(shù)據(jù)的信噪比也有一定要求，信噪比較低掃描結(jié)果誤差較大。有的學(xué)者利用地震波的傳播路徑計算地震反射傾角，但通過計算方式直接獲取地層傾角的研究目前較少[13-17]。

為了提高疊前時間偏移成像的計算效率，很多學(xué)者將目光轉(zhuǎn)向了GPU技術(shù)。李博等[18]將非對稱旅行時疊前時間偏移方法移植到了GPU平臺，取得了較好的加速效果；李肯立等[19]利用GPU對疊前時間偏移進行了移植和效率測試；馬召貴等[20]利用GPU技術(shù)對Kirchhoff疊前時間偏移進行了優(yōu)化，計算效率明顯提升。與常規(guī)疊前時間偏移不同，衰減補償?shù)寞B前時間偏移需要在頻率域完成每個頻率點的補償，計算量巨大，因此GPU加速對衰減補償疊前時間偏移尤為重要。

本文通過引入等效Q參數(shù)實現(xiàn)衰減補償疊前時間偏移，并針對現(xiàn)有技術(shù)求取Q場存在的弊端、偏移孔徑難以確定及偏移方法效率低的難題，給出了針對性的解決方案。在Q場建模方面，通過VSP資料與地面反射地震資料聯(lián)合應(yīng)用求取Q值，并基于衰減合成地震記錄與井旁地震道的匹配情況判斷Q值是否合理；在壓制偏移噪聲方面，通過孔徑自適應(yīng)改變，有效壓制偏移噪聲；在提高計算效率方面，采用GPU加速，有效解決了頻率域吸收補償疊前時間偏移計算效率低的問題。

1 方法原理

1.1 井震聯(lián)合求取Q場

由于在實際生產(chǎn)中不可能通過巖心樣本得到所有地層的真實Q值，對Q值的合理性判斷也就缺乏客觀標(biāo)準，因此根據(jù)Q值的應(yīng)用目的選取判斷標(biāo)準更具有現(xiàn)實意義。如果Q是用于黏滯介質(zhì)疊前時間偏移補償?shù)卣鸩ㄋp，那么Q值的選取應(yīng)由井上的衰減合成地震記錄與井旁原始地震道是否匹配來決定。也有學(xué)者基于“衰減補償”模式，即通過不同Q值的補償結(jié)果與無衰減合成地震記錄是否匹配判斷Q值的合理性，但對于中深層低信噪比地震數(shù)據(jù)，為了壓制高頻補償噪聲而經(jīng)常使用的增益控制因子對頻帶與振幅增益都做了限制，導(dǎo)致中深層補償效果較差，與合成地震記錄匹配效果差，直接影響了中深層Q值的選取[12]。本文利用井上衰減合成地震記錄與井旁原始地震道匹配效果，選擇合理的Q值，運算過程沒有用到增益控制因子，從根本上避免了上述問題的出現(xiàn)。

為了方便研究，引入了等效Q值的概念[21-22]，并假定地下某點的吸收只與該點的等效Q值有關(guān)，而不受上覆地層吸收的影響。等效Q值的定義為

(1)

式中：T代表地震波垂直單程旅行時；ΔTi為地震波在第i層內(nèi)垂直單程旅行時；Qi為第i層的Q值；n為地層層數(shù)。

利用VSP數(shù)據(jù)與地面地震數(shù)據(jù)聯(lián)合求取Q值的思路是：首先利用VSP數(shù)據(jù)求取Q值，然后以該Q值作為初始值，生成衰減合成地震記錄，并根據(jù)該合成地震記錄與井旁道原始地震數(shù)據(jù)匹配情況調(diào)整優(yōu)化Q值，然后采用對數(shù)譜比求導(dǎo)算法求取地震數(shù)據(jù)的Q值，最后將兩種Q值標(biāo)定融合得到整個工區(qū)的Q值。

通過分析VSP數(shù)據(jù)直達波的頻率與振幅屬性，采用質(zhì)心頻移法[21]計算Q值

(2)

式中：τ為地震波在地層內(nèi)的傳播時間；σ2是震源子波方差；fshot和fgeo分別是炮、檢點地震波主頻。

衰減合成記錄的生成分為兩步：首先利用雷克子波生成不考慮衰減的常規(guī)合成記錄，在確定子波主頻時，可以根據(jù)淺層地震數(shù)據(jù)的主頻及目的層對分辨率的需求聯(lián)合確定；然后基于VSP數(shù)據(jù)生成的Q值，利用正Q濾波算法對該合成記錄進行處理，得到了考慮衰減的合成記錄，其中VSP生成的Q值在縱向上是連續(xù)變化的，不是一個常數(shù)。上述合成記錄生成方法的特點是雷克子波主頻在時間方向不發(fā)生變化，而通過Q值的變化描述地震波的衰減。另外還可以通過主頻時變的雷克子波直接生成吸收衰減合成記錄，但面臨著主頻不易準確求取的巨大困難，因此不建議使用該方法。

根據(jù)衰減合成記錄與井旁道原始地震數(shù)據(jù)的相關(guān)性開展局部Q值優(yōu)化與調(diào)整，并迭代更新吸收衰減合成記錄，再次進行相關(guān)性計算，直至衰減合成記錄與井旁道原始地震數(shù)據(jù)達到最佳相關(guān)時停止迭代計算。

在利用地面地震數(shù)據(jù)求取Q值時，采用了對數(shù)譜比求導(dǎo)的計算方法[22]，與常規(guī)對數(shù)譜比法相比，該方法克服了地震數(shù)據(jù)的薄層調(diào)諧效應(yīng)的不利影響。首先采用Q掃描方法對參考點(假設(shè)沒有衰減)的地震記錄進行衰減，每個Q對應(yīng)的時域波形短時傅里葉變換模為FQ(ω)，地下樣本點對應(yīng)的時域波形短時傅里葉變換模為F(ω)，則譜比頻率導(dǎo)數(shù)指標(biāo)為

(3)

當(dāng)δ絕對值為零或最小時，表明FQ(ω)與F(ω)最接近，便認為對應(yīng)的Q是最合理的。這種Q值求取方法與以往相比有一定改進與不同[11]。前人是對地下樣本點進行“Q掃描+衰減補償”處理，需要對疊前數(shù)據(jù)進行多輪次成像，計算量大，又由于增益控制因子的作用，低信噪比數(shù)據(jù)Q值求取穩(wěn)定性差，另外在求取Q值過程中沒有用到井信息，只利用了單一的疊前地震數(shù)據(jù)。而本文方法只對參考點進行“Q掃描+衰減”處理，對于疊后地震剖面而言，參考點選取若干即可，計算量大大降低。雖然該方法是在疊后開展的，但基于地震資料求取的Q值在空間上的趨勢是有保證的，最后利用VSP數(shù)據(jù)計算的Q值對地面地震數(shù)據(jù)生成的Q場進行標(biāo)定融合，從而又保證了Q值的數(shù)值精度。

1.2 衰減補償疊前時間偏移

與常規(guī)Kirchhoff積分法疊前時間偏移不同，衰減補償疊前時間偏移不僅考慮了波場走時和球面擴散效應(yīng)，還考慮了地震波傳播時的吸收效應(yīng)。

本文假設(shè)Q值與頻率無關(guān)或隨頻率緩慢變化。單道地震數(shù)據(jù)在頻率域可用f(ω)表示，基于深度偏移的相移法，檢波點波場的深度延拓可表示為

P(kx,ω,z+Δz) =P(kx,ω,z)exp(-ikzΔz)

=f(ω)exp(-ikzΔz)

(4)

將頻散關(guān)系

(5)

代入式(4)，有

P(kx,ω,z+Δz)

(6)

式中c(ω)為復(fù)速度，與實速度v(ω)關(guān)系為

(7)

其中ωc為高截頻。當(dāng)頻率趨向于ωc時，實速度v(ω)將趨近于常值，本文用主頻ω0代替高截頻ωc，因為主頻對應(yīng)的實速度更容易獲得。經(jīng)過變換消去ωc項，并代入式(7)，可得

(8)

將式(8)代入式(6)，簡化可得

P(kx,ω,z+Δz)

(9)

P(kx,ω,T)≈f(ω)×

(10)

式中：vi為第i層主頻對應(yīng)的實速度；ΔTi=Δzi/vi。定義等效速度為

(11)

根據(jù)式(1)和式(11)，有

(12)

將式(12)代入式(10)，做空間逆傅里葉變換，可得

(13)

式中px=kx/ω為x方向的射線參數(shù)。應(yīng)用穩(wěn)相點原理求得上式漸近解[23]

(14)

上式為檢波點到成像深度的反向傳播表達式。同樣可得炮點到成像深度點的正向傳播波場

(15)

(16)

將成像點所有的成像結(jié)果進行累加便可得到黏滯介質(zhì)中衰減補償疊前時間偏移的成像結(jié)果[24]。

1.3 傾角域自適應(yīng)偏移孔徑

偏移孔徑是影響疊前時間偏移成像信噪比的重要因素，但實際應(yīng)用中，受限于偏移算法，通常采用統(tǒng)一的偏移孔徑。由于地層傾角是空變的，對某一成像點適合的孔徑，對另一成像點就可能過大或過小。為了解決這個難題，本文采用傾角自適應(yīng)的偏移孔徑。

本文采用輸入道的方式實現(xiàn)衰減補償疊前時間偏移，可以根據(jù)地層傾角確定偏移脈沖響應(yīng)在任一成像深度偏移孔徑的范圍，這種成像方式使偏移孔徑時變與空變成為可能。根據(jù)圖1，以下關(guān)系成立

(17)

(18)

式中：t是成像點對應(yīng)的時間深度，為單程旅行時；L為點M與點P之間的水平距離，也是時間深度t處地層傾角β對應(yīng)的偏移孔徑。聯(lián)立式(17)與式(18)，求解可得

(19)

由上式可以看出，當(dāng)vrms、t、h三個參數(shù)固定時，偏移孔徑L隨著地層傾角的增大而增大。圖2a是一個四層地質(zhì)模型[13]，其中①和③表示兩個傾斜界面，②為一個水平界面；圖2b為該地質(zhì)模型數(shù)據(jù)對應(yīng)的傾角成像道集，其中①、②與③分別代表三個同相軸，與圖2a中的界面編號一一對應(yīng)。在傾角成像道集上，界面③的擬雙曲線頂點對應(yīng)的角度便是真實地層傾角，可以確定傾角成像區(qū)為[A，B]，而常規(guī)方法為了實現(xiàn)構(gòu)造的正確成像，偏移孔徑[C，C′]過大而不夠精確，引入了額外的偏移噪聲，降低了地震數(shù)據(jù)的成像品質(zhì)。

圖1 偏移孔徑與地層傾角的關(guān)系

m為炮點s與接收點r的中點，設(shè)定輸入地震數(shù)據(jù)的半炮檢距為h，β為深度t處的地層傾角，P為深度t處最遠成像道對應(yīng)的成像點，p為成像點P在地面的投影，M為過點m的垂線與過P點的水平線的交點，α是地震波在成像點處的最大入射角

圖2 地質(zhì)模型(a)及其對應(yīng)的傾角成像道集(b)

1.4 處理流程及GPU加速實現(xiàn)

1.4.1 計算流程

吸收衰減補償?shù)寞B前時間偏移計算步驟如下：

(1)常規(guī)預(yù)處理；

(2)在炮檢距初疊剖面上選取若干CDP樣本；

(3)生成樣本CDP對應(yīng)的傾角成像道集，并確定樣本CDP對應(yīng)的傾角成像區(qū)；

(4)利用樣本CDP對應(yīng)的傾角成像區(qū)插值出所有CDP對應(yīng)的傾角成像區(qū)；

(5)在CDP對應(yīng)的傾角成像區(qū)內(nèi)，對疊前數(shù)據(jù)進行衰減補償成像處理；

(6)重復(fù)步驟(5)，得到所有CDP最終的成像結(jié)果。

與常規(guī)疊前時間偏移事先將不同炮檢距組對應(yīng)的偏移孔徑計算出來不同，文中給出的疊前時間偏移方法無需事先計算孔徑，只需對傾角成像區(qū)內(nèi)的地震數(shù)據(jù)進行補償成像即可，在保證構(gòu)造正確成像的同時，最大限度地壓制了偏移噪聲。另外實際成像計算區(qū)域與常規(guī)孔徑對應(yīng)區(qū)域相比有明顯減少，計算效率也得到了一定提高。

與常規(guī)疊前時間偏移相比，衰減補償疊前時間偏移最大的不同在于成像時根據(jù)Q值及地震走時逐頻率點進行衰減補償，計算量巨大，計算效率低。

1.4.2 GPU加速實現(xiàn)

業(yè)界對利用GPU進行加速的研究已經(jīng)十分成熟，GPU的加速性能體現(xiàn)在單個GPU設(shè)備有很多計算核心，并且可以同時發(fā)起多個計算線程，而單個CPU設(shè)備則不具備。對于疊前衰減補償時間偏移而言，實現(xiàn)GPU加速，需要在計算熱點、計算粒度與存儲等方面重點考慮。

GPU設(shè)備適合并行非邏輯性運算，比如常規(guī)的加減乘除，而選擇、判斷等邏輯運算更適合用CPU設(shè)備去處理，因此利用GPU設(shè)備實現(xiàn)程序加速，不是簡單的將程序加載到GPU設(shè)備上運算，而是將程序中計算量最大、耗時最長的部分，在這里稱之為計算熱點，放到GPU設(shè)備上進行加速，因此計算熱點的選取直接決定了整體的加速效果。本文將頻率點的衰減補償運算作為計算熱點，這是因為每個成像點都要在有效頻帶內(nèi)逐頻點做衰減補償運算，計算量與耗時在整個程序中是最大的，并且各頻率點的運算互不影響，沒有先后順序，更適合并行運算。計算粒度是用來描述計算單元大小的指標(biāo)，一個線程可以計算一個頻率點的衰減補償，也可以計算整個成像空間所有頻率點的衰減補償，本文將前者稱之為細粒度計算，后者稱之為粗粒度計算。在細粒度計算中，線程并發(fā)量很大，每個線程的計算量卻很小，線程的數(shù)據(jù)讀取與交換用時比有效計算時間要大得多，總體加速效果不理想。在粗粒度計算中，線程的并發(fā)量太小，GPU中處于計算激活狀態(tài)的單元少，設(shè)備中很多計算單元沒有充分利用起來，負載率低，最終的加速效果也不好。借鑒這兩種計算粒度的不足，本文對成像空間進行區(qū)塊劃分，每個線程只負責(zé)單個成像點的衰減補償運算，既兼顧了線程計算量，又考慮了GPU的線程并發(fā)量，可以獲得最佳的加速效果。

衰減補償?shù)寞B前時間偏移在成像過程中需要疊加速度場、地層傾角及Qeff場等參數(shù)文件，對存儲空間需求較大，每一個成像點計算時都需要讀取這些信息，訪問量大，訪問頻率高。雖然GPU上存儲器種類多，但存儲能力相差較大，本文將這些參數(shù)文件放于存儲速度較慢但空間較大的全局存儲器上。一般情況下，穩(wěn)定性增益控制函數(shù)是采用表驅(qū)動的方式進行運算的，事先將補償值計算出來存于表中，在頻率點補償運算時，直接從表中讀取相應(yīng)的補償值，避免了大量重復(fù)運算，這種存儲換計算的策略應(yīng)用在CPU上時，對計算效率提高貢獻很大，但在GPU上時，會出現(xiàn)存儲讀取問題。該補償表占用內(nèi)存空間較大，運算時所有并發(fā)線程會同時從補償表中讀取補償值，進而出現(xiàn)多個線程排隊讀取同一位置數(shù)值的情況，給計算效率帶來了不利影響。針對該問題，穩(wěn)定性增益函數(shù)不再采用表驅(qū)動的方式，而是以補償因子的形式實時參與運算，這種策略表面上增大了線程計算量，但避免了多線程排隊讀取內(nèi)存的問題，整體計算效率仍有大幅提高。

基于上述計算策略，實現(xiàn)了基于GPU計算架構(gòu)的疊前衰減補償時間偏移(圖3)，與只采用CPU相比，計算效率提高了40倍左右，加速效果十分明顯。

圖3 GPU計算架構(gòu)

2 數(shù)據(jù)應(yīng)用

2.1 模擬數(shù)據(jù)

圖4為二維黏彈模型。基于雙程黏性聲波方程利用頻率域有限差分方法模擬地震數(shù)據(jù)，共51炮，炮點位置從447.75m移動到3147.75m，炮點距為54m。每炮有200道，道間距為18m。

圖5為模擬數(shù)據(jù)常規(guī)疊前時間偏移剖面，應(yīng)用了增益顯示。圖6為衰減補償疊前時間偏移剖面。對比圖5、圖6可知，二者都可以得到正確的成像結(jié)果，但圖6的成像結(jié)果分辨率更高。圖7為x=690m處單道偏移結(jié)果的時頻譜對比，可以看出，衰減補償偏移剖面頻帶寬度從淺到深一致，然而常規(guī)偏移剖面的頻帶卻隨著深度的增加逐漸變窄。

圖4 二維模型

圖5 常規(guī)疊前時間偏移剖面

圖6 衰減補償疊前時間偏移剖面

圖7 x=690m處常規(guī)(a)和衰減補償(b)疊前時間偏移單道時頻譜對比

2.2 實際數(shù)據(jù)

南堡凹陷老爺廟火成巖分布廣泛，導(dǎo)致中深層地震信號能量衰減嚴重[25]。首先利用雷克子波生成了不考慮衰減效應(yīng)的常規(guī)合成記錄(圖8a)；然后基于VSP數(shù)據(jù)求取的Q值，利用正Q濾波算法對該合成記錄進行處理，得到考慮衰減效應(yīng)的合成記錄(圖8b)，將該衰減合成記錄與井旁道原始地震數(shù)據(jù)疊合顯示(圖8c)，根據(jù)二者的匹配結(jié)果對VSP生成的Q值進行局部優(yōu)化與調(diào)整(圖8d)，并迭代更新衰減合成記錄，直至吸收衰減合成記錄與井旁道原始地震數(shù)據(jù)達到最佳匹配。

對工區(qū)內(nèi)其他VSP數(shù)據(jù)重復(fù)該過程，最后利用所有優(yōu)化調(diào)整后的Q值去標(biāo)定地震數(shù)據(jù)求取的Q場，得到全工區(qū)的Q場(圖9)。根據(jù)地層分布與Q值的數(shù)值結(jié)構(gòu)，將Q值曲線劃分為三層(圖8d上①、②、③標(biāo)識的橫線)，可以對地震數(shù)據(jù)求取的Q場起到層控約束作用(圖9上①、②、③標(biāo)識的層控約束線)。

基于井震聯(lián)合方法求取的Q場，應(yīng)用衰減補償疊前時間偏移方法對老爺廟實際數(shù)據(jù)進行成像，并與常規(guī)疊前時間偏移及其反Q濾波結(jié)果進行對比(圖10)。反Q濾波沒有考慮地震波實際的傳播路徑，對中高頻信號成分的補償力度要明顯高于低頻信號，這就會造成低頻信號的相對丟失。另外有效信息的恢復(fù)完全取決于輸入的疊后數(shù)據(jù)，并沒有在補償過程中引入新的信息。本文方法在成像過程中圖11是圖10中矩形框局部放大，可以看出衰減補償疊前時間偏移剖面在成像效果上改善較大，細節(jié)信息明顯增多。

圖8 Q值求取過程

圖9 地震資料求取的Q場VSP標(biāo)定結(jié)果對不同傳播路徑的地震信號進行了相應(yīng)補償，有效恢復(fù)了地震信號的高頻能量信息，拓寬了有效頻帶，改善了地震數(shù)據(jù)整體的頻譜形態(tài)，因此在偏移剖面上地質(zhì)信息變得更加真實豐富。另外衰減補償疊前時間偏移剖面的信噪比較高，偏移噪聲與補償噪聲均被控制在一個較低水平。

圖12是三者的頻譜對比，衰減補償疊前時間偏移剖面的頻譜信息更加豐富，低頻與高頻信息均得到了有效拓展，而反Q濾波后雖然高頻信息得到一定恢復(fù)，但低頻信息有所損失，有效頻帶拓寬能力不如本文方法。

圖10 常規(guī)疊前時間偏移剖面(a)及其反Q濾波結(jié)果(b)與衰減補償疊前時間偏移剖面(c)對比

圖11 常規(guī)疊前時間偏移剖面(a)及其反Q濾波結(jié)果(b)與衰減補償疊前時間偏移剖面(c)局部細節(jié)對比

圖12 常規(guī)疊前時間偏移剖面(藍色)及其反Q濾波結(jié)果(綠色)與衰減補償疊前時間偏移剖面(紅色)的頻譜對比

3 結(jié)論

本文提出的衰減補償疊前時間偏移方法，可以在疊前偏移中沿波場傳播路徑補償黏性介質(zhì)的吸收效應(yīng)，提高地震分辨率，且偏移前道集預(yù)處理與常規(guī)偏移方法一致，可以方便地融入現(xiàn)有處理技術(shù)流程。在Q場建模、壓制偏移噪聲和提高計算效率方面，均采用了針對性的技術(shù)策略。通過模型數(shù)據(jù)驗證了方法的有效性與可行性。對于實際數(shù)據(jù)應(yīng)用本文方法得到了分辨率更高、成像更清晰的剖面，與常規(guī)疊前時間偏移剖面及其反Q濾波結(jié)果相比，波組特征更好、地質(zhì)信息更加豐富，因此本文方法在中深層信號衰減強烈且信噪比較低的地區(qū)有著廣闊的應(yīng)用前景。對于Q場的求取，雖然文中給出了一套切實可行的解決方案，但Q的精度依然有較大提升空間，這將是下一步研究的重點。