張 俊， 熊章強(qiáng)，2， 張大洲，2， 陳 杰， 劉云昌

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083；2.中南大學(xué) 有色金屬成礦預(yù)測教育部重點(diǎn)實(shí)驗(yàn)室，長沙 410083)

0 前言

地震干涉法[1](Seismic Interferometry，SI)被用來計(jì)算兩個接收點(diǎn)之間的近似格林函數(shù)，其基本假設(shè)是其中一個點(diǎn)作為虛擬源，計(jì)算空間內(nèi)不同源點(diǎn)下兩接收點(diǎn)之間記錄的互相關(guān)，然后將不同源點(diǎn)的互相關(guān)結(jié)果疊加就可得到標(biāo)準(zhǔn)的干涉格林函數(shù)。上述所涉及到的源可以是人工源[2]、天然地震源[3]、隨機(jī)噪聲源[4-6]等。由于地震干涉法的部分震源所產(chǎn)生的地震波在傳播過程中會相互干涉產(chǎn)生偽影，從而影響計(jì)算得到的格林函數(shù)的準(zhǔn)確性。因此，計(jì)算出一個更加準(zhǔn)確的格林函數(shù)對于地震勘探，尤其是微動勘探的數(shù)據(jù)處理就顯得非常重要。

通常為了得到兩個接收點(diǎn)之間準(zhǔn)確的格林函數(shù)，一般要求檢波器接收到來自一個封閉曲面內(nèi)的震源，震源類型包括單極子和偶極子源。自2004年以來，Snieder、Weaver、Wapenaar等[4，7-10]的研究結(jié)果表明：兩個接收點(diǎn)之間所計(jì)算的格林函數(shù)僅由部分震源所產(chǎn)生的地震波來獲得。對于在菲涅爾帶內(nèi)震源采用固定相方法做近似積分計(jì)算格林函數(shù)可以得到較好的結(jié)果。而處于菲涅耳帶外的震源則會產(chǎn)生快速發(fā)散的能量從而破壞干涉能量。對于整個區(qū)域而言，震源在覆蓋不完整或者沒有偶極子源的情況下，恢復(fù)得到的格林函數(shù)質(zhì)量就相對較差。然而在實(shí)際的工作中，很少使用偶極子源，源分布通常也是不完整的。因此，目前在提高格林函數(shù)質(zhì)量時，重點(diǎn)在于地震波傳播時所到達(dá)的時間(延時)，而不是恢復(fù)傳播過程中的真實(shí)振幅。

為了提高干涉格林函數(shù)計(jì)算的準(zhǔn)確度，解決源覆蓋不完整的問題，Gabriela Melo等[10]使用奇異值分解方法對疊前互相關(guān)譜進(jìn)行分解，重構(gòu)一個秩為“1”的矩陣后計(jì)算得到近似格林函數(shù)，解決某些情況下震源不完整覆蓋下的問題，取得了較好的效果。筆者使用奇異值分解[10-12](Singular Value Decomposition，SVD)并結(jié)合主成分分析[13-15](principal component analysis，PCA)兩種方法計(jì)算近似格林函數(shù)，通過理論分析、數(shù)值模擬以及實(shí)測數(shù)據(jù)處理，分析對比格林函數(shù)在計(jì)算過程中部分源產(chǎn)生干擾影響，使用兩種方法結(jié)合來壓制干擾、提高格林函數(shù)計(jì)算的準(zhǔn)確性。

1 地震干涉理論與算法

1.1 地震干涉法(SI)

(1)

式(1)中包含兩點(diǎn)間所有單極子源和偶極子源的互相關(guān)。為了簡化處理流程，首先在積分運(yùn)算中減少一個互相關(guān)因子，其次只考慮單極子源。假設(shè)介質(zhì)是均勻的(速度c和密度ρ為常量)，其外部邊界為S，且沒有來自邊界外部的能量，則方程簡化為式(2)

(2)

(3)

式中α(x)是x和法向量在邊界S射線之間的角度。假設(shè)S足夠大，則有α(x)≈0和cos[α(x)]≈1，則公式(2)可簡化為式(4)。

(4)

由于在公式推導(dǎo)過程中做了相應(yīng)的簡化，使計(jì)算得到的格林函數(shù)丟失了真實(shí)振幅，但相位沒有損失，且式(4)適用于大多數(shù)地震干涉。為了提高地震波旅行時間的準(zhǔn)確度，邊界從連續(xù)到離散情況下，忽略方程(4)的振幅因子2/ρc。因此干涉格林函數(shù)可以寫成各源得到的互關(guān)聯(lián)之和

