鄭浩，蔡杰雄，王靜波

(1.中國石化石油物探技術(shù)研究院,江蘇 南京 211103; 2.中國石化勘探分公司,四川 成都 610041)

0 引言

近年來隨著“兩寬一高”地震采集和層析反演方法取得突破進(jìn)展，地震成像分辨率得到有效拓寬[1-2], 高分辨率層析反演速度建模方法取得了長足的進(jìn)步。傳統(tǒng)的深度域射線層析反演速度建模方法的分辨率存在中波數(shù)段信息缺失問題[3]，其可以較好地恢復(fù)地下速度的背景低波數(shù)分量，但難以準(zhǔn)確反演中波數(shù)段的高精度速度模型[4-5]。學(xué)者們提出波動(dòng)方程層析反演方法[6]，將速度反演分辨率向中高波數(shù)端擴(kuò)展。該類方法基于波動(dòng)理論，避免了高頻近似假設(shè)引入的誤差，反演結(jié)果精度高，但該類方法存在計(jì)算效率問題，且反演精度非常依賴于數(shù)據(jù)品質(zhì)，當(dāng)數(shù)據(jù)較差時(shí)存在反演穩(wěn)定性問題[7-8]，在實(shí)際應(yīng)用中推廣較少。為了降低射線類層析矩陣的稀疏性，同時(shí)克服波動(dòng)方程層析的計(jì)算量與穩(wěn)定性問題，諸多學(xué)者提出介于常規(guī)射線理論和波動(dòng)方程理論之間的方法。早在1995年，Vasco[9]便提出一種介于射線和波動(dòng)方程之間的胖射線層析方法，該方法通過將射線加寬，提高射線路徑的投影范圍，從而提高層析精度，但該方法缺乏嚴(yán)格的理論推導(dǎo)，僅是對波動(dòng)理論感性認(rèn)識(shí)，應(yīng)用較為局限。菲涅爾體層析[10]及高斯束層析[11-12]的出現(xiàn)實(shí)現(xiàn)了對波動(dòng)方程線性化近似的“束”層析。菲涅爾體層析利用射線周圍的波動(dòng)性質(zhì)來擬合走時(shí)和波場振幅，相比于射線層析的高頻近似假設(shè)，該方法引入了菲涅爾體的有限頻效應(yīng)，更接近真實(shí)的波場傳播，這降低常規(guī)射線層析方法中矩陣的稀疏性，提高計(jì)算穩(wěn)定性[13]。高斯束層析[14]實(shí)際上是將臨近區(qū)域空間內(nèi)的多條射線合并為一條高斯束射線從而建立層析方程，對于觀測系統(tǒng)中的一個(gè)炮檢對，通過提高射線覆蓋范圍來提高反演精度及分辨率，降低層析方程的病態(tài)性，且該方法相對于波動(dòng)方程層析，計(jì)算速度快、計(jì)算結(jié)果穩(wěn)定。但早期的高斯束層析理論推導(dǎo)不夠嚴(yán)格，蔡杰雄[15]從高斯束偏移角度道集出發(fā)，在波動(dòng)方程一階Born近似和Rytov近似下，推導(dǎo)出了基于高斯束算子的成像域走時(shí)層析方程及其敏感度核函數(shù)，完善了高斯束層析反演的理論基礎(chǔ)，推動(dòng)了高斯束層析反演的實(shí)用化。

