Radial basis function-based vector field algorithm for wildfire boundary tracking with UAVs

Abstract

This paper tackles the problem of dynamic wildfire boundary tracking with UAVs. Wildfire boundary is treated as the zero-level set curve of an implicit function and is approximated with radial basis functions. Its propagation is modeled with the Hamilton–Jacobi equation with an arbitrary initial boundary as the input. To navigate UAVs to the wildfire boundary, an analytical velocity vector field, whose integral curves converge to the wildfire boundary, is constructed on the basis of the typical radial basis function thin-plate spline. Computer simulations with a single UAV and multiple UAVs have been conducted for the evaluation of the proposed solution, and numerical results show that the proposed algorithm can ensure the successful tracking of an arbitrarily shaped wildfire boundary.

1 Introduction

Wildfire, also known as forest fire or bushfire, is an area of combustible vegetation that has caught fire and requires significant amounts of resources for its suppression. Wildfires are becoming more and more prevalent throughout the world and cause considerable ecological and economical loss to the wildland environment (Ambroz et al. 2019). In the study of combating wildfires, one main concern is the live wildfire boundary monitoring, but it is risky and costly to adopt manned aerial wildfire monitoring, especially when dealing with uncontrolled fires. Given these situations, it is essential to adopt automated approaches to monitor wildfire frontiers, which can help improve significantly the wildfire suppression efficiency, risk management, and policy implementation.

Unmanned aerial vehicles (UAVs), due to the ease of deployment, low maintenance cost, high mobility, and flexibility, have gained considerable attention in numerous military and civilian applications, such as ground traffic surveillance (Kanistras et al. 2013), geological surveying (Adams and Friedland 2011), reconnaissance, and tracking (Obermeyer et al. 2010), to name a few. In the context of forestry, autonomous UAVs have been proven to be an effective tool for studying environmental features and surveying hazardous phenomena-prone places (Tang and Shao 2015). A survey in Yuan et al. (2015) analyzes the sensing hardware and algorithms for wildfire detection. The work in Twidwell et al. (2016) presents a comprehensive analysis of the usage of UAVs in fire management operations. Survey in Bailon-Ruiz and Lacroix (2020) reviews the state of the art on wildfire remote sensing with UAVs from the perspective of autonomy.

In recent years, UAV-based monitoring has become a prominent study area and plays a more integral role in providing situational awareness (SA) (Stanton et al. 2001), which refers to perception of elements of interest in the environment within a volume of time and space, and the comprehension of their influences on the evolution of the event state over time. Although a great number of monitoring algorithms (Bertozzi et al. 2005; Susca et al. 2008; Cassandras et al. 2011; Lan and Schwager 2013; Smith et al. 2011; Smith and Rus 2010) have emerged, fire frontier monitoring is a largely under investigated domain. The focus of this paper is to develop a strategy for a UAV to track the boundary of an ongoing wildfire event, with the aim of providing wildfire fighters with the up-to-date information of active wildfire frontiers. The major components of this paper include using an easy algorithm to reconstruct the wildfire boundary based on initial information from sensing instruments, such as airborne sensors, and to efficiently solve the wildfire propagation equation. Based on the reconstructed wildfire boundary, develop a velocity vector field for UAV guidance control. In line with the main research points, the main novelties of this paper include: (1) The investigation into applying radial basis functions to construct a level set function from sensor data to describe the wildfire boundary, (2) the application of radial basis functions to estimate the propagated wildfire boundary modeled by Hamilton–Jacobi equation, and (3) the design of an analytical velocity vector field from radial basis functions to guide the UAV to track the wildfire boundary.

To the best knowledge of the authors, there are no other publications presenting such strategies to track dynamic wildfire boundaries. The remainder of the paper is organized as follows. Section 2 surveys the previous literature and establishes the current state of similar research. Section 3 describes the problem statement, assumptions, and models used in this work. Section 4 describes the strategies for wildfire boundary construction and evolution and presents a guidance methodology. Section 5 is devoted to the presentation of simulation results collected to assess the effectiveness and efficiency of the proposed algorithm. Section 6 closes the paper with a discussion of the results and future work.

2 Related work

In this section, we review current state of similar research. The dynamic nature of wildfire propagation poses challenges to keep a continuous and exhaustive surveillance of the advancing fire boundary. To effectively track the dynamic wildfire boundary, it is imperative to understand the mechanism of wildfire growth, and accurately predict its propagation. In the existing literature, a great deal of methods have been developed to estimate the propagation of wildfire boundaries, such as Huygen’s principle (Richards 1990), Lagrangian (direct) approach (Balažovjech and Mikula 2011; Dziuk 1999; Hou et al. 1994), and the Eulerian (level set) approach (Ambroz et al. 2019; Osher and Fedkiw 2006; Sethian 1999). However, when it comes to the issue of tracking wildfire boundary, the propagation of wildfire boundary is seldom considered. Only the work Kumar et al. (2011) encloses the points on the fire front with an elliptical envelop model and deployed UAVs to persistently covering the fire front under the UAV footprints. This paper concludes that the fire fronts follow the Huygen’s principle to propagate in an elliptical fashion when assuming a uniform fuel distribution, uniform landscape and weather, and constant wind direction. However, the initial shapes of wildfire boundaries are not constricted to an elliptical envelop. Therefore, the model presented in Kumar et al. (2011) is not suitable for an initial fire front with an arbitrary shape. The level set approach, devised by Osher and Sethian (1988), is highly robust and accurate to solve the problem of curve evolution under complex conditions. If we can approximate the initial wildfire boundary with a highly accurate level set function, the level set algorithm is a good choice for computing the advancing wildfire boundary.

The trajectory control problem, defined as making a vehicle follow a pre-determined course in space, is essential for safe and efficient operation of UAVs. The problem of tracking wildfire boundary requires that the UAV should reach the wildfire frontier and then follow the wildfire boundary as accurately as possible. In the past decades, a number of tracking algorithms (Bertozzi et al. 2005; Susca et al. 2008; Cassandras et al. 2011; Lan and Schwager 2013; Smith et al. 2011; Smith and Rus 2010) have emerged, whereas the wildfire boundary tracking is a largely unexplored problem. In the literature, boundary tracking is also named as level set tracking (Matveev et al. 2012), curve tracking (Malisoff et al. 2017), and isoline tracking (Dong et al. 2020). In these studies, the environment is always mapped with a scalar concentration function to represent any spatially varying variable, such as density of combustible material, temperature, etc. UAVs are required to track the constant level sets (boundary) determining the extent of the spatial field. These problems have been widely studied in the application of environmental exploration, for instance, tracking an oil or chemical spill in the ocean (Li et al. 2014a; Fahad et al. 2015; Jiang and Li 2018; Wang et al. 2019) and tracking ash cloud (Dong and You 2021). Roughly speaking, existing control algorithms for boundary tracking can be categorized as the gradient-based approach and the gradient-free approach, depending on whether the gradient information is available or not.

