曾 妍，王 迪※，趙小娟

（1. 中國農(nóng)業(yè)科學院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所/農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)遙感重點實驗室，北京100081； 2. 青海省農(nóng)牧業(yè)遙感中心，西寧810008）

0 引言

小麥在中國是僅次于水稻、玉米的第三大糧食作物［1］，小麥估產(chǎn)可為有關部門制定政策和經(jīng)濟計劃提供依據(jù)，在糧食生產(chǎn)、流通和儲備等各環(huán)節(jié)的宏觀調控中發(fā)揮作用。中國是全球最大的小麥生產(chǎn)國和消費國，常年產(chǎn)量約占全球17%［2］，中國小麥產(chǎn)量變化對全球小麥產(chǎn)量的預測有重要影響，進而波及到國際糧價的穩(wěn)定，因此，及時準確估測中國的小麥產(chǎn)量對保證全球糧食安全具有重要意義。

衛(wèi)星遙感能及時高效地獲取大范圍地物信息，廣泛應用于農(nóng)業(yè)監(jiān)測領域，其中產(chǎn)量估算更是農(nóng)情遙感監(jiān)測的重要內(nèi)容。目前，基于遙感信息的作物產(chǎn)量估測方法主要有3種：經(jīng)驗模型或統(tǒng)計模型、機理模型和半經(jīng)驗半機理模型［3-4］。機理模型能夠反映作物產(chǎn)量形成過程，但輸入?yún)?shù)多且校正困難，導致其應用受限。相較于機理模型，半經(jīng)驗半機理模型雖有簡化，但模型輸入中需要的2個關鍵變量，光能利用效率和作物收獲指數(shù)，都難以在時空分布中進行準確地定量模擬［5］。統(tǒng)計模型方法盡管對作物產(chǎn)量形成的機理解釋性不強，但其結構簡單、可操作性強，契合作物估產(chǎn)的業(yè)務化應用需求，因此仍然是當前作物遙感估產(chǎn)的常用方法。

經(jīng)驗模型以數(shù)學統(tǒng)計分析方法為基礎，包括線性模型和非線性模型。而糧食產(chǎn)量受眾多不確定因素的影響，是一個復雜的非線性系統(tǒng)，所以非線性估產(chǎn)模型的精度通常高于線性模型［6］。支持向量回歸（Support Vector Regression，SVR）方法是近年來廣泛應用的一種非線性回歸模型，其基本思想是，找到一個回歸平面，讓一個集合內(nèi)所有數(shù)據(jù)到該平面的距離最近［7］。相較于神經(jīng)網(wǎng)絡等機器學習方法，它的擬合和泛化推廣能力更優(yōu)異［8-9］，在解決糧食估產(chǎn)等非線性、小樣本的實際問題上具有較大 優(yōu)勢。

構建估產(chǎn)模型多采用與糧食產(chǎn)量具有高度相關性的參數(shù)，包括植被指數(shù)、太陽光合有效輻射、生物量等［10］，其中應用最廣泛的是能夠反映植被綠度狀態(tài)的歸一化植被指數(shù)（Normalized Difference Vegetation Index，NDVI）［11-13］。然而，除植被綠度狀態(tài)以外，影響農(nóng)作物生長的土壤水分、氣溫等指標也與產(chǎn)量的形成緊密關聯(lián)。為此，文章綜合考慮作物生長過程中的水分脅迫指標、作物生長狀態(tài)指標，選取條件植被溫度指數(shù)（Vegetation Temperature Condition Index，VTCI）、葉面積指數(shù)（Leaf Area Index，LAI）作為模型輸入變量，采用支持向量回歸方法建立遙感估產(chǎn)模型，并對該方法的精度進行評價，為及時、準確的區(qū)域小麥估產(chǎn)提供新的解決途徑。

1 材料與方法

1.1 研究思路

回歸的基本思路是利用有限的樣本數(shù)據(jù)，建立輸入和輸出之間隱含的函數(shù)關系。該文在MATLAB環(huán)境下，通過LIBSVM3.23工具包［14］，采用支持向量回歸方法，基于所選的核函數(shù)和最優(yōu)參數(shù)構建產(chǎn)量估測模型，找到各影響因素對產(chǎn)量最佳的映射關系，利用實際產(chǎn)量數(shù)據(jù)對模型模擬結果進行精度評價，技術路線如圖1所示。

1.2 研究區(qū)概況

