BLOCK ADJUSTMENT OF MULTISPECTRAL IMAGES WITHOUT GCP AIDED BY STEREO TLC IMAGES FOR ZY-3 SATELLITE

Multispectral images are the main data source of optical satellite remote sensing application. Consistent geometric accuracy is the basis of image registration and fusion. But the multispectral images usually collected by the nadir imaging camera are of weak geometric intersection, which will lead to the rank defect of the adjustment equation when performing block adjustment (BA) directly, making the solution unstable. Thus, a planar BA aided by additional digital surface model (DSM) which can overcome this weak geometry is often used to improve the geometric accuracy and consistency of regional planar images. However, the inevitable elevation error and the indispensable ground control points (GCPs) make the method limited in practical application. In this paper, a new method aiming to the BA of planer multispectral images of ZY-3 satellite without use of GCPs is presented. This method introduces the constraint of stereo images to assist the BA of planar images. By configuring appropriate weights of different observations, the integrated optimization of positioning models of planar and stereo images can be achieved together. The effectiveness of the method was verified by the planar multispectral images and stereo three linear camera (TLC) images collected by ZY-3 satellite. The satisfactory results indicated the rationality and effectiveness of the presented method.


INTRODUCTION
The multispectral images of ZY-3 satellite contain four spectral bands (red, green, blue, and near-infrared) and have a spatial resolution of 5.8 m (Li, 2012). The approximate nadir pushbroom imaging collects the images with a swath width of 52km, thus, the maximum intersection angle of adjacent images is about 5.8 degrees, which is an extremely weak stereo geometry. Consistent geometric accuracy is the basis of subsequent application of multispectral images. Although ZY-3 satellite has been calibrated in orbit, the attitude error and some random influence on the satellite still cause significant accuracy differences between multispectral images and between multispectral and panchromatic images, attenuating the reliability in subsequent applications. Aiming to eliminating this accuracy differences and ensuring the geometric consistency among the ZY-3 images, we developed an effective block adjustment (BA) method to improve the geometric accuracies of multispectral and panchromatic images together.
According to the intersection geometry of the images, BA method for satellite remote sensing images can be divided into stereo BA (Grodecki and Dial, 2003;Rottensteiner, et al., 2009) and planar BA (Teo, et al., 2010;Pi, et al., 2019). The former is developed for such images with sufficient intersection geometry (base-height ratio greater than 0.6 or intersection angle greater than 30 degrees), such as the three linear camera (TLC) images of ZY-3 satellite, the HR images of SPOT5 satellite; from which an adequate BA model with stereo measuring function could be established simply (Yang, et al., 2017). However, the method is not suitable for the case when satellite images are acquired in an approximately nadir viewing mode. In such case, the intersection * Corresponding author geometry among these images is relatively weak due to the characteristics of high orbit and narrow field of view. High resolution nadir images and multispectral images usually have such geometric characteristics, and these images with high spatial resolution and spectral resolution are the main data sources for the follow-up satellite remote sensing applications. The latter planar block adjustment is always adopted for the geometric process of these planar images, in which an available digital surface model (DSM) is usually used as elevation constraint to overcome the weak intersection geometry, suppressing the instability in calculation. However, a qualified DSM reference data is not always available readily in practical application, especially in those region blank of mapping. Some open DSM data at global scale, such as SRTM and ASTER, is restrained by limited sample size and geometric accuracy, and invalid in improving the BA accuracy sometimes, especially in those region of complex terrain. The reason is that in the planar BA, the elevation error will cause the accumulation of planar error in the block, resulting in the decline of the BA accuracy. Furthermore, due to the initial geometric positioning error of satellite images, it is necessary to use ground control points (GCPs) to ensure the registration of satellite image and DSM. Otherwise, the elevation error caused by the dislocation will still reduce the accuracy and reliability of BA. Therefore, this method is not suitable for the BA without use of GCPs and limits its application.
In this paper, we present a new BA method for planar multispectral images without use of GCPs, inspired by the fact that a desirable stereo observing model with sufficient intersection can be derived from stereo satellite images. Through reasonable tie points (TPs) matching, the introduced stereo images can effectively improve the weak intersection of block, enabling stable and accurate estimation. In this method, we adopt the unified rational function model (RFM) as basic imaging model for both planar and stereo images, and establish the BA model by adding a correction model in the imagery space of RFM. On this basis, this paper describes a series of key technologies involved in this method, including establishment of BA model, weight setting of various observations and optimal estimation of unknown parameters. Compared with the traditional planar BA method, this method can not only overcome the weak intersection geometry of planar images, but also will not introduce elevation interpolation error since its elevation constraint directly comes from stereo satellite images. Additionally, the elimination of DSM makes it suitable for the BA without the use of GCPs and more practical.

