数学建模之长江水质

1、数学建模之长江水质监测问题长江水质监测摘要本文解决的是长江水质的评价与监测问题，通过分析过去十年不同监测站 收集到的长江水质数据，运用不同的理论建立不同的模型，对长江过去十年的 水质情况作出评价，然后再预测未来十年长江水质的变化情况。针对问题一：考虑到问题一中需要对长江水质情况作出定量的评价，并分 析各地区水质的污染状况，为此，建立模糊综合评价模型确定了其隶属度函数， 建立评判因子的权重矩阵，求得最终结果为：水质最差的地方是江西南昌滁 槎（15号），其次水质差的地方为四川乐山岷江大桥 （8号）、湖南长沙新港（12 号）以及四川泸州沱江二桥（10号），此四处水质污染严重；水质最好的地方是 湖北丹

2、江口胡家岭（11号）。针对问题二：根据长江的降解系数，可得到污染物随时间的变化量。由于 污染源的污染物排放量等于本地区污染物的流量与上游流下的污染物流量之 差。因此，建立污染物流量随时间变化的微分方程模型。最后求得：高锰酸钾 指数和氨氮的污染源主要集中在宜昌至岳阳之间。针对问题三：根据已知的过去10年的主要统计数据，建立了灰色预测模 型。 在相对误差较小的情况下对未来10年的水质情况作出了预测，分析得出 结论：未来10年可饮用水所占的比例越来越低，排污量有明显的上升趋势。针对问题四：在问题四中建立多元线性回归方程，利用最小二乘法求解系 数，在满足问题四要求的前提下，求出未来 10年的允许最大相

3、对排污量，继而 求得未来10年每年的相应排污量，后者与前者的差值与未来 10年的长江水总 流量的乘积，求得最终结果如下表：未来10年预处理的排污量年代2005200620072008200920102011201220132014预处理排污量（亿吨）71.2483.1194.98106.86118.73130.60142.48154.35166.22178.09针对问题五：分析总结前几个问题的结果，找出水质污染的根本原因。结 合考察团的调查结果，给出合理的建议和意见。最后，对模型中运用的方法进行了优、缺点评价，在模型的推广中提出了 可以建立类似模型解决生活中的一类问题。关键词：模糊评价 微分方

4、程 灰色预测 线性回归问题重述1.1问题背景长江是我国第一、世界第三大河流，其水质的污染程度日趋严重，已引起 了相关政府部门和专家们的高度重视。2004年10月，由全国政协与中国发展研 究院联合组成“保护长江万里行”考察团，从长江上游宜宾到下游上海，对沿 线21个重点城市做了实地考察，揭示了一幅长江污染的真实画面，其污染程度 让人触目惊心。为此，专家们提出“若不及时拯救，长江生态10年内将濒临崩溃”。长江水是许多人赖以生存和发展的资源，保护水资源就是保护我们自己， 就长江近年来的水质情况，采取合理的保护和治理措施刻不容缓。1.2解决问题为了制定出合理的治理长江水质污染的方案，根据长江地区近两年

5、多主要 水质指标的检测数据，现讨论以下几个问题：(1) 对长江近两年多的水质情况做出定量的综合评价，并分析各地区水质 的污染状况。(2) 研究、分析长江干流近一年多主要污染物高锰酸盐指数和氨氮的污 染源的主要分布地区。(3) 假如不采取更有效的治理措施，依照过去 10年的主要统计数据，对 长江未来水质污染的发展趋势做出预测分析。(4) 根据预测分析，如果未来10年内每年都要求长江干流的W类和V类 水的比例控制在20%以内，且没有劣V类水，那么每年需要处理多少污水？(5) 发表对解决长江水质污染问题切实可行的建议和意见。2问题分析问题需要对长江水质作出评价和预测，根据附件中已知的数据，建立相应