層析反演速度建模的另一個(gè)重要研究內(nèi)容是各種約束條件的利用，即在反演過程中加入包括層位、構(gòu)造產(chǎn)狀、井速度等多種已知信息約束以實(shí)現(xiàn)正則化，提高反演穩(wěn)定性及分辨率。目前，層析反演中的構(gòu)造約束方式主要有兩類，一類是在速度層析過程中引入地質(zhì)層位信息，通過層約束提高反演精度，但這類方法需要解釋人員進(jìn)行地質(zhì)構(gòu)造的精細(xì)化解釋，對解釋結(jié)果有較高的精度要求，若解釋結(jié)果出現(xiàn)錯(cuò)誤，會(huì)使得層析結(jié)果誤差較大。該類方法大都基于主觀認(rèn)知強(qiáng)行加入反演過程中，比如層位約束層析反演就是在本輪的反演中不修改層位而僅更新層位內(nèi)速度值，待完成本輪反演后重新偏移成像，在新剖面上重新進(jìn)行層位拾取，之后再加入到新一輪的層析速度更新中，通常我們稱這種約束為“強(qiáng)”約束[16-17]。強(qiáng)約束的反演具有雙刃性，它一方面可能快速促進(jìn)收斂，但另一方面可能引入解釋人員的主觀錯(cuò)誤信息，導(dǎo)致反演結(jié)果出現(xiàn)誤差；另一類方法通過引入含有構(gòu)造信息的正則化算子，用于導(dǎo)引速度更新方向，避免了解釋人員的主觀錯(cuò)誤信息，但構(gòu)造約束精度取決于正則化算子的精度。我們稱這種約束為“軟”約束。目前有多種預(yù)條件正則化算子，包括Sobel算子、Reborts算子、基于斜率的平滑濾波算子、基于擴(kuò)散方程解析解的構(gòu)造導(dǎo)向?yàn)V波算子等。這些方法通過提取地震偏移圖像的構(gòu)造特征，從而構(gòu)建預(yù)條件算子，實(shí)現(xiàn)構(gòu)造約束層析反演。

本文從高斯束層析核函數(shù)出發(fā)，通過引入含有構(gòu)造信息的預(yù)條件模型正則化算子，推導(dǎo)得到構(gòu)造導(dǎo)向?yàn)V波約束下的高斯束層析速度建模方法。在預(yù)條件模型正則化算子的設(shè)計(jì)上，本文引入結(jié)構(gòu)張量算法，用于提取偏移圖像中的平坦、邊緣及角點(diǎn)區(qū)域，通過擴(kuò)散方程解析解實(shí)現(xiàn)構(gòu)造導(dǎo)向平滑濾波，建立層析反演正則化算子，實(shí)現(xiàn)反演過程的自動(dòng)“軟”約束，提高高斯束速度層析反演的精度與穩(wěn)定性。

1 成像域高斯束層析方法

頻率域高斯束成像條件可以表示為：

IGBM(x,θ,φ,w)=S(x,pS,w;xS)R*(x,pR,w;xS)

(1)

式中:IGBM(x,θ,φ,w)表示共成像點(diǎn)道集,S(x,pS,w;xS)表示來自激發(fā)點(diǎn)的高斯束下行波場,R(x,pR,w;xS)表示檢波點(diǎn)接收的高斯束上行波場,x=(x,y,z)表示共成像點(diǎn)位置的坐標(biāo)，pS和pR分別代表激發(fā)點(diǎn)和接收點(diǎn)位置的慢度矢量，θ表示反射張角，φ表示成像點(diǎn)的方位角，w=2πf表示角頻率，*表示共軛。

從角度域高斯束偏移成像條件(1)出發(fā)，在波動(dòng)方程的一階Born近似下，可得到剩余時(shí)差與速度修正量的表達(dá)式，再引入Rytov近似可進(jìn)一步簡化得到剩余時(shí)差與速度修正量的線性關(guān)系式[18]：

(2)

其中，Δt(x,θ,φ,ω)為高斯束偏移的剩余時(shí)差，KT表示高斯束層析核函數(shù)，其可以表示為兩部分：

基于式(2)表示的層析核函數(shù)，就可以建立類似于射線層析的方程:

AΔm=Δt。

(4)

圖1 深度域高斯束層析反演核函數(shù)Fig.1 Gaussian beam tomography kernel function in imaging domain

