基于VIS-NIR 光譜的濱海濕地土壤碳氮比預測建模分析

張清文，楊曉蕪，楊 睿，尹 軒，寇財垚

（1.華北理工大學 礦業工程學院，河北 唐山 063210；2.國家知識產權局專利局專利審查協作廣東中心，廣東 廣州 510700）

引言

土壤有機碳（Soil organic carbon，SOC）和土壤全氮（Soil total nitrogen，STN）作為土壤質量和健康的主要指標，在保持土壤質量方面發揮著重要作用[1]。SOC 與STN 的比值（C/N）也是土壤質量和肥力的重要指標，反映了SOC 與STN 之間的相互作用和耦合關系[2]。此外，C/N 是影響土壤微生物群落的主要因素，因此在陸地碳氮循環中起著關鍵作用[3]。研究土壤C/N 對于土壤管理、生態環境監測和揭示生態系統碳氮耦合關系有重要意義。

傳統的化學分析方法測定土壤C/N 需要分別測定土樣的SOC 含量和STN 含量，再進行計算C/N，雖然精度較高，但是土壤樣品采集勞動強度大，實驗室測定人力財力消耗多[4]。近年來，高光譜技術的發展對于快速，低成本且高精度測定土壤各類理化性質提供了新的思路。目前，很多學者已經使用實測高光譜對SOC 含量和STN 含量進行精確的估測[5,6]，陳秋實利用實測土壤光譜反射率結合互花米草不同入侵年限探究了互花米草入侵下的濱海濕地SOC 含量，模型決定系數（Coefficient of determination，R2）最高達到0.81[7]，盧艷麗對農田耕層土壤的STN 含量進行建模分析，模型R2最高達到0.863[8]。但是將土壤光譜反射率作為自變量，探究土壤的SOC 含量和STN 含量的比值關系相關研究較少。

高光譜數據波段信息豐富，每個波段可以看作是土壤的一種特征，這些特征可以捕捉不同波段下土壤的化學和物理性質，對土壤理化性質的估測有重要意義[9]，但波段間反射率信息具有強共線性，即部分波段信息存在冗余。前人多使用Pearson 相關系數確定敏感波段[10]，李焱等選取了相關系數最大的10 個波段作為自變量，成功預測了新疆地區不同地表利用類型的STN 含量[11]，于雷等比較了全波段與敏感波段建立的預測模型精度，發現敏感波段建立的模型精度略低于全波段建立模型的精度，但全波段模型的復雜度要遠高于敏感波段預測模型的復雜度[10]。目前，降維和特征篩選是處理波段間的多重共線性有效方法[12]，基于這兩種思想，機器學習中的偏最小二乘回歸（Partial least square regression，PLSR）和基于決策樹的方法得到研究人員廣泛重視。

濱海濕地土壤含水量大，可達性差，野外采樣難度大。本文基于有限野外采樣點的SOC 和STN含量數據和實測高光譜數據，并對原始光譜進行光譜變換以挖掘土壤信息，通過降維和特征篩選的方法確定與土壤C/N 相關的敏感波段，最后基于PLSR 和隨機森林回歸 （Random forest regression，RFR）建立濱海濕地土壤C/N 的光譜預測模型。以期為濕地碳氮循環提供幫助。

1 材料與方法

1.1 土壤數據采集分析

研究區位于江蘇省鹽城市大豐區東北部沿海，該地區年降水量為980～1070mm，多集中于夏季，年平均氣溫在13.7～14.6℃之間[13]。研究區蘆葦與堿蓬等本土植物逐漸被外來種互花米草替代，已形成面積廣泛的互花米草單一優勢植物群落，受互花米草入侵與擴散的影響，原有泥灘逐漸演變為鹽沼濕地[14]。本研究考慮到該區范圍內互花米草海向入侵特征，垂直于海岸線設計采樣方案，每隔40 米進行土壤樣本采集，記錄每個采樣點的經緯度坐標，土壤類型。每個樣點處采集3 個表層（0～30cm）土壤樣品，充分混合，形成復合樣品，共獲取33 個土樣，樣點分布如圖1 所示。

圖1 研究區概況圖

為減少外界因素對土樣的影響，土樣室內自然風干，剔除異物，研磨過2mm 篩，分為3 份，分別用于凱氏蒸餾法測定STN 含量[15]，稀釋熱重鉻酸鉀容量法測定SOC 含量[16]，高光譜數據采集。土壤C/N統計特征如表1 所示，數據集偏度為-18.312%，為小偏度數據集且數據分布左偏。

表1 土樣C/N 描述性統計

1.2 高光譜測定及預處理

利用ASD FieldSpec 4 便攜式地物光譜儀測定土壤光譜反射率，其測量范圍是350～2500nm，采樣間隔為1.4nm（350～1000nm）和1.1nm（1001～2500 nm），重采樣間隔為1nm，共輸出2151 個波段。將過篩的土壤樣品放置于深2cm，長寬為5cm 的樣品盒中，均勻平鋪，接觸式探頭獲取該土樣的光譜反射率。在每次測量前進行白板校正，以減小光線漫反射影響，每個土樣采集10 條光譜曲線，ViewSpecPro計算10 條曲線平均值得到該土樣的光譜反射率。

