鄭浩，蔡杰雄，王靜波

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

0 引言

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

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

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

1 成像域高斯束層析方法

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

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

(1)

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

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

(2)

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

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

AΔm=Δt。

(4)

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

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

Δpz=Average(pS+pR) 。

(5)

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

Δt=ΔpzΔz。

(6)

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

Δm=Du。

(7)

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

ADu=Δt。

(8)

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

DTATADu+εu=DTATΔt,

(9)

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

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

2 基于結構張量構造導向濾波正則化技術

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

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

(10)

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

(11)

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

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

(12)

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

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

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

得到每點的結構張量后，通過合理調整λu、λv、λw，從而實現(xiàn)沿層位和斷層的構造導向濾波處理，平滑效果如圖4所示。

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

圖4 構造導向濾波前(a)后(b)對比，可同時加強層位與斷層特征Fig.4 The comparison of (a) before and (b) after structural-guided filtering,which can enhance the structure and fault characteristics

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

D=(I+TT)-1,

(13)

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

3 模型及實際資料測試

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

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

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

a—理論模型;b—初始模型;c—常規(guī)高斯束層析反演結果;d—構造導向濾波高斯束層析反演結果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 理論地質模型(速度值)及反演結果Fig.5 Velocity model and the inverted results

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

a—抽取模型中心道;b—真實值、常規(guī)平滑正則化約束層析反演結果及構造約束正則化層析反演結果單道對比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 反演結果單道對比分析Fig.7 Single trace comparison

深度域高斯束層析反演技術可以緩解層析矩陣的病態(tài)問題，反演得到速度模型的中波束分量，提高地震勘探的成像精度。基于此，通過引入基于結構張量的構造導向濾波算法用于約束速度更新量的計算，能反演出高精度的地下速度模型，增強斷層突變位置的反演效果，彌補地震成像的中波數(shù)段缺失，提高地震成像精度。

3.2 實際地震資料應用

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

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

圖8 某實際數(shù)據(jù)深度域常規(guī)高斯束層析建模(a)及構造導向濾波高斯束層析建模(b)結果對比Fig.8 Comparison of conventional Gaussian beam tomography (a) and structure-guided filter Gaussian beam tomography (b) in depth domain

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

圖9 速度更新量(a)與剩余譜(b)聯(lián)合質控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是對應于圖8速度模型范圍的偏移結果，顯然，采用相同的成像方法，構造導向濾波下的高斯束層析反演對應的偏移結果中斷點更加連續(xù)，復雜構造區(qū)信噪比更高，整體成像質量明顯優(yōu)于常規(guī)方法所得結果，這對后續(xù)的構造落實有較好的指導意義。通過該實際資料處理過程可以認識到，基于構造導向濾波的高斯束層析反演速度建模能夠反演得到速度的中高波數(shù)分量，所得成像結果細節(jié)更加豐富，信噪比高，成像質量好。

a—常規(guī)高斯束層析反演對應的偏移結果;b—構造導向濾波高斯束層析反演對應的偏移結果a—the imaging results corresponding to conventional Gausssian beam tomography;b—the imaging results corresponding to structural-guided filter Gaussian beam tomography圖11 某實際數(shù)據(jù)PSDM成像結果Fig.11 PSDM imaging results

4 結論

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

進一步可考慮將基于結構張量的構造導向濾波算法應用于各向異性層析建模，通過更高精度且自動化地計算構造的方位角及傾角推進各向異性高斯束層析方法的實用化。