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

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

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

哨兵-2號(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)測研究，如SPOT[3-9]、ASTER[10]、Landsat[11-15]以及航空影像[16-17]等。與ASTER、Landsat系列等傳統(tǒng)的光學(xué)數(shù)據(jù)相比，哨兵-2號衛(wèi)星星座具有中高空間分辨率(最高10 m)、較短重訪周期(5 d)及數(shù)據(jù)免費(fèi)等優(yōu)勢，因此在地表形變監(jiān)測中的應(yīng)用潛力將更加明顯[18]。由于目前距哨兵-2號升空時(shí)間較短，其數(shù)據(jù)積累較少，所以目前關(guān)于它的形變監(jiān)測系統(tǒng)誤差改正和應(yīng)用研究還較少，還尚未有學(xué)者對哨兵-2號影像中存在的各種系統(tǒng)誤差源和改正方法進(jìn)行系統(tǒng)的研究。此外，該衛(wèi)星星座可以獲取4個(gè)10 m空間分辨率的波段影像，何種波段更適合地表形變監(jiān)測以及哨兵-2號與Landsat8兩種光學(xué)遙感影像地表形變監(jiān)測精度之間的差異等問題，目前還尚未有分析比較。以上內(nèi)容對于哨兵-2號衛(wèi)星數(shù)據(jù)的應(yīng)用和推廣都具有重要意義。

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

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

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

哨兵-2號多光譜成像儀(MSI)采用推掃式成像模式，該衛(wèi)星星座含13個(gè)通道，其中10個(gè)通道是波長為433～955 nm的可見光近紅外譜段(VNIR)，其他3個(gè)通道是波長為1360～2280 nm的短波紅外譜段(SWIR)。哨兵-2號衛(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號光學(xué)影像光譜波段信息簡介

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

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

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

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

影像編號采集時(shí)間太陽天頂角/(°)太陽方位角/(°)含云量/(%)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)測完整的數(shù)據(jù)處理流程(如圖3所示)，主要包括如下4個(gè)步驟：

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

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

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

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

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

1.3 誤差分析及糾正

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

注：影像形變場處于WGS-84坐標(biāo)系，分別以朝東和朝北方向?yàn)樾巫冎档恼较颉D4 哨兵-2號影像對(2016-12-15—2016-12-25)東西向和南北向影像形變場系統(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é)采集的影像對應(yīng)的太陽角(包括太陽高度角和方位角)不同會(huì)使形變場中某些區(qū)域產(chǎn)生地形陰影。

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

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

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

1.3.2 軌道誤差

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

φorbit=a0+a1x+a2y+a3xy

(1)

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

圖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 條帶誤差

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

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

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

為了消除CCD線陣列引起的條帶誤差，本文采用傳統(tǒng)的“均值相減法”去除條帶誤差[1,16]。本文選用哨兵-2號2016-12-25—2017-01-04震后影像對進(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號影像對(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 衛(wèi)星姿態(tài)角誤差

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

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

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

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

為了去除形變場中衛(wèi)星姿態(tài)角抖動(dòng)引起的誤差，傳統(tǒng)采用與去除條帶誤差相同的方法(均值相減法)。但是本文考慮到沿衛(wèi)星交叉軌道方向地形起伏變化、植被覆蓋等因素的影響，沿衛(wèi)星交叉軌道方向的衛(wèi)星姿態(tài)角誤差是變化而非固定不變的，直接采用傳統(tǒng)方法并不能取得良好的誤差去除效果。因此,本文結(jié)合哨兵-2號衛(wèi)星推掃式成像儀焦平面上共有12個(gè)交錯(cuò)排列的CCD線陣列探測器的特點(diǎn)，提出改進(jìn)后的“均值相減法”去除衛(wèi)星姿態(tài)角誤差。該改進(jìn)方法的主要思想為：首先定義一個(gè)量η用于定量的衡量形變場中的衛(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號2016-12-25—2017-01-04影像對為例，利用改進(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號光學(xué)影像中共有4個(gè)10 m空間分辨率的可見光近紅外波段(Band 2/3/4/8)。雖然這4個(gè)光譜波段具有相同的空間分辨率，但是其波段寬度、SNR值以及輻射特性卻不盡相同[33]，詳見表1。本文利用10對哨兵-2號震后影像對(如圖2所示)對應(yīng)的4波段影像數(shù)據(jù)進(jìn)行試驗(yàn)分析。由于這些震后影像對的時(shí)間基線較短(10 d)，地表一般不會(huì)發(fā)生明顯的形變。因此，通過分析不同波段影像獲取的影像形變場中形變值的誤差水平，可以反映出不同波段影像地表形變監(jiān)測能力的差異。本文利用4個(gè)不同波段的影像對組合進(jìn)行地表形變監(jiān)測的最佳波段選取分析。具體的數(shù)據(jù)處理流程主要包括以下2個(gè)步驟：

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

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

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

2.2 精度評定

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

圖9 哨兵-2號影像對(2016-12-25—2017-01-04)東西向和南北向影像形變場中衛(wèi)星姿態(tài)角誤差去除前后對比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影像對編號東西向(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 新西蘭凱庫拉Mw7.8地震同震形變監(jiān)測

3.1 地震簡介

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

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

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

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

Tab.4Parametersoftheopticalimagesforthecoseismicdisplacementsmonitoring

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

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

3.3 結(jié)果比較分析

如圖11所示，此次地震的地表破裂帶主要發(fā)生在Kekerengu、Jordan thrust及Papatea斷層。其中Kekerengu斷層(右旋走滑)的地表形變并非對稱分布，該斷層上盤(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等)的地表破裂長度以及最大水平滑移量均明顯的小于東北方向的主要斷層(如Kekerengu斷層等)。如圖11中近斷裂帶區(qū)域A,D(哨兵-2號)及區(qū)域B,E(Landsat8)對比所示，與Landsat8相比，整個(gè)哨兵-2號影像形變場中異常值(Nan值)更少，在近場區(qū)域尤為明顯。圖11(c)、(f)分別為哨兵-2號與Landsat8的東西向和南北向形變場差值圖，其形變值的均值在近場(區(qū)域C、F)和遠(yuǎn)場區(qū)域(近場以外的形變區(qū)域)均較為接近。其中近場區(qū)域東西向形變值差值的均值約為0.1 m，南北向約為-0.2 m；而遠(yuǎn)場區(qū)域東西向形變值差值的均值約為0.2 m，南北向約為-0.1 m。但是形變值差值的均方差(RMSE)在近場和遠(yuǎn)場區(qū)域有一定的差異，其中近場區(qū)域的RMSE約為0.6 m，而遠(yuǎn)場區(qū)域的RMSE約為1.1 m。這種差異可能是因?yàn)榻鼒龅匦纹鸱^小、植被覆蓋較少，從而使得光學(xué)影像形變監(jiān)測信噪比較高，而遠(yuǎn)場山區(qū)地形更加復(fù)雜，植被也更為茂盛，所以信噪比較低。根據(jù)COSI-Corr軟件基于亞像素的相關(guān)性匹配算法用來探測同震形變(近場)的匹配精度一般可達(dá)1/20個(gè)像素[7]，即理論上哨兵-2號和Landsat8的同震形變監(jiān)測精度最高分別可以達(dá)到0.5 m和0.75 m，這表明本文所計(jì)算的精度水平(RMSE)較為合理。

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

4 結(jié) 論

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

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

圖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 新西蘭地震形變場形變值中誤差統(tǒng)計(jì)Fig.12 The displacement standard deviation of the New Zealand earthquake deformation field