基于Bedmap2與冰雷達數據的南極局部冰蓋三維建模

劉 春 蘇小崗， 孫 波 李巍岳 陳 昀，

1 同濟大學測繪與地理信息學院，上海市四平路1239號，200092

2 中國極地研究中心，上海市浦東區金橋路451號，200136

在南極冰蓋物質平衡、冰芯研究、冰下水文環境探測以及冰蓋模式等研究中，冰蓋內外部地形數據是不可缺少的特征參數［1］。冰蓋三維模型能夠將冰表面、冰下基巖地形和冰蓋內部結構統一于一體，在南極冰蓋變化研究中發揮重要作用。隨著衛星遙感、冰雷達和深冰鉆探技術在南極的成功應用，南極冰蓋表面以及冰下地形探測取得很大進展［2］。在冰蓋表面地形探測方面，2009年Bamber等［3］利用ICESat數據和ERS－1測高數據構建了空間分辨率為1km 的冰蓋表面DEM，其覆蓋范圍涵蓋南極內陸大部分區域；2010年Tsutomu 等［4］利用InSAR與ICESat數據在南極毛德皇后地東部沿海冰架部分制作了分辨率為50 m 的冰架表層DEM，在精度與分辨率上都有較大提高。在冰下地形探測方面，英國南極局在2001年和2013年先后發布南極Bedmap1 和Bedmap2（Bedrock Mapping Project）數據庫，其中包括整個南極地區的冰蓋厚度信息、冰蓋表層DEM 和冰下基巖DEM，其空間分辨率均為1km［5－6］；2010年崔祥斌等［7］利用中國第21和第24次南極科學考察期間獲取的冰雷達數據，在Dome A 區域中心30km×30km 范圍內構建了網格分辨率為140.5m 的冰厚分布和冰下地形DEM，獲取了較為詳細的Dome A 冰下地形特征。可以看出，多種探測數據的綜合利用，已經可以獲取南極絕大部分區域的冰蓋表層和冰下基巖DEM；在小范圍區域可以獲得較高分辨率的冰蓋地形模型，但在冰蓋內部結構的獲取與冰蓋內外部地形的統一方面有所不足。本文綜合使用包括南極Bedmap2數據和中國第29次南極科學考察（2012－11）期間獲取的冰雷達數據，構建包括冰蓋表層DEM、冰蓋內部等時層和冰下基巖DEM 的南極局部冰蓋三維模型，并對模型建立過程中的數據處理流程作詳細介紹。

1 研究方法

1.1 數據源

冰雷達（ice radar）是基于電磁波理論，通過雷達回波研究冰雪介質特征的地球物理探測方法［8］，20世紀60年代被引入南極冰蓋的探測研究中［2］。本次研究選取中國第29次南極科學考察期間獲取的冰雷達數據。數據采集首次使用中國科學院電子學研究所自主研發的深冰雷達探測系統，采用車載冰雷達方式。雷達波在冰下距離向分辨率為1m，探測深度大于3 000m。

Bedmap2數據在Bedmap1 基礎上整合應用許多最新的探測數據，如激光衛星測高數據、冰雷達探測數據和最新的衛星遙感數據等，形成了覆蓋南極全部區域的冰蓋及冰下基巖地形數據庫［9］。Bedmap2包含冰蓋表面高程、冰厚和冰下基巖高程3類柵格數據，空間分辨率均為1km［5］。建模過程中主要使用Bedmap2中的冰表面柵格影像，其數據來源主要包括3部分：內陸大部分平坦地區使用Bamber等制作的分辨率為1km 的冰表面DEM［3］；在地形較為復雜的山脈區域，使用俄亥俄州州立大學DEM（OSU DEM）；沿海冰架部分則使用Griggs等利用ICESat數據和ERS－1數據制作的冰蓋DEM［9］。Bedmap2數據可直接在網上下載（http：／／www.antarctica.ac.uk／／bas＿research／data／access／bedmap／）。

1.2 技術路線

冰蓋三維建模主要分為冰蓋表面DEM 獲取和冰蓋內部信息提取兩部分。冰蓋表層信息主要通過對南極Bedmap2冰表面柵格影像處理獲得，冰蓋內部信息包含冰下基巖地形和冰蓋內部等時層，主要通過冰雷達數據處理獲得。為了使冰蓋內外部地形數據相互統一，在處理過程中需要將冰雷達數據與Bedmap2數據置于同一坐標系下，最后結合冰蓋表面、等時層、冰下基巖DEM 構建南極局部冰蓋三維模型。圖1為冰蓋三維建模技術路線圖。

