基于净初级生产力变化的草地退化监测研究

吴 红，安 如，李晓雪，曲春梅，陆 玲，杨仁敏，龚天宇

(河海大学地球科学与工程学院，江苏 南京 210098)

草地退化是由于人为活动或不利自然因素所引起的草地(包括植物及土壤)衰退，生产力、经济潜力及服务功能降低，环境变劣以及生物多样性或复杂程度降低，恢复功能减弱或失去恢复功能的过程[1]，其实质是草地生产力的下降。据调查，三江源区的草地已呈现全面退化的趋势。其中，中度以上退化草地面积达120 000 km2，占本区可利用草地面积的58%。黄河源区作为三江源区重要生态区域，近年来草地生态系统的破坏非常严重。玛多县地处黄河源头，是黄河重要的水源涵养地，其生态作用对整个三江源区域有着直接重要的影响。在20世纪70年代，玛多县由于草地资源丰富，畜牧业发达,曾是有名的富裕县，但是随着自然因素及人为因素的变化和破坏，草地资源急剧减少，畜牧业受到很大影响，导致现在的玛多沦为全国有名的贫困县[2-4]。对该区域进行草地退化遥感监测，为相关部门提供有效的草地资源变化数据，对于整个生态环境建设和人民生活水平的提高都具有重要意义。近年来，草地资源的监测方法得到很大的改进，遥感影像的应用使草地资源的监测范围更大，也更具时效性。李昀等[5]以TM影像为数据源，利用归一化植被指数(NDVI)构建草地产草量模型，评价了景泰县草地资源现状；李文娟等[6]利用MODIS植被指数数据建立天然草地资源产量反演模型，估算了青海省草地草地资源载畜量。

目前基于遥感的草地退化信息提取主要方法有：基于草地退化监测指标与植被指数的相关性，建立退化指数的反演模型；基于遥感影像的草地退化信息直接目视综合判读法。前者受到遥感可反演参数与草地退化指标不能直接对应的限制，信息提取的精度差；后者对遥感影像的选取、时相的选择和后期处理具有严格的要求，而且要求解译人员对研究区草地退化状况有比较完整的认识，分析精度较高[7]。对于草地退化指标的选取，李辉霞和刘淑[8]基于NDVI建立了西藏自治区草地退化评价模型。高清竹等[9]选择草地的植被盖度为草地退化的遥感监测指标，研究了藏北地区草地退化的时空分布特征。王根绪和程国栋[10]研究发现，江河源区6县(达日、玛多、玛沁、治多、曲麻莱、杂多)草地退化面积占草地总面积的34.34%，其中草地退化及重度草地退化多发生在玛多和达日，表明黄河源区草地退化比较严重。草地退化的实质是草地生产力的下降，而净第一性生产力是表征草地生产力状况的一项重要指标。本研究采用基于草地净初级生产力的方法，对玛多地区的草地退化信息进行提取。

植被NPP的估算模型主要有：气候统计模型、过程模型和光能利用率模型3种。气候统计模型通过建立NPP与地面气候资料(如年降水量、年平均温度及年蒸散量等)的统计关系来估计NPP；过程模型从机理上对植物的生物物理过程以及影响因子进行分析和模拟，更能揭示生物生产过程，但模型复杂，涉及的参数多而且不易获得；光能利用率模型是基于资源平衡观点来估算NPP。模型形式简单且可以直接利用遥感数据，克服了气候统计模型以点带面的缺点，比较适合快速、大范围地监测草地NPP。本研究以光能利用率模型CASA(Carnegie-Ames-Stanford Approach)[11]为基础，考虑温度和大气水汽对植被的影响，构建基于TM数据的草地NPP的遥感反演模型。在此基础上通过计算草地潜在NPP作为未退化草地的NPP，获取草地NPP减少的相对百分数，以此作为监测草地退化的指标，对研究区的草地退化状况进行监测。

1 研究区概况

