李 明， 秦 凱， 趙寧博， 田 豐， 趙英俊

核工業(yè)北京地質(zhì)研究院遙感信息與圖像分析技術(shù)國家級(jí)重點(diǎn)實(shí)驗(yàn)室， 北京 100029

引 言

土壤是植被生存的基礎(chǔ)， 即是人類賴以生存的基礎(chǔ)。 土壤中含有大量的養(yǎng)分， 這些養(yǎng)分能直接或經(jīng)轉(zhuǎn)化后被植物根系吸收的礦質(zhì)營養(yǎng)成分， 一般包括氮、 磷、 鉀、 鈣、 鎂、 硫、 鐵、 硼、 鉬、 鋅、 錳、 銅和氯等元素[1]， 與傳統(tǒng)的土壤元素化學(xué)測(cè)定方法相比， 航空及航天高光譜能夠高效地對(duì)大面積地區(qū)進(jìn)行檢測(cè)， 其通過定量研究不同波段光譜與元素含量的關(guān)系進(jìn)行建模反演， 能夠獲取大面積研究區(qū)某元素的含量分布圖。 國內(nèi)外學(xué)者在可見-近紅外波段對(duì)土壤的水分[2-3]、 有機(jī)質(zhì)[4-5]、 含鹽量[6-8]、 養(yǎng)分元素[9-11]等的定量反演做了大量研究工作， 這一波段范圍已經(jīng)無法滿足諸多學(xué)者的研究熱情， 因此中紅外波段成為了新的土壤理化參數(shù)定量反演研究熱點(diǎn)。 阿爾達(dá)克·克里木等[12]利用4種變換后的熱紅外發(fā)射率光譜特征通過多元回歸的方法建立其與土壤含鹽量的模型； 夏軍等[13]通過分析125個(gè)土壤樣本熱紅外波段與土壤含鹽量的相關(guān)性以及通過偏最小二乘與逐步回歸建模， 得到土壤熱紅外發(fā)射率估算鹽分含量的模型， 并得出偏最小二乘法定量反演土壤含鹽量效果比逐步回歸好這一結(jié)論， 模型預(yù)測(cè)的R2達(dá)到0.958， RMSE為1.911%； 楊永民等[14]基于熱紅外數(shù)據(jù)， 對(duì)土壤含水量四種估算方法進(jìn)行對(duì)比分析， 使用ASTER數(shù)據(jù)估算了黑河流域中游地區(qū)的土壤水分狀況。 由此可見， 發(fā)射率數(shù)據(jù)與土壤眾多理化指標(biāo)有較為緊密的關(guān)系。 而發(fā)射率數(shù)據(jù)對(duì)于土壤K元素含量的反演研究較少且模型精度較差， 本文旨在探究提高模型精度的新方法。

1 實(shí)驗(yàn)部分

1.1 數(shù)據(jù)采集與處理

熱紅外航空成像光譜儀(thermal airborne hyperspectral imager， TASI)是加拿大研制的先進(jìn)的機(jī)載熱紅外高光譜設(shè)備。 該設(shè)備在8～11.5 μm范圍內(nèi)有32波段， 波段間隔為0.109 5 μm， 半高寬為0.054 8 μm， 總視場(chǎng)角為40°(表1)。 本文利用核工業(yè)北京地質(zhì)研究院遙感信息與圖像分析技術(shù)國家級(jí)重點(diǎn)實(shí)驗(yàn)室的TASI航空高光譜成像系統(tǒng)， 在東北黑龍江省海倫地區(qū)附近獲取了高空間分辨率的高光譜熱紅外遙感數(shù)據(jù)， 在經(jīng)過大氣校正等數(shù)據(jù)預(yù)處理后， 采用楊杭等[15]改進(jìn)的TES(temperature-emissivity separation)分離方法進(jìn)行溫度與發(fā)射率的分離， 獲取研究區(qū)發(fā)射率數(shù)據(jù)。

表1 TASI各通道中心波長

土壤樣品采集地點(diǎn)為東北黑龍江省海倫地區(qū)， 該地區(qū)土壤類型主要為黑土， 土壤腐殖質(zhì)層較厚， 約30～60 cm， 有機(jī)質(zhì)含量在2.5%～4.5%之間， 粘粒含量在40%～60%之間， 屬于粘土， 西南部旱田土壤基本為中性， 中部及東北部旱田土壤為酸性。 在工作區(qū)不同位置共采集土壤樣本40個(gè)(圖1)。 測(cè)區(qū)表層為黑色腐殖質(zhì)層， 當(dāng)天同步飛行采集表層0～20 cm的土樣， 剔除大的植物殘茬、 石礪等雜物， 置于實(shí)驗(yàn)室風(fēng)干研磨， 過0.15 mm篩選用于含量測(cè)定。 全鉀含量采用X射線熒光光譜法測(cè)定(表2)。

