Accuracy Evaluation on Geolocation of the Chinese First Polar Microsatellite (Ice Pathﬁnder) Imagery

: Ice Pathﬁnder (Code: BNU-1), launched on 12 September 2019, is the ﬁrst Chinese polar observation microsatellite. Its main payload is a wide-view camera with a ground resolution of 74 m at the subsatellite point and a scanning width of 744 km. BNU-1 takes into account the balance between spatial resolution and revisit frequency, providing observations with ﬁner spatial resolution than Terra/Aqua MODIS data and more frequent revisits than Landsat-8 OLI and Sentinel-2 MSI. It is a valuable supplement for polar observations. Geolocation is an essential step in satellite image processing. This study aims to geolocate BNU-1 images; this includes two steps. For the ﬁrst step, a geometric calibration model is applied to transform the image coordinates to geographic coordinates. The images calibrated by the geometric model are the Level1A (L1A) product. Due to the inaccuracy of satellite attitude and orbit parameters, the geometric calibration model also exhibits errors, resulting in geolocation errors in the BNU-1 L1A product. Then, a geometric correction method is applied as the second step to ﬁnd the control points (CPs) extracted from the BNU-1 L1A product and the corresponding MODIS images. These CPs are used to estimate and correct geolocation errors. The BNU-1 L1A product corrected by the geometric correction method is processed to the Level1B (L1B) product. Although the geometric correction method based on CPs has been widely used to correct the geolocation errors of visible remote sensing images, it is difﬁcult to extract enough CPs from polar images due to the high reﬂectance of snow and ice. In this study, the geometric correction employs an image division and an image enhancement method to extract more CPs from the BNU-1 L1A products. The results indicate that the number of CPs extracted by the division and image enhancements increases by about 30% to 182%. Twenty-eight images of Antarctica and ﬁfteen images of Arctic regions were evaluated to assess the performance of the geometric correction. The average geolocation error was reduced from 10 km to ~300 m. In general, this study presents the geolocation method, which could serve as a reference for the geolocation of other visible remote sensing images for polar observations.


