王濤，劉立波，張鵬*，王曉麗

1.寧夏大學(xué)信息工程學(xué)院，銀川 750021

2.中國(guó)農(nóng)業(yè)科學(xué)院農(nóng)業(yè)信息研究所，北京 100081

3.國(guó)家農(nóng)業(yè)科學(xué)數(shù)據(jù)中心，北京 100081

4.中國(guó)農(nóng)業(yè)科學(xué)院國(guó)家南繁研究院，海南三亞 572024

引 言

估產(chǎn)數(shù)據(jù)集是表征農(nóng)作物長(zhǎng)勢(shì)和產(chǎn)量的重要信息，被廣泛用于作物長(zhǎng)勢(shì)分析、產(chǎn)量預(yù)測(cè)等研究。然而，傳統(tǒng)的估產(chǎn)數(shù)據(jù)集制作方法采用人工實(shí)地測(cè)量土壤墑情、作物干重、氣象條件等參數(shù)，從而構(gòu)建數(shù)據(jù)集[1-3]。該方法速度慢、工作量大，且過度依賴于實(shí)測(cè)數(shù)據(jù)。

近年來，遙感影像因其低成本、覆蓋面廣等優(yōu)勢(shì)，常被用于提取歸一化植被指數(shù)（Normalized Difference Vegetation Index，NDVI）和增強(qiáng)型植被指數(shù)（Enhanced Vegetation Index，EVI）。如 Quarmby[4-6]等利用紅外波段和近紅外波段計(jì)算得到了NDVI，Bolton等[7-8]利用紅外、近紅外和藍(lán)光波段計(jì)算得到了EVI；此外，滿衛(wèi)東等[9]基于AVHRR影像制作了基于遼寧省1982至2009年AVHRR歸一化植被指數(shù)數(shù)據(jù)集（https://www.osgeo.cn/data/wef57）。更進(jìn)一步，Aerial Intelligence公司2017年發(fā)布的美國(guó)冬小麥產(chǎn)量預(yù)測(cè)數(shù)據(jù)集（https://aerialintel.blob.core.windows.net/recruiting/datasets/wheat-2013-supervised.csv 和 https://aerialintel.blob.core.windows.net/recruiting/datasets/wheat-2014-supervised.csv）不僅提取遙感影像中的NDVI和EVI，而且加入了氣象、地理位置等特征，提升了估產(chǎn)特征的多樣性。但在這類數(shù)據(jù)集中，NDVI、EVI等手工特征僅用少許的波段計(jì)算得到，忽略了其余波段重要信息，且特征的選取過度依賴人工經(jīng)驗(yàn)，具有一定局限性。

因此，針對(duì)以上問題，本文摒棄手工特征制作的方法，制作了一種基于MODIS高光譜遙感影像的、多波段和多時(shí)相融合的寧夏縣級(jí)尺度枸杞估產(chǎn)數(shù)據(jù)集，用于卷積神經(jīng)網(wǎng)絡(luò)特征的自動(dòng)提取，簡(jiǎn)化特征提取操作的同時(shí)，進(jìn)一步增強(qiáng)了特征豐富度。

1 數(shù)據(jù)采集和處理方法

1.1 數(shù)據(jù)來源

本文采用的實(shí)驗(yàn)數(shù)據(jù)由寧夏回族自治區(qū)遙感影像、枸杞種植區(qū)域矢量圖和年際枸杞產(chǎn)量3類數(shù)據(jù)組成。其中，遙感影像采用 MODIS高光譜影像數(shù)據(jù)，來源于 EARTHDATA網(wǎng)站（https://search.earthdata.nasa.gov/），行列號(hào)為h26v04，時(shí)間范圍為枸杞生長(zhǎng)季內(nèi)的每年第97天至第297天，包括MOD09A1、MOD13A1、MYD11A2、MCD15A2H 4種類型的MODIS產(chǎn)品數(shù)據(jù)（表1）；枸杞種植區(qū)域矢量圖由寧夏農(nóng)林科學(xué)院研究人員實(shí)地記錄枸杞種植區(qū)域經(jīng)緯度制成，為shp文件；年際枸杞產(chǎn)量數(shù)據(jù)來源于寧夏回族自治區(qū)統(tǒng)計(jì)局，包括2010-2019年寧夏16縣（縣級(jí)市區(qū)）枸杞種植面積和實(shí)際產(chǎn)量。其中，2019年產(chǎn)量數(shù)據(jù)如表2所示，枸杞種植總面積為27960公頃，總產(chǎn)量為94843噸，平均產(chǎn)量為3.39噸/公頃。其中，同心縣、紅寺堡區(qū)等地區(qū)枸杞產(chǎn)量均高于9000噸，屬于高產(chǎn)地區(qū)；金鳳區(qū)、平羅縣等地區(qū)產(chǎn)量均在1000噸左右，屬于中產(chǎn)地區(qū)；其余為低產(chǎn)地區(qū)。

