0 引言

经过50多年的发展,地理信息系统已经发展成为地理信息科学和地理信息应用的综合高技术领域[1]。由于GIS对空间数据有输入、输出、分析、查询、运算等一系列的功能,GIS已经成为空间数据挖掘的重要工具。GIS模拟空间数据变化的操作简便、处理空间数据速度快、可视化程度高,因此,使用 GIS建模来模拟事物的变化发展也成为各个学科领域的热点。

在冰川学的研究中,探索冰川未来的变化一直被广泛关注,且通常使用动力学模型来模拟冰川长期变化[2-4]。动力学模型所涉及的物理参数较多,模拟过程十分复杂,且处理速度较慢。本文基于GIS栅格操作的方法建模,在简化复杂动力学模拟的同时,抓住冰川变化的主要问题,提高了冰川变化研究获取重要几何参数的效率,对冰川长期变化的快速模拟具有参考价值。

1 研究背景

1.1 国内外研究现状

GIS在极地考察中普遍使用,比如通过 GIS数据处理建立高精度的冰川数字高程模型[5],通过 ArcGIS软件发布南极参考高程模型(REMA,网址: https://www.pgc.umn.edu/data/rema/)等。虽然GIS在极地应用已经较多,但是通过GIS方法开展冰川变化模拟还非常少见。目前普遍使用动力学模型模拟冰川变化,张良甫等人使用 Elmer/Ice动力学模型,对南极冰盖进行了研究[6]; 季青原等人使用 PISM 冰盖模式,对 Lambert冰川-Amery冰架系统的流动过程进行了模拟,并讨论了该模拟方法的适用性[7]; 周在明等人使用动力学原理,模拟了天山乌鲁木齐河源 1号冰川的运动[8]。通常情况下,动力学模型模拟冰川运动计算量大,运算周期长。能否通过GIS方法快速模拟冰川长期变化？本文正是利用己有实测数据的北极山地冰川作为研究对象,采用 GIS建模的方式预估冰川长期变化,探讨GIS方法模拟冰川变化的可行性。

1.2 研究区概况

北极斯瓦尔巴群岛地区(74°N—81°N,10°E—35°E)与格陵兰岛相邻,位于北冰洋之上,总面积为62 000 km2,如图1所示,其中约60%被冰川(帽)所覆盖。

中国黄河站位于斯瓦尔巴群岛西北部的新奥尔松地区,是目前纬度最靠北的站区(78°55′N,11°56′E)[9]。受北大西洋暖流影响,气候与同纬度地区相比较为温暖湿润,是北极同纬度气温最高的地区之一,同时也是北极增温最快的地区之一[10]。

Austre Lovénbreen冰川(简称 A冰川)位于新奥尔松附近,距黄河站 6.2 km,周围山峰的最高海拔为876 m[11]。A冰川属于多温型山谷冰川,其主体温度具有接近或处于融点的特点,是暖湿气候条件形成的冰川[12]。

2 研究数据和方法

2.1 研究数据

本文使用的主要研究数据包含A冰川冰面和冰床数字高程模型(Digital Elevation Model,DEM)[5],近期物质平衡实测分析结果[13],以及 Elmer/Ice模式对A冰川的模拟结果。

冰面观测数据是2009年GPS测点数据,该数据通过空间插值可以获取冰川冰面的 DEM。冰床 DEM 综合使用 GPS和探地雷达(Ground Penetrating Radar,GPR)数据获取,利用 GPS获取的冰面高程减去 GPR获取的冰层厚度即可获得冰床DEM,该DEM包含了A冰川整个山谷的地形。

冰川物质平衡(mass balance)是指一段时间内冰川的收入和支出量[14],冰川物质平衡线(equilibrium line)即冰川上的雪线,它是指冰川是积累量与消融量正好相等的点的连线[15]。物质平衡通常指表面物质平衡(SMB),主要包含两个参数: 零物质平衡线高度(ELA)、冰川物质平衡随海拔分布的梯度。根据已有研究成果,本文选取梯度参数为1/200,模拟初始值ELA为470 m。

Elmer/Ice冰流模式是一个建立在Elmer流体力学模式基础上的三维有限元处理软件,它使用三角网格对研究区进行区域划分。相对于使用正则化的有限差分求解的冰盖模式,使用三角网格划分的 Elmer/Ice冰流模式可以更好地刻画小范围的复杂地形[6]。目前Elmer/Ice已经成为冰川模拟的主流方法,其模拟结果在冰川模拟领域获得了广泛的认可。在该实验中,Elmer/Ice模式所使用的参数是如表1所示。

2.2 方法

本文的GIS模拟方法聚焦于提取冰川模拟结果中,冰川体积、面积、主流线长度和冰厚,并未考虑影响冰川形态变化的动力学因素。模拟过程中使用数学方法限定冰川的积累和消融以及其逐年变化。

2.2.1 冰川冰面变化模拟方法

根据 2.1 节介绍的表面物质平衡参数,为了更真实地模拟冰川随气候变化的情况,本文在模拟中增加了ELA随时间的变化率(rate)这一参数。实验中,模拟冰川表面变化时,首先将2009年冰川表面DEM转化为ASCII文件,然后使用JAVA程序实现对每一个像素值(冰面高程)的计算,最后将结果重新转换为DEM。

公式(1)中,H表示某像元的值(即某像素代表的冰面高程值); n为时间方向上迭代计算的步长; ELA即零物质平衡线高度; rate表示物质平衡线高度随时间的变化率; gra即为冰川物质平衡随海拔分布的梯度。

根据公式(1),本文进行了五组实验。结合近期物质平衡实测分析结果,五种实验都将ELA初始值定为470 m,gra定为1/200,在实验5中,当ELA随年际变化时,根据 IPCC发布的 RCP(Representative Concentration Pathway)8.5,rate取值为2 m·a–1。这五种实验参数如表2所示。

2.2.2 冰川特征提取

为了分析冰川模拟的结果,需要从冰面和冰床 DEM中提取冰川轮廓、厚度、体积、主流线长度等几何特征。

本文使用 Model Builder建模提取冰川的上述几何特征。模型构造器(Model Builder)是ArcGIS提供的构造地理分析、处理工作流和脚本的图形化数据建模工具[16]。Model Builder由输入数据、空间处理工具和输出时局3个基本结构构成[17],其提供的图形化建模工具能够继承三维、空间分析和地理统计等多种空间分析处理工具[18]。在这个模型中,分别定义不同的图形代表输入数据、输出数据、空间处理工具,他们以流程图的形式组合并且可以指向空间分析操作功能[19]。

Model Builder主要完成冰厚的提取和冰川外轮廓的提取,其过程如下:

1)使用冰床DEM,掩膜提取冰面DEM。

