AUTOMATIC SPATIO-TEMPORAL FLOW VELOCITY MEASUREMENT IN SMALL RIVERS USING THERMAL IMAGE SEQUENCES

An automatic spatio-temporal flow velocity measurement approach, using an uncooled thermal camera, is proposed in this paper. The basic principle of the method is to track visible thermal features at the water surface in thermal camera image sequences. Radiometric and geometric calibrations are firstly implemented to remove vignetting effects in thermal imagery and to get the interior orientation parameters of the camera. An object-based unsupervised classification approach is then applied to detect the interest regions for data referencing and thermal feature tracking. Subsequently, GCPs are extracted to orient the river image sequences and local hot points are identified as tracking features. Afterwards, accurate dense tracking outputs are obtained using pyramidal Lucas-Kanade method. To validate the accuracy potential of the method, measurements obtained from thermal feature tracking are compared with reference measurements taken by a propeller gauge. Results show a great potential of automatic flow velocity measurement in small rivers using imagery from a thermal camera. * Corresponding author


INTRODUCTION
Automatic and accurate measurement of flow velocities is of large importance for hydrology research.As a non-intrusive and continuous implementation, traditional image-based methods usually use particles or dye on the free surface as distinguishable tracers to track.These particles may be tracked by PTV (particle tracking velocimetry) techniques.Using neutrally buoyant particles, 3D-PTV systems have been realized mainly for laboratory applications (e.g.Maas et al., 1993).While PTV requires detecting and tracking individual particles, PIV (particle imaging velocimetry) employs area-based matching strategies to track groups of particles in densely seeded plows (Adrian, 1991).However, the insertion of particles to mark a flow is often not desired.Besides pollution aspects, a drawback of floating particles is their tendency to agglomerate (Weitbrecht et al. 2002).Alternatively, laserinduced fluorescence (LIF) is a technique based on dyes, which are added to a flow at a very low concentration and visualized by a thin laser lightsheet.The technique can be extended to 3D-LIF by illuminating a volume in a flow (Maas and Gruen, 1995).A drawback herein is the need of a powerful laser light source, which limits outdoor use of the technique.Furthermore, the quality of tracking results in optical imagery is severely affected by illumination conditions.
A possible solution to these mentioned challenges is the use of thermal tracers, which can simply be applied by pouring warm (or icy) water onto a river surface.Taking thermal camera image sequences of the water surface and tracking warm water patterns in thermal imagery is not only environmentally friendly but is also completely independent on illumination.Therefore, flow velocity measurements using thermography has achieved more and more attention.Lima and Abrantes (2014) suggested a combined tracer (heated dye) to measure the overland and rill flow velocity using leading edge estimation in laboratory set-ups.Lima et al. (2015) extended their technique to various field experiments.Tauro and Grimaldi (2017) utilized ice dices as thermal tracers to monitor the flow velocity of a river in Trento, Italy.However, these techniques can only measure the mean flow velocity of the whole region or very sparsely-distributed velocity field.Therefore, our research focuses on accurate, automatic, and densely-distributed flow velocity measurement using thermal imagery in field studies.
Regarding the usage of thermal signatures in this study, the grid-attached measurements with PIV are less practical due to unevenly distributed textures.Instead, we chose PTV which tracks each of the detected features individually using the pyramidal Lucas-Kanade method (Bouguet, 2001).In order to avoid unreliable tracking results in texture-less regions, a tracking-quality-assessment feature is introduced in this study.In addition, in order to reduce the degree of supervision with respect to the detection of tracking features, an object-based unsupervised classification method is implemented to continuously detect thermal features.

DATA ACQUISITION
The data were acquired at the small river Wernersbach in the Tharandter Wald in Saxony, Germany.The basic system to measure flow velocity consists of an uncooled thermal camera FLIR Ax65.Focal length is 13 mm and pixel size 17 µm.The measurement temperature range is -25°C to +135°C.Maximum acquisition rate is 30 Hz.The camera is positioned above the river, with its optical axis almost perpendicular to the water surface (Figure 1b), which helps to maximize the texture details.The measurement campaign has been performed on August 9th 2017.During the campaign four customized ground control points (GCPs), clearly visible on thermal imagery, were installed.The GCPs are made of PVC with metal nails located in the middle (shown as S in Figure 1a), which makes them distinguishable in thermal imagery because of emissivity differences between different materials.The points were measured with sub-mm accuracy using a total station, which are then used to orientate the thermal image sequence.
The thermal tracer in this study is heated water.It has been heated above 60°C in the laboratory and stored into thermos flask.The tracer has been added in small portions upstream close to the measurement area.To assess the accuracy of the estimated flow velocities based on thermal imagery, reference measurements are necessary.In this study, a conventional propeller gauge has been used to achieve independent velocity values (Figure 2).The device is a propeller that is positioned as close to the surface as possible in flow direction.Thereby, rotations of the propeller are counted and afterwards converted into flow velocity using a calibration function.Measurements have been performed at different positions within the area of interest.Afterwards, the accurate location of the propeller gauge has been determined using a Structure-from-Motion (SfM) tool processing several images taken by a RGB camera during the measurement from varying perspectives.These referenced velocity measurements are subsequently referred to as check points (M) (Figure 1a).As an invasive method, the propeller measurements would interfere with the image-based measurement; therefore reference measurements were not taken simultaneously, and a stationary flow has to be assumed in order to compare the measurements.

