基于贝叶斯累加回归树模型的溶解氧影响因素分析及预测研究——以长江中游为例

Analysis and Prediction of Factors Affecting Dissolved Oxygen Based on a Bayesian Additive Regression Tree Model—A Case Study of the Middle Reaches of the Yangtze River

摘要: 本文以长江中游三个河流断面每隔4小时的水质数据和气象指标为研究对象进行溶解氧的影响因素分析和预测研究。首先，通过模拟数据对BART、SVR、RF和XG模型进行实验对比，分析了BART模型的稳健性优势。其次，运用BART模型变量重要性(变量包含比例)来筛选重要变量，筛选出影响长江中游水质溶解氧浓度的关键变量为pH、水温、电导率、总磷、高锰酸盐指数、露点温度、气温。最后，选择筛选出的最优输入变量组合来构建BART、SVR、RF和XGB模型来对溶解氧进行预测研究。结果表明，BART模型的预测性能最佳，在100次预测统计中R²、RMSE、MAE、MAPE的中位数(M_e)分别为0.936、0.480、0.333、4.36%。

Abstract: This study focuses on the analysis and prediction of dissolved oxygen (DO) in the middle reaches of the Yangtze River, using water quality data and meteorological indicators measured at three river cross-sections every 4 hours. First, experiments were conducted to compare the performance of four models: BART, SVR, RF, and XGB, using simulated data. The robustness advantage of the BART model was highlighted. Next, the variable importance of the BART model was utilized to select key influencing factors, identifying the critical variables affecting DO concentration in the middle reaches of the Yangtze River. These key variables include pH, water temperature, electrical conductivity, total phosphorus, permanganate index, dew point temperature, and air temperature. Finally, the optimal combination of input variables was selected to build the BART, SVR, RF, and XGB models for DO prediction. The results show that the BART model performs best, with median values of R², RMSE, MAE, and MAPE of 0.936, 0.480, 0.333, and 4.36%, respectively, from 100 prediction trials.

1. 引言

1.1. 研究背景

水质监测对于社会经济发展有着至关重要的影响。随着社会及工业的迅猛进展和人口的不断增长，水资源污染日益严重，水质监测作为水污染控制工作中的基础性工作，其意义和作用也变得愈加重要。长江流域的中下游人口密度、经济发展程度在改革开放以来均位居我国前列，随之而来的环境压力也不容忽视[1]。伴随着长江流域中下游人口和用水量的激增，给水处理系统造成巨大压力，也意味着流域内大量污水排入长江[2]。随污水排放的耗氧污染物导致的水环境危机在流域经济快速发展的过程中已有显现，DO较低、CODMn和NH₃-N负荷高为水体恶化的标志[3] [4]。

1.2. 研究现状

在大数据与人工智能的时代背景下，传统的统计方法正逐渐与先进的计算技术融合，从而极大地扩展了数据分析的边界和深度。溶解氧(DO)是评价地表水水质质量的重要指标之一[5]，过低的溶解氧会降低地表水的自我净化能力，还会导致河流呈现还原环境，增大重金属等溶解释放强度[6]。杨明悦等[7] (2022)以深圳湾为例，运用Pearson相关性分析，变量重要性评分和随机森林方法构建了溶解氧实时测试模型，证明了用随机森林预测溶解氧具有可行性。黄燏等[8] (2023)基于长江流域21个水质检测断面2008~2018年的pH、溶解氧、高锰酸钾指数和氨氮数据，采用Mann-Kendall趋势检验、相关性分析和层次聚类分析，并使用绝对主成分回归方法解析污染物来源。

杨春艳等[9] (2022)分析2013~2020年泸沽湖溶解氧的变化规律，通过Spearman相关分析发现溶解氧含量月变化与水温、pH显著正相关。陈湛峰等[10] (2023)结合水质指标(水温、pH、浊度、电导率、溶解氧)、氨氮、高锰酸钾指数及总氮和气象指标包括温度、湿度、气压、风速、风向和降雨量对珠江口总磷含量进行预测，结果表明，气象指标对水质会有一定的影响。郝玉莹等[11] (2021)用随机森林对影响水质指标的特征进行筛选，然后用LSTM模型进行预测，结果表明用随机森林进行变量选择提高了预测的精度。Bleich等[12] (2014)对比了BART、LASSO、随机森林变量重要性选择的性能，发现BART具有优势。

