陳 濤 周世健 陶 歡 侯藝璇
(1. 東華理工大學 測繪工程學院, 江西 南昌 330013; 2. 南昌航空大學, 南昌 330063;3. 中國科學院地理科學與資源研究所, 北京 100101)
油菜是我國最主要的油料作物之一,其種植面積占總油料作物種植面積的50%,產量大約占總油料產量的40%以上[1-3]。目前我國關于油菜種植分布與面積的數據是基于各地區縣域統計上報的結果,更新周期較長。遙感技術的特點市快速、無破壞和大面積監測優勢,對油菜種植面積監測提供了新的方法。油菜不同生長階段對應的光譜反射率不同,表現在遙感影像上呈現出光譜信息的差異,為基于多時相的油菜遙感監測提供理論基礎。目前已有一些研究結合農作物物候創建時間序列數據,以及植被指數,并將兩者應用到農作物提取。在大尺度制圖研究選取的影像大部分為中分辨率成像光譜儀(Moderate-resolution Imaging Spectroradiometer,MODIS)數據,王凱等[4]利用2008—2013年75個時相的MODIS-NDVI(歸一化植被指數,Normalized Difference Vegetation Index)時序數據,將實地樣本數據與油菜物候信息相結合,建立油菜種植面積提取模型,通過多次閾值設定提取了2009—2013年湖北省油菜種植面積。尤慧等[5]基于250 m空間分辨率的MODIS-EVI遙感數據,TM數據作為野外采樣數據與MODIS EVI數據之間的過渡數據,最終,經過閾值設定提取2014—2015年間江漢平原油菜種植面積分布。在小尺度制圖方面,基于高空間分辨率影像數據,從單一影像向時間序列影像數據發展。柴振剛等[6]為提高Sentinel-1 SAR(合成孔徑雷達,Synthetic Aperture Radar)數據作物種植分布提取精度,以湖北省江陵縣為研究區域,基于資源三號衛星CCD(電荷耦合元件,Charge-Coupled Device)融合數據(空間分辨率為2 m),對田間邊界進行提取,采用SAR數據進行閾值分類、通過SAR不同時期的影像數據進行田間信息提取的比較,最終提取油菜種植分布信息。肖善才等[7]基于Landsat-8和資源3號影像,利用決策樹方法,提取南京市2016年油菜種植分布信息,并進行驗證。韓濤等[8]利用Sentinel-2A、Landsat-8影像以及野外調查數據,基于影像的光譜特征與植被指數信息利用不同分類方法提取南京市高淳區的油菜種植面積。
由上述研究可知,在遙感數據的選擇上,從單一影像到時間序列影像的發展。基于單一影像農作物提取的方法,該方法在確定關鍵物候期時選取監督分類的方法對農作物進行提取;基于時間序列影像數據農作物提取方法包括:有基于單一特征參量的方法(閾值法-NDVI或EVI)、基于多特征參量的方法,也常用閾值法,除對農作物NDVI或EVI時間序列曲線構建之外,還應該考慮各個波段的影響[9-10]。因油菜種植結構較為復雜,必須對生長期的狀況與NDVI進行實時監測,所以本文選取時間分辨較高,覆蓋范圍較大的MODIS影像。
2015年,湖南省油菜種植統計面積為1.31×106hm2、總產量210.8×104t,分別占全國比重為17.4%、14.12%,居全國前列[11]。本研究以湖南省為研究區,研究的內容包括:(1)利用MODIS-NDVI時序數據和油菜實地采樣調查數據,獲取油菜像元MODIS-NDVI的時序圖,構建湖南省油菜生長期的NDVI時序標準曲線;(2)采用最小二乘法計算湖南省時序影像的每個像元與標準曲線的相似度,相似性越高則認為該像元是油菜的概率越高;(3)基于湖南省各區縣統計上報的油菜種植面積結果,對油菜種植面積分布的相似度閾值進行優化,獲得湖南省油菜種植的空間分布,并估算在該最優閾值下油菜種植的空間分布誤差。
湖南省(24°38′~30°08′N, 108°47′~114°15′E)地處云貴高原向江南丘陵和南嶺山脈向江漢平原過渡的地帶,地勢呈現三面環山、包括平原、河湖等地貌,主要以丘陵為主,是我國油菜種植的優勢產區。全省年平均氣溫20℃左右,冬季最冷月(1月)平均溫度在4℃以上。春、秋兩季平均氣溫大多在16~19℃之間。夏季平均氣溫比春、秋氣溫高10℃左右,年降水量大約在1 500 mm,屬亞熱帶季風氣候。油菜播種氣溫與之相差不大,因此研究區非常適合油菜的種植。
本文采用的數據來自美國NASA官網提供的2017年MOD13A1 16天合成的空間分辨率為500 m的產品。將所獲得的30景影像利用MRT軟件進行拼接、選取,把所有影像的數學基礎都統一為GCS_WGS_1984的地理坐標系,投影到WGS_1984_UTM_Zone_50N投影坐標系上,最終進行重采樣成500 m。土地利用類型數據使用的清華大學宮鵬團隊發布的2017年全球10 m的土地類型數據(http:∥data.ess.tsinghua.edu.cn/)[12]。在ArcGIS中,將13景影像拼接、裁剪,并重采樣至500 m分辨率,最終提取出耕地類型圖,而湖南省縣域油菜種植面積通過統計年鑒獲得。
油菜NDVI標準物候曲線的獲取通過湖南省長沙市、常德市、衡陽市綜合實驗站獲取GPS經緯度坐標,調查地點的主要分布在長沙市、衡陽市、常德市、懷化市。樣點選取原則是盡量保證油菜端元是純凈的,對油菜品種具有代表性的樣點進行選取[13]。在經過Google Earth確定農田信息侯,通過ArcGIS得到矢量數據并投影成與MODIS數據坐標相一致,最終從遙感數據中提取出對應像元的NDVI作為油菜的端元,確定湖南省油菜NDVI標準曲線。
2.2.1植被指數NDVI
構建NDVI[14],它在遙感影像進行植被覆蓋與植被物候研究得廣泛應用,對植被生長狀態、植被覆蓋度都有著較高的監測效果,并在消除大氣、地形以及其他輻射干擾具有很好的作用。
2.2.2 Savitzky-Golay濾波
在光譜分析中平滑濾波是常見的預處理方法之一。利用Savitzky Golay方法的優勢在于可以降低噪音的干擾、提高光譜的平滑性[15-16]。而且S-G濾波最大的特點在去除噪聲之后還能保證信號的形狀和寬度不變。
2.2.3與標準曲線判定相似性
通過綜合實驗站的實地數據提取相應遙感像元,形成標準的油菜端元NDVI時序曲線圖,隨后利用ENVI(完整的遙感圖像處理平臺,The Environment for Visualizing Images)中的bandmath進行相似性計算,一次來判斷每個像元的NDVI時序曲線與標準曲線的相似性,確定與標準像元相似的程度。如計算公式(1):
(1)
式中:xi為10個時相對應的NDVI值;yi為標準曲線中NDVI的值;min的值表明像元NDVI曲線與標準NDVI時序曲線在一定值域范圍內的擬合程度。
2.2.4相關系數及指標
為評估最終提取結果的精度與效果,本文利用皮爾遜相關系數(r)、決定系數(R2)和均方根誤差(RMSE)進行評估,具體公式參考Carlos Antonio da Silva Junior[17]。
本文的研究思路:利用清華大學宮鵬團隊發布的2017年全球10 m的土地類型數據[12]從中提取出湖南省耕地類型數據與湖南省的時序,MODIS-NDVI數據進行掩膜提取,得到耕地上的NDVI;再與實地的采樣點數據(濾波之后)進行匹配;最終將每個像元的曲線與標準曲線之間進行最小二乘計算,尋找最優相似度閾值的范圍,從而得到相應的油菜面積。具體的技術路線如圖1所示。
通過實地采樣數據的油菜端元進行取平均值,得到油菜主要生長期內的NDVI時間序列曲線(圖2)。NDVI等于0.4一般作為植被生長期的開始[18]。圖2中可以看出油菜的NDVI的開始值在0.43,冬油菜NDVI時序曲線整體趨勢,1~2月處于油菜的蕾臺期,油菜將會停止生長,NDVI變化很小;3~4月正值油菜開花期葉綠素增加,NDVI曲線快速上升;5月處于成熟收獲期,NDVI曲線開始下降,而在儒略日081(4月7日)時NDVI最高,是油菜提取最佳時相。根據Savitzky-Golay濾波后的圖像,剔除異常值使曲線變得更加的平滑。

