洪水預報中的序列到序列模型及其可解釋性擴展

前些天發現了一個巨牛的人工智能學習網站，通俗易懂，風趣幽默，忍不住分享一下給大家，覺得好請收藏。點擊跳轉到網站。

1. 引言

洪水預報是水文科學和災害管理中的重要課題，準確的洪水預報可以顯著減少人員傷亡和經濟損失。近年來，深度學習技術在時間序列預測領域取得了顯著進展，特別是序列到序列(Seq2Seq)模型在各種預測任務中表現出色。本文將詳細介紹三種用于洪水預報的模型：基礎Seq2Seq模型、Seq2Seq-LRP模型和Seq2Seq-LRP-Attention模型，并分析它們的實現細節和性能特點。

2. 基礎Seq2Seq模型實現與問題分析

2.1 Seq2Seq模型架構

Seq2Seq模型由編碼器和解碼器兩部分組成，通常采用RNN、LSTM或GRU作為基礎單元。在洪水預報任務中，我們的輸入是歷史水文時間序列數據，輸出是未來一段時間的水位或流量預測。

import torch
import torch.nn as nn
from torch.autograd import Variableclass Encoder(nn.Module):def __init__(self, input_size, hidden_size, num_layers=1):super(Encoder, self).__init__()self.hidden_size = hidden_sizeself.num_layers = num_layersself.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True)def forward(self, x):# x shape: (batch_size, seq_len, input_size)outputs, (hidden, cell) = self.lstm(x)return hidden, cellclass Decoder(nn.Module):def __init__(self, input_size, hidden_size, output_size, num_layers=1):super(Decoder, self).__init__()self.hidden_size = hidden_sizeself.output_size = output_sizeself.num_layers = num_layersself.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True)self.fc = nn.Linear(hidden_size, output_size)def forward(self, x, hidden, cell):# x shape: (batch_size, 1, input_size)output, (hidden, cell) = self.lstm(x, (hidden, cell))prediction = self.fc(output.squeeze(1))return prediction, hidden, cellclass Seq2Seq(nn.Module):def __init__(self, encoder, decoder, device):super(Seq2Seq, self).__init__()self.encoder = encoderself.decoder = decoderself.device = devicedef forward(self, source, target, teacher_forcing_ratio=0.5):# source shape: (batch_size, seq_len, input_size)# target shape: (batch_size, output_seq_len, output_size)batch_size = target.shape[0]target_len = target.shape[1]target_dim = target.shape[2]# 初始化輸出張量outputs = torch.zeros(batch_size, target_len, target_dim).to(self.device)# 編碼器處理hidden, cell = self.encoder(source)# 第一個解碼器輸入decoder_input = source[:, -1, :].unsqueeze(1)  # 使用最后一個輸入作為初始解碼器輸入for t in range(target_len):# 解碼器一步預測output, hidden, cell = self.decoder(decoder_input, hidden, cell)# 存儲預測結果outputs[:, t, :] = output# 決定是否使用教師強制use_teacher_forcing = True if random.random() < teacher_forcing_ratio else Falseif use_teacher_forcing and t < target_len - 1:# 使用真實值作為下一個輸入decoder_input = target[:, t, :].unsqueeze(1)else:# 使用預測值作為下一個輸入decoder_input = output.unsqueeze(1)return outputs

2.2 當前實現的問題分析

在當前的Seq2Seq實現中，我們發現了幾個可能導致輸出目標錯誤的問題：

維度不匹配問題：解碼器的輸出維度可能沒有正確映射到目標維度。在洪水預報中，我們通常需要預測多個水文站點的水位，輸出維度應與目標站點數量一致。
教師強制策略問題：當前的教師強制策略可能在序列末端處理不當，導致預測偏差累積。
初始化解碼器輸入問題：簡單地使用最后一個輸入作為解碼器初始輸入可能不適合洪水預報場景，因為水位變化具有連續性特征。
缺失歸一化處理：水文數據通常需要適當的歸一化，當前實現中缺少這一關鍵步驟。