在模型上，李玲玲等[13] (2020)用决策树对城市黑臭水体进行分析，Shi Z B等[14] (2014)利用一种结合小波分析的ARIMA模型，通过小波变换将水质监测数据分解为高频和低频部分，并分别应用于ARIMA模型中，以预测未来的溶解氧浓度。李晓瑛等[15] (2024)选取了长江口六个重点检测站点，使用改进的支持向量机回归、人工神经网络和随机森林三种模型预测溶解氧含量，随机森林的总体平均误差为0.19，效果最好。

2. 数据收集与分析

2.1. 数据收集

本文以长江中游三个重点断面二年的水质作为研究对象，收集的是2021年1月1日到2022年10月28日三个断面(湖南长沙新港、湖北武汉宗关和江西南昌滁槎)的水质数据，采集频率为4 h，长江中游三个监测断面位置如图1所示。本文以溶解氧为目标变量，其余11项指标作为预测变量，包括水温、pH、溶解氧、电导率、浊度、高锰酸盐指数、氨氮、总磷和总氮8项水质指标和平均气温、平均露点温度、降水量这3项气象指标。考虑到站点维护时存在异常值和缺失值，本研究将维护站点的数据删除后，再通过线性插补进行补全，得到共8679条数据。水质数据来源于国家水质自动监管平台，气象数据来源NOAA气象数据。

Figure 1. Distribution map of the middle reaches of the Yangtze River

图1. 长江中游断面分布图

2.2. 描述性分析

表1统计了水质和气象指标的特征统计，分别计算出溶解氧和预测变量的平均值、标准差、变异系数和中位数，来展示数据分布情况。其中溶解氧的平均值为8.22 mg/L，最大值达到16.87 mg/L。平均水温为20.68℃，最高能达到39.75℃；平均pH为7.42，水质大体上呈现中性。在所有变量中，氨氮的变异系数最大，浊度次之，说明这些变量的波动程度相对较大；而总磷的变异系数最小，氨氮次之，说明其波动程度较小。

Table 1. Statistics of water quality and meteorological characteristics

表1. 水质和气象特征统计

指标	符号	平均值	标准差	变异系数	最小值	中位数	最大值
溶解氧(mg/L)	y	8.22	1.90	23.12	2.06	8.06	16.87
水温(℃)	x1	20.68	7.56	36.53	6.53	20.80	39.75
pH	X2	7.42	0.52	7.00	5.26	7.34	9.83
电导率($\text{μS}/\text{cm}$ )	X3	250.53	75.56	30.16	0.01	235.90	460.66
浊度(NTU)	X4	43.93	49.01	111.57	2.63	29.61	498.47
高锰酸盐指数(mg/L)	X5	2.35	1.01	43.20	0.25	2.24	20.68
氨氮(mg/L)	X6	0.13	0.16	123.94	0.03	0.08	7.72
总磷(mg/L)	X7	0.07	0.04	66.07	0.01	0.06	0.59
总氮(mg/L)	X8	2.04	1.11	54.37	0.05	1.89	11.16
气温($˚\text{F}$ )	X9	66.66	16.55	24.82	26.00	67.60	95.00
露点温度(Pa)	x10	56.75	15.93	28.08	2.90	58.20	81.90
降雨量(mm)	X11	0.15	0.42	270.30	0.00	0.00	4.53

2.3. 相关性分析

Pearson相关系数衡量两个定量变量之间的线性相关程度，且要求数据服从正态分布。它的值范围从−1到+1，其中+1表示完全正线性相关，−1表示完全负线性相关，0表示无线性相关。公式为：

$\gamma =\frac{{\displaystyle \sum \left(X_i-\overline{X}\right)\left(Y_i-\overline{Y}\right)}}{\sqrt{{\left(X_i-\overline{X}\right)}^2{\displaystyle \sum {\left(Y_i-\overline{Y}\right)}^2}}}$ (1)

其中，$X_i$ 与$Y_i$ 是变量观测值，$\overline{X}$ 与$\overline{Y}$ 是变量的均值。

