An Emergency Georeferencing Framework for GF-4 Imagery Based on GCP Prediction and Dynamic RPC Refinement

GaoFen-4 (GF-4) imagery has very potential in terms of emergency response due to its gazing mode. However, only poor geometric accuracy can be obtained using the rational polynomial coefficient (RPC) parameters provided, making ground control points (GCPs) necessary for emergency response. However, selecting GCPs is traditionally time-consuming, labor-intensive, and not fully reliable. This is mainly due to the facts that (1) manual GCP selection is time-consuming and cumbersome because of too many human interventions, especially for the first few GCPs; (2) typically, GF-4 gives planar array imagery acquired at rather large tilt angles, and the distortion introduces problems in image matching; (3) reference data will not always be available, especially under emergency circumstances. This paper provides a novel emergency georeferencing framework for GF-4 Level 1 imagery. The key feature is GCP prediction based on dynamic RPC refinement, which is able to predict even the first GCP and the prediction will be dynamically refined as the selection goes on. This is done by two techniques: (1) GCP prediction using RPC parameters and (2) dynamic RPC refinement using as few as only one GCP. Besides, online map services are also adopted to automatically provide reference data. Experimental results show that (1) GCP predictions improve using dynamic RPC refinement; (2) GCP selection becomes more efficient with GCP prediction; (3) the integration of online map services constitutes a good example for emergency response.


Introduction
GF-4 (GF is a Chinese abbreviation for high-resolution) is China's first geostationary earth-observation satellite [1,2] enabled with gazing mode, when continuous imagery series of the same area can be captured at intervals of seconds, providing an ideal observation approach for emergency response when facing disasters such as flood.The ground sample distance (GSD) of its panchromatic image is 50 m.However, only poor geometric accuracies (about 4000 m) can be expected from the rational polynomial coefficient (RPC) provided.Even with a GSD of 50 m, its geometric accuracy will not meet the ideal geolocation error requirement of products; for example, Mouillot [3] carried out an analysis of user needs for burned area (BA) products and pointed out that on average an ideal geolocation accuracy requirement would be 1100 m, with the best GSD (500 m) of the discussed dataset being much lower than GF-4.This geometric accuracy is far from satisfying if emergency decisions are to be made based on the provided coordinates, especially when higher geometric accuracies are required, making ground control points (GCPs) necessary to achieve better accuracies.problem.The main objective is to provide a fast and reliable GCP selection method for the urgent georeferencing of GF-4 imagery in emergency response.The key feature is dynamic GCP prediction, which is able to predict even the first GCP and the prediction is refined as the selection goes on.This is being done by two techniques: (1) GCP prediction using RPC parameters and ( 2) dynamic RPC refinement.Additionally, a third technique, the integration of online map services, is adopted to automatically provide reference data, which will save a great deal of time as long as there is Internet access.
This paper is organized as follows: Section 2 gives detailed information about how to predict and refine RPC when selecting GCPs and the core principles of online map services; experimental verification is presented in Section 3; a discussion is provided and a brief conclusion drawn in the last two sections, respectively.

Materials and Methods
The novel emergency georeferencing framework for GF-4 imagery proposed is mainly featured by GCP prediction, dynamic RPC refinement, and online map service integration.The overall framework is shown in Figure 1 with the three key steps colored in green.
Remote Sens. 2017, 9, 1053 3 of 17 In this paper, a novel emergency georeferencing framework for GF-4 imagery based on RPC refinement and GCP selection assisted by integrated online map services is proposed to solve this problem.The main objective is to provide a fast and reliable GCP selection method for the urgent georeferencing of GF-4 imagery in emergency response.The key feature is dynamic GCP prediction, which is able to predict even the first GCP and the prediction is refined as the selection goes on.This is being done by two techniques: (1) GCP prediction using RPC parameters and (2) dynamic RPC refinement.Additionally, a third technique, the integration of online map services, is adopted to automatically provide reference data, which will save a great deal of time as long as there is Internet access.
This paper is organized as follows: Section 2 gives detailed information about how to predict and refine RPC when selecting GCPs and the core principles of online map services; experimental verification is presented in Section 3; a discussion is provided and a brief conclusion drawn in the last two sections, respectively.

Materials and Methods
The novel emergency georeferencing framework for GF-4 imagery proposed is mainly featured by GCP prediction, dynamic RPC refinement, and online map service integration.The overall framework is shown in Figure 1 with the three key steps colored in green.
Figure 1.This shows the overall framework.The ground control points (GCPs) can then be used to carry out georeferencing using tools such as ArcGIS or ENVI, which is rather mature and thus not covered in this paper.
The first part of this section will focus on GCP prediction and dynamic RPC refinement.Some core principles of online map services will be discussed in the second part of this section.

GCP Prediction and RPC Refinement
Major time-consuming issues when carrying out GCP selection include reference data collecting, GSD difference between reference data and the data to be georeferenced, and GCP locating (find the exact corresponding pixel from another image is rather difficult and tedious without prediction).Thus, the proposed GF-4 emergency georeferencing framework introduces three features to improve efficiency: online map services have made image gathering unnecessary, and the remaining two issues will be settled in this section by zoom level estimation, GCP prediction, and dynamic RPC refinement.

