基于GF-1WFV数据的16 m分辨率反照率产品算法与验证

陆彦蓉，李霞，杨凯祥，刘强，2，闻建光，李秀红，2

1.北京师范大学 全球变化与地球系统科学研究院,北京 100875;

2.北京师范大学 遥感科学国家重点实验室,北京 100875;

3.中国科学院空天信息创新研究院,北京 100094

1 引言

地表反照率定义为地表对反射太阳辐射的通量与入射的太阳辐射通量的比值，量化了大气与地表之间的辐射相互作用。它在陆地表面气候、水文和生物圈模型中起着至关重要的作用（Dickinson，1983；Dirmeyer 和Shukla，1994）。随着城市化进程的加快，在对植被覆盖、城市化、区域生态建设的研究中，都需要利用地表反照率作为指标因子或驱动因素（Liang等，2010）。

在不同的时空分辨率下，气候、地球化学、水文和天气预报模型都需要具有绝对精度为0.02—0.05 反射单位的地表反照率（Sellers 等，1994）。所以早期对于地表反照率就有研究，后续为了更好地服务于气候科学与建模领域，研究学者已经从运行卫星上生成了一致的高质量全球陆地反照率数据集，但是监测陆地表面需要高质量高精度数据提供支持，所以高分辨率地表反照率需求也在快速增加。

近年来，极端天气以及异常天气逐渐增多，也同时促使在环境与气候研究中对于高质量的观测数据提出了极大的需求，在此背景下，地表反照率产品得到了发展和较为广泛的应用。目前常用的全球范围地表反照率数据集均为中等空间分辨率，例如美国发布的MODIS 地表反照率产品（MCD43 系列）、欧洲发布的全球地表反照率产品（GLOBALBEDO）（Lewis 等，2011），以及中国生产和发布的GLASS 反照率产品等。但是由于对反照率产品的应用需求日益增长，如在区域尺度的生态水文研究中，现有产品仍然存在着时空分辨率不足的问题，尤其是空间分辨率不足以支持不同领域的研究。

近年来高分辨率卫星遥感技术迅速发展，中国也通过一系列高分卫星获取了大量的高分辨率遥感数据，但是由于高分辨率遥感数据缺乏多角度观测信息，大部分传感器的波段范围也有限，目前主流的地表反照率定量遥感算法无法简单地应用于高分辨率遥感数据，所以制约了反照率估算的精度。

许多研究学者及团队提出了不同的算法开发高分辨率地表反照率产品（Wang 等，2013；He 等，2018）。一类算法是在高分辨率反照率反演中结合中/低分辨率的遥感数据或BRDF 信息。Shuai 等（2014）提出一个算法是用MODIS 500 m 分辨率BRDF产品提供地表各向异性信息，利用30 m分辨率Landsat 地表反射率数据反演得到30 m 分辨率无雪反照率。该算法首先对Landsat 图像进行处理并聚合到500 m 分辨率；从MODIS 500 m 反照率产品MCD43A 中提取纯像元的BRDF 参数，然后由提取的BRDF参数计算光谱反照率和反射率的比；最后由Landsat 地面反射率、BRDF 参数和30 m 分辨率Landsat 分类图，计算黑、白空反照率及宽波段反照率。You 等（2015）提出了一种基于核驱动模型的土地覆盖线性BRDF 解混算法（LLBU），该算法将机载高光谱传感器CASI的反射率与MODIS的日反射率产品相结合，通过对MODIS反射率进行线性分解提取表征其各向异性的反射率分布函数，从CASI反射率中导出各向同性系数，联合计算反照率。Gao 等（2014）对MODIS 反照率产品的AMBRALS算法进行改进，使算法可以考虑地形的影响，并且能够用到从环境灾害监测小卫星HJ-1A/B CCD数据中反演的地表反照率。张虎等（2015）利用BRDF 原型和单方向反射率数据估算了地表反照率，研究以MODIS BRDF 产品统计得到的BRDF 原型为先验知识，结合高分辨率HJ-1 数据反演得到30 m 地表反照率，反演得到的反照率与地面实测反照率有较好的一致性。即使在中/低分辨率下，BRDF 估算也常常不稳定，并且其尺度效应非常复杂，这导致以上算法推广应用的困难。另一类算法是高分辨率数据直接与中低分辨率的反照率产品结合，例如He 等（2014）选择北美大陆作为研究区域，使用多尺度树融合3种不同尺度的反照率产品（MODIS 500 m 反照率产品，MISR 1.1 km 反照率产品和Landsat ETM+30 m 反照率产品），生成多尺度时空连续的反照率产品。但该算法忽略了低分辨率像元的空间响应，且MODIS和MISR 反照率产品存在大量缺失，影响了高分辨率反照率产品的纹理的保真度。

