Antarctic Surface Ice Velocity Retrieval from MODIS-Based Mosaic of Antarctica ( MOA )

The velocity of ice flow in the Antarctic is a crucial factor to determine ice discharge and thus future sea level rise. Feature tracking has been widely used in optical and radar imagery with fine resolution to retrieve flow parameters, although the primitive result may be contaminated by noise. In this paper, we present a series of modified post-processing steps, such as SNR thresholding by residual, complex Butterworth filters, and triple standard deviation truncation, to improve the performance of primitive results, and apply it to MODIS-based Mosaic of Antarctica (MOA) datasets. The final velocity field result displays the general flow pattern of the peripheral Antarctic. Seventy-eight out of 97 streamlines starting from seed points are smooth and continuous. The RMSE with 178 manually selected tie points is within 60 m·a−1. The systematic comparison with Making Earth System Data Records for Use in Research Environments (MEaSUREs) datasets in seven drainages shows that the results regarding high magnitude and large-scale ice shelf are highly reliable; absolute mean and median difference are less than 18 m·a−1, while the result of localized drainage suffered from too much tracking error. The relative differences from manually selected and random points are controlled within 8% when speed is beyond 500 m·a−1, but bias and uncertainty are pronounced when speed is lower than that. The result through our accuracy control strategy highlights that coarse remote-sensed images such as Moderate Resolution Imaging Spectrophotometer (MODIS) can still offer the capability for comprehensive and long-term continental ice sheet surface velocity mapping.


Introduction
Accurate surface flow velocity measurements are fundamentally important for glaciologists since they play a vital role in the regulation of the ice sheet dynamic system.Based on the assumption of steady state, the ice flux can be integrated by ice thickness and balance velocity.The balance velocity is depth-averaged vertically but commonly takes a ratio with observed surface velocity for simplicity.On the other hand, velocity can be used as initial input or boundary condition in numerical flowing models to invert basal and thickness parameters and predict the mass balance state of ice sheet, and thus the contribution to future sea-level rise [1,2].Considering a large population along the coastal area can be influenced, especially in the background of climate change, investigating how ice sheets and ice shelves evolve in response to increasing greenhouse gases(GHGs) in atmosphere is necessary [3].
According to the Intergovernmental Panel on Climate Change's Fifth Assessment Report (IPCC-AR5) [4], both the mass of input to and loss from Antarctic ice sheets still suffer from large uncertainty [5].On the one hand, in terms of "feeders," cyclone systems over the Southern Ocean have a direct influence on high latitude precipitation [6].On the other hand, however, there is high confidence that the largest ice losses occurred along the northern tip of the Antarctic Peninsula, where the collapse of several ice shelves in the last two decades triggered the acceleration of outlet glaciers [7].In the West Antarctic, a recent study revealed that flow accelerations across the grounding lines of the Amundsen Sea Embayment, the Getz Ice Shelf, and the Marguerite Bay on the western Antarctic Peninsula account for 88% of discharge increase in the last 7 years [8].The Pine Island Glacier [9] together with the Thwaites Glacier [10] both attract scientists' attention for their retreating and calving behaviors.Those two glaciers were found to speed up continually, although the acceleration rates have been decreasing to 1% and 8% [11].Consequently, the accelerated ice mass loss in the coastal regions of West Antarctica continues to outpace the gains made in the East Antarctic regions.Even in the East Antarctic, some glaciers will endure instability, such as the Byrd Glacier [12] and the Totten Glacier [13].Previously studies confirm a couple of factors contributing to this status, such as marine ice-sheet instability [14], ocean thermal forcing [15], and surface melting [16].
Satellite measurements of Antarctic glacier velocities have been carried out since the 1980s [17].Since the NASA Making Earth System Data Records for Use in Research Environments (MEaSUREs) InSAR-Based Antarctica Ice Velocity Map revealed a complex flow pattern from the coastal region to the deep interior in the Antarctic for the first time [18,19], ice surface velocity mapping at the regional to continental scale has progressed rapidly [20], not only in the Antarctic updated yearly [21], but also in Greenland Ice Sheet [22,23], Central Asia [24], Cirucum-Arctic [25], Alaska [26], Canada [27], and Patagonia [28].The improvements in contemporary remote sensing technology, either in optical or synthetic aperture radar (SAR), provided by new-launched earth observation satellites or constellations, such as Radarsat-2 [29], ALOS-PALSAR [27], TerraSAR-X and TanDEM-X [30], Landsat-8 [31,32], Sentinel-1 [22], and Sentinel-2 [33], make these projects possible.Based on a quick re-visit of, and the high-quality imagery from, these integrated missions, some of them have already been systematized into near real-time operational retrieval and distribution portals to provide insights to both cryospheric research community and the public [34], taking Global Land Ice Velocity Extraction from Landsat-8 (GoLIVE-V1), for example.
As an excellent data source commonly used in various cryospheric research fields, the Moderate Resolution Imaging Spectrophotometer (MODIS) demonstrates its capabilities in surface melting detection [35], water vapor estimation [36], digital elevation model calibration [37], etc.Perhaps as the result of the coarse spatial resolution, however, MODIS's potentiality to track ice sheet surface motion is seldom implemented.Haug et al. investigated Larsen-C ice shelf's velocity from 2002 to 2009 via MODIS, and the difference from Landsat's retrieval was within 30 m•a −1 [38].Coincidentally, the Larsen-C ice shelf was chosen to be the study area again by Chen et al. [39], and it was found to have accelerated continually since 2000.In fact, Wang et al. [40] and Li et al. [41] have exploited Declassified Intelligence Satellite Photography (DISP) images acquired during the Cold War period to trace back to earlier dynamic hints of the Larsen-C ice shelf.Most images that have been used to derive velocity in glaciology have been high or medium spatial resolution images, such as orthorectified airborne photos [42], SPOT [43], Landsat [44], ASTER [45], or Sentinel-2 [33].In addition, the technology innovation makes commercial cubesat constellation [46] and unmanned aerial vehicles (UAVs) [47] as the efficient sensors to monitor glacier behavior.Low-resolution images cover large section in one scene, making them suitable for investigating the velocity of large areas such as Antarctic ice shelves [48].Therefore, we speculate that MODIS might underperform in the glacier flow's velocity retrieval, since in most cases researchers have better choices.
While MODIS enjoys a shorter re-visit time and a wider swath in the high latitude region, the probability of cloud-free scenes increases tremendously, which means a continental mosaic can be compiled in a relatively shorter time slot.Hulbe et al. [49] has already tested MODIS-based Mosaic of Antarctic (MOA)'s performance in retrieving surface velocity of Ross Ice Shelves, the result of which provides external forcing for the long-term dynamic diffusion ice sheet model.However, regarding the basic tracking techniques, it inevitably yields gross errors in the primitive displacement field [50], especially where high variation dominates or surface melting and snow events transform former features.To enhance the primitive output but avoid laborious human sorting out, additional constraints are needed, such as a crude local statistics filter [51] or a physical assimilation mechanism [52].By expanding Hulbe C's operation to the whole continent, the aim of the present study is to explore the possibility of coarse resolution images over a long time baseline in accurately extracting a macro-scale velocity field from MOA datasets without much human interaction.

