基于稀疏表示的陣列聲波測井儀數(shù)據(jù)無損壓縮傳輸方法

李 明,尹時(shí)松,張 寧,李波宏,莊獻(xiàn)華,李江山

(中國石油集團(tuán)測井有限公司華北分公司,河北 任丘 062552)

作為一種重要的成像測井儀器，交叉多極子陣列聲波測井儀[1]將單極發(fā)射、接收陣列和偶極發(fā)射、接收陣列交叉組合；其中，單極子陣列包括關(guān)于井軸對稱、全方位的2個(gè)發(fā)射換能器和8個(gè)接收換能器，偶極子陣列包含正交放置的2個(gè)發(fā)射換能器和16個(gè)接收換能器。各個(gè)接收換能器記錄聲源發(fā)射后一定時(shí)間段內(nèi)的聲波波列，通過時(shí)間-慢度相關(guān)方法[2]計(jì)算地層的縱波、橫波和斯通利波慢度，可用于識別巖性、計(jì)算巖石力學(xué)參數(shù)、識別氣層和裂縫，而由偶極陣列可得到地層快、慢橫波速度和方位，進(jìn)行地層各向異性分析，已廣泛應(yīng)用于碳酸鹽巖油氣層[3]、煤層氣[4]和致密砂巖氣層[5]的勘探開發(fā)。

陣列聲波測井儀聲波發(fā)射、接收換能器眾多，且采集聲波波列，因此儀器上傳地面的數(shù)據(jù)量很大。目前儀器僅有約42 kbit/s的最大電纜上傳速率，一旦提高儀器的測井速度，傳輸隊(duì)列中當(dāng)前地層的測量數(shù)據(jù)在未上傳到地面系統(tǒng)的情況下，就會(huì)被新地層的測量數(shù)據(jù)覆蓋。因此，在井下儀器采用無損數(shù)據(jù)壓縮上傳是在采集數(shù)據(jù)不失真的同時(shí)，增加傳輸帶寬、提升儀器測井時(shí)效的有效方法。

目前主流的測井?dāng)?shù)據(jù)無損壓縮方法是LZW(Lempel Ziv Welch,藍(lán)波-立夫-衛(wèi)曲)算法[6-9]，但LZW算法在編解碼過程中，需要存儲(chǔ)編碼字典，因此當(dāng)非平緩變化數(shù)據(jù)流的信息冗余較少時(shí)，壓縮率較低；文獻(xiàn)[10]針對陣列聲波測井儀提出一種預(yù)測編碼方法，但沒有考慮聲波波列的數(shù)值時(shí)序分布規(guī)律，預(yù)測誤差編碼碼位較長，導(dǎo)致預(yù)測器性能不佳；文獻(xiàn)[11]提出一種基于BWT(Burrows-Wheeler Transform,伯羅斯-惠勒變換)和PPM(Prediction with Partial Matching,局部匹配預(yù)測)算法的陣列聲波數(shù)據(jù)無損壓縮傳輸方法，雖然考慮了聲波波列中數(shù)據(jù)分布規(guī)律，但算法復(fù)雜程度較高，要求實(shí)際應(yīng)用中使用高性能的DSP硬件，且壓縮率最高只有51%，不能達(dá)到當(dāng)前測井作業(yè)效率指標(biāo)。

以上的壓縮方法都是基于傳統(tǒng)的信息熵減少理論，將聲波波列數(shù)據(jù)盡可能地轉(zhuǎn)換為平緩變化的數(shù)據(jù)流。如果從信號處理的角度考慮，將聲波波列從時(shí)間域轉(zhuǎn)換到某一變換域，且在該變換域，聲波波列只需極少的量值表征，那么井下上傳這些表征量值到地面系統(tǒng)，反變換后，就可恢復(fù)原始聲波波列。因此，從信號稀疏表示的角度出發(fā)，通過稀疏變換將聲波波列變換到其稀疏域。選擇恰當(dāng)?shù)南∈枳儞Q矩陣，可使稀疏表示系數(shù)中非零元素個(gè)數(shù)遠(yuǎn)小于波列長度。這里，采用K-SVD(K-means Singular Value Decomposition,K均值奇異值分解)算法，通過歷史聲波波列數(shù)據(jù)訓(xùn)練稀疏變換矩陣。數(shù)據(jù)上傳時(shí)，將非零稀疏表示系數(shù)和重構(gòu)波列與原始波列的誤差進(jìn)行壓縮編碼上傳。在地面系統(tǒng)采用相同的稀疏變換矩陣進(jìn)行反變換，再疊加誤差信息，恢復(fù)原始波列。通過HB油田陣列聲波仿真實(shí)驗(yàn)、實(shí)測數(shù)據(jù)壓縮率比對和實(shí)際作業(yè)的應(yīng)用測試，驗(yàn)證方法在提升數(shù)據(jù)壓縮率和提高測井作業(yè)時(shí)效方面的有效性。方法的設(shè)計(jì)和開發(fā)有望提升陣列聲波測井的施工速度，為油田勘探的提質(zhì)提效提供有力的技術(shù)支撐。

