基于高光譜成像的煙田土壤pH估測

張 恒,梁太波,宋效東,江 鴻,郭文孟,戴華鑫,翟 振,馮長春,張艷玲

(1. 中國煙草總公司鄭州煙草研究院,鄭州 450001; 2. 中國科學院南京土壤研究所,南京 210008; 3. 四川省煙草科學研究所,成都 610041)

【研究意義】土壤pH是評估土壤質量的重要指標[1]。土壤的物理、化學和生物特性均會受到土壤pH的影響[2],土壤pH不僅影響土壤礦物質和其他養分的生物有效性,也影響微生物活動,改變土壤有機質分解速度,進而影響土壤養分釋放。因此,土壤pH的準確測定對煙田土壤質量監測和土壤改良至關重要。傳統的土壤pH實驗室測定方法費時費力,且易受測定條件的影響,難以滿足智慧農業對田間管理的實時性需求[3]。高光譜遙感技術的發展為土壤成分快速檢測提供了新方法,高光譜技術因其快速、無損且可獲得被檢測物體較高精度的信息,目前已被廣泛用于定量估測土壤中的水分含量[4]、養分狀況[5-6]、有機質(碳)含量[7-8]以及土壤重金屬[9-10]等土壤屬性信息。因此可通過構建土壤pH高光譜定量估測模型實現土壤pH的準確快速測定,為煙田土壤精準施肥提供技術支撐?！厩叭搜芯窟M展】近年來,國內外學者在構建土壤pH和高光譜反射率定量估測模型方面進行了一些研究。如蔡海輝等[11]以棉田土壤為研究區,分析9種單一預處方法處理后的光譜數據與pH的相關性,發現光譜全波段經二階微分處理后結合隨機森林算法建模是估測棉田土壤pH的最佳模型;沈從旺等[12]以重慶涪陵區土壤為研究對象,利用相關分析法篩選特征波段構建pH估測模型,發現支持向量機是估測水稻土pH的最佳建模方法,偏最小二乘回歸是估測紫色土pH的最佳建模方法;Ji等[13]以浙江省水稻田土壤為研究對象,結果發現原始光譜全波段經Savitzky-Golay平滑處理后,使用偏最小二乘回歸建立的土壤pH估測模型效果較好;Yang等[1]以長江中下游平原水稻土為研究對象,在4種建模方法下,對比了使用光譜全波段和遺傳算法篩選特征波段對構建土壤pH估測模型的影響。前人選擇的研究區域多集中于單個縣(區),研究對象多局限于單一土壤類型,所取土壤樣品類型和空間分布差異較小,在更大區域尺度上運用高光譜技術對土壤pH進行建模估測的研究較少;且較多研究在估測不同地區土壤pH時,多使用光譜全波段構建土壤pH估測模型,關于主成分分析對光譜數據降維以及競爭自適應重加權采樣(Competitive adaptive reweighted sampling, CARS)算法在篩選土壤pH特征波段方面的研究鮮有報道。此外,由于高光譜成像技術可獲得被掃描樣品區域內所有像素點的光譜數據[14],且能同時掃描分析多個樣品,所以高光譜成像技術比使用地物光譜儀基于點獲取的高光譜數據更快更準確,然而到目前為止,利用高光譜成像技術估測煙田土壤pH的研究鮮見報道?！颈狙芯壳腥朦c】以我國四川省煙田土壤為研究對象,利用高光譜成像技術獲取土壤樣品的高光譜數據,采用12種光譜預處理方法對原始光譜進行處理,結合原始光譜全波段、主成分分析降維得到的特征光譜以及CARS篩選的特征波段,運用4種建模方法分別構建土壤pH估測模型?！緮M解決的關鍵問題】旨在探索區域尺度下運用高光譜成像技術估測土壤pH的可行性,為四川省煙田土壤pH快速測定提供理論依據和方法參考。

1 材料與方法

1.1 土壤樣品采集與處理

