Sajid Ghuffar
email: sghuffar@gmail.com

Tobias Bolch

Ewelina Rupnik

Atanu Bhattacharya

A Pipeline for Automated Processing of Declassified Corona KH-4 (1962-1972) Stereo Imagery

Keywords: Panoramic Cameras, Spy Satellites, Corona KH-4, SuperGlue, DEM, Glacier Changes, Epipolar Resampling I

The Corona KH-4 reconnaissance satellite missions acquired panoramic stereo imagery with high spatial resolution of 1.8-7.5 m from 1962-1972. 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 the Corona Stereo Pipeline (CoSP): A pipeline for processing of Corona KH-4 stereo panoramic imagery. CoSP utilizes deep learning based feature matcher SuperGlue to automatically match features point between Corona KH-4 images and recent satellite imagery 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. Using the entire frame of the Corona image, bundle adjustment with well-distributed GCPs results in an average standard deviation or σ0 of less than two pixels. We evaluate fiducial marks on the Corona films and show that preprocessing the Corona images to compensate for film bending improves the 3D reconstruction accuracy. The distortion pattern of image residuals of GCPs and y-parallax in epipolar resampled images suggest that film distortions due to long-term storage likely cause systematic deviations of up to six pixels. Compared to the SRTM DEM, the Corona DEM computed using CoSP achieved a Normalized Median Absolute Deviation of elevation differences of ≈ 4 m over an area of approx. 4000 km 2 after a tilebased fine coregistration of the DEMs. We further assess CoSP on complex scenes involving high relief and glacierized terrain and show that the resulting DEMs can be used to compute long-term glacier elevation changes over large areas.

highly valuable in numerous applications [START_REF] Galiatsatos | Assessment of the Corona series of satellite imagery for landscape archaeology: A case study from the Orontes valley, Syria[END_REF][START_REF] Altmaier | Digital surface model generation from Corona satellite images[END_REF][START_REF] Schmidt | The use of Corona satellite images for generating a high resolution digital elevation model[END_REF]. 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 initiated in 1958 and there were 102 Corona missions between 1960 and 1972. These missions acquired more than 800,000 images, which were declassified in 1995 through a presidential order (Executive Order 12951) [START_REF] Usgs | Declassified intelligence satellite photographs, fact sheet 2008-3054[END_REF]. 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 [START_REF] Mcdonald | Corona between the Sun & the Earth: the first NRO reconnaissance eye in space[END_REF][START_REF] Dashora | A repository of Earth resource information -Corona satellite programme[END_REF]. In the KH-1 to KH-3 missions, the satellite imaging payload consisted of a single panoramic camera with a ground resolution between 7.5-12 m. The KH-4 series, which began in 1962, consisted of a dual panoramic camera system with fore (forward) and aft (backward) looking cameras having a 15 • 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 • 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 [START_REF] Nro | Corona technical information, volume II[END_REF].

The main difference in the KH-4, KH-4A and KH-4B cameras was the utilization of a rotating lens/scan arm in the KH-4B camera system instead of an oscillating scan arm in the KH-4/4A camera systems [START_REF] Madden | The Corona camera system: Itek's contribution to world stability[END_REF]. This provided greater 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 m. 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.

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 Corona missions have undergone continuous design improvement since their inception. The cameras in the Corona missions were designed by Itek Corporation and details are given in declassified documents available from National Reconnaissance Office (NRO) [START_REF] Nro | Index to the declassified collection of CORONA, ARGON and LANYARD records declassified by the National Reconnaissance Office on november 26[END_REF][START_REF] Nro | Description and operational manual: The J-3 panoramic camera system[END_REF][START_REF] Nro | Mural system performance evaluation report[END_REF]. Table I summarizes the main specifications of the KH-4 series missions.

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:

• We devise a scheme to automatically generate GCPs between Corona scenes and modern satellite imagery of different spatial resolutions, such as PlanetScope and Landsat-7 ETM+ images. • We experimentally show that the rigorous camera model with the time-dependent collinearity equations can achieve a Standard Deviation (SD) of better than two 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 Image Motion Compensation (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, Panoramic Geometry (PG) stripes). We show that the 3D accuracy can be improved by film bending correction. • We show that the entire Corona stereo pair can be rectified with a residual y-parallax of less than one pixel but suspected film distortions lead to larger systematic yparallax of up to six pixels in certain parts of the image; • We further show a large-scale application of this pipeline by applying it to 24 Corona images, over a footprint of 160x200 km.

