阮小敏，陳明春，劉振東，王志輝，徐志伍，陳 淼，張新港

1 中國地質科學院，北京100037

2 中國地質調查局中國地質科學院 地球深部探測中心，北京100037

3 中石化石油工程地球物理有限公司 南方分公司，成都610016

0 引 言

地震勘探是地下結構探測和研究的關鍵方法，其中主動源是人工激發(fā)震源，其激發(fā)時刻和位置是精確已知的，信噪比高，但施工復雜，勘探成本高.而另一類被動源探測方法，無需人工激發(fā)震源，而是通過固定位置的地震儀臺站或檢波器等，在一段時間內連續(xù)接收天然噪聲信號，然后利用地震干涉從背景噪聲記錄中提取波場傳播信號（Campillo and Paul, 2003; Shapiro et al., 2005; Wapenaar and Fokkema, 2006）.傳統(tǒng)上，被動源地震是開展大尺度地殼和巖石圈結構成像、反演地球深部物性的主要手段（Bensen et al., 2007）.近年來，被動源地震方法以其綠色環(huán)保、采集成本低、數(shù)據(jù)重構靈活性高和低頻信息豐富等優(yōu)點，在勘探地震領域受到越來越多的關注（Chamarczuk et al., 2022; Cheraghi et al., 2017）.

體波成像和面波成像是利用被動源地震勘探的基本方法.總的來說，由于面波傳播距離遠、衰減慢、能量大，利用噪聲干涉法提取面波更容易，因此能夠反演地下橫波速度結構的面波勘探方法獲得了更廣泛應用（Chen et al., 2022; Hu et al., 2022;Le Feuvre et al., 2015; Xia, 2014）.體波信號相對較弱，從背景噪聲中提取體波具有一定挑戰(zhàn)性，但是體波具有更高的頻率，考慮體波成分可以提高地下結構成像的分辨率（Nakata et al., 2015）.從被動源數(shù)據(jù)中提取體波記錄進行反射波地震勘探，由于其成本相對較低，且能克服工區(qū)無法放炮問題和易于實施等優(yōu)勢，近些年成為了研究熱點，并在地質環(huán)境復雜、散射較為嚴重的工區(qū)獲得了較好效果，比如金屬礦集區(qū)的勘探（Chamarczuk et al., 2022;Cheraghi et al., 2015; 劉國峰等, 2021; Olivier et al.,2015）、CO封存點四維動態(tài)監(jiān)測（Boullenger et al., 2015; Cheraghi et al., 2017）、 油 氣 勘 探（Draganov et al., 2009, 2013)和煤田勘探（Ning et al., 2021）等領域.

反射地震勘探中的高品質成像離不開高品質的原始炮集，不同于主動源激發(fā)的高信噪比炮集數(shù)據(jù)，基于被動源連續(xù)噪聲信號構建的虛擬炮集往往質量不高，在實際應用中面臨諸多技術挑戰(zhàn).比如：地下實際震源數(shù)量有限且分布不均勻，使得干涉法計算得到的虛擬炮記錄中出現(xiàn)虛假反射；地表檢波器接收到的實際背景噪聲場往往主要受控于面波能量，反射體波則屬于弱信號；大量檢波器長時間觀測的海量數(shù)據(jù)導致的計算和存儲瓶頸等.針對這些挑戰(zhàn)，近年來人們對被動源反射波勘探成像開展方法研究并取得一些進展，主要集中在不同干涉算法如何更好地重構虛擬炮集，從面波主導的背景噪聲記錄中有效甄別和提取反射信號，以及被動源原始數(shù)據(jù)預處理、虛擬炮集多次波壓制和直接偏移成像等方面.最后通過幾個實例了解被動源反射地震勘探的應用現(xiàn)狀，并對將來的研究前景進行簡單展望.

1 地震干涉法重構虛擬炮集

地震干涉（seismic interferometry）是指通過對原始地震信號進行干涉合成新記錄，新記錄保留了原始信號中相關性高的目標響應、并且能反映原始信號中所不具有的一些重要特征，比如地震干涉被廣泛用在背景噪聲面波成像、震源定位、多次波成像等諸多方面（Shapiro et al., 2005; Snieder and Vrijlandt, 2005; Wapenaar et al., 2004）.通過對被動源噪聲數(shù)據(jù)進行地震干涉計算，可得到類似主動源激發(fā)的炮集記錄，稱之為被動源的虛擬炮集或擬炮集，而虛擬炮集的質量直接決定著最終反射剖面的成像效果，因此地震干涉法對被動源反射地震勘探而言十分關鍵.

1.1 地震干涉法的提出與發(fā)展

地震干涉技術的基本原理是通過對兩個檢波點接收到的地震記錄做一定的數(shù)學運算（互相關、反褶積、互相干等），提取這兩個檢波點間的格林函數(shù)，得到其中一個檢波點為新的虛擬震源，另一個檢波點接收的模擬地震記錄（圖1）.該方法由Claerbout最早提出（Claerbout,1968），他闡述了通過地下埋藏震源產生的地震記錄的自相關來求取地表格林函數(shù)的過程.Cole（1995）則嘗試通過計算實測被動源數(shù)據(jù)來合成反射記錄.Rickett和Claerbout（1999）將相關法推廣到多維模型，并且稱之為日光成像方法.Schuster在2001年的EAGE會議上，正式將該方法命名為地震干涉技術，并對被動源地震數(shù)據(jù)進行虛擬炮記錄合成、偏移成像等，且取得了較好的效果（Schuster, 2001; Schuster et al., 2004）.Wapenaar等（2004）、Wapenaar和 Fokkema（2006）用積分理論和互易定理進行了嚴格的理論推導，證明了地震干涉技術在不同的非衰減介質和不同的震源條件下均成立.

圖1 地震干涉示意圖（修改自Schuster, 2009）Fig.1 Seismic interferometry schematic diagram (modified from Schuster, 2009)

圖2 封閉區(qū)域模型的地震干涉（Wapenaar and Fokkema,2006）.格 林函數(shù)(xA,xB,ω)可由觀測點 x A、xB和沿著邊界? D 的 震源 x積分的互相關獲得Fig.2 Seismic interferometry model enclosed by boundary(Wapenaar and Fokkema, 2006) Green's function(xA,xB,ω)can be obtained by cross-correlating the observations at x AandxB and integrating along source coordinate x at?D

實際應用時通常對式（1）進行近似：（1）?

D

邊界之外是均勻介質；（2）遠場近似；（3）圍繞的介質參數(shù)均勻變化（Wapenaar and Fokkema,2006）.經(jīng)過近似后，所得的重構兩個檢波點之間的格林函數(shù)公式為：

當?

D

是一個球面，并且半徑非常大時，公式（2）是比較精確的，不過當?

D

是有限的時候，這些近似就會造成振幅上的誤差（Ramírez and Weglein, 2009），但是并不影響相位，所以公式（2）在誤差可接受的范圍內仍適用于地震干涉方程（Thorbecke and Draganov, 2011）.該方程也由Derode等（2003）通過物理推理的方式得到.

地震干涉法的提出深化了對地震波傳播規(guī)律的認識，很多方法都借助于地震干涉技術而提出或完善，大大豐富了地震學的理論與應用研究.一是基于背景噪聲的研究.Shapiro等（2005）通過互相關方法提取出美國加州地震臺站的瑞利波群速度頻散曲線，發(fā)現(xiàn)反演得到的速度結構模型與研究區(qū)的構造特性對應關系較好，開啟了背景噪聲無源成像的應用時代.隨后，利用背景噪聲成像研究巖石圈速度結構呈爆發(fā)式增長（Bensen et al., 2007; Ritzwoller, 2009）.二是基于大地震尾波干涉研究.Snieder（2004, 2006）在總結前人工作基礎上提出了大地震尾波干涉理論，進一步完善了地震干涉法理論.他們通過對地震臺站記錄的大地震尾波進行互相關，用來推斷介質參數(shù)隨時間的變化情況，即可實現(xiàn)面向目標介質的連續(xù)動態(tài)監(jiān)測.尾波干涉法在震源定位（Snieder and Vrijlandt, 2005）、地震活動（Peng and Ben-Zion, 2006）、火山監(jiān)測（Brenguier et al., 2008）、冰川活動（Jonsdottir et al.,2013）、地震波速動態(tài)變化（Wegler and Sens-Sch?nfelder, 2007）等研究領域中都有著廣泛的應用.

