- Email: [email protected]
A methodology based on Particle image velocimetry for river ice velocity measurement
Cold Regions Science and Technology 89 (2013) 36–47
Contents lists available at SciVerse ScienceDirect
Cold Regions Science and Technology journal homepage: www.elsevier.com/locate/coldregions
A methodology based on Particle image velocimetry for river ice velocity measurement Anik Daigle a,⁎, Francis Bérubé b, Normand Bergeron b, Pascal Matte b a b
Cégep Garneau, 1600 de l'Entente, Québec (Québec), Canada, G1S 4S3 Institut national de la recherche scientifique, 490, rue de la Couronne, Québec (Québec), Canada G1K 9A9
a r t i c l e
i n f o
Article history: Received 25 October 2012 Accepted 18 January 2013 Keywords: River ice Particle image velocimetry Remote sensing
a b s t r a c t The in-field use of Particle image velocimetry (PIV), or Large-scale PIV, may be complicated by implementation issues related to the movement of the camera, the variations in luminosity and the distortion of the images due to the inclination of the shore-based camera. The correction for the distortion of the images, or orthorectification, is a critical issue that usually requires the georeferencing of control points located on the river surface. In drifting ice conditions, control points cannot be identified directly on the river surface, and water levels are often unstable. A robust, flexible and relatively inexpensive methodology based on Large-scale PIV for the measurement of river ice velocity is presented. Two orthorectification techniques are proposed, that do not require georeferencing of control points directly on the river surface. The proposed methodology also includes image enhancement, correction for camera movements and a segmentation technique that allows dynamical quantitative estimations of the ice cover extent. The methodology is tested on the St. Lawrence and Montmorency rivers, two rivers of very different sizes and ice conditions. The correction for camera movements proves to be quite reliable and robust to illumination variations and to the presence of shade patterns. The two orthorectification techniques allow to estimate pixel physical dimensions with 6% and b 1% precisions, respectively. Velocities obtained for the St. Lawrence case are compared to velocities simulated by a numerical model. There is strong agreement in velocity values in moderate ice coverage conditions. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Understanding the processes of ice cover formation and ice transport is essential to the prediction of ice jams and flooding, and for the continuation of maritime transport in the winter. In Canada, nearly 44 million tons of goods travel through the St. Lawrence waterway each year (Shaw and Kaczkowski, 2009). Since 1996, the Canadian Coast Guard monitors the ice cover of the St. Lawrence River at various strategic locations using cameras, sonars and radars. This observation protocol allows experts to produce periodic ice condition reports. This system, however, works through considerable human and material resources that small municipalities struggling with seasonal ice jams usually cannot afford. In most cases, the analysis of ice movement is based on a qualitative visual assessment. The quantitative monitoring of river ice movement is difficult and usually relies on acoustic Doppler profilers providing velocity values at one point on the river (Chave et al., 2004). These are 0-D measurements, while the understanding of the formation of ice cover,
⁎ Corresponding author. Tel.: +1 418 654 2596; fax: +1 418 654 2600. E-mail addresses: [email protected] (A. Daigle), [email protected] (F. Bérubé), [email protected] (N. Bergeron), [email protected] (P. Matte). 0165-232X/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.coldregions.2013.01.006
transport, as well as the estimation of the amount of flowing ice, requires knowledge of the whole velocity field. Moreover, these are expensive immersed instruments that can be damaged by the presence of ice, and there may be important installation and operation safety concerns in flood conditions and near structures like bridges or hydroelectric plants. These limitations call for the establishment of a remote and twodimensional measurement method, two conditions that are met by Particle image velocimetry (PIV). PIV is an imaging technique that was adapted in the 90s for the measurement of the two-dimensional flow velocity field in laboratory experiments (Adrian, 2005; Bradley et al., 2002; Ettema et al., 1997). The stream is firstly seeded with small floating particles (such as plastic beads or rice grains) that act as surface flow tracers. The surface flow is then filmed, with a video camera or by taking successive images with a still camera. The velocity field is obtained by identifying subareas of the image, also called Interrogation Areas (IA, after Muste et al. (2009)), in two successive frames: the correlation coefficients between a given IA of the image acquired at time t and some neighboring area in the image acquired at time t + Δt are calculated. The new t + Δt coordinates of the IA are selected as those corresponding to the most highly correlated t + Δt subarea, and the velocity vector is obtained by dividing the displacement of the IA by Δt. This procedure, repeated on IAs covering the entire image and for all time steps, results in a dynamical and two-dimensional velocity vector field.
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
PIV thus allows the measurement of the spatial distribution of the surface velocity of a seeded flow in time, information that is crucial in a variety of practical applications such as the dispersion of pollutants, bank erosion, sedimentation, the local effects of flooding, and the interaction between the flow and riverine structures. Knowledge of the spatial distribution of water and/or ice velocity also gives a fundamental insight into the hydrodynamics of rivers, providing reliable data for calibration and validation of numerical models. The in-field use of PIV (also known as Large-scale PIV or LSPIV) was initiated by Fujita et al. (1998), who presented a series of adjustments required to move from laboratory and small-scale applications to in-field and larger scale uses. Among other parameters to be accounted for, the authors mention the variable outdoor luminosity and distortion of the images due to the inclination of the shore-based camera. LSPIV has since been used for the measurement of water surface velocities in various laboratory and field conditions (Bradley et al., 2002; Kantoush et al., 2011; Muste et al., 2009), including during floods (Dramais et al., 2011). The potential of LSPIV for field measurement of river ice velocity has however practically not been exploited. A method for the spatial description of ice velocity from video images was presented by Ferrick et al. (1992), but this method involved manual monitoring of individual pieces of floating ice (particle tracking velocimetry, or PTV). The first study explicitly involving LSPIV for the measurement of ice velocity was presented by Ettema et al. (1997). LSPIV was used in a laboratory model of the confluence of the Mississippi and Missouri rivers, to estimate the velocity field and the flow of ice, simulated by polypropylene particles. The first field LSPIV experiment applied to river ice was presented by Jasek et al. (2001), who estimated the flow of the Yukon River during an ice jam near Dawson City, using the ice floes as tracers of the flow. Tardif (2001) presented a method to extract ice velocity from images using a method similar to LSPIV. The author however pointed out several implementation issues, including the reduced visibility during snowfalls, the low contrast in the images at sunrise and sunset, and the presence of reflections on the water surface complicating the recognition of IAs from one image to another. Similar problems, as well as unavoidable movements of the camera, were reported by Bourgault (2008), who attempted to use LSPIV to obtain ice velocity fields on the St. Lawrence River in Québec City. The evaluation of ice flow depends on three measurements: the drifting ice velocity field, the ice cover extent, and the ice thickness. This paper presents a robust, flexible and relatively inexpensive methodology based on LSPIV and image segmentation for the two first measurements. Results obtained for two case studies are presented. The LSPIV methodology is presented with particular focus on three major issues pointed out by preceding authors for the in-field implementation of LSPIV for the measurement of ice velocity: the correction for camera movements; the enhancement of images in order to account for low contrast, variable luminosity and reflections; and correction for the image distortion due to the inclination of the camera, also called orthorectification. This latter critical issue usually requires land surveying equipment and expertise, and involves the georeferencing of control points on the surface to be orthorectified. In drifting ice conditions, control points cannot be identified directly on the river surface, and water levels are often unstable. Two alternative orthorectification methods are presented: one that uses control points that do not need to be defined directly on the river surface and thus may be georeferenced a posteriori, and another that does not require georeferencing at all and is particularly suited for small rivers. A segmentation technique that allows dynamical quantitative estimations of the ice cover extent is also presented. The technique is based on the classification of the image pixels in the two “ice” and “not ice” categories, based on image texture and color information.
37
2. Study sites and data acquisition Two study sites were selected for the development and testing of the methodology. They were chosen to be different in sizes and ice conditions, in order to test the flexibility of the proposed method. The St. Lawrence River near Québec City is about 1-km wide, but narrows to ~ 600 m between the Québec and Pierre-Laporte bridges. Tides have a maximum amplitude of ~ 5 m, and water velocity can reach 3 m s −1 in both directions. Ice forms on the river every winter. The waterway and maritime transport are maintained by the Canadian Coast Guard throughout the winter, but a steady stream of ice floes generally occurs from December to March. Images were acquired from the South-shore, at an elevation of ~35 m above the water level (Fig. 1a). This site was chosen because of its accessibility, a relatively high elevation above the water level allowing a more perpendicular line of sight, and the presence of permanent structures in the field of view (e.g. the bridge pillars) that facilitate image georeferencing.
Fig. 1. St. Lawrence River, between the Québec and Pierre-Laporte bridges. (a) The three control regions used for the correction of camera movements, as found in a reference image (dashed boxes). (b) Locations of control region 1 in the image to be corrected (solid line box), and in the reference image (dashed box). (c) Control region 1 in the transformed (translated and rotated) image.
38
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
The Montmorency River flows into the St. Lawrence River about 15 km north-east of Québec City. Its drainage basin is about 1100 km 2. An ice cover forms from December yearly, and the river has experienced several ice jams in its history (Leclerc et al., 2001). The images were acquired from a bridge crossing the river about 10 km from its mouth, where the river is about 50-m wide. All images were acquired with a Canon Rebel T2i DSLR camera with 18 Mpix resolution mounted on a regular tripod. In burst mode and using a 30 Mb s −1 memory card, sequences were recorded at a rate of 3.4 images s −1. 3. Methods All development was conducted using Matlab R2011b. 3.1. Image stabilization Small camera movements due to wind and/or oscillations of the structure on which the camera is based can rarely be completely avoided in the field. The effect of such movements is a small continuous variation of the field of view of the camera creating shifts in the coordinates of the IAs to be tracked by the PIV algorithm. These shifts are not necessarily the same throughout the image, since the camera may experience small rotations in addition to translations. To ensure that the shift of a given IA is solely due to a real displacement of the water surface, every individual image of the sequence must be framed and scaled identically. To achieve this, one frame of the sequence is selected as the reference image and all other images are translated, rotated, and reframed according to this reference. This image transformation procedure, also called image registration, is performed for each image by the multiplication of a transformation matrix. A transformation matrix is a 3 × 3 matrix that specifies the required translation, rotation, shear and/or scaling operations for an image to be aligned to the reference image. Such a matrix can be inferred from a minimum of two control points pairs identified in the reference image and in the image to be transformed. Control points have to be fixed recognizable locations, such as buildings or reference poles. Since the measurement of ice velocity may often be needed in natural areas where no buildings are found, the proposed method does not use control points but control regions in the images that may be parts of the natural landscape such as tree trunks, shore ground or rocks. These regions must be as still as possible (unaffected by the flow or the wind) and located at opposite sides of the image in order to detect a possible rotation of the image (Fig. 1). Once control regions are identified on the reference image by the user (Fig. 1a), they are detected in all images of the sequence by an algorithm based on the correlation of the regions with some neighboring area — very much as in the PIV algorithm. The size of the control regions is not critical, but a larger region is more likely to be uniquely positioned in each image. A rather textured region also helps getting a precise and unique location in every image. The control points are then defined as the center of the detected regions. One transformation matrix is defined for each individual image, and each image is transformed accordingly. The result is a sequence for which the control regions maintain the same coordinates, with the expense of small peripheral losses (Fig. 1c). 3.2. Image enhancement Other unavoidable issues in outdoor implementation of PIV are the variation in luminosity and contrast due to clouding, reflections, and sometimes shade, on the water surface. There are several ways to enhance the quality of an image, often involving image filtering. The tool developed in this study makes use of an adaptive low-pass (wiener) filter for the subtraction of noise, and of an adjustment of
the intensity distribution of the image pixels for the improvement of contrast. The variation in luminosity due to the passage of clouds is dampened by the subtraction, for every image, of the mean pixel intensity value. The major effects on the images of this study were however due to reflections and shades on the water surface (e.g. Fig. 1a). An IA travelling through the limits of a reflection or a shade pattern is more difficult to track with precision. Reflections can be substantially reduced by the use of a polarizing filter. Remaining reflections and shades can be further reduced by the subtraction of a local background containing the shade and reflection patterns. For details about the image enhancement techniques that were used, see the wiener2 function (adaptive low-pass filtering of noise), the imadjust function (for the adjustment of the image intensity distribution) and the imtophat function (for the calculation of a local background) of the Matlab Image Processing Toolbox (The Mathworks Inc., 2012). 3.3. Image orthorectification PIV in itself provides velocity in pixels per second. To convert to physical space velocity, the true physical dimensions of the image pixels must be known. In oblique images, such as images taken from the shore, the pixel physical dimensions vary throughout the image. Correction for the oblicity of an image is called orthorectification. Typically, orthorectification is done by georeferencing control points in the image that define a transformation matrix. This matrix assigns new coordinates to every pixel of the image so that they correspond to the physical space. This procedure is relatively straightforward when the user has access to land surveying equipment and expertise, and when fixed georeferenced points can be found directly on the surface to be orthorectified. In the present application however, the surface is liquid, flowing and may even be moving up and down because of rapid water level changes due to local ice accumulations or tides. We thus developed two different orthorectification techniques adapted to the in-field implementation of PIV. The first one is based on an above-surface georeferencing, and makes use of available georeferenced points that do not need to be on the water surface. The other does not use georeferencing at all, and rather relies on the interpolation of several scale measurements provided by an object of known dimensions that is introduced into the flow during the image acquisition. 3.3.1. Above-surface georeferencing Georeferencing can often be facilitated by the presence of permanent fixed structures. Images of the St. Lawrence River section located between the Québec and Pierre-Laporte bridges, in Québec City, include bridges pillars and streetlights (Fig. 2). These are fixed points that are easily georeferenced with regular surveying equipment. Orthorectification based on above-surface georeferecing can be performed following these steps: 1. Georeference fixed points above the water surface level, and one at the water surface level. In the example shown in Fig. 2, 15 points were georeferenced (crosses). Point 15 was georeferenced at the water surface level, on one of the Québec bridge pillar. 2. Determine height scales at many locations in the field of view. In this study, height scales were determined using available information about the dimensions of the bridge pillars bricks and heights of the streetlight poles. 3. Relocate the georeferenced points on the river surface, using the local scale heights and the position of the water surface level. In Fig. 2, circled crosses mark the positions of the 15 georeferenced fixed points, when relocated on the river surface.
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
39
Fig. 2. Fixed points (crosses) georeferenced above the St. Lawrence surface and relocated (circled crosses) at the water surface level.
4. Use the relocated points as control points for the definition of the transformation matrix, and proceed to the image orthorectification. 3.3.2. Interpolation of local dimension references In most situations, especially in remote areas, permanent fixed points do not exist. One possibility is to install landmarks such as metal posts near the water surface and to proceed with the regular georeferencing and orthorectification procedure. In many situations, however, the installation of landmarks can be time consuming and difficult: for the transformation matrix to be reliable, the markers must be positioned near the limits of the field of view, and thus on each sides of the river. In flood conditions and/or with the presence of ice, one side of the river may be difficult to access and the shores can be unstable. In such urgent situations as the monitoring of a possible ice jam, these issues related to orthorectification may thus prevent the use of PIV. Another way of inferring pixel dimensions is to introduce a floating object of known dimensions within the field of view of the camera during image acquisition (Fig. 3). In the example presented here, the object is a floating square plate. The method can be summarized by these steps: 1. The plate is photographed or filmed at many locations in the field of view (Fig. 3); 2. Coordinates of the plate corners are identified for each location, defining squares of known dimensions and a local orthogonal reference system, and providing physical dimensions to the pixels of the squares (Fig. 4); 3. The local horizontal and vertical pixel dimensions obtained at each plate location are interpolated on the river surface. 3.4. Estimation of the ice cover extent The estimation of the ice cover extent is necessary for the quantitative evaluation of the amount of flowing ice, a variable crucial for the prediction of ice jams (Tardif, 2001). The PIV methodology developed here includes an evaluation of the ice cover for every image of a sequence, by the automatic recognition of ice pixels in the image. In image processing, the separation of an image into different parts is called segmentation. There are many different ways to perform segmentation, such as edge detection, texture and color classification
(Jain, 1989). The choice of the segmentation technique depends on the problem that is addressed. In a single river ice image, the ice can present very different colors and textures. In Fig. 5a, it is apparent that the texture of large snow covered ice floes in the foreground differs from that of the smaller ice floes. This difference is further increased by the effect of perspective. A segmentation based on color is not more reliable, with the shade created by clouds or elevated objects (here the Pierre-Laporte Bridge) causing abrupt color changes, and reflections making parts of the free water surface appear almost as pale as ice floes. For these reasons, segmentation based solely on image texture or pixel colors cannot provide reliable estimates of the percentage of ice cover. The combination of both, however, resulted in a clear separation of free water/ice pixels and proved to be especially robust when tested on different images (Fig. 5b). The procedure adopted uses an artificial neural network (ANN) to classify each pixel of the river surface in the two “ice” and “not ice” classes, based on texture and color information assigned to each individual pixel. It can be summarized as follows: 1. The image to be segmented is orthorectified; 2. The texture information is given by a measure of the image local entropy: for each pixel of the image, the entropy is computed on a 9 × 9 pixel neighborhood; 3. The color information is represented by the three-element RGB color array of each pixel; 4. The user identifies 25 “ice” pixels and 25 “not ice” pixels by clicking on the image; 5. These 50 pixels are used to calibrate an ANN that learns to separate “ice” and “not ice” pixels; 6. The calibrated ANN classifies all pixels on the river surface; 7. The ice cover percentage is computed as the sum of “ice” pixels over the total number of pixels on the river surface. The ANN is a three-layer perceptron that classifies each pixel of the river surface based on the four-element [entropy-RGB] inputs assigned to each individual pixel. The user identifies the 50 calibration pixels only once for a given sequence, and the segmentation is performed automatically for every frame, resulting in a dynamical estimation of the ice cover. Illumination variations are fairly well supported by the segmentation procedure, as an average background
40
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
Fig. 3. Montmorency River. A square plate of known dimensions is attached to a fishing line and photographed at several locations in order to provide the scaling information.
is subtracted to each image previous to the ANN classification. Large changes in the illumination or weather conditions however require a new calibration for the segmentation to remain reliable. 4. Results 4.1. Image stabilization and enhancement Image stabilization and enhancement were performed on the Montmorency and St. Lawrence images using the same procedure. Conditions near the St. Lawrence River are most of the time windy, and sequences acquired at the St. Lawrence site were more affected by camera movements than sequences acquired at the Montmorency site. An example of the effect of the stabilization routine is shown in Fig. 6. It shows the correlation coefficients between a control region as found in the reference image and in the unprocessed images (gray line), as well as in the stabilized images (black line). In this example sequence, 1530 images were acquired at 3.4 Hz between 15:34 and 15:41
on February 19th, 2012. Three control regions were defined as shown in Fig. 1a, and the reference image was chosen half-way through the sequence. Correlation coefficients for control region 1 are shown in Fig. 6a. The smallest coefficients between the reference image and the unprocessed images (gray line) correspond to the largest camera shifts. Correlation coefficients between control region 1 as found in the reference image and in the stabilized images are all >0.99. Correlation coefficients for control region 2 are shown in Fig. 6b. Here again, the coefficients increase drastically after the stabilization, but decrease progressively on both sides of the reference image time. This is due to the fact that no background subtraction was performed for this example, and that control region 2 is defined in an area that is affected by the moving of the shade of the Pierre-Laporte Bridge (Fig. 1a). Because control region 2 is well textured, it could still be reliably detected by the stabilization routine throughout the sequence despite this variable shade pattern. Control regions should nevertheless be chosen so as to minimize such alterations, or the performance
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
41
Fig. 4. Local orthogonal system allowing the inference of the local x- and y-pixel dimensions on the Montmorency River surface.
of the stabilization routine could be affected despite image enhancement efforts. The stabilization of one 18Mpix image takes ~ 33 s using Matlab running on an Intel Core Duo CPU [email protected] GHz with a RAM of 1.94 Go. The increase in the precision of the IA displacement measurement due to each of the image enhancement steps was measured by the mean increase in the signal-to-noise ratios of the IA correlation matrices. Our well-contrasted and low-noise images do not contain enough noise for the wiener filter, or the improvement of the contrast, to have a noticeable impact on the SNR ratios. The background subtraction, on the other hand, allows a substantial reduction of the contribution of reflections and shades to the images background and doubles the average SNR of the IA correlation matrices.
The six locations of the floating plate that were used to infer pixel sizes on the river surface are shown in Fig. 7. Two horizontal and two vertical reference scales are defined for each of the plate locations and interpolated on the river surface. Typically, on images acquired from the shore, pixels dimensions vary
4.2. Interpolation of local dimension references: the case of the Montmorency River The case of the Montmorency River was considered mostly to test the orthorectification technique based on local dimension references.
Fig. 5. Segmentation of the St. Lawrence surface between the Québec and Pierre-Laporte Bridges, based on the combination of texture and pixel color information. Image was taken on March 17th, 2011, at 14:42.
Fig. 6. (a) Correlation coefficients between the control region 1 as found in the reference image of the St. Lawrence River and in the unprocessed images (gray line) and the stabilized images (black line), as a function of time. Position of the reference image is indicated. (b) Correlation coefficients for control region 2, as a function of time. Coefficients for control region 3 are not shown. No image enhancement has been applied for this example.
42
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
Fig. 7. Locations of the floating plate that were used to infer pixel sizes on the Montmorency River surface. The reference lengths used to compute the local x- and y-pixel scales are indicated.
mostly in the image foreground-background direction (direction y in Fig. 7). Pixels horizontal and vertical dimensions were interpolated in the x and y directions by distinct linear and quadratic functions (Fig. 8). Performances of both interpolation types are evaluated by measuring the four edges of the plates at each of the six locations. The 24 lengths obtained are compared to the actual 61 cm length of the plate edges. Performance is evaluated by calculating the root mean square error (RMSE): rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 XN 2 RMSE ¼ ðL −LÞ i¼1 i N
ð1Þ
where L is the actual edge length (=61 cm), Li are the lengths as measured on the image based on the pixel size interpolation and N is the number of measurements (24). The RMSE achieved using a linear interpolation is 5.1 cm, while the quadratic interpolation achieves a RMSE = 3.7 cm. Velocities in cm s−1 can now be obtained on the water surface. Fig. 9 shows the average velocity vector field obtained on a 5-minute sequence acquired on December 10th, 2011. The instantaneous velocity components averaged on the subarea delimited in Fig. 9 are plotted as a function of time in Fig. 10. Ice cover of the whole river surface as a
function of time is also shown in Fig. 10. The right hand side column of Fig. 10 shows the distribution of the instantaneous velocity components values and estimated ice cover values during the 5-minute sequence. Velocity measured on the 5-minute sequence has a mean modulus of 30 cm s−1, with a 1.8 cm s −1 standard deviation. 4.3. Above-surface georeferencing: the case of the St. Lawrence River The orthorectified version of Fig. 2 is shown in Fig. 11, where the 15 control points are plotted at their orthorectified coordinates. Points 1–2 and 12–13–14–15 now appear superposed, since these share approximately the same coordinates on the river surface plane (see Fig. 2). The orthorectified coordinates of each point were compared to their georeferenced coordinates in order to evaluate the error associated to the orthorectification. The error on the orthorectified coordinates is evaluated by the RMSE and rRMSE, or relative RMSE, which is given by: rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 XN C i −C i i¼1 rRMSE ¼ N C max −C min
ð2Þ
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
43
whole river surface. Velocity measured for this 7-minute sequence has a mean modulus of 115 cm s−1, with a 5 cm s−1 standard deviation. 4.4. Uncertainty evaluation
Fig. 8. Quadratic interpolation of (a) the x- and (b) y-pixel dimensions of the Montmorency River image.
where Ci and C i are the ith georeferenced and estimated coordinates, and Cmin and Cmax are the minimum and maximum coordinates of the 15 points considered. This relative performance measure allows the rescaling of the coordinate values and a more representative description of the relative estimation precision. In the example shown in Fig. 11, the RMSEs obtained for the 15 transversal and longitudinal coordinates are 5.2 m and 0.63 m, respectively, which represent rRMSEs of 0.7% and 0.2%. The average velocity vector field obtained from a 7-minute sequence acquired on February 19th, 2012 is shown in Fig. 12. The instantaneous velocity components averaged on the subarea delimited in Fig. 12 are plotted as a function of time in Fig. 13, along with the ice cover of the
Velocities measured on the 5- and 7-minute sequences acquired on the Montmorency and on the St. Lawrence (Figs. 10 and 13) show respective coefficient of variation of 6% and 4%. It is difficult to evaluate which part of this variability is real, and which is attributable to measurement error. Several sources of uncertainty are associated to LSPIV that have been discussed in previous studies (Dramais et al., 2011; Hauet et al., 2008). The main sources of uncertainty identified for the particular applications presented here are (1) the orthorectification, and (2) the measurement of the IA displacement. Errors associated to both orthorectification techniques are systematic and thus affect the mean velocity values and not their variability in time. The evaluation of the IA displacement, on the other hand, is made for every time step. As stated earlier, this measurement is based on the computation of a correlation matrix between subareas of an image taken at time t and another taken at time t + Δt. This matrix is typically highly peaked at the position of the IA in the t + Δt frame. To achieve a sub-pixel evaluation of this position, the autocorrelation matrix is modelled by a two-dimensional Gaussian profile. The precision of the displacement measure is thus dependant on the quality of this fit, which in turn depends on the IA size, time step, and on the image contrast and SNR. The measurement of the IA displacement also depends on the image stabilization, also based on the computation and modeling of autocorrelation matrices. The image transformation, which involves an interpolation, also adds to the IA displacement uncertainty. In our 18 Mpix resolution well-contrasted and low-noise images, the uncertainty associated to the measurement of the IA displacement is evaluated to be ±1 pixel. The measured IA displacements being of the order of 10 pixels for both the Montmorency and St. Lawrence cases, the associated relative error can thus reach 10%, which is higher than the observed coefficients of variations. 4.5. Comparison of the St. Lawrence River velocities to H2D2 simulations Direct measurements of ice velocity are difficult and no direct validation of the obtained velocity values could be performed. The velocity values obtained using the PIV methodology on image
Fig. 9. Average velocity vector field obtained on a 5-minute image sequence of the Montmorency River acquired on December 10th, 2011.
44
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
Fig. 10. x- and y-velocity components averaged on the subarea delimited on Fig. 9, and ice coverage estimated on the whole river surface, plotted as a function of time (on the left), and as histogram distributions (on the right).
sequences of the St. Lawrence were however compared to velocities simulated by a numerical model. A 2D high-resolution numerical model of the St. Lawrence fluvial estuary (Matte et al., 2009) was used to simulate the hydrodynamic conditions prevailing during the analysis periods. The finite-element grid underlying the model has an average spatial resolution of 50 m, with refinements down to 1 m around engineering structures such as bridge pillars. The 2D vertically-integrated shallow-water equations are solved using the H2D2 software (Secretan, 2010). The model was calibrated under average flow conditions (~12,000 m 3 s −1) for a neap-spring tide cycle with an accuracy of b5 cm in water levels at most stations. However, simulated water velocities were not verified against measured data, so that deviations between modeled and measured lateral velocity distributions may appear, mainly where cross-sectional variations in bed friction are significant. The accuracy of results for other flow conditions was not verified either. For the purpose of this study, the domain consists of a 30-km reach of the river, with boundaries positioned at the upstream and downstream tide gauges closest to the Quebec Bridge, respectively located at Neuville and Quebec. The model is forced at the boundaries with water levels
measured at the tide gauges. Main tributaries, also influenced by tides, are not included in the model due to a lack of data at their outlets; this may affect the modeled velocities only within the radius of influence of the tributaries as their discharges are quite small compared to the St. Lawrence. Simulations were performed without ice since no information on the extent, thickness and roughness of the ice cover was available for the system at the simulation times. Velocities simulated by the numerical model and measured by the PIV methodology on the same transect on February 19th, 2012, show strong agreement, as shown in Fig. 14a. The same comparison was performed on the same transect of the St. Lawrence but for a sequence acquired on March 17th, 2011. Results are shown in Fig. 14b. This time, velocities simulated by H2D2 are more than twice as high as the velocities measured by the PIV methodology.
5. Discussion The main processing difference in the two case studies presented is in the orthorectification techniques that are proposed.
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
Fig. 11. Orthorectified image of the St. Lawrence River between the Quebec and Pierre-Laporte bridges.
The orthorectification based on the interpolation of local dimension references is proposed as an alternative to the traditional georeferencing of control points located on the surface to be orthorectified. This alternative technique is based on the introduction
45
of a floating object of known dimension in the flow, and is thus particularly suited for situations where land surveying material is not available, and/or when deployment time is a critical issue. This technique does not require georeferencing, is particularly affordable, can handle rapid water level variations and takes almost no time to deploy. It can be handled from one shore only if the river width so permits. In the Montmorency River case that is presented, best results are obtained with a quadratic interpolation of the local dimension references. This small nonlinear trend, more pronounced in the horizontal direction, is probably due to the fact that no correction was made to account for the camera lens distortion. This technique is more suited to small to average sized rivers: large rivers involve large distances, over which the position of the floating object is more difficult to control. The use of a floating object on large rivers would also require the object to be large enough for it to remain visible on the images even when located near the shore opposite to the camera. A floating object of proper dimension attached to a boat could allow the use of this technique on large rivers. Above-surface georeferencing differs from traditional georeferencing by the fact that it is based on available fixed points that are not necessarily positioned on the river surface. Georeferencing may thus be performed after the image acquisition, a significant advantage in urgent situations, or when the water level is changing. The methodology corrects for the height of the control points given known height references in the image, and then proceeds to a regular image registration. The validation of the obtained velocity values remains a challenging problem. Direct measurements (for example, ADCP measurements) could provide reliable validation values. H2D2 simulations of the St. Lawrence velocities nevertheless provide an insightful comparison basis. Velocities obtained by both the PIV-based method and H2D2 agree very closely for the 2012 sequence, but show an important but quasi-constant difference for the 2011 sequence. An interesting hypothesis is that this discrepancy is due to the very extensive ice cover present on the river on March 17th, 2011: ice cover values varied between 75% and 80% with large ice floes in contact (Fig. 5), while the ice cover on the February 2012 sequence did not exceed 50% (Fig. 11). Surface velocity on March 17th, 2011 was thus probably substantially lower than the water-column averaged velocity. The fact that the numerical model does not account for the presence of
Fig. 12. Average velocity vector field obtained on a 7-minute image sequence of the St. Lawrence River acquired on February 19th, 2012.
46
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
Fig. 13. x- and y-velocity components averaged on the subarea delimited on Fig. 12, and ice coverage estimated on the whole river surface, plotted as a function of time (on the left), and as histogram distributions (on the right).
ice also certainly explains part of the observed discrepancy. More comparisons on more sequences showing a range of ice cover values would be needed to verify this hypothesis. The proposed method includes an estimation of the ice cover as a function of time based on image segmentation. The ice cover values obtained on the Montmorency sequence and presented in Fig. 10 show a number of sudden jumps to ~ 75%. This is due to the important presence of steam near the water surface: the water vapor takes a coloration and texture similar to the ice under particular illumination conditions, and this gets to trick the segmentation routine for a few frames at a time. These singular ice cover values can be discarded visually based on the histogram distribution, or automatically filtered. Other specific conditions, such as a snowfall, would probably have similar consequences. Additional tests in more challenging illumination conditions (e.g. during a snowfall or in foggy weather) would allow to better evaluate the effects of image enhancement on the reliability of the obtained velocity and ice cover values. A river discharge can be obtained by multiplying the area of a given cross section and the average velocity of the water through this section. By providing the surface velocity information, PIV could allow a real-time evaluation of river discharge, given the bathymetry
information for the observed transect and an independent measurement of the water level elevation. The segmentation procedure presented here is a promising first step towards a dynamical evaluation of the amount of flowing ice using a PIV-based methodology. The conversion from ice cover to ice volume requires an estimation of the ice thickness, that cannot be obtained from optical images such as the ones presented here. Studies however showed that the ice thickness can be measured by an upward looking sonar instrument (Jasek and Marko, 2008; Morse et al., 2003) or be inferred from thermal images (Emond et al., 2011). 6. Conclusion A methodology based on PIV for the measurement of river ice velocity was presented and tested for two rivers of very different sizes. It was developed to be a simple and affordable implementation procedure allowing organizations and users with limited resources to monitor drifting ice velocity and ice cover. Frequently mentioned implementation problems such as camera movements, variable luminosity, reflections, shades and orthorectification are taken into account. The methodology allows for obtaining velocity fields and ice cover estimations as a function
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
47
References
Fig. 14. Velocities obtained on the same transect of the St. Lawrence River by the PIV methodology (gray line) and the H2D2 numerical simulation (black line) on (a) February 19th, 2012 (~14 h30) and (b) on March 17th, 2011 (~14 h45).
of time, and could be used to define ice jam alarm conditions based on the ice deceleration and/or ice cover.
Acknowledgements The authors wish to thank M. Rousseau and M. Cauchon for their help on the field. This work was funded by the Fonds de recherche Nature et technologies du Québec (FQRNT).
Adrian, R.J., 2005. Twenty years of particle image velocimetry. Experiments in Fluids 39, 159–169. http://dx.doi.org/10.1007/s00348-005-0991-7. Bourgault, D., 2008. Shore-based photogrammetry of river ice. Canadian Journal of Civil Engineering 35, 80–86. Bradley, A.A., Kruger, A., Meselhe, E.A., Muste, M.V.I., 2002. Flow measurement in streams using video imagery. Water Resources Research 38, 1315. Chave, Lemon, Fissel, Dupuis, Dumont, 2004. Real-time measurements of ice draft and velocity in the St. Lawrence River. Presented at Oceans ’04, MTS/IEEE/TechnoOceans'04, Kobe, Japan, November 2004. Dramais, G., Le Coz, J., Camenen, B., Hauet, A., 2011. Advantages of a mobile LSPIV method for measuring flood discharges and improving stage–discharge curves. Journal of Hydro-Environment Research 5, 301–312. http://dx.doi.org/10.1016/ j.jher.2010.12.005. Emond, J., Morse, B., Richard, M., Stander, E., Viau, A.A., 2011. Surface ice observations on the St. Lawrence River using infrared thermography. River Research and Applications 27, 1090–1105. http://dx.doi.org/10.1002/rra. Ettema, R., Fujita, I., Muste, M., Kruger, A., 1997. Particle-image velocimetry for whole-field measurement of ice velocities. Cold Regions Science and Technology 26, 97–112. Ferrick, M.G., Weyrick, P.B., Hunnewell, S.T., 1992. Analysis of river-ice motion near a breaking front. Canadian Journal of Civil Engineering 19, 105–116. Fujita, I., Muste, M., Kruger, A., 1998. Large-scale particle image velocimetry for flow analysis in hydraulic engineering applications. Journal of Hydraulic Research 36 (3), 397–414. Hauet, A., Creutin, J.-D., Belleudy, P., 2008. Sensitivity study of large-scale particle image velocimetry measurement of river discharge using numerical simulation. Journal of Hydrology 349, 178–190. Jain, A.K., 1989. Fundamental of Digital Image Processing. Prentice Hall, Englewood Cliffs, NJ. Jasek, M., Marko, J., 2008. Acoustic detection and study of frazil ice in a freezing river during the 2007–2008 winter. Proceedings of the 19th IAHR Symposium on Ice, July 2008. Vancouver, British Columbia. Jasek, M., Muste, M., Ettema, R., 2001. Estimation of Yukon River discharge during an ice jam near Dawson City. Canadian Journal of Civil Engineering 28, 856–864. Kantoush, S.A., Schleiss, A.J., Sumi, T., Murasaki, M., 2011. LSPIV implementation for environmental flow in various laboratory and field cases. Journal of HydroEnvironment Research 5, 263–276. http://dx.doi.org/10.1016/j.jher.2011.07.002. Leclerc, M., Morse, B., Francoeur, J., Heniche, M., Boudreau, P., Secretan, Y., 2001. Analyse de risques d'inondations par embâcles de la rivière Montmorency et identification de solutions techniques innovatrices. INRS-Eau R577 and Université Laval Reaserach Report (118 pages). Matte, P., Secretan, Y., Morin, J., 2009. A 2D-transient hydrodynamic model of high spatial resolution of the St. Lawrence fluvial estuary: description and calibration. Poster at the 62nd Annual Congress of the Canadian Water Resources Association, 9–12 June, 2009, Quebec, Canada. Morse, B., Hessami, M., Bourel, C., 2003. Characteristics of ice in the St. Lawrence River. Canadian Journal of Civil Engineering 30, 766–774. Muste, M., Fujita, I., Hauet, A., 2009. Large-scale particle image velocimetry for measurements in riverine environments. Water Resources Research 44, W00D19. http://dx.doi.org/10.1029/2008WR006950. Secretan, Y., 2010. H2D2 Software. http://www.gre-ehn.ete.inrs.ca/H2D2 (accessed September 5, 2012). Shaw, G.C., Kaczkowski, V., 2009. Voie maritime du Saint-Laurent. L'Encyclopédie canadienne. Fondation Historica (http://www.thecanadianencyclopedia.com, accessed November 2009). Tardif, P.M., 2001. Ice images processing interface for automatic features extraction. Medical Imaging 2001: Image processing. The Mathworks Inc., 2012. Image Processing Toolbox User's Guide.
Contents lists available at SciVerse ScienceDirect
Cold Regions Science and Technology journal homepage: www.elsevier.com/locate/coldregions
A methodology based on Particle image velocimetry for river ice velocity measurement Anik Daigle a,⁎, Francis Bérubé b, Normand Bergeron b, Pascal Matte b a b
Cégep Garneau, 1600 de l'Entente, Québec (Québec), Canada, G1S 4S3 Institut national de la recherche scientifique, 490, rue de la Couronne, Québec (Québec), Canada G1K 9A9
a r t i c l e
i n f o
Article history: Received 25 October 2012 Accepted 18 January 2013 Keywords: River ice Particle image velocimetry Remote sensing
a b s t r a c t The in-field use of Particle image velocimetry (PIV), or Large-scale PIV, may be complicated by implementation issues related to the movement of the camera, the variations in luminosity and the distortion of the images due to the inclination of the shore-based camera. The correction for the distortion of the images, or orthorectification, is a critical issue that usually requires the georeferencing of control points located on the river surface. In drifting ice conditions, control points cannot be identified directly on the river surface, and water levels are often unstable. A robust, flexible and relatively inexpensive methodology based on Large-scale PIV for the measurement of river ice velocity is presented. Two orthorectification techniques are proposed, that do not require georeferencing of control points directly on the river surface. The proposed methodology also includes image enhancement, correction for camera movements and a segmentation technique that allows dynamical quantitative estimations of the ice cover extent. The methodology is tested on the St. Lawrence and Montmorency rivers, two rivers of very different sizes and ice conditions. The correction for camera movements proves to be quite reliable and robust to illumination variations and to the presence of shade patterns. The two orthorectification techniques allow to estimate pixel physical dimensions with 6% and b 1% precisions, respectively. Velocities obtained for the St. Lawrence case are compared to velocities simulated by a numerical model. There is strong agreement in velocity values in moderate ice coverage conditions. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Understanding the processes of ice cover formation and ice transport is essential to the prediction of ice jams and flooding, and for the continuation of maritime transport in the winter. In Canada, nearly 44 million tons of goods travel through the St. Lawrence waterway each year (Shaw and Kaczkowski, 2009). Since 1996, the Canadian Coast Guard monitors the ice cover of the St. Lawrence River at various strategic locations using cameras, sonars and radars. This observation protocol allows experts to produce periodic ice condition reports. This system, however, works through considerable human and material resources that small municipalities struggling with seasonal ice jams usually cannot afford. In most cases, the analysis of ice movement is based on a qualitative visual assessment. The quantitative monitoring of river ice movement is difficult and usually relies on acoustic Doppler profilers providing velocity values at one point on the river (Chave et al., 2004). These are 0-D measurements, while the understanding of the formation of ice cover,
⁎ Corresponding author. Tel.: +1 418 654 2596; fax: +1 418 654 2600. E-mail addresses: [email protected] (A. Daigle), [email protected] (F. Bérubé), [email protected] (N. Bergeron), [email protected] (P. Matte). 0165-232X/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.coldregions.2013.01.006
transport, as well as the estimation of the amount of flowing ice, requires knowledge of the whole velocity field. Moreover, these are expensive immersed instruments that can be damaged by the presence of ice, and there may be important installation and operation safety concerns in flood conditions and near structures like bridges or hydroelectric plants. These limitations call for the establishment of a remote and twodimensional measurement method, two conditions that are met by Particle image velocimetry (PIV). PIV is an imaging technique that was adapted in the 90s for the measurement of the two-dimensional flow velocity field in laboratory experiments (Adrian, 2005; Bradley et al., 2002; Ettema et al., 1997). The stream is firstly seeded with small floating particles (such as plastic beads or rice grains) that act as surface flow tracers. The surface flow is then filmed, with a video camera or by taking successive images with a still camera. The velocity field is obtained by identifying subareas of the image, also called Interrogation Areas (IA, after Muste et al. (2009)), in two successive frames: the correlation coefficients between a given IA of the image acquired at time t and some neighboring area in the image acquired at time t + Δt are calculated. The new t + Δt coordinates of the IA are selected as those corresponding to the most highly correlated t + Δt subarea, and the velocity vector is obtained by dividing the displacement of the IA by Δt. This procedure, repeated on IAs covering the entire image and for all time steps, results in a dynamical and two-dimensional velocity vector field.
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
PIV thus allows the measurement of the spatial distribution of the surface velocity of a seeded flow in time, information that is crucial in a variety of practical applications such as the dispersion of pollutants, bank erosion, sedimentation, the local effects of flooding, and the interaction between the flow and riverine structures. Knowledge of the spatial distribution of water and/or ice velocity also gives a fundamental insight into the hydrodynamics of rivers, providing reliable data for calibration and validation of numerical models. The in-field use of PIV (also known as Large-scale PIV or LSPIV) was initiated by Fujita et al. (1998), who presented a series of adjustments required to move from laboratory and small-scale applications to in-field and larger scale uses. Among other parameters to be accounted for, the authors mention the variable outdoor luminosity and distortion of the images due to the inclination of the shore-based camera. LSPIV has since been used for the measurement of water surface velocities in various laboratory and field conditions (Bradley et al., 2002; Kantoush et al., 2011; Muste et al., 2009), including during floods (Dramais et al., 2011). The potential of LSPIV for field measurement of river ice velocity has however practically not been exploited. A method for the spatial description of ice velocity from video images was presented by Ferrick et al. (1992), but this method involved manual monitoring of individual pieces of floating ice (particle tracking velocimetry, or PTV). The first study explicitly involving LSPIV for the measurement of ice velocity was presented by Ettema et al. (1997). LSPIV was used in a laboratory model of the confluence of the Mississippi and Missouri rivers, to estimate the velocity field and the flow of ice, simulated by polypropylene particles. The first field LSPIV experiment applied to river ice was presented by Jasek et al. (2001), who estimated the flow of the Yukon River during an ice jam near Dawson City, using the ice floes as tracers of the flow. Tardif (2001) presented a method to extract ice velocity from images using a method similar to LSPIV. The author however pointed out several implementation issues, including the reduced visibility during snowfalls, the low contrast in the images at sunrise and sunset, and the presence of reflections on the water surface complicating the recognition of IAs from one image to another. Similar problems, as well as unavoidable movements of the camera, were reported by Bourgault (2008), who attempted to use LSPIV to obtain ice velocity fields on the St. Lawrence River in Québec City. The evaluation of ice flow depends on three measurements: the drifting ice velocity field, the ice cover extent, and the ice thickness. This paper presents a robust, flexible and relatively inexpensive methodology based on LSPIV and image segmentation for the two first measurements. Results obtained for two case studies are presented. The LSPIV methodology is presented with particular focus on three major issues pointed out by preceding authors for the in-field implementation of LSPIV for the measurement of ice velocity: the correction for camera movements; the enhancement of images in order to account for low contrast, variable luminosity and reflections; and correction for the image distortion due to the inclination of the camera, also called orthorectification. This latter critical issue usually requires land surveying equipment and expertise, and involves the georeferencing of control points on the surface to be orthorectified. In drifting ice conditions, control points cannot be identified directly on the river surface, and water levels are often unstable. Two alternative orthorectification methods are presented: one that uses control points that do not need to be defined directly on the river surface and thus may be georeferenced a posteriori, and another that does not require georeferencing at all and is particularly suited for small rivers. A segmentation technique that allows dynamical quantitative estimations of the ice cover extent is also presented. The technique is based on the classification of the image pixels in the two “ice” and “not ice” categories, based on image texture and color information.
37
2. Study sites and data acquisition Two study sites were selected for the development and testing of the methodology. They were chosen to be different in sizes and ice conditions, in order to test the flexibility of the proposed method. The St. Lawrence River near Québec City is about 1-km wide, but narrows to ~ 600 m between the Québec and Pierre-Laporte bridges. Tides have a maximum amplitude of ~ 5 m, and water velocity can reach 3 m s −1 in both directions. Ice forms on the river every winter. The waterway and maritime transport are maintained by the Canadian Coast Guard throughout the winter, but a steady stream of ice floes generally occurs from December to March. Images were acquired from the South-shore, at an elevation of ~35 m above the water level (Fig. 1a). This site was chosen because of its accessibility, a relatively high elevation above the water level allowing a more perpendicular line of sight, and the presence of permanent structures in the field of view (e.g. the bridge pillars) that facilitate image georeferencing.
Fig. 1. St. Lawrence River, between the Québec and Pierre-Laporte bridges. (a) The three control regions used for the correction of camera movements, as found in a reference image (dashed boxes). (b) Locations of control region 1 in the image to be corrected (solid line box), and in the reference image (dashed box). (c) Control region 1 in the transformed (translated and rotated) image.
38
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
The Montmorency River flows into the St. Lawrence River about 15 km north-east of Québec City. Its drainage basin is about 1100 km 2. An ice cover forms from December yearly, and the river has experienced several ice jams in its history (Leclerc et al., 2001). The images were acquired from a bridge crossing the river about 10 km from its mouth, where the river is about 50-m wide. All images were acquired with a Canon Rebel T2i DSLR camera with 18 Mpix resolution mounted on a regular tripod. In burst mode and using a 30 Mb s −1 memory card, sequences were recorded at a rate of 3.4 images s −1. 3. Methods All development was conducted using Matlab R2011b. 3.1. Image stabilization Small camera movements due to wind and/or oscillations of the structure on which the camera is based can rarely be completely avoided in the field. The effect of such movements is a small continuous variation of the field of view of the camera creating shifts in the coordinates of the IAs to be tracked by the PIV algorithm. These shifts are not necessarily the same throughout the image, since the camera may experience small rotations in addition to translations. To ensure that the shift of a given IA is solely due to a real displacement of the water surface, every individual image of the sequence must be framed and scaled identically. To achieve this, one frame of the sequence is selected as the reference image and all other images are translated, rotated, and reframed according to this reference. This image transformation procedure, also called image registration, is performed for each image by the multiplication of a transformation matrix. A transformation matrix is a 3 × 3 matrix that specifies the required translation, rotation, shear and/or scaling operations for an image to be aligned to the reference image. Such a matrix can be inferred from a minimum of two control points pairs identified in the reference image and in the image to be transformed. Control points have to be fixed recognizable locations, such as buildings or reference poles. Since the measurement of ice velocity may often be needed in natural areas where no buildings are found, the proposed method does not use control points but control regions in the images that may be parts of the natural landscape such as tree trunks, shore ground or rocks. These regions must be as still as possible (unaffected by the flow or the wind) and located at opposite sides of the image in order to detect a possible rotation of the image (Fig. 1). Once control regions are identified on the reference image by the user (Fig. 1a), they are detected in all images of the sequence by an algorithm based on the correlation of the regions with some neighboring area — very much as in the PIV algorithm. The size of the control regions is not critical, but a larger region is more likely to be uniquely positioned in each image. A rather textured region also helps getting a precise and unique location in every image. The control points are then defined as the center of the detected regions. One transformation matrix is defined for each individual image, and each image is transformed accordingly. The result is a sequence for which the control regions maintain the same coordinates, with the expense of small peripheral losses (Fig. 1c). 3.2. Image enhancement Other unavoidable issues in outdoor implementation of PIV are the variation in luminosity and contrast due to clouding, reflections, and sometimes shade, on the water surface. There are several ways to enhance the quality of an image, often involving image filtering. The tool developed in this study makes use of an adaptive low-pass (wiener) filter for the subtraction of noise, and of an adjustment of
the intensity distribution of the image pixels for the improvement of contrast. The variation in luminosity due to the passage of clouds is dampened by the subtraction, for every image, of the mean pixel intensity value. The major effects on the images of this study were however due to reflections and shades on the water surface (e.g. Fig. 1a). An IA travelling through the limits of a reflection or a shade pattern is more difficult to track with precision. Reflections can be substantially reduced by the use of a polarizing filter. Remaining reflections and shades can be further reduced by the subtraction of a local background containing the shade and reflection patterns. For details about the image enhancement techniques that were used, see the wiener2 function (adaptive low-pass filtering of noise), the imadjust function (for the adjustment of the image intensity distribution) and the imtophat function (for the calculation of a local background) of the Matlab Image Processing Toolbox (The Mathworks Inc., 2012). 3.3. Image orthorectification PIV in itself provides velocity in pixels per second. To convert to physical space velocity, the true physical dimensions of the image pixels must be known. In oblique images, such as images taken from the shore, the pixel physical dimensions vary throughout the image. Correction for the oblicity of an image is called orthorectification. Typically, orthorectification is done by georeferencing control points in the image that define a transformation matrix. This matrix assigns new coordinates to every pixel of the image so that they correspond to the physical space. This procedure is relatively straightforward when the user has access to land surveying equipment and expertise, and when fixed georeferenced points can be found directly on the surface to be orthorectified. In the present application however, the surface is liquid, flowing and may even be moving up and down because of rapid water level changes due to local ice accumulations or tides. We thus developed two different orthorectification techniques adapted to the in-field implementation of PIV. The first one is based on an above-surface georeferencing, and makes use of available georeferenced points that do not need to be on the water surface. The other does not use georeferencing at all, and rather relies on the interpolation of several scale measurements provided by an object of known dimensions that is introduced into the flow during the image acquisition. 3.3.1. Above-surface georeferencing Georeferencing can often be facilitated by the presence of permanent fixed structures. Images of the St. Lawrence River section located between the Québec and Pierre-Laporte bridges, in Québec City, include bridges pillars and streetlights (Fig. 2). These are fixed points that are easily georeferenced with regular surveying equipment. Orthorectification based on above-surface georeferecing can be performed following these steps: 1. Georeference fixed points above the water surface level, and one at the water surface level. In the example shown in Fig. 2, 15 points were georeferenced (crosses). Point 15 was georeferenced at the water surface level, on one of the Québec bridge pillar. 2. Determine height scales at many locations in the field of view. In this study, height scales were determined using available information about the dimensions of the bridge pillars bricks and heights of the streetlight poles. 3. Relocate the georeferenced points on the river surface, using the local scale heights and the position of the water surface level. In Fig. 2, circled crosses mark the positions of the 15 georeferenced fixed points, when relocated on the river surface.
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
39
Fig. 2. Fixed points (crosses) georeferenced above the St. Lawrence surface and relocated (circled crosses) at the water surface level.
4. Use the relocated points as control points for the definition of the transformation matrix, and proceed to the image orthorectification. 3.3.2. Interpolation of local dimension references In most situations, especially in remote areas, permanent fixed points do not exist. One possibility is to install landmarks such as metal posts near the water surface and to proceed with the regular georeferencing and orthorectification procedure. In many situations, however, the installation of landmarks can be time consuming and difficult: for the transformation matrix to be reliable, the markers must be positioned near the limits of the field of view, and thus on each sides of the river. In flood conditions and/or with the presence of ice, one side of the river may be difficult to access and the shores can be unstable. In such urgent situations as the monitoring of a possible ice jam, these issues related to orthorectification may thus prevent the use of PIV. Another way of inferring pixel dimensions is to introduce a floating object of known dimensions within the field of view of the camera during image acquisition (Fig. 3). In the example presented here, the object is a floating square plate. The method can be summarized by these steps: 1. The plate is photographed or filmed at many locations in the field of view (Fig. 3); 2. Coordinates of the plate corners are identified for each location, defining squares of known dimensions and a local orthogonal reference system, and providing physical dimensions to the pixels of the squares (Fig. 4); 3. The local horizontal and vertical pixel dimensions obtained at each plate location are interpolated on the river surface. 3.4. Estimation of the ice cover extent The estimation of the ice cover extent is necessary for the quantitative evaluation of the amount of flowing ice, a variable crucial for the prediction of ice jams (Tardif, 2001). The PIV methodology developed here includes an evaluation of the ice cover for every image of a sequence, by the automatic recognition of ice pixels in the image. In image processing, the separation of an image into different parts is called segmentation. There are many different ways to perform segmentation, such as edge detection, texture and color classification
(Jain, 1989). The choice of the segmentation technique depends on the problem that is addressed. In a single river ice image, the ice can present very different colors and textures. In Fig. 5a, it is apparent that the texture of large snow covered ice floes in the foreground differs from that of the smaller ice floes. This difference is further increased by the effect of perspective. A segmentation based on color is not more reliable, with the shade created by clouds or elevated objects (here the Pierre-Laporte Bridge) causing abrupt color changes, and reflections making parts of the free water surface appear almost as pale as ice floes. For these reasons, segmentation based solely on image texture or pixel colors cannot provide reliable estimates of the percentage of ice cover. The combination of both, however, resulted in a clear separation of free water/ice pixels and proved to be especially robust when tested on different images (Fig. 5b). The procedure adopted uses an artificial neural network (ANN) to classify each pixel of the river surface in the two “ice” and “not ice” classes, based on texture and color information assigned to each individual pixel. It can be summarized as follows: 1. The image to be segmented is orthorectified; 2. The texture information is given by a measure of the image local entropy: for each pixel of the image, the entropy is computed on a 9 × 9 pixel neighborhood; 3. The color information is represented by the three-element RGB color array of each pixel; 4. The user identifies 25 “ice” pixels and 25 “not ice” pixels by clicking on the image; 5. These 50 pixels are used to calibrate an ANN that learns to separate “ice” and “not ice” pixels; 6. The calibrated ANN classifies all pixels on the river surface; 7. The ice cover percentage is computed as the sum of “ice” pixels over the total number of pixels on the river surface. The ANN is a three-layer perceptron that classifies each pixel of the river surface based on the four-element [entropy-RGB] inputs assigned to each individual pixel. The user identifies the 50 calibration pixels only once for a given sequence, and the segmentation is performed automatically for every frame, resulting in a dynamical estimation of the ice cover. Illumination variations are fairly well supported by the segmentation procedure, as an average background
40
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
Fig. 3. Montmorency River. A square plate of known dimensions is attached to a fishing line and photographed at several locations in order to provide the scaling information.
is subtracted to each image previous to the ANN classification. Large changes in the illumination or weather conditions however require a new calibration for the segmentation to remain reliable. 4. Results 4.1. Image stabilization and enhancement Image stabilization and enhancement were performed on the Montmorency and St. Lawrence images using the same procedure. Conditions near the St. Lawrence River are most of the time windy, and sequences acquired at the St. Lawrence site were more affected by camera movements than sequences acquired at the Montmorency site. An example of the effect of the stabilization routine is shown in Fig. 6. It shows the correlation coefficients between a control region as found in the reference image and in the unprocessed images (gray line), as well as in the stabilized images (black line). In this example sequence, 1530 images were acquired at 3.4 Hz between 15:34 and 15:41
on February 19th, 2012. Three control regions were defined as shown in Fig. 1a, and the reference image was chosen half-way through the sequence. Correlation coefficients for control region 1 are shown in Fig. 6a. The smallest coefficients between the reference image and the unprocessed images (gray line) correspond to the largest camera shifts. Correlation coefficients between control region 1 as found in the reference image and in the stabilized images are all >0.99. Correlation coefficients for control region 2 are shown in Fig. 6b. Here again, the coefficients increase drastically after the stabilization, but decrease progressively on both sides of the reference image time. This is due to the fact that no background subtraction was performed for this example, and that control region 2 is defined in an area that is affected by the moving of the shade of the Pierre-Laporte Bridge (Fig. 1a). Because control region 2 is well textured, it could still be reliably detected by the stabilization routine throughout the sequence despite this variable shade pattern. Control regions should nevertheless be chosen so as to minimize such alterations, or the performance
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
41
Fig. 4. Local orthogonal system allowing the inference of the local x- and y-pixel dimensions on the Montmorency River surface.
of the stabilization routine could be affected despite image enhancement efforts. The stabilization of one 18Mpix image takes ~ 33 s using Matlab running on an Intel Core Duo CPU [email protected] GHz with a RAM of 1.94 Go. The increase in the precision of the IA displacement measurement due to each of the image enhancement steps was measured by the mean increase in the signal-to-noise ratios of the IA correlation matrices. Our well-contrasted and low-noise images do not contain enough noise for the wiener filter, or the improvement of the contrast, to have a noticeable impact on the SNR ratios. The background subtraction, on the other hand, allows a substantial reduction of the contribution of reflections and shades to the images background and doubles the average SNR of the IA correlation matrices.
The six locations of the floating plate that were used to infer pixel sizes on the river surface are shown in Fig. 7. Two horizontal and two vertical reference scales are defined for each of the plate locations and interpolated on the river surface. Typically, on images acquired from the shore, pixels dimensions vary
4.2. Interpolation of local dimension references: the case of the Montmorency River The case of the Montmorency River was considered mostly to test the orthorectification technique based on local dimension references.
Fig. 5. Segmentation of the St. Lawrence surface between the Québec and Pierre-Laporte Bridges, based on the combination of texture and pixel color information. Image was taken on March 17th, 2011, at 14:42.
Fig. 6. (a) Correlation coefficients between the control region 1 as found in the reference image of the St. Lawrence River and in the unprocessed images (gray line) and the stabilized images (black line), as a function of time. Position of the reference image is indicated. (b) Correlation coefficients for control region 2, as a function of time. Coefficients for control region 3 are not shown. No image enhancement has been applied for this example.
42
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
Fig. 7. Locations of the floating plate that were used to infer pixel sizes on the Montmorency River surface. The reference lengths used to compute the local x- and y-pixel scales are indicated.
mostly in the image foreground-background direction (direction y in Fig. 7). Pixels horizontal and vertical dimensions were interpolated in the x and y directions by distinct linear and quadratic functions (Fig. 8). Performances of both interpolation types are evaluated by measuring the four edges of the plates at each of the six locations. The 24 lengths obtained are compared to the actual 61 cm length of the plate edges. Performance is evaluated by calculating the root mean square error (RMSE): rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 XN 2 RMSE ¼ ðL −LÞ i¼1 i N
ð1Þ
where L is the actual edge length (=61 cm), Li are the lengths as measured on the image based on the pixel size interpolation and N is the number of measurements (24). The RMSE achieved using a linear interpolation is 5.1 cm, while the quadratic interpolation achieves a RMSE = 3.7 cm. Velocities in cm s−1 can now be obtained on the water surface. Fig. 9 shows the average velocity vector field obtained on a 5-minute sequence acquired on December 10th, 2011. The instantaneous velocity components averaged on the subarea delimited in Fig. 9 are plotted as a function of time in Fig. 10. Ice cover of the whole river surface as a
function of time is also shown in Fig. 10. The right hand side column of Fig. 10 shows the distribution of the instantaneous velocity components values and estimated ice cover values during the 5-minute sequence. Velocity measured on the 5-minute sequence has a mean modulus of 30 cm s−1, with a 1.8 cm s −1 standard deviation. 4.3. Above-surface georeferencing: the case of the St. Lawrence River The orthorectified version of Fig. 2 is shown in Fig. 11, where the 15 control points are plotted at their orthorectified coordinates. Points 1–2 and 12–13–14–15 now appear superposed, since these share approximately the same coordinates on the river surface plane (see Fig. 2). The orthorectified coordinates of each point were compared to their georeferenced coordinates in order to evaluate the error associated to the orthorectification. The error on the orthorectified coordinates is evaluated by the RMSE and rRMSE, or relative RMSE, which is given by: rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1 XN C i −C i i¼1 rRMSE ¼ N C max −C min
ð2Þ
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
43
whole river surface. Velocity measured for this 7-minute sequence has a mean modulus of 115 cm s−1, with a 5 cm s−1 standard deviation. 4.4. Uncertainty evaluation
Fig. 8. Quadratic interpolation of (a) the x- and (b) y-pixel dimensions of the Montmorency River image.
where Ci and C i are the ith georeferenced and estimated coordinates, and Cmin and Cmax are the minimum and maximum coordinates of the 15 points considered. This relative performance measure allows the rescaling of the coordinate values and a more representative description of the relative estimation precision. In the example shown in Fig. 11, the RMSEs obtained for the 15 transversal and longitudinal coordinates are 5.2 m and 0.63 m, respectively, which represent rRMSEs of 0.7% and 0.2%. The average velocity vector field obtained from a 7-minute sequence acquired on February 19th, 2012 is shown in Fig. 12. The instantaneous velocity components averaged on the subarea delimited in Fig. 12 are plotted as a function of time in Fig. 13, along with the ice cover of the
Velocities measured on the 5- and 7-minute sequences acquired on the Montmorency and on the St. Lawrence (Figs. 10 and 13) show respective coefficient of variation of 6% and 4%. It is difficult to evaluate which part of this variability is real, and which is attributable to measurement error. Several sources of uncertainty are associated to LSPIV that have been discussed in previous studies (Dramais et al., 2011; Hauet et al., 2008). The main sources of uncertainty identified for the particular applications presented here are (1) the orthorectification, and (2) the measurement of the IA displacement. Errors associated to both orthorectification techniques are systematic and thus affect the mean velocity values and not their variability in time. The evaluation of the IA displacement, on the other hand, is made for every time step. As stated earlier, this measurement is based on the computation of a correlation matrix between subareas of an image taken at time t and another taken at time t + Δt. This matrix is typically highly peaked at the position of the IA in the t + Δt frame. To achieve a sub-pixel evaluation of this position, the autocorrelation matrix is modelled by a two-dimensional Gaussian profile. The precision of the displacement measure is thus dependant on the quality of this fit, which in turn depends on the IA size, time step, and on the image contrast and SNR. The measurement of the IA displacement also depends on the image stabilization, also based on the computation and modeling of autocorrelation matrices. The image transformation, which involves an interpolation, also adds to the IA displacement uncertainty. In our 18 Mpix resolution well-contrasted and low-noise images, the uncertainty associated to the measurement of the IA displacement is evaluated to be ±1 pixel. The measured IA displacements being of the order of 10 pixels for both the Montmorency and St. Lawrence cases, the associated relative error can thus reach 10%, which is higher than the observed coefficients of variations. 4.5. Comparison of the St. Lawrence River velocities to H2D2 simulations Direct measurements of ice velocity are difficult and no direct validation of the obtained velocity values could be performed. The velocity values obtained using the PIV methodology on image
Fig. 9. Average velocity vector field obtained on a 5-minute image sequence of the Montmorency River acquired on December 10th, 2011.
44
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
Fig. 10. x- and y-velocity components averaged on the subarea delimited on Fig. 9, and ice coverage estimated on the whole river surface, plotted as a function of time (on the left), and as histogram distributions (on the right).
sequences of the St. Lawrence were however compared to velocities simulated by a numerical model. A 2D high-resolution numerical model of the St. Lawrence fluvial estuary (Matte et al., 2009) was used to simulate the hydrodynamic conditions prevailing during the analysis periods. The finite-element grid underlying the model has an average spatial resolution of 50 m, with refinements down to 1 m around engineering structures such as bridge pillars. The 2D vertically-integrated shallow-water equations are solved using the H2D2 software (Secretan, 2010). The model was calibrated under average flow conditions (~12,000 m 3 s −1) for a neap-spring tide cycle with an accuracy of b5 cm in water levels at most stations. However, simulated water velocities were not verified against measured data, so that deviations between modeled and measured lateral velocity distributions may appear, mainly where cross-sectional variations in bed friction are significant. The accuracy of results for other flow conditions was not verified either. For the purpose of this study, the domain consists of a 30-km reach of the river, with boundaries positioned at the upstream and downstream tide gauges closest to the Quebec Bridge, respectively located at Neuville and Quebec. The model is forced at the boundaries with water levels
measured at the tide gauges. Main tributaries, also influenced by tides, are not included in the model due to a lack of data at their outlets; this may affect the modeled velocities only within the radius of influence of the tributaries as their discharges are quite small compared to the St. Lawrence. Simulations were performed without ice since no information on the extent, thickness and roughness of the ice cover was available for the system at the simulation times. Velocities simulated by the numerical model and measured by the PIV methodology on the same transect on February 19th, 2012, show strong agreement, as shown in Fig. 14a. The same comparison was performed on the same transect of the St. Lawrence but for a sequence acquired on March 17th, 2011. Results are shown in Fig. 14b. This time, velocities simulated by H2D2 are more than twice as high as the velocities measured by the PIV methodology.
5. Discussion The main processing difference in the two case studies presented is in the orthorectification techniques that are proposed.
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
Fig. 11. Orthorectified image of the St. Lawrence River between the Quebec and Pierre-Laporte bridges.
The orthorectification based on the interpolation of local dimension references is proposed as an alternative to the traditional georeferencing of control points located on the surface to be orthorectified. This alternative technique is based on the introduction
45
of a floating object of known dimension in the flow, and is thus particularly suited for situations where land surveying material is not available, and/or when deployment time is a critical issue. This technique does not require georeferencing, is particularly affordable, can handle rapid water level variations and takes almost no time to deploy. It can be handled from one shore only if the river width so permits. In the Montmorency River case that is presented, best results are obtained with a quadratic interpolation of the local dimension references. This small nonlinear trend, more pronounced in the horizontal direction, is probably due to the fact that no correction was made to account for the camera lens distortion. This technique is more suited to small to average sized rivers: large rivers involve large distances, over which the position of the floating object is more difficult to control. The use of a floating object on large rivers would also require the object to be large enough for it to remain visible on the images even when located near the shore opposite to the camera. A floating object of proper dimension attached to a boat could allow the use of this technique on large rivers. Above-surface georeferencing differs from traditional georeferencing by the fact that it is based on available fixed points that are not necessarily positioned on the river surface. Georeferencing may thus be performed after the image acquisition, a significant advantage in urgent situations, or when the water level is changing. The methodology corrects for the height of the control points given known height references in the image, and then proceeds to a regular image registration. The validation of the obtained velocity values remains a challenging problem. Direct measurements (for example, ADCP measurements) could provide reliable validation values. H2D2 simulations of the St. Lawrence velocities nevertheless provide an insightful comparison basis. Velocities obtained by both the PIV-based method and H2D2 agree very closely for the 2012 sequence, but show an important but quasi-constant difference for the 2011 sequence. An interesting hypothesis is that this discrepancy is due to the very extensive ice cover present on the river on March 17th, 2011: ice cover values varied between 75% and 80% with large ice floes in contact (Fig. 5), while the ice cover on the February 2012 sequence did not exceed 50% (Fig. 11). Surface velocity on March 17th, 2011 was thus probably substantially lower than the water-column averaged velocity. The fact that the numerical model does not account for the presence of
Fig. 12. Average velocity vector field obtained on a 7-minute image sequence of the St. Lawrence River acquired on February 19th, 2012.
46
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
Fig. 13. x- and y-velocity components averaged on the subarea delimited on Fig. 12, and ice coverage estimated on the whole river surface, plotted as a function of time (on the left), and as histogram distributions (on the right).
ice also certainly explains part of the observed discrepancy. More comparisons on more sequences showing a range of ice cover values would be needed to verify this hypothesis. The proposed method includes an estimation of the ice cover as a function of time based on image segmentation. The ice cover values obtained on the Montmorency sequence and presented in Fig. 10 show a number of sudden jumps to ~ 75%. This is due to the important presence of steam near the water surface: the water vapor takes a coloration and texture similar to the ice under particular illumination conditions, and this gets to trick the segmentation routine for a few frames at a time. These singular ice cover values can be discarded visually based on the histogram distribution, or automatically filtered. Other specific conditions, such as a snowfall, would probably have similar consequences. Additional tests in more challenging illumination conditions (e.g. during a snowfall or in foggy weather) would allow to better evaluate the effects of image enhancement on the reliability of the obtained velocity and ice cover values. A river discharge can be obtained by multiplying the area of a given cross section and the average velocity of the water through this section. By providing the surface velocity information, PIV could allow a real-time evaluation of river discharge, given the bathymetry
information for the observed transect and an independent measurement of the water level elevation. The segmentation procedure presented here is a promising first step towards a dynamical evaluation of the amount of flowing ice using a PIV-based methodology. The conversion from ice cover to ice volume requires an estimation of the ice thickness, that cannot be obtained from optical images such as the ones presented here. Studies however showed that the ice thickness can be measured by an upward looking sonar instrument (Jasek and Marko, 2008; Morse et al., 2003) or be inferred from thermal images (Emond et al., 2011). 6. Conclusion A methodology based on PIV for the measurement of river ice velocity was presented and tested for two rivers of very different sizes. It was developed to be a simple and affordable implementation procedure allowing organizations and users with limited resources to monitor drifting ice velocity and ice cover. Frequently mentioned implementation problems such as camera movements, variable luminosity, reflections, shades and orthorectification are taken into account. The methodology allows for obtaining velocity fields and ice cover estimations as a function
A. Daigle et al. / Cold Regions Science and Technology 89 (2013) 36–47
47
References
Fig. 14. Velocities obtained on the same transect of the St. Lawrence River by the PIV methodology (gray line) and the H2D2 numerical simulation (black line) on (a) February 19th, 2012 (~14 h30) and (b) on March 17th, 2011 (~14 h45).
of time, and could be used to define ice jam alarm conditions based on the ice deceleration and/or ice cover.
Acknowledgements The authors wish to thank M. Rousseau and M. Cauchon for their help on the field. This work was funded by the Fonds de recherche Nature et technologies du Québec (FQRNT).
Adrian, R.J., 2005. Twenty years of particle image velocimetry. Experiments in Fluids 39, 159–169. http://dx.doi.org/10.1007/s00348-005-0991-7. Bourgault, D., 2008. Shore-based photogrammetry of river ice. Canadian Journal of Civil Engineering 35, 80–86. Bradley, A.A., Kruger, A., Meselhe, E.A., Muste, M.V.I., 2002. Flow measurement in streams using video imagery. Water Resources Research 38, 1315. Chave, Lemon, Fissel, Dupuis, Dumont, 2004. Real-time measurements of ice draft and velocity in the St. Lawrence River. Presented at Oceans ’04, MTS/IEEE/TechnoOceans'04, Kobe, Japan, November 2004. Dramais, G., Le Coz, J., Camenen, B., Hauet, A., 2011. Advantages of a mobile LSPIV method for measuring flood discharges and improving stage–discharge curves. Journal of Hydro-Environment Research 5, 301–312. http://dx.doi.org/10.1016/ j.jher.2010.12.005. Emond, J., Morse, B., Richard, M., Stander, E., Viau, A.A., 2011. Surface ice observations on the St. Lawrence River using infrared thermography. River Research and Applications 27, 1090–1105. http://dx.doi.org/10.1002/rra. Ettema, R., Fujita, I., Muste, M., Kruger, A., 1997. Particle-image velocimetry for whole-field measurement of ice velocities. Cold Regions Science and Technology 26, 97–112. Ferrick, M.G., Weyrick, P.B., Hunnewell, S.T., 1992. Analysis of river-ice motion near a breaking front. Canadian Journal of Civil Engineering 19, 105–116. Fujita, I., Muste, M., Kruger, A., 1998. Large-scale particle image velocimetry for flow analysis in hydraulic engineering applications. Journal of Hydraulic Research 36 (3), 397–414. Hauet, A., Creutin, J.-D., Belleudy, P., 2008. Sensitivity study of large-scale particle image velocimetry measurement of river discharge using numerical simulation. Journal of Hydrology 349, 178–190. Jain, A.K., 1989. Fundamental of Digital Image Processing. Prentice Hall, Englewood Cliffs, NJ. Jasek, M., Marko, J., 2008. Acoustic detection and study of frazil ice in a freezing river during the 2007–2008 winter. Proceedings of the 19th IAHR Symposium on Ice, July 2008. Vancouver, British Columbia. Jasek, M., Muste, M., Ettema, R., 2001. Estimation of Yukon River discharge during an ice jam near Dawson City. Canadian Journal of Civil Engineering 28, 856–864. Kantoush, S.A., Schleiss, A.J., Sumi, T., Murasaki, M., 2011. LSPIV implementation for environmental flow in various laboratory and field cases. Journal of HydroEnvironment Research 5, 263–276. http://dx.doi.org/10.1016/j.jher.2011.07.002. Leclerc, M., Morse, B., Francoeur, J., Heniche, M., Boudreau, P., Secretan, Y., 2001. Analyse de risques d'inondations par embâcles de la rivière Montmorency et identification de solutions techniques innovatrices. INRS-Eau R577 and Université Laval Reaserach Report (118 pages). Matte, P., Secretan, Y., Morin, J., 2009. A 2D-transient hydrodynamic model of high spatial resolution of the St. Lawrence fluvial estuary: description and calibration. Poster at the 62nd Annual Congress of the Canadian Water Resources Association, 9–12 June, 2009, Quebec, Canada. Morse, B., Hessami, M., Bourel, C., 2003. Characteristics of ice in the St. Lawrence River. Canadian Journal of Civil Engineering 30, 766–774. Muste, M., Fujita, I., Hauet, A., 2009. Large-scale particle image velocimetry for measurements in riverine environments. Water Resources Research 44, W00D19. http://dx.doi.org/10.1029/2008WR006950. Secretan, Y., 2010. H2D2 Software. http://www.gre-ehn.ete.inrs.ca/H2D2 (accessed September 5, 2012). Shaw, G.C., Kaczkowski, V., 2009. Voie maritime du Saint-Laurent. L'Encyclopédie canadienne. Fondation Historica (http://www.thecanadianencyclopedia.com, accessed November 2009). Tardif, P.M., 2001. Ice images processing interface for automatic features extraction. Medical Imaging 2001: Image processing. The Mathworks Inc., 2012. Image Processing Toolbox User's Guide.