地物光谱仪测量中的校正方法研究*

柯衍航

（泉州华中科技大学智能制造研究院，福建泉州 362000）

0 引言

近年来随着科学技术的不断发展，高光谱技术被广泛应用在多种领域，美国ASD 公司的地物光谱仪以良好的精度与多选择的测量波长范围引领相关技术进入光谱分析时代，特别是在地质勘测、农业、植被遥感、海洋遥感、环境监测等方面的勘测应用。目前在中国，有大量的地物光谱仪在使用，遍布高校、研究所等科研机构，涉及地质、环境、海洋、气象、水利、农业、林业、生态等领域，有大量的用户群体和多数的应用案例[1-3]。

光谱仪应用中，必须对波长进行标定。国内外许多研究学者提出了波长标定的方法，例如：Youngquist 等利用白光干涉仪制造等距离的谱线；Perret 等利用法布里-珀罗干涉光片制造出等距离光线的谱线[2]。目前最常用的波长标定方法是将CCD 像元序数和波长之间的关系作线性近似，即认为光谱仪的线色散率是线性的[3]。本文主要采用的是最小二乘法的3 阶拟合，从而对波长进行标定，选取的标定点数达到了8 个，大大提高了标定精度，本文所采用的算法数组求解，还可以根据波长范围扩大，增大标定点数，从而达到提高标定精度的目的。

暗电流会产生散粒噪声影响系统的信噪比，另外，暗电流的非均匀性会影响成像质量。针对暗电流的问题，主要的解决方法有使用数值方法进行暗电流的扣除[4-5]、通过制冷抑制暗电流[6]、通过改变CCD 驱动方式来抑制暗电流[7]等。本文通过制冷技术，抑制大部分暗电流，且在控制制冷温度稳定的前提下，再通过数值的方法，通过上位机算法的编写来扣除暗电流，这样通过制冷使暗电流减小，同时控温下又保证了暗电流的稳定性，从而使得在相同积分时间下与无制冷的光谱仪相比，扣除暗电流后的光强数值大大提高。

CCD 光谱仪探测单元本身的非线性导致了CCD 光谱仪探测值的非线性。具体来说，对于不同的光照强度，测量得到的响应并不是线性的。尤其是在强度较大的情况下，响应会进入饱和区。而问题是标准光源的强度是无法调节的，无法通过直接读取光强来分析非线性。为此，不少研究学者提出了CCD 非线性矫正的方法，例如两点多段矫正法、查表法、软件矫正法、多因子加权法等。本文主要采用最小二乘法多项式拟合，进行非线性矫正[8]，解决光谱仪提前进入饱和区的问题。

光谱仪得到的响应曲线时常会出现振荡等波形变形。从理论上来说，光谱出现的振荡完全可以看作是一种“噪声”，只不过这种噪声不同于一般的电子噪声，因为光谱中的电子噪声幅值不会这么大。波形矫正主要是为了更加准确地测量色温，目前主要的解决方法有黑体轨迹Chebyshev 近似法、内插法以及三角垂足插值法等[9]。本文主要采用最小二乘法多项式拟合，计算出黑体辐射响应与理论响应的关系，进行波形矫正。

1 实验原理

CCD 的数据采集是以像元为单位，一个像元对应一个数据，对于光谱波长标定实质上就是波长对应CCD 序号像元之间的关系[10]。

光谱仪中CCD 的像元与其对应的波长之间采用如式（1）所示的3次多项式拟合。

式中：p为CCD像元的序数（本文是从第0个像元算起）；λ为像元对应的波长；c0、c1、c2及c3是4个待定系数。

假设测得n条光谱线λ（ii=0，1，2，3，…，n-1）对应的像元序数pi（i=0，1，2，3，…，n-1）可以得到下列方程：

根据方程组（2）求出4 个待定系数c0、c1、c2及c3，从而得到波长与像元的转换关系，进行波长的标定。

暗电流其实是从CCD 探测器自身引入的，由于半导体内部热运动产生的载流子填充电势阱，导致在无光的情况下，CCD 输出端仍有电流形成，称之为暗电流噪声[11]。暗电流的噪声强度与温度、电荷存储时间的关系如式（3）[12]所示。

式中：σdark为暗电流噪声强度；Id为暗电流；tint为积分时间。

根据式（3），势阱中存储时间tint的长短和暗电流噪声有关，tint越长，暗电流噪声越大。由于本文中使用的CCD 对温度不敏感，所以通过测量计算推导得到暗电流的噪声强度与积分时间的线性关系，求出每一个像元所对应的斜率ki（i=0，1，2，3，…，2 047）和截距bi（i=0，1，2，3，…，2 047），即可对暗电流进行扣除。

