Confidence-based Self-Corrective Learning: An Application in Height Estimation Using Satellite LiDAR and Imagery

1 Introduction

Climate action is an important component of the 2030 Sustainable Development Goals (SDGs) of the United Nations [UN, 2022]. Human activities have led to widespread transformations of the Earth's environment and climate [P¨ortner et al., 2019]. Observational evidence suggests that the climate system at high latitudes including the Arctic is undergoing significant shifts with faster increases in air temperatures

∗Corresponding author.

[Richter-Menge et al., 2019], which have led to changes including the mass loss of Greenland Ice Sheet and glaciers, permafrost thaw, decreasing snow and ice cover on land, and increasing vegetation cover.

Satellite remote sensing systems play a key role in monitoring essential climate variables in the polar regions, which are hard to access for traditional observation. Indeed, advances in sensing instruments have continued to demonstrate promising potentials. The recent launch of ICESat-2 satellites provides a revolutionary approach to measure the surface heights at high precision to help scientists quantify and better understand the ongoing changes from a third dimension. However, in order to obtain high-resolution heights at large scale with altimeter instruments, ICESat-2 has a very narrow footprint (i.e., width of a scan line) as shown in Fig. 1(a), leaving large empty gaps between scan lines. Filling the gaps will substantially improve the representativeness of the dataset.

We aim to expand the spatial footprints of the height data by combining the single-line height observations with largefootprint satellite imagery. This height estimation task has the following challenges. First, there are often only one or a few scan lines available in each satellite image tile for the corresponding timestamp, and all height observations are spatially constrained to the single or few lines, making it difficult for traditional interpolation methods to expand spatial footprints. Second, a large proportion of pixels in satellite imagery present variations in spectral characteristics (e.g., colors) that are irrelevant to height variations, which brings challenges in building the connections between the height data and spectral imagery. Finally, an image only reflects height variations in the covered region, and the base height can be highly uncertain to set. Different landscapes may also exhibit heterogeneous patterns over space [Xie et al., 2021].

While image-based height estimation has been widely studied, most existing works rely on high-contrast characteristics from man-made structures in urban areas with aerial images, which are not available in non-urban regions at large scale. Or they often require specific image types such as stereo images that are not available for most satellites [Qin, 2019]. Traditional spatial interpolation methods (e.g., Kriging) designed to fill data gaps cannot well capture complex nonlinear relationships between spectral features and heights. While their deep learning extensions (e.g., Kriging convolutional networks [Appleby et al., 2020]) enhance such ability, their performances suffer when observations are highly localized in space. More details are discussed in Sec. 3.

We propose a spatial self-corrective learning framework to expand the coverage of surface height data in high latitudes including polar regions with the following contributions:

• We propose a confidence-based pseudo-interpolation framework to correct lower-confidence predictions using higher-confidence predictions from spatially-adjacent areas. This helps suppress erroneous predictions caused by irrelevant feature variations from the input imagery.
• We propose a dynamic refinement strategy with a recurrent structure to allow the pseudo-interpolation to iteratively correct low-confidence predictions, allowing adaptive self-correction during test.
• We integrate an inverse-confidence-based ridge regression layer to utilize the limited observations from the scan line (i.e., a line crossing an image) to perform truth-based correction. We further improve the robustness of the correction by using it recurrently with random ensembles. The network parameters are trained to collaborate with the recurrent regression correction to better generalize over space.
• We use base-height augmentation to reduce the overfitting effects and improve the model's usability in new test sites.

Our experiments in high-latitude regions show that the proposed self-corrective learning can stably improve height estimation compared to baseline methods on various scenarios.

2 Problem Statement

We first introduce several basic concepts for the problem.

