礦區(qū)土壤重金屬Cd高光譜建模與模型適應(yīng)性分析

蘭 淼，楊 斌，宋 強，陳弘揚，周鵬飛，莊紅娟，方 兵，張世文*

(1.安徽理工大學(xué) 地球與環(huán)境學(xué)院，安徽 淮南 232001；2.安徽理工大學(xué) 空間信息與測繪工程學(xué)院，安徽 淮南 232001)

土壤重金屬的侵襲與累積是一種十分普遍的現(xiàn)象，人類活動將重金屬帶到土壤中，致使土壤中重金屬含量明顯高于背景值，并且造成現(xiàn)存或潛在的土壤質(zhì)量退化、生態(tài)環(huán)境惡化的現(xiàn)象，而礦區(qū)土壤在常年礦產(chǎn)開采加工以及運輸過程中更易受到重金屬污染。鎘(Cd)在自然界中普遍存在，其含量的增加可能會導(dǎo)致鋁、銅、鋅、錳、鉛等元素的高環(huán)境風(fēng)險，通過潛在的地表或地下水的浸出，或通過這些元素更高的生物可給性；Cd的毒性較大，且難以自然降解、清除，容易累積在土壤中；此外重金屬Cd可通過植物根部吸收，對生態(tài)系統(tǒng)及人體健康構(gòu)成威脅，因此，監(jiān)測礦區(qū)土壤中Cd的含量勢在必行。

傳統(tǒng)重金屬檢測大多采用野外采樣-室內(nèi)分析的方法,耗時耗財，無法較好地獲取空間上重金屬的污染情況。近年來高光譜以其動態(tài)、高效等優(yōu)勢已廣泛應(yīng)用于環(huán)境、地質(zhì)、土壤等多個領(lǐng)域，為快速獲取土壤理化信息提供了新的思路。國內(nèi)外學(xué)者在土壤反射率與土壤參數(shù)之間運用多種方法建立了數(shù)學(xué)估算模型。Kemper等對Aznalcollar受污染的土壤利用可見-近紅外光譜實現(xiàn)對6種金屬元素的回歸分析預(yù)測；解憲麗等選擇江西貴溪銅冶煉廠污染區(qū)，分析了9種重金屬元素與可見-近紅外光譜之間的相關(guān)性；龔紹琦等利用光譜一階微分變換，采用逐步回歸方法建立了3種重金屬預(yù)測模型并得到較高精度；Meng Xiangtian等利用14張高光譜圖像采用離散小波變換對原始反射率和光譜一階微分進(jìn)行分解重構(gòu)，采用隨機森林、支持向量機和BP神經(jīng)網(wǎng)絡(luò)算法，建立區(qū)域尺度有機碳預(yù)測模型；Shen Qiang等對湖北大冶鐵礦采用間接反演的方法實現(xiàn)復(fù)墾土壤重金屬Cu的空間分布；陳元鵬等采用偏最小二乘與粒子群算法相結(jié)合的方法,對工礦復(fù)墾區(qū)土壤重金屬反演，主要障礙因子Cd取得較高精度；Hong Yongsheng等利用連續(xù)小波變換結(jié)合隨機森林模型反演土壤有機碳，并得到局地尺度上土壤有機碳空間分布模式；Wei Lifei等采用穩(wěn)定競爭自適應(yīng)加權(quán)采樣算法結(jié)合逐次投影算法解決光譜數(shù)據(jù)冗余，使用重組青蛙跳躍算法優(yōu)化后的徑向基神經(jīng)網(wǎng)絡(luò)模型得到較高預(yù)測精度，以上都說明高光譜對重金屬預(yù)測具有可行性。

前人在土壤重金屬估測方面做了大量的研究并取得了較好的預(yù)測結(jié)果，但多集中于重金屬含量的估測，并未分析模型在不同濃度下的預(yù)測能力。研究針對礦山開采過程中造成或潛在造成的土壤重金屬Cd污染，以粵北南嶺某礦區(qū)為例，通過傳統(tǒng)檢測方法測定研究區(qū)土壤Cd含量；利用ASD Field Spec 4型便攜式高光譜儀進(jìn)行土壤反射率測定，分析不同光譜指標(biāo)與Cd含量的相關(guān)性；篩選敏感波段，采用偏最小二乘和隨機森林建立礦區(qū)土壤重金屬Cd含量估算模型，分析不同Cd含量區(qū)間對建模精度的影響，探討利用高光譜遙感技術(shù)快速檢測土壤Cd含量的可行性，為礦區(qū)土壤重金屬高光譜反演提供方法和理論支持。

