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

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

0 引言

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

1.2 數(shù)據(jù)與預(yù)處理

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

表1 遙感數(shù)據(jù)源選擇Tab.1 Selection of remote sensing image

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

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

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

1.3 研究方法

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

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

1.3.1數(shù)據(jù)融合

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

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

(1)

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

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

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

m——類別數(shù)

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

(2)

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

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

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

(3)

Vti=a(ti-tj)+Vtj

(4)

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

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

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

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

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

(5)

1.3.2構(gòu)建水稻標(biāo)準(zhǔn)時序EVI曲線

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

為準(zhǔn)確地提取水稻的時序VI曲線，從灌區(qū)選取17個地面水稻樣點，保證樣點100～200 m范圍內(nèi)均為水稻，應(yīng)用融合生成的多時相VI數(shù)據(jù)提取各個樣點的EVI曲線。應(yīng)用所有樣點EVI的平均值來構(gòu)造標(biāo)準(zhǔn)時序EVI曲線，其結(jié)果如圖3所示。

圖3 水稻標(biāo)準(zhǔn)時序EVI曲線Fig.3 Standard series EVI curve of rice

1.3.3光譜耦合技術(shù)提取水稻

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

光譜耦合技術(shù)(Spectral matching technique, SMT)是指通過分析多光譜曲線與已知曲線的相似程度來對目標(biāo)對象進(jìn)行分類的技術(shù)[14, 28]?？梢詰?yīng)用該方法通過判斷目標(biāo)類別的時序EVI曲線與地面樣點構(gòu)建的標(biāo)準(zhǔn)時序EVI曲線的相似性來對目標(biāo)類別進(jìn)行分類。光譜耦合技術(shù)中應(yīng)用光譜相似度(Spectral correlation similarity, SCS)來表示多光譜曲線之間的相似程度，本文中用其來表示目標(biāo)類別與標(biāo)準(zhǔn)水稻時序EVI曲線之間的相似程度。相關(guān)計算式為

(6)

(7)

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

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

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

xi——i時刻的目標(biāo)類別EVI

yi——i時刻的標(biāo)準(zhǔn)曲線EVI

2 結(jié)果與分析

2.1 提取結(jié)果

在融合生成高時空分辨率的EVI數(shù)據(jù)后，用掩膜提取研究區(qū)域中的農(nóng)田，采用FCM將土地利用類型中的農(nóng)田分為20類，計算20個類在各個時相的平均EVI，根據(jù)EVI的平均值構(gòu)建20個類別的時序EVI變化曲線。構(gòu)建20個類的平均EVI曲線和標(biāo)準(zhǔn)時序EVI曲線的相似性矩陣。根據(jù)相似矩陣對類別進(jìn)行合并識別來提取水稻種植面積，提取結(jié)果如圖4所示。

圖4 研究區(qū)作物提取結(jié)果Fig.4 Crop extraction result of study area

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

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

2.2 提取精度

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

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

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

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

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

2.3 分辨率對提取的影響

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

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

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

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

3 結(jié)論

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

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

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