　　多源图像融合可类比为多源信息融合，即通过获取源数据原始显著特征，然后通过适当的融合方法将这些优势特征集成到单个图像内。经过多年各个领域的研宄和突破，这些融合方法取得了非凡的融合性能，并广泛应用于许多重要领域。军事领域最先利用多源信息融合研宄成果。早期美国防部采用信息融合技术对声呐信号进行处理，由此信息融合技术开始被广泛关注。现在，世界各国在遥感领域都投入大量人力物力并取得长足发展，较易获得高质量且多元化的数据，因此，多源遥感数据融合，例如红外图像与可见光图像之间融合，高光谱图像与雷达数据融合，成为当前备受关注的研究课题。

　　随着深度学习兴起，基于深度学习的融合算法被大量提出。通过卷积神经网络提取图像基本表征，特征重构后获得融合结果。在这些基于深度网络的融合策略中，只有最后一层的特征被用作源数据的主要特征，显然，在对图像处理期间，中间层丰富有效的特征被丢失。根据现状调研，多源遥感图像融合主要目的之一是为了保留光谱信息的同时，提升图像空间分辨率。通过融合空间分辨率较高的图像，来得到高空间分辨率的高／多光谱图像。融合结果可综合多种数据源的优势，达到增强多源数据协同解译能力的目的。

　　目前图像融合算法有很多种，通常可划分成以下三类：

基于传统方法图像融合，其原理是针对图像的各个波段，使用变换、数值运算等容易实施的方法处理，如线性加权法、高通滤波法、IHS变换法主成分分析法等。
基于多尺度图像融合，其基本理论是把图像解析为不同频率特性、不同分辨率等子信号，然后融合相应分解级别的子信号，将融合后子信号重构得到融合图像，如小波变换融合算法。
基于模型图像融合，该融合方法理论前提是假定空间分辨率较低的图像经过分辨率较高的图像下采样或者其他降低空间分辨率的算法得到。根据上述假设，通过构建能量泛函表征，并且构建高、低分辨率图像之间的模型映射，最后通过求解优化模型得到融合图像。

二、方法

1 gdal_pansharpen.py脚本实现图像融合

脚本位置：

#根据自己项目路径里找
#D:\python\Anaconda3.7\envs\ttt\Lib\site-packages\osgeo\scripts

gdal_pansharpen.py的主要参数

Usage: gdal_pansharpen [--help-general] pan_dataset {spectral_dataset[,band=num]}+ out_dataset
                       [-of format] [-b band]* [-w weight]*
                       [-r {nearest,bilinear,cubic,cubicspline,lanczos,average}]
                       [-threads {ALL_CPUS|number}] [-bitdepth val] [-nodata val]
                       [-spat_adjust {union,intersection,none,nonewithoutwarning}]
                       [-verbose_vrt] [-co NAME=VALUE]* [-q]
主要参数解释如下：
pan_dataset：全色波段数据集，作为输入的第一项数据
spectral_dataset：多光谱波段数据，具有一个或多个光谱波段的数据集。
out_dataset：融合后要输出的数据
-of out_format：输出文件格式，默认GeoTiff格式
-b band:从输入光谱波段中选择波段进行输出。从1开始编号，编号顺序为指定的光谱波段，当这个参数不做特别设置时，即默认所有输入光谱波段进行输出。
-w weight：指定计算伪全色值的权重。有多少个输入波段，就设置几个w
-r：重采样算法，默认为cubic
-n nodata_value：为波段指定nodata值

准备好的待融合的数据，实验一下，下面两张图是配准好的待融合的影像

　　利用上面介绍的pansharpen.py来实现上述影像数据的融合操作，具体实现如下。下图为融合后的效果，看上去还可以

os.chdir(r'D:\Anaconda3\envs\RS\Scripts')
%run gdal_pansharpen.py -of GTiff D:\gdal_data\pan.tif D:\gdal_data\multi.tif D:\gdal_data\\pansharpen.tif

　　再试一下设置不同权重融合后是什么效果，结果就是各波段的值会随着权重的不同而发生变化。不设置W时，默认的每个波段的权重为0.25（验证了）。

　　同样，也为了实践一下gdal_pansharpen.py的效率，选择一整景的高分1号数据测试一下融合的耗时，具体实现如下：

