李福元 韋成龍 鄧桂林 張寶金 張 衡 楊 力
(①自然資源部海底礦產資源重點實驗室,廣東廣州 510760;②自然資源部中國地質調查局廣州海洋地質調查局,廣東廣州 510760)
震源信號子波是反射地震勘探中重要的參數,在地震資料處理及后續很多環節中都起著關鍵作用,包括地震正演模擬、偏移成像、全波形反演、疊前疊后屬性反演、多次波預測和多源數據合成處理等[1-8]。更為一般地說,震源信號是地震勘探信號系統的激勵函數,也是地震波彈性波動方程中的源函數[9]。
研究震源信號的核心問題是獲取震源信號子波。除了采集現場直接觀測記錄外,目前主要的方法有: ①利用近場記錄擬合[10-12]; ②通過專用軟件模擬[13-17]; ③從地震數據本身提取[18]。地震數據采集受現場環境、儀器裝備等因素影響,從現場直接觀測震源子波十分困難。利用近場記錄擬合時,因近場記錄中常伴有干擾,擬合結果往往存在誤差; 在很多情況下地震數據中并沒有近場記錄,故該方法的使用受到限制。軟件模擬方法受軟件所使用理論模型、震源類型及環境參數等因素限制,模擬子波與實際子波存在差異。如果慮及檢波器特性、儀器因素的影響,從地震數據中直接提取震源信號子波應更合理。隨著數字信號處理理論、非線性理論、優化理論的發展,從反射地震數據中提取震源信號子波的方法從理論到實踐都取得了巨大進展[19-20],如何消除反射地震數據中包含的地層信息對提取結果的影響是這類方法成敗的關鍵。
反射勘探地震數據中的直達波,除了表層信息外,不含地下地層信息。如果表層介質是均勻的,則直達波沿直線路徑傳播。在海洋地震勘探中,震源的激發和數據的接收都是在水中進行,海水具有良好的均質性,地震數據中的直達波僅與震源信號、海面反射系數及海水聲波速度有關,因此利用直達波提取震源信號子波是可行的,且無須考慮地層反射信息。
Oliveira[21]提出從地震資料直達波估算震源信號子波的方法,并與近場子波擬合的震源信號做分析對比; 伍忠良等[22]從地震系統的濾波效應、直達波與鬼波傳播路徑及時間差、原始數據的處理對比和理論計算,分析了直達波和鬼波的綜合效應; Koo等[3]研究了一種在拉普拉斯域全波形反演中利用全牛頓法直接從直達波估計震源子波的方法,可用于實際資料震源子波估計; 陳禮等[23]基于深水可完整記錄直達波的特點,在測線近道上應用通過直達波與海底反射波的組合得到統計子波的方法壓制氣泡效應; 鄭江龍等[24]在單道地震記錄中利用從直達波提取地震子波進行相關濾波的方法有效地壓制了干擾波; 任婷等[25]基于深水地震資料利用直達波提取子波對數據進行確定性子波反褶積處理,以有效壓制氣泡效應。這些都從不同角度和層面針對從直達波中提取震源信號子波進行了一定程度的研究及應用。但如何建立震源信號與直達波之間通用的數學關系,進而給出從直達波提取震源信號子波的解析方法,是本文嘗試深入、系統地研究的課題。
海洋地震資料中的直達波除了震源直達波能量外,還包括震源能量的水面反射波,后者即所謂震源鬼波,可看作是位于震源關于水面鏡像點的虛震源的直達波。由于這兩個直達波的傳播路徑不同,加之鬼波還經過水面的反射作用,二者的疊加效應隨炮—檢關系的變化而變化。另外,一般情況下炮檢距遠大于震源和電纜的沉放深度,直達波在接收點是大角度入射,研究直達波必須考慮地震記錄系統的組合效應。顯然,直達波不能直接作為震源信號子波使用。
本文依據氣泡振蕩理論和組合理論,通過推導震源及其虛震源的直達波時距曲線,考慮了震源和接收點組合情況,得到直達波與震源信號之間的關系式,并給出了頻率域利用直達波計算震源信號遠場子波的解析解。當地震數據接收位置滿足震源遠場條件時,利用地震數據中的直達波通過計算提取震源信號遠場子波具有明顯的可操作性和實際意義。實際地震資料處理效果表明,本文方法具有適用性強、抗干擾、計算精度高等特點。利用本文方法提取的地震子波,完全可用于后續的氣泡壓制、虛反射(鬼波)消除、多次波預測、波形反演、偏移成像等眾多處理環節。
海洋地震勘探最常用的是氣槍陣列震源,氣槍在水下激發產生的氣泡快速膨脹形成彈性波能量。根據Keller等[26]提出的氣泡振蕩基礎理論,氣泡在水中被激起后形成一個反復膨脹←→收縮、幅度逐漸衰減的阻尼性振蕩過程; 同時,氣泡本身也受水的浮力而自由上升,直至水面破裂。這種振蕩過程產生一系列脈沖,其中第一次激發膨脹產生的脈沖叫初始脈沖或主脈沖,隨后的多次振蕩膨脹產生數次氣泡脈沖,脈沖的周期與水的密度、靜水壓力、震源能量有關。設計不同脈沖周期的多個震源進行組合激發,使每個單獨震源的初始高能脈沖正相干、而氣泡脈沖為非正相干,就可達到增強初始脈沖、削弱氣泡脈沖的效果。Ziolkowski等[27]和Parkes等[28]研究了氣槍陣列特性,分析了單槍氣泡之間存在的調諧、相干、連通三種關系,指出氣泡在相干半徑外其能量的互相作用呈現線性關系。對于某個檢波點,記錄到的多個單震源的調諧組合子波僅是每個震源子波振幅的簡單疊加,即
(1)
式中:u(t)表示位于檢波點的震源信號的疊加信號;ni(t)表示單震源信號;ri表示檢波點與單震源間的距離;c表示水中聲波速度。
氣泡激發的脈沖能量以彈性波方式在水中傳播,遇到海面(海水與空氣的界面)后形成反射,反射振幅滿足Zoeppritz方程描述的AVO關系。鄭曉東[29]對Aki等[30]提出的AVO近似公式進行了適當變換,利用射線參數統一了不同相態介質界面上入射波與反射波、透射波之間的關系公式。水與空氣界面上的反射系數為
(2)
式中:p是射線參數;ρw、ρa是水和空氣的密度;a是空氣中聲波速度;θ是入射角。由于ρwc?ρaa,則R0絕對值約等于1。受風浪等因素影響,實際海面是不平整的,地震波在粗糙海面除了反射還會發生散射,散射波場只有一部分與反射波相干疊加,因此實際反射系數絕對值應小于1。
Jovanovich等[31]根據波浪譜的高斯分布及散射理論,推導出粗糙海面均方根波高幅值與海面反射系數之間的關系
(3)
式中:σ是波高的均方根幅值;ω是地震波角頻率;α是波的掠射角。由于掠射角與入射角互余,式(3)可寫成包含射線參數的形式
R(p,ω)=R0exp[-2(ωσp)2]
(4)
可見粗糙海面反射系數是地震波射線參數和角頻率的函數。
由上述分析可知,多個單震源的虛震源在檢波點的調諧組合子波可表示為
(5)

