基于稀疏表示的非零井源距VSP波場分離方法

趙彥青 蕭蘊詩

（同濟大學電子與信息工程學院，上海201804）

基于稀疏表示的非零井源距VSP波場分離方法

趙彥青*蕭蘊詩

（同濟大學電子與信息工程學院，上海201804）

隨著垂直地震剖面（VSP）技術的不斷發展，高保真的波場分離成為VSP資料處理與應用的關鍵之一。引入欠定盲源分離思想，利用VSP波場的極化特征及多分量資料在時頻域的稀疏性，提出一種基于稀疏表示的波場分離方法。該方法分為兩步：第一步將地震信號變換到時頻域，并估計極化矩陣；第二步在時頻域分離各單一波場，再反變換回時間域。在第二步波場分離過程中，針對時頻域波場稀疏性不足的情況，提出一種改進的最短路徑分解法，在一定程度上解決了時頻域局部信號混疊的問題。模型試算與實際資料處理結果表明，與傳統方法相比該方法能夠更準確地分離上、下行P-P波、P-SV波，分離后的波場保幅性較好，不存在混波、假頻以及邊界問題。

非零井源距VSP 多分量 波場分離 欠定盲源分離 稀疏表示

1 引言

垂直地震剖面（VSP）是一種地面激發、井中觀測的地震勘探技術。與地面地震資料相比，VSP資料信噪比和分辨率高，具有更明顯的地震波運動學與動力學特征。由于觀測方式的特殊性，VSP資料波場非常豐富，不同的波場能夠解釋不同的地球物理現象，因而高保真的波場分離是VSP資料處理的關鍵之一。

傳統的VSP波場分離方法有很多，按照數據類型可分為兩類：一類是針對單分量數據的波場分離，其中有的方法根據不同波場視速度差異進行數學變換，在變換域內完成波場分離再反變換，如f-k濾波、Radon變換等［1］，還有的方法對某一類型波場進行衰減或增強，達到波場分離的目的，如中值濾波、均值濾波等；另一類是基于多分量數據的波場分離，如根據不同波場極化特征的極化濾波方法、變速視慢度法等［2-5］。這些方法都有各自的優勢，但也存在一定的缺陷。常規f-k濾波、Radon變換法都會引起波場的畸變；中值濾波在分離上行P-P波、P-SV波時會造成有效波場損傷，并且存在邊界問題；常規極化濾波雖然具有較好的保幅性，但采用直達波傳播方向確定極化角僅在直達P波與上行P-P波極化方向完全正交時分離效果較好。針對傳統方法的不足，Wang等［6］提出了利用檢波點位置的速度與密度等參數在τ-P域實現多分量地震數據P波與SV波的分離；Blias［7］提出一種 wave-by-wave波場分離方法，根據振幅與視速度的不同實現逐波分離；王成禮［8］根據不同觀測方式下VSP的波場特征提出了不同的波場分離方法；杜婧等［9］通過估計局部斜率參數分離VSP波場，克服了傳統方法的濾波閾值鑲邊問題；崔汝國等［10］提出了浮動坐標系極化濾波，建立極化方向隨深度變化的濾波坐標系，提高了極化濾波的精度；孫文博等［11］提出一種基于參數反演的矢量波場分離方法，在分離波場的同時反演出速度、極化角等參數；馬志霞等［12］進一步對多種VSP波場分離方法進行了對比，取得了良好的應用效果；Jiang等［13］利用高精度拋物線Radon變換分離P-P波與P-SV波波場，提高了波場分離的精度。

盲源分離（BSS）理論是現代信號處理領域的一個嶄新的研究方向。這里的“盲”有兩層含義：一是源信號不能被觀測；二是傳輸信道不可知或者難以估計。盲源分離算法對系統的先驗知識要求甚少，具有廣闊的應用前景，近年被引入地震資料波場分離領域［14］。Van der Baan［15］提出了在τ-P域利用獨立分量分析（ICA）分離PP波與PS波，僅利用二者統計上的區別，而無須知道速度、密度等參數；Vrabie等［16］提出高階奇異值分解（HOSVD）和單峰獨立分量分析子空間算法，在波場分離過程中具有較好的抗噪性；Liu等［17］提出了基于信息最大化原則的地震信號盲源分離方法，利用高階統計量分離多次波。

