基于隨機變異—Kennard—Stone和偏最小二乘法的土壤重金屬鎘含量反演
——以雄安新區西南部為例

黃照強，倪斌

中國冶金地質總局礦產資源研究院，北京, 101300

內容提要: 土壤重金屬污染對人類身體健康造成了嚴重威脅，作為快捷、高效、無損監測分析土壤重金屬含量的方法之一，高光譜遙感反演土壤重金屬含量的方法正逐步得到發展。本文以雄安新區西南部作為研究區，針對潛在生態風險最大且活性最強并易被植物吸收的重金屬鎘，開展高光譜反演研究。將采集來的426件土壤樣品除雜、風干、過篩在實驗室測得重金屬含量，并用SVC便攜式地面光譜儀測量樣品350～2500 nm范圍光譜。采用Savitzky—Golay卷積平滑方法進行光譜降噪平滑處理。由于粒徑大小而不是化學成分差異可能會導致基線效應和漂移現象，為了增強光譜差異和光譜曲線形狀，將數據進行標準正態變量變換(SNV)、一階微分(FD)、二階微分(SD)、多元散射校正(MSC)等數學變換，并分析變換后光譜與Cd含量的相關性。本文提出一種集成隨機變異法—Kennard—Stone法—偏最小二乘回歸的方法。針對變換后的光譜集和Cd含量集，第一步采用隨機變異法—Kennard—Stone法將樣本集分為70%訓練集和30%驗證集，使樣本數隨性質分布均勻并覆蓋整個樣本空間；第二步用偏最小二乘回歸法結合交叉驗證建立回歸模型，用確定系數R2、均方根誤差RMSE、偏差比值RPD、誤差范圍比值RER等參數開展模型評價，如果沒有達到預期效果，則回到第一步，迭代反復選擇，直至達到最優效果。結果表明，適用的最優反演重金屬鎘含量的技術方法是采用FD變換后，不斷迭代集成RM—KS樣本選擇和PLS偏最小二乘法回歸建立的模型，其驗證綜合效果最好，建模主成分數為11個，確定系數R2達到0.909，RMSE為0.604，RPD為2.696，RER達為15.516。成果可為類似區域快速、無損的土壤重金屬Cd含量反演提供技術支撐。

