Rigorous Co-Registration of KOMPSAT-3 Multispectral and Panchromatic Images for Pan-Sharpening Image Fusion.

KOMPSAT-3, a Korean earth observing satellite, provides the panchromatic (PAN) band and four multispectral (MS) bands. They can be fused to obtain a pan-sharpened image of higher resolution in both the spectral and spatial domain, which is more informative and interpretative for visual inspection. In KOMPSAT-3 Advanced Earth Imaging Sensor System (AEISS) uni-focal camera system, the precise sensor alignment is a prerequisite for the fusion of MS and PAN images because MS and PAN Charge-Coupled Device (CCD) sensors are installed with certain offsets. In addition, exterior effects associated with the ephemeris and terrain elevation lead to the geometric discrepancy between MS and PAN images. Therefore, we propose a rigorous co-registration of KOMPSAT-3 MS and PAN images based on physical sensor modeling. We evaluated the impacts of CCD line offsets, ephemeris, and terrain elevation on the difference in image coordinates. The analysis enables precise co-registration modeling between MS and PAN images. An experiment with KOMPSAT-3 images produced negligible geometric discrepancy between MS and PAN images.


Introduction
KOMPSAT-3, a Korean earth observing satellite, provides four multispectral (MS) bands (i.e., blue, green, red, near infrared) and one panchromatic (PAN) band. The image fusion or the pan-sharpening combines MS and PAN images into a single image that increase the spatial resolution while simultaneously preserving the spectral information in the MS images. This leads to high spectral resolution and high spatial resolution for various applications, such as visual inspection to identify the textures and shape of the various objects, feature extraction, and map updating. A number of pan-sharpening techniques have been studied, including intensity-hue-saturation (IHS), high-pass filtering (HPF), principal component analysis (PCA), Brovey and wavelet transforms, and Contourlet transform [1].
Co-registration is the most critical pre-processing requirement for the fusion of the MS and PAN images, and its accuracy is an important parameter of the fusion product quality [2]. Because MS and PAN CCD (the charge-coupled device) lines are placed with certain offsets, each CCD line captures the same ground target at different times. Under the influence of satellite with the ephemeris and terrain effects remaining, the MS and PAN images cannot be precisely co-registered to each other by a simple image shift. The inaccurate co-registration also impacts further applications like data fusion, change detection, and spectral-signature-based classification.
Conventional image co-registration has been studied out in two ways. First, the image-matching approach performs sub-pixel matching between MS and PAN images, providing tie points that enable Figure 1a shows the configuration of KOMPSAT-3 AEISS camera. Blue, PAN1, PAN2, green, red, and near-infrared (NIR) channels are aligned in a unifocal camera of focal length, 8.56215983 meters. Figure 1b depicts the design of PAN and MS CCD alignments. The gaps between sensors in the focal plane correspond to the differences in the projection centers [13]. The camera system is designed to provide a PAN resolution of 0.7 m and an MS resolution of 2.8 m, as presented in Table 1

Physical Sensor Modeling
KOMPSAT-3 AEISS is a pushbroom sensor based on a non-linear projection form from a given ground point in an earth-centered, earth-fixed (ECEF) coordinate frame to a point in the body coordinate frame, as shown in Equations (1) and (2) [13]. The ephemeris data are provided from KARI, and they are interpolated for the exterior orientation parameters (EOPs) for a given time. Figure 2 depicts the MS and PAN CCD lines acquiring the same ground target with an offset in the corresponding projection centers. This produces row coordinate differences in MS and PAN images (approximate equation equal to the corresponding ground distance of the offset/spatial resolution in KOMPSAT-3). Note that the amount of row coordinate difference depends on the different view directions.

