Distributed cooperative Kalman filter constrained by advection–diffusion equation for mobile sensor networks

In this paper, a distributed cooperative filtering strategy for state estimation has been developed for mobile sensor networks in a spatial–temporal varying field modeled by the advection–diffusion equation. Sensors are organized into distributed cells that resemble a mesh grid covering a spatial area, and estimation of the field value and gradient information at each cell center is obtained by running a constrained cooperative Kalman filter while incorporating the sensor measurements and information from neighboring cells. Within each cell, the finite volume method is applied to discretize and approximate the advection–diffusion equation. These approximations build the weakly coupled relationships between neighboring cells and define the constraints that the cooperative Kalman filters are subjected to. With the estimated information, a gradient-based formation control law has been developed that enables the sensor network to adjust formation size by utilizing the estimated gradient information. Convergence analysis has been conducted for both the distributed constrained cooperative Kalman filter and the formation control. Simulation results with a 9-cell 12-sensor network validate the proposed distributed filtering method and control law.

1 Introduction

Natural phenomena, such as forest fires, hurricanes, and ocean eddies, are complex spatial–temporal processes that are influenced by physical parameters like wind speed, humidity, flow direction, flow speed, and temperature (Eshghi and Schmidtke, 2018; Wei et al., 2019; Chen and Dames, 2020; Park et al., 2018). By understanding these phenomena in real-time, precision countermeasures such as rerouting firefighters/AUVs to the growing edge of a forest fire, monitoring and evacuating local areas that are at risk of hurricane damage, or decreasing the search area of floating plane-crash survivors by measuring and estimating ocean eddies can be effectively deployed (Radmanesh et al., 2021; Elston et al., 2010). Partial differential equations (PDEs) have been proven to be effective in modeling many of these natural phenomena, and solving these PDEs can help us gain a better understanding of the spatial–temporal variations in these phenomena.

One typical PDE is the advection–diffusion equation, which has been applied to model a wide range of physical phenomena, including the propagation of chemicals, particles, energy, and other physical quantities inside media like water and air Ghez (2001). The advection–diffusion equation accounts for both spatial and temporal variations, making it applicable to many physical phenomena with varying parameters. Real-world applications of this equation include modeling oceanic oil spills like the Deepwater Horizon disaster in 2010 (Boufadel et al., 2022) and airborne chemical dispersion (Gao and Acar, 2016) resulting from train derailments, such as the one that occurred in Ohio in 2023. In both cases, understanding the behavior of the advection–diffusion equation can help predict the spread of pollutants and aid in the implementation of effective countermeasures.

Obtaining explicit analytic solutions for PDEs, such as the advection–diffusion equation, can be challenging in many cases. In the literature, static sensor networks have been proposed as a solution to estimation problems, where a large number of sensors are deployed over the domain of interest (Krstic and Smyshlyaev, 2008). However, static sensors have limited accuracy in areas with high spatial–temporal variations. To address this issue, mobile sensor networks have emerged as a more flexible and efficient alternative, requiring fewer sensors and enabling adaptive monitoring of areas of interest (Ge et al., 2019; Kumar et al., 2019; Wang et al., 2019).

In mobile sensor networks, measurements are collected along trajectories while the sensors move in groups within the field of interest. These mobile sensing data exhibit coupling between space and time, making it necessary to convert them into a map of spatial–temporal field estimates, instead of obtaining separate decoupled spatial and temporal maps Ge et al. (2019). To improve the accuracy of estimates, filters can be applied to combine multiple sensor measurements while accounting for measurement noise (Demetriou, 2021; Mishra et al., 2020; Hu et al., 2022).

Our previous work focused on solving the field estimation problem for mobile sensor networks in a centralized manner (Zhang and Leonard, 2010; You et al., 2022). Specifically, we derived information dynamics to model the spatial–temporal variations along the trajectory of the formation center of a cell of mobile sensors. We also proposed a cooperative Kalman filter to give estimations of the field value at formation center based on those derived information dynamics. In cases where the spatial–temporal field can be modeled by a known PDE, the PDE can be discretized as a constraint for the estimates of the cooperative Kalman filter (Wu et al., 2021; You et al., 2022). In our recent work (Zhang et al., 2022), we proposed a constrained cooperative Kalman filter and a framework to decouple the information dynamics and the PDE model.

In Zhang et al. (2023), we generalized the constrained cooperative Kalman filter to distributed mobile sensor networks with rigid formations in a spatial–temporal field modeled by the Poisson equation. In this approach, the mobile sensors are organized into cells and have limited communication capabilities restricted to intra-cell and adjacent cells. The estimated information at each cell center is shared only among neighboring cells according to the derived information dynamics.

The Kalman-consensus filter (Olfati-Saber, 2009) and its variants (Ji et al., 2017; Lian et al., 2020; Song et al., 2019; He et al., 2020) have been widely applied for distributed filtering, where a consensus term is added to the conventional Kalman filter for estimating states and reaching consensus with neighbors on the estimation. Different from these consensus-based filters, our proposed distributed filtering strategy does not focus on achieving consensus between neighboring sensors on the state estimation of targets. Instead, the information dynamics considered in this paper model the spatial–temporal variation along trajectories of cell centers and are defined based on the differences between neighbors. These differences establish the foundation of the proposed distributed filtering.

In this paper, we focus on solving the estimation problem for distributed mobile sensor networks in a field modeled by the advection–diffusion equation (Ghez, 2001). Distributed filtering over sensor networks has been studied with applications in target tracking, environmental monitoring, etc. Instead of maintaining all-to-all communication, sensors are only required to share information with their local neighbors, resulting in greater robustness and scalability compared with centralized filtering and estimation. In real-world scenarios such as underwater environments, low bandwidth, power limitations, high deployment expenses, and communication losses can make it challenging to maintain all-to-all communications among mobile sensor networks (Alfouzan, 2021; Omeke et al., 2021). This proposed distributed formulation is welcome in such real-world applications, as sensor power can be better distributed spatially. With limited communication, mobile sensors maintain a formation of organized cells with time-varying bounded relative distances while maneuvering in the field. For each cell, the field value and gradient information at the cell center are estimated using the sensor measurements from this cell.

The estimated information generated by the distributed filters is required to satisfy the advection–diffusion equation, which is used as a constraint for the estimation. Proper discretization and approximation of the continuous PDE are crucial in obtaining the estimation constraint. Finite volume methods are commonly applied to obtain PDE discretizations within each cell over the mesh grids composed of mobile sensors (Hermeline, 2000; Droniou, 2014). In the context of distributed mobile sensor networks, each node can be regarded as a grid point in a mesh grid that provides measurements of the field value at the grid point. The discretization of the advection–diffusion equation models the inclusion of sensor measurements in a single cell, with the information shared among neighboring cells. Such discretization is incorporated as constraints in the distributed cooperative Kalman filter and reflects the weakly coupled relationships among neighboring cells.

In this work, the desired relative positions between sensors are no longer constant and are determined by the estimated gradient information from the cooperative filters. In order to maintain the desired distances between sensors, a formation controller has been designed. Under the proposed formation control, the networked sensors are capable of spreading out with increased cell sizes when exploring fields with small gradients and forming a tighter size when the field has larger gradients. With adaptive formation sizes, the mobile sensors can explore larger areas in less time with minimal sacrifices to quality.

The major contributions of this paper are summarized as follows: 1) developing a distributed cooperative Kalman filter under constraints induced by the advection–diffusion equation for a mobile sensor network. 2) Applying the finite volume method for the PDE approximation. 3) Proving convergence for both the formation control and the distributed constrained cooperative Kalman filter.

The rest of the paper is organized as follows. Section 2 introduces the Problem formulation. Section 3 reviews the information dynamics and the measurement equations. Section 4 derives the approximation at each cell center using finite volume methods. Section 5 presents the distributed constrained cooperative Kalman filter for state estimation. Section 6 presents the formation control design for distributed mobile sensor networks. Section 7 provides the Convergence analysis. Simulation results are given in Section 8. Conclusion and future works follow in Section 9.

2 Problem formulation

Consider a spatial–temporal field in d-dimensional space, where $\mathrm{d}\in {\mathbb{Z}}_{+}$ and d ≥ 2. Denote the spatial domain of interest as $\mathrm{\Omega }\subseteq {\mathbb{R}}^{\text{d}}$, location as r ∈ Ω, and time as $t\in {\mathbb{R}}_{+}$.

2.1 Environmental model and mobile sensors

We assume that the field can be described by the following advection–diffusion equation:

$\frac{\partial z\left(\mathit{r},t\right)}{\partial t}=\theta \mathrm{\Delta }z\left(\mathit{r},t\right)+{\mathit{v}}^{⊺}\nabla z\left(\mathit{r},t\right),(1)$

where $z(\mathit{r},t):{\mathbb{R}}^{\text{d}}\times {\mathbb{R}}_{+}\mapsto \mathbb{R}$ is the field value that depends on both location r and time t, θ > 0 is a known constant diffusion coefficient, ∇ is the gradient operator, Δ = ∇² represents the Laplacian operator, and $\mathit{v}\in {\mathbb{R}}^{\text{d}}$ is a known constant vector representing flow velocity.

Eq. 1 has the initial condition z(r, 0) = $\mathit{r}\in \mathrm{\Omega }$ and the boundary condition z(r, t) = z_b(r, t) for r ∈ ∂Ω, where z₀(r) and z_b(r, t) are arbitrary initial condition and Dirichlet boundary condition, respectively.

We suppose a number of distributed mobile sensors are taking discrete measurements of the field z(r, t), which conforms to the advection–diffusion equation described in Eq. 1. In a limited communication environment, these mobile sensors can only share information with their local neighbors and are unable to communicate globally. In this manner, the mobile sensors form communication cells with the ability to communicate intra-cell as well as with adjacent cells. No communication occurs between cells not sharing a boundary. The communication network of the mobile sensors is illustrated in Figure 1.

FIGURE 1

FIGURE 1. Example of five cells C₁, …, C₅ with eight sensors (left), where the blue dots denote mobile sensors and the solid lines represent boundary edges of each cell. Intra-cell communication of cell C₁ (right), where the dashed lines represent communication between any two sensors from the same cell. Note that cells C₁ and C₄ do not share a boundary, and thus no communication occurs between these cells.

Assumption 2.1: The communication graph formed by the mobile sensors is predetermined and fixed, but the formation size changes while the sensors move in the field, as shown in Figure 1. The spatial domain covered by the mobile sensor network is defined as Ω_r with boundary ∂Ω_r, which both vary while the sensors move in the field.

Assumption 2.2: For each communication cell C_j, the mobile sensors maintain distance from each other so that the area covered by the cell denoted by $\mathcal{A}(C_j)$ will not shrink to zero size, i.e., $\mathcal{A}(C_j)\ne 0$ for any j = 1, …, N.

Assumption 2.3: The mobile sensors communicate with neighboring sensors, and there is no overlap between any communication cells. The intersection between two cells is the edge connecting two vertices of the communication graph.

Definition 2.4: Two communication cells are called neighboring cells if they have a shared edge on the cell boundaries, e.g., C₁, C₂ in Figure 1.

Assumption 2.5: Every communication cell has at least one neighboring cell. The sensors maintain all-to-all information sharing within each communication cell. The information sharing among different communication cells only happens between neighboring cells, and the shared information is the corresponding estimated information at cell centers.

