王宇翔， 楊鍇， 楊小椿， 薛冬， 陳寶書

1 同濟(jì)大學(xué)海洋地質(zhì)國家重點實驗室， 上海　200092 2 Broadgeo, Inc, Houston, TX 770994 3 中海石油研究中心, 北京　100027

?

基于梯度平方結(jié)構(gòu)張量算法的高密度二維立體層析反演

王宇翔1,2， 楊鍇1*， 楊小椿3， 薛冬3， 陳寶書3

1 同濟(jì)大學(xué)海洋地質(zhì)國家重點實驗室， 上海200092 2 Broadgeo, Inc, Houston, TX 770994 3 中海石油研究中心, 北京100027

摘要射線參數(shù)水平分量是立體層析數(shù)據(jù)空間中最為重要的參數(shù)信息，梯度平方結(jié)構(gòu)張量算法是一種針對圖像的邊緣檢測快速算法.本文將疊前地震數(shù)據(jù)集視為圖像，將梯度平方結(jié)構(gòu)張量算法用于射線參數(shù)水平分量的提取，提高了立體層析數(shù)據(jù)空間的準(zhǔn)備效率.高效率的數(shù)據(jù)空間提取使得高密度立體層析反演成為可能，數(shù)據(jù)空間加密后的立體層析反演精度以及疊前深度成像質(zhì)量相比常規(guī)流程得到明顯提高.基于二維理論數(shù)據(jù)與南海某二維深水?dāng)?shù)據(jù)的嚴(yán)格測試證實了該方法的有效性與穩(wěn)健性.

關(guān)鍵詞立體層析反演； 數(shù)據(jù)空間提??； 梯度平方結(jié)構(gòu)張量算法； 偏移速度建模

1引言

立體層析反演是層析反演類速度建模方法中極具特色的一種方法.其將一個局部相干同相軸在炮道集與檢波點道集內(nèi)的射線參數(shù)水平分量(Psx,Prx)納入到反演的數(shù)據(jù)空間之中，不僅使得數(shù)據(jù)空間的準(zhǔn)備相比傳統(tǒng)反射層析更為簡便，也使得立體層析反演成為運動學(xué)層析反演方法中唯一一種可以同時反演速度結(jié)構(gòu)與反射點位置的方法(Billette and Lambare, 1998；Lambare, 2008；倪瑤等，2013；任麗娟等，2014；李振偉等，2014).近年來亦有嘗試將立體層析反演與有限頻理論相結(jié)合，以期獲得更好的反演結(jié)果(Ni et al.，2012).

在立體層析反演中，數(shù)據(jù)空間可以寫為d=[Sx,Rx,T,Psx,Prx].其中Sx,Rx為炮檢點坐標(biāo)；T為總走時；Psx,Prx為炮檢點處的射線參數(shù)水平分量.注意對有效數(shù)據(jù)點的確定完全是根據(jù)局部波包的Psx,Prx是否足夠精確、可靠作為標(biāo)準(zhǔn).換言之，數(shù)據(jù)點的所有信息中，Psx,Prx首先得到確定，之后再確定該局部波包對應(yīng)的炮檢點坐標(biāo)Sx,Rx和走時T.因此，射線參數(shù)水平分量的精度問題是立體層析反演數(shù)據(jù)準(zhǔn)備過程中首要關(guān)心的問題.

由于射線參數(shù)的確定實質(zhì)就是同相軸的局部傾角估計問題，前人大多使用局部傾斜疊加或平面波摧毀濾波器來獲得這類信息.本文引入圖像處理界中發(fā)展成熟的邊緣檢測算法——結(jié)構(gòu)張量算法來估計同相軸的局部傾角.石油工業(yè)界中已有人將該類算法用于地震解釋中的自動層位追蹤或者提高信噪比的構(gòu)造預(yù)測濾波(Wu and Hale, 2013; Hale, 2009).不同于上述兩個文獻(xiàn)，本文使用的是基于梯度平方的結(jié)構(gòu)張量算法.實際測試表明該算法具有很好的抗噪性，將其用于實際地震數(shù)據(jù)的局部傾角估計是非常合適的.引入該算法的另一個重要理由是它的計算效率遠(yuǎn)高于剛才提到的兩種常規(guī)算法.該算法的引入將使得立體層析數(shù)據(jù)點的分布密度大幅提高而無需擔(dān)心計算成本的大幅上升.這是實施高密度立體層析反演的一個基本前提保證.

本文首先簡要回顧立體層析反演的方法原理并闡明射線參數(shù)水平分量對立體層析數(shù)據(jù)空間中的重要地位.然后詳細(xì)介紹梯度平方結(jié)構(gòu)張量算法，并將其與局部傾斜疊加和平面波摧毀濾波器的計算效率和性能進(jìn)行對比，討論如何將該算法應(yīng)用于疊前地震數(shù)據(jù)、從而實現(xiàn)智能化的自動傾角拾取.在此基礎(chǔ)上，作者提出了一種基于物理規(guī)律的嚴(yán)格質(zhì)量監(jiān)控(QC)方法來過濾非物理的、冗余的數(shù)據(jù)點.需要指出，正是由于梯度平方結(jié)構(gòu)張量算法的高效以及合乎物理規(guī)律的嚴(yán)格QC的引入，使得適用于實際數(shù)據(jù)的高密度立體層析反演成為可能.基于二維理論數(shù)據(jù)與南海某二維深水地震數(shù)據(jù)的測試表明，上述策略的使用不但使得數(shù)據(jù)點的數(shù)量大幅度增加，同時也保證了數(shù)據(jù)點的有效性和可靠性.這對于降低層析矩陣的病態(tài)性、提高反演的精度是至關(guān)重要的.

