基于連續(xù)小波變換下的土壤有害元素砷含量估測

王雪梅，玉米提·買明，黃曉宇，李 銳，劉 東

1.新疆師范大學地理科學與旅游學院, 新疆 烏魯木齊 830054 2.新疆維吾爾自治區(qū)重點實驗室“新疆干旱區(qū)湖泊環(huán)境與資源實驗室”, 新疆 烏魯木齊 830054

引 言

隨著城市化進程的加劇，工業(yè)和生活污水的排放以及農藥等大量使用導致土壤重金屬污染日益嚴重[1]。近幾十年來，新疆石油化工產業(yè)和綠洲農業(yè)的快速發(fā)展，使得土壤中的鉻、汞、鉛、銅、鎘和砷等重金屬元素在土壤環(huán)境中普遍存在。渭干河—庫車河三角洲綠洲是新疆南疆重要的農業(yè)灌溉區(qū)和石油化工園區(qū)，土壤受到重金屬污染的風險較大。在自然界中，有害元素砷(As)常以離子形態(tài)存在于土壤溶液中，對動植物和人體健康危害極大，而以膠狀物質存在于土壤中的砷容易被植物吸收并經食物鏈最終到達人體，危害生命健康[2]。開展土壤重金屬As含量的檢測對防治土壤重金屬污染、保障人體健康，以及構建健康生態(tài)環(huán)境具有十分重要的現實意義[3]。

檢測土壤有害元素As的傳統(tǒng)方法主要采用分光光度法、化學分析法以及原子熒光光譜法等實驗手段進行，雖檢測精度高，但儀器設備和人力耗費大且不易推廣，隨著高光譜遙感技術的快速發(fā)展，為土壤有害元素砷的遙感估測提供了有力的技術手段[4]。土壤樣品的高光譜采集受眾多環(huán)境因素的影響，光譜數據中存在許多細節(jié)噪聲。研究最初嘗試對土壤原始光譜數據進行倒數、倒數對數以及一階微分和二階微分等多種數學變換處理，從而消除土壤樣品中的環(huán)境噪聲，增強土壤光譜數據中的有效信號[5-6]。隨著研究的深入，發(fā)現單純采用傳統(tǒng)的數學變換方法對原始光譜數據進行處理，雖在一定程度上提高了模型的估測精度，卻難以有效去除白噪聲[7]。相關研究顯示，經過連續(xù)小波變換(continuous wavelet transformation，CWT)后的光譜反射率與土壤屬性具有更強的關聯性，可實現光譜信號的近似特征和細節(jié)特征的有效分離，構建的模型估測精度更高[8-9]。Zhang[10]、肖艷[11]等研究認為光譜數據的連續(xù)小波變換對有效波段的挖掘和反演模型的估測都明顯優(yōu)于傳統(tǒng)的數學變換。為了探討在連續(xù)小波變換下進行高光譜遙感估測干旱區(qū)綠洲土壤As含量的最佳小波基函數及最優(yōu)分解尺度，本研究基于野外調查的98個土壤樣品的光譜數據以及有害元素砷含量，通過4種小波基函數對原始光譜采用10個尺度進行連續(xù)小波變換處理，篩選出與土壤As含量密切相關的敏感小波系數。以敏感小波系數為自變量，土壤As為因變量，采用偏最小二乘回歸(partial least squares regression，PLSR)、支持向量機回歸(super vector machine regression，SVMR)、BP神經網絡(back propagation neural network，BPNN)和隨機森林回歸(random forest regression，RFR)方法構建干旱區(qū)綠洲土壤有害元素As含量的反演模型，為有害元素As的準確估算提供依據。

1 實驗部分

1.1 土樣采集與制備

渭干河—庫車河三角洲綠洲位于新疆塔里木盆地北緣，是一個典型的干旱區(qū)扇形平原綠洲，轄區(qū)包括新疆維吾爾自治區(qū)阿克蘇地區(qū)的庫車市、沙雅縣和新和縣，地理位置為北緯39°30′—42°40′，東經81°27′—84°07′(圖1)。該區(qū)多年平均降水量為51.6 mm，多年平均蒸發(fā)量可達2 123 mm，屬干旱與極干旱地區(qū)。綠洲植被主要以小麥(TriticumaestivumL)、棉花(Gossypiumspp)和玉米(ZeamaysLinn)等農作物，以及檉柳(Tamarixramosissima)、鹽爪爪(Kalidiumfoliatum)、鹽穗木(Halostachyscaspica)、駱駝刺(Alhagisparsifolia)、蘆葦(PhragmitesaustraliasTrin)和胡楊(Populuseuphratica)等荒漠植被為主。土壤類型包括棕漠土、灌淤土、沼澤土和鹽土等多種類型。以遙感影像為參考，在室內提前合理布點，2019年7月中到下旬進行野外調查和土壤樣品采集。土壤剖面深度為0～20 cm，共采集98個土壤樣品，并分別采用GPS系統(tǒng)對樣點進行準確定位。采集的樣品裝入做好標記的土樣袋中。土壤樣品在實驗室經風干處理，挑出雜物，研磨過1 mm篩，送至新疆維吾爾自治區(qū)分析測試研究院，由工作人員依據GB/T 22105.2—2008第2部分測定標準，采用原子熒光光譜法進行土壤總砷的測定。

