郑凯旋， 林兴稳， 闻建光， 郝大磊

1. 浙江师范大学 地理与环境科学学院, 金华 321004;

2. 中国科学院空天信息创新研究院 遥感科学国家重点实验室, 北京 100094;

3. 中国科学院大学, 北京 100049

1 引 言

地表反照率（Albedo）定义为地表向各个方向反射的太阳短波辐射（0.3—3.0 μm）与太阳入射总辐射的比值（Liang等，2010），是表征地表反射太阳辐射能量强弱的物理量。地表反照率是辐射强迫的直接驱动因子之一，陆表覆盖变化也是通过改变地表反照率，进而影响辐射强迫，最终导致区域或全球气候环境发生变化（Dickinson，1995；Trenberth 等，2009），因此，研究地表反照率的时空分布特征及其变化趋势对研究全球或者区域气候变化具有十分重要的意义。

遥感卫星以其空间覆盖广的优势，为区域乃至全球尺度地表反照率的研究提供重要的数据支撑（Qu等，2015）。遥感卫星辅以不同的地表反照率反演算法，构建了长时间序列的地表反照率遥感产品（Liang等，2013）。在低分辨率尺度（公里级或百米级），地表反照率反演算法主要有基于地表二向反射分布函数（BRDF）模型算法（Schaaf 等，2002）、直接反演算法（Liu 等，2013）、静止卫星反演算法（孙越君 等，2018）。基于以上这些算法，生产了全球尺度的MODIS（Schaaf等，2002）、GLASS（Liang等，2013）和Copernicus（Lacaze等，2013）产品等；在高分辨率尺度（十米级），卫星重访周期长，短时间范围内难以累积足够多角度的地表反射率，面临着角度信息量不足难以满足BRDF 建模需求（Gao 等，2014）的问题。目前，高分辨率地表反照率反演算法主要是利用低分辨率的BRDF 产品提供先验知识，如基于MODIS BRDF 产品（MCD43A1）构建的“pre-MODIS era”算法（Shuai 等，2011）、“MODIS-concurrent”算法（Shuai 等，2014）、利用Landsat 的光谱响应函数和MODIS BRDF 数据库构建辐射传输模型的直接反演方法等（He等，2018）。上述算法在均一性较好的地表取得了较好的效果（Lin 等，2018a），但在山区地表下，地形改变了太阳—地形—传感器之间的几何关系，导致山区地表太阳辐射传输过程较平坦地表更为复杂（林兴稳 等，2020）。地形坡度、坡向、光照、阴影等效应，导致同种地表类型在不同坡度和坡向下，地表反射率变化各异（Combal等，2000；Schaaf等，1994）以及地表BRDF 呈现出不对称性（Wu 等，2019；Yin 等，2018）。目前BRDF 模型没有考虑地形的影响，反演的地表反照率产品在山区地表存在一定的误差（Lin 等，2018b）。林兴稳（2018b）将山地辐射传输模型和直接反演算法结合，初步构建了适用于山区地表反照率遥感反演方法，并初步应用于黑河地区的高分辨率环境卫星的地表反照率的反演，但是山区云、雪多，云污染、季节性积雪等导致反演的地表反照率产品的缺失值较多，产品在空间分布上的存在缺失和时间序列上的存在不连续问题，严重制约山区地表反照率应用。

低分辨率遥感反照率产品具有良好的时间连续性，但其空间分辨率过低（Holben，1986；Justice 等，1985）；高空间分辨率卫星空间分辨率高，重访周期长（Coops 等，2006）。许多学者将低空间分辨率遥感数据的高频率观测优势与高空间分辨率遥感数据的空间分辨率高的优势相结合，构建了基于不同理论的时空融合算法（Li 等，2020）。这些算法大体上归纳为5 类；（1）基于权重函数的方法，如STARFM（Spatial and Temporal Adaptive Reflectance Fusion Model）系列模型，这类方法能够准确地捕捉物候变化，但精确性受限于景观的特征斑块大小，且采用简单的线性模型，会影响山区地表参量反演（Gao等，2006）。（2）基于混合像元分解的方法，如MMT模型（Zhukov等，1999）、STDFA 模型（Wu 等，2012）等。该类模型采用混合像元分解模型，利用长时间序列的低空间分辨率影像提取地表反射率的非线性时间变化特征，进而与高空间分辨率影像进行融合（吕长春 等，2003）。但模型假设像元反射率不随时间发生变化，与实际情况存在不符。（3）基于机器学习的方法，如SPSTFM 机器学习模型（Huang 和Song，2012），该方法有较好的稳定性与可信度，但是算法过程中需要Landsat 与其对应的MODIS 影像，导致算法在实际应用中的难度增加（蔡家骏，2019）。（4）基于贝叶斯的方法，如TBDF 模型，该方法将时间相关信息融入到图像时间序列中，将融合问题转化为估计问题，通过最大后验估计得到融合后的图像，较为适合于异质性地表（Xue等，2017）。（5）混合类方法，将基于权重函数、贝叶斯理论或者像元分解等方法相结合，得到更好的结果，如FSDAF 模型混合使用了像元分解及权重函数的方法（Zhu 等，2016）。上述方法大多未考虑山区地形，且受到景观斑块大小的影响，较难适用于山区地表时空连续的高空间分辨率地表反照率反演。卡尔曼滤波算法（Kalman，1960），以分析误差的最小方差为最优标准，通过递归处理，可以保持过程结果，使之始终处于最优状态；Evensen（1994）提出了集合卡尔曼滤波资料同化的基本理论框架，克服了卡尔曼滤波仅限于处理线性问题的弱点，不再需要切线性模型和伴随模式，解决了计算量巨大的问题。集合卡尔曼滤波EnKF（Ensemble Kalman Filter）算法通过集合预报与资料同化，在给定观测和预报模式下去描述参量的概率分布及其发展（官元红 等，2007），不受限于景观的斑块大小，故在山区也有很大的应用空间。

本文首先基于改进的直接反演算法获取山区高分辨率地表反照率，利用集合卡尔曼滤波方法，构建山区长时间序列高空间分辨率反照率产品，解决反演的反照率在山区地表下时空缺失值较多的情况。最后，将反演的反照率产品与地面站点实测数据进行对比，评价构建的时空连续地表反照率产品的精度和不确定性。

2 研究区与数据源

2.1 研究区域

黑河流域是中国西北地区第二大内陆流域，位于河西走廊农牧交错带上，横跨3 种不同的自然环境单元，即上游地区主要为祁连山区、中游地区是平原走廊、下游地区则以低山山地为主，整个流域地区是典型的生态敏感区（赵晓冏，2012）。本研究区位于黑河流域南部，整个研究区高程范围为889—5218 m，太阳辐射强烈，地表覆盖类型包含农田，林地，草地、湿地、裸地、冰雪覆盖区等，主导地表类型为裸地及草地（图1）。该地区因独特的地理位置及气候条件，是寒区和旱区生态、水文陆表过程研究的理想区域。大量科学研究与实验已经积累非常丰富的数据资料，已开展的科学实验包括“非均匀下垫面地表蒸散发的多尺度观测试验”（Liu等，2018）、“无人机物联网中继实现无公网偏远地区环境监测”（Zhang 等，2020）、“黑河流域生态—水文过程综合遥感观测联合试验（HiWater）”等（李新 等，2012）。其中，HiWater 实验是黑河流域的多尺度综合观测试验，构建了长时间序列的水文气象观测网络。该水文气象观测网包含有不同下垫面的地面观测站点，提供多种观测数据集，包含风速、气压、四分量辐射、地表辐射温度、光合有效辐射、土壤热通量、土壤水分、土壤温度等，每10 min 记录一次并对外发布。

图1 黑河流域研究区高程及土地利用类型图Fig. 1 Elevation and land cover maps of Heihe Watershed area

2.2 数据源

2.2.1 地面数据

地面验证数据采用青藏高原科学数据中心发布的2016 年和2017 年“黑河流域生态—水文过程综合遥感观测联合试验”地面辐射观测数据集（Liu 等，2018）。由于MODIS 反照率已订正到太阳正午时刻的值（Wu等，2017）。考虑到黑河地区的正午时刻约为下午13：00时左右，故选取12：30 pm至13：30 pm 之间观测的地表上行短波辐射值与太阳下行短波辐射值（观测波段范围：0.3—3.0 μm），计算比值并求平均值作为每天的反照率观测值，计算式为式（1）：

式中，ar表示日均当地正午时刻地表短波反照率，ri表示每10 min短波上行辐射，表示每10 min短波下行辐射，n＝7。为验证算法在不同地表的适用情况，选取了4 个不同下垫面的站点进行算法检验；考虑到MODIS反照率产品及GF-4数据像元大小，统计站点周围2 km 范围内的坡度分布范围，站点信息见表1。

表1 辐射通量观测站点信息Table 1 information of radiation flux observation site

2.2.2 卫星数据

GF-4 卫星是中国第一颗地球同步轨道遥感卫星，可见光和多光谱空间分辨率优于50 m，红外谱段空间分辨率优于400 m，实现对中国及周边地区的观测（孙越君 等，2018），GF-4 传感器各波段参数见表2，本文所用GF-4 反射率数据见表3。经过辐射定标等预处理的GF-4 反射率数据用于坡面反照率反演。

表2 GF-4卫星传感器参数Table 2 Observation parameters of GF-4 satellite sensor

表3 GF-4影像数据Table 3 Information of GF-4 data

MCD43A3 C6 产品是日尺度16 d 合成的500 m地表反照率产品。是通过筛选近16 d质量较好的像元，组成多角度观测数据集，进而拟合BRDF 得到地表反照率（Schaaf等，2002）；该产品经过广泛测试使用，精度已得到普遍认可（冯智明 等，2018）。

地表高程数据采用国家青藏高原科学数据中心发布的ASTER 卫星90 m 地表高程数据（https：//data.tpdc.ac.cn/zh-hans/data［2020-08-03］），地表覆盖类型采用清华大学发布的2017年30 m 地面分类数据（http：//data.ess.tsinghua.edu.cn/［2020-08-03］），地表覆盖类型共分为10 类，包括农田、森林、草地、灌木地、湿地、水、苔原、不透水面、裸地和冰雪（图1）。

3 方法构建

3.1 算法流程

算法流程图如图2所示，首先利用山区地表反照率反演方法得到GF-4 坡面地表反射率（Lin 等，2018b）；然后通过研究区坡度、坡向、地表高程数据、土地利用数据对MCD43A3 数据进行逐像元处理，利用高空间分辨率地表覆盖数据提取纯像元，并将MCD43A3 纯像元的观测值作为时空连续地表反照率反演的背景场；在此基础上代入地面站点观测数据和GF-4 坡面地表反照率，通过EnKF 算法得到长时间序列高分辨率反照率数据；最后进行反演结果的精度验证。

图2 长时间序列山区地表反照率遥感反演流程图Fig. 2 The flow chart of long time series land surface albedo retrieval over rugged terrain

3.2 山区地表反照率的直接反演算法

Lin 等（2018b）基于单一角度反照率反演算法，初步构建了山区地表反照率反演算法。该方法假设高分辨率像元内地表具有单一的坡度和坡向，考虑树木向地性生长；在斜坡上构建坡面坐标系，将太阳入射和卫星观测几何格网化，构建相对太阳入射天顶角、相对太阳入射方位角、相对卫星观测天顶角和相对卫星观测方位角的格网组合。在坡面坐标系的每个格网中，利用核驱动模型，构建不同地表类型下的地表反射率和宽波段地表反照率之间的查找表（Liu 等，2013；孙长奎 等，2013），进而在坡面坐标系中反演坡面地表反照率。

式中，ρdd（is，φs，iv，φv）、ρhd（iv，φv）、ρdh（is，φs）、ρhh分别表示为坡面像元方向—方向反射率、半球—方向反射率、方向—半球反射率和半球—半球反射率。根据定义可将ρhd（iv，φv）、ρdh（is，φs）、ρhh展开如下：

将式（3）—（5）分别带入到式（2）中，构建关于ρdd（is，φs，iv，φv）的一元二次方程，利用各地表类型的MODIS纯像元提供BRDF先验知识，在坡面坐标系下，对坡面方向—方向反射率进行积分获取坡面半球—方向、方向—半球和半球—半球反射率，进而求解关于ρdd（is，φs，iv，φv）的一元二次方程，获取坡面地表反射率。利用构建的坡面地表反射率和宽波段地表反照率的查找表，直接反演得到地表反照率（式（6））

式中，αABT表示高分辨率坡面反照率，is，φs，iv，φv分别表示坡面坐标系上相对的太阳天顶角、太阳方位角、观测天顶角和观测方位角；ρi（is，φs，iv，φv）表示第i波段的地表反射率；Ci（is，φs，iv，φv）表示坡面坐标系下的回归系数。

3.3 山区高分辨率地表反照率时空填补方法

本文采用改进的集合卡尔曼滤波方法对反演的山区地表反照率遥感产品进行处理，以MODIS地表反照率产品和地面实测数据为先验知识，构建时空连续的山区地表反照率产品。集合卡尔曼滤波方法常用来进行资料同化，充分有效的利用已有的各种信息预测未知的状态，关键步骤是计算过程中的状态更新（朱琳，2007）；应用在山区高分辨率地表反照率反演中，有两点需要考虑，一是背景场的构建，二是集合卡尔曼滤波方法构建。

3.3.1 先验知识背景场的构建

MODIS 反照率数据未考虑地形影响，将其作为背景场时，须先对其进行像元尺度的处理，增强其在山区的适用性。本文首先以GF-4 影像空间分辨率和坐标系统为基准，对MODIS 反照率产品进行重采样与投影变换；然后，在GF-4 影像上建立搜索窗口，搜索与窗口中心像元相似的像元，确定每个像元周围的相似像元位置，主要从两方面进行相似像元的选取：一是与窗口中心像元具有相同地物覆盖类型、相似高程及坡度坡向的像元视为相似像元；二是计算中心像元与邻近像元的反照率差异值，通过阈值判断差异是否合理，在合理区间内的像元被视为相似像元。阈值可以通过高分辨率图像中像元总体的标准差和估计的土地覆盖类别数来确定（式（7））（Gao等，2006）；

式中，Fi、Ci分别表示最邻近时刻第i个相似像元位置对应的高、低空间分辨率的反照率值，cov（.）表示求协方差，D（Fi）、D（Ci）分别表示Fi、Ci的方差，R的范围［-1，1］，R值越大表示越相似。第i个相似像素与中心像素的地理距离di可由式（9）计算，窗口内相似像素的距离范围［1，1+20.5］（Zhu等，2010）。

结合光谱相似性和距离，可以计算出光谱和地理距离结合的综合指数D（式（10））：

式中，D越大，表示相似像元的贡献越小，利用倒数归一化形式（式（11）），得到权重W表示为

式中，Wi的范围为0 到1，所有相似像素的总权重为1；故可计算tp时刻搜索窗口中心像元的高分辨率反照率（式（12））：

式中，N为搜索窗口内通过式（7）得到的近似像元个数，包含中心预测像元；C为相似像元位置对应的重采样后MODIS 产品反照率值；vi为第i个相似像元的转换系数，是通过线性回归分析窗口内GF-4 反照率数据和MODIS 反照率数据的相似像元，建立回归模型，取其斜率得到。

3.3.2 集合卡尔曼滤波算法

集合卡尔曼滤波算法根据背景场以及观测值的特征误差分布，对背景场和观测值进行分析，得到一组分析值，再用这组分析值的差异作为分析误差的统计样本进行分析误差协方差的估计；对这组分析值作一个短期预报后，得到一组预报值。把这组预报值的差异作为背景场误差的统计样本进行误差协方差的估计，以此递归处理整个系统（Li 等，2012）。当测量及观测误差已知时，集合卡尔曼滤波器可以从这些数据中估计系统的动态变化，其动态模型一般表示为

式中，X表示动态模型状态向量，本文中为高空间分辨率的地表反照率，Xk‒1，k表示预测的高空间分辨率反照率；Mk‒1，k是k‒1到k时刻的地表反照率的状态转移矩阵，考虑前后两日正午时刻地表反照率之间互不影响，且变化较小，此处将M设为1。Bk‒1，k是k‒1 到k时刻的系统控制矩阵，Uk表示k时刻的状态控制量，本研究是直接获取地表反照率，没有控制量，所以Uk=0，Bk不需要考虑。vk-1为k‒1 时刻的观测误差，符合正态分布。通过式（14）对k时刻的高分辨率地表反照率Xk进行预测：

式中，Kg为卡尔曼增益矩阵；Zk表示集合卡尔曼滤波算法中的实际测量值，在本研究中即为地面站点实际测量反照率值；Hk是测量系统的参数，将状态向量转化为观测类型的值，在本研究中，观测值与状态向量均为地表反照率，即已直接获得反照率的测量结果，故H设为单位矩阵。

式中，Zk是k时刻经过处理的地面站点实际测量反照率值。Kgk可以通过式（16）表示：

式中，R是测量噪声协方差矩阵，HT表示H的转置矩阵。由上一预测时刻的误差协方差Pk‒1和过程噪声矩阵Q预测新的误差协方差矩阵Pk-1，k，可表示如下：

式中，Pk‒1，k是Xk对应的协方差矩阵，Q是系统过程的噪声协方差，MT表示M的转置矩阵。此时得到了k状态下最优的估算值X。为使卡尔曼滤波器不断的运行下去直到系统过程结束，还需更新k状态下X的协方差Pk：

式中，I为单位矩阵。

在本研究中，以3.1 节得到的结果为背景场。考虑到研究区域内并不是每个像素都具有地面测量值，故计算区域内测量误差的协方差矩阵R时，基于站点误差构建了均值为0，方差为0.01 的正态分布矩阵，作为区域内像素的测量误差；即根据站点像素的测量误差设置区域的测量误差，进而计算区域内测量误差的协方差矩阵R；同理，在预测新的误差协方差矩阵P的过程中（公式（17）），将系统过程噪声设置为均值为0，方差为0.01 的矩阵，求得系统过程噪声协方差矩阵Q，进而得到新的误差协方差矩阵P（Zhang等，2019）。

4 结果与验证

4.1 山区时空连续地表反照率反演结果

本文算法利用日尺度MODIS 反照率产品构建先验知识背景场，利用地面观测数据和集合卡尔曼滤波算法，结合数幅高质量的GF-4 影像，获得日尺度50 m 的反照率产品。结果显示，本文的方法可以构建时空连续的高分辨率地表反照率产品，且随着高质量GF-4 影像的增加，同化效果得到提升。为体现算法的普适性，选取四个不同地表覆盖类型的站点进行结果验证。根据各站点首次获得最优高分辨率反照率产品为起始日期，下图展示了各站点地表反照率随时间的变化趋势（图3）；由于大满及花寨子站点距离较近，在图3（c）中用一组图展示。

图3 2017年研究区内不同类型地表反照率反演结果图Fig. 3 Shortwave albedo maps generated by EnKF for different land surface types of study area in 2017

图3（a）为阿柔研究区2017 年第100—300 D地表反照率变化图，其下垫面为亚高山山地草甸，地表反照率变化最为明显，尤其是中下部分。中部颜色由褐色逐渐变为黄色，下半部分深绿色逐渐增加，表明随着时间变化，该区域地表反照率在下降；中部下垫面以草甸为主，此外还存在部分农田，随着草甸回绿及农作物的生长，地表反照率不断下降至0.2 左右；南部山区人迹罕至，下垫面包含森林、草地和裸地，图3（a）中青色部分地表反照率较小，该区域为森林覆盖，反照率较周围草地、裸地等要小。图中红色部分地表高程高，下垫面以草地和裸地为主，部分地区山顶常年积雪，反照率值逼近0.6，但随着季节变化，大部分山顶冰雪消融，地表反照率值下降；图3（b）为张掖站点，深绿色的上半部分为芦苇湿地，海拔较低，在该时间段生态情况较为稳定，变化浮动较小，稳定在0.17 左右。下半部分下垫面为裸地，常年地表反照率稳定保持在0.24左右；图3（c）中部为大满气象站点所在地，下垫面为玉米地，该区域在2017 年前100—240 d 左右地表反照率一直处于下降状态，从0.23 左右降到0.15 附近，在图中体现为中间部分绿色一直不断加深，对应玉米的灌浆期、拔节、封垄期地表反照率的变化；接着反照率出现上升，应是因为在9月底左右农作物处于收割阶段，大量农作物的收割导致该现象的出现；图3（c）下方为花寨子站点，黄色部分的下垫面为山前荒漠，该区域在时间尺度上变化不大，反照率稳定在0.23 左右；左下角绿色区域不断加深，为植被区（姚彤和张强，2014）。图4 展示了阿柔站点2017年第300 d地表高程与反演反照率的局部放大对比图，从中可以直观的发现山区坡度坡向变化引起的反照率变化；在图4中，高程较低的区域下垫面以草甸为主，此外还存在部分农田，为地表反照率值较大，随着高程的增加，阳坡地表类型以草地位置，地表反照率较大，阴坡地表以针叶林为主，地表反照率较阳坡地表低。在地表高程更高的地方存在积雪区域，反照率较高，部分地区的地表反照率值要大于0.4。

图4 地表反照率反演结果随地表高程变化Fig. 4 The distribution of shortwave albedo retrieved by EnKF following the change of elevation

图5 为阿柔、张掖、大满、花寨子研究区2017 年时间段内反照率及与地面实测地表反照率随时间变化趋势图，从中可知，除了张掖湿地站点外，在其余3个站点下，MODIS产品较地面实测反照率存在低估现象，但在时间序列上，MCD43A3产品的变化趋势与地面站点相同，能较好的反映地表的季节性变化特征。张掖站点的MCD43A3 产品较地面观测值偏高，该站点的下垫面为湿地，地面站点观测值可能包含来自水体的反射信号。MCD43A3 产品的空间分辨率为500 m，在该站点下两者之间的尺度效应明显。花寨子站点的地面观测数据波动性较大，通过地面实测数据的时间序列分析发现，花寨子站点前后相邻日期的地表反照率变化较大，可能是站点观测仪器稳定性、观测条件等因素影响，导致地面反照率降低，进而影响了实际测量，但整体局势符合实际情况。在这4个站点上，GF-4影像较少，利用GF-4卫星直接反演的地表反照率与地面站点的一致性较好。利用本文算法得到的时空连续的地表反照率与地面站点的一致性也比较好，比MODIS 产品更能反映地表反照率的时空变化，准确性也更高，同时也可以很好的保持时间上的连续性。比较本文算法反演的2017 年4 个站点的时空连续的地表反照率与站点实测地表反照率可知，4个站点的R2分别为0.8701、0.9865、0.9882、0.4466，大多高达0.8；RMSE分别为0.0116、0.0036、0.0036、0.0123，皆小于0.02。