CAMERA CALIBRATION
In order to be able to derive accurate measurements from thermal camera image data, the camera has to be calibrated both concerning geometric and radiometric aspects.Geometric calibration has to be performed to establish a precise geometric relation between pixel coordinates and corresponding object points.The goal of radiometric calibration in the application is mainly to remove the non-uniformity in thermal image sequences, while the aspect of converting each pixel value to an actual temperature value is of less importance here.

Radiometric Calibration
In terms of radiometric calibration, two concerns must be taken into consideration: Spatial non-uniformity due to fabrication variations and temporal non-uniformity derived by changing sensor temperature.On the one hand, a vignetting effect which belongs to spatial non-uniformity severely degrades the image quality.On the other hand, temporal non-uniformity, which stems from the fact that the camera output depends not only on the object irradiance but also on the time-variant sensor temperature, increases the measurement instability.Therefore, a radiometric calibration procedure (Lin et al., 2017a) is carried out to address these two non-uniformity problems.
Firstly, a shutter-less temporal non-uniformity correction (NUC) (Lin et al. 2017b), which takes advantage of polynomial curve fit models and real-time sensor temperatures, is applied to stabilize the unstable outputs.Next, multi-point correction (Tempelhahn et al. 2016) is implemented to get rid of vignetting effect and to improve the contrast in the thermal imagery.Finally, radiant temperature could be retrieved using a Planck curve after temporal and spatial NUC.However, kinetic temperature, which refers to internal or true temperature, would still be unknown unless the emissivity of the object is known (Vollmer and Möllmann, 2010).

Geometric Calibration
In order to ensure high geometric accuracy of the measurement, the thermal camera (FLIR Ax65) was calibrated before the flow velocity measurement using self-calibrating bundle adjustment supported by a RGB camera (Canon 1200D).The geometric calibration was performed with a 3D calibration field which was designed to permit the circular targets to be clearly visible on both RGB and thermal imagery.The thermal targets consist of circles made of silver foil placed on black velour foil (Westfeld et al. 2015).The calibration field placed outdoor warrants for good contrast, as the silver circular targets made of heat protective foil reflect most of the global irradiance and emit little due to their good specular reflector characteristic, while the black velour foil with excellent absorbing property emits almost all of the absorbed irradiance (Figure 3).
Coded targets and additional calibrated reference bars (which are only sufficiently detectable on the RGB images) were used to find the object coordinates of the uncoded thermal targets.Ellipse fitting is used to measure the center image coordinates of the targets.Afterwards, the interior orientation parameters (including focal length, principle point, lens distortion parameters) of the thermal camera could be determined by selfcalibrating bundle adjustment.

ACTUAL DATA PROCESSING
The core of actual thermal image sequence processing is tracking features through subsequent images.In order to translate the image measurements from pixel-space to metric surface-related measurements, the exterior orientation of the camera above the water surface has to be determined.This is realized by a spatial resection on the basis of ground control points (GCPs).We developed an automatic detection procedure of GCPs and tracking points on thermal image sequences.The main role is to detect new tracking points every several images automatically.Besides the use of orientating the thermal image sequence, another benefit for automatic detection of GCPs in each image is to detect possible camera motion and correct the velocity vectors for this effect.

Image Segmentation
The primary role of performing image segmentation is to extract interest regions (control point regions and hot water regions) matched objects, which helps to automate detection of GCPs and tracking points.Multi-resolution segmentation (MRS) has proven its ability to capture meaningful image objects, which maximize internal homogeneity while preserving external discontinuity (Benz et al. 2004).Thus, an unsupervised classification approach, based on MRS, is implemented on the corrected images after temporal and spatial NUC.

