嘉陵江铊污染事件数学建模.pdf

1 嘉陵江铊污染事件 摘要 本文研究的主要是基于嘉陵江铊元素水质污染事件发生后的最优控制策略。依据事 件相关数据的统计分析，建立只考虑铊元素纵向扩散的一维水质模型，以快速估测出水 质污染后的重点监测区域及时间，并能快速查找出相应污染源类型及大致位置。 针对问题一，首先仅考虑铊元素纵向的分散作用，通过研究铊元素的浓度C，江水 的流速u，沿江流动扩散的距离x等参数之间的关系，建立铊元素在稳态条件下的一维 水质模型S x c E x c u t c x 2 。得到浓度变化曲线 x i i eCC 0 ， )9 , 2 , 1 (i，分析得 出下游水质铊元素浓度随污染源的距离的增加呈现出指数型递减。根据建立的模型和变 化规律求解下游的污染重点监控时间段为污染事故发生后h vv C 21. 9ln 0 。重点监控区 域为距离污染源km C 21. 9 ln 0 范围内。 针对问题二，根据下游各个观测点的相对位置距离 i x和监测到铊元素的浓度数据 i C，根据问题一中所建立的水质中铊元素浓度一维水质模型和各个参数之间的规律，将 所收集到的数据带入计算分析。 （1）若经过计算所得到的数据解满足辨别点污染源关系 9 1 )( i i xxdist 则可以判别出该污染类型为点污染源，且污染源位置位于川陕界上游x （2）若不满足满足辨别点污染源关系式，则可判定其类型为面污染源。其边界端点 01 x ，02 x 分别为区域02010 xxX 上下游边界。且污染源区域位于川陕界上游长度为0 X 的范围内。 针对问题三，收集相关数据，利用 Excel 进行相关统计分析，建立基于最小二乘法 思想的数据分析模型，分析铊元素浓度变化相关规律，得出本次铊污染中铊元素浓度在 不同时间段、区域的变化特点。在利用 Excel 相关统计分析基础上，利用 SPSS 软件进 行回归分析，得出各监测站铊元素浓度与时间的显著性概率均小于 0.005，则铊元素浓 度与时间存在着一定的函数关系。通过 MATLAB 进行最小二乘法曲线拟合得出本次各监 测站浓度随时间变化不同， 经过比较分析得出铊元素浓度在不同时间， 区域的变化特点。 针对问题四，首先根据问题三中收集的本次铊元素污染的相关数据，基于问题一中 所建立的稳态条件下的一维水质模型，得出本次嘉陵江铊污染事件应监测的重点时间为 污染发生后 4.25h内，重点监测区域为距离污染源所在地下游 12km范围内。并求解出 本次嘉陵江铊污染类型为点污染源， 污染源的大致位置为嘉陵江距离川陕界上游 7.2km 处。与真实污染源位置燕子砭镇（6.2km）相错 1km，误差概率为 0.16。 关键词：一维水质模型铊元素浓度变化方程最小二乘法拟合绝对值距离 2 一、问题重述 1.1 问题的背景 2017 年 5 月 5 日 18 时，四川广元市环境监测中心站监测发现嘉陵江入川断面水质 异常，西湾水厂水源地水质铊元素超标 4.6 倍。初步判定污染源为川陕界上游输入型、 一次性污染团。此次事件再一次警示民众水污染监控及预警的重要性。 水污染可主要分为点源污染与面源污染两种类型，当出现水质重金属污染超标 时，管理部门需要能尽早快速诊断出污染类型及污染源位置，这样能尽快切断污染源， 减少污染对环境、民众工作生活等带来的损失。 1.2 问题的描述 根据以上信息，请利用数学模型回答以下问题： 1、假设当出现水质铊元素点源污染或面源污染时，建立模型分析下游水质铊元素 浓度变化趋势，由此估计下游的污染重点监控区域及时间。 2、当在下游水质监测点检测到铊元素超标及相关数据，能否利用数学模型快速估 计出上游污染类型及污染源大致位置？ 3、收集本次广元铊元素事件的相关数据，分析本次铊污染的主要特点。 4、根据第三问数据，进一步回答问题 1 和问题 2，并由此验证问题 1 和问题 2 的合 理性。 二、问题分析 此题主要是研究重金属元素（铊元素）在水质中浓度变化问题。首先根据实际问题 建立合理的模型假设， 其次根据建立的描述铊元素水质中浓度变化模型分析设立参数的 变化规律；建立优化模型对重点监测区域，监测时间，污染源的类型及大致位置进行求 解。现在对各个问题的具体分析如下： 2.1 问题一的分析 针对问题一，首先仅考虑纵向的分散（弥散）作用，而不考虑横向和竖向（水深方 向）扩散作用，查阅相关资料，通过研究铊元素的浓度，江水的流速，沿江流动扩散的 距离等参数之间的关系建立铊元素在稳态条件下的一维水质模型。得到浓度变化曲线， 分析下游水质铊元素浓度变化趋势。根据建立的模型和变化规律求解下游的污染重点监 控区域及时间。 2.2 问题二的分析 针对问题二，根据下游各个观测点的相对位置距离和监测到铊元素的浓度数据，根 据问题一中所建立的铊元素浓度变化的一维水质模型和各个参数之间的规律，将问题转 化为利用已知参数求解未知参数的问题。最后将所收集到的数据带入计算分析，快速估 计出上游污染类型及污染源大致位置。 2.3 问题三的分析 针对问题三，首先收集本次嘉陵江铊污染事件中相关数据，根据收集铊污染发生后 广元沿江地段各监测断面某时刻铊浓度相关数据，利用 Excel 进行数据预处理，而后整 理得出某时刻各监测断面所对应铊元素浓度统计表及各监测断面铊元素浓度随时间而 改变的变化趋势图。然后根据铊元素浓度的变化趋势图进行相关分析，可以得出铊元素 浓度的变化与时间相关。最后基于最小二乘法思想，建立数据分析模型，由此分析本次 铊污染的主要特点。 3 2.4 问题四的分析 针对问题四，首先根据问题三中收集的本次铊元素污染的相关数据，基于问题一中 所建立的稳态条件下的一维水质模型，进行相关分析。在已知各监测站之间的距离、水 流速度的条件下，得出本次嘉陵江铊污染发生后所应监测的重点区域及相应时间。而后 利用，精确求解出本次嘉陵江铊污染法污染类型及污染源的大致位置。最后将模型的预 测值与真实值做对比，由此验证出模型的合理性。 三、模型假设 假设 1：河流中河段均匀、恒定连续排污和水文条件稳定 假设 2：河流纵向长度远大于横向竖向长度 假设 3：污染物在河流中纵向扩散作用远大于横向和竖向作用 假设 4：各个监测站所监测到江中铊元素的含量数据真实可靠 假设 5：各个监测站所在地的水流速度数据真实可靠 四、符号说明 符号符号说明 x E纵向分散系数 S附加项（包括衰减项、源和漏） W污染物的输入量 A断面面积 0 C污染源铊元素的初始浓度 i C监测站监测到铊元素的浓度 i u监测站站处测得的的江水流速 i 监测站浓度计算公式所对应常数因子 v位于重点监控范围内监测站处江水流速 i x各个监测站在干流一维转化图中的位置坐标 x中心点解（污染源到川陕界监测站的距离） i x监测站到川陕界监测站的距离 0 X面污染源在干流一维转化图中的区域 01 x面污染源在干流一维转化图中的区域的上游端点 02 x面污染源在干流一维转化图中的区域的下游端点 中心相似度阈值 4 五、问题一模型的建立与求解 5.1.稳态条件下的一维水质模型建立 一维水质模型 4 是仅考虑纵向的分散（弥散）作用，而不考虑横向和竖向（水深方 向）扩散作用的水质模型，推导的基本数学模型如下： S x c E x c u t c x 2 式中： x E纵向分散系数， 12 TL， S附加项（包括衰减项、源和漏）。 如果只考虑衰减项，则KCS ,K为衰减系数： t c 浓度随时间的变化率； x c 浓度随空间（水流方向）的变化率； 当河流中河段均匀、恒定连续排污和水文条件稳定时，断面积A、流速u、污染物 输入量W和E,都不随时间而变化。此时，断面的污染物浓度C是稳定的，0 t c ，则上 式变成 S x c E x c u x 2 2 对于河底无渗漏、忽略侧向输入，对不受抄袭影响的稳态河流 3 ，纵向分数（弥散） 作业很小，易降解的污水污染物水质模型可变为 CK dx dc u 1 此方程的通解为： x AeC 得 u k 对于0x，污染物的降解系数在本文中取 2k，0)(tX， 0 CC 带入边界条件得到 x u k eCC 0 其k的单位为：/1天，u的单位为：sm/，我们化u单位为/km天，则 uuu k 432 1 606024 10002 . 0 5.2 模型的求解 5.2.1 求解浓度与距离的关系式 （1）数据的收集与处理 对于嘉陵江下游，广元市的各个监测点，我们沿江流方向，从上到下依次将各个监 测点的位置在嘉陵江干流示意图中标注 5 图 1：嘉陵江广元段各监测站示意图 我们对嘉陵江干流的九个主要监测点从 2016 年 4 月到 2017 年 4 月的流速进行了统计， 得到了各地的平均流速。 表格 1：嘉陵江广元段各监测站之间距离 川陕界清风峡沙河镇千佛崖苴国村上石盘昭化古镇张王乡八庙沟 3.098762.948572.690852.564741.959842.213332.094771.984051.96923 （2）浓度变化方程的建立 将上表九个监测点的流速分别带入式得到下游各个监测点浓度的变化方程求解的 结果 x x x x x x x x x eCC eCC eCC eCC eCC eCC eCC eCC eCC 00118. 0 09 00117. 0 08 00111. 0 07 00105. 0 06 00118. 0 05 00090. 0 04 00086. 0 03 00079. 0 02 00075. 0 01 5.2.2 模型求解的结果 题目要求所求下游的污染重点监控区域，即为铊浓度超标区域。将污染源处排放的 废渣中铊浓度记为LmgC/ 0 ，超标区域初始铊浓度大于Lmg /0001. 0，即 LmgeCC x /0001. 0 0 6 通过解出上式得出 21. 9 ln 0 C x 因此，发生铊元素水污染后，距离污染源的下游区域km C 21. 9 ln 0 范围内均为污染重 点监控区域。记位于污染重点监控区域内的检测站流域水速为v，则污染事故发生后重 点监测时间段为h vv C 21. 9ln 0 5.3 利用 MATLAB 进行模型检验与评价 依据表 1 中的相关数据，利用 MATLAB 编程将基于稳态条件下的一维水质模型进 行检验，则各监测断面处铊元素浓度随距离变化曲线图，具体结果如下： 图 2：监测断面处铊元素浓度随距离变化曲线图 六、问题二的模型建立与求解 6.1 一维坐标转化模型的建立 由问题一中的一维水质模型我们能够得到嘉陵江下游广元市的各个监测站的污染 物浓度变化模型。我们将嘉陵江干流部分示意图沿着水流方向，建立一维坐标轴，将沿 岸各个监测点记做坐标轴上的点， 记沿江流方向各个监测点川陕界， 清风峡， 沙河镇 八 里庙分别为 1 x， 2 x， 3 x 9 x由下图轴（1）可知污染源应处于嘉陵江干流上游位置， 即下图 中 1 x的左端区域。 （1） 如果污染源的类型为点污染源， 可记上游点污染源在一维坐标转化图的位置为 0 x， 0 x位于 1 x的左端。即如下图轴（2）所示。 （2） 如果污染源的类型为面污染源， 则将上游面污染源在一维坐标检测站中的位置 区域为 0 X，其左右端点界限分别为 01 x， 02 x。即如下图轴（3）所示。 7 图 3：嘉陵江广元段监测站一维坐标示意图 6.2 数据的统计和处理 收集到的下游各个观测点的相对位置距离和监测到铊元素的浓度数据进行分析 （详 细数据见附录 2） 表格 2：嘉陵江沿程各段距离 上游点污染源在一维坐标转化图的位置为0 x 到第一个监测点川陕界的距离记做 0 x ，可以得到各个监测点的坐标值01 xx ， 112 xxx ， 223 xxx ， ， 889 xxx 。 统计这九个监测站在m天之内每天所测得铊元素的均值，每个监测站得到的监测数 据即为一个m维向量， imijiii ccccC, 21 ，9 , 2 , 1i 根据统计的数据结果得到 9 组m维向量，记做下表 检测地点 ji xx 符号 i x距离km 川陕界至清风峡 21 xx 20 清风峡至沙河镇 32 xx 11.4 沙河镇至千佛崖 43 xx 10.5 千佛崖至苴国村 54 xx 5 苴国村至上石盘 65 xx 9 上石盘至昭化古镇 76 xx 3.8 昭化古镇至张王乡 87 xx 23.7 张王乡至八庙沟 98 xx 55.6 8 表格 3：各监测站对应的 9 组m维向量 川陕界清风峡沙河镇千佛崖苴国村上石盘昭化古镇张王乡八庙沟 1 C 2 C 3 C 4 C 5 C 6 C 7 C 8 C 9 C 6.3 模型的求解 根据所收集到的下游各个观测点的相对位置距离和监测到铊元素的浓度数据 （详细 数据见附录）进行分析，在下游水质监测点检测到铊元素超标及相关数据带入问题一中 的铊元素浓度关系式中，可以得到 9m组的数据解 )1( x， )2( x， )3( x， )4( x， )5( x， )6( x， )7( x， )8( x， )9( x 记九组数据的中心点解为 9 )9()2()1( xxx x 并设立中心相似度阈值2 . 0对所得的九组数据进行分析： （1）若经过计算所得到的数据解 )1( x， )2( x， )3( x， )4( x， )5( x， )6( x， )7( x， )8( x， )9( x 满足辨别点污染源关系式 9 1 )( i i xxdist 则可以判别出该污染类型为点污染源，且污染源位置位于川陕界上游x （2）若经过计算所得到的 )1( x， )2( x， )3( x， )4( x， )5( x， )6( x， )7( x， )8( x， )9( x 不满足满足辨别点污染源关系式， 则可判定其类型为面污染源。 其边界端点01 x ，02 x 分别为 )9()2()1( 101 ,maxxxxxx )9()2()1( 102 ,minxxxxx 02010 xxX 则可以判别出该污染类型为面污染源， 且污染源区域位于川陕界上游 01 x到 02 x长度为 0 X 的范围内。 七、问题三模型的建立与求解 7.1 模型的准备 7.1.1 理论基础 最小二乘法进行拟合逼近原始数据。 常见的曲线拟合方法: 1.使偏差绝对值之和最小 m i ii m i i yx 11 min 2.使偏差绝对值最大的最小 9 iii yx maxmin 3.使偏差平方和最小 2 11 2 min m i ii m i i yx 按偏差平方和最小的原则选取拟合曲线，并且采取二项式方程为拟合曲线 5 的方法,称 为最小二乘法。 推导过程： 1.设拟合多项式为： k kx axaay 10 2.各点到这条曲线的距离之和，即偏差平方和如下： 2 1 110 2 n i k iki xaxaayR 3.为了求得符合条件的a值，对等式右边求 i a偏导数，因而我们得到了： 02- 1 10 n i k kx axaay 02- 1 10 xxaxaay n i k k . 02- 1 10 k n i k k xxaxaay 4.将等式左边进行一下化简，将得到的等式表示成矩阵 6 的形式，就可以得到系数矩 阵 A： n i i k i n i ii n i i k n i k i n i k i n i k i n i k i n i i n i i n i k i n i i yx yx y a a a x xxx xxx xxn 1 1 1 1 0 1 2 1 1 1 1 1 1 2 1 11 5.将这个范德蒙得矩阵化简后可得到: 10 nk k nn k k y y y a a a xx xx xx 2 1 1 0 22 11 1 1 1 7.2 数据的预处理 7.2.1 背景知识 于 5 月 5 日 18 时， 广元市环境监测中心站监测发现嘉陵江入川断面水质异常， 西湾 水厂水源地水质铊元素超标 4.6 倍(国家地表水环境质量标准为 0.0001 /L)。而后于 广元段分别对川陕界、八庙、清风峡、沙河镇、千佛崖、上石盘、苴国村、昭化古镇、 张王乡、亭子口、沙溪、工农镇供水站、工农贺连利供水站等 13 个监测站对嘉陵江内 铊元素的浓度进行观测。 则下表为水质铊污染发生后 5月6日的部分监测结果及评价表： 表 4监测结果及评价表 监测断面采样时间铊 （mg/L）评价结果 上次监测 铊浓度 (mg/L) 川陕界5 月 6 日 8 时0.00004达标 0.00006 (5月6日2 时) 八庙沟（嘉陵江）5 月 6 日 8 时0.00006达标0.00008 清风峡（嘉陵江）5 月 6 日 8 时0.00008达标0.00012 沙河镇（嘉陵江）5 月 6 日 8 时0.00016超标超标 0.60.6 倍倍0.00018 千佛崖（嘉陵江）5 月 6 日 8 时0.00026超标超标 1.61.6 倍倍 0.00030 (5月6日7 时) 上石盘（嘉陵江）5 月 6 日 8 时0.00019超标超标 0.90.9 倍倍0.00016 苴国村（白龙江， 清洁对照断面） 5 月 6 日 7 时未检出达标 昭化古镇（嘉陵 江） 5 月 6 日 7 时0.00006达标0.00003 张王乡（嘉陵江）5 月 6 日 8 时未检出 亭子口（嘉陵江）5 月 6 日 8 时未检出达标未检出 沙溪（嘉陵江）5 月 6 日 8 时未检出达标未检出 工农镇供水站5 月 6 日 8 时未检出达标未检出 工农贺连利供水 站 5 月 6 日 8 时未检出达标未检出 11 收集这 13 个监测点每时段的监测结果，利用 Excel 进行统计分析（因为版面限制只展 示部分数据，详细请见附录 2）得： 表格 5：5 月 6 日各断面铊元素监测结果 2 时8 时10 时12 时14 时16 时18 时20 时22 时24 时 川 陕 界 0.00 006 0.00 004 0.00 007 0.00 003 0.00 003 0.00 003 八 庙 沟 0.00 008 0.00 006 0.00 005 0.00 005 0.00 004 清 风 峡 0.00 012 0.00 008 0.00 006 0.00 006 0.00 005 0.00 005 0.00 005 沙 河 镇 0.00 018 0.00 016 0.00 013 0.00 012 0.00 009 0.00 008 0.00 007 0.00 006 千 佛 崖 0.00 030 0.00 026 0.00 021 0.00 015 0.00 014 0.00 013 0.00 011 0.00 011 0.00 010 0.00 009 上 石 盘 0.00 016 0.00 019 0.00 018 0.00 021 0.00 024 0.00 025 0.00 026 0.00 026 0.00 026 0.00 028 昭 化 古镇 0.00 003 0.00 006 0.00 009 0.00 007 0.00 008 0.00 008 0.00 008 0.00 007 0.00 007 0.00 010 根据事故发生后各断面的铊元素浓度监测结果发现，各监测断面铊元素的浓度随着 时间的变化不断的发生变化，则利用 Excel 进行统计分析做出 5 月 6 日 2 时至 5 月 11 日 12 时各监测断面铊元素浓度变化趋势图： 图 4：各监测断面铊元素浓度变化趋势图 由图 4 可以看出浓度c与时间t存在着一定的关系。 7.3 基于数据统计分析模型的求解 12 7.3.1 利用 SPSS 软件进行回归分析 利用 SPSS 软件对本次嘉陵江铊污染事件的相关数据进行分析。从中选取其中一个 监测断面（上石盘）5 月 6 日 2 时至 5 月 11 日 12 时的铊元素浓度c与时间t相关关系， 具体数据参见附录 3。 图 5：上石盘监测断面铊元素浓度散点图 结论：由散点图可知，浓度c与时间t存在着一定的函数关系，则对其进行线性回 归，进一步检测浓度c与时间t的函数关系。 进行曲线估计后即如下图所示： 图 6：上石盘监测断面铊元素浓度随时间变化曲线估计图 13 7.3.2 利用 MATLAB 进行最小二乘法拟合 依据 MATLAB 编程（详细程序附录三）求出各监测站铊元素浓度与时间相关系数， 得出各监测站铊元素浓度变化趋势，发现沙河镇，上石盘、千佛崖、张王乡等地铊元素 变化趋势明显； 依据各监测站铊元素随时间变化梯度判断出本次铊元素污染的主要特点 为污染源出现在嘉陵江广元段上游。 将统计所得的数据通过 Matlab 程序进行最小二乘法拟合， （详细程序见附录四）得 到各监测断面铊元素浓度拟合曲线，如下图所示： 图 7：各监测断面铊元素浓度拟合曲线 八、问题四的模型建立与求解 8.1 模型的求解 8.1.1 对于问题一的解答 根据问题一中所建立的模型即下游各个监测点浓度的变化方程利用问题三中所收集 到广元市各个监测站初始时刻检测到铊元素浓度9, 2 , 1iCi x x x x x x x x x eCC eCC eCC eCC eCC eCC eCC eCC eCC 00118. 0 09 00117. 0 08 00111. 0 07 00105. 0 06 00118. 0 05 00090. 0 04 00086. 0 03 00079. 0 02 00075. 0 01 所求下游的污染重点监控区域， 即为铊浓度超标区域。 由查找到相关数据可以得知， LmgC/52. 1 0 ，超标区域初始铊浓度大于Lmg /0001. 0，即 LmgeC x /0001. 052. 1 通过解出上式得出 14 9 , 2 , 121. 9 ln max 0 i C Sx i 因此经计算上式得到kmS838.12。发生铊元素水污染后，距离污染源的下游区域km13 范围内均为污染重点监控区域。记位于污染重点监控区域内的检测站流域水速为 hkmv/0568. 3，则污染事故发生后重点监测时间段为 4.25h. 8.1.2 对于问题二的解答 上游点污染源在一维坐标转化图的位置为0 x 到第一个监测点川陕界的距离记做 0 x ，可以得到各个监测点的坐标值01 xx ， 20 2 xx ， 4 .31 3 xx ， 139 9 xx 。 统计这九个监测站在m天之内每天所测得铊元素的均值，每个监测站得到的监测数 据即为一个m维向量， imijiii ccccC, 21 ，9 , 2 , 1i 根据统计的数据结果得到 9 组m维向量。 x x x x x x x x x eCC eCC eCC eCC eCC eCC eCC eCC eCC 00118. 0 09 00117. 0 08 00111. 0 07 00105. 0 06 00118. 0 05 00090. 0 04 00086. 0 03 00079. 0 02 00075. 0 01 在下游水质监测点检测到铊元素超标及相关数据带入上式，可以得到九组的数据解 )1( x， )2( x， )3( x， )4( x， )5( x， )6( x， )7( x， )8( x， )9( x 经过计算可以得到九组数据的中心点解为 2 . 7 9 )9()2()1( xxx x 并设立中心相似度阈值2 对所得的九组数据进行分析，若经过计算所得到的 )1( x， )2( x， )3( x， )4( x， )5( x， )6( x， )7( x， )8( x， )9( x 满足辨别点污染源关系式 213. 1 9 1 )( i i xxdist 则可以判别出该污染类型为点污染源，且污染源位置位于川陕界上游2 . 7x 8.2 模型的检验与评价 根据收集本次广元铊元素事件的相关数据，我们利用问题一中的一维水质模型，对 嘉陵江铊元素污染事件进行研究。 铊元素污染发生后 4 小时 20 分钟内为重点监测时间， 距离污染源所在地下游 12km范围内为污染重点监测区域。利用问题二中的模型，带入 15 各个监测站中监测到铊元素的浓度和各个监测站沿江所在的相对位置，我们可以判断出 此次嘉陵江铊元素污染类型为点源污染，且污染源应位于嘉陵江川陕交界处（广元市与 宁强县交界）上游沿江 7.2km处。 通过此次嘉陵江铊元素污染事件的调查显示，此次污染事件 1 的污染源位于嘉陵江 上游宁强县燕子砭镇， 污染类型为点污染源。 且污染源位于嘉陵江川陕交界处 6.2km处。 将所求得的预测值与实际值相比较，污染源位置相差 1km，相对误差为 0.16。从而验证 出模型的建立合理性，与实际值基本一致。 九、模型评价与改进 9.1 模型的优缺点 9.1.1 模型的优点 本文采用了一维稳态水质模型对嘉陵江铊元素污染事件重点监控区域、时间进行预 测和分析，具有严谨的推导过程，且模型适用性强。 在问题三中利用最小二乘法对铊元素浓度与时间相关数据进行曲线拟合，得出铊元 素浓度在不同时间、不同区域的变化特点。 9.1.2 模型的缺点 问题一中所建立的一维稳态水质模型在均匀河段上定常排污条件下，河段横截面、 流速。流量、污染物输入量和弥散系数都不随时间变化。仅仅考虑了污染物纵向扩散作 用，但并未考虑污染物在水中的横向扩散作用及垂向扩散作用。 9.2 模型的改进与推广 9.2.1 模型的改进 1、问题一模型的改进 问题一所建立的一维稳态水质模型仅仅考虑污染物的纵向扩散作用，适用于预测位 于均匀混合端面下游的敏感点。在此基础上可以对一维稳态水质模型加以改进，建立二 维水质模型，不仅考虑污染物的纵向扩散作用，还考虑横向扩散的作用，用于预测敏感 点位于均匀混合端面的上游，以及污染带分布及超标水域的界定。 2、问题二模型的改进 问题二中依据问题一所建立的一维稳态水质模型，求解出铊元素浓度变化方程，预 测出污染源类型和大致位置，对于数据处理中所出现的误差较大。此题还可以还可以利 用遗传算法进行改进。 在一维稳态水质模型建立的基础上利用遗传算法，构建水污染源遗传算法辨识模 型，追踪上游污染源大致位置并判别相应污染类型。利用数学分析法在获得地表水控制 方程基础上采用优化法反演出水污染事故的污染源。其中遗传算法具备全面搜索的能 力，非常适合应用到水污染源反演模型中去。 9.2.1 模型的推广 本文主要应用一维稳态水质模型及转化方程，快速预测出重点监控区域、污染类型 及污染源大致位置，对本次嘉陵江铊元素污染事件发生后做出最优控制策略。一维稳态 水质模型不仅适用于污染物监测领域，还可广泛应用与其他领域，突出表现为基于一维 水质模型估算河流动态纳污能力，预测不停污水排放量对指定河段水质的不同影响。同 时还可以根据一维水质模型，进行对多参数识别的反演优化算法，预测水库水质变化趋 势等应用。 16 十、参考文献 1数学建模长江水质污染分析与预测， qq-pf-to=pcqq.group，2017.7 2广元市环境保护局， 3广元市嘉陵江水质污染调查监测报, ml?id=13410，2017.7 4韩中庚著.数学建模方法及其应用.北京.高等教育出版社,2009.6 5李明著.详解 MATLAB 在最优化计算中的应用.北京.电子工业出版社,2011.5 6司守奎孙玺菁著.数学建模算法与应用.北京.国防工业出版社.2014.9 17 十一、附录 附录 1 关于铊浓度随距离的变化函数曲线的 matlab 实现 问题一通过对 matlab 中 plot()函数的调用，绘制出铊浓度随距离的变化曲线，并且 利用 line()，ginput()函数绘制在LmgC/00001. 0条件下 0 /CCy 曲线，同时求得交 点坐标，matlab 实现程序如下： x=linspace(0,20000); y1=exp(-0.00075*x); y2=exp(-0.00079*x); y3=exp(-0.00086*x); y4=exp(-0.00090*x); y5=exp(-0.00118*x); y6=exp(-0.00105*x); y7=exp(-0.00111*x); y8=exp(-0.00117*x); y9=exp(-0.00118*x); plot(x,y1); hold on plot(x,y2); plot(x,y3) plot(x,y4) plot(x,y5) plot(x,y6) plot(x,y7) plot(x,y8) plot(x,y9) hold on line(0,10000,6.5785e-5,6.5785e-5) x,y = ginput(1) 附录 2 监测数据统计 %因页面限制，部分统计数据展示如下： 各检测断面铊浓度检测结果汇总（部分） 监测时间监测点 川陕界 八庙沟 清风 峡 沙河镇 千佛崖 上石 盘 昭化古 镇 张王 乡 江口镇 5 月 6 日 2 时68121830163 5 月 6 日 8 时4681626196 5 月 6 日 10 时561321189 5 月 6 日 12 时7561215217 5 月 6 日 14 时914248 18 5 月 6 日 16 时345813258 5 月 6 日 18 时11268 5 月 6 日 20 时35711267 5 月 6 日 22 时610267 5 月 6 日 24 时3592810 5 月 7 日 2 时82812 5 月 7 日 4 时3472912 5 月 7 日 6 时63314 5 月 7 日 7 时716 5 月 7 日 8 时34633 5 月 7 日 9 时618 5 月 7 日 10 时636 5 月 7 日 11 时20 5 月 7 日 12 时346343 5 月 7 日 13 时21 5 月 7 日 14 时34 5 月 7 日 15 时18 5 月 7 日 16 时530 5 月 7 日 17 时18 5 月 7 日 18 时305 5 月 7 日 19 时 5 月 7 日 20 时45635 5 月 7 日 21 时19 5 月 7 日 22 时 5 月 7 日 23 时13 5 月 7 日 24 时3457 5 月 8 日 1 时 5 月 8 日 2 时35 5 月 8 日 3 时 5 月 8 日 4 时3433 5 月 8 日 5 时19 5 月 8 日 6 时346 5 月 8 日 7 时22 5 月 8 日 8 时24431 5 月 8 日 9 时24 5 月 8 日 10 时30 5 月 8 日 11 时24 5 月 8 日 12 时234328 5 月 8 日 13 时12 19 5 月 8 日 14 时29 5 月 8 日 15 时6 5 月 8 日 16 时22424 5 月 8 日 17 时 5 月 8 日 18 时227 5 月 8 日 19 时7 5 月 8 日 20 时02520 5 月 8 日 21 时5 5 月 8 日 22 时23 5 月 8 日 23 时4 5 月 8 日 24 时23188 5 月 9 日 1 时6 5 月 9 日 2 时17 5 月 9 日 3 时10 5 月 9 日 4 时3419 5 月 9 日 5 时 5 月 9 日 6 时2113 5 月 9 日 7 时13 5 月 9 日 8 时3418 5 月 9 日 9 时 5 月 9 日 10 时 5 月 9 日 11 时13 5 月 9 日 12 时1714 5 月 9 日 13 时7 5 月 9 日 14 时18 5 月 9 日 15 时3 5 月 9 日 16 时331714 5 月 9 日 17 时 5 月 9 日 18 时16 5 月 9 日 19 时3 5 月 9 日 20 时15 5 月 9 日 21 时 5 月 9 日 22 时143 5 月 9 日 23 时 5 月 9 日 24 时23153 5 月 10 日 1 时2316 5 月 10 日 2 时144 5 月 10 日 3 时 5 月 10 日 4 时13314 20 附录 3 利用最小二乘法对采集数据二次多项式拟合的 matlab 实现 问题三中利用拟合多项式的系数对数据进行分析，利用 polyfit()函数拟合多项式并 得到各项系数，matlab 实现如下： y1=64733333343222032200 2; x1=0610141822263034424654586266869495 102 118 126; 5 月 10 日 5 时 5 月 10 日 6 时123 5 月 10 日 7 时 5 月 10 日 8 时027132 5 月 10 日 9 时 5 月 10 日 10 时12 5 月 10 日 11 时 5 月 10 日 12 时134 5 月 10 日 13 时 5 月 10 日 14 时1034 5 月 10 日 15 时 5 月 10 日 16 时290124 5 月 10 日 17 时 5 月 10 日 18 时90 5 月 10 日 19 时 5 月 10 日 20 时100 5 月 10 日 21 时 5 月 10 日 22 时80 5 月 10 日 23 时 5 月 10 日 24 时008095 5 月 11 日 1 时 5 月 11 日 2 时 5 月 11 日 3 时 5 月 11 日 4 时86 5 月 11 日 5 时 5 月 11 日 6 时 5 月 11 日 7 时 5 月 11 日 8 时24861010 5 月 11 日 9 时 5 月 11 日 10 时 5 月 11 日 11 时 5 月 11 日 12 时1010 21 p1=polyfit(x1,y1,2)%川陕界 y2=86554; x2=0681014; p2=polyfit(x2,y2,2)%八庙沟 y3=128 665554445434322233; x3=06810141822243034424650545862667074 78; p3=polyfit(x3,y3,2)%清风峡 y4=18 16 13 12 9876; x4=0681012141820; p4=polyfit(x4,y4,2)%沙河镇 y5=30 26 21 15 14 13 11 11 10 987676666565 444453443332204; x5=068101214161820222426282930313234 3842465054586266707478869495102 110 118 126; p5=polyfit(x5,y5,2)%千佛崖 y6=16 19 18