雙層非負矩陣分解的分形圖像壓縮算法

方美東，王 輝，張愛華

南京郵電大學 理學院，南京 210023

如果一個對象的部分與整體以某種規則聯系起來，通過某種變換使之對應，那么可將它看成分形。自1988年Barnsley首次提出分形圖像壓縮編碼的概念以來，分形圖像壓縮編碼已獲得廣泛認可[1-2]，并成為近年來最受歡迎的現代圖像壓縮方法之一[3]。之后，基于塊的分形圖像編碼方法被Jacquin[4]和Barnsley[5]提出。分形塊編碼方法[6]是在自然圖像具有豐富仿射冗余的情況下，發現一個構造規則來構造近似圖像。

分形圖像壓縮將圖像分成Range塊和Domain塊，然后在Domain塊池中搜索與之最佳匹配的Range塊。該方法在編碼過程中只記錄匹配過程中必要的變換參數[7-9]。在分形圖像壓縮過程中，關鍵性問題是如何優化匹配方案，減少壓縮時間，提高重建圖像的質量。Huang等[10]提出了一種基于分形特征的圖像檢索方法。

NMF本質上是一種無監督的特征提取算法，其思想最早可追溯到Paatero等[11]在1994年提出的正矩陣分解的概念。Lee和Seung于1999年在Nature上正式提出了NMF方法和概念[12]。近年來，NMF在人臉識別[13-15]及特征提取[16]等領域得到廣泛應用。

NMF收斂效果不夠明顯，Hoyer首次提出將稀疏約束的思想加入NMF[17]，接著Hoyer提出了一種介于L1范數和L2范數之間的稀疏度量[18]；Du等提出了L3/2范數的稀疏約束NMF[19]和基于L3/2范數的圖正則NMF[20]；Yang等在NMF基礎上增加了投影約束限制[21]。本文在投影非負矩陣分解的基礎上加入L3/2范數的稀疏約束，構建了雙層非負矩陣分解，提取圖像的稀疏特征灰度值，然后對特征進行K-means聚類，接著進行基于分類下的稀疏分解，相較于其他算法，本文改進算法在分形圖像壓縮的過程中提高了重建圖像的效率。

1 理論基礎

1.1 非負矩陣分解

非負矩陣分解是一種非常有效的低秩逼近方法，它通過對非負矩陣進行非負因子分解得到數據的潛在特征。對于給定數據為X∈R M×N，非負矩陣分解的目的是把它分解為基因子U∈R M×K與低維表示因子V∈R K×N乘積的形式：

應用兩種測度方式：Frobenius范數（簡稱F范數）和Kullback-Liebler（KL）散度來度量式（1）逼近程度，基于歐式距離的目標函數為：

本文采用F范數求得分解后的非負基矩陣，針對目標函數（2），Lee和Seung給出了乘性迭代規則：

上述迭代算法（4）可以看成步長自學習的梯度下降算法。

1.2 基本分形編碼

在基本分形編碼中，首先將N×N原始圖像分割成不重疊的Range塊R i(i=1,2,…,Nr)和重疊的Domain塊D j(j=1,2,…,Nd)，其中通常取X=2B，δ是滑動窗口的長度。

對于每個Domain塊D j(j=1,2,…,N)，采用四領域像素平均，得到B×B像素塊，用S（空域收縮算子）表示這種運算，即：

圖1 范圍塊和域塊之間的縮進仿射變換的過程Fig.1 Process of contractive affine transformation between range block and domain block

然后，在灰度匹配過程中對Domain塊進行灰度變換?；叶茸儞Q是由亮度調整s和亮度偏移o兩部分組成的線性變換。因此，可以得到Range塊R與Domian塊D之間的灰度變換：

這里，R、D分別是Range塊和Domain塊的灰度值，|s|＜1，s與o都為實數，I是灰度均為1的大小為B×B的常值塊。根據最小誤差函數E(R,D)計算Range塊與Domain塊之間的匹配誤差：

其中‖·‖是L2范數。通過最小二乘法可得參數s和o如下：

其中·,·為內積，和表示Range塊和Domain塊的

解碼是一個相對簡單的迭代過程，基于壓縮變換的不動點定理和分形碼進行圖像重構。

1.3 稀疏分解