II. 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 [START_REF] Bayram | Coastline change detection using Corona, SPOT and IRS 1D images[END_REF], lakes [START_REF] Saruulzaya | Thermokarst lake changes in the southern fringe of Siberian permafrost region in Mongolia using Corona, Landsat, and ALOS satellite imagery from 1962 to 2007[END_REF] and urban areas [START_REF] Noaje | Environmental changes analysis in Bucharest City using Corona, SPOT HRV and IKONOS images[END_REF]. Most of the earlier work on Corona imagery, however, has focused on applications in archaeology and glaciology. Corona imagery has been used for identification and interpretation of archaeological features such as crop marks and built-up structures [START_REF] Galiatsatos | Assessment of the Corona series of satellite imagery for landscape archaeology: A case study from the Orontes valley, Syria[END_REF][START_REF] Fowler | Detection of archaeological crop marks on declassified Corona KH-4B intelligence satellite photography of southern england[END_REF][START_REF] Casana | Stereo analysis, DEM extraction and orthorectification of Corona satellite imagery: archaeological applications from the Near East[END_REF][START_REF] Watanabe | Uti-lization of structure from motion for processing Corona satellite images: Application to mapping and interpretation of archaeological features in Liangzhu Culture, China[END_REF][START_REF] Casana | The Corona atlas project: Orthorectification of Corona satellite imagery and regionalscale archaeological exploration in the Near East[END_REF], computation of long-term glacier area changes [START_REF] Bhambri | Glacier changes in the Garhwal Himalaya, India, from 1968 to 2006 based on remote sensing[END_REF][START_REF] Schmidt | Changes of high altitude glaciers from 1969 to 2010 in the Trans-Himalayan Kang Yatze Massif, Ladakh, northwest India[END_REF][START_REF] Robson | Spatial variability in patterns of glacier change across the Manaslu range, Central Himalaya[END_REF][START_REF] Racoviteanu | Spatial patterns in glacier characteristics and area changes from 1962 to 2006 in the Kanchenjunga-Sikkim area, Eastern Himalaya[END_REF], estimation of glacier volume change and mass balance [START_REF] Bolch | Planimetric and volumetric glacier changes in the Khumbu Himal, Nepal, since 1962 using Corona, Landsat TM and ASTER data[END_REF][START_REF] Robson | Spatial variability in patterns of glacier change across the Manaslu range, Central Himalaya[END_REF][START_REF] Goerlich | Glacier mass loss during the 1960s and 1970s in the Ak-Shirak range (Kyrgyzstan) from multiple stereoscopic Corona and Hexagon imagery[END_REF][START_REF] King | Six decades of glacier mass changes around Mt. Everest are revealed by historical and contemporary images[END_REF], assessing area elevation and surface velocity changes of rock glaciers [START_REF] Bolch | Occurrence, evolution and ice content of ice-debris complexes in the Ak-Shiirak, Central Tien Shan revealed by geophysical and remotely-sensed investigations[END_REF][START_REF] Kääb | Inventory and changes of rock glacier creep speeds in Ile Alatau and Kungöy Ala-Too, northern Tien Shan, since the 1950s[END_REF] and identification/mapping of glacial lakes [START_REF] Bolch | Identification of potentially dangerous glacial lakes in the northern Tien Shan[END_REF][START_REF] Raj | Remote sensing-based hazard assessment of glacial lakes in Sikkim Himalaya[END_REF]. Bhattacharya et al. [START_REF] Bhattacharya | High Mountain Asian glacier response to climate revealed by multi-temporal satellite observations since the 1960s[END_REF] 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 [START_REF] Pieczonka | Region-wide glacier mass budgets and area changes for the Central Tien Shan between˜1975 and 1999 using Hexagon KH-9 imagery[END_REF][START_REF] Dehecq | Automated processing of declassified KH-9 Hexagon satellite images for global elevation change analysis since the 1970s[END_REF][START_REF] Maurer | Acceleration of ice loss across the Himalayas over the past 40 years[END_REF][START_REF] King | Glacial lakes exacerbate Himalayan glacier mass loss[END_REF]. 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-9 m) as well as image acquisitions earlier in time. This can extend 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 [START_REF] Kääb | Inventory and changes of rock glacier creep speeds in Ile Alatau and Kungöy Ala-Too, northern Tien Shan, since the 1950s[END_REF][START_REF] Bhattacharya | High Mountain Asian glacier response to climate revealed by multi-temporal satellite observations since the 1960s[END_REF].

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. [START_REF] Sohn | Mathematical modelling of historical reconnaissance Corona KH-4B imagery[END_REF] evaluated three different mathematical models for KH-4B panoramic imagery. Their first model used modified collinearity equations, which included panoramic geometry transformation along with scan positional distortion and 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 timedependent 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 model used the Rational Polynomial Coefficients (RPC) [START_REF] Fraser | Sensor orientation via RPCs[END_REF][START_REF] Tao | A comprehensive study of the rational function model for photogrammetric processing[END_REF] based generic sensor model to represent 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 m using GCPs covering an area of around 560 km 2 (approx. 20% area of entire stereo overlap), which was better than the height accuracy achieved using the RPC based model. Shin and Schenk [START_REF] Shin | Rigorous modeling of the first generation of the reconnaissance satellite imagery[END_REF] 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 m on the check points spread over entire Corona stereo pair covering Ohio, US. The NASA Ames Stereo Pipeline (ASP) [START_REF] Beyer | The ames stereo pipeline: Nasa's open source software for deriving and processing terrain data[END_REF] has implemented KH-4B panoramic camera model based on the model of Shin and Schenk [START_REF] Shin | Rigorous modeling of the first generation of the reconnaissance satellite imagery[END_REF] and motion compensation from Sohn et al. [START_REF] Sohn | Mathematical modelling of historical reconnaissance Corona KH-4B imagery[END_REF]. Lauer [START_REF] Lauer | Exploiting space-based optical and radar imagery to measure and model tectonic deformation in continental areas[END_REF] has used a fisheye camera model for Corona imagery and obtained an average height accuracy of 12.5 m for the control points. Jacobsen [START_REF] Jacobsen | Calibration and validation of Corona KH-4B to generate height models and orthoimages[END_REF] used the panoramic transformation terms in the perspective frame camera model and observed an average standard deviation of 11.4 pixels using GCPs over multiple consecutive Corona stereo pairs. Multiple earlier studies [START_REF] Bolch | Planimetric and volumetric glacier changes in the Khumbu Himal, Nepal, since 1962 using Corona, Landsat TM and ASTER data[END_REF][START_REF] Goerlich | Glacier mass loss during the 1960s and 1970s in the Ak-Shirak range (Kyrgyzstan) from multiple stereoscopic Corona and Hexagon imagery[END_REF][START_REF] Bhattacharya | High Mountain Asian glacier response to climate revealed by multi-temporal satellite observations since the 1960s[END_REF] have used the Remote Sensing Software Package Graz (RSG) developed by Joanneum Research Graz, for processing of Corona imagery to create DEMs for glacier volume change estimation. The RSG Corona camera model is also based on panoramic coordinates transformation in the collinearity equations similar to the first model of Sohn et al. [START_REF] Sohn | Mathematical modelling of historical reconnaissance Corona KH-4B imagery[END_REF] 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.

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. [START_REF] Sohn | Mathematical modelling of historical reconnaissance Corona KH-4B imagery[END_REF] 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 whole image. There are also differences in the model parameterization in different studies. For example, in RSG as well as in the work of Jacobsen [START_REF] Jacobsen | Calibration and validation of Corona KH-4B to generate height models and orthoimages[END_REF], the Corona 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 [START_REF] Brown | V/H image motion in aerial cameras (image motion in focal plane of aerial camera result-ing from forward motion of vehicle)[END_REF]. 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. Jacobsen [START_REF] Jacobsen | Calibration and validation of Corona KH-4B to generate height models and orthoimages[END_REF] has shown that the bending of the film has a significant magnitude and it should be corrected for achieving higher accuracy.

