Method and device for calculating river surface flow velocity based on variational principle

Abstract

A method and device for calculating a river surface flow velocity are provided based on a variational principle, which are used to capture and process the images of an objective area, and to obtain the flow velocity field data of the objective area with high precision in a non-contact manner. The method and device include 3 steps: (1) preparation before initial flow measurement; (2) capturing a video of the river by an image acquisition device, converting a motion of a pixel flow field of the fluid in a captured image sequence into solving an energy functional optimization problem, and solving partial differential equations to obtain data of pixel flow field distribution; and (3) obtaining space coordinates of the pixel point in a world coordinate system and calculating the flow velocity according to the data obtained in the step 2 and the transformation relationship determined in the step 1.

Claims

1. A method for calculating a river surface flow velocity based on a variational principle, comprising the following steps: step 1: preparation before initial flow measurement: setting up a plurality of calibration boards on left bank and right bank slopes of a river according to characteristics of a site environment, and measuring and recording coordinates of the plurality of calibration boards using a measuring equipment; setting up an image acquisition device, and adjusting a shooting angle and a focal length such that a field of view is capable of covering the plurality of calibration boards within a range and a motion of a fluid in the river; determining a transformation relationship from three-dimensional coordinates of an image sequence to coordinates of a pixel point of an image based on positions of the plurality of calibration boards; step 2: capturing a video of the river by the image acquisition device, and converting a motion of a pixel flow field of the fluid in a captured image sequence into solving the following energy functional extreme value problem; $\begin{array}{cc}E\left(\omega \left(i,j\right)\right)=\int \Phi \left(\frac{\partial I}{\partial t}+\nabla I^T\omega -D\Delta I\right)+\varphi \left(\nabla \omega \right)d\omega & \left(1\right)\end{array}$ wherein i represents a set of pixels in an i-th row in the image, j represents a set of pixels in a j-th column in the image, t is a shooting time of the image, ω(i,j) represents a pixel flow field distribution to be solved, I=I(i,j,t) represents the captured image sequence, and D is a diffusion coefficient; ∇ is a gradient operator, Δ is a Laplacian operator, and α is a weight coefficient; Φ, ϕ are functions respectively for a flow field and a gradient; a general form of an objective functional is: $\begin{array}{cc}E\left(\omega \left(i,j\right)\right)=\int {\left(\frac{\partial I}{\partial t}+\nabla I^T\omega -D\Delta I\right)}^2+\alpha \nabla {\omega }^T\nabla \omega d\omega & \left(2\right)\end{array}$ the pixel flow field distribution ω(i,j) is solved; and step 3: obtaining space coordinates of the pixel point in a world coordinate system according to data of the pixel flow field distribution ω(i,j) calculated in the step 2 and the transformation relationship determined in the step 1, and further obtaining a flow velocity $v=\frac{d}{\Delta t},$  wherein d is a length of a flow velocity vector in ω projected to world coordinates, and Δt is a time difference between frames of calculated images; wherein the flow velocity is used for measuring a river flow and hydrometry.

2. The method for calculating the river surface flow velocity based on the variational principle according to claim 1, wherein: in the step 1, first, a flow measurement site is investigated such that a selected position has a sufficient field of view to observe 80% or more of an area of the river; and a power supply infrastructure required for video recording and video transmission is configured.

3. The method for calculating the river surface flow velocity based on the variational principle according to claim 1, wherein: in the step 1, the coordinates of the pixel point in the image are marked as (i,j), the coordinates of a corresponding point in a three-dimensional space are marked as (x,y,z), and an n-th undetermined parameter of a camera calibration is marked as l.sub.n; a calculation is carried out according to the following transformation formula: $\begin{array}{cc}i=\frac{l_{1}x+l_{2}y+l_3{z}_w+l_4}{l_{9}x+l_{10}y+l_{11}{z}_w+1},j=\frac{l_{5}x+l_{6}y+l_7{z}_w+l_8}{l_{9}x+l_{10}y+l_{11}{z}_w+1}& \left(3\right)\end{array}$ mapping of pixel coordinates projected to the world coordinates is represented by $\mathrm{Proj}\left(\left[\begin{array}{c}i\\ j\end{array}\right]\right).fwdarw.\left[\begin{array}{c}x\\ y\\ z\end{array}\right],$  mapping of the world coordinates projected to the pixel coordinates is represented by ${\mathrm{Proj}}^{-1}\left(\left[\begin{array}{c}x\\ y\\ z\end{array}\right]\right).fwdarw.\left[\begin{array}{c}i\\ j\end{array}\right],$  and the above equation is solved based on position coordinates corresponding to the plurality of calibration boards, which are preset, to obtain the transformation relationship from the three-dimensional coordinates of the image sequence to the coordinates of the pixel point of the image.

