遥感图像融合TIF算法：Python与MATLAB实战指南

简介：图像融合TIF算法（Transform Invariant Fusion）提供Python与MATLAB两种语言的实现代码，适合图像处理初学者与进阶开发者学习。该算法基于变换不变性设计，能在图像发生平移、缩放或旋转时仍保持稳定融合效果，有效保留边缘、纹理等关键特征，减少信息损失。压缩包共477个文件，主体为471张JPG测试图像，另有2张TIF原始图、Python脚本、MATLAB脚本、Markdown说明文档及TXT注释，容量仅13.32MB，结构清晰便于查阅。Python版本基于OpenCV编写，MATLAB版本使用内置图像处理函数，两者均附带可直接运行的测试脚本，读者可对比输出结果，掌握不同编程环境下的算法实现技巧。该资源已有3373人学习，适用于课堂实验、毕业设计或课题预研，可直接运行代码理解融合流程，也可基于现有实现进行二次开发或算法改进，是图像融合方向不可多得的参考资料。 做遥感图像处理的朋友，对“图像融合TIF”这几个字应该都不陌生。我最早接触这个需求，是拿到一组国产卫星的多光谱和全色波段数据，老板要求出2米分辨率的融合影像，而手里的工具只有ArcGIS，跑一次流程看着挺简单，但批次多了之后慢得让人抓狂。后来我干脆把算法搬到了Python和MATLAB里，自己控制流程，想怎么调参数就怎么调。这一篇文章就围绕“图像融合TIF算法”的Python和MATLAB版本代码展开，讲讲我是怎么选方案、写代码、踩坑、优化的，适合正在做遥感影像预处理的同学、做毕业设计的学生，以及想把融合流程自动化落地的工程师参考。

这套代码解决的核心问题很明确：把高分辨率全色影像的空间细节和低分辨率多光谱影像的光谱信息合到一张TIF里，输出一个兼顾纹理和色彩的成果。放到实际场景中，就是制作正射影像底图、变化检测前处理、目标识别前的数据融合。全文不搞花架子，直接给可跑的代码和参数思路。

1. 需求拆解与方案选型

1.1 为什么图像融合偏偏绕不开TIF

因为TIF（尤其是GeoTIFF）几乎是遥感行业默认的交付格式。相比JPG、PNG，它可以无损保存16位甚至32位辐射分辨率，支持多波段（R、G、B、NIR等），还能在文件头里写入地理坐标、投影信息、仿射变换矩阵这些关键元数据。换句话说，一张TIF不只是“一张图”，它还是一份自带空间参考的地理数据。用Python或MATLAB做融合时，如果读进来的是普通图像格式，地理信息很可能就丢了；处理完再想落回GIS平台里叠加、出图，还得重新做地理配准，非常麻烦。所以我在代码里优先使用能够保留GeoTIFF元数据的读写方式，这也是标题里强调TIF的原因。

1.2 Python和MATLAB怎么选

这个选择没有绝对的对错，我自己的习惯是看场景分：

• Python版强的点在于免费、生态好、适合批量处理和工程部署。配合rasterio、GDAL、numpy这几个库，能把TIF读写成带地理信息的数组，融合算法部分可以透明地看到每一步计算，跟深度学习、机器学习结合也很方便。
• MATLAB版强在矩阵运算语法简洁、内置的Image Processing Toolbox和Wavelet Toolbox好用，做原型验证、课堂实验、算法对比时非常顺手。

我的做法是：原型快速验证用MATLAB，正式批量处理和交付用Python。两者都要写，因为同一个算法在不同语言下的边界处理、数据类型转换、文件读写细节差异很大，能避开很多隐藏的坑。下面我分别讲两套代码的具体实现。

2. 核心算法原理与选型

2.1 像素级融合：从加权平均说起

图像融合分三个层次：像素级、特征级、决策级。做遥感TIF融合，最常用的是像素级，因为它保留的原始信息最多，实现也直观。像素级里又有很多算法，最简单的就是加权平均：

[ F(i,j) = w_1 \times A(i,j) + w_2 \times B(i,j) ]

