利用SBAS-InSAR技術(shù)研究印尼Agung火山的形變特征與巖漿房參數(shù)*

季靈運(yùn)，ZHONG Lu，王慶良，劉瑞春，秦姍蘭

0 引言

Agung火山位于印度尼西亞的巴厘島東部，是該島的最高峰，最高處海拔3 142 m，構(gòu)造上位于爪哇海溝以北約300 km處，屬于印度—澳大利亞板塊向東南亞大陸下方俯沖的弧后火山。Agung火山體呈錐狀，頂部存有520 m×375 m的漏斗狀破火山口，無(wú)植被生長(zhǎng)，火山體的側(cè)面陡峭，植被較茂盛 (Dilmy，1965)。歷史上，Agung火山曾多次噴發(fā)，在過(guò)去的300年里，有記錄的噴發(fā)事件分別發(fā)生在1808、1821(待考證)、1843、1963～1964年 (Seach，2012)。其中，1963～1964年的布里尼式噴發(fā)奪去了至少1 000人的生命，爆炸性噴發(fā)所產(chǎn)生的巨大的火山灰云隨后籠罩整個(gè)地球，造成了巨大的災(zāi)難 (Frederic，1970)。噴發(fā)歷史表明，Agung火山屬于具有巨大噴發(fā)危險(xiǎn)的當(dāng)代活動(dòng)火山，一旦噴發(fā)，將會(huì)給周圍居民乃至全球氣候帶來(lái)巨大災(zāi)難。然而，由于當(dāng)?shù)厣鐣?huì)經(jīng)濟(jì)條件落后等方面的原因，Agung火山目前缺乏地表形變監(jiān)測(cè)資料。

近年發(fā)展起來(lái)的InSAR技術(shù)因其自身的優(yōu)勢(shì)已在火山區(qū)地表形變監(jiān)測(cè)方面得到了廣泛應(yīng)用(Massonnet et al，1995;Hooper et al，2004，2007;Peltier et al，2010)。常規(guī)InSAR技術(shù)容易受到時(shí)間、空間失相干的影響，Berardino等 (2002)提出了小基線算法 (small baseline subsets，簡(jiǎn)稱SBAS)，該算法僅選取短基線干涉對(duì)進(jìn)行時(shí)間序列形變求解，削弱了空間失相干的影響，同時(shí)大大降低了地形因素引入的誤差和大氣延遲相位。Jung等 (2008)對(duì)小基線算法在4個(gè)方面進(jìn)行了精化，提高了時(shí)間序列形變結(jié)果的可靠性:(1)選取高質(zhì)量的干涉圖 (不含相位解纏誤差，且相干性好)估計(jì)外部DEM誤差和線性形變相位;(2)對(duì)形變時(shí)間序列和誤差進(jìn)行迭代計(jì)算;(3)利用有限差分平滑方法削弱時(shí)間域的噪聲;(4)對(duì)參考點(diǎn)本身存在的誤差進(jìn)行了校正。目前小基線算法已在火山區(qū)緩慢地表形變監(jiān)測(cè)方面表現(xiàn)出了極大的潛力和優(yōu)勢(shì)(Lundgren et al，2004;Pepe et al，2008;Casu et al，2008;季靈運(yùn)等，2011)。

1 InSAR數(shù)據(jù)處理與結(jié)果

Agung火山位于熱帶地區(qū)，植被覆蓋茂盛 (圖1)，本文選用了對(duì)植被具有較強(qiáng)穿透能力的日本ALOS PALSAR影像 (波長(zhǎng)23.6 cm)，獲取時(shí)間從2007年7月至2009年1月，選取的影像提供了1.5 a的地表形變觀測(cè)資料。

1.1 InSAR數(shù)據(jù)處理

