謝 歡，顧振雄，劉 俊，童小華

（1.同濟(jì)大學(xué) 空間信息科學(xué)及可持續(xù)發(fā)展應(yīng)用中心，上海200092；2.同濟(jì)大學(xué) 測(cè)繪與地理信息學(xué)院，上海200092）

南極冰雪的動(dòng)態(tài)平衡對(duì)于維持全球氣候穩(wěn)定極其重要，全球氣候的變化也會(huì)對(duì)南極冰雪平衡產(chǎn)生巨大影響，其中最直接的可觀測(cè)影響之一就是南極冰蓋表面的高程以及水平位置變化.包括南極冰川在內(nèi)的冰雪消融是引起海平面上升的因素之一［1］，Amery冰架是南極三大冰架之一，臨近我國(guó)南極科考中山站，對(duì)于研究東南極物質(zhì)平衡有著重大的意義［2］.

ICESat衛(wèi)星發(fā)射的的科學(xué)目的之一是監(jiān)測(cè)南極和格陵蘭島冰蓋的高程變化［3］，對(duì)其激光測(cè)高數(shù)據(jù)常用的分析方法包括重復(fù)點(diǎn)分析、交叉點(diǎn)分析和重復(fù)軌分析等.Zwally［4］用交叉點(diǎn)法對(duì)格陵蘭島區(qū)域的Geosat和Seasat雷達(dá)測(cè)高數(shù)據(jù)進(jìn)行高程變化分 析，Brenner 等［5］用 交 叉 點(diǎn) 法 對(duì) 極 地ERS－2、Envisat雷達(dá)測(cè)高數(shù)據(jù)和ICESat激光測(cè)高數(shù)據(jù)進(jìn)行高程變化分析，史紅嶺等［6］人也曾用ICESat交叉點(diǎn)分析恩德比地區(qū)的高程變化.

雖然重復(fù)點(diǎn)法和交叉點(diǎn)法能夠非常精確地分析出某一位置的高程變化，但是它們對(duì)于ICESat數(shù)據(jù)的利用率太低，一般都低于1%［7］.為了提高數(shù)據(jù)的利用率，有學(xué)者提出了重復(fù)軌分析法.Abdalati和Bamber等［8］早在1995 至2002 年間就已分別在加拿大和挪威用機(jī)載激光測(cè)高儀對(duì)局部區(qū)域進(jìn)行了高程測(cè)量，并通過(guò)比較前后兩次的高差來(lái)推算局部區(qū)域的質(zhì)量平衡.重復(fù)軌法還被用于探測(cè)西南極Ronne冰架接地線、探測(cè)Whillans冰流等地區(qū)的冰下湖以及估算格陵蘭和南極冰蓋的高程變化率［8］.Fricker等［9］用ICESat重復(fù)軌數(shù)據(jù)，基于潮汐對(duì)延伸到海洋的冰架具有抬升作用的原理來(lái)探測(cè)South Ronne冰架的接地線，以及基于類(lèi)似的原理來(lái)探測(cè)冰下湖［10］.

然而這些研究普遍沒(méi)有考慮重復(fù)軌道之間的坡度變化對(duì)高程變化的影響.在劉學(xué)軍等人對(duì)多種計(jì)算任意方向數(shù)字高程模型（digital elevation model，DEM）坡度的方法研究的基礎(chǔ)上［11］，本文使用基于DEM 坡度改正的重復(fù)軌分析法［7］，分析Amery 冰架及附近區(qū)域2003～2009年的年際高程變化，并對(duì)不同海拔區(qū)域的冰蓋表面高程變化進(jìn)行了分析.

1 研究區(qū)域與數(shù)據(jù)

1.1 研究區(qū)域

