華祖林，韓愛秋

(1．河海大學(xué)淺水湖泊綜合治理與資源開發(fā)教育部重點(diǎn)實(shí)驗(yàn)室，江蘇南京 210098;2．河海大學(xué)水資源高效利用與工程安全國家工程研究中心，江蘇南京 210098;3．河海大學(xué)環(huán)境學(xué)院，江蘇南京 210098)

湖泊參照狀態(tài)指湖泊受到人類活動(dòng)影響最小或環(huán)境系統(tǒng)可達(dá)到的最佳狀態(tài)[1]，它是制定水質(zhì)標(biāo)準(zhǔn)的科學(xué)依據(jù)，也是評(píng)估、預(yù)防、控制和管理湖泊富營養(yǎng)化的基礎(chǔ)，更是水污染治理和生態(tài)修復(fù)期望達(dá)到的“終極”狀態(tài)。壓力-響應(yīng)模型是最適用于確定淺水湖泊參照狀態(tài)的模型之一。Lamon等[2]探索性地以貝葉斯多層次線性回歸模型確立了美國不同類型湖泊壓力變量和響應(yīng)變量的響應(yīng)關(guān)系;US EPA[3]介紹了利用線性回歸、多元線性回歸和非參數(shù)拐點(diǎn)分析等建立壓力-響應(yīng)模型的方法;Haggard等[4]利用線性回歸和分類回歸樹相結(jié)合的方法確定了美國紅河流域的營養(yǎng)物基準(zhǔn)值;Qian等[5]構(gòu)建了連續(xù)變量貝葉斯網(wǎng)絡(luò)建?？蚣埽摽蚣芙Y(jié)合了貝葉斯網(wǎng)絡(luò)模型和傳統(tǒng)經(jīng)驗(yàn)統(tǒng)計(jì)模型，并確定了美國俄亥俄州溪流的營養(yǎng)物基準(zhǔn)值;在中東部平原湖泊生態(tài)區(qū)、云貴高原湖泊生態(tài)區(qū)等8個(gè)湖泊營養(yǎng)物一級(jí)生態(tài)區(qū)的基礎(chǔ)上[6-7]，Huo等[8-11]引入壓力-響應(yīng)模型來制訂中國受人類活動(dòng)影響較嚴(yán)重湖泊的參照狀態(tài)，利用線性回歸壓力-響應(yīng)模型制訂了東部平原湖泊生態(tài)區(qū)和云貴湖泊生態(tài)區(qū)的營養(yǎng)物基準(zhǔn)，并在考慮季節(jié)影響水質(zhì)的情況下采用分類回歸樹與拐點(diǎn)分析相結(jié)合的方法建立壓力-響應(yīng)模型深入研究湖泊生態(tài)分區(qū)參照狀態(tài)的確定;吳超等[12]采用分類回歸樹和線性回歸分別建立了壓力-響應(yīng)模型，并通過相互驗(yàn)證確定了太湖TN、TP等的參照狀態(tài);張亞麗等[13]利用貝葉斯拐點(diǎn)分析等3種非線性方法確定了中國9個(gè)典型湖泊的參照狀態(tài);霍守亮等[14-15]利用分類回歸樹、拐點(diǎn)分析法、貝葉斯層次線性回歸建立了壓力-響應(yīng)模型，針對(duì)不同湖泊生態(tài)分區(qū)進(jìn)行參照狀態(tài)的確定和案例研究，并進(jìn)行了歸納總結(jié)。

