基于多源遙感數(shù)據(jù)的青海省高寒植被覆蓋時(shí)空變化分析

趙健赟，丁圓圓，杜 梅，楊 靜，宋 芊

（1.青海大學(xué) 地質(zhì)工程系，青海 西寧 810016； 2.青海省青藏高原北緣新生代資源環(huán)境重點(diǎn)實(shí)驗(yàn)室，青海 西寧 810016； 3.青海省剛察縣氣象站，青海 剛察 812300）

青藏高原特殊的地理位置和地貌條件孕育形成了濕潤、干旱、熱帶和寒帶等多種氣候類型，在區(qū)域水源涵養(yǎng)、生物多樣性保護(hù)、水土保持等方面發(fā)揮著不可替代的作用。 青藏高原生態(tài)系統(tǒng)敏感而脆弱，受全球氣候變化和人類活動(dòng)影響，青藏高原地區(qū)草地退化、水土流失等問題日益突出，認(rèn)識(shí)青藏高原生態(tài)安全功能變化與成因，是保障青藏高原生態(tài)安全的迫切需求［1-5］。高寒植被作為高原生態(tài)系統(tǒng)的重要組成部分，在物質(zhì)能量交換過程中扮演著重要角色，植被覆蓋變化可以在一定程度上反映生態(tài)環(huán)境的變化，因此植被覆蓋變化研究對(duì)于青藏高原生態(tài)環(huán)境變化規(guī)律的探索具有重要意義［6-9］。

歸一化植被指數(shù)NDVI可以準(zhǔn)確反映植被綠度、光合作用強(qiáng)度，對(duì)植被生長(zhǎng)量和長(zhǎng)勢(shì)十分敏感，是反映地表植被生長(zhǎng)狀況的重要指標(biāo)。 目前，已有學(xué)者利用NDVI對(duì)青藏高原開展了較多植被覆蓋變化方面的研究。 Li 等［10］指出1982—2008年青藏高原NDVI總體呈上升趨勢(shì)，季節(jié)變化趨勢(shì)不明顯；Yang 等［2］研究1982—1999年青藏高原草地覆蓋地區(qū)的植被指數(shù)變化發(fā)現(xiàn)，NDVI在生長(zhǎng)季顯著增大，春季NDVI增長(zhǎng)率最大；趙健赟等［11］對(duì)青海高原植被覆蓋時(shí)空變化特征及其驅(qū)動(dòng)因子進(jìn)行了分析和研究。 目前時(shí)間序列NDVI數(shù)據(jù)主要來自AVHRR 和MODIS 這兩種傳感器，獲取的數(shù)據(jù)在時(shí)間尺度和空間尺度上均存在較大差異，因此大多數(shù)研究集中于單一尺度區(qū)域植被覆蓋變化分析，對(duì)于集成兩種NDVI數(shù)據(jù)開展更長(zhǎng)時(shí)間序列的變化分析，以及不同海拔、不同類型植被覆蓋變化的差異性等方面有待深入研究。 本文以青藏高原東北部青海省為例，利用GIMMS NDVI、MODIS NDVI 等多源遙感數(shù)據(jù)，通過自適應(yīng)降噪、回歸建模等方法減弱數(shù)據(jù)間的差異，進(jìn)而分析青海省高寒植被覆蓋時(shí)空變化特征與規(guī)律。

1 研究區(qū)概況

青海省面積約72.12 萬km2，平均海拔在3000 m以上，地勢(shì)起伏西部較小，東部較大，且地貌類型多樣，差異性大。 青海省年降水量為200 ～500 mm，主要集中在6—9月，總體降水偏少且分布不均，降水量由東南向西北遞減。 青海省年平均氣溫在0 ℃以下，中部地區(qū)晝夜溫差較大，月均氣溫7月最高、1月最低，屬高原大陸性氣候區(qū)。 青海省土地資源類型多樣，土地分布具有顯著高原地帶性，主要植被類型有高寒草原、高寒草甸、高山植被、灌叢等，研究區(qū)主要植被類型空間分布見圖1。

圖1 研究區(qū)主要植被類型空間分布

