基于非线性规划的莱州湾营养盐环境容量计算*

虞兰兰,陈艳丽,张海宁,杨　波,江文胜，虞　洋

(1.国家海洋局 北海信息中心，山东 青岛 266061；2.中国海洋大学 物理海洋实验室，山东 青岛 266100；3.海军北海舰队海洋水文气象中心，山东 青岛 266003)



基于非线性规划的莱州湾营养盐环境容量计算*

虞兰兰1,陈艳丽1,张海宁1,杨波2*,江文胜2，虞洋3

(1.国家海洋局 北海信息中心，山东 青岛 266061；2.中国海洋大学 物理海洋实验室，山东 青岛 266100；3.海军北海舰队海洋水文气象中心，山东 青岛 266003)

摘要:在前人研究的基础上提出了基于非线性规划的排海通量最优化法(简称非线性规划法)，以进行海域内非保守物质环境容量计算，并选取莱州湾这个水动力较强、面积较大的半封闭海域进行模拟、计算。为更好地保护生态环境，对于氮(N)、磷(P)等非保守物质(以无机氮(DIN)、无机磷(DIP)为例进行分析)，选取其月最高浓度场数据，依次计算出莱州湾海域每个月的营养盐响应系数场，并进一步使用非线性规划法求解出营养盐的月环境容量。最终计算结果显示， 21世纪初莱州湾的营养盐环境容量为：在保持莱州湾DIP排放现状不变的前提下，其DIN年环境容量为21 843 t/a，剩余年环境容量为-8 264 t/a，需要减少点源DIN的排放量；而在保持莱州湾DIN排放现状不变的前提下，其DIP年环境容量为1 794 t/a，剩余年环境容量达1 485 t/a，可继续容纳部分DIP的排放。

关键词：营养盐；莱州湾；非保守物质；环境容量；非线性规划

对海域污染物的环境容量的研究是由对河流、湖泊污染物环境容量的研究拓展到海洋领域所形成的。环境容量最初只应用于国外海域污染物的排海通量控制[1-4]，后被我国学者借鉴和应用，并不断得以发展[5-8]。污染物在海洋环境中发生的物理、化学、生物等自净过程共同决定了其环境容量的大小。目前，在我国应用比较广泛的环境容量计算方法是，综合考虑各种自净过程、基于线性规划的排海通量最优化法(简称线性规划法)[9-11]。该方法的计算原理是通过优化各污染源的排海通量分配率，使各排污口在满足一定等级海水质量标准的前提下，其污染负荷排放量总和达到极大值。

对于氮(N)、磷(P)营养盐、石油烃等非保守物质，线性规划法的前提假设不能得到满足[12]，经改进目前针对非保守物质，较为可靠的环境容量计算方法为基于非线性规划的排海通量最优化法(简称非线性规划法)[12-13]。本文主要将非线性规划法应用于莱州湾海域N、P营养盐的环境容量计算。

1莱州湾简介

渤海是我国唯一的封闭性内海，莱州湾则位于渤海南部的浅海湾，湾内地势平坦，总面积达6 966.93 km2，水深一般为10～15 m，最深处约为20 m(图1)。黄河、小清河、虞河、界河等河流携带了大量的无机营养盐流入莱州湾，有利于浮游植物的繁殖与生长，使莱州湾形成了黄、渤海水产资源的主要产卵场、栖息地和多种渔业的传统渔场，是我国北方重要的渔业资源基地[13-14]。

近年来，随着莱州湾沿岸工、农业的发展和人口的增长，入海河流携带的营养盐不断增多，引起了湾内生态环境的恶化，这已经成为莱州湾目前面临的重要问题。莱州湾海域生态环境的恶化表现为：湾内海水的富营养化不断加重、赤潮频发，产卵场受到破坏，主要经济鱼类和渤海对虾资源严重衰退。为更好地控制污染物的排放量，减小污染物对莱州湾海域生态系统的危害，对其进行污染物环境容量的计算很有必要。

图1　莱州湾的营养盐点源及水质浓度控制点分布图Fig.1　The distribution of main nutrients sources points and water quality control points in Laizhou Bay

2非保守物质环境容量计算的简要说明