土地資源是人類賴以生存的最基本自然資源(朱永官等，2015)。隨著全球經濟的迅猛發展，工業化、農業現代化進程中所引起的工業污染和農田施肥農藥造成了污染越來越嚴重，特別是重金屬污染物越來越多的進入到農業土壤中。農業土壤中的重金屬離子較易被農作物吸收并形成富集，會對人類身體健康產生威脅。因此，了解土壤中的重金屬濃度將是非常重要的，以便有效進行農用地土壤風險篩查和管制。土壤重金屬含量和農田土壤環境質量評價及治理一直是國內外研究的難點和熱點問題，是環境科學的重要研究領域之一(Pkasads，1995；陳玉娟等，2005；崔邢濤等，2014)。全球工業化導致大量潛在有毒化合物釋放到生物圈中，其中包括微量元素，如鎘、汞、鉛、砷、鋅和鎳等，這些元素通常被稱為重金屬(Lasat et al., 1998)。在過去的50 a里，全球重金屬釋放量達到了22 kt Cd、939 kt Cu、783 kt Pb和1350 kt Zn(Singh et al.,2003)。2014年全國土壤污染狀況調查公報結果顯示，我國耕地土壤環境質量堪憂，耕地土壤點位超標率達19.4%，主要以鎘、鎳、銅、砷、汞、鉛等重金屬污染為主，其中鎘(Cd)污染最為嚴重，點位超標率達7.0%(環境保護部和國土資源部，2014)。農田土壤Cd主要來源為大氣沉降、畜禽糞便(有機肥)、化肥/農藥、污灌及其他(磷肥、餅肥等)等,全國范圍農田土壤中Cd的年輸入通量約為0.725g/(hm2·a), 河北省農田年輸入通量則達到2.153 g/(hm2·a), 為全國農田土壤中Cd年輸入通量的2.97倍(陳世寶等，2019)。近年來，國內外眾多學者對土壤重金屬污染的狀況和成因進行了研究，通過多目標地球化學調查，系統分析土壤重金屬污染特征, 采用內梅羅綜合污染指數法、富集因子、地累積指數、潛在生態危害指數和健康風險評估模型對土壤重金屬污染程度、潛在生態風險和健康風險進行評估(樊燕等，2008；崔邢濤等，2014；張素榮等，2020；郭志娟等，2021)。然而，傳統的獲取土壤重金屬含量和空間分布的方法是采用網格化采樣和費時的實驗室分析，這具有破壞性，然后再采用地統計學插值得到結果(Leenaers et al.,1990; Kemper and Sommer,2002)。反射光譜法是一種快速無損表征不同材料的方法(Ben-Dor et al.,1994， 2009; Clark, 1999)。隨著高光譜遙感技術的發展，其快捷、高效、無損、易獲取等特點將為土壤遙感監測提供有力支撐。可見光—短波紅外分析方法獲得的精度有助于快速土壤表征和遙感應用(Ben-Dor and Banin, 1994)。金屬在近紅外區域沒有吸收特征，但可被檢測，因為它們與有機物絡合并與可檢測的部分如氫氧化物、硫化物、碳酸鹽或氧化物相關聯，或被具有光譜吸收特征的黏土吸附(Malley and Williams, 1997)。Kemper等基于土壤光譜與重金屬間的相關關系，采用多元線性逐步回歸和人工神經網絡方法建立反演模型高精度的預測出礦區土壤中As、Fe、Hg、Pb、S、Sb等的含量(Kemper and Sommer,2002)。Kooistra等采用偏最小二乘回歸和反射光譜建立了萊茵河流域土壤中的鎘、鋅含量的預測模型(Kooistra et al., 2001)。Sun Wei等用反射率偏態分布校正和偏最小二乘法估算礦區重金屬Cu濃度(Sun Wei et al., 2019)。Wu Yunzhao等利用反射光譜法研究了農業土壤中的汞污染，通過地面采集的光譜數據做反演模型研究，發現汞濃度分別與針鐵礦和粘土在496 nm和2210 nm處的吸收特征相關(Wu Yunzhao et al.,2005)。Chen Wei等利用紫外—可見/熒光光譜結合相關分析，研究了SOM與Pb2+和Cd2+相互作用過程中的相干結構變化(Chen Wei et al., 2020)。眾多研究者都是利用多種光譜變換如多元散射校正、一階微分、二階微分、連續統去除等用于提取可識別的光譜敏感波段和/或光譜轉換指數，然后采用各種線性/非線性的分析方法如逐步多元線性回歸、偏最小二乘回歸和神經網絡等從這些變換后的光譜來建模預測土壤重金屬含量(Liu Meiling et al.,2012; Wang Fenghe et al.,2018; Riedel et al., 2018; 鄒濱等，2019)。逐步多元線性回歸方法是選擇幾個特征波段進行建模且特征波段數不能大于樣本數，而提取高光譜數據中的幾個特征波段對于高光譜數據豐富的光譜信息是極其有限，其重要信息易丟失(Lee et al.,2004; 劉克等，2012)，故在分析中需盡量保持高光譜數據的豐富內涵。偏最小二乘回歸則特別適合當變量比觀測樣本值更多，以及變量中具有強共線相關性、噪聲等問題的情況(Wold et al.,2001)。

因此本文以雄安新區西南部作為研究區，針對潛在生態風險最大且活性最強并易被植物吸收的重金屬鎘Cd(張素榮等，2020)，開展土壤中重金屬Cd的含量和空間分布研究，通過采用多種光譜曲線變換方法，消除噪聲等不利因素，提高變換后光譜特征與重金屬鎘Cd含量之間的相關性，然后提出一種集成隨機變異—Kennard—Stone —偏最小二乘回歸(Random Mutation — Kennard—Stone — Partial Least Squares，RM—KS—PLS)的方法，不斷循環迭代選擇合適樣本訓練集和測試集，使得樣本能均勻覆蓋整個樣品空間，建立預測模型，并開展精度對比研究，探索適用的最優反演重金屬鎘含量的預測模型。

1 研究區與實驗分析

1.1 研究區概況