圖2 油菜NDVI時間序列曲線
經過與標準曲線進行最小二乘計算,圖3是湖南省耕地中與標準油菜曲線的相似性分布圖,每取一次概率閾值就會得到相應的油菜分布,將其與統計年鑒數據進行比較,確定最優的取值范圍為0.2~0.7對應的油菜分布如圖4所示,此時的相關性最高。研究區采用的是空間分辨率為500 m的MODIS數據,因此一個油菜像元的面積對應的是地面上的500 m×500 m區域,只需要統計湖南省油菜像元數量即可計算油菜的種植面積,計算公式如下:
(2)
式中:S為提取面積,單位hm2;N為研究區作物提取像元數量;R為MODIS數據空間分辨率(500 m)。
利用ENVI中的分區統計工具對提取結果進行計算,得出共有137 732個油菜像元,根據公式(2)計算,遙感提取湖南省油菜種植面積為3.87×106hm2。
為了確定各個市級縣區油菜提取面積對應相關性,分別將每個市級區縣矢量圖與湖南省提取結果做掩膜處理并進行分區統計。然后將分區統計的各市級縣區的結果與統計年鑒數據做相關性分析。經Pearson相關系數計算,14個市中3個市的縣級相關系數在0.95以上。3個市的縣級相關系數在0.90~0.95之間,3個市的縣級相關系數在0.75~0.85之間,5個市的縣級相關系數在0.35~0.7之間。可得出,大部分市級油菜提取面積與統計年鑒數據的相關性較高,表明各市能達到一定精度要求。

