散亂數(shù)據(jù)插值方法及其在背景速度建模中的應(yīng)用

王咸彬,吳成梁

(1.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;2.波現(xiàn)象與反演成像研究組WPI,同濟大學(xué)海洋與地球科學(xué)學(xué)院,上海200092)

散亂數(shù)據(jù)插值方法及其在背景速度建模中的應(yīng)用

王咸彬1,吳成梁2

(1.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103;2.波現(xiàn)象與反演成像研究組WPI,同濟大學(xué)海洋與地球科學(xué)學(xué)院,上海200092)

地震全波形反演(FWI)主要利用低頻長偏移距透射波估計背景速度,有效的初始背景速度場對于透射波FWI和反射波FWI的有效應(yīng)用都十分重要。因此,根據(jù)散亂的速度和界面深度信息,對不同來源、不同精度的散亂速度場插值產(chǎn)生背景速度模型是速度建模的核心環(huán)節(jié)之一。在總結(jié)了三角剖分法、自然鄰域法、反距離加權(quán)插值法和徑向基函數(shù)插值法以及克里金插值法等方法特點的基礎(chǔ)上,采用三維克里金插值方法將實際測井?dāng)?shù)據(jù)插值成三維速度場，然后利用克里金方差融合策略和多尺度空間波數(shù)域Gabor變換兩種方法對分別來自測井?dāng)?shù)據(jù)和來自偏移速度分析的速度場進行融合,融合后的速度場的精度得到了提高,成像質(zhì)量明顯改善。最后,對基于散亂數(shù)據(jù)支撐的特征速度場的描述問題進行了討論,提出了多尺度特征速度場的表達策略。

散亂數(shù)據(jù)插值;克里金插值;速度融合;速度場特征表達;克里金方差;多尺度Gabor變換

隨著計算機技術(shù)的快速發(fā)展和“兩寬一高”采集技術(shù)的進步,地震全波形反演(FWI)技術(shù)研究也得到快速發(fā)展,有效促進了油氣勘探的發(fā)展。在以估計“全/寬”波數(shù)帶的彈性參數(shù)為目標(biāo)的FWI中,利用低頻長偏移距透射波估計背景速度占有重要位置。但是,在目前的陸上地震數(shù)據(jù)中,缺乏有效的低頻成分。有效的初始背景速度場對于透射波FWI和反射波FWI克服周期跳(Cycle-skipping)現(xiàn)象、提高收斂效率非常重要。充分利用地震數(shù)據(jù)包含的信息,根據(jù)散亂的速度和界面深度信息特點,對不同來源、不同精度的散亂速度場插值產(chǎn)生有效背景速度模型是速度建模的核心環(huán)節(jié)之一。

勘探地震某種程度上屬于信息處理的學(xué)科,抽象地看,地震數(shù)據(jù)處理中經(jīng)常碰到的各種類型的數(shù)據(jù)基本上都是散亂采樣(規(guī)則采樣可以視為散亂數(shù)據(jù)采樣的特殊形式)。對于來自不規(guī)則空間位置的測井?dāng)?shù)據(jù),有必要利用基于散亂數(shù)據(jù)插值的方法獲得全空間的、規(guī)則的參數(shù)場數(shù)據(jù)體。尤其對速度反演和建模而言,參數(shù)場(主要是速度場)的描述也希望用盡可能少的特征點支撐起整個參數(shù)場。利用非規(guī)則(散亂)的數(shù)據(jù)表達連續(xù)的函數(shù)或其一階和二階導(dǎo)數(shù)有著重要的應(yīng)用。剖分方式和基函數(shù)的選擇是其中最關(guān)鍵的兩點。Haar條件保證了一維多項式基函數(shù)插值有唯一的解,也發(fā)展了一套有效的插值方法,但高維情況下解不唯一,高維離散數(shù)據(jù)插值,尤其是高維散亂數(shù)據(jù)插值,是一個正在發(fā)展中的問題。

在不同的應(yīng)用領(lǐng)域,根據(jù)實際需要,不同的學(xué)者提出了不同的插值方法。氣象學(xué)及地質(zhì)學(xué)家SHEPARD[1]在研究降水量和高程時提出了Shepard插值法;SIBSON[2]在研究自然鄰近法的基礎(chǔ)上發(fā)展了自然鄰域插值法;HARDY[3]在處理飛機外形的曲面擬合問題時提出了以多二次曲面(Multiquadric)函數(shù)為代表的徑向基函數(shù)插值法;DUCHON[4]從樣條彎曲能量最小出發(fā)推導(dǎo)出多元問題的薄板樣條函數(shù);南非礦業(yè)工程師KRIGE[5]在礦產(chǎn)估計過程中首先提出了克里金插值技術(shù);后來MATHERON[6]提出地質(zhì)統(tǒng)計學(xué)概念,將克里金插值方法引入地質(zhì)統(tǒng)計學(xué),使得克里金插值方法得到迅速發(fā)展;SCHUMAKER[7]在1976年提出了散亂數(shù)據(jù)的兩步擬合方法。兩步擬合法是指將一個復(fù)雜的問題分解為兩個簡單的問題來處理,成為了目前應(yīng)用較多的散亂數(shù)據(jù)擬合方法。

根據(jù)實際的地下介質(zhì)特點,描述速度場的方法一般用界面加層速度的方式,界面的描述基于解釋出的散亂點,模型的描述需要對散亂數(shù)據(jù)連續(xù)化。更關(guān)鍵的是解波場方程的導(dǎo)函數(shù)也需要離散數(shù)據(jù)(散亂數(shù)據(jù))的逼近。非規(guī)則數(shù)據(jù)的插值和逼近表達連續(xù)函數(shù)以及任意點的導(dǎo)函數(shù)在地震勘探中成為一個廣為關(guān)注的議題。高維疊前數(shù)據(jù)規(guī)則化方法已經(jīng)有很多的文獻進行了討論[8-10],速度場的特征描述也逐漸得到研究。HALE[11]對自然鄰域法進行修改,發(fā)展了一種混合鄰域法,并用成像剖面作為約束,提出了一種成像指導(dǎo)的混合鄰域插值法;DOUMA等[12]將混合鄰域插值法應(yīng)用到反演中的低頻背景速度建模。我們認(rèn)為,速度場的特征描述與基于此特征描述的波傳播方法研究可以為多尺度參數(shù)反演提供良好的基礎(chǔ)。

