魏 峣，陳卓爾，楊明方，佟洪金，劉曉聰

(1.四川省生態(tài)環(huán)境科學(xué)研究院，成都 610041；2.四川省環(huán)保科技工程有限責(zé)任公司，成都 610041)

前 言

水體溶解氧(Dissovled Oxygen，DO)是水生生物生理代謝和魚(yú)類生長(zhǎng)繁育的限制因子之一，是衡量水體自凈能力的一個(gè)重要指標(biāo)，也是水體水質(zhì)的重要指標(biāo)之一[1～3]。常溫下水中DO含量一般為8～14mg/L，DO低于5 mg/L時(shí)便對(duì)水生生物的生存產(chǎn)生不利影響，低于2 mg/L時(shí)大多數(shù)魚(yú)類可能死亡[4]。水體DO的變化受到光合作用、呼吸作用、大氣復(fù)氧過(guò)程和沉積物耗氧等過(guò)程的綜合影響[5～7]。其中，光合作用是水體溶解氧最重要的來(lái)源。歐陽(yáng)瀟然[8]對(duì)太湖研究發(fā)現(xiàn)，藍(lán)藻水華易發(fā)生和堆積的太湖梅梁灣，浮游植物光合作用制氧是梅梁湖水中最重要的來(lái)源；浮游植物的呼吸與死亡耗氧是溶解氧的最大的消耗者。與此相對(duì)應(yīng)，水體中DO的變化受溫度、pH、高錳酸鹽、懸浮物、營(yíng)養(yǎng)物質(zhì)、生物量等多種因素影響[9]。

Logistic回歸分析是研究因變量是分類、定性變量的常用方法，常用于流行病學(xué)危險(xiǎn)因素分析[10-11]，可用于判別引發(fā)事件的因素[12]，預(yù)測(cè)事件發(fā)生的概率[13]。近年來(lái)在自然災(zāi)害[14-15]、環(huán)境[16-18]等領(lǐng)域也有廣泛應(yīng)用。采用Logistic回歸模型分析水質(zhì)、水文、氣象三方面因素對(duì)釜溪河碳研所斷面發(fā)生低DO影響程度，識(shí)別引發(fā)低DO的危險(xiǎn)因素，分析低DO成因，對(duì)開(kāi)展河道水生態(tài)治理修復(fù)具有重要作用。

1 研究區(qū)域

自貢市地處四川盆地南部，位于北緯28°55′37″～29°38′25″、東經(jīng)104°02′57″～105°16′11″之間。地勢(shì)西北高東南低，地形分為低山、丘陵、平壩，其中丘陵分布最廣，占全市總面積的80%以上。屬四川盆地亞熱帶季風(fēng)濕潤(rùn)氣候區(qū)，多年平均氣溫17.9℃，多年平均降雨量1 079.6mm，多年平均相對(duì)濕度80%，多年平均蒸發(fā)量1 010.3mm。全市土地面積4 381km2，常住人口290萬(wàn)人，地區(qū)生產(chǎn)總值1 312億元。

釜溪河是長(zhǎng)江二級(jí)支流、沱江一級(jí)支流。發(fā)源于四川省內(nèi)江市威遠(yuǎn)縣兩母山東北麓，由旭水河和威遠(yuǎn)河在自貢雙河口處匯成干流，流經(jīng)自貢市主城區(qū)，在富順李家灣匯入沱江(圖1)。釜溪河流域面積3 472km2，長(zhǎng)度190km，平均比降2.97‰，多年平均流量42.4m3/s，多年平均徑流量13.4億m3/年。國(guó)控碳研所斷面位于釜溪河自貢市主城區(qū)下游約1.5km。該斷面水質(zhì)長(zhǎng)期較差，被列為《長(zhǎng)江流域劣Ⅴ類國(guó)控?cái)嗝嬲螌ｍ?xiàng)行動(dòng)工作方案》(2019年)中四川省兩個(gè)劣Ⅴ類國(guó)控?cái)嗝嬷弧３钡?、總磷、化學(xué)需氧量等超標(biāo)外，DO濃度是釜溪河流域最低點(diǎn)，春末及夏季容易出現(xiàn)低DO現(xiàn)象，影響水生態(tài)健康。

圖1 研究區(qū)域Fig.1 Study area

2 數(shù)據(jù)來(lái)源與方法

2.1 數(shù)據(jù)來(lái)源

