尼格拉·吐爾遜, 蘇磊·乃比, 高 健, 沈江龍, 鄭江華*, 余丹林
1. 新疆大學資源與環境科學學院, 新疆 烏魯木齊 830046 2. 新疆大學數學與系統科學學院, 新疆 烏魯木齊 830046 3. 新疆林業科學院現代林業研究所, 新疆 烏魯木齊 830063 4. Department of Earth and Environmental Studies, Montclair State University, New Jersey 07043, USA
葉綠素含量是棗樹光合作用能力、 生長狀況、 營養狀況的指示劑[1], 通常采用便攜式葉綠素計(SPAD-502)測定植物葉片SPAD值來直接表征植物葉綠素含量的相對大小, 但使用過程中需要將葉片反復插入測量, 難以用于大范圍的葉綠素檢測, 研究表明SPAD值能與無損、 無污染、 價格低的高光譜遙感數據準確對應, 近年來成為葉綠素含量估算的強有力工具[2]。
20世紀90年代, Pinar[3]和 Blackburn[4]等研究得到葉綠素與高光譜波段之間的相關關系。 隨后, 許多學者在高光譜估算葉綠素模型方面開展了大量的研究, 杜華強基于高斯核函數變換的偏最小二成回歸模型建立了馬尾松針葉葉綠素含量與光譜反射率及9個特征參數之間的預測模型, 其精度遠大于傳統線性回歸模型[5]。 劉京等用實例證實了支持向量機具有更好的SPAD值反演效果[6]。 馮海寬等基于特征光譜參數, 利用隨機森林模型較好的估算了蘋果葉片葉綠素[7]。 李曉麗等證實了最小二乘支持向量機(least sqares support vector regression, LS-SVR)在植物參數估算方面具有較好效果[8]。
上述研究常選用相關系數較高的波段或者植被指數建模使得變量選擇隨機、 單一、 缺乏定量化, 模型估算能力低下。 本文通過CP統計量在預測角度選擇重要性較高的自變量, 篩選重要程度高的特征波段(characteristic band, CB)。 其次, 以往的高光譜估算應用廣泛的多元線性回歸(multiple linear regression, MLR)、 支持向量機(support vector regression, SVR)、 LS-SVR模型較多, 并沒有考慮到地理位置可能對葉片SPAD值產生的影響。 2017年Hwang和Shim對于LSSVM模型加入地理位置影響, 提出了地理加權最小二乘支持向量機模型[9](GWLS-SVM), 證實了其估計精度顯著高于傳統的GWR、 LS-SVR模型。 本研究對于GWLS-SVR模型是否適用于葉綠素含量估算, 能否在紅棗樹葉片葉綠素估算中得到較好效果還需要進一步的驗證。 用CP統計量選擇特征波段, 計算若羌紅棗樹葉片SPAD值的全局莫蘭指數, 分析紅棗樹葉片SPAD值分布是否與空間位置有關, 再運用GWLS-SVR模型, 將建模結果與傳統模型進行對比分析, 檢驗并比較模型的擬合效果。
研究區位于中國新疆若羌縣, 范圍在東經87°00′—89°0′、 北緯38°40′—39°30′之間, 屬暖溫帶大陸性荒漠干旱氣候, 是新疆名牌產品“若羌紅棗”種植區[10]。 于若羌紅棗果實成熟期2019年9月28日—10月2日采樣, 為了保證實驗結果的全面性和精確性, 在去除野外數據異常值后最終保留均勻覆蓋若羌縣的67個棗林樣點, 在預先設計的棗林內確定代表性棗樹1~3棵進行數據采集, 再通過手持GPS記錄地理位置信息, 共采集219條紅棗樹葉片高光譜數據和219個棗樹葉片SPAD值數據, 研究區位置和采樣樣點地理位置分布情況如圖1。

