陕西彬长矿区NDVI3g(1982—2013)变化趋势及气候响应

马雯思,马 超,刘玮玮

(1.河南理工大学 测绘与国土信息工程学院,河南 焦作 454000; 2.河南理工大学 矿山空间信息国家测绘与地理信息局重点实验室，河南 焦作 454000; 3.中国科学院西北生态环境资源研究院 冰冻圈科学国家重点实验室，甘肃 兰州 730000; 4.中国科学院大学，北京 100049)

煤炭资源的开发利用造成了地表沉陷、土地挖损、煤矸石堆存占地等无法弥补的损失，进而造成水土流失加剧、耕地退化、土地荒漠化等更为严重的生态环境问题[1]。矿区生态环境的破坏主要体现在对水、土壤、植被大气的影响，绿色植被作为生态环境最敏感和最主要的环境要素，是矿区生态环境恶化的“指标器”。遥感植被指数可以定量分析植被覆盖及其生长状况，其中归一化差值植被指数(NDVI)为最常用的一种。

利用遥感数据动态监测自然地表植被覆盖变化，国内外已有诸多研究，如Myneni,Fang,Tucker,Los等利用卫星遥感数据对全球范围内植被覆盖变化，研究表明:过去20 a北半球高纬度地区陆地植被活动在显著增强[2-7]。方精云等利用AVHRR/NDVI时间序列数据(1982—1999年)研究中国地区NDVI变化，结果显示:全国平均归一化植被指数(NDVI)年增加7.4%，中国西北、陕西地区植被指数增加高于6.5%[8]，NING等运用归一化植被指数及气候相关因素对中国北方黄土高原研究表明:过去几十年里，年均NDVI最大值上升速率达0.09/10 a，荒漠化面积由54.7%(1999年)降至8.9%(2012年)，气温、降水变化是植被生长状况改变的主要诱因[9]。

针对矿区人为影响植被监测的研究相对较少:RASIM等运用两种植被指数(归一化植被指数NDVI和归一化湿润指数NDWI)对加拿大阿萨巴斯卡油砂区进行检测，研究证实:1984—2012年，NDVI平均减少18.6%，归一化水体指数(NDWI)平均减少31.0%，是气候变化和人为因素的干扰结果[10]。PROSPER等利用遥感和地理信息系统技术评估佤邦东部矿区土地覆盖变化:植被出现下降趋势，NDVI像元均值从0.48降至0.11[11]。HUANG等利用时空自适应反射融合模型对山西大同矿区进行NPP和生物量损失估算表明:2001—2010年，矿区年NPP平均减少24.71 gC/m2，年生物量损失率达33.48 gC/m2，气候变化与平均NDVI值呈正相关，采煤对植被的干扰导致生物量的损失及植被吸收能力的减弱[12]。徐占军等利用Landsat5 TM数据对徐州矿区研究表明，采矿活动对NPP 的影响大于气候变化对NPP的影响，NPP变化对采矿活动具有敏感性[13];郝成元等基于EOS/MODIS卫星遥感数据研究分析，研究时段内(2001—2006年)潞安矿区NPP呈显著减少趋势，其时间异质性主要与煤炭开采、农业耕作等类活动联系密切[14];马超等采用SPOT2/4卫星遥感数据研究采矿扰动区NDVI变化趋势，潞安矿区8个试验工作面在开采前NDVI相对稳定，开采后出现下降，2004—2007年期间，年均降幅达14.83%，认为采矿扰动是矿区地表NDVI变化的主要诱因[15]。

鉴于国外对非自然生态区植被遥感监测研究较少，且国内的多数研究时序间断或过短，缺乏空间对比参考，也没有考虑全球变化的影响，难以量化开采扰动对矿区植被的影响程度。笔者采用时间分辨率较高的GIMMS AVHRR NDVI3g(1982—2013年)长时间序列数据，针对缺乏长期监测的煤炭开采区，设置多重对比实验区(10 km缓冲区、20 km缓冲区、自然生态校验区)，结合气候驱动因子(气温、降水)，分别从时间相关性、空间相关性、气候相关性3方面揭示矿区开采前、后地表植被覆盖的变化特征及植被生长期的变化规律，通过对比分析开采区相对于非开采区的植被变化差异，甄别自然变化和人类活动的影响，旨在为采煤区植被自然恢复、人工修复及其周边地区生态综合治理提供借鉴。

