山東省冬小麥單產(chǎn)監(jiān)測與預報方法研究

郭 銳 朱秀芳 李石波 侯陳瑤

(1.北京師范大學地表過程與資源生態(tài)國家重點實驗室， 北京 100875；2.北京師范大學地理科學學部遙感科學與工程研究院， 北京 100875；3.中國地質(zhì)大學土地科學技術學院， 北京 100083)

0 引言

遙感農(nóng)業(yè)監(jiān)測具有觀測范圍大、時間分辨率高、數(shù)據(jù)客觀可靠的優(yōu)點，比工作量大、成本高、效率低的傳統(tǒng)估產(chǎn)方式節(jié)省了時間和成本，為農(nóng)作物估產(chǎn)提供了科學有效的手段[1-2]。

利用遙感進行作物產(chǎn)量估算的方法主要包括：遙感統(tǒng)計估產(chǎn)模型、干物質(zhì)-產(chǎn)量模型和作物模型模擬。遙感統(tǒng)計估產(chǎn)模型通過建立遙感變量與產(chǎn)量之間的關系表達式來進行產(chǎn)量估算[3-5]；干物質(zhì)-產(chǎn)量模型先基于遙感數(shù)據(jù)估算作物的地上生物量，再通過收獲指數(shù)轉(zhuǎn)換成作物的經(jīng)濟產(chǎn)量[6-7]；遙感作物模型模擬是將遙感數(shù)據(jù)作為模型校正的數(shù)據(jù)源之一，對作物模型進行參數(shù)本地化后，在氣象、土壤、作物種植信息等數(shù)據(jù)的驅(qū)動下進行作物生長模擬和產(chǎn)量的估算[8-10]。3種方法中第1種方法最簡單，對數(shù)據(jù)的要求最低，后兩種方法機理性更強，要求大量輸入數(shù)據(jù)，操作更復雜。

在遙感統(tǒng)計估產(chǎn)模型中，使用最多的輸入變量是植被指數(shù)。植被指數(shù)可以反映植被的生產(chǎn)力和健康狀況，研究表明，植被指數(shù)與作物產(chǎn)量之間高度相關[11-12]。以往研究中用到的植被指數(shù)包括歸一化植被指數(shù)(NDVI)[13-15]、增強型植被指數(shù)(EVI)[16]、葉面積指數(shù)(LAI)[17]、垂直植被指數(shù)(PVI)[18]、植被條件指數(shù)(VCI)[19]、光合有效輻射(APAR)[20]等。盡管利用植被構建遙感估產(chǎn)模型取得了成功，但也存在一些問題，比如NDVI在植被高覆蓋區(qū)容易飽和，對災害(如干旱、病蟲害)的響應滯后，不能很好地反映農(nóng)業(yè)管理和科技進步帶來的產(chǎn)量增加趨勢[21]。

另外，隨著遙感數(shù)據(jù)產(chǎn)品的進一步豐富，發(fā)布了大量環(huán)境變量遙感產(chǎn)品(如溫度、土壤濕度、降水量和蒸散發(fā)量)，有研究開始綜合使用環(huán)境變量和植被指數(shù)遙感產(chǎn)品建立統(tǒng)計估產(chǎn)模型。例如， PRASAD等[22]利用NDVI、VCI和條件溫度指數(shù)(TCI)等遙感數(shù)據(jù)對植被的健康狀況和土壤生產(chǎn)力進行評估，并構建了作物估產(chǎn)模型。DABROWSKA-ZIELINSKA等[23]利用早春和初夏兩個時期的VCI和TCI進行波蘭作物生長狀況的評估與作物產(chǎn)量的估計。上述研究都表明，加入環(huán)境變量對作物估產(chǎn)具有積極作用。

本文以山東省冬小麥為例，綜合使用技術產(chǎn)量、植被指數(shù)和環(huán)境變量構建統(tǒng)計估產(chǎn)模型，同時在監(jiān)測模式和預報模式下進行模型的應用和驗證，以期為冬小麥實時產(chǎn)量預報提供方法參考。

1 研究區(qū)及數(shù)據(jù)獲取

1.1 研究區(qū)概況

