長白山天池火山2015—2022年InSAR變形與活動狀態分析

熊國華 季靈運* 劉傳金

1)中國地震局第二監測中心,西安 710054 2)國家遙感中心地質災害研究部,北京 100036

0 引言

長白山天池火山位于吉林省東南部的中朝邊境線上,是一個休眠了300多年、仍具有噴發危險的復合式活火山(劉若新等,1992)。目前從外形上看,長白山天池火山的錐體形貌完整,破火山口完好,從地貌上可以清晰地辨別火山造錐的4個階段。迄今為止,長白山天池火山可查證的爆炸式噴發活動至少有3次,規模最大的一次噴發發生于公元946年,曾造成巨大損失。長白山天池火山被認為是中國最危險且最具有潛在噴發危險的火山之一(Xuetal.,2013),其周圍分布多條斷裂,斷裂附近、長白山天池火山與望天鵝火山周圍地震活動頻繁。據統計,2002—2005年期間,長白山天池火山發生異常擾動,其間地表出現隆升現象,地震活動亦變得頻繁(Xuetal.,2012)。其后的2010—2020年間,地震發生頻率有所降低,但仍處于活動狀態,火山活動會引發火山灰、地震及海嘯等多種災害,且天池火山的蓄水量高達20.4億立方米,一旦噴發,水體可能會與巖漿發生作用形成巖漿-氣水噴發,并容易引發火山泥石流,造成更大危害。

跟蹤天池火山的變化、明確火山活動情況及規律對認識火山噴發的危險性具有非常重要的意義。大地測量作為一種重要的技術手段,可通過地表形變信息反映地殼內部巖漿的狀態,反演獲取巖漿房的位置和深度,進而使我們了解火山活動變化情況(牛玉芬,2020)。傳統的地表形變監測手段通常是基于單點進行的局部監測(劉俊清等,2012; 顧國華等,2019; 郭明瑞等,2022),具有精度高的優勢,但無法實現廣域監測。遙感技術作為一種非接觸式的監測手段,可實現大范圍監測,且合成孔徑雷達干涉測量技術(Interferometric Synthetic Aperture Radar,InSAR)在監測緩慢蠕動形變方面具有獨特優勢,其可利用歷史存檔數據獲取地表變化區域的位置、變形趨勢及特征,在火山形變監測中也已得到廣泛應用(季靈運等,2013; 郭倩,2020; 楊云飛,2021; 劉媛媛等,2022; Liuetal.,2023)。不同學者采用InSAR技術對長白山天池火山開展了地表形變監測工作。1992—1998年,Kim等(2004)基于InSAR技術監測到長白山天池火山周圍呈現速度為3mm/a的緩慢隆升變形。2004—2010年期間,InSAR技術探測到在長白山天池火山的北坡和西坡有一定程度的地表變形。北坡的監測結果顯示,2004—2006年,垂直形變量≤5mm,形變量級較小; 2006—2009年,地表出現明顯的隆升趨勢,在破火山口邊緣監測到12mm的隆升變形; 2009—2010年,開始發生明顯下沉,量級達10mm。西坡的監測結果表明,2004—2008年,地表呈現隆升態勢,之后整體發生緩慢下沉(季靈運,2012; 何平等,2015)。2018—2020年,Trasatti等(2021)發現長白山天池火山的東南側和西南側均發生明顯變形,東南側的形變速率可達20mm/a,西南側的形變速率達-20mm/a。上述監測結果表明長白山天池火山持續發生著地表變形,火山下方可能持續存在巖漿活動。為了解長白山天池火山的活動狀態,本文采用時序InSAR技術,利用2015—2022年間的Sentinel-1A/B升、降軌數據,對長白山天池火山開展了動態監測,并利用Mogi點源模型反演火山巖漿房的深度與體積變化,以期為天池火山的危險性研究提供更多數據支撐。

1 技術原理

1.1 時序InSAR技術

