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

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

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

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

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

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)

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

1.2 數(shù)據(jù)

所使用的數(shù)據(jù)包括3個來源。

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

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

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

2 研究方法

2.1 技術(shù)路線

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

圖1 技術(shù)路線圖Fig.1 Technical flowchart

2.2 各類指標(biāo)計算

2.2.1 農(nóng)業(yè)氣象類指標(biāo)計算

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

2.2.2 種植情況變化指標(biāo)計算

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

2.2.3 遙感植被指數(shù)計算

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

2.2.4 趨勢單產(chǎn)的計算

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

2.3 估產(chǎn)分區(qū)

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

2.4 指標(biāo)重要性計算與篩選

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

步驟1對參與建模的指標(biāo)在黑龍江整省內(nèi)針對3種不同的作物進(jìn)行重要性排序,進(jìn)行指標(biāo)的初步篩選。

步驟2利用初步篩選的指標(biāo),基于估產(chǎn)分區(qū)的結(jié)果,在各區(qū)內(nèi)針對3種作物進(jìn)行重要性排序,逐一將變量加入隨機(jī)森林中建立估產(chǎn)模型,評估所建模型的估產(chǎn)精度。

2.5 估產(chǎn)模型構(gòu)建與精度驗證

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

(1)

3 結(jié)果與分析

3.1 分區(qū)結(jié)果

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

3.2 指標(biāo)篩選結(jié)果

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

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

表1 水稻、大豆、玉米的各分區(qū)初篩指標(biāo)的重要性排序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指數(shù)、累積降水、降水距平、生長度日、高溫度日、有效積溫和積溫距平圖4 各估產(chǎn)分區(qū)內(nèi)不同類型初篩指標(biāo)的重要性排序均值Fig.4 The importance ranking mean of different types of initial screening indicators within each yield estimation zone

3.3 估產(chǎn)模型及其精度驗證結(jié)果

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

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

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

4 結(jié)論與討論

4.1 結(jié)論

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

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

(2)在全區(qū)基于隨機(jī)森林的變量重要性評價結(jié)果顯示在溫度類指標(biāo)中積溫距平最優(yōu),在降水類指標(biāo)中降水距平、單站降水Z指數(shù)、累積降水指標(biāo)表現(xiàn)相當(dāng)。分區(qū)進(jìn)一步評價變量重要性的結(jié)果顯示趨勢單產(chǎn)的重要性>NDVI距平的重要性>積溫距平的重要性>Z指數(shù)的重要性。

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

4.2 討論

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

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

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

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