基于GPR探测的长江源地区冰川与冻土厚度研究

周黎明,张 杨

(长江科学院 水利部岩土力学与工程重点实验室,武汉 430010)

0 引 言

长江源地区地处青藏高原腹地,素有“中华水塔”的美誉,是中国重要的生态安全屏障,同时也是气候变化的敏感响应区和生态环境脆弱区。长江源地区平均海拔4 760 m,流域面积约13.82万km2。长江源地区水系包括北源楚玛尔河水系、正源沱沱河水系、南源当曲水系以及干流通天河水系。近50 a来,在全球气候变暖和人类活动的双重影响下,长江源地区雪线上升、冰川退缩、水土流失、荒漠化和草地退化等问题日益凸显,已对长江源区域乃至整个长江流域的生态安全和民生经济安全造成威胁。冰川作为响应高原气候变化的指示器,其厚度的精确探测和淡水资源储量的准确计算对于研究长江源地区气候变化和推测长江径流量具有重要意义。获取精确的冰下地形对于推断冰川的发育和运动学过程具有重要作用。另外,探测长江源地区的永久冻土位置对于研究长江源地区水土保持状况和淡水资源储量具有重要价值,同时也为青藏高原地区公路、铁路和水电工程的安全建设提供参考资料。

近年来,随着探地雷达(Ground Penetrating Radar, GPR)设备精度和便携性能的不断提高,我国利用探地雷达在多个冰川展开了探测工作,其中多数为地形崎岖的山地冰川,并取得了良好的效果。2008年,王璞玉等[1]对天山托木尔峰地区的青冰滩72号冰川进行了GPR多剖面的厚度测量,并使用差分全球定位系统(Global Positioning System,GPS)对探测点位进行精准定位。2009年,吴利华等[2]对天山博格达峰地区的四工河4号冰川进行了GPR测厚工作,并在地理信息系统(Geographic Information System,GIS)技术的支持下,采用Co-Kriging插值算法结合理想塑性体理论对冰川非测厚区域的厚度进行了重建。2014年,王玉哲等[3]使用GPR对祁连山老虎沟12号冰川的中流线及多条横剖面进行了厚度测量,获取了准确的冰下地形信息。2018年,李亚楠等[4]结合雷达冰川测厚数据和遥感信息对东昆仑山煤矿冰川冰储量进行了估算并获得冰床地形图。2018年,靳胜强等[5]利用GPR对西藏阿里地区喀喇昆仑山脉南部的嘎尼冰川进行了冰川厚度测量,并结合差分GPS数据和遥感影像数据绘制了冰下地形及冰厚分布图。

土壤层、冻土活动层和多年冻土层的介电性质存在明显差异,为GPR研究冻融界面及多年冻土层上限提供有效地球物理前提。20世纪70年代,GPR技术在高寒地区的浅层地质勘察中得以应用,随后在国内外逐渐推广[6]。杜二计等[7]在祁连山老虎沟和观山河源头等地开展了GPR冻土调查,获得了多年冻土上限的准确位置,并分析了冻土上限分布与海拔和气候变化的关系。武小鹏等[8]将GPR应用在青藏高原共玉高速公路多年冻土工程地质勘察中,获得了准确的多年冻土上限位置以及沿公路设计线路的分布情况,为高原公路施工提供了参考依据。2018年,单波等[9]利用GPR技术为青海玉树、新疆天山等高寒地区的输电线路工程地质勘察提供准确的冻土物探解译结果。

在国外,探地雷达已广泛应用于冰川和冰冻物质的研究,包括冰川和冻土的含水量、冰川地下的排水路线、冰川和冰冻物质的内部结构、冻土和冰川冰的破裂和变形、冰川中沉积物等[10]。Navarro等[11]基于探地雷达冰厚探测数据,对斯瓦尔巴群岛韦德尔-雅尔斯贝格冰川进行了冰储量估计,研究结果表明对冰川、冷冰和温带冰层采用不同的速度进行雷达成果时间-深度转换与采用不同的波速估计冰储量的误差范围是一致的。

冰川厚度探测方法还有合成孔径雷达技术[12],冻土厚度探测技术有高密度电法[13]、航空电磁技术[14]、地震折射波法、地震映像法、瑞雷面波法[15]、瞬变电磁法[16]等。合成孔径雷达技术能较好适用于冰川厚度探测,但不适用于冻土厚度探测。现有冻土厚度探测方法存在如下弊端:航空电磁技术成本较高;不能单独利用地震折射波数据检测山地多年冻土;地震映像法抗干扰能力弱,勘探深度有限,对冻土上限确定的准确度略低;瑞雷面波法地表介质不均匀或激发条件差时,勘探深度减小;瞬变电磁法能较好适用于冻土厚度探测,但不适用于冰川厚度探测。综合其他无损检测技术,GPR能用于冰川及冻土厚度探测,具有高效、高精度、能非常直观地反映出地质体的几何位置及形态的优点。

