格爾木長時間序列遙感生態指數變化監測

汪東川,陳 星,孫志超,辛 燕,王海慶,柴 華,王鴻藝

1 天津城建大學地質與測繪學院,天津 300384 2 天津城建大學/天津市土木建筑結構防護與加固重點實驗室,天津 300384 3 天津綠茵景觀生態建設股份有限公司,天津 300384

資源型城市依托資源產業的開發而興建,由于資源的不可再生性,過度開發在帶來社會貢獻與經濟效益的同時對生態環境造成破壞[1],監測生態環境質量分布與變化對生態環境保護有重要作用。生態系統為人類生產生活提供了重要的物質基礎[2],也影響著居民的身心健康、社會的文明進步、人才的吸引和經濟發展的長期穩定性[3]。對生態環境進行準確、快速的評價與分析,能為城市環境治理與綠色發展提供可靠的科學依據[4]。

遙感技術的發展使長時序對地觀測研究得以實現[5—7],利用多時相遙感數據獲取土地利用類型或地表覆蓋因子等信息,對地表動態變化特征進行量化表達,快速、準確、多尺度的分析變化規律。在長時間序列遙感的實際應用中,學者常選擇3—10年相同間隔年份進行分析[8],時間選擇缺少指向性,容易遺漏發生重要變化的年份。有研究通過觀察多年數據規律進行人為分段[9]或基于梯度下降法自動分段擬合劃分斷點位置[10],利用主觀判斷或數學統計分析方法選取特征時間節點,對長時間序列進行分段。本研究引用Mann-Kendall突變檢驗法對格爾木RESI均值進行非線性檢驗,根據突變點對長時間序列數據分段。

目前利用遙感數據研究生態環境多通過地表覆蓋類型[8]或植被覆蓋指數[11]、不透水面指標[14]、地表溫度[12]、地表濕度[13]等生態指標來表示或量化生態環境的優劣情況[15]。然而生態環境是一個綜合動態系統,單一因子無法全面反映其變化,因此如何構建綜合評價指標是監測生態環境質量研究的重要問題。近年來,有研究者計算多種指數均值[16]或以面積權重法構建綜合模型[17],也有利用變異系數法[18]、層次分析法[19]等地學統計方法構建綜合模型,但結果往往受人為主觀因素影響。2013年徐涵秋[20]綜合多方面生態因子提出新型遙感生態指數(Remote Sensing Ecological Index,RSEI),以主成分分析法綜合綠度、濕度、熱度和干度4個指標,根據各指標對主分量的貢獻值確定權重,綜合、客觀的定量評價生態環境狀況。近年來,利用RSEI對不同空間尺度的生態環境進行監測與評價的研究證明RSEI具有可擴展性、可視化和可比性[21—23]。由于RSEI指標選取全面,在長時間序列研究中數據量多,處理工作重復繁瑣。谷歌地球引擎(Google Earth Engine,GEE)強大的數據運算功能可以批量、快速的處理影像數據集[24],因此使用GEE云平臺預處理并計算各項指標,不僅提高了效率,還可以進行最小云量篩選以獲取質量最優的影像進行研究[10,25]。

格爾木擁有豐富的礦產資源,以資源開采為主的人類活動造成生態環境嚴重破壞。近年來,學者們對格爾木地區的生態環境進行了各方面的研究[26—27],但對格爾木全區綜合生態指數的研究較為少見。本文以格爾木為研究區,通過GEE計算1990—2019年遙感生態指數,進一步分析格爾木生態環境時空變化并探討其影響因素。

1 研究區概況、數據來源

1.1 研究區概況

格爾木地處青藏高原,位于青海省西部,是青藏高原的第三大城市,由柴達木盆地中南部和唐古拉山鎮兩塊互不相連的區域組成,唐古拉山區域內有一級河流長江源頭沱沱河(圖1)。隨著全球氣候變化與人類活動的影響,區內土地荒漠化日益嚴重,生態系統脆弱敏感。受地質構造作用,格爾木擁有眾多礦產資源,對礦產資源的粗獷式開采加劇了資源的破壞和環境的惡化[28]。格爾木是新疆、甘肅、四川、西藏、青海五省(自治區)的中心樞紐地帶,是青海省兩個國家級交通樞紐站之一[29],該地區生態環境可持續發展關系到青藏地區的經濟繁榮和社會穩定。

圖1 研究區位置

1.2 數據來源及預處理