GCP Prediction
Traditional GCP selection is extremely tedious and slow mainly because prediction is not available, which is the fact for traditional GCP selection operations using tools such as ENVI or ArcGIS.
Figure 1.This shows the overall framework.The ground control points (GCPs) can then be used to carry out georeferencing using tools such as ArcGIS or ENVI, which is rather mature and thus not covered in this paper.
The first part of this section will focus on GCP prediction and dynamic RPC refinement.Some core principles of online map services will be discussed in the second part of this section.

GCP Prediction and RPC Refinement
Major time-consuming issues when carrying out GCP selection include reference data collecting, GSD difference between reference data and the data to be georeferenced, and GCP locating (find the exact corresponding pixel from another image is rather difficult and tedious without prediction).Thus, the proposed GF-4 emergency georeferencing framework introduces three features to improve efficiency: online map services have made image gathering unnecessary, and the remaining two issues will be settled in this section by zoom level estimation, GCP prediction, and dynamic RPC refinement.

GCP Prediction
Traditional GCP selection is extremely tedious and slow mainly because prediction is not available, which is the fact for traditional GCP selection operations using tools such as ENVI or ArcGIS.
However, Level 1 data (whose radiometric and systematic distortions have been removed but is not georeferenced) with RPC parameters are able to provide rough geographic coordinates, which can be used for prediction.To predict a GCP, the geographical coordinate of the selected point must be derived using its pixel coordinate (row number and column number) or the other way around.This is carried out using the RPC model, a detailed description of which is given by Zhang [14].Since the method used has been rather mature, there is no need to go into detail in this paper.
With the derived geographical coordinate (typically as latitude and longitude), we will be able to navigate reference data to the predicted GCP or the other way around (see Figure 2).With this feature, GCP selection will be well facilitated.However, we can make more precise predictions as the selection goes on.
However, Level 1 data (whose radiometric and systematic distortions have been removed but is not georeferenced) with RPC parameters are able to provide rough geographic coordinates, which can be used for prediction.To predict a GCP, the geographical coordinate of the selected point must be derived using its pixel coordinate (row number and column number) or the other way around.This is carried out using the RPC model, a detailed description of which is given by Zhang [14].Since the method used has been rather mature, there is no need to go into detail in this paper.
With the derived geographical coordinate (typically as latitude and longitude), we will be able to navigate reference data to the predicted GCP or the other way around (see Figure 2).With this feature, GCP selection will be well facilitated.However, we can make more precise predictions as the selection goes on.
Figure 2.This shows the process of GCP prediction.A potential feature/point is selected from the GF-4 Level 1 image, and the reference data is then navigated to the predicted location (the junction point of the red cross located at the point of the green pushpin in the reference image).

Dynamic RPC Refinement
The main idea of dynamic RPC refinement is to refine RPC parameters as soon as a new GCP is confirmed with the purpose of carrying out a more precise prediction for the next GCPs.
However, there are dozens of parameters included in an RPC parameter set, while in the very first place there is only one GCP available to refine the dozens of parameters in the beginning and only a few more available later.One way to solve this problem is to simulate more observations from one or a few more GCPs.
To do this, an RPC refinement method similar to those used by Zhang [14] as well as Jinshan & Xiuxiao [15] is adopted.This is a terrain-independent RPC model, which means terrain elevation will not be needed to derive geographical coordinates from pixel coordinates.The core idea is to construct a controlling grid consisting of a certain number of layers, each of which is laid at a different height above/below the same terrain area.In this way, the grid will cover a reasonable 3D space range to drive new RPC parameters, so the elevation variations of the ground area covered will not be able to effect the refinement.Once we refine the coordinates of all the junction points in this controlling grid with as few as one GCP, we have enough observations to refine the RPC parameter set.
A sketch map of the controlling grid is shown in Figure 3.
Figure 2.This shows the process of GCP prediction.A potential feature/point is selected from the GF-4 Level 1 image, and the reference data is then navigated to the predicted location (the junction point of the red cross located at the point of the green pushpin in the reference image).

