当前位置：首页 > news >正文

水质指数预测模型R²偏低的原因分析与优化策略

news 2025/6/27 7:55:58

问题背景

在最近的水质指数预测项目中，我们尝试了多种机器学习模型，包括线性回归、Ridge回归、Lasso回归、XGBoost、梯度提升回归和随机森林。然而，所有模型的R²值都不理想，最佳模型Ridge回归的R²仅为0.4058。本文将分析R²偏低的可能原因，并提出优化策略，同时展示各算法的核心代码实现。

模型性能概览

模型性能比较:
- Ridge回归: R² = 0.4058
- Lasso回归: R² = 0.4031
- 线性回归: R² = 0.3881
- 优化后随机森林: R² = 0.3312
- 梯度提升回归: R² = 0.3291
- 随机森林: R² = 0.2981
- XGBoost: R² = 0.2011

各算法实现代码

1. 线性回归

from sklearn.linear_model import LinearRegression# 多元线性回归regression_model = LinearRegression()regression_model.fit(X_train_pca, y_train)y_pred = regression_model.predict(X_test_pca)# 性能指标mse = mean_squared_error(y_test, y_pred)r2 = r2_score(y_test, y_pred)print(f'多元线性回归 - MSE: {mse:.4f}, R²: {r2:.4f}')

2. Ridge回归

from sklearn.linear_model import Ridge# Ridge回归ridge_model = Ridge(alpha=1.0)ridge_model.fit(X_train, y_train)y_pred_ridge = ridge_model.predict(X_test)r2_ridge = r2_score(y_test, y_pred_ridge)print(f'Ridge回归 - R²: {r2_ridge:.4f}')

3. Lasso回归

from sklearn.linear_model import Lasso# Lasso回归lasso_model = Lasso(alpha=0.1)lasso_model.fit(X_train, y_train)y_pred_lasso = lasso_model.predict(X_test)r2_lasso = r2_score(y_test, y_pred_lasso)print(f'Lasso回归 - R²: {r2_lasso:.4f}')

4. XGBoost

import xgboost as xgb# XGBoost模型xgb_model = xgb.XGBRegressor(random_state=15)xgb_model.fit(X_train, y_train)y_pred_xgb = xgb_model.predict(X_test)r2_xgb = r2_score(y_test, y_pred_xgb)print(f'XGBoost - R²: {r2_xgb:.4f}')

5. 梯度提升回归

from sklearn.ensemble import GradientBoostingRegressor# 梯度提升回归gbr_model = GradientBoostingRegressor(random_state=15)gbr_model.fit(X_train, y_train)y_pred_gbr = gbr_model.predict(X_test)r2_gbr = r2_score(y_test, y_pred_gbr)print(f'梯度提升回归 - R²: {r2_gbr:.4f}')

6. 随机森林

from sklearn.ensemble import RandomForestRegressor# 随机森林rf_model = RandomForestRegressor(random_state=15)rf_model.fit(X_train, y_train)y_pred_rf = rf_model.predict(X_test)# 性能指标mse_rf = mean_squared_error(y_test, y_pred_rf)r2_rf = r2_score(y_test, y_pred_rf)print(f'随机森林 - MSE: {mse_rf:.4f}, R²: {r2_rf:.4f}')

7. 随机森林参数优化

from sklearn.model_selection import RandomizedSearchCV# 参数搜索空间param_dist = {'n_estimators': [100, 300, 500, 700, 1000],'max_depth': [10, 20, 30, 40, 50, None],'min_samples_split': [2, 5, 10, 15],'min_samples_leaf': [1, 2, 4, 8],'max_features': ['sqrt', 'log2', None]}# 随机搜索random_search = RandomizedSearchCV(estimator=RandomForestRegressor(random_state=15),param_distributions=param_dist,n_iter=15,cv=5,verbose=1,random_state=17,n_jobs=-1)random_search.fit(X_train, y_train)best_params = random_search.best_params_# 使用最佳参数建模best_rf_model = RandomForestRegressor(**best_params, random_state=15)best_rf_model.fit(X_train, y_train)y_pred_best_rf = best_rf_model.predict(X_test)r2_best_rf = r2_score(y_test, y_pred_best_rf)

R²偏低的原因分析

1. 数据复杂性与样本量限制

水质系统是一个极其复杂的生态环境系统，受多种因素影响：

自然因素：季节变化、气温、降雨量、地质条件等

人为因素：工业排放、农业活动、城市污水等

生物因素：水生植物、微生物活动等

当前数据集可能无法完全捕捉这些复杂的相互作用，特别是在样本量有限的情况下。

2. 特征不足或相关性不强