If the target boundary can be described as a geometric curve rather than a function of time, its tracking can be handled by path following algorithms. There are already many existing path following approaches, and some of them can be applied to deal with tracking dynamic boundaries. Among different path-following algorithms, vector field method is shown to achieve high path-following accuracy while requiring low control efforts (Sujit et al. 2014). As a typical gradient-based method, the gradient vector field algorithm decouples the guidance law into convergence and circulation two terms and weighs them by scalars to influence field strength and direction (Wilhelm and Clem 2019). The resulted vector field has the feature that its integral curve approaches the target path asymptotically, and it can ensure convergence of any trajectory that does not reach the “critical” locations where the vector filed is degenerate (Goncalves et al. 2010; Kapitanyuk et al. 2017). The study of vector field is extensive. The earliest application of vector field approach to the problem of isoline tracking can be found in work (Marthaler and Bertozzi 2003; Triandaf and Schwartz 2005), where UAVs are controlled to track a certain level of an explicitly formulated isoline. Work (Wilhelm and Clem 2019; Lim et al. 2014; Nelson et al. 2007) developed plane vector field path-following control laws for straight-line paths and circular arcs and orbits. To navigate a single mobile robot along a general n-dimensional smooth curve, the work Goncalves et al. (2009) provided a universal algorithm to compute an artificial vector field and extended it in Goncalves et al. (2010), where the implicit functions that define the target curve were used as the convergence term in the designed guidance law. Inspired by the idea of constructing a general vector field, work (Kapitanyuk et al. 2017; De Marina et al. 2017) designed vector field-based controllers to drive the UAV to track closed smooth curves described by explicit equations. In the application of vector field algorithms, the mathematical expression of the desired curve is essential for the construction of the vector field controller. In aforementioned studies, they all assumed that the explicit representation of the target curve is available by default. In gradient-based algorithms, the knowledge of the spatial field gradient is highly desirable. If the gradient information is not explicitly available, some work estimates gradient and their spatial divergence for designing guidance laws. Work (Li et al. 2014a; Fahad et al. 2015; Jiang and Li 2018; Wang et al. 2019) studied the problem of tracking a dynamic ocean plume whose dispersion is modeled by an advection–diffusion equation. As the tracked dynamic concentration-level curves have a real physical meaning, robots are assumed to be equipped with onboard sensors like fluorometer to obtain the measurement of the chemical concentration. By data processing, the estimated gradient of the concentration and the corresponding divergence of the concentration gradient are available. Thus, a control law only depending on concentration measurements was successfully constructed for tracking the chemicals poured in a marine environment.

If it is impossible to obtain the gradient information directly or by estimation, some studies turn to gradient-free algorithms. “Bang-Bang” approach is one of the most common and straightforward gradient-free methods and has been adopted to solve some boundary tracking problems (Jin and Bertozzi 2007; Nelson et al. 2007), however, the Bang-Bang type control switches between alternate steering angles if current measurement is above or below a predefined threshold, and results in unavoidable zigzagging behavior in the tracking process. Some work employed the sliding mode control technique to the tracking problem. In Matveev et al. (2012), the UAV equipped with a onboard sensor was limited to measure the field value of the environmental level set at its current position and was driven to track the static spatial level set with a sliding mode controller. A suboptimal sliding mode algorithm, designed with the density measurement of the contaminant from an onboard sensor of the autonomous vehicle, was proposed in Menon et al. (2014) to steer the single autonomous vehicle along the smooth contour of the static spatial contaminant phenomena. The work (Dong et al. 2020; Dong and You 2021) also assumed that the concentration measurement at the current location is in a point-wise fashion (Matveev et al. 2012) with onboard sensors, and presented a gradient-free controller in a PI-like form for a Dubins vehicle to track a desired isoline of an unknown scalar field using only the concentration feedback. It validated the effectiveness of the controller via simulations using a fixed-wing UAV on the real dataset of the environmental pollution level. For above-mentioned level curve tracking problems solved with gradient-free algorithms, the concentration value of the level curve has a physical meaning and must be available at the UAV’s current position with an onboard sensor.

Regardless of gradient-based or gradient-free controllers, most above-mentioned studies focus on tracking static level curves with a given explicit expression, rather than time-varying curves. As for variable curve tracking based on an advection–diffusion pollution dispersion model, it should be noticed that the concentration is meaningful and the UAV must measure the chemical concentration value at its current position for implementing the guidance law. In this paper, the wildfire boundary is modeled with level set function and its propagation is described by the Hamilton–Jacobi equation; nevertheless, the “concentration” value at the UAV’s position is unavailable as we do not have sensors that can measure a meaningful “concentration”. The only UAV data available are its current position and velocity measured by its onboard sensors. As for the tracked wildfire boundary, only the positions of waypoints on the wildfire boundary are available. To apply the idea of vector field algorithm for wildfire boundary tracking, it is necessary to deal with the scattered waypoints on the wildfire boundary with an effective and efficient approximation algorithm, thus obtaining an explicit level set function describing the wildfire boundary. What is required then is the construction of the vector field using the approximated level set function. In this paper, we consider tracking a time-varying wildfire boundary whose dependence on time cannot be neglected. Different from existing studies which approximate the wildfire boundary with an elliptical envelope, we apply a meshless numerical interpolation approach to reconstruct the wildfire boundary, regardless of the shape of the boundary. Also, the introduced meshless algorithm is incorporated with the Hamilton–Jacobi equation to exploit the propagation of wildfire boundary. In addition, to keep the integral of the velocity vector field converges to the desired curve, it will be derived from the approximated level set function by obeying the principle of gradient-based methods.

3 Problem formulation

