目錄

基於結構化狀態空間擴散模型與自適應固定秩克利金法之未知區域時空預測

碩士論文

摘要

隨著感測與遙測技術迅速發展,交通、氣象與環境監測資料呈現高頻且跨區域的時空特性。然而,實務資料常伴隨觀測缺漏、採樣不均與未觀測地點推估等問題,使時空資訊之重建與預測面臨挑戰。本研究提出時空預測框架 $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK,結合基於 S4 層之結構化狀態空間擴散模型(Structured State Space Diffusion with S4 Layers, $\text{SSSD}^\text{S4}$)與多解析度薄板樣條基底函數(Multi-Resolution Thin-Plate Spline Basis Functions, MRTS),以同時處理時間序列之長期依賴關係與空間結構估計。模型將 MRTS 作為額外空間條件資訊引入 S4 層,使 $\text{SSSD}^{\text{S4+MRTS}}$ 能夠萃取已觀測地點為條件之時空動態特徵,並於推論階段結合以 MRTS 為空間基底之自適應固定秩克利金(Resolution Adaptive Fixed Rank Kriging, AFRK),在維持預測場域空間一致性的同時,對未觀測位置進行插值。實驗分別以 Weather2K 、 MERRA-2 與第二屆大型資料集空間統計競賽資料集進行驗證,並與 Temporal Fusion Transformers (TFT) 、 Vector Autoregression (VAR) 、 Stochastic Variational Gaussian Process (SVGP) 及 Spatio-temporal DeepKriging (STDK) 等基準方法比較。結果顯示,$\text{SSSD}^{\text{S4+MRTS}}$ 在多數實驗設定下均能改善平均平方預測誤差(MSPE),尤其於未觀測地點之未來預測表現最為顯著。綜上所述,本研究展示時間動態模型與空間統計方法之整合,可有效提升缺失資料與空間外插情境下之時空預測能力。

關鍵字:結構化狀態空間擴散模型、多解析度薄板樣條基底函數、自適應固定秩克利金、時空預測、空間插值、缺失資料補值

1 緒論

隨著感測技術與資料基礎建設的快速發展,交通監測、水文分析、氣象觀測及衛星遙測等領域已能在高時間頻率、廣泛空間範圍以及長時間尺度下蒐集大規模資料,因而形成具有豐富時空結構的資料集 (Gneiting, 2002; N. Cressie and Wikle, 2011) 。然而,真實世界的觀測資料經常伴隨觀測缺漏、採樣頻率不一致、測站分布稀疏以及部分區域長期缺乏觀測等問題,使得資料呈現高度不完整性 (Little and Rubin, 2002; Decorte et al., 2024) 。此類缺失不僅降低統計推論的可靠性,也使得模型在未來預測與跨空間推估時的穩健性與泛化能力下降。

在實際應用中,資料缺失可能源自多種機制,包括隨機缺失(Missing at Random, MAR)以及因長期空間空缺或整段時間缺乏觀測所造成的結構性缺失。當研究目標由單純的缺失補值進一步擴展至未觀測地點未來狀態的預測時,問題本質上同時涉及時間外插(temporal extrapolation)與空間外插(spatial extrapolation)。前者要求模型能有效捕捉非線性且具長期依賴性的時間動態,後者則要求模型能維持空間連續性與一致的空間結構。現有方法多半僅能處理其中一個面向。傳統時間序列模型如自迴歸整合移動平均模型(Autoregressive Integrated Moving Average, ARIMA)與向量自迴歸模型(Vector Autoregression, VAR)雖能描述時間動態,但缺乏顯式的空間結構;空間統計方法如克利金法(Kriging)與固定秩克利金法(Fixed Rank Kriging, FRK)能刻畫空間依存關係,卻不適合長序列時間建模;深度學習模型如長短期記憶網路(Long Short-Term Memory, LSTM)與時序融合變換器(Temporal Fusion Transformer, TFT)雖具有良好的非線性表徵能力,但在長序列、高缺失率或跨空間預測情境下,仍面臨穩定性不足與外插能力受限等問題 (Lim et al., 2020; Wu et al., 2022) 。

近年來,狀態空間模型(State Space Models, SSMs)獲得廣泛關注,其中以結構化狀態空間序列模型(Structured State Space Sequence Model, S4)及其延伸之結構化狀態空間擴散模型(Structured State Space Diffusion, SSSD)(Gu, Dao, et al., 2020; Gu, Goel, and Ré, 2022; Alcaraz and Strodthoff, 2023) 為代表,在長序列建模與缺失資料補全方面展現優異表現。然而,此類模型本質上仍屬於時間導向架構,缺乏顯式空間結構,因此在空間外插時較難維持地理連續性。另一方面,自適應固定秩克利金法(Resolution Adaptive Fixed Rank Kriging, AFRK)(N. Cressie and Johannesson, 2008; Tzeng and Huang, 2018) 提供具計算效率的空間低秩近似方法,能有效描述大尺度空間結構,但其設計目的並非用於建模長期時間動態或複雜非線性序列。

上述現象指出現有方法難以同時兼顧長序列時間動態建模、空間結構一致性與高度缺失資料之處理能力。為回應此一問題,本研究提出一個整合結構化狀態空間擴散模型(SSSD)與自適應固定秩克利金(AFRK)之統一時空預測框架。此架構的核心概念在於利用 SSSD 的長序列建模能力克服高缺失率與強時間依賴情境下的預測困難,同時透過 AFRK 的空間低秩結構補足純時間模型缺乏空間連續性的限制。藉由結合以多解析度薄板樣條基底函數(Multi-Resolution Thin-Plate Spline Basis Functions, MRTS)所建構之空間表示方法與擴散式時間建模機制,所提出之框架能同時捕捉長期時間動態、空間結構資訊以及不完整觀測資料的特性,進而提供一套適用於不完整觀測情境下之空間預測與長期時間預測的統一時空分析框架。

2 相關研究

時空預測方法通常可歸納為時間序列模型、空間統計模型,以及整合時間與空間資訊之時空模型。雖然各類方法皆已獲得相當程度之發展,然而如何同時兼顧時間動態特徵與空間相依結構,仍為時空預測領域的重要研究課題。以下將分別介紹各類代表性方法,並探討其優勢、限制及適用情境。

2.1 時間序列模型

2.1.1 向量自迴歸模型

向量自迴歸模型(Vector Autoregression model, VAR)係由 Sims (1980) 提出,為單變量自迴歸模型(AR)之擴展,旨在捕捉多個隨時間變化變數間的動態相互依賴關係。與單變量模型不同,VAR 將系統中所有變數均視為內生變數(endogenous variables),並透過聯立方程式組描述變數間的滯後影響。

對於一個包含 $k$ 個變數之 $p$ 階向量自迴歸模型 $\text{VAR}(p)$,其數學表達式為: \begin{align} \boldsymbol{y}_t = \boldsymbol{c} + \sum_{i=1}^{p} \boldsymbol{\Phi}_i \boldsymbol{y}_{t-i} + \boldsymbol{\epsilon}_t, \end{align} 其中, $\boldsymbol{y}_t \in \mathbb{R}^k$ 為 $t$ 時刻之觀測向量; $\boldsymbol{c}$ 為截距項向量; $\boldsymbol{\Phi}_i \in \mathbb{R}^{k \times k}$ 為滯後算子矩陣,用以衡量 $t-i$ 時刻之變數對當前狀態的影響強度; $\boldsymbol{\epsilon}_t \sim \mathcal{N}(0, \boldsymbol{\Omega})$ 為雜訊向量。

VAR 模型之優勢在於其能透過脈衝響應函數(Impulse Response Function)與變異數分解(Variance Decomposition)分析變數間的動態交互影響 (Sims, 1980) 。然而,隨著變數數量 $k$ 或滯後階數 $p$ 增加,參數數量將呈平方增長,易導致過度擬合問題,且傳統線性 VAR 模型難以捕捉複雜的非線性結構與參數的時變性 (Primiceri, 2005) 。

2.1.2 時序融合變換器

鑒於傳統線性多變量模型在捕捉複雜非線性動態與高維特徵關聯上之侷限,近期研究多轉向利用深度學習架構以提升預測效能。為克服傳統循環架構在處理長距離依賴時之困境,並有效整合包含靜態背景資訊、已知未來變數與歷史觀測值在內的多樣化資訊來源, Lim et al. (2020) 提出了時序融合變換器(Temporal Fusion Transformers, TFT)。該模型係一種專為多步時序預測設計的深度學習架構,其設計重點在於透過專門的網路組件對不同性質的輸入變數進行權重分配,從而提升預測結果的解釋能力,補足傳統統計模型在處理大規模異質資料時之不足。

TFT 之實現原理依賴於門控殘差網路(Gated Residual Networks, GRN),該組件透過門控線性單元(Gated Linear Units, GLU)控制資訊流動,使模型能根據資料特性自動調整非線性轉換之深度。對於輸入向量 $\boldsymbol{a}$ 與選用的背景向量 $\boldsymbol{c}$,其運算過程如下: \begin{align} \text{GRN}_\omega(\boldsymbol{a}, \boldsymbol{c}) = \text{LayerNorm}(\boldsymbol{a} + \text{GLU}_\omega(\boldsymbol{\eta}_1)), \end{align} 其中 $\boldsymbol{\eta}_1$ 為經過權重矩陣轉換後之特徵向量。此機制確保了模型在處理不同複雜度的序列時具備高度彈性,有效避免了深層網路在簡單資料集上的過度擬合問題。

針對包含大量外部因子之時空資料, TFT 導入了變數選擇網路(Variable Selection Networks, VSN),旨在從眾多輸入特徵中篩選關鍵變數。透過為每個特徵分配權重 $\nu_t^{(i)}$ ,模型得以自動忽略冗餘資訊並聚焦於具影響力之因子,其整合後之特徵向量表示如下: \begin{align} \tilde{\boldsymbol{\xi}}_t = \sum_{i=1}^{m} \nu_t^{(i)} \tilde{\boldsymbol{\xi}}_t^{(i)}, \end{align} 其中 $\tilde{\boldsymbol{\xi}}_t^{(i)}$ 為處理後之特徵向量。此設計大幅提升了模型在面對多維度特徵輸入時的穩健性,並使研究者能直觀地量化各類變數對預測結果之貢獻。

在時間序列關係的捕捉上, TFT 採用改良之時序自我注意力機制(Temporal Self-Attention)來處理長期依賴關係。相較於標準 Transformer 架構, TFT 在注意力層中加入門控層進行殘差連接,並透過解碼器整合歷史與未來的時空背景資訊,進而找出對當前預測最具影響力的時間點。 TFT 不僅在預測精確度上表現優異,更賦予了深度學習模型在特定時間步或特定特徵上之解釋能力 (Lim et al., 2020) 。

2.1.3 狀態空間模型

狀態空間模型(State Space Model, SSM)是一類透過潛在狀態(latent state)向量來描述動態系統或序列資料的數學模型。其最初由 Kalman (1960) 在控制理論與濾波領域提出,用以解決線性動態系統的最佳濾波與預測問題。隨後, Gu, Goel, and Ré (2022) 等人將該概念推廣至深度學習架構中,用以處理長序列的時間序列建模,證明 SSM 相較於傳統 RNN 與 LSTM 在捕捉長期依賴性(long-range dependencies)與維持穩定梯度方面表現更優 (Gu, Goel, and Ré, 2022) 。