研究區矢量數據來自中國科學院資源環境科學與數據中心(https://www.resdc.cn/)。遙感數據利用GEE平臺編程(JavaScript API)在線調用經過輻射校正、大氣校正的Landsat 5/8地表反射數據集,對6—9月影像進行云量篩選,鑲嵌合成格爾木近30年(1990—2019年)最小云量影像。通過歸一化水體指數(Normalized Difference Water Index,NDWI)閾值法對水體掩膜處理,以減少水域對濕度指標的影響[30]。在GEE平臺進行波段運算,得到各年份綠度、干度、熱度與濕度四種指數,通過主成分分析得到遙感生態指數(RSEI),再利用ArcGIS與MATLAB等軟件進行統計分析,本文技術路線(圖2)。

圖2 技術路線

2 研究方法

2.1 指數計算

根據徐涵秋提出的遙感生態指數(RSEI)[20],選擇歸一化差值植被指數(Normalized Difference Vegetation Index,NDVI)表示綠度指標,反映區域內植被生長情況及其分布密度[11]；建筑指數(Index-based Built-up Index,IBI)與土壤指數(Soil Index,SI)求均值得到干度指標(Normalized Difference Soil Index,NDSI),反映土地荒漠化與城鎮化光譜響應[31]；單窗算法反演得到地表溫度(Land Surface Temperature,LST)表示熱度指標[32],反映不透水面熱平衡[33]；纓帽變化計算濕度指標(Wet),通過土壤含水量反映土地退化情況[34—35]。各因子計算公式見表1：

表1 指標計算方法

由于各指標在數值單位和大小上有差異,因此要進行標準化處理[20],再通過主成分分析計算初始遙感生態指數(RSEI0)。標準化公式如下：

(1)

式中，N表示標準化后的指標值,I表示該指標在某像元的值,Imin為該指標的最小值,Imax為該指標的最大值。

基于主成分分析法計算初始遙感生態指數(RSEI0)：

RSEI0=1-PC1[f(NDVI,NDSI,LST,Wet)]

(2)

式中，PC1為第一主成分；f是對各指標進行標準化處理。同樣為便于指標的度量和比較也要對RSEI0進行歸一化處理,得到遙感生態指數(RSEI),值越大表示生態質量越好,反之生態質量越差。

2.2 生態指數動態變化分析

2.2.1Mann-Kendall突變檢驗

Mann-Kendall突變檢驗廣泛應用于水文序列數據與氣象站點數據,是常用的非參數統計方法,不受異常值干擾[36]。區域遙感生態指數均值常用來判斷區域生態總體變化趨勢,因此本研究對格爾木30年遙感生態指數均值序列進行非參數檢驗,以獲取生態環境突變點對長時間序列數據進行分段研究。

在Mann-Kendall檢驗中,對長度為n的時間序列R(x1,x2,…,xn)構造秩序列Sk：

(3)

假定時間序列隨機獨立的情況下,定義標準化統計量UFk：

(4)

式中,E(Sk)與var(Sk)為Sk的期望和方差。

再構建逆時間序列,計算統計量定義為UBk,滿足UBk=-UFk。當UFk>0時,序列為上升趨勢；UFk<0,序列為下降趨勢,UFk、UBk在置信水平上有交點,則該時刻為發生突變時間[37]。

2.2.2變化軌跡法

變化軌跡法利用代碼完整的表征時間維度的變化過程[37],有利于深入探討某一屬性動態變化規律,在土地利用變化的研究中有十分廣泛的應用[38]。本研究基于遙感生態指數等級,引用變化軌跡法對格爾木RSEI變化軌跡進行可視化。對變化軌跡編碼計算公式進行簡化：

(5)

式中,Tt為變化軌跡代碼,研究節點共有n個,xi為第i個節點對應的RSEI等級代碼(1為差、2為較差、3為中等、4為良、5為優)；從左至右,變化軌跡編碼首位數字為初始年份RSEI等級,第二位數字為第二節點年份RSEI等級,之后以此類推。

2.2.3Sen+Mann-Kendall趨勢分析

Sen+Mann-Kendall顯著性趨勢檢驗是長時間序列趨勢判斷的重要方法[39],由于Sen氏趨勢分析本身不能判斷序列趨勢的顯著性,可以通過Mann-Kendall檢驗Sen氏趨勢分析結果實現顯著性分析。Sen氏趨勢分析通過計算時間序列的中值,能較好的消除異常值的干擾,計算公式如下[40]：

(6)

式中,QRSEI表示時間序列數據趨勢的變化；mean為取中值函數；10.0005時趨勢上升,當|QRSEI|≤0.0005時趨勢穩定不變,當QRSEI<-0.0005時趨勢下降。

Mann-Kendall檢驗統計量S的計算過程如下：

(7)

(8)

當n≥10時,S近似服從標準正態分布,使用檢驗統計量Z進行趨勢檢驗：

(9)

其中var表示方差。采用顯著性水平α=0.05進行顯著性檢驗,即|Z|≥1.96時發生顯著變化,|Z|<1.96時發生不顯著變化。

3 結果與分析

3.1 生態環境整體分析

根據格爾木RSEI主成分分析結果(表2)可知：(1)第一主成分(PC1)平均貢獻率達到73%以上,說明PC1集中了4個指標的大部分特征信息,能夠表示區域內生態環境情況。(2)綠度與濕度的載荷值為正,干度與熱度的載荷值為負。符合綠度和濕度越高,地表植被覆蓋越高,土壤水分越充足,生態環境質量越好；干度和熱度越高,土壤沙化、巖石裸露等問題越嚴重的自然特征[41]。(3)對比PC1載荷值的變化,LST的絕對值遠大于其他指標,且NDSI與LST載荷值的絕對值之和始終大于NDVI與Wet之和,說明格爾木土地荒漠化與不透水面熱平衡對生態環境的破壞作用大于植被和土壤濕度的優化作用。綜上所述,RSEI指數可以綜合反應格爾木生境質量。

表2 主成分分析結果(第一主成分的載荷和貢獻率)

圖3為1990—2019年各指標及RSEI均值變化情況,RSEI在0.31與0.43之間波動變化,RSEI由0.38下降至0.33,整體降低13%。從單一指標來看綠度變化保持平穩,濕度明顯下降30%,干度下降17%,而熱度上升12%。綜合可知各指標對格爾木RSEI的影響呈下降趨勢,這與格爾木市RSEI變化趨勢相一致。

圖3 生態指標及遙感生態指數均值變化

3.2 格爾木市遙感生態指數突變分析

對1990—2019年的RSEI均值進行Mann-Kendall突變檢驗,結果如圖4。UFk曲線從2001年之后持續低于0,表明格爾木生態質量從2001年開始有變差趨勢；UFk曲線與UBk曲線在2006年存在有效突變點,且統計量小于0；UFk曲線在2015年之后超出顯著性水平。這與年均RSEI變化趨勢基本相同,因此根據突變分析的結果,本研究決定以1990、2001、2006、2015、2019年為時間節點,對格爾木長時間序列數據進行分段。

圖4 RSEI突變檢驗(1990—2019)

對RSEI進行等級劃分,并統計各生態等級占比(圖5)。由圖可知：等級為優的區域所占比例較為穩定,且主要集中在區內冰川積雪地帶；等級為良好的部分占比較少,從1990年的2.60%降至2019年的1.50%,主要位于唐古拉山區域,從均勻分布到只集中在唐古拉山地區西北部,中部的分布明顯減少；中等級的區域主要分布在柴達木盆地南部山區和唐古拉山區域,在2001年大幅度減少后占比基本穩定；格爾木RSEI等級以較差為主,分別達到48.96%、46.91%、49.07%、59.04%和45.72%；而生態環境等級差的區域主要分布在柴達木盆地北部,1990年僅占12.27%,之后基本穩定在25%上下浮動。

圖5 RSEI等級分布(1990—2019)

通過變化軌跡法得到軌跡代碼,表示RSEI等級變化過程。RSEI等級未發生變化的區域占格爾木總面積的61.94%,其中穩定不變的軌跡占總面積的35.93%。根據軌跡代碼面積,保留發生變化的前15種軌跡代碼,占格爾木總面積的30.26%(圖6)。其中RSEI等級向好發展的區域僅占總面積的3.55%,主要為較差轉為中等,其分布范圍在柴達木盆地區域南部和唐古拉山地區東北部。RSEI等級降低的區域占總面積的27.25%,其中2001年發生變化后保持穩定的區域占14.41%,結合等級變化百分比(表3),說明格爾木RSEI等級在2001年下降最為明顯；由較差轉向差的區域占13.50%,主要分布在柴達木地區北部；由中等轉為較差的區域占13.19%,主要分布在唐古拉山地區。根據相關研究,格爾木平均年氣溫在研究時段內有增長趨勢[42],導致冰川融化速度加快,并且降水量上升。短期內河流徑流量增加,周邊植被覆蓋增加,長期發展導致凍土地區植被退化,荒漠化加重,生境質量下降[43]。工業發展對環境的污染也不可忽視,但結合表3結論,RSEI在2015年之后下降幅度有所減緩,說明格爾木城市發展在改善區域生態環境工作中做出了積極的貢獻。

表3 RSEI等級變化百分比/%

圖6 RSEI等級變化軌跡

3.3 格爾木市遙感生態指數趨勢分析

通過Sen+Mann-Kendall趨勢分析良好的反映了1990—2019年格爾木RSEI變化趨勢的空間分布特征。表4統計可知,格爾木生態環境質量以下降為主,占區域總面積的89.73%,上升區域僅占3.24%,穩定不變即沒有發生顯著變化的區域為7.03%。

表4 1990—2019年RSEI變化狀況統計

根據圖7變化趨勢分布特點選擇典型區域,對其變化原因展開討論。格爾木RSEI顯著上升區域僅有0.29%。近年來國家加快柴達木循環經濟試驗區建設,對格爾木老工業基地進行高新技術改造,促進產業綠色發展,這一措施減緩了工業發展對地表溫度的影響,區域熱度指標有所下降。另一方面,三北防護林工程、退耕還林、防沙治沙等項目的加速發展對土地荒漠化治理已有成效,干度指標下降。在城市內部推進林草間作治理、建設農田網格防護林等措施在顯著提高區域植被覆蓋度的同時也降低了土壤蒸散量[44],使綠度與濕度指標同時上升。因此格爾木市區及烏圖美仁鄉RSEI主要受人類活動的積極影響顯著上升。RSEI不顯著上升區域占2.95%。瑪曲鄉位于可可西里自然保護區,根據李睿等[45]的研究,這一地區氣溫升高、降水量增加使氣候向暖濕化轉變,有利于植被生長,區域綠度指標上升明顯。而人類活動,如牧業生產和青藏公路、鐵路修建等,對該地區土地景觀也有一定影響,但整體影響較弱,干度指標微弱下降。在氣候變化與人類活動的共同影響下該區域RSEI呈不顯著上升。大灶火溝位于東昆侖造山成礦帶,主要以黑色及有色金屬冶煉產業發展為主,2005年國務院批復建立柴達木地區循環經濟試驗區[46],促進資源高效、循環利用實現區域可持續發展,有利于熱度指標下降,對RSEI有積極影響,使得該區域RSEI不顯著上升。RSEI不顯著下降區域占43.05%。唐古拉山鎮處于長江源頭沱沱河流域,氣候暖濕化使得植被覆蓋逐步改善,減緩草地退化趨勢,但綠度指標依然在下降。加之該地區土壤保持功能較差,降雨量增加也加劇了這一區域的土壤侵蝕[47],導致干度指標上升,因此氣候變化是唐古拉山鎮RSEI不顯著下降的主要原因。RSEI顯著下降區域占46.68%。近年來,察爾汗鹽湖工業區發展迅速,并提出了一系列環保措施,但鹽田開采面積增大影響植被生長[48],不可避免地導致綠度下降。柴達木盆地區域近30年降雨量變化不明顯,但降水日數略有增加,年平均溫度有升高趨勢,導致蒸發量增加,土地沙漠化面積增加,整體干度上升[49]。在氣候與人類活動共同的負面影響下造成了察爾汗鹽湖工業區RSEI顯著下降。因氣候變暖,格拉丹東冰川地區導致冰川消融增加[50],易引發山洪等地質災害,植被受到嚴重破壞,綠度下降。氣候變暖導致次生災害,使該地區RSEI顯著下降。

圖7 RSEI變化趨勢 (1990—2019年)

4 結論

本研究基于長時間序列遙感影像,構建格爾木市1990—2019年遙感生態指數(RSEI)。通過檢測突變點對時間序列數據進行分段研究,結合趨勢分析對格爾木生態環境進行時空可視化分析,主要結論包括：

(1)1990—2019年間,各數據的第一主成份(PC1)貢獻率平均值達到73%以上,綠度與濕度指標對生態環境有積極影響,干度與熱度指標有消極影響,與自然環境特征一致,且熱度指標對RSEI的影響最大,載荷值關系與RSEI均值變化趨勢相一致,說明RSEI可以綜合評價格爾木生態環境質量。30年間格爾木RSEI均值從0.38下降至0.33,表明格爾木生態環境整體有退化趨勢。

(2)根據突變檢驗結果,以1990、2001、2006、2015、2019年為節點對格爾木長時間序列數據分段研究。將RSEI劃分為5個等級,發現格爾木遙感生態指數等級分布以較差為主,占總面積的45%以上,主要分布在柴達木盆地北部地區。1990—2019年格爾木RSEI等級未發生變化的區域占61.94%,發生變化的軌跡中以等級下降為主,其中2001年下降最為明顯,而2015年之后下降幅度有所減緩。

(3)受不同程度的氣候與人類活動影響,格爾木遙感生態指數變化以下降為主。顯著上升區域集中分布在人類活動密集的格爾木市區及鄉鎮所在地；不顯著上升區域分布于唐古拉山中部瑪曲鄉及柴達木盆地南部成礦帶；不顯著下降區域大面積分布在唐古拉山鎮以及柴達木盆地南部地區；顯著下降區域主要分布在柴達木北部鹽湖工業區與唐古拉山南部冰川地帶。