何頌華，陳 橋*，段 江

1. 深圳職業(yè)技術(shù)學(xué)院傳播工程學(xué)院，廣東 深圳 518055

2. 西南財經(jīng)大學(xué)經(jīng)濟信息工程學(xué)院，四川 成都 610075

基于有約束非負(fù)矩陣分解的原稿基色色料光譜預(yù)測方法

何頌華1，陳 橋1*，段 江2

1. 深圳職業(yè)技術(shù)學(xué)院傳播工程學(xué)院，廣東 深圳 518055

2. 西南財經(jīng)大學(xué)經(jīng)濟信息工程學(xué)院，四川 成都 610075

針對直接在光譜反射率空間，對原稿顏色樣本光譜的主成分分析會導(dǎo)致特征向量的數(shù)目超過真實物理維度(原稿所用基色色料)的數(shù)量，以及特征向量和對應(yīng)系數(shù)存在負(fù)值，不能直接表示原稿基色色料的光譜特性和對應(yīng)濃度等情況。創(chuàng)新性的提出需根據(jù)原稿色料的光學(xué)特性建立一個完全線性的光譜空間，并在該空間中使用帶約束條件的非負(fù)矩陣分解實現(xiàn)對原稿基色數(shù)量和光譜形狀進(jìn)行預(yù)測的方法。對此，首先設(shè)計了一個實現(xiàn)對原稿基色色料光譜預(yù)測方法的總體研究方案和實現(xiàn)步驟，再以透明色料原稿為例，研究如何選擇和構(gòu)建一個符合其光學(xué)特性的光譜線性空間，然后再在基本非負(fù)矩陣分解(BNMF)基礎(chǔ)上提出針對基色色料光譜預(yù)測的有約束非負(fù)矩陣分解算法(SCNMF)。針對BNMF算法會出現(xiàn)多重最優(yōu)解，為了提高預(yù)測精度，使矩陣分解結(jié)果有明確的物理意義，所提出的SCNMF算法需要滿足4個約束條件：非負(fù)性約束；全加性約束；平滑性約束；稀疏性約束。建立了滿足約束條件的目標(biāo)函數(shù)和迭代算法。預(yù)測結(jié)果表明本文提出的新方法能有效的實現(xiàn)對原稿基色物理維度和基色色料光譜的準(zhǔn)確預(yù)測。

基色色料光譜預(yù)測；光譜顏色復(fù)制；非負(fù)矩陣非解；線性混合空間；主成分分析

引 言

光譜顏色復(fù)制技術(shù)使用光譜數(shù)據(jù)編碼即光譜反射率或透射率數(shù)據(jù)來描述顏色信息，因此可使同色異譜程度達(dá)到最小[1]。并且光譜顏色復(fù)制技術(shù)在對顏色光譜進(jìn)行重建時必須采用多基色成色方式，這就大大增加了顏色疊加的自由度，擴大了顏色再現(xiàn)的范圍，因此特別適用于高保真顏色復(fù)制[2-3]。要實現(xiàn)對原稿的光譜顏色復(fù)制，首先必需確定一組基色色料，以在多光譜空間中進(jìn)行原稿色彩的表示與再現(xiàn)。這組基色色料的成功確定必須包括兩個方面：對基色物理維度的估計，即確定原稿顏色是由幾種基色色料合成得到；對基色色料光譜的估計。目前，用于基色色料光譜預(yù)測的方法不多，主要有基本主成分分析法PCA[4-8]，以及Tzeng提出的基于旋轉(zhuǎn)主成分分析的光譜預(yù)測方法[9]。將PCA法用于光譜特征分析時，會出現(xiàn)特征向量和其對應(yīng)權(quán)重系數(shù)出現(xiàn)負(fù)值，因而在物理上沒有任何意義的情況；基于旋轉(zhuǎn)主成分分析的光譜估計方法預(yù)測精度取決于旋轉(zhuǎn)初值的選取，具有很大的不確定性和不準(zhǔn)確性。本研究將結(jié)合非負(fù)矩陣分解提出一種新的基色色料預(yù)測技術(shù)，可較為精確的實現(xiàn)對透明色料原稿基色物理維度和色料光譜的估計。

1 原稿基色色料光譜預(yù)測方法的總體研究方案設(shè)計

