郭云開 ,謝曉峰,謝 瓊,張思愛
(1.長沙理工大學 交通運輸工程學院,湖南 長沙 410076;2.長沙理工大學 測繪遙感應用技術研究所,湖南 長沙 410076;3.湖南工程職業(yè)技術學院 測繪地理學院,湖南 長沙 410151;4.長沙理工大學 公路地質(zhì)災害預警空間信息技術湖南省工程實驗室,湖南 長沙 410114)
耕地是人類生產(chǎn)生活的基本元素。隨著工業(yè)化的快速發(fā)展,耕地土壤面臨的重金屬污染問題日益嚴重。在土壤中,重金屬鉻(Cr)主要以Cr3+、Cr6+的形式存在。Cr6+屬于強致突變物質(zhì),可誘發(fā)肺癌和鼻咽癌。Cr3+主要以氫氧化物的形式附著于土壤組分中,不易遷移且毒性較低[1-2]。傳統(tǒng)監(jiān)測土壤重金屬含量的方法主要是以點代面采集樣本,帶回室內(nèi)進行化學分析。這種方法雖然精度高,但是費時費力且無法在大區(qū)域內(nèi)進行實時監(jiān)測。隨著高光譜遙感技術的發(fā)展,高光譜遙感技術以其快速獲取地面信息的能力和光譜分辨率高的特點為土壤重金屬污染監(jiān)測提供可能[4-5]。目前,在土壤重金屬污染監(jiān)測領域,很多專家已經(jīng)采用高光譜遙感技術展開研究。文獻[6]在對原始反射率光譜進行倒數(shù)對數(shù)、一階微分、二階微分和包絡線去除處理等光譜變換,發(fā)現(xiàn)基于二階微分處理后的數(shù)據(jù)建立的偏最小二乘模型(PLS)的預測精度高于其他方法。袁中強[7]等對比原始地表反射率、原始地表反射率的一階微分和倒數(shù)對數(shù)等光譜變換方法,發(fā)現(xiàn)基于原始地表反射率的一階微分值建立的回歸模型對重金屬鉻的預測效果最好。龔紹琦[8]等在對土壤光譜反射率進行一階微分、倒數(shù)對數(shù)和連續(xù)統(tǒng)去除等預處理的基礎上建立逐步回歸模型預測重金屬鉻含量,得出一階微分為鉻含量反演的最佳光譜變換。
通過微分處理等傳統(tǒng)光譜變換方法,提取特征波段,建立預測模型,在一定程度上可以有效預測土壤中重金屬鉻含量。但是土壤中重金屬的含量屬于微量級,其光譜信息十分微弱。此外,微分處理反射光譜的方法在增強微弱信號的同時也會擴大噪聲,降低光譜的信噪比[9]。在對微弱光譜信號的處理中,小波變換是一種較為理想的提升信噪比的方法[10]。因此,本研究采用Db4小波對原始光譜進行分解與重構,提取特征光譜,利用特征光譜與重金屬含量的相關性建立PLS模型預測重金屬鉻的含量。
研究對象選取湖南省岳陽縣某地區(qū)(見圖1)。該地區(qū)地勢平坦,土地肥沃,是湖南省糧食主產(chǎn)區(qū)。此外,該地區(qū)有很多中小企業(yè)。隨著企業(yè)的快速發(fā)展,工業(yè)廢水、廢渣存在處理不當?shù)默F(xiàn)象,這使得該地區(qū)面臨耕地污染等環(huán)境挑戰(zhàn)。
圖1 岳陽縣某地區(qū)衛(wèi)星圖
根據(jù)農(nóng)田分布情況布設采樣點,每個樣本點采樣范圍為10 m×10 m。在采樣范圍內(nèi),采用5點采樣法采集0~20 cm范圍內(nèi)的表層土壤,取約1 kg土壤裝袋保存。將采集回來的土壤去除石塊、雜草等明顯雜質(zhì),放置在陰涼的通風室進行風干處理。風干后的土壤進行磨碎、過100目尼龍篩后將樣本均勻分為兩份。取其中一份進行重金屬含量的測定。重金屬含量測定采用電熱板消化法,通過電感耦合等離子發(fā)射光譜儀測定。
使用Field-Spec AvaField-3高精度地物波普儀采集室內(nèi)波譜。該儀器波段范圍為:300~2 500 nm。室內(nèi)光譜采集過程中,以聚四氟乙烯白板為標準白板,以1 000 W鹵素燈作為光源,控制視場角為5°,光源與樣本距離約30 cm,照射方向與垂直方向成15°。光譜采集前開機預熱10 min,每個樣品重復采集多條光譜曲線,取平均作為該樣品的原始光譜。
樣本采集總量為40個,剔除兩個光譜異常的樣本。此外,為保證數(shù)據(jù)精度,本研究去除了光譜中信噪比低的340~380 nm和2 401~2 515 nm光譜[11]。由于原始光譜分辨率高,波段之間存在大量冗余,數(shù)據(jù)量龐大。為降低數(shù)據(jù)量,在保證能充分反映土壤光譜特性的前提下,對原始光譜進行6 nm光譜重采樣。
小波分解與重構是能夠同時在時域和頻域內(nèi)進行局部化信號分析的數(shù)學方法。它具有多分辨率分析的特點,在時、頻兩域都具有表征信號細節(jié)特征的能力,適合分析非平穩(wěn)信號和提取信號的局部特征[12-13]。土壤的原始光譜信號屬于非平穩(wěn)信號,適合采用小波變換提取光譜特征。
小波分解是將光譜信號分解為不同的層,各層所包含的信息與土壤不同成分有關。低頻系數(shù)反映決定土壤光譜基本形狀的物質(zhì),高頻系數(shù)反映光譜采集過程的各種噪聲及土壤微量成分在原始光譜的特征[14]。土壤中重金屬含量屬于微量元素,主要體現(xiàn)在高頻系數(shù)中。小波重構是把分解的系數(shù)還原的過程,是小波分解的逆過程。小波單層重構是基于小波分解后單層小波系數(shù)進行重構,獲得特定頻率的重構信號[15]。
目前,小波基的選擇沒有統(tǒng)一的標準,需要結合信號本身的特點和小波函數(shù)的數(shù)學特性綜合考慮。本研究選擇的Db4小波基函數(shù)具有正交性、雙正交性、緊支性和不對稱性的特點。同時,Db4小波基能夠解決線性相位與正交性要求的矛盾,提取出與重金屬含量具有較高相關性的局部高光譜信息[16]。
分解尺度的確立是小波分解的關鍵。分解尺度可以根據(jù)母小波函數(shù)和光譜信號長度共同確定[17]。確定小波分解的最佳分解尺度n的方程為:
n=1b(N).
(1)
其中N為高光譜信號的長度。
此外,最佳分解尺度n也可根據(jù)母小波的濾波器長度來計算確定[17]。本研究以Db4作為小波基,根據(jù)母小波濾波器長度計算得最佳分解尺度為4層。但是該理論計算得到的最佳分解尺度只是一個概率估計,可能存在偏差。因此,綜合小波分解的特點及本研究數(shù)據(jù)特點,確定小波分解尺度為5層。
以波段與重金屬含量的顯著相關程度作為判斷特征波段的依據(jù)。在軟件Matlab中編寫程序實現(xiàn)小波分解和重構,使用重構后的信號與重金屬含量進行相關性分析提取特征波段,提取的主要過程為:
1)采用Db4小波基對原始土壤光譜進行5層離散分解,獲得5個高頻系數(shù)和1個低頻系數(shù)。
2)在Matlab中采用Wrcoef函數(shù)分別對高頻系數(shù)和低頻系數(shù)進行單層系數(shù)重構,獲得不同的重構信號。
3)匯總各層重構信號,將1~5層高頻信息分別命名為B1、B2、B3、B4、B5。由于屬于痕量級的土壤重金屬主要體現(xiàn)在高頻系數(shù)中,同時為充分利用高頻信息,研究將與重金屬含量存在顯著相關的各層高頻信號進行相加,將得到的數(shù)據(jù)命名為Bn~Bm加和。
4)將得到的各組數(shù)據(jù)分別與重金屬鉻含量進行皮爾森相關性分析,通過0.05雙側顯著性檢驗的波段作為特征波段。
如前述,有效樣本總量為38個,隨機選取27個樣本作為建模集,11個樣本作為驗證集。將B1~B5各層重構信號、B1~B4加和、一階微分處理后的土壤光譜分別與重金屬鉻含量進行皮爾森相關性分析,結果如圖2—圖8所示。
圖2 B1層波段與Cr相關系數(shù)
圖3 B2層波段與Cr相關系數(shù)
圖4 B3層波段與Cr相關系數(shù)
圖5 B4層波段與Cr相關系數(shù)
圖6 B5層波段與Cr相關系數(shù)
圖7 B6層波段與Cr相關系數(shù)
圖8 一階微分處理后全波段與Cr相關系數(shù)
統(tǒng)計各組數(shù)據(jù)的特征波段個數(shù),結果如表1所示。單層高頻系數(shù)重構信號中,B4層重構信號中特征波段的個數(shù)最多,達到22個。其次是B3層重構信號,特征波段的個數(shù)為13個。B5層重構信號不存在特征波段,說明在第5層重構信號中沒有重金屬鉻的信息,可以視為噪聲去除?;谝浑A微分處理后的土壤光譜提取特征波段,共計11個。統(tǒng)計分析B1~B4各層重構信號與重金屬鉻含量都存在顯著相關,因此選擇B1~B4層的重構信號進行加和運算,獲得B1~B4層加和數(shù)據(jù)?;贐1~B4層加和數(shù)據(jù)提取特征波段,結果存在24個顯著相關波段。
表1 重金屬含量與各層重構信號存在相關的波段個數(shù)
選擇B1~B4層加和、B4層及一階微分處理后數(shù)據(jù)分別建立多元線性模型。建模結果如表2所示。
表2 不同的小波重構信號建立多元線性模型的結果
基于B1~B4層加和建立的多元線性預測模型的決定系數(shù)為0.916?;贐4層重構信號建立的多元線性模型的決定系數(shù)為0.339?;谝浑A微分處理光譜信息建立的多元線性模型的決定系數(shù)為0.521。結合表2發(fā)現(xiàn),B4層重構信號的特征波段的個數(shù)大于一階微分處理后特征波段的個數(shù),而基于B4層重構信號建立的多元線性模型的決定系數(shù)卻小于基于一階微分處理后的光譜信息所建立模型的決定系數(shù)。出現(xiàn)這樣的結果,可能是因為作為自變量的特征波段之間存在嚴重的共線性問題。
方差膨脹因子(VIF)是衡量多元線性回歸模型中共線性嚴重程度的一種度量。一般情況下方差膨脹因子大于10,表明模型中自變量之間存在很強的共線性問題[18]。分別對B1~B4層加和、B4層、一階微分處理土壤光譜提取的特征光譜進行共線性診斷,結果如表3所示。分析表3可得,B4層、B1~B4層加和的特征波段,方差膨脹因子(VIF)遠大于10,基于一階微分處理提取的特征波段存在多個波段的VIF大于10。因此,基于B4層、B1~B4層加和及一階微分處理后光譜數(shù)據(jù)提取的特征波段之間都存在共線性問題,其中,基于B4層、B1~B4層加和提取的特征波段之間共線性尤為嚴重。
表3 方差膨脹因子匯總
偏最小二乘回歸模型是一種多因變量對多自變量的回歸建模方法,同時包含多元線性回歸分析、典型相關分析和主成分分析的特點,能夠較好的解決變量之間存在多重相關性或樣本點數(shù)過少的回歸難題的方法[19-20]。為消除特征波段之間的共線性影響,建立可靠的預測模型,選擇PLS作為建模方法。基于B1~B4層加和、B4層和一階微分處理后的三組數(shù)據(jù)分別建立PLS預測模型,建模結果如表4所示。
分析表4可知,基于一階導數(shù)處理后的數(shù)據(jù)建立PLS預測模型,在8個潛在因子、自變量累積解釋能力為99.4%的前提下,PLS模型的決定系數(shù)R2為0.519?;贐4層重構信號進行PLS建模,在潛在因子為7個、自變量累積解釋能力100%的前提下,模型的決定系數(shù)R2為0.481。基于B1~B4層加和建模結果為:在潛在因子為9個、自變量累積解釋能力99.5%的前提下,模型的決定系數(shù)R2為0.605?;贐1~B4層加和建模結果比基于B4層重構信號建模的R2提升了0.124,比基于一階導數(shù)處理后的數(shù)據(jù)建模R2提升0.086。結果說明,直接在原始光譜的基礎上進行小波分解與重構,提取特征波段,利用特征波段與重金屬含量之間的相關性建立PLS模型,可以有效預測重金屬鉻的含量。同時,與常用的一階微分變換的方法相比,基于小波分解與重構提取的特征波段建立PLS模型的預測精度有所提升。
表4 不同小波重構信號建立PLS模型結果
根據(jù)建模集所建立的PLS模型系數(shù),得到采用一階微分光譜變換所提取的特征波段與重金屬鉻含量相關性建立的PLS模型為:
y=-635.527×χ1+3 546.074×χ2-
1 147.016×χ3-52.217×χ4+259.237×χ5-
3 093.044×χ6+1 036.373×χ7+
3 599.764×χ8-1 426.175×χ9+
1 923.961×χ10-41.321×χ11+11.406.
(2)
利用B1~B4層加和提取的特征波段與重金屬鉻含量的相關性建立的PLS模型為:
y=-654.654×χ1+646.323×χ2+
2 615.525×χ3-3 940.077×χ4+
4 670.818×χ5-1 957.401×χ6-
1 797.477×χ7-2 956.227×χ8+
3 810.497×χ9+1 645.746×χ10-
220.139×χ11+626.171×χ12-
1 107.56×χ13-982.036×χ14-
3 342.157×χ15-921.935×χ16
4 046.777×χ17+465.223×χ18-
3 928.37×χ19+1 141.567×χ20-
5 354.198×χ21-710.444×χ22+
3 132.64×χ23-1 518.24×χ24+10.416.
(3)
基于驗證集數(shù)據(jù)分別驗證以上兩個模型,計算結果如表5所示。
分析表5可知,在驗證集中,基于B1~B4層加和建立PLS模型的決定系數(shù)R2為0.511,比基于一階微分建立的PLS模型的決定系數(shù)(R2)0.351提高0.16。結果表明:利用小波變換處理原始光譜建立偏最小二乘模型具有較高的精度。
表5 預測集建模精度匯總
本研究采用基于土壤重金屬原始光譜進行小波分解與重構的方法,提取特征波段,建立PLS預測建模,對土壤重金屬鉻的含量進行定量反演,并與基于常用的一階微分處理土壤光譜建立PLS模型進行比較。結果表明:
1)利用小波分解與重構提取的有用信息可以有效預測土壤重金屬鉻的含量;
2)對比基于常用的一階微分處理提取的特征波段和基于小波變換處理所提取的特征波段,發(fā)現(xiàn)利用小波變換處理所提取的特征波段與重金屬鉻的含量相關性建立的PLS模型具有較高精度。
雖然小波分解與重構的方法可以有效提取光譜中有用信息,但本研究是基于室內(nèi)光譜開展的,對于復雜的室外環(huán)境,小波分解與重構提取特征波段的方法是否適用仍需在今后研究中進行探討。