Optical Flow Inversion for Remote Sensing Image Dense Registration and Sensor’s Attitude Motion High-Accurate Measurement

It has been discovered that image motions and optical flows usually become much more nonlinear and anisotropic in space-borne cameras with large field of view, especially when perturbations or jitters exist.The phenomenon arises from the fact that the attitude motion greatly affects the image of the three-dimensional planet. In this paper, utilizing the characteristics, an optical flow inversion method is proposed to treat high-accurate remote sensor attitude motion measurement. The principle of the new method is that angular velocities can be measured precisely by means of rebuilding some nonuniform optical flows. Firstly, to determine the relative displacements and deformations between the overlapped images captured by different detectors is the primary process of the method. A novel dense subpixel image registration approach is developed towards this goal. Based on that, optical flow can be rebuilt and high-accurate attitude measurements are successfully fulfilled. In the experiment, a remote sensor and its original photographs are investigated, and the results validate that the method is highly reliable and highly accurate in a broad frequency band.


Introduction
For the remote sensors in dynamic imaging, one important technology is image motion compensation.Actually, to determine image motion velocity precisely is a very hard problem.In [1,2], optical correlators are utilized to measure image motion in real time based on a sequence of mild smeared images with low exposure.This technique is appropriate to the situations in which the whole image velocity field is uniform.Some other blind motion estimation algorithms in [3][4][5] have been used to image postprocessing, which can roughly detect inhomogeneous image motion, but lack real-time performance because of complexity.As for space imaging, in order to avoid motion blurring, image motion velocity needs to be computed in real time according to the current physical information about spacecraft's orbit and attitude motion, which can be obtained by the spaceborne sensors, such as star trackers, gyroscopes, and GPS.
Wang et al. developed a computational model for image motion vectors and presented error budget analysis in [6].They focused on the small field of view (FOV) space cameras which are used in push-broom imaging with small attitude angles.In that situation, the nonlinearity of image motion velocity field does not appear significantly.However, for others with larger FOV, image motion velocity fields are definitely nonlinear and anisotropic because the geometry of the planet will greatly modulate the moving images.Under the circumstances, the detectors need to be controlled separately to keep the time series synchronized with the instantaneous image velocities.
The time-phase relations between the photos belonging to different detectors are affected by optical flows, which are uniquely determined by the behavior of image velocity field in a specific period.Some phenomena about moving image variation and distortion due to optical flow have been discovered [7][8][9][10].References [7,8] reported the camera on 2 Mathematical Problems in Engineering Mars Reconnaissance Orbiter (MRO) in High Resolution Imaging Science Experiment (HiRISE) missions of NASA.It takes pictures of Mars with resolutions of 0.3 m/pixel resolving objects.Fourteen staggered parallel CCDs are overlapped with 48 pixels at each end to fulfill the entire field of view.Although adjacent detectors overlap with equal physical pixels, yet their lapped image pixels are not equal and varying with time, because spacecraft jitters cause undulating optical flows within the interlaced areas [8].In addition, we found that when large FOV remote sensors perform stereoscopic imaging with large pitch angles, the lapped images belonging to marginal detectors are bound to exceed or lose several hundred pixels compared to their physical overlaps.Furthermore, the unexpected quantity decreases significantly for the detectors mounted at the central region of the focal plane.
Although nonuniform optical flow brings many troubles in image processing such as registration, resample, interconnection, and geometrical rectification, it permits us to measure the spacecraft attitude motion with very high accuracy in a broad bandwidth.It is nearly impossible for conventional space-borne sensors to realize the target.Precision attitude motion measurement is very useful for remote sensing image processing, especially for image restoration from motion blurring as studied in [11,12].Associating the measurement and optical flow models, the dynamic point spread functions (PSF) are able to be estimated to be set as the convolution kernels in nonblinded deconvolution algorithms.
The behavior of optical flow characterizes the entire twodimensional flow field for an image's motion and variation.In [13] optical flow estimation based on image sequences of the same aurora to determine the flow field will provide access to the phase space, the important information for understanding the physical mechanism of the aurora.For the purpose to improve the accuracy of optical flow estimation, a two-step matching paradigm for optical flow estimation is applied in [14]; firstly the coarse distribution measuring of motion vectors is done with a simple frame-to-frame correlation technique also known as the digital symmetric phaseonly filter (SPOF), and after that subpixel accuracy estimation result is achieved by using the sequential tree-reweighted max-product message passing (TRW-S) optimization.Similarly, Sakaino overcame the disadvantages in optical flow determination when moving objects with different shapes and sizes move against a complicated background; the image intensity between frames may violate the common situation image brightness constancy to image brightness 5 change models as constraints in regular situations [15].However, unlike continuous image sequences, if we merely obtained several images of the identical moving objects captured by different detectors with long intervals, the former techniques do not work well for optical flow estimation for lacking of the information of imaging process of the instrument.
In this paper, a new optical flow inversion method is proposed for precise attitude measurement.Unlike the situations in [13][14][15], the image sequences of video do not exist for the transmission type remote sensors, instead of the image pairs of the same earth scene which are captured by different TDI CCD detectors in push-broom fashion.The time intervals between the independent image formations corresponding to the overlapped detectors are much more than the interval between sequential frames in video, for which, the frame rates usually exceed tens of frames per second (fps).However, we can model optical flows based on the working mechanism of the instrument and image processing techniques rather than estimating from frame sequences of a specific detector.The contents of this paper are organized as follows: in Section 2, an analytical model of image motion velocity field is established, which is applicable to dynamic imaging for three-dimensional planet surface by large FOV remote sensors.The phenomenon of moving image deformation due to optical flow is investigated in Section 3. Based on rough inversion of optical flow, a novel method for dense image registration is developed to measure the subpixel offsets between the lapped images which are captured by adjacent detectors.In Section 4, an attitude motion measuring method based on precise optical flow inversion is studied, and the results of the experiment support the whole theory perfectly.

