Multimodal Remote Sensing Benchmark Datasets for Land Cover Classification with A Shared and Specific Feature Learning Model

In plain words

Satellites and aircraft record the same ground in different ways. One sensor sees fine colour detail, another sees radar reflections, a third sees height. The usual way to combine them is to paste the readings side by side and feed the result to a classifier. That can drown out what each sensor is good at, and it struggles when the sensors are very different from each other.

This paper proposes a method, called S2FL, that splits the learned features of all sensors into two parts. One part is shared: what the sensors agree on for the same kind of ground. The other part is specific: what only one sensor can tell you. A graph that pulls together pixels of the same land-cover class guides the shared part. The authors also release three new benchmark datasets that combine hyperspectral data with multispectral, radar or height data, so that such methods can be compared. On those datasets the method classifies land cover more accurately than the earlier methods it is compared with.

TL;DR

S2FL learns a linear projection per sensor set that separates a shared component from modality-specific components, fitting class labels through a regression matrix. Both projections are kept orthonormal and are solved by an alternating scheme that uses ADMM for the shared part. It is compared against JDR-PCA, USMA, SMA, l2-CoSpace and l1-CoSpace on three datasets, with a nearest-neighbour classifier. The results are author-reported, and the official code covers only the Houston2013 setting.

Problem

Multimodal remote-sensing data differ in sensor, resolution and content, so a single compact and discriminative representation is hard to build §1. Most earlier feature methods either concatenate modality features, which can inadequately fuse information and damage each modality's components §1, or align a common feature space across modalities §2.2. Alignment-only methods capture the shared part but not the modality-specific part §2.2. The field also lacked open multimodal benchmarks, especially with three modalities §1.

Method

Overview. S2FL decomposes multimodal data into a modality-shared subspace, obtained with one projection Θ0, and modality-specific subspaces, obtained with one projection Θk per modality. The shared projection is pulled toward a manifold alignment term, and all projections are orthonormal §3.1. The learned features are then classified with a nearest-neighbour rule §4.2.1.

Notation. Each modality k has d_k channels over N pixels; X_k is d_k × N, Y is the C × N one-hot label matrix, d_s is the feature dimension, and P is a C × d_s regression matrix §3.2 tab. 1. The block matrix X̃ places each modality's pixels in its own rows and zeros elsewhere, and Θ = Θ0 + [Θ1, ..., ΘK] §3.3. The graph Laplacian is L = D − W over the 2N samples used in the two-modality case §3.2.

Objective. S2FL minimises the three-term problem below eq. 1. The first term is the least-squares fit of labels to the projected features. The second term is a Frobenius penalty on P, and the third is a graph-alignment penalty on the shared part only.

\min_{P,\Theta_0,\{\Theta_k\}}\; \tfrac12\lVert\tilde Y-P\Theta\tilde X\rVert_F^2 +\tfrac\alpha2\lVert P\rVert_F^2 +\tfrac\beta2\operatorname{tr}\!\left(\Theta_0\tilde X L (\Theta_0\tilde X)^\top\right) \quad\text{s.t. }\Theta_k\Theta_k^\top=I,\; k=0,\dots,K eq. 1

Graph. The adjacency W has K² blocks §3.3. Intra-modality blocks use k-nearest-neighbour edges with a Gaussian weight exp(−‖x_i − x_j‖²/σ²) eq. 2 §3.3 §4.2.3. Inter-modality blocks link pixels of the same class with weight 1/N_C, where N_C is the number of samples in class C eq. 3. The code builds the graph on training pixels only demo.m.

Optimisation. The problem is non-convex, so the method alternates over P, Θ0 and the Θk, each step being a convex subproblem §3.4 alg. 1.

1. Update P in closed form, a ridge regression on the projected features §3.4.1 eq. 5.
2. Update the shared projection Θ0 with the ADMM solver of alg. 2 (below), using the graph term §3.4.2 eq. 6 eq. 7 eq. 8.
3. Update each specific projection Θk with the same solver, with the graph term removed §3.4.3.
4. Stop when the relative change of the objective falls below ζ = 1e-4 alg. 1.

