INFLUENCE OF IMAGE INTERPOLATION ON IMAGERY-BASED DETECTION AND COMPENSATION OF SATELLITE JITTER

Satellite jitter is a common and complicated phenomenon that degrades the geometric quality of high-resolution satellite images. Imagery-based detection and compensation of satellite jitter have recently been widely concerned. However, most of the existing studies overlook the issue of image interpolation in this topic involving subpixel measurements. In this study, the influence of image interpolation on imagery-based detection and compensation of satellite jitter is investigated. Four different interpolators are separately applied in dense least squares matching for jitter detection based on parallax observation and in intensity resampling for jitter distortion compensation. Experiments were carried out using ZiYuan-3 dataset to compare and analyze the results in the case of different image interpolation. The experimental results demonstrate the influence of image interpolation on imagery-based jitter processing. Inferior interpolators can induce pixel locking effect in subpixel matching and position-dependent systematic bias after image correction, which deteriorate the performance of jitter detection and compensation. To ensure the reliability, sophisticated interpolation algorithms with smaller phase errors are preferable in imagery-based jitter detection and compensation. * Corresponding author


INTRODUCTION
Precise sensor orientation plays a crucial role in the geometric processing of high-resolution satellite images. Spacecraft position and absolute attitude can be determined by the onboard global navigation satellite system receivers, star trackers and gyros with post-processing strategy (Tang et al., 2015). The sensor misalignment due to camera mounting and sensor distortions can be corrected by in-orbit geometric calibration. Besides, satellite jitter that causes the periodic instability of satellite attitude is another significant error source of sensor orientation (Tong et al., 2014). Many high-resolution satellites suffer from the issue of satellite jitter that is possibly induced by both external environment variations and internal mechanical operations (Mattson et al., 2010;Tong et al., 2017;Zhu et al., 2019). A subtle variation on attitude angle will cause a considerable geometric error on the ground. Therefore, detection and compensation of satellite jitter is an essential step to ensure the geometric quality of high-resolution satellite images.
The attitude measurement data is the most direct source to extract and estimate the attitude jitter information, if the sampling rate of the on-board attitude sensors is larger than double frequency of jitter and the measurement precision is enough high Zhang, Guan, 2018). However, in the cases when the performance of attitude sensors is inadequate, it is necessary to realize the jitter detection task with additional manners. For linear array pushbroom cameras, each image line scans the ground at different time, and the dynamic micro-vibration of satellite platform will lead to different distortions between adjacent imaging lines . Therefore, satellite imaging sensor provides an alternative solution to detect and estimate satellite jitter since the distorted images record information about satellite jitter. The commonly used method based on imagery is to obtain the satellite jitter information from the dense matching results of images with parallax observation such as multispectral bands and staggered arrays (Teshima, Iwasaki, 2008). This method is convenient to be implemented without the need of external reference, and has been widely applied and extended for jitter detection of many sensors and satellites, such as ASTER, MRO, Pléiades-HR, ZiYuan-3, Mapping Satellite-1 and Gaofen-1 (Teshima, Iwasaki, 2008;Amberg et al., 2013;Sun et al., 2015;Tong et al., 2015b;Pan et al., 2017;Sutton et al., 2017;Zhu et al., 2019). In addition, to remove the negative effect of satellite jitter, the estimated jitter information can be compensated in the image space. This correction can be performed by directly resampling the original distorted images (Tong et al., 2014) or reprojection with sensor models . The corrected jitter-free images are beneficial to the following applications, such as band-to-band registration, change detection and RFM replacement (Ayoub et al., 2008;Wang et al., 2016;Tong et al., 2019).
Imagery-based detection and compensation of satellite jitter have been received a lot of attention in the study area of spaceborne photogrammetry. Since the magnitudes of jitter distortion in the image space are normally less than certain pixels, fractional measurement is required and subpixel-level accuracy is highly pursed in this topic. Therefore, image interpolation that transforms the discrete matrix to continuous image representation is indispensable in both jitter detection and jitter compensation. However, few previous studies have attempted to consider the influence of different image interpolation algorithms. In this study, we focus exclusively on this issue, investigating the influence of image interpolation on imagery-based detection and compensation of satellite jitter, and intend to reveal the internal explanation. Jitter detection method based on parallax observation and jitter compensation method based on direct image resampling are adopted. In this case, the performance of four types of interpolation algorithms are compared and analyzed in dense least squares matching and intensity mapping using the experiments conducted with ZiYuan-3 multispectral images. It is found that image interpolation algorithm is an affecting factor that determines the final results of jitter detection and compensation. The possible reason for the diverse effects and systematic errors of different interpolation algorithms is further discussed.