研究區(qū)域為山東省(圖1)，屬暖溫帶季風氣候，氣候條件穩(wěn)定，冬季寒冷干燥，夏季高溫多雨，光照資源豐富，適合農(nóng)業(yè)種植，主要的種植作物有冬小麥、棉花、玉米。山東省的小麥種植面積常年保持在20萬hm2，小麥總產(chǎn)量約占全國小麥產(chǎn)量的18%，是我國僅次于河南省的小麥生產(chǎn)省。干旱是影響山東省冬小麥產(chǎn)量的一大因素，但隨著種植方法與灌溉技術的不斷進步，山東省冬小麥抵抗干旱、病蟲害等災害的能力逐漸增強，整體小麥產(chǎn)量呈現(xiàn)穩(wěn)中有升的趨勢。

圖1 研究區(qū)域圖Fig.1 Study area map

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

本文使用的遙感數(shù)據(jù)為500 m分辨率地表反射率(Surface reflectance)8 d合成產(chǎn)品MOD09A1和500 m分辨率全球陸地蒸發(fā)蒸騰(Global terrestrial evapotranspiration)8 d合成產(chǎn)品MOD16A2，均下載自LAADS DAAC(https:∥ladsweb.modaps.eosdis.nasa.gov/)，取水平第27景、垂直第5景數(shù)據(jù)，位于102.2°～128.8°E、29.3°～40.0°N之間，得到包括河南、山東等省份在內(nèi)的影像。數(shù)據(jù)時間為2007—2017年，時間間隔為8 d/景，共計308幅影像。

MOD09A1產(chǎn)品提供了MODIS傳感器500 m分辨率1～7波段表面反射率8 d合成的數(shù)據(jù)產(chǎn)品，投影為正弦曲線投影，數(shù)據(jù)集中的每個像素包含了8 d內(nèi)盡可能準確的觀測值。由于MODIS數(shù)據(jù)產(chǎn)品中未包含500 m分辨率的8 d合成EVI數(shù)據(jù)，為了統(tǒng)一遙感數(shù)據(jù)的時間與空間分辨率，本文基于MOD09A1數(shù)據(jù)計算得到研究中所使用的500 m分辨率的8 d合成增強型植被指數(shù)(EVI)數(shù)據(jù)。

MOD16A2產(chǎn)品數(shù)據(jù)集中包含了實際蒸騰和潛在蒸騰，ET和PET表示500 m分辨率8 d內(nèi)單位面積通過蒸騰散失水分(0.012 5 kg/(m2·d))的總和。其中ET表示不同植被覆蓋度條件下植被區(qū)域和非植被區(qū)域的加權平均蒸騰量；PET表示在假設水分供應不受限制的情況下，某一固定下墊面可能達到的最大蒸騰量。

使用的非遙感數(shù)據(jù)包括中國多時期土地利用土地覆被遙感監(jiān)測數(shù)據(jù)集(CNLUCC)[24]、山東省矢量邊界數(shù)據(jù)、山東省作物物候觀測數(shù)據(jù)和山東省歷史產(chǎn)量統(tǒng)計數(shù)據(jù)等。使用的耕地掩膜數(shù)據(jù)由中國多時期土地利用土地覆被遙感監(jiān)測數(shù)據(jù)集(CNLUCC)中一級類型為耕地的像元(一級類型編號為1)提取并二值化后獲得。數(shù)據(jù)的預處理主要包括遙感影像的投影轉(zhuǎn)換、數(shù)據(jù)裁剪、像元篩選等。

2 研究方案與方法

2.1 研究方案

本文的研究方案包括：①計算增強型植被指數(shù)(EVI)和作物水分脅迫指數(shù)(CWSI)。利用山東省2007—2017年冬小麥生育期內(nèi)的MODIS傳感器地表反射率產(chǎn)品(MOD09A1)計算各縣區(qū)平均累計增強型植被指數(shù)(EVI), 利用同時期蒸騰散發(fā)產(chǎn)品(MOD16A2)計算各縣區(qū)平均累計水分脅迫指數(shù)(CWSI)。②確定技術產(chǎn)量因子。技術產(chǎn)量反映了該地區(qū)的平均農(nóng)業(yè)生產(chǎn)技術水平以及抗災能力。結合山東省2007—2017年縣級歷史產(chǎn)量統(tǒng)計數(shù)據(jù)，利用時間序列趨勢方法計算技術產(chǎn)量。③建立估產(chǎn)模型。使用最小二乘線性回歸法分別進行山東省市級尺度和省級尺度的冬小麥產(chǎn)量預測模型的構建。④模型的應用及精度驗證。分別在監(jiān)測模式和預報模式下進行估產(chǎn)模型的應用及精度驗證。技術路線見圖2。

