季灵运，ZHONG Lu，王庆良，刘瑞春，秦姗兰

0 引言

Agung火山位于印度尼西亚的巴厘岛东部，是该岛的最高峰，最高处海拔3 142 m，构造上位于爪哇海沟以北约300 km处，属于印度—澳大利亚板块向东南亚大陆下方俯冲的弧后火山。Agung火山体呈锥状，顶部存有520 m×375 m的漏斗状破火山口，无植被生长，火山体的侧面陡峭，植被较茂盛 (Dilmy，1965)。历史上，Agung火山曾多次喷发，在过去的300年里，有记录的喷发事件分别发生在1808、1821(待考证)、1843、1963～1964年 (Seach，2012)。其中，1963～1964年的布里尼式喷发夺去了至少1 000人的生命，爆炸性喷发所产生的巨大的火山灰云随后笼罩整个地球，造成了巨大的灾难 (Frederic，1970)。喷发历史表明，Agung火山属于具有巨大喷发危险的当代活动火山，一旦喷发，将会给周围居民乃至全球气候带来巨大灾难。然而，由于当地社会经济条件落后等方面的原因，Agung火山目前缺乏地表形变监测资料。

近年发展起来的InSAR技术因其自身的优势已在火山区地表形变监测方面得到了广泛应用(Massonnet et al，1995;Hooper et al，2004，2007;Peltier et al，2010)。常规InSAR技术容易受到时间、空间失相干的影响，Berardino等 (2002)提出了小基线算法 (small baseline subsets，简称SBAS)，该算法仅选取短基线干涉对进行时间序列形变求解，削弱了空间失相干的影响，同时大大降低了地形因素引入的误差和大气延迟相位。Jung等 (2008)对小基线算法在4个方面进行了精化，提高了时间序列形变结果的可靠性:(1)选取高质量的干涉图 (不含相位解缠误差，且相干性好)估计外部DEM误差和线性形变相位;(2)对形变时间序列和误差进行迭代计算;(3)利用有限差分平滑方法削弱时间域的噪声;(4)对参考点本身存在的误差进行了校正。目前小基线算法已在火山区缓慢地表形变监测方面表现出了极大的潜力和优势(Lundgren et al，2004;Pepe et al，2008;Casu et al，2008;季灵运等，2011)。

1 InSAR数据处理与结果

Agung火山位于热带地区，植被覆盖茂盛 (图1)，本文选用了对植被具有较强穿透能力的日本ALOS PALSAR影像 (波长23.6 cm)，获取时间从2007年7月至2009年1月，选取的影像提供了1.5 a的地表形变观测资料。

1.1 InSAR数据处理

常规的D-InSAR技术容易受到大气延迟相位的干扰。Zebker等 (1997)的研究表明，云的相对湿度在空间或时间域内变化20%将会导致10～14 cm的形变误差。对于Agung火山区，大气相位延迟在某些干涉图中表现的比较明显。图2列出了3幅干涉图 (已解缠)，可以看出图2a和图2b中的颜色分布相反，若差分干涉相位不包含误差，则两者反映了相反的形变态势，即在图2a干涉对的时间间隔内，火山区域发生了明显的下沉，最大下沉达7.5 cm;而在图2b干涉对的时间间隔内，火山区域则发生了明显的隆升，最大隆升达3.8 cm;图2c却未反映出较大的相位差信息。2007年7～10月火山并无喷发，如此短的时间内 (92 d)形变态势发生转换，不符合火山区形变的物理解释。图2a，b中的相位分布与地形具有较强的相关性，即靠近火山体顶部，相位的绝对值较大，顶部和侧面的相位差达到3～4弧度，符合对流层水汽垂直分层效应对雷达信号延迟干扰的特征。由此可以得出，图2a，b中的较大相位差均由2007年8月24日这景影像携带的大气延迟相位所导致。

图1 Agung火山周围SRTM DEM(右上角插图显示了巴厘岛及其周围岛屿的位置，虚线框表示干涉图的范围)Fig.1 Digital Elevation Map from Shuttle Radar Topography Mission surrounding Agung volcano(the illustration in the right top corner shows the location of Bali island(small rectangle)and its surrounding islands，the dashedline frame shows the range of interferogram)

图2 带有显著大气延迟相位的解缠干涉图(a)2007－07－09～2007－08－24;(b)2007－08－24～2007－10－09;(c)2007－07－09～2007－10－09Fig.2 Unwrapping interferograms with atmospheric delayed phase

考虑到大气延迟相位的干扰，本文选择小基线算法求解时间序列形变。给定垂直基线小于1 000 m，时间间隔小于6个月，自由组合干涉生成33幅干涉图。对于每个干涉图，多视因子给定4∶9，即每个像元约为 (30×30)m2。对含有基线误差的干涉图，基于已有DEM，利用最小二乘方法进行了基线精化 (Rosen et al，1996;Lu，2007)。选取12幅高信噪比的干涉图，组成了一个子集 (图3)，进行小基线算法解算。