2.3 改進的Seq2Seq實現

針對上述問題，我們提出以下改進：

class ImprovedSeq2Seq(nn.Module):def __init__(self, encoder, decoder, device, output_dim):super(ImprovedSeq2Seq, self).__init__()self.encoder = encoderself.decoder = decoderself.device = deviceself.output_dim = output_dimself.scaler = None  # 用于數據歸一化def set_scaler(self, scaler):self.scaler = scalerdef forward(self, source, target=None, teacher_forcing_ratio=0.5, predict_mode=False):batch_size = source.shape[0]if target is not None:target_len = target.shape[1]else:target_len = 24  # 默認預測24小時# 初始化輸出張量outputs = torch.zeros(batch_size, target_len, self.output_dim).to(self.device)# 編碼器處理hidden, cell = self.encoder(source)# 改進的初始解碼器輸入：使用歷史數據的加權平均weights = torch.linspace(0.1, 1.0, steps=source.shape[1]).to(self.device)weighted_input = (source * weights.view(1, -1, 1)).sum(dim=1) / weights.sum()decoder_input = weighted_input.unsqueeze(1)for t in range(target_len):output, hidden, cell = self.decoder(decoder_input, hidden, cell)outputs[:, t, :] = outputif predict_mode:# 預測模式下不使用教師強制decoder_input = output.unsqueeze(1)else:# 訓練時使用教師強制use_teacher_forcing = True if random.random() < teacher_forcing_ratio else Falseif use_teacher_forcing and target is not None:decoder_input = target[:, t, :].unsqueeze(1)else:decoder_input = output.unsqueeze(1)return outputsdef predict(self, source, steps):self.eval()with torch.no_grad():if self.scaler is not None:source = self.scaler.transform(source)source = torch.FloatTensor(source).unsqueeze(0).to(self.device)prediction = self.forward(source, predict_mode=True)if self.scaler is not None:prediction = self.scaler.inverse_transform(prediction.cpu().numpy())return prediction.squeeze(0)

3. Seq2Seq-LRP模型實現

3.1 LRP(層級相關性傳播)原理

LRP是一種解釋神經網絡決策的方法，通過將輸出相關性反向傳播到輸入層，顯示每個輸入特征對預測結果的貢獻程度。在洪水預報中，這可以幫助我們理解哪些歷史水文特征對當前預測最重要。

3.2 Seq2Seq-LRP模型實現

