聯(lián)合結(jié)構(gòu)張量與運動學(xué)反偏移的立體層析數(shù)據(jù)空間提取與反演策略研究Ⅰ:理論

楊 鍇,熊 凱,王宇翔,汪小將

(1.同濟大學(xué)反射地震學(xué)科組,上海200092;2.阿美遠(yuǎn)東(北京)商業(yè)服務(wù)有限公司國際研發(fā)中心,北京100102;3.中海石油研究中心,北京100027)

聯(lián)合結(jié)構(gòu)張量與運動學(xué)反偏移的立體層析數(shù)據(jù)空間提取與反演策略研究Ⅰ:理論

楊 鍇1,熊 凱1,王宇翔2,汪小將3

(1.同濟大學(xué)反射地震學(xué)科組,上海200092;2.阿美遠(yuǎn)東(北京)商業(yè)服務(wù)有限公司國際研發(fā)中心,北京100102;3.中海石油研究中心,北京100027)

提出了一種聯(lián)合結(jié)構(gòu)張量與運動學(xué)反偏移的立體層析數(shù)據(jù)空間提取與反演策略。首先在數(shù)據(jù)域搜索得到初始數(shù)據(jù)空間,可獲得初始反演結(jié)果與初始成像數(shù)據(jù)體,然后基于初始成像數(shù)據(jù)體內(nèi)的共偏移距成像剖面人工拾取出比較可靠的數(shù)據(jù)點位置,在這些位置搜索構(gòu)造傾角與剩余曲率,最后利用運動學(xué)反偏移將上述成像域運動學(xué)信息反偏移到數(shù)據(jù)空間并實施校正,獲得更為可靠與均勻的立體層析數(shù)據(jù)空間。無論是第一步的數(shù)據(jù)域搜索還是第二步的成像域搜索都是基于高效率、高精度的結(jié)構(gòu)張量算法的實現(xiàn),可認(rèn)為是一種聯(lián)合結(jié)構(gòu)張量與運動學(xué)反偏移的兩步法立體層析數(shù)據(jù)空間提取與反演策略。同時還探討了當(dāng)射線參數(shù)水平分量信息引入層析反演的數(shù)據(jù)空間后,實施反演算法的各種可能性。理論數(shù)據(jù)證實了上述策略的可靠性,為后續(xù)的實踐奠定了理論基礎(chǔ)。

立體層析;結(jié)構(gòu)張量;運動學(xué)反偏移;數(shù)據(jù)空間提取

立體層析反演是層析反演類速度建模方法中極具特色的一種方法。立體層析將一個局部相干同相軸在炮道集與檢波點道集內(nèi)的射線參數(shù)水平分量(Psx,Prx,以下簡稱地表觀測P參數(shù))納入到數(shù)據(jù)空間之中,不僅使得數(shù)據(jù)空間的準(zhǔn)備相比傳統(tǒng)反射層析更為簡便,也使得立體層析成為運動學(xué)層析反演方法中唯一一種可以同時反演速度結(jié)構(gòu)與反射點位置的方法[1-4]。李振偉等[5]將該方法首次應(yīng)用于中國南海實際二維地震數(shù)據(jù)處理,獲得了可靠的反演結(jié)果。反演得到的速度模型能夠滿足疊前深度偏移的需要。王宇翔等[6]利用高效率的結(jié)構(gòu)張量實現(xiàn)了在疊前數(shù)據(jù)域提取P參數(shù)的快速算法,使得基于結(jié)構(gòu)張量的二維高密度立體層析反演成為現(xiàn)實。楊鍇等[7]基于非降階漢密爾頓在直角坐標(biāo)系下實現(xiàn)了三維立體層析反演。XING等[8]將該算法首次用于南海某三維實際數(shù)據(jù)反演,同樣利用結(jié)構(gòu)張量獲取三維立體層析數(shù)據(jù)空間,獲得了理想的速度建模效果。

在實踐中我們體會到,立體層析數(shù)據(jù)空間的正確提取是保證其反演成功的最關(guān)鍵因素。在BILLETTE等[3]和李振偉等[5]的工作中,立體層析數(shù)據(jù)空間提取通過傾斜疊加能量譜的交互分析拾取獲得,這是一種穩(wěn)健但效率較低的數(shù)據(jù)空間提取方式。該提取方式在二維數(shù)據(jù)中還可以使用,但對于三維數(shù)據(jù)其應(yīng)用潛力十分有限。此外,BILLETTE等[1]認(rèn)為單次反射與單次繞射都可以用于立體層析反演。這個觀點值得商榷。我們認(rèn)為,繞射波不適用于基于運動學(xué)信息的絕大多數(shù)層析成像算法,至少對立體層析而言,真正有用的僅是一次反射波的運動學(xué)信息。這一觀點在本文的數(shù)值算例中得到了證實。因此,在提取立體層析數(shù)據(jù)空間時應(yīng)盡量避開非一次反射波(如繞射波、多次波和側(cè)面反射)的干擾,保證提取到的運動學(xué)信息都來自一次反射波。

然而,無論是基于傾斜疊加能量譜的交互拾取還是基于結(jié)構(gòu)張量的自動拾取,都無法完全避開上述非一次反射波的干擾,例如:基于傾斜疊加能量譜的交互拾取對多次波的識別或許有效,但是對繞射波與側(cè)面波則無能為力;基于結(jié)構(gòu)張量的自動拾取盡管效率十分令人滿意,但對于繞射波、多次波和側(cè)面反射均無法自動識別?？傊?無論使用何種快速算法,直接從疊前數(shù)據(jù)體中提取數(shù)據(jù)空間必然遇到無法辨識多次波、繞射波以及側(cè)面反射的問題;其次,由于信噪比和能量強度的原因,自動拾取往往導(dǎo)致淺層的數(shù)據(jù)點獲取密度很高,但是深層的數(shù)據(jù)點信息則大幅減少,無法保證數(shù)據(jù)空間的均勻性。這對于中深部目標(biāo)體的反演相當(dāng)不利。

