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

摘 要 為尋求高效的水稻產(chǎn)量估算方法，以2017年長春市九臺和德惠地區(qū)的采樣點為樣本，遙感數(shù)據(jù)和氣象數(shù)據(jù)為特征變量，通過對產(chǎn)量與特征變量間的相關(guān)性分析與特征變量之間的主成分分析和袋外數(shù)據(jù)(out-of-data,OOB)變量的重要性分析對特征變量進(jìn)行選擇，以選擇后的特征變量為輸入變量建立水稻產(chǎn)量估算的隨機(jī)森林回歸(RFR)模型。結(jié)果表明：特征變量優(yōu)選后的RFR模型對水稻產(chǎn)量估算的精度更高，決定系數(shù)R2和平均相對誤差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。

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

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

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

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

注：GS(1-5)表示水稻5個生育期，例B3GS(1-5)表示5個不同生育期的紅光波段反射率數(shù)值，其余變量的下角標(biāo)GS(1-5)同理。FPARmax與FPARmin分別為0.950和0.001，表示植被覆蓋度最大和無植被覆蓋時的FPAR值，其取值與植被類型無關(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)作物生長的關(guān)鍵因素并且氣候條件的好壞直接影響著農(nóng)作物產(chǎn)量的高低，這是因為農(nóng)作物的生長完全或基本在自然條件下，作物的生長和發(fā)育階段受氣象因素的影響，適宜的溫度降水和輻射條件才能促進(jìn)作物生長，本研究選取水稻不同生育期的溫度、降水以及輻射數(shù)據(jù)作為影響作物產(chǎn)量的特征變量。表2為氣象變量及其單位，共計28個。從中國氣象數(shù)據(jù)網(wǎng)上獲取氣象站點數(shù)據(jù)，包括最高溫度、最低溫度、平均溫度、平均降水量以及日太陽輻射，將日太陽輻射進(jìn)行不同生育期的累積，獲得每個生育期太陽輻射累積值，將5個生育期的降水量進(jìn)行累積獲取生育期總降水量，根據(jù)日平均溫度計算8月平均溫度，并利用ARCGIS軟件通過反距離權(quán)重的插值方法插值出與遙感數(shù)據(jù)分辨率相同的30 m分辨率的矢量，并從矢量中提取研究區(qū)內(nèi)各采樣點的氣象信息。

1.2.3實測數(shù)據(jù)

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

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

注：GS(1-5)表示水稻5個生育期，例平均降水量GS(1-5)表示每個不同生育期內(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水稻采樣點分布圖
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ān)系數(shù)來體現(xiàn)產(chǎn)量與特征變量之間的線性相關(guān)程度，相關(guān)系數(shù)的公式為：

(1)

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

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

(2)

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

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

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

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

1.3.3模型參數(shù)確定

RFR模型對研究區(qū)水稻產(chǎn)量的估算借助于R語言中的Random Forest程序包，在該模型中主要有2個參數(shù)需要確定：決策樹個數(shù)以及隨機(jī)選擇的變量個數(shù)m。回歸樹個數(shù)將直接影響預(yù)測結(jié)果的誤差，但當(dāng)決策樹的個數(shù)為一個合適的數(shù)值時，袋外誤差的變化將趨于恒定不變，本研究RFR模型中決策樹的個數(shù)均根據(jù)決策樹個數(shù)與誤差的關(guān)系圖確定，如圖3所示。隨機(jī)選擇的變量個數(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)用九臺德惠地區(qū)的采樣點進(jìn)行建模，由于模型輸入變量不同水稻產(chǎn)量估算結(jié)果也不盡相同，針對變量選擇的不同分別建立了不同水稻產(chǎn)量估算的隨機(jī)森林回歸模型。首先用全部73個變量作為模型的輸入變量，建立水稻產(chǎn)量估算的RFR1模型；其次分析全部特征變量與產(chǎn)量之間的相關(guān)性，提取15個相關(guān)性較高的變量(其相關(guān)性>0.6)建立水稻產(chǎn)量估算的RFR2模型；其提取的15個變量與產(chǎn)量的相關(guān)性表3所示；對該15個相關(guān)性高的變量進(jìn)行主成分分析，提取3個主成分分析結(jié)果，累計貢獻(xiàn)率為86.040%，建立水稻產(chǎn)量估算的RFR3模型；在RFR2的基礎(chǔ)上，對15個相關(guān)性高的變量進(jìn)行了重要性排序分析，剔除變量重要性排序低的變量(%IncMSE為負(fù)值)，將剩余變量重新作為輸入變量建立了RFR4模型。特征變量重要性排序圖如圖4所示。與此同時，對全部原始變量進(jìn)行主成分分析，提取了10個主成分，累計貢獻(xiàn)率為96.670%，以這10個主成分作為模型輸入變量，建立水稻產(chǎn)量估算的RFR5模型；對10個主成分與水稻產(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 模型精度評價

1.4.1留一法交叉驗證

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

1.4.2決定系數(shù)

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

1.4.3平均相對誤差

平均相對誤差(MRE)是多個樣本測量值與估算值之間相對誤差的平均值，用來作為評價水稻產(chǎn)量估算的結(jié)果與實測產(chǎn)量間的誤差的一個標(biāo)準(zhǔn)，其計算公式為

(3)

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

2 結(jié)果與分析

2.1 模型精度分析

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

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

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

2.2 模型適用性分析

由上分析可知RFR3模型水稻估算精度更好，為了進(jìn)一步評判RFR3模型的適用性，在原本研究區(qū)樣本點的基礎(chǔ)上加入農(nóng)安地區(qū)的6個樣本點進(jìn)行建模，農(nóng)安地區(qū)在地形、土壤類型、以及氣候類型上與研究區(qū)較為接近，并利用了SPSS軟件對RFR3模型輸入變量與產(chǎn)量間進(jìn)行了多元逐步回歸，逐步回歸結(jié)果輸入變量為第三生育期的EVI指數(shù)，移除了其余變量，關(guān)系式為：y=9 596.123×EVIGS3+980.356,將多元逐步回歸結(jié)果與RFR3模型進(jìn)行對比，表4為精度對比結(jié)果，結(jié)果表明應(yīng)用RFR3模型對農(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ā)展對于農(nóng)作物產(chǎn)量估算方面的需求。

3 結(jié) 論

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

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

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

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

4 討 論

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