基于MODIS NDVI數(shù)據(jù)的遼寧省玉米種植面積提取研究

宮詔健,田景仁,陳 杰,朱明宇,宋文錦,朱 紅,黃薈羽,張運(yùn)芝,王 瑩,賀 鑫,劉志邦

(1.遼寧省營口經(jīng)濟(jì)技術(shù)開發(fā)區(qū)氣象局,遼寧 營口 115007;2.遼寧省莊河市氣象局,遼寧 莊河 116400;3.遼寧省氣象裝備保障中心,遼寧 沈陽 110166;4.遼寧省遼陽縣氣象局,遼寧 遼陽 111200;5.遼寧省營口市氣象局,遼寧 營口 115007)

0 引言

玉米是遼寧省第一大種植作物,利用遙感資料準(zhǔn)確、及時(shí)地提取玉米種植面積信息,對提高玉米綜合管理水平、預(yù)測區(qū)域性糧食產(chǎn)量、優(yōu)化農(nóng)作物布局和確保糧食安全等具有重要意義[1-5]。

傳統(tǒng)的農(nóng)作物種植面積獲取采用農(nóng)業(yè)調(diào)查隊(duì)統(tǒng)計(jì)抽樣或參照行政級別逐級上報(bào)的方式,這兩種傳統(tǒng)面積獲取方式工作量大,人工成本高,耗費(fèi)時(shí)間長,且很難分析大尺度空間作物種植分布情況[6-7]。遙感提取作物種植面積研究發(fā)展的時(shí)間雖然不長,但其具有覆蓋區(qū)域廣、獲取頻率高、精確、快速、經(jīng)濟(jì)等優(yōu)勢而逐漸成為農(nóng)作物種植面積提取研究的重要手段[8]。目前農(nóng)作物種植面積遙感提取研究多集中在高光譜、高空間分辨率的數(shù)據(jù)源上[9-10],如Grace K等使用World View高分影像遙感提取了海地和危地馬拉作物的種植面積[11]; Mirik M等使用TM影像提取了德克薩斯州冬小麥的種植面積[12];李丹等使用Hyperion高光譜影像遙感獲取了廣州市荔枝的種植信息[13];余凌翔等使用HJ1/CCD影像遙感提取了西雙版納橡膠的分布信息[14]。高分辨率遙感數(shù)據(jù)具有較高的作物分類識別精度,但對于區(qū)域性作物種植面積提取而言,存在選取范圍內(nèi)遙感數(shù)據(jù)覆蓋不全、獲取難、成本高、時(shí)效差等諸多壁壘,中、低空間分辨率的多光譜數(shù)據(jù)在遙感提取區(qū)域性作物種植面積研究方面具有明顯優(yōu)勢[15-18]。常用于農(nóng)作物監(jiān)督分類、種植信息提取的植被指數(shù)有歸一化植被指數(shù)(Normalized difference vegetation index, NDVI)、比值植被指數(shù)(Ratio vegetation index, RVI)、垂直植被指數(shù)(Perpendicular vegetation index, PVI)、增強(qiáng)植被指數(shù)(Enhanced vegetation index, EVI)等[19]。其中NDVI植被指數(shù)能夠很好地反映出農(nóng)作物的生長信息,體現(xiàn)不同農(nóng)作物的物候特征,是目前衛(wèi)星遙感監(jiān)督分類作物應(yīng)用最穩(wěn)定、最普遍的參數(shù)之一[20-22]。使用擬合、濾波法重構(gòu)植被指數(shù)時(shí)序曲線提取農(nóng)作物種植面積是一種較為常用的方法。對于遼寧省春玉米而言,基于非對稱高斯函數(shù)擬合法重構(gòu)的NDVI時(shí)序曲線保真性能好,基于Savitzky-Golay濾波方法重構(gòu)的NDVI時(shí)序曲線平滑性能好,使用這兩種方法重構(gòu)的NDVI時(shí)序曲線都與玉米實(shí)際的生長狀況接近[23]。

因此,本研究選定遼寧省春玉米種植面積為研究對象,利用S-G濾波法和AG擬合法分別重構(gòu)春玉米NDVI時(shí)序曲線,選用動態(tài)振幅閾值法、拐點(diǎn)法和最大值法提取了遼寧省春玉米的關(guān)鍵物候期,并聯(lián)合決策樹分類法對遼寧省春玉米分布進(jìn)行了監(jiān)督分類,以期為使用衛(wèi)星遙感數(shù)據(jù)提取遼寧省春玉米種植面積提供一種有效的方法。