2)一一对应,比较冰面DEM和冰床DEM每个像素的高程值,提取较大的高程值,生成新的DEM,记为 surface。

3)surface和冰床 DEM的每个像素值对应相减,获取冰川厚度,只选取大于0.1 m的厚度,生成最终的冰川厚度栅格图。

4)提取最终冰川厚度栅格图的边界,并且使用边界圆滑工具,获得最终的冰川轮廓。

根据最终获得的冰川轮廓,使用 ArcGIS里面的测量工具可以测得冰川主流线的长度; 根据最终获取的冰川厚度,使用 ArcGIS中的识别工具,可以获取任意一点的冰川厚度,本文中,将该点记为 O点; 根据最终获得的冰川轮廓,掩膜提取冰面和冰床DEM,使用ArcGIS中的功能性表面工具,分析掩膜提取后的冰面和冰床 DEM,可以获取最终冰川的体积。

3 结果及讨论

3.1 冰川特征模拟

本文提取了 60—100年的模拟结果,对冰川主流线长度、某点冰厚、冰川表面积和冰川体积做了初步对比。

图2(a)中,冰川主流线长度随年份的推移逐年递减。实验4的结果在每个年份都是最大的,实验4模拟中,冰川的长度从约3 300 m递减到约3 100 m,其次是实验 3,然后是实验 2,最小的是实验 1的情况,但此三种情况冰川的长度比较相近,都是从约3 100 m递减到2 800 m,这是由于迭代的次数越多,GIS模拟中冰川的消融量就越大,因此,在时间上的迭代步长越小(n的取值越小),模拟结果的值会越小。在实验5中,由于ELA逐年增加 2 m,冰川的积累区逐年减小,消融区逐年增加,这导致冰川的退缩速度会大大快于ELA不变情况下的退缩速度,因此冰川长度从60年处就已经仅剩约3000 m,到了100年处则不到2 000 m,这就是情景二模拟中,50年退缩1 000 m的情况。

