顧及時空譜規律的城市水體時序變化檢測方法

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

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

0 引言

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

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

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

1 研究方法

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

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

1.1 遙感時間序列建模

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

假設有p×q行列的m景影像，遍歷每景影像中所有像元的NDWI值，依據各影像觀測時間先后順序對同一位置像元建立時間序列，每個像元均可構建一個Wi={(w)i，t1)，(w)i，t2)，(w)i，t3)，…，(w)i，tm}的子序列，其中wi表示像元i點的NDWI指數值，tj表示第j景影像的觀測時間，最終m景影像構建得到的時間序列可表示為N={((x1，y1)，W1)，((x1，y2)，W2)，…，((xp，yq)，Wp×q)}，其中(xa，yb)為像元的行列坐標。

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

1.2 時間相似度分析

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

1.3 空間相關性分析

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

圖2 空間相關性分析原理圖

1.4 精度評價指標

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

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

2 結果與分析

2.1 研究區及數據

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

圖3 研究區內Sentinel-2影像集

2.2 降噪重建結果與分析

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

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

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

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

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

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

表1 時間序列重建效果精度統計

2.3 變化檢測結果與分析

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

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

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

圖7 不同變化檢測方法的變化檢測結果

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

表2 不同變化檢測方法精度統計

3 結束語

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