Figure 2. Pearson correlation coefficients among the variables

图2. 各变量间的Pearson相关系数

通过计算溶解氧与各预测变量的Pearson相关系数，得到的相关系数矩阵如图2所示。其中，水温、气温、露点温度与溶解氧呈强负相关性($\gamma >-0.7$ )，pH与溶解氧呈现中度正相关性($0.3<\gamma <0.4$ )，都具有统计学意义($\text{P}<0.05$ )。其它变量与溶解氧呈现弱相关性($-0.17<\gamma <0.27$ )。此外，图2结果还显示出了预测变量间存在中度以及强的多重共线性。例如，水温和气温、露点温度，pH和电导率，气温和露点温度都呈现强相关性($\gamma >-0.7$ )，总磷和总氮、高锰酸盐指数、浊度呈现中度相关性($0.3<\gamma <0.5$ )。

3. 研究方法

3.1. 模型原理

3.1.1. 随机森林

随机森林是一种基于决策树的机器学习方法，它是利用bootstrap重抽样技术，在原始样本中抽取K个样本，分别训练K颗决策树组成随机森林[16]，用于回归问题，通过计算平均值得出预测值[17]，其原理如图3 (a)所示。

3.1.2. 支持向量回归

支持向量回归(SVR)是一个用于处理回归数据的机器学习方法，它源自于一个用于处理分类数据的模型支持向量机(SVR) [18]。SVR的目标是最小化样本点与最优超平面的总偏差，如图3(b)所示。

3.1.3. XGBoost

Figure 3. Random Forest principle (a), Support vector regression principle (b) and XGBoost principle (c)

图3. 随机森林原理 (a)、支持向量回归原理 (b)及XGBoost原理 (c)示意图

极限提升树(extreme gradient boosting, XGBoost)是多棵CART树组成的集成学习算法，是在梯度提升树(gradient boosting decision tree, GBDT)的基础上改进而来的[19]。XGBoost的核心是对目标函数进行二阶泰勒展开，利用不断分裂的新树拟合上次预测的残差以保证模型的准确性，其原理如图3(c)所示。

3.1.4. BART模型

贝叶斯累加回归树(BART)模型分为3个主要部分：树和模型、正则化先验和Backfitting MCMC算法[20]。

1) 树和模型

对于样本$\left(x,y\right)$ 考虑BART模型

$Y=f\left(x\right)+\epsilon ={\displaystyle \sum _{i=1}^{m}g\left(x;T_j,M_j\right)}+\epsilon ,\epsilon ~N\left(0,{\sigma }^2\right)$ (3)

其中，$\left(T_j,M_j\right)$ 表示第j棵回归树模型，$T_j$ 表示第j棵回归树，$M_j=\left\{{\mu }_{1j},{\mu }_{2j},\cdot \cdot \cdot {\mu }_{bj}\right\}$ 是$T_j$ 的$b_j$ 个终节点对应的参数，$g\left(x;T_j,M_j\right)$ 是将${\mu }_{ij}\in M_j$ 分配给x的函数，$\epsilon ~N\left(0,{\sigma }^2\right)$ 为随机误差项，$\sigma$ 表示BART模型的标准差。

2) 正则化先验

对参数做正则化先验，一方面是为了防止某棵树生长过大而占据预测的主导地位，从而失去建模的意义，另一方面是为了便于计算。由式(2)可见，给定树的数量m，树和模型Y由$\left(T_1,M_1\right),\cdot \cdot \cdot ,\left(T_m,M_m\right),\sigma$ 决定，即模型的参数有$T_1,M_1,\cdot \cdot \cdot ,T_m,M_m,\sigma$ 。

假设每棵树之间是相互独立的，且独立于误差项。同时假设，在给定树的具体结构条件下，树的叶子结点的预测值之间是相互独立的。因此：

$\begin{array}{l}P\left(\left(T_1,M_1\right),\cdot \cdot \cdot ,\left(T_m,M_m\right),\sigma \right)=P\left(\left(T_1,M_1\right),\cdot \cdot \cdot ,\left(T_m,M_m\right)\right)P\left(\sigma \right)\\ =\left\{{\prod }_{j}P\left(T_1,M_1\right)\right\}P\left(\sigma \right)=\left\{{\prod }_{j}P\left(M_j|T_j\right)P\left(T_j\right)\right\}P\left(\sigma \right)\\ =\left\{{\prod }_j\left({\prod }_{j}P\left({\mu }_{ij}|T_j\right)\right)P\left(T_j\right)\right\}P\left(\sigma \right)\end{array}$ (3)

