重慶市南川區(qū)土壤錳元素遙感反演

徐 天,李 敬,劉振華*

1. 華南農(nóng)業(yè)大學(xué)資源環(huán)境學(xué)院,廣東 廣州 510000 2. 廣東省土地信息工程技術(shù)研究中心,廣東 廣州 510000

引 言

錳是植物、 動物和人體不可缺少的微量營養(yǎng)元素,在生物化學(xué)過程中起著十分重要的作用。 但錳的濃度過高將導(dǎo)致環(huán)境污染,影響植物的生長發(fā)育,使得作物產(chǎn)量減少,農(nóng)產(chǎn)品質(zhì)量下降,威脅生物多樣性,人體長期接觸高濃度錳則會出現(xiàn)組織損傷和錳中毒癥狀[1]。 傳統(tǒng)的土壤錳含量監(jiān)測方法是田間調(diào)查、 定點監(jiān)測、 室內(nèi)樣品實驗分析。 盡管該方法在監(jiān)測點上可獲得精確的土壤錳含量信息,但需要大量的人力、 財力和物力,效率低。 且就區(qū)域尺度而言,因受監(jiān)測點網(wǎng)格和插值方法影響,非監(jiān)測點的精度不高。 隨著遙感技術(shù)的發(fā)展,已有學(xué)者利用遙感技術(shù)監(jiān)測土壤錳含量。 目前基于遙感技術(shù)的土壤錳含量監(jiān)測方法主要集中于利用土壤光譜指標(biāo)反演土壤錳含量,可分為兩方面: 一是利用高光譜技術(shù)反演土壤錳含量,二是利用多光譜遙感技術(shù)反演土壤錳含量。 用高光譜技術(shù)反演土壤錳含量是基于土壤高光譜數(shù)據(jù),獲取土壤錳的光譜響應(yīng)指標(biāo),構(gòu)建光譜響應(yīng)指標(biāo)與土壤錳之間的關(guān)聯(lián)模型,反演土壤錳含量。 如: 潘勇等利用地物光譜儀采集礦區(qū)土壤光譜譜數(shù)據(jù),并經(jīng)主成分分析法降維后結(jié)合不同的建模方法對土壤錳元素進(jìn)行分級與評價,結(jié)果表明模型預(yù)測效果良好[2]。 Pandit等通過探索土壤錳與反射光譜之間的響應(yīng)關(guān)系,建立了基于偏最小二乘法的土壤錳反演模型并得到了較高的預(yù)測精度[3]。 田燁等基于大量野外實測土壤光譜數(shù)據(jù),探討了不同光譜變換方法和預(yù)測模型對土壤錳含量反演的影響,并進(jìn)一步探索了利用模擬衛(wèi)星傳感器反射率反演土壤錳含量的可行性,取得了較為理想的效果。

利用多光譜遙感技術(shù)反演土壤錳含量是基于多光譜影像,獲取土壤錳的光譜響應(yīng)指標(biāo),構(gòu)建光譜響應(yīng)指標(biāo)與土壤錳含量之間的關(guān)聯(lián)模型,反演土壤錳含量。 如: 梁遠(yuǎn)玲等利用實測高光譜模擬Landsat8影像進(jìn)行反演土壤錳含量的反演,應(yīng)用最優(yōu)模型得到的反演結(jié)果較實測插值結(jié)果的空間分布更為細(xì)致,但是反演結(jié)果只能反映基本的區(qū)域土壤錳含量的分布情況。 袁濤等利用World View-3多光譜影像數(shù)據(jù)中提取的單波段反射率與裸地的土壤樣點實測重金屬含量建立回歸模型,在裸地的建模結(jié)果中模型精度均較差,模型誤差指數(shù)較高,無法進(jìn)行高質(zhì)量的區(qū)域土壤錳含量反演。

