唐新閏, 劉彥彤, 張 巖, 趙玉瑩, 關(guān)正昊

(1.東北石油大學(xué) 計(jì)算機(jī)與信息技術(shù)學(xué)院, 黑龍江 大慶163318; 2.大慶油田勘探開(kāi)發(fā)研究院 信息研究室, 黑龍江 大慶163000)

0 引 言

巖心是油氣田勘探開(kāi)發(fā)中重要的基礎(chǔ)地質(zhì)資料, 對(duì)其觀察描述在確定巖性、 生儲(chǔ)蓋組合研究以及推斷沉積環(huán)境中具有重要作用。 隨著巖心圖像采集設(shè)備的廣泛使用, 巖心樣本被掃描成數(shù)字形式存儲(chǔ), 在提升油氣田勘探開(kāi)發(fā)水平方面起到重要作用, 是實(shí)現(xiàn)油田數(shù)字化建設(shè)必不可少的內(nèi)容。 由于多年的累積和不斷新取心, 導(dǎo)致巖心數(shù)據(jù)量極其巨大。 因此, 研究適合巖心圖像特點(diǎn)的壓縮算法是非常有意義的[1]。 通過(guò)對(duì)大量典型巖心圖像分析發(fā)現(xiàn), 特殊的地質(zhì)環(huán)境和復(fù)雜的地質(zhì)變遷導(dǎo)致巖心圖像普遍具有紋理信息豐富、 對(duì)比度較弱的特點(diǎn), 傳統(tǒng)壓縮方法受Nyquist 采樣定理的限制, 在低碼率下影響邊緣、 輪廓等高維奇異區(qū)域的保持效果[2-3]。 Donoho[4]和Candès[5]提出了壓縮感知( CS: Compressed Sensing)理論與稀疏表示方法結(jié)合, 同時(shí)進(jìn)行壓縮與采樣操作, 在信號(hào)具有可壓縮性或稀疏性先驗(yàn)的條件下, 通過(guò)采集少量信號(hào)的投影值, 在接收端通過(guò)稀疏性約束方法就可以獲得信號(hào)準(zhǔn)確的重構(gòu), 在某種意義上突破了Nyquist 采樣定理的瓶頸, 該方法對(duì)研究巖心等紋理信息豐富的圖像壓縮與處理具有重要意義[6-7]。

在實(shí)際應(yīng)用領(lǐng)域直接采用CS 方法重構(gòu)整幅圖像將導(dǎo)致巨大的計(jì)算量。 分塊壓縮感知( BCS:Block Compressed Sensing) 方法在2007 年由Gan[8]提出, 該方法降低了存儲(chǔ)和計(jì)算成本。 基于此,Mun 等[9]提出了BCS-SPL(Smooth Projected Landweber)算法, 將圖像分塊進(jìn)行觀測(cè)采樣, 通過(guò)維納濾波結(jié)合Landweber 迭代實(shí)現(xiàn)重構(gòu), 使用離散余弦變換、 離散小波變換等正交基作為稀疏基, 改善了BCS的重構(gòu)效果, 但由于對(duì)各個(gè)圖像塊采用相同的采樣率, 造成資源分配不合理, 且重構(gòu)過(guò)程中存在一定塊效應(yīng)。 Fowler 等[10]提出多尺度分塊變采樣率壓縮感知算法(MS-BCS-SPL), 對(duì)每一層小波分解設(shè)定不同的采樣率, 改善了較低碼率情況下的塊效應(yīng)現(xiàn)象。 但由于不同子塊中的邊緣、 紋理信息可能不同, 對(duì)每一層子塊采用相同的采樣率并不合理。 因此國(guó)內(nèi)外專家與學(xué)者針對(duì)該問(wèn)題進(jìn)行了研究。李玉等[11]和 高東紅[12]分別提出了根據(jù)圖像邊緣結(jié)構(gòu)和紋理信息自適應(yīng)分配采樣率的方法, 在一定程度上減少了采樣數(shù)量, 提高了重建性能。 但其只對(duì)小波分解后的低頻系數(shù)進(jìn)行重建, 而對(duì)于細(xì)節(jié)較復(fù)雜的圖像重建效果不理想。 張學(xué)全等[13]根據(jù)圖像塊的內(nèi)部活動(dòng)性將圖像塊分為紋理塊和平坦塊, 通過(guò)鄰塊邊緣自適應(yīng)加權(quán)濾波方法消除了圖像重構(gòu)的塊效應(yīng), 提高了圖像邊緣的保持, 但是對(duì)平坦塊的檢測(cè)比較耗時(shí)導(dǎo)致算法處理速度較慢。 此外當(dāng)前圖像壓縮感知重構(gòu)采用的稀疏表示基函數(shù)通常有Fourier[14]、 Wavelet[15]等正交變換。 Fourier 是一種全局變換, 不具有局部化分析的能力, 容易產(chǎn)生吉布斯現(xiàn)象; Wavelet 是一種時(shí)頻局域性變換, 但在表示二維奇異線條情況下, 存在不光滑和模糊現(xiàn)象。基于K-SVD(K-Singular Value Decomposition)超完備字典學(xué)習(xí)的圖像稀疏表示方法由Bryt 等[16]提出,學(xué)習(xí)型超完備字典能根據(jù)圖像本身的特點(diǎn), 自適應(yīng)的調(diào)整變換基函數(shù), 該方法在圖像處理領(lǐng)域得到了廣泛應(yīng)用。

