数学建模-嘉陵江铊污染事件解题全过程文档及程序
嘉陵江铊污染事件
原题再现：
某年5月5日18时，四川广元市环境监测中心站监测发现嘉陵江入川断面水质异常，西湾水厂水源地水质铊元素超标4.6倍。初步判定污染源为川陕界上游输入型、一次性污染团。此次事件再一次警示民众水污染监控及预警的重要性。水污染可主要分为"点源污染"与"面源污染"两种类型，当出现水质重金属污染超标时，管理部门需要能尽早快速诊断出污染类型及污染源位置，这样能尽快切断污染源，减少污染对环境、民众工作生活等带来的损失。根据以上信息，请利用数学模型回答以下问题：
问题1：假设当出现水质铊元素点源污染或面源污染时，建立模型分析下游水质铊元素浓度变化趋势，由此估计下游的污染重点监控区域及时间。
问题2：当在下游水质监测点检测到铊元素超标及相关数据，能否利用数学模型快速估计出上游污染类型及污染源大致位置？
问题3：收集本次广元铊元素事件的相关数据，分析本次铊污染的主要特点。
问题4：根据问题3数据，进一步回答问题1和问题2，并由此验证问题1和问题2的合理性。
整体求解过程概述(摘要)
河流水质不仅影响周边环境，也影响居民的日常生活。当水质受到污染时，应及时对水质进行检测，分析污染类型和来源，以便更好地解决水污染问题，确保居民的水安全。以嘉陵江铊污染问题为例，建立了偏微分模型，利用神经网络和回归分析对其进行了分析。
对于问题1，经计算，判断污水混合长度远小于嘉陵江总长，满足建立一维模型的条件。在此基础上，建立了一维扩散方程模型，将污水混合类型简化为瞬间在断面上完全混合均匀。污染点分为点源扩散型和面源扩散型。将非点源扩散模型视为泥沙通过径流流入河流的污染，建立了不同的偏微分方程模型。重点监测范围为5月5日22:00后37小时，区域位于千佛崖下游30km。
对于第二个问题，基于第一个问题的一维对流扩散微分方程模型，建立高斯扩散模型进行反演。未给出该方程的初始解，但已知满足该方程的特定数据集，并根据该数据反推该方程的初值。同时，由于河流中沉积物中铊的浓度不同，通过分析沉积物中的铊浓度，进一步确定污染位置，即千佛崖上游15公里处为污染源。通过BP神经网络预测和污染类型分析，沉积物中铊的浓度在一段时间后会稳定在较低的水平，铊浓度的最终波动值将保持稳定的幅度，波动值低于标准线，可以判断为点源污染类型。
针对第三个问题，为了分析嘉陵江铊污染数据的特征，利用matlab统计工具，从广元市嘉陵江水污染调查监测报告中统计了部分监测站的检测结果。对嘉陵江各区域水质中铊浓度数据进行方差分析和回归分析，并做出可靠的拟合曲线，分析铊污染特征。结果表明，铊浓度的变化与时间密切相关。拟合了各观测点铊浓度随时间的函数，得到了铊浓度与时间的函数。同时，从各观测点铊浓度的时变图像可以看出，铊浓度随时间先升高后降低，存在延迟。下游的铊浓度在上游开始增加之前会增加一段时间。
对于问题4，分析问题3中的数据以验证问题1和问题2。问题3中的结果满足问题1和问题2中的条件。虽然远离输入点的位置会发生一些变化，但它仍然满足条件。
问题重述：
问题一：结合河道污染情况，通过建立数学模型分析下游水质铊元素浓度变化趋势，从而预估重点污染区域和时间。
问题二：结合问题一下游水质铊元素超标情况，以预估上游水质情况为目的建立模型，求解污染类型以及污染源大致位置。
问题三：分析本次嘉陵江水质铊元素超标的主要特点。
问题四：基于问题三数据，在问题一问题二原有预估不变的前提下，进一步分析和验证问题一以及问题二的合理性。
问题分析：
对于嘉陵江铊元素超标，我们利用matlab统计嘉陵江各个时段的千佛崖、上石盘、张王乡、昭化古镇，沙河镇、沙溪、亭子口等13个监测站的铊元素含量的变化曲线。同时，做出水流速度一定、污水在断面处瞬间混合等理想假设来简化模型，从而建立一维模型。一维模型主要应用于某个河流的较长河段，查阅嘉陵江有关数据经过计算得出嘉陵江河段符合模型条件。
基于第一问的一维对流扩散微分方程模型，我们可以通过建立差分方程反演模型。方程的初始解并未给出，但已知具体的满足方程的一组数据，根据数据反推出方程的初始值。假设污染物扩散速度较慢以高斯扩散方式扩散。同时，因为河流中的泥沙含铊元素浓度不同，用BP神经网络对泥沙中铊浓度进行预测和分析来进一步确定污染类型。
为了分析嘉陵江铊元素污染数据的特点，我们采用了matlab统计工具从广元市嘉陵江水质污染调查监测快报上统计了部分监测站的检测结果数据，并对各地区嘉陵江水质中铊元素的浓度数据做出了回归分析，做出了可靠的拟合曲线对铊元素污染特征进行分析。
对问题三中数据进行分析，从而验证问题一与问题二。问题三中结果满足问题一与问题二的情况，离输入点较远的位置虽然会有一定的变化，但仍然满足条件。
模型的建立与求解整体论文缩略图
部分程序代码：
x1=1:7:630;y1=y_unit(x1);fori=1:87P(1,i)=y1(i);endfori=1:87P(2,i)=y1(i+1);endfori=1:87P(3,i)=y1(i+2);endfori=1:87T(i)=y1(i+3);end%归一化处理
[P,Pmin,Pmax,T,Tmin,Tmax]=premnmx(P,T);%神经网络
net=newff(minmax(P),[5,1],{'tansig','purelin'});%输入隐含层有两层，分别有5个和1个神经元，隐含层传递函数为logistic，输出层为线性
net.trainFcn='trainbfg';%创建网络
%设置训练参数
net.trainParam.show=50;%每隔50步显示一次结果
net.trainParam.lr=0.05;%学习速率
net.trainParam.epochs=10000;%最大学习步数
net.trainParam.goal=1e-5;%训练
[net,tr]=train(net,P,T);%仿真
A=sim(net,P);a=postmnmx(A,Tmin,Tmax);T=postmnmx(T,Tmin,Tmax);%优化后输入层权值和阙值
inputWeights=net.IW{1,1};inputbias=net.b{1};%优化后网络层权值和阙值
layerWeights=net.LW{2,1};layerbias=net.b{2};%画图输出
newk=a(1,:);t=4:90;figure(1)plot(t,newk,'r-o',t,T,'b--*')xlabel('年份')ylabel('浓度')legend('预测浓度','实际浓度')Pnew=[y1(88);y1(89);y1(90)];OldNum=zeros(15,1);fori=1:30SamNum=size(Pnew,2);Pnewn=tramnmx(Pnew,Pmin,Pmax);HiddenOut=tansig(inputWeights*Pnewn+repmat(inputbias,1,SamNum));anewn=purelin(layerWeights*HiddenOut+repmat(layerbias,1,SamNum));anewn=postmnmx(anewn,Tmin,Tmax);Pnew(1:3,:)=[Pnew(2:3,:);anewn];OldNum(i)=anewn;end%画图输出
t1=4:120;figure(2)NUM=[newk,OldNum'];plot(t1,NUM,'r--o')xlabel('时间')ylabel('浓度')legend('预测浓度')