以上研究成果有值得商榷和改進(jìn)之處，如線性回歸模型要求響應(yīng)關(guān)系是線性的、因變量的誤差呈正態(tài)分布、采集的數(shù)據(jù)樣本獨(dú)立等，而實(shí)際上，光照、水深等環(huán)境因素會(huì)影響變量間的響應(yīng)關(guān)系，導(dǎo)致其難以呈現(xiàn)出線性;此外，湖泊實(shí)際監(jiān)測數(shù)據(jù)來源復(fù)雜，且為時(shí)間序列，很難滿足線性回歸分析設(shè)定的正態(tài)性和獨(dú)立性等假設(shè)理論。貝葉斯層次回歸模型雖然能有效緩解數(shù)據(jù)測量誤差的問題，但需要在調(diào)查時(shí)以明顯分層的方式采集原始數(shù)據(jù)，以便分層建立模型。分類回歸樹壓力-響應(yīng)模型適用于識(shí)別對(duì)模型結(jié)果變化有顯著貢獻(xiàn)的壓力變量，但在數(shù)據(jù)量較少時(shí)缺乏穩(wěn)健性，壓力變量很小的變化可能引起模型結(jié)果較大的變化，使結(jié)果的準(zhǔn)確性得不到保障;拐點(diǎn)分析法需要進(jìn)一步研究低于營養(yǎng)物拐點(diǎn)對(duì)應(yīng)的響應(yīng)變量數(shù)值是否能夠支持水體的指定用途[14-15]。這些限制性條件會(huì)導(dǎo)致推斷的參照狀態(tài)濃度過高或過低。壓力-響應(yīng)模型方法的改進(jìn)一直是學(xué)者們孜孜不倦探索的課題。

本文針對(duì)原有壓力-響應(yīng)模型需要滿足設(shè)定條件等特點(diǎn)，考慮到響應(yīng)關(guān)系應(yīng)是非線性的，以及異常點(diǎn)會(huì)降低模型穩(wěn)健性等情況，嘗試性地采用非參數(shù)LOWESS穩(wěn)健回歸來改進(jìn)淺水湖泊的壓力-響應(yīng)模型。改進(jìn)模型能消除異常點(diǎn)對(duì)估計(jì)結(jié)果的影響，對(duì)于非線性、非齊次問題有較好的效果。太湖作為東部平原湖泊生態(tài)區(qū)淺水湖泊的典型代表，參照狀態(tài)的確定尤為重要[16-17]，以太湖為研究對(duì)象進(jìn)行湖泊參照狀態(tài)的確定計(jì)算。

1 非參數(shù)LOWESS穩(wěn)健回歸模型的構(gòu)建

非參數(shù)LOWESS穩(wěn)健回歸是擬合難以用具體函數(shù)描述響應(yīng)關(guān)系數(shù)據(jù)的有效方法，其基本思想是先用局部多項(xiàng)式進(jìn)行擬合，然后定義穩(wěn)健的權(quán)數(shù)并進(jìn)行光滑處理，重復(fù)運(yùn)算數(shù)次后就可以得到自變量與因變量穩(wěn)健的擬合曲線[18]。具體步驟如下：

步驟1：將壓力變量TP和響應(yīng)變量Chl-a關(guān)系曲線表示為m(x)，利用局部多項(xiàng)式進(jìn)行擬合，得到擬合值 ^m(xi)：

式中：xi為TP的觀測值，mg/L;β^(xi)為對(duì)局部參數(shù)β(xi)的最小二乘加權(quán)回歸。

步驟2：定義并計(jì)算穩(wěn)健權(quán)數(shù)δk：

其中 ek=m(xk)-m^(xk)

式中：B為二次核函數(shù)，定義為B(x)=(1-x2)2I(I為符號(hào)函數(shù)，當(dāng)|x|＜1時(shí)I=1;當(dāng)|x|≥1時(shí)I=0);ek為殘差;xk為加權(quán)的xi;s為殘差絕對(duì)值的中位數(shù)，s=median(e1， e2，…， en)。

步驟3：利用δkwk(xi)對(duì)m(x)進(jìn)行局部加權(quán)最小二乘估計(jì)，得到新的模型殘差ek。其中