受上述研究的啟發(fā), 針對(duì)傳統(tǒng)方法進(jìn)行巖心圖像壓縮感知重構(gòu)時(shí), 在低碼率下容易產(chǎn)生細(xì)節(jié)丟失的問(wèn)題, 筆者提出了基于K-SVD 的重構(gòu)算法, 并在BCS 算法的基礎(chǔ)上, 采用高斯隨機(jī)矩陣觀測(cè)相應(yīng)層級(jí)的圖像塊, 用MMSE 估計(jì)近似解獲得初始解, 采用K-SVD 超完備學(xué)習(xí)字典稀疏表示巖心圖像, 結(jié)合Landweber 迭代與自適應(yīng)閾值方法實(shí)現(xiàn)巖心圖像的壓縮和重構(gòu)。 通過(guò)仿真實(shí)驗(yàn)結(jié)果證明, 筆者方法重構(gòu)的巖心圖像能更好的保留巖心圖像的紋理細(xì)節(jié), 提高了圖像的整體質(zhì)量。

1 分塊壓縮感知重建算法

根據(jù)壓縮感知理論, 假設(shè)x∈RN是長(zhǎng)度為N 的原始信號(hào), y∈RM是長(zhǎng)度為M 的采樣信號(hào), 且M?N,則觀測(cè)值y 可用

表示。 其中Φ 是一個(gè)具有采樣率為r=M/ N 的M×N 維的觀測(cè)矩陣。

由分塊壓縮感知方法, 假設(shè)一幅大小為N×N 的圖像x 被分成大小為B×B 的圖像塊, 第i 個(gè)巖心圖像塊向量化后表示記為xi, 在觀測(cè)矩陣ΦB下得到的觀測(cè)值

其中i=1,…,n, n =N/ B, B 的尺寸根據(jù)巖心圖像重構(gòu)的效率決定: 當(dāng)B 較小時(shí), 內(nèi)存占用相對(duì)少且計(jì)算速度快; 當(dāng)B 較大時(shí), 圖像重構(gòu)的質(zhì)量相對(duì)高。 ΦB是大小為MB×B 的正交觀測(cè)矩陣, MB=(M×B) / N, M為整幅圖像的觀測(cè)采樣值數(shù)目。 整個(gè)圖像的觀測(cè)矩陣Φ 為對(duì)角矩陣, 如下所示。

