被動源面波和體波成像在內(nèi)蒙古淺覆蓋區(qū)勘探應(yīng)用

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

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

0 引言

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

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

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

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

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

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

相關(guān)計算是被動源反射波勘探成像的計算基礎(chǔ)，其可通過互易性定理，基于聲波方程推導(dǎo)而來.

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

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

(1)

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

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

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

(2)

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

(3)

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

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

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

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

(4)

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

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

(5)

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

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

(6)

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

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

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

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

圖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

在上述理論和模型驗證的基礎(chǔ)上，在內(nèi)蒙古淺覆蓋區(qū)對被動源成像方法進行了數(shù)據(jù)試驗.

2 內(nèi)蒙古淺覆蓋區(qū)被動源勘探應(yīng)用

2.1 被動源數(shù)據(jù)的采集

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

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

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

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

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

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

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

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

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

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

(7)

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

(8)

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

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

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

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

圖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 利用高頻面波反演覆蓋層厚度

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

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

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

2.4 被動源反射波勘探成像

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

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

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

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

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

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

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

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

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

3 結(jié)論和建議

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