震源信號及其水面反射的地震波傳播至接收電纜,以壓力場變化的形式被檢波器記錄。為了壓制噪聲,每道由多個壓力檢波器以一定間隔平行連接,構成檢波器組合陣列(圖1),組合后的地震波是各單獨檢波器記錄地震波的簡單疊加。根據檢波器組合理論,地震道中記錄的信號s(t)可表示為
(6)
式中hk(t)表示組合陣列中的第k個檢波器記錄的信號。

圖1 電纜檢波器組合示意圖
直達波是指從震源出發沿介質直接傳播到接收點的地震波。如果傳播介質是常速的,其傳播路徑就是直線。在海洋地震勘探中,由于聲波在水層中的傳播速度近似常數,其直達波沿直線傳播。如圖2所示,震源能量在海面的反射波被稱為震源鬼波,可看作是從虛震源(S′,即震源點S關于水面的鏡像點)發出的直達波。圖中G為檢波點,ds為震源沉放深度,dg為檢波點沉放深度,x為炮檢距,θ為海面反射的入射角,r表示震源到檢波點的直達波傳播距離,r′為虛震源到檢波點的直達波傳播距離。從該圖可得如下關系式

圖2 直達波傳播路徑示意圖
(7)
(8)
則震源直達波的時距曲線方程為
(9)
虛震源直達波的時距曲線方程為
(10)
從式(9)、式(10)可見: 當震源深度與檢波點相同時,震源直達波時距曲線是直線; 否則,時距曲線為雙曲線。虛震源直達波(水面反射波)的時距曲線是雙曲線。射線參數可表示為
(11)
海洋地震記錄中的直達波h(t)是震源與虛震源直達波的疊合,即
h(t)=u(t)+u′(t)
(12)
實際地震采集中,震源和接收點都是由組合陣列組成。對于檢波器陣列中的每一個檢波器,它記錄的震源信號來自于震源陣列中的每一單個震源。設震源由m個互不相干的單震源組合而成,檢波點由l個檢波器并聯組合而成,每炮地震數據有n道記錄;用i表示震源陣列中單個震源序號、j表示接收道序號、k表示地震道組合中檢波器序號;sj(t)表示第j道的直達波記錄,hjk(t)表示第j道中第k個檢波器的記錄,ni(t)表示震源陣列中第i個單震源的震源信號。由式(1)、式(5)、式(12)推知,單檢波器記錄的直達波疊加信號為
(13)
代入式(6)可得
(14)
變換到頻率域
(15)
對于有n個地震記錄道的單炮記錄,直達波與震源信號子波之間的關系可記成矩陣方程
S=AX
(16)
其中
陣列震源的遠場是相對于近場而言,遠場距離的定義及計算公式[32-33]為
(17)
式中:f為激發子波頻率;λ為激發子波最小波長;D為震源陣列組合的空間尺寸。例如,對于組合尺寸等于30m的震源陣列,觀測數據中保留不高于250Hz震源信號的遠場距離約為150m。這一條件對于海洋多道反射地震勘探而言通常是滿足的,即地震采集的最小炮檢距大于震源陣列的遠場距離,因此可認為地震記錄接收到的震源信號子波是震源陣列的遠場子波。
在遠場條件下,震源陣列可被視為位于震源陣列中心的點震源,則震源信號相位譜與炮檢距不相關,式(16)中向量X也就退化為描述遠場子波的單變量,記為w,矩陣方程退化為
(18)
且有
(19)
此即為直達波計算震源遠場子波解析解公式。
實際地震數據處理中可利用疊加方法增強遠場子波解的抗干擾能力,即

