Mapping high-resolution basal topography of West Antarctica from radar data using non-stationary multiple-point geostatistics (MPS-BedMappingV1)

The subglacial bed topography is critical for modelling the evolution of Thwaites Glacier in the Amundsen Sea Embayment (ASE), where rapid ice loss threatens the stability of the West Antarctic Ice Sheet. However, mapping of subglacial topography is subject to uncertainties of up to hundreds of metres, primarily due to large gaps of up to tens of kilometres in airborne ice-penetrating radar flight lines. Deterministic interpolation approaches do not reflect such spatial uncertainty. While traditional geostatistical simulations can model such uncertainty, they become difficult to apply because of the significant non-stationary spatial variation of topography over such large surface area. In this study, we develop a non-stationary multiple-point geostatistical (MPS) approach to interpolate large areas with irregular geophysical data and apply it to model the spatial uncertainty of entire ASE basal topography. We collect 166 high-quality topographic training images (TIs) of resolution 500 m to train the gap-filling of radar data gaps, thereby simulating realistic topography maps. The TIs are extensively sampled from deglaciated regions in the Arctic as well as Antarctica. To address the non-stationarity in topographic modelling, we introduce a Bayesian framework that models the posterior distribution of non-stationary TIs assigned to the local line data. Sampling from this distribution then provides candidate training images for local topographic modelling with uncertainty, constrained to radar flight line data. Compared to traditional MPS approaches that do not consider uncertain TI sampling, our approach results in a significant improvement in the topographic modelling quality and efficiency of the simulation algorithm. Finally, we simulate multiple realizations of high-resolution ASE topographic maps. We use the multiple realizations to investigate the impact of basal topography uncertainty on subglacial hydrological flow patterns.

1 Introduction