Physical Sensor Modeling
KOMPSAT-3 AEISS is a pushbroom sensor based on a non-linear projection form from a given ground point in an earth-centered, earth-fixed (ECEF) coordinate frame to a point in the body coordinate frame, as shown in Equations (1) and (2) [13]. The ephemeris data are provided from KARI, and they are interpolated for the exterior orientation parameters (EOPs) for a given time. Figure 2 depicts the MS and PAN CCD lines acquiring the same ground target with an offset in the corresponding projection centers. This produces row coordinate differences in MS and PAN images (approximate equation equal to the corresponding ground distance of the offset/spatial resolution in KOMPSAT-3). Note that the amount of row coordinate difference depends on the different view directions.  The body coordinates can be converted from the image coordinates using the CCD line alignment information parameters in Equation (3), which is a second-order polynomial. The variable is the column coordinate in pixels, because the KOMPSAT-3 AEISS sensor is a pushbroom sensor of single CCD line vertical to the flight direction. The parameters are determined based on the precise camera calibration before the launch as well as the in-flight calibration with a number of ground control points [14][15][16].
where 0 2 0 2 ,ã a b b are CCD line alignment parameters, and samp is the column (sample) coordinate in pixels. Figure 3 shows the flowchart of the rigorous MS and PAN image co-registration in this study. From PAN and MS images, a grid of conjugate points is generated. Note that the conjugate points are not generated using image matching, but their coordinates are computed based on sensor modeling. Given a flat terrain, iterative physical sensor modeling is carried out between PAN and each of the four MS images. For example, p' on MS in Figure 2 is projected onto the ground surface and then projected into PAN for p". Similarly, grid points on MS can be projected into PAN for conjugate points. The coordinate differences between the points on MS and PAN are due to the differences in the sensor alignment, ephemeris effects, and terrain elevation. These effects can be modeled as follows.

Rigorous Co-Registration of MS and PAN Images
First, the differences in conjugate point coordinates are modelled to estimate the constant image offsets and linear coordinate differences. This model compensates for the coordinate differences mainly from CCD offsets between MS and PAN. Second, the coordinate differences from satellite The body coordinates can be converted from the image coordinates using the CCD line alignment information parameters in Equation (3), which is a second-order polynomial. The variable is the column coordinate in pixels, because the KOMPSAT-3 AEISS sensor is a pushbroom sensor of single CCD line vertical to the flight direction. The parameters are determined based on the precise camera calibration before the launch as well as the in-flight calibration with a number of ground control points [14][15][16].
where a 0 ∼ a 2 , b 0 ∼ b 2 are CCD line alignment parameters, and samp is the column (sample) coordinate in pixels. Figure 3 shows the flowchart of the rigorous MS and PAN image co-registration in this study. From PAN and MS images, a grid of conjugate points is generated. Note that the conjugate points are not generated using image matching, but their coordinates are computed based on sensor modeling. Given a flat terrain, iterative physical sensor modeling is carried out between PAN and each of the four MS images. For example, p' on MS in Figure 2 is projected onto the ground surface and then projected into PAN for p". Similarly, grid points on MS can be projected into PAN for conjugate points. The coordinate differences between the points on MS and PAN are due to the differences in the sensor alignment, ephemeris effects, and terrain elevation. These effects can be modeled as follows.

Rigorous Co-Registration of MS and PAN Images
First, the differences in conjugate point coordinates are modelled to estimate the constant image offsets and linear coordinate differences. This model compensates for the coordinate differences mainly from CCD offsets between MS and PAN. Second, the coordinate differences from satellite attitude changes are analyzed and compensated. Because satellite attitude changes with PAN and MS data acquisition time, the effect of the attitude should be modeled. Third, the effect of terrain elevation variation on the difference in row coordinates is investigated and modelled. We use the sensor modeling via the ground projection between PAN and MS images. The terrain elevation variation introduces Sensors 2020, 20, 2100 5 of 18 coordinate differences that are fixed. Finally, image resampling for the co-registration between PAN and MS is carried out for pan-sharpening, considering all the aforementioned effects.
Sensors 2020, 20, x FOR PEER REVIEW 5 of 19 attitude changes are analyzed and compensated. Because satellite attitude changes with PAN and MS data acquisition time, the effect of the attitude should be modeled. Third, the effect of terrain elevation variation on the difference in row coordinates is investigated and modelled. We use the sensor modeling via the ground projection between PAN and MS images. The terrain elevation variation introduces coordinate differences that are fixed. Finally, image resampling for the coregistration between PAN and MS is carried out for pan-sharpening, considering all the aforementioned effects.

