一種面向打印的低維光譜空間構造法

王 瑩, 王忠民, 翟社平

(西安郵電大學 計算機學院， 陜西 西安 710121)

一種面向打印的低維光譜空間構造法

王 瑩, 王忠民, 翟社平

(西安郵電大學 計算機學院， 陜西 西安 710121)

基于Kubelka-Munk混色理論模型，對其光譜吸收散射比空間進行線性修正，對光譜反射率空間進行正態化處理，并在兩種空間之間建立經驗變換，之后利用非線性優化技術和主成分分析方法，得到一個保持打印基色光譜特性的低維空間，將光譜反射率數據變換至該空間，以實現針對彩色打印的多光譜數據降維。實驗結果表明，新方法能使低維空間數據保持高維光譜空間數據的最主要信息，克服傳統主成分分析法導致重建光譜超出光譜數據范圍的缺陷。和主成分分析法相比，新方法的色度精度、光譜精度以及同色異譜指數三個評價指標都有提高。

多光譜圖像；Kubelka-Munk混色理論模型；非線性優化；主成分分析

對場景高質量的圖像采集與輸出日益成為科技和商用追求的目標。隨著多通道相機與超四色打印機的發展和普及，多光譜圖像越來越受到人們的重視。利用多通道相機采集場景圖像，通過光譜反射率重建技術即可獲得多光譜圖像，因此多光譜圖像的數據為源場景的光譜反射率。該類圖像主要用于源場景在不同觀察環境下的精確再現，目前在藝術作品分析與修復[1-2]、醫學[3-4]、軍事目標圖像采集與分析等領域已獲得廣泛應用。

多光譜圖像不同于傳統RGB三通道的色度圖像，它是利用窄帶濾色鏡，在可見光波長范圍內對場景反射光進行采樣，通常在380nm到730nm之間間隔幾納米，因此圖像數據維數高，這導致在進行顯示器顯示、打印機輸出等硬拷貝設備再現時，對圖像的處理難度高、計算耗時長、所需存儲空間大。所以建立一個低維空間，將高維光譜反射率圖像數據降維至該空間，然后對降維后的圖像進行處理，成為多光譜圖像硬拷貝設備再現的重要一環。

目前針對多光譜圖像進行降維的方法主要有非負主成分分析法[5]、WSPCAplus法[6]、主成分分析降維法[7-8]、LabPQR空間法[9]等，這些方法都是在主成分分析(Principal Component Analysis, PCA)法的基礎上，將光譜反射率降維至一個6維空間。但光譜反射率空間為非正態空間，使用主成分分析法導致降維不穩定，低維數據并不能真正表示高維光譜的最主要信息。此外，在使用經典主成分分析法降維時，從低維數據重構高維數據往往導致高維光譜反射率超出范圍[0,1]。當光線照射到物體表面時，反射輻通量與入射輻通量之比就是光譜反射率，因此圖像的光譜反射率是一個比值，且為正，這使不在區間[0,1]的反射率數據失去了物理意義。

本文針對多光譜圖像的打印應用，在Kubelka-Munk混色模型基礎上，提出面向打印的低維光譜空間構造法，用以實現多光譜數據的降維。新方法對Kubelka-Munk模型中的吸收散射系數比空間進行修正，并對源光譜反射率空間進行正態變換，用以提高降維空間的線性性及正態性；并用修正后的空間進行降維變換，以提高多光譜圖像的降維精度。

1 理論概述

1.1 Kubelka-Munk混色理論模型

場景物體受到光線照射后會產生反射與吸收，反射輻通量與吸收輻通量相關[10-11]。

假設在可見光波段范圍內，總的采樣數為h,場景物體在第i(i=1,2,…,h)個采樣波長上對光的光譜吸收系數為ki，對應的光譜散射系數為si，則光譜反射率在各波長上的反射系數可表示為光譜吸收系數ki與散射系數si的比值的函數，且吸收與散射系數比和基色顏料的濃度成近似線性關系。另設混合色在第i個采樣波長上的光譜吸收系數為Ki，光譜散射系數為Si，那么，由相應的基色顏料混合形成混合色的混色公式可表示為

(1)

其中cj表示基色顏料的濃度，n表示基色數。光譜反射率與Φi的轉換關系可由Kubelka-Munk理論得出。

針對不透明介質作變換

其中ri,∞是不透明介質的厚度為無窮大時的光譜反射率在第i個采樣波長上的光譜反射系數，即厚度再增加也不會影響物體的光譜反射率。

材料表面光譜反射率與透明色膜及不透明承印介質底層之間在光學上的相互關系可描述為

