被動源面波和體波成像在內蒙古淺覆蓋區勘探應用

劉國峰， 劉語， 孟小紅*， 劉瀾波， 蘇維俊， 王永志， 張致付

1 中國地質大學(北京)地球物理與信息技術學院， 北京 100083 2 美國康涅狄格大學地質系， 美國康涅狄格州 06269 3 吉林大學地球探測科學與技術學院， 長春 130061

0 引言

隨著礦產勘探程度的不斷深化，礦集區“透明化”立體探測成為國家深地資源勘探領域的重要研究方向.其對地球物理探測的重要需求是在3000 m深度范圍內實現高精度、高分辨率成像.反射地震勘探具有深部高分辨探測能力，但是三維大規模應用于礦產勘探存在成本高、采集困難等問題.目前在礦集區立體探測中，一種典型方案是采用網狀分布的二維反射地震，建立勘探區框架結構模型，或以解釋的結構和物性成果為約束，輔助面積性的重、磁反演，完成礦集區三維立體結構或者物性填圖.該方案利用了反射地震的勘探深度和分辨率優勢，但當研究區結構橫向變化大，地質結構復雜時，有限的地震剖面無法控制區域構造變化.因此，在礦集區實施經濟、有效的三維反射地震勘探是目前立體探測的待突破瓶頸之一(呂慶田等，2019).

被動源地震是指不依賴主動源激發地震波，而是通過固定位置的地震儀臺站/檢波器等，在一段時間內連續接收天然噪聲信號進行探測研究的技術方法.它是地震學進行大尺度地殼和巖石圈結構成像，反演地球深部物性的主要手段(Bostock et al.，2001;Moschetti et al.,2007;余大新等，2016).若將臺陣更密集的擺放，可完成更精細尺度的探測應用，例如密集臺陣和城市環境中進行的面波勘探(Wapenaar，2006;Wapenaar et al.，2010；王爽等，2018；Wang et al.,2020).相對于上述應用，從被動源中提取體波記錄進行反射波地震勘探的應用，在近年成為了研究熱點.

Claerbout早期曾指出：“通過地表兩點放置的兩個檢波器接收的天然震動信號進行相關計算，可以構建以一點為檢波器，另外一點為震源的主動源波場記錄.”(Rickett and Claerbout，1999).同時提出該計算思路的還有Cole(Cole，1995)，不同之處在于，Cole是指提取面波，而Claerbout指的是提取體波.基于天然噪聲提取體波的核心環節為相關計算(Interferometry)(Schuster，2001；Schuster et al.,2004;Schuster and Zhou，2006)，其由超聲計算領域發展而來(Weaver and Lobkis，2001,2002)，主要內容為：相鄰的兩個檢波器信號進行互相關計算等于兩個檢波器間激發接收格林函數與隨機震動信號的褶積.

隨著節點檢波器的發展和普及(史子樂等，2019)，連續接收被動源地震記錄成本大幅度降低，被動源反射地震探勘從理論走到了實際應用.自2007年以來，陸續有試驗性應用文章發表，其普遍特點是在主動源地震采集時，進行重合測線的被動源采集，通過主動源和被動源的單炮及成像剖面對比，測試被動源反射波勘探在不同領域應用的可行性，這其中包括在油氣勘探領域的應用(Draganov et al.,2009,2013)、二氧化碳存儲四維動態監測的應用(Xu et al.,2012；Boullenger et al.,2015；Cheraghi et al.,2017)以及金屬礦勘探的應用等(Cheraghi et al.,2015；Olivier et al.,2015；Roots et al.,2017).同時，正演計算的開展也在理論與實踐應用之間搭建了論證紐帶，分析了震源分布、震源長度、震源數量等對被動源成像質量的影響(Thorbecke and Draganov，2010；朱恒等，2012；張盼等，2015).同時，更精細化的被動源勘探地震數據處理方法也不斷提出，以增強弱體波信號，進而提高反射波成像質量，這其中包括照明分析方法(Vidal et al.,2011)以及在原始道集數據進行的各種濾波處理(Girard and Shragge，2019a，b).

