融合目的

多光谱影像具有高光谱分辨率和低空间分辨率，全色影像具有高空间分辨率和低光谱分辨率，将这两种影像进行融合就可以得到同时具有高空间分辨率和高光谱分辨率的影像

PCA主成分变换影像融合

PCA主成分分析原理

参考博客：主成分分析（PCA）原理详解

融合步骤

代码

from osgeo import gdal, gdalconst, gdal_array
import numpy as np
from score import *

if __name__ == \'__main__\':
    pantif = \'./sub_image/SV1-03_20180407_L2A0001046618_1012102060010001_01-PAN.tif\' # 全色影像
    muxtif = \'./sub_image/SV1-03_20180407_L2A0001046618_1012102060010001_01-MUX.tif\' # 多光谱影像
    # 读取全色波段影像
    gdal.AllRegister()
    pands = gdal.Open(pantif, gdalconst.GA_ReadOnly)
    # 宽高
    panX = pands.RasterXSize
    panY = pands.RasterYSize
    # 坐标系
    pansrs = gdal.Dataset.GetSpatialRef(pands)
    # 读取多光谱影像
    muxds = gdal.Open(muxtif, gdalconst.GA_ReadOnly)
    muxsrs = gdal.Dataset.GetSpatialRef(muxds)
    # 数据类型
    muxBand = gdal.Dataset.GetRasterBand(muxds, 1)
    muxType = muxBand.DataType
    # 重采样
    options = gdal.WarpOptions(width=panX, height=panY, srcSRS=muxsrs,
                               dstSRS=pansrs, outputType=gdalconst.GDT_Int32, resampleAlg=gdalconst.GRIORA_Bilinear, creationOptions=[\'COMPRESS=LZW\'])
    gdal.Warp(\"./pca.tif\", muxds, options=options)
    # 主成分分析影像融合
    rsData = gdal_array.LoadFile(\"./pca.tif\")
    nBands, xSzie, ySize = rsData.shape
    rsArray = None
    for i in range(nBands):
        bArray = rsData[i]
        bArray = np.reshape(bArray, (-1, 1))
        if i == 0:
            rsArray = bArray
        else:
            rsArray = np.hstack((rsArray, bArray))
    rsArrayAvg = np.mean(rsArray, axis=0)
    rsArray = rsArray - rsArrayAvg
    u, s, v = np.linalg.svd(rsArray, full_matrices=False)
    rsPCA = rsArray@v.T
    rsPCA1 = rsPCA[:, 0]
    rsPCA1Avg = np.mean(rsPCA1)
    rsPCA1Std = np.std(rsPCA1)
    # 计算全色波段的均值和标准差
    panBand = gdal.Dataset.GetRasterBand(pands, 1)
    panArray = gdal.Band.ReadAsArray(panBand)
    panArray = np.reshape(panArray, (-1, 1))
    panStd = np.std(panArray)
    panAvg = np.mean(panArray)
    # 计算变换系数k1,k2
    k1 = rsPCA1Std/panStd
    k2 = rsPCA1Avg-k1*panAvg
    # 对全色波段进行变换使其与第一主成分的标准差和均值一直
    panArray = panArray*k1+k2
    # 替换第一主成分
    rsPCA[:, 0] = panArray[:, 0]
    # 主成分逆变换
    rsArray = rsPCA@v+rsArrayAvg
    rsds = gdal.Open(\"./pca.tif\", gdalconst.GA_Update)
    for i in range(nBands):
        bArray = rsArray[:, i]
        bArray = np.reshape(bArray, (ySize, xSzie))
        iBand = gdal.Dataset.GetRasterBand(rsds, i+1)
        gdal.Band.WriteArray(iBand, bArray)
    del rsds, pands, muxds
    print(\'ok!\')

融合效果

多光谱影像

全色影像

融合影像

细节对比