基于隨機(jī)森林回歸方法的水稻產(chǎn)量遙感估算

楊北萍 陳圣波 于海洋 安 秦 (吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130012)

摘 要 為尋求高效的水稻產(chǎn)量估算方法，以2017年長(zhǎng)春市九臺(tái)和德惠地區(qū)的采樣點(diǎn)為樣本，遙感數(shù)據(jù)和氣象數(shù)據(jù)為特征變量，通過(guò)對(duì)產(chǎn)量與特征變量間的相關(guān)性分析與特征變量之間的主成分分析和袋外數(shù)據(jù)(out-of-data,OOB)變量的重要性分析對(duì)特征變量進(jìn)行選擇，以選擇后的特征變量為輸入變量建立水稻產(chǎn)量估算的隨機(jī)森林回歸(RFR)模型。結(jié)果表明：特征變量?jī)?yōu)選后的RFR模型對(duì)水稻產(chǎn)量估算的精度更高，決定系數(shù)R2和平均相對(duì)誤差MRE分別為0.950和0.060；并將該模型應(yīng)用到農(nóng)安地區(qū)，以多元逐步回歸模型作為比較模型，表明RFR模型的水稻產(chǎn)量估算精度明顯優(yōu)于多元逐步回歸模型，RFR模型的R2和MRE分別為0.730和0.090，多元逐步回歸模型的R2和MRE分別為0.530和0.120。

水稻是我國(guó)的主要糧食作物之一，水稻產(chǎn)量的多少直接關(guān)系著農(nóng)民的生活保障以及農(nóng)業(yè)經(jīng)濟(jì)的發(fā)展進(jìn)程，隨著農(nóng)業(yè)經(jīng)濟(jì)的發(fā)展，水稻產(chǎn)量估算費(fèi)時(shí)費(fèi)力的方法已經(jīng)難以滿足農(nóng)業(yè)經(jīng)濟(jì)發(fā)展的需求，所以找尋水稻產(chǎn)量及時(shí)、準(zhǔn)確的估算方法也越來(lái)越重要。遙感作為一種新興的探測(cè)技術(shù)，能夠準(zhǔn)確快速的收集信息，目前，越來(lái)越多的研究學(xué)者利用遙感手段獲取信息來(lái)進(jìn)行農(nóng)作物產(chǎn)量估算方面的研究。Hamar等[1]基于Landsat影像獲取的不同遙感植被指數(shù)與玉米和小麥的產(chǎn)量建立了線性回歸模型；洪雪[2]基于高光譜遙感數(shù)據(jù)的植被指數(shù)建立了水稻的產(chǎn)量關(guān)系模型，湯斌等[3]基于遙感和氣象數(shù)據(jù)對(duì)江蘇省水稻面積監(jiān)測(cè)和估產(chǎn)進(jìn)行了研究，但限于不同模型的選擇雖采用了遙感技術(shù)但在產(chǎn)量估算精度上仍存在較大差異。

隨機(jī)森林回歸(Random forest regression，RFR)模型是一種基于機(jī)器學(xué)習(xí)的統(tǒng)計(jì)方法，它能處理很高維度的數(shù)據(jù)并且估算結(jié)果具有較高的準(zhǔn)確率。據(jù)此已有國(guó)內(nèi)外研究學(xué)者將該方法應(yīng)用于農(nóng)業(yè)遙感估算上來(lái)。王麗愛(ài)等[4]基于RFR模型結(jié)合多種植被指數(shù)對(duì)小麥葉片SPAD值進(jìn)行了遙感估算，結(jié)果表明RFR模型對(duì)比支持向量回歸模型和BP神經(jīng)網(wǎng)絡(luò)模型表現(xiàn)了最強(qiáng)的學(xué)習(xí)能力；張偉等[5]基于遙感影像數(shù)據(jù)和樣地調(diào)查數(shù)據(jù)利用RFR模型對(duì)江山市公益林生物量進(jìn)行了估算，估算精度較高。Jeong等[6]應(yīng)用RFR模型通過(guò)多種變量的輸入對(duì)全球農(nóng)作物的產(chǎn)量進(jìn)行了估算，結(jié)果顯示RFR模型是一種高效、可靠地估算方法。然而在目前的研究中應(yīng)用RFR模型對(duì)農(nóng)作物單點(diǎn)產(chǎn)量進(jìn)行估算的研究還相對(duì)較少，并且大部分學(xué)者應(yīng)用RFR模型進(jìn)行估算時(shí)，均采用特征變量直接輸入的方式，并未對(duì)特征變量進(jìn)行分析。