本文提出相关算法将500 m分辨率GLASS反照率产品与16 m 分辨率高分一号反照率进行融合，在保持GLASS 反照率产品绝对辐射精度的前提下提高了分辨率，丰富了细节纹理信息。由于该算法相对简单，且输入的高分一号数据和GLASS 反照率产品均为中国公开发布数据，因此具有业务化生产16 m分辨率反照率产品的潜力。

2 实验区域与数据

2.1 黑河流域概况

黑河流域是中国西北地区第二大内陆流域，位于欧亚大陆中部。全球独特的以水为纽带的“冰雪—冻土—森林—河流—湖泊—绿洲—荒漠”多元自然景观，使黑河流域成为寒区和旱区生态、水文陆表过程研究的理想区域。黑河流域地表类型丰富，四季变化分明，为高分辨率产品研究提供了理想区域。

“黑河综合遥感联合试验”主要由中国科学院寒区旱区环境与工程研究所、中国科学院空天信息创新研究院、北京师范大学地理科学学部和中国林业科学研究院资源信息研究所组织实施，由中国科学院西部行动计划（二期）项目“黑河流域遥感—地面观测同步试验与综合模拟平台建设”与国家重点基础研究发展计划（973 计划）项目“陆表生态环境要素主被动遥感协同反演理论与方法”共同资助，是在流域尺度上开展的，以水循环及与之密切联系的生态过程为主要研究对象的大型航空、卫星遥感与地面同步观测科学试验。总体目标是为发展流域科学积累基础数据；发展能够融合多源遥感观测的流域尺度陆面数据同化系统，为实现卫星遥感对流域的动态监测提供方法和范例；发展尺度转换方法，实现对地表生态和水文变量的主被动遥感协同反演。

通过在黑河流域的遥感联合试验，获取了一套高质量的、多尺度的、经过规范化处理并发布共享的流域综合数据集，该数据集包括大量航空遥感、卫星遥感、地基观测等数据。通过运用这些数据，研究学者已在水文气象观测、生物物理参数反演、积雪参数提取、尺度转换等方面都取得了重要的成果，特别是在一些参数遥感反演方法、估算精度方面获得重要进展（李新 等，2012a，2012b）。黑河遥感实验也促进了定量遥感的发展及其与陆地表层系统科学的深度结合（李新 等，2023）。黑河流域的观测环境以及地面数据对于本算法的研究具有良好的代表性，故本研究选择黑河流域为测试算法的研究区。图1为研究区土地覆盖图。

图1 黑河研究区土地覆盖及站点Fig.1 Land cover and sites in Heihe research area

2.2 数据源介绍

2.2.1 高分数据

GF-1 卫星是中国第一颗作为“国家高分辨率地球观测系统”发射的系列卫星，主要用于陆地监测、环境监测等，为地理测绘、海洋和气候气象观测、水利和林业资源监测等领域提供数据支持。GF-1 卫星搭载了两台2 m 分辨率全色/8 m 分辨率多光谱相机（PMS），4 台16 m 分辨率多光谱中分辨率宽幅相机（WFV）。

GF1-WFV 传感器幅宽高达800 km，其包含蓝、绿、红、近红外等4个波段，重访周期为两天，空间分辨率为16 m，数据实现了高空间分辨率和高时间分辨率的完美结合。在2019年11月16日，中国国家航天局推出“高分卫星16 m 数据共享服务平台”，宣布正式将中国高分16 m数据对国外开放共享，国内用户可通过中国资源卫星应用中心网站进行下载（http：//www.cresda.com/CN［/2021-03-15］），数据下载和使用服务方便。本文用GF1-WFV作为反演16 m反照率的主要数据源。

