高光谱遥感处理

引言

高光谱遥感在可见光到短波红外范围内以几纳米到十几纳米的间隔采样，一次成像能获得上百个连续波段。多光谱只有几个到十几个宽波段，只能区分大类地物；高光谱的连续光谱曲线携带了矿物成分、植被生化参数、水质参数等精细信息，能识别具体矿物种类、区分作物品种、探测伪装目标，也能在单个像元内分辨出多种地物的混合比例。

工程上的第一道门槛是数据量。一景 1000 乘 1000 像元、200 波段、16 位的高光谱影像有 4 亿字节，相当于同等空间分辨率的多光谱影像的几十倍。第二道门槛是维度灾难：波段数远多于可用标注样本时，分类器会过拟合，需要降维。第三道门槛是混合像元：高光谱空间分辨率往往较粗，一个像元里可能同时有植被、土壤与水体，直接分类会把混合光谱判成不存在的类别，必须做解混。

本文按「数据到处理到应用」展开。先讲波段特性与数据量，再覆盖噪声与条带校正，然后进入降维与波段选择，接着是光谱解混、分类与亚像元目标检测，最后落到工程实现与质量评估。高光谱与多光谱的对比基础见 遥感影像基础与元数据 ，入射辐亮度到反射率的转换见 辐射定标与大气校正 。

目录

1. 高光谱数据的特性与数据量
2. 波段特性与地物光谱曲线
3. 噪声估计与条带校正
4. 降维：PCA 与 MNF
5. 波段选择与稀疏表示
6. 光谱解混与端元提取
7. 分类：SAM 与光谱库匹配
8. 亚像元目标检测
9. 工程实现与质量评估

1. 高光谱数据的特性与数据量

高光谱传感器分两类。成像光谱仪逐行推扫，空间维与光谱维同时成像，波段数可达几百；快照式传感器一次曝光获得完整数据立方体，帧率高但光谱分辨率较低。星载代表有 EnMAP、PRISMA、GF-5 的 AHSI，机载代表有 AVIRIS-NG。

数据立方体的三个维度是行、列、波段，存储时按波段顺序或按行顺序排列。波段顺序存储（BSQ）便于逐波段处理，行顺序（BIL）与像元顺序（BIP）便于逐像元处理。选哪种取决于处理方式：逐像元做光谱匹配用 BIP 更高效，逐波段做校正用 BSQ。

数据量估算（16 位，单精度处理时乘 2）
传感器规格              单景大小        说明
400x400x224 (AVIRIS)    72 MB           机载小场景
1000x1000x200           400 MB          典型星载
5000x5000x200           10 GB           大区域，需分块

工程上必须分块处理。按行分块（strip）最自然，因为推扫式数据本身就是按行获取的，且光谱维完整保留，逐像元的光谱运算不受影响。分块时要保证块间有重叠，避免空间滤波类操作在边界失真。

内存布局对性能影响很大。Python 里把数据组织成形状为「行、列、波段」的三维数组，逐像元访问用 cube[y, x, :] 得到连续内存，速度快；若按波段优先存储则逐像元访问会跨大步长，缓存命中率低。用内存映射（np.memmap）打开大文件能避免一次性载入。

2. 波段特性与地物光谱曲线

高光谱的价值在于光谱曲线的形状。不同地物的诊断性吸收特征位于特定波长：

诊断性吸收特征（微米）
特征位置    成因              用于识别
0.43~0.45   叶绿素吸收          植被
0.67~0.68   叶绿素红光吸收      植被、胁迫
0.97, 1.20  液态水吸收          植被含水量
1.40, 1.90  大气水汽吸收        需剔除的坏波段
2.20        黏土矿物 Al-OH      高岭石、伊利石
2.33        碳酸盐 CO3         方解石、白云石
1.00, 2.10  铁氧化物            赤铁矿、针铁矿

两个纪律。第一，大气水汽吸收带（1.35 到 1.42 微米、1.80 到 1.95 微米）信噪比极低，几乎不含地表信息，标准流程要把这些波段剔除或做水汽校正，否则它们只贡献噪声。第二，波段要按诊断性特征分组使用，把 2.2 微米附近的窄波段组合起来算矿物指数，比用全波段分类更可解释也更抗噪。

红边区域（0.68 到 0.75 微米）是植被研究的核心。高光谱能分辨红边的位置与斜率，红边位置偏移反映叶绿素含量与胁迫程度，这是多光谱做不到的。

