基于D-InSAR的樱岛火山喷发 形变场解算与反演分析

刘媛媛 赵振宇 晏 霞 陈人杰

1 东华理工大学测绘工程学院，南昌市广兰大道418号，330013 2 自然资源部环鄱阳湖区域矿山环境监测与治理重点实验室，南昌市广兰大道418号，330013

火山喷发是一种重大自然灾害。据统计，2020年全球共有24个国家的89座火山发生喷发，其中Ⅲ级以上预警的火山达21座[1]。火山喷发时会释放出大量的火山灰和火山气体，这些物质会随大气环流被带至高空，从而影响气候变化[2]。樱岛火山(31°35′N，130°39′E)位于日本鹿儿岛市，东接大隅群岛，西隔锦江湾，距离鹿儿岛市中心约4 km，南北面皆为鹿儿岛湾。火山由北岳(海拔1 117 m)、中岳(海拔1 160 m)、南岳(海拔1 140 m)组成，面积约77 km2，是一座活跃的成层火山[3]，包含南岳火山口、昭和火山口及中央火山锥两侧的多个寄生火山口[4]。樱岛火山在过去的100 a里发生过几次不同类型的喷发。根据日本气象厅公布的樱岛火山活动情况，2020-06下旬开始观测到山体膨胀的缓慢变动，2020-07下旬膨胀停止，随后在2020-08-09 05:41:53监测到南岳火山口发生喷发，当日发生的A型地震达68次，火山喷发产生的烟雾高度达5 km，樱岛部分地区的降灰量达300 g/m2。

火山喷发往往是由地壳内部岩浆活动引起的，根据地表形变可以有效反演出地下岩浆的活动特征，对预测火山喷发及灾害防御有重要意义。传统的监测手段无法及时获取大范围的地表形变信息，而合成孔径雷达差分干涉D-InSAR技术作为一种非接触测量手段，精度可达cm级甚至亚cm级，具有全天时、视域广、密度大等特点[5]，已被广泛应用于各类火山喷发引起的地表形变监测中[6]。2015年樱岛火山岩脉入侵事件中，相关学者使用多视角InSAR数据进行D-InSAR测量，解算岩脉入侵事件导致的三维形变场[7]。为监测樱岛火山的最新活动状态，本文拟采用D-InSAR技术，利用覆盖日本鹿儿岛市樱岛火山的2景Sentinel-1A数据研究2020-08-09樱岛火山喷发事件，并利用点源Mogi模型反演火山下方岩浆源的特征。本文对研究火山喷发的性质和灾害分布具有重要意义。

1 技术原理

1.1 D-InSAR技术

InSAR技术的重要应用领域之一是地表形变监测[8]，而D-InSAR技术是InSAR技术的拓展，可以理解为SAR卫星沿着相同的轨道重复运行，获取同一地区不同时间的影像。将2幅不同时段主、副影像的相位值进行差分处理，即可得到差分干涉相位图。此时干涉图中包含的相位有参考椭球面相位、地形相位、形变相位、大气相位和噪声相位，其中占比较大的相位为参考椭球面相位、地形相位和形变相位，D-InSAR主要提取的是形变相位。

1.2 Mogi模型

对此次火山喷发事件进行多次模拟实验后，本文最终选用点源Mogi模型[9]来恢复地表形变。Mogi模型的参数分别为源中心坐标在地表的投影(x0，y0)，源中心的深度d以及岩浆源体积变化量ΔV。假设地面坐标系为(x，y，z)，岩浆源的空间位置为(x0，y0，d)，根据Mogi模型拟合地面点三维位移与源模型参数的公式为：

(1)

(2)

(3)

r2=(x-x0)2+(y-y0)2

(4)

式中，Ux、Uy、Uz为坐标原点的位移分量，r为Mogi源模型中心与地面点在地表投影坐标间的距离。

2 数据处理与模型反演

