简介面向图像融合学习者的TIF变换不变融合算法双版本实现包内含Python 3.8基于OpenCV库与MATLAB两套可直接运行的代码帮助理解变换不变性融合的核心思路适用于希望掌握多源图像特征保留、信息整合技术的初学者及进阶开发者。压缩包共477个文件以jpg测试图像为主471张另有少量tif示例图、py与m格式的算法脚本以及md/txt说明文档整体大小13.32MB结构清晰便于按需取用。已有3373人学习下载。通过运行示例可直观对比两种语言在图像读取、预处理、融合及保存环节的实现差异同时借助内置的多元场景测试图检验算法在不同光照、人物、夜景等图像上的融合效果既能夯实理论基础也能快速迁移到自己的图像处理项目中。 两年前我第一次做TIF影像融合时差点被一个很基础的问题劝退两幅影像的分辨率差了五倍融合出来的图糊成一片。后来我才意识到图像融合算法的代码其实只占工作量的一小半真正麻烦的是TIF格式背后的地理坐标、波段位深和重采样。这篇文章就是基于我后来整理的一套Python和MATLAB版本融合脚本把设计思路、核心算法和踩过的坑都摊开讲适合刚接触遥感影像处理、或者想把手头融合实验快速跑通的同学参考。1. 拿到两幅TIF之后先想清楚融合到底要解决什么问题图像融合不是简单地把两张图叠加而是要解决传感器物理条件带来的信息瓶颈。遥感或者航拍场景里最常见的情况是全色影像空间分辨率高但光谱信息弱多光谱影像光谱信息丰富但空间分辨率低红外影像能暴露温度特征但细节模糊可见光影像纹理清楚却受光照影响大。融合的目标就是用算法把高空间分辨率的结构细节“灌”进高光谱信息的影像里或者把红外特征与可见光纹理整合到一张图里让后续目视判读和自动解译都更可靠。这个目标听起来很直白但落到TIF数据上就复杂了。TIF尤其是GeoTIFF不只是存像素还带了一整套地理参考信息投影坐标系、像元大小、仿射变换六参数、角点坐标等等。融合如果把这些信息弄丢了后面的GIS叠加、矢量勾绘就全乱套。我见过不少人用OpenCV的imread把两幅TIF读进来融合完用imwrite一存看起来是张图放到ArcMap里却“飞”到了海里。原因就是地理参考信息在读写过程中被丢弃了。另外还要先判断融合需求属于哪一类全色锐化Pan-sharpening同源传感器多光谱全色目标是增强空间分辨率。红外与可见光融合异源数据一个看温度特征一个看纹理细节常用于目标检测与监视场景。多聚焦融合同一相机不同焦距拍摄主要用于显微和摄影领域工程上相对少见。判断清楚场景再选算法因为不同场景对“光谱保真”的要求完全不一样。全色锐化如果光谱扭曲地物分类就会出问题红外与可见光融合如果细节取舍不当目标就被淹没在背景里。这篇文章的代码以拉普拉斯金字塔融合和多尺度分解为主因为它对这几种场景的适配性最好原理也直观便于读者改成自己的版本。2. 融合前的数据对齐分辨率、范围和波段位深都是暗坑很多人以为融合算法是核心结果一跑代码就报维度不匹配或者融合出来半边黑半边亮。这些绝大多数不是算法问题而是预处理没做干净。我自己的处理流程里预处理占的时间比写融合函数还长。2.1 分辨率重采样先统一到同一个像元尺寸两幅TIF的分辨率不一致是最常见的情况。比如多光谱影像是15米分辨率全色影像是2米分辨率直接做像素级融合矩阵shape完全不同。常规做法是把高分辨率影像作为基准用最近邻、双线性或三次卷积把低分辨率影像重采样到高分辨率网格上。用Python的rasterio处理时看一下两个文件的shapeimport rasterio with rasterio.open(multispectral.tif) as src: ms src.read() ms_profile src.profile print(MS shape:, ms.shape, bounds:, src.bounds) with rasterio.open(panchromatic.tif) as src: pan src.read() pan_profile src.profile print(PAN shape:, pan.shape, bounds:, src.bounds)如果形状不一致可以用rasterio.warp.reproject重采样import numpy as np from rasterio.warp import reproject, Resampling # 把多光谱数据重采样到全色影像的空间范围 dst_transform pan_profile[transform] dst_height pan_profile[height] dst_width pan_profile[width] ms_resampled np.empty((ms.shape[0], dst_height, dst_width), dtypenp.float32) reproject( sourcems.astype(np.float32), destinationms_resampled, src_transformms_profile[transform], src_crsms_profile[crs], dst_transformdst_transform, dst_crspan_profile[crs], resamplingResampling.bilinear )MATLAB这边用imresize也可以但没有原生的经纬度感知能力需要先手动读取地理信息再处理。这是我后来坚持用Python做预处理的其中一个原因MATLAB更适合做算法验证。2.2 波段位深别在归一化之前就丢信息TIF影像常见的有8位和16位两种位深。很多遥感产品是16位整型动态范围非常大。如果一开始就astype(np.uint8)高光溢出、暗部细节丢失基本不可避免融合效果再好也救不回来。稳妥的做法是先把数据转成float32再做逐波段归一化def normalize_band(band): band band.astype(np.float32) vmin, vmax np.percentile(band, 2), np.percentile(band, 98) band np.clip((band - vmin) / (vmax - vmin 1e-8), 0, 1) return band用2%和98%分位数做截断而不是直接min/max可以避免个别极亮像元把整个影像压暗。这个细节我是在处理夜间红外影像时踩坑后加上的前后效果差异非常明显。2.3 无效值区域和地理范围裁剪许多TIF边缘含有NoData值如果不处理融合时这些NoData会被当成0值参与计算产出大量异常暗色区域。建议在读取时获取nodata值并做mask# 读取nodata值 nodata src.nodata if nodata is not None: mask (band nodata) band[mask] 0同时要保证两幅影像覆盖的地理范围一致。最简单的方式是统一用目标影像的bounds做窗口裁剪或者在重采样时把目标网格直接对齐。范围不一致导致融合结果边缘有偏移这个问题在后期精度评价时很难解释清楚。3. 核心融合算法怎么选IHS、PCA和多尺度分解的适用边界图像融合算法多到能写一本书但工程上真正高频使用的其实就那么几类。选算法之前先理解每种方法在做什么才有底气根据数据去调参。3.1 IHS变换和PCA替换经典全色锐化思路IHS变换的思路很直观把RGB三波段转到亮度、色相、饱和度空间然后用高分辨率全色影像替换亮度分量再反变换回RGB。好处是计算量小视觉效果强烈坏处是光谱保真度一般尤其当地物颜色很丰富时彩色畸变明显。PCA主成分分析替换稍微高级一点。对多光谱波段做PCA第一主成分通常是信息量最大、也最接近全色影像灰度分布的成分用全色影像替换它后再反变换。这个方法在很多商业遥感软件里是默认的全色锐化方案。但它同样对波段间相关性有要求异源数据用了容易出伪影。3.2 多尺度分解融合质量更稳的主流方案多尺度分解的核心思想是把影像拆成低频基础层和高频细节层。低频包含整体辐射信息高频包含边缘、纹理和结构。融合时低频做加权平均或能量保持高频选择细节更丰富的分量最后重建。这个思路在红外与可见光融合里表现得尤其好。红外影像低频能反映热分布可见光高频能保留轮廓。拉普拉斯金字塔是小波变换的一种雏形实现简单效果可靠小波变换进一步引入了方向选择性红外与可见光边缘更容易被保留。下面我给的示例代码都以拉普拉斯金字塔为主因为它的逻辑最短读者改成DWT也不难。3.3 方法对比不同场景下怎么选方法适用场景优点缺点IHS变换全色锐化、快速预览实现简单速度快光谱畸变明显PCA替换多光谱与全色信息量集中通用性好无法处理异源数据差异拉普拉斯金字塔红外与可见光、通用融合结构清晰光谱保真较好对配准误差较敏感小波变换红外与可见光、多聚焦方向信息保留强参数选择复杂度高深度学习类特定场景专项优化效果上限高需要训练数据工程成本大如果只是做实验验证拉普拉斯金字塔基本不会错。如果做生产级全色锐化我建议PCA和金字塔都跑一遍用定量指标选最优而不是拍脑袋定算法。4. Python版用rasterio numpy实现拉普拉斯金字塔融合Python版本我选择rasterio numpy的组合。rasterio负责读写TIF并保留地理参考numpy负责矩阵运算。算法部分的核心是金字塔分解、融合策略和金字塔重建。4.1 金字塔分解函数高斯金字塔的每一层由上一层做高斯模糊下采样得到拉普拉斯金字塔则保存高斯金字塔相邻两层的差分信息。因为差分值包含了细节纹理融合时可以根据细节强度自适应选取。import numpy as np from scipy.ndimage import gaussian_filter from skimage.transform import resize def gaussian_pyramid(img, levels): pyramid [img] current img for _ in range(levels - 1): current gaussian_filter(current, sigma1.0) current current[::2, ::2] pyramid.append(current) return pyramid def laplacian_pyramid(gauss_pyr): lap_pyr [] for i in range(len(gauss_pyr) - 1): size (gauss_pyr[i].shape[0], gauss_pyr[i].shape[1]) expanded resize(gauss_pyr[i 1], size, modereflect) lap_pyr.append(gauss_pyr[i] - expanded) lap_pyr.append(gauss_pyr[-1]) return lap_pyr def reconstruct(lap_pyr): current lap_pyr[-1] for i in range(len(lap_pyr) - 2, -1, -1): size (lap_pyr[i].shape[0], lap_pyr[i].shape[1]) current resize(current, size, modereflect) lap_pyr[i] return current这个实现里levels一般取3到5。层数太少细节分离不充分层数太多底层信息过于平滑融合结果容易发虚。4.2 融合策略低频取能量高频取最大融合策略是整个算法的灵魂。我的经验是低频层用加权平均权重根据两幅影像的全局亮度统计来定高频层用绝对值取大策略因为细节信号的强弱直接反映边缘清晰度。def fuse_pyramids(lap1, lap2, weight0.5): fused [] for i in range(len(lap1)): if i len(lap1) - 1: # 基础层加权平均 fused.append(weight * lap1[i] (1 - weight) * lap2[i]) else: # 细节层绝对值取大 mask np.abs(lap1[i]) np.abs(lap2[i]) fused_layer np.where(mask, lap1[i], lap2[i]) fused.append(fused_layer) return fused这里有个细节值得注意如果两幅影像的平均辐射水平差异太大加权平均会导致整体亮度产生跳变。我在实际项目中会先对两幅影像做直方图匹配把亮度分布拉齐再做融合效果会稳很多。4.3 主流程与TIF写出最后是完整主流程的代码。读取两幅TIF先重采样对齐再归一化逐波段做金字塔融合最后把地理参考写入新TIF。import rasterio from rasterio.transform import from_origin def fuse_tif(ms_path, pan_path, output_path, levels4): with rasterio.open(pan_path) as src: pan src.read(1).astype(np.float32) profile src.profile transform src.transform crs src.crs with rasterio.open(ms_path) as src: ms src.read().astype(np.float32) ms_transform src.transform ms_crs src.crs # 将ms重采样到pan的网格 from rasterio.warp import reproject, Resampling resampled np.empty((ms.shape[0], pan.shape[0], pan.shape[1]), dtypenp.float32) for b in range(ms.shape[0]): reproject( sourcems[b], destinationresampled[b], src_transformms_transform, src_crsms_crs, dst_transformtransform, dst_crscrs, resamplingResampling.bilinear ) fused_bands [] for b in range(resampled.shape[0]): band1 normalize_band(resampled[b]) band2 normalize_band(pan) gp1 gaussian_pyramid(band1, levels) gp2 gaussian_pyramid(band2, levels) lap1 laplacian_pyramid(gp1) lap2 laplacian_pyramid(gp2) fused_lap fuse_pyramids(lap1, lap2, weight0.6) fused_bands.append(reconstruct(fused_lap)) fused np.stack(fused_bands, axis0) fused np.clip(fused * 255, 0, 255).astype(np.uint8) profile.update(dtyperasterio.uint8, countfused.shape[0], transformtransform, crscrs) with rasterio.open(output_path, w, **profile) as dst: dst.write(fused)输出路径里建议避免在文件名中出现中文和特殊字符否则部分GIS软件读取时会出编码问题。另外写TIF前记得用profile.update保留原来影像的driver、压缩选项等信息。5. MATLAB版适合快速验证的融合脚本写法MATLAB版本更适合快速原型验证不需要管文件路径里的地理信息细节内置的矩阵运算和图像处理工具箱让算法改写非常方便。我通常用它来验证新融合策略验证完再翻译到Python做生产。5.1 用impyramid做高斯金字塔MATLAB的impyramid函数直接支持高斯金字塔的分解和重建。reduce对应降采样expand对应放大。直接用这两个函数就能搭出拉普拉斯金字塔。function fused fuse_tif_matlab(ms_bands, pan_band, levels) % ms_bands: H x W x C多光谱波段 % pan_band: H x W全色波段 fused zeros(size(ms_bands)); for c 1:size(ms_bands, 3) band1 double(ms_bands(:, :, c)); band2 double(pan_band); band1 (band1 - min(band1(:))) / (max(band1(:)) - min(band1(:))); band2 (band2 - min(band2(:))) / (max(band2(:)) - min(band2(:))); fused(:, :, c) pyramid_fuse(band1, band2, levels); end end function fused pyramid_fuse(img1, img2, levels) pyr1 cell(levels, 1); pyr2 cell(levels, 1); curr1 img1; curr2 img2; for i 1:levels pyr1{i} curr1; pyr2{i} curr2; curr1 impyramid(curr1, reduce); curr2 impyramid(curr2, reduce); end % 从最小层开始重建 fused (pyr1{levels} pyr2{levels}) / 2; for i levels - 1 : -1 : 1 lap1 pyr1{i} - imresize(impyramid(pyr1{i}, reduce), size(pyr1{i})); lap2 pyr2{i} - imresize(impyramid(pyr2{i}, reduce), size(pyr2{i})); fused imresize(fused, size(pyr1{i})) max(abs(lap1), abs(lap2)) .* sign(lap1 lap2); end end注意这里重建的时候细节层取绝对值更大的一方。sign(lap1 lap2)是为了让选出来的细节保持原方向避免两张都取最大值之后出现边缘方向反转的假纹理。这个细节是我对比了好多组实验后加上的对边缘质量影响很大。5.2 TIF读取与写出保留地理信息MATLAB的geotiffread和geotiffwrite可以保留TIF的地图坐标信息。如果只用imread和imwrite地理参考会被丢弃融合结果进不了GIS流程。% 读取 [X, R] geotiffread(multispectral.tif); [Y, R2] geotiffread(panchromatic.tif); % 假设需要对齐时可以用 georefcells 对 X 进行重采样 % 融合 fused fuse_tif_matlab(X, Y, 4); % 写出保留坐标系信息 geotiffwrite(fused_output.tif, uint8(fused * 255), R);MATLAB的geotiffwrite要求传入的数据范围匹配栅格参考对象。如果源TIF是16位整型建议在写出前将融合结果缩放到0到65535再转uint16否则亮部会普遍过曝。5.3 用批处理脚本加快参数调试调试阶段经常要遍历不同的金字塔层数和融合权重。我习惯把参数提成结构体写外层循环批量跑results struct(); levels_list [3, 4, 5]; weights_list [0.4, 0.5, 0.6]; idx 1; for lv levels_list for w weights_list fused fuse_tif_matlab(X, Y, lv); results(idx).level lv; results(idx).weight w; results(idx).image fused; idx idx 1; end end然后一次性计算所有结果的客观指标选最优组合。这个过程在MATLAB里写起来最快比Python反复开文件省时间。6. 两个版本我都踩过的坑零值边界、浮点误差和数据范围裁剪这部分是整篇文章最想分享的实操经验。很多问题看起来像算法不行实际是数据处理细节没处理好。6.1 零值边界会导致融合影像整体发黑TIF边缘经常有大块无效区域。归一化之后这些区域变成0高频细节层在这些位置会产生强烈的负差分融进结果后就形成一圈黑边。我当时调试红外与可见光融合时目标的轮廓没出来黑边倒是很显眼。解决办法有两种一是用形态学操作把无效区扩大一圈生成一个可靠的融合掩膜二是融合完成后对无效区做中值模糊修复。第一种更干净from scipy.ndimage import binary_dilation valid1 ~(band1 0) valid2 ~(band2 0) valid_mask valid1 valid2 valid_mask binary_dilation(valid_mask, iterations5) # 融合完成后将无效区置为原影像或0 fused np.where(valid_mask, fused, 0)但保险起见先用NoData信息生成mask再做融合效率和质量都比事后修复高。6.2 浮点误差会让金字塔重建出现负值拉普拉斯金字塔分解时涉及大量差值和resize操作重建时不可避免产生轻微负值。如果直接把负值截成0暗部细节会损失。我习惯在最终写TIF之前做np.clip(fused, 0, 1)却还遇到过整体偏暗的问题后来发现是金字塔重建过程中能量没有完全恢复。排查下来问题出在resize的插值方式与金字塔分解时的下采样方式不一致。解决方式很粗暴但有效重建完成后做一次直方图匹配把融合结果的分布对齐到原多光谱影像的分布def hist_match(source, template): src_sort np.sort(source.ravel()) tmpl_sort np.sort(template.ravel()) src_cdf np.cumsum(src_sort) / src_sort.sum() tmpl_cdf np.cumsum(tmpl_sort) / tmpl_sort.sum() map_values np.interp(source.ravel(), src_sort, tmpl_sort) return map_values.reshape(source.shape)这个方法能在不改变纹理结构的前提下把整体亮度和色彩拉回正常范围。6.3 大文件融合时内存不够和速度过慢的优化一景10厘米分辨率的无人机TIF动辄上亿像素单波段float32就要几百MB。如果一次读入整个影像再堆叠多个金字塔层内存立刻爆掉。我常用的优化策略是分块处理。把图像切成若干个重叠瓦片每个瓦片独立融合再按权重拼回去。重叠区域通常取金字塔最大模糊半径的两倍以上避免瓦片边界可见。另外金字塔层数不一定要从头到尾全算可以只计算高频层底层用两张图的加权平均直接作为基础层能省下不少内存。MATLAB版本同样面临内存压力。如果workers数量够可以并行处理各个波段的融合parfor c 1:size(ms_bands, 3) fused(:, :, c) pyramid_fuse(ms_bands(:, :, c), pan_band, levels); end需要先parpool开启并行池否则parfor就是普通循环。7. 效果验证和批量处理该怎么落地跑通融合流程只是第一步。实际项目中必须用客观指标验证融合效果不能只靠眼睛看。我常用的定量指标有三个熵信息量、平均梯度细节清晰度和光谱偏差与原多光谱影像的差异。三个指标要一起看因为它们之间经常互相冲突——细节提升了光谱却偏了。下面是Python里计算平均梯度和光谱相关系数的函数def average_gradient(image): gx np.diff(image, axis1) gy np.diff(image, axis0) return np.sqrt((gx ** 2 gy ** 2).mean()) def spectral_correlation(fused, original_ms): corr_sum 0 for b in range(fused.shape[0]): f fused[b].ravel() o original_ms[b].ravel() corr_sum np.corrcoef(f, o)[0, 1] return corr_sum / fused.shape[0]评价时至少选择一块地物丰富的区域和一块平坦区域分别计算才能代表整体效果。只用全局指标容易被大面积均质区域拉低差异。批量处理时先把单景融合封装成函数再遍历文件夹。需要注意输出路径的组织方式我一般按“输入数据日期/融合方法/参数”建目录方便后续效果回溯。融合实验会跑很多轮如果不保存参数三个月后回头根本想不起某张结果图是用哪组参数生成的。建议在输出TIF的文件里顺手写入融合参数作为metadatawith rasterio.open(output_path, w, **profile) as dst: dst.update_tags(fusion_methodlaplacian_pyramid, levelsstr(levels), weightstr(weight)) dst.write(fused)这个习惯救了我很多次尤其在给甲方交付时能直接解释清楚每张成果图是怎么来的。最后再分享一个小技巧融合算法的参数往往跟传感器特性强相关。同一套参数换了一台无人机、换了一个卫星传感器效果就可能明显变差。所以每次拿到新数据先用一小块代表性区域跑参数扫描确定最优组合后再全图跑能节省大量时间。本文还有配套的精品资源点击获取