王洪偉，王 波，紀(jì) 童，徐 君，劇 鋒，王彩玲

1. 武警工程大學(xué)，陜西 西安 710086 2. 鹽池縣草原實(shí)驗(yàn)站，寧夏 鹽池 751506 3. 甘肅農(nóng)業(yè)大學(xué)草業(yè)學(xué)院，甘肅 蘭州 730070 4. 西安航空學(xué)院，陜西 西安 710077 5. 中華人民共和國(guó)銀川海關(guān)，寧夏 銀川 750000 6. 西安石油大學(xué)，陜西 西安 710065

引 言

隨著人類物質(zhì)生活水平的提高和工業(yè)化的發(fā)展，水污染已經(jīng)成為當(dāng)今社會(huì)普遍存在的問題，其監(jiān)測(cè)與治理也備受關(guān)注。生化需氧量(biochemical oxygen demand,BOD)是水體中的好氧微生物在一定溫度下將水中有機(jī)物分解成無機(jī)質(zhì)，這一特定時(shí)間內(nèi)的氧化過程中所需要的溶解氧量，是監(jiān)測(cè)水中有機(jī)物染物的一個(gè)綜合指標(biāo)[1]，是地表水、生活污水及絕大多數(shù)工業(yè)廢水的必測(cè)指標(biāo)之一。BOD值越高表明水中溶解氧會(huì)被自身微生物消耗的數(shù)值越高，造成許多的生態(tài)問題[2]?！拔迦张囵B(yǎng)法”為現(xiàn)下普遍的接受測(cè)定BOD的方法，但測(cè)定時(shí)間長(zhǎng)、不能及時(shí)反映水質(zhì)變化，不適合現(xiàn)場(chǎng)監(jiān)測(cè)。

自20世紀(jì)70年代以來，隨著遙感技術(shù)的快速發(fā)展，高光譜技術(shù)已成為現(xiàn)代遙感技術(shù)的重要組成部分[3]，利用高光譜技術(shù)反演水質(zhì)指數(shù)早有研究，劉彥君等[4]利用多光譜數(shù)據(jù)，進(jìn)行線性與非線性模型反演研究，對(duì)浙江農(nóng)林大學(xué)東湖水體的總磷(TP)、濁度(SS)、懸浮物濃度(TUB)進(jìn)行了反演。林劍遠(yuǎn)等[5]利用水質(zhì)化驗(yàn)數(shù)據(jù)和光譜反射率進(jìn)行相關(guān)性分析，建立了浙江省嘉興市河網(wǎng)化學(xué)需氧量(CODcr)、生化需氧量(BOD5)、總磷(TP)、總氮(TN)的反演模型。周亞東等[6]利用GF-1號(hào)WFV遙感影像，通過多元線性回歸和RBF神經(jīng)網(wǎng)絡(luò)模型建立了武漢市周圍水域綜合營(yíng)養(yǎng)狀態(tài)指數(shù)模型。這些成果有效解讀了水體光譜特征規(guī)律，為遙感監(jiān)測(cè)水質(zhì)，生產(chǎn)生活提供了理論支撐與技術(shù)指導(dǎo)。

原始光譜反射數(shù)據(jù)有著數(shù)據(jù)量大，指標(biāo)彼此高度相關(guān)的特性； 原始指標(biāo)高度相關(guān)的特性經(jīng)常會(huì)導(dǎo)致多重共線性問題的產(chǎn)生，從而導(dǎo)致模型失真[7]，因此如何對(duì)大量光譜數(shù)據(jù)進(jìn)行處理和挑選一直是光譜反演模型的重點(diǎn)。主成分分析法(PCA)與偏最小二乘法(PLS)作為常用降維方法在遙感上應(yīng)用廣泛[8]，許多研究結(jié)果也表明應(yīng)用主成分分析與偏最小二乘法篩選的主成分參數(shù)可以更好的反演各自的指標(biāo)。楊國(guó)范等[9]利用比值線性回歸模型與最小二乘支持向量機(jī)，對(duì)鐵嶺清水河庫(kù)葉綠素a濃度與Landast OLI衛(wèi)星數(shù)據(jù)分析，并建立了葉綠素濃度a的反演模型。何金成等[10]利用近紅外光譜數(shù)據(jù)結(jié)合偏最小二乘法回歸建立了BOD預(yù)測(cè)模型。