光谱仪的非线性是制作CCD 器件的单晶材料不可避免存在的缺陷。一般在理想状况下，CCD 输出的响应随曝光时间线性增加。由于响应越大越容易产生溢出的饱和失真，所以取饱和时的响应进行非线性矫正。通过改变积分时间，固定光源的输出功率，来观察CCD 探测器的非线性，有利于减小光源本身不稳定性的影响[13]。利用最小二乘法多阶拟合进行处理，得到拟合系数。假设计算得到的理论响应为Z=（z1,z2,z3,…,zn），测量得到的实际响应为F=（f1,f2,f3,…,fn），通过计算得到理想响应与实际响应的比值C=（c1,c2,c3,…,cn），由此得到理想状况下的线性曲线。

黑体是指一种能吸收投射到其表面的所有波长光源（即a=1）的物体，也就是说投射到黑体表面的所有光，全部被吸收[14]。绝对黑体只是一种理想化的产物，在现实生活中是不存在的[15]。因为对于绝大多数的光源，色品坐标并不能刚好落在黑体辐射线上。所以在测量色品坐标前，先利用最小二乘法多阶拟合进行处理，得到拟合系数，再利用归一化，通过理想响应与黑体辐射的比值关系，运用黑体辐射线作为参考的谱线，进行波形矫正[16-20]。

2 实验结果与分析

2.1 光谱仪波长定标

光纤线将汞-氩灯与光衰减器、光谱仪分别连接在一起，测量出所要的像元序数所对应的响应，再运用算法进行计算推导，实验框图如图1所示。

图1 实验框图

光谱仪的相关参数：2 048×14 像元背照式CCD 探测器；16 位的A/D 转换器；可测量光谱范围为300～1000 nm；CCD积分时间为1～2 ms。

本文采用多项式拟合的方法对光谱仪的波长进行标定。利用汞-氩灯的特征谱线对光谱仪进行波长定标，汞-氩灯的特征谱线有：313.16、365.01、435.84、546.08、696.54、772.40、800.62、912.30 nm。用曲线y=ax3+bx2+cx+d来拟合波长和像元的关系，其中x为像元，y为波长，得到拟合系数a=-3.22×10-9，b=1.78×10-5，c=0.39，d=284.90，通过得到的拟合系数，计算得到标定的波长。由表1 所示的实际波长与计算所得波长的比较结果，可以观察到像元转换为波长的对应值，以及实际波长与计算所得波长的最大偏差略大于0.3。

表1 实际波长与计算所得波长比较

图2 所示为波长定标后的谱图。根据上述算法，得到选定特征谱线：313.16、365.01、435.84、546.08、696.54、772.40、800.62、912.30 nm，与图2 中峰值所对应的波长基本吻合，达到了最终的标定要求。本方法还可适用于波长范围增大后，通过增加标定点个数来提高波长标定精度。

图2 波长定标后的谱图

2.2 光谱仪暗电流扣除

暗电流测量时，将光谱仪的进光口用金属帽密封，且将光谱仪内的制冷片设置控温在20 ℃，从而降低暗电流的数值，使得相同积分时间下与无制冷的光谱仪相比，扣除暗电流后的光强数值得到很大的提高，同时也保证了暗电流的稳定性，大大减小了暗电流数值的波动[20-22]。

光谱仪的相关参数：2 048×14 像元背照式CCD 探测器；16位的A/D转换器；测量的光谱范围为300～1 000 nm；CCD积分时间为1～1 000 ms，间隔为1 ms。

由图3可以观察到暗电流会随着积分时间的变化而变化。图3 中每条曲线代表不同像元（像元：0～2 047，选取其中几组数据）。由图4可以看到不同像元暗电流随积分时间变化的规律还是不尽相同的，但是可以观察到变化的规律大致是线性，每一个暗电流随积分时间的变化可表示为一个线性函数y=kx+b，其中：x为积分时间，y为响应。

图3 暗电流矫正前

图4 暗电流矫正后

经过定标，根据响应、积分时间计算出斜率k和截距b的值，得到如图4所示的线性关系。这个函数的斜率和截距均可存储到仪器中。定标结束后，根据得到的k，b值线性方程，已知积分时间，求解暗电流强度，由此扣除暗电流。

通过上述的线性关系，将计算得到的函数斜率k和截距b存储到仪器中，测得暗电流矫正前后的对比如图5 所示。由图可知，矫正后的谱图就是在矫正前的基础上扣除了暗电流的基底得到的。

图5 暗电流矫正前后对比

2.3 光谱仪非线性矫正

本文的光谱仪波长范围为可见光范围（300～1 100 nm），故可以选用卤钨灯。非线性矫正的实验装置如图6 所示。实验主要是将卤钨灯放入积分球中，光纤线将积分球与光谱仪（以300～1 100 nm 为例）连接在一起，测量出所要的积分时间对应的响应，再运用算法进行计算推导。