圖1 技術路線Fig.1 Technology roadmap

1.3 研究區域

東南極冰蓋中山站至Dome A 斷面是國際橫穿南極計劃（international trans－antarctic scientific expedition，ITASE）的核心斷面之一，沿途經過Lambert冰川東側上游、Gamburtsev 冰下山脈和Dome A 等南極科學研究熱點區域，其中Gamburtsev山脈被認為是南極冰蓋的發源地之一［10］。1996年以來，中國南極科學考察隊沿中山站－Dome A 路線進行多次考察，獲得冰面地形、淺部雪層特征、雪積累率、冰芯樣品等科學數據［11］，并建立昆侖站（2009－01）和泰山站（2014－02）。本次研究選取中山站至Dome A 斷面的Gamburtsev山脈地區的一塊區域，首次對其進行三維模型建立。如圖2，研究區域位于昆侖站與泰山站之間，區域面積11.3km×11.5km，距昆侖站約150km，距中山站約1 050km。

2 數據處理與建模

2.1 冰蓋內部信息提取

冰蓋內部結構和冰下基巖地形主要通過處理冰雷達數據取得。由于本次探測首次采用了中國自主研制的冰雷達系統，雷達本身的采集存儲系統、現場測量方式與以往所使用的國外冰雷達系統有所差別，所以在數據讀取、成像等方面與以往的數據處理方式略有不同。數據處理主要流程包括預處理、常規處理和時深轉換3部分。

圖2 研究區域Fig.2 The studied area

2.1.1 數據預處理

冰雷達采集存儲的數據分為道頭文件與數據文件，道頭文件存儲GPS數據與設備參數，數據文件存儲雷達反射信號的電壓值。數據文件分多個單文件存儲，每個文件名記錄該段冰雷達數據采集的開始時間。以WGS84 橢球為基準，將GPS點的經緯度投影到南極極方位立體投影坐標系下，得到冰雷達測線分布圖（圖3）。為了清晰分辨冰雷達圖像，將測線分為9個測段。由于GPS系統與冰雷達系統的數據采集時間和采集頻率并不一致，需要根據冰雷達數據采集時間對GPS數據的起止時間作相應處理，得到與每個冰雷達數據文件對應的GPS數據。然后，采用三次樣條函數插值法對GPS數據進行插值，得到每道冰雷達數據對應的經緯度信息。

圖3 冰雷達測線圖Fig.3 Iceradar surveying route

2.1.2 常規處理

冰雷達數據文件為二進制格式，可以被MATLAB讀取，其橫向單位為雪地車行進的時間，每一列稱為一道，代表冰雷達發射并接收一次電磁波；縱向單位是從天線發射到接收電磁波的時間差。本次采用的冰雷達系統每隔0.016s采集一次數據，一般研究并不需要如此高的探測精度。在數據讀取時，對原始數據橫向按每128道數據抽取一道，縱向則每10列抽取1列，根據雪地車的行進速度計算，處理后的數據橫向分辨率為10m。在完成數據的初步讀取和處理后，需要合并整理處于同一測段的圖像。由于冰雷達系統工作與雪地車行進狀態不一致，雷達數據剖面影像中會產生冗余道，數據處理時要剔除產生冗余道的數據，同時刪除與其對應的GPS數據，以保證二者相互匹配。

在冰雷達探測過程中，由于系統自身、外界干擾及冰蓋內部起伏等原因，使得接收到的回波信號存在各種誤差或錯誤。主要處理方法包括增益控制、濾波去噪和偏移歸位［12］。圖像處理中主要使用了探地雷達處理軟件Reflex－Win 4.5。

首先，由于電磁波信號從冰下深淺不同層面先后到達地面的振幅相差很大，為了能均勻地記錄并顯示回波信號，需要對圖像進行增益控制［13］。使用手動增益控制方式，通過人為判斷手動調節增益大小，使不同層位的信號強度處于一個合適的值。圖4（a）為原始冰雷達Z－scope剖面影像，圖4（b）為經過增益控制后的剖面影像。

其次，為了減少噪聲、雜波等影響，需將增益后的圖像進行濾波處理，通常采用的濾波算法有帶通濾波、預測反褶積和背景去噪等。在冰雷達數據采集過程中，噪聲頻率一般處于一個較穩定的范圍，可以通過剔除特定頻率的回波達到去噪的目的。圖4（c）為經過帶通濾波去噪后的冰雷達剖面圖像，噪聲點得到有效控制。

