黑龍江大宗作物估產模型: 以水稻、大豆、玉米為例

楊子毅, 朱秀芳, 代佳佳, 姬忠林, 潘耀忠

(1.北京師范大學地理科學學部遙感科學與工程研究院, 北京 100875; 2.青海師范大學地理科學學院, 西寧 810008; 3.聊城大學地理與環境學院, 聊城 252000)

隨著全球氣候變化和人類活動的不斷加劇,各種自然災害頻繁發生,其中旱澇災害嚴重威脅著農業生產、生態環境[1-2],并且擁有發生頻率高、覆蓋面積廣的特點。根據《中國水旱災害防御公報2020》[3],2020年全國因旱澇災害造成農作物受災面積1 544萬hm2,其中絕收面積201萬hm2,僅干旱造成的糧食損失達1 230.4萬t、經濟作物損失169.81億元。黑龍江作為糧食單產多年保持第一的大省,對其大宗作物單產進行預估有著重要的意義。

進入21世紀以來,遙感技術在農業領域方面的應用發展迅速,其時空分辨率的提升更加有利于農作物的監測和估產[4-6]。目前已有許多學者針對不同作物研究了相應的估產模型,主要的估產模型包括有4種:經驗統計模型、光能利用率模型、作物生長模擬模型和耦合模型。其中經驗統計模型因操作簡單、計算方便,而被廣泛應用?；跉庀髷祿墓喇a模型屬于經驗統計模型的一種,它通過建立各類氣象數據與實際單產之間的關系來對作物單產進行估計[7-8]。遙感數據在經驗統計模型的應用包括兩個方面:一是以遙感波段作為自變量進行估產建模。二是利用遙感數據計算植被指數,以植被指數作為自變量進行建模[9-10]。常用來進行作物估產建模的植被指數有歸一化植被指數(NDVI)、增強型植被指數(enhanced vegetation index,EVI)、綠色植被指數(green vegetation index,GVI)、土壤調節植被指數(soil-adjusted vegetation index,SAVI)等。近年來,部分學者在經驗統計模型中將氣象數據和遙感數據相結合進行估產建模,例如,朱秀芳等[11]利用綜合氣象災害指標、遙感植被指數作為自變量,小麥單產作為因變量,利用隨機森林回歸方法建立了中國五大麥區的小麥單產非線性估算模型。

根據《黑龍江統計年鑒—2021》[12],水稻、大豆、玉米3種作物的總產量為黑龍江農作物總產量的98%,但針對該省3種作物的估產模型研究較少,建模時使用的樣本為單一的遙感數據、氣象數據或者社會統計、生產力數據[13-14],沒有考慮黑龍江省內氣候差異對糧食單產的影響。因此,現首先根據黑龍江省的農業氣候和作物種植情況進行估產分區,其次在各區內以遙感植被指數、氣象指標、趨勢單產作為因變量,以黑龍江的3種大宗作物的實際單產作為自變量,利用隨機森林重要性評價對指標進行篩選,最終選擇精度和重要性排序最高的指標作為模型的輸入變量進行建模,以期為黑龍江大宗作物的估產模型研究提供參考。

1 研究區與數據

1.1 研究區

黑龍江省介于東經121°11′～135°05′,北緯43°26′～53°33′,東北部與俄羅斯接壤,西部與內蒙古相鄰,南部與吉林省相鄰,地勢大致呈西北、北部和東南部高,東北、西南部低,由山地、臺地、平原和水面構成。其屬于溫帶大陸性季風氣候,多年平均降水量為400～800 mm,多年平均氣溫為-4～6 ℃,擁有豐富的黑土資源、較高的農業機械化水平,是全國耕地面積最多的省份,也是全國糧食單產最高的省份[15]。與此同時,黑龍江位于高緯度地區,也是對氣候變暖響應最敏感的地區之一,主要表現為氣候變干、降水量減少等[16],使得農作物種植面臨著極大挑戰。

1.2 數據

所使用的數據包括3個來源。

(1)來自中國氣象局氣象信息中心的中國國家級地面氣象站的2000—2021年的中國氣象要素站點觀測逐日數據,該數據包括日照、降水、氣溫等不同類型的日觀測數據。

(2)來自國家統計局黑龍江省調查總隊的2016—2021年玉米、水稻、大豆作物播種面積、空間分布數據以及2014—2021年的3種作物的歷史單產數據。

(3)來自美國航空航天局的2004—2021年的8 d合成的 MOD09Q1.06 數據,空間分辨率為 250 m。

2 研究方法

2.1 技術路線

