基于時間序列分割算法的雅魯藏布江流域NDVI(1985—2018)變化模式研究

王 塞,王思詩,樊風雷,2,*

1 華南師范大學地理科學學院, 廣州 510631 2 西藏大學高原地表環境遙感監測聯合實驗室, 拉薩 850000 3 天津大學環境科學與工程技術學院, 天津 300350

植被是聯結水分、大氣和土壤的重要紐帶,也是評價生物多樣性和生態服務的重要指標,在全球碳平衡和氣候變化中扮演重要角色[1-3]。雅魯藏布江流域生態環境脆弱,對氣候變化和人類活動的響應顯著[4]。近30年來,流域所屬的青藏高原年平均氣溫每10年上升0.7℃,增溫趨勢總體呈現由南向北逐漸降低,降水量受復雜地形影響在不同區域差異顯著,流域氣候暖干化趨勢明顯。由此,干旱程度加劇導致的植被稀疏化對流域脆弱的生態系統產生不利影響。此外,土地過度開墾、過度放牧和城市化進程等人類活動也加劇了流域植被退化[5-7]。因此,雅魯藏布江流域的植被變化及其影響因素倍受學界關注。目前研究的方法與途徑主要是,以遙感影像為數據源,通過對不同類型的植被覆蓋指標與氣溫、降水等氣候因子進行耦合測度,進而判別該流域內生態環境狀況[6,8-10]。開展高原高寒植被物候期變化研究,討論氣候變化對植被生長的影響與作用規律[11]。然而,對雅魯藏布江流域植被變化模式研究鮮有相關報道。