1.2 基于SBAS-InSAR技术的Agung火山区形变时间序列

利用前述精化的SBAS-InSAR技术提取了Agung火山的形变时间序列形变场 (图4)。从整体上看，从2007年7月9日到2009年1月11日，Agung火山地区表现出了较强的地表隆升趋势，与时间呈明显的线性正相关，且形变从火山体侧部到顶部逐渐增大，552 d内累积最大形变超过12 cm。为了表现火山不同部分的时间序列形变特征，图5给出了火山顶部 (CO)、火山体侧部 (FL)以及参考点附近地区 (NR)的形变时间序列曲线。从图中可以看出，火山顶部和火山体侧部表现出了相似的隆升形变特征，其中火山顶部的累积形变约为13 cm，火山体侧部的累积形变约为5 cm。

图3 InSAR干涉组合和基线图Fig.3 InSAR interference combinations and baselines

2 时间序列累积形变场模拟与分析

图5 火山体不同部位的形变时间序列 (火山体不同部位见图4a中黑色三角)Fig.5 InSAR LOS deformation time series for different parts of volcano(the different parts of volcano marked by the black triangles in Fig.4a)

图4 基于SBAS－InSAR技术的LOS方向时间序列形变场 (2007年7月至2009年1月)(无颜色的地区表示雷达信号失相干)Fig.4 Time series deformation fields derived by SBAS－InSAR teohnology in LOS direction from Jul.，2007 to Jan.，2009(decorrelation of radar signal are uncolored)

火山区的地表形变一般由于其下方岩浆房或者岩浆系统的扰动等所引起的 (Dvorak，Dzurisin，1997)。为了解释2007～2009年Agung火山区的隆升形变态势，研究其下方可能存在的岩浆房或岩浆系统的动力学机制，本文选取552 d的累积形变场，分别利用埋置于均匀弹性半空间中的点源Mogi模型(Mogi，1958)和竖直椭球体模型 (Yang et al，1988)进行了模拟。考虑到可能残留的轨道误差，模拟过程中引入了两个方向的线性趋势项 (Massonnet，Feigl，1998)。采用下山单纯形算法和蒙特卡洛搜索算法估计最优参数及其误差 (Press et al，1992)，以残差干涉图 (原始干涉图与模拟干涉图之差)的均方根误差作为迭代结束的阈值。

图6列出了模拟结果，其中图6a为552 d的累积形变场，为了削弱噪声，将观测形变场 (图4h)进行了均值滤波，图6b和d分别为Mogi点源模型和竖直椭球体模型模拟的结果，图6c和e为相应的残差图。从整体上看，两模型均可以很好地模拟原始观测形变场。表1列出了两模型模拟的Agung火山区岩浆房的几何参数，两模型模拟的岩浆房的中心位置较一致，均位于火山口的下方，但是竖直的椭球体模型模拟的结果相对较深。为了定量直观地比较两模型的模拟结果，沿火山锥体做了一条剖面 (图6a中黑线，图6f为模拟结果)，从图6f中可以看出，Mogi点源模型能够很好地拟合火山体侧面的形变场，而对火山体顶部的形变场拟合较差，而竖直的椭球体模型对火山体各个部位的形变场均能很好地模拟。由此表明，竖直的椭球体模型更加接近于真实的岩浆房几何形状。

图6 基于Mogi、竖直椭球体模型的形变场模拟结果(a)时间序列累积形变场 (图4h)的均值滤波图 (黑线表示剖面位置);(b)Mogi模型模拟的形变场;(c)Mogi模型模拟的残差图;(d)竖直椭球体模型模拟的形变场;(e)竖直椭球体模型模拟的残差图;(f)模拟结果(黑线表示地形，红线表示观测值，绿线和蓝线分别表示Mogi、竖直椭球体模型模拟结果)Fig.6 Deformation field simulation results based on Mogi and vertical spheroid models(a)mean filter of time series accumulation deformation field in Fig.4(black line represents the position of the profile);(b)deformation field simulated by the Mogi model;(c)residual simulated by Mogi model;(d)deformation field simulated by the vertical spheroid models;(e)residual simulated by vertical spheroid models;(f)simulation results(the black line represents the terrain，the red line represents the observed value，the green and blue lines show the simulation results based on Mogi and vertical spheroid models respectively)

综合火山区的形变时间序列与基于点源Mogi模型、竖直的椭球体模型模拟实验，笔者认为Agung火山区在2007～2009年发生了比较明显的隆升形变，岩浆房中心位于火口下方约5 km处(表1列出了竖直的椭球体模型得到的深度约为平均海平面以下2 km，火山体最高处约为3 km)，具备再次喷发的危险，需加强跟踪监测。

表1 Agung火山岩浆房几何参数 (误差为95%置信区间)Tab.1 Magma chamber geometric parameters of Agung volcano(the error is at 95%confidence interval)

3 结论