1 材料與方法

1.1 研究區(qū)域概況

遼寧省(東經(jīng)119°～126°，北緯39°～43°)位于中國東北南部,北與吉林省、內(nèi)蒙古自治區(qū)毗連,西與河北省相接,東有鴨綠江,南鄰渤海、黃海。遼寧省總面積為14.8萬km2,地勢為北高南低,東西兩側(cè)丘陵向中部平原下降。遼寧省氣候較為宜人,四季分明,屬溫帶大陸性季風(fēng)氣候區(qū),日照豐富,積溫偏高,溫差較大,降水不均,雨熱同期。

遼寧省是中國重要的玉米種植區(qū)之一,省內(nèi)絕大多數(shù)地區(qū)以種植春玉米為主,年均種植面積約為161.5萬hm2,約占全國總種植面積的6.47%;年均產(chǎn)量約為1069萬t,約占全國玉米總產(chǎn)量的8.61%[24]。

1.2 數(shù)據(jù)材料

1.2.1 遙感類數(shù)據(jù)來源 本研究使用2017年MODIS NDVI數(shù)據(jù)中的3級陸地標(biāo)準(zhǔn)數(shù)據(jù)產(chǎn)品MOD13Q1(MODIS/Terra Vegetation Indices 16-Day L3 Global 250 m SIN Grid),其時(shí)空分辨率為16 d、250 m,從美國國家航空航天局(NASA)官方網(wǎng)站下載。

1.2.2 統(tǒng)計(jì)類數(shù)據(jù)來源 選用遼寧省氣象局11個(gè)農(nóng)業(yè)氣象站(昌圖、沈北、建平、新民、燈塔、建昌、綏中、岫巖、寬甸、瓦房店、莊河，其分布見圖1)的觀測數(shù)據(jù)作為遼寧省春玉米物候期遙感提取的驗(yàn)證數(shù)據(jù),由遼寧省氣象信息中心提供;選用2018年遼寧省統(tǒng)計(jì)年鑒提供的統(tǒng)計(jì)數(shù)據(jù)作為春玉米種植面積遙感提取的驗(yàn)證數(shù)據(jù)。

圖1 遼寧省玉米實(shí)驗(yàn)區(qū)分布

1.3 研究方法

1.3.1 數(shù)據(jù)預(yù)處理 選用NASA研制的MODIS重投影工具(MRT)對2017年MOD13Q1數(shù)據(jù)產(chǎn)品進(jìn)行格式轉(zhuǎn)換、基準(zhǔn)投影、波段選擇和柵格拼接等預(yù)處理工作,其具體參數(shù)如表1所示。

表1 MRT投影參數(shù)設(shè)置

選用ENVI軟件和遼寧省行政邊界矢量對拼接后的數(shù)據(jù)進(jìn)行剪裁處理,得到2017年遼寧省春玉米發(fā)育期的NDVI柵格數(shù)據(jù)(圖2)。

圖2 遼寧省MODIS NDVI柵格數(shù)據(jù)(2017年5～9月)

以遼寧省氣象局11個(gè)農(nóng)業(yè)氣象站的地理坐標(biāo)制作美國標(biāo)準(zhǔn)信息交換代碼(ASCII)文本,使用ENVI Classic中的ROI Tool模塊提取所需的2017年5～9月實(shí)驗(yàn)研究點(diǎn)的NDVI柵格信息。

1.3.2 時(shí)間序列遙感數(shù)據(jù)的重構(gòu)方法 MODIS NDVI數(shù)據(jù)中的3級陸地標(biāo)準(zhǔn)數(shù)據(jù)產(chǎn)品的主要合成算法是角度歸一化合成方法,能夠有效地減少太陽高度角、暗影和殘?jiān)频挠绊?科學(xué)地提高合成數(shù)據(jù)的一致性和質(zhì)量,但該數(shù)據(jù)產(chǎn)品仍受水汽、臭氧、降水等因素的干擾,使其產(chǎn)生不可忽略的、不合理的波動[25-27]。因此,本研究使用TIMESAT軟件提供的以下2種重構(gòu)方法去除時(shí)序曲線的噪聲。