被動源地震勘探無需主動源激發，其采集成本不及主動源十分之一，是目前最具潛力的在活動礦區進行三維勘探的替代選擇.但同時，被動源地震勘探也面臨體波信號弱，現有可參考案例少等問題.基于此，我們在內蒙古集寧淺覆蓋區進行了被動源地震勘探試驗.本文推導了基于聲波方程的被動源反射信息提取理論.在實際數據驗證階段，介紹了數據的采集、頻率域面波和體波分離以及對擬炮集記錄的處理等關鍵環節.面波反演和體波成像結果表明，被動源勘探在內蒙古淺覆蓋區的應用具有廣泛的前景.

1 被動源地震勘探擬炮集記錄提取的理論基礎

相關計算是被動源反射波勘探成像的計算基礎，其可通過互易性定理，基于聲波方程推導而來.

點源激發下的頻率域聲波方程可表達如下：

=-jωρδ(x-xA)，

(1)

其中，G(x,xA,ω)指震源在x處，檢波點在xA,處的波傳播格林函數的頻率域表示，ρ(x)和v(x)分別是傳播介質的密度和速度，ω是角頻率，δ(x-xA)是脈沖震源.

根據互異性定理(Wapenaar,2006)，同屬接收點的A,B兩點間格林函數G(xA,xB,ω)為:

-(?iG(xA,x，ω)*)G(xB,x，ω))nid2x,

(2)

式中?D是任意封閉區間，ni=(n1,n2，n3)表示該區間的法向向量，R{}為取實部計算.公式(2)是A，B兩點間格林函數的準確表達，其中包含了非均勻介質中傳播的直達波、散射、反射以及多次反射等.但在實際應用中，往往對其進行簡化后應用，其中包括假設應用介質為均勻介質、密度和速度為常量、?D邊界處介質光滑等(Wapenaar and Fokkema,2006),而公式(2)也隨之簡化為：

(3)

上述簡化過程雖然會為波場帶來異常振幅，但公式(3)中波場相位沒有變化，因此，依然可以用來進行結構性成像計算.

被動源勘探中，震源為隨機噪聲，若其震源函數為N(x,ω)，則A，B兩點的接收該震動的記錄為：

u(xA,ω)=∮?DG(xA,x，ω)N(x,ω)d2x,

u(xB,ω)=∮?DG(xB,x，ω)N(x,ω)d2x,

(4)

若N(x,ω)為隨機信號，在時間和空間上均互不相關，則滿足：

〈N(x,ω)N(x′,ω)*〉=δ(x-x′)S(ω),

(5)

其中，S(ω)是隨機噪聲的功率譜.〈·〉代表空間的加權平均.

將公式(4)和(5)代入公式(3)中，可得：

(6)

式中，左端項是指接收點A，B之間的格林函數與噪聲震源的褶積.右端項是檢波器A，B記錄數據的相關計算，其中，〈*〉代表不同互相關結果的疊加.

被動源反射波勘探中，依照公式(6)，固定位置A,將A處的震動記錄與移動的B點震動記錄分別進行固定窗口的互相關和疊加計算，可獲得形如炮點位置在A處，接收點在變化的系列B處的炮集記錄.因為A處并無實際的主動源激發，因此該記錄稱為擬炮集記錄.按設計的觀測系統變化A的位置以及與其互相關計算的移動B的范圍，最終可獲得勘探區的系列擬單炮數據集.圖1表示了被動源體波成像計算的原理和步驟，其中圖1a是被動源記錄生成擬單炮記錄的原理示意圖，被動源在A處的記錄與其經界面反射后的B處的記錄褶積，可獲得A處激發，B處接收的反射記錄.圖1b示意了A，B兩點進行相干計算的過程.選定窗口長度，進行對應時間的相關和疊加，最終形成如圖1c的擬單炮記錄.將炮集記錄按主動源反射地震流程處理，可獲得反射成像剖面.

圖1 被動源記錄相干計算生成擬炮集記錄示意圖Fig.1 Schematic diagram of cross-coherent calculation of passive seismic records to generate virtual-shot records

