3、高光谱数据表示与读取：常见数据格式（如ENVI、HDF5、MAT）、使用Python读取和可视化高光谱数据立方体

高光谱图像（Hyperspectral Image, HSI）与普通RGB图像的本质区别在于其“数据立方体”结构：空间维度（行、列）与光谱维度（波段）共同构成一个三维数组。在实际工程中，数据以多种格式存储于磁盘，不同传感器厂商、不同处理软件（ENVI、MATLAB、Python）各有其约定。本章将系统梳理高光谱数据的底层表示逻辑，深入解析ENVI、HDF5、MAT三种主流格式的存储结构，并给出基于Python（numpy、h5py、scipy、rasterio）的读取与可视化完整方案。理解数据格式的差异，是避免“读取出错”“波段顺序颠倒”“反射率定标错误”等高频问题的第一道防线。

上图展示了从数据立方体到三种存储格式，再到Python读取与可视化的完整链路。实际项目中，数据格式的选择往往由传感器或处理软件决定，但读取后的内存布局（波段在前还是在后）需要开发者自行统一。下面逐一展开。

3.1 高光谱数据立方体的底层表示

核心原理：三维数组与波段语义

高光谱数据立方体在数学上是一个三维张量，通常表示为 \( X \in \mathbb{R}^{H \times W \times B} \)，其中 \( H \) 为空间行数（高度），\( W \) 为空间列数（宽度），\( B \) 为光谱波段数。每个空间像元（像素）对应一条长度为 \( B \) 的光谱曲线，记录了该位置在不同波长处的反射率或辐射值。例如，AVIRIS传感器采集的图像通常为 \( 512 \times 512 \times 224 \)，而国产高分五号（GF-5）的可见短波红外数据为 \( 60 \times 330 \times 330 \) 左右（波段数、空间尺寸因产品级别而异）。

在内存中，数据可以按两种主要顺序存储：波段顺序（BSQ, Band Sequential）和波段交叉（BIL, Band Interleaved by Line）。BSQ将每个波段作为一个完整的二维矩阵连续存储，即先存波段1的 \( H \times W \) 矩阵，再存波段2，以此类推。BIL则按行交叉存储，先存第1行的所有波段数据，再存第2行。ENVI格式默认支持BSQ，也支持BIL和BIP（Band Interleaved by Pixel）。理解存储顺序对读取性能至关重要——如果按错误顺序解析，数据将完全错乱。

判断方法：如何识别数据格式与布局

识别数据格式的第一步是查看文件扩展名和伴随文件。ENVI格式通常包含一个ASCII头文件（.hdr）和一个二进制数据文件（.dat、.img或无扩展名）。HDF5文件以.h5或.hdf5结尾，内部采用树状结构存储数据集和属性。MAT文件以.mat结尾，由MATLAB保存，内部包含变量名和数组。此外，还可以通过Python的file命令或读取文件头部字节来初步判断：HDF5文件头包含“HDF5”签名，MAT文件头包含“MATLAB 5.0 MAT-file”文本。

对于波段顺序，可通过头文件中的interleave字段（ENVI）或数据集的shape属性（HDF5）来判断。ENVI头文件中interleave = bsq/bil/bip明确标注；HDF5数据集的shape若为(B, H, W)则波段在前，若为(H, W, B)则波段在后。MAT文件则需查看变量维度，通常MATLAB保存时习惯为(H, W, B)，但需确认。

操作流程：从文件到numpy数组的标准步骤

1. 读取头文件/元数据：提取波段数、行数、列数、数据类型（uint8/uint16/float32）、波长列表、定标系数。
2. 确定存储布局：根据interleave或shape信息决定是否需要转置或重排。
3. 读取二进制数据：使用numpy的fromfile或memmap按指定dtype和偏移量读取。
4. 重塑为三维数组：将一维数据按BSQ/BIL/BIP规则reshape为(H, W, B)或(B, H, W)。
5. 验证数据：检查数值范围、是否存在全零波段、空间维度是否与头文件一致。

常见误区与复盘

• 误区1：忽略字节序。某些传感器数据为big-endian，而x86平台默认little-endian，直接读取会得到乱码。需在fromfile中指定byteorder='>'。
• 误区2：混淆反射率与DN值。DN值（Digital Number）是原始记录，反射率是定标后的物理量。若直接对DN值做分类，结果可能受光照影响。需乘以定标系数或使用反射率产品。
• 误区3：忽略坏波段。水汽吸收波段（如1350-1450nm、1800-2000nm）噪声极大，应在预处理中剔除，否则会干扰分类器。

3.2 ENVI格式详解与读取

核心原理：头文件+二进制数据

