1992—2009 年格拉丹東冰川及冰前湖面積變化的遙感研究

周文明，李志偉，李佳，崔志勇，汪長城

(中南大學 地球科學與信息物理學院，湖南 長沙，410083)

世界范圍內冰川的加速退縮是全球氣候不斷變暖的最明顯信號[1]。目前，全球氣候觀測系統(global climate observing system，GCOS)將冰川作為研究氣候變化的基本監測對象之一[2]。同時，冰川也是地球上重要的固體淡水資源，與人類的生產生活息息相關。處于世界“第三極”的青藏高原蘊藏了豐富的冰川資源。在全球變暖的大背景下，青藏高原大部分冰川呈明顯的退縮趨勢[3]。青藏高原冰川的變化監測對于研究全球變化和評估區域內徑流量的變化都具有重要意義。位于長江發源地的格拉丹東雪山發育了大規模冰川，這些冰川能深刻影響長江源頭生態環境和中下游社會發展，且規模較大，是研究氣候變化的理想標本。但是，由于該地區環境惡劣，地勢險峻，常規實地測量不能滿足大范圍、高空間分辨率、經濟和安全的冰川變化監測要求，而遙感技術憑借其面上的、無接觸式的觀測特點可以滿足這些要求[4]。魯安新等[5]利用人工目視判讀遙感分類方法研究了格拉丹東地區冰川在小冰期、1969 年和2000 年這3 個時期的面積變化，但沒有分析引起冰川變化的原因；Ye 等[6]利用非監督遙感分類和人工修整方法得到該地區冰川在1969—2002 期間的分布變化，并分析了冰川變化與沱沱河站和安多站的氣象數據之間的聯系，但其采用的2002年ETM+影像有季節性積雪覆蓋，雖然冰川邊界經過了人工處理，但結果仍然受到一定的影響，并且沒有研究該地區冰前湖的時空變化規律及與冰川變化之間的關系。上述研究只涉及了該地區冰川2002 年之前的變化，而未研究2002 年之后的變化情況。完善和驗證目前的冰川變化預測模型以及探究區域的氣候變化機制迫切需要最新的冰川變化參數[7]。針對目前格拉丹東冰川的研究狀況，本文作者綜合利用遙感和地理信息系統平臺，采用分類精度更高的方法得到格拉丹東冰川近20 a 來的變化規律，并分析引起這些變化的原因，同時揭示冰前湖的演化規律，發現冰前湖的演化與冰川變化密切相關，而且由于采用的分類方法不同，本文揭示的冰川退縮速度略大于文獻[6]中的結果。

1 研究區與數據源

1.1 研究區概況

格拉丹東雪山位于青藏高原腹地的唐古拉山西段山區(90.0°—91.5°N，33.0°—30.7°E)，海拔5 200～6 621 m，是長江源區的最高峰[5]。該地區冰川規模較大而且集中，冰川的分布形式多為星狀分布，冰川面積大約為1 000 km2(圖1 中虛線矩形框標注)，冰川類型為冷性大陸性冰川。格拉丹東雪山處于青藏高原腹地，受季風氣候影響很小，具有氣溫低和降水少的特點。該地區氣象站資料顯示年平均氣溫約為-4 ℃，每年5～9月份在0 ℃以上；年降水量在500 mm 以內，90%集中在5～9 月份。

1.2 數據源

圖1 研究區域Fig.1 Location of study area

Landsat 衛星的TM/ETM+影像覆蓋整個地球的時間序列較長，并且影像數據的空間分辨率和光譜分辨率均適合冰川監測，所以，TM/ETM+影像成為監測冰川變化的主要數據源，并被廣泛應用于世界冰川編目中[1-2,7-10]。通常作為冰川研究的遙感影像要滿足2 方面要求：1) 無云覆蓋或影像中的冰川區無云覆蓋；2) 影像獲取時間處于或接近冰川的消融末期[1,8-9]。

綜合考慮云、陰影及季節性積雪的影響，本文選取了1992 年、1999 年、2006 年和2009 年中多幅TM/ETM+數據，見表1。影像來源于美國地質調查局(United States Geological Survey, USGS)及全球土地覆蓋數據數據庫(global land cover facility,GLCF)。所有影像均為正射校正影像[7]。

