基于HHT的氣槍震源信號時頻分析與瞬時屬性研究

鄭露露　林建民　龐立臣　祝捍皓

摘要：

通過Hilbert-Huang變換對淺水環境下氣槍震源激發子波特征進行了研究，結果表明：氣槍震源首先激發高頻壓力脈沖，其頻帶相對較寬，優勢頻率為40～80 Hz，主頻隨時間降低；隨后高壓氣體在水下形成的氣泡振蕩激發低頻氣泡脈動，其頻帶相對較窄，優勢頻率為4～12 Hz，且水深較淺處的氣泡振蕩激發的氣泡脈動主頻相對較低。研究結果有助于提高對淺水環境下氣槍震源激發的復雜近場信號的認識。

關鍵詞：Hilbert-Huang變換；氣槍震源；高頻壓力脈沖；低頻氣泡脈動

中圖分類號：P3153文獻標識碼：A文章編號：1000-0666（2016）03-0466-07[HJ]

0引言

地震波是唯一能夠穿透地球內部的振動，迄今為止，關于地球內部的結構、組成、過程和狀態的研究等多來自地震波的信息（陳颙，朱日祥，2005）。震源是產生地震波的源頭，其所產生的地震波信號的品質將直接影響地震探測的效果。氣槍震源自1964年由美國Bolt公司的卡爾米思克發明以來，逐漸被發展成為具有綠色環保、可重復性好、可控性強、子波可測等卓越性能的水中震源，而被廣泛應用于海洋地震勘探（其中以地下油氣、礦藏等資源的局部尺度探測為主）（Calvert，2004）。近年來，氣槍震源也逐漸被應用于區域尺度的地震探測，出現了利用海洋激發氣槍、陸地接收信號的科學嘗試，如美國的LARSE（Los Angeles Region Seismic Experiment）計劃（Lutter et al， 1999；Fuis et al， 2003）、新西蘭的SIGHT（South Island Geophysical Transect）計劃（Okaya et al， 2002；Stern et al， 2002；Avendonk et al， 2004）以及國內的南海北部海陸聯測試驗（趙明輝等，2004；丘學林等，2007）。由于水中炸藥爆破存在可控性差、激發可重復性低以及作業費用高等缺點，氣槍震源亦開始逐漸成為陸地水域激發地震波的選擇。國內在陸上水庫進行的氣槍震源激發實驗表明：大容量氣槍震源是一種優良的水中震源，具有激發能量大、地震波轉換率高、傳播距離遠等特點，是陸上激發地震波的一種新方法（陳颙等，2007a，b；Chen et al， 2008；林建民等，2008a，b；楊微等，2013；劉自鳳等，2015；蘇金波等，2015；王彬等，2015）。

通過大量的理論研究與實驗觀測，人們已經很好地掌握了深海氣槍源的震源特性和以此組建氣槍陣列的技術（Mayn，Quay， 1971；Caldewll， Dragoset， 2000；Dragoset， 2000；羅桂純等，2007）。深海中氣槍作業水深多在30 m以上，如LARSE 計劃中氣槍的沉放水深大于100 m（Fuis et al，2003），周圍水體較大，震源相對比較獨立，激發信號受水底、水面反射等影響較小，激發過程相對比較簡單，這使得震源特性幾乎完全由氣槍本身決定。但水深較淺的陸地水域水體相對較小，使用氣槍震源必須綜合考慮氣槍震源與水域的相互作用、耦合等諸多問題，其激發地震波受激發條件和激發環境多種因素（水深、沉放深度、水底地形、邊界及沉積等）的影響，因此，其激發特性與規律和深海氣槍源存在很大差異。盡管水庫氣槍激發作為陸地震源的有效性已經得到驗證（陳颙等，2007b），但至今針對淺水環境下氣槍震源激發特性及激發地震波規律的研究甚少。林建民等（2010）曾重點分析了沉放深度、工作壓力等參數對水庫氣槍震源激發波形的影響，但尚未從氣槍震源信號本身特征出發加以分析。

