2025.3.23机器学习笔记：文献阅读

题目信息

• 题目： Enhancement of Hydrological Time Series Prediction with Real-World Time Series Generative Adversarial Network-Based Synthetic Data and Deep Learning Models
• 期刊： Scientific Reports
• 作者： Ana Dodiga, Vladimir Stankovic, Lina Stankovic, Milan Stojkovic
• 发表时间： 2023
• 文章链接：https://papers.ssrn.com/sol3/papers.cfm?abstract_id=5022039

摘要

河流水质直接影响人类健康、生态系统、生物多样性以及社会和工业需求。准确预测溶解氧（DO）和水温（WT）等关键水质参数，能为决策者提供预警。但是高质量的连续时间序列水质数据稀缺且获取的成本高，这导致了数据存在不确定性和缺失值，限制了机器学习的预测技术的实施。传统时间序列数据生成模型难以捕捉数据复杂依赖关系，利用生成对抗网络（GAN）进行年度水时间序列数据生成并提高预测模型准确性的研究较少。本论文旨在探究利用水流量数据辅助生成水质和水量数据，并分析生成数据对深度学习预测模型的影响，以提高河流水质时间序列预测的准确性和可靠性。

Abstract

River water quality directly impacts human health, ecosystems, biodiversity, and social and industrial needs. Accurate prediction of key water quality parameters such as dissolved oxygen (DO) and water temperature (WT) can provide early warnings for decision-makers. However, high-quality, high-resolution continuous time-series water quality data are scarce and costly to obtain, leading to uncertainties and missing values in the data, which limit the implementation of machine learning prediction techniques. Traditional time-series data generation models struggle to capture the complex dependencies in the data, and there is limited research on using Generative Adversarial Networks (GANs) for generating annual water time-series data and improving the accuracy of prediction models. This thesis aims to explore the use of water flow data to assist in generating water quality and quantity data, and to analyze the impact of the generated data on deep learning prediction models, with the goal of enhancing the accuracy and reliability of river water quality time-series forecasting.

创新点

以往研究多基于现有水质数据预测，未考虑易获取的水量数据，且少有使用GAN进行年度水时间序列数据生成以提高预测模型准确性的研究。本研究利用RTSGAN生成水质和水量的合成数据，并结合水质和水量数据进行预测。

网络架构

在数据输入前，作者对数据进行了预处理，预处理的方式分别为KNN与Bootstrap
KNN:

Bootstrap:
放回抽样，比如上述例子中，有可能第一次抽到 ①，第二次也抽到 ①，第三次还是抽到 ①，对三个样本 X = 1 , 2 , 3 ，放回抽样 1000 次，用这 1000 次抽样的结果分别计算出各自的均值 $\bar{x}$ , 然后再利用这 1000 个值，计算其整体的均值，方差等统计量，这就是 Bootstrap 方法
比如，有放回随机抽样，生成B个子样本（B为采样次数，如B=100）。每次抽样保持样本量不变，但数据点可能重复或缺失。

模型结构如下图所示：

1. TRAINING DATA为训练数据，包含三个站点的2013-2020共8年每日的水文数据（水流量、含氧量、水温），其用于RTSGAN的输入，用于学习时间序列的分布和模式。

2. RTSGAN是一个生成对抗网络，用于生成合成时间序列数据。其包括编码器和解码器
   编码器将输入时间序列数据如，水流量、含氧量、水温编码为固定维度的潜在向量，并捕捉时序和全局特征；解码器则是从潜在向量自回归生成时间序列。生成器是MLP，生成潜在向量；判别器是三层全连接网络，激活函数为LeakyReLU，区分真实和生成数据。
   优化目标为：
   $$\underset{G}{\min }\underset{D}{\max }{\mathbb{E}}_{r\sim \mathrm{encoder}(X,y)}[D(r)]-{\mathbb{E}}_{z\sim p(z)}[D(G(z))]$$
   关于Wasserstein与RTSGAN可以看我之前的博客解释：https://blog.csdn.net/Zcymatics/article/details/145981612?spm=1001.2014.3001.5501