We suppose the mobile sensors form multiple communication cells C₁, C₂, …, C_N. (The communication cells will be denoted as cells for simplicity in the rest of this paper.) Denote the ith mobile sensor from cell C_j as (i, j) and the location of this mobile sensor at the kth time step as ${\mathit{r}}_{i,j}^k$. Then, the corresponding noisy measurement $p({\mathit{r}}_{i,j}^k,k)$ taken by this sensor can be written as

$p\left({\mathit{r}}_{i,j}^k,k\right)=z\left({\mathit{r}}_{i,j}^k,k\right)+n_{i,j}^k,(2)$

where $n_{i,j}^k\in \mathbb{R}$ is an i.i.d. Gaussian noise with zero mean. Define the center location of sensors from cell C_j at the kth time step as ${\mathit{r}}_{c_j}^k$ and ${\mathit{r}}_{c_j}^k=\frac{1}{|C_j|}{\sum }_{i\in C_j}{\mathit{r}}_{i,j}^k$, where |C_j| is the number of mobile sensors belonging to cell C_j.

2.2 Formation control

We assume that the desired distance $d_{i_1,i_2}^k$ at the kth time step between any two communicating sensors i₁, i₂ from the same communication cell is non-constant and can be determined by gradient information as follows:

$d_{i_1,i_2}^k=d\left(\nabla z\left({\mathit{r}}_{i_1}^k,k\right),\nabla z\left({\mathit{r}}_{i_2}^k,k\right)\right),(3)$

where function $d(\cdot ):{\mathbb{R}}^{\text{d}}\mapsto \mathbb{R}$ describes the relationship between the desired distance and gradient information at the cell center ${\mathit{r}}_{c_j}$ at the kth time step.

The velocity of sensor i from communication cell C_j is described by

${\dot{\mathit{r}}}_{i,j}^k=\mathit{u}\left(z\left({\mathit{r}}_{i,j}^k\right),\nabla z\left({\mathit{r}}_{i,j}^k\right),{\mathit{d}}_i^k\right),(4)$

where u is a control law that relies on the measurement z(r_i), gradient ∇z(r_i) of the field function, and the desired distance vector ${\mathit{d}}_i^k={[d_{i,i^{\prime }}^k]}_{i^{\prime }\in {\mathcal{N}}_i}$ containing all desired distance information of sensor i.

2.3 Design goals

By utilizing the discrete measurements taken by mobile sensors, we would like to solve the estimation problem for distributed mobile sensor networks in a field modeled by the advection–diffusion equation. While the sensors are deployed in the field collecting measurements, the estimated information will be applied to maintain desired formations and guide them moving in the field.

The objectives of this paper can be summarized as follows: 1) estimate field and gradient information at cell centers along trajectories by incorporating discrete sensor measurements. 2) Develop a formation controller such that the sensors will move in desired varying formation sizes by incorporating the estimated gradient information while moving in the field of interest. 3) Provide detailed convergence analysis of the estimation under varying formations. 4) Apply formation control to the distributed sensor network to explore the field utilizing estimated information.

3 Preliminaries

In this section, we will review the results of information dynamics and measurement equations at each cell center. For more detailed derivation, interested readers can refer to our previous papers (Zhang et al., 2022) for information dynamics and (Zhang et al., 2023) for distributed approximation using the finite volume method.

3.1 Information dynamics

The relationship between total derivative $\frac{\mathrm{d}z}{\mathrm{d}t}$ and partial derivative $\frac{\partial z}{\partial t}$ of ${\mathit{r}}_{c_j}^k$ at the kth time step can be approximated using the chain rule as follows:

${\left.\left(\frac{\mathrm{d}z}{\mathrm{d}t}-\frac{\partial z}{\partial t}\right)\right|}_{\left({\mathit{r}}_{c_j}^k,k\right)}\approx \frac{1}{\delta t}{\left({\mathit{r}}_{c_j}^k-{\mathit{r}}_{c_j}^{k-1}\right)}^{⊺}\nabla z\left({\mathit{r}}_{c_j}^{k-1},k-1\right),(5)$

where δt is the sampling time interval. By applying the finite difference method, we can obtain $\frac{\mathrm{d}z}{\mathrm{d}t}$ and $\frac{\partial z}{\partial t}$ as

${\left.\frac{\mathrm{d}z}{\mathrm{d}t}\right|}_{\left({\mathit{r}}_{c_j}^k,k\right)}\approx \frac{z\left({\mathit{r}}_{c_j}^k,k\right)-z\left({\mathit{r}}_{c_j}^{k-1},k-1\right)}{\delta t},{\left.\frac{\partial z}{\partial t}\right|}_{\left({\mathit{r}}_{c_j}^k,k\right)}\approx \frac{z\left({\mathit{r}}_{c_j}^k,k+1\right)-z\left({\mathit{r}}_{c_j}^k,k\right)}{\delta t},(6)$

which leads to

${\left.\left(\frac{\mathrm{d}z}{\mathrm{d}t}-\frac{\partial z}{\partial t}\right)\right|}_{\left({\mathit{r}}_{c_j}^k,k\right)}\approx \frac{1}{\delta t}\left(2z\left({\mathit{r}}_{c_j}^k,k\right)-z\left({\mathit{r}}_{c_j}^{k-1},k-1\right)-z\left({\mathit{r}}_{c_j}^k,k+1\right)\right).(7)$

From Eqs 5, 7, we can build the following relationship:

$\left(2z\left({\mathit{r}}_{c_j}^k,k\right)-z\left({\mathit{r}}_{c_j}^{k-1},k-1\right)-z\left({\mathit{r}}_{c_j}^k,k+1\right)\right)={\left({\mathit{r}}_{c_j}^k-{\mathit{r}}_{c_j}^{k-1}\right)}^{⊺}\nabla z\left({\mathit{r}}_{c_j}^{k-1},k-1\right).(8)$

Since $z({\mathit{r}}_{c_j}^k,k)=z({\mathit{r}}_{c_j}^{k-1},k)+{({\mathit{r}}_{c_j}^k-{\mathit{r}}_{c_j}^{k-1})}^{⊺}\nabla z({\mathit{r}}_{c_j}^{k-1},k)$, we should have

$\begin{array}{cc}\hfill z({\mathit{r}}_{c_j}^k,k+1)=& 2z({\mathit{r}}_{c_j}^k,k)-z({\mathit{r}}_{c_j}^{k-1},k-1)\hfill \\ \hfill & -{({\mathit{r}}_{c_j}^k-{\mathit{r}}_{c_j}^{k-1})}^{⊺}\nabla z({\mathit{r}}_{c_j}^{k-1},k-1)\hfill \\ \hfill =& 2z({\mathit{r}}_{c_j}^{k-1},k)+2{({\mathit{r}}_{c_j}^k-{\mathit{r}}_{c_j}^{k-1})}^{⊺}\nabla z({\mathit{r}}_{c_j}^{k-1},k)\hfill \\ \hfill & -z({\mathit{r}}_{c_j}^{k-1},k-1)\hfill \\ \hfill & -{({\mathit{r}}_{c_j}^k-{\mathit{r}}_{c_j}^{k-1})}^{⊺}\nabla z({\mathit{r}}_{c_j}^{k-1},k-1).\hfill \end{array}(9)$

By utilizing a similar approach, we can obtain

$\begin{array}{ll}\hfill {\left.\frac{\mathrm{d}z}{\mathrm{d}t}\right|}_{\left({\mathit{r}}_{c_j}^k,k+1\right)}& \approx \frac{z\left({\mathit{r}}_{c_j}^k,k+1\right)-z\left({\mathit{r}}_{c_j}^{k-1},k\right)}{\delta t},\hfill \\ \hfill {\left.\frac{\partial z}{\partial t}\right|}_{\left({\mathit{r}}_{c_j}^k,k+1\right)}& \approx \frac{z\left({\mathit{r}}_{c_j}^k,k+1\right)-z\left({\mathit{r}}_{c_j}^k,k\right)}{\delta t},\hfill \end{array}(10)$

leading to

$\begin{array}{ll}\hfill {\left.\left(\frac{\mathrm{d}z}{\mathrm{d}t}-\frac{\partial z}{\partial t}\right)\right|}_{\left({\mathit{r}}_{c_j}^k,k+1\right)}& \approx \frac{1}{\delta t}\left(z\left({\mathit{r}}_{c_j}^k,k\right)-z\left({\mathit{r}}_{c_j}^{k-1},k\right)\right)\hfill \\ \hfill & =\frac{1}{\delta t}{\left({\mathit{r}}_{c_j}^k-{\mathit{r}}_{c_j}^{k-1}\right)}^{⊺}\nabla z\left({\mathit{r}}_{c_j}^{k-1},k\right),\hfill \end{array}(11)$

which means that

$z\left({\mathit{r}}_{c_j}^k,k\right)-z\left({\mathit{r}}_{c_j}^{k-1},k\right)={\left({\mathit{r}}_{c_j}^k-{\mathit{r}}_{c_j}^{k-1}\right)}^{⊺}\nabla z\left({\mathit{r}}_{c_j}^{k-1},k\right).(12)$

We define the state variable $\mathit{x}(j,k)=\left[z({\mathit{r}}_{c_j}^{k-1},k-1),\nabla z({\mathit{r}}_{c_j}^{k-1},k-1),\right.{\left.z({\mathit{r}}_{c_j}^{k-1},k),\nabla z({\mathit{r}}_{c_j}^{k-1},k)\right]}^{⊺}$, and we can get the following state equation based on Eqs 9, 12:

$\mathit{x}\left(j,k+1\right)=\mathit{A}\left(j,k\right)\mathit{x}\left(j,k\right)+\mathit{U}\left(j,k\right)+\mathit{e}\left(j,k\right),(13)$

where $\mathit{A}(j,k)\triangleq \left[\begin{array}{cccc}0& 0& 1& {({\mathit{r}}_{c_j}^k-{\mathit{r}}_{c_j}^{k-1})}^{⊺}\\ 0& {\mathit{I}}_{\text{d}\times \mathrm{d}}& 0& 0\\ -1& -{({\mathit{r}}_{c_j}^k-{\mathit{r}}_{c_j}^{k-1})}^{⊺}& 2& 2{({\mathit{r}}_{c_j}^k-{\mathit{r}}_{c_j}^{k-1})}^{⊺}\\ 0& 0& 0& I\end{array}\right]$, $\mathit{U}(j,k)\triangleq \left[\begin{array}{c}0\\ {\nabla }^{2}z({\mathit{r}}_{c_j}^{k-1},k-1)({\mathit{r}}_{c_j}^k-{\mathit{r}}_{c_j}^{k-1})\\ 0\\ {\nabla }^{2}z({\mathit{r}}_{c_j}^{k-1},k)({\mathit{r}}_{c_j}^k-{\mathit{r}}_{c_j}^{k-1})\end{array}\right]$, and e(j, k) is the noise term.

3.2 Measurement equation

The field concentration can be locally approximated by the Taylor series up to the second order as

