韩 鹏 秦京伟

（中国寰球工程有限公司 北京 100012）

1 研究背景

北京大学未名湖区是北京市第四批文物保护单位，也是北大学子们学习之余休闲、散心之处。未名湖作为整个湖区的中心，在北大人心中占有重要地位。

北京大学的地表水系主要由18个大小湖泊组成[1]。其中未名湖是面积最大的水域。

20世纪60年代以前，未名湖是活水，但是近几十年由于北京市地下水的过量开采，无论是地表水还是地下水都几近枯竭，没有自然水来补充，未名湖成为死水。这样一来，污染物得不到有效清除，未名湖的水质状况逐渐恶化，富营养化日趋严重，总磷、总氮和BOD严重超标，影响了周围师生的身心健康和湖区的观赏价值。

2 研究目的

建立未名湖BOD/DO模型，为未名湖BOD的治理提供支持。

3 研究方法和步骤

3.1 模型构建

BOD为生化需氧量的缩写，是表示水中有机物等需氧污染物质含量的一个综合指标。

DO是指溶解在水里氧的量，用每升水里氧气的毫克数表示。常温差压下，纯水里大约溶解氧8.66mg/L。当水中的溶解氧值降到5mg/L时，一些鱼类的呼吸就发生困难。

当水体受到有机物污染时（即富营养化），耗氧严重，溶解氧得不到及时补充，水体中的厌氧菌就会很快繁殖，有机物因腐败而使水体变黑、发臭。因此有必要建立BOD/DO模型对水中的BOD、DO进行模拟，寻找解决富营养化的途径或检验湖泊治理措施的好坏。

未名湖面积很小，只有36018平方米[1]，因此需要建立一个集总（零维）的水质模型。

根据质量守恒原理，某时段任何污染物含量的变化等于该时段流入总量减去流出总量，再减去元素降解或沉淀等所损失的量。其统一表达为：

其中M表示湖内污染物总量，minmout表示污染物流入、流出湖泊的速率，k表示污染物衰减率，t表示时间。

对于未名湖的BOD来说，假设湖泊水总体积不变，入湖流量等于出湖流量，可得：

其中V表示湖的体积，c表示BOD浓度，Q表示入湖流量，W表示非点源途径入湖速率，kr表示BOD削减速度：

kd表示因好氧有机物氧化分解而削减的速率，ks表示因沉降作用而削减的速率。

对于未名湖来说，由于水并不是常流动的，所以在计算一般状态时可以取Q值为零，在未名湖补水时，Q才有正的值。

对于氧气来说，有：

其中D表示氧负，即氧饱和浓度与当前浓度之差，D=os-o，ka表示复氧速率。其余变量与（3）式相同。

（2）、（4）两式即为未名湖BOD/DO模型的核心，可以看到，模型本身是十分简单的，关键是各个参数的收集和校正。模型的校正可使用一次补水后的BOD、DO浓度随时间变化的曲线，以及补水后经过足够长的时间后BOD、DO的平均值。

3.2 数据收集

（1）湖泊体积V

未名湖的面积在2002年左右为36018平方米[1]，鉴于未名湖的水位变化范围很小，故认为现在的面积也为36018平方米。未名湖平均水深可以使用船配合测深杆、测深锤，在湖面上等距测深。取25米间距的格点，共45个。如图2中的黑点所示。

图2 未名湖等距（25米）网格（底图：Google Earth）

（2）入湖流量

入湖流量数据需要向学校管理部门咨询，包括补水时间、频率、来源、每次持续时间以及水流量等。

（3）入湖BOD浓度

补水来源分为两种情况，一是使用万泉河水补给，二是使用自来水补给（包括直接使用抽取的地下水补给）。目前尚不清楚补水来源，所以BOD浓度的测试需要等到补水来源信息获取后进行。

BOD的测定有多种方法，如标准稀释法、微生物电极法、活性污泥曝气降解法、测压法等。根据我国颁布的水质分析方法GB7488-87，采用便准稀释法，即将水样稀释至一定浓度后，在20摄氏度恒温下培养5天，测量培养后水中的溶解氧量，便可计算出5日生化需氧量（BOD5），也可以根据BOD5计算出总BOD的值。[2]这种方法也是国际上约定俗成的分析方法。

