陳國棟，梁圣豪，孟子淇，朱家亨

（蘇州科技大學(xué) 地理科學(xué)與測繪工程學(xué)院，江蘇 蘇州 215009）

隨著人類社會的發(fā)展，人類活動對地球氣候產(chǎn)生了巨大影響，其中，海平面上升是最受關(guān)注的環(huán)境問題之一，嚴(yán)重影響沿海地區(qū)的生產(chǎn)生活安全。格陵蘭冰蓋是世界第二大陸地冰蓋，儲存了全球十分之一的冰量，完全融化可以使全球海平面上升約6.5 m。進(jìn)入21世紀(jì)以來，格陵蘭冰蓋消融的速度明顯增加，對海平面上升的貢獻(xiàn)已經(jīng)遠(yuǎn)遠(yuǎn)超過南極冰蓋[1]，因此，對格陵蘭冰蓋冰雪消融狀況進(jìn)行大范圍、長時間的監(jiān)測是十分必要的。

由于冰蓋地理位置和自然環(huán)境特殊，一般的地面沉降手段[2]在冰雪覆蓋的極地條件下難以實施，因此，極地冰蓋高程變化的監(jiān)測主要通過衛(wèi)星測高技術(shù)實現(xiàn)，Envisat、ICESat等多顆測高衛(wèi)星的數(shù)據(jù)都在格陵蘭冰蓋的監(jiān)測中貢獻(xiàn)了許多重要成果[3]。Cryosat-2是由歐空局（European Space Agency，ESA）研制的雷達(dá)測高衛(wèi)星，于2010年4月8日發(fā)射，相比于過去的雷達(dá)測高衛(wèi)星，Cryosat-2采用了新一代的雷達(dá)高度計SIRAL，不僅大大提高了沿軌分辨率，還具有三種不同的測量模式，以應(yīng)對不同特點的地表類型[4]。因此，Cryosat-2衛(wèi)星測高數(shù)據(jù)在兩極海冰厚度監(jiān)測、冰蓋及冰川質(zhì)量變化、重力場研究等方面都取得了令人矚目的成果[5-7]。近年來，國內(nèi)也有不少學(xué)者利用這顆衛(wèi)星的觀測數(shù)據(jù)開展了一系列的研究[8-11]。

筆者利用2011年1月至2017年12月間Cryosat-2衛(wèi)星測高數(shù)據(jù)對格陵蘭冰蓋高程及體積變化趨勢進(jìn)行了確定，詳細(xì)分析了冰蓋各流域的變化情況，并將結(jié)果與2003—2009年間ICESat測高數(shù)據(jù)獲得的冰蓋變化情況進(jìn)行了對比，結(jié)果表明格陵蘭冰蓋高程下降趨勢正在加劇。

1 數(shù)據(jù)介紹及預(yù)處理

Cryosat-2衛(wèi)星是ESA地球探測計劃的組成任務(wù)之一，主要科學(xué)目標(biāo)是為了監(jiān)測地球海洋冰蓋和大陸冰蓋厚度的變化[3]。為了更好地應(yīng)對不同地表形態(tài)，Cryosat-2搭載的新一代雷達(dá)高度計SIRAL采用了三種測量模式。其中在格陵蘭冰蓋表面采用了兩種觀測模式：在地形平坦的格陵蘭島內(nèi)陸地區(qū)采用低分辨率模式（Low Resolution Mode，LRM）；在地形復(fù)雜的沿海區(qū)域則采用合成孔徑雷達(dá)干涉測量模式（Interferometry Synthetic Aperture Radar，InSAR）。獲得的數(shù)據(jù)經(jīng)過不同類型的處理后，分為Level 1和Level 2兩種級別，各級數(shù)據(jù)按照模式和處理方式及包含信息的不同又分為23類產(chǎn)品。筆者采用的數(shù)據(jù)為融合處理后的SIR_GDR測高數(shù)據(jù)，使用的數(shù)據(jù)時間跨度為2011年1月至2017年12月。該數(shù)據(jù)在處理過程中已經(jīng)過了粗差剔除、波形重跟蹤處理等各項誤差改正，并將三種測量模式的結(jié)果以統(tǒng)一的形式進(jìn)行存儲，數(shù)據(jù)的坐標(biāo)系統(tǒng)為WGS-84系統(tǒng)，時間系統(tǒng)為國際原子時。

為避免海洋、無冰覆蓋的陸地等地表類型上的數(shù)據(jù)對格陵蘭冰蓋高程變化的計算產(chǎn)生影響，文中利用1 km分辨率格陵蘭島地表類型格網(wǎng)[12]對數(shù)據(jù)進(jìn)行了篩選，只保留冰蓋和冰川上的測高數(shù)據(jù)。此外，SIR_GDR產(chǎn)品中仍然存在一些粗差數(shù)據(jù)，為保證參與計算的觀測值的質(zhì)量，根據(jù)產(chǎn)品中給出的Height_error、SIN_x-track_angle_error和Calibration_Warning三項質(zhì)量標(biāo)記將相應(yīng)的粗差數(shù)據(jù)進(jìn)行剔除[13]。完成數(shù)據(jù)編輯后的格陵蘭冰蓋表面Cryosat-2測高數(shù)據(jù)分布如圖1所示。