class Seq2SeqLRP(nn.Module):def __init__(self, encoder, decoder, device):super(Seq2SeqLRP, self).__init__()self.seq2seq = Seq2Seq(encoder, decoder, device)self.device = devicedef forward(self, source, target=None, teacher_forcing_ratio=0.5):return self.seq2seq(source, target, teacher_forcing_ratio)def lrp_forward(self, x):# 保存所有層的輸入輸出self.activations = {}# 編碼器部分lstm = self.seq2seq.encoder.lstmh_0 = torch.zeros(lstm.num_layers, x.size(0), lstm.hidden_size).to(self.device)c_0 = torch.zeros(lstm.num_layers, x.size(0), lstm.hidden_size).to(self.device)# LSTM前向傳播seq_len = x.size(1)hidden_seq = []for t in range(seq_len):xt = x[:, t, :]gates = lstm.weight_ih @ xt.T + lstm.weight_hh @ h_0 + lstm.bias_ih.unsqueeze(1) + lstm.bias_hh.unsqueeze(1)i_gate = torch.sigmoid(gates[:lstm.hidden_size])f_gate = torch.sigmoid(gates[lstm.hidden_size:2*lstm.hidden_size])g_gate = torch.tanh(gates[2*lstm.hidden_size:3*lstm.hidden_size])o_gate = torch.sigmoid(gates[3*lstm.hidden_size:])c_1 = f_gate * c_0 + i_gate * g_gateh_1 = o_gate * torch.tanh(c_1)self.activations[f'encoder_{t}_i'] = i_gateself.activations[f'encoder_{t}_f'] = f_gateself.activations[f'encoder_{t}_g'] = g_gateself.activations[f'encoder_{t}_o'] = o_gateself.activations[f'encoder_{t}_c'] = c_1self.activations[f'encoder_{t}_h'] = h_1h_0, c_0 = h_1, c_1hidden_seq.append(h_1.T)hidden_seq = torch.stack(hidden_seq, dim=1)return hidden_seqdef compute_relevance(self, x, target_time_step=0):# 前向傳播收集激活值self.lrp_forward(x)# 初始化相關性R = torch.zeros_like(x)# 解碼器部分的反向傳播decoder_lstm = self.seq2seq.decoder.lstmdecoder_fc = self.seq2seq.decoder.fc# 從目標時間步開始反向傳播for t in reversed(range(x.size(1))):# 獲取當前時間步的激活值h = self.activations[f'encoder_{t}_h']c = self.activations[f'encoder_{t}_c']i = self.activations[f'encoder_{t}_i']f = self.activations[f'encoder_{t}_f']g = self.activations[f'encoder_{t}_g']o = self.activations[f'encoder_{t}_o']# 計算LSTM門的相關性R_h = torch.ones_like(h)  # 初始化# 反向傳播到輸入W_ih = decoder_lstm.weight_ihW_hh = decoder_lstm.weight_hh# 計算輸入和隱藏狀態的相關性R_x = (x[:, t, :].unsqueeze(1) * (W_ih @ R_h.T)).TR_h_prev = (h.unsqueeze(1) * (W_hh @ R_h.T)).T# 累積相關性R[:, t, :] = R_x# 更新隱藏狀態相關性R_h = R_h_prevreturn R

3.3 LRP在洪水預報中的應用

在洪水預報場景中，LRP可以幫助我們：

識別關鍵輸入特征：確定哪些歷史時間點的水位、降雨量等對當前預測影響最大。
模型可信度評估：當預測結果與LRP分析的關鍵特征不一致時，可以懷疑模型可能存在問題。
異常檢測：當LRP顯示不尋常的特征重要性分布時，可能表明輸入數據存在異常。

4. Seq2Seq-LRP-Attention模型實現

4.1 注意力機制與LRP的結合

注意力機制可以動態地為輸入序列的不同部分分配權重，而LRP可以解釋這些權重如何影響最終預測。結合兩者可以創建既強大又可解釋的洪水預報模型。

4.2 模型實現