收集了2017～2018年釜溪河自貢水文站日均流量(Discharge，DI)；自貢氣象站日均氣溫(Air Temperature，AT)；碳研所水質(zhì)自動(dòng)站每日6時(shí)點(diǎn)(0、4、8、12、16、20時(shí))DO、總氮(Total Nitrogen，TN)、總磷(Total Phosphorus，TP)小時(shí)數(shù)據(jù)，并計(jì)算得到日均值(部分缺失時(shí)點(diǎn)數(shù)據(jù)采用插值方式補(bǔ)齊)。由于部分天數(shù)水質(zhì)數(shù)據(jù)缺失，研究共收集到444日數(shù)據(jù)(2017年共283日，2018年共161日)。

2.2 研究方法

Logistic回歸模型多采用二分類變量，設(shè)P為某事件(Y=1)發(fā)生的概率，取值范圍為[0，1]，(1-P)為該事件不發(fā)生(Y=0)的概率。將P/(1-P)取自然對(duì)數(shù)lnP/(1-P)，記為logitP，則logitP的取值范圍為(-∞，+∞)。以P為因變量，m個(gè)危險(xiǎn)因素為自變量，建立線性模型(式1)。

(1)

式(1)可變式為概率模型形式(式2)。

(2)

模型中β0為常數(shù)，表示在不接觸任何危險(xiǎn)因素條件下，結(jié)果發(fā)生與不發(fā)生的概率之比的對(duì)數(shù)值(基線值)。系數(shù)βi是某危險(xiǎn)因素Xi存在(Xi=1)與不存在(Xi=0)相比，事件發(fā)生的優(yōu)勢(shì)比(Odds Ratio，OR)的對(duì)數(shù)值。

βi值>0，eβi>1，表示有危險(xiǎn)因素對(duì)事件發(fā)生的優(yōu)勢(shì)大于無(wú)危險(xiǎn)因素；

βi值<0，eβi<1，表示有危險(xiǎn)因素對(duì)事件發(fā)生的優(yōu)勢(shì)低于無(wú)危險(xiǎn)因素；

βi值=0，eβi=1，表示無(wú)論危險(xiǎn)因素是否出現(xiàn)，與事件發(fā)生與否無(wú)關(guān)。

他們披頭散發(fā)，蓬頭垢面，日復(fù)一日，年復(fù)一年，在暗無(wú)天日的井洞里，頑強(qiáng)地堅(jiān)持著。這樣的作業(yè)空間，這樣危險(xiǎn)的工作環(huán)境，每一個(gè)采礦者的生命，隨時(shí)都會(huì)從陽(yáng)間一步踏過(guò)奈何橋，走向死亡。生命之于他們，只是在陽(yáng)世間的短時(shí)間寄存。每一天清晨他們進(jìn)入礦洞，便是生死兩茫茫，晚上，不知道還能不能回到棲息的工棚。

3 結(jié)果與分析

3.1 基本特征

收集的444個(gè)日均值數(shù)據(jù)統(tǒng)計(jì)結(jié)果(表1)，DO平均濃度4.0±1.9 mg/L，最低濃度僅0.5 mg/L，出現(xiàn)在2017年6月。TP平均濃度0.384±0.246 mg/L，最高濃度達(dá)1.612 mg/L，出現(xiàn)在2018年3月。TN平均濃度7.1±2.9 mg/L，最高濃度21.5 mg/L，出現(xiàn)在2018年2月。平均AT為15.6±6.9℃，最高溫31.6℃，出現(xiàn)在2018年8月。平均DI為3.4±14.1 m3/s，最低DI為0.1 m3/s，出現(xiàn)在2017年11月。

表1 各指標(biāo)特征Tab.1 Characteristics of variables

分析各指標(biāo)相關(guān)性(表2)，DO與TP、TN、AT呈顯著負(fù)相關(guān)關(guān)系(P<0.001)，表明較高的TN、TP濃度及AT對(duì)應(yīng)較低的DO。DO與DI(P=0.008)呈顯著正相關(guān)關(guān)系，表明較低的DI對(duì)應(yīng)較低的DO。

表2 各指標(biāo)相關(guān)性分析(n=444)Tab.2 Correlation analysis of variables (n=444)

3.2 危險(xiǎn)事件及危險(xiǎn)因素判別

以地表水Ⅲ類水質(zhì)限值為判別標(biāo)準(zhǔn)，將DO<5mg/L定義為發(fā)生低溶解氧(Y=1)，DO≥5mg/L定義為未發(fā)生低溶解氧(Y=0)。在此判別標(biāo)準(zhǔn)下有65.5%的樣本發(fā)生低溶解氧(表3)。

表3 危險(xiǎn)事件判別標(biāo)準(zhǔn)(n=444)Tab.3 Criteria for risk events