盡管基于遙感技術(shù)的土壤錳含量的監(jiān)測研究取得了很大進(jìn)展,但現(xiàn)有研究僅利用土壤光譜反演求取土壤錳元素含量,然而由于土壤組分的復(fù)雜性,信號干擾大,造成土壤錳的光譜響應(yīng)指標(biāo)獲取困難,難以獲取精準(zhǔn)反演土壤錳的模型,鑒于土壤中微量元素直接影響植被的長勢,導(dǎo)致植被光譜的變化,因此從植被光譜角度出發(fā),探討植被光譜與土壤錳之間的關(guān)聯(lián)性,構(gòu)建土壤錳含量反演模型,實現(xiàn)區(qū)域尺度的土壤錳含量反演。

1 實驗部分

1.1 研究區(qū)概況

研究區(qū)位于重慶市南川區(qū)(東經(jīng)106°54′—107°27′,北緯28°46′—29°30′),幅員面積2 602 km2,境內(nèi)屬亞熱帶濕潤季風(fēng)氣候,氣候溫和,雨量充沛,既無嚴(yán)寒,又無酷暑,四季分明,霜雪稀少,無霜期,熱量豐富,年均氣溫16.6 ℃。 南川區(qū)境內(nèi)多山,地形以山地為主,地勢呈東南向西北傾斜,最高點金佛山風(fēng)吹嶺海拔2 251 m,最低點騎龍魚跳巖海拔340 m。 重慶市南川區(qū)主要土壤類型為紫色土、 黃壤和石灰(巖)土,有機質(zhì)含量低,磷、 鉀豐富,石灰(巖)土質(zhì)地粘重,偏堿性,土壤陽離子交換量和鹽基飽和度均高。

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

基于重慶地理數(shù)據(jù)庫隨機選取694個樣點,用梅花五點采樣法采集土樣,均勻混合后作為該點的土壤樣品,提取樣點的土壤錳含量數(shù)據(jù),土壤錳含量區(qū)間為20～4 462 mg·kg-1。 此外,在幾何校正后的Landsat影像上,提取錳含量有關(guān)的光譜數(shù)據(jù),獲取樣點的植被指數(shù)。 選取其中的555個樣點數(shù)據(jù)進(jìn)行建模與測試,用隨機抽樣法按3: 1的比例對樣點劃分了建模集[圖1(a)綠點]和測試集[圖1(b)黃點]。 為了驗證基于植被光譜響應(yīng)指標(biāo)的土壤錳含量制圖精度,將剩余的139個樣點作為制圖精度驗證點[圖1(c)紅點]。

圖1 研究區(qū)及樣本點分布Fig.1 Study area and sampling distribution

