A Sensor-Shift Image Motion Compensation Method for Aerospace Remote Sensing Cameras Based on Image Motion Velocity Field Calculations

The calculation and analysis of the image motion velocity (IMV) field holds significant importance in the image motion compensation (IMC) of aerospace remote sensing cameras (ARSCs). Thus, this article puts forward a method for calculating the IMV field based on the rigorous imaging model, which takes into account the camera distortion characteristics. The proposed technique is applied to both a virtual and a real remote sensor to analyze the spatial and temporal features of the IMV field and the influence of camera distortion on it. Our experiments revealed that the additional IMV caused by camera distortion should not be disregarded when calculating IMV field due to its relative magnitude and the decrease in IMC performance caused by it. Additionally, we discussed the selection of the sensor-shift IMC strategy and found that for the real remote sensor considered in this article, the 2-D local compensation is already sufficient to achieve the desired compensation effect.


A Sensor-Shift Image Motion Compensation Method for Aerospace Remote Sensing Cameras Based on Image Motion Velocity Field Calculations Yan Bai and Zhonghe Jin
Abstract-The calculation and analysis of the image motion velocity (IMV) field holds significant importance in the image motion compensation (IMC) of aerospace remote sensing cameras (ARSCs). Thus, this article puts forward a method for calculating the IMV field based on the rigorous imaging model, which takes into account the camera distortion characteristics. The proposed technique is applied to both a virtual and a real remote sensor to analyze the spatial and temporal features of the IMV field and the influence of camera distortion on it. Our experiments revealed that the additional IMV caused by camera distortion should not be disregarded when calculating IMV field due to its relative magnitude and the decrease in IMC performance caused by it. Additionally, we discussed the selection of the sensor-shift IMC strategy and found that for the real remote sensor considered in this article, the 2-D local compensation is already sufficient to achieve the desired compensation effect.
Index Terms-Aerospace remote sensing cameras, distortion model, image motion compensation, image motion velocity field, sensor-shift.

I. INTRODUCTION
A EROSPACE remote sensing cameras (ARSCs) exhibit characteristics that necessitate consideration of the image motion problem during their usage. Due to the satellite orbital motion, Earth rotation, and satellite attitude adjustment, image motion refers to the movement of the ground target's image, captured via the camera optical system and projected onto the focal plane of the ARSC. Image motion may degrade the quality of remote sensing imagery, in terms of both geometry and radiation, thereby necessitating image motion compensation (IMC) [1], [2], [3].
The most common approach used for IMC is the time delay integration (TDI) technique, [4] whereby images are acquired line-by-line using a linear array sensor, with successive lines generated by the satellite's orbital motion, rendering it a "pushbroom" mode. Alternatively, for imaging process executed by a planar array sensor, another method called "staring" mode requires rotating the camera's optical axis during the exposure time, via the satellite attitude adjustment or fast steering mirror, to maintain alignment with a specific point on the ground, which puts high requirements on the satellite's maneuverability and rotating mechanism's control precision [5], [6], [7], [8]. However, the sensor-shift IMC method remediates the image motion directly on the focal plane by simply moving the sensor, thereby offering a consistent compensation velocity over the entire image plane, with no rotating mechanism demand and avoiding additional optical axis fluctuation. To achieve effective IMC, calculating the magnitude and direction of the image motion velocity (IMV) on the image plane is a crucial first step. Researchers have developed diverse methods for this purpose. By measuring the satellite speed and multiplying it by the scale factor, the speed-to-height ratio model is simple to calculate, but it neglects many practical factors, especially the vertical IMV caused by the deflection angle, and therefore is only suitable for situations where the forward IMV dominates, such as the aerial photography. [9], [10] Considering IMV due to the Earth rotation, Wang et al. [11] used the homogeneous coordinate transformation method to calculate the IMV vector at the subsatellite point, while Yang et al. [12] completed the same calculation using the vector decomposition method. However, they did not consider the dynamic imaging of satellites with certain attitude and only calculated the IMV at the center of the image plane, whereas knowing the IMV distribution over the entire image plane is crucial for cameras with high resolution or wide field of view (FOV). Based on the homogeneous coordinate transformation approach, Li et al. [13] proposed a method considering the case of scroll imaging, while Han et al. [14] combined the idea of ray tracing with vector projection to compute the velocity field on the ground over the Arctic and Antarctic. However, none of them considered the IMV resulting from the satellite attitude maneuver. Therefore, Xu et al. [15] proposed a mathematical model of the IMV field to analyze the effect of circular scanning imaging. So far, all of the aforementioned methods overlooked the fact that Earth is an ellipsoid. Although the methods proposed by Wang et al. [16] and Wu et al. [17] have addressed the aforementioned shortcoming, they have not taken into account the distortion characteristics of ARSCs, potentially introducing errors in IMV calculation.
Therefore, based on the principle of ray tracing, this article proposes a new approach that combines coordinate transformation and vector operations to calculate the IMV field while considering the distortion characteristics of ARSCs. The proposed approach is applied to a virtual remote sensor to examine the spatial and temporal characteristics of the IMV field before it is applied to identify the optimal sensor-shift compensation strategy using the actual remote sensor's parameters. Furthermore, the significance of considering the distortion characteristics of ARSCs is demonstrated by computing the additional IMV field and comparing the IMC results with and without accounting for the camera distortion.

