中国陆地区域GF-1 WFV地表反射率产品

佘文庆，张兆明，3，彭燕，何国金，3，龙腾飞，3，王桂周

1.中国科学院空天信息创新研究院,北京 100094;2.中国科学院大学,北京 100049;3.海南空天信息研究院海南省地球观测重点实验室,三亚 572029

1 引言

高分一号卫星（以下简称GF-1）是中国高分辨率对地观测系统的首颗卫星，于2013年4月26日在酒泉卫星发射中心成功发射。GF-1 卫星搭载两台2 m分辨率的全色、8 m分辨率的多光谱相机和4台16 m分辨率的宽幅相机。其中4台宽幅相机影像拼接合成后，幅宽可达到800 km（白照广，2013）。目前GF-1 卫星数据已广泛应用于林业、农业和生态环境等领域，在一定程度上弥补了中国高分辨率卫星数据的应用短板。

地表反射率是表征物体表面对太阳辐射反射的能力，是许多生态遥感参数反演的输入数据，例如生物量和叶面积指数等。目前国内外公开共享的地表反射率产品主要有250 m、500 m、1 km分辨率的MODIS地表反射率产品（Liang等，2002；Vermote 等，2002）、30 m 分辨率的Landsat 系列地表反射率产品（Roy 等，2010；Vermote 等，2016；彭燕 等，2020）和10 m、20 m、60 m 分辨率的Sentinel-2地表反射率产品（Main-Knorn等，2017）。国产GF-1 WFV数据已面向公众开放共享，但是仅提供L1 级标准化产品。如果能在GF-1 WFV L1 级数据的基础上生产16 m 分辨率的地表反射率产品，即可与Landsat系列、Sentinel系列的地表反射率产品协同形成中国区域近每日的密集时间序列中高分辨率地表反射率产品。

地表反射率反演的关键步骤是大气校正，大气校正方法主要分为统计模型和物理模型（李正强 等，2018）。统计模型方法是利用地表变量和遥感数据之间建立的统计关系模型进行大气校正，但由于区域之间存在差异性，统计模型通常适用于局部地区（郑伟和曾志远，2004；申茜 等，2021）。物理模型方法的基本原理是利用电磁波在大气中的辐射传输原理建立起来的模型对遥感影像进行大气校正，其大气校正精度高且应用范围广（郑伟和曾志远，2004）。基于物理模型的遥感影像大气校正的一个核心问题是如何获取成像时的气溶胶光学厚度AOD（Aerosol Optical Depth）、大气水蒸汽含量WV（Water Vapor）和臭氧含量TO（Total Ozone）等大气参数。通常的做法是采用相关的大气产品或利用自身数据进行反演。如Hu 等（2014）和Peng 等（2016）分别将MOD04（AOD）、MOD05（WV）和MOD07（TO）应用于Landsat 5和Landsat 8 影像的大气校正，地面实测地表反射率数据验证表明了算法的可行性。Martins等（2018）将MCD19A2（AOD 和WV）和MOD08D3（TO）应用于CBERS-4 MUX 影像的大气校正，生产南美洲的CBERS-4 MUX 地表反射率产品，利用Landsat 8 OLI 地表反射率和AERONET 站点数据验证了产品的精度，验证结果表明产品具有较高的可靠性。美国WELD（Web Enabled Landsat Data）团队利用MODIS 大气产品（MOD04、MOD05 和MOD07）生产美国区域的Landsat 地表反射率产品（Roy 等，2010）。这些研究表明了将相关的MODIS大气产品应用于大区域遥感影像大气校正的可行性。在大气参数反演方面，大气水蒸汽含量反演需要数据具有水汽吸收通道（Main-Knorn 等，2017），不适用于GF-1 WFV 数据；气溶胶光学厚度反演常用的算法为暗目标（Dark Target，DT）法，该算法是利用短波红外波段受大气影响较小的特性，根据暗目标区域短波红外波段与蓝、红波段地表反射率之间的经验线性关系来反演AOD（Kaufman 和Sendra，1988），但由于GF-1 WFV 数据缺少短波红外波段，导致DT 方法应用于GF-1 WFV 数据存在较大困难（Richter等，2006）。