The topography beneath the Greenland and Antarctic ice sheets is essential for nearly every ice sheet investigation, including modelling subglacial hydrology (e.g. De Fleurian et al., 2018; MacKie et al., 2021a; Siegert et al., 2016; Sommers et al., 2018), interpreting geologic conditions (Bingham and Siegert, 2009; King et al., 2009; Rippin et al., 2014; Holschuh et al., 2020; Alley et al., 2021), estimating ice volume and sea level rise contributions (e.g. Fretwell et al., 2013; Morlighem et al., 2020) and ice sheet modelling for sea level rise projections (e.g. Le clec'h et al., 2019; Schlegel et al., 2018; Seroussi et al., 2017). The characterization of subglacial topography is particularly important for Thwaites Glacier in the Amundsen Sea Embayment, which is experiencing accelerating ice loss (Rignot et al., 2019) that could destabilize the West Antarctic Ice Sheet (Joughin et al., 2014). Subglacial topography is predominantly measured with airborne ice-penetrating radar along flight lines separated by up to tens of kilometres (e.g. Bingham et al., 2017; Fretwell et al., 2013; Herzfeld et al., 1993). Large gaps in data remain, which are generally interpolated deterministically using methods such as kriging (Herzfeld et al., 1993), the ArcGIS Topogrid algorithm (Fretwell et al., 2013), spline interpolation (Lytheand Vaughan, 2001; Holt et al., 2006) or ice sheet model inversions (Farinotti et al., 2017; Huss and Farinotti, 2012; Morlighem et al., 2017, 2020). These approaches produce topography that is unrealistically smooth and provide limited morphological information. Furthermore, deterministically interpolated topography does not sample the uncertainty space, making it difficult to quantify uncertainty in ice sheet models with respect to topographic uncertainty. These issues have previously been addressed with two-point geostatistical simulation, such as fast Fourier transforms (Goff et al., 2014; Graham et al., 2017; MacKie et al., 2020) or sequential Gaussian simulation (SGSIM) (MacKie et al., 2021a). The objective of geostatistical simulation is to generate multiple realizations of phenomena that reproduce the spatial variability of observations, as modelled by variogram or spatial covariance and can be used to quantify uncertainty (e.g. Deutsch and Journel, 1998). Goff et al. (2014) also conducted a conditional simulation of Thwaites Glacier. To improve the modelling quality, the channelized structures and the abrupt transition between lowland and highland are individually handled. The method has the advantage to ensure the continuity of fjord-like channels beneath the glacier. Geostatistical simulation has also been applied in Antarctica and Greenland to quantify uncertainty in subglacial hydrology (MacKie et al., 2020, 2021a; Zuo et al., 2020).

However, spatial variation over very large areas is inherently non-stationary. For example, the Greenland and Antarctic ice sheets are thousands of kilometres in length and contain a wide range of topographic and geologic settings. This means that the nature of spatial variation changes significantly and possibly in complex ways over the domain of interest. Traditional geostatistical ways of dealing with non-stationary data are through the modelling of trend functions (e.g. Pyrcz and White, 2015) or using covariates (e.g. Almeida and Journel, 1994; MacKie et al., 2021a). However, such approaches typically model the variation in the mean (trend) or some degree of correlation (co-simulation). Another approach is using a non-stationary spatial covariance model (Schmidt and O'Hagan, 2003). Such an approach becomes exceedingly difficult to apply over large areas because of the use of Markov chain Monte Carlo in its Bayesian inference. Regardless, most two-point geostatistical approaches are limited in expressing non-stationarity in terms of a mean or covariance function only.

The non-stationary bed topography is measured using high-resolution remote sensing data such as satellite imagery but only in deglaciated areas (Porter et al., 2018). Deglaciated topographic images reveal glaciated morphologies resembling the topography beneath the contemporary ice sheets (King et al., 2009; Margold et al., 2015; Spagnolo et al., 2017). They therefore bear significant morphological information on the subglacial topography. Exposed topography has previously been used to perform deterministic interpolations (Clarke et al., 2009). However, deglaciated topography has not been used to stochastically simulate subglacial topography, until the very recent alpine glacier study by Neven et al. (2021) using multiple-point geostatistics.

Recent developments in multiple-point geostatistics (MPS) has shown great potential in using high-resolution training images (e.g. satellite images) to fill remote sensing gaps (e.g. Gravey and Mariethoz, 2020; Mariethoz et al., 2012; Yin et al., 2017; Zakeri and Mariethoz, 2021; Zuo et al., 2020). MPS approaches use the training images (TIs) as explicit prior models to generate realistic topographical models and quantify spatial uncertainty. The simulation of non-stationary and morphologically complex topography can also be achieved with MPS (Hoffimann et al., 2017, 2019; Mariethoz and Caers, 2014). Compared to alternative machine learning or deep learning approaches (Laloy et al., 2018; Mo et al., 2020), MPS has a flexible conditioning capability and can accommodate sparse and non-uniform sampling in space. It can generate multiple topographic model realizations conditioned to the radar line observations, without requiring a large amount of training data.

We briefly review three categories of approaches to build non-stationary geospatial models using MPS. The first way is to divide non-stationary TI or simulation grid into several stationary subareas. Each stationary simulation area has its specified stationary TI (Honarkhah and Caers, 2012; Strebelle, 2002; Wu et al., 2008; Zhou et al., 2014). However, the zonation may create artefacts in areas where two zones meet. Therefore, the second way is most commonly used. It incorporates spatially continuous non-stationary maps (named as “auxiliary variables”) with weighted aggregation or so-called “ad-hoc weighting” (Chugunova and Hu, 2008; Mariethoz et al., 2010; Oriani et al., 2014; Zuo et al., 2020). Such auxiliary variables determine which TI patterns should fill which location in the simulation domain in a spatially smooth manner. The limitation is that the ad-hoc weights do not scale to the complexity of bed topography. The determination of weights is also subjective. Moreover, auxiliary variables are very difficult to obtain in subglacial topographic modelling. Another challenge in the non-stationary modelling is how to choose training images (Tahmasebi, 2018). This is particularly important as the MPS modelling relies on the spatial information provided by the training images. In the third option, Hoffimann et al. (2019) introduced an approach to generate time-series training images to model the spatial and temporal evolutions of geomorphology, which is similar to Pirot et al. (2014, 2015). A training image transitional model in time was proposed to reproduce the non-stationary geomorphologic evolutions. However, in subglacial topographic modelling, there is limited training imagery because subglacial topographic measurements are only made along flight lines. Graham et al. (2017) provided a synthetic bed elevation terrain for Antarctica. Satellite-based observations from deglaciated areas in the Arctic offer a potential source of training imagery. But these Arctic training images would be non-stationary due to the natural variability of the landscape. Furthermore, the Arctic provides a vast amount of deglaciated topographic data, which present a significant computational burden on MPS simulation algorithms. We therefore need a strategy to explicitly specify which training images or patterns should go where when filling the radar line gaps.

In this paper, we generalize a geospatial modelling framework to fill irregular geophysical data gaps in large areas. We use MPS to address the non-stationary topographic modelling by probabilistically selecting non-stationary training images. We first collect a large number of TIs from the deglaciated areas in the Arctic and Antarctica. Then we develop a probability aggregation method to estimate each TI's probability of being assigned to different local radar lines. This probabilistic TI selection scheme avoids the use of auxiliary variables with arbitrary ad-hoc weighting. We demonstrate our method using the entire Amundsen Sea Embayment (ASE) in West Antarctica. This region has alternating areas of sparse and dense measurements with a variety of radar line orientations. We show that the training image sampling process accommodates a range of data configurations. We generate realistic non-stationary topographic realizations that reflect the subglacial topographic uncertainty in ASE. We then use the topographic simulations to model subglacial hydrologic flow in order to investigate the impact of topographic uncertainty on hydrologic uncertainty.

2 Radar data set and training images

The topographic data for the ASE include seafloor bathymetry measurements from the International Bathymetric Chart of the Southern Ocean (IBCSO) (Arndt et al., 2013), subaerial topography from the Reference Elevation Model of Antarctica (REMA) (Howat et al., 2019), and ice-penetrating radar measurements of subglacial topography (Blankenship et al., 2001; Gogineni, 2012; Holschuh et al., 2020; Holt et al., 2006; Vaughan et al., 2006; Young et al., 2016). The data are gridded at a 500 m spatial resolution by averaging the measurements within each grid cell (Fig. 1). The swath bathymetry data (Arndt et al., 2013) and subglacial swath radar data (Holschuh et al., 2020) provide some training imagery. To have more extensive representations of the subglacial topography, we increase the available training data with deglaciated subaerial topography from ArcticDEM V3 (Porter et al., 2018). The Arctic and much of North America were formerly covered by the Laurentide and Cordilleran ice sheets and share morphological similarities with Antarctic subglacial topography (King et al., 2009; Margold et al., 2015). While the seafloor and subaerial topography may have experienced additional erosional and depositional processes after deglaciation, any topographic alterations are likely minimal at a 500 m resolution. We sampled a total 166 candidate training images to capture a variety of geological settings (Fig. 2). Each training image has a size of 100 km × 100 km. The training image data repository is publicly accessible on the Zenodo repository (https://zenodo.org/record/5083715#.YQT2JI5Kiiw (last access: 8 July 2021), https://doi.org/10.5281/zenodo.5083715, MacKie et al., 2021b).

Figure 1 Radar line surveys of the Thwaites and Pine Island glaciers in the Amundsen Sea Embayment of West Antarctica. The red box implies the location of the study area in Antarctica. Black lines indicate boundaries for Thwaites Glacier, ice shelves, and the grounding line (the point where the ice detaches from the bed and achieves flotation). The two topography patches in the centre of Thwaites Glacier were measured using swath radar (Holschuh et al., 2020).

Figure 2 Geographical locations of the 166 training images in (a) ArcticDEM (red boxes encompass TI regions) and (b) Antarctica's Amundsen Sea Embayment region. (c) Examples of the 166 training images.

3 Methodology

3.1 Multiple-point geostatistics

3.1.1 Overview

Multiple-point geostatistics (Journel and Zhang, 2006; Mariethoz and Caers, 2014; Srivastava, 2018; Strebelle, 2002) is the field of study that focuses on the digital representation of physical reality by reproducing high-order statistics inferred from training images. The emphasis in MPS lies on capturing higher order (hence multi-point statistics) from training images that have been selected to be representative for a specific area of study. In that sense, it differs from spatial covariance-based (variograms) methods (e.g. Gaussian process regression or kriging) (Matheron, 1963; Williams and Rasmussen, 1996) that are based on spatial correlation (two values at a time). Both MPS and covariance-based methods have the ability to interpolate data exactly at the data locations. Exact interpolation, if desired, is also where geostatistics differs from machine learning or computer vision methods, where such exact interpolation is not usually considered important.

Several MPS simulation algorithms (e.g. Gravey and Mariethoz, 2020; Hoffimann et al., 2017; Mariethoz et al., 2010; Strebelle and Journel, 2001) have been developed that use training images to generate multiple realizations that interpolate the data exactly. The algorithm used in this work is direct sampling (DS) (Mariethoz et al., 2010; Mariethoz and Renard, 2010), which will be introduced in the following section. These algorithms do not address the challenge of selecting the training images themselves. For example, if an area of the simulation grid contains dense data, few training images may be compatible with those data. On the other hand, an area with sparse data may have many compatible training images. Finally, training images selected for two adjacent areas are not necessarily independent from each other.

To date, there has not been any attempt to use MPS to interpolate ice-penetrating radar measurements of topography at the scale of the Amundsen Sea Embayment. In doing so, additional challenges occur that are not present in smaller study areas. The challenges may include limited number of training images, non-stationarity over the ASE, and running time cost when generating high-resolution topographic maps. Before moving to the methodology that addresses these challenges, we first introduce the direct sampling method and a probabilistic framework for representing training images in metric spaces.

3.1.2 Direct sampling

Direct sampling is a widely used MPS approach for achieving spatial modelling and gap filling (Mariethoz, 2018; Mariethoz et al., 2012; Zuo et al., 2020). Figure 3 provides a simplified example of DS in the context of flight lines. The values in the grid indicate the elevation. In general, there are two major components within DS. First, the algorithm visits an unknown location in the simulation grid and collects neighbouring observed points as conditioning data. For example, in Fig. 3a, three conditioning points are detected near an unknown location (marked with “?”). DS records the values and relative locations of known points. Second, a searching programme is launched to find the similar structure in TI. The similarity within DS is defined by a certain distance metric (e.g. Hamming distance for categorical variable and Euclidean distance for continuous variable). As Fig. 3d shows, the programme finds a matching structure. The centre of the similar instance is pasted into the simulation grid. Thus, the value of an unknown point is predicted. The preceding simulation programme is repeatedly performed until there is no unknown point in the grid.

Figure 3 Conceptual example of the DS point simulation. (a) Radar lines on the simulation grid; (b) three known points (values: 37, 80, 86) constitute a conditioning data pattern; (c) a mismatch pattern in TI; (d) a similar pattern in TI.

Based on the explanation above, there are three main key parameters within DS. (1) The first parameter is the number of conditioning points n. In a continuous simulation scenario, n≥30 is suggested to extract complex patterns from TI as well as the simulation grid (Bruna et al., 2019; Meerschman et al., 2013). (2) Then there is the distance threshold t. It is possible that there is no completely matching structure in the TI. Therefore, the algorithm accepts a training pattern whose distance with the conditioning pattern is lower than t. When many suitable patterns exist in TI, the first pattern found will be accepted. The value of t has a significant influence on the DS performance. A small value will improve modelling quality but may result in a significant computational cost. In most cases, t=0.1 is generally recognized as the upper bound (Meerschman et al., 2013; Zuo et al., 2020). (3) Lastly, there is the fraction of scanned TI f. Repeated morphological structures can be common in TI. With the aim of saving time, we can scan only a fraction of TI. For example, f=0.1 implies that the computer only inspects 10 % TI. According to the investigation conducted by Mariethoz and Caers (2014), a recommended value of f ranges from 0.1 to 0.5.

3.1.3 A metric space for training images

A metric space expresses the relationship between objects by using a distance function defining the similarity between any two objects. In metric spaces, we do not know the exact coordinates of objects but only how far objects are apart. Metric spaces are therefore useful in representing high-dimensional objects, such as training images. In this paper, we employ metric spaces for two purposes: (1) to visualize the difference between training images and (2) to estimate probabilities of training images to occur over some areas.

To define a meaningful distance between any two training images, we create a set of representative patterns for each training image. The TI morphological features are mainly concerned when creating the representative patterns. This requires first removing the effects of the original TI elevations. To do so, we rescale each TI to a range between 0 and 1 by min–max normalization (Han et al., 2012). Then, like other MPS approaches (Honarkhah and Caers, 2010; Strebelle, 2002), we extract all the spatial patterns from each TI using a fixed template. Next, we use agglomerative hierarchical clustering (Romary et al., 2015) to divide the spatial patterns of each TI into a finite number of groups. The number of groups is determined by a distance threshold between the clustered groups in agglomerative hierarchical clustering. As mentioned in Sect. 3.1.2, we set the distance threshold as 0.1 since it is commonly used to distinguish two patterns in DS (Meerschman et al., 2013). Therefore, a TI with more complex spatial patterns will have more clusters and thus more representative patterns. The medoid pattern of each group is taken as the representative pattern of the TI. Figure 4 shows a few representative patterns. The distance used in the clustering is the normalized Euclidean distance between the pixel-wise evaluations.

Figure 4 Calculating the distance between any two training images (TI_A and TI_B) using modified Hausdorff distance. There are three key steps. (1) The first is to extract training patterns with a fixed template. (2) The representatives are selected by a hierarchical clustering method. In this example, the method found 16 important patterns from TI_A, and 21 patterns are from TI_B. The number of representatives is dependent on the complexity of morphology. (3) The third step is to calculate the modified Hausdorff distance between two pattern sets. The output distance becomes an indicator of similarity between two TIs.

After clustering and medoid selection, TIs can now be expressed by a set of representative patterns. We define the difference between any two training images as the difference between their sets of representative patterns. To do this, we use the modified Hausdorff distance (Dubuisson and Jain, 1994; Huttenlocher et al., 1993). This distance is commonly used to quantify the difference between shapes of high-dimensional objects. In detail, if we call the set of representative patterns for training image A as 𝒜 and for training image B as ℬ, then the modified Hausdorff distance is

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{dis}\left({\mathrm{TI}}_A,{\mathrm{TI}}_B\right)& =max\left\{{\displaystyle \frac{\mathrm{1}}{\left|\mathcal{A}\right|}}\sum _{x_{\mathcal{A}}\in \mathcal{A}}mind(x_{\mathcal{A}},\mathcal{B}),\right.\\ & \left.{\displaystyle \frac{\mathrm{1}}{\left|\mathcal{B}\right|}}\sum _{y_{\mathcal{B}}\in \mathcal{B}}mind(y_{\mathcal{B}},\mathcal{A})\right\},\end{array}\end{array}$$

where x_𝒜 any one of the representative patterns in 𝒜, and y_ℬ is any pattern in ℬ. d is the Euclidean distance between any two representative patterns. |A| and |B| are respectively the sizes of 𝒜 and ℬ. In essence, the modified Hausdorff distance represents the maximum of expected minimum distances between the two TIs' representative patterns. Once a distance is defined, we can visualize the metric space in low-dimensional Cartesian space using multi-dimensional scaling or MDS (Scheidt et al., 2018). The main idea of MDS is to project objects from a high-dimensional space into a 2-D Cartesian space to visualize the similarity between all TIs. Figure 5 shows the projection of 166 training images in 2-D; each dot represents a TI. Similar training images map close to each other in the MDS scatterplot.

Figure 5 Visualization of the metric space using multi-dimensional scaling (MDS) into a two-dimensional Cartesian space. Each dot on the plot represents a TI. It shows TIs with similar morphology are close in this metric space.

3.2 Illustration case and overview of the mapping strategy

To illustrate the proposed methodology, we focus on a small area of the ASE overlapping Pine Island Glacier (see Fig. 6). In this area, we observe a variety of radar line geometries and densities, as well as elevation changes. This smaller area is divided into four subareas. Strategies for such subdivision will be discussed later in the application to the entire domain. Direct sampling, by construction, avoids any artefact boundary between the radar line subareas, because the data template is not limited by subareas borders. With this strategy, two problems now need to be addressed. First, we need to find training images that are consistent with radar data within a selected subarea. There could be multiple such training images. Second, we need to model training image cross-correlation between different subareas. Training images of two adjacent subareas are not necessarily independent because of spatial correlations between the subareas. Our approach is to model the posterior TI distribution of each area through a probability aggregation problem.

Figure 6 A subset of the pine island glacier is used to illustrate the methodology. The red line divided four subareas: A₁, A₂, A₃, and A₄. The black coloured areas highlight line gaps of the selected subset.

Figure 7 Illustration of measuring the distance between training image and radar line data (d) in subarea A₁. We first extract a group of radar data patterns D from the simulation grid with flexible-sized templates. Then the Hausdorff distances between the representative patterns A and radar patterns D are individually computed. Representative pattern x_A has a fixed size of 23×23 pixels, while the size of conditioning data pattern y_D varies.

3.3 Formulation of the problem through probability aggregation

Our goal is to estimate, for each area A₁,…A₄, the posterior distribution of training images, given the flight radar line data $d_{A_i}$. This can be formulated as

$$\begin{array}{}\text{(2)}& P\left(\mathrm{TI}\left(A_i\right)\mathrm{|}{d}_{A_{\mathrm{1}}},d_{A_{\mathrm{2}}},d_{A_{\mathrm{3}}},d_{A_{\mathrm{4}}}\right).\end{array}$$

TI(A_i) is a discrete random variable that has 166 possible outcomes (number of candidate TIs). To obtain the posterior distribution, we first estimate individual conditional probability $P\left(\mathrm{TI}\right(A_i)\mathrm{|}{d}_{A_j}$), then aggregate them into a single estimate for Eq. (2). We will use a simple aggregation model that uses log ratios (Allard et al., 2012), as follows:

$$\begin{array}{}\text{(3)}& r_{ij}=\log \left({\displaystyle \frac{P\left(\mathrm{TI}\right(A_i\left)\mathrm{|}{d}_{A_j}\right)}{\mathrm{1}-P\left(\mathrm{TI}\right(A_i\left)\mathrm{|}{d}_{A_j}\right)}}\right).\end{array}$$

Here, i and j are indices of the subarea, respectively. To aggregate these individual conditional probabilities, the log ratios can be summed relative to the prior:

$$\begin{array}{}\text{(4)}& r_i-r={\sum }_{j=\mathrm{1}}^{\mathrm{4}}\left(r_{ij}-r\right).\end{array}$$

Here, r_i is the log ratio of $P\left(\mathrm{TI}\left(A_i\right)\mathrm{|}{d}_{A_{\mathrm{1}}},d_{A_{\mathrm{2}}},d_{A_{\mathrm{3}}},d_{A_{\mathrm{4}}}\right)$ in Eq. (2). r is the log ratio of the prior. The prior is a uniform distribution over all training images. Thus r is calculated as

$$\begin{array}{}\text{(5)}& r=\log \left({\displaystyle \frac{P\left(\mathrm{TI}\right(A_i\left)\right)}{\mathrm{1}-P\left(\mathrm{TI}\right(A_i)}}\right)=\log \left({\displaystyle \frac{\mathrm{1}/\mathrm{166}}{\mathrm{1}-\mathrm{1}/\mathrm{166}}}\right)=-\mathrm{5.10}.\end{array}$$

Then, we can solve Eq. (4) for r_i and invert the log ratio to get $P\left(\mathrm{TI}\right(A_i)\mathrm{|}{d}_{A_{\mathrm{1}}},d_{A_{\mathrm{2}}},d_{A_{\mathrm{3}}},d_{A_{\mathrm{4}}}$).

However, in summing, we make a conditional independence assumption (Allard et al., 2012). Indeed, summing logarithms is equivalent to making products of the actual probabilities, which entails a form of conditional independence. Assuming conditional independence, when that assumption is untrue in reality, often results in overconfidence and too-small uncertainty. To mitigate this issue, we add an additional weight term w_ij:

$$\begin{array}{}\text{(6)}& r_i-r={\sum }_{j=\mathrm{1}}^{\mathrm{4}}{w}_{ij}\left(r_{ij}-r\right).\end{array}$$

Logically, we would like the weight to account for the correlation between data in different regions. For example, if data of region A_i are highly correlated with data in region A_j, then they are likely redundant with respect to the training image selection. Hence, we will make the weight w_ij function of the correlation structure between different subareas. In the next section, we will detail the subtasks ahead: (1) modelling and estimating $P\left(\mathrm{TI}\right(A_i)\mathrm{|}{d}_{A_j}$) and (2) calculating the weights w_ij.

3.4 Probability of training images given radar line data

3.4.1 Most probable set of training images

A direct estimate of $P\left(\mathrm{TI}\right(A_i)\mathrm{|}{d}_{A_i}$) is challenging because the line data $d_{A_i}$ are very high dimensional. For example, there are 7982 radar measurements in subarea A2. We turn this high-dimensional problem into a low-dimensional space as follows. Using the data $d_{A_i}$ in area A_i, we find those training images that constitute a set of most probable training images, i.e. those images closest to the radar lines in that area in terms of morphological similarities. We term this set as $\widehat{\mathrm{TI}}$. Then given this set, we replace the radar line data with the most probable set:

$$\begin{array}{}\text{(7)}& P\left(\mathrm{TI}\right(A_i\left)\mathrm{|}{d}_{A_i}\right)\cong P\left(\mathrm{TI}\right(A_i\left)\mathrm{|}\widehat{\mathrm{TI}}\right).\end{array}$$

To determine this set, we solve the following optimization problem:

$$\begin{array}{}\text{(8)}& \underset{\widehat{\mathrm{TI}}}{\mathrm{argmin}}\left\{\mathrm{dis}\left({\mathbb{I}}_{\mathit{T}\mathit{I}}\left(\widehat{\mathrm{TI}}\right),d_{A_i}\right)\right\},\end{array}$$

where 𝕀_TI is an indicator function which returns $\widehat{\mathrm{TI}}$, a n-size subset of TI. $\mathit{T}\mathit{I}=\left[{\mathrm{TI}}^{\left(\mathrm{1}\right)},{\mathrm{TI}}^{\left(\mathrm{2}\right)},\mathrm{\dots },{\mathrm{TI}}^{\left(\mathrm{166}\right)}\right]$, which is the total set of training images. We will explain how to determine the size n of $\widehat{\mathrm{TI}}$ via particle swarm optimization (see the Appendix). The distance (dis) in Eq. (8) measures the distance between the radar line data and any given training image. To calculate dis, we rely on the same modified Hausdorff distance approach as in Sect. 3.1.3:

$$\begin{array}{}\text{(9)}& \mathrm{dis}\left(\mathrm{TI},d\right)={\displaystyle \frac{\mathrm{1}}{\left|\mathcal{D}\right|}}\sum _{y_D\in \mathcal{D}}mind(y_{\mathcal{D}},\mathcal{A}).\end{array}$$

𝒟 is now the set of patterns y_𝒟 extracted from the radar lines. By pattern, we mean the radar line data are scanned within a given template. We use flexible-sized templates when scanning the radar lines over each subarea A_i. The template size varies in order to collect 40 neighbouring radar data points but with a maximum radius up to 15 pixels. 𝒜 is the set of TI representative patterns x_𝒜, obtained in Sect. 3.1.3. Figure 7 provides an illustration of this idea.

Figure 8 (a) The estimated set of most probable training images $\widehat{\mathrm{TI}}$ on MDS plots for each subarea (A1, A2, A3, A4). The red dots highlight the estimated $\widehat{\mathrm{TI}}$. (b) Examples of $\widehat{\mathrm{TI}}$ displayed in the topographic modelling space.

We use a particle swarm optimization (PSO) to minimize the distance function $\mathrm{dis}\left({\mathbb{I}}_{\mathit{T}\mathit{I}}\left(\widehat{\mathrm{TI}}\right),d_{A_i}\right)$ (see the Appendix). As a heuristic optimization approach, PSO has its specific advantages in requiring less parameterization, easy implementation, and fast convergence with good accuracy (Rezaee Jordehi and Jasni, 2013; Sengupta et al., 2019). These characteristics make PSO a preferred optimizer for our initial training image selection. In conjunction with PSO, we employ the profile log-likelihood function to find the optimal size n of $\widehat{\mathrm{TI}}$. A detailed explanation of the PSO algorithm and profile log-likelihood function implementations is provided in the Appendix. Figure 8 shows the selected $\widehat{\mathrm{TI}}$ in metric space for each subarea (A1, A2, A3, and A4). In this figure, we also plot examples of $\widehat{\mathrm{TI}}$ in the radar line map grid.

3.4.2 Kernel density estimation of $P\left(\mathrm{TI}\right(A_i)\mathrm{|}{d}_{A_i}$)

We will use the optimal set of training images $\widehat{\mathrm{TI}}$ to infer $P\left(\mathrm{TI}\right(A_i)\mathrm{|}{d}_{A_i}$). We assume that TIs near the $\widehat{\mathrm{TI}}$ on the MDS plot tend to have similarly high probability of being assigned to the radar data subarea. This is because nearby TIs in the MDS metric space (see Fig. 5) show similar morphological patterns. We therefore consider a Gaussian kernel density estimation (KDE) to predict the probability of TI being assigned to a subarea A_i. The probability of each TI is estimated according to its distance with $\widehat{\mathrm{TI}}$ in the MDS plot (Fig. 5):

$$\begin{array}{}\text{(10)}& {\displaystyle }P\left(\mathrm{TI}\left(A_i\right)\mathrm{|}{d}_{A_i}\right)& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{n}}\sum _{k=\mathrm{1}}^{n}K\left({\displaystyle \frac{\mathrm{dis}(\mathrm{TI},{\widehat{\mathrm{TI}}}_k)}{h}}\right)\text{(11)}& {\displaystyle }K\left(z\right)& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}}}\exp \left(-{\displaystyle \frac{z^{\mathrm{2}}}{\mathrm{2}}}\right).\end{array}$$