II. METHODOLOGY
This study considers the distortion characteristics of ARSCs by applying the rigorous imaging model. To calculate the IMV field, the connection between the real image point on the image plane and the optical-conjugate point on the ground is established through the transformation into a series of coordinate systems. All calculations are carried out based on vector relations. The proposed workflow, depicted in Fig. 1, entails two primary steps. The first step involves calibrating the camera in the laboratory or on orbit to construct the distortion model. The second step uses the distortion model to compute the IMV field. The subsequent section will outline each specific step in detail.

A. Camera Calibration
Geometric distortion inside the camera, such as lens optical distortion and charge-coupled device (CCD) distortion, alter the geometric properties of the imaging process and must be eliminated to ensure accurate IMV calculations. Mathematical models that describe this internal geometric distortion can be classified into physical and general models. Physical models are established based on the physical characteristics of the various distortions present during camera imaging [18]. While these models clearly explain and define each distortion error source, they also have several disadvantages, such as having too many parameters, and studies have shown that some parameters suffer from strong correlation under long focal length and narrow FOV [19]. Therefore, the general viewing angle model is now commonly adopted for camera calibration of ARSCs, and is also the distortion model used in this article.
The general viewing angle model describes the relationship between the image coordinates (r, c) of the detector unit and its corresponding viewing angle (ϕ x , ϕ y ) [20]. However, to calculate the IMV, the connection between the spatial coordinates of the real and the intended (or the theoretical) image points of an incident light on the image plane is required. Therefore, the conversion of the image coordinates (r, c) in the image system to the spatial coordinates (x, y) in the image plane system needs to be performed.
Suppose that the row and column number commence from the bottom-left corner of the sensor when viewed along the direction of the incident light. Assuming the absence of any internal geometric distortion in the camera, the incident light emanating from the ground point o ought to be imaged at the detector unit (r, c), with the theoretical image plane coordinates (x it , y it ) being where a p denotes the pixel size and (r 0 , c 0 ) are the image coordinates of the principal point of photography. However, because of the internal geometric distortion, the incident light is realistically imaged at (x ir , y ir ), as described in the general viewing angle model From (1) and (2), we can obtain where and Equation (3) is the distortion formula that incorporates the general viewing angle model. On-orbit or in-laboratory camera calibration provides the parameters in the general viewing angle model, which reveal the spatial relationship between the theoretical and real image points [21].

