面源污染SWAT模型实战：从数据准备到负荷出图的完整路径

简介：这份《面源污染SWAT模型讲义》面向水文、环境与农业管理领域初学者及科研人员，系统讲解分布式水文模型SWAT的原理与应用。资源以PPT形式呈现，压缩包内共1个文件，约564KB，内容涵盖模型简介、水文循环模块划分及泾河流域应用实例。讲义从气候、陆地水文与水面水文三部分拆解水循环过程，说明日降水、气温、太阳辐射、下渗、蒸散发、地表径流与河道汇流等环节的模拟思路，并给出SCS曲线数、水库调蓄等关键方法。应用部分以泾河流域为例，基于1996至2000年平凉站气象数据、2000年土地利用与1:400万土壤信息，将流域划分为92个子流域，比较冬小麦、草地、乔木、灌木等植被类型的水土流失差异，并分析坡度与土壤流失的相关性，指出30度以上坡度呈明显正相关。目前已有985人学习，适合希望快速建立SWAT建模框架、理解非点源污染模拟流程的读者参考。

1. 面源污染 SWAT 模型讲义：从数据到出图，一套能跑通的落地路径

面源污染这件事，麻烦就麻烦在它不像点源那样有个明确的排口。农田里的氮磷、畜禽养殖的径流、城镇地表的冲刷，都是跟着降雨和地形慢慢“渗”进河网的，你很难指着某个位置说这就是污染源。SWAT 模型（Soil and Water Assessment Tool）正是为这种“漫散型”污染设计的分布式水文水质模型，它把流域切成若干子流域和水文响应单元，逐个算产流、产沙、氮磷迁移，最后汇到出口断面。这套讲义要解决的核心问题很具体：拿到一个流域的 DEM、土地利用、土壤和气象数据后，怎么把 SWAT 跑起来、参数怎么定、结果怎么验证、面源污染负荷怎么出图。适合做流域水环境管理、农业面源核算、环评水文专题的从业者，也适合刚接触分布式模型的研究生。下面按“建库—率定—出图—避坑”的顺序，把每一步拆到能照着操作。

2. 建模前的数据准备：四类输入怎么凑齐、怎么对齐

SWAT 对输入数据的格式和坐标系非常敏感，很多新手第一次跑失败不是模型问题，而是数据没对齐。这一章把四类必需输入——地形、土地利用、土壤、气象——的获取、处理和常见坑讲清楚。

2.1 DEM 的填洼与流向计算：别跳过这一步

DEM 是 SWAT 划分河网和子流域的基础。原始 DEM 里通常有洼地（sink），如果不填，水流方向会断，河网就画不出来。常见做法是在 GIS 里先做填洼（Fill Sinks），再算流向（Flow Direction）和流量累积（Flow Accumulation）。

# 以常见 GIS 命令行工具为例，先填洼再算流向 # 填洼：消除 DEM 中的封闭洼地 fill_sinks -i dem_raw.tif -o dem_filled.tif # 计算 D8 流向 flow_direction -i dem_filled.tif -o flow_dir.tif # 计算流量累积，用于提取河网 flow_accumulation -i flow_dir.tif -o flow_acc.tif

逻辑说明：填洼是保证水文连通性的前提，流向决定每个栅格的水往哪走，流量累积值越大说明汇水面积越大，通常取累积值超过某个阈值的栅格作为河道。参数说明：填洼阈值一般默认即可；河网提取阈值需要根据流域面积调，小流域取小值（如 1000），大流域取大值（如 10000），阈值太小河网会碎，太大主河道会断。这一步做完，把填洼后的 DEM 作为 SWAT 的 DEM 输入，不要用原始 DEM。

2.2 土地利用和土壤数据的重分类：索引表是命门

SWAT 不认你原始的土地利用代码，它需要一张索引表（lookup table）把你的分类映射到 SWAT 内置的作物/植被类型和土壤属性。土地利用重分类时，把二级分类合并成 SWAT 能识别的类型，比如水田、旱地、林地、草地、建设用地、水域。

