Hyperspectral Pansharpening Based on Homomorphic Filtering and Weighted Tensor Matrix

: Hyperspectral pansharpening is an e ﬀ ective technique to obtain a high spatial resolution hyperspectral (HS) image. In this paper, a new hyperspectral pansharpening algorithm based on homomorphic ﬁltering and weighted tensor matrix (HFWT) is proposed. In the proposed HFWT method, open-closing morphological operation is utilized to remove the noise of the HS image, and homomorphic ﬁltering is introduced to extract the spatial details of each band in the denoised HS image. More importantly, a weighted root mean squared error-based method is proposed to obtain the total spatial information of the HS image, and an optimized weighted tensor matrix based strategy is presented to integrate spatial information of the HS image with spatial information of the panchromatic (PAN) image. With the appropriate integrated spatial details injection, the fused HS image is generated by constructing the suitable gain matrix. Experimental results over both simulated and real datasets demonstrate that the proposed HFWT method e ﬀ ectively generates the fused HS image with high spatial resolution while maintaining the spectral information of the original low spatial resolution HS image.


Introduction
Depending on the number of acquired bands, remote sensing imaging technology has developed from collecting panchromatic (PAN) and color images to multispectral (MS) images, and it can now capture hyperspectral (HS) images with dozens of hundreds of bands. A PAN image with very high spatial resolution is a single-band grayscale image acquired in the visible range. It is able to obtain the shape feature of objects, but cannot distinguish colors. A color image consists of three bands which are red, green and blue, and displays the colors of objects. However, it is difficult to distinguish the features in similar colors. An MS image not only obtains spatial features, but also obtains spectral information in several bands, which is more capable of distinguishing categories of different features. However, the rough spectral resolution of MS images may not meet the requirements in some applications, and it is hard to realize fine feature detection [1]. An HS image with a higher spectral resolution on the order of nanometers can provide finer classification [2], which has been applied to many fields [3][4][5][6][7] and some practical applications, such as vegetation study [8], precision agriculture [8], regional geological mapping [9], mineral exploration [10], and environment monitoring [11]. Due to technical limitations, the spatial resolution of an HS image is low. Many hyperspectral pansharpening algorithms were developed, among which hyperspectral pansharpening methods using Bayesian and matrix factorization have been proposed in recent years. The Bayesian-based approaches include Bayesian naive Gaussian prior [12], Bayesian sparsity promoted Gaussian prior [13], and HySure [14]. These algorithms utilize the posterior distribution, and are based on maximum a posteriori estimation to fuse LRHS and HRPAN images [15]. The matrix factorization approach generates a fused HRHS image by using the nonnegative matrix factorization (NMF) under some constraints to estimate endmember and abundance matrices [16]. The matrix factorization approach is well represented by the nonnegative sparse coding (NNSC) [17] and constrained nonnegative matrix factorization (CNMF) [18] methods. The main challenge in hyperspectral pansharpening is to effectively improve the spatial resolution while preserving the original spectral information. The Bayesian and matrix factorization approaches are able to achieve good results on this challenge, but have a high computational cost.
Component substitution (CS) and multi-resolution analysis (MRA) approaches are two classical hyperspectral pansharpening approaches which have simple and fast implementation. For the CS class, intensity-hue-saturation (IHS) transform [19,20], principal component analysis (PCA) transform [21,22], Gram-Schmidt (GS) [23], and adaptive GS (GSA) [24] are the most representative methods. The CS class extracts spatial details of the HS image, and replaces the extracted spatial details with the HRPAN image. Regardless of superior spatial performance, the CS class suffers from serious spectral distortion [25]. The typical algorithms of the MRA technique are smoothing filter based intensity modulation (SFIM) [26], Laplacian pyramid [27], modulation transfer function generalized Laplacian pyramid (MTF-GLP) [28], and MTF-GLP with high pass modulation (MTF-GLP-HPM) [29]. The MRA methods generally utilize a multi-resolution decomposition to extract spatial details which are imported into the HS image. Compared with the CS methods, the MRA methods generate less spectral distortion, but usually have a larger computational burden [30]. Recently, several algorithms based on the CS and MRA approaches have been proposed, such as the Sentinel-2A CS and MRA based sharpening algorithm [31], the multiband Filter estimation (MBFE) algorithm [32], and the guided filter PCA (GFPCA) algorithm [33]. Moreover, several intelligent processing-based methods have also been proposed, and examples include deep two-branches convolutional neural network (Two-CNN-Fu) [34], Bidirectional Pyramid Network [35], and 3D-convolutional neural network (3D-CNN) [36]. The CS and MRA approaches mostly extract the spatial information of the HRPAN image and inject it into the LRHS image, but without considering the spatial information of the LRHS image. Due to the incomplete spatial information injection, the CS and MRA approaches may result in distortion. To address this problem, we propose a novel hyperspectral pansharpening method by combining homomorphic filtering with a weighted tensor matrix. An optimized weighted tensor matrix-based method which considers the structure information of the LRHS and HRPAN images is proposed to generate more comprehensive spatial information. In addition, to extract the spatial structure information of the LRHS images, open-closing morphological operation is first used for noise removal, and homomorphic filtering is then introduced to extract the spatial details of each band. Finally, a weighted root mean squared error based method is proposed to obtain the total spatial component of the LRHS image from extracted spatial details of each band, and the Laplacian pyramid networks super-resolution algorithm is adopted to enhance the spatial resolution of the obtained spatial component. Comparative analysis was used to demonstrate the applicability and superiority of the proposed method in both spectral and spatial qualities.
As stated above, a new hyperspectral pansharpening method based on homomorphic filtering and weighted tensor matrix is proposed in this paper. The main novelties of the proposed hyperspectral pansharpening method are concluded in the following aspects.

