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

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

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

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

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

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

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

1.1 研究區(qū)域

Amery冰架位于南大洋普里茲灣（Prydz Bay）的南部，是南極三大冰架之一，而澳大利亞的莫森站（Mawson Station）、戴維斯站（Davis Station）、俄羅斯的進(jìn)度站（Progress Station）和中國的中山站都設(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（極地冰蓋測高數(shù)據(jù)）激光測高數(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á)測高、機(jī)載冰雷達(dá)、地面觀測和地理定位系統(tǒng)（GPS）等多種觀測數(shù)據(jù)，采用多種內(nèi)插方法得到的最高分辨率為200 m 的全南極DEM.

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

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

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

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

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

3 實驗與分析

本文選用了15期數(shù)據(jù)，并合并了第15和16兩期獲取時間較近的數(shù)據(jù)，共計14期數(shù)據(jù).首先分析高程的年際變化，結(jié)果如圖3所示；然后依次對相鄰兩期測高數(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)行高程變化比較的時間段，縱坐標(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

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

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

除了上述總體平均高程變化分析，本文還對其進(jìn)行了按高程分層次分析，即分析不同高程范圍內(nèi)在一段時間內(nèi)的高程變化，具體結(jié)果如圖5所示.以小方塊標(biāo)記的點是某一高程范圍內(nèi)的平均高程變化，有極個別處于高程較低區(qū)域因為數(shù)據(jù)較少而且變化比較劇烈，所以變化值超出了繪圖范圍.細(xì)實線表示該高程區(qū)間的高程變化波動范圍，從而可以直觀地顯示該高程區(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ū)確實以負(fù)增長為主，其原因有待解釋.2009年3月～9月的重復(fù)軌投影點數(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，但相鄰的兩個年內(nèi)高程變化都比一般情況要高，因此可以認(rèn)為2005年～2006年期間該區(qū)域的總體物質(zhì)平衡為正平衡.唯一看似異常的是圖5m 中表示的2008年10月至2009年3月的高程變化，而實際上該結(jié)果是由第15～17期數(shù)據(jù)計算得到的，并且把第15和16期數(shù)據(jù)作了合并處理，其數(shù)據(jù)獲取時間分別為2008－10－04～2008－10－19和2008－11－25～2008－12－17，其數(shù)據(jù)點個數(shù)分別為185 290和257 014.其中第16期數(shù)據(jù)的采集時間是南極氣溫最高的時間，而第17期數(shù)據(jù)的采集時間為2009－03－09～2009－04－11，即將進(jìn)入最寒冷的季節(jié)，因此更傾向于代表從夏季到冬季的冰蓋表面高程變化，所以結(jié)果顯示為正積累也是合理的.

4 結(jié)論與展望

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

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

［1］ 胡明城.關(guān)于海面的長期上升趨勢［J］.測繪學(xué)報，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交叉點分析探測恩德比地冰蓋近年高程變化［J］.武漢大學(xué)學(xué)報：信息科學(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］.地域研究與開發(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］ 鄂棟臣，徐瑩，張小紅.星載激光測高及其在極地的應(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.