現(xiàn)有文獻(xiàn)報(bào)道中，利用光譜估測(cè)水質(zhì)參數(shù)BOD指標(biāo)的報(bào)道較少，基于此試驗(yàn)利用光譜數(shù)據(jù)進(jìn)行水體指標(biāo)BOD的反演，測(cè)定水體樣本光譜數(shù)據(jù)的同時(shí)收集水體樣本并帶回實(shí)驗(yàn)室測(cè)定BOD指標(biāo)，將采集到的光譜數(shù)據(jù)與BOD指標(biāo)進(jìn)行Person相關(guān)性分析，挑選敏感光譜指標(biāo)； 由于光譜指標(biāo)之間的高度相關(guān)，為避免模型失真，在建立反演模型之前，利用主成分分析和偏最小二乘法分別對(duì)光譜指標(biāo)進(jìn)行處理，消除指標(biāo)之間的多重共線性問題，最終建立多元線性回歸模型與偏最小二乘模型，比較兩種建模方法的建模精度與預(yù)測(cè)效果，選出更加適合反演BOD指標(biāo)的建模方法。探索利用高光譜技術(shù)估測(cè)水體BOD值的可行性與最優(yōu)方法，為實(shí)時(shí)診斷水體狀況提供理論基礎(chǔ)和關(guān)鍵技術(shù)，為實(shí)現(xiàn)對(duì)BOD指標(biāo)實(shí)時(shí)監(jiān)測(cè)提供可行的途徑。

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

1.1 試驗(yàn)地概況

于2018年對(duì)西安市地表水環(huán)境進(jìn)行取樣研究，取樣地點(diǎn)集中于渭河(林家村)、浐河(田家灣)、灞河(馬渡王)，共計(jì)60處采樣點(diǎn)，每處采樣點(diǎn)共計(jì)10次重復(fù)。

1.2 方法

1.2.1 光譜數(shù)據(jù)獲取與校正

所用儀器為美國(guó)ASD (Analytica Spectra Devices.,Inc)公司制造的適用于遙感測(cè)量、農(nóng)作物監(jiān)測(cè)等方面的 FieldSpec?4 Hi-ResASD便攜式地物光譜儀，其光譜范圍為300～2 500 nm。

光譜采集選擇干燥、無風(fēng)、晴朗無云或少云的天氣進(jìn)行，并根據(jù)天氣條件及時(shí)進(jìn)行標(biāo)準(zhǔn)白板校正，采集時(shí)間盡量在10：00—14：00之間，此時(shí)光照條件良好。進(jìn)行地面水質(zhì)采樣和水體光譜數(shù)據(jù)等實(shí)驗(yàn)數(shù)據(jù)獲取，光譜采集參數(shù)設(shè)置時(shí)間為100 ms，測(cè)量后及時(shí)進(jìn)行白板校正[11]。每塊樣本選擇2～3個(gè)光譜采樣點(diǎn)進(jìn)行高光譜數(shù)據(jù)采集，每個(gè)樣點(diǎn)每次重復(fù)測(cè)量10次，最后以該樣點(diǎn)的光譜反射率均值制作光譜反射率曲線。

1.2.2 BOD指標(biāo)的測(cè)定

采集水體樣本時(shí)，利用聚乙烯桶采集距離水面10～12 cm的水樣，不使漂浮于水面的物質(zhì)混入，每處試驗(yàn)點(diǎn)共取10次樣本，對(duì)水樣加入保存劑，以便將樣本帶回實(shí)驗(yàn)室，利用標(biāo)準(zhǔn)稀釋法[12]處理水樣，并在20 ℃培養(yǎng)箱中培養(yǎng)，5 d后測(cè)出培養(yǎng)后的溶解氧含量，取平均值作為BOD指標(biāo)原始因變量。樣本BOD參數(shù)變化范圍如表1所示。

表1 水質(zhì)參數(shù)變化范圍Table 1 Variation range of water quality parameters

1.3 數(shù)據(jù)處理

普通的多元線性回歸應(yīng)用中有許多限制，最典型的問題就是自變量之間的多重相關(guān)性。為此，利用主成分分析降維與消除指標(biāo)間多重共線性的特性，篩選多元線性模型的自變量，已期解決多重共線性對(duì)參數(shù)估計(jì)的影響，減小模型誤差。偏最小二乘回歸中開辟了一種有效的技術(shù)途徑，通過對(duì)系統(tǒng)中的數(shù)據(jù)信息進(jìn)行分解和篩選，提取對(duì)因變量的解釋性最強(qiáng)的綜合變量，辨識(shí)系統(tǒng)中的信息與噪聲，從而能夠更好地克服變量多重相關(guān)性在系統(tǒng)建模中的不良作用。

2 結(jié)果與討論

2.1 Person相關(guān)性