高光谱的波段相关性极高，相邻波段几乎线性相关。这既是冗余（可用降维压缩），也是信息（可用一阶或二阶导数增强吸收特征）。导数光谱能消除基线漂移、放大吸收峰，是矿物识别与植被参数反演的标准预处理。

import numpy as np

def derivative_spectra(cube, order=1):
    # 沿波段维做 Savitzky-Golay 导数，放大吸收特征
    from scipy.signal import savgol_filter
    return savgol_filter(cube, window_length=9, polyorder=2, deriv=order, axis=-1)

3. 噪声估计与条带校正

高光谱的噪声来源有三：光子散粒噪声、探测器读出噪声，以及推扫式传感器各探元响应不一致造成的条带。条带是视觉上最刺眼的伪影，会让同一地物在相邻列呈现明暗交替。

条带校正有两类方法。基于统计的方法（如矩匹配）把每条探元对应的列调整到全局均值与标准差；基于滤波的方法（如低通滤波、小波）分离条带分量并扣除。

import numpy as np

def destripe_moment(image, valid=None):
    # image: (rows, cols)，逐列做矩匹配，抑制探元响应差异
    img = image.astype("float32").copy()
    if valid is None:
        valid = np.ones_like(img, dtype=bool)
    col_mean = np.array([img[valid[:, c], c].mean() for c in range(img.shape[1])])
    col_std = np.array([img[valid[:, c], c].std() + 1e-6 for c in range(img.shape[1])])
    ref_mean, ref_std = col_mean.mean(), col_std.mean()
    for c in range(img.shape[1]):
        img[:, c] = (img[:, c] - col_mean[c]) / col_std[c] * ref_std + ref_mean
    return img

矩匹配假设条带是乘性与加性的线性偏差，对缓变地物效果好，但在地物沿列方向有真实变化时会把真实信号也拉平。稳妥做法是只对均匀区域估计校正系数，再用这些系数应用到全图。

噪声估计用局部方差法：在均匀小窗口内算方差，取所有窗口方差的中位数作为噪声方差估计，或用 MNF 变换后靠后波段的特征值反推。噪声水平决定了后续能做多少降维：信噪比低的波段降维后仍是噪声，需要在降维前先剔除。

信噪比是评估数据质量的核心指标，也是决定能否做定量反演的前提。信噪比低于 100 比 1 时，精细的矿物识别与生化参数反演基本不可靠，只能做粗分类。

4. 降维：PCA 与 MNF

降维是高光谱处理的第一步，目的是压缩冗余、抑制噪声、缓解维度灾难。

PCA 是最基础的方法，对协方差矩阵做特征分解，取特征值最大的若干主成分。它的优点是简单、无参数，缺点是按方差排序，而方差大的方向未必是信息量大的方向，噪声方差大时 PCA 会把噪声排到前面。

MNF（Minimum Noise Fraction）是专为高光谱设计的降维方法，分两步：先用噪声协方差矩阵白化数据，使噪声在各方向方差一致；再对白化后的数据做 PCA。这样排序依据是信噪比而非方差，前面的成分是信噪比高的信息，后面的成分是噪声。MNF 是矿物识别与解混前的标准步骤。

import numpy as np

def mnf(cube, noise=None, n_components=30):
    # cube: (H, W, B)，noise: 各波段噪声标准差估计，缺省用局部方差法
    h, w, b = cube.shape
    X = cube.reshape(-1, b).astype("float64")
    if noise is None:
        noise = _local_noise_std(cube)
    # 噪声协方差并白化
    Cn = np.diag(noise ** 2)
    evals, evecs = np.linalg.eigh(Cn)
    W = evecs @ np.diag(1.0 / np.sqrt(np.maximum(evals, 1e-12)))
    Xw = X @ W
    Cw = np.cov(Xw, rowvar=False)
    ev, evec = np.linalg.eigh(Cw)
    order = np.argsort(ev)[::-1][:n_components]
    proj = W @ evec[:, order]                    # 最终投影矩阵
    return (X @ proj).reshape(h, w, n_components), proj

def _local_noise_std(cube, win=5):
    # 局部方差中位数近似噪声标准差
    from scipy.ndimage import uniform_filter
    mean = uniform_filter(cube, size=(win, win, 1))
    sq = uniform_filter(cube ** 2, size=(win, win, 1))
    var = np.maximum(sq - mean ** 2, 0)
    return np.sqrt(np.median(var, axis=(0, 1)))

