基于DCT的陣列聲波測井?dāng)?shù)據(jù)壓縮算法

林　博， 葉朝輝

(清華大學(xué) 自動化系，北京 100084)

?

基于DCT的陣列聲波測井?dāng)?shù)據(jù)壓縮算法

林博， 葉朝輝

(清華大學(xué) 自動化系，北京 100084)

摘要:陣列聲波測井是目前一種重要的測井方法。但由于它產(chǎn)生的數(shù)據(jù)量較大，加之電纜傳輸速度的限制，測井效率會因此受到影響。為提升測井效率，有必要對原始聲波測井?dāng)?shù)據(jù)進行壓縮處理。通過分析聲波信號的特點，提出了采用DCT+適當(dāng)量化+算術(shù)編碼的壓縮方法，并根據(jù)實際波形的測試結(jié)果對壓縮方法進行了一系列改進。通過對660組真實測井?dāng)?shù)據(jù)的測試，得到了對于672個點的單極聲波數(shù)據(jù)平均4.56的壓縮比，對于1 200個點的偶極聲波數(shù)據(jù)平均9.18的壓縮比，均方根誤差在2 %左右，證實了該方法的有效性。

關(guān)鍵詞：離散余弦變換; 陣列聲波測井; 數(shù)據(jù)壓縮

0引言

聲波測井是一種重要的地球物理測井方法，它主要利用不同介質(zhì)中聲波的速度、衰減等特性的不同來研究井下的地質(zhì)結(jié)構(gòu)。

聲波測井儀器作為一個嵌入式系統(tǒng)被安置在井內(nèi)，通過電纜將測井?dāng)?shù)據(jù)上傳到上位機進行數(shù)據(jù)解析。但是，由于聲波測井需要利用陣列的大量重復(fù)的信息來彌補可能出現(xiàn)的測量誤差[1]，因此，數(shù)據(jù)量較大；而傳輸線纜能夠提供的數(shù)據(jù)傳輸速度又是有限的[2]，這樣就會極大地影響測井效率。因此，有必要對于井下數(shù)據(jù)進行壓縮處理，傳輸壓縮后的數(shù)據(jù)，等數(shù)據(jù)到達地面上位機系統(tǒng)后再進行解壓縮處理。這樣不僅可以節(jié)省傳輸時間，提高測井效率，同時也能減小對于井下儀器存儲空間的需求。

聲波測井信號屬于一維非平穩(wěn)信號，對于這樣一類信號的壓縮處理目前主要采用的思路是變換+量化+熵編碼的方式。

文獻[3]使用了離散余弦變換(discrete cosine transform,DCT) +差分編碼+Huffman編碼的方法，利用了心電圖(ECG)信號中的周期性特點，在百分比均方根誤差 (percent root-mean-square difference,PRD)等于0.61 %的情況下達到了壓縮比 (compress ratio,CR)等于16.5的壓縮率。文獻[4]中作者使用了離散小波變換(discrete wavelet transform,DWT)+神經(jīng)元網(wǎng)絡(luò)量化+算術(shù)編碼對于肌電圖數(shù)據(jù)進行了壓縮研究，在PRD等于1.4 %時能夠達到CR等于5左右。文獻[5]使用了DWT+蘭博—立夫—衛(wèi)曲編碼(Lempel—Ziv—Welch encoding,LZW)的壓縮算法，用于橋梁振動實時檢測數(shù)據(jù)的壓縮處理，在PRD等于0.41 %時能夠達到CR等于10.86的壓縮率。由于不同種類的信號特性不同，最終得到的CR和PRD也是不盡相同的。

而目前對于聲波測井信號的壓縮方法研究并不多，文獻[6]提出了一種改進的SPIHT方法；文獻[7]嘗試使用LZW編碼對于聲波信號進行壓縮研究。

本文以DCT +適當(dāng)量化+熵編碼為主要思路對于660組聲波測井?dāng)?shù)據(jù)進行數(shù)據(jù)壓縮測試,并依據(jù)波形的特點提出改進壓縮效果的方法。

1壓縮方案設(shè)計

1.1聲波信號及其特點

聲波測井信號分為單極信號和偶極信號，分別由單極子聲源和偶極子聲源產(chǎn)生。