Image Velocity Field Analysis
Suppose that a large FOV camera is performing pushbroom imaging to the earth; the scenario is illustrated in Figure 1.The planet's surface cannot be regarded as a local plane but a three-dimensional ellipsoid since it may greatly influence the image motion and time-varying deformation when complicated relative motion exists between the imager and the earth.
In order to set up the model of space imaging, some coordinate systems need to be defined as follows.
(1) I(  − ): the inertial frame of the earth.For convenience, here we choose 2000 frame.The origin   is located at the earth center.
(2) C(  −   1   2   3 ): the frame of camera.Axis     3 is the optical axis, and origin   is the center of the exit pupil.
(3) O(  −  1  2  3 ): the orbit frame.Axis    3 passes through the center of earth, and axis    2 is perpendicular to the instant orbit plane.
(5) P( − ): the frame of photography.The origin  is the center of the photo.Axis  points to the column direction, and axis  points to the row direction.imager holds a large attitude angle.Obviously, the shapes and lengths of  1 and  2 also have notable differences during push broom, which implies that the geometrical structure of the image is time varying as well as nonuniform.Furthermore, it can be discovered later that the deforming rates mainly depend on the planet's apparent motion observed by the camera.
Considering an object point  on the earth, its position vector relative to   is denoted as ⃗   .As a convention in the following discussions, I ⃗   represents the vector measured in frame I, and accordingly C ⃗   is the same vector measured in frame C. We select one unit vector ⃗  which is tangent to the surface of the earth at .Let ⃗ ( 1 ,  2 ,  3 ) be the position vector of  relative to   ; then C ⃗  and C ̇⃗  characterize the apparent motion of .Assume that the image point   is formed on the focal plane with coordinates (  1 ,   2 ,   3 ) in frame C. Generally, the optical systems of space cameras are well designed and are free from optical aberrations and the static PSF is approximate to the diffraction limit [16,17]; thus following [18], we have where   is the effective focal length, the lateral magnification of   :  = (−1) −1 (  /( ⃗  ⋅ ⃗  3 )),  is the number of intermediate images in the optical system, and ⃗   ( = 1, 2, 3) is the base of C.
Let ⃗   be the position vector of satellite relative to   ; then ⃗  = ⃗  − ⃗   .In imaging, the flight trajectory of the satellite platform in I can be treated as Keplerian orbit, as illustrated in Figure 2. According to the orbit elements,  0 , inclination; Ω, longitude of ascending node; , argument of perigee; , semimajor axis; , eccentricity;   , mean anomaly at epoch, we implement Newton-Raphson method to solve (2) and get the eccentric anomaly  from the given mean anomaly   =  0 + ( −  0 ), where  = 2/,  is the orbit period [11]: ( Orbit plane The equatorial plane Perigee In frame O, The coordinate transform matrix between O and I is For simplicity, we write  := cos ,  := sin .
In engineering, the coordinate transfer matrix T OI also can be derived from the real-time measurements of GPS.Since the base vectors of frame O in I: Associating the equation of boresight with the ellipsoid surface of the earth in C yields Here,   = 6378.137km and   = 6356.752km being the length of earth's semimajor axis and semiminor axis.  ( = 1, 2, 3) are the unit vectors of I ⃗ .We write the solution of (7) as I ⃗  = (  )  .Hence, where M is the coordinate transformation matrix from frame B to frame C; it is a constant matrix for fixed installation; A is the attitude matrix of satellite; according to 1-2-3 rotating order, we have in which where   ,   , and   are in order, the real-time roll angle, pitch angle, and yaw angle at moment .The velocity of  in C can be written in the following scalar form: Thus, the velocity of image point of   will be Substituting ( 2)-( 9) into (10), the velocity vector of image point ⃗ can be expressed as the explicit function of several variables; that is, For conciseness, this analytical expression of ⃗ V  is omitted here.
The orbit elements can be determined according to instantaneous GPS data.Besides, they also can be calculated with sufficient accuracy in celestial mechanics [19].On the other hand, the attitude angles   ,   , and   can be roughly measured by the star trackers and GPS.Meanwhile, their time rates φ  , θ  , and ψ  have the following relations: 1 ,  2 , and  3 are the three components of the remote sensor's angular velocity C ⃗   relative to orbital frame O, which is calibrated in frame C.Those can be roughly measured by space-borne gyroscopes or other attitude sensors.
It is easy to verify from ( 11) that the instantaneous image velocity field on the focal plane appears significantly nonlinear and isotropic for large FOV remote sensors, especially when they are applied to perform large angle attitude maneuvering, for example, in sidelooking by swing or stereoscopic looking by pitching and so forth.Under these circumstances, in order to acquire photos with high spatial, temporal, and spectral resolution, image motion velocity control strategies should be executed in real time [20] based on auxiliary data which measured by reliable space-borne sensors [21,22].In detail, for TDI CCD cameras, the line rates of the detectors must be controlled synchronizing to the local image velocity modules during exposure so as to avoid along-track motion blurring; the attitude of remote sensor should be regulated in time to maintain the detectors push-broom direction aiming at the direction of image motion to avoid cross-track motion blurring.

Optical Flow Rough Inversion and Dense Image Registration
Optical flow is another important physical model carrying the whole energy and information of moving images in dynamic imaging.A specific optical flow trajectory is an integral curve which is always tangent to the image velocity field; thus we have Since ( 13) are coupled nonlinear integral equations, we convert them to numerical forms and solve them iteratively, It is evident that the algorithm has enough precision so long as the step-size of time interval Δ is small enough.It can be inferred from ( 13) that strong nonlinear image velocity field may distort optical flows so much that the geometrical structure of image may have irregular behaviors.Therefore, if we intend to inverse the information of optical flow to measure the attitude motion, the general formula of image deformation due to the optical flows should be deduced.

Time-Varying Image Deformation in Dynamic Imaging.
Firstly, we will investigate some differential characteristics of the moving image of an extended object on the earth surface.As shown in Figure 1, considering a microspatial variation of  along ⃗  on the curved surface can be expressed as  ⃗   =  ⃗ .Its conjugated image is We expand the term of : Taking derivatives with respect to variable  for either part of (15), we have According to (16) we know that  β ≈ 0. On the other hand, the variation of ⃗  can be expressed through a series of coordinate transformations; that is Notice that E ⃗  is a fixed tangent vector of earth surface at object point , which is time-invariant and specifies an orientation of motionless scene on the earth. Consequently, where the coordinate transform matrix from frame E to I is Let   be the angular rate of the earth and   the longitude of  on the earth; then the hour angle of  at time  is   () = GST +   +   , in which, GST represents Greenwich sidereal time.
The microscale image deformation of the extended scene on the earth along the direction of ⃗  during  1 ∼  2 can be written as From ( 17), we have According to ( 16), (18), and ( 19) we obtain the terms in (22): Furthermore, if the camera is fixed to the satellite platform, then Ṁ = 0, ̇⃗   = 0. Consequently, (22) becomes For the motionless scene on the earth surface, E ⃗  is a timeindependent but space-dependent unit tangent vector, which meanwhile represents a specific orientation on the ground.Moreover, the physical meaning of function F  (, ⃗ ) is the image deformation of unit-length curve on the curved surface along the direction of E ⃗  in unit time interval.That is the instantaneous space-time deforming rate of the image of the object along E ⃗ .Consequently, in dynamic imaging, macroscopic deformation on the moving image can be derived from the integral of F  (, ⃗ ) in space and time.Referring to Figure 1, let Γ be an arbitrary curve of the extended object on the earth, let Γ  be its image, let two arbitrary points ,  ∈ Γ, and let their Gaussian images   ,   ∈ Γ  .Let E ⃗  = ⃗ T() be a vector-valued function with variable  (the length of the arc) which is time-invariant in frame E and gives the tangent vectors along the curve.
So, the image deformation taking place during  1 ∼  2 is able to be described as in which ]. Now, in terms of ( 24) and ( 25), we can see that the image deformation is also anisotropic and nonlinear which depends not only on optical flow's evolution but also on the geometry of the scene.