圖1 TASI數(shù)據(jù)采集區(qū)域及采樣點(diǎn)分布圖

表2 土壤樣本鉀(K)元素含量信息表

1.2 建模方法及實(shí)現(xiàn)

探究多個(gè)自變量與因變量之間關(guān)系的方法很多， 對(duì)于光譜數(shù)據(jù)， 常用的建模方法有最小二乘法、 多元逐步回歸擬合、 BP神經(jīng)網(wǎng)絡(luò)、 支持向量機(jī)SVM等。 機(jī)器學(xué)習(xí)類方法雖然能較好地訓(xùn)練模型并進(jìn)行預(yù)測(cè)， 但是存在較多的人為干預(yù)調(diào)參問題， 同時(shí)對(duì)于某一元素在研究區(qū)整體的提取具有較大難度。 本研究主要聚焦于偏最小二乘法與多元逐步回歸擬合， 尤其利用多元逐步回歸方法研究發(fā)射率與K元素含量關(guān)系時(shí)創(chuàng)新性地使用了全二次逐步回歸進(jìn)一步提高模型的精度。

相比常規(guī)的n元線性逐步回歸僅有的常數(shù)項(xiàng)和線性項(xiàng)， 多元全二次逐步回歸引入了交叉乘積項(xiàng)和平方項(xiàng)進(jìn)行回歸， 以對(duì)回歸方程中常數(shù)項(xiàng)、 線性項(xiàng)和二次項(xiàng)進(jìn)行的t檢驗(yàn)的p值是否小于等于0.05為判定依據(jù)， 依次引進(jìn)顯著項(xiàng)剔除非顯著項(xiàng)， 可以有效地解決自變量的多重共線性問題。 但參數(shù)不宜引入過多以免數(shù)據(jù)產(chǎn)生過擬合， 其中， 引入了4個(gè)參數(shù)進(jìn)行全二次多元逐步回歸方程如下所示， 由于增加了更多的系數(shù)， 因此能夠更加精確地進(jìn)行回歸模型的建立。 本次研究所有建模方法均利用Matlab編程實(shí)現(xiàn)。 TASI影像數(shù)據(jù)預(yù)處理及溫度與發(fā)射率分離利用ENVI5.3-IDL編程實(shí)現(xiàn)。

式中，yi為響應(yīng)變量(預(yù)測(cè)值)；b0，bi，bij和bii分別為回歸方程常數(shù)項(xiàng)、 線性項(xiàng)、 交叉乘積項(xiàng)、 平方項(xiàng)的系數(shù)；xi為預(yù)測(cè)變量(輸入值)， 本文中為輸入的四個(gè)相關(guān)性強(qiáng)波段值；εi～N表示數(shù)據(jù)服從期望值μ為0的正態(tài)分布。

2 結(jié)果與討論

2.1 土壤熱紅外發(fā)射率光譜曲線特征

物體的發(fā)射率除了取決于其材質(zhì)， 更取決于其所存在的環(huán)境。 對(duì)于土壤， 其熱紅外發(fā)射率光譜曲線形態(tài)主要取決于土壤中所含的礦物種類以及含量、 水分、 有機(jī)質(zhì)含量和溫度等因素。 研究所采集的40個(gè)土壤樣本的發(fā)射率光譜曲線在8～9.6 μm變化趨勢(shì)基本一致(圖2)， 9.6～11.45 μm后形態(tài)有所變化但基本變化趨于平緩， 發(fā)射率值在8.38， 8.6以及9.26 μm處出現(xiàn)了三處非常明顯的波谷， 其中8.38與9.3 μm的兩處波谷出現(xiàn)了類似石英的波譜特征。 夏軍等研究發(fā)現(xiàn)， 土壤中硅酸鹽礦物導(dǎo)致發(fā)射率光譜曲線呈現(xiàn)明顯的Reststrahlen吸收特征， 即不對(duì)稱雙吸收谷， 兩個(gè)吸收谷分別位于8.23和9.27 μm波長附近， 且后一個(gè)吸收谷較深， 寬度較大[13]。 9.6 μm后光譜曲線出現(xiàn)不同波動(dòng)是由于土壤中所含不同礦物的基團(tuán)內(nèi)部振動(dòng)產(chǎn)生的譜帶不同所致， 但總體幅度變化較小。