1.
A novel HS image spatial component extraction strategy is proposed. Open-closing morphological operation and homomorphic filtering are first introduced to remove the noise and extract the spatial details of each band of the HS image, respectively. Then, a weighted root mean squared error-based method is proposed to obtain the total spatial component of the HS image.

2.
An optimized weighted tensor matrix-based method is proposed to integrate the spatial component of the HS image with the spatial component of the PAN image. The weighted structure tensor matrix that represents the structural information of multiple images is applied to hyperspectral pansharpening for the first time. The classical methods which mostly extract the spatial information of the PAN image inject the incomplete spatial information, and may lead to distortion. Unlike there classical methods, the proposed optimized weighted tensor matrix-based method generates the spatial information not only from the PAN image but also from the HS image, and can reduce the distortion caused by the insufficient spatial information.
The remainder of this paper is organized as follows. Section 2 describes the weighted structure tensor matrix and homomorphic filtering. In Section 3, the proposed homomorphic filtering and weighted tensor matrix-based hyperspectral pansharpening algorithm is presented. Experimental results and discussion are provided in Section 4, and conclusions are drawn in Section 5.

Weighted Structure Tensor Matrix
For an image I, the structure tensor matrix M can be decomposed as: where I x = ∂I/∂x and I y = ∂I/∂y are the horizontal and vertical partial derivatives of the image, ∇I = [ I x I y ] T , () T is the transpose operation, v 1 , v 2 , e 1 and e 2 are the two eigenvectors and the corresponding eigenvalues, respectively. As shown in Equation (1), the tensor matrix M which is a symmetric and semi-definite positive matrix has eigen-decomposition, and it has been exploited in some fields, such as texture synthesis [37], image regularization [38], denoising [39], and recognition systems [40]. The eigenvalues obtained by decomposing the tensor matrix are utilized to describe the where M w is the weighted structure tensor matrix, M l is the structure tensor matrix of the lth image, ∂I l /∂x and ∂I l /∂y are the horizontal and vertical partial derivatives of the lth image, respectively. The weighted tensor matrix M w is also a symmetric and semi-definite positive matrix, and can be proceeded the eigen-decomposition.