a = time.time()
%run gdal_pansharpen.py -nodata 0 -of GTiff D:\gdal_data\GF1_PAN_ortho.tiff D:\gdal_data\GF1_MULTI_ortho.tiff D:\gdal_data\GF1pansharpen.tif
b = time.time()  
print(b-a)

　　耗时84秒左右，也就是1分半不到融合一景数据，效率还是可以的。

2 HSV融合

　　HSV融合的步骤，首先是将多光谱多光谱图像重采样至全色图像的尺寸，其次将RGB图像变换HSV颜色空间，接着用高分辨率的图像代替转换后的颜色亮度值波段（V），最后再将转换后的图像变换回RGB颜色空间（以高分1数据为例)：

2.1 数据读取，包括多光谱、全色波段数据读取及相关的数据处理。

# 读取多光谱数据
ds = gdal.Open('./multi.tif')
multi_arr = ds.ReadAsArray()
multi_arr.shape
# 读取全色数据
ds = gdal.Open('./pan_1.tif')
pan_arr = ds.ReadAsArray()
pan_arr.shape
# 定义一个包含4各波段的arr，大小与multi_arr一致，需要注意的是shape的顺序，并依次将blue,green,red,nir波段赋值
rgbn = np.zeros((270,406,4))
for i in range(multi_arr.shape[0]):
    rgbn[:,:,i] = multi_arr[i,:,:]

2.2 实现多光谱多光谱图像重采样至全色图像的尺寸。

# 定义一个待重采样的数组，shape大小等于multi_arr.shape乘以4？？ 原因是高分1数据的多光谱分辨率是8米，全色分辨率是2米，二者之间相差4倍
rgbn_interpolation = np.zeros((multi_arr.shape[1] * 4, multi_arr.shape[2] * 4, 4))
rgbn_interpolation.shape

#利用skimage库里的resacle进行重采样

from skimage.transform import rescale

for i in range(4):
    img = rgbn[:,:,i]
    rgbn_interpolation[:,:,i] = rescale(img, (4,4))
　　#rgbn_interpolation[:,:,i] = cv2.resize(img,multi_arr.shape[2] * 4, multi_arr.shape[1] * 4)

2.3 将重采样后的图像的RGB图像变换HSV颜色空间，利用全色波段替换转换后图像中的V波段

#利用skimage库里的color实现图像空间转换，也可以用opencv来实现同样的操作
hsv = color.rgb2hsv(rgbn_interpolation[:,:,:3])

pan_arr1= cv2.resize(pan_arr,multi_arr.shape[2] * 4, multi_arr.shape[1] * 4)

hsv[:,:,2] = pan_arr1

2.4 将转换后的图像重新变换回RGB颜色空间

hsv_p_img = color.hsv2rgb(hsv)

　　至此就实现了HSV的图像融合，其实整个过程里比较关键的步骤就是多光谱图像重采样至全色图像的大小。重采样的方法多种，比如最近邻、双线性或三次卷积等等，都可以去尝试。文中默认用的是skimage中scale中重采样方式。

　　其实在实现图像重采样的过程中，会出现一些问题，也是容易出错的问题，那就是重采样的图像按照全色、多光谱的空间分辨率的倍数进行操作，操作后极有可能出现重采样图像的shape与全色图像的shape大小不一致，如何解决，大家可以思考一下。提示一下，可以用替代的方式，就是保证shape范围小的，可以解决这个问题。（文中示例的高分样例数据，是我重采样过的，保证了全色和多光谱数据重采样之后的shape大小是一致的，所以不会出现上面说的问题，但是大家在操作的时候可能会出现上述问题。）

　　在利用全色波段替换V波段时，可以对全色波段进行权重的操作（设置W权重）

# 结合全色波段和nir波段进行V波段转换的操作，对输出的结果会有一定影响，这是必然的，因为改变了V波段的赋值
hsv[:,:,2] = pan_arr - w * band_nir
# 其实也可以去尝试结合其他波段进行操作

效果：

3 Gram-Schmidt图像融合

　　Gram-Schmidt方法效果较好，且应用广泛。该方法由CraigA.Laben等人提出，已经被封装到多个遥感图像处理软件中。对于此算法的叙述，国内的李存军写的《两种高保真遥感影像融合方法比较》复述的很清楚，结合原文看清晰易懂。