常規(guī)的D-InSAR技術(shù)容易受到大氣延遲相位的干擾。Zebker等 (1997)的研究表明，云的相對(duì)濕度在空間或時(shí)間域內(nèi)變化20%將會(huì)導(dǎo)致10～14 cm的形變誤差。對(duì)于Agung火山區(qū)，大氣相位延遲在某些干涉圖中表現(xiàn)的比較明顯。圖2列出了3幅干涉圖 (已解纏)，可以看出圖2a和圖2b中的顏色分布相反，若差分干涉相位不包含誤差，則兩者反映了相反的形變態(tài)勢(shì)，即在圖2a干涉對(duì)的時(shí)間間隔內(nèi)，火山區(qū)域發(fā)生了明顯的下沉，最大下沉達(dá)7.5 cm;而在圖2b干涉對(duì)的時(shí)間間隔內(nèi)，火山區(qū)域則發(fā)生了明顯的隆升，最大隆升達(dá)3.8 cm;圖2c卻未反映出較大的相位差信息。2007年7～10月火山并無(wú)噴發(fā)，如此短的時(shí)間內(nèi) (92 d)形變態(tài)勢(shì)發(fā)生轉(zhuǎn)換，不符合火山區(qū)形變的物理解釋。圖2a，b中的相位分布與地形具有較強(qiáng)的相關(guān)性，即靠近火山體頂部，相位的絕對(duì)值較大，頂部和側(cè)面的相位差達(dá)到3～4弧度，符合對(duì)流層水汽垂直分層效應(yīng)對(duì)雷達(dá)信號(hào)延遲干擾的特征。由此可以得出，圖2a，b中的較大相位差均由2007年8月24日這景影像攜帶的大氣延遲相位所導(dǎo)致。

圖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

考慮到大氣延遲相位的干擾，本文選擇小基線算法求解時(shí)間序列形變。給定垂直基線小于1 000 m，時(shí)間間隔小于6個(gè)月，自由組合干涉生成33幅干涉圖。對(duì)于每個(gè)干涉圖，多視因子給定4∶9，即每個(gè)像元約為 (30×30)m2。對(duì)含有基線誤差的干涉圖，基于已有DEM，利用最小二乘方法進(jìn)行了基線精化 (Rosen et al，1996;Lu，2007)。選取12幅高信噪比的干涉圖，組成了一個(gè)子集 (圖3)，進(jìn)行小基線算法解算。

1.2 基于SBAS-InSAR技術(shù)的Agung火山區(qū)形變時(shí)間序列

利用前述精化的SBAS-InSAR技術(shù)提取了Agung火山的形變時(shí)間序列形變場(chǎng) (圖4)。從整體上看，從2007年7月9日到2009年1月11日，Agung火山地區(qū)表現(xiàn)出了較強(qiáng)的地表隆升趨勢(shì)，與時(shí)間呈明顯的線性正相關(guān)，且形變從火山體側(cè)部到頂部逐漸增大，552 d內(nèi)累積最大形變超過(guò)12 cm。為了表現(xiàn)火山不同部分的時(shí)間序列形變特征，圖5給出了火山頂部 (CO)、火山體側(cè)部 (FL)以及參考點(diǎn)附近地區(qū) (NR)的形變時(shí)間序列曲線。從圖中可以看出，火山頂部和火山體側(cè)部表現(xiàn)出了相似的隆升形變特征，其中火山頂部的累積形變約為13 cm，火山體側(cè)部的累積形變約為5 cm。

圖3 InSAR干涉組合和基線圖Fig.3 InSAR interference combinations and baselines

2 時(shí)間序列累積形變場(chǎng)模擬與分析

圖5 火山體不同部位的形變時(shí)間序列 (火山體不同部位見圖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技術(shù)的LOS方向時(shí)間序列形變場(chǎng) (2007年7月至2009年1月)(無(wú)顏色的地區(qū)表示雷達(dá)信號(hào)失相干)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)

火山區(qū)的地表形變一般由于其下方巖漿房或者巖漿系統(tǒng)的擾動(dòng)等所引起的 (Dvorak，Dzurisin，1997)。為了解釋2007～2009年Agung火山區(qū)的隆升形變態(tài)勢(shì)，研究其下方可能存在的巖漿房或巖漿系統(tǒng)的動(dòng)力學(xué)機(jī)制，本文選取552 d的累積形變場(chǎng)，分別利用埋置于均勻彈性半空間中的點(diǎn)源Mogi模型(Mogi，1958)和豎直橢球體模型 (Yang et al，1988)進(jìn)行了模擬。考慮到可能殘留的軌道誤差，模擬過(guò)程中引入了兩個(gè)方向的線性趨勢(shì)項(xiàng) (Massonnet，F(xiàn)eigl，1998)。采用下山單純形算法和蒙特卡洛搜索算法估計(jì)最優(yōu)參數(shù)及其誤差 (Press et al，1992)，以殘差干涉圖 (原始干涉圖與模擬干涉圖之差)的均方根誤差作為迭代結(jié)束的閾值。