1 研究區(qū)概況與研究方法

1.1 研究區(qū)概況

研究區(qū)位于廣東省韶關(guān)市境內(nèi)(112°50′～114°45′E，23°5′～25°31′N)，位置及采樣點如圖1所示。該區(qū)域地處五嶺山脈南麓，北江中上游地區(qū)，全境在大地構(gòu)造上處于華廈活化陸臺的湘粵褶皺帶。韶關(guān)市屬中亞熱帶濕潤型季風(fēng)氣候區(qū)，年平均氣溫18.8～21.6 ℃，年均降雨1 400～2 400 mm，地形以山地丘陵為主，河谷盆地分布其中。平原、臺地面積約占20%，境內(nèi)河流主要屬珠江水系北江流域，以湞江為干流，自北向南貫穿全境。地勢北高南低，土壤類型以紅壤為主。韶關(guān)市作為“中國有色金屬之鄉(xiāng)”，有“中國鋅都”稱號，境內(nèi)礦產(chǎn)資源豐富。該地區(qū)的開采歷史最早可追溯到千年之前，近年來由于大量私人和小集體不合理的開采，使得廢棄后的礦窿成為了地下水的主要污染源，富含重金屬的裸露山體，經(jīng)雨水的沖刷，不斷析出酸性水，最終對土壤造成污染。

圖1 研究區(qū)位置及采樣點圖

1.2 樣品采集與處理

研究土樣采自廣東省韶關(guān)市某金屬礦區(qū)。2020年7月通過GPS精確定位采樣點，共采集樣本點73個。對每個采樣點，按梅花形收集5個子集，再混和成一個組合樣本。將采集的土壤自然風(fēng)干，去除礫石及動植物殘體后，將土樣平均分為兩份，一份用于測定土壤重金屬Cd含量，一份用于采集高光譜數(shù)據(jù)。其中，土壤重金屬Cd含量采用王水提取-電感耦合等離子體質(zhì)譜法(ICP-MS)測定，土壤光譜反射率采集使用美國ASD(Analytical Spectral Devices)公司生產(chǎn)的Field Spec 4便攜式地物光譜儀。該儀器的光譜測量范圍是350～2 500 nm，包含了可見光-近紅外全部范圍，兩次采樣的間隔為1 nm，共有2 150個波段。由于光譜容易受到外界光線的影響，因此實驗選在無光的暗室中進(jìn)行，光源為12 V、50 W的燈泡，數(shù)據(jù)使用25°裸光纖鏡頭獲取。實驗前將儀器先通電預(yù)熱半小時，實驗過程中將土壤樣本均勻平鋪在直徑100 mm、高20 mm的玻璃盛樣皿中，使用黑色絨布為背景墊在盛樣皿下方。光源與樣本的直線距離為60 cm，與水平方向夾角45°。鏡頭位于樣品正上方10 cm，與樣品垂直。為保證數(shù)據(jù)的準(zhǔn)確性，測試之前去除輻射強度中暗電流的影響，然后以白板進(jìn)行定標(biāo)，為防止測量過程中其他因素影響數(shù)據(jù)準(zhǔn)確性，每測量10個樣本進(jìn)行一次白板校正。每個樣本采集20條光譜曲線，剔除噪聲較大的曲線，取平均值作為該土樣的實際反射光譜數(shù)據(jù)。

1.3 光譜預(yù)處理

光譜數(shù)據(jù)獲取過程中，由于外界環(huán)境的影響以及光譜儀在不同波段對能量響應(yīng)上的差異，會導(dǎo)致光譜曲線存在一些噪聲。噪聲主要來自高頻隨機噪聲、基線漂移、樣本不均勻、光線散射等，因此，對光譜數(shù)據(jù)進(jìn)行預(yù)處理就顯得尤為必要。實踐表明，對光譜曲線進(jìn)行平滑，可以去除信號內(nèi)的少量噪聲，得到平滑的光譜波形。如果噪聲的頻率較高且量值不大，用平滑的方法可在一定程度上降低噪聲。研究采用卷積平滑(Savitzky-Golay)方法對光譜曲線進(jìn)行平滑，將得到的結(jié)果作為原始數(shù)據(jù)。