本研究以遙感數(shù)據(jù)以及氣象數(shù)據(jù)作為特征變量，通過(guò)對(duì)產(chǎn)量與特征變量間的相關(guān)性分析、特征變量之間的主成分分析和OOB變量重要性分析，選擇最優(yōu)數(shù)據(jù)作為輸入的特征變量建立水稻產(chǎn)量估算的RFR模型，以期提高利用該模型進(jìn)行水稻產(chǎn)量估算的精度。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)為長(zhǎng)春市九臺(tái)、德惠以及農(nóng)安地區(qū)，是吉林省水稻主產(chǎn)區(qū)，位于125°14′～126°30′E，43°50′～44°55′N。研究區(qū)氣候以溫帶大陸性季風(fēng)氣候?yàn)橹鳎?017年最高氣溫35 ℃，最低氣溫-29 ℃，年平均氣溫6 ℃；日平均降水量6.7 mm，累積降水量為1 423.9 mm，降水主要集中在6—9月；累積太陽(yáng)總輻射達(dá)到7 923.44 MJ/m2。

1.2 數(shù)據(jù)獲取

1.2.1遙感數(shù)據(jù)

水稻從生長(zhǎng)到成熟主要經(jīng)歷出苗、分蘗、抽穗、灌漿和成熟5個(gè)生育期，由于HJ-1A/B和Landsat8衛(wèi)星遙感數(shù)據(jù)空間分辨率相同，本研究從HJ-1A/B和Landsat8衛(wèi)星數(shù)據(jù)中獲取2017年水稻全生育期(growing-season，GS)的遙感圖像，使用ENVI軟件對(duì)HJ-1A/B衛(wèi)星數(shù)據(jù)以及Landsat衛(wèi)星數(shù)據(jù)進(jìn)行預(yù)處理，主要包括正射校正、輻射定標(biāo)、大氣校正。將原始多光譜圖像的DN值轉(zhuǎn)換為輻射亮度值。表1 為2種遙感數(shù)據(jù)來(lái)源以及影像獲取時(shí)間。本研究選擇HJ-1A/B衛(wèi)星和Landsat8衛(wèi)星光譜覆蓋范圍內(nèi)的藍(lán)光波段(B1)、綠光波段(B2)、紅光波段(B3)、近紅外波段(B4)反射率數(shù)據(jù)以及由這4個(gè)波段構(gòu)建的4種植被指數(shù)和影響作物生長(zhǎng)的光合有效輻射的吸收比例(FPAR)指數(shù)作為遙感變量[7-9]，共計(jì)45個(gè)。

表1 遙感變量及其說(shuō)明Table 1 Remote sensing variable and description

注：GS(1-5)表示水稻5個(gè)生育期，例B3GS(1-5)表示5個(gè)不同生育期的紅光波段反射率數(shù)值，其余變量的下角標(biāo)GS(1-5)同理。FPARmax與FPARmin分別為0.950和0.001，表示植被覆蓋度最大和無(wú)植被覆蓋時(shí)的FPAR值，其取值與植被類型無(wú)關(guān)[10]。

Note：GS(1-5)represents 5 growth periods of rice. Example B3GS(1-5)represents the reflectance value of red light at 5 different growth periods, and the lower corner of other variablesGS(1-5)is in the same way. FPARmaxand FPARminare respectively 0.950 and 0.001, which indicates the FPAR value when vegetation coverage is maximum and without vegetation coverage, and the value is independent of vegetation type[10].

1.2.2氣象數(shù)據(jù)

