利用 SBAS-D InSAR技術提取騰沖火山區形變時間序列＊

季靈運 王慶良 崔篤信 胡亞軒 郝 明 李煜航 秦姍蘭

(中國地震局第二監測中心,西安 710054)

利用 SBAS-D InSAR技術提取騰沖火山區形變時間序列＊

季靈運 王慶良 崔篤信 胡亞軒 郝 明 李煜航 秦姍蘭

(中國地震局第二監測中心,西安 710054)

基于 6景 JERS-1 L波段 SAR影像,利用小基線集-合成孔徑雷達干涉測量技術,通過線性形變相位、非線性形變相位、大氣延遲相位以及地形殘差相位的分離,提取了騰沖火山地區 1995—1997年間地表形變時間序列(雷達視線向),與 2003—2004年的 GPS觀測結果對比表明,SBAS-D InSAR技術提取地殼形變的精度可達亞厘米級。時間序列形變顯示膽札-高田斷裂兩側形變差異性顯著,可能與其下方存在的地殼巖漿囊的活動有關。打鷹山地區地表形變揭示其下方可能存在隱伏斷裂。

SBAS-D InSAR;JERS-1;形變時間序列;大氣延遲相位;騰沖火山區

1 引言

騰沖火山區南北長約 90 km、東西寬約 50 km,共有 68座新生代火山[1],其中全新世有過活動的有3座(黑空山、馬鞍山、打鷹山)[2]。前人研究結果表明,騰沖火山區現今仍有較強的活動性,具有再次噴發的危險[3-10]。對騰沖火山的形變監測,云南省地震局分別于 1997年和 2003年共埋設了 94個水準點,并進行了多期觀測,并于 2002、2003、2004年進行了 GPS組網觀測,李成波等[6]基于 GPS水平形變資料計算了地殼面膨脹率和面應變,認為壓力源位于騰沖西南方向。李春光等[7]計算了 1998—1999年的垂直形變,發現打鷹山附近形成了一個明顯的隆起區,而其南面則呈下降趨勢。胡亞軒[8,9]、施行覺[10]等對 1998、1999、2000 3期水準數據進行了分析,結果顯示該區巖漿活動比較活躍,且壓力源中心位置不定,在破火山口周圍遷移,他們還利用Mogi模型對巖漿源參數進行了反演,結果表明火山活動受斷層控制。由前人的形變研究成果不難得出,騰沖火山區近年來存在形變異常。所以,跟蹤騰沖火山區的形變演化,對了解巖漿源的活動特征具有重要的意義。然而,常規的 GPS、水準測量由于存在作業周期長、僅能獲取點的形變信息等缺陷,很難全局把握火山的形變活動規律,近年發展起來的 InSAR技術可以實現全天候、高分辨率、大空間尺度的連續觀測,D InSAR技術已在地表形變監測方面得到了廣泛應用[11-13],SBAS-D InSAR技術[14]克服了常規D InSAR監測中的時空失相干因素影響,在城市緩慢地表形變[15]、斷層蠕動[16]、火山區域緩慢地表形變[17-20]監測等方面表現出了極大的潛力和優勢。所以本文嘗試利用 SBAS-D InSAR技術提取騰沖火山區形變時間序列,探討騰沖火山活動的狀態。

2 SBAS-D InSAR技術原理與數據處理流程

SBAS算法中,差分干涉相位定義為[14]：

其中 x和 r為像元坐標,λ為雷達波長,Δd為雷達視線方向地表形變,B⊥為垂直基線,θ為 SAR視角, Δz為地形殘差,Δφatm為大氣延遲相位,Δφn為其他噪聲相位。顯然,SBAS算法將差分干涉相位分為地表形變相位、地形殘差相位、大氣延遲相位以及其他噪聲相位 4部分。其數據處理流程如圖 1所示。

圖1 SBAS-D InSAR數據處理流程Fig.1 Flow chart of SBAS-D InSAR processing

