鐘旭珍,王金亮,鄧云程,李 杰,吳瑞娟,董品亮

1 云南師范大學地理學部,昆明 650500 2 內江師范學院地理與資源科學學院,內江 641100 3 云南師范大學-西南聯(lián)合研究生院,昆明 650500 4 云南省高校資源與環(huán)境遙感重點實驗室, 昆明 650500 5 云南省地理空間信息工程技術研究中心, 昆明 650500 6 美國北德克薩斯大學地理系, 丹頓 76203

近百年來全球氣溫持續(xù)上升,根據(jù)IPCC發(fā)布的《氣候變化2021:公眾摘要》報告顯示:與19世紀末(工業(yè)革命前)相比,2011年至2020年地球表面的平均氣溫升高了1.1℃,并且比過去12.5萬年的任何時候都高[1]。歷年的IPCC發(fā)布的報告均認為:近百年來,人類過度的開荒及大規(guī)?；顒?導致所排放的溫室氣體有增無減,這是造成全球變暖最主要的原因[2]。人類活動不僅導致了全球氣候的變暖,還誘發(fā)了如極端天氣、植被退化、水土流失等生態(tài)安全問題。生態(tài)環(huán)境是人民生存的基礎條件,人類可持續(xù)發(fā)展的保障。而植被作為生態(tài)系統(tǒng)主要的組成部分,對區(qū)域生態(tài)系統(tǒng)環(huán)境變化有著重要指示,對全球碳循環(huán)和氣候調節(jié)具有非常重要的作用。植被覆蓋度(FVC)是指植被在地面的垂直投影面積占總面積的百分比,它量化了植被的茂密程度,反映了植被的生長態(tài)勢,是描述生態(tài)系統(tǒng)變化的重要基礎數(shù)據(jù)[3]。怒江-薩爾溫江是東南亞最長的自由流動河流和最重要的跨界河流之一,在中國-中南半島經(jīng)濟走廊建設和中緬經(jīng)濟走廊建設和生態(tài)保護中發(fā)揮著不可替代的作用。在過去的50年里,該流域冬季和春季氣溫呈上升趨勢,氣溫上升的幅度隨著海拔高度的升高而增加[4]。然而,目前對該流域關注不夠,尚沒看到其植被覆蓋研究報道。因此,研究分析其植被覆蓋時空演變,并探討其驅動力,可為其生態(tài)環(huán)境保護和社會經(jīng)濟發(fā)展提供重要指導作用[5]。

估算植被覆蓋度(FVC)是植被覆蓋相關研究的重要基礎,有多種估算方法,如地表實測法、經(jīng)驗模型法、植被指數(shù)法、決策樹分類法等,部分方法存在對實測數(shù)據(jù)依賴性強、計算繁瑣、需要大量樣本數(shù)據(jù)等缺點,特別對于偏遠地區(qū)來說很難實現(xiàn)[6],而遙感技術的快速發(fā)展以及遙感大數(shù)據(jù)時代的到來,使得在大時空尺度上精確監(jiān)測偏遠地區(qū)的植被行為成為可能[7],其中MODIS NDVI產(chǎn)品被大量學者用于植被覆蓋度的計算及其動態(tài)變化的定量分析,而像元二分模型法是常用的FVC估算方法,其操作較為簡單,可直接利用NDVI數(shù)據(jù)來進行計算[8-10]。目前針對植被覆蓋演變規(guī)律的研究主要是運用一些常規(guī)的線性趨勢分析方法,如線性回歸分析方法[11-12]、相關性分析法[13-14]、Theil-Sen斜率與Mann-Kendall檢驗結合的趨勢分析法[15],很少有學者關注其非線性變化特征。然而,植被覆蓋的時空分布與變化受多方面的因素影響,使得其演變規(guī)律也具有復雜性,并不僅僅呈現(xiàn)為簡單的線性變化,而BFAST法可用于分析不同類型的時間序列,對長時間序列過程中趨勢的突變檢測具有很好的應用效果,已有學者將其用于植被的非線性特征[15-17],以及氣溫變化[18]、氣溶膠變化[19]等領域研究中,并顯示了其較好的適用性[20-21]。對于南北跨度極大的怒江-薩爾溫江流域,其不同緯度區(qū)域的植被變化特征及原因也會存在差異,因此有必要運用BFAST法來更加準確地探究其植被變化的非線性特征及其未來發(fā)展趨勢。

對于植被未來可持續(xù)性的研究,Hurst指數(shù)法是廣大學者常用的方法[22-23]。耦合BFAST法與Hurst指數(shù)可以實現(xiàn)植被未來趨勢的更加精確表達。驅動植被變化的因素,有人為因素和自然因素,以往學者大多采用偏相關分析、回歸分析等方法研究植被覆蓋度變化的影響因素[24-25],這些方法假設了植被與各環(huán)境因子存在線性關系。然而植被在環(huán)境中的生長過程影響因素復雜,可能不存在嚴格的線性關系[26]。而由王勁峰等學者提出的地理探測器模型(Geographical Detector,簡稱Geodetector),不僅可用于探測空間分異性并識別其主要的影響因子,還能實現(xiàn)不同因子間的交互作用探測[27]。此方法無線性假設,可以避免多重共線性問題,在計算驅動因子的貢獻率時結果較為準確[28]。其改進的模型OPGD模型(Optimal Parameters-Based Geographical Detector,OPGD)加入了參數(shù)優(yōu)化模塊,可根據(jù)自變量的特征進行多種方式和多種斷點數(shù)目的離散化處理,選取出最優(yōu)的離散化方式和最優(yōu)的斷點數(shù)目,提高各類別間的顯著差異性[29]。因此,研究采用該方法進行植被非線性特征的驅動因子探測。

