基于辐射不变点的高光谱遥感传感器交叉定标方法

李映潭 杨辽 王杰

(1 西华师范大学 地理科学学院,四川南充 637002)(2 中国科学院新疆生态与地理研究所,乌鲁木齐 830011)

对于遥感成像测控系统,数据处理不仅要关注地物的纹理和颜色等几何特征,还需要考虑物体的反射和辐射特性[1]。辐射定标是进行大气校正的关键一环,因此获取逐波段的辐射标定系数(增益与偏置)是光学影像预处理的必要工作。光学遥感影像辐射定标是将原始的像元亮度(DN)值转化为辐射值的过程[2]。卫星发射前会进行大量定标实验和飞行实验,并在升空后进行在轨校准。这些步骤有助于建立影像DN曲线与地表辐射值曲线之间的线性关系,最终将原始的DN值精确转换为辐射值。受图像传感器的非线性响应、太阳辐照度的不稳定以及传感器属性影响[3],导致遥感传感器的测量值难以准确描述真实地物的辐射特性。因此,为了正确描述地表辐射特征,获得准确的辐射值,需要对标定系数进行精确计算。高光谱遥感影像波段众多,出现了图谱合一的特性[4-7],因此,对其辐射标定的精度要求较高。有研究指出资源一号02D卫星高光谱数据传感器中的120个大气窗口波段在辐射定标方面平均存在6.18%的误差,其中有12个波段的误差超过了7%,误差满足遥感卫星应用的精度标准[8],但资源一号02D号卫星提供的标定系数仅包括增益系数,缺乏偏置系数。在实际应用中,使用默认定标参数完成影像辐射定标后,经大气校正的影像反射光谱曲线可能出现一些偏差。

遥感传感器交叉定标的基本原理是在待定标卫星和已定标卫星在轨道运行过程中,对同一区域进行遥感成像。通过比较两者的成像结果并建立关系,实现对待定标卫星传感器的辐射定标[9-12]。对于交叉定标的研究,文献[13]最早基于水卫星(Aqua)与美国国家海洋和大气管理局16号卫星(NOAA-16)数据提出了基于同步天底过境(SNO)的交叉标定方法,该方法严格限制了观测时间与成像角度,且SNO的点多为极地区域,实用性与代表性不足;故文献[14]后续提出来拓展的SNO方法,将适用范围扩展到了中低纬度地区;文献[15]结合双向反射分布函数(BRDF)模型与6S辐射传输模型实现了较大时间差异与观测天顶角差异的卫星数据交叉标定;文献[16]通过辐亮度调整因子和反射率调整因子补偿高分一号(GF-1)卫星与参考卫星之间的时间差异、天顶角的差异以及光谱响应差异,从而提高了交叉标定精度;文献[17]以准不变定标场(PICS)交叉点,降低了标定不确定性。

受图像传感器的非线性响应、太阳辐照度不稳定等因素影响,遥感传感器测量值难以准确描述地物真实的辐射特性。因此,精确获取标定系数成为获取地表辐射特征的关键。交叉定标的实施成本低、频次高,已成为验证光学遥感载荷辐射定标系数的常规方法,但多被用于多光谱卫星数据。为解决高光谱卫星传感器定标参数的精确计算,文章通过交叉定标思路提出了一套完备的定标系数获取方法。

1 方法原理

对于常规空间遥感卫星辐射定标系数的获取与校正,需要配合地面光谱仪多点实地测量或专用定标场地,但空间卫星成像时间短,难以实现仪器点测量与遥感面测量的同步。其次,受限于卫星传感器的空间分辨率,点测量的光谱仪由于混合像元效应,与其对应的像元光谱存在偏差。因此,为实现快速标定系数的获取与校正,本文基于实施成本低且频次高的交叉定标思路,引入NASA的地表矿物粉尘源调查(EMIT)高光谱数据,通过在与资源一号02D卫星经过同一区域的轨道交叉点采集数据时,利用空间、光谱与时间分辨率都较高的哨兵二号(Sentinel-2)数据寻找辐射不变点,逐波段建立关系,并得到一套新的定标系数。

选取辐射不变点是两种数据进行交叉定标的首要工作,通过目视解译选取辐射不变点费时费力,且目视解译还存在主观性的偏差,同时两种数据分别为DN值与辐射值,无法从光谱维度获取信息,存在信息上的缺失。基于以上因素,本文使用Sentinel-2卫星的二等A级数据(L2A)数据来选取辐射不变点。

在遥感应用中,光谱角(Spectral angle)是一种比较两个光谱相似度的指标,通常用于遥感图像处理和分类[18-20]。给定两个光谱向量A和B,可以表示为包含多个波段反射率值的向量。光谱角(θ)可以通过以下公式计算

(1)

