陳占國 曾昭翰 楊心超

(中國石化石油物探技術(shù)研究院，江蘇南京 211103)

0 引言

隨著地震各向異性理論的發(fā)展和應(yīng)用，工業(yè)界基于VTI介質(zhì)的地震建模和成像技術(shù)日趨成熟； TTI各向異性建模和成像技術(shù)在實際中得到較廣泛的應(yīng)用[1]； 學術(shù)界已開始研究如基于正交各向異性理論等的更復(fù)雜應(yīng)用技術(shù)。

地震各向異性參數(shù)的獲取一直是各向異性應(yīng)用的關(guān)鍵。與地面地震相比，VSP測量作為各向異性參數(shù)獲取的一種途徑，具有較高的精度和準確的深度優(yōu)勢。前人研究表明，應(yīng)用 VTI模型可以更加精確地描述長排列上的地震波旅行時特性[2]； 基于 VTI假設(shè)，給出了不同形式的時差關(guān)系：Byun等[3-4]應(yīng)用一種雙參數(shù)速度分析對VSP合成數(shù)據(jù)進行速度分析； Hake等[5]導出了VTI介質(zhì)時差的3項泰勒展開式； Tsvankin等[6-7]對Hake的公式進行了改進，并引入漸近擬合使新公式在大炮檢距上更逼近真值； Alkhalifah等[8-9]提出用兩個參數(shù)表達的非雙曲時差關(guān)系描述 VTI介質(zhì)的地震波旅行時，已在地面地震數(shù)據(jù)各向異性動校正等處理中得到了廣泛的應(yīng)用[10-12]； 白海軍等[13]和陳沅忠等[14]將Alkhalifah的方法應(yīng)用于VSP數(shù)據(jù)，均取得了較好的效果。

以上方法主要是基于VTI或TTI介質(zhì)的縱波非雙曲線時差進行速度分析和旅行時反演，難以獲得更復(fù)雜的任意各向異性參數(shù)。VSP能夠直接測量地震波的慢度和極化信息，且不受上覆地層影響。極化和慢度矢量是表征Christoffel方程的直接參數(shù)，理論上只要獲得足夠多的慢度和極化信息，就可以反演介質(zhì)的所有各向異性參數(shù)。根據(jù)這一特點，Jech等[15]提出了各向異性介質(zhì)的一階擾動方法。Horne 等[16]和Gomes等[17]基于擾動理論，提出了應(yīng)用慢度和極化矢量求解介質(zhì)弱各向異性參數(shù)的方法，并得到了一些成功的應(yīng)用。

在前人研究的基礎(chǔ)上，本文系統(tǒng)地介紹了基于慢度和極化信息反演各向異性參數(shù)的方法原理，并給出了Walkaway VSP慢度和極化數(shù)據(jù)的獲取方法，提出了基于加權(quán)迭代的任意弱各向異性參數(shù)反演算法的詳細實現(xiàn)過程，應(yīng)用一組TTI模型合成數(shù)據(jù)和一口實際井資料進行處理試驗，取得了理想的效果。

1 方法原理

1.1 Christoffel 方程及擾動方法

介質(zhì)中傳播的平面波遵循Christoffel方程[18]

(1)

式中：Γjk=aijklninl為Christoffel矩陣，其中a是密度歸一化彈性張量，n為波的傳播方向矢量，i、j、k、l均取1～3，表示坐標方向；m=1、2、3，分別代表各向異性介質(zhì)中傳播的qS1、qS2、qP波；Gm是Christoffel矩陣的三個特征值；g(m)為三種波的極化矢量。式(1)中使用了愛因斯坦求和約定。

由于實際介質(zhì)大多是弱各向異性[19-20]，基于擾動理論，Christoffel傳播方程可近似表述為[21]

(2)

式中：線“～”表示在各向同性參考介質(zhì)中的基本量； “Δ”表示在弱各向異性介質(zhì)中的擾動量。

1.2 qP波慢度和極化矢量正反演公式

只考慮式(2)中的qP波(m=3)，忽略其中二階項，并引入了三個相互垂直的單位矢量[22]

(3)

(4)

(5)

式中：p和g分別是qP波的慢度和極化矢量；α和β分別表示參考介質(zhì)中qP波和qS波的相速度，K′=1、2；B為與弱各向異性參數(shù)相關(guān)的弱各向異性矩陣[25]，可表示為