Homomorphic Filtering
Homomorphic filtering which is a type of frequency domain filtering can compress the image brightness range and enhance the image contrast. Homomorphic filtering has been applied to some image processing problems [41][42][43] and is based on an image imaging model: where f represents an image, f H represents the high frequency reflectance component, and f L represents the low frequency illumination component. Homomorphic filtering aims to reduce the low frequency component of an image. Logarithmic transformation is utilized to separate the two components: After applying the Fourier transform: where F, F H and F L denote the Fourier transform of ln( f ), ln( f H ) and ln( f L ), respectively. Then, the high-pass filter H is applied to Equation (5) as where S is the filtered result. The final image is obtained by the inverse Fourier transform and the exponential operation: where f h f denotes the homomorphic filtered image, s denotes the inverse Fourier transform of S, and −1 denotes the inverse Fourier transform.  The open-closing operation which belongs to the mathematical morphology operation is an effective denoising processing operation [44,45]. The effects of denoising using the open operation and closed operation alone are usually not very good, since they may cause amplitude deflection. By contrast, the open-closing operation has the better denoising effect. The open operation is first applied on the image, and the selected structure element is larger than the noise size to remove the background noise. Then, the closed operation is utilized to remove the noise of the image obtained in the previous step. The open-closing denoising operation is suitable for the images which have less small details. Since the LRHS image has the low spatial resolution, the fine spatial details are few. The open-closing operation is applicable to removing noise with high interference in the LRHS image. The open-closing morphological operation is applied as:   The open-closing operation which belongs to the mathematical morphology operation is an effective denoising processing operation [44,45]. The effects of denoising using the open operation and closed operation alone are usually not very good, since they may cause amplitude deflection. By contrast, the open-closing operation has the better denoising effect. The open operation is first applied on the image, and the selected structure element is larger than the noise size to remove the background noise. Then, the closed operation is utilized to remove the noise of the image obtained in the previous step. The open-closing denoising operation is suitable for the images which have less small details. Since the LRHS image has the low spatial resolution, the fine spatial details are few. The open-closing operation is applicable to removing noise with high interference in the LRHS image. The open-closing morphological operation is applied as:

Hyperspectral Image Preprocessing
band of the LRHS image and the denoised LRHS image, respectively, and S 1 and S 2 are the structure elements. Here, • represents the opening operation which first uses the erosion operation and then the dilation operation, and • denotes the closing operation which does in reverse. The erosion and dilation operations obtain the local minimum and maximum of the image, respectively. Equation (8) can be expressed in detail as: Remote Sens. 2019, 11, 1005 6 of 18 for k = 1, 2, . . . , B, where Θ and ⊕ denote the erosion and dilation operations, respectively.

Hyperspectral Image Spatial Information Extraction
Homomorphic filtering is a filtering method that transforms the nonlinear problem into a linear problem. It transforms the nonlinear multiplicative mixed problem into an additive model by the logarithmic transformation, and then uses linear filtering to process it. The homomorphic filtering suppresses the low frequency illumination component and enhances the high frequency reflectance component. For an HS image, the high frequency component of each band is considered as the spatial component for each band. To obtain the spatial information for each band, we apply homomorphic filtering to each band of the denoised LRHS image. Through the use of homomorphic filtering, the low frequency component of each band of the denoised LRHS image is suppressed, and the high frequency component is extracted. Therefore, in this research, homomorphic filtering is applied to each band of the denoised LRHS image to extract the spatial component of each band. The homomorphic filtering processing is based on the following image imaging model: for k = 1, 2, . . . , B, where X RNH LR_H represents the high frequency component of the denoised LRHS image, X RNH LR_L represents the low frequency component, and (X RNH LR_H ) k and (X RNH LR_L ) k represent the kth band of X RNH LR_H and X RNH LR_L , respectively. Based on Equation (4)-(6), Logarithmic transformation, Fourier transform, and high-pass filtering operations are applied to Equation (11): for k = 1, 2, . . . , B, where S LR is the high-pass filtered imageh, (S LR ) k is the kth band of S LR , represents Fourier transform, and H is the high-pass filter, defined as: where D 0 is the cut-off frequency, D is the distance between (x, y) and the center, β H and β L are the high and low frequency gains. Figure 3 shows the 3-D mesh of the high-pass filter. Since homomorphic filtering aims to reduce the low frequency component and extract the high frequency component, β H is greater than 1 and β L is smaller than 1. By adjusting the value of the cut-off frequency D 0 , the sharpness of the transition between β L and β H can be controlled. In practice, the values of these parameters are generally determined empirically. In this paper, empirically, β H , β L , and D 0 are set to 2, 0.25, 40, respectively. S LR is the high-pass filtered image in which the low frequency component has been weakened. Then, the spatial component of each band is obtained by applying the inverse Fourier transform and the exponential operation to S LR .
for k = 1, 2, . . . , B, where X I LR denotes the spatial component of each band of the denoised LRHS image, (X I LR ) k denotes the kth band of X I LR , and −1 denotes the inverse Fourier transform. After introducing homomorphic filtering to obtain the spatial component of each band of the denoised LRHS image, a weighted root mean squared error (RMSE)-based method is presented to extract the spatial intensity information of the HS image. Let I LR = B k=1 λ k (X I LR ) k denote the total spatial information of the LRHS image, where [λ 1 , λ 2 , . . . , λ B ] is the weighted vector. To determine the values of the weighted vector, we utilize the RMSE index to measure the deviation of two images. A smaller value of RMSE indicates a better result, and the optimal value is 0. In the RMSE-based method, the spatial information of the PAN image is considered, and the RMSE value between the total spatial information I LR and the PAN image X PAN is calculated. The smallest value of RMSE is computed to obtain the optimal values of the weights [λ 1 , λ 2 , . . . , λ B ]: where T = m × n represents the total pixel number of one band of the LRHS image, ↓ represents down-sampling operation, ↓ X PAN represents that the PAN image is down-sampled to the size of k , respectively. The laplacian pyramid networks (LapSRN) [46] super-resolution method can effectively improve the spatial resolution of an image, and has the advantages of parameter sharing, local skip connections, and multi-scale training. So it is adopted to super-resolve the spatial information of the LRHS image I LR for an I HR with super-resolution spatial information.