其中 (w_1 + w_2 = 1)，实际操作里很多初学者直接取0.5和0.5，结果往往是图像发灰、对比度下降。原因很简单：多光谱影像和全色影像的灰度范围、均值、标准差都不一样，直接加权等于把两堆不同分布的数据混在一起，自然会把原本的清晰层次“稀释”掉。所以我在实际代码里，会先对两幅影像做直方图匹配，让全色影像的均值和方差对齐到多光谱影像，再做融合，效果会好很多。

2.2 PCA、小波与拉普拉斯金字塔

加权平均只是基础，生产环境里我更偏爱两类算法：

• PCA（主成分分析）替换法：核心思想是对多光谱波段做主成分分析，第一主成分（PC1）集中了最大的空间方差，把这个分量用高分辨率全色影像替换，然后再反变换回RGB等波段。好处是光谱保真相对好，缺点是线性变换，对强非线性地物会有失真。
• 多尺度分解（拉普拉斯金字塔、小波变换）：把图像分解成低频和高频，低频代表光谱信息，高频代表空间细节。全色影像的高频细节替换或增强到多光谱影像对应频段里去。这类方法效果好，但计算量略大。

下面两章就是这两种思路的具体实现，我分别给Python和MATLAB的代码。

3. Python版本实现

3.1 环境准备与依赖安装

我的Python环境用的是3.9，建议至少3.8以上。依赖库就那么几个：

pip install rasterio numpy opencv-python pywavelets

如果机器上没有GDAL，也可以另外装，但rasterio已经自带打包了大部分GDAL功能，日常读写TIF够用。至于环境装不上、GDAL报错的问题，可以参考对应操作系统的预编译轮子来解决。我建议用conda创建虚拟环境，因为conda对geospatial库的依赖管理比pip省心：

conda create -n imagefusion python=3.9 conda install -c conda-forge rasterio geopandas conda install numpy opencv pywavelets

3.2 核心代码实现与参数说明

3.2.1 TIF读取与数据检查

我先封装一个读取函数，返回影像数组和元数据：

import numpy as np import rasterio from rasterio.transform import from_origin def read_tif(path): """读取GeoTIFF，返回三维数组和元信息字典。""" with rasterio.open(path) as src: data = src.read() # 形状为 (波段数, 高度, 宽度) meta = src.meta.copy() return data, meta def check_same_size(data1, data2): """检查两幅影像尺寸是否一致，不一致先报警。""" if data1.shape != data2.shape: raise ValueError(f"影像尺寸不一致：{data1.shape} vs {data2.shape}")

这一段看着简单，但实际项目里尺寸不一致的情况特别常见。多光谱波段数可能是4个，全色波段只有1个，如果直接做逐波段的矩阵运算会立刻报错。我一般会用OpenCV的resize或rasterio的reproject先把全色影像裁剪或重采样到跟多光谱一样的大小，再做融合。

3.2.2 直方图匹配

为了让全色影像与多光谱影像的灰度分布对齐，我用了一个实用的直方图匹配方法：

def histogram_matching(source, target): """ 把source影像的灰度分布匹配到target的灰度分布。 输入输出都是二维数组（单波段）。 """ source = source.astype(np.float32) target = target.astype(np.float32) # 计算分位数映射 src_vals = np.sort(source, axis=None) tgt_vals = np.sort(target, axis=None) # 构建LUT src_idx = np.linspace(0, len(src_vals) - 1, 256).astype(np.int32) tgt_idx = np.linspace(0, len(tgt_vals) - 1, 256).astype(np.int32) src_quantiles = src_vals[src_idx] tgt_quantiles = tgt_vals[tgt_idx] # 线性插值映射 lut = np.interp(source, src_quantiles, tgt_quantiles) return lut.astype(np.float32)

直方图匹配的原理不提太深，核心就是把原始影像的累计分布函数映射到目标影像的累计分布函数上。好处是融合前两幅影像已经有相近的亮度分布，后续加权平均不会把影像搞灰。

3.2.3 加权平均融合（带防溢出）

def weighted_fusion(img1, img2, w1=0.5, w2=0.5, meta=None, out_path=None): """两幅影像加权融合，支持3D数组或2D数组，自动按波段处理。""" if img1.dtype != img2.dtype: # 统一转float32避免溢出 img1 = img1.astype(np.float32) img2 = img2.astype(np.float32) fused = w1 * img1 + w2 * img2 # 防止超出原数据范围的溢出 if np.issubdtype(img1.dtype, np.integer): fused = np.clip(fused, np.iinfo(img1.dtype).min, np.iinfo(img1.dtype).max).astype(img1.dtype) else: fused = fused.astype(img1.dtype) if out_path: write_tif(out_path, fused, meta) return fused def write_tif(out_path, data, meta): """保留地理参考写出TIF，注意dtype和波段数变更。""" with rasterio.open(out_path, 'w', **meta) as dst: dst.write(data)