本文首先對以下幾種散亂數(shù)據(jù)插值方法進行了分析和總結(jié),主要包括三角剖分法、自然鄰域法、反距離插值法、徑向基函數(shù)插值法以及克里金插值法,分析了上述各種方法的優(yōu)缺點。在此基礎(chǔ)上,將三維克里金插值方法應(yīng)用到不同信息來源的速度場融合中,首先采用克里金插值方法將實際測井?dāng)?shù)據(jù)插值成三維速度場,然后分別采用克里金方差融合策略和多尺度空間波數(shù)域Gabor變換2種方法融合速度場。融合后的速度場精度提高了,偏移成像質(zhì)量也得到明顯改善。最后,討論了基于散亂數(shù)據(jù)支撐的特征速度場的描述問題,提出了多尺度特征速度場的表達策略。

1 散亂數(shù)據(jù)插值方法

根據(jù)Haar插值條件,高維離散數(shù)據(jù)插值結(jié)果是不唯一的。散亂數(shù)據(jù)插值的核心包括兩點：①網(wǎng)格體系;②基函數(shù)空間。當(dāng)然,散亂數(shù)據(jù)插值也存在無網(wǎng)格、全局的插值。但是,一般而言局部的、有網(wǎng)格的插值更符合函數(shù)變化的本質(zhì),即相鄰點的函數(shù)值之間有更強的關(guān)系。

1.1 散亂數(shù)據(jù)

散亂數(shù)據(jù)插值的定義一般描述為：假設(shè)在d維歐式空間中,存在n個散亂、無規(guī)則的數(shù)據(jù)點X={x1,x2,…,xn}?Rd,每個數(shù)據(jù)點的值為：F={f1,f2…,fn}?R。散亂數(shù)據(jù)插值就是要構(gòu)建一個函數(shù)f(x):Rd→R,使得在每個已知的數(shù)據(jù)點上滿足:

f(xi)=fii=1,2,…,n

(1)

1.2 常見散亂數(shù)據(jù)插值方法

根據(jù)基函數(shù)支撐的空間范圍可將基函數(shù)分為局部支撐和全局支撐,其中局部支撐是指基函數(shù)只受周圍幾個點的影響;而全局支撐是指基函數(shù)受整個空間數(shù)據(jù)的影響?；瘮?shù)為局部支撐的插值方法叫做局部插值法,常見的局部插值方法有:三角剖分法、自然鄰域法等;基函數(shù)為全局支撐的插值方法叫做全局插值法,常見的全局插值方法有:反距離加權(quán)法、徑向基函數(shù)法以及克里金法等。

1.2.1 三角剖分插值法

三角剖分插值法是一種局部插值方法,該方法首先將整個區(qū)域分解成不同的三角形區(qū)域,在每個小區(qū)域上進行插值,然后將不同的區(qū)域拼接在一起,得到整個區(qū)域的插值結(jié)果。三角剖分插值過程分為兩步:①對區(qū)域進行三角剖分;②構(gòu)造插值基函數(shù),求取插值系數(shù)。

在步驟①中,當(dāng)數(shù)據(jù)點大于等于4時,會出現(xiàn)幾種不同的剖分形式,在實際應(yīng)用中為了滿足特定的需要,不同學(xué)者提出了不同的三角剖分判別準(zhǔn)則[13]。其中一種比較常用的判別準(zhǔn)則為最大外接圓最小的優(yōu)化剖分,又稱德勞內(nèi)(Delaunay)三角剖分[14]。Delaunay三角剖分最早可追溯到1907年提出的沃羅諾伊(Voronoi)圖[15],前人研究認(rèn)為任何Delaunay三角剖分都是Voronoi圖的對偶剖分。以20個數(shù)據(jù)點為例,它們形成的Voronoi圖和Delaunay三角剖分見圖1[16]。三角剖分是基于二維平面上的概念,當(dāng)數(shù)據(jù)擴展到三維時,三角剖分就變成了四面體剖分。四面體的剖分思想和三角剖分一致,只是在剖分的實現(xiàn)上比三角剖分復(fù)雜。

圖1 由20個點組成的Voronoi圖(a)和Delaunay三角剖分(b)

在步驟②中,可以直接選用三角剖分的3個頂點進行插值,此時構(gòu)造的插值函數(shù)可表示成:

(2)

式中:i=1,2,3是三角形的頂點;fi是頂點的值;φi(x)是構(gòu)造的插值基函數(shù),它決定了輸入數(shù)據(jù)fi的加權(quán)大小。但是此時的插值只能達到C0連續(xù)(函數(shù)連續(xù)),在三角形的頂點和邊上不連續(xù)。為了保證C1連續(xù)(一階導(dǎo)數(shù)連續(xù)),需要利用其它信息構(gòu)造更高階的插值方法。三次九參數(shù)法是改進插值連續(xù)性的一種重要方法,該方法除了利用頂點上的值,還利用了在三角形頂點的2個沿三角形邊的方向?qū)?shù)值來保證連續(xù)性,但是此時存在一個自由的參數(shù),在求解方程組時,需添加C1的連續(xù)性條件來保證三角形拼接的連續(xù)性,但是使C1連續(xù)性滿足的條件非常難以達到,有時會使得方程組無解,此時通常采用最小二乘法獲得一個與C1連續(xù)的類似曲面當(dāng)做插值曲面[17]。Clough-Tocher法[18]給出了改進插值連續(xù)性的策略,該策略在三角形片的邊界上再增加一個穿越邊界方向的方向?qū)?shù),從而滿足C1的連續(xù)性條件。在Clough-Tocher法的基礎(chǔ)上,Powell-Sabin方法[19]不僅構(gòu)造出C1連續(xù)的曲面,而且降低了分片多項式的次數(shù)。

三角剖分插值是一種局部插值方法,不需要求解大型的線性方程組,即使對大量的數(shù)據(jù)點,插值效率也非常高;當(dāng)添加新的數(shù)據(jù)點時,不需全部重新計算插值系數(shù),應(yīng)用比較靈活。對于存在速度間斷面、具有簡單構(gòu)造的地下速度場,可以將速度場劃分成一系列不規(guī)則的塊體,在塊體內(nèi)部,速度變化光滑緩慢,在塊體邊界上,速度變化劇烈。此時可以依照塊體的形狀進行三角剖分,剖分簡單,可得到比較好的插值效果,但是在邊界處可能會出現(xiàn)不連續(xù)現(xiàn)象。

三角剖分插值要求必須在插值之前完成剖分。對于簡單的區(qū)域,剖分比較方便,但是當(dāng)區(qū)域構(gòu)造比較復(fù)雜時,剖分就變得比較困難,尤其是高維情況。三角剖分插值還要求散亂數(shù)據(jù)點為凸集,且對散亂數(shù)據(jù)點的分布有一定的要求,當(dāng)散亂數(shù)據(jù)點共線(或接近共線)時,會形成狹長狀的三角形(如圖1中的邊界部分),使得插值結(jié)果變差。此外,三角剖分形式不唯一,不同剖分可能會產(chǎn)生截然不同的結(jié)果,無法確保最優(yōu)的剖分。