圖1 研究區位置和采樣點分布圖
紅棗樹的葉片光譜反射率在晴朗無風無云條件下于北京時間11:00—17:00使用PSR-3500便攜式地物光譜儀在野外測定, 波段范圍是350~2 500 nm, 每隔1 nm輸出一個數據, 一共2 151個光譜通道。 在選擇的代表性1~3顆棗樹上、 中、 下層各隨機采集3片葉片。 為減少誤差, 每次光譜測定之前均進行白板標定, 同時用干燥紙巾去除葉片表面浮塵, 測量時將葉片鋪平放置在反射率近似為零的黑板上, 將光纖探頭垂直固定于葉片上方約5 mm, 每個葉片樣本避開葉脈重復測量3次, 取光譜曲線的算術平均值作為該樣點的原始葉片光譜反射率。 為減少噪聲影響, 剔除1 050~2 500 nm噪聲較大波段, 并利用Origin軟件平滑去噪[11]。 另外, 導數光譜可以反映植被中生化物質的吸收引起的波形變化還能夠揭示光譜峰值的內在特征進而估算植被內部葉綠素含量信息[2]。 因此, 對原始光譜反射率(raw reflectance, RR)求光譜一階導數(first derivative of reflectance, FD)。
使用葉綠素計(SPAD-502Plus, Konica Minoita, Japan)對現場采集的多個棗樹葉片SPAD值進行測定, 測量時避開葉脈部分, 從葉柄至葉尖分段隨機測量3次, 將多個葉片測定結果取算術平均值作為該樣點SPAD值。 SPAD值測定時間與葉片光譜測定同步進行, 測定位置與葉片光譜保持一致。
本工作采用CP統計量進行變量選擇。CP統計量可以通過預測的角度選擇重要性較高的自變量。 其原理為由部分變量預測的均方誤差可能比利用所有變量進行預測的均方誤差更小, 故可以去除重要程度不是很高的變量。 其計算方式如式(1)
(1)

f(xi,Ui)=ωTφ(xi)+bi
(2)
設給定x與Ui下的權重矩陣為Wi, 則可以將回歸模型轉化成如式(3)優化問題
(3)
其中, C>0為懲罰參數, wij為用于表示Ui和Uj之間的距離的權重函數。

(4)

圖2是不同范圍SPAD值的紅棗樹葉片平均光譜反射率曲線圖。 由圖可知, 不同范圍SPAD值的紅棗樹葉片平均反射率曲線變化趨勢基本相同。 總體上, 350~750 nm波段內反射率比750~1 050 nm波段低。 在350~675 nm波段內隨著SPAD值的升高, 紅棗樹葉片平均光譜反射率降低, 光譜差異較明顯, 其中, 在500~551 nm波段范圍內反射率緩慢上升, 551 nm附近出現反射峰, 675 nm附近出現吸收谷; 675~750 nm處平均光譜反射率隨著波長呈現快速上升趨勢, 750~1 050 nm范圍內, 隨著SPAD值的升高, 平均光譜反射率升高。 紅棗樹的長勢狀態直接決定了SPAD值的大小, SPAD值也會影響紅棗樹葉片的反射率。