为实现排海通量最优化法的计算，综合应用了2种数值模式，包括HAMSOM水动力模式和基于PIC(Particle-In-Cell)方法的三维生态模式。首先，通过HAMSOM水动力模式模拟出莱州湾的水动力场，包括海流、水位、温度、盐度等要素值；再将各水动力要素值保存并导入三维生态模式中，用来模拟莱州湾的N、P营养盐和浮游植物等的分布情况；利用三维生态模式依次模拟出各点污染源的响应系数场，进而计算出莱州湾N、P营养盐的环境容量。

三维生态模式的空间分辨率为1′×1′，时间分辨率为30 min，模式在垂直方向上采用Z坐标的形式，垂向共分为7层，1～6层层厚均为3 m，第7层层厚5 m。本文各项环境容量工作的开展均基于该生态模式,模式包括营养盐、浮游植物、浮游动物以及碎屑四个模块，涉及的生态过程与生态参数设置等及其验证模拟最终结果可参考文献[12]。为验证模式的稳定性设计了一个数值试验，将河口和开边界的COD浓度都设为1.2 mg/L，并将COD的初始场也设为1.2 mg/L。模拟100 d以后，全场COD的浓度仍然保持1.2 mg/L不变(图略)，这表明生态模式是质量守恒的。

N、P营养盐之间相互影响、此消彼长，N排放量的增大会相应导致P被浮游植物更多地吸收和消耗，从而减小P在海域中的浓度值，反之亦然。故当N(或P)的浓度处于不同状态时，相应的P(或N)的环境容量并不相同。由于N、P营养盐之间相互作用较为复杂，而本研究旨在探讨非保守物质环境容量的计算方法，所以采用简化的手段，将二元非线性问题简化为一元非线性问题：保持一种营养盐按照现状排放不变，计算另一种营养盐的环境容量。

3非线性规划法的应用前提

由于N、P营养盐具有明显的季节变化特征，随时间是不断改变的，不可能随时间的变化逐渐趋于一个稳定不变的状态，所以无法得出最终稳定的响应系数场，应根据实际需求选取响应系数场。为尽可能避免对湾内营养盐环境容量的高估，及更好地保护生态环境，本文采用每月最高浓度场来计算响应系数场[12]。而在使用非线性规划法之前，首先应该验证该方法的前提条件是否能够得到满足，即点污染源是否满足线性叠加。将采用湾内营养盐的月最高浓度场，按月验证非保守物质是否满足点污染源的线性叠加。由于模拟、计算时选用21世纪初(2000-2008年)的气候态平均数据，并且4月份是浮游植物全年生物量的高峰期，营养盐分布情况最复杂，所以选取21世纪初的4月份进行气候态特征的分析。

3.1点污染源线性叠加的验证

点源的响应系数场是单独改变某点源的污染物排放量计算得出的，但环境容量最优化计算时针对的是所有点源同时改变的情况，所以各点源之间必须满足线性叠加关系，即各点源营养盐排放量依次单独改变与同时改变的效果等价，才可将各点源单独排放所得出的响应系数场进行线性叠加，最终得出各点源同时改变所形成的浓度场。

已有研究直接使用线性规划法计算非保守物质的环境容量，并未对各点源之间是否满足线性叠加关系进行验证，本研究将以DIN为例补充完成这一工作。保持各点源的DIP排放现状不变，根据莱州湾DIN浓度的变化分析各点源之间的线性叠加性，取浮游植物生长的峰值月份4月进行研究(图2)。

图2　4月份5个点源采用不同方式改变DIN排放量时莱州湾DIN浓度场改变量的结果对比Fig.2　The comparison results of the DIN variation of concentration field in Laizhou Bay with different way of DIN discharge in five source points in April

由4月份5个点源采用不同方式改变DIN排放量时莱州湾DIN浓度场改变量的结果对比(图2)显示可知，各点源同时改变得到的莱州湾营养盐浓度场改变量与各点源依次单独改变所形成的营养盐浓度场改变量之和的分布在整个莱州湾基本一致，只在南部湾底的小部分海域其线性叠加性稍弱。由此可见，莱州湾各点源之间基本符合线性叠加关系，满足使用非线性规划法的前提条件。故可以使用非线性规划法进行营养盐的环境容量计算。

3.2点污染源响应系数场的计算

