EOF分解原理及Python实现
原创EOF分解原理及Python实现
原创
大三上课第一次听到EOF可以把时空变化场分解为空间函数和时间函数时，我很迷惑，为什么可以把空间和时间分离开，那么空间随时间的偏导数怎么处理？正交说明时间和空间函数没有直接的依赖关系，而空间在随时间变化，空间与时间是彼此关联的，为什么可以正交呢？为什么可以降维呢？总之，就是不理解这样操作的可行性和物理意义。故从原理入手，去理解该分解。
分解原理概述
EOF分解的核心思想是将一个随时间变化的变量场（如气象要素场）分解为两部分：一部分是不随时间变化的空间函数（V），它概括了要素场的空间分布特点；另一部分是只依赖时间变化的时间函数（Z），它表征了典型场随时间的演变特征。
这种分解是通过数学上的线性代数方法实现的，具体涉及计算数据的协方差矩阵，并求解其特征值和特征向量。特征向量对应的是空间函数，而特征向量对应的时间权重系数则构成了时间函数。
降维并不是说把时间和空间在原有空间上完全割裂开，而是通过变换到新空间，从而最大限度的保留原有信息，保留下来的信息并不直接就表示原来的物理实际意义，但仍与原信息相对应。
对正交的理解
在数学上，两个向量（或函数）正交意味着它们的点积（或内积）为零。在EOF分解中，空间函数（特征向量）是正交的，这意味着它们之间没有线性相关性。
EOF分解得到的空间函数和时间函数是相互正交的，这意味着它们之间没有直接的依赖关系。然而，这并不意味着空间分布不随时间变化，而是这些变化已经被有效编码在时间函数中了。
EOF分解的目的是找到数据集中最重要的变化模式，这些模式在空间和时间上都是独立的。通过正交性，我们可以确保每个空间函数都捕捉了数据中的一个独特变化模式，而不会与其他模式重叠。这使得我们能够更容易地理解和解释数据中的变化。
EOF分解中的正交性是一种数学工具，用于提取数据中的独立变化模式。虽然它可能给人一种空间和时间完全独立的印象，但实际上它们是通过时间函数相互关联的。这种分解方法使我们能够更容易地理解和解释数据中的时空变化。EOF分解只是一种统计方法，其得到的分解结果本身没有任何的物理意义。
EOF与PCA
EOF类似于PCA，本质是降维，通过PCA/EOF，提取出信号的主要变化特征，数学原理上相似，它们都涉及协方差矩阵的计算和特征向量的提取。EOF在某些文献中也被称为经验正交函数（EOF）主成分分析（PCA）。
EOF主要应用于大气和海洋科学，主要处理具有空间和时间变化的数据；PCA应用范围更加广泛，主要用于数据降维和特征提取，可以处理高维数据。
PCA主成分分析的数学原理
以下内容参考：CodingLabs - PCA的数学原理
主成分分析PCA以及特征值和特征向量的意义_主成分分析特征值-CSDN博客
PCA通过线性变换将原始数据变换为一组各维度线性无关的表示，可用于提取数据的主要特征分量，常用于高维数据的降维。
降维意味着信息的丢失，因为数据本身存在相关性，我们可以想办法在降维的同时将信息的损失尽量降低。
怎样降维不丢失信息呢？例如：
接下来，从原理出发，重新“发明”一遍PCA。自己从原理上认识一遍该过程比直接套用算法理解更加深刻。
向量与基变换
多个同类数据 → 一组向量
向量内积（即点乘）
基
代数上，我们用线段的终点坐标表示向量，如上图的向量(4,2)，但这样表示的前提是：以x轴和y轴上正方向长度为1的向量为标准。
便引出“基”的概念。
上图中(1,0)和(0,1)叫做二维空间中的一组基。
要描述向量，首先要确定基，不过一般默认(1,0)和(0,1)为一组基。事实上，任意两个线性无关的二维向量都可以是一组基，一般取模为1。
基变换
①原基向量(1,0)和(0,1)
②构建新的基向量坐标空间，例如：
③如何将原本在(1,0)和(0,1)为基向量的坐标空间的向量变换到新的空间呢？即该向量的坐标在新空间中是多少？
两个矩阵相乘的意义是将右边矩阵中的每一列列向量变换到左边矩阵中每一行行向量为基所表示的空间中去。一个矩阵可以表示一种线性变换。
降维投影
原数据分布在二维空间，如果要降维，如何变换到一维空间，还仍然保留原有信息呢？
即在二维平面中选择一个方向，将所有数据都投影到这个方向所在直线上，用投影值表示原始记录。
那么如何选择这个方向（基）才能保留最多的原始信息呢？一种直观的看法是：希望投影后的投影值尽可能分散。
上图如果向x轴做投影，则会导致多点重叠，严重数据丢失。
那么，如果用数学的方式去描述这种“投影值尽量分散”，即分散程度呢？
（1）方差：衡量离散程度
把二维降为一维，找到方差最大的方向。数据离散性越大，表明在所投影的维度上具有越高的区分度。
（2）协方差：表示两个变量之间的线性相关性
考虑三维降到二维问题。同样，先找到一个方向使投影后方差最大，这样就完成了第一个方向的选择，继而寻找第二个投影方向。
如果第二投影方向还是单纯只选择方差最大的方向，很明显，这个方向将与第一个方向“几乎重合在一起”，这样的维度是没用的，因此，还应该有其他约束条件。
从直观上说，让两个字段尽可能表示更多的原始信息，我们是不希望它们之间存在（线性）相关性的，因为相关性意味着两个字段不是完全独立，必然存在重复表示的信息。
故约束条件为协方差为0，选择第二个基时只能在与第一个基正交的方向上选择。
这就解释了我开头的疑惑，正交并不是为了把时间和空间分离开而一股脑的不管时间和空间之间的关系，恰恰相反，而是尽可能的保留了原数据的最多信息，防止数据重叠。
Cov(a,b)=\frac{1}{m}\sum_{i=1}^m{a_ib_i}，每个字段均值已处理为0。
至此，我们得到了降维问题的优化目标：将一组N维向量降为K维（K大于0，小于N），其目标是选择K个单位（模为1）正交基，使得原始数据变换到这组基上后，各字段两两间协方差为0，而字段的方差则尽可能大（在正交的约束下，取最大的K个方差）。
（3）协方差矩阵：便于将上述两点的原理转化为具体算法
协方差矩阵对角线上的元素是两个字段的方差，其它元素是两个字段的协方差，两者被统一到了一个矩阵。
（4）协方差矩阵对角化：距离发明“PCA”仅一步之遥
即除对角线外的其它元素化为0（即协方差化为0），并且在对角线上将元素按大小从上到下排列（即对方差排序）。
特征值和特征向量：
方阵的特征值对应的特征向量就是理想中想取得的变换后的的坐标轴，特征值就等于数据在变换后的坐标上对应维度的方差。
某一特征值除以所有特征值的和的值就为：该特征向量的方差贡献率（方差贡献率代表了该维度下蕴含的信息量的比例）。
EOF分解方法
以下内容参考：经验正交分解EOF详解及案例_eof分解-CSDN博客
气象要素的时空数据集F大多是3维，包括2维的空间场V以及1维的时间序列Z。
分解之后，根据方差贡献率的大小对空间模态进行排序，方差贡献率最大的称为第一空间模态，方差贡献率第二大的称为第二空间模态。
所有空间模态(或时间系数)所对应的方差贡献率加起来等于 1，所以，当方差贡献率最高的几个空间模态(或时间系数)加起来得到的累计方差贡献率较大时(比如达到80%)，我们就认为这几个模态可以表示原始数据中所包含的主要特征。
此刻我有一个新的疑问，找到了这样一个可以保留原数据大部分信息的空间，如何还原出原数据背后隐含的物理意义呢？
答：需要深入理解EOFs和PCs的物理含义，分析它们的时间和空间变化。
EOFs：代表了数据集中的主要空间变化模式，通常按照它们对数据方差的贡献程度进行排序，EOF1代表数据中的主要空间变化模式，EOF2代表次要的空间变化模式，依此类推。
EOFs可以用于理解数据的空间分布特征，例如气候模式的空间分布、海洋温度场的空间结构等。
PCs：是与EOFs相对应的时间序列，描述了EOFs随时间的变化情况。
PCs可以帮助我们理解数据中的时间变化特征，例如气候模式的季节变化、年际变化等。
EOFs代表数据中的空间变化模式，而PCs则描述了这些空间模式随时间的变化情况。
Python实现EOF分解
EOFs代表数据的空间变化模式，而PCs是与EOFs相对应的时间系数，代表这些模式在时间上的变化。
import xarray as xr
import datetime
import numpy as np
from eofs.xarray import Eof
import matplotlib.pyplot as plt
import proplot as pplt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter, LatitudeLocator, LongitudeLocator)
import cmapsimport warningswarnings.filterwarnings('ignore')
data = xr.open_dataset("C:/Users/lenovo/Desktop/sst1991-2020.nc")
print(data)
sst=data['sst']
lon, lat = data['longitude'], data['latitude']
lon2d, lat2d = np.meshgrid(lon, lat)
months = sst['time'].dt.month
jan_indices = (months == 1)
jan_data = sst.isel(time=jan_indices)
feb_indices = (months == 2)
feb_data = sst.isel(time=feb_indices)
mar_indices = (months == 3)
mar_data = sst.isel(time=mar_indices)
print(jan_data)
def mapart(ax):
# 指定投影为经纬度投影，并指定中心经度为180°
projection = ccrs.PlateCarree(central_longitude=180)
ax.set_extent([0, 360, 20, 90], crs=ccrs.PlateCarree(central_longitude=180))
# 设置经纬度标签
ax.set_xticks(np.arange(-180,180+10,30), crs=projection)
ax.set_yticks([20,30,40,50,60,70,80,90], crs=projection)
# 给标签添加对应的N,S,E,W
lon_formatter = LongitudeFormatter()
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
# 添加海岸线
ax.coastlines(color='k', lw=0.5)
# 添加陆地
ax.add_feature(cfeature.LAND, facecolor='white')
coslat = np.cos(np.deg2rad(data.coords['latitude'].values))
wgts = np.sqrt(coslat)[..., np.newaxis]
solver = Eof(jan_data, weights=wgts, center=True) #center=True自动距平化
eofs = solver.eofs(neofs=10, eofscaling=2) #取前10个
pcs = solver.pcs(npcs=10, pcscaling=1)
variance = solver.varianceFraction(neigs=10)
print(eofs,pcs,variance)
plt.figure()
plt.plot(np.arange(1, 11), variance[:10] * 100, 'o-')
plt.xlabel('EOF modes')
plt.ylabel('Variance explained (%)')
plt.title('Variance explained by EOF modes')
plt.grid(True)# 在每个点上添加文本标签
for i, v in enumerate(variance.values[:10] * 100):
plt.text(i + 1, v + 0.9, f'{v:.2f}', fontsize=8, color='black')plt.show()
# 计算特征值和NorthTest的结果
eigenvalues = solver.eigenvalues(neigs=10)
errors = solver.northTest(neigs=10)
# 绘制Error Bar图
plt.figure()
plt.errorbar(np.arange(1, 11), eigenvalues, yerr=errors, fmt='o-', capsize=5)
plt.xlabel('EOF modes')
plt.ylabel('Eigenvalues')
plt.title('NorthTest(with error bars)')
plt.grid(True)plt.show()
var = (variance*100).values
projection = ccrs.PlateCarree(central_longitude=180) # 指定投影为经纬度投影，中心经纬度为180°fig, ax = pplt.subplots(proj=projection,proj_kw={'central_longitude':0})#cmap = pplt.Colormap('ColdHot',cut=0)p = ax.contourf(lon2d,lat2d,eofs[1], levels=np.linspace(-1, 1, 21), cmap=cmaps.BlueWhiteOrangeRed,transform=ccrs.PlateCarree(),extend='both')ax.colorbar(p,loc='bottom')ax.set_title('EOF2 (%s' % (round(var[1], 2))+"%)", loc='left')mapart(ax)plt.show()
years = np.arange(1991, 2021)
fig, ax = pplt.subplots()
b = ax.bar(years, pcs[:, 1], color='r')# 对时间系数值小于0的柱子设置为蓝色
for bar, height in zip(b, pcs[:, 1]):
if height < 0:
bar.set(color='blue')
# 为第二个子图添加标题
ax.set_title('PC2' % (round(var[1], 2)), loc='left')
plt.show()
'''
jan_data_mean = jan_data - jan_data.mean(dim='time')
def mapart(ax):
# 指定投影为经纬度投影，并指定中心经度为180°
projection = ccrs.PlateCarree(central_longitude=180)
ax.set_extent([0, 360, 20, 90], crs=ccrs.PlateCarree(central_longitude=180))
# 设置经纬度标签
ax.set_xticks(np.arange(-180,180+10,30), crs=projection)
ax.set_yticks([20,30,40,50,60,70,80,90], crs=projection)
# 给标签添加对应的N,S,E,W
lon_formatter = LongitudeFormatter()
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
# 添加海岸线
ax.coastlines(color='k', lw=0.5)
# 添加陆地
ax.add_feature(cfeature.LAND, facecolor='white')
projection = ccrs.PlateCarree(central_longitude=180)
fig, ax = pplt.subplots(proj=projection,proj_kw={'central_longitude':0})
cmap = pplt.Colormap('ColdHot',cut=0)
p = ax.contourf(lon2d,lat2d,jan_data_mean[0,:,:]/100,cmap=cmap)
ax.colorbar(p,loc='bottom',label='Mean Sea Level Pressure (1991.1)')
mapart(ax)
plt.show()
'''
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
腾讯云开发者