圖4 冰雷達Z－scope剖面影像處理Fig.4 Ice radar Z－scope profile Image processing

最后，由于在冰下地形較為復雜的冰巖界面，特別是在雷達剖面中的巖性突變點會產生繞射波，使雷達記錄中的反射點偏移其原來的位置。為了糾正偏移點的位置，使用探地雷達中比較通用的繞射掃描疊加偏移法，使雷達反射波自動偏移歸位到空間真實位置［14］。

2.1.3 時深轉換與冰下地形獲取

電磁波在冰層內的傳播時間和傳播速度反映了冰層的厚度信息。冰蓋內部介質以單一冰體為主，電磁波在冰蓋內部的傳播速度按國際普遍采用的統計值1.68×108m／s［8］。電磁波在基巖界面處的反射功率達到極大值，利用這一特點可以在Reflex－Win軟件中自動跟蹤并提取冰下基巖界面［15］。

南極冰蓋冰層主要由雪轉變而成，每次降雪都會在冰蓋表面形成一層新的物質。從表層向冰蓋內部，隨著深度的增加，其物理特征則顯示為層狀有序的變化，通常認定冰蓋內部連續層是“等時”的［16］。將冰蓋內部等時層數據與深冰芯數據結合，通過數值模擬可以克服直接觀測在時間和空間上的限制［17］。本次研究利用處理后的冰蓋剖面影像提取平均埋深在950 m 處的冰蓋內部等時層。在提取過程中，為了保證不同測線的剖面圖中提取的等時層處于同一連續的等時層，需要根據不同測線之間的首尾連接關系，確定不同測線的等時層位置。圖5為其中一個測段內的等時層和冰巖界面提取線。

由GPS現場實測冰面高程與埋深作差，可得到冰下地形高程值。然后，通過ANUDEM 算法，插值得到分辨率為100m 的冰下基巖與等時層DEM。ANUDEM 算法能夠利用粗糙度懲罰函數對結果進行插值修正，相比于其他插值算法，獲得的DEM表面更加平滑，與實際更加相符［18］，如圖6。

圖5 剖面線提取Fig.5 The extraction of profile line

圖6 冰下基巖與等時層與DEMFig.6 The DEM of bedrock and isochronous layer

2.2 冰蓋表層DEM 獲取

本次研究使用Bedmap2 中的冰表層柵格數據提取冰蓋表層DEM。首先，利用ArcGIS10.2中的掩膜提取工具獲取Bedmap2在研究區的冰表面柵格影像。然后，柵格轉點得到1km×1km的高程分布點。為了保證冰蓋內外部地形結構的統一，同時呈現更加詳細的冰表面地形，同樣利用ANUDEM 插值算法，將1km 分辨率的高程點數據插值為100 m 分辨率，最終得到的冰蓋表層DEM 如圖7所示。

2.3 三維建模

冰蓋三維模型的建立使用三維可視化軟件Voxler 3［19］。通過對冰雷達數據的處理，獲取100m 分辨率的冰下基巖DEM 與等時層DEM。然后，從Bedmap2數據中提取100m 分辨率的冰蓋表面DEM。3層數據統一于南極極方位立體投影坐標系下。由于3層DEM 的分布范圍和分辨率都相同，可以從中提取得到如式（1）所示的三維坐標數據：

式中，（Xn，Yn）為平面坐標，Zn1、Zn2、Zn3為相對應的基巖、等時層與冰蓋表面高程值。將數據以txt格式導入Voxler軟件，最終的模型建立結果如圖8所示。

圖7 冰蓋表面DEMFig.7 The DEM of ice sheet surface

圖8 冰蓋三維模型Fig.8 The 3D model of ice sheet

3 模型檢驗與分析