Hyperspectral Image Spatial Information Extraction
Homomorphic filtering is a filtering method that transforms the nonlinear problem into a linear problem. It transforms the nonlinear multiplicative mixed problem into an additive model by the logarithmic transformation, and then uses linear filtering to process it. The homomorphic filtering suppresses the low frequency illumination component and enhances the high frequency reflectance component. For an HS image, the high frequency component of each band is considered as the spatial component for each band. To obtain the spatial information for each band, we apply homomorphic filtering to each band of the denoised LRHS image. Through the use of homomorphic filtering, the low frequency component of each band of the denoised LRHS image is suppressed, and the high frequency component is extracted. Therefore, in this research, homomorphic filtering is applied to each band of the denoised LRHS image to extract the spatial component of each band. The homomorphic filtering processing is based on the following image imaging model: represents Fourier transform, and H is the high-pass filter, defined as:

Panchromatic Image Preprocessing and Total Spatial Information Acquisition
To make the spatial information of the PAN image clearer, the Laplacian of Gaussian (LOG) [47] image enhancement algorithm is applied to the PAN image, which uses a Gaussian filter to reduce noise followed by a Laplace operator for enhancement. Let I PAN s represent the enhanced PAN image. The HS and PAN images contain the different and complementary information for a scene. To acquire the total spatial information, we should consider simultaneously the spatial structure details of these two images, and then propose an optimized weighted tensor matrix-based method. I HR and I PAN s include the spatial structure information of the HS and PAN images, respectively. Based on Equation (2) where () T is the transpose operation, v w1 and v w2 are the two eigenvalues, v w1,p and v w2,p are the two eigenvalues at pixel p, e w1,p = [ e w11,p e w12,p ] T and e w2,p = [ e w21,p e w22,p ] T are the eigenvectors corresponding to the two eigenvalues at pixel p, respectively. The two eigenvalues generally have a larger value and a smaller value. We assume that v w1 is the larger eigenvalue. When v w1 ≈ v w2 ≈ 0, v w1 > v w2 ≈ 0, and v w1 > v w2 > 0, the structure region of this pixel are flat area, edge area, and corner, respectively. We test on some images to study the eigenvalues of weighted tensor matrix. Figure 4 shows the two eigenvalues at each pixel of the weighted tensor matrix for the Salinas scene data. It can be seen that for many pixels, the smaller eigenvalues shown in Figure 4d are approximately 10 −5 , and are very small. By experimenting on other numerous images, we have also discovered that the smaller eigenvalues are mostly very small. Thus, the approximation of M HP w,p is expressed as: where M HP w,p is the approximation of M HP w,p . Based on Equation (1), the structure tensor matrix satisfies that M = ∇I·∇I T . The weighted gradient G w at pixel p satisfies that M HP w,p = G w,p ·(G w,p ) T = v w1,p e w1,p (e w1,p ) T . So, G w,p is deduced as G w,p = √ v w1,p ·e w1,p . Since the direction of the eigenvector corresponding to v w1 is not unique, the direction of the weighted gradient G w,p is also not unique. We specify the direction of the weighted gradient G w,p as the gradients average of the individual multiple source images [I HR , I PAN s ]: where ∇I HR,p = [ ∂(I HR,p )/∂x ∂(I HR,p )/∂y ] and ∇I PAN s,p = [ ∂(I PAN s,p )/∂x ∂(I PAN s,p )/∂y ], ·, · represents the inner product of two vectors, and sign(·) represents the sign function. Once the weighted gradient G w,p is acquired from the multiple images [I HR , I PAN s ], an optimization model is proposed to obtain the total spatial information I HP T as: where ∇I HP T = [ ∂I HP T /∂x ∂I HP T /∂y ], ∂I HP T /∂x, and ∂I HP T /∂y denote the x and y partial derivatives of I HP T . Equation (20) is an unconstrained optimization problem, and we solve it by the conjugate gradient method. Equation (20) can effectively ensure that the total spatial information I HP T contains the spatial structure details of both the HS and PAN images.  The two eigenvalues generally have a larger value and a smaller value. We assume that 1 w v is the larger eigenvalue. When > , the structure region of this pixel are flat area, edge area, and corner, respectively. We test on some images to study the eigenvalues of weighted tensor matrix. Figure 4 shows the two eigenvalues at each pixel of the weighted tensor matrix for the Salinas scene data. It can be seen that for many pixels, the smaller eigenvalues shown in Figure 4d are approximately -5 10 , and are very small. By experimenting on other numerous images, we have also discovered that the smaller eigenvalues are mostly very small.
Thus, the approximation of HP , w p M is expressed as:

Fused High Spatial Resolution Hyperspectral Image Generation
The LRHS image X HS LR is interpolated to the scale of the HRPAN image. By constructing a suitable gain matrix R, the total spatial information I HP T is injected into the interpolated HS image to generate the fused HRHS image X HS HR . For the gain matrix R, it is beneficial to preserve the ratio of each HS pair band unchanged to reduce the spectral distortion. Thus, R should satisfy that R k ∝ where X HS IN is the interpolated HS image, (X HS HR ) k and R k are the kth band of X HS IN and R, respectively. Then, a tradeoff parameter ε is defined to regulate the amount of the injected details to reduce the spatial distortion. This process can be expressed as: for k = 1, 2, . . . , B, where X HS HR is the fused HRHS image, and (X HS HR ) k is the kth band of X HS HR .

Datasets and Experimental Setup
In order to evaluate the effectiveness of the proposed HFWT hyperspectral pansharpening method (named as HFWT), experiments were performed on two simulated hyperspectral datasets which were a Washington DC and a Salinas scene, and one real hyperspectral dataset, the Hyperion dataset. The Salinas scene hyperspectral dataset was collected by Airborne Visible Infrared Imaging Spectrometer (AVIRIS) [48], and the Washington DC dataset was acquired by the Spectral Information Technology Application Center of Virginia. The used real dataset is provided by the EO-1 spacecraft. The EO-1 spacecraft has a Hyperion instrument which provides the real LRSH images and an Advanced Land Imager (ALI) instrument acquires the HRPAN images [48]. Table 1 lists the characteristics of each dataset. The proposed HFWT method is compared with several state-of-the-art hyperspectral pansharpening methods: Gram-Schmidt (GS) [23], guided filter principal component analysis (GFPCA) [33], coupled nonnegative matrix factorization (CNMF) [18], Bayesian sparsity promoted Gaussian prior (Bayesian) [13], and HySure [14]. Four typical quantitative evaluation indexes are adopted: cross correlation (CC) [49], spectral angle mapper (SAM) [50], root mean squared error (RMSE), and erreur relative globale adimensionnelle desynthèse (ERGAS) [51]. The CC and SAM measure the spectral and spatial distortion, respectively. The larger value of CC and the smaller value of SAM indicate the better fusion result. The RMSE and ERGAS are the global indexes that measure both the spatial and spectral performance, and their value ranges are all (0,1), with 0 being the optimal value.
In order to perform the objective fusion evaluation of the simulated hyperspectral datasets, the available HS image is used as the reference HS image. The simulated LRHS image and HRPAN image are generated according to Wald's protocol [52,53]. The reference HS image is blurred and down sampled 4 times to obtain the simulated HS image. The simulated PAN image is obtained by averaging the visible light band of the reference HS image.
For the real datasets, the real LRHS and HRPAN images are available, and the reference high resolution HS image is not available. In order to test the objective quality of the real hyperspectral images, the real LRHS image is served as the reference image. The real LRHS and HRPAN images available are degraded, and the two degraded images are fused to obtain a fused image. This fused image is compared to the real LRHS image to evaluate the objective quality.
In the proposed HFWT method, we define a tradeoff parameter ε which regulates the amount of the injected details and ensures spatial performance. In practice, the optimal value is determined based on experience. By adjusting different values of ε, the optimal value can be determined by the fusion result. In this paper, by experience, the values of tradeoff parameter ε are set as 0.25, 0.05 and 0.2 for the Washington DC, Salinas scene and Hyperion dataset, respectively.

