高 瀟，張 偉，徐錦承，楊 輝，邵理陽，潘 權(quán)，宋章啟

(1.中國科學技術(shù)大學 地球和空間科學學院，合肥 230000;2.南方科技大學 a.地球與空間科學系，b.電子與電氣工程系，c.創(chuàng)新創(chuàng)業(yè)學院，深圳 518055)

0 引言

隨著地震勘探規(guī)模的擴大，海量的地震數(shù)據(jù)需要被存儲以供后期處理使用。盡管這個問題可以通過使用大容量的地震存儲設(shè)備來解決，但是在實際處理中計算機的IO瓶頸又會限制處理速度，從而出現(xiàn)計算機IO能力與計算能力不匹配的現(xiàn)象。因此，將數(shù)據(jù)壓縮之后保存，不僅能解決海量數(shù)據(jù)的存儲問題，同時也能減輕IO負擔，加快處理速度。

目前的數(shù)據(jù)壓縮方法通常被分為兩大類：①無損壓縮;②有損壓縮。對于無損壓縮，數(shù)據(jù)壓縮之后能完全恢復(fù)，不丟失任何信息，但隨之而來的問題就是無法實現(xiàn)較高的壓縮率。而有損壓縮的壓縮過程一般都分為變換-量化-編碼三個過程[1-3]。相對無損壓縮，有損壓縮增加了變換、量化這兩個步驟：變換可以有效地消除數(shù)據(jù)間的相關(guān)性，數(shù)據(jù)間相關(guān)性越大，數(shù)據(jù)的冗余度就越大；量化是將主要變換系數(shù)和次要變換系數(shù)區(qū)分開來，主要系數(shù)包含數(shù)據(jù)的概貌部分，是必須保留的，否則無法重構(gòu)數(shù)據(jù)，次要系數(shù)重構(gòu)數(shù)據(jù)的細節(jié)部分，是對重構(gòu)數(shù)據(jù)質(zhì)量的進一步提高，有損壓縮就是通過舍棄部分次要系數(shù)來實現(xiàn)高壓縮率的。雖然壓縮前、后數(shù)據(jù)會出現(xiàn)一些差異，但通過設(shè)置合適的壓縮參數(shù)，這些差異就會比較微弱，并且丟失的是數(shù)據(jù)中的非主要信息，在很多實際應(yīng)用中這種誤差是完全可以忽略的。

對于地震數(shù)據(jù)的壓縮算法研究，早期開始于傳統(tǒng)的無損壓縮算法(如Huffman編碼壓縮[4]，LZ77編碼壓縮[5]或LZW編碼壓縮[6]等)，之后為了規(guī)范地震數(shù)據(jù)壓縮后的交換格式，國際數(shù)字地震臺網(wǎng)聯(lián)合會(The International Federation of Digital Seismograph Networks, FDSN)，將Steim壓縮算法確定為地震數(shù)據(jù)交換標準格式的核心算法，Steim算法目前常用的主要是兩種:①SEED數(shù)據(jù)格式所采用的Steim1算法;②miniSEED數(shù)據(jù)格式所采用的Steim2 算法，但這些無損壓縮算法用于地震數(shù)據(jù)能實現(xiàn)的最大壓縮率有限，不能滿足實際海量數(shù)據(jù)的存儲需求。而有損壓縮算法根據(jù)研究目的不同，用戶可以自行設(shè)定壓縮參數(shù)從而使壓縮率達到10乃至超過100[7-10]。正因為如此，在過去幾十年，許多地球物理學家都致力于尋找一種高效的地震數(shù)據(jù)有損壓縮算法。當前比較流行的地震數(shù)據(jù)壓縮算法大都是基于變換的壓縮算法，如Candes等[11]提出的多尺度、多方向的各向異性變換方法——Curvelet變換；Pennec等[12]提出的Bandelet變換；Guo等[13]提出的Shearlet變換以及Fomel等[14]提出的Seislet變換等，而目前使用最多的兩個典型代表就是：離散余弦變換[9,15-18]和小波變換[6,10,19-26]。