III. METHODOLOGY

The fore and aft looking camera of Corona KH-4 missions performs a 70 • scan in the across track direction with a rotat- . ing 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 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. [START_REF] Sohn | Mathematical modelling of historical reconnaissance Corona KH-4B imagery[END_REF]. 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. Our developed pipeline, CoSP consists of several different components (Figure 2):

A. Recovery of the film geometry a) Film scans stitching: The physical length of the film (across track direction) containing one panoramic scan is around 745 mm, while the film width (flight direction) is 70 mm. Due to its large size, U.S. Geological Survey (USGS) scans the film in four parts (labeled: a, b, c, d) with an overlap between the successive scans. USGS conducts film scanning at 7 µm or 14 µm resolution. In this work we have used scans at 7 µ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 areas of the scans and estimating a transformation (2D rotation and translation) to align them. We first stitch a with b and c with d, then the parts ab and cd are merged together.

b) 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, holes of 75 microns were carved in the film restraining rail. A light source on the scan arm shone through these holes, which produced a sequence of spots with approx. 1 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, the Panoramic Geometry (PG) stripes/lines on either side of the film format were produced by lamps mounted on the lens, which aid in estimating the axis of the lens rotation. The films of earlier Corona missions i.e. KH-4 and KH-4A did not contain rail holes or PG stripes and instead contained shrinkage markers and center of format indicator, which were created from a v-shaped cutout on the film guide track [START_REF] Nro | Mural system performance evaluation report[END_REF]. In addition to the fiducial marks the KH-4 films also contained titling data and 200 cycles per second time marks (Figure 3).

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 reconstructing the panoramic geometry and correcting the film bending and distortions. A certain number of positions of the PG stripes and the rail holes were manually extracted along the length of the film for evaluation purpose. These positions were then interpolated for each pixel along the length of the film using spline interpolation. The positions of the rail holes along the two sides of the film unexpectedly deviate from each other and from the respective PG stripes (see Figure 4). This may be due to the fact that the rail holes are closer to the bottom/top edges 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 format. 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.

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). For the KH-4B images, rail holes and the center of format markers are extracted using normalized cross correlation with a template of the known size of the markers (Figure 3). The outliers are filtered by the Random Sample Consensus (RANSAC) algorithm using known positions of these markers. The extracted positions of these rail holes lead to the estimation of the rotation to align the scanned image with the axes of the image format. For film bending correction, PG stripes are extracted using a Laplacian of Gaussian kernel (size 9 × 1) averaged along the length of the PG stripe. The extrema of the filtered image indicates the center of the PG stripes. Using RANSAC an 8 th degree polynomial is fitted to these extrema locations to estimate the position of the PG stripes along the length of the film. A correction grid along the y-axis based on the difference of the PG stripes position with respect to their y-coordinates at the center of the image format is computed. The images are then resampled using this grid, so that the PG stripes become straight and parallel to each other. Finally, we clip the 1.5 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 [START_REF] Nro | The non-specified area of J-3 format[END_REF]. The aft camera images are rotated 180 • to align them with fore camera images. For KH-4 and KH-4A missions, the four v-shaped markers are extracted using normalized cross correlation with the v-shape template. The positions of these markers are used to clip and align the exposed image area of the film.

B. GCP Generation

In the earlier studies, the extraction of GCPs was done manually and has been one of the most time consuming component of the Corona imagery processing pipeline [START_REF] Jacobsen | Calibration and validation of Corona KH-4B to generate height models and orthoimages[END_REF][START_REF] Bhattacharya | High Mountain Asian glacier response to climate revealed by multi-temporal satellite observations since the 1960s[END_REF]. Automatic feature matching using SIFT [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF] or other traditional feature detectors is challenging because of physical changes in the scene since image acquisition. The radiometric differences and variations in cast shadows in high-relief terrain makes feature matching particularly difficult. Recent progress 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 [START_REF] Maiwald | An automatic workflow for orientation of historical images with large radiometric and geometric differences[END_REF][START_REF] Zhang | Feature matching for multi-epoch historical aerial images[END_REF].

We use SuperGlue [START_REF] Sarlin | Superglue: Learning feature matching with graph neural networks[END_REF] a deep learning based architecture for feature extraction and matching to compute point matches between the Corona and reference imagery. The SuperGlue architecture consists of a Convolutional Neural Network (CNN) based interest point detector and descriptor called SuperPoint [START_REF] Detone | Superpoint: Self-supervised interest point detection and description[END_REF], which uses VGG (Visual Geometry Group) [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF] style encoder. It uses a graph neural network to perform matching of the features extracted by the CNN and implicitly learn 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. SuperGlue, like any other CNN based architecture, is typically trained on relatively low resolution images. The default image resolution for SuperGlue is 640 × 480 pixels and we found in our experiments that it provides reasonable results up to three times the default resolution i.e, 1920 × 1440. The Corona image has a much larger size i.e. ∼ 106000 × 8000. So, we divide it in to patches of size 1920 × 1440 and extract similar geographic location patch from reference satellite imagery to match the feature points. As the resolution of KH-4B images is around 2-3 m, we use a mosaic of orthorectified PlanetScope images [START_REF] Team | Planet application program interface: In space for life on earth[END_REF], which has a similar spatial resolution (3 m) as reference satellite imagery to match feature points. As reference satellite imagery with a high resolution of 1-5m may not be available for larger areas, we also experimented with feature matching on freely available Landsat-7 ETM+ panchromatic band (15 m). For this, we resample the Corona imagery to have similar spatial resolution of 15 m. Specifically, we divide Corona image into tiles of ∼ 10600 × 8000 and resize it to 1920 × 1440 for input to SuperGlue. We used the original weights of the SuperGlue model as provided by the authors [START_REF] Sarlin | Superglue: Learning feature matching with graph neural networks[END_REF]. Our experimentation revealed that SuperGlue was successful in matching points between Corona KH-4 and satellite imagery like the Landsat series or very high resolution satellite imagery. Hence, other satellite imagery may also be used for the purpose of feature matching and GCPs generation.

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 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. We adopt a coarse-tofine strategy, whereby 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.

