基于鐵氧化物特征光譜和改進遺傳算法反演土壤Pb 含量

張 霞，王一博,2，孫偉超，黃長平，張 茂,2

（1. 中國科學院空天信息創新研究院，北京 100101；2. 中國科學院大學，北京 100049）

0 引 言

土壤為自然生態系統的運轉和人類活動提供了基本保障。隨著工業化的不斷發展，通過污水排放、大氣降塵等途徑進入土壤的重金屬元素在土壤中持續累積，造成土壤重金屬污染[1]。Pb 是土壤中常見的一種重金屬污染物，滲透進土壤后移動性較差，殘留時間長，且不易被微生物分解，其治理和修復難度大，不僅影響農作物的生長，還會通過食物鏈威脅人類健康[2]。因此，廣泛開展土壤Pb 含量調查和監測，對保障農業生產、保護人類健康和生態系統安全具有重要意義。近年來高光譜遙感技術[3-4]不斷發展，其具有光譜分辨率高、波段多且連續性強等特點[5]，能夠快速高效地反演土壤重金屬含量[6-8]。另外，隨著高光譜圖像性能的不斷提高，土壤重金屬含量的大范圍制圖以及周期性監測成為可能。

對土壤光譜進行分析時，國內外學者常常使用全譜段來對土壤重金屬的含量進行預測[9-13]。一般情況下，自然環境中的土壤重金屬含量較低，其本身微弱的光譜響應信號又受到多種成像因素的干擾，難以直接通過分析土壤重金屬元素的光譜特征來估算其含量[14]。但土壤光譜活性物質（主要指黏土礦物、鐵氧化物和有機質等）對重金屬元素具有吸附或賦存作用[15]，因而土壤重金屬元素對土壤光譜的影響間接反映在這些光譜活性物質的特征光譜波段上。Sun 等[16-17]使用有機質和黏土礦物的特征波段組合建立Zn 和Ni 的反演模型，與全譜段建模相比精度有顯著提升。賀軍亮等[18]基于土壤有機質敏感波段對應的多種光譜變換指標，建立了Cd 的高光譜間接反演模型。雖然現有研究在間接反演土壤重金屬元素含量方面取得了一定進展，但關于提取土壤 Pb 相關光譜活性物質的特征光譜進行Pb 含量反演的可行性還有待驗證。

在建模算法的選擇方面，遺傳算法-偏最小二乘回歸（Genetic Algorithm-Partial Least Squares Regression，GA-PLSR）[16-17,19-20]具有特征優選的能力，且能對具有強相關性的自變量進行建模，常被用于高光譜數據模型構建。由于高光譜數據波段眾多、相關性強且冗余性高，使用隨機尋優的遺傳算法對參與建模的波段進行篩選能夠減少自變量的個數，降低數據間相關性，提取出光譜的有效信息，最終提高模型估算精度[21-22]。但是遺傳算法仍然存在“過早收斂”的問題[23]，算法在迭代時只接受適應度更高的解，所以在開始迭代后便立即向初始種群的局部解收斂，從而容易陷入局部解區域，難以做到全局尋優，對最終模型的精度產生了一定的影響。

以往研究多以礦區作為研究對象，而礦區以外的一般農作區也應是監測的重要方面。本文以河北雄安一般農作區農田采集的土壤樣本數據為研究對象，提取土壤中對Pb 起主導吸附作用的土壤光譜活性物質的特征譜段進行Pb 含量反演，以避免全譜段建模的數據冗余問題。為解決GA 算法的“過早收斂”問題，提出改進遺傳算法（Improved Genetic Algorithm，IGA），使算法在運行前期能夠跳出局部解區域，尋找更優解，并結合偏最小二乘法[24]進行回歸建模反演土壤Pb 含量。

1 研究數據與預處理

1.1 研究數據

