王 珺，楊長(zhǎng)春，李波濤

(1．石油大學(xué) 信息與控制工程學(xué)院，東營(yíng)257061;2．中科院地質(zhì)與地球物理研究所，北京100029;3．中國(guó)科學(xué)院 蘭州地質(zhì)研究所，蘭州730000)

壓縮三維走時(shí)表提高克?；舴蚱朴?jì)算效率

王 珺1，楊長(zhǎng)春2，李波濤3

(1．石油大學(xué) 信息與控制工程學(xué)院，東營(yíng)257061;2．中科院地質(zhì)與地球物理研究所，北京100029;3．中國(guó)科學(xué)院 蘭州地質(zhì)研究所，蘭州730000)

地震波走時(shí)是克?；舴虔B前深度偏移的重要參數(shù)。三維克?；舴蚱瞥上裼捎谛桀l繁讀取大數(shù)據(jù)量的走時(shí)表文件，所以計(jì)算效率不高。這里提出一種通過壓縮三維射線追蹤走時(shí)表，來提高克?；舴蚱朴?jì)算效率的方法:首先規(guī)劃一個(gè)具有規(guī)則網(wǎng)格控制點(diǎn)并覆蓋所有走時(shí)表點(diǎn)集的最小長(zhǎng)方體區(qū)域，然后以三維三次B樣條函數(shù)為插值基函數(shù)進(jìn)行最小二乘法曲面體擬合，求出規(guī)則網(wǎng)格控制點(diǎn)的數(shù)值并以數(shù)組形式存儲(chǔ)入內(nèi)存，采用稀疏化存儲(chǔ)進(jìn)一步節(jié)省了內(nèi)存空間。在偏移成像時(shí)，再由這些規(guī)則網(wǎng)格控制點(diǎn)的數(shù)值，使用線性插值公式解編出走時(shí)表。實(shí)際資料算例驗(yàn)證了該走時(shí)表壓縮方法不僅近似精度高，計(jì)算穩(wěn)定度高，計(jì)算效率高，而且由于省去了頻繁進(jìn)行大數(shù)據(jù)量走時(shí)表文件的讀寫操作，所以克?；舴蚱频挠?jì)算效率提高了二倍以上。

克?；舴蚱?計(jì)算效率;三維走時(shí)表壓縮

0 前言

用射線追蹤方法計(jì)算地震波走時(shí)是克希霍夫疊前深度偏移的基礎(chǔ)。在三維地震克希霍夫偏移成像過程中，三維射線追蹤的走時(shí)表往往具有相當(dāng)大的數(shù)據(jù)量，通常作為專門的文件來存儲(chǔ)，成像時(shí)需頻繁地讀取此走時(shí)表文件，這樣勢(shì)必造成成像計(jì)算效率不高。如果能夠?qū)⒋髷?shù)據(jù)量的三維射線追蹤的走時(shí)表壓縮存儲(chǔ)，用有限個(gè)數(shù)值點(diǎn)來代替，直接以數(shù)組的形式存儲(chǔ)在內(nèi)存中，就省去了頻繁地進(jìn)行較大數(shù)據(jù)量文件的讀寫操作，自然會(huì)提高克希霍夫偏移成像的計(jì)算效率。對(duì)射線追蹤走時(shí)表進(jìn)行的科學(xué)合理的壓縮存儲(chǔ)方法，還可用于提高正演模擬以及層析成像等的計(jì)算效率，具有重大的理論和現(xiàn)實(shí)意義。目前，國(guó)內(nèi)、外有一些關(guān)于常規(guī)地震資料的壓縮方法的研究工作(國(guó)內(nèi)有郭洪升，陳俊良［1]、張軍華，仝兆岐［2]、王喜珍等［3]、余平［4];國(guó)外有Averbuch、Meyer［5]、Bernasconi［6]、Wu［7]、Giancarlo［8]、E．J．Candès［9]、Cheng and Zhang［10])。

上述關(guān)于地震數(shù)據(jù)壓縮的研究工作可大體分為兩類:

(1)一類是有損失數(shù)據(jù)壓縮，實(shí)際上就是采用基于變換的方法，比如小波變換方法，正交變換方法等。

(2)第二類是無損失數(shù)據(jù)壓縮，如將地震資料轉(zhuǎn)換為文本形式，然后對(duì)文本進(jìn)行壓縮。其中，有損失壓縮方法壓縮比較高，能滿足信號(hào)傳輸?shù)乃俣纫蟆?/p>

