基于D-InSAR的櫻島火山噴發(fā) 形變場(chǎng)解算與反演分析

劉媛媛 趙振宇 晏 霞 陳人杰

1 東華理工大學(xué)測(cè)繪工程學(xué)院，南昌市廣蘭大道418號(hào)，330013 2 自然資源部環(huán)鄱陽(yáng)湖區(qū)域礦山環(huán)境監(jiān)測(cè)與治理重點(diǎn)實(shí)驗(yàn)室，南昌市廣蘭大道418號(hào)，330013

火山噴發(fā)是一種重大自然災(zāi)害。據(jù)統(tǒng)計(jì)，2020年全球共有24個(gè)國(guó)家的89座火山發(fā)生噴發(fā)，其中Ⅲ級(jí)以上預(yù)警的火山達(dá)21座[1]。火山噴發(fā)時(shí)會(huì)釋放出大量的火山灰和火山氣體，這些物質(zhì)會(huì)隨大氣環(huán)流被帶至高空，從而影響氣候變化[2]。櫻島火山(31°35′N(xiāo)，130°39′E)位于日本鹿兒島市，東接大隅群島，西隔錦江灣，距離鹿兒島市中心約4 km，南北面皆為鹿兒島灣。火山由北岳(海拔1 117 m)、中岳(海拔1 160 m)、南岳(海拔1 140 m)組成，面積約77 km2，是一座活躍的成層火山[3]，包含南岳火山口、昭和火山口及中央火山錐兩側(cè)的多個(gè)寄生火山口[4]。櫻島火山在過(guò)去的100 a里發(fā)生過(guò)幾次不同類(lèi)型的噴發(fā)。根據(jù)日本氣象廳公布的櫻島火山活動(dòng)情況，2020-06下旬開(kāi)始觀測(cè)到山體膨脹的緩慢變動(dòng)，2020-07下旬膨脹停止，隨后在2020-08-09 05:41:53監(jiān)測(cè)到南岳火山口發(fā)生噴發(fā)，當(dāng)日發(fā)生的A型地震達(dá)68次，火山噴發(fā)產(chǎn)生的煙霧高度達(dá)5 km，櫻島部分地區(qū)的降灰量達(dá)300 g/m2。

火山噴發(fā)往往是由地殼內(nèi)部巖漿活動(dòng)引起的，根據(jù)地表形變可以有效反演出地下巖漿的活動(dòng)特征，對(duì)預(yù)測(cè)火山噴發(fā)及災(zāi)害防御有重要意義。傳統(tǒng)的監(jiān)測(cè)手段無(wú)法及時(shí)獲取大范圍的地表形變信息，而合成孔徑雷達(dá)差分干涉D-InSAR技術(shù)作為一種非接觸測(cè)量手段，精度可達(dá)cm級(jí)甚至亞cm級(jí)，具有全天時(shí)、視域廣、密度大等特點(diǎn)[5]，已被廣泛應(yīng)用于各類(lèi)火山噴發(fā)引起的地表形變監(jiān)測(cè)中[6]。2015年櫻島火山巖脈入侵事件中，相關(guān)學(xué)者使用多視角InSAR數(shù)據(jù)進(jìn)行D-InSAR測(cè)量，解算巖脈入侵事件導(dǎo)致的三維形變場(chǎng)[7]。為監(jiān)測(cè)櫻島火山的最新活動(dòng)狀態(tài)，本文擬采用D-InSAR技術(shù)，利用覆蓋日本鹿兒島市櫻島火山的2景Sentinel-1A數(shù)據(jù)研究2020-08-09櫻島火山噴發(fā)事件，并利用點(diǎn)源Mogi模型反演火山下方巖漿源的特征。本文對(duì)研究火山噴發(fā)的性質(zhì)和災(zāi)害分布具有重要意義。

1 技術(shù)原理

1.1 D-InSAR技術(shù)

