可解譯深度網絡的多光譜遙感圖像融合

余典，李坤，張瑋，李對對，田昕*，江昊

1.武漢大學電子信息學院，武漢 430072; 2.中國資源衛星應用中心，北京 100094

0 引 言

受限于傳感器光學系統設計的局限性，遙感圖像空間分辨率和光譜分辨率之間存在相互約束的關系。例如，多光譜圖像(multispectral image, MSI)具有較高的光譜分辨率和較低的空間分辨率，而全色圖像(panchromatic image, PAN)則具有較高的空間分辨率和單一光譜信息。多光譜圖像融合，又稱Pansharpening，將同一時區、同一目標的低空間分辨率MSI和高空間分辨率PAN，按照一定算法進行信息提取和運算結合，可以融合生成高空間分辨率的MSI，從而在光譜特性和空間結構方面均可以達到對成像目標更為準確地識別，極大提升了定量分析的精度和可靠性。近半個世紀以來，多光譜圖像融合可在同時搭載MSI和PAN的遙感衛星中直接應用(如IKONOS、QuickBird、Gaofen-2等)，在自然災害監測、軍事偵探、土地測量和農業分析等多個領域發揮出巨大應用價值(王海榮 等，2021；Pohl和van Genderen，1998；焦姣和吳玲達，2019)。

1 相關工作

目前代表性的多光譜圖像融合方法主要包括：成分替換方法(component substitution, CS), 多分辨率分析方法(multi-resolution analysis, MRA)、變分模型優化方法(variational model optimization, VO)以及深度學習方法(deep learning, DL)(胡鑫，2021)。CS方法首先通過矩陣分解將MSI分離出空間和光譜信息，進一步將PAN與MSI的空間信息疊加后進行相應逆變換得到融合MSI。Tu等人(2001)首先提出了一種基于IHS(intensity hue-saturation)空間映射的多光譜圖像融合的方法。進一步地，非線性IHS變換(Ghahremani和Ghassemian，2016)和非線性主成分分析(principle component analysis, PCA)變換(Licciardi等，2012)分別用于成分替換，從而有效提升融合性能。雖然CS方法可以有效增強空間細節，但是往往存在較強的光譜失真。MRA方法則是通過對MSI和PAN進行多尺度分解，在不同尺度下融合相應的分解系數以得到融合圖像。例如，采用非下采樣輪廓基函數方法進行多尺度分解，克服了小波基下采樣過程中帶來的失真問題(da Cunha等，2006)。其他常用的多尺度分解方法包括Laplacian金字塔分解方法(Do和Vetterli，2002)、Morphological filter方法(Aiazzi等，2006)和基于注入系數的方法(Garzelli等，2018；Restaino等，2020)。相對于CS方法，MRA方法可以獲得更佳的光譜信息，但是其空間細節清晰度相對較差。

VO方法首先建立融合圖像與源圖像MSI和PAN之間的函數關系(又稱先驗能量函數)，然后結合不同先驗約束構建融合模型，從而將圖像融合問題轉化為數學優化問題，并通過迭代優化獲得最優的融合結果。相對于CS方法和MRA方法，VO方法可以獲得更好的高分辨率空間與光譜信息。Ballester等人(2006)首先提出了一種P+XS的VO方法，假設PAN是融合圖像多波段線性組合的結果，而MSI由融合圖像下采樣形成，在此基礎上構建融合模型。為進一步提升融合性能，Fang等人(2013)采用了Guided filter描述圖像中的相似結構。M?ller 等人(2012)使用小波變換基作為函數約束項，提出了一種變分小波多光譜圖像融合方法。相對于通過強度構建PAN與融合圖像之間的關系而言，利用細節相似性進行描述可以獲得更佳的結果，例如，Chen等人(2014)、Li等人(2017)和Tian等人(2020)通過PAN與融合圖像間的稀疏梯度關系構建了多種融合模型。Chen等人(2015)采用L2,1范數描述PAN與融合圖像間的組梯度稀疏關系，在融合過程中有效考慮了配準問題，形成配準與融合統一模型。相比于全局梯度，局部梯度能更好地描述相似性，因此，Fu等人(2019)在融合模型構建上采用了局部梯度約束。

