随着我国现代化进程的不断加速, 非点源污染与点源氮污染问题愈发严重。我国河流面临富营养化加剧、蓝藻水华灾害频繁、水质恶化和水生生态系统破坏等一系列生态问题^[1-2]。作为陆地与海洋之间物质迁移与能量交换的重要通道, 河流不仅承载了密集人类活动带来的环境压力, 也将这些压力逐步传递至海洋系统^[3-4]。因此, 建立河流入海氮污染模型, 对保障近海生态环境安全、推动流域与海洋环境协同治理具有重要意义^[5-6]。根据《2023中国海洋生态环境状况公报》显示, 中国沿海氮污染问题主要集中在途经人口密集或农业发达区域的河口地带^[7]。因此, 构建入海总氮通量预测模型, 识别影响入海总氮通量变化的关键驱动因素, 对于陆海污染协同控制和近海环境管理具有深远的意义。

目前用于入海水质预测的模型有两类, 一类是传统的物理机制模型, 如DELFT-3D^[8]、SWMM-EFDC^[9]、MIKE^[10]、SWAT^[11]模型等, 这类模型需要大量输入数据和参数设置, 并依赖长时间的模型预热过程, 难以捕捉水质影响因素间复杂的非线性关系^[12]。另一类是数据驱动的水质模型, 尤其是深度学习方法, 可以在不依赖物理过程知识的前提下, 从数据中自动学习输入与输出数据之间的非线性特征, 进而以简单的建模流程准确地实现对影响因素的分析^[13-14]。卷积神经网络(Convolutional Neural Networks, CNN)作为最具代表性的数据驱动模型之一, 在处理图像和其他多维数据方面表现出卓越的能力^[15]。CNN模型以人类视觉神经系统为设计灵感, 广泛应用于图像分类、目标识别、医学影像分析和物理过程建模等多个领域^[16-17]。在遥感、水资源等与水环境密切相关的领域, 深度学习技术也逐渐成为主流方法^[18-19]。尽管在2018年前水文学领域并未广泛使用深度学习模型^[20], 近年来, 随着模型结构与训练方法的不断优化, CNN模型在水文领域的应用逐渐增多^{[7, 21-22]}。例如, Pyo等^[23]基于CNN模型学习EFDC模型生成的水质图像, 成功预测了河流中的蓝藻浓度; Quang等^[24]利用韩国汉江10 a的水质监测数据, 基于不同深度学习架构对蓝藻浓度进行建模, 实现了高精度预测。这些研究验证了CNN模型在水文水质建模中的潜力和有效性, 它在处理图像等多维数据和多维特征方面的强大能力, 为精确预测入海水质提供了一种有效的解决方案。然而, 现有的水质预测模型仍存在一些显著不足: 多数模型依赖于监测点获取的低维水质数据(如氨氮、总氮、溶解氧等), 未直接考虑影响流域水质的多种潜在因素; 且河流设置的监测点有限, 导致预测模型的空间泛化能力较弱^[25]。此外, 现有模型大多基于时间序列或单点空间观测, 未能有效融入流域尺度上的空间结构及多种环境变量的交互作用^{[23-24, 26]}, 这限制了模型捕捉总氮通量驱动机制的能力与精度。

此外, 深度学习模型普遍存在“黑箱”问题, 即模型的内部机制难以解释, 影响了其科学解释性与结果可用性。近年来, 从博弈论中Shapley值扩展而来的SHAP(SHapley Additive exPlanations, SHAP)算法为深度学习提供了可解释工具。该方法已被广泛应用于机器学习与深度学习模型的特征重要性分析中。例如, Wang等^[27]利用SHAP识别了小清河流域产流过程的关键驱动因子; Yang等^[28]对滇池水质参数进行解释性建模; Rahmani等^[29]则基于SHAP对CNN模型进行可解释性分析, 从而识别影响地面沉降的关键因素。这些研究表明, SHAP方法可以有效提升CNN模型在环境建模中的透明度与可解释性, 增强模型的实际应用价值。