2立體層析反演方法原理介紹

首先介紹立體層析反演的模型空間、數(shù)據(jù)空間以及層析矩陣的建立過程.為簡明起見，以二維為例.圖1a顯示了一根從炮點S出發(fā)、到地下反射點X反射、回到地表R的射線.從透射的角度，不妨將其理解為從X出發(fā)、分別以入射角θs與反射角θr發(fā)射至炮點和檢波點的兩根透射射線.當(dāng)速度正確時，兩根透射射線分別以正確的出射方向、正確的走時達(dá)到正確的炮點和檢波點位置.這樣構(gòu)成的立體層析數(shù)據(jù)空間各分量為：d=[Sx,Rx,T,Psx,Prx].其中Sx,Rx為炮檢點坐標(biāo)；T為走時；Psx,Prx為炮檢點處的射線參數(shù)水平分量.地下的模型空間m=(X0,Z0,θs,θr,V)，其中X0,Z0為地下反射點X的坐標(biāo)；θs,θr分別代表從炮點一側(cè)與檢波點一側(cè)出射的地下張角,V代表地下介質(zhì)速度信息.

圖1　二維立體層析數(shù)據(jù)空間各分量與模型空間各分量(a) 立體層析的模型與數(shù)據(jù)分量; (b) 同相軸斜率與射線理論的慢度矢量.

注意圖1b中所示的Psx可以在共檢波點道集內(nèi)搜索同相軸的局部傾角得到，Prx可以在共炮點道集內(nèi)搜索同相軸的局部傾角得到.原理如下：在圖1b所顯示的炮記錄中，如果考察其中某一局部相干同相軸，由幾何關(guān)系以及射線理論中慢度矢量的定義得

(1)

式中，Pslope為局部相干同相軸的斜率，Prx為檢波點處慢度矢量的水平分量.式(1)表明共炮數(shù)據(jù)上同相軸的斜率對應(yīng)于檢波點處慢度矢量的水平分量Psx，由炮檢互易原理可知，共檢波點數(shù)據(jù)上同相軸的斜率必對應(yīng)于炮點處慢度矢量的水平分量.

Billette和Lambare(1998)指出，正是由于Prx,Psx(共炮道集和共檢波點道集數(shù)據(jù)內(nèi)同相軸的斜率)被引入層析數(shù)據(jù)空間，使得在層析反演中反射/散射波在炮檢點處的局部傳播方向得到強(qiáng)有力的約束.具體的優(yōu)勢體現(xiàn)在立體層析反演中無需拾取連續(xù)的反射同相軸，只要選擇連續(xù)性好、信噪比高的局部波包即可.同時，將反射問題化為透射問題，也回避了傳統(tǒng)反射層析中速度與反射層深度耦合的問題.由于在數(shù)據(jù)空間中引入了Prx,Psx，模型空間必然也隨之?dāng)U大.因此在反演速度的同時必須還要反演反射點的深度位置以及反射層的傾角.這使得立體層析反演成為運動學(xué)層析反演方法中唯一一種可以同時反演速度、反射點位置以及局部傾角的層析反演方法.而常規(guī)反射層析反演需要在速度更新之后再單獨更新反射層深度，這通常是一個需要分兩步實施的過程.

公式(2)展示了二維立體層析反演的層析矩陣(或FRECHET偏導(dǎo)數(shù)矩陣)，其中σ為針對層析矩陣每一行物理量的數(shù)量級不同設(shè)置的均衡系數(shù).可以看到常規(guī)層析反演由于其模型空間一般僅由速度V構(gòu)成，其數(shù)據(jù)空間一般僅考慮走時T，可認(rèn)為其層析矩陣僅相當(dāng)于公式(2)中的右上角那一部分.可以看出立體層析矩陣無論從規(guī)模還是其稀疏性方面都超過了常規(guī)的走時層析矩陣.如果能夠?qū)崿F(xiàn)高密度、高質(zhì)量的數(shù)據(jù)空間提取，對于改善層析矩陣的條件數(shù)，降低求解時對規(guī)則化的要求，進(jìn)而提高求解精度顯然都是有意義的.

F=

(2)

Lambare(2008)指出，選擇立體層析數(shù)據(jù)空間的標(biāo)準(zhǔn)就是優(yōu)先考慮連續(xù)性好、信噪比高的一次反射局部波包.從數(shù)據(jù)空間的特點不難看出，一旦選擇好了Psx和Prx參數(shù)，其他參數(shù)如炮檢點橫坐標(biāo)和走時都可以輕易隨之得到.如前所述，Psx和Prx參數(shù)估計等價于同相軸局部傾角估計.換言之，如何準(zhǔn)確、快速地估計共炮點道集與共檢波點道集內(nèi)同相軸的局部傾角是立體層析反演中最為重要的數(shù)據(jù)準(zhǔn)備工作.

3梯度平方結(jié)構(gòu)張量算法原理及其特點

工業(yè)界兩種最常用的局部傾角提取手段當(dāng)屬局部傾斜疊加(Ottolini, 1983)與平面波摧毀濾波器(Fomel, 2002).前者在時空域或者頻率空間域針對目標(biāo)同相軸實施局部疊加，通過疊加能量最大找到最合適的局部傾角；后者主要通過對局部平面波方程進(jìn)行有限差分分解，在頻率域內(nèi)構(gòu)建預(yù)測算子，然后再轉(zhuǎn)換到時間域時對頻率域相移算子時進(jìn)行近似得到.這兩種算法都被工業(yè)界大量使用.但是如果從計算成本的角度來考量，圖像處理領(lǐng)域的梯度平方結(jié)構(gòu)張量算法則具有巨大的優(yōu)勢.

