简介这份《长江水质的评价和预测》数学建模文档系统展示了针对长江水质评价与预测问题的完整建模思路。文档基于国标地表水环境质量标准围绕五个子问题构建模糊综合评价模型、污染源判别模型与灰色系统预测模型对长江近两年水质进行定量评价识别高锰酸盐和氨氮主要污染源地区并对未来十年水质河长比例与所需污水处理量做出预测。资源共1个doc文档压缩包约595KB内容包含问题重述、模型假设、符号说明、模型求解与结果表格结构完整。适合数学建模竞赛备赛者、环境类专业学生及水质数据分析人员参考从数据预处理到模型构建与结果解读均有清晰推导可作为完整赛题方案学习模糊综合评价与灰色预测的实际应用。已有348人学习下载。1. 模糊综合评价长江水质定级的第一步不是查表长江水质的评价和预测本质上是一条「定级-溯源-预测-配额」的数据链路先用模糊综合评价把 17 个观测站点的水质类别定出来再用一维稳态水质模型反推高锰酸盐和氨氮的主要排放源然后用 GM(1,1) 灰色系统预测未来 10 年各级水质河长比例最后把比例换算成每年必须处理的污水总量。这套思路适合正在做环境数据分析、有小样本时序预测需求或者想复现数学建模完整流程的工程师——它解决的不只是「水质属于几类」的分类问题而是从监测数据到治理决策的完整推算过程。我拆这份资源时最大的感触是水质定级这类边界模糊的问题硬规则往往不如概率分布可靠而模糊综合评价恰好是处理这种不确定性的成熟工具。2. 模糊综合评价的隶属函数设计与权重标定2.1 为什么排除 PH 值再做模糊评价国标 GB3838-2002 给出了 4 个主要项目标准限值乍看可以直接按浓度查表定级。但仔细看限值表会发现一个反直觉的事实PH 值在Ⅰ类到劣Ⅴ类的标准限值全部是 6~9没有任何区分度。而原始数据里 17 个站点 28 个月的 PH 监测值全部落在 6~9 区间内也就是说无论水质好坏这一项都不会把水质等级向下拉。所以数据预处理阶段就直接把 PH 从评价因子集里拿掉只保留溶解氧DO、高锰酸盐指数CODMn和氨氮NH3-N三个指标。这里的经验是做多指标综合评价时先看每个指标的等级区分度——如果某个指标在所有等级上的限值几乎一致它只会稀释权重不会贡献判别信息应当剔除。这比无脑套四指标模型更贴近实际决策场景。另外要注意 DO 和其他两个指标的优化方向相反DO 是越大越好以最大值为最优CODMn 和 NH3-N 是越小越好以最小值为最优。这个方向差异在后续权重计算和隶属函数设计里要分别处理不能用一个公式通吃。2.2 权重归一化与正逆向指标处理权重的设计逻辑是「超标越多权重越大」。对以最小值为最优的指标CODMn、NH3-N用实测值与该指标各级标准最大值之比作为权重基础对以最大值为最优的指标DO则反过来用最大值与实测值之比。公式化表达为正向指标越大越好w_i Smax_i / S_i逆向指标越小越好w_i S_i / Smax_i其中 S_i 是第 i 个指标的实测均值Smax_i 是该指标各级标准中的最大标准值。算完三个权重后做归一化让权重之和等于 1。注意 DO 的公式里实测值在分母所以 DO 浓度越低权重越大符合「溶解氧低说明污染重」的直觉。以四川攀枝花站为例DO 实测均值 9.15CODMn 为 2.43NH3-N 为 0.18。代入公式再归一化后权重向量大体会落在 W [0.21, 0.43, 0.36] 附近具体数值随 Smax 取值略有浮动。氨氮浓度虽低但因为标准限值区间很窄归一化后权重并不小。2.3 用 Python 实现梯形隶属函数隶属函数是模糊综合评价的核心。水质等级之间不是硬边界比如Ⅲ类和Ⅳ类之间没有一条「断崖」而是渐变过渡。这里采用梯形隶属函数把每个指标在每个等级上的隶属度算出来组成模糊关系矩阵 R。import numpy as np # 国标限值DO / CODMn / NH3-N单位 mg/L # 行索引 0..5 对应 I, II, III, IV, V, 劣V 类 limits { DO: [7.5, 6, 5, 3, 2, 0], CODMn: [2, 4, 6, 10, 15, np.inf], NH3-N: [0.15, 0.5, 1.0, 1.5, 2.0, np.inf] } def assign_level(value, name): 对单个指标值做等级分配用于后续隶属度计算的初判 lv limits[name] if np.isinf(lv[-1]) and value lv[-2]: return len(lv) - 1 for i in range(len(lv) - 1): if value lv[i]: return i return len(lv) - 1 def trapezoid_membership(value, low, high): 梯形隶属度落在 [low, high] 内为 1两侧线性过渡到 0 if value low: return max(0.0, 1 - (low - value) / max(low, 1e-6)) if value high: return max(0.0, 1 - (value - high) / max(high, 1e-6)) return 1.0这段代码的要点assign_level 先做初判确定该指标最优落在哪个等级区间trapezoid_membership 负责在区间边界做线性过渡。参数 low 和 high 对应国标里相邻等级的标准限值比如 CODMn 在Ⅱ类对应区间是 2~4实测 3.2 时对Ⅱ类的隶属度就是 1.0同时它对Ⅲ类4~6的隶属度按线性衰减从而保留「接近边界」的模糊信息。np.inf 用在劣Ⅴ类上限避免数值溢出。有了隶属函数后对 17 个站点、3 个指标分别计算 6 个等级的隶属度组成 3×6 的模糊关系矩阵 R再与归一化权重向量 W 做模糊合成得到多指标模糊概率综合评价矩阵 M W ∘ R。合成算子常见取加权平均型def fuzzy_compose(weights, R): 加权平均型模糊合成M[j] sum(W[i] * R[i][j]) return np.dot(weights, R)逻辑说明weights 是长度为 3 的一维数组R 是 3×6 的二维数组np.dot 对每个等级 j 做加权求和输出长度为 6 的 M 向量。M[j] 表示该站点水质归属于第 j 类的综合模糊概率。定级时取最大隶属度原则——M 中最大值对应的索引就是水质类别同时计算综合质量系数 M 与等级向量 [1,2,3,4,5,6] 的内积这个系数越高说明水质越差用于跨站点横向排名。2.4 17 个站点的模糊概率与综合质量系数按上述流程跑完全部站点抽取几个有代表性的结果站点Ⅰ类Ⅱ类Ⅲ类Ⅳ类Ⅴ类劣Ⅴ类评价等级综合质量系数排名四川攀枝花0.9100.8800.0000.0000.0000.000Ⅰ2.6688湖北丹江口胡家岭1.0000.0000.0000.0000.0000.000Ⅰ1.00017四川乐山岷江大桥0.0000.2901.0000.0000.0000.000Ⅲ3.5802江西南昌滁槎0.3700.7200.3900.0000.1100.170Ⅳ4.5301江苏南京林山1.0000.1400.0000.0000.0000.000Ⅰ1.28016注意两个细节一是丹江口和南京林山虽然都定级为Ⅰ类但综合质量系数不同说明系数提供了等级之外的排序信息适合做细粒度污染排名二是江西南昌滁槎对不同等级都有非零概率最大隶属度原则把它判为Ⅳ类但它的Ⅱ类概率其实有 0.72这种「跨界」信号在传统查表法里会被直接丢弃。这也是模糊综合评价相对硬分类的核心优势——保留不确定性而不是强行二值化。从结果统计看17 个站点中 23.53% 为Ⅰ类、70.59% 为Ⅱ类、5.88% 为Ⅲ类没有站点落入Ⅳ类以下。整体水质偏好但支流站点乐山岷江大桥、南昌滁槎明显比干流差这个结论直接为后续污染溯源提供了方向。3. 污染源判别模型一维稳态水质方程与排放量反推3.1 浓度沿程衰减与质量流量换算模糊综合评价能回答「哪里水质差」但不能回答「污染物是谁排的」。要定位污染源需要利用长江干流 7 个观测站点的距离、水流量和水流速数据建立污染物沿程衰减关系。这里采用一维稳态河流水质模型污染物浓度随距离按指数规律衰减即c_x c_0 · e^(-K · x / v)其中 c_0 是上游起点浓度K 是降解系数单位 1/天x 是站点间距离kmv 是江水流速km/天。原始资源里提到 K 通常介于 0.1~0.5取 0.2 作为默认值。这里的物理含义是污染物在流动过程中通过物理降解、化学降解和生物降解自然消减从上游流到下游时浓度已经变小所以下游站点测到的污染物不全是本地排的有相当一部分是上游残留。质量流量换算也很关键。站点浓度单位是 mg/L要换算成每秒通过断面的污染物质量g/s需要乘上断面流量M_k c_k · Q_k其中 Q_k 是第 k 个站点的水流量m³/s换算时注意单位统一1 mg/L 1 g/m³所以 M_k c_k (mg/L) × Q_k (m³/s) 的结果单位就是 g/s。这一步不做后面所有排放量比较都是空中楼阁。3.2 由站点总量分离本地排放量核心假设是第 k 个站点的污染物总量 M_k 等于「上游流入本站的残留量」加「本站自身排放量」。而上游流入本站的残留量可以用第 k-1 个站点的质量流量按一维稳态模型衰减折算到本站位置M_k^up M_{k-1} · e^(-K · Δx_k / v_k)于是本站自身排污量m_k M_k - M_k^up需要注意边界条件起始站点四川攀枝花没有上游所以它的自身排放量就等于 M_1 本身。这就是假设 2「干流上起始站点所含的污染物只是其本身排放的污染物」的含义。用 Python 实现如下import numpy as np # 7 个干流站点基本数据示意结构实际以附件3为准 stations [攀枝花, 朱沱, 宜昌南津关, 岳阳城陵矶, 九江河西水厂, 安庆皖河口, 南京林山] distances [0, 422, 712, 358, 132, 164, 300] # km相邻站间距 flow_speed 1.2 # m/s干流平均流速简化取值 K 0.2 # 1/day降解系数 Q [4200, 6800, 7200, 9100, 9800, 10500, 12000] # m^3/s示意流量 CODMn_conc [2.1, 2.3, 2.8, 3.1, 2.9, 3.4, 3.6] # mg/L def compute_local_discharge(conc, distances, K, v, Q): v_km_day v * 86.4 # m/s - km/day M np.array(conc) * np.array(Q) # g/s断面质量流量 local np.zeros_like(M) local[0] M[0] for i in range(1, len(M)): delta_x distances[i] residual M[i-1] * np.exp(-K * delta_x / v_km_day) local[i] M[i] - residual return local discharge compute_local_discharge(CODMn_conc, distances, K, flow_speed, Q)逻辑说明v_km_day 把流速从 m/s 换算成 km/天因为 K 的单位是 1/天x 的单位是 km必须保证量纲一致。np.exp 计算的是上游污染物衰减到本站的剩余比例local[i] 就是本站实际排入长江的污染物质量流量。如果 local 出现负值通常是浓度测量误差或 K 取值偏大导致衰减过多这时要把负值按 0 处理并检查该站点的数据可信度。3.3 排放量结果与污染源判定按上述流程计算高锰酸盐和氨氮在 7 个站点上的自身排放量站点高锰酸盐排放量 (g/s)氨氮排放量 (g/s)四川攀枝花66974464重庆朱沱2915827337湖北宜昌南津关45284268湖南岳阳城陵矶5435180845江西九江河西水厂3593459950安徽安庆皖河口5628173298江苏南京三江营155769899判定规则是高锰酸盐排放量超过 40 kg/s、氨氮排放量超过 3 kg/s 的站点列为对应污染物的主要污染源。按此规则高锰酸盐的主要污染源依次是安徽安庆、湖南岳阳、江西九江、江苏南京、湖北宜昌氨氮的主要污染源依次是湖南岳阳、安徽安庆、江西九江、重庆朱沱、湖北宜昌。这里有个单位转换的坑要提醒模型计算结果是 g/s判定阈值是 kg/s两者差 1000 倍很多初做这个项目的人会在这里对不上数。能源头数据发现攀枝花作为上游起点虽然浓度不是最低但因为流量小绝对排放量反而是 7 个站点里最少的。这说明污染源判别必须用质量流量而非浓度——浓度高的小支流在总污染负荷里可能远不如浓度中等的大干流。这个思路可以平移到任何「多个监测断面、需要区分本地排放和上游输入」的水环境场景。4. GM(1,1) 灰色预测级比检验与小样本时序建模4.1 为什么小样本选灰色模型而不是 ARIMA问题三是预测 2005 年到 2014 年各级水质河长比例但手头只有 1995 到 2004 年一共 10 个年度数据点。对 ARIMA 或 LSTM 这类统计/深度模型10 个样本根本不够支撑参数估计而灰色系统理论的核心优势就是「少数据、贫信息」——它不需要大样本也不要求数据满足正态分布或平稳性假设。GM(1,1) 的基本思路把原始数列做一次累加生成1-AGO把原来波动较大的数据序列变成单调递增序列然后用一阶微分方程拟合这条累加序列再通过后减运算还原出预测值。这听起来有点绕但本质上是「先累加消噪再拟合趋势再差分还原」。4.2 级比检验与平移变换在建模之前必须先做级比检验否则模型可能不适用。设原始数据列为 X⁰ [x⁰(1), x⁰(2), ..., x⁰(n)]级比为λ(k) x⁰(k-1) / x⁰(k)如果所有级比都落在可容覆盖区间 (e^(-2/(n1)), e^(2/(n1)))则数据序列可以直接用 GM(1,1)。n10 时可容覆盖区间约为 (0.833, 1.199)。原始资源里对 1995~2004 年枯水期各级水质河长比例做检验时原始序列的级比范围是 0.89~1.18勉强落在区间内但余量很小部分类别如Ⅰ类水比例接近 0 的那些年份会越界因为比例数据大量集中在低位相邻年份比值波动剧烈。处理办法是平移变换对所有数据加一个常数 C使得新序列的级比落入可容区间。这里取 C50。平移后枯水期级比范围变成 0.89~1.18与可容区间 (0.833, 1.199) 对照全部落在区间内。注意平移变换不会改变序列的相对趋势只改变幅度所以预测后再减回 C 即可得到真实比例。4.3 最小二乘求解与后减还原GM(1,1) 的完整 Python 实现如下import numpy as np def gm11(x0, n_pred10): x0: 原始数据列长度 n n_pred: 预测步数 返回预测值序列含原始区间内的拟合值 x0 np.asarray(x0, dtypefloat) n len(x0) # 1-AGO 累加生成 x1 np.cumsum(x0) # 构造 B 矩阵和 Y 向量 B np.column_stack([-0.5 * (x1[:-1] x1[1:]), np.ones(n - 1)]) Y x0[1:] # 最小二乘估计参数 [a, u] theta np.linalg.inv(B.T B) B.T Y a, u theta[0], theta[1] # 累加序列的拟合值 x1_fit np.zeros(n n_pred) x1_fit[0] x1[0] for k in range(1, n n_pred): x1_fit[k] (x0[0] - u / a) * np.exp(-a * k) u / a # 后减运算还原 x0_fit np.zeros(n n_pred) x0_fit[0] x0[0] x0_fit[1:] np.diff(x1_fit) return x0_fit, a, u参数含义a 是发展灰数反映序列的增长或衰减趋势a 的绝对值越大趋势越强u 是内生控制灰数相当于系统的外部驱动项。B 矩阵的每一行是累加序列相邻两项的均值取负这是灰色模型对微分方程离散化的标准做法Y 向量是原始序列从第二项开始的值。np.linalg.inv(B.T B) B.T Y 就是最小二乘闭式解不需要迭代。调用时要注意原始数据如果有零值累加序列会出现连续相等段导致 B 矩阵奇异。这时要用平移变换或对零值做小幅常数偏移避免矩阵求逆报错。另外后减运算用的是 np.diff它会把累加序列的相邻差作为原始序列的拟合值这一步丢失了第一个点的还原所以手动把 x0_fit[0] 赋值为 x0[0]。4.4 残差与级比偏差双重验证单有预测值不足以说明模型可靠还要做两类检验。残差检验计算相对误差 ε(k) (x⁰(k) - x̂⁰(k)) / x⁰(k)要求多数点 ε(k) 0.1级比偏差检验先由原始数据算级比 λ(k)再结合发展系数 a 算理论级比两者偏差小于 0.1 认为合格。以枯水期Ⅰ类水比例为例2005 到 2014 年的相对误差全部落在 0.00~0.07 之间级比偏差绝对值最大 0.15大部分在 0.1 以内模型通过验证。预测得到的水文年各级水质河长比例年份Ⅰ类Ⅱ类Ⅲ类Ⅳ类Ⅴ类劣Ⅴ类20050.027.9117.266.313.217.4020060.028.7515.155.442.727.4220070.029.6013.054.572.247.4320080.030.4410.973.701.757.4520090.031.308.912.831.277.4620100.032.156.881.970.797.4820110.033.014.861.110.317.4920120.033.872.860.260.007.5120130.034.730.890.000.007.5220140.035.600.000.001.877.54数据背后两个信息一是Ⅱ类水比例逐年上升说明总体水质在变好二是劣Ⅴ类比例从 7.40 缓慢爬升到 7.54说明重污染河段在顽固扩张。前者解释了为什么干流站点评价结果偏好后者支撑了「需要严格控制污水排放」的结论。如果只做均值或线性回归很难同时抓住这两种不同速率的趋势这是灰色模型在这个场景里真正起作用的地方。5. 从河长比例到治污配额污水量测算与模型边界5.1 六类水压缩成三类处理量公式化问题四给的条件是未来 10 年内Ⅳ类和Ⅴ类水的比例控制在 20% 以内且没有劣Ⅴ类水。原始资源把六类水按用途压缩成三类——A 类可饮用水含Ⅰ、Ⅱ、Ⅲ类、B 类可用水含Ⅳ、Ⅴ类、C 类废水劣Ⅴ类。需要处理的污水量由两部分组成B 类超过 20% 的部分全部处理C 类全部处理。写成公式L max(0, P_B - 20%) × Q_total P_C × Q_total其中 P_B、P_C 是 B、C 类水质占河长的百分比Q_total 是当年长江总水量。用灰色预测分别得到 2005~2014 年的总流量和各级水质比例代入即得每年处理量。实测结果是逐年递增2005 年 298.9 亿吨2014 年 470.8 亿吨十年涨了 57%。这个单调上升曲线直接说明在不加干预的前提下污水治理压力只增不减。5.2 预测结果如何落入管理指标这个测算过程的价值在于把「水质类别」这种定性指标转化成了「每年必须处理多少亿吨污水」的定量配额。对环保部门来说决策依据不是某段河属于Ⅳ类而是「要在 20% 红线内明年至少新增 xx 亿吨处理能力」。类似的逻辑也能迁移到别的场景——工业园区排放总量控制、污水处理厂扩容规划、流域横向生态补偿金额计算本质都是「先预测比例再按比例配额」的结构。5.3 模型的适用边界复盘时有三点边界值得注意。第一GM(1,1) 的外推步数最好不要超过样本量的一半用 10 年数据外推 10 年已经接近极限预测区间越往后可信度越低尤其是 2014 年的预测值只能当趋势参考不能当精确值用。第二降解系数 K 取常数 0.2 是工程简化实际中 K 随水温、季节、流量变化很大如果要做更精细的溯源需要分段标定 K 值。第三模糊综合评价的最大隶属度原则在相邻等级概率接近时会误判比如江西南昌滁槎同时对Ⅱ类和Ⅳ类都有较高隶属度这时只看评级会丢失信息建议配合综合质量系数使用。把这三点记在项目文档里比直接用模型结论写报告要稳妥得多。本文还有配套的精品资源点击获取