上述應用隨機噪聲記錄生成擬炮集記錄計算的過程可以應用圖2的模型，采用正演計算進一步說明.模型含三個不同速度和密度的層位，其中8000個地下震源隨機分布.正演計算采用聲波方程的有限差分計算完成(Thorbecke and Draganov，2011).地表記錄隨機噪聲的檢波點間距為50 m，記錄長度為4 s，圖3a為地表記錄的隨機噪聲記錄，圖3b為震源在模型正中，所有檢波點參與計算生成的擬炮集記錄.若變換震源位置，可以獲得按特點觀測系統采集的炮集記錄，進行常規處理可獲得反射波成像結果.

圖2 含有三個層位的模型，其中8000個震源隨機分布Fig.2 A model with 3 layers, there are 8000 random sources in the model

圖3 正演獲得的隨機噪聲記錄以及生成的擬炮集記錄Fig.3 The random noise obtained by forward modeling and the generated virtual shot gather

在上述理論和模型驗證的基礎上，在內蒙古淺覆蓋區對被動源成像方法進行了數據試驗.

2 內蒙古淺覆蓋區被動源勘探應用

2.1 被動源數據的采集

本文試驗所用的被動源數據采集區域行政區劃隸屬于內蒙古集寧豐鎮市官屯堡鄉，如圖4所示.該地區富含鉬等多種金屬礦產，屬巖漿熱液型礦床.在前期地質調查項目的支持下，該區進行了主動源二維地震采集工作.該次被動源試驗選擇與前期該區一條二維地震測線重合，可對成像結果進行互相驗證、對比.

圖4 數據采集區域位置圖，圖中紅線為本次被動源采集的2 km測線Fig.4 Position of the passive seismic line. The red line is our 2 km length passive seismic line

被動源數據采集采用節點式單分量檢波器，型號為SmartSolo的IGU-16.該節點檢波器自然頻率為5 Hz，1 ms采樣下可連續采集25天.本次共擺放檢波器100個，檢波點間距是20 m.采集時設置采樣間隔為1 ms，采集時長為連續的10天.為了能夠避免地表風噪，將檢波器埋于地下，埋藏深度約20 cm.

該測區無明顯強活動震源.向北2.4 km為S24省道，其他方向交織村級公路.白天省道車流量相對較大，多有運煤卡車經過.采集時值秋季收割季節，有收割機在農田作業，附近有農民刨挖土豆.上述是地表可見的自然震源.附近礦場處在停工狀態，無礦區活動.

圖5是近村級公路500 m處檢波器10月5日和6日兩天的震動記錄，總體白天震動強度要大于晚上，主要以公路交通和人文活動產生的噪聲為主.

圖5 某路邊檢波器記錄的兩天的被動源記錄Fig.5 Two days of passive seismic records by a roadside receiver

2.2 基于頻率域信噪比計算的面波和體波分離

被動源體波成像的理論計算中，假設天然震源在地下半空間均勻分布，但實際應用中，地表附近震源(以公路交通為主)無論從強度還是數量上均占優，導致面波在被動源信號中的能量強，而相對而言，體波則屬于弱信號.將體波從面波中分離出來，有助于提高體波信號成像信噪比.在重建的擬單炮記錄中，因為體波信號頻率低，與面波在位置和頻帶范圍都存在重合，濾波分離效果不佳.若從原始數據出發，在互相關前對時窗內的數據進行面波和體波成分的甄別和分離，再進行相干和疊加，則更為有效.

分析公共交通產生的面波特點，可見其具有與體波信號不同的頻譜形態.圖6a,b分別是面波和體波占優數據段的信號的頻譜.面波占優時，頻譜呈梳狀分布，而體波占優時，頻譜形態正常.這種差別可通過頻率域的信噪比量化表示.面波占優時時窗數據的信噪比比體波占優時低，可基于此特點，對面波和體波數據進行甄選.本文采用張軍華等(2009)頻率域信噪比的計算公式，對參與互相關計算的時窗內信號信噪比進行計算，再根據選定的信噪比閾值，判斷相關疊加結果分配到體波數據中還是面波數據中.公式(7)是頻率域計算信噪比的公式，圖7是其參數在功率密度譜上的表示.