研究表明，通過不同的數(shù)學(xué)變換可以減少土壤母質(zhì)、成土條件、質(zhì)地、表面粗糙度、微聚體、濕度等土壤本身屬性的影響及大氣溫度、濕度、組分和電磁特性等外界因素干擾，有效地消除了基線和背景干擾，提高了部分波段的分辨率和靈敏度，使重疊樣本得到分離，進(jìn)一步提取了原始數(shù)據(jù)中差異不顯著的光譜信息，使得光譜特征更加明顯。研究主要對土壤反射率(Reflectance，R)進(jìn)行以下變換：倒數(shù)變換(Reciprocal Transform，RT)、對數(shù)變換(Logarithm Transform，LT)、一階微分(Frist Derivative，F(xiàn)D)、二階微分(Second Derivative，SD)、歸一化變換(Normalization Transform，NT)、倒數(shù)對數(shù)變換(吸光率，Absorbance Transformr，AT)、倒數(shù)對數(shù)一階微分(ATFD)及倒數(shù)對數(shù)二階微分(ATSD)。

(1)微分技術(shù)。光譜微分技術(shù)是一種在遙感數(shù)據(jù)處理中特別有應(yīng)用前景的分析方法。光譜微分技術(shù)對不同的背景、噪聲有去除作用，特別是比較容易去除以“加”的形式混入光譜信號中的噪聲，還可以消除基線和其他背景的干擾，分辨重疊峰，提高分辨率和靈敏度。一般認(rèn)為，可用一階微分處理去除部分線性或接近線性的背景干擾，二階微分可消除平方項噪聲的影響，因而其在實際應(yīng)用中較為有效，其公式如下：

式中,

λ

為每個波段的波長；

R

(

λ

)和

R

(

λ

-1)分別為波長

λ

和

λ

-1處的光譜反射率；

R

(

λ

)和

R

(

λ

-1)分別為波長

λ

和

λ

-1處的一階微分光譜；

R

(

λ

)為波長

λ

處的二階微分光譜；

Δλ

為波長

λ

-到

λ

的間隔，視波段波長而定，波長

λ

-2到

λ

的間隔及波長

λ

+2到

λ

的間隔為2

Δλ

。

(2)初等變換。采用初等函數(shù)對光譜數(shù)據(jù)進(jìn)行變化處理叫做初等變換，也稱簡單變換。實際應(yīng)用表明，初等變換一般不能有效提高變換后光譜數(shù)據(jù)與研究對象間的相關(guān)性，但為研究組合變換方法提供了基礎(chǔ)，可根據(jù)具體問題，通過對比試驗獲取最佳組合變換方法。

1.4 模型建立與精度檢驗

經(jīng)過以上數(shù)據(jù)預(yù)處理，從73個土壤樣本中隨機選擇52個作為建模樣本，21個作為檢驗樣本用來檢驗?zāi)Ｐ途取?/p>

(1)偏最小二乘算法。偏最小二乘(Partial Least Squares Regression，PLSR)是一種基于主成分分析的多變量建模方法，通過將光譜數(shù)據(jù)進(jìn)行分解，去除其中無效的噪聲干擾，同時分解重金屬含量數(shù)據(jù)，消除其中無用信息。在分解光譜數(shù)據(jù)的同時考慮了重金屬含量數(shù)據(jù)的影響，將數(shù)據(jù)分解與回歸并為一步。在計算每一個主成分前將光譜數(shù)據(jù)的得分矩陣與重金屬含量數(shù)據(jù)的得分矩陣交換，使得光譜數(shù)據(jù)主成分直接與重金屬含量數(shù)據(jù)關(guān)聯(lián)。