# 用 Python 做土地利用重分类的示意 import rasterio import numpy as np # 读取原始土地利用栅格 with rasterio.open('landuse_raw.tif') as src: landuse = src.read(1) profile = src.profile # 定义映射：原始代码 -> SWAT 代码 # 假设原始代码 1=水田, 2=旱地, 3=林地, 4=草地, 5=建设用地, 6=水域 mapping = {1: 'AGRL', 2: 'AGRR', 3: 'FRST', 4: 'PAST', 5: 'URBN', 6: 'WATR'} # 这里仅示意逻辑，实际需按 SWAT 索引表逐类赋值 reclass = np.zeros_like(landuse) for raw_code, swat_code in [(1, 1), (2, 2), (3, 3), (4, 4), (5, 5), (6, 6)]: reclass[landuse == raw_code] = swat_code # 写出重分类结果 profile.update(dtype=rasterio.int16, nodata=0) with rasterio.open('landuse_swat.tif', 'w', **profile) as dst: dst.write(reclass, 1)

逻辑说明：这段代码把原始土地利用栅格的像元值按映射关系替换成 SWAT 能识别的代码。参数说明：映射关系必须和 SWAT 安装目录下的landuse.lookup或你自定义的索引表一致，否则模型会报“未定义土地利用类型”。土壤数据同理，需要准备土壤类型栅格和对应的土壤属性数据库（.sol 文件），属性包括土层厚度、容重、有效含水量、饱和导水率等。常见坑是土壤属性单位不统一，SWAT 默认用公制，如果你拿到的容重是英制，必须换算。

2.3 气象数据的格式与插补：缺测不是小事

SWAT 需要逐日气象数据：降水、最高温、最低温、太阳辐射、风速、相对湿度。前三个是必需，后三个可以用天气发生器生成。气象站数据经常有缺测，直接删掉会导致时间序列不连续，模型报错。

# 气象数据整理成 SWAT 要求的 .txt 格式 # 每行格式：站号 年 月 日 降水 最高温 最低温 # 缺测值用 -99 标记，SWAT 会自动用天气发生器插补 cat pcp_swat.txt | head -5 # 输出示例： # 101 1980 1 1 0.0 5.2 -3.1 # 101 1980 1 2 2.4 4.8 -2.5

逻辑说明：SWAT 的天气发生器（WXGEN）会根据多年统计特征自动填补缺测，但前提是你把缺测标记为 -99，而不是留空或写 0。参数说明：降水单位是 mm，温度单位是摄氏度。如果研究区有多个气象站，建议全部纳入，SWAT 会用泰森多边形或反距离权重法插值到每个子流域。注意气象数据的年份范围要覆盖模拟期和预热期，预热期一般设 2-3 年，让模型内部状态（如土壤含水量）达到平衡。

3. 模型率定与验证：参数怎么调、指标怎么看

SWAT 有几十个参数，但对面源污染负荷敏感的主要是那几个。这一章讲清楚率定的顺序、敏感参数的范围和验证指标的门槛。

3.1 敏感参数筛选：先做一次全局敏感性分析

不要一上来就手动调所有参数。先用 SWAT-CUP 或类似工具做一次敏感性分析，找出对径流和氮磷负荷影响最大的参数。常见做法是拉丁超立方抽样，跑几百次，看 t-stat 和 p-value。

逻辑说明：这张表是面源污染模拟中最常调的参数。CN2 调大，径流增大，泥沙和吸附态污染物输出增加；ALPHA_BF 调小，基流衰减慢，地下径流占比高，氮的淋失路径改变。参数说明：率定时先调水文参数（CN2、ALPHA_BF、GW_DELAY），让径流过程线对上，再调水质参数（SOL_N、PPERCO），让氮磷浓度对上。顺序反了会互相干扰。

