Geometric Accuracy Assessment and Correction of Imagery from Chinese Earth Observation Satellites (hj-1 A/b, Cbers-02c and Zy-3)

The Chinese satellites HJ-1 A/B, CBERS-02C and ZY-3 have been recently launched and are considered as the main space platforms on orbit to acquire optical images for monitoring the Earth for various applications in China. The commercially distributed products (Level 1 or 2) of those satellites usually lack sufficient information (about platform, sensor and ephemeris) that is the key to geometrically correct the acquired images. It is therefore always a challenging issue and the first step to assess the geometric accuracy, which is a key part of qualities in spatial data, of the images from those satellites before generation of geometrically accurate image products. This paper first describes an operational methodology to assess the geometric accuracy of those satellite images. The methodology automatically collects dense and spatially well distributed ground control points (GCP) against reference imagery and then fits those GCPs to the given geometric math model. The geometric accuracy of an image can then be assessed from the overall fitness of those GCPs and their distribution of geometric errors along and across track. The residual mean square (RMS) parameter is used to indicate the degree of overall fitness of the GCPs to the photogrammetric system. The distribution of geometric errors may be random or approximated by a second or higher order polynomial functions; the latter case is generally considered as a systematic error that was not removed completely in the Level 1 or 2 data product. In order to draw solid conclusions, a significant number of samples are selected for each of those satellites by taking variations of landscapes into consideration. The assessment experiments demonstrate that the accuracy of HJ-1 A/B is often very poor, that of CBERS-02C is better than the situation of HJ-1 A/B but records poor accuracy for most samples, and that of ZY-3 is the best among all satellites under investigation and has few samples with poor accuracy. According to the assessment results, this paper suggests an operational correction methodology to improve the accuracy for those satellites, particularly for the HJ-1 A/B and CBERS-02C. Operational production proves that the proposed correction methodology is capable of achieving much higher accuracy than traditional ones and the achieved accuracy meets high standard product requirements for such applications as mapping.