2021年在四川省4個地區(涼山州、攀枝花、宜賓、瀘州)18個縣(區、縣級市)共采集土壤樣品296份,采樣點分布見圖1。具體采集方法:采用5點取樣法采集煙田0～20 cm耕層土壤樣品,四分法保留2 kg,在實驗室土壤樣品風干箱內風干,剔除雜物,研磨后過2 mm篩,分成兩部分,一部分用于土壤高光譜數據的采集,另一部分用于測定土壤pH。用于采集高光譜數據的土壤樣品,在采集前平攤于通風條件良好的實驗室內72 h,以降低土壤水分對高光譜數據采集過程中的不利影響。土壤pH測定參照《土壤農業化學分析方法》[15],采用水浸提電位法測定pH,水土比為2.5∶1.0。

圖1 采樣點分布Fig.1 Distribution of sampling points

1.2 室內高光譜圖像采集與數據提取

高光譜圖像采集裝置主要包括消色差鏡頭(HSIA-OLE23)、可見—近紅外高光譜相機(GaiaField-V10E-AZ4)、電動載物臺、計算機和Spec View圖像采集軟件等。光譜采集方法:將處理好的土壤樣品放置于直徑10 cm、深2.0 cm被黑布包裹的培養皿中,刮平土壤表面后放置在電動載物臺上,通過操作電腦中的Spec View圖像采集軟件,使電動載物臺以1.35 cm/s的速度帶動土壤樣本進入圖像采集箱內,在暗箱環境下完成對土壤樣品高光譜圖像的采集。光譜數據采集范圍為390～2561 nm。在390～1030 nm范圍內,光譜分辨率為2.6 nm,采集250個波段;在967～2561 nm范圍內,光譜分辨率為5.4 nm,采集288個波段。高光譜圖像采集后,需要對采集的高光譜圖像進行黑白版校正,以消除儀器本身在采集過程中帶來的光譜噪聲。由公式(1)對原始高光譜圖像進行校正。

(1)

式中,I0為采集的原始高光譜圖像,W為采集反射率為99%的標準校正白板獲取白板標定圖像,B為采集反射率為0%的內置黑板獲取黑板標定圖像,I為校正后原始高光譜圖像。

對于每個土壤樣本校正后的高光譜圖像,選擇培養皿內土壤樣本區域作為感興趣區域(Region of interest,ROI)。提取土壤樣本ROI內所有像素點的光譜數據,然后對ROI內所有像素點的光譜反射率值進行平均,以此生成平均光譜作為此樣品的原始高光譜反射率數據。

1.3 高光譜數據的預處理

由于每段光譜數據的頭部和尾部均存在噪聲,因此每個土壤樣本只使用408～1007和1012～2500 nm的數據作為進一步分析和建模的高光譜數據。采用最大最小歸一化(Max-min scaling, MMS)、多元散射校正(Multiplicative scatter correction, MSC)和標準正態分布(Standard normal variate, SNV)3種單一預處理方法對原始光譜(Reflectance,R)進行預處理[16]。同時為比較不同組合光譜預處理方法對建模精度的影響,在一階導數(First derivative, D1)、Savitzky-Golay平滑(Savitzky-Golay smooth, SG)和趨勢校正(DT)基礎上,對光譜數據再分別進行MMS、SNV和MSC預處理。

1.4 高光譜數據特征波段提取與降維

為減少高光譜數據的冗余,降低模型的復雜性和計算量。研究中使用CARS算法篩選特征波段和主成分分析(Principal component analysis,PCA)對光譜數據進行降維處理。PCA主要通過某種線性投影將高維數據映射到低維的空間中,可將眾多具有相關性的指標重新組合一組新的無相關的綜合指標來代替,以達到對原始特征進行降維的目的[17]。CARS是通過適應重加權采樣和指數衰減函數篩選出偏最小二乘模型中回歸系數絕對值大的波段,去掉權重較小的波段,利用交互驗證選出交互驗證均方根誤差(Root mean square error of cross validation,RMSECV)最小時的波段子集,從而有效尋找出最優變量組合。研究中CARS算法設置采樣次數為50次[18]。