Introduction
Visible remote sensing plays an important role in earth observations by providing super-width and high spatial resolution visual images. Along with its advantages, it has a wide range of applications in environmental surveying and mapping, disaster monitoring, resource investigation, vegetation monitoring, etc. [1][2][3][4][5]. In polar regions, visible remote sensing provides comprehensive observations of features on the earth's surface, and thus it is a supplement to limited field observations. With climate warming, dramatic changes have taken place in the polar regions where glaciers have retreated [6,7], ice flow has accelerated [8,9], and sea-ice has shrunk rapidly [10]. However, many of the rapid changes occurring in polar regions are difficult to monitor due to the trade-off between the temporal and spatial resolutions of existing satellite sensors (fine spatial resolution with a long revisit period; coarse resolution with a short revisit period) [4,11]. For example, the Moderate-Resolution Imaging Spectroradiometer (MODIS) sensors aboard the Terra/Aqua satellites can provide daily observations that facilitate the capture of rapid surface changes [4], but the coarse spatial resolution (250-1000 m) of MODIS sensors is often inadequate for monitoring the collapse of small glaciers or the disintegration of small icebergs. In contrast, the Landsat-8 OLI/Sentinel-2 MSI sensor has a higher spatial resolution (30 m/10 m) than MODIS, providing more details of the snow and ice surface changes. However, the 16-day/10-day revisit period of the Landsat-8 OLI/Sentinel-2 MSI sensor limits its application in the study of time-sensitive events, such as sea ice drift, which may evolve rapidly in a few days. Therefore, a sensor that provides high-resolution remote sensing data on a daily frequency or satellite constellations are needed for observing the rapid changes in polar regions.
Launched on 12 September 2019 and developed through the collaboration between Beijing Normal University, Sun Yat-sen University, led by Shenzhen Aerospace Dongfanghong HIT Satellite Ltd., Ice Pathfinder (Code: BNU-1) is the first Chinese polar-observing microsatellite. It is in a sun-synchronous orbit (SSO) with an altitude of 739 km above Earth's surface, a semi-major axis of 7,116,914.419 m, an inclination of 98.5238 Degrees, and an eccentricity of 0.000220908. Weighing only 16 kg, BNU-1 carries an optical payload with a panchromatic band and four multispectral bands. The spatial resolution at the sub-satellite point is approximately 74 m from the ground. The wide swath of BNU-1 (744 km) provides a 5-day revisit period of polar regions up to 85 • latitude. BNU-1 takes into account the balance between spatial resolution and revisit frequency, providing observations with finer spatial resolution than Terra/Aqua MODIS data and more revisit frequency than Landsat-8 OLI and Sentinel-2 MSI, benefiting the environmental monitoring of the polar regions. Also, the low cost of BNU-1 makes it financially feasible to construct a constellation observation system [12]. A five-satellite constellation system provides the ability to observe polar environmental elements on a daily basis with a spatial resolution finer than 100 m.
Image geolocation is an essential process prior to the application of satellites. However, geolocation errors are commonly found in visible images. For example, the images from MODIS have a geolocation error of 1.3 km in the along-track direction and 1.0 km in the along-scan direction without correction [13]. Geolocation errors need to be corrected because they cause uncertainty in satellite data and have a serious impact on the applications of satellite data for environmental monitoring [14]. The geolocation errors are usually corrected by parametric and non-parametric correction models [13,15,16]. Both these models correct the errors of a satellite image by matching the CPs obtained from the target image (the image with geolocation errors) and the corresponding points from the reference image with high geolocation accuracy [17]. The parametric correction model corrects the errors by optimizing the inner and external orientation parameters in the geometric calibration model based on the differences between the CPs from the target and reference images [13,15,18], while the non-parametric model is performed by establishing the coordinates transformation model between coordinates of the target image and coordinates of the reference images based on the CPs [8,11,12].
Both the parametric and non-parametric correction models are highly dependent on the amount of CPs [17,19,20]. Various methods have been used to increase the number of CPs extracted from the images, such as image division [21,22] and histogram equalization [22], and GCP sampling optimization [17]. Other methods have been used to eliminate the mis-matched CPs such as random sample consensus (RANSAC) [17,19], etc. These Both the parametric and non-parametric correction models are highly dependent on the amount of CPs [17,19,20]. Various methods have been used to increase the number of CPs extracted from the images, such as image division [21,22] and histogram equalization [22], and GCP sampling optimization [17]. Other methods have been used to eliminate the mis-matched CPs such as random sample consensus (RANSAC) [17,19], etc. These methods are commonly used for images with rich textures at low-and mid-latitudes. However, due to the high reflectance of ice and snow surfaces at high latitudes, the texture of images in polar regions is rarely observed. It is necessary to explore methods for correcting the geolocation errors of the images of polar regions.
This study aims to develop a geolocation method for polar images from BNU-1. The BNU-1 images for several regions of Antarctica and Greenland were used to demonstrate the effectiveness of this proposed geolocation method to deliminate the geolocation errors. This paper is organized as follows. Section 2 describes the data and the study area. Section 3 describes the geolocation method of the BNU-1 images in detail. Section 4 shows the performance of the geolocation method. The discussion of the results is shown in Section 5. Conclusions are given in Section 6.

Data
BNU-1 Imagery. BNU-1 has obtained more than 6000 images covering Antarctica and Greenland since it was launched. It provides the observations in panchromatic, blue, green, red, and red-edge spectral bands. Twenty-eight images of Antarctica and fifteen images of Greenland in the panchromatic band were collected for geolocation and accuracy evaluation. As shown in Figure 1, the images of Antarctica are distributed over the Amery Ice Shelf, Victoria Land, Dronning Maud Land, and Pine Island Glacier. The images of Greenland cover the west and north of Greenland. MODIS Imagery. MODIS is a key instrument onboard the Terra and Aqua satellites, which were launched on 18 December 1999, and 4 May 2002, respectively, providing global coverage every one to two days. Since the MODIS sensor has high geolocation accuracy (50 m for one standard deviation) [4] and a daily revisit capability, we used MODIS Imagery. MODIS is a key instrument onboard the Terra and Aqua satellites, which were launched on 18 December 1999, and 4 May 2002, respectively, providing global coverage every one to two days. Since the MODIS sensor has high geolocation accuracy (50 m for one standard deviation) [4] and a daily revisit capability, we used MOD02QKM and MYD02QKM products as the reference images for the geometric correction of the BNU-1 images in this study. The geolocation error of MOD02QKM (MYD02QKM) is 50 m or better, which is finer than the pixel size of the MODIS image [13,23,24]. It is reasonable to use MODIS images as the reference data in this study since the spatial resolution of BNU-1 images is 80 m.
Coastline dataset. We used the high-resolution vector polylines of the Antarctic coastline (7.4) [25] of 2021 from the British Antarctic Survey (BAS). We also used the MEaSUREs MODIS Mosaic of Greenland (MOG) 2005, 2010, and 2015 Image Maps, Version 2 [26] from the NASA National Snow and Ice Data Center (NSIDC) to obtain the Greenland coastline. We evaluated the geolocation error of the BNU-1 image visually by comparing the coastline dataset with the geolocation of the BNU-1 image.