Validity Discussion of the Open-Closing Denoising Operation
To verify the effectiveness of the open-closing HS image denoising operation, the proposed HFWT method was conducted on the Washington DC dataset with different denoising processing. The compared denoising algorithms contain average filtering, Gaussian filtering, open operation, closed operation, and open-closing operation. Table 2 shows the fusion performance of different image denoising processing. As outlined in Table 2, the HFWT method without HS image denoising had the worst fusion results. The proposed methods with each HS image denoising preprocessing have the better fusion results compared with the proposed method without HS image denoising, which demonstrates that the HS image denoising preprocessing is significative and effective. By contrast, the proposed HFWT method with the open-closing operation achieves the best fusion performance, and it demonstrates that the open-closing operation is an effective HS image denoising preprocessing.  Figure 5 shows the fusion experimental results for the Washington DC dataset, where Figure 5(a1) shows the reference HS image, and the subjective fused HS images of each method are displayed in Figure 5(b1-g1). Moreover, Figure 5(a2) shows the enlarged subareas of the reference HS image, and the two enlarged subareas of each fused image are shown in Figure 5(b2-g2). The reference error image is shown in Figure 5(a3), and the error images between each fused HS image and the reference HS image are reported in Figure 5(b3-g3). Except for the first column, each column in Figure 5 shows the experimental results corresponding to each method. By visually comparing the fused HS images with the reference HS image, the fused result of the GS method suffers from serious spectral distortion. For example, the two enlarged subareas of the GS method are distorted seriously. The GFPCA approach generates fuzzy spatial details in some regions, such as the enlarged subareas shown in Figure 5(c2). This is because the spatial information of the GFPCA approach is injected insufficiently. As depicted in Figure 5(d1,d2), the spatial information of the fused image is well enhanced using the CNMF method, but some slight spectral distortion is appeared in the roof of the buildings. A closer inspection revealed that the HySure method seems to generate some distortion in the circular building in the upper left corner. By contrast, the fused HS images obtained by the Bayesian and HFWT methods achieve superior performance in terms of both spectral and spatial aspects. In order to further compare the performance of each fusion method, the third row of Figure 5 shows the error images of different methods. The error image is the difference (absolute value) of pixel values between the fused HS images and the reference HS image. We can see that, the GS, GFPCA and CNMF methods have larger differences, the HySure and Bayesian approaches generated relatively smaller differences, and the proposed HFWT approach shows the smallest differences in most areas, which demonstrates the excellent fusion capacity of the proposed method.

