采用Landsat8产品算法流程的高分一号数据大气校正

张晓月，李琳琳，王 莹，张 琪，李国春

采用Landsat8产品算法流程的高分一号数据大气校正

张晓月1,2,3，李琳琳2,3，王 莹2,3，张 琪2,3，李国春4

（1.中国气象局沈阳大气环境研究所，沈阳 110166； 2. 辽宁省生态气象和卫星遥感中心，沈阳 110166； 3.高分辨率对地观测系统辽宁数据与应用中心，沈阳 110166； 4.沈阳农业大学农学院，沈阳 110866）

高分一号（GF-1）卫星搭载的传感器，实现了高分辨率和宽幅成像能力的结合，使其在精准农业方面应用发挥重要作用。该文在6S大气辐射模拟模型基础上，参照LaSRC大气校正流程，设计了GF-1卫星WFV/MSS数据从算法原理分析到编码实现的大气校正。算法应用大气总传输率、水汽透过率和大气后向半球反照率等参数和高程、大气可降水以及臭氧含量等值，循环计算使GF-1像元红蓝通道值比等于来自MODIS数据的红蓝通道值比，求得GF-1像元气溶胶，再将包含当日大气可降水和臭氧含量的辅助数据文件等代入6S模型，得到大气校正后的地表反射率。试验表明，该大气校正方法在有农田和林木等植被覆盖的中低纬度大气校正效果较好，对稀疏植被的荒漠裸地和建筑地表大气校正效果相对稍差。比较GF-1 WFV/MSS数据与Landsat8(LC8) OLI数据的基于6S模型LaSRC流程算法的大气校正结果，GF-1 WFV/MSS各传感器与LC8 OLI大气校正结果的相关系数为0.825～0.972，2种卫星数据大气校正结果相关性高，其中WFV相较于MSS显示出与LC8 OLI更相近的大气校正结果。结果表明，应用6S模型原理参照LaSRC校正流程设计的自行估计数据逐像元水平气溶胶参数的GF-1卫星数据大气校正方法应用方便、可操作性强，适合生长季农林监测等陆面应用。

卫星；遥感；GF1 WFV/MSS；LC8 OLI；6S；LaSRC；大气校正

0 引 言

高分一号卫星具备时间分辨率和空间分辨率双高的特性，因其轨道稳定性和数据可持续性等诸多优势在多领域取得了较好应用效果[1]。高分一号数据不仅局限在对地表地貌的观测，而是越来越多用于定量应用遥感，农业方面主要在农作物种植区划、农业用地精准提取、农作物分类识别、农田面积估算、农作物产量预报等领域发挥着重要作用[2]。但众所周知卫星传感器接收到的地面辐射受大气分子、气溶胶以及臭氧、水汽等气体的散射、吸收影响会造成信号失真，因此在应用高分一号卫星遥感数据进行地表参数定量反演时，应尽可能将大气噪声从地面辐射信息中剥离，保留地面物体真实反射率信息。大气校正是解决这一问题的重要手段，按照是否具备明确物理意义分为相对校正模型和绝对校正模型[3]。目前较为常用的大气校正方法有图像特征法、经验线性法和大气辐射传输理论法[4]。图像特征法是一种相对校正方法，针对图像特征消除大气影响，常见的有直方图均衡法[5]、暗目标法[6-7]、不变目标法[8]、内部平均法（IARR，internal average relative reflectance）[9]、平场域法（FF，flat field）[10]等。经验线性法结合野外试验，依赖于地面同步定标点实测数据[11]。基于辐射传输理论（RT，radiative transfer）的大气校正方法具备明确的物理意义，该类方法以辐射传输理论为基础，考虑大气参数等影响，校正精度较高、应用范围广泛。

目前大气校正处理应用较为广泛的6S模型是美国马里兰大学地理系在5S[12]模型基础上研发的[13]，该模型在气体透过率、瑞利散射及气溶胶光学厚度计算等方面进行了改进，在MODIS产品大气校正计算方面表现尤为突出[14]。2005年6SV1版本公布，该版本较之前的模型主要增加了辐射极化的计算[15]，2015年更新的6SV2.1是最新版本。学者们应用6S模型进行了多种遥感数据的大气校正研究，推进了模型的改进和应用。刘佳等[16]研究了基于6S模型的高分一号卫星16 m分辨率数据大气校正，并将批量处理的效率提高了75.0%以上。孙林等[17]提出了地表反射率支持的GF-1 PMS气溶胶光学厚度反演及大气校正方法，通过对北京、太湖2个AERONET站点观测的气溶胶光学厚度验证，得到反演地表反射率与实测地表反射率误差低于0.015的结果，算法具备较高精度。孙元亨等[18]应用6S模型对GF-1 WFV数据进行了大气校正，开展了GF-4 PMS与GF-1 WFV地表反射率及NDVI一致性的分析，校正前后相关系数分别为0.74和0.77，二者表现出较好的一致性。薛现光等[19]提出在6S模型基础上，以气溶胶厚度为自变量对大气参数进行拟合的算法，对GF-1进行大气校正，与FLAASH对比具有较好的一致性。胡勇等[20]利用MODIS大气产品和6S模型对GF-1 WFV数据进行了大气校正处理，不同地表覆盖类型GF-1 WFV和MODIS反射率数据的交叉对比结果显示出较好的一致性。