Amery冰架位于南大洋普里茲灣（Prydz Bay）的南部，是南極三大冰架之一，而澳大利亞的莫森站（Mawson Station）、戴維斯站（Davis Station）、俄羅斯的進(jìn)度站（Progress Station）和中國(guó)的中山站都設(shè)立在Amery冰架附近.本文的研究區(qū)域定在Amery冰架及其附近經(jīng)緯度范圍為40°E～95°E、68.5°S～81°S的扇形區(qū)域，如圖1所示.

圖1 研究區(qū)域Fig.1 Study Area

1.2 ICESat數(shù)據(jù)介紹

本文中使用的數(shù)據(jù)包括ICESat GLA 12（極地冰蓋測(cè)高數(shù)據(jù)）激光測(cè)高數(shù)據(jù)［12］和南極DEM 數(shù)據(jù)，ICEsat數(shù)據(jù)用于高程的變化分析［13］，DEM 數(shù)據(jù)用于ICESat重復(fù)軌之間的坡度改正.

1.3 DEM 數(shù)據(jù)介紹

本文采用的是Liu等［14］用不同的內(nèi)插方法得到的200m 分辨率的南極DEM，一般記為RAMPv2，如圖1所示.它是基于紙質(zhì)地形圖、雷達(dá)測(cè)高、機(jī)載冰雷達(dá)、地面觀測(cè)和地理定位系統(tǒng)（GPS）等多種觀測(cè)數(shù)據(jù)，采用多種內(nèi)插方法得到的最高分辨率為200 m 的全南極DEM.

2 基于DEM 坡度改正的ICESat重復(fù)軌分析

重復(fù)軌分析法流程圖如圖2a所示.圖中，x，y分別為投影平面的橫坐標(biāo)，縱坐標(biāo)；B，L分別為緯度，徑度；H為大地高；T為時(shí)間.首先對(duì)ICESat激光測(cè)高點(diǎn)進(jìn)行過(guò)濾［15］、橢球轉(zhuǎn)換、投影變換等預(yù)處理和升降軌分離，然后在重復(fù)軌上計(jì)算臨近點(diǎn)的高程投影，并利用DEM 進(jìn)行坡度改正，最后得到高程變化量.

本文使用二次曲線對(duì)分離后的每一條軌道進(jìn)行最小二乘擬合，然后在升軌和降軌中分別搜索重復(fù)軌.重復(fù)軌組中重復(fù)軌的數(shù)量一般是2條或2條以上.

本文研究區(qū)域的ICESat重復(fù)軌之間的平均距離為100m，而該區(qū)域的平均坡度為0.3%，因此必須考慮坡度對(duì)重復(fù)軌之間高程點(diǎn)投影時(shí)的影響后才能進(jìn)行高程變化分析.如圖2b所示，以圖中軌道C上的C2點(diǎn)為例，首先計(jì)算C2到軌道D上所有點(diǎn)的距離，如果最小的兩個(gè)距離都小于某個(gè)給定閾值，則認(rèn)為軌道D上存在相鄰的激光測(cè)高點(diǎn)D1和D2，此時(shí)可以獲取C2到D1和D2連線上的垂足點(diǎn)DREF2，然后分別計(jì)算投影點(diǎn)C2投影到軌道D后的高程HCREF2＝HC2＋dHDEM2和由D1、D2得到的內(nèi)插高程HDREF2，其中HC2為C2點(diǎn)的高程，dHDEM2即為根據(jù)已有DEM 計(jì)算的坡度而得到的投影高程改正［11］.那么HCREF2與HDREF2之間的差別即可代表該投影點(diǎn)處兩個(gè)時(shí)期之間的高程變化，進(jìn)而可以用于分析整個(gè)研究區(qū)域的高程變化［7］.

在確定了投影點(diǎn)C2和內(nèi)插點(diǎn)DREF2之后，首先計(jì)算它們的中點(diǎn)在DEM 上兩個(gè)正交方向的偏導(dǎo)數(shù)，再根據(jù)它們連線的方位角計(jì)算方向?qū)?shù)，即該方向的坡度［11］.

