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

Abstract

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.

1. 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–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 space-borne 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–10]. References [7, 8] reported the camera on 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 two-dimensional 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 phase-only 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–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.

2. Image Velocity Field Analysis

Suppose that a large FOV camera is performing push-broom 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.

Figure 1 

The analysis of dynamic imaging for the three-dimensional planet.

In order to set up the model of space imaging, some coordinate systems need to be defined as follows.(1)ℐ(𝑂_𝑒−𝑋𝑌𝑍): the inertial frame of the earth. For convenience, here we choose 𝐽2000 frame. The origin 𝑂_𝑒 is located at the earth center.(2)𝒞(𝑜_𝑠−𝑥^′₁𝑥^′₂𝑥^′₃): the frame of camera. Axis 𝑜_𝑠𝑥^′₃ is the optical axis, and origin 𝑜_𝑠 is the center of the exit pupil.(3)𝒪(𝑜_𝑠−𝑢₁𝑢₂𝑢₃): the orbit frame. Axis 𝑜_𝑠𝑢₃ passes through the center of earth, and axis 𝑜_𝑠𝑢₂ is perpendicular to the instant orbit plane.(4)ℬ(𝑜_𝑠−𝑥_𝑏𝑦_𝑏𝑧_𝑏): the body frame of the satellite.(5)𝒫(𝑜−𝑥𝑦𝑧): 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.(6)ℱ(𝑜^′−𝑥^′𝑦^′𝑧^′): the frame of focal plane. Axes 𝑜^′𝑥^′ and 𝑜^′𝑦^′ lie in the focal plane. They are, respectively, parallel to 𝑜_𝑠𝑥^′₂ and 𝑜_𝑠𝑥^′₁. Axis 𝑜^′𝑧^′ coincides with the optical axis.(7)ℰ(𝑂_𝑒−𝑥_𝑒𝑦_𝑒𝑧_𝑒): the frame of Terrestrial Reference Frame (TRF). Axis 𝑂_𝑒𝑦_𝑒 points to the North Pole, and axis 𝑂_𝑒𝑥_𝑒 passes through the intersection of Greenwich meridian and the equator.

According to Figure 1, 𝜕₀ is the ground track of the satellite and 𝜕₁ and 𝜕₂ are the ground traces corresponding to two fixed boresights in FOV, which are far away from 𝜕₀ if the imager holds a large attitude angle. Obviously, the shapes and lengths of 𝜕₁ and 𝜕₂ 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, ^ℐ𝜌_𝑝 represents the vector measured in frame ℐ, and accordingly ^𝒞𝜌_𝑝 is the same vector measured in frame 𝒞. We select one unit vector 𝜏 which is tangent to the surface of the earth at 𝑝. Let 𝑟(⃗𝑥₁,𝑥₂,𝑥₃) be the position vector of 𝑝 relative to 𝑜_𝑠; then 𝑟⃗^𝒞 and 𝑟̇^𝒞⃗ characterize the apparent motion of 𝑝. Assume that the image point 𝑝^′ is formed on the focal plane with coordinates (𝑥^′₁,𝑥^′₂,𝑥^′₃) in frame 𝒞. 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 𝑥^′_𝑖=𝛽𝑥_𝑖,(),𝑥𝑖=1,2_′3=𝑓^′, (1) where 𝑓^′ is the effective focal length, the lateral magnification of 𝑝^′: 𝛽=(−1)^𝑚−1(𝑓^′/(𝑟⃗⋅𝑒₃)), 𝑚 is the number of intermediate images in the optical system, and 𝑒_𝑖(𝑖=1,2,3) is the base of 𝒞.

Let 𝑟_𝑠 be the position vector of satellite relative to 𝑂_𝑒; then 𝑟⃗=𝜌−𝑟_𝑠. In imaging, the flight trajectory of the satellite platform in ℐ can be treated as Keplerian orbit, as illustrated in Figure 2. According to the orbit elements, 𝑖₀, 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 𝑀_𝑡=𝑀₀+𝑛(𝑡−𝑡₀), where 𝑛=2𝜋/𝑃, 𝑃 is the orbit period [11]: 𝑀_𝑡−(𝐸−𝑒sin𝐸)=0. (2) In frame 𝒪, ^𝒪𝑟_𝑠=𝑎()0,cos𝐸−𝑒𝑏sin𝐸^𝒪𝑣_𝑠=0𝑛−𝑎sin𝐸𝑏cos𝐸1−𝑒cos𝐸. (3) The coordinate transform matrix between 𝒪 and ℐ is 𝐓_𝒪ℐ=𝐶𝜔𝐶Ω−𝑆𝜔𝐶𝑖₀𝑆Ω−𝑆𝜔𝐶Ω−𝐶𝜔𝐶𝑖₀𝑆Ω𝑆𝑖₀𝑆Ω𝐶𝜔𝑆Ω+𝑆𝜔𝐶𝑖₀𝐶Ω−𝑆𝜔𝑆Ω+𝐶𝜔𝐶𝑖₀𝐶Ω−𝑆𝑖₀𝐶Ω𝑆𝜔𝑆𝑖₀𝐶𝜔𝑆𝑖₀𝐶𝑖₀. (4) For simplicity, we write 𝐶𝛼∶=cos𝛼,𝑆𝛼∶=sin𝛼.