雄安新區是中國新規劃的數字經濟創新發展試驗區，為關注該地區的環境建設，將其作為研究區。2018 年9 月20 日至9 月23 日，于雄縣和安新縣的農田開展土壤樣本采集與野外原位光譜測量試驗，采集農田土壤表層（0～20 cm）作為土壤樣本。野外光譜采集使用SVC HR1024i 地物光譜儀，該光譜儀波長范圍為 350～2 500 nm，測量時探頭離地面距離大約5 cm。采集土壤光譜時每個樣點測3～5 條光譜取平均作為樣點光譜，共采集了土壤樣本及光譜73 組。

土壤樣本送實驗室進行 Pb 含量測試。測量方法為HNO3-HF-HClO4 消煮法，將樣本在微波消解儀進行消解，然后用電感耦合等離子體質譜ICP-MS 測定Pb 元素含量，每個土壤樣本測量3 次后取平均值。

1.2 數據預處理

光譜數據處理主要包括光譜維去噪、異常光譜剔除、光譜特征增強及鐵氧化物特征譜段提取等。

野外光譜測量時，大氣中的水汽在1 400 和1 900 nm處存在強吸收，且1 900～2 500 nm 區間的光譜存在較大噪聲[25]。為去除大氣水汽吸收波段和受噪聲影響的波段，并保留盡可能多的光譜區間，將1 800 nm 之后的光譜波段從土壤野外光譜中剔除。為進一步去除野外光譜在 350～1 800 nm 的隨機噪聲，使用包絡線去除及分段Savitzky-Golay（SG）濾波對光譜曲線進行平滑處理。SG濾波中890～1 020 和1 330～1 520 nm 采用窗口大小為15的二次多項式，其余區間采用窗口大小為7 的二次多項式。為排除光譜異常和重金屬含量異常對試驗結果的影響，通過光譜曲線篩選和z-score 異常點檢測方法，去除3 個異常樣本，最終保留70 個有效樣本用于反演建模。

高光譜數據相鄰兩波段具有較強的相關性，使得整個波段的反射率具有強多重共線性，給反演帶來極大困難。因此采用一階差分的方式，減弱各波段之間的共線性，增強光譜的特征[4]。差分公式見式（1）

式中Bi′為第i段光譜的一階差分，Bi為第i段光譜反射率，n為光譜總數。

2 研究方法

2.1 鐵氧化物光譜特征提取

Covelo 等[26]的重金屬吸附和解析試驗表明：土壤對重金屬Pb 的吸附作用最強，且土壤有機質、鐵氧化物、錳氧化物均對重金屬Pb 具有很強的吸附作用，其中鐵氧化物在重金屬 Pb 的吸附作用中占主導地位，而且隨著土壤中不同重金屬之間對吸附位置的競爭加劇，鐵氧化物在重金屬Pb 的吸附作用中的重要性也隨之上升。Wu等[27]的研究表明，鐵氧化物與重金屬含量相關性在所有土壤光譜活性物質中占主導地位，并試驗得出鐵氧化物的吸收特征在以500 和950 nm 為中心的吸收峰。因此，本文提取500 和950 nm 處的吸收特征譜段預測土壤Pb含量，提取波段區間為[Bm-W/ 4 ,Bm+W/4]，其中Bm是最大的吸收波段，W是吸收區域的寬度。

綜上，從70 條有效土壤光譜中提取500 和950 nm為中心的鐵氧化物吸收峰，提取出的光譜波段范圍為450.7～523.2 以及 9 14.2～1 027.9 nm。

2.2 改進的遺傳波段選擇算法（IGA）

遺傳算法[28-30]仿照物種演化機制，模仿了遺傳進化中發生的基因突變和自然選擇現象，通過隨機選擇、交叉和變異等遺傳操作，產生適應度高的個體作為最終解。傳統遺傳算法中存在最嚴重的問題是“過早收斂”：運行初期會直接向局部解收斂而難以跳出。為解決該問題引入模擬退火算法中的Metropolis 準則[31-32]：若新解優于當前解，則接受新解，否則以某一概率接受新解，且此概率隨適應度增大而降低。Metropolis 準則使得算法在運行初期有一定的概率接受較差解而跳出當前局部解區域，從而做到全局尋優。IGA 算法具體流程如圖1。