1.2.2 自然鄰域插值法

自然鄰域插值法又叫Sibson插值法,由SIBSON[2]在1981年提出。Sibson插值法是一種局部插值方法,是最近鄰點插值法的一種改進。最近鄰點插值法,又稱泰森多邊形插值法,其做法是首先將插值區(qū)域按照Voronoi圖進行剖分,生成鄰近區(qū)域,然后將鄰近區(qū)域內(nèi)的值等于區(qū)域內(nèi)部插值數(shù)據(jù)點的值。所以最近鄰點插值結(jié)果不連續(xù),是分段常量的。SIBSON[2]修改了最近鄰點插值法,對鄰域區(qū)域賦予一個加權(quán)系數(shù),插值時使用鄰點的加權(quán)平均來確定插值點的值。Sibson插值法可歸納為兩步:①生成Voronoi圖;②構(gòu)造插值函數(shù)。

在步驟①中,首先將所有散亂數(shù)據(jù)點按照Voronoi圖進行剖分,如圖2a;然后將待插值點x插入剖分網(wǎng)格中,形成新的Voronoi圖,如圖2b。新加入的插值點x形成的區(qū)域(圖2b中灰色部分)與已知的數(shù)據(jù)點形成的區(qū)域中會存在一條公共的邊(點),我們把與x形成的區(qū)域有公共邊(點)的已知散亂數(shù)據(jù)點稱為自然鄰域點。一般來說,在一個區(qū)域中自然鄰域點最少有d+1個(d是空間的維度),最多有n-1個(n是已知點的個數(shù))。對于二維情況,自然鄰域點的個數(shù)不超過6個。

圖2 散亂數(shù)據(jù)點形成的Voronoi圖a 所有已知點; b 插值點x加入后

在步驟②中,構(gòu)造的插值函數(shù)寫成:

(3)

式中:N代表自然鄰域點個數(shù);fi是已知散亂數(shù)據(jù)點值;φi(x)是加權(quán)系數(shù),可由下式確定:

(4)

式中:αi(x)是在x區(qū)域中與自然鄰域點相鄰的面積,如圖2b中αi所示。

我們對Sibson插值方法中的加權(quán)系數(shù)進行修改。在Sibson插值中,公式(4)中的αi(x)由剖分形成的自然鄰域面積決定,采用一種新的加權(quán)系數(shù)定義方式,可以將αi(x)修改成:

(5)

式中:si(x)是插值點x與自然鄰域點形成的區(qū)域邊長;hi(x)是插值點x到自然鄰域點距離的一半,如圖2b 所示。

SIBSON[2]和FARIN等[20]指出:除了在已知點上,自然鄰域插值法在整個區(qū)域上連續(xù)。為了保持連續(xù)性,SIBSON在插值多項式中添加了導(dǎo)數(shù)約束項使得插值函數(shù)在每個點上也連續(xù),但此方法在實際數(shù)據(jù)應(yīng)用中沒有得到實用化發(fā)展,主要原因是導(dǎo)數(shù)一般未知。WATSON[21]闡述了幾個不同的估計梯度的計算方法。當(dāng)插值點在點集外時,自然鄰域法無法計算,此時可以添加一些鬼點來構(gòu)造新的凸集。

自然鄰域法是一種局部插值方法。一般來說,與插值點x形成的區(qū)域有接觸的自然鄰域點非常少,即使當(dāng)n非常大時,自然鄰域點也遠(yuǎn)遠(yuǎn)小于n,所以構(gòu)建插值函數(shù)比較簡單;自然鄰域法加權(quán)系數(shù)取決于周圍散亂數(shù)據(jù)點與插值點形成的區(qū)域大小;當(dāng)散亂數(shù)據(jù)點的分布密度明顯不均勻時,基于區(qū)域插值的自然鄰域法可以補償數(shù)據(jù)密度的變化,效果優(yōu)于對數(shù)據(jù)敏感、基于距離插值的方法。可見,自然鄰域法是一種穩(wěn)健的、有界的插值方法,并且具有線性精度。但是,自然鄰域法需要預(yù)先進行Voronoi剖分,在三角剖分中遇到的剖分問題在自然鄰域法中同樣會碰到。

1.2.3 反距離加權(quán)插值法

反距離加權(quán)插值法又稱Shepard法,最早由SHEPARD[1]在1968年提出。其基本思想是利用散亂數(shù)據(jù)點與插值點距離不同進行插值,加權(quán)系數(shù)大小與距離成反比。反距離加權(quán)插值法可歸納為兩步:①構(gòu)造插值函數(shù);②求取加權(quán)系數(shù)。

在步驟①中,可以構(gòu)建如下插值函數(shù):

(6)

式中:wi(x)為加權(quán)系數(shù);fi為已知的散亂數(shù)據(jù)點值;n為已知點個數(shù)。加權(quán)系數(shù)可寫成:

(7)

式中:di(x)是x和xi間的歐式距離,可表示成di(x)=‖x-xi‖2;k為指數(shù)常量,一般k≥1。當(dāng)k=1,2時,插值函數(shù)變成SHEPARD[1]和GORDON等[22]給出的形式。

由公式(6)和公式(7)可看出,當(dāng)插值點x很靠近xi時,hi(x)趨于無窮大,f(xi)趨于fi,為了保證插值函數(shù)的連續(xù)性,在xi處可令f(xi)=fi,此時雖然保證了插值函數(shù)連續(xù),但是插值函數(shù)卻在散亂數(shù)據(jù)點附近的結(jié)果表現(xiàn)為平面(原因是其各階導(dǎo)數(shù)為零)。一個相應(yīng)的修改方案就是修改成帶導(dǎo)數(shù)條件的Shepard插值,把已知點的導(dǎo)數(shù)信息加進去。此時公式(6)可修改成:

(8)

式中:αi為xi處的各階導(dǎo)數(shù),但是由于在一般情況下散亂數(shù)據(jù)點的導(dǎo)數(shù)未知,所以需要利用其它方法去獲得導(dǎo)數(shù)信息。

Shepard插值是一種全局插值方法,為了適應(yīng)局部性,可以對其進行修改。其中一種簡單的修改就是限制hi(x)的影響范圍,比如可以將公式(7)中的hi(x)修改成[17]:

(9)

式中:r是選取的影響范圍半徑。[r-di(x)]+滿足條件:

(10)