研究區位于雄安新區西南部，行政區劃歸屬于安新縣，包括何橋鄉老河頭鎮和蘆莊鄉，屬暖溫帶季風型大陸性半濕潤半干旱氣候(朱江等，2017)。區內交通十分便利，分布有331省道、336國道等交通要道。該區地勢相對平坦，地形地貌為沖積洼地平原，地勢低平，土壤母質主要為黃土母質和洪沖積物母質為主，主要為潮土、及少量褐土和沼澤土，唐河和清水河流經該區域。農業主產玉米、小麥等大宗農作物。區內冶煉、有色金屬等工業產業發達。在一段時間內存在數以百計規模大小不等的多金屬冶煉廠和無庫尾礦渣堆，大氣沉降、降雨淋溶、大風吹揚等導致大量有毒有害物質進入土壤耕作層，造成部分耕地土壤環境Cd、Cu、Zn、Pb等重金屬元素不同程度超標(崔邢濤等，2014；張素榮等，2020；郭志娟等，2021)。

1.2 樣品采集

采樣點主要在研究區的農田，少數在建筑用地和尾礦渣堆(圖1)。在研究區采集0～20 cm耕作層土壤樣品，為提高樣品的代表性，在每個采樣點50 m范圍內采用梅花形布點法布設5處子采樣點，去除石塊、雜草、樹根等雜質，然后采用四分法合成一個混合樣，樣品重量大于1 kg，樣品裝在聚乙烯樣品袋中，并獲取WGS84坐標系統的GPS經緯度定位坐標。同時，實地調查和詳細記錄了采樣點土壤性狀和周邊環境。

圖1 雄安新區西南部土壤重金屬采樣點分布圖Fig. 1 Distribution of soil heavy metal sampling sites in southwest of Xiong’an new district

1.3 樣品測試

樣品的測試分析由中國冶金地質總局地球物理勘查院測試中心完成。樣品帶回測試中心后，經風干、研磨、過篩(200目)，然后采用電感耦合等離子質譜儀NexIoN 350x型測定重金屬的含量。根據國家標準《硅酸鹽巖石化學分析方法 第30部分：44個元素量測定》(GB/T 14506.30-2010)測定了土壤重金屬Cd的含量，中國土壤背景值Cd為0.097 μg/g(中國環境監測總站，1990)，冀中平原背景值Cd為0.152 μg/g(張素榮等，2020)，研究區的土壤pH平均值約為8.0，偏堿性，因此其農田土壤污染風險篩選值Cd為0.6 μg/g，管控值為4 μg/g(生態環境部和國家市場監督管理總局，2018)。

1.4 樣品光譜測量

眾多研究學者在研究類似問題時，都是首先采用實驗室測量和建模，在模型成熟后，再應用到實際航空高光譜和衛星遙感數據中。本文也是采取此思路開展模型研究。土壤樣品光譜測量采用已經過風干、研磨和過篩(200目)的樣品，用美國SVC Hr1024便攜式地面光譜儀進行測量，鹵素燈作為光譜測量光源，光譜采樣范圍是350～2500 nm(VNIR)，共1024個波段，其中350～1000 nm范圍的光譜分辨率為3.5 nm，1000～1850 nm范圍的光譜分辨率為9.5 nm，1850～2500 nm范圍的光譜分辨率為6.5 nm。測量前用白色參考板輻射定標，每個樣品測量5次，然后取其算術平均值。受環境影響，小于400 nm紫外光譜的噪聲較大，因此去除小于400 nm的波段，只采用400～2500 nm范圍的光譜，共924個波段。

1.5 光譜處理

預測質量在很大程度上取決于反射光譜的預處理步驟和方法(Kemper and Sommer,2002)。測量的光譜曲線常存在噪聲，為了有效地消除光譜噪聲，同時保持光譜表達的化學信息，采用Savitzky—Golay卷積平滑方法(S—G，5個點平滑)進行光譜平滑處理，該方法是一種在時域內基于局域多項式最小二乘擬合的濾波方法(Savitzky and Golay, 1964)。另外由于粒徑大小而不是化學成分差異可能會導致基線效應，因此為了獨立于強度標度的數據比較，增強吸收特性和光譜曲線形狀(Mackin et al.,1991)，將數據進行標準正態變量(SNV)、一階微分(FD)、二階微分(SD)、多元散射校正(MSC)等變換。