3.2 率定指标：R²、NSE 和 PBIAS 的及格线

率定结果不能只看 R²，水文模型更看重 NSE（纳什效率系数）和 PBIAS（百分比偏差）。

# 计算 NSE 和 PBIAS 的示意 import numpy as np def nse(obs, sim): # 纳什效率系数：越接近 1 越好 return 1 - np.sum((obs - sim)**2) / np.sum((obs - np.mean(obs))**2) def pbias(obs, sim): # 百分比偏差：越接近 0 越好，正偏表示模拟偏大 return 100 * np.sum(sim - obs) / np.sum(obs) # 假设 obs 和 sim 是逐日径流序列 obs = np.array([...]) # 实测值 sim = np.array([...]) # 模拟值 print('NSE =', nse(obs, sim)) print('PBIAS =', pbias(obs, sim))

逻辑说明：NSE 衡量模拟值对实测值方差的解释能力，大于 0.5 算及格，大于 0.65 算好。PBIAS 衡量总量偏差，绝对值小于 25% 算可接受，小于 10% 算好。参数说明：率定期和验证期要分开，通常按 7:3 或 6:4 划分。验证期指标可以比率定期略差，但不能差太多，否则说明过拟合。注意面源污染的验证数据往往比水文数据少，很多流域只有几个断面的水质监测，这时候要结合源解析或文献值做合理性判断。

3.3 水质率定的特殊处理：氮磷的滞后与累积

面源污染的氮磷输出不是线性的，一场大雨可能带走几个月累积的养分。率定时要关注累积效应和滞后效应。

# 在 SWAT 的 .hru 文件中查看氮磷输出 # 关注 ORGN、ORGP（有机氮磷）、NSURQ、PSURQ（地表径流中的氮磷） grep -E "ORGN|ORGP|NSURQ|PSURQ" output.hru | head -20

逻辑说明：有机氮磷主要吸附在泥沙上，跟泥沙输出同步；溶解态氮磷跟径流同步。如果模拟的氮磷峰值比实测提前，可能是 CN2 偏大导致径流响应太快；如果峰值偏低，可能是 SOL_N 或 PPERCO 设置偏小。参数说明：水质率定通常比水文率定更难，因为监测数据少、误差大。常见做法是先用文献值定初始范围，再用少量实测数据微调，最后做不确定性分析。

4. 面源污染负荷出图：从模型输出到空间分布图

模型跑完只是第一步，把结果变成能看的图才是交付。这一章讲怎么从 SWAT 输出文件里提取氮磷负荷，做成空间分布图和时间序列图。

4.1 从 output.rch 提取断面负荷

SWAT 的output.rch文件记录了每个子流域出口断面的逐月/逐日负荷。

# 读取 output.rch 并提取总氮负荷 import pandas as pd # 跳过前几行说明，按固定宽度读取 df = pd.read_csv('output.rch', skiprows=9, sep=r'\s+', header=None) # 列名参考 SWAT 文档：RCH, MON, AREA, FLOW_IN, FLOW_OUT, ... ORGN, ORGP, ... df.columns = ['RCH', 'MON', 'AREA', 'FLOW_IN', 'FLOW_OUT', 'EVAP', 'TLOSS', 'SED_IN', 'SED_OUT', 'SEDCONC', 'ORGN', 'ORGP', 'NSURQ', 'SOLP', 'SEDP', 'ORGN_OUT', 'ORGP_OUT', 'NO3', 'NH4', 'NO2', 'SOLP_OUT', 'MINP', 'ORGP_OUT2', 'TOTN', 'TOTP'] # 按子流域汇总总氮 tn_summary = df.groupby('RCH')['TOTN'].sum().reset_index() print(tn_summary.head())