(2)隨機森林算法。隨機森林(Random Forest Regression，RF)是Breiman 2001年提出的一種利用多棵決策樹對樣本進(jìn)行訓(xùn)練并預(yù)測的機器學(xué)習(xí)算法，是眾多決策樹的集合。采用Bootsrap重抽樣方法隨機抽樣構(gòu)建不同的分類模型，再用它們構(gòu)成一個多分類模型系統(tǒng)，以多數(shù)投票法確定最終分類結(jié)果。對異常值和噪聲具有較高的容忍度，通過引入“隨機性”來處理“過擬合”，因此具有較高的泛化能力。

(3)精度檢驗。模型的驗證主要選用決定系數(shù)(Coefficient of Determination，R)和均方根誤差(Root Mean Squard Error，RMSE)，計算公式如下：

2 結(jié)果與分析

2.1 土壤Cd的經(jīng)典統(tǒng)計學(xué)分析

對本次測定的73個土樣的重金屬Cd含量進(jìn)行統(tǒng)計如表1所示。將檢測結(jié)果與國家相關(guān)標(biāo)準(zhǔn)及廣東省背景值對比，發(fā)現(xiàn)超過72.6%的土壤樣點Cd含量高于農(nóng)用地土壤污染風(fēng)險管控標(biāo)準(zhǔn)，含量最大值點位超過背景值145倍。研究區(qū)Cd含量最大值達(dá)到4.95 mg/kg，最小值為0.04 mg/kg，平均值為1.11 mg/kg。從土樣的平均值來看，Cd含量超過背景值30倍。根據(jù)單因子指數(shù)法測算，該區(qū)域Cd的污染指數(shù)為重度污染。Cd含量變異系數(shù)為101.14%，通常認(rèn)定變異系數(shù)反應(yīng)離散程度，且當(dāng)100%<變異系數(shù)時，為強變異性，可能是由于采樣點受采礦區(qū)擾動程度不同，土壤Cd含量差異較大。因此，該地區(qū)應(yīng)加強土壤重金屬動態(tài)監(jiān)測，為發(fā)現(xiàn)并控制土壤重金屬污染提供依據(jù)。

表1 土壤樣本Cd含量描述性統(tǒng)計

2.2 土壤光譜特征分析

圖2 不同Cd含量土壤光譜曲線和連續(xù)統(tǒng)去除光譜曲線

2.3 土壤Cd與光譜反射率的相關(guān)性分析

Cd含量與不同數(shù)學(xué)變換后的相關(guān)系數(shù)圖如圖3所示。由圖3可知，由于土壤中Cd含量較低，Cd含量與土壤光譜總體相關(guān)性較低，但都表現(xiàn)出與土壤光譜相似的相關(guān)性，分別在700～800 nm、1 380～1 410 nm、1 550～1 600 nm、1 800～1 950 nm、2 320～2 360 nm范圍內(nèi)有較強的相關(guān)性。Cd含量與原始光譜絕大部分呈負(fù)相關(guān)，在可見光部分相關(guān)性先降低再升高，在900 nm附近逐漸趨于穩(wěn)定，在1 900 nm附近相關(guān)性降低后逐漸升高，并在2 212 nm附近達(dá)到顯著相關(guān)。LT、RT相關(guān)性與R相似，均在400～560 nm范圍內(nèi)降低后逐漸升高并趨于穩(wěn)定。R、RT、AT、NT相關(guān)系數(shù)最大值均出現(xiàn)在可見光波段，說明在進(jìn)行微分變換之前，在可見光波段的土壤Cd含量探測能力要比近紅外波段強。經(jīng)微分變換之后的4種光譜指標(biāo)與土壤Cd含量相關(guān)性明顯提高，相關(guān)性在正負(fù)之間波動很大，且近紅外波段的相關(guān)性高于可見光波段，部分波段達(dá)到了0.01的極顯著水平。FD與Cd含量相關(guān)性最高的波段為1 943 nm(

r

=-0

.

496)，SD與Cd含量相關(guān)性最高值略低于FD但總體相差不大，最高的波段為1 945 nm(

r

=-0

.

467)。ATFD與ATSD的相關(guān)性沒有FD和SD那么突出，變化趨勢與FD和SD相似，均在正負(fù)之間波動，相關(guān)性最高的波段分別為758 nm(

r

=0

.

466)和2 387nm(

r

=-0

.

487)。綜上所述，在9種微分指標(biāo)中一階微分變化與Cd含量相關(guān)性最好，為最優(yōu)光譜指標(biāo)。