Dynamic RPC Refinement
The main idea of dynamic RPC refinement is to refine RPC parameters as soon as a new GCP is confirmed with the purpose of carrying out a more precise prediction for the next GCPs.
However, there are dozens of parameters included in an RPC parameter set, while in the very first place there is only one GCP available to refine the dozens of parameters in the beginning and only a few more available later.One way to solve this problem is to simulate more observations from one or a few more GCPs.
To do this, an RPC refinement method similar to those used by Zhang [14] as well as Jinshan & Xiuxiao [15] is adopted.This is a terrain-independent RPC model, which means terrain elevation will not be needed to derive geographical coordinates from pixel coordinates.The core idea is to construct a controlling grid consisting of a certain number of layers, each of which is laid at a different height above/below the same terrain area.In this way, the grid will cover a reasonable 3D space range to drive new RPC parameters, so the elevation variations of the ground area covered will not be able to effect the refinement.Once we refine the coordinates of all the junction points in this controlling grid with as few as one GCP, we have enough observations to refine the RPC parameter set.
A sketch map of the controlling grid is shown in Figure 3.In this paper, the image space is equally divided into a 25 × 25 chessboard.Ten layers of such chessboards are constructed to cover an elevation range from −5000 to 5000 m.The elevation range is manually set to this range not because local elevation varies as such, but because (1) globally the elevation variation of common areas on earth should fall at a relatively "middle" position of this range and (2) the interval between controlling grid layers should be large enough to help keep a stable calculation.
After the controlling grid is set, an image space adjustment model is derived using the confirmed GCPs to fit the distribution of the pixel coordinate adjustments for all the GCPs.The model used differs according to the number of GCPs: a translation model for one GCP, a conformal transformation for two GCPs, and a 1st order polynomial transformation model for three or more GCPs.The notion "adjustment" here is used to represent the difference between the original pixel coordinate of a certain GCP and the pixel coordinate derived through an RPC model from its manually confirmed geographical coordinate.
Basically, the equation used to derive the pixel coordinates from the geographic coordinate of the junction point is shown in Equation (1), where (P, L, H) is the normalized geographic coordinate derived from (Lat, Lng, height) of a junction point and (rn, cn) is the normalized pixel coordinate derived from row and column (r, c).All the coefficients needed to carry out the normalization mentioned as well as the construction of the rational polynomials, i.e., NumL, DenL, NumS and DenS, come from the RPC parameters as is shown in Equation (2).In Equation ( 2), all the coefficients come from the RPC parameters, and (P,L,H) is derived by Equation (3) of a certain junction point.LAT_OFF, LAT_SCALE, LONG_OFF, LONG_SCALE, HEIGHT_OFF, HEIGHT_SCALE are also contained in the RPC parameters used.
(P, L, H) (P, L, H) (P, L, H) (P, L, H) Figure 3.This shows the structure of the controlling grid.Note that the image columns and rows are not necessarily parallel to meridians or the equator and the layers are not necessarily divided into 5 × 4 chessboards as showed.
In this paper, the image space is equally divided into a 25 × 25 chessboards.Ten layers of such chessboards are constructed to cover an elevation range from −5000 to 5000 m.The elevation range is manually set to this range not because local elevation varies as such, but because (1) globally the elevation variation of common areas on earth should fall at a relatively "middle" position of this range and (2) the interval between controlling grid layers should be large enough to help keep a stable calculation.
After the controlling grid is set, an image space adjustment model is derived using the confirmed GCPs to fit the distribution of the pixel coordinate adjustments for all the GCPs.The model used differs according to the number of GCPs: a translation model for one GCP, a conformal transformation for two GCPs, and a 1st order polynomial transformation model for three or more GCPs.The notion "adjustment" here is used to represent the difference between the original pixel coordinate of a certain GCP and the pixel coordinate derived through an RPC model from its manually confirmed geographical coordinate.
Basically, the equation used to derive the pixel coordinates from the geographic coordinate of the junction point is shown in Equation ( 1), where (P, L, H) is the normalized geographic coordinate derived from (Lat, Lng, height) of a junction point and (r n , c n ) is the normalized pixel coordinate derived from row and column (r, c).All the coefficients needed to carry out the normalization mentioned as well as the construction of the rational polynomials, i.e., NumL, DenL, NumS and DenS, come from the RPC parameters as is shown in Equation ( 2).In Equation ( 2), all the coefficients come from the RPC parameters, and (P,L,H) is derived by Equation ( 3 Remote Sens. 2017, 9, 1053 The image space adjustment is derived as follows: For each GCP, its pixel coordinate (r 0 , c 0 ) is logged when it is selected, and the geographical coordinate (Lat 0 , Lng 0 ) can be derived using the RPC parameters.The GCP is then manually confirmed to obtain a new geographic coordinate (Lat 1 , Lng 1 ), which can be used to derive a new pixel coordinate (r 1 , c 1 ) using the same set of RPC parameters.In this way, an image space adjustment (v r , v c ) can be calculated for this GCP as is shown in Equation ( 4).
The geographical coordinates of the junction points can then be refined using the image space adjustment model.This is done in three steps: (1) derive the pixel coordinate of each junction point using the original RPC parameters and its geographical coordinate; (2) derive a new pixel coordinate for each GCP by applying the adjustment model to its pixel coordinate; (3) derive a new geographical coordinate for each junction point using its adjusted pixel coordinate and the original RPC parameters.
With the refined geographical coordinates of the junction points, refined RPC parameters can be derived using least square estimation.To derive new RPC parameters, Equation ( 1) is transformed into Equation ( 5) so that the error equation can be derived as Equation (6), where V is the residual matrix, B the partial derivative matrix, x the new RPC parameter matrix, and l the observation matrix.Using the least square technique, x can be derived.
A brief description of the refinement method is given as follows: (1) Create a controlling grid along the image plane and imagine that there are several copies of this grid at different elevation levels.(2) Compute pixel coordinates offset coefficients using selected GCPs and derive the image space adjustment model.(3) Calculate the latitude/longitude of all the junction points of the grid at all heights using both RPC parameters and the image space adjustment model.(4) Derive a new set of RPC parameters using the latitude/longitude of all the junction points.

Online Map Services
With GCP prediction and dynamic RPC refinement, we will be able to design an emergency georeferencing framework for GF-4 imagery once the reference data can be easily obtained.The most ideal way to do this is to integrate online map services as reference data.
The basic hypothesis behind online map services acting as a substitution for local reference data is that approximate geometric accuracy can be provided by the former.For public online map services, an accuracy evaluation [6] was carried out across China against online map services including Google, Bing, and Mapbox, showing that all results are better than 5 m.These accuracies will not satisfy needs for fundamental GIS/mapping products, but will be helpful for departments and small corporations for whom high-accuracy reference data are not available or under urgent circumstances where low-accuracy GIS/mapping products can also work, as in the GF-4 imagery case.Meanwhile, fundamental GIS/mapping production will also find online map services meaningful, since high-accuracy reference data can be utilized to make a customized online map service.

Technical Framework
The basic function of online map services is to provide georeferenced map data to end-users through network.Original data can be stored and transferred as well as clips (tiles) derived from original data abide predefined rules.The tiling strategy is adopted by most of current online map service standards for convenience and low cost.
The tiling strategy defines a set of zoom levels to present global data in different scales, and the global data will be resampled to specified GSD and clipped into tiles of the same size within each zoom level based on a predefined latitude/longitude grid.As a user pans through the map, tiles will be requested according to the zoom level and coordinate where the user is viewing and transferred to the app, where tiles will be organized back into the right part of a global map.The principle of the online map services is shown in Figure 4.
According to Han [11], online map services can be divided into three layers from a functional view as is shown in Figure 5: the data layer, which provides different kinds of data; the service layer, which mainly covers transferring issues and consists of standards and specifications; and finally the application layer, which is in charge of data requesting, displaying, and other user-interface operations.These three layers can be summarized as data, transmission and application, respectively.

Standards
Organizations such as OGC (Open Geospatial Consortium) and OSGeo (Open Source Geospatial Foundation) have been developing series of standards and specifications to reduce possible developing obscurities or chaos concerning online map services, including but not limited to Web Map Service (WMS), Tile Map Service (TMS), and Web Map Tile Service (WMTS).

Standards
Organizations such as OGC (Open Geospatial Consortium) and OSGeo (Open Source Geospatial Foundation) have been developing series of standards and specifications to reduce possible developing obscurities or chaos concerning online map services, including but not limited to Web Map Service (WMS), Tile Map Service (TMS), and Web Map Tile Service (WMTS).

Standards
Organizations such as OGC (Open Geospatial Consortium) and OSGeo (Open Source Geospatial Foundation) have been developing series of standards and specifications to reduce possible developing obscurities or chaos concerning online map services, including but not limited to Web Map Service (WMS), Tile Map Service (TMS), and Web Map Tile Service (WMTS).
WMS was developed by OGC and first published in 1999, with version 1.3.0(ISO 19128) [16] issued on 15 March 2006.Standards, terms and definitions, basic service elements, and interfaces are covered in the specifications.WMS is characterized mainly by the following: (1) both text and images can be returned through WMS interfaces; (2) images of any reasonable range can be returned by the WMS server; (3) two basic operations, GetCapabilities and GetMap, are provided, with a third optional one, GetFeatureInfo.The transportation of WMS is not based on tiles of a fixed size, which means the data clipping can only be carried out when a request is posted, making the data transfer rather complex.
To make data preparing and transferring more efficient, tiling strategies have been introduced to new standards including Web Map Service Tile Caching (WMS-C by OGS) and TMS (by OSGeo).The most important character of the new standards is that original data is clipped into fixed-sized tiles according to geographic coordinates and pre-defined zoom levels.Since data clipping is no longer needed when responding data requests and tiles can even be cached on user-ends, data transferring becomes more efficient.
However, WMS-C itself becomes too complicated to support an old WMS and a new tiling strategy at the same time.For example, coding by tile index is not supported in WMS, which makes the tiling strategy of WMS-C rather clumsy.On the other side, although TMS is based on a tiling strategy in the very beginning, it is not an official standard.As a result, WMTS was proposed by OGC to make better use of the tiling strategy and to adopt new features such as the GetFeatureInfo interface to provide descriptive information about the data/feature to be transferred.

Coordinates Computation
Before actually predicting GCPs, a proper zoom level of the online map services should be estimated.The purpose is to find out the best zoom level of the online map service where the GSD is closest to the Level 1 data.First, the relationship between the zoom level and GSD must be made clear as is shown in Equation (7), where the zoom level is denoted by n, the average radius of Earth is denoted by R earth , and latitude is denoted by φ.An obvious deduction is that GSD doubles as the zoom level is decreased by 1.
Traverse Equation (7) gives Equation ( 8) to derive the closest zoom level n to the given GSD at given latitude φ.
Note that these equations may also differ for different standards.
The exact relationship between geographic coordinates and zoom level and tile index should also be given to display tiles correctly, which is shown as Equation ( 9), where x is the column index, y the row index, λ the longitude, φ the latitude, and n the zoom level ([0, maxLevel]).
The adverse process is to covert the column and row number of a pixel in a certain tile to its corresponding geographic coordinate, which is most important if we were to select GCPs from online map services.This computation can be carried out through the following steps: (1) Derive global pixel coordinate (x p , y p ) of the current pixel from its position in the tile and tile index (x, y), as is shown in Equation (10), where 256 is the tile size.
(2) Derive lat/lon from (x p , y p ) as is shown in Equation (11) where the total global column number and row number of pixels at the current zoom level are denoted by x sz and y sz , respectively.
Note that different standards may use different tile sizes, i.e., 256 × 256 or 512 × 512, and that the start point of y may be the top-left corner or the bottom-left corner-we only need to carry out minor adjustments to the given conversions.
With the estimated zoom level and coordinate, the GCP prediction process (Figure 2) should actually be Figure 6.
(1) Derive global pixel coordinate (xp, yp) of the current pixel from its position in the tile and tile index (x, y), as is shown in Equation (10), where 256 is the tile size.
(2) Derive lat/lon from (xp, yp) as is shown in Equation (11) where the total global column number and row number of pixels at the current zoom level are denoted by xsz and ysz, respectively.
Note that different standards may use different tile sizes, i.e., 256 × 256 or 512 × 512, and that the start point of y may be the top-left corner or the bottom-left corner-we only need to carry out minor adjustments to the given conversions.
With the estimated zoom level and coordinate, the GCP prediction process (Figure 2) should actually be Figure 6.As is shown in Figure 6, when a point is selected from the Level 1 data, the row and column of the point is logged.Together with the RPC parameters, an RPC computation can be carried out to derive the geographical coordinate of the selected point.On the other hand, with the GSD of the Level 1 imagery (when the imagery source is specified, a rough GSD is known), zoom level estimation can be carried out to evaluate the best zoom level of the online map service used.With the derived geographical coordinate and zoom level, we will be able to navigate the online map service to the predicted position at the best zoom level.

Results
To validate the proposed framework, a typical scenario was chosen GCPs were selected for a GF-4 Level 1 PAN image.These islands are located in Philippines with an elevation variation from 0 to less than 3000 m.This is a typical area for the proposed framework not because of its elevation variation but because it is very hard for operators to identify features as potential GCPs since (1) the Figure 6.This shows the actual process of GCP prediction.GSD is basically constant when the data source is known.For GF-4 PAN imagery, the GSD is about 50 m.
As is shown in Figure 6, when a point is selected from the Level 1 data, the row and column of the point is logged.Together with the RPC parameters, an RPC computation can be carried out to derive the geographical coordinate of the selected point.On the other hand, with the GSD of the Level 1 imagery (when the imagery source is specified, a rough GSD is known), zoom level estimation can be carried out to evaluate the best zoom level of the online map service used.With the derived geographical coordinate and zoom level, we will be able to navigate the online map service to the predicted position at the best zoom level.

Results
To validate the proposed framework, a typical scenario was chosen GCPs were selected for a GF-4 Level 1 PAN image.These islands are located in Philippines with an elevation variation from 0 to less than 3000 m.This is a typical area for the proposed framework not because of its elevation variation but because it is very hard for operators to identify features as potential GCPs since (1) the islands only occupy a small portion of the area covered, (2) the islands are surrounded by oceans, and (3) parts of them are covered by clouds in the image.
First of all, the Level 1 data was georeferenced using only its original RPC parameters to evaluate its geometry accuracy, using checkpoints selected from the online map service.GCPs were then selected using the proposed framework.After that, the original image was georeferenced using the selected GCPs.At last, the accuracies of the two georeferenced results were compared.
The online map service used was a hybrid map provided by Bing for both checkpoints' selection and GCP selection.The ground range of the GF-4 image is shown in Table 1.An overview of the range shows that most of the area is covered by ocean, with a smaller part being islands, as is shown in Figure 7.
Remote Sens. 2017, 9, 1053 11 of 17 islands only occupy a small portion of the area covered, (2) the islands are surrounded by oceans, and (3) parts of them are covered by clouds in the image.First of all, the Level 1 data was georeferenced using only its original RPC parameters to evaluate its geometry accuracy, using checkpoints selected from the online map service.GCPs were then selected using the proposed framework.After that, the original image was georeferenced using the selected GCPs.At last, the accuracies of the two georeferenced results were compared.
The online map service used was a hybrid map provided by Bing for both checkpoints' selection and GCP selection.The ground range of the GF-4 image is shown in Table 1.An overview of the range shows that most of the area is covered by ocean, with a smaller part being islands, as is shown in Figure 7.

GCP Prediction
The GCP prediction result for the first GCP is shown in Figure 8, where the right part is GF-4 Level 1 data (wrap image), with a detail window and global window under the main window, and on the left, an online map service (reference image).Within a rather large ground range (we can say a worldwide range), the GCP was successfully predicted and the zoom level estimated.Note that

GCP Prediction
The GCP prediction result for the first GCP is shown in Figure 8, where the right part is GF-4 Level 1 data (wrap image), with a detail window and global window under the main window, and on the left, an online map service (reference image).Within a rather large ground range (we can say a worldwide range), the was successfully predicted and the zoom level estimated.Note that this was the first GCP, meaning there were no other GCPs available for other possible GCP prediction methods.It is obvious that the framework proposed was able to give a rather good estimation of the zoom level as well as the location even for the first GCP.
Remote Sens. 2017, 9, 1053 12 of 17 this was the first GCP, meaning there were no other GCPs available for other possible GCP prediction methods.It is obvious that the framework proposed was able to give a rather good estimation of the zoom level as well as the location even for the first GCP.
Figure 8.This is a GCP prediction result.In the main window of the wrap image (top-right window), the selected point is denoted by a green cross mark, while the predicted point is denoted by a green pushpin in the main window of the online map service (top-left window).The predicted GCP was located in the same small island where the Level 1 image was clicked (see the two adjacent detail windows).

Dynamic RPC Refinement
Detail windows of the wrap image and reference image show a rather large discrepancy (actually, the distance is about 4 km) between the predicted location and its true location, which is shown in Figure 9.By prediction discrepancy, we mean the distance between the predicted point and its true location.To better show the prediction discrepancy, the predicted location (the red cross in the left view) was also manually drawn in the right view (the yellow cross), thus the distance between the yellow cross and the middle of the right view (the center of the big red cross) would be the prediction discrepancy, which is represented by a green line.Since dynamic RPC refinement was adopted, prediction discrepancy does not remain the same throughout the selecting process.To validate the dynamic RPC refinement, all discrepancies for each GCP selected were logged, as is shown in Table 2.  .This is a GCP prediction result.In the main window of the wrap image (top-right window), the selected point is denoted by a green cross mark, while the predicted point is denoted by a green pushpin in the main window of the online map service (top-left window).The predicted GCP was located in the same small island where the Level 1 image was clicked (see the two adjacent detail windows).

Dynamic RPC Refinement
Detail windows of the wrap image and reference image show a rather large discrepancy (actually, the distance is about 4 km) between the predicted location and its true location, which is shown in Figure 9.By prediction discrepancy, we mean the distance between the predicted point and its true location.To better show the prediction discrepancy, the predicted location (the red cross in the left view) was also manually drawn in the right view (the yellow cross), thus the distance between the yellow cross and the middle of the right view (the center of the big red cross) would be the prediction discrepancy, which is represented by a green line.
Remote Sens. 2017, 9, 1053 12 of 17 this was the first GCP, meaning there were no other GCPs available for other possible GCP prediction methods.It is obvious that the framework proposed was able to give a rather good estimation of the zoom level as well as the location even for the first GCP.
Figure 8.This is a GCP prediction result.In the main window of the wrap image (top-right window), the selected point is denoted by a green cross mark, while the predicted point is denoted by a green pushpin in the main window of the online map service (top-left window).The predicted GCP was located in the same small island where the Level 1 image was clicked (see the two adjacent detail windows).

Dynamic RPC Refinement
Detail windows of the wrap image and reference image show a rather large discrepancy (actually, the distance is about 4 km) between the predicted location and its true location, which is shown in Figure 9.By prediction discrepancy, we mean the distance between the predicted point and its true location.To better show the prediction discrepancy, the predicted location (the red cross in the left view) was also manually drawn in the right view (the yellow cross), thus the distance between the yellow cross and the middle of the right view (the center of the big red cross) would be the prediction discrepancy, which is represented by a green line.Since dynamic RPC refinement was adopted, prediction discrepancy does not remain the same throughout the selecting process.To validate the dynamic RPC refinement, all discrepancies for each GCP selected were logged, as is shown in Table 2.  Since dynamic RPC refinement was adopted, prediction discrepancy does not remain the same throughout the selecting process.To validate the dynamic RPC refinement, all discrepancies for each GCP selected were logged, as is shown in Table 2.The trend line of the prediction discrepancy is shown in Figure 10.The trend line of the prediction discrepancy is shown in Figure 10.The prediction discrepancy was 4829 m for the first GCP and decreased to 1994 m for the second GCP since the former GCP was available for RPC refinement.When RPC was refined using the confirmed two GCPs, the prediction discrepancy for the 3rd GCP turned to 84 m, and remained within [40,276] for the following GCPs, with most of them smaller than 100 m (about 2 pixels).This means, given GCPs, dynamic RPC refinement will help carry out stable and precise predictions.
Prediction discrepancies became stable after only two GCPs (except for the minor ascending trends from three to five GCPs and from 8 to 9 GCPs) because the image space adjustment model adopted changes as the number of GCPs changed, as described in Section 2.1.2.Given three or more, the type of the model used will be fixed, and new GCPs only bring updates to the coefficients of the model.In fact, the difference in the model used for two GCPs and that used for three GCPs is so small that two GCPs will suffice to stabilize a prediction.Thus, with one GCP, it is difficult to simulate a higher order distribution of prediction discrepancy, and the prediction discrepancy will not remain stable for the first two GCPs.In fact, if the second GCP falls near the 1st GCP, the prediction discrepancy will be almost 0 m.In contrast, if it falls far from the 1st GCP, the discrepancy will increase.With three or more GCPs available for RPC refinement, a 1st order polynomial transformation model can be established to better fit the global distribution of the offsets, which means the following prediction discrepancies will be relatively stable compared with the previous ones.
However, there was a minor ascending trend from three to five GCPs.The reason for the trend is that, when there were three or more GCPs, a 1st order polynomial transformation model was adopted; however, three or four GCPs were not enough to cover the whole image area.Therefore, when a new GCP was selected in an uncovered area, the discrepancy would likely increase since the model was subjected to a rather large update.This also explains the slight increase when the 9th GCP was selected, which was also located in an uncovered area.In conclusion, the limited number of GCPs and their distribution are responsible for the ascending trend from three to five GCPs.
After that, six or more GCPs were enough to maintain a good distribution along the entire image area, and the model was finally stabilized, which shows that, with as few as six GCPs, the proposed framework would be able to give steady and precise GCP predictions.
The combination of GCP prediction and dynamic RPC refinement also saved us a great deal of time.However, it was also not appropriate to carry out a quantitative comparison since the experience of the operator matters here.However, we were still able to draw a qualitative conclusion that these two features did help a great deal, according to our experience of selecting GCPs.The prediction discrepancy was 4829 m for the first GCP and decreased to 1994 m for the second GCP since the former GCP was available for RPC refinement.When RPC was refined using the confirmed two GCPs, the prediction discrepancy for the 3rd GCP turned to 84 m, and remained within [40,276] for the following GCPs, with most of them smaller than 100 m (about 2 pixels).This means, given GCPs, dynamic RPC refinement will help carry out stable and precise predictions.
Prediction discrepancies became stable after only two GCPs (except for the minor ascending trends from three to five GCPs and from 8 to 9 GCPs) because the image space adjustment model adopted changes as the number of GCPs changed, as described in Section 2.1.2.Given three or more, the type of the model used will be fixed, and new GCPs only bring updates to the coefficients of the model.In fact, the difference in the model used for two GCPs and that used for three GCPs is so small that two GCPs will suffice to stabilize a prediction.Thus, with one GCP, it is difficult to simulate a higher order distribution of prediction discrepancy, and the prediction discrepancy will not remain stable for the first two GCPs.In fact, if the second GCP falls near the 1st GCP, the prediction discrepancy will be almost 0 m.In contrast, if it falls far from the 1st GCP, the discrepancy will increase.With three or more GCPs available for RPC refinement, a 1st order polynomial transformation model can be established to better fit the global distribution of the offsets, which means the following prediction discrepancies will be relatively stable compared with the previous ones.
However, there was a minor ascending trend from three to five GCPs.The reason for the trend is that, when there were three or more GCPs, a 1st order polynomial transformation model was adopted; however, three or four GCPs were not enough to cover the whole image area.Therefore, when a new GCP was selected in an uncovered area, the discrepancy would likely increase since the model was subjected to a rather large update.This also explains the slight increase when the 9th GCP was selected, which was also located in an uncovered area.In conclusion, the limited number of GCPs and their distribution are responsible for the ascending trend from three to five GCPs.
After that, six or more GCPs were enough to maintain a good distribution along the entire image area, and the model was finally stabilized, which shows that, with as few as six GCPs, the proposed framework would be able to give steady and precise GCP predictions.
The combination of GCP prediction and dynamic RPC refinement also saved us a great deal of time.However, it was also not appropriate to carry out a quantitative comparison since the experience of the operator matters here.However, we were still able to draw a qualitative conclusion that these two features did help a great deal, according to our experience of selecting GCPs.

Georeferencing Accuracy
Eleven GCPs were selected with an RMSE of 12.3 pixels.The 11 GCPs are listed in Table 3, with Figure 11 showing the distribution.With the selected GCPs, we were able to georeference the original image using both the RPC and GCPs to get a Level 2B product.After that, the georeferencing accuracy was checked using 10 checking points, with an RMSE of 86.7 m.The original georeferencing accuracy was also checked, with the RMSE being 3894.4 m.The checking points are listed in Table 4 and the distribution is shown in Figure 12.Note that the GCPs and checking points were manually selected in two independent processes: GCP selection and precision checking.Compared with the original RPC parameters, a new set of RPC parameters given by the proposed framework decreased the RMSE of the georeferencing from 3894.4 (about 78 pixels) to 86.7 m (less than 2 pixels).This was a magnificent improvement and would definitely be quite convenient if it were under emergency response, not to mention the time saved by using GCP prediction, dynamic RPC refinement, and online map services.Compared with the original RPC parameters, a new set of RPC parameters given by the proposed framework decreased the RMSE of the georeferencing from 3894.4 (about 78 pixels) to 86.7 m (less than 2 pixels).This was a magnificent improvement and would definitely be quite convenient if it were under emergency response, not to mention the time saved by using GCP prediction, dynamic RPC refinement, and online map services.
Additionally, the RMSE of the selected GCPs was 12.3 pixels (see Table 3) and about 615 m.However, the accuracy of the final Level 2B product was improved to better than 100 m, which is about 2 pixels.This means that the instability of man-selecting can be compensated by the dynamic RPC refinement to some extent.

Discussion
As has been demonstrated in Section 3, the proposed framework is able to provide GCP prediction throughout the entire selection process, which will greatly improve efficiency.We actually did not carry out any quantitative evaluations about how much time was saved either by using online maps services as reference data or by using the combination of GCP prediction and dynamic RPC refinement to assist GCP selection.However, our experience indicates that, normally, at least 10-20 min will be needed to select GCPs for a GF-4 imagery, while, with the help of the proposed framework, the work can be done within 1 min.Note that there is no need to collect reference images anymore.In fact, the cost of collecting reference images cannot be simply measured by the time spent, the money spent, nor the personnel involved, since, in practice, a failed or even delayed collection of reference data may lead to unimaginable loss.It is also unreasonable to log how long it takes to select GCPs without GCP prediction, because the speed of GCP selection, for the same operator, varies from one operator to another, from one image to another, and even from one time to another.However, we do know that, when using an online map service as a reference image, only Internet access is needed, thus saving time, money, and personnel and even reducing the potential loss of property and human lives.
Aside from the proposed framework, there are also other ways to carry out GCP prediction and to refine the predictions.For example, geographical coordinates can be directly derived from pixel coordinates using a 2nd order polynomial transformation determined by the confirmed GCPs, just like our group has been using for years.However, the major problem lies in the first three GCPs that are needed to determine the transformation: this method is not capable of predicting the first three GCPs.Note that, during manual GCP selection, this is exactly when GCP prediction is most needed, since the first few GCPs are hardest to select, especially when there is a rotation between Level 1 imagery and the reference image or when the area covered is too large.
We may go further based on the proposed framework towards a full ability of automation by introducing an image matching technique.Since no matter how precisely the GCP is predicted, its final confirmation depends on human intervention.In fact, the reliability of the image matching technique will increase if we are able to carry out a good estimation of the potential matching area.This is precisely what we have done here.We carried out trials using SIFT and least squares matching (LSM): once the GCP was predicted, we clipped the image area around the selected point as well as the image area around the predicted point and tried to match them using the mentioned methods.The results were not satisfying, but this still remains a promising direction in which to proceed.
The proposed framework can also be applied to other Level 1 satellite imagery, and its application is not limited to emergency response as long as there are RPC parameters available.In addition, orthoimagery and digital line graphics (DLGs) with low geometric accuracy will also find this framework helpful, with a minor change to the GCP prediction method and the dynamic RPC refinement if needed.

Conclusions
An emergency georeferencing framework for GF-4 imagery based on RPC refinement and GCPs selection assisted by integrated online map services is proposed in the paper.The experiment has demonstrated that this framework will introduce an obvious improvement to the efficiency of GF-4 emergency georeferencing, which is often time-consuming and tedious because of the reference data collecting and lacking of the coordinate prediction.This will be especially meaningful in emergency response when reference data is hard to access, image matching is no longer reliable, and time is limited.
In addition, online map services have been demonstrated to be ideal providers of reference data, and their integration constitutes a good example of remote sensing and photogrammetry applications in emergency response.Online map services, once used as reference data, can be useful in other potential applications, such as the urgent georeferencing of unmanned aerial vehicle (UAV) images based on feature matching.

Figure 3 .
Figure 3.This shows the structure of the controlling grid.Note that the image columns and rows are not necessarily parallel to meridians or the equator and the layers are not necessarily divided into 5*4 chessboards as showed.

Figure 4 .
Figure 4.This shows the principle of online map services.Note that some standards do not tile the original data.

Figure 5 .
Figure 5.This shows the functional layers of online map services.

Figure 4 .Figure 4 .
Figure 4.This shows the principle of online map services.Note that some standards do not tile the original data.

Figure 5 .
Figure 5.This shows the functional layers of online map services.

Figure 5 .
Figure 5.This shows the functional layers of online map services.

Figure 6 .
Figure 6.This shows the actual process of GCP prediction.GSD is basically constant when the data source is known.For GF-4 PAN imagery, the GSD is about 50 m.

Figure 7 .
Figure 7.An overview of the experimental area.Most of the area is covered by ocean.

Figure 7 .
Figure 7.An overview of the experimental area.Most of the area is covered by ocean.

Figure 9 .
Figure 9.This shows the GCP prediction discrepancy.

Figure 8
Figure8.This is a GCP prediction result.In the main window of the wrap image (top-right window), the selected point is denoted by a green cross mark, while the predicted point is denoted by a green pushpin in the main window of the online map service (top-left window).The predicted GCP was located in the same small island where the Level 1 image was clicked (see the two adjacent detail windows).

Figure 9 .
Figure 9.This shows the GCP prediction discrepancy.

Figure 9 .
Figure 9.This shows the GCP prediction discrepancy.

Figure 10 .
Figure 10.This shows the trend line of the prediction discrepancies.The discrepancies (Y-axis) are measured in meters.

Figure 10 .
Figure 10.This the trend line of the prediction discrepancies.The discrepancies (Y-axis) are measured in meters.

Figure 11 .
Figure 11.This shows the distribution of the selected GCPs.

Figure 11 .
Figure 11.This shows the distribution of the selected GCPs.

Figure 12 .
Figure 12.This shows the distribution of the checking points.

Figure 12 .
Figure 12.This shows the distribution of the checking points.

Table 1 .
This table shows the ground range of the GF-4 image.

Table 1 .
This table shows the ground range of the GF-4 image.

Table 2 .
This table shows all the prediction discrepancies logged.The discrepancies are measured in pixels as well as in meters.

Table 2 .
This table shows all the prediction discrepancies logged.The discrepancies are measured in pixels as well as in meters.

Table 2 .
This table shows all the prediction discrepancies logged.The are measured in pixels as well as in meters.

Table 3 .
This table shows the list of the selected GCPs.

Table 3 .
This table shows the list of the selected GCPs.

Table 4 .
This table shows the checking points list.

Table 4 .
This table shows the checking points list.