Virtual radar ice buoys – a method for measuring fine-scale sea ice drift

Here we present an algorithm for continuous ice drift estimation based on coastal and ship radar data. The ice drift is estimated for automatically selected ice targets in the images. These targets are here called virtual buoys (VBs) and are tracked based on an optical flow method. To maintain continuous ice drift tracking new VBs are added after a given number of VBs have been lost; i.e. they can not be tracked reliably any more. Here we also apply the algorithm to data of three test cases to demonstrate its capabilities and properties. Two of these cases use coastal radar data and one ship radar data. Ice drift velocity and direction information derived from the VB motion are computed and compared to the prevailing ice and weather conditions. Also a quantity measuring the local divergence or convergence is computed for some VBs to demonstrate the capability to estimate derived kinematic sea ice parameters from VB location time series. The results produced by the algorithm can be used as an input for estimation of the dynamic properties of sea the ice field, such as ice divergence or convergence, shear, vorticity, and total deformation.


Introduction
Sea ice motion is an important parameter, because ice dynamics has a major effect on the nature of sea ice.Ice motion can cause ice pressure, which in turn contributes to ice deformation, or diverging ice motion can cause opening of ice (cracks, leads).The main goal of this study was to develop and test an algorithm suitable for continuous operational ice drift monitoring based on radar data and to demonstrate it for a few test data sets from coastal and ship radars.
If necessary, radar platform (ship) motion can be compensated based on the geoinformation (GPS (Global Positioning System) position of each radar frame).However, large ship motion is not desirable because the radar signal is rapidly attenuated as a function of the range or, in the worst case, the two radar images with a given time gap between them are not overlapping any more.If ship motion is not compensated then the ice drift with respect to the ship is measured, like in the ship radar test case in this study.
Many publications on ice drift from Synthetic Aperture Radar (SAR) imagery and other satellite-borne Earth Observation (EO) data have been published.The motion estimation is based on detecting the same features in two adjacent images (e.g.Fily and Rothrock, 1987;Kwok et al, 1990;Sun, 1994;Thomas et al., 2004;Haarpaintner, 2006;Thomas et al., 2008;Karvonen, 2012).Motion vector estimation for weather radar data has also been studied, e.g. in Peura and Hohti (2004).Sea ice drift and object tracking from coastal radars have been studied earlier: e.g. in Okhotsk Sea in Tabata et al. (1980), by matching of prominent features preserved from image to another; in the Chukchi Sea near Barrow, Alaska, in MV et al. (2013), using Lucas-Kanade optical flow algorithm (Lucas and Kanade, 1981) for features detected by Harris corner/edge detection (Harris and Stephens, 1988) algorithm; and in the Baltic Sea in Karvonen (2013a), using a combination of phase correlation and normalized cross-correlation.Coastal radar data for sea ice analysis in the Canadian Arctic have been utilized in Shapiro (1976), Shapiro and Metzner (1989), Mahoney et al. (2007), Druckenmiller et al. (2009), and Mahoney et al. (2015).
One practical restriction for all the methods is that they can only operate where distinguishable objects exist, and in the featureless areas either no estimates are given or the es-Published by Copernicus Publications on behalf of the European Geosciences Union.
timates are inter/extrapolated values.The motion detection methods are computationally intensive but in principle easy to parallelize due to the local nature of the computations.
In our earlier studies (Karvonen, 2013a;Karvonen et al., 2013) we have used data with longer temporal differences (10-30 min) than here and an algorithm based on crosscorrelation techniques.This approach was directly adapted from techniques used for ice drift estimation from SAR imagery (Thomas et al., 2004(Thomas et al., , 2008;;Karvonen, 2012).However, with a finer temporal resolution sub-pixel resolution would be desirable, and it can not be achieved using crosscorrelation techniques without suitable interpolation of the results.This interpolation needs to be nonlinear, thus also requiring a larger support area than a simple linear interpolation.Optical flow algorithm in turn inherently gives the estimates with a sub-pixel resolution.For this reason we found the optical flow approach better suitable for sea ice drift analysis from radar imagery of fine temporal resolution.Compared to our earlier algorithms (Karvonen, 2013a;Karvonen et al., 2013) we also get a better (sub-pixel) spatial and temporal resolution of the ice drift estimates.
This optical flow algorithm also enables continuous operational ice drift monitoring.The algorithm automatically adds virtual buoys (VBs) as their number decreases under a given number and the monitoring can be continued for the ice season without any human intervention.Based on the automatically performed ice drift analysis near-real-time information, e.g. on local convergence or divergence (closing or opening ship track), can be delivered for navigation in the area.Also the radar image data can automatically be stored and transmitted with a shorter time interval when the ice is moving or deforming for further analysis.During periods of no ice drift or insignificantly small ice drift radar, image data storage or transmission is not necessary.

Radar image capturing and transmission
Marine radars typically operate at 10 GHz (wavelength λ ≈ 3 cm) and 3 GHz (λ ≈ 10 cm), i.e. at X and S bands.The radar resolution is defined by means of two resolutions: bearing and range resolutions.Bearing (a.k.a. as azimuth or angular) resolution is the ability of a radar system to separate objects at the same range but at slightly different bearings.The bearing resolution depends on radar beam width and the range of the targets.Range resolution is the ability of a radar system to distinguish between two or more targets on the same bearing but at different ranges.Range resolution depends on the radar pulse length, unless radar pulse compression techniques (Cohen, 1987) are used.
There are presently about 60 coastal radars along the Finnish coast, administered by the Finnish Traffic Agency.Also the Finnish navy and coastal guard have complementary coastal radar networks of their own.The radars are typically located 20-50 m above the sea level (the Tankar radar used in two cases presented here about 30 m).The environmental monitoring can be realized by instrumenting the radar with an independent radar server, in our case the Image Soft radar server designed and manufactured by the Finnish company Image Soft Ltd.We presently have radar servers in Marjaniemi and Tankar coastal radars in the Bay of Bothnia and in the Uto coastal radar facing the northern Baltic Proper.Further installations have been planned.The radar server is a LINUX server equipped with a radar image capturing hardware card.The radar servers capture the analog signals of the radar and rasterize a PPI (plan position indicator) image from the radar signal, the radar triggering pulse and the radar antenna pulse for each radar revolution (or an image per a user-defined time interval).PPI is the most common type of radar display: the radar antenna is represented in the centre of the display, so the distance from it can be presented as concentric circles.The sampling rate of the image digitization on the image server is 20 MHz.
The digital-analog conversion produces 12 bit raw radar data values.For the PPI imagery these data are quantized to 8 bits per pixel.According to our experience this quantization is adequate for sea ice tracking.Many of the radar server processing parameters can be adjusted for the user to be suitable for the radar and application; more details can be found in the radar server technical manual (Image Soft, 2014) All the rasterized images are stored on the server hard disk while a subset of preprocessed images are sent to FMI.The preprocessing performed on the radar servers is a temporal median filtering of 15-20 s.The data are transmitted via a stand-alone GSM (Global System for Mobile Communications) mobile link, for details see Fig. 1.Due to the limited bandwidth of the GSM modems a preprocessed image is also transmitted in every 2 min, which has proved to be a suitable time interval for continuous ice drift monitoring.Our plans are to install the presented VB tracking software on the radar servers, leading to a reduced amount of required data transmission.In practice only VB motion, or VB motion data with radar imagery in the case of significant ice motion, then needs to be transmitted.

Data sets used in the study
We have used three data sets to test the new algorithm.Two of the data sets were coastal radar data from the Tankar coastal radar, located at 63.95 • N, 22.84 • E: the first data set period was 25 February 2011 from 03:00 to 16:58 UTC (total time period of about 14 h), and the second data set period was 8 February 2012 from 00:00 to 23:58 UTC (total time period was about 1 day).The February 2011 Tankar coastal radar data set will be referred as case A and the February 2012 Tankar coastal radar data set as case B. The data sets were selected such that they include significant ice drift.The temporal resolution, i.e time interval between two successive radar images of the two coastal radar data sets, was 2 min.The third data set was a longer period data set and collected on-board RV Lance during the period from 21 January 2015, 11:30 UTC, to 18 February 2015, 10:20 UTC.The time difference between successive radar images of the RV Lance data set was 10 min.For demonstration purposes we selected a 1-day period (8 February 2015 from 00:00 to 24:00 UTC) of the RV Lance data set.This data set will be referred as case C here.During this time period set, significant and heterogeneous ice motion with respect to the ship also occurred within the radar coverage.During this 1-day period the location of RV Lance was north of Svalbard, approximately 82.5 • N, 18 • E. The total range of the Tankar coastal radar data sets was 20 km and the image size was 1200 × 1200 pixels; i.e. the nominal resolution was 33.3 m.For the RV Lance radar the image size in pixels was the same as for the Tankar data, but the range was only 7.5 km, resulting to a nominal resolution of 12.5 m.The radar image area of the Tankar radar was shifted 10 km to the west and the images were rotated 49.93 • to the counterclockwise direction on the radar server with respect to the radar location.This kind of operations can be performed on the radar server by giving the displacement and rotation parameters.The translation and rotation were performed to minimize the land area in the imagery to get a larger cover of the sea area.
To reduce radar artifacts (e.g.due to weather phenomena, radar noise or clutter, electromagnetic interference) temporal median filtering in the beginning of each minute (11 images, corresponding to about same amount of time in seconds, assuming that the radar rotation frequency is about 1 Hz) was performed.During this time period the ice motion is neglectable and only the noise and possible artifacts are reduced by the filtering.The images were also processed by homomorphic filtering (HMF) to reduce the signal attenuation as a function of the range (Lensu et al., 2014).This processing mainly makes the visual analysis of the data eas-ier.According to some performed tests it does not have significant effect on the tracking of objects.
We also applied a rough land mask to the Tankar imagery.Our land mask was derived from the GSHHG (Global Self-Consistent Hierarchical High-resolution Geography database) from the NOAA (National Oceanic and Atmospheric Administration) coastline data (Wessel and Smith, 1996).
For the RV Lance data we were unable to perform the temporal median filtering because we only had data with 10 min temporal sampling at our disposal, i.e.only one unfiltered image every 10 min.

Weather and ice conditions at the test sites
The air temperature on 25 February 2011 (corresponding to case A) around the Tankar lighthouse and radar station was from about −15 • C in the morning to about −3 • C in the afternoon.The previous day was colder with a daily maximum temperature of about −15 • C. The wind direction was 150-180 • , and wind speed varied in the range 6-8 m s −1 .In the eastern parts of the area there was landfast ice; west of the fast ice there was a zone of very open ice (concentration 10-30 %), and west of this zone there were very close drift ice (concentration 90-100 %).The ice thickness in the area was 20-55 cm.The ice information were extracted from the FMI ice charts.
On 8 February 2012 (case B) the air temperature around the Tankar lighthouse and radar station was from −11 • in the morning to −20 • C in the evening; also, the previous day the temperatures were relatively cold, below −10 • C. The wind direction was 90-180 • , and wind speed in the range 2-6 m s −1 .According to the FMI ice charts there was a fast ice zone in the eastern parts of the area, a zone of new ice to north and east of the fast ice, and very close drift ice (ice concentration 90-100 %) further in the west.The ice thickness in the area was 5-30 cm.
On 8 February 2015 (case C) the air temperatures around RV Lance were cold, about −30 • C; the wind direction was 300-330 • and wind speed around 8 m s −1 .The ship was drifting with a speed of approximately 0.2 m s −1 , first to the south and later to the southeast.According to the met.norway ice charts the ice in the area was very close drift ice (ice concentration 90-100 %), and according to the operational Nansen Environmental and Remote Sensing Center (NERSC) Topaz ice model (Sakov et al., 2012) the ice thickness was 100-120 cm in the area.As the weather was cold, the opening leads were also frozen relatively fast and the ice thickness in frozen leads was typically around 20 cm.

Virtual buoys and tracking algorithm
In the first phase a filtering to reduce the signal attenuation as a function of the range is performed, which is described in more detail in Sect.5.1.After this we used an edge and corner detection to locate the VBs in the first radar image of a radar image sequence and also when adding new VBs after the number of VBs has reduced to a predefined level (an adjustable parameter).The tracking algorithm is also based on the optical flow (e.g.Horn and Schunck (1981) and Beauchemin and Barron (1995)) between successive image pairs.
The schematic flow diagram of our algorithm has been presented in the diagram of Fig. 2. In the first phase the VBs are initialized based on features (edges and corners) detected by using absolute local binary patterns (ALBPs; described in more detail in Sect.5.2), and then one iteration of motion tracking between the first two images of the image sequence is performed to prune the VBs due to noise amplification by the image filtering: all the images fed into the algorithm are first filtered by the HMF to reduce the range dependence of the radar signal.After this initialization a list of the automatically selected VBs (including the locations and crosscorrelations between the matched windows) are fed to the continuous VB tracking algorithm.At each iteration a new filtered radar image is put into the tracking system and the optical flow tracking between the previous and the novel image is performed.After each tracking iteration the number of resulting VBs is compared to a predefined threshold (T N ), and if the number of the remaining VBs is less than T N , new VBs are added starting from near range until the original number of VBs has been reached.The new VB location are defined based on ALBP features with the limitation that they are not allowed to be added closer than a given radius R x (we have applied R x = 15 pixels here, but this parameter can be defined by the user) from the existing VBs.After the VB adding step the tracking is continued with the updated list of VBs.If still more VBs exist than the defined threshold T N , then the VB tracking is continued without updating the list of VBs.

Radar image preprocessing by HMF
The temporal median filtering is applied already on the radar server before transmitting the data.After receiving the data HMF is applied to reduce the attenuation of the signal as a function of the distance from the radar.
HMF has its background in optical image processing (Pitas and Venetsanopoulos, 1990).HMF intensity I (r, c) at image location (r, c) for an optical image is presented as a product of the illumination L I and reflectance R I , where R I cab be considered as a quantity describing interesting objects in the scene and L I results from the lighting conditions, i.e. (1) The applicability of the method for coastal radar imagery is readily conceived, as the images make the impression that the radar tower "illuminates" the ice field where the sea ice features reflect this "light".To proceed, a logarithm transform is first applied to make the components of Eq. (1) additive and linearly separable.The variations of illumination are then considered as noise that can be reduced by applying a highpass filter.We assume that the reflectance, which is of interest, is represented more by the high-frequency components while the low-frequency components relate to the illumination.In other words, the lighting condition is assumed to vary slowly across the image (the physical attenuation of the radar power) and reflectance is changing faster (due to deformed ice areas such as ridges).In practice the frequency filtering has here been implemented by applying the 2-D fast Fourier transform (FFT) to the image, then performing the high-pass filtering in the Fourier domain, and finally performing the inverse FFT.Because FFT requires the size of the input to be a power of 2, we extend our image to the nearest power of 2 larger than the image size by mirroring with respect to the image boundaries (in our case a 1200 × 1200 pixel image is extended to 2048 × 2048 pixel image).In the frequency (FFT) domain the low-pass coefficients are attenuated by multiplying by a factor f < 1.0; in our case we have used f = 0.0, i.e. totally zeroing the low-frequency components.A schematic presentation of the HMF principle is shown in Fig. 3.An example of HMF for a radar image of the case.A time series is shown in Fig. 4.

Local binary patterns and VB selection
We have used LBPs (Ojala et al., 1996) here for the edge and corner detection.LBSs are computed locally over the image area around each image pixel.LBSs can be used as a tex-  ture measure or for detecting certain imagery features such as edges and corners.We used a step of /4 in direction corresponding to 8 bit binary patterns; the radius R LBP (distance from the centre pixel, see Fig. 5) used here was 2. A (8 bit) local binary pattern is defined as where g c is the grey tone of the centre pixel and the values of g k are the grey tones in the 8 pixels within the given radius R LBP around the centre pixel (see Fig. 5) and s 1 (x) defined as Here x is a generic argument for s 1 ; for LBP, x = g k −g c .We have used a variant which we here call the ALBP: and s 2 (x) defined as T is a threshold value and we have used value T = 10 in this study.Rotational invariance can be achieved by using the minimum among all the (eight) cyclic shifts of the ALBP.This is denoted here by ALBP r .The value 15 of ALBP r corresponds to an edge point, value 31 corresponds to a corner point, and value 63 corresponds to a sharp corner point.The ALBP 15 corresponds to the pattern of four adjacent points G i (i = 0. ..7) with the value of 1 (see Fig. 5), 31 to a pattern of five adjacent points G i with the value of 1, value 63 to a pattern of six adjacent points with a value of 1, and all the other values G i of 0. We also first locate all the corner and sharp corner points over the radar image and, based on the local densities of corner and sharp corner points, search for suitable objects for VBs.The idea is to select such areas of the radar image where there are many corner and sharp corner points, i.e. local maxima of their densities.The corners are searched to locate locally unique features containing nonlinear edges, because linear edges often are similar along the edge and can lead to similarization errors in the tracking process.
We select the points which locally (within a given radius R b ) maximize the complexity function F c : where N c (r, c, R s ) is the number of corner points within a search radius R s from the location described by the column and row coordinates (r, c), and N sc (r, c, R s ) is the number of sharp corner points within the same area.To avoid assigning VBs too close to each other, we only perform the search R s or more outside the already assigned VB locations.We have used values R b = 30 and R s = 15 pixels here, but these parameters can be defined by the user.
VBs are added always when the number of VBs becomes less than a given threshold T N .T N can be defined as an absolute value or relative to the number of original VBs.In the experiments presented here we have used relative T N values of 75-90 % of the number of the original number of the VBs.
www.the-cryosphere.net/10/29/2016/The Cryosphere, 10, 29-42, 2016 Because we apply the HMF prior to the VB selection, we also get relatively many VBs in the far range, because the filtering also amplifies the noise in the far range.For this reason, to properly initialize the VBs we perform one VB tracking iteration between the first and second images of the sequence to remove the VBs due to this noise amplification.This first iteration removes the VBs generated by random fluctuation (amplified radar noise).
In the areas of no distinguishable ice features (low signalto-noise ratio) in the imagery the tracking can not be reliably performed.For example in the far radar range this is typically the case.Typically there are enough ice features available and suitable for tracking ice drift in the drift ice field; within drift ice deformation continuously occurs, leading to ice features visible in radar (e.g.ice ridges).The algorithm only adds and tracks VBs in the areas where such traceable targets exist.

Optical flow and algorithm implementation
Optical flow (Horn and Schunck, 1981;Beauchemin and Barron, 1995) is a method used for estimating motion in image sequences such as in digital video.In optical flow we assume an intensity at a location (x, y) in a digital image at time to be moving such that Using the Taylor expansion and assuming small motion, where HOT stands for higher-than-first-order terms.Dividing by t we get the optical flow equation: where I x , I y , and I t indicate the partial derivatives of the image signal with respect to x, y, and t.The changes of I at (x, y) in x and y directions and change of I in time can be estimated from an image pair.To perform the estimation additional conditions are needed.One practical approach is the Lucas-Kanade method (Lucas and Kanade, 1981) where it is assumed that optical flow equation holds for a block of N pixels p k (k = 1, . .., N ): This corresponds to a (overdetermined The block of N pixels consists of pixels within a roundshaped window around a centre pixel (VB centre).This linear overdetermined system of equations can be solved (least squares solution) (Penrose, 1955) as M = (A T A) −1 A T is known as the Moore-Penrose pseudoinverse.
We have used the following discrete estimates for I x , I y , and I t (Horn and Schunck, 1981): where I 1 and I 2 are the first and second (in this temporal order) image of an image pair, and (r, c) refers to the row and column coordinates which are used here instead of x and y.
The image pixel values with fractional coordinates are computed using bilinear interpolation.For numerical stability it is essential that the estimates for I x , I y , and I t are computed at the same spatiotemporal location.
The optical flow method is best suitable for short motion corresponding to short time differences, e.g. for our coastal radar data with a relatively short time difference (in our case 2 min).In practice some image smoothing at sharp edges is recommendable before the optical flow computation, because optical flow assumes continuity of the signal.For this reason we perform a spatial Gaussian smoothing of the images before the optical flow estimation.The Gaussian smoothing is combined with the original image data by a linear combination to get the smoothed pixel value I (r, c) from the original pixel value I (r, c) and the smoothed pixel value G(r, c) (G refers to the image convolved with a predefined Gaussian kernel): We used a Gaussian kernel with standard deviation σ = 15.0 (pixels) in the Gaussian smoothing, i.e. convolving the original signal with the Gaussian kernel, and for the factor f = 0.8, i.e. the original image pixel value has a weight of 0.2 and the smoothed pixel value a weight of 0.8 in the linear combination used.
In computation of the optical flow we used a spherical window with a radius of 11 pixels (resulting to N = 377 pixels involved) for the coastal radar data and a a spherical window with radius of 21 pixels (N = 1373) for the ship radar (with longer time difference between the successive images) in the optical flow estimation.
If the cross-correlation between two matched windows is less than a given threshold T cc the object is not tracked any more, i.e. the VB is lost.We have used T cc = 0.9 in our experiments presented here.We also studied the use of the coefficient of determination of the linear fits as a measure of the matching quality, but it seems to have insignificant correlation e.g. with respect to the cross-correlation, which seems to be a more useful measure.The first tracking iteration (between the first and second images) is used to prune the unreliable VBs produced by noise amplification due to the HMF.The same cross-correlation thresholding is applied for this purpose.
6 Experimental results

Coastal radar image sequences
The initial radar images (first images in the image sequences) of the two coastal radar image sequences of cases A and B, with the locations of the initial VBs indicated, used in our experiments are shown in Fig. 6.The location of the radar indicated by the green dot and the initial VBs are indicated by red dots and the location of the radar by green dots.
The first and last images of case A and case B are shown in Figs.7 and 8 respectively.From these images we can see the change of the ice field during the whole study period.In case A a large ice field is torn off and drifting away from the land and landfast ice zone.In case B a smaller part of  the ice is also torn off and floating away from the coast.The trajectories resulting from the Lucas-Kanade optical flow algorithm for the two test cases are shown in Fig. 9a and b.We can see that for case A the direction of the motion was rather uniform over the whole drift ice area, but for case B the direction of the motion was less uniform.However, this kind of trajectory plots do not show the ice drift velocity evaluation as a function of time.
The VB trajectories computed by the algorithm correspond to the visual interpretation during the test periods.The visual inspection was performed using animations of the image sequences with the VBs indicated by coloured circles over the radar imagery.The information derived from VBs gives us possibilities to estimate different parameters related to the ice drift.We have computed some features for three selected VBs on both the images.Three VBs were selected to initially be close each other such that they form a triangle.Based on the triangles formed by the VB triplets we could compute the evolution of the area of the triangles as a function of time, indicating local divergence or convergence.The trajectories of the selected VB triplets are indicated in red in Fig. 9.For case A a large part of the ice was torn off the fast ice and drifting to north/northwest (Figs. 7 and 9a).For this case we applied a threshold T N = 75 % of the original number of VBs.The number of VBs as a function of time for this case is shown in Fig. 10a.We adjusted T N this high just to demonstrate the adding of VBs, in practice a lower value could be used.For case B some smaller ice floes were torn off and drifting to the northwest and west (Figs. 8 and 9b).For this case we used T N = 90 % of the original number of VBs.The number of VBs as a function of time is shown in Fig. 10b.Because the total time was longer for this case the VBs were added twice during the whole time period of 24 h.
It is also straightforward to compute the velocity and direction (here given as the compass direction with 0 • to the north, 90 • to east, and so on) by just dividing the displacement converted to metres by the time step of 2 min for each successive image pair.The velocities and direction for the selected three VBs of case A are shown in Fig. 11.It can be seen that the three VBs are moving rather coherently.The velocity is first very slow and the direction rather ambiguous.Then the velocity is accelerated rather rapidly (the average acceleration a can be estimated by dividing the velocity differences between the two adjacent velocities by the length of the time step: a ≈ 3.0 × 10 −5 m s −2 ).It can also be seen that as the direction changes for a short time, the velocity is reduced (around time 600 min).In Fig. 12 the velocity and direction are shown for case B. Also in this case the ice is quite stable for a while in the beginning of the time period and then the velocity accelerates (in the beginning, a ≈ 1.5 × 10 −5 m s −2 ; later the average acceleration in general is lower but exhibits temporal variation) to a top value of about 0.3 m s −1 achieved near the end of the period.For case A the top velocity was a little lower and achieved about in the middle of the time period.For case B the velocities of the three studied VBs were less uniform than for case B. This also results to the larger divergence for this case compared to 25 February 2011 case.The ratio of the area of a triangle formed by a selected triplet of points to the area of the triangle formed by the same triplet in the beginning of the study time period are shown in Fig. 13 for the cases A and B. These figures indicate the local ice convergence (decreasing value) or divergence (increasing value).
For case A we can see divergence after the motion starts (after around 100 min) to around 300 min, and then about 100 min of convergence and then divergence again.The divergence is increased around 600 min, and finally after about 700 min the VBs move uniformly for the rest of the time.The VB velocity accelerates from 0 to about 0.2 m s −1 during the period from approximately 100 to 300 min; also during this period the relative area increases, indicating divergence.During the period from approximately 300 to 400 min the VB velocity remains about similar, but there is some variability in the velocities and directions of the single VBs leading to convergence during this period.After this, during the time period from about 400 to 600 min, the VBs move quite uniformly while their velocity first increases and then decreases.After 600 min the velocity of all the VBs accelerates from about 0.05 m s −1 to about 0.2 m s −1 and during this period we also see differences in the velocity of the different VBs and due to these differences also the relative area increases (indicating divergence).
For case B, there is a period of divergence from about 200 min to about 800 min, and after that the VBs move quite uniformly to the end of the period.During the diverging period, there seem to be one VB (the northernmost one) moving faster than the two other VBs, mostly explaining the divergence.The relative area changes were from 1.0 to 1.2 for case A and from 1.0 to 1.5 for case B. In case B the ice concentration in later parts of the period was quite low, but in case A the VBs were part of a larger drifting ice field.
The ice drift estimated by VBs was also compared to the simple free drift model (Lepparanta, 2009).According to the model the ice drift velocity is the (surface) wind speed multiplied by a factor N a , also known as the Nansen number.A typical value of N a for the Baltic Sea is 0.025 (Lepparanta and Myrberg, 2009) and 0.017 for the Arctic (Lepparanta, 2009).The free drift model ice drift direction in the Northern Hemisphere is the wind direction rotated clockwise by 20-40 • (Lepparanta, 2009).For case A, the wind speed was in the range 6-7 m s −1 for the first 4.5 h of the study period and then increased to 8-9 m s −1 for the 4.5 h.The last 5 h of the period the wind was mainly in the range 7-8 m s −1 .According to the free drift model the ice drift velocities during these three periods would be about 0.16, 0.21, and 0.19 m s −1 .These values correspond to the VB velocities in Fig. 11a rather well.In the beginning of the period there was no significant drift until the ice was torn and started to drift.The wind direction was in the range 150-180 • and the ice was moving approximately to the north.This is also in ac-  cordance with the free drift model ice direction, taking into account that the Finnish coastline in the east was restricting the eastward drift component.
For case B the wind speed was approximately in the range 4-6 m s −1 during the first 8 h of the study period.During the second 8 h period the wind speed was in the range 2-4 m s −1 , and during the last 8 h period it increased to approximately to 4-5 m s −1 .The wind direction during the period varied from about 180 • to about 90 • and then back to about 150 • .According to the free drift model the ice drift velocities in the three 8 h periods would have been approximately 0.12, 0.08, and 0.11 m s −1 .There was, just like in case A, a static period in the beginning of the study period, but after the ice motion started the VB velocities were somewhat larger than the values based on the free drift model.One reason for this may be that the ice was relatively thin and drift ice concentration was low as the ice proceeded faster further from the coast.The ice drift direction of Fig. 12b corresponded to the free drift model quite well.

Ship radar image sequence
The RV Lance ship radar data were collected during the period from 21 January to 18 February, and the temporal difference between each image pair was 10 min, which is quite long for an optical flow algorithm.The ship was drifting in ice and no ship motion correction based on ship GPS has been performed; i.e. the detected ice drift is the drift with respect to the ship.The ship motion could be compensated based on the ship GPS, but we did not have this information available when this study was performed.
We computed the tracking for the whole period but here we only show the results for a 1-day period with some interesting motion.During many of the days the motion with respect to the ship was not significant.Because the ship radar had a shorter range and higher resolution there were sea ice details visible over the whole image.
For demonstration we selected the 8 February 2015 data for the whole day (case C) with significant ice motion with respect to the ship in some areas within the ship radar range.
Also in this case we selected three adjacent VBs, indicated in red in the trajectory image of Fig. 15, for which we computed their velocities, directions, and the divergence based on the area of the triangle formed by these three VBs.The first image of the 1-day time period (8 February 2015, 00:00 UTC) and the last image of the period (9 February 2015, 00:00 UTC) are shown in Fig. 14.It can be seen that two larger leads are opening, one in the northeastern part and another in the southern part of the images.This is also indicated by the VB trajectories shown in Fig. 15.The drift velocity, direction, and divergence for the selected VB triplet www.the-cryosphere.net/10/29/2016/The Cryosphere, 10, 29-42, 2016   are shown in Figs.16 and 17.In this case the direction changes quite suddenly around the 300-400 min mark.However, there is no significant reduction in the computed drift velocity.This may be due to the relatively sparse time step of 10 min; with 2 min time steps the observed speed change might have been larger.It can also be seen that some convergence occurs for the ice area described by the selected triplet of VBs.The area is decreased to about 80 % of the original area.This is because the VBs belong to the ice field east of the opening lead and the ice is converging in this   For case C, the wind direction also varied in the range 300-330 • , and wind speed was around 8 m s −1 .According to the free drift model this would result in drift direction of about 90-120 • and speed around 0.015 m s −1 .These values approximately correspond to the ship drift speed and direction.As the VB drift was estimated with respect to the ship we were not able to apply the free drift model to the estimated VB velocity and direction of Fig. 16.We can see that in the beginning as the ice drift with respect to the ship is slow the ice drift direction is approximately the same as the ship drift direction (ice is moving slightly faster than the ship to approximately same direction), but after the non-homogeneous dynamic ice drift events begin there are some relatively rapid changes in the VB velocities and directions.These relative changes were due to the local surface current and wind deviations.The compression in the triangle described by the triplet of VBs started around 1200 min and continued to the end of the period.During this period the velocity also decreased proportionally to the decrease of the triangle area.This is also an expected result as the ice velocity tends to reduce in compressing ice.

Evaluation of the estimation error
We performed a study in which we compared some visually selected points within some target windows (VB area) of a moving VB and manually defined the location after each 200 min steps (100 2 min time steps) in the four Tankar coastal radar VB sets.The selected points were not necessarily in the middle of the window (which corresponds the VB location in the algorithm), but such points which could be visually located in the time series of radar images.If we assume that the VB area remains unchanged (acts as a rigid object) then the difference between the VB centre and the studied point should remain constant.We computed the standard deviation of the difference for our selected points and the VB centres for the row and column coordinates and got values σ r = 1.41 and σ c = 1.24.The manually defined locations were always integers; i.e. the locations were defined visually from the image by pointing the pixel by a mouse and recording the pixel coordinates with a precision of one pixel.There are three sources of error in the manual estimation: human estimation error, rounding error to an integer pixel (on average 0.25 pixels in both directions), and changes due to possible internal transformations within the VB area.Based on this analysis we can only estimate the error in VB position to roughly be less than or equal to 1 pixel (33 m) in both coordinate directions, because it is difficult to estimate the exact contributions of the error sources present.
In Fig. 18 we show the locations of the VBs (red) and locations of the manually tracked points within the VB area for one case A time series.Also the corresponding VB and some ice area around it at the time instants 0, 200, 400, 600, and 800 min are shown in the figure.A thorough evaluation using this method would in practice be labour intensive.An estimation error of roughly 1 pixel in both coordinate directions in 2 min seems to overestimate the actual error; in addition, we also studied other methods of estimating the error.
Better accuracy for the estimation error can be reached by estimating the differences of both the coordinates when the ice is not moving based on visual interpretation.Such a situation is e.g. at the beginning of case B. We selected a 160 min time period from the beginning of case B time se- ries and computed the means and standard deviations of the successive (with a time difference of 2 min) location differences of r and c coordinates (row and column) for three VBs of this time series.We got for the values µ r = −0.020and µ c = −0.035for the means and σ r = 0.126, and σ c = 0.120 Tankar radar pixels (33 m) for the standard deviations.We also performed a similar study with the RV Lance case by selecting a VB which was visually static during the whole 1day study period and computed the values µ r = −0.00076and µ c = −0.000393for the means and σ r = 0.145, and σ c = 0.142 RV Lance radar pixels (12.5 m).The standard deviations can be considered as estimates of the location estimation error in the row and column directions.Possibly there still were some minor ice motion (not easily visible by eye) in the Tankar radar case as the values of µ r and µ c were significantly higher than for the RV Lance data.
For comparison, the error of GPS is 15m or less for 95 % of the time (Hofmann-Wellenhof et al., 2008).Some additional error in the GPS location is caused by the water level changes.This gives an idea of the location error of buoys equipped with GPS positioning.
The accuracy of the derived quantities such as velocity, direction, and area of a triangle formed by three VBs can be estimated based on the multivariate Taylor expansion; if we assume the errors are relatively small, we can also assume the higher-order terms to be neglectable and estimate the error by the first-order terms.If the derived quantity Z is a function of the original quantities X i , i = 1, . .., n, i.e.Z = F (X 1 , . .., X n ), then the error U at (X 1 , . .., X n ) can be estimated as, see e.g.(Steward, 1996), We can compute the error estimates for the velocity v and direction θ from their formulas v = s/t, and θ = arctan(dr/dc), where s is the drift of the VB between the 2 or 10 min time period, t is the time (2 or 10 min), dr is the row direction drift, and dc is the column direction drift within a given time period.Also the corresponding error estimate for the triangle area A can be estimated based on the equation used for the triangle area A: where (r 1 , c 1 ), (r 2 , c 2 ), and (r 3 , c 3 ) are the corner row and column coordinates of the triangle, i.e. coordinates of the three VBs forming the triangle.Making some simplifying assumptions, e.g. that the error in time measurement is neglectable, we get the estimates v ≈ 0.05 m s −1 for the Tankar cases (A and B) and v ≈ 0.004 m s −1 for the RV Lance case (C).The derived error estimation formula for the ice drift direction has the ratios of dr 2 +dc 2 in the denominator, and for this reason it will give unreliable estimates for small values of drift in the case both dr and dc are small.According to our experiments the estimates of the direction are not very reliable (include a lot of variation) for velocities lower than approximately 0.1 m s −1 .This can be seen in Figs.11b and 12b in the beginning of the two periods as the ice motion is small, or does not exist at all, as large variations of the direction estimates.According to our evaluation the relative error in the computed triangle (formed by VB triplets) size is about 0.5-0.6 % of the triangle area for case A, 0.8-1.1 % for case B, and 1.3-1.5 % for case C.

Discussion and conclusions
An algorithm for ice tracking was developed to track virtual ice buoys.The algorithm enables continuous ice drift tracking by adding VBs after a certain absolute or relative amount of VBs has been lost; this number or ratio can be defined by the algorithm user.The location of the VBs are initialized automatically based on the local information content of the image favouring locations with well-distinguishable features.The VBs can also be initialized manually if some specific case studies are made.This method is best suitable for relatively short temporal steps and giving sub-pixel resolution, unlike our earlier algorithm based on cross-correlation (Karvonen, 2013a).From the VB tracks many kinds of derived quantities, such as divergence, shear, vorticity, and total deformation, can be computed and information on the nature of the ice dynamics extracted.For example it is possible to distinguish between deformed ice and level ice or open water at least in the near range area, and opening of new open water channels can also be identified.
The main advantages of this algorithm are the capability of sub-pixel resolution and capability to monitor sea ice alone and continuously by adding VBs when needed.The estimation cannot be performed in the areas where no distinguishable ice features are present in the radar imagery.However, usually there are such features available in the radar images of drift ice and VB tracking is possible.Also the imaging geometry imposes some restrictions, especially in the far range.The typical ice ridges in the Baltic are at maximum several metres high and are not very well visible in the far range as their backscattering towards radar is less than that of larger (higher) targets.However, their shadowing effect is also less, even though for example a typical 2 m high ridge has a shadow of approximately 1 km in the range of 15 km, assuming the radar to be at a 30 m altitude from the sea surface and (locally) flat earth geometry.This may reduce the number of targets suitable for a VB in the far range.
The lost targets were typically lost as they had drifted away from the radar and signal-to-noise ratio had become lower.In general the number of traceable targets depends on the density of good scatterers (deformed ice, e.g.ridges) within the radar range.The algorithm radius parameters have been selected such that adjacent VBs are not very close to each other.If there are images of poor quality in the radar image sequence, e.g.due to weather conditions or radio-frequency signal interference from other sources (e.g.other radars), more VBs can be lost and the tracking can temporally be interrupted, but it will be restarted (by adding new VBs) as soon as the radar image quality has recovered.
We applied the developed algorithm to three test cases.Here we only computed some relatively simple derived quantities related to the ice dynamics, i.e. velocity, direction, and the relative area defined by triplets of VBs.The relative VB triplet areas give information on the divergence and convergence of the ice and also to some degree on the compression in the ice.Some other and more sophisticated ways to analyse VB data have been presented e.g. in Karvonen et al. (2013).The main purpose of this study was the development and testing of the VB tracking system.
The VB drift results for the test cases were evaluated based on visual inspection of the animations of the image sequences with overlaid VB positions.There were no real buoy data available.The visual analysis showed that the algorithm results correspond to the visual interpretation and same targets were tracked throughout the radar image sequences.This could also be visually verified by extracting single frames of the radar image time series with a time difference of e.g.2-3 h and with some given VB positions indicated on the radar images.Also, according to this verification the algorithm tracking results and visual inspection were in good agreement.Also some attempts to estimate the estimation error numerically were made for the test cases.
This VB tracking software complemented some basic VB data analysis software tools will be included in the radar servers, making automated ice dynamics analysis in real time or near-real time possible.The analysis tools will include the analysis of VB drift velocity and direction as well as divergence based on triangles formed by VB triplets.These will be computed within the convex polygon defined by the outer VBs in the image in a given grid and for a given time step (multiple of the basic time step).This integration is under construction as part of the project "Harnessing Coastal Radars for Environmental Monitoring Purposes" (HARD-CORE) funded by the Baltic Sea Research and Development Programme (BONUS).

Figure 1 .
Figure 1.Coastal radar and data transmission.The temporal median filtering is performed on the radar server; a PPI radar image is transmitted every 2 min by the GSM modem through the GSM base transceiver station (BTS), GSM base station controller (BSC), and mobile switching centre (MSC) and further through the Internet and to the FMI server performing the radar motion analysis.

Figure 2 .
Figure 2. Flow diagram of the continuous VB tracking.

Figure 3 .
Figure3.The principle of homomorphic filtering.A logarithmic transform is applied to make the illumination and reflection terms additive, a high-pass filter is then applied, and an exponential transform is applied to invert the log transform.

Figure 4 .
Figure 4.One temporal median filtered radar image (25 February 2011, left) and the same image after the homomorphic filtering.The land areas according to our rough land mask are masked off (white area in the figures).The associated scale is in kilometres.

Figure 5 .
Figure 5. Geometry of an LBP with an angular step of /4.The circle has a radius of R, which is one LBP parameter.
linear system Av = b, where v = [v x v y ] T and A =   I x (p 1 ) I y (p 1 ) . . .I x (p N ) I y (p N )

Figure 6 .
Figure 6.The initial VBs of the February 2011 test data set (a), and the February 2012 data set (b) drawn over the first radar images of the sequences.The coordinates are in kilometres, and the radar is located at (row, column) = (20, 30).The land area according to our rough land mask is indicated in green.The associated scale is in kilometres.

Figure 7 .
Figure 7.The first and last image of the February 2011 test data set.The coordinates are in kilometres, and the radar is located at (row, column) = (20, 30).The land area according to our rough land mask is indicated in green.The associated scale is in kilometres.

Figure 8 .
Figure 8.The first and last image of the February 2012 test data set.The coordinates are in kilometres, and the radar is located at (row, column) = (20, 30).The land area according to our rough land mask is indicated in green.The associated scale is in kilometres.

Figure 9 .
Figure 9. Trajectories of the VBs which survived the whole February 2011 and February 2012 test periods.The three VBs whose properties are studied in more detail in both the cases are indicated in red.The starting point of a trajectory is indicated by an open circle and the end point by a closed circle.The coordinates are in kilometres, and the radar is located at (row, column) = (20, 30).The land area according to our rough land mask is indicated in green.The associated scale is in kilometres.

Figure 10 .
Figure 10.Number of VBs at each time step for the 2011 case (a) and for the 2012 case (b).The threshold to add new VBs was 75 and 90 % of the original number of VBs respectively.

Figure 11 .
Figure 11.Velocity (a) and direction (b) for the 2011 case selected three VBs.The direction is in degrees.

Figure 12 .
Figure 12.Velocity (a) and direction (b) for the 2012 case selected three VBs.The direction is in degrees.

Figure 13 .
Figure 13.Ratio of the area of the triangle formed by the selected three VBs with respect to the area in the beginning of the period for the 2011 case (a) and for the 2012 case (b).

Figure 14 .
Figure 14.The first and last radar images of the RV Lance 1-day test period.The ship radar is located in the middle of the image.The associated scale is in kilometres.

Figure 15 .
Figure 15.Trajectories of the VBs which survived the RV Lance 1-day February 2015 test periods.The three VBs whose properties are studied in more detail in both cases are indicated in red.The starting point of a trajectory is indicated by an open circle and the end point by a closed circle.The associated scale is in kilometres.

Figure 16 .
Figure 16.Velocity (a) and direction (b) for the RV Lance 2015 case selected three VBs.The direction is in degrees.

Figure 17 .
Figure 17.Ratio of the area of the triangle formed by the selected three VBs with respect to the area in the beginning of the period for the 2015 RV Lance case.

Figure 18 .
Figure 18.An example of manual VB pixel location estimation and VB location defined by the algorithm with a time step of 200 min (case A).The locations are shown in (a): the red dots indicate the VB centre location, the blue dots the selected and tracked feature location, and the corresponding pairs have been marked with the green rectangles.The VB areas (red circles) in the radar imagery with some background included at the given time instants are shown in time order from left to right in (b)-(f).