基于時(shí)空相關(guān)性的NDVI時(shí)序重建方法研究

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

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

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

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

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

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

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

1.1 基于S－G濾波的時(shí)域噪聲去除

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

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

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

1.2 年序列異常值濾波

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

如果待檢測(cè)點(diǎn)，第t年NDVI值NDVIt如果滿足：

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

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

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

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

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

然后計(jì)算異常度L：

一般情況下，如果有L＞3，認(rèn)為該點(diǎn)為異常點(diǎn).

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

2 數(shù)據(jù)準(zhǔn)備與實(shí)驗(yàn)

2.1 數(shù)據(jù)準(zhǔn)備

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

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

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

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

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

2.2 時(shí)序重建

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

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

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

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

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

3 結(jié)果與討論

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

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

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

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

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

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

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

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

4 結(jié)語(yǔ)

本文從地理數(shù)據(jù)的時(shí)空相關(guān)性出發(fā)，建立了一種長(zhǎng)時(shí)間序列NDVI數(shù)據(jù)重建方法，并利用11年時(shí)間序列的SPOT VGT NDVI旬合成數(shù)據(jù)進(jìn)行了實(shí)驗(yàn).首先利用分段S－G濾波對(duì)長(zhǎng)時(shí)間序列NDVI數(shù)據(jù)時(shí)域噪聲去除；基于年際異常信息判斷，對(duì)長(zhǎng)時(shí)間云覆蓋等因素造成的連續(xù)異常值進(jìn)行修復(fù)；然后利用超限中值濾波濾除了空間域離散噪聲.結(jié)果表明，該方法可以有效地去除時(shí)空域噪聲，為全球變化定量分析研究提供高質(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í)間序列數(shù)據(jù)濾波重建算法發(fā)展綜述.遙感學(xué)報(bào)，2009，13（2）：335－341.

［5］ 顧 娟，李 新，黃春林.NDV I時(shí)間序列數(shù)據(jù)集重建方法述評(píng)［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)，宋孟強(qiáng)，等.MODIS植被指數(shù)時(shí)間序列Savitzky－Golay濾波算法重構(gòu)［J］.遙感學(xué)報(bào)，2010，14（4）：725－741.

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