圖3 Cd含量與土壤光譜的相關(guān)系數(shù)圖

2.4 土壤Cd含量高光譜模型分析

選擇與土壤Cd含量相關(guān)性較大的波段作為自變量，采用偏最小二乘模型和隨機森林方法，分別建立9種光譜指標(biāo)(R、RT、LT、FD、SD、NT、AT、ATFD、ATSD)的礦區(qū)土壤重金屬預(yù)測模型，模型的預(yù)測能力和穩(wěn)定性由

R

和

RMSE

檢驗，結(jié)果如表2所示。

表2 不同土壤反射率變換形式模型回歸結(jié)果

從模型的回歸效果看，兩種模型的預(yù)測能力差異較大，偏最小二乘的

R

介于0

.

14～0

.

61，

RMSE

介于0

.

69～1

.

11；隨機森林的

R

介于0

.

13～0

.

86，

RMSE

介于0

.

40～1

.

09。不同預(yù)處理變換方法對模型的預(yù)測結(jié)果有較大的影響，與相關(guān)性結(jié)果類似，經(jīng)過初等變換的光譜曲線預(yù)測能力較差，所建立模型的

R

均小于0

.

50，只能較為粗略地估計Cd含量的高低，無法準(zhǔn)確預(yù)測含量值；經(jīng)過微分變換的光譜曲線預(yù)測能力有較高的提升，所建立模型的

R

最大值超過0

.

80，能夠較為準(zhǔn)確地預(yù)測土壤Cd含量。擬合效果最好的是FD-RF所建立的模型，

R

為0

.

85，

RMSE

為0

.

40，NT變換效果最差?；贔D-RF法建立的高光譜模型的預(yù)測值與實際值的散點圖如圖4所示。由圖4可見，樣本點與1∶1的線很近，說明用此方法建模預(yù)測能力高，對土壤Cd含量具有較好的解釋能力。

圖4 基于FD-RF法土壤Cd含量實測值與預(yù)測值散點圖

2.5 模型適應(yīng)性分析

為探究Cd濃度對模型預(yù)測能力的影響，以最優(yōu)反演模型為建模手段，建立不同濃度區(qū)間的Cd含量預(yù)測模型。將73個樣本分為7個濃度區(qū)間(Cd～Cd)，相鄰兩個區(qū)間的濃度平均值相差約0.20 mg/kg，為排除樣本個數(shù)與Cd濃度值大小對模型精度評價的影響，每個濃度區(qū)間選擇40個樣本，將同一濃度區(qū)間的樣本按照8:2的比例隨機分成建模集與驗證集，同時以

R

與

RMSE

作為最終的評價指標(biāo)，研究結(jié)果如表3所示。

表3 不同Cd濃度區(qū)間回歸模型統(tǒng)計

從模型的回歸效果看，模型對不同Cd濃度的預(yù)測能力不同，其中預(yù)測能力較高的濃度區(qū)間為Cd～Cd，平均值為0.80～1.20 mg/kg，

R

均超過0

.

80，

RMSE

均小于0.10 mg/kg；預(yù)測能力較差的濃度區(qū)間為Cd，平均值為1.60 mg/kg，

R

=0

.

64，只能粗略地估計Cd含量值高低。綜上所述，當(dāng)樣本濃度平均值變化小于20%時，模型的預(yù)測能力較為穩(wěn)定；當(dāng)樣本濃度平均值變化介于20%～40%時，模型的預(yù)測能力下降，當(dāng)超過40%時，只能粗略估計含量高低，無法準(zhǔn)確預(yù)測。

3 討論

礦區(qū)土壤環(huán)境復(fù)雜，常年的煤炭開采活動對土壤重金屬的空間分布產(chǎn)生較大的影響，不同的土地利用方式使得土壤重金屬Cd在一定的區(qū)域內(nèi)發(fā)生變異性較大的現(xiàn)象，過高的Cd含量可能對人體造成潛在的危險，因此快速監(jiān)測土壤Cd含量對于礦區(qū)土壤治理具有重要意義。高光譜遙感近年來以其快速、無損、低成本等特點已廣泛應(yīng)用于土壤屬性預(yù)測。研究以廣東省韶關(guān)市某礦區(qū)為研究對象，利用偏最小二乘與隨機森林算法實現(xiàn)土壤Cd含量預(yù)測，并對模型進(jìn)行了適應(yīng)性分析。

