伊犁地区植被生产力的时空分异及其与地形因子关系研究

张振东，昝 梅

(1.新疆师范大学地理科学与旅游学院，乌鲁木齐 830054；2.山东师范大学地理与环境学院，济南 250358)

陆地生态系统是碳元素重要的源和汇，在全球碳循环中扮演了至关重要的角色，而植被作为陆地生态系统的主体在保证全球陆地生态系统中的物质循环和能量流动的动态平衡中发挥了重要的作用[1].植被净初级生产力(net primary productivity，NPP)是指绿色植物通过光合作用产生的全部有机物量扣除自身生长所必须的有机量后的剩余部分[2].由于NPP能够以统一的标准反应生态系统生产力，可以很好的反映植被对自然资源的利用能力[3-4].因此，从19世纪90年代开始植被净初级生产力的研究开始得到各国的广泛关注，遥感技术也开始迅速发展.遥感技术在获取地面植被动态变化方面具有时间分辨率高、精确度高的优势，很快成为研究植被净初级生产力的主流方式.CASA(carnegie ames stanford approach)模型也随着遥感技术的进步而兴起，CASA模型的优势是利用遥感数据进行估算，参数少，稳定性强，能够实现动态更新植被净初级生产力[5].

伊犁地区是丝绸之路经济带的重要通道，地域优势明显.但随着气候变化和人类经济活动影响的增加，绿洲冷岛效应减弱、冰川面积退缩、草场退化和荒漠化等生态环境问题日趋严重，威胁着绿洲生态安全和生产的可持续发展[1，6].研究伊犁地区的植被净初级生产力有助于更好的掌握伊犁地区植被和绿洲生态系统的变化情况，同时为以NPP为基础的其他有关伊犁地区的科学研究提供理论基础[3].目前已经有很多学者利用不同模型对中国各类生态系统的NPP进行反演.如朱文泉等采用改进的最小二乘法对中国典型植被的最大光利用率进行了系统的模拟，并针对不同植被分类精度可能带来的误差对最大光利用率进行了敏感性分析[7]；艾则孜提约麦尔·麦麦提利用MODIS数据采用GWR(地理加权回归)建模法对2000年—2014年博斯腾湖流域NPP时空变化特征及影响因子进行了分析[8]；另外有多个地区和省份利用CASA模型求取植被的净初级生产力，如东北三省、四川、河南、福建、和新疆等[1，4，9-10].但是，某些在社会、经济和生态方面具有重要地位的热点区域的植被生产力及影响因素的研究却未被充分重视，如伊犁地区.因此，本文利用CASA模型，以伊犁地区为研究区，利用MODIS数据、日照数据和DEM数据运用CASA模型对伊犁地区NPP进行估算，进而对伊犁地区2012年—2014年植被NPP时空格局变化特征进行探究；分析伊犁地区不同植被类型NPP的特征；并阐明不同地形条件下伊犁地区NPP变化特征；最后通过对比验证分析CASA模型在研究区的可用性和估算结果的可信程度[11].

1 材料与方法

1.1 研究区概况

伊犁哈萨克自治州地处新疆维吾尔自治区天山北部的伊犁河谷内(如图1).位于80°09′42″～84°56′50″E，42°14′16″～44°53′30N″.地形三山环绕，地貌类型复杂多样，主要划分为山地、丘陵、平原和山间谷地四类；海拔高程为481～6 308 m .气候属于寒温带半干旱大陆性气候；植被类型主要以草地、农作物、稀疏植被为主，分别占研究区总面积的60.02%、27.01%、7.37%，其他植被类型仅占5.62%.植被覆盖度变幅大，受降水和气温的影响，山地区域植被覆盖度呈上升趋势，伊犁地区是新疆植被覆盖度最高的地区，平均覆盖度约90%.除雪线以上的常年冰雪带和雪线以下的寒冻裸岩带外，其余地面均有植被.低山带为草地，中山带为林地，高山带为夏季优良草场[12].

图1 研究区概况图Fig.1 Overview of the study area

1.2 数据来源