圖1為原始光譜DN值與水體BOD含量的相關(guān)系數(shù)圖，因圖中波段1 023～2 500 nm與水體BOD含量無顯著相關(guān)性，因此圖中只展示了350～1 023 nm波段范圍內(nèi)的相關(guān)系數(shù)，由圖可知BOD指標(biāo)與光譜在350～900 nm呈負(fù)相關(guān)，960～100 nm為正相關(guān)，350～490與920～1 000 nm與BOD指標(biāo)無顯著相關(guān)性，BOD敏感波段大體分布于600～900 nm，其中758 nm處為相關(guān)系數(shù)絕對(duì)值最大值0.418，根據(jù)相關(guān)系數(shù)大小與顯著性原則，篩選出了35個(gè)與BOD指標(biāo)極顯著相關(guān)的原始光譜指標(biāo)，作為多元線性回歸模型與偏最小二乘模型的自變量，篩選指標(biāo)相關(guān)系數(shù)絕對(duì)值由大到小分別為： 758，759，853，809，1000，810，890，813，851，1 012，807，893，618，864，816，806，782，787，785，888，796，808，924，845，663，530，887，724，863，889，757，683，628，909和689 nm。主成分分析要求建模數(shù)據(jù)量高于變量數(shù)，偏最小二乘法允許在樣本點(diǎn)個(gè)數(shù)少于變量個(gè)數(shù)的條件下進(jìn)行回歸建模，Person相關(guān)系數(shù)法共篩選出35個(gè)光譜變量，因此將60組樣本數(shù)據(jù)分為建模組(40)與檢驗(yàn)組(20)。

2.2 主成分分析

主成分分析結(jié)果如圖2所示。

圖1 相關(guān)系數(shù)圖

圖2 主成分分析碎石圖注： 橫坐標(biāo)是主成分，縱坐標(biāo)為解釋程度Fig.2 Principal component analysis lithotripsy

經(jīng)分析共有10個(gè)主成分。其中主成分1方差貢獻(xiàn)率為94.9%，主成分2方差貢獻(xiàn)率為1%，而主成分3～10累積方差貢獻(xiàn)率不足10%，且主成分2到主成分3，斜率開始趨于平緩，因此剔除主成分3～10，只保留主成分1和2(Z1和Z2)，這2個(gè)主成分既能達(dá)到降維的目的，又能反映原始數(shù)據(jù)95.9%的信息。

圖3直觀展現(xiàn)了各植被指數(shù)在主成分1和主成分2中的分布情況。橫縱坐標(biāo)分別代表第一主成分與第二主成分以及各自的貢獻(xiàn)率，Z1和Z2累計(jì)貢獻(xiàn)率高達(dá)95.9%，可以解釋原有變量中的大部分信息，4個(gè)BOD含量分組中0～0.2與0.4～0.6 mol·L-1在4組中彼此獨(dú)立，可以明顯區(qū)分，0.2～0.4與0.6～0.8 mol·L-1彼此交叉分組不明顯。

將特征向量代入主成分公式中，得到主成分Z1和Z2的表達(dá)式

Z1=0.168x1+0.168x2+0.17x3+…+0.166x34+0.169x35

Z2=0.382x1+0.353x2+0.287x3+0.28x4+0.151x5+0.122x7-0.311x26-0.264x27-0.208x29-0.364x30-0.196x32-0.271x34-0.13x35

將2個(gè)主成分分別代入多元線性回歸中，得到的方程

YBOD=-0.000 004 468z1+0.000 059 19z2+9.217

(R2=0.656，RMSE=0.007)

多元回歸模型中BOD與主成分?jǐn)M合方程R2較大，RMSE值較小，說明利用主成分Z1和Z2通過多元線性回歸，可以很好的擬合水體BOD指標(biāo)。

圖3 主成分分析效果圖

2.3 偏最小二乘模型的構(gòu)建

由于自變量與因變量之間的量綱與數(shù)值都是不同的，現(xiàn)將BOD值與篩選的光譜指標(biāo)進(jìn)行標(biāo)準(zhǔn)化處理，利用R語(yǔ)言PLS偏最小二乘函數(shù)包建立水質(zhì)BOD含量的估測(cè)模型，各主成分貢獻(xiàn)率結(jié)果見表2。

表2 主成分貢獻(xiàn)率Table 2 Contribution rate of principal component

如表2可知，當(dāng)主成分為3時(shí)，解釋率逐漸趨于平穩(wěn)，因此選取comps=3時(shí)建立模型。y=0.015 703x1+0.124 092x2+0.423 545x3-0.181 04x4-0.255 47x5+…+0.331 165x34-0.189x35(R2=0.896，RMSEP=0.7469)。

使用函數(shù)包中jack.test函數(shù)對(duì)回歸系數(shù)進(jìn)行顯著性檢驗(yàn)結(jié)果見表3。

表3 jack.test函數(shù)顯著性檢驗(yàn)Table 3 Significance test of jack.test function

