Ice Velocity in Upstream of Heilongjiang Based on UAV Low-Altitude Remote Sensing and the SIFT Algorithm

: In river management, it is important to obtain ice velocity quickly and accurately during ice flood periods. However, traditional ice velocity monitoring methods require buoys, which are costly and inefficient to distribute. It was found that UAV remote sensing images combined with machine vision technology yielded obvious practical advantages in ice velocity monitoring. Current research has mainly monitored sea ice velocity through GPS or satellite remote sensing technology, with few reports available on river ice velocity monitoring. Moreover, traditional river ice velocity monitoring methods are subjective. To solve the problems of existing time-consuming and inaccu-rate ice velocity monitoring methods, a new ice velocity extraction method based on UAV remote sensing technology is proposed in this article. In this study, the Mohe River section in Heilongjiang Province was chosen as the research area. High-resolution orthoimages were obtained with a UAV during the ice flood period, and feature points in drift ice images were then extracted with the scale-invariant feature transform (SIFT) algorithm. Moreover, the extracted feature points were matched with the brute force (BF) algorithm. According to optimization results obtained with the random sample consensus (RANSAC) algorithm, the motion trajectories of these feature points were tracked, and an ice displacement rate field was finally established. The results indicated that the average ice velocities in the research area reached 2.00 and 0.74 m/s, and the maximum ice velocities on the right side of the river center were 2.65 and 1.04 m/s at 16:00 on 25 April 2021 and 8:00 on 26 April 2021, respectively. The ice velocity decreased from the river center toward the river banks. The proposed ice velocity monitoring technique and reported data in this study could provide an effective reference for the prediction of ice flood disasters.


Introduction
Rivers north of 30° N in China generally develop ice hazards at the stages of winter river closure and spring river opening [1]. With increasing ice concentration and ice velocity in a given river section, the impact on hydraulic structures in the river increases. Upon river thawing in spring, specific river sections can become blocked by large amounts of ice, such as narrow river courses, shoals, sharp and continuous bends, and the front edges of unfrozen ice sheets. Ice accumulates to form ice jams or ice dams, which often severely block water flow sections and notably raise upstream water levels, resulting in ice jam flooding. Under long-term severe ice conditions, the impact of ice drift in spring can cause varying degrees of damage to hydraulic structures such as piers and bank revetments, threaten the lives and property of people, and further increase the difficulty of ice flood prevention [2]. Therefore, feasible and accurate ice velocity monitoring technology constitutes a major problem urgently requiring a solution in the field of water disasters. Ice velocity monitoring is not only the main priority in the fields of major natural disaster monitoring and mitigation, but is also one of the major problems related to water engineering safety, the national economy, livelihoods of people, and social development [3].
At present, some achievements had been made regarding the velocity measurements of sea ice and glaciers. Sea ice velocity monitoring methods are mainly divided into two categories: field observation methods and remote sensing monitoring methods [4]. Of these categories, field observation methods mainly include marker monitoring methods and global positioning system (GPS)-based monitoring methods [5,6], whereas remote sensing monitoring methods are largely categorized into optical and microwave remote sensing monitoring methods [7,8]. The methods applied in river ice velocity monitoring are mostly field observation methods [9], for example, visual tracking methods [10]. These methods rely on a simple principle but require highly experienced observers. Moreover, field observation methods exhibit notable subjectivity and high time and labor costs, and these methods are unsuitable for large-scale and long-term tracking. In the GPS monitoring method [11], large ice blocks with a good shape are selected, GPS sensors are installed atop these ice blocks, and the instantaneous speed of the ice blocks equipped with GPS sensors is monitored in real time with a high accuracy through the positioning mode involving a satellite, GPS sensors, and base station. For example, Song Chunshan et al. [12] measured the ice drift velocity in the upper reaches of the Heilongjiang River with the GPS monitoring method in 2019 and studied the drift ice blocking phenomenon in each reach during the ice flood period based on the velocity variation. However, the GPS tracking method can only obtain the velocity of a single ice block but cannot obtain the ice velocity field in large-scale areas, and the observation process (sensor placement and recovery) results in a high time cost and difficulty. Deng Xiao et al. [13] set up a camera on the shore, manually selected two images at a certain time interval in the obtained monitoring footage, and determined the moving distance of drift ice to calculate the ice drift velocity. However, due to the influence of side view, images typically suffer the problem of image geometric distortion, which affects the monitoring accuracy. In addition, it takes a long time to process data manually.
With the continuous improvement of unmanned aerial vehicles (UAVs) and sensor technology, small UAVs have been widely deployed in low-altitude remote sensing measurements [14][15][16][17]. Compared to the visual and GPS-based tracking methods, the multisensor UAV platform exhibits the characteristics of high efficiency, low cost, flexible operation, and strong suitability in complex river environments. This technique not only overcomes the difficulty of image data acquisition but also avoids the image distortion phenomenon attributed to side shooting from the shore, which improves the accuracy of the original data and provides a new technical means for ice velocity monitoring within a large range.
To solve the above problems, this paper proposes the combination of the SIFT and RANSAC algorithms to track feature points and researches the velocity vector of drift ice on this basis. To verify the feasibility of the proposed method, the Mohe section located in the upper reaches of Heilongjiang Province is selected as the observation area. An UAV is employed to obtain high-resolution video data during the ice flood period. Images are processed in OpenCV software that are automatically extracted at 1-s intervals in video. Subsequently, the ice plane velocity field in the observation area is extracted to perform ice dynamic analysis, thus providing a suitable technique and reference data for ice flood disaster early warnings and prediction.