式中：W為距離權(quán)函數(shù);hi為光滑參數(shù)，hi= xi-xj，xj(j=1，2，…，n)為 xi相鄰的 TP 觀測值，hi為計(jì)算得出的r個(gè)數(shù)中最小的一個(gè)。

步驟4：重復(fù)第2步和第3步，最后第N次迭代得到的擬合值產(chǎn)生TP和Chl-a的響應(yīng)關(guān)系曲線。

穩(wěn)健權(quán)數(shù)可將異常值排除在外，并且初始?xì)埐畲?小)的觀測值在下一次加權(quán)最小二乘法中的權(quán)數(shù)就小(大)，重復(fù)幾次后就可把異常值不斷地排除在外，最終得到穩(wěn)健的響應(yīng)關(guān)系曲線。Cleveland[18]推薦迭代次數(shù)N=3。

非參數(shù)LOWESS穩(wěn)健回歸壓力-響應(yīng)模型構(gòu)建流程見圖1。

圖1 壓力-響應(yīng)模型構(gòu)建流程

2 模型數(shù)據(jù)與壓力變量的確定

2.1 模型數(shù)據(jù)

利用LOWESS穩(wěn)健回歸方法建立壓力-響應(yīng)模型推斷太湖參照狀態(tài)，數(shù)據(jù)來源于文獻(xiàn)[19]。文獻(xiàn)[19]記錄了圖2太湖湖泊生態(tài)系統(tǒng)國家野外科學(xué)觀測研究站(簡稱太湖站)全太湖32個(gè)野外站點(diǎn)中8個(gè)站點(diǎn)自1991—2006年的逐月觀測數(shù)據(jù)?？紤]到參照狀態(tài)是受人類活動(dòng)影響最小或環(huán)境系統(tǒng)可達(dá)到的最佳狀態(tài)，所以本文采用位于湖心區(qū)域受人類活動(dòng)影響相對(duì)較小的7號(hào)和8號(hào)站點(diǎn)的數(shù)據(jù)。華祖林等[20]采用季節(jié)分解的非參數(shù)法發(fā)現(xiàn)1999—2001年太湖湖心區(qū)域的TP與Chl-a的觀測值比1995—2006年的觀測值更適合用于推斷太湖的參照狀態(tài)。本文使用7、8號(hào)站點(diǎn)1999—2001年共3年108個(gè)觀測值進(jìn)行分析。

圖2 太湖站監(jiān)測點(diǎn)圖

2.2 壓力變量確定和數(shù)據(jù)處理

首先對(duì)模型數(shù)據(jù)進(jìn)行分析。1999—2001年，太湖湖心TN、TP以及Chl-a質(zhì)量濃度平均值、中位數(shù)、均方差等見表1。1999—2001年太湖湖心觀測值穩(wěn)定性按從大到小順序?yàn)門P、TN、Chl-a。在峰度和偏度檢驗(yàn)中，TN、TP和Chl-a均存在正偏態(tài)現(xiàn)象，偏重程度按從大到小順序?yàn)镃hl-a、TN、TP，且三者均偏離正態(tài)分布。綜合考慮數(shù)據(jù)穩(wěn)定性及峰度和偏度檢驗(yàn)結(jié)果，在選擇壓力變量時(shí)應(yīng)優(yōu)先考慮TP。

在研究湖泊氮磷與Chl-a關(guān)系時(shí)，TN與TP質(zhì)量濃度比值(ρ(TN)/ρ(TP))是一個(gè)關(guān)鍵參數(shù)，常根據(jù)ρ(TN)/ρ(TP)判別藻類受營養(yǎng)鹽限制的類型。ρ(TN)/ρ(TP)比較大(＞17)時(shí)藻類受磷限制，ρ(TN)/ρ(TP)比較小( ＜10)時(shí)藻類受氮限制，ρ(TN)/ρ(TP)中等(10～17)時(shí)藻類受二者共同限制。表1顯示太湖1999—2001年TN質(zhì)量濃度平均值為1．69 mg/L，TP 質(zhì)量濃度平均值為 0．08 mg/L，ρ(TN)/ρ(TP)=21．12，大于 17。 可見，1999—2001年太湖為磷限制型湖泊。張波等[21]研究表明，太湖水體存在一定的固氮作用，固氮速率為每小時(shí)1．53 ng/L，導(dǎo)致水體中TN與Chl-a響應(yīng)關(guān)系不明顯。因此分析太湖TP和Chl-a的響應(yīng)關(guān)系更符合壓力-響應(yīng)模型的釋義。