圖1 改進的遺傳算法流程圖Fig.1 Process of improved genetic algorithm

算法的改進核心為生成新種群，其中具體的Metropolis準則為：當新個體目標函數值小于父個體目標函數值時，用新個體取代父個體；否則以 exp （ -k? Δf/E）的概率接受新個體為父個體。其中k為大于0 的實數，E為種群目標函數值的方差。令n為算法總迭代次數，當E大于0.01/n時，k=E+ 0 .01/n；否則k= 0 .01/n。k保證了E大于0.01/N時算法有足夠概率能夠跳出局部解區域，E小于0.01/N時算法能夠保持收斂狀態。E定義如式（2）

式中N為種群個體數，fi為個體i的目標函數，f為種群目標函數的平均值。

2.3 偏最小二乘建模與模型評價指標

偏最小二乘法（Partial Least Squares Regression，PLSR）能較好地解決高光譜反演中自變量之間的強多重共線性問題，并且能夠實現多對多線性回歸建模，特別是在自變量存在多重相關性而且個數較多，而樣本數目又遠小于自變量的情況下PLSR 仍適用[33-35]。因此，本研究采用PLSR 建模。

使用擬合優度R2，均方根誤差（Root Mean Square Error，RMSE）和相對偏差（Relative Percent Difference，RPD）評價模型的優劣。較好的模型通常擁有較高的RPD 和R2及較低的RMSE。模型的評價參考現有的土壤屬性含量高光譜估算的評價標準[19,35]：出色模型，R2> 0.9，RPD>3.0；良好模型，0.9 >R2> 0.82，3.0 > RPD > 2.5；近似模型，0.82 >R2>0.65，2.5 > RPD > 2.0；具有一定估算能力，0.65 >R2> 0.50，2.0 > RPD > 1.5；不具備估算能力，0.50 >R2，RPD < 1.5。

3 結果與分析

3.1 樣本集劃分

將 70 個樣本分為訓練集和測試集。依據重金屬 Pb的含量從小到大進行排序，對樣本集按照2∶1 的比例進行分層抽樣，獲得47 個樣本的訓練集和23 個樣本的測試集，樣本集的數據統計見表1。

表1 樣本集中Pb 含量數據統計Table 1 Data statistics of Pb content in the sample set

3.2 模型構建

為驗證 IGA 波段選擇算法改進的有效性，本文對GA-PLSR 和IGA-PLSR 這2 種算法進行建模。2 種算法的參數設置相同，個體數為 20，個體編碼方式為二進制編碼，代際差為90%。將RMSE 作為波段選擇算法的目標函數，即算法目標是最小化 RMSE。將選擇算法運行結果中RMSE 最小的個體作為當前最優解輸出，并將其作為PLSR 的自變量建立最終的回歸預測模型。最后使用測試集數據計算模型評價指標。

3.3 土壤Pb 含量估算模型評價及分析

根據2.1 節的理論與方法，提取出鐵氧化物特征譜段（450.7～523.2 nm 和914.2～1 027.9 nm）作為輸入，使用GA-PLSR 及其改進算法IGA-PLSR 對土壤Pb 含量進行反演建模。為對比兩算法的運行效果，選取迭代次數為100、250、500 和1 000，每種迭代次數下運行5 次算法，評價參數及運行時間取5 次運算的平均值。

