基于CNN-S-GPR的寧夏枸杞高光譜影像估產方法

劉立波 王 濤 張 鵬

(寧夏大學信息工程學院， 銀川 750021)

0 引言

枸杞作為寧夏特產，是國家地理標志產品，其產量是重要的經濟信息，預測枸杞年產量對于種植管理和政府決策是一項難度較大但必不可少的研究。傳統的農作物估產采用人工區域調查方法，從農學、氣象學等不同角度建立作物估產模式。該方法速度慢、工作量大、成本高，且不利于時空動態監測。近年來，遙感技術因其覆蓋范圍廣、重返周期短、獲取成本低等優勢，被廣泛運用于農作物估產，成為遙感與農業交叉的研究重點[1]。

目前，遙感估產方法按模式不同主要分為2類：基于機器學習和基于深度學習的估產方法[2]。基于機器學習的方法將歸一化植被指數(Normalized difference vegetation index, NDVI)、增強型植被指數(Enhanced vegetation index, EVI)等植被指數作為表征作物產量的重要信息，進而采用機器學習方法建立模型預測農作物產量。文獻[3-5]利用NDVI、EVI等植被指數，基于Lasso、支持向量回歸(Support vector regression, SVR)等回歸算法構建了產量預測模型，達到了比較高的估產精度。但這類方法NDVI、EVI等指數僅利用近紅外波段和紅外波段2個波段計算得到，忽略了其余波段重要信息，不利于估產精度的提升，對作物長勢反映不理想。所以，文獻[6-10]提出了作物生長模型(Aquacrop、WheatSM、WOFOST等)，根據氣象條件、土壤條件、葉面積指數(Leaf area index, LAI)以及作物栽培管理措施，定量描述作物生長、發育、籽粒形成等動態過程。這類模型雖然模擬作物生長機理性強、估產精度高，但過分依賴實測數據和人工經驗，使工作量和成本大幅度上升，具有一定局限性。基于深度學習的估產方法，利用神經網絡非線性映射能力擬合樣本與觀測值之間的關系，從而預測出產量。如高云[11]采用DNN(深度神經網絡)準確預測了春小麥產量。針對DNN等前饋神經網絡難以處理時序數據的問題，文獻[12-13]采用長短期記憶網絡(Long short-term memory, LSTM)對長期信息依賴進行處理，大大提升了估產精度，但仍需大量手工時間序列特征參與計算。文獻[14-15]采用卷積神經網絡(Convolutional neural networks, CNN)直接從時序影像中提取特征，摒棄了原有手工制作特征，進而利用全連接網絡預測了中國北方冬小麥產量，并取得了很好的效果，但是仍存在以下問題：①基于CNN的遙感估產方法雖然簡化了特征提取操作，卻未對通道間依賴關系進行建模，即無法自動獲取每個特征通道的重要程度，以強調有效信息并抑制無效信息。②CNN回歸模型僅能提取遙感影像時間序列特征，卻無法整合特定位置空間下對作物生長有影響的特征，因此具有一定局限性。

綜上所述，本文在CNN模型的基礎上，提出一種基于CNN-S-GPR的高光譜影像估產模型。通過CNN提取影像特征，融合通道注意力機制SENet(Squeeze-and-Excitation network)對CNN卷積層通道進行加權，以對通道間的依賴關系建模。接著在網絡中引入高斯過程回歸(Gaussian process regression, GPR)方法整合影像特征和空間位置特征，以寧夏回族自治區中寧縣、賀蘭縣等16個縣為估產單元，基于多級Modis高光譜影像數據，對各縣域枸杞產量進行預測。

1 研究區概況

研究區域位于寧夏回族自治區(35°14′～39°23′N， 104°17′～107°39′E)，地處黃河水系，屬北溫帶大陸性季風氣候區，是全國主要枸杞種植地[16]，研究區概況如圖1所示。寧夏枸杞種植主要分布在中南部地區，北部地區種植面積較少。此外，枸杞每年的休眠期從當年11月至次年4月，生長期為7個月，即萌芽期為4月上旬，開花初期為5月上旬，果熟期為6月中旬，連續開花，連續結果，直至10月下旬落葉。

圖1 研究區概況圖Fig.1 Overview of study area

