曾 群，殷守敬，田禮喬，陳曉玲

（1.華中師范大學學報編輯部，武漢430079；2.華中師范大學城市與環(huán)境科學學院，武漢430079；3.環(huán)境保護部衛(wèi)星環(huán)境應(yīng)用中心，北京100094；4.武漢大學測繪遙感信息工程國家重點實驗室，武漢430079）

通過遙感技術(shù)對地球表面的長時間、連續(xù)觀測，已經(jīng)形成了海量的時序遙感數(shù)據(jù)，為全球地表過程和環(huán)境變化研究提供了重要數(shù)據(jù)源.歸一化差分植被指數(shù)（NDVI）作為一個重要遙感參數(shù)，針對NOAA／AVHRR（GIMMS、Pathfinder）、MODIS、SPOT VEGETATION等傳感器都有專門的NDVI產(chǎn)品，形成了長達12～21年的數(shù)據(jù)積累，在全球和區(qū)域環(huán)境變化、地表覆蓋動態(tài)監(jiān)測、植被生長過程分析等領(lǐng)域得到廣泛應(yīng)用［1－2］.

大區(qū)域覆蓋、長時間跨度的NDVI數(shù)據(jù)集由于受到成像天氣條件、氣候變化、傳感器衰減等因素的影響包含了大量噪聲，盡管數(shù)據(jù)集經(jīng)過了最大值合成、云檢測等處理，但是數(shù)據(jù)產(chǎn)品仍然不同程度上受到噪聲的影響，使得對數(shù)據(jù)的分析可能會導致錯誤的結(jié)論［2］，大大限制了其在全球變化定量分析研究中的應(yīng)用.因此，對NDVI時間序列數(shù)據(jù)集進行去噪處理，剔除噪聲對影像質(zhì)量的影響，顯得尤其重要.

針對NDVI時序數(shù)據(jù)去噪問題，目前已發(fā)展出許多算法，對各類算法優(yōu)缺點也有過系統(tǒng)研究［3－5］.MVC（Maximum Value Composition）法［6］操作簡單但殘差保留過多，主要應(yīng)用于NDVI數(shù)據(jù)集初級產(chǎn)品的生產(chǎn)和對定量化分析要求不高的研究；BISE（Best Index Slope Extraction）法［7］主要用于去除NDVI突降值，但閾值設(shè)置受人為影響較大；非線性高斯函數(shù)擬合法［8］適用于周期及趨勢分析，但算法復雜不適合大數(shù)據(jù)量運算，異常值較多時會造成局部變化趨勢錯判；傅里葉變換［9］會對NDVI曲線過平滑造成變化信息的丟失，并且對NDVI數(shù)據(jù)中的偽極值比較敏感造成結(jié)果產(chǎn)生偏移；S－G（Savitzky－Golay）濾波［10］能夠清晰描述NDVI長期變化趨勢，同時保留局部的突變信息，該方法關(guān)鍵因素在于滑動窗口的大小和平滑多項式的階數(shù)，其取值需要對NDVI數(shù)據(jù)分析后確定.

NDVI作為一種典型的地理數(shù)據(jù)，具備了地物時空相關(guān)的基本特性.現(xiàn)有的NDVI時序數(shù)據(jù)重建方法，往往只注重了其時間連續(xù)性，只對時域噪聲進行了處理，而忽略了空間域噪聲的去除.本文從時空相關(guān)性出發(fā)，提出一種綜合利用S－G濾波、異常值檢測和空間域中值濾波，適用于長時間序列NDVI數(shù)據(jù)的噪聲去除方法.

1 基于時空相關(guān)的NDVI時序重建方法

1.1 基于S－G濾波的時域噪聲去除

通過對前述NDVI時域噪聲處理方法的對比，從NDVI數(shù)據(jù)的趨勢變化分析和定量化研究等需求出發(fā)，綜合考慮各種方法的噪聲去除效果、數(shù)據(jù)保真性、參數(shù)自動化程度、計算效率等因素，選取S－G濾波方法進行時域噪聲的去除.

基于S－G濾波的NDVI時序數(shù)據(jù)處理方法從提出［11］后被廣泛應(yīng)用［12－13］.該方法能夠利用NDVI長期變化趨勢（半年或一年）有效修正異常值，