為消除外界環境因素造成的噪聲干擾，本研究采用Savitzky-Golay 算法（多項式階數為2，窗口點數為20）對光譜反射率曲線進行平滑處理。為充分挖掘光譜信息，對平滑后的原始光譜（R）進行數學變換，包括倒數（1/R），對數（LGR），一階微分（R′）和二階微分（R"）[17,18]。

1.3 數據分析方法

1.3.1 偏最小二乘回歸

PLSR 是一種多元回歸分析方法，與多元線性回歸類似，但不同之處在于PLSR 在建立模型時同時考慮了自變量和因變量之間的關系，即通過構建主成分來最大化自變量和因變量之間的協方差[10]。對原始數據進行正交變換構建的主成分是原始自變量的線性組合，且主成分之間是線性無關的，所以PLSR 適用于處理自變量間多重共線性的數據。PLSR 以提取的主成分代替原始變量參與回歸建模，以達到對高維度自變量數據降維的目的[19]。綜上，PLSR 可以解決波段間的共線性問題和樣本數少于光譜波段數的問題，適用處理高光譜數據。本研究依據留一法（Leave-one-out，LOO）交叉驗證計算的預測均方根誤差 （Root mean square error of prediction，RMSEP） 確定主成分數量，RMSEP 值隨主成分個數的增加先遞減，達到最低點后隨主成分個數的增加出現微小上升或波動，為防止模型過擬合，提取最小的RMSEP 對應的主成分作為最終建模變量。PLSR 模型使用R 語言中的“pls”包構建。

1.3.2 Boruta 特征篩選

Boruta 特征篩選的內核是隨機森林算法。首先Boruta 算法會先根據原始數據集創建一個影子數據集，影子數據集與原始數據集有相同因變量，但每一個自變量被隨機打亂，使用兩個數據集同時進行隨機森林建模，比較模型中的變量重要性，如果同一自變量的重要性在兩個模型中有顯著差異，則說明此變量為預測土壤C/N 的敏感變量，反之，則為非敏感變量，最后輸出所有敏感變量，所有的敏感變量被用于RFR 建模[2]。Boruta 特征篩選使用R語言中的“Boruta”包完成。

1.3.3 隨機森林回歸

RFR 采用了集成學習的思想，由多個決策樹組成，每個決策樹是一個弱學習器，通過對這些決策樹的集成，隨機森林對于高維數據和復雜關系的建模效果良好，能夠捕捉非線性、交互作用和噪聲[20]。首先隨機森林從訓練數據集中進行有放回的隨機抽樣，生成多個不同的訓練數據子集。隨機森林在每個訓練數據子集上隨機選擇一個特征子集，隨機森林對每個訓練數據子集和特征子集的組合構建一個決策樹。在得到所有決策樹之后，對于新的輸入數據，每個決策樹產生一個預測結果，取所有決策樹預測結果平均值，作為最終的模型預測結果。所以RFR 對數據異常值不敏感，具有較高的預測準確性和穩定性。隨機森林建模中需要定義兩個參數： 每次拆分時隨機抽樣作為候選特征的數量（mtry）、森林中生成的決策樹數量（ntree），RFR 模型基于R 語言的“caret”包使用網格搜索法調整最優參數，采用均方根誤差（Root mean square error，RMSE）最小的參數組合構建最終模型。

1.4 精度評價指標

本研究建立模型使用LOO 交叉驗證以防止過擬合問題，LOO 交叉驗證使用n-1 個樣本作為訓練集（n 為樣本數），沒有參與訓練的1 個樣本作為測試集，幾乎用到了所有樣本信息，共循環n 次，保證每個樣本都被用作測試集，其充分地利用了每個樣本，更適用于小樣本數據集，驗證結果更具代表性。模型的驗證指標采用預測值和實測值的R2，RMSE 和平均絕對誤差（Mean absolute error，MAE）。R2越接近于1，RMSE 和MAE 越小說明模型精度越高。計算公式為：

2 結果與分析

2.1 不同光譜變換的曲線特征和敏感變量

對R 反射率進行數學變換，各變換后的光譜反射率如圖2 所示。各樣本R 反射率曲線波形基本相似，在可見光波段（350～780nm）光譜反射率迅速上升，在近紅外波段（780～1900nm）上升趨勢減緩，在2100nm 處反射率達到峰值，之后平緩下降。在1410nm 和1920nm 處出現明顯的水分吸收谷。經過對數變換后，曲線形態與變換前基本一致，但相比于R 反射率曲線，各樣本之間的曲線更加聚攏。經倒數變換后的光譜反射率變化趨勢與R 反射率曲線完全相反且更加聚集，但曲線趨勢發生改變的節點未發生變化。R′反射率反映了原始光譜反射率的增減性，經一階微分變換后，放大了原始光譜曲線的上升和下降等特征，曲線波動劇烈。二階微分反映了原始反射率曲線的凹凸性，二階微分光譜反射率在0 周圍浮動。