Generally speaking, the main problem studied in this paper is to establish a control law for the UAV, such that it converges to, and follows the boundary of wildfires, in a specified direction such as clockwise or counterclockwise. In this section, first, the kinematic model of the UAV used in this paper is presented, and then, the details of the wildfire boundary reconstruction and propagation are provided. The following is a formal statement of the control problem.

3.1 UAV model

In the existing literature, it is assumed that the dynamic controller has a quick and precise response, thus considering a simpler kinematic model and dealing with the “outer” kinematic control. The kinematic model of the UAV used in this paper is described by the Dubins car model in 2D plane

$$\begin{aligned} {\dot{x}}_r= & {} \nu _r \ \text{ cos }\theta _r \nonumber \\ {\dot{y}}_r= & {} \nu _r \ \text{ sin }\theta _r \nonumber \\ {\dot{\theta }}_r= & {} \omega _r, \end{aligned}$$

(1)

where ${\varvec{x}}_r = [x_r,y_r]^\mathrm{T}$ and $\theta _r$ represent the global Cartesian coordinates and the heading angle of the UAV, and $v_r$ and $\omega _r$ are the corresponding translational and rotational velocities, respectively.

The position ${\varvec{x}}_r$ and the orientation $\theta _r$ in dynamics (1) cannot be simultaneously stabilized by a continuous static (time-invariant feedback) (Brockett 1983). To avoid this problem, for a positive constant $l_0$, a new variable ${\varvec{x}}$ is defined

$$\begin{aligned} {\varvec{x}} = [x,y]^\mathrm{T} = [x_r + l_0 \text{ cos }\theta _r, y_r+l_0\text{ sin }\theta _r]^\mathrm{T}, \end{aligned}$$

(2)

as the output of the system and the control objective is set based on the position of that output (Sahin et al. 2007). Using feedback linearization, the kinematic model (1) can be transformed into the single-integrator model (Li et al. 2014b)

$$\begin{aligned} \dot{{\varvec{x}}} = {\varvec{u}}, \end{aligned}$$

(3)

where the new control variable ${\varvec{u}} = [u_1,u_2]^\mathrm{T}$ is defined as

$$\begin{aligned} {\varvec{u}} = \begin{bmatrix} \text{ cos }\theta _r &{} -l_0 \ \text{ sin }\theta _r\\ \text{ sin }\theta _r &{} l_0\ \text{ cos } \theta _r \end{bmatrix}\begin{bmatrix} \nu _r\\ \omega _r. \end{bmatrix} \end{aligned}$$

(4)

Note that the determinant of matrix $\begin{bmatrix} \text{ cos }\theta _r &{} -l_0\text{ sin }\theta _r\\ \text{ sin }\theta _r &{} l_0\text{ cos }\theta _r\\ \end{bmatrix}$ is always $l_0$ no matter the value of the heading angle $\theta _r$; therefore, the matrix’s inverse always exists and it is given by $\begin{bmatrix} \text{ cos }\theta _r &{} \text{ sin }\theta _r\\ -\frac{\text{ sin }\theta _r}{l_0} &{} -\frac{\text{ cos } \theta _r}{l_0} \end{bmatrix}$. After coordinate transformation, designed control input ${\varvec{u}}$ can be directly applied to the integrator model (Sahin et al. 2007; Jiang et al. 2020). In the following content, we only consider the single-integrator model (3).

3.2 Radial basis function

Wildfire boundary can be regarded as a series of scattered points $\Gamma (t)$. With the feature of level set algorithm, fire boundary is implicitly defined as a zero level set of a smooth time-dependent function $\Phi ({\varvec{x}},t): {\mathbb {R}}^2 \rightarrow {\mathbb {R}}$, namely

$$\begin{aligned} \Gamma (t) = \{{\varvec{x}}\in {\mathbb {R}}^2 | \Phi ({\varvec{x}},t)=0 \}. \end{aligned}$$

(5)

The initialization of $\Phi$ usually uses signed distance function, given by