1.首先预处理数据，计算多光谱影像和全色波段重叠区域，得到裁剪后的多光谱影像和全色波段。

2.随后模拟产生低分辨率的全色波段影像用于作为GS变换的第一分量。通常是将低分辨率的多光谱影像根据光谱响应函数按一定权重wi进行模拟，得到模拟的全色波段灰度值。或者把全色波段影像模糊，缩小到与多光谱影像相同大小。

这里我们最终对多光谱影像，按波段计算了平均值，来模拟全色波段。

3.接着就是重头环节。GS变换--施密特正交化，具体原理可以百度，这里给出修改后的施密特正交化公式。其中h()是计算矩阵内积，然后做除法。以模拟波段为第一波段，多光谱影像所有波段为后续波段，做GS变换。

　　　　　　　　　　　　施密特正交化

　　　　　　　　　　GS融合正变换

4.接着根据GS第一分量，即模拟波段的mean和var，对全色波段进行修改。

5.然后把修改后的全色波段作为第一分量，进行GS逆变换，输出n+1个波段，去除第一个波段，就是融合后的结果。

最后分析一下具体编码步骤：

1）overlay，求重叠区域图像的函数

2）resample，重采样把多光谱影像重采样到全色波段的形式

3）simulate，模拟全色波段的函数

4）GS正变换

5）modify函数，修改全色波段作为GS第一分量

6）GS逆变换

4 基于IHS变换的图像融合方法

　　IHS方法是将原始多光谱图像从RGB空间变换到IHS空间，然后用高分辨率图像或用不同投影方式得到的待融合图像替代I分量。在IHS系统中，三种成分的相关性比较低，这使得我们能够对这三种分量分别进行处理。I成分主要反映地物辐射总的能量以及空间分布，即表现为空间特征。而H,S则反映光谱信息。

HIS为：亮度（I ）、色调(H)、饱和度(S)；

强度表示光谱的整体亮度大小，对应于图像的空间分辨率；

传统的IHS图像融合方法基本思想是将IHS空间中的低分辨率亮度用高分辨率的图像的亮度成分所代替；

首先是正变换（RGB->IHS)

逆变换（IHS->RGB)

实验代码

def IHS(data_low,data_high):
    """
    基于IHS变换融合算法
    输入：np.ndArray格式的三维数组
    返回：可绘出图像的utf-8格式的三维数组
    """
    A = [[1./3.,1./3.,1./3.],[-np.sqrt(2)/6.,-np.sqrt(2)/6.,2*np.sqrt(2)/6],[1./np.sqrt(2),-1./np.sqrt(2),0.]] 
    #RGB－>IHS正变换矩阵
    B = [[1.,-1./np.sqrt(2),1./np.sqrt(2)],[1.,-1./np.sqrt(2),-1./np.sqrt(2)],[1.,np.sqrt(2),0.]] 
    #IHS－>RGB逆变换矩阵
    A = np.matrix(A)
    B = np.matrix(B)
    band , w , h = data_high.shape
    pixels = w * h
    data_low = data_low.reshape(3,pixels)
    data_high = data_high.reshape(3,pixels)
    a1 = np.dot(A , np.matrix(data_high))#高分影像正变换
    a2 = np.dot(A , np.matrix(data_low))#低分影像正变换
    a2[0,:] = a1[0,:]#用高分影像第一波段替换低分影像第一波段
    RGB = np.array(np.dot(B , a2))#融合影像逆变换
    RGB = RGB.reshape((3,w,h))
    min_val = np.min(RGB.ravel())
    max_val = np.max(RGB.ravel())
    RGB = np.uint8((RGB.astype(np.float) - min_val) / (max_val - min_val) * 255)
    RGB = Image.fromarray(cv2.merge([RGB[0],RGB[1],RGB[2]]))
return RGB
 
def imresize(data_low,data_high):
    """
    图像缩放函数
    输入：np.ndArray格式的三维数组
    返回：np.ndArray格式的三维数组
    """
    band , col , row = data_high.shape
    data = np.zeros(((band,col,row)))
    for i in range(0,band):
            data[i] = smi.imresize(data_low[i],(col,row))