C. Camera Model

During image acquisition, the scene is exposed to the film through a slit, which rotates along with the lens. The scan angle α is the angle of the slit and is given by

α = 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 timedependent parameters

X 0t = X 0 + X 01 • t, (2)
Y 0t = Y 0 + Y 01 • t, (3)
Z 0t = Z 0 + Z 01 • t, (4)
ω 0t = ω 0 + ω 01 • t, (5)
ϕ 0t = ϕ 0 + ϕ 01 • t, (6)
κ 0t = κ 0 + κ 01 • t. (7)
Here, X 01 , Y 01 , Z 01 represent the translation of the camera and ω 01 , ϕ 01 , κ 01 represent the angular motion along the three axis during the panoramic scan. It is assumed that the translational and rotational velocity during the scan time is constant. The duration of the scan was typically around 0.36 s and here, we assume a normalized time i.e. t ∈ [0, 1] . Considering the width of the slit to be small, the geometry of the image formed at any position of the slit is given by,

  0 y p + y IM C -f   = sR α R   X -X 0t Y -Y 0t Z -Z 0t   , (8)
where R α is,

R α =   cosα 0 sinα 0 1 0 -sinα 0 cosα   , (9)
Here, y p is the panoramic photo coordinate perpendicular to x p and X,Y,Z are the coordinates of the object point and the IMC term is given as,

y IM C = - V f Hδ sin(α)cos(ω 0 ), (10)
where, V is the satellite velocity, H is altitude of the satellite and δ is the scan angular velocity i.e. the change of scan angle α with time, which was around 3.3 rad s -1 . The image motion caused by a moving platform is directly related to V/H ratio [START_REF] Kawachi | Image displacements on a tilted panoramic photograph[END_REF][START_REF] Itek-Laboratories | Panoramic progress-part II. Technical analysis of panoramic photography[END_REF][START_REF] Brown | V/H image motion in aerial cameras (image motion in focal plane of aerial camera result-ing from forward motion of vehicle)[END_REF]. The velocity V of the satellite in low Earth orbit is around 7.7 km s -1 and the orbital elevation or height H of Corona satellites typically varied between 150 km to 400 km. 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 [START_REF] Nro | Description and operational manual: The J-3 panoramic camera system[END_REF]. The image motion compensation produces an S-shape distortion across the image length [START_REF] Itek-Laboratories | Panoramic progress-part II. Technical analysis of panoramic photography[END_REF]. By rearranging the terms of Equation ( 8) we obtain,

R T α   0 y p + y IM C -f   = sR   X -X 0t Y -Y 0t Z -Z 0t   , (11)
and finally,  

f sinα y p + y IM C -f cosα   = sR   X -X 0t Y -Y 0t Z -Z 0t   . (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α = - N x N z , (16)
y p + y IM C = -f cosα N y N z , (17)
The final image panoramic coordinates (x p , y p ) being,

x p = f tan -1 (- N x N z ) (18)
y p = -y IM C -f cosα N y N z (19)
There are 13 unknown camera parameters in the model: six parameters for camera position X 0 , Y 0 , Z 0 , X 01 , Y 01 , Z 01 , six parameters for camera orientation ω 0 , ϕ 0 , κ 0 , ω 01 , ϕ 01 , κ 01 and one parameter for the IMC i.e. V Hδ . We keep the focal length fixed in the adjustment due to the correlation with the camera altitude parameter [START_REF] Zhou | Simulation and analysis of photogrammetric UAV image blocks-influence of camera calibration error[END_REF]. The approximate values for ω 0 were set to {-15 • , 15 • } for the fore and aft looking cameras respectively, κ 0 as 0 • , camera altitude Z 0 is set to 170 km and the time-dependent parameters are initialized as zero. These parameters are then optimized in the bundle adjustment using GCPs as well as tie points between the stereo pairs. The bundle adjustment is performed in an Earth Centered Cartesian frame.

D. Stereo Rectification

We follow the stereo rectification algorithm given in [START_REF] Deseilligny | Epipolar rectification of a generic camera[END_REF] and implemented in MicMac [START_REF] Pierrot-Deseilligny | MicMac, apero, pastis and other beverages in a nutshell[END_REF]. This epipolar rectification is generic in the sense that it makes no assumption about 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 stereopair (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 [START_REF] Lowe | Distinctive image features from scaleinvariant keypoints[END_REF] or SuperGlue [START_REF] Sarlin | Superglue: Learning feature matching with graph neural networks[END_REF]. Note that to obtain unambiguous rectifications using this variant, the 3D scene should not be flat [START_REF] Deseilligny | Epipolar rectification of a generic camera[END_REF]. Here, we use a 4th degree polynomial to approximate the epipolar curves.

E. 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) [START_REF] Hirschmuller | Stereo processing by semiglobal matching and mutual information[END_REF] in MicMac, which is based on weighted multiscale image matching. 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. The unreliable pixels in the dense matching were estimated using left-right consistency check and filtered out.

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:

R 1 + tan(α)R 3 R 2 + yp+y IM C f cos(α) R 3   X Y Z   = R 1 + tan(α)R 3 R 2 + yp+y IM C f cos(α) R 3 )   X 0t Y 0t Z 0t   .
(20) Here, R 1 , R 2 and R 3 are the 1 st , 2 nd and 3 rd rows of the rotation matrix respectively. Given the four equations (two equations for each image), 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. To compensate for these systematic deviations, we perform a fine coregistration of Corona DEM and the reference DEM in a post-processing step. We divide the Corona DEM into smaller tiles and estimate a 3D affine transformation using Least Square Matching (LSM) [START_REF] Ressl | Applying 3-D affine transformation and least squares matching for airborne laser scanning strips adjustment without GNSS/IMU trajectory data[END_REF]. We use a tile size of approx. 400 km 2 and estimate the parameters of the transformation individually for each tile. To ensure smooth transition, we keep a two third overlap between the between neighbouring tiles. The LSM approach uses a robust weighting function to filter the outliers in the least squares adjustment process.