$$\begin{aligned} \Phi ({\varvec{x}},0) = {\left\{ \begin{array}{ll} \begin{aligned} &{}-d_{\Gamma (t)}({\varvec{x}}) &{}&{}{\varvec{x}}\ \text{ inside } \text{ the } \text{ fire } \text{ boundary } \\ &{}0 &{}&{}{\varvec{x}}\ \text{ on } \text{ the } \text{ fire } \text{ boundary } \\ &{}d_{\Gamma (t)}({\varvec{x}}) &{}&{}{\varvec{x}}\ \text{ outside } \text{ the } \text{ fire } \text{ boundary } \\ \end{aligned} \end{array}\right. }, \end{aligned}$$

(6)

where $d_{\Gamma (t)}({\varvec{x}})$ is the distance between ${\varvec{x}}$ and its nearest point on the wildfire boundary (Zhai et al. 2020).

Radial basis functions (RBFs) are a set of real-valued axial symmetric functions, namely the value only relies on the distance to the origin. To build a geometric model from multivariate scattered data, the radial basis function approach is a popular and effective tool due to its ease of implementation with a simple form and good approximation behavior. There are many types of popular radial-based functions and thin-plate splines are used in this paper to construct the level set function for wildfire boundary.

The thin-plate spline can be written as

$$\begin{aligned} g_i({\varvec{x}}) = \parallel {\varvec{x}} - {\varvec{x}}_i\parallel ^2 \text{ log }(\parallel {\varvec{x}} - {\varvec{x}}_i \parallel +1), \end{aligned}$$

(7)

where the ${\varvec{x}}_i = [x_i,y_i]^\mathrm{T} \in {\mathbb {R}}^2$ are referred to the centers of radial basis function, $\parallel \cdot \parallel$ represents the Euclidean norm operator. The points on wildfire boundary can be approximated by a N number of thin-plate splines, centered at N fixed centers, which can be written as

$$\begin{aligned} \Phi ({\varvec{x}},t) = \sum _{i=1}^{N}\lambda _i(t)g_i({\varvec{x}}) + p({\varvec{x}},t), \end{aligned}$$

(8)

where coefficients $\lambda _i(t)$ are real numbers, $p({\varvec{x}},t)$ is a first-order polynomial, determined by the selected RBF (Koushki et al. 2020) and changed with time to account for the linear and constant portion of $\Phi ({\varvec{x}},t)$ and to ensure positive definiteness of the solution (Morse et al. 2005).

The polynomial $p({\varvec{x}},t)$ is not required for positive definite RBFs; however, it is necessary for semi-positive RBFs to guarantee singularity (Buhmann 2003). The adopted thin-plate spline RBF is semi-positive; therefore, we take the polynomial part as $p({\varvec{x}},t) = c_1(t) + c_2(t)x +c_3(t)y$ for the 2D problem (Xie and Mirmehdi 2007; Wei et al. 2018). To ensure a unique solution of RBF interpolation of the level set function, the expansion coefficients in Eq. (8) must be subjected to the following orthogonality conditions (Morse et al. 2005):

$$\begin{aligned} \sum _{i=1}^{N}\lambda _i(t) = \sum _{i=1}^{N}\lambda _i(t)x_i = \sum _{i=1}^{N}\lambda _i(t)y_i = 0. \end{aligned}$$

(9)

The function (8) combined with the constraints in (9) can be rewritten in a matrix form (Xie and Mirmehdi 2007; Wei et al. 2018)

$$\begin{aligned} {\varvec{G}}\varvec{\alpha } = {\varvec{f}}, \end{aligned}$$

(10)

where ${\varvec{G}} = \begin{bmatrix} {\varvec{A}}&{}{\varvec{P}}\\ {\varvec{P}}^\mathrm{T}&{} {\varvec{0}}_{3\times 3} \end{bmatrix}$, ${\varvec{A}} = \begin{bmatrix} g_1({\varvec{x}}_1) &{} \ldots &{}g_N({\varvec{x}}_1)\\ \vdots &{} \ddots &{}\vdots \\ g_1({\varvec{x}}_N) &{} \ldots &{}g_N({\varvec{x}}_N)\\ \end{bmatrix}$, ${\varvec{P}} = \left[ \begin{matrix} 1 &{} x_1&{}y_1\\ \vdots &{} \vdots &{}\vdots \\ 1 &{} x_N&{}y_N\\ \end{matrix}\right]$, $\varvec{\alpha } = \begin{bmatrix}\lambda _1&\ldots&\lambda _N&c_1&c_2&c_3\end{bmatrix}^\mathrm{T}$ and ${\varvec{f}} =\begin{bmatrix}\Phi ({\varvec{x}}_1,t)&\ldots&\Phi ({\varvec{x}}_N,t)&0&0&0\end{bmatrix}^\mathrm{T}$.

3.3 Wildfire model

Existing wildfire propagation models can be summarized into physical models, empirical models, and quasi-empirical models (Sullivan 2009a, b; Sullivan and Andrew 2009). Wildfires involve complex combustion chemistry, heat transfer, and fluid dynamics, making many data is unmeasurable for physical models. Empirical models derive statistical correlations between observations and the wildfire, while they are restricted by the lack of physical basis in different environments. Quasi-empirical models use some forms of physical framework for statistical modeling and they often include simplified equations with less number of parameters (Karouni et al. 2014). Miyanishi (2001) proposed a semi-empirical fire spread model focusing on front velocities at the head, the rear, the flanks, and elsewhere between them (Fig. 1). The fire spread velocities primarily depend on the wind velocity.

Fig. 1

A schematic of wildfire propagation during the time interval $\Delta t$ under a wind of magnitude $v_w$ in the direction $\varphi = 0$

In this paper, homogeneous topography, vegetation, and meteorology are assumed in the propagation of wildfire boundary. Without loss of generality, a simplified Fedell’s model (Mallet et al. 2009) which is easier to tune, either via direct trials or with systemic methods for parameter estimation, is adopted in this paper. The rate of spread (RoS) at a point ${\varvec{x}}$ on the fire boundary is given by

$$\begin{aligned} U({\varvec{x}},t) = {\left\{ \begin{array}{ll} \begin{aligned} &{}\epsilon _0 + a\sqrt{\nu _w\text{ cos}^n_u\varphi }&{} \ \text{ if } \ |\varphi |\le \frac{\pi }{2}\\ &{}\epsilon _0[\alpha + (1-\alpha )|\text{ sin }\varphi |]&{} \ \text{ if } \ |\varphi |> \frac{\pi }{2}\\ \end{aligned} \end{array}\right. }, \end{aligned}$$

(11)

where $\nu _w$ is the wind velocity, $\varphi$ is the angle (measured counterclockwise) between the wind direction and the outward normal to the wildfire boundary at point ${\varvec{x}}$, $\epsilon _0$ and a are determined by fuel, and $\alpha \in [0,1]$ is the ratio between RoS at the flanks ($\varphi = \frac{\pi }{2}$) and that at the rear ($\varphi = \pi$). In this paper, we assume that the parameters of RoS are known.

In this work, level set algorithm is implemented to describe the propagation of wildfire boundary. Although $\Phi ({\varvec{x}},t)$ has no physical meaning, the evolution of its zero level set reveals the propagation of wildfire boundary. Assuming that the wildfire boundary only propagates in its normal directions with RoS $U({\varvec{x}},t)$; therefore, the evolution of function $\Phi ({\varvec{x}},t)$ is governed by the Hamilton–Jacobi equation (Alessandri et al. 2021)

$$\begin{aligned} \frac{\partial \Phi {({\varvec{x}},t)}}{\partial t} + U({\varvec{x}},t) \parallel \nabla \Phi ({\varvec{x}},t) \parallel = 0, \end{aligned}$$

(12)

where $\nabla \Phi ({\varvec{x}}) = [\frac{\partial \Phi ({\varvec{x}})}{\partial x }, \frac{\partial \Phi ({\varvec{x}})}{\partial y} ]^\mathrm{T}$ is the gradient of the level set function. Although Eq. (12) can be solved with the conventional finite difference methods proposed in Peng et al. (1999), Mallet et al. (2009), the parameterized level set function based on radial basis functions Wei et al. (2018) is adopted here to obtain smooth wildfire boundaries. Therefore, the original time-dependent Hamilton–Jacobi partial differential equation (12) is discretized into a system of coupled ODEs governing the motion of the dynamics interfaces, which can be regarded as a collocation formulation of the general method of lines (Wei et al. 2018)

$$\begin{aligned} {\varvec{G}}\frac{\mathrm{d}\varvec{\alpha }}{\mathrm{d}t} + {\varvec{B}}(\varvec{\alpha },t) = 0, \end{aligned}$$

(13)

where ${\varvec{B}}(\varvec{\alpha },t) = \begin{bmatrix} U({\varvec{x}}_1,t)|\nabla {\varvec{g}}({\varvec{x}}_1) \varvec{\alpha }|&\ldots&U({\varvec{x}}_1,t)|\nabla {\varvec{g}}({\varvec{x}}_N) \varvec{\alpha }|&0&0&0 \end{bmatrix}^\mathrm{T}$ and ${\varvec{g}}({\varvec{x}}) = \left[ g_1({\varvec{x}})\ \ldots \ g_N({\varvec{x}})\quad 1 \right. \left. x\quad y \right] ^\mathrm{T}$. The system (13) is coupled nonlinear ODEs which can be solved by Euler methods. Thus, it is more computationally efficient to evolve boundaries by updating a set of parameterized coefficients than finite difference methods.

3.4 Control objective

We assume that the velocity of UAV is much faster than that of an advancing wildfire boundary and the information of wildfire boundary is available to each UAV. All deployed UAVs can evaluate their positions and speeds and have information on nearby UAVs. According to the definition of the level set function, it can be regarded as the tracking error. We define the notion of error distance as $|\Phi ({\varvec{x}},t)|$. The goal is to design a controller for the UAV to eliminate the tracking error $\lim \nolimits _{t\rightarrow \infty } | \Phi ({\varvec{x}},t)| \rightarrow 0$, bringing thus the UAV to the wildfire boundary. Upon reaching the boundary of wildfire, the UAV should patrol along the boundary. The control objective is formally stated as follows.

Problem: For the dynamic wildfire boundary, whose propagation is modeled by the Hamilton–Jacobi equation (12), design a control law to drive the UAV, which is subject to the kinematic constraint (3), to track the propagated wildfire boundary $\{{\varvec{x}}\in {\mathbb {R}}^2 | \Phi ({\varvec{x}},t)=0 \}$ with a desired velocity $\nu _d$, eliminating the tracking error $\lim \nolimits _{t\rightarrow \infty } | \Phi ({\varvec{x}},t)| \rightarrow 0$.

4 Trajectory tracking control

In this section, we implement the principle of vector field technique to design a guidance law for the single-integrator modeled UAV to track the wildfire boundary.

The accurate expression of the wildfire boundary is vital for designing the vector field-based controller. Thus, we first use the meshless thin-plate radial basis function to approximate the wildfire boundary with an explicit function before vector field construction. To avoid the meaningless trivial solution of Eq. (10), additional points outside and inside the closed wildfire boundary are appended. A common practice to extend points $\Gamma (t)$ is generating off-surface points in their normalized external and internal normal directions (Carr et al. 2001), which are defined as

$$\begin{aligned} {\varvec{x}}_i^+ = {\varvec{x}}_i + \delta {\varvec{n}}_i,\ {\varvec{x}}_i^- = {\varvec{x}}_i - \delta {\varvec{n}}_i,\ i = 1,\ldots ,n, \end{aligned}$$

(14)

where ${\varvec{x}}_i^+$ and ${\varvec{x}}_i^-$ are off-surface points in external and internal normal directions, respectively. ${\varvec{n}}_i$ is the normalized normal vector at point ${\varvec{x}}_i$ on the wildfire boundary. $\delta$ is a small step size and its value can be selected follow the rule in Zhu and Wathen (2015), Cuomo et al. (2017) .To make the reconstructed boundary is relatively insensitively to the projection distance mentioned in Eq. (6), care must be taken when projecting off-surface points ${\varvec{x}}_i^+$ and ${\varvec{x}}_i^-$ along the normals to ensure that they do not intersect other parts of the surface. Thus, the closest points to these new constructed points are the corresponding base points generated them (Carr et al. 2001). The set of points

$$\begin{aligned} {\tilde{\Gamma }} = \Gamma \cup \Gamma ^+ \cup \Gamma ^-, \end{aligned}$$

(15)

with $\Gamma ^+ = \{{\varvec{x}}_i^+|{\varvec{x}}_i \in \Gamma (t) \}$, $\Gamma ^- = \{{\varvec{x}}_i^-|{\varvec{x}}_i \in \Gamma (t) \}$ is the extended dataset for reconstructing wildfire boundary. The value of extended dataset is given according to Eq. (6).

Applying the new dataset as the input to Eq. (10), the resulting zero level set of its solution is the approximation of the wildfire boundary. After obtaining the approximation level set function of wildfire boundary, the Hamilton–Jacobi equation will be solved with a RBF-based parameterized level set method (Wei et al. 2018) to estimate the propagated wildfire boundary.

According to the problem definition, the studied tracking problem can be decomposed into two subproblems: (i) convergence to the wildfire boundary and (ii) patrol along the wildfire boundary. Each subproblem will be solved separately by designing corresponding vector fields. Then, the composition of individual vector field for two subproblems will ensure both convergence and patrol. The total vector field can be represented in the form as

$$\begin{aligned} {\varvec{u}} = {\varvec{V}}_{\text{ con }} + {\varvec{V}}_{\text{ pat }}, \end{aligned}$$

(16)

where ${\varvec{V}}_{\text{ con }}$ produces vectors that converge to the path and ${\varvec{V}}_{\text{ pat }}$ produces vectors that patrol along the path.

The wildfire boundary is the zero level set, and the value $\Phi ({\varvec{x}},t)$ can be regarded as a signed distance function from the multirotor to the desired path, or the tracking error (Kapitanyuk et al. 2017; Kapitanyuk and Chepinsky 2013). Inspired by this property of level set function, we define a Lyapunov function $V_{p}$ as

$$\begin{aligned} V_{p} = \frac{1}{2}\Phi ({\varvec{x}},t)^2, \end{aligned}$$

(17)

which achieves a minimum at the boundary of the wildfire along the gradient descent direction $-\nabla V_{p}$. Therefore, the vector field for convergence is modeled as

$$\begin{aligned} {\varvec{V}}_{\text{ con }} = -k_n \cdot \Phi ({\varvec{x}},t) \frac{\nabla \Phi ({\varvec{x}},t)}{\parallel \nabla \Phi ({\varvec{x}},t)\parallel }, \end{aligned}$$

(18)

where $k_n$ is a scalar weight, the normal vector $\nabla \Phi ({\varvec{x}},t)$ is given by

$$\begin{aligned} \begin{aligned} \frac{\partial \Phi }{\partial x} = c_2 + \sum _{i=1}^{N}\lambda _i(x - x_i)\left( 2\text{ log }(\parallel {\varvec{x}} - {\varvec{x}}_i \parallel +1) +\frac{\parallel {\varvec{x}} - {\varvec{x}}_i \parallel }{\parallel {\varvec{x}} - {\varvec{x}}_i \parallel +1}\right) \\ \frac{\partial \Phi }{\partial y} = c_3 + \sum _{i=1}^{N}\lambda _i(y - y_i)\left( 2\text{ log }(\parallel {\varvec{x}} - {\varvec{x}}_i \parallel +1) +\frac{\parallel {\varvec{x}} - {\varvec{x}}_i \parallel }{\parallel {\varvec{x}} - {\varvec{x}}_i \parallel +1}\right) \end{aligned}. \end{aligned}$$

(19)

When the UAV reaches the boundary of wildfire, the convergent item becomes zero and the UAV is required to patrol along the wildfire boundary. The direction of the UAV’s velocity should be aligned with the orientation of the vector field. Therefore, an additional behavior constraint related to the velocity vector $\dot{{\varvec{x}}}$ and the tangent vector $\frac{({\varvec{E}} \cdot \nabla \Phi ({\varvec{x}},t))^\mathrm{T}}{\parallel \nabla \Phi ({\varvec{x}},t)\parallel }$ can be described as

$$\begin{aligned} \frac{({\varvec{E}} \cdot \nabla \Phi ({\varvec{x}},t))^\mathrm{T}}{\parallel \nabla \Phi ({\varvec{x}},t)\parallel } \dot{{\varvec{x}}} = k_t, \end{aligned}$$

(20)

where ${\varvec{E}} =\left[ \begin{matrix} 0 &{} 1\\ -1&{}0 \end{matrix}\right]$ is an orthogonal rotation matrix and determines the moving direction of the UAV along the wildfire boundary, and $k_t$ is a scalar weight. According to Eq. (20), the vector field for patrol is defined as

$$\begin{aligned} {\varvec{V}}_{\text{ pat }} = k_t \cdot \frac{ {\varvec{E}} \cdot \nabla \Phi ({\varvec{x}},t)}{\parallel \nabla \Phi ({\varvec{x}},t)\parallel }. \end{aligned}$$

(21)

Combining the two terms defined in (18) and (21), the vector field-based controller is represented as

$$\begin{aligned} {\varvec{u}} = -k_n\cdot \Phi ({\varvec{x}},t) \frac{ \nabla \Phi ({\varvec{x}},t)}{\parallel \nabla \Phi ({\varvec{x}},t)\parallel }+ k_t\cdot \frac{ {\varvec{E}} \cdot \nabla \Phi ({\varvec{x}},t)}{\parallel \nabla \Phi ({\varvec{x}},t)\parallel }. \end{aligned}$$

(22)

To ensure that the velocity of the UAV is the desired value $\nu _d$, Eq. (22) is multiplied with a coefficient $\frac{\nu _d}{\sqrt{k_n^2\Phi ({\varvec{x}},t)^2+k_t^2}}$, and thus, the controller is

$$\begin{aligned} {\varvec{u}} = -\frac{k_n\nu _d\cdot \Phi ({\varvec{x}},t)}{\sqrt{k_n^2\Phi ({\varvec{x}},t)^2+k_t^2}}\cdot \frac{ \nabla \Phi ({\varvec{x}},t)}{\parallel \nabla \Phi ({\varvec{x}},t)\parallel }+ \frac{k_t\nu _d}{\sqrt{k_n^2\Phi ({\varvec{x}},t)^2+k_t^2}}\cdot \frac{ {\varvec{E}} \cdot \nabla \Phi ({\varvec{x}},t)}{\parallel \nabla \Phi ({\varvec{x}},t)\parallel }. \end{aligned}$$

(23)

In the constructed vector field (23), the values of $k_n$ and $k_t$ do not influence the velocity magnitude of the UAV, as the magnitude of the constructed vector field is $v_d$. However, they influence the time taken by the UAV to converge to the tracked trajectory. The increase of $k_n$ or the decrease of $k_t$ will raise the portion of convergence controller term, thus increasing the UAV’s rate of convergence to the target curve.

For the case of multiple UAVs, the potential conflicts among UAVs are avoided by introducing the repelling field. We define a function $V_r$ as

$$\begin{aligned} V_{r} = \sum _{i=1}^m k_r e^{-\frac{\parallel {\varvec{x}} - {\varvec{x}}_{oi} \parallel }{l_r}}, \end{aligned}$$

(24)

where ${\varvec{x}}_{oi}$ is the position of the $i\text {th}$ obstacle, and $k_r$ and $l_r$ are positive coefficients. The repulsive vector can be defined as

$$\begin{aligned} {\varvec{V}}_{\text{ rep }} = -\nabla V_r = \sum _{i=1}^m \frac{k_r}{l_r}e^{-\frac{\parallel {\varvec{x}} - {\varvec{x}}_{oi} \parallel }{l_r}}\cdot \frac{{\varvec{x}} - {\varvec{x}}_{oi}}{\parallel {\varvec{x}} - {\varvec{x}}_{oi} \parallel }. \end{aligned}$$

(25)

Therefore, the final controller for the UAV is the combination of Eqs. (23) and (25)

$$\begin{aligned} {\varvec{u}}= & {} -\frac{k_n\nu _d\cdot \Phi ({\varvec{x}},t)}{\sqrt{k_n^2\Phi ({\varvec{x}},t)^2+k_t^2}}\cdot \frac{ \nabla \Phi ({\varvec{x}},t)}{\parallel \nabla \Phi ({\varvec{x}},t)\parallel }+ \frac{k_t\nu _d}{\sqrt{k_n^2\Phi ({\varvec{x}},t)^2+k_t^2}}\cdot \frac{ {\varvec{E}} \cdot \nabla \Phi ({\varvec{x}},t)}{\parallel \nabla \Phi ({\varvec{x}},t)\parallel } \nonumber \\&+ \sum _{i=1}^m \frac{k_r}{l_r}e^{-\frac{\parallel {\varvec{x}} - {\varvec{x}}_{oi} \parallel }{l_r}}\cdot \frac{{\varvec{x}} - {\varvec{x}}_{oi}}{\parallel {\varvec{x}} - {\varvec{x}}_{oi} \parallel }. \end{aligned}$$

(26)

The scheme of the wildfire boundary tracking mission can be explained in these steps:

• The initial map of the environment with wildfire locations is provided by sensing instruments

• Generate extended data in the internal and external normal directions of sampled waypoints on the wildfire boundary

• Apply thin-plate spline RBFs to reconstruct a level set function to approximate the wildfire boundary with the extended dataset

• Based on available information about wind, topography, and vegetation, the above level set function is then used to estimate the propagation of wildfire boundary by solving the Hamilton–Jacobi equation with RBFs

• If there are new sensing data, reconstruct a new level set function with the updated data, otherwise using the existing level set function to generate a velocity vector field for designing guidance law

• Update the trajectory for the UAV with the generated guidance law.

The constructed velocity vector field is analytical and tolerant to sensed changes, as the characteristics of radial basis functions ensure that any pop-up changes on the wildfire boundary can be modeled and updated in the approximated wildfire boundary, thus refreshing the constructed vector field. It means that once new sensing data of wildfire boundary are available from other instruments, the wildfire boundary will be renewed with the approximation of radial basis functions. Thus, the wildfire boundary will expand with the renewed boundary. Correspondingly, the reconstructed velocity field is refreshed.

5 Simulation results

In this section, numerical simulations produced in MATLAB are used to evaluate the performance of the proposed controller. UAVs are first employed to track a simulated advancing wildfire boundary, and then follows the tracking of a real dynamic wildfire boundary.

5.1 Simulated wildfire boundary application

The initial position of an arbitrarily shaped wildfire boundary is described in a polar coordinate whose origin is $[x_\mathrm{o},y_\mathrm{o}]^\mathrm{T}$

$$\begin{aligned} r = r_\mathrm{o} + k_\mathrm{s}(\text{ sin }6\beta + \text{ sin }3\beta ), \end{aligned}$$

(27)

where $k_\mathrm{s}$ is the parameter that shapes the curve, $r_\mathrm{o}$ is a constant and $\beta \in [0,2\pi ]$ is the angular coordinate. We consider the positions of waypoints on wildfire boundary in x and y directions; thus, the polar coordinate will be converted to Cartesian coordinate. In this simulation, parameters for initial boundary are defined as $x_\mathrm{o}=15$, $y_\mathrm{o}=15$, $r_\mathrm{o} = 5$, and $k_\mathrm{s}=0.8$. The parameters for RoS model (11) are set to be the same as corresponding coefficients in Mallet et al. (2009), namely $\epsilon _0 = 0.2$, $a = 0.5$, $n_u = 3$, and $\alpha = 0.5$. The magnitude of wind velocity is set as $\nu _w = 0.1 \text{ m/s }$, which is smaller than the velocity of the UAV. With the RoS model, the Hamilton–Jacobi equation (12) is numerically solved with the radial basis function-based approach (Wei et al. 2018) to generate the propagation of wildfire boundary. An example of the propagated wildfire with initial boundary described in (27) is visualized in Fig. 2.

Fig. 2

Propagated wildfire boundary by solving Hamilton–Jacobi equation

At each time the wildfire boundary is updated, radial basis functions are applied to reconstruct the boundary curve to update the velocity vector field. First, an extended dataset should be constructed in both internal and external normal directions with a small step $\delta = 0.01$ of sampled waypoints on wildfire boundary. Then follows the approximation with radial basis functions and the construction of velocity vector field. The extended dataset is shown in Fig. 3, which indicates the sampled waypoints on the initial wildfire boundary and their corresponding external and internal normal vectors. Applying RBF algorithm to the new extended dataset, we obtain an approximated surface and its intersection with the zero level set, as shown in Fig. 4. On the basis of the approximated level set function $\Phi$, a vector field for velocity (Fig. 5) is generated. Obviously, the constructed vector field can ensure convergence to the desired boundary.

Fig. 3

Sampled waypoints on initial wildfire boundary with internal (blue) and external (red) waypoint normals indicated

Fig. 4

Approximated surface whose intersection with $\Phi = 0$, i.e., the approximated surface’s zero level set is the desired wildfire boundary

Fig. 5

Vector field constructed with RBF approximation method

5.1.1 Single UAV scenario

In this section, we present the performance of the proposed algorithm with a single UAV and also compare it with the tracking strategy presented in Kumar et al. (2011). The work Kumar et al. (2011) formulated the problem of tracking the fire front as an optimization of a utility function, and designed a controller for a double integrator system on the basis of the gradient descent information of the utility function. In this paper, the square of the approximated level set function is regarded as the utility function. Its gradient descent term and a damping force relating to UAV’s current velocity are combined as the input control.

Simulation with a single UAV is set as follows. A single UAV with initial condition ${\varvec{x}} = [25\text{ m },9\text{ m}]^\mathrm{T}$ is used to track the dynamic wildfire boundary. The desired velocity is $\nu _d = 10 \text{ m/s }$ and control parameters are set as $k_n = 15$, $k_t = 10$. Simulation time and time step are $t_f = 30 \text{ s }$ and $\Delta t = 0.01 \text{ s }$. The wildfire boundary is updated every $1 \text{ s }$. Corresponding simulation results generated from the application of the presented controller in this paper and the controller in Kumar et al. (2011) are shown in Figs. 6 and 7. As shown in Fig. 6a, although the wildfire boundary propagates with time, the single UAV’s trajectory is adaptively adjusted to patrol along the dynamic wildfire boundary with the proposed guidance law. While in Fig. 6b, the controller in Kumar et al. (2011) guides the UAV to track the wildfire boundary only in the descending direction of the designed utility function. The change of tracking error at UAV’s position is shown in Fig. 7. Regardless of the algorithm used [i.e., the proposed algorithm or the comparison algorithm (Kumar et al. 2011)], the computed tracking error $|\Phi |$ has fluctuations at the critical time that the wildfire boundary is updated. The reason for these fluctuations is that when the UAV is commanded the updated wildfire boundary, it receives a command that is is a sudden change from the current boundary to the new boundary. Under the guidance of the proposed controller, it takes some time for the UAV to move across from the current tracking boundary to the updated boundary. However, as shown in Fig. 7, we can see that until such time the updated wildfire boundary is received, the proposed algorithm performs really well, rapidly decreasing the tracking error towards zero. By comparison, the proposed algorithm can ensure that the UAV can converge to and patrol along the wildfire boundary, while the controller (Kumar et al. 2011) only ensures the convergence to the boundary and the trajectory tracking in the outdoor expanding direction. In addition, if the fire boundary is not dynamic, the control law (Kumar et al. 2011) would have the issue of local minimum, where the gradient of utility function is zero. This is because the dynamic nature of the propagated wildfire boundary ensures that the local minimum changes, and that the UAV is not stuck in the local minimum for long.

Fig. 6

Tracking trajectory of a single UAV resulted from the proposed controller (a) and the controller in Kumar et al. (2011) (b); filled square represents the initial position of UAV, filled star represents the final positions of UAV, red line is the trajectory of UAV, and black dot-dash line is wildfire boundary

Fig. 7

The time profile of wildfire boundary tracking error of a single UAV with the proposed controller (a) and the controller in Kumar et al. (2011) (b), the time at which the impulse happens is the critical time for updating wildfire boundary

5.1.2 Multi-UAV scenario

For multi-UAV scenario, we also present a comparison of the proposed guidance law with the method in Kumar et al. (2011). Both algorithms use the same controller for collision avoidance. Three UAVs initially located at $[26\,\text{ m },8.5\,\text{ m}]^\mathrm{T}$, $[25.5\,\text{ m },11\,\text{ m}]^\mathrm{T}$ and $[24.5\,\text{ m },7.5\,\text{ m}]^\mathrm{T}$, are employed to track the dynamic wildfire boundary. Control parameters for collision avoidance are $k_r = 10$, $l_r = 2$, other simulation parameters are same as the single case. Simulation results with three UAVs are displayed in Figs. 8, 9, and 10. Figure 8 shows the dynamic tracking trajectories to the expanding wildfire boundary with three UAVs. The presented algorithm ensures that all UAVs converge to the boundary and then patrol along, as it shown in Fig. 8a. In Fig. 8b, the controller in Kumar et al. (2011) also ensures that all UAVs converge to the boundary, but all UAVs move in the expanding direction of the boundary after reaching the wildfire boundary. Although the three UAVs are deployed at very close positions at the start time, both controllers guarantee that the overall trend of distance between any two UAVs increases along with time, as shown in Fig. 9, indicating that no occurrence of collisions among UAVs. The tracking errors resulted from the application of the presented algorithm and the algorithm in Kumar et al. (2011) are shown in Fig. 10. The tracking error $|\Phi |$ of each UAV displayed in Fig. 10a indicates that although it has a sharp increase at the critical time that the wildfire boundary is updated, it reduces dramatically and converges to zero along with time. In comparison, the tracking error of each UAV resulted from the controller in Kumar et al. (2011) (Fig. 10b) decreases slower to zero at the beginning and it also has fluctuations arising from updating of wildfire boundary.

Fig. 8

Tracking trajectories of three UAVs resulted from the proposed controller (a) and the controller in Kumar et al. (2011) (b); filled squares represent the initial positions of UAVs, filled stars represent the final positions of UAVs, red dot-dash line is the trajectory of UAV 1, blue dash line is the trajectory of UAV 2, and pink line is the trajectory of UAV 3

Fig. 9

The time profiles of distance between any two UAVs with the proposed controller (a) and the controller in Kumar et al. (2011) (b); red dot-dash line is the distance between UAV 1 and UAV 2, blue dash line is the distance between UAV 1 and UAV 3, and pink line is the distance between UAV 2 and UAV 3

Fig. 10

The time profiles of wildfire boundary tracking errors of three UAVs with the proposed controller (a) and the controller in Kumar et al. (2011) (b); the time at which the impulse happens is the critical time for updating wildfire boundary, red dot-dash line represents UAV 1, blue dash line represents UAV 2, and the pink line represents UAV 3

5.2 Actual wildfire boundary application

We also apply the proposed algorithm to track actual wildfire boundaries $^1$ extracted using the image processing technique. An example of extracted wildfire boundaries at different time instants is shown in Fig. 11. The initial position of the single UAV is $[45\text{ m },18\text{ m}]^\mathrm{T}$; other parameters are same as previous. The resulted tracking trajectory and tracking error are displayed in Figs. 12 and 13, respectively. From Fig. 12, we can find that the single UAV can adapt its trajectory to closely follow along the wildfire boundary. Although the tracking error (Fig. 13) has fluctuations when the UAV is commanded to track a renewed wildfire boundary from its current position, this error can be reduced rapidly when the UAV moves across the current boundary to converge to and patrol along the updated wildfire boundary under the guidance of the proposed controller.

To summarize, the proposed algorithm works well no matter deploying a single UAV or multiple UAVs to track the dynamic wildfire boundary. Different from the the controller which only ensures the track of the dynamic wildfire boundary in its expanding direction, the presented controller in this paper ensures both convergence to and patrolling along the wildfire boundary with a constant velocity. In addition, the presented guidance law is robust to deal with abrupt changes in wildfire boundary, as the introduction of RBF-based approximation algorithm ensures that any change in the wildfire boundary can be modeled and reflected in the constructed vector field. Thus, the abrupt position differences between the UAV and the expanding wildfire boundary can rapidly reduce and converge to zero along with time.

Fig. 11

Extracted real wildfire boundaries at different time instants

Fig. 12

Tracking trajectory of a real dynamic wildfire boundary

Fig. 13

The time profile of wildfire tracking error of a single UAV, the time at which the impulse happens is the critical time for updating wildfire boundary

6 Conclusion

In this paper, we presented a guidance strategy for UAVs to track an arbitrarily shaped wildfire boundary using the level set approach. Radial basis functions are used to approximate the wildfire boundary and construct the vector field for tracking dynamic wildfire fronts. Numerical simulations with single and multiple UAVs demonstrate that the proposed guidance law is robust and ensures successful tracking of an expanding wildfire boundary, regardless of the shape of the boundary. In this paper, only the kinematic model of the UAV is considered, and the vector field-based control law for more realistically modeled UAV will be studied. In addition, we assumed that the wildfire boundary is closed and the open boundaries were not considered. However, the open boundaries can be artificially closed and the same method can be applied to track those open boundary segments.