Here, ${\widehat{\mathrm{TI}}}_k$ is the kth selected TI using PSO. n is the size of the set $\widehat{\mathrm{TI}}$. $\mathrm{dis}(\mathrm{TI},{\widehat{\mathrm{TI}}}_k$) is modified Hausdorff distance between a TI and ${\widehat{\mathrm{TI}}}_k$. K is the Gaussian kernel function (Eq. 11). The bandwidth h is the variance of the Gaussian kernel. We calculate the optimal bandwidth h by following Silverman's rule of thumb (Silverman, 1981). Figure 9 shows the KDE-estimated probability of each TI for subarea A1.

Figure 9 (a) Estimated $\widehat{\mathrm{TI}}$ for subarea A1 in MDS space. The red dots are the $\widehat{\mathrm{TI}}$, while blue points represent other TIs. (b) Kernel density smoothing assigns likelihoods (densities) to the total set of training images by using $\widehat{\mathrm{TI}}$.

Figure 10 Probability distribution of final aggregated TI probability in each radar line subarea.

3.5 Aggregation by weighting log ratios

Next, we aggregate the KDE-estimated probabilities by weighting the log ratios to obtain the posterior $P\left(\mathrm{TI}\left(A_i\right)\mathrm{|}{d}_{A_{\mathrm{1}}},d_{A_{\mathrm{2}}},d_{A_{\mathrm{3}}},d_{A_{\mathrm{4}}}\right)$. The weights w_ij required in the log-ratio aggregation of Eq. (6) are used to quantify the spatial correlation between the radar line subareas. We use a variogram-based approach proposed by Fouedjio (2020) to measure spatial correlation between any two areas. In detail, the variogram dissimilarity is calculated as the sum of absolute values of all direct and cross variograms between the two subareas.

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& \mathrm{dissim}\left(A_i,A_j\right)=\\ & {\displaystyle \frac{{\sum }_{\mathrm{l},{\mathrm{l}}^{\prime }}{K}_{\mathrm{E}}\left(\left({\mathit{x}}_i,{\mathit{x}}_j\right),\left(x_{\mathrm{l}},x_{{\mathrm{l}}^{\prime }}\right)\right)\times {\left(z\left(x_{\mathrm{l}}\right)-z\left(x_{{\mathrm{l}}^{\prime }}\right)\right)}^{\mathrm{2}}}{\mathrm{2}{\sum }_{\mathrm{l},{\mathrm{l}}^{\prime }}{K}_{\mathrm{E}}\left(\left({\mathit{x}}_i,{\mathit{x}}_j\right),\left(x_{\mathrm{l}},x_{{\mathrm{l}}^{\prime }}\right)\right)}},\end{array}\end{array}$$