• Satellite footprints: Earth observation satellites collect surface data following pre-defined orbits or scan lines. The data thus center around the scan lines. Satellites with different sensors (e.g., multispectral instrument, altimeter) cover different widths along the lines, resulting in different sizes of footprints. Large or wide footprints can cover all gaps between adjacent scan lines, whereas narrow footprints leave large empty gaps.
• Large-footprint spectral imagery: Satellite imagery contains multiple spectral bands (e.g., 13 bands for the Sentinel-2A product) that can range from visible (i.e., RGB) to non-visible (e.g., near-infrared) wavelengths. Spectral imagery is often collected with large footprints.
• Large-footprint surface heights: Surface height data can also be collected with large footprints (e.g., SRTM collected by the Space Shuttle Endeavour [Farr and Kobrick, 2000]). However, this kind of data has very limited availability at large scale or has low update frequency due to the high cost (e.g., SRTM data were last collected in 2000).
• Narrow-footprint surface heights: The most recent satellite constellation ICESat-2 takes a different approach by using narrow footprints at high resolution. The satellites can also continuously fly along the orbits to provide updates throughout the many-year mission. However, the narrow footprints leave large empty gaps in-between scan lines (e.g., 17m swaths vs. 3km gaps).

Figure 1: Examples of inputs (spectral image Xi and narrowfootprint height Zi) and output (large-footprint height Yi). Areas in black color between scanning lines in Zi are empty gaps.

• Tile vs. patch: Satellite images are often provided as large tiles. For example, Sentinel-2 tiles are aligned with UTM zones, where each zone is a 6◦ × 6◦ region that can correspond to an image of size 10000 × 10000 (∼600MBs). As such tiles are often too large as training samples for deep networks, in practice they are split into small patches (e.g., 400×400, 200×200) during preprocessing, with or without overlapping.

In the problem, we are given the following three data sources for training: (1) A collection of large-footprint spectral imagery tiles that are split into small patches X = {Xi ∈ Rm×m×d } where m is the patch size and d is the number of spectral bands; (2) Narrow-footprint surface heights along scan lines that spatially and temporally intersect with Xi ∈ X. We align the height data with the image patches and preprocess them into the same image patch representations Z = {Zi ∈ Rm×m×1 }, where pixels intersecting the scan lines contain the elevation values and others are set to 0 with a mask; and (3) The large-footprint surface heights {Yi ∈ Rm×m×1 } for patches Xi ∈ X.

The prediction problem takes X and Z as inputs and aims to predict height values that are unavailable in Z due to gaps in the scan lines. Y is used as the ground-truth labels for training. Overall, Y is only available for very limited image tiles in practice, so the goal is to use the trained model to expand height coverage for tiles where height is only available as Z. Fig. 1 shows examples of the inputs and outputs.

3 RelatedWork

Pixel-level prediction with imagery: To derive height information at the pixel level, traditional methods mostly focus on stereo matching [Soergel et al., 2009], which uses two stereo images (i.e., near-simultaneous images at different angles) to reconstruct 3D information. For example, [Liu et al., 2017] estimates building heights using high-resolution stereo images in urban areas. However, stereo images are unavailable in most of the satellite platforms. Additionally, nonurban areas often do not have the high-contrast corresponding points needed by the methods. Deep semantic segmentation networks such as U-Net [Ronneberger et al., 2015], DeepLabV3+ [Chen et al., 2018] and so on can also be utilized for single-image-based height prediction. [Amirkolaee and Arefi, 2019] used a multi-scale encoder-decoder to estimate height using single aerial images. These methods perform well in urban areas. However, they rely on highresolution aerial imagery and do not have the capacity to incorporate additional information from narrow-footprint surface heights.

Traditional interpolation: Spatial interpolation methods estimate unknown values using distance-weighted aggregation of observations. Many interpolation methods have been developed, such as inverse distance weighting [Shepard, 1968], spline [McKinley and Levine, 1998], Kriging [Cressie, 2015], etc. Among the methods, Kriging [Cressie, 2015] is the most widely adopted approach, which is a nonparametric model based on Gaussian processes. Ordinary Kriging [Wackernagel and Wackernagel, 2003] performs estimation directly using known labels (or targets), whereas universal Kriging [Caballero et al., 2013] can also incorporate feature variables. However, traditional Kriging cannot model complex non-linear relationships and its performance relies on the availability of spatially-nearby observations.