近些年人們又在常規(guī)的地震干涉法基礎上發(fā)展出了地震超越干涉法（Beyond Seismic Interferometry），該方法突破了傳統(tǒng)地震干涉法對于虛擬震源必須置于檢波器位置的限制，可實現(xiàn)虛擬震源或檢波器在地下介質的任意位置分布.Broggini和Snieder（2012）首先將地震超越干涉法的Marchenko成像引入地球物理領域，其本質是通過虛源點與地表之間的直達波記錄和地表反射響應反演得到上下行格林函數(shù).Broggini等（2012）、Wapenaar等（2013）將該方法擴展到二維和三維，能夠在無需了解上覆介質構造的情況下，基于一個宏觀的非精確模型就可對目標區(qū)域的介質進行精細成像.Marchenko成像方法在多次波壓制（Meles et al., 2015）、數(shù)據(jù)重構（靳中原等，2017）等方面有諸多應用，該方法也可被擴展到被動源數(shù)據(jù)，直接面向目標成像，避免檢索炮集中的層間多次波影響（Jin, 2020）.

1.2 地震干涉法在勘探中的應用

地震干涉法的提出不僅極大促進了天然地震學的研究，對勘探地震領域也有重大理論和實際意義.比如在主動源地震勘探中，復雜的近地表如風化帶、推覆體、小斷塊等通常制約著地下構造的精確成像，常規(guī)的勘探技術手段已經(jīng)很難滿足實際生產的需要.而地震干涉法可以通過逆時重構將井下檢波器變?yōu)樘摂M震源從而避免復雜近地表的影響，間接地達到面向目標成像的目的，提高了地震資料的信噪比和分辨率.由于這一過程中，不需要震源或地震子波的相關信息，該方法又稱作“虛源法”（圖3）.2004年，Calvert首次提出該方法（Calvert et al., 2004），隨后被陸續(xù)用于實際的VSP剖面（Bakulin and Calvert, 2006）、OBC數(shù)據(jù)（Gaiser and Vasconcelos, 2010）和 OBS數(shù)據(jù)（Carrière and Gerstoft, 2013）等勘探實踐.

圖3 虛源法試驗（曹輝等，2012）.（a）地表激發(fā)、井中接收；（b）虛源構建后的觀測系統(tǒng)，消除了近地表復雜條件對波場的透射影響Fig.3 Virtual source method experiments (Cao et al., 2012).(a) The sources are fired on the surface, and the receivers are located in the borehole; (b) Acquisition geometry after virtual source construction, removing the effects of the complicated near surface

對于因為環(huán)保、成本等因素限制而無法實施主動源地震激發(fā)的勘探區(qū)域，被動源方法可以較好地填補這些區(qū)域地下結構成像的空白.被動源勘探采集的背景噪聲信號粗看或許雜亂無章，但地震干涉法能夠變無序為有序，即能夠從背景噪聲中重構有效的地震波場并使之成像，進而推斷地下地質構造或監(jiān)測儲層變化等.自從 Wapenaar 等（2004）、Wapenaar和Fokkema（2006）對被動源地震干涉法在3D非均勻介質下進行嚴格的理論推導以來，被動源地震干涉技術不斷取得進步，目前利用被動源環(huán)境噪聲的面波勘探應用比較廣泛，比如在城市環(huán)境中開展的工程勘察及各種無損檢測（Chen et al., 2022; Cheng et al., 2015; Hu et al., 2022; 汪利民等, 2022; Xi et al., 2021）.但面波勘探方法分辨率相對較低，無法滿足更高精度的勘探需求，因此被動源體波勘探逐漸興起，這幾年基于地震干涉技術的反射波提取方向成為了地震干涉領域新的熱點.因為基于提取的被動源反射體波信息，不僅能得到速度結構，而且能獲得高頻率的反射剖面（Chamarczuk et al., 2022; Cheraghi et al., 2015,2017; 劉國峰等, 2021; Olivier et al., 2015）.

1.3 被動源虛擬炮集重構

基于地震干涉法對被動源噪聲數(shù)據(jù)實際計算虛擬炮集的時候，為了提高信噪比，通常采用分段時窗互相關＋疊加來處理記錄到的連續(xù)噪聲記錄.比如對于記錄長度為

T

的

N

個臺站數(shù)據(jù)，固定一檢波點不動作為虛擬炮點，與剩余

N

-1個檢波點在選取的小時窗

t

內作互相關，共進行

T

/

t

次互相關計算，最后將所有時窗內的互相關結果進行疊加就得到一個完整的虛擬炮集，再依次滾動計算生成其他位置的虛擬炮集記錄（方捷等，2022）.計算可知，重構所有虛擬炮集需要求解

T

×

N

×(

N

-1)/2的互相關函數(shù)，對于長時間觀測（一周以上）的大量檢波器（上百甚至上千臺）的被動源數(shù)據(jù)，總的計算和存儲量都非常巨大.完成所有虛擬炮集重構之后，再借鑒常規(guī)反射地震方法進行處理，即可獲得被動源反射成像剖面.由此可見，基于地震干涉法高效重構高質量的虛擬炮集是被動源反射波成像的關鍵.

1.3.1 多種干涉算法重構虛擬炮集

對不同檢波器記錄的連續(xù)波形數(shù)據(jù)進行多種干涉運算，比如經(jīng)典的互相關（Schuster et al., 2004）、反褶積（Snieder, 2006）和互相干（Prieto et al.,2009）三種方法，都可以提取這兩個檢波點間的格林函數(shù)，進而獲得地震虛擬炮集記錄.

假設在檢波器

r

點處的頻率域波動波場

u

(

r

,

s

,ω)，是非擾動格林函數(shù)

G

(

r

,

s

,ω)和擴散波格林函數(shù)

G

(

r

,

s

,ω)分別與

s

處的震源函數(shù)

W

(

s

,ω)褶積求和得到的：

式子（3）中的

W

(

s

,ω)是一復雜的頻率函數(shù)，隨震源函數(shù)

s

的變化而變化； ω為頻率；

G

為背景波場；

G

為擴散波場.

u

、

G

、

G

滿足聲波、彈性波、可衰減的波，或者是包含高頻擴散的非均勻波.如果忽略頻率 ω對震源函數(shù)

W

的影響，令

G

=

G

+

G

，則（3）式可簡化為：

因此，頻率域中A、B兩點的波場互相關可表示為（Vasconcelos and Snieder, 2008）：

而反褶積干涉成像消除震源特性的影響，其在頻率域可表示為（Snieder, 2006)：

同樣，互相干地震干涉也消除了震源項的影響，相當于空間自相關，頻率域的互相干干涉可表示為（Prieto et al., 2009)：

從公式（5）可知，互相關方法需要知道震源功率譜，但該方法穩(wěn)定性高.而公式（6）和（7）中已經(jīng)證明，反褶積和互相干方法不需要知道震源功率譜，且消除了震源項的影響，提高了分辨率.

很多學者曾對三種經(jīng)典干涉算法的效果進行研究對比，并分析各方法的優(yōu)缺點.Liu等（2020）用模型數(shù)據(jù)進行虛擬炮集重構試驗，證實反褶積和互相干方法獲得的虛擬炮集記錄信號更接近于尖脈沖信號，分辨率和信噪比也更高（圖4）.Snieder等（2009）認為反褶積方法提取的反射波響應可能會在參考檢波點引入虛假反射.Cao和Askari（2019）比較不同方法效果后認為互相干方法是開展二氧化碳儲存時移監(jiān)測的更好選擇.方捷等（2022）對內蒙地區(qū)實際數(shù)據(jù)進行分析對比后，也建議采用互相干方法，因為較之反褶積方法，它包含的波形信息更加豐富，較之互相關方法，其面波同相軸延伸更短，分布范圍較小，中間無強噪聲分布（圖5）.