从特征重要性分析可以看出：

最重要的5个特征:

1. DO_COD_Ratio: 0.2361

2. Dissolved_Oxygen_mg_L: 0.1616

3. Ammonia_N_mg_L: 0.1583

4. COD_mg_L: 0.1278

5. Total_Nitrogen_mg_L: 0.0567

虽然已经构建了一些派生特征（如DO_COD_Ratio），但现有特征可能仍不足以全面解释水质指数的变化。特征与目标变量之间的相关性可能不够强，导致预测能力有限。

3. 数据质量问题

水质监测数据可能存在以下问题：

测量误差：不同监测站点、不同时间的测量精度可能存在差异

缺失值：通过均值填充处理的缺失值可能引入偏差

异常值：尽管进行了异常值处理，但可能仍有未被识别的异常数据

4. 时空因素未充分考虑

水质具有明显的时空特性：

空间相关性：相邻监测站点的水质往往相互关联

时间序列特性：水质指标通常具有季节性、周期性变化

上下游关系：河流系统中的上游污染会影响下游水质

当前模型可能未能充分捕捉这些时空关系。

5. 聚类分析结果的启示

聚类分析将数据分为3个类别：

聚类分析总结:

聚类 0: 359个样本 (35.9%), 平均水质指数: 52.81

聚类 1: 234个样本 (23.4%), 平均水质指数: 58.66

聚类 2: 407个样本 (40.7%), 平均水质指数: 53.03

不同聚类的水质指数差异不大，这可能表明：

数据中的水质状况相对均质

现有特征难以区分不同水质状况

可能需要针对不同类型的水体分别建模

6. 模型选择与复杂度匹配

线性模型（Ridge、Lasso）表现优于非线性模型（随机森林、XGBoost），这可能表明：

水质指数与特征之间可能存在较强的线性关系

非线性模型可能过拟合了有限的数据

复杂模型的参数调优可能不够充分

优化策略

1. 数据增强

扩充数据来源：

# 整合多源数据def integrate_multiple_data_sources():# 整合气象数据weather_data = pd.read_csv('weather_data.csv')# 整合水文数据hydrology_data = pd.read_csv('hydrology_data.csv')# 整合土地利用数据land_use_data = pd.read_csv('land_use_data.csv')# 按时间和地点合并数据merged_data = pd.merge(df, weather_data, on=['Date', 'Province', 'Monitoring_Station'], how='left')merged_data = pd.merge(merged_data, hydrology_data, on=['Date', 'Province', 'Monitoring_Station'], how='left')merged_data = pd.merge(merged_data, land_use_data, on=['Province'], how='left')return merged_data

增加时间序列特征：

# 添加时间序列特征def add_temporal_features(df):# 确保日期列为datetime类型df['Date'] = pd.to_datetime(df['Date'])# 提取时间特征df['Year'] = df['Date'].dt.yeardf['Month'] = df['Date'].dt.monthdf['Day'] = df['Date'].dt.daydf['DayOfWeek'] = df['Date'].dt.dayofweekdf['Quarter'] = df['Date'].dt.quarter# 添加滞后特征for lag in [1, 3, 7, 14, 30]:for col in ['Dissolved_Oxygen_mg_L', 'COD_mg_L', 'Ammonia_N_mg_L']:df[f'{col}_lag_{lag}'] = df.groupby('Monitoring_Station')[col].shift(lag)# 添加移动平均特征for window in [7, 14, 30]:for col in ['Dissolved_Oxygen_mg_L', 'COD_mg_L', 'Ammonia_N_mg_L']:df[f'{col}_rolling_{window}'] = df.groupby('Monitoring_Station')[col].transform(lambda x: x.rolling(window=window, min_periods=1).mean())return df

2. 增强特征工程

添加交互特征：

# 创建交互特征def create_interaction_features(df):# 溶解氧与COD的交互df['DO_COD_Interaction'] = df['Dissolved_Oxygen_mg_L'] * df['COD_mg_L']# 氨氮与总磷的交互df['Ammonia_Phosphorus_Interaction'] = df['Ammonia_N_mg_L'] * df['Total_Phosphorus_mg_L']# 季节与溶解氧的交互for season in range(1, 5):df[f'DO_Season_{season}'] = df['Dissolved_Oxygen_mg_L'] * (df['Season'] == season)# 创建多项式特征from sklearn.preprocessing import PolynomialFeaturespoly = PolynomialFeatures(degree=2, include_bias=False)# 选择主要特征进行多项式转换main_features = ['Dissolved_Oxygen_mg_L', 'COD_mg_L', 'Ammonia_N_mg_L', 'Total_Phosphorus_mg_L']poly_features = poly.fit_transform(df[main_features])# 将多项式特征添加到数据框poly_feature_names = [f'poly_{i}' for i in range(poly_features.shape[1])]poly_df = pd.DataFrame(poly_features, columns=poly_feature_names)return pd.concat([df, poly_df], axis=1)