圖2 不同范圍SPAD值紅棗樹葉片平均光譜反射特征
為了明確紅棗樹葉片SPAD值相對應的敏感波段, 將紅棗樹葉片SPAD值和原始光譜、 光譜一階導數反射率波段做皮爾遜相關性分析。 由圖3可知, 紅棗樹葉片SPAD值和原始光譜反射率及光譜一階導數反射率緊密相關, 且都存在著極顯著相關。 對原始光譜來說, 在570~620及690~700 nm間達到相關系數峰值, 通過了0.01的顯著性水平, 相關系數分別達到-0.578及-0.561, 此波段范圍受葉綠素吸收的影響, 相關系數呈負相關, 選擇這兩組波段的原始光譜反射率作為估測棗樹葉片SPAD值的敏感波段區間。 SPAD值與光譜一階導數呈正負相關, 相關性極顯著的波段分布在400~750 nm區間內, 最高值出現在688 nm處。 與原始光譜相比, 在492~510, 542~543, 642~652, 657~670和682~692 nm區間內的SPAD相關性有所提高, 且分別達到-0.655, -0.662, -0.697, 0.709和-0.749, 也說明了這些波段的光譜反射率與棗樹葉片SPAD值相關性好, 適合用于敏感波段的挑選。 綜上所述, 紅棗樹葉片反射率光譜做一階導數處理后與SPAD的相關性有較顯著的提高。
結合圖3, 在原始光譜570~620 nm范圍內選擇了相關性高的581, 590, 595和602 nm波段, 690~700 nm波段范圍內選擇695和696 nm共6個特征波段進行CP統計量的計算; 基于光譜一階導數與SPAD值相關性高低, 在492~510, 542~543, 642~652, 657~670和682~692 nm共5個波段內分別選擇相關性達到區間內最高的495, 543, 649, 664和688 nm共5個特征波段計算出其不同組合統計量, 表1為波段的相關系數表。

圖3 SPAD值與光譜反射率之間的相關性

表1 波段的相關系數表
表2為CP統計量計算結果表, 考慮到所有變量組合方式數目較大, 且大部分組合方式的CP統計量都遠高于表2中的幾種組合方式, 只列出CP統計量值靠前的組合。CP統計量越低, 代表該種變量選擇方式重要性程度越高。 且由表2可知, 原始光譜選擇在570~620和690~700 nm范圍內分別選擇595和696 nm時CP統計量絕對值最低, 因此將595與696 nm原始光譜重要程度最高的兩個變量作為建模的特征波段。 光譜一階導數變換后688 nm波段CP統計量絕對值最低, 因此光譜一階導數的特征波段定為688 nm。 原始光譜特征波段696 nm和光譜一階導數特征波段688 nm都處于紅邊波段[13], 說明紅邊與植被的各種理化參數是緊密相關的, 是描述植物色素狀態和健康狀況的重要的指示波段。

表2 特征波段組合及CP統計量計算結果
不難發現, 對于同一個區間的波段組合總有單波段的CP統計量低于多波段組合的CP統計量。 說明相近波段組合建模會使得誤差增大, 這可能是相近波段之間較強共線性造成的, 故每個敏感波段區間只選取一個波段進行建模是合理的。
運用CP統計量選出的3個特征波段以及實測葉片SPAD值建立MLR, SVR和GWLS-SVR模型。 相比較而言, GWLS-SVR主要的優勢是變量系數隨著地理位置而變化, 具有較強的靈活性。 為了明確紅棗樹葉片SPAD值分布是否與地理位置有關, 對其進行Moran’s Ⅰ的計算結果為0.125 8(p<0.1), 呈空間正相關, 說明棗樹葉片SPAD值的分布有顯著的空間聚集性, 適合運用GWLS-SVR模型來建模。
原始光譜(RR)與光譜一階導數(FD)分別基于MLR, SVR以及GWLS-SVR擬合的MSE與R2如圖4所示。 從建模效果來看, 基于原始光譜建立的三種模型中, MLR與SVR的R2低于0.8, MSE也較高, 說明這兩種模型的穩定性較差, 預測效果不理想; GWLS-SVR的R2為0.915, MSE低至3.679, 表明GWLS-SVR的穩定性及估算能力優于MLR與SVR模型。 光譜一階導數變換后的三種模型精度較原始光譜均有所提升, 且MSE整體上都有所降低, 表明數據變換后模型的穩定性和精度有了一定的提高; 而GWLS-SVR在光譜一階導變換后均顯著優于其余兩種模型, 模型的R2提高到了0.975, MLR與SVR的MSE均比GWLS-SVR高約20倍, 綜合上述可得GWLS-SVR模型不僅擬合精度高, 其估計偏差與方差綜合看來均低于其余兩個模型。