這樣就可以將全局的插值法修改成局部插值,插值函數(shù)的影響也被控制在一定的范圍內(nèi)。

反距離加權(quán)插值是一種基于距離的插值方法,對數(shù)據(jù)的分布比較敏感,該方法會產(chǎn)生平滑效應(yīng),在精度和效率上比一些基于三角剖分的插值方法差一點,但是它不需要剖分,很容易擴展到高維空間,而且加權(quán)系數(shù)的選擇比較靈活,同時還可修改距離范數(shù),方便應(yīng)用于各向異性情況。

1.2.4 徑向基函數(shù)插值法

徑向基函數(shù)插值法[3-4]是一種全局插值方法,其基本思想是首先確定一個隨距離變化的基函數(shù),然后把各個散亂數(shù)據(jù)點進行線性加權(quán)求和,得到插值點的值。徑向基函數(shù)的插值過程可歸納為以下兩步:

①選擇基函數(shù),構(gòu)造插值函數(shù)。一般情況下,徑向基函數(shù)插值可寫成如下形式:

(11)

式中:n是已知散亂數(shù)據(jù)點個數(shù);λi是插值系數(shù);φ(‖x-xi‖)是隨距離變化的基函數(shù),基函數(shù)的選擇有多種形式。

②構(gòu)建插值線性方程組,求取加權(quán)系數(shù)。將所有已知散亂點的值代入到構(gòu)造的插值函數(shù)中,形成的插值線性方程組為:

(12)

式中:f為已知的散亂數(shù)據(jù)點的值形成的向量;Φ是基函數(shù)φ(‖x-xi‖)展開形式的矩陣;λ是求取的插值系數(shù)向量。求解上述線性方程組,可以得到系數(shù)向量:

(13)

插值過程中也可利用導(dǎo)數(shù)信息以獲得更好的插值效果,此時可在公式(11)的基礎(chǔ)上,添加一個多項式約束,構(gòu)造的插值函數(shù)就變成:

(14)

式中:αi為m階多項式系數(shù);m表示已知的插值曲面斜率為0的離散點個數(shù)。對公式(14)求導(dǎo)并且令導(dǎo)數(shù)等于0,可得到:

(15)

將已知的n個散亂數(shù)據(jù)點和m個導(dǎo)數(shù)值分別代入公式(14)和公式(15),可得到n+m個方程,方程中未知數(shù)個數(shù)為n+m。

為了構(gòu)造更一般地具有m階多項式精度的插值函數(shù),導(dǎo)數(shù)的約束項可換成具有m階的多項式基函數(shù),插值函數(shù)變成:

(16)

式中:q(‖x-xi‖)為m階多項式基函數(shù);λi,αi為插值系數(shù),它們滿足方程組:

(17)

對于徑向基函數(shù)插值方法而言,基函數(shù)的選擇至關(guān)重要,目前有許多種不同形式的基函數(shù),幾種常見的基函數(shù)形式如表1所示[23-24]。這些徑向基函數(shù)已經(jīng)在不同的領(lǐng)域中得到廣泛應(yīng)用[25-26]。在數(shù)學(xué)上,FRANKE等[27]做了大量的各種散亂數(shù)據(jù)插值方法的實例對比[28],結(jié)果發(fā)現(xiàn)徑向基函數(shù)插值最讓人滿意。許多學(xué)者也曾利用徑向基函數(shù)去求解各種偏微分方程問題。比如,POLLANDT[29]利用徑向基函數(shù)近似求解非線性橢圓偏微分方程中的多維邊界元方法;SONAR[30]利用薄板樣條基函數(shù)去求解雙曲型守恒律方程。

表1 常見的徑向基函數(shù)

從徑向基函數(shù)插值過程中可以看出,徑向基函數(shù)插值可歸結(jié)為一個求解線性方程組問題。但是該線性方程組的系數(shù)矩陣不是稀疏的,并且由于徑向基函數(shù)是一種全局插值方法,從而隨著數(shù)據(jù)點的增加,矩陣的規(guī)模越來越大,使得儲存量和計算量顯著增加,而且在求解過程中還可能出現(xiàn)病態(tài)問題,這都給計算帶來了巨大的挑戰(zhàn)。常規(guī)的解決方案是利用迭代法求解線性方程組,可以從目前大量的迭代法中尋找一個合適的迭代方法進行求解;另一個方法是引入緊支撐正定徑向基函數(shù)[31]。緊支撐正定徑向基函數(shù)對任何給定的連續(xù)性條件及任何給定維數(shù)的自變量空間,可以找到一個截斷多項式,使得由其產(chǎn)生的徑向函數(shù)在給定的維數(shù)空間正定且滿足給定的連續(xù)性條件,由此導(dǎo)出的線性方程組是一個稀疏矩陣。

1.2.5 克里金插值法

克里金插值法是在地質(zhì)統(tǒng)計學(xué)的基礎(chǔ)上,利用變差函數(shù)理論進行最小方差無偏估計(Best Linear Unbiased Evaluation,又稱BLUE估計)的一種全局插值方法?？死锝鸩逯捣譃閮刹?①利用區(qū)域變差函數(shù)理論統(tǒng)計區(qū)域的變差函數(shù);②利用最小方差無偏估計建立并求解關(guān)于插值系數(shù)的方程組。

變差函數(shù)又稱變異函數(shù),是地質(zhì)統(tǒng)計學(xué)的基本理論。在一維情況下,當(dāng)空間點x變化時,設(shè)區(qū)域化變量在x和x+h處的值為Z(x)和Z(x+h),則兩值之差的方差之半定義為區(qū)域化變量在x軸方向上的變差函數(shù),記為:

(18)

式中:Z(x)和Z(x+h)分別為變量在x和x+h處的值;γ(x,h)為變差函數(shù)。

(19)

二階平穩(wěn)是指隨機變量Z(x)的數(shù)學(xué)期望是一個常數(shù),并且每一對隨機變量Z(x)和Z(x+h)之間協(xié)方差函數(shù)僅依賴于兩點之間的差向量h。即有:

C(h)=C(x,x+h)=E{{Z(x+h)-

E[Z(x+h)]}{Z(x)-E[Z(x)]}}

=E[Z(x+h)Z(x)]-m2

(20)

式中:C代表協(xié)方差函數(shù)。在滿足上述假設(shè)條件下,可得到變差函數(shù)與協(xié)方差函數(shù)的關(guān)系:

(21)