玛多县位于青海省的南部，三江源的中东部(图1)，位于96°50′～99°32′ E，33°52′～35°39′ N，全县土地总面积约24 500 km2，行政区划上分为3个乡(扎陵湖乡、黄河乡和黑河乡)和1个镇(花石峡镇)。玛多县属高寒半干旱气候区，年均气温-4℃，全县平均海拔在4 000 m以上，地形起伏不大，地势平坦，多湖泊、沼泽地、平坦地[12]。多年平均降水量321 mm，但年际变化大，最多的年份434.18 mm，最少的年份84.0 mm。玛多县主要土地覆被类型为草地，全县草地面积约22 526 km2。草地类型包括高寒草甸、高寒沼泽草甸和高寒草原类草地。沼泽化草甸具有典型的高原湿地生态特征，是区域内最具有生物多样性的草地生态系统之一[13]；高山草甸亚类草地分布在该县山地阳坡、阴坡、圆顶山、滩地和河谷高阶地上。高寒沼泽草甸亚类草地分布在排水不畅的平缓滩地、山间盆地、碟形洼地、湖畔、河滨、高山鞍部、山麓溢出带。高寒草原类草地常分布在滩地、宽谷、高原湖盆外缘、河流高阶地、剥蚀高原面以及干旱山地[14-16]。

图1 三江源区玛多县位置示意图

2 研究数据与方法

本研究采用玛多县2008年夏季(草场生长最好时期)Landsat TM影像共5景(影像获取时间为2008年6月12号)，并辅以玛多县草地类型分布图、土地覆被类型图、玛多县地形图和该区域的气象数据(水汽压、日照百分率及地表温度)。首先将遥感影像进行预处理：分别对5景影像进行假彩色合成，利用玛多县1∶10万地形图对影像进行几何精度纠正，在ENVI中利用FLASSH模块进行大气辐射校正，然后将影像拼接再根据玛多县行政边界进行剪裁，得到玛多县整个行政区域影像。利用光能利用率模型CASA模型进行草地NPP反演。中华人民共和国国家标准[17]中提出的草地退化评价依据：“未退化草地以监测点附近相同水热条件草地自然保护区中合理利用示范区相同草地类型的植被特征与地表、土壤状况为基准”，实质上就是以监测区域附近相同水热条件、相同草地类型中长势最好的草地状况为基准。潜在NPP是指草地在水热等生长条件适宜、没有自然或人为因素破坏情况下的净第一性生产力，它是评价草地退化程度的基准与依据。因此，可以把监测区域附近相同水热条件、相同草地类型草地的最大NPP作为该区域的潜在NPP，即未退化草地的NPP值[18]。

2.1草地退化指标及等级划分 草地退化遥感监测指标选取，既要考虑草地退化以及遥感原理和数据特点,又要有科学性、系统性、相对独立性、可行性和可操作性。在参考文献[7,18-19]及综合研究区可获取数据的基础上，采用草地NPP减少的相对百分数作为评价草地退化的指标：

(1)

式中，D为草地NPP减少的相对百分数；NPP潜在为草地的潜在NPP，即草地没有发生退化时的NPP。草地退化指数D取值为0～1，值越大说明草地退化越严重。参照中华人民共和国国家标准[17]提出的草地退化程度分级标准，将退化指数值分为4个等级：未退化0～10%、轻度退化10%～20%、中度退化20%～50%和重度退化>50%[18]。

2.2草地NPP的反演 国内许多学者利用CASA模型进行了NPP估算研究，如朴世龙等[19]利用CASA模型估算了1997年中国的NPP及其分布；陈利军等[20]使用1990年每旬的NOAA/AVHRR数据，采用光能利用率模型对中国该年每旬的NPP进行了估算，然后进行累加得到了1990年的NPP分布，并分析了其空间分布规律；朴世龙和方精云[21]利用1982-1999年NOAA/AVHRR数据及其相应的气象数据，结合植被和土壤质地等资料，应用CASA模型研究了青藏高原植被的NPP及其动态变化。

