基于時序InSAR的地下核試驗爆后形變監測與分析

劉智勇，王娟娟

（1. 武漢大學遙感信息工程學院，湖北 武漢 430072；2. 66444部隊，北京 100094）

天然地震、火山活動等自然災害，工程爆破、地下核武器試驗等人類活動都會引起地面震動。通常可根據震源深度來區分該事件是自然災害還是人為活動[1]：地震發生于中上地殼的脆弱區，離地表較深[2]；地下核試驗或工程爆破引起的地面震動位于淺層地表。朝鮮民主主義人民共和國（朝鮮） 在2006-10-09—2017-09-03 共進行6 次被證實的地下核試驗，其中最后一次核試驗威力遠超前5 次，一次又一次的地下核試驗引發了國際社會的廣泛關注[3]。因其引起的地面震動與地震類似，很多學者采用監測地震的方法對地下核試驗進行監測，如Gaber H[4]等采用震源深度與位置、頻率成分、復雜度、譜比值、Ms-mb 震級比等判別方法區分地震和爆炸，成功區分了朝鮮2016-09-09 的核試驗與2000-09-08 的4.9 級地震；Ichinose G A[5]等基于朝鮮第2~5 次核爆炸的表面波和體波對核試驗位置進行了相對定位，分析了核試驗發生的位置區域；謝小碧[6]等利用地震學方法對朝鮮6 次地下核爆炸展開研究，結果表明朝鮮地下核試驗表現出典型淺源爆炸的特征。綜上所述，基于地震學方法可確定核爆的相對位置與當量，但受限于地震臺網的數量與密度，核試驗的詳細位置與爆炸范圍估算存在較大的不確定性。

20世紀70年代遙感技術開始出現于地震災害監測研究領域。隨著高分辨率的合成孔徑雷達（SAR）遙感衛星的發射升空，對于植被覆蓋度高且人員無法實地到達的區域，合成孔徑干涉雷達（InSAR）技術能以高精度、高分辨率、全天候的典型優勢進行觀測。Vincent P[7]等采用InSAR對地下核試驗引起的表面位移進行觀測，干涉圖時間序列顯示震后沉降形變受區域地質和水文環境影響呈現不同的變形速率和大小；Wang T[8]等利用TerraSAR-X SAR影像分析了第6次朝鮮核試驗引起的地表位移，并建立彈性模型模擬山體的上升和倒塌，估計爆炸的位置、深度和能量；Xu X[9]等基于沉降坑形成過程的縮放規律，利用Wang T[8]估計的位置、深度和能量，開發了一種爆炸模型裝置，驗證了沉降帶和坍塌坑的半徑接近經驗公式和InSAR 觀測值；Yocky D A[10]等對比了極化SAR 測量爆炸引起的表面形變與加速計記錄的地面運動預測，結果表明極化SAR表面相干空間范圍與測量和預測的表面形變速率有較強的相關性；Zhu M[11]等利用Mogi 模型研究了朝鮮第4 次地下核試驗引起的地面位移，得出校正電離層延遲能有效減小干涉圖中大尺度相位斜坡影響的結論。本文以朝鮮核試驗場（朝鮮東北部咸境北道吉洲郡）為研究區域，基于33 景Sentinel-1 雷達降軌影像，通過SBAS-InSAR 技術建立研究區LOS 向地表沉降模型，獲取朝鮮核試驗區域山谷、山脊處2017-09-10—2018-12-22 的累積形變量，利用時序InSAR 技術為非天然地震所引起的地表形變監測提供可靠的數據。技術路線見圖1。

圖1 技術路線圖

1 核試驗區概況與數據來源

朝鮮核試驗區位于朝鮮東北部咸鏡北道吉州郡豐溪里地區，地下核試驗發生在山脈中，海拔高程約為500~2 300 m，地形起伏較大，核試驗場地下基巖主要由花崗巖等硬巖組成。萬塔山山體結構穩定，樹林密布，大部分區域植被覆蓋度較高、隱蔽性強，成為地下核試驗的適宜區域，本文選取41.275°～41.31°N、129.055°～129.125°E范圍作為研究區。

由于SAR衛星覆蓋范圍大、空間分辨率高且具備全天候運行能力，已廣泛用于地表位移監測研究。歐空局的Sentinel-1 A/B衛星安裝了適用于大規模成像的SAR傳感器，具有短回訪周期（6 d）和相當高的空間分辨率。Sentinel-1 A 衛星的工作模式主要包括波模式（WV）、極寬幅模式（EW）、干涉寬幅模式（IW）和條帶模式（SM）4種[12]（表1）。本文采用IW模式數據進行地表形變監測分析。

表1 Sentinel-1 A不同成像模式具體參數

2 SBAS-InSAR技術測量核試驗后地表形變

2.1 SBAS-InSAR技術原理與方法

本文采用SBAS-InSAR 技術對研究區進行沉降監測。SBAS 技術是由早期小基線差分干涉演變而來的一項反演時序地表形變的數據處理方法[13]。

