长江水质评价和预测的数学模型.doc

高教社杯全国大学生数学建模竞赛承 诺 书我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B/C/D中选择一项填写）： B 我们的参赛报名号为（如果赛区设置报名号的话）： 所属学校（请填写完整的全名）： 参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)： 日期： 年 月 日赛区评阅编号（由赛区组委会评阅前进行编号）：高教社杯全国大学生数学建模竞赛编 号 专 用 页赛区评阅编号（由赛区组委会评阅前进行编号）：赛区评阅记录（可供赛区评阅时使用）：评阅人评分备注全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：长江水质评价和预测的数学模型摘要水质评价和预测是实施水污染控制规划的重要基础。本文结合模糊数学和灰色系统预测理论，构建长江水质评价和预测的数学模型并利用数学软件对模型进行了求解，解决了如何综合评价长江水质以及如何预测水质的污染情况问题。(1)问题一,结合水质评价的特点, 采用了直接确定多指标权重和有限样本隶属度的模糊分析法，对长江近两年多的水质情况做出了定量的综合评价。用模糊综合评价法对各断面水质的综合评价结果如表：四川攀枝花重庆朱沱湖北宜昌湖南岳阳城陵矶江西九江安徽安庆江苏南京四川乐山四川宜宾四川泸州湖北丹江口湖南长沙湖南岳阳楼湖北武汉江西南昌江西九江江苏扬州(2)问题二，把长江干流的7个观察站分成6个类似的独立系统来研究，通过分析水流流入和流出各个系统的过程，建立江水污染物浓度随时间变化含参变量的微分方程模型，在江水污染浓度恒定和自然净化效率成线性关系情况下，精确算出刚进入各个观察站的上游污水的浓度，利用关系式观察值=上游污水浓度+本地污染造成的浓度，代入数据得出结论：CODMn的主要污染源是湖南岳阳城陵矶、湖北宜昌和四川攀枝花；NH3-N的主要污染源是重庆朱沱和湖南岳阳城陵矶。(3)问题三，根据灰色系统预测法来对未来10年长江水质污染的发展趋势做出预测；首先建立GM（1，1）模型，再用灰色系统理论应用软件对数据进行处理得到未来10年的长江水质报告表，然后用1996年到2004年的模拟值、残值和相对误差对报告表进行检验，经检验可知一次预测10年的水质情况是存在较大误差，最后我们用一次预测两年对水质报告表进行修正得到比较合理的长江水质报告表。(4)问题四，根据问题三中得出的类和类的百分比之和，当两者之和小于20%时，建立约束条件和目标函数进行目标规划，得出未来10年内每年最大排污量，从附件四中10内的数据中建立浓度与时间的回归模型，得出其相关函数并预测出未来10内每年的排污量及每年需要处理的污水量。 最后我们根据实际对模型做了评价和适当改进，使之更符合实际运用。关键词 水质评价和预测 模糊综合评价法 微分方程 灰色系统预测法目标规划一 问题重述长江水质的污染程度日趋严重，已引起了相关政府部门和专家们的高度重视。附件3给出了长江沿线17个观测站（地区）近两年多主要水质指标的检测数据，以及干流上个观测站近一年多的基本数据（站点距离、水流量和水流速）。通常认为一个观测站（地区）的水质污染主要来自于本地区的排污和上游的污水。一般说来，江河自身对污染物都有一定的自然净化能力，即污染物在水环境中通过物理降解、化学降解和生物降解等使水中污染物的浓度降低。反映江河自然净化能力的指标称为降解系数。事实上，长江干流的自然净化能力可以认为是近似均匀的，根据检测可知，主要污染物高锰酸盐指数和氨氮的降解系数通常介于0.10.5之间，比如可以考虑取0.2(单位：1/天)。附件4是“19952004年长江流域水质报告”给出的主要统计数据。下面的附表是国标(GB3838-2002) 给出的地表水环境质量标准中4个主要项目标准限值，其中、类为可饮用水。要求：（1）对长江近两年多的水质情况做出定量的综合评价，并分析各地区水质的污染状况。（2）研究、分析长江干流近一年多主要污染物高锰酸盐指数和氨氮的污染源主要在哪些地区?（3）假如不采取更有效的治理措施，依照过去10年的主要统计数据，对长江未来水质污染的发展趋势做出预测分析，比如研究未来10年的情况。（4）根据预测分析，如果未来10年内每年都要求长江干流的类和类水的比例控制在20%以内，且没有劣类水,那么每年需要处理多少污水？ （5）对解决长江水质污染问题提出切实可行的建议和意见。二 问题分析问题一 由于水体污染指标之间的相容性(不存在传递性、绝对化) , 且水质评价本身是一种实践性、时空性、技术性均很强的多属性、多指标决策,同时评价过程中所涉及的各指标权重。本文针对水质评价的特点、要求以及求解过程中的难点, 采用了直接确定多指标权重和有限样本隶属度的模糊分析法，对问题1进行了模糊综合分析法。根据国家水质标准（表2）来综合评价。问题二 中把长江干流的七个观察站分成六个类似的独立系统来研究，因各个系统中某种污染物的浓度都是时刻在变化的，通过分析水流流入和流出该系统的过程，建立江水污染物浓度随时间变化含参变量的微分方程模型。通过模型算出刚进入各个观察站的上游污水的浓度，利用（观察值=上游污水浓度+本地污染造成的浓度）的关系，算出各个观察站所在地区在某特定时间对长江造成的污染浓度。问题三 首先对附件4的数据进行分析然后整理成为9个表，然后建立灰色系统预测法来对未来十年长江水质污染的发展趋势做出预测，再用灰色系统理论应用软件对数据进行处理得到未来十年的长江水质报告表，最后用1996年到2004年的模拟值、残值和相对误差对报告表进行检验，由此可知一次预测十年的水质情况是不准确的，我们用一次预测两年对水质报告表进行修正得到比较合理的长江水质报告表。问题四是根据问题三中得出的类和类的百分之和，当两者之和小于20%时，建立约束条件和目标函数，得出未来10年内每年最大排污量，从附件四中10内的数据中建立浓度与时间的回归模型，得出其相关函数并预测出未来10内每年的排污量，减去每年至多的排污量，则为每年需要处理的污水量。 三 各问题假设、参数说明、模型建立与求解问题一的建模与求解针对水质评价的特点、要求以及求解过程中的难点, 采用了直接确定多指标权重和有限样本隶属度的模糊分析法，对问题1进行了模糊综合分析。根据国家水质标准（表2）来综合评价。 符号说明 因子的实测浓度值均值； 因子各级水质标准的算术均值； 因子的权重值； 因子对j级水质的隶属度； 因子的实测浓度值； 因子第j级水质标准；1、 计算出17个观察站在28个月中的各个参数的平均值（表1），并建立评价对象的因素，根据国家水质标准（表2）来综合评价因素是参与评价的评价指标，即问题中的实际测定浓度的模糊子集。根据我国地表水环境质量标准（GB3838-2002），并考虑长江水污染的特点，选择四个参数作为评价因子=（PH，DO，CODMn,NH3-N）。表 1 17个观察站在28个月中时间上的平均值PHDOCODMnNH3-N四川攀枝花龙洞8.2560719.1542862.4321430.182857重庆朱沱 7.8935718.8960712.3071430.334643湖北宜昌南津关 7.7507148.5053572.8750.264286湖南岳阳城陵矶 7.8167868.6832143.7857140.330000江西九江河西水厂 7.4242867.7535712.4285710.160357安徽安庆皖河口 7.0846437.0942862.5107140.911429江苏南京林山 7.6500007.4910712.0928570.127857四川乐山岷江大桥 7.4957145.5585715.2428570.924286四川宜宾凉姜沟 8.0757148.9760712.7357140.430357四川泸州沱江二桥 7.6796436.9325003.2964290.815714湖北丹江口胡家岭 7.8771439.2910711.9535710.092143湖南长沙新港 7.0846437.0942862.5107140.911429湖南岳阳岳阳楼 7.7285718.3150004.1928570.385714湖北武汉宗关 7.9475007.4214293.3250.197500江西南昌滁槎 7.1103575.6982142.3239294.633214江西九江蛤蟆石 7.6450007.8946433.7321430.287143江苏扬州三江营 7.6817868.1378573.0214290.287143地表水环境质量标准（GB38382002）中4个主要项目标准限值 单位：mg/L表 2序号  分 类            标准值  项 目类类类类类劣类1溶解氧(DO)          7.5（或饱和率90%）6532 02高锰酸盐指数(CODMn) 24610153氨氮（NH3-N）      0.150.51.01.52.04PH值（无量纲）6-92、 立评价集 V是与U评价因子相应的评价评价标准集合。在水质评价中，V是各个污染因子相应的环境质量标准等级的集合。由于水体污染程度是一个模糊概念，就以问题中的分类为评价等级，把在表水分为6个等级。即评价集为：V，劣3、 定评价因素的模糊权向量通常各因子的重要程度不同，因此对每个因子赋于一个相应的权重（i=1,2,3,4）构成权重集A。的确定方法如下： （1）其中是因子的实测浓度值，是因子各级水质标准的算术均值，是因子的权重值。上式中四项因子指标，分别计算出权重后，组成一个模糊权重集，且。根据以上步骤，可计算出各单项参数权重值，得到模糊权重集A。表 3各项参数各断面点名称PHDOCODMnNH3-N四川攀枝花龙洞0.310.550.10.05重庆朱沱0.30.530.090.09湖北宜昌南津关0.30.520.110.07湖南岳阳城陵矶0.280.50.140.86江西九江河西水厂0.310.530.10.05安徽安庆皖河口0.260.410.090.24江苏南京林山0.340.530.090.04四川乐山岷江大桥0.330.40.230.03四川宜宾凉姜沟0.30.460.160.03四川泸州沱江二桥0.280.460.10.11湖北丹江口胡家岭0.270.440.120.21湖南长沙新港0.310.50.080.27湖南岳阳岳阳楼0.250.410.090.23湖北武汉宗关0.320.470.130.05江西南昌滁槎0.140.220.040.66江西九江蛤蟆石0.30.480.140.08江苏扬州三江营0.30.50.150.084、 隶属度的确定从一个出发进行评价，以确定评价对象对评价集元素的隶属程度亦称为单因素模糊评价。属于第j级水质的隶属函数为： （2）其中是因子对j级水质的隶属度，是因子的实测浓度值，是因子第j级水质标准。对第I个因子评价的结果组成单因素模糊评价集。根据以上步骤，可计算出监测断面的模糊关系矩阵R。得单因素评价模糊矩阵 表4给出四川攀枝花龙洞等断面各参数的隶属度。表 4等级PH值DOCODMnNH3-N四川攀枝花0.1610.7850.91430.8400.2150.0857000000000000劣0000重庆朱沱0.393310.950.48570.606700.050.5143000000000000劣0000湖北宜昌南津关0.510.560.68570.500.440.3143000000000000劣0000湖南岳阳城陵矶0.453310.1050.48570.546700.8950.5143000000000000劣0000江西九江河西水厂0.7210.7850.97140.2800.2150.0286000000000000劣0000安徽安庆皖河口0.70670.97330.710.77140.29330.02670.290.2286000000000000劣0000江苏南京林山0.56670.99330.95510.43330.00670.0450000000000000劣0000四川乐山岷江大桥0.66670000.33330.560.620.1600.440.380.8400000000劣0000四川宜宾凉姜沟0.2810.630.20.7200.370.8000000000000劣0000四川泸州沱江二桥0.54670.53330.3300.45330.46670.670.380000.6200000000劣0000湖北丹江口胡家岭0.41331110.5867000000000000000劣0000湖南长沙新港0.94670.740.75500.05330.260.2450.160000.8400000000劣0000湖南岳阳岳阳楼0.5133100.31430.486700.9050.6857000.095000000000劣0000湖北武汉宗关0.36670.94670.3350.85710.63330.05330.6650.1429000000000000劣0000江西南昌滁槎0.926700.8400.07330.70.16000.30000000000劣0001江西九江蛤蟆石0.586710.130.60.413300.870.4000000000000劣0000江苏扬州三江营0.546710.490.60.453300.510.4000000000000劣00005、 综合评价根据模糊评价的原理，考虑所有因子的影响，将模糊权向量A与单因素模糊评价矩阵R复合，得到四川攀枝花龙洞断面的综合评价向量B。因子的模糊综合评价向量B=（0.5371,0.4628,0,0,0,0）即可得到评价指标，因此，四川攀枝花龙洞断面评价结果为：对级，级，级，级，级，级的隶属度依次是0.5371,0.4628,0,0,0,0,可见对级水的隶属度最大，根据模糊综合评价的最大隶属原则，可知四川攀枝花断面为级水。同理可得其他断面综合评价向量B，结果见下表。表 5 长江各17个监测断面隶属度水质级别计算表级别项各断面名称劣级别四川攀枝花龙洞0.53710.46280000重庆朱沱 0.2504062300湖北宜昌南津关 0.34000.6601000湖南岳阳城陵矶 0.49410.505800江西九江河西水厂 0.16570.83430000安徽安庆皖河口 0.20740.79250000江苏南京林山 0.21530.78470000四川乐山岷江大桥 0.10390.18240.15620.346200需加权四川宜宾凉姜沟 0.31660.68340000四川泸州沱江二桥 0.17050.35790.37140.261400需加权湖北丹江口胡家岭 0.61200.38800000湖南长沙新港 0.16040.26570.5738000湖南岳阳岳阳楼 0.02200.47360.5043000湖北武汉宗关 0.41690.53920000江西南昌滁槎 0.10470.63460.1374000.0934江西九江蛤蟆石 0.407500.5925000江苏扬州三江营 0.34010.65990000根据模糊综合评价的最大隶属度原则，由上表可知四川乐山岷江大桥、四川泸州沱江二桥的水质分别为、级，但由于其他等级隶属度原则。本论文加权均解决了解这一矛盾，加权平均原则的思想是：将事物的等级看作一种相对位置，使其连续化，为了能定量处理，用“1，2，3，m”依次表示各等级，并称其为各等级的秩。然后用B中对应分量将各等级的秩加权求和，得到被评事物的相对位置。其计算公式为：。由上可知式中隶属于第j级水的隶属度，K为待定系数（1或2），目的是控制较大的所起的作用，由上表可知四川乐山岷江大桥和四川泸州二桥断面需用加权平均法进行最后的综合评判。计算得到四川乐山岷江大桥和四川泸州二桥断面的分别为1.23和3.21,所以四川乐山岷江大桥水质级别为：级，略偏向于级；四川泸州二桥断面水质级别为：级，稍偏向于级。由上表可得，长江水各监测断面水质基本都有在国家GB3838-2002V级地面水标准，只有个别断面水质较差，可见长江地面水污染相当严重，这些污染主要来自纺织、化工轻、轻工、造纸等行业及城镇居民生活污水的直接排放。结合模糊评价的最终结果，可知湖南岳阳楼断面水质污染最严重，处于湖北丹江口胡家岭 的水质污染最轻。具体原因是前者处于工业区和人口密集的地方，后者则是处在工业较少和人口密度较小的地方，评价结果与实际情况吻合较好。从28个月中的数据中随机抽取4个，做出如下图像，可以大概知道其趋势。图 2可见原模型的综合分析还是贴近实际的。问题二的建模与求解1、模型假设（1）每两个水文站观察点之间的江面看成一个单流入、单流出的系统，不考虑支流水对干流某种化学物质的影响。（2）某地点的污染物能以很快的速度与干流中的水均匀混合。（3）参与模型的变量是连续变化的，并且充分光滑。（4）干流中每两个水文站观察点之间水体的体积在一定的时间保持定常，时间间隔为一个月。（5）考虑物理降解、化学降解、生物降解在水体自净过程中的作用，综合的水体自净系数为。（6）每个水文站观察点的所在的本地污染源与观察点无限的接近，近似看成一个点。（7）对于每个研究的江面在每个特定的时间内水流的流速不变，包括两个端点。（8）每两个观察站之间没有污染源。（9）系统的体积保持不变，假设有降水等原因引起水的流入量与水的增发、渗透所造成的失量相抵消。2、问题分析 对于本问题，把它细分到在某个特定的时间来进行研究，时间间隔为一个月。首先我们把长江干流的七个观察站分成六个类似的独立的系统来研究。相对于每个系统，系统内污染物的浓度只要是受该系统污染源的污染物输入、系统本身的自净系数共同影响的。水流的动态流程图如下： 系统 图 3由上图可以清楚地知道，系统里面污染物量的变化由上一个观察站的观察值（实际为上游进入观察点污染物的量和该观察点地区的污染物输入的综合）、系统排除的浓度与系统的自净系数共同决定，因此我们可以列出系统污染物变化的微分方程。其次，根据微分方程模型，算出进入后面六个观察点的污染物的浓度。再利用观察值、上游污水浓度、本地污染程度之间的关系，就可以得出本地的污染浓度。3、参数说明 时刻流出系统的污染物的浓度； t流入系统的浓度（水文站的观察值）； 时刻流入系统的速度； 时刻流出系统的速度； 时刻某截断干流中污染物的浓度； 水流的流速； 某系统的水流的体积； 系统的自净系数； 生物降解的速率； 时刻系统的体积； 该地区排进系统的排污浓度。4、模型建立 由假设3，参与模型的变量都是连续而且充分光滑的。由假设9可知=（常数）。在不考虑系统自净能力的情况下可以得到方程： （1） 由于为常数，故有=，流出的污染物应与系统中的污染物有相同的浓度=。根据假设7，进一步我们可以得出=。这样，我们可以得到： （2）在上面的基础上，根据本问题的特点，我们考虑流进系统的自净能力。系统的污染物流入量仍然为，而流出系统的污染物流出量则变成自净速率和流出量两者的和。所以我们可以得到下面方程： （3）由自净系数与自净速率的关系可以得出： （4）由（3）、（4）可得： （5）分离变量得：两边积分：（为常数，同是）即：两边取对数有：化简得：当时， 代入方程得即方程可以变为 （6）由本地区排污能够与江水瞬间混合并且不考虑排入的水量，我们可以得出本地区的排污量的表达式： （7） 由方程（6）我们可以算出刚好进入后面六个观察站某种污染物的浓度。再代入方程（7）就可以解得后面六个观察站本地某种污染物污染的浓度。5、模型求解根据题目的附件（3），在某个特定的研究时间内，这里是一个月 。用EXCEL可以算出水流从一个观察站流到相邻观察站的时间（也就是从系统的头端到系统的尾端的时间）、系统的体积、水流的速度。根据方程（6）和附件5的数据，用MATLAB编程可以得出刚好进入后面六个观察站水流中某种污染物的含量。如下列各表为(单位：mg/L)：表 6刚进入重庆朱沱污染物浓度刚进入湖北宜昌污染物浓度CODMnNH3-NCODMnNH3-N2004.040.001740.00011332004.040.000357592.14554E-052004.050.003225.247E-052004.050.0002061542.53007E-052004.060.001862.972E-052004.060.0002978731.78724E-052004.070.001632.723E-052004.070.0003209081.75041E-052004.080.004470.00077112004.080.0002660192.26116E-052004.090.004025.935E-052004.090.0003699461.34526E-052004.10.000545.446E-052004.10.0001575132.26425E-052004.110.003447.382E-052004.110.0003004924.05664E-052004.120.003730.00018652004.120.0001733096.79904E-052005.010.002190.00012782005.010.0002823650.000117292005.020.001780.0002972005.020.0002315697.07571E-052005.030.001990.00046922005.030.0001559034.43092E-052005.040.002230.00020252005.040.0003404040.0001072272005.050.032850.00246442005.050.0034600430.000588979表 7刚进入湖南岳阳污染物浓度刚进入江西九江污染物浓度CODMnNH3-NCODMnNH3-N2004.048.3388E-054.52E-062004.046.52803E-055.1433E-062004.057.69682E-057.06E-062004.058.89314E-057.905E-062004.060.0001054321.29E-052004.060.0001332951.3329E-052004.070.0001526661.53E-052004.070.0001821721.5615E-052004.080.0001797081.61E-052004.080.0001471151.2137E-052004.099.25828E-054.63E-062004.090.0001031447.9342E-062004.10.0002008791.23E-052004.10.0001787381.8385E-052004.110.0001782141.52E-052004.119.1611E-051.2685E-052004.120.0002620521.63E-052004.120.0001383051.225E-052005.010.0001556411.97E-052005.010.0001655421.3633E-052005.020.0001394939.3E-062005.020.0001066034.3308E-062005.039.21253E-051.07E-052005.030.0001548221.586E-052005.040.0001036769.33E-062005.047.79708E-051.1292E-052005.050.0018228250.0001532005.050.001633530.0001505表 8 刚进入安徽安庆污染物浓度刚进入江苏南京污染物浓度CODMnNH3-NCODMnNH3-N2004.048.58724E-057.995E-062004.046.9032E-056.42711E-062004.059.67704E-050.00011682004.056.60916E-055.16341E-062004.060.0001167976.7383E-062004.065.33885E-055.67253E-062004.070.0001276127.853E-062004.076.21044E-058.03704E-062004.089.49524E-059.9081E-062004.087.16886E-056.03694E-062004.095.01385E-056.5972E-062004.096.28626E-065.50048E-062004.10.000109121.0912E-052004.18.5369E-056.91082E-062004.110.0001106421.1064E-052004.116.80038E-053.7093E-062004.120.000110983