1
Comparing SAR based short time-lag
cross-correlation and Doppler derived sea ice drift
velocities
Thomas Kræmer, Harald Johnsen, Camilla Brekke, Geir Engen
Abstract—This paper shows initial results from estimating
Doppler radial surface velocities over Arctic sea ice using the
Sentinel-1A satellite. Our study presents the first quantitative
comparison between ice drift derived from the Doppler-shifts
and drift derived using time series methods over comparable time
scales. We compare the Doppler-derived ice velocities to global
positioning system tracks from a drifting ice station as well as
vector fields derived using traditional cross-correlation between
a pair of Sentinel-1A and Radarsat-2 images with a time lag of
only 25 minutes. A strategy is provided for precise calibration
of the Doppler values in the context of the the Sentinel-1A level-
2 ocean radial surface velocity product. When comparing the
two methods, root mean squared errors (RMSEs) of 7 cm/s were
found for the EW4 and EW5 swaths while the highest RMSE of
32 cm/s was obtained for the EW1 swath. Though the agreement
is not perfect, our experiment demonstrates that the Doppler
technique is capable of measuring a signal from the ice if the ice
is fast moving. However, for typical ice speeds, the uncertainties
quickly grow beyond the speeds we are trying to measure. Finally,
we show how the application of an antenna pattern correction
reduces a bias in the estimated Doppler offsets.
I. INTRODUCTION
With the launch of the Sentinel-1 synthetic aperture radar
(SAR) satellites by the European space agency (ESA), new
possibilities have emerged for monitoring sea ice motion
from space. Through precise estimates of the azimuth (along-
track) center frequency or Doppler centroid it is possible
to obtain a near-instantaneous measurement of the motion
of surface scatterers parallel to the pointing direction of the
radar antenna. Following the Sentinel-1 level-2 ocean (OCN)
product naming convention this line-of-sight (LOS) speed is
referred to as a radial surface velocity (RVL) [1]. Doppler-
derived radial surface velocities were originally studied in
the context of ocean wind and surface current retrieval [2]
and later demonstrated over sea ice [3], [4] using data from
the advanced synthetic aperture radar (ASAR) instrument
on-board the Envisat satellite with encouraging results. The
ASAR instrument was not designed with this product in mind,
however, and results using stripmap data were degraded by
antenna gain problems [4], [5]. In constrast, the Sentinel
satellites are constructed to provide very fine control over the
antenna, orbit and attitude enabling very high measurement
precision.
Ice buoys with high time resolution are the most natural
source of data for calibration and validation of radial surface
velocities over ice, but the spatial coverage of drifting ice
buoys is low as they only provide point measurements. An at-
tractive alternative is to use drift fields derived from frequently
available satellite data, where a single SAR scene may cover
a swath width in the order of 400–500 km. Two-dimensional
ice drift fields are regularly estimated by cross-correlating
similar image patterns between pairs of satellite images [6],
[7]. A problem with comparing Doppler measurements derived
from a single image to cross-correlation (CC) drift vectors
derived from a pair of images is that the time separation
between scenes can be large (traditionally 1–4 days). Over
such time scales, the motion of the ice may be highly non-
linear, which prevents a direct quantitive comparison between
the two methods [3], [4]. However, due to the increased
number of SAR missions orbiting the Earth, we are now at
a point where images from multiple satellites can be used
together to reduce the time separation between acquisitions to
minutes and hours rather than days. This takes us closer to a
valid assumption of linear drift between the scenes.
In this paper we estimate the 2-D ice displacement field
between a pair of scenes from the Radarsat-2 (RS2) and
Sentinel-1A (S1A) satellites with a time spacing of only 25
minutes. We use global positioning system (GPS) tracks from
a drifting ice station to check that the ice movement was
approximately linear between scenes. The derived CC drift
field is then projected onto the antenna LOS and compared
with the Doppler RVL drift showing good agreement. We
present the S1A RVL product in the context of S1A extra
wide swath (EW) mode data, but the algorithm is general and
can be applied to any appropriately prepared SAR data (see
section III-A for details).
The paper is organized as follows. We first introduce the
data selected for the experiments in section II. Section III
provides the theoretical background on Doppler frequency esti-
mation from S1A data. Section IV details the cross-correlation
algorithm. Results obtained using the two methods are then
compared in Section V. Section VI summarizes our findings
and provides recommendations for futures studies of RVL for
sea ice drift measurements.
II. DATASET
For our experiment we used three sources of ice drift
information; GPS positions from a drifting ice station, Doppler
derived velocities and cross-correlation displacement measure-
ments. The following gives a brief overview of the dataset.
1) Drifting ice station: In the first half of 2015 the Norwe-
gian Polar Institute (NPI) conducted the Norwegian Young
Sea Ice Cruise (N-ICE) whose objective was to increase
2
Table I: Synthetic aperture radar data used in the experiment
Scene Sensor
Sensing time [UTC]
Mode Polarization
Pass
direction
Angle of
incidence
Nominal looks
(Ra x Az)
S1 SENTINEL-1A 2015-03-16 06:05:34 Extra Wide HH / HV Descending 19–47
◦
6 x 2
R1 RADARSAT-2 2015-03-16 06:30:13 ScanSAR Wide HH / HV Descending 20–49
◦
4 x 2
S2 SENTINEL-1A 2015-12-27 10:03:58 Extra Wide HH / HV Ascending 19–47
◦
6 x 2
S3 SENTINEL-1A 2016-01-12 09:55:46 Extra Wide VV / VH Descending 19–47
◦
6 x 2
understanding of the effects of decreasing ice thickness on
ice dynamics, energy fluxes and associated local and global
climate variables [8]. In late December 2014 the research
vessel Lance was frozen into an ice floe North of Svalbard
to become an ice station passively flowing with the drifting
sea ice towards Fram Strait. When the floe broke up or the ship
exited the ice the ship moved back into the ice to freeze into
another ice floe. The ship continually logged its GPS position
with 10 s intervals which we use to check displacement fields
derived using CC.
2) Cross-correlation drift: Two-dimensional ice drift vector
fields are regularly estimated from pairs of remote sensing
images by cross-correlating image patches (technical details
given in section IV). For this purpose, operational services
normally prefer spaceborne sensors with wide geographical
coverage [7]. Traditionally passive microwave instruments
have provided rapid revisit times and wide coverage, but with
poor resolution in the order of several kilometers. SAR sensors
provide a good compromise between wide coverage and high
resolution by electronic steering of the antenna in elevation,
which periodically illuminates a set of swaths. For RS2 this
normally means the ScanSAR Wide (SCW) mode which
covers an area of ∼500 km x 500 km with a square ground
range pixel spacing of 50 m and a resolution of ∼100 m in
each dimension [9]. This mode uses four beams which cover
incidence angles ranging from 20 to 49 degrees. For S1A, the
EW medium resolution product covers an area of ∼400 km x
400 km with a square pixel spacing of 40 m and a resolution
of ∼90 m in each dimension [1]. EW mode images are
acquired using five beams (EW1–EW5) with incidence angles
in the range 19–47 degrees. In contrast with the ScanSAR
mode used on RS2, the S1A EW product implements the
terrain observation by progressive scans (TOPS) mode which
electronically sweeps the antenna in azimuth in addition to
stepping in elevation.
From archives of RS2 and S1A scenes we selected data
based on the following requirements: A spatial overlap of at
least 40% between images was desired to obtain a reasonably
large 2-D drift field using CC. At the same time we wanted
the time spacing between images to be as small as possible,
while still allowing the ice to be displaced sufficiently to be
measured by pattern matching. Furthermore, the search was
limited to image pairs where the ship was located within both
scenes which allows comparison with GPS positions. Because
sea ice drift speeds are small compared to surface wind speeds
we also wanted to have reasonably high drift speeds to increase
chances of having detectable Doppler shifts. Therefore, we
also included the ice speed (estimated from the ship’s GPS) in
the search. The image pair (S1, R1) in Table I stood out as an
excellent candidate. The other scenes are used for calibration
investigations.
3) Doppler-derived drift: The Doppler estimation algo-
rithm requires a full-bandwidth processed single-look complex
(SLC) image as input. This is not a standard S1A product and
we therefore require that the raw unfocused (Level-0) data
is available so a custom SLC can be created without using
window functions, thereby retaining the full bandwidth of the
data. Raw data was not available for RS2 and hence Doppler-
derived velocities were calculated for S1A scenes only. We
therefore use Doppler anomalies from homogeneous parts of
S2 to calibrate the Doppler anomalies in S1.
All the scenes had two polarization channels; horizon-
tal transmit/horizontal receive (HH) and horizontal trans-
mit/vertical receive (HV). In our experiments we have focused
on the HH polarization only. This has long been the preferred
channel for many sea ice applications, however the algorithms
are not limited to use a particular polarization. For CC drift
estimation it has been shown that use of both channels may
be beneficial [7]. Although the Doppler estimation algorithm
presented in section III does not assume a particular polariza-
tion, it should be noted that the signal-to-noise ratio (SNR)
over ice and water is often lower compared to HH which will
lead to larger uncertainties in the RVL estimates.
III. DOPPLER-DERIVED RADIAL SURFACE VELOCITY
MEASUREMENT AND CALIBRATION
The SAR imaging process can be formulated as a convo-
lution of the transmitted signal modulation with the ground
reflectivity, weighted by the antenna directivity pattern [10].
High resolution is achieved by pulse compression in both the
across track and along track direction by proper modeling
and matched filtering of the target phase history φ = −2kR
[rad] where R is the sensor–target range, k = 2π/λ is the
wavenumber and λ is the wavelength. As the satellite moves,
the relative range between the antenna and the ground changes
at a rate
˙
R, introducing Doppler shifts in the signal. The
angular Doppler centroid in [rad/s]
$
dc
= −2k
˙
R(τ
0
) = −2kv
r el
·
ˆ
r (1)
is the frequency offset corresponding to the time τ
0
when the
target is in the beam center. Here, v
r el
= v
t
−v
s
is the relative
velocity between the sensor (v
s
) and target (v
t
) and
ˆ
r is a unit
vector pointing from the sensor to the target. Note that only the
projection of the relative velocity onto the LOS matters, and
hence any along track target motion is not observable using
Doppler measurements.
3
Angular frequencies in [rad/s] are related to their linear
counterpart in Hz by a factor 2π. The pulsed nature of the
SAR system limits the observable $
dc
values to the baseband
region [−$
prf
/2, $
prf
/2], where $
prf
is the angular pulse
repetition frequency (PRF). In general the Doppler centroid,
$
0
dc
, may be expressed as $
0
dc
= $
dc
+ M $
prf
where
$
dc
is the fractional Doppler centroid and M is an integer
referred to as the Doppler ambiguity [10]. It is common for
SAR satellites to follow a yaw-steering law which adjusts the
antenna pointing direction as a function of latitude to provide
M = 0. In this paper we therefore only consider the fractional
part ($
0
dc
= $
dc
), but in general the ambiguity would have
to be estimated as well. Estimation of both $
dc
and M is
covered in standard texts on SAR such as [10, ch. 12]).
At any position we can model the measured Doppler cen-
troid as a linear combination of (i) Doppler shifts due to the
relative sensor–target motion as predicted by eq. (1) and (ii)
Doppler shifts due to antenna effects,
$
dc
≈ $
geom
+ $
phys
+ $
em
(2)
where $
geom
is the contribution due the relative motion
between the satellite and a stationary target on the surface
of the rotating Earth, $
phys
is the geophysical Doppler shift
due to the LOS motion of surface scatterers relative to the
rotating Earth and $
em
is a bias introduced by the antenna
electronic mispointing [11].
1) Geometric Doppler, $
geom
: The geometric contribution
can be calculated by solving the range–Doppler equations
taking into account the sensor attitude. This is explained in
great detail in standard textbooks on SAR (see e.g. [10, chap.
12]). S1A uses total zero-Doppler steering which combines
yaw-steering with an additional pitch-steering to provide a
nominally zero geometric Doppler ($
geom
≈ 0) across the
entire swath [12].
2) Electronic mispointing, $
em
: The most commonly used
antennas for spaceborne SARs are phased array systems which
have the ability to electronically steer the beam in both
azimuth and elevation as well as the freedom to shape the
antenna pattern by varying the amplitude and phase of each
transmit/receive module (TRM). Over time, the characteristics
of the TRMs change due to drift in the electronics or physical
damage to the antenna. For a given elevation angle, these
deviations may cause the maximum gain to occur at an
azimuth angle slightly offset from the nominal pointing angle
which introduces an unintentional squint. The effect, known
as electronic mispointing, contributes to an offset $
em
which
is a function of elevation angle only. If the embedded row
patterns, error matrix and excitation coefficients are available,
the full antenna pattern can be simulated and the mispoint-
ing estimated directly using the antenna model presented in
[13]. However, publicly available auxiliary calibration files
for S1A (referred to as AUX_CAL in ESA documentation
[14]) only provide two slices of the antenna and this dyadic
approximation hence does not capture the range variation of
the mispointing. Alternatively $
em
can be estimated from
data over stationary areas of homogeneous backscatter. This
is discussed further in section V where mispointing profiles
predicted by the antenna model are compared with estimates
from rainforest data.
3) Geophysical Doppler, $
phys
: By explaining away con-
tributions from the motion of the Earth and antenna effects we
can invert eq. (2) to obtain the parameter of interest, $
phys
;
the Doppler shift due to the geophysical motion of scatterers
on the surface:
$
phys
≈ $
dc
− $
geom
− $
em
. (3)
Referring to eqs. (1) and (3), the target LOS speed u
r
can be
calculated as
u
r
= −
$
phys
2k
. (4)
In [2], Chapron et al. interpreted the target speed u
r
using a
simple geometrical model,
u
r
=
h(u
g
sin θ
i
− u
v
cos θ
i
)σ(θ
0
i
)i
hσ(θ
0
i
)i
(5)
where h·i denotes ensemble averaging over the local incidence
angles θ
0
i
, θ
i
is the angle of incidence at the center of the
estimation cell, u
g
is the target speed tangential to the surface,
u
v
is the target speed normal to the surface, and σ is the signal
intensity. Following [3], assuming homogeneous backscatter
and no vertical motion in the central ice pack (u
v
= 0), we
can approximate the ground range surface velocity u
g
as
u
g
≈
u
r
sin θ
i
. (6)
This is of course not always a good assumption, and we
discuss this further in section V. The rest of this section details
the baseband Doppler centroid estimation algorithm.
A. Preprocessing
SLC data for the EW mode is not a standard product
delivered by ESA. We therefore take the unfocused raw data
(Level-0) as our starting point and focus each burst to an SLC
using an ω–K algorithm adapted for TOPS mode data [15].
No window functions are applied during focusing, thereby
retaining the full bandwidth of the data.
Let I
(m)
(t) denote the m-th complex-valued SLC burst
where t = (t, τ) denotes range time (t) and azimuth time (τ ).
While aquiring each burst, the TOPS mode sweeps the antenna
in azimuth, changing the beam center angle at a rate k
ψ
[rad/s]. The one-to-one relationship between the beam center
angle and the Doppler centroid means that a linear sweep of
the antenna introduces an approximately linear change in the
Doppler centroid at a rate k
a
[Hz/s] given by
k
a
≈ −
2v
s
λ
k
ψ
where v
s
is the platform speed. This phase ramp in the SLC
data must be removed before traditional Doppler estimators
can be used [16]. A deramped SLC burst, I
(m)
d
(t), can be
produced by multiplying the burst with a chirp
I
(m)
d
(t) = I
(m)
(t)e
−jk
t
(τ −τ
(m)
c
)
2
(7)
where j =
√
−1, k
t
= k
a
/α and τ
(m)
c
denotes the azimuth
time of the m-th burst center.
4
As described in [16], the range-dependent factor α given by
α = 1 +
k
a
|k
r
|
(8)
can be interpreted as an antenna scaling factor relating a
physical antenna operating in TOPS mode to a mathematically
equivalent scaled antenna operating in stripmap mode. Here,
k
r
≈ −
2v
2
eff
λr
0
(9)
is the range-dependent Doppler rate, v
eff
=
√
v
s
v
g
is the
range-varying effective sensor speed (see e.g. [10, p. 127]),
v
g
is the gound speed of the antenna footprint, and r
0
is the
range at the time of closest approach.
After deramping the bursts are merged onto a common
grid to obtain a connected SLC image I(t) per swath using
the procedure outlined in appendix A. From each swath, the
Doppler centroid is estimated blockwise on a regular grid
where the block size is 298×228 pixels in range and azimuth
direction (∼ 4.2 × 9.5 km), respectively, and the step size is
25% of block side lengths. The Doppler estimation procedure
is described in the following section.
B. Doppler estimation and side-band correction
Doppler centroid estimators exploit the observation that the
azimuth power spectrum of the data, P ($; t), is related to
the two-way azimuth antenna directivity pattern D($), where
$ denotes the azimuth-direction (Doppler) frequency [17].
The observed power spectrum estimated from a block of data
centered at time t is given by
P ($; t) =
Z
dt
0
Z
dτ
0
I(t
0
)h(t
0
− t)e
−j$τ
0
2
(10)
where h(t
0
) = h
r
(t
0
)h
a
(τ
0
) is a dyadically constructed win-
dow function satisfying
Z
h
r
(t
0
)dt
0
=
Z
h
a
(τ
0
)dτ
0
= 1 and
the integration is over a data block of size (B
r
, B
a
) in range
and azimuth, respectively. Estimation of the Doppler centroid
therefore amounts to finding the Doppler frequency $
dc
that
provides the best fit between the observed power spectrum
and a model spectrum based on the expected antenna pattern
described below.
The antenna pattern extends well beyond the PRF which
means that energy in the side bands, i.e., signal components
outside the main band [−$
prf
/2, $
prf
/2] will fold (alias) into
the main band resulting in azimuth ambiguities (also referred
to as ghost images) [10]. This is equivalent to energy from
neighboring geographical areas influencing the spectrum of the
estimation area. The degree of aliasing depends on the PRF
location relative to the antenna pattern, but typically there is
only significant energy coming from the first side band. In
the case of high SNR and homogeneous data, this aliasing
process results in a power spectrum that is well modeled by
a raised cosine [18]. This is the motivation for time-domain
estimators which exploit the Fourier pair relationship between
the autocorrelation function and the power spectrum (Wiener–
Khinchin theorem) [18]. In inhomogeneous areas, backscatter
registered through the side bands of the antenna may contribute
significantly to the total power. This typically happens close
to shore, where the beam center covers the ocean which may
have low backscatter, while the side bands cover neighboring
mountain areas with high backscatter. A consequence of the
pulsed operation of the SAR system is that the spectrum is
periodic. For a stripmap system, this periodicity follows the
raw data PRF, $
prf
, but for TOPS we need to take into
account the antenna scaling factor α making the data periodic
with a separation of $
∆
= $
prf
/α.
Building on work by Madsen [18] and Bamler [19], En-
gen and Johnsen [20] therefore modeled the expected power
spectrum P ($; t) as the sum of the frequency folded antenna
pattern weigted by the average intensity σ(t
l
) within each side
band and a white noise component b(t) capturing the thermal
and quantization noise.
By introducing the normalized azimuth frequency η =
$/$
∆
, the model for the azimuth spectrum can be expressed
as
P (η; t) = W (η − ς
0
)
"
X
l
σ(t
l
)D
l
(η −ς(t
l
)) + b(t)
#
(11)
where
D
l
(η) = D((η + l)$
prf
) (12)
is the two-way antenna gain pattern of the l-th side band, W is
a window covering the critical bandwidth $
∆
, centered on ς
0
,
the normalized Doppler centroid used during focusing and ς
is the normalized Doppler centroid we wish to estimate. The
time t
l
= t + ∆t
l
is the position of the l-th ghost image
where the range component of ∆t
l
is the range migration and
the azimuth component is l · f
prf
/k
r
(∆t
l
= 0 for l = 0).
The corresponding autocorrelation coefficients (index by n)
given by
p
n
(t) =
Z
dηP(η; t)e
j2πnη
(13)
are then
p
n
(t) =
X
l
f
n
(t
l
)d
n
(l + ∆l(t
l
)) + b(t)δ
n
(14)
where f
n
(t) = σ(t)e
−j2πnς
(t) are the side-band corrected
autocorrelation coefficients, ∆l(t) = ς
0
− ς(t) is the offset
between the Doppler centroid used during focusing and the
true Doppler centroid, δ is the Kroenecker delta function (δ
0
=
1, δ
n6=0
= 0) and
d
n
(λ) =
Z
dηD(η)W (η − λ)e
j2πnη
. (15)
The coefficients f
1
and f
0
are of special interest as the first
coefficient provides an estimate of the average signal intensity
and the true Doppler centroid, and the zeroth order coefficient
can be used to estimate the additive noise level. We can invert
eq. (14) solving for f
1
by introducing
(t) =
X
l
f
1
(t
l
)(d
1
(l + ∆l(t
l
)) − d
1
(l)) (16)
expressing the first correlation coefficient as
p
1
(t) =
X
l
f
1
(t
l
)d
1
(l) + (t) (17)
5
Figure 1: Effect of antenna element pattern (AEP) correction
of raw data on Doppler centroids estimated from the EW3
swath of the S3 rainforest scene. If the correction is not
applied, the estimated Doppler centroids contain an approx-
imately linear trend in azimuth within each burst.
and exploiting the Fourier shift property
˜p
1
(ω) =
X
l
˜
f
1
(ω)e
−jω·∆t
l
d
1
(l) + ˜(ω) (18)
from which we can obtain
f
1
(t) = F
−1
˜p
1
(ω) − ˜(ω)
P
l
d
1
(l)e
−jω·∆t
l
(19)
Here, ˜x(ω) = F{x(t)} denotes the Fourier transform of x(t)
with corresponding inverse transform x(t) = F
−1
{˜x(ω)} and
the implicit dependency on f
1
through can be solved by fix
point iteration of eq. (19) starting with the assumption that
ς = ς
0
( = 0). The Doppler centroid is then $
dc
= ς$
prf
.
If the described Doppler estimation algorithm is applied as
stated to the SLC data a trend in the Doppler can be observed
within each burst (see fig. 1). This can be explained by the
elevation direction antenna element pattern (AEP) envelope
which weighs the total phased array beam pattern and biasing
the beam center slightly. The effect can be mitigated by either
calculating the resulting Doppler offset and including another
correction term in eq. (5) or by dividing the raw data by the
element pattern to flatten the data before Doppler estimation.
Antenna pattern profiles are provided in the AUX_CAL auxil-
iary files made available at https://qc.sentinel1.eo.esa.int/. Both
options provide trade-offs. On the one hand, post-processing
the estimated Doppler shifts is still sensitive to the small
radiometric discontinuity at burst overlaps which can introduce
a scalloping pattern with significant harmonics in the estimated
Doppler values. On the other hand, applying a gain correction
to the data will color the noise and is therefore strictly not
in agreement with the proposed model which assumes white
noise. However, the gain correction is very small and exper-
iments where both methods were tested with rain forest data
showed that the gain correction method gave similar results to
the Doppler post-processing method, while providing cleaner
estimates in burst overlap zones. We therefore recommend the
gain correction method. Uncertainties related to calibration of
the Doppler offsets are further discussed in section V.
IV. CROSS-CORRELATION DRIFT ESTIMATION
Motion estimation algorithms are often categorized into
pixel based and feature based algorithms. Pixel based algo-
rithms (e.g., [7], [21], [22]) use the pixel information directly
to maximize a measure of similarity between images, while
feature methods (e.g., [23]–[25]) first detect interest points and
match derived features. Among pixel based methods, cross-
correlation algorithms are the most popular and they have
been used for a long time for motion estimation from SAR
[6]. Using drifting ice buoys as reference, studies on the
accuracy of CC methods with SAR have reported root mean
squared error (RMSE) values as low as 300 m when using
buoy data [22]. Hollands and Dierking [26] obtained RMSE
values in the order of 400-560 m using manually drawn vectors
as a reference. For comparison, a feature matching method
was recently presented with a reported RMSE of 202 meters
when compared to manually drawn vectors [25]. However, the
accuracy will vary depending on the time separation between
the images as longer time separation increases the chance of
image pattern decorrelation. Thus, it is expected that studies
combining multiple satellites like Sentinel-1 A and B will
perform well just due to the increased time resolution. The
specific algorithm used in our example is as follows.
Given two detected images A
1
and A
2
and a set of lati-
tude/longitude positions (θ
lat
(k), θ
lon
(k)), k = 1, . . . , K, we
find the pixel corresponding to the geographical reference
point in each of the two images, p
1
(k) = (x
1
(k), y
1
(k))
and p
2
(k) = (x
2
(k), y
2
(k)). Around each point we extract
a square block of data with side length w = 129 (∼ 6.5 km)
and compute the normalized cross-correlation (NCC)
ρ(s
x
, s
y
) =
X
x
X
y
b
1
(x, y)b
2
(x + s
x
, y + s
y
) (20)
where the image blocks b
1
and b
2
have been normalized by
subtracting the mean and dividing by the standard deviation.
The offset that maximizes the correlation is taken as an
estimate of the displacement of the reference point in pixels
ˆ
s = (ˆs
x
, ˆs
y
) = arg max
(s
x
,s
y
)
ρ(s
x
, s
y
) (21)
where coordinates are relative to the block center. Thus, if we
extract a block from A
2
centered on p
0
2
= p
2
+
ˆ
s it should look
similar to b
1
. Blockwise correlation is repeated independently
for each geographical refrence point, forming a collection of
point correspondences {(p
1
(k), p
0
2
(k))}
K
k=1
.
The algorithm blindly maximizes the correlation which
may be low over e.g. open ocean. Therefore, some of the
estimated vectors will likely be incorrect. It is common to
reduce incorrect vectors by thresholding the NCC under the
assumption that low NCC values indicate incorrect matches.
In our case all vectors with an NCC value of less than 0.3
were discarded. However, this is often not enough to filter out
all incorrect matches. Therefore, several algorithms employ a
two-pass strategy where the first pass matches reference points
from A
1
to A
2
and the second pass takes the matched points