Logistic回歸分析需要將自變量進(jìn)行二分類取值，即選取某一界限值將自變量分為0和1兩組，其中取1表示為危險(xiǎn)因素，取0表示非危險(xiǎn)因素。根據(jù)各指標(biāo)在Y為0和1兩種情形下的累積頻率分布或水質(zhì)標(biāo)準(zhǔn)，劃分自變量危險(xiǎn)因素(表4)。

表4 各指標(biāo)危險(xiǎn)因素判別界限取值Tab.4 Criteria for risk factors of variables

經(jīng)檢驗(yàn)，各指標(biāo)方差膨脹因子(Variance inflation factor，VIF)均小于10(表5)，表明選取的指標(biāo)不存在明顯的多重共線性問(wèn)題[19]。

表5 各指標(biāo)方差膨脹因子Tab.5 VIF of variables

3.3 模型建立

采用方差分析(Analysis of Variance，ANOVA)對(duì)比各指標(biāo)分組間的差異(表6)。TP、TN、AT和DI危險(xiǎn)因素組與非危險(xiǎn)因素組的DO存在顯著差異(P<0.001)。

表6 各自變量方差分析Tab.6 ANOVA of variables

采用Logistic模型進(jìn)行多因素回歸分析，得到各危險(xiǎn)因素的回歸系數(shù)(βi)估計(jì)值，建立回歸模型。模型Hosmer-Lemeshow統(tǒng)計(jì)量P=0.459，大于0.05，表明模型預(yù)測(cè)值與樣本之間不存在顯著差異[20](表7)，即模型擬合較好。各回歸系數(shù)的Wald檢驗(yàn)P值均小于0.001，表明模型回歸系數(shù)顯著(表8)。

表7 Hosmer-Lemeshow 檢驗(yàn)結(jié)果(n=444)Tab.7 Result of Hosmer-Lemeshow test (n=444)

表8 各指標(biāo)回歸系數(shù)估計(jì)值(n=444)Tab.8 Regression coefficients of variables(n=444)

建立的Logistic回歸模型見(jiàn)式(3)。

LogitP=-4.259+2.401TP+1.472TN+3.627AT+1.086DI

(3)

根據(jù)預(yù)測(cè)結(jié)果，模型對(duì)DO<5mg/L的樣本預(yù)測(cè)正確率86.6%，對(duì)DO≥5 mg/L的樣本預(yù)測(cè)正確率73.2%，總體預(yù)測(cè)正確率82%(表9)。

表9 預(yù)測(cè)結(jié)果(n=444)Tab.9 Predicted results(n=444)

3.4 危險(xiǎn)因素識(shí)別

各指標(biāo)OR值均>1，表明TP≥0.2mg/L、TN≥6.4mg/L、AT≥14.2℃、DI≤0.3 m3/s均是發(fā)生低溶解氧的危險(xiǎn)因素(表10)。在控制其他變量后，發(fā)生低溶解氧的概率，AT≥14.2℃是<14.2℃的37.607倍，較高的AT是發(fā)生低溶解氧的首要危險(xiǎn)因素；TP≥0.2mg/L是<0.2 mg/L的11.031倍，TN≥6.4mg/L是<6.4 mg/L的4.358倍，較高的N、P濃度是發(fā)生低溶解氧的重要危險(xiǎn)因素；DI≤0.3 m3/s是>0.3 m3/s的2.962倍，較低的流量是發(fā)生低溶解氧的危險(xiǎn)因素。

表10 危險(xiǎn)因素分析(n=444)Tab.10 Analysis of risk factors(n=444)

3.5 各月DO的主要危險(xiǎn)因素

計(jì)算各指標(biāo)月平均值并按危險(xiǎn)因素分類(圖2)。各危險(xiǎn)因素中，較高AT(≥14.2℃)分布在3～10月，較高TP(≥0.2mg/L)分布于全部12個(gè)月，較高TN(≥6.4mg/L)分布在1～5月和12月，較低DI(≤0.3m3/s)分布在3、5月。時(shí)間尺度上影響最廣泛的是較高TP和較高AT。

圖2 影響各月DO的主要因素Fig.2 Main factors affecting DO of each month

月均DO<5mg/L的月份為2～7月和9～11月，共9個(gè)月。2月主要危險(xiǎn)因素是較高的TP和較高的TN；其余8個(gè)月危險(xiǎn)因素是較高的TP和較高的AT，其中，3、5月還受低流速影響。

4 討 論

