基于小波變換的土壤重金屬鉻含量定量反演

郭云開(kāi) ,謝曉峰,謝 瓊,張思愛(ài)

(1．長(zhǎng)沙理工大學(xué) 交通運(yùn)輸工程學(xué)院，湖南 長(zhǎng)沙 410076；2．長(zhǎng)沙理工大學(xué) 測(cè)繪遙感應(yīng)用技術(shù)研究所，湖南 長(zhǎng)沙 410076；3．湖南工程職業(yè)技術(shù)學(xué)院 測(cè)繪地理學(xué)院，湖南 長(zhǎng)沙 410151;4．長(zhǎng)沙理工大學(xué) 公路地質(zhì)災(zāi)害預(yù)警空間信息技術(shù)湖南省工程實(shí)驗(yàn)室，湖南 長(zhǎng)沙 410114)

耕地是人類生產(chǎn)生活的基本元素。隨著工業(yè)化的快速發(fā)展，耕地土壤面臨的重金屬污染問(wèn)題日益嚴(yán)重。在土壤中，重金屬鉻(Cr)主要以Cr3+、Cr6+的形式存在。Cr6+屬于強(qiáng)致突變物質(zhì)，可誘發(fā)肺癌和鼻咽癌。Cr3+主要以氫氧化物的形式附著于土壤組分中，不易遷移且毒性較低[1-2]。傳統(tǒng)監(jiān)測(cè)土壤重金屬含量的方法主要是以點(diǎn)代面采集樣本，帶回室內(nèi)進(jìn)行化學(xué)分析。這種方法雖然精度高，但是費(fèi)時(shí)費(fèi)力且無(wú)法在大區(qū)域內(nèi)進(jìn)行實(shí)時(shí)監(jiān)測(cè)。隨著高光譜遙感技術(shù)的發(fā)展，高光譜遙感技術(shù)以其快速獲取地面信息的能力和光譜分辨率高的特點(diǎn)為土壤重金屬污染監(jiān)測(cè)提供可能[4-5]。目前，在土壤重金屬污染監(jiān)測(cè)領(lǐng)域，很多專家已經(jīng)采用高光譜遙感技術(shù)展開(kāi)研究。文獻(xiàn)[6]在對(duì)原始反射率光譜進(jìn)行倒數(shù)對(duì)數(shù)、一階微分、二階微分和包絡(luò)線去除處理等光譜變換，發(fā)現(xiàn)基于二階微分處理后的數(shù)據(jù)建立的偏最小二乘模型(PLS)的預(yù)測(cè)精度高于其他方法。袁中強(qiáng)[7]等對(duì)比原始地表反射率、原始地表反射率的一階微分和倒數(shù)對(duì)數(shù)等光譜變換方法，發(fā)現(xiàn)基于原始地表反射率的一階微分值建立的回歸模型對(duì)重金屬鉻的預(yù)測(cè)效果最好。龔紹琦[8]等在對(duì)土壤光譜反射率進(jìn)行一階微分、倒數(shù)對(duì)數(shù)和連續(xù)統(tǒng)去除等預(yù)處理的基礎(chǔ)上建立逐步回歸模型預(yù)測(cè)重金屬鉻含量，得出一階微分為鉻含量反演的最佳光譜變換。

通過(guò)微分處理等傳統(tǒng)光譜變換方法，提取特征波段，建立預(yù)測(cè)模型，在一定程度上可以有效預(yù)測(cè)土壤中重金屬鉻含量。但是土壤中重金屬的含量屬于微量級(jí)，其光譜信息十分微弱。此外，微分處理反射光譜的方法在增強(qiáng)微弱信號(hào)的同時(shí)也會(huì)擴(kuò)大噪聲，降低光譜的信噪比[9]。在對(duì)微弱光譜信號(hào)的處理中，小波變換是一種較為理想的提升信噪比的方法[10]。因此，本研究采用Db4小波對(duì)原始光譜進(jìn)行分解與重構(gòu)，提取特征光譜，利用特征光譜與重金屬含量的相關(guān)性建立PLS模型預(yù)測(cè)重金屬鉻的含量。