其中ri,g為不透明承印介質的光譜反射率在第i個采樣波長上的光譜反射系數，X為色層的厚度。

1.2 主成分分析法

主成分分析法能將高維度數據變換至一個低維特征空間并盡最可能地保持源高維數據的方差，即能用少數變量來表示源數據，同時保持這組數據的最重要信息，是常用的降維方法。

設有由m個h維列矢量形成的矩陣

R=(R1，R2，…，Rm),

其協方差矩陣Σ的秩為u。根據PCA方法，首先求出Σ的所有非零特征值λ1,λ2,…,λu及其對應的特征向量α1,α2,…,αu。利用求得的特征向量形成降維矩陣

V=(α1,α2,…,αu)，

實現源數據的降維，得降維后數據

B=VTR=(β1,β2,…,βu)T，

其中β1稱為第一主成分，能反映出源數據的最主要信息，β2稱為第二主成分，能反映源數據次主要信息，依此類推，βu為第u主成分。

為實現從低維數據重構高維數據,可取

其中β1α1是源數據R的最大近似；β2α2為源數據R的第二大近似；依此類推，βuαu為源數據R的第u近似。

采用主成分分析法降維本質上是使源數據和降維重構數據盡可能地接近[12]，即使得兩者的誤差盡可能地小。

2 低維光譜空間構造法

文[13]為預估構成一幅多光譜圖像的顏料基色的譜密度，從測量得到的多光譜圖像的光譜反射率中，通過PCA法得到三個特征向量，對其進行變換獲得三個非負向量，以此代表顏料基色的譜密度。一旦計算出基色顏料的譜密度，則可以生成光譜圖像的任何譜密度，且其色度誤差和光譜誤差均較小。重建的譜密度進而可以變換為光譜透射率或反射率。由此得到啟示，相對于采用主成分分析法進行降維時以光譜反射率空間作為降維空間，Kubelka-Munk模型中以Φi為特征形成的空間可能更適于作為降維空間，因為光譜吸收散射系數比Φi與基色顏料的濃度成近線性關系，在用有限種基色去合成混合色的應用中，一個線性的高維混色空間更易用一個低維空間去近似。以下將Kubelka-Munk模型中具有m個混合色，且以Φi為特征形成的混色空間記為Φ。

2.1 混色空間與特征向量空間的變換

由PCA法得到的特征向量之間線性無關，這就導致從低維的特征向量空間重構出的高維光譜空間數據有正有負，會出現超出范圍[0,1]的現象，這時重構數據并不能真正反映物理顏料的光譜特性；同時，由特征向量空間變換到高維光譜反射率空間的系數矩陣元素也有正有負，不能用于表示真實顏料的濃度。由于顏料的吸收系數與散射系數在各波長上均為正值，所以Φi=Ki/Si的值也為正值，而基色顏料的單位吸收系數在各波長上也均為正，所以對應各基色濃度必然非負。由Kubelka-Munk混色模型可知，任何打印色均是由不同濃度的基色顏料混合而成，式(1)即為混色公式。

根據式(1)可將Φ表示為h×m矩陣

Φ=φC,

(2)

其中h×u矩陣φ是基色顏料的單位吸收散射系數比的矩陣，u×m矩陣C是重建Φ的濃度矩陣。為推導出特征向量矩陣V與矩陣φ之間的相互變換關系，對表征矩陣Φ施加PCA分析，得到其特征值矩陣B和對應的特征向量矩陣V。由

Φ=VB=φC

看出特征向量矩陣V與矩陣φ之間是近線性關系，可表示為

φ=VBC-1=VN,

(3)

其中C-1為C的廣義逆。若已知φ，則基色濃度矩陣

C=φ-1VB。

(4)

推導出特征向量與顏料基色的光譜特性之間的變換關系后，可采用非線性優化技術獲得基色的光譜特性。

2.2 替代空間的建立

使用PCA法降維，當樣本集為非正態分布時，得到的低維空間并不能代表樣本的光譜特性。由于不確定樣本如何分布，所以PCA法忽略了那些信息貢獻率小，但對重構光譜反射率很重要的特征值所對應的特征向量。若樣本集為多元正態分布，則可保證從樣本集中計算出的特征值與特征向量能用來表征整體采樣空間的特征值與特征向量，這等同于在采樣空間中進行最優采樣，使所采樣本的特征向量張成的空間能表征整體樣本空間。若一個空間中的樣本是多元正態分布，則其中99%的樣本應分布于2維橢圓、3維橢圓體和多維超橢球內。在特征向量坐標系內，橢球的軸向長度與其對應的特征值成正比。若空間中的樣本分布為非橢球體，則降維后的數據并不能很好地重構出源光譜數據。

