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

劉媛媛 趙振宇 晏 霞 陳人杰

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

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

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

1 技術原理

1.1 D-InSAR技術

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

1.2 Mogi模型

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

(1)

(2)

(3)

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

(4)

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

2 數據處理與模型反演

2.1 數據處理

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

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

圖1 LOS向形變

2.2 數據誤差估計

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

γ(h)=nugget+

(5)

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

圖2 InSAR數據估計非變形區域的誤差特征

2.3 模型反演與分析

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

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

表1 Mogi模型參數的約束信息及反演結果

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

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

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

圖4 壓力源體積變化與火山灰的關系[11]

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

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

表3 M源參數

圖5 壓力源平面位置

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

圖6 地下壓力源與巖漿供應

3 結 語

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

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

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