得益于深度學習方法強大的非線性特征表征能力，DL方法通過深度學習進行多光譜圖像融合逐漸成為近年來的研究熱點。Huang等人(2015)首先將深度神經網絡應用到圖像融合領域，通過不同分辨率的全色圖像進行訓練，生成全色圖像的退化模型，并假設多光譜圖像的退化過程也遵循這一模型，從而使用訓練得到的深度神經網絡重建融合圖像。Masi等人(2016)和Yuan等人(2018)分別采用卷積神經網絡(convolutional neural networks, CNN)和多尺度多深度卷積神經網絡提升融合性能。為了保留更多空間信息，Yang等人(2017)提出了一種名為PanNet(deep network architeeture for pansharpening)的融合方法。PanNet利用殘差網絡獲得易丟失的高頻信息，再注入融合圖像中，從而在保持光譜信息的同時有效提升空間細節信息的清晰度。He 等人(2019)基于CS方法的思想，將細節替換過程映射為深度神經網絡，有效減少了光譜失真。Shen等人(2019)使用深度殘差CNN構建VO方法中的梯度特征，從而提升VO方法融合模型的準確性。

從以上分析中可以看出，雖然VO方法和DL方法在多光譜圖像融合領域中取得了較好的效果，但還有如下難題需要解決：1)如何尋求準確的先驗構建有效的融合成像模型對于VO方法是非常重要的，傳統方法往往通過線性特征(例如梯度)構建先驗建立PAN與融合圖像的關系，難以描述自然場景復雜非線性關系，導致成像模型準確性依然有限；模型參數對VO方法影響巨大，傳統方法需要手動調參尋求最優模型參數，其調參過程往往過于耗時，且最優模型參數難尋。2)雖然DL方法可以通過非線性關系建立融合圖像之間的聯系，從而克服VO方法使用線性融合模型的局限性，但是，傳統DL方法往往將融合過程當做黑盒，忽略了真實物理成像意義，融合性能依然有待突破。Tian等人(2022)將VO方法中的融合模型求解過程映射為網絡架構，提出了一種物理可解譯的深度學習多光譜圖像融合方法，為解決上述問題提供了一種新的思路，但其利用L1范數構建先驗約束導致融合模型準確性依然存在提升空間。

本文提出了一種基于可解譯神經網絡的多光譜圖像融合方法。通過深度學習構建深度先驗建立融合圖像與PAN之間的聯系，由于該先驗由學習驅動所得，因此，可以有效提升融合模型先驗的準確性；從物理成像機理出發，將MSI看做是下采樣的融合圖像構建數據保真項，并結合上述深度先驗，形成一種新的多光譜圖像融合模型，此時，圖像融合問題轉化為上述融合模型的優化求解問題；通過近端梯度下降法(proximal gradient descent, PGD)對上述優化問題進行求解，進一步將優化求解步驟映射為一種可解譯深度網絡架構，通過訓練數據學習生成最優求解結果?？梢钥闯?，該方法的優點主要體現在：1)該網絡架構的每個模塊對應于求解過程的每個步驟，具有明確的物理可解譯性；2)非線性模型先驗和融合模型參數由學習所得，有效提升了成像模型的準確性，同時極大降低了傳統VO方法最優模型參數選擇的難題。

本文的主要貢獻為：

1)通過數據驅動形成深度學習映射，從而構建一種新的非線性學習先驗，有效提升模型的準確性；

2)基于上述先驗構建一種新的融合模型，進一步將融合模型求解過程映射為可解譯深度網絡架構，避免傳統VO方法的調參難題；

3)分別從主觀視覺和客觀量化分析兩個方面，在仿真和真實數據集上對上述方法的優越性進行了有效驗證。

2 本文方法

2.1 融合模型建立

所提出融合模型可以表示為下述優化問題

(1)

式中，J(F)和K(F)分別代表數據保真項和先驗約束項。