氣候條件是影響農(nóng)作物生長(zhǎng)的關(guān)鍵因素并且氣候條件的好壞直接影響著農(nóng)作物產(chǎn)量的高低，這是因?yàn)檗r(nóng)作物的生長(zhǎng)完全或基本在自然條件下，作物的生長(zhǎng)和發(fā)育階段受氣象因素的影響，適宜的溫度降水和輻射條件才能促進(jìn)作物生長(zhǎng)，本研究選取水稻不同生育期的溫度、降水以及輻射數(shù)據(jù)作為影響作物產(chǎn)量的特征變量。表2為氣象變量及其單位，共計(jì)28個(gè)。從中國(guó)氣象數(shù)據(jù)網(wǎng)上獲取氣象站點(diǎn)數(shù)據(jù)，包括最高溫度、最低溫度、平均溫度、平均降水量以及日太陽(yáng)輻射，將日太陽(yáng)輻射進(jìn)行不同生育期的累積，獲得每個(gè)生育期太陽(yáng)輻射累積值，將5個(gè)生育期的降水量進(jìn)行累積獲取生育期總降水量，根據(jù)日平均溫度計(jì)算8月平均溫度，并利用ARCGIS軟件通過(guò)反距離權(quán)重的插值方法插值出與遙感數(shù)據(jù)分辨率相同的30 m分辨率的矢量，并從矢量中提取研究區(qū)內(nèi)各采樣點(diǎn)的氣象信息。

1.2.3實(shí)測(cè)數(shù)據(jù)

2017年在研究區(qū)內(nèi)選擇了多個(gè)水稻采樣點(diǎn)，記錄水稻的有效株數(shù)和采集水稻真實(shí)樣品，并進(jìn)行脫粒、烘干和稱重等操作。根據(jù)水稻的有效株數(shù)和重量來(lái)計(jì)算水稻的實(shí)際產(chǎn)量。本研究九臺(tái)德惠地區(qū)共收集16個(gè)水稻采樣點(diǎn)，農(nóng)安地區(qū)6個(gè)水稻采樣點(diǎn)，圖1 為研究區(qū)2017年水稻采樣點(diǎn)分布圖。

表2 氣象變量及其單位Table 3 Meteorological variables and their units

注：GS(1-5)表示水稻5個(gè)生育期，例平均降水量GS(1-5)表示每個(gè)不同生育期內(nèi)的平均降水量，其余變量的下角標(biāo)GS(1-5)同理。

Note：GS(1-5)refers to the five growth periods of rice, 平均降水量GS(1-5)refers to the average precipitation in each different growth period, and the lower cornerGS(1-5)of the other variables is the same.

圖1 研究區(qū)2017水稻采樣點(diǎn)分布圖
Fig.1 Distribution map of 2017 rice sampling sites in the research area

1.3 研究方法

隨機(jī)森林回歸模型原理流程圖如圖2所示。

圖2 RFR模型原理流程圖
Fig.2 Flowchart of model principle

1.3.1變量分析

相關(guān)性分析，通過(guò)相關(guān)系數(shù)來(lái)體現(xiàn)產(chǎn)量與特征變量之間的線性相關(guān)程度，相關(guān)系數(shù)的公式為：

(1)

主成分分析，是將一組相關(guān)變量通過(guò)線性變換轉(zhuǎn)成另一組不相關(guān)的變量，提取的主成分變量最大的包含原變量的所有信息，達(dá)到降維的目的并使得變量間相互獨(dú)立。

袋外數(shù)據(jù)(out-of-bag data，OOB)重要性分析主要基于OOB數(shù)據(jù)，袋外數(shù)據(jù)是模型進(jìn)行中對(duì)訓(xùn)練集做有放回隨機(jī)抽樣時(shí)每次未被抽到的樣本點(diǎn)組成的數(shù)據(jù)集，通過(guò)袋外誤差增長(zhǎng)百分率來(lái)衡量特征變量的重要性，針對(duì)一個(gè)決策樹，將OOB數(shù)據(jù)對(duì)應(yīng)變量打亂前打亂后分別帶入決策樹[11]，計(jì)算其誤差的增長(zhǎng)百分率(IncMSE%)，假設(shè)森林中有N棵樹，對(duì)于第K顆樹的誤差增長(zhǎng)百分率為：

