Python | 垂直模态分解（phython垂直输出）

垂直模态分解概述

七八月的金陵像一口大铁锅，学校里面施工的锤子敲得叮当作响，连空气都带着钢筋水泥的味道。暑期的太阳倒像疯了似的，把地面烤得能煎蛋。开完题歇了几周，脑子里的科研也差点长了锈。乘着暑假把之前的烂摊子收拾收拾，往前推一推，不然就跟学校的脚手架一样，永远停在那儿。顺便写写代码，免得这夏天白热了一场。

★

EOF分析是气象学中提取主要变化模态和识别主导空间模式的重要方法。

1. EOF 基本思想

EOF 分析（即 PCA）用于将高维数据投影到按方差大小排序的正交模态上，提取主要变率。

常见的使用是针对三维数据 [time, lat, lon] ，可以得到空间模态EOF以及对应的时间序列PC。比如说冬季海温异常的信号。

1.1 数据矩阵

1.2 协方差矩阵

EOF的本质通过计算协方差矩阵的特征值和特征向量

1.3 特征分解

1.4 主成分（PC）

表示第 k EOF 模态的时间系数。

1.5 重构

2. 数据预处理：距平

3. SVD

对距平矩阵 做 SVD：

4. 垂直模态分解

以上思路同样的逻辑，那么针对 [time, level, lat, lon] 数据，就可以用来提取典型垂直廓线。

具体来说，如果你输入的变量是垂直高度或者非绝热加热时，

同理，计算PC：

实际数据的处理和EOF类似，或者了解聚类分析的话，只不过将垂直层数作为第一维。从数据维度上来说，即：

[time, level, lat, lon] -> [level, time x lat x lon]

这样应该会很清晰

4.1 重排维度

其中 m = level 数，N = time × lat × lon。

4.2 计算SVD

这里numpy和scipy都有现成的函数，直接调用即可

4.3 方差

实操

输入数据格式

这里使用的数据是ERA5的垂直速度，monthly mean 资料，单位是 Pa/s ，选择热带区域15°S-15°N，高度为：1000-100 hPa

'w' (time: 516, level: 0, latitude: 31, longitude: 360)> Size: 0B
[0 values with dtype=float64]
Coordinates:
  * longitude  (longitude) float32 1kB -180.0-179.0-178.0 ... 178.0179.0
  * latitude   (latitude) float32 124B -15.0-14.0-13.0 ... 13.014.015.0
  * level      (level) int32 0B 
  * time       (time) datetime64[ns] 4kB 1979-01-011979-02-01 ... 2021-12-01
Attributes:
    units:          Pa s**-1
    long_name:      Vertical velocity
    standard_name:  lagrangian_tendency_of_air_pressure

代码逻辑

w(time, level, lat, lon)
↓ subset (lat=-15~15, level=1000–100 hPa, time=1979–2021)
↓ 去除气候态季节循环，得到 w_anom
↓ 转置 -> (nlev, nt, nlat, nlon)
↓ 展平 (time × lat × lon) -> (nlev, N_samples)
↓ 去缺测列 -> (nlev, N_valid)
↓ SVD -> U(nlev, nlev), s(nlev), Vh(nlev, N_valid)
↓ EOF = U (列为垂直模态)
↓ PC = U.T @ X = Σ Vh
↓ eigval = s^2 / N_valid
↓ 计算方差贡献率

EOF模态的气象学解释

1. EOF1 : 通常代表最主要（方差贡献最大）的垂直运动模态，可能对应深对流或大尺度垂直环流
2. EOF2 : 常表现为垂直结构的第二主要变化，可能反映不同层次间的耦合
3. 高阶模态 : 反映更复杂的垂直结构和局地过程

核心代码

def compute_eof(self, data: np.ndarray) -> Dict:
        """
        Compute EOF analysis using SVD.
        
        Parameters:
        -----------
        data : np.ndarray
            Input data matrix (variables x observations)
            
        Returns:
        --------
        Dict
            Dictionary containing EOF results
        """
        # Handle missing values
        valid_mask = ~np.any(np.isnan(data), axis=0)
        data_valid = data[:, valid_mask]
        
        # Perform SVD
        u, s, v = np.linalg.svd(data_valid, full_matrices=False)
        
        # Calculate EOF patterns and PC time series
        eof_patterns = u.T  # EOF patterns
        pc_series = np.dot(eof_patterns, data_valid)  # Principal components
        
        # Calculate eigenvalues and explained variance
        nt = data_valid.shape[1]
        eigenvalues = s**2 / nt
        explained_variance = eigenvalues / np.sum(eigenvalues) * 100
        
        # Estimate degrees of freedom for North test
        L = 1# one-lag autocorrelation
        phi_L, phi_0, dof = self._estimate_dof(data_valid, L)
        eigenvalue_errors = explained_variance * np.sqrt(2 / dof)
        
        return {
            'eof_patterns': eof_patterns,
            'pc_series': pc_series,
            'eigenvalues': eigenvalues,
            'explained_variance': explained_variance,
            'eigenvalue_errors': eigenvalue_errors,
            'degrees_of_freedom': dof,
            'valid_mask': valid_mask,
            'phi_0': phi_0,
            'phi_L': phi_L
        }

展示结果

掩膜陆地，只有海洋

对于热带海洋（15°S-15°N），前两个EOF模态方程贡献率可以达到84.2%，基本上可以代表了热带海洋上的主要垂直结构。 一般来说，

第二个 EOF 中垂直运动的符号在 650 hPa 左右切换。

掩膜海洋，只有陆地

全球区域，无掩膜

拓展read

此外，前人研究表明，热带对流云中降水产生的效率很大程度上取决于沿着垂直方向发展的高强度对流。例如，强降雨事件往往发生在深对流区域，那里的云顶通常可以增长到对流层顶，例如西太平洋暖池。另一方面，小雨事件通常与经常发生在信风逆温地区的浅层对流有关（例如，亚热带海洋东部地区）。

热带太平洋的两种典型结构模态：

通过选择两个区域方框，可以使用上述SVD分解得到其主要的垂直廓线，可以清晰的发现两个方框中截然不同的垂直结构：

参考

★

Bui, H. X., J.-Y. Yu, and C. Chou, 2016: Impacts of Vertical Structure of Large-Scale Vertical Motion in Tropical Climate: Moist Static Energy Framework. Journal of the Atmospheric Sciences, 73, 4427–4437, https://doi.org/10.1175/JAS-D-16-0031.1.

Back, L. E., Z. Hansen, and Z. Handlos, 2017: Estimating Vertical Motion Profile Top-Heaviness: Reanalysis Compared to Satellite-Based Observations and Stratiform Rain Fraction. J. Atmos. Sci., 74, 855–864, https://doi.org/10.1175/JAS-D-16-0062.1.