The use of high-resolution satellite images has expanded due to developments in satellite imaging technology, which has raised demand for high-precision geospatial data for Earth observation. To conduct land cover analysis and object change detection, accurate geographic information is considered essential in various fields, such as environmental monitoring, urban development, and disaster prevention and management. The initial positional errors in raw Level-1 Radiometric (L1R) satellite images must be corrected to generate Level-2R (L2R) images for enhancing geometrical accuracy. Geometric correction of such high-resolution satellite images is typically performed using ground control points (GCPs) (Yoon et al., 2018; Park et al., 2020; Son et al., 2021). However, the collection of GCPs and the geometric correction for satellite images require significant time and cost. This process poses challenges for constructing large volumes of high-resolution satellite images rapidly. Additionally, obtaining the precise ground coordinates of GCPs is subject to data security or policy issues, as well as geographical constraints. These regulations may restrict access to certain areas, making it difficult to obtain GCPs (Oh et al., 2022). Consequently, it is often inevitable to encounter images that are either not geometrically corrected or have suboptimal correction performance.

To overcome these limitations, recent studies have been conducted on methods for performing geometric corrections without relying on GCPs (Ma et al., 2017; Yang et al., 2017). Among these, relative geometric correction methods utilizing bundle block adjustment have been proposed recently (Ban and Kim, 2024). Bundle block adjustment corrects multiple images simultaneously to reduce positional errors and allows relatively accurate position estimation without GCPs (Yang et al., 2018). However, it becomes difficult to construct stable weight values when images with inconsistent initial positional errors are used. In particular, it is difficult to determine the weight parameters in bundle block adjustment due to the initial positional differences between the L1R images of the KOMPSAT-3/3A satellites before geometric correction processing (Shan et al., 2014). Even if the relative geometric correction between L1R images is performed normally, the actual ground position of the images deviates still.

In this study, we propose a novel method to efficiently perform geometric corrections on uncorrected KOMPSAT-3/3A L1R images using orthorectified reference images. This method aims to address the challenges of performing absolute geometric corrections on individual images. At the same time, we aim to improve the initial position accuracy of relative geometric correction across multiple images. To achieve this, we utilized accessible reference data with high accuracy, such as freely available orthorectified CAS500-1 images. The bundle adjustment was performed using L1R/L2R images generated from KOMPSAT-3A and CAS500-1 Level-1/2 Geometric (L1G/L2G) images. This procedure validated the inverse orthorectification experiment and ensured consistency in the data used.

The main idea of the proposed approach is to use reference image-based bundle block adjustment to correct a large number of uncorrected images with different initial position errors. In this process, we improve the overall absolute geospatial accuracy of the dataset while maintaining positional consistency between images without using GCPs. We have confirmed that stable geometric correction is possible even without the precise GCPs required by existing methods. Consequently, this study presents an effective alternative for generating practical and reliable relative geometric correction images, even in regions where acquiring GCPs is challenging.

The proposed method is carried out in two main steps: The first step is generating rational polynomial coefficients (RPCs) and simulated images by inverse orthorectification and the second step is precise bundle adjustment based on RPCs. First of all, inverse orthorectification is performed using the reference L2G orthoimage. Through reference images, L2R simulated images and RPCs information are generated. This allows other uncorrected images to be referenced in the same coordinate system as the reference image. It is possible to adjust relative positioning mistakes using the precise coordinate data from the reference image. Subsequently, the RPCs correction coefficients of satellite images through bundle block adjustment are re-estimated to provide precise positional accuracy for multiple satellite images.

2.1. Generation of RPCs and Simulated Images by Inverse Orthorectification

The CAS500-1 L2G images provided for free by the National Geographic Information Institute (NGII) have high positional accuracy (National Geographic Information Institute, 2024). The images are provided with precise geometric correction and orthorectification. However, there is a possibility of positional errors occurring when georeferencing L2R images that include geometric distortions, as only orthoimages with known ground coordinates and corrected geometric distortions are provided. Thus, we performed inverse orthorectification on the given CA500-1 orthoimages in this study for use as input data. This process was conducted to generate L2R-simulated images that replicated the geometric features at the moment of capture. Through this procedure, inverse projection based on the collinearity requirement was performed using RPCs to reflect the relationship between ground and image coordinates. The digital elevation model (DEM) was used for the required elevation values in the inverse projection process.