MOA
MOA datasets, together with a grain size map, are compiled by the National Snow and Ice Data Center (NSIDC) and University of New Hampshire [53].Their source images were collected from December of the prior year to February of the next year in the periods of 2003-2004 and 2008-2009, respectively, from the MODIS instruments on-board the NASA EOS Aqua and Terra satellites.The two digital maps achieve a better spatial resolution of about 150 m by leveraging pixel cumulating schema.Furthermore, the 125 m grid images have a spatial coverage identical to the Radarsat Antarctic Mapping Project Antarctic Mapping Mission 1 (RAMP AMM-1) 125 m mosaic.They provide a nearly-perfect cloud-free view of the ice sheet, ice shelves, and outcrop surfaces for land areas north of 60 • S, allowing us to investigate its potentiality in determining the displacement between two acquisition times.
This and the following datasets are all released and mapped as Antarctic 71 • S Polar Stereographic Projection (EPSG:3031).

Landsat-8 Rock Outcrop
An updated Landsat-8-based rock outcrop dataset of Antarctic was released in 2016 [54], which had become part of the Antarctic Digital Database (ADD).It used Landsat-8 multispectral data, including thermal bands, to generate a new automated methodology for delineation of rock outcrop in Antarctica.Overcoming the limitations of previous techniques, it enjoys several improvements, such as calibration overestimate areas and more consistent georeferencing.In areas south of 82 • 40 S, where Landsat-8 does not cover, or islands lacking suitable cloud-free images, the outcrops are derived from the existing ADD rock outcrop dataset to supplement.Thanks to the higher accuracies than previously achieved (proven as high as 85 ± 8%), here we utilize it as a stationary area mask to evaluate our velocity result's residual.
It is worth noting that another Antarctic rock outcrop dataset was published recently; the authors improved the detection techniques in the shadow area and achieved greater accuracy [55].We compared the two datasets in the experimental phase and found that Kang's result is featured with finer detail.The rock polygons smaller than the final velocity pixel with postings of 750 m are commonly found.Meanwhile, the spatial scale of the ADD dataset is suitable for following medium resolution analysis.