表1 研究區域的TM/ETM+影像Table 1 TM/ETM+images of the study area

2 研究方法

目前基于光學遙感影像的冰川邊界提取方法有人工目視判讀、監督與非監督分類、波段比值法、NDSI(雪蓋指數法)[11-16]。對于表面無冰磧物覆蓋的冰川，波段比值法被認為是精確、快速、穩健的冰川提取方法，并廣泛應用于冰川編目中[2,7-10]。波段比值法的基本原理是利用冰和雪在可見光近紅外波段(波長為0.4～1.2 μm)的高反射特性和在短波紅外(波長為1.4～2.5 μm)的強吸收特性[10]。對于Landsat TM/ETM+影像，波段比值法主要應用f(TM3)/f(TM5)和f(TM4)/(TM5)這2 個參數(其中f(TM3)，f(TM4)和f(TM5)分別代表波段3，4 和5 的原始灰度)。在陰影區和薄巖屑覆蓋的冰川區域，f(TM3)/f(TM5)比f(TM4)/f(TM5)識別冰川能力更強，所以，比值f(TM3)/f(TM5)在冰川目錄的編譯中應用更廣泛[1-2,7-8]。但是，采用波段比值法對于冰川識別來說有2 個缺陷：冰前湖在遙感圖像上的光譜反射特性與冰雪的反射特性相似，使冰前湖容易被誤分為冰川；冰面上有巖屑覆蓋的部分與周圍裸露的巖石的光譜特性相似，使得這部分冰川不能識別[8]。

本文首先通過波段比值法得到二值圖，對二值圖進行矢量化，得到每個時期冰川的初始邊界矢量圖。有的影像無季節性積雪覆蓋，但部分區域受到云的影響。本文通過當年中的多期影像相互輔助，人工修整被云覆蓋的區域；通過將冰川矢量邊界疊加在假彩色合成的遙感影像上后人工修整冰川的邊界，消除冰前湖對冰川解譯結果影響[8,10]。另外，結合研究區影像的假彩色合成圖，采用人工目視判讀的方法提取該地區冰川湖邊界。在遙感分類精度評價中，一般采用高分辨率或相同分辨率人工解譯的結果作為參考數據。由于該地區缺少成像時間相同的高分辨率影像，因此，本文沒有定量評價冰川和冰川湖信息提取的精度。但是，本文采用波段比值和人工目視判讀相結合的方法提取冰川信息，采用人工目視判讀提取冰川湖信息，不僅兼顧了時間效率，同時也保證了冰川信息提取的精度。圖像處理的基本流程如下：

1) 將4 個時期的影像配準，本次實驗選取1992年TM 影像為主影像，其他影像都配準到這幅影像上，配準誤差都在1 個像元以內(低于30 m)。

圖2 1992 年研究區假彩色圖及冰川分類結果Fig.2 TM false-color composite and glacier mask of study area in 1992

2) 選取TM/ETM+影像的部分波段合成假彩色圖(R,G 和B 分布對應第5，3 和1 波段)，如圖2(a)所示。根據前面提到的波段比值f(TM3)/f(TM5)得到影像二值圖，二值圖中冰川像元值為1，非冰川區像元值為0。

3) 對冰川二值圖進行窗口為3×3 的中值濾波處理，這樣可以減少陰影中的噪聲，同時可以消除一些季節性雪塊[2,9,13]。濾波完成后將二值圖(圖2(b))進行矢量化處理，得到冰川的矢量邊界。

4) 將冰川的矢量邊界疊加在同時期的假彩色合成遙感圖像上，通過人工判讀檢查冰川邊界。修整陰影區的冰川邊界，同時刪除被誤分為冰川的水體(尤其是冰前湖)(見圖2(a)中的c 和f)；有云覆蓋的區域可以通過當年不同時期的遙感影像綜合判別來修改其邊界；在當年不同時期的遙感影像中選擇最小的冰川面積作為當年的統計冰川面積，可以最大限度地消除季節性雪的影響。

5) 統計各個時相冰川的面積參數，并計算該地區冰川1992—2009 年的分布變化數據。

3 結果與分析

3.1 冰川變化及分析

