A Novel Adaptively Optimized PCNN Model for Hyperspectral Image Sharpening

: Hyperspectral satellite imagery has developed rapidly over the last decade because of its high spectral resolution and strong material recognition capability. Nonetheless, the spatial resolution of available hyperspectral imagery is inferior, severely affecting the accuracy of ground object identiﬁcation. In the paper, we propose an adaptively optimized pulse-coupled neural network (PCNN) model to sharpen the spatial resolution of the hyperspectral imagery to the scale of the multispectral imagery. Firstly, a SAM-CC strategy is designed to assign hyperspectral bands to the multispectral bands. Subsequently, an improved PCNN (IPCNN) is proposed, which considers the differences of the neighboring neurons. Furthermore, the Chameleon Swarm Optimization (CSA) optimization is adopted to generate the optimum fusion parameters for IPCNN. Hence, the injected spatial details are acquired in the irregular regions generated by the IPCNN. Extensive experiments are carried out to validate the superiority of the proposed model, which conﬁrms that our method can realize hyperspectral imagery with high spatial resolution, yielding the best spatial details and spectral information among the state-of-the-art approaches. Several ablation studies further corroborate the efﬁciency of our method.


Introduction
Hyperspectral (HS) images are acquired by sampling the spectrum range into a large number of spectral channels, rendering them as enhanced multispectral (MS) images with multiple bands, a narrow spectral range, and abundant spectral information.Since the abundance of spectral bands enables HS imagery to identify ground cover types accurately, HS imagery is broadly applied in environmental monitoring [1], agricultural product assessment [2], geology [3], and mineralogical mapping [4].However, because of the signalto-noise ratio (SNR) constraints of satellite sensors, the spatial resolution and the spectral resolution of HS imagery would have to be inevitably compromised in a unique acquisition.Hence, HS images are characterized by high spectral resolution but low spatial resolution, which limits their applications in higher precision remote sensing (RS) interpretation.RS image sharpening is a cost-effective way to generate imagery with simultaneously rich spectral information and high spatial resolution by injecting spatial details into either HS or MS imagery, which could enhance the spatial resolution of the RS images.RS image sharpening has been extensively studied over the last four decades.Nonetheless, most of these methods aim to fuse the MS imagery and the panchromatic (PAN) imagery, commonly known as pansharpening.As more HS satellites have been launched more recently (such as Earth Observing-1, ZY-1-02D, AVIRIS, and GF-5), HS sharpening has become increasingly prominent.
To solve the HS sharpening problem, many traditional pansharpening methods can still be applied.Pansharpening methods fall into two broad categories, namely the component substitution (CS) method and the multiresolution analysis (MRA) method [5].The CS method involves projecting the original MS image onto the transform domain and substituting the spatial component with the PAN image [6][7][8][9][10][11][12][13].Since the relatively large differences in spectral ranges between the PAN component and the replaced MS spatial component, the fusion result suffers from significant spectral distortion in the CS method.On the contrary, the MRA method preserves spectral information by injecting spatial details from the PAN imagery into the upsampled MS image, which is achieved by multiscale spatial filtering [14][15][16][17][18][19][20].The MRA method is able to conserve the spectral information well, whereas it cannot get impressive spatial details.
HS sharpening aims to produce HS imagery with high spatial resolution by fusing low-resolution HS imagery and high-resolution MS imagery.In order to tackle the HS sharpening issue, traditional CS and MRA pansharpening methods are typically utilized by simply replacing both the PAN imagery and the MS imagery with the MS imagery and the HS imagery simultaneously [21].Gomez et al. are the pioneers to apply the pansharpening approach to HS sharpening [22] and used the 2D-wavelet transform to render the MS and HS images fused in the same wavelength range.Chen et al. [23] present a generic HS sharpening framework based on the pansharpening method, which is the primary inspiration for later transmigration methods from pansharpening to HS sharpening.More recently, some methods for HS sharpening have also been proposed.Picone et al. addressed the band-assignment problem of the HS sharpening [24].Lu et al. propose a spectral modulation hyper-sharpening framework [25], which mitigates the problem of large spectral distortion after fusion.Yokoya et al. [26] propose a coupled nonnegative matrix factorization method (CNMF), which generates the end members and abundance matrices by alternately unmixing HS and MS images with the NMF algorithm.However, the three-dimensional structure of the HS image is hardly reserved by the matrix factorization.Thus, tensor factorization is utilized in place of the matrix factorization [21].Dian et al. present a nonlocal sparse tensor factorization method (NLSTF_SMBF) for the semi-blind fusion of HS and MS images [27], which primarily constructs full-band patches (FBPs), and similar FBPs can share the same dictionary.
In recent years, deep learning (DL) methods have been extensively applied in the field of HS sharpening [28,29].Zhang et al. [30] proposed an unsupervised deep learning network architecture for the simultaneous optimization of the HS super-resolution and the degradation estimation.Qu et al. [31] proposed an unsupervised, unaligned Mutual HS super-resolution Dirichlet-Net, which effectively improves the robustness in the face of alignment errors.
Although promising fusion results are realized by various HS sharpening methods, it is expected to further improve the HS fusion accuracy by neural network algorithms.Nonetheless, due to the independent spectral features of the unique RS image, it is difficult to establish a universal training database for the RS fusion applications.Besides, it is time-consuming to train the network.A pulse-coupled neural network (PCNN) is a kind of biologically inspired neural network without training, which exhibits the characteristics of the pulse-synchronous phenomenon.The property of pulse-synchronous allows the synchronously stimulated pixels to generate the segmentation results, which turns out to be in accordance with the human visual mechanism [32].In order to improve the accuracy of image fusion, it is useful to adopt different fusion strategies for different segmentation regions that are consistent with human perception.Traditional PCNN fusion models include simplified PCNN (SPCNN) [33], dual-channel PCNN (DCPCNN) [34], and shuffled frog-leaping PCNN [35].Particularly, Panigrahy et al. propose adaptive DCPCNN for multi-aggregation and medical image fusion [36,37].However, traditional PCNN fusion approaches are applicable only to either the medical image fusion or the multi-focus image fusion, which cannot be applied to the RS image fusion.Recently, Li et al. presented a modified PCNN model for RS image fusion [32], but it only works on MS imagery.
In this paper, a novel adaptively optimized PCNN model for HS sharpening is proposed.Concretely, we summarize the main contributions of the paper as follows: (1) a SAM-CC band assignment method is proposed to group the HS bands with the MS band; (2) an improved PCNN (IPCNN) model is proposed to obtain irregular injection region of spatial details. (3) A Chameleon Swarm Optimization (CSA)-based IPCNN parameter optimization method is designed to achieve the optimal sharpening imagery.The comparative experiments were carried out on three datasets captured by the ZY-1-02D and GF-2 satellites, and the results substantiate the effectiveness of the proposed method.
The rest of the article is organized as follows.The description of the PCNN and CSA principles is given in Section 2. Section 3 presents the details of our proposed fusion approach.Experimental results and discussions are provided in Section 4. Section 5 contains the conclusion of the paper.