MEaSUREs Ice Velocity
MEaSUREs InSAR-Based Antarctica Ice Velocity Map, assembling 1384 tracks optical and SAR images from multiple sensors, takes full advantage of offset tracking and interferometric techniques (InSAR) [18].The well-corrected error of 17 m•a −1 , together with high spatial resolution and accurate geolocation, wholly make it a nearly perfect benchmark for future investigation.Source data are largely acquired during the International Polar Year (IPY) 2007 to 2009, as well as between 2013 and 2016.Additional data acquired between 1996 and 2016 are used as needed to maximize spatial coverage.It is set as the validation reference in this study, although it is worth pointing out that the data collection time period spans 20 years from 1996 for ERS-1 to 2016 for Radarsat-2, which is quite different from the acquisition time of MOA.

Feature Tracking
Coregistration accuracy of two MOA products was checked by selecting well-dispersed nunataks or outcrops from Landsat-8 Rock Outcrop.Considering the offsets of all of them are within half of a pixel (∼20 m, see Section 4.3 for detail), it can be safely accepted that the whole image's mis-registration is within 1 pixel as well.Feature tracking based on a normalized cross-correlation (NCC) algorithm has already been implemented with dozens of software applications, in which images from two different times are compared using correlation techniques to derive glacier displacement over the time period [56].Here we employ COSI-Corr [57] to acquire the raw displacement field because it outperforms other similar operators in most cases [58].This tool is originally written in IDL as an extension of ENVI software, and it embeds modified phase correlation, Frobenius norm, and singular value decomposition techniques [59] to achieve a more robust field.Based on the ICESat-1 Antarctic drainage system published by the Goddard Space Flight Center [60] and their respective flow characteristics, the whole continent was divided into seven drainages accordingly as shown in Figure 1, namely Amundsen, Dronning Maud, Filchner-Ronne, Lambert, Peninsula, Ross and Wilkes Drainage.Each one of seven drainages was tested iteratively to obtain optimal but slightly different tracking parameters.By adjusting reference and search window size, step, iteration times, and mask threshold, the output binary file comprises the Easting displacement component in the X dimension, the Northing displacement component in the Y dimension, and the signal-to-noise ratio (SNR) as the primitive displacement field.

Post-Processing
The primitive displacement fields cannot be analyzed directly since they are often contaminated by severe noise from spurious matching.Hence it is indispensable to correct the initial output before inferring valuable knowledge.In the following part, an accuracy improvement strategy consisting of three steps is developed.
Usually, the data below the certain SNR threshold, which denotes poor tracking quality, are sweeped out.Considering the spatial heterogeneity in an immense ice sheet, it is obviously unsuitable to use single criteria.To make the threshold adaptive, each drainage has its own threshold based on the residual velocity's decreasing pattern within an outcrop mask (see Figure 5).Here in this first post-processing step, the SNR threshold, which approaches the asymptote, was tuned and determined manually by visual inspection.The automated method could be developed in the future, allowing for efficiently filtering as many outliers as possible while retaining a majority of valid measurements.However, we found that the first procedure is still not always helpful in getting rid of all mis-correlation.
In order to improve the initial output further, the performances of various classical filter configurations were tested before the low-pass Butterworth filter was chosen.Assuming that the shear strain of ice is relatively low and that there is little velocity variation in the adjacent pixels, a low-pass filter acts as a second step of post-processing, like a smoother.To make full use of computational resources and speed up the numerical fast Fourier transform (FFT), the real and imaginary parts were assigned with the X and Y dimensions in the flow field, respectively.Afterwards, the filtered complex array was separated into real and imaginary parts again, as the smoothed X and Y dimensions in the flow field in subsequence of the inverse FFT transformation.
Considering the coarse resolution and long time span of the MOA dataset, a mask was applied in which both the X and Y dimension velocity magnitudes are less than 30 m•a −1 (about 1/5 pixel of MOA image) or larger than 1500 m•a −1 in the MEaSUREs dataset.On top of that, the statistics of the difference between MOA-based velocity and MEaSUREs velocity in each drainage are shown in Figure 6.To obtain a reliable final result, they were used to exclude any measurements whose difference with respect to the MEaSUREs dataset is beyond triple standard deviation recursively, as the third post-processing step.
The above three post-processing steps, including SNR thresholding, the Butterworth filter, and the triple standard deviation truncation, excluded most outliers in the raw tracking field.In the end, the results in seven drainages were resampled to 750 m and mosaicked into a complete Antarctic velocity dataset.In some cases, spatial interpolation is the optional choice to fill the data void when necessary (see Figure 8).All these processing and the following analysis steps are summarized as a flowchart in Figure 1.

