Deep Gradient Projection Networks for Pan-sharpening

In plain words

Satellites often take two pictures of the same place: a sharp black-and-white picture and a blurrier colour picture with several bands. Pan-sharpening combines them into one colour picture that is sharp everywhere.

Most neural networks for this task learn the mapping from examples and ignore how the satellite made the two pictures. This paper starts from the physics. The blurry colour picture is a blurred, shrunk copy of the true sharp colour image. The black-and-white picture is a mix of the colour bands. The authors write one optimisation problem for each of these two facts, and a learned rule stands in for what a good image looks like.

Each problem is solved by a standard step-by-step method: check how well the current guess matches the observation, correct the guess, then clean it up. The authors turn each step into a small network block and stack the blocks into one network. The result is a network whose stages can be read as the steps of that method, and in their experiments it beats the earlier methods they compare against on three satellites.

TL;DR

GPPNN is a model-based pan-sharpening network. It unrolls a gradient projection solver for two observation models, one for the low-resolution multispectral (LRMS) image and one for the panchromatic (PAN) image, into alternating MS and PAN blocks. Use it as a reference design for physics-informed fusion networks, or as a baseline on Landsat8, QuickBird or GaoFen2 pan-sharpening. All results are author-reported. The official code has no licence file and some choices differ from the paper.

Problem

Pan-sharpening fuses a low-resolution multispectral image with a high-resolution PAN image §1. Deep networks for it need many training samples and do not model how the two images are formed, so they are hard to interpret §1. Hyperspectral pan-sharpening networks such as MHNet rely on a low-rank property that is weak when there are no more than 10 bands, and the authors report that MHNet performs badly on pan-sharpening §1 §2.3.

Method

1. Observation models §3 §3.1. The LRMS image L is the true HRMS image H, blurred and downsampled. The PAN image P is H projected onto one band through a spectral response S. The ratio is r = M/m = N/n. Notation: B bands, D downsampling, K blur, S spectral response §3.1.
2. Two problems instead of one §3.1. A single penalised problem with one prior h is first written eq. 2. The paper then splits it into an LRMS-aware problem and a PAN-aware problem, each with its own deep prior eq. 3a eq. 3b. The ablation favours the split §4.4, configuration III.
3. Gradient projection §3.2 §3.3. Each problem is solved by a gradient step with step size ρ followed by a proximal operator eq. 4. The data gradients are −(DK)ᵀ(L − DKH) for the LRMS problem and −(P − HS)Sᵀ for the PAN problem §3.2 §3.3.
4. MS block §3.2. The LRMS gradient step is split into four steps eq. 5a–5d. Each step becomes a network operation eq. 6–9: a learned cascade replaces the blur-and-downsample operator, a second learned cascade replaces its transpose, and a learned cascade replaces the proximal operator. The two cascades are not forced to be transposes of each other §3.2. Section 4.4 tests the forced version.
5. PAN block §3.3. The PAN gradient step is split into four steps eq. 10 eq. 11a–11d and mapped to a network eq. 12a–12d. Because the PAN image is a linear combination of the bands, the convolutions that simulate and lift the PAN image use kernel size 1 §3.3.
6. Stacking §3.4. The estimate starts as the bicubic upsampling of the LRMS image. The network then applies K pairs of blocks, each an MS block followed by a PAN block eq. 13 eq. 14. The last estimate is the output. Each pair has its own weights; the ablation that shares them is worse §4.4.
7. Training §3.5. The loss is L1 between the output and the ground truth. Training pairs are made with the Wald protocol §3.5. Optimiser: Adam, 100 epochs, learning rate 5×10⁻⁴, batch size 16 §3.5.

Key equations (from the paper):

\min_{\bm{H}}\ \tfrac12\|\bm{L}-\bm{DKH}\|_2^2+\lambda h_l(\bm{H}) \quad [\text{eq. 3a}]