目前,针对长江源地区冰川和冻土的探测和研究仅限于遥感影像和地质情况调查,对于长江源头冰川厚度和长江源湿地多年冻土上限的GPR探测尚缺少研究。2022年和2023年7月,长江水利委员会长江科学院开展了自2012年以来的第12次、13次江源科考,并以江源“冰和碳”为这两次科考的主要内容,采用GPR技术进行冰川厚度探测以及冻土上限调查。其中,冰川探测位置选取长江正源沱沱河的发源地格拉丹东主峰冰川,冻土探测位置选择在查旦湿地。本文设计了多种典型冰川和冻土地质体模型并进行GPR波场数值模拟,研究了冰川和冻土探地雷达信号响应特征。分析了在长江源地区获取的实测冰川和冻土雷达数据,结合地质模型的GPR波场模拟结果,增加了实测数据的可靠性。对比分析了2022年和2023年冰川及冻土厚度的变化,为今后开展更大规模的长江源地区冰川冻土科考奠定了基础。

1 研究地区及其物性条件

1.1 格拉丹东主峰

在长江科学院历年的江源科考成果中,已有包括对冰川地区地理特征、地质特征、冰面高程变化监测等多方面的考察研究基础。2022年和2023年科考,在以往的基础上,登上冰川表面对冰川厚度进行精准探测,以掌握长江源地区冰川的分布状况。探测位置选择在长江正源沱沱河的发源地格拉丹东雪山,位于青海省西南部青藏边境唐古拉山中段,由一大片南北长达50余千米、东西宽约30余千米、共50余条山岳冰川群所组成。主峰格拉丹东峰位于91°E、33.5°N,海拔6 621 m,为唐古拉山脉的主峰(图1)。格拉丹东雪山冰川的伸缩变化明显地揭示了当地乃至青藏高原和全球气候的变迁,极具科研价值。

图1 格拉丹东主峰冰川地理位置及冰川卫星照片

1.2 查旦湿地

2022年和2023年江源科考,还进行了长江源地区冻土调查研究工作。调查地点选择在地下水含量丰富的查旦湿地,也是长江南源当曲的发源地。为了验证冻土探测的准确性,2022年和2023年科考冻土探测测线选择在查旦乡一个常年设置的地下水监测井附近,如图2所示。

图2 冻土探测测线位置及地下水监测井

1.3 长江源地区常见介质物性条件

长江源地区常见介质的电性特征见表1[17]。

表1 长江源地区常见介质的电性特征[17]

根据电磁波在介质中的波速和旅行时间可以计算界面深度为

h=vt/2 。

(1)

式中:h表示反射界面深度;t表示传播时间;v表示电磁波在介质中的传播速度,即

v=c/(εr)1/2。

(2)

式中:c为电磁波在真空中的传播速度,c=0.3 m/ns;εr为介质的相对介电常数。根据所探测区域介质的不同选择适当的电性特征参数(表1)可以获取目标体准确的位置信息。

电磁波在传播过程中,遇到不同的阻抗界面(介质分界面)时将产生反射波和透射波,其反射与透射遵循反射与透射定律,反射波能量大小取决于反射系数。反射系数r的数学表达式为

(3)

式中ε1、ε2为反射界面两侧的相对介电常数。

2 冰川和冻土GPR波场模拟

2.1 模型建立

本文采用GprMax[18]软件,基于时域有限差分法,设置满足稳定性条件的适当参数并结合PML[19]衰减边界条件,模拟GPR电磁波在各向同性无限半空间地质模型中的传播过程。本文设计地质模型尺寸均为10 m×20 m,网格间距为0.01 m,GPR在地表采用剖面法观测,天线间距1.0 m、探测点距0.2 m,相关介质参数见表1。图3中0～2 m深度的界面为直达波。

图3 冰川模型及GPR正演记录

图3(a)为2层水平层状冰川模型,其中沿深度方向0.0～9.0 m为冻结冰层、9.0～20.0 m为冰下基岩,且冰下基岩地形呈水平状。图3(b)为2层向右倾斜状冰川模型,其中冰岩界面的深度从左到右,从9.0 m变化到11.0 m。图3(c)为3层水平层状冰川模型,其中沿深度方向0.0～5.0 m为含水的融化冰层、5.0～9.0 m为冻结冰层、9.0～20 m为冰下基岩。冰下基岩地形呈水平状;冻结冰层的相对介电常数ε3=3;融化冰层的相对介电常数ε4=80;基岩的相对介电常数ε5=8。

2.2 GPR波场模拟结果

