陳 永，陶美風(fēng)，陳 錦

(1.蘭州交通大學(xué) 電子與信息工程學(xué)院，蘭州 730070；2.甘肅省人工智能與圖形圖像處理工程研究中心，蘭州 730070)

敦煌莫高窟是世界上現(xiàn)存規(guī)模最宏大、內(nèi)容最豐富的佛教石窟壁畫寶庫，其內(nèi)所存的壁畫、經(jīng)卷等具有珍貴的研究價值。然而，由于自然風(fēng)化的破壞及人為因素的影響，窟內(nèi)壁畫出現(xiàn)了地仗脫落、劃痕、褪色、裂紋等嚴重的災(zāi)害，亟待保護[1]。因此，研究病態(tài)敦煌壁畫的修復(fù)極其重要。但人工修復(fù)存在風(fēng)險大、耗時長且修復(fù)不可逆等問題，將數(shù)字化虛擬修復(fù)應(yīng)用于古代壁畫的保護，是近年來國內(nèi)外文化遺產(chǎn)保護的研究熱點[2]。

目前，傳統(tǒng)圖像修復(fù)方法主要分為3類：基于偏微分方程、基于樣本塊和基于稀疏表示的修復(fù)方法。基于偏微分方程的主要修復(fù)模型有TV(total variation)模型[3]和CDD(curvature-driven diffusion)模型[4]，該類方法可擴散完成小區(qū)域破損圖像的修復(fù)[5]?；跇颖緣K的修復(fù)方法以Criminisi算法[6]為代表，該類算法可以完成較大區(qū)域破損圖像的修復(fù)，但易出現(xiàn)填充錯誤的問題。陶兆勝等[7]在Criminisi算法的基礎(chǔ)上，提出了基于邊緣特征和像素結(jié)構(gòu)相似度的圖像修復(fù)算法，減少了像素誤匹配率。

基于稀疏表示的修復(fù)方法利用字典中極少量原子及其系數(shù)的線性組合，完成對破損圖像的重構(gòu)[8]。Zhang等[9]提出了組稀疏圖像修復(fù)方法，但該方法未考慮圖像全局信息，會導(dǎo)致修復(fù)結(jié)果出現(xiàn)結(jié)構(gòu)傳播錯誤和線條斷裂的現(xiàn)象。Zha等[10]提出將圖像塊與組稀疏聯(lián)合稀疏表示的方法，但該方法也僅對局部固定大小區(qū)域進行聯(lián)合塊組字典建立。Ghorai等[11]提出一種考慮局部圖像塊統(tǒng)計信息和幾何特征的多尺度金字塔稀疏表示修復(fù)方法。Zhao等[12]通過采用RPCA分解的方法，提出了一種結(jié)合閾值選擇和中值濾波器的圖像去噪算法。Qiang等[13]通過形態(tài)成分分析(morphological component analysis,MCA)方法，將圖像分解為結(jié)構(gòu)層和紋理層的修復(fù)方法，并應(yīng)用于云南劍川石寶山壁畫修復(fù)中。Hu等[14]提出一種基于圖像分解的唐卡圖像修復(fù)方法，后采用模糊集和局部二值模式完成修復(fù)。然而，上述稀疏表示算法字典設(shè)計單一，僅對局部固定大小區(qū)域構(gòu)建字典，導(dǎo)致壁畫圖像修復(fù)結(jié)果易出現(xiàn)結(jié)構(gòu)不連貫和模糊塊效應(yīng)等問題，此外在圖像分解修復(fù)時，還存在顏色紋理光學(xué)屬性分離不徹底的問題[15-16]。

針對上述問題，本文提出了一種基于塊核范數(shù)的RPCA分解與熵權(quán)類稀疏的壁畫修復(fù)方法。主要工作有：1)首先，采用提出的基于塊核范數(shù)的RPCA圖像分解算法，通過引入凸先驗信息塊核范數(shù)，利用塊核范數(shù)進行紋理矯正操作，將待修復(fù)圖像分解為結(jié)構(gòu)層和紋理層，更精確地完成結(jié)構(gòu)和紋理信息的分離。2)然后，對結(jié)構(gòu)層采用提出的熵權(quán)類稀疏表示進行修復(fù)，通過熵加權(quán)k-means算法對相似圖像塊進行聚類，構(gòu)建稀疏子類字典，完成結(jié)構(gòu)層圖像的重構(gòu)，克服了傳統(tǒng)稀疏表示字典單一的問題。3)最后，對紋理層采用雙三次插值進行修復(fù)，融合得到修復(fù)后的圖像。通過對敦煌壁畫的修復(fù)實驗結(jié)果表明，本文方法較對比算法，獲得了較好的主客觀評價效果。