將不同類型的VSP波場看作彼此獨立的源信號，多分量VSP波場分離問題屬于典型的欠定盲源分離（UBSS）問題，即觀測信號（VSP分量）個數小于源信號（VSP波場）個數。目前，基于稀疏表示（SR）的方法已經成為欠定盲源分離主要方法。本文利用VSP波場的極化特征，提出一種基于稀疏表示的VSP波場分離方法，該方法能夠同時分離出上、下行P-P波、P-SV波，既保留了極化濾波保幅性較好的優勢，又不存在二維濾波產生的波場畸變問題。

2 欠定盲源分離問題

假設有N個未知的源信號s1（t），s2（t），…，s N（t），經過一個未知的混合系統后，得到M個觀測信號x1（t），x2（t），…，x M（t），不考慮噪聲的影響，則瞬時線性混合模型的數學描述為

式中aij為混合系數。式（1）可以寫成矩陣形式

其中：X是M×T維的觀測信號矩陣；S為N×T維的源信號矩陣；A表示M×N維混合矩陣。盲源分離的主要任務就是在僅通過觀測信號X來獲取源信號S的估計。對于M=N的情況，經典的盲源分離方法通過尋找一個分離矩陣W來完成源信號的估計，即

對于非零井源距VSP（Offset VSP），定義三分量檢波器三個分量接收的記錄為r x、r y、rz，在各向同性介質中，SH波經過檢波器定位后完全被旋轉到y分量（即與炮檢平面垂直的方向），P-P波與P-SV波則旋轉到x-z平面（即炮點檢波點所在的平面）。因為SH波較容易分離，本文主要考慮P-P波與P-SV波的分離。為了書寫方便，仍用r x表示檢波器定位后的x分量。

圖1 檢波器定位后的VSP波場（x-z平面）

檢波器定位后的x-z平面如圖1所示。在某一時刻t，定義檢波器接收的上行P-P波為uP，上行P-SV波為uSV，下行P-P波為dP，下行P-SV波為dSV，水平分量為r x，垂直分量為rz，θuP、θuSV、θdP、θuSV分別為上行P-P波、上行P-SV波、下行P-P波、下行P-SV波傳播方向與x軸的夾角，即各個單一波場的極化角，存在如下關系

上式可以簡寫為

式中：r為t時刻x、z兩分量記錄；p為極化矩陣，其元素即為各個單一波場的傳播方向分別在x、z軸上的單位矢量投影；s為t時刻各上、下行P-P波、P-SV波單一波場。可以看出，在p未知、s也未知的情況下，僅憑r來估計s是典型的欠定盲源分離問題，極化矩陣p即為盲源分離問題中的混合矩陣，s為源信號。

3 基于稀疏表示的VSP波場盲分離

稀疏表示的數學模型與盲源分離模型一致，也可以用式（5）表示，稀疏表示的目標就是估計恰當的p與s，使s盡可能稀疏。稀疏性有兩層含義：其一是指源信號本身服從稀疏分布，大多數時刻取值為零或接近于零；其二是在某一時刻，只有很少的源信號取值較大。如果信號本身的稀疏性并不理想，則需要進行某種線性變換T，使信號在變換域內更具有稀疏性。由于T是線性變換，極化矩陣在變換前后保持不變，由式（5）可得

由于VSP波場復雜，單一波場在時間域的稀疏性較差，因而首先采用短時傅里葉變換（STFT）或小波變換將地震信號從時間域變換到時頻域，在時頻域進行波場分離后，再將分離的波場變換回時間域。

3.1 極化矩陣的估計