Result and Analysis
It is shown in Figure 2 that the flowing pattern extending into the interior of the ice sheet through ice stream channels.By taking processing steps mentioned above in Section 3, we found most valid data located in peripheral ice shelves, which means that the dominant inland ice sheet's velocity is well below the 30 m•a −1 in both the X and Y dimensions.Although the ice stream may firstly manifest itself as a roughly 15-km-wide area of enhanced flow starting within 100 km of the ice divide [61], this characteristic is indistinct because its motion signal upstream is too weak to detect in such scale.Only in large ice shelves such as the Ross or Filchner-Ronne does the dendritic system develop further upstream through tributaries.

Streamline
Being the particle's displacement path in the case of steady flow, streamline (also referred to as flowline for the steady-state velocity field in some literature) is constructed to verify the local orientation of the velocity field [62].Usually the hypothetical particle tracks are constructed as the streamline by interpolating from the current velocity field.Setting the seed point's time marker as t(0), the location of t(i + 1) can be computed iteratively [38].If the spatial variation of ice flows is abnormal, the streamline can terminate abruptly rather than ending up in the ocean.Usually, the velocity magnitude's profiles along the streamlines show a monotonous increasing trend as the elevation decreases.Dozens of start points are settled in Larsen-C, Filchner-Ronne, Ross, and other smaller ice shelves.The streamlines align quite well with the flow features on the ice surface, implying that there has been no or little directional change in ice flow during the last several decades.Most of them can move along a local orientation, forming smooth curves, as shown in Figure 3, although a few are trapped in local extrema, a data void, or a drastic shear zone and thus cannot wind their way to the ocean or converge to the mainstream.Specifically, among 97 start points, 19 streamlines break off halfway, and most of them, such as the Lambert Glacier in Figure 3b, are trapped in upstream narrow bottlenecks through visual inspection.Even so, it is difficult to detect the onset of tributary flow from this result since successful tracking depends both on the absolute velocity magnitude and on topography and morphology.

Tie Points
Since the in-situ measurements are too scarce, ground truth data can hardly be found for appropriate places and times.To validate the final result, we selected 178 pairs of tie-points manually, which are evenly distributed in the final velocity field, both spatially and numerically.Almost all of these points are distinct features in the ice surface, including crevasses, moraine, deformation, and ice-shelf front.Velocities from tie points are compared with those extracted from the MOA by automated tracking in Figure 4. Showing a strictly linear relationship between them, the velocity result from tracking MOA agrees quite well with the tie point velocity, not only in the X and Y dimensions but also in the magnitude and direction domains.Statistical comparison shows that determination coefficients are all close to 1, and RMSEs are roughly 50 m•a −1 , the same order as the original spatial resolution of the MODIS image taking five years time span, although the largest outlier goes to the extreme reaching 200 m•a −1 .The mean direction difference is less than 1 degree with an RMSE of 8.4 degree.A variety of dynamic parameter can be derived based on the velocity field, among which the strain rate is of particular significance [63].The above performance suggests that such parameters in the ice shelf area can be roughly estimated from our velocity result since the spatial distribution of errors is uniform and unbiased.However, noting that the tie-point pairs mentioned above are mostly derived from distinct features, those RMSEs are presumed to be underestimated since the majority of tracking pixels from indistinct features are not included in this comparison.Hence, the statistics can only represent areas with highly recognizable surface features.The different performance between tie-points and randomly sampled points in comparison with MEaSUREs is followed in Section 5 later.

