熊國(guó)華 季靈運(yùn)* 劉傳金

1)中國(guó)地震局第二監(jiān)測(cè)中心,西安 710054 2)國(guó)家遙感中心地質(zhì)災(zāi)害研究部,北京 100036

0 引言

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

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

1 技術(shù)原理

1.1 時(shí)序InSAR技術(shù)

本文采用SBAS InSAR(Small Baseline Subset InSAR)技術(shù)對(duì)長(zhǎng)白山天池火山開展監(jiān)測(cè)。它是一種針對(duì)分布式散射目標(biāo)的時(shí)序InSAR技術(shù)(Berardinoetal.,2002),使用多干涉對(duì)組合進(jìn)行高相干點(diǎn)的相位分析,以獲取可靠的連續(xù)面狀結(jié)果,對(duì)于相干性比較弱的區(qū)域也能夠獲取較好的監(jiān)測(cè)效果。SBAS InSAR基于多主影像的SAR數(shù)據(jù)處理模式,可以在很大程度上保證2幅干涉影像之間的相干性,增加冗余觀測(cè)量并提高影像利用率。其基本原理是假設(shè)區(qū)域有N張SAR影像,根據(jù)小基線原理,對(duì)滿足要求的影像進(jìn)行干涉。在去除平地效應(yīng)和地形相位后生成差分干涉圖,采用最小費(fèi)用流算法對(duì)每幅差分干涉圖進(jìn)行相位解纏。第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采集時(shí)間tB、tA時(shí)的相位(tA

1.2 Mogi點(diǎn)源模型

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

2 數(shù)據(jù)及處理流程

2.1 數(shù)據(jù)

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

表1 Sentinel-1A數(shù)據(jù)參數(shù)列表

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

2.2 處理流程

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

圖2 干涉圖的時(shí)空基線圖

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

3 火山形變場(chǎng)

3.1 2015—2022年InSAR形變速率

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

圖3 形變速率圖

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

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

圖4 觀測(cè)值間的殘差分布直方圖

3.2 2015—2022年InSAR形變時(shí)間序列

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

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

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

圖6 剖面累計(jì)形變圖

圖7 長(zhǎng)白山火山地震數(shù)據(jù)統(tǒng)計(jì)圖(改自Meng et al.,2022)

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

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

4 巖漿房參數(shù)反演

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

表2 Mogi點(diǎn)源模型反演的幾何參數(shù)

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

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

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

圖9 殘差統(tǒng)計(jì)圖

5 討論

5.1 火山巖漿房收縮機(jī)理分析

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

5.2 火山形變與巖漿活動(dòng)時(shí)間序列

長(zhǎng)白山天池火山被認(rèn)為是中國(guó)最具危險(xiǎn)性的活火山(洪漢凈等,2007)?；鹕阶冃慰梢越沂净鹕綆r漿是否發(fā)生擾動(dòng),根據(jù)火山地表形變研究,長(zhǎng)白山天池火山表現(xiàn)出持續(xù)的巖漿活動(dòng)。針對(duì)長(zhǎng)白山天池火山的變形監(jiān)測(cè),主要包括水準(zhǔn)測(cè)量、GPS、InSAR 3種方式。火山周邊共建有15個(gè)GPS站,覆蓋面積達(dá)1200km2,測(cè)量始于2000年。精密水準(zhǔn)測(cè)量是從2002年起在北坡和西坡進(jìn)行的。InSAR技術(shù)利用歷史影像追溯地表變形情況,最早于1992年開始監(jiān)測(cè),隨著目前影像種類和數(shù)量的增多,在長(zhǎng)白山天池火山變形監(jiān)測(cè)中起到了關(guān)鍵的作用。

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

圖10 巖漿活動(dòng)時(shí)間序列示意圖

6 結(jié)論

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

(1)基于InSAR的時(shí)間序列變化結(jié)果表明,長(zhǎng)白山天池火山在監(jiān)測(cè)期間,地表變形表現(xiàn)為下沉,近7a間火山口附近的累計(jì)形變可達(dá)-40mm,形變速率約為-4～-2mm/a。

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

(3)長(zhǎng)白山天池火山于1992—2022年期間經(jīng)歷了平靜—擾動(dòng)—平靜的巖漿活動(dòng)過程。1992—1998年,長(zhǎng)白山火山?jīng)]有明顯形變,火山巖漿活動(dòng)較弱; 于2002—2005年監(jiān)測(cè)到該區(qū)出現(xiàn)明顯的地表隆升變形,巖漿房體積顯著膨脹,之后巖漿活動(dòng)逐漸減弱。