Two-Dimensional Flow Field Measurement Method for Sediment-Laden Flow Based on Optical Flow Algorithm

Abstract

This paper proposes a novel particle image velocimetry (PIV) technique to generate an instantaneous two-dimensional velocity field for sediment-laden fluid based on the optical flow algorithm of ultrasound imaging. In this paper, an ultrasonic PIV (UIV) system is constructed by integrating a medical ultrasound instrument and an ultrasonic particle image velocimetry algorithm. The medical ultrasound instrument with a phased sensor array is used to acquire acoustic echo signals of particles in water and generate two-dimensional underwater ultrasound images. Based on the optical flow field, the instantaneous velocity of the particles in water corresponding to the pixels in the ultrasonic particle images is derived from the grayscale change between adjacent images under the L-K local constraint based on the optical flow field, and finally, the two-dimensional flow field is obtained. Through multiple sets of experiments, the proposed algorithm is verified. The experimental results are compared with those of the conventional cross-correlation algorithms. The results show that the L-K optical flow method can not only obtain the underwater velocity field accurately, but also has the advantages of good smoothness and extensive suitability, especially for the flow field measurement in sediment-laden fluid than conventional algorithms.

1. Introduction

The measurement of the depth-averaged velocity in sediment-laden flow is a prerequisite for many research works on river engineering, such as river-bank erosion and sediment transportation. Two-dimensional quantitative information of the depth-averaged velocity field helps investigate dynamic flow structures in complicated turbulent phenomena. The information enables us to compute vorticity, deformation, and other quantities, which are directly related to the dynamics of coherent flow structures [1]. However, traditional measuring instruments, such as hot-wires and laser-Doppler anemometers, are of one-point measurements, not capable of illustrating an instantaneous spatial flow structure. The conventional Doppler algorithm has a problem on the angulation error, which confines the measurement only to the velocity component along ultrasonic beams. Parallel alignment of the ultrasonic beam to the flow direction is required, which is not always feasible in many applications. It is challenging research to scientists and hydraulic engineers for a long time [2].

Optical PIV is now widely adopted as a reliable method that can obtain quantitative information on a spatial flow structure [3]. However, light cannot penetrate through turbid water because it rapidly attenuates when travelling in turbid flow, which affects its application in sediment-laden flow.

To overcome those disadvantages and meet the requirements of high-precision measurements, ultrasonic imaging techniques have been introduced into the field of PIV for their distinct advantages in terms of cost and non-transparent fluid penetration. For example, the ultrasonic PIV (UIV) technique has been used to measure the time-average velocity field in silt carrying flow under moderately concentrated flow conditions [4]. In addition, UIV can measure the instantaneous flow velocity field of single-phase turbulent pipe flow [5]. In this paper, an implementation based on ultrasound imaging is developed. A UIV system is built by integrating a medical ultrasound instrument and a measurement algorithm of particle velocimetry with ultrasound images. The combination of ultrasound imaging technique and a digital PIV method can generate a new measurement technique of two-dimensional, two-component, depth-averaged velocity field appropriate for sediment-laden flow conditions [6]. Compared with Doppler-based techniques, an advantage of the UIV is its ability to measure the flow velocity components in both parallel and perpendicular directions to the ultrasonic beam. In general, it has excellent performance on simplicity, accuracy, and accessibility, while it does not have the limitations of the Doppler-based algorithm [7].

Even with ultrasonic particle image velocimetry, flow velocity measurement algorithms still face a major challenge. When the concentration of sandy water is high enough, algorithms based on the PTV method of tracking individual particles can no longer meet the measurement requirements, because there is not only In-Plane Motion but also Out-of-Plane Motion in the acoustic plane of ultrasound images, which brings instability to the calculation, as shown in Figure 1. On the other hand, there are a large number of particles in turbid water. Their motion trajectories cross each other, which causes uncertainty in the analysis of particle motion and greatly reduces the tracking robustness.

PIV matching algorithms based on cross-correlation method are widely used for simplicity. For improving the accuracy and efficiency of the flow field measurement, the multigrid deformation algorithm is used to calculate the flow velocity with ultrasonic particle images of sandy water streams with volume sediment concentration less than 5‰ [4]. Cross-correlation is also applied in photoacoustic flowmetry for deep tissue measurement of blood flow velocity at physiological speeds ranging from 3 to 25 mm/s [8]. Based on cross-correlation technology, fibre elimination is realized by calculating the coordinate transformation parameters between separated cotton images, and the Kalman filter tracking is utilized to predict the estimated location to improve efficiency [9]. Although PIV matching algorithms based on the cross-correlation method are continuously improved, they usually require high image quality, good correlation between images, and long computation time [10].

