POD (Proper Orthogonal Decomposition)

POD（Proper Orthogonal Decomposition，本征正交分解）

数据预处理

考虑现在有一个二维速度场 $q=(u,v)$ ，然后通过仿真或者实验获得这个二维速度场的数据

$$q_{ij}=\left( u_{ij},v_{ij} \right) =\left( u\left( \xi _i,\eta _j \right) ,v\left( \xi _i,\eta _j \right) \right)$$

其中 $( \xi _i,\eta _j )$ 是空间离散网格点坐标，$i=1,\cdot \cdot \cdot ,m_{\xi}$ ，$j=1,\cdot \cdot \cdot ,m_{\eta}$ ，下图为该二维速度场均匀网格点上的离散速度数据示意图

假设在 $t=t_k$ 时刻，得到了一组数据并将其堆叠成一个列向量 $x(t_k)$ ，并把这个大小为 $2m_{\xi}m_{\eta}\times1$ 的列向量设为 $m\times1$ 的列向量

$$\left[ \begin{array}{c} \begin{matrix} u_{11}& \cdots& u_{m_{\xi}1}\\ \vdots& \ddots& \vdots\\ u_{1m_{\eta}}& \vdots& u_{m_{\xi}m_{\eta}}\\ \end{matrix}\\ \begin{matrix} v_{11}& \cdots& v_{m_{\xi}1}\\ \vdots& \ddots& \vdots\\ v_{1m_{\eta}}& \cdots& v_{m_{\xi}m_{\eta}}\\ \end{matrix}\\ \end{array} \right] _{t=t_k}\underset{unstack}{\overset{stack}{\begin{array}{c} \rightharpoonup\\ \leftharpoondown\\ \end{array}}}\ \ \left[ \begin{array}{c} \left( \begin{array}{c} u_{11}\\ \vdots\\ u_{1m_{\eta}}\\ \end{array} \right)\\ \vdots\\ \left( \begin{array}{c} u_{n_{\xi}1}\\ \vdots\\ u_{m_{\xi}m_{\eta}}\\ \end{array} \right)\\ \left( \begin{array}{c} v_{11}\\ \vdots\\ v_{1m_{\eta}}\\ \end{array} \right)\\ \vdots\\ \left( \begin{array}{c} v_{m_{\xi}1}\\ \vdots\\ v_{m_{\xi}m_{\eta}}\\ \end{array} \right)\\ \end{array} \right] _{t=t_k}\equiv x\left( t_k \right) \in R^m$$

现收集到一系列时间快照（snapshot）($\begin{matrix} t_1& t_2& \cdots& t_n\\ \end{matrix}$) 的流场数据，并把每一时刻对应的数据整理成上述的列向量，从而构造出一组数据矩阵 $x$

$$x=[ \begin{matrix} x\text{(}t_1\text{)}& x\text{(}t_2\text{)}& \cdots& x\text{(}t_n\text{)}\\ \end{matrix} ] \in R^{m\times n}$$

通常情况下，我们只会考虑流场的速度脉动，所以矩阵 $X$ 需要减去一个时间均值，即

$$\bar{x}\left( t \right) =\frac{1}{n}\underset{i=1}{\overset{n}{\varSigma}}x\left( t_i \right)$$

$$X=x-\left[ \begin{matrix} \bar{x}\left( t \right)& \bar{x}\left( t \right)& \cdots& \bar{x}\left( t \right)\\ \end{matrix} \right]$$

分解算法

POD分解把上述矩阵分解为空间分量和时间分量

$$X_{m\times n}=\underset{r}{\varSigma}a_r\left( t \right) \phi _r\left( \xi \right)$$

其中 $a(t)$ 表示时间系数，$\phi(\xi)$ 表示空间模态，$\xi$ 表示空间向量，m 是参数的个数，n 是时间快照的个数

也可以用矩阵表示

$$X_{m\times n}=\varPhi \varSigma \varPsi ^T$$