Deep-learning-based interpolation: Integrations between Kriging and deep neural networks have been developed to enhance the ability to capture non-linear relationships. Kriging convolutional networks [Appleby et al., 2020] approximate Kriging interpolation using a graph neural network (GNN), which enables learning information propagation between data points. Inductive Graph Neural Network Kriging (IGNNK) model [Wu et al., 2021] trains a GNN to reconstruct information on random subgraph structures, allowing it to learn to generalize to unseen nodes/graphs for spatiotemporal interpolation. However, these methods are computationally expensive and they are intended for data with a smaller number of points (e.g., thousands), which is not suitable for satellite data that often have millions of pixels (points) per image. A recent work [Liu et al., 2022] proposed an image-based interpolation technique that parameterizes Kriging with deep neural networks to learn single-pixel and neighborhood embeddings for interpolation. However, these methods share a common limitation with traditional Kriging, i.e., their performances degrade when points with known labels are highly localized in space.

4 Spatial Self-Corrective Learning

We decompose the spatial self-corrective learning framework into four components: (1) A confidence-based pseudointerpolation method to correct low-confidence predictions in local neighborhoods; (2) A recurrent structure that uses the confidence-based pseudo-interpolation as a sub-routine and self-refines the height in an iterative manner; (3) An inverseconfidence-based regression layer that further corrects the predictions using limited observations (from narrow-footprint height data) in each recurrent step; and (4) A height-based augmentation to reduce overfitting.

4.1 Spatial Correction with Confidence-Based Pseudo-Interpolation

One challenge facing height prediction from spectral imagery is the high volume of irrelevant variations. For example, spectral features may change by surface properties that are not highly correlated with changes in height, such as types of soils or rocks, surface moisture, color differences between vegetation, etc. Such irrelevant variations in signals can be both smooth or sharp. They introduce additional confusion during the training process, and will likely lead to undesired or distorted variations in predicted heights. Since the deep network F still needs to make predictions at all locations, these hard-to-fix errors caused by the local variations may further lead to uninformative gradient updates that increase errors at other easier-to-predict locations.

To address the problems introduced by irrelevant spectral variations, we propose to reduce their impact using an "ignore-and-interpolate" correction strategy. The idea is that, rather than having the variations pollute the results, we can let the network "ignore" predictions at the locations being affected by variations and re-fill their values by doing interpolation with nearby predictions within a spatially-adjacent neighborhood S.

Denote the deep network as F, and FH and FC as two sub-branches (may have shared layers) that predict preliminary heights ˆYi and confidence scores ˆCi, respectively. The correction layer of the network updates height values at lowconfidence locations using high-confidence values to produce the output heights ˆY∗i . In the following, we explain the details of the key steps on confidence estimation and confidencebased pseudo-interpolation.

Confidence estimation: The goal of confidence estimation is to help the model locate predictions with potential high-errors. We estimate the confidence using maximum likelihood estimation (MLE). The normal-based MLE can be estimated by assuming each of the predictions (i.e., each pixel from an input image patch) follows a normal distribution, i.e, p(xj) ∼ N(ˆ yj , (1/ˆ c)2 ), where predicted heights ˆ yj ∈ ˆYi and the inverses of confidences ˆ cj ∈ ˆCi are used as means and standard deviations, respectively; and ˆYi and ˆCi are in Rm×m×1 , which have the same shape as inputs Xi as specified in the problem definition (Sec. 2). Denote ΘH and ΘC as network parameters for FH and FC, respectively, and so we have ˆY = FH(Xi,ΘH) and ˆC = FC(Xi,ΘC). Using pixels from each input image patch (xj , yj ) ∈ (Xi, Yi) as examples, the optimization is then:

where p(ˆ yj|xj) = 1 √2π/FC(xj,ΘC) · exp(- (yj-FH(xj,ΘH))2 2·(1/FC(xj,ΘC))2 ) and M = m2 is the number of pixels per image patch. The

first term shows that the errors on pixels with lower confidences (i.e., higher variances) will be scaled down, and the second term naturally constraints the variances.