圖4 不同地震干涉方法計算的模擬數(shù)據(jù)虛擬炮集（修改自Liu et al., 2020）.（a）互相關；（b）反褶積；（c）互相干Fig.4 Virtual shot gathers using noise-free data as input; generated by (a) cross-correlation, (b) deconvolution, and (c) cross-coherence (modified from Liu et al., 2020)

圖5 不同地震干涉方法計算的實際數(shù)據(jù)虛擬炮集（方捷等，2022）.（a）互相關；（b）反褶積；（c）互相干Fig.5 Virtual shot gathers using low S/N data as input; generated by (a) cross-correlation, (b) deconvolution, and (c) cross-coherence (Fang et al., 2022)

在二維或三維介質中，公式（9）是病態(tài)方程，可通過最小二乘反演來求解：

其中*表示轉置和復共軛，

W

是對角加權矩陣，

I

是單位矩陣， ∈是保證最小二乘反演穩(wěn)定性的一個參數(shù).這個矩陣反演的過程等同于在時間域沿著上邊界?

D

自由表面多次波響應與下傳波場的反褶積.該方法也被稱為多維反褶積.不規(guī)則分布震源（圖6a）的地震干涉數(shù)值計算示例表明，多維反褶積方法的結果（圖6d）較互相關計算結果（圖6c）改善明顯.但值得注意的是，多維反褶積方法需要同時對所有檢波器的記錄進行求逆運算，相較于傳統(tǒng)干涉方法，在穩(wěn)定性上略有不足（Minato et al., 2009）

圖6 不規(guī)則分布震源的地震干涉數(shù)值計算示例（修改自Wapenaar et al., 2008）.（a）橫向均勻的水平層狀模型，不規(guī)則分布震源依次發(fā)出瞬態(tài)信號；（b）單個地下震源的地震記錄響應；（c）互相關重構炮集（黑色點折線）與源在x=0處的直接模擬炮集（紅線）對比；（d）多維反褶積炮集（黑色點折線）Fig.6 Numerical simulation of seismic interferometry with irregular sources (modified from Wapenaar et al., 2008).(a) Configuration with a horizontally layered target below a homogeneous overburden and a free surface.The irregularly distributed sources below the target emit transient signals sequentially.(b) Response of one of the sources.(c) Interferometry from cross-correlation (red traces) compared with the directly modeled response of a source at x=0 (black dashed traces).(d) Interferometry from multidimensional deconvolution

1.3.2 影響虛擬炮集重構的其它因素

除了不同干涉方法對重構虛擬炮集的影響之外，噪聲震源也是影響干涉重構效果的關鍵因素，很多學者探究了如震源類型、震源分布、震源長度、震源數(shù)量等對重建虛擬炮集記錄和最終成像的影響（Thorbecke and Draganov, 2011; 王德利等, 2012;張盼等, 2015; 朱恒等, 2012）.

首先，地下震源數(shù)量對干涉成像影響較大，較多的地下震源數(shù)量對虛擬炮集成像更有利.圖7是Thorbecke和Draganov（2011）進行的一個被動源單炮模擬數(shù)值試驗，對于震源在地下空間500～4 000 m隨機分布的地下五層模型（圖7a），將被動源重構虛擬炮集與模擬的主動源激發(fā)單炮進行比較（圖7b），其中8 000個隨機噪聲震源重構的單炮反射較為清楚（圖7c），1 000個震源重構單炮成像次之（圖7d），50個震源重構單炮效果最差（圖7f），且有很多假象.

其次，非相干的震源一般來說有助于地震干涉反射信號的重構，因此震源性質和位置分布也影響著虛擬炮集的成像.圖7d、7f、7g的震源在數(shù)量上雖然都相同，不過震源之間的非相干性依次降低，其中圖7d采用的是隨機噪聲震源，非相干性最強.圖7f和7g的所有地下震源均采用相同的Ricker子波，但圖7f的激發(fā)時間和震源位置都是隨機的，仍能部分保證震源之間的非相干性，因此反射信號還可以較好地被重構；而圖7g的震源均分布在2 700 m深度且同時激發(fā)，非相干性最低，因此反射數(shù)據(jù)檢索最差，幾乎都是假象.可見，非相干的地下震源對于構建高質量的虛擬炮集反射波十分重要.

圖7 不同震源參數(shù)模擬的地震干涉結果（修改自Thorbecke and Draganov, 2011）.（a）地下震源隨機分布的10×4 km數(shù)值模型；（b）地表中點激發(fā)的主動源單炮；（c-e）分別用8 000、1 000、50個隨機位置地下噪聲震源激發(fā)構建的被動源虛擬炮集；（f）1 000個隨機位置雷克子波震源激發(fā)的重構虛擬炮集；（g） 1 000個固定2 700 m位置深度的雷克子波震源同時激發(fā)的重構虛擬炮集Fig.7 SI results from modeling by different sources (modified from Thorbecke and Draganov, 2011).(a) The 10×4 km model with sources positioned at random locations, visible as black dots.(b) A directly modeled reference result for an actual source at x=0 m.SI results for noise signature sources for varying numbers of (c) 8 000, (d) 1 000, and (e) 50.(f) SI results from 1 000 sources using the Ricker wavelet in random locations.(g) SI results by 1 000 sources using the Ricker wavelet at z=2 700 m;the sources are started simultaneously

基于被動源單炮數(shù)值模擬試驗可知，噪聲源分布和噪聲源類型對于虛擬炮集的反射波提取起著關鍵作用，而數(shù)值試驗本身也可為實際被動源數(shù)據(jù)的采集部署提供一定依據(jù).比如較多的地下震源數(shù)量對虛擬炮集成像更有利，但是實際野外采集數(shù)據(jù)時，震源的數(shù)量人為并不可控，在構造穩(wěn)定的區(qū)域往往缺少足夠多的地下噪聲源，而地下構造相對活躍的區(qū)域（比如地下正在開采的礦區(qū)等）一般來說更適于被動源體波反射成像.此外，適當增加記錄時長通常能增加震源數(shù)量，可有效改善最后重構單炮記錄的效果.但是記錄時長和重構單炮記錄的效果并不是呈簡單線性關系，更多的震源可能意味著來自地面的震源占據(jù)大多數(shù)，從而使被動源信號中面波占據(jù)主導，即更多的噪聲源分布在反射波干涉的非穩(wěn)相區(qū)，使本來能量就很弱的體波被淹沒，增加后期面波和體波分離的難度.因此被動源勘探的野外采集部署時應注意采集時間和最終成像效果兩者的平衡，目前被動源反射波勘探采集的記錄時長一般為幾天（劉國峰等，2021）到幾十天（Chamarczuk et al., 2022）.

2 被動源數(shù)據(jù)中弱體波反射信號的提取

由于實際環(huán)境噪聲場往往主要受控于面波能量，地表附近震源往往無論從數(shù)量還是強度上均占主導，導致面波在被動源信號中的能量強（特別是檢波器通常在地表布置），相對而言，體波反射則屬于弱信號.而想用地震干涉法重建體波，噪聲源必須分布在適當?shù)钠椒€(wěn)區(qū)域，其中大部分需要埋藏在地球內部（Shapiro, 2019）.Vidal等（2014）的數(shù)值模擬結果直觀地反映了這一點（圖8）.如果對原始噪聲記錄不加甄別和提取，直接對所有記錄（包括面波能量占優(yōu)和體波能量占優(yōu)）的噪聲窗口進行相關求和的虛擬炮集將以面波為主導（圖8a），體波反射能量則被淹沒其中，而對體波能量進行有效提取并單獨重構的虛擬炮集反射則十分清晰（圖8b），除了頻率較低之外，幾乎可與主動源激發(fā)的單炮媲美（圖8c）.