METHODOLOGY
The bundle BA method is used to perform the mixed adjustment of multispectral planar images and TLC stereo images, and the RFM model of single image was used as the basic adjustment unit (Tao, et al., 2001). The tie points (TPs) were identified among the adjacent images automatically, and the virtual control points (VCPs) (Yang, et al., 2012) were generated to overcome the freedom of the block result from the lack of GCPs. According to the adjustment theory, the weight of the TPs, VCPs and elevation constraint were determined by prior knowledge and updating according to the posteriori accuracies. With the TPs, VCPs, original RPCs and estimated ground coordinates of TPs as inputs, the BA model then could be established, therefore the calculation of the BA parameter was conducted. The specific approach of the BA method without use of GCPs is shown in Figure. 1.

Matching TPs
Establishing modified normal equation  Figure 1. Approach of the BA without use of GCPs

Block adjustment model
Geometric imaging model is the fundamental mathematic model in BA for optical satellite images, which establishes the relationship between each image point ( , ) ls and its objectspace counterpart ( , , ) B L H . The RFM, fitted from the rigorous physical model (Madani, 1999), is widely used in the geometric process of satellite images due to its simple and unified form. Therefore, we also adopt the RFM for the presented BA method. Through attaching a suitable error correction model ( , ) ls  into the image space of RFM, the basic BA model can be established, as in Eq.(1).
The choice of error correction model depends on the estimation of remaining geometric error existing in the test images. The commonly used models are translation model, affine model and quadratic polynomial model. For the multispectral images and TLC images of ZY-3 satellite, the highorder geometric errors are compensated well by on-orbit geometric calibration in commissioning phase (Wang, et al., 2014;Zhang, et al., 2014), thus an affine model is usually sufficient to correct the remaining errors.
The essence of BA without use of GCPs is the optimization based on a rank defect model. In general, supplementary constraints need to be developed to improve the BA model, and ensure the stable solution of BA. Here, we introduced the VCPs generated with the initial RPCs of images as weighted observations to improve the condition of BA model. The VCPs could be regarded as the GCPs with smaller weights. Therefore, the observations in our method include not only the TPs, but also the VCPs. By linearizing the adjustment model Eq. (1), we established the error equations for the VCP and TP, respectively. The difference between them is that the unknown parameters corresponding to the TP include not only the coefficients of the correction model, but also the corresponding ground point coordinates. As in Eq. (2) where t is the correction for the correction model coefficients,

Weight strategy
In BA, the contribution from each observation is controlled by the weight matrix. According to the classical least-square adjustment theory, each observation should be assigned to a reasonable weight according to its variance to ensure the optimal estimation of unknown parameters. However, in most cases, the accurate variance of observation is unavailable, and an initial weight is set empirically and roughly in terms of some prior knowledge. In the combined error equation Eq. (2), the observations involved in presented BA were divided into three groups, which are the VCPs of TLC images, VCPs of multispectral images and the TPs. The observations in different groups are independent of each other. Supposing that the observations in the same group have the same variance, thereby only three kind of weights need to be determined.
The weights of the TPs could be estimated according to the matching accuracy by high-precision matching operators. In general, the matching accuracy of TPs among remote sensing images is better than one pixel, thus the initial weight of TP can be set to 1.0 directly. The weights of VCPs are difficult to determine due to the lack of prior information, which is critical for the optimal estimation of the BA parameters. If the weights are set too large, the power of the TPs will be weakened, resulting in relative error among images. Conversely, the weights are set too small, the freedom of the adjustment parameters fail to be constrained desirably, leading to poor convergence and unstable estimation. For ZY-3 images, a great number of studies have demonstrated that their initial RPCs provided by the vendors has a geopositioning accuracy of about 20m after geometric calibration (Tang, et al., 2012;Cao, et al., 2015), thus the initial weights of VCPs could be set to 0.0025 according to the relationship 2 1/ P   ( 2  is variance) between weight and accuracy. However, this weight is not constant, because with the iterative solution of BA, the image accuracy continues to improve, if the weight of the VCPs is still too small, it cannot play a role in constraining the deformation and error accumulation within the block, so it is necessary to update the weight according to the current imagery posteriori accuracy. It should be noted that this is only the initial weights of observations, and the final weight need to be balanced according to the number of TPs and VCPs.
The weight of the elevation constraint depends on the intersection condition of the TPs. When the maximum intersection angle of the corresponding light rays in a TP is greater than 30 degrees, it is considered to be a good stereo intersection. Thus, the elevation constraint can be ignored. When the intersection is weak, a larger weight should be configured to ensure the stability of the elevation estimation. In this method, the elevation constraint weight of the TP with less than 0.1 arc intersection angle is set as 0.01, and that with larger angle is set as 0.005.

Estimation of block adjustment
Based on the established error equation, we adopted the least square method to calculate the adjustment parameters iteratively. For all TPs, the current RPC parameters were used to determine the corresponding ground 3D coordinates by forward intersection, and for the TPs of weak intersection, the ground 3D coordinates were estimated under the weighted elevation constraint. Then, the estimated coordinates were introduced into the BA model as the initial values to solve the unknown parameters. In the BA model, there are two kinds of unknown parameters (coefficients of correction model and the ground coordinates of TPs) need to be estimated, but they are not calculated together. In the least square solution, the unknown ground coordinates with lager number were eliminated to construct modified normal equation with only parameters of the error model. In the solution of the modified equation, the ternary storage structure based on sparse matrix was used to reduce the cost of storage and the complexity of data organization, and the conjugate gradient method was used to improve the efficiency of solution. With the estimated parameters, the RPC model were updated and then the ground coordinates were re-estimated. The adjustment calculation is an iterative process performed until the difference between two successive results is less than a predefined tolerance criterion.
Additionally, a stepwise solution strategy from low-order to high-order parameter was developed for the estimation of correction model parameters. In the initial solution, when the weight of VCPs is small, only the translation error is calculated, and then other high-order errors are further estimated with the updating of the weight. The purpose of this is to prevent the deformation and error accumulation in the image block caused by the unequal distribution of TPs when the weight of VCPs is small.

