CALIBRATION AND VALIDATION OF CORONA KH-4B TO GENERATE HEIGHT MODELS AND ORTHOIMAGES

At 48 years of age CORONA-KH-4B images are important for actual urban planning in Bangladesh, where no old maps or aerial images exist, indicating locations of former water courses that are no longer visible, but causing problems for the stability of building ground. CORONA KH-4B images are available for very low cost or even free of charge. A forward and a backward-looking panoramic camera with a nadir angle of approximately 15° enabled a stereoscopic coverage. Taken from a height of 154km, an image covers approximately 220km x 14km up to 17km with a height to base relation of 1:1.85. The ground sampling distance (GSD) varies with the incidence angle from 1.77m in the centre of the image up to 2.18m in Y-direction and 2.69m in the X-direction at the ends of the images. These nice conditions of the old images are affected by geometric problems of the panoramic images. A geometric bending can be improved by a correction based on the available sidelines, nevertheless deformations dominantly in the longitudinal direction only can be determined based on ground control points. A group of neighboured images have similar deformations, allowing the determination of a correction grid, describing the systematic image errors. This improved the geometry strongly, but could not eliminate all individual and local geometric problems. The high morphologic quality of the CORONA images, but limited absolute height accuracy has been improved by merging with the highly accurate TDM90 height model. For ortho images a fitting of neighboured images was required.


INTRODUCTION
The required geometric handling of the now declassified panchromatic CORONA images, which can be downloaded from the USGS Earth Explorer, is poorly described in the literature due to the priority of CORONA for reconnaissance and not for the generation of height models and ortho images. (Schneider et al. 2001) and (Dowman et al. 2012) explain the method for using a part of the CORONA image, but no literature is available for the precise handling of the full images. (Auelmann 2015) gives detailed information about the configuration and imaging condition of the CORONA satellite. Also (Dashora, Lohani, Malik, 2007) describes the image configuration, while McDonald explains the CORONA program. The original thin layer film of CORONA images with a usable size of approximately 75cm x 5.5cm is available scanned with 7µm pixel size in 4 overlapping parts. The pixel size leads from a flying elevation of 254km with the focal length of 610mm to the nominal GSD of 1.77m in the centre and 2.18m x 2.69m at the film ends ( figure 1). An analysis of the image quality showed an effective ground resolution being approximately 10% larger as the nominal resolution (Jacobsen 2008), only the first and last approximately 2cm of the original film is clearly worse. The combination of the image parts is not a problem, but the long film format suggests that the image will be deformed. With 2 sidelines the film curvature can be determined and respected, but there are still larger systematic image errors that dominate in the longitudinal direction. With the common additional parameters that are used for aerial and terrestrial cameras, the complex systematic errors of CORONA images cannot be described. In addition the standard image overlap in flight direction of CORONA is limited to 6% in the centre and 20% at the end of images (figure 2), and dynamic imaging with a rotating optics causes systematic errors that may change from image to image. The overlap of the stereo images does not solve geometric problems due to the use of two cameras having different geometric deformations. Therefore, ground control points (GCPs) had to be used for the geometric analysis and improvement. GCPs have been extracted from Google Earth which has at least a relative accuracy of ~ 1.5m. It is difficult to identify GCPs in images with a time difference of 48 years. This causes bigger problems than the accuracy of Google Earth. On average, 72 GCPs were identified in each CORONA image. The even distribution of points was limited by large streams that moved a few kilometres in 48 years and some other areas with major changes, as well as fact that the streets are now covered by trees. The required geometric transformations and the method of determining the systematic errors are described, as well as the improvement of generated height models and the remaining part of the systematic image.