In recent years, particle image velocimetry algorithms based on the optical flow method have received increasing attention from researchers [11,12]. The optical flow method, which considers the motion vector field as a whole, is suitable for the calculation and analysis of image sequences with highly concentrated particles and continuous deformation. With a dense optical flow algorithm, its resolution can reach pixel scale, which provides a useful tool for accurate flow field measurements [13].

In summary, an ultrasound particle image velocimetry algorithm is proposed to effectively measure the two-dimensional flow field distribution of highly concentrated sand-bearing water streams. The 2D flow field measurement method combined with the UIV measurement system can significantly improve the environmental tolerance and accuracy of flow field measurement, and reduce the computational consumption by pyramidal hierarchical algorithm. Experiments have been conducted to verify the performance of the algorithm proposed in the paper. The results demonstrate that UIV is a promising technique to generate instantaneous two-dimensional fields of flow velocities in sediment-laden fluids [14]. The quantitative performance assessment of the algorithm has been made using ultrasound images acquired in seeded flow. The algorithm accurately measured the 2D flow field in turbid water with 10‰ volume concentration and outperformed the cross correlation method.

2. Ultrasound Image Acquisition

B-mode ultrasound imaging system, which adopts the ultrasound pulse-echo technique, is widely used. It transmits ultrasound waves $s\left(t\right)$ into the medium with ultrasound probes. The waves are reflected, refracted and scattered when meeting objects, and then echo signals $r\left(t\right)$ with location information are received by the ultrasound probes. Signals are processed with AD, time gain compensation (TGC), then beamforming is used to focus the multi-element signals. Finally, ultrasound images representing the distribution of objects inside the medium are obtained, which can be stored and displayed after image processing. The principle is shown in Figure 2.

As shown in Figure 3, a UIV measurement system is proposed in the paper. For this system, an acoustic array are motivated to emit acoustic beams. Those beams point to a specific direction at every launch time [15]. After travelling in the water, echoes reflected from tracer particles or suspended sediment are received by the acoustic array. The intensity and flight duration of the acoustic echoes are recorded in the system. Echo intensity reflects the existence and material of the particles in the fluid, and flight duration reflects target distances from the array. Based on that information, the fluid morphology on a measuring line on depth-direction can be achieved. Repeat above operations, and a two-dimensional fan-shaped ultrasound image is generated.

According to the motion characteristics of high concentration sand-bearing flow, the main influences on the motion of sediment particles are the resistance caused by relative motions between particles, the upward force cased by shearing actions of the water flow, and the additional mass force caused by variable speed motion. These types of forces will reduce the correlation of the B-mode ultrasound image sequences acquired in the presence of high concentrations of sandy water flow.

Because of the high concentration of tracer particles, it will make the correlation map appear with multiple peaks, as shown in Figure 4. From the normalized cross-correlation formula, it is known that the cross-correlation algorithms are based on statistical calculation, and the velocity vectors obtained by those algorithms are only the displacement vectors with the highest confidence in probability.

3. B-Mode Ultrasound Velocimetry Based on Optical Flow

3.1. Constraint Equation for Optical Flow Algorithm

Multiple stationary two-dimensional underwater particle images are captured in a period by the measurement system. We need to analyze the motion state, $(u,v)$, of the underwater particles in the observed profile at a given moment, where $u(x,y)$ and $v(x,y)$ are the velocity components in two dimensions at pixel point $(x,y)$ in underwater images.

In the case of a light source with a constant intensity, the instantaneous change in the grayscale of an image pixel is caused by the motions of particles. The image sequence is expressed by $I(x,y,t)$, $(x,y)$ is a pixel, t is a time point. As mentioned above, the instantaneous rate of the change in two dimensions at a specific point in the underwater profile is expressed by $(u,v)$, which is referred to as the optical flow vector. Every optical flow vector corresponds to the instantaneous velocity vector at a point, and velocity vectors of all pixels in 2D underwater image constitute the optical flow field.

Our task is to determine the optical flow field by the variation and correlation of the gray intensity of the pixels in the image sequence in the time domain, so as to obtain the motion of each pixel position, then obtain the 2D optical flow field.

The optical flow calculation usually uses the brightness constancy model (BCM), of which there are two assumptions: The brightness is constant between adjacent frames; The objects’ motion between adjacent frames is relatively small.

Assuming that $I(x,y,t)$ is the gray level of pixel $(x,y)$ at the moment t, so the BCM can be expressed by the following equation:

$$I(x+dx,y+dy,t+dt)=I(x,y,t)$$

(1)

where, $dx,dy,dt$ are the increments of $x,y,t$ respectively.

Expanding $I(x+dx,y+dy,t+dt)$ with the first-order Taylor formula, we can obtain following expansion formula:

$$I(x+dx,y+dy,t+dt)=I(x,y,z)+\frac{\partial I}{\partial x}dx+\frac{\partial I}{\partial y}+\frac{\partial I}{\partial t}+e$$

(2)

Substituting Equation (2) into Equation (1), the constraint equation for optical flow algorithm is obtained as below:

$$I_{x}u+I_{y}v+I_t=0$$

(3)

where $u(x,y),v(x,y)$ is the components of optical flow at pixel $(x,y)$, $I_x,I_y,I_t$ are the partial derivatives of I with respect to $x,y,t$, respectively.

3.2. L-K Smoothness Constraint Equation

Equation (3) has two variables to be solved, so the solution is not unique. Constraints should be imposed to get the appropriate optical flow vector.

A predominant constraint method is the L-K algorithm with the consistency constraint assumption. The constraint of the L-K algorithm is that the pixel motion is consistent in the local region. Assuming that the motion vector remains constant in the spatial neighborhood of the pixel, weighted least-squares is used to calculate optical flow by assigning different weights to the residual values of different pixel points. So Equation (3) is rewritten to the following formula:

$$e(x,y)=\sum _{(x,y)\in \mathsf{\Omega }}{W}^2(x,y){(I_{x}u+I_{y}v+I_t)}^2$$

(4)

where $W(x,y)$ is the weight assigned to the residual value of a pixel point, $(x,y)$.

In order to find the optimal value of optical flow to make $e(x,y)$ minimal, let the partial derivative of $e(x,y)$ with respect to V($V=(u,v)$), be 0. According to the method, the optical flow vector, $V={(u,v)}^T$, can be derived:

$$V={\left(A^T{W}^{2}A\right)}^{-1}{A}^T{W}^{2}b$$

(5)

where, A is the gradient matrix of the image $I(x,y)$:

$$A=\left[\begin{array}{cccc}\nabla I(x_1,y_1)& \nabla I(x_1,y_2)& \dots & \nabla I(x_1,y_N)\\ \nabla I(x_2,y_1)& \nabla I(x_2,y_2)& \dots & \nabla I(x_2,y_N)\\ ⋮& ⋮& \ddots & ⋮\\ \nabla I(x_N,y_1)& \nabla I(x_N,y_2)& \dots & \nabla I(x_N,y_N)\end{array}\right]$$

(6)

W is the weight matrix:

$$W=\left[\begin{array}{cccc}w(x_1,y_1)& w(x_1,y_2)& \dots & w(x_1,y_N)\\ w(x_2,y_1)& w(x_2,y_2)& \dots & w(x_2,y_N)\\ ⋮& ⋮& \ddots & ⋮\\ w(x_N,y_1)& w(x_N,y_2)& \dots & w(x_N,y_N)\end{array}\right]$$

(7)

b is the partial derivative of $I(x,y)$ with respect to t:

$$b=-\left[\begin{array}{cccc}{I}_t(x_1,y_1)& I_t(x_1,y_2)& \dots & I_t(x_1,y_N)\\ I_t(x_2,y_1)& I_t(x_2,y_2)& \dots & I_t(x_2,y_N)\\ ⋮& ⋮& \ddots & ⋮\\ I_t(x_N,y_1)& i_t(x_N,y_2)& \dots & I_t(x_N,y_N)\end{array}\right]$$

(8)

According to Equation (5), we can calculate the optical flow vectors of the pixels in the acoustic images, from which the two dimensional flow velocity field will be inferred.

3.3. Pyramidal L-K Algorithm

The constraints of the L-K algorithm are small motion speed, constraint brightness, and consistent region, which are stronger assumptions and not easily satisfied in many applications. When particles in water are moving fast, the assumptions above do not hold, and computational errors will arise, which would gradually accumulate, making subsequent calculation more deviant and eventually failing to obtain valid and reliable optical flow values.

Because fast motion speed will generate large errors of the algorithm, it is desirable to reduce the motion speed of particles in the ultrasound images to meet the algorithm requirement [16]. An intuitive way to achieve this is to shrink the size of the images. Suppose that when an image’s size is $400\times 400$, a pixel’s flow velocity is 16 pixels per frame; then when the image is downscaled to $200\times 200$, the pixel’s flow velocity is 8 pixels per frame. If ultrasound images do not satisfy the motion speed assumption, L-K algorithm can still be applied after downscaling the images several times.

In this way, a multi-scale representation of images, image pyramid, can be used in the L-K optical flow algorithm, which generates pyramid images of the ultrasound images with different sizes, to solve the problem mentioned above layer by layer. The basic approach is to obtain the UIV image sequence, $I(x,y,t)$, and downsample it separately until it reaches the specified level. In practice, a pyramid of three or four layers is usually sufficient to meet the requirements of efficiency and accuracy. Figure 5 shows the pyramid images.

