王志偉,岳廣陽,吳曉東

1 貴州省農(nóng)業(yè)科學(xué)院草業(yè)研究所,貴陽 550006 2 中國科學(xué)院西北生態(tài)環(huán)境資源研究院 冰凍圈科學(xué)國家重點(diǎn)實(shí)驗(yàn)室 青藏高原冰凍圈觀測研究站, 蘭州 730020 3 美國德克薩斯州大學(xué)圣安東尼奧分校地質(zhì)學(xué)系,美國圣安東尼奧 78249

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

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

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

1 研究區(qū)概況

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

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

2 數(shù)據(jù)來源與研究方法

2.1 數(shù)據(jù)來源及預(yù)處理

本研究中涉及的ASAR (Advanced Synthetic Aperture Radar,ASAR)數(shù)據(jù)有28景。因基線和多普勒質(zhì)心差的原因,其中4景影像沒有參與配對,剩余24景參與運(yùn)算的影像獲取年月日(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數(shù)據(jù)集屬于單視復(fù)數(shù)影像(Single Look Complex, SLC),為C波段、IS2模式遙感數(shù)據(jù),入射角23°。研究數(shù)據(jù)獲取方式為降軌,數(shù)據(jù)來源于歐空局(European Space Agency, ESA)。方位向分辨率和距離向分辨率分別為22.6 m和4.05 m。

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

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

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

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

2.2 地表形變反演方法

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

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

(1)

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

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

(2)

獲取具有物理意義的形變時(shí)間序列結(jié)果,其計(jì)算過程如下：

假設(shè)vk是i到j(luò)時(shí)刻的平均相位速度,則其相位和時(shí)間關(guān)系滿足如下公式：

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

(3)

式中,φ指相位信息。

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

(4)

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

Av=Δφ

(5)

式中,A是系數(shù)矩陣；速度v是速度矢量,可以通過系數(shù)矩陣A計(jì)算獲取。因?yàn)樵赟BAS-InSAR處理過程中會(huì)有多個(gè)主影像,這就有可能導(dǎo)致系數(shù)矩陣A秩虧。通常利用SVD方法來處理,首先會(huì)生成一個(gè)逆矩陣,然后獲得速度矢量的最小范數(shù)解,最終通過各個(gè)時(shí)間段內(nèi)的速度積分獲得各個(gè)時(shí)間段的形變量。

2.3 生態(tài)指標(biāo)計(jì)算方法

已有研究[28]指出,遙感生態(tài)指數(shù)包括綠度、濕度、熱度和干度4個(gè)指標(biāo),分別可以用植被指數(shù)NDVI、濕度分量Wet、地表溫度LST 和土壤指數(shù)NBSI來代表。

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

濕度指標(biāo)利用TM近似波段的MODIS 09A1產(chǎn)品計(jì)算獲得,其計(jì)算公式如下：

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波段的反射率,分別對應(yīng)TM數(shù)據(jù)的第1、2、3、4、5和7波段數(shù)據(jù),如表1所示。

熱度指標(biāo)在此直接利用MODIS 11A2 LST產(chǎn)品。

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

最后,為了后續(xù)的統(tǒng)計(jì)分析,除了干度指標(biāo)不做考慮外,其他3個(gè)指數(shù)都實(shí)行歸一化(Normalized Differential)操作,具體計(jì)算公式如下：

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

(7)

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

3 結(jié)果與分析

3.1 研究區(qū)地表形變時(shí)間序列反演結(jié)果

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

3.2 研究區(qū)地表形變率結(jié)果

研究區(qū)在2005年至2010年間的地表形變率如圖3所示,所有位置的形變率位于15 mm/a之內(nèi)。圖中無數(shù)據(jù)區(qū)域是形變率為0 mm/a和數(shù)據(jù)質(zhì)量較差的像元,其余區(qū)域地表形變率絕大部分都在±8 mm/a內(nèi)。如圖4所示,其中形變率在±4 mm/a以內(nèi)的像元占所有形變像元的89.24%,變形率為正的區(qū)域占57.70%,地表形變?yōu)樨?fù)的區(qū)域占42.30%。

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

研究區(qū)24個(gè)時(shí)間點(diǎn)在不同植被類型條件下的地表形變結(jié)果如圖5所示,高寒草地整體表現(xiàn)地表下沉的現(xiàn)象,同全球變暖凍土退化進(jìn)而導(dǎo)致冰融化為水造成土層體積減少的規(guī)律一致。

3.4 不同植被類型條件下生態(tài)指數(shù)結(jié)果

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

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

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

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

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

3.5 地表形變與生態(tài)指數(shù)的相關(guān)性

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

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

4 討論

4.1 SBAS-InSAR技術(shù)的優(yōu)勢及局限性

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

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

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

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

4.2 生態(tài)遙感指數(shù)對不同高寒草地類型地表形變的影響分析

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

5 結(jié)論

以上研究利用成熟的SBAS-InSAR技術(shù),對青藏高原多年凍土區(qū)的五道梁研究區(qū)進(jìn)行了地表形變特征分析,準(zhǔn)確的反演出研究區(qū)高分辨率的地表形變信息。同時(shí)借助生態(tài)遙感指標(biāo)對研究區(qū)高寒草甸和高寒草原兩種不同的植被類型進(jìn)行分析,發(fā)現(xiàn)高寒草原的地表沉降狀況要嚴(yán)重于高寒草甸。以上結(jié)果反映出在青藏高原多年凍土區(qū),植被生長越差的區(qū)域,其生態(tài)系統(tǒng)的脆弱性更為強(qiáng)烈。