圖6列出了模擬結(jié)果，其中圖6a為552 d的累積形變場(chǎng)，為了削弱噪聲，將觀測(cè)形變場(chǎng) (圖4h)進(jìn)行了均值濾波，圖6b和d分別為Mogi點(diǎn)源模型和豎直橢球體模型模擬的結(jié)果，圖6c和e為相應(yīng)的殘差圖。從整體上看，兩模型均可以很好地模擬原始觀測(cè)形變場(chǎng)。表1列出了兩模型模擬的Agung火山區(qū)巖漿房的幾何參數(shù)，兩模型模擬的巖漿房的中心位置較一致，均位于火山口的下方，但是豎直的橢球體模型模擬的結(jié)果相對(duì)較深。為了定量直觀地比較兩模型的模擬結(jié)果，沿火山錐體做了一條剖面 (圖6a中黑線，圖6f為模擬結(jié)果)，從圖6f中可以看出，Mogi點(diǎn)源模型能夠很好地?cái)M合火山體側(cè)面的形變場(chǎng)，而對(duì)火山體頂部的形變場(chǎng)擬合較差，而豎直的橢球體模型對(duì)火山體各個(gè)部位的形變場(chǎng)均能很好地模擬。由此表明，豎直的橢球體模型更加接近于真實(shí)的巖漿房幾何形狀。

圖6 基于Mogi、豎直橢球體模型的形變場(chǎng)模擬結(jié)果(a)時(shí)間序列累積形變場(chǎng) (圖4h)的均值濾波圖 (黑線表示剖面位置);(b)Mogi模型模擬的形變場(chǎng);(c)Mogi模型模擬的殘差圖;(d)豎直橢球體模型模擬的形變場(chǎng);(e)豎直橢球體模型模擬的殘差圖;(f)模擬結(jié)果(黑線表示地形，紅線表示觀測(cè)值，綠線和藍(lán)線分別表示Mogi、豎直橢球體模型模擬結(jié)果)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)

綜合火山區(qū)的形變時(shí)間序列與基于點(diǎn)源Mogi模型、豎直的橢球體模型模擬實(shí)驗(yàn)，筆者認(rèn)為Agung火山區(qū)在2007～2009年發(fā)生了比較明顯的隆升形變，巖漿房中心位于火口下方約5 km處(表1列出了豎直的橢球體模型得到的深度約為平均海平面以下2 km，火山體最高處約為3 km)，具備再次噴發(fā)的危險(xiǎn)，需加強(qiáng)跟蹤監(jiān)測(cè)。

表1 Agung火山巖漿房幾何參數(shù) (誤差為95%置信區(qū)間)Tab.1 Magma chamber geometric parameters of Agung volcano(the error is at 95%confidence interval)

3 結(jié)論

本文基于ALOS PALSAR影像，利用SBAS-In-SAR技術(shù)，提取了位于印度尼西亞巴厘島的Agung火山2007～2009年的地表形變時(shí)間序列，彌補(bǔ)了該火山區(qū)一直以來(lái)無(wú)形變監(jiān)測(cè)結(jié)果的現(xiàn)狀。結(jié)果表明，Agung火山在2007～2009年發(fā)生了較明顯的隆升形變，且與時(shí)間呈正相關(guān)，具有再次噴發(fā)的危險(xiǎn)?；贛ogi點(diǎn)源模型以及豎直橢球體模型，利用552 d的累積形變場(chǎng)模擬得到了巖漿房的幾何參數(shù)，結(jié)果表明，Agung火山的巖漿房位于火山體下方約5 km處。

Agung火山的SBAS-InSAR技術(shù)試驗(yàn)再次表明了InSAR技術(shù)在提取區(qū)域地表緩慢形變方面的巨大潛力和優(yōu)勢(shì)。InSAR技術(shù)的應(yīng)用，為Agung火山區(qū)的現(xiàn)今巖漿活動(dòng)狀況和地表形變態(tài)勢(shì)分析提供了重要依據(jù)。我國(guó)境內(nèi)分布有多座活動(dòng)火山，受氣候、當(dāng)?shù)厣鐣?huì)經(jīng)濟(jì)等條件的限制，還有一部分處于無(wú)監(jiān)測(cè)或少監(jiān)測(cè)資料的現(xiàn)狀。本文基于SBAS－InSAR技術(shù)對(duì)Agung火山開展的研究工作為查明和研究我國(guó)火山的現(xiàn)今活動(dòng)特征起到重要的借鑒作用。

季靈運(yùn)，王慶良，崔篤信，等.2011.利用SBAS-DInSAR技術(shù)提取騰沖火山區(qū)形變時(shí)間序列［J］．大地測(cè)量與地球動(dòng)力學(xué)，31(4):149－154.

Berardino P，F(xiàn)ornaro 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，F(xiàn)eigl 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.