表1 影像及矢量數(shù)據(jù)表Table 1 The table of remote sensing images and vector data

表2 2019年產(chǎn)量數(shù)據(jù)表Table 2 The table of yield data in 2019

1.2 數(shù)據(jù)處理

1.2.1 重投影與重采樣

為了保證遙感影像數(shù)據(jù)空間位置的一致性，首先，利用 MODIS影像重投影工具（MODIS Reprojection Tool， MRT）工具將MODIS遙感影像和枸杞種植區(qū)域掩膜數(shù)據(jù)重投影為基于WGS-84橢球體的UTM投影；然后，將MYD11A2影像和枸杞矢量數(shù)據(jù)重采樣為500 m，使MYD11A2影像和枸杞矢量數(shù)據(jù)與其余MODIS數(shù)據(jù)產(chǎn)品的空間分辨率相互統(tǒng)一。

1.2.2 時(shí)間序列補(bǔ)充

由于MOD13A1的時(shí)間分辨率為16天，其余MODIS數(shù)據(jù)產(chǎn)品為8天（表1），為了保證時(shí)間序列的完整性，根據(jù)式（1）采用上下影像求平均的方法對(duì)枸杞生長(zhǎng)季內(nèi)缺失影像進(jìn)行補(bǔ)充，如將MOD13A1第97天和第113天影像的均值作為第105天的影像數(shù)據(jù)。

式中，Ii、Ii-8和Ii+8分別為MOD13A1第i天、第i-8天和第i+8天影像數(shù)據(jù)，i的取值范圍為[105, 289]，時(shí)間間隔為8天。

1.2.3 波段融合和時(shí)間序列融合

為了提高波段信息的豐富度，分別提取了同一景MOD09A1、MOD13A1等遙感影像中的band1-band7、NDVI、EVI等13個(gè)波段（表1），并對(duì)其進(jìn)行了融合。隨后，為了進(jìn)一步融合影像不同時(shí)相上的枸杞生長(zhǎng)信息，對(duì)波段融合結(jié)果在時(shí)間維度上進(jìn)行融合，最終形成了10張波段數(shù)為13，時(shí)相為26的時(shí)間序列影像。

1.2.4 枸杞種植區(qū)域提取

遙感圖像是基于像素點(diǎn)的，包括了各種地面覆蓋類型，為了使各縣年鑒統(tǒng)計(jì)產(chǎn)量與其遙感影像相對(duì)應(yīng)，利用了不同縣域枸杞種植矢量圖提取上述融合數(shù)據(jù)中的種植區(qū)域，以中寧縣枸杞種植區(qū)域提取為例。首先，采用ArcGIS軟件從寧夏行政區(qū)劃矢量圖提取出中寧縣行政區(qū)劃矢量圖；然后，根據(jù)中寧縣行政區(qū)劃矢量圖左上角和右下角經(jīng)緯度裁剪得到中寧縣融合影像（圖1）；接著，調(diào)用GDAL（Geospatial Data Abstraction Library）庫(kù)中的warp函數(shù)以實(shí)現(xiàn)枸杞中寧枸杞種植區(qū)域矢量圖（圖2）對(duì)圖1的裁剪，從而得到中寧市枸杞種植區(qū)域圖（圖3）；最后，以此類推，提取出其余縣枸杞種植區(qū)域。

1.2.5 直方圖降維