由于月最高浓度场计算出的响应系数场包含非线性部分，所以必须选用非线性规划法来进行非保守物质环境容量的计算。由于非线性规划中约束条件的限制不能太强，否则数值上算不出最优解，所以不能采用线性规划选取9～16个水质控制点的方式。本研究在每个点源附近距离1个网格的所有水点中，弃掉浓度极高或者极低的水质点，选取剩余水质点中浓度相对较高的点作为水质控制点，以避免过高或过低估计污染物的环境容量。

如图1所示，在莱州湾5个点污染源(r1到r5)附近对应选取的5个水质控制点依次用P1，P2，P3，P4，P5表示。分别对各水质控制点进行各点源DIN响应系数的函数拟合，具体拟合过程以图3为例进行说明。图3中自变量ΔF代表点源的DIN排放改变量，其单位F0为各点源相应的初始排放量，F0的值对于不同的点源并不相同，R1为一阶拟合的相关系数，R2为二阶拟合的相关系数。优先考虑相关性最好的多项式拟合；若高阶拟合与低阶拟合的相关系数比较接近时，则选用低阶拟合，以降低非线性最优化求解的难度。

黄河、小清河、潍坊直排口、虞河和界河DIN排放量的改变量分别用r1，r2，r3，r4和r5表示；对于某水质控制点(P1)，由各点源引起的DIN浓度场的改变量则分别用C1，C2，C3，C4和C5表示由各点源引起的DIN浓度场改变量；将浓度场的改变量统一用ΔC进行标示。

图3　4月份P1处DIN浓度改变量ΔC相对于各点源DIN排放改变量ΔF的相应系数场拟合Fig.3　The functional simulation of ΔC and ΔF of the water quality control point P1 in April

4月份水质控制点P1处选取的各点源DIN响应系数拟合函数：

C1=27.331r1+3.112;

(1)

(2)

C3=0.008r3+3.907;

(3)

C4=0.019r4+3.907;

(4)

C5=3.907。

(5)

按照上述方法，依次按月对各水质控制点进行各点源DIN响应系数场的拟合，结果显示，基本上每个水质控制点的5个拟合函数中都至少有1个函数的自变量为二次项(即非线性项)，甚至有些函数的自变量出现了三次项。由此亦可再次证明使用非线性规划法进行莱州湾营养盐的环境容量计算的必要性。

4非线性规划法对营养盐环境容量的计算

本研究将选用基于营养盐月最高浓度场数据的非线性规划法计算N、P等非保守物质的环境容量，这样虽然会低估海域的环境容量，然而低估部分可以作为安全保证额，使莱州湾海域对于营养盐具备更好的缓冲能力，从而更好地维持海域的生态平衡[12]。

4.1非线性规划法的提出与应用

根据渤海近岸海域主要海洋功能区划的要求，莱州湾内5个点源附近海域都应满足二类水质要求，(http:∥www.soa.gov.cn/zwgk/fwjgwywj/gwyfgwj/201211/t20121105_5255.html)，即选取的各水质控制点处，DIN和DIP浓度都必须满足二类以上水质标准。根据国家海水水质标准的规定，DIN 的二类水质标准浓度为21.429 mmol/m3，DIP的二类水质标准浓度为0.968 mmol/m3。

可将莱州湾营养盐的环境容量计算可归纳为一个非线性规划问题[9-12]。以保持莱州湾DIP排放不变、求解湾内DIN的环境容量为例进行说明，其目标函数为：Max(r1+r2+r3+r4+r5);由于Matlab要求用极小值表示目标函数，故上式改为：Min-(r1+r2+r3+r4+r5)。

边界约束条件为：-1≤r1≤0;-1≤r2≤0;r3≥0；r4≥0;r5≥0。

该非线性规划不具备等式约束条件。不等式约束条件为水质控制方程，在不同的月份，不等式约束条件是不同的。以4月份黄河点源附近水质控制点P1的DIN浓度不等式约束为例：

(0.019r4+3.907)+3.907≤21.429。

(6)

选用Matlab提供的非线性最优化求解函数fmincon来求解非线性规划问题，从而得出各点源的最优排放量。

