青藏高原多年凍土區不同高寒草地類型地表形變特征

王志偉,岳廣陽,吳曉東

1 貴州省農業科學院草業研究所,貴陽 550006 2 中國科學院西北生態環境資源研究院 冰凍圈科學國家重點實驗室 青藏高原冰凍圈觀測研究站, 蘭州 730020 3 美國德克薩斯州大學圣安東尼奧分校地質學系,美國圣安東尼奧 78249

青藏高原作為“世界第三極”,不僅在亞洲季風系統中發揮著重要作用,對全球氣候變化也異常敏感[1- 2]。它不僅是我國氣候系統的重要組成因子,還影響到全球尺度的氣候變化[3- 4]。被稱為“氣候變化指示器”的多年凍土[5]廣泛分布于青藏高原[6],面積約為1.06×106km2[7]。多年凍土指溫度能夠維持在零攝氏度以下狀態兩年及兩年以上的近地表土壤或巖石層,活動層是多年凍土區位于多年凍土之上,夏季融化、冬季凍結的地表部分。伴隨全球氣候變暖[8],多年凍土在逐年退化[9],活動層厚度相應增加[10],影響到土壤特別是近地表的水、熱循環過程和生態環境[11]。多年凍土的熱融效應增長會促使大團聚體破碎成小團聚體,釋放大量有機碳、硝態氮等物質,進而造成地表植被發生改變,影響到地表的一系列特征,如反照率、降水的滲透速度、土壤中的蒸騰和蒸散、以及土壤侵蝕等,從而打亂水文和氣候系統的循環速率,造成高寒地區的生態環境惡化。

青藏高原多年凍土區蘊含世界上海拔最高、類型最為獨特的高寒草地生態系統[12],其分布面積分別占青藏高原多年凍土區和全國草地總面積的80%[13]和38%[14],具有防風固沙、涵養水源、固氮儲碳、調節碳循環及氣候變化、維護生物多樣性等諸多生態服務功能[15- 16],是國家生態安全的重要屏障[14]。已有研究表明,高緯度和高海拔地區的生態系統對氣候變暖的響應更加敏感[17]。而青藏高原的高寒草地作為高原氣候和凍土環境雙重屬性的生態系統,其對全球氣候變化的響應更迅速、更超前[18],其生態系統的穩定性也更脆弱[19]。綜上所述,對青藏高原多年凍土區高寒草地生態系統地表環境變化的特征分析極為重要,可以為揭示高寒草地生態系統在全球變化中的生態價值和貢獻提供重要的理論基礎和科學依據[12]。

隨著星載合成孔徑雷達干涉測量技術(Interferometric Synthetic Aperture Radar, InSAR)的發展[20],可以實現對多年凍土區近地表凍融循環變化特征的大范圍、高精度和高分辨率監測[21],為寒區生態環境的安全、平穩和可持續發展提供了十分有效的監測工具和技術手段[22]。近些年利用InSAR技術分析多年凍土區的相關研究日益成熟[23- 24],Liu等[25]在阿拉斯加多年凍土區利用ERS 1/2 SAR數據成功探測到厘米級的地表沉降量,之后[26]進一步反演該區域的活動層厚度。Chen等[21]使用ALOS PALSAR和Envisat ASAR數據,獲得青藏高原北麓河的地表形變量,結果顯示部分區域每年的變化量甚至為±2 cm。Zhao等利用Envisat ASAR數據[24]對青藏高原當雄至羊八井的凍土區地表形變進行研究,發現自然地表的形變量在3.6—5.0 cm,而鐵路和公路的形變量則在2.8—3.7 cm。上述研究多只是單純針對地表形變(Ground surface deformation, GSD)和活動層厚度進行探討[26],而發掘不同高寒草地生態環境條件下地表形變特征變化的研究還較為薄弱[27]。鑒于此,本研究以青藏高原多年凍土區五道梁區域的不同類型高寒草地為研究對象,通過分析其地表形變量和遙感生態指數(Remote Sensing Based Ecological Index, RSEI)的關系(本研究中的生態指數包括綠度、濕度和熱度三項[28]),試圖揭示：(1)不同高寒草地類型條件下的地表形變量變化特征如何？(2)地表形變量變化與綠度、濕度和熱度指標存在何種關系？旨在對不同高寒草地類型地表形變特征進行挖掘,以探討青藏高原多年凍土區不同高寒草地生態系統對氣候變化的響應機制。

1 研究區概況

