# A pipeline for automated processing of Corona KH-4 (1962-1972) stereo
imagery

Sajid Ghuffar
School of Geography and Sustainable Development
University of St Andrews, St Andrews, UK
& Tobias Bolch
School of Geography and Sustainable Development
University of St Andrews, St Andrews, UK
Ewelina Rupnik
LASTIG, Univ. Gustave Eiffel
ENSG, IGN, F-94160 Saint-Mande, France
&Atanu Bhattacharya
Department of Remote Sensing and GIS
JIS University, Kolkata, India
Corresponding author, Email<EMAIL_ADDRESS>Also affiliated with
Department of Space Science, Institute of Space Technology, Islamabad,
Pakistan

###### Abstract

The Corona KH-4 reconnaissance satellite missions from 1962-1972 acquired
panoramic stereo imagery with high spatial resolution of 1.8-7.5 m. The
potential of 800,000+ declassified Corona images has not been leveraged due to
the complexities arising from handling of panoramic imaging geometry, film
distortions and limited availability of the metadata required for
georeferencing of the Corona imagery. This paper presents Corona Stereo
Pipeline (CoSP): A pipeline for processing of Corona KH-4 stereo panoramic
imagery. CoSP utlizes a deep learning based feature matcher SuperGlue to
automatically match features point between Corona KH-4 images and recent
satellite imagery using SuperGlue to generate Ground Control Points (GCPs). To
model the imaging geometry and the scanning motion of the panoramic KH-4
cameras, a rigorous camera model consisting of modified collinearity equations
with time dependent exterior orientation parameters is employed. The results
show that using the entire frame of the Corona image, bundle adjustment using
well-distributed GCPs results in an average standard deviation (SD) or
$\sigma_{0}$ of less than $2$ pixels. We evaluate fiducial marks on the Corona
films and show that pre-processing the Corona images to compensate for film
bending improves the accuracy. We further assess a polynomial epipolar
resampling method for rectification of Corona stereo images. The distortion
pattern of image residuals of GCPs and y-parallax in epipolar resampled images
suggest that film distortions due to long term storage as likely cause of
systematic deviations. Compared to the SRTM DEM, the Corona DEM computed using
CoSP achieved a Normalized Median Absolute Deviation (NMAD) of elevation
differences of $\approx$ $4\text{\,}\mathrm{m}$ over an area of approx.
$4000\text{\,}{\mathrm{km}}^{2}$. We show that the proposed pipeline can be
applied to sequence of complex scenes involving high relief and glacierized
terrain and that the resulting DEMs can be used to compute long term glacier
elevation changes over large areas.

_Keywords_ Declassified Imagery $\cdot$ Panoramic Cameras $\cdot$ SuperGlue
$\cdot$ DEM $\cdot$ Glacier Changes $\cdot$ Epipolar Resampling

##  1 Introduction

The US Corona reconnaissance program consisted of a series of low-Earth orbit
satellite missions. These missions were designed by the Central Intelligence
Agency (CIA) for acquiring photographic images of primarily the former Soviet
Union for strategic surveillance during the Cold War era [1]. Panoramic
imagery from Corona missions is the first Earth observation dataset, which
covers large part of the Earth’s landmass. The combination of high spatial
resolution of up to $1.8\text{\,}\mathrm{m}$ along with stereoscopic coverage
makes this data highly valuable in numerous applications [2, 3]. In spite of
that, the potential of this data has remained largely unexploited due to the
difficulties in modeling of the panoramic geometry combined with the presence
of film distortions and the limited availability of image metadata.

The Corona spy satellite program was approved in the year 1958 and there were
102 Corona missions between 1960 and 1972. The data from these missions
consisting of more than 800,000 images, were declassified in the year 1995
through a presidential order (Executive Order 12951). The code names for
Corona series are KH-1, KH-2, KH-3, KH-4 KH-4A, KH-4B, where KH means Key-
Hole, a designation used for spy satellites [4, 1]. In the KH-1 to KH-3
missions, the satellite imaging payload consisted of a single panoramic camera
with ground resolution between 7.5-12 m. The KH-4 series consisted of a dual
panoramic camera system with first mission in the year 1962. The dual camera
system consisted of fore (forward) and aft (backward) looking cameras with a
$15^{\circ}$ off nadir view angle resulting in a convergence of angle of 30∘
and baseline to height ratio of 0.54. These cameras consisted of a rotating
lens and scan arm with a narrow slit that sequentially exposes a stationary
film through a $70^{\circ}$ rotation of the scan arm (Figure 1). The lens
system of the fore and aft looking cameras rotated in opposite direction to
minimize perturbations on the spacecraft.

The Corona panoramic images were captured on photographic film that was
returned to the Earth using a re-entry capsule and recovered by the US Air
Force. The technical design details of the camera systems on board Corona
mission are given in the declassified documents available from National
Reconnaissance Office (NRO) [5]. The Corona missions have undergone continuous
design improvement since their inception and the documents and reports
available from NRO provides valuable insights in to the design and working of
these camera systems [6, 7].

The cameras in the Corona missions were designed by Itek Corporation and the
main difference in the KH-4, KH-4A and KH-B cameras was the utilization of a
rotating lens/scan arm in KH-4B camera system instead of an oscillating scan
arm in KH-4/4A camera systems [8]. This provided more stability to the
platform and allowed faster scanning rates, which enabled KH-4B satellites to
be operated at a lower altitude and achieve a ground resolution of up to
$1.8\text{\,}\mathrm{m}$. There were 17 KH-4B missions with the first launch
in 1967 and the final launch in 1972. Table 1 summarizes the main
specifications of the KH-4 series missions. There were further variations in
the image motion compensation system, camera filters, slit width and panoramic
geometry reference data exposed on the film in the KH-4 series missions.

Table 1: Corona KH-4 mission and camera specifications [6, 9, 10]. Due to a non circular orbit, the satellite altitude varies, which also changes the ground resolution and footprint. | KH-4 | KH-4A | KH-4B
---|---|---|---
$\\#$ of Missions | 26 | 52 | 17
Years | 1962-1963 | 1963-1969 | 1967-1972
Recovery Vehicles | 1 | 2 | 2
Orbit Perigee (km) | 200 | 180 | 150
Camera Name | Mural | J-1 | J-3
Lens/Scan Arm Motion | Reciprocating | Reciprocating | Rotating
Focal Length (mm) | 609.6 | 609.6 | 609.6
Stereo Angle (deg) | 30 | 30 | 30
Panoramic FOV (deg) | 70 | 70 | 70
Best Ground Res. (m) | 7.5 | 2.75 | 1.8
Ground Footprint (km) | 20$\times$280 | 17$\times$232 | 14$\times$ 188

In this work we present CoSP: A pipeline for automated processing of Corona
stereo imagery, which implements a rigorous Corona camera model and automates
the processing over the entire Corona image swath and in multi-image
configurations. In addition to that, this paper contains following
contributions:

  * •

We devise a scheme to automatically generate GCPs between Corona scenes and
modern satellite imagery of different spatial resolutions, such as the images
from Planet Labs nano-satellites and Landsat-7;

  * •

We experimentally show that the rigorous camera model with the time-dependent
collinearity equations can achieve a SD of better than 2 pixels across the
entire image footprint. We also show that estimated camera parameters are
consistent with the Corona camera characteristics and the orbital parameters
and different IMC mechanisms in KH-4A and KH-4B are also observed in the
estimated parameters.

  * •

We perform an investigation of the internal camera calibration through the
reference data exposed to the film (e.g., rail holes, PG stripes). We propose
a correction scheme and show that by removing the observed deviations, the
accuracy in 3D can be improved;

  * •

We further show that the entire Corona stereo pair can be rectified with a
residual y-parallax less one pixel but suspected film distortions lead to
larger systematic y-parallax of up to 5 pixels in certain parts of the image;

  * •

We show a large-scale application of this pipeline by applying it to 24 Corona
images, over a footprint of 160x200 km.

##  2 Related work

Due to the availability of images from the 1960s along with high spatial
resolution, Corona images have enormous value in detecting long-term changes
on the Earth’s surface, such as changes of coastlines [11], lakes [12] and
urban areas [13]. Most of the earlier work on Corona imagery, however, has
focused on applications in archaeology and cryosphere. Corona imagery has been
used for identification and interpretation of archaeological features such as
crop marks and built-up structures [2, 14, 15], computation of long term
glacier area changes [16, 17, 18, 19], estimation of glacier volume change and
mass balance [20, 21, 22], assessing area elevation and surface velocity
changes of rock glaciers [23, 24] and identification/mapping of glacial lakes
[25, 26]. Bhattacharya et al. [27] have used DEMs generated from Corona KH-4
imagery for long term glacier mass balance estimation for several areas in
High Mountain Asia. However, due to the manual work involved in their Corona
image processing workflow, their work focused on certain smaller areas within
the whole region. In this context, it should be mentioned that declassified
images from KH-9 Hexagon (1971-1984) mapping cameras have been used for
mapping of large areas and computation of glacier mass balance [28, 29, 30,
31, 32, 33]. In comparison to Corona, the processing of Hexagon KH-9 mapping
camera imagery is less complicated because of frame camera geometry and
availability of reseau grid for estimation of film distortions. The advantage
of Corona imagery is higher spatial resolution in comparison to Hexagon
(6-$8\text{\,}\mathrm{m}$) as well as image acquisitions earlier in time [31].
This can extent the glacier mass balance time series, as well map the dynamics
of slow flowing landforms such as rock glaciers at a high resolution over
longer time period [24].

