• Data assimilation (DA)
• NWP | DA | Kalman Filter
• NWP | DA | 4D-Var

---

Data assimilation (DA)

Data assimilation (DA) (數據同化，或稱資料同化) is a technique by which numerical model data and observations are combined to obtain an analysis that best represents the state of the atmospheric phenomena of interest

History

• In 1901 Cleveland Abbe founded the United States Weather Bureau. He suggested that the atmosphere followed the principles of thermodynamics and hydrodynamics
• In 1904, Vilhelm Bjerknes proposed a two-step procedure formodel-based weather forecasting. First, a analysis step of data assimilation to generate initial conditions, then a forecasting step solving the initial value problem.
• In 1922, Lewis Fry Richardson carried out the first attempt to perform the weather forecast numerically.
• In 1950, a team of the American meteorologists Jule Charney, Philip Thompson, Larry Gates, and Norwegian meteorologist Ragnar Fj ̈ortoft and the applied mathematician John von Neumann, succeeded in the first numerical weather forecast using the ENIAC digital computer.
• In September 1954, Carl-Gustav Rossby’s group at the Swedish Meteorological and Hydrological Institute produced the first operational forecast (i.e. routine predictions for practical use) based on the barotropic equation. Operational numerical weather prediction in the United States began in 1955 under the Joint Numerical Weather Prediction Unit (JNWPU), a joint project by the U.S. Air Force, Navy, and Weather Bureau.
• In 1959, Karl-Heinz Hinkelmann produced the first reasonable primitive equation forecast, 37 years after Richardson’s failed attempt. Hinkelmann did so by removing high-frequency noise from the numerical model during initialization.
• In 1966, West Germany and the United States began producing operational forecasts based on primitive-equation models, followed by the United Kingdom in 1972, and Australia in 1977.

Polynomial Fitting of Panofsky

Panofsky, R. A. (1949). Objective weather-map analysis. Journal of Atmospheric Sciences, 6(6), 386-392.

\[ \begin{equation} \begin{split} p(x,y) = \sum_{i,j} a_{ij} x^i y^j, \qquad i+j \leq 3 \end{split} \end{equation} \tag{2} \]

• A field of 10 observations can be fitted accurately by a third-degree polynomial, with no smoothing of data.
• The amount of smoothing will depend on the nature of the variable analyzed.

Successive Correction Method (SCM) of Cressman

Cressman G P, 1959. An operational objective analysis system[J]. Mon Wea Rev, 87: 367-374.

• Purpose:
  • SCM is designed to improve the accuracy of gridded meteorological data by integrating observations from various sources, such as weather stations, into a coherent grid. This method is particularly useful in situations where observational data is sparse or unevenly distributed.
• Basic Principle:
  • The method operates on the principle of iteratively correcting a background field (initial guess) using available observational data. Each iteration refines the grid values based on nearby observations, effectively merging the observational data with the background field.
• Weighting Function:
  • Cressman's SCM employs a specific weighting function that determines how much influence each observation has on the grid point being corrected. This function typically decreases with distance from the observation point, allowing closer observations to have a greater impact on the grid value than those further away.
• Computational Efficiency:
  • One of the advantages of SCM is its computational simplicity and efficiency, making it suitable for real-time applications in weather forecasting and analysis. It does not require complex background error covariance matrices or assumptions about error distributions, which simplifies its implementation compared to other assimilation methods like Optimal Interpolation (OI) or Variational methods.

Inverse Problem

DA consists of using the actual result of some measurements to infer the values of the parameters that characterize the system.

Random Variable

Bayesian

Mathematics

Algorithms

The specific data assimilation algorithms can be divided into two main categories:

1. Continuous Data Assimilation Algorithms: This includes Three-Dimensional Variational (3DVar) and Four-Dimensional Variational (4DVar) methods, which optimize the model trajectory within the assimilation window to make the model results closer to the observed data.
   1. \[x^a = x^b+BH^T(R+HBH^T)^{-1}(y-Hx^b) \qquad (1)\]
   2. \[J(x)=\frac12(x-x^\mathrm{b})B^{-1}(x-x^\mathrm{b})^\mathrm{T}+\frac12(y-Hx^\mathrm{b})R^{-1}(y-Hx^\mathrm{b}) \qquad (2)\]