圖2 土樣的原始與數學變換后的光譜反射率曲線

通過PCA 確定PLSR 建模使用的主成分個數，選擇RMSEP 最小時對應的主成分個數，R，LGR，1/R，R′和R"的主成分個數分別為2，2，2，4 和1。表明經簡單數學變換（LGR，1/R），主成分個數未曾發生變化，經過一階微分后，主成分個數提高為4 個，表明原始光譜中一些微弱的變化，在經微分后得到更強烈的信號，導致主成分分析時提取更多的主成分來解釋這些變化，即一階微分變換挖掘出了更多土壤信息。二階微分反射率在0 周圍浮動，且浮動幅度類似，僅一個主成分即可概括絕大多數自變量信息。

Boruta 特征選擇篩選出用于RFR 建模的變量如表2 所示，R 和LGR 篩選出的特征波段除2450 nm 和2486nm 外完全一致，說明R 進行對數變換對識別敏感波段并未有明顯改善。進行倒數變換后，敏感波段從19 個增加至39 個，敏感波段所處的波段范圍基本一致，均集中在1900～2500nm。經一階微分變換后，敏感波段的數量從21 個提升至38 個，二階微分后，篩選出的敏感波段數量減少到15 個。微分變換后，敏感波段的位置發生較大改變，由未變換前短波紅外范圍擴散至可見光范圍，短波紅外范圍的敏感波段減少，且分布較分散，說明經過微分變換凸顯了可見光范圍的光譜特征，挖掘出更多可見光波段的光譜信息。

表2 Boruta 篩選的敏感波段數量及位置

2.2 模型對比分析

結合PLSR，RFR 和不同光譜變換的波段反射率構建的土壤C/N 預測模型評價指標如表3 所示，不同模型和不同變換形式對模型精度有顯著影響。兩種機器學習算法中，未經微分變換（R，LGR 和1/R）的光譜反射率建立的模型有相似精度，PLSR 模型 中R，LGR 和1/R 的R2分 別 為0.698，0.697 和0.694。表明經對數和倒數變換的光譜反射率對預測模型精度沒有顯著提升，甚至模型精度有輕微降低，RFR 模型中也發現相同規律。經微分變換的光譜反射率建立的模型精度得到提高，PLSR 中，相比于R 建立的預測模型，R′，R" 建立的預測模型R2分別提升了42.6%和33.0%，RFR 中，相比于R 建立的預測模型，R′，R"建立的預測模型R2分別提升了42.9%和42.0%，表明一階微分能挖掘潛在的土壤信息，一定程度上消除噪聲，提高預測模型精度，提高幅度相似證明PLSR 中的主成分提取和Boruta 篩選效果相似。

表3 土壤C/N 不同模型精度對比

對比PLSR 和RFR 模型，PLSR 建立的預測模型精度高于RFR 模型，PLSR 通過提取的主成分解釋波段與C/N 的線性關系，而RFR 對非線性關系有更強的解釋能力，PLSR 算法相比于RFR 算法更適用于建立高光譜數據與土壤C/N 之間的關系。R及其4 種變換后的光譜反射率建立的PLSR 模型比RFR 模型有不同程度地提高，R2提高23.7～27.6%，即每種變換的提高幅度差異較小。預測結果與實測值之間散點圖如圖3 所示，表明預測值與實測值之間有很好的相關性，預測結果是可靠的，PLSR 模型的預測值與實測值的擬合直線相比與RFR 模型的擬合直線更接近于1:1 線，表明PLSR模型預測效果比RFR 模型更好。

圖3 不同光譜變換和不同建模方法預測的土壤C/N與實測值比較

3 結論

高光譜技術依靠其數據獲取快捷，波段信息豐富等優勢被廣泛應用于估測土壤理化性質，但野外環境復雜，受到土壤水分、土壤雜質和大氣條件等影響，背景噪聲大，采集土樣實驗室內測定高光譜數據可以獲得較為穩定的光譜信息，結果相對理想，相比復雜的化學測定土壤C/N 方法，又省時省力無污染。本研究通過對土樣高光譜曲線的重采樣、平均和平滑預處理，并對原始光譜反射率曲線進行4 種數學變換，使用主成分分析、Boruta 特征篩選、PLSR 和RFR 對土壤C/N 分析研究，得到的具體結論如下：

（1）原始光譜R，對數變換LGR 和倒數變換1/R，Boruta 篩選出的敏感波段范圍相似，集中在1900～2500nm，微分變換后篩選出的敏感波段均勻分布在整個波段區間，微分變換能有效挖掘光譜信息。

（2）在PLSR 和RFR 模型中，對R 進行對數和倒數變換并不能提高模型精度，對R 進行微分變換，在兩種算法建立的模型精度均得到了顯著提高，R2最高提升42.9%。

（3）兩種算法均能有效預測土壤C/N，PLSR 模型精度比RFR 模型精度高23.7～27.6%，PLSR 更適用于小樣本數據集建模預測，基于PLSR 和R’構建 的 最 優 模 型，R2為0.995，RMSE 為0.216，MAE為0.165。