Experiments on Simulated Hyperspectral Datasets
aspects. In order to further compare the performance of each fusion method, the third row of Figure  5 shows the error images of different methods. The error image is the difference (absolute value) of pixel values between the fused HS images and the reference HS image. We can see that, the GS, GFPCA and CNMF methods have larger differences, the HySure and Bayesian approaches generated relatively smaller differences, and the proposed HFWT approach shows the smallest differences in most areas, which demonstrates the excellent fusion capacity of the proposed method.
Remote Sens. 2019, 11, 1005 12 of 18  Figure 6a4, and the SAM images of enlarged subarea obtained by each method are shown in Figures 6b4-g4. The spectral distortion caused by the GS method is very obvious, and the degree of spatial enhancement is also not acceptable for the GS method, as depicted on Figure 6b1 and b2. Compared with the GS approach, the GFPCA method performs better in terms of the spectral quality. However, the fused HS image obtained by the GFPCA approach shows an indistinct area in the left region of Figure 6c1. Despite having a preeminent spatial quality, the CNMF method generates significant spectral distortion in the triangle region in the lower half of Figure 6d1. From the visual analysis, the HySure, Bayesian and HFWT methods effectively improve spatial performance while maintaining spectral information, and the HFWT method shows better spectral quality compared with the HySure and Bayesian methods in some regions, such as the upper area of the enlarged subarea. The SAM images and the SAM images of the enlarged subareas of different approaches are shown in the third and fourth rows of Figure 6, to further verify the fusion performance of the proposed method. It can be seen that the proposed HFWT method yields the lowest SAM values for most regions. These results demonstrate that the proposed HFWT algorithm performs well in both the spatial and spectral aspects. Similar to the previous experiments, for the Salinas scene dataset, the fused results are shown in Figure 6. Figure 6(a1-a3) show the reference HS image, the enlarged subareas of the reference HS image, and the reference SAM image, respectively. Figure 6(b1-g1) in the first row show the subjective pansharpened results of each algorithm, and Figure 6(b2-g2) in the second row display each enlarged subarea. The SAM images of each approach are shown in Figure 6(b3-g3). The reference SAM image of enlarged subarea is shown in Figure 6(a4), and the SAM images of enlarged subarea obtained by each method are shown in Figure 6(b4-g4). The spectral distortion caused by the GS method is very obvious, and the degree of spatial enhancement is also not acceptable for the GS method, as depicted on Figure 6(b1,b2). Compared with the GS approach, the GFPCA method performs better in terms of the spectral quality. However, the fused HS image obtained by the GFPCA approach shows an indistinct area in the left region of Figure 6(c1). Despite having a preeminent spatial quality, the CNMF method generates significant spectral distortion in the triangle region in the lower half of Figure 6(d1). From the visual analysis, the HySure, Bayesian and HFWT methods effectively improve spatial performance while maintaining spectral information, and the HFWT method shows better spectral quality compared with the HySure and Bayesian methods in some regions, such as the upper area of the enlarged subarea. The SAM images and the SAM images of the enlarged subareas of different approaches are shown in the third and fourth rows of Figure 6, to further verify the fusion performance of the proposed method. It can be seen that the proposed HFWT method yields the lowest SAM values for most regions.
These results demonstrate that the proposed HFWT algorithm performs well in both the spatial and spectral aspects. quality, the CNMF method generates significant spectral distortion in the triangle region in the lower half of Figure 6d1. From the visual analysis, the HySure, Bayesian and HFWT methods effectively improve spatial performance while maintaining spectral information, and the HFWT method shows better spectral quality compared with the HySure and Bayesian methods in some regions, such as the upper area of the enlarged subarea. The SAM images and the SAM images of the enlarged subareas of different approaches are shown in the third and fourth rows of Figure 6, to further verify the fusion performance of the proposed method. It can be seen that the proposed HFWT method yields the lowest SAM values for most regions. These results demonstrate that the proposed HFWT algorithm performs well in both the spatial and spectral aspects.  In addition to visual inspection, the performance of each algorithm for the Washington DC and Salinas scene datasets is analyzed quantitatively in Table 3, where the best results for each quantitative index are marked in bold. As can be seen from Table 3, the objective quantitative results are roughly consistent with the subjective qualitative effects. Same as the subjective results, the GS and GFPCA algorithms produce worse objective performance compared with other algorithms. The HySure approach obtains the best RMSE value for the Washington DC dataset and the optimal ERGAS value for the Salinas scene dataset. Most of the quality indexes generated by applying the proposed HFWT method are the best, in which the SAM, CC and ERGAS values are the best for Washington DC, and the RMSE, SAM and CC indexes are ranked first for the Salinas scene. In addition to visual inspection, the performance of each algorithm for the Washington DC and Salinas scene datasets is analyzed quantitatively in Table 3, where the best results for each quantitative index are marked in bold. As can be seen from Table 3, the objective quantitative results are roughly consistent with the subjective qualitative effects. Same as the subjective results, the GS and GFPCA algorithms produce worse objective performance compared with other algorithms. The HySure approach obtains the best RMSE value for the Washington DC dataset and the optimal ERGAS value for the Salinas scene dataset. Most of the quality indexes generated by applying the proposed HFWT method are the best, in which the SAM, CC and ERGAS values are the best for Washington DC, and the RMSE, SAM and CC indexes are ranked first for the Salinas scene.  Figure 7 shows the pansharpened images of each method for the Hyperion dataset to confirm the fusion performance of the proposed HFWT method in the real dataset. Figure 7a-c show the real HS, real PAN, and interpolated HS images, respectively. The GS method shown in Figure 7d generates obvious spectral distortion, especially in the wharf area. In spite of good spatial improvement, the spatial details shown in the fused images obtained by using the GS and HySure approaches are too sharp. Spectral quality of the GFPCA and Bayesian methods seems be acceptable, but the GFPCA and Bayesian methods perform poorly from the spatial aspect. By contrast, the subjective effects of the CNMF and HFWT approaches are the best, and the HFWT method yields better spatial capacity compared to the CNMF method. The objective quality evaluation for the Hyperion dataset are presented in Table 4. As reported in Table 4, the HFWT method provides the best quantitative evaluation results in terms of the RMSE, SAM, CC, and ERGAS indices, which indicates that the HFWT method successfully maintains the spectral information of the original LRHS image and improves the spatial resolution.  Figure 7 shows the pansharpened images of each method for the Hyperion dataset to confirm the fusion performance of the proposed HFWT method in the real dataset. Figures 7a, b, and c show the real HS, real PAN, and interpolated HS images, respectively. The GS method shown in Figure 7d generates obvious spectral distortion, especially in the wharf area. In spite of good spatial improvement, the spatial details shown in the fused images obtained by using the GS and HySure