We have implemented the camera model and bundle adjustment in MATLAB. Implementation of dense matching, tie

IV. DATA AND STUDY AREAS

We applied and evaluated CoSP using three different case studies (see Table II). We first focus on the evaluation of the camera model and the corresponding Corona DEM. We searched a study site for which well-distributed GCPs can be automatically extracted. We selected an area bordering Mongolia and China (Figure 7) as this site has only shallower slopes, distinctive texture and had minimal natural or anthropogenic changes over the past decades. We chose two Corona stereo pairs from mission no. 1112 (23-11-1970) and 1117 (30-05-1972) covering this area. The images were downloaded from the USGS Earth Explorer. We use SRTM 1 Arc-Second Global DEM v3 (DOI: 10.5066/F7PR7TFT) to extract the elevation values, which ranges from 300-2000 m in the area. 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 m resolution reference DEM (for elevation coordinates).

Our second focus is the evaluation of the proposed pipeline for estimating glacier volume change. To do so, we chose an area around Mount Everest in the Himalaya because of the availability of a high resolution contemporary Cartosat DEM (10 m) as well as Corona DEM from earlier studies computed using RSG [START_REF] Bolch | Multi-decadal mass loss of glaciers in the Everest area (Nepal Himalaya) derived from stereo imagery[END_REF][START_REF] King | Six decades of glacier mass changes around Mt. Everest are revealed by historical and contemporary images[END_REF]. 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 their hydrological importance and accelerating mass loss [START_REF] King | Six decades of glacier mass changes around Mt. Everest are revealed by historical and contemporary images[END_REF]. We selected a Corona stereo pair from mission no. 1108 dated 18-12-1969 (Table II) to create Corona DEM. Here, we use 30 m ALOS World 3D DEM (AW3D30) [START_REF] Tadono | Precise global DEM generation by ALOS PRISM[END_REF] as a reference DEM because of the bias in SRTM DEM due to the penetration of radar signal into snow and ice [START_REF] Gardelle | Impact of resolution and radar penetration on glacier elevation changes computed from DEM differencing[END_REF]. 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 II), covering a footprint of approx. 160 × 200 km. This sequence of images are suitable for large scale glacier change estimation due to late summer acquisition and cloud free coverage.

V. RESULTS

A. Accuracy of the Camera Model (Mongolia Data) a) Accuracy of Check Points: SuperGlue extracted a large number of features points between the Corona images and the high resolution satellite imagery (see Figure 8). Fewer points were detected in areas without texture as well as in 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 SD (σ 0 ) of bundle adjustment is ∼ 2 pixels for both image pairs (Table III). The σ 0 shows an improvement of 10-20% when bending correction is applied. The RMSE of the checkpoints for the two image pairs are [4.85, 3.81, 5.79] and [6. 63, 4.23, 8.42] (m) (Table III). As with the σ 0 , the RMSE decreases when film bending correction is performed.

b) 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. c) Camera parameters: Overall, the estimated camera parameters correspond relatively well to the expected values (see Table IV). 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 in a northerly direction i.e. Y 01 , whose estimated value ranges from 2.59-3.10 km. The expected value of motion along flight direction for a scanning time of 0.36 s is 2.8 km (velocity = 7.7 km 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 relatively high differences between some of the estimated time-dependent orientation parameters of the fore and aft cameras can perhaps be attributed to the uncompensated systematic errors.

The specifications of the J-3 camera state 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 [START_REF] Nro | Description and operational manual: The J-3 panoramic camera system[END_REF]. This may explain why: 1) The estimated value of the IMC term V Hδ is significantly smaller than the expected value of 0.014 (V= 7.7 km s -1 , H = 170 km, δ= 3.3 rad s -1 ) as both image pairs are from KH-4B missions and 2) A value between 0.76-1.09 deg for the parameter ω 01 because the rotation to compensate camera motion will be observed in ω 01 . Given a scan time of 0.36 sec, the 7.7 km 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 ω 01 is consistent with the rotation required to compensate the image motion. To further assess the effect of the IMC mechanism on the estimated parameters, we estimated the camera parameters for different KH-4A scenes using GCPs generated from the scheme introduced in Section III-B. In the majority of the images, the estimated IMC term has a value closer to the value of 0.014 (Table V).

B. Epipolar Resampling and the Accuracy of the DEM (Mongolia Data)

The y-parallax in the epipolar resampled images is computed using 2D correlation algorithm implemented in MicMac [START_REF] Rosu | Measurement of ground displacement from optical satellite image correlation using the free open-source software MicMac[END_REF]. The SD of the y-parallax for the two rectified stereo pairs is 0.89 and 0.99 pixels. The y-parallax over the majority of the image area is less than ±1 pixels (Figure 10). However, certain regions with higher y-parallax up to six 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 (cf. Figure 9).

The Normalized Median Absolute Deviation (NMAD) of elevation differences between the two Corona DEMs and the SRTM DEM (30 m) are 7.26 m and 6.69 m (Figure 11). The elevation differences show systematic deviations of up to 25 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 (cf. Section III-E). This alignment of the DEMs reduces the NMAD of the elevation differences to 3.32 m and 4.15 m for the two Corona DEMs (Figure 11). The size of each tile is approximately 400 km 2 , which is large enough to avoid fitting the 3D affine transform to the temporal surfaces changes that we are interested in estimating.