本研究以研究区2008年的地面气象数据(气温、日照百分率、地面水汽压等)和TM遥感数据作为数据源，气象数据是从中国气象科学数据共享服务网(http://cdc.cma.gov.cn/)下载。NPP可以由植物吸收的光合有效辐射(APAR)与光能转化率(ε)2个变量来确定。

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

(2)

式中，t表示时间，x表示空间位置；APAR(x,t)表示像元x在t时间吸收的光合有效辐射；ε(x,t)表示像元x在t时间的实际光利用率。

2.2.1光合有效辐射 植被所吸收的光合有效辐射取决于太阳总辐射和植被对光合有效辐射的吸收比例，用下列公式表示：

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

(3)

式中,SOL(x,t)是t月份象元x处的太阳总辐射量(MJ/m2)；FPAR(x,t)为植被层对入射光合有效辐射(PAR)的吸收比例；常数0.5表示植被所能利用的太阳有效辐射(波长为0.4～0.7 μm)占太阳总辐射的比例。光合有效辐射比例与比值植被指数(VI)存在以下关系：

(4)

式中，FPARmax=0.950，FPARmin=0.001；VImax和VImin分别表示草地的最大和最小比值植被指数，分别取值为1.08和4.41[22]。

由于所获取的地面气象数据没有太阳总辐射SOL(x,t)，需要对其进行计算，公式如下：

SOL(x,t)=Q0(a+bS)

(5)

式中，a和b是常数，分别取0.248和0.752；S为日照百分率，由气象数据获得；Q0表示最大晴天总辐射量，由地理纬度、海拔高度及平均地面水汽压估算。

Q0=0.418 675(C0+C1φ+C2H+C3e)

(6)

式中，φ为地理纬度；H为海拔高度(m)；e为地面水汽压(hPa)，由玛多气象站地面观测数据可获取；7月份：C0=22 397.100、C1=26.478、C2=0.876、C3=-11.659，4个回归系数值在不同月份其取值也不相同[23]。

2.2.2光能转化率(ε)的确定 光能转化率是指植被把所吸收的入射光合有效辐射转化为有机碳的效率。一般认为在理想条件下植被具有最大光能转化率，而在现实条件下光能转化率主要受温度和水分的影响，用式(7)表示：

ε(x,t)=σE×σT×ε*

(7)

式中，σE和σT分别为大气水分含量和温度对植物生长的影响系数[24]，ε*是理想条件下的最大光能转化率。CASA模型中采用的最大光能利用率0.389 g/MJ，并不适用于中国的植被[25]。Running等[26]根据生态生理过程模型BIOME-BGC模拟了不同的植被覆盖类型的最大光能利用率，其模拟结果更符合实际情况。本研究中采用Running等对草地的模拟结果0.608 g/MJ作为草地的最大光能利用率。温度对植物生长的影响系数，可以通过下式计算[24]。

式中，Ts为近地层气温或地表温度，由气象数据直接获得。

大气水分含量对植物生长的影响系数[24]：

σE=1.2exp(-0.35Ev)-0.2

(9)

式中，σE为大气水汽压差对草的光能利用率影响系数；Ev为大气水汽压差,由下式计算：

式中,Ts为近地层气温或地表温度(K)；Td为近地层露点温度(K)，可根据整层大气水汽含量w(g/cm2)来计算：

(11)

w=α0+α1e+α2e2

(12)

式中，α0、α1、α2为经验系数;e为地面水汽压(hPa)。对于我国大部分地区,经验系数α2≈0，α0、α1的计算公式因海拔高度和地理纬度的不同而不同,玛多地区α0=-0.005 6、α1=0.167 7[27]。根据2.2介绍的方法，得到玛多县2008年夏季草地NPP分布图(图2)。

图2 玛多县2008年夏季草地NPP

2.3草地潜在NPP 本研究把监测区域附近相同水热条件、相同草地覆被类型草地的最大NPP值作为该区域的潜在NPP。对于子区域的划分，经研究发现，中国草地类型图在制作过程中考虑水热条件同等性因素，并且其中型一级草地面积较小、属于同种草地类型。因而实际应用中可直接根据草地类型图进行子区域的划分[18]。

具体步骤为：1)为消除由于遥感数据噪音，或者混有其他类型植被等因素造成的NPP异常值，需对所得NPP图像进行5×5窗口中值滤波，以去除异常值[18]；2)本研究中采用的草地分类图具体到型一级， 在ERDAS中利用Vector模块将草地覆被类型图与草地NPP图叠加，计算草地类型分布图中每个多边形的最大NPP值，得到潜在NPP分布矢量图；3)将潜在NPP矢量图转换成与TM一致的30 m×30 m分辨率的栅格图，即是草地潜在NPP。

3 结果与分析

3.1草地退化结果分析 玛多县草地退化等级分布如图3所示。2008年玛多县全区草地退化状况仍然比较严峻，全区草地以中度退化(D值20%～50%)为主，面积为9 535.13 km2，占全区草地面积50.13%；轻度退化(10%～20%)草地面积占22.34%；重度退化面积3 256.34 km2，达17.18%，主要分布于北部地区(表1)。

