基于GF-1與MODIS時空融合的南方丘陵區水稻提取研究*

王建勛，華 麗※，鄧世超，王惠東， 2，陳家慧

(1.華中農業大學水土保持研究中心，湖北武漢 430070； 2.武漢大學測繪遙感信息工程國家重點實驗室，湖北武漢 430079)

0 引言

水稻作為我國重要的糧食作物，其分布十分廣泛。研究水稻空間種植分布可以為水稻精準估產與農業政策制定提供技術支持，對保障國家農業發展和糧食安全具有重要意義[1-3]。傳統方法對水稻面積的研究多基于統計數據，該方法難以獲得精確的空間分布信息。隨著遙感技術的飛速發展，基于遙感技術的水稻提取研究取得了一定成果[4-8]。但這些研究區域多集中在地勢平坦、地物類型單一的平原區，針對南方丘陵區地形復雜、地類斑塊破碎的特點，如何提高水稻分布面積提取效率和精度成為當前研究的難點和熱點。

當前水稻識別與分類主要包括兩大類方法：基于單一高分辨率影像和基于多源影像數據的水稻分類提取。高分辨率影像主要包括IKONOS和SPOT等，其優勢在于此類影像地物紋理、色調與形狀等特征較為豐富，不同地類之間差異大且容易識別。國內外一些學者利用高分辨率數據對小區域的水稻種植區進行面積提取，達到了一定的提取精度[9， 10]。但是影像的時間分辨率低，無法充分利用作物識別所必需的物候信息。而利用多源數據進行水稻的提取，主要采用中高分辨率的Landsat影像與低分辨率的MODIS影像通過數據融合構建時空序列數據集進行水稻提取[11， 12]。該類方法優勢在于影像的開源性且易于獲取，但僅適合大尺度、地類單一區域的水稻提取，難以適用于地形復雜、地類斑塊破碎的南方丘陵區。

自Gao F等[13]提出自適應遙感影像時空數據融合模型(STARFM)以來，許多學者利用該模型完成TM和MODIS數據融合，并進行地物分類與提取研究[14， 15]。但該模型存在無法解決混合像元的缺陷，Zhu等[16]針對這一缺陷提出了ESTARFM融合模型。該模型假定一段時間內地物反射率為線性變化，混合像元由不同地表覆蓋線性組合而成，相較于STARFM方法利用單時相高空間分辨率影像與高時間、中低空間分辨率影像融合，ESTARFM利用兩景高空間分辨率影像參與融合，可以減少融合過程中產生混合像元的問題，從而提高融合精度。該方法更能好地構建中高分辨率時空數據集并提供滿足用戶使用的多幅影像，從而為研究南方丘陵區水稻分布提供解決方法。

文章以湖南省衡陽縣為研究區域，采用ESTARFM方法將MODIS分別與GF-1WFV和TM數據進行融合，構建高時空分辨率的NDVI時空數據集，再利用TIMESAT軟件[17]對獲取的NDVI數據集進行濾波處理并提取地物的物候參數，最后通過決策樹分類方法提取水稻種植面積并對兩類數據融合后的提取結果進行精度評價。

1 材料與方法

1.1 研究區概況

衡陽縣地處湖南省衡陽市西北部，湘江中游，西起東經111°59′，東至東經112°46′，南起北緯26°51′，北至北緯27°26′，多為低山丘陵地貌，總面積為25.59萬hm2(圖1)。該縣屬亞熱帶季風氣候，年均氣溫17.9℃左右，溫暖濕潤的氣候成為水稻種植的天然優勢。作為全國糧食生產先進縣，衡陽縣2016年水稻播種面積9.57萬hm2，產量57萬t，是湖南省主要水稻種植區。旱地作物與水稻交錯分布現象普遍，水稻種植模式是單、雙季稻皆有，以雙季稻為主。單季稻每年5月開始播種育秧， 9月收割。雙季稻中早稻4月播種， 7月中旬收割，接著播種晚稻， 10月中旬收割。

圖1 研究區

1.2 數據及預處理

1.2.1 Landsat 8數據與GF-1WFV數據