但是，在壓縮前后資料的誤差較大，這對(duì)后期地震資料的精確處理不利。無損失壓縮方法雖然在壓縮前后數(shù)據(jù)具有較小的誤差，但壓縮比很低，最高才能達(dá)到2，無法滿足信號(hào)存儲(chǔ)及傳輸?shù)目臻g和速度要求。此外，上述地震資料壓縮方法的壓縮目的，均為減少地震資料數(shù)據(jù)體以節(jié)省存儲(chǔ)空間或者利于遠(yuǎn)距離傳輸。通過壓縮走時(shí)表來提高地震偏移計(jì)算效率方面的研究，尚未見到報(bào)道。

在三維克希霍夫偏移過程中，計(jì)算得到的射線追蹤走時(shí)表屬于大規(guī)模散亂數(shù)據(jù)。散亂數(shù)據(jù)指的是在二維平面上或三維空間中，無規(guī)則的、隨機(jī)分布的抽樣數(shù)據(jù)點(diǎn)。散亂數(shù)據(jù)的曲面擬合，是指通過這一系列無規(guī)則的抽樣數(shù)據(jù)點(diǎn)構(gòu)造一個(gè)光滑的曲面，以實(shí)現(xiàn)在對(duì)這些散亂數(shù)據(jù)進(jìn)行分析時(shí)，不但可以知道抽樣位置的數(shù)值，還可以知道任意位置的數(shù)值。鑒于三次B樣條函數(shù)具有連續(xù)二階導(dǎo)數(shù)的光滑性、計(jì)算的快速性，控制節(jié)點(diǎn)的穩(wěn)定性，在滿足計(jì)算精度的前提下，對(duì)于插值和曲面擬合具有最優(yōu)性。作者在本文提出了采用三次B樣條函數(shù)曲面體最小二乘擬合，對(duì)三維射線追蹤走時(shí)表進(jìn)行壓縮存儲(chǔ)，以提高克?；舴蚱瞥上裼?jì)算效率的方法。由于本算法最小二乘法線性方程組的系數(shù)矩陣為稀疏矩陣，且是非零元素規(guī)則排列(為25個(gè)條帶，帶寬均為5的帶狀結(jié)構(gòu))，所以采用稀疏矩陣存儲(chǔ)方法對(duì)矩陣進(jìn)行存儲(chǔ)節(jié)省了內(nèi)存空間。實(shí)際三維地震資料驗(yàn)證了本算法具有很好的適用性，可以保證壓縮擬合的計(jì)算效率和計(jì)算精度，提高了三維克希霍夫疊前深度偏移的計(jì)算效率。

1 Kirchhoff偏移基本原理

根據(jù)Claerbout成像原理，介質(zhì)中的成像點(diǎn)r處的偏移成像，可以由該點(diǎn)處的散射波場(chǎng)與入射波場(chǎng)的比值，在頻率域積分得到。即

其中 成像點(diǎn)r處的散射波場(chǎng)為從檢波點(diǎn)處反傳播到介質(zhì)中，使能量集中于成像點(diǎn)r處的反傳播場(chǎng)Ub(r，ω)，入射波場(chǎng)為來自震源的源場(chǎng)Usrc(r，ω)。

根據(jù)半空間壓力場(chǎng)的Kirchhoff積分公式，反傳播場(chǎng)見式(2)。

入射波場(chǎng)見式(3)。

其中r、rs、rd分別為成像點(diǎn)、源點(diǎn)和檢波點(diǎn)的坐標(biāo)(x，y，z);U(rd，ω)為頻率為 ω 時(shí)檢波點(diǎn)rd處的壓力場(chǎng);n?為檢波點(diǎn)rd所在曲面的法向量;G*(rd，r，ω)為成像點(diǎn)r到檢波點(diǎn)rd的格林函數(shù)的傅氏變換的復(fù)共軛，G*(rd，r，ω)=A(rd，r)exp(－iωT(rd，r));A(rd，r);T(rd，r)為沿著從成像點(diǎn)r到檢波點(diǎn)rd的射線的振幅和走時(shí);F(ω)為源函數(shù)的傅氏變換;G(rs，r，ω)為成像點(diǎn)r到震源rs的格林函數(shù)的傅氏變換，G(rs，r，ω)=A(rs，r)exp(iωT(rs，r))，A(rs，r)、T(rs，r)為沿著從源點(diǎn)rs到成像點(diǎn)r的射線振幅和走時(shí)。上述振幅和走時(shí)，均可以通過射線追蹤得到。由上述討論可見，成像點(diǎn)處的走時(shí)數(shù)據(jù)是計(jì)算格林函數(shù)的必要參數(shù)，因而也是克希霍夫偏移的重要參數(shù)。但三維射線追蹤的走時(shí)表往往具有相當(dāng)大的數(shù)據(jù)量，需要作為一個(gè)文件來存儲(chǔ)，在成像時(shí)需頻繁地讀取這個(gè)大數(shù)據(jù)量的走時(shí)表文件，從而導(dǎo)致偏移成像的計(jì)算效率不高。