因此利用MODIS 大气产品辅助GF-1 WFV 大气校正成为了一个重要途径。然而MODIS 大气产品在部分区域存在无效值，尤其是AOD 产品，在一定程度上会影响地表反射率反演的精度。针对此问题，通常的做法是直接对单一AOD 产品采用传统克里金法插值填充无效值区域（Hu等，2014；胡勇 等，2018）。该方法虽然能在一定程度上解决无效值问题，但是其计算复杂，且对大量无效值区域的插值会出现数据失真的现象。目前已发布和共享了多种不同分辨率的MODIS AOD 产品，若能将这些多源产品在空间上进行融合，就可以减少大范围无效值的情况。

本文在参考相关工作的基础上，开展多源MODIS AOD 产品协同的GF-1 WFV 大气校正方法研究。基于多源MODIS AOD 空间融合和动态查找表，设计和实现一种针对GF-1 WFV 的6S 大气校正算法，并利用该算法生产和共享中国陆地区域2020 年地表反射率产品。该产品可为全国植被生长监测、地表覆盖分类和重要生态遥感参数反演等研究提供基础数据支撑。

2 数据源介绍

2.1 GF-1 WFV数据

GF-1 卫星传感器由2 台多光谱相机（PMS）和4台宽视场相机（WFV）组成。本文所用的数据为GF-1 WFV，其参数如下表1 所示。本文选取2020 年覆盖中国陆地区域的无云或者少云的GF-1 WFV 正射数据，总共850 景。所有数据来源于中国资源卫星应用中心（http：//www.cresda.com/［2022-10-22］）和对地观测数据共享计划（http：//ids.ceode.ac.cn/［2022-10-22］）。

表1 GF-1 WFV数据参数Table 1 GF-1 WFV data parameters

2.2 MODIS大气产品

GF-1 WFV 大气校正所使用的MODIS AOD 产品包括MOD08D3、MOD09CMA 和MCD19A2。其中MOD08D3提供3种全球逐日1° AOD产品（Wei等，2019），本文采用DTB（DT和DB（Deep Blue，深蓝）结合）算法获取的AOD；MOD09CMA 提供全球逐日0.05° AOD 产品（Vermote 和Kotchenova，2008）；MCD19A2提供全球逐日1 km AOD（Lyapustin等，2018）。大气水蒸汽含量和臭氧含量分别来自MCD19A2 和MOD09CMG 产品。上述产品均可在MODIS产品网站下载和使用（https：//ladsweb.modaps.eosdis.nasa.gov/archive/allData/［2022-10-22］）。

2.3 其他数据

其他数据包括全球30 m 分辨率 ASTER GDEM（Abrams 等，2020）和重采样到2.5 nm 步长的GF-1 WFV光谱响应函数。

3 数据处理算法

中国陆地区域2020 年GF-1 WFV 地表反射率产品的生产技术流程如图1 所示，主要包括GF-1 WFV 预处理、多源MODIS AOD 产品空间融合、动态查找表构建和快速大气校正等步骤。

图1 GF-1 WFV 地表反射率产品生产技术流程图Fig.1 Flow chart of GF-1 WFV surface reflectance product generation

3.1 GF-1 WFV数据预处理

在对GF-1 WFV L1 级数据进行正射校正的基础上，利用辐射定标系数和式（1）将卫星影像DN值转换为辐射亮度值，然后利用式（2）计算表观反射率TOA（Top of Atmosphere reflectance）。

式中，Le为辐射亮度，Gain和Bias分别为辐射定标系数增量和偏移量。本文中的辐射定标系数来源于包头辐射定标场获取的数据（Ma等，2020）。

式中，ρTOA为表观反射率、Le是辐射亮度、d为日地相对距离，近似取值为1（刘佳 等，2015）、ESUN为大气外层太阳辐照度，θ为太阳天顶角。

3.2 多源MODIS AOD产品空间融合

多源MODIS AOD 产品的空间融合通常是在以确定某一AOD 产品为主数据，其他AOD 产品为副数据的基础上，可采用直接替换法、线性替换法或泛克里金法UK（Universal Kriging）进行替换。直接替换法针对主数据未覆盖而副数据覆盖的地区，直接用副数据的值进行替换。线性替换法针对主数据和副数据的相同覆盖区域，建立线性回归模型；对主数据未覆盖而副数据覆盖的区域，利用线性回归模型计算的值进行替换。UK 方法是普通克里金法的改进版本，它是一种考虑存在漂移的无偏线性估计的地统计方法（Armstrong，1984）。UK 法替换针对主数据和副数据的相同覆盖区域，在泛克里金的漂移项内考虑副数据的影响因素，构建泛克里金模型。对于主数据未覆盖而副数据覆盖的地区，利用泛克里金模型计算的值进行替换。Chatterjee 等（2010）曾利用UK方法将MODIS AOD 和MISR AOD 数据进行融合研究，表明了UK方法在AOD融合上的可行性和有效性。

