基于小波變換和BP神經(jīng)網(wǎng)絡(luò)的水稻冠層重金屬含量反演

李旭青 李 龍 莊連英 劉瑋琦 劉湘南 李 杰

(1.北華航天工業(yè)學院計算機與遙感信息技術(shù)學院， 廊坊 065000;2.河北省航天遙感信息處理與應(yīng)用協(xié)同創(chuàng)新中心， 廊坊 065000;3.中國地質(zhì)大學(北京)信息工程學院， 北京 100083)

0 引言

食品安全是關(guān)系國家安全穩(wěn)定與民生的重要問題。近年來，隨著城市化進程的不斷推進，工廠大量擴張、農(nóng)業(yè)用地大量占用使得農(nóng)作物被污染程度越來越嚴重，污水直接灌溉更是農(nóng)作物受重金屬污染的主要原因之一[1]。根據(jù)2018年最新研究結(jié)果，我國主要糧食產(chǎn)區(qū)耕地土壤重金屬的點位超標率為21.49%，整體以輕度污染為主，其中輕度、中度和重度污染比例分別為13.97%、2.50%和5.02%[2]。重金屬是一類不能被微生物植物等分解的物質(zhì)，但其卻可以通過生物鏈的層級作用而進入人體，并且在人體內(nèi)產(chǎn)生大量聚集，最后可能引發(fā)人體的重金屬慢性中毒[3]。因此，迫切需要一種能夠進行大面積監(jiān)測土地污染情況的方法，為國家食品安全與國計民生的穩(wěn)定發(fā)展提供依據(jù)。利用遙感手段監(jiān)測耕地土壤重金屬污染狀況是一種解決方案。

傳統(tǒng)監(jiān)測方法是實地采樣與實驗室化驗相結(jié)合，測量精度雖高，卻耗時耗力，不具備速度快、范圍廣、成本低的優(yōu)勢，且無法較好地獲取空間上的連續(xù)分布信息。相比較而言，遙感作為空間技術(shù)具有視野開闊、信息豐富、可快速獲取數(shù)據(jù)的特點，為遙感技術(shù)在土壤重金屬污染調(diào)查與監(jiān)測方面發(fā)揮更大作用提供了可能，目前已經(jīng)成為解決農(nóng)田面源重金屬污染問題的最優(yōu)方法之一[4-7]。一些學者開始對“自然農(nóng)田生態(tài)系統(tǒng)中”耕地土壤及作物重金屬污染狀況開展遙感監(jiān)測方面的研究，目前涉及到的方法有：探求敏感波段或光譜指數(shù)與作物重金屬含量的關(guān)系[8-10]，利用多維光譜指數(shù)空間[9]、小波變換、分形等方法選取受重金屬脅迫后光譜奇異性特征，從而構(gòu)建模型[11-12]。這些研究都是從不同角度探索受重金屬脅迫后作物在生理生態(tài)參數(shù)上的變化，通過構(gòu)建光譜指數(shù)提取變化信息從而構(gòu)建模型，說明基于遙感技術(shù)建立耕地土壤及作物重金屬含量反演模型在機理和技術(shù)上可行。

本文利用小波分析法(Db-5)進行小波變換，從中選取具有異常光譜特征的奇異點，利用奇異點對應(yīng)波段(716、745、766 nm)的光譜反射率構(gòu)建BP神經(jīng)網(wǎng)絡(luò)模型，對水稻冠層4種重金屬含量進行反演，并對利用模型得到的預(yù)測值與實測值進行相關(guān)性分析。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

研究區(qū)位于湘江下游地段的湖南省株洲市，地處湖南省的東部。地理坐標位于北緯26°03′05″～28°01′07″，東經(jīng)112°57′30″～114°07′15″。株洲市四季分明，夏季高溫多雨，冬季溫和潮濕，光線充足，雨量充沛，是亞熱帶季風性濕潤氣候。土壤的主要類型是一種有機質(zhì)質(zhì)量分數(shù)(2%～3%)非常大的紅壤，是適合作物生長的土壤類型。該地區(qū)年平均無霜期達286 d以上，年平均氣溫16～18℃，適合多種農(nóng)作物的種植與栽培，是湖南省的主要糧食產(chǎn)區(qū)，同時也是國家重要的商品糧食基地。主要種植的水稻類型是雜交稻博優(yōu)9083，其生長期大約為100 d，通常移栽時間為6月上旬，收獲時間為9月中旬。主要生長期為分蘗期、拔節(jié)孕穗期、抽穗揚花期和成熟期4個時期。株洲市是一個重工業(yè)基地，礦產(chǎn)開發(fā)、冶煉等造成流域水、土、作物污染，嚴重影響了農(nóng)業(yè)生產(chǎn)[13]。湘江及其支流的工業(yè)污水被用來直接灌溉農(nóng)田，造成稻田不可逆轉(zhuǎn)性的重金屬污染[14]。