通過jack.test函數(shù)進(jìn)行顯著性檢驗(yàn)，表3中“*”代表極顯著影響，“**”代表顯著影響，x1—x35代表波長(zhǎng)按由小到大排列的原始光譜變量，由表3可知對(duì)水體BOD含量有顯著影響的光譜指標(biāo)有628，889和893 nm，其中對(duì)BOD有正向影響的光譜指標(biāo)為628與889 nm，對(duì)BOD有負(fù)向影響的光譜指標(biāo)為893 nm。

2.4 最優(yōu)模型篩選

比較多元線性回歸模型與偏最小二乘法模型，依據(jù)R2最大RMSE最小原則，最終采用偏最小二乘法模型y=0.015 703x1+0.124 092x2+0.423 545x3-0.181 04x4-0.255 47x5+…+0.331 165x34-0.189x35(R2=0.896，RMSEP=0.746 9)。

2.5 模型精度檢驗(yàn)

將檢驗(yàn)組的20組BOD與光譜數(shù)據(jù)代入模型進(jìn)行模型精度檢驗(yàn)，檢驗(yàn)結(jié)果見表4。

表4 模型精度檢驗(yàn)Table 4 Model accuracy test

從表4可以看出，偏最小二乘模型，其均方根誤差較低為0.12，且估測(cè)精度R2較高。說明利用偏最小二乘法可以建立精度較好的BOD反演模型。

3 結(jié) 論

利用多元線性回歸與偏最小二乘法建立水質(zhì)BOD指標(biāo)的反演模型。在進(jìn)行光譜單波段與BOD指標(biāo)相關(guān)性分析時(shí)發(fā)現(xiàn)在350～500 nm波段相關(guān)系數(shù)偏低，350～500 nm原始光譜反射曲線雜亂，說明該波段可能受其他水質(zhì)參數(shù)影響，波段敏感性較差，不能作為模型預(yù)測(cè)波段。隨著波長(zhǎng)增加相關(guān)系數(shù)于758 nm達(dá)到最高值，且光譜最優(yōu)反演波段大多分布在600～900 nm處，與林劍遠(yuǎn)等[5]得到的高光譜遙感數(shù)據(jù)與BOD指標(biāo)敏感波段750～900 nm有一定不同但也有相似之處，主要原因有以下幾點(diǎn)：

(1)水體光譜受自然條件與人為干擾，使邊緣波段噪聲很大，導(dǎo)致350～500 nm波段敏感性較差，與BOD指標(biāo)的相關(guān)性較低。

(2)水體光譜易受時(shí)間空間等影響，導(dǎo)致光譜區(qū)別較大，但光譜趨勢(shì)整體相似，且受其他水質(zhì)指標(biāo)影響，其光譜也會(huì)隨之變化。

高光譜具有分辨率高，波段連續(xù)性強(qiáng)的特點(diǎn)，但光譜信息冗雜，數(shù)據(jù)的篩選與模型的簡(jiǎn)化一直是光譜模型研究的重點(diǎn)[13]，通過主成分分析和偏最小二乘法綜合篩選的光譜指標(biāo)，建立了多元線性回歸模型與偏最小二乘法，結(jié)果表明主成分分析與偏最小二乘法可以有效降低數(shù)據(jù)維度，綜合篩選指標(biāo)特性，提高光譜數(shù)據(jù)與BOD參數(shù)的相關(guān)性與模型精度，其中偏最小二乘模型模型精度遠(yuǎn)高于多元線性回歸模型，因?yàn)槠钚《朔ㄊ欠謩e從因變量與自變量中提取成分因子，保證成分因子盡可能多的反應(yīng)變量的變異信息，同時(shí)也保證了兩者之間相關(guān)性最大[14]，且試驗(yàn)中樣本個(gè)數(shù)與變量個(gè)數(shù)大致一致，適用于偏最小二乘法模型。

在擬合偏最小二乘模型時(shí)，利用jack.test函數(shù)得出對(duì)水體BOD含量有顯著影響的光譜指標(biāo)有628，889與893 nm，說明628，889與893 nm可以作為反演BOD指標(biāo)的敏感波長(zhǎng)。林劍遠(yuǎn)等[5]以高光譜數(shù)據(jù)研究是城市河網(wǎng)BOD指標(biāo)中發(fā)現(xiàn)波段565 nm為單波段與BOD指標(biāo)相關(guān)系數(shù)(0.44)最佳波段，689/667 nm為組合波段與BOD指標(biāo)相關(guān)系數(shù)(0.84)最佳波段，與本工作篩選的敏感波段有所不同，但有所相近。

以上試驗(yàn)結(jié)果為水質(zhì)BOD指標(biāo)的快速估算提供了依據(jù)，也為水體質(zhì)量評(píng)估提供更便利的方案。