(6)

其中Δa為密度歸一化彈性張量的擾動量。

反過來，通過慢度矢量和極化矢量，可以求解弱各向異性介質(zhì)的彈性參數(shù)[23-24]

(7)

B33=-2α4piΔpi

(8)

(9)

式(7)和式(8)是關(guān)于15個弱各向異性參數(shù)的線性公式。

1.3 qP波慢度矢量和極化矢量的計算

1.3.1 極化矢量

在數(shù)學上，質(zhì)點的運動軌跡可以認為是空間點的集合。對于三分量的VSP記錄，這種離散采樣點的質(zhì)點運動軌跡具有相關(guān)性，可以通過三分量記錄的協(xié)方差矩陣統(tǒng)計其空間運動規(guī)律。根據(jù)矩陣的特征值和特征向量可定量地確定極化橢球(長、中、短)各軸的大小和橢球的空間位置及其標準橢球方程，并能求出質(zhì)點運動軌跡的極化參數(shù)。協(xié)方差矩陣定義為

(10)

求解協(xié)方差矩陣的特征值和對應(yīng)的特征向量，可以確定極化方向，其中最大特征值對應(yīng)的特征向量就是qP波的極化矢量。

1.3.2 慢度矢量

一般情況下，可通過鄰近兩個檢波器的時差dt和深度差dh求取慢度的垂向分量

(11)

慢度p是一個局部量，不依賴于上覆地層速度結(jié)構(gòu)。根據(jù)地震波走時的互易性，如果介質(zhì)是橫向均勻的，則炮點到接收點的走時等于接收點到炮點的走時。當介質(zhì)橫向非均勻性比較弱時可以將接收點與炮點互易，利用共接收點的旅行時計算慢度矢量的水平分量。

旅行時曲線一般情況下都不光滑，其一階導數(shù)不連續(xù)甚至不存在。為了得到更穩(wěn)定的慢度矢量，本文采用多點平滑和相鄰五點三次多項式求導的方法計算慢度矢量。

1.4 Walkaway VSP反演弱各向異性參數(shù)方法

通過將BKL代入式(7)和式(8)，可以建立弱各向異性參數(shù)與慢度及極化矢量的關(guān)系式。如果慢度矢量的擾動Δp的三個分量和極化矢量g(3)的三個分量都已知，可以單獨使用式(7)或式(8)建立一個關(guān)于未知數(shù)為弱各向異性參數(shù)的線性方程組，或者聯(lián)合式(7)和式(8)建立方程組，求解該線性方程組即可得弱各向異性參數(shù)。

(12)

式中：N是數(shù)據(jù)點個數(shù)；y表示式(7)和式(8)的觀測參數(shù)和已知量；w為觀測參數(shù)的權(quán)系數(shù)；x為待反演的弱各向異性參數(shù)，其維數(shù)為Npar；J=1,2,…,Npar；A為N×Npar階系數(shù)矩陣。弱各向異性參數(shù)的估計可以通過奇異值分解求解，即

x=VΛ-1UTd

(13)

式中：U和V是系數(shù)矩陣ATwA兩個正交的N×N奇異矩陣；Λ是由ATwA的奇異值構(gòu)成N×N對角矩陣；d=ATwy。

權(quán)系數(shù)w由數(shù)據(jù)的初始權(quán)重wd和更新后的測量值與正演計算值殘差的權(quán)重wr組成

(14)

數(shù)據(jù)的初始權(quán)重wd由兩部分構(gòu)成，一部分是參與運算的數(shù)據(jù)類型的質(zhì)量估計權(quán)系數(shù)wq，另一部分權(quán)系數(shù)由慢度或極化數(shù)據(jù)的可靠性決定，即

(15)

權(quán)系數(shù)wr為反演迭代時慢度或者極化矢量殘差的再加權(quán)[26]

(16)

式中：drI為第I個慢度或者極化矢量的觀測值與正演計算值的殘差；σmed=0.67449，為高斯噪聲的中值；ω是dr的標準偏差截斷因子，一般選擇3～6間的數(shù)； drmed定義為

drmed=med(|drI-med(drI)|)

(17)

其中med(?)表示取數(shù)列的中值。

求解式(13)后，每個模型參數(shù)的最小二乘估計誤差EJ為[26]