Residual in Rock Outcrop
Mismatches or outliers are identified and removed using a threshold value of SNR, which is calculated by the correlation peak and background in every tracking's computation as the matching confidence indicator.The choice of the threshold is a compromise between removing most of the mismatches and retaining the interesting information at the same time [24].Assuming that the coregistration of MOA 2004 and MOA 2009 is accurate enough and outcrops are presumed to be stationary in this scale during the observational period, the residual velocity field in the stable rock area can be utilized to determine the appropriate SNR threshold.Figure 5 shows the distributions of velocity measurement's SNR in each drainage and their mean residual speed in rock mask with respect to SNR.The actually thresholds we used in the processing step are shown as vertical dashlines.The corresponding thresholds of the SNR filter were detected as the tipping point where residual is going to slash abruptly.Due to various flow characteristics, the climate regime, and the observation conditions such as the illuminating condition and the imaging quality in seven drainages, there are slightly different SNR thresholds in different regions.It is quite clear that the residual speed decreases to zero as the SNR threshold approaches one, both in our result and MEaSUREs dataset.Furthermore, in most cases, the residual of MEaSUREs (blue line) exceeds that of our result (red line).Since the random error originating from the stable rock area of MOA-based velocity by feature tracking is in no way better than that from the MEaSUREs dataset by InSAR methods, this residual difference indicate that the mask area of the Landsat-8 Rock Outcrop over-extends to the non-stationary ice sheet.The only exception occurs in the Wilkes Land Drainage, where the Transantarctic Mountains occupy a large amount of area of Victoria Land, while only a few small-scale inactive ice areas exist.For the same reason, the largest residual also appears in the Wilkes Land drainage for our result.However, for the MEaSUREs dataset, the largest residual occurs in Amundsen, where that of MOA-based velocity seems to remain invariable.By visual inspection, both overspread rock mask, resulting from snow coverage and from the shadow effect [54] in the ice sheet, and too much data void for the MOA-based velocity in this drainage contribute to this phenomenon.Generally, the residual error in the rock mask is well within 15 m•a −1 when the SNR is set beyond 0.985.Thus the systemic errors can be removed if the residual was accessed properly in a small basin.However, this value cannot be assigned to a constant in a large division such as our drainage because it is spatially dependent as well.