圖1 采樣點分布圖

1.2 光譜數據的測定與預處理

采用ASDFieldSpec3光譜儀進行研究區(qū)土壤樣品的高光譜數據采集。在野外測定光譜時，為避免溫差干擾數據的采集，首先要對光譜儀進行預熱和白板校正。由于環(huán)境因素會使光譜數據在采集時產生干擾噪聲，故要遠離高大物體。按照標簽順序將土壤樣品均勻攤開并完全覆蓋在50 cm×50 cm的牛皮紙上，探頭與采樣點的垂直上方距離為15 cm。為了使采集的光譜數據更為準確，每個樣品需要重復觀測10次，取平均值作為該樣品的光譜曲線值。為了消除儀器兩端產生的較大噪聲，剔除350～399和2 451～2 500 nm的光譜數據，采集的有效光譜范圍為400～2 450 nm。進一步對98個樣品的光譜反射率進行9點Savitzky-Golay濾波平滑處理，通過觀察98條光譜曲線，將受水分吸收和大氣影響的1 341～1 400和1 811～1 950 nm波段的異常數值進一步剔除。

1.3 基于光譜數據的連續(xù)小波變換

連續(xù)小波變換(continuous wavelet transformation，CWT)也稱為積分小波變換(integral wavelet transform，IWT)，常應用于各類信號的處理，具有較高的分辨性和適應性，其變換公式見式(1)和式(2)

(1)

(2)

式中：a為尺度因子；b為平移因子；Wf(a，b)為小波系數，采用二維矩陣表示；λ為波長；f(λ)為土壤光譜反射率；Ψa，b為小波基函數。CWT作為一種信號處理方法，可有效去除白噪音，然而在光譜分解過程中，小波基函數的選擇一直沒有統(tǒng)一的說法。相關研究發(fā)現，采用不同的小波基函數對光譜數據進行各種尺度的分解，會直接導致實驗結果不一致，目前多以經驗和反復實驗比較進行選擇[12]。研究中使用4種小波基函數bior1.3，db4，gaus4與mexh對原始光譜反射率R進行10個尺度的光譜分解，選取敏感小波系數作為自變量構建干旱區(qū)土壤有害物質砷含量的高光譜估測模型。

1.4 模型的建立與檢驗

采用偏最小二乘回歸(PLSR)、支持向量機回歸(SVMR)、BP神經網絡(BPNN)和隨機森林回歸(RFR)4種建模方法對研究區(qū)土壤有害物質砷含量進行高光譜反演，其中PLSR模型的最佳潛在變量(Latent variables)按照小波系數的數量而定；SVMR模型以徑向基函數radial作為核函數，參數Cost、Epsilon和Gamma可通過反復訓練和調試最終確定最佳值；BP神經網絡模型的訓練函數選擇newff，訓練過程中最大迭代次數為500，訓練目標最小誤差為0.001，學習率為0.01，通過不斷改變隱含層的神經元個數訓練該模型以獲得最佳效果；隨機森林回歸(RFR)采用randomForest函數，通過指定決策樹的數目ntree來反復多次訓練尋找最優(yōu)的模型參數mtry。模型的精度和預測能力通過決定系數(R2)、均方根誤差(RMSE)和相對分析誤差(RPD)進行評價。R2和RPD值越大，RMSE值越小，則模型的精度越高，估測能力越強。

2 結果與討論

2.1 有害物質砷含量與光譜特征分析

對研究區(qū)98個土壤樣本的有害物質砷含量進行基本統(tǒng)計(表1)，有害元素砷含量的最大值為22.8 mg·kg-1，最小值為4.35 mg·kg-1，平均水平為11.67 mg·kg-1，屬于中等空間變異。由偏度系數和峰度系數發(fā)現，砷含量總體偏高，且呈高聚集狀態(tài)。根據建模需要將總樣本隨機分為68個訓練樣本和30個驗證樣本。與新疆土壤砷的背景值和國家標準比較后發(fā)現，研究區(qū)As含量的平均水平已超過新疆土壤背景值(11.2 mg·kg-1)，但還未超出國家標準(40 mg·kg-1)，說明As已過量沉積形成輕度污染。分析認為隨著研究區(qū)工業(yè)化進程的加劇以及農藥的過量使用，使耕層土壤中的有害物質As含量在土壤中不斷蓄積，成為影響農業(yè)可持續(xù)發(fā)展和危及人體健康的重要因素。