通过独立假设，可将问题简化为只需分别设计$P\left(T_j\right)$ ，$P\left({\mu }_{ij}|T_j\right)$ 和$P\left(\sigma \right)$ 的先验。

3) Backfitting MCMC算法

BART模型参数训练使用的是贝叶斯Backfitting MCMC算法，该算法的思想为：通过贝叶斯推断确定后验分布$P\left(\left(T_1,M_1\right),\cdot \cdot \cdot ,\left(T_m,M_m\right),\sigma |y\right)$ ，基于确定的后验分布进行吉布斯抽样并建立统计推断。

3.2. 变量筛选原理

BART模型的主要输出是对响应变量y的一系列预测值。尽管这些预测值可以用于描述模型的整体拟合情况，但它们并不能直接用于评估每个预测变量的相对重要性。为了解决这一问题，Chipman等[20]在其研究中开始探索每个预测变量的“变量包含比例”以作为BART模型的变量重要性度量，旨在衡量每个变量对整体模型的贡献程度，有助于选择最具影响力的预测变量。本文运用BART模型中通过评估在构建树过程中每个特征被选择的频率，即“变量包含比例”。其中“变量包含比例”表示在树和模型的后验样本中，变量被选为分裂规则的次数占所有分裂规则中的比例[20]。

$v_i=\frac{1}{k}{\displaystyle \sum _{k=1}^k{z}_{ik}}$ (4)

其中，$z_{ik}$ 为x的第i个分量的分裂比例。

定义$RMS{E}_{change}$ 为相对变化率，计算公式如下所示：

$RMS{E}_{change}=\left(\frac{RMS{E}_{n-1}-RMS{E}_n}{RMS{E}_{n-1}}\right)\times 100\%$ (5)

其中$RMS{E}_n$ 代表第n个变量组合构建模型的$RMSE$ ,设定$RMS{E}_0=\infty$ 。

变量筛选的算法步骤：

(1) 首先通过BART对所有预测变量计算变量包含比例，并以降序方式排序；

(2) 选择序列中的第一个预测变量对响应变量构建BART模型，进行k折交叉验证计算模型性能，得到$RMS{E}_1$ ；

(3) 若$RMS{E}_{change}=\left(\frac{RMS{E}_{n-1}-RMS{E}_n}{RMS{E}_{n-1}}\right)\times 100\%$ > 阈值(阈值默认设定为3%)，则认为模型的性能

有显著减小，将第n个变量添加入最优变量组合中；否则执行步骤(5)。

(4) 依次添加序列中排名后一位的预测变量得到新的变量组合，对响应变量y重新构建BART模型，进行k折交叉验证计算模型性能得到$RMS{E}_n$ (n为预测变量数目)，返回执行步骤(3);

(5) 输出模型的最优变量组合。

3.3. 评价指标

在评估模型性能上，使用到的指标是，决定系数R²，均方根误差$RMSE$ ，平均绝对误差$MAE$ 和平均绝对百分比误差$MAPE$ 。

$R^2=1-\frac{{\displaystyle {\sum }_i{\left({\widehat{y}}_i-y_i\right)}^2}}{{\displaystyle {\sum }_i{\left(\overline{y}-y_i\right)}^2}}$ (6)

$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left({\widehat{y}}_i-y_i\right)}^2}}$ (7)

$MAE=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|{\widehat{y}}_i-y_i\right|}$ (8)

$MAPE=\frac{100\%}{n}{\displaystyle \sum _{i=1}^n\left|\frac{{\widehat{y}}_i-y_i}{y_i}\right|}$ (9)

其中n样本数，$y_i$ 为因变量真实值，${\widehat{y}}_i$ 为因变量预测值，$\overline{y}$ 为因变量真实值的平均值。

4. 数值模拟实验

4.1. 实验设计