图6 非线性矫正的实验装置

光谱仪的相关参数：2 048×14像元背照式CCD 探测器；16位的A/D转换器；测量的光谱范围为300～1 000 nm；CCD积分时间为1～51 ms，间隔为1 ms。

针对光谱仪的非线性问题，先测量得到矫正前的数据（共有2 048组数据），观察到2 048组数据大致是线性变化的，可表示为一个线性函数。通过拟合得到线性方程y=kx+b，其中x表示积分时间，y表示理论响应，计算得到2 048 组k、b，进而可通过得到的线性方程计算得到2 048组理论响应。

取趋于饱和时积分时间的实际响应（51 ms），计算得到实际光强与理论光强的比值；再根据7 阶多项式拟合计算得到y=ax7+bx6+cx5+dx4+ex3+fx2+gx+h非线性方程，得到拟合系数a=9.78×10-32，b=-2.03×10-26，c=1.68×10-21，d=-7.06×10-17，e=1.60×10-12，f=-1.86×10-8，g=9.37×10-5，h=0.93。图7 所示为非线性拟合曲线，图中纵坐标为实际响应与理想响应的比值，当测量得到一个响应（实际响应）时，就可以知道相应的实际响应与理论响应的比值。再根据实际值与比值的关系，就能得到矫正后的理论值，实现光谱仪的非线性矫正。

图7 非线性拟合曲线

根据上述算法，进行光谱采集，表2 所示为非线性矫正后数据，取积分时间分别为5、1、15、20、25 ms中随机5组数据进行验证。从表3数据间的比值可以发现5 组数据均实现积分时间为10 ms 与5 ms 的比值均在2 左右，积分时间为15 ms与5 ms的比值均在3左右，积分时间为20 ms 与5 ms 的比值均在4 左右，积分时间为25 ms 与5 ms 的比值均在5左右，验证了算法是可行的。非线性矫正后的曲线（5～25 ms）如 图8 所 示，通过非线性矫正，随着积分时间的增加，响应也随之线性增加，两者呈线性关系，实现了光谱的非线性矫正。

表2 非线性矫正后数据（5～25 ms）

表3 非线性矫正后积分时间之间数据的比值

图8 非线性矫正后（5～25 ms）

2.4 光谱仪波形矫正

波形矫正的实验装置与非线性矫正一样。光谱仪的相关参数：2 048×14 像元背照式CCD 探测器；16 位的A/D 转换器；测量光谱范围：300～1 000 nm；CCD 积分时间为1～32 ms，间隔为1 ms。

波形矫正前如图9 所示，为改善响应曲线时常会出现振荡等波形变形问题，采用以下解决方法。首先，波形矫正前面的处理步骤与非线性矫正一致，先得到线性方程y=kx+b，其中x表示积分时间，y表示理论响应，得到所需的理论响应；利用归一化原理，计算得到理论响应与黑体频谱亮度（选取所需的波段为300～1 000 nm）的比值。再根据计算得到y=ax5+bx4+cx3+dx2+ex+f拟合非线性方程。在波形矫正中，取接近中间值的理论响应，根据中间值的理论响应（40 000左右），得到所需的比值，根据黑体辐射响应与比值的关系，就可以得到矫正后的值（理论值），实现光谱仪的波形矫正。

图9 波形矫正前

波形矫正是为了后续色温测量做准备的，所以利用黑体频谱亮度的谱图作为参考谱图。波形矫正后与参考曲线的对比如图10所示，得到的相对光谱分布与参考光谱的趋势基本一致，这是利用了归一化，实现了光谱的波形矫正。

图10 波形矫正后与参考曲线的对比

3 结束语

本文利用最小二乘法多项式拟合定标、归一化、查表法、多因子加权法、计算黑体辐射响应与理论响应的关系，通过对公式的反复推导以及多次测量、反复计算提高定标精度，集多种校正方法于一体，在多项式拟合方法的基础上，通过提高定标个数来实现大范围波长段定标的精度，使波长测量精确到0.3 nm。通过制冷技术与扣除暗电流数值相结合，使暗电流数值降低且稳定，波动小，在相同积分时间下与无制冷的光谱仪做对比，大大增强了扣除暗电流后的光强度。通过最小二乘法多项式拟合，使非线性拟合精确到1.029，解决了光谱仪提前进入饱和区的问题，大大提高了光谱强度范围。通过最小二乘法多项式拟合，使波形矫正达到与理想变化趋势基本一致的结果。这些算法在地物光谱仪的实际应用中，取得了较好的结果。地物光谱仪在地球资源研究中具有重要意义，它可以用来检测地表物质的组成，以及气象状况，可以用来进行地质调查和岩石鉴定，也可以用来识别地下矿产资源，在地球资源研究中发挥着不可替代的作用。