基于以上分析可知,盡管已經(jīng)開展了許多關于植被時空演變與驅動因素的研究,但目前還沒有學者采用BFAST法研究其歷史非線性趨勢的同時結合Hurst指數(shù)表達其未來更加精確的可持續(xù)性,也就是考慮了植被非線性特征的未來趨勢預測。另外,怒江-薩爾溫江是東南亞一條重要的國際河流,對南亞、東南亞的氣候調節(jié)具有重要作用,但由于其南北跨度大,涉及中國、緬甸和泰國,對其進行相關研究尤其數(shù)據(jù)收集較為困難,目前還沒有學者對該流域進行植被覆蓋的相關研究。因此,研究以該流域為研究對象,基于MODIS NDVI數(shù)據(jù),在遙感與GIS技術的支持下,采用像元二分模型估算其近22年的植被覆蓋度,并借助 BFAST01模型及 Hurst 指數(shù)等,對區(qū)域內植被覆蓋的時空分布、非線性變化趨勢和未來可持續(xù)性進行分析,并結合坡度、海拔、降雨、氣溫、人口、GDP、土地類型等自然和人類活動因子,引入OPGD對影響研究區(qū)植被覆蓋度的主要因素進行定量探究,以期為研究區(qū)的生態(tài)環(huán)境保護和可持續(xù)經(jīng)濟發(fā)展提供科學依據(jù)。

1 研究區(qū)概況及數(shù)據(jù)來源

1.1 研究區(qū)概況

怒江-薩爾溫江發(fā)源于青藏高原上的唐古拉山脈,經(jīng)中國云南(保山、臨滄)流入緬甸,最后注入印度洋的安達曼海,經(jīng)過中國、緬甸和泰國(圖1)。該流域形狀極為狹長,河流全長約3200 km,流域總面積為32.5萬km2,它在中國的部分被稱為怒江,其余的部分稱為薩爾溫江,是東南亞最長的自由流動河流,也是最重要的跨境河流之一。楊帆等的研究中,將嘉玉橋以上劃分為上游,嘉玉橋到中緬邊界為中游,中國邊界以南為下游[30]。由于該流域南北跨度和海拔高差都很大,橫跨多個氣候帶和生態(tài)區(qū),其氣候多樣,地形復雜,不同區(qū)域的地形、氣候、土地利用類型都存在很大的空間差異性,本研究根據(jù)流域海拔特征并結合行政區(qū)劃界限,將云南和西藏交界線以北區(qū)域劃分為上游,屬于青藏高原,為典型的高原山地氣候區(qū),氣候變化較為敏感;云南西藏交界線到中緬邊界劃分為中游,主要為高山峽谷區(qū),即狹窄又陡峭,氣候垂直差異明顯,流經(jīng)三江并流保護區(qū);中國邊界以南為下游,海拔相對北部較低,地勢較為平坦,屬于熱帶季風氣候,以此分區(qū)作為后期相關研究和分析的基礎。

1.2 數(shù)據(jù)源及預處理

研究數(shù)據(jù)主要包括MODIS NDVI 數(shù)據(jù),地形數(shù)據(jù),土地利用數(shù)據(jù),氣象數(shù)據(jù)以及人類活動數(shù)據(jù),研究區(qū)流域邊界和行政區(qū)劃數(shù)據(jù)等。詳述如下:

1.2.1NDVI 數(shù)據(jù)