2 數(shù)據(jù)來源與研究方法

2.1 數(shù)據(jù)來源

GIMMS NDVI 數(shù)據(jù)采用NOAA 系列氣象衛(wèi)星AVHRR 傳感器的GIMMS NDVI 第三代（NDVI3g）全球覆蓋產(chǎn)品數(shù)據(jù)集，GIMMS NDVI 數(shù)據(jù)進(jìn)行了軌道偏移、大氣水汽、輻射、氣溶膠、除云、幾何畸變等校正與處理［12］。 收集研究區(qū)1982—2015年的GIMMS NDVI 數(shù)據(jù)，其空間分辨率為8 km，時(shí)間分辨率為15 d。MODIS 是美國EOS 系列衛(wèi)星的重要傳感器，擁有36個(gè)波段，1～2 d 覆蓋全球一次，收集研究區(qū)2000—2018年MODIS 的MOD13A2 NDVI V6 產(chǎn)品數(shù)據(jù)，其空間分辨率為1 km，時(shí)間分辨率為16 d。 1 ∶100 萬《中國植被圖集》、中國植被區(qū)劃數(shù)據(jù)、青海省行政區(qū)劃數(shù)據(jù)來源于中國科學(xué)院資源環(huán)境科學(xué)與數(shù)據(jù)中心。 DEM 數(shù)據(jù)來源于美國航天航空局，共8 景數(shù)據(jù)。 氣象數(shù)據(jù)來源于國家氣象科學(xué)數(shù)據(jù)中心中國氣象數(shù)據(jù)網(wǎng)。 對(duì)GIMMS NDVI、MODIS NDVI 數(shù)據(jù)進(jìn)行格式轉(zhuǎn)換、拼接等處理后，利用行政區(qū)劃數(shù)據(jù)裁切和投影，獲得研究區(qū)GIMMS NDVI、MODIS NDVI、1 ∶100 萬植被圖集、植被區(qū)劃數(shù)據(jù)和DEM 數(shù)據(jù)。

青海省不同海拔高度氣候差異顯著，高程2000 m以下的河湟谷地為暖區(qū)，高程3000 m 的柴達(dá)木盆地等地區(qū)為次暖區(qū)，高程4000 m 以上的祁連山區(qū)和青南高原地區(qū)為冷區(qū)。 考慮到2000 m 以下河湟谷地面積極小，將研究區(qū)分為3000 m 以下、3000 ～4000 m和4000 m 以上3 個(gè)高程區(qū)劃帶。

2.2 研究方法

研究采用的技術(shù)路線見圖2，NDVI序列數(shù)據(jù)降噪采用基于SG 與小波融合的自適應(yīng)降噪處理方法，算法流程見圖3，其中ci為濾波系數(shù)，di為NDVI與總趨勢(shì)線間的距離，dmax為NDVI與總趨勢(shì)線間的最大距離。

圖2 技術(shù)路線

為評(píng)估所采用NDVI序列數(shù)據(jù)降噪方法與傳統(tǒng)高斯函數(shù)、SG 算法、Hants算法等降噪方法的優(yōu)劣，構(gòu)建加權(quán)綜合評(píng)價(jià)指標(biāo)WCEI對(duì)其進(jìn)行評(píng)價(jià)［13］。

式中：RMSEi為均方根誤差；SNRi為信噪比；NDVIi為原始數(shù)據(jù)序列；NDVI*i為降噪后的數(shù)據(jù)序列；n為降噪序列數(shù)據(jù)的總數(shù)量。

WCEI的權(quán)重主要依據(jù)NDVI真值無法準(zhǔn)確獲取以及信噪比SNR越大降噪效果越好的基本事實(shí)確定?；诮翟牒蟮腘DVI數(shù)據(jù)集，計(jì)算逐像元逐年NDVI最大值：

利用回歸分析法建立GIMMS NDVI 與MODIS NDVI 之間的轉(zhuǎn)換關(guān)系：

式中：β0、β1為回歸系數(shù)；εi為殘差；為GIMMS NDVI 數(shù)據(jù)；為MODIS NDVI 數(shù)據(jù)。