3 實(shí)驗(yàn)與分析

本文選用了15期數(shù)據(jù)，并合并了第15和16兩期獲取時(shí)間較近的數(shù)據(jù)，共計(jì)14期數(shù)據(jù).首先分析高程的年際變化，結(jié)果如圖3所示；然后依次對(duì)相鄰兩期測(cè)高數(shù)據(jù)進(jìn)行分析，總體結(jié)果如圖4所示，具體高程變化分析如圖5所示.

3.1 總體高程變化分析

圖3顯示的是研究區(qū)域從2003年～2009年相鄰年度2月和10月的平均高程變化結(jié)果，圖4顯示的是該區(qū)域從2003～2009年連續(xù)的平均高程變化結(jié)果.圖中的橫坐標(biāo)均代表進(jìn)行高程變化比較的時(shí)間段，縱坐標(biāo)為研究區(qū)域以cm 為單位的平均高程變化.圖3表明，本文研究區(qū)域的平均高程變化一般在－3.2～6.2cm，只有2004～2005年的高程變化很大，為10cm 以上.

圖5 2003年2月～2009年3月高程變化Fig.5 Elevation change from February，2003to March，2009

類(lèi)似地，圖4表明，在2005年2～10月期間，研究區(qū)域的平均高程變化為上升10.6cm.這段時(shí)間正好覆蓋了南極的整個(gè)冬季，從2月冬季剛剛開(kāi)始到10月冬季即將結(jié)束，所以在寒冷的氣溫作為主導(dǎo)因素的作用下，研究區(qū)域的物質(zhì)平衡有顯著的正積累.圖4中的2006年2～10月、2007年3～10月、2008年2～10月有同樣的結(jié)果，只是量級(jí)上略有差別.至于2003年2～9月和2004年2～10月的變化可能是由于其他因素代替氣溫因素占主導(dǎo)而引起的.又如圖4中2005年10月至2006年2月這段時(shí)間覆蓋的是南極夏季開(kāi)始到夏季結(jié)束，是南極一年中氣溫相對(duì)最高的時(shí)間，物質(zhì)平衡應(yīng)該是負(fù)積累，與本文的結(jié)果上 升1.8cm 不 一 致.但 是 將2005 年10 月 至2006年2月這段時(shí)間的高程變化和與其相鄰的兩個(gè)結(jié)果放在一起分析，即2005年2～10 月上升10.6 cm 和2006年2～10月上升5.2cm 這兩個(gè)結(jié)果，可以看到從2005年10月至2006年10月的高程變化都以顯著上升為主，因此可推測(cè)這段時(shí)期的物質(zhì)平衡以正積累為主.

3.2 不同高程范圍的高程變化分析

除了上述總體平均高程變化分析，本文還對(duì)其進(jìn)行了按高程分層次分析，即分析不同高程范圍內(nèi)在一段時(shí)間內(nèi)的高程變化，具體結(jié)果如圖5所示.以小方塊標(biāo)記的點(diǎn)是某一高程范圍內(nèi)的平均高程變化，有極個(gè)別處于高程較低區(qū)域因?yàn)閿?shù)據(jù)較少而且變化比較劇烈，所以變化值超出了繪圖范圍.細(xì)實(shí)線表示該高程區(qū)間的高程變化波動(dòng)范圍，從而可以直觀地顯示該高程區(qū)間高程變化強(qiáng)度.從圖5a—5m中，首先可以看到Amery冰架及其附近區(qū)域海拔較低的沿海區(qū)域的高程變化比海拔較高的內(nèi)陸區(qū)域更加劇烈.

