李可寒 ,楊俊峰 ,宋克柱
(1.中國科學技術(shù)大學核探測與核電子學國家重點實驗室,合肥,230026;2.中國科學技術(shù)大學近代物理系,合肥,230026)
地震勘探是當前油氣勘探中一種應(yīng)用廣泛的方法[1],它對數(shù)據(jù)采集系統(tǒng)的要求有多通道、高精度、數(shù)據(jù)量大和實時處理等[2]。數(shù)據(jù)采集系統(tǒng)在采集到地震數(shù)據(jù)之后,將采集的數(shù)據(jù)匯聚至控制與紀錄模塊供以后數(shù)據(jù)分析使用?;诘卣鹂碧綌?shù)據(jù)采集系統(tǒng)的實時、多通道等要求,數(shù)據(jù)采集系統(tǒng)需要實時采集和傳輸大量數(shù)據(jù)。
在已有的地震數(shù)據(jù)采集系統(tǒng)中,通常使用拖纜中長度為100 m 的傳輸線進行數(shù)據(jù)傳輸。由于傳輸距離較長,如果傳輸速率過快將導致誤碼率增高。因此,在傳輸線質(zhì)量無法提升的情況下,可以使用數(shù)據(jù)壓縮來降低需要的傳輸速率,從而降低數(shù)據(jù)傳輸?shù)恼`碼率,同時也能降低系統(tǒng)的功耗。
現(xiàn)有的地震數(shù)據(jù)壓縮算法[3-4]利用地震波的物理特征采用小波變換、Dreamlet 變換等方式能夠起到很好的壓縮效果。但這些算法的計算過程較復雜,而且需要輸入整個地震數(shù)據(jù)文件才能完成壓縮,很難用于基于現(xiàn)場可編程門陣列(Field programmable gate array,FPGA)的實時采集處理系統(tǒng),因此需要一種針對地震數(shù)據(jù)流的壓縮算法。
現(xiàn)有的地震數(shù)據(jù)流壓縮算法[5]是一種針對數(shù)據(jù)流的無損壓縮算法,根據(jù)地震數(shù)據(jù)的數(shù)值將3 字節(jié)的數(shù)據(jù)壓縮至1~4 字節(jié),是一種熵編碼。
本文考慮地震波的物理特征,同時使用熵編碼,實現(xiàn)了一種針對地震數(shù)據(jù)流的無損壓縮算法。
海上地震數(shù)據(jù)采集系統(tǒng)[6-7]分為水下部分和水上部分。水下部分包括采集和傳輸數(shù)據(jù)的拖纜,系統(tǒng)可以有多條拖纜;水上部分包括數(shù)據(jù)處理機箱和室內(nèi)控制顯示工作站,數(shù)據(jù)處理機箱負責轉(zhuǎn)發(fā)控制命令和匯總拖纜數(shù)據(jù)。由于水下部分的傳輸距離更長,因此數(shù)據(jù)傳輸?shù)钠款i主要在水下部分。
地震數(shù)據(jù)采集系統(tǒng)支持若干條拖纜,其中單條水下拖纜的結(jié)構(gòu)如圖1 所示。每段拖纜長100 m,里面包含16 個傳感器,道間距為6.25 m,采用24 位的ADC,最大采樣率為2 ksps,最大要求支持150 段拖纜。兩段拖纜之間由傳輸節(jié)點相連,每一個傳輸節(jié)點除了負責收集本級采集數(shù)據(jù)和狀態(tài)數(shù)據(jù)外,還負責接收后級傳輸節(jié)點數(shù)據(jù)并轉(zhuǎn)發(fā)。根據(jù)系統(tǒng)結(jié)構(gòu),越往前級數(shù)據(jù)量越大,其中滿負荷情況下首節(jié)點的凈數(shù)據(jù)量可以達到150×16×24×2 000=115.2 (Mb/s)。如果考慮協(xié)議開銷,首節(jié)點需要的傳輸速率則會更大。
圖1 水下拖纜整體結(jié)構(gòu)Fig.1 Integral structure of underwater towline
當拖纜段數(shù)增加時,首節(jié)點需要更大的傳輸速率。為了控制誤碼率,可以考慮對地震數(shù)據(jù)進行壓縮,減少需要傳輸?shù)臄?shù)據(jù)量。另外,傳輸速率的減小能夠降低單個節(jié)點的功耗,從而使系統(tǒng)能夠支持更多段拖纜。
現(xiàn)有的大部分地震數(shù)據(jù)壓縮算法針對的是一個地震數(shù)據(jù)文件,用于減小地震數(shù)據(jù)占用的存儲空間,無法在數(shù)據(jù)采集的時候使用以提高傳輸效率?,F(xiàn)有的地震數(shù)據(jù)流壓縮算法[5]針對的是單個地震數(shù)據(jù),利用地震數(shù)據(jù)的分布特征,將3 字節(jié)的地震數(shù)據(jù)壓縮至1~4 字節(jié)。該算法沒有考慮數(shù)據(jù)之間的相關(guān)性,故而壓縮性能有待提升。
本文將根據(jù)地震數(shù)據(jù)的時空特征,同時借鑒語音壓縮算法,提出一種易于用FPGA 實現(xiàn)的24 位地震數(shù)據(jù)流無損壓縮算法。
任意兩段拖纜之間由一個傳輸節(jié)點相連,每個傳輸節(jié)點負責收集附近16 個傳感器采集到的地震數(shù)據(jù)。為了提高壓縮效果,傳輸節(jié)點通常需要緩存若干個采樣時刻的數(shù)據(jù)。圖2 為隨機截取的連續(xù)16 個通道、16 個采樣時刻的地震數(shù)據(jù),即為待壓縮的數(shù)據(jù)。
從圖2 可以看出,待壓縮的16×16 地震數(shù)據(jù)具有如下特點:(1)傳感器的數(shù)值隨采樣時刻的變化不大;(2)傳感器的數(shù)值隨傳感器號的變化不可預知。
圖2 地震數(shù)據(jù)Fig.2 Seismic data
傳感器收集的是地震波數(shù)據(jù),地震波本質(zhì)上是機械波,理想情況下應(yīng)該是一個類正弦波信號。激發(fā)人工地震波的震源頻率約為幾十赫茲,而傳感器的采樣周期通常為毫秒量級(對應(yīng)頻率為千赫茲量級),即采樣頻率遠大于振動頻率。因此,在不長的時間里,某一個傳感器的數(shù)值不會有太大的變化。而在同一時刻,相鄰傳感器采集的數(shù)據(jù)可能相近,也可能相差很大,這取決于地層結(jié)構(gòu)。由于地層結(jié)構(gòu)無法預知,所以很難事先得知同時刻不同傳感器數(shù)據(jù)之間的關(guān)系。
因此,本文中的壓縮主要考慮針對同一個傳感器在若干個采樣時刻的地震數(shù)據(jù),這樣的數(shù)據(jù)類似于機械波(聲波),故可以借鑒語音壓縮的相關(guān)算法。
對地震數(shù)據(jù)的壓縮可以借鑒語音壓縮,一種常見的語音壓縮算法分為預測編碼和熵編碼兩步,其過程如圖3 所示。
圖3 語音壓縮過程Fig.3 Process of speech compression
預測編碼是根據(jù)離散信號之間存在著一定關(guān)聯(lián)性的特點,利用前一個或多個信號預測下一個信號,然后對實際值和預測值的差(預測誤差)進行編碼。如果預測比較準確,誤差就會很小。在同等精度要求的條件下,就可以用比較少的比特進行編碼,達到壓縮數(shù)據(jù)的目的。
常見的預測編碼有線性預測編碼,即現(xiàn)在值可用若干個過去值的加權(quán)線性組合來逼近,加權(quán)系數(shù)即為預測系數(shù),預測值為
為使線性預測誤差
最小,可以根據(jù)實際采樣值讓線性預測誤差的平方達到最小值,從而確定最佳預測系數(shù)。
由于線性預測編碼比較復雜,尤其是預測系數(shù)的計算在硬件上實現(xiàn)起來難度很大,因此考慮一種更加簡單的線性預測編碼,即差分編碼。
對于原始數(shù)據(jù)X1,X2,…,Xn,一階差分編碼為X1,X2-X1,…,Xn-Xn-1,一階差分編碼是預測階數(shù)p=1,預測系數(shù)A1=1 的線性預測編碼。預測值=Xn-1,即估計后一個數(shù)據(jù)和前一個數(shù)據(jù)差不多。
二階差分編碼為X1,X2-X1,X3-2X2+X1,…,Xn-2Xn-1+Xn-2,二階差分編碼是預測階數(shù)p=2,預測系數(shù)A1=2,A2=-1 的線性預測編碼。預測值=2Xn-1-Xn-2,即用前兩個點擬合一次函數(shù),并預測下一個點會出現(xiàn)在該直線上。
以此類推,可得高階差分編碼。
差分編碼只需要作若干次減法,易于硬件上的實現(xiàn)。
由于地震數(shù)據(jù)及其各階導數(shù)隨時間的變化不是很明顯,因此地震數(shù)據(jù)在差分編碼之后,其絕對值將變得很?。ㄏ啾仍瓉淼臄?shù)據(jù))?,F(xiàn)在就要想辦法用較少的比特來存儲這些絕對值很小的數(shù)據(jù)。
熵編碼又稱為統(tǒng)計編碼,是根據(jù)消息出現(xiàn)概率的分布特性而進行的無損數(shù)據(jù)壓縮編碼(差分編碼之后,絕對值小的數(shù)據(jù)出現(xiàn)的概率更大)。
2.4.1 指數(shù)哥倫布編碼[8]
指數(shù)哥倫布編碼刪去了數(shù)據(jù)高位多余的0(負數(shù)刪去多余的1),并添加上若干位的標記位。表示非負整數(shù)N的k階指數(shù)哥倫布編碼生成方法如下:
(1)將N用二進制碼表示,去掉低位的k個比特,然后加1;
(2)計算留下的比特數(shù),將此數(shù)減1,并記作m;
(3)將第(1)步中去掉的k個比特補回串尾,并在串頭添加m個0。
例如,對1 階指數(shù)哥倫布編碼,0 到13 的編碼如下:
0——0000——10
1——0001——11
2——0010——0100
3——0011——0101
4——0100——0110
5——0101——0111
6——0110——001000
7——0111——001001
8——1000——001010
9——1001——001011
10——1010——001100
11——1011——001101
12——1100——001110
13——1101——001111
由此可見,如果數(shù)據(jù)出現(xiàn)小幅度的概率比出現(xiàn)大幅度的概率大,則通過指數(shù)哥倫布編碼可以較好地壓縮數(shù)據(jù)。但如果數(shù)據(jù)不滿足上述特征,則壓縮后的數(shù)據(jù)可能會比原來更大。
k階指數(shù)哥倫布編碼包含m位的0 + 1 位的1 + (m+k)位的數(shù)據(jù),若要解碼,對于某串k階指數(shù)哥倫布編碼的數(shù)據(jù),先計算數(shù)據(jù)頭部0 的個數(shù),并記作m;則除去這m個0,該數(shù)據(jù)的有效位數(shù)為m+k+1,將該數(shù)據(jù)減去2k,即可得到指數(shù)哥倫布編碼數(shù)據(jù)的解碼,即原始數(shù)據(jù)。
由編碼和解碼的過程可知,解碼能過完全恢復出編碼之前的數(shù)據(jù),不會損失數(shù)據(jù)的精度,因此指數(shù)哥倫布編碼能夠?qū)崿F(xiàn)無損壓縮。
單個數(shù)據(jù)的指數(shù)哥倫布編碼很容易用硬件實現(xiàn),然而要壓縮的數(shù)據(jù)往往有很多個,在指數(shù)哥倫布編碼之后通常要將這些數(shù)據(jù)拼接起來。由于這些數(shù)據(jù)往往不是整數(shù)字節(jié),而且長度不一,這會導致數(shù)據(jù)的拼接需要消耗大量的邏輯資源。因此,指數(shù)哥倫布編碼并不是熵編碼的最佳選擇。
2.4.2 整數(shù)字節(jié)熵編碼
指數(shù)哥倫布編碼之后的數(shù)據(jù)不是整數(shù)字節(jié),而且長度也不確定,這給后續(xù)處理帶來了很大的麻煩,下面考慮一種將差分編碼后的數(shù)據(jù)壓縮成整數(shù)個字節(jié)的熵編碼。
對于24 bit 地震數(shù)據(jù),一階差分后的數(shù)據(jù)可用25 bit 表示,其大小位于[-16 777 216 , 16 777 215]區(qū)間,并且差分數(shù)據(jù)中的大多數(shù)數(shù)據(jù)數(shù)值都很小,因此可按照如下方式[5]壓縮數(shù)據(jù):
(1)如果差分數(shù)據(jù)位于[-64 , 63]區(qū)間,可用1 個字節(jié)A[7:0]表示該數(shù)據(jù)。A[6]為差分數(shù)據(jù)的符號位,A[7]=0,A[6:0]為差分數(shù)據(jù)的有效值。
(2)如果差分數(shù)據(jù)位于 [-8 192 , 8 191]區(qū)間,可用2 個字節(jié)B[15:0]表示該數(shù)據(jù)。B[13]為差分數(shù)據(jù)的符號位,B[15:14]=10,B[13:0]為差分數(shù)據(jù)的有效值。
(3)如果差分數(shù)據(jù)位于[-1 048 576 , 1 048 575]區(qū)間,可用3 個字節(jié)C[23:0]表示該數(shù)據(jù)。C[20]為差分數(shù)據(jù)的符號位,C[23:21]=110,C[20:0]為差分數(shù)據(jù)的有效值。
(4)無論如何,可用 4 個字節(jié) D[31:0]表示該數(shù)據(jù)。D[24]為差分數(shù)據(jù)的符號位,D[31:25]=1 110 000,D[24:0]為差分數(shù)據(jù)的絕對值。
以上方法簡單易實現(xiàn),可將25bit 差分數(shù)據(jù)壓縮為1~4 字節(jié)的數(shù)據(jù)。
由于壓縮后的數(shù)據(jù)由若干個字節(jié)組成,因此如要解壓,首先要判斷哪幾個字節(jié)(1 到4 字節(jié))對應(yīng)一個差分數(shù)據(jù)。首先從緩存區(qū)中讀取一個字節(jié)的數(shù)據(jù)A[7:0],如果A[6]=0,該數(shù)據(jù)為1 個字節(jié);如果A[6:5]=10,該數(shù)據(jù)為 2 個字節(jié);如果 A[6:4]=110,該數(shù)據(jù)為 3 個字節(jié);如果 A[6:4]=111,該數(shù)據(jù)為 4 個字節(jié)。
先通過以上方法將差分后的數(shù)據(jù)分離出來,轉(zhuǎn)換成25 bit 的數(shù)據(jù);再進行差分逆變換,得到24 bit 的原始地震數(shù)據(jù)。
由編碼和解碼的過程可知,解碼能過完全恢復出編碼之前的數(shù)據(jù),不會損失數(shù)據(jù)的精度,因此整數(shù)字節(jié)熵編碼能夠?qū)崿F(xiàn)無損壓縮。
如果地震數(shù)據(jù)不是24 bit,或者采用的差分編碼不是一階差分編碼,也可以使用類似的整數(shù)字節(jié)熵編碼。
為了減少傳輸節(jié)點需要上傳的數(shù)據(jù)量,可以利用傳輸節(jié)點中的FPGA 對采集到的24 位地震數(shù)據(jù)進行壓縮;當壓縮后的數(shù)據(jù)到達船上時,再利用數(shù)據(jù)處理機箱中的FPGA 對數(shù)據(jù)解壓。每個傳輸節(jié)點只負責壓縮附近16 個傳感器的地震數(shù)據(jù)。
因為對FPGA 的功耗有一定的要求,所以很難使用比較復雜的壓縮算法。本文實現(xiàn)了一種用一階差分編碼和整數(shù)字節(jié)熵編碼來完成地震數(shù)據(jù)壓縮的算法。
使用一階差分編碼和整數(shù)字節(jié)熵編碼來完成地震數(shù)據(jù)的壓縮。在每一個采樣周期內(nèi),需要輸入24 bit 數(shù)據(jù),然后輸出1~4 字節(jié)的數(shù)據(jù)。數(shù)據(jù)壓縮模塊的輸入輸出端口如圖4 所示。
數(shù)據(jù)的輸入和輸出均使用AXI-Stream 接口。輸入的數(shù)據(jù)是24 bit 的地震數(shù)據(jù),i_last 表示某一次壓縮的最后一個數(shù)據(jù);輸出的數(shù)據(jù)是壓縮后的數(shù)據(jù),壓縮后的數(shù)據(jù)以8 bit 為單位,o_last 表示壓縮后數(shù)據(jù)的最后一個字節(jié)。數(shù)據(jù)壓縮過程分差分編碼和熵編碼兩步來實現(xiàn)。
3.1.1 差分編碼
差分編碼由一個24 位寄存器和一個24 位減法器實現(xiàn),寄存器保存上一個數(shù)據(jù),減法器利用當前輸入的數(shù)據(jù)減去上一個數(shù)據(jù)得到差分數(shù)據(jù),如圖5 所示。圖中:i_last 控制寄存器的清零(當最后一個數(shù)據(jù)到來時);i_valid 和i_ready 控制寄存器的使能端(僅i_valid 和i_ready 都有效時才進行差分編碼)。
圖4 數(shù)據(jù)壓縮模塊Fig.4 Module of data compression
圖5 差分編碼器Fig.5 Differential encoder
3.1.2 熵編碼
整數(shù)字節(jié)熵編碼可用一個組合邏輯電路來實現(xiàn),輸入25 位數(shù)據(jù)data_dif 和有效的字節(jié)數(shù)byte_cnt,輸 出 4 字 節(jié) 數(shù) 據(jù) data_out(1~4 字 節(jié) 有效),Verilog 代碼如圖 6 所示。
有效的字節(jié)數(shù)byte_cnt 可根據(jù)輸入數(shù)據(jù)data_dif 高位多余的 0(負數(shù)為 1)得到。
在熵編碼之后,需要將熵編碼器輸出的數(shù)據(jù)data_out[31:0]從 o_data[7:0]端口輸出,根據(jù) byte_cnt 的不同需要 1~4 個時鐘周期。在輸出完成之前,將i_ready 拉低,暫時停止前級的輸入,直到輸出結(jié)束才將i_ready 拉高。
圖6 整數(shù)字節(jié)熵編碼器的Verilog 實現(xiàn)Fig.6 Verilog implementation of integer byte entropy encoder
數(shù)據(jù)解壓在船上數(shù)據(jù)處理機箱中的FPGA 上進行。數(shù)據(jù)解壓模塊也是使用AXI-Stream 接口,輸入8 bit 壓縮數(shù)據(jù),輸出24 bit 地震數(shù)據(jù)。數(shù)據(jù)解壓過程分反熵編碼和反差分編碼兩步。
3.2.1 反熵編碼
反熵編碼的第1 步是根據(jù)每個壓縮數(shù)據(jù)的第一個字節(jié)來判斷該數(shù)據(jù)占用了幾個字節(jié),然后將這幾個字節(jié)的數(shù)據(jù)當作反熵編碼器的輸入,并開始解析下一個數(shù)據(jù),該過程可以通過一個狀態(tài)機來實現(xiàn)。
反熵編碼器則可以通過一個簡單的組合邏輯電路實現(xiàn),輸入4 字節(jié)數(shù)據(jù)data_in(1~4 字節(jié)有效)和有效的字節(jié)數(shù)byte_cnt,輸出25 bit 數(shù)據(jù)data_dif,Verilog 代碼如圖 7 所示。
有效的字節(jié)數(shù)byte_cnt 通過反熵編碼第一步的狀態(tài)機得到。
3.2.2 反差分編碼
類似于差分編碼,反差分編碼可用一個25 位寄存器和一個25 位加法器實現(xiàn),寄存器保存之前所有數(shù)據(jù)的和值(即上一個原始數(shù)據(jù)),加法器將當前差分數(shù)據(jù)和上一個原始數(shù)據(jù)相加得到當前原始數(shù)據(jù),如圖8 所示。圖中:i_last 控制寄存器的清零(當最后一個數(shù)據(jù)到來時);i_valid 控制寄存器的使能端(僅i_valid 有效時才進行反差分編碼)。
圖7 反熵編碼器的Verilog 實現(xiàn)Fig.7 Verilog implementation of inverse entropy encoder
圖8 反差分編碼器Fig.8 Inverse differential encoder
壓縮率指數(shù)據(jù)壓縮之后的大小與數(shù)據(jù)壓縮之前的大小之比。利用本文中所述的一階差分編碼和整數(shù)字節(jié)熵編碼,表1 給出了3 組實測地震數(shù)據(jù)的壓縮結(jié)果。
每一組地震數(shù)據(jù)均由2 336 個傳感器進行12 288 次采樣得到。每次壓縮1 個傳感器64 次采樣得到的地震數(shù)據(jù),壓縮前數(shù)據(jù)長度為64×3=192 字節(jié)??倝嚎s次數(shù)為2 336×12 288/64=448 512 次。剔除掉全零數(shù)據(jù)(采集系統(tǒng)未開始工作),需要壓縮的次數(shù)記作總幀數(shù)。使用本文中所述的壓縮算法壓縮這些數(shù)據(jù),記錄這些數(shù)據(jù)在壓縮之后的總長度(記作總字節(jié)數(shù))。壓縮后的平均長度=總字節(jié)數(shù)/總幀數(shù)。壓縮率=平均長度/192。3 組地震數(shù)據(jù)的平均壓縮率=(101+105+105)/3/192=54%。
數(shù)據(jù)壓縮在傳輸節(jié)點中的FPGA 上進行,相應(yīng)FPGA 的型號是Altera 公司的Cyclone V 5CEBA7F23C7,壓縮模塊占用的資源如圖9 所示。壓縮模塊需要的邏輯資源為54 個ALMs,不到總量的千分之一。
數(shù)據(jù)解壓在船上數(shù)據(jù)處理機箱中的FPGA 上進行,F(xiàn)PGA 的型號為Xilinx 公司的Kintex-7 xc7k325tffg900-2,解壓模塊占用的資源如表2 所示。解壓模塊需要的邏輯資源為79 個LUT,不到總量的萬分之五。
由此可見,不管是壓縮模塊,還是解壓模塊,占用FPGA 的資源都可以忽略不計。
表1 壓縮率測試結(jié)果Table 1 Test result of compression ratio
圖9 壓縮模塊占用資源Fig.9 Utilization of compression module
表2 解壓模塊占用資源Table 2 Utilization of decompression module
圖10 是壓縮模塊的仿真結(jié)果。在圖10 中,采集的24 位數(shù)據(jù)依次是{580,600,630,900,-900(}3 字節(jié)帶符號整數(shù)),壓縮后的數(shù)據(jù)依次是{82,44,14,1E,81,0E,E0,F(xiàn)F,F(xiàn)8,F(xiàn)8(}十六進制)??梢园l(fā)現(xiàn),相鄰兩個數(shù)據(jù)之間的差值越小,則數(shù)據(jù)壓縮之后的長度也會越小。同時,由于壓縮一個數(shù)據(jù)需要多個時鐘周期(1~4 個),因此輸入的數(shù)據(jù)之間應(yīng)該間隔若干個時鐘周期。
圖11 是解壓模塊的仿真結(jié)果,直接利用壓縮模塊的仿真結(jié)果當作解壓模塊的輸入,輸入數(shù)據(jù)依此為{82,44,14,1E,81,0E,E0,F(xiàn)F,F(xiàn)8,F(xiàn)8(}十六進制),輸出的數(shù)據(jù)是{580,600,630,900,-900(}3 字節(jié)帶符號整數(shù))??梢钥闯?,解壓得到的數(shù)據(jù)正是壓縮之前的數(shù)據(jù)。
圖10 壓縮模塊仿真結(jié)果Fig.10 Simulation result of compression module
圖11 解壓模塊仿真結(jié)果Fig.11 Simulation result of decompression module
海洋地震勘探方法是目前海洋油氣勘探的主要方法,更深層和更精確的海洋勘探需要更長的勘探拖纜,這將導致系統(tǒng)的功耗和傳輸速率大幅度增加。本文提出了一種可以提高傳輸效率的24位地震數(shù)據(jù)流無損壓縮算法,可以大幅度降低采集系統(tǒng)需要的傳輸速率,同時略微降低采集系統(tǒng)的功耗。地震數(shù)據(jù)的壓縮率平均可達54%,數(shù)據(jù)壓縮算法很容易用FPGA 實現(xiàn),占用FPGA 的邏輯資源極少,無須占用或只需占用較少的存儲資源,數(shù)據(jù)壓縮可在幾個時鐘周期之內(nèi)完成。簡而言之,本文所述的壓縮算法可以只花費極少的代價就將地震數(shù)據(jù)采集系統(tǒng)需要傳輸?shù)牡卣饠?shù)據(jù)量減少一半左右。