季靈運，ZHONG Lu，王慶良，劉瑞春，秦姍蘭

0 引言

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

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

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

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

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

常規(guī)的D-InSAR技術容易受到大氣延遲相位的干擾。Zebker等 (1997)的研究表明，云的相對濕度在空間或時間域內(nèi)變化20%將會導致10～14 cm的形變誤差。對于Agung火山區(qū)，大氣相位延遲在某些干涉圖中表現(xiàn)的比較明顯。圖2列出了3幅干涉圖 (已解纏)，可以看出圖2a和圖2b中的顏色分布相反，若差分干涉相位不包含誤差，則兩者反映了相反的形變態(tài)勢，即在圖2a干涉對的時間間隔內(nèi)，火山區(qū)域發(fā)生了明顯的下沉，最大下沉達7.5 cm;而在圖2b干涉對的時間間隔內(nèi)，火山區(qū)域則發(fā)生了明顯的隆升，最大隆升達3.8 cm;圖2c卻未反映出較大的相位差信息。2007年7～10月火山并無噴發(fā)，如此短的時間內(nèi) (92 d)形變態(tài)勢發(fā)生轉換，不符合火山區(qū)形變的物理解釋。圖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火山區(qū)形變時間序列

利用前述精化的SBAS-InSAR技術提取了Agung火山的形變時間序列形變場 (圖4)。從整體上看，從2007年7月9日到2009年1月11日，Agung火山地區(qū)表現(xiàn)出了較強的地表隆升趨勢，與時間呈明顯的線性正相關，且形變從火山體側部到頂部逐漸增大，552 d內(nèi)累積最大形變超過12 cm。為了表現(xiàn)火山不同部分的時間序列形變特征，圖5給出了火山頂部 (CO)、火山體側部 (FL)以及參考點附近地區(qū) (NR)的形變時間序列曲線。從圖中可以看出，火山頂部和火山體側部表現(xiàn)出了相似的隆升形變特征，其中火山頂部的累積形變約為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月)(無顏色的地區(qū)表示雷達信號失相干)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)的擾動等所引起的 (Dvorak，Dzurisin，1997)。為了解釋2007～2009年Agung火山區(qū)的隆升形變態(tài)勢，研究其下方可能存在的巖漿房或巖漿系統(tǒng)的動力學機制，本文選取552 d的累積形變場，分別利用埋置于均勻彈性半空間中的點源Mogi模型(Mogi，1958)和豎直橢球體模型 (Yang et al，1988)進行了模擬。考慮到可能殘留的軌道誤差，模擬過程中引入了兩個方向的線性趨勢項 (Massonnet，F(xiàn)eigl，1998)。采用下山單純形算法和蒙特卡洛搜索算法估計最優(yōu)參數(shù)及其誤差 (Press et al，1992)，以殘差干涉圖 (原始干涉圖與模擬干涉圖之差)的均方根誤差作為迭代結束的閾值。

圖6列出了模擬結果，其中圖6a為552 d的累積形變場，為了削弱噪聲，將觀測形變場 (圖4h)進行了均值濾波，圖6b和d分別為Mogi點源模型和豎直橢球體模型模擬的結果，圖6c和e為相應的殘差圖。從整體上看，兩模型均可以很好地模擬原始觀測形變場。表1列出了兩模型模擬的Agung火山區(qū)巖漿房的幾何參數(shù)，兩模型模擬的巖漿房的中心位置較一致，均位于火山口的下方，但是豎直的橢球體模型模擬的結果相對較深。為了定量直觀地比較兩模型的模擬結果，沿火山錐體做了一條剖面 (圖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)

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

表1 Agung火山巖漿房幾何參數(shù) (誤差為95%置信區(qū)間)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年的地表形變時間序列，彌補了該火山區(qū)一直以來無形變監(jiān)測結果的現(xiàn)狀。結果表明，Agung火山在2007～2009年發(fā)生了較明顯的隆升形變，且與時間呈正相關，具有再次噴發(fā)的危險?；贛ogi點源模型以及豎直橢球體模型，利用552 d的累積形變場模擬得到了巖漿房的幾何參數(shù)，結果表明，Agung火山的巖漿房位于火山體下方約5 km處。

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

季靈運，王慶良，崔篤信，等.2011.利用SBAS-DInSAR技術提取騰沖火山區(qū)形變時間序列［J］．大地測量與地球動力學，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.