基于時空數據融合的縣域水稻種植面積提取

牛海鵬 王占奇 肖東洋

(1.河南理工大學測繪與國土信息工程學院， 焦作 454000；2.河南理工大學黃河生態(tài)文明與高質量發(fā)展研究院， 焦作 454000)

0 引言

水稻是世界第三大糧食作物[1]。中國是世界上最大的水稻生產國，水稻產量在中國農業(yè)生產中占有重要地位[2]。黃河兩岸地區(qū)是河南省水稻種植區(qū)。隨著自然、社會和經濟的發(fā)展，原陽縣水稻種植面積從約30 000 hm2下滑至約18 000 hm2。因此，適時、準確監(jiān)測水稻種植區(qū)，對糧食政策制定、糧食安全以及農業(yè)發(fā)展具有重要意義[3]。

隨著衛(wèi)星技術的發(fā)展，利用遙感技術監(jiān)測水稻種植面積已經得到廣泛應用。目前最常用的提取水稻等農作物空間分布的遙感衛(wèi)星為Landsat數據和MODIS數據[4-7]。文獻[8]使用多時段MODIS影像提取了南亞和東南亞地區(qū)的水稻分布信息，文獻[9]基于多時相MODIS數據提取了河南省水稻種植區(qū)的分布。文獻[10]利用Landsat數據并結合農戶調查，提取了翻陽湖平原水稻主產區(qū)的水稻種植區(qū)。Landsat數據具有較高的空間分辨率，但是時間分辨率較低，并且容易受云雨天氣影響，無法獲取作物關鍵生長期的時序數據[11]，因此只根據Landsat數據進行農作物區(qū)提取會受到數據缺失的限制。MODIS數據具有36個光譜波段，時間分辨率較高，但其空間分辨率低，存在大量混合像元，對地面的細節(jié)表征能力不強[12]，使基于MODIS數據的農作物提取只能在大尺度區(qū)域上取得較好的效果。

近年來，國內外學者針對遙感影像時空融合展開了研究，并取得一些成果[13-15]。文獻[16]根據像元的空間相關性和光譜相似性，提出一種自適應遙感影像融合模型(Spatial and temporal adaptive reflectance fusion model, STARFM)，該模型能較準確預測出影像的反射率。文獻[17]提出了反射率變化的時空自適應融合算法(Spatial temporal adaptive algorithm for mapping reflectance change, STAARCH)。文獻[18]在STARFM模型的基礎上提出ESTARFM (Enhanced spatial and temporal adaptive reflectance fusion model)模型，該模型能通過時間和光譜分離理論分析反射率的變化趨勢，提高融合模型預測Landsat數據的精度。文獻[19]基于類別反射率和類內像元反射率的時間變化特征，提出了STDFM (Spatial and temporal data fusion model)模型。文獻[20]對STDFM模型進行改進，提出了ESTDFM(Enhanced spatial and temporal data fusion model)模型。文獻[21]對5種主流遙感數據時空融合模型進行對比實驗，結果表明，ESTARFM模型的融合效果最好，更能反映地物的細節(jié)特征。

利用時空融合模型融合Landsat與MODIS數據對水稻種植面積進行提取的方法值得研究。本文以河南省沿黃稻區(qū)原陽縣為例，基于ESTARFM模型，對Landsat和MODIS影像進行融合，構建完整的時序Landsat數據，計算得到完整的時序NDVI曲線，采用決策樹分類的方法進行水稻種植區(qū)域的提取，以期在縣域尺度上應用ESTARFM模型提取水稻種植面積取得理想的結果，為綜合利用高空間分辨率Landsat數據和高時間分辨率MODIS數據進行農作物遙感識別和監(jiān)測提供一定的實驗基礎。

1 研究區(qū)與數據

1.1 研究區(qū)概況