其中 $\varPhi$ 为空间模态，$\varSigma$ 为幅值，$\varPsi$ 为 temporal dynamics 描述了时间上的演变

---

POD 分解的目的就是找到一组最优的基向量能够最好地描述这个数据矩阵，这个组基向量就是 $\phi$

并且 $\phi$ 就是求协方差矩阵 $XX^T\in R^{m\times m}$ 的特征向量

准确地讲，协方差矩阵是 $\frac{1}{n}XX^T$ 或者是 $\frac{1}{n-1}XX^T$ ，但前面的系数并不会影响特征向量的求解，系数会作用在特征值上

$$XX^T\phi=\lambda \phi$$

且 $\varPhi$ 是单位正交矩阵，即 $\phi$ 模长唯1，两两正交

特征值 $\lambda_i$ 表示每个时刻的观测量在特征向量 $\phi_i$ 上的投影点到原点距离的平方和，再开根号

The eigenvalues $\lambda_i$ convey how well each eigenvector $\phi_i$ captures the original data in the L2 sense (scaled by n)

L2 是向量的2范数，是向量元素平方和的平方根，用于衡量向量的长度或大小

这里可以看一下第一个视频

对于一个速度矩阵，特征值 $\lambda$ 可以衡量动能，$\lambda =\sigma ^2$ 奇异值 $\sigma$​ 对应着幅值，且是一个速度矩阵，幅值的平方就可以衡量能量

通常为了数据降维，只保留能量较大的几个模态，剩余的皆可舍去，这里假设前 3 个模态能量占主要，即

$$\frac{\underset{i=1}{\overset{3}{\varSigma}}\lambda _i}{\underset{i=1}{\overset{n}{\varSigma}}\lambda _i}=\frac{\underset{i=1}{\overset{3}{\varSigma}}\sigma _{i}^{2}}{\underset{i=1}{\overset{n}{\varSigma}}\sigma _{i}^{2}}=99\%\approx 1$$

只需要留下前三个模态即可

Spatial (Classical) POD Method 和 Method of Snapshots

这两种方法没有太大区别，本质上就是比较 $XX^T$ 和 $X^TX$ 两个矩阵的大小

如果 $m<n$ 则矩阵 $XX^T\in R^{m\times m}$ 维度小，就使用 $XX^T$ 矩阵求模态和特征值 $\lambda$ ，是 Spatial (Classical) POD Method

如果 $m>n$ 则矩阵 $X^TX\in R^{n\times n}$ 维度小，就使用 $X^TX$ 矩阵求模态和特征值 $\lambda$ ，是Method of Snapshots

这就是求矩阵 $X$ 奇异值分解的步骤，因此可以直接使用奇异值分解

对 POD 分解的每一个矩阵的讨论

$$X_{m\times n}=\varPhi \varSigma \varPsi ^T=U\varSigma V^T$$

这里为方便表示，使用 $U$ 表示模态，$V$ 为 temporal dynamics

$$X=\left[ \begin{matrix} u_1& u_2& \cdots& u_m\\ \end{matrix} \right] \left[ \begin{matrix} \sigma _1& 0& \cdots& 0& \cdots& 0\\ 0& \sigma _2& \cdots& 0& \cdots& 0\\ \vdots& \vdots& \ddots& \vdots& \ddots& \vdots\\ 0& 0& \cdots& \sigma _r& \cdots& 0\\ \vdots& \vdots& \ddots& \vdots& \ddots& \vdots\\ 0& 0& \cdots& 0& \cdots& 0\\ \end{matrix} \right] \left[ \begin{array}{c} v_1^T\\ v_2^T\\ \vdots\\ v_n^T\\ \end{array} \right]$$