Jitter Detection and Compensation
The overall workflow of the adopted imagery-based detection and compensation of satellite jitter is illustrated in Figure 1. The distortions caused by satellite jitter in the image space are retrieved from the relative disparities between parallax observation, and accordingly corrected on separate image by intensity interpolation.

Jitter Detection:
For a parallax observation system that contains multiple linear arrays that are imaged with a slight time difference t ∆ , satellite jitter is transmitted to the focal plane and will induce relative registration errors between different linear arrays. The relationship between the image distortions caused by satellite jitter ( ) f t and the relative disparity ( ) g t can be expressed by (Teshima, Iwasaki, 2008) In order to estimate the jitter distortion ( ) f t , the relative disparity ( ) g t is firstly obtained by high-density area-based matching. In this study, image matching is implemented on a dense grid. For each grid position, the initial disparity is provided by means of quadric fitting of normalized cross correlation similarity measure within corresponding patches, and least squares matching is further adopted to refine subpixel accuracy (Gruen, 1985).
In least squares matching, considering the affine mapping function and linear radiometric difference, the relationship between the reference patch ( , ) r x y and target patch ( , ) t x y is given by After expanding the nonlinear form with a first-order Taylor series, the residual function is given by where dx is the incremental parameter vector, C is a coefficient matrix that records the partial derivatives (gradient) to each unknown parameter, and L is the reduced observation. The incremental parameter vector can be solved using least-squares adjustment, and the unknows are iteratively updated by until reaching a certain termination condition (Gruen, 1985). The translation term 0 0 ( , ) a b is used to determine the subpixel shift between patches.
The disparity maps in the cross-track and along-track directions can be obtained after high-density matching with least squares matching implemented on patches distributed densely. The matching outliers are filtered out by checking if the normalized cross correlation is higher than a predefined threshold (e.g., 0.6).
The filtered disparity maps are further compensated with affine model to attenuate the relative linear bias resulted from other influence factors except jitter . Finally, the relative disparity ( ) g t can be obtained by averaging the disparity map at each line with three-sigma limits, and is then smoothed by a low-pass function.
According to the relationship in Equation (1), the absolute jitter distortion ( ) f t is retrieved from ( ) g t in the frequency domain based on the Fourier shift property (Tong et al., 2015b): where i is the first solution to the equation 2 , and  and 1 −  denote the Fourier transform and the inverse Fourier transform, respectively.

Jitter Compensation: The derived jitter distortions
( ) x f t in the cross-track direction and ( ) y f t in the along-track direction in the image space of multispectral sensor can be transformed to the variations of roll and pitch angle relative to the orbit coordinate system based on the principal length and off-nadir angle of multispectral sensor (Liu et al., 2016).
With these attitude angle variations, the image correction to eliminate the influence of satellite jitter can be realized through two ways: (1) directly resampling the image needed to be compensated after determining the imaging time range and the corresponding jitter distortions in the specific image space; (2) based on the sensor model, projecting the image to a nominal surface with jittery attitude and re-projecting back to the image space using jitter-free attitude Zhang et al., 2018). Both two ways can effectively realize jitter compensation on imagery and involve intensity interpolation. In this study, the former one without the use of sensor model is adopted for simplicity to straightforward show the influence of different intensity interpolation. With the determined time (line) range t  , the corrected image I  can be mapped from the distorted image x y