本文采用模拟数据对BART模型的稳健回归性能做了探讨研究，并与常见三种机器学习方法作比较，包括支持向量回归(SVR)、随机森林(RF)和XGBoost (XGB)模型。本文模型实验采用Friedman [21]提出的包含一次项、二次项、交叉项和周期函数的经典回归模型：

$Y=10\sin \left(X_1{X}_2\right)+20{\left(X_3-0.5\right)}^2+10X_4+5X_5+\epsilon ,X_n~U\left(0,1\right),n=1,2,\cdot \cdot \cdot ,5$ (10)

其中，Y为目标变量，$X_n$ ，$\left(n=1,2,\cdot \cdot \cdot ,5\right)$ 是与Y具有函数关系的变量，$\epsilon$ 为随机误差项。

本文采用Daniel VandenHeuvel [22]提出的一种对数据进行随机攻击的污染模式：

$y_{i,a}=\left(1+\frac{s}{100}\right)y_i,s~N\left(\mu ,{\sigma }^2\right)$ (11)

随机攻击通过$1+s\%$ 因子随机地缩放数据，其中$y_i$ 为原始数据中受到污染的数据值，$y_{i,a}$ 为$y_i$ 受到污染后缩放的数据值。同时，设定P为污染率，以一定比例P对原始数据进行污染。

本文设定了不同μ的污染情形，μ符号正负代表缩放情形，μ为正时，代表数据被放大；μ为负时，代表数据被缩小。同时，考虑到不同攻击的强度，设计了$\left|\mu \right|=5,15,25,\sigma =\left|\mu \right|/10$ 的情形，分别代表小强度，中强度和高强度的随机攻击情形。综合各种情形，最终设定了$\mu =-25,-15,-5,5,15,25$ 这6种情形，每种情形均设定了污染率$\text{P}=10\%,20\%,30\%,40\%$ 来对数据进行污染。

本文模拟实验均采用随机产生数据的方式独立进行100次，每次模拟生成样本量n = 200 的独立同服从均匀分布的数据，随机划分训练集和测试集的比例均为50%，在训练集上进行随机攻击污染数据。训练集数据用于训练生成BART、RF、SVR和XGB模型，在未受污染的测试集数据上测试和比较模型的预测性能。

4.2. 实验结果分析

表2给出了$\mu =-5,-15,-25$ 在$\text{P}=10\%,20\%,30\%,40\%$ 的污染情形下BART、SVR、RF和XGB模型在100次试验中测试集上三个评价指标RMSE、MAE、MAPE的下四分位数(Q₁)、中位数(M_e)和上四分位数(Q₃)。

Table 2. Prediction evaluation results of each model under the conditions of $\mu =-5,-15,-25$ and $\text{P}=10\%,20\%,30\%,40\%$

表2. $\mu =-5,-15,-25$ ，$\text{P}=10\%,20\%,30\%,40\%$ 情形下各模型的预测评价结果

	P	Method	RMSE			MAE			MAPE
			Q1	Me	Q3	Q1	Me	Q3	Q1	Me	Q3
−5	10%	BART	1.633	1.780	1.858	1.287	1.390	1.475	10.67%	11.78%	13.10%
		SVR	1.860	2.030	2.163	1.446	1.545	1.628	11.97%	13.74%	15.84%
		RF	2.663	2.875	3.023	2.146	2.310	2.451	18.31%	22.33%	23.66%
		XGB	2.625	2.785	2.968	2.053	2.221	2.386	16.65%	17.97%	20.03%
	20%	BART	1.693	1.795	1.881	1.324	1.414	1.484	10.65%	11.49%	12.99%
		SVR	1.875	2.004	2.198	1.435	1.523	1.658	11.53%	12.96%	15.30%
		RF	2.767	2.901	3.089	2.233	2.359	2.510	19.27%	21.35%	23.57%
		XGB	2.670	2.857	3.026	2.155	2.277	2.429	16.31%	18.03%	20.09%
	30%	BART	1.638	1.781	1.897	1.305	1.412	1.518	10.70%	11.82%	12.93%
		SVR	1.803	2.007	2.164	1.406	1.514	1.652	11.86%	13.00%	14.80%
		RF	2.693	2.869	3.018	2.152	2.332	2.479	19.03%	20.80%	23.66%
		XGB	2.653	2.778	3.015	2.109	2.210	2.424	16.17%	17.83%	19.26%
	40%	BART	1.630	1.747	1.887	1.308	1.377	1.481	10.48%	11.68%	12.79%
		SVR	1.868	2.051	2.201	1.442	1.576	1.674	11.93%	13.23%	15.08%
		RF	2.742	2.881	3.031	2.210	2.340	2.456	19.40%	20.76%	23.52%
		XGB	2.699	2.909	3.075	2.186	2.305	2.453	16.60%	18.37%	19.50%
