哨兵-2號(hào)光學(xué)影像地表形變監(jiān)測(cè)：以2016年MW7.8新西蘭凱庫(kù)拉地震為例

賀禮家，馮光財(cái)，馮志雄，高 華

1. 中南大學(xué)地球科學(xué)與信息物理學(xué)院雷達(dá)遙感研究室，湖南 長(zhǎng)沙 410083； 2. 東華理工大學(xué)測(cè)繪工程學(xué)院，江西 南昌 330013

雖然光學(xué)數(shù)據(jù)易受云霧、太陽(yáng)照度等因素的影響，限制了其在地表形變監(jiān)測(cè)中的應(yīng)用，但光學(xué)影像可以利用地表紋理信息變化快速地監(jiān)測(cè)大范圍、大量級(jí)的地表形變(如地震、冰川、滑坡等)，并獲取較高的形變監(jiān)測(cè)精度。在氣候條件允許下，它能很好地彌補(bǔ)和克服GPS技術(shù)空間分辨率低、點(diǎn)位覆蓋稀疏以及InSAR技術(shù)易受相位失相干和最大形變梯度等限制因素影響[1-2]。

哨兵-2號(hào)(Sentinel-2)是哥白尼計(jì)劃(Copernicus Programme)下多光譜成像任務(wù)中的光學(xué)衛(wèi)星星座，該星座包含兩顆重復(fù)軌道衛(wèi)星——哨兵-2A及哨兵-2B衛(wèi)星，它們分別于2015年6月和2017年3月發(fā)射升空。此前，已經(jīng)有許多學(xué)者利用各種光學(xué)數(shù)據(jù)進(jìn)行地表形變監(jiān)測(cè)研究，如SPOT[3-9]、ASTER[10]、Landsat[11-15]以及航空影像[16-17]等。與ASTER、Landsat系列等傳統(tǒng)的光學(xué)數(shù)據(jù)相比，哨兵-2號(hào)衛(wèi)星星座具有中高空間分辨率(最高10 m)、較短重訪周期(5 d)及數(shù)據(jù)免費(fèi)等優(yōu)勢(shì)，因此在地表形變監(jiān)測(cè)中的應(yīng)用潛力將更加明顯[18]。由于目前距哨兵-2號(hào)升空時(shí)間較短，其數(shù)據(jù)積累較少，所以目前關(guān)于它的形變監(jiān)測(cè)系統(tǒng)誤差改正和應(yīng)用研究還較少，還尚未有學(xué)者對(duì)哨兵-2號(hào)影像中存在的各種系統(tǒng)誤差源和改正方法進(jìn)行系統(tǒng)的研究。此外，該衛(wèi)星星座可以獲取4個(gè)10 m空間分辨率的波段影像，何種波段更適合地表形變監(jiān)測(cè)以及哨兵-2號(hào)與Landsat8兩種光學(xué)遙感影像地表形變監(jiān)測(cè)精度之間的差異等問題，目前還尚未有分析比較。以上內(nèi)容對(duì)于哨兵-2號(hào)衛(wèi)星數(shù)據(jù)的應(yīng)用和推廣都具有重要意義。

因此，本文首先介紹了哨兵-2號(hào)光學(xué)數(shù)據(jù)監(jiān)測(cè)地表形變的數(shù)據(jù)處理流程,并分析了其影像形變場(chǎng)中系統(tǒng)誤差源的組成成分及改正方法；并在傳統(tǒng)方法的基礎(chǔ)上，采用了改進(jìn)的均值相減法以去除衛(wèi)星姿態(tài)角抖動(dòng)誤差。另外，還對(duì)哨兵-2號(hào)影像中4個(gè)10 m空間分辨率的波段(Band 2/3/4/8)獲取地表形變信息的能力進(jìn)行了分析和比較，通過對(duì)不同波段影像獲取的形變結(jié)果進(jìn)行精度評(píng)定，選取和驗(yàn)證了地表形變監(jiān)測(cè)的最佳波段。接著，分別利用哨兵-2號(hào)和Landsat8影像獲取了2016年11月14日MW7.8新西蘭凱庫(kù)拉地震的同震形變場(chǎng)，并對(duì)兩種光學(xué)影像的觀測(cè)結(jié)果進(jìn)行了對(duì)比分析。最后，本文總結(jié)了哨兵-2號(hào)光學(xué)影像中的系統(tǒng)誤差源去除方法和精度分析，展望了該數(shù)據(jù)在地表形變監(jiān)測(cè)中的應(yīng)用前景。

1 系統(tǒng)誤差改正

1.1 數(shù)據(jù)介紹