1 數(shù)據(jù)的獲取與處理

1.1 研究區(qū)域概況

研究對(duì)象選取湖南省岳陽(yáng)縣某地區(qū)(見(jiàn)圖1)。該地區(qū)地勢(shì)平坦，土地肥沃，是湖南省糧食主產(chǎn)區(qū)。此外，該地區(qū)有很多中小企業(yè)。隨著企業(yè)的快速發(fā)展，工業(yè)廢水、廢渣存在處理不當(dāng)?shù)默F(xiàn)象，這使得該地區(qū)面臨耕地污染等環(huán)境挑戰(zhàn)。

圖1 岳陽(yáng)縣某地區(qū)衛(wèi)星圖

1.2 數(shù)據(jù)采集及預(yù)處理

根據(jù)農(nóng)田分布情況布設(shè)采樣點(diǎn)，每個(gè)樣本點(diǎn)采樣范圍為10 m×10 m。在采樣范圍內(nèi)，采用5點(diǎn)采樣法采集0～20 cm范圍內(nèi)的表層土壤，取約1 kg土壤裝袋保存。將采集回來(lái)的土壤去除石塊、雜草等明顯雜質(zhì)，放置在陰涼的通風(fēng)室進(jìn)行風(fēng)干處理。風(fēng)干后的土壤進(jìn)行磨碎、過(guò)100目尼龍篩后將樣本均勻分為兩份。取其中一份進(jìn)行重金屬含量的測(cè)定。重金屬含量測(cè)定采用電熱板消化法，通過(guò)電感耦合等離子發(fā)射光譜儀測(cè)定。

使用Field-Spec AvaField-3高精度地物波普儀采集室內(nèi)波譜。該儀器波段范圍為：300～2 500 nm。室內(nèi)光譜采集過(guò)程中，以聚四氟乙烯白板為標(biāo)準(zhǔn)白板，以1 000 W鹵素?zé)糇鳛楣庠矗刂埔晥?chǎng)角為5°，光源與樣本距離約30 cm，照射方向與垂直方向成15°。光譜采集前開(kāi)機(jī)預(yù)熱10 min，每個(gè)樣品重復(fù)采集多條光譜曲線，取平均作為該樣品的原始光譜。

樣本采集總量為40個(gè)，剔除兩個(gè)光譜異常的樣本。此外，為保證數(shù)據(jù)精度，本研究去除了光譜中信噪比低的340～380 nm和2 401～2 515 nm光譜[11]。由于原始光譜分辨率高，波段之間存在大量冗余，數(shù)據(jù)量龐大。為降低數(shù)據(jù)量，在保證能充分反映土壤光譜特性的前提下，對(duì)原始光譜進(jìn)行6 nm光譜重采樣。

2 小波分解與重構(gòu)

小波分解與重構(gòu)是能夠同時(shí)在時(shí)域和頻域內(nèi)進(jìn)行局部化信號(hào)分析的數(shù)學(xué)方法。它具有多分辨率分析的特點(diǎn)，在時(shí)、頻兩域都具有表征信號(hào)細(xì)節(jié)特征的能力，適合分析非平穩(wěn)信號(hào)和提取信號(hào)的局部特征[12-13]。土壤的原始光譜信號(hào)屬于非平穩(wěn)信號(hào)，適合采用小波變換提取光譜特征。

小波分解是將光譜信號(hào)分解為不同的層，各層所包含的信息與土壤不同成分有關(guān)。低頻系數(shù)反映決定土壤光譜基本形狀的物質(zhì)，高頻系數(shù)反映光譜采集過(guò)程的各種噪聲及土壤微量成分在原始光譜的特征[14]。土壤中重金屬含量屬于微量元素，主要體現(xiàn)在高頻系數(shù)中。小波重構(gòu)是把分解的系數(shù)還原的過(guò)程，是小波分解的逆過(guò)程。小波單層重構(gòu)是基于小波分解后單層小波系數(shù)進(jìn)行重構(gòu)，獲得特定頻率的重構(gòu)信號(hào)[15]。