原陽縣隸屬河南省新鄉(xiāng)市，地處豫北平原，位于113°36′～114°15′E，34°55′～35°11′N之間，區(qū)域總面積1 022 km2。南鄰黃河，北面是余河通道，地勢從西南至東北呈逐漸降低趨勢，地貌屬于黃河沖積平原。農業(yè)以小麥、水稻、玉米為主，其中水稻種植以一季稻為主，每年5月上旬育苗播種，10月中上旬收割。

1.2 數據與預處理

空間分辨率為30 m的Landsat數據(表1)來源于美國地質調查局(United States Geological Survey，USGS)，下載2015年品質較好、云量覆蓋小的影像。Landsat8 OLI影像已經基于地形數據進行幾何校正處理，因此本次實驗中不再對Landsat數據進行幾何校正。投影坐標系為UTM-WGS84坐標系，利用ENVI 5.3軟件對影像進行輻射定標后，使用FLAASH模型進行大氣校正。

表1 遙感數據類型及獲取日期Tab.1 Remote sensing data types and acquisition date

全年共23期的空間分辨率為500 m的MODIS13Q1數據(表1)來源于NASA網站(https://ladsweb.nascom.nasa.gov/search)。對MODIS影像利用MRT工具(MODIS reprojection tool)重投影為UTM-WGS84投影，把HDF格式轉換為GEOTIFF格式，采用雙線性內插法(Bilinear)將影像重采樣為30 m分辨率。處理后，影像與Landsat具有相同的投影方式、分辨率和像元尺寸。對Landsat與MODIS的對應波段按順序一一對應，如表2所示。

表2 Landsat8 OLI與MODIS數據波段設置及對應關系Tab.2 Band settings of MODIS and Landsat8 images

2 研究方法

本文采用ESTARFM模型將Landsat8 OLI數據與MODIS13Q1數據進行融合，從而重構出完整時序的融合Landsat數據；計算得到融合影像的NDVI數據，通過TIMESAT軟件中savitzky-Golay(S-G)濾波法進行平滑降噪；得到6種主要地物的NDVI曲線，根據6種地物的曲線特征建立專家決策樹進行分類，進而提取出水稻種植區(qū)域。具體工作流程如圖1所示。

圖1 技術流程圖Fig.1 Flow chart of technical structure

2.1 ESTARFM時空數據融合模型

本文采用的ESTARFM模型是在STARFM模型基礎上改進的。假設在理想情況下：Landsat和MODIS影像之間的反射率差異只存在系統偏差，并且2個時間的影像沒有產生較大的差異，此時ESTARFM模型預測融合影像的公式為

F(x,y,tp,B)=F(x,y,t0,B)+a(C(x,y,tp,B)-
C(x,y,t0,B))

(1)

式中F——Landsat影像的反射率

(x,y)——像元位置坐標

B——影像參與計算的波段

C——MODIS影像的反射率

tp、t0——影像獲取的時間

a——進行線性回歸時的轉換系數，由傳感器之間的系統偏差確定

由于MODIS影像中含有大量的混合像元，此時式(1)不成立。文獻[22]提出假設將混合像元的反射率建模表示為該像元中存在的不同土地覆蓋成分反射率的線性組合，其權重由其像元內混合土地覆蓋的區(qū)域覆蓋率決定。此時模型預測公式為

F(x,y,tp,B)=F(x,y,t0,B)+v(x,y)(C(x,y,tp,B)-
C(x,y,t0,B))

(2)

式中v(x，y)——對Landsat和MODIS影像同一像元進行線性回歸得到的轉換系數

考慮到相鄰的同類像元具有相似的反射率變化，文獻[16]提出一種移動窗口的方法，移動窗口用于搜索窗口內的相似像元，然后將相似像元的信息集成到所需預測的Landsat反射率計算中，此時模型預測公式為

F(xw/2,yw/2,tp,B)=F(xw/2,yw/2,t0,B)+

(3)

式中 (xw/2,yw/2)——搜索窗口中心坐標

N——包括中央預測像元的相似像元的數目

Wi——第i個相似像元的權重

Vi——第i個相似像元的轉換系數

(xi,yi)——第i個相似像元的坐標

選擇2個不同時刻m和n的Landsat和tp時刻的MODIS影像進行預測融合，用m時刻的Landsat影像和在P時刻的MODIS影像預測P時刻的Landsat影像反射率，用n時刻的觀測數據預測P時刻的Landsat影像反射率。結果分別記為Fm(xw/2,yw/2,tp,B)和Fn(xw/2,yw/2,tp,B)。通過2種預測結果的加權組合，可以獲得更精確的P時刻反射率。因此在這種情況下，為Landsat影像輸入設置更大的時間權重是合理的。因此模型公式為

(4)

最終模型預測P時刻的Landsat影像反射率計算公式為

F(xw/2,yw/2,tp,B)=TmFm(xw/2,yw/2,tp,B)+
TnFn(xw/2,yw/2,tp,B)

(5)

2.2 融合影像數據的時序NDVI曲線

歸一化植被指數(NDVI)是用于表達植被生長狀況、覆蓋情況的參數。本文利用2015年已有的質量較好的Landsat影像和ESTARFM融合預測出其他對應時間的Landsat融合影像，通過ENVI 5.3的波段運算，可以得到全年23期的NDVI數據。由于原始NDVI時序數據受噪聲和空值的影響較大[23]，造成NDVI時序曲線出現不規(guī)則的波動，所以需要對NDVI時間序列數據進行平滑降噪處理。本實驗利用TIMESAT軟件對構建的NDVI時序數據進行濾波，該軟件有非對稱高斯函數擬合、雙Logistic函數擬合和S-G濾波法，在本次研究中，通過多次實驗對比發(fā)現，S-G濾波方法基于滑動窗口的平均濾波，以局部擬合為主，有較強的細節(jié)擬合能力，能有效地去除時間序列數據中的噪聲，保留原始時序曲線的有效信息[24-25]。因此本文采用S-G濾波方法對23期的NDVI數據進行濾波得到時序NDVI曲線。

圖2 真實Landsat影像與ESTARFM模型融合影像以及STARFM模型融合影像對比Fig.2 Comparisons of real Landsat image with ESTARFM model fusion image and STARFM model fusion image

3 結果分析

3.1 ESTARFM時空融合模型結果與精度評價

根據實驗數據，對ESTARFM模型運行結果進行精度評價，并與STARFM模型運行結果進行對比。把第209天和第257天的Landsat數據和MODIS13Q1數據、第225天的MODIS13Q1數據作為ESTARFM模型的輸入數據；把第257天的Landsat數據和MODIS13Q1數據、第225天的MODIS13Q1數據作為STARFM模型的輸入數據；融合預測出第225天的ESTARFM融合數據和STARFM融合數據；以此融合數據與實際第225天的Landsat數據做對比，以顯示ESTARFM模型的預測精度。影像數據以R(Blue)G(Red)B(NIR)假彩色合成。

選取其中3塊對比明顯且具有代表性的區(qū)域(建筑物、水體、植被、農作物等)進行對比，如圖2所示。在目視解譯的基礎上，ESTARFM模型時空融合的效果相比STARFM模型效果更好，各種地物大體上紋理清晰，輪廓特征比較準確，與原始影像在空間上基本保持一致。在圖2中標記區(qū)域可以看出，相對于STARFM模型，ESTARFM模型的融合結果更加準確，與真實Landsat影像更具有一致性。只有小部分紋理比較粗糙，這是由于在構建融合數據時，使用不同時刻Landsat數據和MODIS數據存在時相的差異，這會造成融合影像與真實影像間的偏差。

為進一步驗證ESTARFM模型的融合預測效果，選取融合出的第225天的影像，分3個波段與原始影像進行對比，獲得散點圖(圖3)，從圖3看出散點基本都在1∶1趨勢線的周圍，計算真實影像與融合影像在藍波段、紅波段和近紅外波段的決定系數R2分別為0.824 6、0.920 3、0.940 2。決定系數R2都在0.82以上，且在紅波段和近紅外波段達到0.92以上，表明真實反射率與融合預測的反射率相關程度較高，因此證實了ESTARFM模型融合預測精度較高，能取得較好的融合效果，可以作為研究區(qū)水稻面積提取的數據依據。

3.2 時序NDVI曲線

結合實地調研與Google Earth影像目視解釋，將原陽縣的主要地物分為6類：建筑用地、裸地、水體、林地、小麥-水稻以及小麥-玉米。對計算得到的2015年內23個時相的NDVI數據，利用TIMESAT軟件中的S-G濾波函數對時序NDVI數據進行平滑去噪處理。利用實地調研樣點，得到每種地物的NDVI時序曲線，如圖4所示。

圖4 主要地物時序NDVI曲線Fig.4 Time series NDVI curves of main features

從圖4可以看出：主要作物為小麥-水稻輪作與小麥-玉米輪作，一年有2個波峰，在5月中下旬小麥成熟，此時收割小麥造成在儒略日第145天左右NDVI出現下降至波谷，并且由于水稻處于生長初期，稻田中含有大量的水，因此稻田的NDVI相比玉米的NDVI呈現出更低的現象；隨后水稻和玉米逐漸生長，在7月下旬至8月上旬NDVI呈現波峰。因此說明ESTARFM模型融合影像得到的時序NDVI曲線與實際作物NDVI變化情況相一致，能夠較好地反映地物的季節(jié)性變化情況，為后續(xù)提取水稻及其他地物信息提供數據基礎。

3.3 主要地物類別的NDVI閾值設置

根據各個地物NDVI曲線顯示不同的特征，本文利用設置NDVI閾值的方法，建立分層決策樹的分類方法對主要地類進行提取。具體算法如圖5所示。

圖5 分類決策樹Fig.5 Classification decision tree

(1)由圖4的NDVI時序曲線可知，植被和非植被的曲線特征差別比較明顯。非植被地物全年趨于平緩，且全年的NDVI相對于植被的NDVI明顯呈現較低的水平。在第97天時植被與非植被差異最明顯，具有較好的區(qū)分性。因此利用第97天的NDVI小于0.45作為植被與非植被用地的閾值。

(2)水體的NDVI一年中大部分都為負值，雖然有少數水體像元在第145天呈現出正值，但都小于0.10，因此用第145天的NDVI小于0.10作為水體和建筑用地、裸地的閾值。

(3)建筑用地和裸地的時序NDVI曲線大體趨勢相似，都呈現為較低值的水平。但是裸地在春夏交替時期，呈現出增長的趨勢，原因是裸地在這個時期會生長出一些野生綠色植物。因此可以用第161天的NDVI是否小于0.35作為判別建筑用地和裸地的閾值。

(4)林地與農作物的NDVI曲線相比，林地的峰值持續(xù)時間較長，呈現峰值后，NDVI沒有出現較大的波動，與水稻、玉米等農作物的時序NDVI曲線有明顯的區(qū)別，主要農作物水稻和玉米的NDVI在5月下旬至6月上旬呈現最低值，在7月下旬至8月上旬呈現峰值，因此利用波峰與波谷的差值區(qū)分開林地與農作物。分別選取100個林地與農作物的樣點，求取NDVI平均值，計算第209天與第161天的NDVI差值，經反復實驗，對比發(fā)現，利用第209天與第161天的NDVI差值小于0.22時分類效果最好。

(5)由圖4可以看出，水稻、玉米在生長初期的NDVI較低，并且水稻在生長初期田間有大量的水，NDVI比玉米更低，在7月下旬至8月上旬同時呈現峰值，且水稻的NDVI峰值比玉米略高，并且在9月下旬至10月上旬，水稻和玉米成熟，NDVI迅速下降，但下降后，水稻的NDVI比玉米的高，玉米的NDVI基本都在0.20以下，因此利用第225天與第289天NDVI的差值大于0.56作為判別玉米和水稻的閾值，第225天與第289天NDVI的差值小于0.20的判別為其他作物。

3.4 分類結果及精度驗證

根據研究區(qū)內水稻與其他主要地物的時序NDVI曲線特征，利用專家決策樹分類的方法，對研究區(qū)進行土地利用分類，提取出水稻種植區(qū)域。具體的分類結果如圖6所示。

圖6 決策樹分類結果Fig.6 Decision tree classification results

利用實地采集和在Google Earth上選取的各個地類的樣點共2 560個，通過分類后的混淆矩陣對分類結果進行定量分析(表3)，總體分類精度為92.23%，Kappa系數為0.904 3。水稻、玉米、建筑用地的制圖精度和用戶精度均在90%以上，且水稻的制圖精度達到96.73%，用戶精度也達到了93.51%。由于林地分布比較分散且與其他地類相連緊密，造成林地的混合像元問題嚴重，制圖精度和用戶精度都小于90%。

表3 分類結果精度評價Tab.3 Classification result accuracy evaluation

根據原陽縣2015年統計年鑒對分類結果進行定性分析，并且向當地農戶訪問，原陽縣大面積種植水稻的區(qū)域集中在祝樓鄉(xiāng)、太平鎮(zhèn)和葛埠口鄉(xiāng)，與分類結果一致。統計年鑒記錄原陽縣2015年的水稻種植面積為20 468 hm2，分類結果中水稻面積為18 556.38 hm2，提取水稻面積與統計數據的一致性為90.66%。

4 討論

不同地物的時序NDVI曲線具有不同的特征，因此可以利用水稻的時序NDVI將其與其他地物區(qū)分。中高空間分辨率的遙感影像在農作物識別提取中作為重要的數據源，但是由于受云雨天氣的影響，經常無法獲取連續(xù)的中高空間分辨率影像，因此利用ESTARFM模型等數據融合技術，可以解決遙感數據缺失的問題，作為一種對農作物識別提取的有效手段。但是，時空融合模型構建數據的品質受基準Landsat和MODIS影像數據的影響，獲取相近時刻并且高品質的Landsat和MODIS影像是構建高品質時序融合數據的關鍵，但是這在部分地區(qū)是比較困難的。對于零碎分散地塊，融合得到的30 m分辨率的Landsat數據依然存在明顯的混合像元問題，導致提取精度不高。另外本次研究只選取了單一的時序NDVI指數進行分類提取，可能對有些地物反映不明顯。

因此，在以后的研究中，利用ESTARFM模型融合更高空間分辨率影像和更高時間分辨率影像，進一步解決混合像元嚴重的問題，并且能監(jiān)測出短期內地物發(fā)生快速變化的情況。另外，可以開展多種指數共同對地物提取，如EVI/NDWI指數等。探討多種指數對水稻和其他農作物識別能力的適用性及準確性。

5 結論

(1)ESTARFM時空融合模型能較好地融合出缺失時相的Landsat影像，結果顯示，空間細節(jié)信息比較清晰，在Red波段和NIR波段，與真實Landsat的Red波段和NIR波段一致性程度較高，決定系數R2分別達到0.920 3和0.940 2，因此NDVI的一致性同樣較高。說明時空融合模型能夠較好地重構高空間分辨率數據，有效解決多云多雨地區(qū)受天氣影響的Landsat數據缺失問題。

(2)利用TIMESAT軟件對研究區(qū)時序影像進行S-G濾波平滑處理，得到主要地類的時序NDVI曲線，判別出水稻以及其他主要地類的不同特征，按照一定的閾值范圍設計決策樹進行分類，有較好的分類效果。總體分類精度為92.23%，Kappa系數為0.904 3；水稻的制圖精度為96.73%，用戶精度為93.51%；與統計數據的一致性為90.66%。分類精度較高，在技術上具有可行性。