圖6 被動源信號不同頻率成分的功率密度譜特征Fig.6 Power density spectrum characteristics of different frequency components of passive seismic signal

(7)

根據不同工區數據的統計和分析，確定甄選體波和面波的信噪比閾值ε，甄選和分流數據，分別相干計算得面波和體波為主的擬炮集記錄，過程概括為公式(8)，

(8)

其中，Gk(r1,r2,t)是r1和r2兩個檢波器位置震動記錄的相關疊加結果，t是記錄時間，k表示波類型，當k為0時表示為面波，為1時表示為體波.n是疊加的總窗口數，ai(r1,r2,t)是某個時段的互相關結果.

圖7 頻率域計算信噪比參考圖Fig.7 Spectrum for signal-to-noise ratio calculation in frequency domain

被動源數據量大，生成擬炮集過程的計算強度也很大.以本文采集為例，100個檢波器，10天接收時間，1 ms采樣接收的記錄超過500 GB.公式(6)中的互相關在頻率域計算比時間域具有更高的計算效率，輔以圖形處理器GPU上的批量并行傅里葉變換，可實現更快的急速計算.本文介紹的基于頻率域信噪比的面波和體波分離，可插入到上述傅里葉變換流程中進行，無需額外計算.

應用上述方法對本文采集的被動源數據進行擬單炮生成計算，計算時，互相關計算時窗長度選擇為5 s，并通過上述的頻率域信噪比方法進行面波和體波的分離，其中,各頻率值為f1=0,fL=2,fH=15,f2=40,fc=60，信噪比閾值測試后選擇為0.76,小于0.76的歸于面波疊加，大于0.76的部分歸于體波疊加.圖8a是同一炮點計算獲得的以面波為主的單炮，圖8b是分流面波后的突出體波的記錄，雖然其中依然有面波出現，但是相比于圖8a而言，有明顯的反射波能量出現.

圖8 頻率域分離計算后的擬炮集記錄(a) 以面波為主的擬炮集記錄; (b) 分離面波后突出體波的擬炮集記錄，紅色箭頭為反射波信號.Fig.8 Virtual shot gathers after frequency domain separation calculation(a) Shot gather with surface wave as the main energy; (b) Shot gather with body wave as the main energy after surface wave separation， the red arrows point to reflection wave.

2.3 利用高頻面波反演覆蓋層厚度

被動源勘探處理過程能獲得的輔助信息遠少于主動源地震，因此，需要充分發掘不同信息成分的表征能力.在主動源地震勘探實踐中，面波主要用來進行淺層結構探測，被動源重建的擬單炮記錄面波，同樣可以進行近地表結構探測應用.而在本文的勘探區域，面波反演結果可用來進行覆蓋層厚度的估算.

面波反演的過程分為面波的識別，頻散曲線的提取和橫波速度反演三步.本文所用面波即為圖8a所示計算的面波記錄.圖9是單炮面波記錄和從中提取的頻散曲線，從中可以看出，低階面波頻散曲線聚焦性、連續性都很好，易于識別，且可見相對連續的高階面波.和主動源面波反演一樣,一旦從被動源面波中提取了頻散曲線,接下來就可對頻散曲線進行反演以獲取橫波速度結構.本文采用遺傳算法反演,遺傳算法作為一種非線性全局優化方法,不要求苛刻的初始速度和深度模型,得到的最佳擬合模型穩定、可靠(趙東等，1995).反演中，僅應用了低階面波，采用一維反演計算.對測線上所有單炮記錄進行上述反演，綜合成圖，可得到圖10所示的橫波速度剖面.該剖面中，約100 m深度處可見明顯的速度結構縱向變化界面，可判斷其為覆蓋層分界面，這與已知的地質符合良好(嚴昊偉等，2017).覆蓋層厚度的準確計算對區域結構認識，以及面積重、磁的反演等都可提供約束.

圖9 面波炮集記錄及其頻散圖Fig.9 Shot gather with surface wave and its frequency-velocity spectrum

2.4 被動源反射波勘探成像