依據上述步驟獲得格拉丹東雪山冰川的變化情況，見表2。從表2 可以看出：1992—2009 年，格拉丹東冰川的總面積減少了66.68 km2(減小率為7.37%)。其中，在1992—1999 年，該地區冰川面積共減少19.90 km2(減小率為2.20%)，平均每年減少2.84 km2；在1999—2006 年，冰川總面積退縮了39.95 km2(減小率為4.52%)，平均每年減小5.71 km2，與1992—1999 年相比，每年減小的面積增加2.87 km2，冰川的退縮速度加快；在2006—2009 年，總面積減少6.83 km2(減小率為0.81%)，每年減小量為2.28 km2，冰川的退縮速度有所減緩。通過上述數據可知，自1992年以來，格拉丹東冰川總體處于連續退縮之中。

表2 1992—2009 年格拉丹東地區冰川的面積變化Table 2 Glacier area change in Geladandong from 1992 to 2009

文獻[5]和[6]獲取了格拉丹東雪山1969—2002 年的冰川面積變化信息。本文增加該地區2006 年和2009年的冰川面積參數，同時部分覆蓋了文獻[5]和[6]中的監測時段，因此，可以更好地揭示該地區冰川的時空變化規律。本文綜合波段比值法與人工解譯獲取的1999 年面積與文獻[5]中手動解譯的2000 年面積一致，表明本文解譯的結果精度很高。由于各自采用的分類方法對處于特殊位置的積雪與冰川的辨別能力不同，本文得到的1992—1999 年時間段內的冰川年面積變化率與文獻[6]得到的1992—2002 年時間段內的面積變化率略有差別，但并不影響對冰川變化趨勢的判定[7]。

根據該地區冰川空間分布特點，本文將研究區劃分為A 區和B 區(如圖2(a)所示)。通過計算A 和B 這2 個區域面積參數(見表3)發現：A 和B 這2 個區域冰川均出現退縮，但B 區的冰川比A 區退縮更快。其原因可能是B 區的冰川面積基數更大；A 區3 個時段的冰川面積年平均變化率比較穩定，約為0.45 km2/a，B區則波動較大；1999—2006 年期間退縮速度最快，2006—2009 年期間退縮速度有所減緩。

從整體來看，冰川退縮現象廣泛存在，而不是局限于某幾條冰川，冰川退縮主要集中于冰舌區(消融區)。雪山東側的冰川比西側的退縮更嚴重，局部地形以及冰川中心線的走向可能是導致變化不一致的原因。其中坡度會影響冰川流速，在冰川系統中主要靠冰川運動將積累區的冰雪運送到消融區，每條冰川冰下地形的坡度與冰川流速有直接聯系，進而影響冰川的物質平衡，使得冰川出現不同的退縮狀況；冰川中心線的走向直接影響冰川接收太陽輻射量，而冰川接收太陽輻射量與冰川表面溫度直接相關，進而影響冰川消融速度。另外，冰川大多位于山谷，當冰川兩邊的山峰高度較低時，冰川接收更多的太陽光能量，消融會更加劇烈。冰川的加速消融也導致某些復合冰川的末端已經開始分裂。圖3 所示為圖2(a)中e 冰川的末端在4 個時期的變化情況。對比圖3(a)～(d)中的黑色圓圈部分，很清晰地看出在1992—2009 年間冰川消融加快，末端嚴重退縮，面積不斷減少。圖4 所示為圖2(a)中a，b，c 和d 這4 條冰川的末端在不同時間的位置線。對比這4 條冰川末端位置可以發現：a，b和c 這3 條冰川處于逐年退縮狀態；1992—1999 年，d 冰川末端有所前進，而在1999—2009 年間又退縮，但2009 年的冰川面積仍然比1992 年的大，冰川的加速運動或者某些年份的降水增多可能是導致這種現象的原因。

表3 1992—2009 年格拉丹東雪山A 區和B 區冰川的面積變化Table 3 Glacier area change of part A and part B in Geladandong from 1992 to 2009

3.2 冰前湖變化及分析