傳統的極化濾波方法認為直達P波的極化方向與上行P-SV波極化方向相同，而上行P-P波與下行P-SV波的極化方向相同，因而利用直達P波計算極化角，通過極化濾波分離上行P-P波與上行P-SV波。這種做法有很大缺陷：一是時直達P波與上行P-SV波極化方向不同；二是隨著地層深度的增加，各個波場的極化角是變化的，利用單一的直達P波極化角對整道數據進行極化濾波顯然不合理，因此必須對極化角進行時變估計。

首先將地震記錄r變換到時頻域，r的STFT可以表示為

式中：h為窗函數；t為時間；ω為角頻率。由于STFT是線性變換，式（6）可以寫為

在時頻域對估計極化角做如下假設。

（1）在以某一時刻tj為中心的時頻窗（tj，ωk）內僅有單一波場Si發育，即

式中：Si是時間域源信號s i的STFT；j為時間采樣點；k為頻率采樣點。顯然存在

式中：Rz（tj，ωk）、Rx（tj，ωk）分別為z分量和x分量的STFT變換；αi為各單一波場的極化角。

（2）在以t j為中心開始的連續K個時頻窗內，各單一波場的極化角保持不變，即

計算K個時窗內的極化角均值與方差，即

如果在K個時頻窗內存在var（αi）=0，則說明此時僅存在單一波場Si，且其極化角為（t j，ωk）。逐一計算掃描整個時頻域數據，得到所有的單源時頻點及相應的極化角簇（l=1，2，…，L）。采用K-means聚類算法，將按照波場類型聚成4類φ=｛φ1，φ2，φ3，φ4｝，假設每類中的元素數為H=｛H1，H2，H3，H4｝，類心為，同一類之間的緊密程度定義為

3.2 單一波場的估計

對于欠定盲源分離，即使已知混合矩陣，源信號依然是多解的。稀疏表示的最終目標是得到最稀疏的s，也就是說使s中盡可能多的元素取值為零，數學上可以歸結為L0范數優化問題

式中：為s的估計；為L0范數稀疏測度。雖然L0范數最能體現稀疏性的含義，但是求解對噪聲太敏感、魯棒性差，因而難以有效實現，所以Donoho等［19］建議采用L1范數稀疏測度代替L0范數，即

單一波場的估計仍在時頻域進行，按照稀疏性理想程度分以下幾種情形進行分析。

（1）如果在時頻域單一波場源信號足夠稀疏，即在每一個以t j為中心的時頻窗（t j，ωk）內僅有唯一的單一波場Si發育，則最小L1范數解能夠非常好地逼近單一波場源信號，有

式中αi（t j，ωk）為當前時頻窗內單一波場Si的極化角。

（2）如果在每一個以tj為中心的時頻窗（tj，ωk）內僅有兩種波場發育，此時可以用最短路徑分解法對波場進行估算。圖2a為最短路徑分解法示意圖，在Rx-R z平面中，O為原點，R為觀測信號在時頻域中的某一點，其矢量方向與R x軸的夾角為θ，即當前時刻的極化角可由式（10）計算。a1和a2為離θ最近的兩種波場的基矢量，也是極化矩陣的列向量。沿著a1和a2的方向做平行四邊形，分解向量R，則有

（3）如果在每一個以t j為中心的時頻窗（tj，ωk）內不少于三種波場發育，由于多源干涉，最短路徑分解法對源信號的估計效果會較差。下面以三種單一波場源信號為例，提出一種改進的最短路徑分解法。

如圖2b所示，a1和a2仍為離θ最近的兩種波場的基矢量，a3為離θ較遠的波場基矢量，有

圖2 最短路徑分解法及其改進方法示意圖

首先通過以上方法估計離tj最近時刻t′j的單源頻率點。由VSP真振幅恢復的補償公式可推導出與的關系

式中n為補償因子，可利用tj與t′j段的時間—能量對進行最小平方擬合獲得。

根據以上假設條件，基矢量a3在一定時間段內保持不變，即

將向量R與a3方向的向量作矢量相減，有

然后再利用最短路徑法分解R′，估計出與，就完成了全部單一波場的估計。