Standard PCNN Principle
PCNN belongs to the third generation of artificial neural networks, which is a kind of neural network model proposed by Johnson [38] on the basis of the observation of pulse delivery experiments in the cerebral cortex of cats and monkeys.The individual neuron in the PCNN model is partitioned into the receptive field, the modulation field, and the pulse delivery field according to its functions.Mathematical formulae of the receptive field are given in Equations ( 1) and (2).The mathematical description of the modulation field can be found in Equation (3), and the mathematical formulae of the pulse delivery field can be defined in Equations ( 4) and ( 5).
where the subscripts ij and kl refer to the positions of the current neuron and the neighboring neurons, respectively.n denotes the nth iteration.I is the input image, whereas I ij corresponds to the external input of the neuron ij.F and L indicate the feeding input and the linking input, respectively, while the difference between them is that the linking input L can only receive the local pulse stimulation from the neighboring neuron via the synaptic matrix W. On the contrary, F receives not only the local stimulation but also the external stimulation from I. β denotes the linking strength, which ranges from 0 to 1.A larger value of β indicates that the linking relationships are greater, i.e., the current neuron is more susceptible to neighborhood neurons.M and W represent the connection coefficient matrixes, which are typically calculated by the Euclidean distance between the current neuron and the neighborhood.U refers to the internal activity, and E represents the dynamic threshold.
U and E jointly determine whether or not the current neuron is stimulated in the current iteration, i.e., the boolean variable Y ij becomes equal to one whenever U ij is greater than E ij .α F , α L , and α E are time decay constants, which are utilized to adjust the decay rate of F, L and E. V F , V L and V E indicate the normalization constants.The matrix Y[n] denotes the binary output of PCNN in each iteration n, whereas Y ij is 0 (when the neuron ij is unstimulated) and 1 (when stimulated), respectively.In addition, as illustrated in Figure 1, each PCNN neuron corresponds to a pixel in the image when applying PCNN to image processing applications.Therefore, the current neuron ij will release a pulse as long as U ij is greater than E ij .If the current pixel is stimulated in the current iteration, E will automatically spike the value V E ; otherwise, E will decrease gradually as the iteration number n increases.After the nth iteration of PCNN, the neurons with either similar grayscale values or the adjacent positions will be synchronously stimulated, while the output Y will eventually form a binary ignition map.denotes the binary output of PCNN in each iteration n, whereas Yij is 0 (when the neuron ij is unstimulated) and 1 (when stimulated), respectively.In addition, as illustrated in Figure 1, each PCNN neuron corresponds to a pixel in the image when applying PCNN to image processing applications.Therefore, the current neuron ij will release a pulse as long as Uij is greater than Eij.If the current pixel is stimulated in the current iteration, E will automatically spike the value VE; otherwise, E will decrease gradually as the iteration number n increases.After the nth iteration of PCNN, the neurons with either similar grayscale values or the adjacent positions will be synchronously stimulated, while the output Y will eventually form a binary ignition map.