氣槍震源信號是典型的非線性、非平穩信號，即統計量（相關函數、功率譜等）隨時間而變化。Huang等（1998）提出了一種分析非平穩非線性時間序列的新方法，即Hilbert-Huang變換（Hilbert-Huang Transform， 簡稱HHT）。該方法通過經驗模態分解（Empirical Mode Decomposition， 簡稱EMD），根據信號本身的時間尺度特征，將信號分解為有限個可以視為窄帶信號的本征模態函數（Intrinsic Mode Function， 簡稱IMF），并通過各IMF的Hilbert變換獲得具有明確定義和物理意義的瞬時頻率，建立了以本征模態函數為基本時域信號、以瞬時頻率作為基本量表征信號變化的新時頻分析方法。本文將利用HHT方法研究水庫氣槍震源激發信號時頻特征，并獲得瞬時頻率等瞬時地震屬性，提高對淺水環境下氣槍震源激發的復雜近場信號的認識。

11經驗模態分解

經驗模態分解的關鍵是將一個非線性非平穩時間序列分解為有限個IMF分量和一個趨勢項。由于每個時間序列都是由多個時間尺度的振蕩波構成，因此，該方法從本質上是將一個信號序列進行平穩化的過程，其處理過程是將信號中的這些固有的、內在的波動或趨勢項按照不同的時間尺度逐級分解開來，產生一系列具有不同時間尺度的IMF。

該方法可通過包絡線方法進行逐級篩選，實現時間序列的經驗模態分解，求出全部有限個IMF，主要步驟包括：

（1）對于原始序列X （t）的局部極大值和極小值，用三次樣條函數插值形成極大值包絡Xmax（t）和極小值包絡Xmin（t）；

（2）求取極大值包絡與極小值包絡的瞬時均值線

m1（t）=12[Xmax（t）+Xmin（t）]. （1）

（3）考察原始序列X（t）與瞬時均值m1（t）差值，即去掉低頻的新序列h1（t）=X（t）-m1（t），判定其是否滿足以下IMF要求，如果滿足，則將h1（t）作為第一個IMF分量C1；如果不滿足，則將h1（t）作為原始序列重復以上篩選過程，依次篩選，直至第k次篩選結果符合要求，即得到第一個IMF分量C1。

這樣得到的C1應是時間序列中特征時間尺度最小的高頻IMF分量。然后將C1從原序列 X（t）中分離出去，得到一個剩余序列r1=X（t） -C1，由于r1中可能含有比 C1的特征時間尺度大的IMF，因此又將r1作為新的序列實施上述（1）～（3）篩選過程，以求得第二個IMF分量C2。如此重復，求出一系列IMF，直至剩余序列rn（趨勢項）成為一個單調函數，從中再也不能分解出IMF分量時，整個分解過程結束。此時，原始序列可表示成為n個IMF和1個趨勢項之和，即

X（t）=∑nj=1Cj（t）+rn（t）.（2）

其中，IMF分量Cj需滿足以下條件：

（1）在整個資料集里，極值點的數目和穿零點的數目必須相等或者最多相差1個；

（2）由局部極大值所構成的包絡線和由局部極小值所構成的包絡線的平均值為零。

從上述EMD過程可以看出，該方法的優點在于：直觀、直接、基于數據本身特征以及自適應性強。

其中，P表示Cauchy主值。數學上，式（3）表示X（t）與1t的卷積，因此，Hilbert變換強調X（t）的局部性。在1t的基礎上，構造解析信號Z（t）為

Z（t）=X（t）+iY（t）=a（t）eiθ（t）. （4）

式中，a（t）為瞬時振幅，θ（t）為瞬時相位，對其求導可得瞬時頻率ω（t）為

ω（t）=dθ（t）dt.（5）

為使瞬時頻率ω（t）具有實際意義，信號X（t）必須為窄帶信號，其函數圖像關于零均值線是局部對稱的，并且穿零點和極值的數目是一樣的。若以上條件不滿足，所求出的瞬時頻率就可能出現異常值或負數。

通過EMD所分解得到的IMF為窄帶信號，符合作Hilbert變換求取瞬時頻率的要求，因此，原始信號可以用其各個IMF分量的瞬時頻率ωj（t）和瞬時振幅aj（t）近似地表示成

H（ω，t）=RE∑nj=1aj（t）eiωj（t）dt

.（6）

式（6）即為Hilbert譜，其中，RE表示求取實部運算。ωj（t）和 aj（t）都隨時間而變化，體現了非平穩序列的本質特征，突破了傳統傅立葉分析中定常振幅和固定頻率的限制。Hilbert譜將信號的振幅表示成為瞬時頻率和時間的函數，描述了其在整個頻率域上隨時間和頻率的變化規律。

2氣槍震源信號時頻特征及瞬時屬性