稀疏分解的目的是在海量的數據中選取一小部分重構新的數據[23]。在數學上，如果列向量B=[b1,b2,…,b K](b i∈R d)是有限維空間R中的一組基，則定義任意列向量X∈R d有：

其中S∈R K是向量X的系數和每個基向量上的投影系數。一般取B中的元素作為一組標準正交基，即任意兩個基向量b i、b j滿足：

如果基向量的個數等于空間維數，即K=d，那么這組基是完備的。稀疏分解的思想是利用少量的基向量來線性近似一個目標向量。因此，實現稀疏分解，需要基底完整，并且具有一定的冗余性?？梢酝ㄟ^基底數量遠遠大于基向量的維數來實現，即K?d。在稀疏分解中存在重構誤差，一般用稀疏近似，其形式為：

在保證重構精度的前提下，使系數向量達到稀疏最大化，有兩種解決方案，分別是匹配追蹤和正交匹配追蹤。OMP[24]（orthogonal matching pursuit）是一種改進的MP[25]（matching pursuit），在OMP中殘差與之前的線性部分保持完全正交，因此OMP具有更好的收斂性。

2 雙層非負矩陣分解的分形圖像壓縮算法

2.1 投影非負矩陣分解

PNMF在原矩陣空間和特征空間上建立線性映射，不僅可以降低求解的復雜度，還可以得到更稀疏的解。對于一個秩為K的對稱半正定矩陣P∈R M×M，存在非負矩陣U∈R M×K，使得P=UUT，將原始的非負矩陣X∈R M×N進行非負線性映射得到：

因此矩陣X可以表示為：

PNMF有兩種測度方式，分別是Frobenius范數和Kullback-Liebler（KL）散度。下面僅介紹基于歐式距離的目標函數：

根據目標函數式（14），更新規則如下：

相比NMF方法，PNMF方法將系數矩陣V∈R K×N替換為UTX，引入投影限制，并且在非負分解的過程中要求投影矩陣P接近單位矩陣。說明PNMF方法隱含了基矩陣的正交性，這樣可以獲得更稀疏的特征矩陣。

2.2 提出的非負矩陣分解

NMF將原始數據矩陣分解為兩個非負矩陣U和V的乘積，分解后的矩陣僅包含非負因子，并且基向量U可以很好地保持數據的局部結構。但是NMF得到的解不夠稀疏，收斂速度慢。PNMF可以將原始矩陣投影到特征空間，能很好地表示數據隱藏的屬性，對于正交限制的隱藏可以獲得更稀疏的矩陣。

對于范數稀疏約束方法而言，主要有：L0范數、L1/2范數、L1范數、L3/2范數、L2范數。

Lq是一組范數，根據q的變化，范數有著不同的性質。L0范數是NP難問題，而L1是L0的最優凸近似，相對L0范數要容易優化求解。L2范數可以實現正則化約束中的每一個元素都接近于零，而L1范數趨于產生稀疏的特征，其他元素都是零。通過L2范數實現正則化的同時在一定程度上避免了過擬合，從而提升泛化能力，還可以使得優化求解變得穩定和快速，一般用于目標函數求損失函數。

近年來，分數范數Lq(0＜q＜1)被很多學者提出，當分數q取不同值時，得到比L1范數更稀疏的解。在文獻[26]中，Krishnan等證明L1/2范數和L2/3范數正則化約束是非常高效的約束。在文獻[19]中，Du等表明Lq和L1/q具有相同屬性，當q＞1時，L3/2范數的正則化項是一個凸的光滑范數，驗證了L3/2范數約束在稀疏表示上具有很好的約束效果。

因此本文將PNMF與L3/2范數的稀疏約束相結合，提高矩陣分解的魯棒性和稀疏性，目標函數如下：

其中?、β是正則化參數（在仿真實驗中α、β均取10-3），目標函數可寫成如下形式：

下面對拉格朗日函數關于U、V求偏導：

根據KKT條件φik U ik=0，φkj V kj=0，使得：

進一步化簡，可以得到關于U、V的更新迭代公式：

使用梯度下降的方法對U、V進行如下更新：

其中ηik、?kj為步長參數，計算可得：

為了證明目標函數O2的收斂性，需要證明式（23）和（24）在更新過程中是不增的。由于U的證明過程和V的證明過程類似，這里只對V進行證明。將遵循類似于文獻[27]的描述過程。證明過程需要用到輔助函數，類似于期望最大化算法[28]中使用的函數。