降维维度怎么定？看 MNF 特征值曲线，在特征值接近 1（噪声本底）之前截断。经验上 200 波段的高光谱压到 20 到 40 个 MNF 成分能保留绝大部分信息，后续处理量与内存大幅下降。

5. 波段选择与稀疏表示

降维有两类：特征提取（PCA、MNF）把波段线性组合成新特征，波段选择则从原始波段中挑子集。前者压缩率高但新特征失去物理含义，后者保留物理波段，适合需要可解释性或需要保留诊断性吸收特征的任务。

波段选择的准则：

• 信息量：选方差大或熵高的波段。
• 可分性：选类间距离大（如 Jeffries-Matusita 距离）的波段组合。
• 冗余度：用相关系数或互信息去掉高度相关的冗余波段。
• 稀疏性：用稀疏表示或正则化方法选最少而信息最足的波段。

import numpy as np

def jm_distance(mu1, cov1, mu2, cov2):
    # Jeffries-Matusita 距离，衡量两类可分性，值域 0~2
    cov = (cov1 + cov2) / 2.0
    dmu = mu1 - mu2
    bhat = 0.125 * dmu @ np.linalg.pinv(cov) @ dmu + 0.5 * np.log(
        np.linalg.det(cov) / np.sqrt(np.linalg.det(cov1) * np.linalg.det(cov2) + 1e-12) + 1e-12)
    return 2 * (1 - np.exp(-bhat))

工程上常用贪心前向搜索：从空集开始，每次加入使类间可分性提升最大的波段，直到提升饱和。计算量可控，结果可解释。用相关系数矩阵先做粗筛能把候选波段数砍掉一大半，再做精细搜索。

稀疏表示是较新的路线，用少量波段或少量训练样本的线性组合表示待分类像元。它的优势是样本需求少，在标注稀缺的高光谱场景有吸引力，代价是计算量大、参数敏感。

6. 光谱解混与端元提取

混合像元是高光谱的核心问题。粗空间分辨率下，一个像元常包含多种地物，观测光谱是各地物光谱（端元）按面积比例的线性混合。线性混合模型假设：

x = sum_i (a_i * e_i) + n，约束 a_i >= 0 且 sum_i a_i = 1

其中 x 是观测光谱，e_i 是端元光谱，a_i 是丰度。解混分两步：先提取端元，再估计丰度。

端元提取方法：

• N-FINDR：在特征空间中找体积最大的单纯形，顶点即端元。
• VCA（顶点成分分析）：迭代投影找极值点，速度快，是常用方法。
• PPI（像元纯度指数）：投影多次，取极端像元，简单但需要人工定端元数。

import numpy as np

def vca(X, n_endmembers, seed=0):
    # X: (N, B) 像元光谱，返回 (n_endmembers, B) 端元
    rng = np.random.default_rng(seed)
    n, b = X.shape
    Xm = X - X.mean(0)
    U, S, _ = np.linalg.svd(Xm, full_matrices=False)
    Ud = U[:, :n_endmembers]
    idx = []
    for _ in range(n_endmembers):
        w = rng.standard_normal(n_endmembers)
        w /= np.linalg.norm(w)
        v = w @ Ud.T
        proj = Xm @ v
        idx.append(int(np.argmax(np.abs(proj))))
    return X[idx]

丰度估计用全约束最小二乘（FCLS），同时满足非负与和为 1 两个约束。无约束最小二乘会给出负丰度或丰度和不为 1 的解，物理上不成立。

import numpy as np
from scipy.optimize import nnls

def fcls(x, E):
    # 全约束最小二乘：非负且和为 1，用增广矩阵把等式约束并入 nnls
    n_em = E.shape[0]
    delta = 1e-6
    A = np.vstack([E.T, delta * np.ones((1, n_em))])
    b = np.concatenate([x, [delta]])
    a, _ = nnls(A, b)
    return a

端元数怎么定？用 Harsanyi 虚拟维度或特征值拐点估计，通常比真实地物数略多。端元数取多了会出现重复端元，取少了会强迫一个端元表达多种地物。工程上建议先用 VCA 提一批候选，再人工或按光谱库筛掉不合理端元。

解混的验证用重构误差：用丰度与端元重构光谱，与原始光谱的残差反映模型是否合适。残差大的像元往往是端元库没覆盖的地物，或非线性混合（如植被与土壤的紧密混合），此时线性模型失效。