2.1 小波基選擇及分解尺度的確定

目前，小波基的選擇沒(méi)有統(tǒng)一的標(biāo)準(zhǔn)，需要結(jié)合信號(hào)本身的特點(diǎn)和小波函數(shù)的數(shù)學(xué)特性綜合考慮。本研究選擇的Db4小波基函數(shù)具有正交性、雙正交性、緊支性和不對(duì)稱性的特點(diǎn)。同時(shí)，Db4小波基能夠解決線性相位與正交性要求的矛盾，提取出與重金屬含量具有較高相關(guān)性的局部高光譜信息[16]。

分解尺度的確立是小波分解的關(guān)鍵。分解尺度可以根據(jù)母小波函數(shù)和光譜信號(hào)長(zhǎng)度共同確定[17]。確定小波分解的最佳分解尺度n的方程為：

n=1b(N).

(1)

其中N為高光譜信號(hào)的長(zhǎng)度。

此外，最佳分解尺度n也可根據(jù)母小波的濾波器長(zhǎng)度來(lái)計(jì)算確定[17]。本研究以Db4作為小波基，根據(jù)母小波濾波器長(zhǎng)度計(jì)算得最佳分解尺度為4層。但是該理論計(jì)算得到的最佳分解尺度只是一個(gè)概率估計(jì)，可能存在偏差。因此，綜合小波分解的特點(diǎn)及本研究數(shù)據(jù)特點(diǎn)，確定小波分解尺度為5層。

2.2 特征波段提取

以波段與重金屬含量的顯著相關(guān)程度作為判斷特征波段的依據(jù)。在軟件Matlab中編寫程序?qū)崿F(xiàn)小波分解和重構(gòu)，使用重構(gòu)后的信號(hào)與重金屬含量進(jìn)行相關(guān)性分析提取特征波段，提取的主要過(guò)程為：

1)采用Db4小波基對(duì)原始土壤光譜進(jìn)行5層離散分解，獲得5個(gè)高頻系數(shù)和1個(gè)低頻系數(shù)。

2)在Matlab中采用Wrcoef函數(shù)分別對(duì)高頻系數(shù)和低頻系數(shù)進(jìn)行單層系數(shù)重構(gòu)，獲得不同的重構(gòu)信號(hào)。

3)匯總各層重構(gòu)信號(hào)，將1～5層高頻信息分別命名為B1、B2、B3、B4、B5。由于屬于痕量級(jí)的土壤重金屬主要體現(xiàn)在高頻系數(shù)中，同時(shí)為充分利用高頻信息，研究將與重金屬含量存在顯著相關(guān)的各層高頻信號(hào)進(jìn)行相加，將得到的數(shù)據(jù)命名為Bn～Bm加和。

4)將得到的各組數(shù)據(jù)分別與重金屬鉻含量進(jìn)行皮爾森相關(guān)性分析，通過(guò)0.05雙側(cè)顯著性檢驗(yàn)的波段作為特征波段。

3 實(shí)驗(yàn)與分析

3.1 相關(guān)性分析

如前述，有效樣本總量為38個(gè)，隨機(jī)選取27個(gè)樣本作為建模集，11個(gè)樣本作為驗(yàn)證集。將B1～B5各層重構(gòu)信號(hào)、B1～B4加和、一階微分處理后的土壤光譜分別與重金屬鉻含量進(jìn)行皮爾森相關(guān)性分析，結(jié)果如圖2—圖8所示。

圖2 B1層波段與Cr相關(guān)系數(shù)

圖3 B2層波段與Cr相關(guān)系數(shù)