當冰川融化加快時，大量的冰融水以及脫落的冰塊就會在冰川的前端、冰面或冰下匯集形成冰川湖，如圖4(c)所示。該冰前湖在1992 年、1999 年、2006年和2009 年這4 個時期內的面積分別為0.19，0.23，0.28 和0.29 km2，呈連續擴大的趨勢。圖5 所示為研究區4 個時期的冰川和冰川湖面積參數。由于研究區能滿足冰川監測要求的遙感數據數量有限，本文只獲取研究區4 個時期的面積。由于實驗數據過少，當利用函數擬合冰川和冰川湖的面積變化時，會導致擬合函數存在很大的誤差，無法統計擬合誤差。但是，通過分析圖5 發現：冰川面積不斷減少，而冰川湖面積在增加，冰川和冰川湖的面積隨時間變化呈負相關，并且冰川面積和冰川湖面積呈反比例函數關系，因此，認為冰前湖面積與冰川消融面積呈非線性正相關。另外，冰川消融和冰川前端的冰塊脫落速度加快是導致冰前湖水量增加的主要因素。夏季溫度升高，冰川消融加強，融水注入冰川湖后使湖水上漲，伸入湖中的冰川前端底部壓力和浮力隨之增大，冰川冰開始脫落形成浮冰，從而使冰川物質轉移至冰川湖中。為了補充冰舌損失的物質，冰體脫落和冰舌消融反過來又會加速冰川運動，使更多的上游物質抵達至此，導致冰川湖面積持續擴大。冰川湖水量激增至冰壩所能承受的最大值后就引發冰湖潰決，冰川洪水會嚴重威脅流域內人們的生命財產安全。另外，湖水增加會影響當地的生態平衡和農牧業，所以，冰川湖變化監測十分重要。

圖3 冰川末端在4 個不同時期的退縮變化Fig.3 Examples of recession of glacier tongue in four different periods

圖4 4 個時期冰川邊界變化(底圖為2009 年TM 影像)Fig.4 Terminus boundary of a glacier(TM images in 2009 is used as base map)

圖5 研究區冰川和冰川湖在1992 年、1999 年、2006 年和2009 年的面積Fig.5 Areas of glaciers and glacial lake in study area in 1992,1999,2006 and 2009

3.3 冰川變化與氣候的關系

冰川的變化與冰川所在地區的氣溫和降水有密切聯系[3-4,17]。本文依據距離格拉丹東雪山最近的沱沱河(34°13′N,92°25′E)氣象站資料[6]，整理1992—2009 年中夏季平均氣溫和夏季降水量，如圖6 所示。從圖6可見：在1992—1999 年、1999—2006 年和2006—2009年這3 個時間段內，該地區夏季月平均氣溫分別為5.3，5.6 和5.8 ℃，夏季降水量分別為230，300 和332 mm。夏季月平均氣溫一直處于升高狀態。在冰川區域，夏季氣溫的升高會導致3 個方面的變化[17]：1) 增加外部環境融化冰和雪的能量；2) 增加液態水的比例，從而減少了冰川積累區濕雪或干雪的體積；3) 降低了冰川表面的反照率，增加冰川吸熱量，加強冰、濕雪和干雪的消融。當冰川周圍環境溫度升高時，會導致冰川內部的溫度緩慢升高，從而引起冰川內部消融。這些變化最終導致冰川表面高程降低，冰川末端的后退和雪線位置升高，使冰川出現負的物質平衡。格拉丹東地區冰川屬于“夏季積累型冰川”，降水主要集中在夏季(5～9 月份)。數據顯示，這3 個時間段內夏季降水量呈現增加的趨勢。在冰川系統中，氣溫主要影響冰川的消融，降水主要影響冰川的積累。若冰川消融量與積累量相等，則冰川的物質平衡不會發生變化，冰川的總面積也不會發生變化。對于格拉丹東冰川，1992—2009 年期間冰川總面積一直減少，冰川一直處于負的物質平衡狀態，說明該冰川的消融量大于積累量。雖然該區域夏季的降水量增加，使得冰川的物質積累增加，但冰川處于負的物質平衡狀態。夏季升溫導致冰川消融量大于積累量，因此，認為夏季降水量增加量并不能彌補夏季氣溫升高造成的消融增加量[17-19]。通過上述冰川的面積變化和氣象數據，認為夏季氣溫的升高可能是導致格拉丹東冰川消融加快、面積減少的主要因素。

