王 珺，楊長春，李波濤

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

壓縮三維走時表提高克?；舴蚱朴嬎阈?/p>

王 珺1，楊長春2，李波濤3

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

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

克希霍夫偏移;計算效率;三維走時表壓縮

0 前言

用射線追蹤方法計算地震波走時是克?；舴虔B前深度偏移的基礎(chǔ)。在三維地震克?；舴蚱瞥上襁^程中，三維射線追蹤的走時表往往具有相當(dāng)大的數(shù)據(jù)量，通常作為專門的文件來存儲，成像時需頻繁地讀取此走時表文件，這樣勢必造成成像計算效率不高。如果能夠?qū)⒋髷?shù)據(jù)量的三維射線追蹤的走時表壓縮存儲，用有限個數(shù)值點來代替，直接以數(shù)組的形式存儲在內(nèi)存中，就省去了頻繁地進(jìn)行較大數(shù)據(jù)量文件的讀寫操作，自然會提高克?；舴蚱瞥上竦挠嬎阈?。對射線追蹤走時表進(jìn)行的科學(xué)合理的壓縮存儲方法，還可用于提高正演模擬以及層析成像等的計算效率，具有重大的理論和現(xiàn)實意義。目前，國內(nèi)、外有一些關(guān)于常規(guī)地震資料的壓縮方法的研究工作(國內(nèi)有郭洪升，陳俊良［1]、張軍華，仝兆岐［2]、王喜珍等［3]、余平［4];國外有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ù)壓縮，實際上就是采用基于變換的方法，比如小波變換方法，正交變換方法等。

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

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

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

1 Kirchhoff偏移基本原理

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

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

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

入射波場見式(3)。

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

2 三維走時表壓縮

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2．2 曲面體擬合

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

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

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

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

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

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

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

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

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

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

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

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

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

本算法在求解稀疏矩陣線性方程組時，采用非對稱的多波前法。即對于線性方程組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)來指引。多波前法只存儲非零元素的數(shù)值和位置，存儲效率較高。在同一個波前，不管含多少列，其結(jié)構(gòu)都是類似的。因此只需記錄最左邊那列的位置信息，其它列共用這些信息。波前越寬，節(jié)省的空間越多。而索引存儲，鏈接表存儲等，都要記錄每列非零元素的數(shù)值和位置，因此效率要低于多波前法。

4 實際資料算例

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

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

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

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

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

5 結(jié)論

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

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

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

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

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

［4]余平，馬小虎，陳恒金．基于提升小波的地震數(shù)據(jù)壓縮編碼算法［J]．蘇州大學(xué)學(xué)報:工科版，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ī)輔助設(shè)計與圖形學(xué)學(xué)報，2006，18(3):372．

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

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

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

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

P 631.4

A

1001—1749(2011)02—0122—07

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

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

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