黄土高原是我国乃至世界上水土流失最严重的地区。黄河上中游地区人类活动对自然植被的破坏是造成该生态问题的主要原因，尤其是对坡地植被的破坏，大大增加水土流失和洪涝灾害的风险^[1-4]。为了解决这一生态环境问题，中国政府于1999年启动“退耕还林还草”重点生态修复工程^[5]。目前大量研究已表明：近10年内，人类活动对黄土高原下垫面影响巨大，植被覆盖明显增加^[6-8]；同时，对黄土高原的产流产沙也产生了显著的影响^[9-10]。但是，也有相关研究指出，该生态恢复工程的效益被夸大，且存在一定的负影响^[11-13]。因此，对近十几年黄土高原的植被变化进行时空变化分析，有助于全面了解该区域植被变化面积与变化程度，并对科学评价流域生态环境和制定有效的水土保持措施，均具有重要意义。

大量研究表明，植被的结构特征是影响植被水土保持效益的重要因素^[14-15]。叶面积指数 (Leaf Area Index, LAI) 是最重要的植被结构参数之一，也是陆地生态、水文模型模拟中不可或缺的植被参数。LAI是指单位土地面积上，植物叶片总面积占土地面积的倍数，该参数能够同时反映植被水平与垂直结构特征。目前，遥感技术已经成为监测区域尺度植被动态变化的最有效工具。NDVI数据是进行区域植被监测的常用数据源，如AVHRR^[16-17]、SPOT^[18]和MODIS^{[7, 19]}；但NDVI不能直接代表特定植被结构参数，且研究中所选遥感数据空间尺度较大，难以反映地表植被的真实变化情况。目前，针对黄土高原植被LAI的时空变化分析研究还相对较少，笔者首次利用NASA最新发布的高分辨率LAI产品，从区域到像素尺度，对黄土高原植被LAI的时空变化特征进行提取，从不同角度，全面了解2000—2014年黄土高原植被LAI的变化情况，同时对植被LAI的变化原因进行分析。

1 研究区概况

黄河是中国第二长河，是世界第五长河流，也是世界上含沙量最多的河流。黄土高原 (E 100.9°~114.5°，N 33.7°~41.3°) 是世界上水土流失最严重和生态环境最脆弱的地区之一，也是黄河泥沙主要来源区，黄土高原海拔一般在1 000~1 300 m，地貌起伏不平，坡陡沟深，千沟万壑，沟谷面积占40%~50%。黄土高原主要产流产沙的8个子流域分别为皇甫川、窟野河、无定河、延河、北洛河、泾河、渭河和汾河，其地理位置见图 1。

2 研究方法 2.1 数据获取与处理 2.1.1 遥感数据

使用2015年NASA最新发布的MODIS LAI产品 (MOD15A2H)，对2000—2014年黄土高原植被LAI进行时空变化分析研究。其产品是基于地表反射率数据和土地覆盖类型数据，综合利用经验模型和三维辐射传输模型查找表法，对植被真实LAI进行反演得到^[20]，其空间分辨率500 m，时间分辨率8 d。其他数据包括MOD44B产品、ASTER-DEM和MCD12Q1产品。

2.1.2 遥感数据处理

利用MODIS数据处理软件MRT (MODIS Reprojection Tools)，对MODIS产品进行数据预处理。MOD44B产品表示年内植被最大覆盖度，将250 m分辨率的植被覆盖度，上推到500 m分辨率，尺度上推方法采用聚合平均法。对于DEM数据处理，采用基于30 m分辨率的ASTER-DEM，计算得到局地坡度和坡向。

对于MOD15A2H产品，由于在实际成像过程中，受到云与气溶胶的影响，导致遥感影像中存在不可避免的噪声信号，从而不能客观的反映地表真实LAI情况。针对这一问题，笔者采用Chen Jin等^[21]提出的迭代滤波方法，对LAI数据进行滤波处理，某一典型植树造林区域像素的LAI时间序列重构前后对比情况见图 2。

从图 2可知，由于受到云的影响，2003和2008年夏季LAI数据上下波动剧烈。这显然不符合植被的生长规律，而重构算法对受云影响的LAI数据进行平滑，使得重构数据曲线更符合植被物候变化规律，从而可以大大提高数据分析精度。

按照上述方法对MODIS产品处理后，得到2000—2014年时间分辨率为8 d的LAI数据序列。为了反映植被的年内最佳生长状况，对每年的LAI数据序列中最大值进行提取，然后构建2000—2014年LAI年最大值数据序列。

2.2 LAI时空变化特征