圖2　某海洋二維測線原始炮記錄

圖2顯示了南海某二維深水測線中的一個單炮記錄，其信噪比不是很理想.從圖像處理的角度，這個信噪比不高的地震記錄可以視作一張反差不太強(qiáng)烈的二維圖像.而針對一張信噪比不高的地震記錄如何去估算其每一個樣點處的局部傾角和圖像處理中對于反差不夠強(qiáng)烈的圖像如何實施邊緣檢測，本質(zhì)上是同一個問題.梯度平方結(jié)構(gòu)張量算法正是針對這個問題而提出的(Van Vliet and Verbeek, 1995).

對于任何一張圖像，為檢測其邊緣，必須首先計算出每一個樣點處X方向與Y方向的梯度gx,gy.而高效計算每一樣點的梯度在圖像處理界早已有成熟算法，這就是結(jié)構(gòu)張量算法中使用的迭代高斯濾波.需要指出，這正是結(jié)構(gòu)張量算法具有高計算效率的原因所在.關(guān)于迭代高斯濾波的技術(shù)細(xì)節(jié)可以參考相關(guān)文獻(xiàn)(Van Vliet et al.，1998).

針對低信噪比圖像，為了提高方向預(yù)測的穩(wěn)定性，避免梯度矢量在圖像邊緣兩側(cè)互相抵消， Van Vliet和Verbeek (1995)引入了所謂的梯度平方張量矩陣，對圖像中某一樣點而言，其梯度平方張量矩陣的定義如下：

(3)

其中g(shù)x，gy分別表示該樣點處的梯度水平分量與梯度垂直分量.這個梯度平方張量表達(dá)了圖像中某一樣點處，對應(yīng)于單一走向的梯度方向.但是對于低信噪比圖像，為了提高其預(yù)測走向信息的穩(wěn)定性，需要在該樣點附近的鄰域內(nèi)，如對一個5×5的區(qū)域統(tǒng)計25個這樣的梯度平方張量，并將它們加權(quán)疊加獲得一個光滑的梯度平方張量矩陣G′才能獲得關(guān)于該樣點比較可靠的走向信息.這里將G′寫為

(4)

其中〈·〉代表光滑后的梯度平方張量.注意梯度平方張量矩陣的引入使得同一走向但是方向相反的梯度矢量不至于相互抵消，反而可以相互增強(qiáng).

獲得光滑的梯度平方張量矩陣G′之后，針對半正定矩陣G′，求解其特征值與特征向量，具體可由求解特征方程|G′-λ I|=0得到：

(5)

其中,λ1為最大特征值，對應(yīng)于張量能量在第一個特征張量方向V1的能量;λ2為最小特征值，對應(yīng)于張量能量在第二個特征張量方向V2的能量.(λ1-λ2)/λ1表示線性度，反映局部方向的一致性.注意如果基于地震數(shù)據(jù)實施該算法，這個參數(shù)其實表達(dá)了同相性的強(qiáng)弱程度.

特征向量則描述了圖像局部線性結(jié)構(gòu)的方向性.如圖2上的白色箭頭就是作者使用該算法計算中的兩個特征方向，特征向量V1正交于圖像的主結(jié)構(gòu)方向，特征向量V2平行于圖像的主結(jié)構(gòu)方向.圖3則是Van Vliet和Verbeek(1995)對一張信噪比欠佳的巖芯照片實施該算法后的結(jié)果.可以看到各類參數(shù)表達(dá)的圖像信息具有良好的一致性，表現(xiàn)出該算法具有不錯的穩(wěn)定性和抗噪性.

圖4展示了對于圖2所示的實際單炮記錄，應(yīng)用該算法提取的四種屬性參數(shù)剖面.可以看出由于梯度平方結(jié)構(gòu)張量算法的穩(wěn)定和高效，它能夠勝任針對實際地震數(shù)據(jù)的局部傾角估計工作.

對于立體層析而言，似乎只要知道特征方向V2就足夠了.實際上，我們同樣需要知道特征方向V1，以及相應(yīng)的最大與最小特征值λ1，λ2.因為線性度是非常重要的一個參數(shù)，其計算依賴于最大與最小特征值λ1，λ2.在后面的章節(jié)將會提及，線性度的大小正是評價局部波包同相性的一個重要指標(biāo)，也是我們評價數(shù)據(jù)點質(zhì)量的第一個標(biāo)準(zhǔn).

圖5a顯示了針對一個縱橫向樣點數(shù)為200×200的測試數(shù)據(jù)，其中含有少量隨機(jī)噪聲.應(yīng)用局部傾斜疊加、平面波摧毀濾波器與梯度平方結(jié)構(gòu)張量方法估算出的圖像局部傾角信息分別如圖5b,5c,5d所示.三種算法得到的局部傾角信息大致相似，但是計算效率則有明顯不同.表1中是三種算法的計算時間對比.

表1　三種算法的計算時間對比(單位：s)

從表1可見梯度平方結(jié)構(gòu)張量算法的計算效率高于平面波摧毀近一個數(shù)量級，高于局部傾斜疊加達(dá)兩個數(shù)量級.注意梯度平方結(jié)構(gòu)張量算法獲得的局部傾角信息是一個鄰域內(nèi)的平均效應(yīng).因此其變化比較光滑連續(xù)，不像局部傾斜疊加獲得的局部傾角信息似有突變.但是這并不表示其分辨率低于局部傾斜疊加，相反其傾角估計的穩(wěn)定性恰恰來自這種鄰域平均效應(yīng).需要指出，即便考慮了這個鄰域平均效應(yīng)，它在計算成本方面相對于另外兩種常規(guī)算法依然存在巨大優(yōu)勢.這使得它成為立體層析反演實際應(yīng)用中，提取射線參數(shù)水平分量的首選算法.