根據物理成像機制，MS圖像M可視為融合圖像F經模糊和空間下采樣的結果，即M=ΨF+N。Ψ表示模糊和空間下采樣運算符，N表示均值為0的高斯噪聲。因此，為了在融合圖像中有效保持MS圖像中固有的光譜特性，常用數據保真項J(F)可以表示為

(2)

為了將PAN圖像中的空間細節信息有效傳遞至融合圖像中，傳統方法如Chen等人(2015)和Tian 等人(2022)用L1范數構建先驗約束描述PAN與融合圖像間的殘差關系，但由于自然場景特征規律往往較為復雜，難以用僅僅適宜于表征稀疏特性的L1范數去準確描述。為了解決上述問題，擬將其殘差關系通過數據驅動由深度學習函數φ進行描述，在此基礎上通過深度殘差構建先驗約束項K(F)，即

(3)

基于上述分析，所提出融合模型可以表示為

(4)

式中，λ為平衡數據保真項和先驗約束項的參數。

2.2 模型優化求解

為了求解式(4)，采用PGD方法。引入中間變量Z∈Rm×n×b，此時，式(4)可以轉化為如下兩個公式進行迭代求解，即

Zt+1=Ft-μΨT(ΨFt-M)

(5)

(6)

通過式(6)發現，Zt+1可以看做F在第t+1次迭代時的重建結果。因此有如下假設：在第t+1次迭代時，F-Zt+1的殘差服從均值為0、方程為δ2的獨立正態分布(Zhang等，2014)?；谏鲜黾僭O，可以通過與Zhang和Ghanem(2018)相似的推導，得到如下結論：

定理1：假設x1,…,xk∈Rl是均值為0、方差為δ2的獨立同分布變量，X=[x1,…,xk]T。對于任意矩陣A∈Rp×l，可以得到如下結論：

基于定理1，可以得到如下近似

(7)

因此，式(6)可以轉化為

(8)

(9)

可以通過Zhang等人(2022)方法從式(9)中求解出Wt+1，進一步地，可以得到Ft+1的近似求解，即

(10)

基于上述分析，所提出模型的優化求解過程可以總結為：

算法1 基于近端梯度下降的模型優化求解算法。

輸入：M，P，Ψ，Rp，μ，ρ。

初始化：F0。

循環t=0至最大值：

根據式(5)更新Zt+1。

根據式(9)更新Wt+1。

根據式(10)更新Ft+1。

結束。

輸出：F。

2.3 可解譯深度網絡

將算法1所示的優化迭代步驟展開成多個神經網絡模塊，從而形成所需的模型可解譯融合網絡。一方面，利用深度網絡的學習能力挖掘全色圖像和融合圖像之間的結構先驗，使得模型更為逼近真實的物理成像關系；另一方面，網絡中參數優化不僅受訓練數據驅動，同時受物理模型的優化機制所引導，提升訓練模型的融合精度；最后，將模型中涉及的超參數設置嵌入在網絡模塊中，在數據訓練中實現超參數自主尋優，避免復雜的手工設計問題。

圖1 融合網絡中第t個網絡狀態模塊

圖2 模塊2中SwinResUnet的網絡結構圖

網絡模塊2的細節網絡結構圖如圖3所示，受Liu 等人(2021)和Zhang等人(2022)方法的啟發，本文方法采用結合Transformer和CNN架構的網絡模塊SwinResUnet(Zhang等，2022)作為式(9)的去噪器。具體來講，SwinResUnet采用U-Net(Ronneberger等，2015)作為網絡主體架構，主要包含3組不同維度下的編碼器(E1,E2,E3)和解碼器(D1,D2,D3)，中間嵌入特征提取單元F，同時，對應的編碼器(Ei)和解碼器(Di)之間建立殘差短連接結構，增強不同尺度下的特征傳遞效果。因此，編碼器與解碼器單元(Ei與Di)特征維度和卷積通道數均一致。E1，E2，E3和F中的1×1卷積單元通道數分別為64，128，256和512。

