Data Fusion of Scanned Black and White Aerial Photographs with Multispectral Satellite Images

: To date, countless satellite image fusions have been made, mainly with panchromatic spatial resolution to a multispectral image ratio of 1 / 4, fewer fusions with lower ratios, and relatively recently fusions with much higher spatial resolution ratios have been published. Apart from this, there is a small number of publications studying the fusion of aerial photographs with satellite images, with the year of image acquisition varying and the dates of acquisition not mentioned. In addition, in these publications, either no quantitative controls are performed on the composite images produced, or the aerial photographs are recent and colorful and only the RGB bands of the satellite images are used for data fusion purposes. The objective of this paper is the study of the addition of multispectral information from satellite images to black and white aerial photographs of the 80s decade (1980–1990) with small di ﬀ erence (just a few days) in their image acquisition date, the same year and season. Quantitative tests are performed in two case studies and the results are encouraging, as the accuracy of the classiﬁcation of the features and objects of the Earth’s surface is improved and the automatic digital extraction of their form and shape from the archived aerial photographs is now allowed. This opens up a new ﬁeld of use for the black and white aerial photographs and archived multispectral satellite images of the same period in a variety of applications, such as the temporal changes of cities, forests and archaeological sites.


Introduction
Data fusion is the result of using two or more images and incorporating their content information in order for the new composite image to contain more information than can be originally captured by the sensors. Image synthesis methods and techniques are used to create a high spatial resolution image that attempts to maintain the spectral information of the lower spatial resolution original data. In the composite image, the accuracy of the geometric correction can be improved, the intertemporal changes can be better defined and in addition optimal visual interpretation and classification can be allowed. Some wider areas of application are cartography, environment, urban planning, town planning, etc. [1][2][3][4][5][6][7][8].
Black and white (B/W) aerial photographs (film of high sensitivity of 250 lp/mm) of 1987 and 1990 of Sparta and Pyrgos area (Figures 1 and 2) were used respectively, at a scale of 1:20000 (Table  1), which were scanned at 1200 dpi. In both areas, urban and sub-urban areas are included in aerial photographs, while the selection of locations was the result of the first searches of the editorial team for data with small difference in their acquisition date (just a few days), the same year and season. In each study area the aerial photographs belong to one strip and have an overlap of 60%. All aerial photographs were accompanied by camera calibration data (Camera Calibration RMK A 15/23 and Camera Calibration RMK A). Moreover, atmospheric corrected satellite images Landsat 5 of 1987 ( Figure 3) and 1990 were collected in the respective geographic areas (U.S. Geological Survey [20]). Aerial photographs of 1987 were captured 7 days before, while 1990 aerial photographs were taken 1 day after the satellite images were captured in the above geographic areas. Finally, for the realization of the aerial triangulation, Ground Control Points (GCPs) and Check Points (CPs) were used, whose horizontal coordinates were collected from the Hellenic Cadaster [21], ie. the official public provider of cartographic base maps in Greece, with horizontal accuracy of 1 m. The corresponding altitude information was collected from the DTM (5 × 5 m point grid, 1.5 m vertical accuracy) also from the Hellenic Cadastre.    Black and white (B/W) aerial photographs (film of high sensitivity of 250 lp/mm) of 1987 and 1990 of Sparta and Pyrgos area (Figures 1 and 2) were used respectively, at a scale of 1:20000 (Table  1), which were scanned at 1200 dpi. In both areas, urban and sub-urban areas are included in aerial photographs, while the selection of locations was the result of the first searches of the editorial team for data with small difference in their acquisition date (just a few days), the same year and season. In each study area the aerial photographs belong to one strip and have an overlap of 60%. All aerial photographs were accompanied by camera calibration data (Camera Calibration RMK A 15/23 and Camera Calibration RMK A). Moreover, atmospheric corrected satellite images Landsat 5 of 1987 ( Figure 3) and 1990 were collected in the respective geographic areas (U.S. Geological Survey [20]). Aerial photographs of 1987 were captured 7 days before, while 1990 aerial photographs were taken 1 day after the satellite images were captured in the above geographic areas. Finally, for the realization of the aerial triangulation, Ground Control Points (GCPs) and Check Points (CPs) were used, whose horizontal coordinates were collected from the Hellenic Cadaster [21], ie. the official public provider of cartographic base maps in Greece, with horizontal accuracy of 1 m. The corresponding altitude information was collected from the DTM (5 × 5 m point grid, 1.5 m vertical accuracy) also from the Hellenic Cadastre.