Dense Image Registration through Optical Flow Prediction.
As mentioned in the preceding sections, optical flow is the most precise model in describing image motion and timevarying deformation.On the contrary, it is possible to inverse optical flow with high accuracy if the image motion and deformation can be detected.As we know, the low frequency signal components of angular velocity are easier to be sensed precisely by attitude sensors such as gyroscopes and star trackers, but the higher frequency components are hard to be measured with high accuracy.However, actually, perturbations from high frequency jittering are the critical reason for motion blurring and local image deformations, since the influences brought by low components of attitude motion are easier to be restrained in imaging through regulating remote sensors.
Since ( 13) and ( 25) are very sensitive to the attitude motion, the angular velocity is able to be measured with high resolution as well as broad frequency bandwidth so long as the image motion and deformation are to be determined with a certain precision.Fortunately, the lapped images of the overlapped detectors meet the needs, because they were captured in turn as the same parts of the optical flow pass through these adjacent detectors sequentially.Without losing generality, we will investigate the most common form of CCD layout, for which two rows of detectors are arranged in parallel.The time-phase relations of image formation due to optical flow evolution are illustrated in Figure 3, where the moving image elements { 1 ,  2 , . ..} (in the left gap), { 1 ,  2 , . ..} (in the right gap) are captured firstly at the same time since their optical flows pass through the prior detectors.However, because of nonuniform optical flows, they will not be captured simultaneously by the posterior detectors.Therefore, the geometrical structures of photographs will be time varying and nonlinear.It is evident from Figure 3 that the displacements and relative deformations in frame C between the lapped images can be determined by measuring the offsets of the sample image element pairs in frame P.
Let Δ  = Δ  1 , Δ  = Δ  2 be the relative offsets of the same object's image on the two photos; they are all calibrated in C or F. We will measure them by image registration.
As far as image registration method is concerned, one of the hardest problems is complex deformation, which is prone to weaken the similarity between the referenced images and sensed images so that it might introduce large deviations from the true values or even lead to algorithm failure.Some typical methods have been studied in [23][24][25].Generally, most of them concentrated on several simple deforming forms such as affine, shear, translation, rotation, or their combinations instead of investigating more sophisticated dynamic deforming models.In [26][27][28][29][30], some effective approaches have been proposed to increase the accuracy and the robust of algorithms according to the respective reasonable models according to the specific properties of objective images.
For conventional template based registration methods, once a template has been extracted from the referenced image, the information about gray values, shape, and frequency spectrum does not increase since no additional physical information resources would be offered.But actually, such information has changed when the optical flows arrive at the posterior detectors.Therefore, the cross-correlations between the templates and sensed images certainly reduce.So, in order to detect the minor image motions and complex deformations between the lapped images, high-accurate registration is indispensable, which means that more precise model should be implemented.We treat it using the technique called template reconfiguration.In summary, the method is established on the idea of keeping the completion of the information about optical flows.In operating, as indicated in Figure 3, take the lapped images captured by the detectors in prior array as the referenced images and the images captured by posterior detectors as the sensed images.Firstly we will rebuild the optical flows based on the rough measurements of the spaceborne sensors and then reconfigure the original templates to construct the new templates whose morphologies are more approximate to the corresponding parts on the sensed images.With this process, the information about imaging procedures is able to be added into the new templates so as to increase the degree of similarity to the sensed images.The method may dramatically raise the accuracy of dense registration such that the high-accurate offsets between the lapped image pairs are able to be determined.
In the experiment, we examined Mapping Satellite-1, a Chinese surveying satellite operating in 500 km height sun synchronous orbit which is used for high-accurate photogrammetry [31] whose structure is shown in Figure 4.One of the effective payload, three-line-array panchromatic CCD cameras has good geometrical accuracy whose ground pixel resolution is superior to 5 m, spectral range is 0.51 m ∼ 0.69 m, and the swath is 60 km.Another payload is that the high resolution camera is designed possessing Cook-TMA optical system which gives a wide field of view [16,17] and the panchromatic spatial resolution can reach 2 m.
In engineering, for the purpose to improve the image quality and surveying precision, the high-accuracy measurements of jitter and attitude motion are very essential for posterior processing.Thus, here we investigate the images and the auxiliary data of the large FOV high resolution camera to deal with the problem.The experimental photographs were captured with 10 ∘ side looking.The focal plane of the camera consists of 8 panchromatic TDI CCD detectors and there are  = 96 physical lapped pixels between each other.
The scheme of the processing in registering one image element  is illustrated in Figure 5.
Step 1. Set the original lapped image strips (the images which were acquired directly by the detectors and without any postprocessing) in frame C.
Step 2. Compute the deformations of all image elements on referenced template with respect to their optical flow trajectories.
We extract the original template from the referenced image denoted as  1 , which consists of  2 square elements; that is, dim( 1 ) =  × .Let  be its central element and  the width of each element; here  = 8.75 m.Before the moving image was going to be captured by the posterior detector, in terms of (25), their current shapes and energy distribution can be predicted by the optical flow based on the auxiliary data of the remote sensor.
In order to simplify the algorithm, first order approximation is allowed without introducing significant errors.This approximation means that the shape of every image element is always quadrilateral.Linear interpolations are carried out to determine the four sides according to the deformations along the radial directions of the vertexes, as showed in Figure 5.The unit radial vectors are denoted by ⃗   1 ∼ ⃗   4 in frame C: Suppose image point   is the center of an arbitrary element Σ  in  1 .Let Σ be the area element on the earth surface which is conjugate to Σ  .The four unit radial vectors of the vertexes on Σ, ⃗  1 ∼ ⃗  4 are conjugate to ⃗   1 ∼ ⃗   4 and tangent to the earth surface at .From the geometrical relations, we have where E ⃗   is the unit normal vector of Σ at .We predict the deformations along ⃗  1 ∼ ⃗  4 during  1 ∼  2 according to the measurements of GPS, star trackers, and gyroscopes, as explained in Figure 6. 1 is the imaging time on prior detector and  2 is the imaging time on the posterior detector, The shape of deformed image Σ  2 can be got through linear interpolation with Step 3. Reconfigure referenced template  1 according to optical flow prediction, and then get a new template  2 .
Let   1 be the deformed image of  1 computed in Step 2. Let  =  , be the central element of   1 ; integers  and  are, respectively, the row number and column number of  , .The gray value   of each element in    1 is equal to its counterpart in  1 with the same indexes.In addition, we initialize a null template  0 whose shape and orientation are identical to  1 ; the central element of  0 is denoted by  , .Then, we cover  0 upon   1 and let their centers coincide; that is,  , =  , , as shown in Figure 7. Denote the vertexes of   1 as    ( = 1 ∼ 4).Therefore, the connective relation for adjacent elements can be expressed by  1  , =  2 ,−1 =  3  −1,−1 =  4 −1, .Next, we will reassign the gray value ℎ   to  , ( = 1 ⋅ ⋅ ⋅ ,  = 1 ⋅ ⋅ ⋅ ) in sequence to construct a new template  2 .The process is just a simulation of image resample when optical flow arrives at the posterior detector, as indicated in Figure 3.
That is, Weight coefficient   =   / 2 where   is the area of the intersecting polygon of  , with  , .Step 4. Computenormalized cross-correlation coefficients between  2 and the sensed image, and then determine the subpixel offset of  2 relative to the sensed image in frame P. Firstly, for this method, the search space on the sensed image can be contracted so much since the optical flow trajectories for the referenced elements have been predicted in Step 2 Assuming that the search space is   , dim(  ) =  × .When  , moves to the pixel ( 1 ,  2 ) on   , the normalized cross-correlation (NCC) coefficient is given by where  , is the mean gray value of the segment of   that is masked by  2 and ℎ is the mean of  2 .Equation ( 31) requires approximately  2 ( −  + 1) 2 additions and  2 ( −  + 1) 2 multiplications, whereas the complexity of FFT algorithm needs about 12 2 log 2  real multiplications and 18 2 log 2  real additions/subtractions [32,33].
At the beginning, we take  = 101,  = 7 and compute the NCC coefficient.When  is much larger than , the calculation in spatial domain will be efficient.Suppose that the peak value γmax is taken at the coordinate (, ), ,  ∈ Z in the sensed window.Hence we will reduce search space into a smaller one with dimension of 47 × 47 which centered on   (, ).Next, the subpixel registration is realized by phase correlation algorithm with larger  and  to suppress the system errors owing to the deficiencies of detailed textures on the photo.Here we take  = 47,  = 23.Let the subpixel offset between the two registering image elements be denoted as   and   in frame P.
The phase correlation algorithm in the frequency domain becomes more efficient as  approaches  and both have larger scales [28].Moreover, the Fourier coefficients are normalized to unit magnitude prior to computing the correlation so that the correlation is based only on phase information and being insensitive to changes in image intensity [27,29].
Let G(, V) be the 2D Discrete Fourier Transforms (DFT) of the sensed window; then we have Here Cross-phase spectrum is given by where H * is the complex conjugate of H.By inverse Discrete Fourier Transform (IDFT) we have (36) The right side presents the spatial distribution of the normalized cross-correlation coefficients.Therefore, (  ,   ) are able to be measured based on that.In practice, constant  ≤ 1, which tends to decrease when small noise exists and equals unity in ideal cases.
Step 5. Dense registration is executed for the lapped image strips.
Repeating Step 1∼Step 4, we register the along-track sample images selected from the referenced images to the sensed image.The maximal sample rate can reach up to line-by-line.The continuous procedure is shown in Figure 8, in which, the image pairs are marked.
The curves of relative offsets in P are shown in Figures 9  and 10.
Let col  , row  be the column and row indexes of image elements on the referenced image, and let col  , row  be the indexes of the same elements on the sensed image.The total columns of each detector  = 4096 pix, and the vertical distance between the two detector arrays  = 18.4975 mm.According to the results of registration, we get the offsets of images at th gap,   (cross track),   (along track) in frame P, and Δ   , Δ   (mm) in frame F: Four pixels S11, S12, S31, and S32 are examined.Their data are listed in Table 1.S11 and S31 are the images of the same object which was captured in order by CCD1 and CCD2 (Gap 1).S12 and S32 were captured, respectively, by CCD3 and CCD4 (Gap 3).Referring to the auxiliary data, S11 and S31 were captured at same time and S12 and S32 were captured at different time, which means that the along-track speeds of the two moving images were quite different.Moreover, the crosstrack image offsets in Gap 1 and Gap 3 vary so much, which says that the optical flows were also distorted unevenly and deflects away from the along-track direction.On the other