(20)
此次合成了共計36道記錄的直達波數據作為輸入進行測試。使用有52ms延遲、主頻為75Hz的Ricker子波作為震源信號,接收道炮檢距范圍是225.0~662.5m,道間距為12.5m,每個接收道由16個檢波器組合構成,道內檢波器等距組合、間距為0.781m,水速為1540m/s,水面反射系數選用常數-1,震源深度為7m,電纜深度為15m。
分別測試了無干擾、加入10%隨機干擾和加入10%相干干擾(20Hz正弦波)等三種情形。從輸入及輸出結果(圖3)中可見無干擾情況下計算結果與輸入震源信號相一致。由于使用了疊加方法,加入隨機干擾后所得計算結果也具較高信噪比,加入相干干擾的輸入數據對解的影響也很小。
使用本文方法從南海M區二維深水反射實際地震數據中提取遠場子波,并與近場記錄擬合的遠場子波進行對比。部署該測線的目的是獲取地殼深部結構信息。現場采集選用6420in3的大容量震源,采用5排、間隔10m的槍陣構成氣槍陣列,每排槍(序列號依次為G1R1~G1R5)由5對相干槍組成,共計25對相干槍,每對相干槍含兩支單槍。相干槍可被看作是單震源,在其上方1m處布設近場檢波器。表1列出了槍陣組合參數,包括每個相干槍的容量、通道號及x、y坐標。坐標系統以槍陣中心為原點、以船首方向為x軸、船右方向為y軸。采用Sercel公司Sentinel型24位數字固體(接收)電纜,共有480個地震道,每道由8個檢波器組成,表2為檢波器組合參數。

圖3 合成數據測試

G1R5G1R4G1R3G1R2G1R1通道號x/my/m容量in3通道號x/my/m容量in3通道號x/my/m容量in3通道號x/my/m容量in3通道號x/my/m容量in31-6-205006-6-1076011-6050016-61076021-6205002-2-201807-2-1018012-2018017-21018022-22018031-2025081-102501310500181102502312025044-2012094-1080144012019410802442012056-2080106-10120156080206101202562080

