基于LSWI和NDVI時間序列的水田信息提取研究

劉 元 亮，李 艷,吳 劍 亮

(南京大學國際地球系統科學研究所，江蘇省地理信息技術重點實驗室，江蘇 南京 210023)

?

基于LSWI和NDVI時間序列的水田信息提取研究

劉 元 亮，李 艷*,吳 劍 亮

(南京大學國際地球系統科學研究所，江蘇省地理信息技術重點實驗室，江蘇 南京 210023)

利用多時相的Landsat TM/ETM+和DEM數據，分別采用直方圖匹配(HM)和準不變特征點(PIFs)方法對影像序列進行相對輻射校正，減少了影像的光照和大氣條件在時間上的不確定性，提高了歸一化植被指數(NDVI)和地表水分指數(LSWI)的計算精度。根據水稻在不同生育期表現出的生理特征，基于LSWI和NDVI時間序列及高程特征，采用二叉樹方法提取了浙江省金華市水田信息。經過驗證，在空間上水田信息的提取精度達到92.3%，在縣域尺度上提取面積與統計年鑒具有0.97的相關度。

水田；LSWI；NDVI時間序列；金華

0 引言

水稻是世界主要的糧食作物之一，水田面積約占世界農田總面積的11%，其灌溉用水約占全球淡水資源的70%，而且水田是溫室氣體甲烷排放變化的重要影響因素，也是碳匯的重要途徑，因此，及時、準確地了解水田的面積及空間分布狀況，對糧食安全、溫室氣體排放監測、水資源管理等方面具有重要意義[1-4]。目前，部分學者對水田信息的提取進行了相關研究：如楊艷昭等基于TM數據，采用水稻移栽期的地表水分指數(Land Surface Water Index，LSWI)等指標獲取了吉泰盆地的水田分布信息[5]；鐘仕全等基于HJ-1B衛星遙感數據，通過對HJ-1B星CCD數據水稻的光譜反射特性進行分析，建立了水稻遙感信息識別模型，應用決策樹分類方法提取了廣西賓陽縣的水稻作物信息[6]；王力凡等利用水稻成熟期的CBERS-02B衛星遙感影像對比分析了水稻與其他地物的歸一化植被指數(Normalized Difference Vegetation Index，NDVI)值和各個波段DN值的差異，并通過運算擴大數值差異來建立水稻像元提取模型，進而提取了南京市溧水縣的水稻信息[7]。但以上研究均是面向像元的信息提取方法，容易產生由“同物異譜”或“同譜異物”引起的“椒鹽效應”；另外，以上研究多是基于單一時相的遙感影像數據進行水田信息提取，由于水稻和其他植被同一季節的光譜特征類似而易造成誤分現象。

本文以浙江省金華市的Landsat TM/ETM+和數字高程模型(DEM)為基本數據，利用面向對象的分類方法，定量提取了研究區的耕地信息，根據水稻在不同生育時期表現出的光譜特征，綜合水稻的生長發育規律，基于LSWI和NDVI時間序列，在耕地信息的基礎上進一步提取水田信息，為我國作物面積的信息提取提供了新的思路。

1 研究區和數據

金華市位于浙江省中部(東經119°14′～120°46′30″，北緯28°32′～29°41′)，地處金衢盆地東段，為浙中丘陵盆地地區，地勢南北高、中間低。年平均氣溫為16.7℃(浦江)～18.2℃(永康)，年降水量為1 839.9 mm(永康)～2 052.3 mm(浦江)，四季分明，熱量雨量豐富，干濕兩季明顯，屬于亞熱帶季風氣候。金華作為重要的糧食生產基地，素有“浙中糧倉”的美譽。2011年，金華市糧食作物總產量達90.7萬t，目前，水稻仍是其主要種植作物。