研究區位于青藏高原腹地、青海省的西部,大片連續多年凍土區(圖1),地處北緯34.4°—35.9°,東經92.4°—93.9°之間,平均海拔(4600±190) m,屬高原山地氣候,夏季降水多,冬季降水少[29]。2005—2010年期間,年均植被指數為0.085±0.037(-),屬于植被稀疏區,主要以高寒草地植被類型為主,包括高寒草甸和高寒草原兩種。

圖1 青藏高原多年凍土分布及研究區位置Fig.1 Permafrost distribution of Qinghai-Tibet Plateau and location of study area

2 數據來源與研究方法

2.1 數據來源及預處理

本研究中涉及的ASAR (Advanced Synthetic Aperture Radar,ASAR)數據有28景。因基線和多普勒質心差的原因,其中4景影像沒有參與配對,剩余24景參與運算的影像獲取年月日(YYYY-MM-DD)分別為：2005-4-7、2005-9-29、2006-4-27、2006-10-19、2007-5-17、2007-10-4、2008-5-1、2008-10-23、2009-1-1、2009-2-5、2009-3-12、2009-4-16、2009-5-21、2009-6-25、2009-7-30、2009-9-3、2009-10-8、2009-12-17、2010-1-21、2010-2-25、2010-4-1、2010-5-6、2010-6-10和2010-7-15。該ENVISAT ASAR數據集屬于單視復數影像(Single Look Complex, SLC),為C波段、IS2模式遙感數據,入射角23°。研究數據獲取方式為降軌,數據來源于歐空局(European Space Agency, ESA)。方位向分辨率和距離向分辨率分別為22.6 m和4.05 m。

DEM資料為STRM V4版本數據,該數據在赤道的分辨率約為90 m[24]。

本文使用的MODIS產品有09A1、11A2和13A1三類,每類產品選取Terra(MOD產品)衛星結果。09A1產品為500 米8 天合成數據集,選取其1、2、3、4、6和7六個波段的數值(表1)計算生態指數。地表溫度(Land Surface Temperature,LST)11A2和增強型植被指數(Enhanced Vegetation Index,EVI)13A1產品分別為1 km 8天和500 m 16天合成數據集,運算過程中對LST數據執行降尺度處理至500 m。植被類型數據結果從全國植被類型圖中提取[30],重采樣至500 m。

此外,文中所有地圖投影方式都為WGS84坐標系。

表1 MODIS 09A1與TM各波段波長

2.2 地表形變反演方法

地表形變反演,采用短基線集合成孔徑雷達干涉測量技術(Small Baseline Subset Interferometric Synthetic Aperture Radar, SBAS-InSAR)。該技術使用的微波數據,不受或少受云、雪和稀疏植被的影響,不需要光照條件,可全天候、全天時地獲取。SBAS-InSAR技術能夠獲取長時間序列的地表形變特征,成熟應用于青藏高原多年凍土區[21, 24]。其原理是通過構建具有較短時-空基線的影像數據對,計算獲得干涉圖,然后應用奇異值分解(Singular Value Decomposition,SVD)方法生成形變時間序列和平均形變速率結果。SBAS-InSAR干涉時序分析方法計算獲取影像形變結果基本計算模型如下式所示：

φint=φdef+φflat+φtopo+φatmo+φnoise

(1)

式中,φint指干涉相位,φdef指地表形變相位,φflat指平地相位,φtopo指地形相位,φatmo指大氣延遲相位；φnoise指噪聲引起的相位。差分干涉測量主要關注的是地表形變相位,該結果并非借助高程數據(高程信息會作為輔助數據參與運算)反演的高程差,而是相對于第一景影像的形變差(因此第一景影像的形變值為0)。

具體來講,在研究區內參與運算的影像共有N+1景,獲取時間按照順序排列,依次為t0、t1、…、tN。以其中任意一幅影像為例,該影像會與其他影像進行匹配,生成至少一幅差分干涉圖。假設第k幅差分干涉圖由ti和tj時刻獲取的影像生成,則其計算過程如下式所示：

(2)

獲取具有物理意義的形變時間序列結果,其計算過程如下：

假設vk是i到j時刻的平均相位速度,則其相位和時間關系滿足如下公式：

vk×(tj-ti)=Δ(φj-φi)

(3)

式中,φ指相位信息。

考慮到不同的時間分段,T0到Tk時刻的第k幅干涉圖的相位信息可以通過如下公式計算得到：

(4)

式中,Δφk是影響不同時間間隔的相位速度積分信息,改寫成矩陣表達式為

Av=Δφ

(5)