圖6 1992—2009年沱沱河氣象站的夏季(5～9月份)平均氣溫和降水量變化Fig.6 Annual mean summer(May to September)temperate and precipitation at Tuotuohe Station in 1992—2009

4 結論

1) 基于1992 年、1999 年、2006 年和2009 年的LandsatTM/ETM+影像，采用波段比值法和人工解譯結合的分類策略，提取格拉丹東地區冰川與冰前湖面積參數并矢量化，然后對當年中多期遙感影像對應的矢量邊界進行綜合評判消除了陰影、云和季節性積雪的影響，以提高分類的精度。1992—2009 年，冰川面積共減少66.68 km2(減小率為7.37%)，其中在1992—1999 年、1999—2006 年和2006—2009 年分別減少19.90 km2(減小率為2.84 km2/a)，39.95 km2(減小率為5.71 km2/a)和6.83 km2(減小率為2.28 km2/a)，處于連續退縮狀態，而冰前湖面積共增加了0.1 km2(增大率為52.6%)，在相應3 個時段內分別增加了0.04 km2(增大率為0.006 km2/a)，0.05 km2(增大率為0.007 km2/a)和0.01 km2(增大率為0.003 km2/a)，處于連續擴張狀態。冰川與冰川湖兩者緊密相連，其面積呈非線性負相關。

2) 得到的1992—1999 年時間段內的冰川年面積變化率與文獻[6]中1992—2002 年時間段內的冰川年面積變化率略有差別，但并不影響對冰川變化趨勢的判定。冰川退縮主要集中于冰舌區(消融區)，研究區域內幾乎所有冰舌都在退縮，但雪山東側的冰舌比西側退縮更嚴重。沱沱河氣象站資料顯示上述3 個時間段內夏季氣溫和降水量都呈現上升的趨勢。對于夏季積累型的格拉丹東冰川來說，夏季降水量增加量并不能彌補夏季持續氣溫上升造成的消融增加量，因此，認為夏季持續升溫可能是導致格拉丹東冰川連續退縮的主要因素。

3) 對于大范圍的冰川變化監測，在冰川末端的冰磧區，還可考慮高精度的DEM 和遙感熱紅外波段信息綜合解譯冰川冰磧區的邊界變化。目前格拉丹東冰川的實地監測資料缺乏，要更深入地了解該地區冰川的物質平衡變化，還需要利用其他手段監測冰川體積變化和冰川的運動變化。在以后的研究中還將使用其他高精度的光學遙感影像、SAR 影像等綜合研究格拉丹東以及其他區域冰川變化情況，以便更全面地掌握冰川與氣候變化之間的關系。

致謝：

感謝美國地質調查局(USGS)及全球土地覆蓋數據數據庫(GLCF)提供LandSat 影像；感謝中國氣象局提供的氣象數據；感謝中科院寒旱所張華偉博士對本文給予的大量幫助。

[1] Paul F, K??b A. Perspectives on the production of a glacier inventory from multispectral satellite data in Arctic Canada：Cumberland Peninsula, Baffin Island[J]. Annals of Glaciology,2005,42：59-66.

[2] Andreassen L M,Paul F,K??b A, et al.The new landsat-derived glacier inventory for Jotunheimen, Norway, and deduced glacier changes since the 1930s[J]. The Cryosphere, 2008, 2(2)：131-145.

[3] 蒲健辰, 姚檀棟, 王寧練, 等.近百年來青藏高原冰川的進退變化[J]. 冰川凍土,2004,26(5)：517-522.PU Jianchen,YAO Tandong,WANG Ninglian,et al.Fluctuations of the glaciers on the Qinghai—Tibetan Plateau during the past century[J]. Journal of Glaciology and Geocryology, 2004, 26(5)：517-522.

[4] 葉慶華, 陳鋒, 姚檀棟, 等. 近30 年來喜馬拉雅山脈西段納木那尼峰地區的冰川變化的遙感監測研究[J]. 遙感學報,2007,11(4)：511-520.YE Qinghua, CHEN Feng,YAO Tandong, et al. Tupu of glacier variations in the Mt. Naimona’nyi Region, Western Himalayas,in the last three decades[J]. Journal of Remote Sensing, 2007,11(4)：511-520.

