1140812 meeting

發佈於 2025-08-08 更新於 2025-08-09 約 7882 字預計閱讀 16 分鐘次閱讀

前言

由於人類活動與自然氣候變化（例如聖嬰現象）共同作用，全球表面溫度呈現長期上升趨勢。衛星遙測資料（如 GES DISC 的 MERRA-2 tavg1_2d_flx_Nx）雖能提供高頻率的全球表面溫度觀測，但仍受衛星掃描間隔、雲層遮擋、儀器故障與資料傳輸錯誤等因素影響，常出現資料空缺或異常值。為建立可支援政策制定與科學研究的可靠時空溫度場，本研究設計嚴謹的模擬實驗，比較三種時間序列預測方法： SSSD 、 TSMixer 與 RegressionEnsemble ，並結合 autoFRK 進行空間插值與預測，評估不同方法在資料填補與預測上的表現差異。

研究動機

全球表面溫度的持續上升已成為氣候變遷研究與政策制定的核心議題。高解析度且高時頻率的全球溫度資料，對於即時監測極端氣候事件、評估氣候模型準確性及制定減緩與調適策略至關重要。然而，衛星遙測資料在實際應用中，常因觀測限制與技術因素導致資料不完整或失真，若不加以處理，將影響溫度場重建與未來趨勢預測的可靠性。現有的時間序列與空間插值方法雖各有優勢，但對於大規模時空資料的填補與預測效能仍缺乏系統性的比較研究。因此，本研究旨在透過嚴謹的模擬實驗，評估不同方法在全球表面溫度資料重建中的適用性與表現差異，為未來氣候監測與決策提供技術參考。

研究方法

本研究以 GES DISC 所提供之 MERRA-2 tavg1_2d_flx_Nx 資料集為基礎，選取經度 $73^\circ \sim 104^\circ$、緯度 $36^\circ \sim 54^\circ$ 的區域作為研究範圍。該區域涵蓋中國西北部、蒙古國西部，以及哈薩克、吉爾吉斯與烏茲別克部分地區，地形多樣，包括高山、高原、盆地、草原、湖泊與河流。多樣化的地理與氣候條件使得不同觀測點的氣溫變化模式存在差異與關聯性，可作為模型多變量輸入的依據。由於本研究聚焦於相對小範圍的時空資料，能在控制資料異質性的同時提高模型預測的準確度。

資料前處理

資料擷取與合併
- 下載 2024 年全年的 MERRA-2 tavg1_2d_flx_Nx 資料。
- 將各小時的觀測資料按時間順序合併，並檢查時間標記的完整性與一致性。
缺值與異常值檢查
- 檢測原始資料中是否存在缺值或異常值。
資料格式處理
- 將全球全年資料重新塑形（reshape）為三維陣列 $(24\ \text{小時},\ 366\ \text{天},\ 207{,}936\ \text{地點})$。
實驗子集選取
- 從全球資料中擷取目標區域（經度 $73^\circ \sim 104^\circ$，緯度 $36^\circ \sim 54^\circ$）。
- 選取 $(24\ \text{小時},\ 260\ \text{天},\ 1{,}850\ \text{地點})$ 作為實驗樣本。
區分已知與未知區域
- 設定固定偽隨機種子碼（seed = 123），以確保實驗可重現性。
- 在目標區域的 1,850 個地點中，隨機選取 1,500 個地點作為已知地點（約佔 81.5%），用於比較預測後的未來趨勢。
- 剩餘 350 個地點作為未知地點（約佔 19.5%），用於比較過往空間填補效果與預測後的未來趨勢。

使用的程式碼如下：

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
# import modules
import os
import numpy as np
import pandas as pd
import xarray as xr
from tqdm import tqdm
import matplotlib.pyplot as plt
import csv
import cartopy.crs as ccrs
import cartopy.feature as cfeature


# functions
def check_data_folder(folder):
    return os.path.exists(folder) and os.path.isdir(folder)

def generate_date_range(start_date, end_date):
    """
    Generate a list of dates from start_date to end_date.
    """
    return pd.date_range(start=start_date, end=end_date, freq='D').strftime('%Y%m%d').tolist()

def load_data(file_path):
    """
    Load data from a NetCDF file.
    """
    if os.path.exists(file_path):
        return xr.open_dataset(file_path)
    else:
        raise FileNotFoundError(f'File not found: {file_path}')


# main program
## check data folder
data_folder = f'..\\..\\..\\surface_air_temperature\\data2024'
if check_data_folder(data_folder):
    print(f'Data folder found: {data_folder}')