在本文中,光谱向量A与B分别代表两种影像同一坐标位置对应的光谱向量,A·B表示向量A和向量B的点积(内积),‖A‖和‖B‖则分别表示向量A和向量B的范数(或长度)。运算结果在0～1,数值越大,表明两个向量越相似。光谱角的引入使得Sentinel-2影像的辐射不变点选取更具有科学性,且光谱角从角度度量,不受成像时间和太阳辐照度不同引起的亮度和对比度变化的影响。通过编写交互式数据语言(IDL)代码对两景影像逐点求取光谱角,得到一个角度结果矩阵,取矩阵中前100个最大值的坐标作为辐射不变点。相关代码如图1所示。

图1 光谱角匹配的处理代码

将100个辐射不变点放入EMIT与资源一号02D卫星数据进行位置检查,共删除了18个在云区的点,保留82个点用于计算逐波段辐射标定系数。编写IDL代码,基于辐射不变点坐标,以资源一号02D卫星数据的DN值为自变量,EMIT数据的辐射值为因变量,逐波段进行线性回归。标定结果为一个逗号分隔值(CSV)表格,该表格包含了165个波段所对应的偏置、增益相关系数与标准误差。

对处理完后的影像,以IDL代码为基础,结合重新标定的系数,对资源一号02D数据进行逐波段的辐射定标;将得到的辐射值影像导入基于可视化图像环境(ENVI)软件的快速视线大气分析光谱超立方体算法(FLAASH)模块编写的批处理程序进行大气校正处理(相关参数设置见ENVI帮助文档),得到反射率影像;最后在删除误差较大的18个坏波段后,对结果进行3×3模板的光谱平滑。整体的过程框架如图2所示。

图2 交叉辐射标定与大气校正过程框架

2 数据源与数据预处理

资源一号02D卫星于2019年9月12日成功发射,搭载的两个传感器分别提供9波段多光谱数据和166波段高光谱数据[21]。多光谱数据的光谱覆盖范围为0.452～1.047μm,全色波段的空间分辨率为2.5m,其他8个波段的空间分辨率为10m,成像宽幅为115km。高光谱传感器的光谱覆盖范围为0.40～2.50μm,空间分辨率为30m。在可见光-近红外波段,光谱分辨率为10nm,短波红外波段为20nm,成像宽幅为60km。卫星的重访周期为3天,全球覆盖周期为55天。本文使用的是L1级产品数据,该数据级别为未地形校正的DN值影像,其中影像主文件为地理标记图像文件格式(GeoTiff)数据文件,同时附带可拓展标记语言(xml)元数据说明文件、有理多项式系数模型(RPC)校正参数文件、预览图文件、覆盖区域矢量文件、观测几何角度文件以及定标系数文件。其中,影像数据文件按波段分别存为可见近红外谱段(VN)与短波红外谱段(SW)两个GeoTiff文件。

EMIT是部署在“国际空间站”上的传感器,于2022年8月开始测量,其主要任务是监测地表矿物,特别是地球干旱尘埃源覆盖区域,在北纬52°和南纬52°之间的阳光照射区进行矿物测量[22]。EMIT Level 2A数据以网络通用(NC)数据格式存储,包括3个60m空间分辨率的产品:反射率影像(EMIT_L2A_RFL)、反射率不确定性影像(EMIT_L2A_RFLUNCERT)和反射掩膜影像(EMIT_L2A_MASK)。反射率产品包括285个波段,光谱覆盖范围为0.381～2.493μm,光谱分辨率约为7.5nm;反射率不确定性包含了逐像元、逐波段的反射率不确定估计值(见图3)。

图3 资源一号02D影像DN值与EMIT影像辐射值

本文使用的是一景2023年5月29日的EMIT数据与一景2021年6月6日的资源一号02D卫星影像。资源一号02D卫星影像的标准处理流程包括数据打开、辐射定标、大气校正以及正射校正。本文使用该数据的DN值,因此无需进行辐射定标和大气校正:①数据打开,由于ENVI不支持直接读取该卫星数据,需要安装国产卫星支持工具插件,并选择.xml文件打开数据,该软件会自动打开可见光-近红外和短波红外两个图像文件,并自动合并为166个波段的图像文件,同时打开RPC和其他元数据文件,赋予影像基础的位置信息;②正射校正,本文使用ENVI的RPC Orthorectification Using Reference Image工具,以同一覆盖区域的Sentinel-2数据作为参考影像,30m的高级空间热发射与反射辐射计(ASTER)数据作为参考数字高程模型(DEM)。处理后,得到坐标系为通用横轴墨卡托投影(UTM)43N、空间分辨率30m的影像。最后为了统一空间分辨率,将影像重采样到60m。

本文基于IDL程序编写了EMIT数据预处理的批处理程序,包括以下步骤:①NC格式数据读取,通过NC文件ID以及目标参数的ID读取影像的相关参数并存储为变量;②地理查找表(GLT)地理编码,在第一步中读取的数据包括两个波段变量,分别逐像元记录了经度和纬度,通过地理编码将影像投影为经纬度;③增加元数据,最后影像的输出格式为ENVI标准格式,通过添加头文件信息的方式添加影像的中心波长、全波半波宽以及地图信息;④删除坏波段,通过对照光谱曲线,删除误差的波段;⑤影像重投影,将EMIT的经纬度投影转换为Sentinel-2同一投影带的UTM投影,确保二者地理位置完全对应,将两种数据地理位置链接后,通过影像对照发现无需进行额外的影像配准。