2.2.2 GLASS数据

500 m分辨率的GLASS反照率产品使以MODIS数据短波波段观测数据为输入。算法首先采用基于角度网格（Angular Bin）的直接反演算法对遥感观测值进行反演，获得反照率初级产品，然后采用基于先验知识的时空滤波（Statistics based Temporal Filtering）算法进一步加工处理得到最终产品（Liu等，2013；Qu等，2014；Liang等，2013）。基于地面观测资料的验证结果表明，GLASS 反照率产品精度与广泛使用的MCD43 产品相当，而且GLASS 反照率产品也具有时空连续无缺失的优点，因此GLASS 总体质量高，方便数据用户使用（Liu 等，2013；王立钊 等，2014）。然而，500 m分辨率反照率不足以满足区域尺度精细研究的需求，所以我们考虑使用GLASS 反照率数据产品作为均值信息，通过降尺度融合提高分辨率。

2.2.3 地面测量数据

“黑河综合遥感联合试验”获取的站点短波辐射分量测量数据，可用于反照率产品的验证。本文选择黑河上游和中游的8个站点（具体站点信息见表1）的上行和下行短波辐射数据，首先根据下行辐射异常小的一天为阴天的理论，设定阈值判断天气状况。阈值的选取首先是根据辐射值与太阳天顶角余弦的比值进行归一化，从归一化辐射值中统计出95%意义上的最大值，若其他归一化辐射值小于最大值的75%，则认为是阴天。剔除天气状况异常的数据后根据当地正午前后半小时内上、下行辐射的平均值计算正午反照率，用于反照率的验证。各站点周边的均匀程度并不足以支持500 m 像元产品的验证，但各站点仪器架设高度所测量的面积（Siegel 和Howell，1992）与16 m像元的面积是基本匹配的，测量数据具有局部代表性，故可以直接进行16 m 分辨率反照率产品的对比验证。

表1 研究使用的黑河流域观测站点信息Table 1 Information of observation stations in the Heihe River Basin used in the study

3 反照率产品算法介绍

算法基本原理是将GF-1 WFV 传感器16 m 分辨率图像的纹理信息与500 m分辨率GLASS反照率产品的均值信息进行定量的融合得到16 m 反照率产品，包括简化直接反演算法、降尺度融合算法两个主要步骤，除此之外还包括简单的数据预处理过程。具体流程如图2所示。

图2 产品反演流程图Fig.2 Flow chart of product inversion

3.1 数据预处理

本文获得的图像是Level-1 产品，该产品是辐射校正后的大气顶层辐射亮度，并附有一组有理多项式系数（RPC）以进行几何校正。预处理的第一步是根据RPC 将原始图像转换为UTM 投影，称为几何粗校正。几何粗校正通常包含高达数十米甚至数百米的定位误差。因此，处理的第二步是根据地面控制点进行更精确的几何校正，即几何精校正。在山区需要进行正射校正，也就是说，除了获得精确的地面控制点外，还需要校正由地形引起的图像局部变形。为此，采用分辨率为10 m并根据SRTM（Shuttle Radar Topography Mission）全球90 m DEM 产品进行几何正射校正的Sentinel-2-MSI（https：//scihub.copernicus.eu［/2021-03-15］）图像来构成研究区域的几何校正底图。这里使用我们团队开发的软件，软件采用区域分幅、并行计算、金字塔式最大相关匹配算法、校正函数、形变函数方法过程，自动将GF-1WFV 图像与底图匹配，并根据SRTM90 m DEM 进行正射校正。目视检验是为了确保校正后的GF-1WFV 图像与底图之间的匹配度差在1.5 像素以内，这是目视检验不匹配的阈值。如果目视检验失败，则必须手动剔除错误的地面控制点并重做正射校正，直到满足匹配要求。预处理的第三步是大气校正，假定气溶胶类型为大陆型，近地面的可视距离为30 km，使用6SV 代码获得大气光学参数和校正公式将大气层顶辐射亮度校正为地表反射率（Kotchenova 和Eric，2007；刘佳 等，2015）。