Data
The experiment was carried out for Ulaanbaatar, Mongolia with a terrain variation from 1250 to 1950 m. The KOMPSAT-3 image was acquired in 17 July 2016 during 5.17 s with roll, pitch, and yaw angles of −7.7°, −20.0°, and −0.4°, respectively. The satellite azimuth and incidence angles are 145.9° and 23.9°. Table 2

The Effect of CCC Line Offset on the Image Co-Registration
First, a grid of conjugate points was generated in MS and PAN image spaces. Each PAN image coordinate was projected to the ground using a backward model of Equation (1) and then projected to each of the four MS image coordinates using a forward model of Equation (2). For the image projection, the precisely calibrated inner orientation parameters and EOP parameters from the

Data
The experiment was carried out for Ulaanbaatar, Mongolia with a terrain variation from 1250 to 1950 m. The KOMPSAT-3 image was acquired in 17 July 2016 during 5.17 s with roll, pitch, and yaw angles of −7.7 • , −20.0 • , and −0.4 • , respectively. The satellite azimuth and incidence angles are 145.9 • and 23.9 • . Table 2 shows the CCD line alignment parameters of each of the four CCD lines, where slightly different values are observed along the scan direction (difference in b 0 ∼ b 2 ).

The Effect of CCC Line Offset on the Image Co-Registration
First, a grid of conjugate points was generated in MS and PAN image spaces. Each PAN image coordinate was projected to the ground using a backward model of Equation (1) and then projected to each of the four MS image coordinates using a forward model of Equation (2). For the image projection, the precisely calibrated inner orientation parameters and EOP parameters from the ephemeris were used. Mean elevation flat terrain from shuttle radar topography mission (SRTM) was initially assumed for the iterative projection, because the terrain elevation effect is considered later. Figure 4a shows the grid points in MS and PAN with an interval of 250 pixels. Note that PAN image space is shown in MS image scale. Figure 4b depicts the projection of the image center points in MS onto PAN for conjugate points, where large differences in row coordinate were observed.
Sensors 2020, 20, 2100 6 of 18 ephemeris were used. Mean elevation flat terrain from shuttle radar topography mission (SRTM) was initially assumed for the iterative projection, because the terrain elevation effect is considered later. Figure 4a shows the grid points in MS and PAN with an interval of 250 pixels. Note that PAN image space is shown in MS image scale. Figure 4b depicts the projection of the image center points in MS onto PAN for conjugate points, where large differences in row coordinate were observed.  We computed the discrepancy of conjugate points between MS and PAN images, as plotted in Figure 5. Note that the difference is computed by subtracting PAN from MS coordinates. In Figure  5a, BLUE CCD shows that row coordinate difference (i.e., scan direction offset) varies by about 17.1-18.8 pixels with a parabolic shape on the x-axis of the column. Figure 5b-d represents that GREEN, RED, and Near IR cases also show similar patterns in the graph. Interestingly, BLUE is the only CCD placed upper than PAN CCD as seen in Figure 1b. On the right side in Figure 5, the differences in column coordinates are observable up to 15 pixels in straight line shapes while the amounts vary with different CCDs. We computed the discrepancy of conjugate points between MS and PAN images, as plotted in Figure 5. Note that the difference is computed by subtracting PAN from MS coordinates. In Figure 5a, BLUE CCD shows that row coordinate difference (i.e., scan direction offset) varies by about 17.1-18.8 pixels with a parabolic shape on the x-axis of the column. Figure 5b-d represents that GREEN, RED, and Near IR cases also show similar patterns in the graph. Interestingly, BLUE is the only CCD placed upper than PAN CCD as seen in Figure 1b. On the right side in Figure 5, the differences in column coordinates are observable up to 15 pixels in straight line shapes while the amounts vary with different CCDs.
First, we derived CCD line offsets from row coordinates difference in the first image column as shown in Table 3. The offset ranges from 19 pixels in BLUE to −28 pixels in Near IR. The remaining coordinate differences between MS and PAN CCD lines were modelled using Equation (4). The row and column coordinate differences are modelled using second and first order polynomials, respectively. Table 3 shows the estimated transformation parameters for each of the CCD lines.
where A, B, C, D, E are transformation parameters for CCD offsets.    Figure 6 and Table 4 show the residual after applying the offsets and transformation parameters to each of the four MS images. Most coordinate differences are less than one-pixel level. This can be further improved because this study plans to achieve a difference of less than 0.1 pixels. First, we derived CCD line offsets from row coordinates difference in the first image column as shown in Table 3. The offset ranges from 19 pixels in BLUE to −28 pixels in Near IR. The remaining coordinate differences between MS and PAN CCD lines were modelled using Equation (4). The row and column coordinate differences are modelled using second and first order polynomials, respectively. Table 3 shows the estimated transformation parameters for each of the CCD lines.
where , , , , A B C D E are transformation parameters for CCD offsets. Table 3. Estimated parameters of the coordinate differences due to the initial CCD offsets.