1 陣列聲波信號采集和傳輸原理

交叉多極子陣列聲波測井儀在一個(gè)深度點(diǎn)記錄12道單極子波列和32道偶極波列，包括8道單極子聲源T1發(fā)射，8個(gè)接收器采集的全波波列、4道單極子聲源T2發(fā)射，前4個(gè)接收器采集的全波波列和32道正交的兩個(gè)偶極子聲源交替發(fā)射，交叉放置的16個(gè)偶極接收器采集的偶極波列。

表1給出了GBG1井陣列聲波測量數(shù)據(jù)在一個(gè)深度點(diǎn)的存儲(chǔ)格式。由表1可以看到，單極子聲源T1發(fā)射采集的8道全波波列，存儲(chǔ)在曲線TFWV10中，總采樣點(diǎn)數(shù)為5376，每道波列點(diǎn)數(shù)則為672。同樣地，單極子聲源T2發(fā)射采集的4道全波波列，存儲(chǔ)在曲線TNWV10中，總采樣點(diǎn)數(shù)為832，每道波列點(diǎn)數(shù)則為208。而TXXWV10、TXYWV10、TYXWV10和TYYWV10則分別為采集的4個(gè)偶極分量，每個(gè)分量存儲(chǔ)8道全波波列，每道波列采樣點(diǎn)數(shù)為400。

表1 GBG1井陣列聲波測量數(shù)據(jù)存儲(chǔ)格式

圖1給出了GBG1井TXYWV10曲線在3759.94 m深度處，采集的8道偶極波列?？梢钥吹?，各道全波波列具有較高的相似性。

圖1 GBG1井3759.94 m處 TXYWV10波形

如圖2所示，陣列聲波信號的數(shù)據(jù)上傳是通過遙測傳輸短節(jié)實(shí)現(xiàn)的。儀器采集的各道聲波全波列(一維時(shí)間信號)在井下儀器總線控制器的作用下，通過井下儀器總線，傳入遙測短節(jié)內(nèi)進(jìn)行緩存，形成上行數(shù)據(jù)幀。對數(shù)據(jù)幀進(jìn)行PCM編碼后，轉(zhuǎn)換為模擬信號，經(jīng)低通濾波和差分驅(qū)動(dòng)，送入與傳輸電纜相連的模式變壓器，產(chǎn)生電纜中電壓的高低變化。地面系統(tǒng)在檢測到電纜電壓的變化后，經(jīng)過信號濾波處理和模數(shù)轉(zhuǎn)換，以及PCM解碼，可接收到上傳的陣列聲波采集數(shù)據(jù)。由于高頻信號經(jīng)測井電纜長距離傳輸后，幅度大幅衰減且信號嚴(yán)重失真，造成測井電纜傳輸帶寬有限。以陣列聲波數(shù)據(jù)傳輸為例，由表1可知波列每個(gè)采樣點(diǎn)的值為2 B長度的有符號SHORT型，計(jì)算一個(gè)深度點(diǎn)的陣列聲波采樣數(shù)據(jù)幀大小至少為37.125 KB，電纜上傳速率為41.66 kbit/s，則數(shù)據(jù)傳輸時(shí)間約為0.8911 s，加上聲源發(fā)射和陣列接收換能器采樣時(shí)間約為0.2 s，在不計(jì)入其他一維、二維采集數(shù)據(jù)的情況下，一次采集加數(shù)據(jù)上傳的時(shí)間至少為1.0911 s，而陣列聲波的深度采樣間隔為0.1524 m，可知儀器在單位深度地層的工作時(shí)長最小約為7.16 s。因此，在數(shù)據(jù)上傳前，對采集數(shù)據(jù)進(jìn)行有效壓縮編碼，可大幅減少數(shù)據(jù)上傳時(shí)間，變相增加單位時(shí)間內(nèi)傳輸?shù)臄?shù)據(jù)量和電纜上傳帶寬，進(jìn)而減小儀器在單位深度地層的工作時(shí)長。