Interest Regions Extraction
Interest regions extraction (including control point regions and hot water regions) is implemented using the workflow shown in Figure 5. On the basis of image segmentation results, candidate objects, which have higher gray values than the mean value of all the objects, are extracted firstly because the control point regions and hot water regions are recorded brighter than the background.Regarding the control point regions, a border index (BI) feature is used to extract them considering that they are equipped with approximate rectangular shapes.BI is calculated as the ratio between the border length of an image object and the perimeter of its smallest enclosing rectangle, which is expressed in equation ( 1).
(1) Where = border length of image object v = length of the smallest enclosing rectangle = width of the smallest enclosing rectangle BI describes the smoothness and rectangular similarity of an image object.The more similar an image object is to a rectangle, the lower is its BI.Therefore, each image object whose BI is less than the threshold a is classified as a control point region.
The threshold can be determined by quantile statistical analysis shown in equation ( 2), which extracts a specified percentage p of the candidate objects with relatively low BI feature.
(2) Where = a-prior number of control point regions = number of candidate image objects Afterwards, the remaining candidate image objects are classified as hot water regions only if their standard deviation (SD) values are higher than the mean SD value of all the candidate objects, because only regions with rich texture contain well-trackable features.

GCPs and Tracking Features Detection
GCPs show up as strong corner features and are located at the approximate geometric center of the homogeneous control point regions.Therefore, Shi-Tomasi algorithm (Shi and Tomasi, 1994) is implemented in 10 10 pixels search windows within the control point regions to find the strongest corner features as initial positions.Then, these positions are iteratively refined to sub-pixel accuracy using a gradient-based search in OpenCV (Bradski et al. 2005).
Local hot points are taken as interest features to be tracked because of their texture with high thermal gradients.Pixels in a hot water region are determined as interest points when they fulfill the following two conditions: their gray values are local maximum values in a 10 10 pixel spatial window, and their gray values are higher than the average value of hot water region.

Tracking in Thermal Image Sequences
The actual tracking procedure includes four main stages: Adding New Tracking Features: The thermal river surface texture is changing continuously due to mixing of inserted hot water and river water.Therefore, it is necessary to detect new trackable features when previous features lose their contrast.In this study, we added new detected features after every 10 th image.
Tracking Method: The Lucas-Kanade tracker (Equation ( 3)) is used here for tracking detected features between subsequent images.The tracker uses a 6 parameter affine transformation to iteratively minimize the sum of the squares of the gray value differences between candidate patches in every two consecutive images.In order to handle large motions, a pyramidal Lucas-Kanade implementation (Bouguet, 2001) is utilized here.Only the 2 shift parameters of the affine transformation are used for velocity vector determination, the other parameters only serve for patch adaption.Deleting Bad Tracking Features: The accuracy of the tracked features tends to decrease over time because of decreasing contrast between them and their surroundings.A trackingquality-assessment feature, described by minimum eigenvalue of spatial gradient matrix (Tommasini et al., 1998), is calculated as a criterion to remove poor tracking results caused by low contrast.Any tracking point whose feature is lower than a threshold b is eliminated.

Flow Velocity Estimation:
The exterior orientation parameters of the pre-calibrated camera are obtained by a spatial resection using the four GCPs.Final flow velocity measurements (m/sec) are calculated using collinear equation and the assumption of a planar water surface.The orientation parameters are also used to project check points (M) onto the image sequence.Thus, image-based velocity values (pixel/sec) at each check point are calculated by inverse distance weight interpolation of candidate tracked points within a 5 5 pixel spatial window around its location.

RESULTS AND DISCUSSIONS
A pilot study was conducted at a small river to show the performance of the camera calibration (radiometric calibration, geometric calibration) procedure as well as to validate results of actual data processing (interest regions/points extraction, features tracking).

Calibration Performance
Results of the radiometric calibration are shown in Figure 6.As shown in (b), vignetting effect in the original images (a) is removed after temporal and spatial NUC.Corresponding radiant temperature images are shown in (c).The kinetic temperature of inserted hot water can be regarded as similar to its radiant temperature considering that the emissivity of water is close to 1.0 and that reflected radiation under the condition of nadir view is negligible.Thereby, a temperature decrease to around 24 degree is indicated when the hot water passes through the camera measurement area, whereas the original river water reveals a kinetic temperature of about 18 degree.For other objects such as rock and vegetation, the kinetic temperature is unknown because of undiscovered emissivity values.
On the other hand, the interior orientation parameters, calculated by geometric calibration allow for the calculation of undistorted images for accurate feature tracking.Although a low impact of distortion is given for the images (e.g. in terms of the distortion effect for all the check points, maximum change is 1.06 pixel for M4), geometric calibration is still recommended to ensure highest position accuracy, because an inaccurate determination of focal length and principle point may lead to larger position error.