不同区域尺度植被LAI变化发生的空间位置、变化范围和变化强度都是植被变化监测研究的重要内容。

2.2.1 变化面积比例

采用Mann-Kendall (M-K) 非参数检验方法^[22]，对黄河流域植被LAI时间序列变化趋势进行统计分析。若时间序列M-K检验显著性水平α < 0.05，则认为该趋势变化为显著性变化；若显著性水平>0.05，则认为该趋势没有统计学意义，即表示没有发生显著性变化。

LAI变化面积比例 (D) 是指某个区域范围内，LAI显著变化 (P < 0.05) 区域占该区域范围总面积的比重，该指标可反映LAI在某个区域内的变化剧烈程度。计算公式如下：

$D = \frac{a}{s} \times 100\% 。$

(1)

式中：A为LAI显著变化的区域面积，km²；S为某区域范围的总面积，km²。

2.2.2 绝对与相对变化率

采用Kendall倾斜率β表示每年LAI绝对变化率，该指标可反映LAI的年绝对变化幅度，β>0表示上升趋势，β < 0表示下降趋势，其变化趋势显著性水平参考M-K检验结果。计算公式如下：

$\beta = {\rm{Median}}\left( {\frac{{{x_j}-{x_i}}}{{j-i}}} \right), \left( {1-\le i < j \le n} \right)。$

(2)

式中:x_j和x_i分别为第j年和第i年的观测数值，j>i；n为年总数。

R为每年LAI相对变化率，该指标可反映LAI的年相对变化程度，其计算公式如下：

$R = \frac{\beta }{{\frac{1}{n}\sum\limits_{i = 1}^n {{x_i}} }} \times 100\% 。$

(3)

式中:β为LAI绝对变化率; x_i为第i年的观测数值；n为年总数。由于该值利用年均值对β进行标准化，从而使得不同空间像素的趋势变化程度具有可比性^[23]。

2.2.3 区域/流域尺度LAI时空变化分析

为了反映黄土高原以及其内部不同流域间LAI的变化差异，对黄土高原以及其内部流域面积较大，且产流产沙量较大的8个子流域的植被LAI进行时空变化分析。采用绝对变化率和相对变化率2个指标，对区域/流域尺度LAI变化特征进行评估。

2.2.4 像素尺度LAI时空变化分析

为了反映植被LAI变化特征在黄土高原、不同子流域、不同坡度和坡向以及不同植被覆盖类型的差异性，采用变化面积比和相对变化率2个指标，对不同区域范围内的植被LAI变化特征差异进行评估。其中：子流域为所选8个典型流域；坡度分类执行水利部颁布的SL 190—2007《土壤侵蚀分类分级标准》^[24]，将坡度划分为6个等级： < 5°，5°~8°，8°~15°，15°~25°，25°~35°，>35°；坡向分类按照45°夹角，由北开始，顺时针方向旋转分为9个坡向：平地 (-1°)、北坡 (N，0~22. 5°，337.5°~360°)、东北坡 (NE，22.5°~67.5°)、东 (E，67.5°~112.5°)、东南坡 (SE，112.5°~157.5°)、南坡 (S，157.5°~202.5°)、西南坡 (SW，202.5°~247.5°)、西坡 (W，247.5°~292.5°)、西北坡 (NW，292.5°~337.5°)^[25]；植被覆盖类型采用2001年的MOD12Q1土地覆盖产品，提取其中的农田、草地、灌木以及森林类型进行分析。

2.3 LAI时空变化原因

退耕还林还草政策，对黄土高原植被LAI的改善起到至关重要的作用，但不同区域植被LAI变化原因的空间分布规律并没有被探讨。由于在物理意义上，LAI同时反映植被水平覆盖度与垂直维结构，为了分析不同维植被结构特征对植被LAI的影响，首先将植被垂直维结构特征定义为垂直维LAI，指单位植被覆盖面积上，植物叶片总面积占植被覆盖面积的倍数。其具体公式如下：

$\frac{{{\rm{LAI}}}}{{{\rm{VFC}}}} = \frac{{{S_{{\rm{land}}}}{\rm{LAI}}}}{{{S_{{\rm{veg}}}}}} = \frac{{{S_{{\rm{land}}}}\frac{{{S_{{\rm{leaf}}}}}}{{{S_{{\rm{land}}}}}}}}{{{S_{{\rm{veg}}}}}} = \frac{{{S_{{\rm{leaf}}}}}}{{{S_{{\rm{veg}}}}}} = {\rm{LA}}{{\rm{I}}_{{\rm{ver}}}}。$

(4)