(2)

其中i為某一變量，OOBK1對(duì)應(yīng)的袋外誤差，OOBK2對(duì)應(yīng)的打亂后袋外誤差。

對(duì)于N棵樹如果該變量在OOB數(shù)據(jù)上打亂后對(duì)決策樹的結(jié)果沒(méi)什么影響，及打亂后的均方誤差的差值很小，則說(shuō)明該變量不重要[12-13]。

1.3.2選擇最佳分割節(jié)點(diǎn)

在對(duì)決策樹進(jìn)行分割時(shí)，設(shè)每個(gè)觀測(cè)值對(duì)應(yīng)n個(gè)特征，則在每一棵樹的每個(gè)節(jié)點(diǎn)處隨機(jī)從n個(gè)特征中無(wú)放回的隨機(jī)抽取m個(gè)特征(m≤n)，選擇一個(gè)最佳分割屬性作為節(jié)點(diǎn)創(chuàng)建決策樹，對(duì)于回歸模型，最佳分割屬性的評(píng)判標(biāo)準(zhǔn)為使分割后兩部分樣本的均方差結(jié)果達(dá)到最小，然后在分叉的2個(gè)節(jié)點(diǎn)處再利用這樣的準(zhǔn)則，選擇之后的分割屬性，且分割過(guò)程不需要剪枝[14-17]，直到達(dá)到葉子節(jié)點(diǎn)為止。

1.3.3模型參數(shù)確定

RFR模型對(duì)研究區(qū)水稻產(chǎn)量的估算借助于R語(yǔ)言中的Random Forest程序包，在該模型中主要有2個(gè)參數(shù)需要確定：決策樹個(gè)數(shù)以及隨機(jī)選擇的變量個(gè)數(shù)m?；貧w樹個(gè)數(shù)將直接影響預(yù)測(cè)結(jié)果的誤差，但當(dāng)決策樹的個(gè)數(shù)為一個(gè)合適的數(shù)值時(shí)，袋外誤差的變化將趨于恒定不變，本研究RFR模型中決策樹的個(gè)數(shù)均根據(jù)決策樹個(gè)數(shù)與誤差的關(guān)系圖確定，如圖3所示。隨機(jī)選擇的變量個(gè)數(shù)程序包一般默認(rèn)為總變量的1/3[18]。

圖3 決策樹與袋外誤差關(guān)系圖
Fig.3 Relation between decision tree and outside bag error

1.3.4水稻產(chǎn)量估算的RFR模型建立

首先應(yīng)用九臺(tái)德惠地區(qū)的采樣點(diǎn)進(jìn)行建模，由于模型輸入變量不同水稻產(chǎn)量估算結(jié)果也不盡相同，針對(duì)變量選擇的不同分別建立了不同水稻產(chǎn)量估算的隨機(jī)森林回歸模型。首先用全部73個(gè)變量作為模型的輸入變量，建立水稻產(chǎn)量估算的RFR1模型；其次分析全部特征變量與產(chǎn)量之間的相關(guān)性，提取15個(gè)相關(guān)性較高的變量(其相關(guān)性>0.6)建立水稻產(chǎn)量估算的RFR2模型；其提取的15個(gè)變量與產(chǎn)量的相關(guān)性表3所示；對(duì)該15個(gè)相關(guān)性高的變量進(jìn)行主成分分析，提取3個(gè)主成分分析結(jié)果，累計(jì)貢獻(xiàn)率為86.040%，建立水稻產(chǎn)量估算的RFR3模型；在RFR2的基礎(chǔ)上，對(duì)15個(gè)相關(guān)性高的變量進(jìn)行了重要性排序分析，剔除變量重要性排序低的變量(%IncMSE為負(fù)值)，將剩余變量重新作為輸入變量建立了RFR4模型。特征變量重要性排序圖如圖4所示。與此同時(shí)，對(duì)全部原始變量進(jìn)行主成分分析，提取了10個(gè)主成分，累計(jì)貢獻(xiàn)率為96.670%，以這10個(gè)主成分作為模型輸入變量，建立水稻產(chǎn)量估算的RFR5模型；對(duì)10個(gè)主成分與水稻產(chǎn)量間進(jìn)行相關(guān)性分析，發(fā)現(xiàn)只有第二主成分與產(chǎn)量的相關(guān)性較大(相關(guān)系數(shù)為0.638)，以第二主成分為輸入變量建立RFR6模型。