1.3.2.1 賽維迪克·格雷濾波方法(Savitzky-Golay, S-G) S-G濾波方法是移動平均算法的一種改良方法,是一種基于最小二乘卷積的移動平均算法,其權(quán)重由濾波窗口中做最小二乘多項(xiàng)式擬合而得出的次數(shù)決定[28]。濾波公式如下:

(1)

式(1)中:Y為NDVI原值；Y′為擬合值；Ci是第i個(gè)NDVI的權(quán)重值；N為卷積數(shù)目,數(shù)值上為濾波窗口寬度的2n+1。濾波效果由濾波窗口大小和平滑多項(xiàng)式的次數(shù)控制,濾波窗口影響了重構(gòu)后曲線的平滑能力及其斜率變化速率;平滑多項(xiàng)式次數(shù)的高低會使平滑結(jié)果出現(xiàn)過度平滑或平滑次數(shù)不足的現(xiàn)象[29]。

1.3.2.2 非對稱高斯函數(shù)擬合法(Asymmetric Gaussian model, AG) AG擬合法是基于高斯函數(shù)的分段最小二乘擬合算法,在運(yùn)算時(shí)函數(shù)先提取時(shí)序曲線內(nèi)的極大值和極小值,使用獲取的極值將曲線分為多個(gè)區(qū)間,分別在各自區(qū)間內(nèi)進(jìn)行局部擬合,再將各段的擬合曲線進(jìn)行整體拼接。局部擬合函數(shù)公式如下:

f(t)≡f(t;c,d)=c1+c2g(t;d)

(2)

式(2)中:線性系數(shù)c=(c1,c2)控制曲線的振幅和基值;非線性系數(shù)d=(d1,d2,…,di)控制基函數(shù)的形狀g(t;d)。

基函數(shù)公式如下:

(3)

式(3)中:d1為極大、極小值對應(yīng)變量t的位置參數(shù);d2,d3,d4和d5分別決定了極值左右半曲線的寬度和陡峭度。

1.3.3 參數(shù)設(shè)置 在重構(gòu)時(shí)序曲線時(shí),重構(gòu)模型的參數(shù)設(shè)置對重構(gòu)結(jié)果的優(yōu)劣有重要影響。參照前人的研究成果,本研究選用如下參數(shù):有效值域0.00～1.00,包絡(luò)迭代次數(shù)3,擬合峰值0.5,去噪聲模式2:STL replace,滑動窗口大小4,自適應(yīng)強(qiáng)度3.0。

1.3.4 玉米關(guān)鍵物候期的提取方法 使用Origin軟件中的4次擬合平滑功能對重構(gòu)后的曲線進(jìn)行平滑處理,選用動態(tài)振幅閾值法、曲線拐點(diǎn)法和最大值法提取遼寧省春玉米的關(guān)鍵物候期(出苗期、拔節(jié)期、抽雄期)。選用動態(tài)振幅閾值法提取春玉米的出苗期,以春玉米生長季開始前裸地的NDVI值為基值,以距離基值振幅10%的日序?yàn)榇河衩椎某雒缙凇＿x用拐點(diǎn)法提取春玉米的拔節(jié)期,在玉米營養(yǎng)生長過程中,其生長速率為“慢-快-慢”的“S”型曲線,且其拔節(jié)期生長速率最快,因此時(shí)序曲線上升階段的拐點(diǎn)即為春玉米的拔節(jié)期。選用最大值法提取春玉米的抽雄期,玉米抽雄期為從營養(yǎng)生長向生殖生長轉(zhuǎn)變的時(shí)期,此時(shí)玉米的平均葉面積最大,其NDVI時(shí)序曲線達(dá)到峰值,因此以NDVI最大值的日序?yàn)榇河衩椎某樾燮赱6]。

1.3.5 玉米種植面積的提取方法 選用關(guān)鍵期時(shí)序提取法獲取遼寧省春玉米的種植面積。以遙感反演的春玉米關(guān)鍵物候期(出苗期、拔節(jié)期、抽雄期)對應(yīng)日序臨近的2個(gè)NDVI數(shù)據(jù)作為遙感提取目標(biāo)(B1～B4),從遼寧省11個(gè)農(nóng)業(yè)氣象站重構(gòu)平滑后的NDVI數(shù)據(jù)集中檢索提取目標(biāo)的NDVI值，作為決策樹值域。使用ENVI軟件中的決策樹分類模塊對遙感目標(biāo)進(jìn)行分類提取,并通過比值計(jì)算獲得遼寧省春玉米的種植面積。決策樹監(jiān)督分類流程如圖3所示。