Image Interpolation
There are two places in the imagery-based jitter detection and compensation that intensity interpolation is required. The first one is in the process of jitter detection in that the iteration of least square matching needs to calculate the intensity and gradient at the fractional positions. The other one is in the process of jitter compensation that needs to calculate the intensity from the distorted image.
The essence of image interpolation is to reconstruct a 2-D continuous signal s(x, y) from its discrete samples s(i, j), where (x, y) is within real domain and (i, j) is within integer domain (Thévenaz et al., 2000). The ideal interpolator is sinc function. However, it is normally not used in practical due to its infinite support. For the computational speed, symmetrical and separable interpolation kernels with a smaller support is used instead (Lehmann et al., 1999): where h(x) is the interpolation kernel. In this study, four linear methods for image interpolation are adopted including nearest neighbor, bilinear, bicubic, B-spline interpolation. The interpolation kernels of these four interpolators are described in the following.
The easiest way is rounding to the nearest integer. The interpolation kernel for nearest neighbor is a rect function: For separated bilinear interpolation, the values are weighted sum of both direct neighbors. The interpolation kernel follows the triangular function: The bicubic interpolator uses the cubic polynomials to approximate the sinc function, and the interpolation kernel for four-point interpolation is expressed by where a is a free parameter and is set to -0.5 as recommended in Keys (1981). B-spline interpolation is different from the above three interpolator. It belongs to the type of generalized interpolation the replaces discrete samples with interpolation coefficients, and can be solved by a two-step method (Thévenaz et al., 2000). The cardinal kernel of cubic B-spline interpolation is given by Figure 2a shows the interpolation kernels of different interpolation algorithms. It can be seen that cubic B-spline interpolator is closer to the ideal sinc function than others. The frequency responses of various kernels after Fourier transform are plotted in Figure 2b. The frequency analysis indicates the differences of frequency characteristics in both passband and stopband, which causes the different performance of four image interpolation algorithms in jitter detection and compensation.

Experimental details
To investigate the influence of intensity interpolation in imagery-based jitter processing, experiments with Ziyuan-3 multispectral images were conducted. The multispectral sensor was taken as a parallax observation system to detect and estimate attitude jitter of Ziyuan-3 satellite, and the distortions caused by satellite jitter in different bands were corrected by image resampling. The Ziyuan-3 multispectral sensor provides four spectral bands (B1-B4: blue, green, red, and near-infrared) with a ground sampling distance of 5.8 m at a nominal altitude of 505 km and a line integration time of approximately 0.8 milliseconds. As shown in Figure 3, a raw level-0 multispectral image with relatively flaw terrian relief was used as experimental dataset.
Dense image matching with least squares matching was implemented on a dense grid with steps of 2×5 pixels on the first two multispectral bands (B1 and B2). Four types of interpolation algorithms were tested in this study. The results of jitter detection and jitter compensation using each interpolation algorithm were calculated and compared.

Results of Jitter Detection
After dense matching and mismatch removal, the relative disparity curves are obtained by averaging the disparity maps at each line. The results of dense least squares matching with four different interpolators are presented in Figure 4. The periodic trends in the disparity curves and Fourier analysis demonstrate the existence of satellite jitter with frequencies of approximate 0.65 Hz in both directions in this scene.
In addition, the comparison of different interpolators indicates that the matching results of least squares matching with nearest neighbor interpolator obviously deviate from the other three interpolators, while the results of the others have a slight difference. The deviation of nearest neighbor interpolation refers to a systematic error called pixel locking effect (Shimizu, Okutomi, 2005;Tong et al., 2015c). It can be found that the subpixel estimates towards to bias toward integral pixel positions. This leads to the situations that the results of nearest neighbor interpolation are less or higher within the different intervals of 0.5 pixels, and all the results coincide at the approximate multiples of 0.5 pixels. This systematic error is related to the image interpolation in matching algorithms, and also depends on the image content. For this scene, bilinear interpolator is adequate, although spline interpolator would be better. The dense matching is crucial to image-based jitter estimation. Based on the obtained relative disparity, the absolute distortions caused by satellite jitter can be retrieved according to Equation (4), and are then transformed to the time-varying attitude angle variations relative to the orbit coordinate system. Figure 5 shows the results of jitter detection and estimation with different interpolators. Attitude jitter components estimated from the onboard attitude measurement data with a sampling rate of 4 Hz are also plotted in Figure 5 as reference. The jitter distortions in the cross-track direction corresponds to the variations of roll angle, and distortions in the along-track direction corresponds to the variations of pitch angle. Note that the direct-current components of the estimates from imagery and attitude data have been eliminated.
Similar as the results of relative disparity, the derived absolute attitude jitter components using jitter detection method with nearest neighbor interpolator are different from the results from the other three interpolators. Table 1 gives the comparison results with the discrete samples of attitude measurement data. The root-mean-square deviation of the result of nearest neighbor interpolator is 0.362 arc second for roll angle and 0.216 arc second for pitch angle. In contrast, the root-meansquare deviations of the results of the other three interpolators are 56.1% and 68.1%, 54.7% and 66.7%, and 52.7% and 62.5% of the deviation of nearest neighbor interpolator for roll and pitch angle, respectively. It is obvious that the more sophisticated image interpolation leads to the better performance of jitter detection and better consistence with the results of attitude data.  Table 1. Root-mean-square deviation compared with attitude measurement data (Unit: arc second).