本實驗在株洲市設(shè)置了3個實驗區(qū)：蘆淞區(qū)五里墩鄉(xiāng)陳家壩村、株洲縣南陽橋鄉(xiāng)橫江村、天元區(qū)馬加河鎮(zhèn)中路社區(qū)(標記地點為A、B、C)，如圖1所示。3個實驗區(qū)均為典型的水稻種植區(qū)域，實驗田面積為1 km×1 km，根據(jù)GB 15618—2018《土壤環(huán)境質(zhì)量標準》，將3塊實驗田重金屬脅迫等級分別定為：無脅迫、中度脅迫、重度脅迫3個污染水平，如表1所示。實驗區(qū)內(nèi)的溫濕度、降水、氣候等一切自然地理情況相同。同時，水稻灌溉、施肥等對水稻生長狀態(tài)等有影響的各種環(huán)境條件也基本一致。

圖1 研究區(qū)地理位置圖Fig.1 Location map of experiment sites

表1 實驗田基本信息Tab.1 Basic information of test fields

1.2 數(shù)據(jù)采集

于2014年選擇晴朗無風天氣，分別在分蘗期(7月3日)、拔節(jié)孕穗期(7月29日)、抽穗揚花期(8月29日)和成熟期(9月16日)采集水稻光譜數(shù)據(jù)。

采集時間為10:30—14:00。選用FieldSpec 4光譜儀采集光譜數(shù)據(jù)。每次測量前，先打開光譜儀進行預(yù)熱，預(yù)熱時間為10～30 min，然后用標準白板進行校正，其中標準白板的反射率為1，視場角為10°，光譜采集過程中，傳感器探頭始終保持垂直向下，距水稻冠層約1 m處，光譜儀的波段范圍為350～2 500 nm，分辨率為10 nm。350～1 000 nm波段范圍內(nèi)采樣間隔為1.4 nm，1 000～2 500 nm波段范圍內(nèi)采樣間隔為2 nm。采樣點均勻分布于實驗田。每塊樣地選擇30個采樣點，每次測量10條光譜數(shù)據(jù)，并對10條光譜數(shù)據(jù)進行取平均處理，平均值作為該樣本的結(jié)果數(shù)據(jù)，因為水分吸收帶噪聲嚴重，實驗中選取的波段范圍是350～1 000 nm。

在采集高光譜數(shù)據(jù)的同時，采集對應(yīng)水稻植株樣本，本實驗中所有重金屬含量的測定標準均參考國家相關(guān)測定標準，分別為：GB/T 5009.15—2003(鎘)、GB/T 5009.12—2003(鉛)、GB/T 5009.17—2003(汞)、GB/T 5009.11—2003(砷)。

1.3 實驗方法

國內(nèi)外基于小波變換的遙感圖像去噪方法眾多，小波變換在遙感圖像去噪方面已經(jīng)成為學者們研究的熱點，不同學者對小波閾值、小波系數(shù)等進行了改進優(yōu)化[15-20]。小波變換的主要特點是可以表征信號的局部特征，而信號的突變點和奇異點等不規(guī)則部分通常包含重要信息。MAIREG等[21]將小波分析應(yīng)用于反射光譜并利用逐步多元回歸方法反演葉片葉綠素含量，取得了較好的反演精度。小波變換可通過Matlab快速實現(xiàn)。已有研究表明，對實驗區(qū)光譜數(shù)據(jù)在波長350～1 000 nm范圍內(nèi)做Db-5小波5層分解，可使光譜信號的噪聲減弱，奇異信息增強。

自MORLET和ARENS初步提出小波概念后，小波分析被廣泛應(yīng)用于圖像處理、信號處理等眾多領(lǐng)域。小波可以把數(shù)據(jù)分解成頻率不同的組成成分，通過對數(shù)據(jù)時域和頻域的分解來表示數(shù)據(jù)局部特征，本文對光譜數(shù)據(jù)分析是指對高光譜數(shù)據(jù)波長進行小波分析，對母小波ψ(t)進行伸縮和平移，產(chǎn)生一組小波基函數(shù)。

(1)

式中a——小波寬度b——小波位置

ψab——小波基函數(shù)

對任意函數(shù)的小波變換表示為

(2)

式中wf——小波變換函數(shù)

f——反射光譜R——窗口尺寸