在保持莱州湾DIP排放量不变的基础上，湾内DIN年环境容量为21 843 t/a，剩余年环境容量为-8 264 t/a，应减少莱州湾点源的DIN排放量。由图4和图5可见，DIN的月环境容量在夏季和秋季最大，在冬、春季最小；由于黄河的DIN排放量超标太多，导致整个莱州湾的剩余环境容量全年都为负值(即需要削减DIN的排放)，其中秋季需要削减的DIN排放量最大，冬、春季需要削减的量最小。

图4　莱州湾DIN的月环境容量Fig.4　The monthly environment capacity of DIN in Laizhou Bay

图5　莱州湾DIN的剩余月环境容量Fig.5　The monthly surplus environment capacity of DIN in Laizhou Bay

在保持莱州湾DIN排放不变的基础上，湾内DIP年环境容量为1 794 t/a，剩余年环境容量达1 485 t/a，可继续容纳部分DIP的排放。由图6和图7可见，DIP的月环境容量在夏季和秋季最大，在冬、春季最小。由于各点源的DIP排放量较小，莱州湾的剩余环境容量全年都为正值(即不需要削减DIP排放量，湾内仍可继续容纳部分DIP的排放)，并且夏季和秋季的剩余环境容量最大，冬、春季的最小。

图6　莱州湾DIP的月环境容量Fig.6　The monthly environment capacity of DIP in Laizhou Bay

图7　莱州湾DIP的剩余月环境容量Fig.7　The monthly surplus environment capacity of DIP in Laizhou Bay

4.2非线性规划法计算结果——点源排放优化量

以DIN为例对非线性规划计算结果进行说明。为更清楚地了解各点源的DIN排放现状与剩余排放量之间的关系，将5个点源的情况列举在表1中。削减率表示在现状排放的基础上，需要削减的排放量占现状排放量的百分比，正值表示需要削减，负值表示可容纳更多DIN营养盐的排放。

表1　莱州湾各点源DIN年排放现状与剩余年排放量的关系

由表1可知，黄河和小清河的DIN年排放量都已超标，分别需要削减其年排放量的45.8%和41.8%。潍坊直排口、虞河和界河的DIN年排放量未超标，不需要削减其排放量。其中，-5 873.7%表示界河仍可继续容纳的DIN排放量是其排放现状的58.7倍，这个高值的出现主要是由于目前界河的DIN排放量相对较低造成的。

4.3非线性规划法合理性验证

由月最高浓度场求得的营养盐响应系数场包含了非线性部分，必须使用非线性规划法对其进行环境容量的计算。将各点源的DIN营养盐排放量依次设定为非线性规划法计算所得的最优排放量，可得到优化后的各水质控制点处DIN浓度的年变化曲线(图8)。

由图8可知，使用基于月最高浓度场的非线性规划法进行点源DIN的最优排放量规划时，可以保证水质控制点处的营养盐浓度不超标，从而使湾内生态系统对于突发事件具有足够的缓冲能力;而使用线性规划法规划点源的DIN排放量时，在某些时刻水质控制点处将连续出现DIN浓度超标的情况，严重危害了莱州湾海域生态系统的健康[12]。

由此也说明，应该使用基于月最高浓度场的非线性规划法进行营养盐等非保守物质的环境容量计算；线性规划法并不适用于求解非保守物质的环境容量。

图8　莱州湾点源的DIN优化排放时各水质控制点处DIN浓度的年变化曲线Fig.8　The annual variation of DIN with the DIN discharge optimization in Laizhou Bay

5讨论

目前可以进行最优化计算的软件中，比较常见的包括Matlab, SAS, OriginPro, SPSS, GraphPad, TableCurve2D, TableCurve3D, DataFit等。而以上软件中最常用的算法是麦夸特法(Levenberg-Marquardt)或者简面体爬山法(Simplex Method)等，这两种算法均归属于局部最优解法，其计算结果只能保证是局部最优解，并不能完全保证其全局最优性。并且使用局部最优算法时必须给出合适的初始值，若初值选取不恰当，则计算难以收敛，这是局部最优解法难以克服的瓶颈。在今后的工作中应该换用全局最优算法进行最优化计算，以克服局部解法存在的缺陷。