給定一維輸入訊號序列 $\boldsymbol{u}(t)$ 與一維輸出訊號序列 $\boldsymbol{y}(t)$, SSM 的基本形式為 \begin{align} \begin{aligned} \boldsymbol{x}’(t) & = \boldsymbol{A} \boldsymbol{x}(t) + \boldsymbol{B} \boldsymbol{u}(t); \\ \boldsymbol{y}(t) & = \boldsymbol{C} \boldsymbol{x}(t) + \boldsymbol{D} \boldsymbol{u}(t), \end{aligned} \end{align} 其中, $\boldsymbol{x}(t) \in \mathbb{R}^N$ 為 $\boldsymbol{u}(t)$ 映射的 $N$ 維度潛在狀態; $\boldsymbol{x}’(t) = \frac{d}{dt}\boldsymbol{x}(t)$ ; $\boldsymbol{A} \in \mathbb{R}^{N \times N}$ 為狀態矩陣(state matrix); $\boldsymbol{B} \in \mathbb{R}^{N \times 1}$ 與 $\boldsymbol{C} \in \mathbb{R}^{1 \times N}$ 分別為輸入與輸出矩陣,用於描述輸入對狀態的影響以及狀態對輸出的映射; 而 $\boldsymbol{D} \in \mathbb{R}$ 則為前饋矩陣(feedthrough matrix),使輸入向量直接影響輸出結果,通常為零矩陣。在深度學習中,這些參數通常透過梯度下降進行學習。

https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/fig-SSM.svg
圖 1: 典型的狀態空間模型。

為解決 SSM 在實務中,梯度可能會出現隨序列長度呈指數增減的問題, Gu, Goel, and Ré (2022) 引入 HiPPO (High-order Polynomial Projection Operators)矩陣 (Gu, Dao, et al., 2020) ,用以取代原 (4) 中的隨機矩陣 $\boldsymbol{A}$ 。

\begin{align} \boldsymbol{A}_{nk} = - \begin{cases} (2n + 1)^{1/2} (2k + 1)^{1/2}, & \text{if } n > k; \\ n + 1, & \text{if } n = k; \\ 0, & \text{if } n < k. \end{cases} \end{align}

(5) 所示的 HiPPO 矩陣能夠隨時間變化衡量過去每個時間步的重要性,動態調整記憶並記住所有歷史。將原本的隨機矩陣 $\boldsymbol{A}$ 替換為 HiPPO 矩陣後,狀態 $\boldsymbol{x}(t)$ 可以有效記憶輸入序列 $\boldsymbol{u}(t)$ 的歷史資訊,同時避免梯度爆炸或消失。實驗結果顯示,這種設計不僅能提升模型的計算穩定性,也顯著改善長序列預測的性能 (Gu, Goel, and Ré, 2022) 。

2.1.4 結構化狀態空間序列模型

基於前一節的理論, Gu, Goel, and Ré (2022) 提出了結構化狀態空間序列模型(Structured State Space Sequence Model, S4)。該模型旨在將連續時間狀態空間模型離散化,並應用於深度學習框架中以處理長序列資料。

S4 以 SSM 為基礎,將 (4) 離散化,以便將 SSM 應用於離散的輸入序列。定義步長為 $\Delta$ ,則離散化的 SSM 如下: \begin{align} \begin{aligned} \boldsymbol{x}_k & = \overline{\boldsymbol{A}} \boldsymbol{x}_{k - 1} + \overline{\boldsymbol{B}} \boldsymbol{u}_k; \\ \boldsymbol{y}_k & = \overline{\boldsymbol{C}} \boldsymbol{x}_k, \end{aligned} \end{align} 其中, $\overline{\boldsymbol{A}}, \overline{\boldsymbol{B}}, \overline{\boldsymbol{C}}$ 為離散近似形式的 $\boldsymbol{A}, \boldsymbol{B}, \boldsymbol{C}$ 。 \begin{align} \begin{aligned} \overline{\boldsymbol{A}} & = (\boldsymbol{I} - \Delta / 2 \cdot \boldsymbol{A})^{-1} (\boldsymbol{I} + \Delta / 2 \cdot \boldsymbol{A}); \\ \overline{\boldsymbol{B}} & = (\boldsymbol{I} - \Delta / 2 \cdot \boldsymbol{A})^{-1} \Delta \boldsymbol{B}; \\ \overline{\boldsymbol{C}} & = \boldsymbol{C}. \end{aligned} \end{align}

為降低矩陣運算成本, S4 可透過對角化將離散矩陣,表示為不同基底下的等價形式。 \begin{align} \begin{aligned} (\boldsymbol{A}, \boldsymbol{B}, \boldsymbol{C}) \sim (\boldsymbol{V}^{-1} \boldsymbol{AV}, \boldsymbol{V}^{-1} \boldsymbol{B}, \boldsymbol{CV}), \end{aligned} \end{align} 其中 $\boldsymbol{V}$ 為基底變換矩陣。同時,離散化 SSM 可展開為卷積形式,以提高平行運算效率: \begin{align} \boldsymbol{y} & = \overline{\boldsymbol{K}} * \boldsymbol{u}; \\ \overline{\boldsymbol{K}} & \in \mathbb{R}^L := (\overline{\boldsymbol{C}}\overline{\boldsymbol{B}}, \overline{\boldsymbol{C}}\overline{\boldsymbol{A}}\overline{\boldsymbol{B}}, \cdots, \overline{\boldsymbol{C}}\overline{\boldsymbol{A}}^{L - 1}\overline{\boldsymbol{B}}), \end{align} 其中, $\overline{\boldsymbol{K}}$ 為 SSM 卷積核, $L$ 表示卷積長度。

為進一步降低運算複雜度, S4 透過 Normal Plus Low-Rank (NPLR) 方式,將 $\overline{\boldsymbol{A}}$ 矩陣參數化: \begin{align} \begin{aligned} \boldsymbol{A} = \boldsymbol{V \Lambda V}^* - \boldsymbol{PQ}^\top = \boldsymbol{V} (\boldsymbol{\Lambda} - (\boldsymbol{V}^*\boldsymbol{P})(\boldsymbol{V}^*\boldsymbol{Q})^*) \boldsymbol{V}^*, \end{aligned} \end{align} 其中 $\boldsymbol{\Lambda}$ 為對角矩陣,$\boldsymbol{P}, \boldsymbol{Q} \in \mathbb{R}^{N \times r}$ 為低秩(low-rank)矩陣,$\boldsymbol{V} \in \mathbb{C}^{N \times N}$ 為酉(unitary)矩陣。

綜合以上研究, S4 結合 HiPPO 矩陣、離散化、對角化、卷積與 NPLR 參數化,使模型在長序列任務中實現高效運算與優異表現 (Gu, Goel, and Ré, 2022) 。

2.1.5 擴散模型

擴散模型(Diffusion Model)是透過建立正向擴散過程(forward diffusion process)與反向去噪過程(reverse denoising process)的對偶結構,學習資料的生成分佈 (Sohl-Dickstein et al., 2015) 。其模型逐步向原始資料注入高斯雜訊,使其最終接近標準常態分佈,並訓練模型學習該過程的反向映射以重建資料分佈。近年來,該方法被應用於時間序列插補,僅對缺失區段進行擴散與去噪操作,以在條件觀測下恢復完整序列 (Alcaraz and Strodthoff, 2023) 。

令 $\boldsymbol{x}_0 \sim q(\boldsymbol{x}_0)$ 表示原始資料樣本,正向過程被定義為一個固定參數的高斯馬可夫鏈,用以模擬逐步擾動的資料生成機制: \begin{align} \begin{cases} q(\boldsymbol{x}_{1} | \boldsymbol{x}_0) = \prod_{t=1}^T q(\boldsymbol{x}_t | \boldsymbol{x}_{t-1}); \\ q(\boldsymbol{x}_t | \boldsymbol{x}_{t-1}) = \mathcal{N}(\boldsymbol{x}_t; \sqrt{1 - \beta_t} \boldsymbol{x}_{t-1}, \beta_t \boldsymbol{I}), \end{cases} \end{align} 其中, $\beta_t$ 為控制雜訊強度的變異數; $\mathcal{N}$ 為常態分佈。

為了重建資料,模型需學習該過程的反向映射。反向過程定義如下: \begin{align} \begin{cases} p_\theta(\boldsymbol{x}_{0}) = p(\boldsymbol{x}_T)\prod_{t=1}^T p_\theta(\boldsymbol{x}_{t-1}|\boldsymbol{x}_t); \\ p_\theta(\boldsymbol{x}_{t-1} | \boldsymbol{x}_t) = \mathcal{N}(\boldsymbol{x}_{t-1}; \mu_\theta(\boldsymbol{x}_t, t), \Sigma_\theta(\boldsymbol{x}_t, t)), \end{cases} \end{align} 其中, $p(\boldsymbol{x}_T) = \mathcal{N}(\boldsymbol{x}_T; 0, \boldsymbol{I})$ 為標準常態分佈; $\mu_\theta$ 與 $\Sigma_\theta$ 為經參數 $\theta$ 的神經網路所參數化的平均值與共變異數矩陣。

然而,直接對反向過程的均值 $\mu_\theta$ 進行建模在實作上往往難以收斂。因此, Ho, Jain, and Abbeel (2020) 提出了一種稱為去噪擴散機率模型(Denoising Diffusion Probabilistic Models, DDPM)的參數化方式,將 $p_\theta(\boldsymbol{x}_{t-1}|\boldsymbol{x}_t)$ 重新參數化為 \begin{align} \mu_\theta(\boldsymbol{x}_t, t) & = \frac{1}{\sqrt{\alpha_t}} \Big(\boldsymbol{x}_t-\frac{\beta_t}{\sqrt{1-\bar{\alpha}_t}}\, \boldsymbol{\epsilon}_{\theta}(\boldsymbol{x}_t,t)\Big), \\ \Sigma_\theta(\boldsymbol{x}_t, t) & = \sigma_t^2 \boldsymbol{I}, \qquad \sigma_t^2 = \beta_t \text{ or } \sigma_t^2 = \frac{1 - \bar{\alpha}_{t - 1}}{1 - \bar{\alpha}_t} \beta_t, \end{align} 其中, $\alpha_t = 1 - \beta_t$ ; $\bar{\alpha}_t = \prod_{s=1}^t \alpha_s$ 。在此架構下, $\boldsymbol{\epsilon}_{\theta}(\boldsymbol{x}_t,t)$ 是用於估計 $\boldsymbol{x}_t$ 於正向擴散過程中所加入的隨機高斯雜訊,並透過移除 $\boldsymbol{x}_t$ 中該雜訊的估計值重建 $\boldsymbol{x}_{t-1}$ 。此參數化策略避免了對複雜高維資料分佈的直接建模,大幅簡化訓練目標並提升數值穩定性。因此,任意擴散步驟 $t$ 之樣本可表述為原始資料 $\boldsymbol{x}_0$ 與雜訊之線性組合 \begin{align} \boldsymbol{x}_t = \sqrt{\bar{\alpha}_t} \boldsymbol{x}_0 + \sqrt{1 - \bar{\alpha}_t} \boldsymbol{\epsilon}, \qquad \boldsymbol{\epsilon} \sim \mathcal{N}(0, \boldsymbol{I}), \end{align} 使得模型能夠在訓練期間隨機抽樣時間步 $t$ 與雜訊 $\boldsymbol{\epsilon}$ ,而無需迭代計算中間步驟。由於此表示方式將原先對資料生成分佈的直接擬合問題轉化為對高斯雜訊的估計,訓練目標可進一步簡化為 \begin{align} \min_\theta \mathbb{E}_{t, \boldsymbol{x}_0, \boldsymbol{\epsilon}} \left[ \left\| \boldsymbol{\epsilon} - \boldsymbol{\epsilon}_\theta \left(\sqrt{\bar{\alpha}_t} \boldsymbol{x}_0 + \sqrt{1 - \bar{\alpha}_t} \boldsymbol{\epsilon}, t\right) \right\|^2 \right], \end{align} 即僅需最小化預測雜訊與實際注入雜訊之間的平均平方誤差(Mean Squared Error, MSE)。