Comparisons with MEaSUREs
Most MEaSUREs source data were collected between 2006 and 2009 [19], coincidentally falling within the time span of two MOA datasets, which are from the end of 2003 to the beginning of 2009.Assuming that velocity evolves gradually rather than drastically with respect to time (except for surge-type glaciers, which is not typical in Antarctic), a discrepancy resulting from the time difference should be limited except for a few high variable glaciers.To better evaluate the accuracy in the seven drainages, we analyzed the difference between MOA-based velocity and the MEaSUREs dataset, which is resampled to 750 m by the bi-linear interpolation method, pixel to pixel in the remaining valid area.
Figure 6 represents the histograms of the velocity difference in the X and Y components in seven drainages, respectively.Considering the granularity of source images imported into the tracking program, Figure 6c,e shows that the final velocity retrieval from MOA is excellent in Ross drainage, where the mean and median are less than 5 m•a −1 , and the Filchner-Ronne drainage, where the mean and median range from 4 to 13 m•a −1 .Variance is less than 30 m•a −1 in these two drainages as well, roughly half of the RMSE in comparison between the tie-point pairs and our result.This value of 30 m•a −1 is the same as the velocity retrieval difference between MODIS and Landsat, whose outstanding spectral and georeferencing characteristics have already been recognized [39].In drainages where moderate scale ice shelves dominate, such as Dronning Maud, Lambert, and Peninsula, the result is yet acceptable since the mean and median of the difference is still relatively low, and the variance of the difference is the same order of that RMSE.Noting that, along the coastline of Dronning Maud drainage, Fimbul and other small ice shelves flow primarily in the positive Y direction but the Brunt Ice Shelf flows toward the negative X direction, the histogram in the Y direction shows a bimodal pattern.This kind of behavior appears most prominent in the Wilkes Land Drainage since it stretches most widely in longitude.In the worst cases-the Amundsen and Wilkes Land drainages-the discrepancy is considerable, resulting from both surface decorrelation and different acquisition times.Featuring intensive acceleration and a retreating grounding line, the Pine Island Glacier [64] and its neighbor, the Thwaites Glacier [14] are both located in the Amundsen drainage.There are hardly any valid data available in the main trunk of these two dominant glaciers.The corresponding deformation can not only ruin the matching features but also affect the accuracy of nearby areas in filtering processing.Even worse, if taking the observation time of MEaSUREs into account, a great deal of source data from ERS-1 and ERS-2, recorded as early as 1996 in this region, which is more than one decade before MOA's compilation, is integrated.Therefore, we can infer that only a fraction of inconsistent data is caused by tracking outliers.In term of the Wilkes Land drainage, where the vast majority of glaciers span across only dozens of kilometers as a result of steep topography along the coastline, only a few ice tongue or ice shelf motion signals are represented correctly.Six of them are indicated by a digital number from one to six in Figure 2, among which the Mertz Ice Tongue was already broken into two pieces by a floating iceberg B-09B in 2010 [65].
The effectiveness of filter via shrinking within a 3σ range is obvious by comparing the distribution and statistics from each drainage in Figure 6.In both X and Y components, the difference's distribution is all symmetrically bell-shaped roughly, approximating one or multiple Gaussian distribution curves combined.Outliers that exaggerate variance are truncated off, variance decreasing to 1/2-1/8 of the primitive result.Although this method originates from the inherent statistical characteristic of data, the a priori reference dataset is needed to work out.In addition, the statistics of difference between these two velocity datasets implies that the mean and variance both reduce remarkably, but the median, which is not sensitive to extreme outliers, remains of the same order after filtering.Furthermore, the scale of complete and continuous tracking area influences the statistics, namely in the Ross and Filcher-Ronne drainages, and the tracking techniques perform better than that in the Amundsen or Wilkes Land drainage, where only localized ice stream or ice shelves are situated.
To investigate the relative difference from MOA-based velocity with respect to the MEaSUREs dataset, 3000 points were strewn randomly in a valid area.Figure 7 depicts these points' differences between MOA-based and MEaSUREs velocity in magnitude from both tie points and random points.As shown in Figure 7a, our results are systematically biased to some extent, since it is more prone to underestimate than MEaSUREs, even though the mean relative difference is constrained within 8% beyond 500 m•a −1 .Like streamline analysis in Section 4.1, the limiting velocity detection from MOA is about 400 m•a −1 or so, below which the streamline can hardly maintain itself, and the confidence for the raw tracking result is relatively low.For random points, as is evident from the left panel, the underestimate pattern is severe when velocity is less than 400 m•a −1 .Rather, it is not the case for manually selected tie points since they are evenly distributed above and below the zero reference line within a narrow range, regardless of the MEaSUREs velocity.Therefore, selecting tie points with the distinct feature is more or less a matter of "cherry-picking," which reminds us that the interpretation of RMSEs in Section 4.2 might be underestimated.Sometimes a cluster of tracking exceptions can happen when the surface feature changes significantly between 2004 and 2009.For example, in the middle Amery Ice Shelf, intensive supraglacial ponds scattered as dense dark patches occurred in 2004, indicating that the instantaneous environment is warm enough to melt the snow and ice, but they almost disappeared entirely in 2009.As illustrated in Figure 8, the velocity field before filtering repletes with outliers in the red rectangle, and they are restored into a more smooth field afterwards.The evident contrast between the blue and black arrows in the middle part of the Amery Ice Shelf demonstrates the effectiveness of our accuracy constraint strategy.However, some of the erroneous results can still survive even after our correction strategy, since profound image variation like this not only results in failure locally but also propagates incorrect results upstream and downstream.

