胶州湾水体氮磷营养物基准研究❋

徐梦琪，娄 琦，崔晓莹，李正炎,2❋❋

(1.中国海洋大学环境科学与工程学院, 山东 青岛 266100; 2.中国海洋大学海洋环境与生态教育部重点实验室, 山东 青岛 266100)

海湾处于陆地和海洋相互作用强烈的区域，其生态系统较为复杂，极易受到人类活动和陆源污染的影响[1]。由于海湾地区人口分布密集，且快速的城市化和工业化进程使世界上大多数沿海水域都面临着巨大的环境压力[2-3]。过量的人为营养物，如氮(N)和磷(P)，被排放到沿海水域，导致其出现水质恶化和富营养化问题，从而引发有害藻华(Harmful algal blooms, HABs)等生态灾害[4]。《海水水质标准》(GB 3097—1997)是为中国海域的环境监测和风险评估制定的，然而，考虑到不同地区的自然特征差异，目前制定的海水营养物标准还不能为所有海域提供充分有效的保护。为了保护各海湾生态系统免受富营养化和HABs的有害影响，应针对不同海湾等沿海水域制定特定的营养物基准，以开展相应的营养物控制工作。海湾营养物基准定义为对海湾地区没有不利或有害影响的环境营养物参数的最大剂量或浓度[5]。营养物基准的制定和实施对海湾生态系统的保护、恢复和可持续发展具有重要意义。

在营养物基准研究方面，美国在过去几十年中发布了一系列针对不同水体包括河流、湖泊和水库、河口和沿海海域以及湿地[6-9]的技术指导文件。参照条件法、机械建模和压力-响应模型是制定营养物基准的三种典型方法，由美国环境保护局(the United States Environmental Protection Agency, US EPA)推荐并广泛运用[10-12]。除此之外，从水生态保护的角度，欧盟也制定了相应的技术指南使各成员国相对统一地建立营养阈值控制富营养化[13-15]。相比之下，在过去十几年中，中国在营养物基准研究方面也取得一些重要进展。但已有研究主要集中在湖泊和河口区域[16-19]，针对海湾等沿海水域的营养物基准研究仍然十分匮乏。

胶州湾被称为青岛的“母亲湾”，面临来自工业、农业、旅游、港口等各方面的环境压力。随着环湾地区高强度的人类开发活动，近几十年来N、P等营养物不断输入，并呈现不同程度的增加趋势，同时，由于湾内相对较弱的水交换能力，造成了胶州湾水体污染和富营养化等问题[20-21]。近年来，胶州湾水体呈现明显的“高N、低P、低硅(Si)”的营养盐结构，总氮(Total nitrogen, TN)和总磷(Total phosphorus, TP)的年输入通量分别为6 945和160.6 t/a[22]，且目前胶州湾N/P值呈现快速上升趋势，Si/N值则一直远低于Redfield比值[23]。由于这些变化，海湾的叶绿素a(Chlorophylla, Chla)浓度以及浮游植物的生物量都在增加，严重威胁着胶州湾的生态系统安全[22, 24]。目前，胶州湾的营养物基准尚未建立，这在一定程度上影响了胶州湾及其陆域营养物控制的科学性和富营养化风险评估。本文根据2014—2017年的胶州湾现场监测数据，采用频数分布法、分类回归树模型以及非参数拐点分析等方法推导出营养物基准候选值，并基于赤潮敏感藻生态响应，通过室内多藻种实验进一步开展营养物基准校验，最终提出胶州湾水体的营养物基准推荐值。

1 研究区域和数据来源

胶州湾(35°97′N—36°18′N, 120°04′E—120°23′E，见图1)是一个典型的半封闭型海湾，位于山东半岛南部，平均水深7 m，最大水深64 m，出现在湾口附近，面积约360 km2。胶州湾三面被青岛、红岛和黄岛环绕，且有6条主要河流汇入胶州湾，为其带来大量的点源和非点源营养物，显著影响了胶州湾的生态环境[21]。自1990年以来，胶州湾内赤潮频发，记录次数超过30次，这对青岛市社会经济的可持续发展造成了威胁[25-26]。