BP神經(jīng)網(wǎng)絡(luò)的基本思想是梯度下降法，利用梯度搜索技術(shù)，使網(wǎng)絡(luò)的實際輸出值和期望輸出值的誤差均方差為最小。基本BP算法包括信號的前向傳播和誤差的反向傳播兩個過程。即計算誤差輸出時按從輸入到輸出的方向進行，而調(diào)整權(quán)值和閾值則從輸出到輸入的方向進行。基本參數(shù)包括網(wǎng)絡(luò)層數(shù)、網(wǎng)絡(luò)初始權(quán)值和閾值、隱含層層數(shù)、最大訓練步數(shù)、網(wǎng)絡(luò)學習步長和速率等。本研究中，選擇輸入層、隱含層和輸出層建立反向傳播算法的BP 神經(jīng)網(wǎng)絡(luò)，需應(yīng)用newff函數(shù)，在確保模型精度的情況下，神經(jīng)元個數(shù)越少, 越有助于避免過擬合現(xiàn)象, 提高模型的泛化能力。人工神經(jīng)網(wǎng)絡(luò)具有極強的非線性映射能力，非常適合于函數(shù)逼近，即找出2組數(shù)據(jù)之間的關(guān)系[22-26]。人工神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)由輸入層、輸出層、1或2個隱含層構(gòu)成，網(wǎng)絡(luò)的好壞取決于網(wǎng)絡(luò)結(jié)構(gòu)、隱含層個數(shù)、神經(jīng)元個數(shù)、傳遞函數(shù)、訓練函數(shù)等[26]，在實驗中需要通過不斷地嘗試才能進行確定。水稻受到重金屬污染脅迫會有光譜響應(yīng)，水稻重金屬含量和光譜響應(yīng)之間呈復雜的非線性關(guān)系，其復雜的非線性關(guān)系用統(tǒng)計回歸模型難以模擬[23]。

由于在“自然農(nóng)田生態(tài)系統(tǒng)”中，重金屬對作物的影響在光譜上的反映不便于捕捉，而小波分析又具有放大微小信號變化的特性。已有研究表明，Daubechies小波系中的Db-5小波函數(shù)能夠有效檢測重金屬污染脅迫下的水稻光譜異常信號[27]，而且小波變換能夠把其他對水稻原始光譜信號有影響的因素進行有效的排除，如背景影響、設(shè)備工作時產(chǎn)生的噪聲以及大氣的散射、反射和吸收。所以本實驗采用小波變換方法，利用Db-5小波基對水稻原始光譜進行小波變換處理，將敏感波段波譜反射率作為BP神經(jīng)網(wǎng)絡(luò)模型輸入層，建立水稻冠層重金屬含量反演模型。

2 水稻重金屬污染脅迫水平評估模型

2.1 光譜反射率奇異點分析

首先分別對4次測量的若干條光譜數(shù)據(jù)進行均值處理，求得實驗田平均水稻光譜反射率曲線。對原始光譜反射率信號進行Db-5小波基的小波分解變換，分解層級設(shè)置為5，把原始的光譜反射率信號分解為兩部分，即高頻信號與低頻信號。經(jīng)過多層分解之后，即可把原始光譜反射率信號(圖2)轉(zhuǎn)換為多個子信號(圖3)。根據(jù)前人研究以及前文分析可知，原始水稻光譜反射率曲線的異常信息突變點和小波變換后的奇異點相互對應(yīng)。一些作物受到重金屬污染的異常光譜反射率信息可以通過這些奇異點被放大表現(xiàn)出來。對受3種不同重金屬脅迫污染的地區(qū)，進行數(shù)據(jù)奇異點分析，結(jié)果如表2所示。

圖2 A、B、C各時期原始光譜反射率曲線Fig.2 Original spectral curves of A, B and C in different periods

圖3 A、B、C各時期小波變換結(jié)果Fig.3 Wavelet transform result of A，B，C in each period

據(jù)表1可知，實驗區(qū)只有土壤重金屬含量超標而無其他環(huán)境因素異常，在350～1 000 nm范圍內(nèi)光譜奇異性是由重金屬污染脅迫所致。在株洲市的3個實驗田中，在不同的重金屬脅迫水平下生長的水稻，其高光譜數(shù)據(jù)的奇異點所對應(yīng)的波段會表現(xiàn)出一些不同之處，但依舊會在某些特定的波段出現(xiàn)相同的奇異點。通過對4次測量數(shù)據(jù)的對比研究，得出3個實驗田里水稻光譜反射率相同的極大值波段有492、555、605、694、745、798、877、937 nm共8個奇異點對應(yīng)波段，3個實驗田共同的極小值波段有516、578、624、716、766、817、907、963 nm共8個奇異點對應(yīng)波段。處于相同的重金屬脅迫水平下生長的水稻，其對應(yīng)的奇異點的分布情況也與之對應(yīng)相同。反之，處于不同的重金屬脅迫水平下生長的水稻，其對應(yīng)的奇異點的分布情況有一定的差異。同理，不同生長環(huán)境、生長區(qū)域下生長的水稻，對應(yīng)的奇異點的分布情況也不盡相同[28]。所以通過對奇異點進行分析研究，可以快速找到其對應(yīng)波段。