圖5c—5f表明，2005～2008年期間同一年內(nèi)的夏季到冬季，大多數(shù)區(qū)域的高程變化是正值，并且與圖4中的總體結(jié)果一致.在圖4 中2003 年和2004年的年內(nèi)總體高程變化為負(fù)值，從圖5a—5b中更細(xì)化的結(jié)果中也可以看到高程在2 000m 以上的內(nèi)陸地區(qū)確實(shí)以負(fù)增長(zhǎng)為主，其原因有待解釋.2009年3月～9月的重復(fù)軌投影點(diǎn)數(shù)據(jù)較少，約為其他年份的1／3.所以雖然圖5g顯示2009年3月～9月期間低海拔地區(qū)的震蕩較大，但由于高海拔地區(qū)的數(shù)據(jù)量占主導(dǎo)并且變化趨近于零，所以最終總體高程變化結(jié)果趨于穩(wěn)定，與圖4中的總體結(jié)果一致.

圖5h—5l表明，2003 年～2008 年期間冬季到夏季，大多數(shù)區(qū)域的高程變化是負(fù)值.從圖4中也可以得到相同的結(jié)論，雖然2005年10月至2006年2月期間的總體高程變化為1.8cm，但相鄰的兩個(gè)年內(nèi)高程變化都比一般情況要高，因此可以認(rèn)為2005年～2006年期間該區(qū)域的總體物質(zhì)平衡為正平衡.唯一看似異常的是圖5m 中表示的2008年10月至2009年3月的高程變化，而實(shí)際上該結(jié)果是由第15～17期數(shù)據(jù)計(jì)算得到的，并且把第15和16期數(shù)據(jù)作了合并處理，其數(shù)據(jù)獲取時(shí)間分別為2008－10－04～2008－10－19和2008－11－25～2008－12－17，其數(shù)據(jù)點(diǎn)個(gè)數(shù)分別為185 290和257 014.其中第16期數(shù)據(jù)的采集時(shí)間是南極氣溫最高的時(shí)間，而第17期數(shù)據(jù)的采集時(shí)間為2009－03－09～2009－04－11，即將進(jìn)入最寒冷的季節(jié)，因此更傾向于代表從夏季到冬季的冰蓋表面高程變化，所以結(jié)果顯示為正積累也是合理的.

4 結(jié)論與展望

本文用重復(fù)軌分析法對(duì)Amery 冰架及其附近區(qū)域的ICESat測(cè)高數(shù)據(jù)進(jìn)行了冰蓋表面高程變化分析.從18期數(shù)據(jù)中選取了15期數(shù)據(jù)進(jìn)行重復(fù)軌分析，得到了12個(gè)年際高程變化結(jié)果和13 個(gè)相鄰季節(jié)的冰蓋表面高程變化比較結(jié)果.此外，還對(duì)13個(gè)相鄰季節(jié)之間的冰蓋表面高程變化進(jìn)行了按不同高程區(qū)間的分析.海拔較低即靠近海洋的區(qū)域高程變化劇烈，平均高程變化波動(dòng)范圍為±40cm，而海拔較高即內(nèi)陸區(qū)域的高程變化比較緩和，平均高程變化均在±10cm 以?xún)?nèi).冰流從冰原流向冰架，海拔由高到低，因此這些區(qū)域往往伴隨著較大的高程變化.冰原通常比較穩(wěn)定，因此高程變化幅度不會(huì)太大.總體上，從夏季到下一個(gè)冬季的高程變化以增長(zhǎng)為主，從冬季到下一個(gè)夏季的高程變化以降低為主.

由于極地環(huán)境惡劣，地面真實(shí)數(shù)據(jù)稀少，只有將多種獨(dú)立的結(jié)果進(jìn)行比對(duì)，才能更好地證明方法的可靠性.下一步工作將是采用其他數(shù)據(jù)源，如ERS（歐洲遙感衛(wèi)星）雷達(dá)測(cè)高數(shù)據(jù)、ERS InSAR 數(shù)據(jù)對(duì)同一地區(qū)同一時(shí)期和不同時(shí)期的高程變化進(jìn)行分析.另外如果在本文中重復(fù)軌投影時(shí)使用多個(gè)時(shí)期的、精度、分辨率更高的DEM，必定會(huì)使坡度改正更加精確.