圖3 TP與Chl-a冪變換數(shù)據(jù)Q-Q圖

表1 1999—2001年TN、TP與Chl-a統(tǒng)計(jì)特征比較

由于原始逐月觀測數(shù)據(jù)TP與Chl-a存在數(shù)量級(jí)差別，因此需要對(duì)原始數(shù)據(jù)進(jìn)行處理。本文采用Qian[22]提出的冪變換對(duì)TP和Chl-a數(shù)據(jù)進(jìn)行處理。圖3是TP數(shù)據(jù)經(jīng)過對(duì)數(shù)變換后和Chl-a 0．1次方變化后的Q Q圖。從圖3可以明顯地看出，TP數(shù)據(jù)點(diǎn)與標(biāo)準(zhǔn)正態(tài)分布直線擬合效果較差，存在嚴(yán)重的偏離。結(jié)合表1峰度和偏度檢驗(yàn)結(jié)果，說明1999—2001年太湖湖心TP觀測數(shù)據(jù)和冪變換后均不滿足正態(tài)性假設(shè)，所以基于上述假設(shè)建立的壓力-響應(yīng)模型具有不合理性，進(jìn)一步說明建立非參數(shù)方法壓力-響應(yīng)模型的必要性。

在采用壓力-響應(yīng)模型推斷參照狀態(tài)時(shí)，需要給定一個(gè)響應(yīng)變量基準(zhǔn)值。在保證水體使用功能的前提下推斷營養(yǎng)物基準(zhǔn)時(shí)，Chl-a是聯(lián)系營養(yǎng)物濃度的重要響應(yīng)變量[8]。參考 USEPA[23]和霍守亮等[14，24-25]對(duì)Chl-a基準(zhǔn)值的研究，本文在建立壓力-響應(yīng)模型時(shí)設(shè)定4．73 μg/L Chl-a為太湖響應(yīng)變量基準(zhǔn)值。

3 模型驗(yàn)證

采用LOWESS穩(wěn)健回歸建立1999—2001年太湖湖心區(qū)域TP與Chl-a的壓力-響應(yīng)模型，平滑參數(shù)f=0．2，迭代次數(shù) N=3。圖4為太湖湖心區(qū)域1999—2001年Chl-a與TP壓力-響應(yīng)模型及模型殘差。從圖4的散點(diǎn)可以看出，Chl-a與TP并不呈現(xiàn)出簡單的線性響應(yīng)關(guān)系，而是復(fù)雜的非線性響應(yīng)關(guān)系。根據(jù)確定的響應(yīng)變量Chl-a基準(zhǔn)值4．73 μg/L，可以得出壓力變量TP的參照狀態(tài)質(zhì)量濃度為0．018 mg/L。此外，結(jié)合置信區(qū)間和模型殘差評(píng)估模型結(jié)果的準(zhǔn)確性，利用自助法求得1999—2001年太湖TP參照狀態(tài)質(zhì)量濃度的95%置信區(qū)間為0．013～0．030 mg/L。 從圖4(b)可以看出，模型殘差均值為0．0mg/L左右，下25%分位點(diǎn)為-0．06 mg/L，上25%分位點(diǎn)為0．04mg/L，說明LOWESS穩(wěn)健回歸方法可用于建立壓力-響應(yīng)模型進(jìn)行參照狀態(tài)的推斷。