InSAR技術(shù)的重要應(yīng)用領(lǐng)域之一是地表形變監(jiān)測(cè)[8]，而D-InSAR技術(shù)是InSAR技術(shù)的拓展，可以理解為SAR衛(wèi)星沿著相同的軌道重復(fù)運(yùn)行，獲取同一地區(qū)不同時(shí)間的影像。將2幅不同時(shí)段主、副影像的相位值進(jìn)行差分處理，即可得到差分干涉相位圖。此時(shí)干涉圖中包含的相位有參考橢球面相位、地形相位、形變相位、大氣相位和噪聲相位，其中占比較大的相位為參考橢球面相位、地形相位和形變相位，D-InSAR主要提取的是形變相位。

1.2 Mogi模型

對(duì)此次火山噴發(fā)事件進(jìn)行多次模擬實(shí)驗(yàn)后，本文最終選用點(diǎn)源Mogi模型[9]來(lái)恢復(fù)地表形變。Mogi模型的參數(shù)分別為源中心坐標(biāo)在地表的投影(x0，y0)，源中心的深度d以及巖漿源體積變化量ΔV。假設(shè)地面坐標(biāo)系為(x，y，z)，巖漿源的空間位置為(x0，y0，d)，根據(jù)Mogi模型擬合地面點(diǎn)三維位移與源模型參數(shù)的公式為：

(1)

(2)

(3)

r2=(x-x0)2+(y-y0)2

(4)

式中，Ux、Uy、Uz為坐標(biāo)原點(diǎn)的位移分量，r為Mogi源模型中心與地面點(diǎn)在地表投影坐標(biāo)間的距離。

2 數(shù)據(jù)處理與模型反演

2.1 數(shù)據(jù)處理

選取覆蓋櫻島火山區(qū)域2020-07-28(噴發(fā)前)與2020-08-09(噴發(fā)后)21h的2景C波段降軌Sentinel-lA單視復(fù)數(shù)SLC影像數(shù)據(jù)以及對(duì)應(yīng)的精密軌道文件，時(shí)間基線(xiàn)為12 d。同時(shí)，選用分辨率高達(dá)30 m的DEM數(shù)據(jù)，用于減弱外部DEM引起的誤差。

由于火山周?chē)L(zhǎng)的植被會(huì)造成失相干，為抑制數(shù)據(jù)噪聲的影響以及提高相干性，本文采用距離向?yàn)?、方位向?yàn)?的多視比，并使用一種局部自適應(yīng)濾波的方法。在此基礎(chǔ)上，采用基于規(guī)則格網(wǎng)的最小費(fèi)用流算法進(jìn)行相位解纏，相關(guān)系數(shù)閾值設(shè)為0.25，最終通過(guò)地理編碼獲得地理坐標(biāo)系下櫻島火山LOS向一維形變場(chǎng)(圖1)。由圖可見(jiàn)，此次火山噴發(fā)事件造成的形變主要集中在中央火山錐中心，地表沉降主要發(fā)生在火山錐中心地帶，最大沉降量與平均沉降量分別為5.5 cm和2.85 cm；抬升主要發(fā)生在火山邊緣區(qū)域，最大抬升量與平均抬升量分別為5 cm與2.24 cm，其中北部的隆起可能與艾拉火山口下方的巖漿活動(dòng)有關(guān)。

圖1 LOS向形變

2.2 數(shù)據(jù)誤差估計(jì)

通常情況下，SAR衛(wèi)星在對(duì)觀測(cè)區(qū)域進(jìn)行數(shù)據(jù)采集時(shí)會(huì)受到大氣相位的影響，與大氣和地面狀態(tài)相關(guān)的每個(gè)InSAR圖像都具有特定的InSAR數(shù)據(jù)誤差[9]。本文構(gòu)建噪聲的完整協(xié)方差矩陣：首先從原始數(shù)據(jù)中刪除一個(gè)線(xiàn)性斜坡，以校正可能的長(zhǎng)波大氣效應(yīng)；然后估計(jì)去趨勢(shì)數(shù)據(jù)集上的實(shí)驗(yàn)半變異函數(shù)；最后使用nugget模型擬合出指數(shù)函數(shù)。指數(shù)函數(shù)的方差與協(xié)方差矩陣的定義為：

γ(h)=nugget+

(5)