Methods
There are two steps for geolocating the BNU-1 images in this study. The first step is geometric calibration. In this step, a geometric calibration model is constructed to transform the image coordinates to geographic coordinates. The images with geographic coordinates are the BNU-1 Level 1A (L1A) product. The second step is the geometric correction. The geolocation errors of the BNU-1 L1A product are corrected by an automated geometric correction method in this step. This method is designed to correct the geolocation errors of the images of polar regions where surface textures rarely exist. The corrected BNU-1 L1A product, which has high geolocation accuracy, is named the BNU-1 Level 1B (L1B) product.

Description of Geometric Calibration Model
A rigorous geometric calibration model was constructed for transforming the image coordinates to the geographic coordinates for the BNU-1 images. The timing, position, and altitude of satellites and camera position parameters are used as inputs of the model. The model is shown as follows [27]: where t is the scanning time of an imaginary line.
indicates the coordinates of the Global Positioning System (GPS) antenna phase center, which are measured by a GPS receiver on the satellite in the WGS84 coordinate system (derived from ECEF) at t. m is the scale factor determined by the orbital altitude. R WGS84 J2000 t , R J2000 body t and R body camera are the rotation matrix of the coordinate system from J2000 to WGS84 at t, the rotation matrix from the satellite's body-fixed coordinate system to J2000 coordinate system at t, and the rotation matrix from the camera coordinate system to the satellite's bodyfixed coordinate system, respectively. D x D y D z T is the coordinates of the GPS antenna phase center in the satellite's body-fixed coordinate system. d x d y d z T is the translations of the origin of the camera coordinate system relative to the satellite's bodyfixed coordinate system. tan ψ x tan ψ y 1 T represents the value of the coordinates of point (x, y) corresponding to the detector direction angle model composed of the camera's principal point, focal length, charge coupled device (CCD) installation position, and lens distortion. X Y Z T WGS84 represents the ground coordinates of the point (x, y) in the World Geodetic System 1984 (WGS84) coordinate system.
tan ψ x and tan ψ y describe the directional angle of point (x, y) in the camera coordinate system [27][28][29], and this can be calculated by Equation (2), where f is the focal length of the camera. tan This step is conducted on the Windows Server 2016 Standard operating system on the Intel(R) Xeon(R) Gold 5220R CPU @2.20 GHz, 256 GB RAM. It is a whole-day unattended automatic data production system.

Uncertainty Evaluation of Geometric Calibration Model
In addition to systematic errors, the geolocation of acquired images is also affected by random errors. The satellite imaging process is affected by various complex conditions such as attitude adjustment, attitude measurement accuracy, and imaging environment. In addition, as a microsatellite, BNU-1 s low cost and the imaging environment of polar regions limit, to a certain extent, the overall accuracy and stability of the measurement equipment of attitude and position. Moreover, due to the wide swath of BNU-1 s camera, the imaging time of its single-scene image is long. As the satellite attitude is adjusted along the imaging, random errors in attitude measurement cause random geolocation errors in single-scene images and multiple-scene images. Since the measurement error of GPS can be regarded as translation error (Equation (1)), which is equivalent to the satellite rotating at a small angle, we only designed and carried out an experiment to simulate the influence of the satellite's attitude angle change on the geolocation change through Equation 1. The satellite's attitude is determined by the roll angle, pitch angle, and yaw angle. We randomly selected an image and simulated the angles of roll, pitch, and yaw, which changed from 0 • to 0.5 • with a step size of 0.1 • , to obtain 216 (6 * 6 * 6) groups of geolocations. The quintic polynomial method was used to fit the scatter plot.

Automated Geometric Correction Processing Method
Since the space environment is complex and variable during satellite launch and operation [15,28,30,31], the geographic coordinates calculated by geometric calibration models with the pre-launch laboratory measurement parameters usually have geolocation errors of about several hundred meters to several kilometers [32,33]. In addition, the random error of the attitude measurement cannot be eliminated due to the lack of ground control points in polar regions. An automated geometric correction method based on CPs matching was developed to improve the geolocation accuracy of the BNU-1 L1A product. There are three steps involved in the method. Firstly, we selected the reference image with a high geolocation accuracy for the BNU-1 images. Then, the Scale Invariant Feature Transform algorithm (SIFT) [33] was used to extract the CPs from both the BNU-1 image and the corresponding reference image. Finally, geometric correction was conducted on the BNU-1 image based on the CPs. The flow chart is shown in Figure 2. Our experiment was conducted on the operating system of Windows 10 on the Intel(R) Core (TM) i5-5200 CPU @2.20 GHz, 8 GB RAM. We used the programming language of python2.7 to implement the one-stop processing of the automatic geometric correction. In this process, the programming language of MATLAB was used to realize SIFT algorithm, and the software of ArcGIS 10.6 was used to realize data preprocessing and geometric correction.