IMAGE ORIENTATION AND GEOMETRIC CORRECTION
The CORONA satellite was equipped with a forward and a backward looking panoramic camera. The panoramic cameras had a classical type of imaging with rotating optics and film located in a curvature (figure 2). The location of the thin film on a circle, and the optics rotating 360° in 1.25 seconds or imaging the whole used image of 75cm within 0.24 seconds, may cause some bending of the film. For the determination of the bending on both side of the film there are sidelines. Only for the first and the last few cm the side lines cannot be used due to being too small up to disappearing. Before sideline correction, the inner orientation was respected based on one fiducial mark at both ends of the film. For few images the fiducial marks have not been scanned, requiring a reconstruction of the fiducial mark positions based on the imaged part of the film. The sidelines could not be used for this due to the circumstance that they are not visible close to the film ends.   Standard photogrammetric solutions are available for perspective cameras, for this reason the panoramic images have been projected to a tangential plane (figure 2). In addition to the panoramic correction shown in figure 5, during the imaging from one side to the other in 0.24 seconds, the satellite is moving 1850m in the orbit what should lead to a small S-shape deformation of the covered area on ground.

Figure 5. Panoramic correctionvectors enlarged
Nevertheless the image motion would cause a smearing, for this reason it was compensated with forward motion compensation by continuous camera pitch change. This compensates the major influence of the satellite motion to the image footprint during imaging. In addition the exact values of image motion are not known, so no image correction for this reason was introduced, and also no systematic error corresponding to an S-shape deformation could be detected, nevertheless it also should be compensated by the later described correction grids. The scanned sub-images of the whole film are fitted together based on tie points in the overlapping parts. Small rotations caused by the film positioning in the scanner, but no scale differences could be seen (figure 6).

Figure 6. configuration of sub-images
The image orientation has been done in a tangential coordinate system of the earth ellipsoid. This orthogonal coordinate system avoids problems with the earth curvature. In case of a handling in the national coordinate system, second order earth curvature effects and scale changes in X and Y would cause not acceptable discrepanciesthe earth curvature over 110km corresponds to approximately 1km and the scale change in the transversal Mercator systems reaches 0.15‰ or 16m over 110km. The orientation by resection with an average of 72 GCPs leads to a sigma0 of 45µm up to 121µm, with in the average 80µm or 11.4 pixels, respectively 23m on ground.

