哨兵-2號光學影像地表形變監測：以2016年MW7.8新西蘭凱庫拉地震為例

賀禮家，馮光財，馮志雄，高 華

1. 中南大學地球科學與信息物理學院雷達遙感研究室，湖南 長沙 410083； 2. 東華理工大學測繪工程學院，江西 南昌 330013

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

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

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

1 系統誤差改正

1.1 數據介紹

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

表1 哨兵-2號光學影像光譜波段信息簡介

1.2 數據處理

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

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

圖2 哨兵-2號相關影像對采集日期分布Fig.2 Acquisition dates of the Sentinel-2 correlation image pairs

影像編號采集時間太陽天頂角/(°)太陽方位角/(°)含云量/(%)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

光學影像地表形變監測完整的數據處理流程(如圖3所示)，主要包括如下4個步驟：

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

(2) 相關性計算：以COSI-Corr為數據處理平臺，對影像對進行亞像素的相關性匹配(sub-pixel correlation)計算[23-24]。對于哨兵-2號數據，參數設置如下：搜索窗口大小為32×32像素(地面分辨率約為320×320 m)；移動步長為8×8像素(地面分辨率為80×80 m)；掩膜閾值為0.9；迭代次數為2次。

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

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

圖3 光學影像地表形變監測數據處理流程Fig.3 Flow chart of processing the ground deformation obtained from optical images

1.3 誤差分析及糾正

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

注：影像形變場處于WGS-84坐標系，分別以朝東和朝北方向為形變值的正方向。圖4 哨兵-2號影像對(2016-12-15—2016-12-25)東西向和南北向影像形變場系統誤差源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 失相關噪聲

失相關噪聲與地表輻射性能息息相關。當地表輻射性能變化較小時，可能會使影像的紋理特征產生微小變化；而當輻射性能變化較大時，則會引起影像紋理的缺失。失相關噪聲水平可用信噪比值(SNR)進行描述。當影像中某個區域的信噪比值普遍較低時，該區域可能受失相關噪聲的影響[9,16]。失相關噪聲產生的原因可歸納為如下幾點：

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

(2) 地形陰影。同一傳感器在不同季節采集的影像對應的太陽角(包括太陽高度角和方位角)不同會使形變場中某些區域產生地形陰影。

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

(1) 程序自動掩膜掉形變場中SNR低于某閾值(經驗值取0.9)的區域。

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

1.3.2 軌道誤差

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

φorbit=a0+a1x+a2y+a3xy

(1)

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

圖5 哨兵-2號影像對(2016-12-25—2017-01-04)東西向和南北向影像形變場中軌道誤差去除前后對比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 條帶誤差

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

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

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

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

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

圖7 哨兵-2號影像對(2016-12-25—2017-01-04)東西向和南北向影像形變場中條帶誤差去除前后對比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 衛星姿態角誤差

在形變場中，沿交叉軌道方向呈現出周期性均勻平行分布的帶狀信號，即衛星姿態角抖動誤差。該誤差產生的原因可歸納如下：

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

(2) 在影像正射校正過程中，由于衛星姿態角欠采樣，導致無法精確計算衛星姿態角[32]。

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

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

(2)

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

圖8 改進后的均值相減法去除衛星姿態角誤差數據處理流程Fig.8 Flow chart of removing the satellite attitude jitter distortion using the modified mean subtracting method

2 最佳波段分析

2.1 數據處理

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

(1) 相關性計算：利用COSI-Corr分別計算不同波段(Band 2/3/4/8)對應的10個相關性影像對(如圖2所示)的二維地表形變場。計算參數統一為：搜索窗口大小為32×32像素；移動步長為8×8像素；掩膜閾值為0.9；迭代次數為2次。

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

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

2.2 精度評定

本文以整個影像形變場中形變值的標準差作為精度評價指標。由表3中的統計數據可知，在哨兵-2號10個影像對形變場中，無論是東西向還是南北向形變場，Band 8的標準差均表現為最小，其余波段獲取的影像形變場中形變值標準差的統計結果沒有明顯的差異。統計結果表明：就該試驗區域而言，哨兵-2號4個10 m空間分辨率波段中，第8波段(Band 8)影像是地表形變監測的最佳波段。另外，哨兵-2號4個可見光近紅外波段的光譜輻射性能中，第8波段的波段寬度最寬(115 nm)，信噪比值SNR最大(172)(見表1)，這也進一步間接地驗證了本文試驗結論。

圖9 哨兵-2號影像對(2016-12-25—2017-01-04)東西向和南北向影像形變場中衛星姿態角誤差去除前后對比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)

標準差/m影像對編號東西向(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 形變值標準差統計Fig.10 The displacement standard deviation diagram

3 新西蘭凱庫拉Mw7.8地震同震形變監測

3.1 地震簡介

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

3.2 數據處理

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

表4同震形變監測光學影像參數表

Tab.4Parametersoftheopticalimagesforthecoseismicdisplacementsmonitoring

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

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

3.3 結果比較分析

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

總體來說，在近場區域哨兵-2號影像的同震形變監測結果與Landsat8影像的形變結果有較好的一致性，但是哨兵-2號的形變圖更加清晰且形變結果中的異常值更少。如圖12所示，為哨兵-2號與Landsat8地震形變結果的中誤差，結果同樣表明哨兵-2號同震形變結果的精度明顯優于Landsat8影像，這可能與哨兵-2號數據具有更高的時空分辨率密切相關。

4 結 論

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

致謝:本研究所用的哨兵-2號和Landsat8數據來源于USGS(glovis.usgs.gov/)，文中主要采用GMT5.2.1軟件繪制各種圖形，采用美國加州理工學院研發的軟件包COSI-Corr(http:∥www.tectonics.caltech.edu)進行數據處理。

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

圖12 新西蘭地震形變場形變值中誤差統計Fig.12 The displacement standard deviation of the New Zealand earthquake deformation field