HSI数据产品分0~5共6个级别, 包括原始数据产品、辐射校正产品、系统几何纠正产品、几何精纠正产品、高程校正产品和标准镶嵌图像产品。普通用户通过访问中国资源卫星应用中心的网站便可以免费获得HSI数据的2级产品, 即系统几何纠正产品。本研究即在HSI 2级产品基础上进行的。

HSI 2级产品经过辐射校正和系统几何纠正。几何定位精度在1 km左右, 误差分布比较均匀, 纠正后的影像映射到指定的地图投影坐标系下^[1], 其特性如表1所示。

表1

Tab.1

表1(Tab.1)

表1 HSI 2级产品特性说明 Tab.1 Introduction of HSI level 2 data product 参数 产品特性 波长/nm 450~950 波段数/个 115 地面分辨率/m 100 数据类型 16位有符号整型 数据格式 HDF5 辐射校正情况 经过辐射校正 几何纠正情况 经过几何粗纠正

表1 HSI 2级产品特性说明 Tab.1 Introduction of HSI level 2 data product

以吉林省洮南县为中心的Hyperion数据与环境星HSI数据共同覆盖范围为实验区, 实验区中心位置为东经121° 44'31″、北纬45° 53'54″, 其空间位置及遥感数据如图1所示。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 实验区空间位置及遥感数据Fig.1 Position of study area and data

图1遥感数据为HSI第105、70、40波段的假彩色合成图像。实验区由北向南依次覆盖了乌兰浩特市南部、科尔沁右翼前旗东南部、洮南市西部和突泉县东北部的小部分地区。

HSI 2级产品数据集为辐亮度产品, 以无符号整型数据记录, 数值范围为0~65 535。实际上, 由于地物的辐射亮度值很小, 在产品生成时乘了一个扩大系数100, 因此在应用时, 只需要将各波段数值除以系数100, 即可以得到单位为W/(m²· sr· μ m)的绝对辐射亮度值图像。

在HSI数据2级产品的部分波段上可以看到明显的条纹, 这些像元灰度值比较小, 与周围有明显差异, 条纹的存在严重影响了图像的质量和应用。波长较短的波段条纹较清楚, 与周边差异明显, 但数量较少; 波长较长的波段条纹较模糊, 与周边差异不明显, 但数量较多。图3为波长479.6 nm和波长529.6 nm的原始数据产品, 从图中可以明显地看到条纹。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 实验区HSI第10波段(479.6 nm, 左图)和第30波段(529.6 nm, 右图)原始图像Fig.3 The original image of band 10(479.6 nm, left) and band 30(529.6 nm, right)

采用“ 全局去条纹” 方法^[2]去除影像中的垂直条纹。对于HSI数据2级产品, 由于已经过系统几何纠正, 图像发生了旋转, 图像上的条纹也不再呈垂直分布, 因此在进行去条纹之前需要先将图像进行一定的旋转, 使条纹能沿垂直方向分布。

2.2.1 图像旋转

快速、自动地确定图像旋转角度, 减少人工参与, 可以在很大程度上保证处理方法的精确性和普适性, 从而大大提高图像处理效率。通过确定图像两个顶点坐标, 并求得两点所在直线与法线夹角的正切值, 进而确定该夹角为图像旋转的角度, 如图4所示。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 图像旋转角度计算Fig.4 Calculation of image rotation angle

利用图像背景值为0这一特点, 通过求每列像元值的和, 确定第一个和不等于0的列号作为x₁, 将这一列中第一个值不等于0的行号作为y₁。用类似的方法确定x₂和y₂。通过多次试验, 确定图像旋转角度α 的计算公式为

α =arctan[(x₂-x₁)/(y₂-y₁-1)](1)

将各波段数据旋转角度α 后, 图像上的条纹基本上均沿垂直方向分布, 这样即可以采用“ 全局去条纹法” 进行条纹去除。

2.2.2 垂直条纹去除