圖2 40個(gè)土壤樣本的熱紅外發(fā)射率光譜曲線

2.2 土壤元素含量與發(fā)射率相關(guān)性分析

通過40個(gè)土壤樣本K元素含量與發(fā)射率各波段做皮爾森相關(guān)性分析可以看出， K元素與32個(gè)波段發(fā)射率呈負(fù)相關(guān)關(guān)系， 其中相關(guān)系數(shù)大于0.6呈強(qiáng)相關(guān)的波段有5個(gè)， 分別是6， 11， 15， 22和23波段； 相關(guān)系數(shù)介于0.4～0.6中等相關(guān)的波段有17個(gè)(圖3)， 整體具有較強(qiáng)的相關(guān)性。

圖3 鉀元素含量與發(fā)射率相關(guān)性

由于多元全二次逐步回歸引入的參數(shù)數(shù)量直接決定了模型建立引入系數(shù)的多少和復(fù)雜程度， 因此， 為了避免建模時(shí)系數(shù)過多導(dǎo)致過擬合增大無謂的計(jì)算量， 初步選取與K元素含量相關(guān)系數(shù)最高的4個(gè)波段用于模型的建立。 土壤中鉀元素有多種賦存形態(tài)， 大部分以原生或次生的結(jié)晶硅酸鹽狀態(tài)存在于土壤中， 其中云母族礦物參考(白云母、 黑云母)及富鉀長石(正長石)中含鉀元素最多， 白云母、 黑云母以及正長石特征吸收位置選擇6， 11， 15和23波段作為特征波段(表3)。

表3 土壤元素含量與發(fā)射率顯著相關(guān)所對(duì)應(yīng)的TASI波段

2.3 光譜定量預(yù)測(cè)

將40個(gè)土壤樣本隨機(jī)分為兩組， 其中32個(gè)樣本用于含量預(yù)測(cè)模型的建立， 8個(gè)樣本用來測(cè)試模型的精度。 驗(yàn)證K元素含量數(shù)據(jù)符合正態(tài)分布后， 以所選4個(gè)特征波段發(fā)射率數(shù)據(jù)作為自變量， K元素含量為因變量， 以對(duì)回歸方程中常數(shù)項(xiàng)、 線性項(xiàng)和二次項(xiàng)進(jìn)行的t檢驗(yàn)的p值是否小于等于顯著性水平0.05為判定依據(jù)， 依次引進(jìn)顯著項(xiàng)剔除非顯著項(xiàng)， 同時(shí)對(duì)模型總體進(jìn)行F檢驗(yàn)的p值是否小于等于顯著性水平0.05來驗(yàn)證樣本觀測(cè)值與總體假設(shè)值是否存在顯著性差異從而建立模型[圖4(a,b)]。

圖4 逐步回歸模型預(yù)測(cè)值與真實(shí)值擬合效果圖

常規(guī)逐步回歸建立的回歸擬合模型均方根誤差RMSE為0.031， 調(diào)整后的判定系數(shù)R2為0.569， 測(cè)試集的均方根誤差RMSE為0.031， 調(diào)整后的判定系數(shù)R2為0.78[見圖5(a)]； 全二次多元逐步回歸建立的回歸擬合模型均方根誤差RMSE為0.027， 調(diào)整后的判定系數(shù)R2為0.667， 測(cè)試集的均方根誤差RMSE為0.032， 調(diào)整后的判定系數(shù)R2為0.82[見圖5(b)]， 所有指標(biāo)均通過p值小于0.05的顯著性驗(yàn)證(表4)。 通過模型擬合結(jié)果及評(píng)價(jià)指標(biāo)來看， 全二次多元逐步回歸比常規(guī)多元逐步回歸建模精度以及驗(yàn)證精度均有所提高。

表4 逐步回歸模型建模結(jié)果