2 三維走時(shí)表壓縮

在三維情況下，單炮射線追蹤的走時(shí)表用t=f(x，y，z)表示，其中(x，y，z)是三維空間中離散數(shù)據(jù)點(diǎn)的坐標(biāo)，t是對(duì)應(yīng)于此點(diǎn)的走時(shí)值。對(duì)三維系統(tǒng)的散亂數(shù)據(jù)進(jìn)行曲面擬合和數(shù)據(jù)壓縮，實(shí)際上是利用已知的散亂數(shù)據(jù)，描述連續(xù)曲面的過程。在曲面擬合時(shí)，首先尋找一個(gè)光滑的插值基函數(shù)，然后確定散亂數(shù)據(jù)在三維立體中坐標(biāo)(x，y，z)的范圍，規(guī)劃出一個(gè)具有規(guī)則網(wǎng)格控制點(diǎn)并包含所有散亂數(shù)據(jù)點(diǎn)集的最小長(zhǎng)方體區(qū)域，再利用最小二乘法擬合，求出規(guī)則網(wǎng)格控制點(diǎn)的數(shù)值。在存儲(chǔ)時(shí)，只需要保留這些有限個(gè)規(guī)則網(wǎng)格控制點(diǎn)的數(shù)值，就可以代表原有的散亂數(shù)據(jù)。根據(jù)偏移成像計(jì)算的需要，使用線性插值公式就可以解編出原有的散亂數(shù)據(jù)。

2．1 選取插值基函數(shù)

就現(xiàn)有的插值方法而言:

(1)高次函數(shù)插值的計(jì)算量大，有劇烈振蕩，數(shù)值穩(wěn)定性差。

(2)在分段插值中，分段線性插值在分段點(diǎn)上僅連續(xù)但不可導(dǎo)，分段三次埃爾米特插值僅能保證一階導(dǎo)數(shù)連續(xù)，因此光滑程度常不能滿足物理問題的需要。

(3)樣條插值可以同時(shí)解決這二個(gè)問題，其基本思想是:把整個(gè)區(qū)間分段，各段分別用低次多項(xiàng)式來逼近一個(gè)函數(shù)，使整個(gè)函數(shù)成“裝配式”，同時(shí)保證接縫(結(jié)點(diǎn))處具有一定程度的光滑性(直到若干階導(dǎo)數(shù)為連續(xù))［11]。這樣既可避免高次插值所具有的數(shù)值不穩(wěn)定的龍格現(xiàn)象，又可避免一般的分段插值所導(dǎo)致的插值曲線不光滑的缺點(diǎn)。

因此，樣條插值這種分段低次插值能夠以低代價(jià)而獲得較好的效果，在那些要求較高光滑性和收斂性條件的逼近問題中具有重要意義。

B樣條曲線是貝塞爾(Bézier)曲線的拓廣，與貝塞爾曲線相比，其突出優(yōu)點(diǎn)是對(duì)局部的修改，不會(huì)引起樣條形狀變化的遠(yuǎn)距離傳播，能夠更好地實(shí)現(xiàn)曲線形狀的局部控制。也就是說，在修改樣條的某些部份時(shí)，不會(huì)過多地影響曲線到其它部份［12]。

k次B樣條函數(shù)φk是分段k次多項(xiàng)式，其在(－∞，+∞)區(qū)間內(nèi)的定義為:

其中的值可由式(6)來定義。

k次B樣條函數(shù)φk的k－1階導(dǎo)數(shù)連續(xù)，其中k階導(dǎo)數(shù)的間斷點(diǎn)為:

x0、x1、…、xk+1也是k次B樣條函數(shù)φk的節(jié)點(diǎn)。φk雖然在整個(gè)實(shí)數(shù)軸(－∞，+∞)上都有定義，但在區(qū)間［x0，xk+1]，即之外，其值恒為零。