本文采用SBAS InSAR(Small Baseline Subset InSAR)技術對長白山天池火山開展監測。它是一種針對分布式散射目標的時序InSAR技術(Berardinoetal.,2002),使用多干涉對組合進行高相干點的相位分析,以獲取可靠的連續面狀結果,對于相干性比較弱的區域也能夠獲取較好的監測效果。SBAS InSAR基于多主影像的SAR數據處理模式,可以在很大程度上保證2幅干涉影像之間的相干性,增加冗余觀測量并提高影像利用率。其基本原理是假設區域有N張SAR影像,根據小基線原理,對滿足要求的影像進行干涉。在去除平地效應和地形相位后生成差分干涉圖,采用最小費用流算法對每幅差分干涉圖進行相位解纏。第i幅干涉圖中像素(x,y)的解纏相位δφ(x,y)可表示為

δφ(x,y)=δφ(x,y)(tB)-δφ(x,y)(tA)≈δφ(def,x,y)+δφ(e,x,y)+δφ(α,x,y)+δφ(n,x,y)

(1)

式中,δφ(x,y)(tB)與δφ(x,y)(tA)分別表示SAR采集時間tB、tA時的相位(tA

1.2 Mogi點源模型

Mogi點源模型也被稱為火山模型,是反演火山巖漿房參數最常用的模型(陳國滸等,2008; 靳錫波等,2020)。其原理是將火山內部的巖漿壓力源與火山活動引起的地表形變建立聯系的一種數學模型。該模型有2個基本假設(Mogi,1958):1)地殼是一個無限大的彈性體; 2)地表會由于壓力源的壓力或體積變化而產生位移。當地下巖漿發生活動后,會造成工作面上的介質發生變化,破壞平衡狀態,導致地表發生形變?；诖思僭O,使用地表形變作為約束,可反演獲得火山巖漿房的體積變化及位置分布情況。

2 數據及處理流程

2.1 數據

本文圍繞長白山天池火山建立研究區,從歐洲航天局(1)https:∥search.asf.alaska.edu/。獲取了64景升軌Sentinel-1A影像; 2個軌道的降軌SAR影像,其中軌道號為32的數據有172景,軌道號為134的數據有180景?？紤]到長白山冬季積雪會影響InSAR形變監測,選擇影像時主要考慮每年6—9月的數據。數據覆蓋情況如圖1所示,基本參數見表1。同時,采用精密軌道數據去除軌道誤差(2)https:∥s1qc.asf.alaska.edu/aux＿poeorb/。,并采用30m分辨率的SRTM DEM地形數據(3)https:∥e4ftl01.cr.usgs.gov/。模擬并去除地形相位。

表1 Sentinel-1A數據參數列表

圖1 研究區的地理位置及SAR影像覆蓋圖

2.2 處理流程

本文使用GAMMA商業軟件,采用SBAS InSAR技術,在保證空間分辨率的基礎上減少噪聲的干擾,獲取了連續、面狀的形變時間序列監測結果。在處理過程中,為了降低噪聲、提高信噪比,首先設置距離向與方位向之比為5︰1進行多視處理,并設置時間與空間基線閾值,滿足閾值的影像生成干涉對。由于獲取的升軌的數據量較少,設置的時間基線為1500d,空間基線為±250m,最后參與形變計算的干涉對組合如圖2a 所示。設置降軌數據的時間、空間基線分別為500d、±250m,2個軌道的干涉圖基線組合分別如圖2b、2c所示。然后,借助30m分辨率的SRTM DEM模擬并去除地形相位,獲取差分干涉圖。采用Goldstein濾波,設置大小為64、32、16的窗口相繼進行濾波,以提高干涉圖的相干性并降低噪聲。采用最小費用流方法進行相位解纏,設置相干性閾值,剔除大部分由于積雪及植被覆蓋導致相干性較低的像元,分別獲取升、降軌的解纏干涉圖。采用通用型衛星雷達大氣改正系統(Generic Atmospheric Correction Online Service for InSAR,GACOS)(4)http:∥www.gacos.net。去除大氣誤差,通過使用二次多項式擬合的方法去除趨勢誤差。對去除誤差后的解纏干涉圖進行挑選,剔除質量較差的干涉圖。最后,采用相位堆疊算法(Stacking InSAR)獲取年平均形變速率圖,并進行地理編碼,獲得WGS-84坐標系下的形變結果。

圖2 干涉圖的時空基線圖