表3 水稻產(chǎn)量與特征變量相關(guān)性表Table 3 Correlation Table of rice yield and characteristic variables

圖4 特征變量重要性排序圖
Fig.4 Ranking diagram of importance of feature variables

1.4 模型精度評(píng)價(jià)

1.4.1留一法交叉驗(yàn)證

對(duì)于輸入變量不同而建立的模型分別應(yīng)用留一法進(jìn)行交叉驗(yàn)證，每次選出一個(gè)樣本進(jìn)行驗(yàn)證，其他樣本全部作為訓(xùn)練樣本，然后建模并驗(yàn)證一個(gè)測(cè)試樣本的估算精度以及誤差，直至所有樣本均參與了驗(yàn)證[13]，該實(shí)驗(yàn)有16個(gè)水稻樣本點(diǎn)，每個(gè)模型均將進(jìn)行16次的交叉驗(yàn)證。避免了因樣本選擇出現(xiàn)的偶然性，可以有效的對(duì)模型的穩(wěn)定性進(jìn)行評(píng)價(jià)。

1.4.2決定系數(shù)

決定系數(shù)(coefficient of determination,R2)也稱擬合優(yōu)度，是相關(guān)系數(shù)的平方，它的大小決定了相關(guān)的密切程度，R2越接近1，表示2個(gè)數(shù)據(jù)擬合優(yōu)度越好，相反，越接近0，表示擬合結(jié)果越差。

1.4.3平均相對(duì)誤差

平均相對(duì)誤差(MRE)是多個(gè)樣本測(cè)量值與估算值之間相對(duì)誤差的平均值，用來(lái)作為評(píng)價(jià)水稻產(chǎn)量估算的結(jié)果與實(shí)測(cè)產(chǎn)量間的誤差的一個(gè)標(biāo)準(zhǔn)，其計(jì)算公式為

(3)

式中：xm表示水稻產(chǎn)量測(cè)量值，xe表示產(chǎn)量估算值，N表示樣本個(gè)數(shù)，本研究中對(duì)水稻產(chǎn)量的估算，MRE越小表示估算結(jié)果精度越高。

2 結(jié)果與分析

2.1 模型精度分析

利用留一法交叉驗(yàn)證的方法，每次選出一個(gè)樣本進(jìn)行驗(yàn)證，各模型水稻產(chǎn)量估算的訓(xùn)練集和驗(yàn)證集的平均相對(duì)誤差變化如圖5所示，其中，驗(yàn)證集的平均相對(duì)誤差呈現(xiàn)逐漸減小的趨勢(shì)，訓(xùn)練集的平均相對(duì)誤差趨于恒定，研究區(qū)水稻樣本點(diǎn)各模型的水稻產(chǎn)量估算值與實(shí)測(cè)值的對(duì)比如圖6所示，其中RFR3模型產(chǎn)量估算結(jié)果滿足關(guān)系式：y=0.730 3x+1 868.4，R2為0.949，相比于其他模型最高，平均相對(duì)誤差為0.064，對(duì)比于其他模型最??；而剩余模型中RFR1模型對(duì)比于RFR2模型精度較高；RFR2模型對(duì)比于RFR5模型精度較高。RFR4模型以及RFR6模型對(duì)于水稻產(chǎn)量的估算結(jié)果精度較低。

圖5 水稻產(chǎn)量估算誤差變化
Fig.5 Variation of rice yield estimation error

圖6 水稻產(chǎn)量估算與實(shí)測(cè)對(duì)比
Fig.6 Comparison between rice yield estimation and observed results