表1 土壤樣品的基本統(tǒng)計特征

為了分析砷含量與土壤光譜反射率之間的關聯性，將98個樣本按照砷含量由低到高劃分為5個等級，Ⅰ級：As<8 mg·kg-1，Ⅱ級：8 mg·kg-1≤As<10 mg·kg-1，Ⅲ級：10 mg·kg-1≤As<12 mg·kg-1，Ⅳ級：12 mg·kg-1≤As≤14 mg·kg-1，Ⅴ級：As≥14 mg·kg-1。通過分析不同砷含量等級的光譜反射率曲線(圖3)，發(fā)現隨著砷含量的增加，光譜曲線呈現出下降的趨勢，說明土壤中砷元素的含量直接影響到光譜反射率，采用高光譜數據進行土壤砷含量估測具有一定的依據。

圖2 不同砷含量等級的光譜反射率曲線

2.2 CWT下的光譜反射率與有害物質砷含量的相關分析

通過對原始光譜數據分別采用4種小波基函數進行10個尺度的連續(xù)小波變換，變換后的小波系數與As含量進行相關分析并繪制相關譜圖。圖3(a—d)是土壤As含量與不同小波基函數變換下10個尺度的相關系數二維圖，橫坐標為波長(400～2 450 nm)，縱坐標為不同分解尺度(2n，n=1, 2, 3, …, 10)，顏色表示相關系數的大小。從圖3可以看出，有害物質As含量與CWT處理后的光譜反射率存在一定程度的相關性。其中，在可見光400～700nm以及近紅外的1 100～1 700和2 200～2 400 nm附近具有較強的相關關系，且表現出3～8尺度可極大程度提高光譜與土壤屬性之間的相關性，其中4～6尺度具有最佳的光譜分解能力。

進一步分析圖3發(fā)現采用bior1.3小波基函數進行第5尺度分解的小波系數與As含量的相關性最強，能夠通過p<0.01顯著性檢驗的光譜波段有507個，且主要集中分布在可見光400～600 nm與分散分布在近紅外1 100～2 400 nm。為了得到建模所需的敏感小波系數，進一步對第5尺度分解結果通過p<0.001相關顯著性檢驗，篩選出相關系數絕對值的敏感波段共12個，波長范圍主要集中在2 343～2 354 nm，相關系數r最高可達0.687。通過對db4小波基函數的各尺度分解結果分析，發(fā)現第6尺度與As的相關性最強，通過p<0.01顯著性檢驗的光譜波段有203個，主要集中分布在可見光的500～700 nm和近紅外的2 300～2 400 nm附近；通過進一步篩選的敏感波段共有28個(p<0.001)，相關系數r的絕對值最高可達0.686。采用gaus4小波基函數進行10個尺度的光譜分解結果顯示，通過p<0.01顯著性檢驗的光譜波段在第5尺度最多，達到280個，且主要集中分布在400～600 nm以及2 300～2 400 nm；與As具有極顯著相關性且的敏感波段有16個(p<0.001)，相關系數絕對值最高可達0.669。mexh小波基函數的第4尺度的光譜分解能力最強，通過p<0.01顯著性檢驗的光譜波段達到216個，主要集中在400～600，1 100～1 300以及2 300～2 400 nm，通過p<0.001相關顯著性檢驗且的敏感波段有24個，相關系數絕對值最高可達0.689。

圖3 砷含量與小波系數間的相關性

通過上述分析可以看出，采用4種小波基函數對光譜數據進行連續(xù)小波變換后，小波系數與As含量具有較顯著的相關性，表明基于不同小波基函數的CWT能較好地放大各個波段范圍內微弱的土壤特征信息并更有效地挖掘信息。對比4種小波基函數的分解結果發(fā)現，經bior1.3小波處理后與As相關性顯著的小波系數數量遠高于其他3種小波函數(p<0.01)，說明bior1.3小波變換可更大程度上實現光譜信號的近似特征和細節(jié)特征的有效分離，提高光譜與As之間的相關性。進一步進行p<0.001相關極顯著性檢驗發(fā)現，小波基函數db4和mexh的分解效果要優(yōu)于gaus4。通過相關系數絕對值|r|>0.6可篩選出4種小波基函數下的敏感小波系數作為As含量估算的建模變量(見表2)。

表2 篩選的敏感小波系數

2.3 反演模型的建立

