★ FEATURED ARTICLE 简介面向遥感、精准农业、医学影像与地质勘探等领域的MATLAB高光谱数据读取与处理资源包围绕高动态范围HDR格式高光谱图像为需要快速上手的科研人员和工程技术人员提供完整的示例方案。压缩包内共11个文件合计22.2MB包含可运行的MATLAB脚本hsi_read.m、高光谱数据文件001.dat及其头文件001.hdr、ENVI支持文件001.dat.enp、TIFF格式影像以及若干zbak备份文件和说明txt文档文件类型覆盖从数据读取、解析到格式转换的主要环节便于对照理解HDR高光谱数据的存储结构。内容既讲解将高光谱图像组织为行、列、波段三维数组的基本方法也展示了辐射定标、光谱平滑、主成分分析等预处理技术以及imshow、imagesc、hypercube对象等多种可视化手段进一步结合K均值聚类、支持向量机分类和交叉验证评估形成从数据导入到模型验证的完整分析链路。已有47人学习下载特别适合具备一定MATLAB基础、希望系统掌握高光谱数据处理流程的初学者或相关领域研究人员参考。1. HDR高光谱数据集读取为什么总在MATLAB开头就翻车“MATLAB环境下HDR高光谱图像数据集的读取与处理方法”——刚接触这类数据时我第一反应也是直接load直到连续几个数据集读出来全是斜纹和雪花点才意识到高光谱里的HDR和显示器上的HDR完全是两码事。真实的高光谱HDR数据通常是12位甚至16位原始传感器输出一个立方体包含几十到几百个波段动态范围跨度极大它往往被封装成ENVI二进制、多页TIFF或MAT v7.3格式而不是一张用imread就能打开的普通图片。读取时只要interleave顺序、字节序或位深错一个后面的分析全是白搭。这篇文章按我处理高光谱数据的习惯来写怎么判断封装格式怎么把三维立方体正确读进MATLAB怎么让暗区和亮区同时看得见以及做白板校正前必须先处理哪些坑。适合正在做高光谱课题的从业者也适合为MATLAB图像处理大作业选了这个方向、被读取步骤卡住的学生。2. 按封装格式决定读取路径ENVI、多页TIFF与MAT文件怎么读高光谱数据集在外存上最常见的形态有三种ENVI标准格式一个.dat二进制体配一个.hdr文本头、把波段叠成一页一页的多页TIFF以及MATLAB自己的.mat文件大文件常是v7.3。第一种最麻烦因为hdr里藏着字节序和interleave第二种看着简单但imread可能悄悄把位深压成8位第三种如果想着load一把梭内存和变量维度会让你反复翻车。所以我拿到任何数据集后的第一件事永远是先找元数据文件而不是急着写fread。2.1 ENVI格式hdr元数据里的位深、字节序与interleaveENVI格式在遥感领域出现频率极高。它的典型构成是scene.hdr配scene.dathdr里用文本行记录samples、lines、bands、data type、header offset、interleave、byte order等字段。下面这段代码是我常用的取字段写法hdrLines readlines(scene.hdr); % 需要MATLAB R2020b及以上 getField (key) strtrim(extractAfter( ... char(hdrLines(contains(hdrLines, key))), key)); samples str2double(getField(samples )); lines str2double(getField(lines )); bands str2double(getField(bands )); dataType str2double(getField(data type )); byteOrder str2double(getField(byte order )); interleave char(getField(interleave )); headerOffset str2double(getField(header offset ));readlines把hdr每行读成string元素contains滤出包含关键字的行extractAfter截取等号后面的内容strtrim去掉空格。这样写的好处是字段顺序无关hdr文件里samples不一定在第3行不能用固定行号去读。header offset不带时默认0但我见过不少数据集在二进制前塞了512字节的辅助信息这个字段漏掉会直接让数据整体错位。取值时str2double能自动忽略行尾不可见字符ANSI编码的hdr也能正常工作只有遇到带BOM的UTF-8文件第一行偶尔会多一个隐藏字符对数值字段没有实质影响。hdr里的data type是ENVI自己的编号需要转成MATLAB类型。最常用的映射我直接列成表ENVI data type含义MATLAB读取类型单像素字节数18位无符号整数uint81216位有符号整数int1621216位无符号整数uint1621332位无符号整数uint324432位浮点数single4HDR高光谱文件绝大多数是12号也就是16位无符号整数。接下来按元数据打开二进制文件machineFmt ieee-le; % ENVI里byte order0是小端 if byteOrder 1 machineFmt ieee-be; % 大端常见于旧工作站生成的数据 end fid fopen(scene.dat, rb); if fid -1 error(没有找到scene.dat检查文件路径); end fseek(fid, headerOffset, bof); raw fread(fid, samples*lines*bands, *uint16, 0, machineFmt); fclose(fid);fread前面的星号表示按文件原样读出来不转成double。如果数据量是几千乘几千再乘上百波段转double会让内存直接翻四倍这一步坚持用uint16或single。machineFmt决定文件字节序读之前可以先打印一两个值验证比如读完后的最小值往往在0附近最大值接近满量程若数值乱得像随机噪声先怀疑字节序反了。读取后的一维数组要按interleave reshape成三维立方体if strcmpi(interleave, bsq) % 按波段优先每个波段完整存完再存下一个 cube reshape(raw, [samples, lines, bands]); cube permute(cube, [2, 1, 3]); elseif strcmpi(interleave, bil) % 按行交错一行内所有波段的样本连续存放 cube reshape(raw, [samples, bands, lines]); cube permute(cube, [3, 1, 2]); elseif strcmpi(interleave, bip) % 按像素交错每个像元的所有波段连着存放 cube reshape(raw, [bands, samples, lines]); cube permute(cube, [3, 2, 1]); else error(hdr里的interleave不是bsq/bil/bip需要人工确认); end这里有一个必须记住的顺序reshape先按内存排列把一维数组折成“最左边维度最快”的数组permute再交换维度两个步骤顺序反了图像就会像被斜切过一样。BSQ是最常见的布局reshape后permute一次就行BIL和BIP在遥感里是为了IO效率设计的读取时特别容易错。我的习惯是先取hdr里的samples、lines、bands做乘法再用乘法结果和文件大小对一遍做过校验心里就有底了。2.2 多页TIFF与MAT v7.3两种看似简单但容易走样的读法多页TIFF通常把每个波段存成一页读取可以按页索引来info imfinfo(hdr_cube.tif); nBands numel(info); rows info(1).Height; cols info(1).Width; cube zeros(rows, cols, nBands, uint16); for k 1:nBands band imread(hdr_cube.tif, k); cube(:, :, k) band; end预先用zeros分配好整个三维数组比在循环里不断用cat拼接快很多原因是MATLAB对预分配数组做“写时复制”时不用反复搬移大内存块。这个数据集里如果infos里显示BitDepth8但资料说是16位采集那很可能是保存时已经做过色调映射的显示版本而不是原始HDR数据。遇到位深被压成8位的TIFF直接把它当成物理数据参与反演是不靠谱的只能当作预览参考。对于.mat文件旧版本可以直接load新版v7.3本质是HDF5容器推荐用matfile对象按需读取m matfile(hdr_cube_v73.mat); whos(m) cube m.cube; % 文件太大时只取m.cube(:,:,k)单波段使用matfile时要注意下标顺序和保存时一致。如果保存的是(rows, cols, bands)取第50个波段就是m.cube(:,:,50)别和一维数组索引混在一起。v7.3文件超过4GB很常见whos能先看变量名和尺寸避免load半天才发现变量名不是自己猜的那个。2.3 读取前用文件大小校验元数据三个字段对不上就先别往下走写读取代码之前先用文件大小验证元数据能挡掉一半低级错误dirInfo dir(scene.dat); bytesPerSample 2; % 按uint16计算float32就是4 expected samples * lines * bands * bytesPerSample headerOffset; if dirInfo.bytes expected warning(文件小于预期可能截断、位深或samples/lines对不上); elseif dirInfo.bytes expected fprintf(文件比预期多%d字节\n, dirInfo.bytes - expected); disp(多出来的量通常在几十到几百字节优先怀疑header offset没设); end文件大小不等于samples乘lines乘bands乘单像素字节数时五个元数据字段里至少有一个是错的。多出来的字节数如果只是几十常见原因是header offset少了一大截要先怀疑data type或者interleave判断错了而不是急着改reshape。这一步在压缩包数据集里尤其值得做因为很多人重新打包时会在文件头追加注释。3. 读进来却什么都看不见HDR高光谱的显示归一化与量化读进了三维数组用imagesc显示一个波段屏幕上却是一团漆黑或者一片惨白。这不是数据坏了而是HDR数据的动态范围超出了MATLAB默认的显示能力。高光谱一个完整场景里暗背景区域的DN值可能只有几百而太阳耀斑、镜面反射区域能冲到几万整幅图按全局最小值到最大值线性拉伸高光把标尺拉满后地物细节全被压成黑色。这个矛盾其实是信号本身跨了几个数量级需要选择合理的显示尺度。3.1 先看每波段统计范围搞清楚数据是饱和还是本来暗别急着显示先做一轮统计bandStats zeros(bands, 3); for k 1:bands b cube(:, :, k); bandStats(k, :) [min(b(:)), max(b(:)), mean(b(:))]; end figure; plot(bandStats(:, 3)); xlabel(波段序号); ylabel(平均DN值);mean曲线正常情况下是一条平滑的光谱响应曲线如果某个波段突然冲到65535那一带大概率已经饱和。max等于满量程说明拍摄时增益设置过高对应波段里的反射率信息已经截断后期怎么拉伸都救不回来。相邻波段均值如果出现锯齿状跳变则要检查是不是文件里混入了坏波段或者数据做过有损压缩。3.2 单波段显示用百分位拉伸别让一个亮点毁掉整幅图做显示归一化时我一般用2%到98%百分位代替全局min-maxfunction out bandStretch(band) band single(band); lo prctile(band(:), 2); hi prctile(band(:), 98); out (band - lo) / (hi - lo eps); out min(max(out, 0), 1); endprctile取的2%低分位把传感器暗噪声和零星坏点排除掉98%高分位把极亮的镜面反射去掉。用min-max拉伸时一个太阳耀斑就能让整幅图的对比度全部失效百分位拉伸牺牲了最高的2%灰度级换来了地物纹理的清晰度。注意这个函数输出的0到1浮点数只用于显示和预览不应该直接保存成TIFF参与后续反射率计算。3.3 合成RGB显示波段选取与增强的边界需要做一张接近真彩色的预览图时先选三个波段再分别增强idxR 66; idxG 52; idxB 30; % 根据你自己的波段标号选择 R histeq(mat2gray(squeeze(cube(:,:,idxR)))); G histeq(mat2gray(squeeze(cube(:,:,idxG)))); B histeq(mat2gray(squeeze(cube(:,:,idxB)))); rgb cat(3, R, G, B); imshow(rgb);mat2gray把波段映射到0和1之间histeq再做一次直方图均衡目的是让三个通道的对比度都能被眼睛看到。但直方图均衡会改变灰度的相对关系三个通道分别做histeq后色调会漂移这种增强只适合论文里的“示意”图不能当原始数据展示。要在报告里呈现真实性就只用三次bandStretch再合成效果更保守但颜色保真度高。4. 从能看变成能算暗电流、白板与条带噪声的校正流程HDR高光谱数据从DN值变成可用来反演的相对反射率必须按“先暗电流、后白板”的顺序做。采集时通常每个场景会配一组暗帧和一组白板帧暗帧是盖上镜头盖拍的白板帧是对着标准白板拍的。不把这两步做掉直接用原始DN值做波段比值不同波段间因传感器响应不一致混入的误差会直接污染光谱曲线。4.1 暗电流校正先按波段减暗帧再做白板归一暗电流受温度和曝光时间影响而且不同波段的暗电流偏移不一样我一般按波段做去减darkFrame mean(darkCube, 3); % 若有多帧暗帧取平均降低随机噪声 whiteFrame mean(whiteCube, 3); % 白板帧同样取平均 cal zeros(rows, cols, bands, single); for k 1:bands dark darkFrame(:, :, k); white whiteFrame(:, :, k); denom white - dark; denom(denom 1) 1; % 防止近零值把暗区噪声放大 cal(:, :, k) (single(cube(:, :, k)) - dark) ./ denom; end减暗帧去掉了固定模式噪声白板归一化则消除了波段间的响应差异。denom小于1的像元做钳制是必须的否则白板和暗帧在极端情况下数值几乎相等时除法会造出亮度爆炸的伪像。注意减暗帧要在白板归一之前做顺序反了白板帧里自带的暗电流也会被跟着归一化进去最后结果仍然带偏移。4.2 白板参考帧怎么选同曝光优先于平均值数据集里多组白板帧的曝光参数不一致时直接全家平均会得到一条混合了多个增益的光谱响应校正结果容易出现横向条带。我的原则是先看白板帧与目标图像的拍摄时间或曝光时间字段选择同曝光条件下那几帧求平均如果数据集没有给出曝光信息至少要在hdr或者附属文本里找采集说明。没有白板帧的情况下退一步做暗电流减除加列归一化也能得到视觉上可看的图像但这种方法会抹掉真实的光谱差异计算植被指数或矿物指数时需要谨慎。白板选择这类问题很容易被当成“图像处理细节”跳过实际上它决定的是光谱曲线数值上的绝对刻度错一个拍摄条件反演结果全盘偏移。4.3 坏像素与条带噪声快速检测及边界HDR高光谱的传感器上总有坏点表现为单波段图像中孤立的亮斑或黑点medianBand medfilt2(band, [3 3], symmetric); diffMap abs(band - medianBand); thresh 5 * median(diffMap(:)); % 中位数绝对差乘倍 bad diffMap thresh; if sum(bad(:)) 0.001 * numel(band) warning(坏点比例超过0.1%%可能是条带而不是点状坏点); else band(bad) medianBand(bad); end中值滤波窗口取得小只替换孤立坏点不破坏地物边缘。thresh用中位数的倍数而不是均值倍数原因是HDR数据里个别高亮像元会把均值拉高导致阈值过大漏检。如果检测出的bad像元呈竖线分布这不是点坏而是探测器列响应不一致单个像元替换解决不了需要考虑整列校正或条纹去除算法。面对条带噪声我建议先对照暗帧找成因暗电流造成的条带校正暗帧后通常能缓解响应列偏差造成的条带则需要列方向拟合。5. HDR高光谱读取的五个高频坑现象、原因与解决下面这几条是我在真实高光谱数据集上踩过的记录按“现象、原因、解决”写清楚。读过一遍后你能少花一到两天在同样的地方。5.1 打开文件就Out of Memory现象load一个2GB以上的.mat或直接fread整个ENVI文件MATLAB提示内存不足。原因fread默认把数据转成double读进内存16位文件转double后占用的内存直接翻四倍load大mat文件也存在相同问题。解决读取时显式指定*uint16或*single或分块读chunkRows 512; cube zeros(lines, samples, bands, uint16); for startRow 1:chunkRows:lines nRows min(chunkRows, lines - startRow 1); offset headerOffset (startRow - 1) * samples * bands * 2; fseek(fid, offset, bof); tmp fread(fid, nRows * samples * bands, *uint16, 0, machineFmt); tmp reshape(tmp, [samples, nRows, bands]); cube(startRow:startRownRows-1, :, :) permute(tmp, [2, 1, 3]); end这段按行分块一次只读512行内存占用从“整幅立方体”降成“一截片段”。offset计算中的2是uint16的字节数如果文件是float32要改成4。分块读取时最容易错的是行偏移算错一旦图像出现水平断裂先回来检查这个乘数。5.2 明明16位数据读出来最大却只有255现象文档说数据是16位读取后所有波段的值域都在0到255之间。原因读取时用了imread看图或用了fread但没加星号。imread对许多TIFF和PNG会按自己推断的位深转HDR原始文件被当成8位处理fread不加星号则会把数据转成double但值域本身没变看起来像0到255之间的浮点数。解决先确认元数据里的data type再用fread的*uint16。如果文件已经读进来是uint8数组可以用typecast重新解释整段数据但必须确认原始文件确实是16位否则只是把值域拉宽细节并没有变多。5.3 白板校正后满图噪点加负值现象用白板和暗帧校正后图像不但没有变干净反而出现大量极亮极暗的像素很多像元数值是负的。原因白板帧和暗帧的曝光参数不一致或者白板帧接近或达到饱和也可能是denom在某些波段里接近零而没有钳制。负值则说明减去暗帧后的数值小于零常见于热噪声模型不稳或减暗帧之前没有取平均。解决先对white和dark分别做逐波段均值画出一条光谱曲线对比曝光条件然后给denom加保底钳制。我常用的处理是denom(denom 1) 1并在比例校正后把负值置为零。如果校正后的图像依然有很多条带检查whiteFrame那一波段的均值是否已经饱和饱和的白板帧会直接破坏校正结果。5.4 BIL和BSQ顺序搞反图像像被斜切现象单波段图看起来像被斜着切开边缘错位但统计值min、max都正常。原因interleave判断错了。BSQ是按“波段外循环”存的BIL是按“行内波段交错”存的用错reshape顺序就会把一个波段的错位拼到另一个波段上。解决回到第2.1的reshape逻辑。BSQ用reshape(raw, [samples, lines, bands])再permuteBIL用reshape(raw, [samples, bands, lines])再permute。我判断经验是先取单个像元的几十个波段画一条光谱曲线如果相邻波段均值出现周期性的跳跃大概率是BIL与BSQ搞反了。换过一个数据集这两个姿势不同不要照搬上次的顺序。5.5 保存成TIFF后位深被压成8位现象处理完一个HDR高光谱波段用imwrite保存成tif下次读回来发现值域变成了0到255。原因imwrite对不指定位深的TIFF默认按uint8保存16位高光谱数据被截断到8位动态范围损失大量信息。解决用TIFF接口手动写文件显式设置BitsPerSample为16t Tiff(out.tif, w); tag struct(Photometric, Tiff.Photometric.MinIsBlack, ... BitsPerSample, 16, SamplesPerPixel, 1, ... PlanarConfiguration, Tiff.PlanarConfiguration.Chunky, ... ImageLength, rows, ImageWidth, cols); t.setTag(tag); t.write(uint16(band)); t.close();这里tag里每一个字段都不能缺Photometric决定像素解释方式PlanarConfiguration必须和写入时的二维形状一致。写入之前先检查band的数值范围有没有超过65535否则溢出后会静默回绕成锯齿比压成8位更难发现。6. 先读小块校验再整块读入我最后留下的一个习惯处理数据量大的时候我已经很少直接加载完整数据集去调试了。New宏hdr也好、TIFF索引也好我都先取一个小patch跑通全流程确认形状和数值域再去处理完整文件。这个习惯帮我挡掉过很多次低级错误尤其是那些“看起来一切正常、一画图就方块错乱”的场景。我常把下面这个自检函数放在项目公共目录里每次换新数据集先跑一遍function verifyCube(cube) fprintf(尺寸: %d x %d x %d\n, size(cube, 1), size(cube, 2), size(cube, 3)); for k 1:10:size(cube, 3) b cube(:, :, k); fprintf(波段%3d: min%6.1f max%6.1f mean%6.2f\n, ... k, min(b(:)), max(b(:)), mean(b(:))); end end打印出来的min和max能告诉我位深有没有读错一张图上相邻波段的mean若呈周期波动我会马上检查interleave。另一个实用技巧是利用matfile分时读取校验时只取立方体的前64×64×32对比完整读取的结果数值一致再放开整块。遇到已存的v7.3文件我会先用whos看变量名再用m.cube(1:64,1:64,1:32)这一句取样整个过程不占大内存。我到现在仍然记得第一次被BIL布局坑了一整天所有统计值看起来都正常画出来的图却像被折过的纸片一样错位。从那以后我学会了不迷信单个读取函数而是在处理流程的最前面放一个尺寸、文件大小、单波段统计的校验段。任何声称能自动完成HDR读取的方案我都当成半成品看待真正的可靠性来自你先抽样确认过数据的形状和量纲。希望这些方法和踩坑记录能帮你在MATLAB环境下处理HDR高光谱数据集时少走几步弯路少一点玄学多一点可控。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