在三維情況下，B樣條函數(shù)具有自變量獨(dú)立的性質(zhì)，將x、y和z三個(gè)方向的B樣條函數(shù)相乘，就可以得到三維B樣條函數(shù)的表達(dá)式(8)。

三維B樣條函數(shù)構(gòu)成四維空間一個(gè)曲面體(見圖1)。圖1為三維三次B樣條函數(shù)分別在x=0、y=0和z=0三個(gè)平面的切片示意圖。

圖1 三維B樣條函數(shù)在x=0、y=0和z=0三個(gè)平面的切片F(xiàn)ig．1 Slices of 3D cubic B-spline function at x=0，y=0，z=0 respectively

實(shí)際上，三維B樣條函數(shù)在y方向、z方向上的切片，與相應(yīng)位置上x方向的切片在形態(tài)上是完全一致的，這是因?yàn)槿SB樣條函數(shù)本身的構(gòu)建，就是由相互獨(dú)立的三個(gè)方向的一維B樣條函數(shù)相乘而來的，具有x、y、z三個(gè)自變量獨(dú)立的性質(zhì)。而且每一個(gè)平面切面的形態(tài)，類似于二維B樣條函數(shù)。同時(shí)在x=－1平面的切片和x=1平面的切片，在形態(tài)上是完全一致的，這是因?yàn)槊恳粋€(gè)獨(dú)立的一維B樣條函數(shù)都是偶函數(shù)，這樣相乘而來的三維B樣條函數(shù)，具有分別關(guān)于x=0、y=0和z=0三個(gè)平面對(duì)稱的性質(zhì)。