The Pyramid L-K algorithm is an iterative calculation in layers, which continuously improves the resolution of the results. The optical flow value calculated from the L layer image is used as the initial value of the optical flow of the L-1 layer. Then the exact value of the optical flow of the L-1 layer is calculated, which is used as the initial value of the optical flow value of the L-2 layer. This process continues iteratively until the exact value of the optical flow of the original image is calculated.

3.4. Implementation of the Proposed Algorithm

The proposed algorithm of 2D flow field measurement of sediment-laden flow based on Pyramid L-K algorithm has three main steps: building pyramid images, tracking pixels through pyramid images, and optical flow calculation iteration.

3.4.1. Building Pyramid Images

Let $I_0$ be the image of layer 0 of the original image I, which is the highest resolution image in the pyramid image, and the width and height of the image are defined as $W,H$. To build the pyramid images, downsample by $1/2$. Let $L=0,1,2,\cdots$ denote the number of layers of the pyramid: $I^L$ denotes the image of the Lth layer, then the image of the Lth layer, then the image of the Lth layer of the pyramid can be obtained by the following equation:

$$\begin{array}{cccc}{I}^L(x,y)\hfill & =\hfill & \hfill \frac{1}{4}& I^{L-1}(2x,2y)\hfill \\ & & \hfill +\frac{1}{8}& [I^{L-1}(2x-1,2y)+I^{L-1}(2x+1,2y)+\hfill \\ & & & I^{L-1}(2x,2y-1)+I^{L-1}(2x,2y+1)]\hfill \\ & & \hfill +\frac{1}{8}& [I^{L-1}(2x-1,2y-1)+I^{L-1}(2x+1,2y+1)+\hfill \\ & & & I^{L-1}(2x+1,2y-1)+I^{L-1}(2x-1,2y+1)]\hfill \end{array}$$

(9)

In practice, it is usually sufficient to build four layers of pyramid images to meet the application requirements.

3.4.2. Tracking Pixels through Pyramid Images

The purpose of image point tracking is to find the matching pixel point N in the image to be matched for the pixel point M in the specified reference image I, or to calculate the optical flow V of the point M, where $V={(u,v)}^T$. Assuming that the coordinates of the point M in the original image I are $(M_x,M_y)$, and $M_L=(M_x^L,M_y^L)$ is the value of the coordinates of the point M in the Lth level pyramid image, $M^L=M/2^L$, which is known from the principle of image pyramid formation [17].

Let the optical flow prediction in layer L be $g^L={[g_x^L,g_v^L]}^T$, and remaining optical flow value, optical flow residual displacement vector, be $d^L={[d_x^L,d_v^L]}^T$. The relationship between the optical flow vector in layer L and the optical flow vector prediction in layer $L-1$ can be expressed as follows:

$$g^{L-1}=2(g^L+d^L)$$

(10)

Considering the top layer, layer $L_m$, has no reliable optical flow prediction, set $g^{L_m}={[0,0]}^T$, so the optical flow of original image, layer 0, can inferred as follows:

$$\begin{array}{cc}d\hfill & =g^0+d^0\hfill \\ \hfill & =2(g^1+d^1)+d^0\hfill \\ \hfill & =\cdots \hfill \\ \hfill & ={\displaystyle \sum _{L=0}^{L_m}}{2}^L{d}^L\hfill \end{array}$$

(11)

As can be seen from Equation (11), there is an obvious advantage on the pixel tracking of pyramid images, that the residual displacement vector at every level, $d^L$, is always a very small value when calculating the optical flow vector.

3.4.3. Optical Flow Calculation Iteration

The relationship between the matching error ${ϵ}^L$ and the residual value of optical flow $d^L={[d_x^L,d_y^L]}^T$ in the L layer can be expressed as follows:

$${ϵ}^L(d_x^L,d_y^L)=\sum _{x,y}{(I(x,y)-J(x+g_x^L+d_x^L,y+g_y^L+d_y^L))}^2$$

(12)

$$d^L={[d_x^L,d_y^L]}^T=argmin\left\{{ϵ}^L\right\}$$

(13)

Let $A(x,y)=i^l(x,y)$, $B(x,y)=j^l(x,y)$, then Equation (12) can be expressed as:

$$ϵ(d_x^l,d_y^l)=\sum _{x,y}{(A(x,y)-B(x+d_x^l,y+d_y^l))}^2$$

(14)

Solving the above equation using the L-K algorithm, the first-order partial derivative of $ϵ$ with respect to $(d_x^L,d_v^L)$ is zero when $ϵ$ is the minimum, and we have:

$${\left.\frac{\partial ϵ}{\partial d^L}\right|}_{d^L=d_{opt}^L}=[0,0]$$

