问题背景
在最近的水质指数预测项目中，我们尝试了多种机器学习模型，包括线性回归、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.year df['Month'] = df['Date'].dt.month df['Day'] = df['Date'].dt.day df['DayOfWeek'] = df['Date'].dt.dayofweek df['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 PolynomialFeatures poly = 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 lgb lgb_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 StackingRegressor from 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 tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense, Dropout, BatchNormalization
from 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 LSTM
def 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()结论与建议
主要发现
- 数据复杂性是关键挑战：水质系统涉及多种因素的复杂交互，单一的水质参数难以完全解释水质指数的变化。
- 时空特性需要特别关注：水质数据具有明显的时间序列特性和空间相关性，传统回归模型难以充分捕捉这些特性。
- 特征工程的重要性：通过创建交互特征、时间序列特征和非线性变换，可以显著提升模型性能。
- 分层建模的优势：针对不同类型的水体或不同聚类构建专用模型，可以更好地捕捉各自的特性。
- 高级模型的潜力：深度学习和时间序列专用模型在处理复杂水质数据时展现出良好潜力。
实践建议
- 数据收集与整合：
- 增加监测频率，特别是在水质变化敏感时期
- 整合气象、水文、土地利用等多源数据
- 建立完整的上下游监测网络
- 特征工程策略：
- 重视时间序列特征的构建
- 考虑特征之间的交互作用
- 针对不同污染物特性应用适当的非线性变换
- 模型选择与优化：
- 对于小型数据集，优先考虑正则化线性模型
- 数据充足时，尝试集成学习和深度学习方法
- 针对时空数据特性，考虑专用模型（LSTM、GNN等）
- 实施路径：
- 从基础模型开始，逐步增加复杂度
- 持续验证模型在不同条件下的稳定性
- 结合领域知识解释模型预测结果
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