技術路線如圖1所示,利用黑龍江氣象站點數據計算各類氣象指標,并且利用目標作物種植面積和空間分布數據,計算種植情況變化指標,基于SKATER(Spatial “K”luster Analysis by Tree Edge Removal)方法對黑龍江氣候、耕地變化情況進行分區。基于MOD09Q1.061遙感數據和土地利用數據計算黑龍江各縣三種作物范圍內的NDVI平均值和距平值。利用各縣歷史單產數據,采用滑動平均法對趨勢單產進行擬合。利用隨機森林重要性評價方法針對3種作物在黑龍江全境內對農業氣象指標進行初步的篩選。在3種作物的不同氣候區域內,對初步篩選的指標進行重要性計算及排序,采用留一法,并且按變量重要性排序逐個將變量加入模型中建模,評估建立的所有模型的精度,選出使得各年份整體建模精度最高的變量作為最終輸入變量進行最終估產模型的建立。

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

2.2 各類指標計算

2.2.1 農業氣象類指標計算

利用2000—2021年黑龍江氣象站點逐日降水、溫度和日照時數數據計算了14個指標,具體包括:降水Z指數、累積降水、降水距平、降水距平平均值、降水距平變化率、生長度日、高溫度日、有效積溫、活動積溫、積溫距平、生長季開始、生長季結束、生長季長度、平均日照時數[17-20]。其中所有指標參與了估產分區,而降水Z指數(PZ)、累積降水(RA)、降水距平(AP)、生長度日(GD)、高溫度日(KD)、有效積溫(EA)、積溫距平(AT)7個指標作為候選的估產建模指標。要特別進行說明的是:在進行估產分區時,各個站點的各類指標以年為尺度進行計算,并且取多年平均值作為該指標的最終值,而在進行估產建模時,各個站點的各類指標以月為尺度進行計算。

2.2.2 種植情況變化指標計算

基于2016—2021年的農作物空間分布數據,分別統計3種作物的每個區縣內作物連續種植1、2、3、4、5年的像元數,進而計算近5年間只種植1次的像元占比;利用2016—2021年的農作物種植面積,計算各區縣的3種作物的種植面積變異系數(coefficient of variation,CV)值。

2.2.3 遙感植被指數計算

NDVI作為遙感重要的植被指數,可以較好地反映植被的生長狀況以及分布特征[21-23]。以8 d合成的250 m分辨率地表反射率產品MOD09Q1.061為數據源,計算2004—2021年共計17年NDVI數據,采用自適應濾波方法[21]對NDVI進行濾波,對濾波后的NDVI按照自然月逐像元計算17年每個月份的NDVI月合成值,從中剔除最大值和最小值,計算剩余15個NDVI合成值的平均值,進而獲得逐像元的多年的NDVI月平均值,進而逐像元計算6年每個月份的NDVI距平(某年某月的NDVI值與該月歷史的NDVI平均值的差值),最后以縣為單元,分別統計3種作物種植范圍內每個月份的NDVI均值(ND)及NDVI距平均值(NA)。

2.2.4 趨勢單產的計算

實際單產一般由趨勢單產(Yt)、氣象單產、隨機單產三部分組成。趨勢單產一般是指由作物品種、社會生產力等因素影響的單產分量,氣象單產是指由氣象因子影響的波動性單產分量,隨機單產是指由一些隨機因素影響的單產分量,但由于隨機單產的不確定,在進行單產模擬的過程中一般忽略不計。采用滑動平均法[24]擬合趨勢單產。該方法將前兩年實際單產的均值作為本年的趨勢單產,例如計算2014年和2015年的實際單產的均值作為2016年的趨勢單產。

2.3 估產分區

黑龍江的熱量資源和降水變率在空間上均呈現不均勻分布[25-26],且不同年份的作物的種植面積和種植范圍分布變化較大。因此,以縣為基本單位,以縣域內年尺度計算的各類氣候指標、目標作物的種植面積變異系數以及近5年間只種植一次的耕地像元占比為輸入變量,使用SKATER[27]聚類方法對黑龍江各縣進行估產分區。SKATER算法是一種基于最小生成樹的圖聚類算法,通過切除最小生成樹的邊來獲得滿足空間鄰接約束的聚類結果。

2.4 指標重要性計算與篩選

隨機森林算法是由Breiman等[28]在2001年提出來的機器學習算法,不僅具有分類回歸的功能,還可以對數據進行重要性評價[29]。使用袋外誤差法進行重要性評價。袋外誤差指的是未參與建模的數據對構建好的分類器的常規誤差進行無偏估計的結果,由于其高效性,不再需要交叉驗證或者使用單獨的測試集進行無偏估計[30]。具體步驟如下。

步驟1對參與建模的指標在黑龍江整省內針對3種不同的作物進行重要性排序,進行指標的初步篩選。