3.2 直接反演算法

直接反演算法是一种基于BRDF先验知识的地表宽波段反照率估计的半经验方法。该算法早期由Liang 等（1999）；Liang（2003）提出并应用于MODIS 数据估算地表反照率，但是，这两项研究均未考虑陆地表面的各向异性。Liang 等（2005）进一步改进了它，建立了分角度格网的大气层顶反射率与地表短波波段反照率之间的多元线性回归关系，以便校正陆表反射的各向异性，并被用于估算格陵兰冰盖的每日陆地表面反照率。中国生产的GLASS 系列产品中直接反演算法（梁顺林 等，2013）得到进一步发展并用于生产反照率产品（Qu 等，2014）。该算法也被陆续运用到了不同分辨率的遥感数据中，如Landsat数据（He等，2018）和中国的HJ-1CCD卫星数据（He等，2015）。

本文也沿用直接反演算法（梁顺林 等，2013）对数据进行处理，该算法的核心部分是推导出每个角元的经验回归系数的查找表，然后利用查找表中的经验回归系数计算黑、白空反照率。本文算法中与之前的区别就是根据GF 数据的波段设置和角度分布，生成了适合GF 数据的查找表。因为GF-1WFV传感器的波段范围都在1 μm以内，1 μm到3 μm的短波辐射范围内缺少观测，再加上定标和大气校正中存在的不确定性，直接反演的反照率定量精度不高，在此仅把它作为初级产品。

3.3 融合算法

由于GF 数据生产的产品为初级产品，精度不高，而GLASS 反照率产品的绝对精度较高，只是分辨率低，所以我们考虑利用16 m 反照率初级产品中的纹理信息与GLASS 500 m反照率产品的均值信息进行融合，得到最终的目标产品，来提高产品的精度及分辨率。适用于不同空间分辨率定量产品融合的算法有基于小波变换的算法（Sundar 等，2017）、SFIM（基于亮度调节的平滑滤波）算法（李存军 等，2004）以及最小二乘拟合（商荣 等，2015）等。考虑到16 m 和500 m 反照率产品之间的尺度变换是一个空间采样过程，而已有的融合算法往往对像元空间响应函数缺少定量的描述和考虑，造成融合结果的过度锐化或过度平滑，因此本文提出了基于像元空间响应的融合算法，以支持高分16 m和GLASS 500 m反照率的深度融合。

为了便于与500 m GLASS数据做融合，需要把16 m分辨率的反照率初级产品先聚合成为500 m分辨率反照率初级产品。在聚合过程中，考虑了像元的空间响应函数，因此聚合公式为：

式中，Yi为500 m 分辨率反照率初级产品，i=1，…，m，yj为16 m 分辨率反照率初级产品，j=1，…，n，Ωi为低分辨率像元（500 m）的空间响应函数确定的一个邻域，wi，j为高分辨率像元（16 m）yj在Yi中占的权重。这里参考彭菁菁等（2015）的研究，采用高斯椭圆函数来模拟反照率产品的空间响应函数，其空间响应函数如下：

式中，σ为高斯函数的标准差，取值为375 m。

记降尺度融合结果的16 m 分辨率反照率为：zj，j=1，…，n。期望它能够满足以下两个特性：

（1）16 m 分辨率反照率融合结果应尽可能接近500 m 的时空滤波反照率，因此得到反映均值信息第一组方程为（式（3））：

式中，Xi为500 m 分辨率GLASS反照率产品，Ωi为第i个500 m 像元的空间响应范围确定的邻域，z为16 m分辨率反照率。

（2）16 m分辨率融合结果的纹理尽量接近16 m分辨率初级产品的纹理，因此得到反映纹理信息的第二组方程（式（4））：

式中，ε1和ε2分别为均值和纹理信息的不确定因素。

要在严格意义上求解这两组方程，一般来说需要构造代价函数，用迭代优化方法使ε1和ε2极小化，估算最优的融合结果，其计算量非常大。尤其是在低分辨率像元空间响应函数有重叠的情况下，每一个高分辨率像元都会关联到多个低分辨率像元，因此对每一个高分辨率像元反照率的估算都需要考虑到邻域内多个低分辨率像元的综合影响，这使得求解过程成为复杂的全局优化问题，在工程上难以实现。