6、的评价预测模型，分析得出长江水质过去十年以及未来十年的污染状况。2.1问题一的分析问题提供的数据主要反映了三方面的内容即：水质划分等级标准、 17个测 站点在28个月份中主要污染物的浓度和水质等级。而观测点水质等级的确定是 由主要污染物的浓度决定的，因此为了对长江及各观测点作出定量的综合评价 我们需要将主要污染物的浓度归一到一个比较量中。这个比较量要能够快速而 准确的反映出各测站点的水质污染状况，因此我们评价整个长江流域的水质状 况采用模糊综合评价模型。虽然PH值是影响水质等级的因素之一，但是给出的所有PH值都在69这个范围内，即PH值对于水质等级几乎没有影响，所以在以后的问题讨论中不再 考虑

7、PH值的影响。2.2问题二的分析因为污染物在时间和空间上是动态变化的，为简便起见，不考虑支流的因 素，只从纵向的角度，结合本地与上游污水进行水质的污染分析。根据常识可 知，污染物的污染源就是新增污染物较多的地方。因此，考虑到降解系数，求 得污染物随时间的变化量，通过相邻两主要观测点之间每千米的相对排污量判定污染源， 即本地新增污染物的量等于本地区现有的污染物减去上游通过降解后到达该处 的污染物的总量，再比上两地之间的距离，得到相对排污量。建立污染物的量 随时间变化的微分方程模型，最后通过新增污染物的相对量大小关系得出主要 的污染源地区。2.3问题三的分析测第三问根据过去10年长江的总体水污染状

8、况的监测数据，可以看出长江总 体水流量变化不大，但年排污的总量在增加，这使得污染河段比例增加，污染 的严重程度呈现快速增长的趋势，即每年污染情况主要与当年的排污量和总水 流量等因素有关。为此，首先利用回归分析方法确定出可饮用水的比例与总排 污量和总水流量的关系式，然后根据过去 10年排污量，利用灰色预测方法对 未来的年排污量做出预测，最后根据总排污量的增长趋势来推断出可饮用水比 例的变化趋势，从而可以预测出未来 10年长江水质的变化情况。2.4问题四的分析在问题四中建立多元线性回归方程，利用最小二乘法求解系数，在满足未 来10年内每年都要求长江干流的W类和V类水的比例控制在20%以内，且没有劣

9、V类水的前提下，求出未来10年的允许最大相对排污量，继而求得未来10年每年的相应排污量，后者与前者的差值与未来10年的长江水总流量的乘积，即可求得未来10年要处理的排污量。2.5问题五的分析分析总结前几个问题的结果，找出水质污染的原因。结合考察团的调查结 果，给出合理的建议和意见。3模型假设与符号说明3.1模型假设假设一：降解系数在一定时间段固定不变；假设二：长江干流的自然净化能了近似均匀；假设三：干流的相邻两观测站点之间的排污口主要集中在下游处；假设四：干流污染物的富集主要受上游影响，支流的影响忽略不计。3.2符号说明符号符号意义Aui氧气属于第1等级的隶属度Bui咼锰酸盐属于第i等级的隶属