Ahmed[15]在JPEG圖像壓縮算法中首次提出了離散余弦變換的概念；Spanias[17]對比了四種基于變換的壓縮方法對地震數(shù)據(jù)壓縮的效果，結(jié)果表明，DCT是最健壯且有效的；劉財[18]在DCT的基礎(chǔ)上，將游程編碼算法引入到了地震數(shù)據(jù)壓縮中，測試顯示這種技術(shù)對層間相似性較強的數(shù)據(jù)具有不錯的壓縮效果；Dalmau[9]為了解決磁盤的IO限制，將DCT應(yīng)用到了時間域逆時偏移的互相關(guān)處理中。然而，由于該算法自身的特性，即將數(shù)據(jù)分塊處理，當實現(xiàn)較高的壓縮率時，解壓得到的數(shù)據(jù)塊之間存在明顯的邊界，這個即是所謂的“塊狀效應(yīng)”[7]。小波變換是一種用于圖像層次分解的數(shù)學工具[27]。相比于離散余弦變換,由于小波變換引入了不同的母小波，這些母小波可以通過擴張或收縮以及平移來構(gòu)造最適合的子波，從而更易于獲得最佳的壓縮效果。另外，由于不需要對數(shù)據(jù)進行分塊處理，從而使小波變換壓縮得到的數(shù)據(jù)即使在高壓縮率的情況下，也不會產(chǎn)生”塊狀效應(yīng)”[28]。唐向宏等[29]以二維地震數(shù)據(jù)具有較強的分頻帶特性為切入點，提出了一種利用M帶小波變換進行壓縮的算法，該算法能在保證重建數(shù)據(jù)精度滿足后續(xù)處理的前提下，獲得10～20的壓縮比；徐鋒濤等[30]提出了一種基于嵌入式零樹小波編碼(EZW)的地震數(shù)據(jù)壓縮算法，其首先對地震數(shù)據(jù)進行二維小波變換，之后利用系數(shù)之間的相關(guān)性對變換后的系數(shù)采用嵌入式零樹小波進行編碼，最后使用自適應(yīng)算術(shù)編碼對量化后的系數(shù)進行無損熵編碼。但是小波變換不能很好地表示數(shù)據(jù)中的震蕩模式，而地震數(shù)據(jù)又恰好具備這種特性[7]。

近年來，除了上述的一些經(jīng)典地震數(shù)據(jù)壓縮方法外，一些新的方法也被提出。耿瑜等[31]將Dreamlet變換應(yīng)用到了地震數(shù)據(jù)壓縮中；Dreamlet變換中利用到了局部余弦基和頻散關(guān)系，研究結(jié)果表明，對于地震數(shù)據(jù)壓縮而言，Dreamlet變換比Curvelet變換具有更好的性能；Boehm[32]提出了一組專門針對全波形反演中伴隨狀態(tài)方法的波場壓縮策略，包含粗空間網(wǎng)格重新插值，自適應(yīng)選擇多項式階數(shù)和自適應(yīng)選擇浮點數(shù)精度這三個步驟。盡管這種壓縮方法能獲得比較理想的壓縮率，但它其中的一個壓縮策略是專門針對譜元法設(shè)計的，因此如果要在其他方法中使用(有限差分法)，就需要額外的適配工作。Lindstrom[10]將zfp壓縮算法(zfp算法是一個開源的用于壓縮浮點數(shù)組的C/C++函數(shù)庫，由勞倫斯利弗莫爾國家實驗室的Peter Lindstrom, Markus Salasoo和Matt Larsen開發(fā)并維護)用到了三維全空間波形層析成像的散射積分計算實現(xiàn)中，從而大幅減少了數(shù)據(jù)存儲量并顯著提高了IO性能；陳祖斌等[33]將壓縮感知算法應(yīng)用到了地震數(shù)據(jù)壓縮重構(gòu)中，該方案在下機位實時使用混沌伯努利測量矩陣壓縮小波變換后的系數(shù)，同時在上機位結(jié)合使用貝葉斯小波樹結(jié)構(gòu)壓縮感知重構(gòu)算法和馬爾科夫鏈蒙特卡洛推理對數(shù)據(jù)進行解壓；實際數(shù)據(jù)測試表明這種方法具有不錯的壓縮效果；Nuha[34]以神經(jīng)網(wǎng)絡(luò)研究為基礎(chǔ)，提出了一種適用于檢波器地震數(shù)據(jù)的實時壓縮算法，與常規(guī)方法相比，這種方法在壓縮率和重構(gòu)數(shù)據(jù)質(zhì)量上都要更優(yōu)。

