如何用PCA分析并解读300组随海拔变化的温度剖面数据？

我有300组温度-海拔数据，每组是20×300的规模，各样本特征存在差异。这些随海拔变化的温度剖面，让我在以下任务中遇到挑战：

• 数据可视化
• 识别变化幅度最大的样本
• 确定偏离标称温度剖面最小的样本编号
• 归类具有相似模式的温度剖面

我尝试用**主成分分析（PCA）**将数据维度降至2，但不确定如何解读分析结果，也不清楚这个方向是否正确。

复现数据的示例代码 #

import numpy as np
import matplotlib.pyplot as plt


"""  Create 300 Sample data for Altitude and Temperature"""
# Define the number of rows and columns
num_rows, num_columns = 20, 300

# Create an array of altitudes ranging from 100 km to 500 km
altitudes = np.linspace(100, 500, num_rows)

# Repeat the altitude values across all columns
altitude_array = np.tile(altitudes, (num_columns, 1)).T

# Parameters for generating temperature profiles
altitude_midpoint = (altitudes.min() + altitudes.max()) / 2
temperature_amplitude, temperature_frequency = 250, 0.02

# Create a sine wave to represent temperatures with altitude
temperatures = temperature_amplitude * np.sin(temperature_frequency * (altitudes - altitude_midpoint)) + 750

# Create random noise for slight temperature variations
max_variation = 50
noise = np.random.uniform(-max_variation, max_variation, (num_rows, num_columns))

# Add the noise to the temperatures to create temperature profiles
temperature_array = temperatures[:, np.newaxis] + noise


""" Principal Component Analysis"""

from sklearn.decomposition import PCA

run_number = np.arange(1, 301).reshape(300, 1)

pca = PCA(2)  
principal_components = pca.fit_transform(temperature_array.T)
print(temperature_array.shape)
print(principal_components.shape)


# Create a figure with two subplots
fig, axs = plt.subplots(1, 3, figsize=(12, 5))

# First subplot
# Plot temperature vs. altitude for all 300 profiles
for i in range(num_columns):
    sc1 = axs[0].plot(temperature_array[:, i], altitude_array[:, i], label=f'Profile {i+1}', alpha=0.8)

# Plot the nominal temperature profile with marked lines
nominal_temperature = temperature_amplitude * np.sin(temperature_frequency * (altitudes - altitude_midpoint)) + 750
axs[0].plot(nominal_temperature, altitudes, 'k', linewidth=2, label='Nominal Profile', linestyle='--')

# Set plot labels and title
axs[0].set_ylabel('Altitude (km)')
axs[0].set_xlabel('Temperature (K)')
axs[0].set_title('Temperature vs. Altitude Profiles')

# Second subplot
sc2 = axs[1].scatter(principal_components[:, 0], principal_components[:, 1],
                     c=run_number.T, edgecolor='none', alpha=0.5,
                     cmap=plt.cm.get_cmap('Accent'))
axs[1].set_xlabel('component 1')
axs[1].set_ylabel('component 2')
fig.colorbar(sc2, ax=axs[1])

# Third subplot
sc3 = axs[2].plot(np.cumsum(pca.explained_variance_ratio_))
axs[2].set_xlabel('number of components')
axs[2].set_ylabel('cumulative explained variance')


plt.tight_layout()

生成的样本图表 #

PCA方向是否正确？ #

你的方向完全合理：PCA是处理这类高维度剖面数据的理想工具，能有效提取数据中的核心变异模式，完美匹配你需要的可视化、异常识别、聚类等任务。从你的累计方差图可以推断，前2个主成分应该已经解释了绝大多数数据方差（因为你的数据是在正弦标称曲线上叠加随机噪声，主成分会对应噪声的主要变异方向）。

如何解读PCA结果？ #

1. 累计方差图（第三个子图） #

该图展示了前N个主成分能解释的总方差比例。如果前2个主成分的累计方差比例很高（比如超过90%），说明将数据降维到2维几乎保留了所有关键信息，后续分析的可靠性有保障。

2. 主成分散点图（第二个子图） #

每个点对应一个温度剖面样本：

• 相似模式归类：散点图中距离近的点，对应温度剖面的模式更相似；反之则差异更大。你可以直接在这个2维空间里用聚类算法（如K-means）划分相似组别。
• 异常样本识别：远离集群的孤立点，就是变化幅度最大的样本，对应原始数据中偏离标称曲线最显著的剖面。
• 偏离最小样本定位：散点图中心区域的点对应数据的平均模式，和标称曲线最接近，结合颜色映射的样本编号，就能直接找到这些样本。

3. 主成分的物理意义 #

通过查看pca.components_可以理解每个主成分代表的变异模式：

print(pca.components_)

每个主成分是长度为20的向量，对应不同海拔高度的温度变异权重。比如第一个主成分如果在所有海拔上符号一致，说明它代表整体温度的偏移；如果符号交替，说明它代表剖面形状的扭曲（比如正弦曲线的幅度变化）。

针对任务的具体操作建议 #

1. 确定偏离标称最小的样本：计算每个样本与标称曲线的均方误差（MSE），取MSE最小的样本编号：

nominal = temperature_amplitude * np.sin(temperature_frequency * (altitudes - altitude_midpoint)) + 750
mse = np.mean((temperature_array.T - nominal)**2, axis=1)
min_mse_idx = np.argmin(mse) + 1  # 样本编号从1开始
print(f"偏离标称最小的样本编号：{min_mse_idx}")

2. 识别变化幅度最大的样本：取MSE最大的样本，或者在PCA散点图中找到距离中心最远的点对应的编号。

3. 归类相似模式：在PCA的2维结果上运行K-means聚类：

from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=3, random_state=42)
clusters = kmeans.fit_predict(principal_components)

# 可视化聚类结果
plt.scatter(principal_components[:,0], principal_components[:,1], c=clusters, cmap='Accent', alpha=0.5)
plt.xlabel('Component 1')
plt.ylabel('Component 2')
plt.title('Clustered Temperature Profiles')
plt.show()

内容的提问来源于stack exchange，提问作者Mark