在時間序列插補應用中, Alcaraz and Strodthoff (2023) 進一步提出條件式擴散模型(Conditional Diffusion Model),僅對缺失片段施加擴散與去噪操作。於訓練階段,模型接收部分觀測序列作為條件輸入以學習重建完整序列;而於生成階段,藉由固定已觀測部分並執行反向去噪過程,即可在保持時間一致性的前提下補全缺失區段。此方法結合了擴散模型的穩定性與高品質生成能力,能有效處理具長期依賴或結構性缺失的時間序列資料。

2.1.6 基於 S4 層的結構化狀態空間擴散模型

結構化狀態空間擴散模型(Structured State Space Diffusion, SSSD)由 Alcaraz and Strodthoff (2023) 提出,旨在結合 DiffWave 架構擴散模型 (Kong et al., 2021) 的生成穩定性與 SSM 的長期依賴建模能力。該模型以條件擴散(conditional diffusion)的形式應用於時間序列插補任務,僅對缺失區段施加雜訊,而保留已觀測部分不受擾動,藉此避免資料洩漏並維持條件一致性。透過在每一步中學習反向去噪過程, SSSD 能在固定觀測條件下逐步重建缺失資料,展現高品質的插補與生成效果。

在此基礎上, Alcaraz and Strodthoff (2023) 進一步提出基於 S4 層的結構化狀態空間擴散模型(Structured State Space Diffusion with S4 Layers, $\text{SSSD}^\text{S4}$),作為針對時間序列任務的改進變體。 $\text{SSSD}^\text{S4}$ 延續 SSSD 的條件擴散設計,並將原本 DiffWave 架構中的雙向擴張卷積層(Bidirectional Dilated Convolution Layer)替換為 S4 層,以提升模型對長期序列動態的建模能力。實驗結果顯示,該變體在多種缺失模式下均表現出更穩定且準確的插補效果,較以傳統卷積或 Transformer 基礎的擴散模型更穩定且具有更好插補表現 (Alcaraz and Strodthoff, 2023) 。

2.2 空間統計學

2.2.1 克利金法

克利金(Kriging)方法源自地質學中的空間統計學,是依據觀測資料進行線性內插,用以估計未觀測地點的隨機場值 (N. A. C. Cressie, 1993) 。其假設位於 $\boldsymbol{s}$ 的觀測值 $Z(\boldsymbol{s})$ 可表示為: \begin{align} Z(\boldsymbol{s}) = Y(\boldsymbol{s}) + \varepsilon(\boldsymbol{s}), \qquad \boldsymbol{s} \in D \subset \mathbb{R}^D, \end{align} 其中, $Y(\boldsymbol{s}) = \mu(\boldsymbol{s}) + \xi(\boldsymbol{s})$ 是隨空間位置變化的線性均值結構; $\varepsilon(\boldsymbol{s})$ 為零均值且與 $Y(\boldsymbol{s})$ 不相關的隨機雜訊,其共變異數函數為 $C(\boldsymbol{s}, \boldsymbol{s}’) = \text{Cov}(\varepsilon(\boldsymbol{s}), \varepsilon(\boldsymbol{s}’))$,可為非平穩的空間共變異數函數。

傳統克利金方法依賴對完整共變異數矩陣的逆運算,導致當觀測點數量 $n$ 增大時,計算成本呈現指數級增長,形成顯著的計算瓶頸 (N. Cressie and Johannesson, 2008) 。固定秩克利金(Fixed Rank Kriging, FRK)以及其後的自適應變體(Resolution Adaptive Fixed Rank Kriging, AFRK)即於此背景下提出,以有效降低大規模空間資料分析的計算負擔。

2.2.2 固定秩克利金法

為降低克利金法的計算負擔, N. Cressie and Johannesson (2008) 提出固定秩克利金法(Fixed Rank Kriging, FRK),將隨機場以有限基底函數展開,使高維空間隨機效應近似為低維隨機係數: \begin{align} Y(\boldsymbol{s}) = \mu(\boldsymbol{s}) + \boldsymbol{f}(\boldsymbol{s})^\top \boldsymbol{w} + \xi(\boldsymbol{s}), \end{align} 其中, $\boldsymbol{f}(\boldsymbol{s}) = (f_1(\boldsymbol{s}), \dots, f_K(\boldsymbol{s}))^\top$ 為預先指定的 $K$ 維基底函數向量, $K \le n$; $\boldsymbol{w} \sim N(\boldsymbol{0}, \boldsymbol{M})$ , $\boldsymbol{M}$ 為未知的非負定矩陣; $\xi(\boldsymbol{s}) \sim \mathcal{N}(0, \sigma_\xi^2)$ 為細尺度隨機雜訊。此時共變異數矩陣可表示為: \begin{align} \text{Cov}(Y(\boldsymbol{s}), Y(\boldsymbol{s}’)) = \boldsymbol{f}(\boldsymbol{s})^\top \boldsymbol{M} \boldsymbol{f}(\boldsymbol{s}’) + \sigma_\xi^2 \boldsymbol{I}(\boldsymbol{s} = \boldsymbol{s}’). \end{align}

由於 $\boldsymbol{fMf}^\top$ 的秩通常遠小於觀測點數量 $n$, FRK 可有效降低共變異數矩陣逆運算的計算成本,特別適合處理大規模遙測與環境觀測資料。

2.2.3 自適應固定秩克利金法

在 FRK 的基礎上, Tzeng and Huang (2018) 進一步提出自適應固定秩克利金法(Resolution Adaptive Fixed Rank Kriging, AFRK)。其想法是讓基底函數的解析度能夠依據資料分佈自動調整,以更靈活地捕捉空間非均勻性,使模型能針對不同區域分配適當的空間解析度。

AFRK 所使用的基底函數為多解析度薄板樣條基底函數(Multi-Resolution Thin-Plate Spline Basis Functions, MRTS),其係由薄板樣條(Thin-Plate Spline, TPS)所構成。 TPS 為一種常見的平滑樣條方法,透過最小化平方誤差與平滑懲罰項以獲得平滑函數 (Wahba and Wendelberger, 1980; Green and Silverman, 1993) 。在 TPS 的基礎上, Tzeng and Huang (2018) 進一步透過特徵值分解建構出一組具有不同解析度的有序基底函數,並稱其為 MRTS 。為了自動適應資料分佈與空間非均勻性, AFRK 可依據特徵值大小自動選擇使用的基底數量,只保留對大部分變異有解釋力的基底,實現以較少的基底函數捕捉主要空間變異、提高計算效率。

在 AFRK 中, MRTS 定義為: \begin{align} f_k(\boldsymbol{s}) = \begin{cases} 1, & \text{if } k = 1; \\ x_{k-1}, & \text{if } k = 2, \ldots, d+1; \\ \lambda^{-1}_{k-d-1} \times \left \{ \phi(\boldsymbol{s}) - \boldsymbol{\Phi} \boldsymbol{X} (\boldsymbol{X}’\boldsymbol{X})^{-1} \boldsymbol{x} \right \}’ \boldsymbol{v}_{k-d-1}, & \text{if } k = d+2, \ldots, n, \end{cases} \end{align} 其中, $f_k$ 是 $\boldsymbol{f}$ 中第 $k$ 個基底; $\boldsymbol{X} \in \mathbb{R}^{N \times (d+1)}$ 是設計矩陣,其每一列對應觀測位置 $\boldsymbol{s}_i$ 的截距與坐標, $x = (1, \boldsymbol{s}’)’ = (1, x_1, \ldots, x_d)’$ ;而 $\boldsymbol{\Phi}$ 滿足 \begin{align} J(f) = \boldsymbol{\alpha}’ \boldsymbol{\Phi \alpha}, \end{align} 其中, $\boldsymbol{\Phi}$ 是 $n \times n$ 矩陣,其 $(i,j)$ 元素為 $\phi_j(\boldsymbol{s}_i)$ , $\phi(\boldsymbol{s})$ 定義如下: \begin{align} \phi_i (\boldsymbol{s}) = \begin{cases} \frac{1}{12} \| \boldsymbol{s} - \boldsymbol{s}_i \|^3, & \text{if } d = 1; \\ \frac{1}{8\pi} \| \boldsymbol{s} - \boldsymbol{s}_i \|^2 \log(\| \boldsymbol{s} - \boldsymbol{s}_i \|), & \text{if } d = 2; \\ -\frac{1}{8} \| \boldsymbol{s} - \boldsymbol{s}_i \|, & \text{if } d = 3; \end{cases} \end{align} 而 $\boldsymbol{v}_{k}$ 為矩陣 $\boldsymbol{V}$ 的第 $k$ 列,且 $\boldsymbol{V} \operatorname{diag}(\lambda_1, \ldots, \lambda_n) \boldsymbol{V}’$ 為 $\boldsymbol{Q \Phi Q}$ 的特徵分解,其中 $\boldsymbol{Q} = \boldsymbol{I} - \boldsymbol{X} (\boldsymbol{X}’\boldsymbol{X})^{-1} \boldsymbol{X}’$。此方法能依據資料密度調整基底,在資料較密集的區域產生較高解析度的基底,而在資料稀疏的區域則保持較平滑的結構。

藉由保留低秩近似的計算優勢,同時引入資料導向的多解析度基底, AFRK 能有效處理非均勻採樣與局部變異性顯著的空間資料。相較於傳統 FRK , AFRK 能在非均勻與非平穩的空間資料中展現更佳的預測表現。

2.3 時空模型

2.3.1 高斯過程

高斯過程(Gaussian Process, GP)為一種非參數貝氏模型,用於建立輸入空間至輸出值間之隨機函數關係 (Rasmussen and Williams, 2006) 。 GP 由均值函數與共變異數函數(covariance function,或稱核函數 kernel function)所定義,能在貝氏框架下提供預測分佈之均值與不確定性量化。然而,在標準 GP 中,後驗推論需對規模為 $N \times N$ 之共變異數矩陣進行反矩陣運算或 Cholesky 分解,其計算複雜度為 $\mathcal{O}(N^3)$ ,因此難以直接應用於大規模時空資料情境。

為克服計算瓶頸, Hensman, Fusi, and Lawrence (2013) 提出隨機變分高斯過程(Stochastic Variational Gaussian Process, SVGP)。該方法引入一組數量為 $M \ll N$ 之誘導點(inducing points)以近似完整後驗分佈。在 Gardner et al. (2021) 的實作框架中,變分分佈(variational distribution)被定義為具備全共變異數矩陣(full covariance matrix)之多元常態分佈。為了在優化過程中確保共變異數矩陣之正定性(positive definite),實務上將其參數化為均值向量與下三角矩陣(lower triangle),即透過 Cholesky 分解進行參數化 (Hensman, Matthews, and Ghahramani, 2014) 。

在 SVGP 框架下,變分分佈的大小由誘導點數量決定,即變分均值之維度為 $M$,而變分共變異數矩陣之規模為 $M \times M$ 。透過 $M \ll N$ 的誘導點以近似完整後驗,並以變分證據下界(Evidence Lower Bound, ELBO)作為目標函數, SVGP 將計算複雜度降至 $\mathcal{O}(M^3)$ (Gardner et al., 2021) 。

2.3.2 時空深度克利金法

