目的意义¶
掌握主流图像融合方法，包括代数运算法（如加权平均、相乘、Brovey变换）和成分替换法（如IHS、PCA、Gram-Schmidt融合），理解其数学原理、适用条件及优缺点。
利用Python平台实现多种融合算法，处理真实的高分二号（GF-2）卫星多光谱（MSS）与全色（PAN）影像数据。
学习并应用RMSE、CC、UIQI、ERGAS、SAM等定量指标，科学评估融合影像在空间细节增强与光谱保真度之间的平衡。
实习步骤¶
读取 data/GF2-MSS.tif（多光谱影像）和 data/GF2-PAN.tif（全色影像）。 将多光谱影像重采样至全色影像的空间分辨率（建议使用双线性或双三次插值以减少锯齿效应）。
代数运算法（三选一或多种）： 加权平均法：fused = w₁·PAN + w₂·MS_band 相乘法：fused = PAN × MS_band Brovey变换：fused_gs = (PAN × MS_band) / ΣMS_bands
成分替换法（必做IHS与PCA，GS为加分项）： IHS融合：RGB→IHS变换 → 用矩匹配后的PAN替换I分量 → IHS→RGB逆变换。 PCA融合：对重采样MS影像做主成分分析 → 用矩匹配后的PAN替换第一主成分 → 逆PCA变换。 GS融合（可选）：构造辅助波段 → 施密特正交化 → 替换第一分量 → 逆变换。
融合结果保存 将各方法生成的高分辨率多光谱影像保存至 result/ 目录。
融合精度定量评价 至少计算两种评价指标（如RMSE、CC、UIQI、ERGAS、SAM），对比原始多光谱影像与融合影像，评估光谱保真度与空间增强效果。
In [23]:
import rasterio as rio
import numpy as np
import warnings
warnings.filterwarnings("ignore", category=rio.errors.NotGeoreferencedWarning)
def normalize(array):
if array.ndim == 2:
array_min, array_max = array.min(), array.max()
return (array - array_min) / (array_max - array_min)
elif array.ndim == 3:
for i in range(array.shape[0]):
array_min, array_max = array[i].min(), array[i].max()
array[i] = (array[i] - array_min) / (array_max - array_min)
return array
with rio.open("data/GF2-MSS.tif") as src:
mss = src.read().astype(np.float64)
mss = mss[0:3, :, :] # 只取前三个波段（RGB）
mss = mss[::-1, :, :] # 波段顺序调整为RGB
mss = normalize(mss)
mss = np.repeat(np.repeat(mss, 4, axis=1), 4, axis=2) # 重采样到与全色影像相同分辨率
mss_meta = src.meta
with rio.open("data/GF2-PAN.tif") as src:
pan = src.read(1).astype(np.float64)
pan = normalize(pan)
pan_meta = src.meta
from rasterio.plot import show
show(mss, title="original MSS image")
show(pan, title="original PAN image")
Out[23]:
<Axes: title={'center': 'original PAN image'}>
In [24]:
# 加权平均融合
def weighted_average_fusion(mss, pan, alpha=0.5):
fused = np.empty_like(mss)
for i in range(mss.shape[0]):
fused[i] = alpha * mss[i] + (1 - alpha) * pan
return normalize(fused)
fused_wa = weighted_average_fusion(mss, pan)
show(fused_wa, title="Weighted Average Fusion Result")
Out[24]:
<Axes: title={'center': 'Weighted Average Fusion Result'}>
In [25]:
# 相乘法融合
def multiplicative_fusion(mss, pan):
fused = np.empty_like(mss)
for i in range(mss.shape[0]):
fused[i] = mss[i] * pan
return normalize(fused)
fused_mult = multiplicative_fusion(mss, pan)
show(fused_mult, title="Multiplicative Fusion Result")
Out[25]:
<Axes: title={'center': 'Multiplicative Fusion Result'}>
In [26]:
# Brovey变换融合
def brovey_fusion(mss, pan):
fused = np.empty_like(mss)
sum_mss = np.sum(mss, axis=0)
for i in range(mss.shape[0]):
fused[i] = (mss[i] / sum_mss) * pan
return normalize(fused)
fused_brov = brovey_fusion(mss, pan)
show(fused_brov, title="Brovey Fusion Result")
Out[26]:
<Axes: title={'center': 'Brovey Fusion Result'}>
In [27]:
# IHS融合
def rgb_to_ihs(rgb):
r, g, b = rgb[0], rgb[1], rgb[2]
intensity = (r + g + b) / 3
r1 = (2*b - r - g) / (3 * np.sqrt(2))
r2 = (r - g) / np.sqrt(2)
hue = np.arctan2(r2, r1)
saturation = np.sqrt(r1**2 + r2**2)
return np.array([intensity, saturation, hue])
def ihs_to_rgb(ihs):
intensity, saturation, hue = ihs[0], ihs[1], ihs[2]
r1 = saturation * np.cos(hue)
r2 = saturation * np.sin(hue)
r = intensity + (r2 - r1) / np.sqrt(2)
g = intensity - (r1 + r2) / np.sqrt(2)
b = intensity + r1 * np.sqrt(2)
return np.array([r, g, b])
def match_histogram(source, target):
source_mean, source_std = np.mean(source), np.std(source)
target_mean, target_std = np.mean(target), np.std(target)
matched = (source - source_mean) * (target_std / source_std) + target_mean
return matched
def ihs_fusion(mss, pan):
fused = np.empty_like(mss)
ihs = rgb_to_ihs(mss)
ihs[0] = match_histogram(pan, ihs[0])
fused = ihs_to_rgb(ihs)
return normalize(fused)
fused_ihs = ihs_fusion(mss, pan)
show(fused_ihs, title="IHS Fusion Result")
Out[27]:
<Axes: title={'center': 'IHS Fusion Result'}>
In [28]:
# PCA融合
def pca_fusion(mss, pan):
reshaped_mss = mss.reshape(mss.shape[0], -1)
covariance_matrix = np.cov(reshaped_mss)
eigenvalues, eigenvectors = np.linalg.eig(covariance_matrix)
principal_component = eigenvectors[:, np.argmax(eigenvalues)]
pca_image = np.dot(principal_component.T, reshaped_mss).reshape(mss.shape[1], mss.shape[2])
pca_image = match_histogram(pan, pca_image)
fused = np.empty_like(mss)
for i in range(mss.shape[0]):
fused[i] = mss[i] + (pca_image - pca_image.mean())
return normalize(fused)
fused_pca = pca_fusion(mss, pan)
show(fused_pca, title="PCA Fusion Result")
Out[28]:
<Axes: title={'center': 'PCA Fusion Result'}>
In [29]:
# Gram-Schmidt融合
def gs_fusion(mss, pan):
fused = np.empty_like(mss)
g1 = (mss[0] + mss[1] + mss[2]) / 3
g2 = mss[0] - g1
g3 = mss[1] - g1
g4 = mss[2] - g1
g1 = match_histogram(pan, g1)
fused[0] = g1 + g2
fused[1] = g1 + g3
fused[2] = g1 + g4
return normalize(fused)
fused_gs = gs_fusion(mss, pan)
show(fused_gs, title="GS Fusion Result")
Out[29]:
<Axes: title={'center': 'GS Fusion Result'}>
In [30]:
# 结果评价
def resample_to_mss(fused, target_shape):
from skimage.transform import resize
n_bands = fused.shape[0]
resampled = np.empty((n_bands, target_shape[1], target_shape[2]))
for i in range(n_bands):
resampled[i] = resize(fused[i], (target_shape[1], target_shape[2]), mode='reflect', anti_aliasing=True)
return resampled
def rmse(original, fused):
return np.sqrt(np.mean((original - fused) ** 2))
def cc(original, fused):
from scipy.stats import pearsonr
return pearsonr(original.flatten(), fused.flatten())[0]
def uiqi(original, fused):
mean_orig = np.mean(original)
mean_fused = np.mean(fused)
var_orig = np.var(original)
var_fused = np.var(fused)
cov = np.mean((original - mean_orig) * (fused - mean_fused))
return (4 * cov * mean_orig * mean_fused) / ((var_orig + var_fused) * (mean_orig**2 + mean_fused**2) + 1e-6)
def ergas(original, fused):
ratio = 4 # 分辨率比
n_bands = original.shape[0]
sum_term = 0
for i in range(n_bands):
rmse_band = rmse(original[i], fused[i])
mean_orig_band = np.mean(original[i])
sum_term += (rmse_band / (mean_orig_band + 1e-6)) ** 2
ergas_value = 100 / ratio * np.sqrt(sum_term / n_bands)
return ergas_value
def sam(original, fused):
dot_product = np.sum(original * fused, axis=0)
norm_orig = np.linalg.norm(original, axis=0)
norm_fused = np.linalg.norm(fused, axis=0)
cos_theta = dot_product / (norm_orig * norm_fused + 1e-6)
cos_theta = np.clip(cos_theta, -1, 1)
angles = np.arccos(cos_theta)
return np.mean(angles)
methods = {
"Weighted Average": fused_wa,
"Multiplicative": fused_mult,
"Brovey": fused_brov,
"IHS": fused_ihs,
"PCA": fused_pca,
"GS": fused_gs
}
import pandas as pd
pd.set_option('display.float_format', '{:.4f}'.format)
result = []
for name, fused in methods.items():
fused = resample_to_mss(fused, mss.shape)
for i in range(mss.shape[0]):
result.append([name,
i + 1,
rmse(mss[i], fused[i]),
cc(mss[i], fused[i]),
uiqi(mss[i], fused[i]),
ergas(mss, fused),
sam(mss, fused)])
result_df = pd.DataFrame(result, columns=["Method", "Band", "RMSE", "CC", "UIQI", "ERGAS", "SAM"])
print(result_df)
Method Band RMSE CC UIQI ERGAS SAM 0 Weighted Average 1 0.0508 0.9671 0.9553 6.0129 0.0918 1 Weighted Average 2 0.0540 0.9755 0.9577 6.0129 0.0918 2 Weighted Average 3 0.1127 0.9790 0.9193 6.0129 0.0918 3 Multiplicative 1 0.1392 0.9661 0.7281 14.6677 0.0077 4 Multiplicative 2 0.1474 0.9693 0.7331 14.6677 0.0077 5 Multiplicative 3 0.2370 0.9548 0.6826 14.6677 0.0077 6 Brovey 1 0.1429 0.6206 0.5091 14.5465 0.1120 7 Brovey 2 0.1566 0.7080 0.5580 14.5465 0.1120 8 Brovey 3 0.2123 0.8984 0.6105 14.5465 0.1120 9 IHS 1 0.2102 0.7417 0.5437 17.1851 0.1642 10 IHS 2 0.1646 0.7782 0.6170 17.1851 0.1642 11 IHS 3 0.2014 0.8975 0.6442 17.1851 0.1642 12 PCA 1 0.0925 0.8978 0.8618 10.5099 0.1605 13 PCA 2 0.1057 0.9177 0.8513 10.5099 0.1605 14 PCA 3 0.1794 0.9046 0.7490 10.5099 0.1605 15 GS 1 0.2102 0.7417 0.5437 17.1851 0.1642 16 GS 2 0.1646 0.7782 0.6170 17.1851 0.1642 17 GS 3 0.2014 0.8975 0.6442 17.1851 0.1642
In [31]:
# 结果总结
summary = result_df.groupby("Method").mean().reset_index()
summary.drop(columns=["Band"], inplace=True)
summary.sort_values(by="RMSE", inplace=True)
print(summary)
Method RMSE CC UIQI ERGAS SAM 5 Weighted Average 0.0725 0.9739 0.9441 6.0129 0.0918 4 PCA 0.1259 0.9067 0.8207 10.5099 0.1605 0 Brovey 0.1706 0.7423 0.5592 14.5465 0.1120 3 Multiplicative 0.1745 0.9634 0.7146 14.6677 0.0077 2 IHS 0.1921 0.8058 0.6016 17.1851 0.1642 1 GS 0.1921 0.8058 0.6016 17.1851 0.1642
结果与结论¶
加权平均法在所有指标中表现最优，尤其在ERGAS（6.01）和RMSE（最低）方面显著优于其他方法，说明其在保持原始光谱信息的同时较好地引入了空间细节。
相乘法虽然SAM值极低（表明光谱角偏差小），但UIQI和ERGAS指标较差，说明其虽保留了光谱方向一致性，却牺牲了整体结构相似性和亮度保真度。
Brovey变换与IHS/GS融合在光谱保真度方面表现较弱（CC普遍低于0.8，ERGAS >14），尤其IHS和GS因用PAN完全替代强度/主成分，导致光谱失真明显。
PCA融合介于代数法与成分替换法之间，在光谱保持（CC≈0.91）与空间增强之间取得较好平衡。
代数运算法（如加权平均）更适合对光谱保真度要求较高的应用场景，因其操作简单且能有效控制光谱畸变。 成分替换法（如IHS、PCA、GS）虽能显著提升空间分辨率，但往往以牺牲光谱信息为代价，适用于更关注纹理和边缘细节的任务（如城市建筑识别）。 定量评价指标需综合使用：单一指标易产生误导（如相乘法SAM低但UIQI差），应结合RMSE（误差大小）、CC（线性相关性）、UIQI（结构相似性）、ERGAS（整体质量）和SAM（光谱角度）多维度判断。