基于PCA的井間地震偏振分析方法及其應用

金 聰 林 松 鄧小虎 程 飛 程 邈

1 地震預警湖北省重點實驗室，武漢市洪山側路40號，4300712 湖北省地震局，武漢市洪山側路48號，4300713 武漢地震工程研究院有限公司, 武漢市洪山側路40號，4300714 中國地質大學(武漢)海洋地質資源湖北省重點實驗室，武漢市魯磨路388號，430074

井間地震是將震源與檢波器放置在相鄰的兩口井中，在一口井中激發地震，在另一井中接收信號。與常規單分量相比，井間三分量檢波器可采集到包含豐富運動學和動力學信息的矢量波場[1]。但在標準正交坐標系中，各分量極性和幅值會隨檢波器放置狀態(包括方位角和傾角)的變化而變化。由于井下檢波器擺放狀態完全隨機，地震響應水平分量隨機分布在X、Y分量各地震道上，因井斜程度不同和介質各向異性，Z分量各地震道上的響應幅值也會存在不同程度的損失，這對后期利用三分量資料開展巖性研究造成一定干擾[2]。因此，在預處理中需要對三分量檢波器進行方向校正，在VSP中稱之為檢波器重定向，是井間地震數據處理的關鍵[3]。

目前設計的全數字三分量檢波器利用傾斜度和重力測試可實現垂直分量傾斜校正，利用姿態檢測技術可自動進行水平分量數據校正[4]，但在實際生產中還未大規模應用?，F階段國內外主要利用直達P波偏振特征來估算檢波器的方位角與入射傾角[5-6]。三分量地震記錄的偏振分析一般是在時間域實現，常被應用于隨機噪聲壓制、波場識別與分離、介質屬性分析等方面[7]。具體分析方法有兩大類：一是基于最大能量準則分析法，在時窗內計算不同分量中波場的振幅或能量，使之達到最大時的角度即為偏振角[8]；二是基于主成分分析(PCA)法，在時窗內構造協方差矩陣，對該矩陣的特征值和特征向量進行求解，從而獲得相應的偏振信息。相比較而言，后者原理簡單，實現方便，所求偏振角精度較高[9]。

目前，在跨孔井間彈性波波速測試中，關于波形識別的研究較少，如果能將三分量地震偏振分析應用到初至波形起跳的識別中，可以大幅提升在復雜地層中采集的原始資料的質量。本文基于PCA方法進行偏振分析，構建層狀介質模型進行理論試算，分析利用協方差矩陣求取偏振角的有效性，并將該方法應用于跨孔井間彈性波波速測試，通過偏振旋轉可使P波和S波能量得到不同程度的加強，使初至起跳的拾取更加高效和準確。

1 三分量地震偏振分析

1.1 PCA方法求解偏振角

某一波場延續時間T內，質點振動的3個分量X、Y、Z可用矩陣來描述[10]：

(1)

式中，N=T/Δt+1為時窗內采樣點數，C為3×N矩陣，規定其中各分量在延續時間內的均值為0，則協方差矩陣可通過式(2)計算：

(2)

協方差矩陣S為3×3實對稱矩陣，矩陣中元素為三分量振動的自相關和互相關數值，即各分量之間的協方差。矩陣S存在3個非負實特征值λ1、λ2、λ3，且滿足λ1≥λ2≥λ3，令特征向量為u1、u2、u3，則存在如下關系式：

(S-λ2I)u=0

(3)

式中，I為3×3單位矩陣，0為元素值均為0的列向量。

橢球極化參數由特征向量u和特征值λ來表示，特征向量表示3個互相垂直的極化方向，而特征值為地震相平均能量的度量。如果信號為線性偏振(如體波)，在理想情況下協方差矩陣僅有一個非零特征值，此時存在λ1>>λ2、λ3，λ1對應的特征向量代表極化方向；在波場互相干擾或存在噪聲的情況下，可得到3個特征值，此時為橢球偏振，對應的3個特征向量可確定3個軸的取向，其中最大特征值對應的特征向量代表主偏振方向[11]。

主特征向量u1歸一化后的3個標量元素(u11,u21,u31)分別為在直角坐標系Z、X、Y方向的位移標量，因此可求解水平方向方位角θ和垂直方向傾角α，表達式分別為[12]：

(4)

1.2 偏振旋轉

在求得水平方向和垂直方向旋轉角后，對采集到的原始三分量記錄進行坐標旋轉。首先在水平面內進行旋轉，如圖1(a)所示，水平分量H0與X分量正向的夾角為θ，根據幾何關系可知，變換公式為：

H0=Xcosθ+Ysinθ

T=Ycosθ-Xsinθ

(5)

旋轉后P波和SV波的水平分量主要存在于H0分量中，而T分量中主要為SH波能量。圖1(b)為垂直方向旋轉，變換公式為：

