基于多源數據融合模型的水稻面積提取

魏永霞 楊軍明 吳 昱 王 斌 SHEHAKK M 侯景翔

(1.東北農業大學水利與土木工程學院， 哈爾濱 150030； 2.農業部農業水資源高效利用重點實驗室， 哈爾濱 150030； 3.東北林業大學林學院， 哈爾濱 150040； 4.黑龍江農墾勘測設計研究院， 哈爾濱 150090)

0 引言

作物的空間分布信息是農業估產等作物研究的前提和基礎，是糧食問題和水資源分析管理等的重要依據[1]。與歐美等國家相比，我國的作物種植結構復雜，地塊小且分散[2]，這使得遙感影像的空間分辨率制約著作物空間分布信息提取結果的精度[3]，且灌區等中小尺度的農業研究需要高空間分辨率的作物空間分布信息來研究農田系統的復雜變化。

Landsat和MODIS是水稻空間分布信息提取中最常用的遙感衛星[4-5]，Landsat雖然空間分辨率較高，但是時間分辨率低，容易受云和陰雨等天氣的影響，造成作物關鍵生長發育期無衛星覆蓋，以和平灌區2015—2017年為例，作物生長期(5—9月)可獲取的Landsat數據有9～10幅，但是實際可用的只有2～3幅，且可利用影像時間分布隨機，僅僅依靠可用的Landsat數據進行作物提取將會嚴重受到數據源的制約。MODIS時間分辨率較高，但是其低空間分辨率的特征使得混合像元的數目增加，而純凈像元指數和景觀異質性是影響分類精度的主要因素[6]，這必然造成分類精度下降。IKONOS、QuickBird等商業遙感衛星雖然時間分辨率和空間分辨率均較高，但是應用成本較高，不適合長期的作物檢測。基于線性混合模型的數據融合模型[7-8]和時空自適應反射率融合模型[9]及其改進模型[10-11]等多源數據融合模型的出現為高空間分辨率的作物空間分布提取提供了研究思路和方法。為此，很多研究者[12-14]基于多源遙感數據融合模型的成果對提取作物種植面積進行了研究，以期能夠通過數據融合的手段獲得高精度的高空間分辨率作物空間分布信息。但這類模型是為波段數據融合而開發的，基于線性混合模型以初始或者始末兩期影像為基期影像來提取豐度矩陣，當中間時刻的植被指數(Vegetation index, VI)隨時間變化較劇烈時，存在融合精度較差的問題；時空自適應反射率融合模型及其改進模型以始末兩期影像在一定窗口內的高分辨率相似像元與低分辨率像元的變化為線性的假設為理論基礎，當中間時刻的像元VI變化與初始時刻的相似性較低時，也存在融合精度較低的風險。

隨著遙感分類技術的發展，作物面積提取的方法不斷進步。通過關鍵時相的VI和作物時序VI曲線提取作物面積是目前兩種主要的作物提取方法。與作物時序VI曲線提取作物種植面積不同的是，關鍵時相的VI提取作物種植面積需要充分了解研究區域的作物物候信息和種植結構等，從中找出具有明顯差別的時相來提取作物種植面積，如GUSSO等[15]在對巴西南里奧格蘭德州的多年光譜信息統計分析的基礎上，建立了該地區大豆提取模型。張榮群等[16]在研究曲周縣主要作物的生長發育期歸一化植被指數(Normalized differential vegetation index, NDVI)的基礎上，提出了該區域的作物提取模型。與其他地區相比，我國東北地區的作物種植結構簡單，主要的作物水稻、玉米和大豆均為一年一熟，且主要的生長發育期分布在5—10月。這使得作物在生長期內具有更高的光譜相似性，也一定程度增加了關鍵生育期作物提取的難度。作物VI曲線法通過分析作物時序VI曲線與目標像元的VI曲線的相似性來提取作物，不需要清楚地了解研究區域的物候信息和種植結構及能夠反映作物的動態變化信息等特點而在作物提取中具有獨特的優勢，在各種尺度的作物空間分布提取中得到了廣泛的應用[17-19]。

考慮到現存的數據融合模型在融合VI數據時存在精度可能較低的風險，為解決遙感影像數據源時空分辨率不可調和的矛盾對作物提取的限制，本文根據地物的VI變化特征提出一種多源VI數據融合模型，以常用的Landsat 和MODIS數據為數據源，提出一種適用于中高分辨率影像數據缺失區域的水稻等作物空間分布提取方法。

1 材料與方法

1.1 研究區概況