哨兵-2號(hào)多光譜成像儀(MSI)采用推掃式成像模式，該衛(wèi)星星座含13個(gè)通道，其中10個(gè)通道是波長(zhǎng)為433～955 nm的可見光近紅外譜段(VNIR)，其他3個(gè)通道是波長(zhǎng)為1360～2280 nm的短波紅外譜段(SWIR)。哨兵-2號(hào)衛(wèi)星的成像幅寬為290 km，光譜分辨率為15～180 nm，空間分辨率包括10、20和60 m。本文采用4個(gè)10 m空間分辨率的可見光近紅外波段影像(Band 2/3/4/8)進(jìn)行試驗(yàn)分析，這4個(gè)光譜波段的光譜和輻射信息見表1。

表1 哨兵-2號(hào)光學(xué)影像光譜波段信息簡(jiǎn)介

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

本文選取新西蘭東北部凱庫(kù)拉縣(Kaikoura)為研究區(qū)域(如圖1所示)。該區(qū)域在2016年11月14日發(fā)生了MW7.8級(jí)大地震。該研究區(qū)域地形地貌豐富，地表破裂明顯，非常適合利用哨兵-2號(hào)光學(xué)數(shù)據(jù)進(jìn)行系統(tǒng)誤差分析和同震形變監(jiān)測(cè)研究。本文從美國(guó)地質(zhì)調(diào)查局(United States Geological Survey，USGS)收集了哨兵-2號(hào)震前數(shù)據(jù)2景，震后數(shù)據(jù)15景，詳細(xì)數(shù)據(jù)信息見表2。此外，為了驗(yàn)證和比較哨兵-2號(hào)影像獲得的形變結(jié)果，本文也收集了同時(shí)期的兩景Landsat8影像。如圖2所示，這17景哨兵-2號(hào)影像可構(gòu)成10對(duì)時(shí)間間隔均為10 d的震后影像對(duì)組合，這些影像對(duì)獲取的形變場(chǎng)中均存在各種系統(tǒng)誤差源。本文利用COSI-Corr軟件包[19-20]和交叉頻譜相關(guān)算法[9,21-22]獲取了地表形變場(chǎng)，然后選取了系統(tǒng)誤差比較明顯的兩個(gè)影像對(duì)2016-12-15—2016-12-25及2016-12-25—2017-01-04(Band8影像)，對(duì)形變場(chǎng)中的各種系統(tǒng)誤差源進(jìn)行系統(tǒng)分析。

圖1 2016年MW7.8新西蘭凱庫(kù)拉地震構(gòu)造背景圖Fig.1 Tectonic setting of 2016 MW7.8 Kaikoura, New Zealand earthquake

圖2 哨兵-2號(hào)相關(guān)影像對(duì)采集日期分布Fig.2 Acquisition dates of the Sentinel-2 correlation image pairs

影像編號(hào)采集時(shí)間太陽(yáng)天頂角/(°)太陽(yáng)方位角/(°)含云量/(%)12016-02-0939.194655.45470.042922016-02-1941.824252.1309032016-12-0529.634457.214935.20342016-12-1529.774759.89600.008652016-12-2530.550261.55514.228862017-01-0431.849962.02915.918472017-02-1041.183757.77731.392282017-02-2043.777454.37940.008892017-02-2343.134750.45400102017-03-0546.008346.860456.1255112017-03-1549.000743.312442.0905122017-04-2161.255035.17070132017-04-2461.077232.034216.7647142017-05-0163.966733.24110.0658152017-05-0463.765130.33020.0283162017-05-1466.141529.11740.2846172017-05-2468.122128.38840.6865

光學(xué)影像地表形變監(jiān)測(cè)完整的數(shù)據(jù)處理流程(如圖3所示)，主要包括如下4個(gè)步驟：

(1) 數(shù)據(jù)預(yù)處理：由于哨兵-2號(hào)數(shù)據(jù)與Landsat8數(shù)據(jù)的影像覆蓋范圍存在差異，為方便對(duì)兩類數(shù)據(jù)結(jié)果進(jìn)行對(duì)比分析，需要對(duì)這兩類影像進(jìn)行裁剪。

(2) 相關(guān)性計(jì)算：以COSI-Corr為數(shù)據(jù)處理平臺(tái)，對(duì)影像對(duì)進(jìn)行亞像素的相關(guān)性匹配(sub-pixel correlation)計(jì)算[23-24]。對(duì)于哨兵-2號(hào)數(shù)據(jù)，參數(shù)設(shè)置如下：搜索窗口大小為32×32像素(地面分辨率約為320×320 m)；移動(dòng)步長(zhǎng)為8×8像素(地面分辨率為80×80 m)；掩膜閾值為0.9；迭代次數(shù)為2次。