else:
    raise FileNotFoundError(f'Data folder not found: {data_folder}')

## load data
start_date = '2024-01-01'
end_date = '2024-12-31'
date_list = generate_date_range(start_date, end_date)

## get location and shape
path = os.path.join(data_folder, f'M2T1NXFLX.5.12.4%3AMERRA2_400.tavg1_2d_flx_Nx.{date_list[0]}.nc4.dap.nc4')
nc4_data = load_data(path)
lat = nc4_data['lat'].values
lon = nc4_data['lon'].values
shape = nc4_data['TLML'].shape
total_locations = shape[1] * shape[2]

## combine data
sample_path = os.path.join(data_folder, f'M2T1NXFLX.5.12.4%3AMERRA2_400.tavg1_2d_flx_Nx.{date_list[0]}.nc4.dap.nc4')
sample_data = load_data(sample_path)
shape_per_file = sample_data['TLML'].shape   # e.g. (24, 361, 576)
time_per_file = len(sample_data['time'])

total_samples = len(date_list)
combined = np.empty((total_samples * shape_per_file[0], *shape_per_file[1:]), dtype=np.float32)
time_list = np.empty(total_samples * time_per_file, dtype=sample_data['time'].dtype)

for i, date in enumerate(tqdm(date_list, desc="Combining")):
    path = os.path.join(data_folder, f'M2T1NXFLX.5.12.4%3AMERRA2_400.tavg1_2d_flx_Nx.{date}.nc4.dap.nc4')
    nc4_data = load_data(path)

    start = i * shape_per_file[0]
    end = (i + 1) * shape_per_file[0]

    combined[start:end] = nc4_data['TLML'].values
    time_list[start:end] = nc4_data['time'].values

print(f'Combined data shape: {combined.shape}')

## reshape data
locations = np.stack(np.meshgrid(lon, lat), axis=-1).reshape(-1, 2)
reshaped_data = combined.reshape(combined.shape[0], -1)
pd.DataFrame(reshaped_data)  # 2d data with time as rows and locations as columns
reshaped_df = pd.DataFrame(reshaped_data, columns=[f"({lon}, {lat})" for lon, lat in locations], index=list(time_list))
time_num = len(time_list)
locations_num = len(locations)

## reshape data to (24, day, location)  (24 hours)
reshaped_df.index = pd.to_datetime(reshaped_df.index)
groups = reshaped_df.groupby(reshaped_df.index.time)
stacked = np.stack([group.to_numpy() for time, group in sorted(groups)])
print(stacked)
time_order = sorted(groups.groups.keys())
print(time_order)

## train set
np.random.seed(123)
valid_mask = (locations[:, 0] >= 73) & (locations[:, 0] <= 104) & \
             (locations[:, 1] >= 36) & (locations[:, 1] <= 54)  # 緯度
valid_indices = np.where(valid_mask)[0]


day_num_train = 250
known_locations_num = valid_indices.shape[0] - 350  # 81.1%
unknown_locations_num = 350                            # 18.9%
locations_num = known_locations_num + unknown_locations_num

locations_index = np.random.choice(valid_indices, size=locations_num, replace=False)
known_locations_index = locations_index[:known_locations_num]
unknown_locations_index = locations_index[known_locations_num:(known_locations_num + unknown_locations_num)]
known_locations_index.sort()
unknown_locations_index.sort()
known_locations_choose = locations[known_locations_index, :]
unknown_locations_choose = locations[unknown_locations_index, :]
stacked_train = stacked[:, :day_num_train, known_locations_index]

future_days = 10
known_real_data = stacked[:, :day_num_train + future_days, known_locations_index]
unknown_real_data = stacked[:, :day_num_train + future_days, unknown_locations_index]

print(f'Load data complete!')

接下來，我們會使用以下不同的模型生成時間序列預測值，並進行如下操作。

時間序列預測
- 使用不同模型，以已知地點的時間序列進行訓練。
- 預測已知地點未來 10 日，也就是 $24 \times 10 = 240$ 個時間步。
- 每次被預測的形狀皆為 $(24, 10, 1850)$ 。
空間補值
- 使用 autoFRK 模型以已知地點的時間序列與地點坐標進行訓練，其中包含以其他模型預測的未來資訊。
- 填補未知地點整條時間序列之值。

`SSSD`