式中：LAI为植被叶面积指数；VFC为植被总覆盖度；LAI_ver为垂直维LAI; S_leaf、S_veg和S_land分别为植物叶片总面积、植被覆盖面积和土地面积。

利用MOD44B产品，提取树木覆盖度与草本植被覆盖度，二者相加即为植被总覆盖度，并依据公式 (4)，计算得到垂直维LAI。首先将区域尺度LAI与4种植被参数 (植被总覆盖度、垂直维LAI、树木覆盖度、草本植被覆盖度) 分别进行线性回归分析。为了更进一步分析植被类型变化对LAI的影响，根据公式 (5) 和 (6)，分别计算像素尺度LAI与树木覆盖度和草本植被覆盖度的偏相关系数，并对偏相关系数进行显著性检验。通过偏相关分析空间分布图和LAI趋势变化空间分布图的叠加分析，同时确定每个像素LAI的变化趋势和与LAI变化显著相关的变量，从而得到植被LAI显著变化区域的变化原因空间分布图。

${r_{12, 3}} = \frac{{{r_{12}}-{r_{13}}{r_{23}}}}{{\sqrt {1-r_{13}^2} \sqrt {1-r_{23}^2} }}, $

(5)

$t = \frac{{r\sqrt {n-q-2} }}{{\sqrt {1-{r^2}} }}。$

(6)

式中:r_{12, 3}为偏相关系数; r₁₂、r₁₃、r₂₃分别为不同要素之间的Pearson相关系数；n为样本数，q为阶数；统计量t服从n-q-2个自由度的t分布。若t检验显著性水平α < 0.05，则认为2个要素之间存在显著的相关性。

3 结果与分析 3.1 区域/流域尺度LAI时空变化特征

2000—2014年，整个黄土高原以及8个典型子流域的植被LAI年最大值变化趋势和变化特征统计结果见图 3和表 1，整个黄土高原及其子流域的植被LAI均呈显著上升趋势 (P < 0.01)。整个黄土高原地区，植被LAI的年绝对变化幅度为0.042，年相对变化程度为2.71%。在8个子流域中，汾河流域的LAI年绝对变化幅度最大，达到0.081，从大到小排列顺序依次为汾河>延河>北洛河>泾河>渭河>无定河>皇甫川>窟野河。延河流域的LAI年相对变化程度最大，达到5.37%，从大到小排列顺序依次为延河>皇甫川>无定河>汾河>窟野河>泾河>北洛河>渭河。

3.2 像素尺度LAI时空变化特征 3.2.1 LAI变化总体情况

利用公式 (2) 和 (3)，基于像素尺度计算每个像素的绝对变化率和相对变化率，计算结果见图 4。结果表明：绝对变化率呈现出由西北向东南方向递增的趋势，这与黄土高原气候湿润条件的空间梯度分布趋势较为一致，说明在湿润地区，植被由于能够获取充足的水分，导致其生长较快，LAI绝对变化幅度也较大；相反，干旱地区植被由于缺水，导致其生长较慢，LAI绝对变化幅度也较小。因此，LAI绝对变化率在一定程度上反映了地区湿润条件。反观LAI相对变化率，消除了地区湿润条件对植被变化的影响，其空间分布并没有随着地区湿润条件的变化而变化，可真实表达LAI变化程度，且不同区域的LAI变化程度可直接进行比较。

利用公式 (1) 对整个黄土高原的变化面积比进行统计。结果表明：58.6%的区域，LAI呈现显著增加趋势，增加程度最大的区域，主要集中在河口-龙门区间；仅有0.9%的区域，LAI呈现显著减少趋势，主要分布在渭河、伊洛河和湟水流域零星区域；其余40.5%的区域，LAI保持稳定状态。LAI变化趋势空间分布见图 4。

3.2.2 LAI变化的流域差异

LAI变化的流域差异如图 5所示。结果表明：根据变化面积比和相对变化率的排序比较，延河流域LAI变化剧烈程度最大，渭河流域LAI变化剧烈程度最小。除渭河流域外，其余流域内变化面积比均超过69%。

3.2.3 LAI变化的坡度差异

将LAI重采样到30 m尺度，与30 m尺度坡度数据进行叠加分析。LAI变化的坡度差异见图 6。结果表明：15°~25°和25°~35° 2个坡度的植被LAI变化程度最为剧烈，变化面积比均超过60%，年相对变化率均超过3%；>35°区域的植被LAI变化程度最小，可能是由于坡度较大，导致水分较易流失，植被生长较为缓慢。

3.2.4 LAI变化的坡向差异