2.1 数据处理

选取覆盖樱岛火山区域2020-07-28(喷发前)与2020-08-09(喷发后)21h的2景C波段降轨Sentinel-lA单视复数SLC影像数据以及对应的精密轨道文件，时间基线为12 d。同时，选用分辨率高达30 m的DEM数据，用于减弱外部DEM引起的误差。

由于火山周围生长的植被会造成失相干，为抑制数据噪声的影响以及提高相干性，本文采用距离向为8、方位向为2的多视比，并使用一种局部自适应滤波的方法。在此基础上，采用基于规则格网的最小费用流算法进行相位解缠，相关系数阈值设为0.25，最终通过地理编码获得地理坐标系下樱岛火山LOS向一维形变场(图1)。由图可见，此次火山喷发事件造成的形变主要集中在中央火山锥中心，地表沉降主要发生在火山锥中心地带，最大沉降量与平均沉降量分别为5.5 cm和2.85 cm；抬升主要发生在火山边缘区域，最大抬升量与平均抬升量分别为5 cm与2.24 cm，其中北部的隆起可能与艾拉火山口下方的岩浆活动有关。

图1 LOS向形变

2.2 数据误差估计

通常情况下，SAR卫星在对观测区域进行数据采集时会受到大气相位的影响，与大气和地面状态相关的每个InSAR图像都具有特定的InSAR数据误差[9]。本文构建噪声的完整协方差矩阵：首先从原始数据中删除一个线性斜坡，以校正可能的长波大气效应；然后估计去趋势数据集上的实验半变异函数；最后使用nugget模型拟合出指数函数。指数函数的方差与协方差矩阵的定义为：

γ(h)=nugget+

(5)

式中，γ(h)为每个数据点之间任意给定距离h处的协方差矩阵，range为数据点在空间上的相关距离，nugget为空间独立噪声的级别，sill为当range趋于无穷大时半方差的最大值。图2为根据式(5)计算出的InSAR数据误差特性，由图可见，得到的距离值、噪声级别和方差最大值分别为849.845 8 m、1.612 2×10-6m2和8.964 7×10-6m2。

图2 InSAR数据估计非变形区域的误差特征

2.3 模型反演与分析

本文以D-InSAR处理得到的位移信息为观测数据，使用贝叶斯方法检索樱岛火山喷发的Mogi源最佳参数和不确定性，并使用马尔可夫链蒙特卡洛(MCMC)算法结合M-H算法搜索Mogi 模型各项参数的后验概率[10]。

反演中使用的参考系坐标原点为南岳火山口地理坐标(130.66°E，31.58°N)，X为东西方向的位移，向东为“+”；Y为南北方向的位移，向北为“+”；Z为岩浆源的深度。由图1可见，火山周围大量的植被会造成几何失真产生空白，故选择中央火山锥区域进行模型反演。反演开始前，通过设定模型参数的约束信息来提高反演速度(表1)。此外，采用D-InSAR得到的数据集为面状数据，原始数据量较大，需要在反演的初始阶段运用四叉树算法对其进行降采样处理，采样后的数据点为106个，保证了形变场的特征。

表1 Mogi模型参数的约束信息及反演结果

经过1.0×106次迭代后得到模型的最佳拟合参数，由表1可知，东西向位移为0.402 km，南北向位移为0.069 km，深度为1.061 km，岩浆源位于南岳与昭和火山口之间，比较靠近昭和火山口。模型各项参数的联合概率如图3所示，由图可见，总体残差很小。

①X/m, ②Y/m, ③Z/103 m, ④ΔV/105m3, ⑤视线方向恒定偏移量/m， ⑥轨道残差X/10-6 m， ⑦轨道残差Y/10-6 m，⑧迭代次数/104次图3 联合概率分布

