简介长江水质评价与预测数学建模文档内容完整适合数学建模竞赛参赛者、环境科学与水利工程专业学生及从事水质分析的研究人员使用。文档围绕四个核心问题展开基于模糊综合评价法对长江近两年水质进行定级构建主要污染源判别模型确定高锰酸盐与氨氮的重点排放区域利用灰色系统模型预测未来十年各类水质河长比例并计算为控制Ⅳ、Ⅴ类水比例所需处理的污水总量。全文涵盖模型假设、符号说明、问题重述、求解步骤及结果表格并附有对环保部门的参考性建议既可作为数学建模论文写作范例也可用于学习模糊综合评价、灰色预测等方法的环境应用。资源为1个doc文档压缩包大小595KB已有349人学习。文档结构清晰从问题重述到模型求解与结果展示一应俱全适合需要快速理解水质评价与预测建模思路的读者下载研读。1. 长江水质的评价和预测先把“等级”和“趋势”拆成两件事拿到“长江水质的评价和预测”这份文档的人十有八九第一反应是打开附件看数据然后直接上模型。我第一次做的时候也是这个顺序结果评价出来的等级和常识对不上预测曲线平滑得像画上去的。后来想明白一件事长江水质的评价和预测是两条完全不同的技术链路。评价是一次快照处理同一时间断面上多个指标的合成问题核心矛盾在于权重怎么定、最差指标要不要一票否决预测是时间序列处理同一断面多个指标在几十个月里的演变核心矛盾在于样本量太少、季节性太强、缺测太多。把这两条链路混在一起想就会陷入“用评价的等级序列直接做预测”的误区——等级是离散标签跳变多、信息量低拿它喂给时间序列模型残差一定难看。这套东西适合环境类、水利类的建模竞赛选手、以及要做流域水质报表的工程技术人员只要能拿到断面监测原始表就能按下面的路径一步步复现。2. 把附件里的监测表变成可计算的矩阵断面、月份、指标三轴对齐打开附件通常是几张排版不规整的表格列是断面名称行是监测时间中间夹着溶解氧、高锰酸盐指数、氨氮这些指标有的表头还带单位后缀。真正开始算之前有三件事必须先做完把评价标准落到代码里、把宽表摊成长表、把缺测和检出限处理掉。这三步做歪了后面无论用多花哨的模型结论都是错的。我在第一次做的时候跳过了第二步直接在宽表上循环遍历列名代码又长又容易在下一次换数据格式时崩掉血泪经验。2.1 先把 GB 3838 的五类限值写成一张能查的表评价的基础是《地表水环境质量标准》里的限值分级。不同指标的分界值不一样而且方向不一致——大部分指标越小越好溶解氧恰恰相反越大越好。把限值硬编码进if-else是灾难的开始正确做法是做成结构化字典后续查表统一调用。指标I类II类III类IV类V类优劣方向溶解氧 (mg/L)≥7.5≥6≥5≥3≥2越大越好高锰酸盐指数 (mg/L)≤2≤4≤6≤10≤15越小越好氨氮 (mg/L)≤0.15≤0.5≤1.0≤1.5≤2.0越小越好五日生化需氧量 (mg/L)≤3≤3≤4≤6≤10越小越好总磷河流mg/L≤0.02≤0.1≤0.2≤0.3≤0.4越小越好pH6~9无分级超出即不达标区间型这张表我最常用来做两件事一是把每个断面的实测值映射成 1~5 的类别号二是给综合指数提供分母基准。用水质类别做分母还是用 III 类标准值做分母算出来的指数完全不可比论文里必须写清楚用的是哪一种。2.2 宽表摊平成长表三步 reshape 把结构定死附件里的表天然是宽表行是时间、列是断面或指标这种结构做分组运算很别扭。我一般先把它melt成长表每个指标压成一列这样按“断面 指标”分组就顺了。import pandas as pd import numpy as np raw pd.read_excel(changjiang.xlsx, sheet_name监测数据) # 1) 列名规整去掉空格和单位后缀避免同一指标两个拼法 raw.columns [c.strip().replace(mg/L, ).replace((mg/L), ) for c in raw.columns] # 2) 宽表转长表id_vars 是不参与变形的维度列其余列全部转成“指标/测值”两列 long raw.melt( id_vars[断面, 监测时间], var_name指标, value_name测值 ) # 3) 时间列统一成 datetime 再锚定到月防止 2003-6 和 2003-06 被当成两个月 long[监测时间] pd.to_datetime(long[监测时间], errorscoerce) long[月份] long[监测时间].dt.to_period(M) # 4) 数值清洗0.01、未检出、-- 统一处理 def to_num(v): if pd.isna(v): return np.nan s str(v).strip() if s.startswith(): # 低于检出限取 1/2 检出限是通行做法 try: return float(s[1:]) / 2 except ValueError: return np.nan try: return float(v) except ValueError: return np.nan long[测值] long[测值].map(to_num) # 5) 先做一次体检每个指标的样本量、最小值和最大值负值一眼就能看出来 print(long.groupby(指标)[测值].agg([count, min, max]))id_vars指定的是保持不动的列选错了会把断面名也当成指标搅进去。errorscoerce让无法解析的时间变成NaT而不是直接抛异常方便一次性看到所有格式问题。半检出限取值是最常见的处理方式但如果某一指标的检出限占标准限值的比例很大这个近似会系统性抬高整体浓度必须在限制条件里说明。最后那行print是黑匣子我几乎每个项目都会顺手跑一遍最大值为负、最小值比检出限还低一个量级这些异常在这一步就能拦住。2.3 缺测怎么补按断面插值别用全局均值缺测在断面月序列里非常普遍尤其是枯水期某些断面停测。最常见的错误做法是用全表均值填补这会把时间趋势抹平后面预测出来的曲线自然就是一条水平直线。long long.sort_values([断面, 指标, 监测时间]) # 每个断面、每个指标各自插值时间轴上做线性内插 long[测值_填充] ( long.groupby([断面, 指标])[测值] .transform(lambda s: s.interpolate(methodlinear, limit2, limit_directionboth)) ) # 覆盖率体检整段缺失的指标不要插直接标记为不可评价 coverage long.groupby([断面, 指标])[测值].apply(lambda s: s.notna().mean()) print(coverage[coverage 0.7])limit2表示最多连续补两期超过两期的空洞说明监测体系本身有问题硬补出来的数字支撑不了一个评价结论。limit_directionboth处理序列首尾否则开头结尾的缺口永远是空的。覆盖率阈值 0.7 是经验值低于这个比例我一般直接把这个“断面 × 指标”组合排除出评价矩阵而不是插出一个看起来完整的数。异常值处理上用一个双约束3σ 之外、且超出物理下界浓度不能为负的点先标记再人工判断别直接删掉删掉的可能是真实的污染事件。3. 水质等级怎么定单因子、内梅罗与灰色聚类三条路线的取舍评价环节是整个工作的地基。18 个断面、6 项指标、28 个月最简单的做法是逐月逐断面算出一个类别号但“类别号”本身是粗粒度的做趋势分析时会发现大量并列和跳变。我一般同时跑三条路线单因子看最差项内梅罗看综合水平灰色聚类看整体归类三者结论一致时结论才站得住不一致时反而是论文里最有价值的讨论点。3.1 单因子评价一票否决把每个指标映射成类别单因子法的逻辑很硬取所有指标中最差的那一项作为该断面的类别。它对应的是“水质功能不能因为某一项达标就判定合格”的管理思路代码实现上就是一次查表。import numpy as np import pandas as pd # 越小越优型指标的类界升序排列第一段即 I 类 BOUNDS { 高锰酸盐指数: [2, 4, 6, 10, 15], 氨氮: [0.15, 0.5, 1.0, 1.5, 2.0], 五日生化需氧量: [3, 3, 4, 6, 10], 总磷: [0.02, 0.1, 0.2, 0.3, 0.4], } def class_lower_better(value, bounds): 越界返回 6表示劣于 V 类。 if pd.isna(value): return np.nan # sideleft 对应“ 界值即归入该类”的语义 return int(np.searchsorted(bounds, value, sideleft)) 1 def class_dissolved_oxygen(value): 溶解氧越大越好反向查表。 if pd.isna(value): return np.nan for i, low in enumerate([7.5, 6, 5, 3, 2], start1): if value low: return i return 6np.searchsorted的side参数是这里唯一的玄学点sideleft在值恰好等于界值时返回前一段符合标准里“小于等于”的表述换成right会出现刚好 4.0 的氨氮被判成 III 类的情况。溶解氧必须单独走一条函数因为它是唯一越大越优的指标混进统一循环里一定会反过来。断面类别取所有指标类别的最大值极端情况下一个断面的氨氮劣 V 类其余全是 I 类单因子结论就是劣 V——这不是 bug是这套方法的固有属性写报告时要说明。3.2 内梅罗指数与加权综合为什么权重不能拍脑袋单因子给的是极端信息内梅罗指数给的是整体信息。公式是P sqrt((P_avg² P_max²) / 2)把平均值和最大值放在同等权重下平方合成既照顾整体水平也对最差项给出惩罚。# III 类标准值作为基准所有指标统一到这个尺度上才有可比性 STD3 { 溶解氧: 5.0, 高锰酸盐指数: 6.0, 氨氮: 1.0, 五日生化需氧量: 4.0, 总磷: 0.2, } def nemerow(row, std3STD3): row 是某断面某月的 {指标: 测值}。 ps [] for k, v in row.items(): if pd.isna(v) or k not in std3: continue if k 溶解氧: ps.append(std3[k] / v) # 反向指标取比值倒数 else: ps.append(v / std3[k]) if not ps: return np.nan ps np.array(ps) return float(np.sqrt((ps.mean() ** 2 ps.max() ** 2) / 2))平方合成意味着最大值那一项实际拿到了约 0.5 的权重pH 这类没有分级的指标不要塞进来。最常见的坑是不同断面参与计算的指标个数不一致——某个断面缺了总磷内梅罗值天然会偏小看起来“更干净”。解决办法是固定指标集合缺任何一项就整条记录打上缺失标记宁缺毋滥。3.3 灰色聚类与熵权法让权重从数据里长出来到这一步会遇到一个绕不开的问题六项指标谁更重要。拍脑袋给权重在评审面前站不住纯客观赋权又容易被异常值带偏。我一般的做法是熵权法定权重、灰色聚类做归类两者互为校验。def entropy_weight(mat): mat: 样本 × 指标要求已同向化。 x mat.astype(float) x (x - x.min(0)) / (x.max(0) - x.min(0) 1e-12) # 极差归一化 x x 1e-6 # 避免 log(0) p x / x.sum(0, keepdimsTrue) e -(p * np.log(p)).sum(0) / np.log(len(x)) # 各指标信息熵 w (1 - e) / (1 - e).sum() # 差异越大权重越高 return w同向化必须在归一化之前完成溶解氧取倒数或者反向归一化否则熵值算出来是反的。1e-6是为了避开零值取对数取值大小影响有限但不能省。熵权法对样本量敏感只有十来个断面时权重抖得厉害这时我会把它和 AHP 主观权重各取一半加权平均既保留数据信息又不会因为某一个断面的极端值把权重全吸走。灰色聚类部分用白化函数构造各灰类的隶属度按最大隶属度定级实现方式与模糊综合评价非常接近选哪种主要看报告里想强调“信息不完全”还是“边界模糊”。3.4 从“这段水质差”到“上游排了多少”一维水质模型反演评价只能回答哪里差回答不了差从哪来。要做污染来源分析就得引入一维稳态水质模型污染物从上游断面往下游迁移的过程中按指数衰减浓度变化遵循C(x) C0 · exp(-k·x/u)。def emission_between(Q, C_up, C_down, L_km, u_km_per_day, k_per_day): 两断面之间新增的日排放量返回 kg/d。 Q: 断面平均流量 m3/sC: 浓度 mg/Lk: 综合衰减系数 1/d decay np.exp(-k_per_day * L_km / u_km_per_day) delta C_down - C_up * decay # 扣除自然衰减后真正新增的浓度 return delta * Q * 86400 / 1000k是最需要标定的参数氨氮常见在 0.05~0.25 /d、高锰酸盐指数在 0.02~0.1 /d 区间具体值要用区间内多组上下游浓度反算并取稳健值写死在代码里必翻车。86400/1000是把 m³/s 与 mg/L 换算成 kg/d 的系数量纲错了结果差一千倍。如果delta算出来是负的说明衰减系数取大了、或者区间内有支流稀释不要强行解释成“负排放”那是模型在提醒你边界条件没设对。4. 未来水质怎么走GM(1,1)、ARIMA 和 LSTM 的适用边界预测环节最容易陷入的思维定式是“越复杂的模型越好”。断面月序列通常只有几十个时间步深度学习模型的参数量动辄上万这种数据体量下复杂模型往往输给灰色预测这种看起来朴素的方法。我一般的策略是用 GM(1,1) 打底用 ARIMA 处理季节性用 LSTM 做上界参照最后用滚动回测统一裁决。4.1 GM(1,1)几十个点的小样本先用它打底灰色预测对样本量要求低一般 4 个点以上就能建模代价是对数据光滑度有要求。建模前必须做级比检验落不到可容区间就说明数据不适合直接建模。import numpy as np def gm11(x0): x0 np.asarray(x0, dtypefloat) n len(x0) # 1) 级比检验lambda(k) x(k-1)/x(k) 需落在 (e^{-2/(n1)}, e^{2/(n1)}) 内 lam x0[:-1] / x0[1:] lo, hi np.exp(-2 / (n 1)), np.exp(2 / (n 1)) if not ((lam lo) (lam hi)).all(): # 不满足时做平移保证序列全正且级比落回区间 x0 x0 (lo * x0[1] - x0[0]) 1e-6 # 2) 一次累加生成与紧邻均值 x1 np.cumsum(x0) z1 0.5 * (x1[1:] x1[:-1]) B np.column_stack([-z1, np.ones(n - 1)]) Y x0[1:].reshape(-1, 1) a, b np.linalg.lstsq(B, Y, rcondNone)[0].ravel() # 3) 时间响应式还原出原始序列的拟合值 preds np.array([(x0[0] - b / a) * np.exp(-a * k) * (1 - np.exp(a)) for k in range(1, n 1)]) return a, b, predsa是发展系数b是灰作用量a的绝对值越小说明序列变化越平缓。发展系数绝对值超过 0.3 时预测快速衰减、超过 0.5 时短期预测都不可信这是硬边界。平移操作会改变序列的绝对水平还原时记得再减回去我在这上面栽过一次预测值整体偏高一截。4.2 ARIMA月度数据的季节性和差分阶数怎么定月度序列带有明显的年内周期枯水期和丰水期的浓度差异可能比年际变化还大。定阶的起点是平稳性检验然后按 ACF/PACF 或信息准则挑参数。from statsmodels.tsa.stattools import adfuller from statsmodels.tsa.arima.model import ARIMA def fit_arima(series, p1, q1): s series.dropna() d 1 if adfuller(s)[1] 0.05 else 0 # 不平稳就差分一阶 model ARIMA(s, order(p, d, q)).fit() print(model.summary().tables[1]) # 看 ar/ma 系数是否显著 return modeld不要盲目取大差分会吃掉趋势信息本来在改善的序列差分两次之后就变成噪声了。月度数据只有几十个点季节性差分m12会切掉一整年样本除非数据跨度足够长否则宁可把季节性信息放进外生变量或者环比特征里别硬上 SARIMA。定阶之后一定要看残差Ljung-Box 检验的 p 值大于 0.05 才说明残差里没有再可提取的结构。4.3 LSTM 什么时候值得上样本量红线与滑窗构造神经网络的滑窗构造是标准动作把序列切成定长输入和单点输出。def make_windows(arr, lookback): X, y [], [] for i in range(lookback, len(arr)): X.append(arr[i - lookback:i]) y.append(arr[i]) return np.array(X)[..., None], np.array(y)我的经验红线是时间步少于 100 条时LSTM 只做对照模型不做主模型。lookback一般取 6 或 12对应半年和一年周期隐藏单元数控制在 16~32加一层 Dropout 和早停否则几个 epoch 就把训练集背下来了。归一化一定要用训练段的均值和方差去变换验证段用全序列统计量做归一化等于把未来信息泄露给了模型回测指标会漂亮得离谱。4.4 滚动回测别用“拟合得好”证明预测得准拟合优度不能说明任何问题一个能完美拟合历史的模型可以完全不预测未来。判断标准只有一个滚动向前验证。def rolling_backtest(series, fit_predict, horizon6): errs [] for t in range(len(series) - horizon): train series[:t 1] yhat np.asarray(fit_predict(train, horizon), dtypefloat) true np.asarray(series[t 1:t 1 horizon], dtypefloat) errs.append(np.abs((yhat - true) / true)) return float(np.mean(errs)) # 平均绝对百分比误差滚动回测的每一步只能用当期及之前的数据训练horizon取 6 表示向前预测半年。指标上我一般同时报 MAPE 和 GM(1,1) 的后验差比C S_resid / S_origC 0.35为优、 0.5为合格、 0.65只能算勉强。两个指标方向不一致时以滚动回测为准因为它更接近真实使用场景。方法适用样本量优势主要风险GM(1,1)4~20 期不需要平稳性假设短期稳长期发散级比不满足时误差大ARIMA40 期以上可解释、置信区间明确对结构性突变无能为力LSTM100 期以上能捕捉非线性与多变量耦合小样本过拟合可复现性差5. 避坑与排查五个让结论直接反过来的细节这一段是我在反复重跑这套流程后整理出来的排查清单每一条都真实出现过而且从输出结果上很难一眼看出问题。现象一综合指数算出来上游断面比下游还脏。原因几乎总出在溶解氧的方向上。溶解氧是唯一越大越优的指标如果和氨氮、高锰酸盐指数一起做正向归一化溶解氧越高的断面反而被算成污染越重。解决方式是先做同向化对溶解氧取S/C或者反向极差归一化并且在权重矩阵的注释里写死“已同向化”避免下一次改代码时又被覆盖掉。现象二预测曲线是一条水平直线方差接近于零。多半来自缺测填补用了全表均值。均值填补会消灭时间维度上的方差模型学到的就是“永远等于平均值”。排查方式是打印填充前后的序列标准差如果填充后方差下降超过 30%说明填补策略吃掉了趋势改成按断面分组插值或者干脆把缺测月份从训练集中剔除。现象三GM(1,1) 预测出负浓度或者第十年直接爆炸。前者说明级比检验没做或者平移量给错了后者说明发展系数的绝对值过大模型已经把衰减外推成指数级。解决方式是先跑级比检验不满足就做平移拟合完立刻看a绝对值超过 0.5 就把预测区间压到三年以内并在报告里如实写明模型的短期预测属性。现象四达标率算成 100%但报告结论写水质在恶化。这是评价口径不一致造成的。达标率通常按 III 类标准统计而水质类别序列可能包含大量 IV 类断面。排查方式是统一口径达标率的分母是全部断面月数分子是类别号不超过 3 的记录数同时把类别转移概率一并列出。两者不一致时很可能是部分断面在 III 类和 IV 类之间反复横跳这本身就是值得写进结论的现象。现象五污染源反推的排放量量级差了一千倍。单位换算。浓度单位 mg/L、流量单位 m³/s、目标单位 kg/d中间必须有86400/1000这个系数。另一个隐蔽原因是把衰减项漏掉直接用下游浓度减上游浓度结果把所有自然降解都算成了排放量。排查方式是先做量纲检查再用已知无排放的对照区间验证delta是否接近零这一步能筛掉绝大部分低级错误。6. 用类别转移矩阵给预测结果做交叉验证模型跑完最后一步我习惯加一个不看代码只看结果的验证把每个断面的月度水质类别排成序列统计类别之间的转移频率得到一张 6×6 的马尔可夫转移矩阵。import numpy as np def transition_matrix(sequences, n_class6): sequences: 多条类别序列元素取值 1~n_class P np.zeros((n_class, n_class)) for seq in sequences: for a, b in zip(seq[:-1], seq[1:]): P[int(a) - 1, int(b) - 1] 1 row_sum P.sum(1, keepdimsTrue) row_sum[row_sum 0] 1 # 空行保护避免除零 return P / row_sum P transition_matrix(sequences) init np.array([0.1, 0.2, 0.3, 0.25, 0.1, 0.05]) # 当前类别分布 future np.linalg.matrix_power(P, 12) init # 12 个月后的分布 print(future.round(3))这个矩阵的价值在于它给出了一条完全独立的基线。P的对角线元素是对应类别保持不变的稳定度如果 III 类的自持概率只有 0.3 出头说明断面在类别之间频繁跳动任何单点预测的置信区间都应该放宽。matrix_power(P, 12)给出一年后的稳态倾向把它和 GM(1,1) 或者 ARIMA 的预测方向对照两个方向一致结论可以放心写方向矛盾先去查数据口径和缺测而不是急着换更复杂的模型。矩阵里有明显非对角聚集的行往往对应真实的污染事件月份这些位置值得回到原始表逐条核对我在一次复核中就是从这类聚集里发现某断面的异常值其实是重复录入造成的。做完这一步评价的等级、预测的曲线、转移的稳定度三者在同一张图上对得上整份工作才算闭环。希望帮到你。本文还有配套的精品资源点击获取