[5] 魯安新, 姚檀棟, 劉時銀, 等. 青藏高原格拉丹東地區冰川變化的遙感監測[J]. 冰川凍土,2002,24(5)：559-562.LU Anxin, YAO Tandong, LIU Shiyin, et al. Glacier change in the Geladandong area of the Tibetan Plateau monitored by remote sensing[J]. Journal of Glaciology Geocryology, 2002,24(5)：559-562.

[6] Ye Q H, Kang S C, Chen F, et al. Monitoring glacier variations on Mt. Geladandong, Central Tibetan Plateau,from 1969 to 2002 using remote sensing and GIS technologies[J]. Journal of Glaciology,2006,52(179)：537-545.

[7] Bolch T,Yao T,Kang S,et al.A glacier inventory for the western Nyainqentanglha Range and the Nam CoBasin,Tibet,and glacier changes 1976—2009[J].The Cryosphere,2010,4(3)：419-433.

[8] Bolch T, Menounos B, Wheate R. Landsat-based inventory of glaciers in western Canada, 1985—2005[J]. Remote Sensing of Environment,2010,114(1)：127-137.

[9] Paul F, K??b A, Maisch M, et al. The new remote-sensing derived glacier inventory： I Methods[J]. Annals of Glaciology,2002,34：355-361.

[10] Racoviteanu A E, Paul F, Raup B, et al.Challenges and recommendations in mapping of glacier parameters from space：results of the 2008 Global Land Ice Measurements from Space(GLIMS) workshop, Boulder, Colorado, USA[J]. Annals of Glaciology,2009,50：53-69.

[11] Gjermundsen E, Mathieu R, Kaab A, et al. Assessment of multispectral glacier mapping methods and derivation of glacier area changes, 1978—2002, in the central Southern Alps, New Zealand, from ASTER satellite data, field survey and existing inventory data[J]. Journal of Glaciology, 2011, 57(204)：667-683.

[12] Paul F.Evaluation of different methods of glacier mapping using landsat TM[C]// Proceedings of EARSeL SIG-Work-shop Land Ice and Snow.Dresden,2000：16-17.

[13] Raup B, K??b A, Kargel J S, et al. Remote sensing and GIS technology in the global land ice measurements from space(GLIMS)Project[J].Computes&Geosciences,2007,33(1)：104-125.

[14] 宋波, 何元慶, 龐洪喜, 等. 基于遙感和GIS 的我國季風海洋型冰川區冰磧物覆蓋型冰川邊界的自動識別[J]. 冰川凍土,2007,29(3)：456-462.SONG Bo, HE Yuanqing, PANG Hongxi, et al. Identifying automatically the debris-covered glaciers in China’s monsoonal temperate-glacier regions based on remote sensing and GIS[J].Journal of Glaciology and Geocryology,2007,29(3)：456-462.

[15] 王建. 衛星遙感積雪制圖方法對比與分析[J]. 遙感技術與應用,1999,14(4)：29-36.WANG Jian. Comparison and analysis on methods of snow cover mapping by using satellite remote sensing data[J]. Remote Sensing Technology and Application,1999,14(4)：29-36.

[16] 張世強, 盧健, 劉時銀. 利用TM 高光譜圖像提取青藏高原喀喇昆侖山區現代冰川邊界[J]. 武漢大學學報(信息科學版),2001,26(5)：435-440.ZHANG Shiqiang, LU Jian, LIU Shiyin. Deriving glacier border information on Qinghai Tibet by TM high spectrum image[J].Geomatics and Information Science of Wuhan University, 2001,26(5)：435-440.

[17] Narama C, K??b A, Duishonakunov M, et al. Spatial variability of recent glacier area changes in the Tien Shan Mountains,Central Asia,using Corona(～1970),Landsat(～2000),and ALOS(～2007) satellite data[J]. Global and Planetary Change, 2010,71(1/2)：42-54.

[18] Khromova T E, Osipova G B, Tsvetkov D G, et al. Changes in glacier extent in the eastern Pamir, Central Asia, determined from historical data and ASTER imagery[J]. Remote Sensing of Environment,2006,102(1/2)：24-32.

[19] Paul F, Huggel C, K??b A. Combining satellite multispectral image data and a digital elevation model for mapping debris-covered glaciers[J]. Remote Sensing of Environment,2004,89(4)：510-518.