Squeeze壓縮算法(SZ)，是由美國阿貢國家實驗室的Di和Cappello[35-36]開發(fā)維護的一種誤差約束的有損壓縮算法，即數(shù)據(jù)在壓縮過程中會嚴格遵循用戶定義的誤差。SZ最初的研發(fā)目的是為了提升高性能計算中的check-point性能以及數(shù)據(jù)的后續(xù)處理效率。SZ壓縮算法相比之前傳統(tǒng)的壓縮算法，具有壓縮性能高、擴展性好以及移植性佳的優(yōu)點，并且時間復(fù)雜度只有O(N)。為了對比SZ壓縮算法和其他高性能計算中常用的壓縮算法之間的性能差異，Di和Cappello以及Tao對13個涵蓋了各個學科領(lǐng)域的高性能計算應(yīng)用進行了測試[35-38]。實驗結(jié)果表明對任何一個HPC應(yīng)用得到的數(shù)據(jù)集進行壓縮，SZ的性能都是最佳的。此外SZ支持多種語言，包括C、Fortran以及Java。

SZ的壓縮實現(xiàn)十分簡潔，只需在代碼中添加幾行壓縮函數(shù)調(diào)用，就能實現(xiàn)數(shù)據(jù)的高效壓縮存儲。但目前還沒有SZ在地震數(shù)據(jù)壓縮中的相關(guān)工作。筆者首次將SZ壓縮算法引入到地震數(shù)據(jù)壓縮中，首先對Squeeze的算法實現(xiàn)進行闡述，之后對其與實際工作結(jié)合的編碼進行簡單示例，接著通過三個測試案例比較Squeeze算法與DCT壓縮和Dreamlet壓縮的性能差異，最后給出相應(yīng)的結(jié)論。

1 Squeeze算法的原理

SZ的算法實現(xiàn)包括三個步驟：①基于多層和多維數(shù)據(jù)的預(yù)測模型；②自適應(yīng)誤差約束量化；③熵編碼模型。

1.1 多維和多層預(yù)測模型

與傳統(tǒng)基于變換的壓縮方法不同，SZ基于多項式預(yù)測。SZ算法的早期版本[35]使用的是一維預(yù)測算法，包括常量、線性以及二次預(yù)測模型。之后的版本[37-38]將預(yù)測模型改進成了多維，從數(shù)學意義上來說，稱之為Lorenzo預(yù)測，這一改變使數(shù)據(jù)的預(yù)測命中率相比之前有了顯著提高。

圖1 二維多層預(yù)測模型示例

k1,k2≤n}{i0,j0}

(1)

V(i0-k1,j0-k2)

(2)

f(i0,j0)=2V(i0-1,j0)+2V(i0,j0-1)-

4V(i0-1,j0-1)-V(i0-2,j0)-

V(i0,j0-2)+2V(i0-2,j0-1)+

2V(i0-1,j0-2)-V(i0-2,j0-2)

(3)

類似推導(dǎo)可以得到多維多層的預(yù)測模型，

V(x1-k1,…,xd-kd)

(4)

其中：d為預(yù)測模型中數(shù)據(jù)集的維度；n為預(yù)測所用的層數(shù)。對于一個給定的數(shù)據(jù)集，參數(shù)d是固定的，但是n是可變的。從一方面來說，增加預(yù)測層數(shù)可以獲得比較好的預(yù)測精度，但從另一方面而言，使用的層數(shù)越多，引入的不相關(guān)信息也越多，從而在某種程度上也會降低預(yù)測精度。因此，壓縮性能的好壞與預(yù)測層數(shù)n的選擇有很大關(guān)系，需要針對不同問題測試選擇最優(yōu)的n值。

1.2 自適應(yīng)誤差約束量化

經(jīng)過上面的多維和多層預(yù)測之后，需要將預(yù)測值在誤差約束條件下進行量化。SZ算法的誤差約束是指允許用戶通過配置壓縮參數(shù)來對壓縮結(jié)果的誤差進行控制。在SZ中，用戶可以選擇三種誤差約束方式中的一種用于當前的工作：

1)絕對誤差約束(ABS)。每一個數(shù)據(jù)點的壓縮誤差都被限定在同一個最大允許誤差范圍(Max Error, ME)內(nèi)。例如，設(shè)定誤差約束ebabs=0.000 1，數(shù)據(jù)點的原始值為V，那么該點解壓之后的值必須在[V-0.000 1,V+0.000 1]范圍內(nèi)。