A是利用式(2)層析核函數(shù)計(jì)算得到的敏感度核函數(shù)矩陣，Δm是待反演的速度更新量，Δt是共成像點(diǎn)道集(common image gathers,CIG)的走時(shí)殘差，在反射點(diǎn)局部，可以假定炮點(diǎn)和檢波點(diǎn)的慢度pS和pR是常數(shù)，那么該點(diǎn)的慢度可以表示為：

Δpz=Average(pS+pR) 。

(5)

通過在CIG道集上拾取剩余深度差Δz(RMO)，則可利用式(5)得到走時(shí)殘差Δt：

Δt=ΔpzΔz。

(6)

通常，式(4)中的矩陣A是極度稀疏的，直接求解穩(wěn)定性較差，且多解性較強(qiáng)，往往需要引入正則化項(xiàng)，加入一些先驗(yàn)約束，使解在約束條件下收斂。這里采用預(yù)條件模型正則化的方式加入先驗(yàn)信息，待反演的模型擾動(dòng)量可表示為：

Δm=Du。

(7)

其中算子D表示預(yù)條件算子，即引入的先驗(yàn)信息，這樣可以將Δm的求解問題轉(zhuǎn)化為算子u的求解問題，將式(7)代入到式(4)中，可得：

ADu=Δt。

(8)

式(8)即為構(gòu)建的預(yù)條件層析方程，其阻尼最小二乘方程可表示為：

DTATADu+εu=DTATΔt,

(9)

求解上式得到u，代入到式(7)中即可得到速度更新量Δm。

通過上述推導(dǎo)，實(shí)現(xiàn)了波動(dòng)方程線性化的深度域走時(shí)層析反演方法，再與疊前深度偏移相結(jié)合進(jìn)行迭代：速度誤差帶來成像不聚焦誤差，表現(xiàn)為成像道集走時(shí)差；再利用成像誤差反過來修正速度，實(shí)現(xiàn)兩者的迭代收斂，獲得更精確的速度場和成像結(jié)果。

2 基于結(jié)構(gòu)張量構(gòu)造導(dǎo)向?yàn)V波正則化技術(shù)

從式(9)中可以看出，預(yù)條件模型正則化的高斯束層析反演關(guān)鍵是算子D的構(gòu)建，傳統(tǒng)方法常選擇空間平滑算子，主要是通過減小層析系數(shù)矩陣的條件數(shù)，緩解反問題病態(tài)性，在矩陣求解的穩(wěn)定性方面有一定作用，但該方法容易模糊構(gòu)造邊界，其反演結(jié)果對于邊界的刻畫較差。本文采用各向異性擴(kuò)散方程解析解[20]構(gòu)建預(yù)條件算子D，其隱式表達(dá)式為如下的偏微分方程：

g(x)-α·T(x)g(x)=f(x) ,

(10)

其中:f(x)是待濾波的地震數(shù)據(jù)，g(x)是濾波后的地震數(shù)據(jù)，α是表示平滑力度強(qiáng)弱的常數(shù)；T(x)是從地震剖面中計(jì)算得到的結(jié)構(gòu)張量。結(jié)構(gòu)張量實(shí)際上相當(dāng)于沿構(gòu)造梯度的一個(gè)平滑解，它可以用于構(gòu)造的傾角掃描，且平滑過程沿構(gòu)造進(jìn)行，達(dá)到保邊界的效果。對于一個(gè)三維數(shù)據(jù)體而言，其某點(diǎn)的結(jié)構(gòu)張量T(x)可以通過一個(gè)3×3的對稱正定矩陣表示：

(11)

其中:dx、dy、dz分別為數(shù)據(jù)沿x，y，z三個(gè)方向的梯度。通過矩陣的特征值分解，上式可以表示為:

T(x)=λuuuT+λvvvT+λwwwT,

(12)

其中:λu、λv和λw為分布在0～1的正實(shí)數(shù)，表示矩陣T(x)的特征值，u、v和w是對應(yīng)的特征向量。通常定義λu≥λv≥λw≥0，這種情況下，特征向量u表示梯度最大的方向，即表示垂直于構(gòu)造傾角的方向。特征向量v和w分別表示inline和crossline平行于構(gòu)造的方向，λ越大意味沿對應(yīng)特征方向的平滑尺度越大。

