数据同化相关反问题的平滑性及其他超参数估计

数据同化相关反问题的平滑性及其他超参数估计

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 迭代¹。后者被称为 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