−15	10%	BART	1.726	1.810	1.964	1.363	1.475	1.542	11.20%	12.02%	13.32%

续表

		SVR	1.924	2.028	2.169	1.466	1.576	1.671	12.00%	13.13%	16.13%
		RF	2.727	2.897	3.044	2.207	2.331	2.468	18.92%	21.07%	24.97%
		XGB	2.756	2.931	3.072	2.179	2.304	2.436	17.01%	18.29%	20.22%
	20%	BART	1.756	1.880	2.035	1.386	1.498	1.600	11.06%	12.23%	13.66%
		SVR	1.985	2.137	2.351	1.532	1.643	1.774	12.37%	14.10%	16.64%
		RF	2.780	2.978	3.129	2.217	2.407	2.542	19.31%	21.26%	24.56%
		XGB	2.797	3.010	3.242	2.233	2.405	2.576	17.49%	18.98%	20.96%
	30%	BART	1.805	1.923	2.072	1.450	1.530	1.643	11.51%	12.45%	13.43%
		SVR	2.094	2.223	2.416	1.595	1.747	1.873	12.83%	14.50%	15.80%
		RF	2.765	2.993	3.206	2.225	2.419	2.592	19.05%	20.94%	23.28%
		XGB	2.889	3.079	3.284	2.277	2.452	2.623	17.47%	18.94%	20.28%
	40%	BART	1.919	2.086	2.191	1.516	1.631	1.767	11.89%	13.17%	14.24%
		SVR	2.173	2.367	2.546	1.711	1.833	1.981	12.99%	14.48%	16.78%
		RF	2.921	3.075	3.287	2.359	2.492	2.619	19.47%	21.44%	24.63%
		XGB	2.921	3.151	3.375	2.364	2.524	2.686	17.99%	19.54%	21.15%
−25	10%	BART	1.791	1.894	2.034	1.400	1.500	1.593	11.09%	12.24%	14.00%
		SVR	1.966	2.082	2.234	1.522	1.591	1.716	12.11%	13.43%	15.69%
		RF	2.747	2.899	3.123	2.178	2.367	2.504	18.03%	20.41%	23.43%
		XGB	2.872	3.024	3.204	2.254	2.425	2.524	17.11%	18.21%	20.03%
	20%	BART	1.940	2.069	2.220	1.548	1.628	1.773	12.28%	13.14%	14.49%
		SVR	2.170	2.321	2.530	1.658	1.772	1.920	13.16%	14.50%	16.73%
		RF	2.841	3.033	3.276	2.304	2.462	2.618	19.00%	20.97%	23.53%
		XGB	3.005	3.154	3.402	2.414	2.499	2.683	17.81%	19.06%	20.65%
	30%	BART	2.121	2.292	2.468	1.702	1.829	1.984	12.84%	13.81%	15.04%
		SVR	2.437	2.561	2.815	1.899	1.990	2.148	14.01%	15.97%	17.34%
		RF	3.053	3.197	3.412	2.438	2.608	2.754	19.73%	21.47%	23.34%
		XGB	3.255	3.437	3.690	2.589	2.729	2.941	19.28%	20.68%	21.72%
	40%	BART	2.337	2.463	2.641	1.847	1.984	2.124	13.89%	14.66%	16.29%
		SVR	2.641	2.888	3.080	2.097	2.239	2.428	15.05%	16.17%	18.11%
		RF	3.144	3.297	3.621	2.508	2.663	2.875	19.17%	20.71%	22.99%
		XGB	3.384	3.645	3.866	2.679	2.912	3.101	19.54%	20.70%	22.56%

