长白山天池火山2015—2022年InSAR变形与活动状态分析

2024-01-20熊国华季灵运刘传金

地震地质 2023年6期

熊国华季灵运* 刘传金

1)中国地震局第二监测中心,西安 710054 2)国家遥感中心地质灾害研究部,北京 100036

0 引言

长白山天池火山位于吉林省东南部的中朝边境线上,是一个休眠了300多年、仍具有喷发危险的复合式活火山(刘若新等,1992)。目前从外形上看,长白山天池火山的锥体形貌完整,破火山口完好,从地貌上可以清晰地辨别火山造锥的4个阶段。迄今为止,长白山天池火山可查证的爆炸式喷发活动至少有3次,规模最大的一次喷发发生于公元946年,曾造成巨大损失。长白山天池火山被认为是中国最危险且最具有潜在喷发危险的火山之一(Xuetal.,2013),其周围分布多条断裂,断裂附近、长白山天池火山与望天鹅火山周围地震活动频繁。据统计,2002—2005年期间,长白山天池火山发生异常扰动,其间地表出现隆升现象,地震活动亦变得频繁(Xuetal.,2012)。其后的2010—2020年间,地震发生频率有所降低,但仍处于活动状态,火山活动会引发火山灰、地震及海啸等多种灾害,且天池火山的蓄水量高达20.4亿立方米,一旦喷发,水体可能会与岩浆发生作用形成岩浆-气水喷发,并容易引发火山泥石流,造成更大危害。

跟踪天池火山的变化、明确火山活动情况及规律对认识火山喷发的危险性具有非常重要的意义。大地测量作为一种重要的技术手段,可通过地表形变信息反映地壳内部岩浆的状态,反演获取岩浆房的位置和深度,进而使我们了解火山活动变化情况(牛玉芬,2020)。传统的地表形变监测手段通常是基于单点进行的局部监测(刘俊清等,2012; 顾国华等,2019; 郭明瑞等,2022),具有精度高的优势,但无法实现广域监测。遥感技术作为一种非接触式的监测手段,可实现大范围监测,且合成孔径雷达干涉测量技术(Interferometric Synthetic Aperture Radar,InSAR)在监测缓慢蠕动形变方面具有独特优势,其可利用历史存档数据获取地表变化区域的位置、变形趋势及特征,在火山形变监测中也已得到广泛应用(季灵运等,2013; 郭倩,2020; 杨云飞,2021; 刘媛媛等,2022; Liuetal.,2023)。不同学者采用InSAR技术对长白山天池火山开展了地表形变监测工作。1992—1998年,Kim等(2004)基于InSAR技术监测到长白山天池火山周围呈现速度为3mm/a的缓慢隆升变形。2004—2010年期间,InSAR技术探测到在长白山天池火山的北坡和西坡有一定程度的地表变形。北坡的监测结果显示,2004—2006年,垂直形变量≤5mm,形变量级较小; 2006—2009年,地表出现明显的隆升趋势,在破火山口边缘监测到12mm的隆升变形; 2009—2010年,开始发生明显下沉,量级达10mm。西坡的监测结果表明,2004—2008年,地表呈现隆升态势,之后整体发生缓慢下沉(季灵运,2012; 何平等,2015)。2018—2020年,Trasatti等(2021)发现长白山天池火山的东南侧和西南侧均发生明显变形,东南侧的形变速率可达20mm/a,西南侧的形变速率达-20mm/a。上述监测结果表明长白山天池火山持续发生着地表变形,火山下方可能持续存在岩浆活动。为了解长白山天池火山的活动状态,本文采用时序InSAR技术,利用2015—2022年间的Sentinel-1A/B升、降轨数据,对长白山天池火山开展了动态监测,并利用Mogi点源模型反演火山岩浆房的深度与体积变化,以期为天池火山的危险性研究提供更多数据支撑。

1 技术原理

1.1 时序InSAR技术

本文采用SBAS InSAR(Small Baseline Subset InSAR)技术对长白山天池火山开展监测。它是一种针对分布式散射目标的时序InSAR技术(Berardinoetal.,2002),使用多干涉对组合进行高相干点的相位分析,以获取可靠的连续面状结果,对于相干性比较弱的区域也能够获取较好的监测效果。SBAS InSAR基于多主影像的SAR数据处理模式,可以在很大程度上保证2幅干涉影像之间的相干性,增加冗余观测量并提高影像利用率。其基本原理是假设区域有N张SAR影像,根据小基线原理,对满足要求的影像进行干涉。在去除平地效应和地形相位后生成差分干涉图,采用最小费用流算法对每幅差分干涉图进行相位解缠。第i幅干涉图中像素(x,y)的解缠相位δφ(x,y)可表示为

