遙感圖像融合

一、簡介

　　遙感圖像融合的定義是通過將多光譜低分辨率的圖像和高分辨率的全色波段進行融合從而得到信息量更豐富的遙感圖像。常用的遙感圖像融合方法有Brovey\PCA\Gram-Schmidt方法。其中Gram-Schmidt方法效果較好，且應用廣泛。

　　多源圖像融合可類比為多源信息融合，即通過獲取源數據原始顯著特征，然后通過適當的融合方法將這些優勢特征集成到單個圖像內。經過多年各個領域的研宄和突破，這些融合方法取得了非凡的融合性能，並廣泛應用於許多重要領域。軍事領域最先利用多源信息融合研宄成果。早期美國防部采用信息融合技術對聲吶信號進行處理，由此信息融合技術開始被廣泛關注。現在，世界各國在遙感領域都投入大量人力物力並取得長足發展，較易獲得高質量且多元化的數據，因此，多源遙感數據融合，例如紅外圖像與可見光圖像之間融合，高光譜圖像與雷達數據融合，成為當前備受關注的研究課題。

　　隨着深度學習興起，基於深度學習的融合算法被大量提出。通過卷積神經網絡提取圖像基本表征，特征重構后獲得融合結果。在這些基於深度網絡的融合策略中，只有最后一層的特征被用作源數據的主要特征，顯然，在對圖像處理期間，中間層豐富有效的特征被丟失。 根據現狀調研，多源遙感圖像融合主要目的之一是為了保留光譜信息的同時，提升圖像空間分辨率。通過融合空間分辨率較高的圖像，來得到高空間分辨率的高／多光譜圖像。融合結果可綜合多種數據源的優勢，達到增強多源數據協同解譯能力的目的。

　　目前圖像融合算法有很多種，通常可划分成以下三類：

1. 基於傳統方法圖像融合，其原理是針對圖像的各個波段，使用變換、數值運算等容易實施的方法處理，如線性加權法、高通濾波法、IHS變換法主成分分析法等。
2. 基於多尺度圖像融合，其基本理論是把圖像解析為不同頻率特性、不同分辨率等子信號，然后融合相應分解級別的子信號，將融合后子信號重構得到融合圖像，如小波變換融合算法。
3. 基於模型圖像融合，該融合方法理論前提是假定空間分辨率較低的圖像經過分辨率較高的圖像下采樣或者其他降低空間分辨率的算法得到。根據上述假設，通過構建能量泛函表征，並且構建高、低分辨率圖像之間的模型映射，最后通過求解優化模型得到融合圖像。

二、方法

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進行重采樣
for i in range(4):
    img = rgbn[:,:,i]
    rgbn_interpolation[:,:,i] = rescale(img, (4,4))

2.3 將重采樣后的圖像的RGB圖像變換HSV顏色空間，利用全色波段替換轉換后圖像中的V波段

#利用skimage庫里的color實現圖像空間轉換，也可以用opencv來實現同樣的操作
hsv = color.rgb2hsv(rgbn_interpolation[:,:,:3]) 
hsv[:,:,2] = pan_arr 

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實現：

分析：

我們先看下公式，  是什么？

　　 的實現直接用相鄰或者差一個位置的像素相減就可以實現。咱多看點 ，這個式子的意思就是每個像素與周圍八個像素的梯度之和。

逐個對圖像 中的每點取 8 鄰域點與之相減.先求 8 個差值的加權 和，權的大小取決於距離.距離遠.則權小.如 45° 和 135° 方向的差值需乘以需要乘以   。再將所有點所得值相加 。

　　再看一眼整體公式 

　　就是遍歷所有的像素點之后求均值。

實現：

卷積去解決  ，然后我們再求矩陣的均值

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

公式:

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 角二階矩

公式：

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. 峰值信噪比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