CHAURIS等[9]揭示了偏移速度分析中地下成像域運動學(xué)信息與地表觀測P參數(shù)之間的定量關(guān)系,給出了在速度不正確的情況下,基于疊前成像空間中提取的運動學(xué)屬性參數(shù)(它們分別是共偏移距成像剖面上的構(gòu)造傾角與共偏移距成像點道集內(nèi)的剩余曲率),通過特定的修正公式,計算得到正確的地表觀測P參數(shù)。基于上述思想,GUILLAUME等[10]提出了通過“運動學(xué)反偏移+校正”獲得立體層析數(shù)據(jù)空間的流程,這個工作的意義在于,即便對于不正確的初始速度獲得的不正確的成像結(jié)果,通過高效的運動學(xué)信息提取和有效的校正手段,依然可以獲得準(zhǔn)確的地表觀測P參數(shù)。這個流程的應(yīng)用價值在于它提供了一個在成像域確定數(shù)據(jù)點位置的選擇。也就是說,我們可以根據(jù)初始成像道集,結(jié)合自己的地質(zhì)認(rèn)識,避開原始數(shù)據(jù)上可能存在的繞射、多次波、側(cè)面波等各種非一次反射波的干擾,優(yōu)選出最重要、最關(guān)鍵的數(shù)據(jù)點,服務(wù)于后續(xù)的斜率層析或者其它需要P參數(shù)的各種快速成像算法[11-12]。需要指出的是,基于疊前數(shù)據(jù)體的高效結(jié)構(gòu)張量自動拾取與利用反偏移重建數(shù)據(jù)空間并不矛盾,這恰恰構(gòu)成了立體層析在實際應(yīng)用過程的兩步法:先通過結(jié)構(gòu)張量算法在數(shù)據(jù)域自動拾取獲得初始數(shù)據(jù)空間,基于初始數(shù)據(jù)空間反演得到一個很好的初始速度模型;然后在觀察初始成像數(shù)據(jù)體后,通過人工識別去除與繞射波、多次波、側(cè)面反射有關(guān)的數(shù)據(jù)點,提煉出更為合理的數(shù)據(jù)空間,使得立體層析反演精度得到進一步提升。

值得注意的是,當(dāng)射線參數(shù)水平分量和射線出射到地表的坐標(biāo)也被納入到數(shù)據(jù)空間后,會導(dǎo)致層析反演有多種可能的實現(xiàn)路徑。本文首先分析了立體層析的理論基礎(chǔ)——反射地震學(xué)中重要的“走時一一映射”條件,針對立體層析數(shù)據(jù)空間的特點分析了各種可能的實現(xiàn)策略,并指出本文的實現(xiàn)策略與所有這些可能策略的關(guān)系;然后詳細(xì)介紹結(jié)構(gòu)張量以及“運動學(xué)反偏移+校正”獲得立體層析數(shù)據(jù)空間的方法原理;最后利用一個典型的深水理論數(shù)據(jù)驗證了上述算法設(shè)計與實現(xiàn)策略的正確性,證實其可以獲得更可靠的數(shù)據(jù)空間,獲得更合理的反演結(jié)果。

1 立體層析的數(shù)據(jù)空間、模型空間以及“走時一一映射條件”

首先以二維情況為例介紹立體層析反演的模型空間與數(shù)據(jù)空間。如圖1a所示,反射射線可以認(rèn)為是從反射點出發(fā),分別以不同的散射角θs,θr向炮點位置和檢波點位置出射兩根射線構(gòu)成的。當(dāng)射線到達(dá)觀測面時,兩條透射射線的出射方向、走時之和以及炮點和檢波點位置構(gòu)成立體層析數(shù)據(jù)空間。二維立體層析的數(shù)據(jù)分量(d)和模型分量(m)為:

式中:s,r分別為炮點和檢波點橫坐標(biāo);t為雙程走時;x0,z0分別為反射點的橫、縱坐標(biāo);θs,θr分別為炮點和檢波點兩側(cè)的散射角;v(x,z)為模型空間中的速度信息。

公式(1a)表達(dá)了拾取到的ndata個數(shù)據(jù)點,每個數(shù)據(jù)點對應(yīng)于一個反射過程。公式(1b)則表達(dá)了模型空間信息。

圖1 二維立體層析的模型與數(shù)據(jù)分量a 立體層析的模型與數(shù)據(jù)分量; b 同相軸斜率與慢度矢量之間的關(guān)系

立體層析數(shù)據(jù)空間中最重要的數(shù)據(jù)分量是射線參數(shù)水平分量(或慢度矢量水平分量)。圖1中所示Psx可在共檢波點道集內(nèi)搜索局部同相軸斜率得到,Prx可在共炮點道集內(nèi)搜索局部同相軸斜率得到。當(dāng)考察在圖1b所顯示的炮記錄中某一局部同相軸時,由幾何關(guān)系及射線參數(shù)水平分量的定義可得:

(2)

式中:Pslope為局部同相軸的斜率;Prx為檢波點處射線參數(shù)水平分量。公式(2)表明共炮記錄內(nèi)同相軸斜率必對應(yīng)于檢波點處慢度矢量水平分量Prx。由炮檢互易原理可知,共檢波記錄內(nèi)同相軸斜率必對應(yīng)于炮點處慢度矢量的水平分量Psx。只要能夠正確提取P參數(shù),其對應(yīng)的旅行時和坐標(biāo)信息均可輕易得到。

對應(yīng)于二維立體層析的Fréchet偏導(dǎo)數(shù)矩陣為:

(3)

式中:σ為層析矩陣每一行物理量的數(shù)量級不同設(shè)置的均衡系數(shù);Sx，Rx分別為地表出射炮點和檢波點位置;T為總旅行時?？梢钥吹?常規(guī)的旅行時層析由于其模型空間僅包含速度v,其數(shù)據(jù)空間僅包含走時,因此,其Fréchet偏導(dǎo)數(shù)矩陣僅相當(dāng)于公式(3)中右上角的子矩陣。立體層析Fréchet偏導(dǎo)數(shù)矩陣通過在背景介質(zhì)中進行動力學(xué)射線追蹤,然后利用射線擾動理論[13]建立數(shù)據(jù)空間殘差各分量與模型空間殘差各分量之間的線性關(guān)系后計算得到。