定義1G(v,v′)是F(v)的輔助函數，如果G(v,v′)≥F(v)，則有G(v,v)=F(v)。

根據如下引理，輔助函數將非常有用。

引理1如果G是F的輔助函數，F在如下迭代中是不增的：

證明

現在，將證明式（24）中v的更新滿足式（29）中帶有輔助函數的更新。目標函數O2可表示為：

這里用Fab表示O2中只與v ab相關的部分，對于V中任意元素v ab，有：

由于V的更新本質上是元素的更新，這足以表明每一個F ab在更新步驟（24）下是不增的。

引理2函數：

是O2的部分輔助函數，只與V ab相關。

證明由于G(v,v)=Fab(v)是顯而易見，這里僅需要驗證為此比較Fab(v)的泰勒級數展開：

由于

并且

2.3 正交稀疏分形編碼過程

稀疏分解的編碼過程可以看成線性分解求系數矩陣的過程，因此，需要找到每個Range塊在虛擬碼本中的若干個基向量的最優線性組合。根據稀疏分解，一個新的灰度變換定義[29]如下：

通過進一步研究E(R,D)，可以發現：

根據以上推算過程可以得到R-Rˉ·I和D^能夠用來計算匹配誤差E(R,D)。通過直接計算R-Rˉ·I,D^2的最大值作為Domain塊的最優匹配塊?；诖朔椒?，基本分形編碼中的灰度變換公式（6）可以定義為如下形式：

其中sj(j=1,2,…,k(k≤K))是尺度參數。正交稀疏分形編碼方法直接計算尺度參數s，減少了計算的復雜度和迭代次數。

2.4 算法分析

根據上述計算，可以得到基于PNMF方法與L3/2稀疏約束結合的更新迭代過程。在進行分形圖像壓縮時，為了盡可能得到更稀疏的圖像特征，本文對圖像的灰度值矩陣進行了雙層非負矩陣分解，使分解結果盡可能的稀疏，獲得原始圖像的特征向量。結合DLNMF、Kmeans和OMP稀疏分解的思想，提出了一種雙層非負矩陣分解的分形圖像壓縮編碼算法。

本文采用二分K-means聚類[30]，可以有效地避免K-means聚類中在第一步中隨機選取初值中心敏感帶來的影響，保證了每一步誤差最小并且能夠克服Kmeans收斂于局部最優。利用二分K-means聚類將虛擬碼本Ω進行分類，根據稀疏分解的思想在相應的類別塊里進行正交稀疏分形編碼，不僅可以提高Range塊與Domain塊的匹配搜索效率，而且可以提高匹配重建圖像的質量。

2.5 算法實現

基于以上算法分析，下面給出所提出算法的具體過程。

算法1雙層非負矩陣分解的分形圖像壓縮算法

輸入：原始圖像I。

（1）圖像分割。把原圖分割成互不重疊的大小為n×n的固定Range塊，記為R塊。

（2）碼本構成。對同一幅圖像，在縱橫方向上按滑窗步長均為δ（一般的，取δ=2n）個像素來生成尺寸為2n×2n的Domain塊，記為D塊。對于每個D塊，采用4-鄰域像素值平均得到n×n的圖像塊，并考慮8種等距變換，這樣的子塊集合構成碼本Ω。

（3）計算圖像特征向量。對原始圖像的灰度值矩陣進行DLNMF。

①根據非負雙重奇異值分解初始化[31]得到初始矩陣U0和V0。

②首先，進行第一層的非負矩陣分解，根據目標函數：

得到基矩陣U1。然后進行第二層的非負矩陣分解，由目標函數：

得到基矩陣U2。

③對雙層非負矩陣分解的基矩陣U2進行歸一化，得到圖像的特征向量U3。

（4）根據原始圖像的特征向量U3進行K-means聚類。利用特征聚類索引將圖像的虛擬碼本Ω劃分成C類。

（5）用OMP算法進行稀疏灰度匹配處理。

對于每一個R塊R i，設置初始殘差e0=R i-Rˉi·I，索引集Λ0=?，初始化D塊B0=?，迭代次數t從1到K。

①找出殘差e t-1與相應類別字典D中某一列D^j的內積最大值對應的下標λ，即：