圖4 MLR, SVR和GWLS-SVR對實測值與預測值間的擬合圖
從擬合效果來看, GWLS-SVR在原始光譜與光譜一階導數的擬合曲線比起其他兩種模型真實值與預測值均勻分布在1∶1直線周圍, 表明GWLS-SVR的擬合效果較好, 且在光譜一階導數變換后的擬合效果更佳。
為了檢驗三種模型的擬合效果差異, 對于原始樣本利用Bootstrap再抽樣方法進行100次有放回隨機抽樣, 每次抽取67個樣本。 之后對于隨機生成的樣本利用上述三種模型分別建模計算100組MSE的均值(mean of MSE)和方差(varionce of MSE), 并利用配對t檢驗進一步比較原始光譜和光譜一階導數基于GWLS-SVR與其他兩個模型的MSE之間的差異的顯著性。
100次Bootstrap再抽樣并基于三種模型建模后的MSE箱線圖如圖5所示。
由表3及圖5可以看出, 整體上GWLS-SVR的100組MSE在原始光譜與光譜一階導變換后均為最低, 且波動也比較小, 說明100次Bootstrap再抽樣后GWLS-SVR相比于傳統的MLR及SVR模型預測的精度較高且發揮比較穩定, 其中基于光譜一階導數建立的GWLS-SVR模型的MSE最小, 且波動最小, 說明基于光譜一階導數建立的GWLS-SVR模型的模型預測精度最佳且穩定。

表3 Bootstrap再抽樣結果

圖5 100次建模的MSE箱線圖
為了進一步檢驗GWLS-SVR的MSE是否顯著小于其他兩個模型的MSE, 以下分別做GWLS-SVR與其他兩個模型的單邊配對t檢驗。 所得T統計量與p值如表4所示。

表4 t檢驗結果
由表4可見, 4組單邊配對t檢驗的t統計量的絕對值都比較大, 且p值均非常接近0。 所以GWLS-SVR預測的MSE小于其他兩個模型的MSE這一假設在統計學上是高度顯著的。
利用野外實測67個樣點, 219條紅棗樹葉片高光譜數據和棗樹葉片SPAD值數據, 對SPAD值與高光譜波段進行相關性分析、CP統計量特征波段選取、 建立基于特征波段的SPAD值估算模型, 結果表明:
(1)光譜一階導數起到了對原始光譜數據的去噪、 突出高光譜信息的作用, 尤其是在492~510, 542~543, 642~652, 657~670和682~692 nm區間內明顯提高了與SPAD值的相關性。
(2)根據統計量計算發現: 對于同一個敏感波段區間的波段組合總有單個波段的CP統計量低于多個波段組合的CP統計量, 臨近分布的波段之間的存在的較強共線性可能導致這些波段的組合誤差增大。
(3)基于實地采樣數據進行地統計分析若羌縣棗樹SPAD值與地理位置的關聯性, 發現若羌縣存在空間聚集性, 全局莫蘭指數為0.125 8(p<0.1), 表明地理加權最小二乘支持向量機方法適用于估算若羌縣棗樹葉片SPAD值。
(4)基于光譜一階導數的特征波段建立的GWLS-SVR模型的估算能力(R2為0.975, MSE為1.082)優于基于原始光譜特征波段建立的GWLS-SVR模型(R2為0.915, MSE為3.679), 且由結合Bootstrap再抽樣方法與t檢驗的結果來看, 基于光譜一階導數的加入地理位置信息的GWLS-SVR模型為最優的棗樹葉片SPAD值估算模型, 能夠為快速無損的監測紅棗樹生長狀況提供參考。
致謝:感謝若羌縣委的一貫支持, 感謝若羌縣委辦、 縣農業農村局、 自然資源局、 交通運輸局對本項野外調查工作的具體幫助。 感謝縣委辦戶亮亮同志對本項工作的協調和幫助。