式中,A是系數矩陣；速度v是速度矢量,可以通過系數矩陣A計算獲取。因為在SBAS-InSAR處理過程中會有多個主影像,這就有可能導致系數矩陣A秩虧。通常利用SVD方法來處理,首先會生成一個逆矩陣,然后獲得速度矢量的最小范數解,最終通過各個時間段內的速度積分獲得各個時間段的形變量。

2.3 生態指標計算方法

已有研究[28]指出,遙感生態指數包括綠度、濕度、熱度和干度4個指標,分別可以用植被指數NDVI、濕度分量Wet、地表溫度LST 和土壤指數NBSI來代表。

考慮到青藏高原多年凍土區植被稀疏的特點[31],本研究中使用MODIS 13A1的EVI產品替代NDVI作為綠度指標。

濕度指標利用TM近似波段的MODIS 09A1產品計算獲得,其計算公式如下：

Wet(MODIS)=0.0315b3+0.2021b4+0.3102b1+0.1594b2-0.6806b6- 0.6109b7

(6)

式中,b3、b4、b1、b2、b6和b7分別為MODIS 09A1的第3、4、1、2、6和7波段的反射率,分別對應TM數據的第1、2、3、4、5和7波段數據,如表1所示。

熱度指標在此直接利用MODIS 11A2 LST產品。

因研究對象為青藏高原多年凍土區的高寒草地,在此有植被區域不適合用裸土指數。同時,研究區內較少有建筑存在,故在本研究中也不考慮土壤指數。因此本研究中未對干度指標進行計算和分析。

最后,為了后續的統計分析,除了干度指標不做考慮外,其他3個指數都實行歸一化(Normalized Differential)操作,具體計算公式如下：

NDRSEI=(RSEIx- RSEImin)/(RSEImax-RSEImin)

(7)

式中,NDRSEI、RSEIx、RSEImin和RSEImax分別對應綠度、濕度和熱度3個指標的歸一化值、當前值、最小值和最大值,計算后的結果包括歸一化綠度(Normalized differential green index,NDGI)、歸一化濕度(Normalized differential wet index,NDWI)和歸一化熱度(Normalized differential land surface temperature index,NDLI)。

3 結果與分析

3.1 研究區地表形變時間序列反演結果

地表形變的長時間序列數據結果主要覆蓋2005—2010年時間段,如圖2所示。利用SBAS-InSAR方法反演研究區24個時間段(共采集28景影像,4景沒有參與配對的影像無結果)的形變結果。其中,于2005年4月7日獲取的第一景運算影像被設置為參考影像(即t0時刻影像),該影像所有數值皆為0,故未在圖2中展示。此外,圖2中紅色區域代表地表存在抬升現象；綠色則相反,代表地面的沉降現象。

3.2 研究區地表形變率結果

研究區在2005年至2010年間的地表形變率如圖3所示,所有位置的形變率位于15 mm/a之內。圖中無數據區域是形變率為0 mm/a和數據質量較差的像元,其余區域地表形變率絕大部分都在±8 mm/a內。如圖4所示,其中形變率在±4 mm/a以內的像元占所有形變像元的89.24%,變形率為正的區域占57.70%,地表形變為負的區域占42.30%。

3.3 不同植被類型條件下地表形變結果

研究區24個時間點在不同植被類型條件下的地表形變結果如圖5所示,高寒草地整體表現地表下沉的現象,同全球變暖凍土退化進而導致冰融化為水造成土層體積減少的規律一致。

3.4 不同植被類型條件下生態指數結果

利用公式7中的方法,計算歸一化生態指數(Remote sensing based ecological indexes,NDRSEI)獲取結果如圖6所示。圖中歸一化熱度(Normalized differential land surface temperature index,NDLI)體現高寒草甸類的NDLI值大于高寒草原類,但是利用LST計算獲取的NDLI結果中并沒有當即體現出2008年極端氣溫低年,而是在2009年3月后開始出現一次極端低地溫現象。歸一化綠度(Normalized differential green index,NDGI)顯示高寒草甸類的NDGI值要大于高寒草原類,特別是在2008年下半年時,存在NDGI顯著變小的現象。歸一化濕度(Normalized differential wet index,NDWI)不僅較好的體現出濕度季節性變化的特點,而且曲線在2005年4月至2010年7月期間呈現一種濕度穩定增長的狀態。

圖3 研究區2005年到2010年地表形變率圖 Fig.3 Deformation rate of ground surface in the study area from 2005 to 2010 紅色為抬升區域,綠色為沉降區域

圖4 研究區地表形變率頻率分布直方圖 Fig.4 Frequency distribution histogram of ground deformation rate in study areas