δφ(x,y)=δφ(x,y)(tB)-δφ(x,y)(tA)≈δφ(def,x,y)+δφ(e,x,y)+δφ(α,x,y)+δφ(n,x,y)

(1)

式中,δφ(x,y)(tB)与δφ(x,y)(tA)分别表示SAR采集时间tB、tA时的相位(tA

1.2 Mogi点源模型

Mogi点源模型也被称为火山模型,是反演火山岩浆房参数最常用的模型(陈国浒等,2008; 靳锡波等,2020)。其原理是将火山内部的岩浆压力源与火山活动引起的地表形变建立联系的一种数学模型。该模型有2个基本假设(Mogi,1958):1)地壳是一个无限大的弹性体; 2)地表会由于压力源的压力或体积变化而产生位移。当地下岩浆发生活动后,会造成工作面上的介质发生变化,破坏平衡状态,导致地表发生形变。基于此假设,使用地表形变作为约束,可反演获得火山岩浆房的体积变化及位置分布情况。

2 数据及处理流程

2.1 数据

本文围绕长白山天池火山建立研究区,从欧洲航天局(1)https:∥search.asf.alaska.edu/。获取了64景升轨Sentinel-1A影像; 2个轨道的降轨SAR影像,其中轨道号为32的数据有172景,轨道号为134的数据有180景。考虑到长白山冬季积雪会影响InSAR形变监测,选择影像时主要考虑每年6—9月的数据。数据覆盖情况如图1所示,基本参数见表1。同时,采用精密轨道数据去除轨道误差(2)https:∥s1qc.asf.alaska.edu/aux＿poeorb/。,并采用30m分辨率的SRTM DEM地形数据(3)https:∥e4ftl01.cr.usgs.gov/。模拟并去除地形相位。

表1 Sentinel-1A数据参数列表

图1 研究区的地理位置及SAR影像覆盖图

2.2 处理流程

本文使用GAMMA商业软件,采用SBAS InSAR技术,在保证空间分辨率的基础上减少噪声的干扰,获取了连续、面状的形变时间序列监测结果。在处理过程中,为了降低噪声、提高信噪比,首先设置距离向与方位向之比为5︰1进行多视处理,并设置时间与空间基线阈值,满足阈值的影像生成干涉对。由于获取的升轨的数据量较少,设置的时间基线为1500d,空间基线为±250m,最后参与形变计算的干涉对组合如图2a 所示。设置降轨数据的时间、空间基线分别为500d、±250m,2个轨道的干涉图基线组合分别如图2b、2c所示。然后,借助30m分辨率的SRTM DEM模拟并去除地形相位,获取差分干涉图。采用Goldstein滤波,设置大小为64、32、16的窗口相继进行滤波,以提高干涉图的相干性并降低噪声。采用最小费用流方法进行相位解缠,设置相干性阈值,剔除大部分由于积雪及植被覆盖导致相干性较低的像元,分别获取升、降轨的解缠干涉图。采用通用型卫星雷达大气改正系统(Generic Atmospheric Correction Online Service for InSAR,GACOS)(4)http:∥www.gacos.net。去除大气误差,通过使用二次多项式拟合的方法去除趋势误差。对去除误差后的解缠干涉图进行挑选,剔除质量较差的干涉图。最后,采用相位堆叠算法(Stacking InSAR)获取年平均形变速率图,并进行地理编码,获得WGS-84坐标系下的形变结果。

图2 干涉图的时空基线图

获取InSAR的形变速率后,基于Mogi点源模型反演火山的岩浆房参数,得到岩浆房的深度及体积变化。在计算过程中,假设岩浆压力源的参考点位置为(42.0075°N,128.0559°E),使用贝叶斯方法反演模型参数,设置迭代次数为106次,以确定岩浆房的最佳参数及不确定性。采用马尔科夫链-蒙特卡洛算法结合Metropolis-Hastings算法计算参数的后验概率分布。反演输入的约束数据为InSAR视线向形变速率,在开始反演计算之前,采用均匀降采样方法对数据进行重采样,以提高反演效率。然后通过设置岩浆房的中心位置x、y,深度z,体积变化率ΔV的上、下限及相应的步长,约束模型参数的反演。

3 火山形变场

3.1 2015—2022年InSAR形变速率