INTRODUCTION
China has recently launched a large number of Earth observation (EO) satellites from the late 2000's in response to the strong demand of timely spatial information over the territories of China from both public and private sectors (Guo, 2012) (Li & Townshend, 2013).HJ-1, CBERS-02C1 and ZY-3 are three major platforms among those EO satellites on orbit and offering commercially available imagery (CRESDA, 2014).
The constellation of HJ-1 consists of two small optical satellites A and B (launched in September 2008), and one small SAR satellite C (launched in November 2012), dedicated to monitor environment and mitigate disasters; CBERS is an international cooperation space program between Chinese and Brazilian governments that have launched a series of satellites, and CBERS-02C (launched in December 2011) is the most recently launched one; ZY-3 (launched in January 2012) is the first three-line CCD array satellite with high resolution for stereo mapping in China.See the web source (CRESDA, 2014) and literatures (Wang, et al., 2010) (Li, 2012) for detail of the parameters of all these Chinese satellites.
The agencies, including the China Centre for Resources Satellite Data and Application -CRESDA, of these satellites have made their image data commercial.However, the distributed data does not include sufficient information about platform, sensor and ephemeris to build an accurate geometric relationship between an image and the ground, which makes it difficult to use the imagery in applications when geometric accuracy is highly concerned.
Assessment of the quality of satellite image data, particularly the geometric accuracy, is always the fundamental step toward finding solutions and then processing the data.Many literatures (Wang, et al., 2014) (Jiang, et al., 2014) (Zhang, et al., 2014) (Tang, et al., 2013) (Jiang, et al., 2013) (Xu, et al., 2012) focus on data quality analysis in Level 0 products from the three satellite platforms, which are not commercially available; very few literatures dedicate their effort in assessing the Level 1 products from the three satellite platforms or the literatures if available limit their work in the products from a single satellite platform.Particularly, there are a few papers (Pan, et al., 2013) (Liu, et al., 2012) for assessing ZY-3 Level 1 product data quality, a single paper (Hu & Tang, 2011) for HJ-1 A/B Level 1 product data quality, and none for CBERS-02C product at all.Besides making contribution to assessing the Level 1 products of HJ-1 A, HJ-1 B, CBERS-02C and ZY-3, this paper will use an operational methodology to assess the data quality from the four satellites.The operational assessment methodology is able to simplify the data analysis procedure and be more efficient in the time-consuming analysis work.
The operational methodology proposed in this paper for Level 1 data quality assessment includes the following steps: 1. Automatically collect ground control points (GCP) from orthorectified reference imagery; 2. Fit the GCPs to a given photogrammetric math model; 3. Analyse the residual errors.
The phase correlation technology (Kuglin & Hines, 1975) is deployed here to perform automatic GCP collection because of its high reliability (high tolerance for noise and illumination change) and sub-pixel accuracy in matching (Stone, et al., 2001) (Hoge, 2003).The given photogrammetric math model is either a simulated physical model, i.e., the 3D parametric model (Toutin, 2003), or a rational polynomial function if the distributed Level 1 data includes its own rational polynomial coefficients (RPC), which will be referred as to the RPC model hereafter.The residual errors are analysed to evaluate the fitness of a given math model with the collected GCPs.
The data sets for testing the operational assessment methodology are selected carefully with key attention on minimizing or at least lowering the probability of wrong matches when collecting GCPs and representing variable landscapes if possible.As such, the image should be cloud/snow/ice free or almost free and is covered mostly by land since most of water and cloud/snow/ice areas generally have rare texture information and may lead to wrongly matched GCPs.
The detailed analysis demonstrates that the HJ-1 A/B Level 1 data product generally has a systematic distortion in its image internally that cannot be removed properly by regular approaches; the CBERS-02C Level 1 data product often has lens distortion in its cameras that should be calibrated on board and/or when processing Level 0 data into Level 1 one; the ZY-3 Level 1 data product has the best quality overall and almost free of any visually noticeable distortion from lens or other physical reasons.
Having the above analysis, the paper proposes an approach in addition to the operational assessment methodology to correct those images, particularly the HJ-1 A/B and CBERS-02C data, to meet the accuracy criteria set for high-accuracy-demanding applications (e.g., mapping).The operational approach uses the collected GCPs which are spatially evenly distributed and dense to compute a new RPC model.The computed RPC model is better than the finite element approach (Goshtasby, 1988) (Hu & Tang, 2011) because it takes elevation into consideration (Toutin, 2004).The same data sets for the assessment experiments are processed by the new approach to generate orthorectified products.Manual verification, assuming error free when up to sub-pixel accuracy is concerned in applications, is conducted over the orthorectified products against reference imagery and illustrates that the proposed correction methodology can improve the accuracy significantly in comparison with a traditional one (e.g., using the simulated 3D parametric model (Toutin, 2003) or the original RPC).

Photogrammetry Math Models
A photogrammetry math model establishes the relationship between a 2-D image space and a 3-D ground system for an image acquired by a satellite (Slama, et al., 1983).What degree the model fits to the relationship determines the final accuracy of a corrected image.From the processing point of view and regarding how to model the distortions in satellite imagery, Toutin (Toutin, 2004) provides a thorough review of various popular photogrammetric math models, among which the 3D parametric model (referred as to the Toutin model hereafter) and the RPC model have more advantages than the other models and are often applied in production systems.

Toutin Model
Toutin model is a set of equations that are rigidly derived from the classic collinearity equations.See (Toutin, 2003) for the equations.This math model simplifies the collinearity equations for a push broom image, which has many scan lines and each scan line corresponds to one projection center, and hence is implementable for production.Each parameter in the equations represents one or a combination of physical realities of the full viewing geometry.It has been proved to be robust for various optical and radar data (Toutin, 2003) (Toutin, 2004).To achieve pixel-level accuracy, the Toutin model requires in theory a certain number of GCPs.