本研究中使用的胶州湾环境监测数据来自中国海洋大学未公开的现场调查结果和国家海洋环境监测中心的水质监测报告。现场调查在2014—2017年的3—10月进行，采样站点共22个(见图1)，监测指标包括盐度(S)、pH、溶解无机氮(Dissolved inorganic nitrogen, DIN)、活性磷酸盐(Soluble reactive phosphorus, SRP)、TN、TP、氨氮(NH3-N)、溶解氧(Dissolved oxygen, DO)和Chla。所有监测指标的样品采集、保存、运输、预处理方法、化学分析和记录规范均按照《海洋监测规范》(GB17378.4—2007)等相关标准进行。

图1 胶州湾研究区域及采样点位置

2 研究方法

2.1 营养物基准指标筛选

一般来说，营养状态变量涉及营养物指标和生物指标。营养物指标，如DIN、SRP、TN和TP，属于压力变量，能够表征水体的富营养化状况；藻类对营养物质的变化很敏感并能做出直接的反应，藻类生物量(例如Chla)和DO等生物指标则属于响应变量，反映了营养物与生物初级生产力的关系[16]。

本研究结合胶州湾2014—2017年3—10月的调查监测数据，运用相关性分析(Correlation analysis, CA)和主成分分析(Principal component analysis, PCA)，综合考虑我国海水水质评价的主要监测指标类别，进行胶州湾营养物基准指标的筛选。其中，CA是分析连续变量之间线性相关程度，根据变量数据的正态性检验，计算皮尔森相关系数(Pearson correlation coefficient)评估相关性。运用CA确定变量之间的相关性后进行PCA。PCA在尽量保证原信息量不丢失的基础上，按照方差最大化的原则，通过某种线性投影使所投影的维度上数据信息量最大，转换后的变量为主成分。CA和PCA均运用SPSS 26软件进行数据处理和统计分析。

2.2 营养物基准推导方法

2.2.1 频数分布法 频数分布法是确定水体参照状态时常用的统计分析方法，并在国内外研究中广泛应用于制定特定营养状态的营养物基准。频数分布图用于分析DIN、SRP、TN、TP、Chla和DO监测指标数据的正态性，如果数据不呈现正态分布，则采用ln进行数据转换。分析营养状态变量的频数分布曲线时，如果沿海水域受到富营养化的损害程度较小，通常提出将所有采样点数据频数分布曲线的下25%或参照点数据频数分布曲线的上25%作为营养物基准值。如果大多数采样点的水体环境因富营养化而严重退化，则应选择将监测数据频数分布曲线的上5%或下5%作为营养物基准值[6, 8]。频数分布法的分析及作图均使用SPSS 26软件进行处理。

2.2.2 压力响应模型 当无法确定受富营养化影响最小的地点及其历史数据时，可以根据现场调查的监测指标数据来分析养分浓度和初级生产力之间的压力响应关系，通过建立压力-响应模型来推导营养物基准值[10]。该模型能够定量描述营养物压力变量与生物响应变量之间的关系，尤其适用于受人为活动影响明显水域的营养物基准。为了探索压力变量和响应变量之间复杂的非线性关系，分类回归树模型(Classified and regression tree, CART)和非参数拐点分析法(Nonparametric changepoint analysis, nCPA)是两种常见的的压力-响应模型方法[19, 27]，本研究基于胶州湾2014—2017年3—10月的调查监测数据，应用CART模型和nCPA两种方法分析并进行胶州湾营养物基准候选值的推导。

CART模型可以揭示一系列响应和压力变量之间可能存在的非线性关系[28]。CART模型通过基于基尼系数的二叉树结构对压力变量进行分类，压力变量的每一次分类都是根据其最佳分类值来完成，分组后按照响应变量的最大降低目标变量的分布或变异性将其分成两个独立的节点。重复这个分类过程，直到不能进行节点分类并形成终节点为止[30]。通常，分类回归树的第一个根节点是最重要的压力指标，对响应变量有显著影响。