(15)

Substituting Equation (14) into Equation (15) yields:

$$\frac{\partial ϵ}{\partial d^L}=-2\sum _{x,y}(A(x,y)-B(x+d_x^L,y+d_y^L))\left[\frac{\partial B}{\partial x}\frac{\partial B}{\partial y}\right]$$

(16)

Expanding Equation (16) with the first-order Taylor formula yields:

$$B(x+d_x^L,y+d_y^L)=B(x,y)+\left[\frac{\partial B}{\partial x}\frac{\partial B}{\partial y}\right]·{\left[d_x^L,d_y^L\right]}^T+o(x,y)$$

(17)

Substituting Equation (17) into Equation (16) yields:

$$\frac{\partial ϵ}{\partial d^L}\approx -2\sum _{x,y}(A(x,y)-B(x,y)-\left[\frac{\partial B}{\partial x}\frac{\partial B}{\partial y}\right]·d^L)\left[\frac{\partial B}{\partial x}\frac{\partial B}{\partial y}\right]$$

(18)

The value of $A(x,y)-B(x,y)$ can be considered as the derivative of the ultrasound image at the point $(0,0)$ with respect to time, denoted as:

$$\delta I(x,y)=A(x,y)-B(x,y)$$

(19)

The gradient vector of the ultrasound image can be defined as:

$$\nabla I=\left[\begin{array}{c}{I}_x\\ I_y\end{array}\right]={\left[\frac{\partial B}{\partial x}\frac{\partial B}{\partial y}\right]}^T$$

(20)

The gradient vector $[I_x,I_y]$ of the image $A(x,y)$ can be calculated from the neighborhood of the image at the pixel M. The gradient of the image can be expressed in terms of the central difference as:

$$\begin{array}{c}\hfill I_x(x,y)=\frac{\partial A(x,y)}{\partial x}=\frac{A(x+1,y)-A(x-1,y)}{2}\\ \hfill I_y(x,y)=\frac{\partial A(x,y)}{\partial y}=\frac{A(x,y+1)-A(x,y-1)}{2}\end{array}$$

(21)

Substituting Equations (20) and (21) into Equation (18), we get:

$$\begin{array}{ccc}\frac{1}{2}\frac{\partial ϵ}{\partial d^L}\hfill & \approx \hfill & (\nabla I^T·d^L-\delta I)\nabla I^T\hfill \\ & =\hfill & {\displaystyle \sum _{x,y}}\left(\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]-\left[\begin{array}{c}\delta I·I_x\\ \delta I·I_y\end{array}\right]\right)\hfill \\ \hfill & =\hfill & G·d^L-b\hfill \end{array}$$

(22)

where,

$$G=\sum _{x,y}\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right],b=\left[\begin{array}{c}\delta I·I_x\\ \delta I·I_y\end{array}\right]$$

(23)

When the error $ϵ$ achieves its minimum value, the optimal optical flow residual displacement vector $d_{opt}^L$ of the L layer can be derived as follows:

$$d_{opt}^L=G^{-1}b$$

(24)

To further solve the optical flow residual displacement vector $d^L$ in the L layer accurately, the Newton–Raphson method can be used to improve the exact solution of the optical flow.

Let the optical flow value after $k-1$ iterations at the L layer be ${\zeta }^{k-1}={[{\zeta }_x^{k-1},{\zeta }_y^{k-1}]}^T$. Then, when the kth iteration is performed, the image after superimposing a translation of the vector ${\zeta }^{k-1}$ is:

$$B_k(x,y)=B(x+{\zeta }_x^{k-1},y+{\zeta }_y^{k-1})$$

(25)

The objective of the kth iteration is to calculate the optical flow residual displacement vector ${\lambda }^k={[{\lambda }_x^k,{\lambda }_y^k]}^T$ at the L layer, which minimizes the error function represented as follows:

$$ϵ({\lambda }_x^k,{\lambda }_y^k)=\sum _{x,y}{(A(x,y)-B(x+{\lambda }_x^k,y+{\lambda }_y^k))}^2$$

(26)

Solving the above equation, we get:

$${\lambda }^k=G^{-1}{b}_k$$

(27)

where, G remains constant in the iteration, and

$$\begin{array}{c}b={\displaystyle \sum _{x,y}}\left[\begin{array}{c}\delta I_k·I_x\\ \delta I_k·I_y\end{array}\right]\delta I_k(x,y)=A(x,y)-B_k(x,y)\hfill \end{array}$$

(28)

From the residual displacement vector, ${\lambda }^k$, calculated from Equation (26), the displacement, ${\zeta }^k$ of the ultrasound image at the kth iteration can be deduced as:

$${\zeta }^k={\zeta }^{k-1}+{\lambda }^k$$

(29)

The process is continued until the calculated residual displacement vector is less than a given threshold, or the maximum number of iteration is reached.

Suppose the convergence condition of Equation (29) is reached after k iterations, then the final optical flow of the k layer can be expressed as:

$$d^L={\zeta }^N=\sum _{k=1}^K{\lambda }^k$$

(30)

The corresponding optical flow vectors are obtained by using the above method in all layers of the pyramid images. Substituting them back into Equation (11), the total optical flow vectors can be obtain.

4. Experimental Setup

4.1. Water Tank Experimental Environment and Facilities

As shown in Figure 6, a variable-slope recirculating water tank is used for the experiments. The sink is $300\mathrm{cm}$ long, $50\mathrm{cm}$ wide, and $50\mathrm{cm}$ deep, and the inner cavity is divided into left and right parts by a vertical glass plate. A grid is equipped in the middle of the flume to rectify and reduce the fluctuation of the water surface. An inlet pipe and a drainage pipe are equipped at one end of the tank, which are connected to a pump. The high concentration sand-containing water flow is driven by the pump, and enters the tank from a reservoir through the inlet pipe. The water flow circulates and then discharges back to the reservoir through the drainage pipe. The volume concentration of sand-bearing water in experiments is 10‰, i.e., the mass concentration is about $10{\mathrm{kg}/\mathrm{m}}^3$, Reynolds number is 2200, and the average diameter of sand particles is $100\mathsf{\mu }\mathrm{m}$.

The experimental system controls the speed of the pump with a frequency converter to provide a stable flow environment for the experiment. The reservoir with a total area of about $3{\mathrm{m}}^2$ is built at the end of the tank to regulate the flow pattern and water volume for the high concentration sand-bearing water flow experiments, while a flexible rubber tube is used to connect the tank to the pump motor to reduce the impact of the pump motor vibration on the measurement.

4.2. Ultrasonic Water Flow Imaging System

In the experiments, a medical B-type ultrasound instrument, SIUI APOGEE 1200, is deployed. The ultrasound probe is fixed just touching the water surface to avoid affecting the flow field, as shown in Figure 6. Emission frequency of the probe is 5.0 MHz, the temporal resolution is 60 frames/s. In order to monitor the motion of all particles in the imaging plane, the probe is placed in the same direction as the sand-containing water flow. The ultrasound instrument is connected to a computer, which converts the echo signals collected by the probe into ultrasound images and transmits them to the ultrasound particle velocimetry software in the computer for analysis, processing, and display of flow field information.

To calculate the actual particle flow velocity, it is necessary to convert the pixel units of the ultrasonic particle image into actual distance units. We place standard metal parts of different widths at different locations in the water, identify their pixel widths through image processing, and compare them with their standard physical widths to obtain a calibrated rate relationship. According to the rate relationship, the flow velocity field obtained by the ultrasonic image velocimetry algorithm can be calibrated to the actual flow velocity field. The schematic diagram of the calibration experiment is shown in Figure 7, and the software interface of the calibration algorithm we wrote is shown in Figure 8.

4.3. Water Flow Experimental Conditions

Since the wavelength of ultrasound, ${\lambda }_u$ is in the range of 0.15∼1.5 mm, and the wavelength of visible light, ${\lambda }_o$ is in the range of 380∼780 nm, ${\lambda }_u\gg {\lambda }_o$, the spot of the particle image in ultrasound imaging is larger than it in optical imaging, as shown in Figure 9. For this reason, the particle images in ultrasound imaging are more likely to saturate. In general, as shown in Figure 9d, ultrasound images will basically reach saturation when the volume concentration of sand content in the water stream reaches 10‰, namely, mass concentration is about $10{\mathrm{kg}/\mathrm{m}}^3$.

In order to ensure the smoothness of the experiment, the measurements are conducted in a constant sand-bearing water flow with a volume concentration of 10‰, the water depth is about 50 cm, and the range of width-to-depth ratio is 1.2∼1.5, which ensures that the side walls of the flume will not produce obvious secondary turbulence. The experimental observation area is about 2 m away from the inlet of the flume, which ensures that the high concentration of sand-bearing water has been fully spread when arriving at the observation area.

The experimental steps are as below. Firstly, adjust the inclination of the sloping water tank and the speed of the pump motor, so that the water depth meets the measurement requirements. Then adjust the direction of the ultrasonic probe with a laser collimator to be the same as the flow direction of the sand-containing water. Add appropriate amount of sand particles as tracer particles, and wait about 10 min so the various conditions of the whole measurement system are stable. At this time, ultrasonic imaging of the water flow is performed to obtain ultrasound particle images. Finally, the ultrasonic particle images were analyzed using the algorithm proposed in this paper to obtain a two-dimensional flow field of the turbid water stream.