如圖1和圖2所示分別是單極陣列聲波測井?dāng)?shù)據(jù)和偶極陣列聲波測井?dāng)?shù)據(jù)的波形圖。單極聲波數(shù)據(jù)共有8個通道，從上到下依次為1～8通道，每個通道672個數(shù)據(jù)點，每個數(shù)據(jù)點是從A/D轉(zhuǎn)換得來的16 bit的數(shù)據(jù)。偶極聲波數(shù)據(jù)同樣有8個通道，每個通道有1 200個16 bit的數(shù)據(jù)點。(圖3是偶極數(shù)據(jù)的前300個數(shù)據(jù)點)。

圖1　單極信號波形圖Fig 1　Mono pole signal waveform

圖2　偶極信號波形圖Fig 2　Dipole signal waveform

圖3　偶極信號的前300個點Fig 3　First 300pts of dipole signal waveform

每個通道的單極信號依次分為縱波、橫波和斯通利波。而偶極信號除了縱波、橫波和撓曲波以外，還有很長的一段信號非常弱的遠探測波段。

聲波信號有這樣一些特點：信號隨時間變化而變化，沒有什么周期性。頻帶較豐富，有一定噪聲。8通道數(shù)據(jù)之間有一定的相關(guān)性，但又不完全相同。

信號有一定的噪聲，也就決定了聲波信號的數(shù)據(jù)并不十分適合使用通用無損壓縮的技術(shù)進行處理。因為噪聲信號帶來的隨機性會使得無損壓縮效果變壞；同時由于有了噪聲的干擾，也就意味著即使無誤差還原原始信號，意義也并不是很大。所以,壓縮方案設(shè)計首先對于聲波信號做變換編碼處理。

1.2DCT

變換編碼的基本思想是：通過映射變換，將原始數(shù)據(jù)投射到一個更為緊湊的表示空間去。在變換域空間重新表示原始數(shù)據(jù)，可以得到一系列更為有效的系數(shù)[8]。從某種意義上說，這些系數(shù)集中了更多的能量，減少了不確定性，因而，對于這些系數(shù)進行編碼需要的總比特數(shù)會更少。

這其中，在均方誤差度量準(zhǔn)則下，失真最小的一種變換方式是K—L變換，但因為該變換包含求解特征向量這一較為困難的操作，致使它并不具有實用性。1974年，Ahmed N等人提出了DCT成為替代K—L變換的準(zhǔn)最佳變換。DCT可以得到近似的最優(yōu)解，同時它又有快速的算法，這使得它在語音和圖像處理中都得到了廣泛應(yīng)用[8]。

DCT矩陣A具有以下形式

Ai,j=

通過Y=AX即可求得DCT后的系數(shù)矩陣Y。而由于A是一個正交矩陣，因此,反離散余弦變換(IDCT)的變換矩陣

A-1=AT

基于DCT的良好特性，使用它作為變換編碼方式是非常合適的。

1.3通用無損壓縮方法

變換后的系數(shù)具有較大的冗余度，可以使用通用無損壓縮方法進一步壓縮。這其中，LZW字典編碼和熵編碼較為常用。

但是，LZW變換由于是通過生成字典的方式用較短的字符替代較長的字符，因此,它對于字符的相對位置非常敏感。而聲波數(shù)據(jù)總體來說隨機性較強，即使經(jīng)過DCT后也較難保證有良好的長字符重復(fù)性，所以,LZW變換并不適合。

而熵編碼則不涉及到字符的相對位置，它只是利用了字符概率這一信息來減小平均碼長。可以證明，當(dāng)概率為pi的字符碼長Li等于它的自信息量-lgpi時，平均碼長l最小，最小值

(3)

即字符概率分布的信息熵。

算術(shù)編碼屬于熵編碼的一種，它采用化整為零的思想，將各符號的概率值映射為區(qū)間的大小，通過區(qū)間大小與碼字長度的一個對應(yīng)關(guān)系(區(qū)間越大，表示該區(qū)間的碼字長度越短)來實現(xiàn)碼字長度對自信息量的逼近。算術(shù)編碼的效果非常好，可以很好地逼近信息熵。因此，本文采用算術(shù)編碼作為通用無損編碼方法。

