NDVI相似性分区下天山地区草地总产草量遥感估算

刘 艳，聂 磊，杨 耘

(1.中国气象局乌鲁木齐沙漠气象研究所,新疆 乌鲁木齐 830002; 2.中亚大气科学研究中心,新疆 乌鲁木齐 830002; 3.武汉大学测绘遥感信息工程国家重点实验室,湖北 武汉 430079; 4.长安大学地质工程与测绘学院,陕西 西安 710054)

新疆现有牧草地总面积5 116.07万hm2，占全国牧草地总面积的19.52%[1]。北疆地区天然草地有效面积约2 377.52万hm2，占全疆46.47%，是新疆草原畜牧业的重要基地[2]。建立草原产草量监测模型以准确、及时地估算草原产草量对确定载畜量、合理安排畜牧业生产和维护草原生态平衡等具有重要意义。新疆牧草产量估算研究主要为90年代初开展的NOAA/AVHRR卫星数据和天然草地牧草产量关系模型研究[2-7]；地面测产数据和同步MODIS遥感数据集合的天山北坡乌鲁木齐草甸、草原和荒漠草原草地牧草估产模型[8]；遥感分类结果和产量估测模型集合的阜康市1990-2008年近20年草地总面积和总产量趋势变化[9]；1991-1996年新疆天山北坡阜康市内不同草地类型实测草地可食产量、环境与遥感资料等综合的多重相关分析和遥感估产技术研究[10]。这些研究主要集中在1990-2000年；研究区域有限且集中分布在天山北坡的乌鲁木齐(南郊)、阜康和阿勒泰；模型采样点少且时空分布不均，致使模型代表性不足；模型主要针对一个观测年进行，时间序列性不强。仅仅使用一个观测年内高寒草甸、山地草甸、山地草原等7种不同牧草类型逐月产量观测数据，进行产量变化特征分析和实现单一牧草类型产量计算。新疆草地分布错综复杂，上述模型代表性和时效性都非常有限。NOAA/NDVI植被指数产品在空间和光谱分辨率上劣于MODIS/NDVI植被指数产品，相对获取不便，特别是近年来的植被指数产品[11]。MODIS/NDVI植被指数产品空间分辨率有250 m、500 m和1 km，共3种，有单日、16日和月合成3类数据，产品时空分辨率能很好地反映草地植被时空变化和估算草地生物量[12-18]。因此，根据前人已有研究成果不足，考虑到如何有效合理地利用有限样点数据，本研究选取MOD13Q1归一化植被指数(NDVI)16 d合成产品为试验数据，从遥感监测植被指数分布特征入手，基于计算机图像直方图识别手段定量评价遥感监测植被指数的区域相似性，以县(市)为研究单元，计算多年植被指数均值直方图巴氏距离，定量评价各单元及其临近单元NDVI频数直方图相似性以得到有效遥感分区，以此分区为基础，建立NDVI-草地产草量统计模型，将其推广至临近具有相似性NDVI的区域，从而获取整个天山山区牧草(鲜重)总产草量遥感估算模型。

1 数据来源与研究方法

1.1 研究区概况

按草地植被型组归纳新疆草地有草原、荒漠、草甸、沼泽4个类组。草原草地除阿尔泰山南麓冲积扇上部和天山北麓山前倾斜平原存在平原荒漠草原片段外都发育在山地。荒漠草地在新疆分布广、面积大，包括温性草原化荒漠、温性荒漠和高寒荒漠类。草甸草地分布面积仅次于荒漠和草原草地，包括低地草甸、山地草甸和高寒草甸3种类型。沼泽草地分布面积不大，分别占新疆草地面积和可利用土地面积的0.47%和0.50%，主要分布在南北疆各大河流下游或河口、河漫滩、湖泊周围、山前潜水溢出带山间低位盆地的湖滨和山地高位盆地等，植被是由湿中生多年生根茎禾草和莎草类组成。由表1和图1可见，本研究区域草地在类型上具有一定的典型性和完整性，在空间分布上具有一定的连续性，可作为草地总产草量遥感估算的典型试验区域。

表1 研究区4类典型草地统计面积及比例Table 1 Area and proportion of four kinds of typical grasslands in the study area

图1 研究区4类典型草地空间分布图Fig. 1 Spatial distribution of four kinds of typical grasslands in the study area

1.2 数据来源