圖3　梯度平方張量算法應(yīng)用效果(Van Vliet and Verbeek，1995)(a) 輸入圖像; (b) 特征值λ1剖面; (c) 特征值λ2剖面; (d) 局部走向剖面; (e) 線性度(同相性)剖面.

圖4　梯度平方張量算法在某單炮記錄上的應(yīng)用效果(a) 輸入圖像; (b) 特征值λ1剖面; (c) 特征值λ2剖面; (d) 局部走向剖面; (e) 線性度(同相性)剖面.

圖5　三種傾角估計算法結(jié)果比較(a) 理論沉積數(shù)據(jù)模型; (b) 使用局部傾斜疊加獲得的局部傾角剖面; (c) 使用平面波摧毀獲得的局部傾角剖面; (d) 使用梯度平方結(jié)構(gòu)張量算法獲得局部傾角剖面.

4基于嚴(yán)格質(zhì)量監(jiān)控的自動拾取流程

圖6　基于梯度平方結(jié)構(gòu)張量算法的自動拾取流程

梯度平方結(jié)構(gòu)張量算法的高效、穩(wěn)定使得它成為自動化拾取流程中的首選.其設(shè)計思路如下：首先，對于共偏移距剖面中線性度比較大的樣點位置，利用某些準(zhǔn)則判斷其是否位于波峰位置.如果確實位于波峰位置，則讀取其局部傾角，然后換算為射線參數(shù)水平分量信息Pmx；其次，根據(jù)走時信息與炮檢點坐標(biāo)位置，計算該局部同相軸在炮道集與檢波點道集中的射線參數(shù)水平分量Prx，Psx.然后基于三者間應(yīng)滿足的定量關(guān)系進(jìn)行質(zhì)量監(jiān)控(QC)(圖6).該定量關(guān)系見4.2節(jié)詳細(xì)介紹.

4.1基于局部線性度的波峰拾取

如前所述，對共偏移距剖面進(jìn)行梯度平方結(jié)構(gòu)張量計算后得到了每個點的線性度與結(jié)構(gòu)向量.注意線性度實質(zhì)上可以理解為局部波包的同相性指標(biāo)，可以通過該指標(biāo)判斷數(shù)據(jù)點與相鄰樣點的相關(guān)度.從方便拾取的角度出發(fā)，我們同時希望立體層析所需的數(shù)據(jù)點都采集在波峰上.為了判斷波峰位置，進(jìn)一步通過已經(jīng)獲得的方向梯度(gx,gy)與梯度平方結(jié)構(gòu)張量向量V1來計算沿同相軸垂直方向的梯度d：

(6)

獲得垂直于同相軸方向的梯度信息d之后，就可以對每個點進(jìn)行波峰判斷：是否該樣點前后d的符號相反？是否該樣點處振幅大于剖面的絕對值振幅的平均值？是否該樣點處線性度(λ1-λ2)/λ1的值比較大？滿足上述標(biāo)準(zhǔn)的數(shù)據(jù)點將得到保留.圖7展示了利用上述策略獲得的數(shù)據(jù)點的分布情況.這里使用了0.7作為線性度準(zhǔn)則的門檻值，可以看到波峰判斷準(zhǔn)則結(jié)合線性度大于0.7的判別準(zhǔn)則,能夠準(zhǔn)確獲得地震剖面波峰上的數(shù)據(jù)點位置.

4.2基于物理規(guī)律的質(zhì)量監(jiān)控

4.1節(jié)中數(shù)據(jù)點的選擇策略可以認(rèn)為是基于圖像信息本身的特點實施了一次初選.但是疊前地震數(shù)據(jù)具有自己固有的物理波動規(guī)律.注意立體層析所需要的是反射或繞射信息，僅僅根據(jù)數(shù)據(jù)本身的線性度(或同相性)作為入選標(biāo)準(zhǔn)，在實際應(yīng)用中幾乎可以肯定會有不少非反射、非繞射的信息入選.這些數(shù)據(jù)點對立體層析而言缺乏物理意義，屬于應(yīng)該剔除的無效信息.

我們制定的物理QC準(zhǔn)則是：對通過初選的每個局部反射同相軸，首先統(tǒng)計其Pmx信息，再根據(jù)走時信息、炮檢點坐標(biāo)信息計算其在相關(guān)的炮道集、檢波點道集中的射線參數(shù)信息Prx與Psx.然后與共偏移距剖面中的斜率Pmx進(jìn)行比較，如果滿足Psx+Prx=Pmx(推導(dǎo)見附錄)的精度準(zhǔn)則，則為有效拾取，保存拾取結(jié)果，不滿足則剔除.這個物理QC準(zhǔn)則對實際數(shù)據(jù)的反演尤為重要.在實際應(yīng)用中，如果沒有上述物理準(zhǔn)則的嚴(yán)格QC，必然會有相當(dāng)數(shù)量的不滿足物理規(guī)律的、冗余的數(shù)據(jù)點入選.這些數(shù)據(jù)點對反演算法而言完全屬于噪聲，將會對反演結(jié)果造成負(fù)面影響.

5二維理論數(shù)據(jù)算例

二維理論數(shù)據(jù)基于圖8a所示的一個六層鹽丘模型得到.我們采用了三角剖分方式建模，基于反射射線正演的方法獲得炮數(shù)據(jù).共計模擬401炮，每炮401道，道間距10 m，炮間距20 m.圖8a顯示了原始速度模型以及某一炮對應(yīng)的射線路徑.圖8b顯示了反演所用的初始模型.這是一個常規(guī)的垂直梯度模型.