本文获得的长白山天池火山的InSAR形变速率如图3所示。其中,正值表示形变靠近卫星飞行方向,负值表示形变远离卫星飞行方向。天池火山在中国境内的区域植被较为茂密,干涉图的相干性较差,靠近朝鲜的区域植被较少,且火山口附近裸岩较多,获取的干涉点主要分布在火山口附近。一般而言,岩浆活动引起的地表变形表现为整体的趋势性形变,局部形变信号通常不是由火山运动引起的,可能是由人类活动或其他因素导致。

图3 形变速率图

在2015—2022年期间,升、降轨监测到的年平均形变速率为-8～6mm/a。升轨的形变比降轨的量级更小,主要是由于SAR的几何特征及地形差异所致。升、降轨的形变方向较为一致,说明长白山天池火山的形变主要发生在垂直方向上。升轨视线向形变速率较小(图3a),主要的形变分布在破火山口处,形变分布均匀,地表变形主要为沉降,形变速率为-2～-1mm/a。2个降轨数据的形变分布表现较为一致(图3b,c),在火山口附近出现整体的趋势性形变,在远离火山口的区域存在局部变形。火山口附近表现为整体缓慢下沉,量级不大,视线向形变速率为-4～-2mm/a。局部形变主要分布在火山口东侧的2个区域,用A、B表示,其中A区域的最大形变速率为-8mm/a,B区域的最大形变速率为-6mm/a。根据形变速率结果可知,天池火山总体表现为沉降状态,火山口形变较大,向四周逐渐减小,呈扩散状,表明天池火山可能存在深部活动。

由于未获得实地测量数据,本文采用内符合精度评定的方法评价InSAR监测结果的精度。结合SAR影像的局部入射角,求取了3个轨道的SAR数据在垂直方向上的形变速率,分别统计不同轨道间公共点的垂直形变的差值,计算残差的偏差及标准差(图4)。经统计,残差总体均呈正态分布,升轨(ASC54)与降轨(DSC32)垂向形变的残差偏差为0.78mm/a(图4a),但在正值区间存在一定残差。降轨垂向形变(DSC32、DSC134)的残差偏差为-1.74mm/a(图4b),这可能是由于参考点的差异导致结果存在较小的系统偏差,但残差分布更为集中。不同轨道数据间垂向形变的残差均为毫米级,达到InSAR监测的精度要求,说明本文使用InSAR技术监测长白山天池火山的结果具有较好的可靠性。

图4 观测值间的残差分布直方图

3.2 2015—2022年InSAR形变时间序列

为进一步了解长白山天池火山各部位形变的时空变化特征及趋势,使用SBAS InSAR技术进行计算,获取了时间段为2014年11月—2021年12月、共计172景降轨(DSC32)SAR影像的形变时间序列。图5为降轨数据(DSC32)在雷达视线方向上的累计形变时间序列。从图5可以看出,长白山天池火山在监测期间存在持续变形。2014年11月—2017年6月,火山总体形变速度较慢,从2018年6月起形变加速; 2020年7月—2021年12月期间,形变加速明显。在近7a的监测时间内,累计形变量达-40mm。总体而言,近7a以来,长白山天池火山地表持续发生缓慢变形,主要发生在火山口附近,且于近期有加速趋势。

图5 降轨(DSC32)形变时间序列

提取横穿火山口的2条剖线,定量分析长白山天池火山的地表变化情况。剖线方向分别为NW-SE和NE-SW,具体位置见图5中的2021年12月11日累计形变图。获取的剖线高程及形变时间序列如图6所示。长白山天池火山水面的位置相对较低,水面高约2200m,高度由中间向两侧增加,在2条剖线的时间序列监测结果中,天池火山口处存在明显变形。

图6 剖面累计形变图

图7 长白山火山地震数据统计图(改自Meng et al.,2022)

根据2条剖线的累计形变时间序列可知,在天池火山口的南侧监测到较大程度的变形。根据NW-SE剖线的形变时间序列可知,在火山口东南侧存在显著形变,视线向累计形变量达-38mm,发生在天池火山口的内侧。在火山口的北西侧存在小幅度地表隆升。火山附近总体形变主要表现为沉降,距离火山口的位置越远,形变变化趋势越趋于稳定。形变与地形表现出相关性,在地形起伏剧烈的区域,对应的形变也呈现出一定程度的起伏。根据NE-SW剖线的形变时间序列可知,火山口西南侧的形变比东北侧量级更大,量级最大的形变出现在火山口山体外侧的西南侧斜坡中部,视线向形变量达-38mm。