進一步地，圖3(a)展示了所采用的編碼器和解碼器內部網絡細節。如圖3(a)所示，編碼器(E1,E2,E3)包含黑色虛線之間的網絡結構，特征提取單元F包含藍色虛線之間的網絡結構，而解碼器(D1,D2,D3)包含紅色虛線之間的網絡結構。三者結構共享特征提取單元F部分。其中，編碼單元在F結構末端加入了大小為2×2(步長為2)的卷積操作實現特征尺寸的下采樣；而解碼單元在F結構前端加入了大小為2×2(步長為2)的轉置卷積操作實現特征尺寸的上采樣，以此實現數據的多層語義與結構信息提取融合。對于特征提取單元F，本文方法引入了RConv所表示的殘差卷積模塊；此外，考慮到全局感受野下不同區域特征的聯系，本文方法在RConv單元并行加入了Swin Transformer(SwinT)網絡單元。得益于其內部的多頭自注意力機制(multi-head self attention, MSA)以及窗口遷移等操作，該結構能夠在保證局部光譜和結構特性的情況下，從圖像全局層面增強不同區域間空間和光譜信息的交互和融合。在E1,E2,E3和F中，RConv分別采用通道數為32，64，128和256的3×3卷積核對等分后的特征進行處理，SwinT則采用全連接層進行隱式表征，并分別將結果級聯后還原為通道數為64，128，256和512的特征。同時，本文方法在每個編碼器和解碼器的后端和前端分別設計了空間下采樣和上采樣卷積操作(SConv和TConv)，以此實現不同層級的特征維度變換，獲取兼顧高層語義特征和底層細節特征的信息。其中，SConv和TConv分別使用步長和卷積核大小均為2的卷積和轉置卷積實現。圖3(b)表示的是SwinT和RConv單元網絡結構(WMSA(window based multi-head self-attension),SWMSA(shifted WMSA))，可以看到，SwinT網絡單元由兩個相同的結構單元組成，均包含一個基于窗口的MSA基本單元，以及由兩個線性層和介于二者之間的GELU(Gaussian error linear units)非線性層構成的MLP(multilayer perception)基本單元。同時在MSA和MLP前側使用LN(LayerNorm)層進行標準化并進行殘差跳躍連接。圖3(b)中RConv則包含兩個3×3卷積層以及介于二者之間的ReLU非線性層，同樣使用了殘差跳躍連接結構。

圖3 編碼器和解碼器內部網絡細節

本文采用L1損失函數對網絡進行約束并驅動參數優化，即

(11)

3 實驗結果與分析

3.1 實驗設置

實驗主要從仿真實驗、真實實驗和消融分析3個方面展開。仿真實驗主要基于Wald仿真協議：將高分辨率的MSI和PAN分別經過MTF(modulation transfer function)濾波和空間下采樣，生成低分辨率的MSI和PAN用于多光譜圖像融合實驗，將該高分辨率的MSI作為真值圖像用于對融合結果進行定性和定量分析。MTF濾波器與理想低通濾波器相似，在Nyquist頻率處具有截止幅值。對比真值圖像與融合圖像的差異，可以首先從主觀視覺方面對不同方法的特點進行定性分析。進一步地，可以采用全參考圖像質量評價指標對算法的有效性進行定量評估，具體評價指標包括全局相對無量綱誤差(relative dimensionless global error in synthesis, ERGAS)、光譜角映射(spectral angle mapping, SAM)、全局綜合評分Q2n、結構相似度(structural similarity index, SSIM)、均方根誤差(root mean square error, RMSE)、相對平均光譜誤差(relative average spectral error, RASE)、通用圖像質量指數(universal image quality index, UIQI)和峰值信噪比(peak signal-to-noise ratio, PSNR)。ERGAS和Q2n是多光譜圖像融合綜合性能評價指標。SAM通過計算融合圖像與真值圖像對應的兩個矢量間絕對角度來測量融合圖像的光譜失真情況。SSIM用于衡量圖像之間的結構信息相似程度。RASE用于評估融合圖像的全局光譜質量。PSNR和RMSE主要從像素差異角度衡量融合結果與真實結果之間的偏差。UIQI用于評價融合圖像與參考圖像的結構失真程度。