图3(d)为图3(a)所示模型的GPR正演记录,其中在深度方向9.0 m处出现明显反射信号水平状同相轴。由表1和式(3)可知,由于冻结的冰层与冰下基岩相对介电常数差异较大,冰、岩分界面为电磁波强反射界面,在探地雷达正演记录中,反射波同相轴出现的位置及形状与图3(a)所示模型介质分界面一致。

图3(e)为图3(b)所示模型的GPR正演记录,其中在深度方向7.0 m处开始出现明显反射信号向右倾斜状同相轴,反射波同相轴出现的位置及形状与图3(c)中所示模型介质分界面基本一致。图3(c)和图3(e)所示倾斜状冰下地形为冰川槽谷常见地形,多呈“V”形或“U”形。

图3(f)为图3(c)所示模型的GPR正演记录,其中在深度方向5.0 m和15.0 m处各出现强反射信号水平状同相轴,冰层(融化)和冰层(冻结)界面能较好对应理论模型深度,而冰层(冻结)和基岩界面与理论模型深度偏差较大,其原因为冰层(融化)的电磁波传播速度较小,而进入冰层(冻结)波速较大,导致计算结果深度偏大。

图4(a)为2层水平层状冻土模型,其中沿深度方向0.0～2.0 m为含水土壤层(即季节性融土层)、2.0～20 m为冻土层(多年)。

图4 冻土模型及GPR正演记录

图4(b)为图4(a)所示模型的探地雷达正演模拟记录,其中在深度方向2.0 m处出现明显的水平状反射同相轴。由表1和式(3)可知,由于融土层与冻土层之间存在介电性质差异,因此在探地雷达正演记录中,反射波同相轴出现的位置及形状与图4(a)中所示模型介质分界面一致。在冻土模型中,融土层与冻土层之间的电磁波反射界面振幅强弱还与两层介质中含水率有关。

3 长江源地区冰川和冻土探测

3.1 探测仪器设备及基本数据处理流程

冰川、冻土科考所使用设备为加拿大sensor &software公司制造的EKKO PRO型探地雷达仪,天线主频为100 MHz。GPR数据处理[20]要经过零点校正、背景噪声消除、直流滤波、Laplace滤波、归一化、调节增益、时深转换等处理步骤,得到最终突出地下介质特征的波形信息。其中,零点校正主要消除天线离地和天线延时的影响,背景噪声消除主要去除背景中含有的电磁噪声信号;直流滤波主要切除电磁波中含有的直流成分;Laplace带通滤波主要是消除地表环境产生的噪声干扰;归一化是为了在信号强度大的地方(通常为近地表地区)采用小增益,在信号衰减的地方(通常为深部)采用高增益;调节增益主要补偿电磁波在地下介质中传播时的衰减能量;时深转换是将时间剖面转换为深度剖面。GPR的地质解释是在数据处理后得到的雷达剖面上,根据反射波组的波场特征,通过同相轴追踪,确定反射波组的地质意义。

3.2 格拉丹东主峰冰川探测

格拉丹东主峰冰川探测所设探地雷达参数如表2所示。综合考虑了采集的原始数据的质量和人员设备安全等因素,设计了与格拉丹东主峰冰川中流线夹角呈61.52°的横剖面测线一条(图5),测线参数如表3所示。

表2 冰川探测探地雷达参数设置

表3 冰川探测测线参数

图5 冰川探测测线起点位置

格拉丹东主峰冰川GPR探测现场工作由多人完成,分别负责仪器搬运、点位标定、距离测量、雷达操作、现场记录和安全保障(图6)。

图6 冰川探测工作现场

2022年GPR探测数据处理与解译如图7所示。由图7可看出:从冰面深度d=0.0 m至冰下深度d=3.0 m处雷达信号表现为振幅较强,频带较低,同相轴连续性好,在d=3.0 m处存在振幅突变界面,结合第2节中图3(f)所示冰川模型GPR波场模拟记录,推断d=0.0～3.0 m的范围内为融化冰层;沿地表测线方向从x=0.0～9.5 m对应在深度d=11.0～12.5 m处存在振幅突变界面,振幅明显增强,在该界面上方无明显雷达反射信号,符合冰介质对电磁波信号弱衰减作用的特征,该界面下方存在复杂反射信号及多组线性同相轴(多次波),符合冰下基岩的GPR信号响应特征,推断从融化冰层底界面至该强振幅界面中间为冻结冰层,该强振幅分界面为冰、岩分界面。

图7 2022年冰川探测探地雷达数据解译

从图7解译结果中读取的冰厚信息由测线起点(远离冰川中流线一侧)的11.0 m逐渐变化为测线终点(靠近冰川中流线一侧)的12.5 m。冰下基岩坡度由式(4)计算得到。