1.2.1 MODIS的NDVI数据 本文选用的MODIS数据来源于美国国家航空航天局网站(http://ladsweb.nascom.nasa.gov/).下载了2012年—2014年的MOD13A1产品，共144幅图像.空间分辨率为500 m，时间分辨率为16 d.采用最大值合成法(MVC)得到以月为周期的伊犁地区NDVI数据.

1.2.2 MODIS的土地覆被数据 本文选用MODIS土地覆盖类型产品(land cover data)MCD12Q1，空间分辨率为500 m，时间分辨率为1 a.根据伊犁地区植被类型实际分布状况，最终选择植物功能型方案作为研究区土地覆被的分类方案将伊犁地区的植被类型分为针叶林(常绿针叶林、落叶针叶林)、阔叶林(落叶阔叶林)、草地、农作物(谷类作物、阔叶作物)、稀疏植被五种类型.通过以上方法我们分别获得了研究区2012年、2013年和2014年的植被类型空间分布图(如图2).

1.2.3 气象数据 从中国气象数据网(http://data.cma.cn/)下载得到伊犁地区周围11个气象台站2012年—2014年月平均气温、月降水量数据和日照时数的日值数据.通过克里格(Kriging)插值法生成500 m空间分辨率的月均气温和月降水量的序列栅格数据.将日值日照数据累加，计算得到每个站点逐月日照百分率后采用Kriging插值法获得研究区逐月日照百分率的空间时序栅格数据，再利用太阳辐射能计算公式得到2012年—2014年逐月的太阳总辐射能序列栅格数据[13].

1.3 CASA模型的构建

本文选用的CASA模型是由Potter等基于Monteith的理论提出的.该模型是基于植被的机理过程建立的，利用APAR和ε来模拟植被的NPP[14-16]，所需参数相对较少且易于获取.有效的避免了由于参数缺乏而人为简化或估算产生的误差[19].相对于其他模型的优势是利用遥感数据进行估算，覆盖的范围广；数据更新快；时间分辨率高；能够实现动态更新植被净初级生产力[5]，该模型的具体公式为：

NPP=APAR(x，t)×ε(x，t)，

(1)

式中，NPP表示植被净初级生产力，APAR(x，t)表示植被吸收的光合有效辐射，ε(x，t)表示光能利用率.x表示像元的空间位置，t表示时间(月).

APAR主要由SOL和FPAR两部分决定[5].APAR(x，t)的估算公式为(2)：

APAR(x，t)=SOL(x，t)×FPAR(x，t)×0.5，

(2)

式中，SOL(x，t)表示像元x在t月份的太阳总辐射能；FPAP(x，t)表示植被层对入射光合有效辐射的吸收比例；其中0.5是常数，它表示植被能够吸收利用的太阳有效辐射量.

图2 2012年—2014年伊犁地区植被类型空间分布图(a为2012年、b为2013年、c为2014年)Fig.2 Spatial distribution map of vegetation types in Yili area from 2012 to 2014 (a is 2012，b is 2013，and c is 2014)

相对于地球辐射来说，SOL的波长较短，故又称为短波辐射.SOL是利用日照时数采用经验公式进行的估算.公式如下[13]：

(3)

式中a，b是计算太阳总辐射量的参数，其中，a=0.353，b=0.543[2].n/N为日照百分率，Ra为天文辐射.Ra由日地距离d，纬度φ，太阳日落角ω，太阳常数(0.082 MJ·m-2·min-1)，太阳赤纬σ和天数t通过公式(4)计算得到：

(ω·sinφ·sinσ+cosφ·cosσ·sinω)，

(4)

其中

ω=arccos(-tanφ·tanσ)，

根据上述公式可以计算出2012年—2014伊犁地区11个台站逐月天文辐射.再利用ArcGIS软件得到空间插值后的研究区天文辐射逐月栅格数据集.

根据文献资料可以得到FPAR与植被类型和NDVI存在一定的关系.通过Kumar和Monteith等人的验证，提出利用FPAR和NDVI建立的线性关系来估算FPAR，估算结果具有较高的精度[14，17]，具体公式为：

(FPARmax-FPARmin)+FPARmin，

(5)

式中，NDVI(x，t)是指像元x在t月份的NDVI值.NDVI(i，max)是指i种植被类型NDVI的最大取值.NDVI(i，min)是指i种植被类型NDVI的最小值.FPARmax和FPARmin是指FPAR的最大值和最小值，一般为常数，即FPARmax=0.950，FPARmin=0.001[13].其中，NDVI(i，max)和NDVI(i，min)与植被的类型有关，是把不同植被类型每月NDVI的值从小到大排序累积后选取5%和95%处的值分别作为该月的NDVI最小和最大值，详见表1.

表1 各年份每种植被类型NDVI最大值最小值统计Tab.1 Statistics of maximum and minimum values of vegetation

光能利用率指植被层吸收入射光合有效辐射并将其转化为有机物的能力，与温度和水分等因子有关[15-16]．

ε(x，t)=T1(x，t)×T2(x，t)×W(x，t)×εmax，

(6)

式中，T1(x，t)、T2(x，t)是温度胁迫因子，W(x，t)是水分胁迫因子.温度和水分胁迫因子可以参考文献研究结果得到[18-20].εmax是最大光能转化率，即在理想状态下植被吸收光合有效辐射并转化为有机碳的最大光能转化率[21]，一般与植被类型有关，本文主要参考朱文泉和冯益明等模拟的εmax结果求取平均值后获得研究区主要植被类型的εmax[7，22-25]，详见表2.

表2 各植被类型最大光能利用率Tab.2 Maximum light energy utilization rate of each vegetation type

本文利用MODIS数据(MOD13A1，MCD12Q1)、DEM数据和气象数据，依照上述模型构建方法对数据进行处理，生成2012年—2014年月尺度NDVI数据、太阳辐射数据以及年尺度植被分类数据，然后进行CASA模型APAR、FPAR、SOL、ε等参数的估算，以此来研究伊犁地区植被生产力的时空分异及其与地形因子的关系，模拟流程如图3所示.

图3 植被净初级生产力模拟流程Fig.3 Simulation process of vegetation net primary productivity

2 结果与分析

2.1 伊犁地区2012年—2014年NPP空间变化特征

模拟结果表明，研究区NPP的空间分布具有较大的差异.由图4可以看出，2012年—2014年伊犁地区的植被NPP的年平均最大值为868.5 gC·m-2·a-1，年平均最小值为10 gC·m-2·a-1，总体平均值为200 gC·m-2·a-1.研究区2012年—2014年植被NPP的分布总体呈现出东部大于西部的特征.NPP的高值区(>500 gC·m-2·a-1)主要在伊犁东部及北部和南部的山地边缘地区呈半环状分布.中值区(200～500 gC·m-2·a-1)主要集中在中部河流流经的水资源较为丰富的地区.低值区(<200 gC·m-2·a-1)主要集中在伊犁中部植被覆盖率低、城镇较为集中的地区以及海拔较高的山区.

此外，结合2012年—2014年伊犁地区植被类型图(图2)可以看出NPP的高值区在针叶林、阔叶林和草地分布的地区较为集中，中值区在农作物分布的地区较为集中，低值区在稀疏植被分布的地区较为集中.以上分布特征也说明CASA模型估算得到的伊犁地区NPP的结果较为合理.

图4 2012-2014年NPP空间分布(a为2012年、b为2013年、c为2014年)Fig.4 NPP spatial distribution from 2012 to 2014(a is 2012，b is 2013，and c is 2014)

2.2 伊犁地区2012年—2014年不同植被类型NPP月平均变化特征

由2012年—2014年不同植被类型逐月平均NPP变化趋势(图5)可以看出，研究区植被NPP随着季节的变化呈现出周期性变化特征.不同植被类型的NPP随时间变化的周期和大小不同，具有较好的区分度.

不同植被类型的逐月平均NPP呈单峰型曲线特征.各种植被类型的生长季均集中在每年的3～11月.各种植被类型的平均NPP基本都在2月达到最小值，7月达到最大值.另外，不同植被类型NPP的月变化趋势基本相同，即3～7月NPP处于增长期，其中3～6月增长速率较快，6～7月增长速率较慢.各种植被类型NPP在7～12月为下降期，7～11月下降速率较快.

通过比较发现2012年—2014年不同植被类型逐月NPP平均值的变化幅度存在较大的差异.阔叶林月均NPP的变化速率和变换幅度最大.农作物、针叶林和草地月均NPP的变化速率和变化幅度差距不大，处于中等水平.稀疏植被变化速率和变化幅度最小.一般而言林地的固碳能力在同等水热条件下大于草地和农作物，所以导致林地(阔叶林、针叶林)的月均NPP较其他植被类型大.而农作物和草地的NPP分别受到人为灌溉和物候、气候环境的影响更大.由图2研究区植被类型空间分布图可知，稀疏植被主要分布在伊犁地区的南部高海拔山区，受水热限制导致该种植被类型的月均NPP普遍较小.

图5 2012年—2014不同植被类型逐月平均NPP变化趋势Fig.5 Trends of monthly average NPP of different vegetation types from 2012 to 2014

2.3 不同植被类型对NPP的贡献率

结合研究区植被类型空间分布图(图2)可知，伊犁地区针叶林和阔叶林主要分布在伊犁地区的南部山地边缘，约占植被分布总面积的5.62%.草地面积占植被分布总面积的60.02%.农作物则主要集中在伊犁河和喀什河流经的水资源较为丰富的地区，占植被分布总面积的27.01%.草地和农作物所占面积超过植被分布总面积的87.02%.本文利用伊犁地区不同植被类型年平均NPP和不同植被类型分布面积的乘积在整个伊犁地区所有植被类型年平均NPP总和中所占的比例作为不同植被类型对该地区NPP的贡献率.由表3可以看出，研究区草地NPP的贡献率最大为58.59%；第二是农作物，对伊犁地区NPP的贡献率为31.45%；针叶林和阔叶林对伊犁地区NPP的贡献率基本相同，分别为3.75%、3.60%.最小的是稀疏植被，对伊犁地区NPP的贡献率仅为2.16%.

表3 不同植被类型对NPP值的贡献率Tab.3 Contribution rates of different vegetation types to NPP values

2.4 地形因子与伊犁地区植被NPP的关系

2.4.1 海拔高度与伊犁地区植被NPP的关系 利用伊犁地区DEM数据对伊犁地区NPP的空间分布进行分析.由图6可以看出，伊犁地区年均NPP总体上呈现随海拔增加而先增加后减少的特征，年均NPP的中值区分布在低海拔地区(<1 500 m)；高值区分布在中海拔地区(1 500～3 000 m);低值区分布在高海拔地区(>3 000 m).具体来说：一方面2012年—2014年均NPP随海拔变化的梯度性明显，年均NPP为0～100 gC·m-2·a-1主要集中在3 500～5 500 m；年均NPP为100～200 gC·m-2·a-1主要集中在海拔504～1 000 m和3 500～4 000 m；年均NPP为200～400 gC·m-2·a-1主要集中在海拔1 000～1 500 m和3 000～3 500 m；年均NPP为500～600 gC·m-2·a-1主要集中在海拔1 500～3 000 m.另一方面：年均NPP为400～500 gC·m-2·a-1在海拔504～3 500 m范围内分布较均匀，在不同海拔梯度内所占的比例差距较小，分别为27.15%(500～1 000 m)、20.74%(1 000～1 500 m)、28.31%(1 500～2 000 m)、23.77%(2 500～3 000 m)、20.93%(3 000～3 500 m)；年平均NPP为0～100 gC·m-2·a-1在海拔高于3 500 m时占不同海拔梯度内的比例均为最高，分别为52.18%(3 500～4 000 m)、93.43%(4 000～4 500 m)、99.47%(4 500～5 000 m)、100%(5 000～5 500 m)，这与研究区年均NPP随海拔的变化情况基本一致.另外，由图6可以看出，多个年均NPP梯度(500～600 gC·m-2·a-1、600～700 gC·m-2·a-1、700～800 gC·m-2·a-1)在海拔2 000～2 500 m处均呈现出增长的趋势.出现这种现象的原因主要是由于海拔较低的地区人口稠密，天然植被生长受人类影响强烈，破坏严重，但农作物大部分生长于此，光热充足并且膜下滴管提供充足水分，年平均NPP处于中值水平.随着海拔的升高，在中山带形成降雨带，天然植被生长茂盛，年平均NPP最高.但随着海拔的持续升高降水减少，温度降低，土层稀薄，土壤肥力差，年平均NPP降到最低，基本保持在0～100 gC·m-2·a-1范围内.

图6 2012年—2014年平均NPP随高程的变化Fig.6 Average NPP changes with elevation in 2012-2014

2.4.2 坡度与伊犁地区植被NPP的关系 图7为2012年—2014年平均NPP随坡度变化图，可以看出年平均NPP呈现随坡度增加而先增加后减小的特征.当坡度小于35°时，研究区植被年均NPP 400～500 gC·m-2·a-1主要集中在坡度0～5°；年均NPP为600～700 gC·m-2·a-1的植被主要集中在5～25°坡度范围内；年均NPP为0～300 gC·m-2·a-1的植被则主要分布在坡度大于35°范围内，另外，年均NPP 500～600 gC·m-2·a-1在0～35°坡度范围内分布较均匀，在不同坡度梯度内所占的比例差距较小，分别为29.84%(0～5°)、29.82%(5～15°)、28.53%(15～25°)、24.27%(25～35°).由于坡度较低的地区地势平坦，水热条件好，适宜农作物的生长，导致植被年均NPP保持在中值水平；随着坡度的增加，人类干扰减少，林地面积增大，年均NPP增大；当坡度大于35°后，土壤侵蚀加剧，土层薄，土壤肥力差，不利于植被的生长，导致植被年均NPP降低.

图7 2012年—2014年平均NPP随坡度的变化Fig.7 Average NPP changes with slope in 2012-2014

3 讨论

3.1 结果验证

通过查阅文献可知MODIS的NPP产品已经被许多学者验证过结果较好，而且部分地区结果精度较高[16].用MOD17A3 NPP数据对研究区利用CASA模型的估算的植被NPP结果进行对比验证，结果如图8所示.选用R2作为趋势线拟合程度的指标，一般在0～1范围内取值，值越大拟合程度越高.通过对2012年—2014年CASA模型NPP估算值和MOD17A3模拟数据进行比较发现2012年—2014年的R2值分别为0.87、0.86、0.95.说明利用CASA模型估算伊犁地区的NPP与MOD17A3的结果一致性较好，也说明利用CASA模型对伊犁地区植被NPP估算实用性较强.

图8 2012年—2014年CASA模型估算值与MODIS产品NPP值比较Fig.8 Comparison of 2012-2014 CASA model estimates with NPP values of MODIS products

3.2 植被NPP的时空格局

植被NPP的空间分布受地形、土壤、水分和植被类型等因素的影响，具有很强的空间异质性.通过CASA模型模拟结果表明，伊犁地区的植被NPP分布总体呈现出东部大于西部的特征.综合伊犁地区地形及植被类型分布情况分析发现(图2)，伊犁东部多为农作物、林地、草地，该区域位于伊犁河下游，区域水热条件好，草地质量好，面积广阔，并且地形较平坦，利于人为灌溉农田系统的植被生长，加上伊犁北部和南部山地的边缘地区林地面积广阔，同等水热条件下林地的固碳能力要大于其他植被类型；从而导致上述地区的植被NPP较大；中部地区植被类型主要为农作物，植被NPP受水分条件影响大，该区域河流广布，灌溉水源丰富，水土条件优良，因而植被NPP处于中值水平；而南部地区受地形、土壤类型的限制，以山地和荒漠为主，植被稀疏，导致该区域植被NPP较低；伊犁西部地势相对较低，水热条件相对不足，导致该区域植被类型单一的草地NPP明显很低.研究区不同植被类型的年内月均NPP呈单峰型曲线主要是受该区域植被物候周期影响所致.

3.3 植被NPP与地形因子的关系

海拔和坡度是主要的地形指标，伊犁地区年均NPP呈现随海拔的增加而先增加后减少的特征.年均NPP的中值区分布在低海拔地区(<1 500 m)；高值区分布在中海拔地区(1 500～3 000 m);低值区分布在高海拔地区(>3 000 m).研究区受山盆结构分异造成的水热条件限制明显，植被年均NPP也呈现出随坡度增加先增加后减小的特征.植被分布集中在高程小于3 000 m且坡度在0～35°的区域内，此区域水热条件也适合植被的生长，因此在此区域植被NPP年均值最高.而在海拔大于3 000 m并且坡度大于35°的区域内土壤侵蚀加剧，土层薄，土壤肥力差，加之水热条件不足，不利于植被的生长，因此植被年均NPP在此区域出现最低值.

4 结论

本文利用美国国家航空航天局网站上以月为周期的NDVI数据集MODIS13A1、以年为周期的植被类型数据集MCD12Q幅以及各种气象数据，利用CASA模型估算得到了以月为周期，空间分辨率为500 m的伊犁地区植被NPP.进而分析了研究区2012年—2014年植被NPP的时空分布特征并探讨了不同植被类型对NPP的贡献率.主要结论如下.

1) 不同植被类型的月均NPP呈现3～7月增加，7月到11月下降的趋势.2月各种植被类型的NPP达到最小值，7月达到最大值；但不同植被类型月均NPP的变化速率和变换幅度存在较大差异.阔叶林月均NPP的变化速率和变换幅度均为最大.农作物、针叶林及草地月均NPP的变化速率和变化幅度差距不大，在所有植被类型中变化速率和变化幅度居中.稀疏植被月均NPP变化速率和变化幅度均最小.

2) 不同植被类型对伊犁地区NPP的贡献率不同.其中草地的贡献率最大为58.59%；其次是农作物，贡献率都为31.45%；针叶林和阔叶林对伊犁地区NPP的贡献率基本相同，分别为3.75%、3.60%.最小的是稀疏植被，对伊犁地区NPP的贡献率仅为2.16%.

3) 伊犁地区年均NPP随海拔和坡度的增加均呈现出先增加后减少的特征；通过验证说明，利用CASA模型对伊犁地区植被NPP估算的结果与MOD17A3产品一致性较好.

气候尤其是水热因子对植被净初级生产力的影响极为显著，植被本身的物理特性对植被NPP的影响同样不可忽略，但本研究仅仅针对地形因子与植被净初级生产力的关系进行研究，具有一定的片面性；另外，本研究缺乏实地植被NPP的采样验证.这些问题都有待于进一步深入探讨.