宋艷若，莊啟智，陳輝，4，朱大明，程亮，4，5

(1.昆明理工大學(xué) 國土資源工程學(xué)院，昆明 650093；2.南京大學(xué) 地理與海洋科學(xué)學(xué)院，南京 210023；3.南京大學(xué) 江蘇省地理信息技術(shù)重點實驗室，南京 210023；4.南京大學(xué) 中國南海研究協(xié)同創(chuàng)新中心，南京 210023；5.南京大學(xué) 自然資源部國土衛(wèi)星遙感應(yīng)用重點實驗室，南京 210023)

0 引言

遙感變化檢測技術(shù)為全面、快速掌握城市水資源演化趨勢提供了有效手段，對維護城市水域生態(tài)格局和可持續(xù)發(fā)展具有重要的現(xiàn)實意義[1-2]。對于中分辨率影像，遙感時序變化檢測可以克服雙時相變化檢測存在的高質(zhì)量影像要求、低影像誤差容錯率、缺乏時間演化信息利用等問題[3]，弱化“同物異譜”或“同譜異物”現(xiàn)象造成的影響，擁有較高的檢測精度和效率，被廣泛應(yīng)用于城市擴張、土地覆蓋、異常災(zāi)害等各個領(lǐng)域的動態(tài)監(jiān)測[4-5]。

近年來，遙感時序變化檢測方法取得了長足發(fā)展，但依然面臨著眾多挑戰(zhàn)。首先，時間序列降噪重建研究有待加強。時序分析方法的側(cè)重點不同，對應(yīng)的預(yù)處理需求隨之不同，降噪方法的選擇要充分考慮時序數(shù)據(jù)的特點?，F(xiàn)今時間序列的重建方法大多基于歸一化水體指數(shù)(normalized differnce vegetation index，NDVI)、增強型植被指數(shù)(enhanced vegetation index，EVI)等植被指數(shù)[6-7]展開研究，針對歸一化水體指數(shù)(normalized difference water index，NDWI)時間序列降噪重建研究較少。其次，當(dāng)前遙感時序變化檢測模型對影像時空譜信息的挖掘與利用不夠充分。現(xiàn)已有一些成熟的遙感時序變化檢測模型廣泛應(yīng)用于各個領(lǐng)域。如Verbesselt等[8-9]提出的BFAST(breaks for additive season and trend)算法減弱了噪聲和季節(jié)性趨勢對檢測結(jié)果的影響，檢測出小范圍的森林?jǐn)_動。Zhu等[10]提出了連續(xù)變化檢測和分類算法(continuous change detection and classification，CCDC)，并以時序模型參數(shù)作為特征進行土地覆蓋分類。但這些模型均以像元為分析單位，忽略了像元空間和上下文特征。如何兼顧影像的時空譜特征來提升檢測精度已然成為研究的熱點問題，一些學(xué)者已經(jīng)開展了一些研究。Lin等[11]結(jié)合Landsat數(shù)據(jù)的時間特征識別斷點，利用鄰域信息提升檢測精度。Jualea等[12]利用時空數(shù)據(jù)挖掘技術(shù)，結(jié)合空間分布的一致性，提取植物生長信息和土地覆蓋變化信息。但這些方法多應(yīng)用于土地覆蓋領(lǐng)域，且處于基礎(chǔ)方法研究階段，在水資源檢測領(lǐng)域研究較少。

根據(jù)以上分析，本文基于對城市水體變化規(guī)律的理解與表達，提出一種顧及時空譜特征遙感時序變化檢測方法。該方法結(jié)合水體光譜特征建立遙感時間序列，并探討多種降噪重建算法的適用性；依據(jù)變化像元時序在時間維的獨有特性確定疑似變化像元子序列；并利用像元空間鄰近性原則，融合建筑物指數(shù)[13]，獲取最終檢測結(jié)果。本文以杭州錢塘江下沙段為典型實驗區(qū)，利用Sentinel-2影像驗證該方法的穩(wěn)定性和有效性。

1 研究方法

本文提出的顧及影像時空譜規(guī)律的遙感時序變化檢測方法技術(shù)流程如圖1所示。首先，計算預(yù)處理后Sentinel-2影像集的NDWI指數(shù)，建立像元級時間序列，分別采用傅里葉變換、小波變換、中值濾波、S-G濾波進行降噪處理，并對重建結(jié)果評價選優(yōu)；其次，采用動態(tài)時間規(guī)整算法衡量待檢測像元時間序列與樣本像元時間序列的相似度，認(rèn)定結(jié)果小于閾值的像元為疑似變化像元；最后，通過空間相關(guān)性分析濾除偽變化像元，并將最終檢測結(jié)果二值分類，評價分析。