表2 檢波器組合參數
選取其中一炮的計算結果進行分析和說明。如圖4所示,1~25通道的地震信號分別對應25個近場檢波器記錄的震源激發所產生的近場信息,其中用紅色線條顯示的第3、9、18號相干槍未振動、第7號狀態異常,其他21對相干槍工作正常。利用近場記錄擬合的遠場子波和利用本文方法提取的遠場子波分別用綠色和黑色線條同時顯示在35通道號位置,可見兩條曲線幾乎重合,所反映的氣泡振蕩周期完全相同,曲線正、負波峰峰值基本一致。二者的差異主要在高頻部分。分析可能的原因: 一方面與近場記錄的高頻干擾、擬合算法誤差有關; 另一方面,本文方法使用了疊加運算,當計算參數存在誤差時,算法存在一定的高截濾波作用。

圖4 近場記錄(1~25道)及利用其擬合的遠場子波(35道,綠色)、采用本文方法提取的遠場子波(35道,黑色)
為了測試從地震數據直達波提取震源信號遠場子波方法的有效性和實用性,選取2015年5月采集于中國南海北部陸坡的S區二維地震測線的部分單炮數據進行實際資料計算。
現場采集使用總容量6400in3的40槍組合陣列作為震源,槍陣由四排子陣構成,每排子陣陣列長度為15m,子陣間距為10m(圖5);槍陣沉放深度為8.5m,電纜沉放深度為16m; 采用Sercel公司Sentinel型24位數字固體電纜; 每個地震道由8個檢波器組成,檢波器靈敏度為19.7 μV/μbar(20°C),檢波器組合參數參見表2; 地震接收道炮檢距范圍是225.0~6212.5m,道間距為12.5m,數據記錄長度為8s,采樣間隔為2ms,480道接收。在原始數據中,未記錄近場信號。

圖5 槍陣組合結構圖
地震數據預處理首先是根據導航信息定義精確的觀測系統,明確地震數據炮點、檢波點的水平和垂直位置,計算炮檢距; 然后進行一些必要的去噪處理,去除數據中的野值干擾、交流電噪聲,消除直流漂移量,濾除低頻涌浪干擾等,但不做任何改變信號子波或改變振幅強度等方面的處理; 最后從地震數據中分離出直達波,該直達波數據長度從其初至時間開始計算,故應大于所求子波長度。基于直達波時距曲線,由式(9)和式(10)計算出不同炮檢距的初至時間。
為使測試更具針對性,規避直達波與反射波混合的情形,選取數據的海底深度約為1600m,海底反射波出現在2s以后(圖6)。由槍陣組合尺寸、地震數據Nyquist頻率及聲波速度(1500m/s),利用式(17)可算出震源信號遠場距離為150m,地震數據記錄位置的最小炮檢距滿足遠場接收條件。另外,從式(12)可見,震源直達波與虛震源直達波相位相反、二者時差隨炮檢距增大而減少,因此二者的疊加能量隨炮檢距的增大會快速減弱。慮及環境噪聲的存在,遠炮檢距數據信噪比也會變得很低,因此所選數據炮檢距不宜太大。選取近炮檢距36個單炮作為測試輸入數據,其炮檢距范圍約為225~663m(圖6矩形框)。圖7為其中10炮用于提取遠場子波的直達波數據,截取1.5s數據長度以確保求取子波具有1s時長。

圖6 經過預處理的單炮記錄左上角白色矩形框圈出測試計算用數據范圍
有了直達波數據后首先預設一個聲波速度,利用式(9)計算直達波的初至時間,調整聲波速度直至所計算的初至時間與地震數據的初至時間基本一致。用該方法可確定測試數據的聲波速度是1540m/s。從現場采集作業班報獲知,本次測試數據采集時恰遇一定風浪,因此采用平均反射系數,但其值通過迭代計算確定,細述如下。
首先,依據式(20)在地震數據處理軟件Omega系統上編制了在頻率域用直達波提取震源信號遠場子波的流程;預設R=-1,逐炮逐道計算;對于每一道數據,根據檢波器組合參數算出每個檢波器的位置,每一炮的子波是其36道數據計算的36個子波的疊加。然后,用計算出的子波預測直達波,逐步調整、加大反射系數R直至預測的直達波與實際地震數據的直達波基本一致。通過迭代,最后確定測試數據的平均反射系數R是-0.865。此時的子波計算結果就是本文方法最終結果,即利用直達波所提取的震源信號遠場子波(圖8)。還可用檢波器組的靈敏度將子波從電壓單位V轉換成通用的壓力單位mbar。