圖8 體波提取前后的地震干涉模擬結果（修改自Vidal et al., 2014）.（a）利用所有記錄重構的虛擬炮集；（b）提取體波后的重構虛擬炮集；（c）主動源激發(fā)單炮Fig.8 SI result for selection of body wave retrieval (modified from Vidal et al., 2014).(a) Retrieved virtual common-source panel with all recorded data.(b) Retrieved virtual common-source panel after the selection of body wave retrieval.(c) Directly modeled reflection response

因此如何從大量背景噪聲連續(xù)記錄中更好地甄別和提取弱體波信號，就成為被動源反射波勘探成像的核心問題.如果對已重建完成的虛擬單炮記錄再進行處理，因為被動源體波與面波在空間位置和頻帶范圍都存在重合，F(xiàn)K等傳統(tǒng)濾波方法效果往往并不理想.若是從原始數(shù)據(jù)出發(fā)，在互相關之前就對時窗內的數(shù)據(jù)先開展面波和體波成分的甄別和分離，再進行相干和疊加，則更有助于提高體波反射信號的信噪比（Chamarczuk et al., 2022; Cheraghi et al., 2015; 劉國峰等, 2021）.目前的研究主要基于被動源數(shù)據(jù)面波和體波的多方面差異特征，包括信噪比、速度以及方位角等差異來識別和提取弱體波反射信號.

2.1 基于信噪比差異提取被動源體波

被動源野外勘探實際采集的時候，不同時段的噪聲主要來源往往有所差異，如白天時段人文交通等面波噪聲更為突出集中.人們通過分析野外實際被動源采集數(shù)據(jù)時發(fā)現(xiàn)，面波占優(yōu)時段與體波占優(yōu)時段的信號頻譜特征往往差異較大.例如圖9a、9b分別是面波和體波占優(yōu)數(shù)據(jù)段的信號的頻譜.面波占優(yōu)時，頻譜呈梳狀分布，而體波占優(yōu)時，頻譜形態(tài)正常（劉國峰等，2021）.這種差別可通過頻率域的信噪比量化表示，見式（11）：

圖9 被動源信號不同頻率成分的功率密度譜特征（劉國峰等，2021）.（a）面波能量占優(yōu)；（b）體波能量占優(yōu)Fig.9 Power density spectrum characteristics for different frequency components of the passive seismic signal (modified from Liu et al., 2021).(a) Energy dominated by surface wave; (b) Energy dominated by body wave

其中，SNR為所述噪聲信號的信噪比，|

X

(

f

)|為頻率值

f

對應的功率值，頻率值

f

＜

f

＜

f

＜

f

＜

f

＜

f

，其中

f

和

f

為功率譜密度計算的最小和最大頻率值，

f

和

f

為有效信號的頻帶范圍，

f

和

f

則是包括有效信號和隨機噪聲的頻帶范圍（圖10）.

圖10 頻率域計算信噪比參考圖（張軍華等，2009）Fig.10 Spectrum for the signal-to-noise ratio calculation in the frequency domain (Zhang et al., 2009)

劉國峰等（2021）等利用該方法對內蒙古淺覆蓋區(qū)礦區(qū)的數(shù)據(jù)進行了處理，首先計算獲得不同時窗內的噪聲信號的信噪比，并基于信噪比差異將噪聲信號分選為面波信號和體波信號.如果計算確定某段時窗內屬于面波能量為主的信噪比范圍，就可以在合成虛擬炮集之前去除這段時窗內的記錄，最終可得到體波能量占比更高的虛擬炮集.圖11a是未進行面波體波分離處理的單炮，圖11b是分離處理后的記錄，相比于圖11a有較明顯反射波能量出現(xiàn).

圖11 頻率域分離計算后的虛擬炮集記錄（劉國峰等，2021）.（a）以面波為主的虛擬炮集記錄；（b）分離面波后突出體波的虛擬炮集記錄，紅色箭頭為反射波信號Fig.11 Virtual shot gathers after frequency domain separation calculation (Liu et al., 2021).(a) Shot gather with surface wave as the main energy; (b) Shot gather with body wave as the main energy after surface wave separation.Red arrows indicate the reflection wave

2.2 基于速度差異提取被動源體波

除了信噪比特征差異外，面波速度一般要明顯低于體波，人們也可用速度差異來識別和提取體波.例如Vidal等（2014）采用傾斜疊加照明分析的方法來區(qū)分二維測線上的低速面波和高速體波，該研究設計了復雜的地下彈性模型，并基于地表震源和地下震源數(shù)據(jù)分別構建虛擬炮集，再用不同的射線參數(shù)對虛擬炮集記錄進行線性掃描疊加.這類似于對數(shù)據(jù)體進行Tau-

p

變換，其中Tau為截距，

p

為射線參數(shù)，然后求取Taup域上的極值對應的

p

值，即可知原始數(shù)據(jù)窗口內最強能量對應的速度值.地表震源激發(fā)的虛擬單炮對應著高絕對

p

值（低速），地下震源激發(fā)的虛擬單炮對應著低

p

值（高速），分別對應著面波和體波能量（圖12）.

圖12 傾斜疊加照明分析示意圖（修改自Vidal et al., 2014）.（a）彈性地下模型，黃色倒三角地表檢波器，黑點代表地下震源；（b）淺源數(shù)據(jù)（灰色星）重構的虛擬炮集；（c）虛擬炮集圖（b）對應的慢度；（d-e）深源數(shù)據(jù)（青色星）的虛擬炮集和慢度，綠線表示面波和體波慢度的區(qū)分界限；（f）地下深源數(shù)據(jù)傾斜疊加照明分析的慢度分布圖Fig.12 Slant-stack illumination diagnosis (modified from Vidal et al., 2014).(a) Elastic subsurface model with yellow triangles as receivers and black dots as sources.(b) The virtual source retrieved from the shallow source is marked by the grey star.Red lines indicate slownesses.Green lines are pre-defined limits between body and surface wave slownesses in the virtual source function.(c) Slowness representation of the virtual source function from (b).(d) As in (b), but for the source represented by the cyan star in (a).(e) As in (c), but from (d).(f) Illumination diagnosis, consisting of the results for the correlated common source panels from sources in (a) (black dots), with the panel from (b) (grey star) discarded, and the panel from (d) (cyan star) included

Chamarczuk等（2019）把該方法拓展到三維數(shù)據(jù)，研究分析了地表震源S1和地下震源S2對于5個檢波器X1 至 X5組成的十字排列的走時響應，S1和S2震源對于X4、X2、X5組成的inline方向測線時距曲線相同，也就意味著難以從inline方向加以區(qū)分能量來自地表還是地下震源.但是兩個震源對X1、X2、X3組成的crossline方向測線到時卻有很大差別，其中地表S1震源是線性時距曲線，而地下S2震源是雙曲時距曲線（圖13）.因此不單單掃描inline測線方向的速度，同時也掃描crossline方向的速度，就能更好區(qū)分震源是來自于地表面波還是地下體波的能量.該研究表明，實際勘探中若能同時部署crossline測線，將更有助于從速度差異特征上準確識別和提取弱體波信號（Chamarczuk et al., 2019）.

圖13 三維觀測系統(tǒng)更好區(qū)分面波和體波示意圖（修改自Chamarczuk et al., 2019）.（a）五個檢波器 X1 至 X5 記錄的來自地表震源S1和地下震源 S2 的直達走時.兩源同時激發(fā)，且假定介質傳播速度恒定5 km/s.（b）由五個檢波器形成的十字排列觀測系統(tǒng).（c）由三條檢波線組成的三維觀測系統(tǒng)Fig.13 The diagram of event discrimination from surface- or body-waves (modified from Chamarczuk et al., 2019).(a) Direct arrivals from sources S1 and S2 recorded at five receivers X1 to X5 forming a cross-shaped array.We assume a constant propagation velocity of 5 km/s.(b) Configuration of the cross-shaped array formed by the five receivers.(c) 3D survey formed by the three receiver lines

2.3 基于方位角差異被動源體波