Iguchi等[11]研究樱岛昭和火山口2007～2011年火山喷发的构造活动，在研究中假设一个位于南岳火山口下方4 km处的压力源，利用倾斜仪和引伸仪的径向应变变化，估计压力源体积的月变化量与火山灰重量之间的关系(图4)。结果表明，昭和火山口可能通往南岳下方的岩浆储层或与南岳火山口的管道相连。由图4可见，Mogi源的月体积变化主要集中在±0.5×106m3范围内，而根据Mogi源反演得到的此次火山喷发的体积参数为-0.139×106m3，与Iguchi等[11]得到的月体积变化量较为接近。

图4 压力源体积变化与火山灰的关系[11]

2016年Hotta等[12]使用GNSS、应变与倾斜数据构建一个3Mogi源模型，用于计算以往多次火山喷发事件的体积变化量。表2为文献[12]与本研究计算的M源体积变化量，表3为M源参数，图5为压力源平面位置。由表2可见，以往火山喷发事件导致的M源体积最大变化量与最小变化量分别为-0.33×106m3和-0.017×106m3，此次火山喷发事件的体积变化量大约是以往M源体积最小变化量的8倍和最大变化量的0.42倍。由图5可见，Mogi源分别位于艾拉火山口下方(A源)9.6 km、北岳下方(K源)3.3 km和南岳下方(M源)0.7 km处。此次火山喷发事件反演得到的Mogi源位置与Hotta等[12]得到的M源位置较为接近，M1源由本文计算得到，M2源为Hotta等[12]计算得到。M1位于南岳火山口东侧，M2位于昭和火山口东侧，二者在位置上存在一定差距，可能是由于SAR卫星在监测过程中容易受到误差的影响，虽然本文在实验过程中对误差进行了改正和削弱，但无法完全消除。尽管如此，从最终结果来看，本文反演的参数结果仍然具有较好的参考价值。

表2 M源体积变化量[12]

表3 M源参数

图5 压力源平面位置

火山活动的强弱会随时间变化。樱岛火山的常见喷发模式有2种：间歇性喷发和连续性喷发。6月底至7月底樱岛火山无喷发事件，但监测到山体膨胀，最终于2020-08-09发生火山喷发，因此可以将其视作一次短期的间歇性喷发。在此期间火山下方的岩浆由A源入侵至K源，经过一段时间的累积后由K源迁移至M源，最终在地表喷发。樱岛的岩浆管道系统由一个主要岩浆储层、一个次要岩浆储层和一个连接管道组成，在火山喷发后，岩浆各个储层都会降低，这可能会导致岩浆撤退和火山体积的收缩[13]。此次火山喷发事件反演得到的收缩源与Hotta等[12]得到的南岳源位置较为接近，由此可知，樱岛的岩浆管道系统组成为：艾拉火山口下方约10 km的岩浆储层连接北岳下方3.3 km的岩浆储层，经过南岳下方1.016 km的浅层储层，通过管道连接至山顶火山口(图6)，其中北岳下方的K源位置由Hotta等[12]提供。

图6 地下压力源与岩浆供应

3 结 语

1)最大沉降区域为火山锥中心地带，最大沉降量与平均沉降量分别为5.5 cm和2.85 cm；抬升主要发生在火山边缘区域，最大抬升量与平均抬升量分别为5 cm和2.24 cm。

2)根据Mogi模型反演得到的岩浆源深度和体积收缩量分别为1.016 km和-0.139×106m3，岩浆源位于南岳下方，与南岳火山口的喷发活动有关。

樱岛火山是一座活火山，近年来南岳火山口和昭和火山口处于活跃状态，同时也伴随着密集的地震与火山喷发现象，间歇性与连续性的喷发影响着居民的正常生活，因此许多学者致力于研究此火山的构造活动。然而在中央火山锥两侧还存在寄生火山口，因此仍然无法对具体的岩浆供应路径作出详细解释。为了更好地理解火山活动的全过程，还需要利用更加全面的数据进行研究与讨论。