冰蓋三維模型的檢驗主要是對3層DEM 的檢驗。冰蓋表層DEM 誤差主要來自Bedmap2數據。所選研究區域處于南極內陸，其表層DEM由Bamber等根據衛星測高數據制作而成，高程誤差在±23m 左右，中誤差小于1m［5］。冰下基巖與等時層DEM 由冰雷達數據所得，誤差來源主要包括冰雷達數據本身存在的誤差與數據處理過程中產生的誤差，這些誤差最終造成DEM 模型與實際地理位置的差異。目前，對于此類誤差的檢驗方法中比較常用的是交叉點分析驗證法［20］。對于不同測線提取的冰蓋剖面線，如果處在同一連續等時層或基巖界面，則不同測線的交叉點分析誤差應為零，即交叉點的埋深值應該相同。本次研究中共有9條不同的測線、18個交叉點，圖9為等時層與基巖數據的交叉點分析結果。經過統計，基巖數據的平均誤差為32.67m，大部分交叉點誤差小于50 m，相比于第24次中國南極考察隊獲得的冰下地形誤差要小。由于內部等時層地形起伏較小，誤差小很多，經統計得等時層交叉點的平均偏差為24.56 m，相對于km 級的冰蓋該誤差屬于可接受范圍。

圖9 基巖與等時層數據交叉點誤差分析Fig.9 Cross－point analysis results of the bedrock and isochronous layer

本次研究建立的南極局部冰蓋三維模型覆蓋范圍為11.3km×11.5km，網格分辨率為100 m，冰厚變化介于900～2 000m 之間。冰蓋表層海拔范圍為3 679～3 745m，整體較為平坦，海拔由南向北逐漸降低。冰下基巖地形起伏相對劇烈，海拔范圍為1 729～2 718m。可以看出，冰下存在一條由北向南延伸的槽谷，東部可能存在另一個山谷，這與中國南極考察隊在Dome A 探測的冰下地形較為相似［7］，均為典型的冰川作用地貌。相比于冰下基巖，所提取的冰蓋內部等時層地形起伏度較小，但與基巖在起伏變化上具有一致性。

4 結 語

利用冰雷達數據在冰蓋內部地形探測中的優勢，結合南極Bedmap2冰表面數據，通過對兩種數據的綜合處理，有效地將冰蓋內外部地形結構統一于一體，在中山站到Dome A 斷面上的一塊區域首次構建了覆蓋范圍11.3km×11.5km 的南極局部冰蓋三維模型。模型包含了100 m 分辨率的冰下基巖、冰蓋內部等時層和冰蓋表層DEM，詳細展示了冰蓋內外部地形結構特征。在往后的研究中，如何進一步提高冰蓋三維模型精度，尤其是對冰下地形DEM 的網格分辨率與精度的提高將是極地研究的重要內容。在中山站到Dome A 斷面，尤其是冰下地形數據較為稀少的東南極冰蓋和南極冰穹制高點——Dome A 區域，是未來進行冰雷達探測與冰蓋三維模型建立的重要區域［2］。將冰蓋三維模型與冰蓋熱力－動力耦合模式相結合，進行冰層與冰流動力學研究，推演冰蓋演化和深冰芯鉆孔位置選定等，也是未來極地研究的重要方向。

致謝：感謝中國極地研究中心提供中國第29次南極科學考察獲取的冰雷達數據，感謝課題組各位老師、同學的建議和幫助，感謝極地研究中心老師對本次工作的支持與指導。

［1］唐學遠，孫波，李院生，等.南極冰蓋研究最新進展［J］.地球科學進展，2009，24（11）：1 210－1 218（Tang Xueyuan，Sun Bo，Li Yuansheng，et al.Some Recent Progress of Antarctic Ice Sheet Research［J］.Advances in Earth Science，2009，24（11）：1 210－1 218）

［2］張勝凱，鄂棟臣，周春霞，等.南極數字高程模型研究進展［J］.極地研究，2006，18（4）：301－309（Zhang Shengkai，E Dongchen，Zhou Chunxia，et al.Progress on the Antarctic Digital Elevation Model［J］.Chinese Journal of Polar Research，2006，18（4）：301－309）

［3］Bamber J L，Gomez－Dans J L，Griggs J A.A New 1km 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－111

［4］Tsutomu Y，Koichiro D，Kazuo S.Combined Use of In－SAR and GLAS Data to Produce an Accurate DEM of the Antarctic Ice Sheet：Example from the Breivika－Asuka Station Area［J］.Journal of Polar Science，2010，4（1）：1－17

［5］Fretwell P，Pritchard H D，Vaughan D G，et al.Bedmap2：Improved Ice Bed，Surface and Thickness Datasets for Antarctica［J］.The Cryosphere，2013，7（1）：375－393

［6］Lythe M B，Vaughan D G.BEDMAP：A New Ice Thickness and Subglacial Topographic Model of Antarctica［J］.Journal of Geophysical Research：Solid Earth（1978－2012），2001，106（B6）：11 335－11 351