基于辐射传输理论基础，近年来相关机构相继研发了一系列大气校正流程化平台和系统，使大气校正处理流程更便捷，科研人员利用这些算法软件开展了一系列大气校正相关研究。LOWTRAN模型1971年由美国空军剑桥研究所研制，已发展至LOWTRAN7版本[21]。MODTRAN模型从LOWTRAN 模型发展而来，由美国空军地球物理实验室研制，将光谱分辨率由LOWTRAN的20提高到2 cm-1[22-23]。杨贵军等[24]利用MODTRAN4和MODIS产品，对环境与减灾小卫星HSI高光谱数据进行了大气校正。郑盛等[25]利用MODTRAN4模型对HJ-1卫星CCD数据进行逐像元大气校正，并从光谱曲线、MODIS地表反射率产品、NDVI的3个方面比较了大气校正效果。熊世为等[26]利用MODTRAN模型对HJ/CCD进行了大气校正，校正结果与高精度MODIS地表反射率产品在植被、居民地、水体3类地物一致性较高。美国空军实验室1998年开始以MODTRAN为基础发起了FLAASH的研发[27]，专业遥感数据处理软件ENVI大气校正模块集成了FLAASH模型[28]，李卫国等[29]利用FLAASH进行了基于薄云雾去除的ETM+数据大气校正。ATCOR[30]模型是以MODTRAN4为基础开发的大气校正模型，ERDAS软件集成了ATCOR模型。常用的还有商用大气校正软件包ACORN[31]主要针对0.4～2.5m高光谱和多光谱数据进行大气校正。

本文分析了基于6S模型和LaSRC（Landsat Surface Reflectance Code）流程进行GF-1数据大气校正的适用性，进行了通道相关等改进处理和重新规划，设计实现了针对GF-1数据的大气校正算法和软件平台。LaSRC方法依赖于外部辅助数据，需以待校正数据同日的臭氧和水汽等大气影响气体产品作为初值，考虑到对时效性的要求和实时资料有时难以获得，平台提供了利用近6 a臭氧和水汽逐日值作为替代和使用实时值的2种方案。从而保证在无外源数据支持下仍然能够提供大气校正产品。鉴于中国境内暂无LaSRC产品，且LaSRC移植到RSD（Remote Sensing Desktop）平台处理LC8 L1T数据能够代表原LaSCR产品，本文以同地同时次LC8 OLI L1T的RSD大气校正产品和GF1 WFV/MSS数据大气校正处理结果对比分析二者大气校正效果和精度差异，探索对国产高分卫星数据的高精度大气校正处理平台软件工程化实现的实用方法。

1 原理方法

6SV模型是基于6S大气辐射传输模型基础上发展的极化矢量版本，美国USGS EROS项目针对Landsat8（LC8）数据开发了LaSRC大气校正程序。该程序运行在Linux操作系统，用于美国境内LC8数据地表反射率产品反演，对于光谱波段较窄的Operational Land Imager(OLI)传感器具有较好的适用性[32]。应用待校正LC8 OLI数据反演像元水平气溶胶，以及水汽、臭氧含量从MODIS产品或NCEP数据获取是LaSRC流程算法求算地表反射率的2个主要特点。以该流程对LC8 OLI数据进行大气校正时，首先利用查找表和大气总传输率、气体透过率、大气内反射率等6S模型参数从预定义的22个气溶胶光学厚度查找表估算地表反射率初值；从MODIS产品提取红蓝通道地表反射率比值，根据该比值与NDVImir线性关系，由LC8 OLI地表反射率初值和线性斜率及截距估算LC8 OLI红蓝通道比，在3种埃氏系数等级和22个气溶胶光学厚度下循环计算，找到使MODIS红蓝通道比和LC8 OLI红蓝通道比残差最小时的反射率，求得气溶胶；再利用MODIS产品通过空间插值得到LC8 OLI像元水平水汽、臭氧含量，利用DEM数据及气压与海拔高度关系计算并空间插值得到LC8 OLI像元水平大气压强，通过卫星参数文件得到像元几何参数；将表观反射率，气溶胶以及相关参数输入6S模型最终求得地表反射率（图1）。

图1 算法原理流程图

1.1 6S模型算法

假设地表为均匀朗伯面，大气中分子与气溶胶散射及气体吸收水平一致，卫星传感器接收到的辐射由大气路径辐射、透过大气的目标地物直接反射辐射、经由大气散射的目标地物反射辐射、背景地物散射辐射以及背景地物多次地表散射辐射构成，卫星传感器接收到的反射率表示为[33]

式中TOA为表观反射率；为大气内反射率，为地表反射率，为大气后向半球反照率，()为下行辐射传输率，()为上行辐射传输率，为太阳天顶角，为传感器天顶角。

大气内反射率又称大气路径辐射，大气中同时存在大气分子及气溶胶散射和气体吸收，简单将大气分子散射和气体吸收的作用隔离开来分别考虑显然与实际情况不符，则假设存在一种平均状况，为便于计算认为有一半的水汽混合在气溶胶中参与散射，综合考虑大气内各分子及气体耦合作用，将大气内反射率计算如下[34]

式中为瑞利散射反射率，Tg/2为半程水汽透过率。将式（1）带入表达式并改进如下：

式中为大气透过率，Tg为臭氧和其他气体的透过率。上述式中，反射率、传输率、透射率无量纲，角度单位rad。

将式（4）带入式（2），即可进行地表反射率的求算，地气解耦是实现大气校正或气溶胶光学厚度反演的重要一环。应用该方法对LC8 OLI数据进行大气校正时，气溶胶的模拟由预先建立的查找表及MODIS产品估算。

1.2 光谱响应分析

从美国USGS网站及中国资源卫星应用中心网站下载LC8 OLI、GF1 WFV/MSS光谱响应数据，制作LC8 OLI、GF1 WFV/MSS可见光及近红外波段光谱响应曲线（图2），用于计算初始的分子光学厚度参数、臭氧传输参数、水汽传输的两组参数和其他气体的3组参数。

图2 LC8 OLI、GF-1 WFV/MSS光谱响应曲线

比较二者光谱响应差异，GF-1的4个WFV传感器及2个MSS传感器光谱响应曲线相近，LC8 OLI红光及近红外波段波谱范围相对较窄，蓝绿波段光谱响应函数与GF-1 WFV/ MSS差异较小，虽然Landsat8 OLI与 GF-1 WFV/MSS对应通道并非完全相同，但波段设置较为相似的设计，使应用LC8 OLI适用的大气校正算法进行GF1 WFV/MSS数据大气校正成为较好的选择。

