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

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

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

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

1 研究區(qū)概況

青海省面積約72.12 萬km2，平均海拔在3000 m以上，地勢起伏西部較小，東部較大，且地貌類型多樣，差異性大。 青海省年降水量為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ù)進行了軌道偏移、大氣水汽、輻射、氣溶膠、除云、幾何畸變等校正與處理［12］。 收集研究區(qū)1982—2015年的GIMMS NDVI 數(shù)據(jù)，其空間分辨率為8 km，時間分辨率為15 d。MODIS 是美國EOS 系列衛(wèi)星的重要傳感器，擁有36個波段，1～2 d 覆蓋全球一次，收集研究區(qū)2000—2018年MODIS 的MOD13A2 NDVI V6 產(chǎn)品數(shù)據(jù)，其空間分辨率為1 km，時間分辨率為16 d。 1 ∶100 萬《中國植被圖集》、中國植被區(qū)劃數(shù)據(jù)、青海省行政區(qū)劃數(shù)據(jù)來源于中國科學院資源環(huán)境科學與數(shù)據(jù)中心。 DEM 數(shù)據(jù)來源于美國航天航空局，共8 景數(shù)據(jù)。 氣象數(shù)據(jù)來源于國家氣象科學數(shù)據(jù)中心中國氣象數(shù)據(jù)網(wǎng)。 對GIMMS NDVI、MODIS NDVI 數(shù)據(jù)進行格式轉(zhuǎn)換、拼接等處理后，利用行政區(qū)劃數(shù)據(jù)裁切和投影，獲得研究區(qū)GIMMS NDVI、MODIS NDVI、1 ∶100 萬植被圖集、植被區(qū)劃數(shù)據(jù)和DEM 數(shù)據(jù)。

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

2.2 研究方法

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

圖2 技術路線

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

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

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

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

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

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

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

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

基于SG 與小波融合的自適應降噪處理方法對GIMMS NDVI 與MODIS NDVI 數(shù)據(jù)進行降噪處理，并對其降噪效果進行對比分析和評價；計算獲得兩種數(shù)據(jù)逐年最大值合成數(shù)據(jù)，再將MODIS NDVI 數(shù)據(jù)重采樣為8 km，并選取重采樣后2000—2015年的數(shù)據(jù)與對應的GIMMS NDVI 數(shù)據(jù)進行逐像元回歸分析，同時進行回歸效果評價。 依據(jù)各個像元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ù)中隨機選擇9 個像元，利用傳統(tǒng)高斯函數(shù)、SG 算法、Hants 算法，以及本文提出的SG 與小波融合自適應降噪處理方法對其進行處理，利用式（1）～式（3）計算各方法的加權綜合評價指標WCEI，結(jié)果見表1。 由表1可知，基于SG 與小波融合的自適應降噪處理方法的WCEI值最大，說明該方法降噪效果最為理想，并且在降噪處理過程中不用按照年度設置平滑窗口和確定去噪閾值，計算效率較高。

表1 NDVI 序列數(shù)據(jù)降噪方法評價

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

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

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

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

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

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

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

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

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

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

分別提取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ù)在描述植被覆蓋變化過程中具有較強的一致性。 因此，基于8 km 尺度的1982—2018年NDVI數(shù)據(jù)，利用Mann-Kendall 非參數(shù)檢驗方法對研究區(qū)植被覆蓋時空變化趨勢進行分析，通過置信度為99%檢驗的像元存在顯著變化趨勢，通過置信度為90%檢驗的像元呈增加或減少趨勢，見圖7。 由圖7可以看出，近37 a 研究區(qū)北部植被覆蓋增加，南部及柴達木盆地部分地區(qū)植被覆蓋減少。 植被覆蓋呈減少趨勢的面積占12.3%，其中顯著減少的占5.3%；植被覆蓋呈增加趨勢的面積占52.1%，其中顯著增加的占35.7%。

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

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

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

4 討論

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

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

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

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

5 結(jié)語

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