加入高程因子的航空高光谱影像大气辐射校正

伊丕源， 李瀚波， 童鹏， 赵英俊， 张川， 田丰， 车永飞， 吴文欢

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

0 引言

在山区，由于地形的存在造成山区各部位地面目标接收的太阳辐射能量有很大差异，也改变了太阳—地面目标—传感器之间的几何位置关系，从而导致传感器得到的测量值与地物真实的光谱反射率或辐射亮度之间存在差异，形成地形效应。许多研究表明，地形效应会影响甚至扭曲地表的反射特性，进而会影响地表反照率的估算，并带来较大的不确定性[1-3]。地形效应是影响遥感定量化分析与应用的主要障碍之一[4]，为提高遥感影像的利用价值，必须消除地形影响，即开展地形辐射校正。地形辐射校正研究是当前遥感科学与技术发展的重要学科前沿之一[5]，国内外学者对此开展了大量的研究，但目前的地形辐射校正方法大多是针对卫星遥感影像，着重于消除地形角度的影响[6-13]。然而，地形因子主要包括高程因子和地形角度因子，对于卫星遥感而言，由于卫星轨道通常在对地距离100 km以上，所以地形高程起伏导致的辐射传输路径变化相对较小，但对于航空高光谱遥感而言，由于其对地距离近，地形起伏造成传感器与地物之间的辐射传输路径变化较大，因此除地形角度外，地形高程变化的影响亦不可忽略。

本文主要针对地形高程变化的影响，在假定每一个高程点为水平朗伯体的前提下，计算分析了地形起伏导致的辐射传输距离变化对航空高光谱遥感成像造成的影响，并实现了加入高程因子的大气辐射校正。

1 实验数据

实验数据为2012年9月12日在青海省雪鞍山地区，采用国王飞机同机搭载的CASI(compact airborne spectrographic imager)成像光谱仪和ALTM Gemini激光雷达设备获取。设备性能参数如表1所示。

依据测区实际状况与设备性能，设置航飞高度为6 900 m，采集的CASI航空高光谱数据波谱范围： 380～1 050 nm(可见光—近红外波段)，共36个波段，光谱分辨率为19 nm，预处理后插值生成的影像数据空间分辨率为0.75 m，高光谱影像像元值为辐射亮度，单位为10-3μWcm-2·sr-1·nm-1，即1 000DN=1 μWcm-2·sr-1·nm-1。航空激光雷达数据采集使用33 kHz脉冲频率，视场角为±21°，平均点密度为0.3 个/m2，预处理后生成空间分辨率为5 m的数字高程模型(digital elevation model，DEM)。图1(a)为实验选择的2012年9月12日下午15： 08采集的单航带航空高光谱影像(采集时间600 s)，图1(b)为其对应的DEM，高程范围为3 918～5 166 m，落差达1 248 m。

(b) DEM

2 实验方法

鉴于本实验重点在于分析高程因子变化对航空高光谱成像的影响，因此研究思路为： 假定每一个高程点为水平均一朗伯体，然后选择相应的高程范围，分别解算出不同高程对应的大气辐射传输参量并对比分析。

2.1 基于MODTRAN辐射传输分量计算

在假定地表为水平均一朗伯体的情况下，传感器接收的辐射亮度可表示为[14]

(1)

式中：Lpath为大气上行辐射；S为大气半球反照率；Fd为下行总辐射；T为地面像元与传感器之间的大气透过率；ρt为地表反射率。

采用MODTRAN软件，选择中纬度夏季(middle latitude summer)大气模式、乡村(rural，visible=23 km)气溶胶模式； 设置传感器垂直观测。

选择3 500 m，4 000 m，4 500 m，5 000 m，5 500 m共计5个地面高程值。以高程3 500 m为例，在MODTRAN4.0软件中输入设置的反射率和传感器高度参数(如表2所示)，分别模拟计算传感器在航飞高度(6 900 m)、地面高度(3 501 m)以及不同地表反射率数值(0，0.1，0.2)情况下接收的辐射亮度。然后将5组模拟计算的结果代入式(1)，解算出每个波段对应的Lpath，S，Fd和T。

表2 模拟输入参数

参照以上参数设置，分别解算出其他高程值对应的辐射传输分量。进一步对5个高程值对应的大气辐射传输分量对比分析，以确定地形高程变化对传感器接收信号的影响。

2.2 加入高程因子的大气辐射校正

依据式(1)，ρt可表示为[15]

。

(2)

若经2.1节结果分析，高程变化对大气辐射传输分量有较大影响，则依据式(2)开展航空高光谱影像大气辐射校正须考虑到各项参数受高程变化造成的影响。