Remote Sensor Attitude Motion Measurement
In this section, the attitude velocity of the remote sensor is going to resolved by using optical flow inversion method.The results of dense registration are applied to produce conditions of fixed solution for optical flow equations.

The Principle of Optical Inversion.
For clarity, in frame C, the two coordinate components of image displacement of th sample element belonging to th lapped strip pair are written as Δ   , Δ   .From ( 13) and (25), it is easy to show that the contributions to optical flow owing to orbital motion and earth's inertial movement are of very slightly varying in short term such that the corresponding displacements can be regarded as piecewise constants   ,   .
Let  , ,  , be in order the two sequential imaging time of the th image sample on the overlapped detectors in th gap.They are usually recorded in the auxiliary data of the remote sensor.Hence, for every image element, the quantity of discrete status in optical flow tracing will be where  is the amount of CCD gaps,  is the amount of sample groups, and Δ is the time step.We set samples with same  index into the same group, in which the samples are captured by the prior detectors simultaneously.We expand (11), substitute it into ( 14) and ( 13), and then arrange the scalar optical flow inversion equations in terms of the three axial angular velocity components  1 ,  2 , and  3 (the variables in the inverse problem), yielding the linear optical flow equations.For the th group samples, Suppose that the sample process will stop until  groups have been founded.The coefficients are as follows:

Mathematical Problems in Engineering
Here (41) As for the algorithm, to reduce the complexity, all possible values for the coefficients are stored in the matrixes Ξ  and Λ  .The accuracy is guaranteed because the coefficients for the images moving into the same piece of region are almost equal to an identical constant in a short period, which is explained in Figure 11.
It has been mentioned that the optical flow is not sensitive to satellite's orbit motion and earth rotation in a short term; namely, the possible values are assigned by the following functions: =   (, ,  0 , Ω, ,    ,    , Δ) ,   =   (, ,  0 , Ω, ,    ,    , Δ) , Here N is the number of constant-valued segments in the region encompassing all the possible optical flow trajectories.The orbital elements and integral step size Δ are common to all functions.Furthermore, when long term measurements are executed, Ξ  and Λ  only need to be renewed according to the current parameters.
The coefficient matrix of the optical flow equations for th (1 ≤  ≤ ) group can be written as where   = max{ 1 , . . .,   }.
Consequently, as we organize the equations for all groups, the global coefficient matrix will be given in the following form: C is a quasidiagonal partitioned matrix; every subblock has 2 rows.The maximal columns of C are  max = max{ 1 , . . .,   }.
The unknown variables are as follows: The constant are as follows: Δu has been measured by image dense registration.s can be determined by auxiliary data of sensors.The global equations are expressed by As for this problem, it is easy to be verified that conditions.
(1) 2 > 3 max , (2) rank(C) = 3 max easily meet well in practical works.To solve (44), well-posedness is the critical issue for the inverse problem.Strong nonlinearity and anisotropy of optical flow will greatly reduce the relevance between the coefficients in C; meanwhile it increases the wellposedness of the solution.The least-square solution of (47) can be obtained: The well-posedness is able to be examined by Singular Value Decomposition (SVD) to C. Consider the nonnegative definite matrix C  C whose eigenvalues are given in order: where U 2×2 and V 3 max ×3 max are unit orthogonal matrices and the singular values are   = √  .The well-posedness of the solution is acceptable if condition number (C) =  1 / 3 max ≤ .
Associating the process of inverse problem solving in Section 4 with the process of preliminary information acquisition in Section 3, the whole algorithm for remote sensor's attitude measurement is illustrated in the flow chart in Figure 12.