B. Calculation of IMV Field
To avoid any confusion, it is essential to initially define some of the symbols used in this article: S P -The image plane coordinate system: the origin P is the intersection of the camera's optical axis and the focal plane, which is the principal point of photograph, X P is along the flight direction, Y P is perpendicular to the flight direction. If there are no interior errors in the camera, P should be centered at the sensor, X P should be along the column direction of the pixel arrangement, and Y P should be opposite to the row direction. S C -The camera coordinate system: the origin C is the rear nodal point of the objective, Z C is the optical axis pointing to Earth, the X C CY C plane is perpendicular to the optical axis and X P and X C (Y P and Y C ) are parallel to each other. S S -The satellite coordinate system: the origin S is at the satellite centroid. S O -The orbit coordinate system: the origin O is on orbit, Z O points to the geocentric, X O is along the orbital forward direction, and they are both in the orbit plane. The three-axis attitude angle of the satellite is the three-axis attitude angle of the S-system relative to the O-system. S E -The Earth coordinate system: the origin E is the geocentric, X E is in the equator plane pointing to the vernal equinox, and Z E points to the north pole. × -The origin of × system, e.g., point E is the origin of the Earth coordinate system. o -The object point on the ground. i t -The theoretical image point, which is the location where the light imaged at the actual image point should theoretically be imaged in the absence of camera distortion. i r -The real image point, which is the point on the image plane where we want to calculate the IMV.
The coordinates of the vector from point j to point k in the L-system, e.g., V C Co is the coordinates of the vector from the origin of S C to the object point o, which is also the . T JK -The rotation transformation matrix from K-system to J-system, e.g., the camera's mounting matrix on the satellite T CS is the rotation transformation matrix from S-system to C-system. X -The derivative of X with respect to time. The process of calculating IMV comprises of four sequential steps, which are as follows.
The first step involves the selection of interested point and the correction of sight vector. We opt for a point on the image plane, where IMV needs to be calculated and then apply the distortion model that was obtained earlier to convert from the real image point i r to the corresponding theoretical image point i t . Consider the coordinates of i r in S P to be (x P ir , y P ir ), and the coordinates of its corresponding i t in S P can be calculated through (3).
In the second step, we solve for the position of the corresponding ground point according to the spatial relationship in the instantaneous rigorous imaging model. Assuming that the object point o, optical-conjugate to i t , has coordinates V C Suppose the coordinates of o in , and using the WGS84 Earth ellipsoid model, we have Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
To solve for V C Co and V E Eo , wherein there are six unknowns but three equations only, we need to find the connection between them to obtain more equations. V C Co and V E Eo refer to the same point in different coordinate systems, addressing the issue may be resolved by the transformation from S E to S C , which is carried out as follows.
1) Transformation From S E to S O : To obtain the O-system, we can turn the E-system around the Z axis by Ω, then around the X axis by (i 0 − π/2), then around the Y axis by −(π/2 + γ(t) + ω), and finally translate the E-system along the Z axis by −l EO (t). Therefore, we have and and Once the camera is mounted on the satellite, its mounting matrix T CS and the position of the camera photography center on the satellite V S SC can be obtained in the laboratory. Then, we have By incorporating (9), (12), and (13), we will find Simultaneously solving (6), (7), and (14), we can get V C Co and V E Eo . It is important to note that there may be two sets of solutions, which correspond to the two intersections of the camera's view axis and the Earth's surface. Therefore, as per the physical interpretation, the solution related to the point at the back of the Earth should be disregarded.
In the third step, the object point is backprojected onto the image plane in the dynamic rigorous imaging model, and the time-varying relationship of involved physical quantities are used to further determine the position-time relationship of the theoretical image point. According to the assumptions in this article, the P -system can be obtained by translating the C-system along the Z axis by −f . Combining this with the collinearity condition equations, we have Differentiate both sides of the (15) with respect to time t, we haveV To resolve the (16), we differentiate both sides of the (14) with respect to t, givinġ Typically, after the camera is mounted on the satellite, the mounting matrix and position of the photography center remain constant (i.e.,Ṫ CS = 0;V S SC = 0). Thus, we simplify (17) aṡ ( 2 4 ) By solving (16)-(24) together, we can retrieve the position-time relationship of the theoretical image point, i.e.,V P P it . In the final step, we solve for the position-time relationship of the real image point and calculate the IMV there. Since our goal is to acquire the IMV at the real image point, we need to use the distortion model again to convert the theoretical image point to the real image point. The derivative of (3) with respect to time t can be obtained as By substituting V P P it andV P P it into (25), we can obtain the position-time relationship of the real image point (i.e.,V P P ir ), thus enabling us to calculate the IMV and deflection angle ⎧ ⎨ To sum up, in the second and third steps, we combine the coordinate transformation and vector operation methods to determine the position-time relationship of the image point based on the rigorous imaging model of ARSCs. In the first and final steps, we strive to consider the camera distortion to acquire a more precise IMV output. We proceed by selecting several other points on the image plane and repeating the above steps to calculate the IMV and deflection angle at those points. Finally, we obtain the IMV field distributed on the whole image plane.