Chameleon Swarm Optimization Algorithm
The CSA algorithm is a meta-heuristic method proposed by Malik [39].The model imitates the socially intelligent synergistic behavior of chameleons when foraging and capturing food near woods, swamps, and deserts.It is a bio-optimization algorithm for finding the global optimum of nonlinear, nonconvex, and other complex problems, which may prevent entrapment in the local optimum.
The principle of CSA is shown in Figure 2. The algorithm mathematizes the behavioral stages of a chameleon in looking for food, which include initializing the starting position, tracking the prey from a distance, locating the prey by eye rotation, and catching the prey with a high-velocity sticky tongue.

Chameleon Swarm Optimization Algorithm
The CSA algorithm is a meta-heuristic method proposed by Malik [39].The model imitates the socially intelligent synergistic behavior of chameleons when foraging and capturing food near woods, swamps, and deserts.It is a bio-optimization algorithm for finding the global optimum of nonlinear, nonconvex, and other complex problems, which may prevent entrapment in the local optimum.
The principle of CSA is shown in Figure 2. The algorithm mathematizes the behavioral stages of a chameleon in looking for food, which include initializing the starting position, tracking the prey from a distance, locating the prey by eye rotation, and catching the prey with a high-velocity sticky tongue.

Start
Set parameters and initialize the population.
Set fitness function, calculated fitness values of the initial chameleon group according to f, and record the optimal value and its position.
Execute the part of r i ≥P p , and update the position.
Update the position of eyes after rotation.
Update the position of the attacked prey.

Whether r i ≥P p
Execute the part of r i ＜P p , and update the position.
Calculated fitness values of the current chameleon group according to f, and record the optimal value and its position.
Output global optimal value and its position.[39] for all variable definitions.Figure 2. Flowchart of CSA.See reference [39] for all variable definitions.

Proposed Method
In order to tackle the problem of spectral distortion and fuzzy spatial details in HS sharpening, we propose a novel adaptively optimized PCNN fusion algorithm.The flowchart of the proposed HS sharpening algorithm is shown in Figure 3, which contains the following modules:

SAM-CC Band Assignment Block
In contrast to the multi-to-single band assignment in the pansharpening problem, HS sharpening is a multiband-to-multiband fusion.Thus, as shown in Figure 4, correctly assigning each HS band to the MS band is an essential procedure before the actual fusion operation, i.e., the choice of which band in the MS imagery can be used to sharpen the corresponding HS bands plays a crucial role.Classical band selection methods include the minimum spectral distance (MSD) assignment algorithm, the maximum cross-correlation (CC) assignment algorithm, and the minimum spectral distortion (SAM) assignment algorithm [24].Since the optimal band selection is challenging for any single criterion, we propose a joint band selection algorithm using both the SAM and CC indices (SAM-CC) in order to group the HS bands more accurately.

SAM-CC Band Assignment Block
In contrast to the multi-to-single band assignment in the pansharpening problem, HS sharpening is a multiband-to-multiband fusion.Thus, as shown in Figure 4, correctly assigning each HS band to the MS band is an essential procedure before the actual fusion operation, i.e., the choice of which band in the MS imagery can be used to sharpen the corresponding HS bands plays a crucial role.Classical band selection methods include the minimum spectral distance (MSD) assignment algorithm, the maximum cross-correlation (CC) assignment algorithm, and the minimum spectral distortion (SAM) assignment algorithm [24].Since the optimal band selection is challenging for any single criterion, we propose a joint band selection algorithm using both the SAM and CC indices (SAM-CC) in order to group the HS bands more accurately.
Let H = {H} h = 1,. .., Nh , stands for the HS images and M = {M} m = 1,. .., Nm , represents the MS images, where h and m are the band numbers of HS and MS, respectively.N h and N m refer to the total number of bands of HS and MS.The proposed SAM-CC criterion is represented as Equation ( 6).