2)相對誤差約束(REL)。也稱為基于數(shù)據(jù)范圍的相對誤差約束，會在考慮數(shù)據(jù)集范圍的基礎(chǔ)上，得到一個誤差值。例如，設(shè)定ebrel=0.01，原始數(shù)據(jù)集是[100，101，102，…，108，109，110]，那么，由于數(shù)據(jù)集范圍是110-100=10，每一個數(shù)據(jù)點的最大誤差就是10*0.01=0.1。

3)逐點相對誤差約束(PW_REL)。類似REL，但是它不考慮整體的數(shù)據(jù)范圍，而是將當前點的數(shù)值作為一個數(shù)據(jù)范圍。例如，設(shè)定ebpw=0.01，那么對每一個點的數(shù)據(jù)ME就是0.01*V。

這三種壓縮設(shè)置都有各自的優(yōu)缺點，實際處理中選用哪一種需要結(jié)合具體的研究問題而定，并且用戶可以自行調(diào)優(yōu)以獲得最適合的壓縮方式。在本文中，由于地震數(shù)據(jù)變化范圍較大，使用絕對誤差約束方式和逐點相對誤差方式都不能獲得比較好的壓縮效果，因此我們使用第二種誤差約束——相對誤差約束(REL)進行壓縮。

選定誤差約束方式后，就可以根據(jù)相應(yīng)的誤差約束條件計算出最大允許誤差ME值，之后的量化操作以ME值為基礎(chǔ)，其具體又可以分為兩個過程。

1)首先在第一步多層預(yù)測的基礎(chǔ)上，獲得當前點的預(yù)測值，在算法中稱之為“第一階段預(yù)測值”，如圖2中的紅點所示。設(shè)定量化編碼位數(shù)為m，則總共可支持2m-1個量化編碼(編碼“0”有其他用途)。以“第一階段預(yù)測值”為中心，向兩端各增/減1倍ME值，從而獲得一個“基礎(chǔ)預(yù)測范圍”。之后以該“基礎(chǔ)預(yù)測范圍”為單位(2*ME)，向上/向下各線性延拓次，從而共獲得2m-1個“預(yù)測范圍”，從值最小的預(yù)測范圍開始以“1”為起始序號編碼，相應(yīng)的編碼值如圖2右半部分所示。

圖2 在誤差約束下的量化過程示意圖

2)當真實值落入某一個預(yù)測范圍內(nèi)，就用該預(yù)測范圍的中心值表示，稱其為“第二階段預(yù)測值”。不難看出，“第二階段預(yù)測值”與真實值(圖2中藍點)之間的誤差一定小于ME。在實際壓縮中，保存的是“第二階段預(yù)測值”所在區(qū)間對應(yīng)的編碼值。

若當前處理點的真實值能落入到擴展出來的間隔內(nèi)，SZ算法就標記該點是可預(yù)測的，同時用相應(yīng)的“第二階段預(yù)測值”表示解壓縮后的數(shù)值，用相應(yīng)的量化編碼表示壓縮的數(shù)據(jù)。但是如果原始值不能落入擴展出來的間隔內(nèi)，那么SZ算法將標記這個點是不可預(yù)測的，并將該點的量化編碼標記為“0”，之后采用另外一套壓縮算法進行處理——二進制表示分析壓縮算法。如前所述，這個量化過程共需要m位二進制來編碼所有的數(shù)據(jù)。SZ算法允許用戶自定義m的大小，若在壓縮過程中不顯式指定，算法自身也會根據(jù)數(shù)據(jù)集自適應(yīng)地選擇合適的m值進行編碼。

1.3 熵編碼

經(jīng)過以上預(yù)測和量化兩個步驟，數(shù)據(jù)轉(zhuǎn)換成了采用位數(shù)有限的量化編碼表示，實現(xiàn)了有損壓縮。為了進一步提高壓縮效果，SZ算法在最后步驟引入了熵編碼。這是一種無損壓縮算法，基本思想就是將出現(xiàn)次數(shù)越多的數(shù)據(jù)項用位數(shù)越少的編碼項代替。在具體實現(xiàn)中，使用的是經(jīng)過優(yōu)化的哈夫曼編碼算法。