獲取InSAR的形變速率后,基于Mogi點源模型反演火山的巖漿房參數,得到巖漿房的深度及體積變化。在計算過程中,假設巖漿壓力源的參考點位置為(42.0075°N,128.0559°E),使用貝葉斯方法反演模型參數,設置迭代次數為106次,以確定巖漿房的最佳參數及不確定性。采用馬爾科夫鏈-蒙特卡洛算法結合Metropolis-Hastings算法計算參數的后驗概率分布。反演輸入的約束數據為InSAR視線向形變速率,在開始反演計算之前,采用均勻降采樣方法對數據進行重采樣,以提高反演效率。然后通過設置巖漿房的中心位置x、y,深度z,體積變化率ΔV的上、下限及相應的步長,約束模型參數的反演。

3 火山形變場

3.1 2015—2022年InSAR形變速率

本文獲得的長白山天池火山的InSAR形變速率如圖3所示。其中,正值表示形變靠近衛星飛行方向,負值表示形變遠離衛星飛行方向。天池火山在中國境內的區域植被較為茂密,干涉圖的相干性較差,靠近朝鮮的區域植被較少,且火山口附近裸巖較多,獲取的干涉點主要分布在火山口附近。一般而言,巖漿活動引起的地表變形表現為整體的趨勢性形變,局部形變信號通常不是由火山運動引起的,可能是由人類活動或其他因素導致。

圖3 形變速率圖

在2015—2022年期間,升、降軌監測到的年平均形變速率為-8～6mm/a。升軌的形變比降軌的量級更小,主要是由于SAR的幾何特征及地形差異所致。升、降軌的形變方向較為一致,說明長白山天池火山的形變主要發生在垂直方向上。升軌視線向形變速率較小(圖3a),主要的形變分布在破火山口處,形變分布均勻,地表變形主要為沉降,形變速率為-2～-1mm/a。2個降軌數據的形變分布表現較為一致(圖3b,c),在火山口附近出現整體的趨勢性形變,在遠離火山口的區域存在局部變形?；鹕娇诟浇憩F為整體緩慢下沉,量級不大,視線向形變速率為-4～-2mm/a。局部形變主要分布在火山口東側的2個區域,用A、B表示,其中A區域的最大形變速率為-8mm/a,B區域的最大形變速率為-6mm/a。根據形變速率結果可知,天池火山總體表現為沉降狀態,火山口形變較大,向四周逐漸減小,呈擴散狀,表明天池火山可能存在深部活動。

由于未獲得實地測量數據,本文采用內符合精度評定的方法評價InSAR監測結果的精度。結合SAR影像的局部入射角,求取了3個軌道的SAR數據在垂直方向上的形變速率,分別統計不同軌道間公共點的垂直形變的差值,計算殘差的偏差及標準差(圖4)。經統計,殘差總體均呈正態分布,升軌(ASC54)與降軌(DSC32)垂向形變的殘差偏差為0.78mm/a(圖4a),但在正值區間存在一定殘差。降軌垂向形變(DSC32、DSC134)的殘差偏差為-1.74mm/a(圖4b),這可能是由于參考點的差異導致結果存在較小的系統偏差,但殘差分布更為集中。不同軌道數據間垂向形變的殘差均為毫米級,達到InSAR監測的精度要求,說明本文使用InSAR技術監測長白山天池火山的結果具有較好的可靠性。

圖4 觀測值間的殘差分布直方圖

3.2 2015—2022年InSAR形變時間序列

為進一步了解長白山天池火山各部位形變的時空變化特征及趨勢,使用SBAS InSAR技術進行計算,獲取了時間段為2014年11月—2021年12月、共計172景降軌(DSC32)SAR影像的形變時間序列。圖5為降軌數據(DSC32)在雷達視線方向上的累計形變時間序列。從圖5可以看出,長白山天池火山在監測期間存在持續變形。2014年11月—2017年6月,火山總體形變速度較慢,從2018年6月起形變加速; 2020年7月—2021年12月期間,形變加速明顯。在近7a的監測時間內,累計形變量達-40mm?？傮w而言,近7a以來,長白山天池火山地表持續發生緩慢變形,主要發生在火山口附近,且于近期有加速趨勢。

