简介高光谱图像以数十至数百个连续波段记录地物光谱信息在农业、地质、环境监测等领域应用广泛。这份面向遥感初学者与MATLAB使用者的压缩包提供了典型的高光谱TIFF样本及配套读取程序可直接上手体验多波段数据的读取与分析流程。包内共6个文件约11.11MB包含tif格式的示例影像、m格式的读取脚本以及hdr、dat、enp等辅助数据文件可帮助理解数据组织与元数据结构。已有3143人学习下载。通过阅读multibandread()等自编函数源码用户不仅能掌握高光谱图像的多波段读取、数据转换与可视化方法还能学习到辐射校正、光谱分析及PCA统计等常见处理思路适合作为MATLAB遥感数据处理入门的实践素材。高光谱图像和MATLAB读取程序每年都有不少刚接触高光谱的师弟师妹来找我开口第一句话几乎都是同一句师兄这个高光谱数据怎么在MATLAB里打开有人双击.raw文件发现是一堆乱码有人用imread读.tif报错还有人加载.mat之后发现变量结构跟想象中完全不一样。高光谱图像的读取这个环节听起来只是预处理的第一步实际上却卡住了大量刚入门的研究者。这个问题的本质在于高光谱图像和普通图像根本不是同一种数据组织方式。RGB图像只有3个波段而高光谱图像动辄几百个波段每个波段都是一张完整的空间灰度图数据被包装成各种各样的格式散落在不同厂商的传感器里。MATLAB本身不是遥感专用软件它不会自动识别每一种高光谱格式所以你需要自己写读取程序或者借助现成工具把数据翻译成MATLAB能理解的多维数组。这篇文章我会从一个实际处理过高光谱数据的研究者角度把手写的读取思路、不同格式的适配方法、常见坑位一次讲清楚。适合那些已经拿到高光谱数据、正准备在MATLAB里做波段分析、分类、异常检测的本科生、研究生和刚转行的工程师无论你的数据来自航空遥感、无人机载高光谱相机还是实验室内的成像光谱仪思路基本通用。1. 先搞清楚高光谱数据到底长什么样1.1 数据维度从一张图变成一个立方体普通照片是二维的宽乘高乘3个颜色通道。高光谱图像则是三维的宽度、高度、波段数。假设传感器有200个波段每个波段的图像分辨率是1000乘1000像素那你手里的数据就是一个1000乘1000乘200的三维数组。在MATLAB里这个数组通常被存储为data\\(h, w, bands\\)的形式h是行数图像高度w是列数图像宽度bands是波段数。这个立方体结构意味着你不仅能看到地物在空间上的样子还能看到它在连续光谱上的反射曲线。每一像素点对应一条光谱曲线而读取程序的任务就是把这个立方体原原本本地从磁盘文件中还原出来。一个很容易让新手困惑的点在于数据在文件里的排列方式和你在MATLAB里建立的数组维度顺序往往不是一一对应的。这就像同样一本书有人按章保存有人按页保存读取时你得知道它的翻页顺序。1.2 为什么不能直接双击打开很多高光谱原始文件.raw、.bip、.bil、.bsq本质上是纯二进制数据流里面只存着数字不包含任何格式说明。哪个字节属于哪个像素、哪个波段范围是多大、每个数值用多少位存储这些信息全都放在另一个文本文件比如.hdr头文件里。没有头文件数据就是一堆无意义的数字有头文件你才能按规则把数字矩阵还原成三维图像。这就好比一个没有目录的档案柜你知道里面全是文件但不知道哪个格子放什么、每份文件有多厚。读取程序做的事情其实是读取头文件里面的目录信息再按照目录去档案柜里取数据。2. 拿到数据之后的第一件事查看元数据2.1 头文件里藏着所有关键参数在写任何读取代码之前先打开头文件看一眼。ENVI格式的.hdr文件是纯文本用记事本就能打开里面会写明samples图像的宽度列数lines图像的高度行数bands波段数量data type数据类型1代表8位无符号整型2代表16位有符号整型4代表32位浮点型12代表16位无符号整型interleave数据排列方式可能是BSQ、BIL或BIPbyte order字节序0表示小端1表示大端其中interleave是最容易忽略也最容易踩坑的参数。它决定了你读数据时输入用哪种维度顺序来reshape。2.2 BIP、BIL、BSQ三种排列方式的区别用生活化的方式来理解这三种排列方式。假设你有一个4乘4像素、3个波段的图像每一个像元在3个波段各有一个值BSQ波段顺序格式先把第一个波段的完整图像全部存完再存第二个波段再存第三个。相当于按波段打包档案每个波段一个完整文件夹。BIL波段按行交叉格式先存第一行3个波段的值再存第二行3个波段的值。好比每一行档案里把不同波段的材料按行叠放。BIP波段按像素交叉格式一个像素3个波段的值连续存放然后下一个像素。相当于每一个像元的完整信息连着放。MATLAB读取时如果interleave没匹配对读出来的图像会呈现条纹状或者完全撕裂。下面这段代码用于读取ENVI格式高光谱图像并处理三种排列方式% 读取ENVI格式高光谱图像的示例程序 filename hyperspectral_image; % 不带扩展名 hdrfile [filename .hdr]; % 解析头文件参数 samples 1000; % 根据实际头文件修改 lines 1000; % 根据实际头文件修改 bands 200; % 根据实际头文件修改 data_type 4; % 4表示32位浮点型 interleave BSQ; % 可选 BSQ BIL BIP % 打开二进制文件 fid fopen([filename .raw], rb); if data_type 4 precision float32; elseif data_type 2 precision int16; elseif data_type 1 precision uint8; elseif data_type 12 precision uint16; else error(未支持的数据类型); end switch upper(interleave) case BSQ data fread(fid, [samples, lines * bands], precision); data reshape(data, [samples, lines, bands]); data permute(data, [2 1 3]); case BIL % BIL: 每个波段逐行交叉存储 data fread(fid, [samples * bands, lines], precision); data reshape(data, [samples, bands, lines]); data permute(data, [3 1 2]); case BIP % BIP: 每个像素的所有波段连续存储 data fread(fid, [bands, samples * lines], precision); data reshape(data, [bands, samples, lines]); data permute(data, [3 2 1]); end fclose(fid); % 去掉单例维度并显示第50个波段的图像 img_band_50 squeeze(data(:, :, 50)); figure; imagesc(img_band_50); colormap(gray); colorbar;提示permute这一步很多人不理解。它的作用是把维度从按文件存储顺序读入的维度调整成MATLAB习惯的(行, 列, 波段)。如果跳过这一步图像会旋转90度或者看起来像被压扁了。2.3 用哪一列作为波段数还有一个常见误区头文件里bands和samples容易反着填。有个简单的验证方法读入数据后随机取一个波段的imagesc看看如果图像显示出来后高宽比例正常、地物轮廓清晰说明行列顺序对了如果图像看起来像被人横着拉长了或者干脆是扭曲的条纹基本就是samples和lines写反了或者permute的维度顺序不对。3. 几种主流高光谱数据格式的读取方案3.1 直接存成 .mat 文件实验室自用场景下最常见的方式是在采集软件里直接把三维数组导出成.mat文件。这种格式对MATLAB最友好直接load就行。但也要注意有些软件的变量名比较随性加载之后建议第一时间用whos查看变量名和维度确认size的维度顺序。load(captured_data.mat); whos;如果变量名不直观可以直接复制给一个标准变量名方便后续操作hsi_cube your_variable_name; clear your_variable_name;3.2 多波段TIF文件很多高光谱相机出厂时提供TIF输出尤其是显微高光谱和部分无人机载设备。TIF文件的一大好处是软件支持好但注意普通imread只能读取前面几个波段或者需要你用imread指定索引逐个读取。推荐做法是用imread循环读取或者用tiff对象直接读取% 循环读取多波段TIF tif_file multiband_image.tif; info imfinfo(tif_file); num_bands numel(info); for k 1:num_bands band imread(tif_file, k); if k 1 [rows, cols] size(band); hsi_cube zeros(rows, cols, num_bands, like, band); end hsi_cube(:, :, k) band; end fprintf(读取完成维度: %d x %d x %d\n, size(hsi_cube));注意imfinfo会返回所有波段的信息对于波段数特别多的TIF文件imfinfo会读取得比较慢但通常可以接受。3.3 MATLAB自带的多光谱/高光谱读取函数MATLAB在较新版本中加强了多光谱图像支持特别是R2021a之后hypercube函数可以直接读取部分格式的高光谱数据并返回一个hypercube对象配合colorize、endmembers等函数使用。% 尝试用MATLAB原生函数读取ENVI格式 hcube hypercube(hyperspectral_image.raw, hyperspectral_image.hdr); disp(hcube); % 查看元数据 info hcube.Metadata; % 提取数据立方体 data_cube hcube.DataCube;如果你的数据是ENVI格式且有较新的MATLAB版本建议优先试一下hypercube。它能自动处理BSQ/BIL/BIP排列省去手动fread的麻烦。不过不是所有高光谱自定义格式都被支持很多实验室自己定义的二进制格式依然只能手写读取代码。3.4 读取HDF5或NetCDF格式部分卫星高光谱数据如EOS系列产品使用HDF格式存储探头数据也有一些用HDF5。MATLAB里读取HDF5用h5read关键是要先查看文件内部数据结构info h5info(satellite_data.h5); disp(info.Datasets); % 假设数据集名称为 /HDFEOS/SWATHS/Radiance/Data Fields/Radiance radiance h5read(satellite_data.h5, ... /HDFEOS/SWATHS/Radiance/Data Fields/Radiance);读取NetCDF类似用ncinfo查看结构再用ncread读取变量。这两类格式的优点在于自描述性强基本不会出现参数不知道的问题缺点是数据结构组织相对复杂需要一层层找进去。4. 写读取程序时的关键细节与坑位4.1 数据类型和字节序决定一切高光谱数据的存储类型多种多样常见的有data type 值含义MATLAB precision18位无符号整型uint8216位有符号整型int16332位有符号整型int32432位浮点型float32564位浮点型float641216位无符号整型uint16这里特别提醒一个隐蔽问题如果写入数据的处理器和你读取数据的电脑字节序不一致会让所有数值变得非常大或非常小。比如本该是反射率0.5的地方读出来是1.68e-38多半就是字节序搞反了。遇到这种情况用memmapfile读取时指定MachineFormat为b或l可以快速切换字节序。4.2 内存不足别急着一次性全部读入高光谱图像的数据量大得惊人。一台1000乘1000像素、200波段的传感器如果用32位浮点存储单文件就是约800MB。有些数据甚至达到几个GB。在MATLAB里一次性fread全部读入非常容易碰到内存不足。工程上常用的处理方案是分块读取。比如你先读取头文件里的行列数然后按波段分块读入% 分块读取BSQ格式数据一次只读一个波段 rows 1000; cols 1000; bands 200; fid fopen(large_image.bsq, rb); % 跳过前n个波段 n_start 1; n_end 50; for k n_start:n_end offset (k-1) * rows * cols * 4; % 4字节/像素 fseek(fid, offset, bof); band fread(fid, [rows, cols], float32); % 处理band... end fclose(fid);这个策略的优势在于你只需要保留当前处理的波段数据大大减少内存压力。处理完成后再决定是逐波段分析还是把压缩后的特征写回磁盘。4.3 暗电流和反射率定标读取之后不要直接用原始高光谱图像通常记录的是DN值数字量化值并不是物理意义上的反射率。使用前一般要做两步预处理暗电流扣除和反射率定标。暗电流是指传感器在没有光照时依然产生的信号需要将原始数据减去暗电流值反射率定标则是利用标准参考板的已知反射率将DN值转换为反射率百分比。MATLAB里做定标比较简单% 假设dark_current是平均暗电流white_ref是参考板数据 calibrated (data_cube - dark_current) ./ (white_ref - dark_current); % 截取到0~1范围 calibrated(calibrated 0) 0; calibrated(calibrated 1) 1;很多高光谱数据读取后直接做分类或匹配效果不佳问题往往出在这一步——数据没有定标不同时间、不同光照条件下采集的数据无法直接对比。4.4 反射率数据可视化怎么用colormap和波段组合读取完数据第一件事通常是看看图像长什么样。除了前面说的用imagesc显示单个灰度波段你还可以挑选三个波段作为R、G、B通道生成伪彩色图像接近真实的自然色或假彩色效果% 假设红波段在60绿波段在35蓝波段在20 red_band squeeze(data_cube(:, :, 60)); green_band squeeze(data_cube(:, :, 35)); blue_band squeeze(data_cube(:, :, 20)); rgb_img zeros(rows, cols, 3); rgb_img(:, :, 1) mat2gray(red_band); rgb_img(:, :, 2) mat2gray(green_band); rgb_img(:, :, 3) mat2gray(blue_band); figure; imshow(rgb_img);这里mat2gray会把每个波段线性拉伸到0到1范围避免因数值范围不一致导致图像过暗或过曝。如果是反射率数据本身已经归一化到0~1则可以直接用。4.5 检查读取是否正确光谱曲线的快速验证读取完成后建议立即检查光谱曲线是否合理。如果数据的空间范围里有植被、水体或裸地你可以用几个典型像素验证。比如植被在红光波段约650nm有强吸收在近红外波段约800nm反射率陡增。曲线如果完全平坦或充满突变尖峰大概率是读取参数有问题需要回头检查头文件参数。% 检查中心像素的光谱曲线 x_pixel round(cols / 2); y_pixel round(rows / 2); spectral_curve squeeze(data_cube(y_pixel, x_pixel, :)); figure; plot(spectral_curve); xlabel(波段数); ylabel(DN值 / 反射率); title(sprintf(像素 (%d, %d) 的光谱曲线, y_pixel, x_pixel)); grid on;5. 读取程序的工程化封装从一次性脚本到可复用函数5.1 把读取逻辑封装成函数我见过很多人的读取脚本是零散的一大段每次换数据都要从头改参数。工程上更推荐把它封装成一个独立的函数输入头文件路径返回三维数据立方体。function [data, meta] read_envi(hdr_path) % 读取ENVI标准格式的高光谱图像 % 输入: hdr_path - .hdr文件的完整路径 % 输出: data - 三维数据 [rows, cols, bands] % meta - 头文件信息的结构体 [path, name, ~] fileparts(hdr_path); raw_file fullfile(path, [name .raw]); if ~exist(raw_file, file) raw_file fullfile(path, name); % 有些文件没有扩展名 end % 解析头文件 meta parse_envi_hdr(hdr_path); fid fopen(raw_file, rb); precision meta.precision; switch upper(meta.interleave) case BSQ data fread(fid, [meta.samples, meta.lines * meta.bands], precision); data reshape(data, [meta.samples, meta.lines, meta.bands]); data permute(data, [2 1 3]); case BIL data fread(fid, [meta.samples * meta.bands, meta.lines], precision); data reshape(data, [meta.samples, meta.bands, meta.lines]); data permute(data, [3 1 2]); case BIP data fread(fid, [meta.bands, meta.samples * meta.lines], precision); data reshape(data, [meta.bands, meta.samples, meta.lines]); data permute(data, [3 2 1]); end fclose(fid); end头文件解析部分可以独立成函数用文本扫描提取关键字段function meta parse_envi_hdr(hdr_path) fid fopen(hdr_path, r); content fread(fid, *char); fclose(fid); % 用正则表达式提取关键参数 meta.samples str2double(regexp(content, samples\s*\s*(\d), tokens, once)); meta.lines str2double(regexp(content, lines\s*\s*(\d), tokens, once)); meta.bands str2double(regexp(content, bands\s*\s*(\d), tokens, once)); data_type_val str2double(regexp(content, data type\s*\s*(\d), tokens, once)); meta.interleave upper(regexp(content, interleave\s*\s*(\w), tokens, once)); type_map containers.Map(... {1, 2, 3, 4, 5, 12}, ... {uint8, int16, int32, float32, float64, uint16}); meta.precision type_map(data_type_val); end封装之后读取任意ENVI数据只需要一行代码[data, meta] read_envi(path/to/your/file.hdr);5.2 读取数据与后续分析的接口设计比较好的实践是读取函数只负责返回原始三维数组和一个元数据结构体不做过多的预处理。预处理放在单独的脚本里方便后续灵活调整。比如你这次想做波段选择下次想做异常检测它们对数据的预处理要求不同把读取和预处理解耦代码复用率会明显更高。5.3 容错处理数据不完整时怎么办采集过程中偶尔会出现文件损坏或数据缺失。建议在读取代码里加一道校验判断读到的大小是否符合预期expected meta.samples * meta.lines * meta.bands; actual numel(data); if actual ~ expected warning(数据元素数量不匹配期望 %d实际 %d, expected, actual); end注意如果文件确实不完整强行reshape会报错所以最好在reshape之前做这个检查。6. 这些坑我都踩过你直接绕开6.1imagesc显示出来为什么是反的MATLAB的imagesc显示图像时第一维是行从上到下第二维是列从左到右。如果你的数据在读取时没有正确做permute图像就会发生旋转或者镜像翻转。我自己的经验是每次读取完数据先看一眼单波段的imagesc确认地面物体的轮廓方向是否正常。如果发现上下颠倒用flipud修正左右颠倒用fliplr修正如果是90度旋转还是回检查permute的顺序更靠谱。6.2 直接imread读取TIF多波段慢得离谱imread在读取多波段TIF时如果你不知道具体波段数用imfinfo查询是合理的但有些TIF的波段间压缩方式不同imread需要逐波段解压速度会非常慢。这时候建议改用tiff对象的readRGBAImage或者直接考虑转成raw格式再读。读取大量波段时用一个进度提示能让你心里有底for k 1:num_bands if mod(k, 50) 0 fprintf(已读取 %d / %d 波段\n, k, num_bands); end band imread(tif_file, k); % ... end6.3 处理整数型数据时别忘了转成double再参与计算很多相机的原始DN值是uint16直接做减法和除法时MATLAB会按整数规则运算结果容易被截断成0。比如(data - dark) ./ (white - dark)如果white - dark比分子大结果会被四舍五入为0整幅定标图像就会变成黑色。解决办法是先转类型% 整数类型先转成double再做定标 z double(data_cube); white double(white_ref); dark double(dark_current); calibrated (z - dark) ./ (white - dark);6.4 文件名不该带扩展名的时候别硬带ENVI格式里这两个文件有严格约定数据文件如.raw、.dat或没有扩展名和头文件.hdr。很多新手写代码时会写成fopen(xxx.hdr, rb)这当然读不出图像数据因为头文件是文本真正的图像数据在另一个文件里。一个稳妥的做法是用fileparts从.hdr路径剥离出主文件名在同一目录下寻找同名但扩展名为.raw/.dat的文件找不到就直接用无扩展名的主文件名。6.5 不要在循环里反复fopen有的人写按波段读取时每读一个波段就fopen一次读完了再fclose这样效率极低。正确做法是在循环外打开文件一次循环内用fseek跳转循环结束后再关闭。类似道理批量读取多文件时也是先按顺序打开所有文件句柄处理完再统一关闭。7. 读取之外给你的数据配一个快速预览流程7.1 建立“秒看”脚本读取高光谱数据之后我强烈建议你固定一个快速预览流程。对我个人而言这个流程包括显示数据维度、绘制中心像素光谱曲线、显示RGB合成图、显示三五个典型波段的灰度图。把这些代码组合成一个预览脚本以后拿到新数据直接跑一遍预览几秒钟内就能判断数据是否正常。function quick_preview(data_cube, wavelengths, bands_rgb) % 快速预览高光谱数据 % bands_rgb: [红波段索引, 绿波段索引, 蓝波段索引] fprintf(数据维度: %d x %d x %d\n, size(data_cube)); [rows, cols, ~] size(data_cube); % RGB合成 rgb_img zeros(rows, cols, 3); for k 1:3 rgb_img(:, :, k) mat2gray(squeeze(data_cube(:, :, bands_rgb(k)))); end figure(Name, RGB合成); imshow(rgb_img); % 中心像素光谱 figure(Name, 中心像素光谱); plot(squeeze(data_cube(round(rows/2), round(cols/2), :))); grid on; end7.2 把读取结果和地理位置关联起来如果你的高光谱数据来自机载或星载平台通常还伴随地理定位信息像元对应的经纬度。在实践项目中很多分析需要把光谱信息与空间位置结合比如异常目标定位、作物长势空间分布反演。读取程序如果能把经纬度数组一并读出来并且保证每个像元的光谱曲线和坐标一一对应后续做地图叠加时会省很多事。7.3 关于MATLAB版本的提醒不同MATLAB版本对高光谱的支持差异很大。R2021a之前的版本没有hypercubeR2023b之后hypercube又增加了部分格式支持。如果你的MATLAB比较老最好不要依赖原生函数手写fread方案最稳。反过来如果你是为了省事升级到新版本也要注意hypercube返回的对象和普通数组的操作方式不太一样部分函数比如size和squeeze用法没有变化但涉及逐波段操作时偶尔需要先取出DataCube字段再继续。8. 写在最后的一点个人体会做高光谱数据处理的头一个礼拜我就因为读取格式的问题耗了整整三天。反复看头文件、改读取顺序、试不同的permute最后发现只是字节序写反了。从那以后我养成了一个习惯拿到任何高光谱数据第一件事不是急着分析而是先写一个最小的读取脚本确认数据形状、数值范围、光谱曲线合理性。这步做扎实了后面所有分析才有基础。希望这篇文章能帮你跳过那些我当年踩过的坑把时间花在真正有价值的光谱分析上。本文还有配套的精品资源点击获取
MATLAB读取高光谱图像：格式解析、读取程序与工程实践
简介高光谱图像以数十至数百个连续波段记录地物光谱信息在农业、地质、环境监测等领域应用广泛。这份面向遥感初学者与MATLAB使用者的压缩包提供了典型的高光谱TIFF样本及配套读取程序可直接上手体验多波段数据的读取与分析流程。包内共6个文件约11.11MB包含tif格式的示例影像、m格式的读取脚本以及hdr、dat、enp等辅助数据文件可帮助理解数据组织与元数据结构。已有3143人学习下载。通过阅读multibandread()等自编函数源码用户不仅能掌握高光谱图像的多波段读取、数据转换与可视化方法还能学习到辐射校正、光谱分析及PCA统计等常见处理思路适合作为MATLAB遥感数据处理入门的实践素材。高光谱图像和MATLAB读取程序每年都有不少刚接触高光谱的师弟师妹来找我开口第一句话几乎都是同一句师兄这个高光谱数据怎么在MATLAB里打开有人双击.raw文件发现是一堆乱码有人用imread读.tif报错还有人加载.mat之后发现变量结构跟想象中完全不一样。高光谱图像的读取这个环节听起来只是预处理的第一步实际上却卡住了大量刚入门的研究者。这个问题的本质在于高光谱图像和普通图像根本不是同一种数据组织方式。RGB图像只有3个波段而高光谱图像动辄几百个波段每个波段都是一张完整的空间灰度图数据被包装成各种各样的格式散落在不同厂商的传感器里。MATLAB本身不是遥感专用软件它不会自动识别每一种高光谱格式所以你需要自己写读取程序或者借助现成工具把数据翻译成MATLAB能理解的多维数组。这篇文章我会从一个实际处理过高光谱数据的研究者角度把手写的读取思路、不同格式的适配方法、常见坑位一次讲清楚。适合那些已经拿到高光谱数据、正准备在MATLAB里做波段分析、分类、异常检测的本科生、研究生和刚转行的工程师无论你的数据来自航空遥感、无人机载高光谱相机还是实验室内的成像光谱仪思路基本通用。1. 先搞清楚高光谱数据到底长什么样1.1 数据维度从一张图变成一个立方体普通照片是二维的宽乘高乘3个颜色通道。高光谱图像则是三维的宽度、高度、波段数。假设传感器有200个波段每个波段的图像分辨率是1000乘1000像素那你手里的数据就是一个1000乘1000乘200的三维数组。在MATLAB里这个数组通常被存储为data\\(h, w, bands\\)的形式h是行数图像高度w是列数图像宽度bands是波段数。这个立方体结构意味着你不仅能看到地物在空间上的样子还能看到它在连续光谱上的反射曲线。每一像素点对应一条光谱曲线而读取程序的任务就是把这个立方体原原本本地从磁盘文件中还原出来。一个很容易让新手困惑的点在于数据在文件里的排列方式和你在MATLAB里建立的数组维度顺序往往不是一一对应的。这就像同样一本书有人按章保存有人按页保存读取时你得知道它的翻页顺序。1.2 为什么不能直接双击打开很多高光谱原始文件.raw、.bip、.bil、.bsq本质上是纯二进制数据流里面只存着数字不包含任何格式说明。哪个字节属于哪个像素、哪个波段范围是多大、每个数值用多少位存储这些信息全都放在另一个文本文件比如.hdr头文件里。没有头文件数据就是一堆无意义的数字有头文件你才能按规则把数字矩阵还原成三维图像。这就好比一个没有目录的档案柜你知道里面全是文件但不知道哪个格子放什么、每份文件有多厚。读取程序做的事情其实是读取头文件里面的目录信息再按照目录去档案柜里取数据。2. 拿到数据之后的第一件事查看元数据2.1 头文件里藏着所有关键参数在写任何读取代码之前先打开头文件看一眼。ENVI格式的.hdr文件是纯文本用记事本就能打开里面会写明samples图像的宽度列数lines图像的高度行数bands波段数量data type数据类型1代表8位无符号整型2代表16位有符号整型4代表32位浮点型12代表16位无符号整型interleave数据排列方式可能是BSQ、BIL或BIPbyte order字节序0表示小端1表示大端其中interleave是最容易忽略也最容易踩坑的参数。它决定了你读数据时输入用哪种维度顺序来reshape。2.2 BIP、BIL、BSQ三种排列方式的区别用生活化的方式来理解这三种排列方式。假设你有一个4乘4像素、3个波段的图像每一个像元在3个波段各有一个值BSQ波段顺序格式先把第一个波段的完整图像全部存完再存第二个波段再存第三个。相当于按波段打包档案每个波段一个完整文件夹。BIL波段按行交叉格式先存第一行3个波段的值再存第二行3个波段的值。好比每一行档案里把不同波段的材料按行叠放。BIP波段按像素交叉格式一个像素3个波段的值连续存放然后下一个像素。相当于每一个像元的完整信息连着放。MATLAB读取时如果interleave没匹配对读出来的图像会呈现条纹状或者完全撕裂。下面这段代码用于读取ENVI格式高光谱图像并处理三种排列方式% 读取ENVI格式高光谱图像的示例程序 filename hyperspectral_image; % 不带扩展名 hdrfile [filename .hdr]; % 解析头文件参数 samples 1000; % 根据实际头文件修改 lines 1000; % 根据实际头文件修改 bands 200; % 根据实际头文件修改 data_type 4; % 4表示32位浮点型 interleave BSQ; % 可选 BSQ BIL BIP % 打开二进制文件 fid fopen([filename .raw], rb); if data_type 4 precision float32; elseif data_type 2 precision int16; elseif data_type 1 precision uint8; elseif data_type 12 precision uint16; else error(未支持的数据类型); end switch upper(interleave) case BSQ data fread(fid, [samples, lines * bands], precision); data reshape(data, [samples, lines, bands]); data permute(data, [2 1 3]); case BIL % BIL: 每个波段逐行交叉存储 data fread(fid, [samples * bands, lines], precision); data reshape(data, [samples, bands, lines]); data permute(data, [3 1 2]); case BIP % BIP: 每个像素的所有波段连续存储 data fread(fid, [bands, samples * lines], precision); data reshape(data, [bands, samples, lines]); data permute(data, [3 2 1]); end fclose(fid); % 去掉单例维度并显示第50个波段的图像 img_band_50 squeeze(data(:, :, 50)); figure; imagesc(img_band_50); colormap(gray); colorbar;提示permute这一步很多人不理解。它的作用是把维度从按文件存储顺序读入的维度调整成MATLAB习惯的(行, 列, 波段)。如果跳过这一步图像会旋转90度或者看起来像被压扁了。2.3 用哪一列作为波段数还有一个常见误区头文件里bands和samples容易反着填。有个简单的验证方法读入数据后随机取一个波段的imagesc看看如果图像显示出来后高宽比例正常、地物轮廓清晰说明行列顺序对了如果图像看起来像被人横着拉长了或者干脆是扭曲的条纹基本就是samples和lines写反了或者permute的维度顺序不对。3. 几种主流高光谱数据格式的读取方案3.1 直接存成 .mat 文件实验室自用场景下最常见的方式是在采集软件里直接把三维数组导出成.mat文件。这种格式对MATLAB最友好直接load就行。但也要注意有些软件的变量名比较随性加载之后建议第一时间用whos查看变量名和维度确认size的维度顺序。load(captured_data.mat); whos;如果变量名不直观可以直接复制给一个标准变量名方便后续操作hsi_cube your_variable_name; clear your_variable_name;3.2 多波段TIF文件很多高光谱相机出厂时提供TIF输出尤其是显微高光谱和部分无人机载设备。TIF文件的一大好处是软件支持好但注意普通imread只能读取前面几个波段或者需要你用imread指定索引逐个读取。推荐做法是用imread循环读取或者用tiff对象直接读取% 循环读取多波段TIF tif_file multiband_image.tif; info imfinfo(tif_file); num_bands numel(info); for k 1:num_bands band imread(tif_file, k); if k 1 [rows, cols] size(band); hsi_cube zeros(rows, cols, num_bands, like, band); end hsi_cube(:, :, k) band; end fprintf(读取完成维度: %d x %d x %d\n, size(hsi_cube));注意imfinfo会返回所有波段的信息对于波段数特别多的TIF文件imfinfo会读取得比较慢但通常可以接受。3.3 MATLAB自带的多光谱/高光谱读取函数MATLAB在较新版本中加强了多光谱图像支持特别是R2021a之后hypercube函数可以直接读取部分格式的高光谱数据并返回一个hypercube对象配合colorize、endmembers等函数使用。% 尝试用MATLAB原生函数读取ENVI格式 hcube hypercube(hyperspectral_image.raw, hyperspectral_image.hdr); disp(hcube); % 查看元数据 info hcube.Metadata; % 提取数据立方体 data_cube hcube.DataCube;如果你的数据是ENVI格式且有较新的MATLAB版本建议优先试一下hypercube。它能自动处理BSQ/BIL/BIP排列省去手动fread的麻烦。不过不是所有高光谱自定义格式都被支持很多实验室自己定义的二进制格式依然只能手写读取代码。3.4 读取HDF5或NetCDF格式部分卫星高光谱数据如EOS系列产品使用HDF格式存储探头数据也有一些用HDF5。MATLAB里读取HDF5用h5read关键是要先查看文件内部数据结构info h5info(satellite_data.h5); disp(info.Datasets); % 假设数据集名称为 /HDFEOS/SWATHS/Radiance/Data Fields/Radiance radiance h5read(satellite_data.h5, ... /HDFEOS/SWATHS/Radiance/Data Fields/Radiance);读取NetCDF类似用ncinfo查看结构再用ncread读取变量。这两类格式的优点在于自描述性强基本不会出现参数不知道的问题缺点是数据结构组织相对复杂需要一层层找进去。4. 写读取程序时的关键细节与坑位4.1 数据类型和字节序决定一切高光谱数据的存储类型多种多样常见的有data type 值含义MATLAB precision18位无符号整型uint8216位有符号整型int16332位有符号整型int32432位浮点型float32564位浮点型float641216位无符号整型uint16这里特别提醒一个隐蔽问题如果写入数据的处理器和你读取数据的电脑字节序不一致会让所有数值变得非常大或非常小。比如本该是反射率0.5的地方读出来是1.68e-38多半就是字节序搞反了。遇到这种情况用memmapfile读取时指定MachineFormat为b或l可以快速切换字节序。4.2 内存不足别急着一次性全部读入高光谱图像的数据量大得惊人。一台1000乘1000像素、200波段的传感器如果用32位浮点存储单文件就是约800MB。有些数据甚至达到几个GB。在MATLAB里一次性fread全部读入非常容易碰到内存不足。工程上常用的处理方案是分块读取。比如你先读取头文件里的行列数然后按波段分块读入% 分块读取BSQ格式数据一次只读一个波段 rows 1000; cols 1000; bands 200; fid fopen(large_image.bsq, rb); % 跳过前n个波段 n_start 1; n_end 50; for k n_start:n_end offset (k-1) * rows * cols * 4; % 4字节/像素 fseek(fid, offset, bof); band fread(fid, [rows, cols], float32); % 处理band... end fclose(fid);这个策略的优势在于你只需要保留当前处理的波段数据大大减少内存压力。处理完成后再决定是逐波段分析还是把压缩后的特征写回磁盘。4.3 暗电流和反射率定标读取之后不要直接用原始高光谱图像通常记录的是DN值数字量化值并不是物理意义上的反射率。使用前一般要做两步预处理暗电流扣除和反射率定标。暗电流是指传感器在没有光照时依然产生的信号需要将原始数据减去暗电流值反射率定标则是利用标准参考板的已知反射率将DN值转换为反射率百分比。MATLAB里做定标比较简单% 假设dark_current是平均暗电流white_ref是参考板数据 calibrated (data_cube - dark_current) ./ (white_ref - dark_current); % 截取到0~1范围 calibrated(calibrated 0) 0; calibrated(calibrated 1) 1;很多高光谱数据读取后直接做分类或匹配效果不佳问题往往出在这一步——数据没有定标不同时间、不同光照条件下采集的数据无法直接对比。4.4 反射率数据可视化怎么用colormap和波段组合读取完数据第一件事通常是看看图像长什么样。除了前面说的用imagesc显示单个灰度波段你还可以挑选三个波段作为R、G、B通道生成伪彩色图像接近真实的自然色或假彩色效果% 假设红波段在60绿波段在35蓝波段在20 red_band squeeze(data_cube(:, :, 60)); green_band squeeze(data_cube(:, :, 35)); blue_band squeeze(data_cube(:, :, 20)); rgb_img zeros(rows, cols, 3); rgb_img(:, :, 1) mat2gray(red_band); rgb_img(:, :, 2) mat2gray(green_band); rgb_img(:, :, 3) mat2gray(blue_band); figure; imshow(rgb_img);这里mat2gray会把每个波段线性拉伸到0到1范围避免因数值范围不一致导致图像过暗或过曝。如果是反射率数据本身已经归一化到0~1则可以直接用。4.5 检查读取是否正确光谱曲线的快速验证读取完成后建议立即检查光谱曲线是否合理。如果数据的空间范围里有植被、水体或裸地你可以用几个典型像素验证。比如植被在红光波段约650nm有强吸收在近红外波段约800nm反射率陡增。曲线如果完全平坦或充满突变尖峰大概率是读取参数有问题需要回头检查头文件参数。% 检查中心像素的光谱曲线 x_pixel round(cols / 2); y_pixel round(rows / 2); spectral_curve squeeze(data_cube(y_pixel, x_pixel, :)); figure; plot(spectral_curve); xlabel(波段数); ylabel(DN值 / 反射率); title(sprintf(像素 (%d, %d) 的光谱曲线, y_pixel, x_pixel)); grid on;5. 读取程序的工程化封装从一次性脚本到可复用函数5.1 把读取逻辑封装成函数我见过很多人的读取脚本是零散的一大段每次换数据都要从头改参数。工程上更推荐把它封装成一个独立的函数输入头文件路径返回三维数据立方体。function [data, meta] read_envi(hdr_path) % 读取ENVI标准格式的高光谱图像 % 输入: hdr_path - .hdr文件的完整路径 % 输出: data - 三维数据 [rows, cols, bands] % meta - 头文件信息的结构体 [path, name, ~] fileparts(hdr_path); raw_file fullfile(path, [name .raw]); if ~exist(raw_file, file) raw_file fullfile(path, name); % 有些文件没有扩展名 end % 解析头文件 meta parse_envi_hdr(hdr_path); fid fopen(raw_file, rb); precision meta.precision; switch upper(meta.interleave) case BSQ data fread(fid, [meta.samples, meta.lines * meta.bands], precision); data reshape(data, [meta.samples, meta.lines, meta.bands]); data permute(data, [2 1 3]); case BIL data fread(fid, [meta.samples * meta.bands, meta.lines], precision); data reshape(data, [meta.samples, meta.bands, meta.lines]); data permute(data, [3 1 2]); case BIP data fread(fid, [meta.bands, meta.samples * meta.lines], precision); data reshape(data, [meta.bands, meta.samples, meta.lines]); data permute(data, [3 2 1]); end fclose(fid); end头文件解析部分可以独立成函数用文本扫描提取关键字段function meta parse_envi_hdr(hdr_path) fid fopen(hdr_path, r); content fread(fid, *char); fclose(fid); % 用正则表达式提取关键参数 meta.samples str2double(regexp(content, samples\s*\s*(\d), tokens, once)); meta.lines str2double(regexp(content, lines\s*\s*(\d), tokens, once)); meta.bands str2double(regexp(content, bands\s*\s*(\d), tokens, once)); data_type_val str2double(regexp(content, data type\s*\s*(\d), tokens, once)); meta.interleave upper(regexp(content, interleave\s*\s*(\w), tokens, once)); type_map containers.Map(... {1, 2, 3, 4, 5, 12}, ... {uint8, int16, int32, float32, float64, uint16}); meta.precision type_map(data_type_val); end封装之后读取任意ENVI数据只需要一行代码[data, meta] read_envi(path/to/your/file.hdr);5.2 读取数据与后续分析的接口设计比较好的实践是读取函数只负责返回原始三维数组和一个元数据结构体不做过多的预处理。预处理放在单独的脚本里方便后续灵活调整。比如你这次想做波段选择下次想做异常检测它们对数据的预处理要求不同把读取和预处理解耦代码复用率会明显更高。5.3 容错处理数据不完整时怎么办采集过程中偶尔会出现文件损坏或数据缺失。建议在读取代码里加一道校验判断读到的大小是否符合预期expected meta.samples * meta.lines * meta.bands; actual numel(data); if actual ~ expected warning(数据元素数量不匹配期望 %d实际 %d, expected, actual); end注意如果文件确实不完整强行reshape会报错所以最好在reshape之前做这个检查。6. 这些坑我都踩过你直接绕开6.1imagesc显示出来为什么是反的MATLAB的imagesc显示图像时第一维是行从上到下第二维是列从左到右。如果你的数据在读取时没有正确做permute图像就会发生旋转或者镜像翻转。我自己的经验是每次读取完数据先看一眼单波段的imagesc确认地面物体的轮廓方向是否正常。如果发现上下颠倒用flipud修正左右颠倒用fliplr修正如果是90度旋转还是回检查permute的顺序更靠谱。6.2 直接imread读取TIF多波段慢得离谱imread在读取多波段TIF时如果你不知道具体波段数用imfinfo查询是合理的但有些TIF的波段间压缩方式不同imread需要逐波段解压速度会非常慢。这时候建议改用tiff对象的readRGBAImage或者直接考虑转成raw格式再读。读取大量波段时用一个进度提示能让你心里有底for k 1:num_bands if mod(k, 50) 0 fprintf(已读取 %d / %d 波段\n, k, num_bands); end band imread(tif_file, k); % ... end6.3 处理整数型数据时别忘了转成double再参与计算很多相机的原始DN值是uint16直接做减法和除法时MATLAB会按整数规则运算结果容易被截断成0。比如(data - dark) ./ (white - dark)如果white - dark比分子大结果会被四舍五入为0整幅定标图像就会变成黑色。解决办法是先转类型% 整数类型先转成double再做定标 z double(data_cube); white double(white_ref); dark double(dark_current); calibrated (z - dark) ./ (white - dark);6.4 文件名不该带扩展名的时候别硬带ENVI格式里这两个文件有严格约定数据文件如.raw、.dat或没有扩展名和头文件.hdr。很多新手写代码时会写成fopen(xxx.hdr, rb)这当然读不出图像数据因为头文件是文本真正的图像数据在另一个文件里。一个稳妥的做法是用fileparts从.hdr路径剥离出主文件名在同一目录下寻找同名但扩展名为.raw/.dat的文件找不到就直接用无扩展名的主文件名。6.5 不要在循环里反复fopen有的人写按波段读取时每读一个波段就fopen一次读完了再fclose这样效率极低。正确做法是在循环外打开文件一次循环内用fseek跳转循环结束后再关闭。类似道理批量读取多文件时也是先按顺序打开所有文件句柄处理完再统一关闭。7. 读取之外给你的数据配一个快速预览流程7.1 建立“秒看”脚本读取高光谱数据之后我强烈建议你固定一个快速预览流程。对我个人而言这个流程包括显示数据维度、绘制中心像素光谱曲线、显示RGB合成图、显示三五个典型波段的灰度图。把这些代码组合成一个预览脚本以后拿到新数据直接跑一遍预览几秒钟内就能判断数据是否正常。function quick_preview(data_cube, wavelengths, bands_rgb) % 快速预览高光谱数据 % bands_rgb: [红波段索引, 绿波段索引, 蓝波段索引] fprintf(数据维度: %d x %d x %d\n, size(data_cube)); [rows, cols, ~] size(data_cube); % RGB合成 rgb_img zeros(rows, cols, 3); for k 1:3 rgb_img(:, :, k) mat2gray(squeeze(data_cube(:, :, bands_rgb(k)))); end figure(Name, RGB合成); imshow(rgb_img); % 中心像素光谱 figure(Name, 中心像素光谱); plot(squeeze(data_cube(round(rows/2), round(cols/2), :))); grid on; end7.2 把读取结果和地理位置关联起来如果你的高光谱数据来自机载或星载平台通常还伴随地理定位信息像元对应的经纬度。在实践项目中很多分析需要把光谱信息与空间位置结合比如异常目标定位、作物长势空间分布反演。读取程序如果能把经纬度数组一并读出来并且保证每个像元的光谱曲线和坐标一一对应后续做地图叠加时会省很多事。7.3 关于MATLAB版本的提醒不同MATLAB版本对高光谱的支持差异很大。R2021a之前的版本没有hypercubeR2023b之后hypercube又增加了部分格式支持。如果你的MATLAB比较老最好不要依赖原生函数手写fread方案最稳。反过来如果你是为了省事升级到新版本也要注意hypercube返回的对象和普通数组的操作方式不太一样部分函数比如size和squeeze用法没有变化但涉及逐波段操作时偶尔需要先取出DataCube字段再继续。8. 写在最后的一点个人体会做高光谱数据处理的头一个礼拜我就因为读取格式的问题耗了整整三天。反复看头文件、改读取顺序、试不同的permute最后发现只是字节序写反了。从那以后我养成了一个习惯拿到任何高光谱数据第一件事不是急着分析而是先写一个最小的读取脚本确认数据形状、数值范围、光谱曲线合理性。这步做扎实了后面所有分析才有基础。希望这篇文章能帮你跳过那些我当年踩过的坑把时间花在真正有价值的光谱分析上。本文还有配套的精品资源点击获取