数据来源包含地面监测和卫星遥感监测数据。地面监测主要包括草地总产草量观测，观测集中在2009-2015年7月底至8月初(图2)。样地选择在相应群落典型地段且生境条件、植物群落种类组成、群落结构、利用方式和利用强度等具有相对一致性；样地间具有异质性，每个样地能控制的最大范围内地貌、植被等具有同质性，即地貌及植被生长状况相似。草原植被样地面积不小于100 hm2。矮小草本和小半灌木草地样方面积为草原1 m2，荒漠草原2 m2，草甸0.5 m2；株丛小的灌木和高大草本测产面积为4 m×4 m，灌丛或半灌木、高大草本分布均匀的地段设20 m×5 m或10 m×10 m样方。设置一个灌木样方，灌木样方内设置一个草本样方。每个样方要求测定植物名称、植物盖度和草地总产草量。遥感监测植被指数数据选用美国国家宇航局250 m的MOD13Q1归一化植被指数(NDVI)16 d合成产品，时间是2009年1月至2015年12月，轨道号h23v04、h23v05、h24v04和h24v05，共计644景影像。

1.3 研究方法

1.3.1植被指数处理 利用MODIS数据重投影工具(MODIS reprojection tools，MRT)对MOD13Q1数据集进行批量HDF-TIF格式和定义WGS84投影，应用GIS平台计算2009-2015年7月底至8月初多年NDVI均值和最大值合成法(maximum value composition，MVC)合成2009-2015年7月底至8月初最大NDVI，以表征对应年草地总产草量实测期内牧草平均和最好长势情况，基于GPS测点信息提取地面采样点对应NDVI。

图2 2009-2015年7月底至8月初牧草总产量实测位置分布图Fig. 2 Spatial distribution of the measured locations of the total production of the herbage from the end of July to the beginning of August during the years 2009-2015

1.3.2巴氏距离计算 以县(市)为研究单元，统计2009-2015年7月底至8月初时段多年NDVI均值频数直方图，利用巴氏距离(Bhattacharyya distance)测量各县(市)及其临近县(市)NDVI直方图分布的相似性。巴氏距离定义为[19-21]：

(1)

式中：H1与H2为县(市)尺度多年NDVI均值频数直方图，N为频数直方图组数，I为频数直方图组编号，H1(I)为县(市)频数直方图中组数I对应的频数。d(H1,H2)越接近1说明两个直方图越相似，即所在县(市)被划分为一个遥感分区。本研究中，定义d>0.5，即为相似区域。

1.3.3模型精度验证 应用留一交叉验证(leave-one-out cross validation，LOOCV)对草地总产草量遥感模型模拟结果进行验证[22]。假设有n个样本，从中选择一个观测值作为验证数据，其他n-1个样本作为训练样本来建立回归模型，用回归模型模拟选出的单个验证数据的预测值来检验模型精度，如此重复n次，用n个结果平均值来衡量模型模拟精度。每个模型预测能力由均方根误差(RMSE)及拟合系数R2决定。RMSE用来量化模型精度，R2用来评估模型的准确性。

(2)

式中：Fw(yi)表示第i个样方草地总产草量(鲜重)实测值，yi为模型模拟第i个预测值，n是观测样本总数。RMSE越小表明回归模型越精确，R2越接近于1表示模型精度越高。

2 基于巴氏距离的NDVI均值相似性分析和建模分区

各县(市)间NDVI均值巴氏距离计算结果如表2所列，将巴氏距离大于0.5的县市定义为相似县(市)，研究区域被分为7个建模区域(图3)。

各建模区域所含县(市)多年NDVI均值-频数直方图如图4所示，可见在地形、气候条件、土壤环境等方面具有很大相似性的区域应该具有相似的牧草生物量，从而使得该区域NDVI值也相近。而对于地形、气候条件、土壤环境等方面空间差异较大的地区，其牧草生物量也呈现较大的差异性，对应区域NDVI值也具有大的离散度。因此，本研究依据NDVI值的局部相似性准则进行牧草生物量建模分区是合理的。

对各个分区分别采用指数、线性、多项式以及幂函数4类函数建立植被指数与草地总产草量回归模型(表4)。结果显示：Ⅰ、Ⅲ、Ⅵ分区选择幂指数函数、Ⅱ、Ⅳ、Ⅶ分区选择指数函数、Ⅴ分区采用多项式进行建模具有最小的拟合误差。这可能是因为不同分区，气候、地形、牧草生长环境等因子存在差异，使得不同分区NDVI空间分布差异较大，从而使得产草量随NDVI的变化速率及趋势也不同。因此，利用NDVI进行建模分区，不同分区采用不同的拟合模型才能更精确地估算整个研究区的牧草产草量。

