基于Logistic回归模型的低溶解氧风险因素与成因分析——以釜溪河(自贡段)为例

魏 峣，陈卓尔，杨明方，佟洪金，刘晓聪

(1.四川省生态环境科学研究院，成都 610041；2.四川省环保科技工程有限责任公司，成都 610041)

前 言

水体溶解氧(Dissovled Oxygen，DO)是水生生物生理代谢和鱼类生长繁育的限制因子之一，是衡量水体自净能力的一个重要指标，也是水体水质的重要指标之一[1～3]。常温下水中DO含量一般为8～14mg/L，DO低于5 mg/L时便对水生生物的生存产生不利影响，低于2 mg/L时大多数鱼类可能死亡[4]。水体DO的变化受到光合作用、呼吸作用、大气复氧过程和沉积物耗氧等过程的综合影响[5～7]。其中，光合作用是水体溶解氧最重要的来源。欧阳潇然[8]对太湖研究发现，蓝藻水华易发生和堆积的太湖梅梁湾，浮游植物光合作用制氧是梅梁湖水中最重要的来源；浮游植物的呼吸与死亡耗氧是溶解氧的最大的消耗者。与此相对应，水体中DO的变化受温度、pH、高锰酸盐、悬浮物、营养物质、生物量等多种因素影响[9]。

Logistic回归分析是研究因变量是分类、定性变量的常用方法，常用于流行病学危险因素分析[10-11]，可用于判别引发事件的因素[12]，预测事件发生的概率[13]。近年来在自然灾害[14-15]、环境[16-18]等领域也有广泛应用。采用Logistic回归模型分析水质、水文、气象三方面因素对釜溪河碳研所断面发生低DO影响程度，识别引发低DO的危险因素，分析低DO成因，对开展河道水生态治理修复具有重要作用。

1 研究区域

自贡市地处四川盆地南部，位于北纬28°55′37″～29°38′25″、东经104°02′57″～105°16′11″之间。地势西北高东南低，地形分为低山、丘陵、平坝，其中丘陵分布最广，占全市总面积的80%以上。属四川盆地亚热带季风湿润气候区，多年平均气温17.9℃，多年平均降雨量1 079.6mm，多年平均相对湿度80%，多年平均蒸发量1 010.3mm。全市土地面积4 381km2，常住人口290万人，地区生产总值1 312亿元。

釜溪河是长江二级支流、沱江一级支流。发源于四川省内江市威远县两母山东北麓，由旭水河和威远河在自贡双河口处汇成干流，流经自贡市主城区，在富顺李家湾汇入沱江(图1)。釜溪河流域面积3 472km2，长度190km，平均比降2.97‰，多年平均流量42.4m3/s，多年平均径流量13.4亿m3/年。国控碳研所断面位于釜溪河自贡市主城区下游约1.5km。该断面水质长期较差，被列为《长江流域劣Ⅴ类国控断面整治专项行动工作方案》(2019年)中四川省两个劣Ⅴ类国控断面之一。除氨氮、总磷、化学需氧量等超标外，DO浓度是釜溪河流域最低点，春末及夏季容易出现低DO现象，影响水生态健康。

图1 研究区域Fig.1 Study area

2 数据来源与方法

2.1 数据来源

收集了2017～2018年釜溪河自贡水文站日均流量(Discharge，DI)；自贡气象站日均气温(Air Temperature，AT)；碳研所水质自动站每日6时点(0、4、8、12、16、20时)DO、总氮(Total Nitrogen，TN)、总磷(Total Phosphorus，TP)小时数据，并计算得到日均值(部分缺失时点数据采用插值方式补齐)。由于部分天数水质数据缺失，研究共收集到444日数据(2017年共283日，2018年共161日)。

2.2 研究方法

Logistic回归模型多采用二分类变量，设P为某事件(Y=1)发生的概率，取值范围为[0，1]，(1-P)为该事件不发生(Y=0)的概率。将P/(1-P)取自然对数lnP/(1-P)，记为logitP，则logitP的取值范围为(-∞，+∞)。以P为因变量，m个危险因素为自变量，建立线性模型(式1)。

(1)

式(1)可变式为概率模型形式(式2)。

(2)

模型中β0为常数，表示在不接触任何危险因素条件下，结果发生与不发生的概率之比的对数值(基线值)。系数βi是某危险因素Xi存在(Xi=1)与不存在(Xi=0)相比，事件发生的优势比(Odds Ratio，OR)的对数值。