式中，γ(h)為每個(gè)數(shù)據(jù)點(diǎn)之間任意給定距離h處的協(xié)方差矩陣，range為數(shù)據(jù)點(diǎn)在空間上的相關(guān)距離，nugget為空間獨(dú)立噪聲的級(jí)別，sill為當(dāng)range趨于無(wú)窮大時(shí)半方差的最大值。圖2為根據(jù)式(5)計(jì)算出的InSAR數(shù)據(jù)誤差特性，由圖可見(jiàn)，得到的距離值、噪聲級(jí)別和方差最大值分別為849.845 8 m、1.612 2×10-6m2和8.964 7×10-6m2。

圖2 InSAR數(shù)據(jù)估計(jì)非變形區(qū)域的誤差特征

2.3 模型反演與分析

本文以D-InSAR處理得到的位移信息為觀測(cè)數(shù)據(jù)，使用貝葉斯方法檢索櫻島火山噴發(fā)的Mogi源最佳參數(shù)和不確定性，并使用馬爾可夫鏈蒙特卡洛(MCMC)算法結(jié)合M-H算法搜索Mogi 模型各項(xiàng)參數(shù)的后驗(yàn)概率[10]。

反演中使用的參考系坐標(biāo)原點(diǎn)為南岳火山口地理坐標(biāo)(130.66°E，31.58°N)，X為東西方向的位移，向東為“+”；Y為南北方向的位移，向北為“+”；Z為巖漿源的深度。由圖1可見(jiàn)，火山周?chē)罅康闹脖粫?huì)造成幾何失真產(chǎn)生空白，故選擇中央火山錐區(qū)域進(jìn)行模型反演。反演開(kāi)始前，通過(guò)設(shè)定模型參數(shù)的約束信息來(lái)提高反演速度(表1)。此外，采用D-InSAR得到的數(shù)據(jù)集為面狀數(shù)據(jù)，原始數(shù)據(jù)量較大，需要在反演的初始階段運(yùn)用四叉樹(shù)算法對(duì)其進(jìn)行降采樣處理，采樣后的數(shù)據(jù)點(diǎn)為106個(gè)，保證了形變場(chǎng)的特征。

表1 Mogi模型參數(shù)的約束信息及反演結(jié)果

經(jīng)過(guò)1.0×106次迭代后得到模型的最佳擬合參數(shù)，由表1可知，東西向位移為0.402 km，南北向位移為0.069 km，深度為1.061 km，巖漿源位于南岳與昭和火山口之間，比較靠近昭和火山口。模型各項(xiàng)參數(shù)的聯(lián)合概率如圖3所示，由圖可見(jiàn)，總體殘差很小。

①X/m, ②Y/m, ③Z/103 m, ④ΔV/105m3, ⑤視線(xiàn)方向恒定偏移量/m， ⑥軌道殘差X/10-6 m， ⑦軌道殘差Y/10-6 m，⑧迭代次數(shù)/104次圖3 聯(lián)合概率分布

Iguchi等[11]研究櫻島昭和火山口2007～2011年火山噴發(fā)的構(gòu)造活動(dòng)，在研究中假設(shè)一個(gè)位于南岳火山口下方4 km處的壓力源，利用傾斜儀和引伸儀的徑向應(yīng)變變化，估計(jì)壓力源體積的月變化量與火山灰重量之間的關(guān)系(圖4)。結(jié)果表明，昭和火山口可能通往南岳下方的巖漿儲(chǔ)層或與南岳火山口的管道相連。由圖4可見(jiàn)，Mogi源的月體積變化主要集中在±0.5×106m3范圍內(nèi)，而根據(jù)Mogi源反演得到的此次火山噴發(fā)的體積參數(shù)為-0.139×106m3，與Iguchi等[11]得到的月體積變化量較為接近。

圖4 壓力源體積變化與火山灰的關(guān)系[11]