②更新索引集Λt=Λt-1?{}λt，更新所選字典子列構成的集合B t=[B t-1,D^λt]。

③利用最小二乘法得到K階逼近：

④更新殘差e t=e0-B t s^t。如果‖ ‖e t2＜emin，則迭代停止；否則，繼續步驟①。

（6）解碼。根據分形碼重建圖像。

3 實驗結果和分析

為了測試文中提出方法的效果，選擇標準測試圖像Woman2、Aerial、Barbara、Boat、Goldhill、Zelda、Milkdrop、Plane和Lake（512×512，8 bit量化）為實驗對象，在仿真實驗中，實驗平臺為運行Windows10家庭版的Intel?CoreTM（2.60 GHz CPU/16.0 GB內存），程序用MATLAB編寫。測試性能的參數為編碼時間T（單位：s）和峰值信噪比PSNR（單位：dB）。實驗中采用一致方塊分割，并選取R塊大小為8×8，生成D塊池滑窗步長是8。

3.1 評價標準

一般采用峰值信噪比PSNR來評價原始圖像和重建圖像間的相似度：

3.2 參數的選擇

在改進的方法中，OMP算法下的稀疏度K、OMP算法下的重構誤差閾值emin和K-means聚類過程中的聚類數C，這三個參數在編碼的過程中影響重建圖像的質量（PSNR）和編碼時間（Time）。因此，需要進一步確定算法中的三個參數得到最佳的實驗結果。

參數K是正交稀疏匹配過程中灰度值的稀疏程度，K值越大表示匹配獲得的灰度值越多，重建圖像的質量越好，但是由于數據變多，編碼過程中耗費的時長也會增加。根據圖2綜合峰值信噪比和編碼時間可以看出，當K=10時，獲得的編碼結果相對更好。

圖2 參數K對分形圖像壓縮編碼的影響Fig.2 Influence of parameter K on fractal image compression coding

參數emin是正交匹配搜索的終止條件，其目的是進一步區分出Range塊的紋理特征。對于紋理復雜的Range塊需要通過增加虛擬碼本的數量來提高重構質量。emin的值越大表示正交匹配的誤差閾限要求越低，導致重建的圖像質量降低，編碼時間減少。觀察圖3可以發現，當emin∈( )20,30時，PSNR和Time的變化比較快，本文中emin取25。

圖3 參數e min對分形圖像壓縮編碼的影響Fig.3 Influence of parameter e min on fractal image compression coding

參數C是K-means算法的聚類數，編碼時相應的圖像塊劃分的類別數越多，編碼的復雜度會增加，從而影響PSNR和Time。通過圖4觀察參數C取不同值時，圖像的重建質量有所下降，編碼時間也隨之增加，綜合PSNR和Time兩個因素，C=10時，重建圖像效果較好。

圖4 參數C對分形圖像壓縮編碼的影響Fig.4 Influence of parameter C on fractal image compression coding

3.3 結果和分析

使用Woman2、Aerial、Barbara、Boat、Goldhill、Zelda、Milkdrop、Plane和Lake圖像將本文算法分別與基本分形算法[32]、快速稀疏分形算法[33]和小波分形算法[8]進行比較，稀疏度K=10，重建誤差閾限emin=25，類別數C=10，實驗結果見表1～表3，表中用加粗字體標出了本文算法在重建圖像時峰值信噪比提高的值和編碼時間減少的值。

表1 基本分形算法和本文算法算法實驗對比Table 1 Experimental comparison between basic fractal algorithm and algorithm of this paper

表3 小波分形算法和本文算法實驗對比Table 3 Experimental comparison between wavelet fractal algorithm and algorithm of this paper

基本分形算法采用全搜索方法求解，雖然可以得到質量較好的重建圖像，但是耗費大量的時間。提出的方法在快速稀疏分形編碼的基礎上，利用NMF高效的無監督提取灰度值特征的性質，提出了DLNMF有效縮減碼本信息，并且根據提取的灰度特征將圖像塊進行了分類處理，提高編碼效率。從表1中可以看出改進方法在重建圖像質量上和編碼時間上都有很大的提高，其中最為明顯的是在基本分形算法中，Plan和Lake的重建圖像的峰值信噪比非常低，改進方法提高了30.61 dB和24.65 dB。