针对条纹垂直分布的特征, 采用全局去条纹法^{[2, 6]}, 通过像元的列平均值、标准差和波段平均值、标准差之间的差异对像元进行分波段线性化修正, 消除垂直条纹的影响, 即

DN'_{i, j, k}=g_{i, k}· DN_{i, j, k}+b_{i, k}(2)

式中, DN_{i, j, k}和DN'_{i, j, k}分别为原始和消除垂直条纹影响后的像元值。增益g_{i, k}和偏移量b_{i, k}计算公式为

g_{i, k}=Stdv_k/Stdv_{i, k}(3)

b_{i, k}= DNk¯-g_{i, k}· DNi, k¯(4)

式中, Stdv_k和 DNk¯分别为第k波段像元值的标准差和平均值; Stdv_{i, k}和 DNi, k¯分别为第k波段第i列像元值的标准差和平均值。

垂直条纹去除前后效果对比如图5所示。可以看出, 原有的垂直条纹已经被去除掉, 处理后图像的灰度值在空间上过渡平滑。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 第30波段垂直条纹去除前(左)后(右)效果对比Fig.5 Comparison between before(left) and after(right) destripe in band 30

将去除条纹后的图像按之前计算的角度进行反向旋转, 即完成图像的条纹去除。

2.2.3 条纹去除效果评价

通过对原始图像和去除条纹后的图像进行最小噪声分离变换(Minimum Noise Fraction, MNF)来识别和分离数据中的噪声^[7]。经过变换, 高光谱数据中的噪声基本被分离出来了, 通过对比去条纹前后MNF特征图像, 可以非常明显地看出条纹去除的

效果。图6比较了去条纹前后的MNF特征图像第10波段。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 去除条纹前(左)后(右)第10波段MNF特征图像比较Fig.6 Comparison between before(lift) and after(right) destripe in MNF(band 10)

2.3.1 光谱吸收特征波段的确定

本研究采用软件FLAASH(Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes)进行大气校正。由于软件中并没有与HSI数据相对应的传感器类型, 因此需要自定义一个描述光谱吸收特征波段的文件(Spectrograph Definition File)用以进行大气校正^[8]。

大气顶层的辐射亮度具有随电磁辐射波长增加而减弱的趋势, 这是因为, 在传感器端的太阳辐射是乘性因子。通过计算各波段大气顶层反射率, 可以去除太阳辐射的影响; 由大气顶层反射率曲线可以确定820 nm附近的水汽吸收谷波段位置。将相关的波长、FWHM以及波段号等信息输入FLAASH模型可以计算水汽透过率。由于HSI数据的波长范围最长只到950 nm, 在FLAASH中无法进行气溶胶反演, 因此只能根据影像情况设定初始能见度为40 km。大气顶层反射率计算公式为

ρ _p= πLλd2ESUNλcosθs(5)

式中, ρ _p为大气顶层反射率; L_λ 为大气顶层辐射亮度; d为日地距离; ESUN_λ 为各波段大气顶层的平均太阳辐照度; θ _s为太阳天顶角。

2.3.2 大气校正

首先, 在FLAASH软件中输入相关参数, 包括图像成像时间、地理位置、海拔、大气模式以及气溶胶模式, 定义吸收特征波段描述文件, 然后进行大气校正。典型地物(植被)大气校正前后光谱曲线对比如图7所示。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 HSI数据大气校正前后植被反射光谱曲线Fig.7 Spectral curves before and after the atmospheric correction

从图7可以看出, 经过大气校正, 近红外波段范围内由于水汽吸收而形成的反射率低谷明显消失, 植被反射光谱曲线的特征基本保持。

2.3.3 光谱曲线锯齿的去除

经过大气校正后的图像, 部分像元的光谱特征曲线存在光谱抖动, 这是由于在辐射定标等数据预处理过程中不可避免地会引入一些误差和噪音, 因此需要对大气校正后的图像进行处理。