The panoramic imaging geometry of the Corona camera systems differs from the
frame camera geometry. As a result, different camera models and
parameterization have been used in earlier studies to model the Corona
panoramic geometry. Sohn et al. [34] evaluated three different camera models
for KH-4B panoramic imagery. Their first model used modified collinearity
equations, which included panoramic geometry transformation along with scan
positional distortion and Image Motion Compensation (IMC) terms to model the
difference between the frame camera geometry and the panoramic geometry. The
scan positional distortion is caused by the movement of camera during the
scan, while IMC compensates for the resulting image motion to avoid blurring.
Their second model included the time dependent exterior orientation parameters
in the collinearity equations to model the movement of the satellite over time
as the rotating lens scans the scene. Their third camera model used the
Rational Polynomial Coefficients (RPC) [35, 36] to model the relationship
between the image and the object coordinates. Their results showed that the
first two models with modified collinearity equations resulted in a height
RMSE of around $4\text{\,}\mathrm{m}$ using GCPs covering an area of around
$560\text{\,}\mathrm{k}\mathrm{m}^{2}$ (approx. 20% area of entire stereo
overlap), which was better than the height accuracy achieved using the RPC
based model. Shin and Schenk [37] also used a camera model with time dependent
exterior orientation parameters but only included the motion along the flight
direction as an additional parameter. Their results showed a height RMSE of
$12.4\text{\,}\mathrm{m}$ on the check points spread over entire Corona stereo
pair. The NASA Ames Stereo Pipeline (ASP) [38] has implemented KH-4B panroamic
camera model based on the model of [37] and motion compensation from [34].
Lauer [39] has used a fisheye camera model for Corona imagery and obtained an
average height accuracy of $12.5\text{\,}\mathrm{m}$ for the control points.
Jacobsen [40] used the panoramic transformation terms in the perspective frame
camera model and observed an average standard deviation of
$11.4\text{\,}\mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l}\mathrm{s}$
using GCPs over multiple consecutive Corona stereo pairs. Multiple earlier
studies [20, 21, 27] have used RSG software (Remote Sensing Software Package
Graz) developed by Joanneum Research Graz, for processing of Corona imagery to
create DEMs for glacier volume change estimation. RSG Corona camera model is
also based on panoramic coordinates transformation in the collinearity
equations similar to the first model of Sohn et al. [34] but doesn’t include
IMC or scan positional distortion terms. These studies using RSG have reported
on average 2.5 - 3 pixel image residuals using manually extracted GCPs.

Figure 1: Geometry of the Corona panoramic stereo imaging system (Reproduced
and modified from [6])

Analyzing the earlier work on Corona image processing shows that further
investigations are still required to assess the accuracy of the camera model
over the entire panoramic image. Although Sohn et al. [34] presented a
rigorous approach for modeling of Corona imagery, their evaluation was only
performed over a relatively small area in comparison to the whole Corona image
footprint. Thus, additional evaluation is required to assess if similar
accuracy can be achieved over the entire Corona image. There are also
differences in the model parameterization in different studies. For example,
RSG Corona camera model as well as in the work of Jacobsen [40], the camera
model did not include scan motion distortion and the IMC terms. This is
partially grounded by the fact that the IMC is designed to counteract the
image motion. In a perfect scenario the effect of either is cancelled out. In
reality, the image motion is compensated only to a certain degree, and the
remaining residual causes systematic errors [41].

In addition to the camera model, the incorporation of panoramic geometry
reference data exposed on the film during the image acquisition has not been
thoroughly evaluated in previous works. [40] has shown that the bending of the
film has a significant magnitude and it should be corrected for achieving
higher accuracy.

##  3 Methodology

The fore and aft looking camera of Corona KH-4 missions performs a
$70^{\circ}$ scan in the across track direction with a rotating slit
sequentially exposing a static film. The entire Corona image is not captured
at one time instant as the platform is moving while the camera scans the scene
in the across track direction. We therefore decided to use the modified
collinearity equations with time dependent parameters to model the imaging
geometry of Corona panoramic cameras based on the work of Sohn et al. [34].
This model has the advantage of modelling the motion of the platform during
the scan with an explicit parameterization of the time dependent exterior
orientation parameters. Main components of the Corona stereo pipeline are
shown in the Figure 2. The first steps involve recovery of the film geometry
(including the individual scans’ stitching and the reconstruction of the
panoramic geometry) (Section 3.1), followed by GCPs and tie points extraction
(Section 3.2), the estimation of the camera model through bundle adjustment
(Section 3.3), epipolar resampling (Section 3.4), dense matching, 3D
intersection and fine coregistration (Section 3.5).

Figure 2: Workflow of the Corona Stereo Pipeline

###  3.1 Recovery of the film geometry

##### Film scans stitching

The physical length of the film (across track direction) containing one
panoramic scan is around $745\text{\,}\mathrm{mm}$, while the film width
(flight direction) is $70\text{\,}\mathrm{mm}$. Due to its large size, the
U.S. Geological Survey (USGS) scans the film in four parts with an overlap
between the successive scans. The USGS allows film scanning at
$7\text{\,}\mathrm{\SIUnitSymbolMicro m}$ or
$14\text{\,}\mathrm{\SIUnitSymbolMicro m}$ resolution. We stitch the
individual film scans together to generate the image of the whole film for
further processing. The stitching operation is performed automatically by
finding feature points in the overlapping regions of the scans and estimating
a transformation (2D rotation and translation) to align them. The four scans
are labeled by the USGS as a,b,c and d. We first stitch a with b and c with d,
then the parts ab and cd are merged together.

Figure 3: Schematic diagram of the film and the panoramic geometry reference
data in KH-4B missions, (source: [42])

.

##### Fiducial marks

During the image acquisition scan, additional data was exposed on the film to
allow reconstruction of the camera internal geometry and assessment of film
distortions. In the KH-4B missions the J-3 camera consisted of holes on the
rails supporting the film, which produce a sequence of spots with a
$1\text{\,}\mathrm{cm}$ spacing on both sides of the film length (Figure 3).
The center of format included two rail holes on both sides of the film. In
addition to the rail holes, there are Panoramic Geometry (PG) stripes running
along the length of the film on both sides that can be used to reconstruct the
panoramic geometry.

We reconstructed the position of the rail holes and the PG stripes along the
length of the film to assess the usability of this reference data for
correcting the film bending and distortions. A certain number of positions of
the PG stripes and the rail holes were extracted along the length of the film.
These positions were then interpolated for each pixel along the length of the
film. The position of the rail holes and the stripes in two Corona scanned
images are shown in Figure 4. We can observe that the positions of the rail
holes along the two sides of the film unexpectedly deviate from each other and
from the respective PG stripes. This may be due to the fact that the rail
holes are more towards the ends of the film where distortions are higher in
magnitude. This shows that the deviations in the rail holes are not
necessarily representative of the film distortions towards the center of the
image. On the other hand, the PG stripes which are comparatively closer to the
imaged area reveal similar bending effects on either side of the film.
Therefore, we use the location of the PG stripes to resample the image,
imposing that the PG stripes become straight and parallel to each other.

Figure 4: Position of the stripes and the rail holes above and below the
image. These positions are plotted with reference to the respective locations
corresponding to the center of the image.

In CoSP, we extract and align the exposed image area from the entire scan as
the scan also includes reference metadata and images from horizon cameras
(Figure 3), which are not relevant for further processing. We first find an
adaptive threshold to estimate exposed image pixels as the dark background has
a low gray-value. Based on the estimated exposed image pixels, we find the
centroid and principal components directions to rotate the image and align the
length and width of the film scan along the x,y-axes of the image. The PG
stripes are then extracted by finding the location of the edge pixels
extracted using a canny edge detector at top and bottom of the image. The
outliers are filtered out by median filtering the location of the edge pixels.
The images are then resampled so that the PG stripes become straight and
parallel to each other. Finally, we clip the $1.5\text{\,}\mathrm{cm}$ of the
film on both ends along the length of the film due to focus variations. This
focus anomaly is caused by the clearance required to avoid any collision
between the rotating scan-head and the film rollers [43]. The aft camera
images are rotated $180^{\circ}$ to align them with fore camera images. Note
that the rail holes and PG stripes are not available on the films of KH-4 and
KH-4A missions and the above-mentioned correction can not be performed for
images of these missions.

###  3.2 GCP Generation

In the earlier studies, the extraction of GCPs has been one of the most time
consuming component of the Corona imagery processing pipeline. These GCPs were
typically extracted by manually identifying features in Corona images and some
reference satellite imagery and extracting the corresponding elevation values
from a reference DEM. Automatic feature matching using SIFT [44] or other
traditional features detectors is challenging because of physical changes in
the scene due to long temporal difference. The radiometric differences and
variations in cast shadows in high relief terrain makes feature matching
particularly difficult. Recent progresses in deep learning based feature
matching has shown promising results in matching features with high
perspective and radiometric difference. Such detectors have already been used
for matching historical images [45, 46].