4. The method for calculating the river surface flow velocity based on the variational principle according to claim 1, wherein: in the step 1, a number of the plurality of calibration boards is not less than 6.

5. The method for calculating the river surface flow velocity based on the variational principle according to claim 1, wherein: in the step 1, the river is filmed by a camera, the video is transmitted to a PC terminal through bridging equipment for preliminary review, the video is checked for a number of frames and damaged frames, a satisfactory video sequence is transferred to a preprocessing stage, and an unsatisfactory video sequence is deleted.

6. The method for calculating the river surface flow velocity based on the variational principle according to claim 1, wherein: in the step 2, a Euler-Lagrange equation corresponding to a flow velocity field is: $\begin{array}{cc}\nabla I\left(\frac{\partial I}{\partial t}+\nabla I^T\omega -D\Delta I\right)-\mathrm{\alpha \Delta \omega }=0& \left(4\right)\end{array}$ for an image I, a coordinate in a horizontal direction is represented by u, and a coordinate in a vertical downward direction is represented by v, and for a numerical discretization, discretization schemes of the gradient of the I are: $\begin{array}{cc}{\left(\frac{\partial I}{\partial u}\right)}_{i,j}=\{\begin{array}{ccc}{I}_{i,j+1}-I_{i,j}& ,& j<w\\ 0& ,& j=w\end{array}& \left(5\right)\end{array}$ $\begin{array}{cc}{\left(\frac{\partial I}{\partial v}\right)}_{i,j}=\{\begin{array}{ccc}{I}_{i+1,j}-I_{i,j}& ,& i<h\\ 0& ,& i=h\end{array}& \left(6\right)\end{array}$ wherein w represents a width of the image, and h represents a height of the image; and a set of partial differential equations shown in formula (4) is iteratively solved through formula (7) until convergence, and a motion direction and a velocity of the pixel point of the river are inverted: $\begin{array}{cc}{\omega }^{k+1}={\omega }^k-\theta \left\{\nabla I^k\left(\frac{\partial I^k}{\partial t}+\nabla {\left(I^k\right)}^T{\omega }^k-D\Delta I^k\right)-{\mathrm{\alpha \Delta \omega }}^k\right\}& \left(7\right)\end{array}$ wherein θ is a step parameter, and k represents a count of a current round of an iteration.

7. A device for calculating a river surface flow velocity based on a variational principle, the device is configured: to acquire and record coordinates of calibration boards by an image acquisition device whose field of view is capable of covering the calibration boards within a range and a motion of a fluid in a river, and then determine a transformation relationship from three-dimensional coordinates of an image sequence to coordinates of a pixel point of an image based on positions of the calibration boards; to capture a video of the river by the image acquisition device, and convert a motion of a pixel flow field of the fluid in a captured image sequence into solving the following energy functional extreme value problem; $\begin{array}{cc}E\left(\omega \left(i,j\right)\right)=\int \Phi \left(\frac{\partial I}{\partial t}+\nabla I^T\omega -D\Delta I\right)+\varphi \left(\nabla \omega \right)d\omega & \left(1\right)\end{array}$ wherein i represents a set of pixels in an i-th row in the image, j represents a set of pixels in a j-th column in the image, t is a shooting time of the image, ω(i,j) represents a pixel flow field distribution to be solved, I=I(i,j,t) represents the captured image sequence, and D is a diffusion coefficient; ∇ is a gradient operator, Δ is a Laplacian operator, and α is a weight coefficient; ϕ and ϕ are functions respectively for a flow field and a gradient; a general form of the objective functional is: $\begin{array}{cc}E\left(\omega \left(i,j\right)\right)=\int {\left(\frac{\partial I}{\partial t}+\nabla I^T\omega -D\Delta I\right)}^2+\alpha \nabla {\omega }^T\nabla \omega d\omega & \left(2\right)\end{array}$ the pixel flow field distribution w(i,j) is solved; and to obtain space coordinates of the pixel point in a world coordinate system according to a solving result of the solving part and the transformation relationship determined by the device, and further obtain a flow velocity $v=\frac{d}{\Delta t},$  wherein d is a length of a flow velocity vector in ω projected to world coordinates, and Δt is a time difference between frames of calculated images; wherein the flow velocity is used for measuring a river flow and hydrometry.

8. The device for calculating the river surface flow velocity based on the variational principle according to claim 7, wherein: the coordinates of the pixel point in the image are marked as (i,j); the coordinates of a corresponding point in a three-dimensional space are marked as (x,y,z), and an n-th undetermined parameter of a camera calibration is marked as l.sub.n; a calculation is carried out according to the following transformation formula: $\begin{array}{cc}i=\frac{l_{1}x+l_{2}y+l_3{z}_w+l_4}{l_{9}x+l_{10}y+l_{11}{z}_w+1},j=\frac{l_{5}x+l_{6}y+l_7{z}_w+l_8}{l_{9}x+l_{10}y+l_{11}{z}_w+1}& \left(3\right)\end{array}$ mapping of pixel coordinates projected to the world coordinates is represented by $\mathrm{Proj}\left(\left[\begin{array}{c}i\\ j\end{array}\right]\right).fwdarw.\left[\begin{array}{c}x\\ y\\ z\end{array}\right],$  mapping of the world coordinates projected to the pixel coordinates is represented by ${\mathrm{Proj}}^{-1}\left(\left[\begin{array}{c}x\\ y\\ z\end{array}\right]\right).fwdarw.\left[\begin{array}{c}i\\ j\end{array}\right],$  and the above equation is solved based on position coordinates corresponding to the calibration boards, which are preset to obtain the transformation relationship from the three-dimensional coordinates of the image sequence to the coordinates of the pixel point of the image.

9. The device for calculating the river surface flow velocity based on the variational principle according to claim 7, wherein: a Euler-Lagrange equation corresponding to a flow velocity field is: $\begin{array}{cc}\nabla I\left(\frac{\partial I}{\partial t}+\nabla I^T\omega \right)-D\Delta I)-\mathrm{\alpha \Delta \omega })=0& \left(4\right)\end{array}$ for an image I, a coordinate in a horizontal direction is represented by u, and a coordinate in a vertical downward direction is represented by v, and for numerical discretization, discretization schemes of the gradient of the I are: $\begin{array}{cc}{\left(\frac{\partial I}{\partial u}\right)}_{i,j}=\{\begin{array}{ccc}{I}_{i,j+1}-I_{i,j}& ,& j<w\\ 0& ,& j=w\end{array}& \left(5\right)\end{array}$ $\begin{array}{cc}{\left(\frac{\partial I}{\partial v}\right)}_{i,j}=\{\begin{array}{ccc}{I}_{i+1,j}-I_{i,j}& ,& i<h\\ 0& ,& i=h\end{array}& \left(6\right)\end{array}$ wherein w represents a width of the image, and h represents a height of the image; and a set of partial differential equations shown in formula (4) is iteratively solved through formula (7) until convergence, and a motion direction and a velocity of the pixel point of the river are inverted: $\begin{array}{cc}{\omega }^{k+1}={\omega }^k-\theta \left\{\nabla I^k\left(\frac{\partial I^k}{\partial t}+\nabla {\left(I^k\right)}^T{\omega }^k-D\Delta I^k\right)-{\mathrm{\alpha \Delta \omega }}^k\right\}& \left(7\right)\end{array}$ wherein θ is a step parameter, and k represents a count of a current round of an iteration.

10. The device for calculating the river surface flow velocity based on the variational principle according to claim 7, the device is further configured: to allow a user to input an operation instruction and display the solving result.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is a flow chart of a method for calculating a river surface flow velocity based on a variational principle according to the present disclosure;

(2) FIG. 2 is a schematic diagram of camera calibration according to the present disclosure; and

(3) FIG. 3 is a schematic diagram of a pixel flow field according to the present disclosure.

DESCRIPTION OF THE EMBODIMENTS

(4) Specific implementations of the method and device for calculating the river surface flow velocity based on the variational principle according to the present disclosure will be described below in detail in conjunction with the accompanying drawings.

(5) As shown in FIG. 1, the method for calculating the river surface flow velocity based on the variational principle provided by this embodiment may be divided into three parts: early-stage preparation for calibration before initial flow measurement, river video capture, and river surface flow velocity calculation. The method specifically includes the following steps:

(6) Step 1: Preparation Before Initial Flow Measurement.

(7) 1.1) Early-stage preparation before initial flow measurement. First, a flow measurement site is investigated such that a selected position has a sufficient field of view to observe 80% or more of an area of the river and avoids environments that may easily damage the equipment. A power supply and transmission infrastructure required for video recording and video transmission is configured. The infrastructure includes high-definition cameras, power supplies, fixed carriers, etc.; bridging equipment for transmitting data; and remote PC equipment, and destination boards and measuring instruments for calibration. Before the flow measurement, the measuring instruments such as a total station are used to carry out related calibration and positioning processing.

(8) 1.2) Not less than 6 calibration boards are set up on left bank and right bank slopes of a river according to characteristics of a site environment. The calibration boards are designed with high-contrast signs. After the calibration boards are placed properly, coordinates of the calibration boards are measured and recorded using a measuring equipment such as a total station.

