g k =Kk DB k M k Z + h 0Jc I + n k (2) 
where I is the N^xl unit vector, h l k and \ k are 
respectively the gain and offset of the photometric parameters 
which balance the relative effects of sun zenith angle and 
atmosphere condition between different observed images. 
3. JOINT MAP RECONSTRUCTION MODEL 
It has been proved that the image registration (or motion 
estimation) and SR reconstruction are interdependent (Hardie, 
Barnard et al. 1997; Segall, Molina et al. 2003; Shen, Zhang et 
al. 2007), and a desirable solution is to simultaneously 
implement the image registration and SR reconstruction in a 
joint framework. Let the full set of P LR images, registration 
parameter vectors be denoted by g = {^,£2>••••£/>} » 
m = {m v m^,....m P } respectively. Employing the maximum a 
posteriori (MAP) framework, the desired HR image and 
registration parameters can be computed as 
z,m = arg max {p(z, m | g)} (3) 
Z,m 
Applying Bayes’ rule yields, (3) becomes: 
p(g\z,m)p(z,m)} 
P(S) J (4 ) 
P(g\z,m)p(m)p(z) 
P(g) 
Since p(z,m\g) is independent of g , equation (4) can be 
rewritten as 
z,« = aigmaxjj~J[/>(& \z,m k )p{m k )]p{z)^ (5) 
Assuming the model noise in (2) to be Gaussian with mean zero 
and variancecr 2 , the total probability function p(g k \z,m k ) can 
be given by 
P(g k \z,m k ) = ~e\p< 
\ k DB k M k Z h 0k I 
2(7 2 
- arg max 
z,nt 
z,m - arg max 
z,m 
(6) 
where O k is a diagonal matrix that denotes which pixels are 
outliers that strongly deviate from the observation model. In 
order to preserve the edge and detailed information in the 
reconstructed image, a Huber-Markov prior is employed 
for p(z). The Huber-Markov prior has the form as (Schultz and 
Stevenson 1996) 
P(z) = — exp 
P x,y ceC 
(7) 
where C is a constant value, ¡5 is a control parameter, C is a 
clique within the set of all image cliques C , the quantity 
^c( z x ,y) ls a spatial activity measure to pixel Z which is 
x,y 
often formed by first-order or second-order differences, 
a nd /?(■) is the potential Huber function which is defined as 
P( 0 = 
\ï 1*1 ^ P 
\2p\i\-/S \i\>ju 
(8) 
where p is a threshold parameter separating the quadratic and 
linear regions (Schultz and Stevenson 1996). As for 
the d c (z x y ) , we compute the following finite second-order 
differences in four adjacent cliques for every location (x,>')in 
the SR image 
d c ( Z x,y) Z x-l,y ^- Z x,y +Z x+l,y 
(9) 
d c ( Z X,y) = Z x,y-1 ~ ^ Z x,y + z x,y+l 
(10) 
d c ( Z x,y ) ~ ~Jj\_ Z x-\,y-\ ~ ^ Z x,y + Z x+1,>-+1 J 
(11) 
d c ( Z x,y ) = ^J^\_ Z x-l,y+l - ^ Z x,y + Z x+l,y-l] 
(12) 
Generally, the choice of the prior statistical model 
registration parameters p(m k ) is highly application 
for the 
specific 
(Hardie, Barnard et al. 1997), and it depends on the assumed 
motion model. In this paper, we assume the following global 
affine motion model between the observed images 
x 2 =a 0 + a x x x + a 2 y x 
(13) 
y 2 =K + b \X\ +b 2 yi 
(14) 
where (Xj, ) and (x 2 , y 2 ) are respectively the pixel locations 
of the referenced and un-referenced images, 
anda 0 ,a [ ,a 2 ,bQ,b ] ,b 2 are the affine parameters. In this case, the 
registration parameters are significantly overdetermined by the 
data and a prior is not necessary to yield a useful solution 
(Hardie, Barnard et al. 1997). Therefore, we select ppt^) as a 
constant value in this paper. Thus, substituting (6) and (7) in (5), 
after some manipulations, the following minimization cost 
function is obtained. 
¿map = ar g m in[£(z,/«,0)] (15) 
with the objective function E(z,m,0) is 
E(z,m,0) = £|| O t ( gi -l ht DB l ,M k z-KJ)t 
* (16) 
+aXIM(a.,)) 
x,y ceC 
It should be noted that the motion matrix M k embodies the 
registration parameter vector m k . 
4. OPTIMIZATION METHOD 
Figure 1. The optimization procedure 
In the cost function (16), there are three large sets of unknowns: 
the HR image z , the registration parameters m k and the