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

劉媛媛 趙振宇 晏 霞 陳人杰

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

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

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

1 技術(shù)原理

1.1 D-InSAR技術(shù)

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

1.2 Mogi模型

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

(1)

(2)

(3)

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

(4)

式中，Ux、Uy、Uz為坐標(biāo)原點的位移分量，r為Mogi源模型中心與地面點在地表投影坐標(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ù)以及對應(yīng)的精密軌道文件，時間基線為12 d。同時，選用分辨率高達30 m的DEM數(shù)據(jù)，用于減弱外部DEM引起的誤差。

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

圖1 LOS向形變

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

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

γ(h)=nugget+

(5)

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

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

2.3 模型反演與分析

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

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

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

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

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

Iguchi等[11]研究櫻島昭和火山口2007～2011年火山噴發(fā)的構(gòu)造活動，在研究中假設(shè)一個位于南岳火山口下方4 km處的壓力源，利用傾斜儀和引伸儀的徑向應(yīng)變變化，估計壓力源體積的月變化量與火山灰重量之間的關(guān)系(圖4)。結(jié)果表明，昭和火山口可能通往南岳下方的巖漿儲層或與南岳火山口的管道相連。由圖4可見，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)建一個3Mogi源模型，用于計算以往多次火山噴發(fā)事件的體積變化量。表2為文獻[12]與本研究計算的M源體積變化量，表3為M源參數(shù)，圖5為壓力源平面位置。由表2可見，以往火山噴發(fā)事件導(dǎo)致的M源體積最大變化量與最小變化量分別為-0.33×106m3和-0.017×106m3，此次火山噴發(fā)事件的體積變化量大約是以往M源體積最小變化量的8倍和最大變化量的0.42倍。由圖5可見，Mogi源分別位于艾拉火山口下方(A源)9.6 km、北岳下方(K源)3.3 km和南岳下方(M源)0.7 km處。此次火山噴發(fā)事件反演得到的Mogi源位置與Hotta等[12]得到的M源位置較為接近，M1源由本文計算得到，M2源為Hotta等[12]計算得到。M1位于南岳火山口東側(cè)，M2位于昭和火山口東側(cè)，二者在位置上存在一定差距，可能是由于SAR衛(wèi)星在監(jiān)測過程中容易受到誤差的影響，雖然本文在實驗過程中對誤差進行了改正和削弱，但無法完全消除。盡管如此，從最終結(jié)果來看，本文反演的參數(shù)結(jié)果仍然具有較好的參考價值。

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

表3 M源參數(shù)

圖5 壓力源平面位置

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

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

3 結(jié) 語

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

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

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