Computational Complexity Analysis and Time Comparisons
The proposed HFWT algorithm contains simple sequential statements, several loop statements without nesting, and a two-layer loop statements with nesting (In the proposed HFWT algorithm, the program statement of Equation (19) which are applied on each pixel is a two-layer loop statements). The simple sequential statement is naturally O(1) time, and the loop statement without nesting is O(n) time. The two-layer loop statement is O(n 2 ). According to the summation rule of algorithm complexity, the total algorithm complexity is O(n 2 + n + 1) = O(n 2 ). The proposed HFWT algorithm which is O(n 2 ) time belongs to the polynomial time, and is considered a fast algorithm. The computing time (in seconds) of each method for three datasets is shown in Table 5. The experiments in this paper were all performed using MATLAB R2015b, and tested on a PC with an Intel Core i5-7300HQ CPU @ 2.50 GHz and 8 GB of memory. The GS and GFPCA methods are very efficient, but the fusion performance of the GS and GFPCA methods is unsatisfactory. The proposed HFWT method is faster than the CNMF algorithm, and takes much less computing time than the HySure and Bayesian algorithms. The time cost of the proposed HFWT is acceptable.

Conclusions
This paper presents a novel hyperspectral pansharpening method based on the merger of the homomorphic filtering and weighted tensor matrix. The proposed HFWT algorithm introduces the open-closing morphological operation and homomorphic filtering to remove noise and extract spatial information of each band of an HhS image, respectively. Moreover, we propose a weighted RMSE-based method to obtain the total spatial information of the HS image. In order to generate the adequate spatial information from both the HS and the corresponding PAN images, an optimized weighted tensor matrix based method is proposed. Specifically, the weighted tensor matrix, eigenvalues and eigenvectors are deduced and analyzed to obtain the weighted gradient, and an optimization model is presented to acquire the integrated spatial information. Compared with the state-of-the-art methods, experiments performed on the Washington DC, Salinas scene and Hyperion datasets demonstrate the proposed method performs superiorly in terms of both subjective and objective assessment.