研究區(qū)為陜西省中部關中平原（圖2a）的5個市，西安市、寶雞市、咸陽市、渭南市、銅川市，地處106°18′E～110°35′E，33°35′N～35°52′N（圖2b），總面積5.53萬km2。 關中平原三面環(huán)山，渭河橫跨平原，地勢西高東低，中部較為平坦寬闊，平均海拔約400 m。地屬大陸性季風氣候，多年平均氣溫13.9℃，多年平均降水量595.3 mm。該區(qū)域氣候溫潤、水系縱橫、土質疏松，是我國重要的糧食產(chǎn)區(qū)，也是重要的麥、棉產(chǎn)區(qū)，小麥占其耕地面積的50%左右。當?shù)馗髦贫葹閮赡耆欤魑锓N植模式主要是在灌溉區(qū)內(nèi)進行冬小麥、夏玉米輪作，在旱作區(qū)域進行夏季休閑式的冬小麥 連作［15］。

圖2 研究區(qū)地理位置Fig.2 Location of study area

1.3 數(shù)據(jù)源

該文采用的遙感數(shù)據(jù)是研究區(qū)內(nèi)2011—2016年冬小麥主要生育期（3月上旬至5月下旬）的Aqua-MODIS日地表反射率產(chǎn)品（MYD09GA）、日地表溫度產(chǎn)品（MYD11A1）、葉面積指數(shù)產(chǎn)品（MCD15A3）。所用的產(chǎn)量數(shù)據(jù)來源于2011—2016年的《陜西省統(tǒng)計年鑒》。將各市的小麥種植面積以及總產(chǎn)量匯總，表1中列出了計算后得到的對應年份各市冬小麥單產(chǎn)數(shù)據(jù)。

表1 2011—2016年研究區(qū)5市冬小麥單產(chǎn)Table 1 Yield per hectare of winter wheat of research area in five cities from 2011 to 2016 kg/hm2

1.4 研究方法

1.4.1 VTCI計算

條件植被溫度指數(shù)（VTCI）是基于歸一化植被指數(shù)（NDVI）和地表溫度（LST）特征空間呈三角形區(qū)域分布的特點提出的，能體現(xiàn)植物生長過程中受到的水分脅迫，是一種近實時的干旱監(jiān)測方法，計算方法為［16］：

式（1）～（3）中，LSTmax（NDVIi）、LSTmin（NDVIi）分別被稱作熱、冷邊界，表示在研究區(qū)內(nèi)，當NDVIi值為某一特定值時，所有像素地表溫度的最大值和最小值；LST（NDVIi）表示某一像素的NDVI值為NDVIi時的地表溫度；a，b，a′，b′是待定系數(shù)，由研究區(qū)內(nèi)的NDVI和LST散點圖近似得到。

計算VTCI采用的遙感數(shù)據(jù)是研究區(qū)內(nèi)2011—2016年冬小麥主要生育期（3月上旬至5月下旬）的Aqua-MODIS日地表反射率產(chǎn)品（MYD09GA）、日地表溫度產(chǎn)品（MYD11A1）。利用MODIS重投影工具MRT（MODIS Reprojection Tools）對原始數(shù)據(jù)進行拼接、重采樣、投影轉換等預處理操作后，得到研究區(qū)日LST和日地表反射率數(shù)據(jù)。利用日地表反射率數(shù)據(jù)計算日NDVI，用最大值合成方法分別生成逐像素的旬LST和旬NDVI最大值合成數(shù)據(jù)，并按照VTCI定義計算出旬VTCI。

將冬小麥的主要生育期劃分為4個生育時期，結合關中平原冬小麥越冬后的生長情況，將返青期定為3月上旬至3月中旬，拔節(jié)期為3月下旬至4月中旬，抽穗—灌漿期為4月下旬至5月上旬，乳熟期為5月中旬至5月下旬［17］。計算每個生育期內(nèi)各旬VTCI的平均值作為該生育期的VTCI值，計算研究區(qū)2011—2016年5市各生育期的VTCI。

1.4.2 LAI計算

采用2011—2016年冬小麥主要生育期內(nèi)的葉面積指數(shù)產(chǎn)品（MCD15A3）計算LAI，這一產(chǎn)品每4 d合成1次，時間分辨率較高，空間分辨率為500 m，有利于農(nóng)作物長勢監(jiān)測。同樣對影像進行投影轉換、裁剪等預處理，利用S-G濾波器消除云和大氣等噪聲，改善數(shù)據(jù)質量。然后采用取最大值的方法逐像素地生成旬LAI，取各生育期內(nèi)所包含的多旬LAI的均值分別作為該生育期LAI，實現(xiàn)LAI、VTCI時間分辨率的 統(tǒng)一。

1.4.3 SVR模型構建

①支持向量回歸