Methodology and Processing of Data/Products
After collecting the necessary data, the aerial triangulation-s of aerial photographs (for the production of orthophoto mosaics) and geometric correction-s of satellite images are conducted in the study area. Following are the image fusions of the orthophoto mosaics of the aerial photographs with the Landsat 5 orthoimages, as well as the checking of the quality (correlation tables) of the fused spectral information of the Multispectral (MS) satellite. Finally, the following classifications, both on the composite as well as on the MS ortho images, will be checked for their capability to provide accurate area measurements of both built and open surfaces.

Aerial Triangulation
Aerial triangulation was performed with Leica Photogrammetry Suite (LPS) of Erdas Imagine © . Using the Camera Calibration files, the interior orientations of the aerial photographs were restored, 13 GCPs and 5 CPs were selected in the two study areas (Figures 4a and 5a). Digital Surface Model (DSM) of 5 m spatial resolution and the orthophoto mosaics (Figures 4b and 5b) of the study area with spatial resolution of 0.5 m were produced.

Methodology and Processing of Data/Products
After collecting the necessary data, the aerial triangulation-s of aerial photographs (for the production of orthophoto mosaics) and geometric correction-s of satellite images are conducted in the study area. Following are the image fusions of the orthophoto mosaics of the aerial photographs with the Landsat 5 orthoimages, as well as the checking of the quality (correlation tables) of the fused spectral information of the Multispectral (MS) satellite. Finally, the following classifications, both on the composite as well as on the MS ortho images, will be checked for their capability to provide accurate area measurements of both built and open surfaces.

Aerial Triangulation
Aerial triangulation was performed with Leica Photogrammetry Suite (LPS) of Erdas Imagine © . Using the Camera Calibration files, the interior orientations of the aerial photographs were restored, 13 GCPs and 5 CPs were selected in the two study areas (Figures 4a and 5a). Digital Surface Model (DSM) of 5 m spatial resolution and the orthophoto mosaics (Figures 4b and 5b) of the study area with spatial resolution of 0.5 m were produced.   After the completion of the aerial triangulations, the CPs were used to calculate the estimates of the differences x Δ and y Δ and the standard deviations X σ and Υ σ ( Table 2).   After the completion of the aerial triangulations, the CPs were used to calculate the estimates of the differences x Δ and y Δ and the standard deviations X σ and Υ σ ( Table 2). After the completion of the aerial triangulations, the CPs were used to calculate the estimates of the differences∆x and∆y and the standard deviationsσ X andσ Y ( Table 2).
where δx i the differenced of CPs in the X axis between the orthoimage and the actual values, x ORTHO,i the values of CPs in the X axis in the orthoimage, x CP,i the actual values of CPs in the X axis, and n the number of observations (=5).

Geometric Correction of Satellite Images
For the collection of GCPs, in order to perform the geometric correction of satellite images (image resection using Erdas Imagine © ), the products of the aerial photography were utilized, namely the orthophoto mosaics and DSMs, were used. Substantially the image (of the satellite) was recorded on another image (orthophoto mosaic of aerial photographs). In each study area, 10 GCPs and 5 CPs were used. After the geometric corrections, the CPs were used to calculate the estimates of the differenceŝ ∆x and∆y and the estimates of the standard deviationsσ X andσ Y ( Table 2).