2 研究方法

2.1 數據準備及預處理

2.1.1數據準備

本文采用的實驗數據由寧夏回族自治區遙感影像、枸杞種植區域矢量圖和年際枸杞產量3類數據組成。其中，遙感影像采用Modis高光譜影像數據，來源于EARTHDATA網站，行列號為h26v04，其波段信息豐富，光譜分辨率可達納米級，因而可提取農作物的反射峰、吸收谷等特征，更好地刻畫農作物長勢及信息差異；同時，Modis遙感數據具有高時間分辨率特性，可獲得充足的多時相估產數據。根據枸杞每年4月上旬萌芽，10月下旬落葉，選擇每年第97天至第297天，生長季總計201 d，26幅時相影像(時間分辨率為8 d)來構建CNN-S-GPR網絡模型訓練樣本。每個時相影像包括MOD09A1、MOD13A1、MYD11A2、MCD15A2H 4類Modis產品數據，共13個波段，影像及矢量數據如表1所示。MOD09A1為地表反射率數據，其1～7波段可反映作物生長環境及土壤含水率；MOD13A1為植被指數產品數據，其中NDVI和EVI 2個波段可用于作物產量預測[17-18]；MYD11A2為地表溫度數據，包含白天和夜間地表溫度波段，與作物冠層溫度密切相關[19]；MCD15A2H為葉面積指數和光合有效輻射產品數據，LAI和FPAR(光合有效輻射)是大量作物生長模型的基礎，與作物產量之間的關系更為直接[20]。此外，枸杞種植區域矢量數據來源于寧夏農林科學院，為shp文件，由工作人員于2017年實地考察記錄枸杞種植區并繪制而成；年際枸杞產量數據由寧夏回族自治區統計局提供，包括2010—2019年寧夏22個縣(縣級市區)枸杞種植面積和實際產量，因統計數據的不連續性，篩除西吉、隆德等6個縣后，保留了同心、中寧等16個縣作為估產區域(圖1)，2019年產量數據如表2所示，研究區總面積為27 960 hm2，總產量為94 843 t，平均產量為3.39 t/hm2。

表1 影像及矢量數據Tab.1 Image and vector data

表2 2019年產量數據Tab.2 Yield data of 2019

2.1.2數據預處理

數據預處理過程主要分為4步：①利用MRT工具將Modis遙感影像和枸杞種植區域掩膜數據重投影為基于WGS-84橢球體的UTM投影，并將MYD11A2影像和枸杞矢量數據重采樣的空間分辨率為500 m，保證數據空間位置的一致性。②由于MOD13A1時間分辨率為16 d，其余Modis數據產品為8 d，為了保證時序的完整性，采用上下影像求算術平均值的方法對影像進行補充，如將MOD13A1第97天和第113天影像的均值對第105天進行補充。③根據寧夏行政區劃圖分別將枸杞種植區域掩膜和Modis遙感影像數據裁剪為16個縣枸杞種植區域掩膜及216 320幅遙感影像，并提取各縣中心經緯度。④利用不同縣域枸杞種植矢量圖提取Modis高光譜遙感數據的感興趣區域，以中寧縣枸杞種植感興趣區域提取為例，圖2為中寧縣枸杞種植區域矢量圖，綠色區域為枸杞種植地塊。圖3為經步驟③處理所得的中寧縣Modis高光譜影像MOD09A1產品數據，調用GDAL(Geospatial data abstraction library)庫中的warp函數以實現圖2對圖3的裁剪，從而得到中寧縣枸杞種植感興趣區域圖(圖4)，并對同一天不同產品進行波段提取和融合處理后得到16 640幅包含13個波段的影像。

圖2 中寧縣枸杞種植區域矢量圖Fig.2 Wolfberry planting area in Zhongning

圖3 中寧縣Modis高光譜影像MOD09A1圖像Fig.3 Modis HSI MOD09A1 in Zhongning

圖4 中寧縣枸杞種植感興趣區域圖Fig.4 ROI for wolfberry planting in Zhongning

2.2 數據集構建

2.2.1直方圖降維與歸一化