本文選用的遙感影像由2008年、2009年和2010年4-11月份時相組合而成，空間分辨率為30 m，軌道號為119/40和118/40(http://www.gscloud.cn/)。由于Landsat 7 ETM+機載掃描校正器(SLC)發生了故障，造成2003年6月以后的圖像出現條帶而丟失了25%左右的數據，因此，在數據下載時利用地理空間數據云的在線條帶修復模型(多影像局部自適應回歸分析模型)進行了條帶修復，最大限度減小了由條帶缺失對信息提取結果造成的影響。

2 研究方法

本文首先基于決策樹的分類方法定量提取研究區的耕地信息(即先區分原始影像中的植被與非植被地類，再區分植被地類中的林地與耕地)，然后利用特征時間段遙感影像的LSWI以及水稻生長期的NDVI時間序列，提取研究區水田的空間分布信息。

2.1 數據預處理

首先對遙感影像進行輻射定標和大氣校正，將DN值轉換為地表反射率；其次對各期遙感影像之間進行幾何校正，使其誤差在0.5個像元以內，以減小由影像之間的位置偏差對結果造成的影響；然后對遙感影像進行拼接裁剪，得到研究區的遙感影像；最后，因為遙感影像在不同季節的成像條件存在差異，從而會使各期遙感影像的成像效果有所不同，所以需要對各期影像進行相對輻射校正，確保各期影像之間的成像效果差異最小。研究中以7月份影像為參考影像，對同季相的8月、9月份的遙感影像分別采用直方圖匹配的方法進行相對輻射校正[8]，對異季相的其他月份的遙感影像則分別采用“準”不變特征點(PIFs)法[9，10]，即選取影像之間相同地區地物類型沒有變化的像元集合，構造兩幅影像間灰度統計值的線性變換函數，進行相對輻射校正。圖1展示了從各月份影像中選取的不變地物在相對輻射校正前后LSWI值和NDVI值與參考影像值的對比狀況，數據點越靠近對角線說明與參考值越接近，即精確度越好。從圖1可以看出，經過相對輻射校正，各月份同名不變地物點LSWI值和NDVI值與參考數據的差異得到縮小。

圖1 相對輻射校正前后對比

Fig.1 Contrast before and after relative radiometric correction

2.2 研究區水稻信息物候歷

旱地作物的光譜特征和水稻類似，會對水田信息的提取產生影響。基于遙感的水田信息提取主要依據水稻與其他旱地作物的物候期差異進行。本文利用水稻生長期影像的植被指數，通過其在不同生長期表現的物候特征對水田進行提取。通常，水田的動態變化一般包括灌水與秧苗移栽期、秧苗生長期、成熟收割期、休耕期。根據研究區農業氣象觀測數據可知水稻生長發育的時間分布特征，單季稻一般6月中下旬進行灌水移栽，10月上旬進入成熟收割期，其中8月為水稻的生長旺季；雙季早稻則在4月下旬至5月上旬進行移栽，7月下旬至8月上旬成熟，其中6月份進入生長旺季；雙季晚稻則在8月上旬灌水移栽，10月下旬收割，其中9月份生長最旺盛。

在水稻移栽期間，需要大量灌溉，其水深通常在2～15 cm，影像中水田表現出水體和秧苗的混合光譜特征，圖2展示了典型水體、秧苗及移栽期含水體秧苗的光譜特征。數據分別采樣自7月、8月、6月的樣本。

旱地農作物(如棉花)不具備在第5波段反射率下降這一特征[11]。因此本研究首先采用對土壤濕度和植被含水量較為敏感的地表水分指數(LSWI)提取大部分水田信息，計算公式為：

LSWI=(ρnir-ρswir)/(ρnir+ρswir)

(1)

圖2 典型地物光譜特征

Fig.2Thespectralcharacteristicsoftypicalobjects

式中:ρnir是近紅外波段(TM或ETM+的第4波段)的反射率，ρswir則是短波紅外波段(TM或ETM+的第5波段)的反射率。由于6月和8月分別是單季稻和雙季晚稻的灌水移栽期，與生育期相同的旱地作物相比，地表濕度較大，此時相水田的LSWI值相對于這些旱地的LSWI值高，而在其他時相會由于植被含水量的原因使得旱地和水田不易區分，因此，可以通過設置6月份和8月份恰當的LSWI閾值將大部分的水田信息提取出來[5]。圖3(見封2)是研究區耕地樣本6月份和8月份LSWI變換圖像，圖3a中藍色地物為單季稻移栽期水田，綠色地物為旺盛期雙季早稻田，粉紅色地物為旱地；圖3b中藍色地物為雙季晚稻移栽期水田，綠色地物為旺盛期單季稻田，粉紅色地物為旱地；圖3c和圖3d分別為圖3a和圖3b的LSWI變換圖像，從圖中可以看出，旱地的LSWI值較低，因此可以通過設置6月份和8月份的LSWI閾值剔除旱地信息，保留水田信息。但是，某些旱地作物在6月和8月生長較為茂盛，受植被含水量的影響其LSWI值較大；加之6月和8月恰逢我國南方多雨時節，降水會導致某些旱地表面被雨水覆蓋，表現出和水田類似的光譜特征，LSWI值相對于普通旱地作物會偏高。因此，基于LSWI提取的水田可能含有部分旱地信息。

NDVI時間序列能夠反映植被隨季節變化的規律，由于不同植被在各季節表現出的生理特征不同，其植被指數時間序列差異較大，因此，利用植被指數時間序列進行植被信息提取，能夠區分同一時段內無法區分的植被，減少因為季節變化而造成的誤分現象[12]。因此,本文在利用LSWI提取水田的基礎上，又根據水稻在不同生育期表現出的生理特征，利用NDVI在整個水稻生育期的變化狀況，進一步提取水稻信息。NDVI的計算公式為：

NDVI=(ρnir-ρred)/(ρnir+ρred)

(2)

式中:ρnir是近紅外波段(TM或ETM+的第4波段)的反射率值，ρred是紅光波段(TM或ETM+的第3波段)的反射率值。

根據水稻在各個生育期的物候特征可知，單季稻在6月份進行移栽，由于地表水體的影響，NDVI值較低，隨后稻苗開始返青，NDVI值逐漸升高，到8月份進入生長旺季，NDVI值達最大，隨后進入成熟收割期，NDVI值開始下降；雙季早稻則在4月、5月份進行移栽，NDVI值較低，隨后開始升高，到6月達最大；雙季晚稻則在7月下旬進行移栽，9月進入生長旺季，NDVI值達最大。根據實測數據和谷歌地球高清數據，各選取6個雙季稻和單季稻樣本點，分析相對輻射校正前后水稻在其生育期內NDVI值的動態變化，如圖4和圖5所示。

由圖4可以看出，相對輻射校正后雙季稻在生育期內(5-10月份)具有兩個NDVI峰值，分別出現在6月份和9月份，單季稻在生育期內(6-10月份)具有單個峰值，出現在8月份，校正后水稻的NDVI值的動態變化狀況可以基本反映水稻在生育期的生長狀況，從而可以作為區分水田和旱地的依據。

圖4 單季稻(左)和雙季稻(右)NDVI曲線(校正后)

Fig.4 The NDVI value of single cropping rice(left) and double cropping rice(right) in different months (after correction)

圖5 單季稻(左)和雙季稻(右)NDVI曲線(校正前)

Fig.5 The NDVI value of single cropping rice(left) and double cropping rice(right) in different months (before correction)

由圖5可以看出，相對輻射校正前雙季稻在生育期內雖然也具有兩個NDVI峰值，但在其他月份的NDVI值變化不盡一致，單季稻NDVI峰值出現的時間具有不確定性，因此，相對輻射校正前水稻的NDVI值隨時間的變化與實際情況不相符，由此提取的水田信息必然存在較大偏差。

2.3 分類方法

利用eCognition軟件，將分割尺度設為30，形狀參數設為0.15，緊湊度設為0.5，對相對輻射校正后的遙感影像進行分割，得到含有語義信息的影像對象。采用二叉樹分類體系:二叉樹第一層兩個葉片為植被和非植被，根據7月份NDVI值，通過設定恰當的閾值T1將植被提取出來。二叉樹第二層葉片為林地和耕地，其光譜特征類似，但是植被指數隨時間變化的規律差異較大，林地的NDVI值長期較高且保持穩定，而耕地的NDVI值變化相對較大，因此，利用5-10月份的NDVI值大于某一數值T2可剔除林地；此外，耕地一般分布在海拔較低的平原或山谷地帶，所以通過設定恰當的高程閾值T3可進一步提取耕地信息。經過試驗，用于分類提取的閾值最終設置如下：T1為0.25，T2為0.65，T3為330。經過影像分割、分類規則集構建、多閾值分類及后處理等過程，提取耕地信息(圖6a)。二叉樹第三層葉片為水田和旱地，以提取得到的耕地信息為基礎，將LSWI的閾值設置為T4=0.06，T5=0.08，提取規則如下:

其中：g(x,y)=1代表此對象為水田，g(x,y)=0則代表此對象為旱地，LSWI06和LSWI08分別為6月份和8月份的LSWI。

根據以上規則提取了混有旱地的水田，然后計算提取的影像對象在各個時相的NDVI值。通過對單、雙季稻樣點各時相NDVI統計分析，采用以下規則提取水田信息(圖6b)：1)在水稻生長期內NDVI有兩個峰值，分別出現在6月份和9月份；2)在水稻生長期內NDVI有一個峰值，出現在8月份。本文主要采用基于對象水平的二次差分方法對單峰值和雙峰值進行識別，具體過程如下，