Experimental Results and Validation.
In the experiment, 72940 samples on 7 image strip pairs were involved.Considering maintaining the values in Ξ and Λ nearly invariant, we redistributed these samples into 20 subspaces and solved out the three axial components of the angular velocity.According to Shannon's sampling theorem, the measurable frequency   is expected to reach up to the half of line rates of TDI CCD.For the experiment,   ≈ 1.749 KHz.The   ∼  curves of 0 s ∼ 0.148 s are shown in Figure 13.
In this period,  2 max = 0.001104 ∘ /s,   were perturbing the remote sensor; besides, compared to the signals of  1 () and  2 (), the low frequency components in  3 () are higher in magnitude.Actually, according to the remote sensor, satellite yaw angle is needed to be regulated in real time to compensate for the image rotation on the focal plane such that the detectors can always scan along the direction of image motion.Based on the auxiliary data, the image motion velocity vector ⃗ V of the central pixel in FOV can be computed.So the optimal yaw motion in principle will be The mean value of  * 3 (),  * 3 = 0.01198 ∘ /s.We attribute Δ * 3 =  3 −  * 3 = 0.00554 ∘ /s to the error of satellite attitude control.
In order to validate the measurement, the technique of template reconfiguration was implemented again to check the expected phenomenon that, based on the high-accurate information, the correlations between the new templates and   should be further improved.In addition, the distribution of γ near γmax is going to become more compact, which is easy to be understood since much more useful information about remote sensor's motion is introduced into template reconstructions and increases the similarities between the lapped images.
Unlike the processing in image dense registration, in the validation phase, larger original templates are selected.Let  1 be the referenced image template which centered at the examining element,  2 the new template reconfigured by rough prediction of optical flow, T2 the new template reconfigured based on precision attitude motion measurement, and   the template on sensed image which centered at the registration pixel.For all templates,  =  = 101.The distributions of the normalized cross-correlation coefficients corresponding to the referenced template centered on the sampled selected in .1000row belonging to .7 CCD with sensed image belonging to .8 CCD are illustrated in Figure 14.
where  max and  max are, respectively, the column and row number of the peak-valued location.
Judging from the results, the performances in case (c) are better than those in case (b) and much more better than those in case (a), since the precise attitude motion measurements enhance the precision of optical inversion so as to improve the similarities between the new templates and sensed images.Note that, although in case (b) the variance decreases slightly, as we have analyzed in Section 3.2, compared to case (a), the offsets of centroids from the peaks have been corrected well by the use of the rough optical flow predictions.