2. Sequential Data Assimilation Algorithms: Represented by the Ensemble Kalman Filter (EnKF), these methods dynamically update the model state to integrate observed data in real time.
   1. Sequential? 在後續的模式積分過程中，我們不光需要積分模式，也需要將分析誤差協方差矩陣 "積分" (誤差協方差矩陣也進行"更新")，以獲得下次觀測時刻的背景誤差協方差。
   2. Extended Kalman filter (EKF) \[ \begin{aligned} &\boldsymbol{x}_{k} =\boldsymbol{Mx}_{k-1}, \qquad (3) \\ &\boldsymbol{y}_k =\boldsymbol{Hx}_k, \qquad \quad (4)\\ \hline &\boldsymbol{K}=\boldsymbol{P}_k^\mathrm{f}\boldsymbol{H}^\mathrm{T}\left[\boldsymbol{H}\boldsymbol{P}_k^\mathrm{f}\boldsymbol{H}^\mathrm{T}+\boldsymbol{R}\right]^{-1} \qquad (5)\\ &\boldsymbol{x}_k^\mathrm{a}=\boldsymbol{x}_k^\mathrm{f}+\boldsymbol{K}\left[\boldsymbol{y}_k-\boldsymbol{H}\boldsymbol{x}_k^\mathrm{f}\right] \qquad \quad (6) \\ &\boldsymbol{P}_k^\mathrm{a}=(\mathbf{I}-\boldsymbol{K}\boldsymbol{H})\boldsymbol{P}_k^\mathrm{f} \qquad \qquad\qquad (7) \\ \hline &\boldsymbol{x}_{k+1}^\mathrm{f}=\boldsymbol{M}\boldsymbol{x}_{k}^\mathrm{a} \qquad\qquad (8) \\ &\boldsymbol{P}_{k+1}^\mathrm{f}=\boldsymbol{M}\boldsymbol{P}_{k}^\mathrm{a}\boldsymbol{M}^\mathrm{T}+\boldsymbol{Q} \quad (9)\\ \hline \end{aligned} \]
   3. Take Lorenz63 as example, Extended Kalman filter (EKF) \[ \begin{aligned} &\frac{dx}{dt}=\sigma(y-x) \qquad (10)\\ &\frac{dy}{dt}=\rho x-y-xz \qquad (11)\\ &\frac{dz}{dt}=xy-\beta z \qquad (12)\\ \hline &\boldsymbol{x}_{k} =\mathcal{M}(\boldsymbol{x}_{k-1}), \qquad (13) \\ &\boldsymbol{y}_k =h(\boldsymbol{x}_k), \qquad \quad (14) \\ \hline &x_{k}^{\mathrm{f}}=\mathcal{M}\big(x_{k-1}^{\mathrm{a}}\big) \qquad (15)\\ &M=\frac{\partial\mathcal{M}}{\partial x}\bigg|_{x_{k-1}^{\mathrm{a}}} \qquad (16)\\ &P_{k}^{\mathrm{f}}=MP_{k-1}^{\mathrm{a}}\boldsymbol{M}^{\mathrm{T}}+\boldsymbol{Q}\qquad (17)\\ \hline &H={\frac{\partial h}{\partial x}}\bigg|_{x_{k}^{\mathrm{f}}} \qquad (18)\\ &K=P_{k}^{\mathrm{f}}H^{\mathrm{T}}\left(HP_{k}^{\mathrm{f}}H^{\mathrm{T}}+R\right)^{-1} \qquad (19)\\ &x_{k}^{\mathrm{a}}=x_{k}^{\mathrm{f}}+K\left[y_{k}-h\left(x_{k}^{\mathrm{f}}\right)\right]\qquad (20) \\ &P_{k}^{\mathrm{a}}=(\mathbf{I}-\mathbf{K}\mathbf{H})P_{k}^{\mathrm{f}}\qquad (21)\\ \hline &\boldsymbol{M}=\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{x}}=\begin{bmatrix}\frac{\partial f_1}{\partial x}&\frac{\partial f_1}{\partial y}&\frac{\partial f_1}{\partial z}\\\frac{\partial f_2}{\partial x}&\frac{\partial f_2}{\partial y}&\frac{\partial f_2}{\partial z}\\\frac{\partial f_3}{\partial x}&\frac{\partial f_3}{\partial y}&\frac{\partial f_3}{\partial z}\end{bmatrix}=\begin{bmatrix}-\sigma&\sigma&0\\\rho-z&-1&-x\\y&x&-\beta\end{bmatrix} \end{aligned} \]
   4. 集合卡尔曼滤波器（Ensemble Kalman Filter, EnKF)
      1. EKF的不足
      2. 有沒有特別的辦法讓 \(B\) 矩陣可以演化，又不需要切線模式 \(M\)，最好又能不佔用龐大的內容來保存 \(B\) 呢？
      3. EnKF 的基本想法是使用機率空間內的點表示狀態變數的一種可能。在這個想法下，EnKF最具革命性的工作是使用集合來估計預報誤差協方差矩陣。
      4. 在EnKF中，流依賴性不是體現在更新協方差裡面，而是更新狀態集合，然後利用不斷變化的狀態集合來用隨時計算協方差。
      5. Lorenz63模式無法提現其中的計算成本。但如前面說的，實際模式中，切線模式的積分是平方倍的計算量，變數數多的時候會變成天文數值。而在EnKF中，其實 N只要有數十個或數百個就行了。還是比較有限的計算量。兩者對比，還是EnKF可行性更高。此外，EnKF不涉及切線模式的截斷誤差，對於非線性模式也沒有太大障礙。
      6. EnKF的不足: 需要透過添加擾動來產生集合