注意,圖1a表達(dá)了反射地震學(xué)成像理論中的一個非常重要的成像條件:走時一一映射條件。事實上走時一一映射條件是立體層析與共角度域成像共同的理論基礎(chǔ)。該條件由TEN KROODE等[13]通過傅里葉積分算子理論予以嚴(yán)格證明。其物理意義是:地震數(shù)據(jù)中任意一個局部相干的同相軸,可以由炮點坐標(biāo)、檢波點坐標(biāo)、同相軸在兩坐標(biāo)處的局部斜率以及同相軸的局部走時來描述,它嚴(yán)格對應(yīng)于從反射點出發(fā)、到達(dá)炮檢點的一個射線對。在速度模型正確的前提下,記錄到的一個局部相干的地震數(shù)據(jù)同相軸將成像于地下唯一的一個深度點。在反射地震學(xué)觀測中,違反這個條件的幾率幾乎為零。因此該條件有廣泛的適用范圍。XU等[14]以此條件為準(zhǔn)則審視地震數(shù)據(jù)中的各種道集,發(fā)現(xiàn)共炮道集、共偏移距道集均違反此條件的要求,只有共角度道集才滿足這個條件。在獲得這個認(rèn)識后,最終在最小二乘Kirchhoff真振幅反演成像框架下導(dǎo)出了輸出真振幅角度域成像道集的成像條件。

走時一一映射條件的重要性不僅在于正確的理解角度域成像道集,還在于它同時揭示了對層析成像而言一些非常重要的事實:①既然炮檢點處出射方向確定的一個射線對,唯一對應(yīng)于地下的一個反射點,那么在層析反演中,就可以將反射問題轉(zhuǎn)化為透射問題,從而回避傳統(tǒng)反射層析中反射層深度與速度耦合性;②方向信息可以更好地約束射線對,減少其模糊性。方向信息明確的射線對可以給速度分析提供額外的約束信息,使得在復(fù)雜地質(zhì)結(jié)構(gòu)下的層析成像成為可能。正是這些重要的認(rèn)識促使BILLETTE等[1]提出了一種充分利用地震波場所有運動學(xué)信息的全新層析方法,即立體層析反演。

體現(xiàn)這些優(yōu)勢的關(guān)鍵在于同相軸的斜率信息(即P參數(shù))必須進入層析反演的數(shù)據(jù)空間。通過一個簡單的概念實驗即可說明同相軸斜率對速度模型的約束作用。如圖1b所示,當(dāng)?shù)乇砼邳c位置(S)、檢波點位置(R)、炮點處的慢度矢量水平分量(Psx)和檢波點處的慢度矢量水平分量(Prx)已知時,由炮點至檢波點的反射過程S→X(反射點)→R唯一確定。分別根據(jù)Psx,Prx于炮點、檢波點反向朝地下出射射線至兩條射線相交時停止,當(dāng)速度模型正確時,兩條射線相交在正確的反射點,反之當(dāng)速度模型不正確時,兩條射線的交點為錯誤的反射點。那么如何判斷兩條射線是否相交在正確的位置呢？走時信息提供了另外的約束條件,即射線在正確的速度模型上傳播時,兩條射線的走時之和恰與反射走時相等,反之則會出現(xiàn)殘差。

顯然,共炮點、共檢波點道集上同相軸的斜率為層析反演提供了除走時外的新的運動學(xué)約束信息,即數(shù)據(jù)空間可以由更加豐富的運動學(xué)波場屬性參數(shù)來刻畫,可以認(rèn)為立體層析反演的優(yōu)勢主要得益于同相軸斜率的引入,使得數(shù)據(jù)空間和模型空間都得到了合理的擴展。

2 基于立體層析數(shù)據(jù)空間的反演算法設(shè)計與實現(xiàn)策略分析

從上一節(jié)中的概念出發(fā),即可對基于P參數(shù)的層析反演進行歸納和總結(jié)。當(dāng)我們除了擁有旅行時和坐標(biāo)信息,同時擁有炮點和檢波點處的P參數(shù)信息時,反演方案將會有很多。從殘差所屬域的角度來看:數(shù)據(jù)域有旅行時差、P參數(shù)差、炮檢點坐標(biāo)差;成像域殘差則體現(xiàn)為剩余曲率(RMO)不為零,即所謂的不聚焦現(xiàn)象。數(shù)據(jù)殘差的多樣性決定了反演時可以有多種策略供選擇。注意無論使用哪種算法,其數(shù)據(jù)殘差各分量與模型殘差各分量之間的一階擾動關(guān)系,都可以通過文獻(xiàn)[13]中的射線擾動理論來建立。圖2顯示了聯(lián)合P參數(shù)與旅行時的各種層析反演策略。

第1種方案如圖2a所示:先根據(jù)炮點和檢波點的P參數(shù)信息進行運動學(xué)偏移,即利用Psx和Prx分別從炮檢點往下追射線,當(dāng)ts+tr=Tobs時射線追蹤停止(Tobs為初至旅行時)。由于速度模型錯誤,射線追蹤停止時必然發(fā)生欠聚焦現(xiàn)象,即兩個射線終點水平方向有坐標(biāo)差Xerr。這時即可利用坐標(biāo)差Xerr來反演速度。這種方法的實質(zhì)是,用P參數(shù)和走時一起成像,利用成像域欠聚焦信息反演背景速度。這正是SWORD[16]最早的基于同相軸斜率信息反演背景速度的思路。

第2種方案如圖2b所示:先根據(jù)炮點和檢波點處的P參數(shù)信息進行反向射線追蹤,直到兩支射線在地下相交。相當(dāng)于完全用P參數(shù)成像(即利用P參數(shù)確定反射點的位置),因為速度還不準(zhǔn)確,這時ts+tr必然不等于Tobs,于是就可以用旅行時差反演速度。這種方法的實質(zhì)是,用P參數(shù)成像,而后利用數(shù)據(jù)域旅行時殘差信息反演背景速度。熊凱[17]對這種算法進行了測試,在簡單層狀模型時反演效果較為理想,但是,當(dāng)模型趨于復(fù)雜時必須排除兩根射線夾角過小的情況,否則會出現(xiàn)由于旅行時差計算不準(zhǔn)而導(dǎo)致速度更新量發(fā)生異常的情況。