SSSD（Structured State Space Diffusion）是一種結合了條件擴散模型（conditional diffusion model）與結構化狀態空間序列模型（structured state spaces sequence model, S4）的生成式時間序列補值與預測方法。 S4 層能有效捕捉長期依賴性，而 diffusion 機制則提升補值的靈活度與表現力，對於缺值補全與未來趨勢模擬皆擁有優異表現。

https://raw.githubusercontent.com/Josh-test-lab/website-assets-repository/refs/heads/main/posts/1140812%20meeting/updated_architecture.webp — SSSD 模型結構。

在 diffusion 生成機制中，缺值的補值任務將建構為反向擴散過程，逐步從噪音恢復缺失區域，能處理隨機缺值（random missing, RM）、非隨機缺值（non-random missing, NRM），甚至長段黑洞缺值（blackout missing, BM）等複雜情況。

而 S4 則是構化狀態空間模型，利用 HiPPO 理論初始化狀態矩陣，能高效率且穩定地捕捉時間序列中的長期結構關係，為 SSSD 模型提升在長時依賴與跨變量關聯的建模能力。

SSSD 模型能處理多種不同型態的缺值與預測任務，包括：

隨機缺失（random missing, RM）
缺失點隨機分佈於多條時間序列中，缺失位置彼此獨立，常用於模擬因偶發觀測誤差或資料傳輸異常造成的缺值情境。
隨機區塊缺失（random block missing, RBM）
各條時間序列中存在長度與位置隨機的缺失區塊，可視為單一時間序列非隨機缺失的高維版本，常見於儀器短暫失效或局部遮蔽的情況。
黑洞缺失（blockout missing, BM）
與隨機區塊缺失類似，但所有時間序列在相同時間區段出現缺失，屬於同步性缺失模式，多由大範圍系統性故障（如衛星停運或資料伺服器中斷）引起。
時間序列預測（time-series forecasting, TF）
與黑洞缺失類似，但缺失區段位於所有時間序列的末端，對應於未來時段無法觀測的情境，用於模型的未來趨勢預測能力評估。

https://raw.githubusercontent.com/Josh-test-lab/website-assets-repository/refs/heads/main/posts/1140812%20meeting/plots_merged_005.webp — 各類缺失示意圖。

由於 SSSD 需使用三維輸入，因此此處使用調整後的資料集，形狀為 $(24, 250, 1850)$ 。而由先前的經驗，我們知道將 24 小時的資料作為變數輸入，其效果遠比使用 1,850 地點作為變數輸入還差。因此，此處資料集的輸入形狀將調整為 $(1{,}850\ \text{地點},\ 260\ \text{天},\ 24\ \text{小時})$ ，以求更準確的實驗結果。

以下為 SSSD 模型的參數設定。

模型設定 model.yaml
在模型設定中，input_channels 與 output_channels 均設為 24，代表每筆輸入與輸出資料包含 24 維特徵（24 小時逐時氣溫）。 WaveNet 部分採用 32 層殘差層（residual_layers: 32），每層的殘差通道數為 64，並透過 skip connection 將不同層的輸出匯集，以提升多尺度特徵的整合能力。
在 diffusion 模型部分，設定了 200 個擴散步驟（T: 200）， $\beta$ 值由 0.0001 線性增長至 0.02，以控制噪音添加與移除的強度。
而 S4 模型最大序列長度設定為 250（s4_max_sequence_length），與資料集中的時間序列長度一致，代表 250 天；狀態維度為 64（s4_state_dim），並使用雙向計算（s4_bidirectional）與層正規化（s4_use_layer_norm），確保穩定性與泛化能力。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
wavenet:
  # WaveNet model parameters
  input_channels: 24  # Number of input channels
  output_channels: 24  # Number of output channels
  residual_layers: 32  # Number of residual layers
  residual_channels: 64  # Number of channels in residual blocks
  skip_channels: 64  # Number of channels in skip connections

  # Diffusion step embedding dimensions
  diffusion_step_embed_dim_input: 128  # Input dimension
  diffusion_step_embed_dim_hidden: 512  # Middle dimension
  diffusion_step_embed_dim_output: 512  # Output dimension

  # Structured State Spaces sequence model (S4) configurations
  s4_max_sequence_length: 250  # Maximum sequence length
  s4_state_dim: 64  # State dimension
  s4_dropout: 0.0  # Dropout rate
  s4_bidirectional: true  # Whether to use bidirectional layers
  s4_use_layer_norm: true  # Whether to use layer normalization