由于N、P营养盐相互影响时计算量太大，而现有的计算能力有限，所以本文在保持莱州湾N(或P)的排放现状不变的基础上，采用非线性规划法计算P(或N)的环境容量和剩余环境容量。但是在现实情况下N、P营养盐是相互作用的，不可能保持某一种营养盐排放现状不变，故两者相互作用时的营养盐环境容量有待进一步的研究。

6结论

为获得确定的响应系数场，采用满足排海通量最优化法适用条件的营养盐月最高浓度场数据，依次计算出莱州湾海域每个月的营养盐响应系数场，并进一步使用非线性规划法进行营养盐月环境容量的计算。

N、P营养盐之间相互影响、此消彼长，当N(或P)的浓度处于不同状态时，相应的P(或N)的环境容量是不一样的。本文旨在探讨计算非保守物质环境容量的方法，为简化数值计算量，保持P(或N)的排放现状不变，计算相应的莱州湾N(或P)的环境容量和剩余环境容量。计算结果显示，在保持莱州湾DIP排放现状不变的前提下，其莱州湾DIN年环境容量为2 843 t/a，剩余环境容量为-8 264 t/a，需要减少点源DIN的排放量；而在保持莱州湾DIN排放现状不变的前提下，其莱州湾DIP年环境容量为1 794 t/a，剩余年环境容量达1 485 t/a，可继续容纳部分DIP的排放。

参考文献(References)：

[1]KROM M, HORNUNG H, COHEN Y. Determination of the environmental capacity of Haifa Bay with respect to the input of mercury[J]. Marine Pollution Bulletin, 1990, 21(7): 349-354

[2]MARGETA J, BARIC A, GACIC M. Environmental capacity of the Kastela Bay[J]. Eighth International Ocean Disposal Symposium, 1989, (5): 9-13

[3]NORDVARG L, HAKANSON L. Predicting the environmental response of fish farming in coastal areas of the Aland archipelago using management models for coastal water planning[J]. Aquaculture, 2002, 206(3): 217-243.

[4]TEDESCHI S. Assessment of the environmental capacity of enclosed coastal sea[J]. Marine Pollution Bulletin, 1991, 23(2): 449-455.

[5]WU J, WANG Z J. Study of the seawater exchange and self-purification of Dalian Bay[J]. Marine Sciences, 1983, 6(3):32-35.吴俊，王振基. 大连湾海水交换及自净能力的研究[J]. 海洋科学, 1983, 6(3): 32-35.

[6]ZHANG Y, XIE S Y,YANG F. Research and application of the method for water quality evaluation in the sea area off western Guangdong [J]. Advances in Marine Science, 2012, 30(2):198-204. 张莹，谢仕义，杨锋.粤西海域水质评价方法的研究及应用[J].海洋科学进展，2012，30(2)：198-204.

[7]CHEN C H, ZHAN X W. Assessment on present status of seawater quality in eastern sea area of Xiamen[J]. Journal of Oceanography In Taiwan Strait, 2001, 20(2):157-160.陈朝华，詹兴旺. 厦门东侧海域表层海水质量现状与评价[J]. 台湾海峡, 2001, 20(2):157-160.

[8]WANG X, AN Y, ZHANG J, et al. Contribution of biological processes to self-purification of water with respect to petroleum hydrocarbon associated with No.0 diesel in Changjiang estuary and Jiaozhou Bay[J]. Hydrobiologia, 2002, 469(1-3): 179-191.

[9]GUO L B, JIANG W S, LI F Q. Environmental capacity calculation of COD and PHs in the Bohai Sea[J]. Periodical of Ocean University of China: Natural Sciences, 2007, 37(2):310-316.郭良波，江文胜，李凤岐,等. 渤海COD与石油烃环境容量计算[J]. 中国海洋大学学报：自然科学版，2007, 37(2): 310-316.

[10]LI K Q, WANG X L, YAN J, et al. Calculation on environmental capacities of petroleum hydrocarbon in Jiaozhou Bay[J]. Marine Environmental Science, 2003, 22(4):13-17.李克强，王修林，阎菊,等. 胶州湾石油烃污染物环境容量计算[J]. 海洋环境科学，2003, 22(4): 13-17.