支持向量回歸的基本思想是采用一個非線性映射把數(shù)據(jù)映射到高維特征空間，然后在高維特征空間構造回歸估計函數(shù)，再映回到原空間。通過定義適當?shù)暮撕瘮?shù)k（x，xi）來實現(xiàn)非線性變換，無需得到確切的映射函數(shù)Φ（x），給計算帶來了極大方便［18］。回歸方程為：

式（4）中，ai-a*

i是各支持向量的系數(shù)。

②格式轉換及歸一化

通過運行FormatDataLibsvm.xls，將數(shù)據(jù)格式轉化為LIBSVM規(guī)定的輸入格式。為避免數(shù)據(jù)之間的量級差別，對VTCI、LAI、單產(chǎn)數(shù)據(jù)分別采用歸一化處理，歸一化公式 為：

式（5）中，Xmax和Xmin分別為數(shù)據(jù)中的最大值和最小值，經(jīng)過歸一化消除了模型輸入由于量綱和單位不同造成的影響，使樣本數(shù)據(jù)更加適應回歸建模與分析。

③核函數(shù)選擇

支持向量回歸包括線性和非線性，在線性回歸的基礎上引入核函數(shù)，得到非線性回歸。支持向量回歸將高維空間內(nèi)積運算化簡成了從低維空間進行輸入的核函數(shù)計 算［19-20］，選擇的核函數(shù)的類型決定了支持向量回歸的很多特性，從某種程度上直接影響著回歸模型的擬合能力和預測精度。因此，核函數(shù)的選擇是用支持向量回歸預測糧食產(chǎn)量的一個關鍵步驟。

常用的核函數(shù)有多項式核函數(shù)、徑向基（RBF）核函數(shù)（即高斯核函數(shù)）、Sigmoid核函數(shù)等，選用不同的核函數(shù)可以構造出不同的支持向量回歸模型。

許多實驗研究表明［21-22］，當缺少先驗知識時，選用徑向基核函數(shù)訓練建模的效果較好，所得模型的總體性能較高，而且徑向基核函數(shù)只含一個未知參數(shù)σ，易于進行優(yōu)化，所以該研究選擇的核函數(shù)為徑向基核函數(shù)。公式為：

④參數(shù)尋優(yōu)

在用支持向量回歸解決實際問題時，參數(shù)的選擇對回歸模型的性能影響很大。應用RBF核函數(shù)時，懲罰因子C和參數(shù)σ與回歸模型的學習性能息息相關。懲罰因子C用于控制模型對誤差范圍以外的樣本的懲罰程度，也可以說是對離群點的重視程度，懲罰因子過高會造成過擬合。徑向基核函數(shù)的參數(shù)σ反映了訓練樣本數(shù)據(jù)的范圍或分布，確定局部鄰域的寬度。

該研究選用LIBSVM中的ε-SVR方法，用網(wǎng)格搜索和n-fold交叉驗證的方法進行參數(shù)尋優(yōu)，確定C和σ值。當有不同的C和σ都對應最高的精度時，把參數(shù)C最小的那組C和σ作為最佳參數(shù)，避免懲罰參數(shù)C太高造成過學習狀態(tài)。

1.4.4 模型精度評價

該文選取決定系數(shù)（Coefficient of Determination，R2）、均方根誤差（Root Mean Square Error，RMSE）、平均絕對百分比誤差（Mean Absolute Percentage Error，MAPE）作為模型擬合程度的評價指標，具體公式為：

式（7）～（9）中，i表示第i個樣本點數(shù)據(jù)，Yi為第i個實際產(chǎn)量，單位為kg/hm2；Ei為根據(jù)模型擬合的冬小麥產(chǎn)量估算值，單位為kg/hm2；為實際產(chǎn)量的均值，單位為kg/hm2；為模型估測產(chǎn)量的均值，單位為kg/hm2。

2 結果與分析

2.1 回歸結果及誤差分析

研究以2011—2015年研究區(qū)的數(shù)據(jù)作為訓練集，采用SVR算法構建冬小麥估產(chǎn)模型，圖3為2011—2015年研究區(qū)各市的實際產(chǎn)量和回歸產(chǎn)量，可以看出預測產(chǎn)量和對應的實際產(chǎn)量波動的趨勢大致相符。以實際產(chǎn)量和估測產(chǎn)量分別為橫縱坐標繪制散點圖（圖4），添加趨勢線和45°參考線，趨勢線對比參考線偏離較小，可見模型對測試集的擬合效果較好。計算模型對訓練集樣本的預測精度，回歸模型具有較高的決定系數(shù)R2=0.88，估測產(chǎn)量與實際產(chǎn)量間的RMSE為336.39 kg/hm2，MAPE為6.12%。