標準正態變量 (Standard Normal Variate，SNV) 變換是一種面向行的變換，通過對每條光譜進行集中和縮放來消除因固體顆粒大小、表面散射及光程變化引起的光譜散射效應。SNV使用來自該光譜的數據標準化每個光譜。數據X行中的每個值xk根據以下公式進行轉換：

(1)

其中xk,SNV表示轉換后的值，Mean(X)表示數據X行的平均值，SDev(X)表示數據X行的標準差。

微分(Differential)是光譜應用中廣泛使用的一種經典技術，用于校正光譜中的基線效應，以消除非化學效應并建立穩健的校正模型。一階或二階微分處理后，“隱藏”在光譜中的一些信息可能更容易被揭示。光譜的一階微分(First Differential，FD)只是光譜曲線在每個點的斜率的一個度量。曲線的斜率不受光譜中純加性基線偏移的影響，因此一階導數能最小化反射光譜的背景干擾和基線漂移，是消除此類偏移的非常有效的方法。二階微分(Second Differential，SD)是曲線斜率變化的量度。除了去除純加性偏移外，它不受數據中可能存在的任何線性“傾斜”的影響，因此是從光譜中去除基線偏移和斜率的非常有效的方法，二階微分有助于分辨附近的峰和銳化光譜特征，但對噪聲更敏感。

一階微分公式為：

(2)

y′(i) 指波段i處的一階微分，y(i+1) 為下一個采樣波段的反射率，y(i-1) 是上一個采樣波段的反射率,Δy是采樣間隔。

二階微分公式為：

(3)

y″(i) 指波段i處的二階微分，y'(i+1) 為下一個采樣波段的一階微分，y′(i-1) 是上一個采樣波段的一階微分,Δy是采樣間隔。

多元散射校正(Multiplicative Scatter Correction, MSC) 是一種用于補償光譜數據中加性和/或乘性效應的變換方法，用來消除樣本間的基線平移和漂移現象，增強光譜差異性。首先計算所有樣本的平均光譜，然后將平均光譜作為標準光譜，進行一元線性回歸求得各光譜相對于標準光譜的線性平移量(回歸常數)和傾斜偏移量(回歸系數)，最后在每個樣品原始光譜中減去線性平移量同時除以回歸系數，得到修正后光譜，提高光譜信噪比。多元散射校正公式：

(4)

其中Ai是樣品光譜，mi和bi分別表示各樣品光譜與平均光譜進行一元線性回歸得到的相對回歸系數和平移量。

1.6 模型校準和驗證

在建模階段，需要將化學定量數據與光譜曲線組合來建立化學定量模型。在模型構建前，樣本被分成校準集(訓練集，70%)和驗證集(測試集，30%)。利用隨機變異法和Kennard—Stone法(RM—KS)進行訓練集和驗證集樣本的挑選，可明顯改善樣本數隨性質分布的均勻性并覆蓋整個樣本空間，并在生物信息領域得到應用(Kennard and Stone,1969; Morais et al.,2019b)。KS算法是基于歐幾里德距離的計算，首先選取與所有其它樣本距離最大的樣本，接著選取離所選樣本盡可能遠的樣本，直至達到所選樣本數。然后在KS結果中使用隨機變異因子10%，將訓練集中的一些樣本轉移到測試集中，將測試集中的一些樣本轉移到訓練集中，其中隨機變異因子設為10%是一個常用閾值，用于在隨機度和模型收斂性之間保持平衡(Morais et al.,2019a,b)。

偏最小二乘回歸(Partial Least Squares Regression ，PLSR)(Wold et al.,2001)是一種通過線性多元模型將兩個數據矩陣X和Y關聯起來的方法，能同時對X和Y進行建模，以找到X中最能預測Y中潛在變量的隱藏變量，特別適合當預測矩陣比觀測值有更多X的變量，以及X的值中具有強共線相關性、噪聲。其核心思想是在X與Y中提取出成分t1和u1，滿足t1和u1應盡可能大地攜帶他們各自數據表中的變異信息，以及t1和u1與的相關程度能夠達到最大，則此成分確定；不斷反復迭代，選出相應的成分，直至達到滿意的精度。第一個潛在變量可以通過最大化預測變量和響應變量的協方差cov(t1,u1) 來獲得，公式如下：

max [cov(t1,u1)]=

(5)