C. 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. The difference in the Sun position and the cast shadows further complicates the automatic feature matching across mountainous terrain. SuperGlue matched several points between Corona imagery and Landsat-7 ETM+ mosaic (Figure 12), but 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. As no GCPs were available in this section of the image, the resulting DEM shows large elevation differences in comparison with the AW3D30 (Figure 13). The coregistration of the Corona DEM with the AW3D30 DEM, using off-glacier topography removed the majority of the systematic elevation differences between the Corona and the AW3D30. To filter out the glacierized part of the scene, we used the glacier polygons from the work of King et al. [START_REF] King | Six decades of glacier mass changes around Mt. Everest are revealed by historical and contemporary images[END_REF]. Elevation differences also show errors caused by textureless snow/ice cover in the accumulation zone of the glaciers and areas of steep rockwalls. Further large deviations are due to the presence of clouds in the Corona imagery (Figure 13).

The highest spatial resolution in the KH-4 missions is around 1.8 m at the center of the image and around 3 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 m Corona DEM of the Everest region and compare it with a 10 m Cartosat DEM of 15/12/2018 (Figure 14) from the work of King et al. [START_REF] King | Six decades of glacier mass changes around Mt. Everest are revealed by historical and contemporary images[END_REF]. The elevation decrease over the glacierized terrain are well visible and highlight glacier volume loss and thus a negative glacier mass balance. The elevation differences over the stable terrain indicate 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 erroneous elevation differences (Figure 14). Here we generated a spatially complete Corona DEM by interpolating the data gaps due to filtering of unreliable disparity values in dense image matching.

The NMAD of the elevation differences between CoSP Corona DEM and Cartosat DEM over the stable terrain is 10.18 m, which is slightly better than the NMAD of 12.41 m obtained by Corona DEM generated using RSG [START_REF] King | Six decades of glacier mass changes around Mt. Everest are revealed by historical and contemporary images[END_REF] (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 clearly different pattern, which is perhaps due to differences in dense matching algorithms used in CoSP and RSG. After performing the outlier removal using the method of Pieczonka and Bolch [START_REF] Pieczonka | Region-wide glacier mass budgets and area changes for the Central Tien Shan between˜1975 and 1999 using Hexagon KH-9 imagery[END_REF] and gap filling using the median elevation differences [START_REF] Mcnabb | Sensitivity of glacier volume change estimation to DEM void interpolation[END_REF], the elevation difference from 1969 to 2018 over the glacierized area (shown in Figure 14) is -13.55±3.91 m for CoSP and -10.26±4.13 m for RSG computed Corona DEMs. The uncertainty in the elevation differences is computed using the method given in Gardelle et al. [START_REF] Gardelle | Region-wide glacier mass balances over the Pamir-Karakoram-Himalaya during 1999-2011[END_REF]. The elevation change for the Khumbu Glacier is -16.96 ± 3.91 m for CoSP and -17.20 ± 4.13 m for RSG generated Corona DEMs (using the area-wide uncertainty).

D. Large scale Mapping using CoSP (Pamir Data)

In order to show the applicability for large scale mapping, we apply the proposed methodology to 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 m AW3D30. The DEM differences show that CoSP works on a larger scale and there are few systematic elevation errors 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 [START_REF] Kotlyakov | Monitoring surging glaciers of the pamirs, central asia, from space[END_REF][START_REF] Guillet | A regionally resolved inventory of high mountain asia surge-type glaciers, derived from a multi-factor remote sensing approach[END_REF], can be well identified in the DEM differences.

However, this sequence of Corona imagery also shows a potential complication in large scale processing. The consecutive images of the sequence have an overlap of approximately 6% at the center of the format, which is less than the specified overlap of approximately 8%. While the fore and aft looking cameras have an overlap of approximately 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 × 200 km) can be challenging in presence of unstable terrain.

VI. DISCUSSION

A. Camera Model

The results show that a SD (σ 0 ) of better than two 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 [START_REF] Surazakov | Positional accuracy evaluation of declassified hexagon KH-9 mapping camera imagery[END_REF][START_REF] Maurer | Tapping into the Hexagon spy imagery database: A new automated pipeline for geomorphic change detection[END_REF]. Furthermore, bowing/buckling of the film as well as any deviation from the circular position of the film, 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 [START_REF] Sohn | Mathematical modelling of historical reconnaissance Corona KH-4B imagery[END_REF]. 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 m and 8.42 m for the two Corona stereo pairs (Table III), which is higher than the height RMSE of check points of around 4 m reported in Sohn et al. [START_REF] Sohn | Mathematical modelling of historical reconnaissance Corona KH-4B imagery[END_REF]. This decrease in the accuracy is expected due to the presence of film distortions, when considering the entire Corona stereo coverage (area >3000 km 2 ) in comparison to the subset (area of 561 km 2 ) used by [START_REF] Sohn | Mathematical modelling of historical reconnaissance Corona KH-4B imagery[END_REF]. As Figures 9 and 10 suggest, residuals distortions may have a higher magnitude in certain parts of the image and considering a larger image area may decrease the accuracy of check points due to these distortions. The tile-based fine coregistration with the reference DEM improves the accuracy and reduces the SD (σ) of elevation differences to around 4 m to 5 m (Figure 11). This shows that accuracy similar to [START_REF] Sohn | Mathematical modelling of historical reconnaissance Corona KH-4B imagery[END_REF] 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.

[40] have reported a height RMSE of 12.34 m for checkpoints, which is higher than RMSE obtained in our work. However, they used a scan with 12 µm pixel size, which reduces accuracy in comparison to 7 µm pixel size used in this work. The camera model in [START_REF] Shin | Rigorous modeling of the first generation of the reconnaissance satellite imagery[END_REF] was based on time-dependent collinearity equations but they only modeled motion along flight direction as an additional time-dependent parameter.