對訓練集的擬合學習說明了SVR構建的模型能夠以較小的誤差擬合歷史數(shù)據(jù)。采用構建好的模型預測2016年5市冬小麥單產(chǎn)，表2為預測結果及誤差，可見5市數(shù)據(jù)共同構建的模型進行單產(chǎn)估測時有一定誤差，其中面積較小的銅川市的相對誤差較小，僅有1.28%。分別對2011—2015年5市冬小麥產(chǎn)量數(shù)據(jù)進行建模預測（表3），預測結果的RMSE平均值為232.63 kg/hm2，其中寶雞市、渭南市的MAPE低至2.98%、3.37%，5市MAPE的平均值為5.74%，可見模型精度較高。

圖3 2011—2015年5市冬小麥實際產(chǎn)量與回歸產(chǎn)量對比Fig.3 Comparison of actual yield and regression yield from 2011 to 2015

圖4 訓練集擬合情況Fig.4 Fitting effect of training set

表2 2016年5市冬小麥回歸預測結果及誤差Table 2 Regression prediction results and errors

表3 5市分別構建回歸模型的預測結果及精度Table 3 Prediction results and accuracy of regression models constructed by different municipalities

2.2 參數(shù)敏感性分析

通過觀察不同參數(shù)對預測精度的影響，判斷對預測精度提高幫助最大的參數(shù)。圖5、圖6給出了C和σ值變化時3個誤差評價指標的相應變化。結果顯示，發(fā)生同等程度的變化時，核參數(shù)σ引起的MAPE、RMSE、R2的改變都比懲罰因子C引起的改變更大。因此，核參數(shù)σ比懲罰因子C更重要。對2個參數(shù)進行耦合分析，發(fā)現(xiàn)該研究中當懲罰因子C為2、核參數(shù)σ為0.354時模型精度最高。

圖5 不同參數(shù)C條件下的回歸模型精度變化Fig.5 The accuracy change of regression model under different parameter C

3 結論與討論

3.1 結論

該研究計算了2種與植被生長狀況息息相關的參數(shù)——VTCI、LAI，將冬小麥各個生育期的VTCI和LAI作為輸入變量，劃分訓練集和測試集，用支持向量回歸算法對這些數(shù)據(jù)進行學習和訓練，從而找到其內(nèi)部隱含的函數(shù)關系。尋找最優(yōu)的模型參數(shù)，建立冬小麥各生育期的VTCI、LAI與冬小麥產(chǎn)量之間的支持向量回歸估產(chǎn)模型。

圖6 不同參數(shù)σ條件下的回歸模型精度變化Fig.6 The accuracy change of regression model under different parameter σ

支持向量回歸模型較好地反映了輸入的影響因素VTCI、LAI同冬小麥單產(chǎn)之間的復雜的非線性映射關系，建立的回歸模型決定系數(shù)達到0.88，模型擬合較為理想，此外，模型參數(shù)的選取對模型精度有明顯的影響，徑向基核參數(shù)σ比懲罰因子C對預測結果影響更大。

該文采用支持向量機方法對關中平原的冬小麥產(chǎn)量進行預測，在樣本較少的情況下，建立的回歸模型精度較高，能夠推廣到其他地區(qū)，也可以為其他作物的產(chǎn)量預測研究提供參考，該方法在糧食產(chǎn)量預測領域有良好的應用前景。

3.2 討論

該文僅選用了VTCI和LAI 2個指標，只能在一定程度上描述作物本身的生長狀態(tài)。實際上，影響糧食產(chǎn)量的外部因素眾多，還包括很多非遙感信息，如自然災害、勞動主體的積極性、農(nóng)業(yè)科技水平、農(nóng)業(yè)的基礎設施等。該研究缺少足夠的先驗知識，模型設置過于理想化，并未充分考慮多層次的因素，若要得到更為科學合理的估產(chǎn)結果，仍需繼續(xù)探索與糧食產(chǎn)量密切相關的指標。

該研究核函數(shù)的選定依照前人的經(jīng)驗，而無充分的理論支撐，因此核函數(shù)的選用是未來需要進一步研究的內(nèi)容。隨著新技術的不斷涌現(xiàn)，支持向量回歸和其他理論方法的結合值得深入發(fā)掘，例如參數(shù)尋優(yōu)可采用粒子群優(yōu)化算法［23］、遺傳算法，核函數(shù)可選用混合核函數(shù)［24］等。

此外，該研究采用統(tǒng)計年鑒的產(chǎn)量數(shù)據(jù)進行驗證，但產(chǎn)量數(shù)據(jù)不可避免地存在著差錯或漏報［25］, 因此如何有效驗證模型精度也是有待解決的問題之一。