遥感图像融合 |
您所在的位置：网站首页 › 遥感图像融合实验心得体会 › 遥感图像融合 |
遥感图像融合
|
一、简介
遥感图像融合的定义是通过将多光谱低分辨率的图像和高分辨率的全色波段进行融合从而得到信息量更丰富的遥感图像。常用的遥感图像融合方法有Brovey\PCA\Gram-Schmidt方法。其中Gram-Schmidt方法效果较好，且应用广泛。 多源图像融合可类比为多源信息融合，即通过获取源数据原始显著特征，然后通过适当的融合方法将这些优势特征集成到单个图像内。经过多年各个领域的研宄和突破，这些融合方法取得了非凡的融合性能，并广泛应用于许多重要领域。军事领域最先利用多源信息融合研宄成果。早期美国防部采用信息融合技术对声呐信号进行处理，由此信息融合技术开始被广泛关注。现在，世界各国在遥感领域都投入大量人力物力并取得长足发展，较易获得高质量且多元化的数据，因此，多源遥感数据融合，例如红外图像与可见光图像之间融合，高光谱图像与雷达数据融合，成为当前备受关注的研究课题。 随着深度学习兴起，基于深度学习的融合算法被大量提出。通过卷积神经网络提取图像基本表征，特征重构后获得融合结果。在这些基于深度网络的融合策略中，只有最后一层的特征被用作源数据的主要特征，显然，在对图像处理期间，中间层丰富有效的特征被丢失。 根据现状调研，多源遥感图像融合主要目的之一是为了保留光谱信息的同时，提升图像空间分辨率。通过融合空间分辨率较高的图像，来得到高空间分辨率的高／多光谱图像。融合结果可综合多种数据源的优势，达到增强多源数据协同解译能力的目的。 目前图像融合算法有很多种，通常可划分成以下三类： 基于传统方法图像融合，其原理是针对图像的各个波段，使用变换、数值运算等容易实施的方法处理，如线性加权法、高通滤波法、IHS变换法主成分分析法等。 基于多尺度图像融合，其基本理论是把图像解析为不同频率特性、不同分辨率等子信号，然后融合相应分解级别的子信号，将融合后子信号重构得到融合图像，如小波变换融合算法。 基于模型图像融合，该融合方法理论前提是假定空间分辨率较低的图像经过分辨率较高的图像下采样或者其他降低空间分辨率的算法得到。根据上述假设，通过构建能量泛函表征，并且构建高、低分辨率图像之间的模型映射，最后通过求解优化模型得到融合图像。 二、方法 1 gdal_pansharpen.py脚本实现图像融合脚本位置： #根据自己项目路径里找 #D:\python\Anaconda3.7\envs\ttt\Lib\site-packages\osgeo\scriptsgdal_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 # 其实也可以去尝试结合其他波段进行操作效果：
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变换方法不同的是可以将图像分解成多个成分，对于包含轮廓信息的第一主成分进行替换，有效的提高了多光谱图像的空间分辨率，但是同样因为替换成分之间的低相关性造成了颜色的扭曲。
目录 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实现： 分析： 我们先看下公式， 再看一眼整体公式 就是遍历所有的像素点之后求均值。 实现： 卷积去解决 公式: 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 outskimage.measure里面有现成的函数 3 角二阶矩公式： 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公式： 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. 峰值信噪比PSNRPython实现： 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 |
【本文地址】
今日新闻 |
点击排行 |
|
推荐新闻 |
图片新闻 |
|
专题文章 |