MATLAB实现IHS变换图像融合：原理、代码与工程实践

做遥感图像处理的应该都遇到过这种尴尬手里的多光谱影像色彩信息丰富但空间分辨率拉胯地物边缘糊成一片全色影像清晰锐利可惜只有灰度颜色信息为零。两张图各有各的好就是不能完美合体。图像融合要解决的正是这个事而IHS变换融合是这个领域里最经典、也最值得先搞明白的入门算法之一。这篇文章我用MATLAB把IHS变换图像融合从原理到代码完整过一遍包括矩阵变换怎么写、直方图匹配为什么不能省、融合完怎么评估效果以及这个算法有哪些坑、怎么升级。内容对刚接触遥感图像融合的同学友好也适合已经写过一两个融合算法、想系统回看基础原理的从业者。我尽量把每一步的“为什么”说清楚不光是给代码让人照着抄。1. IHS变换到底在融合什么先从色彩空间说起1.1 RGB空间里做不了“干净”的替换很多人第一次接触融合时会有个直觉既然全色影像是灰度图那我直接拿它去替换多光谱影像的R通道让G、B通道不变不就能把细节带进来了吗真要这么干结果一定惨不忍睹整幅影像会严重偏红地物颜色全部乱掉。原因是RGB三个通道描述的是同一场景的红、绿、蓝三波段辐射强度它们彼此高度相关。在植被、农田、水体这类地物结构单一的区域R、G、B三通道的变化几乎是联动的。直接替换其中一个通道相当于人为打破了通道间的平衡关系颜色自然崩坏。RGB本质上是面向显示器设计的色彩空间它追求的是“给人看的时候颜色准”并不适合做图像分析和光谱处理。图像融合真正需要的空间是能把“亮度/结构信息”和“颜色/光谱信息”干净地拆开。亮度管边缘、纹理、几何结构颜色管地物属性两者互不纠缠。IHS变换恰好提供了这样的分离能力。1.2 I、H、S三个分量各管什么IHS色彩空间把一幅RGB影像拆成三个分量IIntensity亮度等于RGB三通道的加权平均描述整体明暗和空间结构视觉效果上基本就是一幅灰度图。HHue色度描述颜色种类对应色调角比如红、绿、蓝、橙、紫各自对应不同的角度值。SSaturation饱和度描述颜色的鲜艳程度灰蒙蒙的颜色饱和度低浓烈鲜艳的颜色饱和度高。把多光谱影像转到IHS空间后I分量几乎承载了原图绝大部分空间结构信息而H和S分量承载的是地物光谱特征。这正好是融合的突破口既然结构信息都浓缩在I里那么把高分辨率的全色影像替换掉I分量空间细节就上去了H和S保持原样颜色又不会崩。这就是IHS融合的命门所在。要提醒的是MATLAB自带的rgb2hsv函数也能提供H、S、V三个分量但它的H归一化到0到1且色度、饱和度定义与遥感文献中常用的IHS矩阵变换并不完全一致。做遥感融合时建议用矩阵变换那一套方便后续逆变换和分量替换操作。1.3 融合的核心思路是“替换”而不是“叠加”很多没接触过融合的人会下意识认为融合就是把两张图按权重叠加在一起比如全色占0.7、多光谱占0.3。实际上IHS融合的思路是“替换”不是“加权混合”。流程浓缩成四步多光谱影像从RGB转到IHS空间得到I、H、S三个分量。对全色影像做直方图匹配让它的灰度分布贴近I分量的灰度分布。用匹配后的全色影像替换原来的I分量。把新的I分量连同原来的H、S分量一起逆变换回RGB得到融合影像。“替换”而不是“叠加”的好处是空间细节完全来自全色影像细节增益更充分色彩信息完全来自多光谱影像只要H、S不变整体色调就能保住。两张图各干各擅长的事互不抢戏。理解了这一点后面所有操作都是顺理成章的。2. 融合流程全拆解从多光谱转IHS到分量替换2.1 数据准备阶段最容易被忽略的三个细节先说一个容易翻车的前提IHS融合对输入数据质量要求很高很多人代码写对了结果仍然惨不忍睹问题基本都出在数据准备阶段。第一分辨率归一化。多光谱影像的分辨率通常低于全色影像比如常见的多光谱是30m全色是15m。融合前必须把多光谱重采样到与全色影像完全一致的网格。MATLAB里用imresize把每个波段放大到目标尺寸即可。这里别偷懒只缩放一个波段然后拼接要保证三个波段的缩放参数完全一致。第二几何配准。两张影像的坐标必须严格对齐。差两三个像素在高分辨率融合结果里就会出现重影道路边缘变成双重轮廓。遇到这种情况先别怀疑算法回去检查配准误差。用地理配准工具或控制点校正都可以但精度必须达到亚像素级。第三位深与数据范围。遥感影像往往是16位整型DN值范围可能到0到65535而MATLAB的很多图像处理函数默认输入输出是0到1的double。直接用im2double处理16位数据没问题但处理8位数据时要注意im2double是按255归一化还是按65535归一化别把两类数据混在一起算。还有一个容易被忽略的点如果影像含有无效值或NoData区域要先做掩膜处理否则直方图统计会被无效值污染后续直方图匹配结果也会有偏差。为了演示方便我用一组模拟数据代表某区域的融合需求ms.tif是8位RGB三波段多光谱影像pan.tif是8位灰度全色影像两者已精确配准且像元尺寸一致。下面代码可以完整跑通。2.2 正变换的矩阵实现与说明我采用的是一组在遥感融合文献里非常常见的正交IHS变换矩阵| I | | 1/3 1/3 1/3 | | R | | v1 | | -1/√6 -1/√6 2/√6 | | G | | v2 | | 1/√2 -1/√2 0 | | B |变换后I (R G B) / 3H atan2(v2, v1)S sqrt(v1² v2²)我在代码里通常直接用v1和v2参与后续计算只要逆变换时v1、v2保持不变色度信息就不会丢失。没有必要非得转换成H角度和S幅值再算回去多一步转换反而引入三角函数误差。为什么选这个矩阵两个原因一是它是正交矩阵逆变换就是转置数值稳定不会放大误差二是I分量的权重是均匀的1/3不偏向某个波段适用面广。很多老文献里用的是带非均匀权重的变换矩阵比如加权平均I 0.5R 0.3G 0.2B那种设计是针对特定传感器数据做过优化的换一个数据源未必合适。MATLAB正变换代码% 读取影像并转为double ms im2double(imread(ms.tif)); % 多光谱影像 pan im2double(imread(pan.tif)); % 全色影像 R ms(:, :, 1); G ms(:, :, 2); B ms(:, :, 3); % IHS正变换 I (R G B) / 3; v1 (-1/sqrt(6)) * R (-1/sqrt(6)) * G (2/sqrt(6)) * B; v2 (1/sqrt(2)) * R (-1/sqrt(2)) * G;2.3 直方图匹配为什么不能直接拿PAN替换I此时I分量就是一幅灰度影像看起来和全色影像很像。但如果你直接把pan替换到I的位置融合结果会出问题首先是亮度水平不一致。比如全色影像整体偏暗原多光谱的I分量整体偏亮替换后融合影像会明显发黑虽然色度还在但明暗关系全错了。其次是对比度不一致。全色影像的动态范围可能比I分量窄替换后融合影像层次感反而下降高光区和暗部细节都会丢失。直方图匹配要解决的就是这个问题让全色影像的灰度分布逼近I分量的灰度分布但保留全色影像自身的空间纹理和边缘信息。打个比方I分量的灰度直方图像一个人的骨架全色影像的纹理像是肌肉直方图匹配相当于把肌肉塑造成骨架需要的形状但肌肉本身的纤维走向不动。MATLAB实现非常简单有一个自带的函数可以直接用% 直方图匹配把pan的灰度分布匹配到I的灰度分布 pan_matched imhistmatch(pan, I);如果你的MATLAB版本比较老没有imhistmatch可以先手动统计两个直方图再做映射。实在嫌麻烦退而求其次用均值方差归一化也能应付大部分场景% 均值方差归一化 pan_norm (pan - mean(pan(:))) / std(pan(:)) * std(I(:)) mean(I(:));直方图匹配比均值方差归一化更强的地方在于它能修正非线性的灰度差异不只是平移和缩放。比如全色影像在暗部非线性压缩、亮部非线性拉伸均值方差归一化处理不了这种差异直方图匹配可以。我自己的习惯是优先用imhistmatch处理完之后可以用imshow对比一下pan_matched和I两者亮度和对比度应该非常接近但pan_matched的边缘明显更锐利。2.4 逆变换与边界处理用匹配后的全色影像替换I分量v1和v2保持不变然后做逆变换。逆矩阵如下| R | | 1 -1/√6 1/√2 | | I | | G | | 1 -1/√6 -1/√2 | | v1 | | B | | 1 2/√6 0 | | v2 |MATLAB代码I_new pan_matched; R_new I_new - (1/sqrt(6)) * v1 (1/sqrt(2)) * v2; G_new I_new - (1/sqrt(6)) * v1 - (1/sqrt(2)) * v2; B_new I_new (2/sqrt(6)) * v1; fused cat(3, R_new, G_new, B_new);逆变换之后最常遇到的问题就是R_new、G_new、B_new的值超出了0到1的范围。原因很简单I分量被替换后原本的RGB数值平衡被打破了某些像素的RGB值在数学上会越界。最粗暴的处理是直接裁剪fused_clip min(max(fused, 0), 1);直接裁剪会带来一个问题高亮区域的细节被“削顶”变成一片纯白丢失纹理。更稳妥的做法是先做一次线性拉伸把三个通道整体的动态范围映射回0到1再做裁剪for k 1:3 ch fused(:, :, k); ch_min min(ch(:)); ch_max max(ch(:)); fused(:, :, k) (ch - ch_min) / (ch_max - ch_min eps); end fused min(max(fused, 0), 1);我个人的经验是如果数据本身没有严重过曝直接clip就够用了画面最自然线性拉伸容易把暗部噪声放大还得配合中值滤波处理。具体用哪种建议先两种都试一下对比局部细节再决定。3. 融合效果怎么看量化指标与实测翻车现场3.1 我常用的几个量化指标及MATLAB实现融合做完了不能只拿眼睛看还得有量化指标。评估融合影像有两条主线一是空间细节增强了多少二是光谱信息保真了多少。指标太多容易眼花我只挑几个最常用的。相关系数CC计算融合影像各波段与原始多光谱对应波段的相关系数。CC越高说明融合影像和原始多光谱在光谱特征上越一致光谱畸变越小。注意这里对比的是融合影像和原始多光谱不是和真值对比所以CC高代表的是“颜色/光谱没跑偏”。光谱角SAM把每个像元的多个波段看作一个光谱向量计算融合影像与原始多光谱对应像元向量的夹角再对全图取平均。SAM越小光谱保真越好。ERGAS综合了波段数、分辨率比例和各波段RMSE的全局误差指标。这个指标在遥感融合论文里出场率很高适合多波段数据。它把空间分辨率比值也考虑进去比较公平。信息熵和平均梯度这两个指标衡量的是空间细节。信息熵越大说明影像的灰度分布越丰富平均梯度越大说明边缘和纹理越锐利。融合后这两个值应该明显高于原始多光谱。MATLAB里自己写CC和SAM都不难function cc calc_cc(a, b) a a(:); b b(:); a a - mean(a); b b - mean(b); cc sum(a .* b) / sqrt(sum(a.^2) * sum(b.^2)); end function sam calc_sam(orig, fused) num sum(orig .* fused, 3); den sqrt(sum(orig.^2, 3)) .* sqrt(sum(fused.^2, 3)); sam acosd(mean(num(:) ./ (den(:) eps))); end3.2 模拟项目上的实测结果与一次翻车现场我拿一组模拟数据跑过完整的IHS融合流程。多光谱影像是从高分辨率彩色图下采样得到的全色影像由彩图转灰度并适度锐化模拟。这种模拟方式的优点是心里有真值能直接比较融合效果和原始高分辨率彩图。跑下来的结果是相关系数CC在0.95以上光谱角SAM保持在6到8度之间平均梯度比原始多光谱提升了将近一倍。这说明空间细节确实注入进去了同时颜色没有大范围失真。视觉上看融合影像的边缘比多光谱锐利得多地物边界清晰植被区域的绿色也基本保留住了。但我在这组数据上还踩过两个坑写出来给大家提个醒。第一个坑一开始我没做直方图匹配直接把pan替换I。结果是融合影像整体发灰对比度特别低SAM值明显偏高。加上直方图匹配之后颜色恢复正常指标也回到正常区间。自此以后我养成了习惯任何时候都要把pan_matched和I放在一起对比看一眼再继续往下走这一步十几秒就能避免一堆后续麻烦。第二个坑换了一批测试数据后融合结果里道路边缘出现双重轮廓。我当时研究了好半天算法以为是IHS变换有问题后来才发现是输入数据没配准全色影像相对多光谱偏移了两个像素左右。重新配准后再跑重影消失。这件事给我的教训是IHS融合对输入数据的几何质量极其敏感数据预处理的工作量可能比写算法本身还大。4. IHS变换的局限与升级思路4.1 光谱失真的根源三波段假设与全局线性变换IHS变换天然是为三波段RGB设计的。你的多光谱如果正好是RGB三波段那没问题但现代遥感数据经常是四波段甚至更多比如蓝、绿、红、近红外。直接用IHS融合只能挑其中三个波段参与近红外信息直接被扔掉了。对植被监测、水体识别这些特别依赖近红外的应用来说这是不可接受的损失。另一个隐患是IHS变换的“全局线性”特性。整个变换是逐像元线性运算替换I分量时默认原来的H、S可以原样保留。但全色影像的灰度变化与多光谱影像的亮度变化并不是处处一致的。在高亮建筑物、浓密阴影、云边界这些特殊区域全色影像带来的空间细节和原I分量之间可能存在明显矛盾强行替换后局部颜色就会偏。这也是为什么很多论文里会说IHS融合光谱失真严重其实严格说不是“严重”而是“特定区域明显”。地物类别单一、结构简单的场景IHS融合完全够用场景一复杂问题就出来了。4.2 与PCA、Brovey等方法的横向对比做融合方案选型时IHS不是唯一选择有几个算法经常拿来一起比较。PCA融合对多光谱各波段做主成分变换第一主成分PC1通常汇聚了绝大部分方差用全色影像匹配后替换PC1再逆变换回去。PCA的优点是数据自适应不依赖固定的三波段假设波段数量也不受限。缺点是主成分变换依赖整幅影像的统计特性如果影像内容差异很大统计结果会被少数地物主导融合效果不稳定。Brovey融合本质是比值变换速度极快适合做快速预览。但光谱失真比较明显颜色经常发飘一般不建议在正式生产流程里用它。小波融合把全色影像做多尺度分解低频分量和多光谱I分量融合高频分量直接注入。保光谱的能力更强空间细节也不错但计算量更大还可能在高频注入时产生振铃伪影。方法波段数限制光谱保真空间细节计算复杂度适用场景IHS三波段中等强低RGB数据、快速验证PCA无中等偏上强中多波段数据、复杂场景Brovey无较弱中等极低实时预览小波无较强中偏上高对光谱保真要求高的任务从我的实测经验看在地物类别简单、边界清晰的场景里IHS和PCA的结果差别不大场景复杂度上来以后PCA和基于多尺度分解的方法优势就开始显现了。IHS最大的价值在于简单、直观、算得快适合作为处理流程的第一版验证方案。4.3 自适应注入细节的改进方向既然IHS的短板是“全局一刀切替换”那改进思路自然就是“局部自适应注入”。我们不把全色影像整个替换I而是先算出全色影像和I分量的高频差异再根据局部区域的特点决定注入多少细节。一个非常朴素的实现是基于局部方差控制注入强度在高纹理区域多注入细节在平坦区域少注入甚至不注入这样既保空间细节又降低光谱失真。% 计算全色影像与I分量的细节差异 detail pan_matched - I; % 局部方差 localVar stdfilt(I, ones(9)); % 归一化权重 w localVar ./ max(localVar(:)); % 自适应融合亮度分量 I_fused I w .* detail;用I_fused替代I_new再走逆变换效果通常比直接替换更温柔颜色过渡更自然。当然这个方法只是抛砖引玉真正工程化还要考虑噪声放大、云影覆盖、山地阴影等复杂因素。比如局部方差对噪声很敏感平坦区域噪声稍微大一点权重就会被顶高反而把噪声细节注入了进去。一个折中的办法是先对detail做高斯模糊或边缘保持滤波降低噪声的影响。这块内容值得单独写一篇去展开。我目前这套自适应注入的思路在公司内部处理某区域地表覆盖数据时已经替代了最原始的IHS流程SAM指标平均能压低2到3度实用性还是看得见的。这几年我处理遥感影像遇到融合需求时第一反应往往是先跑一遍IHS不是因为它最先进而是因为它的结果稳定、代码容易排查还能快速暴露数据配准和预处理的问题。如果你也想上手建议先用下采样模拟数据把整个流程理顺再换真实数据调试的时候多打印中间结果尤其是直方图匹配前后的影像很多问题一眼就能看出来。后面如果大家感兴趣我可以再聊聊如何把这段流程封装成函数批量处理多景影像或者怎么和深度学习超分的结果做交叉验证。有具体问题评论区见。