基于支持向量機(jī)的長(zhǎng)江口及其鄰近海域葉綠素a濃度預(yù)測(cè)模型?

李修竹， 蘇榮國(guó)??， 張傳松， 石曉勇,2

(1.中國(guó)海洋大學(xué)化學(xué)化工學(xué)院，山東 青島 266100； 2.國(guó)家海洋局海洋減災(zāi)中心, 北京 100194)

隨著沿海地區(qū)經(jīng)濟(jì)的飛速發(fā)展，大量工業(yè)廢水和污染物排放入海，導(dǎo)致近岸海域富營(yíng)養(yǎng)化嚴(yán)重，赤潮頻發(fā)。赤潮對(duì)生態(tài)環(huán)境和人類(lèi)健康造成巨大破壞，打破了海域生態(tài)系統(tǒng)平衡，并給水產(chǎn)養(yǎng)殖造成了巨大損失，引起了國(guó)內(nèi)外的廣泛關(guān)注[1-2]。因此，對(duì)水體富營(yíng)養(yǎng)化進(jìn)行早期預(yù)測(cè)很有必要。眾所周知，海洋生態(tài)系統(tǒng)十分復(fù)雜，各因素之間的相互作用及其動(dòng)態(tài)變化過(guò)程未被完全知曉，制約了傳統(tǒng)的生態(tài)水動(dòng)力模型的發(fā)展[3]。

近年來(lái)，數(shù)據(jù)智能技術(shù)在預(yù)測(cè)模型中興起，主要包括遺傳算法[4]、人工神經(jīng)網(wǎng)絡(luò)[5]、多項(xiàng)式回歸[6]、灰色理論[7]等方法。許多學(xué)者運(yùn)用這些方法對(duì)能反映近岸海域富營(yíng)養(yǎng)化狀況的葉綠素a進(jìn)行預(yù)測(cè)，并取得了一定成果。其中，BP神經(jīng)網(wǎng)絡(luò)(BP-ANN)在葉綠素a濃度預(yù)測(cè)中運(yùn)用最廣，實(shí)例最多。但是神經(jīng)網(wǎng)絡(luò)的經(jīng)驗(yàn)風(fēng)險(xiǎn)最小化是基于樣本夠多的情況，在處理小樣本和多變量數(shù)據(jù)時(shí)，并不能很好的保證模型的泛化能力，存在過(guò)擬合現(xiàn)象，即樣本數(shù)較少時(shí)，預(yù)測(cè)結(jié)果不能保證。

支持向量機(jī)(SVM)是根據(jù)Vapnik[8]建立的統(tǒng)計(jì)學(xué)理論為基礎(chǔ)，遵循結(jié)構(gòu)風(fēng)險(xiǎn)最小化原理建立的一種處理數(shù)據(jù)方法。該方法通過(guò)引用核函數(shù)，實(shí)現(xiàn)了在小樣本、非線性的條件下提高模型預(yù)測(cè)的準(zhǔn)確性，避免了BP-ANN局部最優(yōu)和收斂時(shí)長(zhǎng)等問(wèn)題，彌補(bǔ)了神經(jīng)網(wǎng)絡(luò)的不足。

長(zhǎng)江口海域生態(tài)環(huán)境復(fù)雜多樣[9]，隨著人類(lèi)活動(dòng)的增加，長(zhǎng)江口鄰近海域富營(yíng)養(yǎng)化嚴(yán)重，赤潮頻發(fā)，漁業(yè)環(huán)境質(zhì)量下降，嚴(yán)重影響了江蘇沿岸經(jīng)濟(jì)的發(fā)展，因此對(duì)該海域進(jìn)行環(huán)境預(yù)測(cè)技術(shù)研究很有必要。葉綠素a是反映海水中藻類(lèi)生物量的重要指標(biāo)，也是表征水體富營(yíng)養(yǎng)化程度的重要指示之一[10]。預(yù)測(cè)葉綠素a濃度可以為長(zhǎng)江口及鄰近海域富營(yíng)養(yǎng)化的監(jiān)測(cè)和管理提供支持。