被動源實際觀測的時候，地下震源往往不是均勻分布的，其中某些特定方位的波動可能貢獻了大部分能量，因此可根據(jù)方位角差異來區(qū)分面波和體波.人們利用波束成形（beamforming）方法（Aki and Richards,1980; Cheraghi et al., 2015; Girard and Shragge, 2019, 2020）可以求取被動源數(shù)據(jù)中不同頻率成分的波的方向和速度，再據(jù)此去除干擾面波和分離出有效反射波.計算公式如下（Cheraghi et al., 2015）：

式中，

B

(

v

,θ)為計算出的波束（ θ為傳播方向，

v

為傳播的平均速度），

f

表示信號的頻率，

N

是接收器的數(shù)量，?

f

是頻率帶寬，

f

是噪聲的中心頻率，

s

表示第

i

個檢波器接收信號的傅里葉變換，其檢波器的坐標分別為

x

、

y

.例如圖14是加拿大一個CO存儲監(jiān)測項目實測數(shù)據(jù)的計算波束功率，圓的方位角表示波的來源方向，圓的半徑代表不同的波速，振幅的大小表示該點（每個點都對應著一個速度和方向）的歸一化波束功率值.該例子中，工區(qū)附近240°和 330°方位角之間的發(fā)電廠是環(huán)境噪聲的主要來源，而且速度能量明顯大于3 000 m/s，這部分能量以直達波、折射波、反射波等體波能量為主，有助于后續(xù)的被動源反射體波信號提取（Cheraghi et al., 2015）.

圖14 波束成形的方向性分析（修改自Cheraghi et al.,2017）.正北方向為0°，角度沿逆時針方向增加Fig.14 Directional Beam-forming analysis (modified from Cheraghi et al., 2017).The north had an azimuth of 0°,which increased to the west (i.e., counterclockwise)

人們還可利用波束分析方法開展更精細研究，比如Cheraghi等（2015）對加拿大Lalor礦區(qū)的被動源數(shù)據(jù)細分成不同頻率和時段來研究波動的主要方向和速度.首先將數(shù)據(jù)按小時段劃分，并對各時段數(shù)據(jù)分別進行1～3 Hz、10～13 Hz 和 23～25 Hz三個頻段的波束分析.圖15表示研究區(qū)一天之中每個時段對應的最大波束功率譜及其對應的方位角.可見，三個頻段的最大方位能量主要集中在兩個方向，低頻段能量主要來自于300°方位角（東北方向8 km）的人類活動，高頻段能量主要來自于約80°方位角（正西方向）的地下采礦活動，而10～13 Hz能量則來自于這兩個方向的綜合貢獻.基于波束分析的方位角信息，并結合速度差異特征（大于3.2 km/s的主要為體波能量），可以有效識別和提取弱體波信號，再對虛擬炮集數(shù)據(jù)進行處理獲得最終的反射疊加剖面（圖16b），與主動源剖面（圖16a）對比可知，盡管被動源剖面的分辨率和精度要低，但該礦區(qū)的三套主要反射I1、J2、K3還是清晰可見.

圖15 不同小時段的最大波束功率及其視速度和方位角（修改自Cheraghi et al., 2015）Fig.15 Calculated beamformers for the maximum beam power of each hourly panel and its apparent velocity and azimuth (modified from Cheraghi et al., 2015)

圖16 加拿大Lalor礦區(qū)R141線的傾角時差校正疊加剖面（修改自Cheraghi et al., 2015）.（a）主動源剖面；（b）被動源剖面Fig.16 DMO stacked section along line R141 in Lalor deposit, Canada (modified from Cheraghi et al., 2015).(a) Active-source 3D survey and (b) passive interferometry survey

3 被動源反射地震數(shù)據(jù)處理和成像

3.1 被動源原始數(shù)據(jù)預處理

不同于天然地震學研究領域常用的寬頻帶地震儀（如60 s、120 s長周期地震儀），被動源地震勘探領域目前主要采用成本更為低廉、布設更為便捷、間距更為密集的節(jié)點式地震儀（如5 Hz、10 Hz檢波器），旨在記錄更多高頻的體波反射信號.因此在預處理階段，比如帶通濾波等處理，需根據(jù)檢波器的自然頻率以及研究目標的特征選擇合適的頻帶范圍.由于被動源數(shù)據(jù)往往信噪比較低，尤其是連續(xù)記錄中的一些強振幅噪聲會在虛擬炮集干涉計算過程中產生一些虛假振幅，因此需要在數(shù)據(jù)預處理階段盡可能加以去除.壓制強振幅噪聲一般是在對原始數(shù)據(jù)開展去均值、去趨勢、端點歸零等基本處理之后.其中一種較簡便方法是對數(shù)據(jù)進行歸一化處理，如Oren和Nowack（2017）在時間域使用sign-bit信號歸一化來降低噪聲權重，但該方法無法保留原始數(shù)據(jù)集的動態(tài)范圍；另一種常見策略是利用窗口進行數(shù)據(jù)選擇，即在互相關之前識別并去除（或降低權重）強噪聲的窗口（Panea et al.,2014），該方法需要對所有數(shù)據(jù)窗口進行判斷和選擇，工作量比較大.最近幾年，Girard和 Shragge（2020）、方捷等（2022）等開發(fā)自動識別和處理流程來提高效率，總體的處理原則是保持均方根振幅變化的趨勢.比如計算每個滑動時窗內的振幅平均值，若該窗內的樣點值偏離均值較大，就用均值代替當前的樣點值，或者直接定義振幅的一個閾值來將高振幅記錄降低到閾值以下.

除了連續(xù)強振幅噪聲，還有一些在時間域以孤立的尖峰或毛刺形式出現(xiàn)的噪聲，Girard和Shragge（2019）和方捷等（2022）使用基于L1迭代加權的最小二乘方法（Claerbout, 2014）處理.它采用擬合的思想對原始數(shù)據(jù)進行處理，可對高振幅信號進行壓制，而低振幅信號不受影響.從圖17可見，經(jīng)過處理后強振幅噪聲得以壓制，而數(shù)據(jù)的動態(tài)范圍仍然能夠保持.

圖17 時間域強振幅控制示意圖（方捷等，2022）Fig.17 Strong amplitude control in the time domain (Fang et al., 2022)

在時間域去除強振幅噪聲基礎之上，Girard 和Shragge（2019）又在頻率域壓制檢波點的固定頻率噪聲（如機電噪聲等）.該方法首先對數(shù)據(jù)進行傅里葉變換并計算信號的振幅譜

M

和相位 ?，見式（13）和（14），其中，R和I分別表示信號傅里葉變換的實部和虛部.然后與時間域處理類似，定義一個閾值將高譜幅值降低到閾值以下得到振幅譜

M

.最后根據(jù)公式（15）使用

M

重建信號：

式中，

F

是傅里葉反變換算子.圖18是頻率域強振幅控制的示意圖，可見經(jīng)過頻率域處理后強振幅噪聲得以進一步壓制.

圖18 頻率域強振幅控制示意圖（修改自Girard and Shragge, 2019）.（a）頻率域強振幅控制前（藍色）和后（紅色）的時域信號；（b）頻率域強振幅控制前（藍色）和后（紅色）的部分頻譜Fig.18 Strong amplitude control in the frequency domain(modified from Girard and Shragge, 2019).(a) Timedomain trace before (blue) and after (red) debursting.(b) Partial spectrum before (blue) and after (red)debursting

Girard和Shragge（2019）展示了預處理不同階段的虛擬道集（圖19），分別是基于原始數(shù)據(jù)、窗口篩選后數(shù)據(jù)，以及強振幅壓制（包括時間和頻率域）后的數(shù)據(jù).從圖19共偏移距排列的虛擬道集可見，完成預處理后的數(shù)據(jù)噪聲壓制明顯，其中0.5 s左右的反射同相軸的信噪比相對原始數(shù)據(jù)得以大幅提高.

圖19 互相關+疊加得到的以共偏移距（1.0 km）排列的虛擬道集（修改自Girard and Shragge, 2019）.（a）原始數(shù)據(jù)；（b）窗口選擇數(shù)據(jù)；（c）強振幅壓制數(shù)據(jù)Fig.19 Single extracted cross-correlation + stack trace at a virtual source-receiver offset of 1.0 km (modified from Girard and Shragge, 2019).(a) Raw data.(b) After data selection.(c) After time and frequency debursting