(5)

式中M是源的數(shù)量。

設(shè)xi(1≤i≤M)為源位置，τ為時間域的互相關(guān)滯后時間。根據(jù)SI在頻率域的推導(dǎo)方程，可以由它寫出時間域的方程。因此，互相關(guān)函數(shù)C可以寫成：

(6)

其中C(xi,τ)為在源xi的作用下xA和xB求得的互相關(guān)。A、B兩點(diǎn)間的格林函數(shù)G為各源的互相關(guān)之和，則有

G=G(xB,xA,t)+G(xB,xA,-t)=

(7)

1.2 奇異值分解(SVD)

G=G(xB,xA,t)+G(xB,xA,-t)=eC

(8)

式中：e為m個全為1的列向量；C為m個長度為2n-1的行向量，它的每一行表示各源傳播到兩點(diǎn)接收記錄的互相關(guān)，各源的互相關(guān)記錄構(gòu)成如圖1所示的互相關(guān)譜。

圖1 兩點(diǎn)間源的互相關(guān)譜Fig.1 The cross-correlogram for sources

由奇異值分解法，則式(8)的格林函數(shù)可以寫成：

G=e1×mCm×(2n-1)=

(9)

式中：T為轉(zhuǎn)置;s=eUS為右奇異值向量V的權(quán)系數(shù)，稱之為疊加系數(shù)。而當(dāng)sk=0,M

(10)

將求得的sk按絕對值大小排序，求貢獻(xiàn)率Qj為式(11)。

(11)

在找出一個適合的貢獻(xiàn)率Qj后就可以得到j(luò)，然后根據(jù)一個低秩的矩陣近似等于原矩陣的概念，保持S的前j個奇異值Sj重構(gòu)得到一個低秩的近似矩陣。則重構(gòu)的干涉格林函數(shù)為式(12)。

1≤k≤j

(12)

1.3 主成分分析(PCA)

主成分分析(PCA)是一種分析、簡化數(shù)據(jù)集的方法，可以用來降低數(shù)據(jù)集的維數(shù)，同時保持?jǐn)?shù)據(jù)集中的對方差貢獻(xiàn)最大的特征[14-15]。對重構(gòu)的低秩互相關(guān)譜Cj降維計(jì)算新的格林函數(shù)，假設(shè)一個矩陣為X，維度為m，各維度下具有n個樣本數(shù)。將X進(jìn)行線性變換變成另一個矩陣Y，使得Y的協(xié)方差矩陣為對角矩陣。Y是對原始矩陣X提取主成分后的協(xié)方差矩陣，用它可以對原始矩陣X分析主成分構(gòu)造并實(shí)現(xiàn)降維處理。對于矩陣X，其矩陣大小為m×n，則X可寫成式(13)。

(13)

根據(jù)主成分分析的理論[13]，首先對X進(jìn)行均一化

(14)

(15)

2 實(shí)驗(yàn)結(jié)果與分析

2.1 有效震源分析

如圖2(a)所示，在一均勻各向同性介質(zhì)的區(qū)域內(nèi)，位于A、B兩點(diǎn)設(shè)置兩個檢波器，在S1、S2兩點(diǎn)設(shè)置兩個震源。第一類源S1位于AB直線任意位置但非AB內(nèi)(文中設(shè)置在A左側(cè))，第二類源S2位于AB連線兩側(cè)。當(dāng)僅有第一類源S1時，計(jì)算得到的近似格林函數(shù)G1與真實(shí)格林函數(shù)G一致，如圖2(b)中的G和G1曲線所示，峰值時間剛好為源由A傳至B的旅行時差(延時)。當(dāng)源S1和S2同時存在時，計(jì)算得到的近似格林函數(shù)如圖2(b)中G2所示，可以看出G2與真實(shí)格林函數(shù)G有所不同。根據(jù)干涉原理可知，當(dāng)在S2點(diǎn)激發(fā)的地震波被A、B兩個檢波器接收后進(jìn)行互相關(guān)計(jì)算，其實(shí)質(zhì)是在S2和B點(diǎn)連線上的虛擬接收點(diǎn)A’所接收到的信號的互相關(guān)，此處點(diǎn)S2到點(diǎn)A的距離和點(diǎn)S2到點(diǎn)A’的距離相等。通過以上分析可知，當(dāng)震源位于兩個接收點(diǎn)連線的兩側(cè)時，計(jì)算得到的格林函數(shù)不是嚴(yán)格意義上接收點(diǎn)A和B之間的干涉格林函數(shù)，而且這個計(jì)算得到的干涉格林函數(shù)中信號的到達(dá)時間比真實(shí)的時間要小，也就是在計(jì)算得到的格林函數(shù)中在準(zhǔn)確格林函數(shù)之前還存在一個信號。這主要是因?yàn)樘摂M接收點(diǎn)A’到B點(diǎn)的距離較A點(diǎn)到B的距離較小引起的。