1 相關(guān)理論

稀疏表示就是用極少量原子和系數(shù)的線性組合對信號進行重構(gòu)[17]。稀疏表示的關(guān)鍵在于尋找給定信號的最佳稀疏解，其數(shù)學(xué)定義為

(1)

式中：α為稀疏系數(shù)，T為待觀測信號，D為字典，‖·‖0為l0范數(shù)，用來衡量α中非零元素個數(shù)。一般自然信號存在噪聲或重構(gòu)圖像會出現(xiàn)一定的誤差ε，因此，為了平衡圖像的稀疏性和重構(gòu)誤差，式(1)可以改寫為

(2)

在稀疏表示的圖像修復(fù)中，圖像退化的過程可定義為

Y=JX+n

(3)

式中：Y為退化后的圖像，X為未退化的原始圖像，J為退化算子，n為高斯白噪聲。因此，圖像修復(fù)過程就是對退化圖像Y進行推測估計得到原始圖像X的過程。

根據(jù)式(2)，在稀疏字典D給定的條件下，可得到稀疏表示圖像修復(fù)的正則化模型為

(4)

通過式(4)求解稀疏系數(shù)α，最后通過Y=Dα得到修復(fù)后的圖像。

2 本文算法

2.1 基于塊核范數(shù)的RPCA圖像分解

一般圖像中既包含結(jié)構(gòu)信息，也包含紋理信息，其中結(jié)構(gòu)信息主要是圖像邊緣和輪廓等變化劇烈的信息部分；紋理信息主要是圖像中亮度或灰度值變化緩慢或周期性變化區(qū)域，具有重復(fù)性、規(guī)則性和方向性，其反映了圖像中同質(zhì)現(xiàn)象的一種視覺特征[15,18]。RPCA[19-20]是一種高維數(shù)據(jù)降維分析方法，通過RPCA可以將圖像M(M∈Rm×n)分解為兩個矩陣的線性相加，一個是結(jié)構(gòu)矩陣L(L∈Rm×n，rank(L)≤m,n)，主要用來表征圖像的結(jié)構(gòu)信息；另一個是紋理矩陣S(S∈Rm×n)，主要是圖像的紋理細節(jié)信息，其圖像分解模型定義為

(5)

式中：‖·‖0表示l0范數(shù)，γ為平衡圖像的結(jié)構(gòu)和紋理分解的因子。

以圖1敦煌莫高窟第92窟“銜花雙鹿”的局部壁畫為例進行RPCA分解說明。其中，圖1(a)為原始壁畫，圖1(b)和1(c)分別為RPCA分解后的結(jié)構(gòu)層和紋理層，圖1(b)結(jié)構(gòu)層中出現(xiàn)了模糊現(xiàn)象，圖1(c)紋理層也出現(xiàn)了分離不徹底的現(xiàn)象。從圖1實驗可以發(fā)現(xiàn)：采用RPCA分解后紋理層中仍含有部分結(jié)構(gòu)信息。這是因為采用RPCA空間變換分解后，破壞了紋理層的低秩性，導(dǎo)致結(jié)構(gòu)信息仍殘留在紋理層[21]。

圖1 RPCA算法分解結(jié)果Fig.1 Decomposition results of RPCA algorithm

針對以上RPCA圖像分解不徹底的問題，本文提出了一種基于塊核范數(shù)的RPCA圖像分解方法，以更好地分解結(jié)構(gòu)與紋理信息。以圖2紋理圖像為例，圖像塊2(a)和2(b)表現(xiàn)出全局不同但局部相似的特點，此時不滿足全局重復(fù)、相似的低秩性要求。因此，為了將紋理層恢復(fù)到低秩性，通過對局部紋理圖像塊進行矯正操作，使其恢復(fù)為低秩矩陣。低秩矩陣的恢復(fù)是一個非凸連續(xù)函數(shù)優(yōu)化問題，是一個NP難題，可以將求解最小秩的問題轉(zhuǎn)化為核函數(shù)最小化的問題，即核函數(shù)為秩的凸包絡(luò)，是矩陣秩的凸近似[16]。

圖2 紋理的局部相似及全局不相似性示意Fig.2 Schematic of local similarity and global dissimilarity of texture

塊核范數(shù)‖H‖BNN是紋理部分的局部圖像塊，經(jīng)過方向矯正和重疊取塊操作后形成的矩陣奇異值之和，其定義為

(6)

圖3 塊核范數(shù)示意Fig.3 Schematic of block nuclear norm