根据前人的研究可知,长白山天池火山的地表形变在历史监测中表现出周期性变化。同时,结合地震活动性可以了解火山的活动状态。1992—1998年,该区地表没有发生明显变形,地震活动弱,火山岩浆活动并不活跃; 2002—2006年,地表发生约6.5cm的隆升变形(Jietal.,2021),同时根据得到的地震数据表明,此时段内中强地震活动性增强,最大震级为4.4级,地震发生的频率约为600次/a,表明该时段内火山岩浆活动较为活跃(Panetal.,2021; Mengetal.,2022)。在接近4a的隆升变化之后,2008—2011年观测到天池火山处地表开始发生沉降,并且地震活动减弱,发生地震的频率为约70次/a,此阶段火山岩浆活动较弱(Panetal.,2021)。本文的研究表明,2015—2022年期间,火山口附近的地表视线向年平均形变速率约为-4～-2mm/a,累计形变量可达-40mm。同时,统计得到的震中距火山口30km内的地震数量显示,2015—2020年共发生地震195次,表明该时段内地震活动性弱,最大地震的震级为2.2级,说明在此期间火山岩浆活动并不活跃。但值得注意的是,2020年12月22日出现了由38次小地震组成的火山震群事件,地震活动有增多的趋势(盘晓东等,2022)。通常情况下,火山扰动被认为是周期性的岩浆活动。目前,长白山天池火山虽然没有近期喷发的危险,但仍需要持续开展监测。

4 岩浆房参数反演

上文通过近7a的SAR影像获取的InSAR形变监测结果,反映了长白山天池火山的地表动态变化过程,为研究其岩浆活动变化奠定了基础。本文使用Mogi模型,联合2个降轨形变场数据反演了岩浆房参数。在InSAR结果中2个降轨数据在形变分布特征上表现出较好的一致性,表明监测结果的可靠性更高。同时,由于使用升轨数据计算形变结果时SAR影像数据量少,干涉对的时空基线较长,可能导致干涉图之间的相干性相对较弱,使结果中存在一定误差。因此,为保证结果的准确性,在进行参数反演时,主要利用2个降轨数据获取的地表形变结果作为约束,通过多次反演实验,得到拟合残差最小的岩浆房参数(表2)。

表2 Mogi点源模型反演的几何参数

根据反演得到的参数可以看出,岩浆房的位置处于火山口下偏西的区域(41.9928°N,128.0361°E),深约6km,岩浆房的年体积变化率为-3.3×105m3/a,表明长白山天池火山的岩浆在近几年变化缓慢,体积变化表现为收缩。不同学者对长白山天池火山的岩浆房空间位置及活动性开展了研究,并给出了不同结果。地球物理探测的相关研究者发现在长白山天池火山的周围存在明显的P波低速异常,认为在下地壳或更深的位置存在岩浆房,且处于活动状态,并进一步估计了岩浆房的位置和深度。其中,基于三维层析成像技术推测出的长白山天池火山岩浆房位于火山口西侧位置(杨卓欣等,2005); 利用小地震精定位数据反演结果推测得到的岩浆房位于天池火山地下约5km处(吴建平等,2007); 通过钻孔测量技术并基于长白山天池火山地温梯度判断的岩浆房位于深5～8km处(潘波等,2017; 钱程等,2022)。大地电磁观测的反演结果表明,长白山天池火山下方6～8km处存在岩浆通道,在火山口的北部7km处存在岩浆囊(仇根根等,2014)。大地测量的相关学者通过地表监测数据反演获取了长白山天池火山的岩浆房分布及岩浆体积的变化速率,结果表明,岩浆房并非位于火山的正下方,岩浆活动会引起地表的形变。其中,InSAR反演结果显示长白山天池火山的岩浆房位于火山口下偏西的区域,深6～7km(何平等,2015),与本文反演的结果相近。GNSS与InSAR结合的反演结果显示,长白山天池火山的岩浆房位于深6.9km处(陈国浒等,2008)。以上研究表明,长白山天池火山区下方存在岩浆房,深度为5～8km,由于观测方式、数据处理手段、观测时间段的不同及不同空间范围内的岩浆扰动均会导致反演的岩浆房深度出现偏差,因此,对于岩浆房所处的空间位置与展布形态目前还不够清晰。本文的观测结果与GNSS、地球物理探测等技术得到的结果相比,能够更加全面地覆盖天池火山靠近朝鲜一侧的区域,更有利于说明岩浆房的空间位置与活动状态。

本研究通过反演的参数正演出长白山天池火山的地表形变,并计算模拟值与观测值的残差,以验证模拟的准确性。图8给出了各数据的观测值、模拟值与残差值。可以看出,在反演结果中,趋势性形变均被较好拟合,说明岩浆房活动变形引起的地表变化可被反演获取,表明反演结果具有可靠性。同时,我们对反演的残差进行了统计(图9),发现残差均分布在零值附近,DSC32残差的偏差为0.13mm/a,DSC134残差的偏差为0.23mm/a。残差结果在一定程度上说明了数据拟合较好,但也可以看到在负值区间存在一定残差,这可能是Mogi模型反演时无法拟合到局部形变所致。