圖2 陣列聲波測井儀數(shù)據(jù)傳輸原理框圖

2 陣列聲波采集數(shù)據(jù)無損壓縮方法

2.1 信號稀疏表示原理

由信號稀疏表示理論可知，如果原始信號是稀疏的，則可以將其用一稀疏正交基和非零個(gè)數(shù)遠(yuǎn)小于信號序列長度的稀疏系數(shù)序列線性組合表示，即稀疏變換[12]；而作為稀疏變換矩陣的稀疏正交基，稱為稀疏字典，可拓展為超完備的冗余字典。對原始信號f進(jìn)行稀疏變換：

α=Df

(1)

式中:f∈RN×1為N個(gè)采樣點(diǎn)的等間隔采樣原始信號；D∈RM×N為稀疏變換矩陣；α∈RM×1為f進(jìn)行稀疏變換的稀疏表示系數(shù)，如果M=N，則D為稀疏字典；如果M?N，則D為超完備冗余字典。反變換有f=DHα，DH為D的共軛轉(zhuǎn)置陣，滿足DHD=1。

(2)

2.2 基于稀疏表示的聲波全波列數(shù)據(jù)壓縮方法

基于稀疏表示的聲波全波列數(shù)據(jù)壓縮方法原理如圖3所示。

圖3 聲波全波列數(shù)據(jù)無損壓縮方法原理框圖

首先，提取一道的聲波全波列數(shù)據(jù)，去除基值后，根據(jù)采用K-SVD算法預(yù)先構(gòu)建的過完備冗余字典B和正交匹配追蹤算法，計(jì)算稀疏表示系數(shù)向量A，對向量A中的數(shù)值取整后，再根據(jù)稀疏反變換，重構(gòu)波列，得到原始波列與重構(gòu)波列的誤差向量err。將取整后的向量A和誤差向量err中的非零值、索引和基值Base整合，形成二進(jìn)制數(shù)據(jù)流，采用自適應(yīng)算術(shù)編碼后上傳。在地面系統(tǒng)，根據(jù)解碼后的稀疏表示系數(shù)和字典B進(jìn)行稀疏反變換，得到波列估計(jì)值，再加上解碼后的誤差和基值，就可以完全恢復(fù)出原始的聲波全波列。由于原始波列值為16位整型數(shù)，因此向量A和誤差向量err中的非零值和索引值同樣為16位整型。

算術(shù)編碼是一種信息保持型編碼，無需向LZW算法傳送編碼字典，因此壓縮率有了明顯提升。算術(shù)編碼符合概率匹配原則：出現(xiàn)概率較大的符號壓縮為較短的碼字，出現(xiàn)概率較小的符號取較長碼字。自適應(yīng)算術(shù)編碼根據(jù)編碼時(shí)符號出現(xiàn)的頻繁程度動(dòng)態(tài)修改信源符號的概率。當(dāng)出現(xiàn)某個(gè)符號在整個(gè)符號流中全局概率較小，但在符號流某一位置集中分布的情況時(shí)，宜采用自適應(yīng)編碼；由于符號出現(xiàn)的概率隨著符號輸入不斷變化，輸入相同符號數(shù)越多，該符號的概率被更新得越大，則該符號壓縮后的碼字要明顯小于需要遍歷的整個(gè)符號流，統(tǒng)計(jì)每個(gè)符號全局概率的固定模式算術(shù)編碼。因此，采用編碼過程中動(dòng)態(tài)估計(jì)字符概率的自適應(yīng)算術(shù)編碼不僅在壓縮算法時(shí)效上優(yōu)勢明顯，而且壓縮率也優(yōu)于固定模式編碼。