變差函數(shù)存在不同的理論模型,其中常見的有線性模型、指數(shù)模型、高斯模型以及球狀模型等,它們表征了區(qū)域變量的變化形式。在計算變差函數(shù)時,無法直接計算理論的變差函數(shù),需要計算一個實驗變差函數(shù)來擬合理論變差函數(shù),從而求出理論模型的系數(shù)。由于區(qū)域變量符合二階平穩(wěn)假設(shè),可以將一對數(shù)據(jù)[Z(x),Z(x+h)]看成是一對隨機變量的實現(xiàn),所以實驗變差函數(shù)可以寫成:

(22)

式中:N(h)為數(shù)據(jù)對[Z(x),Z(x+h)]的個數(shù);γ*(h)為實驗變差函數(shù)。

得到區(qū)域變差函數(shù)后,可構(gòu)建插值函數(shù):

(23)

式中:xi為已知散亂數(shù)據(jù)點;Z(xi)為散亂數(shù)據(jù)點上的值;n為已知點個數(shù);λi為插值系數(shù);Z(x)為插值結(jié)果?？死锝鸩逯挡捎脽o偏性和最小方差來保證插值的效果,得到:

(24)

式中:Z(x)為期望理論結(jié)果;Z*(x)為實際插值得到的結(jié)果。利用拉格朗日乘數(shù)法可構(gòu)造目標(biāo)函數(shù)F:

(25)

式中:μ為拉格朗日乘數(shù)。為了滿足無偏性和最小方差,對目標(biāo)函數(shù)求導(dǎo),并使其導(dǎo)數(shù)等于零,得到:

(26)

再利用協(xié)方差與變差函數(shù)的關(guān)系,得到克里金插值方程組為:

(27)

利用計算所得到的插值系數(shù)和拉格朗日乘數(shù),可得到插值點的克里金估計方差:

(28)

對比克里金插值和徑向基函數(shù)插值過程可以看出,兩者可歸納到一個框架中,都是選擇一個基函數(shù),構(gòu)造插值函數(shù),然后形成插值線性方程組,求解插值系數(shù)。不同之處在于:在選擇基函數(shù)時,克里金插值采用的是變差函數(shù),變差函數(shù)的形式由實際數(shù)據(jù)驅(qū)動擬合得到,而徑向基函數(shù)插值中基函數(shù)的選擇可以看作是不同形式的協(xié)方差函數(shù),根據(jù)前述的協(xié)方差函數(shù)和變差函數(shù)的關(guān)系可以看出,兩者是等價的;在求取插值系數(shù)時,克里金插值額外添加了2個約束條件,本質(zhì)上還是一種線性加權(quán)方法。STEIN[32]曾證明了克里金插值過程與高斯隨機過程的等價性。

當(dāng)已知數(shù)據(jù)點比較多時,克里金方程組比較龐大,求解就變得較困難,此時需要考慮用局部的少量點進行插值。為了提高計算效率,許多學(xué)者對如何高效選點進行過研究,比如:在二維情況下,HESSAMI等[33]利用Delaunay三角剖分加快選點;在三維情況下,YAO等[34]利用Delaunay四面體剖分加快選點,并在GPU上并行運算提高了效率。

克里金插值存在諸多變種,比如:泛克里金法、協(xié)克里金法以及貝葉斯克里金法等等。雖然克里金插值用到了統(tǒng)計特性,但是克里金插值是一種確定性的地質(zhì)統(tǒng)計學(xué)方法,它只得到一個確定性的模型。后來也有學(xué)者發(fā)展出了一系列隨機地質(zhì)統(tǒng)計學(xué)方法[35-36]。

1.3 小結(jié)

散亂數(shù)據(jù)插值包括局部插值方法和全局插值方法。常見的局部插值方法有三角剖分法、自然鄰域法等;全局插值方法有反距離加權(quán)法、徑向基函數(shù)法和克里金插值法等。所有的散亂數(shù)據(jù)插值方法實現(xiàn)過程都可歸結(jié)為以下2個步驟,其中局部插值法包括:①對區(qū)域進行剖分;②構(gòu)造插值基函數(shù),求取插值系數(shù),進行插值。全局的插值法包括:①選擇或構(gòu)造合適的基函數(shù);②求取加權(quán)系數(shù),利用插值函數(shù)進行插值。

局部插值方法計算簡單,插值效率高,應(yīng)用靈活,需事先對區(qū)域進行剖分,剖分結(jié)果的好壞直接影響插值結(jié)果,插值結(jié)果連續(xù)性差,邊界性保存效果好,比較適合存在斷層、尖滅等具有速度間斷面的速度場建模,所以對于三角剖分、自然鄰域插值這類局部插值方法可以比較方便地應(yīng)用到層析成像或FWI中模型構(gòu)造約束的正則化過程中,以及用于提取構(gòu)造、巖體刻畫、地震解釋等諸多方面。另外局部插值方法可以考慮局部的特征,在精細(xì)化建模、成像等方面有非常大的應(yīng)用潛力。

全局插值方法不需要剖分,對數(shù)據(jù)的分布要求比局部插值要求低,很容易擴展到高維,插值精度高,可以根據(jù)實際數(shù)據(jù)特點靈活選擇基函數(shù),物理含義明確,而且構(gòu)建的插值函數(shù)具有平滑效應(yīng),比較適合在低波數(shù)的大尺度背景速度建模中應(yīng)用。

在全波形反演中,由于初始模型距離反演的真實的參數(shù)模型太遠(yuǎn),使得FWI無法得到預(yù)計的理論目標(biāo),甚至在實際應(yīng)用中無法收斂。有效的背景速度建模對于透射波FWI和反射波FWI克服周期跳現(xiàn)象、提高收斂效率顯得異常重要。結(jié)合局部和全局散亂數(shù)據(jù)插值方法特點和低波數(shù)的背景速度建模需求,限于篇幅原因,本文僅以克里金插值方法為例論述其在低波數(shù)背景速度建模中的應(yīng)用。

2 克里金插值在背景速度建模中的應(yīng)用

在地震速度建模過程中,經(jīng)常會遇到少量高精度的測井?dāng)?shù)據(jù),如何利用這些少量高精度數(shù)據(jù)指導(dǎo)速度建模是一個值得研究的課題。本文利用時移加窗的三維各向異性克里金插值方法插值離散(散亂)的測井?dāng)?shù)據(jù),在此基礎(chǔ)上采用2種融合策略對多源速度進行融合,獲得更高精度的速度建模結(jié)果。