圖5 降軌(DSC32)形變時間序列

提取橫穿火山口的2條剖線,定量分析長白山天池火山的地表變化情況。剖線方向分別為NW-SE和NE-SW,具體位置見圖5中的2021年12月11日累計形變圖。獲取的剖線高程及形變時間序列如圖6所示。長白山天池火山水面的位置相對較低,水面高約2200m,高度由中間向兩側增加,在2條剖線的時間序列監測結果中,天池火山口處存在明顯變形。

圖6 剖面累計形變圖

圖7 長白山火山地震數據統計圖(改自Meng et al.,2022)

根據2條剖線的累計形變時間序列可知,在天池火山口的南側監測到較大程度的變形。根據NW-SE剖線的形變時間序列可知,在火山口東南側存在顯著形變,視線向累計形變量達-38mm,發生在天池火山口的內側。在火山口的北西側存在小幅度地表隆升?；鹕礁浇傮w形變主要表現為沉降,距離火山口的位置越遠,形變變化趨勢越趨于穩定。形變與地形表現出相關性,在地形起伏劇烈的區域,對應的形變也呈現出一定程度的起伏。根據NE-SW剖線的形變時間序列可知,火山口西南側的形變比東北側量級更大,量級最大的形變出現在火山口山體外側的西南側斜坡中部,視線向形變量達-38mm。

根據前人的研究可知,長白山天池火山的地表形變在歷史監測中表現出周期性變化。同時,結合地震活動性可以了解火山的活動狀態。1992—1998年,該區地表沒有發生明顯變形,地震活動弱,火山巖漿活動并不活躍; 2002—2006年,地表發生約6.5cm的隆升變形(Jietal.,2021),同時根據得到的地震數據表明,此時段內中強地震活動性增強,最大震級為4.4級,地震發生的頻率約為600次/a,表明該時段內火山巖漿活動較為活躍(Panetal.,2021; Mengetal.,2022)。在接近4a的隆升變化之后,2008—2011年觀測到天池火山處地表開始發生沉降,并且地震活動減弱,發生地震的頻率為約70次/a,此階段火山巖漿活動較弱(Panetal.,2021)。本文的研究表明,2015—2022年期間,火山口附近的地表視線向年平均形變速率約為-4～-2mm/a,累計形變量可達-40mm。同時,統計得到的震中距火山口30km內的地震數量顯示,2015—2020年共發生地震195次,表明該時段內地震活動性弱,最大地震的震級為2.2級,說明在此期間火山巖漿活動并不活躍。但值得注意的是,2020年12月22日出現了由38次小地震組成的火山震群事件,地震活動有增多的趨勢(盤曉東等,2022)。通常情況下,火山擾動被認為是周期性的巖漿活動。目前,長白山天池火山雖然沒有近期噴發的危險,但仍需要持續開展監測。

4 巖漿房參數反演

上文通過近7a的SAR影像獲取的InSAR形變監測結果,反映了長白山天池火山的地表動態變化過程,為研究其巖漿活動變化奠定了基礎。本文使用Mogi模型,聯合2個降軌形變場數據反演了巖漿房參數。在InSAR結果中2個降軌數據在形變分布特征上表現出較好的一致性,表明監測結果的可靠性更高。同時,由于使用升軌數據計算形變結果時SAR影像數據量少,干涉對的時空基線較長,可能導致干涉圖之間的相干性相對較弱,使結果中存在一定誤差。因此,為保證結果的準確性,在進行參數反演時,主要利用2個降軌數據獲取的地表形變結果作為約束,通過多次反演實驗,得到擬合殘差最小的巖漿房參數(表2)。

表2 Mogi點源模型反演的幾何參數