非线性变换：

# 应用非线性变换def apply_nonlinear_transformations(df):# 对数变换for col in ['COD_mg_L', 'Ammonia_N_mg_L', 'Total_Phosphorus_mg_L', 'Heavy_Metals_Pb_ug_L']:df[f'log_{col}'] = np.log1p(df[col])  # log1p处理零值# 平方根变换for col in ['COD_mg_L', 'Total_Nitrogen_mg_L']:df[f'sqrt_{col}'] = np.sqrt(df[col])# 倒数变换for col in ['Dissolved_Oxygen_mg_L', 'N_P_Ratio']:df[f'inverse_{col}'] = 1 / (df[col] + 1e-5)  # 添加小常数避免除零return df

3. 高级模型与集成学习

尝试更先进的模型：

# LightGBM模型import lightgbm as lgblgb_train = lgb.Dataset(X_train, y_train)lgb_eval = lgb.Dataset(X_test, y_test, reference=lgb_train)params = {'boosting_type': 'gbdt','objective': 'regression','metric': 'rmse','num_leaves': 31,'learning_rate': 0.05,'feature_fraction': 0.9,'bagging_fraction': 0.8,'bagging_freq': 5,'verbose': 0}lgb_model = lgb.train(params,lgb_train,num_boost_round=1000,valid_sets=lgb_eval,early_stopping_rounds=100,verbose_eval=100)y_pred_lgb = lgb_model.predict(X_test, num_iteration=lgb_model.best_iteration)r2_lgb = r2_score(y_test, y_pred_lgb)print(f'LightGBM - R²: {r2_lgb:.4f}')

构建Stacking集成模型：

from sklearn.ensemble import StackingRegressorfrom sklearn.linear_model import ElasticNet# 基础模型base_models = [('ridge', Ridge(alpha=1.0)),('lasso', Lasso(alpha=0.1)),('rf', RandomForestRegressor(random_state=15, n_estimators=300)),('gbr', GradientBoostingRegressor(random_state=15))]# 元模型meta_model = ElasticNet(alpha=0.1, l1_ratio=0.5)# 构建Stacking模型stacking_model = StackingRegressor(estimators=base_models,final_estimator=meta_model,cv=5)# 训练模型stacking_model.fit(X_train, y_train)y_pred_stacking = stacking_model.predict(X_test)r2_stacking = r2_score(y_test, y_pred_stacking)print(f'Stacking集成 - R²: {r2_stacking:.4f}')

4. 分层建模策略

根据聚类结果为不同类型的水体构建专用模型：

# 分层建模def stratified_modeling(df, features, target):results = {}# 按聚类分组训练模型for cluster in df['Cluster'].unique():# 获取当前聚类的数据cluster_data = df[df['Cluster'] == cluster]X_cluster = cluster_data[features]y_cluster = cluster_data[target]# 划分训练测试集X_train_c, X_test_c, y_train_c, y_test_c = train_test_split(X_cluster, y_cluster, test_size=0.2, random_state=15)# 标准化scaler_c = StandardScaler()X_train_c_scaled = scaler_c.fit_transform(X_train_c)X_test_c_scaled = scaler_c.transform(X_test_c)# 训练模型ridge_model_c = Ridge(alpha=1.0)ridge_model_c.fit(X_train_c_scaled, y_train_c)# 评估模型y_pred_c = ridge_model_c.predict(X_test_c_scaled)r2_c = r2_score(y_test_c, y_pred_c)results[f'Cluster_{cluster}'] = {'model': ridge_model_c,'scaler': scaler_c,'r2': r2_c,'sample_size': len(cluster_data)}return results# 执行分层建模stratified_models = stratified_modeling(df, features, 'Water_Quality_Index')# 输出各聚类模型性能for cluster, result in stratified_models.items():print(f"{cluster} - 样本数: {result['sample_size']}, R²: {result['r2']:.4f}")

5. 深度学习方法