where dissim(A_i,A_j) is the dissimilarity between subareas A_i and A_j. According to Fouedjio (2020), (x_i,x_j) are the spatial centre locations of A_i and A_j, respectively. x_l and $x_{{\mathrm{l}}^{\prime }}$ are the radar data point locations in A_i and A_j. K_E(⋅) is the Epanechnikov kernel function (Fouedjio, 2020). z(x) are the radar data measured values at location x. Using the calculated dissim(A_i,A_j), the weights w_ij are

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& w_{ij}=\mathrm{1}-\\ & {\displaystyle \frac{\mathrm{dissim}\left(A_i,A_j\right)}{max\left(\mathrm{dissim}\left(A_i,A_{\mathrm{1}}\right),\mathrm{\dots },\mathrm{dissim}\left(A_i,A_j\right),\mathrm{\dots },\mathrm{dissim}\left(A_i,A_{\mathrm{4}}\right)\right)}},\\ & \mathrm{where}ij=[\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}].\end{array}\end{array}$$

With w_ij, we can aggregate the probability using the Eq. (4). Figure 10 shows the aggregated posterior probability of TIs for each subarea.

Figure 11 Examples of sampled TIs from the posterior distributions corresponding to the subarea (A1, A2, A3, A4). The TIs are rescaled back to the local radar data range by inverting the min–max normalization.

