季靈運(yùn) 王慶良 崔篤信 胡亞軒 郝 明 李煜航 秦姍蘭
(中國(guó)地震局第二監(jiān)測(cè)中心,西安 710054)
利用 SBAS-D InSAR技術(shù)提取騰沖火山區(qū)形變時(shí)間序列*
季靈運(yùn) 王慶良 崔篤信 胡亞軒 郝 明 李煜航 秦姍蘭
(中國(guó)地震局第二監(jiān)測(cè)中心,西安 710054)
基于 6景 JERS-1 L波段 SAR影像,利用小基線集-合成孔徑雷達(dá)干涉測(cè)量技術(shù),通過(guò)線性形變相位、非線性形變相位、大氣延遲相位以及地形殘差相位的分離,提取了騰沖火山地區(qū) 1995—1997年間地表形變時(shí)間序列(雷達(dá)視線向),與 2003—2004年的 GPS觀測(cè)結(jié)果對(duì)比表明,SBAS-D InSAR技術(shù)提取地殼形變的精度可達(dá)亞厘米級(jí)。時(shí)間序列形變顯示膽札-高田斷裂兩側(cè)形變差異性顯著,可能與其下方存在的地殼巖漿囊的活動(dòng)有關(guān)。打鷹山地區(qū)地表形變揭示其下方可能存在隱伏斷裂。
SBAS-D InSAR;JERS-1;形變時(shí)間序列;大氣延遲相位;騰沖火山區(qū)
騰沖火山區(qū)南北長(zhǎng)約 90 km、東西寬約 50 km,共有 68座新生代火山[1],其中全新世有過(guò)活動(dòng)的有3座(黑空山、馬鞍山、打鷹山)[2]。前人研究結(jié)果表明,騰沖火山區(qū)現(xiàn)今仍有較強(qiáng)的活動(dòng)性,具有再次噴發(fā)的危險(xiǎn)[3-10]。對(duì)騰沖火山的形變監(jiān)測(cè),云南省地震局分別于 1997年和 2003年共埋設(shè)了 94個(gè)水準(zhǔn)點(diǎn),并進(jìn)行了多期觀測(cè),并于 2002、2003、2004年進(jìn)行了 GPS組網(wǎng)觀測(cè),李成波等[6]基于 GPS水平形變資料計(jì)算了地殼面膨脹率和面應(yīng)變,認(rèn)為壓力源位于騰沖西南方向。李春光等[7]計(jì)算了 1998—1999年的垂直形變,發(fā)現(xiàn)打鷹山附近形成了一個(gè)明顯的隆起區(qū),而其南面則呈下降趨勢(shì)。胡亞軒[8,9]、施行覺(jué)[10]等對(duì) 1998、1999、2000 3期水準(zhǔn)數(shù)據(jù)進(jìn)行了分析,結(jié)果顯示該區(qū)巖漿活動(dòng)比較活躍,且壓力源中心位置不定,在破火山口周圍遷移,他們還利用Mogi模型對(duì)巖漿源參數(shù)進(jìn)行了反演,結(jié)果表明火山活動(dòng)受斷層控制。由前人的形變研究成果不難得出,騰沖火山區(qū)近年來(lái)存在形變異常。所以,跟蹤騰沖火山區(qū)的形變演化,對(duì)了解巖漿源的活動(dòng)特征具有重要的意義。然而,常規(guī)的 GPS、水準(zhǔn)測(cè)量由于存在作業(yè)周期長(zhǎng)、僅能獲取點(diǎn)的形變信息等缺陷,很難全局把握火山的形變活動(dòng)規(guī)律,近年發(fā)展起來(lái)的 InSAR技術(shù)可以實(shí)現(xiàn)全天候、高分辨率、大空間尺度的連續(xù)觀測(cè),D InSAR技術(shù)已在地表形變監(jiān)測(cè)方面得到了廣泛應(yīng)用[11-13],SBAS-D InSAR技術(shù)[14]克服了常規(guī)D InSAR監(jiān)測(cè)中的時(shí)空失相干因素影響,在城市緩慢地表形變[15]、斷層蠕動(dòng)[16]、火山區(qū)域緩慢地表形變[17-20]監(jiān)測(cè)等方面表現(xiàn)出了極大的潛力和優(yōu)勢(shì)。所以本文嘗試?yán)?SBAS-D InSAR技術(shù)提取騰沖火山區(qū)形變時(shí)間序列,探討騰沖火山活動(dòng)的狀態(tài)。
SBAS算法中,差分干涉相位定義為[14]:
其中 x和 r為像元坐標(biāo),λ為雷達(dá)波長(zhǎng),Δd為雷達(dá)視線方向地表形變,B⊥為垂直基線,θ為 SAR視角, Δz為地形殘差,Δφatm為大氣延遲相位,Δφn為其他噪聲相位。顯然,SBAS算法將差分干涉相位分為地表形變相位、地形殘差相位、大氣延遲相位以及其他噪聲相位 4部分。其數(shù)據(jù)處理流程如圖 1所示。
圖1 SBAS-D InSAR數(shù)據(jù)處理流程Fig.1 Flow chart of SBAS-D InSAR processing
文獻(xiàn)[14-17]的研究表明,SBAS算法的精度可達(dá)毫米級(jí),與 GPS觀測(cè)結(jié)果也吻合得較好。SBAS算法由于在影像自由組合干涉時(shí)限制了時(shí)間基線和空間基線,保證了每幅干涉圖的高相干性,又利用奇異值分解的方法將多個(gè)干涉圖子集聯(lián)合進(jìn)行最小二乘求解,增加了時(shí)間采樣。SBAS算法將外部DEM誤差產(chǎn)生的相位分離出來(lái),降低了使用外部DEM所引入的誤差。充分考慮到大氣延遲相位對(duì)形變結(jié)果的影響,SBAS算法通過(guò)時(shí)空域?yàn)V波的方法削弱了大氣延遲相位。相對(duì)于永久散射體技術(shù),小基線集技術(shù)獲取到的形變序列在空間上更為連續(xù),從而可以應(yīng)用于監(jiān)測(cè)地殼長(zhǎng)時(shí)間緩慢變形。
3.1 SAR影像的選取
騰沖火山區(qū)地貌復(fù)雜,所處緯度偏低,植被覆蓋茂密,很容易造成 SAR影像失相干,所以本文選用L波段 JERS-1 SAR數(shù)據(jù),L波段數(shù)據(jù)能夠較好地克服植被覆蓋所造成的失相干問(wèn)題。數(shù)據(jù)處理結(jié)果表明,在騰沖火山區(qū),時(shí)間基線仍不超過(guò)一年才具有較好的相干性。6景 JERS-1 SAR影像,獲取時(shí)間分別為 1995年 4月 10日、1995年 8月 20日、1995年 11月 16日、1996年 3月 27日、1996年 8月 6日、1997年 4月 27日,給定時(shí)間基線為 365天,空間垂直基線為 2 km作為限制條件,進(jìn)行自由干涉組合,得到6個(gè)可用的干涉對(duì),時(shí)間跨越了 1995—1997年。可用干涉對(duì)的時(shí)間基線和空間基線關(guān)系如圖 2所示。
圖2 干涉對(duì)基線組合Fig.2 Baseline combination of all the interferomeric pairs
3.2 常規(guī) D InSAR技術(shù)處理
對(duì)每個(gè)符合時(shí)間、空間基線限制條件的干涉對(duì)均利用兩軌法 D InSAR技術(shù)處理,處理平臺(tái)為瑞士GAMMA遙感公司開(kāi)發(fā)的GAMMA軟件。處理步驟主要包括:原始 SAR信號(hào)成像、圖像配準(zhǔn)與干涉、外部DEM模擬、基線估計(jì)與精化、相干系數(shù)計(jì)算、干涉圖自適應(yīng)濾波、相位解纏。外部地形數(shù)據(jù)使用美國(guó)NASA(美國(guó)宇航局)公布的 SRT M DEM。另外,為了削弱噪聲,增加相干性,同時(shí)為了檢測(cè)相對(duì)大尺度的形變信息,對(duì)干涉圖進(jìn)行了多視處理,多視因子為9:18。相位解纏采用最小費(fèi)用網(wǎng)絡(luò)流算法,為了與GPS觀測(cè)資料進(jìn)行比對(duì),采用騰沖縣城附近為解纏起始區(qū)域。
3.3 SBAS技術(shù)處理
對(duì)所有符合條件的 6個(gè)差分干涉對(duì)的解纏相位進(jìn)行小基線集技術(shù)處理。首先,通過(guò)差分干涉相位估計(jì)得到研究區(qū)域的線性形變相位,并轉(zhuǎn)換為 LOS向形變速率,其中 3座全新世火山的線性形變速率分別為黑空山 -2.8 mm/a,打鷹山 -6.6 mm/a,馬鞍山 -0.82 mm/a;同時(shí)得到了殘差地形相位,發(fā)現(xiàn)其值均為 -0.1~0.1 rad,換算為 LOS向的形變約為 -1.4~1.4 mm(對(duì)于 JERS-1)。然后從原始差分干涉相位中減去線性形變相位和殘差地形相位,再對(duì)殘余相位進(jìn)行解纏,利用奇異值分解方法求解非線性形變相位,最后聯(lián)合線性形變相位和非線性形變相位得到時(shí)間序列累積形變相位。對(duì)于大氣延遲相位的估計(jì),首先在時(shí)間域上進(jìn)行均值濾波處理,然后在空間域上采用小波分解的方法,分離出高頻成分,從而實(shí)現(xiàn)了大氣延遲相位的分離。圖 3是1995年 8月 20日的大氣延遲相位圖,從圖 3可以看出,最大達(dá) 1.5 rad。
圖3 大氣延遲相位(1995-08-20)Fig.3 I mage on the atmospheric phase delay(1995-08-20)
3.4 SBAS形變結(jié)果的質(zhì)量評(píng)價(jià)
圖4為L(zhǎng)OS向形變時(shí)間序列。根據(jù)雷達(dá)成像的幾何關(guān)系,地面的上、東、北 3方向形變對(duì) LOS向形變的貢獻(xiàn)可表示為:
其中,θ為雷達(dá)脈沖入射角,φ為衛(wèi)星軌道方位角。
本文選用的 GPS站點(diǎn)的形變速率值基本含蓋了 InSAR數(shù)據(jù)所覆蓋的范圍。將 GPS觀測(cè)的 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觀測(cè)時(shí)間與 SAR成像時(shí)間不同,得到的速率不同,但是兩者的一致性很好,最大差值約 1 cm/a,差值的標(biāo)準(zhǔn)差為 0.65 cm/ a。對(duì)比結(jié)果表明 SBAS技術(shù)對(duì)地表形變的監(jiān)測(cè)精度可達(dá)亞厘米級(jí)。
3.5 騰沖火山區(qū)形變時(shí)間序列分析
圖4直觀地表現(xiàn)了騰沖火山區(qū) 1995—1997年的地表形變 (LOS向)隨時(shí)間的演化過(guò)程。從全局來(lái)看,騰沖火山區(qū)發(fā)生了較大范圍的區(qū)域性地表變形,大部分地區(qū)以下降為主要形變特征,其中下降相對(duì)較大的地方出現(xiàn)在黑空山西南以及騰沖縣城以南地區(qū),累積最大約 8 cm;研究區(qū)的西南部在 1995—1996年發(fā)生了隆升,最大達(dá) 4 cm,1996年后又回落;從局部來(lái)看,個(gè)別斷裂在形變場(chǎng)上有明顯顯示,反映出斷裂兩側(cè)有明顯的差異活動(dòng),例如西南部的北西向膽札-高田斷裂以及中部的南北向固?hào)|-騰沖斷裂。圖 5給出了垂直于膽札-高田斷裂走向的時(shí)間序列形變剖面曲線 (圖 4A-A’剖面),可以看出,斷裂東北側(cè)遠(yuǎn)離斷層逐漸下沉,而西南側(cè)遠(yuǎn)離斷層則表現(xiàn)為隆升,斷裂兩側(cè)存在明顯的形變梯度,是斷裂活動(dòng)的反映。這一結(jié)果與李春光等[7]計(jì)算得到的 1998—1999年形變特征一致,與野外近場(chǎng)數(shù)字化臺(tái)站記錄結(jié)果所顯示的火山地震活動(dòng)主要發(fā)生在打鷹山-馬鞍山一帶以及騰沖縣城西南部地區(qū)[1]也是一致的,表明該斷裂現(xiàn)今仍在活動(dòng)。值得注意的是,膽札-高田斷裂的西南側(cè)在 1996年后發(fā)生了約 4 cm的不規(guī)則形變變化,體現(xiàn)了該斷裂的較強(qiáng)活動(dòng)性[4]。
從圖 4還可以定性地得出:1995—1997年,相對(duì)于騰沖縣城,打鷹山的東北和西南兩側(cè)表現(xiàn)為明顯的差異性形變特征,馬鞍山地區(qū)形變不明顯,黑空山一直處于下降狀態(tài)。圖 6給出了穿過(guò)打鷹山的形變時(shí)間序列剖面曲線(圖 4B-B’剖面),從圖 6可以看出,打鷹山的東北側(cè)表現(xiàn)為上升,而西南側(cè)表現(xiàn)為下降,累積最大差異達(dá) 10 cm,這種巨大的差異性運(yùn)動(dòng)特征預(yù)示著打鷹山下方可能存在北西向隱伏斷裂,且活動(dòng)性較強(qiáng)。這一結(jié)果與 1998—2000年的水準(zhǔn)測(cè)量結(jié)果[8]一致,1998—1999年與 1999—2000年打鷹山的南北兩側(cè)的垂直形變差異性非常明顯,推測(cè)壓力源的位置不是固定不變的,沿破火口周圍遷移。表 2列出了馬鞍山和黑空山兩座火山的時(shí)間序列形變量值,可以看出,兩地區(qū)均表現(xiàn)為緩慢下降,但馬鞍山地區(qū)下降幅度相對(duì)較小,累積不到 1 cm,表明其在 1995—1997年活動(dòng)性不顯著。黑空山附近表現(xiàn)為連續(xù)性下沉,累積最大超過(guò) 3 cm,表明其在 1995—1997年活動(dòng)性不強(qiáng)。
圖4 地表形變(LOS向)時(shí)間序列(1995—1997年)Fig.4 Ti me seriers of defor mation on LOS(1995-1997)
圖5 膽札-高田斷裂形變時(shí)間序列剖面Fig.5 Cross section on time series of defor mation of Danzha-Gaotian fault
圖6 打鷹山形變時(shí)間序列剖面Fig.6 Cross section on time series of deformation of Dayingshan volcano area
表2 兩座全新世活動(dòng)火山的形變時(shí)間序列Tab.2 Ti me series of deformation of two Holocene volcanoes
綜合火山區(qū)的整體形變時(shí)間序列和 3座全新世活動(dòng)火山的形變時(shí)間序列演化特征分析,可以得出:騰沖火山區(qū)現(xiàn)今地殼運(yùn)動(dòng)比較活躍,形變特征復(fù)雜多樣;打鷹山地區(qū)活動(dòng)性比較顯著,下方可能存在隱伏斷裂;膽札-高田斷裂兩側(cè)形變相對(duì)較大,差異性顯著,地球物理資料揭示其下方存在活動(dòng)性較強(qiáng)的巖漿囊可能是導(dǎo)致其差異性形變的主要因素。
1)相對(duì)于常規(guī) D InSAR技術(shù),SBAS-D InSAR技術(shù)分離了大氣延遲相位以及地形殘差相位等干擾因素,從而能夠獲得高精度的地表形變時(shí)間序列 (LOS向),騰沖火山區(qū)的試驗(yàn)結(jié)果表明 SBAS-D InSAR技術(shù)能夠以亞厘米級(jí)的精度檢測(cè)區(qū)域地表形變。
2)膽札-高田斷裂兩側(cè)形變差異性顯著,可能與其下方存在的地殼巖漿囊的活動(dòng)有關(guān)。
3)打鷹山地區(qū)地表形變揭示其下方可能存在隱伏斷裂。
致謝 感謝路中博士提供 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 李大明,李齊,陳文寄.騰沖火山區(qū)上新世以來(lái)的火山活動(dòng)[J].巖石學(xué)報(bào),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 葉建慶,等.騰沖火山區(qū)地震群的活動(dòng)特征[J].地震地質(zhì),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 趙慈平,冉華,陳坤華.由相對(duì)地?zé)崽荻韧茢嗟尿v沖火山區(qū)現(xiàn)存巖漿囊 [J].巖石學(xué)報(bào),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 樓海,等.云南騰沖火山區(qū)上部地殼三維地震速度層析成像[J].地震學(xué)報(bào),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 李成波,等.騰沖火山區(qū)的 GPS形變特征[J].地球物理學(xué)進(jìn)展,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 李春光,等.騰沖火山區(qū)的形變特征[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 胡亞軒,等.騰沖火山區(qū)地表垂直形變分析[J].大地測(cè)量與地球動(dòng)力學(xué),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 胡亞軒,王慶良,趙慈平.騰沖火山區(qū)形變分析[J].國(guó)際地震動(dòng)態(tài),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 施行覺(jué),等.以垂直形變資料反演騰沖火山區(qū)巖漿活動(dòng)性的初步研究 [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作地面沉降監(jiān)測(cè)的試驗(yàn)研究[J].大地測(cè)量與地球動(dòng)力學(xué),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 喬學(xué)軍,郭利民.新疆伽師強(qiáng)震區(qū)的 InSAR觀測(cè)研究[J].大地測(cè)量與地球動(dòng)力學(xué),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 喬學(xué)軍,等.當(dāng)雄Ms6.6地震的 InSAR觀測(cè)及斷層位錯(cuò)反演[J].大地測(cè)量與地球動(dòng)力學(xué),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
國(guó)家自然科學(xué)基金(40974062);地震行業(yè)科研重點(diǎn)專項(xiàng)(200908029)
季靈運(yùn),男,1982年生,博士生,主要研究方向?yàn)?InSAR技術(shù)在地震、火山形變中的應(yīng)用.E-mail:dinsar010@163.com
P225.1
A