1 材料與方法

1.1 樣品采集與預(yù)處理

本文以長(zhǎng)江口鄰近海域作為研究對(duì)象，2015年春季(3月)和夏季(7月)在長(zhǎng)江口鄰近海域(29°77′N(xiāo)～32°25′N(xiāo)，122°00′E～124°00′E)，分別設(shè)置了9個(gè)斷面，60個(gè)站位進(jìn)行現(xiàn)場(chǎng)調(diào)查，共采集了172個(gè)表層和次表層的海水樣品，具體站位分布圖見(jiàn)圖1。具體采樣方法：現(xiàn)場(chǎng)用Niskin采水器根據(jù)站位水深進(jìn)行采水，采水后立即用0.7 μm的GF/F膜過(guò)濾，并用馬弗爐燒過(guò)的錫紙包好冷凍保存，帶回實(shí)驗(yàn)室測(cè)定葉綠素a濃度；過(guò)濾后的水樣部分裝在灼燒后的棕色玻璃瓶(已預(yù)先在400 ℃的馬弗爐中灼燒4 h)中冷凍保存，用于測(cè)定有色溶解有機(jī)物(CDOM)的特征吸收系數(shù)；另一部分水樣冷凍保存于100 mL的聚乙烯瓶中，用于總氮、總磷的測(cè)定。

圖1 2015年3和7月采樣站位示意圖Fig.1 Map of sampling station in March and July 2015

1.2 樣品的測(cè)定

溫度(T)、鹽度(S)和溶解氧(DO)：利用多參數(shù)水質(zhì)儀CTD獲得。

TN和TP：采用磷鉬藍(lán)法和Cr-Cu還原法獲得，即先用含硼酸的堿性過(guò)硫酸鉀氧化消解海水樣品，再用Cr-Cu還原法把氧化后的硝酸鹽轉(zhuǎn)化為亞硝酸鹽，用重氮-偶氮反應(yīng)顯色測(cè)定，無(wú)機(jī)磷酸鹽則采用磷鉬藍(lán)法測(cè)定。

CDOM的特征吸收系數(shù)：采用紫外可見(jiàn)分光光度法測(cè)定，測(cè)樣儀器為島津UV-2550紫外可見(jiàn)分光光度計(jì)，并用5 cm石英比色皿，以Mill-Q水為參比測(cè)定過(guò)濾水樣。

葉綠素a：采用分光光度法獲得，先將冷凍的GF/F濾膜放入離心管，并用10 mL 90%的丙酮溶液振蕩萃取得到上清液，將上清液置于1 cm比色皿中，用島津UV-2550紫外可見(jiàn)分光光度計(jì)，以丙酮作參比,測(cè)定630、647、664、750 nm的其吸光度，再利用Jeffrey-Humphrey的方程計(jì)算Chl-a的濃度[11]：

Chl-a(mg·m-3)=[11.85×(A664-A750)-

1.54×(A647-A750)-0.08×(A630-A750)]×

Ve/(L×Vf)。

其中：L為比色皿長(zhǎng)度；Ve為萃取液的體積；Vf為過(guò)濾水樣的體積；A為吸光度。

1.3 研究方法

支持向量機(jī)(SVM)的基本思想就是利用核函數(shù)將低維輸入空間中線性不可分的點(diǎn)映射成高維特征空間中線性可分的點(diǎn)，并通過(guò)劃分超平面使同類(lèi)樣本之間相似性盡可能的大，即所有的點(diǎn)到分類(lèi)超平面的距離最大化，達(dá)到最大泛化能力[12-13]。