(3) 誤差后處理：為了去除二維形變場(chǎng)中的各種系統(tǒng)誤差源，需要進(jìn)行誤差后處理。依次去除失相關(guān)噪聲、軌道誤差、條帶誤差等，最后再通過3×3像素的中值濾波窗口進(jìn)一步降噪[25]。

(4) 形變值重采樣：為方便GMT(generic mapping tools)繪圖和不同類型數(shù)據(jù)比較，需要進(jìn)行坐標(biāo)轉(zhuǎn)換。對(duì)UTM坐標(biāo)系下獲取的二維形變場(chǎng)中的形變值進(jìn)行重采樣，得到WGS-84坐標(biāo)系下的二維形變場(chǎng)。

圖3 光學(xué)影像地表形變監(jiān)測(cè)數(shù)據(jù)處理流程Fig.3 Flow chart of processing the ground deformation obtained from optical images

1.3 誤差分析及糾正

本文以哨兵-2號(hào)2016-12-15—2016-12-25震后影像對(duì)為例，通過相關(guān)性計(jì)算獲取了新西蘭地震覆蓋區(qū)域的東西向和南北向二維地表形變場(chǎng)(如圖4所示)。通過分析發(fā)現(xiàn)哨兵-2號(hào)光學(xué)影像獲取的形變場(chǎng)中主要存在以下幾種系統(tǒng)誤差源：失相關(guān)噪聲、軌道誤差、條帶誤差以及衛(wèi)星姿態(tài)角誤差等(如圖4所示)。此外，光學(xué)影像形變場(chǎng)中一般存在地形陰影誤差[16]，但是考慮到地形陰影誤差的復(fù)雜性，本文不予分析。后文將針對(duì)上述各種系統(tǒng)誤差的分布特點(diǎn)、產(chǎn)生原因以及去除方法進(jìn)行系統(tǒng)性分析。

注：影像形變場(chǎng)處于WGS-84坐標(biāo)系，分別以朝東和朝北方向?yàn)樾巫冎档恼较?。圖4 哨兵-2號(hào)影像對(duì)(2016-12-15—2016-12-25)東西向和南北向影像形變場(chǎng)系統(tǒng)誤差源Fig.4 System error sources in the EW and NS components of the image deformation field from the (2016-12-15—2016-12-25) image pair

1.3.1 失相關(guān)噪聲

失相關(guān)噪聲與地表輻射性能息息相關(guān)。當(dāng)?shù)乇磔椛湫阅茏兓^小時(shí)，可能會(huì)使影像的紋理特征產(chǎn)生微小變化；而當(dāng)輻射性能變化較大時(shí)，則會(huì)引起影像紋理的缺失。失相關(guān)噪聲水平可用信噪比值(SNR)進(jìn)行描述。當(dāng)影像中某個(gè)區(qū)域的信噪比值普遍較低時(shí)，該區(qū)域可能受失相關(guān)噪聲的影響[9,16]。失相關(guān)噪聲產(chǎn)生的原因可歸納為如下幾點(diǎn)：

(1) 時(shí)間失相關(guān)。自然因素：云、雪、植被、水域等覆蓋區(qū)域的變化；人為因素：新增建筑物等。

(2) 地形陰影。同一傳感器在不同季節(jié)采集的影像對(duì)應(yīng)的太陽(yáng)角(包括太陽(yáng)高度角和方位角)不同會(huì)使形變場(chǎng)中某些區(qū)域產(chǎn)生地形陰影。

由于失相關(guān)噪聲的存在，通過對(duì)影像進(jìn)行窗口匹配不僅難以取得好的相關(guān)性計(jì)算結(jié)果，而且還會(huì)大大降低形變值的測(cè)量精度。為了降低失相關(guān)噪聲的影響，一般采用如下兩種解決方案：

(1) 程序自動(dòng)掩膜掉形變場(chǎng)中SNR低于某閾值(經(jīng)驗(yàn)值取0.9)的區(qū)域。

(2) 人為掩膜掉由云、雪、植被等因素引起的范圍較大的失相關(guān)區(qū)域。

1.3.2 軌道誤差

雖然從USGS上下載的哨兵-2號(hào)L1C級(jí)數(shù)據(jù)是正射產(chǎn)品，但是該數(shù)據(jù)沒有經(jīng)過嚴(yán)格的正射校正和幾何校正，因此通過相關(guān)性計(jì)算得到的影像形變場(chǎng)中依然存在明顯的線性長(zhǎng)波長(zhǎng)誤差，即軌道誤差。本文以哨兵-2號(hào)2016-12-25—2017-01-04影像對(duì)為例，利用一次多項(xiàng)式曲面擬合模型去除軌道誤差[26-29]。計(jì)算公式如下