隨著遙感技術的發展,中分辨率遙感數據已經成為在國家和大陸尺度上進行環境監測的重要數據源[12]。基于中分辨率長時間序列遙感衛星影像是對大范圍區域進行長期監測,分析其土地利用變化趨勢,并對植被變化模式進行判別,已成為評價區域植被變化的重要手段[13-14]。歸一化植被指數(Normalized Difference Vegetation Index,NDVI)在植被變化研究中得到廣泛應用。NDVI反映植被冠層在近紅外與紅光波段的吸收和散射作用,與植被的光合作用與生理生化反應密切相關,因此既能監測植被對全球氣候變化的響應,又能識別城市擴張對植被的影響[15]。中分辨率長時間序列NDVI數據集成以上優勢,已逐漸成為研究植被變化重要的數據源。NDVI時間序列變化監測算法按原理可分為閾值、差異、分割、軌跡分類、統計和回歸六類[16]。Landsat數據源的時間序列分割算法(LandTrendr)在森林植被覆蓋變化檢測與生長模式判別等研究方面取得了豐碩的成果[13,17-18]。Google Earth Engine(GEE)是處理和存儲地理大數據的云計算平臺,已在作物生長監測、土地覆蓋變化和災害預警等領域成功應用[19]。GEE平臺實現的LandTrendr算法(https://emapr.github.io/LT-GEE)已經成為大區域長時間序列變化監測的有利工具[20]。

綜上,本研究以GEE平臺共享的Landsat TM/ETM+/OLI影像為數據源,利用LandTrendr算法提取雅魯藏布江流域NDVI時空變化特征。在系統分析1985—2018年雅魯藏布江流域NDVI變化特征基礎上,按NDVI不同的變化強度和持續時間進行分類,進而對流域內NDVI變化模式進行判別,并進一步探查不同變化模式產生的可能原因,以期為雅魯藏布江流域生態環境保護提供決策支持。

1 資料與方法

1.1 研究區概況

雅魯藏布江流域位于北緯27°50′—31°15′,東經81°57′—97°5′,發源于喜馬拉雅山北麓的杰瑪央宗冰川,流經里孜、派鄉和巴昔卡,最后流入印度,流域總面積25.8萬km2(圖1)。流域內自然環境復雜,海拔范圍在149—7172m之間,多年平均氣溫4.7—8.3℃,年降水量251.7—580.0mm[6]。復雜的自然環境造就了多樣的植被覆蓋類型,上游段植被覆蓋以高山草甸和高山草原為主,中游段主要為灌木草原,下游段以喬木和次生植被為主要類型[21]。雅魯藏布江流域因其獨特的自然環境,造就了該區域多樣的植被覆蓋類型,因其脆弱和敏感的自然環境而備受研究者關注。

圖1 研究區示意Fig.1 Map of study area

1.2 數據與處理

本研究中使用的NDVI時間序列數據由Google Earth Engine地理云計算平臺提供的1985—2018年Landsat TM/ETM+/OLI影像生成時間序列堆棧計算生成,其余數據均來源于中國科學院資源環境科學數據中心[22]。其中Landsat 5 TM影像共1737景,Landsat 7 ETM+影像共1877景,Landsat 8 OLI影像共1077景,并對數據堆棧中數據質量較差區域(云,云陰影和雪)進行掩膜。由于Landsat 7 ETM+相比Landsat 8 OLI數據具有較低的輻射分辨率,本研究利用反射率歸一化最小化傳感器之間的差異[23]。為獲取以年為周期的NDVI時間序列數據,將研究區內同一年份所有NDVI數據進行年最大值合成,以削弱大氣氣溶膠和太陽高度角等因素的影響。

1.3 研究方法

1.3.1NDVI時間序列分割算法

時間序列分割算法(LandTrendr)對以年為間隔的NDVI時間序列進行分割、擬合和平滑,獲取單個像元NDVI值在整個研究時間段內的變化特征[17](圖2)。該分割算法主要分為如下4個步驟：

圖2 LandTrendr分割算法示意Fig.2 Conceptual diagram of LandTrendr segmentation algorithm

(1)數據預處理。去除NDVI時間序列軌跡異常值,設定閾值對軌跡中的突變值進行篩選；去除超過NDVI正常變化趨勢的異常值,并對缺省值進行插值填補。為削弱異常值對變化趨勢的影響,增加信噪比,本研究采用3×3滑動窗口均值進行異常值剔除。

(2)時間序列分割。識別軌跡變化劇烈位置作為潛在分割點,并對潛在分割點進行判斷。對比潛在分割點NDVI值,剔除由云、云陰影和雪等因素造成的偽分割點,并通過角度閾值剔除小于閾值的分割點。

(3)分割點擬合。對比第一和第二分割點間連線和擬合直線與原始NDVI軌跡間的均方誤差(Mean Square Error,MSE),選取誤差最小的擬合方法作為初始連接線,并在保證曲線連續的情況下按照MSE最小原則對NDVI時間序列進行擬合。

(4)擬合模型簡化與遴選。通過迭代最大分割時間序列,利用定義的恢復率和置信度對模型進行篩選。考慮到植被覆蓋恢復期應長于1年,剔除恢復率過高的分割段,使擬合模型具有更好的分割效果

本研究基于GEE平臺實現的LandTrendr算法對NDVI時間序列進行分割與擬合[24]。分割算法需要設置適用于研究區的分割和擬合參數(表1),通過多次參數調整以適應該研究區的NDVI變化特征。

表1 時間序列分割模型參數Table 1 Time series segmentation model parameters

1.3.2NDVI變化概念模型

全球氣候變化和人類活動對雅魯藏布江流域脆弱的生態環境產生顯著影響。溫度和水分在整個生長季內對植被生長產生顯著影響,水分控制作用是影響青藏高原植被生長的最主要因素[25],城市化過程加劇和過度放牧等人類活動促使該區域NDVI下降,一定程度的氣溫和降水量增加以及退耕還林還草等生態保護工程都對NDVI上升產生促進作用[26]。基于前人成果,本研究構建雅魯藏布江流域NDVI變化概念模型(圖3),即,全球暖濕化和生物質累積使年時間分辨率NDVI曲線隨植被生長緩慢上升,人類活動加劇和自然環境惡化導致NDVI曲線下降,生態保護工程或氣候暖濕化等自然環境改善又使NDVI曲線呈現上升趨勢,生態保護工程效果減弱或者保護區域再次被破壞都將導致NDVI曲線的再次下降。影響NDVI曲線趨勢的因素包括植被物候特征、人類活動和遙感影像采集的時間和質量,因此本研究將NDVI曲線波動最劇烈年份作為變化時間點,基于持續時間和變化強度對NDVI變化模式進行判別。

圖3 NDVI變化概念模型Fig.3 Conceptual model for NDVI change

2 結果

2.1 干擾與恢復時空分布特征

2.1.1空間分布

利用時間序列分割算法對雅魯藏布江流域1985—2018年NDVI時間序列數據進行變化特征提取,得出NDVI干擾與恢復空間分布情況(圖4)。結果顯示,自流域上游至下游發生干擾年份逐漸推遲的趨勢明顯,其空間異質性與地形特征存在明顯的相關性。流域上游(貢嘎縣以上)的廣大區域NDVI干擾集中發生于2000年左右。流域中上游的貢嘎縣、昂仁縣和日喀則市NDVI干擾集中發生在1986—1990年間。中游當雄縣等海拔較高區域NDVI干擾主要發生于在2000年之前,而拉薩河谷等海拔較低區域NDVI干擾主要發生在2010年后。流域下游工布江達縣、米林縣和林芝縣海拔較低區域NDVI干擾年份發生在2000年左右,而海拔較高區域干擾多發生在2010年之后。自流域上游至下游恢復發生年份具有明顯提前趨勢,其與海拔、地形特征關聯度較高。流域中上游的廣大區域NDVI恢復年份發生在2010年之后。流域下游的米林縣和林芝縣NDVI恢復發生年份存在明顯的海拔差異,高海拔區域恢復年份發生在2010年之后,低海拔區域恢復發生在2000年左右,但念青唐古拉山及唐古拉山以南的部分高海拔區域恢復發生在1990年之前,不同于以上規律。

圖4 干擾與恢復空間格局Fig.4 Spatial pattern of disturbance and recovery

2.1.2時間分布

NDVI干擾與恢復面積隨時間變化特征能表征植被覆蓋變化活動的強弱,按NDVI時序變化趨勢,流域內18.2萬km2劃分為干擾類別,24.3萬km2劃分為恢復類別(圖5)。流域不同時間段植被覆蓋變化強弱不同。1986—1990年,干擾和恢復面積皆呈先上升后下降趨勢,干擾區域占流域總面積20.97%,恢復區域占總面積43.49%,恢復面積大于干擾二者相差22.52%,流域植被覆蓋總體呈恢復態勢。1991—2000年,干擾和恢復面積也呈現先上升后下降趨勢,分別占總面積30.54%和27.67%,干擾面積略大于恢復,二者相差2.87%,研究區植被覆蓋總體呈下降趨勢。2001—2017年期間,干擾和恢復都呈現波動變化趨勢,變化比例分別占流域總面積19.03%和23.02%,恢復面積略大于干擾,二者相差3.99%,流域植被覆蓋度總體呈現恢復態勢。綜觀1986—2017年間,研究區NDVI干擾和恢復植被變化面積總體呈先升后降趨勢,恢復區域面積大于干擾區域面積,二者相差23.64%。

圖5 1985—2018年干擾與恢復面積年際變化 Fig.5 Annual variation of disturbance and recovery area from 1985 to 2018

2.2 干擾和恢復持續時間與強度

2.2.1持續時間

1985—2018年間雅魯藏布江流域NDVI干擾和恢復持續時間具有明顯的時間序列非平穩性,干擾和恢復持續時間皆呈現先下降后上升的趨勢,分別集中在1—10年和28—32年,其中二者下降持續期(1—10年)干擾略小于恢復,而上升持續期(28—32年)干擾明顯大于恢復(圖6)。NDVI干擾持續時間85%少于10年,干擾平均持續時間4.96年,其中49.08%的干擾持續時間少于等于1年。流域內NDVI干擾持續時間超過20年的區域主要集中在拉薩城區、墨竹工卡北部和仲巴縣以南的部分區域(圖7)。NDVI恢復持續時間50%少于6年,平均恢復時間12.55年。對比地形數據可見,NDVI恢復時間超過20年的區域主要分布在河谷等海拔相對降低的區域(圖7)。

圖6 干擾與恢復面積隨持續時間和強度變化Fig.6 The duration and magnitude of disturbance and recovery area

2.2.2強度

NDVI干擾和恢復強度主要集中在0—0.75之間,恢復強度大于干擾強度(圖6)。干擾活動劇烈區域主要集中在拉薩城區、林芝及其北部區域和里加(圖7)。NDVI干擾強度95%都集中在0—0.42之間,平均下降0.14。NDVI干擾強度接近2時表明由植被轉為水體。NDVI恢復強度95%集中在0—0.4之間,平均上升0.13,流域上游和中游區域NDVI上升強度較低,NDVI上升劇烈區域主要集中在里加和林芝以東區域(圖7)。

圖7 干擾與恢復持續時間和強度空間格局Fig.7 Spatial pattern of duration and magnitude of disturbance and recovery

2.3 干擾與恢復變化模式

1985—2018年NDVI時間序列為判別NDVI變化模式提供基礎,將1987、2000和2016三年NDVI擬合數據進行假彩色合成,進而判別NDVI的不同變化模式。流域內NDVI變化主要分為上升和下降兩種趨勢,NDVI變化模式空間分布如圖8所示,NDVI上升在整個流域內都有分布,而NDVI下降主要集中在河流谷地等海拔較低的區域,共劃分為12種典型變化模式(圖8),其中具有明顯變化模式的面積占流域總面積87.23%。

流域內NDVI總體呈現上升趨勢,其中上升模式E和下降模式E所占面積比例最高(圖8)。下降模式E具有NDVI在1987—2000年間逐漸下降,2000—2016年NDVI快速下降的變化模式,主要集中在拉薩河和米林寬谷；下降模式F次之,表征NDVI在研究時間段內降低減速,主要分布在南木林和墨竹工卡等區域。上升模式E為流域內主要的NDVI增加變化模式,NDVI在1987—2000年快速上升,2000—2016年緩慢上升,主要集中在拉薩市北部和林芝縣以東海拔較高的山區；先上升后下降的上升模式A主要分布在林芝和日喀則以東區域。

圖8 雅魯藏布江流域NDVI變化模式Fig.8 Different patterns of NDVI change in the Yarlung Zangbo River Basin

3 討論與結論

3.1 討論

本研究選擇1985—2018年經過地形校正的Landsat TM/ETM+/OLI Level1數據,通過掩膜剔除云、云陰影和雪的區域,計算時間序列數據集的NDVI值,將一年內所有的NDVI數據進行最大值合成,相比較月或更短時間間隔的NDVI時間序列數據,該方法既能突出植被覆蓋年際變化特征,又能削弱云、大氣和太陽高度角等因素的干擾[13]。時間序列分割算法(LandTrendr)對NDVI時間序列進行分割和擬合,提取NDVI變化特征,廣泛應用于植被覆蓋變化及其變化原因分析[27]。相比于基于統計特征的時間序列變化檢測方法,時間序列分割方法能提取子序列的時間變化特征,進而分析時間序列的變化模式[28]。

雅魯藏布江流域植被對全球氣候變化的響應顯著,氣候變暖和降水量增加導致植被覆蓋增加,NDVI總體呈上升趨勢,但受種間關系和降水時空分布的不均勻性,以及人類活動的影響,不同區域具有不同的變化趨勢,這與已有學者研究結論相一致[29-30]。對比圖4與圖8, 念青唐古拉山及唐古拉山以南的部分高海拔區域恢復在1990年前的區域NDVI變化模式不顯著,考慮到NDVI變化的復雜性,尚需進一步研究。

20世紀80和90年代,受全球氣候變化的影響,流域上游和中游地區高原高寒植被物候期延長趨勢明顯,凈初級生產力與得到顯著提高[31],與研究結果上游和中游NDVI恢復持續時間較短且恢復持續強度較低相對應。上游和中游部分海拔較高區域NDVI干擾活動劇烈,可能與高海拔地區降水和長波輻射等因素對NDVI持續上升的限制因素有關[2]。下游森林覆蓋區受采伐和保護雙重作用,致使當地植被覆蓋快速變化,如20世紀90年代前的森林大規模采伐、1998年后天然林資源保護工程實施[32],森林破壞面積下降,恢復面積緩慢增加,流域內恢復面積大于干擾破壞面積,NDVI干擾和恢復都經歷了劇烈變化,但NDVI總體呈現增加趨[9],與本研究中恢復和干擾在2000年前變化劇烈,而2001年后變化活動逐漸減弱的結果相吻合。影響流域內植被恢復強度受不同植被覆蓋類型存在差異,流域上至中下游存在明顯的恢復強度差異。本文研究結果表明,干擾和恢復持續時間都出現了下降后上升的趨勢,干擾和恢復持續時間主要集中在10年內的較短時期,分布集中在流域的中上游,這與近20年來青藏高原暖濕化,草地生態系統恢復期較短相對應,下游區域因森林經歷大規模采伐,森林覆蓋需較長的時間才能恢復。

NDVI增加模式以先加速上升后減速為主,主要位于雅魯藏布江流域下游植被覆蓋較好的區域,植被覆蓋增加受氣候和人類活動的共同影響,NDVI增加速率降低可能與降水的空間分布與樹種組成結構變化有關[30]。隨著生態保護工程的啟動30多年[33],老年林所占的比例有所增加NDVI增加趨勢逐漸減弱。NDVI下降速率逐漸增加主要集中在人類活動較多的寬谷地區,放牧和城市化過程加劇導致的植被破壞,NDVI下降模式在人口密集的海拔較低區域出現空間聚集[34],這與已有的研究結論相一致。考慮到影響雅魯藏布江流域植被變化的因素眾多,不同影響因素在不同的區域時空分布差異巨大,因此在分析NDVI變化模式特征及其時空分異時,還應分別對不同影響因素間相互作用模式進行耦合研究。

3.2 結論

本研究基于1985—2018年Landsat數據合成的年最大NDVI時間序列為數據源,通過Google Earth Engine平臺實現的LandTrendr時間序列分割方法對NDVI時空變化特征進行提取,進而描述研究區內NDVI變化模式,主要結論如下：

(1)雅魯藏布江流域植被時空變化顯著,總體呈現上升趨勢,僅在局部區域出現下降現象,自上游至下游NDVI變化強度逐漸增加。1985—1990年恢復面積大于干擾,NDVI變化劇烈；1990—2000年干擾和恢復變化面積差異不大；2000—2018恢復面積略大于干擾,NDVI變化活動逐漸減弱。

(2)1985—2018年間流域內NDVI干擾95%集中在0—0.42之間,平均下降0.14,平均干擾時間4.96年。NDVI恢復95%集中在0—0.4之間,平均上升0.13,平均恢復時間12.55年。NDVI干擾和恢復在強度上相近,但平均持續時間和變化面積差異顯著,因此在近30年內NDVI總體呈現上升趨勢,但上升趨勢逐漸減弱。

(3)雅魯藏布江流域NDVI變化模式以上升為主,流域內NDVI主要的干擾變化模式為先緩慢下降后快速下降；NDVI恢復模式為先快速上升后緩慢上升。