圖9顯示了正演得到的共偏移距道集、共炮道集、共檢波點道集以及基于這些記錄利用梯度平方結(jié)構(gòu)張量算法獲得的數(shù)據(jù)點信息(與局部傾角信息聯(lián)合顯示).搜索時，首先利用結(jié)構(gòu)張量算法提供的各種信息，結(jié)合波峰判斷準(zhǔn)則和線性度準(zhǔn)則(這里依然選擇線性度大于0.7的同相軸信息)獲得地震記錄上的有效數(shù)據(jù)點，圖9(a,b,c)展示了基于上述原則獲得的數(shù)據(jù)點和局部傾角信息.可以看到這種方式獲得的數(shù)據(jù)點信息是相當(dāng)準(zhǔn)確的.

圖7　梯度平方結(jié)構(gòu)張量算法獲得的波峰位置(a) 原始共偏移距數(shù)據(jù)(理論數(shù)據(jù)); (b)捕捉到的波峰位置 (線性度大于0.7得到保留).

圖8　基于射線正演的理論數(shù)據(jù)及用于反演的初始模型(a) 真實速度模型及某一炮的射線路徑; (b) 初始速度模型.

圖9　理論數(shù)據(jù)梯度平方結(jié)構(gòu)張量數(shù)據(jù)空間顯示(a) 0 m共偏移距數(shù)據(jù)點與傾角信息聯(lián)合顯示; (b) 共炮道集數(shù)據(jù)點與傾角信息聯(lián)合顯示; (c) 共檢波點道集數(shù)據(jù)點與傾角信息聯(lián)合顯示; (d) 0 m共偏移距數(shù)據(jù)點與傾角信息聯(lián)合顯示(紅點代表目標(biāo)數(shù)據(jù)點); (e) 目標(biāo)數(shù)據(jù)點在對應(yīng)的共炮道集(第1炮的201道)上與傾角信息(藍(lán)色線條)聯(lián)合顯示; (f) 目標(biāo)數(shù)據(jù)點在對應(yīng)的檢波點道集(第1個共檢波點道集的第1道)上與傾角信息(藍(lán)色線條)聯(lián)合顯示.

圖9(d,e,f)顯示了在獲得數(shù)據(jù)點信息之后，如何實施嚴(yán)格物理QC的過程.注意在0 m共偏移距剖面上用紅色標(biāo)注了一個目標(biāo)數(shù)據(jù)點(如圖9d上紅點所示)，其搜索到的Pmx=-0.054300 ms·m-1.根據(jù)道頭信息，可以找到在對應(yīng)共炮道集的相應(yīng)位置(第1炮的第201道，藍(lán)色線條代表傾角條)和對應(yīng)的共檢波點道集(第1個共檢波點道集的第1道，藍(lán)色線條代表傾角條)的位置，利用結(jié)構(gòu)張量算法可以分別計算得到Psx=-0.031207 ms·m-1和Prx=-0.023089 ms·m-1.計算Psx+Prx-Pmx=0.000004 ms·m-1，發(fā)現(xiàn)差值小于事先指定的精度門檻0.00001 ms·m-1，符合精度要求，因此這個數(shù)據(jù)點信息得到了保留.

為了測試梯度平方結(jié)構(gòu)張量算法與物理QC的作用，在這個正演記錄中還加入了一定量的隨機(jī)噪聲.可以看到無論是共偏移距記錄還是共炮點記錄或共檢波點記錄，通過線性度指標(biāo)和物理QC準(zhǔn)則的聯(lián)合作用，梯度平方結(jié)構(gòu)張量算法展示了高精度搜索同相軸局部傾角的能力.更重要的是，搜索這些傾角信息的計算成本相比局部傾斜疊加下降近兩個數(shù)量級.本測線共401炮，采用200 m偏移距的間隔，在48個核的并行計算，僅用了40 min就獲得了所有的數(shù)據(jù)點信息.如果采用局部傾斜疊加，需接近20 h才能完成同等規(guī)模的計算.即便應(yīng)用同樣的線性度標(biāo)準(zhǔn)，局部傾斜疊加獲得的數(shù)據(jù)點密度也明顯少于結(jié)構(gòu)張量算法.

對這條理論二維數(shù)據(jù)，以200 m偏移距為間隔，在不同的偏移距剖面上初選了共3萬個數(shù)據(jù)點用于反演.圖10a顯示了沒有實施物理QC的反演結(jié)果.由于數(shù)據(jù)點密度很大，也獲得了大致合理的速度結(jié)構(gòu).但是由于線性噪聲的加入，使得部分?jǐn)?shù)據(jù)點缺乏物理意義，導(dǎo)致靠近鹽丘頂部的第三、四、五、六層反射層傾角反演結(jié)果出現(xiàn)異常，鹽丘結(jié)構(gòu)發(fā)生了變異.這時加入了物理QC的過濾準(zhǔn)則.給出了非?？量痰木乳T檻(0.00001 ms·m-1).這種嚴(yán)格物理QC的加入使得合乎要求的數(shù)據(jù)點驟降到1萬個左右，但是其反演結(jié)果是相當(dāng)穩(wěn)定的.如圖10b所示，雖然在鹽丘兩側(cè)的傾角信息出現(xiàn)了一些不連續(xù)，這顯然與過于嚴(yán)格的過濾準(zhǔn)則有關(guān).但剩下的數(shù)據(jù)點質(zhì)量很高，反演結(jié)果中沒有出現(xiàn)異常的、違反地質(zhì)規(guī)律的現(xiàn)象.