選擇北緯46°51′7.83″～47°4′8.33″、東經127°18′40.28″～127°45′12.22″的區域為研究對象，對多源遙感數據融合提取作物空間分布進行測試。研究區屬于黑龍江省綏化市，位于呼蘭河左岸的干支流河漫灘及一級階地上，地勢平坦，屬寒溫帶大陸性季風氣候，多年平均降雨量為550 mm,平均氣溫2.5℃。研究區以水稻為主要作物，除此之外，有部分的大豆和玉米等作物種植，還包含草地、林地、水體、道路和建筑物等地物。地表結構復雜，如果以MODIS數據為數據源，存在大量的混合像元。經調查，研究區單塊稻田的面積介于0.10～0.17 hm2之間， 研究區所在的慶安縣水稻種植面積大概占總農田面積的52.7%。

1.2 數據與預處理

選擇Landsat8 OLI 數據和MODIS的二級產品MOD09GA和三級VI產品MOD13Q1進行水稻面積提取，Landsat數據和MODIS數據均來自USGS官網 (https:∥earthexplorer.usgs.gov/)。如表1所示，選擇2017年作物生長期內(4—10月)可利用的Landsat遙感影像數據為數據源，其伽利略日為229、245、277，將生育期外的伽利略日為101的數據作為初始影像數據；選擇對應時期或者對應時期附近日期(VI在幾日內不會發生太大變化)可利用的MOD09GA數據為數據源，其伽利略日如表1所示；選擇作物生長期內(4—10月)MOD13Q1為數據源。

表1 遙感數據源選擇Tab.1 Selection of remote sensing image

Landsat8 OLI已經過幾何校正，應用ENVI5.3對其進行輻射定標和大氣校正。MOD09GA和MOD13Q1數據采用的是正弦坐標系統，為和Landsat8 OLI保持一致，應用MODIS批處理工具(MODIS reprojection tool, MRT)將其投影系統轉換為UTM-WGS84 52N，重采樣為240 m×240 m。根據地面樣點對Google Earth 影像進行校正，采用ENVI5.3提供的坐標轉換工具將其坐標轉換為與Landsat8 OLI一致的UTM-WGS84 52N 坐標系統。并選擇水體、道路拐角等處對Google Earth和Landsat影像位置的一致性進行檢驗，經檢驗位置一致性較好。

土地利用圖的精度對作物提取的精度和數據融合的精度具有重要影響，為提高作物提取和數據融合的精度，本文不采用下載的土地利用數據集，而是選擇2015年9月13日的Landsat8 OLI影像數據，采用支持向量機(Support vector machine, SVM)的方法將該區域的土地利用類型分成水體、草地、林地、農田和不透水層5類，其結果如圖1所示。

圖1 研究區土地利用狀況Fig.1 Land cover map of study area

1.3 研究方法

應用數據融合模型和光譜耦合技術結合的方法來提取灌區的水稻種植面積，具體流程如圖2所示。首先應用土地利用類型圖和模糊C聚類算法(Fuzzy C-means algorithm, FCM)對多期可利用Landsat數據計算得到的增強植被指數(Enhanced vegetation index, EVI)數據進行分層聚類，定義為各個土地覆蓋類型的子類。然后根據Landsat計算得到的子類平均EVI與分解MOD09GA的混合像元得到的地表覆蓋類的平均EVI的轉換系數以及轉換系數的線性變化規律。應用MOD13Q1數據的EVI產品融合生成高時空分辨率的時序EVI數據。應用FCM將融合生成的時序EVI數據分成若干類，然后根據標準水稻時序EVI曲線和各類EVI的平均值的相似性來提取水稻空間分布。其主要包括3個步驟：①數據融合。②構建水稻標準時序EVI曲線。③光譜耦合技術提取水稻種植面積。

圖2 方法流程圖Fig.2 Flow chart of algorithm

1.3.1數據融合

研究區主要作物水稻、玉米和大豆均為一年一熟。其主要生長期為5—9月，該時期內的Landsat遙感影像有9～10幅，但是受云和陰雨天氣等因素的影響，作物生長期內可用的影像資源極少，對于本文的研究區域，2015年可用的Landsat影像數據只有3幅，部分可用的有2幅(部分可用是指影像中有部分區域存在云遮蓋等現象)；2016年可用的Landsat影像數據只有2幅，部分可用的有2幅；2017年可用的Landsat影像數據只有3幅，部分可用的有2幅。且可利用影像的時間分布基本無序，因此僅僅依靠可利用的遙感影像數據進行作物的提取將會受到數據源的嚴重限制，作物提取的關鍵時相無衛星覆蓋的現象普遍存在。

按照線性混合模型的假設，低分辨率影像像元(混合像元)是高分辨率端元的線性組合[20]。基于該假設可認為低分辨影像像元的VI是其所包含的高分辨率影像類別VI的線性組合。土地覆蓋類型可用于提取豐度矩陣[21]，因此，可以將低分辨率影像與高分辨率影像類別的關系表示為