为了验证不同融合方法的效果，以MCD19A2多个AOD 产品的融合为例进行实验。MCD19A2 AOD 产品在一天内包含多个时间的产品。将产品按上午和下午时间段划分并取均值，得到上午数据和下午数据。以上午数据为主数据，下午数据为副数据，分别采用直接法、线性法和泛克里金法进行替换，替换结果如图2所示。结果显示由于气溶胶随时间变化的影响较大，直接替换会造成融合数据出现不连续和一致性差的问题（图2（c））。线性替换由于考虑了两者的线性关系，一致性得到了一定的提升，但是融合后的数据在局部地区依然会出现连续性和一致性较差的问题（图2（d））。泛克里金替换通过结合副数据漂移影响和主数据插值，融合结果整体上具有较好的连续性和一致性（图2（e））。

图2 不同替换方法的效果比较Fig.2 Effects of different replacement methods

因此本文采用泛克里金替换，将MOD08D3、MOD09CMA 和MCD19A2 的AOD 数据进行空间融合，获得三者空间覆盖度的并集产品。首先将MCD19A2 上午和下午的AOD 数据空间融合得到MCD19A2（1 km）。然后将MOD08D3 和MOD09 CMA 不同分辨率的AOD 数据空间融合得到MOD0809（5 km）。最后将MCD19A2（1 km）和MOD0809（5 km）的AOD 数据再进行空间融合，得到MOD080919（1km）的AOD 数据，具体流程如图3所示。

3.3 动态查找表构建

6S（Second Simulation of the Satellite Signal in the Solar Spectrum）模型是由法国大气光学实验室Tanre D.和美国马里兰大学地理系Vermote E.在5S模型基础上发展而来（Vermote 等，1997）。6S 模型对大气辐射传输过程进行模拟和计算，得到用于大气校正的系数xa，xb和xc，利用式（3）可计算得到地表反射率。

式中，ρac为地表反射率，Le为辐射亮度，xa、xb和xc为大气校正系数。

查找表可近似替代真实模型，提高大气校正的效率。传统查找表的维度和步长是预先设定的，不同的维度和步长会影响查找表的精度和构建时间。传统查找表在维度上选择气溶胶光学厚度、高程、大气水蒸汽含量、太阳天顶角和卫星天顶角等（Peng 等，2016）。本文在参考传统查找表基础上，提出动态查找表方法。

通过对2013 年—2019 年4 万余景GF-1 WFV的星下点太阳和卫星天顶角进行统计，发现其星下点太阳天顶角主要在60°以内，星下点卫星天顶角均在50°以内（图4）。一景GF-1 WFV 影像星下点到影像边缘的太阳天顶角和卫星天顶角相差约1°和8°（王中挺 等，2015）。逐像元天顶角与星下点天顶角反演的地表反射率结果差异较小，因此以星下点天顶角代表整景影像的天顶角。

图4 天顶角统计结果Fig.4 Statistics of zenith angles

根据上述分析结果，动态查找表在维度上选取气溶胶光学厚度、高程和大气水蒸汽含量。

动态查找表的步长则是根据GF-1 WFV 影像的成像时间和空间范围提取对应的气溶胶光学厚度、高程和大气水蒸汽影像的数值范围，并分别划分为8 等分、6 等分和4 等分来动态确定。由于每景GF-1 WFV 影像的成像角度不同，根据每景影像的具体成像条件，均会构建一个独立的查找表。

3.4 快速大气校正

在上述GF-1 WFV 数据预处理、MODIS 大气产品预处理和动态查找表构建的基础上，实现了GF-1 WFV快速大气校正，具体流程如下：

（1）根据GF-1 WFV 影像提取对应成像时间和空间范围的MODIS WV 和AOD 产品及DEM 数据。利用泛克里金法将多种不同分辨率的MODIS AOD 产品进行空间融合，得到1 km 分辨率全覆盖的MODIS AOD产品。