In the case of KOMPSAT-3/3A L2G images, the transformation from image coordinates to ground coordinates can be calculated by considering the following three factors: considering the X and Y coordinates of the top-left corner, the ground sampling distance (GSD) in the row and column directions, and the rotation elements. The relationship between image coordinates and ground coordinates for the CAS500-1 is also defined in the same way as for KOMPSAT-3/3A. The metadata provided with the L2G images from the CAS500-1 includes the corner point coordinates and geometric information of the images at the time of capture. Using this, we generate the RPCs information for L2R and establish the rational function model (RFM). The Eq. (1) below represents the transformation relationship between the ground coordinates and the image coordinates described by the RFM. To enhance the stability of model establishment, each coordinate information is normalized. Based on this, an RFM in the form of a cubic polynomial is constructed as shown in Eq. (2).

In the above equation, $X_n^i$, $Z_n^k$ represent the ground coordinates of the normalized satellite images, and a_ijk corresponds to the generated RPCs information. Using the established higher-order polynomials P₁~P₄, the image coordinates Line_n and Sample_n are calculated.

Ray-tracing method are applied to generate a simulated image by restoring the shooting geometry based on the established RFM. From the projection center of the image, a ray is used as a reference to estimate the L2G pixel values corresponding to the L2R pixel values. The estimated L2R values are used to generate the simulated image. In the initial processing, the position where the ray emitted from the sensor intersects with the ground surface is back-projected using the collinearity condition based on an arbitrary altitude value. Afterward, based on the altitude data obtained from the DEM for the corresponding location, the reference plane is iteratively estimated until the ground coordinates converge. Fig. 1 visually shows the process of generating simulated images through ray tracing.

Figure 1. Generation of simulated images through inverse orthorectification.

Since the estimated ground coordinates in the simulated image do not completely match the center points of the image pixels, additional brightness value interpolation was performed. Through this, the center point of the orthoimage pixel corresponding to the estimated ground coordinates is calculated, and the brightness value of that center point is referenced. The bicubic interpolation method is used to preserve the high-frequency components of the ortho image, aiming to provide high-quality and smooth results.

2.2. Rigorous Bundle Block Adjustment

We aim to perform relative geometric correction on multiple images using the generated L2R simulated image as the reference image. The RPCs-based rigorous bundle block adjustment proposed in this study consists of four main stages, as shown in Fig. 2.

Figure 2. Flowchart of the proposed RPCs block bundle adjustment.

The first step is feature point extraction, where the Scale Invariant Feature Transform (SIFT) algorithm is used to extract feature point coordinates at the sub-pixel level, thereby increasing the accuracy of bundle adjustment. In the second step, the extracted correspondences are used to form observation equations and construct the bundle adjustment matrix. In this stage, the image coordinates of the observation equation correspondences are applied as coefficients to form the matrix for bundle block adjustment. Next, the RPCs correction parameters are estimated through bundle block adjustment. In this process, the weights are re-estimated based on the covariance matrix of the residuals, ensuring a more rigorous approach. In the final step, the corrected output image is generated using the estimated RPCs correction parameters and ground coordinates.

2.2.1. Tie Point Extraction

This study uses the SIFT algorithm as the feature point extraction algorithm for the accurate alignment of satellite images. The SIFT algorithm has strengths in extracting unique feature points that are invariant to various scales and rotations within the image, and it can extract corresponding point coordinates at the sub-pixel level, enabling more accurate bundle adjustment (Lowe, 2004; Velesaca et al., 2024).

The overall process of extracting correspondences involves calculating the similarity of each feature point after the feature points have been extracted, and then extracting the initial correspondences. At this point, the points with the top 30% similarity among the extracted correspondences are selected to be set as the initial matching set. As shown in Fig. 3, the RANSAC algorithm is applied to remove noise and outliers that may occur in the matching results, and ultimately only highconfidence correspondences are extracted (Fischler and Bolles, 1981).

Figure 3. Flowchart of the tie point extraction.

It is a method of selecting a random subset from the input matching results, repeatedly generating models based on this. The optimal matching set is then found by determining the number of matching points that correspond to the model. In this process, a homography transformation model was applied to estimate the model error. The final correspondences selected through this process are used as input data for the bundle adjustment stage.

2.2.2. Bundle Matrix Construction