綜上所述，基于稀疏表示的VSP波場盲分離包括以下幾個步驟：

（1）對原始的VSP水平分量做檢波器定位，得到定位后的x分量；

（2）對定位后的x分量與原始z分量做STFT，將信號變換到時頻域；

（3）將時頻窗按時間分段，估算各時間段內單一波場的極化角；

（4）在時頻域內應用改進的最短路徑分解法，估算各時間段內的單一波場；

（5）對時頻域的單一波場做STFT反變換，得到時域的單一波場。

4 數值計算

4.1 模型試算

采用如圖3a所示的三層水平模型，檢波器起始深度為500m，共100級，級間距為20m，井源距為1000m。圖3b、圖3c分別為正演模擬的非零井源距VSP水平分量記錄與垂直分量記錄，可以看出，記錄中包含上、下行P-P波與上、下行P-SV波多種類型的波場。由于極化方向不同，各個波場在水平分量與垂直分量記錄上的強弱關系也不同。

圖3 模型示意圖及正演模擬記錄

圖4為模型數據波場分離結果，可以看出，f-k濾波分離的上行P-P波中仍含有下行波的痕跡（圖4a中橢圓內），這是由于二維傅里葉變換會產生空間假頻問題。中值濾波分離的上行P-P波中仍存在下行波、上行P-SV波的痕跡（圖4b中橢圓內）。一方面中值濾波對排齊要求較高，而對于非零井源距VSP，不同深度界面反射的波場無法同時排齊，造成濾波后仍存在殘余波場能量；另一方面中值濾波存在邊界效應，在邊界處無法完全分離單一波場。由于在時頻域單一波場的稀疏性較好，通過本文方法分離的上行P-P波質量較高，既沒有二維濾波帶來的波場畸變，又考慮了波場的偏振特征，因而具有較好的保幅性。

4.2 實際資料處理

圖5是實際非零井源距VSP資料。由于原始地震記錄上存在噪聲，而本文討論的數學模型是無噪的，因而先對記錄進行隨機噪聲衰減，再進行波場分離。圖6為實際資料波場分離結果。可以看出：f-k濾波分離的波場存在混波現象（圖6a中橢圓內）；中值濾波分離的波場信噪比較低，存在明顯的下行波殘余能量（圖6b中橢圓內），且同相軸連續性較差；本文方法分離的波場同相軸清晰，波場連續，分離效果優于f-k濾波與中值濾波。

5 結論

本文介紹了盲源信號分離的基本原理，在此基礎上提出了一種基于稀疏表示的VSP波場分離方法。該方法首先將VSP地震信號變換到時頻域，然后基于稀疏性假設對各個單一波場的極化角進行估計，進而估算極化矩陣，再根據各個單一波場在時頻域的極化特征，采用改進的最短路徑法對單一波場進行分離。由于極化矩陣在線性變換下保持不變，該方法能夠得到保真性較好的上、下行P-P波、PSV波，處理效果優于傳統的f-k濾波與中值濾波。模型試算與實際資料處理效果驗證了本文方法的有效性。

圖4 模型數據波場分離結果

圖5 檢波器定位后的實際非零井源距VSP資料

［1］ Moon W，Carswell A，Tang R et al.Radon transform wave field separation for vertical seismic profiling data.Geophysics，1986，51（4）：940-947.

［2］ Esmersoy C.Inverse of P and SV waves from multicomponent offset vertical seismic profiles.Geophysics，1990，55（1）：39-50.

［3］ Cho W H and Spencer T W.Estimation of polarization and slowness in mixed wavefields.Geophysics，1992，57（6）：805-814.

［4］ 魏修成，吳律，田子奇.變速視慢度波場分離.石油地球物理勘探，1993，28（4）：418-429.Wei Xiucheng，Wu Lüand Tian Ziqi.Wavefield separation using apparent slowness of variable velocity.OGP，1993，28（4）：418-429.

［5］ Anderson S，Nehorai A.Analysis of polarized seismic wave model.IEEE Transactions on Signal Processing，1996，44（2）：379-386.