nCPA法是一种估计二元关系中变化点位置的方法，为营养物基准提供了自然候选值。从操作上来讲，将响应变量沿着压力梯度排序，并将拐点定义为两组响应变量之间差异达到最大值的点[10, 30]。此外，使用随机自举抽样和90%的置信区间来定义拐点位置的不确定性。CART模型和nCPA两种方法能够消除设置基准值的主观偏差，一般来说，CART模型和nCPA这两种方法得到的基准值可以相互验证，以提高确定营养物基准的准确性[19]。CART模型和nCPA法均使用R语言(R64 3.5.1)软件进行数据处理及作图。

2.3 营养物基准值验证方法

不同于单种浮游藻的验证实验，本研究考虑了浮游藻的种间竞争，更好地模拟胶州湾的水生态环境，从种群的角度推导生态响应值进行营养物基准值的验证。根据营养物水平与藻类生长之间的压力响应关系，通过混合藻种实验验证上述营养物基准候选值。验证过程包括实验数据的获取、压力响应模型的建立和生态响应值的计算。

2.3.1 混合藻种实验 近年来，中肋骨条藻(Skeletonemacostatum)和旋链角毛藻(Chaetoceroscurvisetus)一直是胶州湾的常见优势藻种，同样也是主要的赤潮藻种。胶州湾最常暴发双相型赤潮，S.costatum也最易与其他藻种同发赤潮[31-32]。S.costatum、C.curvisetus是常见的广温广盐藻种，在夏季易引发赤潮，对于这些广泛分布的赤潮藻种，应该重视其形成赤潮的潜在风险。除此之外，有研究表明[22-33]，引发胶州湾赤潮的浮游藻中，甲藻的比例明显增高，并且值得注意的是，塔玛亚历山大藻(Alexandriumtamarense)曾于夏季在胶州湾出现过[34]。由于A.tamarense为有毒赤潮藻，需要对其加强监测。因此，结合室内可培养情况，混合藻种实验选择S.costatum、C.curvisetus和A.tamarense为研究对象。

混合藻种实验进行前，将指数生长期的S.costatum、C.curvisetus和A.tamarense分别于1 700、2 000和1 000 r·min-1转速下离心10、10和5 min，去除上清液并保留下层藻液，然后用盐度为30±1的灭菌人工海水进行冲洗再离心2次，从而减少藻液中残存的营养盐。根据胶州湾水体中DIN和SRP的环境浓度水平，加入NaNO3和NaH2PO4以设置8个不同N、P浓度的实验组(见表1)，其它营养物按照f/2培养基的配方加入灭菌后的人工海水。将离心后的三种藻液以初始密度约为1×104cells/mL置于培养基中共培养。实验分为两组，即固定P浓度的N组和固定N浓度的P组，并为每个浓度水平设置三个平行样。为了使藻在生长过程中有充足的N(P)源，固定N组中N浓度及固定P组中P浓度高于实际监测浓度。不同N(P)浓度梯度则主要依据胶州湾现场监测数据中DIN和SRP的浓度范围，同时考虑到曲线拟合的相关要求进行设置。将锥形瓶置于恒温光照培养箱(GXZ-280B，宁波江南仪器厂)中培养，温度为20 ℃，光照强度为3 000～5 000 Lx，光/暗循环时间为12 h∶12 h。为防止藻类聚集，锥形瓶每天定时摇动三次。并于每天的同一时间，取1 mL藻液于取样管中，在显微镜下用血细胞计数板进行中肋骨条藻和旋链角毛藻的细胞密度计数，用浮游动物计数框进行塔玛亚历山大藻的计数。

表1 实验组中的N和P浓度

2.3.2 数据分析 藻类生长曲线通常呈S形，并可以用各种模型进行拟合，拟合模型包括Logistic、Bertalanffy、Mitscherlich和Brody[35]。在这些模型中，S型生长曲线被最广泛地用来描述生物生长过程。因此，本研究采用Logistic模型来拟合三种藻的指数生长阶段。S-logistic2(1)和S-logistic3(2)分别适用于固定P组和固定N组的生长曲线拟合，具体公式如下：

(1)

(2)

式中：Bt是t时的浮游植物生物量(cells/mL)；B0是初始浮游植物生物量(cells/mL)；Bf是终止生物量(cells/mL)；代表浮游植物种群的最大密度；R是浮游植物的种群增长率(cells·mL-1·d-1)；t是培养时间(d)；b和k是模型参数。