3. SLIDING WINDOW用于分割数据，用于将时间序列数据分割为输入-输出对。
   在数据生成之前，使用滑动窗口方法将水质和数量时间序列数据通过归一化和分割预处理成三维(3D)数据集。分割过程如下：
   （1）选择一个窗口大小为w的窗口，并沿着长度为n的时间序列数据滑动。
   （2）在每一步中，生成一个大小为[n, w]的二维矩阵（其中n表示特征数乘以节点数，因为模型一次为所有三个站点生成数据。给定N作为时间序列的总长度(在我们的例子中，8年的数据为2920天)，N = 9(3个特征和3个节点），窗口大小w =365和滑动步长i= 1。
   （3）这就产生了大小为[n, w] =[9,365]的重叠子序列。生成的序列总数为[(N -w)/i + 1] = 2556。
   这些序列作为RTSGAN的输入，用于合成复制原始数据中时间模式的附加序列。

4. 经过数据增强后的数据用于两个模型中，其中，LSTM为每个站训练一个模型，输入10天数据，预测3天DO和WT；GNN训练一个模型，利用空间依赖性，同时处理三个站的数据，其中节点为三个站点，边权重为站间距离。

5. 最后在测试数据（2年）上评估模型，计算RMSE和NSE。
   $\text{NSE}=1-\frac{{\sum }_{i=1}^n(y_t-y_t^{\prime }{)}^2}{{\sum }_{t=1}^n(y_t-\bar{y}{)}^2}$，其中y_t为时间步长t的观测数据，y_t’为预测数据，y为观测数据的均值，n为观测总数。
   

实验

论文围绕利用RTSGAN生成合成数据结合深度学习模型（LSTM和GNN）进行水质时间序列预测展开实验，具体结果如下：

1. 合成数据生成性能

通过PCA和t - SNE对合成数据和原始数据进行对比，结果显示合成数据在不同降维视角下与原始数据保持相似性，整体模式复制良好，表明数据生成过程有效。

	PCA：捕捉数据的线性结构和主要方差，验证合成数据是否保留了原始数据的全局特性。
	t-SNE：关注数据的局部结构和非线性关系，检查合成数据在细节上与原始数据的相似性。

PCA和t - SNE分析如下图所示：

对溶解氧（DO）、水温（WT）和流量三个关键水质和水量参数的分布进行可视化对比，包括密度图、直方图以及均值、中位数和峰度等统计指标。流量数据方面，合成数据集与原始数据总体吻合，仅存在细微差异；水温数据分布在三个站点均与原始数据高度一致；DO数据分布在三个站点也与原始数据极为相似，合成数据的均值、中位数与原始数据对齐，峰度值接近，说明合成数据能有效用于后续分析和建模。

2. 预测模型性能
在三个水文站对不同LSTM模型。如下图所示：
其中，LSTMsynX，X表示添加到训练集的合成序列年数，syn表示是合成数据进行评估。结果表明，纳入合成数据的模型在各站点均显著优于仅使用原始数据的模型。
在诺维萨德站，LSTMsyn12的DO预测RMSE从LSTMsyn0的0.0601降至0.0575，改善4.32%；
LSTMsyn8的水温预测RMSE从0.0285降至0.0264，改善7.37%。在森塔站，LSTMsyn20的DO和水温预测RMSE分别降至0.0276和0.0133，改善率分别为9.51%和14.74%。
在泽蒙站，LSTMsyn16的水温预测RMSE从0.0195降至0.0174，改善10.77%。
总体而言，合成数据显著提高了LSTM模型在各站点对DO和WT的预测准确性，使用较多合成数据的模型表现更优，改善幅度约为4% - 15%。

对不同GNN模型在三个站点的水质参数预测性能如下图所示：
在诺维萨德站，GNNsyn8在第1天的DO预测中表现最佳，RMSE为0.0495，改善1.20%；
GNNsyn12和GNNsyn16在水温预测上有显著改善，第1天RMSE降至0.0175，改善19.35%。
在森塔站，GNNsyn20在第1天和第2天的DO预测中表现最佳，RMSE分别为0.0269和0.0330，分别提升了17.23%和9.84%；
GNNsyn16在水温预测上改良比较大，第1天RMSE降至0.0120，提升了43.40%。
在泽蒙站，基线模型GNNsyn0在DO和水温预测上表现最佳，合成数据未提升预测准确性。
综合来看，合成数据在某些站点和参数上显著提高了GNN模型的预测准确性，但效果因位置和参数而异。

代码如下：

import numpy as np
import pandas as pd
from sklearn.neighbors import KNeighborsRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_squared_error
import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
from torch_geometric.nn import GATConv
from torch_geometric.data import Data
import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense

# 设置随机种子以确保可重复性
np.random.seed(42)
torch.manual_seed(42)
tf.random.set_seed(42)

# 1. 数据预处理：生成模拟数据并填补缺失值
# 模拟数据：Novi Sad完整，Senta和Zemun缺失DO
dates = pd.date_range("2013-01-01", "2022-12-31", freq="D")  # 3650天
n_days = len(dates)

# 模拟Novi Sad数据
n_flow = np.random.normal(3000, 500, n_days)  # 流量
n_wt = np.random.normal(15, 5, n_days)  # 水温
n_do = np.random.normal(10, 1, n_days)  # DO

# 模拟Senta数据
s_flow = np.random.normal(600, 100, n_days)
s_wt = np.random.normal(15, 5, n_days)
s_do = np.full(n_days, np.nan)  # DO初始化为NaN
s_do[::30] = np.random.normal(10, 1, n_days//30)  # 每月1次

# 模拟Zemun数据
z_flow = np.random.normal(3500, 600, n_days)
z_wt = np.random.normal(15, 5, n_days)
z_do = np.full(n_days, np.nan)
z_do[::30] = np.random.normal(10, 1, n_days//30)

# 合并数据
data = pd.DataFrame({
    "N_Flow": n_flow, "N_WT": n_wt, "N_DO": n_do,
    "S_Flow": s_flow, "S_WT": s_wt, "S_DO": s_do,
    "Z_Flow": z_flow, "Z_WT": z_wt, "Z_DO": z_do
}, index=dates)

# KNN+Bootstrap填补缺失值
def knn_bootstrap_impute(data, target_col, feature_cols, n_bootstrap=100, k=5):
    train_data = data[data[target_col].notna()][feature_cols + [target_col]]
    test_data = data[data[target_col].isna()][feature_cols]
    scaler = StandardScaler()
    X_train = scaler.fit_transform(train_data[feature_cols])
    y_train = train_data[target_col].values
    X_test = scaler.transform(test_data)
    predictions = []
    n_samples = len(X_train)
    for _ in range(n_bootstrap):
        indices = np.random.choice(n_samples, n_samples, replace=True)
        X_boot = X_train[indices]
        y_boot = y_train[indices]
        knn = KNeighborsRegressor(n_neighbors=k, weights="distance")
        knn.fit(X_boot, y_boot)
        pred = knn.predict(X_test)
        predictions.append(pred)
    final_pred = np.mean(predictions, axis=0)
    data.loc[data[target_col].isna(), target_col] = final_pred
    return data

# 填补Senta和Zemun的DO
data = knn_bootstrap_impute(
    data, target_col="S_DO", 
    feature_cols=["N_Flow", "N_WT", "N_DO", "S_Flow", "S_WT", "Z_Flow", "Z_WT"],
    n_bootstrap=100, k=5
)
data = knn_bootstrap_impute(
    data, target_col="Z_DO", 
    feature_cols=["N_Flow", "N_WT", "N_DO", "S_Flow", "S_WT", "Z_Flow", "Z_WT"],
    n_bootstrap=100, k=5
)

# 标准化数据
scaler = StandardScaler()
data_scaled = scaler.fit_transform(data)
data_scaled = pd.DataFrame(data_scaled, columns=data.columns, index=data.index)

# 2. RTSGAN生成合成数据
class Generator(nn.Module):
    def __init__(self, input_dim, hidden_dim, output_dim):
        super(Generator, self).__init__()
        self.net = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.LeakyReLU(0.2),
            nn.LayerNorm(hidden_dim),
            nn.Linear(hidden_dim, output_dim)
        )
    
    def forward(self, z):
        return self.net(z)

class Discriminator(nn.Module):
    def __init__(self, input_dim, hidden_dim):
        super(Discriminator, self).__init__()
        self.net = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.LeakyReLU(0.2),
            nn.Linear(hidden_dim, 1)
        )
    
    def forward(self, x):
        return self.net(x)

def train_rtsgan(data, window_size=365, n_synthetic=2556, epochs=100):
    device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
    data_tensor = torch.FloatTensor(data.values).to(device)
    
    # 第一次滑动窗口
    inputs = []
    for i in range(len(data) - window_size):
        inputs.append(data_tensor[i:i+window_size])
    inputs = torch.stack(inputs)
    
    # 参数
    input_dim = window_size * data.shape[1]  # 365 * 9
    hidden_dim = 128
    output_dim = input_dim
    batch_size = 64
    
    # 模型
    G = Generator(input_dim, hidden_dim, output_dim).to(device)
    D = Discriminator(input_dim, hidden_dim).to(device)
    optimizer_G = optim.Adam(G.parameters(), lr=0.0002, betas=(0.5, 0.999))
    optimizer_D = optim.Adam(D.parameters(), lr=0.0002, betas=(0.5, 0.999))
    
    # 训练
    for epoch in range(epochs):
        for i in range(0, len(inputs), batch_size):
            real_data = inputs[i:i+batch_size].view(-1, input_dim)
            batch_size_actual = real_data.size(0)
            
            # 训练判别器
            z = torch.randn(batch_size_actual, input_dim).to(device)
            fake_data = G(z)
            real_score = D(real_data)
            fake_score = D(fake_data.detach())
            d_loss = -torch.mean(real_score) + torch.mean(fake_score)
            
            optimizer_D.zero_grad()
            d_loss.backward()
            optimizer_D.step()
            
            # 训练生成器
            fake_score = D(fake_data)
            g_loss = -torch.mean(fake_score)
            
            optimizer_G.zero_grad()
            g_loss.backward()
            optimizer_G.step()
        
        if epoch % 10 == 0:
            print(f"Epoch {epoch}, D Loss: {d_loss.item():.4f}, G Loss: {g_loss.item():.4f}")
    
    # 生成合成数据
    synthetic_data = []
    for _ in range(n_synthetic):
        z = torch.randn(1, input_dim).to(device)
        fake = G(z).view(window_size, data.shape[1])
        synthetic_data.append(fake.cpu().detach().numpy())
    
    return np.concatenate(synthetic_data, axis=0)

# 生成合成数据（2556个序列，约8年）
synthetic_data = train_rtsgan(data_scaled, window_size=365, n_synthetic=2556, epochs=100)
synthetic_data = pd.DataFrame(synthetic_data, columns=data.columns)

# 3. 合并数据并进行第二次滑动窗口
combined_data = pd.concat([data_scaled, synthetic_data], axis=0, ignore_index=True)

def create_sliding_window(data, input_window=10, output_window=3):
    X, y = [], []
    target_cols = ["N_DO", "N_WT", "S_DO", "S_WT", "Z_DO", "Z_WT"]
    for i in range(len(data) - input_window - output_window + 1):
        X.append(data.iloc[i:i+input_window].values)  # 10天输入
        y.append(data.iloc[i+input_window:i+input_window+output_window][target_cols].values)  # 3天输出
    return np.array(X), np.array(y)

X, y = create_sliding_window(combined_data, input_window=10, output_window=3)

# 划分训练和测试集
train_size = int(0.8 * len(X))
X_train, y_train = X[:train_size], y[:train_size]
X_test, y_test = X[train_size:], y[train_size:]

# 4. LSTM模型（以Novi Sad站为例）
def build_lstm(input_shape, output_shape):
    model = Sequential([
        LSTM(64, input_shape=input_shape, return_sequences=True),
        LSTM(32),
        Dense(output_shape[0] * output_shape[1], activation="linear")
    ])
    model.compile(optimizer="adam", loss="mse")
    return model

# 准备Novi Sad数据
X_train_n = X_train[:, :, [0, 1, 2]]  # Novi Sad: 流量, WT, DO
y_train_n = y_train[:, :, [0, 1]]  # Novi Sad: DO, WT
X_test_n = X_test[:, :, [0, 1, 2]]
y_test_n = y_test[:, :, [0, 1]]

# 训练LSTM
lstm_model = build_lstm(input_shape=(10, 3), output_shape=(3, 2))
lstm_model.fit(X_train_n, y_train_n.reshape(-1, 3*2), epochs=50, batch_size=32, verbose=1)

# 预测
y_pred_n = lstm_model.predict(X_test_n).reshape(-1, 3, 2)

# 5. GNN模型
# 构建图结构：Novi Sad->Senta, Novi Sad->Zemun
edge_index = torch.tensor([[0, 1, 0, 2], [1, 0, 2, 0]], dtype=torch.long)
edge_weight = torch.tensor([1.0, 1.0, 1.0, 1.0], dtype=torch.float)

class GNNModel(nn.Module):
    def __init__(self, input_dim, hidden_dim, output_dim):
        super(GNNModel, self).__init__()
        self.gat1 = GATConv(input_dim, hidden_dim, heads=2, dropout=0.2)
        self.gat2 = GATConv(hidden_dim*2, hidden_dim, heads=1, dropout=0.2)
        self.fc = nn.Linear(hidden_dim*10, output_dim)
    
    def forward(self, x, edge_index, edge_weight):
        # x: [n_nodes, time_steps, features]
        x = x.view(-1, x.size(-1))  # [n_nodes*time_steps, features]
        x = F.relu(self.gat1(x, edge_index, edge_weight))
        x = self.gat2(x, edge_index, edge_weight)
        x = x.view(-1, 10*32)  # [batch, time_steps*hidden]
        x = self.fc(x)
        return x.view(-1, 3, 6)  # [batch, time_steps, targets]

def prepare_gnn_data(X, y):
    data_list = []
    for i in range(len(X)):
        # 3站×10天×3特征
        x = torch.FloatTensor(X[i].reshape(3, 10, 3))
        y_i = torch.FloatTensor(y[i])  # 3天×6目标
        data = Data(x=x, edge_index=edge_index, edge_weight=edge_weight, y=y_i)
        data_list.append(data)
    return data_list

# 准备GNN数据
train_data = prepare_gnn_data(X_train, y_train)
test_data = prepare_gnn_data(X_test, y_test)

# 训练GNN
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
gnn_model = GNNModel(input_dim=3, hidden_dim=32, output_dim=3*6).to(device)
optimizer = optim.Adam(gnn_model.parameters(), lr=0.001)

for epoch in range(50):
    gnn_model.train()
    total_loss = 0
    for data in train_data:
        data = data.to(device)
        optimizer.zero_grad()
        out = gnn_model(data.x, data.edge_index, data.edge_weight)
        loss = F.mse_loss(out, data.y)
        loss.backward()
        optimizer.step()
        total_loss += loss.item()
    if epoch % 10 == 0:
        print(f"Epoch {epoch}, GNN Loss: {total_loss/len(train_data):.4f}")

# 预测
gnn_model.eval()
y_pred_gnn = []
with torch.no_grad():
    for data in test_data:
        data = data.to(device)
        pred = gnn_model(data.x, data.edge_index, data.edge_weight)
        y_pred_gnn.append(pred.cpu().numpy())
y_pred_gnn = np.array(y_pred_gnn)

# 6. 评估
def calculate_rmse(y_true, y_pred):
    rmse = np.sqrt(mean_squared_error(y_true.reshape(-1), y_pred.reshape(-1)))
    return rmse

# 计算NSE
def calculate_nse(y_true, y_pred):
    y_true_flat = y_true.reshape(-1)
    y_pred_flat = y_pred.reshape(-1)
    y_mean = np.mean(y_true_flat)
    numerator = np.sum((y_true_flat - y_pred_flat) ** 2)
    denominator = np.sum((y_true_flat - y_mean) ** 2)
    return 1 - numerator / denominator

# LSTM评估
rmse_lstm = calculate_rmse(y_test_n, y_pred_n)
nse_lstm = calculate_nse(y_test_n, y_pred_n)
print(f"LSTM RMSE: {rmse_lstm:.4f}, NSE: {nse_lstm:.4f}")

# GNN评估
rmse_gnn = calculate_rmse(y_test, y_pred_gnn)
nse_gnn = calculate_nse(y_test, y_pred_gnn)
print(f"GNN RMSE: {rmse_gnn:.4f}, NSE: {nse_gnn:.4f}")

不足以及展望

研究局限于三个水文站，增加站点时，LSTM方法资源消耗大、可扩展性差；Senta和Zemun站的溶解氧数据仅约每月一次，影响RTSGAN性能评估；合成数据并非在所有情况下都能提升预测效果，生成方法有待优化；训练GAN和深度学习模型需大量计算资源，不利于大规模实时应用。后续可将该工作流程应用于更多水文站，纳入气候数据等环境因素，以丰富训练数据集，提高模型泛化能力。