We use SuperGlue [47] a deep learning based architecture for feature
extraction and matching to compute point matches between Corona imagery and
reference satellite imagery. The SuperGlue architecture consists of a
Convolutional Neural Network (CNN) based interest point detector and
descriptor called SuperPoint [48], which uses VGG [49] style encoder. It uses
a graph neural network to perform matching of the features extracted by the
CNN and implicitly learns their spatial relationship within one image and
across a pair with the help of the attention mechanism. The attention
mechanism serves as a spatial filter and is learned end-to-end. The SuperGlue
like any other CNN based architecture is typically trained on relatively low
resolution images. The default image resolution for SuperGlue is $640\times
480$ and we found in our experiments that it provides reasonable results up to
three times the default resolution i.e, $1920\times 1440$. The Corona image
has a much larger size i.e. $\sim 106000\times 8000$. So, we divide it in to
patches of size $1920\times 1440$ and extract similar geographic location
patch from references satellite imagery to match the feature points. As the
resolution of KH-4B images is around 2-3 $\text{\,}\mathrm{m}$, we use
reference satellite imagery at similar spatial resolution to match feature
points. We used a mosaic of orthorectified PlanetScope images [50] with
$3\text{\,}\mathrm{m}$ spatial resolution for feature matching with the Corona
imagery.

As high resolution reference satellite imagery may not be available for larger
areas (here, we refer to 1-5$\text{\,}\mathrm{m}$ spatial resolution as high
resolution imagery), we also experimented with feature matching on freely
available medium resolution satellite imagery like Landsat-7 ETM+. When using
medium resolution imagery, we resample the Corona imagery to have approx.
similar spatial resolution as the ETM+ panchromatic band, i.e.,
$15\text{\,}\mathrm{m}$. Specifically, we divide Corona image into tiles of
$\sim 10600\times 8000$ and resize it to $1920\times 1440$ for input to
SuperGlue.

We use the geographic locations of the image corners specified in the Corona
image metadata to extract overlapping tiles from the Corona and reference
satellite imagery. The locations of the image corners are derived from the
satellite orbit, and can deviate several kilometers from their real positions
in a typical scenario or up to $100\text{\,}\mathrm{km}$ in certain cases. If
this deviation is larger than the size of the tile, feature matching will not
work as the tiles from Corona and reference satellite imagery will have no
overlap. Therefore, we adopt a coarse-to-fine strategy. We first find the
overlapping zones with low resolution images over an extended area. Then we
use the refined estimates of the Corona footprint to recover precise
correspondences in a localised area using higher resolution images.

###  3.3 Camera Model

Figure 5: The panoramic image coordinate system and the angle $\alpha$.

The scene is exposed to the film through a slit, which rotates along with the
lens. The scan angle $\alpha$ is the angle of the slit and is given by

$\alpha=\frac{x_{p}}{f}$ (1)

where $f$ is the focal length and $x_{p}$ is the panoramic photo coordinate
along the length of film (Figure 5). The exterior orientation parameters are
modeled as first order time dependent parameters

$X_{0t}=X_{0}+X_{01}\cdot t,$ (2)

$Y_{0t}=Y_{0}+Y_{01}\cdot t,$ (3)

$Z_{0t}=Z_{0}+Z_{01}\cdot t,$ (4)

$\omega_{0t}=\omega_{0}+\omega_{01}\cdot t,$ (5)

$\phi_{0t}=\phi_{0}+\phi_{01}\cdot t,$ (6)

$\kappa_{0t}=\kappa_{0}+\kappa_{01}\cdot t.$ (7)

Here, $X_{01},Y_{01},Z_{01}$ represent the motion of the camera and
$\omega_{01},\phi_{01},\kappa_{01}$ represent the angular motion along the
three axis during the panoramic scan with $t\in[0,1]$. Considering that the
width of the slit to be small, the geometry of the image formed at any
position of the slit is given by,

$\begin{bmatrix}0\\\ y_{p}+y_{IMC}\\\
-f\end{bmatrix}=sR_{\alpha}R\begin{bmatrix}X-X_{0t}\\\ Y-Y_{0t}\\\
Z-Z_{0t}\end{bmatrix},$ (8)

where $R_{\alpha}$ is,

$R_{\alpha}=\begin{bmatrix}cos\alpha&0&sin\alpha\\\ 0&1&0\\\
-sin\alpha&0&cos\alpha\end{bmatrix},$ (9)

and the IMC term is given as,

$y_{IMC}=-\frac{Vf}{H\delta}sin(\alpha)cos(\omega_{0}),$ (10)

where, $V$ is the satellite velocity, $H$ is altitude of the satellite and
$\delta$ is the scan angular velocity. The image motion caused by a moving
platform is directly related to V/H ratio [51, 52, 41]. Therefore, as the
orbital altitude varies the observed image motion also varies. In J-3 camera,
the image motion was compensated by rotating the lens in the direction
opposite to the image motion, while in the earlier J-1 system the lens was
moved along the image y-axis to compensate image motion [6]. The image motion
compensation produces an S-shape distortion across the image length [52].

By rearranging the terms of Equation (8) we obtain,

$R_{\alpha}^{T}\begin{bmatrix}0\\\ y_{p}+y_{IMC}\\\
-f\end{bmatrix}=sR\begin{bmatrix}X-X_{0t}\\\ Y-Y_{0t}\\\
Z-Z_{0t}\end{bmatrix},$ (11)

and finally,

$\begin{bmatrix}fsin\alpha\\\ y_{p}+y_{IMC}\\\
-fcos\alpha\end{bmatrix}=sR\begin{bmatrix}X-X_{0t}\\\ Y-Y_{0t}\\\
Z-Z_{0t}\end{bmatrix}.$ (12)

Let,

$N_{x}=r_{11}(X-X_{0t})+r_{12}(Y-Y_{0t})+r_{13}(Z-Z_{0t})$ (13)
$N_{y}={r_{21}(X-X_{0t})+r_{22}(Y-Y_{0t})+r_{23}(Z-Z_{0t})}$ (14)
$N_{z}=r_{31}(X-X_{0t})+r_{32}(Y-Y_{0t})+r_{33}(Z-Z_{0t})$ (15)

where $r_{ii}$ are the elements of the rotation matrix $R$. Division of
$1^{st}$ and the $2^{nd}$ row with the last row cancels out the scale factor
and gives the following equations,

$tan\alpha=-\frac{N_{x}}{N_{z}},$ (16)

$y+y_{IMC}=-fcos\alpha\frac{N_{y}}{N_{z}},$ (17)

The final image panoramic coordinates $(x_{p},y_{p})$ being,

$x_{p}=ftan^{-1}(-\frac{N_{x}}{N_{z}})$ (18)

$y_{p}=-y_{IMC}-fcos\alpha\frac{N_{y}}{N_{z}}$ (19)

There are 13 unknown camera parameters in the model: 6 parameters for camera
position $X_{0},Y_{0},Z_{0},X_{01},Y_{01},Z_{01}$, six parameters for camera
orientation $\omega_{0},\phi_{0},\kappa_{0},\omega_{01},\phi_{01},\kappa_{01}$
and one parameter for the IMC i.e. $\frac{V}{H\delta}$. We keep the focal
length fixed in the adjustment due to the correlation with the camera altitude
parameter [53]. The approximate values for $\omega_{0}$ were set to
$\\{-15^{\circ},15^{\circ}\\}$ for the fore and aft looking cameras
respectively, $\kappa_{0}$ as $0^{\circ}$, camera altitude $Z_{0}$ is set to
$170\text{\,}\mathrm{km}$ and the time dependent parameters are initialized as
zero. These parameters are then optimized in bundle adjustment using GCPs as
well as tie points between the stereo pairs. The bundle adjustment is
performed in an Earth Centered Cartesian frame.

###  3.4 Stereo Rectification

We follow the stereo rectification algorithm given in [54] and implemented in
MicMac [55]. This epipolar rectification is generic in the sense that it makes
no assumption on the image pair camera model, as long as the image-to-ground
projection functions are smooth. It is therefore suitable to uncommon image
geometries such as that of the Corona images. The algorithm consists of four
stages: (1) estimation of the global epipolar directions in either image of
the stereo-pair (see Figure 6 (top)); (2) generation of a set of virtual
correspondences between images using the known projection functions; (3)
estimation of the rectifying polynomials in a rotated image coordinate frame
where the epipolar curves are approximately horizontal; (4) computation of the
final rectifying functions as a composition of the rotations computed in step
1 and the polynomials computed in step 3. If the camera geometry and the
projection functions are unknown (e.g., due to lack of metadata or when
handling challenging camera geometries such as Corona), the image
correspondences can be replaced with image features extracted with image
processing algorithms such as SIFT [44] or SuperGlue [47]. Note that to obtain
unambiguous rectifications using this variant, the 3D scene should not be
flat. For more details, readers are referred to the original publication [54].
Here, we use a 4th degree polynomial to approximate the epipolar curves.

Figure 6: Epipolar curves in the native geometry of a Corona stereo pair (top)
and the resulting epipolar rectification (bottom). These curves are
exaggerated for visualization. Note that the curves are globally vertical and
they will be rotated by $\approx 90^{\circ}$ prior to estimating the
rectifying polynomials.

###  3.5 DEM Generation and DEM Coregistration

We compute the dense correspondences between the stereo rectified image pairs
with an implementation of the Semi Global Matching (SGM) [56]. SGM
approximates the 2D disparity cost as 1D cost along several paths, and has
been widely used for DSM generation from aerial and satellite imagery.

Once the dense point correspondences and the camera parameters are known, 3D
triangulation of the corresponding points is performed. Rearranging Eq. 11.
gives the following two equations for each image observation:

$\begin{bmatrix}R_{1}+tan(\alpha)R_{3}\\\
R_{2}+\frac{y+y_{IMC}}{fcos(\alpha)}R_{3}\end{bmatrix}\begin{bmatrix}X\\\ Y\\\
Z\end{bmatrix}=\begin{bmatrix}R_{1}+tan(\alpha)R_{3}\\\
R_{2}+\frac{y+y_{IMC}}{fcos(\alpha)}R_{3})\end{bmatrix}\begin{bmatrix}X_{0t}\\\
Y_{0t}\\\ Z_{0t}\end{bmatrix}.$ (20)

Here, $R_{1},R_{2}$ and $R_{3}$ are the 1st, 2nd and 3rd rows of the rotation
matrix respectively. Given the four equations, a least squares minimization is
performed to estimate the 3D coordinates for each point. This leads to a dense
3D point cloud, which is then interpolated into a raster DEM at the required
spatial resolution.

In practice, the GCPs computed automatically may not be well distributed
across the image. This leads to a systematic deviation between the Corona DEM
and the reference DEM. Even with well distributed GCPs, film distortions can
lead to systematic deviations between the two DEMs (see Figure 11). In
applications requiring DEM differencing, a fine coregistration between the two
DEMs is necessary to produce unbiased results. Therefore, to compensate for
these systematic deviations, we perform a fine coregistration of Corona DEM
and the reference DEM in a post-processing step. Concretely, we divide the
Corona DEM into smaller tiles and estimate a 3D affine transformation using
Least Square Matching (LSM) [57]. We use a tile size of approx.
$400\text{\,}{\mathrm{km}}^{2}$ and estimate the parameters of the
transformation individually for each tile. To ensure smooth transition between
neighbouring tiles we keep a certain amount of overlap between them.

We have implemented the camera model and bundle adjustment in MATLAB.
Implementation of dense matching, tie point extraction and epipolar resampling
is based on MicMac, while GCP extraction uses SuperGlue’s original inference
and evaluation script, and the DEM coregistration uses LSM module of OPALS
[58]. All these modules are executed from a script written in MATLAB.

##  4 Data and Study Areas

We applied and evaluated CoSP using three different case studies (see Table
2). We first focus on the evaluation of camera model and select a study site
for which well distributed GCPs can be automatically extracted. We searched
for an area with distinctive texture and minimal natural or anthropogenic
changes over the past decades. Furthermore, we restricted our search to areas
with shallower slopes as high relief terrain is associated with larger
elevation errors in the global DEMs, which we use in our study to extract the
GCPs’ elevation values. Our visual search revealed that area bordering
Mongolia and China is a suitable for this purpose. Therefore, we chose two
Corona stereo pairs over an area bordering China and Mongolia for the
evaluation of the camera model and the corresponding Corona DEM (Figure 7).
The two image pairs belong to mission no. 1112 and 1117 and the image
acquisition dates are 23-11-1970 and 30-05-1972, respectively. The images were
downloaded from the USGS Earth Explorer. The elevation range in the subject
area is between 300-2500 $\text{\,}\mathrm{m}$. We use SRTM 1 Arc-Second
Global DEM v3 (DOI: 10.5066/F7PR7TFT) for extracting the elevation values and
comparison with the Corona DEM. A mosaic of orthorectified PlanetScope imagery
of this area is used for extraction of GCPs using SuperGlue. It is pertinent
to mention that the GCPs also have limited accuracy as they are derived from
PlanetScope imagery (for planimetric coordinates) and $30\text{\,}\mathrm{m}$
resolution reference DEM (for elevation coordinates).

Figure 7: Corona scenes and corresponding footprints of the stereo pairs
DS1117-2071DF008 DS1117-2071DA014 and DS1112-1071DF023 DS1112-1071DA029 over
an area bordering Mongolia and China

Our second focus is the evaluation of the proposed pipeline in estimating
glacier change. To do that, we chose an area around Mount Everest in the
Himalaya because of the availability of high resolution Cartosat DEM
($10\text{\,}\mathrm{m}$) and Corona DEM computed using RSG. The latter allows
us to compare the DEMs computed from our rigorous model, and the approximated
variant of the Corona model implemented in RSG. The glaciers in the Himalaya
are of high interest due to accelerating mass loss [22]. We use a Corona
stereo pair from mission no. 1108 dated 18-12-1969 (Table 2) to create Corona
DEM. Here, we use $30\text{\,}\mathrm{m}$ ALOS PRISM DEM (AW3D) [59] as a
reference DEM because of the bias in SRTM DEM due to snow penetration of radar
signal.

Finally, to evaluate the potential of large scale mapping using CoSP, we use
12 successive stereo pairs over Central Pamir in Tajikistan from mission 1104,
dated 18-08-1968 (Table 2), covering a footprint of approx. $160\times 200$
km. This sequence of images are suitable for large scale glacier change
estimation due to late summer acquisition and cloud free coverage.

Table 2: Corona scenes used for each study site and the corresponding reference DEMs and satellite imagery used for GCPs Study | Corona Scenes | Date | Imagery | Reference
---|---|---|---|---
Site |  |  | for GCPs | DEM
Mongolia | DS1117-2071DF008 | 30-05-1972 | PlanetScope | SRTM
DS1117-2071DA014
DS1112-1071DF023 | 23-11-1970
DS1112-1071DA029
Everest | DS1108-2217DA070 | 18-12-1969 | Landsat-7 | ALOS
DS1108-2217DF064 | Cartosat
Pamir | DS1104-2169DF092- | 18-08-1968 | Landsat-7 | ALOS
DS1104-2169DF103
DS1104-2169DA098-
DS1104-2169DA109

##  5 Results

###  5.1 Accuracy of the Camera Model (Mongolia Data)

##### Accuracy on Check Points

SuperGlue extracted large number of features points between the Corona images
and the high resolution satellite imagery (see Figure 8). Fewer number of
points were detected in areas without texture as well as the areas with
surface changes. These feature points along with the corresponding height
values derived from the SRTM DEM form the GCPs and half of them are used as
check points to evaluate the accuracy of the camera model. The average
standard deviation ($\sigma_{0}$) of bundle adjustment is $\sim 2$ pixels for
both the image pairs (Table 3). The $\sigma_{0}$ shows an improvement of
$10-20\%$ when bending correction is applied. The RMSE error of the
checkpoints for the two image pairs are $[4.85,3.81,5.79]$ and
$[6.63,4.23,8.42]$ (m) (Table 3). As with the $\sigma_{0}$, the RMSE decreases
when film bending correction is performed.

Figure 8: Point matches between a Corona images and multi satellite PlanetScope imagery mosaic computed by SuperGlue. The grid type pattern of the feature points is due to tile based matching. Table 3: Accuracy of the bundle adjustment and RMSE of the check points for the two image pairs shown in Figure 7. The RMSE is computed in the corresponding UTM zone. Image Pair | $\sigma_{0}$ | RMSE X | RMSE Y | RMSE Z
---|---|---|---|---
(pixels) | (m) | (m) | (m)
Bending Correction |  |  |  |
DS1112-1071DF023 | 1.71 | 4.85 | 3.81 | 5.79
DS1112-1071DA029
DS1117-2071DF008 | 1.90 | 6.63 | 4.23 | 8.42
DS1117-2071DA014
Without Correction |  |  |  |
DS1112-1071DF023 | 1.89 | 5.84 | 4.09 | 9.92
DS1112-1071DA029
DS1117-2071DF008 | 2.43 | 7.39 | 4.75 | 15.11
DS1117-2071DA014

##### Residual systematic errors

In order to assess the residual distortions in the image, residuals are
computed from back projection of the GCPs and interpolated over a regular grid
(Figure 9). Image residuals show higher magnitude in some parts of the image
especially towards the edges of the film. These image residuals can mainly be
attributed to film distortions that occurred during the mission and long term
storage of the film.

Figure 9: Image residuals of the GCPs after resection are interpolated over a regular grid. These images have been corrected for film bending.  Table 4: The estimated camera parameters are given according to the geographic coordinate system with the position of the camera given as longitude, latitude and elevation. Parameters | DS1112- | DS1112- | DS1117- | DS1117-
---|---|---|---|---
1071DF023 | 1071DA029 | 2071DF008 | 2071DA014
$X_{0}$ (Lon.) | 96.24 | 96.39 | 98.26 | 98.04
$Y_{0}$ (Lat.) | 44.59 | 43.70 | 44.31 | 43.54
$Z_{0}$ (km) | 187.27 | 186.77 | 162.21 | 162.89
$X_{01}$ (km) | 0.06 | 1.20 | -1.53 | 0.41
$Y_{01}$ (km) | -2.75 | -2.86 | -2.59 | -3.10
$Z_{01}$ (km) | -0.41 | -0.03 | 0.15 | -0.14
$\omega_{0}$ (deg) | -15.20 | 15.72 | -15.52 | 14.68
$\phi_{0}$ (deg) | -1.56 | 1.46 | 3.71 | -2.25
$\kappa_{0}$ (deg) | 5.69 | 5.74 | -10.19 | -10.64
$\omega_{01}$ (deg) | 0.83 | 0.94 | 0.76 | 1.09
$\phi_{01}$ (deg) | -0.04 | -0.03 | -0.51 | -0.23
$\kappa_{10}$ (deg) | 0.00 | -0.23 | 0.33 | 0.02
$\frac{V}{H\delta}$ | 0.0025 | -0.0001 | -0.0002 | 0.0002

##### Camera parameters