圖10　理論數(shù)據(jù)反演結(jié)果對比(a) 未實施物理QC； (b) 實施物理QC后.

圖11　共偏移距剖面上基于結(jié)構(gòu)張量算法實現(xiàn)數(shù)據(jù)空間提取(a) 原始數(shù)據(jù)的300 m偏移距的共偏移距道集； (b) 僅利用數(shù)據(jù)信息自動拾取出的位置，綠色小點處為波峰； (c) 物理QC準(zhǔn)則應(yīng)用于圖b之后留下的數(shù)據(jù)點，綠色小點為拾取點，藍(lán)色斜線為該點處局部傾角.

6南海某二維深水實際數(shù)據(jù)反演

6.1數(shù)據(jù)拾取

選擇南海某二維深水測線進(jìn)行立體層析反演處理.整條測線長150 km，共2935炮，最大偏移距8275 m，炮間距50 m，檢波點間距25 m.整條測線包括淺水、深水兩個部分，海水深度變化從200 m到2000 m.反演之前的數(shù)據(jù)做了精細(xì)的前處理：如壓制與海面有關(guān)的多次波、實施預(yù)測反褶積提高分辨率等等.在數(shù)據(jù)拾取之前我們進(jìn)行了10～45 Hz的窄帶濾波器并實施短時窗的自動增益以保證同相軸可以識別，并且切除了直達(dá)波以上的數(shù)據(jù)來減少拾取工作量，使得預(yù)處理后的數(shù)據(jù)有較高的信噪比來保證拾取的可靠性.

由于梯度平方結(jié)構(gòu)張量算法的高效，使得本次反演的數(shù)據(jù)空間提取效率大為提高.反射波拾取范圍劃定在100～6200 m偏移距，時間范圍在8 s以內(nèi).每隔200 m偏移距進(jìn)行一次拾取.圖11展示了在南海某二維深水?dāng)?shù)據(jù)的300 m共偏移距剖面上，應(yīng)用物理QC之后的數(shù)據(jù)點前后的變化.可以看到，在應(yīng)用物理QC過濾準(zhǔn)則后，一部分拾取點被刪除了.由于采取了比較嚴(yán)格的過濾條件，圖11c中展示的數(shù)據(jù)點信息都是非?？煽康?

立體層析的目標(biāo)泛函是基于擬合數(shù)據(jù)域信息建立的，因此拾取數(shù)據(jù)的精度十分重要，這是反演能否成功的最基本保證.結(jié)合梯度平方結(jié)構(gòu)張量算法，實現(xiàn)了針對二維實際數(shù)據(jù)的高密度立體層析反演.在MPI多進(jìn)程并行技術(shù)支持下，采用120進(jìn)程并行計算2.5 h就全自動完成數(shù)據(jù)空間準(zhǔn)備工作.得到了關(guān)于整條測線的14萬個數(shù)據(jù)點的拾取數(shù)據(jù)量.如果使用局部傾斜疊加的方法得到同等的數(shù)據(jù)空間準(zhǔn)備需要6～8天，而且還達(dá)不到目前結(jié)構(gòu)張量算法所能提供的數(shù)據(jù)密度.梯度平方結(jié)構(gòu)張量算法的應(yīng)用意味著在更短的生產(chǎn)周期內(nèi)得到了更密集均勻的數(shù)據(jù)拾取覆蓋率.這對于立體層析反演的實際應(yīng)用能力無疑具有積極意義.

6.2反演策略與結(jié)果分析

反演模型的規(guī)模為橫向150 km，深度為10 km，為了壓縮模型空間，采用了B樣條表達(dá).采用的B樣條系數(shù)的橫向間距和縱向間距均為200 m.即便采用了B樣條表達(dá)，由于模型空間很大，射線覆蓋往往不均勻，層析矩陣依然是稀疏、欠定、病態(tài)的.為了改善層析矩陣的稀疏性，降低規(guī)則化對反演結(jié)果的平滑效應(yīng)，高密度的數(shù)據(jù)空間準(zhǔn)備是非常重要的.

初始模型使用最普通的梯度模型.在初始模型上通過疊前深度偏移可以得到初始深度成像剖面，由于初始模型中對應(yīng)于海水的速度是準(zhǔn)確的，因此在初始深度成像剖面上可以拾取出海底高程信息.因此在反演過程中無需計算數(shù)據(jù)空間對海水速度的Frechet導(dǎo)數(shù).同時針對海水部分給出很強(qiáng)的阻尼規(guī)則化參數(shù)，使得海水層內(nèi)速度保持穩(wěn)定.此外，在規(guī)則化參數(shù)設(shè)置方面，在沒有射線或低射線覆蓋區(qū)域設(shè)置了較高的阻尼系數(shù)，在高射線覆蓋區(qū)域則設(shè)置較低的阻尼系數(shù)，以期在反演的穩(wěn)定與精度之間獲得一種平衡.

為求解這個大規(guī)模稀疏立體層析矩陣方程，我們采用了多波前多線程的稀疏矩陣QR分解算法(Davis, 2011)來求解方程.在MPI環(huán)境下進(jìn)行多進(jìn)程并行計算射線正演以及Frechet導(dǎo)數(shù)求取等.在48個進(jìn)程并行計算的情況下，在4 h內(nèi)即完成24次迭代，得到了滿意的反演結(jié)果.