综上所述, 本文聚焦于小清河流域的入海总氮通量问题, 基于流域的气象、土地利用、土壤类型、点源与非点源氮排放等多源数据, 构建了一个融合丰富空间结构与多维特征信息的三维CNN输入形式, 并设计了相应的CNN预测模型用于模拟总氮通量的时空分布特征。同时, 结合SHAP方法对模型结果进行可解释性分析, 从而识别出入海总氮通量变化的关键影响因素及其空间异质性。本研究旨在为流域-海洋一体化氮污染管控与近海水质精准治理提供科学支撑。

1 材料与方法 1.1 研究区域与数据来源

小清河位于山东省内, 是黄河流域渤海水系的重要支流。其源头位于省会济南市的泉群区域, 依次流经济南市、淄博市、滨州市、东营市、潍坊市, 最终于潍坊的寿光市羊口镇进入莱州湾。小清河主干流全长237 km, 流域面积约10 336 km², 约占全省总面积的6.7%。小清河流域季节变化显著, 气候类型属于温带季风气候。作为莱州湾主要入海河流之一, 小清河在区域氮输入中占据重要地位, 其氮的输送贡献量达到30%, 仅次于黄河^[30]。

本文采用的DEM(Digital Elevation Model)数据来源于SRTM数据集(http://dwtkns.com/srtm30m/), 空间分辨率为30 m, 土地利用的数据来源于全球土地利用数据库(http://www.globallandcover.com/GLC30Download/index.aspx), 土壤数据来源于联合国粮农组织的全球土壤数据库(http://www.fao.org/soils-portal/data-hub/soil-maps-and-databases/harmonized-world-soil-database-v12/en/), 气象数据来源于中国气象同化驱动数据集(http://www.cmads.org/), 选取小清河流域内2009—2024年的日降水(P)、日最高气温(T_max)、日最低气温(T_min)、日平均相对湿度(R)、日平均风速(W)、日太阳辐射(S)这6种日尺度气象要素信息。点源和非点源的氮排放数据来自于济南市、淄博市、潍坊市统计年鉴(2009—2024年)以及全国第二次污染源普查, 其中非点源的排放主要考虑了3种指标: 日平均农村生活污水氮排放(DW)、日平均畜禽养殖氮排放(LW)、日平均农业施肥量(AF)。

鉴于水质数据的稀缺性, 本研究采用SWAT (Soil and Water Assessment Tool)模型对小清河流域的水质数据进行了扩展与补充^[31]。在SWAT模型构建过程中, 集成了小清河的逐日水文监测数据、逐月水质监测数据以及其他必要的地理与环境数据, 将整个流域划分为80个子流域和471个水文响应单元。随后, 利用SWAT的自带插件SWAT-CUP对模型进行了率定与验证: 率定期为2009—2016年, 验证期为2017—2024年。模型的验证基于位于小清河入海口处羊口监测站的月度总氮通量观测数据, 相关的率定与验证结果详见表 1, 总氮实测值与SWAT模拟结果的对比如图 1所示。从图 1可以看出, SWAT模型模拟结果与实测值整体拟合较好, 能够较为准确地反映出总氮通量的季节性变化特征。在完成对月尺度数据的有效率定和验证的基础上, 用构建的SWAT模型对2009—2024年小清河入海总氮通量进行了日尺度的模拟计算。

1.2 数据预处理

构建CNN模型所采用的输入数据主要包括气象数据(6个指标)、地形数据(8个指标)、土地利用(6个指标)、土壤类型(8个指标)、污染物排放数据(包括点源排放1个指标和面源排放3个指标), 模型的输出为小清河入海总氮通量。其中, 地形、土地利用和土壤类型被视为静态空间属性, 不随时间变化; 气象数据以日尺度时间序列的形式输入模型; 面源排放数据为年尺度时间序列, 包含生活和工业污水、畜禽养殖和农业施肥3类污染指标; 点源排放数据采用日平均排放量进行输入, 反映工业及城市污水处理厂的稳定排放过程。

由于输入特征之间存在不同的量纲和尺度差异, 可能会导致模型在权重更新过程中出现收敛困难或训练效率低下的问题。为提高模型训练的稳定性与收敛速度, 本文对所有输入数据进行了归一化处理, 采用Z-Score标准化方法, 将各特征值转换为具有相同均值和标准差的分布形式:

$ X^*=\left(X-\mu_X\right) / \sigma_X, $

(1)

式中: X^*为归一化后的数据; X为原始数据; μ_X为均值; σ_X为标准差。

在模型构建过程中, 原始数据被划分为训练集与测试集。训练集用于模型参数的学习, 测试集用于评估模型的泛化能力与性能表现。若训练集样本量不足, 可能导致模型对噪声过度拟合; 而测试集比例过大则会降低有效训练样本的利用率。结合现有研究经验^[25], 本研究将全部数据按8∶2的比例随机划分为训练集与测试集, 以在模型训练效果与数据利用效率之间取得平衡。

大型流域的河网在经过数字化之后, 其结构与马尔科夫链所描述的状态空间转换过程具有高度相似性。河流的流动过程符合马尔科夫性质, 即系统的下一状态仅依赖当前状态, 而与其历史状态无关。这种特性使得可以将流域河网的结构抽象为一个状态转移过程, 从而实现对流域空间结构的数字化建模。

基于上述思路, 本文将河流简化为1个N×N关系矩阵P。若子流域j是子流域i的最近下游空间单元, 则有H(i, j)=1; 否则, H(i, j)设置为0。如图 2所示, 矩阵P是基于马尔科夫链转移概率矩阵方法构建的河网关系矩阵, 可有效表达由多个子流域组成的上下游水文连接关系。由于流向关系是确定性的, 因此矩阵中的所有转移概率均为0或1。在此基础上, 将流域内的地形、气象、土地利用、土壤类型和氮排放等多源异构数据按照河网关系矩阵中河道的空间位置进行组织与填充, 每类数据对应一个通道, 多通道数据叠加后构成了用于模型输入的流域特征图像, 既保留了流域空间结构的水文连通性, 也实现了多源信息的融合表达。

1.3 CNN模型

卷积神经网络(CNN)是一种深度学习框架, 它主要由卷积层、池化层和全连接层组成(图 3)。卷积层通过卷积操作提取输入的特征, 使用多个卷积核来捕捉输入的局部依赖性和空间层次结构, 每个卷积核负责检测输入中的特定类型特征。池化层(或下采样层)则用于降低数据的空间维度和复杂性, 增强模型的泛化能力。通过堆叠多个卷积和池化层, CNN能够构建复杂的数据表征, 进而通过全连接层对这些特征进行分类或回归分析。

卷积层后的特征图大小计算公式为

$ N=(W-F+2 P) / S+1, $

(2)

式中: N为卷积后产生的特征图大小; W为输入矩阵大小; F为卷积核大小; P为填充大小; S为步长。

池化层后特征图大小计算公式为

$ N=[W+2 P-D(F-1)-1] / S+1, $

(3)

式中: D为通道数。

1.4 SHAP算法

CNN模型是一个具有大量参数的黑盒模型, 算法的透明性和可信度通常要借助于事后可解释方法进行分析。SHAP算法是一种以计算特征对模型输出的边际贡献为核心思想的事后解释的方法, 它利用博弈论中的Shapley值概念, 通过计算每个特征对模型输出结果的贡献程度, 量化每个特征在模型预测中的贡献。

对于模型中特征i的Shapley值的计算公式如下:

$ {\rm{ Shapley }}(i) = \sum\limits_{S \subseteq N\backslash i} {\frac{{S!(N - S - 1)!}}{{N!}}} [f(S\bigcup i ) - f(S)] $

(4)

式中: N表示所有特征的集合; S是不包含特征i的任何一个特征子集; f(S)是在S中具有特征的模型预测, f(S∪i)是具有S中的特征i的模型预测。SHAP值是Shapley值在机器学习模型解释中的实际应用形式, 借助高效算法实现对特征贡献的近似计算。

1.5 模型评价方法

模型预测结果采用纳什效率系数(NSE)和相关系数(R)来评价。这些指标的计算公式如下:

纳什效率系数(NSE):

$ {\rm{NSE}} = 1 - \frac{{\sum\limits_{i = 1}^n {{{\left( {{O_i} - {{\rm{P}}_i}} \right)}^2}} }}{{\sum\limits_{i = 1}^n {{{\left( {{O_i} - \bar O} \right)}^2}} }}. $

(5)

相关系数(R):

$ R = \sum\limits_{i = 1}^n {\frac{{\left( {{O_i} - \bar O} \right)\left( {{P_i} - \bar P} \right)}}{{\sqrt {\sum\limits_{i = 1}^n {{{\left( {{O_i} - \bar O} \right)}^2}} \sum\limits_{i = 1}^n {{{\left( {{P_i} - \bar P} \right)}^2}} } }}} , $

(6)

式中: O_i为实际值; P_i为模型预测值; O和P分别代表实际值和预测值的平均值。由于深度学习模型的参数训练具有一定的随机性, 因此选择分别训练3次后的评价指标的平均值作为最终的模型结果。

1.6 模型结构与超参数设置

本研究采用马尔科夫链转移概率矩阵的方法, 对各个指标的输入数据进行重新整合, 形成了一个尺寸为80×80×32的三维输入数据, 以此作为CNN模型的输入。作为对比, 我们还将输入数据按子流域编号整合为80×32的二维输入方式。在二维输入方式中, 每一行代表一个子流域, 每一列对应一个指标。这种方式不考虑子流域之间的相互关系, 因此每个子流域的指标数据是独立的。模型的主要结构流程如图 4所示。

硬件环境配置为Intel i7-12700H处理器, 配备NVIDIA GeForce RTX3060显卡(6GB显存)和16GB内存。模型使用Python3.9进行编程实现, 采用TensorFlow作为深度学习框架。超参数的优化采用贝叶斯搜索方法, 经过优化后得到的学习率为0.0026, 正则化丢弃率为0.65, 优化算法选用Adamax, 激活函数为ReLU。

2 结果与分析 2.1 模型效果评估

基于两种不同的输入结构, 分别构建了入海总氮通量预测的卷积神经网络(CNN)模型, 评价结果如表 2所示。通过对入海总氮通量的预测, CNN模型在处理大量数据时展现出了良好的拟合效果。其中, 基于马尔科夫链转移概率矩阵的三维输入方法在拟合效果上明显优于二维输入方法。在训练集上, 三维输入方法的纳什效率系数(NSE)和相关系数(R)分别达到了0.98和0.99, 均高于二维输入方法的0.89(NSE)和0.92(R), 这表明三维输入方法在捕捉数据中的空间结构和多层次特征方面表现更佳。具体来看, 三维输入方法相对于普通的二维输入方法, 在训练集上, NSE提升了约10.11%, R提升了约7.61%; 在测试集上, NSE提升了约5.68%, R提升了约4.35%。

与其他相关研究相比, 例如徐洁等^[32]对多个河流断面进行的总氮通量预测研究, 以及武暕等^[33]对辽宁省入海河流总氮浓度的预测研究, 本研究通过采用CNN模型对入海总氮通量的预测, 获得了更高的结果精度。这些提升反映出三维输入方法能够有效利用更丰富的信息进行特征提取, 进而提升了预测性能和模型的泛化能力。

此外, 使用三维矩阵表示数据使得CNN模型能够更好地处理复杂的数据模式, 从而显著提高了模型的效果和精度。图 5给出了小清河入海总氮通量的模拟值与实测值对比。从图中可以看出, 模型预测结果与实测数据整体趋势高度一致, 尤其在高值和突变点处仍能保持较好的拟合度, 进一步验证了CNN模型在入海总氮通量预测中的有效性。这些结果表明, 基于马尔科夫链转移概率矩阵的三维输入结构充分发挥了CNN模型处理多维数据的优势^[34]。

2.2 驱动因素重要性及影响分析

在本研究中, 采用SHAP的DeepExplainer方法来计算每个变量和每个空间位置的SHAP值, 量化不同变量及空间位置对模型预测结果的影响。SHAP值的绝对值越大, 表示该变量对模型输出的贡献越大。

图 6展示了SHAP值不为0的输入变量的重要性。从图中可以看出, 降水(P)是影响入海总氮通量的最主要因素, 其次是风速(W)和最低温度(T_min), 其中降水的特征重要性显著高于其他变量。在非点源污染中, 畜禽养殖氮排放(LW)是最具代表性的影响因子, 其次是农业施肥量(AF)。而地形、土地利用类型、土壤性质等静态变量由于不随时间变化, 其SHAP值均为0。这是因为CNN模型主要依赖输入变量的时序特征进行学习, 而这些静态变量仅在空间维度存在差异, 缺乏动态变化信息, 难以为模型提供有效的时变特征, 从而对预测结果影响有限。气象因素在模型中表现出更高的重要性, 主要得益于其强时效性和高时间分辨率。气象数据具有日尺度的连续变化, 能够更好地捕捉入海总氮通量的短期动态。而非点源污染数据多以年为单位记录, 序列波动较小, 短期变化不明显, 这限制了其在模型中对预测的即时响应能力。

为了进一步验证SHAP分析的结果, 本文计算了各输入变量与入海总氮通量之间的Pearson相关系数(图 7), 以刻画变量之间的线性关系。结果显示, 降水(P)与总氮通量的线性相关性最高, 其次为相对湿度(R)、最低温度(T_min)和风速(W)。在非点源排放中, 农业施肥量(AF)与入海总氮通量的相关性高于畜禽养殖(LW), 而农村生活污水(DW)的相关性最低。图 7还补充了变量之间的相关性结构。我们发现, 最高温度(T_max)和最低温度(T_min)之间存在较强相关性, 可能会导致模型中出现多重共线性问题, 从而影响特征权重的分配。

图 8展示了SHAP值随降水和最高温度变化的交互变化图。横坐标和色带分别代表经过归一化后的降水值和最高温度, 纵坐标表示SHAP值。其中图 8(a)显示, 降水量越大, 其SHAP值越高, 表明其对总氮通量的贡献持续增强, 尤其在极端降水事件下影响更为显著, 这也解释了为何图 6中的降水具有最高的平均SHAP值。图 8(b)则表明, 最高温度对SHAP值的影响没有明显的单调变化趋势, 说明其作用可能更为复杂, 受其他变量调控。

结合SHAP和Pearson相关系数的分析结果, 降水被一致认为是最关键的驱动因素, 这也与大量已有研究结果一致^[35-38]。降水会直接增加地表径流, 促进氮素从农业面源、污染源区域向水体迁移^[39-40]。在强降水条件下, 溶解态氮和颗粒态氮更易被冲刷进河道, 从而提高了入海总氮通量。温度同样是重要影响因子。SHAP分析显示, 最低温度对模型的影响更为显著, 可能是由于生物化学反应速率对最低温度更敏感^{[38, 41]}, 在低温条件下, 氮素的矿化、硝化和反硝化过程可能发生显著变化, 进而影响氮的迁移与转化路径。而Pearson结果中, 最高温度和最低温度的相关性较高(图 7), 这可能导致CNN模型在评估特征重要性时偏向选择具有更非线性特征的最低温度, 从而避免多重共线性带来的不确定性^[42]。在非点源污染中, 农村生活污水的入海总氮通量SHAP权重和Pearson相关性都是最低的, 这主要由于其总排放量相对较小, 对总氮通量的实际贡献有限^[43]。

综上所述, CNN模型在预测入海总氮通量时, 对时间动态性强的变量表现出更高敏感度, 尤其是降水、温度等气象因子和非点源排放的动态数据。而模型在处理静态变量方面能力有限, 这提示我们在未来的建模中, 可以考虑引入动态地表信息(如遥感数据或时序土地利用变化)以提升模型的综合表达能力和预测准确性。

2.3 空间特征影响的SHAP分析

基于采用三维输入方式的CNN模型, 进一步计算了小清河各子流域的空间平均SHAP值, 结果如图 9所示。分析结果表明, 各影响因素在空间尺度上对入海总氮通量的贡献存在显著差异。对入海总氮通量影响程度较大的区域主要集中在中部和西部, 包括淄博市张店区附近的55号与57号子流域, 以及济南市历城区附近的53号与54号子流域。张店区和历城区均为区域内人口集中、城市化程度较高的城区, 城市生活污水和工业废水排放量大, 为氮素输入提供了持续来源; 并且这些区域周边农业种植密集, 农田施肥量较高, 降雨后农业面源污染更易进入水体。

每个子流域对入海总氮通量的影响是多因素共同控制的, 受到多方面的影响。图 10展示了小清河入海总氮通量影响最大的2个气象因素的空间SHAP分布。其中降水的贡献主要体现在流域中部和西部, 风速的贡献主要体现在西部。

图 11展示了小清河3类非点源污染(生活污水、畜禽养殖和农业施肥)的平均日负荷空间分布图, 以及相应的空间SHAP值。从对比结果可以看出, 实际排放负荷的空间分布与其对入海总氮通量影响的SHAP值之间存在显著异质性, 说明非点源排放对总氮通量的影响不仅仅取决于排放强度, 还受到多种环境和过程因素的共同作用。

首先, 非点源污染本身具有分布广泛、排放随机和受环境制约强等特征, 其空间分布往往受到土地利用类型、气候条件、地形地貌等因素的综合影响^[44-45]。其次, 气候变化也是影响氮素迁移与转化的重要驱动因素, 温度升高和降水模式的改变可能通过调节土壤微生物活性而影响氮循环过程^[37]。例如, 极端降雨事件可导致地表径流迅速增强, 进而引发氮素在短时间内大量流失, 显著抬升入海总氮通量^[46]。这种现象也在图 8中得到佐证: 图 8(a)展示了极端降雨下SHAP值显著升高的趋势, 图 8(b)则揭示了降水与温度之间的交互作用。

对比排放负荷与SHAP分布可以发现, 在部分区域, 排放强度较高的子流域也对应了较高的SHAP值, 这说明该区域排放对总氮通量具有显著贡献。例如, 35号子流域位于淄博市中心城区下游, 流经周村区和张店区等人口密集区域, 生活污水排放量大, 氮输入负荷显著; 13号子流域覆盖淄博市高青县, 该县为典型的畜禽养殖大县, 畜禽养殖活动频繁, 氮排放强度高; 49号子流域地处临淄区东部与寿光市西部交界区域, 该地以蔬菜大棚农业为主, 农业施肥强度大, 是典型的农业面源污染高发区。这些子流域因其排放源强大, 导致对应的SHAP值明显升高, 说明其对入海总氮通量的影响不仅在排放量层面上显著, 也在模型中体现出较强的驱动效应。这一结论与李运东^[47]的研究结果相一致, 进一步验证了本研究模型对关键影响因子的准确识别能力。

然而, 并非所有排放量大的区域都对入海总氮通量有显著贡献。如历城区西南部的73、74、75号等上游子流域虽然具有较高的污染源强, 但由于其距离入海口较远, 氮素在向下游输运过程中受到沉降、吸附、植物吸收、反硝化等多种物理和生物过程的影响, 浓度逐步衰减, 从而削弱了其最终对入海总氮通量的影响^[48]。

相反, 部分距离入海较远的区域如济南市中心的54、59和61号子流域, 其空间SHAP值却较高。这种现象主要与当地气象条件密切相关。如图 7与图 8所示, 强降水事件在这些子流域频发, 加剧了地表径流强度, 从而加速了氮素从源区向河道的迁移过程, 即便源强一般, 其输送效率却显著增强^[49]。这表明, 水动力过程与气象条件的空间异质性会改变氮素排放与其实际入海影响之间的传统认知关系。

综上所述, CNN模型借助SHAP方法所揭示的空间影响结构, 不仅反映了排放源强的直接贡献, 还捕捉到了污染物输运路径、气象条件、空间位置等多重因素对入海总氮通量的复合作用机制。这也体现出深度学习模型在非线性特征提取和复杂过程识别方面的独特优势, 能够更全面地揭示非点源污染在流域尺度下的真实影响过程^[50], 为流域污染治理和区域优化配置提供了有力支持。

3 结论

(1) 基于马尔科夫链中生成转移概率矩阵方法创建的三维输入结构, 更适合CNN模型处理复杂的数据模式, 有效提高了CNN模型利用多种入海总氮影响因素对入海总氮通量预测的效果和精度。

(2) 在众多入海总氮通量影响因素中, 降水是影响入海总氮通量变化的最重要的特征。

(3) 对小清河入海总氮通量影响程度较大的区域主要集中在中部淄博市区和西部济南市区。每个子流域对入海总氮通量的影响是多因素共同控制的, 受到污染物排放、气象因素等多方面的影响, 并且多因素的叠加影响会显著改变原有区域间的排放特征。