如何大幅提高遥感影像分辨率（Python+MATLAB）

前言：

算法：NSCT算法（非下采样变换）
数据：Landsat8 OLI 遥感图像数据
编程平台：MATLAB+Python
论文参考：毛克.一种快速的全色和多光谱图像融合算法[J].测绘科学,2016,41(01):151-153+98.DOI:10.16251/j.cnki.1009-2307.2016.01.028.

             左图：未进行融合的多光谱真彩色合成                              右图：利用NSCT变换后的多光谱真彩色合成

算法概述

• （1）上采样低分辨率多光谱图像LRM，使其和全色图像HRP的分辨率相同，记为LRM杠；
• （2）根据第i个多光谱波段LRMi杠来匹配全色图像得到新的图像HRPi杠；
  
• （3）使用NSCT变换两层分解HRPi杠，得到低频分量和细节子带；
• 此处虽然说的是两层分解，但是我找遍了资料，还是不能理解是怎么得出。但是考虑到是低频分量，于是我利用的高斯模糊，方差为20，提取出了全色波段的低频部分。
  

                                          图：全色波段的低频分量

• （4）提高全色高频信息HF，即
  
• (5)按比例投影高频信息HF至LRMi杠中，得到融合图像HRMi
  

matlab批量读取遥感影像并将影像矩阵存入.mat文件中：

clc;clear all;
file_path="C:\Users\稳魂\Desktop\data\";
img_path_list = dir(strcat(file_path,'LC08_L1TP_118038_20130525_20170504_01_T1_B*.TIF'));
img_num = length(img_path_list);
fprintf('正在读取的图像为：\n');
data_base={};
name={};
if img_num > 0 
    for j = 1:img_num 
        if j~=9
                img_name = strcat(file_path,'LC08_L1TP_118038_20130525_20170504_01_T1_B',int2str(j),'.TIF');
                pitch=geotiffread(img_name);
                name{j}=img_name;
                data_base{j}=pitch;
                fprintf('第%02d个：%s\n',j,img_name);
        else
                j=j+1;
        end
    end
end
save("D:\pythonProject\data_base.mat","data_base");

Python主程序代码（NSCT融合算法）

附加解释：
def load_mat：读取matlab的mat文件，对全色波段影像进行均值、方差、低频分量提取操作；
def hsc_data： 对多光谱影像进行上述均值、方差操作，并进行NSCT算法对各多光谱影像进行融合。

# coding=utf-8
from typing import List
import scipy.io as scio
import numpy as np
import cv2 as cv

# matlab//python混编.mat数据传递
# 读取.mat数据
def load_mat(path):
    data_data1 = scio.loadmat(path)
    band8_1 = data_data1['data_base'][0, 7]
    row1, col1 = band8_1.shape  # get band8 size
    print(band8_1.shape)  # output
    miu_p1 = np.nanmean(band8_1)  # 计算均值
    std_p1 = np.nanstd(band8_1)  # 计算标准差
    # 获取图像的低频分量//高斯模糊
    L_p1 = cv.GaussianBlur(band8_1, (0, 0), 20)
    # cv.imwrite("my.tif", L_p)
    # print(miu_p)
    # print(std_p)
    return data_data1, row1, col1, miu_p1, std_p1, L_p1, band8_1

def hsc_data(data1, row1, col1, miu_p1, std_p1, L_p1, band81):
    # band_list = []
    miu_i_list = []  # 均值
    std_i_list = []  # 标准差
    # HRP_dd_list = []  # 各个波段的HRP影像
    HF_list = []  # 各个波段的高频分量
    HRM_list = []  # 融合图像列表

    for i in range(0, 4):
        band = data1['data_base'][0, i]
        band = cv.resize(band, (col1, row1))
        # band_list.append(band)
        miu_i = np.nanmean(band)
        miu_i_list.append(miu_i)
        std_i = np.nanstd(band)
        std_i_list.append(std_i)
        HRP_dd = ((band81 - miu_p1) / std_p1) * std_i + miu_i
        # HRP_dd_list.append(HRP_dd)
        HF = HRP_dd - L_p1
        HF_list.append(HF)
        band_aver = np.mean(miu_i_list)
        HRM = band + (band / band_aver) * HF
        HRM_list.append(HRM)
    return miu_i_list, std_i_list, HF_list, HRM_list


if __name__ == '__main__':
    path1 = r"D:\pythonProject\data_base.mat"
    data_data, row, col, miu_p, std_p, L_p, band8 = load_mat(path1)
    miu_i_l, std_i_l, HF_l, HRM_l = hsc_data(data_data, row, col, miu_p, std_p,
                                             L_p, band8)
    # print(HRM_l[1])
    # cv.imwrite('HRM_4.tif', HRM_l[3])

如下结果是NSCT变换前与变换的图，利用arcmap工具打开：

                                                    图：未进行NSCT变换的原图

                                                   图：进行NSCT变换后的图

matlab真彩色合成

clc;clear all

HSM_2=uint16(imread("HRM_2.tif"));
HSM_3=uint16(imread("HRM_3.tif"));
HSM_4=uint16(imread("HRM_4.tif"));
hsm_img=cat(3,HSM_4,HSM_3,HSM_2);
imwrite(hsm_img,'HSM.tif','tif');
%% 

band2=imread("LC08_L1TP_118038_20130525_20170504_01_T1_B2.TIF");
band3=imread("LC08_L1TP_118038_20130525_20170504_01_T1_B3.TIF");
band4=imread("LC08_L1TP_118038_20130525_20170504_01_T1_B4.TIF");
band_img=cat(3,band4,band3,band2);
imwrite(band_img,'BAND.tif','tif')