基于植被指數(shù)的獼猴桃根域土壤水分反演影響因素研究

張 軍 鄧俊濤 倪國威 牛子杰 潘時佳 韓文霆

(1.西北農(nóng)林科技大學機械與電子工程學院， 陜西楊凌 712100；2.西北農(nóng)林科技大學中國旱區(qū)節(jié)水農(nóng)業(yè)研究院， 陜西楊凌 712100)

0 引言

水資源是農(nóng)業(yè)生產(chǎn)中不可或缺的自然資源。植物水分脅迫狀況反映了其水分需求狀況，從而影響灌溉決策，但現(xiàn)有接觸式測量方法難以大面積準確監(jiān)測植株個體的水分脅迫狀況。針對此問題，現(xiàn)階段對于植物水分參數(shù)的部分研究已經(jīng)確立了其與光譜、植被指數(shù)和溫度的相關性，并建立了較為準確的反演模型。

植物表型監(jiān)測的研究表明植物葉片形態(tài)具有判斷植物脅迫狀態(tài)的潛力，這種形態(tài)包括植物葉片多光譜反射率、葉片表面顏色、葉片熱紅外信息。基于植物葉片形態(tài)的無人機、衛(wèi)星遙感作物監(jiān)測模型也證明了其在遙感反演領域的應用潛力。在對監(jiān)測模型的研究中，人工神經(jīng)網(wǎng)絡(Artificial neural network，ANN)算法的引入更進一步提升了遙感反演植物脅迫狀態(tài)的能力[1-9]。

植物水分參數(shù)反演模型中植物的生理參數(shù)野外測量易受測量環(huán)境、測量方式影響，實驗室測量方式耗費時間，但是對植物水分參數(shù)及其監(jiān)測手段的研究為對土壤水分的反演提供了理論基礎。研究證明在衛(wèi)星遙感尺度上，土壤不同波長反射率也可以反映土壤干旱情況，這種反演方式也為低空無人機遙感反演土壤干旱情況提供了參考。相關研究已證明植物冠層反射率與土壤干旱程度具有一定關系，并建立了基于植被指數(shù)的土壤含水率反演模型，在此基礎上有相關學者提出通過圖像卷積進一步提取植物冠層反射率特征從而提升土壤含水率反演精度，但目前尚未對模型精度與影響因素進行研究[10-14]。

土壤含水率反演模型可以分析遙感圖像上植被指數(shù)的空間異質(zhì)性，對農(nóng)田土壤水分狀況進行評估，但是現(xiàn)有模型與研究依然面臨以下問題：①果樹相比于大田作物體積更大，根系與冠層分布更廣。②冠層中存在不同年齡的葉片，這些葉片在反射率上有著較大的差異使得單個植株冠層植被指數(shù)具有空間異質(zhì)性。③植株冠層之間會發(fā)生相互重疊，這種情況在非規(guī)范化果園中尤為明顯，這使得在進行土壤含水率反演時不能像大田作物連續(xù)采集區(qū)域反射率，在評估植株根域土壤含水率(Root domain soil water content，RSWC)時需要討論其數(shù)據(jù)采集的位置與大小。

針對以上問題，本文對獼猴桃冠層植被指數(shù)與RSWC進行相關性和顯著性分析，并根據(jù)分析結果建立不同采樣尺寸的植被指數(shù)-RSWC反演模型，最后通過繪制各個模型精度與采樣尺寸的回歸曲線研究采樣尺寸對反演精度的影響并分析原因。

1 材料與方法

1.1 試驗區(qū)域與試驗設計

于2021年8月在陜西省寶雞市眉縣西北農(nóng)林科技大學獼猴桃實驗站(107°59′31″N，34°07′28″E，海拔643.22 m)進行試驗，該地區(qū)屬暖溫帶大陸性半濕潤氣候，年平均氣溫12.9℃，平均降水量609.5 mm，平均日照時數(shù)2 015.2 h，無霜期218 d。每年3—5月，回暖較快，秋季受冷空氣影響，晝夜溫差較為明顯，是關中地區(qū)秋雨最多的區(qū)域之一；較為適合獼猴桃生長，但仍需人工對果園進行水分管理。該實驗站種植有多種品種獼猴桃，本文試驗植株品種為徐香，正處于果實膨大期(5—9月)，水分需求量大，樹齡5～6年，株距3 m，行距4 m。試驗地位置及試驗區(qū)域標注如圖1所示。