[11]WANG X L, DENG N N, LI K Q, et al. Petroleum pollution condition and estimation of it environmental capacities summer in Bohai Sea[J]. Marine Environmental Science, 2004, 23(4):14-18.王修林，邓宁宁，李克强,等. 渤海海域夏季石油烃污染状况及其环境容量估算[J]. 海洋环境科学，2004, 23(4): 14-18.

[12]YU L L. The environmental capacity study of nitrogen and phosphorus nutrients in Laizhou Bay based on a three dimensional ecosystem dynamic numerical model[D]. Qingdao: Ocean University of China, 2012:1-131.虞兰兰. 基于三维海洋生态动力学模式的莱州湾氮、磷营养盐环境容量研究[D]. 青岛：中国海洋大学,2012:1-131.

[13]YU Y,PENG C S,YU L L, et al.Environmental capacity calculation of petroleum hydrocarbons in Laizhou Bay—The discharge optimization method based on the nonlinear programming[J].Marine Environmental Science,2014,(2):293-299.虞洋，彭昌盛，虞兰兰，等.莱州湾石油烃海洋环境容量计算——基于非线性规划的排涨通量最优化法[J].海洋环境科学，2014，(2)：293-299.

[14]GAO H W, WU D X, BAI J, et al. Distributions of environmental parameters in Laizhou Bay in summer, 2000[J]. Journal of Ocean University of Qingdao, 2003, 33(2):185-191.高会旺，吴德星，白洁，等. 2000年夏季莱州湾生态环境要素的分布特征[J]. 青岛海洋大学学报, 2003, 33(2): 185-191.

[15]WAN X Q, WU D X, BAO X W, et al. Distribution features of main observation elements in Laizhou Bay in the summertime of 2000[J]. Periodical of Ocean University of China: Natural Sciences, 2004,34(1):7-13.万修全，吴德星，鲍献文，等. 2000年夏季莱州湾主要观测要素的分布特征[J]. 中国海洋大学学报：自然科学版, 2004,34(1):7-13.

Received: May 22,2015

*收稿日期：2015-05-22

作者简介：虞兰兰(1985-)，女，江苏盐城人，工程师，博士，主要从事物理海洋学资料分析方面研究.E-mail ：229857443@qq.com *通讯作者：杨波(1974-)，男，山东莱州人，副教授，博士，主要从事浅海动力学方面研究.E-mail：yang.bo@ouc.edu.cn(王佳实编辑)

中图分类号：X834

文献标识码：A

文章编号：1671-6647(2016)02-0304-09

doi:10.3969/j.issn.1671-6647.2016.02.015

Environmental Capacity Calculation of Nutrients in the Laizhou Bay Based on the Nonlinear Programming

YU Lan-lan1, CHEN Yan-li1, ZHANG Hai-ning1, YANG Bo2, JIANG Wen-sheng2, YU Yang3

(1.NorthChinaSeaData&InformationCenter,SOA, Qingdao 266061, China;2.PhysicalOceanographyLaboratory,OceanUniversityofChina, Qingdao 266100, China;3.NorthSeaFleetMarineHydrometeorologicalCenter, Qingdao 266003, China)

Abstract:In order to calculate the nutrients environmental capacity (EC) of nonconservative substances in the sea, the nonlinear programming method was proposed based on the previous studies. The Laizhou Bay with stronger hydrodynamic force, larger and semi-closed sea area was selected for simulation and calculation. In this study, the data of highest monthly concentration field of nonconservative substances such as nitrogen and phosphorus (DIN and DIP) were selected to calculate the monthly nutrient response coefficient field in Laizhou sea area, and the monthly nutrient EC was obtained by using nonlinear programming. The results showed that the nutrient EC in Laizhou Bay in early 21st century could be as below. In the context of DIP discharge unchanged, the DIN EC in the bay is 2 1843 t/a and DIN SEC is -8 264 t/a, suggesting the DIN discharge of source point should be reduced. In the condition of DIN discharge unchanged, the DIP EC in the bay is 1 794 t/a and DIP SEC is 1 485 t/a, suggesting the discharge of DIP could be partly accommodated.

Key words:Nutrients; Laizhou Bay; nonconservative substance; environmental capacity; nonlinear programming

资助项目：2015年北海分局海洋科技项目——海洋站海冰观测资料的质控方法研究与初步应用(2015B14)