$\begin{array}{cc}\hfill z\left({\mathit{r}}_{i,j}^{k-1},k-1\right)\approx & z\left({\mathit{r}}_{c_j}^{k-1},k-1\right)+{\left({\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{c_j}^{k-1}\right)}^{⊺}\nabla z\left({\mathit{r}}_{c_j}^{k-1},k-1\right)\hfill \\ \hfill & +\frac{1}{2}{\left({\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{c_j}^{k-1}\right)}^{⊺}\mathit{H}\left({\mathit{r}}_{c_j}^{k-1},k-1\right)\left({\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{c_j}^{k-1}\right),\hfill \\ \hfill z\left({\mathit{r}}_{i,j}^k,k\right)\approx & z\left({\mathit{r}}_{c_j}^{k-1},k\right)+{\left({\mathit{r}}_{i,j}^k-{\mathit{r}}_{c_j}^{k-1}\right)}^{⊺}\nabla z\left({\mathit{r}}_{c_j}^{k-1},k\right)\hfill \\ \hfill & +\frac{1}{2}{\left({\mathit{r}}_{i,j}^k-{\mathit{r}}_{c_j}^{k-1}\right)}^{⊺}\mathit{H}\left({\mathit{r}}_{c_j}^{k-1},k-1\right)\left({\mathit{r}}_{i,j}^k-{\mathit{r}}_{c_j}^{k-1}\right).\hfill \end{array}(14)$

Let $\mathit{Z}(j,k)={[z({\mathit{r}}_{1,j}^{k-1},k-1),\dots ,z({\mathit{r}}_{|C_j|,j}^{k-1},k-1),z({\mathit{r}}_{1,j}^k,k),\dots ,z({\mathit{r}}_{|C_j|,j}^k,k)]}^{⊺}$ be the vector of true field values. Define

$\begin{array}{ll}\hfill \mathit{C}\left(j,k\right)& \triangleq \left[\begin{array}{cccc}1& {\left({\mathit{r}}_{1,j}^{k-1}-{\mathit{r}}_{c_j}^{k-1}\right)}^{⊺}& 0& 0\\ ⋮& ⋮& ⋮& ⋮\\ 1& {\left({\mathit{r}}_{|C_j|,j}^{k-1}-{\mathit{r}}_{c_j}^{k-1}\right)}^{⊺}& 0& 0\\ 0& 0& 1& {\left({\mathit{r}}_{1,j}^k-{\mathit{r}}_{c_j}^{k-1}\right)}^{⊺}\\ ⋮& ⋮& ⋮& ⋮\\ 0& 0& 1& {\left({\mathit{r}}_{|C_j|,j}^k-{\mathit{r}}_{c_j}^{k-1}\right)}^{⊺}\end{array}\right],\hfill \\ \mathit{D}\left(j,k\right)\hfill & \triangleq \left[\begin{array}{c}\frac{1}{2}{\left[\left({\mathit{r}}_{1,j}^{k-1}-{\mathit{r}}_{c_j}^{k-1}\right)\otimes \left({\mathit{r}}_{1,j}^{k-1}-{\mathit{r}}_{c_j}^{k-1}\right)\right]}^{⊺}\\ ⋮\\ \frac{1}{2}{\left[\left({\mathit{r}}_{|C_j|,j}^{k-1}-{\mathit{r}}_{c_j}^{k-1}\right)\otimes \left({\mathit{r}}_{|C_j|,j}^{k-1}-{\mathit{r}}_{c_j}^{k-1}\right)\right]}^{⊺}\\ \frac{1}{2}{\left[\left({\mathit{r}}_{1,j}^k-{\mathit{r}}_{c_j}^{k-1}\right)\otimes \left({\mathit{r}}_{1,j}^k-{\mathit{r}}_{c_j}^{k-1}\right)\right]}^{⊺}\\ ⋮\\ \frac{1}{2}{\left[\left({\mathit{r}}_{|C_j|,j}^k-{\mathit{r}}_{c_j}^{k-1}\right)\otimes \left({\mathit{r}}_{|C_j|,j}^k-{\mathit{r}}_{c_j}^{k-1}\right)\right]}^{⊺}\end{array}\right],\hfill \end{array}(15)$

where ⊗ is the Kronecker product. The Taylor expansions (Eq. 14) for all sensors near ${\mathit{r}}_{c_j}^{k-1}$ can be rewritten in a vector form as

$\mathit{Z}\left(j,k\right)=\mathit{C}\left(j,k\right)\mathit{x}\left(j,k\right)+\mathit{D}\left(j,k\right)\mathit{H}\left(j,k\right),(16)$

where H(j, k) is a column vector obtained by rearranging elements of the Hessian $\mathit{H}({\mathit{r}}_{c_j}^{k-1},k-1)$.

Supposing $\widehat{\mathit{H}}(j,k)$ represents the estimate of the vector form Hessian H(j, k) at ${\mathit{r}}_{c_j}^k$, Eq. 2 can be remodeled as

$\mathit{P}\left(j,k\right)=\mathit{C}\left(j,k\right)\mathit{x}\left(j,k\right)+\mathit{D}\left(j,k\right)\widehat{\mathit{H}}\left(j,k\right)+\mathit{D}\left(j,k\right)\mathit{\epsilon }\left(j,k\right)+\mathbf{n}\left(j,k\right),(17)$

where $\mathit{P}(j,k)=\left[p({\mathit{r}}_{1,j}^{k-1},k-1),\dots ,p({\mathit{r}}_{|C_j|,j}^{k-1},k-1),p({\mathit{r}}_{1,j}^k,k),\dots ,\right.{\left.p({\mathit{r}}_{|C_j|,j}^k,k)\right]}^{⊺}$ is the measurement vector, ɛ(j, k) represents the error in Hessian estimation, and n(j, k) is the vector of Gaussian noise n_i in Eq. 2.

If D(j,k)^⊺D(j, k) is invertible, then the least mean square method can be applied to approximate the Hessian vector as follows (You et al., 2022):

$\widehat{\mathit{H}}\left(j,k\right)={\left(\mathit{D}{\left(j,k\right)}^T\mathit{D}\left(j,k\right)\right)}^{-1}\mathit{D}{\left(j,k\right)}^T\left[\mathit{P}\left(j,k\right)-\mathit{C}\left(j,k\right)\mathit{x}\left(j,k\right)\right].(18)$

According to the definition of D(j, k) in Eq. 15, at least d² linearly independent vectors of $({\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{c_j}^{k-1}),({\mathit{r}}_{i,j}^k-{\mathit{r}}_{c_j}^{k-1})$ are needed for invertible D(j,k)^⊺D(j, k), which means 2|C_j|≥d². In the case of d = 2, |C_j| should satisfy that $|C_j|\ge \frac{{\mathrm{d}}^2}{2}=2$. Since |C_j| = 2 means two sensors in one cell leading to $\mathcal{A}(C_j)=0$, which violates $\mathcal{A}(C_j)\ne 0$, at least three sensors are needed to form one cell.

If D(j,k)^⊺D(j, k) is not invertible, then a regularization term ϵ₀I will be added as follows:

$\widehat{\mathit{H}}\left(j,k\right)={\left(\mathit{D}{\left(j,k\right)}^T\mathit{D}\left(j,k\right)+{ϵ}_0\mathit{I}\right)}^{-1}\mathit{D}{\left(j,k\right)}^T\left[\mathit{P}\left(j,k\right)-\mathit{C}\left(j,k\right)\mathit{x}\left(j,k\right)\right],(19)$

where ϵ₀ > 0 is a small positive scalar and $\mathit{I}\in {\mathbb{R}}^{\text{d}\times \mathrm{d}}$ is the identity matrix.

4 Approximation using the finite volume method

Since the measurements taken by the mobile sensors in Eq. 2 are discrete, the continuous advection–diffusion equation in Eq. 1 should be discretized properly. The mobile sensors have already been organized into different cells, and the finite volume method can be applied for PDE discretization within each cell.

The advection–diffusion equation in Eq. 1 can be discretized at ${\mathit{r}}_{c_j}^{k-1}$ at the (k − 1)-th time step as follows:

$\frac{z\left({\mathit{r}}_{c_j}^{k-1},k\right)-z\left({\mathit{r}}_{c_j}^{k-1},k-1\right)}{\mathrm{\Delta }t}=\theta \mathrm{\Delta }z\left.\left({\mathit{r}}_{c_j}^{k-1},k-1\right)+{\mathit{v}}^{⊺}\nabla z\left({\mathit{r}}_{c_j}^{k-1},k-1\right)\right).(20)$

As we can observe that $z({\mathit{r}}_{c_j}^{k-1},k),z({\mathit{r}}_{c_j}^{k-1},k-1),\nabla z({\mathit{r}}_{c_j}^{k-1},k-1)$ are elements in the state variable x(j, k), the only term left is the Laplacian $\mathrm{\Delta }z({\mathit{r}}_{c_j}^{k-1},k-1)$. We consider using the finite volume method to obtain the Laplacian approximation such that the discretized advection–diffusion equation can be utilized as a constraint for the distributed Kalman filter.

For each formation of sensors belonging to the same communication cell, the finite volume method in Hermeline (2000) will be applied to generate the approximation. We will present the general results covering both non-boundary edge and boundary edge cases, as illustrated in Figure 2A. Since all the values are at the same kth time step, we will drop the time index k for simplicity in this approximation.

Consider a shared edge s by C_j and neighboring cell C_j′. Denote the cell centers for C_j, C_j′ as c_j, c_j′, respectively. Denote the two vertices of edge s as i, i′. As shown in Figure 2B, denote ν_jj′ as the unit outward normal vector on edge s connecting i and i′, with τ_jj′ as the corresponding unit counterclockwise tangent vector, and denote ν_ii′ as the unit outward normal vector on edge connecting c_j and c_j′, with τ_ii′ as the corresponding unit counterclockwise tangent vector. Define θ_s to be the angle between ν_jj′ and τ_ii′. Suppose the edge connecting c_j and c_j′ intersects with edge s at point a_s. Denote the length of the edge connecting c_j and a_s as d. Then, the three vectors dτ_ii′, d cos θ_sν_jj′, d sin θ_sτ_jj′ form a right triangle, as shown by the green shaded area in Figure 2B. This leads to the following relationship:

${\mathit{\nu }}_{j{j}^{\prime }}=-\tan {\theta }_s{\mathit{\tau }}_{j{j}^{\prime }}+\frac{1}{\cos {\theta }_s}{\mathit{\tau }}_{i{i}^{\prime }}.(21)$

If there is no such neighboring cell C_j′, then the middle point of edge s can be treated as c_j′, as illustrated in Figure 2C. By taking integrals of Δz over cell C_j, we can obtain $\int {\int }_{C_j}\mathrm{\Delta }z={\sum }_{s\in \partial C_j}{\int }_{S_{i{i}^{\prime }}}\nabla z\cdot {\mathit{\nu }}_{j{j}^{\prime }}$. Since $\mathcal{A}(C_j)\ne 0$, substituting ν_jj′ with Eq. 21 leads to21 leads to

$\begin{array}{cc}\hfill \mathrm{\Delta }z\approx \frac{1}{\mathcal{A}\left(C_j\right)}{\iint }_{C_j}\mathrm{\Delta }z=& \frac{1}{\mathcal{A}\left(C_j\right)}\sum _{s\in \partial C_j}{\int }_{S_{i{i}^{\prime }}}\nabla z\cdot {\mathit{\nu }}_{j{j}^{\prime }}\hfill \\ \hfill =& \frac{1}{\mathcal{A}\left(C_j\right)}\sum _{s\in \partial C_j}\left[-\tan {\theta }_s{\int }_{S_{i{i}^{\prime }}}\nabla z\cdot {\mathit{\tau }}_{j{j}^{\prime }}\right.\hfill \\ \hfill & \left.+\frac{1}{\cos {\theta }_s}{\int }_{S_{i{i}^{\prime }}}\nabla z\cdot {\mathit{\tau }}_{i{i}^{\prime }}\right].\hfill \end{array}(22)$

By the finite difference method, we know that

$\begin{array}{cc}\hfill {\int }_{S_{i{i}^{\prime }}}\nabla z\cdot {\mathit{\tau }}_{j{j}^{\prime }}& \approx z_i-z_{i^{\prime }},\hfill \\ \hfill {\int }_{S_{i{i}^{\prime }}}\nabla z\cdot {\mathit{\tau }}_{i{i}^{\prime }}& \approx \left\{\begin{array}{cc}\hfill & \frac{\Vert {\mathit{r}}_{i,j}-{\mathit{r}}_{i^{\prime },j}{\Vert }_2}{\Vert {\mathit{r}}_{c_j}-{\mathit{r}}_{c_{j^{\prime }}}{\Vert }_2}\left(z_{c_{j^{\prime }}}-z_{c_j}\right),\text{\hspace{0.17em}if\hspace{0.17em}}{S}_{i{i}^{\prime }}\notin \partial {\mathrm{\Omega }}_{\mathit{r}},\hfill \\ \hfill & \frac{\Vert {\mathit{r}}_{i,j}-{\mathit{r}}_{i^{\prime },j}{\Vert }_2}{\Vert {\mathit{r}}_{c_j}-{\mathit{r}}_{c_{j^{\prime }}}{\Vert }_2}\left(\frac{z_i+z_{i^{\prime }}}{2}-z_{c_j}\right),\text{\hspace{0.17em}if\hspace{0.17em}}{S}_{i{i}^{\prime }}\in \partial {\mathrm{\Omega }}_{\mathit{r}},\hfill \end{array}\right.\hfill \end{array}(23)$

where S_ii′ is the edge connecting mobile sensors i, i′, and ∂Ω_r is the boundary of the area covered the entire mobile sensor network. This will lead to

$\begin{array}{cc}\hfill \mathrm{\Delta }z=& -\frac{1}{\mathcal{A}\left(C_j\right)}\sum _{s\in \partial C_j}\tan {\theta }_s\left(z_i-z_{i^{\prime }}\right)\hfill \\ \hfill & +\frac{1}{\mathcal{A}\left(C_j\right)}\sum _{s\in \left(\partial C_j⧵\partial {\mathrm{\Omega }}_{\mathit{r}}\right)}\frac{1}{\cos {\theta }_s}\frac{\Vert {\mathit{r}}_{i,j}-{\mathit{r}}_{i^{\prime },j}{\Vert }_2}{\Vert {\mathit{r}}_{c_j}-{\mathit{r}}_{c_{j^{\prime }}}{\Vert }_2}\left(z_{c_{j^{\prime }}}-z_{c_j}\right)\hfill \\ \hfill & +\frac{1}{\mathcal{A}\left(C_j\right)}\sum _{s\in \left(\partial C_j\cap \partial {\mathrm{\Omega }}_{\mathit{r}}\right)}\frac{1}{\cos {\theta }_s}\frac{\Vert {\mathit{r}}_{i,j}-{\mathit{r}}_{i^{\prime },j}{\Vert }_2}{\Vert {\mathit{r}}_{c_j}-{\mathit{r}}_{c_{j^{\prime }}}{\Vert }_2}\hfill \\ \hfill & \times \left(\frac{z_i+z_{i^{\prime }}}{2}-z_{c_j}\right).\hfill \end{array}(24)$

Note that if cell C_j has no boundary edges, then ∂C_j⧵∂Ω_r = ∂C_j and ∂C_m ∩ ∂Ω_r = ∅. Thus, the approximation in Eq. 24 covers the cases with and without boundary edges. Substituting Δz in Eq. 24 back to the discretized advection–diffusion (Eq. 20) will lead to

$\begin{array}{cc}\hfill & \frac{z\left({\mathit{r}}_{c_j}^{k-1},k\right)-z\left({\mathit{r}}_{c_j}^{k-1},k-1\right)}{\mathrm{\Delta }t}={\mathit{v}}^{⊺}\nabla z\left({\mathit{r}}_{c_j}^{k-1},k-1\right)\hfill \\ \hfill & -\frac{\theta }{\mathcal{A}\left(C_j\right)}\sum _{s\in \partial C_j}\tan {\theta }_s\left[z\left({\mathit{r}}_{i,j}^{k-1},k-1\right)-z\left({\mathit{r}}_{i^{\prime },j}^{k-1},k-1\right)\right]\hfill \\ \hfill & +\frac{\theta }{\mathcal{A}\left(C_j\right)}\sum _{s\in \left(\partial C_j⧵\partial {\mathrm{\Omega }}_{\mathit{r}}\right)}\frac{1}{\cos {\theta }_s}\frac{\Vert {\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{i^{\prime },j}^{k-1}{\Vert }_2}{\Vert {\mathit{r}}_{c_j}^{k-1}-{\mathit{r}}_{c_{j^{\prime }}}^{k-1}{\Vert }_2}\hfill \\ \hfill & \times \left[z\left({\mathit{r}}_{c_{j^{\prime }}}^{k-1},k-1\right)-z\left({\mathit{r}}_{c_j}^{k-1},k-1\right)\right]\hfill \\ \hfill & +\frac{\theta }{\mathcal{A}\left(C_j\right)}\sum _{s\in \left(\partial C_j\cap \partial {\mathrm{\Omega }}_{\mathit{r}}\right)}\frac{1}{\cos {\theta }_s}\frac{\Vert {\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{i^{\prime },j}^{k-1}{\Vert }_2}{\Vert {\mathit{r}}_{c_j}^{k-1}-{\mathit{r}}_{c_{j^{\prime }}}^{k-1}{\Vert }_2}\hfill \\ \hfill & \times \left[\frac{z\left({\mathit{r}}_{i,j}^{k-1},k-1\right)+z\left({\mathit{r}}_{i^{\prime },j}^{k-1},k-1\right)}{2}-z\left({\mathit{r}}_{c_j}^{k-1},k-1\right)\right].\hfill \end{array}(25)$

This can be rewritten in linear form with respect to state variable $\mathit{x}(j,k)={[z({\mathit{r}}_{c_j}^{k-1},k-1),\nabla z({\mathit{r}}_{c_j}^{k-1},k-1),z({\mathit{r}}_{c_j}^{k-1},k),\nabla z({\mathit{r}}_{c_j}^{k-1},k)]}^{⊺}$

${\mathit{G}}^{⊺}\left(j,k\right)\mathit{x}\left(j,k\right)=\mathit{g}\left(j,k\right),(26)$

where $\mathit{G}(j,k)\triangleq {\left[\begin{array}{cccc}\hfill \left(-1+\frac{\theta \mathrm{\Delta }t}{\mathcal{A}(C_j)}{\sum }_{s\in \partial C_j}\frac{1}{\cos {\theta }_s}\frac{\Vert {\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{i^{\prime },j}^{k-1}{\Vert }_2}{\Vert {\mathit{r}}_{c_j}^{k-1}-{\mathit{r}}_{c_{j^{\prime }}}^{k-1}{\Vert }_2}\right)\hfill & \hfill -{(\mathit{v}\mathrm{\Delta }t)}^{⊺}\hfill & \hfill 1\hfill & \hfill 0_{1\times d}\hfill \end{array}\right]}^{⊺}$, and

$\begin{array}{cc}\hfill \mathit{g}(j,k)\triangleq & -\frac{\theta \mathrm{\Delta }t}{\mathcal{A}(C_j)}\sum _{s\in \partial C_j}\tan {\theta }_s\left[z({\mathit{r}}_{i,j}^{k-1},k-1)-z({\mathit{r}}_{i^{\prime },j}^{k-1},k-1)\right]\hfill \\ \hfill & +\frac{\theta \mathrm{\Delta }t}{\mathcal{A}(C_j)}\sum _{s\in (\partial C_j⧵\partial {\mathrm{\Omega }}_{\mathit{r}})}\frac{1}{\cos {\theta }_s}\frac{\Vert {\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{i^{\prime },j}^{k-1}{\Vert }_2}{\Vert {\mathit{r}}_{c_j}^{k-1}-{\mathit{r}}_{c_{j^{\prime }}}^{k-1}{\Vert }_2}z({\mathit{r}}_{c_{j^{\prime }}}^{k-1},k-1)\hfill \\ \hfill & +\frac{\theta \mathrm{\Delta }t}{\mathcal{A}(C_j)}\sum _{s\in (\partial C_j\cap \partial {\mathrm{\Omega }}_{\mathit{r}})}\frac{1}{\cos {\theta }_s}\frac{\Vert {\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{i^{\prime },j}^{k-1}{\Vert }_2}{\Vert {\mathit{r}}_{c_j}^{k-1}-{\mathit{r}}_{c_{j^{\prime }}}^{k-1}{\Vert }_2}\frac{z({\mathit{r}}_{i,j}^{k-1},k-1)+z({\mathit{r}}_{i^{\prime },j}^{k-1},k-1)}{2}.\hfill \end{array}(27)$

FIGURE 2

FIGURE 2. Illustration of ν_jj′, τ_jj′, ν_ii′, τ_ii′, θ_s with given cell centers c_j, c_j′ and given edge s connecting vertices i, i′. Orange triangles represent cell centers. (A) Cell 2 with non-boundary and boundary edges. (B) Illustration of non-boundary edge. (C) Illustration of boundary edge.

5 State estimation using constrained cooperative Kalman filter

In this section, we will present the distributed constrained cooperative Kalman filter, which uses local sensor measurements to estimate information at each cell center and treats the PDE model as a constraint for the estimated information.

Assumption 5.1: We assume that the noises e(j, k), ɛ(j, k), and n(j, k) are i.i.d. Gaussian noise with zero mean and with constant covariance matrix for all j, i.e., E[e(j, k)e^⊺(j, k)] = W, E[ɛ(j, k)ɛ^⊺(j, k)] = Q, and E[n(j, k)n^⊺(j, k)] = R_n.The constrained cooperative Kalman filter can be constructed using six steps:

(1) One-step state prediction

${\widehat{\mathit{x}}}^{-}\left(j,k\right)=\mathit{A}\left(j,k-1\right){\tilde{\mathit{x}}}^{+}\left(j,k-1\right)+\mathit{U}\left(j,k-1\right),(28)$

where ${\tilde{\mathit{x}}}^{+}(j,k-1)$ is the constrained state estimate from the previous time step and ${\widehat{\mathit{x}}}^{-}(j,k)$ is the one-step state prediction.

(2) Error covariance of ${\widehat{\mathit{x}}}^{-}(j,k)$

${\mathit{R}}^{-}\left(j,k\right)=\mathit{A}\left(j,k-1\right){\mathit{R}}^{+}\left(j,k-1\right){\mathit{A}}^{⊺}\left(j,k-1\right)+\mathit{W}.(29)$

(3) Optimal gain

$\begin{array}{ll}\hfill \mathit{K}\left(j,k\right)& ={\mathit{R}}^{-}\left(j,k\right){\mathit{C}}^{⊺}\left(j,k\right)\left[\mathit{C}\left(j,k\right){\mathit{R}}^{-}\left(j,k\right){\mathit{C}}^{⊺}\left(j,k\right)\right.\hfill \\ \hfill & +{\left.\mathit{D}\left(j,k\right)\mathit{Q}{\mathit{D}}^{⊺}\left(j,k\right)+{\mathit{R}}_n\right]}^{-1}.\hfill \end{array}(30)$

The Hessian estimate can be approximated using the one-step state prediction as

$\widehat{\mathit{H}}\left(j,k\right)={\left(\mathit{D}{\left(j,k\right)}^T\mathit{D}\left(j,k\right)\right)}^{-1}\mathit{D}{\left(j,k\right)}^T\left(\mathit{P}\left(j,k\right)-\mathit{C}\left(j,k\right){\widehat{\mathit{x}}}^{-}\left(j,k\right)\right).$

(4) Updated unconstrained state estimate

${\widehat{\mathit{x}}}^{+}\left(j,k\right)={\widehat{\mathit{x}}}^{-}\left(j,k\right)+\mathit{K}\left(j,k\right)\left(\mathit{P}\left(j,k\right)-\mathit{C}\left(j,k\right){\widehat{\mathit{x}}}^{-}\left(j,k\right)-\mathit{D}\left(j,k\right)\widehat{\mathit{H}}\left(j,k\right)\right).(31)$

(5) Error covariance of ${\widehat{\mathit{x}}}^{+}(j,k)$

${\left({\mathit{R}}^{+}\left(j,k\right)\right)}^{-1}={\left({\mathit{R}}^{-}\left(j,k\right)\right)}^{-1}+{\mathit{C}}^{⊺}\left(j,k\right){\left[\mathit{D}\left(j,k\right)Q{\mathit{D}}^{⊺}\left(j,k\right)+{\mathit{R}}_n\right]}^{-1}\mathit{C}\left(j,k\right).(32)$

After running the five steps of the unconstrained filter, each cell obtains an unconstrained state estimate ${\widehat{\mathit{x}}}^{+}(j,k)$, and it will share this information with all neighboring cells to update the constrained estimate.Since each cell will not shrink to zero size and the flow velocity vector v ≠ 0, the following relationship can be always satisfied:

$\begin{array}{ll}\hfill {\mathit{G}}^{⊺}\left(j,k\right)\mathit{G}\left(j,k\right)& ={\left(-1+\frac{\theta \mathrm{\Delta }t}{\mathcal{A}\left(C_j\right)}\sum _{s\in \partial C_j}\frac{1}{\cos {\theta }_s}\frac{\Vert {\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{i^{\prime },j}^{k-1}{\Vert }_2}{\Vert {\mathit{r}}_{c_j}^{k-1}-{\mathit{r}}_{c_{j^{\prime }}}^{k-1}{\Vert }_2}\right)}^2\hfill \\ \hfill & +{\left(\mathrm{\Delta }t\right)}^2\Vert \mathit{v}{\Vert }^2+1\ne 0,\forall j,k.\hfill \end{array}(33)$

(6) Updated constrained state estimate

$\begin{array}{ll}\hfill {\tilde{\mathit{x}}}^{+}\left(j,k\right)& ={\widehat{\mathit{x}}}^{+}\left(j,k\right)-\mathit{G}\left(j,k\right){\left({\mathit{G}}^{⊺}\left(j,k\right)\mathit{G}\left(j,k\right)\right)}^{-1}\hfill \\ \hfill & \times \left({\mathit{G}}^{⊺}\left(j,k\right){\widehat{\mathit{x}}}^{+}\left(j,k\right)-\widehat{\mathit{g}}\left(j,k\right)\right),\hfill \end{array}(34)$

where

$\begin{array}{cc}\hfill \widehat{\mathit{g}}\left(j,k\right)\triangleq & -\frac{\theta \mathrm{\Delta }t}{\mathcal{A}\left(C_j\right)}\sum _{s\in \partial C_j}\tan {\theta }_s\left[p\left({\mathit{r}}_{i,j}^{k-1},k-1\right)-p\left({\mathit{r}}_{i^{\prime },j}^{k-1},k-1\right)\right]\hfill \\ \hfill & +\frac{\theta \mathrm{\Delta }t}{\mathcal{A}\left(C_j\right)}\sum _{s\in \left(\partial C_j⧵\partial {\mathrm{\Omega }}_{\mathit{r}}\right)}\frac{1}{\cos {\theta }_s}\frac{\Vert {\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{i^{\prime },j}^{k-1}{\Vert }_2}{\Vert {\mathit{r}}_{c_j}^{k-1}-{\mathit{r}}_{c_{j^{\prime }}}^{k-1}{\Vert }_2}{\widehat{z}}^{+}\left({\mathit{r}}_{c_{j^{\prime }}}^{k-1},k-1\right)\hfill \\ \hfill & +\frac{\theta \mathrm{\Delta }t}{\mathcal{A}\left(C_j\right)}\sum _{s\in \left(\partial C_j\cap \partial {\mathrm{\Omega }}_{\mathit{r}}\right)}\frac{1}{\cos {\theta }_s}\frac{\Vert {\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{i^{\prime },j}^{k-1}{\Vert }_2}{\Vert {\mathit{r}}_{c_j}^{k-1}-{\mathit{r}}_{c_{j^{\prime }}}^{k-1}{\Vert }_2}\hfill \\ \hfill & \times \frac{p\left({\mathit{r}}_{i,j}^{k-1},k-1\right)+p\left({\mathit{r}}_{i^{\prime },j}^{k-1},k-1\right)}{2}.\hfill \end{array}(35)$

This is an approximation of g(j, k) with ${\widehat{z}}^{+}({\mathit{r}}_{c_{j^{\prime }}}^{k-1},k-1)=\left[\begin{array}{cccc}\hfill 1\hfill & \hfill 0_{1\times \mathrm{d}}\hfill & \hfill 0\hfill & \hfill 0_{1\times \mathrm{d}}\hfill \end{array}\right]{\widehat{\mathit{x}}}^{+}(j^{\prime },k)$, and it uses measurements $p({\mathit{r}}_{i,j}^{k-1},k-1),p({\mathit{r}}_{i^{\prime },j}^{k-1},k-1)$ as approximations for $z({\mathit{r}}_{i,j}^{k-1},k-1),z({\mathit{r}}_{i^{\prime },j}^{k-1},k-1)$. According to the definition$\mathit{G}(j,k)={\left[\begin{array}{cccc}\hfill \left(-1+\frac{\theta \mathrm{\Delta }t}{\mathcal{A}(C_j)}{\sum }_{s\in \partial C_j}\frac{1}{\cos {\theta }_s}\frac{\Vert {\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{i^{\prime },j}^{k-1}{\Vert }_2}{\Vert {\mathit{r}}_{c_j}^{k-1}-{\mathit{r}}_{c_{j^{\prime }}}^{k-1}{\Vert }_2}\right)\hfill & \hfill -{(\mathit{v}\mathrm{\Delta }t)}^{⊺}\hfill & \hfill 1\hfill & \hfill 0_{1\times \mathrm{d}}\hfill \end{array}\right]}^{⊺}$, the updated constrained estimate in Eq. 34 updates the estimates of $z({\mathit{r}}_{c_j}^{k-1},k-1),\nabla z({\mathit{r}}_{c_j}^{k-1},k-1),z({\mathit{r}}_{c_j}^{k-1},k)$ and does not affect the gradient estimate $\nabla z({\mathit{r}}_{c_j}^{k-1},k)$ obtained by Eq. 31.

Remark 5.2: For the six steps of the distributed constrained cooperative Kalman filter described in Eqs 28–34, the first five steps in Eqs 28–32 use the local information in one cell, while the last step in Eq. 34 uses estimated information from neighboring cells.

Remark 5.3: The distributed constrained cooperative Kalman filter proposed in Zhang et al. (2023) updates the field value estimate without affecting the gradient estimation, while this work updates the whole state estimate using the constraint information. In addition, the formation cells can be arbitrary shapes and arbitrary sizes as long as the areas are non-zero indicated by Eq. 33.Define a combined state vector x(k) to include all distributed state vectors as

$\mathit{x}\left(k\right)={\left[\begin{array}{c}\hfill {\mathit{x}}^{⊺}\left(1,k\right),\dots ,{\mathit{x}}^{⊺}\left(N,k\right)\hfill \end{array}\right]}^{⊺},(36)$

which represents the true state value. Similarly, we can have updated unconstrained combined state estimate ${\widehat{\mathit{X}}}^{+}(k)={\left[\begin{array}{c}\hfill {\widehat{\mathit{x}}}^{+}{(1,k)}^{⊺},\dots ,{\widehat{\mathit{x}}}^{+}{(N,k)}^{⊺}\hfill \end{array}\right]}^{⊺}$ and constrained combined state estimate ${\tilde{\mathit{X}}}^{+}(k)={\left[\begin{array}{c}\hfill {\tilde{\mathit{x}}}^{+}{(1,k)}^{⊺},\dots ,{\tilde{\mathit{x}}}^{+}{(N,k)}^{⊺}\hfill \end{array}\right]}^{⊺}$. Then, we can have one constrained cooperative Kalman filter of x(j, k) as

$\begin{array}{cc}\hfill {\widehat{\mathit{X}}}^{-}\left(k\right)& =\mathit{A}\left(k-1\right){\tilde{\mathit{X}}}^{+}\left(k-1\right)+\mathit{U}\left(k-1\right),\hfill \\ \hfill {\mathit{R}}^{-}\left(k\right)& =\mathit{A}\left(k-1\right){\mathit{R}}^{+}\left(k-1\right){\mathit{A}}^{⊺}\left(k-1\right)+\mathit{W},\hfill \\ \hfill \mathit{K}\left(k\right)& ={\mathit{R}}^{-}\left(k\right){\mathit{C}}^{⊺}\left(k\right){\left[\mathit{C}\left(k\right){\mathit{R}}^{-}\left(k\right){\mathit{C}}^{⊺}\left(k\right)+\mathit{D}\left(k\right)Q{\mathit{D}}^{⊺}\left(k\right)+{\mathit{R}}_n\right]}^{-1},\hfill \\ \hfill {\widehat{\mathit{X}}}^{+}\left(k\right)& ={\widehat{\mathit{X}}}^{-}\left(k\right)+\mathit{K}\left(k\right)\left(\mathit{P}\left(k\right)-\mathit{C}\left(k\right){\widehat{\mathit{X}}}^{-}\left(k\right)-\mathit{D}\left(k\right)\widehat{\mathit{H}}\left(k\right)\right),\hfill \\ \hfill {\left({\mathit{R}}^{+}\left(k\right)\right)}^{-1}& ={\left({\mathit{R}}^{-}\left(k\right)\right)}^{-1}+{\mathit{C}}^{⊺}\left(k\right){\left[\mathit{D}\left(k\right)\mathit{Q}{\mathit{D}}^{⊺}\left(k\right)+{\mathit{R}}_n\right]}^{-1}\mathit{C}\left(k\right),\hfill \\ \hfill {\tilde{\mathit{X}}}^{+}\left(k\right)& ={\widehat{\mathit{X}}}^{+}\left(k\right)-\mathit{G}\left(k\right){\left({\mathit{G}}^{⊺}\left(k\right)\mathit{G}\left(k\right)\right)}^{-1}\left({\mathit{G}}^{⊺}\left(k\right){\widehat{\mathit{X}}}^{+}\left(k\right)-\widehat{\mathit{g}}\left(k\right)\right),\hfill \end{array}(37)$

where A(k) = diag(A(1, k), …, A(N, k)), C(k) = diag(C(1, k), …, C(N, k)), $\mathit{D}(k)={\left[\begin{array}{c}\hfill {\mathit{D}}^{⊺}(1,k),\dots ,{\mathit{D}}^{⊺}(N,k)\hfill \end{array}\right]}^{⊺}$, $\mathit{P}(k)={\left[\begin{array}{c}\hfill P^{⊺}(1,k),\dots ,P^{⊺}(N,k)\hfill \end{array}\right]}^{⊺}$, $\widehat{\mathit{H}}(k)={\left[\begin{array}{c}\hfill {\widehat{\mathit{H}}}^{⊺}(1,k),\dots ,{\widehat{\mathit{H}}}^{⊺}(N,k)\hfill \end{array}\right]}^{⊺}$, G(k) = diag(G(1, k), …, G(N, k)), and $\widehat{\mathit{g}}(k)={\left[\begin{array}{c}\hfill \widehat{\mathit{g}}(1,k),\dots ,\widehat{\mathit{g}}(N,k)\hfill \end{array}\right]}^{⊺}$.

Algorithm 1. Distributed constrained cooperative Kalman filter

For each individual cell C_j, j = 1, …, N:

Initialize the constrained state estimate, ${\tilde{\mathit{x}}}^{+}(j,0)$;

Initialize the locations of mobile sensors, ${\{{\mathit{r}}_{i,j}^1\}}_{i=1,\dots ,|C_j|}$;

Initialize time step k = 1;

1: while true do

2: Receive measurements ${\{z_{i,j}^k\}}_{i\in C_j}$ from mobile sensors in cell C_j and generate unconstrained state estimation ${\widehat{\mathit{x}}}^{+}(j,k)$ according to Eqs 28–32 (intra-cell communication);

3: Communicate the unconstrained state estimation ${\widehat{\mathit{x}}}^{+}(j,k)$ to the neighboring cells of C_j and obtain the constrained state estimate ${\tilde{\mathit{x}}}^{+}(j,k)$ using Eq. 34 (inter-cell communication);

4: Move to the next mobile sensor locations, ${\{{\mathit{r}}_{i,j}^{k+1}\}}_{i=1,\dots ,|C_j|}$;

5: k≔k + 1;

6: end while

One challenge encountered by the centralized cooperative Kalman filter is the high computation cost for large-scale systems. Our proposed distributed filtering strategy can avoid the problem of high computation costs when applied in large-scale systems. One possible solution is that we split the mobile sensors into two types—“computing, sensing, and communicating” (CSC) mobile sensors and “sensing and communicating” (SC) mobile sensors. CSC sensors take measurements and compute all the Kalman filter estimations, while SC sensors only collect measurements. Each SC sensor is required to share a cell with one CSC sensor. Each communication cell has at least one CSC sensor and multiple SC sensors.

The communication and data flow within the distributed mobile sensor network can be described as follows (also shown in Figure 3):

1. All sensors (SC and CSC) move and take a measurement of the field.

2. SC sensors (blue dots) in each cell send their data to their connected CSC sensors (green dots) (Figure 3A).

3. The CSC sensor performs the unconstrained cooperative Kalman filter to estimate the measurement values and gradients at the cell center (orange triangles) for each cell it is connected to.

4. The CSC sensors communicate with their “connected” CSC neighbors to get all the necessary sensor positions, measurements, and cell center estimates. (“connected” CSC neighbors are CSC mobile sensors that share a mobile sensor in their connected cells) (Figure 3B).

5. Each CSC sensor computes the constrained cooperative Kalman filter for each connected cell.

6. Cell estimates are pushed back to the connected mobile sensors for feedback control (Figure 3C).

7. Repeat.

FIGURE 3

FIGURE 3. Communication flow of the distributed mobile sensor network during estimation and velocity control. Green dots represent the CSC sensors, while blue dots represent the SC sensors, and orange triangles represent the cell formation centers. Green dashed arrows represent communication or information sharing with direction. Orange and dark red dashed lines represent unconstrained and constrained state estimations, respectively. (A) Intra-cell communication for unconstrained estimation. (B) Inter-cell communication for constrained estimation. (C) Intra-cell communication for velocity control.

6 Formation control design

In this section, we will design a formation controller for the mobile sensor network using the gradient estimation obtained from the distributed constrained cooperative Kalman filter.

Denote the desired distance between two sensors i, i′ as d_ii′. According to the state estimation using the constrained cooperative Kalman filter, we are able to obtain gradient estimation at each cell center. Thus, instead of directly using the gradient of each sensor, we consider utilizing the gradient estimation of cell centers related to the mobile sensor.

If sensor i belongs to cell j, we write it as i ∈ C_j. Define ${\mathcal{J}}_i=\{j|i\in C_j\}$ as the set that contains all the cell indices that sensor i is related to. Then, define $v_i^k$ as a weighted summation of gradient estimations from connected cell centers:

${\mathit{v}}_i^k\triangleq \sum _{j\in {\mathcal{J}}_i}\frac{1}{\Vert {\mathit{r}}_i^k-{\mathit{r}}_{c_j}^k{\Vert }_2}\nabla z\left({\mathit{r}}_{c_j}^k\right),(38)$

which implies that if r_i is closer to cell center ${\mathit{r}}_{c_j}$ than ${\mathit{r}}_{c_{j^{\prime }}}$, then estimation $\nabla z({\mathit{r}}_{c_j})$ has more influence on $v_i^k$ than $\nabla z({\mathit{r}}_{c_{j^{\prime }}})$ does. This describes a dynamical way of gradient approximation determined by distances.

By applying the gradient approximation in Eq. 38, the gradient-based desired distance can be designed as follows:

$d_{i{i}^{\prime }}=\frac{{ϵ}_1}{{‖\frac{1}{2}\left({\mathit{v}}_i^k+{\mathit{v}}_{i^{\prime }}^k\right)‖}_2+{ϵ}_2}+{ϵ}_3\in \left({ϵ}_3,\frac{{ϵ}_1}{{ϵ}_2}+{ϵ}_3\right),(39)$

where ${ϵ}_1,{ϵ}_2,{ϵ}_3\in {\mathbb{R}}_{+}$ are tuning parameters.

The desired distances between two communicating sensors Eq. 39 are designed to be gradient-dependent. Moreover, if the norm of the approximated gradient increases, the desired distance will decrease, which indicates that the formation size will shrink while encountering an area of larger field variation. In this way, a smaller formation will provide more accurate state estimation by collecting closer sensor measurements and can avoid ignoring large variations within a small domain. When the sensors are exploring a field with decreasing gradient norm, larger distances can still ensure the accuracy of state estimation due to the slow field variation. Both cases demonstrate the ability of the proposed gradient-based adaptive formation to be suitable for various scenarios, such as source-seeking and contour mapping (Han and Chen, 2014), where the rigid formation fails to provide accurate estimation for a field with a large variation or the field is too simple for the capability of the deployed sensors. In this paper, we are focusing on the application of source-seeking using distributed mobile sensor networks. Additionally, the tuning parameters in Eq. 39 ensure the boundedness of the desired distances so that the formation will be controllable, not shrinking to zero size or becoming too large.

Since the mobile sensors are taking discrete measurements, we will apply discrete time control written as follows:

${\mathit{r}}_{i,j}^{k+1}={\mathit{r}}_{i,j}^k+{\mathit{u}}_i^k\mathrm{\Delta }t,(40)$

where Δt is the sampling time and is set to be constant. Define $\mathit{w}=\frac{\nabla z+e}{\Vert \nabla z+e{\Vert }_2}$ as the unit length direction vector along the gradient direction ∇z perturbed by some estimation error e.

We approximate the gradient direction at r_i using the gradient estimations at all the cell centers around it as

${\mathit{w}}_i^k=\frac{{\mathit{v}}_i^k}{\Vert {\mathit{v}}_i^k{\Vert }_2}=\frac{{\sum }_{j\in {\mathcal{J}}_i^v}\frac{1}{\Vert {\mathit{r}}_i-{\mathit{r}}_{c_j}{\Vert }_2}\nabla z\left({\mathit{r}}_{c_j}^k\right)}{{‖{\sum }_{j\in {\mathcal{J}}_i^v}\frac{1}{\Vert {\mathit{r}}_i-{\mathit{r}}_{c_j}{\Vert }_2}\nabla z\left({\mathit{r}}_{c_j}^k\right)‖}_2}.(41)$

The control input for sensor i is designed as

${\mathit{u}}_i^k=-a{\mathit{w}}_i^k+b{\left({\mathit{w}}_i^k\right)}^{\perp }+c{\mathit{f}}_i^k,(42)$

where a, b, c are tuning parameters and ${({\mathit{w}}_i^k)}^{\perp }$ is the normal unit vector of ${\mathit{w}}_i^k$. The formation terms ${\mathit{f}}_i^k$ are defined as

${\mathit{f}}_i^k=\sum _{i^{\prime }\ne i}\frac{\Vert {\mathit{r}}_{i^{\prime }}^k-{\mathit{r}}_i^k{\Vert }_2-d_{i{i}^{\prime }}^k}{\Vert {\mathit{r}}_{i^{\prime }}^k-{\mathit{r}}_i^k{\Vert }_2}\left({\mathit{r}}_{i^{\prime }}^k-{\mathit{r}}_i^k\right),(43)$

where $d_{i{i}^{\prime }}^k$ is a desired separation distance between sensors i and i′. The formation term enables the mobile sensors to maintain formations close to the desired ones, by letting sensors move away from each other if the actual distance is larger than the desired one and letting sensors move closer if the actual distance is smaller.

Remark 6.1: The control input ${\mathit{u}}_i^k$ defined in Eq. 42 is designed for source-seeking, and it enables the sensors to follow the gradient direction and reach the source, which is the main application considered in this paper. In addition to source-seeking, mobile sensors can also be deployed for other applications, such as contour mapping, where sensors are required to follow the normal direction of the gradient. The control input in Eq. 42 can be modified based on the speeding-up and slowing-down (SUSD) strategy (Wu and Zhang, 2015; Al-Abri et al., 2018; Zhang et al., 2020).

${\mathit{u}}_i^k=-a\left(p_i^k-z_d\right){\mathit{w}}_i^k+b{\left({\mathit{w}}_i^k\right)}^{\perp }+c{\mathit{f}}_i^k,(44)$

and the only change is the multiplication of $(p_i^k-z_d)$, where $p_i^k$ is the noisy measurement of the field taken by agent i at the kth time step and z_d is the field value of the desired level.

7 Convergence analysis

In this section, convergence analysis for the formation control and the distributed constrained cooperative Kalman filter will be provided.

Theorem 7.1: Under the formation control (Eq. 42) along with (Eq. 43) and (Eq. 39), the mobile sensors can always maintain bounded formations with bounded speed.

Proof: In Oh et al. (2015), the controller defined by Eq. 43 is a well-understood control law that is guaranteed to achieve the formation specified by the desired constant inter-agent distances. Since $d_{i{i}^{\prime }}^k\in ({ϵ}_3,\frac{{ϵ}_1}{{ϵ}_2}+{ϵ}_3)$ is bounded by Eq. 39, the formation can be achieved with bounded error. Then, there exists ${ϵ}_4\in {\mathbb{R}}_{+}$ such that for any two neighboring sensors,

$\Vert {\mathit{r}}_{i^{\prime }}^k-{\mathit{r}}_i^k{\Vert }_2-d_{i{i}^{\prime }}^k\in \left[-{ϵ}_4,{ϵ}_4\right].(45)$

This implies that

$\Vert {\mathit{r}}_{i^{\prime }}^k-{\mathit{r}}_i^k{\Vert }_2\in \left[{ϵ}_3-{ϵ}_4,\frac{{ϵ}_1}{{ϵ}_2}+{ϵ}_3+{ϵ}_4\right],(46)$

where the parameter ϵ₃ can be chosen such that ϵ₃ > ϵ₄.Since the formation is bounded by Eq. 45, the formation term satisfies

$\begin{array}{ll}\hfill \Vert {\mathit{f}}_i^k{\Vert }_2& ={‖\sum _{i^{\prime }\in {\mathcal{N}}_i}\frac{\Vert {\mathit{r}}_{i^{\prime }}^k-{\mathit{r}}_i^k{\Vert }_2-d_{i{i}^{\prime }}^k}{\Vert {\mathit{r}}_{i^{\prime }}^k-{\mathit{r}}_i^k{\Vert }_2}\left({\mathit{r}}_{i^{\prime }}^k-{\mathit{r}}_i^k\right)‖}_2\hfill \\ \hfill & \le \sum _{i^{\prime }\in {\mathcal{N}}_i}{‖\frac{\Vert {\mathit{r}}_{i^{\prime }}^k-{\mathit{r}}_i^k{\Vert }_2-d_{i{i}^{\prime }}^k}{\Vert {\mathit{r}}_{i^{\prime }}^k-{\mathit{r}}_i^k{\Vert }_2}\left({\mathit{r}}_{i^{\prime }}^k-{\mathit{r}}_i^k\right)‖}_2\le \sum _{i^{\prime }\in {\mathcal{N}}_i}{ϵ}_4.\hfill \end{array}(47)$

According to the discrete time control law in Eq. 42, the norm of control input ${\mathit{u}}_i^k$ satisfies

$\begin{array}{ll}\hfill \Vert {\mathit{u}}_i^k{\Vert }_2& =\Vert -a{\mathit{w}}_i^k+b{\left({\mathit{w}}_i^k\right)}^{\perp }+c{\mathit{f}}_i^k{\Vert }_2\le \Vert a{\mathit{w}}_i^k{\Vert }_2+\Vert b{\left({\mathit{w}}_i^k\right)}^{\perp }{\Vert }_2\hfill \\ \hfill & +\Vert c{\mathit{f}}_i^k{\Vert }_2\le |a|+|b|+|c|\sum _{i^{\prime }\in {\mathcal{N}}_i}{ϵ}_4,\hfill \end{array}(48)$

where the last inequality can be obtained from Eq. 47 and $\Vert {\mathit{w}}_i^k{\Vert }_2=\Vert {({\mathit{w}}_i^k)}^{\perp }{\Vert }_2=1$.Therefore, the mobile sensors are capable of maintaining bounded formation with bounded speed by Eqs 46, 48.The information dynamics and measurement equation considered for each individual cell in this paper share the same structure as those in Zhang et al. (2022), where a centralized cooperative Kalman filter has been proposed with provable convergence.

Proposition 7.2: [Propositions VI.4 and VI.6 in Zhang et al. (2022)] The state dynamics (13) with the measurement equation (17) are uniformly completely controllable and uniformly completely observable if the following conditions are satisfied: (Cd1) The covariance matrix W is bounded, i.e., λ₁I ≤W ≤ λ₂I for some constants λ₁, λ₂ > 0. (Cd2) The speed of each agent is uniformly bounded, i.e., $\Vert {\mathit{r}}_{i,j}^k-{\mathit{r}}_{i,j}^{k-1}{\Vert }_2\le {\lambda }_3$ for all i, j, k, and some constant λ₃ > 0. (Cd3) The number of sensors in one cell C_j satisfies $|C_j|>\frac{{\mathrm{d}}^2}{2}$ for all j. (Cd4) The covariance matrices R_n and Q are bounded, i.e., λ₄I ≤ R_n ≤ λ₅I and 0 ≤ Q ≤ λ₆I for some λ₄, λ₅, λ₆> 0. (Cd5) The distance between each agent and the cell center is uniformly bounded from both above and below, i.e.,${\lambda }_7\le \Vert {\mathit{r}}_{i,j}^{k-1}-{\mathit{r}}_{c_j}^{k-1}\Vert \le {\lambda }_8$ for all i, j, k, and some λ₇, λ₈ > 0. (Cd6) There exists a constant time difference τ₂, and for all k > τ₂, there exists a time instance k₁ ∈ [k − τ₂, k], as well as two sets of agents indexed by {i₁, …, i_d}, {i_d+1, …, i_2d}, respectively, such that $({\mathit{r}}_{i_1,j}^{k_1-1}-{\mathit{r}}_{c_j}^{k_1-1}),\dots ,({\mathit{r}}_{i_{\text{d}},j}^{k_1-1}-{\mathit{r}}_{c_j}^{k_1-1})$ are linearly independent, and $({\mathit{r}}_{i_{d+1,j}}^{k_1}-{\mathit{r}}_{c_j}^{k_1-1}),\dots ,({\mathit{r}}_{i_{2\mathrm{d},j}}^{k_1}-{\mathit{r}}_{c_j}^{k_1-1})$ are linearly independent.

Remark 7.3: In Proposition 7.2, (Cd2) and (Cd5) are guaranteed by bounded sensor speed and bounded formation cells under (Eq. 40) along with (Eq. 43).Since the unconstrained cooperative Kalman filter is both uniformly completely controllable and observable (Proposition 7.2), the unconstrained filter for each individual cell is convergent. This means that $\Vert {\mathit{x}}_j-{\widehat{\mathit{x}}}_j\Vert$ is bounded for all j, where x_j represents the true state value.Since the combined filter is constructed by stacking all distributed filters together and the unconstrained filter for each individual cell is convergent, the combined unconstrained filter is also convergent.

$\Vert \mathit{X}-\widehat{\mathit{X}}{\Vert }_2=\sqrt{\sum _{j=1}^N\Vert \mathit{x}-{\mathit{x}}_j{\Vert }_2^2}\le \sum _{j=1}^N\Vert \mathit{x}-{\mathit{x}}_j{\Vert }_2.(49)$

From Theorem 4 in Simon and Chia (2002), since the unconstrained combined cooperative Kalman filter is convergent, the constrained combined cooperative Kalman filter is also convergent. Thus, each distributed constrained cooperative Kalman filter is convergent.

8 Simulation results

In this section, simulation results using a distributed mobile sensor network of 12 sensors will be presented to validate that the proposed distributed algorithm enables mobile sensor networks to estimate information at each cell center along trajectories of a collection of cell centers. According to the 1D advection–diffusion solution in Mojtabi and Deville (2015), we consider a 2D advection–diffusion equation as $z(\mathit{r},t)=-\sin (\pi (k_1{\mathit{r}}_x+v_1{k}_{2}t))e^{-\theta {\pi }^2{k}_{2}t}-\sin (\pi (k_1{\mathit{r}}_y+k_2{v}_{2}t))e^{-\theta {\pi }^2{k}_{2}t}$ and $\mathit{r}={[{\mathit{r}}_x,{\mathit{r}}_y]}^{⊺}$, with diffusion coefficient θ = 2 and flow velocity v = [0.5,1]^⊺, k₁ = 0.05, k₂ = 0.0001, and the SNR for the field is set to be 10. The sample rate is 1 Hz, and the simulation of 12 sensors will be running for 600 s.

At each time step, mobile sensors will take noisy measurements of the field, and the sensor locations are also available. Based on the network structure of the mobile sensors (Figure 4), three or four sensor measurements will be incorporated for estimation within each cell. The estimation for field value and gradient at each cell center is composed of two parts: unconstrained and constrained estimations. In order to complete unconstrained estimation (Eqs. 28–32), noisy measurements, previous and current sensor locations, and previous estimations pertaining to each cell are required. Then, the constrained estimation can be obtained by projecting the unconstrained estimation onto the constrained space (Eq. 26) defined by the PDE field, where information sharing of unconstrained estimation between neighboring cells is required to have the PDE constraint.

FIGURE 4

FIGURE 4. Trajectory of cell centers in a mobile sensor network. (A) Initial 9-cell formation shape of a 12-sensor network. (B) Trajectories of cell centers with initial formation (bottom left smaller one) and final formation (upper left larger one), where the colored dashed lines represent the level curve of the field. For simulation videos, please see https://youtu.be/aHrvXjt8u4M (same initial location), and https://youtu.be/8_gvyacW-5g (different initial locations).

As shown in Figure 4A, the mobile sensors have been organized into four triangle cells and five rectangle cells, where the blue dots represent the mobile sensors, the orange triangles represent cell centers, and the blue lines represent the edges of each individual cell. The shape of each cell is fixed due to the predetermined communication graph by Assumption 2.1, and there is no requirement for the cell shapes to be symmetric. The five rectangles in the initial formation are symmetric, while the four triangles are asymmetric.

In Figure 4B, the dashed contours represent the level curves of the field values, and the colored lines are the trajectories of the nine cell centers. The smaller formation on the bottom left is the initial location, and the larger formation on the upper right is the final location. Driven by Eq. 42, the mobile sensors will adjust the distances between neighbors according to the norms of estimated gradients, which leads to non-rigid gradient-based formation for the mobile sensor networks and also explains the size difference between initial and final formations.

The field estimations generated by the proposed distributed constrained filtering at the cell centers of cell 2 and cell 8 along trajectories are provided in Figure 5. The noisy field values indicated by the blue dashed lines are relatively inaccurate, compared with the true field values marked by purple solid lines. Both the unconstrained and constrained filtering strategies (marked by red and yellow, respectively) are capable of reducing the noise. These four lines in each plot of Figure 5 share the same trajectory. At each time step, we run the distributed constrained cooperative Kalman filter in Eq. 37 for each cell, and then the formation control in Eq. 42 is updated based on the estimated gradient information. At the same time, a distributed unconstrained filter is also running for each cell to obtain an unconstrained state estimation, but such estimation will not be used to update the formation control. The constrained one shows improved performance, which demonstrates that more accurate estimation can be generated when more field information is incorporated by adding PDE constraints. We also observe the differences in estimation performance among cells, as shown by the two example cells in Figure 5. According to Figure 4A, cell 2 is the center cell of the formation, and it has eight neighboring cells, which provide information for updated constrained state estimation, while cell 8 only has three neighboring cells. This explains why when comparing with unconstrained estimations, the constrained estimations in cell 8 show less improvement than the constrained estimations in cell 2.

FIGURE 5

FIGURE 5. Field estimation generated by distributed constrained cooperative Kalman filter (yellow lines) at the cell center of cell 2 and cell 8 along trajectories, compared with field value with noise (blue dashed lines), distributed unconstrained cooperative Kalman filter estimation (red lines) and true field value (purple lines). The x-axis is time t, and the y-axis is the field value z(r, t). (A) Cell 2. (B) Cell 8.

Statistical data of estimation errors at cell centers ${\mathit{r}}_{c_1},\dots ,{\mathit{r}}_{c_9}$ are provided in Table 1 along with the corresponding data of measurement noise. The subscripts n, u, c represent measurement noise, unconstrained estimation error, and constrained estimation error, respectively. The first row denotes the nine cells. Rows 2–4 represent the means of the errors between true field value and noise, unconstrained estimation, and constrained estimation, respectively. The last three rows represent the standard deviations of the errors between true field value and noise, unconstrained estimation, and constrained estimation, respectively.

TABLE 1

TABLE 1. Statistical data of state estimation errors and noise.

According to Table 1, both unconstrained and constrained filtering improved the estimation of z(r, t), and the constrained filtering had a lower standard deviation than unconstrained filtering in all cells. This validates that more accurate estimation can be generated when more field information is provided.

9 Conclusion and future work

In this paper, we proposed a distributed cooperative Kalman filter constrained by the advection–diffusion equation, which solves the estimation problem using a large number of mobile sensors in a distributed way. Both convergence analysis and simulation results demonstrate the effectiveness of our proposed filtering strategy, and the mobile sensors can maintain gradient-based formation. For the next step, we will generalize our method to higher-order PDEs and incorporate parameter identification for the PDE models.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author contributions

ZZ, SM, WW, and FZ conceived the study and formulated the problem. ZZ, WW, and FZ designed the filtering algorithm. SM performed the simulation experiments. ZZ and SM contributed to the formation control design and simulation evaluation. ZZ conducted the FVM approximation and the convergence analysis and wrote the draft of the manuscript. All authors contributed to the manuscript revision, read, and approved the submitted version.

Funding

The research work is supported by ONR grants N00014-19-1-2556 and N00014-19-1-2266; AFOSR grant FA9550-19-1-0283; NSF grants GCR-1934836, CNS-2016582, ITE-2137798, CMMI-1917300, and RINGS-2148353; and NOAA grant NA16NOS0120028.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Al-Abri, S., Wu, W., and Zhang, F. (2018). A gradient-free three-dimensional source seeking strategy with robustness analysis. IEEE Trans. Automatic Control64, 3439–3446. doi:10.1109/tac.2018.2882172

CrossRef Full Text | Google Scholar

Alfouzan, F. A. (2021). Energy-efficient collision avoidance mac protocols for underwater sensor networks: Survey and challenges. J. Mar. Sci. Eng.9, 741. doi:10.3390/jmse9070741

CrossRef Full Text | Google Scholar

Boufadel, M. C., Özgökmen, T., Socolofsky, S. A., Kourafalou, V. H., Liu, R., and Lee, K. (2022). Oil transport following the deepwater Horizon blowout. doi:10.1146/annurev-marine-040821

CrossRef Full Text | Google Scholar

Chen, J., and Dames, P. (2020). “Distributed and collision-free coverage control of a team of mobile sensors using the convex uncertain voronoi diagram,” in 2020 American control conference (ACC) (IEEE), 5307–5313.

CrossRef Full Text | Google Scholar

Demetriou, M. A. (2021). “Employing mobile sensor density to approximate state feedback kernels in static output feedback control of pdes,” in 2021 American control conference (ACC) (IEEE), 2775–2781.

CrossRef Full Text | Google Scholar

Droniou, J. (2014). Finite volume schemes for diffusion equations: Introduction to and review of modern methods. Math. Models Methods Appl. Sci.24, 1575–1619. doi:10.1142/s0218202514400041

CrossRef Full Text | Google Scholar

Elston, J., Argrow, B., Frew, E., and Houston, A. (2010). “Evaluation of UAS concepts of operation for severe storm penetration using hardware-in-the-loop simulations,” in AIAA guidance, navigation, and control conference (Reston, Virigina: American Institute of Aeronautics and Astronautics). doi:10.2514/6.2010-8178

CrossRef Full Text | Google Scholar

Eshghi, M., and Schmidtke, H. R. (2018). An approach for safer navigation under severe hurricane damage. J. Reliab. Intelligent Environ.4, 161–185. doi:10.1007/s40860-018-0066-1

CrossRef Full Text | Google Scholar

Gao, X., and Acar, L. (2016). Using a mobile robot with interpolation and extrapolation method for chemical source localization in dynamic advection-diffusion environment. Int. J. Robotics Automation (IJRA)5, 87–97. doi:10.11591/ijra.v5i2.pp87-97

CrossRef Full Text | Google Scholar

Ge, X., Han, Q.-L., Zhang, X.-M., Ding, L., and Yang, F. (2019). Distributed event-triggered estimation over sensor networks: A survey. IEEE Trans. Cybern.50, 1306–1320. doi:10.1109/tcyb.2019.2917179

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghez, R. (2001). Diffusion phenomena: Cases and studies. Springer Science & Business Media.

Google Scholar

Han, J., and Chen, Y. (2014). Multiple uav formations for cooperative source seeking and contour mapping of a radiative signal field. J. Intelligent Robotic Syst.74, 323–332. doi:10.1007/s10846-013-9897-4

CrossRef Full Text | Google Scholar

He, X., Hu, C., Hong, Y., Shi, L., and Fang, H.-T. (2020). Distributed kalman filters with state equality constraints: Time-based and event-triggered communications. IEEE Trans. Automatic Control65, 28–43. doi:10.1109/tac.2019.2906462

CrossRef Full Text | Google Scholar

Hermeline, F. (2000). A finite volume method for the approximation of diffusion operators on distorted meshes. J. Comput. Phys.160, 481–499. doi:10.1006/jcph.2000.6466

CrossRef Full Text | Google Scholar

Hu, W., Demetriou, M. A., Tian, X., and Gatsonis, N. A. (2022). Hybrid domain decomposition filters for advection–diffusion pdes with mobile sensors. Automatica138, 110109. doi:10.1016/j.automatica.2021.110109

CrossRef Full Text | Google Scholar

Ji, H., Lewis, F. L., Hou, Z., and Mikulski, D. (2017). Distributed information-weighted kalman consensus filter for sensor networks. Automatica77, 18–30. doi:10.1016/j.automatica.2016.11.014

CrossRef Full Text | Google Scholar

Krstic, M., and Smyshlyaev, A. (2008). Adaptive control of pdes. Annu. Rev. Control32, 149–160. doi:10.1016/j.arcontrol.2008.05.001

CrossRef Full Text | Google Scholar

Kumar, D. P., Amgoth, T., and Annavarapu, C. S. R. (2019). Machine learning algorithms for wireless sensor networks: A survey. Inf. Fusion49, 1–25. doi:10.1016/j.inffus.2018.09.013

CrossRef Full Text | Google Scholar

Lian, B., Wan, Y., Zhang, Y., Liu, M., Lewis, F. L., and Chai, T. (2020). Distributed kalman consensus filter for estimation with moving targets. IEEE Trans. Cybern.52, 5242–5254. doi:10.1109/tcyb.2020.3029007

PubMed Abstract | CrossRef Full Text | Google Scholar

Mishra, P. K., Chatterjee, D., and Quevedo, D. E. (2020). Stochastic predictive control under intermittent observations and unreliable actions. Automatica118, 109012. doi:10.1016/j.automatica.2020.109012

CrossRef Full Text | Google Scholar

Mojtabi, A., and Deville, M. O. (2015). One-dimensional linear advection–diffusion equation: Analytical and finite element solutions. Comput. Fluids107, 189–195. doi:10.1016/j.compfluid.2014.11.006

CrossRef Full Text | Google Scholar

Oh, K.-K., Park, M.-C., and Ahn, H.-S. (2015). A survey of multi-agent formation control. Automatica53, 424–440. doi:10.1016/j.automatica.2014.10.022

CrossRef Full Text | Google Scholar

Olfati-Saber, R. (2009). “Kalman-consensus filter: Optimality, stability, and performance,” in Proceedings of the 48h IEEE conference on decision and control (CDC) held jointly with 2009 28th Chinese control conference (IEEE), 7036–7042.

CrossRef Full Text | Google Scholar

Omeke, K. G., Mollel, M. S., Ozturk, M., Ansari, S., Zhang, L., Abbasi, Q. H., et al. (2021). Dekcs: A dynamic clustering protocol to prolong underwater sensor networks. IEEE Sensors J.21, 9457–9464. doi:10.1109/jsen.2021.3054943

CrossRef Full Text | Google Scholar

Park, H., Liu, J., Johnson-Roberson, M., and Vasudevan, R. (2018). “Robust environmental mapping by mobile sensor networks,” in 2018 IEEE international conference on Robotics and automation (ICRA) (IEEE), 2395–2402.

CrossRef Full Text | Google Scholar

Radmanesh, M., Sharma, B., Kumar, M., and French, D. (2021). PDE solution to UAV/UGV trajectory planning problem by spatio-temporal estimation during wildfires. Chin. J. Aeronautics34, 601–616. doi:10.1016/j.cja.2020.11.002

CrossRef Full Text | Google Scholar

Simon, D., and Chia, T. L. (2002). Kalman filtering with state equality constraints. IEEE Trans. Aerosp. Electron. Syst.38, 128–136. doi:10.1109/7.993234

CrossRef Full Text | Google Scholar

Song, W., Wang, J., Zhao, S., and Shan, J. (2019). Event-triggered cooperative unscented kalman filtering and its application in multi-uav systems. Automatica105, 264–273. doi:10.1016/j.automatica.2019.03.029

CrossRef Full Text | Google Scholar

Wang, J., Gao, Y., Liu, W., Sangaiah, A. K., and Kim, H.-J. (2019). An intelligent data gathering schema with data fusion supported for mobile sink in wireless sensor networks. Int. J. Distributed Sens. Netw.15, 155014771983958. doi:10.1177/1550147719839581

CrossRef Full Text | Google Scholar

Wei, C., Tanner, H. G., Yu, X., and Hsieh, M. A. (2019). “Low-range interaction periodic rendezvous along Lagrangian coherent structures,” in 2019 American control conference (ACC) (IEEE), 4012–4017.

CrossRef Full Text | Google Scholar

Wu, W., You, J., Zhang, Y., Li, M., and Su, K. (2021). Parameter identification of spatial–temporal varying processes by a multi-robot system in realistic diffusion fields. Robotica39, 842–861. doi:10.1017/s0263574720000788

CrossRef Full Text | Google Scholar

Wu, W., and Zhang, F. (2015). A speeding-up and slowing-down strategy for distributed source seeking with robustness analysis. IEEE Trans. Control Netw. Syst.3, 231–240. doi:10.1109/tcns.2015.2459414

CrossRef Full Text | Google Scholar

You, J., Zhang, Z., Zhang, F., and Wu, W. (2022). Cooperative filtering and parameter identification for advection–diffusion processes using a mobile sensor network. IEEE Trans. Control Syst. Technol.31, 527–542. doi:10.1109/tcst.2022.3183585

CrossRef Full Text | Google Scholar

Zhang, F., and Leonard, N. E. (2010). Cooperative filters and control for cooperative exploration. IEEE Trans. Automatic Control55, 650–663. doi:10.1109/tac.2009.2039240

CrossRef Full Text | Google Scholar

Zhang, Z., Al-Abri, S., Wu, W., and Zhang, F. (2020). “Level curve tracking without localization enabled by recurrent neural networks,” in 2020 5th international conference on automation, control and Robotics engineering (CACRE) (IEEE), 759–763.

CrossRef Full Text | Google Scholar

Zhang, Z., Mayberry, S. T., Wu, W., and Zhang, F. (2023). “Distributed cooperative kalman filter constrained by discretized Poisson equation for mobile sensor networks,” in 2023 American control conference (ACC) (IEEE). accepted.

Google Scholar

Zhang, Z., Wu, W., and Zhang, F. (2022). “Cooperative filtering and parameter estimation for polynomial pdes using a mobile sensor network,” in 2022 American control conference (ACC) (IEEE), 982–987.

CrossRef Full Text | Google Scholar