Overall, the estimated camera parameters correspond relatively well to the
expected values (see Table 4). The estimated position and the orientation of
the satellite is consistent with the orbital parameters of the respective
acquisitions and stereo geometry. The main component of the estimated
satellite motion during the scan is along north direction i.e. $Y_{01}$, whose
estimated value ranges from 2.59-3.10 $\text{\,}\mathrm{km}$. The expected
value of motion along flight direction for a scanning time of 0.36 sec is
$2.8\text{\,}\mathrm{km}$ (velocity =
$7.7\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$). The expected values of
other time dependent orientation parameters will be close to zero if we assume
a stable satellite attitude during the acquisition.

The specifications of the J-3 camera states that the IMC in KH-4B missions was
implemented by rotating the camera in a direction opposite to the flight
direction, while in the earlier Corona missions IMC was implemented by
translating the lens relative to the film [6]. This may explain why: 1) The
estimated value of the IMC term $\frac{V}{H\delta}$ is significantly smaller
than the expected value of 0.014 (V=
$7.7\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$, H =
$170\text{\,}\mathrm{km}$, $\delta$=
$3.3\text{\,}\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$) as both image pairs are
from KH-4B missions and 2) A value between 0.76-1.09 deg for the parameter
$\omega_{01}$ because the rotation to compensate camera motion will be
observed in $\omega_{01}$. Given a scan time of 0.36 sec, the
$7.7\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ velocity of the satellite
will lead to a change of approx. 0.9 deg angle for an object point at the
center of the format. Thus, the estimated $\omega_{01}$ is consistent with the
rotation required to compensate the image motion. To further assess the effect
of IMC mechanism on the estimated parameters, we estimated camera parameters
for different KH-4A scenes using GCPs generated from the scheme introduced in
Section 3.2. In majority of the images, the estimated IMC term has a value
closer to the value of 0.014 (Table 5).

Table 5: Estimated time dependent orientation parameters and IMC for KH-4A scenes Image | $X_{01}$ | $Y_{01}$ | $Z_{01}$ | $\omega_{01}$ | $\phi_{01}$ | $\kappa_{01}$ | $v/H\delta$
---|---|---|---|---|---|---|---
| (km) | (km) | (km) | (deg) | (deg) | (deg) |
DS1024- | -0.187 | -3.909 | -0.652 | -0.020 | -0.202 | 0.050 | 0.016
1038DF096
DS1024- | 0.037 | -3.167 | -2.307 | 0.028 | -0.06 | 0.019 | 0.012
1038DF095
DS1049- | 0.463 | -3.181 | 0.453 | -0.15 | -0.002 | 0.014 | 0.017
2119DF055

###  5.2 Epipolar Resampling and the Accuracy of the DEM (Mongolia Data)

The y-parallax in the epipolar resmapled images is computed using 2D
correlation algorithm implemented in MicMac [60]. The SD of the y-parallax for
the two rectified stereo pairs is 0.89 and 0.99 pixels. The y-parallax over
majority of the area is less than $\pm 1$ pixels (Figure 10). However, certain
regions with higher y-parallax of around 5 pixels are also observed. This can
be due to the distortions in the film that affect locally certain part of the
images. Such systematic errors correlate with the image residuals as shown in
Figure 9.

Figure 10: The two image pairs and the corresponding y-parallax in the
epipolar resampled image pairs

The NMAD of elevation differences between the two Corona DEMs and the SRTM DEM
($30\text{\,}\mathrm{m}$) are $7.26\text{\,}\mathrm{m}$ and
$6.69\text{\,}\mathrm{m}$ (Figure 11). The elevation differences show
systematic deviations of up to $25\text{\,}\mathrm{m}$ in the respective DEMs.
To eliminate these systematic deviations in the Corona DEM, we further align
tiles of the Corona DEMs with the reference DEMs by estimating a 3D affine
transformation, as explained in Section 3.5. This alignment of the DEMs
reduces the NMAD of the elevation differences to $3.32\text{\,}\mathrm{m}$ and
$4.15\text{\,}\mathrm{m}$ for the two Corona DEMs (Figure 11). The size of
each tile is approx. $400\text{\,}{\mathrm{km}}^{2}$, which is large enough to
avoid fitting the 3D affine transform to the temporal surfaces changes that we
are interested in estimating.

Figure 11: Elevation Differences between the Corona DEM and SRTM DEM for the
two stereo pairs. The DEMs are divided in to patches and coregistered, which
reduces the systematic errors in the Corona DEM.

###  5.3 Case Study: Glacier Change Estimation (Everest Data)

In contrast to the scenes used for evaluating the camera model, the images
over glacierized regions show significant seasonal and long-term variations.
Due to difference in the Sun position and the cast shadows; the automatic
feature matching across the mountainous terrain can be quite challenging. Even
with these challenges, SuperGlue matched some points between Corona imagery
and Landsat-7 ETM+ mosaic (Figure 12). However, no matches were found in the
southwest part of the scene consisting of hilly forested terrain due to low
texture and changes in the terrain and shadows.

Figure 12: Feature points matched between 1969 Corona KH-4 and Landsat-7
panchromatic image mosaic (year 2000) of the Everest region. Due to large area
covered by the images, three locations are zoomed for visual comparison.

As no GCPs were available in one section of the image, the resulting DEM over
that area show large elevation differences in comparison with the ALOS DEM
$30\text{\,}\mathrm{m}$ (Figure 13). The Corona DEM is then coregistered with
the ALOS DEM using smaller tiles, which removed the majority of the systematic
elevation differences between the Corona and the ALOS DEM. The DEM
coregistration is performed using the stable terrain i.e. the area without the
glaciers. To filter out the glacierized part of the scene, we used the glacier
polygons from the work of [22]. The Randolph Glacier Inventory (RGI) [61]
polygons were manually modified to match glacier extent from 1969. The
elevation differences show errors caused by textureless snow/ice cover in the
accumulation zone of the glacier and steep rockwalls. Further large deviations
are due to the presence of clouds in the Corona imagery (Figure 13).

Figure 13: Corona stereo pair and the elevation differences between the Corona
DEM and ALOS DEM ($30\text{\,}\mathrm{m}$). a) and b) show the Corona image
pair c) shows the elevation difference with the ALOS DEM d) shows the
elevation differences after tile based coregistration of the Corona DEM and
the ALOS DEM e) shows the shaded Corona DEM.

The highest spatial resolution in the KH-4 missions is around
$1.8\text{\,}\mathrm{m}$ at the center of the image and around
$3\text{\,}\mathrm{m}$ towards the edges of the film. Hence, DEMs with
relatively high spatial resolution can be generated from the Corona images.
Here, we compute a $10\text{\,}\mathrm{m}$ Corona DEM of the Everest region
and compare it with a $10\text{\,}\mathrm{m}$ Cartosat DEM of 15/12/2018
(Figure 14) from the work of [22]. The elevation differences between the
Corona and Cartosat DEMs show a decrease in the elevation of the glaciers and
thus a negative glacier mass balance. The elevation differences over the
stable terrain shows that the two DEMs are well coregistered. A visual
inspection of the shaded DEMs show similar level of details in the Corona and
Cartosat DEMs over the well textured part of the image. Though the Corona DEM
has higher noise level. The accumulation region of the glaciers’ with no
texture result in large elevation differences (Figure 14). Note that we
generated spatially complete Corona DEM but in practice such unreliable
elevation values can be filtered out to reduce the chances of introducing bias
in the estimation of glacier volume change.

This $10\text{\,}\mathrm{m}$ Corona DEM of Everest region, generated using
CoSP is compared with the Corona DEM generated using RSG as used in [22]. The
NMAD of the elevation differences between CoSP Corona DEM and Cartosat DEM
over the stable terrain is $10.18\text{\,}\mathrm{m}$, which is slightly
better than the NMAD of $12.41\text{\,}\mathrm{m}$ obtained by Corona DEM
generated using RSG (Figure 14). The elevation differences over the textured
part of the glaciers are very similar for CoSP and RSG Corona DEMs. However,
the pattern of outliers in textureless accumulation regions of the glaciers
show a very different pattern, which is possibly due to differences in
regularization parameters of the dense matching algorithms used in CoSP and
RSG. After performing the outlier removal and gap filling, the elevation
difference from 1969 to 2018 over the glacierized area shown in Figure 14 is
$-13.55\pm$3.91\text{\,}\mathrm{m}$$ for CoSP and
$-10.26\pm$4.13\text{\,}\mathrm{m}$$ for RSG computed Corona DEMs. The
elevation change for the Khumbu Glacier is
$-16.96\pm$3.91\text{\,}\mathrm{m}$$ for CoSP and
$-17.20\pm$4.13\text{\,}\mathrm{m}$$ for RSG generated Corona DEMs.

Figure 14: a) Hillshade of $10\text{\,}\mathrm{m}$ Corona DEM b) Hillshade of
$10\text{\,}\mathrm{m}$ Cartosat DEM (Date: 15/12/2018) c) Elevation
differences between Cartosat and Corona DEM computed using the proposed
methodology d) Elevation differences between Cartosat and Corona DEM computed
using RSG [22] e) Histogram of elevation differences given in [c] f) Corona
orthophoto of the Everest Region g) Histogram of elevation differences given
in [d]

###  5.4 Large scale Mapping using CoSP (Pamir Data)