2.3 基于K-SVD的稀疏變換矩陣訓(xùn)練方法

信號的稀疏度與超完備冗余字典(稀疏變換矩陣)密切相關(guān)，選用符合聲波波列特征的稀疏變換矩陣，可使聲波波列在該稀疏變換矩陣下表示更稀疏。因此，采用K-SVD算法[15]，以歷史多道聲波全波列為訓(xùn)練數(shù)據(jù)，預(yù)構(gòu)建更符合聲波波列特征的過完備冗余字典，可在提高聲波波列稀疏度的基礎(chǔ)上，進(jìn)一步提升聲波波列的重構(gòu)精度，減少誤差向量非零值數(shù)目，從而使編碼數(shù)據(jù)量更小。

稀疏變換矩陣D∈RN×L(L>N)，矩陣中列向量為字典原子，訓(xùn)練數(shù)據(jù)集X∈RN×S，X中的列向量即為一道訓(xùn)練用的聲波全波列數(shù)據(jù)，其中L為字典原子個(gè)數(shù)，S為訓(xùn)練用聲波全波列道數(shù)。訓(xùn)練過程中稀疏表示系數(shù)向量集為?！蔙L×S。K-SVD算法的目標(biāo)函數(shù)為

(3)

式中:γi為向量集Γ中第i個(gè)列向量；K為稀疏表示系數(shù)向量中非零分量個(gè)數(shù)的最大值。K-SVD算法的訓(xùn)練為一個(gè)迭代過程，首先給定初始矩陣D，采用正交匹配追蹤算法求解訓(xùn)練數(shù)據(jù)集X的稀疏表示系數(shù)矩陣Γ；再根據(jù)Γ，更新矩陣D。重復(fù)進(jìn)行稀疏表示系數(shù)矩陣求解和稀疏變換矩陣更新兩個(gè)步驟，直至滿足式(3)的迭代條件。稀疏變換矩陣需要逐一更新矩陣D中的列。假設(shè)更新D中的第j個(gè)列向量dj，計(jì)算去掉該列代表的原子成分在所有S個(gè)訓(xùn)練樣本中造成的誤差

Ej=UΔVT

式中:U為N×N維矩陣;V為S×S維矩陣;U和V是酉矩陣，與本身轉(zhuǎn)置矩陣相乘等于單位陣;Δ為N×S維矩陣，其主對角線上的元素為奇異值，且除主對角線上的元素以外全為零值。更新dj為U的第1列，計(jì)算V的第1列和Δ(1,1)的乘積，更新矩陣Γ中與列向量dj對應(yīng)的第j行中非零元素值。

3 實(shí)驗(yàn)結(jié)果分析和現(xiàn)場測試

3.1 仿真實(shí)驗(yàn)結(jié)果和分析

在仿真實(shí)驗(yàn)中采用COMSOL軟件，模擬陣列聲波儀器分別在快速地層(花崗巖)、標(biāo)準(zhǔn)地層(砂巖)和慢速地層(泥巖)的陣列聲波單極子接收換能器響應(yīng)。其中，模擬儀器源距為3.0 m，接收換能器間距為0.15 m，聲源采用中心頻率為10 kHz的Ricker子波。圖4給出了快速地層中儀器8個(gè)接收換能器響應(yīng)波形。

圖4 快速地層中8個(gè)接收換能器響應(yīng)波形

如圖4所示，最上部的虛線波形為聲源發(fā)射波形，源距3～4 m處從上至下分別為8道接收換能器波形。實(shí)驗(yàn)中，共模擬3種地層接收波形176道，采用175道全波列波形構(gòu)建K-SVD算法的訓(xùn)練樣本集。初始過完備冗余字典采用DCT矩陣，原子數(shù)設(shè)為128，訓(xùn)練稀疏度上限為50，訓(xùn)練迭代次數(shù)為100。圖5和圖6分別給出了1道未參與訓(xùn)練的聲波波列采用正交匹配追蹤算法在訓(xùn)練后冗余字典和原始Hadamard字典上的重構(gòu)波形、稀疏表示系數(shù)和與原始波列的誤差。

圖5 仿真實(shí)驗(yàn)訓(xùn)練后冗余字典上的波列重構(gòu)結(jié)果