Figure 12 (a) Two topographical model realizations from using our proposed DS with uncertain TI sampling to fill the radar line gaps. The model realization number corresponds to the TI realization number in Fig. 11. (b) Line gap filling by traditional DS using all the 166 TIs (without TI sampling). Panels (c) and (d) show line gap filling using kriging and SGSIM.

3.6 Direct sampling with uncertain TI sampling

Using the aggregated posterior TI probability $P\left(\mathrm{TI}\right(A_i)\mathrm{|}{d}_{A_{\mathrm{1}}},d_{A_{\mathrm{2}}},d_{A_{\mathrm{3}}},d_{A_{\mathrm{4}}}$), we can now sample training images from the posterior distribution (Fig. 10) in each subarea. Figure 11 plots two realizations of sampled training images on the radar line map. We observe that the sampled topography TIs are different between the realizations. For example, the TIs sampled for A1 tend to have higher elevations and more mountain peaks than A2. A2 and A4 tend to have larger-scale low-elevation valleys, while the TIs in A1 and A3 data have more small-scale valleys. For each realization of training images, we run a DS simulation. Finally, multiple realizations of topographical models are generated, each with multiple realizations of TIs. Figure 12a shows two example realizations of the DS simulated results. We can observe large valleys in A2 and A4 areas, while A1 and A3 areas mainly have high elevation peaks. The non-stationarity of both simulated topography realizations also agreed well with their sampled TIs when compared to the TIs in Fig. 11.