［7］崔祥斌，孫波，田鋼，等.東南極Dome A 冰雷達探測：冰厚分布和冰下地形［J］.科學通報，2010，55（3）：268－273（Cui X B，Sun B，Tian G，et al.Ice Radar Investigation at Dome A，East Antarctica：Ice Thickness and Subglacial Topography［J］.Chinese Sci Bull，2010，55（3）：268－273）

［8］Bianchi C，Forieri A，Frezzotti M，et al.Radio Echo Sounding（RES）Investigations at Talos Dome（East Antarctica）：Bedrock Topography and Ice Thickness［J］.Annals of Geophysics，2003，46（6）：1 265－1 270

［9］陳昀，孫波，劉春，等.南極冰蓋地形數據庫BEDMAP2述評［J］.極地研究，2014，26（2）：255－261（Chen Yun，Sun Bo，Liu Chun，et al.The Analysis of A New Antarctic Topography Database：Bedmap2［J］.Chinese Journal of Polar Research，2014，26（2）：255－261）

［10］崔祥斌，孫波，田鋼，等.東南極冰蓋中山站至Domeb A 斷面冰雷達探測初步結果：冰厚和冰下地形［J］.科學通報，2010，55（19）：1 937－1 943（Cui Xiangbin，Sun Bo，Tian Gang，et al.Preliminary Results of Ice Radar Investigation Along the Traverse Between Zhongshan and Dome A in East Antarctic Ice Sheet：Ice Thickness and Subglacial Topography［J］.Chinese Sci Bull，2010，55（19）：1 937－1 943）

［11］任賈文，秦大河，效存德.東南極冰蓋中山站－Dome A 斷面路線考察的初步結果［J］.冰川凍土，2001，23（1）：51－56（Ren Jiawen，Qin Dahe，Xiao Cunde.Preliminary Results of the Inland Expeditions along a Transect from the Zhongshan Station to Dome A，East Antarctica［J］.Journal of Glacilolgy and Geocryology，2001，23（1）：51－56）

［12］King E C.Ice Stream or not？Radio－Echo Sounding of Carlson Inlet，West Antarctica［J］.The Cryosphere，2011，5（4）：907－916

［13］王甜甜，孫波，關澤群，等.冰雷達探測數據處理方法研究［J］.極地研究，2013，25（2）：197－204（Wang Tiantian，Sun Bo，Guan Zequn，et al.Research on Radio－Echo Sounding Data Processing：a Case Study at Dome A，Antarctica［J］.Chinese Journal of Polar Research，2013，25（2）：197－204）

［14］徐玉增，盧海.偏移繞射技術在探地雷達資料處理中的應用［J］.能源技術與管理，2010（5）：17－18（Xu Yuzeng，Lu Hai.The Application of Offset Diffraction Technique in Ground Penetrating Radar Data Processing［J］.Energy Technology and Management，2010（5）：17－18）

［15］Welch B C，Jacobel R W.Analysis of Deep Penetrating Radar Surveys of West Antarctica，US－ITASE 2001［J］.Geophys Res Lett，2003，30：1 444

［16］Eisen O.Inference of Velocity Pattern from Isochronous Layers in Firn，Using an Inverse［J］.Journal of Glaciology，2008，54（187）：613－630

［17］Nereson N A，Raymond C F，Jacobel R W，et al.The Accumulation Pattern Across Siple Dome，West Antarctica，Inferred from Radar－Detected Internal Layers［J］.Journal of Glaciology，2000，46（152）：75－87

［18］張棟.基于ICESat和冰雷達數椐的南極Lambert冰川流域冰蓋特征提取研究［D］.南京：南京大學，2013（Zhang Dong.Research on Ice Sheet feature Extraction of Lambert Glacier Drainage Basin，Antarctic based on ICESat and Ice Radar data［D］.Nanjing：Nanjing University，2013）

［19］梁慶華，劉明偉，胡玉超.基于Voxler的井下瞬變電磁三維可視化研究［J］.礦業安全與環保，2013，40（5）：21－24（Liang Qinghua，Liu Mingwei，Hu Yuchao.Study on 3D Visualization of Mine Transient Electromagnetic Detection Based on Voxler［J］.Mining Safety ＆Environment Protection，2013，40（5）：21－24）

［20］Rippin D M，Bamber J L，Siegert M J，et al.Basal Topography and Ice Flow in the Bailey／Slessor Region of East Antarctica［J］.Journal of Geophysical Research：Earth Surface（2003－2012），2003，108（F1）：1－11