時空深度克利金法(Spatio-temporal DeepKriging, STDK)為近年提出之結合深度學習與空間統計之非參數化建模方法,主要用於大規模時空資料之內插與機率性預測問題 (Nag, Sun, and Reich, 2023) 。相較於傳統高斯過程需事前假設共變異數函數,且在樣本數增加時面臨 $\mathcal{O}(N^3)$ 計算複雜度之限制,STDK 採用資料驅動(data-driven)方式學習時空相關結構,藉此提升在大規模資料情境下之適用性。

STDK 透過嵌入層(Embedding Layer)將時空位置 $(\mathbf{s}, t)$ 映射至高維特徵空間。令基底函數向量為 \begin{align} \boldsymbol{\phi}(\mathbf{s}, t) = \left[\phi_1(\mathbf{s}, t), \dots, \phi_K(\mathbf{s}, t)\right]^\top, \end{align} 其中 $\{\phi_k(\cdot)\}_{k=1}^K$ 為具多尺度(multi-resolution)特性之基底函數。既有研究多採用具緊支撐性質之溫德蘭函數(Wendland functions)或徑向基底函數(Radial Basis Functions, RBF),以捕捉不同尺度之空間依賴關係 (Nag, Sun, and Reich, 2023) 。

在此基礎上,嵌入後之特徵將輸入深度神經網路以進行非線性映射,其模型形式可表示為 \begin{align} Z(\mathbf{s}, t) = f_{\boldsymbol{\theta}}\big(\boldsymbol{\phi}(\mathbf{s}, t)\big) + \epsilon, \end{align} 其中 $f_{\boldsymbol{\theta}}(\cdot)$ 為深度神經網路所表示之函數; $\boldsymbol{\theta}$ 為模型參數; $\epsilon$ 為隨機誤差項。

在預測方面, STDK 為進行機率性預測(probabilistic forecasting),採用分位數損失函數(quantile loss)作為訓練目標,以估計條件分佈之不同分位數,進而構建預測區間。相較於傳統以均方誤差(MSE)為目標之方法,此作法能提供不確定性量化之資訊。

此外,STDK 透過結合基底函數嵌入(basis function embeddings)與深度神經網路,避免對共變異數矩陣進行顯式分解,因此能有效擴展至大型資料集,並成為大規模時空建模任務中具代表性的基準方法之一。

3 研究方法

本研究提出一套結合深度時間序列建模與空間統計建模之時空預測框架。透過在統一架構下同時納入時間動態與空間相關結構,所提出之方法可用於重建缺失觀測資料,並進一步進行未來時間點之預測。

3.1 時空預測任務

考慮空間域中存在位置集合 $\mathcal{S}$,並可將其劃分為具備觀測資料之已觀測位置集合 $\mathcal{S}_{\mathrm{observed}}$,以及無任何觀測紀錄之未觀測位置集合 $\mathcal{S}_{\mathrm{unobserved}}$。對於已觀測位置 $\boldsymbol{s} \in \mathcal{S}_{\mathrm{observed}}$,其於時間 $t \in \{1, \dots, T\}$ 之目標變數均已完整觀測,將其觀測值記為 $y_t(\boldsymbol{s})$;相反地,未觀測位置 $\boldsymbol{s}^\ast \in \mathcal{S}_{\mathrm{unobserved}}$ 則完全缺乏此時間段內之任何歷史資料。本研究旨在建立一預測模型,利用已觀測位置之歷史時空特徵 $\boldsymbol{Y}_{1:T}(\mathcal{S}_{\mathrm{observed}})$,實現對未觀測位置於未來時間點 $t > T$ 之目標變數值 $\hat{y}_{t}(\boldsymbol{s}^\ast)$ 的時間外插與空間推估。

為利後續模型推導與方法說明,表 1 彙整本研究方法使用之主要符號及其定義。

表 1: 本研究方法所使用之主要符號說明。

符號說明
$\mathcal{S}$空間域(Spatial domain)。
$\mathcal{S}_{\mathrm{observed}}$已觀測位置集合。
$\mathcal{S}_{\mathrm{unobserved}}$未觀測位置集合。
$\boldsymbol{s}$已觀測空間位置。
$\boldsymbol{s}^{\ast}$未觀測空間位置。
$N$已觀測位置數量。
$T$已觀測時間序列長度。
$y_t(\boldsymbol{s})$位置 $\boldsymbol{s}$ 於時間點 $t$ 之觀測值。
$\boldsymbol{Y}_{1:T}(\mathcal{S}_{\mathrm{observed}})$所有已觀測位置於時間區間 $1,\ldots,T$ 之觀測資料集合。
$\tilde{y}^{(T)}_t(\boldsymbol{s})$位置 $\boldsymbol{s}$ 於時間點 $t$ 之時間標準化觀測值。
$\hat{\tilde{y}}^{(T)}_t(\boldsymbol{s})$$\text{SSSD}^{\text{S4+MRTS}}$ 於時間標準化尺度下之預測結果。
$\tilde{y}^{(S)}_t(\boldsymbol{s})$已觀測位置 $\boldsymbol{s}$ 於時間點 $t$ 之空間標準化預測值。
$\hat{\tilde{y}}^{(S)}_t(\boldsymbol{s}^{\ast})$AFRK 於空間標準化尺度下對未觀測位置之預測結果。
$\hat{y}_t(\boldsymbol{s})$已觀測位置 $\boldsymbol{s}$ 於原始尺度下之預測值。
$\hat{y}_t(\boldsymbol{s}^{\ast})$未觀測位置 $\boldsymbol{s}^{\ast}$ 於原始尺度下之預測值。
$\mu(\boldsymbol{s})$位置 $\boldsymbol{s}$ 時間序列之平均值。
$\sigma(\boldsymbol{s})$位置 $\boldsymbol{s}$ 時間序列之標準差。
$\bar{y}_t$時間點 $t$ 所有已觀測位置預測值之空間平均。
$\tau_t$時間點 $t$ 所有已觀測位置預測值之空間標準差。
$\boldsymbol{\phi}(\boldsymbol{s})$位置 $\boldsymbol{s}$ 對應之多解析度薄板樣條基底函數(MRTS)特徵向量。
$\boldsymbol{\Phi}$由所有已觀測位置建構之 MRTS 基底矩陣。
$\tilde{\boldsymbol{\Phi}}$沿時間維度擴展後之 MRTS 基底張量。
$D_{\mathrm{MRTS}}$MRTS 基底函數數量。
$\boldsymbol{m}_t(\boldsymbol{s})$表示位置 $\boldsymbol{s}$ 於時間點 $t$ 是否存在觀測值之二元遮罩。
$\mathbf{c}_t(\boldsymbol{s})$輸入至 S4 層之條件向量。
$\theta$擴散模型之可訓練參數。
$g_\theta(\cdot)$反向去噪過程中使用之雜訊估計網路。
$\boldsymbol{x}_t$擴散步驟 $t$ 之潛在變數。
$\boldsymbol{\epsilon}$正向擴散過程中加入之高斯雜訊。
$\hat{\boldsymbol{\epsilon}}_\theta$由網路 $g_\theta$ 預測之擴散雜訊。
$\alpha_t$擴散步驟 $t$ 之變異數排程參數。
$\bar{\alpha}_t$截至步驟 $t$ 為止之擴散係數累積乘積。
$K$原始特徵表示之通道數。
$C$條件卷積層之輸出通道數。
$\boldsymbol{f}(\boldsymbol{s})$AFRK 使用之空間基底向量。
$\hat{\boldsymbol{w}}_t$時間點 $t$ 估計得到之 AFRK 基底係數。
$\hat{\xi}_t(\boldsymbol{s})$位置 $\boldsymbol{s}$ 之 AFRK 殘差成分估計值。

此研究目標可形式化為建構一映射函數 $f(\cdot)$ \begin{align} \hat{y}_{t}(\boldsymbol{s}^\ast) = f \left (\boldsymbol{Y}_{1 : T}(\mathcal{S}_{\mathrm{observed}}), \boldsymbol{s}^\ast \right ), \qquad \boldsymbol{s}^\ast\in\mathcal{S}_{\mathrm{unobserved}}, ~ t > T, \end{align} 其中, $\boldsymbol{Y}_{1 : T}(\mathcal{S}_{\mathrm{observed}})$ 表示截至時間 $T$ ,由已觀測位置所組成的觀測序列集合; $\hat{y}_{t}(\boldsymbol{s}^\ast)$ 則為對未觀測位置 $\boldsymbol{s}^\ast$ 在未來時間 $t$ 的估計值。

3.2 時空建模方法

為實現上述映射目標,本研究將 $\text{SSSD}^{\text{S4}}$ 模型中的 S4 層結合 MRTS 基底函數,稱為 $\text{SSSD}^{\text{S4+MRTS}}$ ,使時間序列資訊進行卷積訓練時能同時考慮空間相依性;而在推論階段,則將 $\text{SSSD}^{\text{S4+MRTS}}$ 所推論之已觀測地點預測值與未觀測地點之空間坐標輸入 AFRK ,綜合時空雙維度資訊估計未觀測位置之目標值。

3.2.1 基於 $\text{SSSD}^{\text{S4}}$ 之時間建模

為有效捕捉資料於時間序列上的相依結構,本研究採用 $\text{SSSD}^{\text{S4}}$ 作為時間特徵模型。其目的為藉由擴散模型與 S4 層,從不完整或含有雜訊的時間序列中,學習兼具長期依賴與局部動態之潛在時間特徵。

以已觀測位置 $\mathcal{S}_{\mathrm{observed}}$ 之觀測資料作為訓練基礎,對資料進行標準化處理,使得 \begin{align} \tilde{y}^{(T)}_t(\boldsymbol{s}) = \frac{y_t(\boldsymbol{s}) - \mu(\boldsymbol{s})}{\sigma(\boldsymbol{s})}, \qquad \boldsymbol{s} \in \mathcal{S}_{\mathrm{observed}}, \end{align} 其中,$\mu(\boldsymbol{s})$ 與 $\sigma(\boldsymbol{s})$ 分別表示位置 $\boldsymbol{s}$ 在整段時間序列上的平均值與標準差,即 \begin{align} \mu(\boldsymbol{s}) & = \frac{1}{T} \sum_{t=1}^{T} y_t(\boldsymbol{s}), \\ \sigma(\boldsymbol{s}) & = \sqrt{\frac{1}{T}\sum_{t=1}^{T} \left(y_t(\boldsymbol{s})-\mu(\boldsymbol{s})\right)^2 }. \end{align} 隨後,將標準化後的輸入序列 $\tilde{y}^{(T)}_t(\boldsymbol{s})$ 送入 $\text{SSSD}^{\text{S4}}$ 模型,以學習其時間依賴結構。

為將空間結構資訊有效融入時間序列建模過程,這裡將 MRTS 作為條件輸入之一,並與觀測特徵及遮罩資訊共同構成 S4 層之條件通道。對於空間位置 $\boldsymbol{s} \in \mathcal{S}_{\mathrm{observed}}$ ,其 MRTS 特徵向量定義為 $\boldsymbol{\phi}(\boldsymbol{s}) \in \mathbb{R}^{D_{\mathrm{MRTS}}}$ 。所有觀測位置所組成之 MRTS 基底矩陣為 \begin{align} \boldsymbol{\Phi} = \begin{bmatrix} \boldsymbol{\phi}(\boldsymbol{s}_1)^\top\\ \boldsymbol{\phi}(\boldsymbol{s}_2)^\top\\ \vdots\\ \boldsymbol{\phi}(\boldsymbol{s}_N)^\top \end{bmatrix} \in \mathbb{R}^{N \times D_{\mathrm{MRTS}}}, \end{align} 其中, MRTS 特徵維度為 \begin{align} D_{\mathrm{MRTS}} = \min \left( N, \operatorname{round}\left( 10\sqrt{N} \right) \right), \end{align} 其中 $N$ 表示已觀測位置數量。此設定採用 AFRK 相關實作中所建議之經驗法則 (Tzeng, Huang, Wang, Nychka, et al., 2021) 。透過限制基底函數數量隨觀測位置數增加而適度成長,可在空間表示能力與計算效率之間取得平衡,並避免於大規模空間資料中產生過多基底函數。