本文所研究信號為容量2 000 in3氣槍震源（Bolt 1500LL型）在水深18 m的陸上水庫環境下的激發信號，氣槍沉放深度10 m，工作壓力（氣槍中壓縮氣體的壓力）10 MPa，距離氣槍激發點水平距離50 m、深8 m位置處布設水聽器（采樣率200 Hz，頻帶50 Hz～40 kHz）所記錄的近場子波信號（林建民等，2010）。海洋石油勘探通常采用氣槍陣列作為震源，其子波測量點與氣槍源之間的距離通常在幾百米以上，主要是為了使陣列中各單槍激發信號能同時（在一個采樣時間間隔內）到達水聽器，而本文工作中采用單槍激發，且重點研究近場信號特征。

圖1為氣槍激發信號的原始波形及其EMD結果，圖1a是水聽器記錄的實際波形，圖1b是分解得到的7個本征模態函數（IMF1～IMF7），其中IMF7是剩余項。觀察各IMF分量時域波形，不同IMF分量包含了不同的時間特征尺度，分別以不同的分辨率顯示了氣槍信號特征，并且這種分辨率是自適應的。IMF1～IMF7分別是從高頻至低頻分量逐級分解獲得，其中，IMF1和IMF2主要對應氣槍震源信號的高頻壓力脈沖，平穩性相對較差；而IMF3和IMF4主要對應后續的低頻氣泡脈沖，平穩性比較好。

若直接利用原始波形進行Hilbert變換求瞬時頻率，則因其不完全滿足Hilbert變換的前提假設條件而產生異常值，不利于真正利用瞬時頻率掌握源信號頻率的時間變化，如圖3所示。由圖2可知，盡管IMF3的主要頻率范圍對應低頻氣泡脈沖，但是在高頻壓力脈沖位置仍有部分微小能量分布，導致在瞬時頻率曲線上出現高頻成分分布，這也說明不同IMF之間有頻率交叉，并不是上述的簡單頻域帶通濾波。

氣槍震源激發信號分高頻壓力脈沖和低頻氣泡脈沖，前者為高壓氣體瞬間釋放到水中激發壓力波，因具有較高頻率通常用于淺層石油勘探，后者為高壓氣體釋放到水中之后形成后續氣泡振蕩所產生，主頻相對較低，適合長炮檢距的深部結構探測。從圖2可以看出，IMF3在02 s附近瞬時頻率迅速由高頻轉為低頻，對應壓力脈沖向氣泡脈沖過渡，即高壓氣體瞬間釋放到水中激發壓力波之后形成的氣泡在02 s附近開始振蕩。但為掌握氣槍震源信號整體頻率隨時間的變化規律以及信號振幅或能量隨頻率和時間的變化分布情況，則需通過Hilbert變換的譜分析獲得Hilbert時頻譜，將振幅表示成為瞬時頻率和時間的函數。圖4給出了本文所分析氣槍信號的波形、功率譜密度分析及相應的Hilbert時頻譜，能清晰地刻畫出氣槍信號的時頻特征。由于Hilbert時頻譜包含了詳細的IMF分量頻率變化的信息，能準確反應信號的時頻特性，可以清晰地分辨出信號主頻的突變及所對應的時刻，因此，具有更高的時間—頻率域分辨率。

3不同時頻分析方法討論

目前，用于分析地震信號的最常用的、且較為成熟的方法是傅立葉變換。但由于傅立葉變換嚴格要求系統的線性、數據的周期性或平穩性（即傅立葉頻率代表著信號的周期性），在處理非平穩、非線性數據時，其實際應用受到很大的限制。如本文中氣槍信號為非平穩信號，其特點之一就是無周期性，若按傅立葉變換方法定義其頻率進行頻譜分析將缺乏物理基礎。因此，引入具有真實物理意義的瞬時頻率ω（t）對于研究氣槍信號瞬態特性具有重要的意義。

通常可利用Hilbert變換求取信號具有實際物理意義的瞬時屬性，但要求被分析信號是窄帶或平穩的，而且該方法對噪聲很敏感。而氣槍震源信號等實際地震信號既非平穩又含有噪聲，直接對其進行Hilbert變換求取瞬時頻率等屬性參數，將缺乏物理意義甚至失真（圖3）。HHT的根本創新就在于先通過EMD將信號進行平穩化，將信號中固有的、內在的波動成分按照不同的時間尺度逐級分解成窄帶信號分量IMF，使其滿足利用Hilbert變換求取具有實際物理意義瞬時屬性的要求，再分別計算每個分量對應的瞬時頻率等瞬時屬性。瞬時頻率作為時間的函數，能敏銳地識別出待分析信號中的多尺度嵌套結構（圖2）。