return data
 
def gdal_open(path):
    """
    读取图像函数
    输入：图像路径
    返回：np.ndArray格式的三维数组
    """
    data = gdal.Open(path)
    col = data.RasterXSize#读取图像长度
    row = data.RasterYSize#读取图像宽度
    data_array_r = data.GetRasterBand(1).ReadAsArray(0,0,col,row).astype(np.float)#读取图像第一波段并转换为数组
    data_array_g = data.GetRasterBand(2).ReadAsArray(0,0,col,row).astype(np.float)#读取图像第二波段并转换为数组
    data_array_b = data.GetRasterBand(3).ReadAsArray(0,0,col,row).astype(np.float)#读取图像第三波段并转换为数组
    data_array = np.array((data_array_r,data_array_g,data_array_b))
return data_array

def main(path_low,path_high):
    data_low = gdal_open(path_low)
    data_high = gdal_open(path_high)
    data_low = imresize(data_low,data_high)
    RGB = IHS(data_low,data_high)
    RGB.save("IHS.png",'png')
if __name__ == "__main__":
    path_low = 'RGB.tif'
    path_high = 'Band8.tif'
    main(path_low,path_high)

5 基于PCA变换的图像融合方法

　　对于PCA变换的融合方法与IHS变换方法不同的是可以将图像分解成多个成分，对于包含轮廓信息的第一主成分进行替换，有效的提高了多光谱图像的空间分辨率，但是同样因为替换成分之间的低相关性造成了颜色的扭曲。

6 IHS结合小波变换的双输入单输出流程图

7 IHS+NSCT

8 PCA+NSCT

三、利用Python实现遥感图像融合评价指标

1. 点锐度EVA(Python实现) ;
2. 信息熵Entropy(Python实现);
3. 角二阶矩ASM(Python实现)；
4. 光谱角测度SAM(Python实现)；
5. 结构相似度SSIM(Python实现);
6. 峰值信噪比PSNR(Python实现)；
7. 其他的指标Python实现。

无参考的有EVA、ASM、Entropy；有参考的SAM、SSIM、PSNR。

1 点锐度

公式：

　　其中，m、n为图像的长和宽，df 为灰度变化幅值，dx 为像元间的距离增量。

Python实现：

分析：

我们先看下公式， $\left| df/dx \right|$ 是什么？

　　 $\left| df/dx \right|$ 的实现直接用相邻或者差一个位置的像素相减就可以实现。咱多看点 $\sum_{a=1}^{8}\left| df/dx \right|$ ，这个式子的意思就是每个像素与周围八个像素的梯度之和。

逐个对图像中的每点取 8 邻域点与之相减.先求 8 个差值的加权和，权的大小取决于距离.距离远.则权小.如 45° 和 135° 方向的差值需乘以需要乘以 $1/\sqrt{2}$ 。再将所有点所得值相加。

　　再看一眼整体公式 $P = \sum_{i=1}^{m \times n} \sum_{a=1}^{8}\left| df/dx \right| /\ m \times n \\$

　　就是遍历所有的像素点之后求均值。

实现：

卷积去解决 $\sum_{a=1}^{8}\left| df/dx \right|$ ，然后我们再求矩阵的均值

xiejiao = 1 / 2 ** 0.5
sizhou = 1.0

kernel = np.array([[xiejiao, sizhou, xiejiao],
                   [sizhou,      -8, sizhou],
                   [xiejiao, sizhou, xiejiao]])
dst_matrix = cv2.filter2D(image, -1, kernel)
EVA_index = np.mean(dst_matrix)

2. 信息熵Entropy

公式:

$H(U) = E\left[ -\log(1/p_{i},2) \right]= -\sum_{i=1}^{n}{p_{i}\log (p_{i}, 2)} \\$

Python实现：

 def entropy(img):
     out = 0
     count = np.shape(img)[0]*np.shape(img)[1]
     p = np.bincount(np.array(img).flatten())
     for i in range(0, len(p)):
         if p[i]!=0:
             out-=p[i]*math.log(p[i]/count,2)/count
     return out

　　skimage.measure里面有现成的函数

3 角二阶矩

公式：

$ASM = \sum_{i=1}^{k} \sum_{j=1}^{k} {P(i, j)^{2}} \\$