同時(shí)，利用GIMMS NDVI 和MODIS NDVI 重疊的2000—2015年數(shù)據(jù)，通過回歸均方誤差（MSE）對(duì)建立的回歸模型進(jìn)行可靠性評(píng)估：

利用Mann-Kendall 非參數(shù)檢驗(yàn)方法進(jìn)行植被覆蓋時(shí)空變化趨勢(shì)分析，其統(tǒng)計(jì)檢驗(yàn)指標(biāo)如下：

式中：S為檢驗(yàn)統(tǒng)計(jì)量；Xp、Xq為時(shí)序數(shù)據(jù)；l為時(shí)序長(zhǎng)度；Z為標(biāo)準(zhǔn)化檢驗(yàn)統(tǒng)計(jì)量（l ＞8 時(shí)，S為近似正態(tài)分布）；V（S） 為方差。

基于SG 與小波融合的自適應(yīng)降噪處理方法對(duì)GIMMS NDVI 與MODIS NDVI 數(shù)據(jù)進(jìn)行降噪處理，并對(duì)其降噪效果進(jìn)行對(duì)比分析和評(píng)價(jià)；計(jì)算獲得兩種數(shù)據(jù)逐年最大值合成數(shù)據(jù)，再將MODIS NDVI 數(shù)據(jù)重采樣為8 km，并選取重采樣后2000—2015年的數(shù)據(jù)與對(duì)應(yīng)的GIMMS NDVI 數(shù)據(jù)進(jìn)行逐像元回歸分析，同時(shí)進(jìn)行回歸效果評(píng)價(jià)。 依據(jù)各個(gè)像元2000—2015年重采樣后的MODIS NDVI 數(shù)據(jù)和回歸方程，估算2016—2018年的GIMMS NDVI 數(shù)據(jù)，從而建立研究區(qū)1982—2018年8 km 尺度逐年NDVI最大值合成數(shù)據(jù)。

3 結(jié)果與分析

3.1 NDVI 序列數(shù)據(jù)重構(gòu)效果分析

在GIMMS NDVI 與MODIS NDVI 原始數(shù)據(jù)中隨機(jī)選擇9 個(gè)像元，利用傳統(tǒng)高斯函數(shù)、SG 算法、Hants 算法，以及本文提出的SG 與小波融合自適應(yīng)降噪處理方法對(duì)其進(jìn)行處理，利用式（1）～式（3）計(jì)算各方法的加權(quán)綜合評(píng)價(jià)指標(biāo)WCEI，結(jié)果見表1。 由表1可知，基于SG 與小波融合的自適應(yīng)降噪處理方法的WCEI值最大，說明該方法降噪效果最為理想，并且在降噪處理過程中不用按照年度設(shè)置平滑窗口和確定去噪閾值，計(jì)算效率較高。

表1 NDVI 序列數(shù)據(jù)降噪方法評(píng)價(jià)

選取2000—2015年重采樣至8 km 的MODIS NDVI 數(shù)據(jù)與對(duì)應(yīng)的GIMMS NDVI 數(shù)據(jù)進(jìn)行逐像元回歸分析，得到回歸方程：

依據(jù)該回歸方程，利用式（6）與GIMMS NDVI 原始數(shù)據(jù)，計(jì)算回歸均方誤差MSE，回歸值與殘差值見見表2，計(jì)算得到的均方誤差MSE為0.0098。 建立的MODIS NDVI 與GIMMS NDVI 數(shù)據(jù)的回歸方程殘差在±0.010 附近，最大殘差為0.024，該值與2000年MODIS NDVI 原始數(shù)據(jù)缺失有關(guān)，且回歸均方誤差小于0.001，因此利用該回歸方程重構(gòu)的NDVI序列數(shù)據(jù)是可靠的。

表2 NDVI 序列數(shù)據(jù)回歸均方誤差

3.2 研究區(qū)NDVI 時(shí)序變化特征