Kubelka-Munk混色模型在推導時有兩個假設前提:(1) 顏料的膜層在水平方向朝向二維空間蔓延，當光以水平方向通過與之垂直的邊緣時，其能量損失與向上和向下擴展時的損失相比非常小，可忽略不計； (2) 僅考察漫射光在色膜中向上下兩個方向擴展時，其在色膜中傳播時所產生的情況。應用中，往往實際情況與假設前提相違背，這就使得該模型在使用中往往會遇到問題而需要進行改進。

Kubleka-Munk模型本身僅是著色過程的近似，真實的材料是向四面八方散射光，這就背離了Kubleka-Munk模型的假設前提。此外，在采用分光光度計進行光譜測量時，儀器的照明光幾乎為直光，而并非現實環境中的漫射光。對于吸收光較少的厚膜，在光線照射入色膜不太深之前就被全部輻射了，遵從Kubelka-Munk混色模型。但假如色膜很薄以至于光未能進行充足的散射或光尚未發生散射就已經被全部吸收，這將使實際結果與Kubelka-Munk模型的模擬結果有很大不同，致使Φ空間的線性假設不成立，此時基色顏料的濃度與混合色的吸收散射系數比已經不再是線性關系，這些均給式(3)帶來嚴重誤差。

針對Kubelka-Munk混色模型存在的問題，研究者提出了很多辦法來提高混色預測精度，但是這些方法大多較為復雜，在實現時主要依靠優化技術，為使精度達到要求，需設置較多的優化參數。

為了解決Kubelka-Munk混色模型和PCA法自身在使用時所產生的問題，基于Kubelka-Munk混色模型中的Φ空間，建立一個類似的線性光譜吸收散射系數比空間，以該空間替代Φ空間,建立關于基色顏料光譜特性φ的線性變換，將該空間記為Ψ空間。對新空間的推導與對Φ空間的推導一致，不再贅述。針對多光譜降維，需要建立一種針對Ψ空間與光譜反射率空間的變換關系。

(1) 由于樣本正態性是影響光譜降維穩定性及精確度的重要因素，因此該變換要能夠將非正態分布樣本集轉換為正態分布樣本集，而對正態分布樣本集的變換要求仍然能保證樣本的正態分布特性。

(2) 對正態變換后的矩陣進行降維，能夠獲得一個與顏料基色光譜特性相對應的低維空間。

(3) 針對打印輸出，在所建立的低維空間中，標量相乘與矢量相加能夠描述減色成色原理，即打印及印染成色原理。

為滿足上述特性要求，在實驗基礎上總結得到經驗變換

Ψ=a- norm(R),

(5)

其逆變換為

norm(R) =a-Ψ,

(6)

其中Ψ表示新建立的線性空間，該空間數據符合正態分布。a為變換補償矩陣，主要用于減小Kubelka-Munk混色模型的假設前提與實際不符時所產生的線性誤差。a矩陣從樣本集中通過優化方法求得，優化的目標函數為樣本集與通過PCA法得到的低維空間中重構的光譜的均方根誤差。R表示源光譜反射率數據，norm(R)表示對R施加正態變換。

式(5)要滿足前述變換特性，必須具備三個條件。

(1) 補償矩陣a為常量矩陣，norm(R)為對R進行正態變換，結果為正態分布，因此兩者相加結果仍符合正態分布。

(2) 新空間Ψ是吸收散射系數比空間Φ的替代，也是關于φ的線性空間，可將其分解為特征值矩陣和特征向量矩陣的乘積。特征向量空間是一個低維空間，可取其維數與顏料基色的個數相同。

(3) 采用非線性優化技術對Ψ空間降維，將其變換到一個低維空間，其光譜特性與顏料基色的光譜特性一致，因此在利用超四色打印對多光譜圖像輸出時，顏料基色的光譜特性可用Ψ空間降維后形成的低維光譜空間的光譜特性φ近似表示。

2.3 基于替代空間的多光譜圖像降維算法

建立Ψ空間后，可在Ψ空間上實現多光譜圖像的降維。具體算法分為如下3個步驟。

(1) 對于任意一幅多光譜圖像，利用式(5)，首先對圖像的光譜反射率數據進行正態變換，進而變換到Ψ空間。

(2) 通過式(3)，對Ψ空間數據采用主成分分析法，得到特征向量矩陣V及特征值矩陣B，并利用非線性優化技術，得到非負特征向量矩陣φ，作為對顏料基色的光譜特性的近似，并用φ張成低維的光譜空間。