2016年Hotta等[12]使用GNSS、應(yīng)變與傾斜數(shù)據(jù)構(gòu)建一個(gè)3Mogi源模型，用于計(jì)算以往多次火山噴發(fā)事件的體積變化量。表2為文獻(xiàn)[12]與本研究計(jì)算的M源體積變化量，表3為M源參數(shù)，圖5為壓力源平面位置。由表2可見(jiàn)，以往火山噴發(fā)事件導(dǎo)致的M源體積最大變化量與最小變化量分別為-0.33×106m3和-0.017×106m3，此次火山噴發(fā)事件的體積變化量大約是以往M源體積最小變化量的8倍和最大變化量的0.42倍。由圖5可見(jiàn)，Mogi源分別位于艾拉火山口下方(A源)9.6 km、北岳下方(K源)3.3 km和南岳下方(M源)0.7 km處。此次火山噴發(fā)事件反演得到的Mogi源位置與Hotta等[12]得到的M源位置較為接近，M1源由本文計(jì)算得到，M2源為Hotta等[12]計(jì)算得到。M1位于南岳火山口東側(cè)，M2位于昭和火山口東側(cè)，二者在位置上存在一定差距，可能是由于SAR衛(wèi)星在監(jiān)測(cè)過(guò)程中容易受到誤差的影響，雖然本文在實(shí)驗(yàn)過(guò)程中對(duì)誤差進(jìn)行了改正和削弱，但無(wú)法完全消除。盡管如此，從最終結(jié)果來(lái)看，本文反演的參數(shù)結(jié)果仍然具有較好的參考價(jià)值。

表2 M源體積變化量[12]

表3 M源參數(shù)

圖5 壓力源平面位置

火山活動(dòng)的強(qiáng)弱會(huì)隨時(shí)間變化。櫻島火山的常見(jiàn)噴發(fā)模式有2種：間歇性噴發(fā)和連續(xù)性噴發(fā)。6月底至7月底櫻島火山無(wú)噴發(fā)事件，但監(jiān)測(cè)到山體膨脹，最終于2020-08-09發(fā)生火山噴發(fā)，因此可以將其視作一次短期的間歇性噴發(fā)。在此期間火山下方的巖漿由A源入侵至K源，經(jīng)過(guò)一段時(shí)間的累積后由K源遷移至M源，最終在地表噴發(fā)。櫻島的巖漿管道系統(tǒng)由一個(gè)主要巖漿儲(chǔ)層、一個(gè)次要巖漿儲(chǔ)層和一個(gè)連接管道組成，在火山噴發(fā)后，巖漿各個(gè)儲(chǔ)層都會(huì)降低，這可能會(huì)導(dǎo)致巖漿撤退和火山體積的收縮[13]。此次火山噴發(fā)事件反演得到的收縮源與Hotta等[12]得到的南岳源位置較為接近，由此可知，櫻島的巖漿管道系統(tǒng)組成為：艾拉火山口下方約10 km的巖漿儲(chǔ)層連接北岳下方3.3 km的巖漿儲(chǔ)層，經(jīng)過(guò)南岳下方1.016 km的淺層儲(chǔ)層，通過(guò)管道連接至山頂火山口(圖6)，其中北岳下方的K源位置由Hotta等[12]提供。

圖6 地下壓力源與巖漿供應(yīng)

3 結(jié) 語(yǔ)

1)最大沉降區(qū)域?yàn)榛鹕藉F中心地帶，最大沉降量與平均沉降量分別為5.5 cm和2.85 cm；抬升主要發(fā)生在火山邊緣區(qū)域，最大抬升量與平均抬升量分別為5 cm和2.24 cm。

2)根據(jù)Mogi模型反演得到的巖漿源深度和體積收縮量分別為1.016 km和-0.139×106m3，巖漿源位于南岳下方，與南岳火山口的噴發(fā)活動(dòng)有關(guān)。

櫻島火山是一座活火山，近年來(lái)南岳火山口和昭和火山口處于活躍狀態(tài)，同時(shí)也伴隨著密集的地震與火山噴發(fā)現(xiàn)象，間歇性與連續(xù)性的噴發(fā)影響著居民的正常生活，因此許多學(xué)者致力于研究此火山的構(gòu)造活動(dòng)。然而在中央火山錐兩側(cè)還存在寄生火山口，因此仍然無(wú)法對(duì)具體的巖漿供應(yīng)路徑作出詳細(xì)解釋。為了更好地理解火山活動(dòng)的全過(guò)程，還需要利用更加全面的數(shù)據(jù)進(jìn)行研究與討論。