分塊后觀測(cè)矩陣ΦB的維數(shù)會(huì)相對(duì)降低, 顯然式(1)是欠定的方程組, 無(wú)法根據(jù)y 求解x, 當(dāng)具備圖像數(shù)據(jù)x 在某個(gè)變換域C 下系數(shù)α 是稀疏的先驗(yàn), 可以通過(guò)求解L1范數(shù)優(yōu)化問(wèn)題

2018年9月17日，美國(guó)森圖斯能源公司(Centrus Energy)和韓國(guó)斗山重工建設(shè)有限公司(Doosan Heavy Industries and Construction)宣布雙方已簽署一份諒解備忘錄，未來(lái)將在先進(jìn)反應(yīng)堆領(lǐng)域開(kāi)展合作。

求解式(1)。 轉(zhuǎn)化為求解凸優(yōu)化的問(wèn)題, 由于α 是具有稀疏性的先驗(yàn), 可以利用不斷促進(jìn)α 稀疏性的方法求解, 即有

其中λ 是平衡因子, α 是α 的估計(jì)。 在促進(jìn)α 稀疏性的過(guò)程可采用Landweber 冷卻閾值迭代法, 其過(guò)程如下

標(biāo)準(zhǔn)Landweber 算法迭代重構(gòu)圖像亮度值可能是負(fù)數(shù), 在選擇松馳參數(shù)及控制收斂速度方面難度較大, Piana 等[17]為此提出改進(jìn)的Landweber 算法, 在重構(gòu)過(guò)程中利用先驗(yàn), 得到投影Landweber 的方法,然后在傳統(tǒng)Landweber 迭代過(guò)程中利用維納濾波較弱塊效應(yīng)現(xiàn)象, 該方法不僅考慮信號(hào)的稀疏性還具有平滑的作用, 即在空間域內(nèi)進(jìn)行維納濾波操作, 在變換域內(nèi)實(shí)現(xiàn)平滑和閾值控制。

2 基于K-SVD 的壓縮感知重構(gòu)算法

2.1 自適應(yīng)閾值算法

在分塊壓縮感知重構(gòu)的基礎(chǔ)上, 對(duì)巖心圖像壓縮感知重構(gòu)操作中的實(shí)驗(yàn)參數(shù)進(jìn)行研究, 發(fā)現(xiàn)Canny算法[18]等方法的閾值參數(shù)通常依據(jù)數(shù)據(jù)的先驗(yàn)知識(shí)或經(jīng)驗(yàn)設(shè)定, 不能得到自適應(yīng)圖像特點(diǎn)的最佳閾值。全局和局部適應(yīng)閾值法是目前廣泛使用的計(jì)算閾值方法, 其中全局閾值時(shí)間復(fù)雜度較低, 與信號(hào)尺寸對(duì)數(shù)的平方根成正比。 巖心圖像分塊壓縮感知初始解圖像塊的尺寸相對(duì)較小, 克服了因尺寸較大出現(xiàn)的“過(guò)扼殺” 系數(shù)后果, 所以選擇計(jì)算巖心圖像BCS 初始解圖像塊全局閾值, 實(shí)現(xiàn)根據(jù)巖心圖像的自身特征計(jì)算閾值的重構(gòu)操作[19], 計(jì)算公式如下

其中ˉσ 是初始解中不同巖心圖像塊小波系數(shù)標(biāo)準(zhǔn)差的數(shù)學(xué)期望, M,N 為初始解的分解尺度。

依據(jù)式(7)得到閾值T, 另外需要具有篩選作用的閾值控制函數(shù), 閾值控制函數(shù)中硬閾值函數(shù)與軟閾值函數(shù)較為常見(jiàn)。 設(shè)信號(hào)數(shù)據(jù)可以由ω 表示, 閾值用T 表示, 符號(hào)函數(shù)表示為sgn(*), 軟閾值函數(shù)為

由式(8)可見(jiàn), 軟閾值函數(shù)首先比較ω 和T, 然后根據(jù)比較的結(jié)果再向0 進(jìn)行收縮, 相比之下更符合極小與極大準(zhǔn)則的理想值, 因此可以結(jié)合T 與軟閾值控制函數(shù)實(shí)現(xiàn)篩選功能。