在當前常見的壓縮算法中，都會使用熵編碼來提高壓縮率，筆者在之后數(shù)值測試部分所使用的DCT和Dreamelt壓縮算法，在最后的操作步驟也使用了標準JPEG算法所采用的游程編碼和哈夫曼編碼來實現(xiàn)熵編碼操作。

2 在地震波模擬中的應(yīng)用

Squeeze算法相對于其他壓縮算法而言，簡單易用是其一個顯著特點。由于底層實現(xiàn)已經(jīng)完全封裝好，只需要調(diào)用相應(yīng)的接口就能實現(xiàn)數(shù)據(jù)的高效壓縮，地球物理學家不需要花費很多精力去學習壓縮算法原理，可以將精力集中在要解決的科學問題上。

將Squeeze算法與地震波正演模擬程序相結(jié)合，其總體實現(xiàn)流程如表1所示(以Fortran代碼為例)。從表1中可以看出，只需要在工作中添加幾行代碼，就能實現(xiàn)一個比較理想的壓縮方案。

表1 Squeeze算法與地震波模擬程序結(jié)合示例

3 數(shù)值測試

對SZ的壓縮效果進行測試，研究壓縮中一些典型參數(shù)的取值方法。壓縮效果主要通過壓縮率(CR)和壓縮誤差(CE)來進行衡量，定義如下：

CR=原始數(shù)據(jù)大小/壓縮數(shù)據(jù)大小

(5)

CE=abs(原始數(shù)據(jù)數(shù)值-解壓數(shù)據(jù)數(shù)值)

(6)

3.1 二維地震記錄壓縮

筆者先對SZ在地震記錄數(shù)據(jù)上的壓縮效果進行測試，所用數(shù)據(jù)是疊后二維SEG-EAGE鹽丘的地震記錄數(shù)據(jù)(圖3)。同時，與離散余弦變換壓縮和Dreamlet變換壓縮進行效果對比。與SZ壓縮不同，基于變換的DCT壓縮和Dreamlet壓縮不是誤差約束的，它們都是通過控制變換系數(shù)的保留比例來實現(xiàn)數(shù)據(jù)壓縮。一般保留比例用參數(shù)rt表示，具體實現(xiàn)為：原始數(shù)據(jù)經(jīng)過變換后將獲得變換后的系數(shù)，如果沒丟失其中的系數(shù)，數(shù)據(jù)能夠完全變換回原始數(shù)據(jù)，但為了獲得較高的壓縮性能，勢必需要丟棄部分系數(shù)，定義絕對值最大的變換系數(shù)和rt的乘積作為閾值，當系數(shù)的絕對值小于這個閾值時，就會被抹零。

圖3 2維SEG-EAGE疊后鹽丘地震記錄數(shù)據(jù)

測試了這三種方法在不同的壓縮設(shè)置下的壓縮效果，如表2所示。同時為了衡量數(shù)據(jù)集的平均壓縮誤差大小，引入了峰值信噪比系數(shù)(Peak Signal-to-Noise Ratio,PSNR)這一指標，數(shù)學定義為式(7)。

表2 不同壓縮設(shè)置下的離散余弦變換，Dreamlet變換和SZ的壓縮率以及PSNR

(7)

其中：Rx是數(shù)據(jù)范圍，Rx=xmax-xmin；rmse表示均方根誤差，定義為式(8)。

(8)

一般來說，rmse數(shù)值越大或PSNR數(shù)值越小就表示壓縮誤差越大，反之亦然。圖4給出了這三種方法壓縮率和PSNR之間的關(guān)系曲線(部分)。從圖4中可以清楚地看到，在相同壓縮率的情況下，Dreamlet壓縮的PSNR最小，其次是DCT，SZ算法的PNSR最大，從而說明其誤差相對較小。圖5是三種方法在壓縮率接近的情況下，解壓縮之后的數(shù)據(jù)，其中圖5(a)是壓縮設(shè)置為rt=5.0E-4情況下的DCT解壓縮結(jié)果，其壓縮率CR=9.578；圖5(b)是壓縮設(shè)置為rt=3.0E-2情況下的Dreamlet解壓縮結(jié)果，其壓縮率CR=9.710；圖5(c)是壓縮設(shè)置為rel=5.0E-5情況下的SZ解壓縮結(jié)果，其壓縮率為CR=9.778。從圖5中可以看出，Dreamlet變換的解壓縮結(jié)果相對較差，在600左部和1 200右部丟失了許多重要的細節(jié)信息，相比之下，DCT和SZ的解壓縮結(jié)果要好很多，且在視覺上整體差異不大。