（2）提取对应的MODIS WV、泛克里金融合的AOD 和DEM 数据的值域范围以及GF-1 WFV 影像星下点的天顶角，调用6S模型构建动态查找表。

（3）综合考虑大气校正算法的空间成本和时间效率，将大气校正辅助数据重采样到64 m。对所有数据分块，利用动态查找表生成64 m 分辨率的校正系数图像。利用线性插值将校正系数图像重采样到16 m，结合式（3），计算地表反射率。

（4）在多进程支持下，并行逐块执行步骤（3），最后将每一块结果进行拼接得到GF-1 WFV 的地表反射率图像。

在Python 开发环境中按照上述流程编写代码实现GF-1 WFV快速自动化地表反射率反演。

4 结果与分析

利用开发的GF-1 WFV 地表反射率反演算法生产2020 年中国陆地区域的GF-1 WFV 地表反射率产品，共计850景。为了验证产品的可靠性，分别采用大气校正前后影像对比、与已有地表反射率产品交叉验证以及实测数据验证进行产品质量分析和精度验证。

4.1 GF-1 WFV大气校正前后影像对比

将全国陆地区域GF-1 WFV 的表观反射率和地表反射率影像采用相同的拉伸比进行对比，结果如图5 所示。从图5 可以看出，大气校正后影像的清晰度和质量都得到了改善，影像具有更好的辐射一致性。

通过计算大气校正前后的影像直方图，发现对于蓝、绿和红波段，地表反射率的峰值明显低于表观反射率峰值（图6）。这是由于对于蓝、绿和红等较短波长，大气影响主要表现为气溶胶散射，经过大气校正，消除了气溶胶散射的影响，地表反射率直方图峰值整体“左移”。从“左移”的幅度来看，由于蓝光波长最短，受大气散射影响最大，“左移”的幅度最大，其次是绿光，红光波段“左移”的幅度最小。

图6 GF-1 WFV大气校正前后的反射率直方图对比Fig.6 Comparison of GF-1 WFV reflectance histograms before and after atmospheric correction

在GF-1 WFV 影像内采集林地、农田、海水、湖水、裸土和建筑6种典型地物，比较大气校正前后的反射率，如图7所示。经过大气校正，各地物的反射率曲线与标准地物反射率曲线更为一致，例如植被在绿波段有小的反射峰，在红和蓝波段有吸收谷。同时典型地物可见光波段反射率的差异随波长的变短而增大，这也符合较短波长波段更容易受到大气影响的规律。

图7 GF-1 WFV典型地物大气校正前后反射率差异Fig.7 Reflectance difference of typical land cover with atmospheric correction of GF-1 WFV

4.2 GF-1 WFV 与Landsat 8 OLI交叉精度验证

采用USGS 的Landsat 8 OLI Collection1 Level2地表反射率产品（https：//earthexplorer.usgs.gov/［2022-10-22］）作为参考，对GF-1 WFV地表反射率产品开展交叉精度验证。

选取同步过境且覆盖同一地理区域的Landsat 8 OLI 和GF-1 WFV 影像，在影像重叠区域目视解译得到验证样本，样本类型包含植被、裸土和水体（图8），得到40 组Landsat 8 OLI和GF-1 WFV 验证数据（表2）。

图8 交叉验证样区样本空间分布Fig.8 Spatial distribution of validation samples in cross-validation

表2 交叉验证数据对应关系表Table 2 Cross-validation data correspondence table

在每组验证数据的重叠区域目视解译得到约30 个验证样区，40 组验证数据共得到1206 个验证样区，其中裸土样区740 个、植被样区293 个、水体样区173个。计算每个验证样区的均值来代表整个验证样区。选择均方根误差（RMSE）和平均绝对误差（MAE）作为精度评估指标。

式中，n为验证样区总数，yi为Landsat 8 OLI 地表反射率，为GF-1 WFV地表反射率。

考虑到GF-1 WFV 数据观测角度与Landsat 8 OLI 观测角度存在一定差异，利用C 因子BRDF 校正方法（Roy 等，2016）对GF-1 WFV 数据进行BRDF 校正，减小两种卫星数据观测角度差异的影响。为减小两种卫星传感器光谱响应函数差异的影响，本文基于USGS splib07 的光谱测量数据（Kokaly 等，2017），利用光谱匹配因子（SBAF）（Chander 等，2013）将GF1 WFV 的地表反射率转换为等效的Landsat 8 OLI 地表反射率。最终交叉验证结果如图9所示。