Figure 2 

Orbital motion of remote sensor.

In engineering, the coordinate transfer matrix 𝐓_𝒪ℐ also can be derived from the real-time measurements of GPS. Since the base vectors of frame 𝒪 in ℐ: 𝑢₃=−^ℐ𝑟_𝑠/|𝑟_𝑠|, 𝑢₂=(^ℐ𝑣_𝑠×^ℐ𝑟_𝑠)/|𝑣_𝑠×𝑟_𝑠|, and 𝑢₁=𝑢₂×𝑢₃ then 𝐓_𝒪ℐ=𝑢₁𝑢₂𝑢₃⁻¹: ^ℐ𝑟_𝑠=𝐓_𝒪ℐ⋅^𝒪𝑟_𝑠,^ℐ𝑣_𝑠=𝐓_𝒪ℐ⋅^𝒪𝑣_𝑠. (5) Associating the equation of boresight with the ellipsoid surface of the earth in 𝒞 yields 𝑋²+𝑍²𝐴²_𝑒+𝑌²𝐵²_𝑒=1,𝑋−𝑋_𝑠𝑠₁=𝑋−𝑌_𝑠𝑠₂=𝑋−𝑍_𝑠𝑠₃. (6) Here, 𝐴_𝑒=6378.137 km and 𝐵_𝑒=6356.752 km being the length of earth's semimajor axis and semiminor axis. 𝑠_𝑖(𝑖=1,2,3) are the unit vectors of 𝑟⃗^ℐ. We write the solution of (7) as^ℐ𝜌=𝑋𝑌𝑍^𝑇. Hence, 𝑟⃗^ℐ=^ℐ𝜌−^ℐ𝑟_𝑠,𝑟⃗^𝒞=𝐌⋅𝐀⋅𝐓⁻¹_𝒪ℐ⋅^ℐ𝑟⃗, where 𝐌 is the coordinate transformation matrix from frame ℬ to frame 𝒞; it is a constant matrix for fixed installation; 𝐀 is the attitude matrix of satellite; according to 1-2-3 rotating order, we have 𝐀=𝐑_𝜓⋅𝐑_𝜃⋅𝐑_𝜑 (7) in which 𝐑_𝜓=cos𝜓_𝑡sin𝜓_𝑡0−sin𝜓_𝑡cos𝜓_𝑡0,𝐑001_𝜃=cos𝜃_𝑡0−sin𝜃_𝑡010sin𝜃_𝑡0cos𝜃_𝑡,𝐑_𝜑=1000cos𝜑_𝑡sin𝜑_𝑡0−sin𝜑_𝑡cos𝜑_𝑡, (8) where 𝜑_𝑡, 𝜃_𝑡, and 𝜓_𝑡 are in order, the real-time roll angle, pitch angle, and yaw angle at moment 𝑡. The velocity of 𝑝 in 𝒞 can be written in the following scalar form: 𝑥̇_𝑖=^𝒞𝑟̇⃗⋅𝑒_𝑖(𝑖=1,2,3). (9) Thus, the velocity of image point of 𝑝^′ will be 𝑥̇^′_𝑖=𝛽̇𝑥_𝑖+𝛽𝑥̇_𝑖=(−1)^𝑚𝑓^′𝑟⃗⋅𝑒₃²𝑥̇₃𝑥_𝑖+(−1)^𝑚−1𝑓^′𝑟⃗⋅𝑒₃𝑥̇_𝑖(𝑖=1,2). (10) Substituting (2)–(9) into (10), the velocity vector of image point 𝑣^′=(𝑥̇^′₁,𝑥̇^′₂)^𝑇 can be expressed as the explicit function of several variables; that is, 𝑣^′=𝒱𝑖₀,Ω,𝜔,𝑒,𝑀_𝑡0,𝜑_𝑡,𝜃_𝑡,𝜓_𝑡,𝜑̇_𝑡,𝜃̇_𝑡,𝜓̇_𝑡,𝑥^′₁,𝑥^′₂. (11) For conciseness, this analytical expression of 𝑣^′ 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: 𝜔₁𝜔₂𝜔₃=𝐑_𝜓00𝜓̇_𝑡+𝐑_𝜃0𝜃̇_𝑡0+𝐑_𝜑00𝜑̇_𝑡, (12)𝜔₁, 𝜔₂, and 𝜔₃ are the three components of the remote sensor's angular velocity ^𝒞𝜔_𝑠 relative to orbital frame 𝒪, which is calibrated in frame 𝒞. 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.

