高俊杰, 張 進,2,3??, 彭陽陽, 董博藝
(1. 中國海洋大學海洋地球科學學院, 山東 青島 266100;2. 青島海洋科學與技術(shù)試點國家實驗室 海洋礦產(chǎn)資源評價與探測技術(shù)功能實驗室, 山東 青島 266237;3.中國海洋大學 海底科學與探測技術(shù)教育部重點實驗室, 山東 青島 266100)
在海上地震勘探中,震源和檢波器往往沉放在海面以下,這是由海洋地震勘探特殊施工條件所決定的,當震源激發(fā)后,作為強反射層的海面或海底,將會把地震波反射回來并形成虛反射。在海上作業(yè)時,虛反射的廣泛發(fā)育嚴重影響地震資料的質(zhì)量,需加以消除。常規(guī)的資料處理工作中,海面往往被假設(shè)為水平且反射系數(shù)簡單定義為-1,實際上,海面呈現(xiàn)出起伏不平的復(fù)雜變化,導(dǎo)致實際地震記錄中含有復(fù)雜的虛反射,常規(guī)的水平海面假設(shè)將會在反射系數(shù)和延遲時間的計算過程中引入誤差,這就導(dǎo)致了虛反射算子估計不準確,影響地下儲層的最終成像。
很多學者在模擬起伏海面方面做了大量的工作,以期得較為準確描述起伏海面反射系數(shù)。1983年Jovanovich等[1]提出了粗糙海面下的高斯曲面相干系數(shù)表達式,并用該相干系數(shù)來代替海面反射系數(shù);1988年Macaskill[2]提出了粗糙海面下的平面波二維散射公式;2002年Laws和Kragh[3-4]分析時移地震成像觀測結(jié)果后認為起伏海面造成的成像誤差不可忽略;2005年Saenger[5]提出了隨機海面下描述聲波散射近似式;Robert等[6-7]隨后在2006年又提出了用低頻信息重構(gòu)海面的方法;2015年Egorov等[8]提出了近似計算粗糙海面反射系數(shù)的積分表達式;2016年Zhang等[9]將海面反射系數(shù)看成是頻率和波高的函數(shù);2016年P(guān)ark等[10]利用矩量法計算隨風速變化的海面反射系數(shù),并檢驗了模型準確性和有效區(qū)域;2017年Zhang等[11]通過高斯統(tǒng)計公式計算粗糙海面的反射系數(shù)實現(xiàn)粗糙海面虛反射壓制;2021年Lv等[12]提出了一種基于環(huán)境參數(shù)和實時測量數(shù)據(jù)的淺海海面模型。上述方法重構(gòu)的海面海浪形態(tài)不能完全符合海浪的統(tǒng)計學觀測規(guī)律,與真實海面形態(tài)存在一定誤差,從而影響反射系數(shù)和延遲時間的計算精度。
為了更真實地模擬起伏海面,本文引入海浪譜概念,基于文圣常等[13]提出并改進的文氏海浪譜(下文簡稱文氏譜),采用有效波高、海浪成長系數(shù)等各項有效海浪參量,建立符合中國海域海浪起伏特征的起伏海面模型;在滿足差分條件下,使用有限差分交錯網(wǎng)格方法進行正演模擬,通過時頻分析方法描述起伏海面虛反射的振幅特征和頻率特征。
根據(jù)傳播路徑,虛反射可分為直達波虛反射、激發(fā)虛反射、接收虛反射、激發(fā)-接收虛反射[14]四種,如圖1所示。
圖1 一次反射與虛反射傳播路徑示意圖
通常情況下,直達波虛反射在數(shù)據(jù)預(yù)處理過程中就被切除,所以實際地震記錄關(guān)注三種虛反射,分別是激發(fā)虛反射、接收虛反射、激發(fā)-接收虛反射。
用s(t)表示實際地震記錄,則:
s(t)=u(t)+Ru(t-Δts)+Ru(t-Δtr)+R2u(t-Δts-Δtr)。
(1)
式中:u(t)為一次反射波;R為反射系數(shù);Δts、Δtr和(t-Δts-Δtr)分別表示激發(fā)、接收和激發(fā)接收虛反射延遲時間;u(t-Δts)和u(t-Δtr)為激發(fā)虛反射和接收虛反射。通過傅氏變換,式(1)為:
s(t)=U(f)(1+Re-i2πfΔts)(1+Re-i2πfΔtr)=U(f)·G(f)。
(2)
可以發(fā)現(xiàn),含有虛反射的地震記錄S在頻率域相當于一次反射波U(f)和算子G(f)的乘積,其陷波周期主要與虛反射延遲時間Δt相關(guān)。
海面的反射系數(shù)是虛反射模擬以及虛反射壓制處理中的一個重要參數(shù),Jovanovich等[15]結(jié)合經(jīng)典的散射理論并利用Kirchhoff近似推導(dǎo)出了高斯曲面的相關(guān)系數(shù)表達式,并提出該相關(guān)系數(shù)就是起伏海面下的海面反射系數(shù),其表達式如下:
(3)
式中:f為入射波頻率;φ為均方根振幅;θ是入射角;v為海水速度。根據(jù)公式(3),繪制出海面反射系數(shù)隨入射角和頻率的關(guān)系(見圖2),當海面均方根振幅一定時,海面反射系數(shù)絕對值隨入射角減小、地震波頻率增大而減小[15]。在垂直入射條件下,當頻率130 Hz時海面反射系數(shù)為-0.75,誤差可達25%。因此,如果忽略海面起伏不平的復(fù)雜變化,普遍假設(shè)海水面平坦,將在虛反射反射系數(shù)和延遲時間計算中將引入巨大誤差,故海面起伏形態(tài)與海面反射系數(shù)有關(guān),且不可忽略。
圖2 起伏海面反射系數(shù)隨入射角和地震波頻率的變化關(guān)系
早在20世紀50年代初,海浪被認為是由許多振幅、頻率、方向和相位不同的簡單波動的疊加[16],因為這些簡單波的振幅或相位被學者們規(guī)定為隨機量,所以這些簡單波疊加所得的結(jié)果也是隨機量。其中由S(ω)表示角頻率為ω的組成波的能量,則描述有限區(qū)域內(nèi)海浪能量的函數(shù)S(ω)就被稱為海浪譜。
海面起伏形態(tài)一般與水域深度、海面風速、地理位置和海浪狀態(tài)等因素相關(guān),常見的海浪譜主要有Neumann譜(諾依曼譜)、Pierson-Moskowitz譜(P-M譜),ITTC譜和雙參數(shù)譜[17],而上述海浪譜主要來自北大西洋海域,與中國海域有所差異。為彌補上述海浪譜的缺陷,文圣常等結(jié)合中國黃、東、南、渤海資料,對理論風浪頻譜進行改進和更廣泛的檢驗,提出更符合中國海域海面起伏特征的海浪譜:
(4)
圖3 不同風速下的文氏譜(a)和海面起伏模型(b)
為探究起伏海面條件下的虛反射特征,設(shè)置如圖4所示的正演模型,其中海底深度為50 m,海水速度為1 500 m/s,地層速度為2 400 m/s,采用左側(cè)單邊放炮的激發(fā)方式,正演子波為30 Hz雷克子波,時間采樣間隔為0.5 ms,空間采樣間隔為1 m。根據(jù)文氏譜分別設(shè)置不同有效波高的海面進行正演模擬。
圖4 起伏海面正演模型示意圖
本文利用等距差分格式的有限差分方法,建立聲波方程的層狀速度模型進行數(shù)值模擬,并設(shè)置PML層,在兼顧精度與計算效率的同時減少邊界效應(yīng)的影響。
根據(jù)二維均勻介質(zhì)中的彈性波方程可推導(dǎo)出聲波方程[18]為:
(5)
式中:ρ為地層密度;u為應(yīng)力分量;v為速度;vx、vz分別表示在x和z方向上的速度分量。
用交錯網(wǎng)格有限差分方法(見圖5)對式(5)中各項偏導(dǎo)數(shù)進行離散[19-21],如式(6):
圖5 有限差分方程的交錯網(wǎng)格示意圖
(6)
式中:Δx、Δz分別表示 X 和 Z 方向上的空間網(wǎng)格步長;i、j分別表示相應(yīng)的空間網(wǎng)格點;Δt表示時間網(wǎng)格步長;k表示相應(yīng)的時間網(wǎng)格點;用U、V、W來分別表示應(yīng)力、速度分量vx、vz的離散形式。
本文基于文氏譜起伏海面建模和有限差分正演模擬方法,分別探究起伏海面虛反射對地震記錄的改造作用、起伏海面條件下的虛反射時頻譜特征。
為探究起伏海面虛反射對地震記錄的改造作用,保持震源沉放深度為水下3 m,檢波器沉放深度為水下10 m,分別設(shè)置海面有效波高H為0 m(水平海面)、2 m(低起伏海面)、3 m(中起伏海面)和4 m(高起伏海面)進行正演模擬。為符合實際海洋地震勘探現(xiàn)狀,本次模擬中采用30 Hz雷克子波作為震源子波,時間采樣間隔為1 ms。由此觀測方式計算的激發(fā)、接收虛反射的延遲時間差值小于單位時間間隔,激發(fā)虛反射、接收虛反射交混在一起難以分離,故本文分析虛反射的整體響應(yīng)。
由不同有效波高起伏海面模型的地震單炮正演記錄(見圖6)發(fā)現(xiàn):
與水平海面模型的地震記錄(見圖6(a))相比,當海面為起伏形態(tài)時,記錄中各同相軸粗糙程度明顯增加,出現(xiàn)抖動現(xiàn)象。隨著海面起伏增大,地震記錄同相軸更加粗糙,能量分布不均勻程度增加,同相軸難以分辨(見圖6(b)—(d))。
造成地震記錄同相軸“抖動”的原因為:對于直達波,當炮檢距較大時,直達波的入射和反射平行于海面,故與虛反射發(fā)生相對規(guī)則的干涉,海面起伏的影響較小;對于反射波,當海面為起伏狀態(tài)時,相對直達波其入射方向于海面的夾角角度更大,故一次反射波和虛反射發(fā)生干涉的不規(guī)則程度相對增加,受海面起伏的影響更大。
為探究不同海面起伏程度對地下反射層成像效果的影響,將多個單炮記錄進行動校正、疊加等處理,獲得成像剖面(見圖7)如下:
((a)水平海面;(b)低起伏海面(H=2 m);(c)中起伏海面(H=3 m);(d)高起伏海面(H=4 m)。 (a) Horizontal sea surface; (b) Slightly undulating sea surface (H=2 m); (c) Moderately undulating sea surface (H=3 m); (d) Heavily undulating sea surface (H=4 m).)
對比水平、低、中、高起伏海面條件地下反射層成像結(jié)果可見:相對于水平海面,起伏海面記錄有著較為粗糙的同相軸,能量分布不均勻,且有可能在實際記錄中出現(xiàn)虛假同相軸。此例中,當起伏海面H=2 m時,地層成像結(jié)果(見圖7(b))與水平海面條件下的成像結(jié)果(見圖7(a))基本一致;當起伏海面H=3 m時,0.55 s、1.17 s處出現(xiàn)了較為明顯的虛反射同相軸(見圖7(c)),在圖7(d)中H=4 m時,起伏海面導(dǎo)致的虛反射更加嚴重,不僅強烈扭曲了海底一次波反射同相軸,且與一次波同相軸發(fā)生了強烈的干涉作用,一次波位置和能量發(fā)生變化(0.95 s),成像剖面出現(xiàn)假同相軸(0.8 s、1.25 s處)。
為進一步探究起伏海面條件下的虛反射時頻譜特征,分別抽取上述四種不同波高起伏海面的單道虛反射記錄(見圖8(a),(c),(e),(g)),通過短時傅里葉變換方法獲得時頻譜(見圖8(b),(d),(f),(h))。
((a)水平海面(H=0 m)單道記錄;(b)水平海面時頻譜;(c)低起伏海面(H=2 m)單道記錄;(d)低起伏海面時頻譜;(e)中起伏海面(H=3 m)單道記錄;(f)中起伏海面時頻譜;(g)高起伏海面(H=4 m)單道記錄;(h)高起伏海面時頻譜。(a) Single-trace record of horizontal sea level (H=0 m); (b) The spectrum of (a); (c) Single-trace record of slightly undulating sea surface (H=2 m); (d) The spectrum of (c); (e) Single-trace record of moderately undulating sea surface (H=3 m); (f) The spectrum of (e); (g) Single-trace record of heavily undulating sea surface (H=4 m); (h) The spectrum of (g).)
通過對比發(fā)現(xiàn):與水平海面單道記錄的時頻分析結(jié)果相比,直達波虛反射能量強度隨著海面有效波高的增大而相對減小。隨著起伏海面有效波高的增加,有效波虛反射能量逐漸增強,與直達波虛反射響應(yīng)時間差變化劇烈;多次波能量強度增大且分布紊亂程度加劇。之所以造成上述現(xiàn)象,除了受海綿起伏影響較大的一次反射波和虛反射的不規(guī)則干涉,還因為海面起伏造成的不均勻照明,加劇了反射能量的空間不均勻分布。
通過在模型四周設(shè)置PML層,可以得到不含虛反射的有效波記錄。在同一個地層模型、相同的震檢組合下,得到含有虛反射的水平海面正演結(jié)果(見圖9(a)紅)與不含虛反射的水平海面有效波正演結(jié)果(見圖9(a)黑),并分別選取偏移距為30、35和40的單道虛反射記錄進行頻譜分析。同理,得到含有虛反射的有效波高2m的起伏海面正演結(jié)果(見圖9(b)紅)與不含虛反射的水平海面有效波正演結(jié)果(見圖9(b)黑),并對不同偏移距的單道記錄進行頻譜分析,探究不同波高起伏海面條件下的虛反射信號頻譜特征。結(jié)果如圖9所示。
通過對比發(fā)現(xiàn):不同海面的正演結(jié)果頻譜中都存在陷波點,水平海面含虛反射記錄頻譜中陷波點分別在17、55 Hz的位置;起伏海面含虛反射記錄頻譜中陷波點分別在19、52 Hz的位置。相較于水平海面,起伏海面頻譜陷波效應(yīng)明顯增強,低頻和高頻端陷波點處能量明顯增強,主頻幅值相對降低。隨著偏移距的增加,水平、起伏海面條件下陷波點呈現(xiàn)相似規(guī)律,即主頻幅值隨著偏移距增加而逐漸減小,主頻位置向高頻方向移動;但起伏海面變化更為劇烈。
為探究不同觀測系統(tǒng)條件下的虛反射時頻譜特征,設(shè)置有效波高為3 m的起伏海面,保持震源沉放深度水下3 m不變,分別設(shè)置檢波器沉放深度為水下4、6、8和10 m進行正演模擬。抽取單道虛反射記錄(見圖10(a),(c),(e),(g)),通過短時傅里葉變換方法獲得時頻譜(見圖10(b),(d),(f),(h))。
((a)沉放4 m單道記錄;(b)沉放4 m時頻譜;(c)沉放6 m單道記錄;(d)沉放6 m時頻譜;(e)沉放8 m單道記錄;(f)沉放8 m時頻譜;(g)沉放10 m單道記錄;(h)沉放10 m時頻譜。(a) Single-trace record of geophone at depth of 4 m; (b) The spectrum of (a); (c) Single-trace record of geophone at depth of 6 m; (d) The spectrum of (c); (e) Single-trace record of geophone at depth of 8 m; (f) The spectrum of (e);(g) Single-trace record of geophone at depth of 10 m; (h) The spectrum of (g).)
由不同檢波器深度單道記錄的時頻分析結(jié)果發(fā)現(xiàn):與檢波器沉放深度為4 m的虛反射單道記錄的時頻分析結(jié)果(見圖10(b))相比,隨著檢波器沉放深度逐漸增加,深檢波器沉放條件下直達波虛反射記錄的能量強度逐漸減小,而有效波虛反射的能量強度逐漸增加。
在同一個地層模型、相同的偏移距下,正演模擬得到含有虛反射的水平海面正演結(jié)果(見圖11(a)紅)與不含虛反射的有效波正演結(jié)果(見圖11(a)黑),分別對檢波器沉放深度為4、6和8 m的單道虛反射記錄進行頻譜分析。同理,正演模擬得到含有虛反射的有效波高3 m的起伏海面正演結(jié)果(見圖11(b)紅)和不含虛反射的水平海面有效波正演結(jié)果(見圖11(b)黑),并對不同沉放深度的單道記錄進行頻譜分析,探究不同觀測系統(tǒng)條件下的虛反射頻譜特征,結(jié)果如圖11所示。
((a)水平海面含虛反射記錄(紅)、水平海面無虛反射記錄(黑);(b)起伏海面含虛反射記錄(紅)、水平海面無虛反射記錄(黑)。(a) under condition of horizontal sea surface with ghost wave (red), and without ghost wave (black); (b) under condition of undulating sea surface with ghost wave (red), and without ghost wave (black).)
通過對比發(fā)現(xiàn):不同海面的正演結(jié)果頻譜中都存在陷波點,水平海面含虛反射記錄頻譜中陷波點分別在13、64 Hz的位置;起伏海面含虛反射記錄頻譜中陷波點分別在17、64 Hz的位置。相較水平海面,起伏海面頻譜的陷波效應(yīng)明顯增強,且頻寬變窄。隨著檢波器沉放深度的增加,起伏海面條件下頻譜受改造程度加劇,幅值波動明顯增加。根據(jù)震源、檢波器空間互易原理,激發(fā)虛反射與接收虛反射規(guī)律一致。
本文設(shè)置海浪狀態(tài)為充分成長型海浪,該狀態(tài)下由于其內(nèi)部各組成波之間能量的非線性輸送平衡,即使隨著時間的變換,海浪能量的頻率分布仍具有很好的相似性[22]。經(jīng)過實驗,不同時間、相同狀態(tài)的海面正演記錄所獲結(jié)論一致。
本文利用文氏譜,建立了更符合真實海面形態(tài)的起伏海面,并基于交錯網(wǎng)格有限差分方法,正演模擬得到了不同有效波高起伏海面地震記錄,通過短時傅里葉變換方法進一步研究了起伏海面條件下的虛反射特征,在研究過程中,得到了以下的結(jié)論和認識:
(1)當海面起伏時,虛反射常常與地下一次波發(fā)生干涉、疊加。當海面起伏程度較大時,海面虛反射對地震記錄的改造作用強烈,不僅強烈扭曲了海底一次波反射同相軸,還導(dǎo)致地層成像結(jié)果中出現(xiàn)了很多假同相軸,嚴重影響了地震資料的信噪比和分辨率。
(2)當海面起伏時,虛反射和多次波能量強度增大、分布紊亂程度增加且無法忽略。地震記錄受陷波效應(yīng)逐漸增強,主頻幅值逐漸減小,虛反射能量逐漸增加。
(3)隨著偏移距的增加,地震記錄陷波點幅值變化豐富,主頻幅值逐漸減小,主頻位置向高頻方向移動。
(4)隨著震源或檢波器沉放深度逐漸增加,海面起伏對信號陷波點位置和能量的改造明顯。起伏海面有效波虛反射的能量強度逐漸增加,陷波作用加劇,主頻幅值逐漸減小。