圖1 試驗地信息Fig.1 Experimental site information

試驗設置了3組不同梯度(每個梯度包含4個植株)的試驗組，統(tǒng)計不同處理下的土壤含水率空間變異性；并記錄灌溉后5 d的數(shù)據(jù)，統(tǒng)計其時間變異性。3組試驗組僅對土壤含水率進行控制，對應田間持水量的100%(充分灌溉)、90%(輕度虧缺)、80%(中度虧缺)，分別表示為T1、T2、T3，達到對應土壤含水率后使其土壤含水率自然下降，并在下降期間采集其土壤含水率與多光譜數(shù)據(jù)。

1.2 數(shù)據(jù)采集

試驗數(shù)據(jù)采集于2021年8月3、4、6、9、10日的12：00—14：00，該時段獼猴桃果樹較易出現(xiàn)水分脅迫現(xiàn)象且無人機獲取的正射影像陰影較少。

遙感數(shù)據(jù)與地面數(shù)據(jù)采集同時進行，使用大疆Phantom 4 V2.0型無人機搭載Survy3型多光譜相機(MAPIR，美國)在高度15 m采集試驗區(qū)域不同波長反射率的像元亮度(Digital number，DN)值。獲取的圖像使用MAPIR多光譜校準軟件(MAPIR camera controler，MCC)與多光譜標定板校正；圖像地面分辨率為0.007 m，根據(jù)果園的株行距試驗區(qū)域按照600像素×400像素(約為4 m×3 m)進行分割，并作為初始采樣尺寸對區(qū)域內(nèi)的DN值進行計算，如圖2所示。

圖2 采集設備與采集點Fig.2 Acquisition equipment and acquisition points

為了獲得正射影像，無人機在采集試驗區(qū)圖像時會向外拓展一部分距離，保證Pix4D Mapper軟件在進行正射拼接時獲得足夠的重疊率。

土壤含水率采集深度為40 cm，處于水分活躍層(20～80 cm)，受到大氣影響較小，但依然活躍，是獼猴桃根系的主要生長深度[15]。灌溉后在獼猴桃主干附近采集40 cm處的土壤剖面，并使用PR-3001-TRREC-N01型土壤水分檢測儀(仁智測控技術公司)重復測量3次取平均值作為被測植株RSWC，并對異常值進行篩選去除。

1.3 數(shù)據(jù)處理

圖像處理首先使用MCC對圖像進行校正；其次使用Pix4D Mapper軟件拼接試驗區(qū)域正射影像；使用Python的PIL庫對圖像進行不同尺寸的中心剪裁，最后提取不同區(qū)域內(nèi)的DN值并計算植被指數(shù)(Vegetation index，VI)，評估的植被指數(shù)包括綠色指數(shù)(GI)、修改型土壤調(diào)整植被指數(shù)(MSAVI)[16]、綠度歸一化植被指數(shù)(gNDVI)[17]、歸一化植被指數(shù)(NDVI)[18]、優(yōu)化型土壤調(diào)節(jié)植被指數(shù)(OSAVI)[19]、重歸一化植被指數(shù)(RDVI)[20]、土壤調(diào)整植被指數(shù)(SAVI)[21]、簡單比值指數(shù)(SR)[22]。

1.4 統(tǒng)計分析

使用SPSS軟件在95%置信水平(P<0.05)下使用單向單因素方差分析。為了比較不同處理的植株之間的植被指數(shù)，將5 d數(shù)據(jù)分為60組，RSWC為25.4%～30.9%，并使用Pearson系數(shù)對植被指數(shù)與RSWC共計9個變量進行相關性分析。