基于有約束非負(fù)矩陣分解的原稿基色色料光譜預(yù)測方法將完全摒棄由主成分分析所建立的特征向量，而是根據(jù)非負(fù)矩陣分解局部特征構(gòu)成整體的思想，直接對非負(fù)的光譜數(shù)據(jù)集進(jìn)行矩陣分解，并進(jìn)行非負(fù)性約束和其他條件約束，分解得到的結(jié)果基向量不存在負(fù)值，在實際應(yīng)用中具有明確物理意義，并具有純加性的特點。由于在光譜反射率空間，主成分分析得到的顯著維度超出了實際物理維度[9]，因此該方法必須首先要將光譜反射率空間轉(zhuǎn)換到一個能反映實際物理維度的線性混合空間，并在這個線性混合空間中進(jìn)行有約束非負(fù)矩陣分解。

基于有約束非負(fù)矩陣分解的原稿基色色料光譜預(yù)測方法的總體研究方案如圖1所示。

在該研究方案中最為關(guān)鍵的三步為色料光譜線性模型建立、針對原稿色料光學(xué)特性建立線性混合空間并在該空間進(jìn)行有約束非負(fù)矩陣分解預(yù)測原稿色料光譜。為了對色料光譜預(yù)測結(jié)果進(jìn)行驗證，在本研究方案中還設(shè)計了預(yù)測光譜與實際基色色料光譜的評價方法。

圖1 基于有約束非負(fù)矩陣分解的原稿基色色料光譜預(yù)測方法的總體研究方案

Fig.1 Spectral prediction research program for original primary pigments based on constraint NMF

2 色料線性混合模型和線性混合空間建立

不同原稿其色料混合的光學(xué)特性不一樣，因此在建立原稿色料線性混合空間時應(yīng)充分考慮原稿色料的光學(xué)特性。本研究以透明色料原稿為例來分析其色料性線混合模型和線性混合空間的建立。kubleka-munk混介理論給了一個很好的啟發(fā)，即可以通過各色料吸收系數(shù)和散射系數(shù)與色料濃度加權(quán)后的線性相加得到混合后的吸收系數(shù)和散射系數(shù)。但大量的研究文獻(xiàn)已表明對于不透明和半透明色料，其色料濃度與吸收系數(shù)和散射系數(shù)并不是一種線性關(guān)系[10]。

對于透明色料原稿來說，由于色料完全透明，并附著在一高反射基底上的原稿，如相紙，其各色劑層的散射系數(shù)s近似為0，只有吸收系數(shù)k。因此透明色料相對于不透明和半透明色料其光學(xué)特性能較多的符合Kubleka-Munk模型的假設(shè)條件，則多種透明色料混合后的吸收系數(shù)可以由參與混合的各色料濃度線性表示，即滿足

(1)

且在式(1)中各色料濃度之和為1，這樣對于透明色料，光譜吸收系數(shù)空間就是一個完全線性的混合空間。

3 用于基色色料光譜預(yù)測的約束非負(fù)矩陣分解算法(SCNMF)

在確定線性混合空間后，就可以直接利用所提出的光譜約束非負(fù)矩陣算法(spectra constrained non-negative matrix factorization algorithm，SCNMF)進(jìn)行基色數(shù)量和基色光譜的預(yù)測。為了算法的通用性，這里用Ψλ表示線性的光譜空間。

3.1 約束條件的確定

由于本研究方案是直接對線性空間Ψλ的光譜數(shù)據(jù)進(jìn)行矩陣分解以估計基色色料光譜和(濃度)混合百分比，為使分解的結(jié)果有明確的物理意義，首先必須保證矩陣分解時滿足約束條件：

(1)非負(fù)性約束，基色色料光譜和(濃度)混合百分比不可能為負(fù)值；

(2)全加性約束，本改進(jìn)算法由于基于式(1)所建立的色料混合線性模型，該線性混合模型對各色料(濃度)混合百分比除了非負(fù)的約束外還有全加性約束，即要求任意像素其混合色料的各濃度百分比之和為1。

利用最優(yōu)化技術(shù)可歸結(jié)為求解下面的問題

(2)

(3)

(4)

(5)