隨著迭代次數(shù)的增加，F(xiàn)k呈拋物線形狀變化，先是逐漸減小達到最小值而后逐漸增大.因此，F(xiàn)k達到最小值的判斷條件可以定義如下：

1.2 年序列異常值濾波

部分區(qū)域，如南部高海拔山區(qū)，由于長期云覆蓋等條件影響在NDVI序列中連續(xù)出現(xiàn)較多的異常值，使得擬合出來的NDVI變化趨勢偏差較大，S－G濾波后也難以達到理想的結(jié)果.對此，本文提出年際異常值濾波法，利用相鄰年份同一時相的NDVI數(shù)據(jù)對S－G濾波后的值進行異常值判斷和調(diào)整.具體如下：

如果待檢測點，第t年NDVI值NDVIt如果滿足：

即，比前1年和后1年的NDVI值NDVIt－1和NDVIt＋1都低，并且超過閾值，則判定為異常值，用前后年的均值代替，T1、T2表示判定閾值，表示NDVI的年際正常波動幅度，是一個大于零的小數(shù)，一般取0.2.

1.3 基于超限中值平滑濾波的空間域噪聲去除

空間域噪聲去除方法主要是借助地物的空間相關(guān)性，根據(jù)“距離越近的地物其相似性也較高”這一原理，利用鄰域像素信息對影像進行平滑以去除噪聲點.一般的空間平滑濾波在去噪的同時，往往造成邊界模糊和信息失真.對于變化檢測等定量分析研究，類型過渡的邊界區(qū)域往往是研究的重點，在數(shù)據(jù)預處理中，要在最大程度上保留像素的原始值.基于該需求，結(jié)合現(xiàn)有的空間濾波方法，提出超限中值平滑濾波方法（Anomaly based Median Smoothing Filter，AMSF），即只是對超限像素進行中值平滑處理.

由于地物分布的空間異質(zhì)性，用固定的閾值判定像素值是否超限是不合理的.在此我們借助基于統(tǒng)計的異常點檢測方法判斷像素是否超限.方法如下：

設(shè)有數(shù)組｛x｝，首先計算｛d（x）｝＝｛｜x0－m（x）｜，…，｜x1－m（x）｜，…，｜xi－m（x）｜，…，｜xn－m（x）｜｝，m（x）為｛x｝的中值，并計算數(shù)組｛d（x）｝的中值Md和中值絕對偏差MAD：MAD＝1.4826·Md.

然后計算異常度L：

一般情況下，如果有L＞3，認為該點為異常點.

中值濾波具有較好的椒鹽噪聲去除效果，并且具有較好的保邊緣效果.在此，對異常點取窗口像素中值替代中心像素值.這樣既可以有效抑制原始影像中的噪聲，同時可以保證中心像素具備較好的同質(zhì)性和連續(xù)性.

2 數(shù)據(jù)準備與實驗

2.1 數(shù)據(jù)準備

選取江西省1998年4月1日～2009年4月1日共11年的SPOT VGT NDVI旬最大值合成產(chǎn)品S10.該數(shù)據(jù)以8位整型存儲，在產(chǎn)品中自帶了云標記產(chǎn)品，描述了成像質(zhì)量及各像素成像狀況（晴空、陰影或云等）.

對數(shù)據(jù)的預處理步驟包括：1）通過線性變換轉(zhuǎn)化為反射率值；2）利用云標記產(chǎn)品，識別出云覆蓋或陰影區(qū)域；3）利用前后時相無云數(shù)據(jù)線性插值，公式如下：同時可以體現(xiàn)NDVI的短期波動規(guī)律.為避免出現(xiàn)局部擬合不足的同時其它部分局部擬合過度的問題，對長時間序列NDVI數(shù)據(jù)，本文中以一年為一個周期進行分段擬合.具體步驟如下：

1）長期變化趨勢擬合：擬合多項式表示NDVI的長期變化趨勢，取較大的窗口寬度和較小的多項式系數(shù)可以獲取平滑的擬合曲線.對不同的NDVI數(shù)據(jù)類型，需要根據(jù)數(shù)據(jù)時間分辨率等因素對參數(shù)進行設(shè)置.