本文基于ALOS PALSAR影像，利用SBAS-In-SAR技术，提取了位于印度尼西亚巴厘岛的Agung火山2007～2009年的地表形变时间序列，弥补了该火山区一直以来无形变监测结果的现状。结果表明，Agung火山在2007～2009年发生了较明显的隆升形变，且与时间呈正相关，具有再次喷发的危险。基于Mogi点源模型以及竖直椭球体模型，利用552 d的累积形变场模拟得到了岩浆房的几何参数，结果表明，Agung火山的岩浆房位于火山体下方约5 km处。

Agung火山的SBAS-InSAR技术试验再次表明了InSAR技术在提取区域地表缓慢形变方面的巨大潜力和优势。InSAR技术的应用，为Agung火山区的现今岩浆活动状况和地表形变态势分析提供了重要依据。我国境内分布有多座活动火山，受气候、当地社会经济等条件的限制，还有一部分处于无监测或少监测资料的现状。本文基于SBAS－InSAR技术对Agung火山开展的研究工作为查明和研究我国火山的现今活动特征起到重要的借鉴作用。

季灵运，王庆良，崔笃信，等.2011.利用SBAS-DInSAR技术提取腾冲火山区形变时间序列［J］．大地测量与地球动力学，31(4):149－154.

Berardino P，Fornaro G，Lanari R，et al.2002.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms［J］．Institute of Electrical and Electronics Engineers Transactions on Geosciences and Remote Sensing，40(11):2 375 －2 383.

Casu F，Poland M P，Solaro G，et al.2008.Surface deformation dynamics of Mauna Loa and Kilauea volcanoes，Hawaii，revealed by InSAR time series analysis［R］．Sanfrancisco:American Geophysical Union.

Dilmy A.1965.Pioneer plants found one year after the 1963 eruption of Agung in Bali［J］．Pac Sci，19(4):498 － 501.

Dvorak J J，Dzurisin D.1997.Volcano geodesy:the search for magma reservoirs and the formation of eruptive vents［J］．Rev Geophys，35(3):343－384.

Frederic E V.1970.Atmospheric turbidity after the Agung eruption of 1963 and size distribution of the volcanic aerosol［J］．JGR，75(27):5 185－5 194.

Hooper A，Segall P，Zebker H.2007.Persistent scatterer InSAR for crustal deformation analysis，with application to Volcano Alcedo，Galapagos［J］．JGR，112，B07407，doi:10.1029/2006JB004763，2007.

Hooper A，Zebker H，Segall P，et al.2004.A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers［J］．Geophys Res Lett，31(23)，L23611，doi:10.1029/2004GL021737，2004.

Jung H S，Lee C W，Park J W，et al.2008.Improvement of small baseline subset(SBAS)algorithm for measuring time-series surface deformations from differential SAR interferograms［J］．Korean Journal of Remote Sensing，24(2)，165 －177.

Lu Z.2007.InSAR imaging of volcanic deformation over cloud-prone areas－ Aleutian islands［J］．Photogrammetric Engineering ＆ Remote Sensing，73(3):245 －257.

Lundgren P，Casu F，Manzo M，et al.2004.Gravity and magma induced spreading of Mount Etna volcano revealed by satellite radar interferometry［J］．Geophys Res Lett，31(4):L04602.

Massonnet D，Briole P，Arnaud A.1995.Deflation of Mount Etna monitored by spaceborne radar interferometry［J］．Nature，375:567－570.

Massonnet D，Feigl K.1998.Radar interferometry and its application to changes in the Earth's surface［J］．Rev Geophys，36(4):441 － 500.

Mogi K.1958.Relations between the eruptions of various volcanoes and the deformations of the ground surface around them［J］．Bull Earthquake Res Inst Univ Tokyo，36:99 －134.

Peltier A，Bianchi M，Kaminski E，et al.2010.PSInSAR as a new tool to monitor preeruptive volcano ground deformation:Validation using GPS measurements on Piton de la Fournaise［J］．Geophys Res Lett，37(12)，L12301，doi:10.1029/2010GL043846.

Pepe A，Manzo M，Casu F，et al.2008.Surface deformation of active volcanic areas retrieved with the SBAS-DInSAR technique:an overview［J］．Annals of Geophysics，51(1):247 －263.

Press W，Teukolsky S，Vetterling W，et al.1992.Numerical Recipes in C，the Art of Scientific Computing［M］．New York:Cambridge Univ Press.

Rosen P，Hensley S，Zebker H，et al.1996.Surface deformation and coherence measurements of Kilauea Volcano，Hawaii，from SIR-C radar interferometry［J］．J Geophys Rea，101(E10):23 109 － 23 125.

Seach J.2012.(2000－01－01)［2012－11－01］．Agung Volcano-John Seach.Volcanism reference base［EB/OL］．http://www.volcanolive.com/agung.html.Last accessed 29th October 2012.

Yang X M，Davis P，Dieterich J H.1988.Deformation from inflation of a dipping finite prolate spheroid in an elastic half-space as a model for volcanic stressing［J］．JGR，93(B5):4 249 － 4257.

Zebker H A，Rosen P A，Hensley S.1997.Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps［J］．JGR，102(B4):7 547 －7 563.