Pseudo-Interpolation: We use a confidence-based interpolation layer, which is attached after the preliminary height predictions ˆYi, to generate the output heights ˆY∗i using ˆCi. We call it pseudo-interpolation as the approach we use is different from traditional interpolation where known values are used to fill in unknown values based on explicit spatial distances. In our case, all values are known but with different levels of confidences in ˆCi. Additionally, we do not use explicit distance-based interpolation (e.g., calculating distances to all points and then do a weighted average), as that process is often non-differentiable. Since the confidence-based pseudo-interpolation will be used as a sub-routine in the entire framework (Sec. 4.2), we need it to be differentiable to better co-learn confidences and heights for better final outputs. Specifically, our pseudo-interpolation uses local neighborhoods Sj around each pixel, where Sj is a W ×W window. We use 7×7 window for Sj by default. The new values are weighted averages of all values in Sj :

where JW is a W ×W matrix with all ones; ⊙ and ⊘ are Hadamard product and division, respectively; and ∗ is the convolutional operation with the latter as the kernel.

Here ˆY∗i is not the final output, and will be iteratively updated together with the confidence ˆCi in the coming section.

4.2 Recurrent Interpolation for Self-Adaptive Refinement

As after each pseudo-interpolation the confidence values should be changed at each pixel, we propose a recurrent interpolation structure to allow the predictions to continue to refine themselves based on the new heights ˆY∗i . Specifically, here we use the pseudo-interpolation as a sub-routine in a recurrent manner. After each round of executing FH, we feed the resulting ˆY∗i - together with the original features Xi - as inputs back to FH and FC for the next round of refinement. In other words, here we include ˆY∗i as an additional input to FH and FC as compared to the base version in Sec. 4.1:

• Initial inputs: Denote ˆY∗,t i as the output after the tth round. Given the new input structure, we need to provide an initial input for ˆY∗,0 i for the very first round. Here we use the narrow-footprint height Zi as ˆY∗,0 i , which is part of the known inputs to this problem (Sec. 2) and provides a peek of true heights in the patch.
• Number of recurrent steps: In test phase, the recurrent procedure can run till convergence, i.e., when:

where e is a vector of ones, and τ is a user-set tolerance. It can be constrained by a maximum number of iterations (e.g., 5, 10). For training, we fix the number of recurrent steps to 5 to reduce overhead.

• Training: We evaluate the MLE loss LMLE at each recurrent step and update the parameters after the final recurrent step using the average of the losses (Alg. 1).

4.3 Truth-Based Correction with Narrow-Footprint Heights

Although the narrow-footprint heights Zi is used as an initial seed for height expansion in the recurrent framework, the training has not yet explicitly utilized its truth-nature to more effectively correct the predictions ˆY∗i . This section aims to perform such explicit and truth-based correction with an inverse-confidence-based regression layer FR.

FR replaces the last layer of the height-prediction branch FH to perform the correction. To avoid confusion, FR becomes the new final layer of FH, which happens before the pseudo-interpolation (lines 9-10 in Alg. 1) and is included as part of the recurrent process. FR takes the features learned from its previous layer and performs a self-optimized linear combination using weights ΘR to generate the outputs ˆYi. In other words, while FR is part of the whole network, its weights ΘR are obtained by a direct least-squares type of solutions instead of by gradient descent, which has been used as a paradigm in [Bertinetto et al., 2018] for adaptation. We use this paradigm rather than standard gradient descent because: (1) It provides a direct one-step correction to allow the predictions to quickly adapt to the limited ground truth Zi in-network; (2) The least-squares solution can be computed with a closed-form solution during feed-forward, which does not interfere with or introduce major overhead to the overall gradient descent process of FH and FC; and (3) Because it can be seamlessly injected into the network, the network can learn weights to adapt to additional correction from the regression layer FR.

Note that this linear optimization is independent for each input image patch Xi, and as ΘR for FR is not updated in back-propagation, we can also consider weights in ΘR as auxiliary inputs, which are dynamically updated on-the-fly. In addition, the samples only involve the pixels intersecting with the narrow-footprints from the scan lines (Fig. 1, stripes in Zi) where ground-truth heights are known as inputs. For simplicity, denote G as the set of pixels that have narrowfootprint ground-truth given, and |G| as the set cardinality. Further, denote Hi(G) ∈ R|G|×(h+1) as a list containing h learned features from the previous layer of FR for all pixels in |G|; we add an additional "1" value to each row to learn the bias term. Similarly, denote ˆCi(G) and Zi(G) in R|G|×1 as the confidence and known-heights (inputs), respectively. The weights ΘR ∈ Rh+1 are then learned as:

where λ is the scaling factor for ridge regression, and Hi = QiRi is the QR factorization of Hi for numerical stability. We use inverse confidence to weigh pixels in the ridge regression to strengthen the correction on the potentially underperforming ones.

Finally, to make the correction more stable, we perform a random ensemble on the linear optimization step. Specifically, we split G into k overlapping subsets (e.g., 3 subsets where each contains 2/3 of points in G), and perform k separate linear optimizations of ΘR. The ensemble is then:

Alg. 1 illustrates the key steps in Sec. 4.1 to 4.3 using one batch as an example.

Algorithm 1 An iteration in spatial self-corrective learning

• Require: A batch B of {(Xi, Zi, Yi)} with M pixels per image. 1: for (Xi, Zi,Yi) in B do 2: ˆY∗,0 i = Init(Zi) {Recurrent correction. Network notation for convenience: FH = F ′H + FR (last correction layer), with ΘH = Θ′H∪ΘR; and FC has parameters ΘC.} 3: for t = 1 to T do 4: Hti = F ′H(Xi, ˆY∗,t-1 i , Θ′H) 5: ˆCti = FC(Xi, ˆY∗,t-1 i , ΘC ) 6: for j = 1 to k do 7: ΘjR = Linear Opt(Hti, ˆCti, Zi, G) {G are pixels intersecting scan lines} 8: end for 9: ˆYens,t i = 1k · ∑kj=1 FR(Hti,ΘjR) 10: ˆY∗,t i = Conf Pseudo Interpolate(ˆYens,t i , ˆCti) 11: end for 12: LiMLE = ∑Tt=1 LMLE(Yi, ˆY∗,t i , ˆCti)/(T · M) 13: end for 14: Update Θ′H and ΘC using ∑|B| i=1 LiMLE/|B|

4.4 Height Augmentation