In order to show the applicability for large scale mapping, we apply the
proposed methodology on a sequence of 12 Corona stereo pairs over Central
Pamir (Figure 15). We applied all processing steps and differenced the
resultant DEM with the $30\text{\,}\mathrm{m}$ ALOS DEM. The DEM differences
show that CoSP works on a larger scale. The systematic elevation errors are
low in magnitude over the entire area. The large number of surge type glaciers
with either strong elevation loss in the middle reaches of the glacier and an
elevation gain at its tongues or vice versa [62], can be well identified in
the DEM differences.

This sequence of Corona imagery also shows a potential complication in large
scale processing. The consecutive images of the sequence have an overlap of
approx. 6% at the center of the format, which is less than the specified
overlap of approx. 8%. While the the fore and aft looking cameras have an
overlap of approx. 90%. As a result there are gaps in between the consecutive
stereo pairs with greater voids towards the center of the image. These voids
can be filled by considering the small overlap with the adjacent aft image of
the sequence. However, the coregistration of the DEMs with small width (e.g.
$3\times 200$ km) can be challenging in presence of unstable terrain.

Figure 15: DEM differences between the Corona DEM computed using 12 Corona
stereo pairs and the ALOS DEM over Pamir region.

##  6 Discussion

###  6.1  Camera Model

The results show that a SD ($\sigma_{0}$) of better than $2$ pixels in bundle
adjustment can be obtained with a rigorous Corona camera model consisting of
time dependent collinearity equations. The image residuals show a pattern,
which can be attributed to film distortions due to long term storage and large
temperature variations during the mission as also observed for Hexagon KH-9
films [63, 31]. Furthermore, film bowing/buckling, misalignment in the
components of the camera and atmospheric effects may also impart distortions
in the image. It is pertinent to mention that these instruments were
continuously upgraded and technical problems have occurred during the Corona
missions, so there can be other factors associated with the residuals
distortions observed in the Corona images.

The rigorous Corona camera model used here is based on the work of [34]. This
camera model is more flexible as it includes six additional parameters for
time dependent modelling of exterior orientation parameters. The estimated
camera parameters are consistent with the expected values and the variation in
their values is due to the systematic errors caused by film distortions. The
estimated camera parameters show that the difference in the IMC mechanism in
KH-4A and KH-4B is observed in the estimated camera parameters and the
inclusion of IMC term in the camera model may not be necessary for KH-4B
cameras.

The height RMSE of check points obtained here i.e. $5.79\text{\,}\mathrm{m}$
and $8.42\text{\,}\mathrm{m}$ for the two Corona stereo pair (Table 3) is
higher than the height RMSE of check points of around $4\text{\,}\mathrm{m}$
reported in [34]. This decrease in the accuracy is expected due to the
presence of film distortions, when considering the entire Corona stereo
coverage (area $>$$3000\text{\,}{\mathrm{km}}^{2}$) in comparison to using
certain section of the image (area of $561\text{\,}{\mathrm{km}}^{2}$) as done
in [34]. As Figures 9 and 10 suggest, residuals distortions may have a higher
magnitude in certain parts of the image and considering larger image area may
decrease the accuracy of check points due to these distortions. The patch
based fine coregistration with the reference DEM improves the accuracy and
reduces the SD ($\sigma$) of elevation differences to around
$4\text{\,}\mathrm{m}5\text{\,}\mathrm{m}$ (Figure 11). This shows that
accuracy similar to [34] can be obtained over the entire stereo coverage as
fine coregistration compensates the effect of the distortions along the length
of the film to a certain extent.

[37] have reported a height RMSE of $12.34\text{\,}\mathrm{m}$ for
checkpoints, which is higher than RMSE obtained in our work. However, they
used a scan with $12\text{\,}\mathrm{\SIUnitSymbolMicro m}$ pixel size, which
reduces accuracy in comparison to $7\text{\,}\mathrm{\SIUnitSymbolMicro m}$
pixel size used in this work. The camera model in [37] was based on time
dependent collinearity equations but they only modeled motion along flight
direction as an additional time dependent parameter.

These results further show an improvement in accuracy in comparison to the
work of [27], who have reported an RMSE of around $2.5$ pixels for RSG Corona
camera model. In addition, the statistics of the elevations differences with a
$10\text{\,}\mathrm{m}$ Cartosat DEM also show slight improvement in
comparison to the Corona DEM generated using RSG in the work of [22] (see
Figure 14). This demonstrates that CoSP can achieve accuracy similar to the
best reported accuracy in the previous work on Corona imagery, while
considering the entire Corona stereo image coverage.

###  6.2 Corona Data Limitations

####  6.2.1 Film Reference Data

The reference metadata exposed on the film is important for reconstruction of
the panoramic geometry and correction of film distortions. However, as shown
earlier, the PG rail holes show unexpected deviations and are not feasible for
correction of film distortions or reconstruction of film geometry. The KH-4
and KH-4A imagery doesn’t have PG stripes and PG rail holes. Furthermore, we
have observed several KH-4B images in which PG stripes are not fully visible.
So, bending correction as performed here may not be applicable to all Corona
images, which will have an effect on the expected accuracy.

####  6.2.2 Low Texture and Film Saturation

One limitation of Corona imagery, which is especially significant for glacier
mass balance studies is the low texture and contrast over snow covered part of
the glaciers mainly caused by the saturation due to light exposed on the film.
This issue of textureless area has a greater impact in snow covered parts of
the glaciers especially in the accumulation regions, where it may be difficult
to estimate reliable elevation data [21]. This problem is further exacerbated
by lack of cloud free images in the late summer to early winter time period.
Although the problem of low image texture over snow covered area is common to
all optical sensors, the spatial resolution, wavelength range, radiometric
resolution, dynamic range and image acquisition time can significantly impact
the quality of DEM over the snow covered areas.

####  6.2.3 Inconsistent Overlap

There are two overlaps to consider, when creating DEMs from a sequence of
Corona stereo images: First is the overlap between the successive frames of a
Corona camera and the second overlap is between the stereo pair. The overlap
between the successive frames is specified to be approx. 8% at the center of
the format and is directly related to the camera scan rate, which in turn
depends on satellite velocity to height ratio [6]. The second overlap is
between the stereo pair i.e. images of an area from the fore and aft looking
cameras. This overlap is typically around 90% but varies over different
acquisitions. The overlap between the successive images is rather small, which
results in a small overlap between the DEMs from the successive stereo pairs
[21, 40]. Hence, deviations from the specified scanning rate or attitude
perturbations can cause gaps in the DEMs as seen in the image sequence over
Pamir, where the overlap between the successive images was approx. 6%.

####  6.2.4 Scanning Artifacts

The scanned imagery delivered by USGS has been reported to contain scanning
artifacts [39, 64, 65, 30]. These artifacts occur due to error in the scanner
calibration. After processing KH-4 imagery of different missions, we have also
observed such artifacts in the scanned Corona images. These blocking artifacts
can be observed as a misalignment in the overlapping stitched parts of the
scans. To quantify the misalignments, we first align the individual scans by a
rotation and translation using the tie points extracted in the overlapping
region of the neighboring scans. Then dense image matching is performed on the
aligned images. The result of the matching of the three image parts along $x$
and $y$ directions are shown in Figure 16. As the Corona images have been
available since long, these images may have been scanned with a different
scanner or with varying scanner calibration. As a consequence, these artifacts
show varying magnitude in different images. Such blocking artifacts are also
visible in the DEM difference shown in Figures 10 and 11.

Figure 16: Scanning artifacts observed in the aligned overlapping parts of the
scans. The shift along the x and y axis are estimated using 2D cross
correlation in MicMac.

###  6.3 GCP Generation

The main advantage of CoSP is the automated processing of Corona imagery,
which enables large scale mapping using Corona imagery. Automatic GCP
generation is an essential part of this pipeline. Our main focus of
application has been on the glaciers of high mountain Asia and GCP generation
using SuperGlue has largely been successful in generating automatic GCPs over
such terrain. However, automatic GCP generation may still be challenging over
scenes dominated by large forested areas, agricultural land, large water
bodies, deserts and urban areas, which may have changed significantly over a
period of time. Performing multi image bundle adjustment for a sequence of
Corona images can potentially resolve the issue of unavailability of GCPs in
some images of the sequence if GCPs can be extracted for other images of the
sequence. Another possible solution is to utilize DEMs to matching features
for georeferencing of the historical imagery as done in [46, 31].

###  6.4 DEM Coregistration

Here, we used a tile based approach to coregister the Corona DEM with a
reference DEM, which appears to work well in the above test cases. However,
problems in coregistration can appear if large parts of the scene contain
unstable terrain, outliers or data voids. Thus, a coregisteration using larger
area or whole DEM will be more suitable. A linear transformation in 3D was
only sufficient to coregister small parts (approx. 10%) of the scene and a
nonlinear transformation will be required for coregistration of larger areas.

There are other approaches to DEM cogregistration, which apply corrections
functions to a raster DEM. Such approaches have widely been used for
coregistration of DEMs for quantification of glacier elevation change [66, 29,
67, 68]. These approaches typically consist of multiple components such as
correction of DEM shift, elevation dependent bias correction, sensor specific
corrections such as the effect of jitter in ASTER DEMs and correction of
remaining systematic elevation differences by polynomial fittting. Such
approaches may also be investigated for improvement of the DEM coregistration
strategy in CoSP.

##  7 Conclusions