(1)

式中IC(k,ti,B)——低分辨率的混合像元k在ti時刻的平均VI

A(k,C)——像元k的豐度矩陣

ε(k,ti)——混合像元k在ti時刻計算的殘差

m——類別數

但是各種土地覆蓋類型包含不同的地物，各個地物VI各異，且隨時間變化存在不同的變化規律。為使得豐度矩陣中的同類像元具有相似的反射率和反射率變化，用可對多維數組進行聚類的FCM將各土地覆蓋類型聚成若干類，定義為各個土地覆蓋類型的子類，使得子類內部的像元具有相似的反射率和反射率變化。并假設子類內部像元的反射率變化等于子類平均反射率變化，從該假設中可以得到

(2)

式中If(k,tp)——tp時刻子類S中的像元k的VI

If(k,t0)——t0時刻子類S中的像元k的VI

鑒于直接對子類的豐度矩陣進行混合像元分解可能造成較大誤差。為此，假設土地覆蓋類型與其子類的比值Vti在一段時間內呈線性變化。從作物時序VI曲線可以看出，作物時序VI曲線在一段時間內的變化基本為線性。因此，可以假設作物時序植被指數曲線在一段時間內呈線性變化。

(3)

Vti=a(ti-tj)+Vtj

(4)

式中Vti——ti時刻子類S平均VI與其所屬的地表覆蓋類的平均VI比值

Vtj——tj時刻子類與其所屬的地表覆蓋類的平均VI比值

a——常數，是比值Vti隨時間的變化率

應用線最小二乘法可以從式(1)中計算出低分辨率影像的子類在各個時刻的平均VI(IC(k,t0,B)，IC(k,t1,B)，…，IC(k,ti,B))。應用線性回歸模型可以從式(4)中求解出參數a。

結合式(2)和式(3)可知，ti時刻像元k的VI可以表示為

(5)

1.3.2構建水稻標準時序EVI曲線

VI是對植被生長發育狀況簡單、有效的度量參數[22]，被作為特征參數應用于作物面積提取。NDVI是作物面積提取中最常用的VI[14]。考慮到大氣和土壤等對NDVI的影響，很多研究者對NDVI進行了改進，提出了新的VI。如LIU等[23]考慮到大氣和土壤的相互作用，引入了EVI，其因引入了藍光波段，能夠消除背景噪聲和大氣干擾而具有一定的優勢[22,24]，在水稻面積提取中得到了應用[25-27]。本文選擇EVI進行水稻面積提取。

為準確地提取水稻的時序VI曲線，從灌區選取17個地面水稻樣點，保證樣點100～200 m范圍內均為水稻，應用融合生成的多時相VI數據提取各個樣點的EVI曲線。應用所有樣點EVI的平均值來構造標準時序EVI曲線，其結果如圖3所示。

圖3 水稻標準時序EVI曲線Fig.3 Standard series EVI curve of rice

1.3.3光譜耦合技術提取水稻

受物力人力等因素的限制，大范圍的地面樣點調查一般很難實現。為有效識別地物，采用FCM將多時相EVI宏影像聚成若干類。該算法采用模糊數學的思想，在確定聚類中心數后，根據隸屬度來判斷一組多維數對另一組多維數的隸屬程度，使得相同類的相似性最大，不同類之間的相似度最小。其可使聚為一類像元的EVI的相似性最大，不同類之間的相似性最小。

光譜耦合技術(Spectral matching technique, SMT)是指通過分析多光譜曲線與已知曲線的相似程度來對目標對象進行分類的技術[14, 28]。可以應用該方法通過判斷目標類別的時序EVI曲線與地面樣點構建的標準時序EVI曲線的相似性來對目標類別進行分類。光譜耦合技術中應用光譜相似度(Spectral correlation similarity, SCS)來表示多光譜曲線之間的相似程度，本文中用其來表示目標類別與標準水稻時序EVI曲線之間的相似程度。相關計算式為

(6)

(7)

式中SSV——光譜相似度，用來度量兩條光譜曲線之間的相似程度

de——歐幾里德距離，值越大，表示曲線之間的相似度越小

l——時序曲線的時間序列步長

xi——i時刻的目標類別EVI

yi——i時刻的標準曲線EVI

2 結果與分析

2.1 提取結果

在融合生成高時空分辨率的EVI數據后，用掩膜提取研究區域中的農田，采用FCM將土地利用類型中的農田分為20類，計算20個類在各個時相的平均EVI，根據EVI的平均值構建20個類別的時序EVI變化曲線。構建20個類的平均EVI曲線和標準時序EVI曲線的相似性矩陣。根據相似矩陣對類別進行合并識別來提取水稻種植面積，提取結果如圖4所示。

圖4 研究區作物提取結果Fig.4 Crop extraction result of study area