3.7 Comparison with traditional MPS modelling and two-point geostatistical modelling

Our results are compared to the conventional MPS simulation without proposed TI sampling. Here, we use the same DS simulation parameters as our TI sampling approach, except that the training images are different. In the conventional test, we run the DS simulation by scanning all 166 TIs to fill the radar line gaps. Figure 12b shows one realization created by the conventional DS. It is obvious that the conventional approach results in a much noisier topographical model, and there are significant line artefacts that make the model unrealistic. To gain a detailed understanding, we take a cross section A–A^′ across the trunk of Pine Island Glacier and plot the gap-filling result comparison in Fig. 13. We can observe that the DS without TI sampling creates more unrealistic elevation peaks and troughs. In particular, at the main channel of Pine Island Glacier (marked by the dashed box in Fig. 13) where dense radar data are available, unrealistic channels are simulated between the radar data observations. This suggests that, when using all the 166 TIs without proper sampling, the DS finds too-large a set of patterns of which many are incompatible with the sparse data. Our TI sampling approach avoids this problem by limiting the algorithm to a small number of most suitable TIs, thereby improving the result (see Fig. 12a). More importantly, avoiding the channel artefacts is critically important for modelling subglacial hydrological flow (see Sect. 4.2). In terms of running time, the conventional DS approach with 166 TIs took nearly 21 h to simulate one realization. By contrast, our TI sampling approach took less than 1 h. Our initial DS implementation tests are run on a PC with an Intel i9-11900 2.5 GHz processor and 32 GB of RAM.

Figure 13 Cross-section view of the modelled topography maps at line A–A^′. The dashed black box shows main channel of Pine Island Glacier.

Figure 14 Comparison of SGSIM and the proposed method of DS with uncertain TI sampling in a local sparse line area. Red ovals highlight the areas where SGSIM does not to simulate morphologically meaningful channels.

Figure 15 Global variograms of the radar line data and modelled topography maps.

Figure 16 Subareas with corresponding data density, overlapped by the radar line data.

We further compare our approach with the two-point geostatistical modelling methods, kriging, and SGSIM (Fig. 12c and d). We observe that kriging produces the smoothest topographical model. The oversmoothed features are clear from the detailed cross section in Fig. 13. Besides, kriging is a deterministic modelling approach. It cannot generate multiple topographical models to quantify the spatial uncertainty. The SGSIM approach here uses local ordinary kriging; this way, non-stationarity is addressed by limiting the neighbourhood of spatial inference. As a spatial covariance-based approach, SGSIM is limited in its ability to capture complex morphological features, especially when the radar line data are very sparse. This limitation can be clearly observed in Fig. 14, where SGSIM does not produce morphologically meaningful channels in the area with sparse lines when comparing to our proposed approach. In Fig. 15, we further compare the empirical variograms from the modelled topographical maps using the three different approaches. It shows the DS using sampled TIs has reproduced the observed radar data variogram. The SGSIM maps also reproduce the variograms from the observed radar data, because they directly use the radar data variograms for modelling. However, the DS without TI sampling has a nugget (noise) effect. Overall, it shows the TI sampling approach performs the best in terms of improving the modelling speed, simulation quality, and capturing the spatial uncertainty. It is also important to note that there are other TI selection approaches. For example, Pérez et al. (2014) ranked TIs using high-order consistency with conditioning data, while Abdollahifard et al. (2019) used image contours to select compatible TIs. The unique contribution of our approach is that it quantifies posterior probabilities and uses a sampling method. In this sampling method, we may use different TIs for different realizations generated, as shown later in Fig. 17.

