随着图像处理技术的不断发展，多光谱图像越来越多地应用于颜色信息的再现、复制及匹配当中，由于传统的颜色图像再现依赖于同色异谱，并取决于特定观察条件和观察者，一旦这些条件改变，同色异谱现象将被破坏，为了更准确地得到颜色信息，可以采用多光谱图像重建来获取物体表面的光谱反射率，通过光谱反射率映射出物体的颜色信息，从而得到真正的颜色再现。

多光谱反射率重建是多光谱成像技术研究的重要内容，主要包括两个方面：光谱反射率重建算法的研究和训练样本方法的研究。国内外学者对光谱反射率重建的研究做了深入工作，提出了很多重建算法，如线性成像系统下的伪逆法^[1]、维纳估计法^{[2, 3]}、有限维模型法^[4]、非线性成像系统下的多项式回归模型法和有限制的多项式回归模型法^[5]等。训练样本的选取方法主要包括Hardeberg提出的基于最小条件数选取方法^[6]、Cheung与Westland提出的4种基于颜色样本相互间矩最大化的选取方法^[7]、Shen等提出的基于特征向量优化的颜色样本选取方法^[8]，以及Mohammadi等提出的基于聚类的颜色样本选取方法^{[9, 10]}。他们选择训练样本的共同目的都是如何从众多的颜色样本中选出有代表性的颜色，同时保证色彩的广泛性和均匀性。但有的是基于光谱反射率空间挑选的，有的是基于色度空间挑选的，均只考虑了光谱反射率空间或者色度空间。

为了进一步提高光谱重建精度，本文考虑到光谱反射率空间和色度空间的综合特性，提出了一种基于核模糊C聚类的样本优化选取方法。首先基于光谱反射率空间来选取聚类初始点，再基于色度空间进行模糊C聚类，使同一聚类内的样本之间具有较高的颜色相似性。实验结果表明，这种优化样本选取方法综合考虑了光谱反射率空间的广泛性和色度空间的相似性，提高了光谱重建反射率精度和色度精度。在实验中使用11通道宽带多光谱成像系统对本文方法和已有4种方法的性能进行了对比。

1 光谱成像系统模型及光谱重建算法 1.1 多光谱成像原理

通常，多光谱成像系统是由一组滤光片和单色CCD相机组成，采用多光谱成像技术来重建物体的光谱反射率需要得到第i个滤光片和相机的数字响应g_i，对应的响应输出模型为：

$\begin{align} & {{g}_{i}}=\int\limits_{{{\lambda }_{max}}}^{{{\lambda }_{min}}}{l\left( \lambda \right)r\left( \lambda \right){{f}_{i}}\left( \lambda \right)o\left( \lambda \right)s\left( \lambda \right)d\lambda +{{b}_{i}}+{{n}_{i}}} \\ & =\int\limits_{{{\lambda }_{max}}}^{{{\lambda }_{min}}}{{{q}_{i}}\left( \lambda \right)r\left( \lambda \right)d\lambda +{{b}_{i}}+{{n}_{i}}} \\ \end{align}$

(1)

式中，l(λ)、f_i(λ)、o(λ)和s(λ)都为未知量，可将q_i(λ)=l(λ)f_i(λ)o(λ)s(λ)合并为单一的光谱函数，b_i、n_i分别表示第i个滤光片下的暗电流响应和图像噪声。将380～780 nm的可见光波段的光谱反射率通过10 nm间隔采样点表示，可以得到一个N×1维的向量R。相应地，相机通道响应为C×1维向量用g表示。从而(1)可以改写成向量矩阵形式：

$g=QR+b+n$

(2)

式中，Q为C×N的光谱响应函数矩阵，令G=g-b将(2)可以改写为G=QR+n，可以通过盖住镜头的响应来去除暗电流b，校正公式如下：

$G=\frac{{{G}_{data}}-{{G}_{dark}}}{{{G}_{white}}-{{G}_{dark}}}$

(3)

式中，G_data为实际响应，G_dark为盖住镜头响应，G_white为白板响应，若不考虑噪声的影响，其转换公式的基本形式为：

$g=QR$

(4)

1.2 主成分分析(PCA)光谱重建

在实际应用中，任何一个光谱反射率数据集R都可以被近似表示为k个彼此正交的单位基向量b_i的线性组合^[11-13]，即：

$R=Ba=\sum\limits_{i=1}^{k}{{{b}_{i}}{{a}_{i}}}$

(5)

式中，B=[b₁,b₂,…,b_k]为相应的特征向量矩阵，a=[a₁,a₂,…,a_k]^T为相应的转换矩阵。