This work presented CoSP: A pipeline for automated processing of Corona KH-4
series imagery. A rigorous camera model with modified collinearity equations
consisting of time dependent exterior orientation parameters was used to model
the panoramic image acquisition geometry of the Corona cameras. A deep
learning-based feature extractor was used to automatically match feature
points between Corona imagery and medium to high resolution satellite imagery.
These feature points were used as GCPs to evaluate the accuracy of the camera
model. The results show that overall a $\sigma_{0}$ better than 2 pixels can
be achieved, while considering the entire Corona image. The estimated camera
parameters are consistent with the known parameters as well as the IMC
mechanism on the Corona camera system. Systematic image residuals up to six
pixels have been observed in different Corona scenes. These errors are perhaps
due to film distortions that have occurred during the mission and the storage
in the 50-year period. Consequently, the resulting Corona DEMs exhibit
systematic elevation differences of up to $25\text{\,}\mathrm{m}$, which could
be reduced by a tile-based coregisteration of the Corona DEM. A tie point
based epipolar resampling of the Corona images showed residual y-parallax of
$\pm 1$ across majority of the rectified image pair. The accuracy of the
derived 3D data was further improved by performing bending correction using
the PG stripes exposed on the film. Finally, it was shown that the proposed
methodology can be used on Corona scenes over high relief and glacierized
terrain, which enables calculation of multi-decadal glacier mass balance over
large areas.