1 材料与方法

1.1 研究区

1.1.1直接影响区

彬长矿区位于黄陇侏罗纪煤田中部、鄂尔多斯盆地南端、太峪背斜以北、彬县和长武的交界处。坐标为东经107.65°～108.33°、北纬34.98°～35.32°，南北距离约36.5 km，东西距离约60.6 km，总面积978 km2。下辖大佛寺井田、小庄井田、文家坡井田等13座井工矿。研究区属渭北黄土高原沟壑地貌，地貌类型以黄土梁塬、黄土沟谷为主。其中黄土沟谷(49.79%)面积最大，其次为黄土梁(36.59%)，黄土塬(6.02%)面积较小，剩余为河流阶地、漫滩及水体。研究区处于暖温带半干旱大陆性季风气候区，年平均气温 11.1 ℃，年平均降水量为 561.3 mm。区内植被以中覆盖度为主，植被类型较为单调，主要包括刺槐林、白羊草、长芒草和农业植被，研究区土地利用类型以草地(44.63%)和耕地(33.06%)为主，林地(11.66%)次之。

1.1.2间接影响区(缓冲区)和自然生态校验区

煤炭开采不仅对直接影响区域造成一定生态影响，也会影响到周边区域生态系统，为更好的定性分析开采扰动对矿区植被覆盖的影响，笔者增设了10 km缓冲区、20 km缓冲区、自然生态校验区[16](图1)。其中，在选择自然生态校验区时，参照了遥感“伪不变特征区”(Pseudo Invariant Feature areas,PIFs)的设置原则[17]，即:① 必须是植被覆盖区;② 与研究区均处于同一纬度，尽可能远离扰动区，受人为因素影响较少;③ 与研究区高程相近，太阳光照条件相近。两区经差0.64°，纬差0.94°，几何中心相距150.5 km，鉴于两区原始植被覆盖类型并不一致，因此通过NDVI的相对变化程度(即增长率)来进行对比分析。

图1 研究区与生态校验区的地理位置Fig.1 Geographical location of the study area and the CK

1.2 实验数据

1.2.1遥感数据