為使 MRTS 基底能與時間序列模型整合,將 $\boldsymbol{\Phi}$ 視為不隨時間變化之空間基底,並沿時間維度重複展開,使其與時間序列長度 $T$ 對齊。展開後之張量為 \begin{align} \tilde{\boldsymbol{\Phi}} \in \mathbb{R}^{N \times D_{\mathrm{MRTS}} \times T}. \end{align}

將空間特徵與時間標準化後之觀測資料 $\tilde{y}^{(T)}_t(\boldsymbol{s})$ 以及缺失資料遮罩 $\boldsymbol{m}_t(\boldsymbol{s})$ 進行串接(concatenation)後, S4 層之條件通道可由原始形式 $\left[ \tilde{y}^{(T)}_t(\boldsymbol{s}), \boldsymbol{m}_t(\boldsymbol{s}) \right]$ 擴展為 \begin{align} \mathbf{c}_t(\boldsymbol{s}) = \left[ \tilde{y}^{(T)}_t(\boldsymbol{s}), \boldsymbol{m}_t(\boldsymbol{s}), \tilde{\boldsymbol{\Phi}}(\boldsymbol{s}) \right]. \end{align}

此條件向量 $\mathbf{c}_t(\boldsymbol{s})$ 與擴散過程之噪聲輸入共同作為 S4 層之輸入,以進行時間序列之狀態建模。因此, S4 層之條件通道數由原架構的 $2K$ 擴增為 \begin{align} (2 + D_{\mathrm{MRTS}}) K, \end{align} 相應地,條件卷積層(conditional convolution layer)調整為 \begin{align} \mathrm{Conv}\left((2 + D_{\mathrm{MRTS}})K, 2C\right), \end{align} 其中, $K$ 為原始輸入特徵的通道數(input channels)。如此可使 MRTS 空間基底特徵得以融入後續 S4 層之狀態空間轉換過程。

透過條件通道(conditional channel)機制, MRTS 所提供之空間基底資訊得以輸入至 S4 層,使空間資訊能夠參與擴散模型之生成過程,進而在統一生成架構下同時學習時間依賴關係與空間結構特徵。經反向過程後,模型可產生已觀測位置之時間標準化預測結果 \begin{align} \hat{\tilde{y}}^{(T)}_t(\boldsymbol{s}), \qquad \boldsymbol{s}\in\mathcal S_{\mathrm{observed}}, \end{align} 其後再轉換回原始量測尺度。 由於模型訓練與推論皆於標準化資料上進行,因此最終預測結果需透過反標準化(inverse standardization)還原至原始量測尺度: \begin{align} \hat{y}_{t}(\boldsymbol{s}) = \sigma(\boldsymbol{s}) \hat{\tilde{y}}^{(T)}_{t}(\boldsymbol{s}) + \mu(\boldsymbol{s}), \end{align} 其中, $\mu(\boldsymbol{s})$ 與 $\sigma(\boldsymbol{s})$ 分別表示位置 $\boldsymbol{s}$ 時間序列之平均值與標準差。此轉換可將模型輸出還原至原始資料尺度,以利後續分析與空間插值應用。

SSSD 模型採用平均平方誤差(Mean Squared Error, MSE)作為訓練時的損失函數,並透過迭代更新參數 $\theta$ 進行模型學習,以確保所擷取之時間特徵具有穩定且具預測能力之表示。在缺失值重建與多步未來預測等任務中,$\text{SSSD}^{\text{S4}}$ 所學得之時間表示(temporal representations)可提供具有辨識能力的序列特徵,而基於 MRTS 所建構之空間表示(spatial representations)則提供額外的空間資訊,進而有助於提升時空估計之表現。

3.2.2 基於 AFRK 之空間建模

在本研究所提出之整合架構中, AFRK 將 $\text{SSSD}^{\text{S4+MRTS}}$ 所產生之已觀測位置預測結果 $\hat{y}_{t}(\boldsymbol{s})$($\boldsymbol{s}\in\mathcal{S}{\mathrm{observed}}$)作為沿時間維度之輸入資訊。為消除不同位置資料尺度之差異,首先將同一時間點 $t$ 之預測值於已觀測位置間進行空間標準化: \begin{align} \tilde{y}^{(S)}_t(\boldsymbol{s}) = \frac{\hat{y}_t(\boldsymbol{s}) - \bar{y}_t}{\tau_t}, \qquad \boldsymbol{s} \in \mathcal{S}_{\mathrm{observed}}, \end{align} 其中, $\bar{y}_t$ 與 $\tau_t$ 分別為時間點 $t$ 於所有已觀測位置 $\mathcal{S}_{\mathrm{observed}}$ 之空間平均與空間標準差,其定義如下: \begin{align} \bar{y}_t & = \frac{1}{N} \sum_{\boldsymbol{s} \in \mathcal{S}_{\mathrm{observed}}} \hat{y}_t(\boldsymbol{s}), \\ \tau_t & = \sqrt{\frac{1}{N} \sum_{\boldsymbol{s} \in \mathcal{S}_{\mathrm{observed}}} \left(\hat{y}_t(\boldsymbol{s}) - \bar{y}_t\right)^2}. \end{align}

而後,以空間標準化後之 $\tilde{y}^{(S)}_t(\boldsymbol{s})$ 作為 AFRK 之輸入,對未觀測位置 $\boldsymbol{s}^\ast \in \mathcal{S}_{\mathrm{unobserved}}$ 進行條件推估。於空間標準化尺度下,未觀測位置之預測值可表示為 \begin{align} \hat{\tilde{y}}^{(S)}_t(\boldsymbol{s}^\ast) = \boldsymbol{f}(\boldsymbol{s}^\ast)^\top \hat{\boldsymbol{w}}_t + \hat{\xi}_t(\boldsymbol{s}^\ast), \end{align} 其中, $\hat{\boldsymbol{w}}_t$ 與 $\hat{\xi}_t(\boldsymbol{s}^{\ast})$ 分別為 AFRK 根據已觀測位置資訊 $\tilde{y}^{(S)}_t(\boldsymbol{s})$ 所估計之空間條件分布參數。

最後,將空間推估結果反標準化,還原至原始量測尺度 \begin{align} \hat{y}_t(\boldsymbol{s}^\ast) = \tau_t \hat{\tilde{y}}^{(S)}_t(\boldsymbol{s}^\ast) + \bar{y}_t. \end{align} 至此,即完成未觀測位置之空間推估以及整體時空場之重建與預測。

3.3 模型訓練與推論演算法

本節總結上述時空整合模型之運作流程。演算法 1 說明 $\text{SSSD}^{\text{S4+MRTS}}$ 之訓練過程;演算法 2 則展示 $\text{SSSD}^{\text{S4+MRTS}}$ 結合 AFRK 進行推論之完整步驟。

演算法 1 訓練階段
Require: $y_t(\boldsymbol{s}), \boldsymbol{s} \in \mathcal{S}_{\mathrm{observed}}$; Max steps $T$; Learning rate $\eta$; Model parameters $\theta$
1. $\tilde{y}^{(T)}_t(\boldsymbol{s}) = \left(y_t(\boldsymbol{s}) - \mu(\boldsymbol{s})\right) / \sigma(\boldsymbol{s})$
2. $\tilde{\boldsymbol{\Phi}} \in \mathbb{R}^{N \times D_{\mathrm{MRTS}} \times T}$
3. repeat
4. $t \sim \text{Uniform}(\{1, \dots, T\}), \quad \boldsymbol{\epsilon} \sim \mathcal{N}(0, \boldsymbol{I})$
5. $\boldsymbol{x}_t = \sqrt{\bar{\alpha}_t} \tilde{y}^{(T)}_0 + \sqrt{1 - \bar{\alpha}_t} \boldsymbol{\epsilon}$
6. $\mathbf{c}_t = \left[\tilde{y}^{(T)}_t(\boldsymbol{s}), \boldsymbol{m}_t(\boldsymbol{s}), \tilde{\boldsymbol{\Phi}}(\boldsymbol{s})\right]$
7. $\hat{\boldsymbol{\epsilon}}_{\theta} = g_{\theta}(\boldsymbol{x}_t, \mathbf{c}_t)$
8. $\mathcal{L} = \| \boldsymbol{\epsilon} - \hat{\boldsymbol{\epsilon}}_{\theta} \|^2$
9. $\theta \leftarrow \theta - \eta \nabla_{\theta} \mathcal{L}$
10. until convergence
演算法 2 推論階段
Require: $y_t(\boldsymbol{s}), \boldsymbol{s} \in \mathcal{S}_{\mathrm{observed}}$; $\boldsymbol{s}^* \in \mathcal{S}_{\mathrm{unobserved}}$; Trained parameters $\theta$
1. $\tilde{y}^{(T)}_t(\boldsymbol{s}) = \left(y_t(\boldsymbol{s}) - \mu(\boldsymbol{s})\right) / \sigma(\boldsymbol{s})$
2. $\tilde{\boldsymbol{\Phi}} \in \mathbb{R}^{N \times D_{\mathrm{MRTS}} \times T}$
3. $\boldsymbol{x}_T \sim \mathcal{N}(0, \boldsymbol{I})$
4. for $t = T, \dots, 1$ do
5. $\mathbf{c}_t(\boldsymbol{s}) = \left[\tilde{y}^{(T)}_t(\boldsymbol{s}), \boldsymbol{m}_t(\boldsymbol{s}), \tilde{\boldsymbol{\Phi}}(\boldsymbol{s})\right]$
6. $\hat{\boldsymbol{\epsilon}}_{\theta} = g_{\theta}(\boldsymbol{x}_t, \mathbf{c}_t)$
7. $\boldsymbol{x}_{t-1} = \frac{1}{\sqrt{\alpha_t}} \left( \boldsymbol{x}_t - \frac{1-\alpha_t}{\sqrt{1-\bar{\alpha}_t}} \hat{\boldsymbol{\epsilon}}_{\theta} \right) + \sqrt{1-\alpha_t} \mathbf{z}, \quad \mathbf{z} \sim \mathcal{N}(0, \boldsymbol{I})$
8. end for
9. $\hat{\tilde{y}}^{(T)}_t(\boldsymbol{s}) = \boldsymbol{x}_0$
10. $\hat{y}_t(\boldsymbol{s}) = \sigma(\boldsymbol{s}) \hat{\tilde{y}}^{(T)}_t(\boldsymbol{s}) + \mu(\boldsymbol{s}), \quad \boldsymbol{s} \in \mathcal{S}_{\mathrm{observed}}$
11. $\tilde{y}^{(S)}_t(\boldsymbol{s}) = (\hat{y}_t(\boldsymbol{s}) - \bar{y}_t) / \tau_t, \quad \boldsymbol{s} \in \mathcal{S}_{\mathrm{observed}}$
12. $\hat{\tilde{y}}^{(S)}_t(\boldsymbol{s}^*) = \text{AFRK}(\tilde{y}^{(S)}_t(\boldsymbol{s}) \mid \boldsymbol{s}^*), \quad \forall \boldsymbol{s}^* \in \mathcal{S}_{\mathrm{unobserved}}$
13. $\hat{y}_t(\boldsymbol{s}^*) = \tau_t \hat{\tilde{y}}^{(S)}_t(\boldsymbol{s}^*) + \bar{y}_t$
14. return Full field $\hat{y}_t(\mathcal{S}_{\mathrm{observed}} \cup \mathcal{S}_{\mathrm{unobserved}})$