3.2 被動源虛擬炮集的數(shù)據(jù)處理

經(jīng)過預處理和面波、體波分離后構建的虛擬單炮盡管反射依然較弱，但可借鑒主動源反射地震數(shù)據(jù)的處理思路，利用多次疊加技術增強反射能量.被動源數(shù)據(jù)處理的主要目的是盡可能提高炮集反射的信噪比，改善最終成像質量，一般包括靜校正、振幅垂向補償和橫向均衡、隨機和相干噪聲壓制、多次波壓制、速度分析、切除、疊加和反褶積等步驟.

被動源數(shù)據(jù)的靜校正一般采用高程校正為主，將起伏地表虛擬炮集校正為水平地表虛擬炮集.值得注意的是，如果直接在被動源原始數(shù)據(jù)上先應用靜校正再相關合成虛擬炮集，會導致靜校正量發(fā)生改變，需對反射記錄非因果部分進行兩倍炮點校正，因果部分進行兩倍檢波點校正才能疊加生成正確的虛擬炮集（鄭鈞等，2022）.一般來說，通過被動源數(shù)據(jù)預處理過程中的去均值、去趨勢已較好完成單炮能量的重新均衡.在此基礎上，還可采用球面擴散補償和地表一致性振幅補償進行校正.

在虛擬炮集靜校正和振幅均衡基礎上可更好地開展噪聲壓制工作，其中隨機噪聲往往是被動源數(shù)據(jù)中面波之外最嚴重的噪聲，我們可通過原始數(shù)據(jù)篩選、去均值、去趨勢、強振幅控制、頻率濾波等處理進行較好壓制.除了隨機噪聲，虛擬炮集記錄中的錯誤旅行時、干涉噪聲還有多次波的存在都會大大影響后續(xù)地震處理和成像的精度.針對這些問題，Hu等（2020）和Wang等（2021）提出了一種迭代去噪方法，將聚焦變換引入到三維L1范數(shù)稀疏約束反演Closed-loop SRME（SRME：自由表面多次波壓制技術）算法框架，實現(xiàn)被動源虛擬炮集的相干噪聲和多次波一體化壓制.首先，選取估計的一次波構建聚焦變換的反算子，使聚焦變換具有多次波降階的作用，一次波變換為聚焦點，高階多次波變換為低階多次波，利用多次波降階的深層信息糾正，重構被動源虛炮集中錯誤的遠偏移距走時信息.其次，由于聚焦變換能夠將有效信號聚焦到聚焦點，使有效信號與噪聲的分布區(qū)域產生明顯的區(qū)別，通過多次迭代的方式構建聚焦域切除算子逐步壓制噪聲，最終獲得具有正確遠偏移距走時信息與高信噪比的一次波虛擬炮集.該方法克服了虛擬炮集中相干噪聲的限制，大大提高了噪聲型被動源數(shù)據(jù)處理與解釋的準確度.

Wang等（2021）通過對模擬數(shù)據(jù)試驗，驗證了基于聚焦變換和SRME三維稀疏反演方法對于被動源干涉噪聲和多次波壓制的有效性（圖20）.圖20a是用來試驗的一個五層模型，該模型內部存在一個斷層和凹陷式地層.299個隨機震源布設于地下1.2～1.5 km的位置（圖20b），用來模擬天然背景噪聲（圖20c）.圖20d為互相關生成的虛擬炮集，圖20e是利用該迭代去噪方法后的虛擬炮集，圖20f則展示用來對比的SRME處理后的主動源記錄.對比可知，經(jīng)過處理后的虛擬炮集走時信息和同相軸曲率獲得了糾正（白色箭頭），且干涉噪聲和多次波去除效果明顯（黑色箭頭）.經(jīng)過多種方法對不同噪聲壓制之后，還需對被動源炮集記錄兩側的虛假反射以及低速面波殘留能量進行切除和疊加來進一步提高信噪比.

圖20 基于聚焦變換和SRME三維稀疏反演的被動源數(shù)據(jù)去噪試驗（修改自Wang et al., 2021）.（a）數(shù)值試驗速度模型；（b）隨機震源采用的子波；（c）震源分布位置；（d）互相關生成的被動源虛擬炮集；（e）應用該方法處理的虛擬炮集；（f）SRME處理后的主動源炮集Fig.20 Numerical tests of noise removal via focal-denoising closed-loop SRME based on the 3D sparse inversion of passive data(modified from Wang et al., 2021).(a) Model for simulating noise-source data and (b) wavelets of the random sources.(c) Source location.Virtual shot gathers from the noise-source data using (d) cross-correlation and (e) the proposed method.(f) Active-source data using SRME

此外，速度分析對最終疊加或偏移成像也十分重要，一般要對速度譜進行多輪、加密拾取.由于相關合成虛擬炮集過程中，可能會產生不同頻率的干擾波，在生成速度譜時利用帶通濾波只保留頻率能量最集中的部分，這樣有利于速度譜能量團的拾取.圖21是一個實際被動源資料的速度分析例子.速度譜中有明顯的能量團，共中心點道集有相對應的反射同相軸，而局部疊加結果中也可見連續(xù)性的反射（劉國峰等，2021）.

圖21 被動源反射信號速度分析和局部疊加（劉國峰等，2021）Fig.21 Reflection velocity analysis and stacked local passive data (Liu et al., 2021)

3.3 被動源數(shù)據(jù)直接偏移成像

目前被動源成像一般采用間接偏移的方法，即首先將被動源數(shù)據(jù)重構為主動源虛擬炮集，再進一步采用主動源地震的處理和偏移方法來成像.但大部分地震干涉法重構的炮集仍然面臨著數(shù)據(jù)量大、信噪比低以及假設條件過多等情況.另外一種思路是Artman（2006）最早提出的被動源直接偏移成像方法，他采用噪聲源數(shù)據(jù)本身作為邊界條件，在頻率域直接利用被動源傳動波場數(shù)據(jù)采用單程波偏移，得到了地下地質體構造圖像.直接偏移成像方法免去了重構虛擬炮集時產生的二次噪聲，提高了成像的準確性，不過該方法對于速度模型的精度要求比較高.

被動源直接偏移方法本質上是利用了多次波信息進行偏移成像，人們又將該方法擴展到逆時偏移成像上，從而適應更復雜的速度模型和高陡地下構造.Zheng等（2015）將逆時偏移結合多次波成像的邊界條件應用于被動源地震數(shù)據(jù)，直接利用數(shù)據(jù)本身正傳反傳并進行偏移成像.Cai 等（2019）提出了被動源逆時偏移能量歸一化成像條件，并采用波場分離的方法較好消除逆時偏移的低頻噪聲，同時利用并行分段疊加的成像策略提高計算效率.圖22的數(shù)值試驗表明該方法能較好克服被動源震源分布不均勻導致的成像能量比例失調問題，改善了模型深部與邊緣成像的質量.

圖22 被動源地震數(shù)據(jù)直接偏移的數(shù)值試驗（修改自Cai et al., 2019）.（a）Marmousi模型，震源分布在左下方的黑框紅框；（b）隨機截取的一段4 s的被動源記錄；（c）基于全波場分離歸一化的被動源數(shù)據(jù)直接逆時偏移成像Fig.22 Numerical tests for the direct migration of passive data (modified from Cai et al., 2019).(a) Marmousi model with the passive data generated from source locations marked by black and red boxes.(b) Example of a passive seismic record with a length of 4 s.(c) The direct migration image by full wavefields decomposition normalized reverse time migration

Girard和 Shragge（2020）利用加拿大Lalor 湖附近礦區(qū)的被動源實際數(shù)據(jù)開展三維直接偏移成像研究.圖23展示了主、被動源成像結果在礦體中心位置的[

E

,

N

,

Z

] = [2.36, 2.54, 0.8] km提取的水平切片和兩個垂直深度剖面.可見，在0.5 km和1.5 km深度范圍內，被動源三維數(shù)據(jù)體揭示工區(qū)西部是NE傾向反射，而東部則發(fā)育西傾構造，與主動源數(shù)據(jù)體的解釋基本一致.被動源直接成像因為無需重構炮集，因此在諸如隨鉆測井、礦區(qū)現(xiàn)場等處理時效性要求高，背景速度又比較確定的研究場景具有較大應用潛力.