(3) 根據非負特征向量矩陣φ，利用式(4)，將光譜反射率數據變換為打印顏料的基色濃度矩陣，作為光譜反射率的降維數據，從而實現了多光譜圖像面向打印的數據降維。

從低維的濃度數據重構高維光譜反射率數據的過程分為以下3個步驟。

(1) 利用式(2)，計算得到Ψ空間數據。

(2) 利用式(6)，得到正態分布的光譜反射率數據。

(3) 對上述數據實施正態變換的逆運算，獲得重建的光譜反射率數據。

3 實驗驗證

為驗證算法的有效性，實驗中對超四色打印機HP designjet 130nr的色空間進行采樣，生成實驗樣本集。130nr為C(青)、M(洋紅)、Y(黃)、K(黑)、c(淺青)、m(淺洋紅)6通道打印機；光譜測量設備為GretagMacbeth SpectroScanT分光光度計，測量值為光譜反射率，該設備在380 nm到730 nm波長范圍內每隔10 nm進行采樣，可取得36維光譜反射率，實驗選擇400 nm到700 nm譜段內的31維光譜反射率數據。實驗中的樣本集是在130nr的 CMYKcm空間中，每個色軸上均勻劃分5級，將各個色軸上的分級數據進行組合，然后對組合得到的色塊進行打印，采用分光光度計進行測量，得到15 626個樣本。實驗中，非正態分布的樣本集到正態分布樣本集的變換采用Box-cox變換族中的平方根變換。

實驗針對超四色打印色空間中的15 626個樣本，分別采用PCA法及新算法進行仿真。對降維結果分別從色度、光譜及同色異譜指數三方面予以評價。

對色度精度的評價，需要將多光譜圖像變換至色度圖像，在此不妨采用CIELAB均勻顏色空間的標準色差公式ΔEab進行評價，其定義為

其中L、a、b分別表示CIELAB空間色度值的三個分量，可由光譜反射率與光照的光譜功率分布及觀察者色匹配函數積分計算獲得。

光譜精度評價采用均方根誤差ERMS公式，定義為

其中N為高維光譜空間的維度，Δβ(λi)為在波長λi處，由低維空間重構的高維光譜反射率與源光譜反射率的誤差。

同色異譜指數是指同色異譜對在改變光照情況下的色差大小，實驗中將光照從D65變換為A，采用ΔEab進行衡量。

表1給出采用兩種方法，將光譜反射率樣本降維至6維空間的平均誤差統計結果。由表1可知，新算法在三個評價指標上明顯優于經典主成分分析降維法。

表1 PCA法與新算法的實驗結果比較

圖1給出實驗樣本集中的10個樣本采用兩種方法的降維效果比較。圖1(a)為源光譜反射率曲線，圖1(b)為采用PCA法降維后重構的光譜反射率曲線，圖1(c)為采用新算法降維后重構的光譜反射率曲線。從圖1可以看出，采用PCA法重構的光譜反射率，有超出范圍[0, 1]的現象，這使得重構光譜失去了光譜反射率的物理意義。此外，對于波長大于650 nm的光譜，PCA法對源光譜的模擬明顯變差，而新算法基于Kubelka-Munk理論模型，利用光譜吸收散射比空間的線性特性及光譜吸收系數和散射系數均為正數的特點，使重構光譜數據保持在 范圍內，可用于表征光譜反射率的真實特性。對于波長大于650 nm的光譜，通過新算法重構的光譜，對源光譜的模擬遠遠優于PCA法。

(a) 源光譜

(b) PCA法重構光譜

(c) 新算法重構光譜

圖1 源光譜反射率與重構光譜反射率對比

4 總結

針對多光譜圖像維度高導致圖像打印輸出困難的問題，提出一種構造低維光譜空間的方法，并將圖像的光譜反射率變換至該空間，實現了多光譜圖像的降維。

新方法通過采用Kubelka-Munk混色理論模型，并結合主成分分析法，利用光譜吸收散射比空間為顏料基色濃度的線性變換這一特點，用光譜吸收散射比空間代替光譜反射率空間作為降維空間，提高了降維過程的線性性。

針對Kubelka-Munk模型的局限性，建立了一個新的線性空間。新空間通過引入補償矩陣，顯著縮小了PCA及Kubelka-Munk模型實際應用時所帶來的誤差。

在將光譜反射率變換至新空間時，首先對光譜反射率進行正態變換，保證了線性降維的穩定性，使得低維空間數據能保持高維光譜空間數據最主要信息。

此外，新方法克服了傳統主成分分析法導致重建光譜超出光譜數據范圍的問題。與主成分分析法相比，新算法有效降低了多光譜圖像降維的誤差，從而提高了彩色打印的精度。