5. Result and Error Analysis

The algorithm proposed in the paper is applied to the acquired ultrasonic underwater images to calculate the underwater two-dimensional flow field to verify the effectiveness of the algorithm.

The original ultrasound particle image of $640\times 480$ pixels is downsampled by 2:1 to obtain the image of the next layer. The bilinear interpolation method is used between adjacent layers to make the transition smoother. The obtained image sequence of the pyramid structure is shown in Figure 10.

The experimental results are measured by average angle error (AAE) and standard angle error (SAD). The average angle error is calculated as:

$$\mathrm{AAE}=\frac{1}{N}\sum _{i=1}^N{\varphi }_e\left(i\right)$$

(31)

where:

$${\varphi }_e\left(i\right)=arccos\left[\frac{u_i^c{u}_i^e+v_i^c{v}_i^e+k^2}{\sqrt{{\left(u_i^c\right)}^2+{\left(v_i^c\right)}^2+k^2}·\sqrt{{\left(u_i^e\right)}^2+{\left(v_i^e\right)}^2+k^2}}\right]$$

(32)

where N denotes the total number of pixels in a frame of ultrasound particle image, $(u_i^c,v_i^c)$ denotes the ground truth value of the velocity of the ith pixel, i.e., standard optical flow value, $(u_i^e,v_i^e)$ denotes the optical flow vector value of the ith pixel, and k denotes the number of frame intervals between images. ${\varphi }_i\left(i\right)$ denotes the angular error of the ith pixel, and AAE denotes the average angular error of the entire optical flow field.

The average angular error reflects the deviation of the optical flow vector field calculated using the optical flow algorithm from the standard vector field. In the experiments, the frame interval k is set as 1, and the velocity measured by a high precision acoustic doppler velocimetry (ADV) is used as the standard optical flow value $u_i^e,v_i^e$, to evaluate the depth average flow velocity of the flow field obtained by the algorithm in the paper. Since the real underwater flow field is difficult to acquire, this verification is only used as a consistency verification.

The standard angle error can be expressed as:

$$\mathrm{SAD}=\sqrt{\frac{1}{n}\sum _{i=1}^n{({\varphi }_e\left(i\right)-\mathrm{AAE})}^2}$$

(33)

The standard angular error reflects the fluctuation of AAE in the calculated optical flow vector field.

The performance of Pyramid L-K optical flow and the traditional algorithm, Cross-Correlation algorithm, are compared horizontally to further validate the performance of the proposed algorithm. The experiments choose a high concentration ultrasonic particle image sequence of 50 consecutive frames of ultrasonic particle images with a volume concentration of 10‰ and a flow rate of $25\mathrm{cm}/\mathrm{s}$ for experimental analysis, as shown in Figure 11, to acquire a time-averaged flow field. The experiments are conducted with typical templates with sizes of $32\times 32$, $48\times 48$, and $64\times 64$ pixels are selected for matching in the template matching-based mutual correlation algorithm with window deformation technique [4]. The number of stratifications of 0, 1, and 2 are selected for comparative analysis in the Pyramid L-K optical flow-based algorithm, as shown in

Table 1. Performance Comparison of Pyramid L-K, Cross-Correlation algorithms.

Algorithm	Parameter	AAE	SAD	Time Consuming
cross-correlation	32 × 32	12.14	21.04	6.2
	48 × 48	14.85	23.86	5.5
	64 × 64	15.73	24.43	4.7
H-S algorithm	$\alpha =0.25$	3.92	9.11	66
	$\alpha =0.5$	5.31	8.54	62
	$\alpha =0.75$	7.04	7.99	65
Pyramid L-K	0 Layer	4.59	6.85	54
	1 Layer	2.43	5.33	28
	2 Layers	5.12	8.39	17

. Figure 12(a1,b1) shows the velocity field obtained using the Pyramid L-K optical flow, Cross-Correlation algorithm, and Figure 12(a2,b2) shows the velocity cloud maps corresponding to the velocity fields.

