基于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