圖2 技術路線圖Fig.2 Technical flowchart

2.2 平均累計指數(shù)計算

2.2.1平均作物水分脅迫指數(shù)(CWSI)

作物水分脅迫指數(shù)(CWSI)由500 m分辨率全球陸地蒸發(fā)蒸騰8 d合成產(chǎn)品MOD16A2計算得到。作物水分脅迫指數(shù)可以反映植被不同生長狀況下蒸騰量的變化和生長環(huán)境的干旱程度[25]，計算式為

(1)

式中CWSI——作物水分脅迫指數(shù)

ET——實際蒸騰量

PET——潛在蒸騰量

縣級平均CWSI具體計算過程如下：首先，利用耕地掩膜數(shù)據(jù)對CWSI數(shù)據(jù)進行掩膜處理，得到山東省耕地范圍內(nèi)的CWSI數(shù)據(jù)。之后，利用山東省縣級行政邊界矢量數(shù)據(jù)，分別提取每個縣區(qū)對應的共計11年、每年14期的所有耕地像元的CWSI平均值，將求得的平均值作為該縣區(qū)該年該期的脅迫指數(shù)。最后，統(tǒng)計縣級生育期內(nèi)的累計平均脅迫指數(shù)，對缺失3期及3期以下的數(shù)據(jù)(缺失期數(shù)不超過總期數(shù)的20%)使用該縣區(qū)其他年份的多年平均值作為替代，構成完整的生育期脅迫指數(shù)數(shù)據(jù)；對于缺失3期及3期以上的數(shù)據(jù)(缺失期數(shù)超過總期數(shù)的20%)，則判定該地區(qū)該年的數(shù)據(jù)缺失過多，舍棄該年的數(shù)據(jù)，不參與后續(xù)的估產(chǎn)建模。最后求得共135個縣區(qū)、每個縣區(qū)共11年(除數(shù)據(jù)缺失年份)的生育期內(nèi)的累計平均脅迫指數(shù)。

2.2.2平均增強型植被指數(shù)(EVI)

之前有研究對比了NDVI、LAI、EVI、GPP(總初級生產(chǎn)力)和LST(地表溫度)等變量對作物產(chǎn)量估計的影響，指出EVI在多種作物估產(chǎn)的應用中效果最優(yōu)[11]。因此，本研究中利用500 m分辨率表面反射率8 d合成的數(shù)據(jù)產(chǎn)品(MOD09A1)進行了EVI的計算[26]，計算式為

(2)

式中EVI——增強型植被指數(shù)

ρRED——紅波段反射率

G——增益參數(shù)，取2.5

ρNIR——近紅外波段反射率

ρBLUE——藍波段反射率

L——增益參數(shù)，取1

C1——大氣修正紅光校正參數(shù)，取6.0

C2——大氣修正藍光校正參數(shù)，取7.5

然后使用TIMESAT軟件對EVI時間序列進行平滑，去除圖像噪聲[27-28]，計算縣級平均EVI和生育期內(nèi)的累計平均EVI。其計算過程同2.2.1節(jié)中縣級累計平均CWSI的計算過程。

2.3 技術產(chǎn)量計算

技術產(chǎn)量代表研究區(qū)域長期的社會生產(chǎn)力發(fā)展水平所決定的產(chǎn)量分量。本文使用時間序列趨勢分析的算法求取技術產(chǎn)量[29]。該方法將整個時間序列內(nèi)的歷史產(chǎn)量，在某個與滑動步長時間內(nèi)進行線性擬合，形成一條線性函數(shù)的直線。隨著滑動直線不斷向后移動，不斷生成新的擬合直線。在直線滑動完成后，各時間點上均對應有大于或等于1個直線的模擬值，再對各時間點上的模擬值求平均值，即得到技術產(chǎn)量。這種模擬方法既不損失樣本序列的年數(shù)，也避免了主觀假定長時間序列產(chǎn)量變化的曲線類型，是一種較為實用的趨勢模擬方法。