為進(jìn)一步提高建模精度， 利用Matlab回歸診斷的學(xué)生化殘差來進(jìn)行模型改進(jìn)。 通過|Sei|>2來查找遠(yuǎn)離數(shù)據(jù)集中心觀測(cè)點(diǎn)即異常點(diǎn)， 剔除異常點(diǎn)來進(jìn)一步提高模型精度。 同時(shí)， 以32個(gè)波段發(fā)射率數(shù)據(jù)為自變量， K元素含量為因變量進(jìn)行偏最小二乘法建模， 進(jìn)一步對(duì)比三種模型的優(yōu)劣。

(2)

偏最小二乘法建立的回歸擬合模型入選主成分?jǐn)?shù)為2， 均方根誤差RMSE為0.033， 判定系數(shù)R2為0.45， 測(cè)試集的均方根誤差RMSE為0.037， 判定系數(shù)R2為0.51(圖6)。

圖5 鉀元素含量實(shí)測(cè)值與預(yù)測(cè)值散點(diǎn)圖

圖6 鉀元素含量實(shí)測(cè)值與預(yù)測(cè)值散點(diǎn)圖

通過評(píng)價(jià)指標(biāo)分析改進(jìn)后的多元逐步回歸模型發(fā)現(xiàn)， 雖然建模精度有所提高， 但測(cè)試集的精度卻均有所下降。 與改進(jìn)前模型相比， 常規(guī)多元逐步回歸建模樣本的均方根誤差RMSE降低了0.7%， 判定系數(shù)R2提高了0.163； 測(cè)試集的均方根誤差RMSE提高了0.2%， 判定系數(shù)R2降低了0.015； 全二次多元逐步回歸均方根誤差RMSE降低了0.71%， 判定系數(shù)R2提高了0.135； 測(cè)試集的均方根誤差RMSE提高了0.2%， 判定系數(shù)R2降低了0.1， 同時(shí)由于剔除了某些不顯著的變量， 模型再次引入了新的變量， 參數(shù)從7個(gè)增加到了10個(gè)(表5)。 改進(jìn)后的模型訓(xùn)練集精度上升而測(cè)試集精度下降的原因應(yīng)該是訓(xùn)練集數(shù)據(jù)發(fā)生了過擬合， 從實(shí)驗(yàn)結(jié)果分析， 改進(jìn)前的模型泛化能力更強(qiáng)， 更適用于研究區(qū)K元素的反演。 同時(shí)對(duì)比偏最小二乘法建模， 全二次多元逐步回歸各項(xiàng)評(píng)價(jià)指標(biāo)均優(yōu)于其余兩種方法。 由于混合像元影響以及樣本的選擇有所差異， 模型整體擬合精度不是很高， 但本研究提出新的逐步回歸方法有效地提高了模型的精度。

表5 建模結(jié)果對(duì)比

相比常規(guī)多元逐步回歸僅考慮常數(shù)項(xiàng)和線性項(xiàng)， 全二次多元逐步回歸能夠引入更多的變量參與到回歸模型的建立中， 從而提高模型的反演精度。 說明利用TASI數(shù)據(jù)的相關(guān)波段通過全二次多元逐步回歸方法反演元素含量是可行的， 比起傳統(tǒng)的化學(xué)填圖， 遙感反演的方法在損失部分精度的條件下能夠高效大面積地反演某個(gè)地區(qū)元素含量。

3 結(jié) 論

針對(duì)土壤中K元素含量反演， 利用熱紅外航空成像光譜儀TASI數(shù)據(jù)的發(fā)射率數(shù)據(jù)， 創(chuàng)新性地使用了一種新的逐步回歸方法-“全二次多元逐步回歸”建立模型， 相對(duì)于常規(guī)多元逐步回歸， 引入了更多的參數(shù)進(jìn)行模型的建立， 能夠有效提高反演精度。 通過研究發(fā)現(xiàn)， 土壤發(fā)射率數(shù)據(jù)對(duì)于選用有效特征波段對(duì)K元素具有較高的反演精度。 K元素通過多元逐步回歸建模與預(yù)測(cè)的均方根誤差RMSE： 0.027和0.032， 判定系數(shù)R2： 0.667和0.82， 相比于常規(guī)多元逐步回歸建模與預(yù)測(cè)的均方根誤差RMSE： 0.031和0.031， 判定系數(shù)R2： 0.569和0.78與偏最小二乘法建模與預(yù)測(cè)的均方根誤差RMSE： 0.033和0.037， 判定系數(shù)R2： 0.45和0.51評(píng)價(jià)指標(biāo)精度均有所提高， 說明本方法有效提高了利用發(fā)射率數(shù)據(jù)對(duì)K元素的反演精度。