影像數(shù)據(jù)來源于地理空間數(shù)據(jù)云(https: //www.gscloud.cn)的Landsat8 OLI影像。 Landsat-8衛(wèi)星上攜帶兩個傳感器,分別是OLI陸地成像儀(Operational Land Imager)和TIRS熱紅外傳感器(Thermal Infrared Sensor),衛(wèi)星一共有11個波段,波段1～7,9～11的空間分辨率為30 m,波段8為15 m分辨率的全色波段,成像寬幅為185 km×185 km。 本研究獲取的數(shù)據(jù)成像時間為2019年8月13日,云量為0.09%,影像質(zhì)量較高,符合使用要求,為消除傳感器和大氣影響,利用ENVI 5.3軟件對Landsat8影像進(jìn)行輻射定標(biāo)和大氣校正等預(yù)處理。 此外,土地覆蓋數(shù)據(jù)來源于GlobeLand30(http://www.globallandcover.com/)。

利用Origin軟件對694個土壤樣點中土壤錳含量和植被樣點進(jìn)行描述性分析,包括樣點數(shù)據(jù)的取值范圍、 平均值、 標(biāo)準(zhǔn)差和平均絕對偏差等能夠反映數(shù)據(jù)基本情況的各項指標(biāo),結(jié)果如表1。

表1 土壤錳含量及植被樣點描述性統(tǒng)計Table 1 Descriptive statistics of soil Mn content and vegetation sampling sites

1.3 方法

1.3.1 植被光譜指標(biāo)的篩選

土壤中錳含量過少,植被葉綠體結(jié)構(gòu)將受到破壞,葉綠素含量隨之下降,將會導(dǎo)致葉片灰綠色或出現(xiàn)斑點等[4]。 錳含量過高將阻礙植物對鐵、 鈣、 鉬的吸收,經(jīng)常出現(xiàn)缺鉬癥狀。 葉片出現(xiàn)褐色斑點,葉緣白化或變紫,幼葉卷曲等。 大量研究表明植被光譜指標(biāo)能夠表征植被生長健康狀況[5]。 因此選用了11個植被光譜指標(biāo),其計算公式分別如下:

增強植被指數(shù)(enhanced vegetation index,EVI)

2.5×(B5-B4)/(B5+6B4-7.5B2-1)

(1)

比值植被指數(shù)(ratio vegetation index,RVI)[6]

B5/B4

(2)

綠度總和指數(shù)(sum green index,SG)

B5/B3

(3)

三角植被指數(shù)(triangular vegetation index,TVI)[7]

60×(B5-B3)-100×(B4-B3)

(4)

可見光大氣阻抗植被指數(shù)(visible atmospherically resistant index green,VARI)[6]

(B3-B4)/(B3+B4-B2)

(5)

歸一化植被指數(shù)(normailized difference vegetation index,NDVI)[8]

(B5-B4)/(B5+B4)

(6)

綠度歸一化植被指數(shù)(green NDVI,GNDVI)[9]

(B5-B3)/(B5+B3)

(7)

藍(lán)度歸一化植被指數(shù)(blue NDVI,BNDVI)

(B5-B2)/(B5+B2)

(8)

藍(lán)-綠通道歸一化植被指數(shù)(green-blue NDVI)[10]

(B5-B3-B2)/(B5+B3+B2)

(9)

紅-綠歸一化植被指數(shù)(green-red NDVI)[11]

(B5-B3-B4)/(B5+B3+B4)

(10)

葉綠素指數(shù)(modified chlorophyll reflectance index Green,mCRIG)[12]

[(B2)-1-(B3)-1]×B5

(11)

式中B2,B3,B4,B5分別為藍(lán)波段、 綠波段、 紅波段和近紅外波段反射率。

為了獲取最優(yōu)植被光譜響應(yīng)指標(biāo),利用皮爾遜相關(guān)系數(shù)對11個植被指數(shù)進(jìn)行篩選,皮爾遜相關(guān)系數(shù)表達(dá)式如式(12)

(12)

1.3.2 土壤錳含量最佳反演模型的篩選

為了獲取土壤錳含量最佳反演模型,基于最佳植被光譜響應(yīng)指標(biāo),利用偏最小二乘法(partial least-squares regression,PLSR)、 多元逐步回歸(multiple stepwise regression,MSR)和BP神經(jīng)網(wǎng)絡(luò)(back propagation neural network,BPNN)構(gòu)建植被光譜響應(yīng)指標(biāo)和土壤錳元素之間的關(guān)系模型,比較分析三個模型的精度,確定最佳反演模型。

(1)偏最小二乘法

偏最小二乘回歸(PLSR)綜合了典型相關(guān)分析、 多元線性回歸和主成分分析三種方法,可以用于解決很多普通多元回歸無法解決的問題,具有計算簡便、 預(yù)測能力強、 模型穩(wěn)健等特點[14]。 該方法的最簡形式是用線性模型來描述預(yù)測變量x和獨立變量y的關(guān)系,如式(13)所示。

Y=Xβ+ε

(13)

式(13)中,Y為經(jīng)過歸一化的因變量(土壤錳含量);X為經(jīng)過歸一化處理的自變量(最佳植被光譜指標(biāo));β為系數(shù)矩陣;ε為殘差矩陣。

(2)多元逐步回歸

多元逐步回歸(MSR)是回歸分析中一種篩選變量的過程,結(jié)合了前向選擇和后向剔除的優(yōu)點,其本質(zhì)是建立反映變量之間變化關(guān)系的最優(yōu)多元回歸模型。 采用逐步回歸法來建立植被指數(shù)與土壤錳元素含量的統(tǒng)計回歸模型[15]。 回歸方程如式(14)

Y=b0+b1x1+b2x2+…+bpxp+ε

(14)

式(14)中,Y為經(jīng)過歸一化的因變量(土壤錳含量);Xp為經(jīng)過歸一化處理的自變量(最佳植被光譜指標(biāo));b0為截距;bi為回歸系數(shù);ε為誤差項。

(3)BP神經(jīng)網(wǎng)絡(luò)

BP神經(jīng)網(wǎng)絡(luò)是一種由輸入層(input layer)、 隱含層(hide layer)和輸出層(output) 構(gòu)成的前饋式誤差反向傳播神經(jīng)網(wǎng)絡(luò),該網(wǎng)絡(luò)的特點是神經(jīng)信號向前傳遞,誤差沿著網(wǎng)絡(luò)反向傳播,該網(wǎng)絡(luò)包含了神經(jīng)網(wǎng)絡(luò)的精髓內(nèi)容,可以實現(xiàn)由任意的m維到n維的映射關(guān)系。 BP算法的中心思想是調(diào)整權(quán)值使網(wǎng)絡(luò)總誤差最小,通過把學(xué)習(xí)的結(jié)果反饋到中間層次的隱含層單元,改變它們的權(quán)系數(shù)矩陣,從而達(dá)到預(yù)期的學(xué)習(xí)目的[16]。 典型的BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)如圖2。