Step 1: Reference Images Selection
There were three criteria for selecting a reference image. Firstly, we selected the MODIS images with the same spatial coverage as the BNU-1 image as well as up to 5 h of different acquisition times. Secondly, the MODIS images covering the BNU-1 image and its surrounding 10 km area were chosen (one BNU-1 image corresponds to multiple MODIS images). Thirdly, the SIFT algorithm was adopted to extract the CPs from each image pair, where the BNU-1 image and the MODIS image are the target image and the reference image, respectively. The MODIS image with the most CPs was used as a reference image for the geometric correction. If more than one MODIS image has the highest number of CPs, the one whose acquisition time is closer to the BNU-1 image's acquisition time is preferred as the reference image. The reference image used for geometric correction of the BNU-1 image is referred to as the corresponding MODIS image hereinafter.

Step 2: Automatic CPs Extraction
The amount and spatial distribution of CPs are key factors for geometric correction because they have direct impacts on the geometric correction accuracy of the corrected images. In this study, we applied the SIFT algorithm based on MATLAB language to extract the CPs automatically from the BNU-1 L1A and the corresponding MODIS image. Due to the lack of texture features of snow and ice surfaces at high latitudes, CPs extracted from the original image pair are usually not sufficient for correcting the geolocation errors. When the number of CPs extracted from the image pair needs to be increased, image division and image enhancement methods are used to enhance the texture features of satellite images [1,21,22,34].
The combination of an image division method and an image enhancement method was applied to highlight surface features of the BNU-1 and MODIS images in this study. The extraction of CPs was carried out in three steps in Step 2 ( Figure 2). Firstly, we extracted the CPs from the original BNU1-1 image and the corresponding MODIS image by using the SIFT algorithm. The Euclidean distances between each pair of CPs from the image pair were calculated. To avoid mismatches of the points, we eliminated the largest 10% points in the Euclidean distance. Secondly, we extracted the CPs from the image pair after processing by different image division schemes. The paired images were divided into 2 × 2 = 4 (Scheme 1) and 3 × 3 = 9 (Scheme 2) sub-images [22]. Then, an adaptive piecewise linear enhancement consisting of three rules for low, middle, and high reflectance ranges was used to enhance each sub-image [34]. We extracted the CPs from all the pairs of the sub-images again by the SIFT algorithm. The largest 30% of the extracted CPs in Euclidean distance were eliminated in this step. Finally, the CPs extracted in the above two steps were merged and de-duplicated to obtain the CPs with the largest number, which were taken as the final CPs for the geometric correction.

Step 3: Geometric Correction
This study performed the geometric correction of the BNU-1 L1A product with the original spatial resolution (74 m) by using a quadratic polynomial (POLYORDER2) model in ArcGIS 10.6 software. We obtained the BNU-1 L1B product by resampling the geolocation error-corrected image to 80 m spatial resolution using the nearest resampling method.

Geolocation Accuracy Evaluation
To evaluate the geolocation accuracy of the BNU-1 L1A/L1B product, we re-extracted the CPs from the BNU-1 L1A/L1B product by the SIFT algorithm as the verification points and calculated the root mean squared error (RMSE) of the verification points: where n is the number of points participating in the accuracy evaluation, R xi and R yi represent the residuals of the i-th extracted points from the BNU-1 and the reference MODIS image in the X and Y coordinates, respectively.

Geolocation Accuracy of the BNU-1 L1A Product
The BNU-1 L1A images with 50% transparency are superimposed on the corresponding MODIS images in Figure 3. The sub-figures (a), (b), (c), and (d) correspond to the Sample Images A, B, C, and D shown in the red box in Figure 1. The prominent features in the images, such as coastlines, rocks, sea ice, etc., are blurred, indicating the mismatch of the geometric position between the BNU-1 images and the corresponding MODIS images. Obvious geolocation errors are observed in the BNU-1 L1A images. Table 1 shows the geolocation errors of the 42 scene BNU-1 L1A images. The errors of the BNU-1 L1A images range from 3 to 20 km, with an average of about 10 km. The geolocation errors of the sub-graphs in Figure 3a- where n is the number of points participating in the accuracy evaluation, xi R and yi R represent the residuals of the i-th extracted points from the BNU-1 and the reference MODIS image in the X and Y coordinates, respectively.