The coordinate transformation between the KOMPSAT-3/3A satellite images and the ground is based on the RFM and the RPCs coefficients that constitute it. RFM is a model that can effectively describe the nonlinear relationship between pixel coordinates in satellite images and ground coordinates, allowing the conversion of the location information of extracted connection points within each image to ground coordinates. Eq. (3) represents the RFM-based bundle adjustment observation equation used in this study.

In Eq. (3), line and sample represent the row and column coordinates of corresponding points in the image space, and Δline and Δsample refer to the correction functions in the row and column directions.

Additionally, Line and Sample are the row and column coordinates calculated from the RPCs using the ground coordinates of the corresponding points, which are latitude, longitude, and height. In the proposed method, the initial RPCs correction function in the observation equation can be represented as a first-degree polynomial in the form of an affine transformation, as expressed in the following Eq. (4).

The coefficients a₀, a_s, a_l, b₀, b_s, and b_l in Eq. (4) are the coefficients of the correction function. It performs the role of correcting errors caused by various factors all at once.

Ultimately, the observation equation established from the correspondence point observation results to be used for bundle adjustment is expressed as Eq. (5). At this time, i represents the image index, j represents the feature point index, and k represents the ground coordinate index.

2.2.3. Block Bundle Adjustment

Bundle adjustment is a calibration method that simultaneously optimizes the 3D coordinates of multiple images and corresponding points to improve overall positional accuracy (Ban and Kim, 2024; Grodecki and Dial, 2003; Fu et al, 2019). In this study, bundle adjustment is performed by minimizing the reprojection error between the image coordinates of each corresponding point using RFM-based observation equations. The nonlinear observation equations are linearized by applying a first-order Taylor expansion, and the increments of the parameters are estimated in the form of least squares as shown in Eqs. (6) and 7.

The coefficients of the RPCs correction function and the ground projection coordinates of the tie points are selected as adjustment parameters. Each parameter is initially set to 0 and iteratively optimized through the bundle adjustment process. The weight matrix is dynamically recalculated and applied at each iteration step to reduce the error of observations during bundle adjustment, based on the adjusted results. And the covariance matrix of the residuals, thereby applying weights that reflect the reliability of each observation for more rigorous adjustments. The calculation of the covariance matrix for the new weights using the covariance matrix of the residuals is performed as shown in Eq. (8).

2.2.4. Result Image Generation

The correction parameters of each image are calculated through bundle adjustment and the ground coordinates of the corresponding points form the basis for generating the final result image. In this proposed method, a virtual DEM is constructed to generate result images tailored to a complex three-dimensional model rather than a simple plane (Fig. 4). For this construction, the inverse distance weighting (IDW) interpolation method is used to estimate the elevation values for all pixels of the resulting image to be generated based on the ground coordinates of the adjusted corresponding points. The 3D ground coordinates corresponding to each pixel are calculated, and the image coordinates corresponding to those pixels are retrieved from the original image using a backprojection method. Finally, the image resampling method is applied to generate the corrected result image.

Figure 4. The notion of result image generation using virtual DEM creation.

For the experiment, 29 KOMPSAT-3A L1R images of the Gwangju area taken between 2018 and 2022 were used. The K3A images were captured from 10 stripes, and the areas were selected to create overlapping regions between the stripes. A single CAS500-1 L2G image captured within the same area was collected to conduct the experiment. The reference image used was selected from the CAS500-1 that overlapped with as many KOMPSAT-3A images as possible. The L1R images used in the experiment were only subjected to radiometric and sensor corrections, without any geometric or ortho-rectification applied. Fig. 5 shows the locations where the 30 images used were taken, overlaid on each other. Table 1 shows the information for the 29 KOMPSAT-3A images used, and Table 2 shows the information for the one CAS500-1 image.

Figure 5. The 29 KOMPSAT-3A (black) images and 1 CAS500-1 (green) image we used.

Table 1 . Information on KOMPSAT-3A images used in the experiment.