圖2 不同震源位置計(jì)算得到的格林函數(shù)對比Fig.2 Compare and analyze the green function(a)觀測系統(tǒng)圖;(b)格林函數(shù)對比圖

圖3 合成數(shù)據(jù)處理結(jié)果Fig.3 Synthesized data processing results(a)觀測系統(tǒng)；(b)原始互相關(guān)譜；(c)解析解與直接疊加對比;(d)SVD重構(gòu)互相關(guān)譜；(e)解析解與SVD重構(gòu)解對

通過對不同位置震源的近似格林函數(shù)峰值延時進(jìn)行分析發(fā)現(xiàn)，第一類源計(jì)算得到的峰值延時最大，與真實(shí)格林函數(shù)的峰值時間一致。對于格林函數(shù)GAB，源在接收點(diǎn)A、B中線左側(cè)的峰值為正延時，反之為負(fù)延時。在各向同性均勻介質(zhì)中第二類源計(jì)算結(jié)果會對格林函數(shù)產(chǎn)生不同程度的影響，其中靠近AB延長線兩側(cè)(雙曲線范圍內(nèi))的影響較小，它們的峰值延時也必定位于真實(shí)格林函數(shù)峰值時間的正負(fù)值之間，這將會引起計(jì)算得到的近似格林函數(shù)的峰值延時減小，而對于實(shí)測數(shù)據(jù)的處理，選擇一個參考具有重要的意義。

2.2 數(shù)值模擬及數(shù)據(jù)處理

在實(shí)測的地震數(shù)據(jù)處理應(yīng)用中，臺站接收到的信號來自各個方向，臺站連線兩側(cè)的振動信號將影響計(jì)算真實(shí)格林函數(shù)，因此采用SVD和PCA對重構(gòu)更加準(zhǔn)確的格林函數(shù)具有重要作用。我們使用數(shù)值模擬方法，對來自不同方向的地震記錄進(jìn)行處理，再使用PCA的方法分離出主信號(主信號)，求取得到新的格林函數(shù)。由于傳統(tǒng)標(biāo)準(zhǔn)格林函數(shù)計(jì)算中采用直接疊加互相關(guān)的方法得到的格林函數(shù)與真實(shí)格林函數(shù)有差別，其主要影響來自連線兩側(cè)的震源，在使用SVD分解時可以對其起到很好地壓制作用。如圖3(a)所示，在一個1 000 m×1 000 m的均勻各向同性空間區(qū)域內(nèi)不同位置處布置震源，其中在接收點(diǎn)A、B連線左側(cè)布置14個震源，橫向兩側(cè)分別布置5個和8個震源，共27個源。

圖4 加噪后重構(gòu)互相關(guān)譜對比Fig.4 Reconstructed correlogram with noise(a)原始互相關(guān)譜;(b)SVD重構(gòu)互相關(guān)譜;(c)PCA重構(gòu)互相關(guān)譜;(d)SVD+PCA重構(gòu)互相關(guān)譜

由圖3(b)和圖3(d)可知，SVD重構(gòu)得到的低秩互相關(guān)譜，較好地壓制了來自遠(yuǎn)離A、B連線方向震源(S15-S27)的能量，保留來自縱向及其附近震源的能量(S1-S14)。這說明在縱向附近的源對格林函數(shù)的影響較小。結(jié)合圖3(c)和3(e)可以看出，遠(yuǎn)離縱向的這部分能量影響了真實(shí)格林函數(shù)，在SVD方法處理后，影響格林函數(shù)的能量被很好地壓制。因此，更加確定來自A、B兩點(diǎn)間縱向的源對格林函數(shù)產(chǎn)生積極作用。

為了驗(yàn)證在有噪聲的情況下，PCA能夠更好地分離出有效信號。對上面相同情況下模擬的信號添加白噪聲，然后進(jìn)行處理分析(圖4)。這里以信噪比為標(biāo)準(zhǔn)進(jìn)行定量對比分析，定義信噪比為：

(16)

式中：S為有效信號；SN為加噪信號。