圖1 顧及時空譜規(guī)律的遙感時序變化檢測方法

1.1 遙感時間序列建模

1)遙感時間序列構(gòu)建。遙感時間序列是由影像數(shù)據(jù)集和其生成時間構(gòu)成的有序元素集合[14]。本研究將大氣校正、影像配準(zhǔn)等預(yù)處理后的數(shù)據(jù)集進行歸一化處理，生成對應(yīng)的NDWI指數(shù)[15]，并以像元為單位建立遙感時間序列。構(gòu)建方法如下。

假設(shè)有p×q行列的m景影像，遍歷每景影像中所有像元的NDWI值，依據(jù)各影像觀測時間先后順序?qū)ν晃恢孟裨r間序列，每個像元均可構(gòu)建一個Wi={(w)i，t1)，(w)i，t2)，(w)i，t3)，…，(w)i，tm}的子序列，其中wi表示像元i點的NDWI指數(shù)值，tj表示第j景影像的觀測時間，最終m景影像構(gòu)建得到的時間序列可表示為N={((x1，y1)，W1)，((x1，y2)，W2)，…，((xp，yq)，Wp×q)}，其中(xa，yb)為像元的行列坐標(biāo)。

2)時間序列降噪重建。為平衡遙感時間序列在時間分辨率和影像質(zhì)量上的沖突，削弱大氣、云層或其他噪聲導(dǎo)致的誤差，對遙感時間序列降噪重建處理。不同時間序列重建算法對不同噪聲點的去除和趨勢完整性的保持度不同，本文從降噪算法的時頻特性出發(fā)，選取中值濾波[16]、S-G濾波[17]、傅里葉變換[18]、小波變換[19]4種常用的降噪重建方法，通過對比分析在NDWI時間序列中的適用性，選擇最優(yōu)重建結(jié)果進行下一步分析。

1.2 時間相似度分析

利用遙感時序數(shù)據(jù)時間維特征識別變化像元的關(guān)鍵為相似性分析，即計算待檢測像元時間序列Vi與變化樣本像元時間序列Vs的相似度TSis，存在一個閾值α(α>0)，當(dāng)TSis>α?xí)r認(rèn)定時序Vi和Vs是相似的，反之亦然。本研究通過動態(tài)時間規(guī)整[20](dynamic time warping，DTW)算法衡量時序數(shù)據(jù)間的相似度。

1.3 空間相關(guān)性分析

地理學(xué)第一定律[21]揭示了地物空間分布的鄰近性原則，變化檢測中則可體現(xiàn)為變化像元孤立存在時極大可能存在假陽性。遙感時間序列的構(gòu)建是以像元為單位獨立建模，彼此之間互不相關(guān)，忽略了像元間存在的空間和上下文特征，造成很多異常像元的誤提取。同時，為弱化隨建筑物類型改變導(dǎo)致光譜反射率改變的偽變化區(qū)域帶來的影響，本文在時間相似度分析的基礎(chǔ)上，融合歸一化差值建筑用地指數(shù)(normalized difference building index，NDBI)，制定像元空間鄰域分布規(guī)則，濾除偽變化像元控制誤報率。如圖2所示，首先確定疑似變化像元Mi，j位置的NDBI值是否為0，若“是”則進行下一步分析，反之濾除；進而判斷Mi，j像元周圍的Mi-1，j、Mi，j-1、Mi+1，j、Mi，j+1像元中是否存在一個及以上可能發(fā)生變化的像元，若“是”則認(rèn)定為疑似變化像元，并以3×3的窗口遍歷。若窗口內(nèi)除中心位置外的8個像元中存在3個及以上的疑似變化像元，認(rèn)定Mi，j為變化像元，反之是偽變化像元。圖2中綠色斜線填充方框為疑似變化像元，空白方框表示未變化像元。

圖2 空間相關(guān)性分析原理圖

1.4 精度評價指標(biāo)

1)降噪重建評價指標(biāo)。為了定量對比時間序列降噪重建效果，本文選用信噪比(signal-to-noise ratio，SNR)、均方根誤差(root mean square error，RMSE)、互相關(guān)系數(shù)(cross-correlation coefficient，R)等評價指標(biāo)，評估各方法重建時序前后的噪聲去除度、保真性和相似性。其中，RMSE值越小，表示時序保真性越強，SNR、R值越大，表明時間序列重建效果越好。