傅立葉變換是一種純頻域分析方法，是將信號從時域轉換到頻域的整體變換，無法提供時域信息。若要獲得某一時刻的頻率特性，則通常采用時頻分析方法，這也是針對非平穩信號較為常用的處理方法，比較典型的有短時傅立葉變換、小波變換等（李春峰，Linea， 2005；陳雨紅等，2006）。但這些方法基本上都是將信號序列在某種基底函數系上展開，然后分析展開系數及各分量的特征。如短時傅立葉變換實際上是加窗傅立葉變換，通常假設信號在窗函數的有效持續時間內是平穩的，但此條件通常無法滿足或只是近似滿足，因而它的時間和頻率分辨率都很低。圖5為所分析的氣槍震源信號通過短時傅立葉變換所獲得的時頻特征。小波分析雖然通過一系列可伸縮平移的小波函數實現了信號時頻局部化分析目的，但在本質上仍然是具有柔性時頻窗的加窗傅立葉變換（崔錦泰，2006），它也要求小波窗內的信號是短時平穩的，還沒有根本擺脫傅立葉變換的局限性。而且，小波基選擇的多樣性和使用過程中的不可變更性也影響了它在非平穩信號分析中的適應能力。同時，小波基的有限長會造成信號的能量泄漏，因而難以對信號作精確的時頻分析。

與傳統時頻分析方法最大的不同就是不采用固定的基底分解，而是基于信號本身的局部特征時間尺度進行分解獲得IMF分量，大大提高了其自適應性，有效克服了傅立葉變換中測不準原理的限制，在客觀性和分辨率方面具有明顯的優勢，更適用于非平穩、非線性信號處理。該方法能用分解得到的IMF分量及1個剩余項來揭示信號序列的振蕩結構特征和非平穩性，并最終由每個IMF分量經過Hilbert變換得到時域和頻域均具有較高分辨率的Hilbert譜（圖4c），具有明確的物理意義，且可以準確給出原信號序列及其IMF分量的主要振幅變化所對應的頻率和時間，即反映物理過程能量—頻率—時間的分布。

由于Hilbert時頻譜是通過計算IMF分量的瞬時頻率和瞬時幅度而直接得到，因此，其時間和頻率分辨率不同于其他時頻分析方法受Heisenberg不確定原理的限制，而是相互獨立的。相比于圖5所示基于短時傅里葉變換所得的時頻特征，圖4c所示的Hilbert時頻譜具有更好的局部時頻特性表現能力，從圖中可以觀察到氣槍信號能量隨時間和頻率的精細變化趨勢。從012 s左右開始，高頻壓力脈沖激發，由于信號脈沖較陡，其頻帶覆蓋范圍較寬，優勢頻段則主要集中在40～80 Hz。同時，得益于Hilbert譜的高時頻分辨率，還能觀察到壓力脈沖（圖4c中012～018 s之間原始波形中的4個脈沖峰）中后面的脈沖比前面的脈沖主頻低，即壓力脈沖主頻隨時間遞減，從80 Hz降低到40～60 Hz，然后直接過渡到后續的低頻氣泡脈沖。其原因可能是后面的脈沖是由壓力脈沖在水面及水底的反射或壓力脈沖本身的振蕩產生，其主頻在水介質中傳播及兩相界面上反射過程中下降或自身振蕩過程中下降。從02 s附近開始，釋放的高壓氣體在水下形成的氣泡開始振蕩激發低頻氣泡脈沖，持續時間達05 s以上，其優勢頻率集中在4～12 Hz的低頻范圍。從圖4b的分析結果可得，低頻氣泡脈沖在6 Hz和10 Hz附近有2個尖峰，但缺少了時域的分辨信息。而從圖4c中可以推測，10 Hz的尖峰能量可能主要來自于第一次氣泡振蕩，而6 Hz的尖峰能量主要來自于第二次氣泡振蕩。該實際計算結果與林建民等（2010）的分析結果中公式（1）相符合，由于氣泡脈沖激發過程對應著逐漸上浮的水下氣泡振蕩過程，第二次氣泡振蕩水深位置相對于第一次氣泡振蕩更淺，因此，公式（1）中P00更小，且R0更大，導致氣泡振蕩周期T更大，即對應氣泡脈沖主頻更低，所以，后續的氣泡脈沖主頻要低于前面的氣泡脈沖（圖4）。此外，在012 s壓力脈沖到達之前，存在微小振幅的振蕩波形，這可能由采集過程中濾波因素引起，它對后續的Hilbert譜計算可能引入一定的干擾。