Fusion of Images
The fusion of Orthophoto mosaics of the aerial photographs with the Landsat 5 orthoimages, as well as the checking of the quality of the fused spectral information of the Multispectral (MS) satellite, was carried out in Erdas Imagine © . Figures 6 and 7 show the b/w orthophoto mosaics (Figures 6a and 7a) and the multi-spectral satellite images (Figures 6b and 7b) in rectangular sections in the study areas. Initially, Resolution Merge was performed using the Principal Component Analysis method, Bilinear Interpolation and Output 8 bit data reconstruction technique, producing a data-fusion image for each study area (Figure 6c,d and Figure 7c,d). The evaluation of fused image is based on qualitative-visual analysis and quantitative-statistical analysis. The qualitative-visual analysis is subjective and is directly related to the experience of the fused image creator (e.g., are more details recognized in the image or are colors, contrasts preserved? etc.) [7]. The quantitative-statistical analysis is objective and is based on spectral analysis of images. The most commonly used method is the correlation coefficient between the original bands of the MS image and the corresponding bands of the fused image. The correlation coefficient values range from −1 to 1. Usually, the values between the corresponding bands of the two images (MS and fused image) must be from 0.9 to 1, so that the fused image can be used for the e.g., successful classification of earth's surface coverings and objects [7,[22][23][24]. In order to create the correlation tables (Tables 3 and 4) of ortho multispectral satellite images with composite images, composite images were degraded spatially (60 times, image degradation) to obtain the spatial resolution of MS images.
correlation coefficient values range from -1 to 1. Usually, the values between the corresponding bands of the two images (MS and fused image) must be from 0.9 to 1, so that the fused image can be used for the e.g. successful classification of earth's surface coverings and objects [7,[22][23][24]. In order to create the correlation tables (Tables 3 and 4) of ortho multispectral satellite images with composite images, composite images were degraded spatially (60 times, image degradation) to obtain the spatial resolution of MS images.  correlation coefficient values range from -1 to 1. Usually, the values between the corresponding bands of the two images (MS and fused image) must be from 0.9 to 1, so that the fused image can be used for the e.g. successful classification of earth's surface coverings and objects [7,[22][23][24]. In order to create the correlation tables (Tables 3 and 4) of ortho multispectral satellite images with composite images, composite images were degraded spatially (60 times, image degradation) to obtain the spatial resolution of MS images.

Classifications and Area Measurements
The classifications that follow, both on the composite as well as on the MS ortho images, will be checked for the capability of providing correct area measurements of both built and open surfaces, into smaller geographical areas (Figure 6c compared with Figures 8a and 7c with Figure 9a). The use of smaller geographic areas will avoid generalizations more apparent in geographically larger areas. As reference data for and comparison reasons, the results of photo interpretation and manual rendering of complex images will be used.
For the clarification of composite and MS satellite orthoimages [25], the Unsupervised -a Pixel-based technique-and ISODATA as a classifier were used. 35 classes (by estimation) of classification were selected, then grouped into regions of built and open surface (Figure 8c,d and Figure 9c,d) and, finally, their areas were automatically calculated ( Table 5). The above processing was performed with Erdas Imagine © .
In order to determine the actual areas of the built and open surface, composite images were introduced into ArcGIS © , where the digitization (Figures 8e and 9e) and the automated calculations of their areas (

Classifications and Area Measurements
The classifications that follow, both on the composite as well as on the MS ortho images, will be checked for the capability of providing correct area measurements of both built and open surfaces, into smaller geographical areas (Figure 6c compared with Figure 8a and Figure 7c with Figure 9a). The use of smaller geographic areas will avoid generalizations more apparent in geographically larger areas. As reference data for and comparison reasons, the results of photo interpretation and manual rendering of complex images will be used.
For the clarification of composite and MS satellite orthoimages [25], the Unsupervised -a Pixelbased technique-and ISODATA as a classifier were used. 35 classes (by estimation) of classification were selected, then grouped into regions of built and open surface (Figures 8c,d and 9c,d) and, finally, their areas were automatically calculated ( Table 5). The above processing was performed with Erdas Imagine © .
In order to determine the actual areas of the built and open surface, composite images were introduced into ArcGIS © , where the digitization (Figures 8e and 9e) and the automated calculations of their areas (Table 5) were performed.

Results
The wider areas of two Greek cities are studied in the paper, in order to check the effectiveness of the composition process of satellite images with b/w aerial photographs with small difference (just a few days) in the image acquisition date, the same year and season, in areas with diversified coverages.
The geometric accuracy of the produced orthophoto mosaic and orthoimages is presented in Table 2. In particular, the achieved results in terms of standard errors on the X and Y axis range 0.4-1.8 m (estimates of the differences 1.0-2.4 m) in the case of orthophoto mosaics of aerial photographs and 3.6-5.7 m (estimates of the differences 7.2-9.5 m) in the case of satellite orthoimages.
Satellite image providers supply MS images with a spatial resolution four times larger than the panchromatic. This was simulated by producing the orthophoto mosaics of the aerial photographs spatially degraded from 0.5 m to 7.5 m (four times smaller than the spatial resolution of the MS image)

Results
The wider areas of two Greek cities are studied in the paper, in order to check the effectiveness of the composition process of satellite images with b/w aerial photographs with small difference (just a few days) in the image acquisition date, the same year and season, in areas with diversified coverages.
The geometric accuracy of the produced orthophoto mosaic and orthoimages is presented in Table 2. In particular, the achieved results in terms of standard errors on the X and Y axis range 0.4-1.8 m (estimates of the differences 1.0-2.4 m) in the case of orthophoto mosaics of aerial photographs and 3.6-5.7 m (estimates of the differences 7.2-9.5 m) in the case of satellite orthoimages.
Satellite image providers supply MS images with a spatial resolution four times larger than the panchromatic. This was simulated by producing the orthophoto mosaics of the aerial photographs spatially degraded from 0.5 m to 7.5 m (four times smaller than the spatial resolution of the MS image) and then the fusion of the images was carried out. From the Correlation Tables that created, it is shown that the rates of retaining the original spectral information of the ortho MS images in the composite images are similar (at least if not better) to the corresponding percentages of the non-initial degradation of the orthophoto mosaics of the aerial photographs presented in this paper. For this reason, the case of the initial degradation of orthophoto mosaics was not further analyzed in the paper.
In addition, during the fusion of the images (with or without initial degradation of the orthomosaics of aerial photographs) apart from the Principal Component Analysis Method (Transformation Based Fusion), other methods were used, such as Multiplacative and Brovey Transform (Additive and Multiplicative Technique) [6,[26][27][28][29][30][31][32][33][34], which did not provide better results in the maintenance of spectral information, and, therefore, were not further analyzed in the paper.
According to the Correlation Table (Table 3) for the wider area of the city of Sparta, the transmission of the spectral information from the MS satellite orthoimage to the composite image is performed satisfactory by 70-85% (except for NIR that is 53%). Also, for the wider area of the city of Pyrgos, it is found that the transmission of the spectral information from the MS satellite orthoimage to the composite image is, also, performed satisfactory by 74-82%. Consequently, despite the fact that the rates are not from 90-100% (ideally), these are generally satisfactory and composite images may be used to perform further classifications. Table 5 shows that in the case of the composite image, the overall estimate of built surfaces is about 20% underestimated, while the overall estimate of open surfaces is about 6% over-estimated. In the case of the original MS satellite images, the overall estimate of built surfaces is approximately 31% over-estimated, while the overall estimate of open surfaces is about 9% underestimated. The important conclusion is that the calculated areas of both the built and the open surfaces in the case of classification of composite images are much closer to reality. In addition, it is of major importance that in the case of composite image classification, it is possible to automatically extract the geometry/shape of the earth's objects (compare Figure 8c Figure 9e). The above is obviously due to the very good spatial resolution of composite images with respect to the lower spatial resolution of Landsat 5 images.

Discussion
Few publications have been made on data fusion of aerial photography with satellite images. In addition to their scarcity, the acquisition date of aerial photographs differs from the acquisition date of satellite images from 2 to 4 years. This fact clearly results in differences-small or large-in the characteristics and objects of the earth's surface that are documented, triggering difficulties and distortions during their study. Moreover, the acquisition dates are not mentioned, despite their importance, as different acquisition periods create errors in the composition of images. Lastly, three issues emerge, not necessarily at the same time: only the RGB bands are used for satellite image fusion, no additional quantitative controls are conducted and the aerial photographs are recent (21th century) and colorful.
No new technique for image composition is presented in the paper, but for the first time in international literature, the addition of multispectral information from satellite images to b/w aerial photographs with small difference in their image acquisition date (just a few days), the same year and season is presented.
In fact, in several cases of composing modern PAN and MS images of the same satellite system, the acceptable limits of retaining the original MS image spectral information in the composite image are either not satisfied or marginally satisfied. Even more so, if images come from such different sensors as in this paper, it is expected that the above limits could not be met. Therefore, it is of primary interest to determine the percentage of spectral information retained in the composite image, and whether this percentage permits the classification of new images, the digital identification of basic structures of the terrestrial surface (built area and non-built area) and the extraction of quantitative data. The results have allowed the above to reach a satisfactory level of acceptance and, in consequence, the research on the creation of new techniques and methodologies for their improvement can now begin.

Conclusions
The paper highlights a new possibility of using archived b/w aerial photographs, since their spectral information can be improved. It is a prerequisite that the corresponding spatial multi-spectral satellite data have the same date of acquisition or differ for just a few days from the acquisition dates of the aerial photographs. New possibilities of using archived Landsat 5 satellite images are also highlighted.
Flights for b/w aerial photographs, in the past, aimed at studies on Cadastre, Forestry, Town Planning, Spatial Planning, etc. Each flight was accompanied by the acquisition of hundreds of aerial photographs, which, nowadays, will allow the enrichment with multi-spectral data of orthophoto mosaics that map spatially large areas. Therefore, automatic extraction of temporal information, e.g., of built and open areas on hundreds of square kilometers of land, will no longer be the result of continuous photointerpretation process and rendering of orthophoto mosaics, but of their automatic extraction through automatic classification. Finally, a variety of new applications is open, such as the identification of the temporal changes of cities, forests, archaeological sites, etc.
Obviously, the data fusion of high spatial resolution images with low spatial resolution ones, as in our case, is accompanied by the problem of large ratios of their spatial resolutions, which often causes the appearance of spectral distortions. However, as the first results are encouraging, the development of optimizations to minimize these spectral distortions are now possible.