Overview of the Research Area
The Heilongjiang River is located in the northernmost part of China, with a wide drainage area and abundant water resources. The drainage basin is located in the cold temperate zone. The ice thickness can reach 2 m in winter, and the ice surface is covered with snow exhibiting an average thickness of 30 cm [18]. Due to the unique geographical location of Heilongjiang, the upper reaches of the Heilongjiang River flow from low latitudes to high latitudes [19] through mountainous areas with narrow and tortuous river courses. There is abundant rainfall and runoff but insufficient observation equipment and measurement means. During the ice flood period, ice mainly accumulates in the Mohe section to form ice dams [20]. Therefore, the research area is the Mohe section in Heilongjiang Province, and the takeoff and landing point of the UAV is an open area southeast of the Mohe hydrological station (as shown in Figure 1).

Research Method
At 16:00 on 25 April 2021 and 8:00 on 26 April 2021, a DJI quad-rotor small UAV (Mavic 2 Pro) was employed for fixed-point shooting. The UAV is equipped with a Hasselblad L1d-20c (CMOS) aerial camera with an effective resolution of 20 million pixels and can record 4K high-dynamic range (HDR) footage. The weather conditions captured in each aerial photograph were sunny, the wind speed was lower than 10.7 m/s, the flight altitude was set to 400 m, the elevation difference between the take-off point and water surface was 5 m, the lens was oriented vertically downward, and the UAV hovered during shooting.
In this study, the tracking process of ice feature points included four main steps. Step 1. Feature point selection. First, a large number of ice feature points were extracted from the phase I image (a certain frame of image in the observation video) with the SIFT operator, and ice feature points were extracted from the phase II image (the frame of image after 1 s in the observation video) in the same way.
Step 2. Feature point matching. Based on the Euclidean distance, a BF matcher was adopted to track the same feature points in the above two images.
Step 3. Match result filtering. The obtained matching feature point pairs were filtered with the RANSAC algorithm to remove any incorrect matching feature point pairs.
Step 4. Ice velocity determination. The relationship between the UAV height and the actual length represented by a unit pixel was determined through calibration, the distance of the filtered feature point pairs was transformed into the moving distance of drift ice, and, finally, the ice drift velocity was determined.

Feature Point Selection
The SIFT algorithm is a feature point detection algorithm based on the scale space first proposed by Lowe [21] in 1999 and formally put forward after summary and improvement in 2004 [22]. A difference of Gaussian (DOG) scale space is constructed. The scale space describes the original image at different levels, and the different levels are represented by various scale parameters. Koenderink [23] and Lindeberg [24] demonstrated that the Gaussian kernel is the only linear kernel for scale space construction.
The Gaussian function of scale transformation is defined as: where is a variable-scale Gaussian operator,  is the convolution operation, and σ is the mean square deviation of the Gaussian normal distribution, which is denoted as the spatial scale factor. The larger the σ value, the more blurred the image becomes. Conversely, the smaller the σ value, the clearer the image becomes. To ensure that the extreme points in the scale space attain a suitable uniqueness and stability, Lowe [22] proposed the DOG concept, which is the difference between two adjacent scale images. The calculation equation of the differential Gaussian image is as follows: where k is a constant factor of the difference between adjacent scale images, and the k value is related to the number of pictures (S) in each layer of the scale space, i.e., = 2 ⁄ . After construction of the DOG pyramid, it is necessary to select key points from the pyramid scale space, namely, local feature extreme points in the image space. Each sampling point is compared to 8 adjacent points in the current scale layer and 18 points in the upper and lower adjacent scale layers to determine whether the considered sampling point is an extreme point in the current scale space and two-dimensional space [25].
Because the extreme points detected with the above method are distributed in the discrete space, the extreme points in the known discrete space are interpolated to extreme points in the continuous space through subpixel interpolation. To obtain the exact position of an extreme point, the DOG function is expanded according to the Taylor formula as follows: where D is the D (x, y, σ) value at the feature point position and X = (x, y, σ) T is the offset relative to this feature point. The extreme value of D(x) is determined, i.e., = 0, and the position of the extreme point is: The extreme value at the corresponding extreme point is: .03, this feature point is removed as an unstable feature point. Because the DOG operator can produce strong edge responses, it is necessary to apply the Hessian matrix at the feature points to eliminate points with an unstable edge response [26]. The Hessian matrix is defined as: where Dxx, Dxy, and Dyy denote derivations of x or y. Among the accurately located stable feature points, the parameter direction is specified for each feature point based on the gradient direction distribution characteristics of its neighborhood pixels, and the feature vector is then generated relative to this direction so that the SIFT-determined feature points remain invariant to image rotation.
Gradient modulus: Gradient direction: Through the above process, the feature points contain three basic components of information: position, scale, and direction. Moreover, the feature points are invariant to translation, scaling, and rotation. To ensure that the feature points remain unchanged under various conditions, such as solar radiation level and viewing angle, a descriptor is established for each feature point. The feature point is adopted as the center, the coordinate axis of the image is rotated to match the direction along which the feature point is located, the neighborhood near the feature point is divided into 4×4 subwindows, the gradient information in each subwindow is obtained, an 8-dimensional direction histogram is generated, a 128-dimensional feature descriptor is finally established, and the descriptor subvector elements are normalized to improve the resistance to light changes [27].