## References

  * McDonald [1997] Robert A McDonald.  _Corona between the Sun & the Earth: the first NRO reconnaissance eye in space_.  American Society for Photogrammetry and Remote Sensing, 1997.
  * Galiatsatos [2004] Nikolaos Galiatsatos.  _Assessment of the CORONA series of satellite imagery for landscape archaeology: A case study from the Orontes valley, Syria_.  PhD thesis, Durham University, 2004.
  * Altmaier and Kany [2002] Angela Altmaier and Christoph Kany.  Digital surface model generation from Corona satellite images.  _ISPRS Journal of Photogrammetry and Remote Sensing_ , 56(4):221–235, 2002.
  * Dashora et al. [2007] Ajay Dashora, Bharat Lohani, and Javed N. Malik.  A repository of Earth resource information – Corona satellite programme.  _Current Science_ , 92(7):926–932, 2007.
  * NRO [1997] NRO.  Index to the declassified collection of CORONA, ARGON and LANYARD records declassified by the National Reconnaissance Office on november 26, 1997, 1997.  URL https://www.nro.gov/FOIA/Major-NRO-Programs-and-Projects.  accessed 18-Oct-2021.
  * NRO [6/30/1967] NRO.  Description and operational manual: The J-3 panoramic camera system.  Technical report, In National Reconnaissance Office Data Book, 6/30/1967.
  * NRO [1/28/1964] NRO.  Mural system performance evaluation report.  Technical report, In National Reconnaissance Office Data Book, 1/28/1964.
  * Madden [1996] FJ Madden.  The Corona camera system: Itek’s contribution to world stability, self-published.  1996\.
  * NRO [1970] NRO.  Corona technical information, volume I.  Technical report, In National Reconnaissance Office Data Book, 1970.
  * [10] NASA JPL.  Mission and spacecraft library, jpl: Corona.  URL space.jpl.nasa.gov/msl/Programs/corona.html.  accessed 18-Oct-2021.
  * Bayram et al. [2004] B Bayram, H Bayraktar, C Helvaci, and U Acar.  Coastline change detection using CORONA, SPOT and IRS 1D images.  _International Archives of Photogrammetry and Remote Sensing_ , 35(B7):437–441, 2004.
  * Saruulzaya et al. [2016] Adiya Saruulzaya, Mamoru Ishikawa, Yamkhin Jambaljav, et al.  Thermokarst lake changes in the southern fringe of Siberian permafrost region in Mongolia using Corona, Landsat, and ALOS satellite imagery from 1962 to 2007.  _Advances in Remote Sensing_ , 5(04):215, 2016\.
  * Noaje and Sion [2012] Ioan Noaje and Ion Gr Sion.  Environmental changes analysis in Bucharest City using Corona, SPOT HRV and IKONOS images.  _International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences_ , 39:B7, 2012.
  * Fowler and Fowler [2005] Martin JF Fowler and Yvonne M Fowler.  Detection of archaeological crop marks on declassified Corona KH-4B intelligence satellite photography of southern england.  _Archaeological Prospection_ , 12(4):257–264, 2005.
  * Watanabe et al. [2017] Nobuya Watanabe, Shinichi Nakamura, Bin Liu, and Ningyuan Wang.  Utilization of structure from motion for processing Corona satellite images: Application to mapping and interpretation of archaeological features in Liangzhu Culture, China.  _Archaeological Research in Asia_ , 11:38–50, 2017.
  * Bhambri et al. [2011] Rakesh Bhambri, Tobias Bolch, Ravinder Kumar Chaujar, and Subhash Chandra Kulshreshtha.  Glacier changes in the Garhwal Himalaya, India, from 1968 to 2006 based on remote sensing.  _Journal of Glaciology_ , 57(203):543–556, 2011\.
  * Schmidt and Nüsser [2012] Susanne Schmidt and Marcus Nüsser.  Changes of high altitude glaciers from 1969 to 2010 in the Trans-Himalayan Kang Yatze Massif, Ladakh, northwest India.  _Arctic, Antarctic, and Alpine Research_ , 44(1):107–121, 2012.
  * Robson et al. [2018] Benjamin A Robson, Christopher Nuth, Pål R Nielsen, Luc Girod, Marijn Hendrickx, and Svein Olaf Dahl.  Spatial variability in patterns of glacier change across the Manaslu range, Central Himalaya.  _Frontiers in Earth Science_ , 6:12, 2018.
  * Racoviteanu et al. [2015] Adina E Racoviteanu, Yves Arnaud, Mark W Williams, and William F Manley.  Spatial patterns in glacier characteristics and area changes from 1962 to 2006 in the Kanchenjunga–Sikkim area, Eastern Himalaya.  _The Cryosphere_ , 9(2):505–523, 2015.
  * Bolch et al. [2008] Tobias Bolch, Manfred Buchroithner, Tino Pieczonka, and André Kunert.  Planimetric and volumetric glacier changes in the Khumbu Himal, Nepal, since 1962 using Corona, Landsat TM and ASTER data.  _Journal of Glaciology_ , 54(187):592–600, 2008\.
  * Goerlich et al. [2017] Franz Goerlich, Tobias Bolch, Kriti Mukherjee, and Tino Pieczonka.  Glacier mass loss during the 1960s and 1970s in the Ak-Shirak range (Kyrgyzstan) from multiple stereoscopic Corona and Hexagon imagery.  _Remote Sensing_ , 9(3):275, 2017.
  * King et al. [2020] Owen King, Atanu Bhattacharya, Sajid Ghuffar, Alex Tait, Sam Guilford, Aurora C Elmore, and Tobias Bolch.  Six decades of glacier mass changes around Mt. Everest are revealed by historical and contemporary images.  _One Earth_ , 3(5):608–620, 2020.
  * Bolch et al. [2019] Tobias Bolch, N Rohrbach, S Kutuzov, BA Robson, and A Osmonov.  Occurrence, evolution and ice content of ice-debris complexes in the Ak-Shiirak, Central Tien Shan revealed by geophysical and remotely-sensed investigations.  _Earth Surface Processes and Landforms_ , 44(1):129–143, 2019.
  * Kääb et al. [2021] Andreas Kääb, Tazio Strozzi, Tobias Bolch, Rafael Caduff, Håkon Trefall, Markus Stoffel, and Alexander Kokarev.  Inventory and changes of rock glacier creep speeds in Ile Alatau and Kungöy Ala-Too, northern Tien Shan, since the 1950s.  _The Cryosphere_ , 15(2):927–949, 2021.
  * Bolch et al. [2011] Tobias Bolch, Juliane Peters, Alexandr Yegorov, Biswajeet Pradhan, Manfred Buchroithner, and Victor Blagoveshchensky.  Identification of potentially dangerous glacial lakes in the northern tien shan.  _Nat Hazards_ , 59:1691–1714, 2011.
  * Raj et al. [2013] K Babu Govindha Raj, SN Remya, and K Vinod Kumar.  Remote sensing-based hazard assessment of glacial lakes in Sikkim Himalaya.  _Current Science_ , pages 359–364, 2013.
  * Bhattacharya et al. [2021] Atanu Bhattacharya, Tobias Bolch, Kriti Mukherjee, Owen King, Brian Menounos, Vassiliy Kapitsa, Niklas Neckel, Wei Yang, and Tandong Yao.  High Mountain Asian glacier response to climate revealed by multi-temporal satellite observations since the 1960s.  _Nature Communications_ , 12(1):1–13, 2021.
  * Pieczonka et al. [2013] Tino Pieczonka, Tobias Bolch, Wei Junfeng, and Liu Shiyin.  Heterogeneous mass loss of glaciers in the Aksu-Tarim catchment (Central Tien Shan) revealed by 1976 KH-9 Hexagon and 2009 SPOT-5 stereo imagery.  _Remote Sensing of Environment_ , 130:233–244, 2013.
  * Pieczonka and Bolch [2015] Tino Pieczonka and Tobias Bolch.  Region-wide glacier mass budgets and area changes for the Central Tien Shan between~ 1975 and 1999 using Hexagon KH-9 imagery.  _Global and Planetary Change_ , 128:1–13, 2015.
  * Dehecq et al. [2020] Amaury Dehecq, Alex S Gardner, Oleg Alexandrov, Scott McMichael, Romain Hugonnet, David Shean, and Mauro Marty.  Automated processing of declassified KH-9 Hexagon satellite images for global elevation change analysis since the 1970s.  _Frontiers in Earth Science_ , 8:566802, 2020.
  * Maurer and Rupper [2015] Joshua Maurer and Summer Rupper.  Tapping into the Hexagon spy imagery database: A new automated pipeline for geomorphic change detection.  _ISPRS Journal of Photogrammetry and Remote Sensing_ , 108:113–127, 2015.
  * Maurer et al. [2019] Josh M Maurer, JM Schaefer, S Rupper, and A Corley.  Acceleration of ice loss across the Himalayas over the past 40 years.  _Science advances_ , 5(6):eaav7266, 2019.
  * King et al. [2019] Owen King, Atanu Bhattacharya, Rakesh Bhambri, and Tobias Bolch.  Glacial lakes exacerbate Himalayan glacier mass loss.  _Scientific Reports_ , 9(1):1–9, 2019.
  * Sohn et al. [2004] Hong-Gyoo Sohn, Gi-Hong Kim, and Jae-Hong Yom.  Mathematical modelling of historical reconnaissance Corona KH-4B imagery.  _The Photogrammetric Record_ , 19(105):51–66, 2004.
  * Fraser et al. [2006] Clive S Fraser, G Dial, and J Grodecki.  Sensor orientation via RPCs.  _ISPRS Journal of Photogrammetry and Remote Sensing_ , 60(3):182–194, 2006.
  * Tao and Hu [2001] C Vincent Tao and Yong Hu.  A comprehensive study of the rational function model for photogrammetric processing.  _Photogrammetric Engineering and Remote Sensing_ , 67(12):1347–1358, 2001.
  * Shin and Schenk [2008] Sung-Woong Shin and Tony Schenk.  Rigorous modeling of the first generation of the reconnaissance satellite imagery.  _Korean Journal of Remote Sensing_ , 24(3):223–233, 2008.
  * Beyer et al. [2018] Ross A Beyer, Oleg Alexandrov, and Scott McMichael.  The ames stereo pipeline: Nasa’s open source software for deriving and processing terrain data.  _Earth and Space Science_ , 5(9):537–548, 2018\.
  * Lauer [2019] Benjamin Lauer.  _Exploiting space-based optical and radar imagery to measure and model tectonic deformation in continental areas_.  PhD thesis, Université de Paris (2019), 2019.
  * Jacobsen [2020] Karsten Jacobsen.  Calibration and validation of Corona KH-4B to generate height models and orthoimages.  _ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences; volume 5, 1_ , 5(1):151–155, 2020.
  * Brown [1965] Earle B Brown.  V/H image motion in aerial cameras (image motion in focal plane of aerial camera resulting from forward motion of vehicle).  _Photogrammetric Engineering_ , 31:308–323, 1965.
  * NRO [1967] NRO.  The corona KH-4B camera system.  Technical report, In National Reconnaissance Office Data Book, 1967.
  * NRO [8/15/1968] NRO.  The non-specified area of j-3 format.  Technical report, In National Reconnaissance Office Data Book, 8/15/1968.
  * Lowe [2004] D. G. Lowe.  Distinctive image features from scale-invariant keypoints.  _International Journal of Computer Vision_ , 60(2):91–110, 2004.  doi:10.1023/B:VISI.0000029664.99615.94.
  * Maiwald and Maas [2021] Ferdinand Maiwald and Hans-Gerd Maas.  An automatic workflow for orientation of historical images with large radiometric and geometric differences.  _The Photogrammetric Record_ , 36(174):77–103, 2021.
  * Zhang et al. [2021] Lulin Zhang, Ewelina Rupnik, and Marc Pierrot-Deseilligny.  Feature matching for multi-epoch historical aerial images.  _ISPRS Journal of Photogrammetry and Remote Sensing_ , 182:176–189, 2021.
  * Sarlin et al. [2020] Paul-Edouard Sarlin, Daniel DeTone, Tomasz Malisiewicz, and Andrew Rabinovich.  Superglue: Learning feature matching with graph neural networks.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 4938–4947, 2020.
  * DeTone et al. [2018] Daniel DeTone, Tomasz Malisiewicz, and Andrew Rabinovich.  Superpoint: Self-supervised interest point detection and description.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , pages 224–236, 2018.
  * Simonyan and Zisserman [2014] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  _arXiv preprint arXiv:1409.1556_ , 2014.
  * Team [2018–] Planet Team.  Planet application program interface: In space for life on earth, 2018–.  URL https://api.planet.com.
  * Kawachi [1963] Donald Kawachi.  Image displacements on a tilted panoramic photograph.  _Optical Engineering_ , 1(3):010388, 1963.
  * Itek-Laboratories [1962] Itek-Laboratories.  Panoramic progress-part II. Technical analysis of panoramic photography.  _Photogrammetric Engineering_ , 28:99–107, 1962.
  * Zhou et al. [2020] Yilin Zhou, Ewelina Rupnik, Christophe Meynard, Christian Thom, and Marc Pierrot-Deseilligny.  Simulation and analysis of photogrammetric UAV image blocks—influence of camera calibration error.  _Remote Sensing_ , 12(1):22, 2020.
  * Deseilligny and Rupnik [2020] Marc Pierrot Deseilligny and Ewelina Rupnik.  Epipolar rectification of a generic camera.  working paper or preprint, October 2020.  URL https://hal.archives-ouvertes.fr/hal-02968078.
  * Pierrot-Deseilligny et al. [2014] Marc Pierrot-Deseilligny, D Jouin, J Belvaux, G Maillet, L Girod, E Rupnik, J Muller, M Daakir, G Choqueux, and M Deveau.  MicMac, apero, pastis and other beverages in a nutshell.  _Institut Géographique National_ , 2014.
  * Hirschmuller [2007] Heiko Hirschmuller.  Stereo processing by semiglobal matching and mutual information.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 30(2):328–341, 2007.
  * Ressl et al. [2011] Camillo Ressl, Norbert Pfeifer, and Gottfried Mandlburger.  Applying 3-D affine transformation and least squares matching for airborne laser scanning strips adjustment without GNSS/IMU trajectory data.  In _ISPRS workshop laser scanning_ , 2011.
  * Pfeifer et al. [2014] Norbert Pfeifer, Gottfried Mandlburger, Johannes Otepka, and Wilfried Karel.  OPALS–a framework for airborne laser scanning data analysis.  _Computers, Environment and Urban Systems_ , 45:125–136, 2014.
  * Tadono et al. [2014] T Tadono, H Ishida, F Oda, S Naito, K Minakawa, and H Iwamoto.  Precise global dem generation by alos prism.  _ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences_ , 2(4):71, 2014.
  * Rosu et al. [2015] Ana-Maria Rosu, Marc Pierrot-Deseilligny, Arthur Delorme, Renaud Binet, and Yann Klinger.  Measurement of ground displacement from optical satellite image correlation using the free open-source software MicMac.  _ISPRS Journal of Photogrammetry and Remote Sensing_ , 100:48–59, 2015.
  * Pfeffer et al. [2014] W Tad Pfeffer, Anthony A Arendt, Andrew Bliss, Tobias Bolch, J Graham Cogley, Alex S Gardner, Jon-Ove Hagen, Regine Hock, Georg Kaser, Christian Kienholz, et al.  The Randolph Glacier Inventory: A globally complete inventory of glaciers.  _Journal of Glaciology_ , 60(221):537–552, 2014\.
  * Goerlich et al. [2020] Franz Goerlich, Tobias Bolch, and Frank Paul.  More dynamic than expected: an updated survey of surging glaciers in the Pamir.  _Earth System Science Data_ , 12(4):3161–3176, 2020.
  * Surazakov and Aizen [2010] Arzhan Surazakov and Vladimir Aizen.  Positional accuracy evaluation of declassified hexagon kh-9 mapping camera imagery.  _Photogrammetric Engineering & Remote Sensing_, 76(5):603–608, 2010.
  * Galiatsatos et al. [2007] Nikolaos Galiatsatos, Daniel NM Donoghue, and Graham Philip.  High resolution elevation data derived from stereoscopic Corona imagery with minimal ground control.  _Photogrammetric Engineering & Remote Sensing_, 73(9):1093–1106, 2007.
  * Gheyle et al. [2011] Wouter Gheyle, Jean Bourgeois, Rudi Goossens, and Karsten Jacobsen.  Scan problems in digital Corona satellite images from USGS archives.  _Photogrammetric Engineering & Remote Sensing_, 77(12):1257–1264, 2011.
  * Nuth and Kääb [2011] C. Nuth and A. Kääb.  Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change.  _The Cryosphere_ , 5(1):271–290, 2011.
  * Girod et al. [2017] Luc Girod, Christopher Nuth, Andreas Kääb, Robert McNabb, and Olivier Galland.  MMASTER: improved ASTER DEMs for elevation change monitoring.  _Remote Sensing_ , 9(7):704, 2017.
  * Hugonnet et al. [2021] Romain Hugonnet, Robert McNabb, Etienne Berthier, Brian Menounos, Christopher Nuth, Luc Girod, Daniel Farinotti, Matthias Huss, Ines Dussaillant, Fanny Brun, et al.  Accelerated global glacier mass loss in the early twenty-first century.  _Nature_ , 592(7856):726–731, 2021.