近海海域富營(yíng)養(yǎng)化是一個(gè)多因素耦合、多維度協(xié)同作用的結(jié)果，關(guān)系復(fù)雜且具有多維的非線性特征[14-15]，此外，葉綠素a濃度與各影響因素之間也存在典型的非線性關(guān)系，而支持向量機(jī)回歸(SVR)就是將實(shí)際問(wèn)題通過(guò)非線性映射到高維特征空間，并在高維特征空間構(gòu)建線性回歸，從而得到低維空間的非線性回歸效果[16]。支持向量機(jī)回歸(SVR)模型的樣本只有一類(lèi)，所尋求的最優(yōu)平面是使所有樣本點(diǎn)離超平面的“總偏差”最小，樣本點(diǎn)都在兩條邊界之間，相當(dāng)于求最大間隔的問(wèn)題[17]。

支持向量機(jī)回歸模型與人工神經(jīng)網(wǎng)絡(luò)類(lèi)似，主要包括三個(gè)部分：輸入層、中間層和輸出層，具體模型見(jiàn)圖2。

圖2 支持向量機(jī)模型的輸入與輸出Fig.2 The inputs and outputs of support vector machine

其中：x1～xn為低維輸入向量，即xi=(x1,x2, …,xn)；K(xi,x) =〈φ(xi) ·φ(x)〉為核函數(shù)，xi為支持向量，φ向高維空間映射的函數(shù)[13]，通過(guò)對(duì)應(yīng)支持向量的拉格朗日系數(shù)(β1, β2, …, βn)連接輸入向量和核函數(shù)就能得到線性組合函數(shù)f(x)。

(1)

兩個(gè)低維空間中的輸入向量經(jīng)過(guò)某種變化后計(jì)算出其在高維空間中的向量?jī)?nèi)積值即為核函數(shù)，巧妙的避免了向量由低維向高維空間映射時(shí)計(jì)算復(fù)雜的問(wèn)題。目前支持向量機(jī)回歸模型常用的核函數(shù)有線性核函數(shù)、多項(xiàng)式核函數(shù)、Sigmoid核函數(shù)和徑向基核函數(shù)(RBF核函數(shù))[18]。與多項(xiàng)式和Sigmoid核函數(shù)相比，徑向基核函數(shù)參數(shù)少，更適合非線性映射，且具有較寬的收斂域，計(jì)算更為簡(jiǎn)單，因此，徑向基核函數(shù)應(yīng)用更廣[18-19]。本文使用徑向基核函數(shù)，公式如下：

(2)

其中g(shù)為核參數(shù)。

核函數(shù)選定后，需要確定相應(yīng)的最優(yōu)懲罰參數(shù)c和核參數(shù)g，其中c主要表示懲罰系數(shù)，即對(duì)誤差的容忍度，g表示主要影響樣本數(shù)據(jù)子空間分布的復(fù)雜程度。最優(yōu)懲罰參數(shù)c和核參數(shù)g通常通過(guò)K-折交叉驗(yàn)證法(K-fold Cross Validation)得到[20]。具體操作是將數(shù)據(jù)集平均分為K組，輪流將其中的K-1組做訓(xùn)練，剩余的1組做驗(yàn)證，在給定參數(shù)情況下，K次的結(jié)果的均值作為對(duì)模型的評(píng)價(jià)指標(biāo)，對(duì)支持向量機(jī)回歸模型表示為預(yù)測(cè)值和實(shí)際值的均方誤差(MSE)。

2 模型的建立

2.1 參數(shù)選擇及預(yù)處理

有色溶解有機(jī)物(CDOM)是水體中溶解有機(jī)物(DOM)的基本組分，能夠影響控制C、N、P等元素的生物地球化學(xué)循環(huán)和浮游植物進(jìn)行光合作用[21]。CDOM的特征吸收系數(shù)與水體DOM含量及性質(zhì)、水體濁度等密切相關(guān)，是近海海域生態(tài)環(huán)境監(jiān)測(cè)主要指標(biāo)的組成部分[22]。水溫(T)、鹽度(S)、總氮(TN)、總磷(TP)、溶解氧(DO)這些參數(shù)是評(píng)價(jià)海水水質(zhì)的基本指標(biāo)，對(duì)海洋中藻類(lèi)產(chǎn)生直接或間接影響，進(jìn)而影響水體中葉綠素a濃度[23]。因此，本文確定的監(jiān)測(cè)參數(shù)包括水溫(T)、鹽度(S)、總氮(TN)、總磷(TP)、溶解氧(DO)等基本參數(shù)和有色溶解有機(jī)物(CDOM)特征吸收系數(shù)aCDOM(355)和aCDOM(455)。