III. EXPERIMENTS
We examine the spatial and temporal characteristics of the IMV field and investigate the impact of camera distortion by analyzing the calculation results on the virtual remote sensor. Subsequently, we discuss the feasibility of the sensor-shift IMC method by applying the approach proposed in this study to a real remote sensor. Table I displays the parameters of the virtual remote sensor, which are similar to those of the real remote sensor, albeit with slight variations in specific values.   satellite attitude has a more significant effect on the IMV field when compared to the argument of latitude and it is worth noting that the tangential distortion added to the virtual remote sensor results in the nonaxisymmetric distribution of IMV field on the image plane when the camera points to the subsatellite. The four points at the edge of the image plane [numbered as (1,−1), (1,1), (−1,−1), (−1,1)] and the center point [numbered as (0,0)] are selected to calculate the variation curves of the IMV within one track at these positions, and the results are shown in Fig. 3(a). Regular IMV fluctuations are observed at each point, and there are small deviations in the argument of latitude when the IMV reaches the extreme value at different positions. This is mainly caused by the rotation of the Earth, so that the IMV of each point on the image plane reaches a maximum when the corresponding object point passes through the equator, and reaches a minimum when it passes through southernmost and northernmost ends. The four points at the edge of the image plane and the center point are selected to calculate the variation curves of the IMV with roll and pitch angle, and the results are shown in Fig. 3(b) and (c). At these positions, the IMV decreases with an increase in the angle of satellite attitude at the imaging moment. Moreover, the deviation of the IMV between the edge and center of the image plane also decreases. Fig. 4 illustrates the additional IMV field that results from camera distortions, which is calculated by taking the difference of the IMV field results before and after considering the camera distortion. The length of the arrow in the figure shows the absolute value of the additional IMV, while the colors indicate the relative value of the additional IMV compared to the IMV without considering the camera distortion. The results show that the additional IMV at the edge of the image plane is more significant in both absolute and relative levels, which is consistent with the characteristics of radial distortion. The absolute value of the additional IMV due to camera distortion is more substantial when the camera is pointing at the subsatellite, with a maximum of 0.98 mm/s. However, the relative value of the additional IMV is larger at specific satellite attitude, where the maximum can be up to 8.5%.

C. Sensor-Shift Image Motion Compensation Strategy
Now consider the case of a real remote sensor in order to discuss how sensor-shift IMC can be applied. The imagery quality degradation due to image motion can be measured by the image motion modulation transfer function (MTF): where f c denotes the Nyquist frequency and s is the magnitude of image motion. Empirically, MTF greater than 0.95 after compensation is an acceptable result, corresponding to an image motion of 1/3 pixel size. From the previous analysis, we know that the IMV on the image plane reaches its maximum when the subsatellite passes through the equator and the camera points directly on the subsatellite. Therefore, in the next discussion, the imaging conditions are fixed as mentioned above.
In order to implement sensor-shift IMC, high-precision shift stages must be mounted at the focal plane to drive the sensor to move at the same speed as the image during the exposure time. The compensation process primarily considers the IMV along the track induced by the satellite's orbital motion, which is much more substantial than the IMV across the track due to Earth's rotation, but it is also pertinent to carefully analyze whether or not to compensate for the IMV across the track. Employing the calculation result from just one point on the image plane, usually the center point, as the input of IMC can markedly decrease the computational complexity. However, as we have previously noted, the IMV's magnitude and direction varies across the image plane, and if the compensation only targets one point, the compensation residual would affect the imagery quality of other positions. Hence, it is crucial to know whether the IMV field over the entire image plane should to be considered for the compensation process.
Using the parameters of the real remote sensor, the IMV field and MTF distribution on the image plane are calculated. The calculations are conducted on three scenarios: before IMC, after IMC considering the IMV along the track at the center point of image plane (1-D local compensation), and after IMC considering the IMV along the track over the entire image plane (1-D global compensation). The results are depicted in Fig. 5(b) and (c), and compensation residuals are presented in Table II. Because the IMV along the track accounts for a large proportion in the whole combined IMV, the results show a significant  improvement in MTF in most areas of the image plane after 1-D compensation, but in some areas, the decline in MTF is still noticeable. In addition, while compensating for the entire image plane can enhance MTF at those areas, there is still a drop in MTF by 10-13% after compensation, which exceeds the tolerable limit of 5%. Thus, compensating image motion in only one direction (i.e., along the track) is insufficient, and it is necessary to compensate for the complete combined velocity vector.
The compensated MTF has been significantly enhanced by orthogonally decomposing the IMV in the along-track and across-track directions and compensating for them, respectively, as depicted in Fig. 5(e) and (f), and compensation residuals are presented in Table II. Even only considering the IMV at the center point of image plane (2-D local compensation), the compensated MTF has reached above 0.95 over the entire image plane. Although a slight improvement in MTF has been observed after considering the entire image plane (2-D global compensation), it is no longer advantageous over the 1-D compensation. Hence, the optimal strategy is the 2-D compensation considering only the IMV at the center point of image plane.