Summary and Discussions.
In terms of the preceding sections, we can see that, comparing to ordinary NCC, the precision of image registration is greatly improved, since it is attributed to the assistance of the technique of template reconfiguration.Implementing the auxiliary data from the space-borne sensors to optical flow prediction, the relative deformations between the lapped image pairs can be computed in considerable accuracy.Afterwards, it will be used to estimate the gray values of the corresponding parts on sensed images and help us to construct a new template for registration.As we know, the space-borne sensors may give middle and low frequency components of imager's attitude motion in excellent precision.Thus, comparing to the classical direct template based registration algorithms, the similarity between the reconfigured template and sensed images may greatly increase.Furthermore, the minor deformations attributed to high frequency jitters can be detected by using subpixel registration between the reconfigured templates and sensed images.This point of view is the exact basis of high frequency jitters measurement with optical flow inversion.

Conclusion
In this paper, optical flows and time-varying image deformation in space dynamic imaging are analyzed in detail.The nonlinear and anisotropic image motion velocity and optical flows are utilized to strengthen the well-posedness of the inverse problem of attitude precise measurement by optical

( 6 )
F(  −       ): the frame of focal plane.Axes     and     lie in the focal plane.They are, respectively, parallel to     2 and     1 .Axis     coincides with the optical axis.