逻辑说明：output.rch是固定宽度格式，用sep=r'\s+'可以自动分割。列名要根据你用的 SWAT 版本核对，不同版本列顺序可能略有差异。参数说明：TOTN 是总氮，TOTP 是总磷，单位通常是 kg/ha 或 kg。如果要做空间分布图，把tn_summary按 RCH 编号关联到子流域矢量图层，再用 GIS 做分级设色。

4.2 空间分布图：用子流域矢量关联负荷

# 用 geopandas 把负荷关联到子流域并出图 import geopandas as gpd import matplotlib.pyplot as plt # 读取子流域矢量 subbasins = gpd.read_file('subbasins.shp') # 关联负荷 subbasins = subbasins.merge(tn_summary, left_on='Subbasin', right_on='RCH') # 出图 fig, ax = plt.subplots(1, 1, figsize=(10, 8)) subbasins.plot(column='TOTN', cmap='YlOrRd', legend=True, ax=ax, legend_kwds={'label': '总氮负荷 (kg)'}) ax.set_title('子流域总氮负荷空间分布') plt.savefig('tn_spatial.png', dpi=300, bbox_inches='tight')

逻辑说明：这段代码把output.rch里的负荷数据按子流域编号合并到矢量图层，然后用颜色深浅表示负荷高低。参数说明：cmap选YlOrRd适合表示污染负荷，颜色越深负荷越高。如果要做年际变化，把MON列按年汇总后再关联。注意子流域编号要和矢量图层的编号一致，不一致就先做字段匹配。

4.3 时间序列图：看氮磷的月际和年际变化

# 画总氮月际变化 df['DATE'] = pd.to_datetime(df[['YEAR', 'MON']].assign(DAY=1)) monthly_tn = df.groupby('DATE')['TOTN'].sum() fig, ax = plt.subplots(figsize=(12, 5)) monthly_tn.plot(ax=ax, color='steelblue') ax.set_ylabel('总氮负荷 (kg)') ax.set_xlabel('时间') ax.set_title('流域总氮负荷月际变化') plt.savefig('tn_timeseries.png', dpi=300, bbox_inches='tight')

逻辑说明：时间序列图能看出氮磷输出的季节性，通常汛期负荷高，非汛期低。参数说明：如果要做趋势分析，可以加一条滑动平均线。注意 SWAT 输出的时间步长要和你的分析尺度匹配，日输出适合看暴雨事件，月输出适合看季节规律。

5. 避坑与排查：面源污染 SWAT 建模的五个血泪教训

这一章是我自己踩过的坑，每条按“现象—原因—解决”写，希望能帮你省点时间。

5.1 模型跑完但河网是断的

现象：SWAT 成功运行，但生成的河网图层只有几段，不连续。原因：DEM 填洼不彻底，或者河网提取阈值设得太大。解决：重新做填洼，检查填洼后的 DEM 是否还有封闭洼地；把河网提取阈值调小，小流域建议从 1000 开始试。

5.2 率定结果 NSE 死活上不去

现象：调了很多次参数，NSE 始终在 0.3 以下。原因：可能是气象数据有系统偏差，或者子流域划分太粗，或者实测径流数据本身有问题。解决：先检查气象数据的单位和缺测处理，再检查子流域数量是否合理（一般 20-50 个为宜），最后核对实测数据的整编方式是否和模型输出一致。

5.3 氮磷负荷模拟值比实测高一个数量级

现象：总氮模拟值远大于实测。原因：SOL_N 初始值设得太大，或者施肥数据输入错误。解决：检查.hru文件里的初始氮浓度，通常土壤初始氮在 0-100 ppm 量级；核对施肥日期、施肥量和施肥类型，SWAT 的施肥单位是 kg/ha。

5.4 模型报错“未定义土地利用类型”

现象：运行时报错，提示某个土地利用代码没有对应的索引。原因：土地利用重分类时漏了某个类别，或者索引表里没有对应项。解决：用 GIS 统计所有土地利用代码，逐一核对索引表，确保每个代码都有映射。

5.5 输出文件太大，打不开