Geolocation Accuracy of the BNU-1 L1A Product
The BNU-1 L1A images with 50% transparency are superimposed on the corresponding MODIS images in Figure 3. The sub-figures (a), (b), (c), and (d) correspond to the Sample Images A, B, C, and D shown in the red box in Figure 1. The prominent features in the images, such as coastlines, rocks, sea ice, etc., are blurred, indicating the mismatch of the geometric position between the BNU-1 images and the corresponding MODIS images. Obvious geolocation errors are observed in the BNU-1 L1A images. Table 1 shows the geolocation errors of the 42 scene BNU-1 L1A images. The errors of the BNU-1 L1A images range from 3 to 20 km, with an average of about 10 km. The geolocation errors of the subgraphs in Figure 3a- Figure 4 shows the distributions of the geolocation errors of the CPs in the X and Y directions for each image shown in Figure 3. The length and direction of the vectors in Figure 4 represent the magnitude and direction of the CPs' geolocation errors. The directions and the magnitude of geolocation errors for each image are not consistent ( Figure 4). For example, the CPs' geolocation errors of Sample Image A, B, and D are less than 12 km, while most CPs' geolocation errors of Sample Image C are up to 19 km. And the geolocation errors in the middle part of Sample Image A are smaller than the errors in the edges of the image, while Sample Image B shows a quite different distribution of geolocation errors. In addition, the direction of the geolocation errors of the CPs shown in Sample Images B, C, and D are also different from the center to the periphery of the images. The CPs' geolocation errors within an image also vary significantly. For example, geolocation errors in Sample Image C are less than 2000 m in the center-west parts and more than 15,000 m in the east and southwest parts (Figure 4). The results illustrate that the distribution of the CPs' geolocation errors varies in each image and indicates that some local distortions exist in the BNU-1 L1A product.  Figure 5). Under imaging conditions in polar regions, the random error of attitude urement cannot be eliminated due to the lack of ground control points in polar r [28]. To obtain high-precision geolocation products, it is necessary to add CPs to the for geometric correction. Since there are few textures observed on high-reflectance snow surface in polar regions, the extraction of CPs is the key to geometric correct   Figure 5). Under imaging conditions in polar regions, the random error of attitude measurement cannot be eliminated due to the lack of ground control points in polar regions [28]. To obtain high-precision geolocation products, it is necessary to add CPs to the image for geometric correction. Since there are few textures observed on high-reflectance ice and snow surface in polar regions, the extraction of CPs is the key to geometric correction.

Influence of Image Division and Enhancement on the CPs Extraction
Sample Image E is a typical image for polar regions. Most of the features in the im are ice sheets and snow with limited boundary features, and only a few of them are ice with well-defined boundaries. However, due to the high reflectance of ice and sn surfaces in polar regions, the textures of the ice sheet and snow can rarely be observed images. The ice sheet area of Sample Image E is a good case for evaluating the effecti ness of various image enhancement methods for increasing the control points on the sheet. Five image enhancement methods, which are linear enhancement, piecewise lin enhancement, Gaussian enhancement, equalization enhancement, and square root hancement, were applied to enhance Sample Image E. Figure 6 shows the distribution CPs extracted from the images enhanced by different image enhancement methods. T numbers of CPs extracted from the original image and the image stretched by the f enhancement methods were 31, 32, 76, 37, 27, and 22, respectively. By comparing th five enhancement methods, we found that the piecewise linear enhancement meth makes the surface textures in the interior of the ice sheet more distinct, and as a result, most CPs were extracted from the image. Therefore, the piecewise linear enhancem (Figure 6c) is considered to be more suitable for enhancing the images of polar regions