图8 反演得到的观测值、模拟值与残差值图

图9 残差统计图

5 讨论

5.1 火山岩浆房收缩机理分析

长白山天池火山被普遍认为是爆裂式喷发型火山,根据目前多种资料佐证,在长白山天池火山地下的浅地壳中存在岩浆房,且在地下深部存在岩浆储存库(许建东等,2022)。浅部岩浆房为位于地壳中的部分或全部熔融体,其来源多是由深部岩浆储存库的岩浆多次注入。火山喷发通常是浅部岩浆房中的岩浆在浮力的驱动下通过火山通道涌出地表的过程(段政等,2022)。一般情况下,岩浆房处于平衡状态,当有新的岩浆注入时,岩浆房承受的压力发生变化,呈现增加趋势,岩浆房的晶体平衡被打破,从而产生新的反应,使得岩浆处于活动状态。因此,监测岩浆房体积的变化情况有助于了解火山活动的强弱。岩浆房体积变化一般有膨胀与收缩2种类型。火山岩浆房的体积膨胀一般由岩浆入侵或增压导致(Luetal.,2003; Jietal.,2013)。体积收缩的原因较为复杂,可能由以下几种情况导致:1)岩浆体的缓慢冷却和结晶(Luetal.,2005); 2)火山气体由火山喷气孔或温泉释放,导致压力降低(Sepeetal.,2007; Fournier,2008; Jietal.,2018); 3)浅储层的岩浆流出(Cervelletal.,2002); 4)火山强烈活动、喷发而导致岩浆损失(Luetal.,2000)。本研究获取的长白山天池火山在监测时间段内的地表形变主要表现为沉降,部分区域存在小幅度隆升,反演得到的岩浆房体积变化表现为收缩,在监测期间未发生喷发,推测体积的变化很可能是由于岩浆的冷却和结晶所致。

5.2 火山形变与岩浆活动时间序列

长白山天池火山被认为是中国最具危险性的活火山(洪汉净等,2007)。火山变形可以揭示火山岩浆是否发生扰动,根据火山地表形变研究,长白山天池火山表现出持续的岩浆活动。针对长白山天池火山的变形监测,主要包括水准测量、GPS、InSAR 3种方式。火山周边共建有15个GPS站,覆盖面积达1200km2,测量始于2000年。精密水准测量是从2002年起在北坡和西坡进行的。InSAR技术利用历史影像追溯地表变形情况,最早于1992年开始监测,随着目前影像种类和数量的增多,在长白山天池火山变形监测中起到了关键的作用。

多项研究表明,自1992年至今,长白山天池火山地下岩浆房的演化过程为平静—发生扰动—归于平静(图10)。1992—1998年期间,长白山天池火山没有发生明显变形(Jietal.,2021)。1995—1998年间,利用InSAR技术结合点源模型的反演结果显示,火山地下可能存在双岩浆房,深度分别为7.9km、5.5km,并得到岩浆房的体积变化量分别为1.25×106m3、3.6×106m3,体积变化表现为膨胀(陈国浒等,2008)。2002—2010年,长白山天池火山地区的变形监测数据更加丰富。2002—2006年,GPS与InSAR监测到长白山天池火山附近的岩浆发生动荡,岩浆活动频繁,出现明显的异常信号。2002—2005年,采用GPS监测数据结合点源模型的反演结果显示,岩浆房的体积变化率为4.6×106m3/a,分布在深4.4km处,体积变化表现为膨胀(Xuetal.,2012; Trasattietal.,2021)。2006—2009年,InSAR监测结合模型的反演结果显示,岩浆房分布在6.3km深度处,体积变化率为1.11×106m3/a,体积变化仍为膨胀,但变化速率出现明显下降,说明自2006年起岩浆活动变弱,开始变得平静(何平等,2015)。2009—2011年,GPS监测到火山出现明显下沉,结合点源模型的反演结果显示,岩浆房以-1.4×106m3/a的速度发生收缩(Xuetal.,2012; Trasattietal.,2021)。2015—2022年,InSAR监测到火山仍在下沉,但速度较为缓慢,结合点源模型的反演结果显示,深度为6km处的岩浆房以-0.33×106m3/a的速度发生收缩,推测长白山天池火山继2002—2006年的扰动期过后一直处于相对较为平静的状态。