通常用主成分分析法来计算基向量。假设N*M矩阵R=(R₁,R₂,…,R_m)^T中包含M个不同的光谱反射率，对其进行奇异值分解：

$R=US{{V}^{T}}$

(6)

其中，S是由r=rank(R)个由大到小排列的奇异值组成的对角矩阵。

$S=diag({{\sigma }_{1}},{{\sigma }_{2}},\ldots ,{{\sigma }_{i}})$

(7)

N*N单位正交矩阵U由RR^T的特征向量组成，M*M单位矩阵V由R^TR的特征向量组成。

如果选取前k个主成分来表示原光反射率数据的线性组合，可以用主成分贡献率来表示前k个主成分的重要程度，以确定主成分的个数，前k个特征向量的主成分累积贡献率可以通过式(8)定义。

${{\rho }_{k}}=\sum\limits_{i=1}^{k}{{{\sigma }_{i}}}/\sum\limits_{i=1}^{N}{{{\sigma }_{i}}}$

(8)

本试验中取得光谱反射率的前k个累计贡献率为ρ_k>99.95%，因此可以确定主成分个数k=6。这里将式(5)代入式(4)得到：

${{g}_{0}}=QB{{a}_{0}}~$

(9)

式中，Q、B分别是前k个主成分的响应刺激值，通过式(9)得到系数转换矩阵为a₀=[(QB)^T(QB)]^-1(QB)^Tg₀，若Q和B都是可以根据训练样本求出的已知量，则QB是随着数码响应值变化而变化的，对于新的测试样本，其系数转换矩阵可以表示为：

$a={{a}_{0}}g_{0}^{T}{{\left( {{g}_{0}}g_{0}^{T} \right)}^{-1}}g$

(10)

因此，重建的光谱反射率可以表示为：

$\hat{R}=B{{a}_{0}}g_{0}^{T}{{\left( {{g}_{0}}g_{0}^{T} \right)}^{-1}}g$

(11)

2 训练样本优化选取方法

在反射率重建的过程中，训练样本的选取是影响反射率重建精度的一个重要的因素，上述选取样本的方法都是单一的基于反射率空间或者基于色度空间的选取，并没有综合考虑两者特性的结合，选取出的训练样本不能很好地代表整个颜色样本集合。对于一个众多色样的色序系统，如劳尔色卡，有几百个颜色样本，如果直接作为训练样本会产生过饱和现象，同时影响重建精度。有研究证明^[14]，当训练样本达到一定合适数量时，可以用该数量的训练样本作为标准颜色来重建其他样本的反射率。由于实际系统中噪声的影响，训练样本个数增加到一定程度时，光谱重建精度的提升不再显著，但较少的样本数又不能够代表整个色序系统，因此，本文提出了核模糊C均值聚类算法(KFCM)，来优化训练样本选取方法，达到选取合适数量的样本并能够代表反射率空间和色度空间的特性的目的。

为了从众多颜色样本中选取有代表性的颜色，即保证范围广泛且要求均匀，达到少而精的结果。首先，采用MAXSUM法在光谱反射率空间中选取一定数量的样本数L=(l₁,l₂,…,l_c)，将L个样本转换到色度坐标(a^*,b^*)中作为C个分类初始点；其次，在CIELAB色度空间按照某种相似性进行核模糊C聚类。

2.1 Mercer核函数

在CIELAB色度空间中，L代表明度坐标，a、b代表色度坐标，这里忽略L(亮度)对色度空间样本选取的影响，本文只利用a、b两个维度进行聚类，从根本上避免了由于光源照射不均匀所引起的明度值变化的情况，由于(a^*,b^*)是一个二维空间，在进行分类的过程中，并不一定是线性可分的。为了能够将训练样本正确分类，可以通过核函数内核诱导距离度量将原始的二维空间映射到三维特征空间，找到一个合适的划分超平面，使得不同类中的数据能够更容易地进行清晰的划分，基于此，将映射后的特征空间数据按照某种相似性进行聚类，从而达到更好的聚类效果。假设总的颜色样本数据集有P个样本点l_p,p=(1,2,…,P)，其中l_p∈L_ab，设$\phi $是一非线性映射函数，使色度空间L_ab被此函数映射到高维特征空间H得到$\phi $(l₁)，$\phi $(l₂),…,$\phi $(l_p)输入空间的点积形式在特征空间可以用Mercer核表示：

$~K\left( {{l}_{i}},{{l}_{j}} \right)=\phi \left( {{l}_{i}} \right)\cdot \phi ({{l}_{j}})~$

(12)