Influence of Image Division and Enhancement on the CPs Extraction
Sample Image E is a typical image for polar regions. Most of the features in the image are ice sheets and snow with limited boundary features, and only a few of them are sea ice with well-defined boundaries. However, due to the high reflectance of ice and snow surfaces in polar regions, the textures of the ice sheet and snow can rarely be observed in images. The ice sheet area of Sample Image E is a good case for evaluating the effectiveness of various image enhancement methods for increasing the control points on the ice sheet. Five image enhancement methods, which are linear enhancement, piecewise linear enhancement, Gaussian enhancement, equalization enhancement, and square root enhancement, were applied to enhance Sample Image E. Figure 6 shows the distributions of CPs extracted from the images enhanced by different image enhancement methods. The numbers of CPs extracted from the original image and the image stretched by the five enhancement methods were 31, 32, 76, 37, 27, and 22, respectively. By comparing these five enhancement methods, we found that the piecewise linear enhancement method makes the surface textures in the interior of the ice sheet more distinct, and as a result, the most CPs were extracted from the image. Therefore, the piecewise linear enhancement (Figure 6c) is considered to be more suitable for enhancing the images of polar regions. Remote Sens. 2021, 13, x FOR PEER REVIEW 13 Sample Image E was also used to assess the influence of image division on CP traction. The image was divided into four sub-images (Scheme 1) and nine sub-im (Scheme 2) and then each sub-image was individually enhanced with the piecewise li stretching method. Figure 7 shows the distribution of CPs extracted from the origina age by Scheme 1, and by Scheme 2, respectively. The number of extracted points are 435, and 596, respectively. More CPs are extracted in the center and the lower right co of the image (the blue border area) divided by Scheme 2 (Figure 7c) compared to the inal image (Figure 7a) and the image divided by Scheme 1 (Figure 7b). As shown in T 2, the amount of CPs extracted from the image increases by 30% to 182% when the im division and piecewise linear enhancement were applied to the images. However, does not mean we can get more CPs if the image is divided into more sub-images. amount of CPs extracted from the Sample Image A, B, C, and D is less when Scheme applied to divide these images.  Sample Image E was also used to assess the influence of image division on CPs extraction. The image was divided into four sub-images (Scheme 1) and nine sub-images (Scheme 2) and then each sub-image was individually enhanced with the piecewise linear stretching method. Figure 7 shows the distribution of CPs extracted from the original image by Scheme 1, and by Scheme 2, respectively. The number of extracted points are 245, 435, and 596, respectively. More CPs are extracted in the center and the lower right corner of the image (the blue border area) divided by Scheme 2 (Figure 7c) compared to the original image ( Figure 7a) and the image divided by Scheme 1 (Figure 7b). As shown in Table 2, the amount of CPs extracted from the image increases by 30% to 182% when the image division and piecewise linear enhancement were applied to the images. However, this does not mean we can get more CPs if the image is divided into more sub-images. The amount of CPs extracted from the Sample Image A, B, C, and D is less when Scheme 2 is applied to divide these images.
The CPs extracted from Sample Image E were used to correct the geolocation errors of the image. The BNU L1A/L1B image with 50% transparency is superimposed on the corresponding MODIS image in Figure 8. There is a distinct displacement between the BNU-1 L1A image and the MODIS image, while the displacement between the BNU-1 L1B image and the MODIS image can barely be discerned. This result indicates that the geolocation correction method improves the geolocation accuracy of the BNU-1 L1Aproduct.
of the image (the blue border area) divided by Scheme 2 (Figure 7c) compared to the inal image (Figure 7a) and the image divided by Scheme 1 (Figure 7b). As shown in 2, the amount of CPs extracted from the image increases by 30% to 182% when the i division and piecewise linear enhancement were applied to the images. However does not mean we can get more CPs if the image is divided into more sub-images amount of CPs extracted from the Sample Image A, B, C, and D is less when Schem applied to divide these images.    The CPs extracted from Sample Image E were used to correct the geolocation errors of the image. The BNU L1A/L1B image with 50% transparency is superimposed on the corresponding MODIS image in Figure 8. There is a distinct displacement between the BNU-1 L1A image and the MODIS image, while the displacement between the BNU-1 L1B image and the MODIS image can barely be discerned. This result indicates that the geolocation correction method improves the geolocation accuracy of the BNU-1 L1A product. In addition to Sample Image E, Sample Image F was also selected to evaluate the effectiveness of the CPs extraction scheme proposed. Most of the areas of Sample Image F are covered by the ice sheet, and only a few areas are fjords. Figure 9 shows the image after geometric correction of the BNU-1 L1A image using the CPs extracted from the original image and the optimal control point extraction scheme (Scheme 2). The CPs extracted from the original image are few and unevenly distributed. If these points are directly used for geometric correction of the BNU-1 L1A image, the corrected image will be severely distorted (Figure 9a). More, and more evenly distributed CPs are extracted of the image divided by Scheme 2 (Figure 9b) compared to the original image (Figure 9a). The corrected BNU-1 image overlaps well with the MODIS image (Figure 9b). In addition to Sample Image E, Sample Image F was also selected to evaluate the effectiveness of the CPs extraction scheme proposed. Most of the areas of Sample Image F are covered by the ice sheet, and only a few areas are fjords. Figure 9 shows the image after geometric correction of the BNU-1 L1A image using the CPs extracted from the original image and the optimal control point extraction scheme (Scheme 2). The CPs extracted from the original image are few and unevenly distributed. If these points are directly used for geometric correction of the BNU-1 L1A image, the corrected image will be severely distorted (Figure 9a). More, and more evenly distributed CPs are extracted of the image divided by Scheme 2 (Figure 9b)