快速稀疏分形算法采用了正交匹配追蹤的思想，提高了Range塊與Domain塊的匹配效率，但是實驗結果表明該方法有待進一步完善。本文算法在快速稀疏算法的基礎上加入特征提取及分類匹配的方法進行優化，由表2可以發現，重建圖像質量和編碼時間都有所提高。

表2 快速稀疏分形算法和本文算法實驗對比Table 2 Experimental comparison between fast sparse fractal algorithm and algorithm of this paper

小波分析自創立以來，在分形理論上的應用就日漸廣泛。小波分形編碼算法利用小波變換后多級子帶之間具有相似性，將全局匹配變為鄰域匹配，提高編碼速度，但是重建圖像的質量并不理想。觀察表3可以看出，在小波分形算法和改進算法的比較中，可以發現本文算法在重建圖像質量上都有非常大的提升，尤其是Lake的峰值信噪比，提高了20.86 dB，是改進算法的近2.63倍，但是在編碼時間上部分略快一些。

將四種算法的峰值信噪比PSNR展示在圖5上。在圖6中只展示其余三種算法。

圖5 四種算法在9幅圖像壓縮的PSNR對比結果Fig.5 PSNR comparison results of four algorithms in 9 images compression

觀察圖5和圖6可以看出，改進算法在重建圖像質量和編碼效率上都有改善。

圖6 三種算法在9幅圖像壓縮的Time對比結果Fig.6 Time comparison results of three algorithms in 9 images compression

這里選取了紋理復雜度不同的四幅圖像，分別是Milkdrop、Woman2、Lake、Goldhill。將四幅圖像在四種算法（基本分形算法[32]、快速稀疏分形算法[33]、小波分形算法[8]、本文算法）下進行仿真實驗，得到圖7，截取其中一個片段放大處理。

圖7 四種算法的重建圖像Fig.7 Reconstructed images of four algorithms

目前圖像壓縮從大體上看已經有相對良好的水平，但是部分細節還需要進一步完善和優化。從仿真實驗得到的圖像中可以看出，紋理相對簡單的圖像對比效果更明顯。根據四種算法得到的圖像1、圖像2、圖像3和圖像4，分別將它們與原始圖像進行對比，觀察圖8可以看見圖像4與原圖最接近，有的地方甚至比原圖更清晰。在圖像4中實物的邊緣部分更光滑，而圖像3與原圖比較最模糊，在放大的圖像3中只能看見大致的實物結構。圖像2僅次于圖像4，比圖像1清晰。通過仿真圖像的對比可以發現本文的算法得到的圖像效果更好。

圖8 四種算法的重建圖像的部分放大圖Fig.8 Partial enlarged view of reconstructed image of four algorithms

本文根據原始圖像的特征分類索引，將字典分成了相應的類別塊，在相應類別塊中利用正交稀疏方法，記錄其與初始化殘差乘積最大值對應的位置，通過迭代得到滿足重構誤差限下的最佳匹配Domain塊的序號，根據序號找到分形圖像的Range塊?？焖傧∈璺中尉幋a算法每次迭代都需要將整個字典與殘差做內積，因此用改進方法計算內積在很大程度上減少編碼時間，并且在相應類別塊中搜索，可以縮小搜索范圍、提高匹配的準確率。圖像Woman2、Lake、Milkdrop的匹配效率分別是91.18%、93.52%、90.43%。根據數據可以發現改進算法的匹配效率均超過了90%，改進算法利用稀疏思想減少信息冗余，在保證高匹配效率的同時，解決了基本分形算法耗時的問題。

4 結束語

針對分形圖像在壓縮后，重建圖像的質量和時間不理想，本文算法利用NMF無監督特征提取的性質，將PNMF與L3/2范數約束相結合，利用稀疏分解來改進分形壓縮算法，提出了雙層非負矩陣分解的分形圖像壓縮算法，提取的圖像特征具有稀疏性和魯棒性。在仿真實驗中，對比四種算法可以發現本文算法在分形圖像編碼時，顯著改善重建圖像的細節效果。今后可以嘗試將本文算法運用在對于重建圖像質量要求較高的分形圖像壓縮領域，如醫療圖像或者存檔掃描的圖像等有價值圖像的保存和傳輸。除了運用本文算法進行特征提取外，還可以嘗試其他方法提高分形圖像重構的匹配搜索效率，進一步縮短編碼時間。