顾及时间因子的航空高光谱影像辐射归一化校正

伊丕源， 刘原麟， 武 鼎， 赵英俊， 刘洪成

(核工业北京地质研究院， 遥感信息与图像分析技术国家级重点实验室， 北京 100029)

遥感影像采集受到太阳辐射、大气条件、地形等多种因素的影响，导致多时相影像中的同类地物光谱特征存在差异，进而会影响到定量产品的一致性[1-2]，因此须要进行辐射归一化处理。辐射归一化包括绝对辐射校正和相对辐射归一化[3]。绝对辐射校正是将传感器获得像元灰度值转换成地表反射率，通常是通过大气辐射校正实现，主要包括统计学模型和基于大气辐射传输理论的模型。此外，为控制和减少不同卫星影像之间，由于成像条件、传感器响应等因素造成的差异，许多学者对遥感影像的相对辐射归一化开展了研究，并提出了多种校正方法[4-5]，基本原理是针对多时相遥感图像的各个波段分别建立对应的波段之间的校正方程，进行辐射归一化计算处理，从而减少各种影响因素产生的噪声。目前的相对辐射归一化校正处理研究主要是针对多时相、多传感器的卫星遥感影像，然而与卫星影像不同，航空高光谱遥感(尤其是推扫式传感器)每条航带的数据采集需要一定的时间过程，且作业时间内太阳辐射、大气条件处于持续变化之中，而这些变化必定对传感器接收的辐亮度信息和反射率反演产生一定的影响。但针对航空光谱遥感的相对辐射归一化校正方法研究较少。

高精度的反射率光谱反演是基于高光谱图像开展应用研究的前提和基础[6-7]。首先基于常用的统计学模型开展了大气辐射校正，进一步分析了时间因子对反射率反演的影响。然后，提出了顾及时间因子的相对辐射归一化校正方法，从而减少数据采集过程中时间变化对同类地物反射率校正结果的影响。

1 实验数据

实验区位于新疆卡姆斯特地区，2018年10月14日采用塞斯纳飞机搭载航空成像光谱仪(compact airborne spectrographic imager, CASI)完成了高光谱影像采集。航飞高度为海拔2 600 m，采集的CASI航空高光谱数据波谱范围：380～1 050 nm(可见光与近红外波段)，48个波段，空间分辨率为1 m。影像像元值为DN(digital number)值，与辐亮度单位之间的关系为1 000 DN=1 μW/(cm2·sr·nm)。航线敷设如图1所示，每条航线长度为80 km，航带旁向重叠率约15%，航带1的采集作业时间为13：06—13：21，航带2的采集作业时间为13：28—13：43。

2 实验方法

2.1 大气辐射校正方法

高光谱影像绝对辐射校正是将传感器获得的辐亮度值转换为反射率值。目前反射率光谱反演主要有统计学模型和基于大气辐射传输理论的模型等[8]。基于统计学模型的方法，如经验线性法，需要两个以上光谱均一、有一定面积大小的目标。实际数据采集过程中，在每条航带中布设明暗地物是难以实施的，通常是在某一航带中布设明暗地物，并以其校正系数对其他航带进行大气辐射校正。基于大气辐射传输的模型，常用的软件包括6S、Modtran、FLAASH等[9-11]。统计学模型选择常用的经验线性法，并基于大气辐射传输模型对不同时间的大气下行辐射等参数进行计算与对比分析。

2.1.1 经验线性法

假定图像像元DN值与反射率R存在线性关系为

DN=aR+b

(1)

式(1)中:a、b为校正系数。通常需要在图像中两个以上光谱均一、有一定面积的定标地物，分别为暗目标和亮目标，利用线性回归[式(1)]建立起两个目标与对应的图像像元平均值间的关系，即解算出校正系数a、b[12]。

实验过程中，在航带1的北端、航线2的南端都铺设了定标黑白布作为明暗地物，并同步测量了其光谱反射率，然后分别以其为参考，计算对应的校正系数。首先，以航带1的校正系数对航带1、航带2进行校正，然后以航带2的校正系数对航带2自身进行校正；然后选择两航带重叠区域的同名地物，对其反射率计算结果进行对比。

2.1.2 大气辐射传输法

大气辐射传输模型表征了大气对遥感传感器接收到的地表辐射信号的影响。对于可见光-短波红外遥感而言，在假定地表为水平均一朗伯体的情况下，依据大气辐射传输方程，传感器接收的辐射信号L可表示为[13-14]

(2)

式(2)中:Lpath为路径辐射；S为大气半球反照率；Fd为下行总辐射；Tv为地面像元与传感器之间的大气传输透过率；ρ为地表反射率。实验中，选用Modtran软件，基于两条航带各自的获取时间，计算出航空高光谱影像采集过程中不同时间因子对应的大气辐射传输参数，辅助分析时间变化的影响。

2.2 相对辐射归一化校正方法