class Attention(nn.Module):def __init__(self, enc_hidden_size, dec_hidden_size):super(Attention, self).__init__()self.enc_hidden_size = enc_hidden_sizeself.dec_hidden_size = dec_hidden_sizeself.attn = nn.Linear(enc_hidden_size + dec_hidden_size, dec_hidden_size)self.v = nn.Parameter(torch.rand(dec_hidden_size))def forward(self, hidden, encoder_outputs):# hidden shape: (batch_size, dec_hidden_size)# encoder_outputs shape: (batch_size, seq_len, enc_hidden_size)seq_len = encoder_outputs.shape[1]# 重復隱藏狀態以匹配序列長度hidden = hidden.unsqueeze(1).repeat(1, seq_len, 1)# 計算注意力能量energy = torch.tanh(self.attn(torch.cat((hidden, encoder_outputs), dim=2)))energy = energy.permute(0, 2, 1)# 計算注意力分數v = self.v.repeat(encoder_outputs.shape[0], 1).unsqueeze(1)attention = torch.bmm(v, energy).squeeze(1)return torch.softmax(attention, dim=1)class AttnDecoder(nn.Module):def __init__(self, input_size, hidden_size, output_size, encoder_hidden_size, num_layers=1):super(AttnDecoder, self).__init__()self.hidden_size = hidden_sizeself.output_size = output_sizeself.num_layers = num_layersself.attention = Attention(encoder_hidden_size, hidden_size)self.lstm = nn.LSTM(input_size + encoder_hidden_size, hidden_size, num_layers, batch_first=True)self.fc = nn.Linear(hidden_size, output_size)def forward(self, x, hidden, cell, encoder_outputs):# x shape: (batch_size, 1, input_size)# hidden shape: (num_layers, batch_size, hidden_size)# cell shape: (num_layers, batch_size, hidden_size)# encoder_outputs shape: (batch_size, seq_len, encoder_hidden_size)# 計算注意力權重attn_weights = self.attention(hidden[-1], encoder_outputs)# 計算上下文向量context = torch.bmm(attn_weights.unsqueeze(1), encoder_outputs).squeeze(1)# 連接輸入和上下文rnn_input = torch.cat((x.squeeze(1), context), dim=1).unsqueeze(1)# LSTM處理output, (hidden, cell) = self.lstm(rnn_input, (hidden, cell))# 全連接層prediction = self.fc(output.squeeze(1))return prediction, hidden, cell, attn_weightsclass Seq2SeqLRPAttention(nn.Module):def __init__(self, encoder, decoder, device):super(Seq2SeqLRPAttention, self).__init__()self.encoder = encoderself.decoder = decoderself.device = devicedef forward(self, source, target, teacher_forcing_ratio=0.5):batch_size = target.shape[0]target_len = target.shape[1]target_dim = target.shape[2]# 初始化輸出張量outputs = torch.zeros(batch_size, target_len, target_dim).to(self.device)attentions = torch.zeros(batch_size, target_len, source.shape[1]).to(self.device)# 編碼器處理encoder_outputs, (hidden, cell) = self.encoder(source)# 初始解碼器輸入decoder_input = source[:, -1, :].unsqueeze(1)for t in range(target_len):output, hidden, cell, attn_weights = self.decoder(decoder_input, hidden, cell, encoder_outputs)outputs[:, t, :] = outputattentions[:, t, :] = attn_weightsuse_teacher_forcing = True if random.random() < teacher_forcing_ratio else Falseif use_teacher_forcing and t < target_len - 1:decoder_input = target[:, t, :].unsqueeze(1)else:decoder_input = output.unsqueeze(1)return outputs, attentionsdef compute_lrp(self, x, target_time_step=0):# 前向傳播收集激活值和注意力權重outputs, attentions = self.forward(x, None, teacher_forcing_ratio=0)# 初始化相關性R = torch.zeros_like(x)# 獲取目標時間步的注意力權重attn = attentions[:, target_time_step, :]# 計算每個輸入時間步的相關性for t in range(x.size(1)):R[:, t, :] = x[:, t, :] * attn[:, t].unsqueeze(1)return R, attentions

4.3 模型優勢與應用

Seq2Seq-LRP-Attention模型結合了注意力機制和LRP解釋性技術，在洪水預報中具有以下優勢：

動態特征關注：注意力機制可以自動學習不同歷史時間點對當前預測的重要性。
可解釋的注意力：LRP可以解釋注意力權重的合理性，驗證模型是否關注了正確的特征。
多尺度分析：可以同時分析長期和短期水文模式對預測的影響。
不確定性量化：通過分析注意力分布的變化，可以評估預測結果的不確定性。

5. 模型訓練與評估

5.1 數據準備與預處理

洪水預報數據通常包括水位、降雨量、蒸發量等多個時間序列。我們需要進行以下預處理：