圖1 LRM模式（左）和SARin模式（右）數(shù)據(jù)分布

2 確定冰蓋表面高程和體積變化的方法

在利用衛(wèi)星測高數(shù)據(jù)對冰蓋變化的研究中，常用的方法包括交叉點法和重復(fù)軌道法，其中交叉點計算的高程變化精度較高，但交叉點分布稀疏，數(shù)據(jù)利用率低[14]，而重復(fù)軌道法雖然數(shù)據(jù)利用率高，但在計算時需要考慮地面坡度的影響，結(jié)果精度不高[15]。因此，為了能夠更好地確定格陵蘭冰蓋表面高程的變化，文中采用了一種以地面格網(wǎng)為單位，在格網(wǎng)內(nèi)平差確定地面高程變化的方法，將Cryosat-2數(shù)據(jù)中的大地坐標(biāo)以切圓錐投影方式投影到平面，切緯度為71°N，中央經(jīng)度為39°W，并將格陵蘭區(qū)域分割成5 km×5 km分辨率的格網(wǎng)。圖2為一個格網(wǎng)的示意圖，其中的黑點表示Cryosat-2測高數(shù)據(jù)的平面位置。設(shè)格網(wǎng)中心點在起始時刻t0時的高程為h0，而每個Cryosat-2數(shù)據(jù)點高程hi與格網(wǎng)中心點起始高程h0之差包括兩部分：一部分是由地面坡度及平面位置差異引起的空間上的高程差異；另一部分則是由于冰蓋消融等原因造成的時間上的高程變化。因此，每個Cryosat-2測高值可表示為以下函數(shù)

圖2 格網(wǎng)內(nèi)測高點分布示意圖

式中，Δx、Δy分別表示數(shù)據(jù)點與格網(wǎng)中心點在南北和東西方向上的坐標(biāo)差，dh/dx和dh/dy分別表示由坡度引起的高程變化沿x軸和y軸方向上的分量，t表示數(shù)據(jù)采集的時間，dh/dt是高程隨時間的線性變化率。

若每個格網(wǎng)內(nèi)有n個測高點，則每個測高點都可以列出一個形如式（1）的方程，此時將方程改為誤差方程，將每個誤差方程聯(lián)立成方程組

由圖1可知，每個5 km×5 km格網(wǎng)內(nèi)的測高數(shù)據(jù)個數(shù)約為數(shù)百個甚至上千個，顯然式（2）可以通過最小二乘原理獲得最或然值，其中dh/dt的估值即為該格網(wǎng)內(nèi)冰蓋表面高程的變化速率。

為了進(jìn)一步避免觀測值粗差對平差結(jié)果造成影響，文中對平差后獲得的觀測值改正數(shù)按照三倍中誤差的限差進(jìn)行檢驗，在剔除對應(yīng)的觀測值后重新進(jìn)行平差計算，如此迭代，直到所有的改正值都在限差范圍內(nèi)，即可輸出最終結(jié)果。獲得每個格網(wǎng)的高程變化后，將其乘以格網(wǎng)投影面積，即可求得格網(wǎng)內(nèi)的體積變化量，在一定范圍內(nèi)求和就可以獲得相應(yīng)區(qū)域內(nèi)的冰蓋體積變化。

3 結(jié)果及分析

文中利用Cryosat-2測高數(shù)據(jù)計算了2011年1月至2017年12月間格陵蘭冰蓋表面高程變化，并按照5 km×5 km分辨率進(jìn)行格網(wǎng)化，結(jié)果如圖3所示。參考國外相關(guān)研究[16]，采用自助抽樣法對精度進(jìn)行估計后，得到整個格陵蘭地區(qū)冰蓋、冰川的平均高程變化為（-12.3±0.38）cm·a－1，體積變化速率為（-220.38±6.0）km3·a－1。

格陵蘭東南沿海地區(qū)冰川密集，是冰蓋消融現(xiàn)象嚴(yán)重的區(qū)域[17－18]，但由于該區(qū)域地形起伏較大[19]，Cryosat-2數(shù)據(jù)質(zhì)量不高，導(dǎo)致這一區(qū)域數(shù)據(jù)稀疏（如圖1所示），因此，圖3結(jié)果在該區(qū)域出現(xiàn)較多空白，無法體現(xiàn)當(dāng)?shù)乇魇У那闆r，這一問題目前還沒有較好的解決方案[19]，也是今后研究的重點之一。除東南沿海地區(qū)以外，文中在格陵蘭冰蓋的絕大多數(shù)區(qū)域都獲得了有效的冰蓋高程變化結(jié)果。由圖3可知，格陵蘭冰蓋表面高程呈現(xiàn)上升或基本維持不變的區(qū)域主要集中在海拔較高的冰蓋內(nèi)陸地區(qū)；而體積流失現(xiàn)象則主要發(fā)生在冰蓋邊緣，特別是西北沿海、東南沿海以及雅各布港等冰川密集區(qū)域。