目前，对于核函数的选择还没有一个通用的标准，常用的Mercer核函数有多项式核函数、高斯核函数和两层神经网络核函数。高斯核函数所对应的特征空间是无穷的，有限样本在该特征空间中是线性可分的，因此本文中的核函数采用高斯核函数，如式(13)所示：

$K({{l}_{i}},{{l}_{j}})=exp(-\frac{\|{{l}_{i}}-{{l}_{j}}\|}{2{{\sigma }^{2}}})$

(13)

式中，K(l_i,l_j)表示高斯核函数，σ>0是自定义参数。

2.2 FKCM聚类方法

使用高斯核函数将输入色度空间样本映射到特征空间$\phi $(l₁),$\phi $(l₂),…,$\phi $(l_p)，定义特征空间中向量间的欧氏距离为：

$\begin{align} & {{d}_{H}}({{l}_{i}},{{l}_{j}})={{\left\| \varphi \left( {{l}_{i}} \right)-\varphi ({{l}_{j}}) \right\|}^{2}} \\ & \sqrt{=\varphi \left( {{l}_{i}} \right)\cdot \varphi ({{l}_{i}})-2\varphi \left( {{l}_{i}} \right)\cdot \varphi ({{l}_{j}})+\varphi \left( {{l}_{j}} \right)\cdot \varphi ({{l}_{j}})} \\ \end{align}$

(14)

由于核函数是选定的，相似度测量公式相应地可以写成式(15)形式：

${{d}_{H}}({{l}_{i}},{{l}_{j}})=\sqrt{K({{l}_{i}},{{l}_{i}})-2K({{l}_{i}},{{l}_{j}})+K({{l}_{j}},{{l}_{j}})}$

(15)

设数据集U={l_i,i=1,…,P}，u_ij表示样本l_i对类C_i(1,…,C)的隶属度。KFCM算法的目标函数定义为式(16):

${{Q}_{KFCM}}=\sum\limits_{i=1}^{C}{\sum\limits_{j=1}^{P}{u_{ij}^{2}d_{ij}^{2}}}$

(16)

满足式(17)条件：

$\sum\limits_{i=1}^{C}{{{u}_{ij}}}=1$

(17)

C为样本聚类数，d_ij²=‖$\phi $(l_j)-φ(v_i)‖²，v_i是第i类的聚类中心，通过求目标函数Q的最小值，得到KFCM算法的隶属度u_ik和聚类中心v_i的迭代公式(18)：

$\begin{align} & {{u}_{ik}}=\frac{{{(1-K({{l}_{k}},{{v}_{i}}))}^{-1}}}{\sum\limits_{j=1}^{C}{{{(1-K({{l}_{i}},{{v}_{j}}))}^{-1}}}} \\ & {{v}_{i}}=\frac{\sum\limits_{k=1}^{P}{u_{ik}^{2}\cdot K\left( {{l}_{k}},{{v}_{i}} \right)\cdot {{l}_{k}}}}{\sum\limits_{k=1}^{P}{u_{ik}^{2}\cdot K({{l}_{k}},{{v}_{i}})}} \\ \end{align}$

(18)

由以上可得KFCM方法的具体步骤如下:

1. 设定聚类初始点及聚类数C=50，最大迭代次数N=100，高斯核参数σ=250，目标函数值改变量的最小阀值e<1e-6。

2. 将色度空间的(a^*,b^*)值映射到特征空间得到$\phi $(l₁),$\phi $(l₂),…,$\phi $(l_p)。

3. 用模糊C均值聚类得到隶属度矩阵U和类中心矩阵V。

4. 用式(18)依次计算u_ik、v_i。

5. 根据式(16)计算目标函数，如果它小于某个确定的阀值，则算法停止。

6. 若步骤(5)的条件没有满足，则重复上述计算步骤直到满足终止条件。

7 .聚类完成后，选出距离C个聚类中心最近的C个样本数据。

3 实验与讨论

为了验证本文提出的方法的有效性和适用性，采用CIE标准照明体D65光源和多光谱成像系统进行性能考察。实验中，多光谱成像系统由Ocean Optics公司的SpectroCam VIS CCD单色相机(见图 1a)和11个不同波峰中心值(波峰分别为420 nm、460 nm、500 nm、540 nm、580 nm、620 nm、660 nm、680 nm、700 nm、740 nm、780 nm)的宽带滤光片组成，其带宽为20 nm以及滤光片透射率曲线(见图 1b)。实验色卡采用德国RAL(劳尔)色卡，包括213种颜色，通过多光谱成像采集系统拍摄得到RAL色卡样本集在11个波段通道下的图像数据，其光谱反射率由Ocean Optics SpectroSuite分光光度计测得，光谱波长范围380~780 nm，采样间隔10 nm，分别采用Hardeberg方法、Cheung方法、Mohammadi方法、Fvecter方法和本文提出的KFCM聚类优化法进行训练样本的选取。