Description of dataset and test area
We choose the Weifang in Shandong Province of China as the test area to verify the present BA. A total of 16 multispectral images with corresponding TLC stereo image pairs were used in this BA, and the corresponding RPC files were provided. These images were collected between from January 2014 to September 2018. The overlap along and across the track is greater than 15% and 50%, respectively. Based on the high precision matching algorithm (Temizel and Yardimci, 2011), a total of 6649 reliable TPs were identified, the distribution of these images and the identified TP over the test area is shown in Figure.  Based on the identified TPs and generated VCPs, the BA of the multispectral images and stereo TLC images were achieved together under the condition without use of GCPs. All the BA calculations were achieved on a personal computer in one minute.
Detailed information about the test area is listed in Table 1.

Parameter Value
Size of test area About 22,100 km 2

Number of images 64
Number of TPs 6649 Topographic relief About 400 m Table 1. Basic information of the experimental area.

Accuracy assessment and analysis
Since the BA was performed without the use of GCPs, which is unable to improve the absolute geo-positioning accuracy of images, we focused on the assessment of the relative geometric accuracy. The relative geometric accuracies of the images before and after BA were compared in this experiment. Two kinds of quantitative indexes were used to evaluate the relative geometric accuracy of the block. One is the statistical accuracy of the relative residual of the TPs on each image, and another is the relative geometric accuracy between adjacent images. The original RPCs and that refined using the estimated parameters were used to verify the accuracy before and after BA, respectively. The image relative residuals of a TP on an image is the bias between the measured and estimated image points. After determining the relative residuals of all TPs, the root mean square error (RMSE), mean value (Mean), and maximum error (Max) of the residuals could be counted, as listed in Table 2.  Table 2. The statistics of the relative image residuals of TPs in the whole block before and after BA.

Stage
The overall relative error of the block before BA is about 13 pixels according to comprehensive RMSE. The large differences between the mean values and the maximum error indicate that the original geometric accuracies of the images are erratic. After the presented BA, the RMSE of the relative residuals was improved to better than 0.5 pixels, and the mean values of the residuals were almost zero, indicating the presented model and the solution strategy such as weight determination and stepwise estimation is reasonable. Additionally, the small deviations between mean values and maximum errors after BA shows that the images in the block have the consistent geometric accuracy, and further verify the effectiveness of the method.
To further verify the relative geometric accuracy of BA, the splicing accuracies of each adjacent image pair before and after BA were also evaluated. The corresponding points identified in the overlapping area of adjacent images were used as checkpoints to assess the relative geometric accuracy. The root mean square (RMS) of the relative residuals for the checkpoints was used as the index to describe the relative geometric accuracy of each image pair. A total of 10 pairs of images distributed in the whole block were selected for the accuracy assessment. For each pair of images, the relative geometric accuracy between multispectral images (M-M), between panchromatic images (P-P), and between multispectral and panchromatic images (M-P) before and after BA were verified. The distribution of these splicing accuracies is shown in Figure. 3.
It can be observed in Figure. 3 that the distribution of the relative geometric errors between multispectral images, between panchromatic images, and that between multispectral and panchromatic images all are disordered before BA, thus the initial geometric accuracy of the images in the block is extremely inconsistent. After the presented BA, the three kinds of splicing accuracies of the image pairs all are improved to better than 1 pixel, and the distribution of the accuracies is much more consistent, which indicate the effectiveness of the presented BA method. For an intuitive comparison of the improvement in relative geometric accuracy, we performed a vision comparison of splicing between multispectral images ( Figure. 4(a)), between multispectral and panchromatic images ( Figure. 4(b)), and between panchromatic images ( Figure. 4(c)), before and after BA. As shown in Figure. 4, in which improvement of the relative geometric accuracy is apparent.

CONCLUSION
A mixed BA method was presented in this paper aiming to the geometric correction of the planar multispectral images. With the supplement of stereo TLC images of ZY-3 satellite, we overcome the weak intersection geometry of the planar images by a weighted constraint BA model. Through a reasonable weighting strategy and a stepwise solution method, the optimal optimization of relative geometric accuracy among the multispectral images and that between the multispectral images and panchromatic images were achieved under the condition without use of GCPs and additional DSM. Based on our method, the BA of 64 multispectral and TLC images of ZY-3 satellite was performed. The relative geometric accuracy over the whole block is within a pixel, which meets the demand for seamless mosaic. The satisfactory experimental results indicate the rationality and effectiveness of this method.