1.5 建模方法

采用偏最小二乘回歸、嶺回歸、核嶺回歸和支持向量機4種方法構建土壤pH估測模型。偏最小二乘回歸(Partial least square regression,PLSR)是一種廣泛用于土壤高光譜定量分析的線性回歸模型[19],可有效降低自變量間多重共線性問題,避免模型過度擬合以及優化模型估測的精度。嶺回歸(Ridge regression,RR)是常見包含懲罰項的線性學習算法,可減少線性回歸模型的方差[20]。核嶺回歸(Kernel ridge regression,KRR)是嶺回歸(線性最小二乘與L2范數正則化)與核函數結合起來的算法,不同的核函數會導致不同的空間函數,研究中核嶺回歸使用高斯徑向基(RBF)為核函數進行建模。支持向量機(Support vector machine,SVM)是基于統計學習理論的機器學習算法,通過核函數將數據轉換到高維特征空間中去,已被用于解決線性和非線性問題,研究中SVM使用高斯徑向基(RBF)核函數[21]。

建模過程中采用系統網格搜索方法對模型參數進行優化選擇,通過10折交叉驗證,以建模集均方根誤差(Root mean square error of training set,RMSET)最小時確定最佳估測模型。

1.6 模型評估

以決定系數(Coefficient of determination,R2)、建模集均方根誤差(RMSET)、驗證集均方根誤差(Root mean square error of validation,RMSEV)、相對分析誤差(Residual prediction deviation,RPD)等參數對建立的pH估測模型性能進行評估[22]。模型的均方根誤差越小,R2越接近1,表明模型的準確性越高。當RPD≥2時,表明模型可以較準確的對土壤pH進行估測;當1.4≤RPD<2時,表明模型可以粗略對土壤pH進行估測。模型評估參數具體計算公式參考文獻[22]。

1.7 數據分析

高光譜數據的采集和校正均使用 Spec View軟件完成;感興趣區選取、光譜數據提取、光譜數據預處理、特征波段篩選以及模型的建立均使用Python3.9軟件;使用Microsoft Excel 2016軟件制圖。

2 結果與分析

2.1 土壤pH統計特征

為保證建模集和驗證集樣本中pH分布的均勻性,采用梯度質量法[16]對土壤樣本集進行劃分,首先將所有樣本按照土壤pH進行升序排序,然后從小到大每隔3個樣品取1個作為驗證集樣本,把所有樣本以3∶1的比例劃分為建模集和驗證集,如表1所示。土壤總樣本中pH變化范圍為4.360～8.715,均值為6.299。建模集和驗證集的描述統計特征相近,與總樣本各統計特征基本在同一水平,且驗證集樣本pH分布范圍在建模集樣本pH分布范圍之內,說明驗證集樣本在建模集中分布均勻。

表1 土壤樣本pH描述統計

2.2 不同土壤pH高光譜反射率特征

將總樣本按照土壤pH大小升序排列,均分為6組后求每組對應pH和光譜數據的平均值(圖2)。不同pH的土壤光譜反射率曲線雖各有不同,但變化趨勢較為相似。光譜反射率在400～800 nm可見光范圍內迅速上升且光譜反射率值相對較小,在800～2500 nm范圍內光譜反射率曲線波動上升,在1400、1900和2200 nm處均存在明顯的水分吸收谷。此外,在600～2500 nm范圍內,不同pH的光譜反射率曲線差異明顯,土壤pH越高,土壤光譜反射率越低,這種差異的存在是后續運用高光譜建立土壤pH估測模型的基礎。

圖2 不同pH土壤原始光譜反射率曲線Fig.2 Original spectral reflectance curves of soil with different pH

2.3 基于光譜全波段的土壤pH估測模型構建