通过MNF变换将数据分离成包含有用信息且连贯的MNF特征波段和以噪声为主的MNF特征波段, 将含有信息量大的前10个波段进行MNF反变换, 得到光谱优化后的图像。这种方法可以有效去除光谱抖动, 而且不会降低图像的空间分辨率^[6]。图8为光谱优化前后HSI数据光谱特征曲线的比较。从图中可以看到, 曲线上不正常的光谱抖动和锯齿基本上都被消除了。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 HSI数据光谱优化前后植被反射光谱曲线Fig.8 Spectral curves before and after the optimization

以地形图或经过几何精纠正的图像为参考, 通过人工在一景图像上选取25~30个控制点, 利用二次多项式纠正模型和最邻近像元重采样方法对大气校正后的HSI数据进行几何纠正, 误差控制在一个像元以内。

Hyperion作为第一个星载民用成像光谱仪, 已运行近7 a。经过多年的研究, 其数据的预处理流程已十分成熟和完善, 得到的地表反射率也比较可靠。本研究采用相同位置且同期的一景Hyperion数据进行标准化预处理^{[3, 9]}, 得到地表反射率数据, 并以此作为参照, 对比分析HSI数据预处理结果的可靠性, 如图9所示。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 HSI(左)与Hyperion(右)地表反射率数据Fig.9 HSI (lift) and Hyperion data (right) for comparison

Hyperion数据的空间分辨率为30 m, 光谱分辨率在10 nm左右, 其波段设置与HSI数据并不相同。为了使被对比的数据与HSI数据在空间和光谱上具有一致性, 本研究将Hyperion数据的预处理结果与HSI数据的地表反射率结果进行了几何配准, 误差控制在1个像元以内, 并降空间分辨率至100 m, 使之空间尺度保持一致; 将Hyperion数据中与HSI数据相同光谱范围的波段提取出来, 并插值成光谱分辨率为1 nm的数据, 结合高斯函数模拟出HSI光谱响应函数, 得到模拟的HSI数据, 作为参照。

2.5.1 地表反射率一致性分析

对450~950 nm光谱范围内的HSI与Hyperion模拟的光谱反射率进行逐像元相关分析, 其结果如图10所示。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 逐像元的地表反射率相关系数及其直方图Fig.10 Correlation coefficient of two images and its histogram

可以看出, 在整个研究区, HSI与Hyperion模拟数据在整个波段范围内表现为极显著相关。从相关系数及其空间差异上来看, 两种数据的光谱反射率变化表现得非常一致。

2.5.2 差值分析

对Hyperion模拟与真实环境星HSI的地表反射率进行差值计算, 各波段差值图的相关统计信息变化趋势如图11所示。

图11

Fig.11

	Figure Option ViewDownloadNew Window
	图11 Hyperion与HSI地表反射率差值统计Fig.11 Statistics for the differential image of reflection between Hyperion and HSI

从图11来看, Hyperion模拟与真实HSI地表反射率差值的均值基本在0附近, 表明两种数据中大部分区域的地表反射率大致相同; 标准差均远小于其值域分布范围, 说明其分布集中、波动程度小。各波段的差值大体服从正态分布。以第10波段为例(图12), 差值最值约为0.1, 最值像元主要出现在不同覆被类型的边缘、破碎地区和特殊土地覆盖区域

图12

Fig.12

	Figure Option ViewDownloadNew Window
	图12 第10波段逐像元的地表反射率差值及其直方图Fig.12 Differential image and its histogram of reflection between Hyperion and HSI

(如水体), 这主要是由于两种数据在获取时的空间尺度差异造成的, 这种差异无法通过像元分辨率的变化得到完全消除, 并在地表覆被破碎地区表现的更为突出。

从相关性及差值统计分析结果来看, 经过标准化预处理得到的环境星HSI地表反射率与Hyperion模拟的地表反射率一致性较高, 说明针对HSI的数据预处理方法是有效的。