其中，p為Ψλ空間進(jìn)行主成分分析得到的顯著維數(shù)，反映參與混合的基色色料數(shù)量。Ψλ, p為預(yù)測的光譜值，Ψλ為需要驗證目標(biāo)的光譜值，F(xiàn)是最小化目標(biāo)函數(shù)。式(4)和式(5)為約束條件。

上面的最優(yōu)化問題正好與非負(fù)矩陣分解算法思想相一致，對比于基本非負(fù)矩陣分解算法，其約束優(yōu)化條件多個一個全加性約束，該約束條件可在非負(fù)矩陣分解的迭代過程中加入。采用文獻(xiàn)[11]所使用的表示形式，具體做法是引入新的矩陣ψs和Ψs

(6)

其中1T表示一個向量的所有分量都為1，即1=(1，1，…，1)T,δ為權(quán)值，控制總和為1這一約束條件對目標(biāo)函數(shù)的影響。用ψs和Ψs替代ψ和Ψ進(jìn)行算法的迭代，并在每次迭代完ψ之后，從ψs和Ψs中取出原本屬于ψ和Ψ的部分，并重新進(jìn)行一次式(6)的擴展即可。這樣在算法收斂時同時也能進(jìn)行全加性和為一的約束。

對于基色色料光譜的物理意義除了滿足約束條件(1)和(2)外，這里還應(yīng)該考慮到色料光譜自身特點和色料分布特點，需要增加約束條件(3)和(4)：

(3)平滑性約束(smoothness constraint)，色料來自于自然界或人工提取，其光譜本身具有平滑性，曲線形狀變化平緩，一般不存在窄帶光譜；

(4)稀疏性約束(sparseness constraint)，所謂稀疏性是指因為在一般情況下，任何基色色料的(濃度)混合百分比分布都不會充滿整個圖像空間，混合色料中包含的只是幾種基色色料，而不是所有基色色料。即每種基色色料的混合百分比只是分布于空間中的部分區(qū)域，具有一定的稀疏度。對于非負(fù)矩陣分解來說，其中一個矩陣的稀疏性必然導(dǎo)致另一個矩陣是非稀疏性(或稱平滑性)，使得最終的乘積能夠更好的擬合觀測信號。這一點恰好加強了約束條件(3)色料光譜的平滑性特征。

3.2 SCNMF目標(biāo)函數(shù)和迭代公式

由于平滑性約束和稀疏性約束是針對非負(fù)基色色料光譜矩陣ψ和濃度百分比矩陣C本身數(shù)據(jù)分布的約束，一般都是通過對非負(fù)矩陣目標(biāo)函數(shù)引入懲罰項來約束矩陣ψ和C[12]。在非負(fù)矩陣分解的兩種目標(biāo)函數(shù)中，歐氏距離最為直觀，應(yīng)用最為廣泛，由于色料光譜預(yù)測的物理意義與歐氏距離比較符合，故以歐氏距離的情況下進(jìn)行算法推導(dǎo)。加入懲罰項后，目標(biāo)函數(shù)由式(7)變?yōu)槭?8)的形式

(7)

(8)

其中，‖‖F(xiàn)代表矩陣的Frobenius范數(shù)，當(dāng)且僅當(dāng)R=WH時，函數(shù)取得最小值。式(8)中，aJ1(ψ)和βJ2(C)是依賴于應(yīng)用的懲罰項，a和β是規(guī)則化參數(shù)，它們平衡著逼近誤差和約束之間的消長關(guān)系。對于色料光譜矩陣的平滑性約束和濃度百分比矩陣的稀疏性約束，引入懲罰項見式(9)

(9)

這樣修改的目標(biāo)函數(shù)為

(10)

在進(jìn)行算法迭代時，采用基本非負(fù)矩陣算法(basic non-negative matrix factorization algorithm，BNMF)的乘性迭代規(guī)則，目標(biāo)函數(shù)式(10)是單調(diào)非增的。為滿足全加性限制，迭代時分別用ψe和Ψe替代ψ和Ψ，其迭代公式為

(11)

(12)

ψe和C的梯度為

(13)

(14)

將式(13)和式(14)代入式(11)和式(12)中，得

(15)

(16)

定義迭代步長

(17)

(18)