圖12a顯示了立體層析在24次迭代更新之后的反演結(jié)果.注意其中的藍(lán)色線條為模型空間中的局部反射傾角信息，將其貼合到反演速度上聯(lián)合顯示的是立體層析中展示反演成果最常用的一種方式.由于實際數(shù)據(jù)的信噪比在不同的構(gòu)造區(qū)域內(nèi)有比較大的差別，因此拾取數(shù)據(jù)在信噪比高的區(qū)域比較密集，在信噪比低的區(qū)域則相對稀疏.圖12b將圖12a中20～60 km段的剖面放大顯示.可以看到高密度的數(shù)據(jù)點拾取所對應(yīng)的反射點位置在反演過程中自動聚集形成了有效的層位信息，和實際的成像層位位置符合的非常好.需要強(qiáng)調(diào)，類似的高密度反演結(jié)果是利用局部傾斜疊加提取數(shù)據(jù)空間時無法做到的.圖12c展示了基于該模型的全測線疊前深度成像結(jié)果，可以看到整體構(gòu)造形態(tài)得到了精細(xì)的刻畫，有利于后續(xù)的地質(zhì)解釋.

圖12　通過結(jié)構(gòu)張量算法獲得的數(shù)據(jù)空間反演得到的速度模型(a) 完整測線速度模型與模型空間傾角; (b) 局部展示20000～60000 m范圍的速度反演結(jié)果與模型空間傾角; (c) 基于完整測線速度模型偏移得到的疊前深度偏移結(jié)果.

圖13　測線局部(110～150 km)的傾斜疊加數(shù)據(jù)空間與梯度平方結(jié)構(gòu)張量數(shù)據(jù)空間的反演與成像結(jié)果對比(a) 局部傾斜疊加在800 m單偏移距剖面內(nèi)獲得的數(shù)據(jù)空間; (b) 梯度平方結(jié)構(gòu)張量在800 m單偏移距剖面獲得的數(shù)據(jù)空間; (c) 局部傾斜疊加數(shù)據(jù)空間反演的速度模型(110～150 km); (d) 梯度平方結(jié)構(gòu)張量數(shù)據(jù)空間反演的速度模型(110～150 km); (e) 用圖c所示模型實施PSDM后，在模型130 km處的若干CIG顯示; (f) 用圖d所示模型實施PSDM后，在模型130 km處的若干CIG顯示; (g) 局部傾斜疊加數(shù)據(jù)空間的疊前深度偏移結(jié)果(110～150 km); (h) 梯度平方結(jié)構(gòu)張量數(shù)據(jù)空間的疊前深度偏移結(jié)果(110～150 km).

由于局部傾斜疊加準(zhǔn)備全測線的數(shù)據(jù)空間過于耗時，因此這里無法提供全測線的、兩種算法的反演結(jié)果對比.圖13展示了在測線末尾抽取了900炮記錄后在其中分選得到的800 m偏移距剖面.其中圖13a和圖13b分別顯示了局部傾斜疊加與梯度平方結(jié)構(gòu)張量算法提取得到的數(shù)據(jù)空間分步對比.其中黃色線條代表射線覆蓋密度，可以看出后者的數(shù)據(jù)點密度與均勻度均明顯超過前者，形成了更為均勻的射線覆蓋強(qiáng)度，這對于提高反演的穩(wěn)定性顯然是有利的.圖13c,圖13d展示了基于這兩種數(shù)據(jù)空間，最終反演得到的速度模型.顯然結(jié)構(gòu)張量算法提供了更多的有效數(shù)據(jù)點，使得反演得到的結(jié)構(gòu)更為精細(xì).這種差異也體現(xiàn)在圖13e,13f所示的兩張疊前深度成像點道集以及圖13g,13h所示疊前深度偏移剖面上.可以看出基于結(jié)構(gòu)張量算法數(shù)據(jù)空間得到的反演結(jié)果在疊前深度成像之后，無論對于能量的聚焦(CDP 3800-CDP 4000)、或?qū)τ诨?CDP 5200-CDP 5500)的歸位、或?qū)τ跇?gòu)造形態(tài)的正確刻畫(CDP 4100-CDP 4400)，其精度都超過了傾斜疊加數(shù)據(jù)空間反演結(jié)果提供的深度成像剖面.

6結(jié)論與討論

本文將圖像處理中用于邊緣檢測的梯度平方結(jié)構(gòu)張量算法引入立體層析數(shù)據(jù)空間的準(zhǔn)備.實踐表明，該算法能夠高效率、高密度的獲得局部同相軸的傾角信息，這是實施高密度立體層析反演的基本保證.同時，作者提出一種基于物理規(guī)律的質(zhì)量監(jiān)控方法來過濾非物理的冗余數(shù)據(jù)點.由于梯度平方結(jié)構(gòu)張量算法以及合乎物理規(guī)律的嚴(yán)格QC的引入，使得自動化的高密度立體層析反演成為可能.基于二維理論數(shù)據(jù)與南海某二維深水地震數(shù)據(jù)的測試表明，與常規(guī)局部傾角估計算法準(zhǔn)備數(shù)據(jù)空間的反演結(jié)果相比，高密度的有效數(shù)據(jù)點的使用降低了層析矩陣的病態(tài)性、提高了反演的精度，改善了疊前深度成像的品質(zhì).上述特點使得立體層析反演具備了應(yīng)對實際數(shù)據(jù)的能力.同時，由于梯度平方結(jié)構(gòu)張量算法很容易推廣到三維，因此基于三維梯度平方結(jié)構(gòu)張量算法的三維立體層析反演的全自動實現(xiàn)也具有很好的應(yīng)用潛力.

致謝作者感謝科羅拉多礦業(yè)學(xué)院伍新民博士研究生在梯度平方結(jié)構(gòu)張量算法方面的有益討論；作者同時感謝國家自然科學(xué)基金面上項目(41274117、41574098)、國家科技重大專項子課題(2011ZX05025-001-03)、以及海洋地質(zhì)國家重點實驗室自主課題(MG20130304)對于這項研究的支持.