随着营养物浓度的变化，Bf的形状大致与S型曲线相同，并且可以用Boltzmann模型来拟合。曲线的拐点x0是Bf的突变点，可用于确定生态响应值。Boltzmann方程(3)表示如下：

(3)

式中：a,b和c为模型参数；x为营养物浓度(mg/L)。

在此混合藻种实验中，固定N(P)实验组设置为8个浓度梯度，每个浓度水平设有三个平行样。培养7天后，将每个浓度梯度的每个平行样的藻密度值作为一组。根据各组的实验测量数据(t,Bt)，通过Logistic拟合生长曲线得出Bf值。然后，进行Bootstrap随机抽样，从每个浓度梯度中随机抽取Bf值，8个浓度梯度提取的Bf值形成一组数据，总计38组数据。此外，从38组数据中随机抽取200组(x,Bf)数据，并利用Boltzmann模型进行拟合，得到200个x0值。将200个x0值的中值作为生态响应值，并使用95%的置信区间来表示平均响应值的不确定性。在获得各藻种的响应值后，选择最低值作为混合藻种实验的生态响应值，以代表最敏感藻种的响应。使用Matlab 2019b和OriginPro2021软件进行Bootstrap随机抽样和模型拟合。

3 结果

3.1 营养物基准指标的选择

从各监测指标的相关性系数矩阵和显著性检验结果(见表2)中可以直观地看到，各监测指标之间是具有明显相关性的。其次，表3是PCA结果，包括各成分的特征值、成分贡献率和累积贡献率。成分1、2、3的特征值分别为2.997、1.748、1.124，均大于1且其累积贡献率达到73.361%，因此提取前三个成分为主成分。此外，通过各因子旋转形成旋转后的成分矩阵(见表4)来分析与各主成分密切相关的监测指标信息。

表2 各指标的相关系数矩阵与显著性检验结果

表3 特征值、成分贡献率与累积贡献率

表4 旋转后的成分矩阵

PCA结果显示，DIN、SRP、TN和TP为第一主成分相关的指标，第二主成分相关指标为Chla和S，pH和DO则为第三主成分相关指标。DIN和SRP通常用于评估富营养化状况，考虑我国《海水水质标准》(GB 3097—1997)中的主要监测指标类别，将DIN和SRP选为主要的营养物指标，TN、TP为辅助性营养物指标。除此之外，由于S与DIN、SRP、TN和TP的相关性较弱，且pH容易受藻类生长、海水扰动等作用影响而产生较大波动，因此不考虑将S和pH作为营养物基准的指标。而CA结果表明DO与DIN和SRP密切相关，且DO是生态健康的重要指标，同时也是水体营养状况的早期预警指标。除此之外，Chla含量表明浮游植物的生物量，是HABs最具代表性的响应指标[36-37]。随着营养物浓度的升高，浮游植物密度增加，DO含量减少。因此，DIN、SRP、TN和TP被选为压力变量，Chla和DO被选为响应变量。

3.2 频数分布法

我国《海水水质标准》(GB3097—1997)将水质分为四类，通过对胶州湾的监测数据进行统计分析，发现76.7%采样点的DIN和98.3%采样点的SRP属于前三类水质，同时考虑到胶州湾整体的富营养化状况，本研究基于胶州湾2014—2017年3—10月的调查监测数据，选择将所有采样点数据的频数分布曲线的下25%来推导胶州湾营养物基准的候选值。频数分布曲线及其统计值如图2和表5所示，DIN、TN、SRP和TP的营养物基准候选值分别为0.152、0.394、0.003 5、0.008 9 mg/L(见图2和表5)。

图2 胶州湾营养状态指标的频数分布曲线

表5 胶州湾营养状态指标的频数分布统计结果

3.3 分类回归树模型