实验思路为依据一定数量高程值对应的大气辐射传输分量，通过函数拟合与插值计算，并结合高光谱像元对应的高程值，计算出任意高程值对应的各项大气辐射传输参数，然后完成航空高光谱影像的反射率反演计算。

基于IDL语言，在ENVI+IDL编程模型下实现以上过程。具体实现步骤如下：

1)首先，将计算的各项大气辐射传输分量依据标准高程值和波段存储构建为查找表，格式如图2所示。将各项参数按照对应的波段分开，逐波段将大气参数生成矩阵，每一行对应一个高程值，每一列对应一项大气参数。具体到本实验中，包括3 500 m，4 000 m，4 500 m，5 000 m，5 500 m共5行，以及Lpath，T，S和Fd共4列。

图2 查找表结构示意图

2)读取DEM，与高光谱影像进行坐标范围和空间分辨率匹配，并获取每个高光谱像元对应的高程值，生成高程数组。具体IDL代码为

dem=envi_get_data(fid=fid_dem,dims=dims1,pos=0)

envi_enter_data,dem,r_fid=r_fid

resize_doit,fid=r_fid,dims=dims1,pos=0,/in_memory,interp=3,r_fid=r_fid1,rfact=rfact

elevation=envi_get_data(fid=r_fid1,dims=dims1,pos=0)

3)逐波段对各项大气参数进行回归计算，解算出回归系数，进一步与高程数组运算，解算出每个像元对应的大气辐射校正参数。具体IDL代码为

lut=envi_get_data(fid=fid_lut,dims=dims_lut,pos=pos_lut[k])

HeightArr=[3500,4000,4500,5000,5500]

Lup=reform(lut[0,*,*],5,1)

Trans=reform(lut[1,*,*],5,1)

Sr=reform(lut[2,*,*],5,1)

Fd=reform(lut[3,*,*],5,1)

L0=regress(heightArr,Lup,const=L1)

T0=regress(heightArr,Trans,const=t1)

S0=regress(heightArr,Sr,const=s1)

F0=regress(heightArr,Fd,const=f1)

L[mask]=L1+elevation[mask]*L0

T[mask]=T1+elevation[mask]*T0

S[mask]=S1+elevation[mask]*S0

F[mask]=F1+elevation[mask]*F0

4)最后，逐波段、逐像元读取高光谱影像辐射亮度数据，结合大气参数，并基于式(2)进行反射率计算。具体IDL代码为

fork=0,num_bands-1do begin

band=envi_get_data(fid=fid_image,dims=dims_image,pos=pos_image[k])

image[mask,k]=(image[mask,k]-L[mask])/(F[mask]*T[mask]+S[mask]*(image[mask,k]-L[mask]))

endfor

image=reform(image,num_cols,num_rows,num_bands)

3 实验结果与分析

3.1 大气辐射传输分量分析

将3 500 m，4 000 m，4 500 m，5 000 m，5 500 m共5个高程值分别对应的Lpath，T，S和Fd进行对比分析。图3为传感器接收的5个不同高程的大气上行辐射Lpath计算结果对比。

图3 不同高程对应的大气上行辐射

图3中可以看出，高程越低，与传感器之间的距离越大，则Lpath值越大。尤其在可见光波段范围内，5个不同高程(对应不同传输距离值)的Lpath值差异较大，这主要与该波段范围内大气散射强烈有关。经统计，可见光波段范围内，3 500 m与5 500 m对应的Lpath值差异为0.23～1.23 μWcm-2·sr-1·nm-1，由于CASI影像的像元DN值单位为1 000DN=1 μWcm-2·sr-1·nm-1，则Lpath差异造成的信号DN值差异可达230～1 230之间，而且CASI高光谱影像大部分像元辐射亮度值(DN值)在2 000～8 000之间，因此这种差异影响较大。在波长大于800 nm的范围内，Lpath的差异逐渐变小，数值在0.05～0.14 μWcm-2·sr-1·nm-1之间，影响相对较小。由此可见，航空传感器在获取数据过程中，每个地面像元对应的高程不同，其Lpath值均存在一定差异，且这种差异不可忽略。

图4为传感器至5个不同高程间的大气透过率T计算结果对比。

图4 不同高程至传感器的大气透过率

图4中可见，不同高程对应的T均不相同，高程越高，对传感器距离值越小，则T越高，但除第32波段(938.4 nm)和第33波段(957.5 nm)外，其他波段的T差别在10-2数量级。而CASI数据的第32波段(938.4 nm)和第33波段(957.5 nm)处的T较低，主要是由于水汽在940 nm波长具有吸收作用。