7. 分类：SAM 与光谱库匹配

高光谱分类方法分监督与非监督两类，监督方法又分基于光谱曲线与基于统计特征。

光谱角匹配（SAM）把光谱看作高维空间中的向量，用两向量夹角衡量相似度，对亮度变化不敏感，适合同一地物在不同光照下的匹配。

import numpy as np

def sam(x, refs):
    # x: (B,) 待分类光谱，refs: (C, B) 参考光谱，返回每类夹角（弧度）
    xn = x / (np.linalg.norm(x) + 1e-9)
    rn = refs / (np.linalg.norm(refs, axis=1, keepdims=True) + 1e-9)
    cos = rn @ xn
    return np.arccos(np.clip(cos, -1, 1))

监督分类的经典方法是支持向量机，它在样本少、维度高时表现好，是高光谱分类的长期基准。随机森林与深度学习（一维卷积、光谱注意力网络）近年也常用，后者样本需求大但在样本充足时精度更高。

光谱库匹配是矿物识别的特色路线。USGS、JPL 等机构维护了标准矿物光谱库，把像元光谱与库中光谱做 SAM 或匹配滤波，无需标注样本即可识别矿物种类。前提是数据已做大气校正并转为反射率，否则观测光谱与库光谱不可比。

8. 亚像元目标检测

亚像元检测要在单个像元内判断某目标是否存在，目标占比可能只有百分之几。典型应用是矿物勘探、气体泄漏检测、伪装识别。

两个经典算子：

• ACE（自适应余弦估计）：假设背景服从多元高斯，用背景协方差白化后计算目标方向上的投影，对背景分布做了自适应。
• CEM（约束能量最小化）：设计一个滤波器，使目标响应恒为 1，同时最小化背景能量，输出即为目标丰度。

import numpy as np

def cem(X, d, reg=1e-3):
    # X: (N, B) 背景像元，d: (B,) 目标光谱，返回权重向量
    n, b = X.shape
    R = X.T @ X / n + reg * np.eye(b)
    Rinv = np.linalg.inv(R)
    w = Rinv @ d / (d @ Rinv @ d + 1e-12)
    return w                     # 目标丰度 = X @ w

def ace(X, d):
    n, b = X.shape
    mu = X.mean(0)
    Xc = X - mu
    C = Xc.T @ Xc / n + 1e-6 * np.eye(b)
    Cinv = np.linalg.inv(C)
    dmu = d - mu
    denom = np.sqrt((dmu @ Cinv @ dmu) * ((Xc @ Cinv * Xc).sum(1)))
    return (Xc @ Cinv @ dmu) / (denom + 1e-12)

CEM 对目标光谱准确度很敏感，目标光谱有偏差时性能骤降；ACE 对背景统计做了自适应，鲁棒性更好但对非高斯背景仍会退化。实践中常用两者交叉验证，并结合丰度阈值与形态学后处理抑制虚警。

亚像元检测的最大挑战是虚警率。目标占比低意味着信号弱，背景的微小波动就可能超过阈值。控制虚警的手段有：用恒虚警率检测器按背景统计自适应设阈值、用空间上下文要求目标成片出现、用时序或多角度观测做一致性检验。

9. 工程实现与质量评估

高光谱处理的工程瓶颈是内存与计算量。一景 5000 乘 5000 乘 200 的数据立方体是 10 GB，任何全局操作（协方差、特征分解）都要先把数据读进来。做法是分块统计再合并：先分块算各块的均值与协方差，再按块大小加权合并成全局统计量，最后用全局统计量做变换。

import numpy as np

def streaming_stats(reader, block_rows=256):
    # 分块累加均值与协方差，避免整幅载入内存
    n = 0
    s = None
    ss = None
    for block in reader(block_rows):        # block: (r, W, B)
        X = block.reshape(-1, block.shape[-1]).astype("float64")
        if s is None:
            b = X.shape[1]
            s = np.zeros(b); ss = np.zeros((b, b))
        n += X.shape[0]
        s += X.sum(0)
        ss += X.T @ X
    mean = s / n
    cov = ss / n - np.outer(mean, mean)
    return mean, cov

计算密集的步骤（特征分解、解混、分类）适合并行，可以借助 高性能计算 里的多进程与分块调度思路，把瓦片任务分发到多核。逐像元的解混与检测天然并行，几乎没有通信开销。