圖2 典型BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.2 Typical BP neural network structure

1.3.3 精度驗證

參考現(xiàn)有的土壤屬性含量光譜估算的評價標(biāo)準(zhǔn),采用均方根誤差(RMSE)、 相對分析誤差(RPD)和決定系數(shù)(R2)3個評價指標(biāo)評估模型反演效果[17-18]。 決定系數(shù)R2用于檢驗回歸方程對樣本預(yù)測值的擬合程度[式(15)],RMSE用于對模型的穩(wěn)定性進(jìn)行評價[式(16)],RPD能夠反映回歸模型的可靠性[式(17)、 式(18)]。

(15)

(16)

(17)

(18)

2 結(jié)果與討論

2.1 土壤錳的最佳光譜響應(yīng)指標(biāo)

利用Landsat8數(shù)據(jù),按照表1中的公式獲取11種植被光譜響應(yīng)指標(biāo)(RVI、 GVI、 mCRIG、 GBNDVI、 GRNDVI、 GNDVI、 BNDVI、 NDVI、 VARI、 EVI、 TVI),結(jié)果如圖3。

圖3 植被光譜指標(biāo)Fig.3 Vegetation spectral indices

利用IBM SPSS Statistics26對南川區(qū)森林覆蓋區(qū)域進(jìn)行土壤錳和11種植被光譜響應(yīng)指標(biāo)的相關(guān)性分析,結(jié)果如圖4所示,選取相關(guān)性大于0.5且相關(guān)性顯著的植被光譜響應(yīng)指標(biāo)進(jìn)行多重共線性診斷。 獲得3個與土壤錳元素有較高相關(guān)性且VIF均小于10的植被光譜響應(yīng)指標(biāo)(如表2)。 因此,本文選取RVI、 GRNDVI和VARI為最佳光譜響應(yīng)指標(biāo)。

表2 最佳植被光譜響應(yīng)指標(biāo)Table 2 Best vegetation spectral response index

圖4 土壤錳與植被光譜指標(biāo)相關(guān)性Fig.4 Correlation between soil manganese and vegetation spectral indexes

2.2 土壤錳最佳反演模型