CCD Initial Offset [pixels]
A  Table 4 show the residual after applying the offsets and transformation parameters to each of the four MS images. Most coordinate differences are less than one-pixel level. This can be further improved because this study plans to achieve a difference of less than 0.1 pixels.

The Effect of Attitude Changes on the Image Co-Registration
Previous forward and backward projections are carried out based on the assumptions of accurate ephemeris and flat terrain. First, we analyzed the errors from attitude changes between two acquisition times. The attitude changes in the same image were computed with respect to the change in row coordinates. For BLUE CCD in Figure 7a, both roll and pitch angles changes range from −0.3 to +0.3 in pixels and yaw angle changes vary from −0.08 to +0.04 in pixels. For GREEN CCD in Figure  7b, different amounts of changes are observable. Note that RED and Near IR CCD cases are omitted here because they show the similar patterns to GREEN CCD. Table 5 presents the statistics of the attitude changes in pixels.

The Effect of Attitude Changes on the Image Co-Registration
Previous forward and backward projections are carried out based on the assumptions of accurate ephemeris and flat terrain. First, we analyzed the errors from attitude changes between two acquisition times. The attitude changes in the same image were computed with respect to the change in row coordinates. For BLUE CCD in Figure 7a, both roll and pitch angles changes range from −0.3 to +0.3 in pixels and yaw angle changes vary from −0.08 to +0.04 in pixels. For GREEN CCD in Figure 7b, different amounts of changes are observable. Note that RED and Near IR CCD cases are omitted here because they show the similar patterns to GREEN CCD. Table 5 presents the statistics of the attitude changes in pixels.   We also investigated the effect of satellite attitudes on the MS and PAN image co-registration. We plotted the coordinate differences in the selected grid points with respect to the changes in roll, pitch, and yaw angles in Figure 8. The limit of angle variations were derived from the maximum and minimum attitude angles for the entire image. A linear pattern between the column coordinates' differences was observed with the change in roll angle. In addition, the row coordinates' difference was correlated with pitch angle. Yaw angle was associated with both row and column differences, but its effect was smaller than the other angles. It is noteworthy that the roll angle change in the pushbroom sensor was orthogonal to the scan direction, whereas pitch and yaw angles were different from the scan direction viewing angles. The standard deviations in the coordinate differences for each of the four MS sensors are presented in Table 6.  We also investigated the effect of satellite attitudes on the MS and PAN image co-registration. We plotted the coordinate differences in the selected grid points with respect to the changes in roll, pitch, and yaw angles in Figure 8. The limit of angle variations were derived from the maximum and minimum attitude angles for the entire image. A linear pattern between the column coordinates' differences was observed with the change in roll angle. In addition, the row coordinates' difference was correlated with pitch angle. Yaw angle was associated with both row and column differences, but its effect was smaller than the other angles. It is noteworthy that the roll angle change in the pushbroom sensor was orthogonal to the scan direction, whereas pitch and yaw angles were different from the scan direction viewing angles. The standard deviations in the coordinate differences for each of the four MS sensors are presented in Table 6.     From the analysis, we formed the compensation equation, Equation (5) for attitude changes between two image acquisition epochs derived from MS and PAN CCD line offsets. In Equation (5), the column coordinate correction was linearly modelled using roll and yaw angles. The row coordinate correction is modelled using pitch and yaw angles. Table 7 shows the estimated parameters of the compensation model for attitude changes for each of the four MS bands where ∆column, ∆row are the differences of column and row coordinates between sub-images, respectively. α, β, γ, η are correction parameters. ∆roll, ∆pitch, ∆yaw are the viewing direction changes in roll, pitch, yaw angles caused by MS-PAN CCD line offsets. We applied this equation for the differences in the row and column coordinates plotted in Figure 9. The differences become less than our acceptable level, 0.1 pixels. The statistics of each of the four MS sensors are presented in Table 8. We applied this equation for the differences in the row and column coordinates plotted in Figure  9. The differences become less than our acceptable level, 0.1 pixels. The statistics of each of the four MS sensors are presented in Table 8.