2.4 估產(chǎn)模型構建

以2017年之前縣級數(shù)據(jù)作為建模樣本輸入，基于最小二乘法的多元線性回歸方法分別建立了省級和市級估產(chǎn)模型。省級尺度上的估產(chǎn)可以反映山東省整體的農(nóng)作物生產(chǎn)力，對于大范圍監(jiān)測冬小麥長勢與實時估算小麥產(chǎn)量有一定的理論意義。市級估產(chǎn)可以反映不同市域范圍內(nèi)由于農(nóng)作物生產(chǎn)技術、氣候和土壤條件所導致的產(chǎn)量差異，對于各市政府針對性地進行科學生產(chǎn)管理與農(nóng)業(yè)政策制定有積極的作用。

對于山東省市級冬小麥單產(chǎn)估產(chǎn)模型，去除了樣本點低于50個的城市(威海市、日照市、萊蕪市)以及遙感數(shù)據(jù)有嚴重缺失的城市(濱州市)，最終保留了共13個城市的樣本數(shù)據(jù)，完成山東省市級小麥單產(chǎn)估產(chǎn)建模，得到13個市級線性回歸方程。

對于山東省省級冬小麥單產(chǎn)估產(chǎn)模型，首先對全省所有縣區(qū)的數(shù)據(jù)樣本進行篩選，去除數(shù)據(jù)缺失和存在異常值的樣本點，最終保留736個樣本，完成山東省省級小麥單產(chǎn)估產(chǎn)建模，得到1個省級線性回歸方程。

2.5 模型應用與精度驗證

以2017年山東省小麥產(chǎn)量數(shù)據(jù)為驗證樣本，采用監(jiān)測和預報兩種模式進行估產(chǎn)模型的應用與驗證。監(jiān)測模式面向生長季結束后的最終產(chǎn)量估算，需在小麥生長季結束后獲得全生長季完整的遙感數(shù)據(jù)后才能進行；而預報模式是在小麥生長季節(jié)開始后進行實時的產(chǎn)量預測。

對于監(jiān)測模式，首先對各市和全省范圍內(nèi)所包含的縣所對應的2017年縣級累計平均EVI、縣級累計平均CWSI和技術產(chǎn)量數(shù)據(jù)分別求取市級和省級的平均值，以2017年的生長季內(nèi)(第65天到第169天)市級累計平均EVI、市級累計平均CWSI和估算的2017年的各市的技術產(chǎn)量數(shù)據(jù)為輸入，分別代入對應的市級估產(chǎn)模型中進行各市的冬小麥單產(chǎn)估算，將省級累計平均EVI、累計平均CWSI和技術產(chǎn)量代入省級估產(chǎn)模型進行全省的冬小麥的單產(chǎn)估算。

預報模式在小麥生長季節(jié)開始時就進行預測。然而在生長季初期，模型所需的整個生長季的累計平均EVI和CWSI是未知的。為此，對于未知時間段的EVI和CWSI采用多年平均值進行代替。例如，生長季一開始(第65天)就進行產(chǎn)量監(jiān)測時，第65天的EVI和CWSI用2017年的實測數(shù)據(jù)，而后的第73天到第165天用歷史平均值替代，最終得到完整的生長季的累計平均EVI和CWSI。隨著季節(jié)的推進，當前生長季越來越多的觀測值被納入模型，預報模型的結果會越來越接近監(jiān)測模式的結果。本文在生長季中選取3個時間點進行預報模式下的模型應用，分別是返青期結束(第89天)、拔節(jié)期結束(第121天)和乳熟期結束(第145天)時。預報模式下市級與省級模型的應用方法同監(jiān)測模式所述過程。

本文中估產(chǎn)精度評價的指標包括絕對誤差(Absolute error，AE)和絕對相對精度(Absolute relative accuracy，ARA)，計算式為

YAE=|Ye-Ya|

(3)

(4)

式中Ye——山東省2017年冬小麥模型估計產(chǎn)量

Ya——山東省2017年冬小麥統(tǒng)計數(shù)據(jù)的真實產(chǎn)量

3 結果與分析

3.1 自變量與小麥產(chǎn)量的相關性分析