文献[4]中指出：20世纪90年代至2004年，玛多南部草地大范围属于中、轻度持续退化区，玛多县北部草地大范围中度持续退化沙化区。玛多县1987年与1997年两次草地调查统计，20世纪80年代全县草地退化面积为10 702 km2，占全区草地面积46.55%；90年代退化草地面积达16 082 km2占全区草地面积69.95%，其中，重度退化草地面积达9 200 km2呈上升趋势[28]。

将本研究所得2008年玛多县草地退化信息与该县20世纪90年代和80年代两次草地资源调查结果进行对比分析。2008年重度退化草地面积较1997年减少了5 944 km2，中度退化面积增加了3 965 km2，轻度退化草地面积增加了2 933 km2(表2)。可以发现，1997年到2008年，重度退化草地向轻度、重度退化转变，草地退化现象得到缓解。青海省政府于2003年开始，在玛多县划分自然保护区并实行生态移民工程，消除草地超载对草地的危害。玛多县境内共有2个保护区：约古宗列保护区、扎陵湖和鄂陵湖保护区。可见，青海省采取的一系列政策措施很有成效。

3.2草地退化与草地类型 玛多县主要草地类型包括高寒草原类、高寒草甸类(高寒草甸亚类和沼泽类高寒草甸类)沼泽类草地。其中，高寒草原类草地为9 200.66 km2，占草地总面积的48.1%；高寒草甸类草地面积为8 873.56 km2，占草地总面积的46.4%，其中高寒草甸亚类为7 826.41 km2，沼泽化高寒草甸类1 047.15 km2；沼泽类草地所占面积很小，只有48.59 km2。利用不同草地类型对草地退化等级进行分类统计发现(图4)：草地退化等级的分布状况与草地类型存在一定的关系。高寒草原类和高寒草甸类草地均以中度退化为主，沼泽类草地则以轻度退化为主，同时中度退化比例也很高。高寒草原类草地退化较其他类型草地严重，重度退化面积为2 852.03 km2，占高寒草原类草地总面积的31%，而高寒草甸类和沼泽类草地重度退化比例分别为2.7%和0.6%。

图3 2008年玛多草地退化等级分布图

表1 玛多县2008年草地退化面积统计

表2 玛多县3个时期草地退化面积统计

图4 不同草地类型退化等级示意图

高寒草甸类草地包括高寒草甸亚类和沼泽化高寒草甸亚类，这两种亚类高寒草甸草地的退化情况也存在差异。其中高寒草甸亚类以中度退化为主，占该类草地面积的54.8%，而沼泽类高寒草甸亚类以轻度退化为主，占该类草地面积的46.8%(图5)。

4 结论

本研究利用玛多县2008年夏季Landsat TM遥感影像作为数据源，并辅以研究区草地类型图、土地覆被类型图、地面气象观测数据，利用CASA模型反演草地NPP并结合草地类型图计算潜在NPP。以草地NPP减少的相对百分数作为监测草地退化的指标，可操作性强，数据类型要求简单且易获取。利用ArcMap软件平台，叠加草地退化等级分布图与草地类型图，进行对比分析，发现高寒草原类草地退化较其它类型草地严重，重度退化草地达到该类草地面积的31%。针对高寒草原类草地的保护措施有待加强。

图5 两种高寒草甸亚类退化面积示意图

利用该方法所提取草地退化信息与其他学者所研究结果进行对比，结果可信。与玛多县1987年与1997年两次草地调查统计结果进行对比发现：玛多县草地退化现象仍然比较严重。从20世纪80年代至90年代，未退化草地面积急剧减少，中度、重度退化草地面积大幅增加，退化呈现加速的趋势；从20世纪90年代至2008年，重度退化草地面积大幅减少，并向轻度、中度退化转变，草地退化情况得到一定缓解。可能的原因是青海省政府于2003年开始，在玛多县划分自然保护区并实行的生态移民工程，减轻了畜牧超载对草地的危害。同时玛多县保护区的建设也起到一定的积极作用。

[1] 李博.中国北方草地退化及其防治对策[J].中国农业科学,1997,30(6):1-9.