由表2可以看出，BART模型的RMSE、MAE、MAPE的Q₁、M_e和Q₃在相同污染情形下都会远小于SVR、RF和XGB模型，在所有指标中均表现出最好的预测性能。具体而言，以$\mu =-25$ ，$\text{P}=40\%$ 这种大范围高强度攻击下为例，BART模型的MAPE中位数为14.66%，而SVR、RF、XGB模型的MAPE中位数分别为16.17%、20.71%、20.70%，分别比BART模型高出了10.3%、41.27%、41.20%。由BART得出MAPE的下四分位数(Q₁)、上四分位数(Q₃)是最小的，它们的差值也是最小的，即$Q_3-Q_1=2.4\%$ ，会低于其他三个模型(SVR、RF、XGB模型的$Q_3-Q_1$ 的值分别为3.06%、3.08%、3.02%)。同样的，对RMSE和MAE来进行对比来看，也是有类似的结果。

图4给出了在$\mu =-5,-15,-25,15,25$ ，$\text{P}=10\%,20\%,30\%,40\%$ 情形下得到四个模型的平均RMSE折线图，从总体上可以看出，随着污染比例P的提高和攻击强度μ的提高，各模型的RMSE均有所上升，这是由于更大的P和μ意味着引入了更多噪声。从各模型的对比结果来看，BART模型(标圆)的RMSE在所有不同的P和μ的设定下会低于其他模型。这表明，尽管异常数据对于所有模型都具有一定的影响，但BART模型在异常数据下的预测性能相对更为稳健和优越。综上所述，说明BART模型相对其他模型更优，具有更高的预测准确性和更小的预测误差。

Figure 4. Comparison of RMSE across models under the conditions $\mu =-5,-15,-25,5,15,25$ and $\text{P}=10\%,20\%,30\%,$ $40\%$

图4. $\mu =-5,-15,-25,5,15,25$ ，$\text{P}=10\%,20\%,30\%,40\%$ 情形下各模型的RMSE对比

5. 结果和分析

5.1. 预测因子的筛选

运用BART模型对长江中游水质数据计算各预测变量的重要性(变量包含比例)，独立重复进行100次，得到的各预测变量重要性得分箱线图如图5所示。其中，pH是预测溶解氧最重要的变量，随后依次是水温，电导率，总磷等；而变量重要性较小的是降雨量、浊度，这代表降雨量和浊度对预测水体中的溶解氧的影响较小。

将各变量根据变量重要性得分排序后，首先选取变量重要性得分最大的变量为第1个组合，组合1是pH (变量重要性得分最高的)；其次，再依据得分排序顺序逐次加入变量，从而得到一系列不同变量的组合，如图5所示。

本文采用5折交叉验证对各变量组合进行变量筛选，将数据集随机划分为5个子集，其中4个子集用于训练模型，剩余1个子集用于测试。该过程重复5次，每次使用不同的子集作为测试集，其余子集作为训练集。最终，计算在每个测试集上的评价指标RMSE、MAE、MAPE和R²，并对不同组合的模型结果进行对比分析。

Figure 5. Box plots of variable importance scores and combinations of different predictor variables

图5. 变量重要性得分箱线图和不同预测变量的组合

Table 3. Results of variable screening

表3. 变量筛选结果

组合	R^2	RMSE	MAE	MAPE	$RMS{E}_{change}$
1	0.2361	1.6604	1.3299	17.40%	100%
2	0.8319	0.7789	0.5456	7.20%	53.09%
3	0.8883	0.6350	0.4311	5.67%	18.48%
4	0.9044	0.5870	0.3960	5.17%	7.55%
5	0.9122	0.5622	0.3779	4.92%	4.23%
6	0.9181	0.5431	0.3675	4.78%	3.40%
7	0.9274	0.5113	0.3494	4.55%	5.85%
8	0.9304	0.5007	0.3480	4.56%	2.08%
9	0.9345	0.4861	0.3377	4.44%	2.91%
10	0.9359	0.4807	0.3366	4.43%	1.10%
11	0.9373	0.4757	0.3350	4.39%	1.05%