βi值>0，eβi>1，表示有危险因素对事件发生的优势大于无危险因素；

βi值<0，eβi<1，表示有危险因素对事件发生的优势低于无危险因素；

βi值=0，eβi=1，表示无论危险因素是否出现，与事件发生与否无关。

他们披头散发，蓬头垢面，日复一日，年复一年，在暗无天日的井洞里，顽强地坚持着。这样的作业空间，这样危险的工作环境，每一个采矿者的生命，随时都会从阳间一步踏过奈何桥，走向死亡。生命之于他们，只是在阳世间的短时间寄存。每一天清晨他们进入矿洞，便是生死两茫茫，晚上，不知道还能不能回到栖息的工棚。

3 结果与分析

3.1 基本特征

收集的444个日均值数据统计结果(表1)，DO平均浓度4.0±1.9 mg/L，最低浓度仅0.5 mg/L，出现在2017年6月。TP平均浓度0.384±0.246 mg/L，最高浓度达1.612 mg/L，出现在2018年3月。TN平均浓度7.1±2.9 mg/L，最高浓度21.5 mg/L，出现在2018年2月。平均AT为15.6±6.9℃，最高温31.6℃，出现在2018年8月。平均DI为3.4±14.1 m3/s，最低DI为0.1 m3/s，出现在2017年11月。

表1 各指标特征Tab.1 Characteristics of variables

分析各指标相关性(表2)，DO与TP、TN、AT呈显著负相关关系(P<0.001)，表明较高的TN、TP浓度及AT对应较低的DO。DO与DI(P=0.008)呈显著正相关关系，表明较低的DI对应较低的DO。

表2 各指标相关性分析(n=444)Tab.2 Correlation analysis of variables (n=444)

3.2 危险事件及危险因素判别

以地表水Ⅲ类水质限值为判别标准，将DO<5mg/L定义为发生低溶解氧(Y=1)，DO≥5mg/L定义为未发生低溶解氧(Y=0)。在此判别标准下有65.5%的样本发生低溶解氧(表3)。

表3 危险事件判别标准(n=444)Tab.3 Criteria for risk events

Logistic回归分析需要将自变量进行二分类取值，即选取某一界限值将自变量分为0和1两组，其中取1表示为危险因素，取0表示非危险因素。根据各指标在Y为0和1两种情形下的累积频率分布或水质标准，划分自变量危险因素(表4)。

表4 各指标危险因素判别界限取值Tab.4 Criteria for risk factors of variables

经检验，各指标方差膨胀因子(Variance inflation factor，VIF)均小于10(表5)，表明选取的指标不存在明显的多重共线性问题[19]。

表5 各指标方差膨胀因子Tab.5 VIF of variables

3.3 模型建立

采用方差分析(Analysis of Variance，ANOVA)对比各指标分组间的差异(表6)。TP、TN、AT和DI危险因素组与非危险因素组的DO存在显著差异(P<0.001)。

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

采用Logistic模型进行多因素回归分析，得到各危险因素的回归系数(βi)估计值，建立回归模型。模型Hosmer-Lemeshow统计量P=0.459，大于0.05，表明模型预测值与样本之间不存在显著差异[20](表7)，即模型拟合较好。各回归系数的Wald检验P值均小于0.001，表明模型回归系数显著(表8)。

表7 Hosmer-Lemeshow 检验结果(n=444)Tab.7 Result of Hosmer-Lemeshow test (n=444)

表8 各指标回归系数估计值(n=444)Tab.8 Regression coefficients of variables(n=444)

建立的Logistic回归模型见式(3)。

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

(3)

根据预测结果，模型对DO<5mg/L的样本预测正确率86.6%，对DO≥5 mg/L的样本预测正确率73.2%，总体预测正确率82%(表9)。

表9 预测结果(n=444)Tab.9 Predicted results(n=444)

3.4 危险因素识别

各指标OR值均>1，表明TP≥0.2mg/L、TN≥6.4mg/L、AT≥14.2℃、DI≤0.3 m3/s均是发生低溶解氧的危险因素(表10)。在控制其他变量后，发生低溶解氧的概率，AT≥14.2℃是<14.2℃的37.607倍，较高的AT是发生低溶解氧的首要危险因素；TP≥0.2mg/L是<0.2 mg/L的11.031倍，TN≥6.4mg/L是<6.4 mg/L的4.358倍，较高的N、P浓度是发生低溶解氧的重要危险因素；DI≤0.3 m3/s是>0.3 m3/s的2.962倍，较低的流量是发生低溶解氧的危险因素。

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