這里以二維數(shù)據(jù)為例，對特征向量進(jìn)行說明。二維情況下，式(12)只有前兩項(xiàng)，圖2a、b中黃色的橢圓分別是地震切片與剖面上結(jié)構(gòu)張量計(jì)算結(jié)果，實(shí)際上在無構(gòu)造(各向同性)區(qū)域，該點(diǎn)的結(jié)構(gòu)張量形態(tài)是圓形，如圖3a所示，此時(shí)u和v的長度相等；當(dāng)存在構(gòu)造形態(tài)(各向異性)時(shí)，就會(huì)出現(xiàn)橢圓形態(tài)的結(jié)構(gòu)張量，如圖3b所示，其中v表示平行于構(gòu)造方向的特征向量，u表示垂直于構(gòu)造方向的特征向量，這樣我們就可以計(jì)算出所有點(diǎn)的結(jié)構(gòu)張量。

圖2 地震切片(a)及剖面(b)計(jì)算得到的結(jié)構(gòu)張量Fig.2 Structural tensor calculated on seismic slice (a) and section (b)

得到每點(diǎn)的結(jié)構(gòu)張量后，通過合理調(diào)整λu、λv、λw，從而實(shí)現(xiàn)沿層位和斷層的構(gòu)造導(dǎo)向?yàn)V波處理，平滑效果如圖4所示。

圖3 各向同性(a)和各向異性(b)的結(jié)構(gòu)張量算子Fig.3 Isotropic(a) and anisotropic(b) structural tensor operators

圖4 構(gòu)造導(dǎo)向?yàn)V波前(a)后(b)對比，可同時(shí)加強(qiáng)層位與斷層特征Fig.4 The comparison of (a) before and (b) after structural-guided filtering,which can enhance the structure and fault characteristics

引入有限差分近似，預(yù)條件算子D可寫作：

D=(I+TT)-1,

(13)

其中T即為式(12)所求的結(jié)構(gòu)張量，單位矩陣I的作用是在λu、λv、λw均為零時(shí)保持圖像不改變，將式(13)所得結(jié)果代入式(9)中即可實(shí)現(xiàn)基于構(gòu)造導(dǎo)向?yàn)V波的高斯束層析反演。此外，相比于其他的構(gòu)造約束算法，該方法無需指定平滑角度，完全自動(dòng)提取，計(jì)算靈活；且方法的抗噪性較好，即便是在信噪比極低的情況下，構(gòu)造導(dǎo)向?yàn)V波算法也可以得到相對精確的結(jié)果。

3 模型及實(shí)際資料測試

3.1 二維層狀模型數(shù)據(jù)測試

設(shè)計(jì)二維層狀地質(zhì)模型驗(yàn)證本文方法的精度，模型中填充速度屬性值。所設(shè)計(jì)模型為速度沿縱向遞增的層狀模型， 其中模型中部有一高速隆起構(gòu)造(圖5a)，縱向采樣點(diǎn)為601，采樣間隔為10 m，橫向共有1201道，道間距為10 m。利用聲波波動(dòng)方程正演得到單炮記錄(圖6)，觀測系統(tǒng)采用中心放炮，兩邊接收，每炮801道。初始模型采用常梯度模型，見圖5，其中第一層速度給成準(zhǔn)確值。圖5c、5d展示了常規(guī)高斯束層析與構(gòu)造約束預(yù)條件高斯束層析的反演結(jié)果，顯然，常規(guī)的高斯束層析對于淺層的水平層狀構(gòu)造反演結(jié)果尚可，但對于中部隆起構(gòu)造反演結(jié)果較差，分辨率明顯不足，且中深部反演結(jié)果誤差較大；相比較而言，基于構(gòu)造導(dǎo)向?yàn)V波的高斯束層析算法所得結(jié)果精度較高，與真實(shí)模型基本吻合，更有利于精確成像。