The closed-form P update is eq. 5. The shared-projection ADMM introduces auxiliary variables H (for Θ0 X̃) and G (for Θ0), with Lagrange multipliers Λ1 and Λ2 and a penalty μ eq. 8. Its updates are: H in closed form eq. 12, Θ0 in closed form eq. 10, G by an SVD-based orthonormal factor eq. 13 eq. 14 eq. 15, then the multiplier updates eq. 16 and the penalty growth μ ← min(ρμ, μmax) eq. 17.

# S2FL global optimisation  [alg. 1]
# Input:  Y_all (C, 2N) one-hot labels repeated per modality
#         X_all (D, 2N) block-padded features; L (2N, 2N) Laplacian
# Output: P (C, d_s); theta = shared + specific (d_s, D)

def s2fl(Y_all, X_all, L, alpha, beta, max_iter):
    # Start with theta from LPP and a zero shared part  [alg. 1]
    shared = zeros(d_s, D)
    specific = lpp_projection(X_all, L)           # (d_s, D)
    prev_loss = None
    for t in range(max_iter):
        theta = shared + specific                 # (d_s, D)
        features = theta @ X_all                  # (d_s, 2N)
        # Closed-form regression matrix  [§3.4.1] [eq. 5]
        P = (Y_all @ features.T) @ inverse(
            features @ features.T + alpha * eye(d_s))
        # Shared projection by ADMM, with the graph term  [§3.4.2] [alg. 2]
        residual = Y_all - P @ specific @ X_all
        shared = solve_shared(residual, P, X_all, beta, L)
        # Specific projection, same solver without the graph term  [§3.4.3]
        residual = Y_all - P @ shared @ X_all
        specific = solve_specific(residual, P, X_all)
        # Relative objective change below 1e-4 ends the loop  [alg. 1]
        loss = objective(Y_all, P, shared + specific, X_all, L, alpha, beta)
        if prev_loss is not None:
            if abs(loss - prev_loss) / abs(prev_loss) < 1e-4:
                break
        prev_loss = loss
    return P, shared, specific

# ADMM for the shared projection, closed-form updates  [alg. 2] [eq. 9-17]
# Input:  residual (C, 2N); P (C, d_s); X_all (D, 2N); L (2N, 2N); beta
# Output: shared projection (d_s, D)

def solve_shared(residual, P, X_all, beta, L):
    # Graph term of the objective  [eq. 10]
    graph_term = X_all @ L @ X_all.T              # (D, D)
    # Initial values: zero projection, zero auxiliaries and multipliers
    shared = zeros(d_s, D)
    # Auxiliaries stand for shared @ X_all (aux_h) and shared (aux_g)
    aux_h = zeros(d_s, 2N)
    aux_g = zeros(d_s, D)
    # Dual variables Lambda_1 and Lambda_2
    dual_h = zeros(d_s, 2N)
    dual_g = zeros(d_s, D)
    mu = 1e-3
    for it in range(max_iter):
        # Update H in closed form  [eq. 12]
        aux_h = inverse(P.T @ P + mu * eye(d_s)) @ (
            P.T @ residual + mu * shared @ X_all - dual_h)
        # Update the shared projection in closed form  [eq. 10]
        numerator = mu * aux_h @ X_all.T + dual_h @ X_all.T
        numerator = numerator + mu * aux_g + dual_g
        denominator = mu * X_all @ X_all.T + mu * eye(D)
        denominator = denominator + beta * graph_term
        shared = numerator @ inverse(denominator)
        # Update G: SVD of shared - dual_g / mu, keep the orthonormal factor
        # SVD of shared - dual_g / mu  [eq. 14]
        U, _, Vt = svd(shared - dual_g / mu)
        # Orthonormal factor, eye(d_s, D) is the d_s x D identity  [eq. 15]
        aux_g = U @ eye(d_s, D) @ Vt
        # Dual ascent on both constraints  [eq. 16]
        dual_h = dual_h + mu * (aux_h - shared @ X_all)
        dual_g = dual_g + mu * (aux_g - shared)
        # Grow the penalty, capped at mu_max  [eq. 17]
        mu = min(1.5 * mu, 1e6)
        # Stop when both residuals are below 1e-6  [alg. 2]
        r_h = norm(aux_h - shared @ X_all)
        r_g = norm(aux_g - shared)
        if r_h < 1e-6 and r_g < 1e-6:
            break
    return shared