鑒于影像數(shù)據(jù)集的稀疏性，無(wú)法采用端到端方式訓(xùn)練估產(chǎn)模型。因此，本文將影像無(wú)差別劃分為32個(gè)像素區(qū)間，進(jìn)而將影像中每個(gè)波段不同像素值映射至不同區(qū)間，以達(dá)到直方圖降維的目的。

2 數(shù)據(jù)樣本描述

基于MODIS的寧夏縣級(jí)尺度枸杞估產(chǎn)數(shù)據(jù)集（2010-2019）主要由遙感影像數(shù)據(jù)和年鑒統(tǒng)計(jì)產(chǎn)量數(shù)據(jù)組成，具體數(shù)據(jù)樣本描述如下。

圖1 中寧縣高光譜影像Figure 1 Hyperspectral image of Zhongning County

圖2 中寧縣枸杞種植區(qū)域矢量圖Figure 2 Vector map of wolfberry planting area in Zhongning County

圖3 中寧縣枸杞種植區(qū)域圖Figure 3 The area of interest for wolfberry planting in Zhongning County

2.1 遙感影像數(shù)據(jù)描述

遙感影像數(shù)據(jù)以MODI高光譜影像為數(shù)據(jù)源，經(jīng)上述方法處理后，共形成了160個(gè)大小為32×26×13（32為像素區(qū)間數(shù)，26為時(shí)相數(shù)，13為波段數(shù)）的直方圖矩陣，以“hist_年份_地名.npy”的方式進(jìn)行命名，存放于“dataset/hist_data”文件夾中。其中，16個(gè)縣（縣級(jí)市區(qū)）（表2）中的每個(gè)地區(qū)各涵蓋10年直方圖矩陣數(shù)據(jù)，時(shí)間范圍為2010-2019年，直方圖矩陣數(shù)據(jù)示意如圖4所示。

圖4 直方圖矩陣數(shù)據(jù)示意圖Figure 4 Schematic diagram of histogram matrix data

2.2 統(tǒng)計(jì)產(chǎn)量數(shù)據(jù)描述

為了方便用戶對(duì)數(shù)據(jù)的加載及使用，本文提供了與直方圖矩陣年份和地名相對(duì)應(yīng)一致的年鑒統(tǒng)計(jì)產(chǎn)量數(shù)據(jù)。該數(shù)據(jù)為表格數(shù)據(jù)，保存于“dataset/yield_data/ yield_data.csv”文件中，共包含7個(gè)屬性列，分別為年份、省份名、省份編號(hào)、縣名、縣級(jí)編號(hào)、枸杞種植面積和產(chǎn)量，總計(jì)160條產(chǎn)量數(shù)據(jù)，平均每個(gè)縣（縣級(jí)市區(qū)）包含10條2010-2019年的枸杞產(chǎn)量數(shù)據(jù)，部分產(chǎn)量數(shù)據(jù)如表3所示。

表3 部分產(chǎn)量數(shù)據(jù)表Table 3 The table of partial yield data

3 數(shù)據(jù)質(zhì)量控制和評(píng)估

為了驗(yàn)證本文數(shù)據(jù)集的可用性，本文分別從定性和定量?jī)蓚€(gè)角度評(píng)價(jià)本數(shù)據(jù)質(zhì)量和準(zhǔn)確性，具體如下：

3.1 定性評(píng)價(jià)

通過可視化對(duì)比高中低產(chǎn)量所對(duì)應(yīng)的直方圖降維結(jié)果（圖5）可見，在高產(chǎn)、中產(chǎn)和低產(chǎn)中波段1、波段7、NDVI和EVI波段明顯存在視覺差異，表明可以從本數(shù)據(jù)集中提取到有利的特征，用于表征枸杞的長(zhǎng)勢(shì)和產(chǎn)量。此外，在時(shí)間維度上，像素區(qū)間最大值基本集中于第180天左右，正值枸杞植株長(zhǎng)勢(shì)茂盛期，符合枸杞的生長(zhǎng)規(guī)律，說明本文數(shù)據(jù)集質(zhì)量較好。

圖5 可視化結(jié)果圖Figure 5 Visualized result graph

3.2 定量評(píng)價(jià)

