📚 课程目录第4章基于主成分分析的光谱降维与线性判别分析分类

4. 基于主成分分析的光谱降维与线性判别分析分类

任务说明

本案例以一组模拟的近红外光谱分类数据为输入，完成以下闭环任务：使用主成分分析（PCA）将高维光谱特征降至低维子空间，再以线性判别分析（LDA）对降维后的特征进行分类建模，最终给出训练集与测试集的分类准确率，并输出PCA-LDA决策边界可视化所需的投影数据。代码需完整可运行，不依赖外部光谱文件。

实现方案

整体流程分为四步：数据生成（模拟三类光谱，每类含不同吸收峰位置与噪声）、PCA降维（保留前两个主成分）、LDA分类（在二维主成分空间上训练线性判别分类器）、评估与可视化（计算准确率并绘制投影散点图）。PCA负责无监督地压缩光谱波长维度，LDA则利用类别标签寻找最佳判别方向，二者串联构成典型的“先降维后分类”光谱分析管线。

代码实现

import numpy as np
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
import matplotlib.pyplot as plt

# ---------- 1. 生成模拟光谱数据 ----------
np.random.seed(42)
n_samples_per_class = 50
n_wavelengths = 200
wavelengths = np.linspace(900, 1700, n_wavelengths)  # 近红外波段

def generate_spectra(class_id, n_samples):
    """生成三类模拟光谱：每类在特定波段有高斯吸收峰"""
    spectra = np.zeros((n_samples, n_wavelengths))
    for i in range(n_samples):
        base = 0.9 + 0.1 * np.random.randn(n_wavelengths)
        if class_id == 0:
            peak = 0.5 * np.exp(-((wavelengths - 1100) ** 2) / (2 * 30 ** 2))
        elif class_id == 1:
            peak = 0.5 * np.exp(-((wavelengths - 1300) ** 2) / (2 * 30 ** 2))
        else:
            peak = 0.5 * np.exp(-((wavelengths - 1500) ** 2) / (2 * 30 ** 2))
        spectra[i] = base + peak + 0.02 * np.random.randn(n_wavelengths)
    return spectra

X = np.vstack([generate_spectra(c, n_samples_per_class) for c in range(3)])
y = np.hstack([np.full(n_samples_per_class, c) for c in range(3)])

# ---------- 2. 划分训练/测试集 ----------
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, stratify=y, random_state=42
)

# ---------- 3. PCA降维至2维 ----------
pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train)
X_test_pca = pca.transform(X_test)

# ---------- 4. LDA分类 ----------
lda = LinearDiscriminantAnalysis()
lda.fit(X_train_pca, y_train)
y_pred = lda.predict(X_test_pca)

# ---------- 5. 评估 ----------
train_acc = accuracy_score(y_train, lda.predict(X_train_pca))
test_acc = accuracy_score(y_test, y_pred)
print(f"训练集准确率: {train_acc:.4f}")
print(f"测试集准确率: {test_acc:.4f}")
print(f"PCA解释方差比: {pca.explained_variance_ratio_}")

# ---------- 6. 可视化 ----------
plt.figure(figsize=(8, 6))
colors = ['#e74c3c', '#3498db', '#2ecc71']
for cls in range(3):
    mask = y_test == cls
    plt.scatter(X_test_pca[mask, 0], X_test_pca[mask, 1],
                c=colors[cls], label=f'Class {cls}', edgecolors='k', s=60)

# 绘制LDA决策边界
x_min, x_max = X_test_pca[:, 0].min() - 1, X_test_pca[:, 0].max() + 1
y_min, y_max = X_test_pca[:, 1].min() - 1, X_test_pca[:, 1].max() + 1
xx, yy = np.meshgrid(np.linspace(x_min, x_max, 200),
                     np.linspace(y_min, y_max, 200))
Z = lda.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
plt.contourf(xx, yy, Z, alpha=0.2, colors=['#e74c3c', '#3498db', '#2ecc71'])
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title('PCA-LDA Classification on Test Set')
plt.legend()
plt.grid(alpha=0.3)
plt.savefig('pca_lda_result.png', dpi=150)
plt.close()
print("可视化图像已保存为 pca_lda_result.png")

运行与验证

依赖安装：需Python 3.8+，并安装numpy、scikit-learn、matplotlib。执行pip install numpy scikit-learn matplotlib即可。

运行方式：将上述代码保存为pca_lda_spectra.py，在终端执行python pca_lda_spectra.py。

预期输出：由于随机种子固定为42，每次运行结果一致。训练集准确率约为0.99～1.00，测试集准确率约为0.93～1.00（具体值取决于数据划分，但应显著高于随机猜测的0.33）。PCA前两个主成分的累计解释方差比应大于0.95，表明二维子空间保留了绝大部分光谱信息。程序同时生成pca_lda_result.png，图中应显示三个类别在PC1-PC2平面上呈明显分离的簇，且LDA决策边界将平面划分为三个区域。

判定标准：若测试集准确率≥0.90且图像中三类样本无大面积重叠，则案例实现成功。若准确率偏低，可检查PCA降维后类别是否可分（观察散点图），或适当增加模拟光谱中吸收峰的强度差异。