由表2～3可知,與使用原始光譜全波段建模相比,使用單一預處理或組合預處理方法處理后的光譜全波段,結合不同建模方法建立的土壤pH估測模型性能均有不同程度提升。其中使用PLSR方法建模時,基于DT-MMS預處理方法使用全波段建立的土壤pH估測模型精度最高;分別使用SVM、RR以及KRR 3種建模方法時,基于D1-SNV預處理方法使用光譜全波段建立的pH估測模型精度均最高。整體看,不同建模方法得到的土壤pH最優估測模型,均使用光譜組合預處理方法處理后的光譜全波段,說明在相同建模方法下,選擇合適的光譜組合預處理方法對原始光譜進行處理,有助于提高土壤pH估測模型的準確度。原始光譜經D1-SNV預處理方法處理后,結合4種建模方法建立的pH估測模型精度均較高,因此在后續分析中選擇D1-SNV對原始光譜進行預處理,以降低數據噪聲,提高光譜數據有效性。

表2 基于光譜全波段的PLSR和SVM模型構建及驗證

表3 基于光譜全波段的RR和KRR模型構建及驗證

2.4 基于CARS篩選特征波段和PCA降維的pH估測模型構建

2.4.1 CARS算法篩選特征波段 利用CARS算法篩選D1-SNV預處理后的特征波段。由圖3-a可知,隨著運算次數增加CARS算法篩選出的特征波段變量個數逐漸減少。經過前10次的運算,篩選的特征波段變量個數由504個降低至182,呈明顯降低趨勢,隨后特征變量個數呈平穩降低趨勢。由圖3-b可知,在0～50次運算過程中RMSECV波動變化明顯,在0～16運算次數過程中RMSECV整體呈降低趨勢,在16～50運算次數過程中RMSECV呈明顯升高趨勢,說明在16～50運算次數過程中與土壤pH無關的噪聲信息增多。當運算次數為第16次時,RMSECV最小,篩選的特征光譜波段組合最優。結合圖3-a可知當運算次數為16時,對應的特征波段變量數目為93,其中有35個波段分布在可見光400～760 nm范圍內,有58個波段分布在近紅外800～2500 nm。

2.4.2 確定PCA最佳維度數 圖4顯示基于D1-SNV預處理方法,不同建模方法下模型評估參數RPD隨維度數(1～222)增長的變化規律?？梢娽槍Σ煌７椒?運用PCA降維過程中得到的最佳維度數差異明顯。當維度數為60時,估測土壤pH的PLSR模型驗證集RPD達最大值;當維度數是180時,土壤pH的SVM估測模型驗證集PRD達最大值;當維度數為最大值176時,估測土壤pH的RR模型驗證集PRD達最大值;當維度數為最大值222時,估測土壤pH的KRR模型驗證集PRD達最大值。因此選擇60、180、176和222分別作為PLSR、SVM、RR和KRR 4種建模方法在后續構建土壤pH估測模型過程中PCA降維的最佳維度數。

2.4.3 土壤pH估測模型的構建 由表4可知,使用PLSR和RR方法建模時,運用CARS篩選特征波段建立的pH估測模型精度均較高,PLSR模型驗證集R2、RMSEV和PRD分別為0.758、0.555和2.034;RR模型驗證集R2、RMSEV和PRD分別為0.751、0.563和2.002;使用KRR和SVM方法建模時,運用全波段和PCA降維建立的土壤pH估測模型精度均較高,模型驗證集的R2均大于0.750,RMSEV均小于0.570,RPD均大于2.0。

整體看,使用全波段建模情況下,土壤pH估測模型精度排序為KRR>SVM>RR>PLSR,采用KRR模型估測效果最佳(RPD=2.105)。經過CARS篩選特征波段后,pH估測模型精度排序為PLSR>RR>KRR>SVM,采用CARS-PLSR模型估測效果最佳(RPD=2.034),RPD較全波段下的PLSR模型提升了0.29。經過PCA對光譜全波段降維后,pH估測模型精度排序為SVM>KRR>RR>PLSR,采用PCA-SVM模型估測效果最佳(RPD=2.040)?？梢娫谌ǘ?、PCA降維和CARS篩選特征波段建模情況下,得到的pH最優估測模型精度相近。但由于CARS-PLSR模型入選變量數目較少、建模方法解釋性較強和模型精度較高,故在D1-SNV預處理方法下,結合CARS篩選特征波段使用PLSR方法建立的土壤pH估測模型(D1-SNV-CARS-PLSR)效果最佳。