\min_{\bm{H}}\ \tfrac12\|\bm{P}-\bm{HS}\|_2^2+\lambda h_p(\bm{H}) \quad [\text{eq. 3b}]

\bm{H}^{(t)}=\mathrm{prox}_{h_l}\!\big(\bm{H}^{(t-1)}-\rho\nabla f(\bm{H}^{(t-1)})\big) \quad [\text{eq. 4}]

\mathtt{Conv}(\cdot;c^{in},c^{mid},c^{out})=\mathtt{Conv}\big(\mathtt{ReLU}(\mathtt{Conv}(\cdot;c^{in},c^{mid}));c^{mid},c^{out}\big) \quad [\text{eq. 1}]

Network blocks as pseudocode. The code shown is the paper's structure, written in our words, with the code's choices noted in the comments.

# GPPNN forward pass  [§3.4] [code: models/GPPNN.py#L91-L106]
# Input:  lrms (B, bands, m, n), pan (B, 1, M, N), per-sample normalised
# Output: hrms (B, bands, M, N)

def gppnn_forward(lrms, pan, num_pairs):
    # Initial estimate: bicubic upsampling of the LRMS image  [§3.4]
    hrms = resize_bicubic(lrms, pan.shape[-2:])
    for pair in range(num_pairs):  # K = 8 in the paper  [§3.4]
        # LRMS-aware update  [§3.2] [eq. 13]
        hrms = ms_block(hrms, lrms)
        # PAN-aware update  [§3.3] [eq. 14]
        hrms = pan_block(hrms, pan)
    return hrms

# MS block: LRMS-aware gradient projection step  [§3.2]
# [code: models/GPPNN.py#L52-L60]
# Input:  hrms (B, bands, M, N), lrms (B, bands, m, n)
# Output: hrms (B, bands, M, N)

def ms_block(hrms, lrms):
    # Simulate an LRMS image: learned stand-in for D K  [eq. 6]
    # cascade(x) = conv (B -> C), ReLU, conv (C -> B)  [eq. 1]
    simulated_lrms = resize_bicubic(cascade(hrms), lrms.shape[-2:])
    # Residual between observed and simulated LRMS  [eq. 7]
    lrms_residual = lrms - simulated_lrms
    # Lift the residual to HRMS size  [eq. 8]
    # The step size rho of the paper has no parameter in the code
    hrms_residual = resize_bicubic(cascade(lrms_residual), hrms.shape[-2:])
    # Gradient step, then the learned proximal unit  [eq. 9]
    return cascade(hrms + hrms_residual)

# PAN block: PAN-aware gradient projection step  [§3.3]
# [code: models/GPPNN.py#L73-L78]
# Input:  hrms (B, bands, M, N), pan (B, 1, M, N)
# Output: hrms (B, bands, M, N)

def pan_block(hrms, pan):
    # Simulate the PAN image from the bands  [eq. 12a]
    # Kernel size 1; conv_1x1(x, c_in, c_out) is conv, ReLU, conv
    simulated_pan = conv_1x1(hrms, bands, 1)
    # Residual between observed and simulated PAN  [eq. 12b]
    pan_residual = pan - simulated_pan
    # Spread the residual over the bands  [eq. 12c]
    hrms_residual = conv_1x1(pan_residual, 1, bands)
    # Gradient step, then the proximal unit  [eq. 12d]
    # The code also uses kernel size 1 here (see Paper vs code)
    return conv_1x1(hrms + hrms_residual, bands, bands)

Hyperparameters

Values stated in the paper, with the code's value where it differs or adds detail.

Results

Author-reported. The setup is three satellites, with the splits in tab. 1. Metrics are PSNR, SSIM, SAM and ERGAS §4.1. Baselines are six deep methods (MIPSM, DRPNN, MSDCNN, RSIFNN, PanNet, MHNet) and seven classical methods (BDSD, Brovey, GS, HPF, IHS, Indusion, SFIM) §4. The paper does not say whether baselines were retrained. Hardware is an Intel i7-9700K with an RTX 2080 Ti §4.