D1i=NDVIi+1-NDVIi

(3)

D2i=1ifD1i>0; D2i=-1ifD1i≤0

(4)

D3i=D2i+1-D2i

(5)

根據式(3)計算每個對象前后兩個時相的NDVI差值，得到D1序列；根據式(4)判斷規則得到D2序列；最后由式(5)計算D2前后兩元素之差，得到D3序列。在D3序列中，若有兩個元素的值為-2，且位置分別在2和5，則表示在水稻生育期有兩個峰值且分別在6月份和9月份；若只有一個元素的值為2，且位置為4，則表示在水稻生育期內有一個峰值且出現在8月份。

圖6 耕地、水田和旱地信息提取

Fig.6 The spatial distribution of cultivated land,paddy field and dry land

3 結果與分析

基于以上水稻識別規則提取了水稻信息，經過適當的修正處理得到了研究區水稻信息的空間分布狀況(圖6b)。結合實驗區野外實測數據和谷歌地球高清數據目視判讀相結合的方法，對研究區選取的樣點進行相關指標計算，對耕地和水田提取結果進行了精度分析和評價。從表1和表2可見，耕地信息提取的總體精度達93.7%，Kappa系數為0.874；水田信息提取的總體精度達92.3%，Kappa系數為0.84。此外，還用相同方法對相對輻射校正前影像信息提取的結果進行了精度驗證，發現校正前耕地信息的總體精度為86.9%，Kappa系數為0.737；水田信息提取的總體精度為85.6%，Kappa系數為0.712。