以Chla和DO为响应变量，DIN、SRP、TN和TP为压力变量，使用CART模型建立压力响应关系，选择CART根节点的养分浓度作为胶州湾营养物基准候选值，并将DIN和TN的结果保留三位有效数字，以Chla和DO为响应变量的CART分析结果如图3、图4和图5所示。从图3中可以看出，SRP是影响Chla和DO浓度变化的最首要且最重要的压力变量，这与目前胶州湾存在的P限制现象相符合。以Chla为响应变量，得出DIN、TN、SRP和TP的营养物基准候选值分别为0.257、0.410、0.003 2和0.009 8 mg/L(见图4)；将DO作为响应变量，计算出DIN、TN、SRP和TP的营养物基准候选值分别为0.235、0.885、0.010 9和0.035 5 mg/L(见图5)。

(箱图是根据分类数据中Chl a和DO的平均浓度升序排列，表示每个终端节点上Chl a和DO的浓度分布。The boxplots represent the concentration distributions of Chl a and DO at each terminal nodes, arranged according to the average Chl a and DO concentrations in the classified data in an ascending order.)

(箱图是根据分类数据中Chl a的平均浓度升序排列，表示每个终端节点上Chl a的浓度分布。The boxplots represent the concentration distributions of Chl a at each terminal nodes, arranged according to the average Chl a concentrations in the classified data in an ascending order.)

(箱图是根据分类数据中DO的平均浓度升序排列，表示每个终端节点上DO的浓度分布。The boxplots represent the concentration distributions of DO at each terminal nodes, arranged according to the average DO concentrations in the classified data in an ascending order.)

3.4 非参数拐点分析法

以DO为响应变量，DIN、SRP、TN和TP为压力变量，使用nCPA方法建立压力响应关系。将DIN、SRP、TN和TP的拐点作为这些变量对DO的突变点，并将DIN和TN的结果保留三位有效数字。nCPA的结果显示(见图6)，DIN、TN、SRP和TP的营养物基准候选值分别为0.235、0.885、0.010 9和0.035 5 mg/L。由于Chla数据不足，未对该响应变量进行nCPA分析。

(折线表示响应变量在拐点处的阶梯变化，虚线表示拐点90%的置信区间，圆圈表示样本数据。The broken line shows the stepped variation of the response variable at the changepoint, the vertical dashed lines show the 90% confidence intervals of the changepoint, and the circles represent sample data.)

3.5 混合藻种实验

C.curvisetus和A.tamarense的生长曲线与总生物量有很大不同，但S.costatum的生长曲线与总生物量相似(见图7、8)。在固定N和P组中，S.costatum为优势种，占总生物量的大部分，且由于营养物质丰富，S.costatum和总生物量在0～3 d和0～4 d阶段的生长曲线基本呈S形，因此，选择固N组0～3 d和固P组0～4 d的阶段进行Bf值的拟合。S.costatum的DIN和SRP生态响应值分别为0.521和0.044 6 mg/L，总生物量的DIN和SRP生态响应值分别为0.642和0.047 4 mg/L(见表6)。因此，对DIN和SRP最敏感藻种的生态响应值分别为0.521和0.044 6 mg/L。根据胶州湾的监测数据，SRP∶TP和DIN∶TN比值的平均值分别为0.315和0.519，TN和TP的生态响应值则分别为1.004和0.142 mg/L。

图7 固N组中肋骨条藻、旋链角毛藻、塔玛亚历山大藻和总生物量的生长曲线

图8 固P组中肋骨条藻、旋链角毛藻、塔玛亚历山大藻和总生物量的生长曲线

表6 DIN和SRP的生态响应值

3.6 营养物基准推荐值

通过频数分布法、分类回归树模型和非参数拐点分析法得出的营养物基准候选值和通过混合藻种实验计算出的生态响应值如表7所示。鉴于对胶州湾现场监测数据的统计分析，76.7%采样点的DIN和98.3%采样点的SRP属于前三类水质，综合考虑胶州湾目前的水质状况和富营养化水平，选择所有方法得出的基准最低值作为胶州湾的营养物基准推荐值，以此对胶州湾的营养物控制和评价提出更高的标准和要求。因此，本研究建议胶州湾的DIN、TN、SRP和TP营养物基准值分别为0.152、0.394、0.003 2和0.008 9 mg/L。

表7 不同方法推导的营养物基准候选值

4 讨论