其中載荷權重w1,c1分別是對稱矩陣X’YY’X，Y’XX’Y的特征向量，X’、Y’分別是X、Y的轉置矩陣。||w1||=1 表示w1的1范數，也就是絕對值的和是1。第一個潛變量定義為：t1=Xw1。利用最小二乘法對X和Y進行壓縮。根據主成分回歸思想可以把X和Y分別對它們的主成分t1和u1進行回歸建模：

(6)

而本研究中作為自變量的光譜曲線具有幾百乃至上千個波段，往往超過了觀測樣本數，因此正可以利用偏最小二乘回歸法的適合于當預測矩陣比觀測值有更多X變量的特點，以變換后光譜訓練集與重金屬鎘Cd含量之間建立校準模型，用驗證集進行驗證以便獲取最佳預測模型。本研究集成一種RM—KS—PLS方法，用Matlab軟件開發功能模塊，進行不斷循環迭代以便獲取最佳的預測模型。模型建立流程見圖2。

圖2 模型建立流程Fig. 2 Flow of Modeling

2 結果與討論

2.1 土壤樣本重金屬Cd含量分析

通過對研究區土壤樣本初步分析，剔除異常值數據，則選擇了426個樣本劃分為模型訓練集和測試集，按照KS算法+隨機變異因子法，得到訓練集樣本297個和測試集樣本129個。通過上述的土壤樣品采集和實驗室分析測試得到了研究區土壤樣本重金屬Cd含量值，其基本特征統計見表1，重金屬Cd含量平均值為1.31 μg/g高于冀中平原背景值和土壤污染風險篩選值，同時根據單因子污染指數評價法(張素榮等，2020)對重金屬Cd累計污染程度進行評價，計算公式為：

表1 雄安新區西南部土壤重金屬Cd含量統計表Table 1 Heavy metal Cd concentration (Range, Average) of soil in southwestern Xiongan area

(5)

其中Ci是污染物的實測值，Si指重金屬Cd的風險管理目標值即風險篩查值，pi指污染指數，pi>1表明土壤受到污染，越大則重金屬累積污染程度越大，其中污染指數在1～3為輕度污染，3～6.7為中度污染，大于6.7為重度污染。樣本分析結果顯示198個點無污染，輕度污染有158個點，中度有35個點，重度污染有35個點(圖3)。故表明研究區存在土壤重金屬Cd污染。而且數據涵蓋了從清潔到重度污染，為建模提供了可靠的基礎數據和驗證集。

圖3 土壤樣品Cd污染情況分類圖Fig. 3 Soil Samples Cd Contamination classification

2.2 土壤光譜分析

圖4 土壤樣品的SVC光譜曲線圖Fig. 4 SVC Spectral Curve of Soil Samples

SVC光譜反射率數據經過S—G平滑后，采用FD、MSC、FD+MSC、SD、SNV五種方法進行光譜變換，得到的結果與Cd含量進行相關分析。FD變換后光譜與Cd含量的相關系數R較大且顯著相關的有R為0.48在1911 nm波長(水吸收)，-0.42在1721 nm波長(有機質吸收)，-0.37在2365 nm波長(碳酸鹽吸收)，0.34在466 nm波長(鐵錳氧化物吸收)；MSC和FD+MSC變換后光譜與Cd含量相關系數絕對值較大的且顯著相關在0.5～0.6之間，在570(氧化鐵吸收特征)、2200和2300 nm(黏土礦物吸收特征)高度相關；SNV和SD變換后光譜與Cd含量相關系數絕對值較大且顯著相關的在0.4～0.5之間(圖5)。顯著相關指示了變換光譜與Cd含量之間關系建模過程中，這些光譜波長具有很大的貢獻。

圖5 FD、MSC、FD+MSC、SD、SNV變換光譜波長與Cd含量相關關系Fig. 5 Correlation Coefficient between wavelength of FD, MSC, FD+MSC, SD, SNV transformed Spectrum and Cd concentration

2.3 PLS模型預測