表1 耕地精度評價

Table 1 Accuracy evaluation of cultivated land

地物類別驗證樣本耕地其他總計用戶精度(%)耕地其他總計生產精度(%)總精度(%)102510795．393．789210090Kappa系數110972070．87492．794．8

表2 水田精度評價

Table 2 Accuracy evaluation of paddy land

地物類別驗證樣本水田其他總計用戶精度(%)水田其他總計生產精度(%)總精度(%)97710493．392．399610591．4Kappa系數1061032090．8491．593．2

分別對研究區各個縣區進行水田面積驗證，由于研究區缺少2010年土地利用面積統計年鑒，但土地利用的年際變化較小，因此本文利用2007年統計年鑒的土地利用面積數據與本研究結果進行對比。由統計年鑒可知，研究區水田總面積為1 420.41 km2,本文提取的水田總面積為1 404.91 km2，面積一致性達R=0.989。在縣區尺度上，由統計年鑒得到的水田面積和本實驗提取的水田面積具有較好的相關性(圖7)。由圖6b可知，金華地區水田呈現西多東少的態勢，主要集中在河流兩側的河谷地帶及水庫周邊等灌溉條件較好的地區，在灌溉條件較差的山地丘陵地區，水田分布較少。

圖7 面積驗證

Fig.7 The validation of area in county scale

4 結論

本文基于LSWI和NDVI時間序列，利用2010年前后水稻生長期的多時相遙感影像數據，依據水稻在各個生長期與其他作物表現出的生理差異特征，通過一定的識別規則提取了浙江省金華市水田的空間分布信息，得出以下結論：1)6月和8月LSWI值較高；雙季稻NDVI峰值出現在6月和9月，單季稻NDVI峰值出現在8月。2)相對輻射校正使得耕地及水田的提取總體精度分別比校正前提高了6.8%和6.7%。3)二叉樹水稻提取方法綜合利用了遙感及空間數據在時間、空間兩方面的信息優勢，獲得了92.3%的精度，在縣域尺度上與統計年鑒具有0.97的相關度。