Tracking Performance
Examples of reconstructed surface flow velocity fields are displayed in Figure 9.As shown in (a) and (b), most of the tracking results are reliable due to high contrast at the beginning stage of inserting hot water into the river.However, increasingly more unreliable results, which are shown as 'fixed points without movement' in (c) and (d), occur because of decreasing contrast over time.If no criterion is implemented to remove these inaccurate results, outliers may occur in the velocity fields.Therefore, threshold b plays an important role in removing unqualified tracking features.
Reference results provided by propeller gauge and image-based results for different threshold values b are displayed in Figure 10.Furthermore, statistical analysis of comparative results is shown in Table 1.
Taking the results achieved by the propeller gauge as a reference, in general, image-based results are in agreement with reference results measured with the propeller gauge.As shown in Figure 10 and Table 1, in terms of the results provided by 'Measurements under automatic selection of b', most of the deviations between both approaches (except for M3) are below 0.03 m/sec.M3 reveals the largest differences (more than 0.08 m/sec).
In addition, for most of the check points, the accuracy of imagebased measurement increases when rising the threshold b.However, too high thresholds would lead to a lack of results.Thus, an automatic determination method of b is to select the possible maximum value which ensures that all of the check points are measurable (e.g.0.004 shown in Figure 10).Furthermore, maximum cut-off b values of different check points largely depend on the surrounding tracking density.For example, in terms of check points M2, M4 and M8 which are located in a dense tracking environment, their cut-off b values are much higher than others with sparse particles (M3, M7).
Highest accuracies in regard to comparing thermal image tracking results and propeller velocities are achieved at M2 and M4 (Figure 10 and Table 1).The main reason for this is that these 2 CPs are located very near to the position where hot water was inserted and thus very rich textures pass through these two areas.
Deviations between propeller and image based velocity measurements are highest at M3 (Figure 10 and Table.1).Although M3 is also positioned next to the water-added spot, it is believed that something hinders the textures from passing through this area.Thus only very sparse textures pass through the area where M3 is positioned, which leads to few candidate velocity measurements available and clearly decreases the reliability of the measurement.
Another possible explanation for larger deviations here is the location of M3, which is surrounded by lower water level areas due to underlying stones.This might create diverse velocity fields in different vertical depths due to turbulences over varying river bed roughness.Thus, different flow velocities are measured because the propeller gauge has to be placed under the water while the hot water only floats at the water surface.However, this is a common drawback for image-based methods because no vertical velocity profiles can be acquired and thus limiting the transferability of the propeller gauge measurements to the image based results.In general, comparison between image-based results and propeller gauge results should be interpreted with care.

CONCLUSIONS
In this paper, a method for automatic flow velocity measurement using thermal image sequences is introduced.In order to avoid the problem of particle clustering, heated water is preferred over other thermal tracers, such as ice dices.The main advantage is that this method is not affected by natural illumination conditions and could provide dense tracking results automatically.While the main disadvantage is its limitation to close-range applications with small rivers.The introduced approach provides a field suitable method to measure flow velocities with minimal labour and material cost.

Figure 1 .
Figure 1.Data acquisition setup: (a) area of interest: GCPs are numbered as S and positions of propeller gauge measurements are marked as M. (b) Thermal camera on tripod above the small river.

Figure 2 :
Figure 2: Flow velocity measurement with a propeller gauge: Small image illustrates the propeller gauge from Dyck & Peschke (1995).
3D Geometric Calibration Field: (a) RGB Image (b) Thermal Image

Figure 5 .
Figure 4. 1st Image Segmentation Result: (a) Best SP Selection (b) Segmentation Result under Best SP MRS uses a scale parameter (SP) to control the maximum allowed heterogeneity and the size of image objects while SP is calculated by a weighted function of spectral and shape information of objects.However, the automatic selection of the user-defined SP in MRS is still a big problem.Drǎguţ et al.(2010) proposed a hypothesis as follows: When the size of an object grows, its local variance (LV) value increases continuously up to the point that it matches the object in the real world.Thus an ascendant trend is expected when plotting LV against SP graph, with break points shown as optimal SP value.Under this hypothesis, an optimal SP estimation method, which searches for the peaks in the rate of change of LV (ROC-LV) graph inside the initial SP range as the most appropriate segmented manner, is used to find the best SP for each image.Furthermore, considering that SP value mainly indicates the spectral difference between the interest regions and the background, the initial SP range can be obtained by roughly estimating the gray value subtraction between hot water and river water as well as the subtraction between GCPs and river water in the images.All the images, which are used to acquire new tracked features, are segmented with the above workflow.For example, for the first image, LV and ROC-LV values change via SP and segmentation results under best SP are shown in Figure4.Note that the initial SP range is roughly determined by simple manual estimation in the first image (60-90), which is then applied to all of the sequential images.The delineated relevant objects provide the basic results for further interest regions extraction.
Figure 7. Interest Regions Extraction Results of 61 st Image: blue regions refer to control point area, while pink region represents hot water area.

Figure 8 .
Figure 8. Detected GCPs and tracking features on 61 st image: blue points refer to GCPs while red points represent tracking features

Table 1 .
Statistical Analysis of Measurement Results