Geolocation Accuracy of the BNU-1 L1B Product
The verification points used to evaluate the geolocation errors of the BNU-1 L1B product for the Sample Images A, B, C, and D are shown in Figure 10a. We compared the coordinates of verification points in the BNU-1 L1A/L1B images with those in the corresponding MODIS images in Figure 10b, c. The verification points extracted from the BNU-1 L1A product (green dots) are distributed on one side of the 1:1 line (black diagonal line), which means that geolocation errors exist in the BNU-1 L1A product, while the verification points extracted from the BNU-1 L1B product (red dots) are almost scattered on the 1:1 line. We fitted the linear relationships between the coordinates of the verification points from BNU-1 L1A/BUN-1 L1B and the MODIS image. The regression coefficients, intercepts, and determined coefficients of the relationship fitted by the BNU-1 L1B product are significantly better than those fitted by BNU-1 L1A. The coordinates of the points from the BNU-1 L1B product show great consistency with the coordinates from the MODIS images. The geolocation accuracy of the BNU-1 L1B images was improved significantly. After geometric correction, the average geolocation error was reduced from 10,480.31 m to 301.14 m (Table 1).
We obtained the image mosaics of the Amery Ice Shelf and Victoria Land in Antarctica and northern Greenland in the panchromatic band of BNU-1 ( Figure 11). Mismatches in the coastlines were found in the image mosaics from the BNU-1 L1A product. However, the coastlines in the image mosaics from the BNU-1 L1B product are consistent with the existing coastline dataset [25]. Even though the junction of adjacent images in the image mosaic from the BNU-1 L1B product has greatly improved coherence compared to the BNU-1 L1A product, the BNU-1 L1B product still has an average geolocation error of ~300 m. For example, the discontinuous waters at the junction are found in Victoria Land (obvious mismatches in the yellow circle in Figure 11a) from the BNU-1 L1A product, while the mosaic from the BNU-1 L1B product has consistent waters at the junction regions (slight mismatches in the green circle in Figure 11a).

Geolocation Accuracy of the BNU-1 L1B Product
The verification points used to evaluate the geolocation errors of the BNU-1 L1B product for the Sample Images A, B, C, and D are shown in Figure 10a. We compared the coordinates of verification points in the BNU-1 L1A/L1B images with those in the corresponding MODIS images in Figure 10b, c. The verification points extracted from the BNU-1 L1A product (green dots) are distributed on one side of the 1:1 line (black diagonal line), which means that geolocation errors exist in the BNU-1 L1A product, while the verification points extracted from the BNU-1 L1B product (red dots) are almost scattered on the 1:1 line. We fitted the linear relationships between the coordinates of the verification points from BNU-1 L1A/BUN-1 L1B and the MODIS image. The regression coefficients, intercepts, and determined coefficients of the relationship fitted by the BNU-1 L1B product are significantly better than those fitted by BNU-1 L1A. The coordinates of the points from the BNU-1 L1B product show great consistency with the coordinates from the MODIS images. The geolocation accuracy of the BNU-1 L1B images was improved significantly. After geometric correction, the average geolocation error was reduced from 10,480.31 m to 301.14 m (Table 1).
We obtained the image mosaics of the Amery Ice Shelf and Victoria Land in Antarctica and northern Greenland in the panchromatic band of BNU-1 ( Figure 11). Mismatches in the coastlines were found in the image mosaics from the BNU-1 L1A product. However, the coastlines in the image mosaics from the BNU-1 L1B product are consistent with the existing coastline dataset [25]. Even though the junction of adjacent images in the image mosaic from the BNU-1 L1B product has greatly improved coherence compared to the BNU-1 L1A product, the BNU-1 L1B product still has an average geolocation error of 300 m. For example, the discontinuous waters at the junction are found in Victoria Land (obvious mismatches in the yellow circle in Figure 11a) from the BNU-1 L1A product, while the mosaic from the BNU-1 L1B product has consistent waters at the junction regions (slight mismatches in the green circle in Figure 11a).   The two small graphs in each sub-graph are the enlarged version of the orange and blue rectangular areas on the black diagonal line.

