黄河源区高寒草地NDVI格局与梯度变化

赵志平，吴晓莆，李 果，李俊生

(中国环境科学研究院，北京100012)

黄河源区位于青藏高原东部，是黄河流域上游主要产流区、水源涵养区，也是我国重要的生态屏障[1]。由于严酷而独特的自然环境，该地区同时也是我国生态环境脆弱区。近年来该区因高寒草甸退化[2-3]、高寒草原沙化[4-5]、草原鼠害猖獗[6]、沼泽和湿地面积减少[7]以及生物多样性急剧萎缩[8]而成为国内学者研究的热点地区[9-11]。受水、热梯度影响，该区顺黄河源头而下分布着高寒草甸和高寒草原两种草地类型。国内学者在内蒙古典型温带草原的研究证实，草地生长状况与水、热梯度有着密切的关系[12-14]，但在黄河源区国内尚未见相关报道。

归一化植被指数(NDVI)是反映地表植被覆盖状况、生长状况以及生物量的重要参数[15-19]。本研究利用2000-2011年黄河源区NDVI作为草地生长状况的代用指标，分析该区草地状况随水热梯度变化过程，揭示气候环境因子的梯度变化对高寒草原和高寒草甸状况影响程度，为全球气候变化与草地生态系统响应的研究提供科学依据。

1 研究区概况

研究区位于青藏高原东部、巴颜喀拉山北部的黄河流域源头地区，地处32.83°-35.64° N，95.87°-101.79° E。行政区划上该区主体位于果洛藏族自治州，主要涉及曲麻莱、玛多、称多、达日、甘德、玛沁、久治和班玛8县，总面积62 224 km2。该区平均海拔4 393 m，气候上属高原高寒气候，表现为冷热两季交替、干旱两季分明，年平均气温-3.64 ℃，年降水量513 mm，从东南向西北年平均气温和年降水量逐渐降低，海拔逐渐升高[20]。草地是该区的主要植被类型，主要为高寒草原和高寒草甸。土壤类型主要有高山草甸土、高山草原土、灰褐土、栗钙土、沼泽土和风沙土等类型，其中以高山草甸土分布最多。

2 数据和方法