3. 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 𝑥^′₁(𝑇)=^𝑇₀𝑥̇^′₁𝑥^′₁,𝑥^′₂𝑥,𝑡𝑑𝑡,_′2(𝑇)=^𝑇₀𝑥̇^′₂𝑥^′₁,𝑥^′₂,𝑡𝑑𝑡. (13) Since (13) are coupled nonlinear integral equations, we convert them to numerical forms and solve them iteratively, 𝑥^′_𝑖(0)=𝑥^′_𝑖(𝑡)_𝑡=0,𝑥^′_𝑗(𝑛)=𝑥^′_𝑗()+1𝑛−12𝑥̇^′_𝑗𝑥^′₁(𝑛−1),𝑥^′₂()+𝑥̇𝑛−1,𝑛_′𝑗𝑥^′₁(𝑛−1),𝑥^′₂(𝑛−1),𝑛−1Δ𝑡,𝑗=1,2;𝑛∈ℤ⁺. (14)

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.

3.1. 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 𝛿𝑥^′_𝑖=𝛿𝛽𝑥_𝑖+𝛽𝛿𝑥_𝑖. (15) We expand the term of 𝛿𝛽: 𝛿𝛽=(−1)^𝑚𝑓^′(𝑟⃗+𝛿𝑟)⃗⋅𝑒₃−𝑓^′𝑟⃗⋅𝑒₃=(−1)^𝑚−1𝑓^′𝑟⃗⋅𝑒₃^∞_𝑘=1(−1)^𝑘𝛿𝑟⃗⋅𝑒₃𝑟⃗⋅𝑒₃^𝑘≈(−1)^𝑚𝑓^′𝜏⋅𝑒₃𝛿𝑙𝑟⃗⋅𝑒₃². (16) Taking derivatives with respect to variable 𝑡 for either part of (15), we have 𝛿𝑥̇^′_𝑖=𝛿𝛽̇𝑥_𝑖+𝛿𝛽𝑥̇_𝑖+𝛽̇𝛿𝑥_𝑖+𝛽𝛿𝑥̇_𝑖. (17) 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 ^𝒞(𝛿𝑟)⃗=𝛿𝑙𝐌𝐀𝐓⁻¹_𝒪ℐ𝐓^ℰ_ℰℐ𝜏. (18) Notice that ^ℰ𝜏 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, ^𝒞𝛿𝑟̇⃗𝛿𝑙_𝜏=𝐌̇𝐀𝐓⁻¹_𝒪ℐ𝐓_ℰℐ+𝐌𝐀̇𝐓⁻¹_𝒪ℐ𝐓_ℰℐ+𝐌𝐀𝐓̇⁻¹_𝒪ℐ𝐓_ℰℐ+𝐌𝐀𝐓⁻¹_𝒪ℐ𝐓̇_ℰℐ^ℰ𝜏, (19) where the coordinate transform matrix from frame ℰ to ℐ is 𝐓_ℰℐ=cos𝐻_𝑝0−sin𝐻_𝑝010sin𝐻_𝑝0cos𝐻_𝑝. (20) 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 𝑡₁∼𝑡₂ can be written as 𝛿𝑥^′_𝑖^𝑡2_𝜏−𝛿𝑥^′_𝑖^𝑡1_𝜏=^𝑡2_𝑡1𝛿𝑥̇^′_𝑖𝜏𝑑𝑡. (21) From (17), we have 𝛿𝑥̇^′_𝑖𝜏=𝛿𝑙𝛿𝛽𝑥̇𝛿𝑙^𝑖+𝛽̇𝛿𝑥_𝑖𝛿𝑥̇𝛿𝑙+𝛽^𝑖.𝛿𝑙 (22) According to (16), (18), and (19) we obtain the terms in (22): 𝛿𝛽=()𝛿𝑙−1^𝑚𝑓^′𝒞𝜏⋅𝑒₃𝑟⃗⋅𝑒₃²,𝛿𝑥_𝑖𝛿𝑙=𝐌𝐀𝐓⁻¹_𝒪ℐ𝐓_ℰℐ𝑒_𝑖⋅^ℰ𝛿𝑥̇𝜏,_𝑖=𝛿𝑙^𝒞𝛿𝑟̇⃗𝛿𝑙_𝜏⋅𝑒_𝑖+^𝒞𝛿𝑟⃗𝛿𝑙_𝜏⋅𝑒_𝑖̇⃗. (23) Furthermore, if the camera is fixed to the satellite platform, then 𝐌̇=0,𝑒_𝑖̇⃗=0.

Consequently, (22) becomes 𝓕_𝑖(𝜏)=𝛿𝑥̇𝑡,^′𝑖𝜏=()𝛿𝑙−1^𝑚𝑓^′𝒞𝜏⋅𝑒₃𝑟⃗⋅𝑒₃²𝑥̇_𝑖+(−1)^𝑚𝑓^′𝑟̇⃗⋅𝑒_𝑖𝑟⃗⋅𝑒₃²𝐌𝐀𝐓⁻¹_𝒪ℐ𝐓_ℰℐ𝑒_𝑖⋅^ℰ𝜏+(−1)^𝑚−1𝑓^′𝑟⃗⋅𝑒₃𝐌𝐀̇𝐓⁻¹_𝒪ℐ𝐓_ℰℐ+𝐌𝐀𝐓̇⁻¹_𝒪ℐ𝐓_ℰℐ+𝐌𝐀𝐓⁻¹_𝒪ℐ𝐓̇_ℰℐ^ℰ𝜏⋅𝑒_𝑖. (24) For the motionless scene on the earth surface, ^ℰ𝜏 is a time-independent but space-dependent unit tangent vector, which meanwhile represents a specific orientation on the ground. Moreover, the physical meaning of function 𝓕_𝑖(𝑡,𝜏) is the image deformation of unit-length curve on the curved surface along the direction of ^ℰ𝜏 in unit time interval. That is the instantaneous space-time deforming rate of the image of the object along ^ℰ𝜏.