Results of Jitter Compensation
Based on the derived jitter distortions in both directions and the corresponding time range, the multispectral images can be corrected. Firstly, the reliability of the jitter compensation method by direct image resampling was validated using a jitterfree reference image. As shown in Figure 6a, a Ziyuan-3 nadir panchromatic image covering the same region but captured at different date that has been demonstrated to be free from satellite jitter (Tong et al., 2015a) was adopted as the reference image. Two images were co-registered using a feature-based parametric alignment. The disparities between the warped reference image and the multispectral image were calculated separately before and after jitter compensation. After outlier removal and affine model compensation as described in Section 2.1.1, the averaged disparities in the cross-track direction over time were presented in Figure 6b. It can be seen that there exists an apparent periodic trend with a similar frequency of 0.65 Hz before compensation, which is resulted from the satellite jitter of the multispectral image. In contrast, this periodic varying trend basically disappears after jitter compensation on the multispectral image. The results demonstrate the effectiveness of imagery-based correction of satellite jitter.
The band-to-band registration between the corrected B1 and B3 were conducted to investigate the influence of image interpolation on jitter compensation. The dense matching scheme using least squares matching with B-spline interpolator was applied to calculate the registration errors between the corrected bands. Figures 7 and 8 show the registration errors averaged in each line between the compensated bands using different interpolators. Besides, the jitter distortions to be interpolated on B1 and B3 are also plotted in Figures 7 and 8 to reflect the relationship between interpolation artifacts and subpixel positions.  Theoretically, the residuals after jitter compensation should be the constant direct current components. The performance of nearest neighbor interpolator is intuitive. Since it rounds to the nearest integer, the registration errors equal to the original disparities when the subpixel interpolated positions of two images at the same intervals of 0.5 pixels, while appear obvious artifacts at the opposite positions. Note that it worsens the results before compensation if we look along-track direction. Interestingly, the calculated registration errors using other interpolators also present systematic artifacts to some extent. The registration errors all coincides when the interpolated distortion of either image meets the multiples of 0.5 pixels, and the magnitude of registration errors is related to the subpixel interpolated positions of both images (the detailed relationships can be found in the following subsection). Table 2 lists the standard variations of the registration errors before compensation and after compensation with different interpolators in two directions. The bilinear and bicubic resampling produce the similar systematic artifacts, and the Bspline resampling achieves the best results, which is almost 9 times better than nearest neighbor interpolator.

Discussion
According to the experimental results, different image interpolation can definitively affect the performance of imagerybased detection and compensation of satellite jitter. The inferior interpolation algorithms present evident systematic bias dependent on the subpixel positions in the case of both jitter detection and jitter compensation. A possible explanation of these systematic errors lies in the phase error introduced by the linear interpolation filter (Schreier et al., 2000). The magnitude of phase error depends on the frequency content of image represented as the wavelength λ, the fractional position δ, and the transfer function of different interpolator ( , ) H λ δ (Baldi, Bertolino, 2015): = arg[ ( , )] 2 / H φ λ δ δπ λ ∆ − The systematic positional error is accordingly given by = 2 p λ φ π ∆ ∆ (12) Figure 9 shows the theoretical positional errors of bilinear, bicubic and cubic B-spline interpolating functions for λ=3. It can be seen that the positional errors of bilinear and bicubic functions are same, and cubic B-spline function generates smaller errors. The positional errors vary with subpixel position, which cause the pixel locking effect in dense matching and systematic artifacts after image correction. To decrease the systematic errors, using high-order interpolators and reducing the high frequency content of image are recommended.

CONCLUSION
In this study, satellite jitter detection based on parallax observation and jitter compensation based on imagery with four image interpolation algorithms were performed and compared. Image interpolation is applied in iterative least square matching and direct image resampling in the process of jitter detection and jitter compensation, respectively. Experiments conducted using image data and attitude data of ZiYuan-3 satellite indicate that diverse interpolators achieve different performance of jitter detection and compensation. Inferior interpolators can generate larger subpixel position-dependent errors. This induces the pixel locking effect of image matching which biases the jitter detection, and leads to the systematic artifacts after image correction which limits the potential of the compensated images in applications that require high subpixel accuracy, such as band-to-band registration and change detection.