Our results further show an improvement in accuracy in comparison to the work of [START_REF] Bhattacharya | High Mountain Asian glacier response to climate revealed by multi-temporal satellite observations since the 1960s[END_REF], 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 m Cartosat DEM also exhibit slight improvement in comparison to the Corona DEM generated using RSG in the work of [START_REF] King | Six decades of glacier mass changes around Mt. Everest are revealed by historical and contemporary images[END_REF] (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.

B. Corona Data Limitations

1) Film Reference Data: The reference metadata exposed on the film is important for the reconstruction of the panoramic geometry and correction of film distortions. However, as shown earlier, the PG rail holes show unexpected deviations and are not reliable for the 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, PG stripes and rail holes are not fully visible in some images. Hence, image clipping and bending correction as performed here is not applicable to all Corona KH-4 images, which will have an effect on the expected accuracy.

2) Low Texture and Film Saturation: The low texture and contrast over the snow covered part of the glacier surfaces is a notable limitation in the utilization of the Corona imagery for glacier mass balance estimation. This is caused by the saturation due to light exposed on the film as well as limited surface roughness at the spatial resolution of the image. 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 [START_REF] Goerlich | Glacier mass loss during the 1960s and 1970s in the Ak-Shirak range (Kyrgyzstan) from multiple stereoscopic Corona and Hexagon imagery[END_REF]. 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 these areas.

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 [START_REF] Nro | Description and operational manual: The J-3 panoramic camera system[END_REF]. 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 and may vary with orbital elevation. The overlap between the successive images is rather small, which results in a small overlap between the DEMs from the successive stereo pairs [START_REF] Goerlich | Glacier mass loss during the 1960s and 1970s in the Ak-Shirak range (Kyrgyzstan) from multiple stereoscopic Corona and Hexagon imagery[END_REF][START_REF] Jacobsen | Calibration and validation of Corona KH-4B to generate height models and orthoimages[END_REF]. 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%.

4) Scanning Artifacts: The scanned imagery delivered by USGS has been reported to contain scanning artifacts [START_REF] Lauer | Exploiting space-based optical and radar imagery to measure and model tectonic deformation in continental areas[END_REF][START_REF] Galiatsatos | High resolution elevation data derived from stereoscopic Corona imagery with minimal ground control[END_REF][START_REF] Gheyle | Scan problems in digital Corona satellite images from USGS archives[END_REF][START_REF] Dehecq | Automated processing of declassified KH-9 Hexagon satellite images for global elevation change analysis since the 1970s[END_REF]. These artifacts occur mainly due to an error in the scanner calibration. After processing KH-4 imagery of different missions, we have also observed such artifacts in the scanned Corona images (Figure 16). 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. As the Corona images have been available since now two decades, 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 and11. The stitching of individual scans using an affine transformation can allow compensation of any scaling or shear component during the scanning process. However, our experimentation indicated that the presence of blocking artifacts limits a reliable estimation of these components. Hence, we estimate only a rotation and translation to stitch the individual scans.

C. GCP Generation

Automatic GCP generation is an essential part of automated processing of Corona imagery. GCP generation using SuperGlue has largely been successful in generating automatic GCPs not only over flat terrain with high contrast but also over glacierised mountainous 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. 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 use DEMs for matching features for georeferencing of the historical imagery [START_REF] Maurer | Tapping into the Hexagon spy imagery database: A new automated pipeline for geomorphic change detection[END_REF][START_REF] Dehecq | Automated processing of declassified KH-9 Hexagon satellite images for global elevation change analysis since the 1970s[END_REF][START_REF] Zhang | Feature matching for multi-epoch historical aerial images[END_REF].

D. DEM Coregistration

We used a tile-based approach to coregister the Corona DEM with a reference DEM. This worked well in our test cases. However, problems in coregistration can appear if large parts of the scene contain unstable terrain, outliers or data voids. Therefore, the tile-based coregistration strategy may not work on scenes with large clouds as well as scenes of ice-sheets. Thus, a coregistration 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.

Other approaches to DEM coregistration, have been widely used for coregistration of DEMs for quantification of glacier elevation change [START_REF] Nuth | Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change[END_REF][START_REF] Pieczonka | Region-wide glacier mass budgets and area changes for the Central Tien Shan between˜1975 and 1999 using Hexagon KH-9 imagery[END_REF][START_REF] Girod | MMASTER: improved ASTER DEMs for elevation change monitoring[END_REF][START_REF] Hugonnet | Accelerated global glacier mass loss in the early twenty-first century[END_REF]. 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 and correction of remaining systematic elevation differences by polynomial fitting. Such approaches may also be investigated for improvement of the DEM coregistration strategy in CoSP.

VII. 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 timedependent exterior orientation parameters was used to model the panoramic image acquisition geometry of the Corona cameras. SuperGlue was applied to automatically match feature points between Corona imagery and medium to high resolution satellite imagery. The accuracy of the camera model was assessed by using these feature points as GCPs. The results show that overall a SD (σ 0 ) better than two pixels can be achieved, while considering the entire Corona image. Systematic image residuals up to six pixels have been observed in different Corona scenes. These errors are probably due to film distortions that have occurred during the mission and the storage in the intervening 50-year period. Consequently, the resulting Corona DEMs exhibit systematic elevation differences of up to 25 m, which could be reduced by a tile-based coregistration of the Corona DEM. The estimated camera parameters are consistent with the known parameters as well as the IMC mechanism on the Corona camera system. A tie point based epipolar resampling of the Corona images showed residual y-parallax with in ±1 pixel 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. Improvements in the coregistration strategy are required to utilize a larger area or whole Corona DEM to make the process more robust to scenes involving clouds, large number of outliers and unstable terrain. 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.

Fig. 1 .

 1 Fig. 1. Geometry of the Corona panoramic stereo imaging system (Modified from NRO [10])