No.	Orbit no.	Shooting date	Image center coordinate (Latitude, Longitude)	Image GSD (Column, Row)
1	17156	2018.05.04	34.82858029°, 126.82493711°	2.818, 2.696
2	17156	2018.05.04	34.94434339°, 126.79807505°	2.818, 2.696
3	17156	2018.05.04	35.05840244°, 126.77157367°	2.819, 2.697
4	17156	2018.05.04	35.17412625°, 126.74462852°	2.819, 2.697
5	17564	2018.05.31	34.92619426°, 126.69319410°	2.843, 2.493
6	17564	2018.05.31	35.04133984°, 126.66724353°	2.843, 2.493
7	17564	2018.05.31	35.15656655°, 126.64125727°	2.843, 2.493
8	19332	2018.09.25	35.01289916°, 126.93778748°	2.203, 2.200
9	19332	2018.09.25	35.12740934°, 126.90761118°	2.203, 2.200
10	22022	2019.03.22	34.86189533°, 126.48532307°	2.954, 2.531
11	22022	2019.03.22	34.97719495°, 126.45974778°	2.955, 2.532
12	22022	2019.03.22	35.09235398°, 126.43417530°	2.955, 2.532
13	22430	2019.04.18	34.85716336°, 126.63585741°	2.673, 2.416
14	22430	2019.04.18	34.97223208°, 126.60921339°	2.673, 2.415
15	22430	2019.04.18	35.08733444°, 126.58249726°	2.673, 2.416
16	26208	2019.12.24	34.88704441°, 126.74605892°	2.690, 2.930
17	26208	2019.12.24	35.00174973°, 126.71779278°	2.690, 2.930
18	26208	2019.12.24	35.11653393°, 126.68945303°	2.690, 2.930
19	27281	2020.03.04	34.85683172°, 126.87171237°	2.583, 2.971
20	27281	2020.03.04	34.97162974°, 126.84242082°	2.583, 2.970
21	38389	2022.03.09	34.93942529°, 126.84854691°	2.787, 2.829
22	38389	2022.03.09	35.05414486°, 126.82108828°	2.787, 2.830
23	38389	2022.03.09	35.16874358°, 126.79360806°	2.787, 2.830
24	38661	2022.03.27	34.87650292°, 126.88798877°	2.623, 2.996
25	38661	2022.03.27	34.99238290°, 126.85866555°	2.623, 2.996
26	38661	2022.03.27	35.10625536°, 126.82978093°	2.623, 2.995
27	38782	2022.04.04	34.85256082°, 126.91484286°	2.953, 2.633
28	38782	2022.04.04	34.96797849°, 126.88871778°	2.953, 2.633
29	38782	2022.0404	35.08277419°, 126.86268768°	2.954, 2.633

Table 2 . Information on CAS500-1 images used in the experiment.

No.	Orbit no.	Shooting date	Image center coordinate (Latitude, Longitude)	Image GSD (Column, Row)
1	10019	2023.01.10	35.03137775°, 126.77581116°	2.188, 2.099

4.1. Generation Simulated Images

In the process of performing strict bundle adjustment, a simulated image of the CAS500-1 with precise location information was generated to be used as the reference image. However, in the case of the generated simulated images, it is impossible for the NGII to obtain the original images due to the precise location information they possess. To verify the accuracy of the simulated image generation process and to enhance the reliability of the generated simulated images, we initially generated simulated images using KOMPSAT-3A images. Fig. 6 shows the comparison results between the L2R image generated using the KOPSAT-3A L2G image and the original L2R image. During the process of converting ortho images to oblique images, the size of the images decreases compared to the range of ground coordinates estimated during the initial orthorectification. As a result, blank areas that do not exist in the original image are created, but it can be confirmed that the non-blank areas are generated normally.

Figure 6. Overlay results of KOMPSAT-3A L2R original and generated simulated images.

Using the proposed method validated with KOMPSAT-3A images, a simulated CAS500-1 L2R image was generated. The simulated CAS500-1 L2R images can be generated by estimating the RPCs based on the L2G images and the metadata information of the images. These data are provided for free through the NGII. The L2G images used for the simulated images in this study and the generated L2R simulated image results are shown in Fig. 7. In the restored L2R image, it can be observed that the original ortho images has rotated similarly to its geometry at the time of capture.

Figure 7. Overlay results of CAS500-1 L2R original and generated simulated images.

4.2. Bundle Adjust Results

In this study, bundle adjustment was performed using the generated simulation CAS500-1 L2R image in Section 4.1 as a reference. Table 3 presents positional positional errors of the extracted tie point before and after bundle adjustment. Throughout the feature extraction process, matching was omitted for image pairs that shared the same orbit number. This approach aimed to reduce the number of excessively redundant initial feature points before performing the feature extraction processing.

Table 3 . Results of tie point extraction before and after bundle adjustment.