( 7 ) 3 Figure 1 :
Figure 1: The analysis of dynamic imaging for the three-dimensional planet.

Figure 3 :
Figure 3: Nonlinear image velocity field and optical flow trajectories influence the time-phase relations between the lapped images captured by the adjacent overlapped detectors.

Figure 4 :
Figure 4: The structure of Mapping Satellite-1 and its effective payloads.

Figure 5 :
Figure 5: Optical flow prediction and template reconfiguration.

Figure 9 :Figure 10 :
Figure 9: The offsets of lapped images captured by CCD1 and CCD2.

Figure 11 :
Figure 11: Coefficients Determination according to the Current Location of the Image.

) 2 Figure 12 :
photography frame between T 2 and

Figure 13 :
Figure 13: Solutions for the angular velocities of the remote sensor.
(a) shows the situation for  1 and   (b) for  2 and   , and (c) for T2 and   .The compactness of the data is characterized by the peak value  max and the location variances  2  ,  2  :

Table 1 :
The offsets between overlapped images.ishasbeen discovered in Figures9 and 10that the fluctuation of image offsets taking place in Gap 1 is greater in magnitude than in Gap 3.All the facts indicate that the distorted optical flows can be detected from a plenty of image offsets.We will see later that the nonlinear distribution of the data strengthens the well-posedness of optical flow inversion algorithm.
1 max = 0.001194 ∘ /s The signal of  3 () is fluctuating around mean value  3 = 0.01752 ∘ /s.It is not hard to infer that high frequency jitters