Source-linked AI summary
Hyperspectral and Multispectral Image Fusion based on a Sparse Representation
Qi Wei, José Bioucas-Dias, Nicolas Dobigeon, Jean-Yves Tourneret
TL;DR
HS/MS fusion is an ill-posed problem for which conventional pansharpening approaches can be inapplicable or inefficient. The paper learns image-dependent dictionaries and code supports, then solves the resulting sparse-regularized problem by alternating ADMM-based image updates and least-square coding; simulations report smaller spatial error and spectral distortion with manageable computational complexity.
Problem
HS/MS fusion is generally ill-posed, and many pansharpening methods are inapplicable or inefficient because both spatial and spectral information occur in the multi-band observations.
Method
The method learns dictionaries and active-code supports from observed images, then alternates ADMM-based target-image optimization with least-square coding under sparse regularization.
Results
The proposed method offered smaller spatial error and smaller spectral distortion than four state-of-the-art fusion methods, with manageable computation complexity.
Takeaways & Limitations
The specific sparse prior regularized the fusion optimization and was attributed with improving spatial and spectral results while retaining manageable computation complexity.
Abstract
from arXiv · showhide
This paper presents a variational based approach to fusing hyperspectral and multispectral images. The fusion process is formulated as an inverse problem whose solution is the target image assumed to live in a much lower dimensional subspace. A sparse regularization term is carefully designed, relying on a decomposition of the scene on a set of dictionaries. The dictionary atoms and the corresponding supports of active coding coefficients are learned from the observed images. Then, conditionally on these dictionaries and supports, the fusion problem is solved via alternating optimization with respect to the target image (using the alternating direction method of multipliers) and the coding coefficients. Simulation results demonstrate the efficiency of the proposed algorithm when compared with the state-of-the-art fusion methods.
I. INTRODUCTION
HS/MS fusion addresses the need to combine complementary spectral and spatial information, but differs from pansharpening and is generally ill-posed. The paper proposes image-adaptive sparse regularization with learned dictionaries and an alternating optimization strategy.
- HS sensors provide abundant spectral information, but their spatial resolution is generally limited, motivating fusion with MS or PAN images.
- HS/MS fusion differs from pansharpening because both spatial and spectral information are contained in the multi-band images.
- Many pansharpening methods are inapplicable or inefficient for the HS/MS fusion problem.
- Because the fusion problem is generally ill-posed, Bayesian inference can regularize it through an appropriate prior for the scene.
- The proposed method learns dictionaries from observed images and incorporates sparse regularization within a constrained optimization framework.
- Fixed code supports enable alternating optimization of the projected target image and sparse code, using SALSA/ADMM for the image and least squares for coding.
II. PROBLEM FORMULATION
The formulation models HS and MS observations as degraded, noisy versions of a full-resolution target image. Subspace learning reduces the unknown image to lower-dimensional coefficients, yielding a more tractable fusion problem.
- A. Notations and observation model: The HS observation is modeled as a blurred and down-sampled version of the target image, whereas the MS observation is spectrally degraded.
- A. Notations and observation model: The full-resolution target image X has mλ spectral bands and n pixels, while the observed HS and MS images are represented band-by-band as matrices.
- A. Notations and observation model: The observation model includes a band convolution operator B, a down-sampling matrix S with factor d, and HS/MS noise matrices.
- A. Notations and observation model: Each noise column is assumed to follow a band-dependent Gaussian distribution, an assumption used to formulate the likelihood and optimization algorithms.
- B. Subspace learning: The target spectrum at each spatial location is assumed to lie in a subspace whose dimension is much smaller than the number of HS bands.
- B. Subspace learning: Using an orthogonal basis H, the target is written as X = HU and the projected coefficients satisfy U = H^T X.
- B. Subspace learning: Estimating U reformulates the fusion problem through the MS equation YM = RHU + NM.
- B. Subspace learning: The subspace basis H is obtained by principal component analysis of the HS data, reducing computational burden and bypassing possible spectral-dependency singularity.
III. PROPOSED FUSION RULE FOR MULTISPECTRAL AND HYPERSPECTRAL IMAGES
The proposed fusion rule treats hyperspectral–multispectral fusion as a potentially ill-posed inverse problem and regularizes it with an adaptive sparse image representation. The MAP formulation combines sensor-data fidelity with learned, patch-based dictionaries and sparse codes.
- Ill-posed inverse problem: Fusion from a single observation is ill-posed or under-constrained, whereas multi-source fusion may be ill-posed or well-posed depending on subspace dimension and spectral-band count.The paper focuses on the under-determined case, which includes pansharpening.
- MAP optimization: The MAP estimator is obtained by minimizing two data-fidelity terms for the MS and HS observations plus a regularization penalty weighted by λ.λ controls the importance of regularization relative to data fidelity.
- Sparse regularization: Learning dictionaries from observed images is preferred because predefined wavelet, curvelet, or DCT bases may not match natural or remote-sensing image structure.The passages state that image-dependent dictionary learning generally improves representation.
- Sparse regularization: The sparse regularizer represents target-image patches as combinations of a few atoms from learned overcomplete dictionaries.Each band has a dictionary and corresponding code matrix, while overlapping patches are recombined by averaging.
- Sparse regularization: The formulation decomposes each spectral band into overlapping patches and uses the adjoint patch operator so that P(P∗(X)) = X.The operator P averages overlapping patches to reconstruct each band.
C. Dictionary learning step
The paper learns over-complete dictionaries and sparse-code supports from a rough estimate of the latent subspace, then incorporates them into the fusion regularization. A hierarchical Bayesian formulation allows these learned quantities to inform the estimation rather than remaining fixed independently.
- Dictionary learning step: Dictionaries and associated sparse-code supports are learned from a rough estimate of the latent subspace obtained using the observed HS and MS images.The estimate is decomposed into overlapping patches before dictionary training and support identification.
- Patch decomposition: Overlapping patches are used to decompose each estimated band and prevent block artifacts when reconstructing the image representation.The patch matrix contains npat overlapping patches for each band.
- Dictionary learning step: Online dictionary learning trains an over-complete dictionary for each subspace band, motivated by its computational efficiency and empirically relevant representations.The dictionary associated with band Ui is trained separately using ODL.
- Sparse coding: Each patch is represented with at most K atoms, where K is much smaller than the dictionary size, and orthogonal matching pursuit can solve the resulting ℓ0-constrained regression.The positions of non-zero code elements define the supports used later in optimization.
- Hierarchical formulation: Fixed dictionary and code hyperparameters can strongly affect representation accuracy, motivating a second hierarchical Bayesian level for their treatment.The paper notes that a standard fixed choice can reduce the relevance of the regularization term.
- Regularization: Fixing learned dictionaries and supports makes the regularization separable across bands and reduces coefficient optimization to an ℓ2-norm problem over non-zero code elements.The sparse-code supports are enforced by updating only the non-zero elements of each code column.
IV. ALTERNATE OPTIMIZATION
The fusion optimization is a large constrained quadratic problem, so the paper alternates between latent-image and coefficient updates. ADMM/SALSA handles the image update, while least squares updates the non-zero coefficients under fixed supports.
- Problem and strategy: The fusion problem is difficult because of its large dimension and non-diagonalizable operators, motivating an alternating optimization strategy.The alternating scheme is described as a simple block coordinate descent method.
- Latent-image update: Conditioned on the coefficients, the latent image is optimized efficiently with ADMM, whose convergence is established in the convex case.The objective is convex with respect to the latent image and coefficients, though not strictly convex.
- Coefficient update: Conditioned on the latent image, coefficient optimization becomes a least-squares regression over the non-zero elements selected by the support constraint.The support constraint is handled by updating only those non-zero elements.
- Convergence: Adding a very small quadratic term to the coefficients makes the algorithm’s solution unique and yields linear ADMM convergence.The paper reports convergence to a unique point even when the added coefficient term is zero.
- SALSA implementation: SALSA implements the ADMM image-update subproblem using variable splittings and scaled Lagrange multipliers.The split variables enforce the terms involving UB and repeated copies of U.
- SALSA implementation: SALSA has guaranteed convergence, and the latent-image optimization can be solved efficiently in the Fourier domain.The SALSA sub-iterations are summarized in the paper’s second algorithm.
B. Patchwise Sparse Coding
Patchwise sparse coding reformulates the coefficient update through patch extraction and reconstruction operators. With fixed dictionaries and supports, only the active coefficients require optimization.
- Patchwise formulation: The coefficient update conditioned on the latent image is formulated as a constrained optimization over patchwise sparse codes.The support constraint sets coefficients outside each learned support to zero.
- Operator formulation: The patch-to-image operator is linear and satisfies P(P*(X)) = X, enabling an equivalent reformulation of the coefficient problem.This identity connects patch-domain coding with image-domain reconstruction.
- Operator approximation: Using P*(P(·)) ≈ cI(·), the paper approximates the patchwise problem with a simpler optimization problem.Here c is a constant and I is the identity operator.
- Support constraint: Enforcing the support constraint requires updating only the non-zero elements of each coefficient column.Inactive elements remain fixed at zero.
- Implementation: The vectors used in the patchwise formulation can be computed in parallel, while projection operators need to be calculated only once after dictionary and support learning.The projection operator is symmetric because it is a projection.
C. Complexity Analysis
The proposed method is evaluated on synthetic HS/MS observations derived from the University of Pavia scene. Its complexity is reduced by operating in a lower-dimensional subspace, while the main computational costs depend on SALSA iterations and patchwise coding.
- Complexity: O(nit emλn log(emλn)) is the stated SALSA complexity, where nit denotes the number of SALSA iterations.Patchwise sparse coding additionally costs O(Knpnpat emλ).
- Complexity: Operating in an emλ-dimensional subspace instead of the original mλ-dimensional space greatly decreases the complexity of SALSA and sparse coding.The reduction applies to both major computational stages.
- Simulation scenario: The simulated low-spatial-resolution HS image is produced with a 5 × 5 exponential filter and four-pixel down-sampling in both spatial directions.A four-band MS image is generated using IKONOS-like spectral responses.
- Noise model: Both simulated observations contain zero-mean additive Gaussian noise, with band-dependent HS SNRs and a 30 dB SNR for every MS band.The HS SNR is 35 dB for the first 127 bands and 30 dB for the remaining 50 bands.
B. Learning the Subspace, the Dictionaries and the Code Supports
The method learns a low-dimensional subspace, over-complete dictionaries, and sparse code supports from the observed imagery. These learned structures guide reconstruction while balancing patch representation, computational cost, and over-fitting risk.
- Subspace: PCA constructs the transform matrix H from eigenvectors associated with the largest eigenvalues of the empirical covariance matrix Υ.The eigenvalues are ordered from largest to smallest, and the top e_mλ components define the subspace.
- Subspace: For the Pavia image, e_mλ = 5 eigenvectors contain 99.9% of the information.This illustrates the strong dimensionality reduction used for the target image.
- Dictionaries: ODL learns dictionaries from a rough projected-image estimate, using 6 × 6 patches and 256 atoms in this study.Because the number of atoms exceeds the patch dimension, the dictionary is over-complete.
- Dictionaries: Smaller patches approximate more objects, whereas larger patches capture textures and edges but require more atoms and computation.Patch size is selected empirically, so the representation involves a fidelity–cost trade-off.
- Code supports: The learned dictionary atoms capture spatial properties of the target image, with common atoms representing recurring details and other atoms encoding patch-specific details.The corresponding supports show that some atoms are frequently used while others are rarely active.
- Code supports: The sparsity level K controls the maximum number of atoms per patch; overly large K can over-fit when too many parameters are estimated from limited samples.Supports are jointly estimated with the projected image U.
1) RMSE:
The evaluation uses RMSE, SAM, UIQI, ERGAS, and DD to compare the proposed fusion method with four state-of-the-art algorithms on the Pavia dataset. The proposed method is reported as visually competitive and best across the listed quantitative measures.
- RMSE: RMSE measures similarity between the target and fused images, with smaller values indicating better fusion quality.The evaluation also considers SAM, UIQI, ERGAS, and DD.
- Other quality measures: SAM measures spectral distortion between actual and estimated images, and smaller absolute values indicate less distortion.Overall SAM is obtained by averaging pixelwise spectral angles, expressed in degrees.
- Other quality measures: ERGAS measures global spectral distortion relative to the reference image, while smaller ERGAS values indicate less distortion.Its formulation accounts for the MS-to-HS pixel-size ratio and the means of HS bands.
- Implementation and convergence: The SALSA regularization parameter is set to μ = 0.05, while the choice of μ remains an open tuning issue despite convergence for any μ > 0 when a solution exists.The reported fusion performance is not significantly affected by μ as long as it is positive.
- Comparison: The proposed method performs competitively visually and provides the best results across RMSE, UIQI, SAM, ERGAS, and DD in the Pavia comparison.Figure 6 compares the reference, input HS and MS images, three competing methods, and the proposed method.
E. Selection of the regularization parameter λ
The regularization parameter λ is evaluated empirically because no single value is optimal across all quality measures. The proposed method performs favorably across a broad λ range and on additional fusion experiments.
- Parameter selection: No single λ value is optimal for all quality measures.Performance is evaluated as a function of λ in Fig. 7.
- Parameter selection: λ = 5 provides the best fusion results in terms of RMSE for the Section V-D simulation.
- Comparative performance: Across a wide range of λ, the proposed method outperforms the other four methods.
- Additional experiments: The method was also tested on a 128 × 128 × 177 hyperspectral Moffett field image and compared with several fusion methods.The experiments include reference, MS, HS, and fusion-result images, with further quantitative results reported in a table.
- Conclusions: The conclusions report smaller spatial error, smaller spectral distortion, and manageable computational complexity than four state-of-the-art fusion methods.The improvement is attributed to the specific sparse prior used to regularize the inverse problem.
- Limitations and future work: Future work includes estimating λ within the fusion algorithm and jointly updating the dictionary with the target image.
APPENDIX
The appendix models the target pixels and multispectral observations with Gaussian conditional distributions. It approximates the projected target image using a conditional mean derived from spatial and covariance estimates.
- Gaussian model: The conditional distribution of U given YM is Gaussian under the stated Gaussian prior and observation assumptions.
- Variables: The target image U is represented by pixel vectors ui, while YM is represented by multispectral pixel vectors yM,i.
- Approximation: The rough approximation ˜U is defined as the conditional mean of the projected target image U.
- Approximation: Computing the conditional mean requires the prior mean and covariance terms involving ui and yM,i.
- Approximations: The prior mean is approximated from spatially interpolated observed HS imagery, while the MS-image prior mean uses spatial smoothing.
- Approximations: Covariance estimation is performed at the lower resolution of the observed HS imagery.