將塊核范數(shù)引入到RPCA圖像分解模型，則得到如下的分解模型

(7)

為了驗證本文分解方法的有效性，以圖1所選取壁畫為例進行實驗，并與RPCA分解結(jié)果比較，見圖4。其中，圖4(a)、(b)分別為RPCA分解和本文方法得到的結(jié)構(gòu)層圖像，比較發(fā)現(xiàn)：本文方法可以使結(jié)構(gòu)層信息更加清晰；圖4(c)、(d)為紋理層分解結(jié)果，發(fā)現(xiàn)通過本文方法分離的紋理和結(jié)構(gòu)信息更加徹底。

圖4 引入塊核范數(shù)的RPCA分解結(jié)果比較Fig.4 Comparison of RPCA decomposition results with block nuclear norm

2.2 熵加權(quán)k-means結(jié)構(gòu)層圖像聚類

對結(jié)構(gòu)層圖像進行類稀疏修復(fù)時，在k-means聚類算法的基礎(chǔ)上，對歐氏距離進行熵加權(quán)改進作為新的聚類準則，提出了熵加權(quán)的k-means聚類算法，從而得到結(jié)構(gòu)相似的子類圖像。待分類圖像塊A的信息熵E(A)為

(8)

式中：PA(i)為圖像塊A中RGB顏色通道中第i個等級的像素占比率；i表示像素等級，范圍為0～255。定義圖像塊A與聚類中心B的信息熵特征差異性系數(shù)為

qj=|E(A)-E(B)|

(9)

式中0≤qj≤1，qj越小，圖像塊A與B越相似。然后，采用式(10)得到圖像塊A的熵權(quán)重vj

(10)

從而得到圖像塊A與B的熵權(quán)歐氏距離計算準則

(11)

通過式(11)將待分類圖像塊A劃分到熵權(quán)歐氏距離最小的類中，從而構(gòu)造得到結(jié)構(gòu)層的子類圖像。為了驗證熵加權(quán)k-means聚類算法對壁畫聚類的有效性，以圖4(b)結(jié)構(gòu)圖像為例進行實驗，結(jié)果見圖5。其中，圖5(a)為原始圖像，圖5(b)為聚類結(jié)果，圖5(c)、(d)、(e)為各子類圖像，可以發(fā)現(xiàn)，具有相似結(jié)構(gòu)的圖像塊能夠較好地進行聚類，得到的結(jié)構(gòu)層圖像聚類更為準確。

2.3 類稀疏表示的結(jié)構(gòu)層圖像修復(fù)

在結(jié)構(gòu)層子類圖像劃分的基礎(chǔ)上，對各子類進行類稀疏修復(fù)，定義類稀疏修復(fù)模型為

(12)

式中：yk是第k類結(jié)構(gòu)層圖像塊；Dk和αk分別是對應(yīng)的字典和稀疏系數(shù)。然后，采用奇異值分解和分裂Bregman迭代優(yōu)化算法更新求解字典和系數(shù)[9]。將圖像塊yk的估計值rk進行奇異值分解，即

(13)

式中：Uk、Vk分別為yk的左奇異正交矩陣和右奇異正交矩陣,∑k為對角矩陣，則字典Dk中每個原子dki為

dki=ukivki,i=1,2,…,t

(14)

式中uki、vki分別為Uk、Vk的列向量，t為第k類字典的原子數(shù)。所以，第k類自適應(yīng)字典Dk為

Dk={dk1,dk2,dk3,…,dki},i=1,2,…,t

(15)

以圖5結(jié)構(gòu)層圖像聚類結(jié)果為例，采用奇異值分解學(xué)習(xí)更新得到自適應(yīng)子類字典，見圖6。其中，圖6(a)、(b)、(c)分別為圖5子類對應(yīng)的自適應(yīng)字典?？梢园l(fā)現(xiàn)，本文提出的類稀疏表示方法得到的各子類字典更能夠體現(xiàn)圖像的結(jié)構(gòu)特征，且聚類字典訓(xùn)練更加高效。

圖5 熵加權(quán)的k-means算法聚類結(jié)果Fig.5 Clustering results of entropy weighted k-means algorithm

圖6 自適應(yīng)字典Fig.6 Adaptive dictionaries

在完成子類字典的更新后，接著對子類圖像進行類稀疏系數(shù)求解

(16)

采用分裂Bregman迭代優(yōu)化算法對式(16)求解，從而完成破損圖像塊的稀疏編碼。最后，利用各子類圖像的字典和系數(shù)完成結(jié)構(gòu)層各子類圖像的修復(fù)，并按照索引矩陣放回原始結(jié)構(gòu)層圖像，即得到結(jié)構(gòu)層圖像的修復(fù)結(jié)果YS