SAM-CC(H
where E is a unit matrix with size N M × N H .The symbol <•> indicates the inner product operation.The symbol • stands for the l 2 norm.r refers to the spatial resolution ratio of HS to MS. M r j represents the low-pass filtered down-sampled image of M j .
Remote Sens. 2023, 15, 4205 7 of 23 where E is a unit matrix with size NM × NH.The symbol <•> indicates the inner product operation.The symbol ‖ • ‖ stands for the l2 norm.r refers to the spatial resolution ratio of HS to MS. Mj r represents the low-pass filtered down-sampled image of Mj.

Improved PCNN Model
To make the PCNN model more suitable for HS sharpening, we propose an IPCNN model in the paper.In the standard PCNN model, the feeding input F stands for the accumulative influence of the external stimuli.Nonetheless, the human eyes are more sensitive to the edges and orientations rather than individual pixels within an image.Therefore, a new feeding input SF is designed in the IPCNN model, which considers the local neighborhood differences in both the horizontal and vertical directions.The SF is calculated in Equations ( 9)- (11).

Improved PCNN Model
To make the PCNN model more suitable for HS sharpening, we propose an IPCNN model in the paper.In the standard PCNN model, the feeding input F stands for the accumulative influence of the external stimuli.Nonetheless, the human eyes are more sensitive to the edges and orientations rather than individual pixels within an image.Therefore, a new feeding input SF is designed in the IPCNN model, which considers the local neighborhood differences in both the horizontal and vertical directions.The SF is calculated in Equations ( 9)- (11).
where RF ij denotes the row difference of the two neighboring feeding inputs F via the local rectangular window M 1 × M 2 , and CF ij indicates the column difference.After performing a number of experiments, we take the window size to be 5 × 5. Furthermore, I ij is designed to strengthen the impact of the local neighborhood, which is given by Equation (12).
where H ij refers to the pixel value of the upsampled HS in row i and column j.
After modeling the new SF and I , the proposed IPCNN neuron model is illustrated in Figure 5.The IPCNN model has the following advantages compared to the standard PCNN: (1) the IPCNN model utilizes the spatial difference of the neighborhoods to stimulate the feeding input F ij of the neuron, which can describe the local detail features; (2) the external stimuli I ij considers the influence of the surrounding pixels towards the central pixels.
where RFij denotes the row difference of the two neighboring feeding inputs F via the local rectangular window M1 × M2, and CFij indicates the column difference.After performing a number of experiments, we take the window size to be 5 × 5. Furthermore, Iij′ is designed to strengthen the impact of the local neighborhood, which is given by Equation ( 12).
where Hij refers to the pixel value of the upsampled HS in row i and column j.
After modeling the new SF and I′, the proposed IPCNN neuron model is illustrated in Figure 5.The IPCNN model has the following advantages compared to the standard PCNN: (1) the IPCNN model utilizes the spatial difference of the neighborhoods to stimulate the feeding input Fij of the neuron, which can describe the local detail features; (2) the external stimuli Iij considers the influence of the surrounding pixels towards the central pixels.

Automatic Parameter Optimization of IPCNN by CSA
Similar to the standard PCNN, the setting of appropriate parameters (i.e., αF, αL, αE, β) of the IPCNN model is most critical.As shown in Figure 6, different combinations of IPCNN parameters will lead to different segmentation results for the same input image.If there are too many segmentation pieces, the segmentation regions will be too small and will further affect the accuracy and complexity of the statistical computation of subsequent injection weights.Conversely, if the number of image segmentation pieces is too few, it is impossible to distinguish each region or to take advantage of the characteristics of each region.So, different parameters can lead to different segmentation results, and different segmentations will have corresponding impacts on the final fusion results.Most researchers choose to simplify the PCNN model or to use the manually set uniform parameters for the different images.However, fixed parameters for different input images do not lead to the optimal fusion results for all images.Thus, we propose a CSA-based optimization approach for setting IPCNN parameters αF, αL, αE, β, and W, which can adaptively generate their own optimal parameters for various input images.

Automatic Parameter Optimization of IPCNN by CSA
Similar to the standard PCNN, the setting of appropriate parameters (i.e., α F , α L , α E , β) of the IPCNN model is most critical.As shown in Figure 6, different combinations of IPCNN parameters will lead to different segmentation results for the same input image.If there are too many segmentation pieces, the segmentation regions will be too small and will further affect the accuracy and complexity of the statistical computation of subsequent injection weights.Conversely, if the number of image segmentation pieces is too few, it is impossible to distinguish each region or to take advantage of the characteristics of each region.So, different parameters can lead to different segmentation results, and different segmentations will have corresponding impacts on the final fusion results.Most researchers choose to simplify the PCNN model or to use the manually set uniform parameters for the different images.However, fixed parameters for different input images do not lead to the optimal fusion results for all images.Thus, we propose a CSA-based optimization approach for setting IPCNN parameters α F , α L , α E , β, and W, which can adaptively generate their own optimal parameters for various input images.
The IPCNN parameters need to be jointly optimized together, so the CSA is employed to automatically optimize all five IPCNN parameters.For convenience, the connection weight W is denoted as W = [w, 1, w; 1, 0, 1; w, 1, w], where w is optimized rather than the entire matrix W. The flowchart of the automatic IPCNN parameter optimization algorithm based on CSA is shown in Figure 7. Firstly, the chameleon position is initialized with the classical parameter values.Secondly, the fitness function is set as the weighted summation of spectral fidelity SAM [5] and spatial detail representation ERGAS [40].The fitness function for the proposed optimization method is given as follows.
where FH denotes the reference imagery, and F refers to the fusion imagery.The IPCNN parameters need to be jointly optimized together, so the CSA is employed to automatically optimize all five IPCNN parameters.For convenience, the connection weight W is denoted as W = [w, 1, w; 1, 0, 1; w, 1, w], where w is optimized rather than the entire matrix W. The flowchart of the automatic IPCNN parameter optimization algorithm based on CSA is shown in Figure 7. Firstly, the chameleon position is initialized with the classical parameter values.Secondly, the fitness function is set as the weighted summation of spectral fidelity SAM [5] and spatial detail representation ERGAS [40].The fitness function for the proposed optimization method is given as follows.The IPCNN parameters need to be jointly optimized together, so the CSA is employed to automatically optimize all five IPCNN parameters.For convenience, the connection weight W is denoted as W = [w, 1, w; 1, 0, 1; w, 1, w], where w is optimized rather than the entire matrix W. The flowchart of the automatic IPCNN parameter optimization algorithm based on CSA is shown in Figure 7. Firstly, the chameleon position is initialized with the classical parameter values.Secondly, the fitness function is set as the weighted summation of spectral fidelity SAM [5] and spatial detail representation ERGAS [40].The fitness function for the proposed optimization method is given as follows.
( ) , , arccos FH F where FH denotes the reference imagery, and F refers to the fusion imagery.SAMr and ERGASr stand for the ranges of SAM and ERGAS.h and l indicate the spatial resolutions of the MS and HS imagery, respectively.uh indicates the mean value of the hth band of the reference imagery.<⸳> means the inner product.‖⸳‖2 refers to the second-order norm operation.Since there are too many bands in the HS image, directly optimizing the parameters for all bands will be time-consuming.Now that the adjacent HS bands acquired by the SAM-CC band assignment method are usually strongly correlated, it would be more efficient to obtain optimized parameters by only using the representative bands of each group.Wang et al. proposed a HS band selection method with optimal domain reconstruction [41], which can adaptively achieve the optimal neighborhood reconstruction (ONR) to find the subsets of bands that can best represent the HS image.Thus, the ONR band selection method is used before IPCNN parameter optimization.In addition, the method can also reduce the impacts of noisy bands by taking advantage of the characteristics of neighborhood bands.
To be more efficient, the optimal IPCNN parameters are generated with the representative bands of HS images.Other HS bands could utilize these parameters as long as both bands are very similar to each other.Representative bands from each group are chosen by the ONR method.If there are more than five HS bands in each group, five of these are then selected using ONR to compose the optimal subset of bands.Otherwise, all of the bands from the original HS within the group are retained.Equation ( 16) is defined to express the procedure.

Extracting MS Details
The multiscale analysis has demonstrated excellent capability in RS image fusion.The "atrous" wavelet transform is an undecimated multiscale analysis method, which is widely used to address image fusion issues due to its fast decomposition and reconstruction as well as the richer extracted details.For each histogram-matched MS band M hm (matched by HS image), the "atrous" wavelet decomposition is performed to extract the spatial details from the multispectral imagery.After wavelet decomposition, the low-frequency component stands for the approximate image, and the high-frequency component indicates the noise and the local features.After setting the high-frequency components to zero, the low-frequency MS imagery M L is then realized via wavelet reconstruction.Finally, the spatial detail information M d is obtained through Equation (17).

Adaptive Injected Gains
The HS imagery is divided into different irregular segmentation regions by the proposed CSA-based IPCNN segmentation algorithm.The injected gains are, therefore, calculated in each irregular region.The overall steps of the adaptive gain approach are described as follows: where (i, j) stands for the pixel coordinates of the imagery.H k u denotes the upsampled image of HS, and M L indicates the low-resolution version of the MS imagery.Y[n] refers to the activated neurons in the current iteration.R is the correlation ratio between the HS and MS imagery.corr(A, B) represents the correlation coefficient between matrix A and matrix B. cov(A, B) denotes the covariance between matrix A and matrix B. std(A) refers to the standard deviation of A.

Fusion Output
The final fusion result can be calculated by Equation (21).
where FH k denotes the high-resolution HS fusion image.
The fusion pseudo-code of the proposed HS sharpening algorithm (Algorithm 1) is described as follows.
Algorithm 1 AT-AIPCNN ("atrous" transform-adaptive IPCNN) method matrix B. cov(A, B) denotes the covariance between matrix A and matrix B. std(A) refers to the standard deviation of A.

Fusion Output
The final fusion result can be calculated by Equation (21).
where FHk denotes the high-resolution HS fusion image.The fusion pseudo-code of the proposed HS sharpening algorithm (Algorithm 1) is described as follows.

Datasets
Three real datasets captured from the ZY1-02D and GF-2 satellites are utilized to test the validity of the proposed method.The HS sensor AHSI of the ZY1-02D satellite, which works in the spectral range from the visible to the short-wave infrared (with the wavelength from 395 nm to 2500 nm), provides 166 HS bands with a spatial resolution of 30 m.The GF-2 satellite captures four multispectral bands with a spatial resolution of 3.4 m, spanning the spectral range from the visible to the near-infrared (450-890 nm).For HS sharpening application, HS images from the ZY1-02D sensor would need to be sharpened

Datasets
Three real datasets captured from the ZY1-02D and GF-2 satellites are utilized to test the validity of the proposed method.The HS sensor AHSI of the ZY1-02D satellite, which works in the spectral range from the visible to the short-wave infrared (with the wavelength from 395 nm to 2500 nm), provides 166 HS bands with a spatial resolution of 30 m.The GF-2 satellite captures four multispectral bands with a spatial resolution of 3.4 m, spanning the spectral range from the visible to the near-infrared (450-890 nm).For HS sharpening application, HS images from the ZY1-02D sensor would need to be sharpened to the spatial resolution of MS imagery of the GF-2 sensor.Thus, dataset1 consists of imagery taken over the Liujiaxia reservoir in China from the ZY1-02D and GF-2 sensors, which mainly include the lake and the village.Dataset2 stands for the suburban area of Linxia City in China, mainly containing mountainous terrain.The dataset3 is the farmland area of Yongchang City in China, which primarily includes farmland, mountains, and buildings.The original images of the three datasets are shown in Figure 8.

Experimental Setup
In the experiments, most program code is executed with matlab2020a on the Intel CPU Core i7-13700K and NVIDIA GPU GeForce GTX 3090.The number of chameleon populations is set to 20 in the CSA-based IPCNN adaptive optimization, and the maximum number of iterations is set to 30.
Nine classical and competitive methods of different sharpening categories are compared with the proposed approach, i.e., the Gram-Schmidt adaptive (GSA) algorithm [12], the smoothing filtered-based intensity modulation (SFIM) algorithm [42], the generalized Laplacian pyramid (GLP) algorithm [43], the CNMF algorithm [26], the nonlocal sparse tensor factorization (NLSTF) algorithm [27], the NLSTF_SMBF algorithm [27], the HYSURE algorithm [44], the fast fusion based on Sylvester equation (FUSE) algorithm [45], and the UDALN algorithm [46].Among which, GSA, SFIM, and GLP are owned by the classical pansharpening methods, GSA belongs to the CS methods, SFIM and GLP are the MRA method, CNMF belongs to the matrix decomposition methods, NLSTF and NLSTF_SMBF are owned by tensor decomposition categories, while the HYSURE and FUSE belong to Bayesian-based methods.In addition, GSA, SFIM, and GLP do not require any parameter setting.The main parameters for the NLSTF and NLSTF_SMBF methods are set according to Table 1, and the parameters of other comparative approaches are set to be the same as in the literature [47].The UDALN is implemented in PyTorch, and the parameter settings are the same as in the literature [46].to the spatial resolution of MS imagery of the GF-2 sensor.Thus, dataset1 consists of imagery taken over the Liujiaxia reservoir in China from the ZY1-02D and GF-2 sensors, which mainly include the lake and the village.Dataset2 stands for the suburban area of Linxia City in China, mainly containing mountainous terrain.The dataset3 is the farmland area of Yongchang City in China, which primarily includes farmland, mountains, and buildings.The original images of the three datasets are shown in Figure 8.

Experimental Setup
In the experiments, most program code is executed with matlab2020a on the Intel CPU Core i7-13700K and NVIDIA GPU GeForce GTX 3090.The number of chameleon populations is set to 20 in the CSA-based IPCNN adaptive optimization, and the maximum number of iterations is set to 30.
Nine classical and competitive methods of different sharpening categories are compared with the proposed approach, i.e., the Gram-Schmidt adaptive (GSA) algorithm [12], the smoothing filtered-based intensity modulation (SFIM) algorithm [42], the generalized Laplacian pyramid (GLP) algorithm [43], the CNMF algorithm [26], the nonlocal sparse tensor factorization (NLSTF) algorithm [27], the NLSTF_SMBF algorithm [27], the HYS-URE algorithm [44], the fast fusion based on Sylvester equation (FUSE) algorithm [45], and the UDALN algorithm [46].Among which, GSA, SFIM, and GLP are owned by the classical pansharpening methods, GSA belongs to the CS methods, SFIM and GLP are the MRA method, CNMF belongs to the matrix decomposition methods, NLSTF and NLSTF_SMBF are owned by tensor decomposition categories, while the HYSURE and FUSE belong to Bayesian-based methods.In addition, GSA, SFIM, and GLP do not require any parameter setting.The main parameters for the NLSTF and NLSTF_SMBF methods are set according to Table 1, and the parameters of other comparative approaches are set to be the same as in the literature [47].The UDALN is implemented in PyTorch, and the parameter settings are the same as in the literature [46].Cluster scaling parameter: K = 151.The spectral response matrix R was estimated by HYSURE [44].
The performance of the proposed HS sharpening approach is evaluated by eight complementary quantitative indices to verify the spectral and spatial qualities of the fusion imagery, such as the peak signal-to-noise ratio (PSNR) [47], the root mean square error (RMSE) [48], the error relative global adimensionnelle de synthèse (ERGAS) [40], the spectral angle mapper (SAM) [5], the structural similarity index measurement (SSIM) [49,50], the universal image quality index (UIQI) [48], the inter-correlation (CC) [51], and the degree of distortion (DD) [52].The ideal values for RMSE, ERGAS, SAM, and DD are 0, whereas 1 for SSIM, UIQI, and CC.Moreover, the ideal value for PSNR is positive infinity.
PSNR measures the spatial similarity between the fusion imagery and the reference imagery.In general, the larger the PSNR, the more spatially similar the fusion imagery is to the reference imagery.The PSNR is defined as: where FH denotes the reference imagery, and F refers to the fusion imagery.h stands for the band number, while N h indicates the total number of all bands.log10(•) represents the logarithm function with base 10.FH h denotes the hth band of the reference imagery.max(•) denotes the maximum value.L denotes the total number of pixels per band.FH hj denotes the jth pixel of the hth band of the reference imagery, and F hj stands for the jth pixel of the hth band of the fusion image.RMSE refers to the measure of deviation between two images.The RMSE is defined as: where T stands for the total pixel number of the reference imagery.ERGAS is a comprehensive index that reflects both spectral distortions and spatial detail differences.The ERGAS can be computed through Equation (14).SAM calculates the similarity of the spectral vectors.The smaller one indicates that the spectrum is better maintained.The SAM index of two images FH and F is computed as Equation (15).UIQI measures the similarity in terms of the brightness and contrast.For the window, a × b, Q is defined as: where µ x and σ 2 x denote the expectation and variance of x, respectively, and σ 2 x,y represents the covariance of x and y.
Then, the UIQI can be computed through the following formula: SSIM compares the structural similarities between two images, which indicate the luminance distortion, the contrast distortion, and the structural distortion.Its formula is described in Equations ( 26) and (27).
where C 1 and C 2 are infinitesimal small constants to ensure stability.DD is an indicator to verify the quality of the fusion spectra.The smaller the value, the better the spectral retention.If the DD is equal to 0, there is no spectral distortion.Its calculation formula is shown in Equation (28).
where T denotes the total pixel number of the reference imagery.vec(x) represents matrix x as a vector, and • 1 is the one-order norm operation.CC represents the spatial correlation between two images, which measures the geometric distortion of the fusion imagery.It is calculated as follows: where X h refers to the hth band of X, and µ x indicates the mean value of x.

Experimental Results
The first experiment is conducted on the Liujiaxia dataset.The fusion results are shown in Figure 9, in which the local area is enlarged for convenient observation.As can be seen from Figure 9, the color of the GSA, SFIM, GLP, CNMF, HYSURE, UDALN, and FUSE methods appear to be dark.Our proposed method, NLSTF and NLSTF_SMBF methods, perform best in spectral preservation; nevertheless, the details of both the NLSTF and NLSTF_SMBF methods are blurred into chunks.Figure 10 presents the SAM error maps of the fusion images in all 166 bands with the Liujiaxia dataset.It is noted that the spectral distortion for all methods primarily occurs at the reservoir-mountain boundary, of which NLSTF and NLSTF_SMBF exhibit the largest aberration.Besides, UDALN has larger distortions in the reservoir area, and the proposed method, SFIM, and GLP perform better in the SAM error maps.Furthermore, quantitative assessments with the Liujiaxia dataset are shown in Table 2.It is clear from Table 2 that the proposed method obtains the best performance across all quantitative metrics, which demonstrates the superiority of the proposed method with respect to spectral preservation and spatial details.The fusion results with the Linxia dataset are shown in Figure 11.As can be seen, the large color difference and block blur are presented in the NLSTF and NLSTF_SMBF fusion methods.Moreover, the SFIM fusion image appears dark.GSA algorithm, as a kind of CS method, has greater spectral distortion, while MRA methods, such as SFIM and GLP, have less spatial detail.Furthermore, the proposed method, GSA, and UDALN, have fine spatial detail information.In addition, the spectra of the proposed method are the closest to the ground-truth image in subjective visualization.Figure 12 shows the SAM error maps of all the fusion results with the Linxia dataset.As can be seen from Figure 12, spectral distortion is more likely to occur in the ridge area and in the land-lake boundary.Although the SFIM method has less spectral distortion in most regions, it exhibits larger errors in some local areas (top left corner and bottom left of Figure 12d).NLSTF and NLSTF_SMBF appear to have more spectral distortion, whereas the proposed method performs well.Table 3 gives the quantitative evaluation results of the Linxia dataset, which demonstrates that our method achieves the best indicator results.
the ground-truth image in subjective visualization.Figure 12 shows the SAM error maps of all the fusion results with the Linxia dataset.As can be seen from Figure 12, spectral distortion is more likely to occur in the ridge area and in the land-lake boundary.Although the SFIM method has less spectral distortion in most regions, it exhibits larger errors in some local areas (top left corner and bottom left of Figure 12d).NLSTF and NLSTF_SMBF appear to have more spectral distortion, whereas the proposed method performs well.Table 3 gives the quantitative evaluation results of the Linxia dataset, which demonstrates that our method achieves the best indicator results.The best result is in bold, and the second best is underlined.
Figure 13 presents the fusion results of the Yongchang dataset.From Figure 13, we can see that the color error occurred in the fusion result of the CNMF, NLSTF, and NLSTF_SMBF methods.In addition, The SAM error maps of the Yongchang dataset are shown in Figure 14.We note from Figure 14 that the proposed method has the least amount of spectral distortion at the edges, which is also verified by Table 4.The best result is in bold, and the second best is underlined.
Figure 13 presents the fusion results of the Yongchang dataset.From Figure 13, we can see that the color error occurred in the fusion result of the CNMF, NLSTF, and NLSTF_SMBF methods.In addition, The SAM error maps of the Yongchang dataset are shown in Figure 14.We note from Figure 14 that the proposed method has the least amount of spectral distortion at the edges, which is also verified by Table 4.

Ablation Experiments
Three kinds of ablation experiments are carried out to test our three primary fusion modules, i.e., SAM-CC band assignment module, ONR band selection module, and CSA adaptive PCNN parameter module.In addition, different automatic parameter optimization strategies are also investigated.
The fusion results obtained on three different datasets with different band assignment strategies (SAM-CC, SAM, and CC) are presented in Table 5.For the Liujiaxia dataset, our proposed SAM-CC assignment indicator performs best.Furthermore, the SAM-CC achieves better fusion results in most cases with the Linxia dataset and the Yongchang dataset, which illustrates the effectiveness of the SAM-CC band assignment.To validate the efficiency of the ONR band selection module, we compared the fusion results of using ONR and without ONR.In order to verify the effectiveness of ONR band selection, we fused the data sets with and without ONR band selection.Comparisons of the fusion accuracy and the running time are shown in Table 6.As can be seen from Table 6, the running time has been reduced by more than 84% with the ONR band selection module, while the fusion accuracy remains almost unchanged.Overall, this indicates significant time cost savings for the ONR band selection module.The best result is in bold, and the second best is underlined.
In order to examine the rationality of the CSA-based parameter optimization strategy in IPCNN, Table 7 shows the comparison of the CSA fusion metrics with other three parameter optimization strategies, i.e., the sparrow search strategy (SSA) [53], the improved grey wolf optimizer (IGWO) strategy [54] and the enhanced whale optimization strategy (EWOA) [55].As shown in Table 7, the CSA-based fusion algorithm achieves a better fusion performance as well as less time-consuming.The best result is in bold, and the second best is underlined.
Table 8 exhibits the quantitative analysis of the fusion results with and without adaptive CSA optimization to verify the impact of the CSA algorithm.Besides, the traditional IPCNN parameters are set as classical values, i.e., α F = 0.1, α L = 1, α E =0.62, β = 0.1, w = 0.5.In all three datasets, the adaptive IPCNN approach is significantly superior to the traditional parameters in most cases, which indicates that the adaptive parameter optimization effectively improves the quality of the fusion images.The best result is in bold, and the second best is underlined.

Conclusions
In the paper, a novel algorithm for HS sharpening is proposed.Firstly, the band groups are acquired by performing joint SAM-CC band matching of HS and MS assignment method, the band grouping results are obtained, and then ONR is applied to the HS images of each group to improve the efficiency for subsequent adaption.Besides, according to the characteristics of remote sensing image fusion application, an IPCNN model is proposed, which can obtain irregular injection regions of spatial details.In addition, a CSA-based IPCNN parameter optimization method is designed to achieve the optimal sharpening imagery.In summary, the proposed method is simple and easy to implement, which presents good fusion results on a variety of datasets, including reservoir, mountain, town, and river landscapes.Furthermore, several ablation experiments are also conducted to corroborate the efficiency of the proposed method.
(1) the SAM-CC band assignment part; (2) the improved PCNN (IPCNN) model part; (3) the automatic parameter optimization of IPCNN by CSA part; (4) the extracting MS detail part; (5) the adaptive injected gains part; (6) the fusion output part.Remote Sens. 2023, 15, 4205 6 of 23 3. Proposed Method In order to tackle the problem of spectral distortion and fuzzy spatial details in HS sharpening, we propose a novel adaptively optimized PCNN fusion algorithm.The flowchart of the proposed HS sharpening algorithm is shown in Figure 3, which contains the following modules: (1) the SAM-CC band assignment part; (2) the improved PCNN (IPCNN) model part; (3) the automatic parameter optimization of IPCNN by CSA part; (4) the extracting MS detail part; (5) the adaptive injected gains part; (6) the fusion output part.

Figure 3 .
Figure 3. Flowchart of the proposed fusion algorithm.

Figure 3 .
Figure 3. Flowchart of the proposed fusion algorithm.

Figure 4 .
Figure 4. Schematic diagram of HS and MS bands grouping.Let H = {H}h = 1,…, Nh, stands for the HS images and M = {M}m = 1,…, Nm, represents the MS images, where h and m are the band numbers of HS and MS, respectively.Nh and Nm refer to the total number of bands of HS and MS.The proposed SAM-CC criterion is represented as Equation (6).

Figure 4 .
Figure 4. Schematic diagram of HS and MS bands grouping.
) Optimize IPCNN parameters α F , α L , α E , β, and W using the CSA-based IPCNN optimization algorithm.(3) Obtain the irregular segmentation region of the IPCNN model in the current iteration n.And calculate the injected gain G k [n] according to Equation (20).

Figure 8 .
Figure 8.The original images of three datasets.(a) HS image of Liujiaxia dataset.(b) HS image of Linxia dataset.(c) HS image of Yongchang dataset.(d) MS image of Liujiaxia dataset.(e) MS image of Linxia dataset.(f) MS image of Yongchang dataset.

Figure 8 .
Figure 8.The original images of three datasets.(a) HS image of Liujiaxia dataset.(b) HS image of Linxia dataset.(c) HS image of Yongchang dataset.(d) MS image of Liujiaxia dataset.(e) MS image of Linxia dataset.(f) MS image of Yongchang dataset.

Table 2 .
Quantitative evaluation results of the Liujiaxia dataset.

Table 2 .
Quantitative evaluation results of the Liujiaxia dataset.

Table 3 .
Quantitative evaluation results of the Linxia dataset.

Table 3 .
Quantitative evaluation results of the Linxia dataset.

Table 4 .
Quantitative evaluation results of the Yongchang farmland area data set.

Table 4 .
Quantitative evaluation results of the Yongchang farmland area data set.

Table 4 .
Quantitative evaluation results of the Yongchang farmland area data set.

Table 5 .
SAM-CC ablation.The best result is in bold, and the second best is underlined.

Table 6 .
ONR band selection for ablation.