本次克里金插值采用的數(shù)據(jù)是某工區(qū)實測得到的58口井的聲波時差(AC)測井?dāng)?shù)據(jù),測井?dāng)?shù)據(jù)的二維平面和三維立體展示如圖3所示。對AC測井?dāng)?shù)據(jù)做濾波、平滑等處理后進行三維克里金插值?？紤]到地下速度場的分層構(gòu)造,速度在水平方向上變化比較緩慢,而在垂直方向上變化比較劇烈,在計算區(qū)域變差函數(shù)時,在垂直方向上乘以一個各向異性系數(shù),拉伸z軸,然后分方向計算實驗變差函數(shù),采用加權(quán)最小二乘法將實驗變差函數(shù)與理論變差函數(shù)擬合。實際處理過程中,在垂向上進行加窗時移處理,既保證了所需要的精度,又提高了效率。

插值后得到一個三維速度場和一個三維速度誤差數(shù)據(jù)。同時,對工區(qū)內(nèi)采集到的地震數(shù)據(jù)做常規(guī)偏移速度分析(Migration Velocity Analysis,MVA)處理,得到一個低波數(shù)帶的三維速度場。分別從克里金插值速度場和偏移速度分析得到的速度場中隨機抽取x,y,z3個方向的二維剖面,如圖4,圖5,圖6所示。在圖4所示的方向上,從克里金速度誤差數(shù)據(jù)中抽取一個二維剖面,如圖7所示。從圖4至圖6中可看出,克里金插值速度場的大體趨勢和偏移速度分析得到的速度一致,此外,克里金插值結(jié)果中包含很多偏移速度分析中沒有的高波數(shù)成分,這是因為測井?dāng)?shù)據(jù)本身就包含精度較高的高波數(shù)成分。

分析圖3和圖7,可以看出有測井?dāng)?shù)據(jù)的地方方差比較小,而沒有測井?dāng)?shù)據(jù)的地方方差比較大。方差大表示該區(qū)域插值結(jié)果可信度低。從圖4中也可看出,方差大的地方,插值結(jié)果偏離速度分析結(jié)果也比較大?？梢娍死锝鸱讲羁梢宰鳛樵u定插值結(jié)果的一個準(zhǔn)則,它表征了插值結(jié)果偏離理論結(jié)果的程度。

圖3 實際測井?dāng)?shù)據(jù)分布展示a 平面分布; b 立體分布

圖4 沿x軸抽取三維克里金插值結(jié)果顯示a 抽取剖面位置; b偏移速度分析結(jié)果; c克里金插值結(jié)果

圖5 沿y軸抽取三維克里金插值結(jié)果顯示a 抽取剖面位置; b偏移速度分析結(jié)果; c克里金插值結(jié)果

圖6 沿z軸抽取三維克里金插值結(jié)果顯示a 抽取剖面位置; b偏移速度分析結(jié)果; c克里金插值結(jié)果

圖7 圖4a所示位置抽取的二維剖面a 克里金插值結(jié)果; b克里金方差結(jié)果

由于地下0～1000m段基本上沒有測井?dāng)?shù)據(jù),所以克里金插值結(jié)果淺層數(shù)據(jù)不可靠。為了進一步提高速度建模的質(zhì)量,我們提出了2種速度場融合方法。

1) 基于克里金方差融合速度場。

根據(jù)上文可知,克里金方差可以作為評定插值結(jié)果的一個準(zhǔn)則,克里金方差小的區(qū)域,速度精度高;克里金方差大的區(qū)域,速度精度低。我們利用偏移速度分析得到的速度來改善克里金方差大的區(qū)域的速度,提高速度場的精度。克里金方差融合可寫成:

(29)

式中:v(x,y,z)為融合后速度場;vok(x,y,z)為克里金插值速度;vmva(x,y,z)為偏移速度分析得到的速度;α(x,y,z,σ)和β(x,y,z,σ)分別為關(guān)于位置和方差的融合系數(shù)。

分別從克里金插值、偏移速度分析和基于克里金方差融合的速度場中抽取一剖面(剖面選取的方向和位置如圖4a所示),結(jié)果如圖8所示。然后從圖8中的紅線位置抽取一條線,結(jié)果如圖9所示。對圖8所示的速度場進行二維克?；舴蚱?得到的疊加剖面如圖10所示。

2) 空間波數(shù)域多尺度數(shù)據(jù)融合。

偏移速度分析中含有高精度的低波數(shù)成分,測井速度中含有中、高波數(shù)成分。一般來說,在某些未測量深度范圍得不到可靠的速度。我們將偏移速度分析和克里金插值兩個速度場進行多尺度Gabor變換,轉(zhuǎn)換到空間波數(shù)域,對不同的空間和波數(shù)進行融合,提高速度場精度。頻率域和時間域多尺度Gabor函數(shù)公式[37]表示為:

圖8 從圖4a所示位置抽取的速度剖面a 克里金插值結(jié)果; b 偏移速度分析結(jié)果; c 采用克里金方差融合結(jié)果

圖9 從圖8紅線位置抽取的速度場對比

圖10 對圖8的速度場采用克?；舴蚱频玫降寞B加剖面a 偏移速度分析結(jié)果; b 克里金插值結(jié)果; c 采用克里金方差融合結(jié)果

(30)

v(x,y,z,kz)=α(x,y,z)vok(x,y,z,kz)+

β(x,y,z)vmva(x,y,z,kz)

(31)

式中:v(x,y,z,kz)為波數(shù)域融合速度;vok(x,y,z,kz)為波數(shù)域克里金插值速度;vmva(x,y,z,kz)為波數(shù)域偏移速度分析速度;α(x,y,z)和β(x,y,z)為融合系數(shù)。然后利用Gabor反變換將波數(shù)域融合的速度變換到空間域。

分別從克里金插值、偏移速度分析和基于空間波數(shù)域多尺度數(shù)據(jù)融合的速度場中抽取一剖面(剖面選取的方向和位置如圖4a所示),結(jié)果如圖11所示。然后從圖11中的紅線位置抽取一條線,結(jié)果如圖12所示。分別對圖11所示的速度場進行二維克希霍夫偏移,得到的疊加剖面如圖13所示。

圖11 從圖4a所示位置抽取的速度剖面a 克里金插值結(jié)果; b 偏移速度分析結(jié)果; c 采用空間波數(shù)域多尺度數(shù)據(jù)融合的結(jié)果

圖12 從圖11紅線位置抽取的速度場對比

圖13 對圖11速度場采用克?；舴蚱频玫降寞B加剖面a 偏移速度分析結(jié)果; b 克里金插值結(jié)果; c 采用空間波數(shù)域多尺度數(shù)據(jù)融合的結(jié)果

上述處理表明,克里金插值結(jié)果和偏移速度分析結(jié)果可融合在一起得到更高精度的速度場,并且偏移成像結(jié)果得到明顯改善。本文提出的速度融合方法為提高速度建模精度提供了一條實用的方法。