2.2 K-SVD 字典學(xué)習(xí)算法

目前在數(shù)字圖像超完備字典學(xué)習(xí)方法中, K-SVD 表現(xiàn)出較高的稀疏表示效果[20], 在實(shí)際應(yīng)用中具有高效、 靈活的特點(diǎn), 因此筆者采用K-SVD 自適應(yīng)學(xué)習(xí)獲得巖心圖像的稀疏表示基。

K-SVD 主要分兩步操作: 第1 步為稀疏編碼; 第2 步為字典更新。 令分別代表信號(hào)、 字典和稀疏表示向量,是N個(gè)信號(hào)集合,是S解向量集合, K-SVD 算法的目標(biāo)方程如下

通過(guò)K-SVD 訓(xùn)練超完備字典的步驟如下。

1) 設(shè)定D0為超完備字典初始值: 超完備字典D0∈RN×M中的原子可以通過(guò)從信號(hào)樣本x中任選的k個(gè), 令k=0, 即D0為0 矩陣。

2) 稀疏編碼階段: 利用OMP(Orthogonal Matching Pursuit)算法[21]通過(guò)解t.‖y‖0≤k0, 以獲得稀疏表達(dá)^yi(1≤i≤M), 并構(gòu)成稀疏表達(dá)集y={y1,y2,…,yM}。

3) K-SVD 字典學(xué)習(xí)步驟: 令j0=1,2,…,m表示將要更新的原子dj0, 開(kāi)始j0= 1, 從y中找出所使用原子dj0的列, 即在j0行中非零稀疏系數(shù)所在的位置

其中為y(k)矩陣的第j行, 依據(jù)Ωj0從Ej0中選擇對(duì)應(yīng)的列得, 通過(guò)SVD 算法得=UΔVT。 更新字典中的原子dj0=ui與對(duì)應(yīng)的稀疏表示yR j0=Δ[1,1]·vi。

4) 當(dāng)‖x-D(k)y2(k)F‖22≤ε時(shí), 停止更新。

對(duì)巖心圖像, 利用標(biāo)準(zhǔn)差σˉ 和初始解尺度獲得自適應(yīng)閾值, 通過(guò)K-SVD 算法學(xué)習(xí)訓(xùn)練稀疏表示巖心圖像的超完備字典, 得到巖心圖像小波系數(shù)的標(biāo)準(zhǔn)差, 圖1 給出了兩幅巖心圖像與稀疏字典圖示。

圖1 訓(xùn)練巖心圖像與字典原子Fig.1 Training images and dictionaries

3 算法描述

首先基于分塊壓縮感知, 將巖心圖像劃分為尺寸相同的圖像塊, 通過(guò)高斯隨機(jī)矩陣與每個(gè)圖像塊相作用, 獲得對(duì)應(yīng)的觀測(cè)值。 接著利用MMSE 方法得到初始解的估計(jì)值, 然后采用提升小波基對(duì)初始解進(jìn)行分解, 依據(jù)式(7)得到自適應(yīng)閾值。 最后利用K-SVD 字典學(xué)習(xí)方法自適應(yīng)閾值以及維納濾波, 迭代求解, 實(shí)現(xiàn)巖心圖像壓縮編碼和重構(gòu)。 整體算法的具體步驟如下:

1) 讀取原始巖心圖像, 設(shè)定圖像分塊大小和采樣率r等參數(shù), 生成巖心圖像子塊與高斯隨機(jī)矩陣ΦB;

2) 用ΦB作用于每個(gè)巖心圖像塊xi, 進(jìn)行觀測(cè), 得到相應(yīng)的觀測(cè)值yi;

3) 用MMSE 方法獲得巖心圖像初始解的估計(jì), 此時(shí)k=0;

4) 在初始解的基礎(chǔ)上使用3×3Winner 的濾波器進(jìn)行k次平滑操作;