(9) 1.3) River video capture. After a camera is mounted on a fixed platform for shooting, a shooting angle and a focal length are adjusted such that a field of view is capable of covering the calibration boards within a range and a motion of a fluid in the river. The river is filmed by the camera, the video is transmitted to a PC terminal through the bridging equipment for preliminary review, and the video is checked for a number of frames and damaged frames. A satisfactory video sequence is transferred to a preprocessing stage, and an unsatisfactory video sequence is deleted.

(10) 1.4) On-site calibration is carried out. A transformation matrix and a relationship between world coordinates and pixel coordinates are obtained.

(11) In order to determine the relationship between the pixel point coordinates and the three-dimensional coordinates, internal and external parameters of a camera model are determined, and a DLT linear transformation method is introduced. The at least 6 calibration boards are set up on site. The calibration boards should be clearly visible in the field of view, as shown in FIG. 2. Coordinates of the pixel point in the image are marked as (i,j), where i represents a set of pixels in an i-th row in the image, and j represents a set of pixels in a j-th column in the image. The coordinates of a corresponding point in a three-dimensional space are marked as (x,y,z), and an n-th undetermined parameter of a camera calibration is marked as l.sub.n. A calculation is carried out according to the following transformation formula:

(12) $i=\frac{l_{1}x+l_{2}y+l_3{z}_w+l_4}{l_{9}x+l_{10}y+l_{11}{z}_w+1},j=\frac{l_{5}x+l_{6}y+l_7{z}_w+l_8}{l_{9}x+l_{10}y+l_{11}{z}_w+1}.$

