劉智勇，王娟娟

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

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

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

圖1 技術(shù)路線圖

1 核試驗(yàn)區(qū)概況與數(shù)據(jù)來源

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

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

表1 Sentinel-1 A不同成像模式具體參數(shù)

2 SBAS-InSAR技術(shù)測(cè)量核試驗(yàn)后地表形變

2.1 SBAS-InSAR技術(shù)原理與方法

本文采用SBAS-InSAR 技術(shù)對(duì)研究區(qū)進(jìn)行沉降監(jiān)測(cè)。SBAS 技術(shù)是由早期小基線差分干涉演變而來的一項(xiàng)反演時(shí)序地表形變的數(shù)據(jù)處理方法[13]。

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

2.2 第6次朝鮮核試驗(yàn)爆后地表形變提取

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

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

圖2 時(shí)空基線連接

圖3 剔除前后相干系數(shù)矩陣

圖4 地表形變速率圖

3 核試驗(yàn)場(chǎng)地表形變分析

3.1 核試驗(yàn)場(chǎng)整體形變分析

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

圖5 朝鮮核試驗(yàn)場(chǎng)時(shí)間序列形變圖

3.2 核試驗(yàn)場(chǎng)地表形變單點(diǎn)分析

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

圖6 朝鮮核試驗(yàn)場(chǎng)形變單點(diǎn)分析

4 結(jié) 語

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