φorbit=a0+a1x+a2y+a3xy

(1)

其中，φorbit為軌道趨勢(shì)項(xiàng)誤差；x為距離向坐標(biāo)；y為方位向坐標(biāo)；a0、a1、a2、a3為待求參數(shù)，可根據(jù)最小二乘平差原理求解。具體步驟為：首先在形變場(chǎng)中遠(yuǎn)離地表破裂帶的非形變區(qū)域均勻的選取若干像素點(diǎn)，根據(jù)這些點(diǎn)的圖像坐標(biāo)及形變值建立多項(xiàng)式誤差曲面以計(jì)算待求系數(shù)a0、a1、a2、a3，然后根據(jù)這些參數(shù)模擬整個(gè)形變場(chǎng)的軌道趨勢(shì)面，從而獲得去除軌道誤差的形變場(chǎng)。如圖5所示，為軌道誤差去除前后的對(duì)比圖，哨兵-2號(hào)東西向和南北向影像形變場(chǎng)中軌道誤差變化緩慢，以常數(shù)為主。

圖5 哨兵-2號(hào)影像對(duì)(2016-12-25—2017-01-04)東西向和南北向影像形變場(chǎng)中軌道誤差去除前后對(duì)比Fig.5 Before and after removing the orbital errors in the EW and NS components of the image deformation field from the Sentinel-2 image pair (2016-12-25—2017-01-04)

1.3.3 條帶誤差

對(duì)于大多數(shù)推掃式成像衛(wèi)星(如SPOT系列、Landsat系列、ASTER等)而言，沿軌道方向都會(huì)產(chǎn)生明顯且均勻分布的線性信號(hào)，該信號(hào)即條帶誤差。條帶誤差產(chǎn)生的原因可歸納為如下兩點(diǎn)：

(1) 正射影像在內(nèi)定向時(shí)沒有經(jīng)過合理的建模，使得多光譜儀器探測(cè)器中的CCD(charge coupled device，電耦合器件)線陣列錯(cuò)位排列[16,30]。

(2) CCD線陣列抖動(dòng)引起的誤差沒有通過建模消除。相鄰且交錯(cuò)排列的CCD線陣列振動(dòng)會(huì)使采集的影像產(chǎn)生偏移，其偏移量大小取決于CCD陣列振動(dòng)的頻率、振幅及相鄰CCD陣列對(duì)同一地理位置成像的時(shí)間間隔[31]。

為了消除CCD線陣列引起的條帶誤差，本文采用傳統(tǒng)的“均值相減法”去除條帶誤差[1,16]。本文選用哨兵-2號(hào)2016-12-25—2017-01-04震后影像對(duì)進(jìn)行計(jì)算分析，這不僅能保證良好的相干性又可以避免同震形變的干擾。條帶誤差去除的具體流程如圖6所示，結(jié)果如圖7所示。

圖6 均值相減法去除條帶誤差數(shù)據(jù)處理流程Fig.6 Flow chart of removing the stripe artifact using the mean subtracting method

圖7 哨兵-2號(hào)影像對(duì)(2016-12-25—2017-01-04)東西向和南北向影像形變場(chǎng)中條帶誤差去除前后對(duì)比Fig.7 Before and after removing the stripe artifact in the E-W and N-S components of the image deformation field from the Sentinel-2 image pair (2016-12-25—2017-01-04)

1.3.4 衛(wèi)星姿態(tài)角誤差

在形變場(chǎng)中，沿交叉軌道方向呈現(xiàn)出周期性均勻平行分布的帶狀信號(hào)，即衛(wèi)星姿態(tài)角抖動(dòng)誤差。該誤差產(chǎn)生的原因可歸納如下：

(1) 在影像采集過程中，航天器振動(dòng)引起的誤差沒有通過建模消除[31]。

(2) 在影像正射校正過程中，由于衛(wèi)星姿態(tài)角欠采樣，導(dǎo)致無(wú)法精確計(jì)算衛(wèi)星姿態(tài)角[32]。

(3) 在東西向形變圖中，該誤差主要是由未被記錄的roll角變化引起，在南北向形變圖中，該誤差主要是由未被記錄的pitch角變化引起[1,9,32]。