Figure 17 Two realizations of TIs assigned to the entire ASE area.

Figure 18 Realizations of the ASE topography after filling the radar line gaps. The ovals highlight the areas with high spatial variations across the realizations. The black area is the non-study area with no radar data.

4 Application to the entire ASE

4.1 Training image sampling and DS simulation

We apply the proposed methodology to fill the radar line gaps of the entire ASE area to generate high-quality topography maps at a resolution of 500 m. To address the spatial non-stationarity and sample TIs, we first divide the whole ASE area into local subareas. We use the following recursive steps to divide the entire ASE into L subareas based on the line data density:

• Step 1.

  Equally divide the ASE area into four quadrants in the shape of a square.

• Step 2.

  For each subarea, if it has more than N radar data points, continue to divide it into four equally sized areas.

• Step 3.

  Repeat step 2 until the number of data measurements is below the threshold number N. In this case, we set N=10 000 and divide ASE into a total of 56 subareas.

Figure 16 shows the final ASE subareas with the corresponding radar data density. Next, we calculate posterior TI probability for each subarea $P\left(\mathrm{TI}\left(A_i\right)\mathrm{|}{d}_{A_{\mathrm{1}}},\mathrm{\dots },d_{A_L}\right)$. We can then use the posterior probability to sample one single TI for each subarea. Figure 17 plots two realizations of the sampled TIs in the entire ASE space. For each TI realization, we run DS to fill the radar line gaps to generate high-resolution topography maps. To reflect the spatial uncertainty, 20 topography map realizations are simulated using 20 realizations of TI sets. The generated high-resolution bed topography is shown in Fig. 18.

4.2 Uncertainty in subglacial hydrological flow

We use the topographic realizations to investigate the sensitivity of subglacial water routing to topographic uncertainty. A water-routing model was applied to the 20 realizations generated from Sect. 4.1 to model the flow of water at the ice–bed interface. The direction of water flow can be determined by calculating hydraulic potential, ϕ, using the Shreve equation (Shreve, 1972):

$$\begin{array}{}\text{(14)}& \mathit{\varphi }={\mathit{\rho }}_{\mathrm{w}}gh+{\mathit{\rho }}_{\mathrm{i}}gH,\end{array}$$

where ρ_w is the density of water (1000 kg m⁻³), ρ_i is the density of ice (917 kg m⁻³), g is gravitational acceleration 9.8 m s⁻²), h is bed elevation, and H is ice thickness. The hydrological model was implemented using the Antarctic Mapping Tools (Greene et al., 2017) and the FLOWobj function and multiple flow direction (MFD) algorithm from the TopoToolbox package (Schwanghart and Scherler, 2014). These functions use the hydraulic potential gradient to compute flow accumulation, or the number of pixels that flow into another pixel. We use the “multi” setting, where downward flow in different directions is distributed based on a hydraulic gradient. We assume spatially uniform basal melt rates and that the subglacial pressure is equal to the ice overburden pressure. The results are compared to a hydrological model made using BedMachine topography (Morlighem et al., 2020), which was derived using the mass conservation inversion method.

The water-routing models vary across each realization (Fig. 19). In particular, the Thwaites Glacier tributaries flowing towards the grounding zone (the area where the ice meets the ocean and decouples from the bed) show significant differences across each realization. The average of the hydrological models across different realizations is different from the hydrological model made using the mass conservation topography, particularly in the main trunk of Thwaites Glacier (Fig. 19d and e). This demonstrates that deterministic DEMs cannot be used to sample the range of possible flow path locations, which could lead to the misinterpretation of hydrological conditions. In contrast, geostatistical simulation provides a framework for quantifying hydrological uncertainty with respect to topographic uncertainty.

Figure 19 (a) Subglacial water routing from mass conservation; (b, c) water routing using two topographic realizations from our DS simulation with TI sampling; (d) flow accumulation for mass conservation; (e) mean flow accumulation for 20 topographic realizations from our DS simulation.

Some of the modelled tributaries are located over a system of active subglacial lakes – lakes at the ice–bed interface that periodically drain and refill (Hoffman et al., 2020; Smith et al., 2017). These lakes are hypothesized to be hydrologically connected, with a drain and refill cycle that depends on the level of connectivity (Malczyk et al., 2020; Smith et al., 2017). Lake drainage events are sometimes associated with increases in ice velocity (Fricker et al., 2016; Stearns et al., 2008), making it important to characterize the connectivity of active lake systems. The topographic models created by the proposed method can be used to investigate the nature of hydrological drainage at Thwaites and highlight areas that require additional observational constraints.

5 Conclusions

We developed a non-stationary multiple-point geostatistical approach to fill large-scale geophysical data gaps and applied it to map high-resolution (500 m) subglacial topography of the Amundsen Sea Embayment in West Antarctica. The radar data gaps were filled using morphological features learned from high-resolution topographic training images. To reflect the geospatial uncertainty, we modelled multiple realizations of topography maps using 166 high-resolution training images from the Arctic and Antarctica. These training images represent the diversity of subglacial geologic settings. We have placed them in a publicly accessible repository for training subglacial topography models (see Data and code availability section). The TI repository can be further expanded in the future upon the acquisition of additional swath bathymetry and swath radar measurements.

Our major contribution was to show a probabilistic method to model posterior TI probabilities, then sample TIs to model the global non-stationarity in subglacial topography. This was achieved by probabilistically assigning non-stationary TIs from the provided repository to the local radar data. We used the collected 166 topographic training images as prior. The posterior distribution is calculated based on the modified Hausdorff distance between each TI and local radar data. To address the spatial correlation across the global area, we aggregate the TI probability between the local areas based on their spatial correlation. The aggregated posterior TI distribution then enabled us to sample training images. Finally, we ran DS to fill the radar line gaps. Multiple realizations of high-resolution topography maps were generated using multiple realizations of sampled training images. This non-stationary TI sampling framework avoids the use of auxiliary variables and arbitrary ad-hoc weighting. It has significantly improved the topography modelling quality from DS. It also dramatically reduced the DS running time from 21 to 1 h when given a large number of training images. Compared to the traditional deterministic interpolation (kriging) and two-point geostatistical simulation (SGSIM) approaches, our approach was shown to provide more realistic topographic maps for spatial uncertainty quantification, whilst retaining the spatial correlation measured by radar data.