质量评估要报三类指标。第一是辐射质量：信噪比、条带残差、坏波段比例。第二是降维质量：MNF 前若干成分解释的信噪比占比、重构误差。第三是应用精度：分类用混淆矩阵与 Kappa，解混用丰度 RMSE 与重构残差，检测用 ROC 曲线与虚警率。

输出产品要附带完整元数据：传感器与波段定义、坏波段列表、降维方法与前若干成分、端元光谱与其来源、解混约束、分类器与训练样本来源。高光谱产品若无这些信息几乎无法复现，跨团队比较也会失去意义。做定量反演时还要记录大气校正方法与气溶胶假设，因为反射率的绝对精度直接决定反演结果的可信度。

权衡取舍

• PCA vs MNF：PCA 简单但按方差排序会把噪声排前，MNF 按信噪比排序更适合高光谱，代价是需要噪声估计。
• 特征提取 vs 波段选择：特征提取压缩率高但失去物理含义，波段选择保留物理波段但压缩率有限，按是否需要可解释性选择。
• 线性解混 vs 非线性：线性模型简单、可解析，但植被与土壤的紧密混合是非线性的，用线性模型会有系统残差。
• 端元数取多 vs 取少：取多会出重复端元、丰度不稳，取少会强迫一个端元表达多种地物，用虚拟维度估计再加人工筛选。
• SAM vs 统计分类：SAM 无需训练样本但依赖参考光谱质量，统计分类精度高但需标注，有光谱库时优先 SAM。
• ACE vs CEM：CEM 对目标光谱准确度敏感，ACE 对背景统计自适应更鲁棒，二者交叉验证。
• 全波段 vs 精选波段：全波段信息最全但噪声与计算量最大，精选诊断性波段更抗噪、更快，按任务选择。

常见坑清单

• 保留水汽吸收波段：现象是分类与解混精度低，原因是 1.4 与 1.9 微米波段信噪比极低，规避方法是按波长剔除坏波段。
• 不做大气校正直接匹配光谱库：现象是矿物识别全错，原因是观测是辐亮度而非反射率，规避方法是先做大气校正转为反射率。
• 无约束最小二乘解混：现象是出现负丰度、丰度和不为 1，原因是未加物理约束，规避方法是改用全约束最小二乘。
• 端元数取太多：现象是丰度图噪声大、端元重复，原因是高估了地物数，规避方法是用虚拟维度估计并筛选端元。
• 降维过度：现象是细分矿物无法区分，原因是把区分性弱的成分也丢了，规避方法是看 MNF 特征值曲线在噪声本底前截断。
• 条带校正拉到真实信号：现象是沿列方向的真实变化被抹平，原因是矩匹配对全图应用，规避方法是只用均匀区域估计校正系数。
• 用 PCA 前几个成分做矿物识别：现象是关键吸收特征被淹没，原因是 PCA 按方差排序忽略低方差诊断特征，规避方法是改用 MNF 或波段选择。
• 训练与测试样本空间重叠：现象是精度虚高，原因是相邻像元泄漏，规避方法是按空间块划分样本。
• CEM 目标光谱不准：现象是虚警率高，原因是滤波器对目标方向极敏感，规避方法是用实测或库光谱并做敏感性分析。
• 忽略混合像元直接硬分类：现象是边界处出现不存在的类别，原因是混合光谱被当作纯光谱，规避方法是先解混再决策。

小结

高光谱的价值来自连续光谱曲线，代价是数据量大、维度高、噪声重。处理链的顺序是有讲究的：先剔除坏波段、做条带与噪声校正，再用 MNF 降维抑制噪声，然后才进入解混、分类或检测。跳过前置校正直接上算法，噪声与条带会被当成信号，结果不可用且难以排查。

落地的关键是分清任务类型。识别矿物或有光谱库时，走光谱库匹配加 SAM，不需要标注样本；做地物分类时，用 MNF 降维加支持向量机或随机森林，样本少也能跑；探测亚像元目标时，用 ACE 与 CEM 交叉验证，重点控制虚警率。每一步都要输出质量标记，让下游知道哪些像元可信。

下一步可以对照 遥感影像分类 理解多光谱分类的经典流程与高光谱分类的差异，也可以结合 辐射定标与大气校正 补足反射率反演的细节，把计算密集的解混与检测任务交给 高性能计算 的分块并行框架。