Consequently, in dynamic imaging, macroscopic deformation on the moving image can be derived from the integral of 𝓕_𝑖(𝑡,𝜏) 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 ^ℰ𝜏=𝒯(𝑠) be a vector-valued function with variable 𝑠 (the length of the arc) which is time-invariant in frame ℰ and gives the tangent vectors along the curve.

So, the image deformation taking place during 𝑡₁∼𝑡₂ is able to be described as 𝑥^′_𝑝𝑖^𝑡2_𝑡1−𝑥^′_𝑞𝑖^𝑡2_𝑡1=_Γ^𝑡2_𝑡1𝓕_𝑖∘𝒯𝑑𝑡𝑑𝑠 (25) 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.

3.2. 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 time-varying 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 {𝛼₁,𝛼₂,…} (in the left gap), {𝛽₁,𝛽₂,…} (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 𝒞 between the lapped images can be determined by measuring the offsets of the sample image element pairs in frame 𝒫.

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.

Let Δ𝑦^′=Δ𝑥^′₁,Δ𝑥^′=Δ𝑥^′₂ be the relative offsets of the same object’s image on the two photos; they are all calibrated in 𝒞 or ℱ. 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–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–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 space-borne 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.

Figure 4 

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

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.

Figure 5 

Optical flow prediction and template reconfiguration.

Step 1. Set the original lapped image strips (the images which were acquired directly by the detectors and without any postprocessing) in frame 𝒞.

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 𝑇₁, which consists of 𝑁² square elements; that is, dim(𝑇₁)=𝑁×𝑁. 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 𝜏^′₁∼𝜏^′₄ in frame 𝒞: 𝜏^′₁=√22𝑒₁−√22𝑒₂,𝜏^′₃√=−22𝑒₁+√22𝑒₂,𝜏^′₂=√22𝑒₁+√22𝑒₂,𝜏^′₄√=−22𝑒₁−√22𝑒₂. (26) Suppose image point 𝑝^′ is the center of an arbitrary element Σ^′ in 𝑇₁. Let Σ be the area element on the earth surface which is conjugate to Σ^′. The four unit radial vectors of the vertexes on Σ, 𝜏₁∼𝜏₄ are conjugate to 𝜏^′₁∼𝜏^′₄ and tangent to the earth surface at 𝑝. From the geometrical relations, we have ^𝒞𝜏_𝑖=(−1)^𝑚𝑟^′×𝜏^′_𝑖×^𝒞𝑛_𝑝𝑟^′×𝜏^′_𝑖×^𝒞𝑛_𝑝ℰ,𝜏_𝑖=𝐓⁻¹_ℰℐ𝐓_𝒪ℐ𝐀⁻¹𝐌^−1𝒞𝜏_𝑖,_𝒞𝑛_𝑝=𝐌𝐀𝐓⁻¹_𝒪ℐ𝐓^ℰ_ℰℐ𝑛_𝑝, (27) where ^ℰ𝑛_𝑝 is the unit normal vector of Σ at 𝑝. We predict the deformations along 𝜏₁∼𝜏₄ during 𝑡₁∼𝑡₂ according to the measurements of GPS, star trackers, and gyroscopes, as explained in Figure 6. 𝑡₁ is the imaging time on prior detector and 𝑡₂ is the imaging time on the posterior detector, 𝛿𝑥^′_1𝜏^Δ𝑡_𝑘=𝛿𝑥^′₁^𝑡2_𝜏𝑘−𝛿𝑥^′₁^𝑡1_𝜏𝑘,𝛿𝑥^′_2𝜏^Δ𝑡_𝑘=𝛿𝑥^′₂^𝑡2_𝜏𝑘−𝛿𝑥^′₂^𝑡1_𝜏𝑘,(𝑘=1∼4). (28) The shape of deformed image Σ^′_𝑡2 can be got through linear interpolation with 𝛿𝑟^′_𝜏^Δ𝑡_𝑘=𝛿𝑥^′_1𝜏^Δ𝑡_𝑘𝛿𝑥^′_2𝜏^Δ𝑡_𝑘. (29)

Figure 6 

Deformation of single element.

Step 3. Reconfigure referenced template 𝑇₁ according to optical flow prediction, and then get a new template 𝑇₂.
Let 𝑇^′₁ be the deformed image of 𝑇₁ computed in Step 2. Let 𝜒=𝐵_𝑖,𝑗 be the central element of 𝑇^′₁; integers 𝑖 and 𝑗 are, respectively, the row number and column number of 𝐵_𝑖,𝑗. The gray value 𝑙_𝑖𝑗 of each element in 𝑇^′₁ is equal to its counterpart in 𝑇₁ with the same indexes. In addition, we initialize a null template 𝑇₀ whose shape and orientation are identical to 𝑇₁; the central element of 𝑇₀ is denoted by 𝑇_𝑖,𝑗.
Then, we cover 𝑇₀ upon 𝑇^′₁ and let their centers coincide; that is, 𝑇_𝑖,𝑗=𝐵_𝑖,𝑗, as shown in Figure 7. Denote the vertexes of 𝑇^′₁ as 𝑉^𝑘_𝑖𝑗(𝑘=1∼4). Therefore, the connective relation for adjacent elements can be expressed by 𝑉¹_𝑖,𝑗=𝑉²_{𝑖,𝑗−1}=𝑉³_{𝑖−1,𝑗−1}=𝑉⁴_{𝑖−1,𝑗}.
Next, we will reassign the gray value ℎ^′_𝑖𝑗 to 𝑇_𝑖,𝑗(𝑖=1⋯𝑁,𝑗=1⋯𝑁) in sequence to construct a new template 𝑇₂. The process is just a simulation of image resample when optical flow arrives at the posterior detector, as indicated in Figure 3.
That is, ℎ^′_𝑖𝑗=^𝑖+1_{𝑚=𝑖−1}^𝑗+1_{𝑛=𝑗−1}𝜂_𝑚𝑛𝑙_𝑚𝑛. (30) Weight coefficient 𝜂_𝑚𝑛=𝑆_𝑚𝑛/𝑤² where 𝑆_𝑚𝑛 is the area of the intersecting polygon of 𝐵_𝑚,𝑛 with 𝑇_𝑖,𝑗.

Figure 7 

Template reconfiguration.

Step 4. Computenormalized cross-correlation coefficients between 𝑇₂ and the sensed image, and then determine the subpixel offset of 𝑇₂ relative to the sensed image in frame 𝒫.
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 (𝑛₁,𝑛₂) on 𝑇_𝑠, the normalized cross-correlation (NCC) coefficient is given by𝛾𝑛₁,𝑛₂=∑_𝑥,𝑦𝑔𝑥,𝑦−𝑔_𝑥,𝑦ℎ𝑥−𝑛₁,𝑦−𝑛₂−ℎ∑_𝑥,𝑦𝑔𝑥,𝑦−𝑔_𝑥,𝑦²∑_𝑥,𝑦ℎ𝑥−𝑛₁,𝑦−𝑛₂−ℎ^20.5, (31) where 𝑔_𝑥,𝑦 is the mean gray value of the segment of 𝑇_𝑠 that is masked by 𝑇₂ and ℎ is the mean of 𝑇₂. Equation (31) requires approximately 𝑁²(𝑀−𝑁+1)² additions and 𝑁²(𝑀−𝑁+1)² multiplications, whereas the complexity of FFT algorithm needs about 12𝑀²log₂𝑀 real multiplications and 18𝑀²log₂𝑀 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 (𝑘,𝑚),𝑘,𝑚∈ℤ 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 𝒫.
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 𝒢(𝑢,𝑣) be the 2D Discrete Fourier Transforms (DFT) of the sensed window; then we have 𝒢(𝑢,𝑣)=^(𝑁−1)/2_{𝑥=−(𝑁−1)/2}^(𝑁−1)/2_{𝑦=−(𝑁−1)/2}𝑔𝑥,𝑦𝑊_𝑀^𝑢𝑥𝑊_𝑀^𝑣𝑦,ℋ(𝑢,𝑣)=^(𝑁−1)/2_{𝑥=−(𝑁−1)/2}^(𝑁−1)/2_{𝑦=−(𝑁−1)/2}ℎ𝑥,𝑦𝑊_𝑁^𝑢𝑥𝑊_𝑁^𝑣𝑦. (32) Here 𝑊_𝑁=exp−𝑗2𝜋𝑁. (33)
Cross-phase spectrum is given by ℛ()=𝒢()ℋ𝑢,𝑣𝑢,𝑣^∗(𝑢,𝑣)𝒢(𝑢,𝑣)ℋ^∗((),𝑢,𝑣)=exp𝑗𝜙𝑢,𝑣 (34) where ℋ^∗ is the complex conjugate of ℋ. By inverse Discrete Fourier Transform (IDFT) we have 𝛾𝑛₁,𝑛₂=1𝑁^{2(𝑁−1)/2}_{𝑢=−(𝑁−1)/2}^(𝑁−1)/2_{𝑣=−(𝑁−1)/2}ℛ(𝑢,𝑣)𝑊^−𝑢𝑛1_𝑁𝑊^−𝑣𝑛2_𝑁. (35) Suppose that the new peak 𝛾_max appears at (𝑘^′,𝑚^′),𝑘^′,𝑚^′∈ℤ, referring to [27]; we have the following relation: 𝛾_max𝑘^′,𝑚^′≈𝜆𝑁²sin𝜋𝑘^′+𝛿_𝑥sin𝜋𝑚^′+𝛿_𝑦sin(𝜋/𝑁)𝑘^′+𝛿_𝑥sin(𝜋/𝑁)𝑚^′+𝛿_𝑦. (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 𝒫 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 𝒫, and Δ𝑥^′_𝑛,Δ𝑦^′_𝑛 (mm) in frame ℱ∶𝛿𝑛_𝑥=col_𝑟+col_𝑠−𝑄−𝜂_𝑛,Δ𝑥^′_𝑛=Δ𝑥^′_2𝑛=𝛿𝑛_𝑥⋅𝑤,𝛿𝑛_𝑦=row_𝑠−row_𝑟−𝐷𝑤,Δ𝑦^′_𝑛=Δ𝑥^′_1𝑛=𝛿𝑛_𝑦⋅𝑤+𝐷. (37) 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 cross-track 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 hand, it is has been discovered in Figures 9 and 10 that 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.

Table 1 

The offsets between overlapped images.

Figure 8 

Dense image registration for lapped image strips, CCD1 versus CCD2 (Gap 1: the left two) and CCD3 versus CCD4 (Gap. 3 the right two).

Figure 9 

The offsets of lapped images captured by CCD1 and CCD2.

Figure 10 

The offsets of lapped images captured by CCD3 and CCD4.

4. Remote Sensor Attitude Motion Measurement

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

4.1. The Principle of Optical Inversion

For clarity, in frame 𝒞, 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 𝑁_𝑖𝑗=𝑡_𝑖,𝑗−𝜏_𝑖,𝑗Δ𝑡∈ℤ⁺,𝑖=1⋯𝑛,𝑗=1⋯𝑚, (38) 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 𝜔₁, 𝜔₂, and 𝜔₃ (the variables in the inverse problem), yielding the linear optical flow equations.

For the 𝑙th group samples, ^𝑁1𝑙_𝑖=𝑙𝑐^𝑖_1𝑙1𝜔^𝑖₁+𝑐^𝑖_1𝑙2𝜔^𝑖₂+𝑐^𝑖_1𝑙3𝜔^𝑖₃=Δ𝑥^′_1𝑙−𝑠_𝑥1,_{𝑁1𝑙𝑖=𝑙}𝑑^𝑖_1𝑙1𝜔^𝑖₁+𝑑^𝑖_1𝑙2𝜔^𝑖₂+𝑑^𝑖_1𝑙3𝜔^𝑖₃=Δ𝑦^′_1𝑙−𝑠_𝑦1,_{𝑁𝑛𝑙𝑖=𝑙}⋮𝑐^𝑖_𝑛𝑙1𝜔^𝑖₁+𝑐^𝑖_𝑛𝑙2𝜔^𝑖₂+𝑐^𝑖_𝑛𝑙3𝜔^𝑖₃=Δ𝑥^′_𝑛𝑙−𝑠_𝑥𝑛,_{𝑁𝑛𝑙𝑖=𝑙}𝑑^𝑖_𝑛𝑙1𝜔^𝑖₁+𝑑^𝑖_𝑛𝑙2𝜔^𝑖₂+𝑑^𝑖_𝑛𝑙3𝜔^𝑖₃=Δ𝑦^′_𝑛𝑙−𝑠_𝑦𝑛. (39)

Suppose that the sample process will stop until 𝑚 groups have been founded. The coefficients are as follows: 𝑐^𝑖_𝜇𝜈𝜅=Ξ_𝜅𝜇,𝑖−𝜈+1𝑁_𝜇𝜈𝒩,𝑑^𝑖_𝜇𝜈𝜅=Λ_𝜅𝜇,𝑖−𝜈+1𝑁_𝜇𝜈𝒩(𝜅=1,2,3). (40) Here Ξ_𝑘=𝜉_11𝑘𝜉_12𝑘⋯𝜉_1𝒩𝑘𝜉_21𝑘𝜉_22𝑘⋯𝜉_2𝒩𝑘𝜉⋯⋯_𝑛1𝑘𝜉_𝑛2𝑘⋯𝜉_𝑛𝒩𝑘,Λ_𝑘=𝜆_11𝑘𝜆_12𝑘⋯𝜆_1𝒩𝑘𝜆_21𝑘𝜆_22𝑘⋯𝜆_2𝒩𝑘𝜆⋯⋯_𝑛1𝑘𝜆_𝑛2𝑘⋯𝜆_𝑛𝒩𝑘. (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.

Figure 11 

Coefficients Determination according to the Current Location of the Image.

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: 𝜉_𝑖𝑗𝑘=𝜉_𝑘𝑎,𝑒,𝑖₀,Ω,𝜔,𝑥^′_𝑞,𝑦^′_𝑞,𝜆,Δ𝑡_𝑖𝑗𝑘=𝜆_𝑘𝑎,𝑒,𝑖₀,Ω,𝜔,𝑥^′_𝑞,𝑦^′_𝑞,Δ𝑡,𝑖=1∼𝑛,𝑗=1∼𝒩,𝑞=1∼𝒩. (42)

Here 𝒩 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 𝐂_𝑗=𝑐¹_1𝑗1𝑐¹_1𝑗2𝑐¹_1𝑗3⋯𝑐^𝑁1𝑗_1𝑗1𝑐^𝑁1𝑗_1𝑗2𝑐^𝑁1𝑗_1𝑗3𝑑⋯00_11𝑗1𝑑¹_1𝑗2𝑑¹_1𝑗3⋯𝑑^𝑁1𝑗_1𝑗1𝑑^𝑁1𝑗_1𝑗2𝑑^𝑁1𝑗_1𝑗3𝑐⋯00⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯_1𝑞𝑗1𝑐¹_𝑞𝑗2𝑐¹_𝑞𝑗3⋯⋯⋯⋯𝑐^𝑁𝑞𝑗_𝑞𝑗1𝑐^𝑁𝑞𝑗_𝑞𝑗1𝑐^𝑁𝑞𝑗_𝑞𝑗3𝑑¹_𝑞𝑗1𝑑¹_𝑞𝑗2𝑑¹_𝑞𝑗3⋯⋯⋯⋯𝑑^𝑁𝑞𝑗_𝑞𝑗1𝑑^𝑁𝑞𝑗_𝑞𝑗2𝑑^𝑁𝑞𝑗_𝑞𝑗3𝑐⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯_1𝑛𝑗1𝑐¹_𝑛𝑗2𝑐¹_𝑛𝑗3⋯⋯𝑐^𝑁𝑛𝑗_𝑛𝑗1𝑐^𝑁𝑛𝑗_𝑛𝑗2𝑐^𝑁𝑛𝑗_𝑛𝑗3𝑑⋯0_1𝑛𝑗1𝑑¹_𝑛𝑗2𝑑¹_𝑛𝑗3⋯⋯𝑑¹_𝑛𝑗1𝑑¹_𝑛𝑗2𝑑¹_𝑛𝑗3⋯0_{2𝑛×3𝑁𝑞𝑗}, (43) where 𝑁_𝑞𝑗=max{𝑁_1𝑗,…,𝑁_𝑛𝑗}.

Consequently, as we organize the equations for all groups, the global coefficient matrix will be given in the following form: 𝐂𝓒=^{12𝑛×3𝑁𝑞1}0𝐂0⋯⋯⋯⋯0_{22𝑛×3𝑁𝑞2}[𝐂]0⋯⋯⋯0⋱⋯⋯_{2𝑛×3𝑁max}⋱𝐂_{𝑚2𝑛×3𝑁𝑞𝑚}0_{2𝑛𝑚×3𝑁max}. (44)

𝓒 is a quasidiagonal partitioned matrix; every subblock has 2𝑛 rows. The maximal columns of 𝓒 are 𝑁_max=max{𝑁_𝑞1,…,𝑁_𝑞𝑚}.

The unknown variables are as follows: [Ω]_3𝑁max×1=𝜔¹₁𝜔¹₂𝜔¹₃⋯𝜔^𝑁max₁𝜔^𝑁max₂𝜔^𝑁max₃^𝑇. (45)

The constant are as follows: Δ𝐮_2𝑚𝑛×1=Δ𝑥^′₁₁Δ𝑦^′₁₁⋯Δ𝑥^′_𝑛1Δ𝑦^′_𝑛1⋯Δ𝑥^′_𝑚1Δ𝑦^′_𝑚1⋯Δ𝑥^′_𝑛𝑚Δ𝑦^′_𝑛𝑚^𝑇,𝐬_2𝑚𝑛×1=𝑠_𝑥1𝑠_𝑦1⋯𝑠_𝑥𝑛𝑠_𝑦𝑛⋯𝑠_𝑥1𝑠_𝑦1⋯𝑠_𝑥𝑛𝑠_𝑦𝑛^𝑇. (46)

Δ𝐮 has been measured by image dense registration. 𝐬 can be determined by auxiliary data of sensors. The global equations are expressed by 𝓒_{2𝑚𝑛×3𝑁max}⋅[Ω]_3𝑁max×1=Δ𝐮_2𝑚𝑛×1−𝐬_2𝑚𝑛×1. (47) As for this problem, it is easy to be verified that conditions.

(1)  2𝑛𝑚>3𝑁_max,  (2)  rank(𝓒)=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 𝓒; meanwhile it increases the well-posedness of the solution. The least-square solution of (47) can be obtained: [Ω]=𝓒^𝑇𝓒⁻¹𝓒^𝑇(Δ𝐮−𝐬). (48) The well-posedness is able to be examined by Singular Value Decomposition (SVD) to 𝓒. Consider the nonnegative definite matrix 𝓒^𝑇𝓒 whose eigenvalues are given in order: 𝜆₁≥𝜆₂≥⋯≥𝜆_3𝑁max, 𝓒=𝐔[𝜎]𝐕^𝑇, (49) where 𝐔_{2𝑚𝑛×2𝑚𝑛} and 𝐕_{3𝑁max×3𝑁max} are unit orthogonal matrices and the singular values are 𝜎_𝑖=𝜆_𝑖. The well-posedness of the solution is acceptable if condition number 𝜅(𝓒)=𝜎₁/𝜎_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.

Figure 12 

The flow chart of the attitude motion measurement.

4.2. 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 0s~0.148s are shown in Figure 13.

Figure 13 

Solutions for the angular velocities of the remote sensor.

In this period, 𝜔_2max=0.001104^∘/s, 𝜔_1max=0.001194^∘/s The signal of 𝜔₃(𝑡) is fluctuating around mean value 𝜔₃=0.01752^∘/s. It is not hard to infer that high frequency jitters were perturbing the remote sensor; besides, compared to the signals of 𝜔₁(𝑡) and 𝜔₂(𝑡), the low frequency components in 𝜔₃(𝑡) 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 𝑣 of the central pixel in FOV can be computed. So the optimal yaw motion in principle will be 𝜓^∗_𝑡=𝑣_𝑦′𝑣_𝑥′,𝜔^∗₃(𝑡)=𝜓̇^∗_𝑡=𝑣̇_𝑦′𝑣_𝑥′−𝑣_𝑦′𝑣̇_𝑥′𝑣²_𝑥′. (50) The mean value of 𝜔^∗₃(𝑡),𝜔^∗₃=0.01198^∘/s. We attribute Δ𝜔^∗₃=𝜔₃−𝜔^∗₃=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 𝑇₁ be the referenced image template which centered at the examining element, 𝑇₂ the new template reconfigured by rough prediction of optical flow, 𝑇₂ 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 𝑁𝑜.1000 row belonging to 𝑁𝑜.7 CCD with sensed image belonging to 𝑁𝑜.8 CCD are illustrated in Figure 14.

(a)

(b)

(c)

Figure 14 

Normalized cross-correlations comparison. ((a) shows the distribution of 𝛾 by applying direct NCC algorithm. (b) shows the distribution of 𝛾 after template reconfiguration with optical flow prediction. (c) shows the distribution of 𝛾 derived from posterior template reconfiguration with high-accurate senor’s attitude measurement. It can be noticed that the values of 𝛾 tend to be distributed uniformly around the peak value location from left to right.)

(a) shows the situation for 𝑇₁ and 𝑇_𝑠 (b) for 𝑇₂ and 𝑇_𝑠, and (c) for 𝑇₂ and 𝑇_𝑠. The compactness of the data is characterized by the peak value 𝛾_max and the location variances 𝜎²_𝑥,𝜎²_𝑦∶𝜎²_𝑥=∑^𝑀_𝑖=1∑^𝑀_𝑗=1𝛾_𝑖𝑗⋅𝑖−𝑥_max²∑^𝑀_𝑖=1∑^𝑀_𝑗=1𝛾_𝑖𝑗,𝜎²_𝑦=∑^𝑀_𝑖=1∑^𝑀_𝑗=1𝛾_𝑖𝑗⋅𝑗−𝑦_max²∑^𝑀_𝑖=1∑^𝑀_𝑗=1𝛾_𝑖𝑗, (51) where 𝑥_max and 𝑦_max are, respectively, the column and row number of the peak-valued location.

In case (a), 𝛾_max(𝑎)=0.893, standard deviation 𝜎_𝑥(𝑎)=56.53, and 𝜎_𝑦(𝑎)=81.92; in case (b), 𝛾_max(𝑏)=0.918,𝜎_𝑥(𝑏)=48.39, and 𝜎_𝑦(𝑏)=66.86; in case (c),  𝛾_max(𝑐)=0.976; however the variance sharply shrinks to 𝜎_𝑥(𝑐)=3.27,𝜎_𝑦(𝑐)=4.06. In Table 2, some other samples with 1000 rows interval are also examined. The samples can be regarded as independent to each other.

Table 2 

Correlation coefficients distribution for registration templates.

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.

4.3. 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.

5. 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 flow inversion method. For the purpose of determining the conditions of fixed solutions of optical flow equations, information based image registration algorithms are proposed. We apply rough optical flow prediction to improve the efficiency and accuracy of dense image registration. Based on the results of registration, the attitude motions of remote sensors in imaging are measured by using precise optical flow inversion method. The experiment on a remote sensor showed that the measurements are achieved in very high accuracy as well as with broad bandwidth. This method can extensively be used in remote sensing missions such as image strips splicing, geometrical rectification, and nonblind image restoration to promote the surveying precision and resolving power.

Conflict of Interests

The authors declare that they have no financial nor personal relationships with other people or organizations that can inappropriately influence their work; there is no professional or other personal interest of any nature or kind in any product, service, and/or company that could be construed as influencing the position presented in, or the review of, this paper.

Acknowledgments

This work is supported by the National High Technology Research and Development Program of China (863 Program) (Grant no. 2012AA121503, Grant no. 2013AA12260, and Grant no. 2012AA120603) and the National Natural Science Foundation of China (Grant no. 61377012).