步驟2利用初步篩選的指標,基于估產分區的結果,在各區內針對3種作物進行重要性排序,逐一將變量加入隨機森林中建立估產模型,評估所建模型的估產精度。

2.5 估產模型構建與精度驗證

應用留一法進行最終模型的建立和精度驗證。具體步驟為留出某一年的數據不參與建模而作為驗證數據,應用其他年份數據建立估產模型,利用驗證數據進行估產精度驗證。建模過程按照步驟2執行。例如留出2021年的數據不參與建模,利用2016—2020年數據建模,按照重要性排序,逐一將變量加入隨機森林中建立估產模型,并評估所建模型的估產精度。對按照留一法計算出來的逐年的所有模型的驗證精度情況進行綜合分析,選出所有年份平均精度最高的變量作為最終的建模輸入變量,基于最終篩選的建模輸入變量建立的模型即為最終的估產模型,最后利用留一法中基于最終篩選的建模輸入變量建立的模型的平均絕對相對誤差(mean absolute relative precision,MARE)的統計值(最大值、最小值、平均值)來衡量最終估產模型的精度和穩定性。

(1)

3 結果與分析

3.1 分區結果

基于SKATER方法得到估產氣候分區如圖2所示。整體來說,3種作物的分區結果均不同,且玉米比其他兩種作物多出一個分區。另外,針對大豆、玉米,大興安嶺及周邊地區都被單獨分區,這是因為大興安嶺地區年平均氣溫較其他地區寒冷所致。以下將大豆、水稻、玉米的分區簡稱大豆1、2、3區,水稻1、2、3區,玉米1、2、3、4區。

3.2 指標篩選結果

逐月計算2.2.1節中選擇的建模指標,利用隨機森林方法,對生長季每個月的指標進行重要性評價,進一步按類型計算該類型內各月指標的綜合排名,結果如圖3所示。在氣溫類指標中無論哪種作物的積溫距平(編號為0的指標)的重要性排名均最高;3種降水類指標(編號為4～6的指標)在3種作物中的重要性排名差異很小,但Z指數能更好地反映旱澇等級;NDVI距平的綜合排名高于NDVI。因此,在后續分區建模的溫度類、降水類和植被指數類指標中分別選擇積溫距平、Z指數和NDVI距平參與建模。利用初步篩選后的指標,分區分作物類型進行進一步的重要性評價,結果如表1所示,并計算各類指標的綜合平均排名,結果如圖4所示。大多數情況下,趨勢單產的重要性都排在首位。綜合來看,趨勢單產的重要性>NDVI距平>積溫距平>Z指數。

NA、ND、PZ、RA 、AP、GD、KD、EA和AT分別表示NDVI均值和NDVI距平均值、降水Z指數、累積降水、降水距平、生長度日、高溫度日、有效積溫和積溫距平圖3 全區指標重要性排序Fig.3 The importance ranking of indicators for the entire study area

表1 水稻、大豆、玉米的各分區初篩指標的重要性排序Table 1 Importance ranking of primary screening indicators for rice, soybean, and corn in each zone

NA、ND、PZ、RA 、AP、GD、KD、EA和AT分別表示NDVI均值和NDVI距平均值、降水Z指數、累積降水、降水距平、生長度日、高溫度日、有效積溫和積溫距平圖4 各估產分區內不同類型初篩指標的重要性排序均值Fig.4 The importance ranking mean of different types of initial screening indicators within each yield estimation zone

3.3 估產模型及其精度驗證結果

結合留一法,將變量按重要性排序后逐一加入建模并計算MARE(圖5),其中橫坐標為逐一加入變量后參與建模的變量總個數,縱坐標是6個MARE的范圍(用色帶寬度表示)和平均值(用實線表示)。例如,在水稻1區中第一次加入趨勢單產,通過留一法建立了6個模型(基于2016—2021年,共計6年的數據,每次留出1年的數據不參與建模,而作為驗證數據,用其他年份數據建立模型并進行精度驗證,整個過程經歷6遍,對應6個模型),計算所建立的每個模型的MARE和6個MARE的均值;第二次加入趨勢單產、7月Z指數,重復上面的步驟;第三次加入趨勢單產、5月NDVI距平、5月積溫距平,重復上述步驟,以此類推直到加入最后一個指標建模并計算精度。比較隨變量增加,留一法建立的模型的MARE的均值,選擇MARE的均值最小時對應的變量為最終的變量。例如,在水稻2區,加入最后一個變量后利用留一法計算的各個模型的整體精度最高,因此以全部變量作為最終的建模變量;而針對水稻1區,累計加入11個變量后,利用留一法計算的各個模型的整體精度最高,以前11個變量作為建模的最終變量。在各區針對不同作物篩選的最終變量,及由最終變量建立的模型的精度的范圍如表2所示。整體來說,3種作物中水稻的平均估產精度最高,其次是大豆,最后是玉米,除了玉米1區,其他區域的所有作物的MARE均在0.1以下。