由于野外调查和室内实验的环境条件不同，营养物基准候选值与生态响应值相差较大(见表7)。沿海地区N、P的浓度会受陆源输入、S、水停留时间、海水交换、沉积特征、海洋环流和大气沉降等的影响[38-39]。因此，基于现场监测数据的营养物基准候选值也会受到上述因素的影响。基于营养物压力指标和生物响应指标之间的关系，混合藻种验证实验模拟了浮游藻最佳生长条件，为藻类生长提供了固定的温度、S和光照条件。然而，实验室无法完全模拟出海湾自然环境中的水文状况和生物环境，使得营养物质的浓度和组成结构是影响生态响应值的主要因素。胶州湾目前呈现高N、低P的营养盐结构，因此，本研究提出的DIN和SRP营养物基准推荐值的N/P值相应较高。为了避免过高的N/P值可能会导致生态系统失衡的结果，在后续的研究中，需要进一步进行现场微宇宙实验，综合考虑水文条件、浮游动物等胶州湾自然生态环境条件的影响，来验证理论推导方法得出的营养物基准候选值，同时结合较为适合的N/P值进行营养物基准候选值的校正，在实现有效控制胶州湾N、P浓度的前提下，进一步保证胶州湾生态系统的平衡和稳定。

从表8可以看出，DIN和SRP的营养物基准值在不同的研究中差异较大，这可能是由于各个国家和地区沿海海域富营养化状况的差异以及对营养物质的敏感性不同造成的。此外，推导方法选择的不同，也是导致营养物基准值存在一定差异的原因。因早期辽河口水质状况较为良好，除辽河口沿海地区外，胶州湾的DIN基准值低于其他沿海地区。而不同沿海地区的SRP基准值大约是胶州湾的2～9倍，这可能是与胶州湾在2010年以来存在的P限制现象有关[23]。本研究主要针对DIN和SRP指标，而TN和TP仅作为辅助性指标。在现行的《海水水质标准》(GB3097—1997)中，水质监测和评价的营养物指标也主要针对DIN和SRP，而TN和TP主要作为淡水中的营养物基准指标。不同研究中TN和TP的营养物基准值如表9所示，由于不同国家和地区藻类对营养物质的敏感性差异，除长江口外及舟山邻海、佛罗里达群岛外，胶州湾推荐的TN基准值低于表9中其它地区中营养物基准的范围。而胶州湾推荐的TP基准值低于或在其他研究的范围内，但比斯坎湾除外。

表8 不同沿海地区DIN和SRP的参照条件或营养物基准

表9 不同国家和地区中TN和TP的营养物基准

此外，有研究表明胶州湾在20世纪70年代以前，受人为影响很小，N、P等营养物基本呈现背景值且较为稳定，其中，SRP浓度在20世纪60年代约为0.005 4 mg/L。自20世纪80年代起，胶州湾DIN和SRP浓度呈现波动变化趋势，目前胶州湾DIN浓度略高于国家Ⅰ类海水水质标准，而SRP浓度基本符合国家Ⅰ类海水水质标准[22-23]。本研究提出的DIN和SRP营养物基准推荐值均低于我国《海水水质标准》(GB 3097—1997)的Ⅰ类。TN的营养物基准推荐值在我国《地表水环境质量标准》(GB 3838—2002)的Ⅰ类和Ⅱ类范围内，TP的营养物基准推荐值低于地表水水质的Ⅰ类。总体而言，本研究中建议的胶州湾DIN、TN、SRP和TP营养物基准值是相对合理可靠的，能对胶州湾实现有效保护。

5 结论

(1)根据2014—2017年的胶州湾调查监测数据，运用相关性分析和主成分分析选择DIN、SRP、TN和TP作为营养物指标即压力变量，Chla和DO作为生物指标即响应变量。

(2)运用频数分布法、分类回归树模型以及非参数拐点分析三种方法推导出DIN、TN、SRP和TP营养物基准候选值，然后通过混合藻种实验计算出生态响应值开展进一步验证，最终提出胶州湾DIN、TN、SRP和TP的营养物基准推荐值分别为0.152、0.394、0.003 2和0.008 9 mg/L。

(3)通过比较不同国家和地区的营养物基准研究，本研究中建议的DIN、TN、SRP和TP营养物基准值合理可靠，能够为胶州湾的水质监测及评估等环境管理工作提供科学支撑。