经裁剪后,EMIT共285波段,资源一号02D卫星共166波段,为了使两种数据的波段相对应,以资源一号02D卫星数据的光谱为参考,将EMIT数据的光谱重采样至166波段。但两种数据波段覆盖范围存在差异,光谱重采样后的EMIT数据最后一个波段为无值,故完全对应的波段数为165。

于ESA的哥白尼数据中心下载了两景与两种数据时间各自相近的Sentinel-2数据[23],由于L2A数据为反射率产品,采用超级分辨率将所有波段空间分辨率增强至10m,然后将所有像元除以10000.0转换为反射率,再重采样至待定标数据分辨率后裁剪到研究区域,经对照二者无需再进行地理配准,可直接作为光谱角计算输入(见图4)。

图4 两个时间节点的Sentinel-2影像

3 结果与精度评价

本文使用标准误差(Se)和相关系数(R)来评价标定结果,标准误差是用于衡量测量值与真实值之间的偏差程度的指标,数值越小表示标定结果的精度越高;相关系数则用于衡量两个变量之间的线性关系,数值范围在-1～+1,越接近1表示标定结果与真实值之间的关联性越强。165个波段评估结果及选取的辐射不变点如图5所示。

图5 自动选择的样本点以及各波段的标准误差与相关系数

为了展示逐波段辐射标定的精度,绘制了9个常用波段的回归结果,如图6所示,标注了对应的相关系数、标准误差与中心波长。

注:图中R为相关系数;Se为标准误差;λ为中心波长。

最终的成果影像与样例光谱如图7所示。在影像上均匀选取了10个点作为示例点,并展示了它们对应的辐射值曲线与反射率曲线。通过分析光谱曲线可以得知经过该方法定标的影像能够获取高质量的反射率光谱曲线。在影像上均匀的选择示例点既能确保它均匀的分布在整个影像区域,又能充分代表不同地物类型和场景,这样的选择有助于验证交叉定标的普适性和可靠性。光谱曲线的分析揭示了影像的辐射值与反射率之间的关系,对于提高影像数据的准确性和可信度具有重要意义,尤其是在定量遥感研究方面。

图7 经校正后的反射率影像以及各样点辐射值曲线与反射率曲线

通过图6分析可知,有小部分波段存在较大误差且相关性较低,研究表明:前几个波段的误差主要受到大气的紫外吸收效应影响,而在1374nm附近与1845nm附近的波段主要受到水汽吸收的影响。[8]

700～1000nm的近红外区域虽然相关性较高,但标准误差较大,原因在于近红外波段对大气参数具有敏感性,这是大气中的水汽和其他气体对这个波段的辐射吸收特性。这部分数据涉及到资源一号02D卫星数据的VN与SW数据重叠部分,由图8分析可知:以VN与SW两种成像仪于1040nm的散点为例,其分布均满足线性规律,但二者DN值大小存在数量级上的差别,故考虑为两种传感器的设计方案不同导致的该误差,且该范围SW成像仪的回归结果均优于VN成像仪,故在后续处理流程中删除了重叠部分的VN波段(即图5中标准误差曲线峰值部分)。

注:图中R为相关系数;Se为标准误差;λ为中心波长。

4 结束语

由于大气状况的变化以及遥感传感器自身损耗等因素,传感器的标定系数存在不可避免的波动,这种波动直接影响辐射定标的准确性,从而对大气校正产生不良影响,导致反射率出现偏差,最终影响后续数据的处理与分析。因此,定标系数需要及时更新。传统的辐射传输模型获取标定系数需要选择适当的模型,并且需要结合大量地表观测数据与大气参数。另外,通过光谱仪实测光谱来进行回归分析获取定标系数也需要大量实测数据[24],并涉及到多个地点和时间点的全面系统的数据收集,因此需要大量的人力和时间成本。本文提出的交叉辐射标定方法在执行上简单快捷,该方法引入了NASA高精度标定的EMIT数据作为参考,并通过IDL代码实现了大部分流程的自动化,相较于传统方法,这一套流程框架能够在短时间内完成数据处理,快速生成高质量的高光谱反射率产品。迄今为止,EMIT已经成功运行了一年多,并积累了大量的数据资源,丰富的数据储备为交叉辐射标定和后续的科学研究提供了数据支撑。在今后进一步的研究中,将会考虑对不同月份的辐射标定进行详细探讨,时间维度的分析能够更好地了解太阳辐射在不同季节和时间节点上对影像辐射标定系数的影响。通过深入研究这些时空变化,能够更准确地把握辐射标定的动态特性,为数据的高质量应用提供更为科学的依据。这一流程框架不仅为当前的科研数据处理工作提供了新思路,也为未来EMIT数据的更全面利用和优化运营提供了有益的方法参考。在高光谱交叉辐射标定方面的深入研究进一步提升了资源一号02D卫星的数据质量,使其在遥感领域具有更广泛的应用。