算術(shù)編碼的整個過程如圖4所示，假設(shè)要對于00,01,10,11這4個字符進行編碼，首先在[0,1]區(qū)間上根據(jù)概率分布給這4個符號劃定相應(yīng)大小的區(qū)間，然后找到第一個字符10的區(qū)間是[0.5,0.7]。再將這個區(qū)間按照符號的概率分布劃為4個區(qū)間，找第2個字符對應(yīng)的區(qū)間位置，依次類推，最終將得到一個很小的區(qū)間，在這個區(qū)間中隨意找到一個最短的2進制輸出，就是最終的編碼結(jié)果。

圖4　算術(shù)編碼過程示意圖 (圖片來源文獻[9])Fig 4　Process of arithmetic coding

1.4字符統(tǒng)計方法

波形數(shù)據(jù)是16bit的有符號整型，而使用算術(shù)編碼需要定義符號，統(tǒng)計它們的概率，則以多少比特為一個整體統(tǒng)計符號的概率才能達到最好的效果就成為一個值得研究的問題。

為了研究這一問題，不妨比較一下以8bit為一個整體進行字符統(tǒng)計的結(jié)果和以4bit為一個整體分別對高低字節(jié)進行字符統(tǒng)計的結(jié)果。

假設(shè)對于4bit的高字節(jié)的統(tǒng)計概率是[p0,p1,…,p15]，對于其低字節(jié)的統(tǒng)計概率是[q0,q1,…,q15]。而8bit的統(tǒng)計概率是[s0,s1,…,s255]。

以4bit為單位的信息熵為

(4)

以8 bit為單位的信息熵為

(5)

因此，H(4)≥H(8)，且只有當(dāng)si,j=piqj時，即高低位統(tǒng)計概率相互獨立時，等號才成立。

基于以上分析，可以得到：H(16)≥H(8)≥H(4)≥H(2)≥H(1),即一般情況下,以更多的bit為單位進行字符概率統(tǒng)計能達到更好的壓縮效果。但是，由于還要考慮樣本數(shù)量的影響，如果位數(shù)過多造成樣本過少，也會影響信息熵的逼近，因此,這里選擇8 bit為統(tǒng)計單位。

2壓縮效果與改進

2.1初始壓縮效果

如上所述，采用DCT編碼，加上以8 bit為單位的算術(shù)編碼進行壓縮測試。由于低8位的數(shù)據(jù)沒有明顯的統(tǒng)計規(guī)律，因此，考慮舍去。這樣做會帶來一定的誤差，但由于低8位所占數(shù)值較小，所占有的能量不大，因此，誤差并不大。

選取5個深度點的單極聲波數(shù)據(jù)為樣本，CR和PRD的測試結(jié)果如表1所示。

表1　DCT+算術(shù)編碼+舍去DCT系數(shù)低8位后的

可以看到:CR約為4左右，而PRD約為2 %，換算成信噪比(SNR)約為33 dB，這樣的信號誤差是完全可以接受的。

2.2對DCT編碼的分類處理

DCT后的系數(shù)(16 bit)呈現(xiàn)低頻部分系數(shù)大，高頻部分系數(shù)小的統(tǒng)計規(guī)律，如圖5所示。

可以利用這樣的規(guī)律對于DCT系數(shù)進行分割統(tǒng)計，前300個系數(shù)數(shù)值較大，占有的能量較大，但其統(tǒng)計規(guī)律不如后372個明顯。將二者分開可以使后372個系數(shù)達到更好的概率分布。

圖5　DCT系數(shù)的波形Fig 5　Waveform of DCT coefficient

事實上，672個采樣點的時間間隔是8 μs，也就意味著最高頻率為125 kHz，做完DCT后，DCT系數(shù)的頻率間隔為186 Hz。而單極聲波信號的最高頻率只有23 kHz，因此，只有前125個系數(shù)擁有較為主要的能量。

同樣，對于5個樣本進行測試，分割統(tǒng)計后得到的CR如表2所示，可以看到整體CR有了一定的提升。

表2　DCT系數(shù)分開處理的CR

2.38通道間的數(shù)據(jù)關(guān)聯(lián)

陣列聲波的數(shù)據(jù)共有8個通道，這8個通道間的數(shù)據(jù)是不同的聲波接收器采集的，波形數(shù)據(jù)不盡相同，但是這些數(shù)據(jù)又來源于同一個采樣深度點，因此它們又具有一定的聯(lián)系。