將式(17)和式(18)代入式(15)和式(16)之后可得到乘性迭代式(19)和式(20)

(19)

(20)

由于BNMF的解不具有唯一性，即它存在多重最優(yōu)解，這顯然與色料光譜預(yù)測的唯一性相矛盾，為了解決解的非唯一性問題，可通過合理的構(gòu)建目標(biāo)函數(shù)，對分解結(jié)果增加除了非負(fù)性外的其他更多條件限制，因而提出的SCNMF算法針對基色色料光譜的自身特點和分布特點，在解的結(jié)果為非負(fù)性的基礎(chǔ)上增加了全加性約束、平滑性約束、稀疏性約束，將從理論上更容易找到全局最優(yōu)解，增強基色色料光譜預(yù)測的準(zhǔn)確性。

4 結(jié)果與討論

實驗中將對kodak Q60C的IT8.7/2反射色標(biāo)(圖2)進(jìn)行實驗驗證，該色標(biāo)介質(zhì)是相紙，相紙是由一高反射不透明基底加上三種透明染料層(C，M，Y)制得的，因此已知各顏色塊呈色是通過三種透明染料合成得到。首先使用Gretag spectraeye測量色標(biāo)的光譜反射率，并將其轉(zhuǎn)換到吸收系數(shù)K空間，再在這兩個空間分別進(jìn)行PCA分析。

圖2 kodak Q60C的IT8.7/2反射色標(biāo)

表1是分別在光譜反射率空間和吸收系數(shù)空間用PCA得到的前六個特征向量重構(gòu)光譜的貢獻(xiàn)率比較?？擅黠@看出，在光譜反射率空間使用6個特征向量進(jìn)行光譜重構(gòu)時，累積貢獻(xiàn)率達(dá)到99.85%，而在吸收系數(shù)空間只需要使用3個特征向量進(jìn)行光譜重構(gòu)，累積貢獻(xiàn)率就達(dá)到99.95%，在使用5個特征向量重構(gòu)時，累積貢獻(xiàn)率接近100%。因為累積貢獻(xiàn)率反映了多個特征向量對信息的貢獻(xiàn)程度，將累積貢獻(xiàn)率增加到99%以上并且增長速度不再明顯變化時的特征向量數(shù)目確定為基色染料混合的物理維度，因此在吸收系數(shù)空間中可確定染料混合的物理維度為3，與實際情況相吻合。

表1 光譜反射率空間和吸收系數(shù)空間特征向量貢獻(xiàn)率比較

Table 1 The percent variance in both reflectance and absorption space of IT8.7/2 refflection target

第i個特征根光譜反射率R空間吸收系數(shù)K空間特征值貢獻(xiàn)率累積貢獻(xiàn)率特征值貢獻(xiàn)率累積貢獻(xiàn)率11 54600 80890 80899 61450 79100 791020 25860 13710 94601 64400 13530 926330 09190 04870 99470 88930 07320 999540 00470 00250 99720 00310 00030 999750 00230 00120 99850 00190 00020 999960 00200 00100 99950 00080 00010 9999

接下來再分析在兩個空間進(jìn)行光譜重構(gòu)時光譜重構(gòu)精度和色度精度的表現(xiàn)。評價光譜重構(gòu)精度使用均方根誤差sRMS，評價色度精度使用CIELAB ΔE*ab色差公式，計算色差時采用D50標(biāo)準(zhǔn)光源和CIE1931 2度視場的顏色匹配函數(shù)。同樣發(fā)現(xiàn)，在吸收系數(shù)空間，通過三個特征向量重構(gòu)后的光譜精度和色差都遠(yuǎn)好于在光譜反射率空間通過五個特征向量重構(gòu)的光譜。

從表2和表3中都可看出，在吸收系數(shù)空間只需要三個特征向量進(jìn)行光譜重構(gòu)后，其sRMS與反射率空間使用六個特征向量的sRMS接近，但其色差的最小值、最大值、平均值都要好于在反射率空間使用六個特征向量進(jìn)行光譜重構(gòu)后表現(xiàn)。

表2 在光譜反射率空間和吸收系數(shù)空間進(jìn)行光譜重構(gòu)后sRMS比較

表3 在光譜反射率空間和吸收系數(shù)空間進(jìn)行光譜重構(gòu)后色差比較