(17)

2.4 分層圖像融合重構(gòu)

采用雙三次插值(bicubic interpolation)對紋理層破損區(qū)域進行插值修復(fù)，完成圖像紋理細節(jié)信息的重建。雙三次線性插值是通過對待修復(fù)點最近的16個采樣點的加權(quán)平均完成像素點的估計，按下式計算

(18)

式中：f(p,q)為待修復(fù)點坐標，W(p-pi)和W(q-qj)分別為待修復(fù)點的橫坐標權(quán)重和縱坐標權(quán)重。權(quán)重的計算如下

(19)

最后將結(jié)構(gòu)層的修復(fù)結(jié)果YS與紋理層的修復(fù)結(jié)果YT融合，得到最終的修復(fù)圖像Y

Y=YS+YT

(20)

2.5 算法步驟

Step1輸入待修復(fù)壁畫圖像，并獲得其掩膜圖像。

Step2將待修復(fù)壁畫圖像利用塊核范數(shù)的RPCA方法進行分解，得到結(jié)構(gòu)層和紋理層。

Step3對結(jié)構(gòu)層圖像采用熵加權(quán)的k-means算法進行聚類，對每個子類圖像構(gòu)造類稀疏字典，并采用奇異值分解和分裂Bregman進行類稀疏表示修復(fù)。

Step4對紋理層圖像進行雙三次插值修復(fù)。

Step5將結(jié)構(gòu)層和紋理層的修復(fù)結(jié)果融合操作，得到修復(fù)后的壁畫圖像。

3 實驗結(jié)果與分析

實驗運行軟件環(huán)境為Windows 10操作系統(tǒng)，采用MATLAB R2016a軟件，硬件配置為Intel(R)Core i7-9700K CPU@3.6 GHz，16.0 GB RAM,NVIDIA GeForce GTX 1660。修復(fù)結(jié)果采用主觀和客觀兩種方式進行評價，客觀評價使用峰值信噪比(peak signal to noise ratio,PSNR)、結(jié)構(gòu)相似性(structural similarity index measurement,SSIM)和算法時間復(fù)雜度進行評價。此外，為了驗證本文方法的有效性，采用人為添加破損和真實破損的敦煌壁畫進行修復(fù)實驗比較，并與經(jīng)典修復(fù)CDD算法[4]、Criminisi算法[6]、文獻[7]和文獻[9]的修復(fù)結(jié)果進行對比分析。

3.1 人為添加破損壁畫修復(fù)

考慮到壁畫修復(fù)是一個病態(tài)逆問題,為了定量分析修復(fù)算法的有效性，一般對完好壁畫圖像人工添加破損后進行定量分析。首先對壁畫進行人為添加劃痕修復(fù)對比實驗，如圖7所示，從上往下分別為“北壁彌勒經(jīng)變”局部壁畫、“雙笛合奏”局部壁畫和“觀無量壽經(jīng)變”局部壁畫。其中，第1列為原始圖像，第2列為掩膜圖像，第3列為CDD算法修復(fù)結(jié)果，從圖7(c)可以發(fā)現(xiàn)，CDD算法存在修復(fù)不徹底的問題，如第二幅壁畫左側(cè)矩形框內(nèi)出現(xiàn)了明顯的修復(fù)痕跡。第4列為Criminisi算法修復(fù)結(jié)果，出現(xiàn)了塊匹配錯誤和結(jié)構(gòu)傳播錯誤的現(xiàn)象，如圖7(d)第三幅壁畫的左側(cè)眉毛部分和右側(cè)矩形框中均出現(xiàn)了誤匹配現(xiàn)象。第5列為文獻[7]修復(fù)結(jié)果，圖7(e)第三幅壁畫的左側(cè)眉毛的連續(xù)性較好，但第一幅壁畫的修復(fù)結(jié)果中仍存在像素匹配錯誤的問題。第6列為文獻[9]修復(fù)結(jié)果，其修復(fù)結(jié)果較符合人類的視覺感受，但因為該算法僅考慮了圖像局部相關(guān)性，導(dǎo)致修復(fù)結(jié)果中出現(xiàn)了線條斷裂和模糊的現(xiàn)象，如圖7(f)第三幅壁畫的眉毛部分出現(xiàn)了線條不連貫的問題。最后一列圖7(g)為本文修復(fù)結(jié)果，與其他算法比較，本文算法修復(fù)后圖像更加連貫自然，獲得了更好的視覺效果。