D. Effect of Camera Distortions on IMC Performance
The additional IMV field due to camera distortion of the real remote sensor and its effect on IMC performance are illustrated in Fig. 6. Negligible camera distortion from the real remote sensor result in small additional IMV, as seen in Fig. 6(a) and (b), where the maximum additional IMV is 0.76 mm/s and the maximum relative additional IMV is only 3.5%.
MTF results after IMC with and without considering the distortion characteristics are presented in Fig. 6(d) and (c), respectively. The relative improvement in MTF after considering camera distortion is depicted in Fig. 6(e), with the green highlighted part indicating the area where MTF increases from below 0.95 to above 0.95. The camera distortion of the real remote sensor have been minimally controlled, resulting in the small MTF improvement. Nonetheless, the presence of distortions still has a significant effect on IMC performance at the edge of the image plane, and it is still necessary to take the camera distortion into consideration.

IV. DISCUSSION
In this article, we propose a new method to calculate the IMV field, which takes into consideration the distortion of ARSCs. We further apply our methodology on a virtual remote sensor to analyze the spatial and temporal characteristics of the IMV field and the effect of camera distortion on it. The results reveal that the distribution of the IMV field at any given instance is nonlinear, owing to the ellipsoidal shape of the Earth and the camera distortion. Moreover, we show that different imaging poses or positions leads to variations in the magnitude of IMV Fig. 6. Additional IMV field due to distortion (a), (b); MTF distribution after 2-D global compensation with (d) and without (c) considering the camera distortion; relative MTF enhancement after considering the camera distortion(e) @ at-nadir viewing and over the ascending node and no attitude jitters. field and the IMV will increase as satellites get closer to the equator or have smaller attitude angles. In addition, the distribution of the additional IMV field caused by the camera distortion on the image plane is consistent with the distortion characteristics, i.e., the additional IMV at the edge of the image plane is larger than that at the center of the image plane, whether in absolute or relative level. Furthermore, we discuss the selection of sensor-shift IMC strategy of a real ARSC under development, and find that 2-D local compensation can significantly enhance the quality of remote sensing imagery, without the necessity of extensive computational resources. Finally, we quantify the impact of camera distortion, demonstrating that distortion at the edge of the image plane has a significant impact on IMV performance and should not be ignored during the calculation of IMV.
We must mention that the well-controlled distortion of the real remote sensor used in this study makes the calculated additional IMV field and the subsequent IMC performance enhancement from considering it less substantial. Nonetheless, in the case of remote sensors with more evident distortion, a comprehensive IMV calculation that takes into account the distortion characteristics of the ARSC is critical. Besides, our study has selected the IMC strategy solely based on the outcomes of the theoretical case. We have to acknowledge that an on-orbit experiment can further demonstrate the actual compensation performance of IMC, due to the inclusion of multiple practical factors, such as the scarcity of various resources that the microsatellite faced and the precision limit of the shift stages. Moreover, according to the existing literature, the accuracy of the IMV field calculations is difficult to verify precisely. Thus, we hope that a well-designed experiment combined with some new technique may address these concerns as a topic for further research after the satellite is launched.