10、度Cui氨氮属于第i等级的隶属度ai第i种指标的权重Rj各观测点模糊关系矩阵Xi第i种指标的实时均值W平均排污量w相对排污量u断面水流平均流速C某组分子在x处的浓度k降解常数x1 t每年长江总水流量x2 t每年总排污量yi tV类水所占的比例y21劣V类水所占比例4模型的建立与求解4.1问题一的解答4.1.1模型一的建立与求解(1) 确立评判指标：根据地表水环境质量标准(GB383 2002)中4个主要项目标准限值设置表，结合 28个月的各观测点的测量结果，选取溶氧量(DO、高锰酸盐指数(CODMn和氨氮(NH3-N为评价因子，设置评价因子集：u X1,X2,X3 ,求出17个观测点28个月的

11、个评价因子的均值，得到17个点的评价因子集。(2) 建立评价集序 号项目标准值分类I类U类川类W类V类劣V类1溶解氧(DO)>7.5(或饱和率90%653202高锰酸盐(COD Mn)w2461015OO3氨氮(NH3-N <0.150.51.01.52.0OO4PH值(无量纲)6-94个指标将水质分为6个等级。(3) 建立隶属函数，进行单因素评价：由于水质污染程度和水质分级标准都是模糊的，所以用隶属度来描述分级 界限较为合理，先根据各指标的 6级标准，做出6个级别的隶属函数，其中DO 的评价指标以数值大为优，其余两个指标以数值小为优，函数如下：D0( A)以升半梯形分布建立隶属函

12、数x 7.50 x 7.5 or x 5Aui xx 67.5 66 x 7.5x 7.5x 6Au3 xx 6 or x 35x63x50x 34 x5 3x 2x 5or x 23x52x3Au5 x0x 2x2x 3or x 02x30x2x 20x2x 06#COOMn(B) NH3-N(C)以升半梯形分布建立隶属函数:#1x 20x 2 or x 660x10or x 40x15 or xx4x3Bu3 x4x6Bu4 x6x 1064106x6x106x1010x 1510615100x0x10Bu5 x10x15Bu6 x 1x15101x15Bu1 X2x4x 4Bu2 x2x

13、44x6cu1 x1xx0.150.15x0.15 x0.50.50.500.150x1.5 or x0.5cu3xx0.50.5 x 110.5x11x 1.51.510x1.5Cu5xx1.51.5x 221.51x20x0.15 or x 1Cu2 xx0.150.15 x 0.50.50.15x0.50.5 x 110.50x 2or x 1Cu4 xx11 x 1.51.51x1.51.5 x 221.5cu6 x1x有了各指标的隶属函数，就可以进行单因素评价，将各观测点的 3个指标 带入相应的隶属函数，计算出指标的隶属度，得到各观测点模糊关系矩阵Rj（j=1，2, 3. ，17），

14、具体数据见附录一。78(4) 建立评价因素的权重集:由于DO,COOM和NH3-N等污染指标对水质的影响不同，因此对各指标应赋 予不同的权重，根据污染物对水质的污染大、权重的的原则来决定权重的大小。疋义：对于DO越大越优型：Xa ，对于COOM和 NH3-N越小越优型: minai氏，其中，ai为第i种指标的权重，人是第i种指标的实时均值,Xi别为多级浓度标准的最大值和最小值。再将权重归一化处理,Qiai3，aii 1因素权重的集合为A=0.145,0.226,0.629(5) 模糊综合评价：评价因素的权重 集A乘以单因子矩阵Rj得到模糊 综合评价结果，Bj ARj。将6个等级看作一种相对位置

15、使其连续化， 设各等级分别用1, 2, 3 4，5，6表示，并称为各等级的秩，然后用 B中对应的量将各级的秩加权求积， 得到被评等级的相对位置，即：bk kk 16bkk 1式中bk为隶属于第k级的隶属度。具体计算结果见下表:表1 ：各观测站水质等级序号123456789水111111121质级.2991.4116.4038.4569.2404.5000.1648.5000.4450序11111111号012345671水21211411质级.1254.0000.3888.9156.5000.1572.4324.41654.1.2问题一的结果分析由表1可知，各观测站中水质等级最大的为江西南昌滁

16、槎（15号），可知此处水质情况最差，其次水质差的地方为四川乐山岷江大桥（8号）、湖南长沙新港（12号）以及四川泸州沱江二桥（10号），此四处水质污染严重；水质最 好的地方是湖北丹江口胡家岭（11号）。4.2问题二的解答因为污染物在时间和空间上是动态变化的，为简便起见，不考虑支流的因 素，只从纵向的角度，结合本地与上游污水进行水质的污染分析。建立污染物 流量随时间变化的微分方程模型，最后通过新增污染物的流量大小关系得出主 要的污染源地区。4.2.1模型二的建立与求解（1）某一污染物扩散所满足的微分方程是一个抛物线方程，结合实际问题 的假设，常假设其水流近似的处于稳定状态，断面均匀，普通对流扩散方

17、程为C 止 D$ kCt xx上式是解决污染物的一般方程，结合问题二的情况，假设观测站干流的污C染物浓度在一个月内保持不变，则C 0,又因为D 0,所以上式进一步简化t为:kC其中，U为断面水流平均流速，C为某组分子在x处的浓度，k为降解常数取0.2。由此解出长江干流污染物浓度 Cx与据观测站x米处的函数关系：Cieui10（2）对于长江上任一段干流 AB A为第i个观测点，B为第i 1个观测点， 两观测点的距离为dii 1，假设该段干流之间有n个排污口（包括支流入口和直排 口），为计算的方便，假设排污口聚集在一起，其总的排污口的流量为q，平均流量为u，污染物的浓度为c，qc即为AB段的排污量

18、,则CiqcUi 1CbqcQACAeqc C B Q BC AQAe#Ua Ub2对于任意一个江段AB污染物的浓度Ca,Cb，水流量Qa,Qb，流速Ua,Ub 距离dAB和降解系数k都已知，则可以算出AB段污染物的排放量Wj（3）根据所给数据可以算出，对于每个月每个江段的排污区间Wj(j=1,2,.,13个月），取13个月的平均值W11313Wj,表示平均排污量,j 1该江段每千米的排污量w芒为相对排污量，每个江段都有一个相对排污量,以此为指标可以判断长江干流的主要污染源。422模型二的求解按步骤运用MATLABS求解（程序见附录3）,结果如下表:表2： COOM的排污量江段s1 s2S2

19、S3S3 S4S4 S5S5 S6s6 SyW32510.5435266.774600229855.6211104.1531404.31相对排污量34.2245.33116.4659.71105.7560.05排序651423表3： NH3-N的排污量江段S1 S2S2 S3S3 S4S4 S5S5 S6S6 S7W2822.623191.464451.772559.311191.38158.46相对排污量2.974.1011.275.127.260.34排序5413264.2.3问题二的结果分析问题二用长江干流7个观测站点将长江分为6个江段，通过污染物流量随 时间的变化关系建立了微分方程模型

20、，先计算出每月每段的高锰酸盐和氨氮的 量，再对六个污染源近一年所排放的污染物的数量求期望，分析得出污染物高 锰酸盐指数和氨氮的污染源主要集中在宜昌至岳阳之间。4.3问题三的解答4.3.1模型三的建立与求解当系统内部信息和特性是部分已知的，另一部分是未知时，人们往往难以 建立客观的物理原形，内部因素难以辨识，以及相互之间的关系较为隐蔽，难已准确了解这类系统的行为特征，因此，对于这类问题进行定量描述，即建立 模型难度较大，便选择建立灰色系统，对问题求解。1.灰色动态模型附件4中给出了 1995-2004年长江各等级在丰水期，枯水期，水文年的百 分比，为了预测此趋势之下未来十年长江水质的变化，通过对

21、水文年全流域水 质的变化研究来对长江未来水质污染的发展趋势作出预测分析。建立灰色预测 模型，首先舍去异常数据：1995年IV类、劣V类数据，然后从附件4中提取出 水文年全流域随时间变化的原始数据。00,0 r0“x x 1,x2 .x n n 10对于原始数据进行累加处理，即k1 0x k x i k 1,2.10i 1加生成数列为x 1 x 1 1 , x 1 2 .x1 n n 10采用一阶单变量微分方程进行拟合，得到白化方程的GM( 1,1 )模型dx1ax 1 tudt其中式中的a,u为待定系数。综上所述，灰色动态模型为：z1 k 0.5x1 k 0.5x1 k 1 x0 k az1

22、k u2.模型求解取水文年全流域的饮用水预测为例，可通过对其变化研究对长江未来水质 污 染的发展趋势作出预测分析求解步骤：第一步：从附表中提取出数据x 093.1,85.3,80.7,88.4,80.2,74,73.7,76.7,77.5,68第二步：利用公式对x0做累加生成得93.1 178.4 259.1 347.5 427.7501.7 575.4 652.1 729.6 797.6第三步：构造矩阵B和数据向量,x°与x1满足关系丫 nbi?其中13#Yn0 0 Tx 2 .x n n 10T85.3 80.7 88.4 80.2 74 73.7 76.7 77.5 68135

23、.75 1218.75 1121x111Bx2211xn2x1 1 ,1x1 3 ,11 x1 n ,1303.3 1387.6 1464.7 1538.55 1613.75 1690.85 1690.85 1 763.6 1第四步：求解常系数Bt B 1 BTYn0.023188.8436第五步：得出表达式:、"11ua kU0.0231kY k 1X1e3753 e3846aaY1ki所以丫0表示还原后的值。(2) 结果求解按步骤运用MATLABS求解(程序见附录4),运用此程序根据过去十年的数据可以预测出未来10年的排污量和可饮用水如下表：年代20052006200720082

24、00920102011201220132014排污量303.01322.52343.29365.39388.92413.96440.61468.98499.18531.32表4:未来10年排污量年2020202020202020202()1份0506070809101112131可6968666463626058565/1饮.6.0.4.9.4.0.6.2.8.5水127572275表5 :未来10年可饮用水所占比例4.3.2问题三的结果分析从表5可以看出，可饮用水所占的比例越来越低，所以，可以看出水质越 来越差，从2012年起，可饮用水的百分比甚至低于 60%,如此严重的水污染，应 该引起足

25、够重视，采用合理的措施来减少对长江水的污染。4.4问题四的解答4.4.1模型的建立与求解(1) 长江干流的W类和V类水以及劣V类水的比例与每年长江总水流量X, t和总排污量x2 t有关，因此分别以W类和V类水所占的比例y,和劣V类水所占比例y2为响应变量，以历年来的x, t = 9205, 9513, 9171.26, 9513, 9924,8892.8,10210, 9980, 9405 x2 t = 174,179,183,207,234,220.5,256,270,285为解释变量，做出多元线性回归，得到W类和V类水的总和及劣V类水与长江总水流量及总排污量的关系，做出一般的多元线性回归模

26、型：% a/ t y2 a2 X1 tbx21b2x2 tC1C2运用matlab用最小二乘法求解得到回归系数a1,b1,G0.0032, 0.079, 15.9276a2,b2,c20.0077, 0.1485, 56.2805就求得了W类和V类水的总和及劣V类水与长江总水流量及总排污量的关系。问题四要解决问题是：未来10年内每年都要求长江干流的W类和V类水的 比例控制在20%以内，且没有劣V类水。则约束条件为：y10.0077x1 t 0.1485x2 t 56.2805 20y20.0032x1 t0.079x2 t 15.9276 0综上所述，所建立的模型为:16Z2 t 0.0225

27、,将t 11,12,13,14，,20带入乙t的表达式中得到未来10年每年预yi0.0077xi t 0.1485x2 t 56.2805 20s.t.y20.0077 x1 t0.1845x2 t cy0.0032xi t 0.079x2 t 15.9276 0对上述回归方程进行回归检验时，发现回归效果并不是太好，考虑到长江干流的W类和V类水的比例与江总水流量及总排污量不是直接的线性关 系，因此对模型进行改进。(2) 优化模型X2 t长江干流的W类和V类水以及劣 V类水的比例与每年总x1 t排污量和长江总水流量 的比值Z tX2 :相对排污量有关，因此分别以W类和V类水所占的比例yi和劣V类

28、水所占比例 y为解释变量，Z t为响应变量，做多元线性回归。对表四的数据进行观测，得到Zt与时间t成线性关系，因此列出如下线性回归方程：Zi t tZ2 tayi t by2 t c运用matlab用最小二乘法进行求解，得到回归系数分别为,0.0012,0.0165a,b,c 0.0002,0.0005,0.0185计算时剔除了 2003年这个突变点，且检验到回归性较好。得到的回归方程如下：乙 t 0.0012t 0.0165乙 t0.0002% t0.0005y2 t 0.0165将yi t20,y2 t 0带入Z2 t的表达式中，得到最大允许相对排污量测的相对排污量，则污水处理量 w t

29、max乙t Z2 t ,0 L ,因为每年的总水 流量(除特殊年份如1998年)变化不大，故取前10年的总水流量均值为每年的总水流量L 丄10 t10t，得到每年的污水处理量如下表:X1表6：未来10年预处理的排污量年代2005200620072008200920102011201220132014预处理排污量(亿吨)71.2483.1194.98106.86118.73130.60142.48154.35166.22178.094.4.2问题四的结果分析从表6可以看出，未来10年所需处理的污水越来越多，所以，可以看出水 质越来越差，如此严重的水污染，应该引起足够重视，采用合理的措施来减少 对

30、长江水的污染。4.5解决长江水质问题的建议通过问题一的计算可以看出，长江水质的枯水期、丰水期与平水期都由水 流量的不同导致水质情况也不同。从问题三、四的求解结果可以看出，在未来 的十年里，长江水质不断恶化，废水排放量逐年递增，需要处理的污水也越来 越多。我们认为可以从以下几个方面着手解决长江水质问题:根据模型假设，在不发生特大洪涝灾害时长江的水流量是一定的，总结得 出解决长江水质问题的关键，在于减少废水的排放量的同时增强废水的处理能 力。鉴于此，(1)加大处罚力度。(2) 对污染严重地区进行重点治理。(3) 加强宣传，增强民众保护长江的意识。严格进行检测，处理违规企业，严惩腐败。19（5）保护

31、长江立法，引起政府和人民的高度重视5模型的评价、改进与推广5.1模型的评价5.1.1模型的优点（1）问题一中建立了模糊综合评价模型，通过精确的数字手段处理模糊的 评价对象，对蕴藏信息呈现模糊性的资料作出了比较科学、合理、贴近实际的 量化评价；（2）问题二的模型，简化了许多因素，忽略了支流的影响，但这种理想化 还是合理的，较符合实际情况；（3）问题三采用灰色预测模型，具有少数据性、良好的时效性、较强的系 统性和关联性等特征，可以合理的对数据做出预测。5.1.2模型的缺点在求解模型时，为了使问题得到方便的解决，往往采用简化的手段进行求 解，因此求出的结果与真实值会存在一定的偏差。5.2模型的改进在

32、问题三中运用了灰色动态模型对未来10年长江的水质作出了预测，但并未对预测的结果进行检验，不知预测结果的准确度有多高，在此，最好能将模 型三的基础上对灰色模型的预测结果进行检验，以保证预测结果有较高的可信 度。5.3模型的推广模型一中，除了用模糊算法评价外，还可以用BP神经网络来求长江水质的在综合评价而；而在求解问题三和问题四时采用的灰色预测，还可以用回归分 析预测、时间序列预测和神经网络预测等代替。206参考文献11999宋晓秋编者，模糊数学原理与方法，北京，中国矿业大学出版社,2 张俊福等编著，应用模糊数学，北京，地质出版社，20013 岳超源，决策理论与方法，北京：科学出版社， 20034

33、 邓聚龙，灰理论基础，武汉：华中科技大学出版社，2002韩中庚，数学建模方法及其应用，北京：高等教育出版社，2005附录各观测点模糊关系矩阵:1000001000001号0.21610.216100002 号 0.04820.048200000.0940.09400000.51940.519400001000001000003号0.43750.437500004 号 0.89290.892900000.32660.326600000.90.900001000000.97030.970300005号0.21430.214300006 号 0.28750.287500000.02970.02970

34、0000.2280.22800000.99410.994100000 0.55860.55860007号0.04650.046500008 号 0 0.6215 0.62150001000000 0.8486 0.84860001000000.57670.57670009号0.36790.3679000010号 0.6697 0.66970000.80110.8011000000.62360.623600001000000.07330.0733000011号10000012号0.24290.2429000010000000.83280.83280001000 0 00.94760.947600

35、0013号00.09650.096500014 号 0.66250.662500000.67340.673400000.13570.1357000000.69820.698200010000015号0.16200.1620000016 号 0.87150.871500000000100.38970.3897000010000017 号0.51070.510700000.39170.39170000附录二问题一的程序(以17号观测站为例)：clc;cleara=10 0 0 00;0.51070.51070 0 0 0;0.39170.3917000 0；b=0.1450.2260.629;c=b*a;d=0;e=0;fo