diffusion:
  # Diffusion model parameters
  T: 200  # Number of diffusion steps
  beta_0: 0.0001  # Initial beta value
  beta_T: 0.02  # Final beta value

模型訓練設定 training.yaml
在訓練設定中，每批資料大小（batch_size）設為 300 。由於訓練資料中有 1,500 筆已知地點，每一迭代需分成 $1500 \div 300 = 5$ 次批次訓練，確保每次迭代涵蓋完整的訓練集。最大迭代數（n_iters）設定為 3,800 ，學習率（learning_rate）為 0.0005。
訓練策略採用只生成缺失值（only_generate_missing: true），並使用時間序列預測缺失遮罩方式（masking: "forecast"），即將序列末端視為未知進行預測。每次遮罩長度（missing_k）為 10。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
# Training configuration
batch_size: 300  # Batch size
output_directory: "./results/surface air temperature/control"  # Output directory for checkpoints and logs
ckpt_iter: "max"  # Checkpoint mode (max or min)
iters_per_ckpt: 100  # Checkpoint frequency (number of epochs)
iters_per_logging: 100  # Log frequency (number of iterations)
n_iters: 3800  # Maximum number of iterations
learning_rate: 0.0005  # Learning rate

# Additional training settings
only_generate_missing: true  # Generate missing values only
use_model: 2  # Model to use for training
masking: "forecast"  # Masking strategy for missing values
missing_k: 10  # Number of missing values

# Data paths
data:
  train_path: "./datasets/surface air temperature/train/train.npy"  # Path to training data

模型填補設定 inference.yaml
推論設定與訓練設定相似，批次大小同樣為 300（batch_size），並從訓練階段的最佳檢查點路徑（ckpt_path）載入模型。推論重複次數（trials）設為 1，每次同樣只針對缺失部分進行生成。缺失遮罩策略與長度與訓練保持一致，以確保一致性。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
# Inference configuration
batch_size: 300  # Batch size for inference
output_directory: "./results/surface air temperature/inference/control"  # Output directory for inference results
ckpt_path: "./results/surface air temperature/control"  # Path to checkpoint for inference
trials: 1 # Replications

# Additional training settings
only_generate_missing: true  # Generate missing values only
use_model: 2  # Model to use for training
masking: "forecast"  # Masking strategy for missing values
missing_k: 10  # Number of missing values

# Data paths
data:
  test_path: "./datasets/surface air temperature/test/test.npy"  # Path to test data

`TSMixer`

與 SSSD 這類深度生成模型不同， TSMixer（Time-Series Mixer）是一種多層感知器（multilayer perceptron, MLP）架構的時間序列預測模型，旨在高效地捕捉時間序列中的時間與特徵維度關聯，並進行多變量預測。 TSMixer 透過堆疊的 MLP 混合層，交替在時間與特徵維度上進行融合，實現時間序列資料的建模。

TSMixer 的核心是多層的混合層，每層首先進行時間維度的混合，然後進行特徵維度的混合，使得模型能夠有效地捕捉時間序列中的時間依賴性和特徵間的關聯。整個模型僅由 MLP 組成，避免了 RNN 或注意力機制中的計算瓶頸，從而提高了訓練效率和預測速度。

此外，TSMixer 支援多種類型的輔助變數，包括過去已知（past covariates）、未來已知（future covariates）及靜態變數（static covariates），使其在多變量及複雜預測任務中具備高度靈活性。在多個長期預測基準測試中，TSMixer 在計算效率與準確度上皆優於傳統 Transformer 類模型。

TSMixer 適用於：

多變量時間序列預測
能夠處理多個變數之間的關聯，適用於如氣象預測、金融市場分析等領域。
長期預測任務
在長期預測中，TSMixer 能夠提供穩定且準確的預測結果。
資源受限的環境
由於其高效的計算性能，TSMixer 適用於計算資源有限的場景。

以下使用基於 PyTorch 的 darts 模組中的 TSMixerModel，逐日（24小時）訓練並預測多變量時間序列。

在此程式碼中，input_chunk_length=30 指模型觀察過去 30 天的資料，output_chunk_length=10 則為預測未來 10 天。同時訓練迭代達 3,800 次。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
from darts.models import TSMixerModel
from darts import TimeSeries
from datetime import datetime
from tqdm import tqdm
real_data = known_real_data
location_choose = known_locations_choose
locations_index = known_locations_index
inference = np.array([[[0.0] * real_data.shape[2]] * future_days] * real_data.shape[0])
start_time = datetime.now()