We applied our proposed approach to fill the radar line gaps for the entire Amundsen Sea Embayment in West Antarctica. The improved modelling efficiency enabled us to simulate 20 realistic high-resolution topographic maps on a local PC. We then used the 20 topographic realizations to investigate the sensitivity of subglacial water routing to topographic uncertainty. The results reveal significant variabilities in the Thwaites Glacier tributaries across different realizations. These tributaries intersect a system of active subglacial lakes, which are hypothesized to be hydrologically connected and could have the potential to influence ice sheet velocity. The high hydrological uncertainty in this area highlights the need for additional measurement constraints. These findings demonstrate the utility of geostatistically simulating subglacial topography rather than performing deterministic interpolations. Our non-stationary MPS framework provides a path forward for implementing geostatistical simulations at continental scales.

Appendix A: PSO and optimal TI numbers

Figure A1 (a) PSO minimum distance vs. the swarm size n. (b) Profile log likelihood of the curve in panel (a), suggesting the optimal number of training images is three.

We perform particle swarm optimization to minimize the distance function $\mathrm{dis}\left({\mathbb{I}}_{\mathrm{TI}}\left(\widehat{\mathrm{TI}}\right),d_{A_i}\right)$ in Eq. (8). Following the PSO algorithm (Rezaee Jordehi and Jasni, 2013), we start with a random initialization of m selected TIs). Each individual TI selection is regarded as an individual particle (P_i). To find TIs that minimize the distance function, each particle will explore the whole TI space iteratively with a velocity V_i. The position (TI index) of P_i at time step t+1 is determined by its previous position P_i(t) and searching “velocity” V_i(t+1).

$$\begin{array}{}\text{(A1)}& P_{\mathrm{i}}\left(t+\mathrm{1}\right)=P_{\mathrm{i}}\left(t\right)+V_{\mathrm{i}}\left(t+\mathrm{1}\right)\end{array}$$

The velocity V_i(t+1) is determined by the particle's current logged best TI index $P_{\mathrm{i}}^{\mathrm{best}}$ and the best TI index $P_{\mathrm{g}}^{\mathrm{best}}$ for the whole swarm, as

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{V}_{\mathrm{i}}\left(t+\mathrm{1}\right)& =w\times V_{\mathrm{i}}\left(t\right)+c_{\mathrm{1}}\times r_{\mathrm{1}}\times \left(P_{\mathrm{i}}-P_{\mathrm{i}}\left(t\right)\right)\\ & +c_{\mathrm{2}}\times r_{\mathrm{2}}\times \left(P_{\mathrm{g}}-P_{\mathrm{i}}\left(t\right)\right),\end{array}\end{array}$$

where V_i(t) is the velocity from the previous time step. w is the “inertia weight” that controls the contribution of V_i(t) to V_i(t+1). A smaller w means less influence from the previous velocity and thus higher PSO exploration capability. Here, we set w as 0.8 according to the study by Han et al. (2010). r₁ and r₂ are two random numbers for a stochastic update of the velocity. They have a uniform distribution with the interval of [0,1]. c₁ and c₂ are the acceleration parameters that pull the particles towards $P_{\mathrm{i}}^{\mathrm{best}}$ and $P_{\mathrm{g}}^{\mathrm{best}}$. $c_{\mathrm{1}}=c_{\mathrm{2}}=\mathrm{2}$ are recommended for most optimization problems according to Ozcan and Mohan (1999). We adopt the recommended settings. The swarm size m also affects the PSO performance. So far, there are no exact rules for the selection of swarm size (Rezaee Jordehi and Jasni, 2013). Here, we use the size n of $\widehat{\mathrm{TI}}$ to determine the swarm size m. We create $m=\mathrm{10}\times n$ particles in the PSO population to enhance searching ability and running time.

Another important question is how to determine the optimal number n of $\widehat{\mathrm{TI}}$. To specify n, we use a profile log-likelihood approach from Zhu and Ghodsi (2006) and Honarkhah and Caers (2010). Specifically, we expect that the distance between training images and radar line data will decrease as we visit more training images. The PSO optimized distance should decrease dramatically when the optimal n TIs is visited and then start flattening out. Hence, there will be an elbow point corresponding to the optimal number of TIs. Based on the study in Honarkhah and Caers (2010), the elbow is found by maximizing profile log likelihood. We find the optimal number of TIs using the following steps.

1. Run PSO to obtain the minimized distances dis_n with different $\widehat{\mathrm{TI}}$ size n, where $n=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots },N-\mathrm{1}$.

2. For every n, we define two samples of ${\mathit{\phi }}_{\mathrm{1}}=\left\{{\mathrm{dis}}_{\mathrm{1}},{\mathrm{dis}}_{\mathrm{2}},\mathrm{\dots },{\mathrm{dis}}_n\right\}$ and ${\mathit{\phi }}_{\mathrm{2}}=\left\{{\mathrm{dis}}_{n+\mathrm{1}},{\mathrm{dis}}_{n+\mathrm{2}},\mathrm{\dots },{\mathrm{dis}}_N\right\}$. Calculate the log likelihood l_n(n) as

   $$\begin{array}{}\text{(A3)}& {\displaystyle }\begin{array}{rl}{l}_n\left(n\right)& =-n\log \left({\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}}}}\right)\sum _{i=\mathrm{1}}^n{\displaystyle \frac{{\left({\mathrm{dis}}_i-{\mathit{\mu }}_{\mathrm{1}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\\ & +(n-N)\log \left({\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}}}}\right)\sum _{i=n+\mathrm{1}}^N{\displaystyle \frac{{\left({\mathrm{dis}}_i-{\mathit{\mu }}_{\mathrm{2}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\end{array}\text{(A4)}& {\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}={\displaystyle \frac{\left(n-\mathrm{1}\right){\mathit{\sigma }}_{\mathrm{1}}^{\mathrm{2}}+(N-n-\mathrm{1}){\mathit{\sigma }}_{\mathrm{2}}^{\mathrm{2}}}{N-\mathrm{2}}},\end{array}$$

   where μ₁, μ₂ are the means of φ₁ and φ₂, By contrast, σ² is the common scale variance. σ₁ and σ₂ are the sample variances of φ₁ and φ₂.

3. Obtain the optimal (elbow) size $\widehat{n}$ based on the empirical maximum value of l_n(n).

We use the illustration case area A1 as an example to show how to select $\widehat{n}$. Figure A1a plots the PSO distance function between A1 radar lines and TIs with varying TI numbers. We can observe a fast drop of the distance at the beginning, and the distance then drops slowly after n=3. To find out the exact elbow, we calculate the log-likelihood values and plot them in Fig. A1b. Figure A1b clearly indicates that the optimal number of TIs is three.