使用不同采樣尺寸下的植被指數(shù)建模，評估采樣尺寸對建模精度的影響。使用多層感知機(Multi-layer perceptron，MLP)算法進行建模，MLP的隱節(jié)點采用輸入向量與權向量的內(nèi)積作為激活函數(shù)的自變量，激活函數(shù)采用ReLU函數(shù)。各參數(shù)對網(wǎng)絡的輸出具有同等地位的影響，因此MLP是對非線性映射的全局逼近，該模型以冠層植被指數(shù)作為輸入，使用MLP模型預測RSWC。

最后，采集試驗植株新生與成熟葉片葉綠素SPAD值，統(tǒng)計分析同一植株的葉綠素含量方差，討論各個參數(shù)對反演精度的影響原因。

2 結果與分析

2.1 植被指數(shù)與根域含水率變化

于2021年8月3日進行灌溉，灌溉后讓植株自然生長。根據(jù)周邊氣象站數(shù)據(jù)記錄，試驗期內(nèi)日平均空氣相對濕度(Mean relative humidity，MRH)與測量時空氣相對濕度(Relative humidity，RH)總體呈上升趨勢；測量時空氣相對濕度最低值出現(xiàn)在8月6日，為46%，試驗期間空氣相對濕度變化如圖3所示。

圖3 2021年8月日平均空氣相對濕度與采集時間點 空氣相對濕度Fig.3 Daily mean relative humidity in August, 2021 and relative humidity at collection time point

根據(jù)記錄，灌溉后的RSWC自然下降，第3天下降至25.4%，為了防止水分虧缺對獼猴桃植株造成不可逆損傷，在8月6日后對獼猴桃果園進行少量連續(xù)灌溉使其逐漸恢復；之后3組試驗組RSWC持續(xù)回升。T3(中度虧缺)組RSWC在整個試驗期間一直較低，最高為27.4%(8月10日)，同日T1、T2組大致為30.0%。T3組RSWC最低為25.4%，同日T1、T2組RSWC也達到試驗期間最低，大致為28.0%。8月6日和10日，T1和T2組RSWC十分接近，這是由于RSWC受到作物自身生理作用、冠層結構以及氣候因素多重影響，這使得灌溉后自然下降的試驗組間可能出現(xiàn)相近的測量數(shù)據(jù)，這也是果園內(nèi)RSWC時間與空間異質(zhì)性的一種體現(xiàn)。同時，也正是復雜因素導致的水分異質(zhì)性，才有必要進行果園內(nèi)水分精細管理。試驗期間不同組RSWC變化如圖4所示。

圖4 2021年8月不同處理試驗區(qū)域平均RSWC日變化Fig.4 Daily variations of average RSWC in different test areas in August, 2021

根據(jù)不同采樣尺寸下獲取的冠層DN值，計算不同的植被指數(shù)。統(tǒng)計5 d內(nèi)3種不同處理的12個植株，共獲取60組RSWC與植被指數(shù)樣本，每組樣本包含8個植被指數(shù)；60組樣本的植被指數(shù)與RSWC變化如圖5所示。從圖5可以看出，除RDVI、SR、GI外，其他的植被指數(shù)均處于0～0.4之間；在RDVI、SR、GI中，SR與GI的變化幅度較小，分別處于1.0～1.5與0.5～1.0區(qū)間；而RDVI則處于1.0～2.1之間，變化跨度大，且不同處理與時間的樣本之間變化幅度較大，這種現(xiàn)象表明RDVI對土壤含水率的變化較為敏感，將有利于反演模型的建立；植被指數(shù)GI與其他植被指數(shù)變化趨勢相反。

圖5 所有樣本RSWC與8個植被指數(shù)變化曲線Fig.5 Changes of RSWC and eight vegetation indices for all samples

2.2 不同植被指數(shù)與RSWC關系

為了研究不同植被指數(shù)對RSWC反演的影響，使用SPSS軟件對正則化后的不同植被指數(shù)與RSWC進行Pearson相關性分析。

不同植被指數(shù)與RSWC的Pearson相關系數(shù)如表1所示，其中只有RDVI的Pearson相關系數(shù)大于0.6，表示其與RSWC呈強相關性；GI、SR與RSWC的相關系數(shù)絕對值小于0.5，其中SR在所有植被指數(shù)中與RSWC相關系數(shù)最小，呈弱相關性；其余植被指數(shù)的相關系數(shù)相互之間差異不明顯。