RPC Model
For the purpose of maintaining confidentiality of the parameters of satellites, the RPC model has been gaining more and more popularity in the society of satellite data processing and tends to become dominant for most commercial satellites launched recently.
The RPC model uses two rational functions (each is a ratio of two polynomial functions) to approximate the rigorous geometric relationship when an image is acquired.Because the rational polynomial coefficients represent no physical realities explicitly, the clients who receive the data do not know what the physical parameters of a satellite are.The two rational functions are: where  ,  are image coordinates,  ,  and  are ground coordinates, and  ,  ,  ,    ( = 0,1, … ,19) are the coefficients for the four polynomial functions.
To remove bias and compensate high order variations, an affine transformation function is often used as an addition to the rational functions (Fraser, et al., 2001) (Dial & Grodecki, 2005) (Hu & Tao, 2004).It is proved (Fraser, et al., 2001) that the additional polynomial functions can very well remove bias and compensate various variations due to platform drifting etc.The compensated RPC model is where   (, )    (, ) are a first order polynomial function of variables  and .
One distinguishing advantage in the RPC model is its generality for implementation regardless of difference among satellites.

GCP Collection Technologies
GCPs are indispensable to achieve accurate orthorectified products (referred as orthos later on) and assess the accuracy of the orthos.There are many ways to collect GCPs that can be categorized into manual collection and automatic collection.
Manual GCP collection may require either field trip using GPS tools or human visual collection against reference data (assisted by computer or not) or both.Automatic GCP collection is generally referred to automatically matching an image (raw) to another (reference), which is a critical step for automating orthos generation.
There are many algorithms developed for image matching (Zitova & Flusser, 2003) (Wyawahare, et al., 2009).The normalized cross correlation (NCC) method (Gonzalez & Woods, 2008) is popularly applied among them.It requires less computation and can be easily implemented.This algorithm is however very sensitive to image noise and illumination change in an image.For orthorectification of satellite imagery, this sensitivity usually forces users to choose the reference imagery that are close enough to the target ones regarding their satellite platforms, radiometric characteristics, weather conditions as well as acquisition date and time.
As an image matching algorithm in the frequency domain, the phase correlation algorithm (Kuglin & Hines, 1975) can overcome the above limitations.It can even generate sub-pixel accuracy (Stone, et al., 2001).Since this algorithm is reliable due to its much less sensitivity to image noise and illumination change, the successful rate of matching is extremely high and the automatically collected GCPs make the following operational assessment and correction methodology applicable.

Operational Accuracy Assessment Methodology
Many researches were dedicated to assess the geometric accuracy of satellite images (Toutin, 2004).The most commonly used approach for assessment is to use ground control points (Toutin, 2004) (Jacobsen & Passini, 2003) (Zhou & Li, 2000), either manually collected over reference data or being facilitated by GPS tools.Those GCPs are split into two sets: the first set (literally named as GCPs) is used to correct the image and the second set (named as independent check points -ICPs) used to verify the accuracy after correction.This popular approach combines the correction process into the quality assessment procedure.Other approaches include the figure condition (Sertel, et al., 2007) etc.
This paper proposes an operational methodology for geometric accuracy assessment.This methodology spins off the correction process and is designed with three major steps: automatic collection of GCPs, fitting GCPs to a math model, and analysis of residual errors of all GCPs.

Automatic Collection of GCPs
Orthorectified imagery is used as reference.The phase correlation algorithm is deployed for matching a target (raw) satellite image to the reference imagery for collecting GCPs.Since the chosen image matching algorithm has high tolerance for illumination change and low signal/noise ratio over images, more flexibility is allowed to select the reference imagery from various sources, particularly available orthos from the satellites other than the one to be assessed.In order to ensure the collected GCPs are evenly distributed over the horizontal direction, all candidates of GCPs are evenly gridded across the target image; Meanwhile, in order to increase the possibility of distributing GCPs evenly over the vertical direction, the gridded GCP candidates are densified enough to cover the most range of elevations in the ground area covered by the image.