試驗結果如表 2 所示，可以看到隨著迭代次數的增加，2 種算法的平均耗時增加，同時平均精度也在不斷提高，IGA-PLSR 因加入了Metropolis 準則的篩選過程，相同迭代次數下耗時略大于GA-PLSR 算法。由模型精度隨迭代次數的變化可得，2 種算法的最高精度對應的迭代次數均為1 000 次，這與Wang 等[36]應用GA 反演土壤重金屬的迭代次數設置一致。在迭代次數為 1000 次時，GA-PLSR 算法的反演精度R2和 RPD 平均值分別為0.782、2.117，RMSE 的平均值為2.487 mg/kg，滿足了近似模型的條件，能夠近似估計Pb 含量；而IGA-PLSR 的3 種評價參數均要優于改進前的算法，模型的R2和RPD為 0.822、2.377，RMSE 為 2.221 mg/kg，其中R2已經達到了良好模型的標準，表明模型對土壤Pb 含量具有良好的估算性能。因此，IGA 算法相對于傳統的 GA 算法具有更優的波段選擇能力，顯著提升了反演模型的精度。

表2 基于鐵氧化物特征譜段建模的反演精度Table 2 Inversion accuracy based on iron oxide characteristic spectrum modeling

圖2 分別展示了2 種算法在1 000 次迭代時的5 次運算中的最優結果。可以看到IGA-PLSR的R2、RPD為0.837、2.403，RMSE 是 2.063 mg/kg，顯著優于傳統算法 GA-PLSR（R2、RPD 分別為 0.788、2.140，RMSE 為 2.447 mg/kg），且通過散點的分布可以看到IGA-PLSR 模型中的樣本點更加集中于1∶1 線附近。因此，IGA-PLSR 算法的反演能力顯著優于改進前的GA-PLSR 算法。

圖2 最優Pb 含量估算結果Fig.2 Optimal results of Pb content pridiction

3.4 鐵氧化物特征譜段的有效性分析

為對比分析鐵氧化物特征譜段的有效性，基于350～1 800 nm 的全部光譜波段對土壤重金Pb 含量進行建模反演，根據3.3 節模型分析結果，設置迭代次數為1 000 次，表3 展示了試驗結果。

表3 基于全譜段建模的反演精度(迭代次數設為1 000 次)Table 3 Inversion accuracy based on full spectrum modeling(1 000 iterations)

對比表2（鐵氧化物特征譜段）和表3（全譜段）的結果可以看出，2 種建模方法在采用全譜段建模時模型精度較低。GA-PLSR 模型和IGA-PLSR 模型的平均R2分別為 0.292、0.440，平均 RPD 分別為 1.216、1.366。IGA-PLSR模型有較高的精度和反演能力，表明IGA 波段選擇算法在光譜波段較多、冗余度較高的情況下依然能選擇出相對有效的光譜波段。但此時 2 種模型均未達到估算的精度要求，而表 2 中使用鐵氧化物特征譜段建模的結果要顯著優于表 3 全譜段的結果，說明鐵氧化物特征波段的提取能切實有效地提高Pb 含量反演的精度。

為直觀顯示鐵氧化物特征波段對土壤Pb 含量預測的貢獻度，圖3 展示出全譜段IGA-PLSR 模型中IGA 算法的波段篩選結果。波段選擇算法的作用是篩選出能有效反演土壤重金屬含量的波段，可以看到篩選的波段在鐵氧化物的特征譜段（500 和950 nm 吸收峰）分布最為密集，表明對研究區土壤Pb 含量進行建模反演時，鐵氧化物特征譜段相比其他波段更加有效。

以上試驗均說明了提取相關土壤活性物質的光譜特征能夠在一定程度上減少光譜信息的冗余，從而提高Pb含量反演的精度。本試驗為基于鐵氧化物特征波段反演土壤Pb 含量提供了試驗依據，同時為反演土壤中其他重金屬含量提供方法參考。

圖3 全譜段IGA-PLSR 模型波段篩選結果Fig.3 Band selecting results of the IGA-PLSR model with full spectrum

3.5 算法收斂性及不確定性分析