表2 最終各分區的估產模型建模變量和精度Table 2 Modeling variables and accuracy of the final yield estimation model for each zone

圖5 MARE隨變量個數增加的變化Fig.5 Change of MARE with increasing number of variables

4 結論與討論

4.1 結論

綜合使用統計數據、氣象數據和遙感數據,利用SKATER方法在黑龍江省進行氣候分析,利用隨機森林結合留一法,分區、分作物逐步進行建模變量篩選,最終分區建立了水稻、大豆、玉米3種主要作物的估產模型。得到的主要結論如下。

(1)3種作物的分區結果均不相同,其中針對大豆、玉米,大興安嶺地區及周邊地區都被單獨分區。

(2)在全區基于隨機森林的變量重要性評價結果顯示在溫度類指標中積溫距平最優,在降水類指標中降水距平、單站降水Z指數、累積降水指標表現相當。分區進一步評價變量重要性的結果顯示趨勢單產的重要性>NDVI距平的重要性>積溫距平的重要性>Z指數的重要性。

(3)針對3種作物在不同區域篩選的最終指標,使用留一法建立的10個模型的MARE除玉米1區外均小于10%,其中水稻的平均估產精度最高(平均MARE在水稻1區為0.050,在水稻2區為0.071,在水稻3區為0.050),其次是大豆(平均MARE在大豆1區為0.099,在大豆2區為0.057,在大豆3區為0.069),最后是玉米(平均MARE在玉米1區為0.157,在玉米2區為0.059,在玉米3區為0.054,在玉米4區為0.069),所有模型的MARE的方差均小于0.05%,說明估產模型精度表現出良好的穩定性。

4.2 討論

研究的特色主要體現在3個方面。首先,在估產分區時不僅考慮了農業氣候指標,還考慮了種植的波動情況,使用SKATER聚類方法對黑龍江省進行估產分區,有利于精細化估產模型的建立。SKATER不僅考慮了數據的屬性,還考慮了數據的空間鄰接關系,比起傳統的聚類方法(如K-Means、ISODATA),其分區結果在空間上更連續。使用的分區指標和分區方法可以為相關分區研究提供參考。其次,采用了“兩步走”策略篩選建模變量,第一步選擇建模指標的大類,如選用了積溫距平和Z指數作為氣象類指標參與建模,第二步在第一步的基礎上,進一步分區篩選建模變量。這樣做可以在一定程度上避免高度相關的指標、物理意義相近的指標(如降水距平和Z指數)同時被選做建模變量,避免信息冗余。最后,利用隨機森林重要性評價方法和留一法相結合篩選指標,可以保證篩選的指標適用于不同的年景,減少了隨機性引起的波動。

最終選用了積溫距平和Z指數作為氣象類指標參與建模。積溫距平和Z指數可以反映出熱量資源和降水資源的變化,也可以指示農業氣象災害[31-32],如高溫、冷害、洪澇和干旱。然而,影響作物生長的氣象因素不僅僅包括溫度、降水,還有光照、輻射、空氣濕度等并未考慮。此外,除了氣象因素,非氣象因素,例如病蟲害也會影響作物的單產。無論是氣象因素還是非氣象因素,如果影響了作物生長,都可能導致NDVI偏離歷史常年的平均,因此NDVI距平不僅可以反映氣象因子的影響,還可以反映其他脅迫的影響,是一個更加綜合的指標。

在分區后的指標重要性排序中,將玉米、大豆、水稻的不同區域的相同指標的排名分別求平均值,結果顯示趨勢單產對于水稻、大豆、玉米的重要性的綜合平均排名分別為1、1、3名,積溫距平的綜合平均排名分別為7.6、7.6、7.1名,Z指數的綜合平均排名為8.3、9.7、8.2名。不論哪一種作物,排名最高的指標均為趨勢單產,說明趨勢單產對于單產的影響最大,加入趨勢單產建?？梢蕴岣呓５木?。3種作物的氣象類指標中積溫距平的綜合平均排名高于Z指數,這可能是由于氣溫是黑龍江地區農作物生長的主要限制因子。

此外,研究仍然存在一些不足之處:進行趨勢單產模擬時,使用的歷史單產數據集時間為2014—2021年,時間序列有些短,對于最終的建模精度可能存在一定的影響;使用的數據除遙感數據以外,其余數據的空間分辨率較低,例如氣象數據為各觀測站的氣象站點數據,僅可以做到縣級尺度的估產,無法做到更精細尺度的估產建模。