While the spectral variations in each scene may reflect height variations, they tend to be limited in estimating the base height of the scene. More specifically, we can decompose the heights of a scene into two parts Yi = Yvar i + β, where β is a constant representing a base height. Imagine if we increase or decrease β by 100m (Fig. 2, we may obtain very similar spectral variations (e.g., shadows) as the lights from the sun can be considered very much as parallel beams with similar angles at different altitudes.

While the limited observations in Zi help, the network may still overfit to certain observed height values. Thus, we include a simple base-height augmentation, where we randomly increase or decrease the ground-truth height (identical for Yi and Zi) of different samples during training. The random change range is bounded by 500m as over-large altitude changes may no longer keep similar characteristics in some spectra. The augmentation encourages the network to focus more on predicting variations, which can then be combined with information in Zi to generate the height predictions.

5 Experiments

5.1 Datasets

Height data: As the goal is to fill height data gaps and generate large-footprint height data, ideally we also need such data for evaluation. Due to the expensive collection cost, large-footprint height information is in general less available at large scales. Among existing datasets, the most commonly used height data are provided by the Shuttle Radar Topography Mission (SRTM). Thus, we use the SRTM data as a main source for our experiments. As SRTMmainly covers regions between 60◦N and 56◦S latitudes, we additionally include the Interferometric Synthetic Aperture Radar (IFSAR) dataset, which offers elevation measurements in the Arctic regions in Alaska. Finally, we include test data from the most recent ICESat-2 satellites. Since ICESat-2 only provides narrowfootprint height observations (e.g., Fig. 1), in the testing we leave out data from a subset of the scan lines (i.e., 2 out of 6) per image patch for evaluation.

Spectral imagery: For the spectral imagery data, we need to use images that are temporally aligned with the data from SRTM, IFSAR and ICESat-2. Among them, SRTM and IFSAR were collected in 2000 and 2010, respectively, and the ICESat-2 mission started in 2018. Given the temporal distribution, we use multispectral imagery from Landsat-7, which covers all the durations. Especially, for SRTM and IFSAR, Landsat-7 was one of the few high-resolution multispectral satellite missions that overlapped with their timestamps (e.g., other popular missions such as Landsat-8 and Sentinel-2 started in 2013 and 2015, respectively).

Additional data details: As most of the bands in Landsat7 imagery have a spatial resolution at 30m (side-length of a pixel), we use 30m as the resolution in all experiments. All height data have equal or higher resolution, so they can be safely resampled to 30m. We select 5 locations in highlatitude regions with different landscapes: (1) Two locations for SRTM in the northwestern US (i.e., Washington and Idaho in Table 1), which are mountainous regions with larger height variations (e.g., steep mountainsides, deep valleys) and more forests. Each region is about 22,500 km2 in size. (2) Two locations for ISFAR in the Arctic region within northern Alaska, US (Alaska-1 and Alaska-2 in Table 1). One location has a bare-earth type of landscape with limited snow cover, and the other is mostly covered by snow. Each region is about 3,600 km2 in size as ISFAR has smaller footprints. (3) Two locations for ICESat-2 in Alaska, US (Alaska-3) and Yukon, Canada, with bare-earth landscapes. Each location is about 2,000 km2 in size, which follows along a segment of the trajectory of an ICESat-2 satellite. Finally, all the data are split and formatted into image patches of size 500 × 500 (Sec. 2).

Training and testing: We group SRTM and ISFAR together for the first evaluation as we have ground-truth largefootprint height data. We use four train-test splits, each with three locations as training and one as testing. Numbers in Tables 1 and 2 represent results with the location in the corresponding column being the test site. To simulate narrowfootprint height observations at intersecting scan lines (inputs Zi, Fig. 1), we sample 2 sets of 3-lines (ICESat-2 scans three lines simultaneously) in two different directions with equal intervals (3km; same as ICESat-2). During testing, all models are first fine-tuned with Zi before prediction. For ICESat-2, since there is no large-footprint height data, we use a trained model from ISFAR(Alaska-1) and fine-tune it with 4 out of 6 lines per scene and use the rest for testing.

5.2 Methods for Evaluation

The following baseline methods are used in the comparison:

• Base: The base network is a U-Net architecture with 3 encoding blocks and 3 decoding blocks to learn features across multiple scales and gradually combines them to achieve full-resolution height prediction.
• DA: The base network with domain adaptation by using a generative adversarial network for adversarial learning [Goodfellow et al., 2020; Tzeng et al., 2017]. DA learns domain-invariant features, but is limited when landscapes are very different in training and testing [Li et al., 2023].
• FNN: A 5-layer fully-connected network that estimates height values based on spectral features at the pixel-level.
• Kriging: The universal Kriging interpolation [Cressie, 2015] that uses spectral features and heights on the narrowfootprint surface heights along scan lines (Fig. 1(a)) as inputs to predict the unknown height at other locations. We use the closest 50 points as the neighborhood.
• Ridge: The ridge regression using L2-regularization with the scaling factor λ of 0.1. It is limited for non-linear and noisy problems [Bao et al., 2022].
• RF: Random forest regression with 100 trees using spectral features at pixel-levels.
• KCN: The Kriging convolutional network [Appleby et al., 2020] that interpolates height values based on known values on narrow-footprint surface heights along scan lines (Fig. 1(a)) for each image. It relies on Graph Convolutional

Network (GCN) to embed neighborhood information. We use the code provided by the authors in our experiments [KCN, 2019]. We keep the recommended hyperparameters and set the number of neighbors to use for interpolation to 50 (default was 5, which led to lower-quality results).

• KCN-att: KCN that uses graph attention [Veliˇckovi´c et al., 2017] instead of original GCN layers to compute the attention weights of a node's neighbors.
• KCN-sage: KCN that uses GraphSAGE [Hamilton et al., 2017] to learn aggregator functions of nodes.
• SCL: Our proposed spatial self-corrective learning method.

Metrics. We use the symmetric mean absolute percentage error (sMAPE) and correlation coefficients to evaluate the performances of the methods. We use sMAPE instead of MAPE as there are many lower-valued heights (e.g., near 0), which causes instability in MAPE calculation [Xie et al., 2023a]. sMAPE is a standard extension of MAPE, which includes predicted values in the denominators to constrain the range to [0,1].

5.3 Results

Comparative Analysis

The neural-network-based models are trained with the Adam optimizer with an initial learning rate of 10-4 . From Tables 1 and 2, we can see that the proposed spatial SCL approach outperformed the baseline methods for all the scenarios in the experiments. Domain adaptation methods did not improve upon base networks mainly because different landscapes may have different height variation patterns [Li et al., 2023]. Therefore, forcing a similar representation may hurt the ability to adapt to various landscapes. Pixel-wise regressors, such as FNN, ridge regression, and random forest perform poorly, as they can be largely affected by the spectral variations that are irrelevant to height changes. They also have limited ability to adapt to new regions. Kriging performs better than the SCL

Ground Truth

Table 3: Effects of different components in SCL.

pixel-wise regressors in most cases as it exploits the known height values in its neighborhood, which provide necessary guidance for the unlabeled locations. However, Kriging convolutional networks did not show much improvements over, and sometimes performed worse than, Kriging in this problem, as the major limiting factor of both types of methods may be the spatially-constrained distribution of known values, under which the enhancement from the additional nonlinear capacity in KCN's feature learning becomes limited.

Comparing across regions, most methods perform the best in Idaho, which has relatively sharper landscapes and visual features that make the estimation easier. In contrast, Alaska2 and Yukon were more challenging for all methods. This is potentially caused by the more flat landscapes with smaller height variations, which may reduce the amount of useful features (e.g., shadows).

Fig. 3 visualizes several examples of height estimation in several areas. We can see SCL can better capture the high mountains and deeper valleys, and at the same time keep the terrains smoother (i.e., less affected by irrelevant height variations) with the confidence-based interpolation. Fig. 4 shows detailed height profiles along randomly sampled straight lines intersecting the regions. Similarly, we can see that the predictions from SCL show the best alignment with the ground truth compared to the others.

Ablation Studies

We conducted ablation experiments to validate the effectiveness of the proposed recurrent refinement and regressionbased correction. In the first scenario, we removed both recurrent refinement and regression-based correction (No-RR) from the proposed model. The models were tested on the Seattle region and trained on the other three regions. In the second scenario, we replaced the inverse-confidence-based regression layer with a trainable 1 × 1 convolutional layer (No-Reg) for height regression. In the third scenario, we set the number of recurrent steps to 1 (No-Rec) to remove the effect of recurrent refinement while keeping other settings constant. The results in Table 3 show that the removal of either recurrent refinement or regression-based correction leads to a drop in performance.

6 Conclusions

We proposed a spatial self-corrective learning framework to expand the coverage of surface height data in high latitudes by combining narrow-footprint heights and large-footprint satellite imagery. The new framework includes several key components, i.e., confidence-based pseudo-interpolation, recurrent self-correction, and truth-based correction with a regression layer. The proposed method demonstrated consistent improvements over baseline methods in the experiments for different landscapes.

Our future work will evaluate the approach for more types of landscapes and explore the integration of the methods in domain analysis pipelines. We will also develop a benchmark dataset to facilitate future comparisons on this problem. Finally, this work has not considered the issues of cloud coverage in satellite imagery, which may need further exploration especially for the polar regions [Xie et al., 2023b].

Acknowledgments

This material is based upon work supported by the National Science Foundation under Grant No. 2105133, 2126474 and 2147195; NASA under Grant No. 80NSSC22K1164 and 80NSSC21K0314; USGS under Grant No. G21AC10207; Google's AI for Social Good Impact Scholars program; the DRI award and the Zaratan supercomputing cluster at the University of Maryland; and Pitt Momentum Funds award and CRC at the University of Pittsburgh. We also thank Dr. Sinead Farrell for the inputs on ICESat-2 data.