在插值和曲面擬合時(shí)，使用次數(shù)越高的B樣條函數(shù)，其相對(duì)誤差越小。但這并不是說B樣條函數(shù)的次數(shù)越高，其壓縮擬合效果就越好。因?yàn)闃訔l次數(shù)越高，其計(jì)算量越大，高次冪計(jì)算的截?cái)嗾`差也越大，這在邊界上容易出現(xiàn)激烈抖動(dòng)現(xiàn)象，從而造成邊界點(diǎn)的相對(duì)誤差較大。事實(shí)上，三次B樣條函數(shù)插值和曲面擬合的效果已經(jīng)很好，能夠滿足計(jì)算精度。三次B樣條最高冪次為三次，在每一區(qū)間上振動(dòng)不太大，連接點(diǎn)處具有二級(jí)光滑程度，這保證了根據(jù)實(shí)際問題設(shè)置的端點(diǎn)控制的計(jì)算穩(wěn)定性。因此，作者選用三次B樣條函數(shù)對(duì)走時(shí)表數(shù)據(jù)進(jìn)行擬合。

2．2 曲面體擬合

如圖2(見下頁(yè))所示，設(shè)已知N個(gè)散亂數(shù)據(jù)點(diǎn)f(xk，yk，zk)在x、y、z三個(gè)方向上的網(wǎng)格節(jié)點(diǎn)數(shù)目分別為Xres、Yres、Zres，這樣就可以用三維空間中的一個(gè)最小長(zhǎng)方體區(qū)域Ω={(x，y，z)|0≤x≤Xres，0≤y≤Yres，0≤z≤Zres}(圖2中的小長(zhǎng)方體)來描述這些散亂數(shù)據(jù)點(diǎn)集在三維立體中的投影區(qū)域。

圖2 三維三次B樣條函數(shù)插值示意圖Fig．2 3D cubic B－spline function interpolation schematic

為了逼近散亂點(diǎn)集，構(gòu)造一個(gè)連續(xù)均勻的三維三次B樣條曲面體，這一曲面體由覆蓋在矩形區(qū)域Ω上的控制頂點(diǎn)網(wǎng)格D(圖2中外面的大長(zhǎng)方體)來描述。為了消除邊界因素地影響，分別向四周擴(kuò)展了二個(gè)網(wǎng)格單位的邊界。因此，D為共包含(Xres+4)*(Yres+4)*(Zres+4)個(gè)控制點(diǎn)的網(wǎng)格，且所有控制點(diǎn)均落在矩形域中的整數(shù)網(wǎng)格上。將每一個(gè)網(wǎng)格控制點(diǎn)的值表示為Cijl，其中i、j、l分別為x、y、z三個(gè)方向的序號(hào)，并且i=－1、0、…、Xres+2;j=－1、0、…、Yres+2;l=－1，0，…，Zres+2。因?yàn)樗械纳y數(shù)據(jù)點(diǎn)都落在Ω中，所以只要計(jì)算出D上的三維三次B樣條曲面體控制點(diǎn)網(wǎng)格的值，就可以用插值方法來逼近數(shù)據(jù)點(diǎn)，集中的所有散亂數(shù)據(jù)。

因?yàn)槿螛訔l函數(shù)的值不為零的定義域?yàn)?－2，2)，受其限制，每一個(gè)散亂數(shù)據(jù)點(diǎn)最多由其相鄰的5*5*5個(gè)控制點(diǎn)網(wǎng)格唯一確定，其余控制點(diǎn)網(wǎng)格的貢獻(xiàn)為0。因此，一個(gè)散亂數(shù)據(jù)點(diǎn)的線性插值公式可簡(jiǎn)化為式(9)。

其中 φ(xk－xi)、φ(yk－yj)、φ(zk－zl)分別為x、y、z三個(gè)方向的插值基函數(shù)，即三次B樣條函數(shù)。式(9)是一個(gè)欠定方程，有許多組Cijl滿足式(9)的要求。為了得到Cijl點(diǎn)的值，可以構(gòu)造如下最小平方目標(biāo)函數(shù)Obj，設(shè)其為Cijl點(diǎn)對(duì)函數(shù)f在點(diǎn)(xk，yk，zk)處的真實(shí)貢獻(xiàn)與期望值之差的平方和:

為得到最佳擬合參數(shù)Cijl(共(Xres+4)(Yres+4)(Zres+4)個(gè)參數(shù))，令這個(gè)目標(biāo)函數(shù)最小，即令其對(duì)Cijl的偏導(dǎo)數(shù)為零:

這樣即得到一個(gè)(Xres+4)(Yres+4)(Zres+4)階線性方程組，求解此線性方程組，就可得到規(guī)則控制網(wǎng)格點(diǎn)的Cijl。再根據(jù)式(9)就可計(jì)算出各散亂數(shù)據(jù)點(diǎn)的值f(xk，yk，zk)，甚至可估算出原來沒有數(shù)值的點(diǎn)的值。只需存儲(chǔ)有限個(gè)規(guī)則控制網(wǎng)格點(diǎn)的Cijl數(shù)據(jù)，就代表了原有許多散亂走時(shí)點(diǎn)，達(dá)到了散亂走時(shí)數(shù)據(jù)曲面擬合和壓縮的目的。

3 走時(shí)數(shù)據(jù)存儲(chǔ)

如前所述，每一個(gè)散亂數(shù)據(jù)點(diǎn)最多可由其相鄰的5*5*5個(gè)控制點(diǎn)網(wǎng)格唯一確定，其余控制點(diǎn)網(wǎng)格的貢獻(xiàn)為0。因此，當(dāng)對(duì)式(10)中某一個(gè)Cijl求偏導(dǎo)數(shù)時(shí)，并不是所有控制點(diǎn)網(wǎng)格Cijl都參與計(jì)算，最多只有其周圍的5*5*5個(gè)Cijl參與計(jì)算。這樣在由一系列式(11)組成的線性方程組Ax=b中，每一個(gè)方程都有(Xres+4)(Yres+4)(Zres+4)個(gè)Cijl，但最多只有25組共125個(gè)Cijl的系數(shù)不為零。因此，系數(shù)矩陣A實(shí)際上是一個(gè)稀疏矩陣，具有25個(gè)條帶的帶狀結(jié)構(gòu)，每一個(gè)條帶的帶寬為5，如圖3所示。

圖3 三維走時(shí)數(shù)據(jù)曲面擬合的系數(shù)矩陣Fig．3 Coefficientmatrix of3D travel time data surface fitting

在圖3中，外圍的虛線是對(duì)擴(kuò)展二個(gè)單位邊界的網(wǎng)格點(diǎn)求偏導(dǎo)數(shù)的系數(shù)，這些系數(shù)中有些可能值不為零，有些可能整個(gè)一行值都為零或者是很小的值，這樣就可能造成系數(shù)矩陣A不滿秩。為了保證系數(shù)矩陣A滿秩，同時(shí)保證線性方程組A x=b有穩(wěn)定的解，所以在求解之前，需要在稀疏矩陣A的主對(duì)角線加上一個(gè)較小的阻尼，比如0.000 001。

為了能夠利用較少的存儲(chǔ)器空間儲(chǔ)存稀疏矩陣完整的數(shù)據(jù)，通常采用壓縮存儲(chǔ)的方法，即不存儲(chǔ)或盡量少存儲(chǔ)稀疏矩陣的零元素，同時(shí)給出必要的索引信息，以便可以方便地找到所要用的非零元素在矩陣中的位置。當(dāng)運(yùn)算過程中產(chǎn)生新的非零元素時(shí)，有位置能方便地將其填入。

本算法選用CSC(Compressed Sparse Colum)存儲(chǔ)法來存儲(chǔ)如圖3所示的稀疏矩陣A，該方法是對(duì)稀疏矩陣進(jìn)行逐列壓縮存儲(chǔ)。為了存儲(chǔ)n階矩陣A，假設(shè)矩陣A中共有l(wèi)個(gè)非零元素，則需要一個(gè)數(shù)據(jù)類型與矩陣A相同的l維向量x，按照先列后行的順序，依次存放矩陣A中的非零元素;然后用一個(gè)整形的l維向量x(I)，按照同樣的順序依次存放矩陣A中這些非零元素的行號(hào);最后還需要引入一個(gè)n+1維的整形向量x(R)，來指明矩陣A中第i列中第一個(gè)非零元素被存儲(chǔ)在向量x中的位置。

本算法在求解稀疏矩陣線性方程組時(shí)，采用非對(duì)稱的多波前法。即對(duì)于線性方程組Ax=b，通過矩陣A的LU分解來得到最后的解向量。多波前法(Multifrontal)是求解稀疏矩陣的比較穩(wěn)定和高效的算法。其原理是找出并構(gòu)造稀疏矩陣中的密集子塊(Frontal)，F(xiàn)rontal的分解直接調(diào)用BLAS(Basic Linear Algebra Subprograms)［13～15]。Frontal的生成、消去、裝配、釋放等，都是由特定的數(shù)據(jù)結(jié)構(gòu)來指引。多波前法只存儲(chǔ)非零元素的數(shù)值和位置，存儲(chǔ)效率較高。在同一個(gè)波前，不管含多少列，其結(jié)構(gòu)都是類似的。因此只需記錄最左邊那列的位置信息，其它列共用這些信息。波前越寬，節(jié)省的空間越多。而索引存儲(chǔ)，鏈接表存儲(chǔ)等，都要記錄每列非零元素的數(shù)值和位置，因此效率要低于多波前法。

4 實(shí)際資料算例

為了驗(yàn)證本文方法的可行性，作者首先用三次B樣條函數(shù)最小二乘法曲面擬合算法，對(duì)實(shí)際三維地震射線追蹤走時(shí)表進(jìn)行壓縮;然后將規(guī)則網(wǎng)格控制點(diǎn)的數(shù)值存儲(chǔ)入內(nèi)存，以代表原有的散亂走時(shí)數(shù)據(jù);最后在克希霍夫偏移成像時(shí)，直接在內(nèi)存中對(duì)壓縮后的走時(shí)表進(jìn)行插值解編即可用于成像計(jì)算。

如圖4(a)、圖4(b)、圖4(c)分別是中國(guó)西部某地區(qū)三維射線追蹤走時(shí)表，在x、y和z三個(gè)方向上的三個(gè)切片示意圖。原始走時(shí)表共有756 939個(gè)走時(shí)數(shù)據(jù)，使用三次B樣條函數(shù)擬合壓縮后，用14*14*14個(gè)網(wǎng)格點(diǎn)來表示。圖4(d)為擬合相對(duì)誤差分析，可以看出，擬合的相對(duì)誤差絕對(duì)值最大不超過0.001，大多數(shù)點(diǎn)的擬合誤差絕對(duì)值小于0.000 2，滿足精度的要求。圖4(e)是由原始走時(shí)表得到的克?；舴蚱平Y(jié)果在主測(cè)線方向的一個(gè)切片;圖4(f)是由壓縮走時(shí)表方法得到的相同數(shù)據(jù)體的克?；舴蚱魄衅６N不同的走時(shí)數(shù)據(jù)調(diào)用方法得到的偏移結(jié)果差別不大，都能對(duì)重要構(gòu)造進(jìn)行很好地成像，并且都具有較小的偏移噪音。這說明壓縮走時(shí)表對(duì)偏移結(jié)果精度的影響可以忽略不計(jì)。

通過算例可以看出，采用三次B樣條函數(shù)擬合方法，對(duì)三維射線追蹤走時(shí)表壓縮具有較高的壓縮比。同時(shí)，通過合理地規(guī)劃規(guī)則網(wǎng)格節(jié)點(diǎn)，可以保證壓縮擬合的計(jì)算效率和計(jì)算精度。把這些有限個(gè)規(guī)則網(wǎng)格控制點(diǎn)的數(shù)值存儲(chǔ)于內(nèi)存中，在偏移成像時(shí)，在內(nèi)存中使用三次B樣條函數(shù)的線性插值，即可解編出原有的散亂走時(shí)數(shù)據(jù)。由于省去了頻繁讀取大數(shù)據(jù)量的走時(shí)表文件，所以克?；舴蚱瞥上竦挠?jì)算效率提高了二倍以上。此外，由于走時(shí)表壓縮和解編的計(jì)算精度都很高，所以通過壓縮走時(shí)表來得到的偏移結(jié)果，與常規(guī)的直接讀取走時(shí)表文件得到的偏移結(jié)果相比，具有很高的逼近程度。

關(guān)于壓縮比例需要說明的是，根據(jù)Nyquist采樣定律，當(dāng)采樣頻率fs不小于信號(hào)中最高頻率fmax的二倍，即fs≥2fmax時(shí)，離散信號(hào)采樣能無失真地恢復(fù)到原來的連續(xù)信號(hào)。一般在實(shí)際應(yīng)用中，應(yīng)保證采樣頻率為信號(hào)最高頻率的5倍～10倍才能保證合理的精度。由理論模型分析可得，一個(gè)周期至少需要十個(gè)左右的規(guī)則網(wǎng)格點(diǎn)，才能保證擬合相對(duì)誤差的精度。比如對(duì)于一個(gè)模型，模型數(shù)據(jù)波動(dòng)大約1.33個(gè)周期，用大于14*14*14的三維規(guī)則網(wǎng)格點(diǎn)去擬合壓縮，數(shù)據(jù)波動(dòng)在三個(gè)方向的每個(gè)周期中，大約有14/1.33≈10個(gè)規(guī)則網(wǎng)格點(diǎn)，這樣的擬合壓縮是合適的。如果數(shù)據(jù)波動(dòng)的每個(gè)周期中的規(guī)則網(wǎng)格點(diǎn)數(shù)小于10，這樣的擬合壓縮則不能保證合理的精度，自然是不合適的。

圖4 壓縮走時(shí)表偏移實(shí)例Fig．4 Migration example by compressing travel time data

5 結(jié)論

地震波走時(shí)是克希霍夫疊前深度偏移的必要參數(shù)，在偏移過程中需頻繁調(diào)用各成像點(diǎn)的走時(shí)數(shù)據(jù)。但是，由于三維射線追蹤的走時(shí)表屬于大型散亂數(shù)據(jù)，現(xiàn)有的克?；舴蚱瞥上穹椒ㄊ菍⒆邥r(shí)表作為文件來存儲(chǔ)。由于需要頻繁讀取大數(shù)據(jù)量的走時(shí)表文件，所以計(jì)算效率不盡人意。為了提高克?；舴蚱瞥上竦挠?jì)算效率，而不至于對(duì)計(jì)算精度造成不利影響，作者在本文中采用了一種先對(duì)三維射線追蹤走時(shí)表進(jìn)行壓縮，再利用壓縮后的走時(shí)數(shù)據(jù)進(jìn)行偏移的方法。具體做法是:首先規(guī)劃一個(gè)具有規(guī)則網(wǎng)格控制點(diǎn)，并包含所有散亂數(shù)據(jù)點(diǎn)集的最小盛放區(qū)域;然后利用三次B樣條函數(shù)，對(duì)散亂的走時(shí)數(shù)據(jù)進(jìn)行最小二乘曲面擬合，從而求出各規(guī)則網(wǎng)格控制點(diǎn)的數(shù)值;最后只需要把這些有限個(gè)規(guī)則網(wǎng)格控制點(diǎn)的數(shù)值，以數(shù)組形式存儲(chǔ)于內(nèi)存中。在偏移成像時(shí)，再由這些規(guī)則網(wǎng)格控制點(diǎn)的數(shù)值，使用B樣條函數(shù)的線性插值公式解編出走時(shí)表。

在計(jì)算時(shí)充分考慮本算法最小二乘法線性方程組的系數(shù)矩陣為稀疏矩陣，并且具有非零元素規(guī)則排列的特征。因此為了進(jìn)一步節(jié)省內(nèi)存空間，同時(shí)提高存儲(chǔ)效率，在矩陣存儲(chǔ)時(shí)作者采用了稀疏矩陣存儲(chǔ)方法，并使用基于多波前法的LU分解求解稀疏矩陣線性方程組。實(shí)際地震資料算例驗(yàn)證了本三維走時(shí)表壓縮算法的可行性，計(jì)算的快速性以及控制節(jié)點(diǎn)的穩(wěn)定性。對(duì)射線追蹤走時(shí)表進(jìn)行科學(xué)合理的壓縮存儲(chǔ)，省去了頻繁地進(jìn)行大數(shù)據(jù)量文件的讀寫操作，促使三維克?；舴蚱瞥上竦挠?jì)算效率提高了二倍以上，并且對(duì)偏移結(jié)果精度的影響還可以忽略不計(jì)。對(duì)射線追蹤走時(shí)表進(jìn)行的科學(xué)合理的壓縮存儲(chǔ)方法，可用于提高其它地震正演模擬，以及地震層析成像等的計(jì)算效率。

［1]郭洪升，陳俊良．地震數(shù)據(jù)實(shí)時(shí)自適應(yīng)壓縮方法研究［J]．地震學(xué)報(bào)，1989，11(1):68．

［2]張軍華，仝兆岐．用小波變換法定量壓縮地震數(shù)據(jù)［J]．石油大學(xué)學(xué)報(bào):自然科學(xué)版 ，2003，27(5):29．

［3]王喜珍，滕云田，高孟潭，等．基于整數(shù)小波變換的地震數(shù)據(jù)壓縮［J]．地震學(xué)報(bào)，2004，26(增刊):116．

［4]余平，馬小虎，陳恒金．基于提升小波的地震數(shù)據(jù)壓縮編碼算法［J]．蘇州大學(xué)學(xué)報(bào):工科版，2009，29(1):7．

［5]AVERBUCH A Z，METER F，STROMBERG J O，et al．Vassiliou，Low bit－rate efficient compression for seismic data［J]．Institute of Electrical and Electronics Engineers Transactions on Image Processing，2001，10:1801．

［6]BERNASCONIG，VASSALLO M．Efficient data compression for seismic－while－drilling applications［J]．IEEE Transactions on Geoscience and Remote Sensing，2003，41:687．

［7]WU R S，WANG Y．Seismic data compression using adapted local cosine transform and its effect on imaging，61st Meeting［E]．European Association of Geoscientists and Engineers，Extended Abstracts，2003．

［8]GIANCARLO B，VITTORIO R．High－quality compression of MWD signals［J]．Geophysics，2004，69:849．

［9]CABDèSE J，DEMANET L．The curvelet representation of wave propagators is optimally sparse［J]．Communications on Pure and Applied Mathematics，2005，58:1472．

［10]CHENG G，ZHANG B．Compression storage and solution of large sparse matrix［J]．Progress in Geophysics，2008，23:674．

［11]關(guān)履泰，覃廉，張?。脜?shù)樣條插值挖補(bǔ)方法進(jìn)行大規(guī)模散亂數(shù)據(jù)曲面造型［J]．計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2006，18(3):372．

［12]毛可飛，路輝．基于多層B樣條的海底地形生成方法［J]．計(jì)算機(jī)仿真，2005，22(4):222．

［13]夏愛生，胡寶安，申楠公，等．系數(shù)矩陣為帶狀大型稀疏矩陣線性方程組解法研究［J]．軍事交通學(xué)院學(xué)報(bào)，2007，9(1):83．

［14]成谷，張寶金．反射地震走時(shí)層析成像中的大型稀疏矩陣壓縮存儲(chǔ)和求解［J]．地球物理學(xué)進(jìn)展，2008，23(3):674．

［15]陳英時(shí)，吳文勇．采用多波前法求解大型結(jié)構(gòu)方程組［J]．建筑結(jié)構(gòu)，2007，2(9):138．

P 631.4

A

1001—1749(2011)02—0122—07

教育部高等教育博士點(diǎn)專項(xiàng)研究基金資助(200804251502);中國(guó)石油天然氣集團(tuán)公司創(chuàng)新基金資助(07E1019);中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金資助(09CX04013A)

2010－09－10改回日期:2010－12－13

王珺(1973－)，女，山東東營(yíng)人，博士，講師，現(xiàn)于中國(guó)石油大學(xué)信息與控制工程學(xué)院任教。