圖3 實(shí)驗(yàn)區(qū)春玉米監(jiān)督分類流程

選用極值法分別篩選玉米出苗期和拔節(jié)期的決策樹值域(T1和T2～T5)。在春玉米抽雄期易發(fā)生NDVI過飽和現(xiàn)象,因此選用四分位法篩選玉米抽雄期決策樹值域(T6、T7)。決策樹監(jiān)督分類閾值見表2。

表2 決策樹監(jiān)督分類的閾值

2 結(jié)果與分析

2.1 重構(gòu)方法的分析與評價(jià)

春玉米原NDVI時(shí)序曲線異常值波動明顯。本研究選用S-G和AG重構(gòu)法糾正時(shí)序曲線的異常值,使用四次擬合法進(jìn)行平滑處理,重構(gòu)平滑后的時(shí)序曲線(圖4)基本上保持原曲線的形態(tài),且除噪聲效果明顯。

以2017年綏中農(nóng)業(yè)氣象站為例。

本研究使用重構(gòu)曲線與原曲線的擬合度(R2)和均方根誤差(RMSE)來評價(jià)重構(gòu)方法的平滑能力和保真能力，結(jié)果如表3所示。從表3可以看出：基于S-G濾波法的各實(shí)驗(yàn)站點(diǎn)的平均擬合度為0.983，優(yōu)于AG重構(gòu)法的平均擬合度0.956,即S-G濾波法有更好的平滑能力;基于S-G濾波法和AG重構(gòu)法的各實(shí)驗(yàn)站點(diǎn)的平均RMSE分別為0.196和0.197,即這兩種重構(gòu)法的保真能力基本相當(dāng)?；趯ι鲜鰞煞N指標(biāo)的綜合評價(jià),用S-G濾波法重構(gòu)后的時(shí)序曲線綜合性能更優(yōu)。

表3 兩種NDVI時(shí)序曲線重構(gòu)方法的效果評估

2.2 春玉米關(guān)鍵物候期的對比分析

2.2.1 出苗期、拔節(jié)期和抽雄期的對比分析 分別對使用兩種重構(gòu)曲線提取的春玉米出苗期、拔節(jié)期和抽雄期日序與實(shí)際觀測物候期日序進(jìn)行比較，結(jié)果如圖5所示，基于S-G重構(gòu)法提取的各物候期精度普遍低于5 d,精度劣于5 d的樣本為6個(gè)(占18.1%),出苗期、拔節(jié)期與抽雄期的平均誤差分別為3.9、2.9和2.6 d;基于AG重構(gòu)法提取的各物候期精度普遍低于6 d,精度劣于6 d的樣本為4個(gè)(占12.1%),出苗期、拔節(jié)期與抽雄期的平均誤差分別為4.5、2.8和3.2 d。對于春玉米出苗期和抽雄期的提取效果而言,基于S-G重構(gòu)的時(shí)序曲線有更好的提取精度;對于春玉米拔節(jié)期的提取效果而言,基于AG重構(gòu)的時(shí)序曲線的提取精度略好。

圖5 2017年遼寧省春玉米關(guān)鍵物候期提取值與觀測值比較

2.2.2 不同地區(qū)關(guān)鍵物候期的對比分析 從表4可以看出，兩種重構(gòu)時(shí)序曲線在不同地區(qū)的提取精度有明顯差異,其中基于S-G重構(gòu)曲線的關(guān)鍵物候期提取精度在燈塔市實(shí)驗(yàn)區(qū)最優(yōu),其MAE為1.55 d;在寬甸縣實(shí)驗(yàn)區(qū)最劣,其MAE為3.96 d。基于AG重構(gòu)曲線的關(guān)鍵物候期的提取精度在新民縣實(shí)驗(yàn)區(qū)最優(yōu),其MAE為2.31 d;在建平縣實(shí)驗(yàn)區(qū)最劣,其MAE為5.00 d。基于各地的遙感提取效果,S-G重構(gòu)曲線的平均提取精度為3.11 d，優(yōu)于基于AG重構(gòu)曲線的平均提取精度3.45 d。

表4 基于S-G與AG重構(gòu)法的時(shí)序提取精度(MAE)對比 d