Figure 7. Discrepancies of 3 resections at GCPs
ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume V-1-2020, 2020 XXIV ISPRS Congress (2020 edition) In figure 7 the discrepancies at GCPs of the resection of three neighboured afterward CORONA images are shown. The character of the discrepancies (size and direction at corresponding image position) is similar. It has to be respected that the identification of the GCPs is very difficult due to the object changes over 48 years. Mainly road and railway crossings have been used as GCPs, but 48 years ago the road crossings have been made for animal drawn carriage and today for cars. In addition the roads are wider and the centre line has been changed in unknown direction. Also the railroads changed from single to double tracks and it is not clear where the old track is located in relation to the new one. Another problem is the today coverage by trees what was not the case 48 years ago.
After eliminating approximately 5% of the image points due to poor identifications and changed ground locations, based on image coordinates of the joint sub-images, corrected by the sidelines and transformed from panoramic to perspective geometry, resections yielded in Sigma0 values between 45 µm and 121 µm, corresponding to standard deviations in X (SX) between 23m and 42m and in Y (SY) between 9m and 18m, in one case even up to 38m, for the forward and the backward cameras. With regard to ~2m GSD, the discrepancies at GCPs are unsatisfactory and cannot be explained by the limited GCP identification. Despite of the changing size, the shape of the systematic image errors for the afterward camera is similar, so that the image coordinate residuals of all 10 used afterward images were used together to generate an image coordinate correction grid.  10). This common correction grid for the afterward camera improved the resections by approximately 22% to SX=23m and SY=11m. In the central 50% of CORONA images with ~ 1.85m GSD, the accuracy with SX=17m and SY=10m is slightly better than in the outer parts with ~ 2.2m GSD. This is still not optimal, so an additional perspective transformation of the image parts to the GCPs has been used, whereby SX and SY have been improved to approximately 10m in both centre image parts and to approximately 15m at both outside image parts. It was more difficult with the images taken with the forward camera. Only the first 4 images showed a similar deformation (figure 13), which allowed the same procedure as with the afterward images. One of the remaining images also showed a strong curvature in spite of the sideline correction and the other images had a different behaviour. Correction grids based on only one image had to be used (figures 11 and 12), which led to the final accuracy mentioned above, which was achieved with the afterward images. Nevertheless also for the forward camera some basic similarities between the correction grids exists (figures 11 up to 13). The tilted view of the panoramic images causes an angular affine deformation of corresponding windows for matching. This is respected by least squares matching. In addition the grey value variation by the film corn and the missing of sharp corners causes problems for semi global matching, so finally an area based least squares matching was used for the DSM generation. The height models generated by least squares matching have a good relative accuracy in the range of a standard deviation (SZ) of 3m to 4m corresponding to the simple accuracy estimation SZ=h/b * Spx (height to base multiplied with the standard deviation of the x-parallax) with 3.7m for the height to base relation of 1.85 and an accuracy of the x-parallax of ~2m. The morphologic details of the generated DSM corresponded. However the absolute accuracy is still limited to SZ=25m to 30m what can be explained by remaining systematic image errors of the individual images and the influence of limited identification of the GCPs. Even if the morphologic details are more important for the project, the absolute height deviation disturbs. For this reason, the CORONA height models have been improved by TDM90 (Aldosari, Jacobsen 2019). TDM90 is a free version of the TanDEM-X-height model (WorldDEM) with a point spacing of approximately 90m. WorldDEM currently is the world's most accurate height model. For flat areas the height accuracy of WorldDEM is in the range of 1m up to 1.5m (Wessel et al. 2018), (Aldosari, Jacobsen 2019). To improve the CORONA height model, all height differences to TDM90 were computed. The individual CORONA height points were corrected for the average of the height differences in a moving window of 1km x 1km; this lead to an absolute height accuracy of the improved CORONA height model in the range of 3m to 5m without loss of morphologic details.  Figure 14 shows the colour coded eastern part (part 4 in figure  6) of a CORONA DSM based on least squares matching and the orientation using the correction grids (same colour scale as figure 17). Figure 15 presents the colour coded height differences against TDM90 (same colour scale as figures 16 and 18). Between the CORONA DSM and TDM90 there is a bias of -7.70m, a RMSZ of 11.64m and a standard deviation of 8.73m.
Such differences are typical for the whole project. Based on the Hannover program ZFIT the CORONA DSM has been fitted to TDM90. The CORONA DSM has a point spacing of 6m while TDM90 has approximately 90m point spacing. By ZFIT the morphologic details of the CORONA DSM are not influenced (figure 21), but the absolute height values have been changedsee colour coded correction values of figure 16. The corrected CORONA DSM is shown in figure 17 this corresponds to the expectation. The colour coded height differences of the corrected CORONA DSM against TDM90 are shown in figure  18 with the same colour scale as figure 15. Only few small red spots can be seen (close to 20m height differences) caused mainly by trees in the CORONA image which are not anymore available today. After fitting the CORONA DSM to TDM90 a bias of 0.18m and a standard deviation of 2.90m are shown for the difference of the corrected CORONA image and TDM90. The height differences of the uncorrected CORONA DSM and TDM90 (image 15) in general is smaller for both central parts (part 2 and 3 in figure 6) as for the outside parts 1 and 4. The required morphologic details of the CORONA DSM show very well the location of former water courses. If a new building site is generated, the ground is filled up with an approximate 7m thick layer of sand, hiding former water courses. New houses are built approximately 3 years after filling the area with sand when nobody knows something about the location of former water courses which are not a stable building ground. The corresponding images of a CORONA stereo pair are not fitting 100% as it can be seen with the combination of the forward image 60 and the afterward image 66 in figure 19. In addition figure 20 shows the very small overview of the footprints in the centre parts. Due to this situation CORONA height models of neighboured stereo models have gaps in the centre area. This problem has been solved by crossing stereo combinationsin addition to the combination of forward image 60 and the backward image 66 ( figure 19) also the combination of forward image 59 and backward image 66 was used. So the whole area could be covered by DSMs without gaps. Gaps in the height models could not be avoided in areas covered by water, caused by the frequent flooding in Bangladesh. The influence of buildings and trees can be removed by filtering the DSMs to DTMs. This is not a problem due to the fact that obstacles in the DSMs being usually only single or small groups of trees or buildings. Forest is not available in a radius of 200km around Dhaka. Nevertheless also TDM90 is a DSM. For small parts of the project area also WorldDEM DSMs and DTMs have been bought. WorldDEM is the commercial version of the TerraSAR-X height model with 10m point spacing. From WorldDEM DSM the lower resolution TDM90 is generated. It has been shown in the project area that the filtering of WorldDEM DSM was leading to more detailed DTMs as the commercial WorldDEM DTMs. The image matching has been done separately for the four image parts (figure 6) due to the large size of the image parts of approximately 33000 x 12000 pixels.