圖1 坐標旋轉示意圖Fig.1 Schematic diagram of coordinate rotation

HP=H0cosα+Zsinα

R=Zcosα-H0sinα

(6)

通過坐標旋轉后，P波能量主要集中在HP分量中，SV波和SH波能量分別集中在R分量和T分量中，有利于后期資料處理。

2 理論模型試算

2.1 正演模擬

為說明利用協方差矩陣求取偏振角的有效性，構建三維均勻層狀介質模型進行理論試算。模型大小為150 m×150 m×300 m，網格間距為1.0 m，物性參數見表1。圖2(a)為模型XOZ平面示意圖，井位分布見圖2(b)，激發點設置在井1，深度在150 m處，接收井設置在井3，全井接收。

圖2 激發和接收井位分布Fig.2 Distribution of shooting and receiving wells

表1 層狀理論模型物性參數

采用彈性介質交錯網格有限差分法進行數值模擬，選擇時間2階、空間6階精度。具體參數為：時間采樣間隔為0.1 ms，道間距為1.0 m，記錄長度為300 ms，雷克子波主頻為120 Hz，震源類型為脹縮力震源。

通過三維數值模擬得到理論三分量地震記錄(圖3)。在均勻層狀介質模型中，除初至縱波(P)、透射縱波(Tp)和上下界面反射縱波(Rp)外，還存在上下界面處的轉換橫波(Rps、Tps)。X分量上P波能量較強，特別是在震源深度附近幅值最大，也有S波分布；Y分量上P波和S波能量均較弱；Z分量上有S波分布。在界面所在深度100 m和200 m附近存在多種波形相互疊加。

圖3 層狀模型理論三分量記錄Fig.3 Three-component records of theoreticallayered model

2.2 理論記錄偏振分析

選取第4道進行偏振分析，三分量波形見圖4。由激發點與接收點的相對位置和Snell定律可知，水平方向方位角θ=18.434 9°，垂直方向傾角α=-47.922 2°。

圖4 井3接收第4道三分量波形Fig.4 Three-component waveform of channel 4 of well-3

通過PCA方法計算得到θ=18.369 6°，α=-47.681 0°，與理論值偏差分別為0.065 3°和0.241 2°。如果考慮到原始記錄因網格化數值模擬存在的系統誤差，計算得到的θ和α精度較高。

利用式(5)和式(6)進行坐標旋轉，得到HP分量、R分量和T分量波形圖見圖5。與圖4原始三分量波形相比，坐標旋轉后各波形運動學特征無改變，R分量中主要包含S波成分，P波能量較弱；HP分量中P波成分得到加強，S波能量被壓制；T分量中能量幾乎為0，由于該分量方向與射線平面垂直，理論上無SV波存在，試算結果與理論一致。

圖5 井3接收第4道旋轉后三分量波形Fig.5 Three-component waveform after rotation ofchannel 4 of well-3

將圖3整個三分量記錄進行偏振分析，通過PCA方法計算得到垂直方向傾角(圖6)。圖中藍色實線為計算結果，紅色實線為Snell定律計算的理論結果。通過對比可知，對于無波形疊加的記錄道，計算的偏振角與理論值吻合較好；在界面分層深度附近，因多種波形疊加在直達P波上，影響了P波偏振橢圓的極化特征，導致偏振角計算出現誤差，但與理論值相差較小。

圖6 層狀模型垂直方向偏振曲線Fig.6 Vertical polarization curve of layered model

利用式(5)和式(6)進行坐標旋轉得到HP分量和R分量記錄(圖7)?？梢钥闯?，HP分量中P波能量得到加強，S波能量被壓制；R分量中主要為S波成分，幾乎無法看到直達P波，僅在100～200 m深度范圍內存在部分界面反射P波能量，這是由于反射P波偏振方向與該深度處S波偏振方向之間夾角較小，垂直方向坐標旋轉會將部分反射P波能量旋轉到R分量上。

圖7 井3偏振后三分量記錄Fig.7 Three-component records after rotation of well-3

3 應用實例

3.1 現場測試

跨孔井間彈性波波速測試一般采用單孔激發、單孔或兩孔接收的方式，孔中由三分量傳感器接收信號。該方法是一種快速、準確的原位測試技術，能提供高分辨率的P波和S波波速曲線，并可根據波速計算動彈性模量、動剪切模量、動泊松比等動力學參數，為巖土工程設計提供依據。

本次野外數據采集采用IPG5000脈沖發生器和BIS-SH震源探頭，接收采用BGK3信號檢波器。圖8為實際工作示意圖，將震源下放至鉆孔中預定深度，使用充氣氣囊機制與井壁耦合，震源方向可通過地面剛性抗扭管進行調節，將三分量檢波器下放至另一鉆孔中相同深度并緊貼孔壁，激發震源，記錄相應波形圖；檢查無誤后將震源旋轉180°，反向激發并再次記錄相應波形。從上至下(或從下至上)依次按上述步驟進行采集，即可完成整個鉆孔波速測試。