图2(b)中,冰川上O点的厚度也是处于逐年递减的情况,这说明O点处于冰川的消融区,且其海拔在2009年冰面DEM上显示的是低于470 m。实验4的结果仍然高于其他实验1、2、3的模拟情况,其厚度从大于110 m到小于110 m,实验1、2、3的模拟也遵循 n的取值越小,模拟结果越小的规律,但其厚度的变化则是在100 m附近由厚变薄。实验5,由于O点处于消融区,其消融速度由于 ELA的升高快于实验 1—4的速度,因此仅在50年中,O的厚度则减少了70 m。

图2(c)是冰川表面面积,其变化规律与O点厚度的变化规律一致。实验4的结果远高于实验1—3的模拟情况; 实验5中,ELA变化时面积的缩减速度也非常快。

图2(d)中,体积的变化规律也与O点厚度的变化规律一致,但是在研究变化趋势时,可以发现在五个实验中,体积变化的趋势处于一种平滑下降的趋势。据此,本文推测: GIS模拟冰川体积的变化时,可能符合一定的函数规律。实验1—4模拟情况的变化趋势一致,本文选取实验2和实验3这两种情况深入研究; 对于实验 5,由于该实验的 ELA是变化的,因此也特别对此模拟结果进行了深入的研究。

3.2 冰川体积

通过3.1节可知,冰川模拟特征提取中,冰川体积的变化规律性较强,模拟效果较好。因此,本文对于实验2、3、5,模拟了100年的冰川变化结果,并且每隔10年输出模拟的冰川体积。

由图3可知,实验2和实验3模拟的值相差不大,且两者的变化趋势较为平缓; 实验 5中,冰川体积随时间变化的趋势同样比较平缓。为定量研究其变化趋势,本文使用最小二乘法拟合的方式,对模拟结果进行了简单的拟合。

最小二乘法曲线拟合的基本思想是: 使所有数据点与拟合点的误差的平方和达到最小[20]。函数的拟合只是一种近似地求取数值之间关系的分析方法,因此在拟合函数的过程中,根据散点的分布和变化趋势,本文使用了比较简单直观的高次常函数拟合方法。为使得拟合结果不至于过于复杂,在拟合的过程中,尽量选取“确定系数”(R2)接近1,当拟合的函数最高次数是k+1但R2不再变化时,则认为拟合结果的R2收敛于k次; 另外加入自由度调整R2(Adjusted R2)作为参考。R2是通过数据的变化来表征一个拟合好坏的判定值,R2的取值范围是[0,1],越接近 1,表明方程的变量对 y的解释能力越强,这个模型对数据拟合的结果也较好; Adjusted R2在接近 1时,也表示匹配结果较好。因此,本文对比了在拟合过程中,R2和Adjusted R2的变化过程。

根据表3—5可知,对于实验2和实验3这两种情况,在二次函数拟合时R2和Adjusted R2都已经达到 1,因此两者可以高度拟合于二次函数; 对于实验5,在四次拟合时R2和Adjusted R2达到了1(五次函数拟合时由于数据较少,出现了错误警告),因此实验 5的模拟中,体积变化情况可以高度拟合于四次函数。