GENERATION OF ORTHO IMAGES
Not only height models are important for urban planning. Often ortho images show more important details. As reference height model not only CORONA DSM can be used, it is satisfying to use TDM90, or even in some areas a constant terrain height. For the 15° nadir angle of the image centre and the up to 35° half field of view of the panoramic camera, or maximal 38° nadir angle 2.5m height error leads to a positional error of 2m, corresponding to 1 GSD. Nevertheless as mentioned above, the orientation accuracy is influenced by systematic image errors, so TDM90 is accurate enough for the project. As the image matching, the ortho images have been generated for the image parts. Only few image formats can handle the large size of corresponding ortho images with ~33000 pixels across flight direction. As mentioned above, the GSD of the panoramic images varies between 1.77m in the centre and 2.18m x 2.69m at the image ends. As compromise 2m GSD have been used for the ortho images. Of course this leads to a slightly reduced image quality at both extreme film ends, but this can be accepted. The used program estimates the ground coordinates of the image corners to support the optimal corner positions for the ortho images. For the generation of ortho images from the defined ortho pixel with the interpolated height in the national coordinate system, a transformation into geographic, geocentric to the tangential coordinate system was required, followed by transformation into the perspective image and with a panoramic correction, sideline correction, systematic image error correction to the final CORONA pixel. This is time-consuming, so that the exact computation was only carried out for the corners of a selectable patch and within the patch by perspective transformation. No commercial software was available for this topic, requiring an own development.
Due to the limited accuracy of the orientation in the range of 10m to 15m, a fitting of neighboured ortho images was required to avoid misfits in the overlay area of neighboured images.

CONCLUSION
CORONA images are optimal for the reconstruction of the situation 48 years ago for which at least in Bangladesh has no geo-referenced information as base of urban planning to avoid unstable building ground. The geometry of the scanned panoramic images with a total size of approximately 75cm x 5.5cm does not correspond to the accuracy available by modern optical satellites and aerial cameras, requiring a correction based directly or indirectly on ground control points by a correction grid. The footprint size of in total ~220km length demands a handling in an orthogonal coordinate system, realized internally by a tangential plane to the earth ellipsoid. So finally a horizontal standard deviation of 10m to 15m can be reached for the ortho images. The height models based on CORONA stereo pairs have a very good relative accuracy and high morphologic quality, nevertheless the absolute accuracy is still influenced by local film deformations, which was compensated by fitting to the high absolute accuracy of the free available TDM90. With 10 stereo pairs approximately 25% of the area of Bangladesh has been covered and more nearly cloud free CORONA images are available on the USGS Earth Explorer.