為了去除形變場(chǎng)中衛(wèi)星姿態(tài)角抖動(dòng)引起的誤差，傳統(tǒng)采用與去除條帶誤差相同的方法(均值相減法)。但是本文考慮到沿衛(wèi)星交叉軌道方向地形起伏變化、植被覆蓋等因素的影響，沿衛(wèi)星交叉軌道方向的衛(wèi)星姿態(tài)角誤差是變化而非固定不變的，直接采用傳統(tǒng)方法并不能取得良好的誤差去除效果。因此,本文結(jié)合哨兵-2號(hào)衛(wèi)星推掃式成像儀焦平面上共有12個(gè)交錯(cuò)排列的CCD線陣列探測(cè)器的特點(diǎn)，提出改進(jìn)后的“均值相減法”去除衛(wèi)星姿態(tài)角誤差。該改進(jìn)方法的主要思想為：首先定義一個(gè)量η用于定量的衡量形變場(chǎng)中的衛(wèi)星姿態(tài)角誤差水平，可表示為

(2)

式中，m、n分別表示形變圖的最大行數(shù)、最大列數(shù)；a、b分別表示形變圖中某行第j等份中始端、末端像素點(diǎn)的位置參數(shù)；r表示將形變圖每行所含像素點(diǎn)總數(shù)均分成12等份；D表示某一像素點(diǎn)位的形變值；ceil(·)函數(shù)用于求取大于某一數(shù)值的最小整數(shù)。采用該改進(jìn)的方法后，誤差的去除效果有很大的提升(如圖9所示)，而且與傳統(tǒng)方法相比計(jì)算的復(fù)雜程度并沒有太大的增加。本文以哨兵-2號(hào)2016-12-25—2017-01-04影像對(duì)為例，利用改進(jìn)后的均值相減法去除衛(wèi)星姿態(tài)角誤差，具體的處理流程如圖8所示，結(jié)果如圖9所示。

圖8 改進(jìn)后的均值相減法去除衛(wèi)星姿態(tài)角誤差數(shù)據(jù)處理流程Fig.8 Flow chart of removing the satellite attitude jitter distortion using the modified mean subtracting method

2 最佳波段分析

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

哨兵-2號(hào)光學(xué)影像中共有4個(gè)10 m空間分辨率的可見光近紅外波段(Band 2/3/4/8)。雖然這4個(gè)光譜波段具有相同的空間分辨率，但是其波段寬度、SNR值以及輻射特性卻不盡相同[33]，詳見表1。本文利用10對(duì)哨兵-2號(hào)震后影像對(duì)(如圖2所示)對(duì)應(yīng)的4波段影像數(shù)據(jù)進(jìn)行試驗(yàn)分析。由于這些震后影像對(duì)的時(shí)間基線較短(10 d)，地表一般不會(huì)發(fā)生明顯的形變。因此，通過分析不同波段影像獲取的影像形變場(chǎng)中形變值的誤差水平，可以反映出不同波段影像地表形變監(jiān)測(cè)能力的差異。本文利用4個(gè)不同波段的影像對(duì)組合進(jìn)行地表形變監(jiān)測(cè)的最佳波段選取分析。具體的數(shù)據(jù)處理流程主要包括以下2個(gè)步驟：

(1) 相關(guān)性計(jì)算：利用COSI-Corr分別計(jì)算不同波段(Band 2/3/4/8)對(duì)應(yīng)的10個(gè)相關(guān)性影像對(duì)(如圖2所示)的二維地表形變場(chǎng)。計(jì)算參數(shù)統(tǒng)一為：搜索窗口大小為32×32像素；移動(dòng)步長(zhǎng)為8×8像素；掩膜閾值為0.9；迭代次數(shù)為2次。

(2) 誤差后處理：采用上一節(jié)介紹的系統(tǒng)誤差源改正方法，依次去除影像形變場(chǎng)中的各種系統(tǒng)誤差源(失相關(guān)噪聲、軌道誤差、條帶誤差、衛(wèi)星姿態(tài)角誤差等)，再通過3×3的中值濾波窗口進(jìn)一步降噪。

在經(jīng)過以上步驟后，本文對(duì)獲得的形變值的平均值和標(biāo)準(zhǔn)差進(jìn)行統(tǒng)計(jì)分析。統(tǒng)計(jì)結(jié)果顯示不同波段影像獲取的影像形變場(chǎng)中形變值的均值均接近于0 m，這表明圖2中這些震后影像對(duì)的形變場(chǎng)受余震形變影響較小。因此，本文主要對(duì)不同波段獲取的形變場(chǎng)中形變值的標(biāo)準(zhǔn)差進(jìn)行對(duì)比分析，統(tǒng)計(jì)數(shù)據(jù)見表3和圖10。值得注意的是，在統(tǒng)計(jì)形變值標(biāo)準(zhǔn)差時(shí)，本文將形變值為空值(Nan)以及形變值絕對(duì)值較大的異常值不參與統(tǒng)計(jì)。

