數據同化相關反問題的平滑性及其他超參數估計

數據同化相關反問題的平滑性及其他超參數估計

Smoothness and other hyperparameter estimation for inverse problems related to data assimilation

https://arxiv.org/pdf/2602.18328

摘要

本文研究由偏微分方程和隨機偏微分方程控制的動力系統在數據同化中出現的貝葉斯反問題。我們聯合推斷時空相關場以及先驗密度和似然密度的靜態參數。特別關注控制先驗光滑性和正則性的超參數，這對確保問題適定性、塑造后驗結構以及確定預測不確定性至關重要。通常該參數被假定為先驗已知且固定，但在本文中我們將采用分層貝葉斯框架，將光滑性及其他超參數視為未知并賦予超先驗。后驗推斷使用適用于高維的Metropolis-within-Gibbs抽樣進行，其中超參數估計幾乎不增加計算負擔。該方法在Navier-Stokes方程和隨機對流擴散方程的反問題中得到驗證，涵蓋稀疏和密集觀測情形，并使用具有不同協方差結構的高斯先驗。數值結果表明，聯合估計光滑性可顯著降低因光滑性設定錯誤而導致的不確定性量化和參數估計誤差，其性能可與已知真實光滑性的場景相媲美。

1 引言
在氣象學、海洋學與大氣科學、水文學、地球物理學等領域，需要估計各種時空相關的動態量并用于預測；參見[6, 20, 5, 4]的綜述。這些量通常采用流體動力學中的物理模型進行建模，表述為偏微分方程，最近也開始使用隨機偏微分方程，描述諸如速度、溫度、氣壓、風或洋流的湍流、鹽度等大氣量的演變[62, 3, 35]。過去幾十年中，從遙感衛星或原位傳感器獲取的時空數據集規模日益增大，例如[24]。衛星能夠生成地球物理變量的高分辨率圖像，如平流層臭氧、海面高度風等，這些圖像被用于從季節周期估計到長期趨勢識別的多種任務[54, 53]。數據同化提供了一個框架，將這些觀測數據與動力學模型相結合，以準確估計和預測這些狀態。數值天氣預報的核心目標之一——預報，可以通過推斷狀態的初始條件并向前傳播其控制動力學來實現。由于這些模型軌跡會迅速偏離真實狀態，我們需引入觀測數據來修正這種偏差，詳見[46, 4, 6]。除了被稱為動力學變量的狀態外，還存在一些靜態參數，它們決定了模型方程及其屬性，如光滑性、噪聲幅度、旋轉角度等。這些參數也需要從數據中估計，本文將考慮對動態變量和靜態變量進行聯合貝葉斯推斷[1, 7]。

數據同化問題自然地被表述為貝葉斯反問題：這在氣象學、大氣科學與海洋學等領域已很常見，例如參見[6, 20]。受這些應用啟發，我們采用貝葉斯觀點，從含噪觀測中推斷一個定義在合適可分向量空間(V, ∥·∥V)上的向量場v = (v(t, x); t ∈ [0, T ], x ∈ X)。過去十年中，得益于基于模擬的算法的發展和強大計算能力的普及，貝葉斯反問題的應用顯著增長[25]。在此框架下，我們旨在通過抽象貝葉斯法則從數據y學習后驗分布μ：

1.2 主要貢獻

• 我們針對偏微分方程和隨機偏微分方程的反問題，提出了一種高效的 MwG 算法。對于狀態變量，我們使用 pCN 算法；對于超參數，我們采用吉布斯采樣或吉布斯 MCMC 核，與文獻[7]類似。我們強調，考慮超參數所帶來的計算開銷非常小。

• 我們特別關注光滑度參數對狀態推斷和預測的影響，這在先前的文獻中常常被忽略。為此，我們證明了，與依賴固定光滑度假設的標準方法相比，我們的方法在估計和預測的不確定性量化方面均有改進。我們還探討了光滑度對其他參數估計的影響。

本文的結構如下：第2節介紹了問題表述、使用的符號以及將要研究的兩個不同問題。第3節描述了我們在高維設置下考慮的Metropolis-within-Gibbs算法及其在每個案例研究中的具體實現。隨后，第4節展示了數值結果。最后，第5節進行總結并討論未來工作。附錄包含補充材料和額外的數值實驗。

2.2 案例研究 1：推斷確定性 Navier-Stokes 方程的初始條件

2.3 案例研究 2：基于隨機對流擴散方程的數據同化

現在，我們研究一個模型，其中狀態動力學遵循上的線性隨機偏微分方程。

3 用于濾波和反問題的蒙特卡洛方法 3.1 馬爾可夫鏈蒙特卡洛方法