根據反演得到的參數可以看出,巖漿房的位置處于火山口下偏西的區域(41.9928°N,128.0361°E),深約6km,巖漿房的年體積變化率為-3.3×105m3/a,表明長白山天池火山的巖漿在近幾年變化緩慢,體積變化表現為收縮。不同學者對長白山天池火山的巖漿房空間位置及活動性開展了研究,并給出了不同結果。地球物理探測的相關研究者發現在長白山天池火山的周圍存在明顯的P波低速異常,認為在下地殼或更深的位置存在巖漿房,且處于活動狀態,并進一步估計了巖漿房的位置和深度。其中,基于三維層析成像技術推測出的長白山天池火山巖漿房位于火山口西側位置(楊卓欣等,2005); 利用小地震精定位數據反演結果推測得到的巖漿房位于天池火山地下約5km處(吳建平等,2007); 通過鉆孔測量技術并基于長白山天池火山地溫梯度判斷的巖漿房位于深5～8km處(潘波等,2017; 錢程等,2022)。大地電磁觀測的反演結果表明,長白山天池火山下方6～8km處存在巖漿通道,在火山口的北部7km處存在巖漿囊(仇根根等,2014)。大地測量的相關學者通過地表監測數據反演獲取了長白山天池火山的巖漿房分布及巖漿體積的變化速率,結果表明,巖漿房并非位于火山的正下方,巖漿活動會引起地表的形變。其中,InSAR反演結果顯示長白山天池火山的巖漿房位于火山口下偏西的區域,深6～7km(何平等,2015),與本文反演的結果相近。GNSS與InSAR結合的反演結果顯示,長白山天池火山的巖漿房位于深6.9km處(陳國滸等,2008)。以上研究表明,長白山天池火山區下方存在巖漿房,深度為5～8km,由于觀測方式、數據處理手段、觀測時間段的不同及不同空間范圍內的巖漿擾動均會導致反演的巖漿房深度出現偏差,因此,對于巖漿房所處的空間位置與展布形態目前還不夠清晰。本文的觀測結果與GNSS、地球物理探測等技術得到的結果相比,能夠更加全面地覆蓋天池火山靠近朝鮮一側的區域,更有利于說明巖漿房的空間位置與活動狀態。

本研究通過反演的參數正演出長白山天池火山的地表形變,并計算模擬值與觀測值的殘差,以驗證模擬的準確性。圖8給出了各數據的觀測值、模擬值與殘差值?？梢钥闯?在反演結果中,趨勢性形變均被較好擬合,說明巖漿房活動變形引起的地表變化可被反演獲取,表明反演結果具有可靠性。同時,我們對反演的殘差進行了統計(圖9),發現殘差均分布在零值附近,DSC32殘差的偏差為0.13mm/a,DSC134殘差的偏差為0.23mm/a。殘差結果在一定程度上說明了數據擬合較好,但也可以看到在負值區間存在一定殘差,這可能是Mogi模型反演時無法擬合到局部形變所致。

圖8 反演得到的觀測值、模擬值與殘差值圖

圖9 殘差統計圖

5 討論

5.1 火山巖漿房收縮機理分析

長白山天池火山被普遍認為是爆裂式噴發型火山,根據目前多種資料佐證,在長白山天池火山地下的淺地殼中存在巖漿房,且在地下深部存在巖漿儲存庫(許建東等,2022)。淺部巖漿房為位于地殼中的部分或全部熔融體,其來源多是由深部巖漿儲存庫的巖漿多次注入?；鹕絿姲l通常是淺部巖漿房中的巖漿在浮力的驅動下通過火山通道涌出地表的過程(段政等,2022)。一般情況下,巖漿房處于平衡狀態,當有新的巖漿注入時,巖漿房承受的壓力發生變化,呈現增加趨勢,巖漿房的晶體平衡被打破,從而產生新的反應,使得巖漿處于活動狀態。因此,監測巖漿房體積的變化情況有助于了解火山活動的強弱。巖漿房體積變化一般有膨脹與收縮2種類型?；鹕綆r漿房的體積膨脹一般由巖漿入侵或增壓導致(Luetal.,2003; Jietal.,2013)。體積收縮的原因較為復雜,可能由以下幾種情況導致:1)巖漿體的緩慢冷卻和結晶(Luetal.,2005); 2)火山氣體由火山噴氣孔或溫泉釋放,導致壓力降低(Sepeetal.,2007; Fournier,2008; Jietal.,2018); 3)淺儲層的巖漿流出(Cervelletal.,2002); 4)火山強烈活動、噴發而導致巖漿損失(Luetal.,2000)。本研究獲取的長白山天池火山在監測時間段內的地表形變主要表現為沉降,部分區域存在小幅度隆升,反演得到的巖漿房體積變化表現為收縮,在監測期間未發生噴發,推測體積的變化很可能是由于巖漿的冷卻和結晶所致。