Discussion
Retrieving the ice velocity via the tracking algorithm is usually conducted over the high-resolution remote sensing image.Since little research involves coarse imagery such as MODIS, we adopted several tailored procedures to achieve a better final result.In this section, the data quality, processing performance, and potential promotions are discussed.
It is worth pointing out a couple of co-registration bias is simply ignored in the processing steps and analysis section.This can exert a little bias to pixel alignment, such as ocean tidal motion can entail up to a few meters height anomaly to the floating ice shelf [66].However, it has to be taken into account in the interferometry method [67].Negligible effects like this also include atmospheric artifacts, errors in sensor orientation and the ortho-rectification model, and DEM discrepancy through time and image distortion from the altitude effect [68], especially in steep margins of ice sheets.Given that the magnitude of pixel offsets generated by the above sources are usually two orders less than the spatial resolution of original MODIS image, their effects are not taken into account either.
Nevertheless, image coregistration can never be absolutely accurate, and the displacements generated from the misalignments among co-registered images need to be removed from the total composite displacement field.A common strategy to limit such error is to test the displacement field where no movement should be observed, i.e., in rock outcrops.When motions are observed in these areas, the errors can be spatially uniform, or present heterogeneous variations.In the former case, the residual can be used for the removal of co-registration errors.However, for our case, this operation is infeasible since the residual is non-uniform within each drainage due to the large scale across thousands of kilometers, just like between various drainages.Dividing them into a smaller drainage system might help, but that would require either time-consuming work in parameter tuning or locally adaptive matching operators.At the same time, if doing so, the difference distribution in Figure 6 should reveal only one peak in the histogram in a homogeneous flowing regime with a single dominant direction.
In terms of different final performance in each drainage division, the detail varies considerably from one drainage system to another, even within the same basin.The worst case refers to the Amundsen drainage, where warm ocean water beneath the continent's floating margins is eating away at the ice attached to the seabed [69].For example, the Pine Island Glacier, the grounding line of it has been continuously retreating thanks to heavy forcing by the warm ocean from below.It experienced as large as ±150 m•a −1 variation in velocity at a different location along one specific profile from 2009 to 2015 [70].Such drastic variation would no doubt sharply alter its surface feature.As shown in Figure 8, surface melting can be blamed for another reason as well.At least for large ice shelves such as Filchner and Ross, the tracking velocity is believed to reveal the real motion, and the result is comparable to Hulbe's experiment in the Ross Ice Shelf [49].In the Larsen-C ice shelf, our result seems even to achieve a better spatial coverage than the original MODIS images [38].
The longer the time delay between images, the lower the error related to the measured displacement per year, but the larger the error in finding an exact match, thanks to the surface pattern changes through time.In this research, the uniformity of time span in a mosaic from area to area complicated this problem.We treated it with a constant of five years for every pixel between two MOA images without special adjustment during processing.However, for most extreme cases, such as a pixel in MOA2004 observed at the very end of February 2004, and a corresponding pixel in MOA2009 observed at the very beginning of December 2008, the time span is nearly three months less than exactly five years.A relative error of up to 5% could be attributed to the above time difference, giving rise to both directional odds in slow moving areas and a large magnitude disparity in fast moving areas.In addition, while seasonal changes in glacier velocity do occur across the Antarctic, their effect is canceled out during the inherent averaging process as well.According to the "On Thin Ice" project website, the MOA2014 is in preparation now.After its release, another velocity field tracking from MOA2009 and MOA2014 would make monitoring Antarctic's spatial and temporal possible.
To improve the performance of raw tracking displacement, some constraints based on either signal processing or physical mechanisms are indispensable.Filtering of mismatches from a single band image pair is still difficult and remains a challenging topic in post-processing after raw tracking.Usually, the SNR and correlation coefficient are used to detect spurious matches [71].However, as the SNR and correlation coefficient can vary depending on the visual contrast and feature distribution, the user-specified thresholds might be arbitrary and even remove correct ones [58].With respect to this deficiency, it seems that the SNR might not be a perfect indicator to label tracking quality on account of insensitivity.Several studies have designed complicated techniques, including manual extraction of glacial areas [72], prior knowledge of flow direction or maximal displacement [62], and statistics of multi-bands within its neighborhood [64], to erase mismatches and preserve as many good matches as possible, leaving thus only those who meet the criteria.Combining the interferometry and offset tracking, Joughin et al. achieved a seamless velocity mosaic of the Amery Ice Shelf [73], also offering valuable insight into the ice sheet study.Recently, some researchers have tried to merge displacement measurements from multiple sources by the fusion technique [74].
Here in our research, after excluding outliers by SNR thresholding, low-pass Butterworth filtering in the frequency domain was employed to smooth the remaining field, although physical assimilation based on Navier-Stokes equations is worth further investigation [52].To leverage the triple-sigma rule of thumb, our approach requires an a priori velocity field with the same spatial coverage to derive standard deviation statistics.Similarly, an algorithm consisting of three processing steps was developed as the general filter framework for the ice surface velocity field deriving from remote sensing methods [75].The MEaSUREs dataset acts as an excellent reference, although the a priori reference field might be absent in some cases.Apart from outliers, some valid measurements with high strain rates in the shear margin can also be removed inevitably during accuracy control strategy.Although we did not deal with direction directly, the X and Y components in the filtering process inherently include the direction constraint indirectly.In addition, it would be better if no significant rotation or torsion [42] existed in two images, since otherwise some other techniques, such as affine transform and the SIFT operator, might be introduced [76].Such deficiency, such as the aberrant phenomenon in Figure 8, cannot be tackled by any statistical or physical methods independently unless another data source is provided.In addition, featured surfaces with crevasses are usually only present due to the stream confluence to ice shelves, dropping out of the vast upstream inland region.All of the above show that our final result is full of data voids and that other datasets or spatial interpolation are required to complement or to provide more complete spatial coverage, respectively.