应用积分算法对各传感器不同波段光谱响应数据进行运算（式5），通过建立以LC8 OLI为基准的光谱匹配系数（式6），与USGS典型地物光谱库数据进行对比，修订GF-1 WFV/MSS光谱响应函数差异[35]。

式中1和2为波段范围最小值和最大值，()为波长处光谱响应值，()为波长处反射率。计算Landsat8,i和GF1,i，并进行相应波段的比较求得光谱匹配系数。

1.3 像元几何参数

本文进行大气校正时，需要求算每个像元卫星方位角和天顶角、太阳方位角和天顶角。利用随卫星元数据发放的参数文件中给定的仰俯、侧滚、偏航等角度及地理定位等信息计算，对于一条扫描线任意列数据卫星方位角和天顶角与给定参数关系如下：

按照卫星过境时间及像元经纬度进行每一像元太阳方位角和天顶角的计算，需对影像进行精准的几何（正射）校正，理论上认为对影像的几何定位与大气校正无既定处理顺序[36]。本文首先对GF-1进行正射校正，根据校正后的幅宽与日地距离比得到整幅的太阳天顶角变化，然后根据幅宽像元数线性插值，得到逐像元太阳天顶角。太阳方位角对计算反射率影响较小且变化很小，这里使用了景中心数值。

1.4 表观反射率

获取像元表观反射率首先需通过辐射定标系数将像元通道计数值转换为像元表观辐亮度，再由大气层顶太阳辐照度、日地天文距离、太阳天顶角数据将大气层顶表观辐亮度转换为大气层顶表观反射率，计算公式如下：

式中为表观辐亮度，W/(m2·sr·m)；DN为像元通道计数值；scale和offset为辐射定标系数；为日地天文距离；为太阳辐照度，W/(m2·sr·m)；为太阳天顶角，rad。

1.5 大气参数

1.5.1 气溶胶光学厚度

地表反射率与气溶胶的耦合关系使得在进行大气校正时，气溶胶的剥离存在一定难度，在一些大气校正平台中往往是根据观测点地理位置及时间选择一个气溶胶模式，如限定纬度、季节和区域类型等。这种预定义的气溶胶模式较粗糙，即使计算当日气溶胶也多是整幅的气溶胶而不是逐像元的气溶胶值，代表性较差。

本文以LaSRC流程对GF-1数据进行大气校正时，气溶胶的估算采用引入MODIS产品为辅助数据的方法，首先计算MODIS产品红蓝波段地表反射率比值，根据红蓝波段比与NDVImir线性关系，以及NDVImir与NDVI的线性关系[36]，估算GF-1 WFV/MSS红蓝波段地表反射率比值，分别比较埃氏系数为1.0、1.75、2.5这3种情况和22个气溶胶光学厚度时MODIS和GF-1 WFV/MSS数据红蓝波段地表反射率比值的残差，由残差最小时地表反射率求得气溶胶。

以2016年8月26日辽宁境内条带号121、行编号31的LC8数据为试验，分别采集有植被和稀疏植被覆盖地区样本各300个，比较NDVI与NDVImir关系。

图3 NDVI与NDVImir比较分析

图3也可看出NDVI与NDVImir线性关系明显，图中有植被覆盖地区NDVI与NDVImir数值均较稀疏植被覆盖地区大，有植被覆盖地区二者相关系数为0.958，线性回归方程常数项和一次项分别为−0.262和1.097；稀疏植被覆盖地区二者相关系数为0.912，线性回归方程常数项和一次项分别为−0.284和1.128；2个方程均通过0.05水平显著性检验。

由于计算输出结果为波长550 nm处AOT，需根据AOT与波长之间的关系[37]求算其他各波段AOT数值，公式如下：

1.5.2 瑞利光学厚度

在获取瑞利散射反射率时需以瑞利光学厚度（ROD）为输入项，ROD是波长及海拔高度的函数[38-39]，波长一定时海平面处ROD计算如下

式中()为波长处ROD。为海拔高度，m；对于任意海拔高度处ROD计算为

1.5.3 大气影响气体含量

应用MODIS产品获取逐日大气水汽、臭氧含量，该数据空间分辨率5.5 km，将其与待校正数据进行像元尺度的空间匹配。当MODIS产品实时性受限，在所需处理时段数据无逐日MODIS产品数据时，对于臭氧、大气可降水气候值的计算采用2013－2018年近6 a MODIS产品逐日动态滚动平均值进行补充，虽然不是待校正当日数据，但相较于根据经纬度及时间的查找表参数确定方式仍具有较高精度。当平台未检测到与待校正数据同日MODIS产品辅助数据时，提示并自动调用动态平均值作为输入量完成大气校正处理。

1.5.4 大气压强

大气压强数据用于计算大气后向半球反照率、气体透过率、大气总传输率、大气内反射率等参数，影像各像元处大气压强通过海平面气压及海拔高度关系计算

式中为大气压强，hPa；0为海平面气压，hPa。

2 数据及处理

2.1 试验数据

LC8 OLI数据通过美国USGS网站下载，数据级别为L1TP。GF-1 WFV/MSS数据来源于中国资源卫星应用中心网站，数据级别L1A。筛选2016—2017年GF-1数据，在比较GF-1与LC8大气校正结果差异时，选取二者过境时传感器采集数据有较大交叉覆盖范围的同日数据，LC8与GF-1均在地方时午前过境，选取同日数据大气状况差异小，便于对二者的大气校正结果进行比较分析（见表1）。

表1 卫星数据集

2.2 数据处理