The Effect of Terrain Elevation Variation on the Image Co-Registration
We assumed a flat target terrain derived from SRTM elevation data in the previous section. Here, we tested the ground elevation change in the image projections with some grid points in MS and PAN image spaces. As the red dots show in Figure 10, elevation changes led to less than 2 pixels of row coordinate mismatches. It was also reported that 2.8 km of terrain elevation variation creates about one pixel of mismatch in Quickbird image [17] and the accuracy of DEM is not very critical for determining the relative mismatch between the bands [18]. We can identify that row coordinate differences are highly correlated with terrain elevation variation. Specifically, high row coordinate differences are observed in both RED and Near IR bands.
(d) Near IR Figure 9. Differences in rows and columns in coordinates between MS and PAN images after attitude change compensation.

The Effect of Terrain Elevation Variation on the Image Co-Registration
We assumed a flat target terrain derived from SRTM elevation data in the previous section. Here, we tested the ground elevation change in the image projections with some grid points in MS and PAN image spaces. As the red dots show in Figure 10, elevation changes led to less than 2 pixels of row coordinate mismatches. It was also reported that 2.8 km of terrain elevation variation creates about one pixel of mismatch in Quickbird image [17] and the accuracy of DEM is not very critical for determining the relative mismatch between the bands [18]. We can identify that row coordinate differences are highly correlated with terrain elevation variation. Specifically, high row coordinate differences are observed in both RED and Near IR bands. Therefore it is required to compensate the terrain elevation effect using Equation (6) Therefore it is required to compensate the terrain elevation effect using Equation (6), where H min , ∆row H min , H max , ∆row H max can be derived from a linear regression in Figure 10.
where ∆row H min , ∆row H max are the row coordinate differences at the minimum and maximum terrain heights H min and H max , respectively. ∆row H is the row coordinate correction value at the terrain height H. The terrain elevation compensation parameters for the test data were estimated using the minimum and maximum elevation, 785 and 2789 m (H min and H max ), derived from SRTM data. The blue dots in Figure 10 show the coordinates compensation derived from Equation (6). Table 9 shows the estimations that are different for each of the MS CCD. Note that the Near IR shows the largest row coordinate difference. We applied the estimated compensation parameters and obtained row coordinate differences smaller than 0.1 pixels between PAN and MS images, as shown in Figure 11.  Figure 11. Differences in row coordinate between MS and PAN images after the terrain elevation compensation.

Resampling
We estimated the parameters of the coordinate differences between PAN and each of the four MS sensors due to the initial CCD offsets, attitude changes, and terrain elevation variations. The estimated parameters were used for the co-registration image resampling of four MS images into a Figure 11. Differences in row coordinate between MS and PAN images after the terrain elevation compensation.

Resampling
We estimated the parameters of the coordinate differences between PAN and each of the four MS sensors due to the initial CCD offsets, attitude changes, and terrain elevation variations. The estimated parameters were used for the co-registration image resampling of four MS images into a PAN image plane for the pan-sharpening. Figure 12 show the result of pan-sharpening in false color before and after the image co-registration.
(c) RED (d) Near IR Figure 11. Differences in row coordinate between MS and PAN images after the terrain elevation compensation.