［6］ Wang Y B，Singh SC and Barton P J.Separation of P-and SV-wavefields from multi-component seismic data in theτ-Pdomain.Geophysical Journal International，2002，151（2）：663-672.

［7］ Blias E.VSP wavefield separation：Wave-by-wave optimization approach.Geophysics，2007，72（4）：47-55.

［8］ 王成禮.VSP不同觀測方式的波場特征與波場分離方法.石油地球物理勘探，1990，25（2）：148-160.Wang Chengli.Various VSP wave fields generated by different observation ways and the corresponding wave field separation methods.OGP，1990，25（2）：148-160.

［9］ 杜婧，王尚旭，劉國昌等.基于局部斜率屬性的VSP波場分離研究.地球物理學報，2009，52（7）：1867-1872.Du Jing，Wang Shangxu，Liu Guochang et al.VSP-wavefield separation using local slopes attribute.Chinese Journal of Geophysics，2009，52（7）：1867-1872.

［10］ 崔汝國，牟風明，宋維琪等.VSP浮動坐標系偏振濾波.石油地球物理勘探，2010，45（1）：10-14.Cui Ruguo，Mu Fengming，Song Weiqi et al.Polarization filtering with drift coordinates for VSP data.OGP，2010，45（1）：10-14.

［11］ 孫文博，孫贊東.3D3C VSP資料處理中矢量波場分離方法的研究.石油地球物理勘探，2009，44（6）：708-719.Sun Wenbo，Sun Zandong.Studies on vector wave field separation method for 3D3C VSP data.OGP，2009，44（6）：708-719.

［12］ 馬志霞，孫贊東，白海軍等.三維三分量VSP多種波場分離方法對比.石油地球物理勘探，2010，45（2）：219-224.Ma Zhixia，Sun Zandong，Bai Haijun et al.Comparative studies on 3D3C VSP multiple wave field separation methods.OGP，2010，45（2）：219-224.

［13］ Jiang X，Lin J，Ye F et al.Separation of P-P and P-SV wavefields by high resolution parabolic Radon transform.Journal of Applied Geophysics，2015，119（8）：192-201.

［14］ Moni A，Bean C J，Lokmer I et al.Source separation on seismic data：Application in a geophysical setting.IEEE on Signal Processing，2012，29（3）：16-28.

［15］ Van der Baan M.PP／PS wavefield separation by independent component analysis.Geophysical Journal International，2006，166（1）：339-348.

［16］ Vrabie V D，Bihan N L，Mars J I.Multicomponent wave separation using HOSVD／unimodal-ICA subspace method.Geophysics，2006，71（5）：V133-V143.

［17］ Liu K H，Dragoset W H.Blind-source separation of seismic signals based on information maximization.Geophysics，2013，78（4）：V119-V130.

［18］ Bofill P，Zibulevsky M.Underdetermined blind source separation using sparse representations.Signal Processing，2001，81（11）：2353-2362.

［19］ Donoho D L，Elad M.Maximal sparsity representation viat1minimization.Proceedings of National Academy Science，2003，100（3）：2197-2202.

P631

A

10.13810／j.cnki.issn.1000-7210.2017.01.005

趙彥青，蕭蘊詩.基于稀疏表示的非零井源距VSP波場分離方法.石油地球物理勘探，2017，52（1）：27-33，55.

1000-7210（2017）01-0027-07

*上海市嘉定區曹安公路4800號同濟大學電子與信息工程學院，201804。Email：forrest_zyq＠163.com

本文于2016年3月28日收到，最終修改稿于同年11月13日收到。

本項研究受國家自然科學基金項目（40872090）資助。

（本文編輯：宜明理）

趙彥青 博士研究生，1981年生；2002年畢業于青島大學，獲工業自動化專業學士學位；2006年畢業于同濟大學，獲控制理論與控制工程專業碩士學位；現為同濟大學控制理論與控制工程專業博士研究生，主要從事盲信號和地震信號處理及其方法研究。