对于LC8 OLI数据大气校正流程的实现，首先使用数据集提供的质量评价数据对遥感数据进行云区、云阴影、水体的检验，对识别出的像元缺漏和云等无法处理的数据进行标记隔离。通过属性文件获取定标信息进行辐射定标、获取太阳角度进行太阳高度订正，计算出表观反射率。利用查找表估算初始地表反射率，通过NDVImir估算红蓝波段地表反射率比值，与MODIS产品红蓝波段地表反射率比值对比求算AOT，最后结合表观反射率和大气气体含量以及其他参数完成大气校正。对于GF-1 WFV/MSS数据大气校正流程的实现，与LC8 OLI不同的是，GF-1 WFV/MSS为未经正射校正的L1A级数据，使用数据集提供的RPC参数结合全球30 m分辨率DEM数据进行正射校正，使二者进行精确几何配准。其次GF-1 WFV/MSS的定标系数由中国资源卫星应用中心每年更新一次，数据期间可能有变化，应用对应的LC8 OLI数据对定标系数进行交叉验证，提高定标的准确性。同时GF-1 WFV/MSS需要对侧滚、仰俯和偏航角处理得到相关角度数据（式7、8）。经上述处理后，计算得到GF-1 WFV/MSS的表观反射率，再通过NDVI求算红蓝波段地表反射率比值，完成后续大气校正处理。最后将LC8 OLI和GF-1 WFV/MSS数据进行相同区域的裁剪和配准（图4）。

图4 数据处理主要流程

3 结果与分析

3.1 LC8 OLI大气校正

LaSRC的开发是为进行LC8 L1G产品大气校正的，目前USGS共享的数据为L1T级别，由于不同级别数据的像元重定位，使校正结果存在少量偏差。利用同一时次LC8 OLI数据以及其地表反射率产品对算法实现精度及适用性进行分析。从USGS网站下载同时次LC8 OLI数据（数据集名称为：LC08_L1TP_030031_20170831_20180125_01_T1）和地表反射率产品（数据集名称为：LC08_CU_015008_20170831_20180126_C01_V01_SR），对LC8 OIL数据进行大气校正处理并与USGS LC8 OLI地表反射率产品逐波段比较（图5）。

从LC8 OLI数据大气校正结果与USGS LC8 OLI地表反射率产品中以100×100个像元的窗口采集样本点各2 500个，比较二者相关性。从图4中可以看出应用LC8 OIL数据的大气校正结果与USGS LC8 OLI地表反射率产品可见光至短波红外的第1～第7个波段相关关系均较好，决定系数为0.993～0.998，对LC8 OIL数据进行大气校正算法与美国USGS发布的LC8 OLI地表反射率产品拟合度高。本地化LaSRC方法对于LC8 L1TP级别数据的大气校正能够代表USGS LC8 OLI地表反射率产品水平。

3.2 GF1 WFV/MSS大气校正分析

对2016－2017年GF1 WFV/MSS试验数据进行大气校正处理。图6为校正前后2017-09-07日辽宁境内GF-1 WFV/MSS数据1、2、3波段合成真彩色图像，分析大气校正处理对影像的目视清晰度的影响。

注：λ为波长，n为样本容量，R2为决定系数。

图6 GF-1大气校正前后影像

从图6中可以看出，由于大气校正过程对大气分子、气溶胶、水汽、臭氧等噪声的处理，大气校正后影像清晰度提升，尤其是影像下部受大气影响明显区域，大气校正后目视判读性得到改善，这一效果是相对大气校正无法实现的。

应用下垫面地类较为丰富的2017年6月5日GF-1 MSS数据为例，该数据下垫面类型包括：裸土、农田、林地、建筑、淡水、海水，且该时段农林地表植被生长利于进行大气校正前后的地表分析，数据景中心大气层臭氧质量浓度0.642 mg/L，水汽质量浓度0.96 mg/L，气溶胶光学厚度0.35。分析大气校正对通道数据的影响，从该数据提取大气校正前后所有像元，计算逐通道每景影像所有像元混合的表观反射率均值与地表反射率均值的差值，定义为校正量（以符号表示）。

从表2看出，对于GF-1 MSS传感器的蓝、绿、红通道，表观反射率与地表反射率差值逐渐减小，尤其蓝通道降低最明显、绿通道次之、红通道降低最少；相较于可见光波段，近红外波段大气校正前后地表反射率表现恰好相反，各传感器近红外波段地表反射率较表观反射率略有增加，且部分数值增加较多，算法适应性还需进一步研究。大气校正处理对蓝、绿通道影响较大，这与波长较短的蓝、绿光易受大气影响，长波辐射对大气的穿透性较强的规律相一致。对于受大气影响较大的蓝、绿通道，大气校正处理使各传感器地表反射率相较于表观反射率值降低了6.12%～8.2%和2.6%～4.29%，红波段校正量为0.1%～2.18%，由此可见，数据定量应用时缺乏精确的大气校正处理必将引起相应后续误差。

表2 GF-1 WFV/MSS校正量

从2017年6月5日GF-1 MSS数据中采集裸土、农田、林地、建筑、淡水、海水区域大气校正前后试验影像样本像元每种300个，分析大气校正处理对不同地物类型影像光谱影响，比较校正前表观反射率和校正后地表反射率差异（图7）。

对于裸土、农田、林地、建筑4种地表下垫面类型，在可见光波段均表现出随波长增加大气校正处理使反射率差值差异减小，尤其在红波段大气校正处理与否反射率值变化不大的特征，这符合可见光短波波段易受大气影响的规律。无论是否进行大气校正处理，裸土及建筑可见光波段反射率均高于农田及林地、近红外波段反射率均低于农田及林地。大气校正前后，反射率随波长变化曲线趋势较为相似，无植被覆盖的裸土及建筑大气校正前各可见光波段反射率值相近，大气校正后反射率随波长逐渐增加；农田及林地大气校正前可见光波段反射率随波长增加减小，大气校正使波长短的蓝波段反射率明显降低；4个陆表地类近红外波段大气校正后反射率值均略增加。水体反射率整体表现出与陆表相反的特征，两类水体大气校正前可见光波段反射率特征相似，均表现为随波长降低，且海水反射率较淡水偏高，大气校正处理使水体各波段反射率值不同程度降低。