(4)

2022年探测所得冰下基岩坡度为9.08°,地形走势呈向右倾斜状,结合图3(c)所示冰川模型和图3(d)所示GPR波场模拟记录,此次探测数据中解译的冰下地形符合冰川槽谷地形特征。冰川槽谷是冰川长期作用地表的结果,其形态受冰川蠕动、作用时间和基岩岩性等因素的影响。不同的槽谷形态对冰川运动速度的影响很大。通常需根据冰川槽谷形态(如“U”形、“V”形) 对侧向应力进行参数化,从而建立预测冰川运动过程的数值模型。由于2022年科考受地形因素限制,考虑到冰上作业的安全问题,横剖面测线未穿越冰川中流线,尚不能完全掌握格拉丹东主峰冰川冰下地形的整体情况。

2023年GPR探测数据处理与解译如图8所示,其中:对应在深度d=6.0～7.0 m处存在振幅突变界面,振幅明显增强,该界面下方存在复杂反射信号及多组线性同相轴(多次波),符合冰下基岩的探地雷达信号响应特征,推断该强振幅分界面为冰、岩分界面。冰川下部基岩不完整,导致基岩裂缝中有冰层,深度d=7.0～8.0 m范围内可能为冰层,深度d=8.0 m以下才为完整基岩层,导致在深度d=8.0 m的位置也存在振幅突变界面。冰川探测的厚度有待进一步采用实测资料加以佐证。

图8 2023年冰川探测探地雷达数据解译

3.3 查旦湿地冻土探测

查旦湿地冻土探测所设GPR参数如表4所示。 探测测线长度9.0 m。 2022年探测数据处理与解译如图9所示。 为突出界面分层效果采用波形图显示, 其中从土壤表面d=0.0 m至地下d=4.7 m处雷达信号表现为振幅较强, 频带较低, 同相轴连续性较好, 局部波形细密, 在d=4.7 m处存在振幅突变界面, 在d>4.7 m处电磁波频率升高, 振幅变小, 且局部波形细密。 结合第2节中图4(a)所示冻土模型和图4(b)所示GPR波场模拟记录, 推断在d=0.0～4.7 m的范围内为土壤层和季节性融土层; 在d=4.7 m为冻土上限位置。 冻土探测现场地下水监测井中可以清晰看到冻土层表面结冰, 经深度测量得到结冰位置位于d=4.73 m(此处未考虑气温对地下水监测井中的冻土层的影响)。 由此验证了科考采用探地雷达对长江源地区冻土调查的有效性。

表4 冻土探测探地雷达参数设置

图9 2022年冻土GPR数据解译

2023年探测结果如图10所示, 从土壤表面d=0.0 m至地下d=3.6 m处雷达信号表现为振幅较强, 频带较低, 同相轴连续性较好;在d=3.6 m处存在振幅突变界面, 在d>3.6 m处电磁波振幅突然减弱, 几乎消失。 推断在d=0.0～3.6 m的范围内为土壤层和季节性融土层; 在d=3.6 m为冻土厚度上限位置。

图10 2023年冻土探测探地雷达数据解译

4 结论和展望

本文基于时域有限差分方法和完全匹配层的吸收边界条件,通过设置多种冰川和冻土理论模型,对冰川和冻土进行探地雷达正演模拟及江源科考探测分析,采用GprMax软件进行数值模拟,得到以下结论:

(1)不同冰川、冻土模型所对应的GPR信号特征同相轴形态与介质分布形态一致,且反射信号振幅由不同介质之间的介电性质差异决定,其中融冰层与冻结冰层之间为较弱电磁波反射界面,冻结冰层与冰下基岩为强电磁波反射界面,融土层和冻土层之间的电磁波反射界面强度还与两层之间的含水率有关。

(2)探测结果表明,格拉丹东雪山主峰冰川厚度和查旦湿地冻土厚度上限均有不同程度降低,冰川厚度和冻土厚度上限观测是一个常年累月的结果,后续还需要持续进行观测,积累更多数据,分析变化趋势,以估算探测区域内冰储量,研究气候变化对冰川的影响。

本文依托2022年和2023年长江源科考对格拉丹东主峰冰川探测和查旦湿地冻土探测,探测解译结果均与理论模型研究结果符合,获得了冰川测线处准确冰厚信息、冰下地形信息和冻土上限的准确位置,验证了GPR方法在长江源地区进行冰川和冻土精准探测方面的可行性。在未来的江源科考工作中,将进一步提高GPR在长江源冰川探测和冻土调查中的应用范围,获取精确的三维冰川和冻土资料,估算冰储量和地下淡水资源储量,并引进航空探地雷达对受地形因素制约的冰川和冻土地区进行探测试验。