2.1数据 本研究用到的NDVI数据来源于美国MODIS Terra 2000-2011年植被指数产品，空间分辨率为1 km(https://wist.echo.nasa.gov/～wist/api/imswelcome/)。1 km分辨率数字高程模型数据和1∶100万植被类型数据来源于国家自然科学基金委员会“中国西部环境与生态科学数据中心”(http://westdc.westgis.ac.cn)。研究区草地空间分布从1∶100万植被类型数据中提取得到(图1)。气象站点观测数据由中国气象局数据共享中心提供，包括1980-2009年黄河源区及周边气象站点日观测数据，数据项为日平均气温、日最高气温、日最低气温、风速、相对湿度、降水和日照时数。

图1 研究区草地类型空间分布Fig.1 Spatial grassland pattern of research area

2.2方法 研究区为黄河源头地区，流域边界提取是基于1 km分辨率数字高程模型数据，利用ArcGIS软件中水文分析模块，进行集水流域提取而生成的。由于研究区位于青藏高原高寒区，植被生长期较短，因此采用最大合成法(Maximum Value Composite Syntheses，MVC)获得每个像元一年中草地NDVI最大值(NDVImax，记为NDVI)来代表当年植被生长状况：

NDVImax(x,t)=max[NDVI(x,t,i)]

(1)

式中，x表示空间位置，t表示年份，i表示t年中第i个16 d，其范围在1～23。然后采用最小二乘法线性回归方程的斜率来分析2000-2011年研究区NDVI变化趋势，计算公式为：

(2)

式中，NDVI为一年中NDVI最大值，slope为2000-2011年NDVI变化斜率；j为年序号，n=12。

NDVI变异系数(CV)计算公式为：

(3)

Thornthwaite[21]使用湿润指数来指示气候的湿润程度，并提出了以下计算湿润指数(Im)的公式：

(4)

式中，P为年降水量，ET0为潜在蒸散。本研究采用此式计算湿润指数来定量表示研究区湿润程度。潜在蒸散(ET0)采用联合国粮农组织(FAO)1998年对Penman-Monteith模型修订后的版本计算[22-23]：

(5)

式中，Rn为净辐射，G为土壤通量，γ为干湿常数，Δ为饱和水汽压曲线斜率，U2为2 m高处的风速，ea为实际水汽压，es为平均饱和水汽压。净辐射Rn的计算公式如下：

(6)

式中，σ为Stefan-Boltzmann常数(4.903×10-9MJK-4·m-2·d-1)，Tmax,k、Tmin,k分别为绝对温标的最高和最低气温，n为实际日照时数，N为可照时数，RSO为晴天辐射。

利用ANUSPLIN软件将黄河源区及周边气象站点数据插值形成本区空间1 km栅格气候数据[24-25]，包括年平均气温、年降水量和湿润指数[26]，其中湿润指数由公式(4)、(5)和(6)计算得到。基于上述工作，形成研究区年平均气温、降水量、湿润指数和海拔梯度数据。

3 结果

3.1环境梯度变化 利用1980-2009年年平均气温、降水量观测数据和海拔数据，分析该区环境因子梯度变化。近30年来研究区年平均气温在-7～1 ℃，年降水量在334～885 mm，湿润指数在-37～67，海拔在3 515～5 175 m。自西北向东南年平均气温、年降水量和湿润指数呈梯度上升，海拔高度呈梯度下降(图2)。

3.2NDVI梯度变化 研究区近12年来草地NDVI均值变化在0.062 7～0.869 8，平均值为0.557 8。空间上该区从东南至西北沿黄河逆流而上NDVI逐渐下降，该格局明显是受水、热条件制约而形成。NDVI最高的地区为久治和甘德县，其次为达日和玛沁县，玛多、称多和曲麻莱县NDVI最低。扎陵湖和鄂陵湖以北的高寒草原地区为研究区NDVI低值区(图3)。

研究区草地NDVI与年平均气温具有斜率为0.052 7/℃的极显著正相关关系(P<0.01)，表明该区气温对植被生长具有正向作用(图4)。草地NDVI与年降水量具有极显著相关关系(P<0.01)，总体上该区降水量对植被生长具有正向作用，但降水量超过750 mm后，NDVI反而下降。草地NDVI与湿润指数具有极显著相关关系(P<0.01)，总体上该区湿润程度对植被生长具有正向作用，但湿润指数超过40后，NDVI反而下降。草地NDVI与海拔高度具有-0.03·100 m-1的极显著(P<0.01)负相关关系，表明该区海拔高度对植被生长具有负向作用。

图2 环境因子空间格局Fig.2 Spatial pattern of environmental factors

图3 2000-2011年研究区草地NDVI均值空间分布Fig.3 Spatial pattern of NDVI in research area from 2000 to 2011

图4 草地NDVI均值随年平均气温、年降水量、湿润指数和海拔的梯度变化Fig.4 Change of grassland NDVI with annual average temperature, annual precipitation, moisture index and elevation

3.3NDVI变化率 研究区近12年来草地NDVI年变化速率介于-0.049 7～0.031 5，平均值为0.002 4。空间上该区从东南至西北NDVI变化速率上升。NDVI变化速率最高的地区为玛多、曲麻莱和称多县，其次为达日和玛沁县，久治和甘德县NDVI变化速率最低。扎陵湖和鄂陵湖附近地区为NDVI变化速率较高，达日县黄河以南局部地区NDVI变化速率较低(图5)。

NDVI变化率能够反映草地NDVI变化的方向和速率。研究区草地NDVI变化率与年平均气温具有极显著相关关系(P<0.01)，低温地区草地NDVI上升速率较大；与年降水量具有极显著相关关系(P<0.01)，总体上年降水量少的地区草地NDVI上升速率较大；与湿润指数具有极显著相关关系(P<0.01)，总体上湿润程度低的地区草地NDVI上升速率较大；与海拔高度具有极显著相关关系(P<0.01)，总体上高海拔地区草地NDVI上升速率较大(图6)。

3.4NDVI变异系数 研究区近12年来草地NDVI变异系数变化范围为0.008 7～1.094 7，平均值为0.076 7，表明草地NDVI变化幅度较小。空间上该区从东南至西北NDVI变异系数逐渐上升，与NDVI变化率相似。NDVI变异系数最高的地区为玛多、曲麻莱和称多县，其次为玛沁和达日县，甘德和久治县NDVI变异系数最低。扎陵湖和鄂陵湖周围地区为研究区NDVI变异系数高值区，久治县为NDVI变异系数低值区(图7)。

变异系数能够反映草地NDVI变动幅度。研究区草地NDVI变异系数与年平均气温具有-0.01/℃的极显著负相关关系(P<0.01)，即高温地区草地NDVI变异系数较小；与年降水量具有极显著相关关系(P<0.01)，总体上年降水量多的地区草地NDVI变异系数较小；与湿润指数具有极显著相关关系(P<0.01)，总体上湿润程度高的地区草地NDVI变异系数较小；与海拔具有极显著正相关关系(P<0.01)，即高海拔地区草地NDVI变异系数较大(图8)。

图5 研究区草地NDVI变化速率空间分布Fig.5 Spatial pattern of NDVI slope in research area

图6 草地NDVI变化率随年均温、年降水量、湿润指数和海拔的变化Fig.6 Variation of grassland NDVI with annual average temperature, annual precipitation, moisture index and elevation

图7 研究区草地NDVI变异系数空间分布Fig.7 Spatial pattern of NDVI variable coefficient in research area

图8 草地NDVI变异系数随年均温、年降水量、湿润指数和海拔的变化Fig.8 Variation of the coefficient of variation of grassland NDVI with annual average temperature, annual precipitation, moisture index and elevation

4 讨论

本研究表明，研究区草地NDVI分布表现与年平均气温、年降水量和湿润指数之间均呈正相关关系，与海拔呈负相关关系，这是由该区自然环境现状决定的。

草地NDVI变化率格局为从东南至西北上升，一方面可能受到水、热条件变好的影响，另外一方面也与研究区自2005年以来实施的退牧还草、以草定畜和草地退化治理等生态恢复措施有关。近12年来，研究区年平均气温变化率为0.089 ℃，其中中部和西部增温速率较大(图9)；年降水量变化率为11.444 mm(图9)，其中东部和西部降水量增加速率较大；湿润指数年变化率为1.523，其中西北部地区湿润指数增加速率较大，东南部地区湿润指数增加速率较小(图9)。2005年三江源生态保护和建设工程实施以来，果洛藏族自治州年末牲畜存栏数持续下降(图9)，尤其是研究区西部的玛多县，在生态移民与减畜、退化草地治理和人工降水等方面的工作使得该区生态状况明显好转。因此，近12年来研究区草地NDVI变化率格局是气候变化、草地现实载畜量明显下降和生态建设工程实施三方面因素叠加的结果。

草地NDVI变异系数是近12年NDVI标准差与均值的比值。本研究中，草地NDVI变化与其分布格局叠加，导致草地NDVI变异系数从东南至西北逐渐上升。

5 结论

近12年来黄河源区草地NDVI平均值为0.557 8，从东南至西北NDVI逐渐下降，年变化速率均值为0.002 4，变异系数均值为0.076 7，变化幅度较小，但是NDVI变化速率和变异系数上升。

黄河源区气候环境因子自西北向东南存在梯度变化，形成水热梯度，从而导致草地NDVI、NDVI变化率和NDVI变异系数存在梯度变化。年平均气温与NDVI极显著相关，与NDVI变化率和NDVI变异系数极显著负相关。年降水量和湿润指数与NDVI极显著相关，与NDVI变化率和NDVI变异系数极显著负相关。海拔与NDVI具有极显著负相关关系，与NDVI变化率和NDVI变异系数极显著相关。

图9 研究区气候因子变化率和果洛藏族自治州年末牲畜存栏数变化Fig.9 Variation of climate factors and year-end livestock stock in Golog Tibetan Autonomous Prefecture

[1] 赵成章,贾亮红.黄河源区退牧还草工程生态绩效与问题[J].兰州大学学报(自然科学版),2009,45(1):37-41.

[2] 周华坤,赵新全,周立,等.青藏高原高寒草甸的植被退化与土壤退化特征研究[J].草业学报,2005,14(3):31-40.

[3] 马玉寿,郎百宁.建立草业系统恢复青藏高原“黑土型”退化草地[J].草业科学,1998,15(1):5-9.

[4] 曾永年,冯兆东.黄河源区土地沙漠化时空变化遥感分析[J].地理学报,2007,62(5):529-536.

[5] 郄妍飞,颜长珍,宋翔,等.近30a黄河源地区荒漠遥感动态监测[J].中国沙漠,2008,28(3):405-409.

[6] 周雪荣,郭正刚,郭兴华.高原鼠兔和高原鼢鼠在高寒草甸中的作用[J].草业科学,2010,27(5):38-44.

[7] 王根绪,丁永建,王建,等.近15年来长江黄河源区的土地覆被变化[J].地理学报,2004,59(2):163-173.

[8] 杨建平,丁永建,陈仁升.长江黄河源区高寒植被变化的NDVI记录[J].地理学报,2005,60(3):467-478.

[9] 谷源泽,李庆金,杨风栋,等.黄河源地区水文水资源及生态环境变化研究[J].海洋湖沼通报,2002(1):18-25.

[10] 薛娴,郭坚,张芳,等.高寒草甸地区沙漠化发展过程及成因分析——以黄河源区玛多县为例[J].中国沙漠,2007,27(5):725-732

[11] 芦清水,赵志平.应对草地退化的生态移民政策及牧户响应分析——基于黄河源区玛多县的牧户调查[J].地理研究,2009,28(1):143-152.

[12] 李林芝,张德罡,辛晓平,等.呼伦贝尔草甸草原不同土壤水分梯度下羊草的光合特性[J].生态学报,2009,29(10):5271-5279.

[13] 程杰,呼天明,程积民.黄土高原白羊草种群分布格局对水热梯度的响应[J].草地学报,2010,18(2):167-171.

[14] 白永飞,张丽霞,张焱,等.内蒙古锡林河流域草原群落植物功能群组成沿水热梯度变化的样带研究[J].植物生态学报,2002,26(3):308-316.

[15] 王维,王文杰,李俊生,等.基于归一化差值植被指数的极端干旱气象对西南地区生态系统影响遥感分析[J].环境科学研究,2010,23(12):1447-1455.

[16] Hua W,Fan G Z,Zhou D W,etal.Preliminary analysis on the relationships between Tibetan Plateau NDVI change and its surface heat source and precipitation of China [J].Science in China Series D:Earth Sciences,2008,51(5):677-685.

[17] Han L J,Wang P X,Yang H,etal.Study on NDVI-Ts space by combining LAI and evapotranspiration[J].Science in China Series D:Earth Sciences,2006,49(7):747-754.

[18] 王蕊,李虎.2001-2010年蒙古国MODIS-NDVI时空变化监测分析[J].地球信息科学学报,2011,13(5):665-671.

[19] 张宏斌,杨桂霞,吴文斌,等.呼伦贝尔草原MODIS NDVI的时空变化特征[J].应用生态学报,2009,20(11):2743-2749.

[20] 青海果洛自治州人民政府网站[EB/OL].(2011-05-24)[2012-12-20].http://www.guoluo.gov.cn.

[21] Thornthwaite C W.An approach toward a rational classification of climate[J].Geographical Review,1948,38(1):55-94.

[22] 王懿贤.高度对彭曼蒸发公式二因子δ/(δ+γ)与γ/(δ+γ)的影响[J].气象学报,1981,39(4):503-506.

[23] Allen R G,Pereira L S,Raes D,etal.Crop evapotranspiration: guidelines for computing crop water requirements[A].Crop Evapotranspiration-Guidelines for Computing Crop Water Requirements[C].Rome:Food and Agriculture Organization of the United Nations,1998:377-384.

[24] Hutchinson M F.Anusplin Version 4.37 User Guide[M].Canberra:The Australian National University,2007.

[25] Hutchinson M F.Interpolation of rainfall data with thin plate smoothing splines I two dimensional smoothing of data with short range correlation[J].Geographic Information Decision Analysis,1998,2:153-167.

[26] 王英,曹明奎,陶波,等.全球气候变化背景下中国降水量空间格局的变化特征[J].地理研究,2006,25(6):1031-1040.