［1］ 胡明城.關(guān)于海面的長(zhǎng)期上升趨勢(shì)［J］.測(cè)繪學(xué)報(bào)，1992，21（2）：151.HU Mingcheng.On the global long－term sea level rise trends［J］.ActaGeodaeticaEtCartographicaSinica，1992，21（2）：151.

［2］ Williams M J M，Warner R C，Budd W F.Sensitivity of the amery ice shelf，Antarctica，to changes in the climate of the southern ocean［J］.Journal of Climate，2002，15（19）：2740.

［3］ Schutz B E，Zwally H J，Shuman C A，et al.Overview of the ICESat mission［J］.Geophysical Research Letters，2005，32（21）：L21S01.

［4］ Zwally H J.Growth of greenland ice sheet：measurement［J］.Science，1989，246（4937）：1587.

［5］ Brenner A C，DiMarzio J P，Zwally H J.Precision and accuracy of satellite radar and laser altimeter data over the continental ice sheets［J］.IEEE Transactions on Geoscience and Remote Sensing，2007，45（2）：321.

［6］ 史紅嶺，陸洋，鮑李峰，等.利用ICESat交叉點(diǎn)分析探測(cè)恩德比地冰蓋近年高程變化［J］.武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2009，34（4）：440.SHI Hongling，LU Yang，BAO Lifeng，et al.Recent elevation change detection of enderby land ice sheet using ICESat crossover analysis［J］.Geomatics and Information Science of Wuhan University，2009，34（4）：440.

［7］ Moholdta G，Nutha C，Hagena J O，et al.Recent elevation changes of Svalbard glaciers derived from ICESat laser altimetry ［J］.Remote Sensing of Environment，2010，114（11）：2756.

［8］ Bamber J L，Krabill W，Raper V B，et al.Elevation changes measured on Svalbard glaciers and ice caps from airborne laser data［J］.Annals of Glaciology，2005，42（1）：202.

［9］ Fricker H A，Padman L.Ice shelf grounding zone structure from ICESat laser altimetry ［J］.Geophysical Research Letters，2006，33（15）：L15502.

［10］ Fricker H A，Scambos T，Bindschadler R，et al.An active subglacial water system in West Antarctica mapped from space［J］.Science，2007，315（5818）：1544.

［11］ 劉學(xué)軍，任志峰，王彥芳，等.基于DEM 的任意方向坡度計(jì)算方法［J］.地域研究與開(kāi)發(fā)，2009，28（4）：139.LIU Xuejun，REN Zhifeng，WANG Yanfang，et al.Slope model at arbitrary direction derived from grid－based DEM ［J］.Areal Research and Development，2009，28（4）：139.

［12］ Zwally H J，Schutz R，Bentley C，et al.GLAS／ICESat L2 Antarctic and Greenland ice sheet altimetry data V33 ［GLA 12］［DB／OL］.［2012－02－14］.http：／／nsidc.org／data／docs／daac／glas＿icesat＿l1＿l2＿global＿altimetry.gd.html.

［13］ 鄂棟臣，徐瑩，張小紅.星載激光測(cè)高及其在極地的應(yīng)用研究分析［J］.極地研究，2006，18（2）：148.E Dongchen，XU Ying，ZHANG Xiaohong.Preliminary study of application of satellite laser altimetry technology in polar region［J］.Advances in Polar Science，2006，18（2）：148.

［14］ Liu H X，Jezek K C，Li B Y.Development of an antarctic digital elevation model by integrating cartographic and remotely sensed data：ageographic information system based approach［J］.Journal of Geophysical Research，1999，104（B10）：199.

［15］ Banber J L，Gomez－Dans J L，Griggs J A.A new1 km digital elevation model of the Antarctic derived from combined satellite radar and laser data—part 1：data and methods［J］.The Cryosphere，2009，3（1）：101.