[1] Sitnik R, Krzeslowski J, Maczkowski G. Archiving shape and appearance of cultural heritage objects using structured light projection and multispectral imaging[J]. Optical Engineering, 2012, 51(2): 021115-021123.

[2] Berns S R, Chen Tongbo. Practical total appearance imaging of paintings[C]//Proceedings of IS&T. Denmark Copenhagen: Wiley, 2012: 162-167.

[3] Bouzid M, Khalfallah A, Bouchot A, et al. Automatic cell nuclei detection: a protocol to acquire multispectral images and to compare results between color and multispectral images[C]//Proceedings of SPIE, Imaging, Manipulation, and Analysis of Biomolecules, Cells, and Tissues. USA California San Francisco:SPIE, 2013: 85871-85881.

[4] Jakovels D, Kuzmina I, Berzina A, et al. Noncontact monitoring of vascular lesion phototherapy efficiency by RGB multispectral imaging[J]. Journal of Biomedical Optics, 2013, 18(12): 126019-126026.

[5] 王瑩, 曾平, 羅雪梅, 等. 一種用于低維光譜空間構造的非負主成分分析法[J]. 四川大學學報, 2010, 42(2):165-170.

[6] 王瑩, 王忠民, 王義峰, 等. 面向色彩再現的多光譜圖像非線性降維方法[J]. 光學精密工程, 2011, 19(5):1171-1178.

[7] Wu Chongyi, Lee Shunming, Wen Chanhua, et al. Multi-spectral image acquisition system for color spectrum reproduction[C]//Proceedings of CVGIP. Taiwan Kinmen:Mingchuan University Press, 2003: 115-122.

[8] Bakke M A, Farup I, Hardedberg Y J. Multispectral gamut mapping and visualization: a first attempt[C]//Proceedings of SPIE. USA California San Jose: SPIE, 2005, 5667: 193-200.

[9] Derhak W M, Rosen R M. Spectral colorimetry using LabPQR: An interim connection space[J]. Journal of IS&T, 2006, 50(1): 53-63.

[10] Kubela P, Munk F. Ein Beitrag zur Optik der Farbanstriche[J]. Zeitschrift fur Technische Physik (German), 1931,12(5): 593-601.

[11] Kubleka P. New Contribution to the Optics of Intensely Light-Scattering Materials[J]. Journal of the Optical Society of America, 1948, 38(5): 448-457.

[12] 王小剛, 田小平, 吳成茂. 基于壓縮感知的圖像快速重構去噪算法[J]. 西安郵電學院學報,2012, 17(4):11-14.

[13] Ohta N. Estimating Absorption Bands of Component Dyes by Means of Principal Component Analysis[J]. Analytical Chemistry, 1973,45(2): 553-557.

[責任編輯:瑞金]

Construction of low-dimensional spectral space for multi-ink printing

WANG Ying, WANG Zhongmin, ZHAI Sheping

(School of Computer Science and Technology, Xi’an University of Posts and Telecommunications, Xi’an 710121, China)

By correcting coefficient ratio space of the absorption based on scattering of Kubelka-Munk turbid media theory, and exerting a normalization transformation to the spectral reflectance space, a new linear space is created. Then a transformation between the spectral reflectance space and the new space is established. A low dimension spectral space which is consistent with the multi-ink color space is finally achieved by utilizing the nonlinear optimization and principal component analysis. The dimensionality reduction of the multi-spectral image for multi-ink printing is accomplished through the transformation of spectral reflectance to this low dimension space. Experiments show that the data obtained by this new method in the low dimension space can keep the main spectral information of the data in the high dimension spectral space. Furthermore, the new method can overcome the shortcoming brought by the principal component analysis, which is that the spectral reflectance reconstructed from the low dimension space may exceed the range of the spectral data. Compared with the principal component analysis method, this new method can make great improvements in the colorimetric accuracy, spectral accuracy and metamerism index.

multispectral image, Kubelka-Munk turbid media theory, nonlinear optimization, principal component analysis

10.13682/j.issn.2095-6533.2014.06.016

2014-06-05

陜西省自然科學基金資助項目(2012JM8044)；陜西省教育廳科學研究計劃資助項目(12JK0733)；西安郵電大學青年教師科研基金資助項目(ZL2011-09)

王瑩(1977-)，女，博士，講師，從事顏色科學及多光譜圖像處理研究。E-mail:wangyingjsj@xupt.edu.cn 王忠民(1967-)，男，博士，教授，從事智能信息處理研究。E-mail:zhmwang@xupt.edu.cn

TP391.4

A

2095-6533(2014)06-0080-06