在LaSRC大气校正流程算法中，LC8 OLI气溶胶的模拟应用到了近红外和短波红外波段归一化指数。GF-1 WFV/MSS没有短波红外通道，在获取GF-1 WFV/MSS气溶胶时，利用归一化植被指数进行了替代，虽然NDVI与NDVImir具有较好的线性关系，但在稀疏植被覆盖地区差异大于植被覆盖地区。利用前述采集的300个裸土和建筑大气校正结果样本数据代表稀疏植被覆盖地区，以农田和林地代表植被覆盖地区，以同时次对应样本点LC8 OLI大气校正结果为基准，比较稀疏/有植被覆盖地区大气校正效果。

从表3可以看出，对于裸土和建筑样本数据GF-1 MSS大气校正结果与LC8 OLI大气校正结果相关系数分别为0.765～0.860和0.692～0.833，农田和林地样本数据大气校正结果相关系数分别为0.810～0.928和0.850～0.906。该方法在有农田和林木等植被覆盖的中低纬度大气校正效果较好，对稀疏植被的裸地和建筑地表大气校正效果相对稍差。

表3 稀疏/有植被覆盖地区大气校正结果比较

3.3 大气校正精度分析

Vermote等[33]已验证LC8 OLI大气校正结果与MODIS产品相比具有较高精度，而MODIS地表反射率产品普遍认为较为准确。本文以LC8 OLI数据大气校正结果为基准，比较同日同区域GF-1 WFV/MSS各传感器大气校正结果准确性。所用数据为2016 年8月26日WFV2/MSS1/MSS2、2017年3月31日WFV1、2017年6月5日WFV4、2017年9月7日WFV3，这些数据涉及所有GF-1传感器，且GF-1与LC8过境时传感器采集数据有较大交叉覆盖范围便于二者进行比较分析。

将同日GF-1与LC8数据进行地理配准，对两景影像分别出现的云、云阴影区域进行合并剔除，剪裁两景影像交叉部位。根据通道波段范围及中心波长，以GF-1 WFV/MSS通道1～4顺序对应LC8 OLI通道2～5，即将2种卫星传感器蓝、绿、红、近红外波段两两对应，提取GF-1 和LC8两景影像交叉区域所有对应像元，比较GF-1和LC8各景地表反射率算数平均值（）和二者差值（），逐像元平均绝对误差（）、相关系数（）进行大气校正结果的比较。

式中1、2分别为GF-1各传感器及LC8 OLI每通道整景影像表观反射率（TOA）或地表反射率（SR）平均值；1i、2i分别为GF-1各传感器及LC8 OLI各景每像元大气校正结果；为样本容量。

所选取的每种传感器的数据采集时间及下垫面不同，表4中GF-1及LC8各传感器大气校正前各通道表观反射率平均值（GF1,TOA与LC8,TOA）相差较多。但无论各传感器各通道平均大气校正结果还是逐像元大气校正结果比较，GF-1各传感器与LC8 OLI地表反射率差异均较小。大气校正后GF1各传感器每通道地表反射率平均值（GF1,SR）与LC8 OLI对应通道地表反射率平均值（LC8,SR）蓝通道相差0.13%～1.38%、绿通道相差为0.23%～1.84%、红通道相差为0.71%～1.51%、近红外通道相差为0.97%～6.64%；GF-1和LC8大气校正结果逐像元绝对误差，蓝通道为0.81%～1.72%、绿通道为0.81%～1.94%、红通道为0.95%～2.26%、近红外通道为1.83%～4.39%；其中除WFV4近红外通道校正结果相差稍多，其他误差均低于3%。

表4 GF-1与LC8大气校正结果比较

由于比较时进行了质检剔除了云覆盖数据，无损采集GF-1与LC8同日同地交叉地域覆盖的其他所有有效像元，样本容量都超过了千万级别，GF-1各传感器与LC8 OLI逐像元大气校正结果都具有较高相关性，且均通过0.05显著性检验。WFV1与OLI大气校正结果相关系数为0.825～0.963、WFV2与OLI为0.909～0.944、WFV3与OLI为0.905～0.956、WFV4与OLI为0.951～0.972、MSS1与OLI为0.860～0.900、MSS2与OLI为0.875～0.920，总体上相关系数与波长呈正相关，这是由于波长越长对大气穿透性越好，所需进行的校正相对较少，相应的引入误差的可能性越小。相较于WFV系列传感器，MSS 2个传感器大气校正结果与OLI相比相关性较低，是由于MSS空间分辨率8 m在像元重定位空间匹配时存在一定的偏差。从比较结果分析，对于GF-1 WFV/MSS传感器应用LC8 OLI大气校正算法进行大气校正处理是适用的。

4 结 论

1）本文验证了本地化LaSRC流程进行LC8 OLI数据大气校正结果与USGS LC8 OLI地表反射率产品通道1～7决定系数分别为0.997、0.998、0.998、0.997、0.995、0.996、0.993，即相关性较好的基础上，通过对逐像元几何参数的计算、逐像元气溶胶光学厚度的求算、MODIS产品的调用等系统建设，构建了GF-1 WFV/MSS数据大气校正方法流程及应用平台。

2）通过对试验数据进行大气校正处理，GF-1数据大气校正后目视判读性得到改善。研究表明，大气校正使蓝、绿通道地表反射率有所降低，大气校正后的地表反射率能更好反映下垫面情况；陆表地类大气校正结果随下垫面不同表现出差异，GF-1 MSS大气校正结果与LC8 OLI大气校正结果相关系数在裸土和建筑样本区域分别为0.765～0.860和0.692～0.833，农田和林地样本区域分别为0.810～0.928和0.850～0.906。