根据拟合情况,本文求出了拟合之后的函数(如表6所示)。

3.3 与Elmer/Ice模拟结果比较

本文在使用Elmer/Ice模拟冰川变化时,模拟的初始ELA为470 m,且分ELA不随年份推移变化和ELA每年增加2 m这两种情景,并提取了模拟100年的结果。本文分析了这100年中A冰川在五种实验中轮廓变化和体积变化结果(在计算中,发现长期模拟的冰川厚度、表面积、主流线长度与 Elmer/Ice的结果相差较大,大部分都在5%以上,因此下文中不讨论这几项)。

图4是实验 1—4下的模拟,此时 Elmer/Ice模拟的ELA是不变的。如图可知,从形态变化来说,Elmer/Ice模拟与GIS模拟相比,冰川的形态大体相同,特别是与实验 4模拟情况相差最小。但事实上,Elmer/Ice模拟的冰川末端总是长于GIS模拟的结果,这是因为在模拟中,Elmer/Ice考虑到了冰川的流变,高海拔地区的冰川会流动到低海拔处进行补充,因此Elmer/Ice模拟的冰川轮廓在相同年份总是大于GIS模拟结果。

对于实验5的模拟和Elmer/Ice模拟(ELA变化)结果,如图5,60年和70年,冰川形态在两种方法的模拟下基本一致; 80—90年,GIS模拟冰川西侧支流已经断流,但仍然保留其小部分分支,Elmer/Ice模拟下,西侧支流则已经完全消失; 100年时,GIS模拟的冰川东侧支流断流,主流线和上游也已经断开,Elmer/Ice模拟的冰川形态中,主流线和东、西支流以及上游都已经断开。但与图4相同,Elmer/Ice模拟的冰川末端一直长于 GIS模拟结果,这同样是Elmer/Ice模拟所考虑的冰川流变的情况导致。

对于冰川体积,在不同的情景下的 100年模拟中,GIS模拟与Elmer/Ice模拟结果相差率差别较大(如图6)。实验2和实验3中,体积相差率变化趋势一致,且相差率在10—30年、70—80年时最接近于0%(左侧主坐标轴),实验 5中,相差率由负到正,逐渐增加,且10—40年的相差率都接近于 0%(右侧次坐标轴),如果将相差率规定在3%内认为GIS模拟结果有效,那么在实验2和实验3中,前90年的模拟结果都是有效的。实验5中,前60年的模拟都是有效的。

GIS模拟冰川变化时,没有考虑冰川流变,只考虑了物质积累和消融; Elmer/Ice则考虑了冰川的物理性质,其在物质积累和消融的情况下,还加入了高海拔物质流向低海拔的情况。因此,两种模拟方法的物质平衡情况是比较相近的,这使得GIS模拟中,冰川体积在十年到几十年的尺度上,和Elmer/Ice的结果是一致的。

4 结论

本文通过GIS方法模拟了冰川长期变化,提取了冰川主流线、冰川表面积和冰川体积,并且着眼于冰川体积做了比较详细的研究,还与南北极冰川模拟的主流模式 Elmer/Ice的结果相比较,探讨了GIS模拟的准确性和可行性,得到了以下结论。

1.GIS方法可以分不同情景,比较真实地模拟冰川长期变化。提取冰川主流线、冰川表面积、冰川体积的结果进行分析,并研究冰川随时间的发育情况,能快速获取可视化的整条冰川变化特征,对冰川冰厚变化的研究同样有效。

2.GIS方法模拟冰川长期变化时,实验2、3中冰川体积变化符合二次函数变化规律,实验 5中冰川体积变化符合四次函数变化规律。这有利于定量分析冰川体积的变化。