Figures: fig. 1 is a flowchart of the shared and specific subspaces for two modalities; fig. 2 shows an example of the adjacency matrix W.

Hyperparameters

Values for the Houston2013 setting are in tab. 5, and the parameter search grids in §4.2.3. Values for the other datasets are only in tab. 6 and tab. 9.

Results

Author-reported. Setup: the nearest-neighbour classifier of §4.2.1 on features learned by each method, with parameters chosen by 10-fold cross validation on the training set §4.2.3. Each comparison below is taken from the paper's own tables; no value was re-run.

Best baseline is l1-CoSpace in all three tables.

Houston2013 ablation, tab. 7 (OA %): HS only 69.21 (no feature learning); S2FL on HS 76.50; MS only 70.82; S2FL on MS 72.74; shared part only 78.77; specific part only 83.11; both parts without the orthogonality constraint 66.11; full S2FL with both parts and orthogonality 85.07.

Cross-modality, tab. 8 (OA %): trained on HS-MS, tested on HS only: 80.86, against 69.21 for HS features alone.

Augsburg modality combinations tab. 10 (OA %): S2FL on SAR+DSM 69.49; HS+DSM 81.11; HS+SAR 82.34; HS+SAR+DSM 83.36. For the same combinations l1-CoSpace reaches 81.17 on HS+DSM, above S2FL's OA on that pair; S2FL is higher on AA and kappa there.

The convergence of the relative objective per iteration is shown for the three datasets in fig. 3. Classification maps are in fig. 5 (Houston), fig. 7 (Berlin) and fig. 8 (Augsburg). The parameter sensitivity over k, σ, α, β and d on Houston is in fig. 6; the paper says the accuracy is stable over a wide range.

Paper vs code

Official code: danfenghong/ISPRS_S2FL at 7046f139f8d3204f51593e5945183fbd0d2a41eb, linked from the abstract. It contains the Houston2013 demo and the S2FL and CoSpace solvers. It has no Berlin or Augsburg scripts and does not ship the datasets: the demo loads .mat and .tif files that are not in the repo demo.m. The repo's CoSpace baseline files were not read for this card. The README states GPL v3 README.md and the repo has no LICENSE file, so the code is cited by path only.

Matches the paper:

• P update eq. 5 S2FL.m.
• Θ0 and H closed forms eq. 10 eq. 12 Solving_Theta00.m.
• μ0 = 1e-3, ρ = 1.5, μmax = 1e6, inner ε = 1e-6 alg. 2 Solving_Theta00.m.
• Outer ζ = 1e-4 alg. 1 S2FL.m.
• Houston α = 0.01, β = 0.1, d = 30 tab. 5 demo.m.

Discrepancies:

• β for Houston. The text gives the optimum as (q, σ, α, β, d) = (45, 10^0, 10^-2, 10^0, 30) §4.3.1, so β = 1. Table 5 and the code use β = 0.1 tab. 5 demo.m. The card uses 0.1.
• Orthonormality of specific projections. Eq. 1 writes Θk Θk^T = I for each k eq. 1. The code solves one specific projection for all modalities together, with one orthonormality constraint on the stacked matrix Solving_Theta12.m S2FL.m. Taken literally, the per-modality constraint cannot hold when a block has fewer features than d_s: Houston's MS has 8 bands with d = 30 §4.1.1 tab. 5, and Augsburg's SAR has 4 features with d = 20 §4.1.3 tab. 9. The code's joint form avoids this.
• Update order. Alg. 1 updates Θ0 and then each Θk alg. 1. The code updates the specific projection first and then the shared one S2FL.m.
• Initialisation. Alg. 1 sets Θ by LPP and Θ0 to zero alg. 1. The code uses the LPP projection only as the regressor for the first P update demo.m S2FL.m; both Θ parts start at zero S2FL.m Solving_Theta12.m. The inner ADMM also restarts from zero on every outer iteration Solving_Theta00.m.
• Kernel width. Eq. 2 uses exp(−‖x_i − x_j‖²/σ²) eq. 2. The code uses constructW with a heat kernel exp(−d²/(2t²)) where t = σ, on squared distances creatLap.m constructW.m constructW.m. At σ = 1 the two differ by a factor 2 in the exponent.
• Inter-modality weights. Eq. 3 gives weight 1/N_C per class eq. 3. The code gives weight 1 to all same-class pairs, then scales the whole W so that its entries sum to 2N demo.m. Class size is not normalised.
• Symmetry. The code symmetrises the kNN graph with max(W, W') constructW.m. The paper does not say how the kNN graph is symmetrised.
• Hidden preprocessing. Each band is min-max scaled over the whole image before the train/test split demo.m. The paper does not state normalisation, and the scaling uses every pixel, including test and background pixels.
• Combining the two halves. The code sums the projections of the two modality halves of each sample demo.m. This equals projecting the concatenated input, because the projection is linear. The paper does not state how the modalities are combined for classification (inferred from code).
• Hard-coded size. The demo hard-codes 2832 training samples demo.m and builds dense 2N x 2N matrices for W and L demo.m.

No discrepancy found for the NN classifier, which uses one neighbour and Euclidean distance §4.2.1 demo.m.

Pitfalls

• Scale of the inner stopping rule. The inner ADMM stops on absolute Frobenius residuals below 1e-6 alg. 2. The code's inputs are scaled to [0, 1] demo.m, so the threshold has a fixed meaning only for that scaling.
• Restarts are costly. Each outer iteration restarts both ADMM solvers from zero Solving_Theta00.m Solving_Theta12.m. Warm-starting is not stated in the paper. The paper gives the complexity per step, with SVDs and matrix inverses dominating §3.5.
• Memory grows with the square of the sample count. The 2N x 2N dense Laplacian is built in the demo demo.m. For Houston's 2832 training samples, 2N = 5664 and one dense double matrix is about 256 MB. The paper does not discuss memory.
• Non-convex; no proof. The problem is non-convex. The paper does not prove the ADMM's convergence and says a strict proof is still to be done §3.5.
• Maximum iterations. Alg. 1 takes a maxIter argument alg. 1, but the paper does not give a value. The demo uses 1000 demo.m.
• Table 5 vs Table 7. The full S2FL on Houston has AA 86.11 in tab. 5 and 86.01 in tab. 7, with the same OA and kappa.
• Table 8 conflicts with Table 5. The "Only MS" column of tab. 8 gives OA 70.01 and AA 80.46. The MS column of tab. 5 gives OA 70.82 and AA 74.75. The kappa of "Only MS" is printed as 75.09, which is outside the kappa scale used elsewhere. "Only MS" AA equals the S2FL-CML-MS AA. These look like transcription errors; no code for these runs exists in the repo.
• Classifier choice. The paper uses a 1-nearest-neighbour classifier on purpose, so that the feature quality is measured §4.2.1. Its numbers are not comparable with SVM, random-forest or deep-learning pipelines.
• Unstated Berlin and Augsburg settings. q and σ for S2FL on Berlin and Augsburg are not stated tab. 6 tab. 9. The DSM attribute profiles that turn 1 band into 3 are not specified §4.2.3.
• Resampling is described loosely. The Houston HS degradation is written as bilinear to 10 m and then nearest-neighbour upsampling §4.1.1. Implementations may differ in details.
• Labels and licences. Berlin and Augsburg labels come from OpenStreetMap §4.1.2 §4.1.3, and the paper does not measure label noise. The paper also does not state a licence for the released datasets; this card does not assess either point and uses no data.

Lineage

• Builds on: CoSpace (ℓ2) and ℓ1-CoSpace §2.2 §4.2.2: manifold-regularised common subspace learning with a label-based graph. S2FL adds a specific part per modality and orthonormal projections. Also uses LPP for initialisation alg. 1 and an ADMM variant after the inexact augmented Lagrangian method §3.5.
• Compared against: JDR-PCA, USMA, SMA, l2-CoSpace and l1-CoSpace §4.2.2 tab. 5.
• Superseded by: not stated in the paper, and no later card in this wiki.