對(duì)獲得的172個(gè)樣品的溶解氧、鹽度、溫度、TN、TP、aCDOM(355)、aCDOM(455)參數(shù)及測(cè)得的Chl-a濃度進(jìn)行統(tǒng)計(jì)分析，由表1可知，Chl-a的平均值為1.754 9 μg·L-1，變化范圍在0.013 1～18.954 4 μg·L-1之間。其中，長(zhǎng)江口附近海域表層和近岸站位所采集的水樣中Chl-a濃度較高。長(zhǎng)江口近岸海域受陸源輸入影響較大，營(yíng)養(yǎng)物質(zhì)濃度較高，浮游植物生長(zhǎng)旺盛[24]。

表1 各參數(shù)數(shù)值特征Table 1 The numerical characteristic of each parameter

此外，考慮到數(shù)據(jù)中各變量存在量綱和數(shù)量級(jí)的差異，對(duì)數(shù)據(jù)進(jìn)行了歸一化處理以減小數(shù)值差異帶來(lái)的影響[25]。具體處理方法如下：

其中：x′為變量x歸一化后的值；xmax和xmin分別是數(shù)據(jù)的最大值和最小值。

2.2 相關(guān)性分析

對(duì)支持向量機(jī)回歸模型的輸入變量進(jìn)行相關(guān)性分析，去除不相關(guān)或重復(fù)變量，對(duì)于保證模型的合理性以及提高模型的準(zhǔn)確度具有重要意義。本文以葉綠素a濃度作為輸出變量，以水溫、鹽度、總氮(TN)、總磷(TP)、溶解氧和有色溶解有機(jī)物(CDOM)特征吸收系數(shù)aCDOM(355)、aCDOM(455)作為候選輸入變量，分別計(jì)算了各輸入變量與葉綠素a濃度的Pearson相關(guān)系數(shù)，結(jié)果如表2。

由表2可知，在0.01的顯著性水平下，溶解氧(DO)、溫度和吸收系數(shù)aCDOM(355)、aCDOM(455)與Chl-a都具有顯著相關(guān)關(guān)系，相關(guān)系數(shù)在-0.204～0.479之間，而鹽度和TN與Chl-a在0.05的顯著性水平下也呈顯著相關(guān)關(guān)系，其相關(guān)系性系數(shù)分別為-0.192和0.165，這表明本研究所選取的大多參數(shù)與Chl-a之間具有相關(guān)性，考慮到TP與其余參數(shù)在0.01的顯著性水平下有顯著相關(guān)關(guān)系，且TP是海水富營(yíng)養(yǎng)化監(jiān)測(cè)的基本要素，故最終的輸入變量確定為溫度、鹽度、總氮(TN)、總磷(TP)、溶解氧(DO)以及CDOM的紫外特征吸收系數(shù)aCDOM(355)和aCDOM(455)。

表2 輸入變量和葉綠素a濃度的Pearson相關(guān)系數(shù)Table 2 The Pearson correlation coefficient between the input variables and Chl-a

注:**表示P<0.01；*表示P< 0.05。**indicates very significant association；*indicates very significant association.

2.3 支持向量機(jī)回歸模型的建立

從172個(gè)樣品中隨機(jī)抽取112個(gè)樣品作為訓(xùn)練集，剩余的60個(gè)樣品作驗(yàn)證集。以CDOM特征吸收系數(shù)aCDOM(355)和aCDOM(455)以及溫度、鹽度、溶解氧等7個(gè)參數(shù)作為輸入變量，以葉綠素a濃度為因變量。