Dataset	No. of images	No. of pairs	Total number of feature points	Maximum distance error	Minimum distance error	Mean distance error
Before bundle adjustment	30	435	10068	311.9412	0.4263	163.0033
After bundle adjustment	30	435	10068	6.6136	0.0002	1.2611

Prior to the bundle adjustment experiment, the average positional error between tie points was approximately 163.0 pixels. However, after the bundle adjustment, this value was significantly reduced to about 1.3 pixels. Notably, the maximum positional error also decreased substantially, confirming the successful execution of the bundle adjustment process.

Fig. 8 shows the overlay results of the satellite images before and after applying the proposed method to 30 satellite images. Fig. 9 shows zoomed-in regions of the overlay results from Fig. 8, both before and after bundle adjustment. As observed in Fig. 9, the overlay of uncorrected satellite images appears blurry due to significant positional errors. The overlay results of the processed images demonstrate that the blurriness caused by positional errors has been effectively resolved. In the initial overlay results, distinct separations are observed such as rivers appearing as multiple entities. However, ghosting effects are eliminated and the overlay of single objects is accurately achieved in the images after bundle adjustment.

Figure 8. Results of 30 images before and after bundle adjustment.

Figure 9. Zoom-in results of 30 images before and after bundle adjustment.

Fig. 10 is a 3D map showing the results of a virtual DEM generated by applying interpolation to the estimated ground coordinates. The visualization reveals that high-altitude areas such as mountains or tall buildings appear significantly elevated compared to their surroundings. Terrain variations are observed even in flat areas like farmlands, with elevation changes occurring more gradually compared to nearby hills or mountainous regions. This phenomenon is presumed to result from errors introduced during the interpolation process using ground coordinates estimated through bundle adjustment.

Figure 10. Results of generated virtual DEM using interpolated estimated ground coordinates.

Fig. 11 shows the reprojection errors when the estimated ground coordinates of the tie points are reprojected onto the image. Blue points represent reprojection errors extracted before bundle adjustment and red points indicate the errors after bundle adjustment. The graph shows that all extracted tie points exhibit errors below 7 pixels, and the previously scattered distribution of distance errors converges to a consistent value. This confirms that the distance errors between tie points decreased after bundle adjustment, indicating that the model achieved stable estimation.

Figure 11. Results of reprojection error graph for ground coordinates.

This study proposes a method for estimating geometric alignment from L1R satellite images by combining the inverse georeferencing method with relative geometry correction. It performs bundle block adjustment using only tie points which does not require GCPs through the process. The proposed method successfully aligned a large number of satellite images through inverse georeferencing and bundle adjustment. The reprojection errors between the satellite images were reduced to less than 1.3 pixels, as quantitatively confirmed. We demonstrated the feasibility of geometric correction and ground coordinate estimation using multiple images without GCPs. This research has established the applying geometric correction to images without GCPs, or with an insufficient number of them. As a result, location information estimation using only images and metadata is possible for areas that are inaccessible due to physical or institutional reasons. Furthermore, we observed improved results based on reference images regarding the inaccurate absolute position accuracy of the existing relative geometric correction.

When the proposed algorithm is extended and applied, bundle block adjustment and geometric correction for heterogeneous (multi-type) satellite images can be performed. Since this experiment generates L1R-simulated images based on L2G images, KOMPSAT-3A images were also generated through inverse orthorectification to ensure consistency in the data used. Therefore, conducting experiments using original L1R images would allow for validation of the results obtained in this study.

In this paper, we observed improved bundle adjustment results when using simulated images created through inverse georeferencing as the reference image. We conducted experiments in urban areas to confirm the feasibility of bundle block adjustment using reference images for a large number of images. Given the significant improvement in accuracy, it seems feasible to conduct follow-up research targeting mountainous areas and North Korea, where GCPs is lacking. Quantitative analysis is needed for further research comparing the proposed method with existing bundle adjustment techniques. With the successful correction of L1R images, subsequent research on real-time bundle adjustment processing using collected satellite images is possible.

This work is supported by the Korea Agency for Infrastructure Technology Advancement (KAIA) grant funded by the Ministry of Land, Infrastructure and Transport (Grant No. RS-2022-00155763). Additionally, it was carried out with the support of the “Cooperative Research Program for Agriculture Science and Technology Development (Project No. PJ 016233)” supported by the Rural Development Administration, Republic of Korea.

No potential conflict of interest relevant to this article was reported.