4 實驗

本章介紹所提出時空整合框架之實驗設計與設定,旨在評估其於不同觀測條件下對未觀測時空值進行預測之有效性。實驗採用三個主要資料集,並詳細說明模型配置、訓練與推論流程,以及各資料集之特性與資料切分方式,以驗證所提出框架在不同資料條件下之預測表現。

4.1 資料集

本研究採用兩個具代表性的氣象資料集,以及一個模擬資料集作為實驗資料。其內容涵蓋地面觀測站觀測資料(in-situ ground station observations)、全球再分析資料(global reanalysis data),以及可控制之模擬空間過程(synthetic spatial processes)。藉由不同資料來源與觀測型態之設計,本研究得以評估所提出方法於異質觀測情境下之時空建模表現。

4.1.1 Weather2K

Weather2K 是近期提出的多變量地面觀測基準資料集,收錄中國境內數千個氣象站的觀測資料,時間解析度為 3 小時,涵蓋氣溫、氣壓、濕度、風速等近地表氣象變數 (Zhu et al., 2023) 。開源版本 Weather2K-R (下稱 Weather2K)包含 1,866 個氣象站及 13,632 個連續時間步,時間序列完整且定期採樣,無缺失值。資料同時提供位置常數(經緯度與海拔),方便時空模式分析。

Weather2K 資料集呈現空間分佈異質及測站密度不均,涵蓋都市、平原與山區等多種地理環境與氣候條件,提供豐富的時空變化訊號。這些特性使 Weather2K 適合用於時空插值、短期序列預測,以及模型泛化能力與穩健性評估。

本研究選取 Weather2K-R 資料集中自 2021 年 07 月 01 日 00 時至 2021 年 08 月 31 日 21 時之觀測序列,並選取氣溫(Air Temperature)變數作為資料來源。詳細變數說明參閱附錄 A

4.1.2 MERRA-2

MERRA-2(Modern-Era Retrospective Analysis for Research and Applications, Version 2)由 NASA Goddard Earth Sciences Data and Information Services Center(GES DISC)提供,是一套全球大氣再分析資料集,透過數值天氣預報模式與以衛星觀測為主的多源觀測資料,自 1980 年起以資料同化技術重建全球大氣狀態 (GMAO, 2015) 。本研究採用其中逐時、單層、瞬時同化診斷資料產品 M2I1NXASM(Version 5.12.4)作為分析資料集。

MERRA-2 由 NASA Goddard Space Flight Center 所屬之 Distributed Active Archive Center(DAAC)典藏與管理,提供全球尺度且具一致性的再分析資料。其資料同化流程整合衛星、地面與遙測觀測資料,並於氣膠、輻射收支及水文循環等物理過程之模擬能力上較前代 MERRA 有顯著改進,因此廣泛應用於氣候分析、大氣診斷及模式評估。 M2I1NXASM 為逐時瞬時資料產品,提供高時間解析度且具代表性的氣候訊號,適合用於短期時空預測與統計特徵分析。

本研究選取 MERRA-2 之 M2I1NXASM(Version 5.12.4)資料產品,時間範圍涵蓋 2023 年 12 月 11 日 00 時至 2023 年 12 月 31 日 23 時,並以地表溫度(Surface Skin Temperature)作為主要研究變數。詳細變數說明請參閱附錄 B 。實驗區域為自全球再分析場擷取之區域子集,其範圍界定於北緯 $26.0^\circ$ 至 $48.0^\circ$、西經 $70.0^\circ$ 至 $123.0^\circ$ 之間。

4.1.3 KAUST 空間統計競賽資料集

KAUST 空間統計競賽資料集源自 2022 KAUST 空間統計競賽(KAUST Spatial Statistics Competition),該競賽旨在探討大規模空間與時空數據分析中的預測精度與計算效率議題 (Abdulah, Alamri, Nag, et al., 2022) 。該競賽資料集提供了一系列經過嚴格採樣與模擬之空間數據,專門設計用於評估複雜時空模型的插值效能,特別是在面對高維度、大規模且具備複雜空間相關性之數據情境下的表現。

本研究選取競賽資料集的配置 2a-7 進行評估 (Abdulah, Alamri, Ltaief, et al., 2022)。詳細資料集說明參閱附錄 C

4.2 實驗環境與計算資源

為確保所提出之時空整合框架在大規模資料上的訓練與推理實驗具有可行性與可靠性,本研究建立統一的運算環境,並整合多種現有軟體套件以支援模型開發、訓練及評估。

SSSD 已由原作者以 Python 實現 (AI4HealthUOL, 2023) ,可有效捕捉時間序列資料中的長程依賴結構。AFRK 則由 Wen-Ting Wang 以 R 語言實現,命名為 autoFRK (Tzeng, Huang, Wang, Nychka, et al., 2021) ,提供穩定且可擴展的空間插值能力。為整合上述功能,本研究將 autoFRK 重新實現並封裝為 Python 套件 (Tzeng, Huang, Wang, and Hsu, 2025) ,並將其演算法與 SSSD 原始碼結合,在 Python 的 PyTorch 框架下完成完整的時空整合運算流程。

所有實驗,包括模型訓練、驗證與推理,皆在 Taiwania 2 高效能運算平臺 (NCHC, 2018) 上執行。 Taiwania 2 提供高效能 GPU 運算資源,以及大容量記憶體與高速儲存系統,使本研究能夠有效處理高解析度資料與長時間序列,同時確保運算的穩定性與結果的重現性。

4.3 實驗設計

為清楚呈現 $\text{SSSD}^{\text{S4+MRTS}}$ 與 autoFRK 的設定,以下將訓練與推理的超參數整理如表 2

表 2: 模型超參數設定。

模型超參數數值
Training Configuration
Batch size40
Learning rate0.001
Only generate missingTrue
MaskingForecast
Missing $k$
SSSDS4
WaveNetInput channels1
Output channels1
Residual layers32
Residual channels64
Skip channels64
Diffusion step embeddingInput dimension128
Hidden dimension256
Output dimension256
S4Max sequence length
State dimension128
Dropout0.1
BidirectionalTrue
Layer normalizationTrue
DiffusionDiffusion steps ($T$)100
$\beta_0$0.0001
$\beta_T$0.05
MRTSUse MRTSTrue
autoFRK
MethodFast

表 3 彙整本研究各項實驗之參數設定與資料集配置,以比較不同模型之預測表現。在所有實驗中,皆自各資料集之空間範圍內隨機抽取 200 個坐標點,並進一步劃分為 160 個已觀測位置與 40 個未觀測位置。

由於各資料集之時間解析度不同,其中 MERRA-2 提供逐時(hourly)觀測資料,而 Weather2K 提供每 3 小時(3-hourly)觀測資料,因此輸入序列長度與預測區間之設定係以完整天數為單位進行規劃,而非強制於不同資料集間採用完全相同之時間比例切分。基於此原則,各資料集之時間切分方式在維持約 9:1 訓練與測試比例的前提下,盡可能使序列邊界與完整天數對齊,以符合資料原始時間解析度。

此設計可避免產生非完整天數之預測區間,並使不同取樣頻率資料集之實驗結果具有較合理之比較基礎。各資料集之空間抽樣區域如圖 2圖 4圖 6 所示。

表 3: 實驗參數與資料集設定(Weather2K / MERRA-2 / 2a-7)。

參數數值
迭代次數500
已觀測位置160
未觀測位置40
輸入時序長度448 / 456 / 90
預測時序 $k$48 / 48 / 10

https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/Stations%20of%20Weather2k%203-Hourly%20Data.png
圖 2: Weather2K 資料集之氣溫(Air Temperature)變數時空觀測資料分佈圖,呈現自 2021 年 07 月 01 日 00 時起,各觀測站位置於連續 5 個時間步之氣溫變化情形,並展示不同地理區域間異質且不規則之空間配置特性。
https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/Weather2k/sample_locations/sample_locations.png
圖 3: Weather2K 資料集中兩個代表性空間位置之觀測序列。藍色區段表示模型輸入之歷史觀測資料,紅色區段表示預測目標區間,可觀察不同地點間時間動態與變異程度之差異。
https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/Stations%20of%20MERRA-2%20Hourly%20Data.png
圖 4: MERRA-2 資料集之地表溫度(Surface Skin Temperature)變數時空資料分佈圖,呈現自 2023 年 12 月 11 日 00 時起連續 5 個時間步之再分析場資料,以及用於時空重建與預測任務之規則空間網格結構。
https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/MERRA2/sample_locations/sample_locations.png
圖 5: MERRA-2 資料集中兩個空間網格位置之地表溫度時間序列。藍色區段表示歷史觀測資料,紅色區段表示未來預測目標,可觀察受大尺度大氣動力影響之結構化時間演化特徵。
圖 6: KAUST 2a-7 模擬資料集之時空觀測資料分佈圖,呈現前 5 個時間步之模擬場值,以及受控制但具有不規則空間分佈特性的觀測位置配置。
https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/2a7/sample_locations/sample_locations.png
圖 7: KAUST 2a-7 資料集中兩個空間位置之代表性時間序列。藍色區段與紅色區段分別表示歷史觀測資料與預測目標區間,可觀察具有隨機波動且較弱決定性結構之時間動態特徵。

為評估 MRTS 空間表示之貢獻,本研究比較基準模型 $\text{SSSD}^{\text{S4}}$ 與整合 MRTS 條件資訊之 $\text{SSSD}^{\text{S4+MRTS}}$ 。藉此分析 MRTS 是否能增強模型對空間依賴結構之表徵能力,並改善推論階段對未觀測位置之時空推估表現。

此外,本研究亦納入 TFT 、 VAR 、 SVGP 與 STDK 等代表性方法作為比較基準,以涵蓋不同模型架構與建模方法,並全面評估各模型於時空預測任務中的表現差異。

5 實驗結果

本研究採用平均平方預測誤差(Mean Squared Prediction Error, MSPE)作為主要評估指標,以衡量模型之預測準確度。表 4 彙整所提出整合式時空框架之實驗結果,其中以 $\text{SSSD}^{\text{S4+MRTS}}$ 與 $\text{SSSD}^{\text{S4}}$ 訓練之模型,皆於推論階段結合 AFRK 進行空間場推估與未觀測位置之重建。

表中結果為各資料集於不同隨機種子(42 至 71)下進行 30 次獨立實驗後之平均 MSPE。此外,表中亦納入多個基準模型作為比較。由於部分基準方法僅具時間序列預測能力,因此同樣透過 AFRK 對未觀測空間位置進行推估,以確保各模型能在一致之時空預測情境下進行公平比較。

表 4: 各模型於不同資料集與預測情境下之平均 MSPE 表現。

訓練策略$\text{SSSD}^{\text{S4+MRTS}}$$\text{SSSD}^{\text{S4}}$TFTVARSVGPSTDK
空間填補AFRKAFRKAFRKAFRK
Weather2K
未觀測地點未來15.163918.572820.194223.428446.679335.1823
未觀測地點過去7.06737.07447.08777.091034.448327.0631
已觀測地點未來11.085814.565018.343419.256429.902234.9767
MERRA-2
未觀測地點未來12.907413.295716.5911389.56263.0083e4111.0262
未觀測地點過去6.47876.48116.48676.49371.6750e4111.2048
已觀測地點未來7.45127.468511.9694459.52762.8654e4107.9238
2a-7
未觀測地點未來0.89950.89870.89900.92710.97831.1991
未觀測地點過去0.90910.90930.90920.90980.99031.2128
已觀測地點未來1.02871.07521.59842.17030.91651.1642