因此，本文采用一种简单直观的近似解法，认为GLASS 反照率与降低分辨率的初级产品的差值就是初级产品需要的修正量，然后对领域内不同低分辨率像元的修正量进行加权平均，用于修正领域中心点的高分辨率像元，简化算法过程，将综合影响尽可能缩小，如下式（式（5））：

式中，zj为最优估算融合结果，为所有对第j个16 m像元有响应的500 m像元构成的邻域。

根据数值模拟试验的结果，这种简单直观解法获得的图像是迭代优化方法获得图像的很好的近似，同时也将计算量尽可能的缩小，更具有可行性，所以我们在生产16 m 产品时使用式（5）计算。

4 产品验证与分析

4.1 实验区的16 m分辨率反照率产品生产及效果

选取黑河上游与中游2016 年—2017 年两年的500 m 分辨率GLASS 反照率产品和16 m 分辨率高分数据进行融合反演，得到16 m 反照率产品，并结合8个站点地面测量值进行验证。

为了直观说明16 m 分辨率反照率产品的效果，我们选择甘肃张掖市南荒漠地区2016年10月18日的500 m 分辨率GLASS 反照率产品与生产的16 m反照率产品作对比（图3）。图3（a）是张掖市光伏产业园实况图（图片出自新华网新闻无人机队）。张掖市位于河西走廊中部，太阳能资源丰富，近几年张掖市大力开发清洁型新能源，建成了南荒漠上规模较大的光伏产业园，蓝色太阳能板成了戈壁滩上的一道风景。但是光伏发电改变地表的能量分配，其潜在的气候和生态效应需要深入研究。图3（b）是GF-1 假彩色合成影像上光伏产业园的位置，产业园离黑河遥感站较近。图3（c）是光伏产业园周边的500 m 分辨率GLASS 反照率产品，中间区域的蓝色椭圆圈出了光伏产业园的大致位置，可以看出光伏产业园的地表反照率低于周边的地表。由于空间分辨率较低，地面上壮观的人类工程在图像上退化为一个个马赛克，因此500 m 分辨率产品在研究人类活动对环境的影响时显得不足。图3（d）是同样区域的16 m 反照率产品，其信息量远高于500 m 反照率产品，图中光伏产业园位置清晰可见，从图3 中也可以看到园区反照率较低既有太阳能电池板的贡献，也有该区域荒漠土壤本身就比较暗的结果。另外我们看到16 m 反照率产品值的变化范围和趋势与GLASS 反照率产品一致，算法的这个特性能够保证我们在不同尺度的能量平滑研究中获得一致的结果。

图3 甘肃张掖市光伏产业园区的不同分辨率的图像Fig.3 Different resolution images of the Photovoltaic industrial park in Zhangye，Gansul

4.2 基于地面测量数据的验证

以选定研究区内的8个站点的实测反照率为基准，与反演结果进行比较，验证最终反照率的精度。剔除缺失和受降水影响的地面数据，以年为单位做实测数据和遥感反演结果的时间序列图，如图4展示了垭口和张掖湿地站的结果，结合站点2016 年与2017 年结果，总体来说，在夏季融合反照率与实测数据有较好的一致性，且融合反照率比初级反照率更加贴合实测数据。由此可以看出，融合的结果比较理想，在初级反照率的基础上进一步提高了反演反照率的精度。但是，在春冬两季垭口、阿柔、大沙龙、景阳岭4个站点地面观测值波动较大，同时也出现个别初级反照率比融合反照率更接近地面实测值的情况，如图4（a）—（d）所示。这是由于这4个站点海拔较高，冬季积雪导致反照率出现大幅度的波动和异常高值。

图4 研究站点反照率实测值时间序列与反演反照率对比图Fig.4 comparison of the time series of measured albedo values at the research site and the inverted albedo

