B 水库排污问题的数学模型(3).doc

数学建模 学号：3110801236 姓名：诸震亚 班级：数学112 11级数学数学建模与数学实验(实践)任务书一、设计目的通过数学建模与数学实验(实践)实践环节，掌握本门课程的众多数学建模方法和原理，并通过编写C语言或matlab程序，掌握各种基本算法在计算机中的具体表达方法，并逐一了解它们的优劣、稳定性以及收敛性。在熟练掌握C语言或matlab语言编程的基础上，编写算法和稳定性均佳、通用性强、可读性好，输入输出方便的程序，以解决实际中的一些科学计算问题。二、设计教学内容1线性规划（掌握线性规划的模型、算法以及Matlab 实现）。整数线性规划（掌握整数线性规划形式和解法）。2微分方程建模（掌握根据规律建立微分方程模型及解法；微分方程模型的Matlab实现）。3最短路问题（掌握最短路问题及算法，了解利用最短路问题解决实际问题）。行遍性问题（了解行遍性问题，掌握其TSP算法）。4回归分析（掌握一元线性回归和多元线性回归，掌握回归的Matlab实现）。5计算机模拟（掌握Monte-carlo方法、了解随机数的产生；能够用Monte-carlo解决实际问题）。6插值与拟合（了解数据拟合基本原理，掌握用利用Matlab工具箱解决曲线拟合问题）。 摘 要 本文针对水库突发性事故排污问题，首先通过建立二维水质污染物浓度模型，给出了单个水库对干流造成大面积污染的可能性；然后建立两水库排污模型，分析了在另一水库有连续点源污染物排放及水流相互影响的情况下，两水库对干流造成大面积污染的可能性大小；并进一步针对第三种情况的发生，给出在短时间内控制污染的有效措施；且讨论了若污染物具有挥发性，上述各情况造成干流发生大面积污染的可能性大小，为水库事故性排污问题提供了有价值的理论依据。关键词：水库排污；污染物浓度；流量；水流速度一、问题的提出 近年来水库污染问题日益严重，某条江流上有2条支流，每条支流上都兴建了规模相当的水库。由于正处在雨水多发季节，因此两个水库都以一定规模的流量进行泄洪。某天晚上10：00，在其中的一个水库中发生了两船相撞的事故，而其中的一条船装载的p吨化学物质（这里的化学物质可以是具有挥发性的，也可能是极难挥发的）全部泄漏至水库中。当水上航运事故处置中心接获事故报告，立即要求该水库关闭水库泄洪闸，以免化学物质随洪水流入干流，发生更大规模的污染。水库闸门开始关闭时，已经处在事故发生后的1个小时，而水库闸门彻底关闭也需要1个小时的时间。 根据当地环境监测的有关规定，干流大面积污染的危险警戒值设为：三小时内q吨该化学物质发生泄漏。(1) 试建立合理的数学模型，讨论由于此次事故的发生，干流发生大面积污染的可能性；(2) 如果在另外的一水库中有一化工厂违规排放废料。废料中同样含有该化学物质。该工厂为躲避环境监测站的监控，均在晚上9：00-12：00违规进行周期性排放。在这种情形下，讨论由于此次事故的发生，干流发生大面积污染的可能性；(3) 如果以上两个水库间有一条人工修建的水渠相连接，水渠中的水流流向不定，但保证两水库之间的水流能够相互影响。那么上述结果是否会改变？请给出说明，若有改变，则给出修正的模型及结果；(4) 针对第三种情况，试给出在短时间内控制污染模型。二、问题的分析 由题目的背景知道，事故发生时两水库都正在泄洪，因此此时水库中的水流速水库较快。而泄露到水库中的化学物质不论是具有挥发性的，还是急难挥发的，它们对干流污染的情况总是类同的，因此我们总可以认为污染物是易溶急难挥发性物质。为使我们的模型简单，我们可以先假设事故发生在水库1中，污染物在水库中的分布是符合零维迁移模型的，此时流入水库的污染物能以很快的速度与水库中的水均匀混合，水库中任何部份水体的污染状况都是一样的，污染程度与水体在水库中的位置无关。而实际上，污染物在水中达到分布均匀是有一个迁移过程，符合污染物一维迁移方程，但在这样一个突发事件要求短时间内得到控制的问题中，我们总可以用污染物瞬时混合均匀状态模型来代替污染物一维迁移的过程。在解答问题一时，我们只要考虑在事故发生到关闭水库的两个小时内，流出水库的污染物的质量小于吨即可。问题二中提到的情况只是将干流的污染源从一个增加到两个，那么发生大面积污染的可能性就要增大。由于两个水库之间没有联系，我们只需要单独考虑2库的排污情况，然后加上1水库的污染物排放情况，最后综合考虑两个水库所排污染物的总量对干流的影响就可以了。问题三中建立人工水渠就是在问题二的基础上使水库1和水库2发生联系。由于水都是从高水位流向低水位，因此，我们只需考虑从1水库向2水库的流入情况，而不必考虑从2水库向1水库的流入情况。 分析问题四，干流已经发生大面积污染，对比水体污染物处理的各种手段，不论是化学、物理还是生物手段都不可能在短时间内除去污染物，因此我们只能通过稀释原理，建立污染物在干流中迁移迁移模型，使干流水体计算体积元内该污染物的增量为负值时，从而使得在干流水体污染物的浓度低于危险警戒值时的浓度，才能在短时间内达到控制污染。 三、模型的假设1 事故发生时，化学物质在短时间内全部泄漏至水库中，不考虑泄漏的速度和时间2水库水流平稳，水流为匀速，当水库关闭水闸时，由于需要1小时的时间，闸门是缓慢关闭的，对水流的影响不大，水流速度不变3不具挥发性的化学物质全部溶解在河水中，不考虑重金属或不溶性物质出现沉淀的情况4具有挥发性的化学物质溶解在河水中，其挥发情况是均匀、恒定的，即在两小时的时间内，其挥发速率恒定不变 四、符号约定 流量 第个水库的闸门个数 时间 水流速度 第个水库的闸门宽度 第个水库的闸门高度 闸门关闭速度 排污点对污染带内点处浓度贡献值 污染物总质量 敏感点到排污点的纵向距离 敏感点到排染点的横行距离 污染物降解系数 排污口污染物排放速率 排放到干流的污染物质量 五、模型的建立和求解5.1极难挥发的化学物质污染模型 干流大面积污染的危险警戒值为三小时内吨该化学物质发生泄漏，从事故发生到闸门完全关闭共有两小时的时间，则干流大面积污染的危险警戒值可相应转化为：两小时内吨该化学物资发生泄漏。对河流中污染物质量的计算遵循以下计算法则：污染物质量水库总流量污染物浓度5.1.1水库总流量计算 在事故发生后的一小时内，闸门并未开始关闭，在这一小时内，水流速度不变，水库任意时刻的流量流速闸门总的排水面积，水库情况如图所示：VbN (1)在事故发生一小时后，闸门开始关闭，关闭的时间需要一小时，关闭的速度较慢，对水流的影响不大，而闸门的排水面积在逐渐减小，闸门是以恒定速度关闭的，则闸门在时刻的排水面积为：一小时后，水库任意时刻的流量为： (2)从而可推出从事故发生到闸门关闭的两小时内，水库的总流量为与的积分和： (3)5.1.2二维水质浓度模型 当吨化学物质排放到河水中，由于河水的冲释、河床结构、流速等因素的影响，在不同的河段、时刻，污染物的浓度均不相同，计算污染物的浓度是本模型的难点所在。70年前美国Streeter和Phelps的河流一维水质模型，经过后人不断的改进和完善，现已发展起适用于河流、海湾、水库和湖泊的较完善的水质模型，但这些水质模型大多只适用于污染浓度控制计算。张玉清1对适用于污染浓度定量计算的水质模型结构进行了改进，建立了适用于总量控制的二维水质模型，能满足水污染浓度的定量计算。而干流中的污染物总质量均经水库闸口排出，因此在此仅计算水库闸口处排出的污染物总质量即可。在参考文献1中，建立了二维水污染浓度计算模型：方程形式 解析解 (4)为河流在横断面处的平均水深，为河流在横断面处的平均水深0VbN事故此模型可引用到小江小河江心事故性排放的浓度场预测，为污染物降解系数，由于从事故发生到污染物浓度场的预测时间较短，水环境对污染物的降解作用不明显，故水库各点的污染物浓度变化是连续的，由上图所示的水库事故发生后，各点的数据情况可得，在时刻经过水库闸门处这一横断面的污染物浓度为从点到点的积分：污染物质量水库流量污染物浓度由上式可得在时刻闸门处的污染物质量为：前一小时： (5)后一小时： (6)为污染物从排放到水中到水库闸口处开始有污染物的时间 ： 发生大面积污染的可能性 (7)5.2两水库排污模型若水库在晚上9：0012：00排放该化学物质，此时干流发生大面积污染的可能性又如何呢？事故在晚上10：00发生，到水库闸门完全关闭是晚上12：00，在此过程中，水库从9：00开始排污，讨论干流发生大面积污染的时间应改为9：0012：00，若在这3小时内，两水库排放到干流的污染物质量达到吨，则该干流发生大面积污染。水库属于污染物岸边点源连续排放，根据参考文献1中建立的适用于小河岸边点源连续排放浓度场模型为：方程形式 (8)解析解 (9)在3小时内，由于环境的降解作用不明显，故值为0水库并没有关闸，其流速恒定不变，根据模型5.1.1在任意时刻其流量为：则水库排放到干流的污染物质量为： (10)水库的排污时间是从9：00开始，所以对干流的影响时间为3小时，在此情况下干流中的污染物总质量为水库与水库的排污量的总和，干流发生大面积污染的可能性为：发生大面积污染的可能性5.3两水库相互影响模型如果两水库间有一条人工水渠相连接，水库间的水流能够相互影响，上述结果是否会改变？当两水库间有水渠相连，根据扩散原理，分子必定会从浓度高的水库向浓度低的水库扩散，这样一部分的污染物将扩散到水库中，再经过水库冲释到干流中，一方面水库排放到干流中的污染物将减少，而另一方面虽然经水渠排放到水库中的污染物最终也会排放到干流中，但经过此过程，污染物排放到干流的浓度和时间都会减缓。水库水库排污口水渠干流相对于水库来说，水渠的宽度很窄，水流的稳定性较高，水渠内各点的浓度差不明显，当污染物通过水渠扩散到水库后，水渠内各点的浓度是相同的，等于水渠在水库的入口处的浓度。两水库及干流的情况如下图所示： 根据二维水污染浓度模型，由(4)式得，水渠入口处的污染物浓度为： (11)由水渠相连的两个水库，结构原理与连通型管类似，在水渠的作用下，两水库的水流会达到一个平衡状态，假设在水库没有关闸时，两水库及水渠的流速相同，为，当闸门开始关闭时，水位开始升高，根据型管的原理，此时水库总流量的变化值将有通过水渠排放到水库，以达到两水库平衡，则水库在没有水渠的情况下时刻流量的变化值为：在有水渠的情况下，将有的流量变化值排放到水库，则水渠在时刻的流量为：水库在两小时内排放到干流的污染物质量计算如下:前一小时 (12)后一小时 (13)由于水坝的作用，且根据型管原理，水库的水位虽然升高了，但水流在水坝的作用下，水流变化不明显，此时水渠对于水库来说，相当于一个排污口，根据点源排放浓度场模型，它对水库闸口处的浓度贡献值为： (14)水库的污染物浓度值为排污口与水渠两者的叠加，但排污口对干流的作用时间为3小时，而水渠对干流的作用时间为2小时，水库排放到干流的污染物质量计算如下： (15)此时，干流发生大面积污染的可能性5.4控制污染模型 针对第三种情况，要在短时间内控制污染物，则水库也必须关闭，以防止过多的污染物排放到干流中，根据题意，水库要在事故发生1小时后才被通知关闭闸门，那么水库和水库的情况相同，也要在事故发生后1小时才开始关闭闸门，假设两水库闸门关闭的情况相同，闸门关闭的速率同为，水库向干流排放的污染物情况不变，仍用(12)、(13)式计算，水库在事故后1小时关闭，即9：0010：00时段只有一个排污口排放污染物，10：0011：00时段有排污口和水渠排放污染物，闸门没开始关闭，11：0012：00有两个排污处，但闸门以开始慢慢关闭。在二维水污染浓度计算模型和点源连续排放浓度场模型中，由于时间较短，水环境的降解作用不明显，故降解系数，但是我们在控制污染物的过程中，可以通过加入一些针对该种污染物的降解物质，使降解系数增大，从而降低污染物对干流造成大面积的污染。例如：若此污染物为油质性污染物，则可往水库中抛入稻草、布碎，这两种物质对油质性污染物有很高的吸附能力，能够有效地降低污染。对水库的污染物计算模型如下：9：0010：00时段 (16)10：0011：00时段 (17)11：0012：00时段闸门匀速关闭，水流量改变 (18) (19)在(16)、(17)、(19)式中 排放到干流的污染物总量通过控制污染模型，在两水库的闸门开始关闭后，干流污染物的质量将比无控制污染时减少，对模型中的变量进行赋值，作图象(1)，对比控制污染前后，干流污染物的总质量变化情况如下：图(1)由图(1)可以看出，污染控制模型对排放到干流的污染物质量有明显的减弱作用，可见措施的有效性，根据此结果，可运用到实际生活中，详见模型的推广。5.5具有挥发性的化学物质污染模型 该化学物质挥发的速率是恒定的，二维水污染浓度模型可改写为： (20)被排放到干流中的污染物的总质量计算方法仍按(5)、(6)式计算，只需将(20)式代入到(5)、(6)中即可对于水库点源连续排放浓度场模型，由于该化学物质具有挥发性，(9)式改写为： (21)把(21)式代入到(10)式中，可计算出该化学污染物由水库排放到干流的质量，干流的污染物总质量仍由两水库排放的污染物质量的总和。对于两水库相互影响模型，(11)式改写如下： (22)将(20)、(21)、(22)式代入到(12)、(13)、(14)、(15)、(16)、(17)，(19)式中，同理可讨论出在两水库相互影响下，可挥发化学物质使干流发生大面积污染的可能性，以及控制污染模型。 六、模型的推广与应用通过模型5.3和模型5.4可以看出，在两水库间修建人工水渠，使两水库的水流可以相互影响，当突发性事故污染发生时，两水库的水流能够有效地稀释污染物，对突发性事故污染问题有明显的减弱作用，图(1)证明了此结果的正确性及可行性。根据此模型的结论，在实际生活中，我们可以在两个或三个地理条件合适的水库间修建水渠，使彼此独立的水环境相互影响，从而降低污染物对干流造成大面积污染的可能性。另外，针对突发性的事故污染，还可以在污染源附近加入相应的污染物降解物质，使降解系数增大，从而降低污染物对干流造成大面积的污染。该加入何种降解物质，则要根据具体的污染物及相应的化学原理来决定。 参考文献：1 张玉清著，河流功能区水污染物容量总量控制的原理和方法，北京：中国环境科学出版社， 118132，2001。2 姜启源，数学模型(第三版)M，北京：高等教育出版社，2003。3 林芳荣等编译，面污染源管理与控制手册，广州：科学普及出版社广州分社， 1543，1987，1，6。4 方子云主编，水资源保护工作工作手册，河海大学出版社，1998，7。附录遗传算法程序function Geneticmaxgen=200; sizepop=100; fselect=tournament; fcode=float; pcross=0.6; fcross=float; pmutation=0.2; fmutation=float; lenchrom=1 ; bound=0 100000;individuals=struct(fitness,zeros(1,sizepop)； avgfitness=; bestfitness=; bestchrom=; for i=1:sizepop individuals.chrom(i,:)=Code(lenchrom,fcode,bound); x=Decode(lenchrom,bound,individuals.chrom(i,:),fcode); individuals.fitness(i)=Aimfunc(x);endbestfitness bestindex=min(individuals.fitness);bestchrom=individuals.chrom(bestindex,:);avgfitness=sum(individuals.fitness)/sizepop;trace=avgfitness bestfitness; for i=1:maxgen individuals=Select(individuals,sizepop,fselect); avgfitness=sum(individuals.fitness)/sizepop; individuals.chrom=Cross(pcross,lenchrom,individuals.chrom,. sizepop,fcross,i maxgen); individuals.chrom=Mutation(pmutation,lenchrom,individuals.chrom,. sizepop,fmutation,i maxgen,bound); for j=1:sizepop x=Decode(lenchrom,bound,individuals.chrom(j,:),fcode); individuals.fitness(j)=Aimfunc(x); end newbestfitness,newbestindex=min(individuals.fitness); worestfitness,worestindex=max(individuals.fitness); if bestfitnessnewbestfitness bestfitness=newbestfitness; bestchrom=individuals.chrom(newbestindex,:); end individuals.chrom(worestindex,:)=bestchrom; individuals.fitness(worestindex)=bestfitness; avgfitness=sum(individuals.fitness)/sizepop; trace=trace;avgfitness bestfitness; end if ishandle(hfig) figure(hfig);else hfig=figure(Tag,trace);endfigure(hfig);r c=size(trace);plot(1:r,trace(:,1),r-,1:r,trace(:,2),b-);title(适应度曲线 终止代数 num2str(maxgen);xlabel(进化代数);ylabel(适应度);legend(平均适应度,最佳适应度);disp(适应度 变量);x=Decode(lenchrom,bound,bestchrom,fcode);% show in command windowvpa(bestfitness,10);vpa(x,10);disp(bestfitness x);function ret=AimFunc(x)z=quadl(f,0.0001,7200,x);ret=abs(z-49.771090910000000000000000000000);function y=f(t,x)y=7565./sqrt(4.8*pi*t).*exp(-(x-3*t).(2)./(4.8*t); pick=rand(1,sum(lenchrom); bits=ceil(pick-0.5); temp=sum(lenchrom)-1:-1:0; ret=sum(bits.*(2.temp);case grey pick=rand(1,sum(lenchrom); bits=ceil(pick-0.5); greybits=bits; for i=2:length(greybits) greybits(i)=bitxor(bits(i-1),bits(i); end temp=sum(lenchrom)-1:-1:0; ret=sum(greybits.*(2.temp); case float % float coding pick=rand(1,length(lenchrom); ret=bound(:,1)+(bound(:,2)-bound(:,1).*pick;endfunction ret=Cross(pcross,lenchrom,chrom,sizepop,opts,pop) input : probability of crossover input : length of a chromosome input : set of all chromosomes input : size of population input : tag for choosing method of crossover input : current serial number of generation and maximum gemeration output : new set of chromosomeswitch opts case simple for i=1:sizepop pick=rand(1,2); index=ceil(pick.*sizepop); while prod(pick)=0 | index(1)=index(2) pick=rand(1,2); index=ceil(pick.*sizepop); end pick=rand; if pickpcross continue; end pick=rand; while pick=0 pick=rand; end pos=ceil(pick.*sum(lenchrom); tail1=bitand(chrom(index(1),2.pos-1); tail2=bitand(chrom(index(2),2.pos-1); chrom(index(1)=chrom(index(1)-tail1+tail2; chrom(index(2)=chrom(index(2)-tail2+tail1; end ret=chrom; case uniform for i=1:sizepop pick=rand(1,2); while prod(pick)=0 pick=rand(1,2); end index=ceil(pick.*sizepop); pick=rand; while pick=0 pick=rand; end if pickpcross continue; end pick=rand; while pick=0 pick=rand; end mask=2ceil(pick*sum(lenchrom); chrom1=chrom(index(1); chrom2=chrom(index(2); for j=1:sum(lenchrom) v=bitget(mask,j); if v=1 chrom1=bitset(chrom1,. j,bitget(chrom(index(2),j); chrom2=bitset(chrom2,. j,bitget(chrom(index(1),j); end end chrom(index(1)=chrom1; chrom(index(2)=chrom2; end ret=chrom; case float for i=1:sizepop pick=rand(1,2); while prod(pick)=0 pick=rand(1,2); end index=ceil(pick.*sizepop); pick=rand; while pick=0 pick=rand; end if pickpcross continue; end pick=rand; while pick=0 pick=rand; end pos=ceil(pick.*sum(lenchrom); pick=rand; v1=chrom(index(1),pos); v2=chrom(index(2),pos); chrom(index(1),pos)=pick*v2+(1-pick)*v1; chrom(index(2),pos)=pick*v1+(1-pick)*v2; end ret=chrom;endfunction ret=Decode(lenchrom,bound,code,opts) input : length of chromosome input : tag of coding method input : boundary of varibles output: value of variblesswitch opts for i=length(lenchrom):-1:1 data(i)=bitand(code,2lenchrom(i)-1); code=(code-data(i)/(2lenchrom(i); end ret=bound(:,1)+data./(2.lenchrom-1).*(bound(:,2)-bound(:,1);case grey for i=sum(lenchrom):-1:2 code=bitset(code,i-1,bitxor(bitget(code,i),bitget(code,i-1); end for i=length(lenchrom):-1:1 data(i)=bitand(code,2lenchrom(i)-1); code=(code-data(i)/(2lenchrom(i); end ret=bound(:,1)+data./(2.lenchrom-1).*(bound(:,2)-bound(:,1); case float ret=code;endfunction y=f1(t)y=0.58123108938362951707736011098632*erfc(30000-30*t)./sqrt(4.8*t); pick=rand; while pick=0 pick=rand; end index=ceil(pick*sizepop); pick=rand; if pickpmutation continue; end pick=rand; while pick=0 pick=rand; end pos=ceil(pick*sum(lenchrom); v=bitget(chrom(index),pos); v=v; chrom