圖6 仿真實(shí)驗(yàn)原始Hadamard字典上的波列重構(gòu)結(jié)果

這里，換能器接收波列長度為501，原始Hadamard字典的原子數(shù)為512。由圖5和圖6可以看到，重構(gòu)波列和原始波列基本一致，但訓(xùn)練后冗余字典上的波列稀疏度為49，明顯優(yōu)于原始Hadamard字典上的波列稀疏度398。由于Hadamard字典的原子數(shù)大于K-SVD算法的過完備冗余字典原子數(shù)，因此原始Hadamard字典上的波列誤差向量為0，優(yōu)于訓(xùn)練后冗余字典上的誤差向量。模擬數(shù)據(jù)壓縮過程，一道原始波列為501×2 B，而經(jīng)稀疏變換后，稀疏表示系數(shù)的非零值與索引值占49×4 B，加上非零誤差值和索引的6×4 B，再加上波形基值2 B和標(biāo)準(zhǔn)化用到的波形最大和最小值4 B，共226 B，不考慮二進(jìn)制數(shù)據(jù)流最后的自適應(yīng)算術(shù)編碼，則一道原始波列的最小壓縮率可達(dá)到77.45%。仿真實(shí)驗(yàn)驗(yàn)證了采用K-SVD算法訓(xùn)練冗余字典的有效性。

3.2 實(shí)際測井?dāng)?shù)據(jù)實(shí)驗(yàn)結(jié)果和分析

本文采用GBG1井TFWV10曲線中4000～4038 m深度范圍內(nèi)的2000道聲波全波列數(shù)據(jù)構(gòu)建K-SVD算法的訓(xùn)練樣本集，初始過完備冗余字典采用Hadamard矩陣，原子數(shù)設(shè)為1024，訓(xùn)練稀疏度上限為50，訓(xùn)練迭代次數(shù)為100。圖7給出了該訓(xùn)練樣本集上的迭代誤差，這里對每個(gè)訓(xùn)練樣本進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化，范圍為[-1，1]。由圖7可以看到，在保證聲波全波列稀疏變換稀疏度的情況下，隨著字典更新次數(shù)的增加，算法的均方根誤差逐步減少，說明了K-SVD算法的收斂性。

圖7 K-SVD算法迭代誤差

隨機(jī)提取GBG1井TFWV10曲線中在訓(xùn)練用全波波形數(shù)據(jù)深度范圍外的一道聲波波列，進(jìn)行數(shù)據(jù)無損壓縮算法實(shí)驗(yàn)，驗(yàn)證K-SVD算法的有效性。圖8和圖9分別給出了該聲波波列采用正交匹配追蹤算法在訓(xùn)練后冗余字典和原始Hadamard字典上的重構(gòu)波形、稀疏表示系數(shù)和與原始波列的誤差。

圖8 訓(xùn)練后冗余字典上的波列重構(gòu)結(jié)果

圖9 原始Hadamard字典上的波列重構(gòu)結(jié)果

由圖8和圖9可以看到，在兩個(gè)冗余字典上的重構(gòu)波列與原始波列基本一致，原始波列和重構(gòu)波列的誤差向量(圖8(c)和圖9(c))也驗(yàn)證了這一點(diǎn)，統(tǒng)計(jì)平均誤差分別為0.98%和1.65%；雖然兩個(gè)誤差近似，但比較圖8(b)和圖9(b)的稀疏表示系數(shù)向量，發(fā)現(xiàn)在訓(xùn)練后冗余字典上的稀疏度明顯優(yōu)于原始Hadamard字典；在訓(xùn)練后冗余字典上的稀疏表示系數(shù)向量中非零值個(gè)數(shù)為94，而在原始Hadamard字典上的相應(yīng)值為444，驗(yàn)證了筆者提出的無損壓縮算法中對于冗余字典進(jìn)行預(yù)訓(xùn)練步驟的必要性和有效性。模擬數(shù)據(jù)壓縮過程，一道原始波列為660×2 B，而經(jīng)稀疏變換后，稀疏表示系數(shù)的非零值與索引值占94×4 B，加上非零誤差值和索引的37×4 B，再加上波形基值2 B和標(biāo)準(zhǔn)化用到的波形最大和最小值4 B，共530 B，不考慮二進(jìn)制數(shù)據(jù)流最后的自適應(yīng)算術(shù)編碼，則一道原始波列的最小壓縮率可達(dá)到59.85%。