图5为不同高程对应的大气半球反照率S的结果对比。

图5 不同高程对应的大气半球反照率

图5中可见，首先，S的数值在短波波段范围偏高，但数值也低于0.25，在近红外波段范围内数值低于0.05。另外，高程越高，对应的S数值越小。

图6为不同高程接收的下行总辐射Fd的结果对比。

图6 不同高程接收的下行总辐射

图6中可见，高程越低，接收的Fd值越大。经统计，3 500 m与5 500 m对应的大部分波段下行辐射数值差异为1～2 μWcm-2·sr-1·nm-1，只有第32波段(938.4 nm)和第33波段(957.5 nm)，由于水汽影响大气透过率差异较大，从而造成不同高程对应的Fd差异更大，分别为3.1 μWcm-2·sr-1·nm-1和3.8 μWcm-2·sr-1·nm-1。因此，大部分波段下行辐射的DN值差异可达1 000～2 000，假定地物反射率为0.3，大气透过率为0.95，则此部分对传感器信号造成的DN值差异约在285～570范围，其影响亦不可忽略。

综合以上分析，不同高程对应的Lpath，T，S和Fd均存在差异，且不可忽略。因此要实现精确的山地航空高光谱影像大气辐射校正，必须考虑加入高程因子的影响。

3.2 反射率计算结果对比

本文所选实验区为高海拔无人区，因此缺乏地面实测光谱数据为参考验证。为分析加入高程因子校正后的效果，将解算的反射率结果与FLAASH大气校正处理后的结果进行对比分析。FLAASH是遥感图像反射率反演的常用方法，基于MODTRAN辐射传输模型开发[16-17]。在FLAASH大气校正过程中，同样选择中纬度夏季(middle latitude summer)大气模式、乡村(rural，visible=23 km)气溶胶模式，高程参数选择测区平均高程为4 500 m。

选择实验区中变质岩地表、苔藓植被、冰雪以及岩浆岩地表共4类地物(如图7所示)，对反演后的反射率进行对比。

(a) 所选特征地物分布

(b) 变质岩地表(1号区域，高程4 920 m) (c) 苔藓植物(2号区域，高程4 390 m)

(d) 冰雪(3号区域，高程4 950 m) (e) 岩浆岩地表(4号区域，高程4 260 m)

对比图7中所选4类地物的FLAASH校正后与加入高程因子校正后的反射率，结果如图8所示。

(a) 变质岩FLAASH校正结果 (b) 变质岩加入高程因子校正后结果

(c) 苔藓FLAASH校正结果 (d) 苔藓加入高程因子校正后结果

(e) 冰雪FLAASH校正结果 (f) 冰雪加入高程因子校正后结果

(g) 岩浆岩FLAASH校正结果 (h) 岩浆岩加入高程因子校正后结果

从图8可以看出，FLAASH大气校正结果与加入高程因子校正后的结果相比，其反射率谱形基本一致，主要差别在于短波波段，FLAASH校正计算的结果甚至出现负值，无疑是错误的。本文认为如第3.1节所述，不同高程对应的Lpath，T，S和Fd之间均存在差异，且波长越短差异越大，而FLAASH校正在假定地表为均一朗伯体的前提下，参考单一高程值进行校正，从而导致反射率计算结果存在误差，且在短波波段尤为明显。因此，对于山区航空高光谱影像而言必须加入高程因子的影响，才能实现更为精确的大气辐射校正。

4 结论

针对地形高程因子，在假定每一个高程点为水平朗伯体的前提下，计算分析了地形起伏导致的辐射传输距离变化对航空高光谱遥感成像的影响，并实现了加入高程因子的大气辐射校正。主要结论如下：

1)基于MODTRAN的计算结果表明，地形高程的变化对航空高光谱影像的成像过程存在重要影响。

2)加入高程因子的大气辐射校正结果与FLAASH大气校正结果相比，两者在反射率曲线的谱形方面总体一致，但反射率数值存在差异，尤其是在短波波段FLAASH大气校正结果中甚至出现负值，无疑是错误的。

3)无论是基于MODTRAN计算的大气辐射校正参数(大气上行辐射、大气透过率)，还是FLAASH大气校正结果，都可以发现高程因子变化的影响在短波波段更为显著，这主要与短波波段的大气散射辐射强烈有关。

高程因子的变化对山地航空高光谱影像成像过程的影响不可忽略，要实现精确的大气辐射校正就必须消除其影响。未来还需要从地形高程和地形角度2个方面共同开展研究，才能进一步实现高精度的航空影像地形辐射校正。