1、长江水质监测摘要本文解决的是长江水质的评价与监测问题，通过分析过去十年不同监测站收集到的长江水质数据，运用不同的理论建立不同的模型，对长江过去十年的水质情况作出评价，然后再预测未来十年长江水质的变化情况。针对问题一：考虑到问题一中需要对长江水质情况作出定量的评价，并分析各地区水质的污染状况，为此，建立模糊综合评价模型确定了其隶属度函数，建立评判因子的权重矩阵，求得最终结果为：水质最差的地方是江西南昌滁槎（15 号），其次水质差的地方为四川乐山岷江大桥（8 号）、湖南长沙新港（12 号）以及四川泸州沱江二桥（10 号），此四处水质污染严重；水质最好的地方是湖北丹江口胡家岭 （11 号）。针对问题
2、二：根据长江的降解系数，可得到污染物随时间的变化量。由于污染源的污染物排放量等于本地区污染物的流量与上游流下的污染物流量之差。因此，建立污染物流量随时间变化的微分方程模型。最后求得：高锰酸钾指数和氨氮的污染源主要集中在宜昌至岳阳之间。针对问题三：根据已知的过去 10 年的主要统计数据，建立了灰色预测模型。在相对误差较小的情况下对未来 10 年的水质情况作出了预测，分析得出结论：未来 10 年可饮用水所占的比例越来越低，排污量有明显的上升趋势。针对问题四：在问题四中建立多元线性回归方程，利用最小二乘法求解系数，在满足问题四要求的前提下，求出未来 10 年的允许最大相对排污量，继而求得未来 10
3、年每年的相应排污量，后者与前者的差值与未来 10 年的长江水总流量的乘积，求得最终结果如下表：未来 10 年预处理的排污量年代 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014预处理排污量(亿吨)71.24 83.11 94.98 106.86 118.73 130.60 142.48 154.35 166.22 178.09针对问题五：分析总结前几个问题的结果，找出水质污染的根本原因。结合考察团的调查结果，给出合理的建议和意见。 最后，对模型中运用的方法进行了优、缺点评价，在模型的推广中提出了可以建立类似模型解决生活中的一类问题。关键词：模糊
4、评价 微分方程 灰色预测 线性回归11 问题重述1.1 问题背景长江是我国第一、世界第三大河流，其水质的污染程度日趋严重，已引起了相关政府部门和专家们的高度重视。2004 年 10 月，由全国政协与中国发展研究院联合组成“保护长江万里行”考察团，从长江上游宜宾到下游上海，对沿线 21 个重点城市做了实地考察，揭示了一幅长江污染的真实画面，其污染程度让人触目惊心。为此，专家们提出“若不及时拯救，长江生态 10 年内将濒临崩溃”。 长江水是许多人赖以生存和发展的资源，保护水资源就是保护我们自己，就长江近年来的水质情况，采取合理的保护和治理措施刻不容缓。1.2 解决问题为了制定出合理的治理长江水质污
5、染的方案，根据长江地区近两年多主要水质指标的检测数据，现讨论以下几个问题： (1) 对长江近两年多的水质情况做出定量的综合评价，并分析各地区水质的污染状况。(2) 研究、分析长江干流近一年多主要污染物高锰酸盐指数和氨氮的污染源的主要分布地区。(3) 假如不采取更有效的治理措施，依照过去 10 年的主要统计数据，对长江未来水质污染的发展趋势做出预测分析。(4) 根据预测分析，如果未来 10 年内每年都要求长江干流的类和类水的比例控制在 20%以内，且没有劣类水, 那么每年需要处理多少污水？ (5) 发表对解决长江水质污染问题切实可行的建议和意见。2 问题分析问题需要对长江水质作出评价和预测，根据
6、附件中已知的数据，建立相应2的评价预测模型，分析得出长江水质过去十年以及未来十年的污染状况。2.1 问题一的分析问题提供的数据主要反映了三方面的内容即：水质划分等级标准、17 个测站点在 28 个月份中主要污染物的浓度和水质等级。而观测点水质等级的确定是由主要污染物的浓度决定的，因此为了对长江及各观测点作出定量的综合评价我们需要将主要污染物的浓度归一到一个比较量中。这个比较量要能够快速而准确的反映出各测站点的水质污染状况，因此我们评价整个长江流域的水质状况采用模糊综合评价模型。 虽然 PH 值是影响水质等级的因素之一，但是给出的所有 PH 值都在 69 这个范围内，即 PH 值对于水质等级几乎
7、没有影响，所以在以后的问题讨论中不再考虑 PH 值的影响。2.2 问题二的分析因为污染物在时间和空间上是动态变化的，为简便起见，不考虑支流的因素，只从纵向的角度，结合本地与上游污水进行水质的污染分析。根据常识可知，污染物的污染源就是新增污染物较多的地方。因此，考虑到降解系数，求得污染物随时间的变化量，通过相邻两主要观测点之间每千米的相对排污量判定污染源，即本地新增污染物的量等于本地区现有的污染物减去上游通过降解后到达该处的污染物的总量，再比上两地之间的距离，得到相对排污量。建立污染物的量随时间变化的微分方程模型，最后通过新增污染物的相对量大小关系得出主要的污染源地区。2.3 问题三的分析测第三
8、问根据过去 10 年长江的总体水污染状况的监测数据，可以看出长江总体水流量变化不大，但年排污的总量在增加，这使得污染河段比例增加，污染的严重程度呈现快速增长的趋势，即每年污染情况主要与当年的排污量和总水流量等因素有关。为此，首先利用回归分析方法确定出可饮用水的比例与总排污量和总水流量的关系式，然后根据过去 10 年排污量，利用灰色预测方法对未来的年排污量做出预测，最后根据总排污量的增长趋势来推断出可饮用水比例的变化趋势，从而可以预测出未来 10 年长江水质的变化情况。32.4 问题四的分析在问题四中建立多元线性回归方程，利用最小二乘法求解系数，在满足未来 10 年内每年都要求长江干流的类和类水
9、的比例控制在 20%以内，且没有劣类水的前提下，求出未来 10 年的允许最大相对排污量，继而求得未来 10年每年的相应排污量，后者与前者的差值与未来 10 年的长江水总流量的乘积，即可求得未来 10 年要处理的排污量。2.5 问题五的分析分析总结前几个问题的结果，找出水质污染的原因。结合考察团的调查结果，给出合理的建议和意见。3 模型假设与符号说明3.1 模型假设假设一： 降解系数在一定时间段固定不变；假设二：长江干流的自然净化能了近似均匀； 假设三：干流的相邻两观测站点之间的排污口主要集中在下游处；假设四：干流污染物的富集主要受上游影响，支流的影响忽略不计。3.2 符号说明符号 符号意义ui
10、A氧气属于第 i 等级的隶属度uiB高锰酸盐属于第 i 等级的隶属度uiC氨氮属于第 i 等级的隶属度4ia 第 i 种指标的权重jR各观测点模糊关系矩阵ix第 i 种指标的实时均值W平均排污量w相对排污量u断面水流平均流速C某组分子在 处的浓度xk降解常数tx1每年长江总水流量t2每年总排污量ty1V 类水所占的比例t2劣 V 类水所占比例4 模型的建立与求解4.1 问题一的解答4.1.1 模型一的建立与求解（1）确立评判指标：根据地表水环境质量标准（GB38382002）中4 个主要项目标准限值设置表，结合 28 个月的各观测点的测量结果，选取溶氧5量（DO）、高锰酸盐指数(CODMn)和
11、氨氮（NH3-N）为评价因子，设置评价因子集： ，求出 17 个观测点 28 个月的个评价因子的均值，得到 17321,xu个点的评价因子集。（2）建立评价集序号项目标 准值分类类 类 类 类 类 劣类1 溶解氧(DO) 7.5（或饱和率90%）6 5 3 2 02 高锰酸盐（CODMn) 2 4 6 10 15 3 氨氮（NH3-N） 0.15 0.5 1.0 1.5 2.0 4 PH 值（无量纲） 6-94 个指标将水质分为 6 个等级。（3）建立隶属函数，进行单因素评价：6由于水质污染程度和水质分级标准都是模糊的，所以用隶属度来描述分级界限较为合理，先根据各指标的 6 级标准，做出 6
12、个级别的隶属函数，其中 DO的评价指标以数值大为优，其余两个指标以数值小为优，函数如下：DO（A）以升半梯形分布建立隶属函数 012202023 3253520535630 65.75.705.765.7165 43 2xxAxxorA xxorAxxorA xorxxx uu uu uuCOOMn(B)和 NH3-N（ C）以升半梯形分布建立隶属函数：7 xxBxxB xxorBxxxorB xxxorBxxBuu uu uu 115105 150561361614400 4260424165 43 2 xxCxxC xxorxxxor xxxorCxxCuu uu uu 1215210 2
13、52.1.10150.01 151.005.0.01.5. 65 43 21有了各指标的隶属函数，就可以进行单因素评价，将各观测点的 3 个指标带入相应的隶属函数，计算出指标的隶属度，得到各观测点模糊关系矩阵（j=1，2，3.，17），具体数据见附录一。jR8（4）建立评价因素的权重集：由于 DO,COOMn 和 NH3-N 等污染指标对水质的影响不同，因此对各指标应赋予不同的权重，根据污染物对水质的污染大、权重的的原则来决定权重的大小。定义：对于 DO 越大越优型： ，对于 COOMn 和 NH3-N 越小越优型：iimxa，其中， 为第 i 种指标的权重， 是第 i 种指标的实时均值， 分
14、iixnaiai inm,别为多级浓度标准的最大值和最小值。再将权重归一化处理， ，评31iiiaQ价因素权重的集合为 A=0.145,0.226,0.629。（5）模糊综合评价：评价因素的权重集 A 乘以单因子矩阵 得到模糊综合评价结果，jR。将 6 个等级看作一种相对位置使其连续化，设各等级分别用jjARB1，2 ，3 ，4，5，6 表示，并称为各等级的秩，然后用 B 中对应的量将各级的秩加权求积，得到被评等级的相对位置，即：61*kkbB式中 为隶属于第 k 级的隶属度。具体计算结果见下表：kb表 1：各观测站水质等级序号 1 2 3 4 5 6 7 8 99水质级1.29911.411
15、61.40381.45691.24041.50001.16482.50001.4450序号1011121314151617水质级2.12541.00002.38881.91561.50004.15721.43241.41654.1.2 问题一的结果分析由表 1 可知，各观测站中水质等级最大的为江西南昌滁槎 （15 号），可知此处水质情况最差，其次水质差的地方为四川乐山岷江大桥（8 号）、湖南长沙新港（12 号）以及四川泸州沱江二桥（10 号），此四处水质污染严重；水质最好的地方是湖北丹江口胡家岭 （11 号）。4.2 问题二的解答因为污染物在时间和空间上是动态变化的，为简便起见，不考虑支流的因
16、素，只从纵向的角度，结合本地与上游污水进行水质的污染分析。建立污染物流量随时间变化的微分方程模型，最后通过新增污染物的流量大小关系得出主要的污染源地区。4.2.1 模型二的建立与求解（1）某一污染物扩散所满足的微分方程是一个抛物线方程，结合实际问题的假设，常假设其水流近似的处于稳定状态，断面均匀，普通对流扩散方程为kCxDutC2上式是解决污染物的一般方程，结合问题二的情况，假设观测站干流的污染物浓度在一个月内保持不变，则 ，又因为 ,所以上式进一步简化0t0为：kCxu10其中， 为断面水流平均流速， 为某组分子在 处的浓度， 为降解常数uCxk取 0.2。由此解出长江干流污染物浓度 与据观
17、测站 米处的函数关系：xukixie（2）对于长江上任一段干流 AB，A 为第 个观测点，B 为第 个观测点，1i两观测点的距离为 ，假设该段干流之间有 个排污口（包括支流入口和直1idn排口），为计算的方便，假设排污口聚集在一起，其总的排污口的流量为 ，q平均流量为 ，污染物的浓度为 ， 即为 AB 段的排污量,则ucqiuxki QeCCi11BUdkABeqcBUdkABeQCcW2BAU对于任意一个江段 AB，污染物的浓度 ，水流量 ，流速 ，BAC,BAQ,BAU,距离 和降解系数 k 都已知，则可以算出 AB 段污染物的排放量 。ABd jW（3）根据所给数据可以算出，对于每个月每
18、个江段的排污区间（j=1,2,.,13 个月），取 13 个月的平均值 表示平均排污量，jW,13jj则该江段每千米的排污量 为相对排污量，每个江段都有一个相对排污ABdWw量，以此为指标可以判断长江干流的主要污染源。114.2.2 模型二的求解按步骤运用 MATLAB 编求解（程序见附录 3），结果如下表：表 2：COOMn 的排污量江段 21s32s4s5s6s7sW32510.54 35266.77 46002 29855.62 11104.15 31404.31相对排污量 34.22 45.33 116.46 59.71 105.75 60.05排序 6 5 1 4 2 3表 3：NH
19、3-N 的排污量江段 21s2s4s54s65s76sW 2822.62 3191.46 4451.77 2559.31 1191.38 158.46相对排污量 2.97 4.10 11.27 5.12 7.26 0.34排序 5 4 1 3 2 64.2.3 问题二的结果分析问题二用长江干流 7 个观测站点将长江分为 6 个江段，通过污染物流量随时间的变化关系建立了微分方程模型，先计算出每月每段的高锰酸盐和氨氮的量，再对六个污染源近一年所排放的污染物的数量求期望，分析得出污染物高锰酸盐指数和氨氮的污染源主要集中在宜昌至岳阳之间。4.3 问题三的解答4.3.1 模型三的建立与求解当系统内部信息
20、和特性是部分已知的，另一部分是未知时，人们往往难以建立客观的物理原形，内部因素难以辨识，以及相互之间的关系较为隐蔽，难12已准确了解这类系统的行为特征，因此，对于这类问题进行定量描述，即建立模型难度较大，便选择建立灰色系统，对问题求解。1. 灰色动态模型 附件 4 中给出了 1995-2004 年长江各等级在丰水期，枯水期，水文年的百分比，为了预测此趋势之下未来十年长江水质的变化，通过对水文年全流域水质的变化研究来对长江未来水质污染的发展趋势作出预测分析。建立灰色预测模型，首先舍去异常数据：1995 年 IV 类、劣 V 类数据，然后从附件 4 中提取出水文年全流域随时间变化的原始数据。 10
21、.2,1000 nxx对于原始数据进行累加处理，即 10.2,101kixki加生成数列为 10.2,111 nxx采用一阶单变量微分方程进行拟合，得到白化方程的 GM（1,1 ）模型 utaxdt1其中式中的 为待定系数。 ua,综上所述，灰色动态模型为：ukazxkxk10115.0.2. 模型求解 取水文年全流域的饮用水预测为例，可通过对其变化研究对长江未来水质污染的发展趋势作出预测分析。 13(1) 求解步骤： 第一步：从附表中提取出数据 68,5.7,.34,2.80,7.385,1.90x第二步：利用公式对 做累加生成得0 6.79.21.654.7.501.42.371.2594
22、.781931x第三步：构造矩阵 B 和数据向量， 与 满足关系 其中 xabYnTTnnxY 685.7.3742.807.3851200 16.738590.1.3746.801.1,21,31,21nxxB第四步：求解常系数 8436.0211nTYBua第五步：得出表达式： 3846375021.11 kkak euexY14 令 10kkY所以 表示还原后的值。 0kY(2) 结果求解 按步骤运用 MATLAB 编求解（程序见附录 4），运用此程序根据过去十年的数据可以预测出未来 10 年的排污量和可饮用水如下表：表 4：未来 10 年排污量表 5：未来 10 年可饮用水所占比例4.
23、3.2 问题三的结果分析从表 5 可以看出，可饮用水所占的比例越来越低，所以，可以看出水质越来越差，从 2012 年起，可饮用水的百分比甚至低于 60%,如此严重的水污染，应该引起足够重视，采用合理的措施来减少对长江水的污染。4.4 问题四的解答4.4.1 模型的建立与求解（1）长江干流的类和类水以及劣 V 类水的比例与每年长江总水流量 和总排污量 有关，因此分别以类和类水所占的比例 和劣 V 类txtx2 1y年代 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014排污量 303.01 322.52 343.29 365.39 388.92 4
24、13.96 440.61 468.98 499.18 531.32年份 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014可饮水 69.61 68.02 66.47 64.95 63.47 62.02 60.6 58.22 56.87 54.5515水所占比例 为响应变量,以历年来的2y=tx1 9405,8102,.89,249513,6.7,95130 =t2 5,7,.0,28,74为解释变量，做出多元线性回归，得到类和类水的总和及劣 V 类水与长江总水流量及总排污量的关系，做出一般的多元线性回归模型：运用 matlab 用最小二乘法求解得
25、到回归系数9276.15,0.,32.,1cba8.,4.,7.,2就求得了类和类水的总和及劣 V 类水与长江总水流量及总排污量的关系。问题四要解决问题是：未来 10 年内每年都要求长江干流的类和类水的比例控制在 20%以内，且没有劣类水。则约束条件为： 09276.1507.032. 8482121 txtxy综上所述，所建立的模型为： 09276.1507.032. 848. 9276.15070322121 2121 txtxyts cttxxy对上述回归方程进行回归检验时，发现回归效果并不是太好，考虑到22121ctxbtay16长江干流的类和类水的比例与江总水流量及总排污量不是直接的
26、线性关系，因此对模型进行改进。（2）优化模型长江干流的类和类水以及劣 V 类水的 比例与每年总排污量和长江总水流量 的比值 相对排污量 有关，因此分别以txtZ12类和类水所占的比例 和劣 V 类水所占比例 为解释变量， 为响应变量，1y2ytZ做多元线性回归。对表四的数据进行观测，得到 与时间 成线性关系，因此tZ列出如下线性回归方程：ctbytatZ2121运用 matlab 用最小二乘法进行求解，得到回归系数分别为0165.,2.,8.,.,.,cba计算时剔除了 2003 年这个突变点，且检验到回归性较好。得到的回归方程如下：0165.05.02.162tytytZ21将 带入 的表达
27、式中，得到最大允许相对排污量0,21tyt t2，将 带入 的表达式中得到未来 10 年每年05.2Z,.143,tZ1预测的相对排污量，则污水处理量 ，因为每年的总Lttw0,max2水流量（除特殊年份如 1998 年）变化不大，故取前 10 年的总水流量均值为每年的总水流量 ，得到每年的污水处理量如下表：10txLtx2tx117表 6： 未来 10 年预处理的排污量年代 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014预处理排污量(亿吨)71.24 83.11 94.98 106.86 118.73 130.60 142.48 154.3
28、5 166.22 178.094.4.2 问题四的结果分析从表 6 可以看出，未来 10 年所需处理的污水越来越多，所以，可以看出水质越来越差，如此严重的水污染，应该引起足够重视，采用合理的措施来减少对长江水的污染。4.5 解决长江水质问题的建议通过问题一的计算可以看出，长江水质的枯水期、丰水期与平水期都由水流量的不同导致水质情况也不同。从问题三、四的求解结果可以看出，在未来的十年里，长江水质不断恶化，废水排放量逐年递增，需要处理的污水也越来越多。根据模型假设，在不发生特大洪涝灾害时长江的水流量是一定的，总结得出解决长江水质问题的关键，在于减少废水的排放量的同时增强废水的处理能力。鉴于此，我们
29、认为可以从以下几个方面着手解决长江水质问题： （1）加大处罚力度。（2） 对污染严重地区进行重点治理。（3） 加强宣传，增强民众保护长江的意识。（4） 严格进行检测，处理违规企业，严惩腐败。（5） 保护长江立法，引起政府和人民的高度重视。185 模型的评价、改进与推广5.1 模型的评价5.1.1 模型的优点（1）问题一中建立了模糊综合评价模型，通过精确的数字手段处理模糊的评价对象，对蕴藏信息呈现模糊性的资料作出了比较科学、合理、贴近实际的量化评价；（2）问题二的模型，简化了许多因素，忽略了支流的影响，但这种理想化还是合理的，较符合实际情况;（3）问题三采用灰色预测模型，具有少数据性、良好的时效
30、性、较强的系统性和关联性等特征，可以合理的对数据做出预测。5.1.2 模型的缺点在求解模型时，为了使问题得到方便的解决，往往采用简化的手段进行求解，因此求出的结果与真实值会存在一定的偏差。5.2 模型的改进在问题三中运用了灰色动态模型对未来 10 年长江的水质作出了预测，但并未对预测的结果进行检验，不知预测结果的准确度有多高，在此，最好能将模型三的基础上对灰色模型的预测结果进行检验，以保证预测结果有较高的可信度。5.3 模型的推广模型一中，除了用模糊算法评价外，还可以用 BP 神经网络来求长江水质的在综合评价而；而在求解问题三和问题四时采用的灰色预测，还可以用回归分析预测、时间序列预测和神经网
31、络预测等代替。196 参考文献1 宋晓秋编著，模糊数学原理与方法，北京，中国矿业大学出版社，1999 2 张俊福等编著，应用模糊数学，北京，地质出版社，20013 岳超源，决策理论与方法，北京：科学出版社，20034 邓聚龙，灰理论基础,武汉：华中科技大学出版社，20025 韩中庚，数学建模方法及其应用，北京：高等教育出版社，2005 附录一各观测点模糊关系矩阵：1 号2 号0094216 05194.0.823 号4 号0326.0.75 09.0.825 号6 号0029714 028.0.7537 号8 号001465 0846.0.159 号 10 号08013679623067.52
32、011 号12 号001 0832.0.0492.7.313 号14 号 006734095. 01357.0.6215 号 16 号 010012 8. 03897.0.17 号 03917.0.5附录二问题一的程序（以 17 号观测站为例）：clc;cleara=1 0 0 0 0 0;0.5107 0.5107 0 0 0 0;0.3917 0.3917 0 0 0 0;b=0.145 0.226 0.629;c=b*a;d=0;e=0;for i=1:6d=d+c(1,i)*i;endfor i=1:6e=e+c(1,i);endB=d/e21附录三问题二的程序：clc;clearnu
33、m11=xlsread(11.xls);num29=xlsread(29.xls);for i=1:6Ca=num11(i,6);Cb=num11(i+1,6);Qa=num29(2,i+2);Qb=num29(2,i+3);d1=num29(1,i+1);d2=num29(1,i+2);d=(d2-d1)*1000;Ua=num29(3,i+2);Ub=num29(3,i+3);U=(86400*(Ua+Ub)/2;m=exp(-0.2*d/U);W=Cb*Qb-Ca*Qa*m;Wend附录四问题三的程序：clc;clear y=input(请输入数据 );n=length(y); yy=o
34、nes(n,1);yy(1)=y(1); for i=2:n yy(i)=yy(i-1)+y(i);22end B=ones(n-1,2); for i=1:(n-1) B(i,1)=-(yy(i)+yy(i+1)/2; B(i,2)=1; endBT=B; for j=1:n-1 YN(j)=y(j+1); endYN=YN;A=inv(BT*B)*BT*YN;a=A(1); u=A(2); t=u/a; t_test=input(请输入需要预测个数：); i=1:t_test+n; yys(i+1)=(y(1)-t).*exp(-a.*i)+t; yys(1)=y(1); for j=n+
35、t_test:-1:2 ys(j)=yys(j)-yys(j-1); end x=1:n; xs=2:n+t_test;yn=ys(2:n+t_test);23plot(x,y,r,xs,yn,*-b);det=0;for i=2:n det=det+abs(yn(i)-y(i)/xs(i);end det=det/(n-1); disp(百分相对误差为：, num2str(det),%);disp(预测值为： , num2str(ys(n+1:n+t_test);附录五问题四的程序：clc;cleara=9205 9513 9171.26 13127 9513 9924 8892.8 10210 9980 9405;c=mean2(a);z=0.0002*20+0.0185;b=11 12 13 14 15 16 17 18 19 20;for i=1:10Z=0.0012*b(1,i)+0.0165;M=(Z-z)*c;Mend