為進一步分析 IGA-PLSR 算法對“過早收斂”的改進，選擇GA-PLSR 和IGA -PLSR 迭代1 000 次時的模型，繪制代價函數 RMSE 隨迭代次數變化的曲線，如圖 4。該曲線圖將算法尋優過程可視化，以便對比 2 種算法運算過程中的異同。向局部最優解收斂的特點為緩慢逼近，由圖4 可看到GA-PLSR 模型由于僅接受更優解，整體呈單調下降趨勢，在250 次迭代后RMSE 基本穩定，隨后400、600、720 及900 次迭代鄰域中都僅僅在某一次迭代后RMSE 值有極小的降低，表明已陷入局部最優，后續僅是逼近局部最優解；IGA-PLSR 算法在迭代初期有較高的概率接受較差解，RMSE 波動較大，在 300 次后開始收斂，每次迭代都有一定的概率接受較差解，從而跳出當前的局部解區域，如在380 和640 次迭代附近接受較差解后RMSE 都會有大幅降低出現，說明跳出了當前局部解區域從而進一步尋找更優解。RMSE 曲線圖的對比表明了 IGA-PLSR 算法提高反演精度的原因，凸顯了其在全局尋優方面的優勢。

改進前后的GA 和IGA 都是隨機尋優算法，故存在一定的不確定性。主要體現在隨機產生初始解、隨機選擇、交叉和變異等一系列遺傳操作中，同時引入的Metropolis 準則中依概率接受較差解也增加了模型的不確定性。因此在實際應用時，應多次運行算法后選擇最優的模型進行重金屬含量的反演。

圖4 不同算法RMSE 隨迭代次數的變化Fig.4 Changes of RMSE with iterations of different algorithms

4 結論與展望

本文基于河北雄安一般農作區采集的70 個土壤樣本的野外光譜數據，研究具有機理性的特征波段的選取方法對土壤Pb 含量進行反演。首先，為規避全譜段建模的數據冗余和機理性不足問題，提取出土壤中對Pb 具有主導吸附性作用的鐵氧化物的特征波段用于Pb 含量反演，在此基礎上，通過引入Metropolis 準則對GA 波段選擇算法改進并提出改進遺傳算法（Improved Genetic Algorithm，IGA），解決遺傳算法的“過早收斂”的問題，利用IGA 算法篩選鐵氧化物特征譜段中更加有效的光譜波段組合，采用偏最小二乘法（Partial Least Squares Regression，PLSR）構建反演模型。結論如下：

1）IGA-PLSR 模型相對GA-PLSR 擁有更高的反演精度，其最優模型的R2、RPD 分別為0.837、2.403，RMSE為2.063 mg/kg，其中R2達到了良好模型的標準，可用于雄安新區農田土壤Pb 含量的定量估算。

2）鐵氧化物特征譜段的反演精度要顯著優于全譜段的精度，說明了依據土壤Pb 的吸附機理而提取相關土壤光譜活性物質的特征譜段能夠有效減少全譜段建模的冗余信息，從而提高Pb 含量反演精度。

3）通過對2 種算法運行過程中RMSE 的變化曲線分析，表明IGA-PLSR 能夠有效地解決GA 算法的“過早收斂”問題，能夠跳出局部解區域尋找更加有效的光譜波段組合，提高PLSR 的建模精度。

綜上所述，在應用于實際區域土壤Pb 含量快速檢測時，可以從兩方面提高土壤元素含量反演模型的精度和適用性：提取出與土壤元素相關的光譜活性物質的特征光譜代替全譜段進行建模反演；改進現有 GA 算法使用IGA-PLSR 算法進行反演建模。但是，在野外光譜及航空/衛星高光譜圖像獲取過程中，不可避免地會受到土壤粒徑、含水率等環境因素影響，增加土壤重金屬高光譜反演的不確定性，因此，下一步將研究野外光譜乃至像元光譜的環境因素去除方法，以期提高土壤重金屬高光譜遙感反演方法的穩定性。