真實實驗直接使用原始的MSI和PAN進行融合，用原始數據尺寸進行真實數據的視覺比較和定量分析。由于真實實驗中并沒有真值圖像，本文采用無參考評價指標QNR(quality with no reference)，Ds和Dλ。QNR是一種綜合性評價指標，它由空間失真指標Ds和光譜失真指標Dλ所組成(Yang等，2017)，計算為

fQNR=(1-Dλ)α(1-Ds)β

(12)

式中，α和β為系數，默認取值為1。

對比方法包括一種通過23個系數進行多項式內插的插值方法EXP；兩種CS方法：GSA(gram-schmidt adaptive)(Aiazzi等，2007)和PRACS(partial replacement-based adaptive component substitution)(Choi等，2011)；兩種MRA方法：Indusion(induction scaling)(Khan等，2008)和AWLP(additive wavelet luminance proportional)(Vivone等，2015)；兩種VO方法：SIRF(simultaneously registration and fusion)(Chen等，2015)和LGC(local gradient constrants)(Fu等，2019)；4種深度學習方法：PNN(pansharpening by using a convolutional neural network)(Masi等，2016)，PanNet(Yang等，2017)，VPNet(interpretable deep network for variational pansharpening)(Tian等，2022)和本文方法。為了公正比較，實驗所有對比方法均運行在配有英特爾至強核心W-2200CPU@3.70 GHz和11 GB顯存的GeForce GTX 2080Ti GPU的計算機上。所有對比方法均采用原始默認參數，同時深度學習方法如PNN、PanNet和VP-Net等均使用相同的數據處理方式和數據集。本文采用Adam優化器更新參數，初始學習率的值設置為0.000 5，訓練batch大小為16，經過400訓練迭代次數達到了穩定的性能。

3.2 仿真實驗

本文仿真實驗主要采用GeoEye-1和Gaofen-2遙感衛星提供的兩組數據集進行多光譜圖像融合實驗驗證。

3.2.1 Gaofen-2遙感衛星數據仿真實驗

圖4 Gaofen-2仿真數據集融合圖像主觀視覺比較(選用紅、綠、藍3個波段顯示)

根據融合圖像與真值圖像之間的平均絕對誤差生成殘差圖像，可以驗證實驗效果。其結果如圖5所示。由于存在空間模糊或光譜失真，EXP，GSA，PRACS，Indusion和AWLP(圖5(a)—(e))具有較大的殘差。通過構建融合模型并進行求解，VO方法可以獲得更好的空間與光譜信息的平衡，因此，SIRF和LGC相對于上述方法而言殘差較小。通過深度學習，PNN、PanNet、VPNet和本文方法具有較好的融合效果，因此，可以有效地減少殘差。相對于其他方法而言，本文方法具有最小的殘差，從而表明了本文算法在Gaofen-2遙感衛星仿真數據上具有最好的融合性能。

圖5 Gaofen-2仿真數據集融合結果殘差圖像比較

定量評估方面，通過28幅測試圖像進行統計分析，多種評價指標的平均值如表1所示。從表1中可以看出，在GSA，Indusion和AWLP方法中，光譜失真評價指標SAM和RASE表現較差，導致其綜合評價指標ERGAS和Q2n與其他方法尚有差距，例如，在上述方法中ERGAS大于5。PRACS，SIRF和LGC可以在光譜失真和空間細節增強方面達到較好的平衡，所以具有較優的綜合評價指標(ERGAS在4.0左右)?？傮w而言，深度學習方法PNN，PanNet，VPNet和本文方法相對于其他方法而言，性能提升顯著。例如，對比于LGC，PNN的ERGAS可以有效減少14.45%。本文方法在所有指標上均具有最佳的融合效果，例如，相對于性能第2的VPNet方法，ERGAS可以有效減少7.58%，證明了其在定量評估實驗中相對于其他方法的優越性。

表1 基于Gaofen-2仿真數據集28幅測試圖像的定量評估結果

3.2.2 GeoEye-1遙感衛星數據仿真實驗

圖6 GeoEye-1仿真數據集融合圖像主觀視覺比較(選用紅、綠、藍3個波段顯示)