從圖4中可以看出，SVD對于壓制側(cè)向能量的影響效果明顯，而PCA對于提取主信號具有很好的效果。因此在SVD處理后互相關(guān)譜具有很好的同相軸性質(zhì)，再使用PCA進(jìn)行提取同相軸的互相關(guān)譜，能夠快速分離得到有效的格林函數(shù)，其互相關(guān)譜如圖4(d)所示。此外，還使用上述公式計(jì)算加噪后各處理方法得到格林函數(shù)與真實(shí)格林函數(shù)的性噪比，SVD處理后的性噪比為22.3 dB，PCA處理后的性噪比為18.5 dB，SVD和PCA同時處理后的性噪比為36.3 dB。

2.3 實(shí)測數(shù)據(jù)處理

為了驗(yàn)證上述理論分析及模擬結(jié)果的準(zhǔn)確性，在某處地形平坦區(qū)域使用錘擊作為震源，在不同位置進(jìn)行激發(fā)并采集數(shù)據(jù)，觀測系統(tǒng)如圖5所示。源s1～s10分布位于r1左右兩側(cè)4 m處的小圓內(nèi)；源s11～s16位于排列縱向右側(cè)，距離r1最近的s11有5 m，隨后每隔1 m激發(fā)一次；源s17～s33隨機(jī)位于排列左側(cè)橫向兩側(cè)(橢圓內(nèi))。檢波器共有12道，道間距為1 m，總共采集了33個不同位置源點(diǎn)的數(shù)據(jù)。圖5(b)和圖5(c)為1號檢波器和11號檢波器在不同震源激發(fā)下接收到的波形記錄，圖中可見不同位置激發(fā)的地震波先到達(dá)1號檢波器，然后才到達(dá)11號檢波器。

根據(jù)前面的分析可知，實(shí)測數(shù)據(jù)處理可以使用s11～s16中任意一個源的計(jì)算結(jié)果作為參考(選擇s11震源數(shù)據(jù)計(jì)算格林函數(shù)作為參考)，用r1和r11接收點(diǎn)的波形記錄計(jì)算近似干涉格林函數(shù)。在計(jì)算格林函數(shù)時先求得兩點(diǎn)間的互相關(guān)譜。圖6(a)、圖6(b)和圖6(c)分別為直接疊加、SVD及SVD+PCA結(jié)合處理后的互相關(guān)譜，圖6(d)為不同方法計(jì)算近似干涉格林函數(shù)進(jìn)行對比分析。

圖5 兩接收器波形記錄圖Fig.5 The waveform of two receivers(a)數(shù)據(jù)采集布置; (b)r1記錄; (c)r11記錄

圖6 兩道實(shí)測數(shù)據(jù)處理結(jié)果Fig.6 Field data processing results(a)原始互相關(guān)譜;(b)SVD重構(gòu)互相關(guān)譜;(c)SVD+PCA重構(gòu)互相關(guān)譜;(d)四種格林函數(shù)對比

從圖6(a)原始互相關(guān)譜可以看出，s11峰值時間、最大時間軸以及能量集中區(qū)域均在0.022 s附近，而s1～s10峰值位于0.01 s附近產(chǎn)生了較強(qiáng)的同相軸。在經(jīng)過圖6(b)的SVD分解后，橢圓內(nèi)的部分源被不同程度的壓制，但是s1～s10的同相軸仍然清晰可見。經(jīng)過對SVD重構(gòu)的低秩互相關(guān)譜采用PCA方法分離主信號得到圖6(c)所示的結(jié)果，然后在計(jì)算新的格林函數(shù)。圖6(d)為不同方法求得的格林函數(shù)對比圖，由圖6(d)中可知，直接疊加法將0.01 s附近的干擾加到了格林函數(shù)中，引起了0.005 s處相位改變以及0.018 s處的峰值時間右移減小。在使用SVD和PCA結(jié)合處理后，上述兩個位置的影響均得到了顯著改善，使得計(jì)算的格林函數(shù)更加接近真實(shí)情況。

3 結(jié)論

通過計(jì)算分析不同位置的震源對格林函數(shù)的影響以及進(jìn)行數(shù)值模擬與實(shí)測數(shù)據(jù)資料處理，得出如下結(jié)論：

1)在均勻各向同性介質(zhì)中，兩個接收點(diǎn)連線方向上的震源對格林函數(shù)起主要作用，而兩個接收點(diǎn)橫向兩側(cè)的震源對格林函數(shù)的計(jì)算結(jié)果有不同程度的影響。

2)使用SVD和PCA結(jié)合方法重構(gòu)互相關(guān)譜可有效壓制部分干擾震源對格林函數(shù)的影響，而分離提取格林函數(shù)的主要成分，可有效提高計(jì)算近似格林函數(shù)的準(zhǔn)確性。