基于PCA的井間地震偏振分析方法及其應(yīng)用

金 聰 林 松 鄧小虎 程 飛 程 邈

1 地震預(yù)警湖北省重點(diǎn)實(shí)驗(yàn)室，武漢市洪山側(cè)路40號(hào)，4300712 湖北省地震局，武漢市洪山側(cè)路48號(hào)，4300713 武漢地震工程研究院有限公司, 武漢市洪山側(cè)路40號(hào)，4300714 中國地質(zhì)大學(xué)(武漢)海洋地質(zhì)資源湖北省重點(diǎn)實(shí)驗(yàn)室，武漢市魯磨路388號(hào)，430074

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

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

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

1 三分量地震偏振分析

1.1 PCA方法求解偏振角

某一波場延續(xù)時(shí)間T內(nèi)，質(zhì)點(diǎn)振動(dòng)的3個(gè)分量X、Y、Z可用矩陣來描述[10]：

(1)

式中，N=T/Δt+1為時(shí)窗內(nèi)采樣點(diǎn)數(shù)，C為3×N矩陣，規(guī)定其中各分量在延續(xù)時(shí)間內(nèi)的均值為0，則協(xié)方差矩陣可通過式(2)計(jì)算：

(2)

協(xié)方差矩陣S為3×3實(shí)對(duì)稱矩陣，矩陣中元素為三分量振動(dòng)的自相關(guān)和互相關(guān)數(shù)值，即各分量之間的協(xié)方差。矩陣S存在3個(gè)非負(fù)實(shí)特征值λ1、λ2、λ3，且滿足λ1≥λ2≥λ3，令特征向量為u1、u2、u3，則存在如下關(guān)系式：

(S-λ2I)u=0

(3)

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

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

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

(4)

1.2 偏振旋轉(zhuǎn)

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

H0=Xcosθ+Ysinθ

T=Ycosθ-Xsinθ

(5)

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

圖1 坐標(biāo)旋轉(zhuǎn)示意圖Fig.1 Schematic diagram of coordinate rotation

HP=H0cosα+Zsinα

R=Zcosα-H0sinα

(6)

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

2 理論模型試算

2.1 正演模擬

為說明利用協(xié)方差矩陣求取偏振角的有效性，構(gòu)建三維均勻?qū)訝罱橘|(zhì)模型進(jìn)行理論試算。模型大小為150 m×150 m×300 m，網(wǎng)格間距為1.0 m，物性參數(shù)見表1。圖2(a)為模型XOZ平面示意圖，井位分布見圖2(b)，激發(fā)點(diǎn)設(shè)置在井1，深度在150 m處，接收井設(shè)置在井3，全井接收。

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

表1 層狀理論模型物性參數(shù)

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

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

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

2.2 理論記錄偏振分析

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

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

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

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

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

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

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

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

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

3 應(yīng)用實(shí)例

3.1 現(xiàn)場測試

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

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

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

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

(7)

由式(7)可知，影響波速計(jì)算的2個(gè)參數(shù)分別為距離L和初至走時(shí)t?？紤]各鉆孔因各種外界因素可能會(huì)存在偏差和傾斜，導(dǎo)致其不嚴(yán)格與地面垂直或互相平行，因此需要同時(shí)進(jìn)行井斜測試，以便對(duì)激發(fā)點(diǎn)和接收點(diǎn)之間的距離L進(jìn)行校正。參數(shù)t的拾取與原始記錄質(zhì)量密切相關(guān)。

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

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

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

3.2 實(shí)測資料偏振分析

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

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

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

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

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

圖12 實(shí)測資料旋轉(zhuǎn)后三分量記錄Fig.12 Three-component records after rotation of measured data

4 結(jié) 語

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

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

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

3)該方法對(duì)初至P波的信噪比依賴較高，在實(shí)際應(yīng)用中如果未能獲取信號(hào)質(zhì)量較好的初至P波，可能無法得到理想的偏振旋轉(zhuǎn)角，后續(xù)工作將在該方法的抗噪能力上作進(jìn)一步研究。