对于多时相卫星影像的相对辐射归一化校正，主要有基于全局影像统计信息的辐射归一化方法和基于样本集的辐射归一化方法[15]。

基于全局影像统计信息的辐射归一化方法，如直方图匹配、平均值-标准偏差归一化法等，容易造成原始光谱特征的扭曲，不适用于高光谱影像处理[16]。

基于样本集的辐射归一化方法是通过选择不变的像元作为控制集样本进行辐射归一化，如伪不变特征法(pseudo-invariant feature method, PIF)[17]、时相不变集群法(temporally invariant cluster, TIC)[18]、全景回归法(wall-to-wall regression, WWR)[19]、基于加权不变点归一化法(weighted invariant pixels，WIP)[20]等。以上方法都是基于选择的样本集构建线性回归方程，实现辐射归一化计算，主要区别在于选择样本集的方法不同。常用的伪不变特征法是通过选取多景遥感影像中的特征地物点，以参考影像的特征点光谱值为自变量，其他影像的特征点光谱值因变量，通过线性回归求解校正系数[21]。

主要采用基于样本集的归一化的思想，在航空高光谱影像的旁向重叠区域内选取同名地物像点作为样本集，具体步骤如下。

(1)选取两条航带S1、S2重叠区域中的同名地物像元作为样本集，并分别记录S1、S2航带中同名地物像元的反射率数值、成像时间，即P1={V11，V12,…,V1n}，P2={V21，V22,…,V2n}。其中，n为同名地物像元编号，V1n、V2n分别表示第n个同名地物点在航带1和航带2中的光谱值。