ENVI格式是遥感领域最通用的格式之一，由Exelis VIS（现L3Harris）制定。它由两个文件组成：头文件（.hdr）和数据文件。头文件是纯文本，包含关键元数据；数据文件是裸二进制，不含任何自我描述信息。头文件的关键字段包括：samples（列数）、lines（行数）、bands（波段数）、data type（1=uint8, 2=int16, 4=float32, 12=uint16）、interleave（bsq/bil/bip）、byte order（0=little-endian, 1=big-endian）。

操作流程：使用spectral库或纯numpy读取

推荐使用spectral库（需安装：pip install spectral），它内置了ENVI格式读写支持。示例代码如下：

import spectral.io.envi as envi
# 读取头文件和数据
img = envi.open('path/to/file.hdr', 'path/to/file.dat')
# 转换为numpy数组，shape为 (lines, samples, bands)
data = img.load()
print(data.shape, data.dtype)  # 例如 (512, 512, 224) uint16
# 获取波长信息
wavelengths = img.bands.centers  # 单位通常为纳米

若不想依赖spectral库，可手动解析头文件并用numpy读取：

import numpy as np
# 解析头文件（简化版）
with open('file.hdr', 'r') as f:
    header = f.read()
# 提取关键字段（实际需用正则或逐行解析）
samples = 512; lines = 512; bands = 224
dtype = np.uint16  # data type = 12
# 读取二进制数据（BSQ布局）
data = np.fromfile('file.dat', dtype=dtype)
data = data.reshape((bands, lines, samples))
# 转置为 (lines, samples, bands)
data = data.transpose(1, 2, 0)

案例拆解：AVIRIS数据读取与波段筛选

假设我们有一份AVIRIS Classic数据，尺寸为 \( 512 \times 512 \times 224 \)，波长范围400-2500nm。目标是提取特定地物（如植被）的光谱曲线。首先读取数据，然后剔除水汽波段（波段105-115和150-170附近），再计算NDVI（归一化植被指数）并可视化。关键步骤：

1. 读取头文件，确认interleave=bsq，dtype=uint16。
2. 加载数据，检查内存占用：\( 512 \times 512 \times 224 \times 2 \) 字节 ≈ 117 MB，可一次性载入。
3. 剔除坏波段：valid_bands = [i for i in range(224) if not (105 <= i <= 115 or 150 <= i <= 170)]。
4. 计算NDVI：使用波段50（约680nm，红）和波段90（约800nm，近红外）。

3.3 HDF5格式详解与读取

核心原理：分层数据模型与属性

HDF5（Hierarchical Data Format version 5）是一种面向科学计算的自描述格式，支持大规模数据、压缩和并行I/O。高光谱数据常以HDF5存储，尤其是NASA的卫星产品（如HYPXIM、EnMAP模拟数据）。HDF5文件内部结构类似文件系统：组（Group）相当于文件夹，数据集（Dataset）相当于文件，属性（Attribute）存储元数据。例如，一个HDF5文件可能包含组/data，其下有数据集reflectance（shape为\( B \times H \times W \)），以及属性wavelength和scale_factor。

操作流程：使用h5py读取

import h5py
import numpy as np
with h5py.File('hyperspectral.h5', 'r') as f:
    # 查看文件结构
    def print_structure(name, obj):
        print(name, obj.shape if hasattr(obj, 'shape') else '')
    f.visititems(print_structure)
    # 读取数据集（假设路径为 /data/reflectance）
    dataset = f['/data/reflectance']
    data = dataset[:]  # 加载到内存
    # 读取属性
    wavelengths = dataset.attrs['wavelength']
    scale = dataset.attrs.get('scale_factor', 1.0)
    data = data * scale  # 定标

HDF5支持分块读取和内存映射，对于超大立方体（如 \( 10000 \times 10000 \times 200 \)），可使用dataset[1000:2000, :, :]按需加载，避免内存溢出。

案例拆解：EnMAP模拟数据读取与波段子集提取

EnMAP（Environmental Mapping and Analysis Program）模拟数据通常以HDF5格式提供，包含约244个波段，空间分辨率30m。实际处理时，我们往往只需要可见光-近红外（VNIR）部分（波段1-100）。操作流程：

1. 打开文件，定位数据集/EnMAP/reflectance。
2. 读取shape，确认波段顺序为\( B \times H \times W \)。
3. 提取VNIR子集：vnir = data[0:100, :, :]。
4. 转置为\( H \times W \times B \)以便可视化：vnir = vnir.transpose(1, 2, 0)。
5. 保存为npy格式供后续处理。

3.4 MAT格式详解与读取

核心原理：MATLAB变量存储

MAT文件是MATLAB的二进制格式，使用scipy.io.loadmat可读取。高光谱数据在MATLAB中常保存为三维数组，变量名如data或hsi。MAT文件支持v4、v5、v7.3等版本，其中v7.3基于HDF5，读取方式不同。对于v7.3文件，需使用h5py读取。

操作流程：scipy与h5py双方案