為了更加直觀地看出造構(gòu)造導(dǎo)向?yàn)V波高斯束層析的優(yōu)勢，在層析結(jié)果豎直抽取某一道進(jìn)行對比，如圖7a所示。圖7b表示了理論速度值、平滑速度值和反演的速度值的對比曲線。，顯然，基于構(gòu)造導(dǎo)向?yàn)V波的高斯束層析反演方法不僅可以反演得到中波數(shù)分量的結(jié)果，提供更加豐富的速度信息，且能夠達(dá)到保邊界的效果，構(gòu)造邊界刻畫更為清晰，這說明基于構(gòu)造導(dǎo)向?yàn)V波的高斯束層析對于構(gòu)造的刻畫明顯優(yōu)于常規(guī)高斯束層析方法，且反演精度更高，減少了由于人為干預(yù)帶來的誤差，結(jié)果更加合理。

a—理論模型;b—初始模型;c—常規(guī)高斯束層析反演結(jié)果;d—構(gòu)造導(dǎo)向?yàn)V波高斯束層析反演結(jié)果a—velocity model;b—initial model;c—updated velocity model by conventional Gaussian beam tomography;d—updated velocity model by Gaussian beam tomography based on structural-guided filtering圖5 理論地質(zhì)模型(速度值)及反演結(jié)果Fig.5 Velocity model and the inverted results

圖6 模型正演的單炮記錄Fig.6 The modelling common shot gathers

a—抽取模型中心道;b—真實(shí)值、常規(guī)平滑正則化約束層析反演結(jié)果及構(gòu)造約束正則化層析反演結(jié)果單道對比a—the location of the extracted trace in the model;b—comparison of theoretical result,conventional smoothing regularization tomography inversion result and structural-guided regularization tomography result圖7 反演結(jié)果單道對比分析Fig.7 Single trace comparison

深度域高斯束層析反演技術(shù)可以緩解層析矩陣的病態(tài)問題，反演得到速度模型的中波束分量，提高地震勘探的成像精度?；诖?，通過引入基于結(jié)構(gòu)張量的構(gòu)造導(dǎo)向?yàn)V波算法用于約束速度更新量的計(jì)算，能反演出高精度的地下速度模型，增強(qiáng)斷層突變位置的反演效果，彌補(bǔ)地震成像的中波數(shù)段缺失，提高地震成像精度。

3.2 實(shí)際地震資料應(yīng)用

實(shí)際資料選自中國西南部復(fù)雜山地頁巖氣某工區(qū)。該地區(qū)地表起伏大，懸崖峭壁密布，山高谷深，數(shù)據(jù)信噪比非常低，且地下結(jié)構(gòu)復(fù)雜，速度橫向變化劇烈，速度建模難度大，成像困難，具有典型的復(fù)雜地表和復(fù)雜構(gòu)造的“雙復(fù)雜”特性，其對地震資料成像精度要求遠(yuǎn)高于常規(guī)勘探。

首先對該區(qū)域的資料進(jìn)行預(yù)處理、疊前時(shí)間偏移、深度域初始建模及偏移成像等基礎(chǔ)工作，得到高品質(zhì)的CIG道集，在此基礎(chǔ)上，通過基于構(gòu)造導(dǎo)向?yàn)V波的高斯束層析反演得到深度域速度模型，結(jié)合克希霍夫疊前深度偏移得到最終的成像剖面。圖8是某實(shí)際數(shù)據(jù)深度域常規(guī)高斯束層析建模和構(gòu)造導(dǎo)向?yàn)V波高斯束層析速度建模結(jié)果對比，從兩張對比圖上可以看出，構(gòu)造導(dǎo)向?yàn)V波高斯束層析速度建模結(jié)果表現(xiàn)了速度模型更多細(xì)節(jié)，清晰地刻畫了速度縱橫向空間變化，提高了速度建模的精度。