鑒于訓練樣本的稀疏性，無法采用端到端方式訓練估產模型。因此，本文將各縣枸杞種植感興趣區域的高維遙感影像映射到直方圖中，使得產量與不同灰度的像素數量相關，而與影像像素位置無關。然而遙感影像中的像素值是離散的，假設每個波段最多可以取b個不同值，那么對于d個波段的影像將得到具有bd個像素區間的直方圖。如在可見光圖像中，每個波段可以取b=256個不同值，波段數d=3，則直方圖像素區間個數約為1.67×107，顯然將b個不同的取值離散至b個區間不可取，會導致維數災難問題。因此，本文劃定x(x

(1)

式中Hi——直方圖歸一化結果

hi、hj——像素直方圖中每個區間像素數

2.2.2時間序列融合與維度轉換

經上述直方圖降維與歸一化處理后，使得每個影像變成尺寸為32×13的矩陣。為了進一步融合高光譜影像不同時相上的枸杞生長信息，對上述直方圖統計結果在時間序列上進行融合，從而得到時序樣本(x軸為直方圖區間維度，y軸為圖像波段維度，z軸為時間維度)。然而，不同于3D CNN提取3維時序數據樣本特征的方式，本文采用2D CNN獲取3維時序數據樣本的植被長勢信息。為了提取時序樣本中枸杞生長時間特征，需對圖像波段維度和時間維度進行轉換，最終形成32×26×13的矩陣作為網絡的輸入。將高、中、低產量不同波段的數據融合結果(x軸為像素區間，y軸為時間維度，z軸為影像波段維度)進行可視化對比，如圖5所示。從圖5可以看出，在高產、中產和低產中波段1、波段7、NDVI、EVI等波段明顯存在視覺差異，表明特征抽取網絡能夠提取到有效特征。

圖5 數據融合可視化結果Fig.5 Visualization results of data fusion

2.3 枸杞估產模型構建

2.3.1網絡結構設計

高光譜影像波段寬而廣，包含豐富的光譜信息和圖像信息，能較好地反映作物長勢。本研究首先通過構建CNN對影像特征進行提取，接著為了將注意力集中在對枸杞估產具有重要作用的通道上，在CNN后融合了通道注意力(SENet)模塊，以建模通道間的依賴關系；最后引入高斯過程回歸方法(GPR)整合影像特征和空間位置特征，對各縣枸杞年際產量進行預測。綜上所述，采用CNN作為骨干網絡，通過引入SENet模塊和GPR方法，提出一種基于CNN-S-GPR網絡的高光譜影像估產模型用于枸杞估產，網絡結構如圖6所示。其中，第1層為輸入層，為32×26×13的矩陣。第2層為特征提取層，由6個卷積層構成，其中淡青色為步長為1的卷積層，藍色為步長為2的卷積層，邊緣采用1個像素零填充，卷積核數量分別為128、256、256、512、512、512，卷積核大小全為3×3。在每個卷積層上進行批歸一化處理后，采用ReLU函數激活，并加入Dropout層和L2正則化，以避免“梯度彌散”和模型過擬合問題[21]。第3層為SENet模塊，對骨干網絡最后一層進行Squeeze和Excitation 2個操作，以學習到不同通道的權重。最后為全連接層，將提取特征全連接為2 048維向量，并融合位置特征后作為高斯過程回歸模型的輸入。

圖6 CNN-S-GPR網絡結構圖Fig.6 CNN-S-GPR network structure diagram

2.3.2高斯過程回歸

基于預處理后的遙感影像數據集，首先采用CNN網絡對2019年寧夏16個縣的枸杞產量進行預測；然后利用半變異函數將各縣距離與產量絕對誤差相關聯，將計算結果進行可視化，如圖7所示。由圖7可知，當空間上距離越接近時，半變異函數值越小。由此可見，枸杞種植空間位置對枸杞估產精度影響較大。若將遙感影像特征與枸杞種植空間位置結合起來，將對提升枸杞估產精度具有重要意義。

圖7 半變異函數計算結果可視化Fig.7 Visualization of semivariogram results

因此，采用高斯過程回歸預測枸杞產量，不僅整合了影像特征和空間位置特征，而且還解決了小樣本模型訓練問題[22]，其線性模型為

y(x)=h(x)Tβ+
(β～N(w,σwI),～N(0,σ2))

(2)

式中y(x)——觀測值(產量預測值)

h(x)——SENet加權后的CNN特征向量

w——h(x)的權重向量

σw——超參數

I——主對角線為1的方陣

由于觀測值y的絕對誤差與x分布有關，即2個x距離越近，y的絕對誤差越小。所以引入高斯過程f(x)替換了式(2)中的殘差函數。不同于深度學習回歸模型，該方法隱式地將輸入特征映射到不同維度的空間，從而求解到模型最優權重，GPR方法直接對函數建模，可顯式整合影像特征和空間位置特征，公式為

y(x)=h(x)Tβ+f(x)
(f(x)～N(0,k(x,x′)))

(3)

式中f(x)——核函數為k(x,x′)、零均值的高斯過程函數

核函數k(x,x′)采用徑向基函數核，加入高斯噪聲后，其公式為

(4)

式中gloc、g′loc——各縣歸一化后的中心經緯度

rloc、σ、σn——超參數

經上述推導后，整合影像特征和位置特征的GPR公式為

(5)

為了使高斯過程回歸擬合效果達到最優，基于交叉驗證，采用網格搜索進行最優超參數值的搜索，最終將σ設為1，σw為0.01，σn為0.32，rloc為0.5。

2.4 損失函數

在回歸問題中，損失函數常被用于描述模型真實值與預測值的接近程度，常見的損失函數有L1損失和L2損失，具體公式為

(6)

(7)

式中n——樣本總數

yi——第i個樣本的真實值

L1loss——L1損失L2loss——L2損失

L1損失表示真實值和預測值之間絕對誤差的平均值。然而，由于L1損失梯度不變，即使損失很小，梯度也非常大，不利于模型收斂。所以本文采用梯度變化的L2損失作為損失函數，以表示真實值和預測值之間的距離平方和。

2.5 評價指標

采用平均相對誤差(Mean relative error, MRE)、均方根誤差(Root mean square error, RMSE)和決定系數(Coefficient of determination,R2)作為評價指標，以驗證基于CNN-S-GPR模型準確性。

3 實驗結果與分析

3.1 模型訓練

通過對寧夏16個縣2010—2019年Modis遙感影像、枸杞種植區域矢量圖和枸杞產量統計數據進行預處理、直方圖降維、數據融合和維度轉換后，共獲得160個數據樣本。實驗將2010—2016年共112個樣本作為訓練集，2017—2018年32個樣本作為驗證集，2019年16個樣本作為測試集，訓練估產模型，訓練集、驗證集和測試集劃分比例為7∶2∶1。

本文實驗環境為Ubuntu 16.04.12操作系統，深度學習框架選用 Pytorch，軟件環境為Python 3.6.10和Cuda 9.0，GPU 采 用 NVIDIA GeForce RTX 2080 Ti。模型采用分批訓練，每個批次隨機選擇26個樣本訓練模型，共訓練30 000次，并利用Adam優化器優化。學習率初始化為0.001，為防止模型因學習率較大導致跳過全局最優點的問題，采用分段常數衰減策略降低學習率，即當訓練次數分別達到4 000和15 000時，學習率降為原來的10%。同時設置patient值，當驗證集損失累計10次未發生改變時，提前終止訓練。

訓練效果如圖8所示，模型逐漸收斂，直至迭代次數為20 000時，訓練集和測試集損失達到最小并趨于穩定。保存模型并在測試集上進行測試，測試結果如圖9所示，R2達到0.91。同時采用直線回歸擬合真實值與預測值間的關系，直線斜率為0.81，接近于1。由此可見，經直方圖降維和歸一化處理后的高光譜影像信息丟失較少，且通過高斯過程回歸將影像信息與位置信息整合后，能夠準確擬合枸杞高產量和低產量樣本。

圖8 訓練效果Fig.8 Training effect diagram

圖9 測試結果Fig.9 Test result graph

3.2 實驗對比

3.2.1消融實驗

為了驗證SENet模塊和GPR方法對本文所提模型性能的提升效果，在自建數據集上，通過以下3個不同的模型進行消融實驗：①用CNN(baseline)表示原卷積神經網絡模型。②用CNN-SENet表示在模型①中融合通道注意力后的網絡。③用CNN-S- GPR表示在模型②中添加GPR后的網絡，實驗結果如表3所示。

表3 消融實驗結果Tab.3 Results of ablation experiments

從表3可以看出，各個改進后的模型在性能上均有一定程度的提升。其中，僅融合通道注意力后的網絡(CNN-SENet)與基線網絡模型相比，在MRE和RMSE上分別下降了0.26個百分點和15.39 t，且R2達到0.85，說明產量預測值中85%變異可由加權圖像特征解釋，通過通道注意力機制對CNN特征進行加權后，強調了有效信息，能更好地表征枸杞生長機理。融合了通道注意力和高斯過程回歸后的網絡(CNN-S- GPR)，針對影像特征及枸杞種植位置特征進行整合，可避免因缺少空間位置信息產生的估產誤差，MRE、RMSE和R2分別達到13.57%、776.58 t和0.91，比其他2種方法效果更好。與CNN-SENet相比， MRE和RMSE分別下降了0.69個百分點和67.26 t，R2提升了0.06；與基線模型相比， MRE和RMSE分別下降了0.95個百分點和82.65 t，R2提升了0.08。由此驗證了本文方法的有效性。

3.2.2對比實驗

為了使實驗結論更具有準確性與說服力，將本文模型與Lasso、嶺回歸(Ridge regression，RR)、SVR、DNN和 LSTM 5種近年來用于估產的主流回歸模型進行對比。其中，Lasso、RR和SVR為機器學習回歸方法；DNN和LSTM為深度學習回歸方法，均有3個隱藏層和256個神經元。與高斯過程回歸不同，對比的回歸模型輸入均采用手工特征，由高光譜影像感興趣區域的平均NDVI、EVI、地表溫度、LAI和FPAR值序列組成，共計16 640條數據。對比實驗結果如表4所示。

由表4可知，本文方法在自建數據集上，MRE、RMSE和R2均優于近年來用于估產的主流回歸方法。與表中對比機器學習回歸方法不同，本文方法一方面運用時間序列融合方法融合了不同時相上的時間信息，另一方面采用CNN提取了NDVI、EVI等手工特征不具有的其他波段信息，使特征豐富度進一步增強，所以本文R2高于其他方法。雖然LSTM也能很好地處理時間序列，但是忽略了空間位置對枸杞生長的影響，所以GPR方法在整合影像特征和空間位置特征后，MRE和RMSE相對于LSTM下降了0.44個百分點和52.48 t。總之，與上述所提方法相比，本文方法提升了枸杞估產的準確性。

表4 對比實驗結果Tab.4 Comparative experimental results

3.3 各縣估產準確性分析

為了更進一步評估CNN-S-GPR模型在各縣區枸杞估產的準確性，對2019年寧夏16個縣(縣級市區)枸杞統計產量與預測產量進行比較，各縣估產準確性分析如圖10所示。誤差較大的區域主要為靈武市、鹽池縣、大武口區和賀蘭縣，誤差產生的原因可能是沒有足量的一年一度枸杞種植區域矢量圖，枸杞種植面積的年際變化性使得模型預測產量數據偏差較大；此外，靈武市、鹽池縣和大武口區種植面積過小且存在混合像元現象，CNN所提取的枸杞生長信息過少，給模型準確預測帶來了一定的困難。與之相反，模型在中寧縣、同心縣、沙坡頭區等縣表現優秀，統計產量曲線與預測產量曲線幾乎重合。總之，CNN-S-GPR模型在研究區的估產精度較高，可實現寧夏枸杞年際估產。

圖10 各縣估產準確性分析Fig.10 Analysis chart of yield estimation accuracy by counties

4 結束語

提出的CNN-S-GPR高光譜影像年際作物估產模型，在CNN特征提取網絡后融合通道注意力機制，表征不同通道的重要性，增強了枸杞生長機理性模擬，在此基礎上，引入高斯過程回歸方法，解決了小樣本估產問題，同時整合影像特征和空間位置特征，使估產特征更加豐富。與其他估產模型相比，該模型平均相對誤差和均方根誤差下降了0.44～0.95個百分點和52.48～82.65 t，且R2達到0.91，驗證了本文方法的有效性。