表1 不同植被指數(shù)與RSWC的Pearson相關系數(shù)Tab.1 Pearson correlation coefficient between different vegetation indices and RSWC

為了排除植被指數(shù)之間的互相關，通過進行兩兩參數(shù)之間的Pearson相關性分析，獲得Pearson相關性矩陣，如圖6所示。橫縱軸交叉處即為Pearson相關系數(shù)。GI與RSWC及其他植被指數(shù)的Pearson相關系數(shù)小于0，呈負相關；其中GI、SR、gNDVI與RDVI的互相關性相較于其他植被指數(shù)較弱(Pearson相關系數(shù)低于0.85)。

圖6 Pearson相關系數(shù)矩陣Fig.6 Pearson correlation coefficient matrix

通過SPSS的單因素方差分析計算植被指數(shù)與RSWC之間的顯著性，結果如表2所示；除SR外其他植被指數(shù)均與RSWC呈現(xiàn)出較強的顯著性(P<0.05)，其中RDVI與RSWC的顯著性最強，其次是gNDVI、MSAVI、GI，這些植被指數(shù)的P值均小于0.04。

表2 植被指數(shù)與RSWC顯著性Tab.2 Significance result of vegetation index and RSWC

根據(jù)相關性與方差分析結果，考慮到植被指數(shù)與RSWC的相關性與顯著性以及植被指數(shù)之間的互相關性，本文使用RDVI作為反演模型的主要植被指數(shù),使用RDVI與RSWC建立基于多層感知機反演模型。

2.3 不同采樣尺寸與RSWC相關性

使用RDVI作為輸入(兩列輸入向量均為RDVI)建立多層感知機反演模型時，發(fā)現(xiàn)不同冠層采樣尺寸獲得的RDVI數(shù)據(jù)組對同一參數(shù)的人工神經(jīng)網(wǎng)絡模型反演精度會造成影響。在每個植株冠層初始采樣尺寸(600像素×400像素)的基礎上對采樣范圍進行變換，如圖7所示，圖中b為采樣區(qū)域?qū)挾龋琱為采樣區(qū)域高度，L為采樣區(qū)域?qū)蔷€長度，黑框為初始采樣尺寸。采樣區(qū)域?qū)挾扔嬎愎綖?/p>

式中i——比例系數(shù)

bi——經(jīng)過比例系數(shù)i裁剪后的采樣區(qū)域?qū)挾炔煌叽缦碌闹脖恢笖?shù)將被用來訓練同一個MLP，以比較不同植被指數(shù)下對反演模型精度的影響，MLP的參數(shù)分別為：層結構為(2，2，1)，優(yōu)化器使用L-BFGS算法，最大訓練世代為2000代，激活函數(shù)為ReLU函數(shù)，學習率使用動態(tài)更新策略。

圖7 采樣尺寸變換過程示意圖Fig.7 Sampling size transformation process

模型訓練集包括所有不同處理與觀測時間數(shù)據(jù)的60%，剩余40%作為測試集。其中訓練集每一組數(shù)據(jù)由兩個植被指數(shù)與對應土壤含水率組成，用于神經(jīng)網(wǎng)絡模型訓練，測試集由具有未知標簽的植被指數(shù)組成，用于評估模型預測效果。

對使用不同采樣尺寸下的植被指數(shù)建立的反演模型決定系數(shù)R2進行統(tǒng)計，獲得不同采樣尺寸與模型精度，如表3所示。其中，采樣面積A=hb。

表3 采樣尺寸與模型精度Tab.3 Sampling size and model accuracy

從R2在L、A上的分布可以看出，模型結果的決定系數(shù)隨著采樣尺寸的減小，先增后減，但是不同采樣尺寸建模下的模型均方根誤差(RMSE)沒有顯著差別，這說明會存在最佳采樣尺寸使得模型精度達到最大。分別建立采樣區(qū)域?qū)蔷€長度與模型決定系數(shù)二次、三次擬合曲線和采樣區(qū)域面積與模型決定系數(shù)二次、三次擬合曲線L-R2和A-R2，如圖8所示，可以看到L-R2與A-R2散點圖的不同階次多項式擬合曲線獲得的最高決定系數(shù)分別為0.993與0.991，說明采樣尺寸對模型結果的決定系數(shù)有很大影響。