3 討 論

3.1 不同光譜預處理方法對土壤pH建模的影響

構建估測模型前對原始光譜進行適當的預處理,可明顯降低外界環境(溫度、濕度、光線等)和土壤自身(質地、顆粒大小以及土壤類型等)對采集光譜數據產生的干擾,有助于提高模型的估測精度[3]。本研究發現,原始光譜經單一預處理或組合預處理后,建立的土壤pH估測模型精度均有所提升,且不同建模方法得到的估測土壤pH最佳模型,均使用的是光譜組合預處理方法。這與Lin等[23]研究結果一致,其研究表明基于組合預處理方法建立的土壤pH估測模型精度較高。同時研究還發現基于DT-MMS預處理后光譜全波段建立的土壤pH PLSR估測模型均方根誤差較小且精度較高,基于D1-SNV預處理后光譜全波段建立的土壤pH SVM、RR和KRR估測模型精度均較高。原因可能是:趨勢校正(DT)有效校正光譜數據中的傾斜和曲線漂移[24];最大最小歸一化(MMS)有效消除了尺度差異帶來的不良影響[25];一階導數(D1)對重疊光譜進行拆分,增大了光譜曲線中波谷和波峰的特征[22];標準正態分布(SNV)有效降低了因顆粒分布不均或顆粒大小不同產生的散射影響[26];本研究將上述預處理方法組合起來使用,有效消除了原始光譜曲線中的曲線漂移和散射效應,增大了光譜數據的信噪比,從而顯著提高了模型估測土壤pH的準確度和穩定性。

圖3 CARS算法篩選特征波段Fig.3 Feature bands screened by CARS algorithm

n為PCA降維過程中的維度數。n is the number of dimensions in PCA dimension reduction.

表4 不同變量篩選方法的土壤pH估測模型構建及驗證

3.2 CARS篩選的土壤pH特征波段

CARS算法是一種有效的特征波段篩選方法,主要通過逐步去除不重要變量和冗余,從全波段中篩選出用于建模的最優特征波段變量集合。本研究中,與全光譜相比,CARS算法篩選出了93個土壤pH特征波段,占光譜全波段數目的18%左右,其中篩選的特征波段在可見光和近紅外區域內均有分布,主要集中在近紅外波段800～2500 nm。與Wang等[9]研究結果相似,其研究發現CARS算法篩選土壤pH的特征波段主要在1903 nm以下的可見光和近紅外區域;與Vohland等[27]發現CARS算法篩選土壤pH特征波段只限于1915 nm以外的結果相悖。究其原因可能是由于光譜活性土壤成分(如有機質和粘土)的協變性在不同區域土壤中存在顯著差異,導致在不同地區使用同一特征波段篩選方法篩選出的pH特征波段存在差異。

3.3 PCA光譜數據降維對土壤pH建模的影響

主成分分析是最常用的數據降維方法,可以保證在信息量盡量不丟失的情況下對高維高光譜數據進行降維處理,以減少建模過程中的計算量和建模時間。研究中發現與使用光譜全波段建模相比,使用PCA對光譜降維后結合PLSR和RR建模方法建立的土壤pH估測模型性能有所提升,但提升程度并不顯著。應該注意到運用PCA對光譜降維建模雖然模型性能提升不明顯,但可以有效降低建模過程中的運算量,提高模型的擬合效果[12]。此外,使用SVM和KRR 2種建模方法時,基于光譜全波段建立的土壤pH估測模型較PCA 降維建模取得了更好的估測性能,造成這一結果的因素可能是由于PCA降維過程中雖然減少了用于建模的數據量,但也損失了部分重要光譜信息,導致所建土壤pH估測模型效果不理想[28],這一點在本研究中也得到了證實,研究中基于PCA篩選特征波段運用KRR方法建模時,模型評估參數RPD隨維度數增長而增加,但由于降維過程中最大維度數限制導致模型精度低于使用光譜全波段建立的KRR估測模型。

