2022, Vol. 42
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模型拟合地面点三维位移与源模型参数的公式为：
| $ {U_x} = \frac{{\Delta V\left( {x - {x_0}} \right)}}{{{{\left( {{d^2} + {r^2}} \right)}^{\frac{3}{2}}}}} $ | (1) |
| $ {U_y} = \frac{{\Delta V\left( {y - {y_0}} \right)}}{{{{\left( {{d^2} + {r^2}} \right)}^{\frac{3}{2}}}}} $ | (2) |
| $ {U_z} = \frac{{r\Delta V}}{{{{\left( {{d^2} + {r^2}} \right)}^{\frac{3}{2}}}}} $ | (3) |
| $ {r^2} = {\left( {x - {x_0}} \right)^2} + {\left( {y - {y_0}} \right)^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向形变 Fig. 1 LOS deformation |
通常情况下，SAR卫星在对观测区域进行数据采集时会受到大气相位的影响，与大气和地面状态相关的每个InSAR图像都具有特定的InSAR数据误差[9]。本文构建噪声的完整协方差矩阵：首先从原始数据中删除一个线性斜坡，以校正可能的长波大气效应；然后估计去趋势数据集上的实验半变异函数；最后使用nugget模型拟合出指数函数。指数函数的方差与协方差矩阵的定义为：
| $ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\gamma }}(h) = {\rm{ nugget }} + }\\ {({\rm{ sill - nugget }})\left( {1 - \exp \left( {\frac{{ - h}}{{{\rm{ range }}}}} \right)} \right)} \end{array} $ | (5) |
式中，γ(h)为每个数据点之间任意给定距离h处的协方差矩阵，range为数据点在空间上的相关距离，nugget为空间独立噪声的级别，sill为当range趋于无穷大时半方差的最大值。图 2为根据式(5)计算出的InSAR数据误差特性，由图可见，得到的距离值、噪声级别和方差最大值分别为849.845 8 m、1.612 2×10-6 m2和8.964 7×10-6 m2。
|
图 2 InSAR数据估计非变形区域的误差特征 Fig. 2 Error characteristics of InSAR data estimation in non-deformation region |
本文以D-InSAR处理得到的位移信息为观测数据，使用贝叶斯方法检索樱岛火山喷发的Mogi源最佳参数和不确定性，并使用马尔可夫链蒙特卡洛(MCMC)算法结合M-H算法搜索Mogi模型各项参数的后验概率[10]。
反演中使用的参考系坐标原点为南岳火山口地理坐标(130.66°E，31.58°N)，X为东西方向的位移，向东为“+”；Y为南北方向的位移，向北为“+”；Z为岩浆源的深度。由图 1可见，火山周围大量的植被会造成几何失真产生空白，故选择中央火山锥区域进行模型反演。反演开始前，通过设定模型参数的约束信息来提高反演速度(表 1)。此外，采用D-InSAR得到的数据集为面状数据，原始数据量较大，需要在反演的初始阶段运用四叉树算法对其进行降采样处理，采样后的数据点为106个，保证了形变场的特征。
|
表 1 Mogi模型参数的约束信息及反演结果 Tab. 1 Constraint information and inversion results of Mogi model parameters |
经过1.0×106次迭代后得到模型的最佳拟合参数，由表 1可知，东西向位移为0.402 km，南北向位移为0.069 km，深度为1.061 km，岩浆源位于南岳与昭和火山口之间，比较靠近昭和火山口。模型各项参数的联合概率如图 3所示，由图可见，总体残差很小。
|
①X/m, ②Y/m, ③Z/103 m, ④△V/105 m3, ⑤视线方向恒定偏移量/m, ⑥轨道残差X/10-6 m, ⑦轨道残差Y/10-6 m, ⑧迭代次数/104次 图 3 联合概率分布 Fig. 3 Joint probability distribution |
Iguchi等[11]研究樱岛昭和火山口2007~2011年火山喷发的构造活动，在研究中假设一个位于南岳火山口下方4 km处的压力源，利用倾斜仪和引伸仪的径向应变变化，估计压力源体积的月变化量与火山灰重量之间的关系(图 4)。结果表明，昭和火山口可能通往南岳下方的岩浆储层或与南岳火山口的管道相连。由图 4可见，Mogi源的月体积变化主要集中在±0.5×106 m3范围内，而根据Mogi源反演得到的此次火山喷发的体积参数为-0.139×106 m3，与Iguchi等[11]得到的月体积变化量较为接近。
|
图 4 压力源体积变化与火山灰的关系[11] Fig. 4 Relationship between volume change of pressure source and volcanic ash[11] |
2016年Hotta等[12]使用GNSS、应变与倾斜数据构建一个3Mogi源模型，用于计算以往多次火山喷发事件的体积变化量。表 2为文献[12]与本研究计算的M源体积变化量，表 3为M源参数，图 5为压力源平面位置。由表 2可见，以往火山喷发事件导致的M源体积最大变化量与最小变化量分别为-0.33×106 m3和-0.017×106 m3，此次火山喷发事件的体积变化量大约是以往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卫星在监测过程中容易受到误差的影响，虽然本文在实验过程中对误差进行了改正和削弱，但无法完全消除。尽管如此，从最终结果来看，本文反演的参数结果仍然具有较好的参考价值。
|
|
表 3 M源参数 Tab. 3 M source parameters |
|
图 5 压力源平面位置 Fig. 5 Plane position of pressure source |
火山活动的强弱会随时间变化。樱岛火山的常见喷发模式有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 地下压力源与岩浆供应 Fig. 6 Underground pressure source and magma supply |
1) 最大沉降区域为火山锥中心地带，最大沉降量与平均沉降量分别为5.5 cm和2.85 cm；抬升主要发生在火山边缘区域，最大抬升量与平均抬升量分别为5 cm和2.24 cm。
2) 根据Mogi模型反演得到的岩浆源深度和体积收缩量分别为1.016 km和-0.139×106 m3，岩浆源位于南岳下方，与南岳火山口的喷发活动有关。
樱岛火山是一座活火山，近年来南岳火山口和昭和火山口处于活跃状态，同时也伴随着密集的地震与火山喷发现象，间歇性与连续性的喷发影响着居民的正常生活，因此许多学者致力于研究此火山的构造活动。然而在中央火山锥两侧还存在寄生火山口，因此仍然无法对具体的岩浆供应路径作出详细解释。为了更好地理解火山活动的全过程，还需要利用更加全面的数据进行研究与讨论。
| [1] |
孙立影, 盘晓东, 朱大庆. 2020年全球火山活动综述[J]. 城市与减灾, 2021, 138(3): 22-28 (Sun Liying, Pan Xiaodong, Zhu Daqing. Overview of Global Volcanic Activities in 2020[J]. City and Disaster Reduction, 2021, 138(3): 22-28)
(  0) |
| [2] |
李传金, 任贾文, 秦大河, 等. 火山活动的气候影响及其冰芯记录研究进展[J]. 冰川冻土, 2012, 34(4): 863-876 (Li Chuanjin, Ren Jiawen, Qin Dahe, et al. Summary of Research on Climatic Influences from Volcanic Activities and Depositional Records of Volcanic Matters in Ice Cores[J]. Journal of Glaciology and Geocryology, 2012, 34(4): 863-876)
( 0) |
| [3] |
于红梅. 火山分类[J]. 城市与减灾, 2018, 122(5): 12-17 (Yu Hongmei. Classification of Volcanoes[J]. City and Disaster Reduction, 2018, 122(5): 12-17)
( 0) |
| [4] |
Yamamoto K, Sonoda T, Takayama T, et al. Vertical Ground Deformation Associated with the Volcanic Activity of Sakurajima Volcano, Japan during 1996-2010 as Revealed by Repeated Precise Leveling Surveys[J]. Bulletin of the Volcanological Social of Japan, 2013, 58(1): 137-151
( 0) |
| [5] |
朱建军, 李志伟, 胡俊. InSAR变形监测方法与研究进展[J]. 测绘学报, 2017, 46(10): 1 717-1 733 (Zhu Jianjun, Li Zhiwei, Hu Jun. Research Progress and Methods of InSAR for Deformation Monitoring[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(10): 1 717-1 733)
( 0) |
| [6] |
史国强, 何秀凤, 肖儒雅. 基于差分干涉与GPS测量的夏威夷基拉韦厄火山形变监测[J]. 大地测量与地球动力学, 2016, 36(5): 404-407 (Shi Guoqiang, He Xiufeng, Xiao Ruya. Monitoring Surface Deformation of Kilauea Volcano, Hawaii, Based on D-InSAR and GPS Measurements[J]. Journal of Geodesy and Geodynamics, 2016, 36(5): 404-407)
( 0) |
| [7] |
Morishita Y, Kobayashi T, Yarai H. Three-Dimensional Deformation Mapping of Dike Intrusion Event in Sakurajima in 2015 by Exploiting the Right and Left Looking ALOS-2 InSAR[J]. Geophysical Research Letters, 2016, 43(9): 4 197-4 204 DOI:10.1002/2016GL068293
( 0) |
| [8] |
刘媛媛, 张勤, 赵超英, 等. PS-InSAR技术用于太原市地面沉降形变分析[J]. 大地测量与地球动力学, 2014, 34(2): 6-9 (Liu Yuanyuan, Zhang Qin, Zhao Chaoying, et al. Analysis of Ground Subsidence Deformation in Taiyuan City with PS-InSAR Technique[J]. Journal of Geodesy and Geodynamics, 2014, 34(2): 6-9)
( 0) |
| [9] |
Yu J S, Zhao B, Xu W B, et al. Oblique Fault Movement during the 2016 MW5.9 Zaduo Earthquake: Insights into Regional Tectonics of the Qiangtang Block, Tibetan Plateau[J]. Journal of Seismology, 2020, 24(3): 693-708 DOI:10.1007/s10950-020-09930-7
( 0) |
| [10] |
Bagnardi M, Hooper A. Inversion of Surface Deformation Data for Rapid Estimates of Source Parameters and Uncertainties: A Bayesian Approach[J]. Geochemistry, Geophysics, Geosystems, 2018, 19(7): 2 194-2 211 DOI:10.1029/2018GC007585
( 0) |
| [11] |
Iguchi M, Tameguri T, Ohta Y, et al. Characteristics of Volcanic Activity at Sakurajima Volcano's Showa Crater during the Period 2006 to 2011[J]. Bulletin of the Volcanological Social of Japan, 2013, 58(1): 115-135
( 0) |
| [12] |
Hotta K, Iguchi M, Ohkura T, et al. Multiple-Pressure-Source Model for Ground Inflation during the Period of High Explosivity at Sakurajima Volcano, Japan - Combination Analysis of Continuous GNSS, Tilt and Strain Data[J]. Journal of Volcanology and Geothermal Research, 2016, 310: 12-25 DOI:10.1016/j.jvolgeores.2015.11.017
( 0) |
| [13] |
Iguchi M. Magma Movement from the Deep to Shallow Sakurajima Volcano as Revealed by Geophysical Observations[J]. Bulletin of the Volcanological Society of Japan, 2013, 58: 1-18
( 0) |
2. Key Laboratory of Mine Environmental Monitoring and Improving around Poyang Lake, MNR, 418 Guanglan Road, Nanchang 330013, China