文獻[14-17]的研究表明,SBAS算法的精度可達毫米級,與 GPS觀測結果也吻合得較好。SBAS算法由于在影像自由組合干涉時限制了時間基線和空間基線,保證了每幅干涉圖的高相干性,又利用奇異值分解的方法將多個干涉圖子集聯合進行最小二乘求解,增加了時間采樣。SBAS算法將外部DEM誤差產生的相位分離出來,降低了使用外部DEM所引入的誤差。充分考慮到大氣延遲相位對形變結果的影響,SBAS算法通過時空域濾波的方法削弱了大氣延遲相位。相對于永久散射體技術,小基線集技術獲取到的形變序列在空間上更為連續,從而可以應用于監測地殼長時間緩慢變形。

3 SBAS-D InSAR技術在騰沖火山區的應用

3.1 SAR影像的選取

騰沖火山區地貌復雜,所處緯度偏低,植被覆蓋茂密,很容易造成 SAR影像失相干,所以本文選用L波段 JERS-1 SAR數據,L波段數據能夠較好地克服植被覆蓋所造成的失相干問題。數據處理結果表明,在騰沖火山區,時間基線仍不超過一年才具有較好的相干性。6景 JERS-1 SAR影像,獲取時間分別為 1995年 4月 10日、1995年 8月 20日、1995年 11月 16日、1996年 3月 27日、1996年 8月 6日、1997年 4月 27日,給定時間基線為 365天,空間垂直基線為 2 km作為限制條件,進行自由干涉組合,得到6個可用的干涉對,時間跨越了 1995—1997年。可用干涉對的時間基線和空間基線關系如圖 2所示。

圖2 干涉對基線組合Fig.2 Baseline combination of all the interferomeric pairs

3.2 常規 D InSAR技術處理

對每個符合時間、空間基線限制條件的干涉對均利用兩軌法 D InSAR技術處理,處理平臺為瑞士GAMMA遙感公司開發的GAMMA軟件。處理步驟主要包括：原始 SAR信號成像、圖像配準與干涉、外部DEM模擬、基線估計與精化、相干系數計算、干涉圖自適應濾波、相位解纏。外部地形數據使用美國NASA(美國宇航局)公布的 SRT M DEM。另外,為了削弱噪聲,增加相干性,同時為了檢測相對大尺度的形變信息,對干涉圖進行了多視處理,多視因子為9：18。相位解纏采用最小費用網絡流算法,為了與GPS觀測資料進行比對,采用騰沖縣城附近為解纏起始區域。

3.3 SBAS技術處理

對所有符合條件的 6個差分干涉對的解纏相位進行小基線集技術處理。首先,通過差分干涉相位估計得到研究區域的線性形變相位,并轉換為 LOS向形變速率,其中 3座全新世火山的線性形變速率分別為黑空山 -2.8 mm/a,打鷹山 -6.6 mm/a,馬鞍山 -0.82 mm/a;同時得到了殘差地形相位,發現其值均為 -0.1～0.1 rad,換算為 LOS向的形變約為 -1.4～1.4 mm(對于 JERS-1)。然后從原始差分干涉相位中減去線性形變相位和殘差地形相位,再對殘余相位進行解纏,利用奇異值分解方法求解非線性形變相位,最后聯合線性形變相位和非線性形變相位得到時間序列累積形變相位。對于大氣延遲相位的估計,首先在時間域上進行均值濾波處理,然后在空間域上采用小波分解的方法,分離出高頻成分,從而實現了大氣延遲相位的分離。圖 3是1995年 8月 20日的大氣延遲相位圖,從圖 3可以看出,最大達 1.5 rad。

圖3 大氣延遲相位(1995-08-20)Fig.3 I mage on the atmospheric phase delay(1995-08-20)

3.4 SBAS形變結果的質量評價

圖4為LOS向形變時間序列。根據雷達成像的幾何關系,地面的上、東、北 3方向形變對 LOS向形變的貢獻可表示為：

其中,θ為雷達脈沖入射角,φ為衛星軌道方位角。

