张 飞,刘 珩,杨 洁,渠立永,吕 琪,黄 伟,李 浩
(1. 61606部队, 北京 100093; 2.陆军工程大学, 南京 210007)
光谱反射率重建技术在颜色科学、文物保护、农业、医学、艺术、纺织等领域具有重要的科学应用价值[1-8]。光谱反射率是指地物在某波段反射通量与入射通量之比[9],是表征物质特有属性的主要特征之一。通过目标样本周围环境信息、光照条件、相机光谱特征以及输出的多通道图像可获得物体的光谱反射率[10]。但颜色数据在采集-提取-传输-配色过程中容易产生失真的问题,因此需要运用光谱反射率重建技术实现准确校正和复制还原。
对光谱反射率重建技术的研究得到了社会各界的高度重视。Cohen测量了433个标准色卡的光谱反射率,利用主成分分析法(PCA)高度拟合出色卡的光谱反射率[11-12]。Hardeberg等[13]提出主特征向量法,解决了伪逆法在光谱重建中需要相机通道数多的问题。Tominaga[14]利用二色反射模型来描述物体反射的光,利用线性有限模型估计物体表面的光谱反射率。岑奕等[15]利用光谱重建技术对唐卡的矿物颜料进行分析,以确定唐卡主色中各颜料的用量,达到颜料修复的目的;张璐提出Sine-SSA-BP算法,主要将经过Sine混沌映射改进的麻雀搜索算法[16]和BP神经网络相结合以提高光谱重建精度[17];孔祥伟[18]为了对纱线颜色进行测量,提出了基于反馈式多光谱成像方法对纱线颜色进行测量。龙群燕等[19]提出了一种基于R矩阵与正则化多项式相结合的光谱重建算法,能够实现较高精度的壁画颜色还原。王可等[20]运用主成分分析法、R矩阵和正则化R矩阵,利用壁画色块颜色进行验证试验,实验结果表明,使用正则化R矩阵重建的光谱具备较高精度。赵海等[21]在主成分分析法基础上引入加权系数,以及误差校正函数,获得了较高的光谱反射率重建精度。何鹏浩等[22]针对现有单图像超分辨率卷积神经网络存在模型参数过多以及重建失真过大的问题,提出了一种基于动态金字塔结构与子空间注意力模块的轻量级单图像超分辨率网络模型,马子杰等[23]为了得到置信度更高的超分辨率先验模型,实现重建结果在噪声和细节之间的平衡,建立了基于混合稀疏表达框架下的高斯—洛伦兹混合先验模型。研究了该先验模型在超分算法中的应用优势和具体的应用方案。
本文中提出的加权平均的光谱重建方法(TWA),基于数码相机RGB值的背景光谱反射率重建及颜色校正开展研究,自制绿色和沙土色色板,采用多种重建方法进行光谱反射率重建实验,分析对比重建效果,为实现背景颜色从采集-传输-配色-显示全过程高精度校正还原奠定技术基础。
光谱反射率重建技术主要有2种数据集的采集分类方式,第1种是通过分光光度计获取数据集,其对应的主要为直接重建法,包括4种具体方法:伪逆法[24]、维纳估计法[25]、有限维模型法、多项式回归法[26]等。第2种是通过多光谱相机获得数据集,主要对应3种重建方法:插值重建法、 基于人工智能的重建算法[27-28]以及组合重建法。
传统最小二乘法就是将图像的数字信号进行多项式拓展获得多通道的数字信号,然后以多通道的数字信号为基础,利用指定的训练样本按照式(1)求解光谱反射率重建矩阵。
Q=Rtrain(Dtrain)-1
(1)
其约束条件为:
(2)
式(2)中:Q为光谱反射率重建矩阵;Rtrain为训练样本的光谱反射率矩阵;Dtrain表示数码相机的RGB拓展值;“-1”表示伪逆运算,“‖·‖”表示矩阵或向量的范数。获得光谱反射率重建矩阵之后再进行光谱反射率重建。
为了解决训练数据过拟合问题,正则化最小二乘法(RLS)采用Tikhonov正则化方法对光谱反射率重建矩阵的求解进行优化,以达到求解伪逆矩阵系数稳定的目的。如式(3)所示:
Dtrain=UWVT
(3)
W=S+αE
(4)
其中,α为正则化参数,是一个极小值,需要在实验过程中综合分析数据进行求解。
偏最小二乘法将经过多项式拓展的响应值矩阵Dtrain分解为得分矩阵M和载荷矩阵N,如式(5)所示:
Dtrain=MNT+I
(5)
式(5)中:I表示残余矩阵;符号“T”表示矩阵的转置。同样的方法可以将训练样本的光谱反射率矩阵Rtrain分解为得分矩阵H和载荷矩阵K,如式(6)所示:
Rtrain=HKT+F
(6)
式(6)中:F为残余矩阵。偏最小二乘法的实质就是使得得分矩阵M和H之间的协方差最大化,如式(7)所示:
H=MB
(7)
式(7)中:B为回归矩阵。利用迭代的方式求出得分矩阵M之后,即可利用式(8)计算出光谱反射率矩阵,进而计算得出光谱反射率重建矩阵。
QT=D(NTD)-1(MTM)-1MTRT
(8)
主成分分析法通过对训练样本进行主成分分析,通过比较各成分对计算结果的影响程度(通常以百分占比的形式呈现),影响程度最大的即为主成分,其余为次成分。得到训练样本集的光谱反射率矩阵的系数矩阵和特征矩阵。如式(9)所示:
Rtrain=UkM
(9)
式(9)中:Uk表示由前k个特征向量组成的特征矩阵;M为相对应的系数矩阵。对于随机的测试样本,均可由式(10)计算出其对应的特征向量。
(10)
式(10)中:α为随机测试样本在主成分分析特征向量空间中的特征向量;m表示多项式的阶数;al表示各多项式的系数;r、g、b为随机测试样本的数码相机响应值。最后通过式(11)计算出测试样本的光谱数据。
Rtest=Ukα
(11)
根据专业背景颜色特点,对传统的伪逆法进行改进,提出加权平均的光谱重建方法,引入加权函数对每一个训练样本赋予不同的权值。使得训练后更加逼近真实反射率,确定每一块训练样本所对应的权值之后,构建随机测试样本的光谱重建矩阵。
首先,根据式(12)计算测试样本与所有训练样本的色差。
(12)
式(12)中: ΔEi为测试样本与第i个训练样本之间的色差;Ltest、atest、btest为测试样本在CIE1976L*a*b*均匀颜色空间中的L*a*b*值;Li、ai、bi为第i个训练样本的L*a*b*值,n=1,2,3,…,N,N为训练样本的个数。然后,按照式(13)对色差进行归一化处理,并按照式(14)计算权值:
(13)
(14)
式(13)中: ΔEi为测试样本与训练样本的色差;ΔEimax为测试样本与训练样本色差的最大值;mi即为加权系数。构建的加权矩阵为
(15)
获得加权函数之后,对测试样本和训练样本的RGB值进行拓展,以一阶线性拓展为例,拓展之后的响应值可由dexp表示。
dexp=[R,G,B,1]T
(16)
式(16)中:R、G、B表示样本的数码相机响应值,符号“T”表示矩阵的转置。然后利用训练样本的光谱反射率矩阵,响应值拓展矩阵以及加权矩阵,由式(17)计算该测试样本的光谱反射率重建矩阵Q。
Q=RtrainM(Dtrain,expM)+
(17)
式(17)中:Rtrain表示训练样本的光谱反射率矩阵;M表示加权矩阵;Dtrain,exp表示训练样本的响应值拓展矩阵。符号“+”表示求矩阵的伪逆矩阵。最后,由式(18)完成对测试样本的光谱反射率重建。
Rtest=QDtest,exp
(18)
基于数码相机RGB值的光谱反射率重建研究是基于训练的方式进行的,首先需利用数码相机和光谱仪分别获取训练样本的RGB值和光谱反射率数据,再对所提取的RGB值和光谱反射率数据进行训练,获得光谱反射率重建矩阵,最后利用光谱重建矩阵对测试样本响应值进行光谱反射率重建。
实验仪器:NikonD90相机、OHSP-660光谱透反射率测试仪、Az-Pro100全光谱光源。
标准色板:自制沙土色系列色板(115块),规格:5 cm×15 cm;色系列色板(98块),规格:30 cm×30 cm。
实验中的成像设备选用NikonD90数码相机1台,光源选用Az-Pro100全光谱光源,实验环境在暗室中进行,以保证所选用的光源为实验中的唯一光源,排除其他光源的干扰。在暗室内垂直拍摄沙土色色板和绿色色板的数字图像,相机的拍摄参数分别设置为相机的最佳参数,相机距离色板的距离为50 cm。
1) 标定相机的光谱灵敏度,得到相机在最佳参数时的光谱灵敏度;
2) 利用NikonD90数码相机在搭建的实验环境下,获取训练样本数码相机响应值,利用光谱仪获取训练样本的光谱反射率数据;
3) 利用训练样本的数码相机响应值与光谱反射率数据,采用传统伪逆法求解训练样本的光谱反射率矩阵;
4) 采用基于加权平均的光谱反射率重建方法,对测试样本的光谱反射率进行重建;
5) 从光谱评价指标和色度评价指标2方面实验样本的光谱反射率重建精度及色差值进行分析评价。
利用115块自制沙土色和98块绿色色板作为训练样本,随机选取一块色板作为测试样板,分别运用传统最小二乘法(OLS)、正则化最小二乘法(RLS)、偏最小二乘法(PLS)、主成分分析法(PLS)和加权平均法进行光谱反射率重建,并将重建之后的光谱与用光谱仪测试出来的真实光谱反射率比较。首先以115块自制沙土色色板为训练样本进行光谱反射率重建。
3.1.1传统最小二乘法(OLS)
首先利用光谱仪获取训练样本的光谱反射率矩阵R,利用数码相机拍摄训练样本并获取色板的响应值(RGB)矩阵D,之后按照式Q=Rtrain(Dtrain)-1求解出光谱反射率重建矩阵QOLS:
(19)
之后,对于任意的一块色板即可计算出重建之后的光谱反射率,并将重建之后的光谱反射率绘制成曲线见图1。
3.1.2正则化最小二乘法(RLS)
首先利用Matlab将数码相机的响应值矩阵进行奇异值分解,分别求出数码相机响应值矩阵D的奇异值矩阵S、正交矩阵U和V,通过对奇异矩阵S加入一个极小值,来解决过拟合问题,之后求解出光谱反射率重建矩阵QRLS:
(20)
获得光谱反射率重建矩阵之后,对于任意的一块色板分别计算出重建之后的光谱反射率,并将重建之后的光谱反射率绘制成曲线见图1。
3.1.3偏最小二乘法(PLS)
首先将光谱反射率矩阵R和数码相机响应值矩阵D分解为得分矩阵和载荷矩阵,在Matlab中以得分矩阵M和H之间的协方差最大为约束条件,求出得分矩阵M,进而求解出光谱反射率重建矩阵,并将重建之后的光谱反射率绘制成曲线见图1。
图1 各种方法重建之后的沙土色光谱反射率曲线
3.1.4主成分分析法(PCA)
利用Matlab中的princomp函数对训练样本的光谱反射率矩阵进行主成分分析,选取前3个特征值,求解出系数矩阵进而获得重建的光谱反射率。并将重建之后的光谱反射率绘制成曲线见图1。
3.1.5加权平均法(TWA)
首先利用光谱仪获取训练样本的光谱反射率矩阵R,利用数码相机拍摄训练样本并获取色板的响应值(RGB)矩阵D,然后分别将训练样本与测试样本(随机选取一块色板作为测试样本)的颜色值转化为L*a*b*值,并且在CIE 1976L*a*b*颜色空间中计算测试样本与每一块色板之间的色差,共有115个色差值;获得色差之后计算出每一块训练样本对应的权值,利用所有的权值构建一个115×115的加权矩阵M:
(21)
获得加权矩阵之后计算得到一个401×3的光谱反射率重建矩阵QTWA。
(22)
最后完成对测试样本光谱反射率的重建,并将重建之后的光谱反射率绘制成曲线如图1所示,real线表示真实的光谱反射率,其他表示经过各种方法重建之后的光谱反射率,从图中可看出,采用加权平均光谱反射率重建方法(TWA)重建的光谱反射率曲线与真实的光谱反射率曲线最为接近,在接近趋势以及曲线平滑度上均优于现有方法。
绿色(及其他色彩)色板的计算方法与沙土色计算方法一致,所获得的光谱反射率重建矩阵为一个401×98的矩阵。并将重建之后的光谱反射率绘制成曲线见图2。图中real线表示真实的光谱反射率曲线,其他表示经过各种方法重建之后的光谱反射率,从图中可看出,采用加权平均光谱反射率重建方法(TWA)重建的光谱反射率曲线与真实光谱反射率重合程度最高。
图2 各种方法重建之后的绿色光谱反射率曲线
(23)
(24)
式(23)和式(24)中:R1表示由重建的光谱反射率组成的向量;R2表示由真实的光谱反射率组成的向量;N表示在可见光范围内采样点的数量;“T”表示矩阵的转置。
(25)
各方法对实验样本的光谱反射率重建精度及色差值计算结果如表1所示。表中OLS1、OLS2、OLS3和OLS4为分别采用一阶、二阶、三阶和四阶多项式对实验样本进行处理计算。
表1 各方法对实验样本的光谱反射率重建精度计算结果
续表(表1)
运用本文中提出的加权平均法(TWA)与现有的主成分分析法(PCA)等方法对沙土色和绿色色板的光谱反射率进行重建,并对重建效果进行分析比较: