基于Landsat和MODIS数据融合的农牧区NPP模拟

尹小君 祝宏辉 GAO Jerry 高 军 郭丽洁 苟贞珍

(1.石河子大学信息科学与技术学院， 石河子 832003； 2.兵团空间信息工程技术研究中心， 石河子 832003；3.石河子大学经济管理学院， 石河子 832003； 4.圣何塞州立大学计算机工程系， 圣何塞 CA 95192；5.安徽四创电子股份有限公司， 合肥 230094)

0 引言

植被NPP是指绿色植物经过光合作用在单位时间、单位面积所获取的有机物质总量与维持自养呼吸所消耗量的差值[1-3]，它是植被自身生理特性与外界环境共同作用的结果[4-5]，能够直观反映地表植被在外界环境作用下的生产能力[6]，是评判陆地生态系统健康状况以及研究碳源/汇过程的重要指标[7-9]。天山北坡是我国西北干旱区的主体，同时也是我国重要的农牧业发展基地[10]，由于近年来气候变化、人类活动以及土地利用变更等因素的影响，已造成天山北坡区域内草地退化、林地萎缩、地下水位下降等一系列生态环境问题[11-13]。因此，适时进行区域内植被NPP的动态模拟，掌握植被生长健康状况，对研究植被动态变化、合理分配农牧场草地资源至关重要。

随着遥感技术的不断发展，遥感数据以其成像周期短、覆盖范围广、成本低廉等特点，可快速、准确地获取植被NPP的动态变化信息，为植被生长健康状况监测提供了有效途径。目前，基于遥感的植被NPP模拟模型主要分为3类[14]：气候生产力模型、光能利用率模型和植被生理生态过程模型。其中，基于光能利用率的CASA(Carnegie-Ames-Stanford Approach)模型可操作性强，在不同地域、不同研究尺度的NPP模拟中得到了广泛应用[15-17]。自2000年以来，MODIS数据因其时间分辨率高，能够较好地反映区域内NPP时间序列变化信息，而成为研究植被NPP动态变化的主要数据源[18-20]。但MODIS数据受到空间分辨率的限制，导致植被变化的空间细节信息表达不够精细，难以满足在精细尺度下的NPP模拟。为此，学者们开始寻求中分辨率遥感数据对区域植被NPP进行模拟。池源等[21]基于Landsat 8 OLI影像探讨了庙岛群岛北五岛景观格局与植被NPP的关系；管小彬等[22]基于TM/ETM+影像分析了武汉市冬季植被NPP时空变化特征；杨会巾等[23]基于Landsat 8 OLI对玛纳斯河流域植被NPP进行了模拟。然而，由于中分辨率卫星传感器重访周期较长，且易受到云、雨、雪等自然天气的影响，使得获取连续时间序列的NPP数据受到了限制。

为满足大范围、高精度及快速变化的地表信息遥感监测需求，多源遥感数据时空融合技术为解决上述问题提供了有效的途径。目前，国内外学者已将多源遥感数据时空融合技术广泛应用于不同领域[24-26]，因此，利用该方法将不同时间分辨率、空间分辨率和光谱分辨率的遥感数据进行时空融合，获取同时具有高时间分辨率遥感数据的时相信息和中空间分辨率数据的空间纹理信息，可为陆地植被NPP的时空动态变化研究提供有效的数据源。本文以天山北坡中段区域为研究区，探讨利用多源遥感数据时空融合技术获取中空间分辨率遥感影像，模拟研究区2016年的植被NPP，并对融合数据的有效性进行精度检验，以期为多源遥感影像融合技术与光能利用率模型协同模拟NPP提供新的思路。

1 研究区概况及数据获取

1.1 研究区概况

研究区位于新疆天山北坡中段绿洲核心区，地理位置为43°40′～45°33′N，84°72′～86°56′E，海拔260～4 900 m，总面积为2.38×104km2，是新疆典型的农牧业发展基地，如图1所示。研究区地处亚欧大陆腹地，远离海洋，属于典型的温带大陆性气候，区域内四季分明，夏季炎热干燥，温差较大，蒸发强烈，冬季寒冷漫长，降水量稀少，生态系统比较脆弱。

图1 研究区行政区划图Fig.1 Administrative map of study area

研究区植被类型分布主要与海拔有关，因不同的海拔梯度，植被整体呈现垂直结构分布。在研究区低海拔山麓平原区，主要分布荒漠耐旱植被。海拔2 000 m左右主要分布低山草原，是研究区主要的植被类型，也是家庭牧场的主要聚集区，其中森林与草甸相间出现。高海拔区域，主要分布一些高山草甸和苔藓地衣。

1.2 数据及处理

Landsat8是NASA于2013年成功发射的卫星，有2个主要载荷：OLI和TIRS。其中OLI为陆地成像仪，TIRS为热红外传感器。选取2016年研究区内数据质量较好，云层覆盖度小于5%的Landsat 8 OLI影像，空间分辨率为30 m，轨道重访周期16 d。影像行列号为144/29和144/30，数据来源于中国科学院计算机网络信息中心的“地理空间数据云平台”(http:∥www.gsclud.cn)。利用Landsat 8 OLI 近红外波段反射率和红光波段反射率，获取NDVI数据[27]。

MODIS是美国进行对地观测系统战略实施发射的重要光谱传感器，该数据不仅可以免费获取，而且数据时间序列范围广，更新频率快，能够较好地获取快速变化的自然植被等信息。MODIS搭载在Terra和Aqua 2颗卫星上，每天获取最少2次白天和2次黑夜MODIS数据。MODIS陆地标准数据产品具有每日产品、8 d合成产品、16 d合成产品、月合成产品、季度合成产品和年合成产品。本文使用的MODIS数据为MODIS13Q1产品中的NDVI数据，空间分辨率为250 m，时间分辨率为16 d。数据来源于美国陆地过程分布式数据档案中心(https:∥lpdaac.usgs.gov/)。影像获取时间为2016年1月1日—12月18日，共23期，行列号为h23v04、h24v04。数据预处理首先使用MRT(MODIS reprojection tool)软件对影像进行镶嵌、格式转换和重投影等操作，再利用ArcGIS软件对研究区范围进行裁剪。气象数据选取研究区内及周边气象站点的日平均温度、日降水总量和日太阳总辐射数据资料，该数据来源于中国气象数据网(http:∥data.cma.cn/)，将所得日数据信息合成月数据信息，利用澳大利亚科学家Hutchinson基于薄盘样条理论编写的专业气象数据插值工具ANUSPLIN[28]，以经纬度为自变量，高程作为协变量，对研究区内气象数据进行空间插值，获取气象栅格数据。

研究区地形数据采用SRTM 30 m DEM产品，数据来源于“地理空间数据云平台”(https:∥www.gscloud.cn)。土地利用覆被数据来源于国家地球系统科学数据共享网服务平台(http:∥www.geodata.cn/)。所有栅格数据的行列数、像元大小及投影方式均与Landsat 8 OLI数据一致。

于2016年7月18日、7月26日、8月9日、8月18日和8月26日，在天山北坡中段农牧区，即在牧场草地生长期对草地地上生物量进行采集，根据DEM数据掌握当地的地形，分别在阳坡和阴坡设置采样点140个，如图2所示。根据实际情况选取一定数量的样方，样方选取为1 m×1 m，在实验区每隔300 m重复的点采样3次，其平均值作为样本值，共计46个样本。齐地剪取样方地上生物量，剪取样本的同时用手持GPS定位仪，记下单位样方的经纬度，样本装袋并编号统一带回实验室称量。样本数据的详细信息主要包括生物量干质量、湿质量、样点坐标、高程、覆盖度、经纬度、海拔、时间、草地群落类型等。

图2 采样点空间分布图Fig.2 Spatial distribution map of sampling points

2 研究方法

基于Landsat数据与MODIS数据，采用时空融合技术实现对研究区2016年NPP的模拟。首先将研究区遥感数据进行预处理，获取研究区内不同数据源的NDVI数据，然后利用STARFM数据融合算法得到研究区中空间分辨率时间序列的NDVI(30 m)数据，其次结合已处理好的气象数据(温度、降水、太阳总辐射)、土地利用覆被数据采用CASA模型对研究区2016年植被NPP进行模拟，最后根据研究区实测数据对估测NPP结果进行精度验证，其具体流程如图3所示。

2.1 遥感数据时空融合

采用的遥感数据时空融合算法为时空适应性反射率融合算法模型(Spatial and temporal adaptive reflectance fusion model，STARFM)[25]，该方法是基于t0时刻的中空间分辨率影像与低空间分辨率MODIS影像，结合t1时刻低空间分辨率MODIS影像，通过计算像元间的空间分布差异及波谱差异，模拟预测t1时刻中空间分辨率数据的方法。基于t0时刻获取的中空间分辨率Landsat 8 OLI NDVI数据(2016年6月24日)和低空间分辨率MODIS NDVI数据(2016年6月25日)以及t1时刻获取的低空间分辨率MODIS NDVI数据(2016年4月6日和2016年8月26日)，预测t1时刻Landsat 8 OLI NDVI 数据(2016年4月6日和2016年8月26日)，预测过程通过使用滑动窗口法减少低分辨率遥感数据像元边界的影响[29]。其计算公式为

(1)

其中

(2)

式中L(xw/2,yw/2,t1)——预测t1时刻的Landsat 8 OLI NDVI像元值

w——移动窗口尺寸

xw/2、yw/2——窗口中心的像元坐标

M(xi,yi,t1)——窗口位置(xi,yi)处在t1时刻的MODIS像元NDVI值

L(xi,yi,t0)——t0时刻给定位置(xi,yi)中Landsat像元的NDVI值

M(xi,yi,t0)——t0时刻给定位置(xi,yi)中MODIS像元的NDVI值

Wijk——移动窗口内各像元对预测中心像元值权重

cijk——由移动窗口内的中心像元与窗口内其他像元的光谱距离、时间距离和空间距离共同确定的系数

2.2 CASA模型

CASA(Carnegie-Ames-Stanford Approach)模型是由遥感、气象、土地覆被数据共同驱动的光能利用率模型，是目前NPP模拟模型中应用最广且较为成熟的模型[30]。采用朱文泉等[31]改进的CASA模型对天山北坡中段植被NPP进行模拟，计算公式为

NPP(x,t)=APAR(x,t)E(x,t)

(3)

式中NPP(x,t)——像元x在t月的NPP,g/(m2·month)

APAR(x,t)——像元x在t月植被所吸收的光合有效辐射，MJ/(m2·month)

E(x,t)——像元x在t月的光能利用率，g/MJ

植被吸收的光合有效辐射(APAR)取决于太阳总辐射量(SOL)和植被对APAR的吸收比例(FPAR)，计算公式为

APAR(x,t)=0.5SOL(x,t)FPAR(x,t)

(4)

式中SOL(x,t)——像元x在t月的太阳总辐射量，MJ/m2

FPAR(x,t)——植被冠层对光合有效辐射的吸收比例

在一定范围内FPAR与归一化植被指数(NDVI)、比值植被指数(SR)都存在良好的线性关系[32]，可以通过遥感影像提取植被的NDVI来对FPAR进行模拟，计算公式为

(5)

FPAR(x,t)SR=

(6)

(7)

式中NDVIi,max——对应第i种植被类型的NDVI最大值

NDVIi,min——对应第i种植被类型的NDVI最小值

SRmax——对应第i种植被类型的SR最大值

SRmin——对应第i种植被类型的SR最小值

FPARmax、FPARmin——常数，分别为0.950、0.001

FIELD等[33]研究发现通过这2种方法模拟所得的结果各有不同，由NDVI模拟的FPAR结果比实测值偏高，而由SR模拟所得的FPAR结果比实测值偏低，因此，可以利用调整系数0.5，将两者取平均值的方式作为对FPAR的模拟结果，从而使其得到的FPAR模拟误差最小，计算公式为

(8)

天山北坡中段为洪积冲积扇，呈荒漠-绿洲-山地景观格局，植被覆盖稀疏，把调节系统调整为0.52。

光能利用率是指绿色植被把所吸收的光合有效辐射转换为有机碳的效率，它主要与温度及水分胁迫系数有关，其计算公式为

E(x,t)=TE1(x,t)TE2(x,t)WE(x,t)εmax

(9)

式中TE1(x,t)——低温作用下对植被光能转换率的胁迫系数

TE2(x,t)——高温作用下对植被光能转换率的胁迫系数

WE(x,t)——水分条件对植被光能转换率的影响系数

εmax——最大光能转换率

3 结果与分析

3.1 数据融合精度

为验证Landsat与MODIS数据采用STARFM模型融合生成NDVI数据的精度，本文以2016年6月24日的Landsat 8 OLI NDVI数据和2016年6月25日的MODIS NDVI数据为基础，分别与2016年4月6日和2016年8月12日的MODIS NDVI数据进行融合，得到融合后的2016年4月6日与2016年8月12日的NDVI数据，空间分辨率为30 m。然后，将融合后的两期NDVI数据分别与真实的Landsat 8 OLI NDVI数据进行对比分析。如图4所示，通过分析可知基于STARFM融合算法融合的NDVI数据相比于原MODIS NDVI数据，空间分辨率得到显著提高，且具有良好的空间细节信息，能够较好地表达较小地物之间空间差异，与真实的Landsat 8 OLI NDVI数据相比，具有良好的空间分辨率，且空间分布趋势基本一致。

图4 MODIS、Landsat 8 OLI与融合的NDVI对比分析Fig.4 Comparative analysis of NDVI under different data sources

为验证融合数据与真实数据之间的差异，除了通过目视判别融合效果外，还可通过相关系数r、偏差Bias、均方根误差RMSE等指标对影像融合精度进行定量评价[34-35]，其偏差计算公式为

(10)

式中N——像元个数i——像元序号

xi——融合数据像元值

yi——真实数据像元值

利用2016年多期融合的空间分辨率为30 m的NDVI数据，与对应时刻的空间分辨率为30 m的Landsat 8 OLI NDVI真实影像进行对比(表1)，检验结果表明，融合生成的NDVI数据与真实NDVI数据之间的r不小于0.759，Bias在0.006 2～0.009 4之间，RMSE在0.074～0.135之间，因此，说明利用该融合算法生成的NDVI数据精度较高，与真实影像具有较好的一致性，能够用于后续植被NPP的模拟。

3.2 NPP模拟结果及精度验证

3.2.1NPP模拟结果

利用CASA模型模拟研究区2016年植被NPP空间分布，结果如图5所示，基于MODIS NDVI数据模拟的植被NPP平均值为327.28 g/(m2·a)，基于融合的NDVI数据模拟的NPP平均值为305.47 g/(m2·a)，两者空间上整体呈现南高北低的趋势，沿天山山脉一侧植被的NPP平均值最高，原因主要为天山上的融雪积水为植被的生长提供了必备的有机质，相对较高的海拔使得人类活动对植被的干扰强度减弱，因此NPP均值较高。基于不同数据源虽在模拟研究区内植被NPP空间分布状况整体相差不大，但是由于空间分辨率的较大差异，导致在反映植被内部变化信息时差异较为明显，整体上，基于MODIS NDVI数据模拟的NPP结果由于空间分辨率较低，导致研究区内存在大量的混合像元，难以清晰地表达不同地物间的空间纹理细节信息，而基于融合的NDVI数据模拟的NPP结果空间分辨率相对较高，影像内存在混合像元的数量大大减少，空间纹理信息较清晰，能够很好地表达微小地物间的空间差异。基于MODIS NDVI数据模拟的NPP像元大部分集中在240～660 g/(m2·a)之间，变化幅度不明显，而基于融合数据模拟的NPP像元大部分集中在320 g/(m2·a)左右，整体近似于正态分布，能够清晰地表达农田、草地、建筑物、道路、河流等不同地物间的NPP差异。

表1 融合数据与真实数据间的相关系数、偏差及 均方根误差Tab.1 Correlation, bias and RMSE between fused data and true data

图5 研究区2016年NPP模拟结果Fig.5 Simulation of NPP in 2016

3.2.2融合数据模拟NPP精度验证

利用2016年在研究区内采集的46个样本生物量数据，验证基于融合后的NDVI数据模拟研究区植被NPP的精度，由于实测NPP难度较大，因此通常采用实测生物量换算为NPP数据，根据生物量与NPP之间的转换系数0.475，进而得到采样点实测NPP数据[31]。由于采样时间为7、8月，年度植被最为旺盛的时期，其NPP可以代表年度NPP[36]，将实测数据的样本生物量位置信息与CASA模型模拟的NPP位置信息进行对比，如图6所示。 结果表明，利用融合数据模拟的NPP模拟值与NPP实测值基本吻合(R2=0.860 1，P<0.01)。因此利用Landsat与MODIS数据的融合结果能够较好地对研究区植被NPP进行模拟。

图6 NPP实测值与NPP模拟值对比分析Fig.6 Comparison analysis between simulated NPP and observed NPP

4 结论

(1)利用STARFM时空融合算法，将MODIS NDVI数据与Landsat 8 OLI NDVI数据时空信息进行融合，获得中空间分辨率NDVI时间序列数据，既能在时间上较好地保留MODIS数据时间变化信息，又能在空间上反映Landsat数据的空间细节信息，与8个时期真实的NDVI数据相比，相关系数r不小于0.759，偏差在0.006 2～0.009 4之间，均方根误差RMSE在0.074～0.135之间，能够较好地模拟植被NPP。

(2)运用CASA模型模拟NPP，MODIS数据模拟的NPP像元大部分集中在240～660 g/(m2·a)之间，变化幅度不明显，而基于融合数据模拟的NPP像元大部分集中在320 g/(m2·a)左右，整体近似于正态分布，能够清晰表达农田、草地、建筑物、道路、河流等不同地物间的NPP差异。

(3)利用野外实测NPP值进行遥感数据融合与CASA模型协同精度验证，研究区融合后模型植被NPP与实测NPP数据具有较好的相关性，R2=0.860 1(P<0.01)。