本文選用的 GPS站點的形變速率值基本含蓋了 InSAR數據所覆蓋的范圍。將 GPS觀測的 3方向速率分量帶入式 (2),得到其在 LOS向的形變速率,與 InSAR形變速率值的比較如表 1所示。

表1 GPS與 I nSAR獲得的LOS向形變速率值的比較Tab.1 Comparison between deformation velocities(LOS) measured by GPS and I nSAR

由表1可以看出,雖然 GPS觀測時間與 SAR成像時間不同,得到的速率不同,但是兩者的一致性很好,最大差值約 1 cm/a,差值的標準差為 0.65 cm/ a。對比結果表明 SBAS技術對地表形變的監測精度可達亞厘米級。

3.5 騰沖火山區形變時間序列分析

圖4直觀地表現了騰沖火山區 1995—1997年的地表形變 (LOS向)隨時間的演化過程。從全局來看,騰沖火山區發生了較大范圍的區域性地表變形,大部分地區以下降為主要形變特征,其中下降相對較大的地方出現在黑空山西南以及騰沖縣城以南地區,累積最大約 8 cm;研究區的西南部在 1995—1996年發生了隆升,最大達 4 cm,1996年后又回落;從局部來看,個別斷裂在形變場上有明顯顯示,反映出斷裂兩側有明顯的差異活動,例如西南部的北西向膽札-高田斷裂以及中部的南北向固東-騰沖斷裂。圖 5給出了垂直于膽札-高田斷裂走向的時間序列形變剖面曲線 (圖 4A-A’剖面),可以看出,斷裂東北側遠離斷層逐漸下沉,而西南側遠離斷層則表現為隆升,斷裂兩側存在明顯的形變梯度,是斷裂活動的反映。這一結果與李春光等[7]計算得到的 1998—1999年形變特征一致,與野外近場數字化臺站記錄結果所顯示的火山地震活動主要發生在打鷹山-馬鞍山一帶以及騰沖縣城西南部地區[1]也是一致的,表明該斷裂現今仍在活動。值得注意的是,膽札-高田斷裂的西南側在 1996年后發生了約 4 cm的不規則形變變化,體現了該斷裂的較強活動性[4]。

從圖 4還可以定性地得出：1995—1997年,相對于騰沖縣城,打鷹山的東北和西南兩側表現為明顯的差異性形變特征,馬鞍山地區形變不明顯,黑空山一直處于下降狀態。圖 6給出了穿過打鷹山的形變時間序列剖面曲線(圖 4B-B’剖面),從圖 6可以看出,打鷹山的東北側表現為上升,而西南側表現為下降,累積最大差異達 10 cm,這種巨大的差異性運動特征預示著打鷹山下方可能存在北西向隱伏斷裂,且活動性較強。這一結果與 1998—2000年的水準測量結果[8]一致,1998—1999年與 1999—2000年打鷹山的南北兩側的垂直形變差異性非常明顯,推測壓力源的位置不是固定不變的,沿破火口周圍遷移。表 2列出了馬鞍山和黑空山兩座火山的時間序列形變量值,可以看出,兩地區均表現為緩慢下降,但馬鞍山地區下降幅度相對較小,累積不到 1 cm,表明其在 1995—1997年活動性不顯著。黑空山附近表現為連續性下沉,累積最大超過 3 cm,表明其在 1995—1997年活動性不強。

圖4 地表形變(LOS向)時間序列(1995—1997年)Fig.4 Ti me seriers of defor mation on LOS(1995-1997)

圖5 膽札-高田斷裂形變時間序列剖面Fig.5 Cross section on time series of defor mation of Danzha-Gaotian fault

圖6 打鷹山形變時間序列剖面Fig.6 Cross section on time series of deformation of Dayingshan volcano area

表2 兩座全新世活動火山的形變時間序列Tab.2 Ti me series of deformation of two Holocene volcanoes