各種光譜變換后的曲線，采用PLS建模和驗證得到的結果見表2所示。將預測值與化學分析的參考值進行相關分析，常采用擬合優度(Goodness of Fit)來表示回歸直線對觀測值的擬合程度。度量擬合優度的統計量是確定系數R2，其取值范圍是[0, 1]，越接近1，說明回歸直線對觀測值的擬合程度越好，反之則越差。微量元素在自然背景環境中含量很低，多呈對數正態分布，而在人為污染情況下可能大量富集，造成含量的概率分布向高濃度方向偏斜(Zhang and Selinus，1998)。因此，需要計算另外的參數值共同來更好的評價模型精度，即均方根誤差RMSE(越小越好)、偏差比值RPD(Ratio of percent deviation，為參考值標準差與RMSE的比值，RPD>1.4效果較好)、誤差范圍比值RER(Ratio of error range，為預測值的取值范圍與RMSE的比值，RER>10效果較好)(Malley and Williams，1997)等。RPD<1.0表示非常差的模型，不建議使用它們；RPD在1.0和1.4之間表示差的模型，其中只有高值和低值可區分；RPD在1.4和1.8之間表示恰當的能用于評估和相關分析的模型/預測；RPD值在1.8和2.0之間表示模型良好，可以進行定量預測；RPD值>2.0之間表示非常好的定量模型/預測(Viscarra et al.,2006)。

表2 RM-KS-PLS建模和驗證的結果 Table 2 Calibration and Validation Results For RM-KS-PLS

5種光譜變換后建模和驗證的結果見表2，F檢驗的P值都小于0.01，具有顯著相關性。光譜FD變換后，建模主成分數為11個，驗證集(Validation)預測值與驗證集化學實測含量之間的相關關系確定系數R2達到0.909， 為這幾種光譜變換方法中最大的，RMSE最小為0.604，RPD最大為2.696，RER大于10為15.516(表2，圖6)，表明模型驗證很好，但有些高值被低估了。而光譜FD+MSC變換后，建模主成分數為11個，驗證集預測值與化學實測值之間的確定系數R2為0.858，RMSE為0.622，RPD為2.298，RER大于10為21.177，模型驗證良好程度次之。SD變換后，驗證集預測值與化學實測值之間的確定系數R2為0.831，RMSE為0.593，RPD為2.080，RER為12.350；MSC變換后，驗證集預測值與化學實測值之間的確定系數R2為0.775，RMSE為1.238，RPD為1.310，RER為17.922，有些值被高估了；這兩類總體效果良好。如果剔除一些異常值，結果會更好。光譜經過SNV變換后，驗證集預測值與驗證集化學實測含量之間的確定系數R2小于0.680，效果較差些。五種光譜變換后的PLS校準集(訓練集)預測值與參考值之間的一致性都較好，FD中的部分高值被低估，部分低水平樣本的預測值為負值；而FD、MSC、SD、FD+MSC變換后PLS驗證集預測值與參考值之間的一致性較好，SNV一致性相對偏低(圖6)。

總體效果是FD變換后，采用RM—KS方法進行樣本集選擇和PLS方法迭代循環建模得到的預測效果最好，效果排在其次的是FD+MSC變換、SD和MSC變換，最后是SNV變換的效果。

3 結論

通過對雄安新區西南部受到不同程度重金屬污染的區域采集土壤樣品，開展重金屬污染監測，分析光譜與Cd含量的相關關系，采用集成的隨機變異法—Kennard—Stone法(RM—KS)—偏最小二乘(PLS)方法不斷迭代交叉驗證優選訓練集和驗證集以使樣本分布均勻并覆蓋整個樣本空間，建立PLS模型實現土壤重金屬含量反演，探索適用的最優反演重金屬鎘含量的預測模型。得出結論如下：

(1)樣品SVC光譜反射率數據經過S—G平滑后，采用FD、MSC、FD+MSC、SD、SNV五種方法進行光譜變換，消除因粒徑大小導致基線效應，增強了光譜響應吸收特征和光譜曲線形狀。

(2)集成隨機變異法—Kennard—Stone法(RM—KS)—PLS方法，不斷迭代交叉驗證，優選訓練集(70%)和驗證集(30%)，改善樣本選取的均勻性和覆蓋全域，有利于提高模型效果；通過確定系數R2、均方根誤差RMSE、偏差比值RPD、誤差范圍比值RER等參數比較，表明采用FD變換后，RM—KS—PLS方法建模效果最好，為將來利用光譜方法快速、無損實現土壤中重金屬Cd含量的反演提供支撐。