為了進(jìn)一步說明該數(shù)據(jù)集的準(zhǔn)確性，本文以平均相對(duì)誤差（Mean Relative Error, MRE）、均方根誤差（Root Mean Square Error, RMSE）和決定系數(shù)（Coefficient of Determination, R2）為評(píng)價(jià)指標(biāo)，計(jì)算公式如式（2）-式（4）所示，分別基于本文數(shù)據(jù)集和Aerial Intelligence公司發(fā)布的美國(guó)冬小麥產(chǎn)量預(yù)測(cè)數(shù)據(jù)集，設(shè)置了兩個(gè)對(duì)比實(shí)驗(yàn)組，分別是：

（1）以本文數(shù)據(jù)集為數(shù)據(jù)源，構(gòu)建了卷積神經(jīng)網(wǎng)絡(luò)（Convolutional Neural Networks, CNN），用于枸杞長(zhǎng)勢(shì)特征的抽取，進(jìn)而采用全連接網(wǎng)絡(luò)（Fully Connected Network, FCN）、支持向量回歸（Support Vector Regression, SVR）和嶺回歸（Ridge Regression, RR）方法預(yù)測(cè)枸杞年際產(chǎn)量。卷積神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)如圖6所示，第一層為輸入層，為本文數(shù)據(jù)集，訓(xùn)練集、驗(yàn)證集和測(cè)試集劃分比率為7:2:1；第二層為特征提取層，由 6個(gè)卷積層構(gòu)成，卷積核數(shù)量分別是 128、256、256、512、512、512，卷積核大小全為3×3；最后一層為全連接層，將特征圖映射至2048維embedding（嵌入向量）空間，以表征枸杞的長(zhǎng)勢(shì)信息。

（2）與上述方法相比，對(duì)比實(shí)驗(yàn)基于冬小麥產(chǎn)量預(yù)測(cè)數(shù)據(jù)集，僅用單一的FCN、RR和SVR預(yù)測(cè)冬小麥產(chǎn)量。其中，F(xiàn)CN包含5個(gè)隱藏層，維度分別為1024、512、256、128和1，實(shí)驗(yàn)結(jié)果如表4所示。

表4 對(duì)比實(shí)驗(yàn)結(jié)果表Table 4 The table of comparative experiment results

由表4可知，在本文數(shù)據(jù)集上，MRE和RMSE分別為14.52%、859.23噸，且R2達(dá)到了0.83，均優(yōu)于對(duì)比數(shù)據(jù)集。與2017年冬小麥產(chǎn)量預(yù)測(cè)比賽數(shù)據(jù)集相比，本文數(shù)據(jù)集采用波段融合和時(shí)間序列融合方法分別融合了不同波段和時(shí)相上的影像信息，增強(qiáng)了特征豐富度，所以基于本文數(shù)據(jù)集的CNN回歸方法，在 MRE和 RMSE上分別下降了1.48%和89.14噸，R2上升了0.06，說明增加了0.06%的特征可由回歸方法解釋，驗(yàn)證了本文數(shù)據(jù)集的準(zhǔn)確性。

其中， Ti為第i個(gè)縣枸杞統(tǒng)計(jì)產(chǎn)量，Pi為第i個(gè)縣產(chǎn)量預(yù)測(cè)值，為16縣統(tǒng)計(jì)數(shù)據(jù)平均產(chǎn)量，單位為噸。

圖6 卷積神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖Figure 6 Convolutional neural network structure diagram

4 數(shù)據(jù)價(jià)值

目前，公開的遙感估產(chǎn)數(shù)據(jù)集非常有限，且目視解譯紋理特征等相關(guān)提法過度依賴于實(shí)測(cè)數(shù)據(jù)和人工經(jīng)驗(yàn)。因此，本文構(gòu)建了一種基于MODIS的寧夏縣級(jí)尺度枸杞估產(chǎn)數(shù)據(jù)集，具有高光譜、多波段、多時(shí)相等特點(diǎn)，可用于卷積神經(jīng)網(wǎng)絡(luò)特征自動(dòng)提取農(nóng)作物長(zhǎng)勢(shì)特征，為作物產(chǎn)量預(yù)測(cè)研究提供數(shù)據(jù)支撐，同時(shí)為遙感多時(shí)相估產(chǎn)影像數(shù)據(jù)集的制作提供了一定的參考價(jià)值。