基于建立的1982—2018年逐年NDVI最大值合成數(shù)據(jù)，統(tǒng)計(jì)研究區(qū)和各高程分區(qū)的NDVI均值，繪制年際變化曲線及趨勢(shì)線，見圖4。 研究區(qū)1982—2018年NDVI整體呈上升趨勢(shì)，年增長(zhǎng)速率為0.04%；在高程3000～4000 m 區(qū)域，NDVI上升趨勢(shì)最明顯，年增長(zhǎng)速率為0.07%。 由此表明，研究區(qū)37 a 間植被覆蓋整體有改善趨勢(shì)，中高海拔地區(qū)植被覆蓋改善較為顯著。

圖4 研究區(qū)年際總體及各高程分區(qū)NDVI 變化曲線

3.3 主要植被類型NDVI年際變化特征

基于建立的1982—2018年逐年NDVI最大值合成數(shù)據(jù)，統(tǒng)計(jì)研究區(qū)各主要植被類型的NDVI均值，繪制年際變化曲線及趨勢(shì)線，見圖5。 1982—2018年研究區(qū)主要植被類型中，高寒草原植被的NDVI值上升趨勢(shì)最為明顯，年增長(zhǎng)速率為0.08%；灌叢、高寒草甸和高山植被NDVI值的年增長(zhǎng)速率均在0.04%左右。 由此表明，研究區(qū)37 a 間，各主要植被類型的覆蓋均有所改善，高寒草原植被覆蓋改善最為顯著。

圖5 主要植被年際NDVI 變化曲線

3.4 青海省植被覆蓋多尺度時(shí)空變化特征

分別提取8 km 尺度GIMMS NDVI 數(shù)據(jù)與1 km 尺度MODIS NDVI 數(shù)據(jù)重疊年份（2000—2015年）最大值合成數(shù)據(jù)，繪制年際變化過程曲線，見圖6。 由圖6可以看出，雖然兩種數(shù)據(jù)NDVI整體存在0.02 左右的差異，但回歸曲線的斜率接近，即兩種尺度的數(shù)據(jù)在描述植被覆蓋變化過程中具有較強(qiáng)的一致性。 因此，基于8 km 尺度的1982—2018年NDVI數(shù)據(jù)，利用Mann-Kendall 非參數(shù)檢驗(yàn)方法對(duì)研究區(qū)植被覆蓋時(shí)空變化趨勢(shì)進(jìn)行分析，通過置信度為99%檢驗(yàn)的像元存在顯著變化趨勢(shì)，通過置信度為90%檢驗(yàn)的像元呈增加或減少趨勢(shì)，見圖7。 由圖7可以看出，近37 a 研究區(qū)北部植被覆蓋增加，南部及柴達(dá)木盆地部分地區(qū)植被覆蓋減少。 植被覆蓋呈減少趨勢(shì)的面積占12.3%，其中顯著減少的占5.3%；植被覆蓋呈增加趨勢(shì)的面積占52.1%，其中顯著增加的占35.7%。

圖6 8 km 與1 km 尺度下NDVI年際變化曲線

基于建立的2000—2018年1 km 尺度逐年MODIS NDVI 最大值合成數(shù)據(jù)，采用Mann-Kendall 非參數(shù)檢驗(yàn)方法對(duì)研究區(qū)植被覆蓋時(shí)空變化趨勢(shì)進(jìn)行分析，結(jié)果見圖8。 由圖8可以看出，研究區(qū)西北部和中東部地區(qū)植被覆蓋增加，南部、西北部和柴達(dá)木盆地部分地區(qū)植被覆蓋減少。 植被覆蓋呈減少趨勢(shì)的面積占4.5%，其中顯著減少的占1.0%；植被覆蓋呈增加趨勢(shì)的面積占53.8%，其中顯著增加的占26.7%。 同時(shí)，結(jié)合8 km 尺度的研究區(qū)植被覆蓋時(shí)空變化情況，發(fā)現(xiàn)在不同時(shí)空尺度下研究區(qū)植被覆蓋變化的空間分布一致，研究結(jié)果具有一定可靠性。

圖8 1 km 尺度下2000—2018年研究區(qū)植被覆蓋時(shí)空變化趨勢(shì)