綜合火山區的整體形變時間序列和 3座全新世活動火山的形變時間序列演化特征分析,可以得出：騰沖火山區現今地殼運動比較活躍,形變特征復雜多樣;打鷹山地區活動性比較顯著,下方可能存在隱伏斷裂;膽札-高田斷裂兩側形變相對較大,差異性顯著,地球物理資料揭示其下方存在活動性較強的巖漿囊可能是導致其差異性形變的主要因素。

4 結論

1)相對于常規 D InSAR技術,SBAS-D InSAR技術分離了大氣延遲相位以及地形殘差相位等干擾因素,從而能夠獲得高精度的地表形變時間序列 (LOS向),騰沖火山區的試驗結果表明 SBAS-D InSAR技術能夠以亞厘米級的精度檢測區域地表形變。

2)膽札-高田斷裂兩側形變差異性顯著,可能與其下方存在的地殼巖漿囊的活動有關。

3)打鷹山地區地表形變揭示其下方可能存在隱伏斷裂。

致謝 感謝路中博士提供 JERS-1 SAR影像!

1 皇甫崗,姜朝松.騰沖火山研究[M].昆明：云南科技出版社,2000.(Huangfu Gang and Jiang Chaosong.Study on Tengchong volcanic activity[M].Kunming：Yunnan Press of Science and Technology,2000)

2 李大明,李齊,陳文寄.騰沖火山區上新世以來的火山活動[J].巖石學報,2000,16(3)：362-370.(Li Daming, LiQi and ChenWenji.Volcanic activities in the Tengchong volcano area since Pliocene[J].Acta Petrologica Sinica, 2000,16(3)：362-370)

3 葉建慶,等.騰沖火山區地震群的活動特征[J].地震地質,2003,25(增)：128-137. (Ye Jian qing,et al. Characteristics of earthquake cluster activity in Tengchong volcanic area[J].Seismology and Geology,2003,25(supp l.)：128-137)

4 趙慈平,冉華,陳坤華.由相對地熱梯度推斷的騰沖火山區現存巖漿囊 [J].巖石學報,2006,22(6)：1 517-1 528.(Zhao Ciping,Ran Hua and Chen Kunhua.Presentday magma chambers in Tengchong volcano area inferred from relative geother mal gradient[J].Acta Petrologica Sinica,2006,22(6)：1 517-1 528)

5 樓海,等.云南騰沖火山區上部地殼三維地震速度層析成像[J].地震學報,2002,24(3)：243-251.(Lou Hai, et al.Three dimensional seis mic velocity tomography of the upper crust in Tengchong volcanic area,Yunnan province [J].Acta Seismologica Sinica,2002,24(3)：243-251)

6 李成波,等.騰沖火山區的 GPS形變特征[J].地球物理學進展,2007,22(3)：765-770.(Li Chengbo,et al.Research on the character of the GPS defor mation in the Tengchong volcano area[J].Progress in Geophysics,2007,22 (3)：765-770)

7 李春光,等.騰沖火山區的形變特征[J].地震研究,2000, 23(2)：165-171.(Li Chunguang,et al.Deformation feature in Tengchong volcano areas[J].Journal of Seis mological Research,2000,23(2)：165-171)

8 胡亞軒,等.騰沖火山區地表垂直形變分析[J].大地測量與地球動力學,2003,(2)：37-41.(Hu Yaxuan,et al. Analysis of vertical deformation in Tengchong volcanic area [J].Journal of Geodesy and Geodynamics,2003,(2)：37 -41)

9 胡亞軒,王慶良,趙慈平.騰沖火山區形變分析[J].國際地震動態,2008,352(4)：42-47.(Hu Yaxuan,Wang Qingliang and Zhao Ciping.Analysis on the deformation in Tengchong volcanic region[J]. Recent Developments in World Seismology,2008,352(4)：42-47)

10 施行覺,等.以垂直形變資料反演騰沖火山區巖漿活動性的初步研究 [J].地震研究,2005,28(3)：256-261. (Shi Xingjue,et al.Vertical deformation inversion ofmagma activity in Tengchong volcanic area[J].Journalof Seismological Research,2005,28(3)：256-261)