為進一步確定敏感小波系數與土壤As含量之間的定量關系，以各小波基函數CWT處理后篩選出的敏感小波系數為自變量，采用PLSR，SVMR，BPNN和RFR方法構建土壤As含量的高光譜估測模型，并通過決定系數(R2)、均方根誤差(RMSE)和相對分析誤差(RPD)對模型進行精度評價和可靠性檢驗。分析表3反演結果認為，4種模型整體估測效果良好，與其他3種機器學習算法相比，PLSR模型的估測精度相對較低；對比SVMR，BPNN與RFR模型發(fā)現，SVMR和BPNN模型訓練集的R2和RPD值較為接近，驗證集則表現為BPNN模型的估測能力要優(yōu)于SVMR模型，而RFR模型的訓練集和驗證集估測精度均較高，模型更穩(wěn)定。

表3 不同反演模型的估算精度

對比4種小波基函數的估測效果發(fā)現，PLSR模型中mexh函數的估測精度要高于其他3種小波基函數；在SVMR模型中，4種小波基函數的光譜分解能力一般，模型估測效果依次為mexh-24>db4-26>bior1.3-25>gaus4-25；而在BPNN和RFR模型中，以bior1.3小波基函數構建的模型具有更強的估測能力，其建模集和驗證集的R2均在0.6以上，RMSE在1.9 mg·kg-1以下，RPD值均大于1.6，說明bior1.3小波基函數對光譜數據的分解效果較佳，構建的模型具有較好的估測能力。通過對比4種小波基函數對原始光譜數據中有效信號的辨識能力，可看出bior1.3和mexh的變換效果明顯優(yōu)于db4和gaus4函數。通過綜合比較分析認為，bior1.3-25-RFR模型對As的估測精度最高，具有最強的穩(wěn)定性。該模型參數ntree和mtry的最佳值分別為500和8，其建模集和驗證集的R2分別為0.893和0.639，RMSE為1.075和1.651 mg·kg-1，RPD值分別為2.89和1.64，均大于1.6，模型具有較好的估測能力。

分別對4種建模方法下有害物質As的最佳反演結果繪制散點圖[圖4(a—d)]，發(fā)現各模型的估測值與實測值構成的散點絕大多數都分布在1∶1線附近，相對較為集中，說明這4種建模方法對有害物質As均具有良好的解釋能力。對比分析認為bior1.3-25-RFR模型的散點在1∶1線附近分布的更為集中，總體估測誤差最小，說明該模型能更好地反演土壤有害物質As的含量。

圖4 實測值與估測值的散點圖

3 結 論

適當的光譜變換方法有利于土壤信息的提取和光譜數據中有效信息的識別。作為一種光譜數據處理方法，CWT可以有效去除土壤高光譜數據中的噪音，檢測有效信號的效果明顯，可更好地挖掘出光譜數據中的微弱卻有效的隱藏信息，對土壤高光譜反演起著重要的作用[13]。Li等[14]的研究表明，土壤原始光譜反射率經CWT處理在一定程度上提高了與土壤有機質含量間的相關性，說明光譜的CWT變換可消除不確定因素對光譜信息的干擾，從土壤光譜數據中分解出穩(wěn)定的高頻信息。Zhang等[10]對4種元素(Cr，As，Ni和Cd)進行高光譜反演結果顯示，與傳統(tǒng)數學變換方法相比，CWT變換能更好地提取光譜中的有效信息，具有良好的光譜分解能力。Hong等[15]研究顯示，采用CWT變換的隨機森林模型對土壤有機碳具有更好的估算能力。由于土壤性質的復雜性以及土壤類型的多樣性，不同研究區(qū)域的土壤屬性存在較大空間差異，如何采用CWT方法選擇合適的小波基函數和最優(yōu)的分解尺度，以及小波基函數是否穩(wěn)定和具有普遍適用性仍需做進一步探討。

研究基于bior1.3，db4，gaus4和mexh 4種小波基函數進行連續(xù)小波變換處理，結合土壤As含量進行相關分析，通過篩選的敏感小波系數采用PLSR，SVMR，BPNN和RFR模型對干旱區(qū)土壤As含量進行估算，結果表明：(1)在CWT過程中，4種小波基函數在4～6尺度上的光譜分解效果明顯優(yōu)于其他尺度，與土壤As具有顯著相關的小波系數數量更多(p<0.01)；(2)對比4種小波基函數的光譜變換效果，bior1.3和mexh對光譜數據中有效信息的提取能力要強于db4和gaus4，說明這兩種小波基函數更利于挖掘光譜數據中的隱藏信息，可在一定程度上放大土壤中的有效信息；(3)與PLSR，SVMR和BPNN模型相比，RFR模型具有更高的估測精度和較好的穩(wěn)定性；土壤中As的最佳估算模型為bior1.3-25-RFR模型，其訓練集和驗證集的R2均大于0.6，均方根誤差RMSE小于1.7 mg·kg-1，RPD值大于1.6，估測能力較好。