圖4 B3層波段與Cr相關(guān)系數(shù)

圖5 B4層波段與Cr相關(guān)系數(shù)

圖6 B5層波段與Cr相關(guān)系數(shù)

圖7 B6層波段與Cr相關(guān)系數(shù)

圖8 一階微分處理后全波段與Cr相關(guān)系數(shù)

統(tǒng)計(jì)各組數(shù)據(jù)的特征波段個(gè)數(shù)，結(jié)果如表1所示。單層高頻系數(shù)重構(gòu)信號(hào)中，B4層重構(gòu)信號(hào)中特征波段的個(gè)數(shù)最多，達(dá)到22個(gè)。其次是B3層重構(gòu)信號(hào)，特征波段的個(gè)數(shù)為13個(gè)。B5層重構(gòu)信號(hào)不存在特征波段，說(shuō)明在第5層重構(gòu)信號(hào)中沒(méi)有重金屬鉻的信息，可以視為噪聲去除?；谝浑A微分處理后的土壤光譜提取特征波段，共計(jì)11個(gè)。統(tǒng)計(jì)分析B1～B4各層重構(gòu)信號(hào)與重金屬鉻含量都存在顯著相關(guān)，因此選擇B1～B4層的重構(gòu)信號(hào)進(jìn)行加和運(yùn)算，獲得B1～B4層加和數(shù)據(jù)?；贐1～B4層加和數(shù)據(jù)提取特征波段，結(jié)果存在24個(gè)顯著相關(guān)波段。

表1 重金屬含量與各層重構(gòu)信號(hào)存在相關(guān)的波段個(gè)數(shù)

3.2 多元線性建模

選擇B1～B4層加和、B4層及一階微分處理后數(shù)據(jù)分別建立多元線性模型。建模結(jié)果如表2所示。

表2 不同的小波重構(gòu)信號(hào)建立多元線性模型的結(jié)果

基于B1～B4層加和建立的多元線性預(yù)測(cè)模型的決定系數(shù)為0.916?；贐4層重構(gòu)信號(hào)建立的多元線性模型的決定系數(shù)為0.339。基于一階微分處理光譜信息建立的多元線性模型的決定系數(shù)為0.521。結(jié)合表2發(fā)現(xiàn)，B4層重構(gòu)信號(hào)的特征波段的個(gè)數(shù)大于一階微分處理后特征波段的個(gè)數(shù)，而基于B4層重構(gòu)信號(hào)建立的多元線性模型的決定系數(shù)卻小于基于一階微分處理后的光譜信息所建立模型的決定系數(shù)。出現(xiàn)這樣的結(jié)果，可能是因?yàn)樽鳛樽宰兞康奶卣鞑ǘ沃g存在嚴(yán)重的共線性問(wèn)題。

方差膨脹因子(VIF)是衡量多元線性回歸模型中共線性嚴(yán)重程度的一種度量。一般情況下方差膨脹因子大于10，表明模型中自變量之間存在很強(qiáng)的共線性問(wèn)題[18]。分別對(duì)B1～B4層加和、B4層、一階微分處理土壤光譜提取的特征光譜進(jìn)行共線性診斷，結(jié)果如表3所示。分析表3可得，B4層、B1～B4層加和的特征波段，方差膨脹因子(VIF)遠(yuǎn)大于10，基于一階微分處理提取的特征波段存在多個(gè)波段的VIF大于10。因此，基于B4層、B1～B4層加和及一階微分處理后光譜數(shù)據(jù)提取的特征波段之間都存在共線性問(wèn)題，其中，基于B4層、B1～B4層加和提取的特征波段之間共線性尤為嚴(yán)重。

表3 方差膨脹因子匯總

3.3 偏最小二乘模型建立