Fitting GCPs to a Math Model
The commercial products (Level 1 or 2) from the four satellites do not include critical information regarding the instant positions of platforms and orientations of their cameras when scanning images, it is therefore impossible to establish an accurate geometric model between the ground and the image.When a product data set comes with a RPC file, a RPC math model will be built accordingly; when there is no RPC file available, a Toutin model will be simulated by using defaulted values for its unknown parameters, which need to be resolved by using GCPs.All GCPs will then be applied to the built model.The built model is adjusted first by the given GCPs (Toutin, 2004) (Hu & Tao, 2004).The adjusted model is then used to calculate the residual errors of all GCPs that represent how well the built math models the relationship between the ground and the image.If most GCPs do not fit well (i.e., with large residual errors), it can be concluded that either the math model is inappropriate or the image has internal distortion that cannot be modelled by the math model.

Analysis of Residual Errors of All GCPs
Scatter diagrams are used to plot the residual errors and can illustrate the trend of errors with respect to the cross track or along track directions.If there is a trend of changes instead of a random distribution of changes, a polynomial function may be estimated from the scattered data.A displacement map of the residual errors of all GCPs as vectors can then be used to visualize the patterns of changes across the entire image.

Operational Correction Methodology
If well approximating the photogrammetry of an image, a RPC model is able to model satellite images and may achieve subpixel accuracy in orthos (Fraser, et al., 2001) (Hu & Tao, 2004).
If not, the RPC could not correct the image up to a sub-pixel accuracy level.This is either because there is inaccurate camera calibration (Dinguirard & Slater, 1999) or because the field of view (FOV) is too wide to use the affine transformation to compensate an RPC (Xiong & Zhang, 2009).In theory the Toutin model should instead be capable of establishing the accurate geometric relationship but in reality always ends with a simulated model because of the important parameters missing regarding camera calibration.A few literatures (Wang, et al., 2011) (Toutin & Cheng, 2000) (Hu & Tang, 2011) (Toutin, 2004) try to use polynomials, thin plate spline functions or finite element functions for image correction.All those models have one known weakness that they do not take the elevation of the ground into consideration, which certainly cannot guarantee the final accuracy (Toutin, 2004).
This paper proposes using an operational approach for correction in addition to the operational accuracy assessment methodology.It uses all GCPs collected from the accuracy assessment methodology and then computes a new RPC model.A similar approach was discussed but questioned by a few researchers (Toutin & Cheng, 2000) (Hu & Tao, 2002) (Di, et al., 2003) (Xiong & Zhang, 2009) because it requires an extreme large number of GCPs that should be accurate, evenly distributed horizontally, and evenly distributed vertically.Those restrictions made it impossible to be operational under circumstances before.
A successful experiment using the RPC computed from GCPs was reported in the work of (Gianinetto & Scaioni, 2008), which yet could not ensure an even distribution of GCPs due to its choice of GCP collection technology (feature matching in spatial domain), and its application was made to two satellites (i.e., IKONOS, QB) that have high quality in geometry and radiometry.This paper extends the work in (Gianinetto & Scaioni, 2008) with the GCPs collected by the operational accuracy assessment methodology and applies to the four Chinese satellites.The computed RPC has the following equations: (5) where  ′ ,  ′ are compensated image coordinates using the collected GCPs with  ,  and  ground coordinates from Equations ( 3) and ( 4), and   ′ ,   ′ ,   ′    ′ ( = 0,1, … ,19) are the coefficients of four new polynomial functions.Linearization of Equations ( 5) and ( 6) with the method of least squares can resolve the coefficients (Di, et al., 2003) (Hu & Tao, 2002).
Although the computed RPC model could not fully recover the rigorous sensor geometry (Di, et al., 2003) (Xiong & Zhang, 2009), it is good enough to approximate the sensor model when some incompleteness of camera calibration or a large FOV exists (Gianinetto & Scaioni, 2008).