支持向量機(jī)回歸采用臺(tái)灣大學(xué)林智仁開(kāi)發(fā)設(shè)計(jì)的LIBSVM-3.1工具包實(shí)現(xiàn)，在MATLAB平臺(tái)下進(jìn)行建模[26]，設(shè)置核函數(shù)為徑向基核函數(shù)，相應(yīng)參數(shù)為最優(yōu)懲罰參數(shù)c和核參數(shù)g。為確定c和g最佳參數(shù)值，將c和g分別取以2為底的指數(shù)離散值，代入 K-CV交叉驗(yàn)證的算法中，選取這K個(gè)模型中平均驗(yàn)證準(zhǔn)確度最大，即平均驗(yàn)證均方根誤差(MSE)最小的 那組c、g值作為該模型的參數(shù)，該方法被稱為“網(wǎng)格尋優(yōu)法”(GS)。為減少計(jì)算量，把c、g的間隔設(shè)置大一點(diǎn)，再通過(guò)最佳參數(shù)位置范圍逐漸減小其范圍和間隔，進(jìn)行精細(xì)的網(wǎng)格尋優(yōu)，以此確定最終的參數(shù)值。

設(shè)置模型參數(shù)c∈ {2-10,2-9.5, …,210}，g∈ {2-10,2-9.5, …,210}，V=10，進(jìn)行網(wǎng)格尋優(yōu)搜索，結(jié)果如圖3所示。

圖3 網(wǎng)格尋優(yōu)搜索結(jié)果Fig.3 The optimization results of Grid Search

利用網(wǎng)格尋優(yōu)搜索方法得到的最佳參數(shù)值為：c=11.313 4、g=0.5,得到最佳參數(shù)后，核函數(shù)為，按該參數(shù)進(jìn)行設(shè)置后，輸入訓(xùn)練集數(shù)據(jù)，就能得到最終的葉綠素a濃度預(yù)測(cè)模型，該支持向量機(jī)回歸模型的函數(shù)可表達(dá)為：

模型的性能由可決系數(shù)(R2)和均方誤差(MSE)決定，R2表示測(cè)量值與預(yù)測(cè)值之間的相關(guān)性，R2越接近于1，表示樣本的預(yù)測(cè)值對(duì)實(shí)測(cè)樣本的擬合度越好，模型的擬合效果越好[27]。而MSE主要用來(lái)表征樣本數(shù)據(jù)之間的變化程度，MSE的數(shù)值越小，表明預(yù)測(cè)模型對(duì)實(shí)驗(yàn)數(shù)據(jù)的分析具有越好的精確度[28]，該模型的MSE=0.048 7。

2.4 訓(xùn)練和驗(yàn)證結(jié)果

將訓(xùn)練集和驗(yàn)證集數(shù)據(jù)輸入上述預(yù)測(cè)模型中，對(duì)輸出值進(jìn)行反歸一化，得到模擬的葉綠素a濃度，并將實(shí)測(cè)葉綠素a濃度值和預(yù)測(cè)值進(jìn)行對(duì)比，對(duì)比結(jié)果具體見(jiàn)圖4所示。