将LAI重采样到30 m尺度，与30 m尺度坡度数据进行叠加分析。LAI变化的坡向差异见图 7。结果表明，不同坡向植被LAI的变化面积比和相对变化率基本一致。说明黄土高原LAI的变化，在坡向上没有显著差异，在一定程度上反映出植被变化并没有明显的坡向规律。

3.2.5 LAI变化的植被覆盖差异

LAI变化的植被覆盖差异见图 8。结果表明，不同植被类型的LAI变化剧烈程度从大到小排序为农田>草地>灌木>森林。这也反映出退耕还林还草政策的实施效果，使得原来的农田和草地区域植被LAI变化剧烈。

3.3 LAI时空变化原因分析

区域尺度LAI与4种植被参数的线性回归分析结果见图 9。从图中可知，对于整个黄土高原植被总覆盖度与植被LAI相关性较低 (R²=0.21)，垂直维LAI与植被LAI之间存在显著的正相关性 (R²=0.94)，树木覆盖度与植被LAI也存在显著正相关性 (R²=0.82)，说明黄土高原植被LAI的变化与植被垂直维结构的变化更为相关，而树木覆盖度的持续增加，则是植被垂直维结构变化的重要原因之一。

通过偏相关分析可知，某地区的植被LAI变化主要是由木本植被或草本植被的改变所致。将像素尺度LAI与树木覆盖度和草本植被覆盖度偏相关分析空间分布图，与LAI趋势变化空间分布图进行叠加分析见图 10。结果表明，在LAI显著增加区域，42%的区域是由树木覆盖度的增加所致，4.3%的区域是由草本植被覆盖度的增加所致，7.5%的区域是由树木和草本植被覆盖度共同增加所致。在LAI显著减少区域中，9.5%的区域是由树木覆盖度的减少所致，16.8%的区域是由草本植被覆盖度减少所致，4.9%的区域是由树木和草本覆盖度共同减少所致。

4 结论

1) 黄土高原植被LAI的年绝对变化幅度为0.042，年相对变化程度为2.71%。在8个子流域中，汾河流域的LAI年绝对变化幅度最大，达到0.081，从大到小排列顺序为汾河>延河>北洛河>泾河>渭河>无定河>皇甫川>窟野河。延河流域的LAI年相对变化程度最大，达到5.37%，从大到小排列顺序为延河>皇甫川>无定河>汾河>窟野河>泾河>北洛河>渭河。

2) 在空间上，黄土高原58.6%的区域LAI呈现显著增加趋势，仅有0.9%的区域LAI呈现显著减少趋势，其余40.5%的区域LAI保持稳定状态。延河流域LAI的变化面积比最大，达到97%，而渭河流域LAI的变化面积比最小，仅为47.5%。15°~25°和25°~35° 2个坡度的植被LAI变化程度最为剧烈，变化面积比均超过60%，年相对变化率均超过3%。黄土高原LAI的变化在坡向上没有显著差异。不同植被类型的LAI变化剧烈程度，从大到小排序为农田>草地>灌木>森林。

3) 对于整个黄土高原与植被总覆盖度相比，植被垂直维结构与植被LAI的变化更为相关；其中，树木覆盖度的持续增加是植被垂直维结构变化的重要原因之一。在LAI显著增加区域中，42%的区域是由树木覆盖度的增加所致，4.3%的区域是由草本植被覆盖度的增加所致，7.5%的区域是由树木和草本植被覆盖度共同增加所致。在LAI显著减少区域中，9.5%的区域是由树木覆盖度的减少所致，16.8%的区域是由草本植被覆盖度减少所致，4.9%的区域是由树木和草本覆盖度共同减少所致。

5 讨论

在已有研究报道中，均使用低分辨率NDVI或植被覆盖度数据对黄土高原植被动态变化进行分析，但NDVI在物理意义上不能直接代表植被结构参数。笔者选用500 m空间分辨率的LAI数据进行分析，与MODIS 1km^[26]、LTDR 5km^[7]、GIMMS 8km^[27]NDVI数据的分析结果相比，能够更为精细的反映黄土高原植被垂直结构的时空变化特征。另外，与已有成果相比，本研究利用相对变化率指标，对黄土高原植被变化剧烈程度的空间异质性规律进行分析和总结。在此需要指出的是，MODIS LAI数据是全球尺度产品，对于黄土高原这一特定区域必然存在绝对误差，但LAI年际相对变化趋势仍具有一定的可信度。今后可利用多尺度和多源遥感数据，对黄土高原植被LAI变化进行多尺度分析研究。