2)檢測結(jié)果精度評價指標(biāo)。本文利用混淆矩陣和Kappa系數(shù)評價變化檢測結(jié)果精度，將各方法的檢測結(jié)果二值化并與真值對比生成混淆矩陣，計算總體精度(over accuracy，OA)、查全率(true positive rate，TPR)、虛檢率(false detection rate，F(xiàn)DR)和Kappa系數(shù)等指標(biāo)值。其中OA、TPR、Kappa值越大，表明變化檢測結(jié)果越好，F(xiàn)DR值越小，表明檢測結(jié)果準(zhǔn)確性越高。

2 結(jié)果與分析

2.1 研究區(qū)及數(shù)據(jù)

本研究利用Sentinel-2衛(wèi)星遙感影像，以錢塘江下沙段為研究區(qū)開展實驗，研究區(qū)概況及數(shù)據(jù)集如圖3所示。研究區(qū)江北岸建筑密集，江南岸有裸地、耕地、水田等分布，水體提取干擾因素多、土地利用類型復(fù)雜、特征變化明顯，具有較強代表性。數(shù)據(jù)源為2016年12月30日至2020年12月24日的Sentinel-2衛(wèi)星遙感影像，為保證研究區(qū)影像質(zhì)量，篩選得到研究區(qū)波段完整、云量小于30%、數(shù)據(jù)無缺失的40景影像。數(shù)據(jù)預(yù)處理包括大氣校正、影像配準(zhǔn)、影像融合及裁剪等，最終生成930像素×790像素的研究區(qū)影像集，預(yù)處理操作在SNAP軟件、Sen2Cor_v2.8.0插件、ENVI 5.3中完成。

圖3 研究區(qū)內(nèi)Sentinel-2影像集

2.2 降噪重建結(jié)果與分析

時間序列降噪重建精度對于參數(shù)設(shè)置具有較高的敏感性，本研究將原始時序像元劃分為非水體像元、純凈水體像元及變化像元3種類型，根據(jù)NDWI時序數(shù)據(jù)源特點及不同類型像元趨勢規(guī)律，基于不同參數(shù)設(shè)置的多降噪方法進行了大量對比預(yù)實驗，建立參數(shù)設(shè)置規(guī)則分別對時間序列降噪重建。具體規(guī)則為：中值濾波的滑動窗口大小為5；S-G濾波采用的滑動窗口及多項式分別為7、3；傅里葉變換的截斷點數(shù)設(shè)置為7；小波變換采用heursure閾值進行2層分解，并選用db4小波為小波基。

圖4、圖5、圖6分別表示4種方法對非水體像元、變化像元及純凈水體像元時間序列重建前后的效果對比圖。由圖可知，經(jīng)重建后的像元時間序列與原始像元時間序列相比有明顯平滑效果，異常點均可被去除，且重建后曲線仍可保持原有特征和獨有變化趨勢。由圖5可知，4種降噪方法對變化像元時間序列重建結(jié)果基本一致，并可以從時間維度中體現(xiàn)由水體變?yōu)榉撬w的演化趨勢。從圖4、圖6中可以看出，非變化像元時間序列重建后整體值的波動區(qū)間變得更為平緩，突顯出非水體像元值基本小于0、水體像元值基本大于0的特征。對比不同方法的重建效果，傅里葉變換的平滑度最好，但容易造成“高值低谷、低值高估”的情況，小波變換略優(yōu)于S-G濾波，中值濾波平滑度最差；就保持整體趨勢而言，S-G濾波更有優(yōu)勢，中值濾波的效果相對較差。

表1反映不同降噪算法對所有像元時間序列重建前后SNR、RMSE和R值的平均狀況。由表1可知，S-G濾波的SNR和R統(tǒng)計值均遠大于其他方法，RMSE值僅次于小波變換，表明該方法降噪重建前后時間序列相關(guān)性較高，整體特征保持性較好。中值濾波的RMSE值最大、R值最小，相較于其他方法保真性最差。多指標(biāo)綜合評價，S-G濾波效果最優(yōu)，故本研究采用S-G濾波重建后的時間序列進行時空分析。

圖4 非水體像元時間序列重建效果對比圖

圖5 變化像元時間序列重建效果對比圖

圖6 純凈水體像元時間序列重建效果對比圖

表1 時間序列重建效果精度統(tǒng)計

2.3 變化檢測結(jié)果與分析

為充分驗證本文方法的有效性，將本文方法檢測結(jié)果與采用閾值分割法[22]、支持向量機(support vector machine，SVM)法的雙時相變化檢測結(jié)果，以及采用CCDC算法的遙感時序變化檢測結(jié)果進行對比分析。不同方法的變化檢測結(jié)果如圖7所示，其中圖7(a)為真實變化區(qū)域的二值分類結(jié)果圖，白色表示變化區(qū)域，由在實地調(diào)研的基礎(chǔ)上參照高分辨率影像等資料人工勾繪得出。