將最佳植被光譜響應(yīng)指標(biāo)作為自變量,南川區(qū)土壤錳含量作為因變量,利用PLSR、 MSR和BP神經(jīng)網(wǎng)絡(luò)構(gòu)建兩者的關(guān)系模型。 基于PLSR的錳含量反演模型為:y=282.69+419.07RVI-5263.27GRNDVI-752.27VARI(R2=0.66),結(jié)果如圖5(a)。 基于MSR的錳含量反演模型為:y=272.776+417.174RVI-5216.026GRNDVI-748.28VARI(R2=0.66),結(jié)果如圖5(b)。 基于BPNN的錳含量反演模型中設(shè)置了3個隱含層和10個神經(jīng)元節(jié)點網(wǎng)絡(luò)迭代次數(shù)設(shè)置為2 000,學(xué)習(xí)率為0.01,學(xué)習(xí)目標(biāo)為0.01; 輸入層設(shè)置共3個,包括GRNDVI、 RVI和VARI; 輸出層為土壤錳含量,這里設(shè)為1; 模型的建立和運行在MATLAB R2019a上完成,結(jié)果如圖5(c)。

圖5 建模集實測值與估測值散點圖Fig.5 Scatter plots of measured and estimated values of modeling set

由圖5可知,PLSR模型和SMR模型的建模結(jié)果較為相近,在土壤錳含量0～1 500 mg·kg-1區(qū)間具有較好的擬合能力,預(yù)測點能夠較為均勻分布在1∶1線,而對于錳含量1 500以上的預(yù)測偏差值較大。 BPNN模型相對于前兩種模型擬合程度更高,預(yù)測效果最優(yōu),在土壤錳含量0～4 000 mg·kg-1區(qū)間都具有較好的擬合能力。

三種模型的測試精度如圖6所示。 在測試過程中,PLSR模型和MSR模型在土壤錳含量0～1 500 mg·kg-1區(qū)間具有較好的擬合能力,預(yù)測點能夠較為均勻分布在1∶1線,而對于錳含量1 500以上的預(yù)測偏差值較大。 總體而言,基于BPNN模型所建立的反演模型表現(xiàn)出最優(yōu)的預(yù)測效果,BPNN模擬的重金屬含量與實測值比較更接近于1∶1直線。 其中,BPNN模型建模的R2為0.78; 測試的R2為0.71; PLSR模型和MSR模型的擬合精度較為接近,建模R2均小于0.7。 已有研究證明,BPNN在處理復(fù)雜的非線性建模問題中有較好的擬合能力,而PLSR和MSR通常被用于確定自變量和因變量之間的線性關(guān)系。 本研究中,非線性模型(BPNN)的驗證精度明顯高于其他兩種線性模型(PLSR和MSR),其原因可能是土壤錳含量和部分植被指數(shù)之間存在明顯的非線性的關(guān)系。

圖6 測試集實測值與估測值散點圖Fig.6 Scatter plots of measured and estimated values of test set

2.3 土壤錳含量空間分布

由上述結(jié)果可知,BPNN模型明顯優(yōu)于其他模型,可以較好地反映土壤錳與最佳植被指數(shù)之間的映射關(guān)系。 因此選用BPNN模型進(jìn)行土壤錳空間制圖,其結(jié)果見圖7。

圖7 南川區(qū)土壤錳含量空間分布Fig.7 Spatial distribution of soil manganese content in Nanchuan District

土壤錳含量反演能夠充分體現(xiàn)出區(qū)域每個像元位置的土壤錳含量的變化特征。 從土壤錳含量空間分布來看,研究區(qū)的東北—西南方向出現(xiàn)錳含量偏高的條帶,土壤錳含量相對偏高的區(qū)域,是金佛山國家級自然保護(hù)區(qū),該地區(qū)可能存在較為原始的未開發(fā)錳礦區(qū),這也與樣本點實測值的分布情況基本一致。