由圖4可知，以7個(gè)變量構(gòu)建的GS-SVR模型所輸出的葉綠素a濃度預(yù)測(cè)值和實(shí)測(cè)值在變化趨勢(shì)上大致相同。在訓(xùn)練集中，模擬值和實(shí)測(cè)值在0.01的顯著性水平下，Pearson相關(guān)系數(shù)為0.886(p<0.01)，均方誤差MSE為0.024 0；而驗(yàn)證集中，模擬值和實(shí)測(cè)值在0.01的顯著性水平下，Pearson相關(guān)系數(shù)為0.840(p<0.01)，均方誤差MSE為0.041 8。Zhang等[29]利用基于主成分分析(PCA)方法的模糊BP神經(jīng)網(wǎng)絡(luò)模型預(yù)測(cè)中國(guó)東海近岸海域葉綠素a濃度，預(yù)測(cè)結(jié)果與實(shí)測(cè)值具有良好一致性，MSE為0.109；Rocha等[30]通過(guò)多元線性回歸方法預(yù)測(cè)了巴西帕爾杜河葉綠素a濃度，其預(yù)測(cè)結(jié)果與實(shí)測(cè)值的Pearson相關(guān)系數(shù)為0.520；Zheng等[31]利用元胞自動(dòng)機(jī)與支持向量機(jī)結(jié)合(CA-SVM)建立了渤海灣葉綠素a濃度預(yù)測(cè)模型，其預(yù)測(cè)結(jié)果與實(shí)測(cè)值的R2為0.861，均方差MSE為0.190。相較而言，本研究所建立的模型得到的預(yù)測(cè)值和實(shí)測(cè)值具有更好的一致性。

圖4 實(shí)際值與預(yù)測(cè)值對(duì)比Fig.4 Comparison of measured values and predicted values

2.5 輸入變量相對(duì)重要性評(píng)價(jià)

根據(jù)7個(gè)輸入變量建立的支持向量機(jī)回歸預(yù)測(cè)模型，由w=∑s.v.βiφ(xi)求出各變量權(quán)重系數(shù)，并得到各輸入變量對(duì)輸出變量的重要性，將權(quán)重系數(shù)最大的溫度賦值100，可以得到各輸入?yún)?shù)的相對(duì)重要性，具體情況見(jiàn)表3。

表3 輸入變量的權(quán)重系數(shù)和相對(duì)重要性Table 3 Input variables in importance according to their weights and standardized weights

由上述表3中的權(quán)重系數(shù)可知，支持向量機(jī)回歸預(yù)測(cè)模型中對(duì)葉綠素a濃度預(yù)測(cè)影響顯著的輸入變量是溫度和CDOM 特征吸收系數(shù)aCDOM(355)。海水中的浮游植物通過(guò)光合作用進(jìn)行初級(jí)生產(chǎn)，而葉綠素a是浮游植物進(jìn)行光合作用的重要色素。有研究表明，溫度是一切酶促反應(yīng)的控制因子，水溫與浮游植物的初級(jí)生產(chǎn)密切相關(guān)[32]，浮游植物代謝率和光合作用暗反應(yīng)都取決于水溫，當(dāng)光照充足時(shí),光合作用的速度與溫度呈正相關(guān)[33]。

除此以外，鹽度和溶氧對(duì)長(zhǎng)江口鄰近海域葉綠素a濃度預(yù)測(cè)模型影響也較大。鹽度是反應(yīng)近岸海域特別是海域陸源輸入的常用指標(biāo)[41]，在這些區(qū)域鹽度與營(yíng)養(yǎng)鹽之間都有較為明顯的負(fù)相關(guān)關(guān)系，水系混合影響著營(yíng)養(yǎng)鹽的消長(zhǎng)[42]。溶解氧是海洋浮游植物光合作用的產(chǎn)物，也是海洋中影響異養(yǎng)生物活動(dòng)的主要因素[43]，是衡量海水水質(zhì)的基本參數(shù)之一。

3 結(jié)語(yǔ)

基于CDOM 特征吸收系數(shù)aCDOM(355)和aCDOM(455)以及溫度、鹽度、溶解氧、TP、TN等5個(gè)基本水質(zhì)參數(shù)作為輸入變量，利用支持向量機(jī)回歸(SVR)建立了長(zhǎng)江口鄰近海域葉綠素a濃度預(yù)測(cè)模型，預(yù)測(cè)值與實(shí)測(cè)值具有較好的一致性，且溫度和CDOM 特征吸收系數(shù)影響顯著，表明該模型能較好的預(yù)測(cè)長(zhǎng)江口鄰近海域葉綠素a濃度，可為長(zhǎng)江口及鄰近海域富營(yíng)養(yǎng)化監(jiān)測(cè)提供技術(shù)支持。