3）应用同日同地交叉地域GF-1 WFV/MSS与LC8 OLI大气校正结果分析表明，各通道地表反射率差值和逐像元绝对误差较小；二者相关系数WFV1与OLI为0.825～0.963、WFV2与OLI为0.909～0.944、WFV3与OLI为0.905～0.956、WFV4与OLI为0.951～0.972、MSS1与OLI为0.860～0.900、MSS2与OLI为0.875～0.920，相关系数最低为0.825，相关关系较好。因此，对于GF-1 WFV/MSS数据应用LaSRC流程进行大气校正处理能够获得较为可观的精度。

4）与LC8 OLI相比GF-1缺少短波红外通道，为继承LC8 OLI大气校正过程中气溶胶光学厚度通过MODIS产品获取待校正数据像元水平的做法，GF-1数据大气校正处理时，求解气溶胶时用NDVI线性斜率和截距代表过程变量NDVImir，实现了用待校正GF-1数据获取其气溶胶光学厚度的方法。近红外波段大气校正前后地表反射率表现恰好相反，各传感器近红外波段地表反射率较表观反射率略有增加，且部分数值增加较多，并且在GF-1 WFV/MSS与LC8 OLI大气校正结果比较时，近红外通道校正结果也表现出偏差较大，由于LC8 OLI近红外波段光谱响应范围较GF-1相比偏窄，虽然利用光谱匹配系数对一些参数做了调整，难免还存在偏差，后续还需进行算法的改进及参数调整。

[1]邱鹏勋，汪小钦，茶明星，等. 基于TWDTW的时间序列GF-1 WFV农作物分类[J]. 中国农业科学，2019，52(17)：2951－2961.

Qiu Pengxun, Wang Xiaoqian, Cha Mingxing, et al. Crop identification based on TWDTW method and time series GF-1 WFV[J]. Scientia Agricultura Sinica, 2019, 52(17): 2951－2961. (in Chinese with English abstract)

[2]常布辉，王军涛，罗玉丽，等. 河套灌区沈乌灌域GF-1/WFV遥感耕地提取[J]. 农业工程学报，2017，33(23)：188－195.