对于上述5种样本选取方法的光谱重建效果的评价，分别从色度精度和光谱精度两方面进行，色度精度采用不同光源下(CIE标准照明体D65和A，CIE1964标准色度观察者)的CIE2000标准色差公式ΔE_ab进行评价，光谱精度采用光谱均方根误差RMS、适应度系数GFC和光谱匹配偏度指数ISSD进行评价。为了保持和已有样本选取方法的一致，并不失一般性，各种方法选取了同样多的训练样本个数，另外，在本文提出的聚类优化样本选取中，为了确定训练样本个数L，选择出最有代表性的样本，本文得到了5种方法所选取的训练样本个数与重建的平均光谱误差之间的变化关系图，如图 2所示。

从图 2可以看出，本文方法所选训练样本个数在L=50处重建的平均光谱误差达到最小，之后样本个数增加，光谱误差略有升高，分析原因可能是后续选择的部分样本带有较大噪声的影响，总体来讲，5种方法所选取训练样本个数重建的光谱误差分布总趋势均先随着L增大而减小，并在L=50时减小最为明显，随后逐渐趋于平缓。由此可知，在本文中训练样本个数L=50时，对剩余的样本进行光谱重建能够得到较好的重建精度。

图 3为Fvecter法和核模糊C聚类优化法所选训练样本集的色度空间、主成分坐标空间分布图，从图中样本集分布可以看出，Fvecter法所选训练样本相对比较集中，如在色度空间中多集中于黑色区域，所选样本与重建样本的相似性较小，主成分空间中所选样本同样比较集中，与重建样本相关性较大的样本过多，在光谱重建中容易造成过拟合；采用核模糊C聚类优化法所选训练样本的代表性和相关性较强，在色度空间和主成分空间中所选样本都分布比较均匀，且与重建样本具有较强的相似性和连续性，整体也具有较强的代表性和相关性。

表 1为5种样本选取方法的CIE2000色差分析，数据表明，在不同光源(D65和A)下，本文提出的核模糊聚类优化法的最大色差为8.9963，最小色差为0.4698，平均色差为2.9102，较其他4种方法的色差有了提高，且在不同的光照环境下都具有较好的色彩一致性。

表 2为5种样本选取方法重建光谱的光谱误差、适应度系数和光谱匹配偏指数分析。数据表明，本文方法的光谱误差和光谱匹配偏指数的平均值、最大值、最小值均小于其他4种方法，适应度系数的平均值、最大值、最小值均大于其他4种方法。同时，本文方法的GFC平均值在99.01%，非常接近99.5%，表明该方法总体的光谱重建是完全可以接受的，GFC最大值达到了99.9%，达到了非常完美的程度。并且，本文方法的光谱匹配偏度值ISSD的最小绝对值在0.0001，几乎达到了两光谱向量完全相同的程度。

图 4为在本实验宽带光谱成像系统下，使用50个训练样本，通过已有样本选取方法以及本文方法对3个样本的重建结果，其中(a)为8007色样，代表了重建光谱曲线中与实际光谱曲线吻合度较高的样本；(b)为5019色样，代表了重建光谱曲线中与实际光谱曲线吻合度一般的样本；(c)为6038色样，代表了重建光谱曲线中与实际光谱曲线吻合度较差的样本。

结合图 4(a)、图 4(b)和图 4(c)可以看出，本文方法与其他4种方法的重建光谱反射率在低反射率的情况下的对比并不明显，但在高反射率情况下是明显的，该现象可能是恰好所选的3个样本对低反射率重建效果比较好的原因，虽然本文方法的重建光谱反射率与实际测得光谱反射率之间仍存在一些偏差，但是优于其他4种方法。

4 结论

本文针对已有样本选取方法的光谱反射率重建存在的问题，提出了一种核模糊C聚类优化选取方法：首先基于光谱反射率空间选取一定数量的样本初始点，再基于色度空间进行聚类划分，选择离聚类中心最近的样本作为光谱反射率重建的训练样本。这样优化选取的样本既代表了光谱反射率空间的广泛性又体现了颜色空间的相似性。实验结果表明，本文方法在进行光谱反射率重建时具有较高的色度精度与光谱精度，其精度优于已有方法，能在较大程度上满足工业生产中对颜色测量、颜色质量评估等高精度需求。