4 討論

在青藏高原及其局部地區(qū)植被覆蓋變化研究中，諸多學(xué)者將NDVI作為基礎(chǔ)數(shù)據(jù)，而目前NDVI值的來源主要有GIMMS NDVI、SPOT VGT NDVI、MODIS NDVI 和Landsat NDVI，GIMMS NDVI 數(shù)據(jù)時(shí)間序列最長(zhǎng)，SPOT VGT NDVI 與MODIS NDVI 數(shù)據(jù)時(shí)間序列較短，而Landsat NDVI 數(shù)據(jù)雖然具有最高的空間分辨率，但無法構(gòu)建連續(xù)有效的數(shù)據(jù)開展長(zhǎng)時(shí)序分析。 為研究植被的長(zhǎng)期變化趨勢(shì)，需要集成上述不同時(shí)空分辨率的NDVI數(shù)據(jù)集進(jìn)行多源NDVI數(shù)據(jù)重構(gòu)，其方法主要是建立不同數(shù)據(jù)集之間的統(tǒng)計(jì)模型［14-16］。 如Gallo 等［17］建立了灌木、草地、常綠闊葉林等土地覆蓋類型AVHRR NDVI 與MODIS NDVI 數(shù)據(jù)的關(guān)系模型；張宏 斌 等［18］在宏觀尺度上根據(jù)NOAA NDVI 與MODIS NDVI 數(shù)據(jù)之間的相關(guān)關(guān)系建立了線性空間尺度轉(zhuǎn)換模型，取得了較好的應(yīng)用效果；張戈麗等［19］通過建立MODIS NDVI 數(shù)據(jù)與NOAA NDVI 數(shù)據(jù)線性回歸方程實(shí)現(xiàn)了NOAA NDVI 數(shù)據(jù)的插補(bǔ)與延長(zhǎng)。 本研究在SG 與小波融合的自適應(yīng)處理方法對(duì)GIMMS NDVI 與MODIS NDVI 數(shù)據(jù)整體降噪的基礎(chǔ)上，采用線性回歸模型實(shí)現(xiàn)了GIMMS NDVI 數(shù)據(jù)的延長(zhǎng)，且降噪效果優(yōu)于傳統(tǒng)方法，重構(gòu)數(shù)據(jù)的時(shí)間、空間變化一致性較好。

目前采用的植被覆蓋時(shí)空變化分析方法有線性回歸法［14］、變異系數(shù)法［16］、綠度變化率法［20］等，這些方法在獲得基礎(chǔ)參數(shù)的同時(shí)，還需要輔以顯著性檢驗(yàn)才能確定植被覆蓋變化的顯著性變化趨勢(shì)。 本研究所采用的Mann-Kendall 法是一種經(jīng)典的非參數(shù)趨勢(shì)檢驗(yàn)方法，無須樣本遵從特定分布，亦不受少量離群值的干擾，適用于植被覆蓋、降水等非正態(tài)分布數(shù)據(jù)的趨勢(shì)分析。 在研究區(qū)植被覆蓋時(shí)空變化方面，Zhao 等［21-23］利用GIMMS NDVI 和MODIS NDVI 數(shù)據(jù)開展了1982—2013年青海省植被覆蓋變化研究，發(fā)現(xiàn)植被覆蓋整體有改善趨勢(shì)，尤其是高海拔地區(qū)、祁連山東部地區(qū)植被覆蓋改善最為明顯。 本研究結(jié)果表明，青海省1982—2018年植被覆蓋整體有改善趨勢(shì)，中高海拔地區(qū)改善較為顯著，與上述不同時(shí)空尺度下植被覆蓋變化的研究結(jié)果一致，同時(shí)研究發(fā)現(xiàn)盡管各主要植被類型的覆蓋均有所改善，但高寒草原植被覆蓋改善最為顯著。