圖8 某實(shí)際數(shù)據(jù)深度域常規(guī)高斯束層析建模(a)及構(gòu)造導(dǎo)向?yàn)V波高斯束層析建模(b)結(jié)果對比Fig.8 Comparison of conventional Gaussian beam tomography (a) and structure-guided filter Gaussian beam tomography (b) in depth domain

速度模型的更新和剩余速度分析采用速度更新量和剩余速度譜聯(lián)合質(zhì)量控制方法，圖9是質(zhì)控圖,圖10是更新前后的CIG道集。從更新前的道集上看，存在同相軸未校平現(xiàn)象，更新后這些同相軸變平，圖上前頭所指。由于存在各向異性，還有一些同相軸未校平。另外，該區(qū)域存在多次波的影響(如圖中紅框區(qū)域所示)，因此仍有部分道集下拉現(xiàn)象。

圖9 速度更新量(a)與剩余譜(b)聯(lián)合質(zhì)控Fig.9 Joint quality control of velocity update (a) and residual spectrum (b)

圖10 速度更新前(a)后(b)CIG道集Fig.10 The CIGs before (a) and after (b) velocity update

圖11是對應(yīng)于圖8速度模型范圍的偏移結(jié)果，顯然，采用相同的成像方法，構(gòu)造導(dǎo)向?yàn)V波下的高斯束層析反演對應(yīng)的偏移結(jié)果中斷點(diǎn)更加連續(xù)，復(fù)雜構(gòu)造區(qū)信噪比更高，整體成像質(zhì)量明顯優(yōu)于常規(guī)方法所得結(jié)果，這對后續(xù)的構(gòu)造落實(shí)有較好的指導(dǎo)意義。通過該實(shí)際資料處理過程可以認(rèn)識(shí)到，基于構(gòu)造導(dǎo)向?yàn)V波的高斯束層析反演速度建模能夠反演得到速度的中高波數(shù)分量，所得成像結(jié)果細(xì)節(jié)更加豐富，信噪比高，成像質(zhì)量好。

a—常規(guī)高斯束層析反演對應(yīng)的偏移結(jié)果;b—構(gòu)造導(dǎo)向?yàn)V波高斯束層析反演對應(yīng)的偏移結(jié)果a—the imaging results corresponding to conventional Gausssian beam tomography;b—the imaging results corresponding to structural-guided filter Gaussian beam tomography圖11 某實(shí)際數(shù)據(jù)PSDM成像結(jié)果Fig.11 PSDM imaging results

4 結(jié)論

基于波動(dòng)方程線性化近似的深度域高斯束層析速度建模方法具有較高的反演精度與實(shí)用性。在此基礎(chǔ)上，本文引入結(jié)構(gòu)張量算法，推導(dǎo)了構(gòu)造導(dǎo)向?yàn)V波正則化技術(shù)。該技術(shù)通過計(jì)算地震數(shù)據(jù)的結(jié)構(gòu)張量，并利用各向異性光滑算子將其應(yīng)用于高斯束層析反演中，作為正則化約束項(xiàng)，實(shí)現(xiàn)了構(gòu)造導(dǎo)向?yàn)V波約束下的高斯束層析速度建模方法。相比于常規(guī)的高斯束層析速度建模技術(shù)，該方法完全基于數(shù)據(jù)驅(qū)動(dòng)，實(shí)現(xiàn)了對層析反演的“軟約束”，既緩解了層析反演求解的多解性，同時(shí)提高反演分辨率與穩(wěn)定性，得到更加符合地質(zhì)認(rèn)識(shí)的高精度速度模型，有效改善成像效果，提高了高斯束層析反演對復(fù)雜模型的反演分辨率，為高斯束實(shí)用化推廣鋪平了道路。

進(jìn)一步可考慮將基于結(jié)構(gòu)張量的構(gòu)造導(dǎo)向?yàn)V波算法應(yīng)用于各向異性層析建模，通過更高精度且自動(dòng)化地計(jì)算構(gòu)造的方位角及傾角推進(jìn)各向異性高斯束層析方法的實(shí)用化。