The computational cost of variational assimilation lies in the multiple iterations of the optimization algorithm and the potential for non-convergence. In contrast, the cost of sequential assimilation comes from the need to use ensembles to estimate uncertainty, resulting in a significant computational load from integrating multiple model runs.

Disturbation method

• 典型的擾動方式是奇異向量（Singular Vector）擾動，即先計算初始場中的快速增長模態（奇異向量），然後在這些方向上添加擾動。這種方法能夠有效捕捉初始場中不確定性成長最快的模式，進而提高集合預報的離散度。
• 另一個熱門的擾動方式為隨機擾動參數化傾向方案（Stochastically Perturbed Parameterization Tendencies，SPPT），透過在物理參數化方案的傾向項中引入隨機擾動，來表徵模式不確定性。相較於完全隨機的擾動方式，上述方法產生的擾動場通常具有空間和時間上的自相關性，並確保一定的物理合理性。

對於集合預報的結果，集合成員的平均值被視為確定性預報，主要用於評估集合預報的平均表現。確定性預報指標包括平均絕對誤差（MAE），均方根誤差（RMSE）等。而機率預報指標用於評估集合預報的機率分佈特性，衡量預報的可靠性和不確定性，包括 連續排位機率評分（CRPS），可靠性（Reliability），集合離散度（Spread）等。其中集合離散度通常由集合成員的標準差表示，用於衡量集合成員與集合平均值之間的差異。在集合預報中，集合離散度與集合平均的均方根誤差（RMSE）之間的關係被用來衡量集合預報系統的可靠性。理論上，一個可靠的集合預報系統的集合離散度應與集合平均預報誤差大致相當。這種特性在集合預報領域通常被稱為一致性（Consistency）。

The inremental method (增量法)

• Courtier et al. 1994, A strategy for operational implementation of 4D-Var, using an incremental approach

增量法由 Courtier 提出，此方法主要應用在三維或四維變分同化中，特別是在四維變分同化（4DVar）時。在把4DVar業務化的過程中，要考慮同化的時效性，而同化中 \(x\) 向量的維度可達 \(10^7\) 到 \(10^9\) 。4DVar 中伴隨模式後向積分的計算量是傳統數值模式（NWP）計算量的兩倍，從而使4DVar進行一次迭代的計算量是數值模式計算量的三倍。再加上數值模式本身的非線性特點，使得同化中損失函數的收斂速度很慢。增量法理論上是用精確度換時間的方法（a cost-benefit trade-off）。增量法的使用要結合變分求解時所使用的內循環 (inner loop) 和 外循環 (outer loop) 過程。

References

1. 朱国富, 2015. 数值天气预报中分析同化基本方法的历史发展脉络和评述. 气象, 41(8): 986-996. DOI: 10.7519/j.issn.1000-0526.2015.08.008.
2. 朱国富, 2015. 理解大气资料同化的根本概念[J]. 气象, 41(4): 456-463. DOI:10.7519/j.issn.1000-0526.2015.04.008
3. 2014: Data Assimilation in Numerical Weather Prediction models
4. 2018: 中央氣象局全球資料同化系統
5. 数据同化的算法演变：3DVar、EKF、EnKF (推薦)| 海洋学院，沈浙奇，2024年6月
6. 资料同化中的增量法（The inremental method）