第4種方案如圖2d所示:以初始速度的運動學(xué)偏移(即向地下出射射線)初始化地下反射點位置,從地下反射點往地表出射射線,找到炮點和檢波點,建立δPsx,δPrx,δT與δv,δx0,δz0之間的擾動關(guān)系,從而更新背景速度和反射點的位置。這種方法的特點是,運動學(xué)偏移初始化反射點后確定地下反射點與地表的局部同相軸的關(guān)系,再反偏移到地表得到數(shù)據(jù)域殘差δPsx,δPrx,δT來更新背景速度和反射點位置。

第5種方案如圖1a所示:通過運動學(xué)偏移初始化地下反射的位置和出射角,根據(jù)地下出射角從地下點往地表出射射線,不追求找到正確的炮點與檢波點坐標(biāo),因此射線出射到地表后會有δPsx,δPrx,δT,δSx,δRx5種數(shù)據(jù)殘差。這種方法的特點是,運動學(xué)偏移獲得初始反射點位置后,即可確定地下反射點與地表的局部同相軸的關(guān)系,反偏移到地表得到數(shù)據(jù)域殘差δPsx,δPrx,δT,δSx,δRx,然后利用這些數(shù)據(jù)殘差更新背景速度模型。這其實就是傳統(tǒng)立體層析的實現(xiàn)思路。

圖2 聯(lián)合P參數(shù)與旅行時的各種層析反演策略分析示意a 第1種方案 (利用坐標(biāo)殘差Xerr來反演速度); b 第2種方案(用P參數(shù)成像,而后利用數(shù)據(jù)域旅行時殘差信息反演背景速度); c 第3種方案(利用數(shù)據(jù)域δPsx殘差信息反演背景速度); d 第4種方案(利用P參數(shù)與旅行時殘差信息反演背景速度)

本文的思路是將第3種方案略作調(diào)整后與第5種方案相結(jié)合。CHAURIS等[9]提出基于炮道集推導(dǎo)的P參數(shù)校正公式,同時指出,在共偏移距成像道集內(nèi)有等價的校正公式。由于Kirchhoff偏移算法輸出共偏移距成像道集是最常規(guī)選項,因此在共偏移距成像域?qū)崿F(xiàn)“運動學(xué)反偏移+P參數(shù)校正”是一個更為自然的選擇。我們在共偏移距域推導(dǎo)了相應(yīng)的運動學(xué)反偏移與校正公式。此外,在實施校正后獲得了正確的P參數(shù)之后,原則上認(rèn)為坐標(biāo)殘差和旅行時殘差不需要再進入數(shù)據(jù)空間,但考慮到實際數(shù)據(jù)中運動學(xué)反偏移的精度問題以及反演的穩(wěn)定性問題,我們依然將所有的數(shù)據(jù)殘差即δPsx,δPrx,δT,δSx,δRx都納入到數(shù)據(jù)空間用于反演。我們認(rèn)為對實際數(shù)據(jù)而言這是最為穩(wěn)妥的處理方案。

3 “運動學(xué)反偏移+P參數(shù)校正”的算法原理及其實現(xiàn)流程

本節(jié)對CHAURIS等[9]提出的“運動學(xué)反偏移+P參數(shù)校正”的推導(dǎo)過程做一詳細(xì)介紹。不同于原文基于共炮道集的公式推導(dǎo),本文基于同樣的思想,直接利用共偏移距域與共中心點域的P參數(shù)(Pmx和Phx)導(dǎo)出適用于共偏移距成像剖面與共偏移距成像道集的“運動學(xué)反偏移+P參數(shù)校正”公式,更易于為讀者所理解。圖3為運動學(xué)反偏移示意圖。

在Kirchhoff偏移中,成像條件滿足這樣的等時曲線:

(4)

式中:u是地下速度;ts,tr分別是地下成像點(x,z)到炮點、檢波點的單程走時;h=(s-r)/2是半偏移距;m=(s+r)/2是炮點、檢波點的中點;t*(h,m)對應(yīng)于一根射線從炮點S出發(fā)、經(jīng)過地下某點反射回到地表檢波點R的總旅行時。

圖3 運動學(xué)反偏移示意a 地下模型信息與地表數(shù)據(jù)信息的幾何關(guān)系; b 共偏移距成像剖面(在其上拾取構(gòu)造傾角ξ); c 偏移距域共成像點道集(在其上拾取剩余曲率的斜率tan φ)

固定橫向位置x和速度模型u,基于公式(4)對深度z,半偏移距h和中點m做全微分可得:

(5)

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

(6)

由此可得(具體推導(dǎo)見附錄A):

(8)

因為:

將(7a)式代入(10)式,得到:

(11)

如圖3a所示,根據(jù)慢度矢量的定義有:

(12)

式中:θs,θr是射線與垂向的夾角。

(13)

式中:β是反射角,即射線與反射界面法向的夾角;ξ是構(gòu)造傾角。

根據(jù)(12)式和(13)式,得到:

(14)

令α=2ucosβcosξ,將(14)式代入(8)式,得:

(15)

根據(jù)(9)式,得到:

(16)

將(7b)式代入(15)式和(16)式,得到:

公式(17)將時間域射線的出射斜率信息與成像域共成像點道集內(nèi)的運動學(xué)信息建立了聯(lián)系,同時也證明了立體層析中擬合出射斜率與成像域速度分析拉平道集內(nèi)同相軸等價。

從(4)式出發(fā)亦可推導(dǎo)出成像剖面傾角的公式,固定s,r,u,基于公式(4)對x,z求全微分,得:

(18)

化簡(18)式得到:

(19)

公式(17)的應(yīng)用價值在于,利用不正確的初始速度模型和不正確的構(gòu)造成像信息,依然有可能獲得正確的立體層析數(shù)據(jù)空間。畢竟用于疊前深度偏移的初始速度信息不可能完全正確,因此相應(yīng)的初始深度成像結(jié)果也不可能完美。但是初始成像結(jié)果依然有兩個優(yōu)點:①大部分繞射波得到了收斂歸位;②側(cè)面反射與多次波的能量將聚焦到不合理的位置上。此時,處理人員結(jié)合自己的地質(zhì)認(rèn)識即可從共偏移距成像剖面上識別出比較可靠的、同時在橫向和深度上分布也比較均勻的一次反射波的成像位置,然后利用結(jié)構(gòu)張量算法提取這些位置處的構(gòu)造傾角與剩余深度差信息,通過公式(7)實施運動學(xué)反偏移即可獲得可靠的、均勻分布的立體層析數(shù)據(jù)空間。這時再實施立體層析反演將會獲得更為可靠的反演結(jié)果。