圖7 GeoEye-1仿真數據集融合結果殘差圖像比較

定量評估方面，通過40幅測試圖像進行統計分析，多種評價指標的平均值如表2所示。與表1類似，深度學習方法PNN，PanNet，VPNet和本文方法，相對于其他方法而言，各種性能指標提升均較為顯著。例如，PNN相對于GSA，ERGAS可以降低19.94%。

表2 基于GeoEye-1仿真數據集40幅測試圖像的定量評估結果

總體而言，本文方法在該數據集上優勢也較為顯著，例如，相對于性能第2的VPNet方法，ERGAS可以有效減少4.61%，因此，在7個客觀評價指標上均獲得了最佳的效果。

3.3 真實實驗

圖8 QuickBird真實數據集融合圖像主觀視覺比較(選用紅、綠、藍3個波段顯示)

定量評估方面，本文選取了一組典型測試圖像進行統計分析，QNR，Ds和Dλ等無參考評價指標的平均值如表3所示。由于本文方法具有最優的空間清晰度和較優的光譜保真度，所以獲得最佳的空間評價指標Ds和較好的光譜評價指標Dλ，從而產生最佳的綜合無參考融合性能評價指標QNR，有效證明了本文方法在QuickBird真實數據集上相對于其他方法的優越性。

表3 基于QuickBird真實數據集的定量評估

3.4 消融分析

在本文可解譯的融合網絡設計過程中，模型迭代次數T(網絡狀態模塊個數)和超參數ρ等初始化設置均對融合性能產生直接的影響，為此，本小節針對以上重要因素對融合網絡分別進行了消融實驗分析，從而得到最適合的模型架構和網路參數設置，提升多光譜和全色圖像的融合質量。

表4給出了不同模型迭代次數下對應的融合網絡在GeoEye-1數據集上的測試性能，可以看出，隨著T的增加，網絡訓練也隨之增加，且融合指標ERGAS和RMSE在T小于5時下降很快，當T大于5時兩項指標下降幅度較小，即隨著迭代次數T增加，模型性能提升不再明顯，綜合訓練時間成本和融合性能的權衡考慮，本文將模型迭代次數T設置為5。

表4 不同迭代次數下的模型性能對比

此外，式(8)中權重參數ρ對于融合結果質量有重要影響，為此，本文設置了多組不同ρ初始化值，并在GeoEye-1數據集進行相同設置下的訓練，最后將測試結果統計在表5中?？梢园l現，當其初始化值設為0.001時綜合指標ERGAS取得相對最優的測試結果，光譜保真度指標SAM同樣表現最佳。因此，本文方法采用ρ初始值為0.001，從而得到兼顧光譜和空間質量的融合結果。

表5 不同ρ下的融合性能對比

4 結 論

面向多光譜遙感圖像融合，本文提出了一種結合VO和DL方法優點的可解譯深度網絡。首先基于深度先驗構建了一種新的多光譜融合模型。接著為了求解融合模型，通過PGD方法將求解過程映射為多個迭代步驟，進一步通過深度展開技術將上述步驟映射為深度網絡模塊，形成一個具有物理可解譯性的深度網絡架構。其優勢在于：由于深度先驗具有更好的非線性表征能力，所以可以有效提升融合模型的準確性。同時，由于模型求解參數由深度網絡訓練生成，降低了傳統VO方法參數設置難度。

為了有效驗證本文方法相對于傳統方法的有效性，分別在仿真和真實遙感衛星數據集上進行了主客觀對比實驗。從主觀實驗中可以看出，本文方法相對于其他方法而言，在保持光譜信息的同時有效增強了融合圖像的空間細節信息?？陀^實驗分別使用了全參考評價指標ERGAS，SAM，Q2n，SSIM，RMSE，RASE，UIQI及PNSR和無參考評價指標QNR，Dλ和Ds進行評價。大量統計實驗表明了本文方法在上述指標上的優越性。

如何在更多的真實衛星數據上驗證本文方法的有效性是未來需要考慮的問題。另外，將本文方法應用于高光譜與多光譜圖像融合，也是未來研究方向之一。