2）NDVI序列初始化：低于長期變化趨勢的數(shù)據(jù)點判定為噪音.利用長期變化趨勢線擬合出的新值，取代NDVI序列中的噪聲點，從而產(chǎn)生新的NDVI時序數(shù)列.新曲線會接近NDVI數(shù)據(jù)的上包絡(luò)線.

3）NDVI序列迭代平滑：在迭代過程中，仍然利用S－G濾波來擬合長期變化趨勢.但是重復迭代時，窗口寬度取值相對較小，多項式系數(shù)取值相對較大，以盡量體現(xiàn)出短期的變化趨勢.令第k次擬合的擬合效果系數(shù)為Fk：NDV＊t表示時間t處插值后的NDVI值，Δt1、Δt2表示該像素點上后一時相和前一時相中第一個無云日期相隔天數(shù).

2.2 時序重建

時序重建技術(shù)路線設(shè)計如圖1所示.

圖1 NDVI時序數(shù)據(jù)重建技術(shù)路線Fig.1 The technical route NDVI time series

2.2.1 S－G濾波 首先利用S－G濾波，對預處理后的NDVI數(shù)據(jù)以1年為一個時間段分段進行序列重建.S－G濾波參數(shù)設(shè)置時，滑動窗口越大，被平滑的峰谷值也就越多；平滑多項式次數(shù)越低，濾波結(jié)果會越平滑，但是也有可能會保留異常值，次數(shù)越高，就越容易去掉異常值，但是也有可能會因為擬合過度而出現(xiàn)更多的噪聲.

對長期變化趨勢模擬窗口半徑寬度取9，即取6個月的開窗大小，多項式次數(shù)取2；短期變化趨勢模擬窗口半徑寬度取5，即開窗大小為3 mon，多項式次數(shù)取4.

2.2.2 年際異常濾波 對連續(xù)長時間受云覆蓋等因素影響的區(qū)域，會產(chǎn)生連續(xù)異常值無法濾除.因此，利用歷年同一時相的NDVI數(shù)據(jù)組成年序列影像（如1998年～2008年4月1日影像組成長度為11的序列），對每個像素點利用異常值檢測，并對異常值利用公式（5）基于年序列數(shù)組進行線性插值補值.2.2.3空間平滑濾波 利用AMSF空間濾波（窗口大小取5×5）對數(shù)據(jù)逐波段平滑.

3 結(jié)果與討論

如圖2為山區(qū)林地區(qū)域樣本點預處理前后NDVI曲線.對云標記已經(jīng)標出的3個噪聲數(shù)據(jù)點A、M、N，通過線性插值進行了修復，但是對于明顯的噪聲B點，未能檢測出來.并且由于臨近噪聲點B，A點的修復效果也明顯受到影響.

圖2 樣本點NDVI時序數(shù)據(jù)云插值預處理結(jié)果Fig.2 The cloud interpolation of sample points’NDVI time series

如圖3所示，是對樣本點預處理后的NDVI序列先后利用S－G濾波和年序列異常濾波方法處理后的結(jié)果.

從圖3可以看出，通過S－G濾波，對數(shù)據(jù)進行了很好的平滑作用.不僅對噪聲點A、B，原NDVI曲線的突降值（如E點）等大多數(shù)噪聲點也通過多項式擬合檢測出來，并進行了修復，修復后的曲線更加連續(xù)，能更真實地反映植被變化情況.

利用年序列異常濾波，對C、D時間段內(nèi)的低異常值成功檢測出來并進行了修復.5～6月是林地植被茂盛季，該時間段內(nèi)出現(xiàn)了連續(xù)低值并且在之后NDVI值又恢復到峰值，很有可能是由于連續(xù)的多云天氣造成的（可以排除發(fā)生地表覆蓋類型變化、病蟲害、火災等情況），連續(xù)異常值的出現(xiàn)使得利用S－G濾波局部趨勢擬合時出現(xiàn)偏差，未能有效修復.而通過年際同時相影像的異常分析，該時間段內(nèi)數(shù)據(jù)均被檢測為異常低值，并且通過插值達到了對其進行有效修復的目的.