3.5 各月DO的主要危险因素

计算各指标月平均值并按危险因素分类(图2)。各危险因素中，较高AT(≥14.2℃)分布在3～10月，较高TP(≥0.2mg/L)分布于全部12个月，较高TN(≥6.4mg/L)分布在1～5月和12月，较低DI(≤0.3m3/s)分布在3、5月。时间尺度上影响最广泛的是较高TP和较高AT。

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

月均DO<5mg/L的月份为2～7月和9～11月，共9个月。2月主要危险因素是较高的TP和较高的TN；其余8个月危险因素是较高的TP和较高的AT，其中，3、5月还受低流速影响。

4 讨 论

AT对DO的影响主要包括两个方面。一是气体溶解在液体中通常为放热过程，水的饱和溶解浓度随水温升高而降低[21]。二是水温通过影响藻类营养盐摄入、细胞中酶的活性从而影响藻类生产力[22]。水温上升，藻类生长速度加快[23]。赵巧华[24]等研究认为水温是影响太湖春季藻类生长的关键因素之一。水体DO变化与藻类的光合作用、呼吸作用有关[25]。当藻类数量增加到一定数量级时，其数量多少和生命活动的旺盛程度控制了水体中DO的变化[26]。无光照时，藻类因内源呼吸，吸收氧气放出二氧化碳；有光照时，藻类光合作用放出氧气。随气温升高，藻类白天的净释氧量和夜间的净耗氧量均增加，日平均结果表现为DO的净消耗，DO浓度逐渐降低。碳研所断面DO与AT呈显著负相关关系(P<0.001)，DO随AT升高而降低，一方面与AT升高导致饱和溶解氧降低有关。更重要的是，从4月起日均AT均在14.2℃以上，较高的温度一直持续到10月。相应月份的水温多在20℃以上，这有利于藻类繁殖和活性增强，呼吸作用耗氧大于光合作用释氧，导致DO逐渐降低。

N、P营养盐与水体中生物量存在紧密关联[27]。水体中N、P浓度将通过藻类繁殖及光合作用、呼吸作用影响水体DO。过量的N、P输入将导致水体富营养化。富营养化水体在适宜的环境条件下可引起藻类大量繁殖，消耗水体中溶解氧[28]。一般认为水体富营养化的条件为TN≥0.2mg/L，TP≥0.02mg/L[29]。釜溪河自贡城区段TN平均浓度7.1±2.9mg/L， TP平均浓度0.384±0.246mg/L，属于富营养化水体。在适宜的环境条件下可能引发藻类大量繁殖，导致DO降低。

一定河道断面下，流量与流速呈正相关。流速(流量)主要影响水体动力学条件。湖泊或者类湖泊型形态水体一般流速较小，水力停留时间长，给藻类营造一种稳定的水力环境。较快的流速，水体更新周期更快，不容易造成单一藻类蓄积和过量增殖[30]。自贡水文站流量≤0.3m3/s多分布在12月至次年6月的枯平水期，较小的流速有利于藻类的聚集增殖，在一定程度上增加了DO消耗。

5 结 论

5.1 TP≥0.2mg/L、TN≥6.4mg/L、AT≥14.2℃、DI≤0.3 m3/s是碳研所断面发生低溶解氧的危险因素。其中，较高的AT是发生低溶解氧的首要危险因素，较高的N、P浓度是发生低溶解氧的重要危险因素，较低的DI是发生低溶解氧的危险因素。

5.2 碳研所断面发生低DO成因一方面由于气温升高，水体中水的饱和溶解浓度下降。但更主要的是由于水体TP、TN浓度较高，特别是TP浓度长期较高，随着气温升高，藻类大量繁殖，光合作用与呼吸作用增强，且表现为净耗氧，河道DO降低。此外，枯平水期的低DI有利于藻类聚集增殖，在一定程度上增加了DO消耗。

5.3 提升DO，首先应控制TP和TN浓度，特别是TP浓度，控制藻类增殖的物质条件。其次，随着AT升高，可以考虑增加上游下泄流量，提高河道DI和流速以降低水温，降低藻类密度和活性。

5.4 采用Logistic回归模型分析河流低DO危险因素是一种便捷、可行的研究方法。本文研究结论可为釜溪河(自贡段)及我省重点小流域生态修复、恢复水体DO提供可行路径。