圖3 格陵蘭冰蓋高程變化示意圖

根據(jù)冰川流動情況，格陵蘭冰蓋可以劃分為8個冰川流域系統(tǒng)，每個流域又分割為若干個子流域[20]（如圖3中灰色實線所示），為更好地對格陵蘭冰蓋進(jìn)行研究和分析，文中結(jié)合格陵蘭冰蓋冰川流域系統(tǒng)劃分?jǐn)?shù)據(jù)對各流域的變化情況進(jìn)行了分析，其結(jié)果見表1。由表1可知，除去東南沿海數(shù)據(jù)缺失較多的DS3和DS4流域以外，在格陵蘭島的西部DS6、DS7、DS8流域是冰蓋體積流失較多的地區(qū)，也是平均高程變化最大的地區(qū)，下降速度超過20 cm·a－1。而在北部地區(qū)，高程變化相較DS6、DS7、DS8這三個流域來說，冰蓋消融較少，平均高程變化也比較小。

表1 2011—2017年格陵蘭島各流域冰蓋體積變化表

文獻(xiàn)[17]曾利用ICESat激光測高衛(wèi)星對2003—2009年間格陵蘭冰蓋的體積變化進(jìn)行過詳細(xì)的分析，將其結(jié)果與文中獲得的2011—2017年間的體積變化進(jìn)行對比，結(jié)果如圖4所示，由于Cryosat-2在東南沿海區(qū)域有較多數(shù)據(jù)缺失，故DS3、DS4、DS5三個流域未進(jìn)行對比。數(shù)據(jù)表明，相比于2003—2009年的變化情況，2010年后的DS1流域中冰蓋消融速度有所增長，但增長速度并不明顯；而在DS6、DS7流域中，冰蓋消融速度迅速增加，冰川流失情況更加嚴(yán)重。在DS2流域中，2003—2009年冰蓋體積變化為體積增加，而2010年后，冰蓋體積則在減少，即表明體積損失的區(qū)域進(jìn)一步擴大，并延伸至氣候寒冷的格陵蘭北部。

圖4 Cryosat-2與ICESat獲得的格陵蘭冰蓋體積變化對比（單位：km3·a－1）

位于DS7.1流域的雅各布港冰川（圖3中方框所示區(qū)域）是格陵蘭體積損失最嚴(yán)重的冰川，因此，對該冰川所在區(qū)域進(jìn)行進(jìn)一步分析。圖5所示為ICESat獲得的2003—2009年以及Cryosat-2獲得的2011—2017年雅各布港冰川地區(qū)高程變化速度，兩者之差（Cryosat-2結(jié)果減去ICESat結(jié)果）如圖6所示。由圖5可知，雅各布港冰川地區(qū)呈明顯的高程下降趨勢，且下降速度較快，冰川中心區(qū)域超過了-5 m·a－1，對比兩個時間段的結(jié)果可知，高程下降的區(qū)域和速度都有所增加。根據(jù)圖6結(jié)果，其中變化趨勢為正（即2011—2017年的消融速度相較于2003—2009年的速度有所減緩）的區(qū)域僅為3 700 km2，而變化趨勢為負(fù)（即2011—2017年間消融速度相較于2003—2009年的速度有進(jìn)一步增加）的區(qū)域面積則高達(dá)16 300 km2，由此可以明顯看出該冰川及周邊區(qū)域的體積損失在近10年中仍在進(jìn)一步加劇，其范圍也明顯擴大。

圖5 2003—2009年（左）與2011—2017年（右）雅各布港冰川高程變化（單位：m·a－1）

圖6 2003—2009年與2011—2017年雅各布港冰川高程變化差異（單位：m·a－1）

4 結(jié)語

基于Cryosat-2衛(wèi)星測高數(shù)據(jù)，以5 km×5 km分辨率確定了格陵蘭冰蓋在2011年1月至2017年12月間的高程和體積變化，對格陵蘭冰蓋流域進(jìn)行了詳細(xì)分析，并與ICESat激光測高衛(wèi)星獲得的2003—2009年間冰蓋高程和體積變化進(jìn)行了對比。結(jié)果表明，2011—2017年間，格陵蘭冰蓋仍處在明顯的體積損失狀態(tài)，其平均高程變化約為（-12.3±0.38）cm·a-1，對應(yīng)的體積變化則為（-220.38±0.60）km3·a-1。進(jìn)一步對各流域在不同時期的高程和體積變化速度的對比分析則表明，在格陵蘭內(nèi)部地區(qū)，冰川體積的積累速度已經(jīng)有明顯減緩，而在沿海地區(qū)，冰川流失現(xiàn)象則進(jìn)一步加劇。

致謝：感謝歐空局（ESA）提供Cryosat-2數(shù)據(jù)產(chǎn)品、布里斯托大學(xué)Bamber教授提供格陵蘭地區(qū)1 km分辨率地表類型格網(wǎng)數(shù)據(jù)及美國國家航空航天局Zwally教授提供格陵蘭冰蓋冰川流域系統(tǒng)劃分。