The velocity vector fields obtained after processing by each algorithm are shown in Figure 12(a1–c1), from which it can be intuitively seen that the flow spectra obtained by the H-S optical flow and Pyramid L-K optical flow based algorithms are better and the flow fields are smoother with only a few error matches, while the results obtained by the cross-correlation algorithm based on template matching have more error vectors and the velocity vector fields are not smooth enough. On the other hand, Figure 12(a2–c2) qualitatively compare the errors of H-S optical flow, Pyramid L-K optical flow, and Cross-Correlation algorithm. From that it can be seen that the velocity cloud map of the cross-correlation algorithm has the largest number and area of read spot regions, which reflects the largest error of the algorithm. It is due to the high concentration of particles in the sand-bearing water flow, which makes the correlation plane emerge with multiple peaks and leads to a large number of erroneous displacement vectors. There are some small red spots which appear in the lower left corner of the velocity cloud from the H-S optical flow algorithm. These indicate that the error is slightly larger in this area than the others. Because the concentration of sand in this region is high, the secondary reflection of ultrasound in these places is more serious, and the grayscale of the ultrasound image changes more drastically. Under this condition, the basic constraint of the H-S optical flow algorithm, i.e., the assumption of grayscale conservation, is not strictly valid, resulting large errors. The velocity cloud map from the Pyramid L-K optical flow algorithm is the smoothest and with the smallest error. However, the errors are larger at the edges of the observation area than other areas, which is because the L-K optical flow algorithm uses local constraints so that the calculation of the optical flow at each pixel needs to depend on the optical flow values of other pixels in the neighborhood. While the grayscale values at the edge of the image are zero, the optical flow in those regions can not be calculated accurately.

To quantitatively analyze the processing results of various algorithms, Table 1 lists the mean angular error, standard angular error and operation time of the H-S optical flow, Pyramid L-K optical flow and FFT Cross-Correlation algorithms. From the table, it can be seen that the FFT Cross-Correlation is the fastest, but its average angular error and standard angular error are the largest, too. The H-S algorithm has the longest operation time, but its accuracy is better than the FFT Cross-Correlation algorithm. By adjusting the smoothing factor, the average angular error can be reduced to 3.92. The Pyramid L-K algorithm has a better calculation efficiency than the H-S algorithm, and the best accuracy among the three algorithms. It shows that the Pyramid L-K algorithm can achieve the good balance of computational efficiency and accuracy by choosing the appropriate number of layers.

It can also be seen from Table 1 that when the template size increases from $32\times 32$ to $64\times 64$, the calculated average angular error does not change much. The Cross-Correlation algorithm is not very sensitive to its parameter. However, the H-S optical flow and Pyramid L-K optical flow algorithms are particularly sensitive to their parameters. Therefore, it is necessary to choose carefully when using parameters to keep the error from increasing.

For verifying the effectiveness of the Pyramid L-K optical flow algorithm in the case of large flow velocities, experiments are conducted using a sequence of 50 consecutive frames of ultrasonic particle images with a volume concentration of 10‰ and a flow velocity of $35\mathrm{cm}/\mathrm{s}$, as shown in Figure 13, to acquire a time-averaged flow field. Figure 14(a1–c1) shows the optical flow field obtained using the Pyramid L-K optical flow algorithm for the number of stratifications 0, 1, and 2, and Figure 14(a2–c2) shows the velocity clouds corresponding to the optical flow field.

From Figure 14(a1–c1), we can see that the Pyramid L-K optical flow algorithm can calculate the optical flow at each point in the image observation area. However, with the increase in the number of layers, the optical flow calculation at the boundary of the observation area becomes worse. This is caused by the edge effect of the Pyramid L-K optical flow algorithm as discussed previously.

Table 2. Comparison of optical flow calculation error of ultrasound particle image sequence with volume concentration of 10‰ and flow velocity of 35 cm/s.

Layer	AAE	SAD	Time-Consuming
0	4.87	7.13	56
1	3.74	5.64	29
2	5.43	8.67	16

shows the mean angular error, standard angular error and operation time of the optical flow vector field calculated by the proposed ultrasound image velocimetry algorithm. From the results, it can be seen that the mean angular error decreases from 4.87 to 3.74, and the standard angular error decreases from 7.13 to 5.64 when the number of stratification increases from 0 to 1. When the number of stratification increases from 1 to 2, the mean angular error increases from 3.74 to 5.43 and the standard angular error increases from 5.64 to 8.67, which is because the basic constraints of the L-K algorithm are: motion velocity is small and the local pixel motion is consistency.

6. Conclusions

This paper describes a ultrasonic PIV algorithm to measure instantaneous two-dimensional velocity field and depth-average velocity for sediment-laden fluid. Experiments are conducted in turbid and sandy water flow to measure and analyze the flow velocity field with the algorithm proposed in the paper. From the experimental results, it can be seen that the Pyramid L-K optical flow algorithm can satisfy the optical flow constraint equation even at large flow velocities, because of the downsampling of the original image to form a pyramid structure. The algorithm also adopts the “Coarse to Fine” matching, which finds the large motions in the higher level pyramid images first, and then gradually refines the result, correcting the accuracy of the large motions. The measure not only improves the computational accuracy but also reduces the computing time. The experimental results prove that the Pyramid L-K optical flow algorithm has the advantages of high computational accuracy, a wide range of application, and efficient operation speed.