圖4 DCT, Dreamlet和SZ三種壓縮方法壓縮率與PSNR之間的關(guān)系

為了更加確切地比較三者的差異大小，在圖5的基礎(chǔ)上，獲得了它們各自的壓縮誤差，如圖6所示。其中DCT的最大壓縮誤差絕對值為0.002 679，Dreamlet壓縮為0.191 113，而SZ壓縮為0.000 033，這比DCT壓縮誤差小2個數(shù)量級，比Dreamlet壓縮誤差小4個數(shù)量級。為了對比明顯，將圖6的顯示刻度范圍都限定在[-1.5*10-3,1.5*10-3]，從中可以看出，在壓縮率為9.5左右時，SZ和DCT的壓縮誤差相對Dreamelet方法都比較小，并且SZ更略優(yōu)一些。

圖5 DCT、Dreamlet和Squeeze的壓縮數(shù)據(jù)對比

圖6 DCT、Dreamlet和Squeeze的壓縮誤差對比

3.2 三維實際地震道集記錄壓縮

我們采用實際原始道集記錄對Squeeze算法的壓縮能力進行測試，同樣也與Dreamlet和DCT壓縮算法對比測試結(jié)果。測試使用的道集記錄數(shù)據(jù)共有432炮，每一炮都有12條測線接收數(shù)據(jù)，每條測線有檢波器120個，間距為50 m，單個檢波器記錄時間為6 s,采樣間隔為4 ms。

表3展示了壓縮率約為8.5時，三種方法的壓縮參數(shù)設(shè)置。為了對比數(shù)據(jù)誤差，我們在這里隨機選擇第27炮的道集數(shù)據(jù)進行分析，圖7是其所有測線原始道集記錄，圖8分別是DCT、Dreamlet和SZ算法的壓縮重構(gòu)結(jié)果。從圖8中可以直觀地看出，對于DCT和Dreamlet壓縮，它們的重構(gòu)數(shù)據(jù)在數(shù)據(jù)采集的初始階段(例如t<500 ms)以及末尾階段(例如t>4 000 ms)的誤差比較大，同時由于這兩個算法都是分塊操作的，因此重構(gòu)得到的數(shù)據(jù)中會有比較明顯的分塊特征，誤差呈現(xiàn)比較明顯的塊狀分布；SZ重構(gòu)得到的數(shù)據(jù)誤差除了在數(shù)據(jù)采集的初始階段會有比較明顯的誤差外，之后階段都能比較好地保留原始數(shù)據(jù)特征。圖9是它們?nèi)N方法的重構(gòu)誤差對比，DCT和Dreamlet的誤差比較接近，它們在t>1 600 ms區(qū)域丟失了比較多的細節(jié)信息，而SZ誤差相對小很多，未見有明顯的細節(jié)信息丟失，從而驗證了SZ算法用于原始道集記錄壓縮的可行性。

圖9 壓縮率約為8.5DCT、Dreamlet和Squeeze的重構(gòu)誤差對比

表3 壓縮率約為8.5時的離散余弦變換、Dreamlet變換和Squeeze算法的壓縮設(shè)置

圖7 第27炮的所有測線原始地震記錄數(shù)據(jù)(AGC時間窗=500 ms)

3.3 全空間波場數(shù)據(jù)壓縮測試

為了進一步測試Squeeze算法對波場數(shù)據(jù)的壓縮效果，使用有限差分方法[39-41]對Marmousi模型(部分)進行正演波場模擬，之后對獲得的全空間波場數(shù)據(jù)進行壓縮。具體操作為每模擬一個時間步，就將相應(yīng)的全空間三維波場數(shù)據(jù)壓縮并輸出。為了對比壓縮效果，同樣使用DCT，Dreamlet和SZ三種算法對數(shù)據(jù)進行操作。正演模擬相關(guān)參數(shù)如表4所示，相應(yīng)的Marmousi模型(切面)如圖10所示，實際使用的三維模型是該切面沿著y方向擴展得到的。

表4 Marmousi模型的正演模擬參數(shù)設(shè)置

圖10 Marmousi模型(部分)的P波速度切面(y=250 m)

