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

Authors: Chong Wang; Zheng You; Fei Xing; Borui Zhao; Bin Li; Gaofei Zhang; Qingchang Tao
Journal: Mathematical Problems in Engineering
 (2014)

---

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


---

## Body

## 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 Figure1. 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)
ℐ(Oe-XYZ): the inertial frame of the earth. For convenience, here we choose J2000 frame. The origin Oe is located at the earth center.(2)
𝒞(os-x1′x2′x3′): the frame of camera. Axis osx3′ is the optical axis, and origin os is the center of the exit pupil.(3)
𝒪(os-u1u2u3): the orbit frame. Axis osu3 passes through the center of earth, and axis osu2 is perpendicular to the instant orbit plane.(4)
ℬ(os-xbybzb): the body frame of the satellite.(5)
𝒫(o-xyz): the frame of photography. The origin o is the center of the photo. Axis ox points to the column direction, and axis oy points to the row direction.(6)
ℱ(o′-x′y′z′): the frame of focal plane. Axes o′x′ and o′y′ lie in the focal plane. They are, respectively, parallel to osx2′ and osx1′. Axis o′z′ coincides with the optical axis.(7)
ℰ(Oe-xeyeze): the frame of Terrestrial Reference Frame (TRF). Axis Oeye points to the North Pole, and axis Oexe passes through the intersection of Greenwich meridian and the equator.According to Figure1, ∂0 is the ground track of the satellite and ∂1 and ∂2 are the ground traces corresponding to two fixed boresights in FOV, which are far away from ∂0 if the 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 pointp on the earth, its position vector relative to Oe is denoted as ρ→p. As a convention in the following discussions, ρ→pℐ represents the vector measured in frame ℐ, and accordingly ρ→p𝒞 is the same vector measured in frame 𝒞. We select one unit vector τ→ which is tangent to the surface of the earth at p. Let r→(x1,x2,x3) be the position vector of p relative to os; then r→𝒞 and r→˙𝒞 characterize the apparent motion of p. Assume that the image point p′ is formed on the focal plane with coordinates (x1′,x2′,x3′) 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)xi′=βxi,(i=1,2),x3′=f′,
where f′ is the effective focal length, the lateral magnification of p′: β=(-1)m-1(f′/(r→·e→3)), m is the number of intermediate images in the optical system, and e→i(i=1,2,3) is the base of 𝒞.Letr→s be the position vector of satellite relative to Oe; then r→=ρ→-r→s. 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, i0, inclination; Ω, longitude of ascending node; ω, argument of perigee; a, semimajor axis; e, eccentricity; Mt, mean anomaly at epoch, we implement Newton-Raphson method to solve (2) and get the eccentric anomaly E from the given mean anomaly Mt=M0+n(t-t0), where n=2π/P, P is the orbit period [11]:
(2)Mt-(E-esinE)=0.
In frame 𝒪,
(3)𝒪r→s=(a(cos⁡E-e)bsinE0),𝒪v→s=(-asinEbcos⁡E0)n1-ecos⁡E.
The coordinate transform matrix between 𝒪 and ℐ is
(4)T𝒪ℐ=(CωCΩ-SωCi0SΩ-SωCΩ-CωCi0SΩSi0SΩCωSΩ+SωCi0CΩ-SωSΩ+CωCi0CΩ-Si0CΩSωSi0CωSi0Ci0).
For simplicity, we write Cα:=cos⁡α,Sα≔sinα.Figure 2
Orbital motion of remote sensor.In engineering, the coordinate transfer matrixT𝒪ℐ also can be derived from the real-time measurements of GPS. Since the base vectors of frame 𝒪 in ℐ: u→3=-r→sℐ/|r→s|, u→2=(v→sℐ×r→sℐ)/|v→s×r→s|, and u→1=u→2×u→3 then T𝒪ℐ=(u→1u→2u→3)-1:
(5)ℐr→s=T𝒪ℐ·𝒪r→s,ℐv→s=T𝒪ℐ·𝒪v→s.
Associating the equation of boresight with the ellipsoid surface of the earth in 𝒞 yields
(6)X2+Z2Ae2+Y2Be2=1,X-Xss1=X-Yss2=X-Zss3.
Here, Ae=6378.137 km and Be=6356.752 km being the length of earth's semimajor axis and semiminor axis. si(i=1,2,3) are the unit vectors of ℐr→. We write the solution of (7) asℐρ→=(XYZ)T. Hence, r→ℐ=ρ→ℐ-r→sℐ,𝒞r→=M·A·T𝒪ℐ-1·r→ℐ, where M is the coordinate transformation matrix from frame ℬ to frame 𝒞; it is a constant matrix for fixed installation; A is the attitude matrix of satellite; according to 1-2-3 rotating order, we have
(7)A=Rψ·Rθ·Rφ
in which
(8)Rψ=(cos⁡ψtsinψt0-sinψtcos⁡ψt0001),Rθ=(cos⁡θt0-sinθt010sinθt0cos⁡θt),Rφ=(1000cos⁡φtsinφt0-sinφtcos⁡φt),
where φt, θt, and ψt are in order, the real-time roll angle, pitch angle, and yaw angle at moment t. The velocity of p in 𝒞 can be written in the following scalar form:
(9)x˙i=r→˙𝒞·e→i(i=1,2,3).
Thus, the velocity of image point of p′ will be
(10)x˙i′=β˙xi+βx˙i=(-1)mf′(r→·e→3)2x˙3xi+(-1)m-1f′r→·e→3x˙i(i=1,2).
Substituting (2)–(9) into (10), the velocity vector of image point v→′=(x˙1′,x˙2′)T can be expressed as the explicit function of several variables; that is,
(11)v→′=𝒱(i0,Ω,ω,e,Mt0,φt,θt,ψt,φ˙t,θ˙t,ψ˙t,x1′,x2′).
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 φt,θt, and ψt can be roughly measured by the star trackers and GPS. Meanwhile, their time rates φ˙t, θ˙t, and ψ˙t have the following relations:
(12)(ω1ω2ω3)=Rψ{(00ψ˙t)+Rθ[(0θ˙t0)+Rφ(00φ˙t)]},ω1, ω2, and ω3 are the three components of the remote sensor's angular velocity 𝒞ω→s 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(13)x1′(T)=∫0Tx˙1′(x1′,x2′,t)dt,x2′(T)=∫0Tx˙2′(x1′,x2′,t)dt.
Since (13) are coupled nonlinear integral equations, we convert them to numerical forms and solve them iteratively,
(14)xi′(0)=xi′(t)|t=0,xj′(n)=xj′(n-1)+12{x˙j′[x1′(n-1),x2′(n-1),n]+x˙j′[x1′(n-1),x2′(n-1),hhhhhhhn-1x′1(n-1)x2′(n-1),]}Δt,(j=1,2;n∈ℤ+).It is evident that the algorithm has enough precision so long as the step-size of time intervalΔt 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 Figure1, considering a microspatial variation of p along τ→ on the curved surface can be expressed as δρ→p=δlτ→. Its conjugated image is
(15)δxi′=δβxi+βδxi.
We expand the term of δβ:
(16)δβ=(-1)m{f′(r→+δr→)·e→3-f′r→·e→3}=(-1)m-1f′r→·e→3∑k=1∞(-1)k(δr→·e→3r→·e→3)k≈(-1)mf′τ→·e→3δl(r→·e→3)2.
Taking derivatives with respect to variable t for either part of (15), we have
(17)δx˙i′=δβ˙xi+δβx˙i+β˙δxi+βδx˙i.
According to (16) we know that δβ˙≈0. On the other hand, the variation of r→ can be expressed through a series of coordinate transformations; that is
(18)(δr→)𝒞=δl[MAT𝒪ℐ-1Tℰℐτ→ℰ].
Notice that τ→ℰ is a fixed tangent vector of earth surface at object point p, which is time-invariant and specifies an orientation of motionless scene on the earth.Consequently,(19)(δr→˙𝒞δl)τ→=(M˙AT𝒪ℐ-1Tℰℐ+MA˙T𝒪ℐ-1Tℰℐ+MAT˙𝒪ℐ-1Tℰℐ+MAT𝒪ℐ-1T˙ℰℐ)ℰτ→,
where the coordinate transform matrix from frame ℰ to ℐ is
(20)Tℰℐ=(cos⁡Hp0-sinHp010sinHp0cos⁡Hp).
Let ωe be the angular rate of the earth and αp the longitude of p on the earth; then the hour angle of p at time t is Hp(t)=GST+αp+ωet, in which, GST represents Greenwich sidereal time.The microscale image deformation of the extended scene on the earth along the direction ofτ→ during t1~t2 can be written as
(21)[δxi′]τ→t2-[δxi′]τ→t1=∫t1t2(δx˙i′)τ→dt.
From (17), we have
(22)(δx˙i′)τ→δl=δβδlx˙i+β˙δxiδl+βδx˙iδl.
According to (16), (18), and (19) we obtain the terms in (22):
(23)δβδl=(-1)mf′τ→𝒞·e→3(r→·e→3)2,δxiδl=MAT𝒪ℐ-1Tℰℐe→i·τ→ℰ,δx˙iδl=(δr→˙𝒞δl)τ→·e→i+(δr→𝒞δl)τ→·e→˙i.
Furthermore, if the camera is fixed to the satellite platform, then M˙=0,e→˙i=0.Consequently, (22) becomes
(24)Fi(t,τ→)=(δx˙i′)τ→δl=(-1)mf′τ→𝒞·e→3(r→·e→3)2x˙i+(-1)mf′(r→˙·e→i)(r→·e→3)2MAT𝒪ℐ-1Tℰℐe→i·τ→ℰ+(-1)m-1f′r→·e→3(MA˙T𝒪ℐ-1Tℰℐ+MAT˙𝒪ℐ-1Tℰℐ+MAT𝒪ℐ-1T˙ℰℐ)τ→ℰ·e→i.
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 Fi(t,τ→) 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 ofFi(t,τ→) 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 p,q∈Γ, and let their Gaussian images p′,q′∈Γ′. Let ℰτ→=𝒯→(s) be a vector-valued function with variable s (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 duringt1~t2 is able to be described as
(25)[(xp′)i]t1t2-[(xq′)i]t1t2=∫Γ∫t1t2Fi∘𝒯→dtds
in which Fi∘𝒯→=Fi[t,𝒯→(s)].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 {α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 𝒞 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Δy′=Δx1′,Δx′=Δx2′ 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 Figure3, 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 asT1, which consists of N2 square elements; that is, dim⁡(T1)=N×N. Let χ be its central element and w the width of each element; here w=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 Figure5. The unit radial vectors are denoted by τ→1′~τ→4′ in frame 𝒞:
(26)τ→1′=22e→1-22e→2,τ→3′=-22e→1+22e→2,τ→2′=22e→1+22e→2,τ→4′=-22e→1-22e→2.
Suppose image point p′ is the center of an arbitrary element Σ′ in T1. 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 p. From the geometrical relations, we have
(27)𝒞τ→i=(-1)mr→′×τ→i′×n→p𝒞|r→′×τ→i′×n→p𝒞|,τ→iℰ=Tℰℐ-1T𝒪ℐA-1M-1τ→i𝒞,𝒞n→p=MAT𝒪ℐ-1Tℰℐℰn→p,
where ℰn→p is the unit normal vector of Σ at p. We predict the deformations along τ→1~τ→4 during t1~t2 according to the measurements of GPS, star trackers, and gyroscopes, as explained in Figure 6. t1 is the imaging time on prior detector and t2 is the imaging time on the posterior detector,
(28)[δx1′]τ→kΔt=[δx1′]τ→kt2-[δx1′]τ→kt1,[δx2′]τ→kΔt=[δx2′]τ→kt2-[δx2′]τ→kt1,(k=1~4).
The shape of deformed image Σt2′ can be got through linear interpolation with
(29)[δr→′]τ→kΔt=([δx1′]τ→kΔt[δx2′]τ→kΔt).Figure 6
Deformation of single element.Step 3.
Reconfigure referenced templateT1 according to optical flow prediction, and then get a new template T2.
LetT1′ be the deformed image of T1 computed in Step 2. Let χ=Bi,j be the central element of T1′; integers i and j are, respectively, the row number and column number of Bi,j. The gray value lij of each element in T1′ is equal to its counterpart in T1 with the same indexes. In addition, we initialize a null template T0 whose shape and orientation are identical to T1; the central element of T0 is denoted by Ti,j.
Then, we coverT0 upon T1′ and let their centers coincide; that is, Ti,j=Bi,j, as shown in Figure 7. Denote the vertexes of T1′ as Vijk(k=1~4). Therefore, the connective relation for adjacent elements can be expressed by Vi,j1=Vi,j-12=Vi-1,j-13=Vi-1,j4.
Next, we will reassign the gray valuehij′ to Ti,j(i=1⋯N,j=1⋯N) in sequence to construct a new template T2. The process is just a simulation of image resample when optical flow arrives at the posterior detector, as indicated in Figure 3.
That is,(30)hij′=∑m=i-1i+1∑n=j-1j+1ηmnlmn.
Weight coefficient ηmn=Smn/w2 where Smn is the area of the intersecting polygon of Bm,n with Ti,j.Figure 7
Template reconfiguration.Step 4.
Computenormalized cross-correlation coefficients betweenT2 and the sensed image, and then determine the subpixel offset of T2 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 Step2 Assuming that the search space is Ts, dim⁡(Ts)=M×M. When Ti,j moves to the pixel (n1,n2) on Ts, the normalized cross-correlation (NCC) coefficient is given by(31)γ^(n1,n2)=∑x,y[g(x,y)-g-x,y][h(x-n1,y-n2)-h-]{∑x,y[g(x,y)-g-x,y]2∑x,y[h(x-n1,y-n2)-h-]2}0.5,
where g-x,y is the mean gray value of the segment of Ts that is masked by T2 and h- is the mean of T2. Equation (31) requires approximately N2(M-N+1)2 additions and N2(M-N+1)2 multiplications, whereas the complexity of FFT algorithm needs about 12M2log2M real multiplications and 18M2log2M real additions/subtractions [32, 33].
At the beginning, we takeM=101,N=7 and compute the NCC coefficient. When M is much larger than N, the calculation in spatial domain will be efficient. Suppose that the peak value γ^max⁡ is taken at the coordinate (k,m),k,m∈ℤ in the sensed window. Hence we will reduce search space into a smaller one with dimension of 47×47 which centered on Ts(k,m). Next, the subpixel registration is realized by phase correlation algorithm with larger M and N to suppress the system errors owing to the deficiencies of detailed textures on the photo. Here we take M=47,N=23. Let the subpixel offset between the two registering image elements be denoted as δx and δy in frame 𝒫.
The phase correlation algorithm in the frequency domain becomes more efficient asN approaches M 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𝒢(u,v) be the 2D Discrete Fourier Transforms (DFT) of the sensed window; then we have
(32)𝒢(u,v)=∑x=-(N-1)/2(N-1)/2∑y=-(N-1)/2(N-1)/2g(x,y)WMuxWMvy,ℋ(u,v)=∑x=-(N-1)/2(N-1)/2∑y=-(N-1)/2(N-1)/2h(x,y)WNuxWNvy.
Here
(33)WN=exp⁡(-j2πN).
Cross-phase spectrum is given by(34)ℛ^(u,v)=𝒢(u,v)ℋ*(u,v)|𝒢(u,v)ℋ*(u,v)|=exp⁡(jϕ(u,v)),
where ℋ* is the complex conjugate of ℋ. By inverse Discrete Fourier Transform (IDFT) we have
(35)γ^(n1,n2)=1N2∑u=-(N-1)/2(N-1)/2∑v=-(N-1)/2(N-1)/2ℛ^(u,v)WN-un1WN-vn2.
Suppose that the new peak γ^max⁡ appears at (k′,m′),k′,m′∈ℤ, referring to [27]; we have the following relation:
(36)γ^max⁡(k′,m′)≈λN2sin[π(k′+δx)]sin[π(m′+δy)]sin[(π/N)(k′+δx)]sin[(π/N)(m′+δy)].
The right side presents the spatial distribution of the normalized cross-correlation coefficients. Therefore, (δx,δy) 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 Step1~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.
Letcolr,rowr be the column and row indexes of image elements on the referenced image, and let cols, rows be the indexes of the same elements on the sensed image. The total columns of each detector Q=4096 pix, and the vertical distance between the two detector arrays D=18.4975 mm. According to the results of registration, we get the offsets of images at nth gap, δnx (cross track), δny (along track) in frame 𝒫, and Δxn′,Δyn′ (mm) in frame ℱ:(37)δnx=colr+cols-Q-ηn,Δxn′=Δ(x2′)n=δnx·w,δny=rows-rowr-Dw,Δyn′=Δ(x1′)n=δny·w+D.
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 Figures9 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.
Sample
Row no.(pixel)
δ
n
x (pixel)
Δ
x
n
′ (mm)
δ
n
y
(pixel)
Δ
y
n
′ (mm)
S11
258
−25.15
−0.2200625
−5.39
18.4503
S12
423
−23.78
−0.2080750
−7.36
18.4331
S31
266
−12.85
−0.1124375
−7.66
18.4304
S32
436
−12.97
−0.1134875
−6.87
18.4374Figure 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.


## 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 Figure1, considering a microspatial variation of p along τ→ on the curved surface can be expressed as δρ→p=δlτ→. Its conjugated image is
(15)δxi′=δβxi+βδxi.
We expand the term of δβ:
(16)δβ=(-1)m{f′(r→+δr→)·e→3-f′r→·e→3}=(-1)m-1f′r→·e→3∑k=1∞(-1)k(δr→·e→3r→·e→3)k≈(-1)mf′τ→·e→3δl(r→·e→3)2.
Taking derivatives with respect to variable t for either part of (15), we have
(17)δx˙i′=δβ˙xi+δβx˙i+β˙δxi+βδx˙i.
According to (16) we know that δβ˙≈0. On the other hand, the variation of r→ can be expressed through a series of coordinate transformations; that is
(18)(δr→)𝒞=δl[MAT𝒪ℐ-1Tℰℐτ→ℰ].
Notice that τ→ℰ is a fixed tangent vector of earth surface at object point p, which is time-invariant and specifies an orientation of motionless scene on the earth.Consequently,(19)(δr→˙𝒞δl)τ→=(M˙AT𝒪ℐ-1Tℰℐ+MA˙T𝒪ℐ-1Tℰℐ+MAT˙𝒪ℐ-1Tℰℐ+MAT𝒪ℐ-1T˙ℰℐ)ℰτ→,
where the coordinate transform matrix from frame ℰ to ℐ is
(20)Tℰℐ=(cos⁡Hp0-sinHp010sinHp0cos⁡Hp).
Let ωe be the angular rate of the earth and αp the longitude of p on the earth; then the hour angle of p at time t is Hp(t)=GST+αp+ωet, in which, GST represents Greenwich sidereal time.The microscale image deformation of the extended scene on the earth along the direction ofτ→ during t1~t2 can be written as
(21)[δxi′]τ→t2-[δxi′]τ→t1=∫t1t2(δx˙i′)τ→dt.
From (17), we have
(22)(δx˙i′)τ→δl=δβδlx˙i+β˙δxiδl+βδx˙iδl.
According to (16), (18), and (19) we obtain the terms in (22):
(23)δβδl=(-1)mf′τ→𝒞·e→3(r→·e→3)2,δxiδl=MAT𝒪ℐ-1Tℰℐe→i·τ→ℰ,δx˙iδl=(δr→˙𝒞δl)τ→·e→i+(δr→𝒞δl)τ→·e→˙i.
Furthermore, if the camera is fixed to the satellite platform, then M˙=0,e→˙i=0.Consequently, (22) becomes
(24)Fi(t,τ→)=(δx˙i′)τ→δl=(-1)mf′τ→𝒞·e→3(r→·e→3)2x˙i+(-1)mf′(r→˙·e→i)(r→·e→3)2MAT𝒪ℐ-1Tℰℐe→i·τ→ℰ+(-1)m-1f′r→·e→3(MA˙T𝒪ℐ-1Tℰℐ+MAT˙𝒪ℐ-1Tℰℐ+MAT𝒪ℐ-1T˙ℰℐ)τ→ℰ·e→i.
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 Fi(t,τ→) 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 ofFi(t,τ→) 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 p,q∈Γ, and let their Gaussian images p′,q′∈Γ′. Let ℰτ→=𝒯→(s) be a vector-valued function with variable s (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 duringt1~t2 is able to be described as
(25)[(xp′)i]t1t2-[(xq′)i]t1t2=∫Γ∫t1t2Fi∘𝒯→dtds
in which Fi∘𝒯→=Fi[t,𝒯→(s)].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 {α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 𝒞 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Δy′=Δx1′,Δx′=Δx2′ 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 Figure3, 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 asT1, which consists of N2 square elements; that is, dim⁡(T1)=N×N. Let χ be its central element and w the width of each element; here w=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 Figure5. The unit radial vectors are denoted by τ→1′~τ→4′ in frame 𝒞:
(26)τ→1′=22e→1-22e→2,τ→3′=-22e→1+22e→2,τ→2′=22e→1+22e→2,τ→4′=-22e→1-22e→2.
Suppose image point p′ is the center of an arbitrary element Σ′ in T1. 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 p. From the geometrical relations, we have
(27)𝒞τ→i=(-1)mr→′×τ→i′×n→p𝒞|r→′×τ→i′×n→p𝒞|,τ→iℰ=Tℰℐ-1T𝒪ℐA-1M-1τ→i𝒞,𝒞n→p=MAT𝒪ℐ-1Tℰℐℰn→p,
where ℰn→p is the unit normal vector of Σ at p. We predict the deformations along τ→1~τ→4 during t1~t2 according to the measurements of GPS, star trackers, and gyroscopes, as explained in Figure 6. t1 is the imaging time on prior detector and t2 is the imaging time on the posterior detector,
(28)[δx1′]τ→kΔt=[δx1′]τ→kt2-[δx1′]τ→kt1,[δx2′]τ→kΔt=[δx2′]τ→kt2-[δx2′]τ→kt1,(k=1~4).
The shape of deformed image Σt2′ can be got through linear interpolation with
(29)[δr→′]τ→kΔt=([δx1′]τ→kΔt[δx2′]τ→kΔt).Figure 6
Deformation of single element.Step 3.
Reconfigure referenced templateT1 according to optical flow prediction, and then get a new template T2.
LetT1′ be the deformed image of T1 computed in Step 2. Let χ=Bi,j be the central element of T1′; integers i and j are, respectively, the row number and column number of Bi,j. The gray value lij of each element in T1′ is equal to its counterpart in T1 with the same indexes. In addition, we initialize a null template T0 whose shape and orientation are identical to T1; the central element of T0 is denoted by Ti,j.
Then, we coverT0 upon T1′ and let their centers coincide; that is, Ti,j=Bi,j, as shown in Figure 7. Denote the vertexes of T1′ as Vijk(k=1~4). Therefore, the connective relation for adjacent elements can be expressed by Vi,j1=Vi,j-12=Vi-1,j-13=Vi-1,j4.
Next, we will reassign the gray valuehij′ to Ti,j(i=1⋯N,j=1⋯N) in sequence to construct a new template T2. The process is just a simulation of image resample when optical flow arrives at the posterior detector, as indicated in Figure 3.
That is,(30)hij′=∑m=i-1i+1∑n=j-1j+1ηmnlmn.
Weight coefficient ηmn=Smn/w2 where Smn is the area of the intersecting polygon of Bm,n with Ti,j.Figure 7
Template reconfiguration.Step 4.
Computenormalized cross-correlation coefficients betweenT2 and the sensed image, and then determine the subpixel offset of T2 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 Step2 Assuming that the search space is Ts, dim⁡(Ts)=M×M. When Ti,j moves to the pixel (n1,n2) on Ts, the normalized cross-correlation (NCC) coefficient is given by(31)γ^(n1,n2)=∑x,y[g(x,y)-g-x,y][h(x-n1,y-n2)-h-]{∑x,y[g(x,y)-g-x,y]2∑x,y[h(x-n1,y-n2)-h-]2}0.5,
where g-x,y is the mean gray value of the segment of Ts that is masked by T2 and h- is the mean of T2. Equation (31) requires approximately N2(M-N+1)2 additions and N2(M-N+1)2 multiplications, whereas the complexity of FFT algorithm needs about 12M2log2M real multiplications and 18M2log2M real additions/subtractions [32, 33].
At the beginning, we takeM=101,N=7 and compute the NCC coefficient. When M is much larger than N, the calculation in spatial domain will be efficient. Suppose that the peak value γ^max⁡ is taken at the coordinate (k,m),k,m∈ℤ in the sensed window. Hence we will reduce search space into a smaller one with dimension of 47×47 which centered on Ts(k,m). Next, the subpixel registration is realized by phase correlation algorithm with larger M and N to suppress the system errors owing to the deficiencies of detailed textures on the photo. Here we take M=47,N=23. Let the subpixel offset between the two registering image elements be denoted as δx and δy in frame 𝒫.
The phase correlation algorithm in the frequency domain becomes more efficient asN approaches M 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𝒢(u,v) be the 2D Discrete Fourier Transforms (DFT) of the sensed window; then we have
(32)𝒢(u,v)=∑x=-(N-1)/2(N-1)/2∑y=-(N-1)/2(N-1)/2g(x,y)WMuxWMvy,ℋ(u,v)=∑x=-(N-1)/2(N-1)/2∑y=-(N-1)/2(N-1)/2h(x,y)WNuxWNvy.
Here
(33)WN=exp⁡(-j2πN).
Cross-phase spectrum is given by(34)ℛ^(u,v)=𝒢(u,v)ℋ*(u,v)|𝒢(u,v)ℋ*(u,v)|=exp⁡(jϕ(u,v)),
where ℋ* is the complex conjugate of ℋ. By inverse Discrete Fourier Transform (IDFT) we have
(35)γ^(n1,n2)=1N2∑u=-(N-1)/2(N-1)/2∑v=-(N-1)/2(N-1)/2ℛ^(u,v)WN-un1WN-vn2.
Suppose that the new peak γ^max⁡ appears at (k′,m′),k′,m′∈ℤ, referring to [27]; we have the following relation:
(36)γ^max⁡(k′,m′)≈λN2sin[π(k′+δx)]sin[π(m′+δy)]sin[(π/N)(k′+δx)]sin[(π/N)(m′+δy)].
The right side presents the spatial distribution of the normalized cross-correlation coefficients. Therefore, (δx,δy) 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 Step1~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.
Letcolr,rowr be the column and row indexes of image elements on the referenced image, and let cols, rows be the indexes of the same elements on the sensed image. The total columns of each detector Q=4096 pix, and the vertical distance between the two detector arrays D=18.4975 mm. According to the results of registration, we get the offsets of images at nth gap, δnx (cross track), δny (along track) in frame 𝒫, and Δxn′,Δyn′ (mm) in frame ℱ:(37)δnx=colr+cols-Q-ηn,Δxn′=Δ(x2′)n=δnx·w,δny=rows-rowr-Dw,Δyn′=Δ(x1′)n=δny·w+D.
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 Figures9 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.
Sample
Row no.(pixel)
δ
n
x (pixel)
Δ
x
n
′ (mm)
δ
n
y
(pixel)
Δ
y
n
′ (mm)
S11
258
−25.15
−0.2200625
−5.39
18.4503
S12
423
−23.78
−0.2080750
−7.36
18.4331
S31
266
−12.85
−0.1124375
−7.66
18.4304
S32
436
−12.97
−0.1134875
−6.87
18.4374Figure 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 kth sample element belonging to nth lapped strip pair are written as Δxnk′, Δynk′. 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 sx,sy.Letτi,j,ti,j be in order the two sequential imaging time of the jth image sample on the overlapped detectors in jth 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
(38)Nij=[ti,j-τi,jΔt]∈ℤ+,(i=1⋯n,j=1⋯m),
where n is the amount of CCD gaps, m is the amount of sample groups, and Δt is the time step. We set samples with same j 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 thelth group samples,
(39)∑i=lN1lc1l1iω1i+c1l2iω2i+c1l3iω3i=Δx1l′-sx1,∑i=lN1ld1l1iω1i+d1l2iω2i+d1l3iω3i=Δy1l′-sy1,⋮∑i=lNnlcnl1iω1i+cnl2iω2i+cnl3iω3i=Δxnl′-sxn,∑i=lNnldnl1iω1i+dnl2iω2i+dnl3iω3i=Δynl′-syn.Suppose that the sample process will stop untilm groups have been founded. The coefficients are as follows:
(40)cμνκi=Ξκ(μ,⌈i-ν+1Nμν𝒩⌉),dμνκi=Λκ(μ,⌈i-ν+1Nμν𝒩⌉)(κ=1,2,3).
Here
(41)Ξk=(ξ11kξ12k⋯ξ1𝒩kξ21kξ22k⋯ξ2𝒩k⋯⋯ξn1kξn2k⋯ξn𝒩k),Λk=(λ11kλ12k⋯λ1𝒩kλ21kλ22k⋯λ2𝒩k⋯⋯λn1kλn2k⋯λn𝒩k).As for the algorithm, to reduce the complexity, all possible values for the coefficients are stored in the matrixesΞk and Λk. 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:(42)ξijk=ξk(a,e,i0,Ω,ω,xq′,yq′,Δt),iiλijk=λk(a,e,i0,Ω,ω,xq′,yq′,Δt),i=1~n,j=1~𝒩,q=1~𝒩.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 Δt are common to all functions. Furthermore, when long term measurements are executed, Ξk and Λk only need to be renewed according to the current parameters.The coefficient matrix of the optical flow equations forjth (1≤j≤m) group can be written as (43)Cj=(c1j11c1j21c1j31⋯c1j1N1jc1j2N1jc1j3N1j⋯00d1j11d1j21d1j31⋯d1j1N1jd1j2N1jd1j3N1j⋯00⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯cqj11cqj21cqj31⋯⋯⋯⋯cqj1Nqjcqj1Nqjcqj3Nqjdqj11dqj21dqj31⋯⋯⋯⋯dqj1Nqjdqj2Nqjdqj3Nqj⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯cnj11cnj21cnj31⋯⋯cnj1Nnjcnj2Nnjcnj3Nnj⋯0dnj11dnj21dnj31⋯⋯dnj11dnj21dnj31⋯0)2n×3Nqj, where Nqj=max⁡{N1j,…,Nnj}.Consequently, as we organize the equations for all groups, the global coefficient matrix will be given in the following form:(44)C=([C1]2n×3Nq10⋯⋯⋯⋯00[C2]2n×3Nq20⋯⋯⋯0⋱⋯⋯[C]2n×3Nmax⁡⋱[Cm]2n×3Nqm0)2nm×3Nmax⁡.C is a quasidiagonal partitioned matrix; every subblock has 2n rows. The maximal columns of C are Nmax⁡=max⁡{Nq1,…,Nqm}.The unknown variables are as follows:(45)[Ω]3Nmax⁡×1=[ω11ω21ω31⋯ω1Nmax⁡ω2Nmax⁡ω3Nmax⁡]T.The constant are as follows:(46)Δu2mn×1=[Δx11′Δy11′⋯Δxn1′Δyn1′⋯Δxm1′Δym1′⋯Δxnm′Δynm′]T,s2mn×1=[sx1sy1⋯sxnsyn⋯sx1sy1⋯sxnsyn]T.Δ
u has been measured by image dense registration. s can be determined by auxiliary data of sensors. The global equations are expressed by
(47)C2mn×3Nmax⁡·[Ω]3Nmax⁡×1=Δu2mn×1-s2mn×1.
As for this problem, it is easy to be verified that conditions.(
1
)  2nm>3Nmax⁡,  (2)  rank⁡(C)=3Nmax⁡ 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 well-posedness of the solution. The least-square solution of (47) can be obtained:
(48)[Ω]=(CTC)-1CT(Δu-s).
The well-posedness is able to be examined by Singular Value Decomposition (SVD) to C. Consider the nonnegative definite matrix CTC whose eigenvalues are given in order: λ1≥λ2≥⋯≥λ3Nmax⁡,
(49)C=U[σ]VT,
where U2mn×2mn and V3Nmax⁡×3Nmax⁡ are unit orthogonal matrices and the singular values are σi=λi. The well-posedness of the solution is acceptable if condition number κ(C)=σ1/σ3Nmax⁡≤tol.Associating the process of inverse problem solving in Section4 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 fc is expected to reach up to the half of line rates of TDI CCD. For the experiment, fc≈1.749 KHz. The ωi~t 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 ω3(t) is fluctuating around mean value ω-3=0.01752°/s. It is not hard to infer that high frequency jitters were perturbing the remote sensor; besides, compared to the signals of ω1(t) and ω2(t), the low frequency components in ω3(t) 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
(50)ψt*=vy′vx′,ω3*(t)=ψ˙t*=v˙y′vx′-vy′v˙x′vx′2.
The mean value of ω3*(t),ω-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 andTs 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. LetT1 be the referenced image template which centered at the examining element, T2 the new template reconfigured by rough prediction of optical flow, T~2 the new template reconfigured based on precision attitude motion measurement, and Ts the template on sensed image which centered at the registration pixel. For all templates, M=N=101. The distributions of the normalized cross-correlation coefficients corresponding to the referenced template centered on the sampled selected in No.1000 row belonging to No.7 CCD with sensed image belonging to No.8 CCD are illustrated in 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)
(b)
(c)(a) shows the situation forT1 and Ts (b) for T2 and Ts, and (c) for T~2 and Ts. The compactness of the data is characterized by the peak value γmax⁡ and the location variances σx2,σy2:(51)σx2=∑i=1M∑j=1Mγ^ij·(i-xmax⁡)2∑i=1M∑j=1Mγ^ij,σy2=∑i=1M∑j=1Mγ^ij·(j-ymax⁡)2∑i=1M∑j=1Mγ^ij,
where xmax⁡ and ymax⁡ are, respectively, the column and row number of the peak-valued location.In case (a),γ^max⁡(a)=0.893, standard deviation σx(a)=56.53, and σy(a)=81.92; in case (b), γ^max⁡(b)=0.918,σx(b)=48.39, and σy(b)=66.86; in case (c),  γ^max⁡(c)=0.976; however the variance sharply shrinks to σx(c)=3.27,σy(c)=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.
Row number
γ
max
⁡
(
a
,
b
,
c
)
σ
x
~
(
a
,
b
,
c
)
σ
y
~
(
a
,
b
,
c
)
No. 1000
0.893, 0.918, 0.976
56.53, 48.39, 3.27
81.92, 66.86, 4.06
No. 2000
0.807, 0.885, 0.929
87.04, 64.52, 2.13
63.80, 73.42, 5.71
No. 3000
0.832, 0.940, 0.988
49.91, 30.23, 1.55
77.04, 40.16, 1.93
No. 4000
0.919, 0.935, 0.983
50.79, 39.95, 3.61
58.73, 51.55, 3.85
No. 5000
0.865, 0.922, 0.951
59.18, 48.01, 2.37
61.51, 23.71, 2.57
No. 6000
0.751, 0.801, 0.907
125.7, 99.85, 7.89
146.6, 82.13, 2.06
No. 7000
0.759, 0.846, 0.924
116.3, 108.4, 7.14
127.1, 82.67, 4.90
No. 8000
0.884, 0.900, 0.943
81.25, 35.46, 5.42
82.47, 67.70, 2.88Judging 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 Section3.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.


## 4.1. The Principle of Optical Inversion

For clarity, in frame𝒞, the two coordinate components of image displacement of kth sample element belonging to nth lapped strip pair are written as Δxnk′, Δynk′. 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 sx,sy.Letτi,j,ti,j be in order the two sequential imaging time of the jth image sample on the overlapped detectors in jth 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
(38)Nij=[ti,j-τi,jΔt]∈ℤ+,(i=1⋯n,j=1⋯m),
where n is the amount of CCD gaps, m is the amount of sample groups, and Δt is the time step. We set samples with same j 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 thelth group samples,
(39)∑i=lN1lc1l1iω1i+c1l2iω2i+c1l3iω3i=Δx1l′-sx1,∑i=lN1ld1l1iω1i+d1l2iω2i+d1l3iω3i=Δy1l′-sy1,⋮∑i=lNnlcnl1iω1i+cnl2iω2i+cnl3iω3i=Δxnl′-sxn,∑i=lNnldnl1iω1i+dnl2iω2i+dnl3iω3i=Δynl′-syn.Suppose that the sample process will stop untilm groups have been founded. The coefficients are as follows:
(40)cμνκi=Ξκ(μ,⌈i-ν+1Nμν𝒩⌉),dμνκi=Λκ(μ,⌈i-ν+1Nμν𝒩⌉)(κ=1,2,3).
Here
(41)Ξk=(ξ11kξ12k⋯ξ1𝒩kξ21kξ22k⋯ξ2𝒩k⋯⋯ξn1kξn2k⋯ξn𝒩k),Λk=(λ11kλ12k⋯λ1𝒩kλ21kλ22k⋯λ2𝒩k⋯⋯λn1kλn2k⋯λn𝒩k).As for the algorithm, to reduce the complexity, all possible values for the coefficients are stored in the matrixesΞk and Λk. 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:(42)ξijk=ξk(a,e,i0,Ω,ω,xq′,yq′,Δt),iiλijk=λk(a,e,i0,Ω,ω,xq′,yq′,Δt),i=1~n,j=1~𝒩,q=1~𝒩.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 Δt are common to all functions. Furthermore, when long term measurements are executed, Ξk and Λk only need to be renewed according to the current parameters.The coefficient matrix of the optical flow equations forjth (1≤j≤m) group can be written as (43)Cj=(c1j11c1j21c1j31⋯c1j1N1jc1j2N1jc1j3N1j⋯00d1j11d1j21d1j31⋯d1j1N1jd1j2N1jd1j3N1j⋯00⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯cqj11cqj21cqj31⋯⋯⋯⋯cqj1Nqjcqj1Nqjcqj3Nqjdqj11dqj21dqj31⋯⋯⋯⋯dqj1Nqjdqj2Nqjdqj3Nqj⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯cnj11cnj21cnj31⋯⋯cnj1Nnjcnj2Nnjcnj3Nnj⋯0dnj11dnj21dnj31⋯⋯dnj11dnj21dnj31⋯0)2n×3Nqj, where Nqj=max⁡{N1j,…,Nnj}.Consequently, as we organize the equations for all groups, the global coefficient matrix will be given in the following form:(44)C=([C1]2n×3Nq10⋯⋯⋯⋯00[C2]2n×3Nq20⋯⋯⋯0⋱⋯⋯[C]2n×3Nmax⁡⋱[Cm]2n×3Nqm0)2nm×3Nmax⁡.C is a quasidiagonal partitioned matrix; every subblock has 2n rows. The maximal columns of C are Nmax⁡=max⁡{Nq1,…,Nqm}.The unknown variables are as follows:(45)[Ω]3Nmax⁡×1=[ω11ω21ω31⋯ω1Nmax⁡ω2Nmax⁡ω3Nmax⁡]T.The constant are as follows:(46)Δu2mn×1=[Δx11′Δy11′⋯Δxn1′Δyn1′⋯Δxm1′Δym1′⋯Δxnm′Δynm′]T,s2mn×1=[sx1sy1⋯sxnsyn⋯sx1sy1⋯sxnsyn]T.Δ
u has been measured by image dense registration. s can be determined by auxiliary data of sensors. The global equations are expressed by
(47)C2mn×3Nmax⁡·[Ω]3Nmax⁡×1=Δu2mn×1-s2mn×1.
As for this problem, it is easy to be verified that conditions.(
1
)  2nm>3Nmax⁡,  (2)  rank⁡(C)=3Nmax⁡ 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 well-posedness of the solution. The least-square solution of (47) can be obtained:
(48)[Ω]=(CTC)-1CT(Δu-s).
The well-posedness is able to be examined by Singular Value Decomposition (SVD) to C. Consider the nonnegative definite matrix CTC whose eigenvalues are given in order: λ1≥λ2≥⋯≥λ3Nmax⁡,
(49)C=U[σ]VT,
where U2mn×2mn and V3Nmax⁡×3Nmax⁡ are unit orthogonal matrices and the singular values are σi=λi. The well-posedness of the solution is acceptable if condition number κ(C)=σ1/σ3Nmax⁡≤tol.Associating the process of inverse problem solving in Section4 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 fc is expected to reach up to the half of line rates of TDI CCD. For the experiment, fc≈1.749 KHz. The ωi~t 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 ω3(t) is fluctuating around mean value ω-3=0.01752°/s. It is not hard to infer that high frequency jitters were perturbing the remote sensor; besides, compared to the signals of ω1(t) and ω2(t), the low frequency components in ω3(t) 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
(50)ψt*=vy′vx′,ω3*(t)=ψ˙t*=v˙y′vx′-vy′v˙x′vx′2.
The mean value of ω3*(t),ω-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 andTs 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. LetT1 be the referenced image template which centered at the examining element, T2 the new template reconfigured by rough prediction of optical flow, T~2 the new template reconfigured based on precision attitude motion measurement, and Ts the template on sensed image which centered at the registration pixel. For all templates, M=N=101. The distributions of the normalized cross-correlation coefficients corresponding to the referenced template centered on the sampled selected in No.1000 row belonging to No.7 CCD with sensed image belonging to No.8 CCD are illustrated in 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)
(b)
(c)(a) shows the situation forT1 and Ts (b) for T2 and Ts, and (c) for T~2 and Ts. The compactness of the data is characterized by the peak value γmax⁡ and the location variances σx2,σy2:(51)σx2=∑i=1M∑j=1Mγ^ij·(i-xmax⁡)2∑i=1M∑j=1Mγ^ij,σy2=∑i=1M∑j=1Mγ^ij·(j-ymax⁡)2∑i=1M∑j=1Mγ^ij,
where xmax⁡ and ymax⁡ are, respectively, the column and row number of the peak-valued location.In case (a),γ^max⁡(a)=0.893, standard deviation σx(a)=56.53, and σy(a)=81.92; in case (b), γ^max⁡(b)=0.918,σx(b)=48.39, and σy(b)=66.86; in case (c),  γ^max⁡(c)=0.976; however the variance sharply shrinks to σx(c)=3.27,σy(c)=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.
Row number
γ
max
⁡
(
a
,
b
,
c
)
σ
x
~
(
a
,
b
,
c
)
σ
y
~
(
a
,
b
,
c
)
No. 1000
0.893, 0.918, 0.976
56.53, 48.39, 3.27
81.92, 66.86, 4.06
No. 2000
0.807, 0.885, 0.929
87.04, 64.52, 2.13
63.80, 73.42, 5.71
No. 3000
0.832, 0.940, 0.988
49.91, 30.23, 1.55
77.04, 40.16, 1.93
No. 4000
0.919, 0.935, 0.983
50.79, 39.95, 3.61
58.73, 51.55, 3.85
No. 5000
0.865, 0.922, 0.951
59.18, 48.01, 2.37
61.51, 23.71, 2.57
No. 6000
0.751, 0.801, 0.907
125.7, 99.85, 7.89
146.6, 82.13, 2.06
No. 7000
0.759, 0.846, 0.924
116.3, 108.4, 7.14
127.1, 82.67, 4.90
No. 8000
0.884, 0.900, 0.943
81.25, 35.46, 5.42
82.47, 67.70, 2.88Judging 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 Section3.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.



---

*Source: 432613-2014-01-16.xml*