可以考慮利用通道間數(shù)據(jù)的關(guān)聯(lián)進一步提升壓縮效果。具體做法是將8通道的數(shù)據(jù)整合到一起做DCT，這樣由于數(shù)據(jù)之間的關(guān)聯(lián)會使得DCT系數(shù)呈現(xiàn)出更為集中的特性，從而能達到更好的壓縮效果。

5個樣本的測試結(jié)果如表3所示，可以看到CR又有了一定程度的提高。

表3　合并8通道波形后做DCT得到的CR

3最終的結(jié)果

3.1壓縮算法流程

最終的壓縮算法流程如圖6所示，其中虛線框出的是應(yīng)該在上位機上實現(xiàn)的解壓縮步驟。

圖6　壓縮算法流程框圖Fig 6　Flow chart of compression algorithm

3.2660組單偶極數(shù)據(jù)測試

使用660組單偶極數(shù)據(jù)進行測試，測試結(jié)果的平均值如表4所示。

4結(jié)束語

本文著重研究了基于DCT的陣列聲波測井?dāng)?shù)據(jù)壓縮方法，并根據(jù)聲波數(shù)據(jù)自身的特點改進了CR。通過大量的數(shù)據(jù)測試得出了平均的CR和PRD，進而從理論和實踐證實:該方法是行之有效的。

表4　660組數(shù)據(jù)測試結(jié)果

參考文獻:

[1]唐曉明,鄭傳漢,趙曉敏.定量測井聲學(xué)[M].北京：石油工業(yè)出版社,2004.

[2]劉雁飛.基于CAN總線的聲波測井地面數(shù)據(jù)接口與傳輸系統(tǒng)設(shè)計[D].青島：中國石油大學(xué)(青島),2009.

[3]Lee S J,Kim J,Lee M.A real-time ECG data compression and transmission algorithm for an e-health device[J].IEEE Transactions on Biomedical Engineering,2011,58(9):2448-2455.

[4]de A Berger P,de O Nascimento F A,da Rocha A F,et al.A new wavelet-based algorithm for compression of EMG signals[C]∥2007 The 29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society,EMBS 2007,2007:1554-1557.

[5]Zhang Fengyuan,Yang Dong,Gong Xuyan,et al.A wavelet-based data compression algorithm for bridge vibration[C]∥2012 IEEE 14th International Conference on Communication Technology(ICCT),IEEE,2012:334-341.

[6]張偉,師奕兵.聲波測井?dāng)?shù)據(jù)壓縮的一種 SPIHT改進算法[J].電子測量與儀器學(xué)報,2008,22(1):15-19.

[7]葉平凡.聲波測井?dāng)?shù)據(jù)壓縮算法研究及實現(xiàn)[D].成都：電子科技大學(xué),2012.

[8]吳樂南.數(shù)據(jù)壓縮[M].北京：電子工業(yè)出版社,2000：108，112-113.

[9]林福宗.無損數(shù)據(jù)壓縮[DB/OL].[2015—04—01].http:∥source.mastvu.ah.cn/dzjxb/dmtzzsyrj-jm/mxzy/MultiMedia/multimedia/course/course/course1—4—2.html.

Data compression algorithm for DCT-based array acoustic logging

LIN Bo, YE Zhao-hui

(Department of Automation,Tsinghua University,Beijing 100084,China)

Abstract:Array acoustic logging is one of the most important logging methods presently.However,due to its huge amounts of data and the speed limit of data transmission,logging efficiency is restricted.Data compression and processing is necessary to improve logging efficiency.By analyzing characters of acoustic signal,introduce a compression method comprised of DCT,proper quantization and arithmetic coding,further,a series of improvements are made to the method.By test on 660 group real logging data,it is proved that the method is effective,average compression ratio(CR)is 4.56 to the monopole acoustic data of 672 points and 9.18 to dipole acoustic data to 1200 points,with root-mean-square error around 2 %.

Key words:discrete cosine transform(DCT); array acoustic logging; data compression

DOI:10.13873/J.1000—9787(2016)02—0143—04

收稿日期：2015—04—28

中圖分類號：TP 274

文獻標(biāo)識碼：A

文章編號：1000—9787(2016)02—0143—04

作者簡介:

林博(1990-)，男，安徽淮南人，工學(xué)碩士，主要從事檢測與電子技術(shù)方向的研究。