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!')

多光谱影像

全色影像

融合影像

细节对比