但是，本研究中未對未耕種水田信息進行處理,數據成像時間受TM數據的時間分辨率和天氣狀況影響，不容易獲取滿足作物關鍵生長期的影像，今后可以考慮利用更高時間分辨率數據進一步提高作物信息提取的精度。

本文工作受軟件新技術與產業化協同創新中心部分資助，此致謝忱！

[1] XIAO X M,STEPHEN B,STEVE F,et al.Mapping paddy rice agriculture in south and southeast Asia using multi-temporal MODIS images[J].Remote Sensing of Environment,2006,100:95-113.

[2] 張友水,原立峰,姚永慧.多時相水田信息提取研究[J],遙感學報,2007,11(2):282-288.

[3] DENIER VAN DER GON H.Changes in CH4 emission from rice fields from 1960 to 1990s:Impacts of modern rice technology[J].Global Biogeochemical Cycles,2000,14(1):61-72.

[4] 鄭長春,王秀珍,黃敬峰.多時相MODIS影像的浙江省水稻種植面積信息提取方法研究[J].浙江大學學報(農業與生命科學版),2009,35(1):98-104.

[5] 楊艷昭,鄭亞南,姜魯光,等.吉泰盆地水稻熟制遙感識別研究[J].地理與地理信息科學,2014,30(3):16-20.

[6] 鐘仕全,莫建飛,陳燕麗,等.基于HJ-1B衛星遙感數據的水稻識別技術研究[J].遙感技術與應用,2010,25(4):464-468.

[7] 王力凡,潘劍君.基于CBERS-02B衛星影像光譜系想你的水稻種植面積提取方法——以南京市溧水縣為例[J].南京農業大學學報,2013,36(1):87-91.

[8] 潘志強,顧行發,劉國棟,等.基于探元直方圖匹配的CBERS-01星CCD 數據相對輻射校正方法[J].武漢大學學報(信息科學版),2005,30(10):925-927.

[9] 張友水,馮學智,周成虎.多時相TM影像相對輻射校正研究[J].測繪學報,2006,35(2):122-127.

[10] 張鵬強,余旭初,劉 智,等.多時相遙感圖像相對輻射校正[J].遙感學報,2006,10(3):339-344.

[11] 魏新彩.基于HJ衛星影像的水稻種植面積遙感信息提取方法研究[D].武漢:湖北大學,2013.

[12] 朱滿,胡光宇,于之峰.基于融合NDVI和EVI時間序列的遙感影像分類研究[J].遙感信息,2009(5):44-46.

Study on Extraction of Paddy Fields Based on LSWI and Time-Series NDVI

LIU Yuan-liang,LI Yan,WU Jian-liang

(InternationalInstituteforEarthSystemScience,NanjingUniversity,JiangsuProvincialKeyLaboratoryofGeographicInformationScienceandTechnology,Nanjing210023,China)

Rice is one of the most important food crops in China.Timely and accurately acquiring the area and spatial distribution of paddy fields is significant to the crop yield estimates,food security and water resources management and so on.In this paper,the multi-temporal Landsat TM/ETM+ images and DEM are used and histogram matching and pseudo invariant feature spot methods are employed respectively to make relative radiometric correction for the TM/ETM+ image sequence.It reduces the uncertainty of the lighting and atmospheric conditions in the temporal phase for the images,and improves the accuracy of the NDVI and LSWI.Firstly,the cultivated land of study area is extracted quantificationally according to multi-temporal NDVI,then based on the physiological characteristics of rice in different growth stages,the paddy field containing partial dry land is extracted by defining the threshold of LSWI in June and August.On the basis,according to the feature that the NDVI of double cropping rice has two peaks which appear in June and September respectively and the NDVI of single cropping rice has one peak in August,the paddy field is extracted more accurately.The results show that the precisions of paddy field extraction are as high as 92.3%.Besides,the correlation coefficient between extraction area and statistical yearbook is 0.97 in county scale.

paddy field;LSWI;time-series NDVI;Jinhua

2014-11-09；

2015-01-16

國家自然科學基金項目“遙感影像小支撐規范化去卷積快速算法研究”(41371331)

劉元亮(1990-)，男，碩士研究生，主要從事遙感與GIS應用方面的研究。*通訊作者E-mai:liyan@nju.edu.cn

10.3969/j.issn.1672-0504.2015.03.007

TP751

A

1672-0504(2015)03-0032-06