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

王宇翔， 楊鍇， 楊小椿， 薛冬， 陳寶書(shū)

1 同濟(jì)大學(xué)海洋地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 上?！?00092 2 Broadgeo, Inc, Houston, TX 770994 3 中海石油研究中心, 北京　100027

?

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

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

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

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

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

1引言

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

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

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

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

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

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

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

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

(1)

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

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

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

F=

(2)

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

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

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

圖2　某海洋二維測(cè)線原始炮記錄

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

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

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

(3)

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

(4)

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

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

(5)

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

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

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

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

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

表1　三種算法的計(jì)算時(shí)間對(duì)比(單位：s)

從表1可見(jiàn)梯度平方結(jié)構(gòu)張量算法的計(jì)算效率高于平面波摧毀近一個(gè)數(shù)量級(jí)，高于局部?jī)A斜疊加達(dá)兩個(gè)數(shù)量級(jí).注意梯度平方結(jié)構(gòu)張量算法獲得的局部?jī)A角信息是一個(gè)鄰域內(nèi)的平均效應(yīng).因此其變化比較光滑連續(xù)，不像局部?jī)A斜疊加獲得的局部?jī)A角信息似有突變.但是這并不表示其分辨率低于局部?jī)A斜疊加，相反其傾角估計(jì)的穩(wěn)定性恰恰來(lái)自這種鄰域平均效應(yīng).需要指出，即便考慮了這個(gè)鄰域平均效應(yīng)，它在計(jì)算成本方面相對(duì)于另外兩種常規(guī)算法依然存在巨大優(yōu)勢(shì).這使得它成為立體層析反演實(shí)際應(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　三種傾角估計(jì)算法結(jié)果比較(a) 理論沉積數(shù)據(jù)模型; (b) 使用局部?jī)A斜疊加獲得的局部?jī)A角剖面; (c) 使用平面波摧毀獲得的局部?jī)A角剖面; (d) 使用梯度平方結(jié)構(gòu)張量算法獲得局部?jī)A角剖面.

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

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

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

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

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

(6)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

6結(jié)論與討論

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

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

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

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

(A1)

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

(A2)

(A3)

化簡(jiǎn)上式，得到

(A4)

(A5)

根據(jù)半偏移距h的定義h=(s-r)/2和中點(diǎn)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.

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

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

(本文編輯汪海英)

基金項(xiàng)目國(guó)家科技重大專(zhuān)項(xiàng)子課題(2011ZX05025-001-03), 國(guó)家自然科學(xué)基金面上項(xiàng)目(41274117、41574098)以及海洋地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室自主課題(MG20130304) 資助.

作者簡(jiǎn)介王宇翔，男，1988年生，畢業(yè)于同濟(jì)大學(xué)地球物理系，理學(xué)碩士，研究方向?yàn)槠扑俣确治雠c立體層析反演. E-mail：86wyx@#edu.cn*通訊作者楊鍇，男, 1972年生,理學(xué)博士,教授,主要研究方向?yàn)榈卣鸩ǚ囱菖c成像. E-mail:yang_kai@#edu.cn

doi:10.6038/cjg20160122 中圖分類(lèi)號(hào)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é)報(bào)，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.