2.2 精度評(píng)定

本文以整個(gè)影像形變場(chǎng)中形變值的標(biāo)準(zhǔn)差作為精度評(píng)價(jià)指標(biāo)。由表3中的統(tǒng)計(jì)數(shù)據(jù)可知，在哨兵-2號(hào)10個(gè)影像對(duì)形變場(chǎng)中，無(wú)論是東西向還是南北向形變場(chǎng)，Band 8的標(biāo)準(zhǔn)差均表現(xiàn)為最小，其余波段獲取的影像形變場(chǎng)中形變值標(biāo)準(zhǔn)差的統(tǒng)計(jì)結(jié)果沒有明顯的差異。統(tǒng)計(jì)結(jié)果表明：就該試驗(yàn)區(qū)域而言，哨兵-2號(hào)4個(gè)10 m空間分辨率波段中，第8波段(Band 8)影像是地表形變監(jiān)測(cè)的最佳波段。另外，哨兵-2號(hào)4個(gè)可見光近紅外波段的光譜輻射性能中，第8波段的波段寬度最寬(115 nm)，信噪比值SNR最大(172)(見表1)，這也進(jìn)一步間接地驗(yàn)證了本文試驗(yàn)結(jié)論。

圖9 哨兵-2號(hào)影像對(duì)(2016-12-25—2017-01-04)東西向和南北向影像形變場(chǎng)中衛(wèi)星姿態(tài)角誤差去除前后對(duì)比Fig.9 Before and after removing the satellite attitude jitter distortion in the EW and NS components of image deformation field from the Sentinel-2 image pair (2016-12-25—2017-01-04)

標(biāo)準(zhǔn)差/m影像對(duì)編號(hào)東西向(E-W)南北向(N-S)Band 2Band 3Band 4Band 8Band 2Band 3Band 4Band 810.4170.390.3710.3340.4020.3770.3580.32620.4280.4080.4030.3640.410.3930.3890.35930.4370.420.4020.3690.4140.3990.3790.35740.4190.3750.3620.3310.4230.3950.3750.37450.4760.4440.4310.3920.4810.4610.4440.42260.4710.4350.4250.3990.4810.4650.4540.44670.4950.4650.4560.4210.4920.4710.4630.44180.4850.4450.4350.4010.4850.4630.4540.44190.4220.4060.3910.3760.4120.3980.3830.383100.4990.4750.470.4440.4880.4720.4620.447平均值0.4540.4260.4140.3830.4490.4290.4160.400

圖10 形變值標(biāo)準(zhǔn)差統(tǒng)計(jì)Fig.10 The displacement standard deviation diagram

3 新西蘭凱庫(kù)拉Mw7.8地震同震形變監(jiān)測(cè)

3.1 地震簡(jiǎn)介

2016年11月14日凌晨2時(shí)56分，新西蘭境內(nèi)凱庫(kù)拉鎮(zhèn)(Kaikoura)發(fā)生Mw7.8大地震，此次地震是該地區(qū)150多年來(lái)發(fā)生的強(qiáng)度最大的地震。根據(jù)USGS給出的震源機(jī)制解，該地震震中位于42.737°S,173.054°E，距離凱庫(kù)拉鎮(zhèn)西南方向約60 km，震源深度約為15 km；該地震是一個(gè)以右旋走滑為主兼少許逆沖分量的地震。該地震地表破裂始于Hope斷層，之后依次沿著Jordan thrust、Papatea及Kekerengu[34-36]等主要斷層向東北方向傳播。此次地震的地表破裂總長(zhǎng)超過150 km，沿Needles斷層的近海地表破裂長(zhǎng)度約為34 km(如圖1所示)。

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

本文利用以上數(shù)據(jù)處理方法和哨兵-2號(hào)影像獲取2016年11月14日Mw7.8新西蘭凱庫(kù)拉地震的東西向和南北向同震形變場(chǎng)，具體的數(shù)據(jù)處理流程如圖3所示。本文選用的哨兵-2號(hào)數(shù)據(jù)見表4，這兩景數(shù)據(jù)覆蓋了整個(gè)地震破裂帶區(qū)域(如圖2所示)，含云量小于1%，時(shí)間間隔約為1 a，太陽(yáng)照度和數(shù)據(jù)相干性得到了良好的保障，極大地減少了地形陰影誤差和時(shí)間失相關(guān)噪聲對(duì)形變結(jié)果的影響[15]。另外為了對(duì)地震形變結(jié)果進(jìn)行對(duì)比分析和交叉驗(yàn)證，本研究采用“小空間基線法”[16]選取了地震前后15 m空間分辨率的Landsat8全色影像(Band 8)各1景(見表4)，兩景影像采集于相同的季節(jié)，含云量較低，時(shí)間間隔也約為1 a。