[2] 赵新全,周华坤.三江源区生态环境退化恢复治理及其可持续发展[J].中国科学院院刊,2005，20(6):471-476.

[3] 催庆虎,蒋志刚,刘季科,等.青藏高原草地退化原因评述[J].草业科学,2007,24(5):20-26.

[4] 王海英,周长进,董锁成.三江源地区主要生态环境问题与对策[J].自然资源学报,2002,17(6):713-720.

[5] 李昀,俞联平,陈兴荣,等.基于3S技术的景泰县草地资源现状评价[J].草业科学,2010,27(4):64-68.

[6] 李文娟,马轩龙,陈全功.青海省海东、海北地区草地资源产量与草畜平衡现状研究[J].草业学报,2009,18(5):270-275.

[7] 刘纪远,徐新良,邵全琴.近30年来青海三江源地区草地退化的主要时空特征[J].地理学报,2008，63(4):364-376.

[8] 李辉霞,刘淑.基于NDVI的西藏自治区草地退化评价模型[J].山地学报,2003,21(增刊):69-71.

[9] 高清竹,李玉娥,林而达，等.藏北地区草地退化的时空分布特征[J].地理学报,2005,60(6):965-973.

[10] 王根绪，程国栋.江河源区的草地资源特征与草地生态变[J].中国沙漠,2001,21(2):101-107.

[11] Potter C S,Randerson J T,Fiell C B,etal.Terestrial ecosystem production:a process model based on global satellite and surface data[J].Global Biogeochemical Cycles,1993,7:811-841.

[12] 封建民,王涛,齐善忠,等.黄河源区土地沙漠化的动态变化及成因分析[J].水土保持学报,2004,18(3):141-145.

[13] 李春秀,孙海松.青海省西藏嵩草沼泽化草甸资源及保护[J].草业科学,2009,26(9):77-83.

[14] 农业部畜牧兽医总站主编.中国草地资源[M].北京:中国科学技术出版社,1996.

[15] 玛多县志编辑委员会.玛多县志[M].中国县镇年鉴出版社,2001.

[16] 张静,李希来,谢得雄.三江源地区不同退化草地聚类分析[J].草业科学,2008,25(6):8-13.

[17] 中华人民共和国国家质量监督检验检疫总局. 中华人民共和国国家标准GB19377:天然草地退化、沙化、盐泽化的分级指标[P].2003.

[18] 姜立鹏,覃志豪,谢雯.基于单时相MODIS数据的草地退化遥感监测研究[J].中国草地学报,2007(7):39-43.

[19] 朴世龙,方精云,郭庆华.利用CASA模型估算我国植被净第一性生产力[J].植物生态学报,2001,25(5):603-608.

[20] 陈利军,刘高焕,冯险峰.遥感在植被净第一性生产力研究中的应用[J].生态学杂志，2002,21(2):53-57.

[21] 朴世龙,方精云.1982-1999年青藏高原植被净第一性生产力及其时空变化[J].自然资源学报，2002,17(3):373-380.

[22] Sellers P J,Tucker C J,Collatz G J.A golbal 1°by 1°NDVI data set for climate studies.Part 2:the generation of global fields of terrestrial biophysical parameters from the NDVI[J].International Journal Remote Sensing,1994,15(17):3519-3545.

[23] 称德瑜,高安,古建泉.农业气候学[M].第1版.北京:气象出版社,1994:114-115.

[24] Chen L J,Liu G H,Feng X F.Estimation of net primary productivity of terrestrial vegetation in china by remote sensing[J].Acta Botanica Sinica,2001,43(11):1191-1198.

[25] 张娜,于贵瑞,于振良.基于3S的自然植被光能利用率的时空分布特征的模拟[J].植物生态学报,2003,27(3):25-36.

[26] Runing S W,Thornton P E,Nemani R.Global terrestrial gross and net primary productivity from the earth observing system[A].In:Sala Jackson O R,Mooney H.Methods in Eco-System science[C].New York:Springer Verlag,2000:44-57.

[27] 杨景梅,邱金桓.用地面湿度参量计算我国整层大气可降水量及有效水汽含量方法的研究[J].大气科学,2002,26(1):9-22.

[28] 周华坤,周立,刘伟,赵新全,来德珍.青海省玛多县草地退化原因及畜牧业可持续发展[J].中国草地,2003,25(6):63-67.