这里特别说明一下，rasterio的meta里包含dtype、width、height、count、crs、transform等关键信息。融合完后如果dtype变了（比如uint16变成float32），必须手动更新meta里的dtype再写，否则rasterio会报错或输出无效文件。这也是很多新人经常踩的坑。

3.2.4 拉普拉斯金字塔融合

加权平均简单但效果有限，我实际生产里更常用拉普拉斯金字塔。它的思路讲起来不复杂：先把图像一层层高斯降采样，再用下一层与上一层插值放大的差值构建拉普拉斯层，在每一层分别做融合，再反过来逐层重建。这样做能有效保留高频细节，也没有简单的加权平均那么容易发灰。

def gaussian_pyramid(img, levels): """生成高斯金字塔，levels越大塔越高。""" gp = [img.astype(np.float32)] for _ in range(levels): img = cv2.pyrDown(gp[-1]) gp.append(img) return gp def laplacian_pyramid(img, levels): """由高斯金字塔推导拉普拉斯金字塔。""" gp = gaussian_pyramid(img, levels) lp = [] for i in range(len(gp) - 1, 0, -1): size = (gp[i-1].shape[1], gp[i-1].shape[0]) up = cv2.pyrUp(gp[i], dstsize=size) lp.append(gp[i-1] - up) return lp def pyramid_fusion(img1, img2, levels=3): """ 对二维单波段影像做拉普拉斯金字塔融合。 多波段可以循环调用本函数。 """ lp1 = laplacian_pyramid(img1, levels) lp2 = laplacian_pyramid(img2, levels) # 高频层取绝对值大的那种，低频层取平均 fused_lp = [] for l1, l2 in zip(lp1, lp2): fused_lp.append(np.where(np.abs(l1) > np.abs(l2), l1, l2)) # 取最后一层低频做平均 base1 = gaussian_pyramid(img1, levels)[-1] base2 = gaussian_pyramid(img2, levels)[-1] base = (base1 + base2) / 2.0 # 重建 result = base for layer in reversed(fused_lp): size = (layer.shape[1], layer.shape[0]) result = cv2.pyrUp(result, dstsize=size) + layer return result

这个实现的“高频绝对值取大”逻辑，实际上就是选择更加锐利、细节更丰富的边缘信息，融合后的道路、建筑边界会明显清晰很多。低频取平均则保留双方的光谱整体趋势，不会出现色偏。多波段影像就拆成单波段循环处理，最后再stack回去。

3.2.5 小波融合（备选）

如果安装了PyWavelets，小波融合其实代码也简单：

import pywt def wavelet_fusion_band(band1, band2, level=2, wavelet='db2'): coeffs1 = pywt.wavedec2(band1, wavelet, level=level) coeffs2 = pywt.wavedec2(band2, wavelet, level=level) fused_coeffs = [] for c1, c2 in zip(coeffs1, coeffs2): if isinstance(c1, np.ndarray): # 低频近似分量取平均 fused_coeffs.append((c1 + c2) / 2.0) else: # 高频细节分量取模极大值 fused_detail = [] for d1, d2 in zip(c1, c2): fused_detail.append(tuple(np.where(np.abs(a) > np.abs(b), a, b) for a, b in zip(d1, d2))) fused_coeffs.append(tuple(fused_detail)) fused = pywt.waverec2(fused_coeffs, wavelet) return fused

小波融合的好处是能保留更多中高频信息，适用于地物纹理复杂区域；缺点是对小波基和分解层数比较敏感，batch处理时要多测几组参数。我一般固定用db2或db4，层数选2或3，层数太高会导致光谱失真。

3.3 Python版踩坑记录