for i in tqdm(range(real_data.shape[0])):
    train_set = real_data[i, :day_num_train, :]
    train_set = pd.DataFrame(train_set, index=pd.to_datetime(date_list[:day_num_train]), columns=[f"loc_{i}" for i in range(known_locations_num)])
    train_set = TimeSeries.from_dataframe(train_set)

    test_set = real_data[i, day_num_train:, :]
    test_set = pd.DataFrame(test_set, index=pd.to_datetime(date_list[day_num_train:day_num_train + future_days]), columns=[f"loc_{i}" for i in range(known_locations_num)])
    test_set = TimeSeries.from_dataframe(test_set)

    model = TSMixerModel(
        input_chunk_length=30,
        output_chunk_length=10,
        n_epochs=3800,
        dropout=0.0005,
        use_reversible_instance_norm=True,
        random_state=42,
        pl_trainer_kwargs={"accelerator": "gpu"}
    )

    model.fit(train_set)

    forecast = model.predict(future_days)
    inference[i] = forecast.values()

print(f'Inference complete! Time taken: {datetime.now() - start_time}')

`RegressionEnsemble`

RegressionEnsemble 是 darts 模組中的一種集成預測模型，利用迴歸模型（如線性迴歸）融合多個基礎預測模型的輸出，以提升整體預測的準確度。此模型採用「堆疊（stacking）」技術，將多個基礎預測模型的預測結果作為特徵，再訓練一個迴歸模型來學習最佳融合權重。

此模型可同時預測多個變數，並能整合多種輔助資訊，包括過去時間點已知的變數、未來時間點已知的變數，以及靜態不變的變數，提升預測準確度。 RegressionEnsemble 適用於：

多變量時間序列預測
當需要同時預測多個變數時，該模型能夠有效地融合不同模型的預測結果。
長期預測任務
集成多模型結果提升預測穩定度與精準度。
資源受限的環境
由於迴歸融合（Regression Ensemble）具備計算效率高的優點，因此特別適合在計算資源有限的環境中使用。

以下使用基於 PyTorch 的 darts 模組中的 RegressionEnsembleModel，逐日（24小時）訓練並預測多變量時間序列。

在此程式碼中，我們使用了三個基礎預測模型來構建迴歸融合模型。首先，NaiveSeasonal(K=7) 是一個基於週期性假設的簡單季節性模型，假設時間序列每隔 7 天會呈現相似的模式。其次，LinearRegressionModel(lags=30) 是一個線性迴歸模型，利用過去 30 個時間點的資料作為特徵，來捕捉時間序列中的趨勢與變化。最後，NaiveDrift() 模型基於序列的整體漂移趨勢，進行簡單的趨勢預測。

這三個基礎模型的預測結果將作為輸入特徵，通過 RegressionEnsembleModel 中的迴歸融合器來學習如何最佳組合這些預測，提高整體的預測精度。其中，參數 regression_train_n_points=30 指定融合器使用過去30個時間點的預測結果來訓練迴歸模型，使模型能更有效地捕捉時間序列的變化規律並生成準確的預測結果。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
from darts.models import RegressionEnsembleModel, NaiveSeasonal, LinearRegressionModel, NaiveDrift
from darts import TimeSeries
from datetime import datetime
from tqdm import tqdm
real_data = known_real_data
location_choose = known_locations_choose
locations_index = known_locations_index
inference = np.array([[[0.0] * real_data.shape[2]] * future_days] * real_data.shape[0])
start_time = datetime.now()

for i in tqdm(range(real_data.shape[0])):
    train_set = real_data[i, :day_num_train, :]
    train_set = pd.DataFrame(train_set, index=pd.to_datetime(date_list[:day_num_train]), columns=[f"loc_{i}" for i in range(known_locations_num)])
    train_set = TimeSeries.from_dataframe(train_set)

    test_set = real_data[i, day_num_train:, :]
    test_set = pd.DataFrame(test_set, index=pd.to_datetime(date_list[day_num_train:day_num_train + future_days]), columns=[f"loc_{i}" for i in range(known_locations_num)])
    test_set = TimeSeries.from_dataframe(test_set)

    base_models = [
        NaiveSeasonal(K=7),
        LinearRegressionModel(lags=30),
        NaiveDrift()
    ]

    model = RegressionEnsembleModel(
        forecasting_models=base_models,
        regression_train_n_points=30
    )

    model.fit(train_set)

    forecast = model.predict(n=future_days)
    inference[i] = forecast.values()