11 路旭,等.用 I NSAR作地面沉降監測的試驗研究[J].大地測量與地球動力學,2002,(4)：66-70.(Lu Xu,et al.Experiment on land subsidence monitoring by InSAR [J].Journal of Geodesy and Geodynamics,2002,22(4)：66-70)

12 喬學軍,郭利民.新疆伽師強震區的 InSAR觀測研究[J].大地測量與地球動力學,2007,(1)：7-13.(Qiao Xuejun and Guo Limin.Study on Jiashi strong earthquake s warm area by InSAR[J].Journal of Geodesy and Geodynamics,2007,(1)：7-13)

13 喬學軍,等.當雄Ms6.6地震的 InSAR觀測及斷層位錯反演[J].大地測量與地球動力學,2009,(6)：1-7. (Qiao Xuejun,et al. Study on dislocation inversion of Ms6.6 Dangxiong earthquake as constrained by InSAR measurement[J].Journal of Geodesy and Geodynamics, 2009,(6)：1-7)

14 Berardino P,et al.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J].IEEE Transact. Geoscience and Remote Sensing,2002,40(11)：2 375-2 383.

15 Antonio Pepe,et al.On the generation of ERS/ENV ISAT D InSAR time-series via the SBAS technique[J].IEEE Geoscience and Remote sensing letters,2005,2(3)：265-269.

16 Lanaria R,et al.Application of the SBAS-D InSAR technique to fault creep：A case study of the Hayward fault, California[J].Remote Sensing of Environment,2007,109 (1)：20-28.

17 Pepe A,et al.Surface deformation of active volcanic areas retrieved with the SBAS-D InSAR technique：an overview [J].Annals of Geophysics,2008,51(1)：247-263.

18 Chang-Wook Lee,et al.Surface defor mation of Augustine Volcano,1992-2005,from multiple-interferogram processing using a refined SmallBaseline Subset(SBAS)Interferometric Synthetic Aperture Radar(InSAR)approach [A].The 2006 eruption of Augustine Volcano,Alaska：U.S.Geological Survey Professional Paper 1769[C].453 -465.

19 Lundgren P,et al.Gravity and magma induced spreading ofMount Etna volcano revealed by satellite radar interferometry[J].Geophysical Research Letters,2004,31(4)：L04602.

20 Casu F,et al.Surface defor mation dynamics ofMauna Loa and Kilauea volcanoes,Hawaii,revealed by InSAR time series analysis[R].American Geophysical Union,Fall Meeting,2008,V11B-2023.

TI M E SERIES OF DEFORMATI ON IN TENGCHONG VOLCANIC AREA EXTRACTED BY SBAS-D InSAR

JiLingyun,WangQingliang,CuiDuxin,Hu Yaxuan,HaoMing,Li Yuhang and Qin Shanlan
(Second CrustM onitoring and Application Center,CEA,X i’an 710054)

On the basisof the 6 JERS-1 SAR images,the time seriesof deformation in Tengchong volcanic area were obtained by SBAS-D InSAR technique through separating the linear deformation phase,the nonlinear deformation phase,the at mospheric phase screen and theDEM residual errorphase.Comparedwith the GPSmeasurement results,SBAS-D InSAR technique can detect the crustal defor mation at the precision of sub-cm level.The defor mation of the Danzha-Gaotian fault is obviously different bet ween the two sides,which may relate the chamber below. The defor mation pattern of the Dayingshan area indicates that there is probably an buried fault underground.

SBAS-D InSAR;JERS-1;time seriesof deformation;atmospheric phase delay;Tengchong volcanic area

1671-5942(2011)04-0149-06

2011-03-23

國家自然科學基金(40974062);地震行業科研重點專項(200908029)

季靈運,男,1982年生,博士生,主要研究方向為 InSAR技術在地震、火山形變中的應用.E-mail：dinsar010@163.com

P225.1

A