(2)假定单航带影像的采集时间长度为T，开始采集时间为Tstart，设置时间差阈值tthr(0

(3)

依据记录的同名像元成像时间，将样本集P1、P2分为k个样本子集{P11,P12,…,P1k}，{P21,P22,…,P2k}，每个子集的成像时间计为tk,可表示为

(4)

(3)假定选取航带S1为参考航带(可选择明暗地物所在的航带)，则以P1k中的同名地物像元光谱值为参考值，P2k中的同名地物像元光谱值为因变量，通过线性回归求解出线性系数{(a1,b1)，(a2,b2), …, (ak,bk)}，其中a为斜率、b为截距。

(4)以每个子集的成像时间tk为自变量，线性系数为因变量，通过线性回归计算出航带S2的对应的线性校正系数{c,d}。

(5)以航带S2的像元成像时间为自变量，结合线性系数{c,d}，对航带S2进行校正计算。

3 实验结果分析

首先对大气辐射校正参数进行对比，然后在航带中选择4种典型地物，对比其反射率光谱校正的结果，分析确认时间变化的影响。进一步对相对辐射归一化方法的校正结果，与常用的伪不变特征点方法校正结果进行对比。

3.1 大气辐射校正参数对比

实际数据采集过程中，在每条航带中布设明暗地物是难以实施的，通常是在个别航带中布设明暗地物，并以其校正系数对其他航带进行大气辐射校正。基于大气辐射传输模型的校正，也是以航带采集过程中的某一时间点为参考进行大气校正参数解算。然而，航空高光谱遥感(尤其是推扫式传感器)每条航带的数据采集需要一定时间，且作业时间内太阳辐射、大气条件处于持续变化之中，而这些变化必定对传感器接收的辐亮度信息和反射率反演产生一定的影响。

3.1.1 线性回归系数对比

实验中在航线1、航带2的首端都布设了黑白布靶标，航空高光谱影像中记录的黑白布辐亮度DN值如图2所示。随着时间变化，反射率强的白布对应的辐亮度DN值变化较为明显，可见光谱段的DN值差异在220～530[即0.22～0.53 μW/(cm2·sr·nm)]。

图2 黑白布辐亮度DN值

分别基于首端、末端黑白布辐亮度值和实测的黑白布反射率进行经验线性回归计算。表1列出了部分波段对应的线性回归系数。黑白布的反射率是一定的，而不同时间对应的黑白布辐亮度存在差异，因此计算的校正系数也不同。

表1 黑白布1、黑白布2对应的校正系数

3.1.2 大气辐射传输参数对比

基于Modtran软件模拟计算获取了单航带采集时间(15 min，13：06—13：21)，以及两条航带整体采集时间内(37 min，13：06—13：43)对应的大气辐射传输分量，并进行对比。

总体上，大气上行辐射、地面至传感间的大气透过率、大气半球反照率变化非常小。大气上行辐射在可见光谱段稍大，单航带时间内在1.04×10-3～8.88×10-3μW/(cm2·sr·nm)，两航带时间内在1.43×10-3～11.78×10-3μW/(cm2·sr·nm)。大气透过率的变化，只有水汽吸收波段(820、940 nm)在10-3～10-4量级，其他波段在10-5量级。大气半球反照率十分稳定，数值变化在10-5～10-6量级。

地面接收的下行总辐射存在一定变化，但不同波段的变化也不同，如图3所示。可见光谱段(400～760 nm)差值单航带时间(15 min)内在0.30～0.49 μW/(cm2·sr·nm)(即DN值为300～490)，两航带时间(37 min)内为0.50～0.85 μW/(cm2·sr·nm)(即DN值为500～850)；940～960 nm波段受水汽吸收影响，变化较小；其他谱段的单航带时间(15 min)内变化在0.15～0.30 μW/(cm2·sr·nm)(即DN值为150～300)，两航带时间(37 min)内在0.27～0.51 μW/(cm2·sr·nm)(即DN值为270～510)。这一点与图2中白布DN值的差异相对应。

图3 不同时间对应的下行总辐射

3.2 地物反射率光谱对比分析

3.2.1 绝对辐射校正结果分析

选取实验区5种典型地物的反射率反演结果进行对比，分析时间变化对大气辐射传输模型、经验线性模型校正结果的影响。表2为航带重叠区同名地物影像采集时间。如图4～图8所示，分别为黄土、沙土、石英正长岩、盐碱的反射率光谱曲线对比。

表2 航带重叠区同名地物影像采集时间

选取航带重叠区域中不同分布位置的4种典型地物，对比其反射率校正结果。

图4中所选的黄土位于航带的北端。两条航带中的采集时间差为35 min。经验线性法校正结果中，航带1中黄土的反射率计算结果与航带2的结果数值差异较大，差值在0.02～0.96。航带2中分别采用1号黑白布和2号黑白布的校正结果也存在一定差异，但整体一致性较好。

图4 黄土反射率计算结果

图5中所选的沙土地物，两条航带中的采集时间差为29 min。经验线性法校正结果中，航带1中黄土的反射率计算结果与航带2的反射率结果在可见光谱段差值较小，在0.01～0.02，近红外谱段差值较大，在0.03～0.09。航带2中分别采用1号黑白布和2号黑白布的校正结果也存在一定差异，但整体一致性较好。

图5 沙土反射率计算结果

图6中所选的石英正长岩，两条航带中的采集时间差为15 min。15 min时间内，太阳辐照度等的变化并不大。经验线性法的校正结果主要是近红外波段略有差异。此外，在382.5、396.9、411.3 nm 3个波段出现负值，且数值不一致。主要是由于高程变化引起的，1号航带中黑白布铺设地点高程为海拔980 m，2号航带中黑白布铺设地点高程为海拔759 m，石英正长岩地物点的高程为海拔904 m。航空高光谱数据采集过程中，飞行高度一定，地面高程变化造成传感器对地距离的变化，而不同传输距离对应的辐射传输分量在短波波段差异较大[22]。图4、图5中黄土、沙土地物的高程分别为968、1 012 m，与1号黑白布、测区高程中值接近，受高程变化影响较小。

图6 石英正长岩反射率计算结果

图7中测区南端的盐碱地物，两条航带中的采集时间差为9 min，且自身辐亮度(反射率)值较高，校正结果的数值一致性较好。

图7 盐碱反射率计算结果

3.2.2 相对辐射校正结果对比

实验中，设置时间差阈值为5 min，自数据采集时间开始每隔5 min时间内的同名地物作为一个样本子集，共构建3个样本子集，以航带1为参考航带，分别计算出每个样本子集的回归系数，进一步结合成像时间计算航带2的回归系数。

此外，基于选定的同名地物像元，采用伪不变特征法，同样以航带1为参考航带，计算回归系数进行校正，与本文方法进行对比，结果如图8～图11所示。

图8 黄土反射率校正前后对比

图9 沙土反射率校正前后对比

图10 石英正长岩反射率校正前后对比

图11 盐碱反射率校正前后对比

以均方根误差作为归一化校正结果与参考反射率数值的相似性度量，均方根误差的值越小，则表明校正前后的光谱数值越接近。具体统计结果如表3所示。通过对均方根误差的统计对比，所提出的方法校正结果优于常用的伪不变特征点法。

表3 均方根误差统计

4 结论

分析了航空高光谱影像数据采集过程中，太阳辐射、大气传输等随着时间变化，对地物反射率反演带来的影响。得出如下主要结论。

(1)随着时间变化，不同时刻地物接收的下行总辐射存在一定差异，在37 min时间内可达0.27～0.51 μW/(cm2·sr·nm)，不可忽略。对于经验线性校正法而言，不同时刻的明暗参考地物所计算的校正系数不同，会对大气辐射校正结果产生一定的影响。

(2)选取两条航带重叠区域4类同名地物的反射率光谱对比，基于经验线性法的反演结果，在反射率谱形方面基本一致，但是反射率数值存在一定偏差。同名地物采集的时间间隔越长(间隔30 min以上)，校正结果的数值偏差越大。

(3)提出的相对辐射归一化方法以航带间同名地物像元作为不变特征点，并加入了时间因子影响，能够有效地去除航带数据采集过程中太阳辐射、大气变化等因素的影响，提高航带间同名地物像元的光谱一致性。