研究區2017年土地利用狀況如表2所示，研究區中水稻種植面積最大，占總面積的63.24%，占農田面積的75.85%。其次為旱田，種植作物為玉米和大豆，占總面積的20.14%，占農田面積的24.15%。草地、林地、水體和不透水層的面積相對較小，占總面積的16.62%。

表2 研究區2017年各土地利用比例及面積Tab.2 Land use area and proportion of study area in 2017

2.2 提取精度

為有效地對提取結果進行驗證，采用地面樣點和Google Earth影像兩種方法對提取結果進行驗證。85個地面樣點的分類結果如表3所示，從結果可以看出，除草地的分類精度較低之外，其他地物的分類結果精度都相對較高。水稻的提取精度為0.92，宏影像的分類精度為0.91，精度相對較高。

表3 基于地面樣點的分類結果精度評估Tab.3 Accuracy assessment of classification result using ground sample

Google Earth 包含著豐富的高空間分辨率衛星影像數據，其分辨率可以達到亞米級，而且包括較新的影像數據。從Google Earth 影像中不僅可以清楚地區分出農田、城鎮、林地等，而且還可以區分出稻田和旱田，但是很難從Google Earth影像中區分出旱田作物為玉米還是大豆。由于亞米級分辨率影像的缺失使得結果的驗證較難，為此，本文隨機從Google Earth 影像中選擇150個樣點為標準點對提取的結果精度進行驗證。其結果如表4所示。

表4 基于Google Earth影像分類結果精度評估Tab.4 Accuracy assessment of classification result using Google Earth image

從表4可以看出，宏影像的分類精度為0.87，水稻面積提取精度較高，達到了0.94。草地和水體的分類精度相對較低。草地的分類精度只有0.67。從圖4可以看出，灌區的水體主要以河流的形式存在，其在Landsat影像中的寬度只有1～3個像元。這使得其分類容易受到影像幾何校正等的影響，從而將水體分到距離其較近的水田和不透水層等地物中。草地在灌區中分布分散且占的面積比較少，這使得其大多數與其他地物混合存在，從而使其容易被分為其他類。

2.3 分辨率對提取的影響

Landsat像元大小為30 m×30 m，約0.09 hm2，基本小于單塊稻田的面積，故忽略識別誤差的情況下，除邊界的稻田外，大多數稻田基本能被識別出來。MODIS像元遠大于單塊稻田的面積，混合像元中的地物大于一定比例才能夠被識別出來。假如圖5所示為1個包含7×7個高分辨率像元的混合像元，雖然其包含的高分辨像元包括水田、旱田、不透水層和水體等地物，但因水田占較大的比例，混合像元識別結果為水田，識別的水田面積會過大。

圖5 混合像元示意圖Fig.5 Mixed pixel schematic

研究區水田面積占總農田面積的63.24%，假如本文中的模型提取的結果為真值，用包含16×16個Landsat像元的網格(一個MODIS像元大概包含16×16個Landsat像元)去劃分本文中的模型提取結果(圖4)，統計每個網格中農田面積的百分比，其結果如表5所示，從結果可以看出，MODIS像元混合現象明顯。假如水田低于20%的會被識別為其他地物，則4.22%的像元被識別為其他地物，但其包含10%～20%的水田。假如水田面積大于80%的地物會被識別為水田，則14.44%的地物被識別為水田，但其包含10%～20%的其他地物。若采用低分辨率像元，這些誤差是不可避免的，從中可以看出高分辨率像元對提高提取精度具有重要意義。

表5 網格中水田面積百分比區間所占的比例Tab.5 Proportion of percentage interval of rice pixel in grid %

3 結論

(1)應用地面樣點對提取結果的評估表明，水稻提取精度為0.92，宏影像分類精度為0.91；應用Google Earth影像對提取結果的評估表明，水稻提取精度為0.94，宏影像分類精度為0.87；水體和草地等地物因分布分散而使得提取精度較差，這說明地物的分類結果受地物離散程度的影響。

(2)低分辨率像元可能包含幾種不同類型的高分辨率像元，采用低分辨率影像進行作物提取時存在被識別其他地物的像元包含一部分高分辨率水稻像元和被識別為水稻的地物包含一部分其他地物的可能性，對于地表結構復雜的地區，這種現象可能造成較大誤差，故采用高分辨率像元來提取作物對提高精度具有重要意義。

(3)通過關鍵時相的VI(如根據稻田移栽期積水的特性來提取水稻種植面積)提取作物空間分布是一種主要的作物提取方法。如果能夠通過多源數據波段融合或者VI數據融合獲得較高精度的融合數據，也可以融合生成高時空分辨率用于關鍵時相提取作物空間分布，這將會對高精度的作物提取產生重要的意義。