print(f'Inference complete! Time taken: {datetime.now() - start_time}')

`autoFRK`

autoFRK（Automatic Fixed Rank Kriging）是一種基於空間統計理論的高效空間插值與預測方法。它結合了固定秩克里金（fixed rank kriging, FRK）的維度約簡技術，利用多尺度基底函數來捕捉空間數據中的不同變異特徵，並自動化選擇模型參數。

數學上， autoFRK 的模型形式可表示為：

$$ z[t] = \mu + G \cdot w[t] + \eta[t] + e[t], \quad w[t] \sim N(0, M), \quad e[t] \sim N(0, s \cdot D); \quad t = 1, \cdots, T, $$

其中， $z[t]$ 是在 $n$ 個地點觀測到的（部分）資料向量； $\mu$是長度為 $n$ 的常數均值向量； $D$ 是已知的 $n \times n$ 矩陣； $G$ 是已知的 $n \times K$ 矩陣； $\eta[t]$ 是長度為 $n$ 的隨機向量，對應於空間平穩過程； $w[t]$ 是長度為 $K$ 的不可觀測隨機權重向量。

參數透過最大概似（maximum likelihood）以封閉形式（closed-form expression）估計。基底函數矩陣 $G$ 採用有序薄板樣條函數（thin-plate spline functions），基底數量則由赤池資訊準則（Akaike’s information criterion, AIC）選擇。

autoFRK 利用基底函數 $G$ 展開來表示空間隨機場，將高維度的空間資料壓縮至低維度的係數空間，大幅降低計算複雜度，特別適用於大規模的空間資料分析。透過最大化似然函數來估計模型參數，autoFRK 可有效捕捉空間結構並提供準確的空間預測與不確定性評估。

在本研究中，autoFRK 用於將基於時間序列預測模型（如 SSSD、TSMixer、RegressionEnsemble）所產生的時間維度預測結果，進行空間維度的插值與補全。這種時空結合的方式，有助於彌補資料中的空缺，並提升全球表面溫度監測的完整性與精確度。

以下為 autoFRK 在本研究中的應用流程：

輸入時間序列模型產生的逐時間點預測結果。
透過 autoFRK 進行空間基底函數展開與參數估計。
對未知地點進行空間插值與預測，補足資料缺失區域。
結合時間與空間預測結果，提供完整的全球表面溫度時空場。

用於空間填補的程式碼如下：

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
# load data
library(reticulate)
np <- import("numpy")

data_path = "../train_frk.npy"
data = np$load(data_path)

real_path = "../real_data.npy"
real_data = np$load(real_path)
known_data = real_data[ , , 1:1500]
unknown_data = real_data[ , , 1501:1850]

known_loc_path = "../known_location.npy"
known_locs = np$load(known_loc_path)

unknown_loc_path = "../unknown_location.npy"
unknown_locs = np$load(unknown_loc_path)

# autoFRK
library(autoFRK)
library(dplyr)

start_time <- Sys.time()

unknown_inference_shape = c(dim(data)[1:2],  dim(known_locs)[1] + dim(unknown_locs)[1])
unknown_inference <- array(NA, dim = unknown_inference_shape)
mrts_basis <- NULL

locs = rbind(known_locs, unknown_locs)
total_iter <- 24 * dim(data)[2]
pb <- txtProgressBar(min = 0, max = total_iter, style = 3)
iter <- 0

for (time_ in 1:24) {
  hour_data <- data[time_, , ]
  
  for (day_ in 1:dim(hour_data)[1]) {
    iter <- iter + 1
    iter_start <- Sys.time()
    
    day_data = hour_data[day_, ]
    
    if (is.null(mrts_basis)) {
      model = autoFRK(data = day_data, loc = known_locs)
      mrts_basis = model$G
    } else {
      model = autoFRK(data = day_data, loc = known_locs, G = mrts_basis)
    }
    
    pred = predict.FRK(object = model, newloc = locs)
    unknown_inference[time_, day_, ] = pred$pred.value
    
    elapsed <- as.numeric(difftime(Sys.time(), start_time, units = "secs"))
    avg_time <- elapsed / iter
    remaining <- avg_time * (total_iter - iter)
    est_done <- format(Sys.time() + remaining, "%H:%M:%S")
    
    setTxtProgressBar(pb, iter)
    cat(sprintf(" | Est. done at: %s | Remaining: %ds\r", est_done, round(remaining)))
  }
}

close(pb)

end_time <- Sys.time()
cat("\n")
cat("Total time elapsed:", format(end_time - start_time))
cat("\n")