Python实现：

计算灰度共生矩阵；

# 计算距离为1，2和角度为0度，90度的GLCM的ASM
def asm(img, distances, angles):
    
    # 计算GLCM
    g = skimage.feature.greycomatrix(src_matrix, [1,2],[0,np.pi/2])

    # 计算该glcm的ASM
    asm = skimage.feature.greycoprops(g, 'ASM')
    
    return asm

4 光谱角测度SAM

公式：

$D_{SAM}(X) = cos^{-1}(d^{T}x / (\sqrt{d^{T}d} \sqrt{x^{T}x} )) \\$

Python实现：

def sam(src, dst):
    # src, dst均为numpy数组
    val = np.dot(src.T, dst)/(np.linalg.norm(src)*np.linalg.norm(dst))
    sam = 1.0/np.cos(val)
    return sam

5. 结构相似度SSIM

　　详见：https://www.cnblogs.com/vincent2012/archive/2012/10/13/2723152.html

Python实现：

def ssim(src, dst):
    
    # 数据准备
    # 均值mean
    mean_src = np.mean(src)
    mean_dst = np.mean(dst)
    # 方差var
    var_src = np.var(src)
    var_dst = np.var(dst)
    cov = np.cov(src, dst)
    # 标准差std
    std_src = np.std(src)
    std_dst = np.std(dst)
    # 常数c1,c2,c3
    K1 = 0.01
    K2 = 0.03
    L = 255
    c1 = (K1*L)**2
    c2 = (K2*L)**2
    c3 = c2 / 2


    # 计算ssim
    l = (2*mean_src*mean_dst + c1)/(mean_src**2 + mean_dst**2 + c1)
    c = (2*var_src*var_dst + c2)/(var_src**2 + var_dst**2 + c2)
    s = (cov + c3)/(var_src*var_dst + c3)
    
    ssim = l * c * s
    
    return ssim

直接调用skimage库来计算ssim

ssim = skimage.measure.compare_ssim(src, dst, data_range=255)

6. 峰值信噪比PSNR

Python实现：

def mse(src, dst):
    return np.mean((src.astype(np.float64)-dst.astype(np.float64))**2)

def psnr(src, dst, Max= None):
    if MAX is None:
        MAX = np.iinfo(GT.dtype).max

    mse_value = mse(src, dst)
    
    if mse_value == 0.:
    return np.inf

    return 10 * np.log10(MAX**2 /mse_value)

直接调用skimage

psnr = skimage.measure.compare_psnr(src, dst, data_range=255)

7. 其他的指标：

skimage.measure中除了有SSIM、PSNR、Entropy、MSE、NRMS
平均梯度AG、空间频率SF、标准差STD、互信息MI、标准化NMI，详见https://blog.csdn.net/weixin_37143678/article/details/103537374

参考链接：

https://zhuanlan.zhihu.com/p/140954243

https://zhuanlan.zhihu.com/p/144818150

https://www.cnblogs.com/ljwgis/p/12774005.html

https://blog.csdn.net/qq_39630875/article/details/109333550

https://zhuanlan.zhihu.com/p/136013000

posted on 2021-09-09 21:09 行走的蓑衣客阅读(7693) 评论(0) 编辑收藏举报

刷新页面返回顶部

导航

一、简介

二、方法

1 gdal_pansharpen.py脚本实现图像融合

2 HSV融合

2.1 数据读取，包括多光谱、全色波段数据读取及相关的数据处理。

2.2 实现多光谱多光谱图像重采样至全色图像的尺寸。

2.3 将重采样后的图像的RGB图像变换HSV颜色空间，利用全色波段替换转换后图像中的V波段

2.4 将转换后的图像重新变换回RGB颜色空间

3 Gram-Schmidt图像融合

4 基于IHS变换的图像融合方法

5 基于PCA变换的图像融合方法

6 IHS结合小波变换的双输入单输出流程图

7 IHS+NSCT

8 PCA+NSCT

三、利用Python实现遥感图像融合评价指标

1 点锐度

2. 信息熵Entropy

3 角二阶矩

4 光谱角测度SAM

5. 结构相似度SSIM

6. 峰值信噪比PSNR

7. 其他的指标：

参考链接：