需要強調(diào)的是,運動學(xué)反偏移對初始速度模型有一定要求。如果初始速度模型品質(zhì)不夠好,無論在共偏移距成像點道集內(nèi)提取RMO或在共偏移距成像剖面內(nèi)提取構(gòu)造傾角都將難以進行,那么后續(xù)的反偏移和校正也都失去了前提。不過好在生產(chǎn)中的初始深度模型一般來自于疊前時間偏移速度分析的時深轉(zhuǎn)換或者初至波層析。就本文而言,我們的初始模型來自全自動拾取的數(shù)據(jù)域立體層析,無論通過何種方式獲取初始模型其品質(zhì)都不會差,基于該模型的初始成像能夠滿足結(jié)構(gòu)張量算法的拾取要求。

圖4 運動學(xué)反偏移拾取參數(shù)流程

圖5顯示了“運動學(xué)反偏移+校正”在一個相對簡單的6層背斜模型上的測試結(jié)果,主要測試了算法的校正精度。圖5a展示了真實模型以及模擬某一炮的反射射線信息。基于該模型共正演模擬401炮,每炮401道,炮間距20m,道間距10m,記錄時間為4s。圖5b展示了基于該模型模擬獲得的多炮正演記錄中的一炮反射射線信息。將真實模型乘以0.9,并略做光滑后實施Kirchhoff疊前深度偏移(PSDM)后、6.2km處的CIG和2.0km共偏移距成像剖面分別如圖5b和圖5c所示,可以看到,由于速度偏小,出現(xiàn)明顯的欠偏移現(xiàn)象。圖5d展示了將運動學(xué)反偏移+P參數(shù)校正后得到的數(shù)據(jù)點位置與斜率信息貼合到對應(yīng)的2.0km共偏移距同相軸上后的質(zhì)量監(jiān)控顯示結(jié)果,可以看出,其貼合效果相當(dāng)不錯。當(dāng)然,僅從視覺上判斷其校正精度是不夠的,必須有相應(yīng)的反演結(jié)果才更有說服力。圖5e和圖5f分別顯示了利用未做校正的數(shù)據(jù)空間實施反演的結(jié)果和利用校正后的數(shù)據(jù)空間實施反演的結(jié)果,可以看出,反演結(jié)果差別很大,顯示了校正公式(17)的正確性和實施這種校正的必要性。

圖5 運動學(xué)反偏移+P參數(shù)校正在簡單理論數(shù)據(jù)上的精度測試a 真實模型及某一炮的反射射線信息; b 真實模型乘以0.9后實施PSDM后的CIG(6.2km處); c 真實模型乘以0.9后實施PSDM的2.0km共偏移距成像剖面; d 基于真實模型乘以0.9后實施PSDM的2.0km共偏移距成像剖面實施反偏移+P參數(shù)校正后獲得的正確反射點位置以及P參數(shù)信息; e 利用未做校正的數(shù)據(jù)空間實施反演的結(jié)果; f 利用校正后的數(shù)據(jù)空間實施反演的結(jié)果

4 結(jié)構(gòu)張量算法及其在數(shù)據(jù)域和成像域中的應(yīng)用

如前所述,立體層析數(shù)據(jù)空間的提取可以在數(shù)據(jù)域直接進行,也可以在成像域提取構(gòu)造傾角和RMO信息后再實施“運動學(xué)反偏移+校正”后換算得到,因此,有一個高效率、高精度的運動學(xué)信息提取方法至關(guān)重要。梯度平方結(jié)構(gòu)張量是圖像處理的一種高效的邊緣檢測算法[18]。如果將地震數(shù)據(jù)視為一張圖像,結(jié)構(gòu)張量算法是一個用于檢測地震數(shù)據(jù)局部同相軸傾角的極好方法。

二維梯度平方結(jié)構(gòu)張量算法原理可概括如下。對于一張二維圖像,為了檢測其邊緣,須首先計算出每一個樣點處x,y方向的梯度gx,gy。對于低信噪比圖像,為了提高方向預(yù)測的穩(wěn)定性,避免梯度矢量在圖像邊緣兩側(cè)互相抵消,引入梯度平方張量矩陣,使得同一走向但是方向相反的梯度矢量不至于相互抵消,反而可以相互增強。實際二維地震資料可以視為一個低信噪比圖像,這是梯度平方張量算法適用于地震數(shù)據(jù)同相軸檢測的重要原因。

對二維圖像中某一樣點而言,其梯度平方張量矩陣定義如下:

(20)

式中:gx,gy分別表示該樣點處的梯度水平分量與梯度垂直分量。梯度平方張量的含義是圖像中某一樣點處,對應(yīng)于單一走向的梯度方向。為提高低信噪比圖像中走向信息預(yù)測的穩(wěn)定性,在該樣點附近的鄰域內(nèi),比如一個4樣點×4樣點區(qū)域內(nèi)統(tǒng)計16個這樣的梯度平方張量,并將它們加權(quán)疊加得到一個光滑的梯度平方張量矩陣G′將有助于獲得該樣點位置處比較穩(wěn)定的走向信息。權(quán)系數(shù)可以設(shè)置為本區(qū)域內(nèi)定義的一個高斯函數(shù),或可以簡單設(shè)為1即平均權(quán)值,視具體的搜索效果而定。不妨將G′寫為:

(21)

式中:〈·〉代表光滑后的梯度平方。獲得光滑的梯度平方張量矩陣G′之后,針對半正定矩陣G′,通過求解特征方程|G′-λI|=0可求得其特征值與特征向量:

(22)