Feature Point Matching
The BF matcher is a simple matching algorithm. Suppose the set of all feature points in image P1 is {a1, a2, …an}, and the set of all feature points in image P2 is expressed as {b1, b2, …bn}. When the BF matcher is applied to match the feature points in images P1 and P2, the Euclidean distance from the feature vectors of each feature point in image P1 to each feature point in image P2 is calculated. The result of the algorithm is an n × m-dimensional , where dij is the Euclidean distance between the feature vectors of feature points ai and bj [28].

Match Result Filtering
After all feature points are matched, incorrect matching points inevitably occur, which can affect the accuracy of ice velocity calculation. In this paper, the RANSAC algorithm is implemented to further optimize the matching results, and the matching feature points are screened through continuous iteration.
The RANSAC algorithm is an iterative algorithm proposed by Fischler and Bolles [29] to estimate the parameters of mathematical models. This algorithm has been widely applied in line fitting and plane fitting. Correct data contained in the considered data are recorded as inner points, and abnormal data (or noise) are recorded as outer points [30]. The probability of obtaining the correct model through the RANSAC algorithm is: where P is the probability of the inner points, K is the minimum number of data points required to solve the model, and M is the number of iterations. The above equation yields the following: Finally, optimal matching feature point pairs are obtained, i.e., the feature point pairs with the same name between the input phase I ice image and matching phase II ice image.

Ice Velocity Calculation
After all feature point pairs with the same feature descriptor are obtained, the horizontal displacement of ice can be extracted, and the feature point pairs with the same feature descriptor can then be transformed into the ice velocity according to the pixel scale and image interval, as follows: where X1, X2, Y1, and Y2 are the coordinates of the feature point pairs with the same feature descriptor, S is the moving distance of the feature point in the image (px), R is the pixel scale (cm/px), T is the time between the two images (s), and V is the ice velocity (m/s).
To realize quantitative analysis of the ice velocity, first, the UAV is deployed to shoot A1-sized paper panels at different heights to determine the shooting height of UAV remote sensing images and the actual area represented by a unit pixel. After the actual area is squared, a fitting equation between the unit pixel length and actual length in remote sensing images at the different heights is obtained, as shown in Figure 2.