從兩個空間的特征向量累積貢獻(xiàn)率、光譜精度和色度精度的比較分析可得出，將光譜值從反射率空間轉(zhuǎn)化到吸收系數(shù)空間后，需要的特征向量將大為減少。而已知kodak Q60C的IT8.7/2色標(biāo)本身是用相紙制得的，因此各顏色塊呈色是通過三種染料合成得到，其色料混合的線性空間當(dāng)然是三維的，但如果單從光譜反射率空間通過分析特征向量的方式尋找物理色料的數(shù)目，表1證明其結(jié)果與實際不符合，但如果將光譜值從反射率空間轉(zhuǎn)到吸收系數(shù)空間，會發(fā)現(xiàn)進(jìn)行重構(gòu)的特征向量的數(shù)量正好與實際物理色料的數(shù)量一致。

由以上分析可知對于透明色料連續(xù)調(diào)原稿，吸收系數(shù)空間可作為色料混合的線性空間，并通過主成份分析確定物理色料數(shù)為3，從而明確了在吸收系數(shù)空間中進(jìn)行SCNMF算法需要分解的基向量個數(shù)為3。圖3為在吸收系數(shù)空間中采用SCNMF算法對基色色料吸收光譜的估計，其預(yù)測結(jié)果與該相紙實際所用的色料光譜曲線形狀非常接近。在使用SCNMF算法進(jìn)行平滑性限制時α=0.1，β=0，最后的量化結(jié)果是通過取6次實驗的平均值得到的。

圖3 吸收系數(shù)空間中對kodak Q60C的IT8.7/2 色標(biāo)透明色料光譜的預(yù)測

Table 4 Similarity degree between estimated spectrum and real spectrum of IT8.7/2 in absorption space

透明染料光譜曲線相似度BNMF估計SCNMF估計青0 99880 9993品0 99810 9987黃0 99840 9992

由于預(yù)測光譜與實際光譜在比較時進(jìn)行了歸一化，因此不能通過均方根方根誤差等評價方式判斷預(yù)測精度，這里本研究中通過矢量相關(guān)分析比較兩者之間的相似度。相關(guān)系數(shù)越接近于1，預(yù)測光譜與實際光譜越相似。對相關(guān)系數(shù)的計算可利用matlab的corrcoef函數(shù)，三種透明染料在吸收系數(shù)空間中預(yù)測光譜與實際光譜的曲線相似度見表4。表中的數(shù)據(jù)比較也表明通過SCNMF算法預(yù)測的結(jié)果要好于BNMF算法預(yù)測的結(jié)果。

5 結(jié) 論

針對本研究提出的有約束非負(fù)矩陣分解的原稿基色色料光譜預(yù)測方法得出以下結(jié)論：

(1)不同原稿其色料混合的光學(xué)特性不一樣，因此在建立原稿色料線性混合空間時應(yīng)充分考慮原稿色料的光學(xué)特性。對于透明色料連續(xù)調(diào)原稿選用吸收系數(shù)空間作為線性混合空間，該空間是完全線性的。

(2)在對色料混合的物理維度判斷上，可將累積貢獻(xiàn)率增加到99%以上并且增長速度不再明顯變化時的特征向量數(shù)目確定為基色染料混合的物理維度，但該判斷的前提條件是所建立的線性混合空間是一個能真實反映色料混合的完全線性空間。

(3)為了提高預(yù)測精度，使矩陣分解結(jié)果有明確的物理意義，所提出的SCNMF算法需要滿足4個約束條件：非負(fù)性約束、全加性約束、平滑性約束、稀疏性約束。

(4)在評價預(yù)測光譜的預(yù)測精度時，由于對預(yù)測光譜進(jìn)行了歸一化，因此不能通過均方根方根誤差等評價方式判斷預(yù)測精度，可通過矢量相關(guān)分析比較預(yù)測光譜與實際光譜之間的相似度。相關(guān)系數(shù)越接近于1，預(yù)測光譜與實際光譜越相似。

[1] HE Song-hua, CHEN Qiao, DUAN Jiang(何頌華, 陳 橋, 段 江). Spectroscopy and Spectral Analysis(光譜學(xué)與光譜分析), 2015, 35(6): 1459.