$$X=\left[ \begin{matrix} X\text{(}t_1\text{)}& X\text{(}t_2\text{)}& \cdots& X\text{(}t_n\text{)}\\ \end{matrix} \right] =\sigma _1u_1v_{1}^{T}+\sigma _2u_2v_{2}^{T}+\cdot \cdot \cdot +\sigma _ru_rv_{r}^{T}=\underset{i=1}{\overset{r}{\varSigma}}\sigma _iu_iv_{i}^{T}=\underset{i=1}{\overset{r}{\varSigma}}u_i\left( \sigma _iv_{i}^{T} \right)$$

为后续解释方便，将矩阵 $V^T$ 展开写

$$V=\left[ \begin{matrix} v_1& v_2& \cdots& v_n\\ \end{matrix} \right] =\left[ \begin{matrix} V_{11}& V_{12}& \cdots& V_{1n}\\ V_{21}& V_{22}& \cdots& V_{2n}\\ \vdots& \vdots& \ddots& \vdots\\ V_{n1}& V_{n2}& \cdots& V_{nn}\\ \end{matrix} \right] \ \ \ \ V^T=\left[ \begin{array}{c} v_{1}^{T}\\ v_{2}^{T}\\ \vdots\\ v_{n}^{T}\\ \end{array} \right] =\left[ \begin{matrix} V_{11}& V_{21}& \cdots& V_{n1}\\ V_{12}& V_{22}& \cdots& V_{n2}\\ \vdots& \vdots& \ddots& \vdots\\ V_{1n}& V_{2n}& \cdots& V_{nn}\\ \end{matrix} \right]$$

可以得出

$$X\left( t_1 \right) =\sigma _1u_1V_{11}+\sigma _2u_2V_{12}+\cdots +\sigma _ru_rV_{1r}$$

$$X\left( t_2 \right) =\sigma _1u_1V_{21}+\sigma _2u_2V_{22}+\cdots +\sigma _ru_rV_{2r}$$

$$\vdots$$

$$X\left( t_i \right) =\sigma _1u_1V_{i1}+\sigma _2u_2V_{i2}+\cdots +\sigma _ru_rV_{ir}$$

从上式可以看出，$V$ temporal dynamics 矩阵的第 i 行代表在第 i 个时间快照每个模态的“权重”，第 j 列代表模态 $u_j$ 在每个时间快照上“权重”的变化（时间上的形状）， $\sigma$​​ 表示相应的幅值

$\sigma _iv_{i}^{T}$ 对应时间系数 $a_r(t)$ ，时间系数矩阵是 $\varSigma V^T$ ， $u_i$ 对应空间模态 $\phi(\xi)$ 

关于 temporal dynamics 矩阵的应用

假设做实验时获取数据的时间步长是1s，为了更精细描述流场，需要更短的时间步长0.5s，又不想重新做实验，则可以对数据进行插值。

插值有两种选择：

1. 直接对原来的数据矩阵 $X_{m\times n}$ 在时间方向（每一行）上进行插值
2. 对时间矩阵 $\varSigma V^T\ _{m\times n}$ 每一行进行插（或者对 $V_{n\times n}$ 的每一列/ $V^T_{n\times n}$ 的每一行进行插值）

插值完成之后，原矩阵的列数变成对应的时间点数，如 $X_{m\times n}$ --> $X_{m\times 2n}$ / $X_{m\times 3n}$​

两种方法其实没有太大区别，如果保留全部 POD 模态，基于 POD 的时间插值与直接数学插值本质上是等价的；POD 插值真正有意义的地方，在于只保留前几个主模态时，能够降低矩阵维度，减少一些计算量，具体可参考文献【5】

References

1. Principal Component Analysis
2. Strang, G. (2023). Introduction to Linear Algebra (6th ed.). Cambridge: CUP.
3. 2017-Modal Analysis of Fluid Flows: An Overview
4. Peter J. Schmid, Chapter Six - Data-driven and operator-based tools for the analysis of turbulent flows, Advanced Approaches in Turbulence
5. Note on the POD-based time interpolation from successive PIV images