圖4 TP與Chl-a壓力-響應(yīng)模型及模型殘差

表2顯示的是不同地區(qū)TP參照狀態(tài)質(zhì)量濃度。鄭丙輝等[26]應(yīng)用頻率分析法建立了太湖TP參照狀態(tài)，并結(jié)合時(shí)間參照狀態(tài)法和沉積物反演法進(jìn)行驗(yàn)證，綜合分析得出TP參照狀態(tài)為0．03 mg/L;華祖林等[20]利用1999—2001年TP數(shù)據(jù)5%分位點(diǎn)的值作為太湖參照狀態(tài)質(zhì)量濃度，其結(jié)果為0．03 mg/L，并利用幾何分布分塊自助法確定其95%置信區(qū)間為0．025～0．046 mg/L。采用頻率分析法在確立湖泊參照狀態(tài)時(shí)需要確定頻率分布分位點(diǎn)，通常都是選擇上(下)25%或5%，由于25%或5%是人為確定的，因而結(jié)果受到人為影響大。吳超等[12]利用分類回歸樹和線性回歸法建立的壓力-響應(yīng)模型相互驗(yàn)證，確定太湖流域TP質(zhì)量濃度基準(zhǔn)值為0．67mg/L。

表2 不同地區(qū)TP參照狀態(tài)濃度

霍守亮等[14]利用非參數(shù)拐點(diǎn)分析和貝葉斯拐點(diǎn)分析，得到不同湖泊生態(tài)區(qū)壓力變量TP突變點(diǎn)的質(zhì)量濃度范圍為0．015～0．222mg/L，利用簡單線性回歸模型得到中東部湖泊生態(tài)分區(qū)TP參照狀態(tài)質(zhì)量濃度范圍為(0．022±0．007)mg/L;顧莉等[27]將1960年總堿度實(shí)測值代入其建立的改進(jìn)MEI模型得到太湖TP質(zhì)量濃度為0．018mg/L。在20世紀(jì)60年代對(duì)太湖的調(diào)查研究表明，太湖TP質(zhì)量濃度為0．01～0．05 mg/L[28]。 太湖藍(lán)藻暴發(fā)最近30年才出現(xiàn)，60年代太湖的狀態(tài)應(yīng)該比較接近其歷史上受人類活動(dòng)影響較小時(shí)的情況。本文采用季節(jié)分解的非參數(shù)局部線性回歸模型推斷得出更適合建立太湖參照狀態(tài)的數(shù)據(jù)，而且非參數(shù)LOWESS穩(wěn)健回歸不需要設(shè)定分位點(diǎn)，計(jì)算過程受人為影響小，計(jì)算得到的太湖TP參照狀態(tài)質(zhì)量濃度為0．018 mg/L，在線性回歸和拐點(diǎn)方法得出的參照狀態(tài)質(zhì)量濃度范圍內(nèi)。因此，本文得到的結(jié)果可以作為太湖參照狀態(tài)質(zhì)量濃度，并且是可信的。

4 結(jié) 語

以季節(jié)分解的非參數(shù)方法驗(yàn)證的適合推斷太湖參照狀態(tài)的數(shù)據(jù)為基礎(chǔ)，結(jié)合峰度和偏度檢驗(yàn)，并考慮湖泊營養(yǎng)鹽的限制問題，以TP作為壓力變量，構(gòu)建了穩(wěn)健性高、對(duì)非線性問題適應(yīng)性強(qiáng)的LOWESS穩(wěn)健回歸的湖泊壓力-響應(yīng)模型。根據(jù)設(shè)定的Chl-a推斷TP的參照狀態(tài)質(zhì)量濃度為0．018 mg/L;通過自助法得出其95%置信區(qū)間為0．013～0．030 mg/L，同時(shí)從多角度對(duì)所建模型進(jìn)行了驗(yàn)證，表明結(jié)果較為合理可信。