import scipy.io as sio
import numpy as np
# 读取v5及以下版本
mat = sio.loadmat('data.mat')
# 查看变量名
print([k for k in mat.keys() if not k.startswith('__')])
data = mat['hsi']  # 假设变量名为hsi
print(data.shape)  # 通常为 (H, W, B)

# 读取v7.3版本（基于HDF5）
import h5py
with h5py.File('data_v73.mat', 'r') as f:
    # 变量名可能带有转义字符
    data = f['hsi'][:]  # shape可能为 (B, H, W)
    data = data.transpose(1, 2, 0)  # 转置为 (H, W, B)

常见误区与复盘

• 误区1：MATLAB保存时自动转置。MATLAB按列优先存储，而Python按行优先，直接读取可能导致维度顺序错误。需检查shape并转置。
• 误区2：v7.3文件无法用scipy读取。scipy不支持v7.3，需改用h5py，且变量名可能包含#refs#等字符，需处理。
• 误区3：忽略数据类型。MATLAB中double类型在Python中为float64，但若原始数据为uint16，MATLAB可能自动转为double，导致内存翻倍。读取后应检查dtype并转换。

3.5 数据可视化：从立方体到图像与光谱曲线

核心原理：波段组合与光谱曲线绘制

可视化高光谱数据主要有两种形式：空间可视化（展示地物分布）和光谱可视化（展示像元光谱特征）。空间可视化通过选择三个波段分别映射到R、G、B通道，生成伪彩色图像。波段选择需根据目标地物特征：例如植被在近红外波段反射率高，可用近红外、红、绿波段组合（如波段90、50、30）突出植被。光谱可视化则提取单个像元或区域平均光谱，绘制反射率随波长变化的曲线。

操作流程：matplotlib与spectral库

import matplotlib.pyplot as plt
import numpy as np

# 假设data为 (H, W, B)，wavelengths为长度B的数组
# 1. 伪彩色合成
rgb_bands = [90, 50, 30]  # 近红外、红、绿
img_rgb = np.stack([data[:, :, i] for i in rgb_bands], axis=-1)
# 归一化到0-1
img_rgb = (img_rgb - img_rgb.min()) / (img_rgb.max() - img_rgb.min())
plt.imshow(img_rgb)
plt.title('False Color Composite (NIR-R-G)')
plt.axis('off')
plt.show()

# 2. 光谱曲线
row, col = 200, 150  # 指定像元
spectrum = data[row, col, :]
plt.figure(figsize=(8, 4))
plt.plot(wavelengths, spectrum, 'b-', linewidth=1.5)
plt.xlabel('Wavelength (nm)')
plt.ylabel('Reflectance')
plt.title(f'Spectral Signature at Pixel ({row}, {col})')
plt.grid(True)
plt.show()

案例拆解：植被与土壤光谱区分

在农业遥感场景中，需要区分植被与裸土。通过绘制两类地物的平均光谱曲线，可以直观看到：植被在680nm附近有吸收谷（叶绿素吸收），在750nm后反射率急剧上升（红边效应）；而土壤光谱则相对平缓。操作步骤：

1. 在图像上手动选取植被区域（如矩形ROI）和土壤区域。
2. 计算区域平均光谱。
3. 绘制两条曲线对比，并标注红边位置（最大斜率点）。

3.6 格式对比与选型建议

选型建议：若数据来自ENVI或需要与ENVI软件交互，优先使用ENVI格式；若数据量大且需要频繁子集读取，HDF5是最佳选择；若团队以MATLAB为主，则保留MAT格式。在Python中，建议统一转换为numpy数组或HDF5作为中间格式，以简化后续处理。

3.7 实战排障：读取失败与数据异常

常见问题与解决方案

• 问题1：读取ENVI数据时shape错误。原因：头文件中的samples和lines顺序混淆。解决：确认samples为列数，lines为行数，reshape时先按(bands, lines, samples)再转置。
• 问题2：HDF5读取时出现KeyError。原因：数据集路径错误或变量名包含特殊字符。解决：使用f.visititems打印所有路径，确认准确名称。
• 问题3：MAT文件读取后数据全为0。原因：MATLAB保存时使用了稀疏矩阵或逻辑类型。解决：检查mat['hsi'].dtype，若为bool或object，需转换。
• 问题4：内存不足。原因：一次性加载整个立方体。解决：使用np.memmap或HDF5分块读取，仅加载需要的波段或空间区域。

案例复盘：一次ENVI头文件解析失败

某项目从机载高光谱传感器获取数据，ENVI头文件中data type = 12（uint16），但读取后数值范围异常（0-65535），且图像出现条纹。排查过程：

1. 检查头文件，发现byte order = 1（big-endian），而Python默认little-endian。
2. 使用np.fromfile(..., dtype='>u2')重新读取，数值恢复正常。
3. 进一步发现部分波段存在坏像元，需用中值滤波修复。

复盘经验：读取任何二进制数据前，务必确认字节序、数据类型和存储布局，三者缺一不可。