# result
y_inf <- unknown_inference
real <- real_data
train <- known_data
test <- unknown_data

mspe <- function(pred, true) mean((pred - true)^2)
rmspe <- function(pred, true) sqrt(mean((pred - true)^2))
mape <- function(pred, true) mean(abs(pred - true))

mspe_p <- function(pred, true) mean((pred - true)^2 / true)
rmspe_p <- function(pred, true) sqrt(mean((pred - true)^2 / true))
mape_p <- function(pred, true) mean(abs(pred - true) / true)

future_days = 10
future_idx <- (dim(y_inf)[2] - future_days + 1):dim(y_inf)[2]
past_idx <- 1:(dim(y_inf)[2] - future_days)
known_idx <- 1:dim(known_locs)[1]
unknown_idx <- (dim(known_locs)[1] + 1):(dim(known_locs)[1] + dim(unknown_locs)[1])

compute_metrics <- function(pred, true) {
  c(
    MSPE = mspe(pred, true),
    RMSPE = rmspe(pred, true),
    `MSPE%` = mspe_p(pred, true),
    `RMSPE%` = rmspe_p(pred, true),
    MAPE = mape(pred, true),
    `MAPE%` = mape_p(pred, true)
  )
}

result_table <- data.frame(
  row.names = c('MSPE', 'RMSPE', 'MSPE%', 'RMSPE%', 'MAPE', 'MAPE%'),

  `ALL Locs & All Time` = compute_metrics(y_inf, real),
  `Known Locs & All Time` = compute_metrics(y_inf[, , known_idx], train),
  `Unknown Locs & All Time` = compute_metrics(y_inf[, , unknown_idx], test),

  `ALL Locs & Future` = compute_metrics(y_inf[, future_idx, ], real[, future_idx, ]),
  `Known Locs & Future` = compute_metrics(y_inf[, future_idx, known_idx], train[, future_idx, ]),
  `Unknown Locs & Future` = compute_metrics(y_inf[, future_idx, unknown_idx], test[, future_idx, ]),

  `ALL Locs & Past` = compute_metrics(y_inf[, past_idx, ], real[, past_idx, ]),
  `Known Locs & Past` = compute_metrics(y_inf[, past_idx, known_idx], train[, past_idx, ]),
  `Unknown Locs & Past` = compute_metrics(y_inf[, past_idx, unknown_idx], test[, past_idx, ])
)

print(result_table)


# save to .npy

# combined data
output_matrix = unknown_inference
output_matrix %>% dim() %>% print()

# save to .npy
output_matrix = output_matrix %>% r_to_py()
save_path = "../plot_frk.npy"
np$save(save_path, output_matrix)

實驗結果

上述實驗除 RegressionEnsemble 模型外，其迭代次數皆設為 3,800 次，且三種模型之資料格式皆相同，惟 TSMixer 與 RegressionEnsemble 模型僅接受二維資料及作為輸入，故使用迴圈方式進行訓練與預測。

結論

指標 / 模型	SSSD + autoFRK	TSMixer + autoFRK	RegressionEnsemble + autoFRK
MSPE (ALL Locs & Future)	24.97269111	75.99501027	33.41590875
MSPE (Known Locs & Future)	24.97548082	76.24826754	33.23029709
MSPE (Unknown Locs & Future)	24.96073520	74.90962190	34.21138730

RMSPE (ALL Locs & Future)	4.99726836	8.71751170	5.78064951
RMSPE (Known Locs & Future)	4.99754748	8.73202540	5.76457259
RMSPE (Unknown Locs & Future)	4.99607198	8.65503450	5.84905012

MSPE% (ALL Locs & Future)	0.08875877	0.26929858	0.11762572
MSPE% (Known Locs & Future)	0.08882456	0.27033777	0.11704136
MSPE% (Unknown Locs & Future)	0.08847681	0.26484490	0.12013010

RMSPE% (ALL Locs & Future)	0.29792410	0.51893986	0.34296606
RMSPE% (Known Locs & Future)	0.29803450	0.51994016	0.34211308
RMSPE% (Unknown Locs & Future)	0.29745052	0.51463080	0.34659790

MAPE (ALL Locs & Future)	4.06021896	7.58070981	4.61292859
MAPE (Known Locs & Future)	4.05973676	7.59905546	4.59836726
MAPE (Unknown Locs & Future)	4.06228551	7.50208560	4.67533431