研究使用2000-2021年MODIS NDVI數(shù)據(jù),通過GEE(Google Earth Engine,GEE)遙感大數(shù)據(jù)云平臺進行調用、預處理及下載,數(shù)據(jù)原網(wǎng)站為USGS (https://lpdaac.usgs.gov/products/),其時間分辨率為16d,空間分辨率為250m。在GEE平臺中對其進行鑲嵌、投影、轉換和裁剪等預處理,為了消除云、大氣效應、掃描角度和太陽天頂角等對影像質量的影響,研究采用了最大值合成法MVC(Maximum Value Composite)來進行NDVI數(shù)據(jù)的合成[15],該方法被廣泛用于植被動態(tài)研究的預處理中[31]。借助ENVI IDL平臺,選用像元二分模型(DPM)來計算研究區(qū)月度及年度FVC,用于后續(xù)植被覆蓋的相關分析。

1.2.2地形數(shù)據(jù)

研究使用的地形數(shù)據(jù)為2000年采集的SRTM Digital Elevation Data Version 4,分辨率為90m,數(shù)據(jù)生產(chǎn)于 NASA (https://srtm.csi.cgiar.org/),通過GEE平臺對其進行拼接、裁剪和下載,并通過ArcGIS 10.2軟件提取如海拔、坡度、坡向等地形因素作為后續(xù)植被覆蓋驅動力分析的指標。

1.2.3LULC 數(shù)據(jù)

研究使用的土地利用數(shù)據(jù)為ESA WorldCover 10m v100,來自esa (https://esa-worldcover.org/en),也是通過GEE平臺對其進行投影轉換和裁剪等預處理和下載。該產(chǎn)品將土地利用類型分為了11類,比較詳細的表示了不同地表覆蓋類型,我們將其作為植被覆蓋的人類活動驅動因子之一。

1.2.4氣象數(shù)據(jù)

研究使用的氣象數(shù)據(jù)主要包括2000-2021年的氣溫、降水、蒸散發(fā)以及氣候類型。其中氣溫數(shù)據(jù)由ERA5 提供的第五代 ECMWF 全球氣候大氣再分析數(shù)據(jù),通過GEE平臺進行預處理與下載。我們選用的是mean_2m_air_temperature,單位為K;降水量數(shù)據(jù)選用的是Climate Hazards Group InfraRed Precipitation With Station Data 2.0 (https://chc.ucsb.edu/data/chirps), 單位為mm,分辨率為0.05°。蒸散發(fā)數(shù)據(jù)為GLEAM ET產(chǎn)品,空間分辨率為0.25°,下載于gleam (https://www.gleam.eu/#downloads)。蒸散發(fā)可以反映植被的干旱情況,對植被生長狀況有一定影響,因此研究選取該指標作為驅動因子之一。氣候類型數(shù)據(jù)為全球柯本氣候類型空間分布數(shù)據(jù)集,空間分辨率為0.1°,下載于國家地球系統(tǒng)科學數(shù)據(jù)中心(http://www.geodata.cn/),將其作為植被覆蓋驅動力分析因子之一。

1.2.5人類活動等其他數(shù)據(jù)

人口格網(wǎng)數(shù)據(jù)來源于world pop (https://hub.worldpop.org/geodata/listing?id=64)。單位為people/km2,空間分辨率為1km。全球GDP數(shù)據(jù)來源于GitHub (https://github.com/Nowosad/global_population_and_gdp), 分辨率為0.5° × 0.5°。中國區(qū)域人口、GDP數(shù)據(jù)來源于中國科學院資源環(huán)境科學數(shù)據(jù)中心(http://www.resdc.cn)。其他數(shù)據(jù)如研究區(qū)行政區(qū)劃,流域邊界等分別來源于中國科學院資源環(huán)境科學數(shù)據(jù)中心 (https://www.resdc.cn/data.aspx?DATAID=259) 和HydroSHEDS 網(wǎng)站 (https://www.hydrosheds.org/)。數(shù)據(jù)來源如表1,驅動因子指標如圖2。

表1 研究數(shù)據(jù)來源與描述Table 1 Source and description of research data

圖2 FVC驅動因子空間分布圖Fig.2 Spatial distribution of FVC driving factorsET:蒸散發(fā) Evapotranspiration;FVC:植被覆蓋度 Fractional vegetation cover

2 研究方法

2.1 像元二分模型

像元二分模型是一種簡單實用的遙感估算模型,對綠度信息具有很好的敏感性,且對不同氣候區(qū)、植被類型等的植被信息具有很好的識別效果,被廣泛用于計算植被覆蓋度。其計算植被覆蓋度的公式如下:

(1)

式中,NDVIsoil和NDVIveg表示全部為裸土和全部為植被的像元NDVI值。理論上,NDVIsoil為0,NDVIveg為1。然而,實際中植被覆蓋的影響因素復雜,NDVIsoil和NDVIveg的值并不等于0和1,而是受圖像本身質量的影響。參考相關研究[23,32],根據(jù)MODIS NDVI時間序列圖像,對NDVI圖像進行頻率直方圖統(tǒng)計,去除異常值以及歸一化處理,并取累計頻率在2%和98%時的NDVI值分別作為NDVIsoil和NDVIveg參數(shù),然后按照公式(1)計算研究區(qū)FVC。

2.2 Hurst指數(shù)

Hurst 指數(shù)是定量描述時間序列長程依賴性的有效方法[33]。近年來在植被覆蓋變化的未來可持續(xù)性研究中得到廣泛應用[34-35]。研究利用matlab計算出H值即Hurst指數(shù),來判斷植被覆蓋時間序列是否存在持續(xù)性。H值介于0-1間,參考相關研究成果[23,36],其分類有三種:0

2.3 BFAST模型

BFAST是一種迭代算法,它利用分段線性趨勢和季節(jié)性模型將時間序列分解為趨勢、季節(jié)性和殘差組分,并能檢測趨勢和季節(jié)性成分的突變[7,37]。BFAST假設非線性可以用分段線性模型來近似,因此,當一個序列可以更好地用多個線性段表示,而不是單一的單調趨勢時,就可以使用結構變化測試確定一個斷點并確定日期[26]。相比其他突變檢測方法,在時間序列內,BFAST受季節(jié)差異和噪聲的影響較小,可以更快地檢測到變化[19]。Verbesselt等對該方法進行了比較完整的描述和介紹[38-39]。在數(shù)學上,BFAST的趨勢分量和季節(jié)分量通過以下分解模型得到:

Yt=Tt+St+et,t=1,…,n

(2)

式中,Yt為t時觀測到的值,Tt為趨勢組分,St為季節(jié)性組分,et為殘差組分。其中趨勢組分可用下式表示:

Tt=ai+bit(t=1,2,…,n)

(3)

季節(jié)性成分可用下式表示:

(4)

式中,ai和bi為趨勢項系數(shù),γi為振幅,δi為分段數(shù)量,f為頻率,其中振幅、分段數(shù)量為未知項,而頻率為已知項。

季節(jié)組分和趨勢組分中突變點的識別需要確定突變點的數(shù)量和突變點在時間序列中的位置。此處利用最小二乘移動求和從季節(jié)組分和趨勢組分中檢驗判斷是否存在突變點[40]。首先,采用普通最小二乘方法(Ordinary Least Squares, OLS)在一個不斷增加的移動窗口內進行檢驗判斷是否存在斷點,如果根據(jù)統(tǒng)計檢驗,窗口內的普通最小二乘殘差的移動和很大,則被描述為一個斷點。接著,通過最小化所有窗口上斷點之間的最小二乘移動總和(OLS-MOSUM),以此來確定斷點的位置[19,38]。

根據(jù)BFAST原理,本研究采用了BFAST的改進模型BFAST01。BFAST01模型相比BFAST模型的優(yōu)勢在于,它同時考慮了季節(jié)模型和趨勢模型,只檢測時間序列中最具影響的趨勢變化,而不是多個較小的變化,也就是說BFAST01只檢測0個或者1個斷點,其為評估長期季節(jié)性和年際變化的主要轉折點提供了一種有效的、全面的變化檢測方法[26,41]。本研究利用R環(huán)境,從GitHub加載最新的bfast包(https://github.com/),完成所有BFAST01的統(tǒng)計分析[42],提取出植被覆蓋的趨勢顯著性、變化強度、趨勢斷點發(fā)生的時間和次數(shù)、突變類型、突變年份等。參考相關研究[15,19,26],將BFAST01檢測到的植被覆蓋變化趨勢類別分為8類,如表2所示。

表2 BFAST01檢測的FVC變化趨勢類型Table 2 Types of FVC change trends detected by BFAST01

2.4 地理探測器

地理探測器Geodetector,是一種統(tǒng)計工具,可以探測地理現(xiàn)象的空間分異性并揭示其背后驅動力,既適用于點數(shù)據(jù),也適用于面數(shù)據(jù),對于面數(shù)據(jù)一般要先進行重分類和離散化處理[27]。其包括4個探測器,我們選用因子探測器和交互探測器來探測植被覆蓋變化的影響因素。

因子探測器,用于探測因變量Y的空間分異性,及某個因素X對因變量Y空間分布的解釋力,其大小可用q值來度量(0≤q≤1),本研究中,若某影響因子的q值越接近于1,說明該因子對FVC的解釋力越強,反之越弱[27,43-44]。

交互探測器可以用于識別不同因子之間交互對植被覆蓋度變化的影響是增強還是減弱,或是獨立[27]。其交互方式見表3。

表3 地理探測器交互作用方式Table 3 Interaction modes of geodetectors

3 結果與分析

3.1 怒江-薩爾溫江流域植被覆蓋時空特征分析

3.1.1空間分布特征

通過像元二分模型計算得到FVC,并利用GIS空間分析可視化工具得到研究區(qū)2000-2021年平均FVC空間分布圖(圖3),對其進行分析可知,怒江-薩爾溫江流域低植被覆蓋度(0-0.2)、較低植被覆蓋度(0.2-0.4)、中植被覆蓋度(0.4-0.6)、較高植被覆蓋度(0.6-0.8)、高植被覆蓋度(0.8-1)占比依次為:5.41%,5.66%,10.41%,21.46%,57.06%,FVC分布具有明顯的空間異質性,從圖3可以看出,下游和中游地區(qū)植被覆蓋情況明顯優(yōu)于上游地區(qū),總體表現(xiàn)為南高北低。整個流域平均FVC值為0.73,高植被覆蓋度和較高植被覆蓋度地區(qū)面積占比達到78.52%,說明該流域植被覆蓋情況較好。結合研究區(qū)土地利用類型分析發(fā)現(xiàn),低植被覆蓋區(qū)域主要分布在海拔較高的青藏高原,主要土地利用類型為苔蘚、裸地、灌木等;高植被覆蓋區(qū)域主要在中游和下游地區(qū),土地利用類型以林地、草地和耕地為主。

圖3 怒江-薩爾溫江流域2000-2021年平均FVC空間分布圖Fig.3 Spatial distribution of the average FVC in the Nujiang-Salween River Basin from 2000 to 2021

3.1.2時間變化特征

通過統(tǒng)計研究區(qū)年、月、季節(jié)的FVC均值,并制作FVC變化折線圖(圖4)。從圖4的年際變化圖可以看出,怒江-薩爾溫江流域2000-2021年植被覆蓋整體呈波動上升趨勢,FVC最大值出現(xiàn)在2021年為0.750,最小值出現(xiàn)在2015年為0.722,在2007年和2015年出現(xiàn)了明顯的低值拐點,但2016年以來,植被生長狀況出現(xiàn)明顯上升趨勢。從研究區(qū)連續(xù)22年植被覆蓋季節(jié)變化圖可以看出怒江-薩爾溫江流域植被變化具有明顯的季節(jié)規(guī)律性,這與研究區(qū)氣候有關。而FVC月際變化圖表明,3月份FVC值最低,3-9月份FVC值隨時間逐漸上升,在9月份的時候FVC值達到最高,隨后呈下降趨勢;整體上夏秋季植被覆蓋狀況較春冬季植被覆蓋狀況好,不同季節(jié)FVC隨著時間變化呈現(xiàn)較小的波動,春季和夏季波動較秋季和冬季大,但整體較為穩(wěn)定。

圖4 怒江-薩爾溫江流域FVC時間變化Fig.4 Time variation of FVC in the Nujiang-Salween River BasinFVC:植被覆蓋度Fractional vegetation cover

3.2 怒江-薩爾溫江流域植被覆蓋非線性趨勢特征

運用BFAST01算法對怒江-薩爾溫江流域2000-2021年FVC進行突變檢測,得到FVC的非線性趨勢特征(圖5)。從圖中可知,研究區(qū)分為了8種非線性植被覆蓋趨勢類型,其中單調遞增占34.63%,是占比最多的,在流域上、中、下游各區(qū)域均有分布;單調遞減類型占比12.22%,主要分布在流域下游區(qū)域;“單調遞增+”占比2.3%,主要分布在中緬邊界;“單調遞減-”占比最少為0.63%;“中斷-+”占比22.91%,排名第二,主要分布在流域最北端和云南以及緬甸北部;“中斷+-”占比8.19%,主要分布在流域最南角,云南和北部也有些許分布;“反轉+-”占比7.73%,“反轉-+”占比11.39%;通過分析發(fā)現(xiàn),總體呈增加趨勢類型的占比(71.24%)大于總體呈減少趨勢類型的占比(28.76%),表明研究區(qū)植被發(fā)展趨勢向好。從顯著性檢測結果來看,變化呈現(xiàn)顯著的占比為44.29%,不顯著的占比為55.71%;突變時間檢測表明,發(fā)生突變的時間主要集中在2003-2018年,2003年和2018年發(fā)生突變的數(shù)量較多,其他年份較為均勻,“中斷-+”發(fā)生的突變最多,發(fā)生突變占比較少的類型是:“單調遞增+”和“單調遞減-”。需要提醒的是,單調遞增和單調遞減沒有發(fā)生斷點,因此圖5(突變時間)中空白區(qū)域表示沒有發(fā)生突變的區(qū)域。

圖5 BFAST01檢測結果Fig.5 BFAST01 detection results

為了更加詳細的描述植被覆蓋突變情況,研究還對FVC變化的最大突變強度、最大突變強度發(fā)生的時間、首次季節(jié)性斷點發(fā)生的時間以及斷點發(fā)生次數(shù)等信息進行了檢測(圖6)。通過分析發(fā)現(xiàn)整個流域突變發(fā)生最大強度的地區(qū)主要在云南以及青藏高原北部,最大突變強度范圍介于-1.2078-1.0287之間,最大突變強度時間主要在2002-2019年,2006年發(fā)生最大強度突變的比例是最大的,2013年、2015年和2018年最大突變強度數(shù)量也較多,上游地區(qū)最大突變強度斷點發(fā)生時間變化趨勢與此相似,說明上游地區(qū)在全流域FVC的變化趨勢中具有重要作用;而中游和下游區(qū)域最大突變強度斷點發(fā)生時間變化趨勢較上游地區(qū)波動平穩(wěn),中游地區(qū)突變數(shù)量較多的年份為2006年、2010年和2015年,下游地區(qū)在2011年最大強度突變數(shù)量最多,占下游地區(qū)突變總面積的9.90%。各區(qū)域最大強度突變面積從2018年開始大面積減少。此外,根據(jù)研究區(qū)FVC時間變化可知,其植被變化具有很強的季節(jié)性,從圖6可以看出,首次發(fā)生季節(jié)性斷點的時間主要集中在2010-2014年。

此外,怒江-薩爾溫江流域最大突變強度斷點發(fā)生突變的次數(shù)分別有1次、2次和3次(圖7),表示某點在研究時間段內植被發(fā)生非線性突變次數(shù)的空間分布,1次的占比最多為77.02%,3次的最少為4.18%,2次的區(qū)域主要分布在中緬邊界和最北部。上游地區(qū)發(fā)生1次斷點的比重最大,占上游地區(qū)發(fā)生斷點次數(shù)總面積的85.72%,而發(fā)生3次斷點的占比僅占本區(qū)斷點總面積的0.91%;下游和中游地區(qū)各斷點次數(shù)占本區(qū)域斷點總次數(shù)面積比重相差不大,也是以1次為主,占比略小于上游地區(qū),而發(fā)生3次斷點的次數(shù)占比均在7%左右,略大于上游地區(qū)。

圖7 植被覆蓋突變頻率分布Fig.7 Frequency distribution of vegetation cover mutations

3.3 怒江-薩爾溫江流域植被未來可持續(xù)性

為了表示植被覆蓋的未來可持續(xù)性,我們計算了怒江-薩爾溫江流域FVC 的 Hurst 指數(shù),通過統(tǒng)計分析發(fā)現(xiàn),Hurst 指數(shù)小于0.5的區(qū)域占總面積的2.76%,這些區(qū)域具有反持續(xù)性,表示過去增加的趨勢將來總體上減少,反之亦然;Hurst 指數(shù)大于0.5的區(qū)域占比為 94.89%,表明研究區(qū)植被覆蓋的正向持續(xù)性較強;Hurst 指數(shù)為0.5的區(qū)域占2.35%,表明未來植被覆蓋發(fā)展趨勢不明確。

為了進一步揭示植被覆蓋的非線性變化趨勢及其持續(xù)性,將BFAST變化趨勢結果與 Hurst 指數(shù)結果進行疊加,得到非線性變化趨勢與持續(xù)性的耦合信息(圖8和表4)。耦合結果總共為17種情形,其中“持續(xù)性&單調遞增”耦合類型含義表示未來植被與過去發(fā)展趨勢一樣,將持續(xù)增加;“反持續(xù)性&單調遞增”表示未來植被與過去發(fā)展趨勢相反,將從增加變?yōu)闇p少,其它情形解釋相似,比較詳細的展示了植被覆蓋的未來發(fā)展趨勢。未來將呈持續(xù)改善情景的面積占比達到68.99%,持續(xù)改善的趨勢類型分別為:“持續(xù)性&反轉-+”、“持續(xù)性&中斷-+”、“持續(xù)性&單調增加+”、“持續(xù)性&單調遞增”、“反持續(xù)性&單調遞減”、“反持續(xù)性&單調遞減-”、“反持續(xù)性&中斷+-”和“反持續(xù)性&反轉+-”,這些區(qū)域占據(jù)了流域大部分地區(qū);未來發(fā)展將呈持續(xù)退化的情景所占比例為29.09%,這些趨勢類型分別為:“反持續(xù)性&反轉-+”、“反持續(xù)性&中斷-+”、“反持續(xù)性&單調遞增+”、“反持續(xù)性&單調遞增”、“持續(xù)性&單調遞減”、“持續(xù)性&單調遞減-”、“持續(xù)性&中斷+-”和“持續(xù)性&反轉+-”,這些區(qū)域主要分布在流域中游和下游;剩下1.92%的區(qū)域未來變化趨勢無法確定,主要分布在流域上游區(qū)域。持續(xù)退化和未來變化趨勢無法確定的區(qū)域的植被變化狀況需要研究人員繼續(xù)關注。

表4 BFAST01和Hurst指數(shù)耦合類型統(tǒng)計Table 4 Statistics of BFAST01 and Hurst coupling types

圖8 植被覆蓋未來發(fā)展趨勢Fig.8 Future development trend of vegetation coverage

3.4 怒江-薩爾溫江流域植被覆蓋驅動因子探測

研究選取了海拔、坡度、坡向、蒸散發(fā)、氣溫、降水、氣候區(qū)劃等自然因素和GDP、人口、土地利用等人為因素作為植被覆蓋的影響因子,分別為X1-X10,植被覆蓋度FVC為Y,作為地理探測器的數(shù)據(jù)輸入,空間分辨率統(tǒng)一重采樣為250m。由于這些因子及FVC都是連續(xù)變量,而地理探測器的輸入數(shù)據(jù)要求為類別數(shù)據(jù),則根據(jù)王勁峰等提出的數(shù)據(jù)離散化方法[27],利用OPGD在R語言環(huán)境中進行參數(shù)優(yōu)化,為了使結果更加符合實際,研究同時對離散化進行了人為調試,最終得到比較理想的結果,具體離散化方法和類別見表5。并利用ArcGIS創(chuàng)建漁網(wǎng),通過多次調試設置采樣間隔為3km,運用Spatial Analyst-提取分析-采樣工具,輸入柵格為所有重分類的X和Y,采樣點為研究區(qū)漁網(wǎng)點,將自變量和因變量的值提取到點,作為地理探測器輸入數(shù)據(jù)。

表5 地理探測器因子離散化方法與類別Table 5 Factor discretization methods and categories of geographic detectors

利用因子探測和交互探測器,得到全流域以及各子流域植被覆蓋驅動力探測結果(表6和圖9)。除中游地區(qū)的坡向因子p值為0.0791外,其余因子p值均為0,通過顯著性檢驗。根據(jù)表6可知,全流域中海拔的q值是最大的為0.5761,接著影響較大的是氣溫、降水、蒸散發(fā)等氣象因子,其次是土地利用類型、GDP、人口等人為因子,人為因子中以土地利用影響最大。各子流域因子探測結果顯示同一個影響因子在不同區(qū)域對植被覆蓋的影響具有差異性,縱向上分析,如海拔和土地利用因子對上游和中游區(qū)域的影響大于下游區(qū)域;而人口、GDP、蒸散發(fā)等因子對下游區(qū)域的影響則大于上游和中游區(qū)域。橫向上分析,上游區(qū)域影響最大的是土地利用因子,其次是海拔,而人口、GDP等人類活動因子對于較為偏遠的上游地區(qū)影響較小。對于中游區(qū)域,對植被影響最大的是海拔,其次是土地利用,氣溫、降水、蒸散發(fā)等也是影響中游區(qū)域植被覆蓋變化的重要因素,自然因素對中游區(qū)域植被的影響大于人為因素。與中游區(qū)域相反,下游區(qū)域的人口和GDP是對植被影響最大的因素,其次是蒸散發(fā),土地利用和降水影響也較大,而影響最低的是海拔。

表6 怒江-薩爾溫江流域各區(qū)域影響因子q值Table 6 Impact factor q values in the Nujiang Salween River Basin

圖9 怒江-薩爾溫江流域各區(qū)域影響因子交互作用探測結果Fig.9 Detection results of the interaction of influencing factors in the Nujiang-Salween River Basin

交互探測結果顯示(圖9),任何兩個因子交互影響都會增強其影響力,對于全流域來說,發(fā)現(xiàn)其它每個因子和海拔交互都達到最大影響值,說明海拔對植被覆蓋有著很強的影響力,這與整個研究區(qū)海拔差異顯著具有關聯(lián)性。上游區(qū)域和中游區(qū)域與全區(qū)相似,各因子與海拔的交互作用對植被覆蓋的影響增強作用是最顯著的。對于下游區(qū)域來說,各因子交互作用具有差異性,其中氣溫、降水、蒸散發(fā)與坡度的交互作用對植被的影響大于它們與其它因子的交互作用;氣候類型、土地利用、人口與蒸散發(fā)的交互作用對植被的影響大于它們與其它因子的交互作用;說明坡度和蒸散發(fā)對下游區(qū)域的植被覆蓋影響具有很大的可利用空間;而GDP與人口的交互作用也大于GDP與其它任何因子的交互作用,成為所有交互作用中影響最大的,這與因子探測的結果具有一致性。

4 討論

4.1 怒江-薩爾溫江流域植被覆蓋時空總體特征

怒江-薩爾溫江流域區(qū)位特殊,上、中、下游地區(qū)地形海拔和緯度等差異明顯,植被覆蓋空間分布特征也形成了鮮明的對比。從全區(qū)域來看,高植被覆蓋區(qū)域占比大于50%,而其中上游地區(qū)高植被覆蓋區(qū)域占比僅占7%左右,平均FVC值為0.53;中游和下游區(qū)域高植被覆蓋區(qū)域占比達到85%和90%以上,平均FVC值分別為0.86和0.88,說明中游和下游地區(qū)拉高了整個區(qū)域的植被覆蓋水平,這與研究區(qū)地形地貌有關[45-46],這將在驅動力分析中進行探討。

從FVC時間演變特征來看,怒江-薩爾溫江流域植被覆蓋在2015年以前增長趨勢不明顯,但2016年開始發(fā)展趨勢向好。這與李杰、王春雅等的研究結論相似[23,47]。但在2007年、2015年和2019年出現(xiàn)了明顯的低值拐點,結合圖4與圖10,發(fā)現(xiàn)上游地區(qū)FVC年際變化也在2007年和2015年發(fā)生了低值拐點,變化趨勢與全區(qū)FVC非常相似;而下游和中游區(qū)域FVC年際變化增長趨勢明顯,其中下游地區(qū)在2004、2012、2018等年份出現(xiàn)低值拐點,查詢資料發(fā)現(xiàn),2004年中緬邊境發(fā)生了一場大森林火災,2006年的颶風“馬拉”和2007年的臺風“利奇馬”,2012年和2013年緬甸發(fā)生了嚴重干旱,2014年、2015年的厄爾尼諾現(xiàn)象、極端干旱事件等,此外,據(jù)聯(lián)合國人類事務協(xié)調辦公室稱,緬甸在2015年經(jīng)歷了幾十年來最嚴重的洪水,同時2015年“4·25”西藏也發(fā)生了地震災害,以及2018、2019年的極端氣候自然災害如印度西南部的特大洪災等[48],這些事件均嚴重威脅到當?shù)刂脖坏纳L狀態(tài),可能是導致FVC出現(xiàn)低值拐點的原因。對于下游地區(qū)來說,有研究表明森林火災、農業(yè)轉型、公路建設增加等是導致其植被減少的主要原因[49]。而研究區(qū)整體植被變化趨勢向好與生態(tài)工程的實施和圍欄禁牧政策以及其他人類活動等有關[47]。

圖10 怒江-薩爾溫江流域各流域FVC年際變化Fig.10 Interannual variation of FVC in each sub-basin of the Nujiang-Salween River Basin

4.2 怒江-薩爾溫江流域植被覆蓋非線性特征

從研究區(qū)FVC年、月、季節(jié)的變化分析以及相關研究表明,植被變化往往受多種因素的影響使其很難呈現(xiàn)嚴格的線性變化[15]。因此,我們在進行研究時應該考慮其趨勢的非線性特征及其影響因素,特別是大尺度長時間序列的演變分析。研究運用BFAST01方法來分析FVC的演變趨勢以及非線性突變,BFAST方法可以很好地適用于突變檢測[50]。在3.2小節(jié)中,我們詳細分析了怒江-薩爾溫江流域FVC變化的突變類型、突變時間、突變顯著性與突變次數(shù)等特征,將研究區(qū)植被變化詳細分為了8種突變類型,每種突變類型的含義見表2。為了進一步分析各突變類型空間分布特征及驗證分類的準確性,我們結合 Google Earth time-series images來進行突變前后年份影像的比對(圖11)。

圖11 突變年前后FVC及Google Earth imagesFig.11 FVC and Google Earth images before and after the mutation year圖中圓圈4、5、6、8分別代表對應的突變類型;黑色矩形框所指圖像代表突變年份的Google Earth images和FVC圖像;紅色矩形框圖像代表突變年份前后的Google Earth images和FVC圖像

從上、中、下游區(qū)域分別選擇了某局部區(qū)域進行制圖比對,上游區(qū)域選擇的代表像元突變類型為 “單遞減-”,突變時間為2017年,從Google Earth time-series images上我們可以明顯看到該區(qū)域范圍突變時間前后土地利用從草地和稀疏植被逐漸轉變?yōu)榻ㄔO用地,FVC圖像也可以看出其顏色逐漸變深,也就是FVC值在減小,這就說明了該區(qū)域植被覆蓋在減少,這和突變類型“單調遞減-”的特征是相符的。中游區(qū)域選擇的代表像元突變類型為 “中斷+-”,該突變類型有1個明顯突變點,且斷點處值突然增大,即顯著正中斷,然后顯著減小,從Google Earth time-series images和FVC圖像比對來看,也符合其突變特征。下游區(qū)域選擇的兩個突變類型代表分別為“中斷-+”和“反轉-+”,“中斷-+”特征與“中斷+- 相反,從圖像特征上可以看出該局部區(qū)域植被覆蓋在斷點2012年處突然減小,即負中斷,然后增加;而“反轉-+”表示趨勢從顯著減小轉換為顯著增加,即在突變點處FVC顯著減小,也符合Google Earth time-series images和FVC圖像的比對特征。因此,通過圖像比對驗證,發(fā)現(xiàn)用BFAST模型進行長時間序列植被非線性趨勢檢測是適用的,結果也具有一定可靠性,研究可為相關生態(tài)環(huán)境保護措施的制定提供科學參考依據(jù)。

4.3 怒江-薩爾溫江流域植被覆蓋驅動力

為了更有針對性地分析植被覆蓋的影響因素,考慮研究區(qū)地形地貌的差異性,研究選擇了適用于具有非線性特征分析的基于最優(yōu)參數(shù)的地理探測器(OPGD),對研究區(qū)植被覆蓋影響因素進行分區(qū)探討,根據(jù)探測結果進行更有針對性的原因分析。探測結果表明,對整個區(qū)域來說,海拔的影響是最大的,原因在于研究區(qū)海拔南北差異極大,對整個區(qū)域的地形地貌起著宏觀控制作用,進而影響局部地區(qū)的水熱條件,對植被空間分異產(chǎn)生重要影響。此外,影響較大的是氣溫、降水、蒸散發(fā)等氣象因子,這與許多研究以及現(xiàn)實情況也是相符的[51-53],植被受氣候影響較大,特別是怒江-薩爾溫江流域跨越了多個不同氣候區(qū),氣候差異大,因此氣候因子對其影響較顯著;人為因子對全區(qū)植被的影響程度比自然因子較低,其中以土地利用影響最大,原因在于研究區(qū)地形偏遠,人口經(jīng)濟發(fā)展較緩慢,特別是高海拔地區(qū)和高山峽谷區(qū),人類活動很難影響到,而土地利用作為能較好反映人類活動的因素,能很好的表現(xiàn)植被覆蓋的變化,同時也有很大的可利用空間來增強或者減弱植被覆蓋。

對各子流域來說,植被覆蓋的影響因素具有差異性。上游區(qū)域土地利用對植被覆蓋的影響大于海拔,原因在于區(qū)域內部大部分屬于高原,高原內部地勢起伏不大,因此海拔以及坡度、坡向等對植被的影響小于土地利用,人類活動通過改變土地利用結構,特別是耕地面積、耕作類型等從而影響植被覆蓋;在剩下的因子中降水是影響較大的,原因在于上游地區(qū)季風難以到達,降水又是植被生長的重要因素,因此,降水對植被影響較大。由于其較為偏遠,人為因素影響小于自然因素。中游區(qū)域,由于其大部分屬于高山峽谷區(qū),海拔差異顯著,海拔成為了植被覆蓋影響最大的因素,而氣溫、降水、蒸散發(fā)等受海拔差異的影響,對植被覆蓋影響也較大。下游區(qū)域,由于以平原為主,地形較為平坦,人類活動頻繁,因此人口和GDP是對植被影響最大的因素,由于緯度低,氣溫高,蒸散發(fā)也成為了植被變化的重要影響因素,與上游和中游地區(qū)不一樣的是海拔對植被影響微弱。

綜上可知,由于怒江-薩爾溫江流域地理特征的復雜性,在不同環(huán)境條件下,各環(huán)境因子對植被覆蓋的影響也存在顯著的差異性。因此,建議綜合考慮不同環(huán)境條件下FVC的空間分布與各影響因子的空間關聯(lián)性,明確不同區(qū)域的環(huán)境制約因子,因地制宜制定生態(tài)修復措施,規(guī)劃流域的經(jīng)濟建設。對于上游地區(qū)來說,植被覆蓋較低,未來植被呈退化趨勢的占比較多為9.03%,應防止植被覆蓋率進一步下降加強生態(tài)恢復,如優(yōu)化土地利用結構、發(fā)展現(xiàn)代生態(tài)農業(yè)、更改耕作類型、進行輪作和休耕管理、退耕還林、建設自然保護區(qū)和實施生態(tài)移民、以及爭取國家政策的扶持等,注重社會經(jīng)濟和生態(tài)環(huán)境的共同發(fā)展;對于中游地區(qū),自然環(huán)境因子對植被覆蓋影響更大,其經(jīng)過三江源等重要生態(tài)源地區(qū),應該確保生態(tài)環(huán)境不被破壞;對于下游地區(qū),植被覆蓋度較高,但影響植被覆蓋的主要是社會經(jīng)濟因素,且未來植被呈退化趨勢的占比是最多的為17.39%,應切實保護好現(xiàn)有資源,維持生態(tài)系統(tǒng)的穩(wěn)定性,并發(fā)展高效、多元化、可持續(xù)的生態(tài)經(jīng)濟。植被動態(tài)是一個復雜的過程,還應加強對極端事件的預防和控制,以防止氣候變化造成植被的進一步退化。總之,應充分發(fā)揮人類活動對于生態(tài)環(huán)境的積極作用,找準生態(tài)保護與經(jīng)濟發(fā)展的平衡點,促進生態(tài)與經(jīng)濟的協(xié)調可持續(xù)發(fā)展。

4.4 局限和展望

研究參考楊帆等[30]的成果重新將怒江-薩爾溫江流域劃分為了上游、中游和下游區(qū)域進行植被趨勢分析與驅動力探討,我們主要考慮的是海拔和行政區(qū)劃,區(qū)域劃分較為主觀,這是未來研究需要改善的地方,比如利用具體的氣候區(qū)、土壤區(qū)、地貌區(qū)、植被區(qū)等進行分析。此外,研究在進行驅動力指標選取的時候,由于跨境各個國家的數(shù)據(jù)統(tǒng)計標準不一,精度也有差異,使得部分數(shù)據(jù)獲取較為困難,特別是人為因素的數(shù)量和質量,這可能會影響結果的精度,因此,未來應該考慮在分區(qū)探討的同時對不同區(qū)域建立不同的指標體系,以提高數(shù)據(jù)的精度和結果的準確性。

5 結論

研究基于怒江-薩爾溫江流域2000-2021年MODIS NDVI數(shù)據(jù),計算出其植被覆蓋度FVC,探討了研究區(qū)近22年植被覆蓋的時空趨勢變化、非線性特征、未來可持續(xù)性以及驅動力。主要結論為:

(1)近22年,怒江-薩爾溫江流域植被覆蓋總體發(fā)展趨勢向好,2007年和2015年出現(xiàn)較大突變,多年平均FVC值為0.73。受水熱條件、地形地貌等影響,植被分布具有明顯的空間異質性,下游和中游區(qū)域植被覆蓋明顯優(yōu)于上游區(qū)域,總體表現(xiàn)為南高北低,植被覆蓋狀況較好。(2)BFAST趨勢表明,近22年怒江-薩爾溫江流域植被覆蓋改善和退化的區(qū)域面積占比分別為71.24%、28.76%,改善的區(qū)域遠大于退化的區(qū)域,說明研究區(qū)植被得到較好的保護。Hurst指數(shù)顯示,未來植被相比以前將得到改善和將有所退化的區(qū)域占比分別為94.89%、2.76%。BFAST與Hurst二者的耦合表達了植被覆蓋的未來非線性趨勢,也是改善的區(qū)域大于退化的區(qū)域,占比分別為68.99%和29.09%。應特別關注未來發(fā)展趨勢退化的區(qū)域。(3)地理探測器結果表明各區(qū)域植被覆蓋影響因素具有差異性,針對全區(qū)域海拔的影響是最大的,進一步說明了地形與研究區(qū)植被具有密切關系。而針對各子流域,土地利用對上游地區(qū)植被覆蓋的影響比中游和下游區(qū)域顯著,中游地區(qū)海拔的影響力最大,下游地區(qū)以人口、GDP等人為因素影響為主。因此,要合理根據(jù)研究區(qū)不同區(qū)域的環(huán)境制約因子等制定合理的植被保護措施。