EXPERIMENTS
HJ-1 A/B CCD Level 1, CBERS-02C HRC Level 1, and ZY-3 TLC Nadir Level 1 are three popular commercial products from the four Chinese satellites for applications that require accurate orthos (CRESDA, 2014).Detailed information of parameters regarding their platforms and sensors can be found from the web source (CRESDA, 2014).The experiments conducted here select the above three products to prove the soundness of the proposed operational accuracy assessment and correction methodology.
The selection of scenes for each product is managed to meet the following criteria: 1.Each scene should be covered most by lands; 2. Each scene should be less or free of cloud/snow/ice; 3.All scenes for one satellite product should have the typical landscapes such as plains and hills/mountains.

HJ-1 A/B CCD Level 1
A total number of 31 scenes, mixed from both HJ-1 A and HJ-1 B satellites, are selected for testing.The nominal spatial resolution of those scenes is 30 meters.These scenes cover most land territories of China with landscapes of plains, hills/mountains, and mixture of plains and hills/mountains.See Table 1 for detailed information of those scenes.Since the swath of each HJ -1 A/B scene has 360km in width, the province(s) covered by each scene is/are listed in the table particularly.

Table 1 HJ-1 A/B testing scenes and accuracy results
There are about 500 GCPs automatically collected for each scene against a mosaic (in the 2000's) of the 15m Landsat 7 ETM orthos (NASA, 2014) over China.After fitting those GCPs to the simulated Toutin model and the RPC computed from those GCPs, it can be seen that the residual errors with the Toutin model at cross track are large and forms a second-order polynomial function trend, while those with the computed RPC are small and randomly distributed.See Figure 1 for the comparison.The rest of all combinations of RMS of X or Y with respect to the cross/along track show a random distribution of errors and are not shown here in order for shortening the paper.

Figure 1 HJ-1 A/B RMS of X cross track for the scene of HJ1A-CCD2-19-72-20130411. Left --by the simulated Toutin model; Rightby the computed RPC
All 31 scenes are then processed using either the simulated Toutin model or the computed RPC and have orthos generated.The 90m SRTM DEM (CGIAR-CSI, 2014) is used in orthorectification.The accuracy of the orthos is measured using nine human-picked points that generally consist of four points around the four image corners, four points in the centers of four image borders, and one in the middle of the entire image.
Table 1 lists the residual errors for each image.The residual errors of all scenes in this test distinguish the computed RPC from the simulated Toutin model.It can be seen that the proposed operational correction methodology is much more stable than the simulated Toutin model and achieves about one and half pixels (i.e., 45m on ground) accuracy at average with half pixel standard deviation, while the simulated Toutin model is very unstable because half of its orthos have more than 100 pixels error at average and the rest have about three pixels errors at average.

CBERS-02C HRC Level 1
The testing data set consists of 30 scenes (15 pairs of images from left and right cameras) of HRC Level 1 product: 22 scenes with most hilly/mountainous landscapes and 8 scenes with most flat areas.See Table 2 for detailed information of these scenes.These scenes have a nominal 2.36m spatial resolution.

Table 2 CBERS-02C testing scenes and accuracy results
Since those commercial products are distributed with vendorcomputed RPC files (referred as to the original RPC), the models to be compared are an original RPC and a computed RPC.The reference imagery is a mosaic from ZY-3 TLC Nadir orthos (at 2m ground resolution with 3m accuracy at average) generated before.About 500 GCPs are automatically collected for each scene and fit to the original RPC and the computed RPC.Plots of the RMS of each scene show that there are systematic distortions for X and Y values at the cross track when using the original RPC (see Figure 2 as an example), while only random distortion is found when using the computed RPC (see Figure 3 as an example).
When plotting the displacement (the RMS of X and Y as a vector) map of all GCPs using the original RPC model, it can be seen that the distortions across an image are not uniform but form regular patterns.Figure 4 is the displacement map of the same scene as that used in  The geometric accuracy of CBERS-02C HRC orthos using the original RPC is about 4 pixels (10m on ground) at average, while that of the orthos using the computed RPC can achieve one and half pixels accuracy (4m on ground) at average.See Table 2 for RMS of all scenes for comparison.

ZY-3 TLC Nadir Level 1
The testing data set consists of 28 scenes: 13 of them have plains covering most areas of each scene; 15 of them have hills/mountains covering most areas of each scene.Those images were acquired from May 2013 to April 2014.See Table 3 for detailed information of these scenes.These scenes have a nominal 2.1m spatial resolution.
Since those commercial products are distributed with original RPCs, the models to be compared are an original RPC and a computed RPC.The reference imagery is a mosaic from ZY-3 TLC Nadir orthos (at 2m ground resolution with 3m accuracy at average) generated before.

Table 3 ZY-3 testing scenes and accuracy results
With about 500 GCPs automatically collected for each scene, the fitness of those GCPs to either the original RPC or the computed RPC is analysed.The plots for the RMS with respect to the cross or along track directions shows no systematic trend but looks almost like random distributions.Figure 5 is an example of the RMS analysis for one ZY3 scene that shows no significant difference between both methods.The same correction procedure is applied to those scenes and the accuracy of each ortho is measured manually in the same way as that in the experiments for HJ-1 A/B and CBERS-02C.
The manual measurement shows that the computed RPC can achieve half pixel (1.2m) accuracy at average, while the original RPC achieves about one pixel (2.2m) accuracy at average.Both results are however acceptable in application productions.Table 3 lists the RMS for all scenes for comparison.

CONCLUSIONS AND DISCUSSIONS
The proposed accuracy assessment methodology is demonstrated via tests against a significant number of samples to be operational in three folds: first, it is able to draw a reliable assessment conclusion regarding the commercial satellite images for production; second, it requires minimum manual labour, particularly regarding the full automation of GCP collection when compared with traditional methodologies; third, the assessment results can be used to compute a new RPC model that can lead to significant improvement of geometric accuracy of orthos.
The causes for the unfitness of the collected GCPs to a math model in the satellites assessed here are either from incomplete calibration of cameras in the given information (the original RPC from CBERS-02C HRC data) or missing critical information that makes a simulated physical model (simulated Toutin model for the HJ-1 A/B scenes) incapable of modelling images, particularly with a large swath for one image.
A RPC model computed from GCPs that are collected from reference imagery works well to model the image particularly when the given model (either the original RPC or the simulated Toutin model) is unsuitable or information is missing.It improves the accuracy of orthos for production.Particularly, the computed RPC model improves the accuracy of HJ-1 A/B products dramatically (more than ten times accurate than the simulated Toutin model) and that of CBERS-02C significantly (more than twice accurate than the original RPC); the orthos for both HJ-1 A/B and CBERS-02C using the computed RPC model meet requirements for productions (i.e., 2 pixels at average).The accuracy improvement for ZY-3 is small when using the computed RPC model in comparison with the original RPC one, which is however not critical in daily production.
According to daily production results (when the criteria of selecting images for the experiments conducted in this paper is much relaxed), the ratio of good orthos over the total processed scenes is about 80%.Although this ratio is acceptable to clients for their daily production work, there are rooms for enhancing the relevant technologies to achieve higher success rate in the future.
One issue for increasing the success rate was identified that the computation of a new RPC from GCPs becomes unstable and could lead to significant distortion in parts of an ortho when the distribution of GCPs is not even horizontally or vertically.It is yet challenging to make GCPs distributed evenly in both horizontal and vertical directions when full automation is required in daily production.The future work will explore an operational approach to ensure the even vertical distribution while keeping an even horizontal distribution, which is believed to be able to warrantee the stability of the computed RPC model.

Figure 2 .
The pixels and lines of the image are scaled to 5 times smaller while the values of displacement of X and Y are magnified by 10 in the displacement map.There are two major displacement patterns: one on the left of the image has displacement toward right; another one on the right of the image has larger displacement toward the lower right image corner.

Figure 2
Figure 2 RMS by the original RPC for scene ZY02C_HR1_E112.4_N34.8_20131107