式中:λ1為最大特征值,對應(yīng)于張量在特征張量方向v1的能量;λ2為最小特征值,對應(yīng)于張量在特征張量方向v2的能量。(λ1-λ2)/λ1表示線性度,反映局部方向的一致性。對于地震數(shù)據(jù)而言,線性度的含義等價于同相性。特征向量描述了圖像局部線性結(jié)構(gòu)的方向性。特征向量v1正交于圖像的主結(jié)構(gòu)方向,特征向量v2平行于圖像的主結(jié)構(gòu)方向。圖6顯示了南海某實際單炮記錄以及應(yīng)用梯度平方結(jié)構(gòu)張量算法提取的屬性參數(shù)剖面。其中,圖6e所示的線性度剖面其實就是地震學(xué)家所說的同相性剖面,可以看到,線性度剖面(圖6e)和圖6a所示的炮數(shù)據(jù)同相軸有很好的對應(yīng)關(guān)系。注意,這個線性度的計算完全依賴于之前計算的特征值和特征向量,因此,如果圖6b,圖6c 和圖6d的計算有誤,那么線性度將和同相軸失去對應(yīng)關(guān)系,由此可以看出,梯度平方結(jié)構(gòu)張量算法穩(wěn)定,完全可以用于實際地震數(shù)據(jù)的局部同相軸斜率估計。

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

王宇翔等[6]對比了梯度平方結(jié)構(gòu)張量與常規(guī)的傾斜疊加及平面波摧毀算法的計算成本,結(jié)果表明,結(jié)構(gòu)張量算法的計算效率高于平面波摧毀近一個數(shù)量級,高于局部傾斜疊加兩個數(shù)量級。因此,該算法用于立體層析數(shù)據(jù)空間的準(zhǔn)備是一個自然的選擇。本文數(shù)值實驗中涉及的地表觀測P參數(shù)提取以及成像域中的構(gòu)造傾角與RMO信息的提取,都是應(yīng)用梯度平方結(jié)構(gòu)張量算法實現(xiàn)的。圖7a和圖7b分別顯示了結(jié)構(gòu)張量算法在共成像點道集(CIG)上搜索RMO以及在共偏移距成像剖面上搜索構(gòu)造傾角的結(jié)果,可以看出,搜索到的位置與同相軸的波峰疊合非常一致,表現(xiàn)出該算法的優(yōu)秀性能。

圖7 梯度平方結(jié)構(gòu)張量算法在疊前成像域的應(yīng)用效果a 在共成像點道集(CIG)上搜索RMO; b 在共偏移距成像剖面上搜索構(gòu)造傾角

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

二維理論數(shù)據(jù)基于圖8a所示的理論模型得到。該模型針對南海深水某工區(qū)的復(fù)雜斷陷模型設(shè)計?；谠撃Ｐ驼莨灿?01炮,每炮401道,道間距10m,炮間距20m,正演數(shù)據(jù)通過高階有限差分聲波模擬計算得到。反演之前首先利用結(jié)構(gòu)張量算法自動搜索出相干性比較好的同相軸斜率信息,再結(jié)合王宇翔等[6]提出的波峰判斷準(zhǔn)則和線性度準(zhǔn)則篩出地震記錄上的有效數(shù)據(jù)點。我們在這個理論數(shù)據(jù)上每隔200m偏移距做一次拾取,共獲得2.2×104個數(shù)據(jù)點。圖8b顯示了1km偏移距剖面內(nèi)搜索到的數(shù)據(jù)點信息,可以看到,盡管使用了比較高的線性度門檻,4.5s以下依然有不少數(shù)據(jù)點信息,這些數(shù)據(jù)點大多來自于繞射波。繞射是立體層析中的射線正演無法表達(dá)的波現(xiàn)象,數(shù)據(jù)空間中存在繞射波的運動學(xué)信息使得反演難以正常收斂。圖8c顯示了第1次迭代后的結(jié)果,可以看出,這時傾角反演結(jié)果已有異常發(fā)生。圖8d顯示了第20次迭代的結(jié)果,深部甚至出現(xiàn)了接近90°傾角的垂直構(gòu)造,顯然不正常。圖8e顯示了20次迭代的二范數(shù)目標(biāo)泛函下降曲線,可以看到,這個結(jié)果無法令人滿意。

在初次反演的基礎(chǔ)上采用運動學(xué)反偏移+校正的方法對數(shù)據(jù)空間進行提煉,這樣的實施策略有很明確的優(yōu)點。盡管初始疊前深度偏移(PSDM)成像剖面并不是最好的成像結(jié)果,但是已經(jīng)可以辨識出大部分?jǐn)帱c的位置。圖9顯示了運動學(xué)反偏移+校正后的數(shù)據(jù)點的反演結(jié)果。圖9a是基于圖8d所示模型實施PSDM后,用一部分小偏移距剖面疊加獲得的初始成像結(jié)果。盡管成像結(jié)果不是很理想,但是在人工拾取數(shù)據(jù)點的過程中,通過人工辨識依然可以避開絕大多數(shù)的對應(yīng)于繞射波的成像點,因此,第2次選擇的數(shù)據(jù)點都無一例外落在了反射界面上。對于實際應(yīng)用而言,非一次反射波的信息除了繞射波外還可能包含多次波甚至側(cè)面波,因此,這種人工拾取更有意義。

圖9b顯示了對圖9a所示的數(shù)據(jù)點實施運動學(xué)反偏移到共偏移距剖面后的質(zhì)量控制結(jié)果,可以看出,其絕對數(shù)量相比圖8b有明顯減少,事實上,反偏移+校正后僅獲得5000個左右的數(shù)據(jù)點,但是后續(xù)的反演中證明,這些數(shù)據(jù)點反而能提供更為合理的反演結(jié)果。圖9c和圖9d顯示了第1次和第20次迭代后的反演結(jié)果。與圖8d相比,圖9d中無論是速度反演還是傾角反演結(jié)果都要穩(wěn)定得多。圖9e顯示了20次迭代后的二范數(shù)目標(biāo)泛函下降曲線,可以看到,隨著迭代次數(shù)的增多,二范數(shù)目標(biāo)泛函會非常穩(wěn)定地下降,收斂結(jié)果明顯優(yōu)于圖8e。圖10是將圖9 反演中獲得的傾角信息重疊在真實模型上的質(zhì)量控制顯示,可以看到,傾角信息和淺中部的主要構(gòu)造界面疊合得非常一致。圖11顯示了基于初始反演速度和新反演速度的PSDM結(jié)果,可以看到,應(yīng)用新模型獲得的CIG相比用初始模型獲得的CIG有更好的聚焦效果,同時新模型對應(yīng)的最終成像剖面的分辨率更好,充分體現(xiàn)了反偏移后獲得的數(shù)據(jù)點更為合理。如果基于圖11所示的偏移結(jié)果再挖掘出一些深部構(gòu)造的數(shù)據(jù)點,重復(fù)上述流程,可以進一步提高對應(yīng)于深部構(gòu)造的反演精度,提高深層構(gòu)造的成像品質(zhì)。