圖3 S－G濾波和年序列異常濾波效果圖Fig.3 The effect picture of Savitzky－Golay filter and anomaly removal algorithm

如圖4所示，經(jīng)過前述時域處理的影像，通過在空間域中的分析，可以看出數(shù)據(jù)中存在空間離散點噪聲.利用超限中值濾波，可以將空間上不連續(xù)的異常值噪聲檢測出來并修復.從圖中可以看出，濾波后的影像中大部分區(qū)域像素值沒有改變，僅對散點噪聲進行了有效濾除，同時很好地保留了邊緣部分信息.

圖4 超限像素中值平滑濾波效果圖左圖為濾波前NDVI，右圖為濾波后NDVIFig.4 The effect picture of median smoothing filter of overrun pixel

4 結(jié)語

本文從地理數(shù)據(jù)的時空相關(guān)性出發(fā)，建立了一種長時間序列NDVI數(shù)據(jù)重建方法，并利用11年時間序列的SPOT VGT NDVI旬合成數(shù)據(jù)進行了實驗.首先利用分段S－G濾波對長時間序列NDVI數(shù)據(jù)時域噪聲去除；基于年際異常信息判斷，對長時間云覆蓋等因素造成的連續(xù)異常值進行修復；然后利用超限中值濾波濾除了空間域離散噪聲.結(jié)果表明，該方法可以有效地去除時空域噪聲，為全球變化定量分析研究提供高質(zhì)量數(shù)據(jù)集.

［1］ Wang L，Chen J，Gong P，et al.Land cover change detection with a cross－correlogram spectral matching algorithm［J］.International Journal of Remote Sensing，2009，30（12）：3259－3273.

［2］ Bradley B A，Jacob R W，Hermance J F，et al.A curve fitting procedure to derive inter－annual phonologies from time series of noisy satellite NDVI data［J］.Remote Sensing of Environment，2007，106：137－145.

［3］ Jennifer N H，Gregory J M.Noise reduction of NDVI time series：An empirical comparison of selected techniques［J］.Remote Sensing of Environment，2009，113：248–258.

［4］ 李 儒，張 霞，劉 波，等.遙感時間序列數(shù)據(jù)濾波重建算法發(fā)展綜述.遙感學報，2009，13（2）：335－341.

［5］ 顧 娟，李 新，黃春林.NDV I時間序列數(shù)據(jù)集重建方法述評［J］.遙感技術(shù)與應(yīng)用，2006，21（4）：391－395.

［6］ Holben B N.Characteristic of maximum value composite images for temporal AVHRR data［J］.International Journal of Remote Sensing，1986，7（11）：1417－1434.

［7］ Viovy N，Arino O，Belward A S.The best index slope extraction（BISE）：a method for reducing noise in NDVI timeseries［J］.International Journal of Remote Sensing，1992，13（8）：1585－1590.

［8］ Jensson P，Eklundh L.Seasonality extraction by function fitting to time－Series of satellite sensor data［J］.IEEE Transaction on Geoscience and Remote Sensing，2002，40（8）：1824－1832.

［9］ Roerink G，Menenti M，Verhoef W.Reconstructing cloud free NDVI composites using Fourier analysis of time series［J］.International Journal of Remote Sensing，2000，21（9）：1911－1917.

［10］ Savitzky A，Golay M J E.Smoothing and differentiation of data by simplified least squares procedures［J］.Analytical Chemistry，1964，36：1627－1639.

［11］ Chen J，Jonsson P，Tamura M，et al.2004.A simple method for reconstructing a high－quality NDVI time－series data set based on the Savitzky－Golay filter［J］.Remote Sensing of Environment，91（34）：332－344.

［12］ 邊金虎，李愛農(nóng)，宋孟強，等.MODIS植被指數(shù)時間序列Savitzky－Golay濾波算法重構(gòu)［J］.遙感學報，2010，14（4）：725－741.

［13］ 宮 攀，陳仲新，唐華俊，等.基于MODIS溫度／植被指數(shù)的東北地區(qū)土地覆蓋分類［J］.農(nóng)業(yè)工程學報，2006，22（9）：94－99.