圖7 用于提取震源信號遠場子波的直達波數據

圖8 從圖7炮記錄直達波中提取的震源信號遠場子波
利用震源信號子波預測直達波有兩個主要目的:檢驗計算結果是否正確和確定平均反射系數。如前所述,當震源信號和聲波速度確定后,直達波形態只與反射系數有關,通過對比預測直達波與地震記錄即可判斷反射系數的值是否合適。
選取測試數據第241炮記錄中225.0、362.5、512.5、662.5m四個炮檢距做對比顯示(圖9),可見預測的直達波數據(紅色曲線)與實際地震記錄的直達波數據(藍色曲線)基本吻合。

圖9 預測的直達波(紅色)與實際地震記錄(藍色)
提取的震源信號子波常用于氣泡壓制、子波一致性處理等,以改善地震數據質量。
在Omega系統中使用Debubble模塊壓制地震資料中的氣泡,該模塊要求輸入每個單炮的震源信號及期望輸出子波。考慮到鬼波效應,輸入的震源信號是本文方法提取的子波加入了垂直方向的震源和檢波點鬼波,期望輸出子波是將作為輸入的震源信號子波壓制氣泡脈沖后的子波。
在子波中壓制氣泡脈沖的方法是利用脈沖的周期性使用預測反褶積壓縮子波,預測步長與脈沖周期基本一致,可通過子波自相關獲得。經過預測反褶積處理后的子波,氣泡脈沖被壓制,只保留了主脈沖。
由于本文提取震源信號子波的處理流程是在Omega系統中實現的,因此易將提取的計算結果直接應用于Debubble模塊以實現氣泡壓制功能(圖10)。Debubble模塊的另一項功能是震源信號的一致性處理。盡管使用了相同的槍陣參數,氣槍在激發過程中受其穩定性、環境因素等的影響,輸出信號很難保證完全一致,這種不一致性導致地震剖面存在非地質因素的橫向變化。Debubble模塊使用匹配濾波法使其輸出地震數據子波都與期望輸出子波一致,從而消除了震源信號不一致對地震數據質量造成的影響(圖11)。

圖10 氣泡壓制即子波壓縮過程

圖11 對測試數據做子波一致性處理前(a)、后(b)的近道剖面
震源信號是地震數據處理的重要參數。從信號與系統的觀點看,震源信號是主動源地震探測系統的輸入函數,它不依賴于探測目標、只與震源系統有關,因此其提取方法應盡量消除地層信息的影響。在此基礎上,本文提出一種地震子波提取的新思路,即從地震數據中分離出某類波場信息,且該波場傳播的介質信息已知、 并可建立地震子波與該波場之間的關系方程,然后通過解析方法計算得出地震子波。海洋地震探測由于其特殊的采集環境因素,地震數據中的直達波波場完全滿足上述條件,通過推導直達波與震源信號之間的關系式,形成了利用海洋地震資料直達波提取震源子波的方法。
(1)以氣泡振蕩理論、組合理論為基礎,利用地震運動學方法建立了直達波與震源信號之間的矩陣方程,并在遠場條件下將其簡化且導出頻率域利用直達波計算震源信號遠場子波的解析式。
(2)通過合成數據及實際地震數據測試,驗證了方法的正確性,同時也證明該方法具有較高計算精度和較強抗干擾能力。
(3)該方法適用于海洋地震資料的子波提取,無須近場記錄和其他設備,不依賴于褶積模型、不需地層反射信息,對地震子波的相位特性不做任何假設,且支持超長序列地震子波提取。
(4)震源信號子波在虛反射(鬼波)壓制、全波形反演、多次波預測等地震數據處理領域具有廣泛應用; 文中實例說明該方法提取的震源子波可方便地應用于其他后續處理環節。
(5)震源和檢波器的空間位置、海面反射系數、組合參數等的準確與否直接影響計算結果的精度;準確的直達波場分離技術是提高子波提取效果的另一關鍵因素。