Chang Buhui, Wang Juntao, Luo Yüli, et al. Cultivated land extraction based on GF-1/WFV remote sensing in Shenwu irrigation area of Hetao Irrigation District[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(23): 188－195. (in Chinese with English abstract)

[3]Bernardo N, Watanabe F, Rodrigues F, et al. An investigation into the effectiveness of relative and absolute atmospheric correction for retrieval the TSM concentration in inland waters[J]. Modeling Earth Systems and Environment, 2016, 2(3): 1－7.

[4]赵英时. 遥感应用分析原理与方法[M]. 北京：科学出版社，2003.

[5]Richter R. A spatially adaptive fast atmospheric correction algorithm[J]. International Journal of Remote Sensing, 1996, 17(6): 1201－1214.

[6]Chavez P S. An improved dark-object subtraction technique for atmospheric scattering correction of multispectral data[J]. Remote Sensing of Environment, 1988, 24(3): 459－479.

[7]Kaufman Y J, Wald A E, Remer L A, et al. The MODIS 2.1-m channel-correlation with visible reflectance for use in remote sensing of aerosol[J]. IEEE Transactions on Geoscience and Remote Sensing, 1997, 35(5): 1286－1298.

[8]Hall F G, Strebel D E, Nickeson J E, et al. Radiometric rectification: Toward a common radiometric response among multi-date, multisensor images[J]. Remote Sensing of Environment, 1991, 35(1): 11－27.

[9]Kruse F A. Use of airborne imaging spectrometer data to map minerals associated with hydrothermally altered rocks in the northern grapevine mountains, Nevada and California[J]. Remote Sensing of Environment, 1988, 24(1): 31－51.

[10]Nieuwenhove V V, Beenhouwer J D, Carlo F D, et al. Dynamic intensity normalization using eigen flat field in X-ray imaging[J]. Optics Express, 2005, 23(21): 27975－27989.

[11]Smith G M, Milton E J. The use of the empirical line method to calibrate remotely sensed data to reflectance[J]. International Journal of Remote Sensing, 1999, 20(13): 2653－2662.

[12]Tanré D, Deroo C, Duhaut P, et al. Description of a computer code to simulate the satellite signal in the solar spectrum: the 5S code[J]. International Journal of Remote Sensing, 1990, 11(4): 659－668.

[13]Tanré D, Herman M, Deschamps P Y, et al. Atmospheric modeling for space measurements of ground reflectances, including bidirectional properties[J]. Applied Optics, 1979, 18(21): 3587－3594.

[14]Kotchenova S Y, Vermote E F, Matarrese R, et al. Validation of a vector version of the 6S radiative transfer code for atmospheric correction of satellite data. Part I: Path radiance[J]. Applied Optics, 2006, 45(26): 6762－6774.

[15]Kotchenova S Y, Vermote E F, Levy R, et al. Radiative transfer codes for atmospheric correction and aerosol retrieval: intercomparision study[J]. Applied Optics, 2008, 47(13): 2215－2226.

[16]刘佳，王利民，杨玲波，等. 基于6S模型的GF-1卫星影像大气校正及效果[J]. 农业工程学报，2015，31(19)：159－168. Liu Jia, Wang Limin, Yang Lingbo, et al. GF-1 satelliate image atmospheric correction based on 6S model and its effcet[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(19): 159－168. (in Chinese with English abstract)

[17]孙林，于会泳，傅俏燕，等. 地表反射率产品支持的GF-1PMS气溶胶光学厚度反演及大气校正[J]. 遥感学报，2016，20(6)：216－228. Sun Lin, Yu Huiyong, Fu Qiaoyan, et al. Aerosol optical depth retrieval and atmospheric correction application for GF-1 PMS supported by land surface reflectance data[J]. Journal of Remote Sensing, 2016, 20(6): 216－228. (in Chinese with English abstract)

[18]孙元亨，秦其明，任华忠，等. GF-4/PMS与GF-1/WFV两种传感器地表反射率及NDVI一致性分析[J]. 农业工程学报，2017，33(9) ：167－173. Sun Yuanheng, Qin Qiming, Ren Huazhong, et al. Consistency analysis of surface reflectance and NDVI between GF-4/PMS and GF-1/WFV[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(9): 167－173. (in Chinese with English abstract)

[19]薛现光，周扬，胡校飞，等. 基于大气参数拟合的GF-1多光谱影像大气校正[J]. 海洋测绘，2018，38(5)：50－54.

Xue Xianguang, Zhou Yang, Hu Xiaofei, et al. Atmospheric correction of GF-1 multi-spectral image based on fitting of atmospheric parameters[J]. Hydrographic Surveying and Charting, 2018, 38(5): 50－54. (in Chinese with English abstract)

[20]胡勇，仲波，马泽忠，等. 采用MODIS大气产品的高分一号WFV数据大气校正[J]. 遥感信息，2018，33(5)：35－40.

Hu Yong, Zhong Bo, Ma Zezhong, et al. GF-1 WFV atmospheric correction based on MODIS atmosphere produce and 6S model[J]. Remote Sensing Informaton, 2018, 33(5): 35－40. (in Chinese with English abstract)

[21]吴北婴. 大气辐射传输实用算法[M]. 北京：北京气象出版社，1998.

[22]Berk A, Bernstein L S, Anderson G P, et al. MODTRAN cloud and multiple scattering upgrades with application to AVIRIS[J]. Remote Sensing of Environment, 1998, 65(3): 367－375.

[23]Callieco F, Dell-Acqua F. A comparison between two radiative transfer models for atmospheric correction over a wide range of wavelengths[J]. International Journal of Remote Sensing, 2011, 32(5): 1357－1370.

[24]杨贵军，黄文江，刘三超，等. 环境减灾卫星高光谱数据大气校正模型及验证[J] . 北京大学学报：自然科学版，2010，46(5)：821－828.Yang Guijun, Huang Wenjiang, Liu Sanchao, et al. Research on modeling and validating of atmospheric correction for HJ-1A hyperspectral imager data[J]. Acta Scientiarum Naturalium Universitatis Pekinensis, 2010, 46(5): 821－828. (in Chinese with English abstract)

[25]郑盛，赵祥，张颢，等. HJ-1卫星CCD数据的大气校正及其效果分析[J] . 遥感学报，2011，15(4)：709－721. Zheng Sheng, Zhao Xiang, Zhang Hao, et al. Atmospheric correction on CCD data of HJ-1 satellite and analysis of its effect[J]. Journal of Remote Sensing, 2011, 15(4): 709－721. (in Chinese with English abstract)

[26]熊世为，沈安云，李卫国，等. MODTRAN模型在HJ/CCD影像大气校中的应用[J] . 江苏农业学报，2016，32(2)：319－324. Xiong Shiwei, Shen Anyun, Li Weiguo, et al. Application of MODTRAN model in atmospheric correction of HJ/CCD data[J]. Jiangsu Journal of Agricultural, 2016, 32(2): 319－324. (in Chinese with English abstract)

[27]Cooley T, Anderson G P, Felde G W, et al. FLAASH, a MODTRAN4-based atmospheric correction algorithm ,its application and validation[J]. IEEE International Geosciences and Remote Sensing Symposium, 2002(3): 1414－1418.

[28]陈玲，陈理，李伟，等. 基于FLAASH模型的Worldview3大气校正[J]. 国土资源遥感，2019，31(4)：26－31. Chen Ling, Chen Li, Li Wei, et al. Atmospheric correction of Worldview3 image based on FLAASH model[J]. Remote Sensing for Land and Resource, 2019, 31(4): 26－31. (in Chinese with English abstract)

[29]李卫国，蒋楠，王纪华. 基于薄云雾去除的ETM+影像大气校正[J] . 农业工程学报，2013，29(增刊1)：82－88. Li Weiguo, Jiang Nan, Wang Jihua. Atmospheric correction for ETM+ image based on thin cloud removal[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2013, 29(Supp.1): 82－88. (in Chinese with English abstract)

[30]Richter R. A fast atmospheric correction algorithm applied to Landsat TM images[J]. International Journal of Remote Sensing, 1990, 11(1): 159－166.

[31]王正海，段建军，耿欣. 基于波普匹配的Hyperion数据大气校正方法对比研究[J] . 遥感技术与应用，2011，26(4)：432－436. Wang Zhenghai, Duan Jianjun, Geng Xin, et al. Comparative study atmospheric correction methods of Hyperion data based on spectral matching[J]. Remote Sensing for Land and Resource, 2011, 26(4): 432－436. (in Chinese with English abstract)

[32]Vermote E F, Justice C, Claverie M, et al. Preliminary analysis of the performance of the Landsat 8/OLI land surface reflectance product[J]. Remote Sensing of Environment, 2016, 185: 46－56.

[33]Vermote E F, Tanré D, Deuzé J L, et al. Second simulation of the satellite signal in the solar spectrum, 6S: An overview[J]. IEEE Geoscience and Remote Sensing Society, 1997, 35(3): 675－686.

[34]Vermote E F. Landsat8/OLI and sentinel 2 and VIIRS surface reflectance is largely based on MODIS C6(LaSRC) [EB/OL]. Workshop on Uncertainties in Remote Sensing ESA/ESRIN Frascati Italy [2017-10-25]. https://earth.esa.int/documents/ 700255/3264428/3. Vermote_ESRIN-WS-Uncertainties.pdf.

[35]Steven M D, Malthus T J, Baret F, et al. Intercalibration of vegetation indices from different sensor systems[J]. Remote Sensing of Environment, 2003, 88(4): 412－422.

[36]Tan B, Masek J G, Wolfe R, et al. Improved forest change detection with terrain illumination corrected Landsat images[J]. Remote Sensing of Environment, 2013, 136: 469－483.

[37]Ångström A. The parameters of atmospheric turbidity[J]. Tellus, 1964, 16(1): 64－75.

[38]Dutton E G, Reddy P, Ryan S, et al. Features and effects of aerosol optical depth observed at Mauna Loa, Hawaii: 1982－1992[J]. Journal of Geophysical Research Atmospheres, 1994, 99(4): 8295－8306.

[39]Bodhaine B A, Norman B W, Ellsworth G D, et al. On Rayleigh optical depth calculations[J]. Journal of Atmospheric and Oceanic Technology, 1999, 16(11): 1854－1861.

Atmospheric correction method of GF-1 data based on Landsat8 product algorithm flow

Zhang Xiaoyue1,2,3, Li Linlin2,3, Wang Ying2,3, Zhang Qi2,3, Li Guochun4

(1.110166,; 2.110166,; 3.110166,; 4.,,110866,)

GF-1 satellite has the characteristics of high spatial resolution and short revisiting period, serving as the first satellite of the High-resolution Earth Observation System for National Science and Technology Major Project in China. The satellite carries two multi-spectral high-resolution cameras (panchromatic multispectral sensor, PMS) and four multi-spectral medium-resolution camera (wide field of view, WFV). GF-1 captured data play an important role in the identification of the underlying surface, and these data can be obtained free of charge from the website of the China Centre for Resources Satellite Data and Application. An important step for the application of GF-1 satellite data can be the interference removal of atmospheric molecules, aerosols, ozone, water vapor. The spectral response curves from Landsat8 (LC8) operational land imager(OLI) and GF-1 MSS/WFV were then analyzed in the visible and near the infrared bands. The results showed that the spectral range of LC8 OLI in the red- and near the infrared bands was relatively narrower than that of GF-1 MSS/WFV, whereas the spectral response function in the blue- and green bands was slightly different from that of GF-1 MSS/WFV, indicating that it is feasible to transplant LaSRC correction process to GF-1 MSS/WFV data. Since GF-1 satellite lacks the short-wave infrared band compared with LC8, the algorithm was modified to adapt to the characteristics of GF-1 channels. The atmospheric correction project was designed for the GF-1 satellite MSS/WFV data, including the algorithm analysis and code implement based on 6S atmospheric radiation simulation model and C++ programming language. Some parameters were used to estimate initial aerosol, including total atmospheric transmission, gaseous transmission, atmosphere spherical albedo and actual values of digital elevation model, atmospheric precipitation, ozone content. The loop calculation of the aerosol optical thickness(AOT) was carried out until the ratio between the red- and blue bands of GF-1 MSS/WFV data equal to the prescribed ratio of MODIS, according to the relationship between the blue- and the red surface reflectance known from MODIS. The results can be obtained the surface reflectance with the minimum residual error during different Ångström coefficients, and retrieved the aerosols in the pixel level of GF-1 MSS/WFV data. The pixel aerosol and these parameters were then substituted into 6S model to calculate the surface reflectance. Since the project was equipped a data file containing the atmospheric precipitation and ozone content at the current day, the surface reflectance could be obtained when only inputting GF-1 MSS/PMS data. Because the data of atmospheric influence gases, such as ozone and water vapor, on the same day of the data to be corrected were sometimes difficult to identify, two schemes can be provided, one is to use the data at that time, the other is to use the daily values of ozone and water vapor in past six years instead. The experimental results show that the proposed method has a good effect on the atmospheric correction in the middle and low latitudes that covered by vegetation, such as farmland and trees, but not good effect on that of the bare land and building surface that covered by sparse vegetation. Based on 6S model and LaSRC correction process, the correlation coefficient of the atmospheric correction between GF-1 MSS/WFV and LC8 OLI was from 0.825 to 0.972, indicating a high correlation of atmospheric correction results for two satellites. WFV similar spatial resolution to that of LC8 OLI was in good agreement with that of LC8 OLI atmospheric correction compared with that of MSS. The results show that it is convenient and operable for the GF-1 satellite data atmospheric correction method using the self-estimation aerosol parameters in the pixel level based on 6S model and LaSRC process. This promising atmospheric method can be very suitable for the land surface application, such as agricultural and forestry monitoring in growing season. At present, this method has been successfully implemented on Remote Sensing data processing platform Remote Sensing Desktop (RSD) in China.

satellites; remote sensing; GF-1 WFV/MSS; Landsat8 OLI; 6S; LaSRC; atmospheric correction

张晓月，李琳琳，王 莹，张 琪，李国春. 采用Landsat8产品算法流程的高分一号数据大气校正[J]. 农业工程学报，2020，36(1)：182－192.doi：10.11975/j.issn.1002-6819.2020.01.021 http://www.tcsae.org

Zhang Xiaoyue, Li Linlin, Wang Ying, Zhang Qi, Li Guochun. Atmospheric correction method of GF-1 data based on Landsat8 product algorithm flow[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(1): 182－192. (in Chinese with English abstract) doi：10.11975/j.issn.1002-6819.2020.01.021 http://www.tcsae.org

2019-07-02

2019-12-02

中国气象局沈阳大气环境研究所开放基金（2017SYIAE03）

张晓月，高级工程师，主要从事农业气象和遥感数据处理及应用研究。Email：zhangxiaoyue225@163.com

10.11975/j.issn.1002-6819.2020.01.021

TP701

A

1002-6819(2020)-01-0182-11