2.2 模型適用性分析

由上分析可知RFR3模型水稻估算精度更好，為了進(jìn)一步評(píng)判RFR3模型的適用性，在原本研究區(qū)樣本點(diǎn)的基礎(chǔ)上加入農(nóng)安地區(qū)的6個(gè)樣本點(diǎn)進(jìn)行建模，農(nóng)安地區(qū)在地形、土壤類型、以及氣候類型上與研究區(qū)較為接近，并利用了SPSS軟件對(duì)RFR3模型輸入變量與產(chǎn)量間進(jìn)行了多元逐步回歸，逐步回歸結(jié)果輸入變量為第三生育期的EVI指數(shù)，移除了其余變量，關(guān)系式為：y=9 596.123×EVIGS3+980.356,將多元逐步回歸結(jié)果與RFR3模型進(jìn)行對(duì)比，表4為精度對(duì)比結(jié)果，結(jié)果表明應(yīng)用RFR3模型對(duì)農(nóng)安地區(qū)的水稻產(chǎn)量進(jìn)行估算的結(jié)果依然較好，R2達(dá)到0.730，MRE達(dá)到0.090，明顯優(yōu)于多元逐步回歸模型的估算精度。所以應(yīng)用RFR模型估算水稻產(chǎn)量結(jié)果較為可靠、并且精度較高，可以很好地滿足農(nóng)業(yè)發(fā)展對(duì)于農(nóng)作物產(chǎn)量估算方面的需求。

3 結(jié) 論

本研究以遙感數(shù)據(jù)以及氣象數(shù)據(jù)為特征變量，通過(guò)對(duì)產(chǎn)量與特征變量間的相關(guān)性分析、特征變量之間的主成分分析和OOB變量重要性分析選擇了最優(yōu)的特征變量建立水稻產(chǎn)量估算的RFR模型，同時(shí)建立多元逐步回歸模型與優(yōu)選后的RFR模型的估算結(jié)果進(jìn)行比較，進(jìn)一步評(píng)價(jià)RFR模型的估算精度。研究得到以下結(jié)論：

1)應(yīng)用RFR模型對(duì)研究區(qū)水稻產(chǎn)量估算時(shí)需要對(duì)特征變量進(jìn)行選擇，經(jīng)過(guò)優(yōu)選后的RFR模型比未優(yōu)選的估算結(jié)果精度更高。特征變量的選擇明顯改善了模型估算精度。

2)將優(yōu)選后的RFR模型應(yīng)用到農(nóng)安地區(qū)，農(nóng)安地區(qū)的產(chǎn)量估算結(jié)果較好，初步驗(yàn)證了該模型在產(chǎn)量估算上的適用性。

3)優(yōu)選后RFR模型對(duì)水稻產(chǎn)量的估算精度高于多元逐步回歸的估算結(jié)果。說(shuō)明優(yōu)選后RFR模型能很好的估算農(nóng)作物產(chǎn)量，為農(nóng)作物產(chǎn)量估算方法提供新的參考。

4 討 論

通過(guò)單一變量對(duì)產(chǎn)量進(jìn)行估算時(shí)往往誤差較大，通過(guò)圖6可以發(fā)現(xiàn)結(jié)合多種數(shù)據(jù)進(jìn)行分析，進(jìn)而對(duì)農(nóng)作物產(chǎn)量進(jìn)行估算可以達(dá)到很好的效果，本研究的研究區(qū)為九臺(tái)和德惠地區(qū)，僅對(duì)農(nóng)安地區(qū)的采樣點(diǎn)進(jìn)行了初步驗(yàn)證，該模型的適用性還有待于進(jìn)一步驗(yàn)證。本研究的輸入變量為遙感方法獲取反映植被生長(zhǎng)狀態(tài)的一些指數(shù)數(shù)據(jù)和氣象網(wǎng)站獲取的氣象數(shù)據(jù)。然而，影響作物生長(zhǎng)的因素眾多。本研究只考慮了部分遙感數(shù)據(jù)和氣象數(shù)據(jù)，在一定程度上降低了最終的估算精度。