Results
Through feature point selection, matching, and filtration of the two obtained orthophoto images, 1429 pairs of effective feature points (as shown in Figure 3) were extracted from the river section at 16:00 on 25 April 2021, and 2726 pairs of effective feature points (as shown in Figure 4) were extracted from the river section at 8:00 on 26 April 2021. According to the principle of feature point selection, there were sparse or even absent feature extreme points on the river surface, but the distribution of the feature points on the ice surface was relatively dense and stable, which satisfied the basic requirements of ice velocity analysis. At 8:00 on 26 April 2021, due to the large amount of broken ice in the river section, the selected effective feature point density was higher than that at 16:00 on 25 April 2021.   Figures 5 and 6, respectively. The distribution characteristics of the ice velocity are similar. Affected by the river bend, areas with the highest ice velocity were concentrated on the right side of the river center, and the ice velocity tended to decrease toward the banks on both sides. Moreover, there occurred a sudden change in the ice flow velocity near the river banks on both sides, which was likely caused by collision between drift ice and the bank slopes. Figure 7a     To study the distribution characteristics of the ice velocity during the ice flood period in the Mohe section of Heilongjiang more comprehensively, the transverse and longitudinal components of the ice velocity during the different periods were calculated by tracking the displacement components of the effective characteristic points in the corresponding coordinate system. For this calculation, Vx is the transverse component of ice velocity, Vy is the longitudinal component of ice velocity, and the positive and negative symbols indicate the direction of the velocity component. According to the schematic diagram of the observation geographical location, the observation area is located at a river bend. As indicated by the transverse velocity component shown in Figure 9a,c, the ice velocity upstream of the river was biased to the right side of the river, which can easily damage the right bank of the river, form ice accumulation, seriously block the water crossing section, and threaten the lives and property of people. For example, the ice accumulation and damage to the bank slope are shown in Figure 10, which are in line with the actual site conditions. As indicated by the longitudinal velocity component depicted in Figure 9b,d, the ice velocity was oriented along the river channel, and there occurred no ice reflux phenomenon affected by vortexing.  To highlight the significance of this research, the related research results are compared. Deng Xiao [13] carried out a series of studies on the ice velocity in the Heilongjiang River. Under the assumption that the ice flows horizontally, the ice velocity was confirmed by manual processing. It needs to be noted that the average velocity of ice was 1.62 m/s and the image data were recorded at 10:28 on 2 May 2016. The ice movement at different times are shown in Figure 11. Although the ice velocity can be determined, this method can only manually acquire the velocity for a single ice block, and this method typically suffers the problem of image geometric distortion, which affects the monitoring accuracy and is time consuming. In addition, Chave [31] installed a real time data acquisition system in Lake St. Pierre in the St. Lawrence River to measure the ice velocity. The system can work around the clock. Unfortunately, the instruments often work underwater and cannot be moved. Therefore, the cost of installation and maintenance is not economical. However, the proposed method, UAV combined with SIFT and RANSAC algorithms, can fill these gaps. The velocity of ice blocks can be determined simultaneously, and a large number of in situ monitoring instruments are not required. In addition, the proposed method is suitable for ice running at high speed.

Ice Velocity Verification
Upon completion of feature point matching in ice images obtained with the UAV, although the matching results are further screened with the RANSAC algorithm, ice may collide with other ice blocks or bank slopes, thus changing the local extreme points, which can lead to incorrect matching or the absence of feature point pairs in images. To quantitatively evaluate the accuracy of ice velocity determination in the research area, ice was selected in the river section, feature points with the same name were extracted through visual interpretation (as shown in Figure 12), the ice velocity was manually calculated, and the obtained calculation value was compared to the ice velocity within 1 m before and after the section based on the SIFT algorithm, as shown in Figure 13a,b, respectively. The errors were less than 0.02 m/s. It was determined that the ice velocity determined based on the SIFT algorithm was consistent with the actual situation.

Conclusions
To monitor and realize in-depth analysis of ice change information, this paper tracks ice surface feature points based on UAV remote sensing images and the SIFT algorithm, establishes the corresponding ice velocity field, and finally carries out in-depth analysis and interpretation of ice movement based on the established ice velocity field. The following conclusions are obtained: 1. Feature points were selected as local extreme points. Only ice surfaces, edges, and corners can produce dense and stable feature points. The SIFT algorithm provided the advantages of high precision and notable robustness. The SIFT and RANSAC algorithms were jointly employed to track the feature points, which realized ice velocity monitoring. 2. At 16:00 on 25 April 2021, the ice velocity was mostly in the range 2.40~2.50 m/s, accounting for 30.02% of all velocity values. The maximum, minimum, and average ice velocities were2.65 m/s, 1.11 m/s, and 2.00 m/s, respectively. At 8:00 on 26 April 2021, the ice velocity was mainly in the range 0.90~1.00 m/s, accounting for 50.33% of all velocity values. The maximum, minimum, and average ice velocities were1.04 m/s, 0.38 m/s, and 0.74 m/s, respectively. The areas with the highest ice velocity in the studied reach were concentrated on the right side of the river center, and the ice velocity tended to decrease toward the banks on both sides. In addition, based on the ice velocity field, ice accumulation formation on the bank slope was highly related to the ice velocity. 3. Under the influence of illumination changes and ice collisions, the local extreme points varied. The matching results were further screened with the RANSAC algorithm. The error between the ice velocity obtained with the SIFT algorithm and the manually calculated ice velocity was less than 0.02 m/s, which was in line with the actual situation. 4. Compared to traditional ice velocity monitoring methods, the proposed method attained a high precision, and greatly reduced time and labor costs. The research method in this paper provides a new technical means for large-scale ice change information monitoring. In the future, we can perform further research on ice jams and ice dams based on this technique.