3.与Elmer/Ice模拟冰川轮廓结果相比,GIS模拟的实验 4,对冰川轮廓的模拟结果和 Elmer/Ice模拟最相近; 实验 5中,冰川轮廓在 60年和70年与 Elmer/Ice模拟结果相近。误差的产生是由于 Elmer/Ice模式考虑了冰川流变情况下高海拔对低海拔地区的物质补充。

4.关于冰川体积的预测,如果以 Elmer/Ice模式为参照并将同期体积误差控制在3%以内,实验2和实验3最长可以模拟未来90年的冰川变化,实验5最长可以模拟60年。因此,在十年至几十年的尺度上,GIS模拟结果与Elmer/Ice结果比较一致。

参考文献

1 周成虎.全空间地理信息系统展望[J].地理科学进展,2015,34(2): 129–131.

2 ZWINGER T,MOORE J C.Diagnostic and prognostic simulations with a full Stokes model accounting for superimposed ice of Midtre Lovénbreen,Svalbard[J].Cryosphere,2009,3(2): 217—229.

3 RÉVEILLET M,RABATEL A,GILLET-CHAULET F,et al.Simulations of changes to Glaciar Zongo,Bolivia (16°S),over the 21st century using a 3-D full-Stokes model and CMIP5 climate projections[J].Annals of Glaciology,2015,56(70): 89—97.

4 ZHAO L,TIAN L,ZWINGER T,et al.Numerical simulations of Gurenhekou Glacier on the Tibetan Plateau using a full-Stokes ice dynamical model[J].Journal of Glaciology,2014,60(219): 71—82.

5 艾松涛,王泽民,鄂栋臣,等.基于GPS的北极冰川表面地形测量与制图[J].极地研究,2012,24(1): 53—59.

6 张良甫,唐学远,杨树瑚,等.基于Elmer/Ice冰盖模式的南极Gamburtsev山脉Lambert冰流区域的数值模拟研究[J].极地研究,2017,29(3): 391—398.

7 季青原,王帮兵,孙波.PISM冰盖模式对Amery冰架流速场模拟的适用性[J].极地研究,2015,27(3): 229—236.

8 周在明,李忠勤,李慧林,等.天山乌鲁木齐河源区1号冰川运动速度特征及其动力学模拟[J].冰川冻土,2009,31(1): 56—61.

9 孙维君,闫明,艾松涛,等.北极新奥尔松地区Austre Lovénbreen冰川温度变化特征[J].武大学报,2016,41(1): 79—80.

10 程振波,石学法,陈志华,等.2007中国北极黄河站夏季科学考察及黄河站概况[J].海洋科学进展,2008,26(1): 113—114.

11 极地测绘科学国家测绘局重点实验室.南北极地图集[M].北京:中国地图出版社,2009: 159.

12 秦大河.冰冻圈科学辞典[M].第二版.北京:气象出版社,2016: 231.

13 徐明星,闫明,任贾文,等.北极Svalbard地区Austre Lovénbreen和Pedersenbreen冰川表面物质平衡和运动分析特征[J].极地研究,2010,22(1): 10—22.

14 雒园.新疆天山参照冰川物质平衡与气象要素的关系及其模拟[D].兰州: 西北师范大学,2012.

15 鞠远江,刘耕年,张晓咏,等.山地冰川物质平衡线与气候[J].地理科学进展,2004,23(3): 43—49.

16 朱一晗,薛丰昌,唐步兴,等.小麦自然水分亏缺干旱评估模型的GIS建模[J].测绘通报,2017,(11): 101—107.

17 李元杰.基于Model Builder的林业专题图境界线自动提取方法[J].林业调查规则,2014,39(2): 11—13.

18 周春峰.浅谈基于ArcGIS的影像批量裁剪[J].北京测绘,2015 (6): 136—137,146.

19 汤国安,杨昕.ArcGIS地理信息系统空间分析实验教程 [M].第二版.北京:科学出版社,2012: 579.

20 贾小勇,徐传胜,白欣.最小二乘法的创立及其思想方法[J].西北大学学报 (自然科学版),2006,36(3): 507—511.