圖23 加拿大Lalor湖附近礦區(qū)主、被動源水平切片和垂直剖面（修改自Girard and Shragge, 2020）.（a）被動源數(shù)據(jù)的直接偏移成像；（b）主動源數(shù)據(jù)的偏移結果Fig.23 Comparison of the horizontal slice and vertical section of 3D active and passive surveys in Lalor deposit, Canada (modified from Girard and Shragge, 2020).(a) Passive direct migration stack and (b) the active-source migration volume

4 被動源反射波勘探應用實例

隨著便攜節(jié)點地震儀等采集設備的日益普及，以及被動源反射地震處理成像方法的發(fā)展，被動源反射波勘探逐漸從理論探索走向實際應用.目前人們在開展被動源反射勘探實際研究時，往往還進行重合測線的主動源采集，這樣可通過主、被動源成像剖面的對比，進一步分析被動源反射方法的成像精度，以及探討其在不同勘探領域的應用潛力.

4.1 CO2封存場地監(jiān)測研究實例

Cheraghi等（2017）利用被動源地震數(shù)據(jù)對加拿大Saskatchewan地區(qū)的一個CO封存場地開展研究，分析被動源地下結構成像的精度，以及被動源方法用于時移地震監(jiān)測的潛力.該項目沿兩條垂直檢波線共采集兩個時間段（2014年6月的23天、2015年2月的13天）的被動源數(shù)據(jù).圖24展示了L1線23天被動源數(shù)據(jù)的剖面結果，其在0.8～1.5 s目的層范圍內成像精度較高，不僅在一定程度上可與帶通濾波后的主動源地震剖面相比較，甚至與注入井的合成記錄的相關系數(shù)也達到0.5～0.65.特別是一些關鍵層位，比如下科羅拉多頁巖（LC）（0.85 s，約 930 m）、Vanguard 組砂巖（VF）（1.1 s，約 1 100 m）和 Bakken 頁巖（BK）（1.45～1.5 s，約2 100 m）等標志層的反射響應均很明顯.不過受限于信噪比，被動源剖面無法識別更深部的前寒武地層（PC）.

Roach等（2015）研究表明CO地下儲層的監(jiān)測研究需要重復度很高的時移地震數(shù)據(jù)，即不同時段采集的地震振幅變化要控制在10%～20%，這樣才能刻畫CO儲層性質隨著時間推移的變化（預期波阻抗變化范圍約為-8%～-17%）.盡管該項目的被動源成像對關鍵層位的結構信息刻畫比較清楚（圖24），但是其精度還無法滿足時移地震監(jiān)測的需求，因為不同時段的成像結果對比表明，23天采集數(shù)據(jù)獲得的剖面與時隔8月后的13天剖面成像差異很大（圖25），而這顯然不可能是CO儲層性質變化所致.由于兩次被動源成像采用的檢波器和處理流程都相同，推測成像的巨大差異可能由于不同的采集時長、不同噪聲強度以及噪聲功率譜隨時間的變化（Mehta et al., 2008）等原因引起.該研究表明目前被動源反射還只能達到構造成像的精度，尚很難通過被動源的時移反射地震來研究地下儲層的性質變化，但Cheraghi等（2017）仍認為他們的研究可以通過使用更長持續(xù)時間的噪聲記錄來加以改進，且被動源地震由于其低成本優(yōu)勢，如果能開展不間斷的連續(xù)記錄，也許能在地下結構的長期動態(tài)監(jiān)測方面更有作為.

圖24 加拿大Saskatchewan的CO2封存場所的被動源研究（修改自Cheraghi et al., 2017）.（a）被動源剖面（紅色為注入井的合成記錄）；（b）井的合成記錄（子波從被動源數(shù)據(jù)提取）；（c）主動源數(shù)據(jù)剖面；（d）深時轉換的P波速度測井記錄.LC：下科羅拉多頁巖；VF：Vanguard組砂巖；BK：Bakken頁巖；PC：前寒武紀巖石Fig.24 Passive surveys of aquistore CO2 storage site Saskatchewan, Canada.(modified from Cheraghi et al., 2017).(a) Passive section (red curve is the log-based synthetic seismogram).(b) Synthetic seismogram generated from a wavelet extracted using the passive data convolved with borehole log-based reflectivity.(c) Active-source cube.(d) Depth-to-time converted P-wave velocity log.LC, VF, BK, and PC represent the Lower Colorado shale, Vanguard Formation, Bakken shale, and Precambrian,respectively

圖25 CO2封存場所的L1和L2線不同時段的被動源地震成像（修改自Cheraghi et al., 2017）Fig.25 Passive seismic imaging for different time periods at the CO2 storage site along lines L1 and L2 (modified from Cheraghi et al., 2017)

4.2 金屬礦區(qū)應用實例

金屬礦床多位于結晶巖地區(qū)和構造發(fā)育帶，地下地質情況復雜，且金屬礦體形態(tài)復雜，非均質性極強，一直是反射地震勘探的難點，相較于油氣勘探領域，金屬礦的主動源地震技術進展較為緩慢（湯聰?shù)龋?022）.不過對于被動源反射波成像來說，金屬礦集區(qū)特別是那些正在開采的礦區(qū)因為能夠提供十分關鍵的地下噪聲源，反而是比較適合被動源勘探應用的潛在領域.

Chamarczuk等（2022）最近在芬蘭東部的一個結晶巖多金屬（Cu-Co-Zn-Ni-Ag-Au）礦床Kylylahti開展了一項被動源三維地震成像工作.該研究在3.5 km× 3 km的區(qū)域內以200 m線距和50 m點距總共布置了994個垂直檢波器，并持續(xù)記錄30天的被動源噪聲數(shù)據(jù).數(shù)據(jù)采集期間日常采礦活動持續(xù)進行，尤其是每天100～800 m深處的雷管爆破提供了大量地下噪聲源.圖26是基于所采數(shù)據(jù)構建的三維地震虛擬炮集與主動源數(shù)據(jù)的對比（均顯示其中6條檢波線）.由于礦區(qū)地下構造異常復雜，即使在主動源數(shù)據(jù)的單炮中也較難識別有效反射信號，好在開展的是三維地震采集，相鄰檢波線的數(shù)據(jù)可用來輔助追蹤有效反射.總體上，主、被動源單炮的信噪比略有差異外，主動源單炮上的反射信號在被動源虛擬炮集也有所指示（圖26綠色箭頭），但在大偏移距的被動源炮集上則存在著一些假象（圖26黑色箭頭）.

圖26 主、被動源三維共炮集數(shù)據(jù)對比（修改自Chamarczuk et al., 2022）.（a）主動源炮集；（b）30天記錄重構的被動源虛擬炮集；（c）10天記錄重構的被動源虛擬炮.VSG：虛擬炮集Fig.26 Comparison of 3D common-source gathers using active and passive data (modified from Chamarczuk et al., 2022).(a) Activeshot gather, (b) VSGs obtained using 30 d of noise, and (c) VSGs obtained using 10 d of noise.VSG:Virture Shot Gathers

基于重構虛擬炮集開展后續(xù)的三維反射地震數(shù)據(jù)處理，獲得了最終的疊后偏移剖面（圖27）.對比發(fā)現(xiàn)，主、被動源的剖面具有很好一致性，并且與該區(qū)基于鉆孔的地質解釋成果較為吻合，比如圖27中多處紅色箭頭所指的巖體邊界均能看到較明顯反射能量.特別要指出的是，被動源剖面在該礦區(qū)有效彌補了主動源地震由于地形限制無法采集而導致的淺部數(shù)據(jù)缺口部分（圖27藍框部分），其中箭頭4所示的位置清楚揭示了巖體在淺部的走向和邊界.該研究證明了被動源反射地震勘探方法在金屬礦地區(qū)有較大的應用潛力，因為最終成像不僅與基于鉆孔的地質解釋吻合度較高，甚至還提供了主動源地震成像缺失的一些關鍵信息.

