基于SVD和PCA方法重構干涉格林函數

張 俊， 熊章強，2， 張大洲，2， 陳 杰， 劉云昌

(1.中南大學 地球科學與信息物理學院，長沙 410083；2.中南大學 有色金屬成礦預測教育部重點實驗室，長沙 410083)

0 前言

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

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

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

1 地震干涉理論與算法

1.1 地震干涉法(SI)

(1)

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

(2)

(3)

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

(4)

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

(5)

式中M是源的數量。

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

(6)

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

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的行向量，它的每一行表示各源傳播到兩點接收記錄的互相關，各源的互相關記錄構成如圖1所示的互相關譜。

圖1 兩點間源的互相關譜Fig.1 The cross-correlogram for sources

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

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

(9)

式中：T為轉置;s=eUS為右奇異值向量V的權系數，稱之為疊加系數。而當sk=0,M

(10)

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

(11)

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

1≤k≤j

(12)

1.3 主成分分析(PCA)

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

(13)

根據主成分分析的理論[13]，首先對X進行均一化

(14)

(15)

2 實驗結果與分析

2.1 有效震源分析

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

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

圖3 合成數據處理結果Fig.3 Synthesized data processing results(a)觀測系統；(b)原始互相關譜；(c)解析解與直接疊加對比;(d)SVD重構互相關譜；(e)解析解與SVD重構解對

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

2.2 數值模擬及數據處理

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

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

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

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

(16)

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

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

2.3 實測數據處理

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

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

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

圖6 兩道實測數據處理結果Fig.6 Field data processing results(a)原始互相關譜;(b)SVD重構互相關譜;(c)SVD+PCA重構互相關譜;(d)四種格林函數對比

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

3 結論

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

1)在均勻各向同性介質中，兩個接收點連線方向上的震源對格林函數起主要作用，而兩個接收點橫向兩側的震源對格林函數的計算結果有不同程度的影響。

2)使用SVD和PCA結合方法重構互相關譜可有效壓制部分干擾震源對格林函數的影響，而分離提取格林函數的主要成分，可有效提高計算近似格林函數的準確性。