圖8 在數(shù)據(jù)域直接搜索數(shù)據(jù)空間獲得的反演結(jié)果(基于中海油復(fù)雜斷陷模型理論數(shù)據(jù))a 真實模型; b 共偏移距剖面內(nèi)通過結(jié)構(gòu)張量自動拾取獲得的數(shù)據(jù)點信息; c 第1次迭代后結(jié)果; d 第20次迭代后結(jié)果; e 目標(biāo)泛函下降曲線

圖9 基于運動學(xué)反偏移+校正后的反演結(jié)果a 基于初始成像拾取的反射點; b 將反射點信息反偏移到時間域的新數(shù)據(jù)點信息; c 第1次迭代后結(jié)果; d 第20次迭代后結(jié)果; e 目標(biāo)泛函下降曲線

圖10 將新反演結(jié)果中傾角信息疊合到真實模型

圖11 基于初始反演速度的PSDM結(jié)果與基于新反演速度的PSDM結(jié)果對比a 基于初始反演速度PSDM得到的CIG; b 基于新反演速度PSDM得到的CIG; c 基于初始反演速度實施PSDM得到的成像剖面; d 基于新反演速度實施PSDM得到的成像剖面

6 結(jié)論

當(dāng)射線參數(shù)水平分量和射線出射坐標(biāo)都被納入到數(shù)據(jù)空間后,會導(dǎo)致層析反演有多種可能的實現(xiàn)路徑。本文深入分析了基于立體數(shù)據(jù)空間實現(xiàn)層析反演的多種可能性,設(shè)計了一種首先在數(shù)據(jù)域獲取初始數(shù)據(jù)空間并實施初步立體層析,而后在初始成像體內(nèi)搜索運動學(xué)信息、實施運動學(xué)反偏移+校正獲取更可靠的數(shù)據(jù)空間的兩步法策略。該方法應(yīng)用于斷陷模型理論數(shù)據(jù)的結(jié)果證明了這種策略的優(yōu)勢。

雖然初始模型并非完美,但是通過初始成像體上的運動學(xué)信息,利用運動學(xué)反偏移+校正可以獲得在地表觀測得到的正確射線參數(shù)水平分量信息。這些信息正是立體層析所需的數(shù)據(jù)空間。運動學(xué)反偏移+校正獲得數(shù)據(jù)空間的優(yōu)勢在于:通過在初始成像剖面上人工選擇數(shù)據(jù)點,可以避開原始數(shù)據(jù)上可能存在的繞射波、多次波、側(cè)面波等各種非一次反射波的干擾。需要指出,基于疊前數(shù)據(jù)體的高效結(jié)構(gòu)張量自動拾取與利用反偏移+校正重建數(shù)據(jù)空間并不矛盾,這恰恰體現(xiàn)了兩步法的優(yōu)勢:先通過高效率的結(jié)構(gòu)張量算法實施自動拾取獲得初始數(shù)據(jù)空間,基于初始數(shù)據(jù)空間反演得到一個較為理想的初始速度模型;然后在此基礎(chǔ)上通過運動學(xué)反偏移去除與繞射波、多次波、側(cè)面反射波有關(guān)的數(shù)據(jù)點,提煉出更合理的數(shù)據(jù)空間,使得立體層析反演精度得到進一步提高。

在兩步法策略中,結(jié)構(gòu)張量始終扮演著重要的角色。無論是第1步的數(shù)據(jù)域初始搜索還是第2步的成像域搜索都是基于高效率的結(jié)構(gòu)張量算法實現(xiàn),尤其是第2步在初始最小偏移距成像數(shù)據(jù)體內(nèi)人工拾取出比較可靠的數(shù)據(jù)點位置后,這些數(shù)據(jù)點的位置對應(yīng)于一次反射波的關(guān)鍵反射層,這些位置需要沿著成像點道集內(nèi)的同相軸延拓到大偏移距剖面去,這時結(jié)構(gòu)張量事先在成像點道集內(nèi)搜索到的同相性非常重要。在數(shù)據(jù)點位置的延拓完成后,依然需要使用結(jié)構(gòu)張量搜索出構(gòu)造傾角,否則運動學(xué)反偏移將無從談起。本文在典型理論數(shù)據(jù)的嚴(yán)格測試證實了結(jié)構(gòu)張量的穩(wěn)定和高效,為后續(xù)實際應(yīng)用奠定了基礎(chǔ)。

[1] BILLETTE F,LAMBARé G.Velocity macro-model estimation by stereotomography[J].Geophysical Journal International,1998,135(2):671-680

[2] LAMBARé G.Stereotomgraphy[J].Geophysics,2008,73(5):VE25-VE34

[3] BILLETTE F,LE BéGAT S,PODVIN P,et al.Practical aspects and applications of 2D stereotomography[J].Geophysics,2003,68(3):1008-1021

[4] 倪瑤,楊鍇,陳寶書.立體層析反演方法理論分析與應(yīng)用測試[J].石油物探,2013,52(2):121-130 NI Y,YANG K,CHEN B S.Stereotomography inversion method:theory and application testing[J].Geophysical Prospecting for Petroleum,2013,52(2):121-130

[5] 李振偉,楊鍇,倪瑤,等.基于立體層析反演的偏移速度建模應(yīng)用研究[J].石油物探,2014,53(4):444-452 LI Z W,YANG K,NI Y,et al.Migration velocity analysis with stereo-tomography inversion[J].Geophysical Prospecting for Petroleum,2014,53(4):444-452

[6] 王宇翔,楊鍇,汪小將,等.基于梯度平方結(jié)構(gòu)張量算法的高密度二維立體層析反演[J].地球物理學(xué)報,2016,59(1):263-276 WANG Y X,YANG K,WANG X J,et al.A high-density stereo-tomography method based on the gradient square structure tensors algorithm[J].Chinese Journal of Geophysics,2016,59(1):263-276

