长江水质评价与预测的数学模型摘要首先,利用附录3的数据，从时间上和空间上分析长江流域水质，得出长江水质的污染程度有所增加，但不明显.大部分污染比较严重的地区都位于支流上. 在求解主要污染物高锰酸盐指数和氨氮的污染源时先计算出单独一个河段内的排污量，进而求出一个河段内包括降解的污染物总量.计算出长江干流近一年多高锰酸盐指数和氨氮的主要污染源是位于重庆朱沱至湖北宜昌之间（每月排放约230万吨高锰酸盐和21万吨氨氮）、湖北宜昌和湖南岳阳之间（每月排放约206万吨高锰酸盐和20万吨氨氮）的工业带.然后,借鉴马氏链模型，结合过去十年水文年全流域数据，拟合出转移矩阵，预测得到未来十年各类水占河长的百分比(如2009年长江劣V类水占河段长约26.28%，2014年长江劣V类水河段占河段长约36.93%).如果未来10年内每年都要求长江干流的Ⅳ类和Ⅴ类水的比例控制在20%以内，且没有劣Ⅴ类水,2005年～2014年每年需要处理的废水分别为(单位：亿吨）：25.48，24.89，26.59，29.42，32.82，36.48，40.23，44.01，47.76，51.48.最后,针对目前长江水质污染状况，提出了切实可行的建议，具有较强的参考价值. 关键词:马氏链模型;水质评价;水质污染;数据拟合1 问题的提出水是人类赖以生存的资源，保护水资源就是保护我们自己，对于我国大江大河水资源的保护和治理应是重中之重.专家们呼吁：“以人为本，建设文明和谐社会，改善人与自然的环境，减少污染.”长江是我国第一、世界第三大河流，长江水质的污染程度日趋严重，已引起了相关政府部门和专家们的高度重视.2004年10月，由全国政协与中国发展研究院联合组成“保护长江万里行”考察团，从长江上游宜宾到下游上海，对沿线21个重点城市做了实地考察，揭示了一幅长江污染的真实画面，其污染程度让人触目惊心.为此，专家们提出“若不及时拯救，长江生态10年内将濒临崩溃”，并发出了“拿什么拯救癌变长江”的呼唤.研究下列问题：(1)对长江近两年多的水质情况做出定量的综合评价，并分析各地区水质的污染状况.(2)研究、分析长江干流近一年多主要污染物高锰酸盐指数和氨氮的污染源主要在哪些地区?(3)假如不采取更有效的治理措施，依照过去10年的主要统计数据，对长江未来水质污染的发展趋势做出预测分析，比如研究未来10年的情况.(4)根据预测分析，如果未来10年内每年都要求长江干流的Ⅳ类和Ⅴ类水的比例控制在20%以内，且没有劣Ⅴ类水,每年需要处理多少污水？(5)对解决长江水质污染问题切实可行的建议和意见.2 模型的假设2.1 长江干流的自然净化能力均匀；2.2水质只受高锰酸盐指数(CODMn)、氨氮(NH3-N)、溶解氧(DO) 、PH值的影响，与其它因素无关；2.3 污染源均匀分布在河岸两侧；2.4 不同水质的水随机分布在全流域上；2.5 两个监测站之间水流速均匀变化.3 符号的约定W表示第k监测站13个月的CODMn总量(包含上一个监测站对它的影响) 1kW表示第k监测站13个月NH3-N的总量(包含上一个监测站对它的影响) k2W表示第k监测站13个月CODMn的总量(不包含上一个监测站对它的影响) '1k'2k W 表示第k 监测站13个月NH3-N 的总量(不包含上一个监测站对它的影响) 0C 表示综合降解系数ijkρ 表示在第k 个监测站第i 个月第j 种污染物的浓度)(i Q o 第i 年长江废水排放总量i len 第i 个河段的长度ij v 第i 月第j 个监测点的平均水流速 ij M 第i 月第j 个监测点的平均水流量4 模型的建立与求解4.1 长江近两年多水质情况与各地区水质污染状况的定量综合评价 4.1.1 污染物含量的变化规律要看出近两年多长江水质的变化情况，我们首先要找出主要污染物的含量随时间的变化规律.由题目我们可以知道，影响水质的主要污染物有高锰酸盐指数(CODMn)、氨氮(NH3-N)、溶解氧(DO) 、PH 值.把长江全流域看成一个整体，),(j i W 为全流域第i 个月第j 种污染物的总量，ijP 为第i 个月第j 种污染物总量占28个月总污染物总量的百分比，即：)1,2,3,4j 28;,1,2,i ( ),(),(281===∑= i ij j i W j i W P ，其中∑==171*k ijk i ij V W ρ,)17,,2,1( =k ,i V 表示第i 个月的水流量，ijk ρ表示在第k 个监测站的第i 个月第j 种污染物的浓度.这里的ijP 都是百分率，数值上有良好的可比性，能很好的反映出三种污染物含量随着时间(月份)的变化规律.作出28个月长江流域四项监测项目总量所占百分比随时间的变化曲线图，如图1所示：由图1我们可以很直观地看出污染物含量的变化规律：（1） 污染物总量是以年为单位成周期性变化的，且相对总量逐年增加； （2） 污染物总量与水流量的变化趋势基本一致；由附录4我们可以知道，每年长江的枯水期为1月～4月，丰水期为5月～10月平水期为11月～12月，在丰水期，污染物的总量有明显的增长，随着枯水期的到来，污染物的总量也随之减少.（3） 四种污染物的变化趋势基本一致；（4） 主要污染物高锰酸盐和溶解氧含量都逐年增加，氨氮的含量有所下降. 4.1.2 长江整体水质情况及各地区水质的污染状况首先，把17个观测站按地理位置由西到东、由支流到干流的顺序重新排列(1、四川攀枝花龙洞 2、四川乐山岷江大桥 3、四川宜宾凉姜沟 4、四川泸州沱江二桥 5、重庆朱沱 6、湖北宜昌南津关 7、湖南长沙新港 8、湖南岳阳岳阳楼 9、湖南岳阳城陵矶 10、湖北丹江口胡家岭 11、湖北武汉宗关 12、江西南昌滁槎 13、江西九江蛤蟆石 14、江西九江河西水厂 15、安徽安庆皖河口 16、江苏南京林山 17、江苏扬州三江营)，根据附录3算出各观测站28个月四种污染物的平均浓度：28),,(),(281∑==i k j i k j ρρ，)28,,2,1;17,,2,1;4,3,2,1( ===i k j图1其中),,(k j i 表示第i 个月第k 个监测站第j 种污染物的浓度.作出四种主要监测项目(pH*、DO 、CODMn 、NH3-N)的平均浓度与各观测站之间的关系图，如图2所示：由图2我们可以很直观地看出各观测站4种污染物的浓度曲线，得污染物浓度的分布情况：（1） 支流比干流污染严重，大部分重污染地区都位于支流上； （2） 四川乐山岷江大桥的氨氮和高锰酸盐的含量都很高； （3） 江西南昌滁槎的氨氮含量很高；（4） 湖北宜昌南津关、湖南岳阳岳阳楼、湖北武汉宗关的高锰酸盐含量都很高. 由于重庆观测站朱沱位于川、渝省界，尚未进入重庆，未能反映重庆河段的污染情况.综上所述，长江水质的污染程度虽有所增加，但不明显，这是由于长江水量目前还比较大，掩盖了问题的严重性.4.2研究、分析长江干流近一年多主要污染物高锰酸盐指数和氨氮的主要污染源 4.2.1 流过各监测点的高锰酸盐和氨氮总量由附录3的数据可以知道，2004.04～2005.04长江干流主要观测站点的平均水流量及其所对应的高锰酸盐指数和氨氮浓度，可以求出在这13个月内，每个月流过干流各图2监测点的高锰酸盐和氨氮的总量：ikijk jk M W *ρ=，)7,,2,1( =k 所以13个月流过干流监测点污染物为：∑==131i ijkkj W Q )321(含量表示含量，表示N NH j CODMn j -==.得到流经各监测点的污染物量如下表（单位：吨/月）：流经4、5、6三个监测点的高锰酸盐和氨氮的总量最多，但这里没有消除上游河段的影响，jkW 是包含上游全河段污染物的总量，不利于我们分析主要污染物高锰酸盐和氨氮的污染源主要在哪些地区，因而求出各个河段单独的污染物排放量.以下我们考虑消除前面河段的影响.4.2.2 单独一个河段内的排污量由河流污染物一维稳定衰减规律的微分方程ρ0c dxdcV -=积分解得Vx c e 00-=ρρ,其中X 表示测试河段离污染源的距离，V 表示水流在该河段的平均速度，0ρ表示前一节点污染物的浓度，ρ表示前一节点对测试点的影响浓度.2004.04～2005.04长江干流7个观测站点平均水流量、不同时间的高锰酸盐指数和氨氮含量、不同时间的水流数度，可以求出在每一个月份内，在监测点上一河段内所排放的高锰酸盐和氨氮的总量： VCox i i i eW W W ---=1' ，把十三个月的总量相加：;2,17,,2,1,)(131==-=∑=-j k eW W W i VCox ikj ikj kj （程序代码见附录1）得到在长江干流不同河段近一年多主要污染物高锰酸盐和氨氮的排放总量（单位：吨/月）：从表中可知重庆朱沱到湖南岳阳城陵矶之间河段排放的高锰酸盐和氨氮总量最多，其中重庆朱沱到湖北宜昌南津关之间的河段每个月排放的高锰酸盐和和氨氮分别有1.2083×10^6吨、1.0933×10^5吨，湖北宜昌南津关到湖南岳阳城陵矶之间和的河段每个月排放的高锰酸盐和氨氮分别有1.5759×10^6吨、1.5251×10^5吨.干流近一年多主要污染物高锰酸盐指数和氨氮的污染源主要在重庆朱沱到湖南岳阳城陵矶之间河段. 4.2.3 每个河段的总排放量(包括降解的污染物)长江两岸工业分布呈工业带形态，这使对应的污染物排放量也呈带状，比如九江市岸边污染带约8公里长,重庆市在80年代就形成了81公里的岸边污染带，武汉几个大的排污口形成了几公里长的污染带.所以我们引进“单位长度江段污染物排放量”，即每一河段一公里排放的污染物量a .设第i 与第1+i 监测点之间的水流速是均匀变化的，得到单位长度江段水流速变化为i i i len v v )(1-=+α，第i 个河段内的有1-i len 个单位的排放量，第j 个单位河段排放的污染物经过降解流到下一个监测点“剩余”的污染物为αj v j C i ae +-0，1-i len 个单位排放量的“剩余”的总污染物为∑-=+-11i i len j j v Cojaeα.当a 表示高锰酸盐排放量时，这就是我们上面统计的'1k W ，有∑-=+-=111'i i len j j v Cojkk ea W α，解出各段对应得a 值，同理可以算出氨氮排放量b 值.如下所示(单位：吨/月)：1817.8 2994.6 5252.13558.3 4229.0 3545.9 157.80 270.96508.25305.02 310.03 17.89对应的河段排放总量(单位：10^5吨/月):重庆朱沱到湖南岳阳城陵矶之间河段排放的高锰酸盐和氨氮总量都是最多的，其中重庆朱沱到湖北宜昌南津关之间的河段每个月排放的高锰酸盐和和氨氮分别有2.3268×10^6吨、2.1054×10^5吨，湖北宜昌南津关到湖南岳阳城陵矶之间河段每个月排放的高锰酸盐和氨氮分别有2.0693×10^6吨、2.0025×10^5吨.攀枝花重庆朱沱江西九江安徽安庆江苏南京湖南岳阳湖北宜昌a 值：b 值：综合5.2.2得出的数据，我们可以确定干流近一年多主要污染物高锰酸盐指数和氨氮的污染源主要集中在在重庆朱沱到湖南岳阳城陵矶之间的河段. 4.3长江未来10年水质污染的发展趋势 4.3.1 模型的准备马氏链模型及其基本方程 按照系统的发展，时间离散化为 ,2,1,0=n ，对每个n ，系统的状态用随机变量n X 表示，设n X 可以取k 个离散值k X n ,,2,1 =，且记)()(i X P n a n i ==，即状态概率，从i X n =到j X n =+1的概率记)|(1i X j X P p n n ij ===+，即转移概率.如果1+n X 的取值只取决于n X 的取值及转移概率，而与 ,,21--n n X X 的取值无关，那么这种离散状态按照离散时间的随机转移过程称为马氏链.由状态转移的无后效性和全概率公式可以写出马氏链的基本方程为∑==+kj jij i p n a n a 1)()1( k i ,,2,1 =并且)(n a i 和ijp 应满足：()∑==ki i n a 11， ,2,1,0=n ∞≥ij p ， k j i ,,2,1, =∑==kj ijp11， k i ,,2,1 =4.3.2 水质污染的发展趋势根据水环境质量的不同，把水环境分为6类(Ⅰ类、Ⅱ类、Ⅲ类、Ⅳ类、Ⅴ类、劣Ⅴ类)，其中Ⅰ、Ⅱ、Ⅲ类为可饮用水.各类水质间可以相互转化，且下一个状态只取决于上一状态，与以前的状态无关，这正符合马氏链模型.把6类水质分为四种状态：Ⅰ、Ⅱ类同属于状态1；Ⅲ类属于状态2；Ⅳ、Ⅴ类同属于状态3；劣Ⅴ类属于状态4.用随机变量n X 表示第n 年的状态，1=n X 表示水质属于Ⅰ类或Ⅱ类，2=n X 表示水质属于Ⅲ类，3=n X 表示水质属于Ⅳ类或Ⅴ类，4=n X 表示水质属于劣Ⅴ类，n=0,1,2,….用)(n a i 表示第n 年处于状态i 的概率，i=1,2,3，4，即)()(i X P n a n i ==.用ijp 表示已知今年处于状态i ，来年处于状态j 的概率，i,j=1,2,3，4，即)|(1i X j X P p n n ij ===+.第n+1年的状态1+n X 只取决于第n 年的状态，n X 和转移概率ijp ，而与以前的状态 ,,21--n n X X 无关，即状态转移具有无后效性.第n+1年的状态概率可由概率公式得到:⎪⎪⎩⎪⎪⎨⎧+++=++++=++++=++++=+4443432421414434333232131342432322212124143132121111)()()()()1()()()()()1()()()()()1()()()()()1(p n a p n a p n a p n a n a p n a p n a p n a p n a n a p n a p n a p n a p n a n a p n a p n a p n a p n a n a这里的)(n a i 为第n 年处于第i 种状态的水域长度占河长的百分比.对于不同的转移矩阵P ，来年有不同的百分比，结合附录4十年不同水质变化数据，可以建立以下模型：min∑=9121i ix+∑=9122i ix+∑=9123i ix+∑=9124i ix, 其中41431321211111)()()()()1(p i a p i a p i a p i a i a x i ----+=； 42432322212122)()()()()1(p i a p i a p i a p i a i a x i ----+=； 43433323213133)()()()()1(p i a p i a p i a p i a i a x i ----+=； 44434324214144)()()()()1(p i a p i a p i a p i a i a x i ----+=；s.t.∑==411j ijp, i=1,2,3,4； ……(1) 10≤≤ij p ； (2)由模型的实际意义，不同状态水质之间的转化应满足关系： (1)向自身状态转化的概率不小于向其它状态转化的概率； (2)向相邻状态转化的概率不小于向相隔状态转化的概率；0;p -p 0;p -p 0;p -p 141313121211≥≥≥ (3)0;p -p 0;p -p 0;p -p 0;p -p 2423242123222122≥≥≥≥ (4)0;p -p 0;p -p 0;p -p 0;p -p 3433323331323134≥≥≥≥ (5)0;p -p 0;p -p 0;p -p 434442434142≥≥≥ (6)对于不同的水期，不同的流域可以求出对应的转移矩阵P ，然后再递推出长江未来十年各类水质的变化情况.现对水文年的全流域十年数据进行分析，用数学软件lingo8进行求解，得到转移矩阵⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=1000215.0355.0215.0215.00106.0449.0445.00098.0451.0451.0P ,其对应可以求出长江未来十年水文年的各类水质;20...11 ,*)1()(=-=i P i a i a 有：用枯水期全流域的数据求出干流未来十年水质的变化情况：用丰水期全流域的数据求出干流未来十年水质的变化情况：用水文年干流和支流的数据可以求出各自未来十年水质的变化情况（见附录2、3）. 通过对5个表格的分析对比，我们可以知道长江水质变化呈现统一趋势，劣Ⅴ类水河长百分比逐年增加，支流污染比干流污染严重(如水文年劣Ⅴ类水十年变化：干流10.5％～23.6％，支流15.4％～37.6％)，枯水期劣Ⅴ类水河长的变化比丰水期快(枯水期17％～45％，丰水期12.7％～28.8％).在不采取有效治理措施的情况下，劣Ⅴ类水河长百分比会以一个很快的增长速度不断增大，致使长江陷入深度危机，若不及时拯救，10年之内，长江水系生态将濒临崩溃.4.4在满足要求情况下的污水处理方法第i 年排放的污染物为 )(i Q o ，水文年干流第i 年四个状态水量百分比为)(4),(3),(2),(1i a i a i a i a ，有关系e i a d i a c i a b i a a i Q o ++++=)(4*)(3*)(2*)(1*)( 10,,2,1 =i ，其中e d c b a ,,,,为常系数.对过去十年数据，用最小二乘法拟和得到)(i Q o 和)(4),(3),(2),(1i a i a i a i a 的关系为： 104.312;+a4(i)*3.983+a3(i)*2.706+a2(i)*1.612=(i)Q 0用5.3计算得到的未来十个水文年干流四个状态水量百分比，而我们的废水处理目标要求要把0a4(i)=而且20a3(i)≤，即312.1042.0*706.2)(2*1612.1)'(0++=i a i Q 这样可以得到每年要处理的污水量为：)'()(00i Q i Q -，即（单位：亿吨）：25.48，24.89，26.59，29.42，32.82，36.48，40.23，44.01，47.76，51.48 .4.5对解决长江水质污染问题的建议和意见意见：目前长江造成局部污染严重而总体水质良好的原因不在水量和水源缺乏，而在水源浪费，废水处理率低甚至出现“直排”现象，人们环保意识薄弱.建议：（1）抓好对人民环保意识的教育，进一步提高水资源保护意识，增强对水资源保护重要性、紧迫性的认识唤醒民众的环保意识.（2）强化流域水资源保护机构，制定流域管理法规，统一规划和监督管理水环境质量，防治水环境污染.（3）依据水环境功能用途要求和水体稀释自净能力，建立不同类型的保护区，优先保护好生活用水水源，使其不受污染.（4）积极引进和开发无废或少废，不用水或少用水的工业技术，研究适合流域工业特色和自然环境特点的废污水处理利用和资源化技术，加快建设城市污水处理设施和资源化工程，降低产业水耗，提高废水利用率，使有限的水资源发挥更大的经济、社会和环境效益.（5）以水环境质量目标和污染物总量控制目标为导向，合理规划工业布局和调整工业结构，严格执行环保“三同时”制度和乡镇企业污染防治法规，使新增废污水及污染物排放量得到有效控制.（6）从生态环境和社会发展对水环境的需求出发，确定社会和环境可承受的水资源开发程度，合理规划水利工程布局、调节调度水资源水量，使水资源的环境功能得以充分的发挥和利用.5 模型的评价与推广5.1 模型的评价5.1.1 模型的优点：（1）模型能抓住影响水质的主要因素（高锰酸盐指数(CODMn)、氨氮(NH3-N)、溶解氧(DO)），能正确预测长江水质变化，过程清晰明了，结果科学合理;（2）模型具有较好的通用性，实用性强，对现实有很强的指导意义;（3）在求解主要污染物高锰酸盐指数和氨氮的污染源时先计算出单独一个河段内的排污量，进而求出一个河段内包括降解的污染物总量.5.1.2 模型的不足以及需要改进的地方：（1）监测点太少，不能全面反映长江流域水质;（2）没有考虑长江沿岸地形对流速、流量的影响.5.2 模型的推广我们建立模型的方法和思想对其他类似的问题也很适用，适用于其它河流、湖泊水质的评价和预测.参考文献[1] 姜启源等,数学模型[M],北京：高等教育出版社，2004.[2] 李强等,Maple基础应用教程[M],北京：中国水利水电出版社,2004.[3] 宋兆基等，MATLAB6.5在科学计算中的应用[M]：清华大学出版社.2005.[4] 总参谋部测绘局编制，中华人民共和国地图集[M]：星球地图出版社.2000.[5] 董哲仁等.中国江河1000问[M]：黄河水利出版社.2001.附录附录1单独一个河段内的排污量da 是221*6 的矩阵，记录了题目所给附录3的04年4月～05年4月的数据v 是13*7矩阵，记录了题目所给附录3的水流量数据len=[ 0 950 1728 2123 2623 2787 3251];w_Mn=zeros(1,7); % 临时记录7个监测点的高锰酸盐含量w_NH=zeros(1,7); % 临时记录7个监测点的氨氮含量pure1=zeros(1,7); % 记录7个监测点13个月的高锰酸盐含量pure2=zeros(1,7); % 记录7个监测点13个月的氨氮含量pure1(1)=116820 ; pure2(1) =6261; % 四川攀枝花监测点为第一个河段的上点flag=0; flag2=0;for i=1:221temp=mod(i,17);if temp~=0 && temp<=7t=temp; flag=flag+1;m=mod(flag,7);if m==0m=7;endn=(flag-m)/7;w_Mn(t)=da(i,3)*V(n,m);w_NH(t)=da(i,4)*V(n,m);endif mod(i,17)==0flag2=flag2+1;w1=w_Mn;w2=w_NH;vv=v(flag2,:);for j=7:-1:2 ％第j个监测点每一个月的污染物量Co=0.2; ff=1; ％减去了上一河段的影响time=2*(len(j)-len(j-1))/((vv(j)+vv(j-1))*3.6*24);wo=w1(j-1); temp1=wo*exp(-Co*time);wo2=w2(j-1); temp2=wo2*exp(-Co*time);w1(j)=w1(j)-temp1; w2(j)=w2(j)-temp2;pure1(j)=pure1(j)+w1(j);pure2(j)=pure2(j)+w2(j);endendendplot(pure1,'-p')hold onplot(pure2,'-*')grid on附录2用水文年干流的数据可以求出干流未来十年水质的变化情况：附录3用水文年支流的数据可以求出支流未来十年水质的变化情况：。
