长江水质评价与预测：从数据处理到模型选择的完整实践指南

简介：长江水质评价与预测数学建模文档，内容完整，适合数学建模竞赛参赛者、环境科学与水利工程专业学生及从事水质分析的研究人员使用。文档围绕四个核心问题展开：基于模糊综合评价法对长江近两年水质进行定级，构建主要污染源判别模型确定高锰酸盐与氨氮的重点排放区域，利用灰色系统模型预测未来十年各类水质河长比例，并计算为控制Ⅳ、Ⅴ类水比例所需处理的污水总量。全文涵盖模型假设、符号说明、问题重述、求解步骤及结果表格，并附有对环保部门的参考性建议，既可作为数学建模论文写作范例，也可用于学习模糊综合评价、灰色预测等方法的环境应用。资源为1个doc文档，压缩包大小595KB，已有349人学习。文档结构清晰，从问题重述到模型求解与结果展示一应俱全，适合需要快速理解水质评价与预测建模思路的读者下载研读。

1. 长江水质的评价和预测：先把“等级”和“趋势”拆成两件事

拿到“长江水质的评价和预测”这份文档的人，十有八九第一反应是打开附件看数据，然后直接上模型。我第一次做的时候也是这个顺序，结果评价出来的等级和常识对不上，预测曲线平滑得像画上去的。后来想明白一件事：长江水质的评价和预测是两条完全不同的技术链路。评价是一次快照，处理同一时间断面上多个指标的合成问题，核心矛盾在于权重怎么定、最差指标要不要一票否决；预测是时间序列，处理同一断面多个指标在几十个月里的演变，核心矛盾在于样本量太少、季节性太强、缺测太多。把这两条链路混在一起想，就会陷入“用评价的等级序列直接做预测”的误区——等级是离散标签，跳变多、信息量低，拿它喂给时间序列模型，残差一定难看。这套东西适合环境类、水利类的建模竞赛选手、以及要做流域水质报表的工程技术人员，只要能拿到断面监测原始表，就能按下面的路径一步步复现。

2. 把附件里的监测表变成可计算的矩阵：断面、月份、指标三轴对齐

打开附件通常是几张排版不规整的表格，列是断面名称，行是监测时间，中间夹着溶解氧、高锰酸盐指数、氨氮这些指标，有的表头还带单位后缀。真正开始算之前有三件事必须先做完：把评价标准落到代码里、把宽表摊成长表、把缺测和检出限处理掉。这三步做歪了，后面无论用多花哨的模型，结论都是错的。我在第一次做的时候跳过了第二步，直接在宽表上循环遍历列名，代码又长又容易在下一次换数据格式时崩掉，血泪经验。

2.1 先把 GB 3838 的五类限值写成一张能查的表

评价的基础是《地表水环境质量标准》里的限值分级。不同指标的分界值不一样，而且方向不一致——大部分指标越小越好，溶解氧恰恰相反，越大越好。把限值硬编码进if-else是灾难的开始，正确做法是做成结构化字典，后续查表统一调用。

这张表我最常用来做两件事：一是把每个断面的实测值映射成 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["监测时间"], errors="coerce") 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指定的是保持不动的列，选错了会把断面名也当成指标搅进去。errors="coerce"让无法解析的时间变成NaT而不是直接抛异常，方便一次性看到所有格式问题。半检出限取值是最常见的处理方式，但如果某一指标的检出限占标准限值的比例很大，这个近似会系统性抬高整体浓度，必须在限制条件里说明。最后那行print是黑匣子，我几乎每个项目都会顺手跑一遍，最大值为负、最小值比检出限还低一个量级，这些异常在这一步就能拦住。

2.3 缺测怎么补：按断面插值，别用全局均值

缺测在断面月序列里非常普遍，尤其是枯水期某些断面停测。最常见的错误做法是用全表均值填补，这会把时间趋势抹平，后面预测出来的曲线自然就是一条水平直线。

long = long.sort_values(["断面", "指标", "监测时间"]) # 每个断面、每个指标各自插值，时间轴上做线性内插 long["测值_填充"] = ( long.groupby(["断面", "指标"])["测值"] .transform(lambda s: s.interpolate(method="linear", limit=2, limit_direction="both")) ) # 覆盖率体检：整段缺失的指标不要插，直接标记为不可评价 coverage = long.groupby(["断面", "指标"])["测值"].apply(lambda s: s.notna().mean()) print(coverage[coverage < 0.7])

limit=2表示最多连续补两期，超过两期的空洞说明监测体系本身有问题，硬补出来的数字支撑不了一个评价结论。limit_direction="both"处理序列首尾，否则开头结尾的缺口永远是空的。覆盖率阈值 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 # side='left' 对应“<= 界值即归入该类”的语义 return int(np.searchsorted(bounds, value, side="left")) + 1 def class_dissolved_oxygen(value): """溶解氧越大越好，反向查表。""" if pd.isna(value): return np.nan for i, low in enumerate([7.5, 6, 5, 3, 2], start=1): if value >= low: return i return 6

np.searchsorted的side参数是这里唯一的玄学点，side="left"在值恰好等于界值时返回前一段，符合标准里“小于等于”的表述；换成"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, std3=STD3): """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, keepdims=True) 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/s；C: 浓度 mg/L；k: 综合衰减系数 1/d""" decay = np.exp(-k_per_day * L_km / u_km_per_day) delta = C_down - C_up * decay # 扣除自然衰减后真正新增的浓度 return delta * Q * 86400 / 1000

k是最需要标定的参数，氨氮常见在 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, dtype=float) n = len(x0) # 1) 级比检验：lambda(k) = x(k-1)/x(k) 需落在 (e^{-2/(n+1)}, e^{2/(n+1)}) 内 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, rcond=None)[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, preds

a是发展系数，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, p=1, q=1): 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 model

d不要盲目取大，差分会吃掉趋势信息，本来在改善的序列差分两次之后就变成噪声了。月度数据只有几十个点，季节性差分m=12会切掉一整年样本，除非数据跨度足够长，否则宁可把季节性信息放进外生变量或者环比特征里，别硬上 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, horizon=6): errs = [] for t in range(len(series) - horizon): train = series[:t + 1] yhat = np.asarray(fit_predict(train, horizon), dtype=float) true = np.asarray(series[t + 1:t + 1 + horizon], dtype=float) errs.append(np.abs((yhat - true) / true)) return float(np.mean(errs)) # 平均绝对百分比误差

滚动回测的每一步只能用当期及之前的数据训练，horizon取 6 表示向前预测半年。指标上我一般同时报 MAPE 和 GM(1,1) 的后验差比C = S_resid / S_orig，C < 0.35为优、< 0.5为合格、< 0.65只能算勉强。两个指标方向不一致时，以滚动回测为准，因为它更接近真实使用场景。

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_class=6): """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, keepdims=True) 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 的预测方向对照：两个方向一致，结论可以放心写；方向矛盾，先去查数据口径和缺测，而不是急着换更复杂的模型。矩阵里有明显非对角聚集的行往往对应真实的污染事件月份，这些位置值得回到原始表逐条核对，我在一次复核中就是从这类聚集里发现某断面的异常值其实是重复录入造成的。做完这一步，评价的等级、预测的曲线、转移的稳定度三者在同一张图上对得上，整份工作才算闭环。希望帮到你。

本文还有配套的精品资源，点击获取