AT對(duì)DO的影響主要包括兩個(gè)方面。一是氣體溶解在液體中通常為放熱過(guò)程，水的飽和溶解濃度隨水溫升高而降低[21]。二是水溫通過(guò)影響藻類營(yíng)養(yǎng)鹽攝入、細(xì)胞中酶的活性從而影響藻類生產(chǎn)力[22]。水溫上升，藻類生長(zhǎng)速度加快[23]。趙巧華[24]等研究認(rèn)為水溫是影響太湖春季藻類生長(zhǎng)的關(guān)鍵因素之一。水體DO變化與藻類的光合作用、呼吸作用有關(guān)[25]。當(dāng)藻類數(shù)量增加到一定數(shù)量級(jí)時(shí)，其數(shù)量多少和生命活動(dòng)的旺盛程度控制了水體中DO的變化[26]。無(wú)光照時(shí)，藻類因內(nèi)源呼吸，吸收氧氣放出二氧化碳；有光照時(shí)，藻類光合作用放出氧氣。隨氣溫升高，藻類白天的凈釋氧量和夜間的凈耗氧量均增加，日平均結(jié)果表現(xiàn)為DO的凈消耗，DO濃度逐漸降低。碳研所斷面DO與AT呈顯著負(fù)相關(guān)關(guān)系(P<0.001)，DO隨AT升高而降低，一方面與AT升高導(dǎo)致飽和溶解氧降低有關(guān)。更重要的是，從4月起日均AT均在14.2℃以上，較高的溫度一直持續(xù)到10月。相應(yīng)月份的水溫多在20℃以上，這有利于藻類繁殖和活性增強(qiáng)，呼吸作用耗氧大于光合作用釋氧，導(dǎo)致DO逐漸降低。

N、P營(yíng)養(yǎng)鹽與水體中生物量存在緊密關(guān)聯(lián)[27]。水體中N、P濃度將通過(guò)藻類繁殖及光合作用、呼吸作用影響水體DO。過(guò)量的N、P輸入將導(dǎo)致水體富營(yíng)養(yǎng)化。富營(yíng)養(yǎng)化水體在適宜的環(huán)境條件下可引起藻類大量繁殖，消耗水體中溶解氧[28]。一般認(rèn)為水體富營(yíng)養(yǎng)化的條件為TN≥0.2mg/L，TP≥0.02mg/L[29]。釜溪河自貢城區(qū)段TN平均濃度7.1±2.9mg/L， TP平均濃度0.384±0.246mg/L，屬于富營(yíng)養(yǎng)化水體。在適宜的環(huán)境條件下可能引發(fā)藻類大量繁殖，導(dǎo)致DO降低。

一定河道斷面下，流量與流速呈正相關(guān)。流速(流量)主要影響水體動(dòng)力學(xué)條件。湖泊或者類湖泊型形態(tài)水體一般流速較小，水力停留時(shí)間長(zhǎng)，給藻類營(yíng)造一種穩(wěn)定的水力環(huán)境。較快的流速，水體更新周期更快，不容易造成單一藻類蓄積和過(guò)量增殖[30]。自貢水文站流量≤0.3m3/s多分布在12月至次年6月的枯平水期，較小的流速有利于藻類的聚集增殖，在一定程度上增加了DO消耗。

5 結(jié) 論

5.1 TP≥0.2mg/L、TN≥6.4mg/L、AT≥14.2℃、DI≤0.3 m3/s是碳研所斷面發(fā)生低溶解氧的危險(xiǎn)因素。其中，較高的AT是發(fā)生低溶解氧的首要危險(xiǎn)因素，較高的N、P濃度是發(fā)生低溶解氧的重要危險(xiǎn)因素，較低的DI是發(fā)生低溶解氧的危險(xiǎn)因素。

5.2 碳研所斷面發(fā)生低DO成因一方面由于氣溫升高，水體中水的飽和溶解濃度下降。但更主要的是由于水體TP、TN濃度較高，特別是TP濃度長(zhǎng)期較高，隨著氣溫升高，藻類大量繁殖，光合作用與呼吸作用增強(qiáng)，且表現(xiàn)為凈耗氧，河道DO降低。此外，枯平水期的低DI有利于藻類聚集增殖，在一定程度上增加了DO消耗。

5.3 提升DO，首先應(yīng)控制TP和TN濃度，特別是TP濃度，控制藻類增殖的物質(zhì)條件。其次，隨著AT升高，可以考慮增加上游下泄流量，提高河道DI和流速以降低水溫，降低藻類密度和活性。

5.4 采用Logistic回歸模型分析河流低DO危險(xiǎn)因素是一種便捷、可行的研究方法。本文研究結(jié)論可為釜溪河(自貢段)及我省重點(diǎn)小流域生態(tài)修復(fù)、恢復(fù)水體DO提供可行路徑。