附錄A質(zhì)量監(jiān)控公式推導(dǎo)

對原始數(shù)據(jù)集，可看做炮域數(shù)據(jù)和共偏移距域數(shù)據(jù)，關(guān)系如下：

(A1)

對公式(A1)，對炮點、檢波點求導(dǎo)，得

(A2)

(A3)

化簡上式，得到

(A4)

(A5)

根據(jù)半偏移距h的定義h=(s-r)/2和中點m的定義m=(s+r)/2，我們有

(A6)

即

(A7)

(A8)

得到

(A9)

PH=PS-PR.

(A10)

References

Billette F, Lambaré G. 1998. Velocity macro-model estimation from seismic reflection data by stereotomography.GeophysicalJournalInternational, 135(2): 671-690.

Davis T A. 2011. Algorithm 915, Suite Sparse QR: Multifrontal multithreaded rank-revealing sparse QR factorization.ACMTransactionsonMathematicalSoftware(TOMS), 38: 8.

Fomel S. 2002. Applications of plane-wave destruction filters.Geophysics, 67(6): 1946-1960.

Hale D. 2009. Structure-oriented smoothing and semblance. CWP Report 635, Center Wave Phenomena.

Ni Yao, Yang Kai. 2012. Slope tomography assisted by finite-frequency sensitivity kernel. ∥ 82th SEG Expanded Abstract.

Lambare. 2008. Stereotomgraphy.Geophysics, 73(5): VE25-VE34.

Wu X M, Hale D. 2013. Extracting horizons and sequence boundaries from 3D seismic images. ∥ 83th SEG Expanded Abstract.Ottolini R. 1983. Signal/noise separation in dip space.StanfordExplorationProject, 37: 143-149.

Van Vliet L J. Verbeek P W. 1995. Estimators for orientation and anisotropy in digitized images. ∥Proc. of the First Annual Conf. of the Advanced School for Computing and Imaging (ASCI′95), Heijen, NL, May 16-18, 442-450.

Van Vliet L J, Young I T, Verbeek P W. 1998. Recursive Gaussian derivative filters. ∥ Proceedings of the Fourteenth International Conference on Pattern Recognition. Brisbane, Qld: IEEE, 509-514.

Li Z W, Yang K, Ni Y, et al. 2014. Migration velocity analysis with stereo-tomography inversion.GeophysicalProspectingforPetroleum(in Chinese), 53(4): 444-452.

Ni Y, Yang K, Chen B S. 2013. Stereotomography inversion method theory and application testing.GeophysicalProspectingforPetroleum(in Chinese), 52(2): 430-436.Ren L J, Sun X D, Li Z C. 2014. The stereotomography based on eigen-wave attributes. ∥ The expanded abstract of the CPS/SEG Beijing 2014 International Geophysical Conference (in Chinese).

附中文參考文獻(xiàn)

李振偉, 楊鍇, 倪瑤等. 2014. 基于立體層析反演的偏移速度建模應(yīng)用研究. 石油物探, 53 (4): 444-452.

倪瑤, 楊鍇, 陳寶書. 2013. 立體層析反演方法理論分析與應(yīng)用測試. 石油物探, 52(2): 121-130.

任麗娟, 孫小東, 李振春. 2014. 基于特征波屬性參數(shù)的立體層析速度反演方法研究. ∥ CPS/SEG北京2014國際地球物理會議摘要.

(本文編輯汪海英)

基金項目國家科技重大專項子課題(2011ZX05025-001-03), 國家自然科學(xué)基金面上項目(41274117、41574098)以及海洋地質(zhì)國家重點實驗室自主課題(MG20130304) 資助.

作者簡介王宇翔，男，1988年生，畢業(yè)于同濟(jì)大學(xué)地球物理系，理學(xué)碩士，研究方向為偏移速度分析與立體層析反演. E-mail：86wyx@#edu.cn*通訊作者楊鍇，男, 1972年生,理學(xué)博士,教授,主要研究方向為地震波反演與成像. E-mail:yang_kai@#edu.cn

doi:10.6038/cjg20160122 中圖分類號P631

收稿日期2014-09-23，2015-09-27收修定稿

A high-density stereo-tomography method based on the gradient square structure tensors algorithm

WANG Yu-Xiang1,2, YANG Kai1*, YANG Xiao-Chun3, XUE Dong3, CHEN Bao-Shu3

1StateKeyLaboratoryofMarineGeology,TongjiUniversity,Shanghai200092,China2Broadgeo,Inc,Houston,TX770994 3CNOOCResearchCenter,Beijing100027,China

AbstractThe efficiency and effectiveness of the stereotomography is highly dependent on the quality of its data space, the so-called kinematic invariant. The structure tensor is a very robust tool for the slope estimation. In this paper, we present a highly efficient high-density kinematic invariant extraction method based on structure tensors. Compared with the conventional slope search methods, the presented technique can improve the computational efficiency by one or two orders of magnitudes which will greatly enhance the applicability of the stereomography.

KeywordsStereotomography; Data space extraction; Gradient square based structure tensor algorithm; Migration velocity analysis

王宇翔， 楊鍇， 楊小椿等. 2016. 基于梯度平方結(jié)構(gòu)張量算法的高密度二維立體層析反演.地球物理學(xué)報，59(1):263-276,doi:10.6038/cjg20160122.

Wang Y X, Yang K, Yang X C, et al. 2016. A high-density stereo-tomography method based on the gradient square structure tensors algorithm.ChineseJ.Geophys. (in Chinese),59(1):263-276,doi:10.6038/cjg20160122.