(13) For convenience of introduction, mapping of pixel coordinates projected to the world coordinates is represented by

(14) 0$\mathrm{Proj}\left(\left[\begin{array}{c}i\\ j\end{array}\right]\right).fwdarw.\left[\begin{array}{c}x\\ y\\ z\end{array}\right].$
Mapping of the world coordinates projected to the pixel coordinates is represented by

(15) ${\mathrm{Proj}}^{-1}\left(\left[\begin{array}{c}x\\ y\\ z\end{array}\right]\right).fwdarw.\left[\begin{array}{c}i\\ j\end{array}\right].$
The above equation is solved based on the preset at least 6 calibration points to complete the transformation from the three-dimensional coordinates of the image sequence to the coordinates of the pixel point of the image.

(16) Step 2: A video of the river is captured by the image acquisition device, and a motion estimation of a pixel flow field of the fluid in the captured image sequence is converted into solving the following energy functional extreme value problem:

(17) $E\left(\omega \left(i,j\right)\right)=\int \Phi \left(\frac{\partial I}{\partial t}+\nabla I^T\omega -D\Delta I\right)+\varphi \left(\nabla \omega \right)d\omega ,$

(18) where i represents a set of pixels in an i-th row in the image, j represents a set of pixels in a j-th column in the image, t is a shooting time of the image, ω(i,j) represents a pixel flow field distribution to be solved, I=I(i,j,t) represents the captured image sequence, and D is a diffusion coefficient. ∇ is a gradient operator, Δ is a Laplacian operator, and α is a weight coefficient. ϕ and ϕ are functions respectively for a flow field and a gradient, which may be square functions, Huber functions and the like. The following is an introduction in a case where ϕ and ϕ are the L.sub.2 norm, and other similar loss function combinations are similar. At this time, the form of the objective functional is:

(19) $E\left(\omega \left(i,j\right)\right)=\int {\left(\frac{\partial I}{\partial t}+\nabla I^T\omega -D\Delta I\right)}^2+\alpha \nabla {\omega }^T\nabla \omega d\omega .$

(20) A Euler-Lagrange equation corresponding to a flow velocity field is:

(21) $\nabla I\left(\frac{\partial I}{\partial t}+\nabla I^T\omega -D\Delta I\right)-\alpha \Delta \omega =0.$

(22) Taking the numerical value of the gradient descent as an example, the energy functional extreme value problem is solved to obtain the pixel flow field distribution ω(i,j). For an image I, a coordinate in a horizontal direction is represented by u, and a coordinate in a vertical downward direction is represented by v. For numerical discretization, taking finite difference as an example, discretization schemes of the gradient of I are:

(23) ${\left(\frac{\partial I}{\partial u}\right)}_{i,j}=\{\begin{array}{cc}{I}_{i,j+1}-I_{i,j},& j<w\\ 0,& j=w\end{array};$${\left(\frac{\partial I}{\partial v}\right)}_{i,j}=\{\begin{array}{cc}{I}_{i+1,j}-I_{i,j},& i<h\\ 0,& i=h\end{array}.$

(24) The equations are iteratively solved through the following iterative scheme until convergence, where θ is a step parameter. The obtained pixel flow field result is shown in FIG. 3. The motion direction and velocity of the pixel point of the river are inverted.

(25) ${\omega }^{k+1}={\omega }^k-\theta \left\{\nabla I^k\left(\frac{\partial I^k}{\partial t}+\nabla {\left(I^k\right)}^T{\omega }^k-D\Delta I^k\right)-{\mathrm{\alpha \Delta \omega }}^k\right\},$

(26) where θ is the step parameter, w represents a width of the image, h represents a height of the image, and k represents a count of the current round of iteration.

(27) Step 3: Flow velocity calculation.

(28) After the flow velocity field of the pixel in the step 2 is obtained, projective transformation is carried out in conjunction with the projection relationship between the world coordinate system and the pixel coordinate system in the step 1.4, that is,

(29) $\mathrm{Proj}\left(\left[\begin{array}{c}i\\ j\end{array}\right]\right).fwdarw.\left[\begin{array}{c}x\\ y\\ z\end{array}\right]\mathrm{and}\mathrm{Proj}\left(\left[\begin{array}{c}i\\ j\end{array}\right]\right).fwdarw.\left[\begin{array}{c}x\\ y\\ z\end{array}\right].$
Space coordinates of the pixel point in a world coordinate system are obtained, and the flow velocity

(30) $v=\frac{d}{\Delta t}$
may be further obtained, where v is the flow velocity in m/s. d is the length of the flow velocity vector in ω projected to the world coordinates in m. Δt is a time difference between frames of the calculated images in s.

(31) Further, this embodiment further provides a device for calculating a river surface flow velocity capable of automatically implementing the above method. The device includes a transformation relationship acquisition part, a solving part, a flow velocity calculation part, an input and display part and a control part.

(32) The transformation relationship acquisition part is configured to acquire and record coordinates of calibration boards by an image acquisition device whose field of view is capable of covering the calibration boards (calibration points) within a range and a motion of a fluid in the river, and then determine a transformation relationship from three-dimensional coordinates of the image sequence to coordinates of the pixel point of the image based on positions of the calibration boards.

(33) Specifically, coordinates of the pixel point in the image are marked as (i,j), where i represents a set of pixels in an i-th row in the image, and j represents a set of pixels in a j-th column in the image. The coordinates of a corresponding point in the three-dimensional space are marked as (x,y,z), and an n-th undetermined parameter of camera calibration is marked as l.sub.n. Calculation is carried out according to the following transformation formula:

(34) $i=\frac{l_{1}x+l_{2}y+l_3{z}_w+l_4}{l_{9}x+l_{10}y+l_{11}{z}_w+1},j=\frac{l_{5}x+l_{6}y+l_7{z}_w+l_8}{l_{9}x+l_{10}y+l_{11}{z}_w+1}.$

(35) Mapping of pixel coordinates projected to the world coordinates is represented by

(36) 0$\mathrm{Proj}\left(\left[\begin{array}{c}i\\ j\end{array}\right]\right).fwdarw.\left[\begin{array}{c}x\\ y\\ z\end{array}\right],$
and mapping of the world coordinates projected to the pixel coordinates is represented by

(37) ${\mathrm{Proj}}^{-1}\left(\left[\begin{array}{c}x\\ y\\ z\end{array}\right]\right).fwdarw.\left[\begin{array}{c}i\\ j\end{array}\right].$
The above equation is solved based on the preset at least 6 calibration points to complete the transformation from the three-dimensional coordinates of the image sequence to the coordinates of the pixel point of the image.

(38) The solving part is configured to capture a video of the river by the image acquisition device, and convert a motion of a pixel flow field of the fluid in the captured image sequence into solving the following energy functional extreme value problem.

(39) $E\left(\omega \left(i,j\right)\right)=\int \Phi \left(\frac{\partial I}{\partial t}+\nabla I^T\omega -D\Delta I\right)+\varphi \left(\nabla \omega \right)d\omega$

(40) where ω(i,j) represents a pixel flow field distribution to be solved, I=I(i,j,t) represents the captured image sequence, and D is a diffusion coefficient. ∇ is a gradient operator, Δ is a Laplacian operator, and α is a weight coefficient. ϕ and ϕ are functions respectively for the flow field and a gradient. The general form of the objective functional is:

(41) $E\left(\omega \left(i,j\right)\right)=\int {\left(\frac{\partial I}{\partial t}+\nabla I^T\omega -D\Delta I\right)}^2+\alpha \nabla {\omega }^T\nabla \omega d\omega .$

(42) The pixel flow field distribution ω(i,j) is solved. Specifically, a Euler-Lagrange equation corresponding to a flow velocity field is:

(43) $\nabla I\left(\frac{\partial I}{\partial t}+\nabla I^T\omega -D\Delta I\right)-\mathrm{\alpha \Delta \omega }=0.$

(44) For an image I, a coordinate in a horizontal direction is represented by u, and a coordinate in a vertical downward direction is represented by v, and for numerical discretization, discretization schemes of the gradient of I are:

(45) $$\begin{gathered} {\left(\frac{\partial I}{\partial u}\right)}_{i,j}=\{\begin{array}{ccc}{I}_{i,j+1}-I_{i,j},& & j<w\\ 0,& & j=w\end{array} \\ {\left(\frac{\partial I}{\partial v}\right)}_{i,j}=\{\begin{array}{ccc}{I}_{i+j,1},-I_{i,j},& & i<h\\ 0,& & i=h\end{array} \end{gathered}$$

(46) where w represents a width of the image, and h represents a height of the image.

(47) The equations are iteratively solved through the following formula until convergence, and the motion direction and velocity of the pixel point of the river are inverted:

(48) ${\omega }^{k+1}={\omega }^k-\theta \left\{\nabla I^k\left(\frac{\partial I^k}{\partial t}+\nabla {\left(I^k\right)}^T{\omega }^k-D\Delta I^k\right)-{\mathrm{\alpha \Delta \omega }}^k\right\}$

(49) where θ is a step, and k represents a count of the current round of iteration.

(50) The flow velocity calculation part is configured to obtain space coordinates of the pixel point in a world coordinate system according to the solving result of the solving part and the transformation relationship determined by the transformation relationship acquisition part, and further obtain the flow velocity

(51) $v=\frac{d}{\Delta t},$
where d is a length of a flow velocity vector in w projected to the world coordinates, and Δt is a time difference between frames of the calculated images.

(52) The input and display part communicates with the transformation relationship acquisition part, the solving part, the flow velocity calculation part and the control part, and is configured to allow a user to input an operation instruction and display a corresponding result.

(53) The control part communicates with the transformation relationship acquisition part, the solving part, the flow velocity calculation part, and the input and display part, and is configured to control operations thereof.

(54) The above embodiment is merely an example to illustrate the technical solution of the present disclosure. The method and device for calculating a river surface flow velocity based on a variational principle according to the present disclosure are not limited to the content described in the above embodiment, but are subject to the scope defined in the claims. Any modifications, additions or equivalent substitutions made by those skilled in the art of the present disclosure on the basis of this embodiment are within the protection scope of the claims of the present disclosure.