Discussion
Although microsatellites have the advantages of compactness, low cost, and flexibility, their flight attitude may be unstable occasionally and the equipment measuring the attitude and position may be inaccurate or perform poorly, which leads to large geolocation errors in the images [35]. For example, the geolocation errors in UNIFORM-1 s visible images are 50-100 km [35]. The BNU-1 L1A product has smaller geolocation errors, but the average geolocation error can still be up to 10 km, which is close to that of the Luojia 1-01 data [30].
The non-parametric geometric correction methods are widely applied to geolocation correction of images without distinguishing the error sources [16,31], such as HJ-1A/B CCD images [36] and Unmanned Aerial Vehicle (UAV) images [37]. The geolocation accuracy of these geometric correction methods relies mainly on adequate CPs. However, the limited texture features of the ice and snow surface in polar images make it difficult to extract the CPs. Some studies prove that image division and image enhancement have the ability to increase the amount of extracted CPs [1,21,22,34]. However, images used by these previous studies are from low-and mid-latitudes and contain rich land surface features. The correction for polar images with few texture features is rarely documented. This study proposes the geometric correction method to reduce the geolocation errors of the visible images for polar regions. The results indicate that piecewise linear enhancement highlights more surface features of ice and snow surfaces than other image enhancement methods. Some other studies have also proved that piecewise linear enhancement is effective in highlighting more texture features of the ice and snow surfaces in polar images [34]. More CPs can be observed after the image pair is processed by image division and piecewise linear enhancement. Different division schemes can be adopted to obtain more CPs for different image pairs.
In addition, we compared the geolocation accuracy of Sample Image A-F after the correction through the CPs extracted from the original image and the geolocation accuracy after correction through the optimal CPs extraction scheme (Table 3). It was found that the geolocation accuracy of Sample Image A-E was not significantly improved. The geolocation accuracy of Sample Image A-E after geometric correction based on the CPs extracted from the original image was close to the level of 250 m (the pixel size of MODIS image), and it was difficult to further improve by adding CPs on this basis. However, for some images where the ice sheet is widely distributed, such as Sample Image F, the proposed method effectively prevents the distortion of the corrected image caused by the lack of CPs on the ice sheet by adding CPs. The increase in CPs can remarkably improve the geolocation accuracy of such images. Therefore, the automatic geometric correction method proposed in this study is of great significance for the correction of images in polar regions with rare feature points. Although the method presented in this study has some advantages in correcting the geolocation errors of polar images, it also has its limitations. For example, only two division schemes were applied in the BNU-1 images, and the division fractions for different images are still worth further study. Besides, the geolocation accuracy of some sub-images with relatively uniform surface features is difficult to improve by using the proposed method. This indicates that the SIFT algorithm has its limitation for finding more CPs. Therefore, it is necessary to explore some other CPs extraction methods for increasing the number of CPs. Since deep learning methods have been widely used in image registration [38,39], it is worth exploring the possible application of deep learning methods on CPs extraction from the images of polar regions.

Conclusions
In this study, we present the geolocation method for BNU-1 images including two steps. For the first step, a rigorous geometric calibration model was applied to transform the image coordinates to the geographic coordinates for the BNU-1 images. The images geolocated by the geometric calibration model are the BNU-1 L1A product. For the second step, an automated geometric correction method was used to reduce the geolocation errors of the BNU-1 L1A product. The images corrected by the geometric correction method are the BNU-1 L1B product.
The geometric correction method is commonly used for improving the geolocation errors of the visible image. However, the texture features of the ice and snow surfaces are rarely seen in polar images, which makes it difficult to find the CPs. The combination of the image division method and piecewise linear image enhancement method was applied to the BNU-1 L1A product and the corresponding MODIS images, and the results indicate that the CPs extracted increased by 30% to 182%, which can effectively improve the geometric accuracy of the BNU-1 images.
The geolocation method was applied to 28 images of Antarctica and 15 images of Arctic regions. The average geolocation error was reduced from 10 km to~300 m. The coastlines in the image mosaics from the BNU-1 L1B product were consistent with the coastline dataset. These results suggest that the geolocation method has the ability to improve the geolocation errors of BNU-1 images and other satellite images in polar regions.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.