偏最小二乘回歸模型是一種多因變量對(duì)多自變量的回歸建模方法，同時(shí)包含多元線性回歸分析、典型相關(guān)分析和主成分分析的特點(diǎn)，能夠較好的解決變量之間存在多重相關(guān)性或樣本點(diǎn)數(shù)過(guò)少的回歸難題的方法[19-20]。為消除特征波段之間的共線性影響，建立可靠的預(yù)測(cè)模型，選擇PLS作為建模方法?；贐1～B4層加和、B4層和一階微分處理后的三組數(shù)據(jù)分別建立PLS預(yù)測(cè)模型，建模結(jié)果如表4所示。

分析表4可知，基于一階導(dǎo)數(shù)處理后的數(shù)據(jù)建立PLS預(yù)測(cè)模型，在8個(gè)潛在因子、自變量累積解釋能力為99.4%的前提下，PLS模型的決定系數(shù)R2為0.519。基于B4層重構(gòu)信號(hào)進(jìn)行PLS建模，在潛在因子為7個(gè)、自變量累積解釋能力100%的前提下，模型的決定系數(shù)R2為0.481?；贐1～B4層加和建模結(jié)果為：在潛在因子為9個(gè)、自變量累積解釋能力99.5%的前提下，模型的決定系數(shù)R2為0.605?；贐1～B4層加和建模結(jié)果比基于B4層重構(gòu)信號(hào)建模的R2提升了0.124，比基于一階導(dǎo)數(shù)處理后的數(shù)據(jù)建模R2提升0.086。結(jié)果說(shuō)明，直接在原始光譜的基礎(chǔ)上進(jìn)行小波分解與重構(gòu)，提取特征波段，利用特征波段與重金屬含量之間的相關(guān)性建立PLS模型，可以有效預(yù)測(cè)重金屬鉻的含量。同時(shí)，與常用的一階微分變換的方法相比，基于小波分解與重構(gòu)提取的特征波段建立PLS模型的預(yù)測(cè)精度有所提升。

表4 不同小波重構(gòu)信號(hào)建立PLS模型結(jié)果

3.4 模型驗(yàn)證與分析

根據(jù)建模集所建立的PLS模型系數(shù)，得到采用一階微分光譜變換所提取的特征波段與重金屬鉻含量相關(guān)性建立的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層加和提取的特征波段與重金屬鉻含量的相關(guān)性建立的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)

基于驗(yàn)證集數(shù)據(jù)分別驗(yàn)證以上兩個(gè)模型，計(jì)算結(jié)果如表5所示。

分析表5可知，在驗(yàn)證集中，基于B1～B4層加和建立PLS模型的決定系數(shù)R2為0.511，比基于一階微分建立的PLS模型的決定系數(shù)(R2)0.351提高0.16。結(jié)果表明：利用小波變換處理原始光譜建立偏最小二乘模型具有較高的精度。

表5 預(yù)測(cè)集建模精度匯總

4 結(jié)束語(yǔ)

本研究采用基于土壤重金屬原始光譜進(jìn)行小波分解與重構(gòu)的方法，提取特征波段，建立PLS預(yù)測(cè)建模，對(duì)土壤重金屬鉻的含量進(jìn)行定量反演,并與基于常用的一階微分處理土壤光譜建立PLS模型進(jìn)行比較。結(jié)果表明：

1)利用小波分解與重構(gòu)提取的有用信息可以有效預(yù)測(cè)土壤重金屬鉻的含量；

2)對(duì)比基于常用的一階微分處理提取的特征波段和基于小波變換處理所提取的特征波段，發(fā)現(xiàn)利用小波變換處理所提取的特征波段與重金屬鉻的含量相關(guān)性建立的PLS模型具有較高精度。

雖然小波分解與重構(gòu)的方法可以有效提取光譜中有用信息，但本研究是基于室內(nèi)光譜開(kāi)展的，對(duì)于復(fù)雜的室外環(huán)境，小波分解與重構(gòu)提取特征波段的方法是否適用仍需在今后研究中進(jìn)行探討。