圖4 最優參數下的油菜種植分布注:審圖號GS(2019)3266號。
本文利用遙感方法提取的湖南省2017年油菜面積與統計面積作驗證(數據來源于湖南統計年鑒),如圖4所示,根據相似性閾值范圍進行取值,然后利用湖南省統計年鑒數據(縣級)與遙感提取面積數據的相關性分析結果可知,當閾值取0.7時相關性最高,以此進一步確定閾值范圍在0.2~0.7之間。遙感提取面積的相關系數(r)達到0.81,其值域等級為強相關[19],表明提取結果有較高的可靠性。從湖南省縣級油菜遙感提取結果與統計年鑒數據的二維散點圖所示(圖5),R2為0.66,RMSE為31.4。精度分析表明本文方法在單一農作物遙感提取中有良好的精度。

圖5 MODIS-NDVI提取結果與統計年鑒數據二維散點圖

表1 湖南省縣級統計年鑒數據與MODIS遙感提取面積的差值 單位:hm2
從單一影像發展到基于時間序列影像,利用油菜主要的物候期與MODIS時間序列數據相結合[9-10,20-21],通過遙感影像分析獲得的主觀數據和客觀數據之間的比較可能會導致兩種方法之間的差異,造成這些差異的可能是南省油菜種植期多云多雨,并插花種植嚴重,加之地塊破碎,極易產生混合像元的現象[17,22-25]。為此對遙感結果與統計結果差異較大區域進行分析,如表1所示,其中寧鄉市的油菜面積被高估了1.04×105hm2;而在懷化市、長沙市、張家界市、郴州市四個地區中有十個縣區出現低估的現象,特別是在安化縣表現出低估值最大最明顯。
本文利用MODIS-NDVI數據描述了我國冬油菜在湖南省的空間格局分布,結合研究區的物候歷,構建最小二乘的方法,分析了冬油菜在湖南省的空間格局,并提取冬油菜面積,同時利用統計年鑒數據對MODIS遙感提取結果進行精度驗證與相關系數檢驗。得出以下結論:MODIS數據具有高時間分辨率、光譜范圍廣、更新頻率高,覆蓋了農作物的生長期,為大尺度制圖提供了新的數據源。研究結果顯示,與統計年鑒數據相比,縣級相關性系數(r)達到0.81,R2達到0.66。此方法操作流程簡單,因此選擇合適的MODIS時間序列影像數據,為油菜面積遙感提取提供了可行性。