[7] 楊鍇,邢逢源,李振偉,等.基于非降階漢密爾頓算子的三維立體層析反演[J].地球物理學(xué)報,2016,59(9):3366-3378 YANG K,XING F Y,LI Z W,et al.3D stereo-tomography based on the non-reduced Hamiltonian[J].Chinese Journal of Geophysics,2016,59(9):3366-3378

[8] XING F X,YANG K,XUE D,et al.3D stereotomography applied to the deep-sea data acquired in the South China Sea,part II:the real case study[J].Expanded Abstracts of 86thAnnual Internat SEG Mtg,2016:3788-3793

[9] CHAURIS H,NOBLE M S,LAMBARé G,et al.Migration velocity analysis from locally coherent events in 2-D laterally heterogeneous media,Part I:theoretical aspects[J].Geophysics,2002, 67(4):1202-1212

[10] GUILLAUME P,LAMBARé G,LEBLANC O,et al.Kinematic invariants:an efficient and flexible approach for velocity model building[J].Expanded Abstracts of 78thAnnual Internat SEG Mtg,2008:3687-3692

[11] 王華忠,馮波,劉少勇,等.疊前地震數(shù)據(jù)特征波場分解、偏移成像與層析反演[J].地球物理學(xué)報,2015,58(6):2024-2034 WANG H Z,FENG B,LIU S Y,et al.Characteristic wavefield decomposition,imaging and inversion with prestack seismic data[J].Chinese Journal of Geophysics,2015,58(6):2024-2034

[12] 劉少勇,蔡杰雄,王華忠,等.山前帶地震數(shù)據(jù)射線 (束) 疊前成像方法研究與應(yīng)用[J].石油物探,2012,51(6):598-605 LIU S Y,CAI J X,WANG H Z,et al.Technology and application of beam ray prestack imaging method in foothill areas[J].Geophysical Prospecting for Petroleum,2012,51(6):598-605

[13] TEN KROODE A P E,SMIT D J,VERDEL A R.A micro local analysis of migration[J].Wave Motion,1998,28(2):149-172

[14] XU S,CHAURIS H,LAMBARé G,et al.Common angle migration:a strategy for imaging complex media[J].Geophysics,2001,66(6):1877-1894

[15] FARRA V,MADARIAGA R.Seismic waveform modeling in heterogeneous media by ray perturbation theory[J].Journal of Geophysical Research Solid Earth,1987,92(B3):2697-2712

[16] SWORD C.Tomographic determination of interval velocities from picked reflection seismic data[J].Expanded Abstracts of 56thAnnual Internat SEG Mtg,1986:657-660

[17] 熊凱.特征波框架下的背景速度與角度反射系數(shù)反演[D].上海:同濟大學(xué),2017 XIONG K.Background velocity and angle reflectivity inversion based on the characteristic wave theory[D].Shanghai:Tongji University,2017

[18] VAN VlIET L J,Verbeek P W.Estimators for orientation and anisotropy in digitized images[R].Tubingen:Proceeding of First Conference of the Advanced School for Computing & Imaging,1995:442-450

附錄A 共偏移距、共中心點剖面內(nèi)Pmx,Phx參數(shù)與共炮點、共檢波點道集內(nèi)Psx,Prx參數(shù)之間的定量關(guān)系

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

(A1)

基于公式(A1)對炮點、檢波點求導(dǎo),得:

(A2)

化簡得到:

(A3)

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

(A4)

即:

(A5)

得到:

(A6)

(編輯：顧石慶)

Inversionstrategyanddataspaceconstructionforstereo-tomographyusingstructuretensorandkinematicdemigration.Ⅰ:theory

YANG Kai1,XIONG Kai1,WANG Yuxiang2,WANG Xiaojiang3

(1.TongjiUniversity,Shanghai200092,China;2.InternationalResearchandDevelopmentCenter,AramcoFarEastBusinessServicesCompanyLimited,Beijing100013,China;3.CNOOCResearchCenter,Beijing100027,China)

Using a structure tensor and kinematic demigration,we propose a two-step scheme for data space construction and inversion in stereo-tomography.First,the scheme searches a prestack data volume to obtain the initial data space and obtains the initial inversion model and the initial prestack depth migrated (PSDM) image.Next,the scheme selects reliable data point locations in the initial PSDM image manually to search for the corresponding residual moveout (RMO) and structural dip.Then,it inverts the kinematic information from the imaging domain to the data space by kinematic demigration followed by correction,to obtain a reliable and uniform data space for stereo-tomography.Note that both the first step search in the data domain and the second step search in the imaging domain are performed using an efficient structure tensor algorithm with high precision.In addition,the paper discusses the types of possible approaches and strategies to accomplish the tomography when the horizontal component information of the ray parameter is introduced into the data space of the tomographic inversion.Finally,the reliability of the proposed scheme was verified through a synthetic data test,thus laying the foundation for its application in field data.

stereo-tomography,structure tensor,kinematic demigration,data space construction

P631

：A

1000-1441(2017)05-0694-13DOI:10.3969/j.issn.1000-1441.2017.05.010

楊鍇,熊凱,王宇翔,等.聯(lián)合結(jié)構(gòu)張量與運動學(xué)反偏移的立體層析數(shù)據(jù)空間提取與反演策略研究Ⅰ:理論[J].石油物探,2017,56(5):706

YANG Kai,XIONG Kai,WANG Yuxiang,et al.Inversion strategy and data space construction for stereotomography using structure tensor and kinematic demigration.Ⅰ:theory

[J].Geophysical Prospecting for Petroleum,2017,56(5):706

2016-11-15;改回日期：2017-06-28。

楊鍇(1972—),男,教授,博士生導(dǎo)師,主要研究方向為立體層析反演、角度域最小平方偏移與全波形反演。

國家自然科學(xué)基金面上項目(41574098,41630964)、國家科技重大專項子課題(2016ZX05026-001-03)共同資助。

This research is financially supported by the National Natural Science Foundation of China (Grant Nos.41574098,41630964) and the National Science and Technology Major Project of China (Grant No.2016ZX05026-001-03).