5) 對(duì)初始解使用提升小波基進(jìn)行分解, 獲得小波系數(shù)標(biāo)準(zhǔn)差σi, 依據(jù)式(7)得到自適應(yīng)的閾值T;

6) 設(shè)定迭代次數(shù)k和超完備字典D原子的個(gè)數(shù)Dn;

7) 對(duì)每一巖心圖像塊xi, 使用迭代;

8) 采用K-SVD 算法訓(xùn)練yi, 獲得大小B2×Dn的超完備字典D;

9) 依據(jù)軟閾值函數(shù)與T、 字典D重構(gòu)巖心圖像子塊~xi;

11) 通過(guò)Oi+1=‖xi+1-x^i‖2計(jì)算殘差, 其中i=i+1;

12) 當(dāng)Oi-Oi-1＜10-2時(shí), 得到重構(gòu)的最優(yōu)解, 迭代終止, 否則重新返回步驟7)。

4 實(shí)驗(yàn)結(jié)果及分析

實(shí)驗(yàn)環(huán)境采用Windows7 操作系統(tǒng), 軟件平臺(tái)為Matlab R2013b 環(huán)境, 實(shí)現(xiàn)巖心圖像的壓縮與重構(gòu),并與BCS-SPL-DWT 和BCS-SPL-DDWT 算法的實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。 測(cè)試巖心采用尺寸為512×512 像素的圖像, 通過(guò)高斯隨機(jī)矩陣構(gòu)造采樣觀測(cè)矩陣, 整體巖心圖像劃分成尺寸為B×B的子塊, 超完備余字典D的原子個(gè)數(shù)為Dn(實(shí)驗(yàn)中取B=8,Dn=4), 則得到大小為B2×Dn的超完備字典D, 分別用離散小波變換(DWT: Discrete Wavelet Transformation)、 雙樹(shù)復(fù)小波變換(DDWT: Dual-plex Wavelet Transform)和K-SVD字典作為稀疏表示基, DWT 和DDWT 分別采用3 級(jí)分解。 表1 給出2 幅巖心圖像在3 種不同方案下壓縮與重構(gòu)PSNR 值。 可見(jiàn)與其他兩種方案實(shí)現(xiàn)巖心圖像的壓縮與重構(gòu)的PSNR 值相比, 筆者算法的PSNR提高0.1 ～0.8 dB。

表1 3 種不同方案對(duì)巖心圖像壓縮與重構(gòu)的峰值信噪比PSNR Tab.1 PSNR of the core image of compression and reconstruction with three schemes (dB)

圖2 給出了采樣率r=0.25 時(shí)3 種方法的結(jié)果。 由圖2 可見(jiàn), BCS-SPL-KSVD 重構(gòu)的巖心圖像質(zhì)量相較BCS-SPL-DWT 約提高0.3 ～0.6 dB, 與BCS-SPL-DDWT 相比約提高0.1 ～0.4 dB。

圖2 巖心1 的3 種方法對(duì)比(r=0.25)Fig.2 Comparison of three reconstruction methods for core image one

圖3 巖心2 的3 種方法對(duì)比(r=0.25)Fig.3 Comparison of three methods for reconstructing core image two

5 結(jié) 語(yǔ)

筆者在分塊壓縮感知重構(gòu)的框架下, 根據(jù)巖心圖像的特點(diǎn), 以全局閾值作為自適應(yīng)閾值的依據(jù),結(jié)合K-SVD 字典學(xué)習(xí)方法, 實(shí)現(xiàn)巖心圖像的隨機(jī)壓縮采樣觀測(cè)與重構(gòu)。 仿真實(shí)驗(yàn)證明, 筆者算法得到的巖心圖像視覺(jué)效果和峰值信噪比均有所提升, 與BCS-SPL-DWT 和BCS-SPL-DDWT 算法相比, PSNR 值分別平均提高0.6 dB 和0.4 dB, 并且有效保留了巖心圖像的紋理特征, 實(shí)現(xiàn)了對(duì)巖心圖像較高質(zhì)量的重構(gòu)。