简介面向遥感图像处理初学者、课程设计及科研人员这份MATLAB资源围绕HIS/IHS色彩空间变换演示如何把多光谱与全色影像进行融合以提升空间分辨率同时保留光谱信息。压缩包共7个文件包括4幅TIFF格式测试影像含多光谱与全色数据、1个可直接运行的.m主程序、1篇IHS变换与自适应区域特征的算法论文及LICENSE说明整体仅9.38MB结构紧凑。已有812人浏览学习。内容覆盖从算法原理到代码实现的闭环借助主程序和配套影像可完整体验从图像预处理、RGB到HIS转换、特征选择、融合算法实施到结果评估的完整链路论文提供原理推导与参数设定依据尤其适合课程设计、毕业设计或遥感融合算法复现。对照代码与文献能快速理解HIS融合在提升影像细节、改善目视效果方面的实际应用价值。1. HIS不是医院信息系统遥感图像融合里的IHS变换到底干什么搜索“HIS遥感图像处理”大概率先撞上医院信息系统Hospital Information System这在遥感领域里纯属同名误会。遥感图像处理中更准确的拼写是 IHSIntensity-Hue-Saturation即亮度-色度-饱和度变换它是遥感图像融合里最经典、最容易被低估的入口。核心操作一句话用高分辨率全色波段PAN替换多光谱影像MS中代表亮度信息的 I 分量再做一次逆变换从而在保留多光谱色彩的同时获得全色的空间分辨率。这套流程没有神经网络、不需要训练样本纯代数变换就能把 2 米的多光谱“提”到 0.5 米。适合刚接触遥感图像融合的人做第一个可复现基线也适合老手在评估新融合算法时拿来当标尺。2. IHS融合的数学原理与波段对应为什么换的是I不是H也不是S2.1 RGB到IHS的正反变换从三通道到亮度-色度-饱和度IHS色彩空间将RGB三波段重新分解为三个互不相关的分量I亮度反映光谱辐射总强度承载空间纹理信息H色度描述主波长决定色调偏红、偏蓝还是偏绿S饱和度描述颜色与灰度的偏离程度从RGB到IHS有球体变换、圆柱体变换和三角形变换等多种数学模型。遥感图像处理中用得最多的是通用矩阵形式的线性近似变换因为计算效率高且逆变换矩阵稳定。正向变换定义如下[ I ] [ 1/3 1/3 1/3 ] [ R ] [v1] [ 1/√2 -1/√2 0 ] [ G ] [v2] [ 1/√6 1/√6 -2/√6] [ B ]对应的逆变换为R I - (1/√2) * v1 - (1/√6) * v2 G I (1/√2) * v1 - (1/√6) * v2 B I (2/√6) * v2注意这里采用的线性IHS是严格正交线性变换可以无损重建回RGB。遥感图像融合中许多开源库和商业软件如ENVI的标准IHS工具使用的原理与此一致只是矩阵系数的取值略有差异。使用该变换的关键前提是输入波段必须经过辐射定标和大气校正否则I分量会携带云影和气溶胶的噪声污染后续替换结果。import numpy as np def rgb_to_ihs(r, g, b): i (r g b) / 3.0 v1 (g - r) / np.sqrt(2.0) v2 (2.0 * b - r - g) / np.sqrt(6.0) return i, v1, v2 def ihs_to_rgb(i, v1, v2): r i - v1 / np.sqrt(2.0) - v2 / np.sqrt(6.0) g i v1 / np.sqrt(2.0) - v2 / np.sqrt(6.0) b i 2.0 * v2 / np.sqrt(6.0) return r, g, b这段代码是实现IHS遥感图像融合的骨架。符号含义r、g、b为同一像素在三个波段上的反射率值v1、v2是变换后的中间色度分量本身不参与替换只在逆变换时用。正变换输出中只有I分量与全色波段在物理含义上对应这也是后续替换操作的理论依据。2.2 三步融合流程正变换、替换、逆变换IHS融合的流程在遥感图像处理教科书中被总结为标准三步几乎所有后续变体都是在这三步上做文章将多光谱影像的三个波段执行RGB到IHS正变换得到I、H、S此处用v1、v2代替H、S对全色波段做直方图匹配或均值方差拉伸使其与I分量具有相同的一阶矩和二阶矩然后用处理后的全色波段替换I分量以替换后的I、原始v1、v2执行逆变换生成融合后的RGB影像其中第二步是关键设计决策。如果直接用原始PAN替换I而不做任何拉伸融合结果会出现整体偏亮或偏暗的问题因为PAN的均值与I分量的均值往往相差很大。常见做法是使用如下线性拉伸pan_stretched (pan - pan.mean()) * (i.std() / pan.std()) i.mean()这个公式将PAN的均值对齐到I分量的均值把PAN的标准差对齐到I分量的标准差。参数含义pan.mean()和pan.std()是全色波段的全局统计量i.mean()和i.std()是I分量的全局统计量。效果是保留PAN的高频空间细节同时让它的亮度水平与MS影像保持一致减少融合后的色偏。对遥感图像融合而言这一步直接决定最终结果的质量。2.3 为什么替换的是I而不是H或S光谱失真的根源从信息论角度看RGB到IHS的变换本质上是坐标旋转。I分量几乎承载了三个波段的公共辐射能量也就是空间结构信息而v1和v2承载的是波段间的差异即色彩信息。因此I分量是高频纹理和亮度信息的理想载体这和PAN影像的物理特性天然匹配。如果试图替换H或S分量结果会非常糟糕H代表色调替换后植被、水体、裸土的颜色会被彻底打乱S代表饱和度替换后整个影像要么灰蒙蒙要么颜色浓到溢出。只有I分量和PAN之间相关性最强替换带来的光谱扭曲最小。但所谓“最小”只是相对而言线性IHS变换本身无法彻底消除色彩畸变这就是后续章节GS变换、小波变换出现的原因——它们本质上都在寻找一个比I分量更优的替换空间。3. 用PythonGDAL在本地跑通IHS融合的最小命令3.1 环境准备与数据读取遥感图像融合的最小环境是Python加上GDAL、NumPy和Rasterio三个库。GDAL负责读取常用遥感影像格式GeoTIFF、ENVI、HDF4/5Rasterio提供更友好的数组接口NumPy承担矩阵运算。pip install gdal rasterio numpy这里不指定版本号因为GDAL的版本与Python版本严格绑定建议在conda环境中用conda install gdal安装预编译版本。数据准备方面需要一对经过配准的影像PAN影像单波段空间分辨率高如0.5mMS影像至少三个波段空间分辨率低如2m两幅影像必须位于同一投影坐标系范围可以略有出入但重叠区要达到95%以上。以下代码读取影像并将MS重采样到PAN的分辨率import rasterio from rasterio.enums import Resampling import numpy as np with rasterio.open(pan.tif) as src: pan src.read(1) pan_profile src.profile with rasterio.open(ms.tif) as src: ms src.read() # 形状: (band_count, height, width) ms_profile src.profile # 将MS重采样到PAN影像的尺寸 with rasterio.open(ms.tif) as src: ms_resampled src.read( out_shape(src.count, pan.shape[0], pan.shape[1]), resamplingResampling.cubic ) print(PAN shape:, pan.shape) print(MS shape after resample:, ms_resampled.shape)参数说明Resampling.cubic表示三次卷积插值比最近邻法nearest更适合光学影像的多光谱重采样因为保留了波段间的光谱关系。out_shape是重采样后的目标尺寸这里直接对齐到PAN的宽高。读取的ms_resampled是一个三维数组维度顺序是波段优先后面所有处理都沿用这个布局。3.2 正变换、全色替换与逆变换的实现核心融合代码可以直接复用第2章的rgb_to_ihs和ihs_to_rgb函数但需要处理一个实际问题遥感影像的DN值范围不是0-255一般8位数据是0-25516位数据是0-65535归一化到浮点型后变换更稳定。# 取前三个波段做RGB融合默认按R-G-B顺序 r ms_resampled[0].astype(np.float64) g ms_resampled[1].astype(np.float64) b ms_resampled[2].astype(np.float64) # 正变换 i, v1, v2 rgb_to_ihs(r, g, b) # 全色波段均值/方差对齐到I分量 pan_f pan.astype(np.float64) pan_stretched (pan_f - pan_f.mean()) * (i.std() / pan_f.std()) i.mean() # 替换I分量 i_new pan_stretched # 逆变换得到融合后的R、G、B r_fused, g_fused, b_fused ihs_to_rgb(i_new, v1, v2) # 裁剪到输入数据的有效值范围 r_fused np.clip(r_fused, 0, np.iinfo(ms.dtype).max) g_fused np.clip(g_fused, 0, np.iinfo(ms.dtype).max) b_fused np.clip(b_fused, 0, np.iinfo(ms.dtype).max)参数说明i.std() / pan_f.std()是标准差比值控制替换后影像的对比度与I分量一致加上i.mean()是亮度对齐。这里的np.clip用于把超出有效值范围的预测值拉回原影像的数据类型范围否则输出16位影像时会出现溢出噪声。将替代后的三个波段拼回数组并写盘# 合并融合波段并写为GeoTIFF fused np.stack([r_fused, g_fused, b_fused], axis0) out_profile ms_profile.copy() out_profile.update({ height: fused.shape[1], width: fused.shape[2], transform: pan_profile[transform] }) with rasterio.open(fused_ihs.tif, w, **out_profile) as dst: dst.write(fused.astype(ms.dtype))这里把输出文件的投影信息与地理变换参数替换为PAN的元数据因为融合后影像的几何分辨率继承自PAN。ms_profile.copy()保留了波段数、压缩格式、数据类型的原始设置transform字段必须更新不然输出的地理坐标会错位。3.3 输出检查先看波谱曲线再看空间纹理融合完成后不能只靠肉眼观察。遥感图像处理的标准质检流程是同时看两类指标空间纹理是否达到全色分辨率放大到像素级别检查道路边缘、屋顶轮廓是否锐利光谱色彩是否接近原MS影像在未变化区域对比融合前后典型地物的光谱曲线from osgeo import gdal ds gdal.Open(fused_ihs.tif) bands ds.RasterCount print(融合影像波段数:, bands) # 输出单波段统计量作为质量记录 for b in range(1, bands 1): band ds.GetRasterBand(b) stats band.ComputeStatistics(False) print(fBand {b}: min{stats[0]:.2f}, max{stats[1]:.2f}, mean{stats[2]:.2f}, std{stats[3]:.2f})波段数通常为3如果你希望保留更多光谱波段如8波段多光谱需要多次执行IHS融合或者改用面向任意波段数的GS变换这一点在第5章展开。检查统计量的目的在于如果某波段的标准差明显高于原始MS影像说明融合过程中发生了严重的光谱拉伸可能需要回头调整直方图匹配的参数。4. IHS融合参数设置与常见坑波段、重采样与直方图匹配4.1 三个波段怎么选真彩色还是假彩色IHS变换一次只能处理三个波段多光谱只要超过3个波段就面临选波段的问题。不同卫星的波段组合习惯不同必须按应用场景决定应用场景推荐波段组合效果地表生态、植被分析R近红外G红B绿植被呈红色适合做植被指数目视解译城市地物分类R红G绿B蓝接近人眼视觉道路与建筑边界清晰水体/湿地监测R近红外G短波红外B红水体呈深黑色湿地边界锐利选择波段时优先选相关性低的三个波段这样v1和v2的方差贡献更大色彩信息保留更充分。如果三个波段高度相关例如红、绿、蓝的原始反射率在植被区接近饱和融合后色彩区分度会显著下降。4.2 重采样方法与分辨率匹配多光谱影像重采样到全色分辨率这一步最容易踩坑。不同插值方法对光谱保持能力的影响如下表nearest速度快但不产生新值容易造成地物边缘锯齿光谱曲线保留最好bilinear线性加权平滑效果好但会轻微模糊纹理cubic三次卷积纹理锐化明显但容易产生超出原值范围的震荡lanczos窗口更大精度更高但计算量成倍增加对于IHS融合推荐使用cubic或lanczos。理由是PAN波段本身分辨率远高于MSMS重采样后要与PAN逐像素对齐插值误差会被融合过程放大。nearest虽然在光谱上最安全但几何对齐误差会直接导致融合结果出现“色彩重影”也就是地物边缘出现一道亮边。另外要注意重采样的边界处理默认情况下GDAL会在影像边缘外推值这会让边缘几列像素出现异常值建议在重采样前先对影像做适当的对称填充。4.3 直方图匹配与光谱失真的三个典型场景直方图匹配是IHS融合中唯一能主动调节光谱保真度的步骤。直接替换I分量在以下三种场景中会引发严重光谱失真场景一PAN与MS成像时间不同步。例如PAN摄取于6月MS摄取于9月两者均值相差悬殊不做对齐会出现整体色偏。场景二PAN包含云影区域。云的亮度极高替换后I分量被抬升IHS逆变换会生成偏白的地物。场景三水体或阴影区域。PAN在这些区域的信噪比极低替换后I分量出现块状噪声逆变换后的色彩完全偏离。针对以上问题常见做法是引入局部直方图匹配。全局匹配不够用就按分块统计直方图。以256×256像元为窗口计算滑动窗口内的均值与方差对PAN做局部对齐后再进行替换可以显著减少局部辐射差异的影响。def local_hist_match(pan, i, window256): from scipy.ndimage import uniform_filter pan_mean uniform_filter(pan, sizewindow, modereflect) i_mean uniform_filter(i, sizewindow, modereflect) pan_std np.sqrt(uniform_filter((pan - pan_mean)**2, sizewindow, modereflect)) i_std np.sqrt(uniform_filter((i - i_mean)**2, sizewindow, modereflect)) # 避免除零 i_std[i_std 1e-6] 1e-6 return (pan - pan_mean) * (i_std / pan_std) i_mean参数说明window是滑动窗口大小值越小局部辐射校正越强但过小会把空间纹理中的真实反差吞掉一般取PAN影像宽度的1/20到1/50之间。modereflect是边界填充方式防止窗口滑到边缘时数据不足。这个函数替换掉第3章中的全局拉伸公式后融合色彩还原度会有肉眼可见的提升。5. 选型对比IHS与Brovey、GS、小波融合怎么选以及一个验证技巧用IHS完成第一版融合后建议顺手跟其他三种方法做一次对比因为它们分别代表了互补的思路。方法替换空间波段数光谱保真计算成本典型场景IHSRGB色彩空间3中等色彩易偏极低快速预览、真彩色合成BroveyRGB比例3较差高亮区易过曝极低单张快出图Gram-SchmidtGS任意统计空间任意波段高低多光谱全波段融合小波变换频率子带任意波段高高高分影像精细融合GS变换是IHS在数学上的推广它用统计方法生成第一分量与I分量类似再做正交化替换因此能处理4波段以上的多光谱数据光谱保真度也更稳。小波变换则把空间信息分解为高频和低频只替换低频部分避免全色波段全盘替换导致的光谱突变。如果你的数据是WorldView或高分一号这种多波段传感器GS是比IHS更合理的基线。一个实用的验证技巧是计算融合前后影像的相关系数CC和相对全局综合误差ERGAS。相关系数衡量光谱保真度ERGAS衡量整体质量公式如下CC 计算原始MS与融合影像对应波段的相关系数逐波段平均 ERGAS 100 * (h/l) * sqrt(mean( (RMSE_b / mean_b)^2 ))其中h/l是PAN与MS分辨率之比。如果ERGAS值小于3说明融合质量很好大于5则需要检查参数设置。计算代码很轻量def ergas(original, fused, ratio4): bands original.shape[0] mse_list [] for b in range(bands): mse np.mean((original[b] - fused[b]) ** 2) mse_list.append(mse / np.mean(original[b]) ** 2) return 100 * ratio * np.sqrt(np.mean(mse_list))把第3章的fused_ihs.tif和重采样后的原始MS影像喂进这个函数快速就能判断该不该继续调参。这个技巧的价值在于遥感图像融合的方法选型不该靠肉眼也不该靠“新方法一定更好”的直觉用定量的ERGAS数据说话比换十个算法都管用。本文还有配套的精品资源点击获取
简介面向遥感图像处理初学者、课程设计及科研人员这份MATLAB资源围绕HIS/IHS色彩空间变换演示如何把多光谱与全色影像进行融合以提升空间分辨率同时保留光谱信息。压缩包共7个文件包括4幅TIFF格式测试影像含多光谱与全色数据、1个可直接运行的.m主程序、1篇IHS变换与自适应区域特征的算法论文及LICENSE说明整体仅9.38MB结构紧凑。已有812人浏览学习。内容覆盖从算法原理到代码实现的闭环借助主程序和配套影像可完整体验从图像预处理、RGB到HIS转换、特征选择、融合算法实施到结果评估的完整链路论文提供原理推导与参数设定依据尤其适合课程设计、毕业设计或遥感融合算法复现。对照代码与文献能快速理解HIS融合在提升影像细节、改善目视效果方面的实际应用价值。1. HIS不是医院信息系统遥感图像融合里的IHS变换到底干什么搜索“HIS遥感图像处理”大概率先撞上医院信息系统Hospital Information System这在遥感领域里纯属同名误会。遥感图像处理中更准确的拼写是 IHSIntensity-Hue-Saturation即亮度-色度-饱和度变换它是遥感图像融合里最经典、最容易被低估的入口。核心操作一句话用高分辨率全色波段PAN替换多光谱影像MS中代表亮度信息的 I 分量再做一次逆变换从而在保留多光谱色彩的同时获得全色的空间分辨率。这套流程没有神经网络、不需要训练样本纯代数变换就能把 2 米的多光谱“提”到 0.5 米。适合刚接触遥感图像融合的人做第一个可复现基线也适合老手在评估新融合算法时拿来当标尺。2. IHS融合的数学原理与波段对应为什么换的是I不是H也不是S2.1 RGB到IHS的正反变换从三通道到亮度-色度-饱和度IHS色彩空间将RGB三波段重新分解为三个互不相关的分量I亮度反映光谱辐射总强度承载空间纹理信息H色度描述主波长决定色调偏红、偏蓝还是偏绿S饱和度描述颜色与灰度的偏离程度从RGB到IHS有球体变换、圆柱体变换和三角形变换等多种数学模型。遥感图像处理中用得最多的是通用矩阵形式的线性近似变换因为计算效率高且逆变换矩阵稳定。正向变换定义如下[ I ] [ 1/3 1/3 1/3 ] [ R ] [v1] [ 1/√2 -1/√2 0 ] [ G ] [v2] [ 1/√6 1/√6 -2/√6] [ B ]对应的逆变换为R I - (1/√2) * v1 - (1/√6) * v2 G I (1/√2) * v1 - (1/√6) * v2 B I (2/√6) * v2注意这里采用的线性IHS是严格正交线性变换可以无损重建回RGB。遥感图像融合中许多开源库和商业软件如ENVI的标准IHS工具使用的原理与此一致只是矩阵系数的取值略有差异。使用该变换的关键前提是输入波段必须经过辐射定标和大气校正否则I分量会携带云影和气溶胶的噪声污染后续替换结果。import numpy as np def rgb_to_ihs(r, g, b): i (r g b) / 3.0 v1 (g - r) / np.sqrt(2.0) v2 (2.0 * b - r - g) / np.sqrt(6.0) return i, v1, v2 def ihs_to_rgb(i, v1, v2): r i - v1 / np.sqrt(2.0) - v2 / np.sqrt(6.0) g i v1 / np.sqrt(2.0) - v2 / np.sqrt(6.0) b i 2.0 * v2 / np.sqrt(6.0) return r, g, b这段代码是实现IHS遥感图像融合的骨架。符号含义r、g、b为同一像素在三个波段上的反射率值v1、v2是变换后的中间色度分量本身不参与替换只在逆变换时用。正变换输出中只有I分量与全色波段在物理含义上对应这也是后续替换操作的理论依据。2.2 三步融合流程正变换、替换、逆变换IHS融合的流程在遥感图像处理教科书中被总结为标准三步几乎所有后续变体都是在这三步上做文章将多光谱影像的三个波段执行RGB到IHS正变换得到I、H、S此处用v1、v2代替H、S对全色波段做直方图匹配或均值方差拉伸使其与I分量具有相同的一阶矩和二阶矩然后用处理后的全色波段替换I分量以替换后的I、原始v1、v2执行逆变换生成融合后的RGB影像其中第二步是关键设计决策。如果直接用原始PAN替换I而不做任何拉伸融合结果会出现整体偏亮或偏暗的问题因为PAN的均值与I分量的均值往往相差很大。常见做法是使用如下线性拉伸pan_stretched (pan - pan.mean()) * (i.std() / pan.std()) i.mean()这个公式将PAN的均值对齐到I分量的均值把PAN的标准差对齐到I分量的标准差。参数含义pan.mean()和pan.std()是全色波段的全局统计量i.mean()和i.std()是I分量的全局统计量。效果是保留PAN的高频空间细节同时让它的亮度水平与MS影像保持一致减少融合后的色偏。对遥感图像融合而言这一步直接决定最终结果的质量。2.3 为什么替换的是I而不是H或S光谱失真的根源从信息论角度看RGB到IHS的变换本质上是坐标旋转。I分量几乎承载了三个波段的公共辐射能量也就是空间结构信息而v1和v2承载的是波段间的差异即色彩信息。因此I分量是高频纹理和亮度信息的理想载体这和PAN影像的物理特性天然匹配。如果试图替换H或S分量结果会非常糟糕H代表色调替换后植被、水体、裸土的颜色会被彻底打乱S代表饱和度替换后整个影像要么灰蒙蒙要么颜色浓到溢出。只有I分量和PAN之间相关性最强替换带来的光谱扭曲最小。但所谓“最小”只是相对而言线性IHS变换本身无法彻底消除色彩畸变这就是后续章节GS变换、小波变换出现的原因——它们本质上都在寻找一个比I分量更优的替换空间。3. 用PythonGDAL在本地跑通IHS融合的最小命令3.1 环境准备与数据读取遥感图像融合的最小环境是Python加上GDAL、NumPy和Rasterio三个库。GDAL负责读取常用遥感影像格式GeoTIFF、ENVI、HDF4/5Rasterio提供更友好的数组接口NumPy承担矩阵运算。pip install gdal rasterio numpy这里不指定版本号因为GDAL的版本与Python版本严格绑定建议在conda环境中用conda install gdal安装预编译版本。数据准备方面需要一对经过配准的影像PAN影像单波段空间分辨率高如0.5mMS影像至少三个波段空间分辨率低如2m两幅影像必须位于同一投影坐标系范围可以略有出入但重叠区要达到95%以上。以下代码读取影像并将MS重采样到PAN的分辨率import rasterio from rasterio.enums import Resampling import numpy as np with rasterio.open(pan.tif) as src: pan src.read(1) pan_profile src.profile with rasterio.open(ms.tif) as src: ms src.read() # 形状: (band_count, height, width) ms_profile src.profile # 将MS重采样到PAN影像的尺寸 with rasterio.open(ms.tif) as src: ms_resampled src.read( out_shape(src.count, pan.shape[0], pan.shape[1]), resamplingResampling.cubic ) print(PAN shape:, pan.shape) print(MS shape after resample:, ms_resampled.shape)参数说明Resampling.cubic表示三次卷积插值比最近邻法nearest更适合光学影像的多光谱重采样因为保留了波段间的光谱关系。out_shape是重采样后的目标尺寸这里直接对齐到PAN的宽高。读取的ms_resampled是一个三维数组维度顺序是波段优先后面所有处理都沿用这个布局。3.2 正变换、全色替换与逆变换的实现核心融合代码可以直接复用第2章的rgb_to_ihs和ihs_to_rgb函数但需要处理一个实际问题遥感影像的DN值范围不是0-255一般8位数据是0-25516位数据是0-65535归一化到浮点型后变换更稳定。# 取前三个波段做RGB融合默认按R-G-B顺序 r ms_resampled[0].astype(np.float64) g ms_resampled[1].astype(np.float64) b ms_resampled[2].astype(np.float64) # 正变换 i, v1, v2 rgb_to_ihs(r, g, b) # 全色波段均值/方差对齐到I分量 pan_f pan.astype(np.float64) pan_stretched (pan_f - pan_f.mean()) * (i.std() / pan_f.std()) i.mean() # 替换I分量 i_new pan_stretched # 逆变换得到融合后的R、G、B r_fused, g_fused, b_fused ihs_to_rgb(i_new, v1, v2) # 裁剪到输入数据的有效值范围 r_fused np.clip(r_fused, 0, np.iinfo(ms.dtype).max) g_fused np.clip(g_fused, 0, np.iinfo(ms.dtype).max) b_fused np.clip(b_fused, 0, np.iinfo(ms.dtype).max)参数说明i.std() / pan_f.std()是标准差比值控制替换后影像的对比度与I分量一致加上i.mean()是亮度对齐。这里的np.clip用于把超出有效值范围的预测值拉回原影像的数据类型范围否则输出16位影像时会出现溢出噪声。将替代后的三个波段拼回数组并写盘# 合并融合波段并写为GeoTIFF fused np.stack([r_fused, g_fused, b_fused], axis0) out_profile ms_profile.copy() out_profile.update({ height: fused.shape[1], width: fused.shape[2], transform: pan_profile[transform] }) with rasterio.open(fused_ihs.tif, w, **out_profile) as dst: dst.write(fused.astype(ms.dtype))这里把输出文件的投影信息与地理变换参数替换为PAN的元数据因为融合后影像的几何分辨率继承自PAN。ms_profile.copy()保留了波段数、压缩格式、数据类型的原始设置transform字段必须更新不然输出的地理坐标会错位。3.3 输出检查先看波谱曲线再看空间纹理融合完成后不能只靠肉眼观察。遥感图像处理的标准质检流程是同时看两类指标空间纹理是否达到全色分辨率放大到像素级别检查道路边缘、屋顶轮廓是否锐利光谱色彩是否接近原MS影像在未变化区域对比融合前后典型地物的光谱曲线from osgeo import gdal ds gdal.Open(fused_ihs.tif) bands ds.RasterCount print(融合影像波段数:, bands) # 输出单波段统计量作为质量记录 for b in range(1, bands 1): band ds.GetRasterBand(b) stats band.ComputeStatistics(False) print(fBand {b}: min{stats[0]:.2f}, max{stats[1]:.2f}, mean{stats[2]:.2f}, std{stats[3]:.2f})波段数通常为3如果你希望保留更多光谱波段如8波段多光谱需要多次执行IHS融合或者改用面向任意波段数的GS变换这一点在第5章展开。检查统计量的目的在于如果某波段的标准差明显高于原始MS影像说明融合过程中发生了严重的光谱拉伸可能需要回头调整直方图匹配的参数。4. IHS融合参数设置与常见坑波段、重采样与直方图匹配4.1 三个波段怎么选真彩色还是假彩色IHS变换一次只能处理三个波段多光谱只要超过3个波段就面临选波段的问题。不同卫星的波段组合习惯不同必须按应用场景决定应用场景推荐波段组合效果地表生态、植被分析R近红外G红B绿植被呈红色适合做植被指数目视解译城市地物分类R红G绿B蓝接近人眼视觉道路与建筑边界清晰水体/湿地监测R近红外G短波红外B红水体呈深黑色湿地边界锐利选择波段时优先选相关性低的三个波段这样v1和v2的方差贡献更大色彩信息保留更充分。如果三个波段高度相关例如红、绿、蓝的原始反射率在植被区接近饱和融合后色彩区分度会显著下降。4.2 重采样方法与分辨率匹配多光谱影像重采样到全色分辨率这一步最容易踩坑。不同插值方法对光谱保持能力的影响如下表nearest速度快但不产生新值容易造成地物边缘锯齿光谱曲线保留最好bilinear线性加权平滑效果好但会轻微模糊纹理cubic三次卷积纹理锐化明显但容易产生超出原值范围的震荡lanczos窗口更大精度更高但计算量成倍增加对于IHS融合推荐使用cubic或lanczos。理由是PAN波段本身分辨率远高于MSMS重采样后要与PAN逐像素对齐插值误差会被融合过程放大。nearest虽然在光谱上最安全但几何对齐误差会直接导致融合结果出现“色彩重影”也就是地物边缘出现一道亮边。另外要注意重采样的边界处理默认情况下GDAL会在影像边缘外推值这会让边缘几列像素出现异常值建议在重采样前先对影像做适当的对称填充。4.3 直方图匹配与光谱失真的三个典型场景直方图匹配是IHS融合中唯一能主动调节光谱保真度的步骤。直接替换I分量在以下三种场景中会引发严重光谱失真场景一PAN与MS成像时间不同步。例如PAN摄取于6月MS摄取于9月两者均值相差悬殊不做对齐会出现整体色偏。场景二PAN包含云影区域。云的亮度极高替换后I分量被抬升IHS逆变换会生成偏白的地物。场景三水体或阴影区域。PAN在这些区域的信噪比极低替换后I分量出现块状噪声逆变换后的色彩完全偏离。针对以上问题常见做法是引入局部直方图匹配。全局匹配不够用就按分块统计直方图。以256×256像元为窗口计算滑动窗口内的均值与方差对PAN做局部对齐后再进行替换可以显著减少局部辐射差异的影响。def local_hist_match(pan, i, window256): from scipy.ndimage import uniform_filter pan_mean uniform_filter(pan, sizewindow, modereflect) i_mean uniform_filter(i, sizewindow, modereflect) pan_std np.sqrt(uniform_filter((pan - pan_mean)**2, sizewindow, modereflect)) i_std np.sqrt(uniform_filter((i - i_mean)**2, sizewindow, modereflect)) # 避免除零 i_std[i_std 1e-6] 1e-6 return (pan - pan_mean) * (i_std / pan_std) i_mean参数说明window是滑动窗口大小值越小局部辐射校正越强但过小会把空间纹理中的真实反差吞掉一般取PAN影像宽度的1/20到1/50之间。modereflect是边界填充方式防止窗口滑到边缘时数据不足。这个函数替换掉第3章中的全局拉伸公式后融合色彩还原度会有肉眼可见的提升。5. 选型对比IHS与Brovey、GS、小波融合怎么选以及一个验证技巧用IHS完成第一版融合后建议顺手跟其他三种方法做一次对比因为它们分别代表了互补的思路。方法替换空间波段数光谱保真计算成本典型场景IHSRGB色彩空间3中等色彩易偏极低快速预览、真彩色合成BroveyRGB比例3较差高亮区易过曝极低单张快出图Gram-SchmidtGS任意统计空间任意波段高低多光谱全波段融合小波变换频率子带任意波段高高高分影像精细融合GS变换是IHS在数学上的推广它用统计方法生成第一分量与I分量类似再做正交化替换因此能处理4波段以上的多光谱数据光谱保真度也更稳。小波变换则把空间信息分解为高频和低频只替换低频部分避免全色波段全盘替换导致的光谱突变。如果你的数据是WorldView或高分一号这种多波段传感器GS是比IHS更合理的基线。一个实用的验证技巧是计算融合前后影像的相关系数CC和相对全局综合误差ERGAS。相关系数衡量光谱保真度ERGAS衡量整体质量公式如下CC 计算原始MS与融合影像对应波段的相关系数逐波段平均 ERGAS 100 * (h/l) * sqrt(mean( (RMSE_b / mean_b)^2 ))其中h/l是PAN与MS分辨率之比。如果ERGAS值小于3说明融合质量很好大于5则需要检查参数设置。计算代码很轻量def ergas(original, fused, ratio4): bands original.shape[0] mse_list [] for b in range(bands): mse np.mean((original[b] - fused[b]) ** 2) mse_list.append(mse / np.mean(original[b]) ** 2) return 100 * ratio * np.sqrt(np.mean(mse_list))把第3章的fused_ihs.tif和重采样后的原始MS影像喂进这个函数快速就能判断该不该继续调参。这个技巧的价值在于遥感图像融合的方法选型不该靠肉眼也不该靠“新方法一定更好”的直觉用定量的ERGAS数据说话比换十个算法都管用。本文还有配套的精品资源点击获取
核心思想与竞赛应用解析)