表4同震形變監(jiān)測(cè)光學(xué)影像參數(shù)表

Tab.4Parametersoftheopticalimagesforthecoseismicdisplacementsmonitoring

衛(wèi)星波段影像采集時(shí)間含云量/(%)震前震后震前震后時(shí)間間隔/dSentinel-2Band 82016-02-192017-02-2300370Landsat8Band 82015-12-132016-12-154.766.76368

為了讓哨兵-2號(hào)和Landsat8數(shù)據(jù)能獲取相似分辨率的同震形變結(jié)果，本文的計(jì)算參數(shù)設(shè)置如下：搜索窗口的大小均設(shè)置為32×32像素(對(duì)應(yīng)地面分辨率：哨兵-2號(hào)為320×320 m，Landsat8為480×480 m)；移動(dòng)步長(zhǎng)：哨兵-2號(hào)設(shè)置為9×9像素，Landsat8設(shè)置為6×6像素(對(duì)應(yīng)地面分辨率均為90 m)；掩膜閾值均為0.9；迭代次數(shù)均為2次。通過以上相關(guān)性計(jì)算獲取的形變場(chǎng)中包含了各種系統(tǒng)誤差源(失相關(guān)噪聲、軌道誤差、條帶誤差、衛(wèi)星姿態(tài)角誤差等)，因此需要對(duì)上述同震形變結(jié)果進(jìn)行誤差后處理，本文采用第2章介紹的方法去除哨兵-2號(hào)和Landsat8形變場(chǎng)中的各種系統(tǒng)誤差。由于上述Landsat8形變場(chǎng)中沒有明顯的條帶誤差及衛(wèi)星姿態(tài)角誤差，因此本文只采用一次多項(xiàng)式曲面模型去除其形變場(chǎng)中的軌道誤差。最后，本文采用3×3像素大小的中值濾波窗口對(duì)哨兵-2號(hào)和Landsat8的結(jié)果進(jìn)行進(jìn)一步降噪。

3.3 結(jié)果比較分析

如圖11所示，此次地震的地表破裂帶主要發(fā)生在Kekerengu、Jordan thrust及Papatea斷層。其中Kekerengu斷層(右旋走滑)的地表形變并非對(duì)稱分布，該斷層上盤(NW)形變明顯大于下盤(SE)，這說明Kekerengu斷層向北方向傾斜，該斷層水平方向的最大水平滑移量為9～10 m。Jordan thrust斷層(右旋走滑)的最大水平滑移量約為3 m。位于Kekerengu斷層南部以及Jordan thrust斷層?xùn)|南方向存在一個(gè)與這兩個(gè)斷層相交的Papatea斷層，該斷層此前是一個(gè)未知的活躍斷層，該斷層滑動(dòng)兼具左旋走滑(最大水平滑移量約4 m)，位于Clarence河西南方向約20 km處。此外，在Marlborough斷層系統(tǒng)中，沿西南方向的主要斷層(如Jordan thrust等)的地表破裂長(zhǎng)度以及最大水平滑移量均明顯的小于東北方向的主要斷層(如Kekerengu斷層等)。如圖11中近斷裂帶區(qū)域A,D(哨兵-2號(hào))及區(qū)域B,E(Landsat8)對(duì)比所示，與Landsat8相比，整個(gè)哨兵-2號(hào)影像形變場(chǎng)中異常值(Nan值)更少，在近場(chǎng)區(qū)域尤為明顯。圖11(c)、(f)分別為哨兵-2號(hào)與Landsat8的東西向和南北向形變場(chǎng)差值圖，其形變值的均值在近場(chǎng)(區(qū)域C、F)和遠(yuǎn)場(chǎng)區(qū)域(近場(chǎng)以外的形變區(qū)域)均較為接近。其中近場(chǎng)區(qū)域東西向形變值差值的均值約為0.1 m，南北向約為-0.2 m；而遠(yuǎn)場(chǎng)區(qū)域東西向形變值差值的均值約為0.2 m，南北向約為-0.1 m。但是形變值差值的均方差(RMSE)在近場(chǎng)和遠(yuǎn)場(chǎng)區(qū)域有一定的差異，其中近場(chǎng)區(qū)域的RMSE約為0.6 m，而遠(yuǎn)場(chǎng)區(qū)域的RMSE約為1.1 m。這種差異可能是因?yàn)榻鼒?chǎng)地形起伏較小、植被覆蓋較少，從而使得光學(xué)影像形變監(jiān)測(cè)信噪比較高，而遠(yuǎn)場(chǎng)山區(qū)地形更加復(fù)雜，植被也更為茂盛，所以信噪比較低。根據(jù)COSI-Corr軟件基于亞像素的相關(guān)性匹配算法用來(lái)探測(cè)同震形變(近場(chǎng))的匹配精度一般可達(dá)1/20個(gè)像素[7]，即理論上哨兵-2號(hào)和Landsat8的同震形變監(jiān)測(cè)精度最高分別可以達(dá)到0.5 m和0.75 m，這表明本文所計(jì)算的精度水平(RMSE)較為合理。