(18)

式中：TJJ為奇異值和奇異矩陣構(gòu)成的協(xié)方差矩陣T=VΛ-2VT的對角元素；Ω為加權(quán)剩余量的變化

(19)

獨立的弱各向異性參數(shù)共有21個，其中與qP波相關(guān)的有15個：εx、εy、εz、δx、δy、δz、χx、χy、χz、ε15、ε16、ε24、ε26、ε34和ε35。若要反演15個參數(shù)，需要布設(shè)5條以上Walkaway VSP測線； 如果對兩個正交的測線進行聯(lián)合反演，最多可得到9個參數(shù)：εx、εy、εz、δx、δy、ε15、ε24、ε34和ε35； 如果僅使用一條測線，只能反演5個參數(shù)，即εx、εz、δx、ε15和ε35[25]。

2 模型資料測試

模型是一個彈性參數(shù)隨深度線性變化的TTI模型，相速度的各向異性程度(相速度的最大最小值之差與最大值的比值)為4%。地表處的TTI彈性參數(shù)(密度歸一化的)矩陣Θ3和10km深處的彈性參數(shù)矩陣Θ4分別由相應(yīng)的兩個VTI彈性矩陣Θ1和Θ2經(jīng)過繞y軸順時針旋轉(zhuǎn)40°、再繞z軸(向下為正)順時針旋轉(zhuǎn)30°后得到。四個彈性參數(shù)矩陣分別為

(20)

(21)

(22)

(23)

觀測系統(tǒng)由6條等方位間隔的Walkaway VSP測線組成。每條測線布設(shè)32個炮點，井孔兩邊各16個； 井源距范圍是100～3100m，炮間距為200m。井中1000～1550m深度等間距布設(shè)12個檢波器，檢波器間距為50m。

聯(lián)合應(yīng)用式(7)和式(8)，通過6條Walkaway VSP的慢度矢量和極化矢量的3個分量，反演了模型的15個弱各向異性參數(shù)。

以深度1250m的檢波點為例，其他深度結(jié)果類似。圖1展示了方位角為30°/210°、90°/270°和150°/330°測線測量的慢度和極化矢量隨井源距的變化曲線，可見測量的慢度/極化矢量與真值完全吻合。圖2為測量的極化矢量和慢度矢量的相對誤差，可見極化矢量的相對誤差非常小，慢度矢量的水平分量誤差稍大一些，但基本都小于1%。

圖3為6條測線聯(lián)合反演的不同深度15個弱各向異性參數(shù)結(jié)果。反演的弱各向異性參數(shù)(散點)在不同深度上有所不同，但基本都在真值(直線) 的附近，其絕對誤差(圖4)最大為4.05×10-3，遠小于最小的各向異性參數(shù)，說明模型的反演結(jié)果精度很高。

圖5是利用反演的弱各向異性參數(shù)通過式(4)和式(5)正演計算的不同方位測線(30°/210°、 90°/270°和150°/330°)1250m深度的慢度和極化矢量與測量值的對比，可以看出二者高度吻合。圖6是利用式(4)正演的不同方位、不同傾角在1250m深度處接收點的相速度曲線，可見：120°/300°測線的相速度曲線對稱，說明這個方位是TTI對稱平面的垂直方向，那么對稱平面方位角則為30°/210°； 在對稱平面上，-40°傾角時相速度最小，說明這是對稱軸傾角，通過對稱平面上的相速度兩個極值計算的各向異性程度為4%，與TTI模型的設(shè)計值完全一致。表1給出了所有深度點計算的各向異性程度和TTI對稱軸方向，與模型設(shè)計值高度一致。

圖1 3條Walkaway測線測量的慢度、極化矢量與真值對比(a)30°/210°； (b)90°/270°； (c)150°/330°。上標“obs”表示測量值，“r”表示真值

圖2 3條Walkaway測線測量的慢度、極化矢量的相對誤差(a)30°/210°； (b)90°/270°； (c)150°/330°

圖3 弱各向異性參數(shù)反演結(jié)果與真值的對比(a)εx、δx、χx、ε15、ε24； (b)εy、δy、χy、ε16、ε35； (c)εz、δz、χz、ε26、ε34。上標“i”表示反演值，“r”表示真值