图5 2017年4个研究区反演结果与地面实测对比图Fig. 5 Comparison between the inversion results and four observation sites measured data

4.2 预测地表反照率和直接反演地表反照率比较结果

图6—图8展示的是本文算法预测的GF-4地表反照率与同一天的利用当天已有反射率数据采用直接反演算法得到的GF-4 地表反照率的比较结果。张掖、大满及花寨子站点对比日期选取了2017 年6 月12 日，阿柔站点选择了2016 年12 月15 日。从各站点反演结果与反照率产品对比图可以看出，本文预测的地表反照率与直接反演的GF-4 地表反照率的空间分布具有很高的相似性。与直接反演算法得到的GF-4 反照率数据相比，模拟结果在不同的坡度和坡向下特征明显，表现出地形特征效果。在不同下垫面下，尤其是阿柔、张掖等站点存在较大坡度，但站点反演结果与GF-4 反照率数据有较好的可比性，整体反演情况较好；说明本方法在模拟过程中能有效地利用已有数据，对崎岖地表完成较好的模拟。反演结果与反照率产品对比图显示张掖站点的计算结果与同时期GF-4 反照率数据值集中在0.17—0.32，统计结果显示两者之间的RMSE 为0.0045，R2为0.9623；大满及花寨子站点的计算结果与GF-4反照率的对比，RMSE为0.0096，R2为0.7372；阿柔站点的预测的地表反照率与GF-4直接反演的地表反照率的对比，RMSE为0.0163，R2为0.9163。RMSE 皆小于0.02，R2大多高达0.8，这说明了本文算法预测的地表反照率精度较好，可以使用于长时间序列的地表反照率的构建。

图6 张掖站点2017年6月12日GF-4反照率产品与反演结果对比Fig. 6 Comparison between the GF-4 albedo and inversion results of Zhangye site on June 12， 2017

图7 大满及花寨子站点2017年6月12日GF-4反照率产品与反演结果对比Fig. 7 Comparison between the GF-4 albedo and inversion results of DaMan and HuaZhaizi site on December 12， 2017

图8 阿柔站点2016年12月15日GF-4反照率产品与反演结果对比Fig. 8 Comparison between the GF-4 albedo and inversion results of ARou site on December 25， 2016

5 结 论

本文参考前人工作提出一种山区地表时间序列反照率的估算方法，以地面站点观测数据和时空连续的MODIS反照率产品构建先验知识背景场，结合少量的GF-4 数据通过集合卡尔曼率波算法得到黑河流域长时间序列高时空分辨率的反照率产品。通过与地面站点地表反照率测量数据对比，在不同地表坡度和地表类型下，反演的地表反照率与实测值时空一致性较好，RMSE 皆小于0.02，R2高达0.8。与直接反演的GF-4 反照率数据相比，预测的地表反照率与同期的地表反照率产品之间的一致性较好，全景影像的统计结果显示两幅影像的RMSE 小于0.02。考虑到GF-4 卫星高空间分辨率的特点，通过增加GF-4 影像数量，能够得到更精确的预测结果。本文的研究具有以下创新点：第一，考虑了山区地表的地形效应，在通过邻近像元预测反照率时，通过选取相似像元，考虑了坡度坡向、高程、地表分类等信息，增强了背景场在山区的适用性，同时发挥其时间连续性优势，以此作为背景场来反映地表随时间变化的物候规律；第二：利用集合卡尔曼滤波算法，通过使用较为少数的高精度影像，结合MODIS反照率产品，利用算法生成计算出全年的高时空分辨率反照率产品，有效的提高了山区反照率反演精度。

本文能够为其他卫星数据或者生物量的反演提供一种思路，然而具体在其他卫星上的应用效果还需进一步检验；考虑到已发布的地面实测数据及卫星数据，仅使用较为少数的站点进行算法验证，后续应增加站点数量及观测年份，提高算法普适性。