對2007—2017年縣級小麥產(chǎn)量數(shù)據(jù)分別與其對應縣區(qū)生育期內(nèi)的技術產(chǎn)量、累計平均EVI、累計平均CWSI進行相關性分析，結果如表1所示。

各市小麥產(chǎn)量和技術產(chǎn)量因子之間的相關性最高，相關系數(shù)均不小于0.974，顯著性概率均高于0.01水平；其次與產(chǎn)量相關度較高的是累計平均EVI，相關系數(shù)均在0.522～0.867之間，也在0.01水平顯著；累計平均CWSI與產(chǎn)量的相關系數(shù)在0.370～0.650之間，大部分通過0.01水平的顯著性檢驗。

對全省小麥產(chǎn)量數(shù)據(jù)與全省范圍內(nèi)的生育期內(nèi)的各指數(shù)也進行相關性分析。其中與產(chǎn)量相關性最高的同樣是技術產(chǎn)量，相關系數(shù)高達0.993，在0.01水平顯著；累計平均EVI與產(chǎn)量的相關性達到了0.780，在0.01水平顯著；累計平均CWSI與產(chǎn)量的相關性相對較低，為0.388，在0.05水平顯著。

由表1可以看出，技術產(chǎn)量、累計平均EVI、累計平均CWSI與實際產(chǎn)量都有良好的相關性，均可作為對冬小麥單產(chǎn)估計進行建模的因子。

表1 小麥產(chǎn)量與技術產(chǎn)量、累計平均EVI、累計平均CWSI的相關系數(shù)Tab.1 Correlation coefficients between wheat yield and trend yield, cumulative EVI and cumulative CWSI

3.2 小麥估產(chǎn)模型

以各縣技術產(chǎn)量Yt、累計平均EVI和累計平均CWSI為自變量，以相應縣的實際歷史產(chǎn)量為因變量，基于最小二乘多元線性回歸法建立了2007—2017年山東省13個市和1個全省的產(chǎn)量估算模型，如表2所示。估產(chǎn)模型的R2均不小于0.962并且在0.01水平顯著。其中，泰安市小麥單產(chǎn)估算模型R2最高(0.990)，濟南市小麥單產(chǎn)估算模型R2最低(0.962)，全省小麥產(chǎn)量估算模型R2為0.985。

表2 冬小麥估產(chǎn)模型Tab.2 Winter wheat yield estimation model

3.3 小麥估產(chǎn)模型的應用與評價

對山東省各個城市和全省的產(chǎn)量估算模型分別在監(jiān)測模式和預測模式下進行應用，單產(chǎn)估算結果如表3所示，精度驗證結果如表4所示。

表3 2017年山東省小麥產(chǎn)量監(jiān)測和預報結果Tab.3 Monitoring and forecasting results of wheat yield of Shandong Province in 2017 kg/hm2

表4 模型精度驗證結果Tab.4 Model accuracy verification results

監(jiān)測模式下精度驗證結果表明，所構建的模型在市級和省級作物單產(chǎn)估測中都很可靠。全省范圍內(nèi)的估產(chǎn)精度為96.91%。各市的估測精度均不小于89.64%，其中模型精度最高為菏澤市，高達99.31%，估產(chǎn)精度最低為濟寧市，為89.64%。

預報模式下的精度驗證結果表明，在大多數(shù)地區(qū)，隨著生育期內(nèi)時間的推進，現(xiàn)勢遙感數(shù)據(jù)不斷加入模型，預報模式下估產(chǎn)模型的精度逐步提高，越來越逼近監(jiān)測模式下模型的估產(chǎn)結果。全省3個時間點的預報精度分別達到了96.44%、97.13%、96.91%。各市中，青島市預報精度最高，3個時間點的預報精度均達到了99%以上。

4 討論

針對小麥的生長特點和研究區(qū)域特點，結合歷史產(chǎn)量數(shù)據(jù)和遙感數(shù)據(jù)，對山東省小麥產(chǎn)量進行建模，并進行了監(jiān)測模式和預報模式下的應用和驗證。該模型在農(nóng)業(yè)遙感產(chǎn)量估算中有一定的參考價值。分析本文對產(chǎn)量預測建模的方法和結果，可能會對建模結果產(chǎn)生影響的因素有：