[2] KONG Ling-wang, ZHU Yuan-hong, LI Qiong(孔令罔, 朱元泓, 李 瓊, 等). Geomatics and Information Science of Wuhan University(武漢大學(xué)學(xué)報·信息科學(xué)版), 2006, 31(9): 788.

[3] WANG Ying, ZENG Ping, WANG Yi-feng(王 瑩，曾 平，王義峰). Acta Optica Sinica(光學(xué)學(xué)報), 2009, 29(8): 2122.

[4] LI Jin-cheng, LIU Zhen, CHEN Guang-xue(李金城, 劉 真, 陳廣學(xué)). Acta Optica Sinica(光學(xué)學(xué)報), 2009, 29(8): 2122.

[5] Xu Dongbo, Wang Xiangzhao, Bu Yang. Chinese Optics Letters, 2012, 10(12): 121202-1.

[6] XU Fa-qiang, WAN Xiao-xia, ZHU Yuan-hong(許法強, 萬曉霞, 朱元泓). Optics and Precision Engineering(光學(xué) 精密工程), 2008, 16(3): 518.

[7] HE Song-hua, LIU Zhen, CHEN Qiao(何頌華, 劉 真, 陳 橋). Acta Optica Sinica(光學(xué)學(xué)報), 2014, 34(2): 0233001-1.

[8] Wang Ying, Zeng Ping, Luo Xuemei, et al. Journal of Southeast University·English Edition, 2013, 2(4): 486.

[9] Tzeng D Y, Berns R S. Color Research and Application, 2005, 30(2): 84.

[10] He Songhua, Liu Zhen. Procedia Engineering, 2011, 23: 320.

[11] Heinz D C, Chang C I. IEEE Geoscience and Remote Sensing Letters, 2001, 39(3): 529.

[12] Lee D, Seung H. Nature, 1999, 401: 788.

(Received Aug. 3, 2015; accepted Nov. 16, 2015)

*Corresponding author

The Spectral Prediction of Original Primary Pigment Based on Constrained Non-Negative Matrix Factorization

HE Song-hua1，CHEN Qiao1*，DUAN Jiang2

1. School of Communication, Shenzhen Polytechnic, Shenzhen 518055, China

2. Computer Science Department of Southwestern University of Finance and Economics, Chengdu 610075, China

With direct prediction in the spectral reflectance space with principal component analysis, the numbers of eigenvectors will surpass the numbers of real primary pigments while the eigenvectors and the corresponding coefficients have negative value, which can not directly presented original primary pigment spectral characteristics and corresponding concentration. We proposed an innovative spectral prediction method in which a complete linear spectral space was created according to optical properties of originals pigment. A constrained non-negative matrix factorization algorithm to predict the numbers and spectral curve shapes of real primary pigments was used in the space. So, this paper designed an overall research plan and implementation process about spectral prediction method firstly, and studied how to select and establish a spectral linear space which was conformed to optical properties of originals; taking transparent pigments as example, and spectra constrained non-negative matrix factorization (SCNMF) algorithm was established to predict primary pigment spectra based on basic non-negative matrix factorization algorithm (BNMF). Aiming at realizing multiple optimal solution of BNMF and improving the prediction accuracy as well as make the matrix decomposition results to be clearly physically meaningful; the proposed SCNMF needs to satisfy four constraints: non negative constraint, additive constraint, smoothness constraint and sparseness constraint. The objective function and iterative algorithm to meet four constraints were set up. The prediction results show that the proposed method can realize accurate prediction of original primary pigments’ numbers and spectra effectively.

Spectral estimation of primary pigments; Spectral color reproduction; Non-negative matrix factorization; Linear mixing space; Principal component analysis

2015-08-03，

2015-11-16

國家自然科學(xué)基金項目(61108087)，四川省杰出青年學(xué)術(shù)技術(shù)帶頭人資助計劃項目(2012jq0017)資助

何頌華，1977年生，深圳職業(yè)技術(shù)學(xué)院教授 e-mail: hdh1818@szpt.edu.cn *通訊聯(lián)系人 e-mail: qiaochen@szpt.edu.cn

TP391.4; O432.3

A

10.3964/j.issn.1000-0593(2016)10-3274-06