圖8 采樣區(qū)域參數(shù)與模型決定系數(shù)二次與三次擬合曲線Fig.8 Quadratic and cubic fitting curves of sampling area parameters and model determination coefficient

在L-R2散點圖中，模型最大決定系數(shù)出現(xiàn)在L為360像素左右，此時b大致為300像素；在A-R2散點圖的二次擬合曲線中，模型最大決定系數(shù)出現(xiàn)在A為50 000像素左右，三次擬合曲線則A為64 000像素左右，此時b為273～309像素。根據(jù)以上分析，b在273～300像素之間時建模精度達到最大，由于地面分辨率為0.007 m，所以這個范圍的寬度在實際中應該在1.953～2.170 m之間，面積在2.540～3.038 m2之間。

3 討論

獼猴桃作為喜水植物對水分的需求量比其他果樹多[23]，最適宜獼猴桃種植地區(qū)的降雨量要達到800～1 200 mm，而陜西關中地區(qū)的降雨量不能滿足這一要求；在自然條件下，獼猴桃植株會逐漸出現(xiàn)水分脅迫，這個現(xiàn)象在果實膨大期尤為明顯。這種由于土壤含水率下降導致的脅迫可以通過葉水勢來進行評估，其中正午葉水勢反映植株水分最大虧缺程度。眾多關于植被指數(shù)與水分參數(shù)相關性的研究[4,6,24-26]以及葉水勢與土壤含水率相關性的研究表明，植物在水分脅迫時的植被指數(shù)、葉水勢以及土壤含水率三者之間存在較強的相關性[27]，這也使得通過植株冠層對土壤含水率進行反演成為了可能；本文將8種植被指數(shù)與土壤含水率進行了相關性與顯著性分析，確定了RDVI為表征RSWC的最佳指標(相關系數(shù)為0.744，P<0.05)，結果與張智韜等[13]的剔除土壤背景的RSWC研究結論類似，通過冠層植被指數(shù)可以較好地反映出一定深度土壤的含水率；SEO等[14]通過卷積神經(jīng)網(wǎng)絡建立冠層圖像與土壤含水率的研究也說明了這一點，但是針對于不同根系深度的植株還有待研究。

獼猴桃作為藤本果樹與大田作物有著較大的區(qū)別，這主要體現(xiàn)在：①單株經(jīng)濟價值更高。②植株體型龐大，根系更深覆蓋面積更廣。③葉片寬大相互交疊，新老葉片的水分狀態(tài)與顏色差異較大。在對獼猴桃冠層的DN值進行采集時，發(fā)現(xiàn)采集區(qū)域不同將對結果造成較大的影響，從而影響植被指數(shù)與建模精度。對獼猴桃冠層的NDVI圖像反歸一化后進行了假彩色處理，如圖9所示，圖中黑色圓圈代表植株的主干位置。

圖9 NDVI反歸一化后的假彩色圖像Fig.9 False color images after inverse normalization of NDVI

從圖9可以發(fā)現(xiàn)，獼猴桃冠層的NDVI分布并不均勻，如圖9冠層頂部的植被指數(shù)普遍大于冠層邊緣的植被指數(shù)，且呈向中心聚集的趨勢；部分冠層并沒有完全覆蓋土壤，且土壤的NDVI與冠層的植被指數(shù)有較大差異；如圖9b果園內(nèi)部獼猴桃植株之間會發(fā)生藤蔓纏繞與冠層重疊，重疊部分NDVI較高，但是無法判斷屬于哪一個植株。