表2 實驗區(qū)水稻光譜反射率奇異點對應(yīng)波段Tab.2 Band corresponding to singular point of rice spectral reflectance in test fields

通過對株洲市3塊實驗田水稻光譜反射率奇異點對應(yīng)波段(表2)進行統(tǒng)計分析，選出在不同重金屬污染脅迫條件下，光譜反射率奇異點相同的幾個波段作為構(gòu)建水稻冠層重金屬含量反演模型的自變量。這些相同的波段分別是555、578、716、745、766 nm 5個波段。相關(guān)研究表明，農(nóng)作物在重金屬污染脅迫條件下，光譜反射率在750 nm波段附近表現(xiàn)最為明顯[29-31]。因此，選擇716、745、766 nm 3個波段來構(gòu)建水稻冠層重金屬含量反演模型。

2.2 水稻重金屬污染脅迫水平評估模型構(gòu)建

輸入層有3個輸入量，分別為X1、X2、X3；隱含層1有3個神經(jīng)元節(jié)點，分別為L11、L12、L13，隱含層2有兩個神經(jīng)元節(jié)點，分別為L21、L22，隱含層3有1個神經(jīng)元節(jié)點L31；最后為輸出層Y，如圖4所示。采用實驗區(qū)地面實測光譜數(shù)據(jù)作為BP神經(jīng)網(wǎng)絡(luò)模型的訓練集建立水稻冠層重金屬含量反演模型，A、B、C實驗區(qū)分別有16、29、13個樣本。經(jīng)反復實驗，結(jié)果表明，學習函數(shù)為learngdm(梯度下降動量函數(shù))、訓練函數(shù)為traingdm(梯度下降動量BP算法函數(shù))、學習速率為0.2、動量系數(shù)為0.2、隱含層神經(jīng)元傳遞函數(shù)均為tansig、輸出層神經(jīng)元傳遞函數(shù)為logsig(對數(shù)傳遞函數(shù))時，模型精度較高。

選取約2/3的樣本(40個)作為訓練集進行建模，約1/3的樣本(18個)作為驗證集對模型結(jié)果進行檢驗。以決定系數(shù)R2和均方根誤差(RMSE )作為評價指標，對模型結(jié)果的優(yōu)越性進行評價(表3)。將基于BP神經(jīng)網(wǎng)絡(luò)構(gòu)建的水稻冠層重金屬含量反演模型預(yù)測結(jié)果與實際污染水平進行比較，效果較好(表4)。結(jié)果顯示，A區(qū)無重金屬污染脅迫，B區(qū)為中度脅迫，C區(qū)為重度脅迫，這一結(jié)果與研究區(qū)實際情況相符。最終表明，水稻冠層重金屬含量與水稻光譜之間并非簡單的線性關(guān)系，利用BP神經(jīng)網(wǎng)絡(luò)方法評估污染脅迫有較好效果。

表3 基于BP神經(jīng)網(wǎng)絡(luò)的污染脅迫水平估測模型Tab.3 Pollution stress level estimation models based on BP neural network model

表4 BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)模型污染脅迫水平預(yù)測結(jié)果Tab.4 Prediction results of pollution stress level using BP neural network model

3 結(jié)論

(1)Db-5小波變換能夠有效增強作物受到重金屬脅迫污染后的弱光譜信息，選取716、745、766 nm 3個波段構(gòu)建水稻冠層重金屬含量反演模型。在反演模型中，通過光譜反射率對水稻冠層重金屬含量進行預(yù)測，精度較好，相關(guān)性較高，模型具有一定的實用性。

(2)在BP神經(jīng)網(wǎng)絡(luò)模型中，當學習函數(shù)為learngdm(梯度下降動量函數(shù))、訓練函數(shù)為traingdm(梯度下降動量BP算法函數(shù))、學習速率為0.2、動量系數(shù)為0.2、隱含層神經(jīng)元傳遞函數(shù)均為tansig、輸出層神經(jīng)元傳遞函數(shù)為logsig(對數(shù)傳遞函數(shù))時，模型精度最高。