首先根據研究范圍將覆蓋研究區的N+1幅SAR影像按時間順序排列，并通過合適的時間與空間閾值生成SAR影像干涉像對，可在一定程度上避免時空相干的影響，進而使用最小二乘法計算得到每個子集合的形變時間序列[14]；再根據所生成連接像對的關聯結果，依次對所有像對進行處理，由于干涉圖中還存在地形和殘余相位，需引入外部DEM 和選取合適（不能在形變處選點）的GCP控制點來去除和估算；然后對干涉圖進行軌道精煉和重去平，并再次利用最小費用流方法進行相位解纏和濾波；最后利用奇異值分解法計算得到地表位移的時間序列形變信息，并進行地理編碼獲取所需坐標系下（WGS84）的時序地表變形和地表形變速率[15]。

2.2 第6次朝鮮核試驗爆后地表形變提取

地下核試驗的發生必然會導致地下巖層的受力狀態發生變化，對于震級大的核試驗，爆后的持續下沉也十分明顯。已有資料表明，朝鮮第6 次核試驗核爆釋放的能量相當于發生了6.3 級地震，將導致核試驗后地下空腔內的溫度和壓力較長時間高于外界，在重力作用下，山體會逐漸壓實，宏觀上表現為地表下沉。由于核試驗區域地形復雜，山體表面植被覆蓋度高，難以獲取第6 次朝鮮核試驗爆后地表形變信息。在考慮監測對象地表形變時空變化特征的前提下，保持影像干涉對相位高相干性就顯得十分重要。

時序InSAR 處理的時空基線連接圖見圖2，可以清楚發現，很多影像的相干系數均值低于0.4，這主要是由于研究區植被覆蓋度高，C波段的Sentinel數據穿透性不足，導致干涉效果非常差。因此，本文在影像連接的空間基線范圍內對連接的干涉圖進行篩選，剔除相干性差的干涉相對。剔除前后干涉系數的變化情況見圖3，通過保留高相干的干涉像對，提高相干性。地理編碼后的監測結果見圖4，結果表明，2017年9 月—2018年12 月核試驗區在核爆后出現地表下沉，具體表現為靠近山體南側區域的形變速率最快，最大形變達到100 cm/a；山體東西兩側也出現部分滑動，位移速率約為50 cm/a。

圖2 時空基線連接

圖3 剔除前后相干系數矩陣

圖4 地表形變速率圖

3 核試驗場地表形變分析

3.1 核試驗場整體形變分析

由核試驗場時間序列形變結果（圖5）可知，核試驗場東北側存在異常變形，且隨著時間的推移，形變緩慢增大；選取高相干性區域，并進一步篩選測量結果中保留下來的高相干點位，大部分區域分布在無植被覆蓋的裸巖或建筑區域，觀測獲得的時間序列形變表明在朝鮮第6次地下核試驗后的15個多月里，核爆后地表形變向遠離衛星方向移動（LOS地表形變為負）。整體來看，萬塔山山體地表形變出現不同程度的沉降，山坡上位移最明顯，主要原因可能是核試驗發生導致地下爆室的收縮，在山坡上出現緩慢局部滑坡；山頂形變量并不是很大，可能是由于山體構造的基巖主要為花崗巖，結構較穩定，導致形變速率緩慢。此外，萬塔山之外出現了許多正值點位，其原因可能是提取形變時根據高程范圍選取掩膜區域，導致少許噪聲點未去除，該部分點位并不在研究區域內，對研究萬塔山試驗區形變影響不大。

圖5 朝鮮核試驗場時間序列形變圖

3.2 核試驗場地表形變單點分析

為定量分析萬塔山地表形變趨勢，本文在關鍵區域選擇典型點位進行單點時序序列形變分析。由于SBAS-InSAR 技術獲取的是地表LOS 向形變，需轉換為地表沉降值LOS cosθ，本文采用Sentinel數據的平均入射角（約為39°）。圖6a中紅色三角形位置為形變區域山頂上的點（41.3007N、129.07E）；圖6b為該三角形點在2017-09-10—2018-12-22 的形變量變化情況，隨著時間增加，形變量逐漸增大，形變速率變化較平均，最大累積沉降量約為30 mm。圖6c中紅色三角形位置為形變區域山坡上的點（41.3007N、129.07E）；圖6d 為該三角形點在在2017-09-10—2018-12-22的形變量變化情況，隨著時間增加，形變量逐漸增大再減小，出現明顯“兩段”特征，前段形變速率較大，后段形變速率較緩和，定量驗證了核試驗發生后地表下沉由快變慢的現象。

圖6 朝鮮核試驗場形變單點分析

4 結 語

本文以朝鮮第6 次地下核試驗為研究對象，基于長時序Sentinel-1 雷達影像集，開展了核試驗區爆后地表形變監測分析，主要結論為：①時序InSAR 技術可用于爆后核試驗場長時序地表形變監測，為非可達地區核試驗監測提供了重要的技術支撐，監測結果對于爆后核試驗場山體穩定性分析具有重要價值；②基于時序InSAR 技術發現朝鮮第6 次核試驗爆后長時間序列累積地的表形變量在2017年9 月—2018年12 月核試驗場各區域出現不同量級沉降，其中山頂出現約30 mm 的下沉，邊坡的形變量更高，約為70~100 mm。