tab. 1 Dataset information.

tab. 2 PSNR on the validation sets for different K (C fixed at 64) and C (K fixed at 8). Values are as printed.

tab. 3 Test results. Columns per satellite are PSNR (higher is better), SSIM (higher), SAM (lower) and ERGAS (lower). The paper bolds GPPNN as the best value in all twelve columns. Values are as printed.

tab. 4 Ablation on Landsat8. Each row changes one design choice from GPPNN. The paper's column labels are paraphrased here: proximal unit, shared weights, naive-model block, forced transpose.

Figures fig. 3, fig. 4 and fig. 5 are visual comparisons. Each panel's caption gives a PSNR for that one image, so these numbers are not the means in tab. 3. For example, the Landsat8 panel for GPPNN is 24.40 dB fig. 3, while the Landsat8 mean in tab. 3 is 38.99 dB. The paper does not say which test image the figures show.

Paper vs code

Official repository: github.com/xsxjtu/GPPNN, linked from the abstract. Clone at raw/code/2103_04584/, commit d3443120a0ed1c9a04e740e367a60744ef80c84f. The repository has no LICENSE file, so the code is cite-only.

• Step size ρ. The paper keeps ρ in the residual lifts eq. 8 eq. 12c and the gradient steps eq. 4. It gives no value §3.2. The code has no ρ term GPPNN.py [L76]. The units have no bias and use ReLU, so a positive ρ can be folded into the last convolution. The code learns that scale instead of fixing it. This is an inference from the structure, not a stated equivalence.
• Kernel sizes. The paper states kernel size 1 only for eq. 12a and eq. 12c §3.3. The code uses kernel size 1 for all three PAN-block units, including the proximal unit GPPNN.py. The MS-block kernel is not stated; the code uses 3×3 with padding 1 GPPNN.py [L48-L50].
• k_LR and k_PAN. The paper sets k_LR = 3 and k_PAN = 1 §3.5 but never defines them. Neither name appears in the code. The code's kernel sizes are 3 for the MS block and 1 for the PAN block. The most likely reading is "kernel size", but this is inference, not stated.
• Padding and bias. Not stated in the paper. The code pads to keep the size and uses no bias terms GPPNN.py.
• Resizing. The paper says bicubic for the LRMS downsampling §3.2 and does not state the upsampling method. The code uses bicubic with align_corners=True for every resize, including the initial estimate GPPNN.py [L56] [L58] [L100].
• Learning-rate schedule. The paper states a learning rate of 5×10⁻⁴ §3.5. The code builds a StepLR (step 20, factor 0.5) train_GPPNN.py but never calls scheduler.step(). The only step call is optimizer.step() train_GPPNN.py. The rate stays at 5×10⁻⁴.
• Checkpoint selection on the test set. The paper picks K and C on the validation set §4.1 and does not say how the epoch is chosen. The code stores the test PSNR in the variable it uses to select the checkpoint train_GPPNN.py. It saves the best one when that value improves train_GPPNN.py and reports the test metrics from it train_GPPNN.py [L245-L246]. The validation PSNR is logged but not used for selection train_GPPNN.py.
• Normalisation. The paper says each patch is divided by its maximum §4.1. The code divides the LRMS, PAN and ground-truth tensors each by their own per-sample maximum, in training, validation and test train_GPPNN.py [L177-L179] [L194-L196] utils.py. The metrics are computed on these normalised tensors train_GPPNN.py.
• Training pairs. The paper uses the Wald protocol §3.5 and says nothing about test-time inputs. The code's training data uses nearest-neighbour resizing with no mode given utils.py [L118-L128]. The PAN image is resized to the MS grid, and the LRMS input is the ground truth downsampled by r utils.py. For validation and test, the LRMS and PAN inputs are nearest-downsampled from the native MS and PAN images utils.py.
• Augmentation and stride. The data call passes aug_times=1 and stride=32 train_GPPNN.py. Rotations are off and the patches do not overlap. The paper states neither.
• ERGAS ratio. get_metrics_reduced calls ergas without a scale, so the default of 4 applies to all satellites metrics.py [L17]. The paper gives no ratio, but its Table 1 gives PAN/MS = 2 for Landsat8. If the paper's numbers come from this code, the Landsat8 ERGAS is computed with ratio 4. This cannot be checked from the paper.
• SAM units. The code returns SAM in radians utils.py. The paper does not state units. The code also computes CC metrics.py, which the paper does not report.
• Matches. Layers, channels, bands and scale match §3.4 §4.2 tab. 1 train_GPPNN.py [code: models/init.py#L3-L18]. The optimiser, L1 loss, batch size and epoch count match §3.5 train_GPPNN.py. The block order (MS then PAN, per pair) and the separate weights per pair match GPPNN.py [L87-L88].
• Repository contents. Only one sample per split is shipped, and only for Landsat8 Readme.txt. The GaoFen2 and QuickBird data, and the full splits in tab. 1, are not in the repository.
• Link mismatch. The README links arXiv 2103.05946, not this paper's 2103.04584 README.md. I did not check 2103.05946. This card uses the repository linked from the abstract of 2103.04584.

Pitfalls

• Test-set checkpoint selection. The published code's test numbers come from the epoch with the best test PSNR (see Paper vs code). If you reproduce these results, select on validation.
• Constant learning rate. The code keeps 5×10⁻⁴ throughout, not the step schedule it builds. Schedule-based comparisons with the code will differ.
• Per-tensor normalisation. The LRMS and PAN inputs are scaled independently, so their relative brightness is lost. Metrics are computed in the normalised space; the code's inverse-normalisation helper is not used.
• Missing ρ. An implementation that follows the paper literally needs a step size the paper does not give. The code learns the scale instead.
• Interpolation. Bicubic resizing with align_corners=True differs from the common default. Frameworks that default to align_corners=False will not match the code's behaviour without a change.
• Unstated values. The paper does not state: the MS-block kernel size, the upsampling method, the Adam betas, weight initialisation, the training stride or augmentation, the test-time construction of inputs, or the number of training patches. The code answers some of these (see above), not all.
• Patch sizes. "32×32" is the LRMS input. The ground truth in the code is r × 32 (64 for Landsat8, 128 for GaoFen2 and QuickBird).
• Paper text vs its table. The text says C = 64 is the best width §4.2. tab. 2 shows C = 128 slightly ahead on QuickBird validation (30.5979 vs 30.5636). For Landsat8, K = 10 is ahead on validation (40.1336 vs 40.0368), and the paper picks K = 8 as a balance §4.2.
• Single-satellite models. The paper states that a network trained on one satellite does not transfer to another, because each satellite has its own D, K and S §5.
• Cost. Every block works at the PAN resolution, so memory and compute grow with the PAN image size. The paper gives no timings §4.
• Figures are single images. Per-panel PSNR values in fig. 3 to fig. 5 are not comparable with the means in tab. 3.
• Author-reported. All numbers are the authors'. The paper does not say whether baselines were retrained or taken from their papers.

Lineage

• Builds on: model-driven deep learning and algorithm unrolling §2.3. The gradient projection method is cited as [Gradient_projection] §3.2. No cards are used for either.
• Compared against: MIPSM, DRPNN, MSDCNN, RSIFNN, PanNet, MHNet, BDSD, Brovey, GS, HPF, IHS, Indusion and SFIM §4. Nothing from them is used here.
• Related: MHNet, a model-driven hyperspectral pan-sharpening network §2.3. The paper argues it does not transfer well to multispectral pan-sharpening §1.
• Superseded by: not known from content in this wiki.