MAPE% (ALL Locs & Future)	0.01437071	0.02678231	0.01624002
MAPE% (Known Locs & Future)	0.01437764	0.02686171	0.01619796
MAPE% (Unknown Locs & Future)	0.01434101	0.02644200	0.01642026

由上述實驗可以發現，SSSD + autoFRK 在時間序列預測（Future）表現相對穩定，尤其是在未知區域時間序列預測（Unknown Locs & Future）的各項指標（MSPE、RMSPE、MAPE）數值明顯低於 TSMixer + autoFRK，且與已知區域（Known Locs）的差距極小，顯示其在未觀測地點的泛化能力較佳。然而，SSSD 模型的時間成本較高，對於需要即時預測的場景可能不太適合。若能在不顯著增加運算負擔的情況下提升迭代次數，或許能進一步增強其預測能力。

相較之下，TSMixer + autoFRK 在時間序列預測的 MSPE 與 RMSPE 明顯高於其他方法（例如 MSPE 高達 75.99），特別是在未知區域的誤差更大，顯示其在未觀測地點的預測偏差較為嚴重。此模型可能需要針對泛化能力進行調整，例如增加迭代次數、改進訓練策略或引入額外特徵，才能縮小與其他方法的差距。

RegressionEnsemble + autoFRK 的整體表現介於 SSSD 與 TSMixer 之間，在未知區域的 MSPE 與 RMSPE 雖然高於 SSSD，但遠低於 TSMixer，泛化能力屬中等偏佳。此外，RegressionEnsemble 的推論時間明顯低於 SSSD，若能在不顯著增加運算成本的情況下引入更多高效的基礎預測模型，則有潛力超越 SSSD 的預測效果。

結語

https://raw.githubusercontent.com/Josh-test-lab/website-assets-repository/refs/heads/main/posts/1140812%20meeting/To%20be%20continued.jpg — To be continued!

運行環境

本機作業系統：Windows 11 24H2
- 程式語言：Python 3.12.9
計算平臺：財團法人國家實驗研究院國家高速網路與計算中心臺灣 AI 雲
- 作業系統：Ubuntu
- Miniconda
- GPU：NVIDIA Tesla V100 32GB GPU
- CUDA 12.8 driver
- 程式語言：Python 3.10.16 for Linux

延伸學習

我測試此項目的 Github 儲存庫。

參考資料

Global Modeling and Assimilation Office (GMAO)。(2015)。MERRA-2 tavg1_2d_flx_Nx: 2d,1-Hourly,Time-Averaged,Single-Level,Assimilation,Surface Flux Diagnostics (Version 5.12.4) [資料集]。Goddard Earth Sciences Data and Information Services Center (GES DISC)。參考自 https://doi.org/10.5067/7MCPBJ41Y0K6
Juan Lopez Alcaraz 、 Nils Strodthoff（2022）。Diffusion-based time series imputation and forecasting with structured state space models。Transactions on Machine Learning Research。參考自 https://openreview.net/forum?id=hHiIbk7ApW
SSSD（2022）。GitHub。參考自 https://github.com/AI4HealthUOL/SSSD
SSSD_CP（2024）。GitHub。參考自 https://github.com/egpivo/SSSD_CP
Unit8 SA（無日期）。Time Series Made Easy in Python。Darts。參考自 https://unit8co.github.io/darts/index.html
darts（2025）。GitHub。參考自 https://github.com/unit8co/darts
Tzeng, S., & Huang, H. C. (2018). Resolution Adaptive Fixed Rank Kriging. Technometrics, 60(2), 198–208. 參考自 https://doi.org/10.1080/00401706.2017.1345701
autoFRK（2024）。GitHub。參考自 https://github.com/egpivo/autoFRK
Si-An Chen, Chun-Liang Li, Nate Yoder, Sercan O. Arik and Tomas Pfister. (2023). TSMixer: An all-MLP architecture for time series forecasting. arXiv. 參考自 https://arxiv.org/abs/2303.06053

目錄

1140812 meeting

前言

研究動機

研究方法

資料前處理

SSSD

TSMixer

RegressionEnsemble

autoFRK

實驗結果

SSSD + autoFRK

TSMixer + autoFRK

RegressionEnsemble + autoFRK

結論

結語

運行環境

延伸學習

參考資料

相關文章

`SSSD`

`TSMixer`

`RegressionEnsemble`

`autoFRK`

`SSSD` + `autoFRK`

`TSMixer` + `autoFRK`

`RegressionEnsemble` + `autoFRK`