同時在8月10日，使用SPAD502(Konica Minolta Japan)隨機對每個植株的3片葉子的上中下部分進行了3次測量，獲得每個植株3片葉子的平均SPAD值，該測量值表示葉片葉綠素濃度，如 圖10 所示。記錄葉片包括冠層頂部成熟葉片與邊緣新生葉片，統(tǒng)計數(shù)據(jù)表明，不同處理植株之間，葉片SPAD值差距最大為21.67；即使是同樣處理的植株，葉片SPAD值也有很大差異，最大差值也達到21.17，標準差為6.403；同一個植株不同區(qū)域葉片SPAD最大差值為12.26。通過觀察，成熟葉片與新生葉片在顏色，質(zhì)地上也存在較大差異。研究表明植物葉片通過葉綠素和葉黃素等植物色素吸收可見光波段的大部分輻亮度[24]，植物受外界脅迫時會導致葉綠素色素的損害，光合作用效率下降，這也使得葉片吸光度和反射率的變化[28]，改變了自身冠層的反射率格局，使可見光波段反射率增加，近紅外波段反射率降低，從而改變植被指數(shù)。

圖10 所有植株葉片SPAD平均值Fig.10 Mean values of SPAD in leaves of all plants

根據(jù)對假彩色圖像中NDVI的分布與葉片SPAD值的統(tǒng)計，推測以下原因使得采樣尺寸對模型精度造成了影響：①新生與成熟葉片的反射率差別及兩者在冠層上的不均勻分布導致了冠層植被指數(shù)分布不均勻，新生葉片植被指數(shù)較低，而成熟葉片的相對較高，這也使得在變換采樣尺寸時獲得的植被指數(shù)平均值不同；越是向冠層中心集中，植被指數(shù)均值也會更加趨近于成熟葉片均值。②土壤像元植被指數(shù)與冠層像元植被指數(shù)有著較大差異，部分獼猴桃冠層沒有完全遮蓋土地，與同組遮蓋完全的植株相比，提取到的植被指數(shù)均值會更??；越是接近冠層中心，裸露的土地越少，從而也減弱了土壤像元對冠層植被指數(shù)的影響。③獼猴桃植株屬藤本類果樹，果園內(nèi)棚架相連使得獼猴桃藤蔓相互纏繞重疊，在大采樣尺寸下，小冠層植株的植被指數(shù)均值可能包括周圍其他處理植株冠層的植被指數(shù)；采樣范圍接近冠層中心時，將避免其他冠層對目標植株的影響。④當采樣尺寸過小時實際上是對冠層頂部中心的幾片葉子進行植被指數(shù)采集，葉子自身狀態(tài)(例如病蟲害等)將會對植被指數(shù)均值造成極大的影響，這也是在縮小采樣尺寸到一定程度后模型精度下降的原因。

經(jīng)過試驗，果實膨大期徐香獼猴桃的植被指數(shù)采集面積在2.540～3.038 m2之間較為合適，本研究中并沒有考慮采樣形狀與植株品種的影響。

ROMERO等[6]結合機載高光譜獲得的高分辨率圖像與人工神經(jīng)網(wǎng)絡，建立了葡萄園的葉水勢反演模型，該方法取代了傳統(tǒng)的壓力室法測量葉水勢，研究證明了通過無人機高光譜監(jiān)測農(nóng)田含水率的可行性與ANN在遙感數(shù)據(jù)挖掘中的巨大潛力。本文通過MLP神經(jīng)網(wǎng)絡建立了植被指數(shù)-RSWC反演模型，使用2.540～3.038 m2之間采集的RDVI數(shù)據(jù)訓練的模型表明，基于MLP網(wǎng)絡的土壤含水率反演模型可以在所有數(shù)據(jù)擬合上取得較高的精度(R2為0.638，RMSE為0.016)，表明了模型的精確性與泛用性。人工神經(jīng)網(wǎng)絡等算法仍處于發(fā)展階段，隨著算法更替將會使得模型的精度與泛化能力進一步提升，這也將是未來的一個研究方向。

4 結束語

在獼猴桃水分虧缺脅迫時RDVI為指示40 cm處RSWC的最相關與最顯著指標。不同采樣尺寸的模型精度表明，模型反演精度與采樣尺寸有著密切關系，確定了果實膨大期徐香獼猴桃的植被指數(shù)最佳采集面積在2.540～3.038 m2之間。