Fig. 2 .

 2 Fig. 2. Workflow of the Corona Stereo Pipeline

Fig. 4 .

 4 Fig. 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.

Fig. 5 .

 5 Fig. 5. The panoramic image coordinate system and the angle α.

Fig. 6 .

 6 Fig. 6. Epipolar curves in the native geometry of a Corona stereo pair (top). These curves are exaggerated for visualization. Note that the curves are globally vertical and they will be rotated by ≈ 90 • prior to estimating the rectifying polynomials. The resulting epipolar rectified fore and aft looking Corona images (bottom). For compact visualization the rectified images have been rotated by 90 • .

Fig. 7 .

 7 Fig. 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

Fig. 8 .

 8 Fig. 8. Point matches between a Corona images and multi satellite Plan-etScope imagery mosaic computed by SuperGlue. The grid type pattern of the feature points is due to tile-based matching.

Fig. 9 .

 9 Fig. 9. Image residuals of the GCPs after resection are interpolated over a regular grid. These images have been corrected for film bending.

Fig. 10 .

 10 Fig. 10. The two image pairs and the corresponding y-parallax in the epipolar resampled image pairs

Fig. 12 .

 12 Fig. 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.

Fig. 13 .

 13 Fig. 13. Corona stereo pair and the elevation differences between the Corona DEM and AW3D30 (30 m). a) and b) The Corona image pair c) Elevation differences with the AW3D30 d) Elevation differences after tilebased coregistration of the Corona DEM and the AW3D30 e) Shaded Corona DEM.

Fig. 14 .

 14 Fig. 14. a) Hillshade of 10 m Corona DEM b) Hillshade of 10 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 [27] e) Histogram of elevation differences given in [c] f) Corona orthophoto of the Everest Region g) Histogram of elevation differences given in [d]

Fig. 15 .

 15 Fig. 15. DEM differences between the Corona DEM computed using 12 Corona stereo pairs and the AW3D30 over Pamir region.

Fig. 16 .

 16 Fig. 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.

TABLE II CORONA

 II SCENES USED FOR EACH STUDY SITE AND THE CORRESPONDING REFERENCE DEMS AND SATELLITE IMAGERY USED FOR GCPS. THESE CORONA SCENES BELONG TO KH-4B SERIES (MISSION IDS 1101-1117).

	Study	Corona Scenes	Date	Imagery	Reference
	Site			for GCPs	DEM
		DS1117-2071DF008 DS1117-2071DA014	30-05-1972
	Mongolia	DS1112-1071DF023 DS1112-1071DA029	23-11-1970	PlanetScope	SRTM
		DS1108-2217DA070			ALOS
	Everest	DS1108-2217DF064	18-12-1969	Landsat-7	Cartosat
		DS1104-2169DF092-
		DS1104-2169DF103
	Pamir	DS1104-2169DA098-	18-08-1968	Landsat-7	ALOS
		DS1104-2169DA109

TABLE III ACCURACY

 III 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	σ 0 (pixels)	RMSE X (m)	RMSE Y (m)	RMSE Z (m)
	Bending Correction
	DS1112-1071DF023 DS1112-1071DA029	1.71	4.85	3.81	5.79
	DS1117-2071DF008 DS1117-2071DA014	1.90	6.63	4.23	8.42
	Without Correction
	DS1112-1071DF023 DS1112-1071DA029	1.89	5.84	4.09	9.92
	DS1117-2071DF008 DS1117-2071DA014	2.43	7.39	4.75	15.11

TABLE IV THE

 IV ESTIMATED CAMERA PARAMETERS ARE GIVEN ACCORDING TO THE GEOGRAPHIC COORDINATE SYSTEM WITH THE POSITION OF THE CAMERA GIVEN AS LONGITUDE, LATITUDE AND ELEVATION.

			DS1112-	DS1112-			DS1117-	DS1117-
	Parameters	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
	ω 0 (deg)			-15.20	15.72			-15.52	14.68
	ϕ 0 (deg)			-1.56	1.46			3.71	-2.25
	κ 0 (deg)			5.69	5.74			-10.19	-10.64
	ω 01 (deg)			0.83	0.94			0.76	1.09
	ϕ 01 (deg)			-0.04	-0.03			-0.51	-0.23
	κ 10 (deg)			0.00	-0.23			0.33	0.02
	V Hδ			0.0025	-0.0001			-0.0002	0.0002
					TABLE V
	ESTIMATED TIME-DEPENDENT ORIENTATION PARAMETERS AND IMC FOR KH-4A SCENES
	Image		X 01	Y 01	Z 01	ω 01		ϕ 01	κ 01	v/Hδ
			(km)	(km)	(km)	(deg)	(deg)	(deg)
	DS1024-1038DF096	-0.187	-3.909	-0.652	-0.020	-0.202	0.050	0.016
	DS1024-1038DF095		0.037	-3.167	-2.307	0.028	-0.06	0.019	0.012
	DS1049-2119DF055		0.463	-3.181	0.453	-0.15	-0.002	0.014	0.017

This article has been accepted for publication in IEEE Transactions on Geoscience and Remote Sensing. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TGRS.2022.3200151 © 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information. Authorized licensed use limited to: ST ANDREWS UNIVERSITY. Downloaded on August 31,2022 at 10:17:32 UTC from IEEE Xplore. Restrictions apply.

This article has been accepted for publication in IEEE Transactions on Geoscience and Remote Sensing. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TGRS.2022.3200151 © 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information.

ACKNOWLEDGMENT

This study was supported by the Strategic Priority Research Program of Chinese Academy of Sciences (XDA20100300) and the Swiss National Science Foundation (200021E 177652/1) within the framework of the DFG Research Unit GlobalCDA (FOR2630). We are thankful to Planet Labs for providing access to PlanetScope imagery through the Education and Research Program. We also thank Marc Pierrot-Deseilligny for his comments and suggestions during this work.