如图5a所示，分区Ⅰ采用幂指数模型拟合草地总产草量与NDVI之间的函数关系，效果最佳，能体现牧草鲜重随NDVI呈明显的快速递增现象。该分区内NDVI在0.25～0.85变化，表明该分区内大部分区域草地长势较好、产草量较高。图5b所示，分区Ⅱ采用指数模型拟合草地总产草量与NDVI之间的函数关系效果较好，能够很好地反映牧草鲜重随NDVI值增长的变化趋势。该分区内NDVI主要在0.04～0.65间变化。当少数子区域NDVI增加至0.7左右时，对应产草量骤然变大。当NDVI在0.04～0.50内时牧草鲜重随NDVI呈近似线性变化。图5c所示，分区Ⅲ采用幂指数模型拟合草地总产草量与NDVI之间的函数关系效果较好，能够反映牧草鲜重随NDVI的变化趋势，该分区内NDVI主要在0.10～0.47变化，但是，少数子区域产草量超出分区产量均值，异常大。如图5d所示，分区IV采用指数模型拟合草地总产草量与NDVI之间的函数关系，误差最小，该分区内NDVI值大都大于0.25。但是，存在一些子区域的产草量异常大。图5e所示，分区Ⅴ采用多项式拟合模型拟合草地总产草量与NDVI之间的函数关系，误差最小。该分区内NDVI值都大于0.05，牧草鲜重随NDVI值变化速度较慢。图5f所示，分区Ⅵ采用幂指数模型拟察布查尔Chabhar County，伊宁县Yining County，霍城县Huocheng County，伊宁市Yining City，尼勒克县Nilka County，巩留县Gongliu County，特克斯县Tekesi County，新源县Xinyuan County，昭苏县Zhaosu County，乌苏市Wusu City，博乐市Bole City，温泉县Wenquan County，精河县Jinghe County，石河子市Shihezi City，沙湾县Shawan County，玛纳斯县Manas County，呼图壁县Hutubi County，奎屯市Kuitun City，昌吉市Changji City，米泉市Miquan City，阜康市Fukang City，奇台县Qitai County，木垒自治县Mulei Autonomous County，吉木萨尔县Jimsar County，巴里坤县Barlikun County，伊吾县Yiwu County，和静Hejing County.

图3 基于巴氏距离的NDVI时间序列相似性分区示意图Fig. 3 Diagram of similarity zoning of NDVI time serial data based on the Bhattacharyya distance

图4 归一化NDVI均值-频数直方图Fig. 4 The normalized mean-frequency histogram of NDVI data in the study area

合草地总产草量与NDVI之间的函数关系，误差最小。该分区内NDVI值都大于0.1，当NDVI在0.1～0.3牧草鲜重随NDVI呈近似的线性增加。图5g所示，分区Ⅶ采用指数模型拟合草地总产草量与NDVI之间的函数关系是合适的，该分区内NDVI都大于0.3，主要集中在0.3～0.8，牧草鲜重随NDVI上升的速度较快。

3 讨论

3.1 与不分区建模对比

忽略NDVI的空间异质性，采用研究区内所有采样点建立NDVI-草地总产草量回归模型(表5)，与考虑NDVI空间异质性采用相似性分区的回归模型对比，以探究相似性分区回归模型在各个区域建模的提升效果。可以看出，除线性回归模型决定系数R2为0.389，其余回归模型的决定系数R2均在0.570～0.600，最高的幂指数模型为0.600。同时，除线性回归模型的估计误差RMSE大于2 000 kg·hm-2，为2 205.219 kg·hm-2外，其余回归模型的估计误差RMSE均在2 000 kg·hm-2以下，最低为幂指数模型，为1 478.640 kg·hm-2。结合R2与RMSE两个指标，可知适合全区域产草量估算的回归模型为幂指数模型。

表4 NDVI-草地总产草量回归模型及交叉检验结果列表Table 4 List of regression models and the results of cross validation of total production of the herbage-NDVI

图5 各个建模分区草地总产草量-NDVI拟合关系图Fig. 5 The fitting relation curves of the total production of the herbage-NDVI for each modeling zone

回归模型Modeling type回归系数Regression coefficientsR2调整R2 Adjusted correlationcoefficientRMSE/(kg·hm-2)样点数Sample numbery=aebx+ca=57.185; b=0.703; c=425.1560.5970.5961 492.482y=ax+ba=1 330.417; b=-2 619.2600.3890.3892 205.219y=a+bx+cx2a=3 820.722; b=-2 682.155; c=482.3200.5730.5711 680.967y=axb+ca=0.430; b=5.035; c=750.2530.6000.5991 478.640600

对比表4中各分区的最优模型以及表5中不分区情况下的最优模型，不分区回归模型的决定系数R2低于各个分区回归模型。说明采用NDVI相似性分区后，回归模型自变量能解释更多因变量的变异，更适合用于牧草产量的估算。同时，除分区Ⅰ外，其余分区回归模型的估计误差RMSE均在1 000 kg·hm-2以内，最低可到266 kg·hm-2，均优于不分区回归模型。但分区Ⅰ回归模型的估计误差达到了2 951.910 kg·hm-2，这主要是因为分区Ⅰ内实测草地总产量在10 000～30 000 kg·hm-2的样点多；而不分区回归模型估计误差低于分区Ⅰ的原因是研究区内其余区域的实测草地总产草量较低，低于10 000 kg·hm-2，由交叉验证原理可知，较低的实测值可将较高的模型估计误差平均。因此，分区回归模型能在各个区域内提升牧草产量的估算。