将选定的8 个站点2016 年—2017 年的数据做汇总，剔除掉其中数据缺失和受积雪、降水影响的数据，即反照率大于1 的数据或小于0 的数据。制作融合反照率、初级反照率与地面反照率散点图。验证结果显示，图5（a）散点图中的散点分布密集，均方根误差RMSE 为0.02439，R2为0.7028；图5（b）中散点较分散，RMSE 为0.05135，R2为0.1924。对比发现，融合反照率结果与地面实测值具有较好的一致性，而初级反照率结果与地面实测值的一致性较差。

图5 2016年份—2017年站点反照率散点图Fig.5 Site albedo scatter plot in 2016 and 2017

4.3 误差影响因素分析

在我们所提出的算法流程中，GF-1 WFV 数据的定标和大气校正的不确定性将影响初级反照率的估算，但对最终的融合结果几乎没有影响。为了进一步分析16 m 分辨率融合反照率产品误差的原因，对选定研究站点的GLASS 500 m反照率时间序列与地面观测反照率进行对比。从时间序列图（图6）可看出，高海拔站点（垭口站）反照率在春冬两季变化激烈，反映了高海拔地区冬季的反复降雪、融化过程；GLASS 反照率因为是时间序列滤波的产品，其实际时间分辨率较低，序列平滑，所以在春冬两季有较大不确定性。与高海拔站点相比，海拔较低的站点（张掖湿地站）时间序列变化平稳，与GLASS 反照率更为一致。另外，从台站观测反照率曲线的波动程度可以看出，地面反照率波动大的时段GLASS 反照率产品的误差偏大，这是因为台站观测反照率的波动往往是云雨等天气过程造成的，而云雨等天气影响了遥感数据的获取，在GLASS 产品中如果晴空遥感观测数据不足，则产品值为时空滤波的填充值，被标记为低质量产品。而融合的反照率也表现出对GLASS 反照率的依赖性。需要说明的是这里把地面站点观测与低分辨率的遥感产品简单对比，未进行尺度匹配，所以并不能作为GLASS 反照率产品精度验证，仅用于分析16 m 反照率融合产品的误差原因。

图6 研究站点GLASS500 m反照率与实测数据时间序列图Fig.6 Time series diagram of the 500 m resolution GLASS albedo and measured albedo at research site

以上分析表明，积雪对验证精度影响很大。16 m 分照率融合产品的精度很大程度上是由GLASS 反照率产品决定的，地面有雪覆盖时GLASS 反照率产品不确定性较大，就可能出现初级反照率比融合反照率更接近实测值的情况。另外，云对验证精度也会产生影响，在验证中尽可能取晴天数据进行验证，但在高海拔地区云的干扰有时难以避免，云的干扰会使数据出现缺失或异常情况。

5 结论

本文提出了基于GF-1WFV数据的16 m反照率产品算法，算法将高分数据16 m 分辨率的纹理信息与500 m分辨率GLASS反照率产品的均值信息进行定量的融合得到16 m 反照率。有效结合GLASS反照率产品辐射精度高和高分卫星数据空间分辨率高的优势，主要结论与创新有：

（1）结合直接反演算法与降尺度融合算法反演了16 m 反照率产品，其具有与500 m 分辨率的GLASS 反照率产品取值范围和趋势一致的特点。算法流程简单，效果稳定，可快速生成高质量、高精度的16 m分辨率反照率产品。

（2）通过地面实测数据的验证，发现16 m 分辨率反照率与实测反照率有较好的一致性，且在无雪季节融合反照率比初级反照率更接近实测值，在有雪季节融合产品的优势不明显。

（3）通过在张掖市光伏产业园区域的对比可知，最终16 m反照率产品空间分辨率远高于GLASS 500 m反照率产品，影像效果更佳。最终16 m反照率产品也更适用于研究人类活动对环境的影响。

本文所得到的16 m 反照率产品与500 m 分辨率GLASS 反照率产品能保持一致性，而且具有更高的分辨率和精度，是局部区域环境研究的友好数据源。由于研究时间、研究数据下载、研究区天气、海拔等问题，地面数据存在缺失或有无效数据后续研究中，要加强反照率产品在积雪区域的研究，提高数据的准确度，并进行更为严密和全面的精度验证和分析。