總體來(lái)說，在近場(chǎng)區(qū)域哨兵-2號(hào)影像的同震形變監(jiān)測(cè)結(jié)果與Landsat8影像的形變結(jié)果有較好的一致性，但是哨兵-2號(hào)的形變圖更加清晰且形變結(jié)果中的異常值更少。如圖12所示，為哨兵-2號(hào)與Landsat8地震形變結(jié)果的中誤差，結(jié)果同樣表明哨兵-2號(hào)同震形變結(jié)果的精度明顯優(yōu)于Landsat8影像，這可能與哨兵-2號(hào)數(shù)據(jù)具有更高的時(shí)空分辨率密切相關(guān)。

4 結(jié) 論

本文以新西蘭地震覆蓋區(qū)域?yàn)槔?，?duì)哨兵-2號(hào)影像形變場(chǎng)中各種系統(tǒng)誤差源的時(shí)空分布特征及誤差去除方法進(jìn)行了系統(tǒng)分析，并且創(chuàng)新地提出了改進(jìn)的均值相減法對(duì)形變場(chǎng)中的衛(wèi)星姿態(tài)角抖動(dòng)誤差進(jìn)行有效去除。在誤差分析的基礎(chǔ)上，還對(duì)哨兵-2號(hào)影像4個(gè)10 m空間分辨率的波段(Band 2/3/4/8)獲取地表形變信息的能力進(jìn)行了統(tǒng)計(jì)分析及精度評(píng)定，統(tǒng)計(jì)結(jié)果揭示了Band8為地表形變監(jiān)測(cè)的最佳波段。最后利用時(shí)間間隔均約為1年的哨兵-2號(hào)和Landsat8影像對(duì)，獲取了2016年11月14日MW7.8新西蘭凱庫(kù)拉地震的同震形變場(chǎng)。形變監(jiān)測(cè)結(jié)果表明：地表滑動(dòng)主要集中在Kekerengu、Papatea及Jordan Thrust斷層，其中沿Kekerengu斷層的最大水平滑移量為9～10 m，沿Papatea斷層的最大水平滑移量為3～4 m，沿Jordan Thrust斷層的最大水平滑移量約為3 m。另外，兩種光學(xué)數(shù)據(jù)的形變結(jié)果比較表明：雖然在近場(chǎng)哨兵-2號(hào)與Landsat8的同震形變結(jié)果一致性較好，但是在整個(gè)影像中，與Landsat8相比，哨兵-2號(hào)獲取的形變結(jié)果異常值更少，形變圖更清晰，中誤差更小。因此，哨兵-2號(hào)數(shù)據(jù)可以很好地彌補(bǔ)其他數(shù)據(jù)在空間分辨率和時(shí)間分辨率上的缺口，為防震減災(zāi)工作的開展提供一個(gè)很好的數(shù)據(jù)平臺(tái)。

致謝:本研究所用的哨兵-2號(hào)和Landsat8數(shù)據(jù)來(lái)源于USGS(glovis.usgs.gov/)，文中主要采用GMT5.2.1軟件繪制各種圖形，采用美國(guó)加州理工學(xué)院研發(fā)的軟件包COSI-Corr(http:∥www.tectonics.caltech.edu)進(jìn)行數(shù)據(jù)處理。

圖11 哨兵-2號(hào)影像對(duì)(2016-02-19—2017-02-23)和Landsat8影像對(duì)(2015-12-13—2016-12-15)分別獲取的2016年11月14日Mw7.8新西蘭凱庫(kù)拉地震的同震形變場(chǎng)Fig.11 Coseismic deformation fields obtained by Sentinel-2 (2016-02-19—2017-02-23) and Landsat8 (2015-12-13—2016-12-15)

圖12 新西蘭地震形變場(chǎng)形變值中誤差統(tǒng)計(jì)Fig.12 The displacement standard deviation of the New Zealand earthquake deformation field