2005年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎,卢亚军。许超指导教师汪晓银 长江水质评价及预测模型的建立与分析 摘要: 水资源污染日趋严重,本文仅就长江流域若干观测点的最近10年的观测数据进行 建模分析,探讨水资源污染的预测与控制的问题。 对于问题一、我们首先采用模糊综合评判法进行评价,我们以溶解氧、髙锰酸盐指 数、氨氮组成的评价参数集为U={u1,u2,u3},水质等级集为{vl,v2,v3,v4,v5,v6},用隶 属度来刻画长江沿岸17个测点水质在三个参数方面属于不同水质级别的情况,由观测 值和隶属度确定模糊评判矩阵,由最大隶属原则,求出28个时点平均值得到17个测点 的2年多的平均水质情况,其中江西南昌滁槎污染最为严重(测点15),其次是四川乐 山岷江大桥(测点8)。 我们还采用了另外一种一般的综合指数评价方法对水质情况作综合评价。选取Ⅲ等 水的4项指标值PH=7,DO=5,COD=6,NH=1为评价标准。先计算分指数K,然后再把 分指数加权(n)求和得到综合指标爬∑nK。因为17个测点分别有28个时点的监测 数据,我们把这28个月按年分为枯水期、丰水期和平水期。分类标准参照题中附件4 的说明(3)。然后用平均的方法求出按年按期分类的综合指标矩阵。用 Excel分别画出 2003、2004、2005年各测点的水质综合评价柱状图(见附件中图、图、图)。其中 Ⅲ等水的量化指标为0.1,值越低水质越好。最后得出:江西南昌滁槎污染最为严重(测 点15),其次是四川乐山岷江大桥(测点8)。可以看出两种方法评价的结果是一致的。 对于问题二、考虑到江水的流动性和自降解性,测点浓度不能与污染源等同起来, 所以我们把每个监测值看作是由上游对它的影响和本地污染两部分组成的。所以我们考 虑水流速度、测点间距计算出降解次数,进而计算出一个测点对下一个测点的影响度, 通过减去来自上游的影响进而求出每个测点的本地污染度。这样就可以通过比较本地污 染度找出主要污染源了。具体计算出主要污染源有:江西南昌滁槎、四川泸州沱江二桥 四川乐山岷江大桥、湖南长沙新港、湖南岳阳城陵矶、重庆朱沱、湖南岳阳楼等地区 对于问题三、问题要求预测未来十年的水质污染趋势,我们用预测未来10年中六 种水的比例来评价水质污染趋势。具体用灰色系统方法进行预测。又由于每年的总评价 河长不一样,所以不能直接预测,但是比例是不会变的,我们把10年的六种水的比例 用同一个总评价河长(我们假定为10000)转化。再进行预测,得出未来10年的六种水 的河长,然后再转化为比例值。可以看出未来10年长江劣V类水越来多,污染越来越 严重。 对于问题四、我们把6种水重新分为3类:清洁水(I、Ⅱ和Ⅲ)、污染水(Ⅳ和V) 劣质水(V)。用灰色系统方法预测未来10年原6种水的比例,折算出新3类水的比例, 建立废水率与3类水比例的逐步回归模型,确定未来10年每一年废水率与各种水比例 之间的回归方程。然后根据回归方程得到在控制状态下(Ⅳ和Ⅴ类不超过20%、劣Ⅴ类 没有)的废水率,然后与不加控制的预测值进行比较,得到需要处理的废水量(亿吨)。 2005年2006年2007年2008年2009年2010年2011年2012年2013年2014年 79497.51189141.9166.11926231.3261.32948330.1 对于问题五、根据我们所做模型的分析和预测结果,以及目前我国治理长江水质污 染的现状给出:合理制定工业发展规划、污水处理设施的利用、控制排污总量、加大执 法力度等四点建议和意见 关键字:水质评价模糊综合评价综合指数评价灰色系统逐步回归
2005 年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎，卢亚军。许超 指导教师 汪晓银 1 长江水质评价及预测模型的建立与分析 摘要： 水资源污染日趋严重，本文仅就长江流域若干观测点的最近 10 年的观测数据进行 建模分析，探讨水资源污染的预测与控制的问题。 对于问题一、我们首先采用模糊综合评判法进行评价，我们以溶解氧、高锰酸盐指 数、氨氮组成的评价参数集为 U={u1,u2,u3},水质等级集为{v1,v2,v3,v4,v5,v6}，用隶 属度来刻画长江沿岸 17 个测点水质在三个参数方面属于不同水质级别的情况，由观测 值和隶属度确定模糊评判矩阵，由最大隶属原则，求出 28 个时点平均值得到 17 个测点 的 2 年多的平均水质情况，其中江西南昌滁槎污染最为严重（测点 15），其次是四川乐 山岷江大桥（测点 8）。 我们还采用了另外一种一般的综合指数评价方法对水质情况作综合评价。选取Ⅲ等 水的 4 项指标值 PH0=7，DO0=5，COD0=6，NH0=1 为评价标准。先计算分指数 K，然后再把 分指数加权（η）求和得到综合指标 W=∑ηKj。因为 17 个测点分别有 28 个时点的监测 数据，我们把这 28 个月按年分为枯水期、丰水期和平水期。分类标准参照题中附件 4 的说明（3）。然后用平均的方法求出按年按期分类的综合指标矩阵。用 Excel 分别画出 2003、2004、2005 年各测点的水质综合评价柱状图（见附件中图[1]、图[2]、图[3]）。其中 Ⅲ等水的量化指标为 0.1，值越低水质越好。最后得出：江西南昌滁槎污染最为严重（测 点 15），其次是四川乐山岷江大桥（测点 8）。可以看出两种方法评价的结果是一致的。 对于问题二、考虑到江水的流动性和自降解性，测点浓度不能与污染源等同起来， 所以我们把每个监测值看作是由上游对它的影响和本地污染两部分组成的。所以我们考 虑水流速度、测点间距计算出降解次数，进而计算出一个测点对下一个测点的影响度， 通过减去来自上游的影响进而求出每个测点的本地污染度。这样就可以通过比较本地污 染度找出主要污染源了。具体计算出主要污染源有：江西南昌滁槎、四川泸州沱江二桥 、 四川乐山岷江大桥、湖南长沙新港、湖南岳阳城陵矶、重庆朱沱、湖南岳阳楼等地区。 对于问题三、问题要求预测未来十年的水质污染趋势，我们用预测未来 10 年中六 种水的比例来评价水质污染趋势。具体用灰色系统方法进行预测。又由于每年的总评价 河长不一样，所以不能直接预测，但是比例是不会变的，我们把 10 年的六种水的比例 用同一个总评价河长（我们假定为 10000）转化。再进行预测，得出未来 10 年的六种水 的河长，然后再转化为比例值。可以看出未来 10 年长江劣Ⅴ类水越来多，污染越来越 严重。 对于问题四、我们把 6 种水重新分为 3 类：清洁水（I、Ⅱ和Ⅲ）、污染水（Ⅳ和Ⅴ）、 劣质水（Ⅴ）。用灰色系统方法预测未来 10 年原 6 种水的比例，折算出新 3 类水的比例， 建立废水率与 3 类水比例的逐步回归模型，确定未来 10 年每一年废水率与各种水比例 之间的回归方程。然后根据回归方程得到在控制状态下（Ⅳ和Ⅴ类不超过 20%、劣Ⅴ类 没有）的废水率，然后与不加控制的预测值进行比较，得到需要处理的废水量（亿吨）。 2005 年 2006 年 2007 年 2008 年 2009 年 2010 年 2011 年 2012 年 2013 年 2014 年 79.4 97.5 118.9 141.9 166.1 192.6 231.3 261.3 294.8 330.1 对于问题五、根据我们所做模型的分析和预测结果，以及目前我国治理长江水质污 染的现状给出：合理制定工业发展规划、污水处理设施的利用、控制排污总量、加大执 法力度等四点建议和意见。 关键字：水质评价 模糊综合评价 综合指数评价 灰色系统 逐步回归
2005年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎,卢亚军。许超指导教师汪晓银 正文 问题的提出 黄河脏了,长江可以救黄河,长江完了,拿什么救长江?”。我们的生命线一 长江正在倍受煎熬,保护长江、保护水资源就是保护我们自己。而如今长江污染日趋严 重。那么长江现在到底现况如何?污染源在哪里?如果再不采取有效措施,长江的命运 将会如何?最主要的污染物一一污水的排放对长江产生什么样的影响?怎么样规划、采 取怎么样的措施才能使长江在保持“生命力”的前提下达到环境与经济和谐发展?等等 这些都是我们亟待解决的问题 二、问题的分析 2.1对水质情况的定量综合评价 水质的评价是对水中各种物质的浓度进行测定,分析是否超过某个标准的过程。题 目中只给出了4种物质作为评价指标。通过查阅资料我们知道一般分析水质情况的方法 有指数评价法、模糊评价法、灰色评价法、物元分析法、人工神经网络评价法等诸多方 法叫。我们在这里先用模糊评价法作出评价,然后用比较简单常用的指数评价法进行验 22寻找主要污染源 污染源是指污染的源头,即污染物的排放点。但是题目只给出了各个测点的污染物 的浓度值,我们很容易理解一个测点的监测值高不一定说明该点就是主要污染源。因为 江水是流动的,一个测点的数据应该包括两个部分,即来自上游的污染物的影响和本地 的污染物的浓度。各测点本地污染物的浓度才是我们考察的目标。那么上游污染物对下 测点的影响度该怎么计算呢?这就要考虑两测点间距离、水流速度及污染物降解与上游 污染物影响度之间的关系 23水质污染趋势的预测 预测需要有历史数据支持,根据题目所给数据的特点:过去10年每一年中6类水 的比例,以及水流量总量及排污量等不确定因素。我们可以灰色系统方法进行预测。且 用一年中6类水的比例也可充分表明水质状况 24水质污染的控制 在2.3的预测中,我们知道6类水的比例反映水质污染的程度,且题目也是要求我 们通过对污水的处理达到对6类水的比例进行控制的目的,即6类水的比例是我们的控 制目标,污水量是我们的控制点。所以我们要先找到污水量与6类水比例之间的关系 我们采用多元回归对污水排放比例和6类水的比例进行逐步回归分析,这样我们可以得 到一个回归方程,它就是我们要找排污比例与6类水比例之间的关系式,通过它我们就 可以在给定控制标准下的污水排放上限,然后跟我们预测的不加控制的排污量比较就得 到了需要处理的污水量。 、模型的假设 (1)假设溶解氧(DO)浓度越髙水质越好,不考虑过含氧情况。 (2)假设各监测指标之间无相互作用。 (3)假设我们研究的长江是一条平直的河流。 (4)假设所给数据真实可靠。 (5)假设水质状况只与题目给我们的4个项目有关,不考虑其他项目 2
2005 年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎，卢亚军。许超 指导教师 汪晓银 2 正文： 一、问题的提出 “黄河脏了，长江可以救黄河，长江完了，拿什么救长江？”。我们的生命线—— 长江正在倍受煎熬，保护长江、保护水资源就是保护我们自己。而如今长江污染日趋严 重。那么长江现在到底现况如何？污染源在哪里？如果再不采取有效措施，长江的命运 将会如何？最主要的污染物——污水的排放对长江产生什么样的影响？怎么样规划、采 取怎么样的措施才能使长江在保持“生命力”的前提下达到环境与经济和谐发展？等等 这些都是我们亟待解决的问题。 二、问题的分析 2.1 对水质情况的定量综合评价 水质的评价是对水中各种物质的浓度进行测定，分析是否超过某个标准的过程。题 目中只给出了 4 种物质作为评价指标。通过查阅资料我们知道一般分析水质情况的方法 有指数评价法、模糊评价法、灰色评价法、物元分析法、人工神经网络评价法等诸多方 法[1]。我们在这里先用模糊评价法作出评价，然后用比较简单常用的指数评价法进行验 证。 2.2 寻找主要污染源 污染源是指污染的源头，即污染物的排放点。但是题目只给出了各个测点的污染物 的浓度值，我们很容易理解一个测点的监测值高不一定说明该点就是主要污染源。因为 江水是流动的，一个测点的数据应该包括两个部分，即来自上游的污染物的影响和本地 的污染物的浓度。各测点本地污染物的浓度才是我们考察的目标。那么上游污染物对下 测点的影响度该怎么计算呢？这就要考虑两测点间距离、水流速度及污染物降解与上游 污染物影响度之间的关系。 2.3 水质污染趋势的预测 预测需要有历史数据支持，根据题目所给数据的特点：过去 10 年每一年中 6 类水 的比例，以及水流量总量及排污量等不确定因素。我们可以灰色系统方法进行预测。且 用一年中 6 类水的比例也可充分表明水质状况。 2.4 水质污染的控制 在 2.3 的预测中，我们知道 6 类水的比例反映水质污染的程度，且题目也是要求我 们通过对污水的处理达到对 6 类水的比例进行控制的目的，即 6 类水的比例是我们的控 制目标，污水量是我们的控制点。所以我们要先找到污水量与 6 类水比例之间的关系， 我们采用多元回归对污水排放比例和 6 类水的比例进行逐步回归分析，这样我们可以得 到一个回归方程，它就是我们要找排污比例与 6 类水比例之间的关系式，通过它我们就 可以在给定控制标准下的污水排放上限，然后跟我们预测的不加控制的排污量比较就得 到了需要处理的污水量。 三、模型的假设 （1）假设溶解氧(DO)浓度越高水质越好，不考虑过含氧情况。 （2）假设各监测指标之间无相互作用。 （3）假设我们研究的长江是一条平直的河流。 （4）假设所给数据真实可靠。 （5）假设水质状况只与题目给我们的 4 个项目有关，不考虑其他项目
2005年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎,卢亚军。许超指导教师汪晓银 四、号的定义与说明 41模糊评判模型符号定义 符号 符号说明 第i个参数在第j级别上的标准值 Y 第i个参数在某一级别上的监测值的第k个 i参数监测值介于Ay到A,之间的个数 L 第i个参数污染物监测值的个数 Z 参数污染物监测值介于AP到A.之间的1,个 监测值的平均值 z对第j级水质的隶属度 计算得到的i参数的权重 参数而言,介于A.到AJ之间的监测值发生在 J水质下的概率 归一化处理后i参数的权重 参数发在j水质下的模糊概率 水域水体出现j级水质的模糊综合概率 42综合指数评价模型符号定义 符号 符号说明 W 水质综合指标 K 分指数(监测数据与评价标准之比) 评价指标的权重 全流域17个测点COD的监测值 全流域17个测点NH3的监测值 P 全流域17个测点PH的监测值 D 全流域17个测点DO的监测值 支流10个测点COD的监测值 支流10个测点NH3的监测值 干流7个测点COD的监测值 N 干流7个测点NH3的监测值 五、模型的建立与求解 5.1长江水质的综合评价 5.1.1模糊综合评判模型 根据水域情况的质量标准我们把水污染监测浓度看成是一个离散的随机变量,用概 率统计方法进行统计可以得到水域属于某个标准的概率,因为可以拟定不同的水域标 准,评价参数集为U={u1,u2,u3},水质分级集为{v1,V2,v3,V45,v6},其中u1,u2,u3分别表 示为溶解氧,高锰酸盐指数,氨氮(NH3-N),因为PH值对水域影响不大,所以对其不
2005 年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎，卢亚军。许超 指导教师 汪晓银 3 四、 号的定义与说明 4.1 模糊评判模型符号定义 符号 符号说明 Aij 第 i 个参数在第 j 级别上的标准值 Xi,k 第 i 个参数在某一级别上的监测值的第 k 个 li,j i 参数监测值介于 Ai,j-1到 Ai,j之间的个数 Li 第 i 个参数污染物监测值的个数 Zi,j i 参数污染物监测值介于 Ai,j-1到 Ai,j之间的 li,j个 监测值的平均值 ri,j,j zi,j对第 j 级水质的隶属度 wi 计算得到的 i 参数的权重 pi,j i 参数而言，介于 Ai,j-1到 Ai,j之间的监测值发生在 j 水质下的概率 ai 归一化处理后 i 参数的权重 qi,j i 参数发在 j 水质下的模糊概率 qj 水域水体出现 j 级水质的模糊综合概率 4.2 综合指数评价模型符号定义 符号 符号说明 W 水质综合指标 K 分指数（监测数据与评价标准之比） η 评价指标的权重 Cj 全流域 17 个测点 COD 的监测值 Nj 全流域 17 个测点 NH3 的监测值 Pj 全流域 17 个测点 PH 的监测值 Dj 全流域 17 个测点 DO 的监测值 Czj 支流 10 个测点 COD 的监测值 Nzj 支流 10 个测点 NH3 的监测值 Cgj 干流 7 个测点 COD 的监测值 Ngj 干流 7 个测点 NH3 的监测值 五、模型的建立与求解 5.1 长江水质的综合评价 5.1.1 模糊综合评判模型[2] 根据水域情况的质量标准我们把水污染监测浓度看成是一个离散的随机变量，用概 率统计方法进行统计可以得到水域属于某个标准的概率，因为可以拟定不同的水域标 准，评价参数集为 U={u1,u2,u3},水质分级集为{v1,v2,v3,v4,v5,v6},其中 u1,u2,u3 分别表 示为溶解氧，高锰酸盐指数，氨氮（NH3-N），因为 PH 值对水域影响不大，所以对其不
2005年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎,卢亚军。许超指导教师汪晓银 予考虑,-V6分别表示为I类到劣Ⅴ类,设i参数污染物监测值共有L个,其中介 于Ar到A之间的监测值有l个, 高锰酸盐指数,氨氮(NH3-N),的监测值为(i=2,3) 监测值Xk≤A1 l=〈lA1≤Hk≤A=2,3456) 溶解氧的监测值为(i=1) l,1监测值Xk≥A1 A≤Mk≤AA(=2……n) 那么,对于i参数而言,介于A1到A之间的监测值Xk发生在A1下或A下的 概率为p=l/L1 我们用隶属度来刻画水质分级界限: 已知水质等级标准为A1,i参数污染物监测值k介于A到A之间的l个监测值 的平均值为=(∑Xkl(k=1,2.bk为第i个参数污染物位于A到A之间的监 测值的第k个。则对;参数而言,对第j级水质和j-1级水质的隶属度m和r,可由 下列各式求得(其中产=2,3),即: r12:≤Acru=A}4j/A1-A Aij-l<<Aij(Ii,n+1,n=0 ≤A ri,1,1=0z,1=0 b≈00 而仁=1时对于溶解氧的隶属度的求法与上面方法相反 对于评价参数的权重的确定: 对于溶解氧权重按如下确定w=(×0-x1)x0-s)而高锰酸盐指数,氨氮的权重分别 为w=xs1,其中ⅹ-—-第ⅰ种污染物的实测浓度算术平均值,x0-—溶解氧在某条件下 的饱和浓度(标准浓度),s-—-第ⅰ种污染物各级标准的算术平均值。对其进行归 化处理得到a=w/∑;(r1,2,3) 3个参数构成权重矩阵即为A=(a1,a2,a3) 三个指标的权重(见表1):
2005 年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎，卢亚军。许超 指导教师 汪晓银 4 予考虑，v1-- v6 分别表示为Ⅰ类到劣Ⅴ类，设 i 参数污染物监测值共有 Li个，其中介 于 Ai,j-1到 Ai,j之间的监测值有 li,j 个， 高锰酸盐指数，氨氮（NH3-N），的监测值为(i=2,3) li,1 监测值 Xi,k ≤Ai,1 li,j= li,j Ai,j-1≤Xi,k≤Ai,,j(j=2,3,4,5,6) li,n+1 Xi,k ≥Ai,n(k=1,2,… … , Li) 溶解氧的监测值为(i=1)： li,1 监测值 Xi,k≥Ai,1 li,j= li,j-1 Ai,j-1≤Xi,k≤Ai,j(j=2,… … ,n) li,n Xi,k≥Ai,n 那么，对于 i 参数而言，介于 Ai,j-1 到 Ai,j 之间的监测值 Xi,k 发生在 Ai,j-1 下或 Ai,j 下 的 概率为 pi,j= li,j / Li 我们用隶属度来刻画水质分级界限： 已知水质等级标准为 Ai,j, i 参数污染物监测值 Xi,k 介于 Ai,j-1 到 Ai,j 之间的 li,j 个监测值 的平均值为 zi,j=(∑Xi,k)/ li,j(k=1 ,2,…li,j), Xi,k 为第 i个参数污染物位于 Ai,j-1 到 Ai,j 之间的监 测值的第 k 个。则对 i 参数而言，zi,j 对第 j 级水质和 j-1 级水质的隶属度 ri,j,j 和 ri,j,j-1,可由 下列各式求得(其中 i=2,3)，即： ri,1,1=1 zi,1≤Ai,j ri,j,j-1=Ai,j-zi,j／Ai,j- Ai,j-1 Ai,j-1≤zi,j≤Ai,j ri,n+1,n=0 ri,j,j=zi,j 一 Ai,j/／Ai,j 一 Ai,j-1 Ai,j-1≤zi,j≤Ai,j (zi,n+1=0) ri,j,j-1=0 zi,j=0 ri,n+1,n=1 ri,1,1=0 zi,1=0 ri,j,j=0 (zi,n+1> Ai,n 而 i=1 时对于溶解氧的隶属度的求法与上面方法相反 对于评价参数的权重的确定： 对于溶解氧权重按如下确定 w1=(x0-x1)/(x0-s1)， 而高锰酸盐指数，氨氮的权重分别 为 wi=xi/si ，其中 xi---第 i 种污染物的实测浓度算术平均值，x0---溶解氧在某条件下 的饱和浓度（标准浓度） ，si---第 i 种污染物各级标准的算术平均值。对其进行归一 化处理得到 ai= wi/∑wi（i=1,2,3） 3 个参数构成权重矩阵即为 A=（a1,a2,a3） 三个指标的权重（见表 1）：
2005年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎,卢亚军。许超指导教师汪晓银 指标溶解氧 高锰酸盐指数氨氮 0.2523 0.6854 0.053 i参数发生各级水质的模糊概率为qi=F·pu+r1+1·P+(=1,23j=1,23456) 则模糊概率关系矩阵为q3×6于是每个单指标模糊概率评价矩阵为q×6=A*q36(见附表 1)其中q=Σa*q即为水域水体出现j级水质的模糊综合概率,并按最大隶属原则求 得q判断属于第几个级别。 在我们已知条件中我们对28个月中17个长江流域主要城市的水质监测项目求出平均值按以 上说明做出模糊综合评判,由 matlab程序求出17个城市28个月的平均水质情况,得到长江沿岸 17个观测点近两年多的水质平均情况(见表2): 表2 测点1234567891011121314151617 水质等级ⅡIIⅡⅡⅡⅡⅡIvⅢⅢIⅢⅡⅡ 劣 IlⅡ 5.1.2综合指数模型 为了验证我们评价的客观性,我们现在用综合指数评价模型对模糊综合评价模型进 行验证。 ①首先求分指数K 4个指标的K分别为:P/PH;D/DO;c/CODb;N/NH 其中选取第三类等级水的4项标准P4=7,DO=5,COD=6,NH=1为评价标准。由 于D0是逆指标,P是非单调指标(不是值越大表示水质越好或越坏)。所以综合指标计 算公式如下: P-/ P (不加权的) COD NH ②求权数n 权重的确定是参照附件3中的数据,通过统计28组数据中“主要污染指标中”各 指标出现次数,依据次数越多,权重越大的思想计算各权重(见表3) 表3 PH DO COD NHo 1.1% 25.9% 16.5% 56.5% W= n P n n (加权) Nh 因为有28个月,这样我们就得到了一个17×28的综合指标矩阵,数据过于庞大
2005 年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎，卢亚军。许超 指导教师 汪晓银 5 表 1 指标 溶解氧 高锰酸盐指数 氨氮 a 0.2523 0.6854 0.053 i 参数发生各级水质的模糊概率为 qi,j=ri,j,j * pi,j+ ri,j+1,j * pi,j+1(i=1,2,3,j=1,2,3,4,5,6), 则模糊概率关系矩阵为 q3×6,于是每个单指标模糊概率评价矩阵为 q1×6=A* q3×6（见附表 1）,其中 qj=∑ai *qi,j，即为水域水体出现 j 级水质的模糊综合概率，并按最大隶属原则求 得 qj 判断属于第几个级别。 在我们已知条件中我们对 28 个月中 17 个长江流域主要城市的水质监测项目求出平均值按以 上说明做出模糊综合评判，由 matlab 程序求出 17 个城市 28 个月的平均水质情况，得到长江沿岸 17 个观测点近两年多的水质平均情况（见表 2）： 表 2 测点 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 水质等级 Ⅱ Ⅱ Ⅱ Ⅱ Ⅱ Ⅱ Ⅱ IV Ⅱ Ⅱ Ⅰ Ⅲ Ⅱ Ⅱ 劣 Ⅴ Ⅱ Ⅱ 5.1.2 综合指数模型[3] 为了验证我们评价的客观性，我们现在用综合指数评价模型对模糊综合评价模型进 行验证。 ①首先求分指数 K 4 个指标的 K 分别为：Pj / PH0 ； Dj / DO0 ； Cj / COD0 ；Nj / NH0 其中选取第三类等级水的 4 项标准 PH0=7，DO0=5，COD0=6，NH0=1 为评价标准。由 于 DO 是逆指标，PH 是非单调指标（不是值越大表示水质越好或越坏）。所以综合指标计 算公式如下： W= 1 0 - PH Pj + D j DO0 + COD0 C j + NH0 N j （不加权的） ②求权数η 权重的确定是参照附件 3 中的数据，通过统计 28 组数据中“主要污染指标中”各 指标出现次数，依据次数越多，权重越大的思想计算各权重（见表 3） 表 3 PH DO0 COD0 NH0 η 1.1% 25.9% 16.5% 56.5% W=η1 1 0 - PH Pj +η2 D j DO0 +η3 COD0 C j +η4 NH0 N j （加权） 因为有 28 个月，这样我们就得到了一个 17×28 的综合指标矩阵，数据过于庞大
2005年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎,卢亚军。许超指导教师汪晓银 不便直观分析,所以我们按照原题附件4中说明(3)的划分标准把28个月划分为 2003年丰水期、2003年平水期;2004年枯水期、2004年丰水期、2004年平水期;2005 年枯水期、2005年丰水期、2005年平水期。每个时期按所包含月份求平均,这样就大 大缩小了数据量。用直观的柱形图表示如下图(按年评价图见附件图图3): 水质综合指数评价直观图 0.3 0.25 0.2 0.15 口枯水期平均 0.1 彐丰水期平均 0.05 口平水期平均 把题目中分类标准中的前5类标准作为5个样本,也按上述综合指标的计算方法计 算的出6个等级水质的量化标准(见表4) 表4 劣V 0.030.060.10.160. 0.22 从图1中不难看出各个测点在3个时期的水质等级,比如:江西南昌滁槎的水质最 差,在枯水期和平水期已经达到了劣V类标准。在丰水期也已达到Ⅳ类标准 然后我们再把3个时期的平均指标再综合,得出长江最近一两年的水质综合评价指 标(见图2) 3个水期综合评价 0.24 0.06 1234567891011121314151617 17个观测点 图2 根据上表中的量化标准可以得出水质等级评价(见表5)
2005 年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎，卢亚军。许超 指导教师 汪晓银 6 不便直观分析，所以我们按照原题附件 4 中说明（3）的划分标准把 28 个月划分为 2003 年丰水期、2003 年平水期；2004 年枯水期、2004 年丰水期、2004 年平水期；2005 年枯水期、2005 年丰水期、2005 年平水期。每个时期按所包含月份求平均，这样就大 大缩小了数据量。用直观的柱形图表示如下图（按年评价图见附件图[1] ~图[3]）： 水质综合指数评价直观图 0 0.05 0.1 0.15 0.2 0.25 0.3 四川攀枝花龙洞 重庆朱沱 ? 湖北宜昌南津关 ? 湖南岳阳城陵矶 ? 江西九江河西水厂 ? 安徽安庆皖河口 ? 江苏南京林山 ? 四川乐山岷江大桥 ? 四川宜宾凉姜沟 ? 四川泸州沱江二桥 ? 湖北丹江口胡家岭 ? 湖南长沙新港 ? 湖南岳阳岳阳楼 ? 湖北武汉宗关 ? 江西南昌滁槎 ? 江西九江蛤蟆石 ? 枯水期平均 丰水期平均 平水期平均 图 1 把题目中分类标准中的前 5 类标准作为 5 个样本，也按上述综合指标的计算方法计 算的出 6 个等级水质的量化标准（见表 4） 表 4 Ⅰ Ⅱ Ⅲ Ⅳ Ⅴ 劣Ⅴ 0.03 0.06 0.1 0.16 0.22 ＞0.22 从图 1 中不难看出各个测点在 3 个时期的水质等级，比如：江西南昌滁槎的水质最 差，在枯水期和平水期已经达到了劣Ⅴ类标准。在丰水期也已达到Ⅳ类标准。 然后我们再把 3 个时期的平均指标再综合，得出长江最近一两年的水质综合评价指 标（见图 2） 3个水期综合评价 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0.27 0.3 1 2 3 4 5 6 7 8 9 1011121314151617 17个观测点 图 2 根据上表中的量化标准可以得出水质等级评价（见表 5）
2005年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎,卢亚军。许超指导教师汪晓银 表5 测点1234567891011121314151617 水质等级ⅡⅡⅡⅡⅡⅡⅡⅢⅡⅢIⅢⅡⅡ劣VⅢⅡ 比较综合指数模型和模糊综合评判模型得出的结论,发现两者基本一致,所以可以 判定我们模型是客观有效的 52主要污染物COD和NH3-N的污染源分析 5.2.1因为支流互不干扰且不受干流影响,所以支流可以先分离出来单独考虑。 我们先找出长江支流10个测点从2004年4月到2005年4月13个月中COD和NH 两项指标的监测数据,然后求平均得到支流每个测点COD和NH的监测值C3和N(见 附件图), 又由于COD和MH在4项指标中的权重为16.5%和56.5%,所以它们两者之间的权重 可以重新划分为0.23和0.77。我们在把C3和Nz分别进行归一,然后再加权平均得到 综合指标W W=0.23×C3/C2+0.77×N/∑N 用柱形图表示(见下图3) 支流测点本地污染综合评价Wj 0.3 0.2 0.15 0.05 0 图3 我们取W=005为比较标准。 可以看出污染较严重的测点依次是:测点8(江西南昌滁槎)、测点3(四川泸州沱 江二桥)、测点1(四川乐山岷江大桥)、测点5(湖南长沙新港)、测点6(湖南岳阳楼) 因为各个支流之间是互不干扰的,所以我们可以初步判断这5个地区为主要污染源 522然后我们再来考虑干流上的污染源。 我们假设同一个测点上的监测值在几天之内是不变的,也就是说,可以用同一时刻 上上游的监测值经过几天降解后的剩余两表示其现时刻对下游的影响。 ①对干流各个测点间降解天数(次数)m的计算: 各测点水流速度V(13个月平均)。(产1~7)(单位由秒换算成天) V,+1 测点之间水流平均速度:Vm+1 (产1~6) 2 各测点间距离:Sj,州1 7
2005 年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎，卢亚军。许超 指导教师 汪晓银 7 表 5 测点 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 水质等级 Ⅱ Ⅱ Ⅱ Ⅱ Ⅱ Ⅱ Ⅱ Ⅲ Ⅱ Ⅲ Ⅰ Ⅲ Ⅱ Ⅱ 劣Ⅴ Ⅱ Ⅱ 比较综合指数模型和模糊综合评判模型得出的结论，发现两者基本一致，所以可以 判定我们模型是客观有效的 5.2 主要污染物 COD 和 NH3-N 的污染源分析 5.2.1 因为支流互不干扰且不受干流影响，所以支流可以先分离出来单独考虑。 我们先找出长江支流 10 个测点从 2004年 4 月到 2005年 4 月 13 个月中 COD 和 NH 两项指标的监测数据，然后求平均得到支流每个测点 COD 和 NH 的监测值 Czj 和 Nzj（见 附件图[4]）， 又由于 COD 和 NH 在 4 项指标中的权重为 16.5%和 56.5%，所以它们两者之间的权重 可以重新划分为 0.23 和 0.77。我们在把 Czj 和 Nzj 分别进行归一，然后再加权平均得到 综合指标 Wj。 Wj=0.23×Czj /∑Czj+0.77× Nzj/∑Nzj 用柱形图表示（见下图 3） 支流测点本地污染综合评价Wj 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 测点1 2 3 4 5 6 7 8 9 10 、 图 3 我们取 Wj=0.05 为比较标准。 可以看出污染较严重的测点依次是：测点 8（江西南昌滁槎）、测点 3（四川泸州沱 江二桥）、测点 1（四川乐山岷江大桥）、测点 5（湖南长沙新港）、测点 6（湖南岳阳楼）。 因为各个支流之间是互不干扰的，所以我们可以初步判断这5个地区为主要污染源。 5.2.2 然后我们再来考虑干流上的污染源。 我们假设同一个测点上的监测值在几天之内是不变的，也就是说，可以用同一时刻 上上游的监测值经过几天降解后的剩余两表示其现时刻对下游的影响。 ① 对干流各个测点间降解天数（次数）m 的计算： 各测点水流速度 Vj （13 个月平均）。（j=1～7）（单位由秒换算成天） 测点之间水流平均速度： 2 1 , 1 + + + = j j j j V V V （j=1～6） 各测点间距离：Sj，j+1；
2005年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎,卢亚军。许超指导教师汪晓银 所以 m)=S#I/V,,+l ②根据表中数据算出干流各测点COD和N13个月的平均监测值Cg和Ng 因为没有数据,我们在这里不考虑第一个测点上游对它的影响,认为它的监测值就 等于本地污染。从第2点开始计算本地污染(降解): C= Car C1×(1-02 Ngr-1×(1-0.2) 所以本地污染CD和N的综合评价指标形为(先把Cg和Ng进行归一,然后加 权平均): W=023×Cg/ΣCg+0.7×NgΣNg 用柱形图表示(见图4) 干流测点本地污染综合评价W 0.25 0.15 0.05 0 测点 图4 因为干流水质整体情况优于支流,所以在确定干流污染源时不能用支流的比较标 准。我们取W=015。得出干流污染源为:测点4(湖南岳阳城陵矶)、测点2(重庆朱 沱)。 5.3对长江未来水质污染的发展趋势的预测分析 由于河流水质状况存在许多的影响因素,而题中所给数据又十分有限,系统存在很 强的不明确的因素,于是我们考虑用灰色系统法对其进行分析和预测。, 灰色系统皚M(1,1)模型是依据系统中已知的多种因素的综合资料,将此资料的时 间序列按微分方程拟合去逼近上述时间序列所描述的动态过程,进而外推,达到预测的 目的。这种拟合得到的模型是时间序列的一阶微分方程,因此,简记为GM(1,1)模型。 建立GM(1,1)模型 给定原始时间1995-2004年资料列 x=(x(1),x(2),x(3),x(4),x0(5),x(6),x(7),x(8),x(9),x(10),) =(2564,1533,1225,1256,1571,1705,2013,1572,1824,473), 对x0坐AG0生成,有x=AGOx,x"(k)=∑x(m),则 =(x(1),x(2),x(3),x(4),x(5),x(6),x(7),x"(8),x"(9),x0(10), (11),x0(12) 8
2005 年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎，卢亚军。许超 指导教师 汪晓银 8 所以 mj = Sj+1/V j , j +1 ②根据表中数据算出干流各测点 COD 和 NH 13 个月的平均监测值 Cgj 和 Ngj 因为没有数据，我们在这里不考虑第一个测点上游对它的影响，认为它的监测值就 等于本地污染。从第 2 点开始计算本地污染（降解）： Cg * j= Cgj- Cgj-1× mj (1- 0.2) Ngj * = Ngj- Ngj-1× mj (1- 0.2) 所以本地污染 COD 和 NH 的综合评价指标 Wj 为（先把 Cg * j 和 Ngj *进行归一，然后加 权平均）： Wj=0.23×Cg * j /∑Cg * j+0.77× Ngj * /∑Ngj * 用柱形图表示（见图 4） 干流测点本地污染综合评价Wj 0 0.05 0.1 0.15 0.2 0.25 1 2 3 4 5 6 7 测点 图 4 因为干流水质整体情况优于支流，所以在确定干流污染源时不能用支流的比较标 准。我们取 Wj=0.15。得出干流污染源为：测点 4（湖南岳阳城陵矶）、测点 2（重庆朱 沱）。 5.3 对长江未来水质污染的发展趋势的预测分析 由于河流水质状况存在许多的影响因素，而题中所给数据又十分有限，系统存在很 强的不明确的因素，于是我们考虑用灰色系统法对其进行分析和预测。， 灰色系统 GM（1，1）模型是依据系统中已知的多种因素的综合资料，将此资料的时 间序列按微分方程拟合去逼近上述时间序列所描述的动态过程，进而外推，达到预测的 目的。这种拟合得到的模型是时间序列的一阶微分方程，因此，简记为 GM(1，1)模型。 建立 GM(1，1)模型 给定原始时间 1995---2004 年资料列： x (0) =（x (0)(1)，x (0)(2)，x (0)(3)，x (0)(4)，x (0)(5)，x (0)(6)，x (0)(7)，x (0)(8)，x (0)(9)，x (0)(10)，） =（2564,1533,1225,1256,1571,1705,2013,1572,1824,473）， 对 x (0) 坐 AGO 生成，有 x (1) =AGOx(0)，x (1)(k)= å = k m 1 (0) x (m), 则 x (1) =（x (1)(1)，x (1)(2)，x (1)(3)，x (1)(4)，x (1)(5)，x (1)(6)，x (1)(7)，x (1)(8)，x (1)(9)，x (1)(10)， x (1)(11)，x (1)(12)）
2005年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎,卢亚军。许超指导教师汪晓银 =(2564,4097,5322,6578,8149,9854,11867,13439,15263,15726), 对上述x的GM(1,1)参数a,b,按下述算式辨识 (BB)"By 基于x0x0与x,有 )1「-0.5(x"(1)+x0(2) 29308 0.5(x(1)+x(2)1-489155 05(x(1)+x(2)1 68581.51 5(x(1)+x(2) B=|-=0(6)1+|-05(x(1)+x"(2)1 107895 (7)1-05(x(1)+x"(2)1 135943.5 05(x"(1)+x"(2)1-1683695 z"(9)1-05x"(1)+x"(2)1-20484851 0.5(x(1)+x(2)1 2443481 y=[x(2),x(3),x(4),x0(5),x(6),x(7),x(8),x0(9),x0(10)]= (2564,1533,1225,1256,1571,1705,2013,1572,1824,473)t。 将B,y代入辨识算式,有 (BB)B'y 4115.64 =-0.146999 b=4115.64 得GM(1,1)模型为 ①灰微分方程 x0(k)-0.146999z°(k)=4115.64 ②白化方程 x() 0.1469992(k)=4115 ③白化方程的时间响应式 x(+1)=(x0-bm+b=254357010695, (t+1) (t), 1)=2564① 得还原方程:x^(+1)=1515e1469 预测结果(见表6)
2005 年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎，卢亚军。许超 指导教师 汪晓银 9 =（2564,4097,5322,6578,8149,9854,11867,13439,15263,15726）， 对上述 x (0)的 GM（1，1）参数 a,b，按下述算式辨识： =(BT B)-1 B T yN 基于 x (0) (0) x 与 x (1)，有 ú ú ú ú ú ú ú ú ú ú ú ú û ù ê ê ê ê ê ê ê ê ê ê ê ê ë é - - - - - - - - - = ú ú ú ú ú ú ú ú ú ú ú ú û ù ê ê ê ê ê ê ê ê ê ê ê ê ë é - + - + - + - + - + - + - + - + - + + ú ú ú ú ú ú ú ú ú ú ú ú û ù ê ê ê ê ê ê ê ê ê ê ê ê ë é - - - - - - - - - = 244348 1 204848.5 1 168369.5 1 135943.5 1 10789.5 1 86730 1 68581.5 1 48915.5 1 29308 1 0.5( (1) (2)) 1 0.5( (1) (2)) 1 0.5( (1) (2)) 1 0.5( (1) (2)) 1 0.5( (1) (2)) 1 0.5( (1) (2)) 1 0.5( (1) (2)) 1 0.5( (1) (2)) 1 0.5( (1) (2)) 1 (10) 1 (9) 1 (8) 1 (7) 1 (6) 1 (5) 1 (4) 1 (3) 1 (2) 1 (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) (1) x x x x x x x x x x x x x x x x x x z z z z z z z z z B yN=[ x (0)(2)，x (0)(3)，x (0)(4)，x (0)(5)，x (0)(6)，x (0)(7)，x (0)(8)，x (0)(9)，x (0)(10)]T = （2564,1533,1225,1256,1571,1705,2013,1572,1824,473）T。 将 B,yN代入辨识算式，有 . ^ ú û ù ê ë é = b a a =(BT B)-1 B T yN = ú û ù ê ë é 4115.64 - 0.146999 ， a=－0.146999 b=4115.64 得 GM(1，1)模型为 ①灰微分方程 x (0)(k)－0.146999z(1)(k)=4115.64； ②白化方程 dt dx(1) －0.146999z(1)(k)= 4115.64； ③白化方程的时间响应式 = 25433.57 t e 0.1062105 －27997.57， ( 1) ( 1) ( ), (1) (1) ( 1) ^ ( 0 ) ^ (1 ) ^ (1 ) ^ ( 0 ) ^ x t + = x t + - x t x = x =2564 ① 得还原方程 ： ( 1) ^(0) x t + =1515 t e -0.146999 ② 预测结果（见表 6） a b e a b x t x at + = - + - ( 1) ( (1) ) (0) (1) ^
2005年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎,卢亚军。许超指导教师汪晓银 年分 1995199619971998199920002001200220032004 实际值25641533 1571 705201315721824473 预测值256415441523[124014831463184314241705|387 00.0070.243-0.012-0.056-0.141-0.084-0.094-0.065-0.181 实际值4225203025822645121139944112891554|1596410602 预测值42252303291128477931918710641123261427812539 0013401270076-0345-00760057-0208-01050182 m 实际值L24524998435757861073811079211855100212048115725 预测值|2452554152876134909591861042311827134201527 00.10802130060-0153-0148-012001800113-0031 实际值3879741329905288850374784354424585833 预测值|387869 855260147013463340026105318 0-0.1070169-0055009900660276-00400061-0088 「实际值「297191259818618971320187613622462325 预测值297175300205152312401503182122072675 0-008301540102-0197-0060-0198|0602-00170150 实际值03113971801239160823209353839734454 预测值032633 l12148219762635335146835244 00.048-0.01-0.370.1960.2280.131-0.050.1780.177 注:O预测值=x()-x(k),3+om1、40 )x(k) (k) ②因1998年发洪水,数据做了相应调整 ③5为相对误差(单位%) 把预测值转换成比例(见表7): 表7 Ⅱ类 Ⅲ类 Ⅳ类 V类劣V类 20050.0246460.3451580.3112870.1105490.0583910.149968 20060.0206740.3401020.3004550.1085040.0601830.170082 20070.0172840.3338180.2888870.1060830.06180.192128 20080.0143990.3263090.2766230.1032830.0632190.216166 20090.0119350.3176290.2637670.1001360.064360.242173 20100.0098530.3077860.2503720.0966480.0652490.270093 20110.0080950.2968520.2365460.0928510.0658350.299822 20120.0066170.2849160.2223950.0887640.0661030.331206 20130.005380.2720790.2080410.0844290.066040.364031 201440.0043620.2584770.1936030.079890.0656350.398033 未来10年长江水质污染情况趋势用折线图表示(见图5)
2005 年全国大学生数学建模竞赛二等奖获奖论文(华中农业大学) 蔡虎，卢亚军。许超 指导教师 汪晓银 10 表 6 注：①预测值= x (k) x (k) (0) ^ (0) - ，?= % x (k) x (k) x (k) (0) (0) ^ (0 ) - ②因 1998 年发洪水，数据做了相应调整 ③ ξ 为相对误差 （单位%） 把预测值转换成比例（见表 7）： 表 7 年份 I 类 Ⅱ类 Ⅲ类 Ⅳ类 Ⅴ类 劣Ⅴ类 2005 0.024646 0.345158 0.311287 0.110549 0.058391 0.149968 2006 0.020674 0.340102 0.300455 0.108504 0.060183 0.170082 2007 0.017284 0.333818 0.288887 0.106083 0.0618 0.192128 2008 0.014399 0.326309 0.276623 0.103283 0.063219 0.216166 2009 0.011935 0.317629 0.263767 0.100136 0.06436 0.242173 2010 0.009853 0.307786 0.250372 0.096648 0.065249 0.270093 2011 0.008095 0.296852 0.236546 0.092851 0.065835 0.299822 2012 0.006617 0.284916 0.222395 0.088764 0.066103 0.331206 2013 0.00538 0.272079 0.208041 0.084429 0.06604 0.364031 2014 0.004362 0.258477 0.193603 0.07989 0.065635 0.398033 未来 10 年长江水质污染情况趋势用折线图表示（见图 5） 年分 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 实际值 2564 1533 1225 1256 1571 1705 2013 1572 1824 473 预测值 2564 1544 1523 1240 1483 1463 1843 1424 1705 387 I 类 ξ 0 0.007 0.243 -0.012 -0.056 -0.141 -0.084 -0.094 -0.065 -0.181 实际值 4225 2030 2582 2645 12113 9944 11289 15574 15964 10602 预测值 4225 2303 2911 2847 7931 9187 10641 12326 14278 12539 Ⅱ 类 ξ 0 0.134 0.127 0.076 -0.345 -0.076 -0.057 -0.208 -0.105 0.182 实际值 2452 4998 4357 5786 10738 10792 11855 10022 12048 15725 预测值 2452 5541 5287 6134 9095 9186 10423 11827 13420 15227 Ⅲ 类 ξ 0 0.108 0.213 0.060 -0.153 -0.148 -0.120 0.180 0.113 -0.031 实际值 387 974 1329 905 2888 5037 4784 3544 2458 5833 预测值 387 869 1554 855 2601 4701 3463 3400 2610 5318 Ⅳ 类 ξ 0 -0.107 0.169 -0.055 -0.099 -0.066 -0.276 -0.040 0.061 -0.088 实际值 297 191 259.8 186 1897 1320 1876 1136 2246 2325 预测值 297 175 300 205 1523 1240 1503 1821 2207 2675 Ⅴ 类 ξ 0 -0.083 0.154 0.102 -0.197 -0.060 -0.198 0.602 -0.017 0.150 实际值 0 311 339.7 180 1239 1608 2329 3538 3973 4454 预测值 0 326 334 112 1482 1976 2635 3351 4683 5244 劣 Ⅴ 类 ξ 0 0.048 -0.01 -0.37 0.196 0.228 0.131 -0.05 0.178 0.177