圖4 反演的弱各向異性參數(shù)的絕對誤差條形圖

圖5 3條Walkaway測線正演計算的慢度、極化矢量與測量值的對比(a)30°/210°； (b)90°/270°； (c)150°/330°。上標“obs”表示測量值，“cal”表示計算值，圖9、圖10同

表1 不同深度接收點處反演的各向異性程度和TTI對稱軸方向

圖6 不同方位測線速度隨傾角的變化曲線

3 實際資料應(yīng)用

以A井3條Walkaway VSP采集數(shù)據(jù)為例。測線布置如圖7，其中測線1布設(shè)450炮、測線2布設(shè)460炮、測線3布設(shè)470炮。接收點深度范圍為4560～5670m，共75 級，間距為15m。

基于上述3條測線，可反演9個弱各向異性參數(shù)，其中5個參數(shù)(εx、εz、δx、ε15和ε35)結(jié)果比較穩(wěn)定。

圖7 實際資料測線布置示意圖

圖8a和圖8b分別是反演的速度和5個弱各向異性參數(shù)隨深度的變化曲線。在地層分界面附近，尤其是第二和第三界面處，相速度有較大的抖動(圖8a)，反演估計的絕對誤差也較大(圖8c)，可能是由于在地層分界面上、下的接收點時差變化較大，造成垂直慢度的誤差較大。在鹽丘(第三條虛線以下)中各向異性參數(shù)有明顯的變大趨勢，說明鹽丘各向異性很強。從圖8c可發(fā)現(xiàn)，幾個ε參數(shù)的反演誤差都比δx小，這是因為ε只與qP波各向異性有關(guān)，δ不僅與qP波有關(guān)，還與qS波有關(guān)，而本文反演只用到了qP波信息，但反演誤差的最大值為2.296×10-4，仍遠小于弱各向異性參數(shù)值。

圖9和圖10分別是應(yīng)用反演的弱各向異性參數(shù)正演計算的3條測線的第一個接收點的慢度和極化矢量與測量值的對比，可以看出，測線1的慢度和極化矢量的計算值與測量值都吻合的較好，測線2和測線3基本能吻合，說明反演的弱各向異性參數(shù)正確地反映了測量的慢度和極化矢量所攜帶的地層各向異性信息。

圖11是計算的不同深度的各向異性程度及其各層的平均各向異性程度。如圖中所示，測線1的各向異性程度較弱，測線2、測線3較強，其中測線2最強。從深度上看，深度5500m以下鹽丘各向異性明顯增強，其中測線2最大，超過18%。

圖8 反演的速度(a)、弱各向異性參數(shù)(b)及其誤差(c)隨深度的變化曲線水平虛線表示地層界面

圖9 應(yīng)用反演的弱各向異性參數(shù)計算的慢度矢量與測量值的對比(a)測線1； (b)測線2； (c)測線3

圖10 應(yīng)用反演的弱各向異性參數(shù)計算的極化矢量與測量值的對比(a)測線1； (b)測線2； (c)測線3

圖11 3條測線不同深度的各向異性程度曲線(a)與地層平均各向異性程度(b)的對比

4 結(jié)束語

本文系統(tǒng)地介紹了基于擾動理論的qP波的慢度和極化矢量正、反演方法，給出了精度高、穩(wěn)定的慢度和極化矢量計算方法，提出了適用于慢度和極化矢量反演任意弱各向異性參數(shù)的加權(quán)迭代反演算法。

通過設(shè)計6條Walkaway VSP的TTI模型處理試驗，反演的15個弱各向異性參數(shù)精度很高； 應(yīng)用反演的弱各向異性參數(shù)計算的慢度和極化矢量與測量值完全匹配； 正演的相速度結(jié)果表明，模型各向異性程度和對稱軸與理論值完全一致。

實際3條測線Walkaway VSP資料反演的5個弱各向異性參數(shù)取值范圍和各向異性程度都符合地質(zhì)認識，各向異性隨深度變化趨勢與地層分界面對應(yīng)較好； 應(yīng)用弱各向異性參數(shù)計算的慢度和極化矢量與測量值吻合很好； 反演誤差遠小于各向異性參數(shù)的數(shù)量級。

模型和實際應(yīng)用表明，本文方法能夠用于高精度的地震各向異性參數(shù)提取，具備較高的應(yīng)用價值。