2.2.3 玉米種植面積提取分析 基于S-G重構(gòu)法提取的遼寧省2018年春玉米種植面積為2.54萬km2,提取精度為94.35%;基于AG重構(gòu)法提取的春玉米種植面積為2.4萬km2,提取精度為89.29%(圖6)。整體而言,基于S-G重構(gòu)曲線的春玉米種植面積提取精度優(yōu)于AG重構(gòu)曲線的面積提取精度。通過ENVI的三窗功能對比發(fā)現(xiàn),兩者提取差異主要集中在遼寧東、西部的山區(qū),說明相較于AG重構(gòu)曲線,S-G重構(gòu)曲線對海拔高、坡度陡、玉米地塊分散破碎、純凈度不高的像元有更好的提取效果。

a：S-G重構(gòu)法；b：AG重構(gòu)法。藍(lán)色區(qū)域?yàn)橛衩追N植分布。

3 結(jié)論與討論

本研究從春玉米物候信息遙感研究的特殊性出發(fā),使用2017年4～9月遼寧省MOD13Q1數(shù)據(jù)產(chǎn)品中的NDVI波段構(gòu)建時(shí)序曲線,采用S-G重構(gòu)法、AG重構(gòu)法和四次擬合平滑法除噪濾波平滑曲線,應(yīng)用動態(tài)振幅閾值法、拐點(diǎn)法和最大值法提取春玉米關(guān)鍵物候期,并聯(lián)合決策樹監(jiān)督分類法估算了遼寧省春玉米種植面積。研究結(jié)果表明:S-G重構(gòu)法對時(shí)序曲線的平滑能力和綜合重構(gòu)性能更優(yōu),S-G重構(gòu)法和AG重構(gòu)法對時(shí)序曲線的保真能力基本相當(dāng)。基于S-G重構(gòu)法提取的出苗期和抽雄期優(yōu)于基于AG重構(gòu)法的提取結(jié)果,基于AG重構(gòu)法提取的拔節(jié)期略優(yōu)于S-G重構(gòu)法的提取結(jié)果?；赟-G重構(gòu)法對各實(shí)驗(yàn)區(qū)關(guān)鍵物候期的整體提取效果好于基于AG重構(gòu)法的整體提取效果?；赟-G重構(gòu)法的遼寧省春玉米種植面積的提取精度為94.35%,優(yōu)于基于AG重構(gòu)法的面積提取精度89.29%。整體而言,S-G重構(gòu)時(shí)序曲線的綜合性能更優(yōu),對關(guān)鍵物候期的提取效果更好,對春玉米種植面積的提取精度更高,因此聯(lián)合S-G時(shí)序重構(gòu)法和決策樹監(jiān)督分類法可為省級遙感估算春玉米種植面積提供一種新的途徑。

本研究在以下3個(gè)方面還需繼續(xù)改進(jìn):一是春玉米監(jiān)督分類精度有待提高。遼寧省中部為平原地區(qū),春玉米種植基本連片,像元的純凈度較高,提取效果好；但遼寧省東西兩側(cè)為丘陵,春玉米種植破碎,像元易混入雜物,從而降低了春玉米種植面積的提取精度。今后的研究需采用多元數(shù)據(jù)融合技術(shù),以提高融合后的MODIS復(fù)合數(shù)據(jù)影像的空間分辨率,獲取更高精度的春玉米監(jiān)督分類,減少同物異譜。二是春玉米提取面積區(qū)域需要精細(xì)化劃分。由于省域提取面積范圍較大,受地理位置、環(huán)境氣候等影響,各地區(qū)春玉米的種植品種、物候和長勢存在一定的差異。本研究受數(shù)據(jù)資源的限制,沒有分區(qū)域精細(xì)化遙感提取種植面積,導(dǎo)致部分縣(區(qū))玉米種植面積的提取精度下降。因此，今后的研究需構(gòu)建省內(nèi)精細(xì)化的種植面積提取區(qū)域,優(yōu)化種植面積提取模型。三是數(shù)據(jù)產(chǎn)品的內(nèi)在性能有待優(yōu)化。由于NDVI數(shù)據(jù)產(chǎn)品易出現(xiàn)飽和現(xiàn)象,致使春玉米成熟期附近的監(jiān)測效果下降,因此今后的研究需通過對比分析增強(qiáng)型植被指數(shù)(EVI)數(shù)據(jù)產(chǎn)品來克服NDVI數(shù)據(jù)產(chǎn)品的不足。