3.2 与随机分区建模对比

分区建模能够在各个区域内很好提升牧草产量的估算精度。根据采样点的聚集程度进行随机分区(图6)，并在各个随机分区内建立NDVI-草地总产草量回归模型(表5)，以探究相似性分区回归建模对牧草产量估算的提升。可以看出，随机分区Ⅰ的最佳回归模型为幂指数模型，其决定系数R2为0.549，估计误差RMSE为2 508.562 kg·hm-2(表6)；随机分区Ⅱ的最佳回归模型为抛物线模型，其R2为0.524，RMSE为479.584 kg·hm-2；随机分区Ⅲ最佳回归模型为指数模型，其模型决定系数R2达到了0.692，模型的估计误差为604.101 kg·hm-2。

结合表5和表6，随机分区后模型对因变量变异的解释程度提升不大，只有随机分区Ⅲ回归模型决定系数大于不分区回归模型，这表明分区建模并不总是带来好的估算结果。同时，模型估计精度RMSE也只在随机分区Ⅱ与随机分区Ⅲ低于不分区回归模型，表明分区建模对产草量估算的提升具有区域性。结合表4和表6，采用相似性分区后，回归模型在模型决定系数R2、模型估计精度RMSE两个指标上都有显著提升，除了分区相似性分区Ⅰ回归模型的RMSE。

图6 NDVI时间序列随机分区示意图Fig. 6 Diagram of random zoning of NDVI time serial data

随机分区 Random zone回归模型Modeling type回归系数Regression coefficientsR2调整R2 Adjusted correlationcoefficientRMSE/(kg·hm-2)样点数Sample numberⅠy=aebx+ca=333.549; b=0.485; c=-690.4080.5470.5442 516.493y=ax+ba=2 383.587; b=-7.565E+030.4700.4682 951.910y=a+bx+cx2a=4 929.915; b=-3 038.495; c=529.2840.5450.5412 552.954y=axb+ca=5.675; b=3.782; c=459.5520.5490.5452 508.562251Ⅱy=aebx+ca=21.182; b=0.720 6; c=686.1250.5290.522501.783y=ax+ba=366.143; b=-82.7850.4040.400540.456y=a+bx+cx2a=1 530.595; b=-662.416; c=133.3700.5240.517497.584y=axb+ca=0.533; b=4.449; c=776.7870.5300.524500.581149Ⅲy=aebx+ca=0.036 0; b=1.802; c=860.0360.6920.689604.101y=ax+ba=615.016; b=-389.0690.3220.318915.760y=a+bx+cx2a=2 441.925; b=-1 643.751; c=316.8540.5200.515821.791y=axb+ca=8.301E-07; b=11.965; c=882.2260.6890.686610.138200

4 结论

不同类型草地因其种类构成、长势、盖度、产量差异引起卫星植被指数的变化，使得充分利用卫星资料进行大面积周而复始地监测草地生产力动态变化成为可能[5]。遥感估算天山地区草地总产量是一项具有重要实用价值的工作。天山北坡草地垂直分异明显，类型复杂，给遥感识别分类和监测带来了很大的困难，本研究提出的NDVI相似性分区下天山山区草地总产草量遥感估算技术实用可行，为后期高级估算方法开发打下了良好基础，将研究区域分为若干个分区可以减少估产模型的异质性，使得回归模型中的误差项完全符合理想假设(误差独立同分布期望为零，从而提高模型对数据的整体拟合度)。根据NDVI直方图显著的空间差异性特征和邻近区域特征相似性，分区域对2009-2015年内7-8月牧草总产量(鲜重)和NDVI进行单一响应因子建模，发现不同建模区域拟合方程形式也不同，总产草量与植被指数具有线性、幂指数、指数关系、一元二次回归方程多种拟合关系。总体来看，7个建模分区NDVI-草地总产草量拟合相关系数在0.836～0.784。交叉检验除天山北坡西段-伊犁河谷草原畜牧业区RMSE值为2 951 kg·hm-2外，其他分区RMSE值在266～928 kg·hm-2，原因在于伊犁河谷草原畜牧业区实测草地总产量在10 000～30 000 kg·hm-2的样点多。

后续研究中，会在NDVI分区建模的基础上，综合考虑NDVI分布的巴式距离分区和NDVI空间分布(NDVI空间结构信息)，增加确定性自变量的一般回归模型和带有观测误差自变量的扩展回归模型(遥感反演误差到回归统计误差的传播)的比较分析，加入草地长势其他有关解释变量进行建模，设计区域差异、区域效应在回归系数中，如空间滤波和空间面板回归模型里的空间或个体效应变量。