5.2 火山形變與巖漿活動時間序列

長白山天池火山被認為是中國最具危險性的活火山(洪漢凈等,2007)?；鹕阶冃慰梢越沂净鹕綆r漿是否發生擾動,根據火山地表形變研究,長白山天池火山表現出持續的巖漿活動。針對長白山天池火山的變形監測,主要包括水準測量、GPS、InSAR 3種方式?；鹕街苓吂步ㄓ?5個GPS站,覆蓋面積達1200km2,測量始于2000年。精密水準測量是從2002年起在北坡和西坡進行的。InSAR技術利用歷史影像追溯地表變形情況,最早于1992年開始監測,隨著目前影像種類和數量的增多,在長白山天池火山變形監測中起到了關鍵的作用。

多項研究表明,自1992年至今,長白山天池火山地下巖漿房的演化過程為平靜—發生擾動—歸于平靜(圖10)。1992—1998年期間,長白山天池火山沒有發生明顯變形(Jietal.,2021)。1995—1998年間,利用InSAR技術結合點源模型的反演結果顯示,火山地下可能存在雙巖漿房,深度分別為7.9km、5.5km,并得到巖漿房的體積變化量分別為1.25×106m3、3.6×106m3,體積變化表現為膨脹(陳國滸等,2008)。2002—2010年,長白山天池火山地區的變形監測數據更加豐富。2002—2006年,GPS與InSAR監測到長白山天池火山附近的巖漿發生動蕩,巖漿活動頻繁,出現明顯的異常信號。2002—2005年,采用GPS監測數據結合點源模型的反演結果顯示,巖漿房的體積變化率為4.6×106m3/a,分布在深4.4km處,體積變化表現為膨脹(Xuetal.,2012; Trasattietal.,2021)。2006—2009年,InSAR監測結合模型的反演結果顯示,巖漿房分布在6.3km深度處,體積變化率為1.11×106m3/a,體積變化仍為膨脹,但變化速率出現明顯下降,說明自2006年起巖漿活動變弱,開始變得平靜(何平等,2015)。2009—2011年,GPS監測到火山出現明顯下沉,結合點源模型的反演結果顯示,巖漿房以-1.4×106m3/a的速度發生收縮(Xuetal.,2012; Trasattietal.,2021)。2015—2022年,InSAR監測到火山仍在下沉,但速度較為緩慢,結合點源模型的反演結果顯示,深度為6km處的巖漿房以-0.33×106m3/a的速度發生收縮,推測長白山天池火山繼2002—2006年的擾動期過后一直處于相對較為平靜的狀態。

圖10 巖漿活動時間序列示意圖

6 結論

本文利用2015—2022年升、降軌Sentinel-1A/B SAR衛星數據計算獲取了長白山天池火山近7a的形變時間序列及年平均形變速率,并以此地表形變監測結果作為約束,反演獲取了長白山天池火山的巖漿房位置分布與體積變化率,得到以下結論:

(1)基于InSAR的時間序列變化結果表明,長白山天池火山在監測期間,地表變形表現為下沉,近7a間火山口附近的累計形變可達-40mm,形變速率約為-4～-2mm/a。

(2)基于Mogi模型反演獲取的巖漿房參數表明,長白山天池火山的巖漿房處于火山口下偏西的位置,深約6km,巖漿房的體積變化率為-3.3×105m3/a。反演結果表明,在近7a的監測時間內,長白山天池火山的巖漿活動并不活躍,結合地震活動性分析認為目前長白山天池火山的巖漿活動較弱。

(3)長白山天池火山于1992—2022年期間經歷了平靜—擾動—平靜的巖漿活動過程。1992—1998年,長白山火山沒有明顯形變,火山巖漿活動較弱; 于2002—2005年監測到該區出現明顯的地表隆升變形,巖漿房體積顯著膨脹,之后巖漿活動逐漸減弱。