圖27 芬蘭Kylylahti結晶巖礦集區(qū)主、被動源疊后偏移剖面（修改自Chamarczuk et al., 2022）.（a） 主動源剖面；（b）主動源剖面疊合地質解釋；（c）被動源剖面疊合地質解釋KAL沉積（藍色）；OUM巖體（綠色）；S/MS礦化半礦化巖體（紅色）Fig.27 Comparison of post-stack migrated sections obtained from active and passive surveys in Kylylahti polymetallic mine, Finland (modified from Chamarczuk et al., 2022).(a) Active survey; (b) Active survey with the geological model background;(c) Passive survey with the geological model background.Geology in the background is color-coded as follows: Kaleva Sedimentary Belt (KAL): blue; semi-massive to massive sulfide (S/MS) mineralization: red; Outokumpu ultra-mafics (OUM):green

4.3 煤礦地區(qū)應用實例

除了金屬礦勘探方面的應用之外，最近Ning等（2021）還將被動源反射成像技術應用在某煤礦地區(qū)的地下結構成像.該研究在600 m× 800 m 的區(qū)域使用了202個檢波器，間距為50～100 m，連續(xù)觀測記錄27天，期間煤礦開采活動正常進行.值得說明的是，該研究只選取其中一天的數(shù)據(jù)以兩種不同方式開展試驗，一種是首先應用波形偏移處理和tomoDD 地震定位方法 （Miao et al., 2021; Zhang and Thurber, 2003）從連續(xù)波形記錄中檢測到 576個高信噪比的地震事件，然后只利用包含這些高信噪比地震事件的記錄進行虛擬炮集重構；另一種則是將連續(xù)地震記錄以30 min時窗分段再互相關進行虛擬炮集重構.該研究采用Chamarczuk等（2021）的方法對兩種方法重構的虛擬炮集進行相似度定量分析，其中相似度定義的范圍在 -1～1之間變化，值越大說明兩者越相似.實際計算所得的兩組虛擬炮集的相似度高于0.7，說明兩組炮集總體非常相近，也從側面反映該研究提取的體波信號可信度較高.

基于提取體波信號重構虛擬炮集，并進行后續(xù)處理最終獲得研究區(qū)的柯?；舴驎r間偏移剖面（圖28）.兩條不同檢波線的時間偏移成果均有兩個明顯的反射層（黑色虛線表示），對應著區(qū)域的煤層構造.其中較淺的煤層（箭頭“1”表示）連續(xù)性較差，較深的煤層（箭頭“2”表示）連續(xù)性相對較好，反映了該區(qū)較淺的煤層開采程度較高的現(xiàn)狀.另外工區(qū)自西向東，煤層呈現(xiàn)出東傾的趨勢，奧陶系地層處于第2個煤層之下，與該區(qū)已有的地質認識和工作面采掘資料吻合良好.該項研究表明正在開采中的煤礦能產生豐富的地下噪聲源，這非常有助于體波信號的提取，被動源地震反射成像方法可以用來輔助調查區(qū)域內的煤層構造賦存形態(tài)和采空區(qū)范圍.實際工作中還可進一步將被動源反射地震剖面與面波成像的速度反演結果等進行聯(lián)合解釋以提高研究精度.

圖28 不同檢波線的偏移剖面（修改自Ning et al., 2021）.黑色虛線表明煤層成像位置Fig.28 Migration sections of receiver lines (modified from Ning et al., 2021).The black dashed lines show the coal seams imaged using ambient-noise data

5 結論與展望

本文主要介紹了被動源反射波勘探成像的一些關鍵技術，包括地震干涉法重構虛擬炮集、基于多種差異特征分離面波和體波以及被動源原始數(shù)據(jù)預處理、虛擬炮集多次波壓制和直接偏移成像等方法.如何在面波能量所主導的大量噪聲連續(xù)記錄中高效率地甄別和提取弱體波信號，是被動源反射地震成像的核心問題.通常需要在干涉和疊加之前對原始數(shù)據(jù)不同時窗內的能量加以判斷和選擇，數(shù)據(jù)存儲和計算量都十分巨大.為了提高效率和精度，人們已開發(fā)各種自動化流程加以處理，同時大數(shù)據(jù)處理中的圖像識別、機器學習、人工智能等技術已經(jīng)被用于并將進一步在被動源反射地震勘探成像中發(fā)揮更重要作用（Chamarczuk et al., 2019）.

隨著勘探技術的不斷發(fā)展，一些新方法也被應用于被動源反射波成像.比如將地震超越干涉法的Marchenko成像應用到被動源數(shù)據(jù)，該方法能夠直接面向目標成像，避免檢索炮集中自由表面多次波和內部多次波的影響（Jin, 2020）.同時，以全波形反演為代表的高分辨率反演技術也被引進到被動源成像之中，如利用隨機干涉重構（尚旭佳等，2021）和多尺度震源編碼反演（蔡中正，2021）等技術獲得較為準確的地下介質速度模型，再用該模型對被動源數(shù)據(jù)利用逆時偏移方法進行直接成像.而Shiraishi和Watanabe（2022）最近將被動源聲波逆時偏移推廣到彈性波逆時偏移，并獲得了P波和S波的被動源成像.此外，主、被動源數(shù)據(jù)的聯(lián)合反演和成像技術也在不斷發(fā)展，比如可以利用被動源數(shù)據(jù)豐富的低頻信號來降低高速屏蔽層對地震信號的散射和屏蔽作用，通過主、被動源據(jù)的聯(lián)合應用，使其在照明和頻帶上優(yōu)勢互補，能顯著提高成像和反演質量（張盼等，2019）.

除了成像方面研究之外，被動源數(shù)據(jù)的采集也值得關注.噪聲源的分布是影響干涉重構效果的關鍵因素，如有可能盡量選擇那些地下噪聲源豐富的區(qū)域開展被動源反射成像研究，比如正在開采的金屬礦區(qū)和煤田等.同時被動源勘探進行數(shù)據(jù)采集部署時應特別注意人工噪聲源的分布角度，條件允許時應盡量部署十字排列或者三維觀測系統(tǒng).相比單純的二維采集，多出的另一個方位信息有助于更好區(qū)分記錄到的能量主要來自地表還是地下（Chamarczuk et al., 2019）.而在特定目標區(qū)進行主、被動源的混合采集也頗有前景，比如對構造穩(wěn)定、缺少足夠多地下噪聲源的地區(qū)進行被動源勘探時，補充少量主動源數(shù)據(jù)就可以有效控制和改善被動源的成像效果.同時在稀疏炮點的主動源勘探中，有效利用被動源的信息也能夠在成像中增加更多的細節(jié)信息，提高成像質量（張盼等，2015）.

被動源反射地震勘探，最近在地下礦產勘查等結構成像領域已初步取得實效，而在碳封存等需要動態(tài)監(jiān)測儲層性質變化的勘探領域仍在不斷探索之中.盡管目前被動源反射成像與主動源成像的效果尚存在明顯差距，但是作為一種低成本的綠色環(huán)保方法，且被動源勘探的部署和應用方式非常靈活，既可作為主動源地震的補充，又能結合已有的被動源面波速度成果進行聯(lián)合分析，因此該技術的發(fā)展一直持續(xù)受到人們關注.特別是近年來由于節(jié)點檢波器、分布式光纖等技術的日益普及，連續(xù)接收被動源地震記錄成本已大幅度降低，以及高性能計算的大規(guī)模興起，使得海量數(shù)據(jù)的存儲和計算成本快速下降，被動源勘探在成本方面的優(yōu)勢得以進一步凸顯，也在一定程度上刺激著相關處理新技術和成像新方法的不斷發(fā)展，人們應該可以期待被動源反射地震勘探在未來會有更多更好的實際應用.

致謝

感謝三位評審專家提出的修改完善建議，感謝編輯部的大力支持.