3.4 不同建模方法對土壤pH建模的影響

在4種建模方法中,無論使用光譜全波段還是PCA降維后的特征光譜,KRR和SVM 2種建模方法都取得了比PLSR和RR更好的估測準確度。這可能是由于土壤pH與有機質、粘土等光譜活性成分之間關系較密切,且主要存在內在非線性關系[2],而偏最小二乘回歸和嶺回歸屬于線性模型,在解決非線性關系時存在一些不足,故導致模型精度相對偏低。但使用CARS算法篩選的特征波段建模時,PLSR和RR 2種建模方法較KRR和SVM取得了更好的估測精度,可能是由于CARS算法主要通過偏最小二乘回歸模型從光譜全波段數據中篩選出光譜波段變量組合,故CARS篩選的特征波段與土壤pH之間主要存在線性關系。另外,通過對比分析發現基于D1-SNV預處理方法,使用CARS篩選特征波段建立的土壤pH PLSR估測模型精度較高且計算量較小,模型驗證集R2、RMSEV和RPD分別為0.758、0.555和2.034,說明在區域尺度上運用高光譜成像技術估測土壤pH具有可行性。但與蔡海輝等[11]在阿拉爾市對棉田土壤pH建立的最優模型相比,本研究建立的最佳估測模型精度相對偏低,原因可能是由于本研究所選研究區域尺度較大,所取土壤樣品理化性質和土壤類型復雜多樣,故導致建立的pH最佳估測模型精度偏低,但也應該注意到本研究所建立的土壤pH估測模型適用區域較大。

因此,在未來運用高光譜成像技術對土壤化學成分進行估測研究時,需要充分考慮土壤樣本自身的光譜特性和土壤類型,深入分析土壤化學成分與光譜數據間的內在關系,加強對光譜數據預處理方法、特征波段篩選方法以及機器學習建模算法方面的研究[23],以提高高光譜估測土壤pH的準確度。此外,雖然研究中建立的土壤pH估測模型精度相對不太高,但所建模型的適用區域較廣,模型的泛化能力較強,可以有效地對四川省植煙土壤pH進行估測,是否對于其它省份或區域的植煙土壤具有普適性,還需要采集更多具有代表性的植煙土壤樣品對模型進行完善。

4 結 論

(1)研究區域耕層土壤pH變化范圍為4.360～8.715,均值為6.299,變異系數為17.998%,土壤pH呈中等變異性。

(2)使用單一預處理或組合預處理方法處理后的光譜波段建模,有利于提高土壤pH估測模型的精度;其中原始光譜經D1-SNV預處理方法處理后,結合4種建模方法建立的pH估測模型精度均較高。

(3)CARS算法篩選出了93個土壤pH特征波段,篩選的特征波段較均勻的分布在可見光和近紅外區域內,且多數集中在近紅外波段800～2500 nm。原始光譜經D1-SNV預處理后,基于偏最小二乘回歸和嶺回歸兩種線性建模方法建模時,使用PCA降維或光譜全波段建立的土壤pH估測模型精度均低于CARS篩選特征波段建立的估測模型;基于核嶺回歸和支持向量機兩種非線性建模方法建模時,使用PCA降維或光譜全波段建立的土壤pH估測模型精度均高于CARS篩選特征波段建立的估測模型。綜合比較下D1-SNV-CARS-PLSR模型估測土壤pH的效果最好,模型建模集與驗證集的R2分別為0.802和0.758,RMSET和RMSEV分別為0.504和0.555,驗證集RPD為2.034。