表 4 顯示,在結合 AFRK 進行空間場推估後,$\text{SSSD}^{\text{S4}}$ 已優於多數基準模型。當進一步納入 MRTS 時,$\text{SSSD}^{\text{S4+MRTS}}$ 結合 AFRK 在多數預測情境下均能取得較 $\text{SSSD}^{\text{S4}}$ + AFRK 更低的 MSPE,且此現象於 Weather2K 與 MERRA-2 資料集上尤為明顯。此結果顯示,MRTS 所提供之空間表示資訊有助於提升模型之預測準確度,特別是在涉及未觀測位置之預測任務中,能為模型提供額外的空間資訊,進而改善時空推估表現。

表 5: 所提出模型與未整合 MRTS 版本於 30 次獨立實驗之 MSPE 平均值與標準差。

模型$\text{SSSD}^{\text{S4+MRTS}}$ + AFRK$\text{SSSD}^{\text{S4}}$ + AFRK
Weather2K
未觀測地點未來15.1639 $\pm$ 3.043418.5728 $\pm$ 3.7072
未觀測地點過去7.0673 $\pm$ 2.79487.0744 $\pm$ 2.7995
已觀測地點未來11.0858 $\pm$ 0.649414.5650 $\pm$ 1.4134
MERRA-2
未觀測地點未來12.9074 $\pm$ 2.330113.2957 $\pm$ 3.0640
未觀測地點過去6.4787 $\pm$ 1.94636.4811 $\pm$ 1.9348
已觀測地點未來7.4512 $\pm$ 0.62537.4685 $\pm$ 1.0986
2a-7
未觀測地點未來0.8995 $\pm$ 0.07750.8987 $\pm$ 0.0753
未觀測地點過去0.9091 $\pm$ 0.03860.9093 $\pm$ 0.0390
已觀測地點未來1.0287 $\pm$ 0.04131.0752 $\pm$ 0.0684

為進一步評估所提出框架之穩定性,本研究針對 $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 與 $\text{SSSD}^{\text{S4}}$ + AFRK 於 30 次獨立實驗所得之 MSPE 計算平均值與標準差,結果如表 5 所示。由於各次實驗採用不同隨機種子進行空間位置抽樣,因此標準差可反映模型在不同抽樣配置下之預測穩定性;標準差越小,表示模型對於資料抽樣變化之敏感度越低,預測表現越穩定。

表 5 的結果可觀察到, $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 在多數實驗情境下不僅具有較低的平均 MSPE,其標準差亦與 $\text{SSSD}^{\text{S4}}$ + AFRK 相當或更低,顯示 MRTS 所提供之空間表示資訊在提升預測準確度的同時,並未降低模型於不同空間抽樣配置下之穩定性。

以下所有圖表所採用之模型代號定義如下: A 為 $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK;B 為 $\text{SSSD}^{\text{S4}}$ + AFRK;C 為 TFT + AFRK;D 為 VAR + AFRK;E 為 SVGP;F 為 STDK。上述代號將在後續各圖中沿用,以便於模型間之比較。

圖 8圖 9圖 10 顯示 Weather2K 資料集於不同預測任務下,各模型 MSPE 分佈之盒狀圖。在未觀測位置之未來預測任務中, $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 的 MSPE 分佈較為集中且整體誤差較低,其中位數與下四分位數均優於 $\text{SSSD}^{\text{S4}}$ + AFRK 及其他基準模型。較低的預測誤差顯示 MRTS 所提供之空間表示資訊有助於模型學習具有空間依賴性的預測模式。此外,在圖 10 所示之已觀測位置未來預測任務中,$\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 所得到的 MSPE 亦低於 $\text{SSSD}^{\text{S4}}$ + AFRK。此結果顯示,納入 MRTS 除了有助於未觀測位置之預測外,亦能對已觀測位置之未來預測帶來額外效益。

https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/Weather2k/Unobs_and_Future/Unobs_and_Future_SSSD_S4+MRTS+AFRK_SSSD_S4+AFRK_TFT+AFRK_VAR+AFRK_SVGP_STDK_MSPE_BoxGrid1xn_nickname.png
圖 8: Weather2K 資料集於未觀測地點未來預測任務,各模型 MSPE 分佈之盒狀圖。
https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/Weather2k/Unobs_and_Past/Unobs_and_Past_SSSD_S4+MRTS+AFRK_SSSD_S4+AFRK_TFT+AFRK_VAR+AFRK_SVGP_STDK_MSPE_BoxGrid1xn_nickname.png
圖 9: Weather2K 資料集於未觀測地點過去預測任務,各模型 MSPE 分佈之盒狀圖。
https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/Weather2k/Obs_and_Future/Obs_and_Future_SSSD_S4+MRTS+AFRK_SSSD_S4+AFRK_TFT+AFRK_VAR+AFRK_SVGP_STDK_MSPE_BoxGrid1xn_nickname.png
圖 10: Weather2K 資料集於已觀測地點未來預測任務,各模型 MSPE 分佈之盒狀圖。

圖 11圖 12圖 13 所示,在 MERRA-2 資料集與不同預測任務下,各模型的 MSPE 差異更為顯著。 VAR + AFRK 、 SVGP 與 STDK 等模型在 MERRA-2 資料集上表現不佳,尤其在未觀測地點之未來預測任務中,其 MSPE 分佈呈現極端高值,顯示這些模型在處理大尺度且高維度的氣候資料時,無法有效捕捉其複雜的空間結構與動態變化,導致預測結果嚴重偏離真值。相較之下, $\text{SSSD}^{\text{S4}}$ + AFRK 雖然在 MERRA-2 上的 MSPE 較 Weather2K 高,但仍能保持相對穩定的預測表現;而加入 MRTS 之 $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 則進一步降低了 MSPE 的分佈範圍與中位數,顯示 MRTS 在增強模型對大尺度空間結構的捕捉能力方面具有顯著效果,使得模型在面對複雜氣候資料時能夠更準確地進行預測。

https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/MERRA2/Unobs_and_Future/Unobs_and_Future_SSSD_S4+MRTS+AFRK_SSSD_S4+AFRK_TFT+AFRK_VAR+AFRK_SVGP_STDK_MSPE_BoxGrid1xn_nickname.png
圖 11: MERRA-2 資料集於未觀測地點未來預測任務,各模型 MSPE 分佈之盒狀圖。
https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/MERRA2/Unobs_and_Past/Unobs_and_Past_SSSD_S4+MRTS+AFRK_SSSD_S4+AFRK_TFT+AFRK_VAR+AFRK_SVGP_STDK_MSPE_BoxGrid1xn_nickname.png
圖 12: MERRA-2 資料集於未觀測地點過去預測任務,各模型 MSPE 分佈之盒狀圖。
https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/MERRA2/Obs_and_Future/Obs_and_Future_SSSD_S4+MRTS+AFRK_SSSD_S4+AFRK_TFT+AFRK_VAR+AFRK_SVGP_STDK_MSPE_BoxGrid1xn_nickname.png
圖 13: MERRA-2 資料集於已觀測地點未來預測任務,各模型 MSPE 分佈之盒狀圖。

圖 14圖 15圖 16 呈現在 MERRA-2 資料集下, $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 、 $\text{SSSD}^{\text{S4}}$ + AFRK 與 TFT + AFRK 三種模型的 MSPE 分佈之盒狀圖。結果顯示,在各預測任務中, $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 的 MSPE 分佈明顯較其他二者更集中於較低的誤差範圍,且其四分位數優於其他兩個模型,反映 MRTS 在增強模型對空間結構的捕捉能力方面具有顯著效果。並且在已觀測地點之未來預測任務中, $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 的 MSPE 分佈更明顯較 $\text{SSSD}^{\text{S4}}$ + AFRK 與 TFT + AFRK 低,進一步證實了 MRTS 在複雜時空條件下提高精度的有效性。

https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/MERRA2/Unobs_and_Future/Unobs_and_Future_SSSD_S4+MRTS+AFRK_SSSD_S4+AFRK_TFT+AFRK_MSPE_BoxGrid1xn_nickname.png
圖 14: MERRA-2 資料集於未觀測地點未來預測任務, $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 、 $\text{SSSD}^{\text{S4}}$ + AFRK 與 TFT + AFRK 三種模型 MSPE 分佈之盒狀圖。
https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/MERRA2/Unobs_and_Past/Unobs_and_Past_SSSD_S4+MRTS+AFRK_SSSD_S4+AFRK_TFT+AFRK_MSPE_BoxGrid1xn_nickname.png
圖 15: MERRA-2 資料集於未觀測地點過去預測任務, $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 、 $\text{SSSD}^{\text{S4}}$ + AFRK 與 TFT + AFRK 三種模型 MSPE 分佈之盒狀圖。
https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/MERRA2/Obs_and_Future/Obs_and_Future_SSSD_S4+MRTS+AFRK_SSSD_S4+AFRK_TFT+AFRK_MSPE_BoxGrid1xn_nickname.png
圖 16: MERRA-2 資料集於已觀測地點未來預測任務, $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 、 $\text{SSSD}^{\text{S4}}$ + AFRK 與 TFT + AFRK 三種模型 MSPE 分佈之盒狀圖。

圖 17圖 18圖 19 顯示,在 2a-7 資料集與不同預測任務下,各模型的 MSPE 分佈之盒狀圖。在未觀測地點的預測任務中,各模型的 MSPE 分佈相當接近,顯示在空間結構較簡單且規律的模擬資料下, $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 雖然無顯著優勢,但仍能達到與其他空間建模模型相當的準確度。此外,在已觀測地點之未來預測任務中, $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 的 MSPE 分佈略高於 SVGP ,可能是因為在此資料集下,加入 MRTS 後模型對空間結構的捕捉能力提升有限,且可能引入了過度擬合的風險,使得在已觀測地點的預測表現略有下降。

https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/2a7/Unobs_and_Future/Unobs_and_Future_SSSD_S4+MRTS+AFRK_SSSD_S4+AFRK_TFT+AFRK_VAR+AFRK_SVGP_STDK_MSPE_BoxGrid1xn_nickname.png
圖 17: 2a-7 資料集於未觀測地點未來預測任務,各模型 MSPE 分佈之盒狀圖。
https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/2a7/Unobs_and_Past/Unobs_and_Past_SSSD_S4+MRTS+AFRK_SSSD_S4+AFRK_TFT+AFRK_VAR+AFRK_SVGP_STDK_MSPE_BoxGrid1xn_nickname.png
圖 18: 2a-7 資料集於未觀測地點過去預測任務,各模型 MSPE 分佈之盒狀圖。
https://josh-test-lab.github.io/posts/Spatiotemporal%20Prediction%20of%20Unknown%20Areas%20Based%20on%20Structured%20State%20Space%20Diffusion%20and%20Resolution%20Adaptive%20Fixed%20Rank%20Kriging/results/2a7/Obs_and_Future/Obs_and_Future_SSSD_S4+MRTS+AFRK_SSSD_S4+AFRK_TFT+AFRK_VAR+AFRK_SVGP_STDK_MSPE_BoxGrid1xn_nickname.png
圖 19: 2a-7 資料集於已觀測地點未來預測任務,各模型 MSPE 分佈之盒狀圖。

綜合表 4 與各資料集之 MSPE 盒狀圖可知,本研究所提出之 $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 並非在所有資料集與所有預測情境下皆取得最小的 MSPE,但其整體表現仍具高度競爭力,且在多數設定下展現出穩定且優異的預測能力。特別是在 Weather2K 與 MERRA-2 等具複雜空間結構之真實資料中, $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 在多數情境下均能維持較低誤差與較小變異,顯示 MRTS 有助於強化模型對空間依賴關係與跨區域變化的表徵能力。