• 内存爆炸：TIF文件太大时，直接整幅读入内存会爆。比如一张4波段、2万×2万的影像，光numpy数组就往几个GB去了。我的解决思路是分块（block）读取，rasterio支持window参数，按256×256或512×512的窗口循环读出、融合、写回，再把地理变换信息补齐。相关热搜里的“gis tif 文件太大”就是这么来的，文件大不是无解的。
• 类型溢出：uint16虽然能存0到65535，但两幅影像加权相加时如果转成了uint8或者中途中转的numpy数组类型不对，会出现截断、溢出、颜色条带。处理前统一转float32，最后再转回原数据类型，这是最稳的。
• 地理坐标丢失：只调cv2.imread去读TIF一定会丢地理信息。老老实实用rasterio或GDAL，融合结果才能直接拖回ArcMap里定位。

4. MATLAB版本实现

4.1 环境准备

MATLAB做图像融合相对省心，官方工具箱已经把大部分函数封装好了。我常用的工具箱：

• Image Processing Toolbox：提供imread、imwrite、impyramid、imresize等基础操作
• Mapping Toolbox：提供geotiffread、geotiffwrite
• Wavelet Toolbox（可选）：提供wavedec2、waverec2

版本上我最早用的是R2019b，后面换到R2023a，代码基本没大改，兼容性还是不错的。如果手头没有Mapping Toolbox授权，也可以用imread读取TIF，但那样地理坐标信息会丢，写回时还得自己拼坐标，所以我还是建议装上Mapping Toolbox。

4.2 核心代码实现

4.2.1 GeoTIFF读写

function fused_img = matlab_pixel_fusion(path1, path2, out_path) % 读取两幅带地理信息的TIF [img1, R1] = geotiffread(path1); [img2, R2] = geotiffread(path2); % 转为double计算，避免溢出 img1 = double(img1); img2 = double(img2); % 尺寸检查 if any(size(img1) ~= size(img2)) error('两幅影像尺寸不一致，请先配准或重采样'); end % 简单加权融合 w1 = 0.5; w2 = 0.5; fused_img = w1 * img1 + w2 * img2; % 转回uint16 fused_img = uint16(round(fused_img)); % 写入带地理参考的TIF geotiffwrite(out_path, fused_img, R1); end

这里有一个细节：geotiffwrite的输出文件名后缀必须是.tif，否则会报格式不支持。另外，两张TIF的投影坐标系最好一致，R1、R2如果不一致，输出会混乱，需要提前用projfwd或mapproj转换。不过大多数业务场景里，同一批影像的坐标系都一样，所以偶尔才需要处理投影转换的问题。

4.2.2 PCA替换融合

MATLAB做PCA有现成函数pca或princomp，我常用很短的代码实现：

function fused = matlab_pca_fusion(pan, ms) % pan: 全色波段 单波段二维数组 % ms: 多光谱影像 H×W×B [H, W, B] = size(ms); ms_2d = reshape(ms, H*W, B); % 像素×波段 ms_2d = ms_2d'; % 波段×像素 % PCA [coeff, score, ~] = pca(ms_2d'); % score: N×B, 第一列为PC1 % 取PC1并转为影像形状 pc1 = score(:, 1); pc1_img = reshape(pc1, H, W); % 直方图匹配全色到PC1 pan = double(pan); pan_std = (pan - mean(pan(:))) / std(pan(:)) * std(pc1_img(:)) + mean(pc1_img(:)); % 替换score的第一列 score_new = score; score_new(:, 1) = pan_std(:); % 逆变换回波段域 ms_new = score_new * coeff'; fused = reshape(ms_new, H, W, B); fused = uint16(round(fused)); end

PCA融合的核心不在于代码多复杂，而在于理解：多光谱波段之间有很强的相关性，PCA把主要空间结构压缩进PC1，用高分辨率全色影像替换PC1再反变换回原始波段空间，本质上是用高分辨率空间信息去“调制”原有多光谱的光谱分布。这也是红外与可见光图像融合中常见的一类思路。

4.2.3 拉普拉斯金字塔融合

MATLAB里做金字塔融合实在是方便，因为buildpyr工具函数都帮你封装好了：