研究所用的遥感数据为GIMMS NOAA/AVHRRR 的第1,2通道生成的NDVI3g数字影像，该NDVI数据最初是由美国航天局(NASA)全球监测与模型研究组发布的最大合成(Maximum Value Composites,MVC)数据(ftp://ftp.glcf.uniacs.umd.edu/glcf/GIMMS/)。图像分辨率为8 km×8 km，时间分辨率为15 d，时间跨度为1982-01—2013-12，空间范围为70°～140°E，15°～55°N，32 a涵盖了768幅影像。经过辐射校正和几何粗校正的NOAA/AVHRR，再进一步消除云、大气、太阳高度角等的部分干扰，运用国际通用的MVC(最大合成)方法获得半月合成数据。相对于其他NDVI数据源来说，GIMMS NDVI3g时间序列数据适合于长时间、大尺度植被活动变化的分析，用于分析我国的植被活动状况是合适的[8]。

1.2.2气象数据

气候数据为1982—2013年各省、市、自治区气候资料处理部门逐月上报的“地面气象记录月报表”的年均气温和降水资料，由于研究区域内的气象站点较稀疏，为消除极端数据所带来的影响，笔者结合当地气候地貌环境，以矿区几何中心点为圆心，300 km为半径作一包含34个气象站点的圆形插值区进行插值，土地利用图取自国家基础地理信息中心的30 m地表覆盖数据(GlobeLand30)，矢量文件为1∶25万县界数字线划图。

1.3 数据预处理

1.3.1遥感数据预处理

NDVI半月合成数据再处理方法主要有累加处理、平均值处理与最大值处理。累加处理是将分析某时间段内的NDVI值求和，可反映植被在某时间段内的生物量当量累加值;平均值处理是将某时间段内的NDVI值取平均值，该方法可消除气候异常对NDVI值的影响;最大值处理消除了云、雾等自然条件对NDVI值的影响，是将某时间段内的NDVI值取最大值。矿区的地理位置、所在区域的降水量、气温、以及作物生长状况等条件共同决定合成方法的选择。低植被覆盖度的半干旱地区(年均降雨量<400 mm,如神东矿区)，为突出植被覆盖状况，可取全年NDVI的最大值[18];中等植被覆盖度的半湿润地区(400 mm<年均降雨量<800 mm,如抚顺矿区)，可取植被生长期内的累加值[19-20];高植被覆盖度湿润地区(800 mm<年均降雨量,如淮南矿区)，可取全年NDVI的平均值[21]。

(1)累加值处理。研究区属中度植被覆盖半湿润地区，宜采用累加值方法，运用累加值也可估算区域NDVI总量年际变化，为开采扰动对矿区带来的植被变化提供一定依据，计算数据模型为

(1)

式中，ai为第i个像元NDVI半月合成值;m为累加次数，分别为2(月总值)、6(季总值)、24(年总值);bi为第i个像元累加值;n为像元个数。

(2)线性回归方程。STOW等用采用一元线性回归来计算植被的绿度变化率(Greenness Rate of Change,GRC)，GRC被定义为某时间段内NDVI一元线性回归方程的斜率S[22]。该斜率揭示了在一定时间内每个像元所在地区植被指数的变化趋势，公式为

(2)

(3)变化率。在单位时间内的变化量就是变化率。方精云[8]提出用变化率判断区域植被指数增加或减少的程度大小，并约定NDVI的变化率为

(3)

式中，S为32 a 的NDVI累积值线性回归直线的斜率;M为32 a的NDVI累加值均值，该参数用以对比分析开采扰动区与非开采区植被活动的趋势。

(4)植被生长期阈值。利用32 a的NDVI3g月均值数据，计算反映植被物候的3个参数:返青期(SOS)，枯黄期(EOS)和生长期(LOS)。月度曲线突然升高和降低处即为生长期开始或结束临界点，生长期长度为枯黄期开始时间减去返青期开始时间[23-24]。因原始数据月NDVI值时间分辨率为30 d，为进一步精确计算植被生长期，利用Matlab对32 a的月NDVI值以3 d为单位进行插值，可以计算出0.1个月的时间变化，即有效计算精度为±3 d。

1.3.2气象数据预处理

运用反距离加权法(Inverse Distance Weighted Method，IDW),对插值区34个气象站点年平均降水、气温进行插值，得到研究区32 a年均降水及气温(图2)。

(1)矿区年均降水分析。如图2(a)所示，彬长矿区32 a降水量波动较大，1982—2002年降水量整体呈下降趋势，与全国降水量变化趋势一致[25]，2003年后，降水量有上升之势。研究区平均年降水量在1997年降到极小值418.6 mm，1983年出现最高值818.3 mm，1981—2012年均降水量为592.56 mm，降水下降速率为21.55 mm/(10 a)。

(2)矿区年均气温分析。32 a研究区年均气温最低为9.81 ℃(1984年)，最高为12.12 ℃(2006年)和12.05 ℃(1998年)，32 a平均气温为11.11 ℃(图2(b))。矿区升温速率为0.47 ℃/(10 a)，高于陕北地区增温速率0.29 ℃/(10 a)[26]，更高于全国近50 a(1951—2004年)的平均增温速率0.25 ℃/(10 a)[27]，远高于同期全球气温变化速率0.12 ℃/(10 a)[28-29]。

图2 彬长矿区32 a气候变化趋势Fig.2 Trend of average temperature & precipitation for 32 years of Binchang coalfield

2 结果与讨论

2.1 时空相关性分析

2.1.1研究区月度NDVI值比较

(1)高斯拟合

将每个月的NDVI求平均，获得植被NDVI的年内特征曲线，图3为4个区域NDVI月均值拟合曲线，可见月度植被指数动态变化曲线符合高斯分布:

其中，y0为基线偏移;A为曲线下方的积分面积;x0为中央峰值;w=2σ近似于峰值半高宽的0.849，4个区域的拟合度R2分别为0.910 0(矿区)，0.921 8(对比区，也称校验区)，0.950 7(10 km缓冲区)，0.927 5(20 km缓冲区)。

图3 研究区32 a月度NDVI趋势Fig.3 Trend of monthly NDVI for 32 years in the study area

彬长矿区地处彬县与长武县境内，当地以落叶阔叶林和农业植物为主，年内变化曲线整体呈单峰型，1—3月份NDVI值处于低值期，4月份后增加迅速(NDVI值陡升)，5—8月达到并保持峰值状态，9月底至10月初NDVI值大幅度下降，直至2月份到达谷底。此变化特征与陕西省植被年内变化趋势相一致[30]。

(2)生长期分析

由图4可以判断，矿区和缓冲区(10 km,20 km)NDVI基数较高，有些年份12月时仍保持0.25的NDVI值，而且0.3的曲线极不稳定，其返青期和枯黄期宜采用相同阈值0.35;校验区受人为扰动较少，原始植被覆盖类型为草地，其返青期和枯黄期宜采用相同阈值0.20。

图4 研究区时间序列月NDVI晕渲图Fig.4 Shading map of monthly average NDVI in the study area

① 图4的“倒梯形”柱状结构表明，4个研究区返青期均有提前趋势，枯黄期有滞后趋势，生长期有延长趋势。② 对1982—2013年矿区、10 km缓冲区、20 km缓冲区和校验区植被返青期、枯黄期、生长期进行一元线性回归，得到4个区琙植被生长期变化趋势(图5)。32 a四区生长期均有延长趋势，生长期分别为182 d(矿区)，185 d(10 km缓冲区)，188 d(20 km缓冲区)和180 d(校验区)，非开采区生长期延长13±3 d(10 km缓冲区)，8±3 d(20 km缓冲区)，6±3 d(校验区)，矿区生长期延长11±3 d，线性增长率为0.34 d/a，1984年为矿区生长期最小年(156 d)，生长期最大年是2004年(175 d)。

图5 4类区域植被生长期变化Fig.5 LOS change in four types of regions

2.1.2年NDVI累加值比较

对1982—2013年裁剪、去除DN值后的遥感图像运用IDL进行累加值处理，得到32 a累积值遥感图像，进而提取出每年各个区域的NDVI累积值，可反映不同研究区植被生长状况，为矿区植被自然恢复、人工修复提供参考。

图6 32 a NDVI累加值变化趋势Fig.6 Slope of NDVI during 32 years

根据图6，研究区年NDVI总值对比趋势显示:

(1)通过0.01的显著性水平检验，矿区与10 km缓冲区、20 km缓冲区、校验区NDVI的相关系数分别为0.979 1,0.961 2,0.651 3;说明矿区与缓冲区植被响应高度一致，而矿区植被活动变化趋势与校验区相关性不大，这样更加适合相对变化的比较。32 aNDVI总量从高到低依次为20 km缓冲区(358.33)、10 km缓冲区(237.42)、直接影响区(158.43)、生态校验区(121.93)，4区植被覆盖均在波动中呈增加之势，其中矿区植被活动变化幅度最弱。2006年开始，校验区植被累积值开始逐渐缩小与矿区的差距，至2013年两区年际植被总量基本一致，在自然生态环境条件下，校验区植被自然生长，矿区植被受到人类活动影响，导致两区植被总量差距缩小。

(2)对研究区32 a连续时间序列NDVI累积值做一元线性回归，由回归方程(式(2))得到的绿度变化率(S)分别为0.659 2(矿区)、0.809 0(10 km缓冲区)、1.174 8(20 km缓冲区)、0.741 1(校验区)。32 a NDVI总量均值分别为158.43，237.42，358.33，121.93，结合年变化率式(式(3))，可得研究区32 a NDVI增长百分比为10.49%(20 km缓冲区)<10.9%(10 km缓冲区)<13.31%(矿区)<19.45%(生态校验区)，矿区因受采煤影响植被增长率显著低于无人工开采干扰的自然生态区。若以校验区NDVI上升速率为自然增长率，可初步认为亏欠的部分是采矿活动的影响，则开采活动对矿区、10 km缓冲区和20 km缓冲区NDVI增长率的贡献分别为-6.14%，-8.55%，-8.96%，对覆盖度高的区域影响似更为显著。

2.2 气候相关性分析

降水和气温是植被响应气候变化最主要的2个驱动因子，笔者选取32 a气象数据与研究区NDVI做相关性统计分析(图7)。结果表明:矿区和缓冲区植被与气温呈显著性正相关，相关系数分别为0.443 0,0.412 1,0.501 0(p<0.05)，均通过显著性检验。校验区对气温响应敏感性较差。4区与降水呈不显著负相关:-0.094 9,-0.085 9,-0.152 8,-0.085 1，植被对气温的响应高于降水，这与以往研究结论有一致性[27-30]，气候因子对研究区植被做“超前-滑动增量”表明，降水对植被生长有2～3 a的滞后性(相关系数达-0.469,0.587)，气温对于植被生长无明显滞后性，秦超等对陕西省植被响应气候变化研究中也得出相同结论[31]。

图7 研究区AVHRR NDVI与降水(上)/气温(下)相关性曲线Fig.7 Correlation curves of AVHRR NDVI and precipitation (above)/temperature (below) in study area

2.3 讨 论

(1)确定NDVI生长期阈值的方法很多(如滑动平均值法[32-33]，NDVI阈值法[34-36]，温度界限法[37]等)，研究方法不同，得出结论也会有所偏差，笔者所采用的NDVI阈值法适合于我国北方植被生长期变化规律研究。由于植被类型不同，不同植被区域的NDVI阈值也有所不同。气候驱动因子只是植被活动的影响因素的一部分，且本文只研究了气候影响因子里最主要的降水和气温，并未研究蒸散量，湿度，日照强度等因素对植被的作用效果。研究结果表明，32 a来研究区植被活动显著增强，这与MYNENI,FANG,TUCKER,LOS等研究成果一致[4-6，8-9];针对RASIM，PROSPER，HUANG，徐占军、马超等对矿区植被研究表明[10-15]，植被在受到开采扰动后生长受到较大影响，生物量呈减少之势，而本文研究矿区植被NDVI累积量呈增加趋势，但低于自然环境校验区NDVI增长率，仍能说明煤炭开采对植被造成了一定影响。

(2)矿区植被生长期延长趋势高于校验区，一方面，说明植被在全球气候变化条件下，生长期有延长倾向;植被受气候变化影响返青期提前，生长期延长的特点，与农牧交错带、陕西及我国北方等地植被生长期研究结论一致[38-41]。而彬长矿区属粮矿复合区，工矿区工农业生产均会导致植被生长周期发生变化，有如下方面:① 工矿区人类生产生活实践，可能形成矿区热岛，改变局部水热条件[42-45]，使区域植被返青期提前，枯黄期滞后;② 矿区绿化、植树造林、退耕还林和土地复垦实践，都可能引入人工植被、改良生态结构[46]，使区域植被返青期提前，枯黄期滞后;③ 矿区周边为传统农业种植区，粮矿复合区的植被物候受农业生产影响较大，越冬作物或大秋作物都可能使区域植被返青期提前，枯黄期滞后;④ 煤矿开采造成土壤理化性质的变化导致植被生理周期被延长。如土壤氮、磷等养分元素损失，开采塌陷造成土壤表面裂缝生成，土壤持水保肥能力弱，植物根系受损，光合作用能力变差，不利于植物生长及植被恢复，植被生理周期被延长[47]。这一点可以从“矿区植被生长期延长，但NDVI增长率却下降”的结论得到一定的支持。

综上所述，矿区及其缓冲区的植被生长期由于人工干扰，不同程度地被二度延长，凸显了人类活动对生态环境的重大影响。

3 结 论

(1)32 a(1982—2013)彬长矿区月度NDVI变化表明，研究区年内植被变化呈双峰型，5—8月份植被指数较高，8月份达到峰值，9月底植被指数骤降，2月份植被指数达到谷值。年际上，4区平均NDVI总量为20 km缓冲区(358.33)、10 km缓冲区(237.42)、直接影响区(158.43)、生态校验区(121.93)，32 a间NDVI增长百分比为生态校验区(19.45%)>矿区(13.31%)>10 km缓冲区(10.9%)>20 km缓冲区(10.49%)，自然环境生长下的植被活动比煤矿开采区活跃，采煤活动对矿区植被造成了一定影响。

(2)1982—2013年研究区植被生长期研究表明:矿区及缓冲区在4月初进入返青期，9月底进入枯黄期，校验区进入返青期日期为4月上旬，枯黄期日期为10月上旬，四区生长期长度在185 d左右。时间序列拟合曲线表明，四区返青期提前，枯黄期推迟，生长期延长，矿区延长天数高于生态校验区5 d。

(3)植被响应气候因子相关性分析:研究区受全球气候变化影响，气温升高，降水量减少，处于半湿润地区的矿区与缓冲区植被对气温的响应更为敏感;校验区NDVI受其气温影响较小，植被与同期降水量相关性不大，但降水对其有明显滞后性。植被受气候变化影响返青期提前，生长期延长的特点，与农牧交错带、陕西及我国北方等地植被生长期研究结论一致。