Resampling
We estimated the parameters of the coordinate differences between PAN and each of the four MS sensors due to the initial CCD offsets, attitude changes, and terrain elevation variations. The estimated parameters were used for the co-registration image resampling of four MS images into a PAN image plane for the pan-sharpening. Figure 12 show the result of pan-sharpening in false color before and after the image co-registration.

Comparison to Image Matching-Based Approach
We compared the rigorous method to image-matching-based co-registrations. For imagematching, Harris operator and Normalized Cross Correlation (NCC) were used for key points' extraction and image-matching between PAN and MS, as the extracted key points are presented in Figure 13. The result shows that it is not easy to extract points over mountainous areas such as the upper image region. The threshold for NCC matching was set to 0.8 and a total of 750 points were matched after outlier removal. Half of the extracted points were used for co-registration modeling and the other half were used for the co-registration accuracy check. We applied two co-registration modelings: the image shift only and the interpolation [19]. Table 10 shows the accuracy results of coregistration at the check points. The shift-only was not able to achieve sub-pixel accuracy. The interpolation and rigorous methods could achieve sub-pixel accuracy, while the rigorous method show slightly better results. It is notable that the rigorous method may not seem accurate, as the above rigorous analysis showed, but it is likely that the check points are not error-free because they were also extracted using image matching.

Comparison to Image Matching-Based Approach
We compared the rigorous method to image-matching-based co-registrations. For image-matching, Harris operator and Normalized Cross Correlation (NCC) were used for key points' extraction and image-matching between PAN and MS, as the extracted key points are presented in Figure 13. The result shows that it is not easy to extract points over mountainous areas such as the upper image region. The threshold for NCC matching was set to 0.8 and a total of 750 points were matched after outlier removal. Half of the extracted points were used for co-registration modeling and the other half were used for the co-registration accuracy check. We applied two co-registration modelings: the image shift only and the interpolation [19]. Table 10 shows the accuracy results of co-registration at the check points. The shift-only was not able to achieve sub-pixel accuracy. The interpolation and rigorous methods could achieve sub-pixel accuracy, while the rigorous method show slightly better results. It is notable that the rigorous method may not seem accurate, as the above rigorous analysis showed, but it is likely that the check points are not error-free because they were also extracted using image matching.
upper image region. The threshold for NCC matching was set to 0.8 and a total of 750 points were matched after outlier removal. Half of the extracted points were used for co-registration modeling and the other half were used for the co-registration accuracy check. We applied two co-registration modelings: the image shift only and the interpolation [19]. Table 10 shows the accuracy results of coregistration at the check points. The shift-only was not able to achieve sub-pixel accuracy. The interpolation and rigorous methods could achieve sub-pixel accuracy, while the rigorous method show slightly better results. It is notable that the rigorous method may not seem accurate, as the above rigorous analysis showed, but it is likely that the check points are not error-free because they were also extracted using image matching.

Conclusions
We studied the effects of CCD offsets, attitude effects, and terrain elevation variation on the KOMPSAT-3 MS and PAN image co-registration. The KOMPSAT-3 CCD lines acquire the same ground features with certain offsets in the corresponding projection centers. Therefore, we compensated for the offsets based on the rigorous physical sensor modeling. The remaining mismatches were corrected using the first-and second-order equations formed from the coordinate difference analysis. It is noteworthy that the attitudes between each CCD data acquisition change when imaging the same target are modelled with a linear combination of roll, pitch and yaw angle changes for the time gap. Row coordinate differences between MS and PAN images from terrain elevation variation can be modelled using a linear compensation equation. Finally, we compensated for all aforementioned effects on the KOMPSAT-3 MS and PAN image co-registrations with negligible discrepancy, less than 0.1 pixels. A rigorous co-registration approach is more robust and useful with available ephemeris data, whereas image matching-based co-registrations are less reliable over a monotonous terrains such as desert and forest. In further studies, more co-registration methods based on image matching can be analyzed to compare physical sensor model approaches.