圖7(b)和圖7(c)分別是基于雙時相影像，采用閾值分割法和SVM分類法得到的變化檢測結(jié)果，其基本思想為：將雙時相影像通過計算NDWI值進行閾值分割，或采用SVM分類法分類，提取影像水域，并進行數(shù)學(xué)運算得到最終變化檢測結(jié)果。在復(fù)雜的城區(qū)背景，河流中包含更多的雜質(zhì)，且道路、建筑物及陰影等廣泛存在，這些因素都會對檢測結(jié)果精度產(chǎn)生直接影響，特別在閾值分割法中體現(xiàn)得尤為明顯。由于建筑物陰影、道路等光譜特性與水體相近，且含沙量較多的水體區(qū)域易與非水體區(qū)域混淆，僅從水體指數(shù)中難以區(qū)分，從而導(dǎo)致閾值分割法的檢測結(jié)果中存在大量的建筑物陰影、道路及水體含沙量較高等區(qū)域的誤識別，造成較多的虛檢。SVM分類法在訓(xùn)練樣本時獲取了先驗知識并運用于水體分類中，可以弱化城區(qū)人造物和水體雜質(zhì)的影響，一定程度上降低虛檢現(xiàn)象，但漏檢現(xiàn)象卻較為明顯，原因是在分類過程中易將大部分混合像元判定為非水體區(qū)域，最終歸為未變化像元，導(dǎo)致漏檢現(xiàn)象。

圖7(d)為采用CCDC算法得到的遙感時序變化檢測結(jié)果，該方法基于像元NDWI值建立時間序列，通過比較像元時間序列的預(yù)測值與觀測值偏離程度來衡量是否發(fā)生變化。CCDC算法從時間演化趨勢出發(fā)深入挖掘遙感時序的時譜信息，可以有效減少由水體含沙量季節(jié)性變化引起的噪聲干擾，或因“同譜異物”現(xiàn)象造成的誤差，能夠有效識別橋梁、道路等細小地物未發(fā)生變化區(qū)域，相較于圖7(b)、圖7(c)的檢測效果有明顯提升，特別是橋梁、道路及河岸下方等區(qū)域的虛檢現(xiàn)象明顯降低。但該方法忽略了變化像元的空間和上下文特征，導(dǎo)致由建筑物類型改變引起光譜反射率改變的偽變化區(qū)域分布零散，破碎度較高的現(xiàn)象。本文在充分利用像元時間序列時譜信息的基礎(chǔ)上，融合NDBI指數(shù)，兼顧像元的空間特征，提出的遙感時序變化檢測方法的檢測結(jié)果如圖7(e)所示，相較于前3種方法檢測效果最佳，最接近真實變化參考圖。

圖7 不同變化檢測方法的變化檢測結(jié)果

表2為不同變化檢測方法定量對比的評價結(jié)果。由表2可知，基于遙感時間序列進行變化檢測的方法優(yōu)勢明顯，Kappa系數(shù)和OA值均高于雙時相變化檢測方法。本文提出方法的OA值和Kappa系數(shù)是最高的，查全率僅次于CCDC算法，TPR值降低了4%，但FDR值降低了23.45%，這是由于進行空間相關(guān)性分析有效控制誤報率的同時存在小部分誤刪現(xiàn)象。對比簡單閾值法和SVM分類法的雙時相變化檢測，TPR分別提升了3.89%、22.47%，F(xiàn)DR值分別降低了43.51%、23.37%，Kappa系數(shù)提高了0.332 7和0.233 2。綜合TPR、FDR、OA和Kappa系數(shù)等指標(biāo)分析，本文所提出變化檢測方法效果最優(yōu)。

表2 不同變化檢測方法精度統(tǒng)計

3 結(jié)束語

針對傳統(tǒng)遙感時序變化檢測未能考慮像元間的空間鄰近特征，對時空信息挖掘不夠深入等問題，本文提出一種顧及時空譜特征的遙感時序城市水體變化檢測方法。并以Sentinel-2影像為數(shù)據(jù)源，以杭州錢塘江下沙段為典型研究區(qū)開展實驗。結(jié)果表明，在時間序列降噪重建的適用性分析中，S-G濾波對NDWI時間序列的降噪重建效果最佳；相較于雙時相變化檢測，遙感時序變化檢測擁有更高的檢測精度；本文方法精度最高，總體精度達到99%，Kappa系數(shù)達到0.808 6，虛檢率降低了20%～50%，已在實驗區(qū)研究階段論證了其有效性，但仍需加強整個技術(shù)流程的自動化實現(xiàn)，自動設(shè)定閾值是下一步研究的重點。