在對三種壓縮方法各自測試了不同的壓縮設(shè)置后，獲得了三種方法壓縮率大致相同的情況，如表5所列。

表5 相近壓縮率下的離散余弦變換，Dreamlet變換和Squeeze算法的壓縮設(shè)置

為了對比壓縮效果，使用正應(yīng)力Txx分量在t=0.6 s、y=250 m處的波場快照進行分析。原始波場快照如圖11所示，相應(yīng)的壓縮之后的數(shù)據(jù)如圖12所示，從圖12中可以明顯看出，三種方法的壓縮都沒有使原始數(shù)據(jù)發(fā)生比較大的差異。圖13則展示了具體的細節(jié)差異。從圖13中可以直觀看出，Dreamlet的壓縮誤差相對較大，其次是DCT，Squeeze的誤差相對較小。這個結(jié)果同前面測試相同，表明Squeeze算法對地震波場數(shù)據(jù)壓縮也有比較好的適用性。

圖11 正應(yīng)力Txx分量在t=0.6 s、y=250 m處的原始數(shù)據(jù)波場快照

圖12 壓縮率相近DCT、Dreamlet和Squeeze的波場快照壓縮數(shù)據(jù)對比的原始數(shù)據(jù)波場快照

圖13 壓縮率相近DCT、Dreamlet和Squeeze的波場快照壓縮誤差對比

基于變換的DCT和Dreamlet壓縮算法由于在變換過程中涉及到數(shù)據(jù)的多維變換操作，其中包含大量的指數(shù)以及乘法運算，其壓縮過程的時間復(fù)雜度不是O(N)，而SZ算法的時間復(fù)雜度只有O(N)，并且其預(yù)測步驟也只有簡單的一些加法和乘法運算(式(3))，量化步驟也沒有涉及到復(fù)雜運算，因此SZ的計算量以及計算效率從理論上來說會比前兩種壓縮方法好。實際測試結(jié)果也驗證了這一點，在如上4 000時間步正演波場壓縮過程中，DCT平均壓縮耗時為0.42 s/步，Dreamlet平均壓縮耗時為0.45 s/步，而SZ平均壓縮耗時為0.31 s/步。對于解壓縮過程，DCT平均0.29 s/步，Dreamlet平均0.32 s/步，SZ平均0.19 s/步，這說明利用SZ進行數(shù)據(jù)壓縮的處理效率是相對較高的。

為了探究由于壓縮造成的誤差對地震波場研究的影響，我們在此以敏感核[42]計算為例進行分析，因為敏感核計算需要使用到波場數(shù)據(jù)。圖14是利用原始地震波場數(shù)據(jù)計算得到的振幅衰減敏感核，圖15分別是利用DCT、Dreamlet以及SZ算法壓縮得到的波場數(shù)據(jù)計算得到的敏感核，從圖15中可以看出，三種方法得到的敏感核形態(tài)基本與原始數(shù)據(jù)計算得到的相同，圖16分別是它們各自的誤差，從圖16中可以看出，DCT和Dreamlet的誤差比較接近，但是SZ的誤差相對它們兩個要更小一些。這說明SZ能在實現(xiàn)較高壓縮率的同時保證后續(xù)處理的數(shù)據(jù)精度。

圖14 原始波場數(shù)據(jù)計算得到的振幅衰減敏感核

圖15 壓縮率相近DCT、Dreamlet和Squeeze壓縮數(shù)據(jù)計算得到的振幅衰減敏感核對比

圖16 壓縮率相近DCT、Dreamlet和Squeeze壓縮數(shù)據(jù)計算得到的振幅衰減敏感核誤差對比

4 結(jié)論

筆者對Squeeze算法原理進行了簡要概述，并將其應(yīng)用到了疊后數(shù)據(jù)，實際道集記錄數(shù)據(jù)和波場數(shù)據(jù)的壓縮實現(xiàn)中。通過與基于變換的離散余弦變換壓縮和Dreamlet變換壓縮的測試結(jié)果對比，SZ算法在相同壓縮率的前提下，可以在較大程度上保留原始數(shù)據(jù)特征，實現(xiàn)較小的壓縮誤差，同時對后續(xù)處理造成的誤差影響也是相對較小的，并且SZ算法的壓縮和解壓縮時間開銷也很低。另外，SZ算法由于使用的便捷性，其在實際工作中的編碼復(fù)雜度在這三種方法中也是較為簡單的。