然而,在 2a-7 資料集中, $\text{SSSD}^{\text{S4}}$ + AFRK 與 $\text{SSSD}^{\text{S4+MRTS}}$ + AFRK 之間的差距相當有限,顯示當資料本身具備較簡單且規律的空間結構時,MRTS 所帶來的效益較不明顯。另一方面,已觀測地點未來預測任務中,SVGP 在 2a-7 資料集上表現最佳,亦反映不同模型在不同資料分布與任務條件下,可能存在明顯的適用性差異。

整體而言,本研究方法在高維度與高變異資料中展現出較佳的穩健性與泛化能力,而在結構較簡單的資料集上,則呈現與其他模型相當之競爭表現。

6 結論

本研究針對具空間觀測不完整性之時空預測問題,提出結合 $\text{SSSD}^{\text{S4}}$ 、 MRTS 與 AFRK 的整合式框架,目的在於提升模型於未觀測位置與未來時間點之預測能力。實驗結果顯示,該方法在 Weather2K 與 MERRA-2 等真實資料集上,整體而言能有效降低預測誤差,並維持較穩定的 MSPE 分佈,顯示其對複雜時空結構具有良好的建模能力。相較之下,傳統基準模型在高維且具強烈空間異質性的資料中,較容易出現誤差放大與預測結果分佈不穩定的現象,尤其在未觀測地點未來預測任務中更為明顯。

此外,研究結果亦顯示,所提出之整合式架構並非適用於所有時空資料,其效能會隨資料特性而有所差異。當資料具有較強的空間非線性與跨尺度變化時, MRTS 能有效補強 S4 架構於空間表徵上的不足,進而提升整體預測效能。然而,對於空間結構較為平滑或變化相對單純之資料,其所帶來之效益可能有限,甚至於部分任務中未必優於其他方法。此外, AFRK 雖能有效推估未觀測位置之空間資訊,但其補值能力仍建立於既有觀測資料所呈現之空間相關性,因此當觀測資料較為稀疏、缺失比例提高,或觀測資料無法充分反映整體空間變異時,其預測能力仍可能受到影響。

另一方面,本研究與一般資料驅動式模型相同,假設訓練資料與測試資料具有相近之統計特性,因此若未來資料因極端事件、環境變遷或其他因素導致資料分佈產生明顯改變,模型效能亦可能下降,需透過重新訓練或模型調整以維持預測品質。此外,相較於原始 $\text{SSSD}^{\text{S4}}$ ,本研究方法因整合 MRTS 與 AFRK ,增加了模型建構與推論階段之計算成本,於實際應用時仍需綜合考量預測效能與計算資源之間的平衡。因此,模型之選擇仍應依資料本身之空間特性、觀測條件及實際應用需求進行評估,以兼顧預測效能與計算成本。

本研究證實深度狀態空間模型結合多尺度空間表徵與推論階段之空間補值機制,能有效提升具空間觀測不完整性之時空預測準確度與穩健性。未來可進一步探討模型於多變數、不同空間尺度、更高缺失比例及更多元環境資料上的適用性,並發展兼具預測效能與計算效率之時空建模方法,以提升模型於環境監測、氣象預報及地球科學分析等實際應用情境中的泛化能力與實用價值。

參考文獻

附錄

A Weather2K 資料集

Weather2K-R 資料集所包含之變數如下。本研究使用其地表溫度(Air Temperature)變數作為主要分析對象。

表 6: Weather2K 變數表。
變數名稱縮寫單位
Latitudelatdegrees east
Longitudelondegrees north
Altitudealtm
Air PressureaphPa
Air Temperaturet°C
Maximum Temperaturemxt°C
Minimum Temperaturemnt°C
Relative Humidityrh%
Precipitation in 3hp3mm
Wind Directionwddegrees
Wind Speedwsm s$^{-1}$
Maximum Wind Directionmwddegrees
Maximum Wind Speedmwsm s$^{-1}$

B MERRA-2 資料集

MERRA-2 資料集所包含之變數如下。本研究使用其 Surface Skin Temperature 變數作為主要分析對象。

表 7: MERRA-2 變數表。
變數名稱縮寫單位
Longitudelondegrees east
Latitudelatdegrees north
Timetimeminutes since 2024-06-01 00:00:00
2-Meter Air Temperaturet2mK
Total Precipitable Liquid Watertqlkg m$^{-2}$
Total Column Odd Oxygentoxkg m$^{-2}$
2-Meter Eastward Windu2mm s$^{-1}$
Surface PressurepsPa
Tropopause Temperature Using Blended TROPP EstimatetroptK
Northward Wind at 50 Metersv50mm s$^{-1}$
Zero Plane Displacement Heightdisphm
Total Column Ozoneto3Dobsons
Surface Skin TemperaturetsK
10-Meter Air Temperaturet10mK
Tropopause Pressure Based on Thermal EstimatetropptPa
Total Precipitable Ice Watertqikg m$^{-2}$
Sea Level PressureslpPa
Tropopause Pressure Based on Blended EstimatetroppbPa
Total Precipitable Water Vaportqvkg m$^{-2}$
2-Meter Northward Windv2mm s$^{-1}$
Tropopause Specific Humidity Using Blended TROPP Estimatetropqkg kg$^{-1}$
10-Meter Northward Windv10mm s$^{-1}$
Eastward Wind at 50 Metersu50mm s$^{-1}$
10-Meter Eastward Windu10mm s$^{-1}$
2-Meter Specific Humidityqv2mkg kg$^{-1}$
Tropopause Pressure Based on EPV EstimatetroppvPa
10-Meter Specific Humidityqv10mkg kg$^{-1}$

C 第二屆大型資料集空間統計競賽資料集

本資料集源自 2022 KAUST Competition on Spatial Statistics for Large Datasets 所公布之大規模模擬資料,設計目的在於於統一環境下評估各類空間統計方法在預測與參數推估任務中的表現。此系列資料以 ExaGeoStat 高效能統計運算框架生成,透過可重現之高斯過程(Gaussian process, GP)模擬,提供標準化且具比較性的實驗條件 (Abdulah, Alamri, Nag, et al., 2022; Abdulah, Alamri, Ltaief, et al., 2022) 。

競賽中 Sub-competition 2a 與 2b 採用非分離且平穩之 GP 模型,其共變異數結構依據 Gneiting (2002) 所提出之非分離型共變異數形式進行設定。對於任兩空間位置 $s \in [0, 1]^2$ 與時間 $t \in \mathbb{R}$,共變異數函數表示如下: \begin{equation} C(\boldsymbol{h},u ; \boldsymbol{\theta}) = \frac{\sigma^{2}}{a_t |u|^{2\alpha} + 1} \, M_{\nu}\!\left( \frac{\lVert \boldsymbol{h} \rVert / a_s}{\left(a_t |u|^{2\alpha} + 1\right)^{\beta/2}} \right), \end{equation} 其中 $\boldsymbol{h}$ 為空間距離、$u$ 為時間差; $\sigma^{2} > 0$ 為變異數, $\nu>0$ 與 $\alpha \in [0, 1]$ 為平滑參數; $a_s, a_t > 0$ 為空間與時間尺度; $\beta \in (0, 1]$ 描述空間與時間間之交互強度; $M_{\nu}(\cdot)$ 為 Matérn 型相關函數。

依據競賽規格, Sub-competition 2a 與 2b 共產生十八組資料,涵蓋不同尺度配置(弱、中度、強)、不同空間點數(1K 與 10K)、不同時間長度(100 與 1000 時距),並針對預測任務提供三類缺失設定:隨機移除空間位置(RS)、隨機移除空間與時間位置(RST)、以及最後十個時間點系統性缺失(T10)。所有參數組合詳列於 Abdulah, Alamri, Nag, et al. (2022) 所列之 Table 1 。 而 2a-7 即為其中一項具特定空間點數、時間長度與模型參數組合之資料版本。

D 其他模型設定

本章說明實驗中所使用之其他基準模型的主要超參數設定,包含 TFT、 VAR 、 SVGP 與 STDK 等模型,以供參考與重現實驗結果。

程式設計部分請參考 https://github.com/Josh-test-lab/SSSDS4-AFRK

D.1 TFT

本研究使用 Darts 套件實作 Temporal Fusion Transformer (TFT) 模型,表 8 為其主要超參數設定。

表 8: TFT 模型超參數設定。

超參數數值
Batch size40
Learning rate0.001
Epochs500
Add relative indexTrue
Input chunk length80
Hidden size32
LSTM layers1
Dropout0.1
Output size1

D.2 VAR

本研究使用 statsmodels 套件實作 Vector Autoregression (VAR) 模型,表 9 為其主要超參數設定。

表 9: VAR 模型超參數設定。

超參數數值
Maximum lag order240
Information criterionaic
Channel dimension1

D.3 SVGP

本研究使用 https://pypi.org/project/gpytorch/ 之 GPyTorch 套件實作 Sparse Variational Gaussian Process (SVGP) 模型,表 10 為其主要超參數設定。

表 10: SVGP 模型超參數設定。

超參數數值
Learning rate0.01
Epochs500
OptimizerAdam
LikelihoodGaussian
DeviceCUDA / CPU
Kernel typeRBF Kernel
ARD dimensions3
Mean functionConstant Mean
Variational distributionCholesky Variational
Variational strategyLearn inducing points
Inducing points1024
Loss functionVariational ELBO
Likelihood noiseLearned Gaussian noise

D.4 STDK

本研究使用 https://pypi.org/project/da-stdk/ 之 da-stdk 套件實作 Spatio-temporal DeepKriging (STDK) 模型,表 11 為其主要超參數設定。

表 11: STDK 模型超參數設定。

超參數數值
DeviceCUDA
Split methodRandom
Train ratio0.8
Calibration ratio0.2
Calibration split methodRandom
Observation methodSite-wise
Observation patternUniform
Observation ratio0.1
Spatial intensity10.0
Provider splitFalse
Learning rate0.001
Basis learning rate ratio0.05
Epochs500
Batch size40
Weight decay5e-4
Patience50
Gradient clipping10.0
SchedulerCosine
Warmup epochs10
Basis unfreeze epoch10
Basis LR ramp-up epochs10
Number of workers0
Hidden dimensions128, 128, 128
Spatial basis functionWendland
Spatial centers25, 81, 121
Temporal centers10, 15, 45
Spatial learnableFalse
Spatial init methodUniform
Delta reparameterizationFalse
Layer normalizationTrue
Gradient dampingTrue
Damping threshold0.0
Damping strength5.0
Domain penalty weight0.01
Movement penalty weight0.0
Sparsity typeNone
Sparsity L1 weight0.0
Sparsity group weight0.0
Spatial sparsity appliedTrue
Temporal sparsity appliedFalse
Sparsity threshold ratio0.01
Non-crossing lambda0.0
Non-crossing L1True
Input covariates0
Regression typeMulti-quantile
Quantile levels0.5
Conformal modeBoth
Conformal alpha0.1
Conformal center sourceTrained
Conformal cluster min N30
Save quantile predictionsTrue
Target normalizationTrue

參見

1
2
3
4
5
6
7
8
@mastersthesis{Hsu2026SSSD_AFRK,
    author      = {HSU, YAO-CHIH},
    title       = {Spatiotemporal Prediction of Unknown Areas Based on Structured State Space Diffusion and Resolution Adaptive Fixed Rank Kriging},
    school      = {National Dong Hwa University},
    type        = {Master's thesis},
    year        = {2026},
    url         = {https://ndltd.ncl.edu.tw/cgi-bin/gs32/gsweb.cgi?o=dnclcdr&s=id=%22114NDHU0507005%22.&searchmode=basic},
}