from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import train_test_split
import numpy as npdef prepare_flood_data(data, input_len=72, output_len=24, test_size=0.2):"""準備洪水預報數據集data: (samples, features) 形狀的numpy數組input_len: 輸入序列長度(小時)output_len: 輸出序列長度(小時)test_size: 測試集比例"""X, y = [], []for i in range(len(data) - input_len - output_len):X.append(data[i:i+input_len])y.append(data[i+input_len:i+input_len+output_len, 0])  # 假設第一列是目標水位X = np.array(X)y = np.array(y)# 歸一化scaler = MinMaxScaler()X = scaler.fit_transform(X.reshape(-1, X.shape[-1])).reshape(X.shape)y = scaler.fit_transform(y)# 劃分訓練測試集X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=test_size, shuffle=False)return X_train, X_test, y_train, y_test, scaler

5.2 訓練過程

import torch.optim as optim
from torch.utils.data import DataLoader, TensorDatasetdef train_model(model, X_train, y_train, X_val, y_val, scaler, epochs=100, batch_size=32, lr=0.001):device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')model = model.to(device)model.set_scaler(scaler)# 準備數據加載器train_data = TensorDataset(torch.FloatTensor(X_train), torch.FloatTensor(y_train))val_data = TensorDataset(torch.FloatTensor(X_val), torch.FloatTensor(y_val))train_loader = DataLoader(train_data, batch_size=batch_size, shuffle=True)val_loader = DataLoader(val_data, batch_size=batch_size, shuffle=False)# 定義損失函數和優化器criterion = nn.MSELoss()optimizer = optim.Adam(model.parameters(), lr=lr)best_val_loss = float('inf')train_losses = []val_losses = []for epoch in range(epochs):model.train()train_loss = 0for batch_x, batch_y in train_loader:batch_x, batch_y = batch_x.to(device), batch_y.to(device)optimizer.zero_grad()outputs = model(batch_x, batch_y)loss = criterion(outputs, batch_y)loss.backward()optimizer.step()train_loss += loss.item()# 驗證model.eval()val_loss = 0with torch.no_grad():for batch_x, batch_y in val_loader:batch_x, batch_y = batch_x.to(device), batch_y.to(device)outputs = model(batch_x)loss = criterion(outputs, batch_y)val_loss += loss.item()train_loss /= len(train_loader)val_loss /= len(val_loader)train_losses.append(train_loss)val_losses.append(val_loss)# 保存最佳模型if val_loss < best_val_loss:best_val_loss = val_losstorch.save(model.state_dict(), 'best_model.pth')print(f'Epoch {epoch+1}/{epochs}, Train Loss: {train_loss:.4f}, Val Loss: {val_loss:.4f}')return train_losses, val_losses

5.3 模型評估指標

在洪水預報中，我們通常使用以下指標評估模型性能：

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_scoredef evaluate_model(model, X_test, y_test, scaler):device = next(model.parameters()).devicemodel.eval()with torch.no_grad():test_data = TensorDataset(torch.FloatTensor(X_test), torch.FloatTensor(y_test))test_loader = DataLoader(test_data, batch_size=32, shuffle=False)predictions = []true_values = []for batch_x, batch_y in test_loader:batch_x, batch_y = batch_x.to(device), batch_y.to(device)outputs = model(batch_x)if scaler is not None:outputs = scaler.inverse_transform(outputs.cpu().numpy())batch_y = scaler.inverse_transform(batch_y.cpu().numpy())predictions.append(outputs)true_values.append(batch_y)predictions = np.concatenate(predictions, axis=0)true_values = np.concatenate(true_values, axis=0)# 計算評估指標mse = mean_squared_error(true_values, predictions)rmse = np.sqrt(mse)mae = mean_absolute_error(true_values, predictions)r2 = r2_score(true_values, predictions)# 計算Nash-Sutcliffe效率系數(水文模型常用指標)numerator = np.sum((true_values - predictions) ** 2)denominator = np.sum((true_values - np.mean(true_values)) ** 2)nse = 1 - (numerator / denominator)return {'MSE': mse,'RMSE': rmse,'MAE': mae,'R2': r2,'NSE': nse,'predictions': predictions,'true_values': true_values}