由于Cd在土壤中含量較低，導(dǎo)致土壤光譜反射率與Cd含量的相關(guān)性較低。經(jīng)過不同的光譜變換可以有效放大波光譜中差異不顯著的信息，并去除外界因素引起的噪聲，進(jìn)而提升土壤光譜反射率與Cd含量的相關(guān)性，但總體上均沒有超過0.5。本研究中，與偏最小二乘相比，隨機森林具有更好的預(yù)測能力，這可能是由于光譜信息與土壤Cd含量之間不僅存在著線性關(guān)系，還存在非線性關(guān)系，這與陳亦凡的相關(guān)研究結(jié)果相近。土壤光譜特征是土壤系統(tǒng)與外部環(huán)境交互的綜合反映，不同的成土母質(zhì)、粒徑大小、有機質(zhì)含量高低都會影響土壤光譜曲線的形狀。在實際模型驗證中發(fā)現(xiàn)，某些含量較低的樣本的預(yù)測濃度存在虛高現(xiàn)象，這可能是由于Cd含量較低，其光譜特征被其他土壤成分掩蓋，這與彭杰等的研究結(jié)果較為一致。通過對比Cd含量在不同濃度下的預(yù)測能力可以發(fā)現(xiàn)，所建立的模型在不同濃度區(qū)間的預(yù)測能力不同。在Cd～Cd區(qū)間范圍內(nèi)的預(yù)測能力較高，在Cd的預(yù)測能力最低，與全樣本相比，當(dāng)樣本濃度平均值變化小于20%時，模型的預(yù)測能力較為穩(wěn)定；當(dāng)樣本濃度平均值變化超過40%時，模型的參數(shù)需要重新調(diào)整。這可能是由于當(dāng)Cd濃度變化過大時，土壤理化性質(zhì)差異性也變大，影響模型預(yù)測精度。近年來，隨著不同國家與地區(qū)建立的不同尺度的土壤高光譜數(shù)據(jù)庫，積累了大量的土壤光譜基礎(chǔ)數(shù)據(jù)，這為研究提高模型預(yù)測適用性提供了數(shù)據(jù)基礎(chǔ)。研究通過偏最小二乘和隨機森林算法實現(xiàn)了礦區(qū)土壤Cd含量預(yù)測，并將樣本分成不同的濃度區(qū)間，探討了濃度變化對模型的影響。該思路為今后檢測土壤屬性參數(shù)提供了理論支持。

4 結(jié)論

基于本研究實驗數(shù)據(jù)，采用偏最小二乘和隨機森林方法建立廣東某礦區(qū)的土壤重金屬Cd含量的高光譜反演模型，重點研究了不同Cd含量下光譜曲線差異以及不同光譜處理方法對建模精度的差異，探究了不同Cd濃度區(qū)間對建模精度的影響，明確了運用微分處理方法可以有效提升模型預(yù)測精度。結(jié)果表明:

(1)Cd含量與光譜反射率呈負(fù)相關(guān)，Cd含量增加會導(dǎo)致反射率下降，但不會影響特征波段與吸收峰的位置。通過對原始光譜數(shù)據(jù)進(jìn)行不同方式的處理，可以有效地去除噪聲的影響，提高相關(guān)性。其中一階微分效果最佳，在1 406 nm處達(dá)到最大值，相關(guān)系數(shù)

r

=-0

.

503。(2)基于不同回歸方法建立回歸模型，預(yù)測效果最佳的模型為FD-RF，模型的預(yù)測值與實測值的

R

=0

.

83、

RMSE

=0

.

40，對于理化性質(zhì)差異較大的土壤，非線性模型的預(yù)測能力優(yōu)于線性模型。當(dāng)樣本濃度平均值變化小于20%時，模型的預(yù)測能力較為穩(wěn)定；當(dāng)樣本濃度平均值變化介于20%～40%時，模型的預(yù)測能力下降，當(dāng)樣本濃度平均值變化超過40%時，只能粗略估計含量高低，無法準(zhǔn)確預(yù)測。