3 速度場的特征表達及多尺度反演

我們把地下介質(zhì)的分布特征抽象為:空間上廣泛分布的層狀沉積層,加上火山活動,構(gòu)造運動等引起的大、小尺度的速度異常體(或波阻抗變化的異常體)。

此介質(zhì)分布特征基本上可以用反射界面+層速度來描述。更進一步地,可以用空間散亂分布的特征點支撐起速度場的背景變化。這樣的表達對速度場的反演具有十分重要的作用,可以大大減少要反演的模型空間,提高速度反演的穩(wěn)定性。該思想也可以推廣到諸如密度、各向異性參數(shù)和Q值的反演中。

利用B樣條函數(shù)基于離散的解釋數(shù)據(jù)描述反射界面,然后在反射界面間充填層速度值是目前基本的速度模型描述方法[38]。圖14a為原始的Marmousi模型,共750×248個點。圖14b是利用等間隔的30道每道90個點進行B樣條插值,然后再橫向線性插值得到的結(jié)果?？梢娫摻Y(jié)果基本表達了Marmousi模型的背景速度變化。

事實上,與地震波場數(shù)據(jù)一樣,速度參數(shù)場也是有特征的,其特征就是反射界面的幾何形狀。在成像剖面上進行地震解釋(可以是自動的),得到特征散亂點的位置和速度值,利用合適的散亂數(shù)據(jù)插值方法可以用更少的特征點很好地描述速度場的背景變化。這樣的描述方法可以很方便地與Monte Carlo類的速度反演方法結(jié)合。更進一步地,若特征散亂點的剖分與波動方程的數(shù)值模擬結(jié)合在一起,可以發(fā)展一套多尺度的速度場反演方法。

圖14 速度模型描述結(jié)果a 原始的Marmousi模型(采用750×248個點描述); b 利用B樣條插值結(jié)果表達的Marmousi模型(采用30×90個點描述)

4 結(jié)論

散亂數(shù)據(jù)插值(與擬合)目前已經(jīng)在地震勘探中得到廣泛應(yīng)用,特別是在速度(或其它參數(shù))建模中的作用越來越重要。散亂數(shù)據(jù)插值包括局部(有網(wǎng)格)插值方法和全局(無網(wǎng)格)插值方法。前者的代表方法有三角剖分法、自然鄰域法以及樣條函數(shù)法等;后者有反距離加權(quán)法、徑向基函數(shù)法(克里金插值法也屬于徑向基函數(shù)方法)。前者計算簡單,需事先對區(qū)域進行剖分,插值結(jié)果連續(xù)性差,擴展到高維困難;后者不需剖分,容易擴展到高維,但有過強的平滑效應(yīng),插值系數(shù)求解計算量大。

利用三維克里金插值方法將各種信息來源的速度值融合在一起,是信息融合的主要方法。本文將實際測井?dāng)?shù)據(jù)插值成三維速度場,在實現(xiàn)三維克里金插值過程中,考慮到地下速度場在水平方向上變化比較緩慢,而在垂直方向上變化比較劇烈的情況,采用在垂向上乘以各向異性系數(shù)、分方向計算實驗變差函數(shù)方法來適應(yīng)地下速度場的變化特點。在計算過程中,采用加窗時移處理,提高了三維插值效率。然后采用克里金方差融合和頻率空間域融合方法將來自測井的速度場和來自偏移速度分析的速度場融合在一起。實驗結(jié)果表明融合后的速度場有更好的成像質(zhì)量,可以為后續(xù)的速度反演提供更好的初始模型。

散亂數(shù)據(jù)插值(或擬合)是一個反問題,如何根據(jù)信號或圖像的特征來約束插值結(jié)果是目前關(guān)注的重點。另一個關(guān)鍵的問題是在地震參數(shù)反演過程中,將要反演的速度場的特征表達(一般地,可以提為基于散亂點的表達)、基于特征表達的地震波正演模擬和地震波反演方法融合在一起,形成多尺度的速度反演方法。我們認(rèn)為這是散亂數(shù)據(jù)插值(擬合)在地震勘探中最有意義的發(fā)展方向。

[1] SHEPARD D.A two-dimensional interpolation function for irregularly-spaced data[J].Proceedings of the 23rdACM National Conference,1968:517-524

[2] SIBSON R.A brief description of natural neighbour interpolation[C]∥BARNETT V.Interpreting multivariate data.New York:John Wiley & Sons,1981:21-36

[3] HARDY R L.Multiquadric equations of topography and other irregular surfaces[J].Journal of Geophysical Research,1971,76(8):1905-1915

[4] DUCHON J.Splines minimizing rotation-invariant semi-norms in Sobolev spaces[M].Berlin:Springer,1977:85-100

[5] KRIGE D G.A statistical approach to some mine valuations and allied problem on the Witwatersrand[D].Johannesburg:University of Witwaterskand,1951

[6] MATHERON G.Principles of geostatistics[J].Economic Geology,1963,58(8):1246-1266

[7] SCHUMAKER L L.Fitting surfaces to scattered data[M].NewYork:Academic Press,1976:203-268

[8] ABMA R,KABIR N.3D interpolation of irregular data with a POCS algorithm[J].Geophysics,2006,71(6):E91-E97

[9] GAO J,STANTON A,SACCHI M D.Parallel matrix factorization algorithm and its application to 5D seismic reconstruction and denoising[J].Geophysics,2015,80(6):V173-V187

[10] ELY G,AERON S,HAO N,et al.5D seismic data completion and denoising using a novel class of tensor decompositions[J].Geophysics,2015,80(4):V83-V95

[11] HALE D.Image-guided blended neighbor interpolation of scattered data[J].Expanded Abstracts of 79thAnnual Internat SEG Mtg,2009:1127-1131

[12] DOUMA J,NAEINI E Z.Application of image-guided interpolation to build low frequency background model prior to inversion[J].Expanded Abstracts of 76thEAGE Annual Conference,2014:We G106 05

[13] 周曉云,朱心雄.散亂數(shù)據(jù)點三角剖分方法綜述[J].工程圖學(xué)學(xué)報,1993,14(1):48-54 ZHOU X Y,ZHU X X.A survey of triangulation methods for seattered data points[J].Journal of Engineering Graphics,1993,14(1):48-54

[14] DELAUNAY B.Neue darstellung der geometrischen kristallographie[J].Zeitschrift für Kristallographie-Crystalline Materials,1933,84(1):109-149

[15] VORONOI G M.Nouvelles applications des parameters continues a la theorie des formes quadratiques.deuxieme memoire:recherches sur les parallelloedres primitives[J].Journal fur die Reine und Angewandte Mathematik,1908(134):198-287