(1)采用累計平均EVI和累計平均水分脅迫指數(shù)反映環(huán)境因素對冬小麥產(chǎn)量的影響，但在實際生產(chǎn)過程中，小麥產(chǎn)量會受到多種環(huán)境因素的影響，以及各種人為因素造成的產(chǎn)量波動。而這部分影響因素在本文的模型中難以用遙感數(shù)據(jù)表達和預測。

(2)使用特定的作物掩膜可以改善作物估產(chǎn)模型的精度[30]，但由于特定作物的空間分布每年均變化，獲取時間序列內(nèi)每年特定作物的掩膜數(shù)據(jù)工作量大，且存在分類誤差，每年的誤差和分布也會不同，作物分布圖的誤差也會傳遞到后續(xù)的模型中。耕地掩膜相對于冬小麥掩膜更容易獲取。此外，有研究指出使用耕地掩膜替代特定作物的掩膜是可行的[1]。綜上考慮，本文對冬小麥像元篩選時，使用耕地掩膜，可能導致作物的特異性被忽略，參與運算的像元中混雜了其他作物或裸地的信息，從而影響估產(chǎn)模型的精度。

(3)本文技術產(chǎn)量是對估產(chǎn)模型精度具有重要影響的變量。在估算技術產(chǎn)量的過程中，所使用的時間序列長度和計算方法是影響技術產(chǎn)量估算精度的兩個重要因素。只使用了近10年的冬小麥歷史產(chǎn)量數(shù)據(jù)進行技術產(chǎn)量的估算。長時間序列的歷史產(chǎn)量有助于提高技術產(chǎn)量的精度，但在縣、村級精細規(guī)模上，長時間序列的歷史產(chǎn)量數(shù)據(jù)獲取困難。

(4)由于云雨等不良天氣的影響，部分地區(qū)的數(shù)據(jù)存在缺失值。在計算縣級平均EVI和CWSI的過程中，僅選取了在同期數(shù)據(jù)中EVI、CWSI均無缺失情況的像元參與平均計算，減少了參與建模的像元數(shù)量。若采取插值方法對缺失像元進行插補，可以得到全覆蓋的數(shù)據(jù)集，同時插值的誤差也會傳遞到估產(chǎn)模型中，通過插值得到全覆蓋的數(shù)據(jù)是否可以提高建模精度還需進一步討論。

(5)僅針對山東省范圍內(nèi)數(shù)據(jù)進行模型的構建與應用，在預報模式中選取了山東省冬小麥生育期的3個時間節(jié)點進行產(chǎn)量預報，并認為在全省范圍內(nèi)冬小麥的生育期節(jié)點均一致。若有更詳細的，如分縣的生育期數(shù)據(jù)，則可按照各縣具體的生育期時間階段進行模型輸入變量的計算。

5 結論

(1)技術產(chǎn)量、EVI和CWSI與冬小麥產(chǎn)量具有顯著相關性，可以作為冬小麥估產(chǎn)建模的因子。結果表明，技術產(chǎn)量與冬小麥實際產(chǎn)量的相關性最高，在估產(chǎn)模型中為主要影響因子，體現(xiàn)了研究區(qū)域的客觀種植條件和農(nóng)作物生產(chǎn)水平；累計平均EVI與冬小麥實際產(chǎn)量的相關性次之，主要表現(xiàn)了作物的實際生長狀況；累計平均CWSI與冬小麥實際產(chǎn)量的相關性相對最低，主要體現(xiàn)了研究區(qū)域的氣候條件、土壤水分和干旱災害等情況。

(2)基于技術產(chǎn)量、EVI和CWSI建立的省級和市級估產(chǎn)模型的R2均不小于0.962。監(jiān)測模式的驗證結果表明，本研究所構建的模型對市級和省級作物單產(chǎn)的估測都有較高的適用性，采用技術產(chǎn)量與遙感數(shù)據(jù)結合的模式可以實現(xiàn)高精度估產(chǎn)；預報模式驗證結果表明，采用現(xiàn)勢數(shù)據(jù)和歷史數(shù)據(jù)結合的方式可實現(xiàn)冬小麥實時產(chǎn)量預報。本研究可以為作物生長狀況監(jiān)測提供參考依據(jù)。