根据BOD5计算出总BOD的方法见下：

（4）非点源BOD入湖速率

非点源BOD的测定十分困难，这里为了计算简便，做一定简化：假设非点源BOD入湖速率是常量，不随时间变化。首先根据文献，估计出非点源BOD入湖速率的大概值以及可能的范围，在模型校正时再确定此项的敏感性，做出相应调整。

（5）BOD氧化分解速率

Hydroscience在1971年发表的文章中推荐使用下面的经验公式来计算20摄氏度下BOD氧化分解速率（d-1）：[3]

其中H表示水体的深度（ft）。

该式表明深度在2.4米（8英尺）以下，BOD氧化分解速率随着深度增加而减少，在2.4米以上BOD氧化分解速率为恒定值。

当温度为其他温度时，使用下式计算[3]：

其中kd20为20摄氏度的速率，T表示温度，θ≌1.047。

（6）BOD沉降分解速率

未经处理的污水中，悬浮固体的含量很高，故沉降作用很明显，而且水越浅，沉降越快。

BOD沉降分解速率可由下式计算[3]：

其中vs表示BOD的沉降速率，H表示水深。

vs的测定可以根据文献估计，然后在模型校正时进行调整。

（7）复氧速率

停滞水体的复氧速率与风速有很大关系，对于湖体，有很多经验公式，这里选用Banks在1975年提出的一个经验公式（d-1）：[3]

其中Uw表示水面上空10米的风速（m/s），H表示水深（m）。

当温度为其他温度时，使用下式计算[3]：

其中ka20为20摄氏度的速率，T表示温度，θ≌1.024。

风速的确定可以使用北京市平均风速，也可以自己实地测定。不过由于风速的可变性，使用平均值可能使结果不准确，所以在模型校正时要对ka进行校正。

（8）饱和溶解氧

影响饱和溶解氧浓度（os）的三大因素是温度、盐度、气压。未名湖海拔约为50米，盐度尚未测定。可以在测量未名湖深度的同时取样，使用盐度计测定未名湖水的平均盐度。

标准大气压、盐度为零的水的饱和溶解氧计算公式见下（mg/L）：[3]

其中Ta表示绝对温度（K）：

当标准大气压下，盐度不为零时，首先计算，然后使用下式计算：[3]

其中 S表示盐度（ppt）。

当海拔不为零时，使用下式计算：[3]

os1为1个大气压下的饱和溶解氧。

整个模型的关键参数以及取得方法见表1：

表1 模型参数及取得方法

3.3 模型计算及校正

将（2）、（4）两式进行数值求解，为保证计算精度，同时考虑实际需要，时间步长选用1天，数值方法选用4级RK方法，由于BOD的方程独立于DO的方程，所以首先计算BOD，然后计算DO。

由于不知道未名湖是否有BOD/DO的监测数据，故要根据调查结果判断。如果官方有检测数据，那么可以使用已有数据；如果没有监测数据，那么研究就要花费一段时间（半年以上）进行BOD/DO监测。

当拿到监测数据后，就可以根据监测数据进行模型校正。计划使用遗传算法（Genetic Algorithm，GA），关于目标函数的选取将在拿到所有数据后再进行。

3.4 不确定性分析

计划使用MCMC（MarkovchainMonteCarlo）方法进行模型的不确定性分析，以确定各个参数的不确定性对模型输出结果的影响。

结语

模型搭建后，可以分析模型在管理上的作用。首先，要确定管理目标。管理目标可以分为近期目标和远期目标，比如在2年内达到国家IV类水标准，5年内达到国家III类水标准。一般情况下，管理部门要制定不同难度的多个目标。

其次，根据管理目标确定治理措施。治理措施分为工程性措施和非工程性措施。这些措施可以被参数化写入模型，根据模型来确定各种措施的治理效果，以及为达到目标所需措施的强度或种类。比如用污水处理厂的中水来代替万泉河的补水，在模型中体现为Q和cin的变化和延续时间；又如使用各种生物方式来治理湖泊，在模型中体现为W和kd的变化；再比如清淤措施可以表现为ks的变化等等。

可以根据刚刚搭建的模型来建立一个优化模型，使在一定置信水平下达到管理目标，同时成本最小。然后根据成本效益分析来确定管理目标和相应的治理手段。