針對HB油田FY1井等3口井的實(shí)測陣列聲波數(shù)據(jù)，比較本算法與LZW算法、算術(shù)編碼和文獻(xiàn)[11]算法的數(shù)據(jù)壓縮率，結(jié)果如表2所示。其中，算法分別針對TFWV10、TNWV10兩種單極子波列和偶極接收波列，訓(xùn)練3種過完備冗余字典。

由表2可以看到，筆者算法的壓縮率在不同數(shù)據(jù)類型上是最高的，較文獻(xiàn)[11]算法的壓縮率平均提升了約17.3%。其中，隨著數(shù)據(jù)文件變大，LZW算法的壓縮率有所提升，這是由于數(shù)據(jù)流中重復(fù)字符串的增多，使得字典中對應(yīng)字符串編碼縮短導(dǎo)致的。對于筆者算法，可以看到，在TNWV10數(shù)據(jù)上的壓縮率略優(yōu)于在TFWV10和偶極波列數(shù)據(jù)上的壓縮率，這是由于TNWV10數(shù)據(jù)的全波波列長度小于256，因此在壓縮時(shí)，誤差的索引值只占用1 B，導(dǎo)致碼流長度減少，壓縮率增加。同時(shí)，P273井4種方法的壓縮率較其他兩口井在整體上有明顯提升。經(jīng)該井實(shí)驗(yàn)數(shù)據(jù)研究發(fā)現(xiàn)，測量井段對應(yīng)地層中出現(xiàn)各向異性極弱的大段泥巖層，使得單極子接收波列和偶極子接收波列各道全波波列具有高度相似性，則LZW算法構(gòu)造的編碼字典大幅度減小，算術(shù)編碼由于字符出現(xiàn)概率增大，編碼碼字縮短，導(dǎo)致兩種算法的壓縮率明顯增加。由于二進(jìn)制碼流中，字符變化的規(guī)律性增強(qiáng)，也使得文獻(xiàn)[11]算法的壓縮率增加。對于筆者算法，波列的高度相似性導(dǎo)致訓(xùn)練的冗余字典更能表征波列特征，使得壓縮時(shí)稀疏變換的稀疏度和重構(gòu)精度增加，壓縮編碼長度進(jìn)一步降低。

表2 算法壓縮率實(shí)驗(yàn)結(jié)果

3.3 實(shí)際作業(yè)應(yīng)用測試

重新設(shè)計(jì)的SL6680型陣列聲波測井儀，應(yīng)用無損壓縮傳輸方法在HB油田GuGu8井等測井作業(yè)井眼環(huán)境較好的4口井進(jìn)行了陣列聲波實(shí)際作業(yè)應(yīng)用測試。在應(yīng)用測試過程中，對于同一口井分別采用原有的SL6680型儀器和新系統(tǒng)進(jìn)行測井時(shí)效對比測試，結(jié)果如表3所示。

表3 測井時(shí)效對比

由表3可以看出，原有儀器和新系統(tǒng)的測量數(shù)據(jù)平均誤差在測井重復(fù)測量標(biāo)準(zhǔn)規(guī)定的數(shù)據(jù)誤差范圍內(nèi)，說明兩者的測量數(shù)據(jù)幾乎相同。在這種情況下，新系統(tǒng)的最高測井速度整體高于原有儀器，測井作業(yè)耗時(shí)大幅度縮短。經(jīng)數(shù)據(jù)統(tǒng)計(jì)計(jì)算，新系統(tǒng)的測井效率較原有儀器平均提升了約20.2%。

4 結(jié)束語

針對陣列聲波測井儀器，對采集的多道單極子和偶極子全波波列，采用基于信號稀疏表示的數(shù)據(jù)無損壓縮傳輸算法，在數(shù)據(jù)傳輸帶寬受限的情況下，有效提升單位時(shí)間內(nèi)的上傳信息量，最終提高了測井時(shí)效，使陣列聲波測井儀器在保證采集數(shù)據(jù)不失真的條件下，滿足油田勘探提質(zhì)提效的形勢需求。