function fused = matlab_laplacian_fusion(img1, img2, levels) % 使用拉普拉斯金字塔融合 if nargin < 3 levels = 3; end % 转double img1 = double(img1); img2 = double(img2); if size(img1, 3) == 3 % 多波段分别处理 fused = zeros(size(img1)); for b = 1:3 fused(:, :, b) = fuse_single_band(img1(:, :, b), img2(:, :, b), levels); end else fused = fuse_single_band(img1, img2, levels); end fused = uint16(round(fused)); end function out = fuse_single_band(im1, im2, levels) g1 = im1; g2 = im2; % 构建高斯金字塔 for k = 1:levels p1{k} = g1; p2{k} = g2; g1 = impyramid(g1, 'reduce'); g2 = impyramid(g2, 'reduce'); end % 拉普拉斯层 l1 = cell(1, levels); l2 = cell(1, levels); for k = 1:levels e1 = impyramid(p1{k+1}, 'expand'); e2 = impyramid(p2{k+1}, 'expand'); % 尺寸可能不严格对齐，需要裁剪 [r, c] = size(p1{k}); e1 = e1(1:r, 1:c); e2 = e2(1:r, 1:c); l1{k} = p1{k} - e1; l2{k} = p2{k} - e2; end % 顶层低频取平均 top = (g1 + g2) / 2; % 逐层重建 out = top; for k = levels:-1:1 out = impyramid(out, 'expand'); sz = size(l1{k}); out = out(1:sz(1), 1:sz(2)); fused_layer = max(abs(l1{k}), abs(l2{k})); % 高频取最强 out = out + fused_layer; end end

这个实现跟Python版的高频取大是同一个思路，但MATLAB的impyramid内置了高斯滤波和样本抽取，代码量会少一些。我测试过，两张2万×2万的TIF在这个流程下，用时基本在30秒以内，取决于CPU内存带宽。

4.3 MATLAB的边界处理

MATLAB里impyramid在expand时，尺寸会跟上一层差1个像素或更少，所以代码里我做了尺寸裁剪。这个细节如果不处理，矩阵维度对不上会直接报错。另一个坑是geotiffwrite要求输入数据类型必须与原始TIF一致，如果原来是uint16，你给个int16它会提示错误。因此转uint16这一步是必需的。

5. 常见问题与排查技巧

5.1 TIF读取报错

排查思路整理成一张速查表：

5.2 融合后图像发灰

这算是最常见的问题了。我见过很多人融合完一看，图像整体灰蒙蒙的，就像蒙了一层雾。原因主要有两个：

• 没做直方图匹配就直接加权
• 融合后没有做反差拉伸或线性增强

处理技巧是：融合输出的TIF在ArcGIS里如果显示发灰，不要急着改算法，先用拉伸方式看看（比如拉伸到99%）。如果拉伸后效果不错，说明是显示问题，不是算法问题。如果拉伸后依然发灰，再回去检查直方图匹配是否做了、权重分配是否合理。真正的生产流程里，我会在保存前加一个2%线性拉伸：

def linear_stretch(data, low=2, high=98): """按百分比截断做线性拉伸。""" for band in range(data.shape[0]): lo, hi = np.percentile(data[band], [low, high]) data[band] = np.clip((data[band] - lo) / (hi - lo + 1e-8) * 65535.0, 0, 65535) return data

5.3 大文件内存溢出处理

图像融合任务经常遇到超大TIF。我实测一张4波段、单波段10000×10000的TIF，如果用整幅读入再做金字塔，内存峰值可能到3GB以上；如果分块逐窗口融合，内存占用能控制在几百MB。分块实现时会引入块与块之间的边缘伪影，但拉普拉斯金字塔融合的伪影不明显，实测下来可以接受。我建议512×512的窗口最平衡，既不产生过多文件读写开销，也不会撑爆内存。

6. 实操经验与个人心得

我个人的体会是，图像融合TIF这事，代码本身不复杂，真正的难点往往在前期数据检查和后期结果验证。拿到数据先看三样东西：尺寸、波段数、坐标范围。尺寸不一致先重采样，波段顺序别搞混，坐标系不一致就先投影转换。做完融合再看三样东西：道路边缘是否清晰、地物颜色是否偏离实际、地理叠加是否贴合。这三关过了，基本就能交付。

最后分享一个小技巧：无论是Python还是MATLAB版本，我建议都封装成命令行或函数接口，入参是“输入A路径+输入B路径+权重+层数+输出路径”，出参是融合后TIF。这样一旦方案稳定，就可以批量去处理整个目录的影像，而不是一张张手动跑。做遥感的人都知道，批处理才是效率的王道。

本文还有配套的精品资源，点击获取