圖8 跨孔波速測試工作示意圖Fig.8 Schematic diagram of cross-hole wave velocity test

通過改變震源的極性方向，結合正向激發和反向激發的三分量記錄，從中選取振幅能量較大、初至較明顯的分量，分別讀取P波初至與S波初至，根據孔距即可得到地下巖土層P波和S波速度，波速計算公式為：

(7)

由式(7)可知，影響波速計算的2個參數分別為距離L和初至走時t。考慮各鉆孔因各種外界因素可能會存在偏差和傾斜，導致其不嚴格與地面垂直或互相平行，因此需要同時進行井斜測試，以便對激發點和接收點之間的距離L進行校正。參數t的拾取與原始記錄質量密切相關。

由于測試采用水平同步方式進行，理論上激發點與接收點位于同一水平面，Z分量與旋轉后R分量一致，垂直方向傾角為0°，S波主要分布在Z分量上。但在實際操作中，震源與接收探頭下放深度往往存在人為誤差、接收探頭緊貼孔壁時存在偏差、孔壁可能存在傾斜等因素，導致采集的三分量波形中P波和S波能量會根據鉆孔傾斜程度在3個分量上重新分布，不利于初至波拾取。

本文以某港口地層跨孔彈性波波速測試為例，選取其中一對鉆孔采集的實測資料，圖9為其三分量記錄，圖中(a)～(c)為正向激發，(d)～(f)為反向激發。可以看出，P波初至起跳較為明顯，而S波由于地層巖性不同，能量衰減程度存在差異。此時如果鉆孔傾斜或激發點和接收點不在同一水平面，會導致S波能量分散投影到3個分量上，從而影響初至S波的拾取。

圖9 實測資料三分量記錄Fig.9 Three-component records of measured data

3.2 實測資料偏振分析

圖10為第12道三分量波形圖，其中Rx1、Ry1、Rz1分別為正向激發X、Y、Z分量，Rx2、Ry2、Rz2分別為反向激發X、Y、Z分量。

圖10 實測資料第12道三分量波形Fig.10 Three-component waveform of measured data of channel 12

利用初至P波構建協方差矩陣，求取三分量波形的水平方位角和垂直傾角，通過式(5)和式(6)進行坐標旋轉，得到HP分量和R分量(圖11，圖中黑色實線為正向分量，藍色實線為反向分量)。通過與圖10對比可知，HP分量中P波能量得到加強，振幅增大到原始記錄的2倍左右；S波能量主要集中在R分量中，與原始記錄S波振幅相比，同樣存在小幅度增加。通過偏振旋轉可將正向分量與反向分量進行重疊，S波初至更加明顯，從而容易識別和讀取，得到的橫波速度也更加精確。

圖11 實測資料第12道旋轉后三分量波形Fig.11 Three-component waveform after rotation ofmeasured data of channel 12

對每一道三分量記錄重復上述步驟，得到經偏振旋轉后的HP分量和R分量記錄(圖12)，P波初至在HP分量中讀取，S波初至在R分量中讀取。相比原始三分量記錄而言，可減少篩選具有較明顯初至分量的步驟，而且P波和S波能量均存在不同程度的加強，對初至拾取更加有利。從圖中可以明顯看出，深度10 m附近和18～26 m處存在明顯的橫波低速帶。

圖12 實測資料旋轉后三分量記錄Fig.12 Three-component records after rotation of measured data

4 結 語

本文介紹了一種基于PCA方法的井間地震三分量記錄偏振分析方法，通過構建理論層狀介質模型，詳細論證了該方法在偏振分析中的有效性，并將該方法應用于跨孔井間波速測試。通過實測資料計算，獲得以下結論：

1)利用PCA方法求取三分量記錄中直達P波的偏振狀態方便、有效，在無其他波形混疊的情況下，得到的偏振旋轉角精度較高；但在界面深度處存在波形疊加時，會影響初至P波偏振橢圓的極化特征，從而導致偏振旋轉角計算出現誤差，不過該誤差較小。

2)跨孔井間波速測試采集的三分量資料經偏振旋轉后，在HP分量中讀取P波初至，在R分量中讀取S波初至，與原始三分量記錄相比，可減少篩選具有較明顯初至分量的步驟，且P波和S波能量均存在不同程度的加強，可在復雜地層中更加準確地拾取初至起跳，為巖土工程設計提供更加精確的動力學參數。

3)該方法對初至P波的信噪比依賴較高，在實際應用中如果未能獲取信號質量較好的初至P波，可能無法得到理想的偏振旋轉角，后續工作將在該方法的抗噪能力上作進一步研究。