根據(jù)1982—2018年研究區(qū)植被覆蓋時(shí)空變化趨勢(shì)分析結(jié)果，提取1982—2018年青海省有增加或顯著增加趨勢(shì)的氣象站的年平均氣溫、風(fēng)速、年降水量和日照時(shí)數(shù)等氣象因子，與所在像元的NDVI進(jìn)行相關(guān)分析，結(jié)果見表3。 由表3可以看出，1982—2018年研究區(qū)大部分地區(qū)植被覆蓋的增加與氣溫、降水量的增加顯著相關(guān)，且大部分地區(qū)與風(fēng)速、日照時(shí)數(shù)存在負(fù)相關(guān)關(guān)系；在3000 m 以上的部分草原地區(qū)，植被覆蓋與風(fēng)速、日照時(shí)數(shù)具有顯著相關(guān)關(guān)系。 另外，研究區(qū)實(shí)施的退耕還林（草）、退牧還草政策，2000—2001年為試點(diǎn)階段，2002—2006年全面啟動(dòng)，2007—2015年為成果鞏固階段。 2000—2015年實(shí)施退耕還林（草）1934.3 km2（還林1740.1 km2，還草194.2 km2），涉及青海省44 個(gè)縣（市、區(qū)、場(chǎng)）、327 個(gè)鄉(xiāng)鎮(zhèn)、3911 個(gè)行政村，說明實(shí)施的退耕還林（草）、退牧還草等政策有效驅(qū)動(dòng)了研究區(qū)植被覆蓋的增加。 從2000—2018年8 km 尺度和1 km 尺度的植被覆蓋變化研究結(jié)果來看，盡管植被覆蓋變化的空間分布情況基本一致，但數(shù)量存在較大差異，說明尺度效應(yīng)對(duì)研究結(jié)果影響較大。

表3 NDVI 與各氣象因子的相關(guān)系數(shù)

5 結(jié)語

基于不同時(shí)空尺度的GIMMS NDVI 和MODIS NDVI 數(shù)據(jù)，構(gòu)建了青海省1982—2018年長(zhǎng)時(shí)序NDVI數(shù)據(jù)集，對(duì)研究區(qū)不同海拔、不同類型高寒植被覆蓋變化進(jìn)行了分析，得到以下結(jié)論：研究區(qū)1982—2018年NDVI整體呈上升趨勢(shì)，在海拔3000 ～4000 m 區(qū)域，NDVI上升趨勢(shì)最為明顯，年增長(zhǎng)速率為0.07%；研究區(qū)主要植被類型區(qū)劃中，高寒草原的NDVI上升趨勢(shì)最為明顯，年增長(zhǎng)速率為0.08%；灌叢、高寒草甸和高山植被NDVI的年增長(zhǎng)速率均在0.04%左右；1982—2018年8 km 尺度下，研究區(qū)北部植被覆蓋增加，南部及柴達(dá)木盆地部分地區(qū)植被覆蓋減少，植被覆蓋呈減少趨勢(shì)的面積占12.3%，植被覆蓋呈增加趨勢(shì)的面積占52.1%。 2000—2018年1 km 尺度下，研究區(qū)西北部和中東部植被覆蓋增加，南部、柴達(dá)木盆地部分地區(qū)植被覆蓋減少，植被覆蓋呈減少趨勢(shì)的面積占4.5%，植被覆蓋呈增加趨勢(shì)的面積占53.8%。 1982—2018年研究區(qū)植被覆蓋的增加與氣溫升高、降水量增加具有顯著相關(guān)關(guān)系，大部分地區(qū)與風(fēng)速、日照時(shí)數(shù)存在負(fù)相關(guān)關(guān)系，在高程3000 m 以上的部分草原地區(qū)，植被覆蓋與風(fēng)速、日照時(shí)數(shù)的變化存在較為顯著的相關(guān)關(guān)系。退耕還林（草）、退牧還草等政策有效驅(qū)動(dòng)了研究區(qū)植被覆蓋的變化，尺度效應(yīng)對(duì)研究結(jié)果影響較大，實(shí)現(xiàn)有效的尺度轉(zhuǎn)換并構(gòu)建長(zhǎng)時(shí)序、高時(shí)空分辨率的地表覆蓋數(shù)據(jù)是下一步研究的重點(diǎn)。