[16] ISKE A.Multiresolution methods in scattered data modelling[M].Berlin:Springer Science & Business Media,2004:1-182

[17] 吳宗敏.散亂數(shù)據(jù)擬合的模型、方法和理論[M].北京:科學(xué)出版社,2007:1-166 WU Z M.Scattered data fitting model,method and theory[M].Beijing:Science Press,2007:1-166

[18] CLOUGH R W,TOCHER J L.Finite element stiffness matrices for analysis of plates in bending[J].Proceedings of Conference on Matrix Methods in Structural Analysis,1965:515-545

[19] POWELL M J D,SABIN M A.Piecewise quadratic approximations on triangles[J].ACM Transactions on Mathematical Software (TOMS),1977,3(4):316-325

[20] FARIN G.Surfaces over Dirichlet tessellations[J].Computer Aided Geometric Design,1990,7 (1/2/3/4):281-292

[21] WATSON D F.Contouring:a guide to the analysis and display of spatial data[M].New York:Pergamon Press,1992:1-321

[22] GORDON W J,WIXOM J A.Shepard’s method of “metric interpolation” to bivariate and multivariate interpolation[J].Mathematics of Computation,1978,32(141):253-264

[23] POWELL M J D.The theory of radial basis function approximation in 1990[C]∥LIGHT W A.Advances in numerical analysis.Oxford:Clarendon Press,1992:105-203

[24] LIGHT W A.Some aspects of radial basis function approximation[C]∥SINGH S P.Approximation theory,spline functions and applications.Dordrecht:Springer Science+Business Media.1992:163-190

[25] ZALA C A,BARRODALE I.Warping aerial photographs to orthomaps using thin plate splines[J].Advances in Computational Mathematics,1999,11(2/3):211-227

[26] HARDY R L.Theory and applications of the multiquadrics-biharmonic method[J].Computers & Mathematics with Application,1990,19 (8/9):163-208

[27] FRANKE R,NIELSON G M.Scattered data interpolation and applications:a tutorial and survey[C]∥HAGEN H,ROLLER D.Geometric Modeling.Berlin:Springer Berlin Heidelberg,1991:131-160

[28] FRANKE R.Scattered data interpolation:test of some methods[J].Mathematics of Computation,1982,38(157):181-200

[29] POLLANDT R.Solving nonlinear differential equations of mechanics with the boundary element method and radial basis functions[J].International Journal for Numerical Methods in Engineering,1997,40(1):61-73

[30] SONAR T.Optimal recovery using thin plate splines in finite volume methods for the numerical solution of hyperbolic conservation laws[J].IMA Journal of Numerical Analysis,1996,16(4):549-581

[31] WENDLAND H.Piecewise polynomial,positive definite and compactly supported radial basis functions of minimal degree[J].Advances in Computational Mathematics,1995,4(1):389-396

[32] STEIN M L.Interpolation of spatial data:some theory for Kriging[M].New York:Springer Science & Business Media,1999:1-249

[33] HESSAMI M,ANCTIL F,VIAU A A.Delaunay implementation to improve Kriging computing efficiency[J].Computers & Geosciences,2001,27(2):237-240

[34] YAO X M,WANG Q,LIU Z N,et al.Fast 3D Kriging interpolation using Delaunay tetrahedron with CUDA-enabled GPU[J].Expanded Abstracts of 85thAnnual Internat SEG Mtg,2015:1739-1743

[35] DEUTSCH C V.A sequential indicator simulation program for categorical variables with point and block data:BlockSIS[J].Computers & Geosciences,2006,32(10):1669-1681

[36] DELBARI M,AFRASIAB P,LOISKANDL W.Using sequential Gaussian simulation to assess the field-scale spatial uncertainty of soil water content[J].Catena,2009,79(2):163-169

[37] QIAN J L,YING L X.Fast Gaussian wavepacket transforms and Gaussian beams for the Schr?dinger equation[J].Journal of Computational Physics,2010,229(20):7848-7873

[38] LING Y,WANG H Z,LIU S Y.Monte Carlo background velocity inversion method[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:735-738

(編輯:朱文杰)

Scattered data interpolation methods and its application in background velocity model building

WANG Xianbin1,WU Chengliang2

(1.SinopecGeophysicalResearchInstitute,Nanjing,211103,China;2.WavePhenomenaandInversionImagingResearchGroup(WPI),SchoolofOceanandEarthScience,TongjiUniversity,Shanghai200092,China)

Full waveform inversion is the current hot research topic in the field of seismic exploration,mainly using the low frequency and long offset transmitted wave to estimate the background velocity.The lack of effective low frequency component about seismic data is the dominating problem in land seismic exploration.The accuracy and validity of initial background velocity is very important for transmission FWI and reflection wave FWI.According to the scattered velocity data and interface depth information,no matter where they come from,interpolating the scattered data to produce an accurate background velocity model is the core of velocity model building.We firstly summarize some of the scattered data interpolation methods such as triangulation based interpolation,natural neighbor interpolation,inverse distance weighted,radial basis function interpolation and Kriging interpolation.Then,we employ 3D Kriging interpolation method to interpolate the practical log data into 3D velocity field and provide two methods to merge the velocity fields from log data and velocity analysis for migration.One of these methods is utilizing Kriging variance and the other one is making use of the multiscale Gabor transform.After that we apply the merged velocity fields to migration imaging process.Application results show that the accuracy of velocity field can be enhanced and both these merged methods improve the imaging quality.In the end,we discuss the description problem on the characteristics of the velocity field based on scattered data interpolation and put forward the multi-scale expression of characteristic velocity field.

scattered data interpolation,Kriging interpolation,velocity merging,characteristic expression of velocity field,Kriging variance,multiscale Gabor transform

2016-10-08;改回日期：2016-12-04。

王咸彬(1965—),男,高級工程師,博士,現(xiàn)主要從事地球物理方法技術(shù)研究與管理工作。

國家重點基礎(chǔ)研究發(fā)展計劃(973計劃)(2011CB201002)、國家自然科學(xué)基金(41374117)和國家科技重大專項(2011ZX05005-005-008HZ,2011ZX05006-002,2011ZX05023)共同資助。

P631

A

1000-1441(2017)01-0126-15

10.3969/j.issn.1000-1441.2017.01.015

This research is financially supported by the National Key Basic Research Program of China (973 Program)(Grant No.2011CB201002),The National Natural Science Fund of China (Grant No.41374117) and the National Science and Technology Major Project of China (Grant Nos.2011ZX05005-005-008HZ,2011ZX05006-002,2011ZX05023).