如何大幅提高遥感影像分辨率（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~=9img_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);elsej=j+1;endend
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 sizeprint(band8_1.shape)  # outputmiu_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_1def 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_p1HF_list.append(HF)band_aver = np.mean(miu_i_list)HRM = band + (band / band_aver) * HFHRM_list.append(HRM)return miu_i_list, std_i_list, HF_list, HRM_listif __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 allHSM_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')