面波的反演是被動源勘探的傳統應用領域，而基于被動源的反射波勘探則是目前的研究熱點.經面波分離后的擬炮集記錄雖然反射依然較弱，但可應用反射地震資料處理的多次疊加技術增強反射能量.在擬炮集計算過程中，設定炮間距和炮檢距范圍，則可形成測線反射勘探的觀測系統.表1是本文采集數據生成的擬單炮記錄的觀測系統參數.為了提高反射波的信噪比，每個檢波點位置都進行擬單炮記錄的計算，圖11是測線覆蓋次數圖，其中最大覆蓋次數達到100.

圖10 面波反演的橫波速度模型Fig.10 Inverted shear wave velocity model with the surface wave

圖11 擬炮集記錄覆蓋次數圖Fig.11 Fold number of the virtual shot gathers

表1 被動源生成的擬單炮記錄觀測系統主要參數Table 1 The main geometry parameters of the virtual shot gathers from passive seismic

擬炮集數據處理采用主動源反射地震數據的處理流程，主要流程包括野外靜校正，本文以高程校正為主，振幅垂向補償和橫向均衡，去噪，反褶積，速度分析及疊加處理.在去噪中，只進行了異常振幅消除，因為在共中心點道集中，同相軸的連續性差，因此沒進行剩余靜校正處理.考慮到速度對成像結果的影響，在與主動源對比時，只進行疊加剖面的對比.圖12是處理中速度分析和局部疊加結果圖，速度譜中有明顯的能量團，相對應的共中心點道集有相對應的反射同相軸，而局部疊加結果中，可見連續性的反射.

圖12 被動源反射信號速度分析和局部疊加Fig.12 Reflection velocity analysis and local stack for passive seismic data

經過上述處理流程處理，圖13a是獲得的被動源記錄疊加剖面，而圖13b是與之重合的主動源成像結果.被動源成像結果頻率較低，所以以主動源記錄上的明顯反射為標準進行成像結果的對比.總體可見二者明顯反射同相軸都具有較好的一致性.首先A處均無明顯反射，且反射的起始位置一致.B和C處的同相軸形態和位置均較一致.D處被動源和主動源斜率稍有差別，但橫向上的延展一致.E處被動源剖面上有與主動源對應的同相軸，但是斜率上稍有差別，被動源可能因為頻率低的緣故，可見同相軸橫向上的延展和連續.B和C相接處出現錯段，可能解釋為斷層F，這種錯段在被動源剖面上也可以清晰可見.

圖13 被動源(a)和主動源(b)反射波疊加剖面對比Fig.13 Stack sections of (a) passive seismic and (b) active seismic for reflections

綜合上述的分析對比，雖然該區無明顯地下主動源記錄，且分離的面波表明天然噪聲中以面波為主，但經過去面波的擬炮集反射波成像結果依然與主動源記錄有較好的一致性，這為被動源進一步應用提供了參考.但同時，被動源成像剖面的低頻特性一定程度上限制了其單獨應用，但在主動源二維框架的約束下，進行三維被動源勘探，聯合主動源結果對主要結構的解釋，或可是被動源地震在該區的一個應用方向.而對被動源數據進行拓頻處理，可更易于對比和解釋.

3 結論和建議

本文基于被動源地震方法，在內蒙古淺覆蓋區進行了勘探試驗.提出了一種頻率域信噪比計算的分離面波和體波信息方法，分離的面波炮集用于覆蓋層厚度反演.體波炮集經反射波常規地震資料處理實現多次覆蓋疊加后，獲得了與主動源記錄一致性較好的反射地震疊加剖面.該試驗表明被動源勘探方法在內蒙古淺覆蓋區具有較好的應用效果.應用被動源進行反射波勘探不需要主動源激發，其采集成本是主動源方法的十分之一，可低成本的完成三維地震勘探，結合主動源二維框架的約束，可對研究區實施更大范圍的三維勘探.但同時，被動源反射波勘探還存在諸如數據頻率低、信噪比低等問題，需要更多的細致的、具有針對性的處理方法進行精細處理，特別是將體波信息從強面波干擾中分離出來的方法技術，以推動被動源反射地震勘探的生產應用.