圖5 2005年到2010年間研究區不同高寒草地類型地表形變Fig.5 Ground deformation of different alpine grassland types from 2005 to 2010

圖6 2005—2010年間不同高寒草地類型生態指數Fig.6 Ecological indexes of different alpine grassland types from 2005 to 2010

3.5 地表形變與生態指數的相關性

為分析研究區2005年4月至2010年7月3種生態指數對地表形變的貢獻與作用,本文利用研究區地表形變(Ground surface deformation,GSD)和3種歸一化生態指數,包括NDLI(熱度)、NDGI(綠度)和NDWI(濕度)結果,計算獲得圖7所示的相關圖。

圖7 高寒草甸和高寒草原條件下地表形變與3種生態遙感參數的相關性Fig.7 Correlations between ground surface deformation and different remote sensing based ecological index in alpine meadow and steppe

4 討論

4.1 SBAS-InSAR技術的優勢及局限性

InSAR技術提供了一種大范圍地表形變估算方法,特別是SBAS-InSAR技術尤為適用于青藏高原多年凍土區,同時多種SAR影像傳感器的在軌運行,也為長時間序列的地表形變研究提供了必要的數據基礎[26]。同時借助現有的植被類型圖,甚至可以在青藏高原多年凍土區腹地開展不同植被類型條件下的地表形變特征研究。而高原凍土區腹地一般而言多位于海拔高、人類生存環境惡劣的區域。因此,本文嘗試利用現有的SAR影像和植被類型分類結果,利用SBAS-InSAR技術在五道梁地區(青藏高原多年凍土區腹地)開展對不同高寒草地地表形變特征的分析,并得到一些初步結果。

雖然本研究成功的利用InSAR技術、植被類型圖集和生態遙感植被對青藏高原多年凍土區進行了地表形變技術特征分析,并得到一些初步的結論,但仍存在以下幾方面需要繼續加強和改進。

(1)SAR數據過于零散。特別是在圖6NDGI曲線中,2007年下半年的點甚至未表現出植被季節性變化的特點。在以后的研究中,需要增加每一年的數據獲取量,以此來更準確的反映地表形變的客觀規律。

(2)凍土區的研究數據存在時間滯后現象。如圖6NDLI曲線中,LST產品對2008年極端低氣溫年的滯后表現。時間滯后性的研究,在多年凍土區的重要性日益凸顯,在下一步的研究中,將會利用合理、有效的模型、算法來挖掘這種規律。

4.2 生態遙感指數對不同高寒草地類型地表形變的影響分析

如圖2所示,呈現大規模抬升的區域主要位于研究區東北部,該區域包括3個湖泊,從西至東依次為庫塞湖、海丁諾爾湖和鹽湖。研究區東北部的大范圍抬升可能受區域內湖泊影響,而湖泊區域外,大部分區域表現地面沉降的特點,該結果同當雄至羊八井段的地表形變特征一致[10]。此外,通過圖5可知,在2008年附近,地表形變存在短期抬升的現象,可能同2008年我國的極端低溫事件發生有關。同時,圖5中也展示出高寒草原的地表沉降現象明顯強于高寒草甸地區,符合植被條件較好的區域,凍土更加穩定,地下冰融化慢的規律。在圖7中,對于高寒草甸而言,地表形變量和3種生態指數的負相關性略大于高寒草原。其中在高寒草甸條件下地表形變與NDLI的負相關性最大,數值為-0.0285,與NDWI和NDGI的負相關性數值依次為-0.0241和-0.0214。而在高寒草原條件下,地表形變量則與NDWI的負相關性最大為-0.0121,然后為NDLI和NDGI,負相關性數值分別為-0.0111和-0.0107。以上研究結果表明,針對不同草地類型,其主導因子存在差異。高寒草甸的地表形變有可能更多的受限于溫度變化(NDLI熱度指標,由地表溫度LST產品計算獲得),而高寒草原的地表形變則可能更多的由水分條件所影響(NDWI濕度指標,在遙感生態指數中為表征濕度的分量)。

5 結論

以上研究利用成熟的SBAS-InSAR技術,對青藏高原多年凍土區的五道梁研究區進行了地表形變特征分析,準確的反演出研究區高分辨率的地表形變信息。同時借助生態遙感指標對研究區高寒草甸和高寒草原兩種不同的植被類型進行分析,發現高寒草原的地表沉降狀況要嚴重于高寒草甸。以上結果反映出在青藏高原多年凍土區,植被生長越差的區域,其生態系統的脆弱性更為強烈。