該研究選取云量少于10%以下的Landsat 8和GF-1WFV數據，Landsat 8影像獲取時間是2016-03-01，2016-11-28，2016-12-30，該數據由美國地質勘探局提供(https：//earthexplorer.usgs.gov/)。GF-1WFV影像獲取時間分別是2016-02-29， 2016-03-28和2016-05-12，數據來源自遙感集市(http：//www.rscloudmart.com/)。Landsat 8 陸地成像儀(OLI)數據包含藍、綠、紅、近紅外和2個短波紅外波段，GF-1WFV數據包含藍、綠、紅、近紅外4個波段，波段信息如表1所示。對GF-1WFV與Landsat影像的紅波段和近紅外波段進行重投影、大氣校正、幾何校正等預處理操作。

1.2.2 MODIS數據

該研究選取8d合成的250m分辨率的地表反射率產品數據(MOD09Q1)，該數據是由NASA提供(https：//earthexplorer.usgs.gov/)的標準2級產品，包括紅波段(Red)和近紅外波段(NIR)的地表反射率，波段信息如表1所示。

數據預處理過程包括：首先使用MODIS Reprojection Tool工具將MODIS數據重投影為UTM-WGS84坐標系； 其次，將數據重采樣至30m和16m分辨率，分別與Landsat 8和GF-1WFV數據進行融合處理。用于實驗處理的數據特征如表2所示。

表1 Landsat8、GF-1WFV和MODIS影像波段信息

Landsat 8(OLI)GF1-WFVMODIS(Terra)波段號波長寬度(nm)波段號波長寬度(nm)波段號波長寬度(nm)4(Red)630～6803(Red)630～6903(Red)620～6705(NIR)845～8504(NIR)770～8904(NIR)841～876

表2 Landsat8、GF-1WFV和MODIS影像波段信息

Landsat 8(OLI)GF1-WFVMODIS(Terra)獲取時間數據用途獲取時間數據用途獲取時間數據用途2016-03-01融合驗證2016-02-29融合輸入2016-01-012016-11-28融合輸入2016-03-28融合驗證…融合輸入2016-12-30融合輸入2016-05-12融合輸入2016-12-27

1.3 研究方法

1.3.1 NDVI計算

NDVI指數是由Rouse[18]等提出，其本質是利用植物葉片在紅色與近紅外波段的光譜反差組合來設計的，是目前應用最廣泛的植被指數。該指數已被證明可以真實地反映地表植被的生長狀況。NDVI的計算公式為：

NDVI=(Nir-Red)/(Nir-Red)

(1)

式(1)中，Nir和Red各代表影像的紅色與近紅外波段。

1.3.2 ESTARFM模型

在忽略大氣校正和誤差配準影響的前提下，同一分辨率下的MODIS和Landsat數據地表反射率之間存在線性關系(式2)。

L(x，y，tk，B)=a×M(x，y，tk，B)+b

(2)

式(2)中，L、M分別代表Landsat和MODIS數據； (x，y)表示像元位置，tk為獲取影像的時間，B表示影像的波段，a、b為線性方程的系數且會隨著像元位置的不同發生改變。

當存在2個時期的2對影像，且地物類型在2個時期都沒有明顯變化，則根據式(2)可得到式(3)[19]。

L(x，y，tp，B)=L(x，y，t0，B)+a×(M(x，y，tp，B)-M(x，y，t0，B))

(3)

使用ESTARFM方法能有效獲取30m的融合結果，但該分辨率應用在南方丘陵山區的水稻分布提取顯然是不夠的，因此該研究嘗試使用GF-1數據的16m分辨率替換Landsat數據進行融合并構建時空數據集。由于GF-1數據與MODIS數據在波譜響應和Landsat數據存在差異，因此試驗中還加入Landsat融合結果作為對比分析。

ESTARFM模型是一種基于移動窗口的數據融合模型，在進行融合時，考慮從鄰近的光譜相似的像元獲取輔助信息來提高精度。在一個給定的移動的窗口內，中心像元的預測反射率可以表述為：

(4)

式(4)中，N為包括預測像元在內的相似像元的數目；ω為滑動窗口大小； (xi，yi)是第i個相似像元的位置；Wi是第i個相似像元的權重大小，可描述為對預測中心像元的貢獻程度，是依據空間距離、時間距離、光譜距離3項來確定的；Vi是轉換系數。

為提高結果的準確率，可使用兩期數據計算(xω/2，yω/2)點的反射率，公式如(5)(6)。

L(xω/2，yω/2，tp，B)=μm×L(xω/2，yω/2，tm，B)+μn×L(xω/2，yω/2，tn，B)

(5)

(6)

式(5)、(6)中，L(xω/2，yω/2，tm，B)和L(xω/2，yω/2，tn，B)分別為tm和tn時刻的觀測數據融合得到的tp時刻的反射率和二者的權重。

1.4 物候參數提取

植物物候是指植物受生物因子與非生物因子如氣候、水文、土壤等影響而出現的以年為周期的自然現象，它包括各種植物的發芽、展葉、開花、葉變色、落葉等現象[20]。利用時間序列NDVI可以提取出與植被作物長勢有關的物候參數。該研究使用TIMESAT軟件提取6個物候參數以獲取季節性植被生長規律信息：作物生長期開始(SOS)、生長期結束(EOS)、生長期長度(LOS)、生長期NDVI最大值(MON)、NDVI基準值(BON)和NDVI幅值(AON)。其中，SOS、EOS及LOS以天為單位計算，涉及NDVI的MON、BON及AON等參數指標的范圍為-1～1。

1.5 決策樹分類

決策樹分類是一種常見的分類方法，利用二叉樹的邏輯性質對研究對象進行逐步剔除。遙感應用中經常利用該方法對地物進行提取分類[21-24]。該方法最大的特點和優勢是能夠綜合運用遙感特征量和非遙感特征量等多源數據進行研究物的提取分類[25]。主要思想是：按照一定規則把遙感數據集逐層次由上向下細分從而得到具有不同屬性的各個子集。決策樹由3部分組成：根節點、內部節點(簡稱為分支)以及最終節點(簡稱為葉)，在每一個內部節點處，根據定義的規則，可以將該節點處的數據集進行劃分得到滿足條件的兩個子集，如此針對所有需要分類的對象進行相同處理，直至所有數據被分為預期設定的不同類別，即最終節點。這些規則一般基于專家知識和經驗積累進行設置，也可以定義一些算法自動提取[26]。

2 結果與分析

2.1 融合結果

該文采用ESTARFM模型分別對Landsat 8和GF-1WFV數據進行融合構建時空數據集。分別將Landsat 8的2016-11-28和2016-12-30以及GF-1WFV的2016-02-29和2016-05-12四景影像作為輸入數據，與對應分辨率的MODIS數據融合得到1年期的多景融合結果。選取Landsat 8的2016-03-01和GF-1WFV的2016-03-28兩景影像的子區域與對應時刻融合結果的子區域進行對比分析(圖2、圖3)。

圖2 融合影像與真實Landsat8影像比較及對應波段的相關性系數

圖3 融合影像與真實GF-1WFV影像比較及對應波段的相關性系數

對比兩種影像融合結果，使用ESTARFM方法融合得到的Landsat和GF-1WFV數據與真實的紅光、近紅外波段之間仍存有少許差異，且Landsat數據融合結果與真實值之間的相關性系數略高于GF-1數據結果，主要原因是MODIS傳感器與Landsat傳感器的波譜信息較為相似。但通過對兩種影像融合后的NDVI進行相關性分析，發現兩類NDVI之間的差距較小(圖4)。另外針對南方丘陵山區地形，使用16m的GF-1數據在細節表現上優于30m的Landsat數據，這對后續水稻種植區域提取有很大優勢。因此選用質量好的GF-1數據替代Landsat數據進行融合構建時序數據集是可行的。

圖4 融合影像NDVI與真實影像NDVI比較及相關性系數

圖5 不同地物的標準NDVI時序曲線圖

2.2 NDVI擬合及分類結果

分別將預測的44景GF-1和Landsat 8 NDVI數據與原始的NDVI數據構成1年期的時間序列數據導入TIMESAT軟件中，選取S-G濾波算法進行NDVI時序平滑。以GF-1WFV數據為例，該研究將研究區地類分為單季稻、雙季稻、林地、居民地和水體5類，針對不同地類各隨機選取50個樣點，以各類樣點的平均時序曲線作為各地類標準時序曲線，統計了林地、水體、居民地和單季稻的物候信息，以提取的物候信息作為分類參數，并通過決策樹分類法區分5種地類。各地物的NDVI擬合曲線如圖5所示。

根據不同地類的物候參數，利用決策樹方法建立分類規則并進行水稻識別與提取。以GF-1WFV數據為例，選取質量較好的第17景融合后的影像作為輸入影像，從NDVI統計圖可以看出，此時的水體生長期NDVI最大值(MON)處在0值附近，且大多都為負值，故用NDVI<0.15和NIR<1100建立第一個規則提取水體； 再利用統計圖及物候參數MON、AON建立第二個規則NDVI<0.25和MON<0.25和AON<0.2提取居民地； 利用統計圖以及物候參數LOS、MON建立第三個規則： 0.4

表3 基于 GF-1WFV數據的不同地物物候信息參數

SOS(d)EOS(d)LOS(d)BONMONAON林地1153162010.4 0.8 0.3 水體134290155-0.2 0.2 0.3 居民地103268 165 0.1 0.2 0.1 單季稻148 2691210.3 0.6 0.3

圖6 分類結果

表4 兩種數據分類結果面積統計

類別GF-1WFVLandsat 8象元個數面積(萬km2)象元個數面積(萬km2)雙季稻331.734 4849.2493.212 2838.91水體15.636 740.034.542 240.88居民地19.718 850.485.465 649.19單季稻38.707 099.099.615 686.54林地593.812 51 520.16171.497 81 543.48

對兩種分類結果進行像元個數統計得到各地類的面積，結果如表4所示。對比兩種數據源融合提取結果得出：居民地和水體面積統計結果比較接近，而在水稻和林地的統計結果上存在一些差異，GF-1WFV數據和MODIS融合數據提取的水稻面積大于Landsat 8融合結果，而在林地面積上，GF-1 WFV則較小。

2.3 分類精度評價

2.3.1 混淆矩陣精度驗證

結合實地采樣點與Google Earth高分辨率影像樣本點，在ENVI軟件中結合野外調查建立的解譯標志，隨機挑選出感興趣區域ROI作為標準。通過混淆矩陣對感興趣區與分類結果進行精度驗證，由于GF-1WFV數據與Landsat 8數據分辨率不同，結合兩者分辨率比例系數關系，該研究共挑選GF-1WFV驗證圖斑22.418 5萬個，Landsat 8驗證圖斑共5.303 8萬個，比例大致為4： 1關系，與兩者分辨率比例保持一致。其中GF-1WFV數據總體分類精度為86.37%，kappa系數為0.80，其中雙季稻的分類精度為86.70%，單季稻分類精度為62.26%，林地分類精度為93.29%，水體分類精度為82.98%，居民地精度為70.85%，混淆矩陣見表5。Landsat 8數據總體分類精度為80.96%，Kappa系數為0.72，各地類除林地較好外，其他分類結果與GF-1WFV分類均存在差距。根據衡陽市統計局衡陽縣2016年國民經濟和社會發展統計公報的數據[27]： 2016年衡陽縣稻谷播種面積9.566萬hm2，該結果與GF-1WFV數據統計結果十分接近，因此可以表明基于GF-1WFV數據提取的水稻分布面積為9.483萬hm2，較Landsat 8的9.254萬hm2更為準確。

表5 GF-1WFV數據混淆矩陣

類別水體居民地雙季稻單季稻林地總像元數用戶精度(%)水體22 232475987922 89397.11居民地1 03120 2881 828203623 20387.44雙季稻2 5247 90846 1122 4915 60464 63971.34單季稻2011622 8285 6381 50110 33054.58林地8052301 81890099 367103 12096.36總像元數26 79328 63553 1849 056106 517224 185制圖精度(%)82.9870.8586.7062.2693.29 總體精度: 86.37%;Kappa系數: 0.80

表6 Landsat 8數據混淆矩陣

類別水體居民地雙季稻單季稻林地總像元數用戶精度(%)水體5 219421222945 87588.83居民地2223 986616104 28593.02雙季稻7931 5949 4794512 01114 32866.16單季稻15445881118561 82161.39林地53883 04838223 05826 62986.59總像元數6 3026 13313 3981 96625 13953 038制圖精度(%)81.5264.9970.7456.8791.72 總體精度: 80.96%;Kappa系數: 0.72

衡陽縣三面環山，山地區域大面積林地可以被有效提取，研究區是雙季稻的主要種植區，耕地基本以種植雙季稻為主，因此雙季稻提取精度較高。而單季稻提取精度較低則因為單季稻種植面積較少且零星分布于雙季稻周圍，從而與雙季稻混雜在影像上構成混合像元，易被誤分為雙季稻。研究區湖泊與河流較多，水體識別比較容易，但受空間分辨率影響，少許面積較小的坑塘難以進行有效提取。居民地分為主城區和農村居民地，其中主城區能進行有效區分，而農村居民地會因為混合象元效應造成分類效果不理想。

2.3.2 地形復雜度對水稻提取影響

研究區地形復雜多樣，地塊破碎，為水稻提取增加了難度。坡度是評價地形復雜度的一個有效而且簡單的指標。該研究以5°為區分地形復雜度高低的閾值，將研究區分為兩類：一類為地形復雜度較低的區域，該區域地勢較為平坦，地形起伏小，田塊大都呈現團塊聚集且較為規整的特點； 另一類為地形復雜度較高的區域，該類區域主要位于海拔較高或地形起伏較大的低山丘陵區，水田呈現零星分布特點。針對此類特點，對兩種不同地形復雜度區域中的水稻提取精度進行了對比分析。在兩種區域內分別選取300個單季稻和600個雙季稻的驗證像元點進行了精度驗證，統計結果見表7。結果顯示：GF-1 WFV融合后的單季稻在地形復雜度低與地形復雜度高的區域提取精度分別為57.6%與56.3%，略低于Landsat 8融合后的數據的提取精度(63.3%與56.3%)，主要原因是單季稻在衡陽縣分布較少，樣本數量對精度計算有一定的影響； GF-1 WFV融合后的雙季稻在地形復雜度低與地形復雜度高的區域提取精度分別為73.7%與68.7%，明顯優于Landsat 8融合后的提取精度(68.5%與63.2%)，不同地形下單雙季稻的提取精度與混淆矩陣中的提取精度具有一致性。地形復雜度較低的區域水田主要圍繞河流與居民點進行分布，單雙季稻混種比較常見，容易造成分類混淆； 而地形復雜度較高區域的內水田主要位于河流或交通道路一側，周圍遍布灌木、果園和小片林地等其他植被，受混合像元的影響大，因此單雙季稻的提取精度低于平緩區。根據對比研究結果，說明GF-1 WFV數據應用于地形較為復雜的南方丘陵區的作物提取具有一定的優勢。

表7 平緩區與丘陵區單雙季稻提取精度對比統計

GF-1 WFVLandsat 8驗證像元正確像元錯誤像元正確率(%)驗證像元正確像元錯誤像元正確率(%)地形復雜度低單季稻30017312757.630019011063.3雙季稻60044215873.760041118968.5地形復雜度高單季稻30016913156.330016913156.3雙季稻60041218868.760037922163.2

3 結論

該文基于ESTARFM模型用GF-1WFV、Landsat 8兩種數據分別與MODIS數據進行時間序列融合重構，計算2016年衡陽縣NDVI時序數據，并獲取物候參數，通過決策樹分類方法進行水稻種植面積提取。得到以下結論。

(1)在地塊破碎、地形復雜多樣的南方丘陵區，GF-1 WFV數據的空間分辨率高于Landsat 8，故對水稻提取有一定的優勢。采用ESTARFM模型分別做了Landsat 8 與MODIS數據、GF-1WFV與MODIS數據的兩組融合處理，雖然Landsat 8數據融合結果與原始影像的紅、近紅外波段真實值的相關性優于GF-1WFV結果，但經過NDVI計算后的結果相關性差異較小。因此可以使用GF-1WFV與MODIS數據構建NDVI時序數據。

(2)利用ESTARFM模型分別構建16m和30m分辨率NDVI時空數據集，對生成的NDVI數據集進行擬合，消除噪聲干擾的同時提取各地物物候參數。利用決策樹對物候信息進行分類得到兩種融合數據的分類精度。其中16m數據集的總體分類精度為86.37%，單季稻和雙季稻的提取精度分別為62.26%和86.70%，Kappa系數為0.80。而30m數據集的總體分類精度為80.96%，單季稻和雙季稻的提取精度分別為56.87%和70.74%，Kappa系數為0.72。對比兩者分類精度，可以看出16m數據集提取的總體精度和水稻精度均優于30m數據集； 此外，對不同地形復雜度下的水稻提取精度進行對比，發現GF-1 WFV數據在地形復雜的區域提取精度優于Landsat 8數據。因此，針對南方丘陵盆地區域復雜地貌和破碎斑塊特征下的水稻提取，使用GF-1WFV數據進行融合處理提取的結果是可以信賴的，該方法的適用性可以推廣應用于同類典型區域進行作物提取。