Conclusions
This research demonstrates that MOA datasets are capable of generating the flow pattern of Antarctic ice streams, especially in tremendous ice shelves whose flow speed is relatively high enough to be detected by the tracking algorithm while preserving its surface morphology.In general, the accuracy is acceptable with the RMSE less then 60 m•a −1 if appropriate parameters are provided.The magnitude of the X and Y components in whole drainages as well as the local orientation variation agree well both with the surface texture of ice sheets in the original images and with the a priori authoritative MEaSUREs velocity dataset.
However, for small-scale drainages with high spatial and temporal flow variability, such as that at the Pine Island Glacier or the Thwaites Glacier in the Amundsen sector, the signal can barely be interpreted correctly.In addition, the heterogeneity of the melting pattern in the summer from each other limits the tracking ability as well.Considering the inherent noise in the primitive output, the post-processing correction is indispensable for accuracy control.The strategy of the SNR thresholding, the Butterworth smoothing filter, and the iterative triple standard deviation truncation was proven effective and efficient in weeding out outliers, although more robust physical assimilation is worth further studying.
The final velocity result shows that in spite of the coarse spatial resolution compared with the Landsat series, the MOA dataset with rapid revisiting and time consistency shows good performance in terms of the retrieval of ice surface velocity of ice streams and ice shelves.

Figure 1 .
Figure 1.Flowchart of processing and analysis steps.

Figure 2 .
Figure 2. MOA-based Antarctic Surface Ice Velocity with stretched MOA 2009 as background, the seven divisions are combined from Zwally's 2012 Antarctic Drainage System [60] (their names are just for notation in latter analysis, rather than exact demarcation of geographical entities), Landsat-8 Rock Outcrop areas, and all tie-points distribute all over the continent.Brown enlarged cross mark exemplify two neighboring tie-points in surface crevasse of the Larsen-C shelf in the close-up above.The floating velocity field of Drygalski Ice Tongue is enlarged in the close-up below.Digital numbers denote glaciers and ice shelves discussed in Section 4: 1 the Drygalski Ice Tongue, 2 the Cook Ice Shelf, 3 the Mertz Ice Tongue, 4 the Totten Glacier, 5 the Shackleton Ice Shelf, 6 the West Ice Shelf, 7 the Fimbul Ice Shelf, 8 the Brunt Ice Shelf, 9 the Pine Island Glacier, and 10 the Thwaites Glacier.

Figure 3 .
Figure 3. Start points and corresponding streamlines in (a) the Larsen-C Ice Shelf, (b) the Filchner-Ronne Ice Shelf, (c) the Donning Maud Land, (d) the Amery Ice Shelf, and (e) the Ross Ice Shelf, with transparent MOA tacking velocity field overlapping on the enhanced MOA 2009 image as the background.

Figure 4 .
Figure 4. Velocity comparison between MOA-based velocity and manually selected tie points offsets in the (a) X dimension, (b) Y dimension and (c) magnitude, and (d) direction domains.

Figure 5 .
Figure 5.The SNR distribution histogram, the residual speed with the SNR in the rock outcrop, and the corresponding SNR threshold labeled as a vertical dash-line in seven drainages: (a) Amundsen, (b) Dronning Maud, (c) Filchner-Ronne, (d) Lambert, (e) Peninsula, (f) Ross, (g) Wilkes Land.Note that the scale of the left Y axis in each panel differs.

Figure 6 .
Figure 6.Distribution of difference between MOA-based velocity and MEaSUREs velocity, before (upper row) and after (lower row) the Butterworth filter and the 3σ truncation iteration, in seven drainages: (a) Amundsen, (b) Dronning Maud, (c) Filchner-Ronne, (d) Lambert, (e) Peninsula, (f) Ross, and (g) Wilkes Land, respectively.Note that the value of all three statistics, mean, median, and variance, drops after the operations.

Figure 7 .
Figure 7. Relative difference with respect to MEaSUREs in (a) 3000 randomly strewn points and (b) 178 manually selected tie points.The blue margins represent 95% confidence interval, and vertical dash lines represent 500 m•a −1 .

Figure 8 .
Figure 8.The velocity field in the middle Amery Ice Shelf, with the transparent speed of the MEaSUREs dataset overlapping the enhanced MOA2009 image as background.Black and blue arrows represent the primitive tracking result and the final refined result of our accuracy constraint strategy, respectively.Enhanced original 2004 and 2009 MOA images of the same region are shown on top of the figure.The location is shown in the red box of the inset.