现象：output.hru或output.rch文件几个 G，Excel 打不开。原因：输出时间步长设成了逐日，且模拟年限长。解决：在file.cio里把输出步长改成月或年，或者用 Python 分块读取。如果只需要特定子流域，可以在file.cio里指定输出范围。

6. 进阶技巧：用 Python 批量处理多情景和不确定性分析

单情景跑通只是入门，实际项目往往要做多个情景（比如不同施肥方案、不同土地利用变化）的对比，还要做不确定性分析。这一章讲怎么用 Python 把 SWAT 的批量运行和结果分析自动化。

6.1 批量修改参数并运行

SWAT 的参数文件是文本格式，可以用 Python 批量替换。

# 批量修改 CN2 并运行 SWAT import os import shutil def modify_cn2(project_dir, cn2_value): # 读取原始 .mgt 文件，替换 CN2 mgt_file = os.path.join(project_dir, '000010001.mgt') with open(mgt_file, 'r') as f: lines = f.readlines() # 找到 CN2 所在行并替换（具体行号因文件而异） for i, line in enumerate(lines): if 'CN2' in line: lines[i] = f' {cn2_value:.2f} | CN2\n' with open(mgt_file, 'w') as f: f.writelines(lines) # 批量运行不同 CN2 情景 for cn2 in [70, 75, 80, 85]: modify_cn2('project', cn2) os.system('swat.exe') # 调用 SWAT 可执行文件 shutil.copy('output.rch', f'output_cn2_{cn2}.rch')

逻辑说明：这段代码演示了批量修改参数并运行的基本流程。参数说明：CN2 的典型范围是 35-98，具体取值取决于土壤水文组和土地利用。注意每次运行前要备份原始文件，避免参数被覆盖。实际项目中，更规范的做法是用 SWAT-CUP 或 PySWAT 做参数采样和批量运行。

6.2 不确定性分析：GLUE 方法的简化实现

GLUE（广义似然不确定性估计）是水文模型常用的不确定性分析方法，核心思想是给参数随机采样，跑很多次，保留表现好的参数组，看预测区间。

# GLUE 简化实现：随机采样 + 保留 NSE > 0.5 的参数组 import numpy as np np.random.seed(42) n_samples = 1000 cn2_samples = np.random.uniform(35, 98, n_samples) alpha_bf_samples = np.random.uniform(0, 1, n_samples) results = [] for cn2, alpha_bf in zip(cn2_samples, alpha_bf_samples): # 这里省略实际运行 SWAT 的步骤 # 假设 run_swat 返回 NSE nse = run_swat(cn2, alpha_bf) results.append((cn2, alpha_bf, nse)) # 保留 NSE > 0.5 的参数组 good = [(c, a) for c, a, n in results if n > 0.5] print(f'保留 {len(good)} 组参数') # 用这些参数组的模拟结果计算 95% 置信区间

逻辑说明：GLUE 的关键是“行为阈值”，NSE > 0.5 是常用阈值。参数说明：采样次数一般 1000-10000 次，取决于计算资源。保留的参数组越多，不确定性区间越宽。实际项目中，GLUE 通常和 SWAT-CUP 结合使用，后者内置了 SUFI-2 算法，效率更高。

6.3 结果出图的几个细节

出图时注意坐标系和比例尺。子流域矢量要和 DEM 坐标系一致，否则叠不上。分级设色时，用自然断点法（Jenks）比等间距法更能反映空间差异。如果要做多情景对比，用统一图例，否则读者没法直观比较。

我自己做面源污染项目这些年，最大的教训是：不要迷信模型输出。SWAT 是个工具，它的结果取决于输入数据和参数设置，而面源污染的监测数据往往又少又粗。每次出图前，我都会问自己一句：这个负荷值合理吗？和文献比、和实地观察比、和上下游关系比，有没有明显矛盾？如果答案不确定，就回去检查数据。希望帮到你。

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