本文對氣槍震源信號的處理分析結果顯示了HHT方法的有效性，但同時也注意到該方法在應用過程中對噪聲的敏感性以及邊界問題的影響，如圖4中t=0附近出現近微小能量的近100 Hz高頻能量，這個是由邊界效應引起的干擾。如何進一步消除信號中夾雜噪聲及邊界效應的影響，對于進一步提高HHT的應用效果至關重要，這是該方法需要進一步完善的地方，也是在應用過程中需要進一步研究的問題。

4結論

HHT是近年發展起來的處理非線性、非平穩信號的高效時頻分析方法，該方法通過EMD將信號分解為IMF窄帶分量，并計算其具有實際物理意義的瞬時頻率，掌握信號頻率的時間變化特征。在此基礎上，進一步通過Hilbert譜分析獲得在時域和頻域均具較高分辨率的時頻譜，更為清晰地表征信號的時頻分布特征。氣槍震源信號為典型的非平穩信號，若直接利用基于傅立葉變換的各種處理方法，所獲得的頻率信息物理意義不明確，時頻分析結果分辨率不高。本文利用HHT對水庫氣槍震源近場信號實際水聽器記錄數據進行了分析，相對于以往工作，在氣槍震源信號在不同時間尺度上振蕩特征、具真實物理意義瞬時頻率時間變化特征、時頻域高分辨率能量分布特征及原因等方面獲得了更多認識，有助于對淺水環境下氣槍震源激發的復雜近場信號的理解。

感謝中國地震局地球物理研究所水庫氣槍試驗及其他參與單位對本文工作的支持。

參考文獻：

陳颙， 朱日祥. 2005. 設立“地下明燈研究計劃”的建議[J]. 地球科學進展， 20（5）： 485-489.

陳颙， 李宜靖. 2007a. 地震波雷達研究展望： 用人工震源探測大陸地殼結構[J]. 中國科學技術大學學報， 37（8）： 813-819.

陳颙， 張先康， 丘學林等. 2007b. 陸地人工激發地震波的一種新方法[J]. 科學通報， 52（11）： 1317-1321.

陳雨紅， 楊長春， 曹齊放等. 2006. 幾種時頻分析方法比較[J]. 地球物理學進展， 21（4）： 1180-1185.

崔錦泰. 2006. 小波分析導論[M]. 西安： 西安交通大學出版社.

李春峰， Liner C. 2005. 基于小波多尺度分析的奇性指數： 一種新地震屬性[J]. 地球物理學報， 48（4）： 882-888.

林建民， 王寶善， 葛洪魁等. 2008a. 大容量氣槍震源特征及深部介質中傳播的震相分析[J]. 地球物理學報， 51（1）： 206-212.

林建民. 2008b. 基于人工震源的長偏移距地震信號檢測和探測研究[D]. 合肥： 中國科學技術大學.

林建民， 王寶善， 葛洪魁等. 2010. 大容量氣槍震源子波激發特性分析[J]. 地球物理學報， 53（2）： 342-349.

劉自鳳， 蘇有錦， 王寶善等. 2015. 賓川主動源地震波走時變化分析方法研究[J]. 地震研究，38（4）：591-597.

羅桂純， 葛洪魁， 王寶善等. 2007. 氣槍震源激發模式及應用[J]. 中國地震， 23（3）： 225-232.

丘學林， 陳颙， 朱日祥等. 2007. 大容量氣槍震源在海陸聯測中的應用： 南海北部試驗結果分析[J]. 科學通報， 52（1）： 1-7.

蘇金波， 王寶善， 王海濤等. 2015. 利用大容量氣槍震源資料研究北天山地區介質衰減特征[J]. 地震研究，38（4）： 598-605.

王彬， 吳國華， 蘇有錦等. 2015. 賓川地震信號發射臺的選址、建設及初步觀測結果[J]. 地震研究，38（1）： 1-6.

楊微， 王寶善， 葛洪魁等. 2013. 大容量氣槍震源主動探測技術系統及試驗研究[J]. 中國地震，29（4）： 339-410.

趙明輝， 丘學林， 夏戡原等. 2004. 南海東北部海陸聯測地震數據處理及初步結果[J]. 熱帶海洋學報， 23（1）： 58-63.