為了進一步對圖7的修復(fù)結(jié)果做出定量評價，采用PSNR和SSIM進行比較，見表1。PSNR和SSIM是衡量圖像質(zhì)量的重要定量評價指標，其中PSNR值越大，代表失真越少，即修復(fù)效果越好；SSIM是一種衡量兩幅圖像相似度的指標，該值越大，表明修復(fù)后圖像的失真程度越小[18]。從表1中可以發(fā)現(xiàn)，本文算法在PSNR和SSIM上均優(yōu)于其他對比算法，表明本文方法修復(fù)質(zhì)量更好，從而驗證了本文方法對于人工破損壁畫修復(fù)的有效性。

圖7 不同算法對人為添加破損壁畫的修復(fù)結(jié)果對比Fig.7 Comparison of mural inpainting results of artificial damages by different algorithms

表1 不同算法修復(fù)結(jié)果PSNR和SSIM對比Tab.1 Comparison of PSNR and SSIM results of different algorithms

表2 不同算法修復(fù)時間復(fù)雜度對比Tab.2 Comparison of inpainting time complexity of different algorithms

通過表2可以發(fā)現(xiàn)：CDD算法時間復(fù)雜度最低，但CDD算法修復(fù)效果較差，無法徹底完成修復(fù)；文獻[9]組稀疏時間復(fù)雜度最高，修復(fù)時間耗時最長；Criminisi算法和文獻[7]時間復(fù)雜度均低于文獻[9]，但遠高于本文算法，這是因為文獻[7]中果蠅優(yōu)化迭代j值以及Criminisi算法中匹配塊搜索值r遠大于本文劃分類k值。本文算法時間復(fù)雜度僅高于CDD算法，遠小于其他3種比較算法，說明本文在達到較好修復(fù)效果的同時，算法執(zhí)行效率也較高。

3.2 真實破損壁畫修復(fù)

為了進一步驗證本文算法的有效性，采用4組真實破損壁畫圖像進行修復(fù)實驗，見圖8。其中，圖8(a)為真實壁畫圖像，從左往右分別為“第263窟·人字坡·供養(yǎng)菩薩”、“第61窟·東壁維摩詰經(jīng)變·騎隊”、“第285窟·北壁四龕楣·供養(yǎng)菩薩”和“第220窟·飛天”的局部壁畫，圖8(b)為掩膜圖像。圖8(c)為CDD修復(fù)結(jié)果，從第一幅和第三幅可以看出該算法出現(xiàn)了結(jié)構(gòu)模糊的問題。圖8(d)為Criminisi修復(fù)結(jié)果，第二幅和第三幅壁畫右側(cè)矩形框內(nèi)的供養(yǎng)菩薩頭光部分出現(xiàn)了結(jié)構(gòu)傳播錯誤問題，第四幅壁畫出現(xiàn)了塊匹配錯誤。圖8(e)為文獻[7]修復(fù)結(jié)果，同樣在第三幅和第四幅壁畫中存在塊匹配錯誤；圖8(f)為文獻[9]組稀疏修復(fù)結(jié)果，從第一幅壁畫圖像和第四幅壁畫矩形框內(nèi)箭頭所指區(qū)域存在修復(fù)痕跡，第二幅和第三幅壁畫出現(xiàn)了線條斷裂現(xiàn)象。圖8(g)為本文結(jié)果，可以看出，在線條連續(xù)性和清晰性方面本文修復(fù)效果均有所提高。

圖8 不同算法對真實破損壁畫的修復(fù)效果對比Fig.8 Comparison of mural inpainting results of real damages by different algorithms

4 結(jié) 論

1)提出了一種基于塊核范數(shù)的RPCA圖像分解方法，將塊核范數(shù)引入RPCA圖像分解模型，對RPCA分解后的紋理層通過矯正操作，克服了圖像結(jié)構(gòu)層與紋理層分離不徹底的問題。

2)提出了熵權(quán)類稀疏的修復(fù)方法，采用熵加權(quán)的k-means聚類算法得到結(jié)構(gòu)相似的子類圖像，并通過奇異值分解和分裂Bregman的類稀疏表示完成結(jié)構(gòu)層圖像的重構(gòu)，克服了單一字典學(xué)習(xí)出現(xiàn)的結(jié)構(gòu)傳播錯誤和模糊現(xiàn)象。

3)破損敦煌壁畫修復(fù)實驗表明，本文方法取得了較好的視覺效果和較高的定量評價，修復(fù)效率也進一步提高。

4)盡管本文方法具有較好的修復(fù)效果，但未對壁畫圖像的歷史、宗教文化等語義信息進行考慮，后續(xù)將采用語義分析等深度學(xué)習(xí)方法進一步研究。