MCMC 是一種迭代過程，通過模擬一個遍歷的、時間齊次的、以 μ 為不變分布的馬爾可夫鏈的長時間運行，來從目標分布 μ 中采樣。經過一個預熱期后，鏈的樣本構成 μ 的經驗近似。有多種馬爾可夫鏈方法可用于從給定的目標后驗分布 μ 中采樣，例如吉布斯采樣器或 Metropolis-Hastings 算法；概述可參見文獻[59]。

我們的設置需要在希爾伯特空間上定義的分布中進行采樣，并且后驗分布位于無限維狀態空間中。對于關于 v 的反問題，我們不容易利用條件依賴性來開發吉布斯類型的采樣器以更新系數塊并提高效率。在實踐中，我們使用高維離散化，但我們的方法仍然需要對網格細化具有穩健性。在給定 ? , θ , y 的條件下，直接從 v 的滿條件分布中采樣是不可能的，因此我們使用適用于高維的 MCMC 方法。對于具有高斯先驗的目標分布，隨機游走 Metropolis 算法不是有效的選擇，因為隨著維度的增加，它們的譜隙會必然崩潰[22]。文獻[59]提出了對先驗分布保持不變的自回歸提議分布。這些提議分布最近被重新命名為預條件 Crank–Nicolson，并在文獻[10, 55]中引入用于高維反問題的推斷。文獻[22]證實，與隨機游走方法相比，pCN 的譜隙以及因此方法的混合速度不會隨著維度增長而崩潰，并且克服了標準 MCMC 在高維中遇到的一些混合緩慢問題。我們將 ? , θ
固定為當前估計值，并在 v v的滿條件分布上運行 pCN 迭代。偽代碼如算法 1 所示。

可以使用直接從上述每個條件分布中采樣的方法，或者如果無法直接采樣，則可以采用少量針對上述每個密度的 Metropolis 迭代1。后者被稱為 MwG。當可用時，推薦使用基于共軛關系的直接采樣以加速采樣過程。

4 數值結果

4.1 案例研究1：Navier-Stokes逆問題

平穩態

我們現在在圖 3 中繪制估計渦度的方差，以全面評估算法的不確定性量化。從數值范圍可以看出，方差受 α 的選擇影響很大。只有在使用 MwG 時，我們才能恢復出已知真實 α 時的方差量級。當 α 被高估和低估時，不確定性分別被低估和高估。實踐者常常為了改進點估計而犧牲不確定性量化，反之亦然。在我們的框架內，可以通過調整 α 的先驗來實現這一點。但請注意，我們將使用高粘度，這會降低系統的能量；這樣一來，這種平穩狀態使得我們的方差估計變得不那么關鍵，并且所有算法在少數時間步后的預測能力都相似。下面我們將在粘度較低、觀測更頻繁且信息量更大的“混沌”狀態下重復此實驗，以研究其中的差異。

混沌狀態

我們在參數估計中獲得了與平穩狀態類似的好結果，見圖9。圖10所示的渦度點估計也是如此。這些圖可以在附錄C.2節中找到。在圖4中，我們展示了與之前圖3類似的結果。在未觀測位置，渦度估計方差的量級與之前相同，但觀測引入了空間相關性。似然在信息量大的觀測位置附近是尖峭的，這降低了這些位置速度后驗的方差，也降低了后驗對先驗的依賴。對于較高的 α 值，這種空間相關性更為明顯，因為相應的先驗方差被平滑了。因此，較高的 α 導致先驗方差的信息量減少，進一步增加了后驗對似然的依賴。這就是為什么方差像在平穩設置中一樣，僅在未觀測位置受到 α 的影響。從觀測時刻速度估計的軌跡圖中也可以看出這一點。由于混沌狀態具有高能量，軌跡圖顯示方差的差異不會隨時間消失，進而影響估計的預測能力。 α = 3 導致未觀測位置的方差被低估：在觀測之后很短的時間步內，置信區間就無法捕捉真實值。在觀測位置，預測表現良好，這意味著似然足以修正錯誤設定的先驗。至于 α = 1.5
，方差被高估，導致在所有時間點的不確定性量化都很差。在觀測位置，其預測直到最后一次觀測前表現良好，但其預測能力隨著方差的爆炸而迅速下降。那么，使用 α = 1.5
的 pCN 在不確定性量化方面表現不佳，尤其是對于未觀測狀態和預測。這表明，即使在似然影響最大的地方，先驗仍然影響著算法的預測能力。最后，MwG 重現了與假設正確 α α值的算法相似的性能，即使在預測中也表現同樣出色。

5 結論

原文鏈接：https://arxiv.org/pdf/2602.18328