变量筛选结果如表3和图6所示，当逐步添加重要性的变量时，模型的RMSE、MAE、MAPE均有所下降，R²会有所上升。同时，为了考虑模型的复杂度和精度，根据$RMS{E}_{change}$ 来选择最优变量组合，有结果显示，当添加到第8个变量时，$RMS{E}_{change}=2.08\%<3\%$ ，可以认为第8个变量添加对模型的预测性能没有显著提升，可以舍弃这个变量。综上所述，可以得到变量筛选出的最优变量组合为组合7，即选择pH、水温、电导率、总磷、高锰酸盐指数、露点温度、气温这7个变量来对溶解氧(目标变量)进行预测。

Figure 6. Model performance of 11 variable combinations

图6. 11 种变量组合的模型性能

通过变量筛选，将氨氮、总氮、浊度和降雨量这四个变量筛除，在Pearson相关性分析结果中这四个变量与溶解氧的相关性也较低，因此在BART变量筛选中，得到了是与溶解氧(目标变量)具有相关性的变量组合。

5.2. 模型的构建与评估

通过变量筛选选择组合7 (pH、水温、电导率、总磷、高锰酸盐指数、露点温度、气温)来对溶解氧进行预测。在当前的数据集上选取70%数据为训练集，剩余的30%作为测试集，以保充分训练和有效评估模型性能。在训练集上分别构建BART、SVR、RF和XGB模型，并在测试集上计算R²、RMSE、MAE、MAPE来评估模型预测性能。该实验重复进行100次，对每次实验记录并统计R²、RMSE、MAE、MAPE，统计计算得出这四个指标的下四分位数(Q₁)、中位数(M_e)和上四分位数(Q₃)，得到的溶解氧预测结果如表4所示。

Table 4. Predicted results of dissolved oxygen

表4. 溶解氧预测结果

Method	R^2			RMSE			MAE			MAPE
	Q1	Me	Q3	Q1	Me	Q3	Q1	Me	Q3	Q1	Me	Q3
BART	0.933	0.936	0.939	0.469	0.480	0.492	0.327	0.333	0.337	4.28%	4.36%	4.43%
SVR	0.919	0.922	0.925	0.517	0.530	0.541	0.347	0.351	0.359	4.70%	4.76%	4.87%
RF	0.887	0.890	0.894	0.619	0.627	0.640	0.424	0.428	0.436	5.45%	5.53%	5.62%
XGB	0.887	0.890	0.894	0.616	0.628	0.639	0.445	0.451	0.458	5.58%	5.69%	5.78%

根据表4所示，在BART、SVR、RF和XGB模型中的结果对比来看，BART模型的在各指标的表现都会优于其他模型。具体而言，BART模型的RMSE、MAE、MAPE的中位数分别为0.480、0.333、4.36%，都低于其他三个模型；BART模型的R²的中位数为0.936，都高于其他三个模型。从图7结果对比显示来看，BART模型的预测精度会优于SVR、RF和XGB模型，表现出更好的拟合能力和稳定性。

Figure 7. The prediction fitting analysis of the BART, SVR, RF, and XGB models on the test set (the black solid line represents$y=x$ , and the dashed line indicates the 95% confidence interval of$y=x$ )

图7. BART、SVR、RF和XGB模型在测试集的预测拟合分析(黑色线段表示$y=x$ ，虚线段表示$y=x$ 的95%置信区间)

6. 结论

1) 本文基于对长江中游三个监测断面2021年1月1日到2022年10月28日的数据进行分析，通过BART进行变量筛选，从11个水质变量和气象变量，选择出对溶解氧影响最关键的7个变量，分别为pH、水温、电导率、总磷、高锰酸盐指数、露点温度、气温。其中pH和水温对溶解氧的影响较强。

2) 本文对BART、SVR、RF和XGB模型的稳健性进行研究，通过一系列的实验设计和模拟，详细研究探索了BART模型对于SVR、RF和XGB模型的优势。实验结果显示，BART模型相较于SVR、RF和XGB模型在面对异常值情况的预测性能更突出，具有强大的抗干扰性能，显示出强大的鲁棒性。将这四种模型运用于长江中游水质溶解氧的预测研究中，结果显示BART模型的预测性能会优于SVR、RF和XGB模型，说明BART模型适用于环境数据的建模，体现了BART模型具有实际应用价值。