图9 交叉验证结果Fig.9 Results of cross-validation

从图9可以看出，蓝、绿、红和近红外波段的平均RMSE分别为1.79%、2.33%、2.90%和2.83%，最大RMSE 不超过3%，且GF-1 WFV 和Landsat 8 OLI 各波段地表反射率的相关性较好，R2均大于0.9。交叉精度验证结果表明本文生产的GF-1 WFV地表反射率产品具有较好的可靠性。

4.3 实测数据验证

为进一步验证GF-1 WFV 地表反射率产品精度，采用ASD 光谱仪实际观测的地表反射率数据开展直接精度验证。实测数据共有6 组，包括2019年9月13日福建南平的水稻、2019年9月24日甘肃张掖的燕麦、2019 年10 月20 日湖南株洲的水体、2019 年9 月15 日福建平潭的沙滩、2019 年9 月20 日甘肃敦煌的戈壁滩和2019 年9 月22 日甘肃嘉峪关的戈壁滩。

首先将所有实测数据转化为GF-1 WFV 等效的地表反射率，然后挑选与实测数据同天或者相差3 天内的GF-1 WFV 数据，并提取实测数据对应的验证区域，计算验证区域内蓝、绿、红和近红外波段的地表反射率均值，最后与实测数据进行对比验证，验证结果如表3和图10所示。

图10 光谱曲线对比Fig.10 Comparison of spectral curves

表3 实测数据验证结果Table 3 Validation result of ground-measured data

从表3 和图10 可以看出，各实测点各波段的平均RMSE在1.21%—6.14%，各实测点的测量光谱曲线和GF-1 WFV的光谱曲线趋势基本一致。对于植被（水稻和燕麦），实测反射率和GF-1 WFV反演反射率在蓝、绿和红波段较为接近，但近红外波段的反演反射率低于实测反射率。对于水体，两者在蓝、绿和红波段的差异也较小，但在近红外波段反演反射率高于实测反射率。对于裸地（沙滩和戈壁滩），GF-1 WFV 蓝、绿、红和近红外波段的反演反射率与实测值均非常接近。直接验证结果表明本产品具有较好的精度和可靠性。但从验证结果中也可以发现，相较于蓝、绿和红波段，植被和水体在近红外波段的反演误差较大。这可能是由于清洁水体在近红外波段的反射率非常低（接近于0），而健康植被在近红外波段具有较高的反射率（超过50%），但是GF-1 WFV 传感器在近红外波段的波谱范围相对较宽，导致该波段的光谱探测精度和敏感性降低，对于较高或较低反射特性地物，反演误差较大。

5 产品简介

中国陆地区域2020 年GF-1 WFV 地表反射率产品共有850 景，总计1.9 TB。数据格式包括TIF、XML 和json。产品由地表反射率数据（*_AIRAC_MOD080919_SR.TIF）、表观反射率数据（*_TOA.TIF）、元数据信息（*.XML）、动态查找表（*.json）和区域DEM、TO、WV、AOD（*.tif）组成。其中地表反射率和表观反射率数据无量纲，除以10000得到真实的地表反射率和表观反射率。产品面向国内外免费共享，可在ftp：//bigrs-info.com/gf1wfv_sr/2020内下载和使用。

6 结论

GF-1 WFV 是目前应用最为广泛的国产卫星数据之一。本文利用多源MODIS AOD 融合产品和动态查找表，采用分块多进程的方法设计和实现GF-1 WFV 快速大气校正算法，并利用该算法生产和共享中国陆地区域2020 年GF-1 WFV 地表反射率产品。通过定性分析、Landsat 8 OLI地表反射率产品交叉验证和地面实际测量的地表反射率数据对产品进行精度验证和分析。验证结果表明，本文生产的GF-1 WFV 地表反射率产品的质量良好，与Landsat 8 OLI 地表反射率产品具有较高的一致性，基于实测地表反射率数据的直接精度验证结果表明各波段的均方根误差（Root Mean Square Error）位于1.21%—6.14%，产品具有较高的精度和可靠性。