import tensorflow as tffrom tensorflow.keras.models import Sequentialfrom tensorflow.keras.layers import Dense, Dropout, BatchNormalizationfrom tensorflow.keras.callbacks import EarlyStopping# 构建神经网络模型def build_nn_model(input_dim):model = Sequential([Dense(64, activation='relu', input_dim=input_dim),BatchNormalization(),Dropout(0.3),Dense(32, activation='relu'),BatchNormalization(),Dropout(0.2),Dense(16, activation='relu'),Dense(1)])model.compile(optimizer='adam', loss='mse', metrics=['mae'])return model# 训练神经网络nn_model = build_nn_model(X_train.shape[1])early_stopping = EarlyStopping(monitor='val_loss',patience=20,restore_best_weights=True)history = nn_model.fit(X_train, y_train,epochs=200,batch_size=32,validation_split=0.2,callbacks=[early_stopping],verbose=1
)# 评估模型
y_pred_nn = nn_model.predict(X_test)
r2_nn = r2_score(y_test, y_pred_nn)
print(f'深度神经网络 - R²: {r2_nn:.4f}')# 可视化训练过程
plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.plot(history.history['loss'], label='训练损失')
plt.plot(history.history['val_loss'], label='验证损失')
plt.title('模型损失')
plt.xlabel('轮次')
plt.ylabel('损失')
plt.legend()plt.subplot(1, 2, 2)
plt.plot(history.history['mae'], label='训练MAE')
plt.plot(history.history['val_mae'], label='验证MAE')
plt.title('平均绝对误差')
plt.xlabel('轮次')
plt.ylabel('MAE')
plt.legend()
plt.tight_layout()
plt.show()

6. 时间序列专用模型

针对水质数据的时间序列特性，尝试专用模型：

# 准备时间序列数据def prepare_time_series_data(df, target_col, sequence_length=10):# 按监测站点和日期排序df = df.sort_values(['Monitoring_Station', 'Date'])X_seq, y_seq = [], []# 为每个监测站点创建序列for station in df['Monitoring_Station'].unique():station_data = df[df['Monitoring_Station'] == station]if len(station_data) < sequence_length + 1:continue# 提取特征列feature_cols = [col for col in df.columns if col not in ['Date', 'Province', 'Monitoring_Station', 'Cluster', target_col]]# 创建序列for i in range(len(station_data) - sequence_length):X_seq.append(station_data[feature_cols].iloc[i:i+sequence_length].values)y_seq.append(station_data[target_col].iloc[i+sequence_length])return np.array(X_seq), np.array(y_seq)# 构建LSTM模型from tensorflow.keras.layers import LSTMdef build_lstm_model(input_shape):model = Sequential([LSTM(64, return_sequences=True, input_shape=input_shape),Dropout(0.2),LSTM(32),Dropout(0.2),Dense(16, activation='relu'),Dense(1)])model.compile(optimizer='adam', loss='mse', metrics=['mae'])return model# 准备数据并训练模型X_seq, y_seq = prepare_time_series_data(df, 'Water_Quality_Index', sequence_length=10)# 划分训练测试集X_train_seq, X_test_seq, y_train_seq, y_test_seq = train_test_split(X_seq, y_seq, test_size=0.2, random_state=15)# 构建并训练LSTM模型lstm_model = build_lstm_model((X_train_seq.shape[1], X_train_seq.shape[2]))early_stopping = EarlyStopping(monitor='val_loss',patience=20,restore_best_weights=True)lstm_history = lstm_model.fit(X_train_seq, y_train_seq,epochs=100,batch_size=32,validation_split=0.2,callbacks=[early_stopping],verbose=1)# 评估模型y_pred_lstm = lstm_model.predict(X_test_seq)r2_lstm = r2_score(y_test_seq, y_pred_lstm)print(f'LSTM模型 - R²: {r2_lstm:.4f}')

结果对比与分析

实施上述优化策略后，我们可以对比各模型性能：

# 汇总所有模型性能all_models = {'原始Ridge回归': 0.4058,'原始Lasso回归': 0.4031,'原始线性回归': 0.3881,'原始优化随机森林': 0.3312,'特征工程后Ridge': r2_ridge_enhanced,  # 假设这是增强特征后的结果'LightGBM': r2_lgb,'Stacking集成': r2_stacking,'深度神经网络': r2_nn,'LSTM时间序列': r2_lstm}# 可视化对比plt.figure(figsize=(12, 6))models = list(all_models.keys())scores = list(all_models.values())# 按性能排序sorted_indices = np.argsort(scores)[::-1]models = [models[i] for i in sorted_indices]scores = [scores[i] for i in sorted_indices]plt.barh(models, scores, color='skyblue')plt.xlabel('R²值')plt.title('各模型性能对比')plt.xlim(0, 1)plt.grid(axis='x', linestyle='--', alpha=0.7)for i, v in enumerate(scores):plt.text(v + 0.01, i, f'{v:.4f}', va='center')plt.tight_layout()plt.show()