為了驗證BPNN在區(qū)域尺度上估算土壤錳含量的可行性,使用139個樣點作為驗證集,對比分析其實測值與估算值的離散程度,結(jié)果如圖8所示,可以看出,錳含量的實測值與預(yù)測值的R2、 RMSE和RPD分別為0.69,567.64,1.30,制圖效果較佳,可以基本反映土壤錳含量的分布特征。 說明利用BPNN模型進(jìn)行土壤錳空間制圖具有一定的可行性。

圖8 驗證集實測值與估測值散點圖Fig.8 Scatter plots of measured and estimated values of validation set

2.4 討 論

利用遙感影像數(shù)據(jù)進(jìn)行土壤錳含量的反演有助于高效監(jiān)測區(qū)域尺度土壤錳含量。 當(dāng)前已有研究利用土壤光譜數(shù)據(jù)直接對土壤錳含量進(jìn)行反演,并且反演精度較高(R2約為0.70～0.78),可以滿足土壤錳的監(jiān)測需求。 但在植被常年覆蓋的南方地區(qū),難以從衛(wèi)星影像中獲取土壤光譜,此方法難以執(zhí)行。

本研究引入植被光譜響應(yīng)指標(biāo),利用線性(PLSR和MSR)和非線性(BPNN)方法構(gòu)建土壤錳含量反演模型,與PLSR和MSR相比,BPNN能顯著提高土壤錳含量的估算精度,非線性回歸模型在大面積反演土壤錳含量方面具有很大的潛力,這表明土壤錳和植被指數(shù)之間的關(guān)系是較為復(fù)雜的,不能用簡單的線性關(guān)系來表達(dá)。 有學(xué)者研究土壤光譜和錳含量關(guān)聯(lián)時,發(fā)現(xiàn)在400～900 nm波段對土壤錳含量有明顯響應(yīng),會隨著錳濃度增加出現(xiàn)依次降低的趨勢,特別是紅波段響應(yīng)最為明顯,而在綠波段其反射率變化不太大。 這與我們利用植被光譜進(jìn)行土壤錳含量反演研究比較一致。

在利用植被光譜響應(yīng)指標(biāo)反演土壤錳含量時,樣點建模和空間制圖均表現(xiàn)出不錯的精度,說明在區(qū)域尺度進(jìn)行大面積反演土壤錳含量是可行的。 但受限于Landsat8數(shù)據(jù)具有光譜范圍較寬、 波段數(shù)量較少等缺點,反演精度偏低。 在未來研究中,將嘗試?yán)酶吖庾V影像(如Hyperion影像等)反演土壤錳含量。 此外,僅利用BPNN一種非線性方法反演土壤錳含量,這可能會導(dǎo)致植被光譜響應(yīng)指標(biāo)與土壤錳含量響應(yīng)關(guān)系未能被準(zhǔn)確刻畫,未來引入更多非線性算法(隨機森林回歸算法、 決策樹回歸算法等)以期提高反演精度。

3 結(jié) 論

引入了植被光譜響應(yīng)指標(biāo),以重慶市南川區(qū)694個地理數(shù)據(jù)庫樣點為數(shù)據(jù)源,利用BPNN算法進(jìn)行土壤錳含量的預(yù)測,構(gòu)建了土壤錳含量最佳反演模型,反演重慶市南川區(qū)的土壤錳含量。 除EVI、 TVI外,其余9種植被光譜響應(yīng)指標(biāo)均與土壤錳有較高相關(guān)性,說明利用植被光譜響應(yīng)指標(biāo)與土壤錳含量具有較強的關(guān)聯(lián)性。 土壤錳含量的BPNN估算精度R2、 RMSE和RPD分別為0.78,334.24,2.13,制圖驗證點精度R2、 RMSE和RPD分別為0.69,567.64,1.30。 結(jié)果表明利用植被光譜響應(yīng)指標(biāo)反演土壤錳含量具有一定的可行性。