© The Author(s) 2025. Published by Science Press and EDP Sciences.

INTRODUCTION

Recent developments in machine learning algorithms and sensor technologies have facilitated the integration of autonomous systems in real-world marine applications, such as environmental monitoring [1-3], disaster response [4, 5], and industrial automation [6, 7]. Unmanned surface vehicles (USVs) have emerged as indispensable tools, enabling a wide variety of missions [8-10]. Despite the tremendous progress in multi-USV fleet coordination, real-time monitoring remains a dilemma when USVs are noncooperative and do not share data with external monitors. Furthermore, practical deployments impose constraints on real-time efficiency and limited onboard calculation capability, necessitating low-complexity algorithms.

A global positioning system (GPS)-based tracking control system for wheeled mobile robots is introduced in such unmanned systems, which compensates for skidding and slipping effects using real time kinematic (RTK)-GPS, ensuring navigational trajectory tracking [11]. In Ref. [12], a vision-based target detection and localization system is developed, utilizing a cooperative swarm composed of UAVs and unmanned ground vehicles (UGVs), where UAVs employ an optical flow for motion detection and UGVs make individual detection. In Ref. [13], a vehicle monitoring system is developed by integrating an Arduino microcontroller with global system for mobile communications (GSM) and GPS modules. A framework for adaptive learning navigation with nested guidance layers is introduced in Ref. [14] for UAVs, enabling horizontal monitoring and vertical descent in confined landing zones using solely relative position feedback. However, these schemes depend on motion information directly provided by the non-cooperative monitored targets. To address this issue, a target-tracking control system for underactuated autonomous surface vehicles (ASVs) is proposed in Ref. [15], which relies solely on line-of-sight range and angle measurements. Moreover, this system integrates an extended state observer and a single hidden layer neural network to estimate both target dynamics and external disturbances. A monocular camera-based method was proposed in Ref. [16], leveraging optical flow for target localization, and integrating it with an extended Kalman filter (EKF) to estimate motion dynamics. Furthermore, Ref. [17] establishes a hierarchical coarse-to-fine deep reinforcement learning framework for UAV tracking, where a coarse stage initializes the bounding box, and a fine stage refines it to handle aspect-ratio, scale, and occlusion changes. However, these approaches rely on ideal models, which exhibit high computational complexity. To analyze multi-source data, Ref. [18] develops a UAV-based tracking-and-recognition system integrating consensus-based tracking, neural-network detection, and gimbal stabilization, where real-time tracking is achieved via multimodal fusion together with moving-background compensation. A tracking system for USVs is tailored, utilizing an EKF and a visibility-aware control strategy to enhance target detection, positioning accuracy, and trajectory prediction [19]. Additionally, Ref. [20] introduces a target detection method with the assistance of a single shot multibox detector, a support vector machine classifier, and a tracking algorithm. Despite these advancements, few existing studies address the scenario of multi-target coordinated monitoring.

To this end, we design a cooperative method for multi-USV systems. The UAVs are assigned according to USV-UAV pairwise matching to observe relative positions online. An efficient block sparse Bayesian learning algorithm (EBSBL) with low computational complexity is proposed to identify the coordinated dynamics of the multi-USV fleet, leveraging the advantages of sparse Bayesian learning (SBL) over traditional ${\ell }_1$ methods for sparse, high-quality signal recovery, and incorporating structural information for improved performance [21, 22]. Additionally, the unscented Kalman filter (UKF) is employed to facilitate real-time prediction, USV trajectory estimation, and UAV monitoring coordination. In summary, the contributions of this work are two-fold.

(1) Propose a real-time cross-domain monitoring method not requiring motion information provided by the multi-USV fleet.

(2) Propose an EBSBL with theoretically guaranteed feasibility.

The remainder of this paper is organized as follows. Section PRELIMINARIES AND PROBLEM FORMULATION introduces the problem addressed by the paper with necessary preliminaries. Section METHOD develops the monitoring scheme, which includes UAVs assignment, coordinated dynamics learning, and cross-domain coordinated tracking modules. Experiments are conducted on a self-established cross-domain platform in Section NUMERICAL AND EXPERIMENTAL RESULTS to demonstrate both the effectiveness and superiority of the proposed monitoring method. Finally, the conclusion is drawn in Section CONCLUSIONS.

PRELIMINARIES AND PROBLEM FORMULATION

Consider a multi-UAV-multi-USV scenario where $n$ USVs are monitored by $n$ UAVs, as shown in Figure 1. Denote the positions of the $i$-th USV ${\mathit{p}}_i(t)={\left[p_{i[x]}(t)\mathrm{,}{p}_{i[y]}(t)\right]}^{\text{T}}$, $i\in \{\mathrm{1,2,}\dots \mathrm{,}n\}$, where $p_{i[x]}(t)$ and $p_{i[y]}(t)$ represent the position along the $x$- and $y$-axes of the $i$-th USV, respectively. Denote the positions of the $j$-th UAV ${\mathit{q}}_j(t)={\left[q_{j[x]}(t)\mathrm{,}{q}_{j[y]}(t)\mathrm{,}{q}_{[z]}\right]}^{\text{T}}$, $j\in \{\mathrm{1,2,}\dots \mathrm{,}n\}$, where $q_{i[x]}(t)$ and $q_{i[y]}(t)$ represent the positions along the $x$- and $y$-axes of the $j$-th UAV, respectively, and $q_z$ represents the fixed altitude. Assume that the dynamics of USVs ${\mathit{M}}_i(t)$ is governed by the velocity function ${\mathit{f}}_i({\mathit{p}}_1(t)\mathrm{,}{\mathit{p}}_2(t)\dots \mathrm{,}{\mathit{p}}_n(t))$, which is widely applied in cooperative control of USVs as follows[23-25]:$\begin{array}{l}{\mathit{M}}_i(t)={\mathit{f}}_i({\mathit{p}}_1(t)\mathrm{,}{\mathit{p}}_2(t)\mathrm{,}\dots \mathrm{,}{\mathit{p}}_n(t))\mathrm{,}i=\mathrm{1,2,}\dots \mathrm{,}n\mathrm{.}\hfill \end{array}$(1)The dynamics of USVs are represented by the following kinematic model [26]:$\begin{array}{l}{\dot{p}}_i(t)={\mathit{g}}_i({\mathit{M}}_i(t))=\mathit{G}({\rho }_i){\left[{\iota }_i\mathrm{,}{\varrho }_i\right]}^{\text{T}}\mathrm{,}\mathit{G}({\rho }_i)=\left[\begin{array}{cc}\cos {\rho }_i& -\sin {\rho }_i\\ \sin {\rho }_i& \cos {\rho }_i\end{array}\right]\hfill \end{array}\mathrm{,}$(2)where $\mathit{G}({\rho }_i)$ represents the rotation matrix, and ${\rho }_i$, ${\iota }_i$, and ${\varrho }_i$ are the orientation angle, forward velocity, and transverse velocity, respectively. The dynamics of UAVs is modeled as follows:$\begin{array}{l}{\dot{q}}_j(t)={\mathit{\kappa }}_j(t)\mathrm{,}{\dot{\kappa }}_j(t)={\mathit{u}}_j(t)\mathrm{,}j=\mathrm{1,2,}\dots \mathrm{,}n\mathrm{,}\hfill \end{array}$(3)where ${\mathit{\kappa }}_j(t)$ and ${\mathit{u}}_j(t)$ denote the velocity and control input of the $j$-th UAV, respectively.

Figure 1 Diagram of the algorithm for UAVs monitoring (or tracking) USVs with EBSBL, consisting of three stages. Stage 1: Assign UAVs to USVs using auction algorithm and make observations. Stage 2: Identify USV dynamics using efficient block sparse Bayesian learning. Stage 3: Monitor USVs in coordination using identified results and UKF by UAVs.

Note that the USVs do not share their position and velocity information with the UAVs, which can only be observed by the UAVs. More precisely, during the observation period, each UAV can observe any USV, rather than being restricted to a fixed pairwise monitoring scheme. Define the relative position of the $i$-th USV observed by the $j$-th UAV at time $t$ as ${\mathit{r}}_{ij}(t)=\left[r_{ij[x]}(t)\mathrm{,}{r}_{ij[y]}(t)\right]$, where $r_{ij[x]}(t)$ and $r_{ij[y]}(t)$ represent the relative positions along the $x$- and $y$-axes, respectively. Denote $\tilde{t}$ as the time when the observation is not available. The problem addressed by this paper is motivated as below.

Problem 1: Monitor the USVs by identifying the dynamics of USVs ${\mathit{M}}_i(t)$ and predicting their positions ${\mathit{p}}_i(\tilde{t})$ based on the relative observed data ${\mathit{r}}_{ij}(t)$ and the positions ${\mathit{q}}_j(t)$ of UAVs, i.e., ${\dot{p}}_i(\tilde{t})={\mathit{g}}_i({\mathit{M}}_i(\tilde{t}))$.

METHOD

UAVs assignment for tracking USVs

To enable trajectory observation and tracking of USVs, each USV is assigned to a unique UAV at each observation time $t_h$, which inspires to the following pairwise matching optimization problem:$\underset{\alpha }{{\displaystyle \min }}{\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{j=1}^n}\parallel {\mathit{r}}_{ij}(t_h)\parallel {\alpha }_{ij}\mathrm{,}$(4a)$\text{s.t}\mathrm{.}{\displaystyle \sum _{j=1}^n}{\alpha }_{ij}=\mathrm{1,}{\displaystyle \sum _{i=1}^n}{\alpha }_{ij}=\mathrm{1,}{\alpha }_{ij}\in \left\{\mathrm{0,1}\right\}\mathrm{,}$(4b)where ${\alpha }_{ij}$ is a binary variable equal to $1$ if the $i$-th USV is assigned to the $j$-th UAV, and otherwise. The objective function (4a) seeks to minimize the total observation distance, given that the quality of UAV-collected data deteriorates with increasing distance. Furthermore, when the UAV is closer to the target, it can more rapidly follow the trajectory of USV. To solve problem (4), the auction algorithm [27] is employed, which iteratively alternates between a bidding phase and an assignment phase. In the bidding phase, for each unassigned $i$-th USV, i.e., ${\displaystyle {\sum }_{j=1}^n}{\alpha }_{ij}=0$, the reward function is defined as$d_{ij}:=-\parallel {\mathit{r}}_{ij}(t_h)\parallel -z_j\mathrm{,}$(5)where $z_j$ denotes the current price of $j$-UAV, initialized as $0$. The optimal and suboptimal UAVs for the $i$-th USV are determined as$j^{*}=\arg \underset{j}{{\displaystyle \max }}{d}_{ij}\mathrm{,}{j}^{\prime }=\arg \underset{j\ne j^{*}}{{\displaystyle \max }}{d}_{ij}\mathrm{.}$(6)Accordingly, the $i$-th USV submits a bid to $j^{*}$-th UAV given by$o_{i{j}^{*}}=z_{j^{*}}+d_{i{j}^{*}}-d_{i{j}^{\prime }}+ϵ\mathrm{,}$(7)where $ϵ\in {ℝ}^{+}$ is a small positive constant. In the assignment phase, each UAV is allocated to the USV offering the highest bid, i.e.,$i^{*}=\arg \underset{i}{{\displaystyle \max }}{o}_{ij}\mathrm{,}$(8)the price of the $j$-th UAV is then updated as $z_j=o_{i^{*}j}$. If the $j$-th UAV was previously assigned to another USV $i\ne i^{*}$, the earlier assignment is canceled, i.e., ${\alpha }_{ij}=0$, and the new allocation is established with ${\alpha }_{i^{*}j}=1$.

Remark 1. The assignment problem in Eq. (4) imposes the one-to-one matching constraints in Eq. (4a), which ensures that each USV is assigned to exactly one UAV at each observation time $t_h$. When the numbers of UAVs and USVs are equal, these constraints define a matching between the two sets. Combined with the auction-based solution procedure, which iteratively assigns all remaining unassigned USVs, the proposed method guarantees that all USVs are observed at each observation time [27].

Define the total observation time as $\tilde{H}$, the observation number of the $i$-th USV by the $j$-th UAV as ${\tilde{h}}_{ij}$, $\tilde{H}={\displaystyle {\sum }_{j=1}^n}{\tilde{h}}_{ij}$. Building on the allocation matrix $\alpha$ obtained from Eq. (4), the estimated position of the $i$-th USV by its assigned the $j$-th UAV at time $t_h$ is described as${\tilde{p}}_i(t_h)=\left[{\tilde{p}}_{i[x]}(t_h)\mathrm{,}{\tilde{p}}_{i[y]}(t_h)\right]=\left[r_{ij[x]}(t_h)+q_{j[x]}(t_h)\mathrm{,}{r}_{ij[y]}(t_h)+q_{j[y]}(t_h)\right]\mathrm{.}$(9)Thereby, the trajectory of the $i$-th USV can be expressed as ${\tilde{s}}_i=\left[{\mathit{s}}_{i[x]}\mathrm{,}{\mathit{s}}_{i[y]}\right]={\left[{\tilde{p}}_i{(t_1)}^{\text{T}}\mathrm{,}{\tilde{p}}_i{(t_2)}^{\text{T}}\mathrm{,}\dots \mathrm{,}{\tilde{p}}_i{(t_{\tilde{H}})}^{\text{T}}\right]}^{\text{T}}\in {ℝ}^{\tilde{H}\times 2}$, and the velocity ${\tilde{v}}_i=\left[{\mathit{v}}_{i[x]}\mathrm{,}{v}_{i[y]}\right]\in {ℝ}^{\tilde{H}\times 2}$ is approximated by using the Euler method.

Coordinated dynamics with efficient block sparse Bayesian learning

An EBSBL is proposed to identify the coordinated dynamics of USVs. To approximate the unknown velocity function ${\mathit{f}}_i({\tilde{p}}_1(t)\mathrm{,}{\tilde{p}}_2(t)\mathrm{,}\dots \mathrm{,}{\tilde{p}}_n(t))$, we build up the vector of candidate functions $\mathit{\varphi }({\tilde{p}}_1(t)\mathrm{,}{\tilde{p}}_2(t)\mathrm{,}\dots \mathrm{,}{\tilde{p}}_n(t))\in {ℝ}^g$ composed of nonlinear candidate functions, where $g$ denotes the number of functions. Define ${\mathit{\Phi }}_i\in {ℝ}^{H\times g}$ as follows:${\mathit{\Phi }}_i:={\left[\mathit{\varphi }({\tilde{p}}_1(t_1)\mathrm{,}\dots \mathrm{,}{\tilde{p}}_n(t_1))\mathrm{,}\dots \mathrm{,}\mathit{\varphi }({\tilde{p}}_1(t_H)\mathrm{,}\dots \mathrm{,}{\tilde{p}}_n(t_H))\right]}^{\text{T}}\mathrm{,}$(10)where $H$ denotes the current observation number. Define the set of time observations associated with the $j$-th UAV of the $i$-th USV as ${\mathcal{I}}_{ij}$, the vector of weights to be identified as ${\tilde{\theta }}_{ij}=\left[{\mathit{\theta }}_{ij[x]}\mathrm{,}{\mathit{\theta }}_{ij[x]}\right]\in {ℝ}^{g\times 2}$, ${\mathit{v}}_{ij[x]}={\mathit{v}}_{i[x]}\left[{\mathcal{I}}_{ij}\mathrm{,}:\right]\in {ℝ}^{h_{ij}}$, ${\mathit{v}}_{ij[y]}={\mathit{v}}_{i[y]}\left[{\mathcal{I}}_{ij}\mathrm{,}:\right]\in {ℝ}^{h_{ij}}$, ${\mathit{\Phi }}_{ij}={\mathit{\Phi }}_i\left[{\mathcal{I}}_{ij}\mathrm{,}:\right]\in {ℝ}^{h_{ij}\times g}$, one has$\left[{\mathit{v}}_{ij[x]}\mathrm{,}{\mathit{v}}_{ij[y]}\right]={\mathit{\Phi }}_{ij}{\tilde{\theta }}_{ij}+\left[{\mathit{\lambda }}_{[x]}\mathrm{,}{\mathit{\lambda }}_{[y]}\right]\mathrm{,}{\mathit{\lambda }}_{[x]}\mathrm{,}{\mathit{\lambda }}_{[y]}\sim \mathcal{N}\left(\mathit{0}\mathrm{,}\sigma {\mathit{I}}_{h_{ij}}\right)\mathrm{,}$(11)where $h_{ij}$ denotes the current observation number of the $i$-th USV by the $j$-th UAV.

Since the proposed method independently identifies the coordinated dynamics of each USV in both $x$- and $y$- directions, the subscripts $i$, $[x]$, and $[y]$ are omitted for conciseness. For the $x$- and $y$-axes dynamics of the $i$-th USV, define data vector $\mathit{\tau }\in {ℝ}^H$ and dictionary matrix $\mathit{\Psi }\in {ℝ}^{H\times gn}$ stacked from all ${\mathit{v}}_{ij}$ and ${\mathit{\Phi }}_{ij}$, respectively, as follows:$\mathit{\tau }:={\left[{\mathit{v}}_{i1}^{\text{T}}\mathrm{,}{\mathit{v}}_{i2}^{\text{T}}\mathrm{,}\dots \mathrm{,}{\mathit{v}}_{in}^{\text{T}}\right]}^{\text{T}}\mathrm{,}\mathit{\Psi }:=\left[{\mathit{\Psi }}_{jg\text{'}}\right]\mathrm{,}$(12)and the elements of ${\mathit{\Psi }}_{jg\text{'}}\in {ℝ}^{h_{ij}\times n}$ are defined as follows:${\Psi }_{jg\text{'}\mathrm{,}pq}:=\{\begin{array}{ll}{\Phi }_{ij\mathrm{,}pg\text{'}}\mathrm{,}\hfill & \text{if}q=j\mathrm{,}\hfill \\ \mathrm{0,}\hfill & \text{otherwise}\mathrm{,}\hfill \end{array}$(13)where ${\Psi }_{jg\text{'}\mathrm{,}pq}$ and ${\Phi }_{ij\mathrm{,}pg\text{'}}$ denote the elements in the $p$-th row and $q$-th column of ${\mathit{\Psi }}_{jg\text{'}}$ and the $p$-th row and $g\text{'}$-th column of ${\mathit{\Phi }}_{ij}$, respectively, $p\in \left\{\mathrm{1,2,}\dots \mathrm{,}{h}_{ij}\right\}\mathrm{,}q\in \left\{\mathrm{1,2,}\dots \mathrm{,}n\right\}\mathrm{,}g\text{'}\in \left\{\mathrm{1,2,}\dots \mathrm{,}g\right\}$. Define $\mathit{w}\in {ℝ}^{gn}$ stacked from all ${\theta }_{ij}$ as follows:$\mathit{w}={\left[{\mathit{w}}_1\mathrm{,}{\mathit{w}}_2\mathrm{,}\dots \mathrm{,}{\mathit{w}}_g\right]}^T\mathrm{,}{\mathit{w}}_{g\text{'}}=\left[{\theta }_{i\mathrm{1,}g\text{'}}\mathrm{,}{\theta }_{i\mathrm{2,}g\text{'}}\mathrm{,}\dots \mathrm{,}{\theta }_{in\mathrm{,}g\text{'}}\right]\in {ℝ}^n\mathrm{,}$(14)where ${\theta }_{ij\mathrm{,}g\text{'}}$ denotes the $g\text{'}$-th element of ${\mathit{\theta }}_{ij}$. As a result, one has$\mathit{\tau }=\mathit{\Psi }\mathit{w}+\mathit{\epsilon }\mathrm{,}\mathit{\epsilon }\sim \mathcal{N}\left(\mathit{0}\mathrm{,}\sigma {\mathit{I}}_H\right)\mathrm{.}$(15)

Block prior is introduced as follows:$p\left(\mathit{w}|\mathit{\gamma }\right)={\displaystyle \prod _{g\text{'}=1}^g}p\left({\mathit{w}}_{g\text{'}}|{\mathit{\gamma }}_{g\text{'}}\right)\mathrm{,}p\left({\mathit{w}}_{g\text{'}}|{\mathit{\gamma }}_{g\text{'}}\right)={\displaystyle \prod _{j=1}^n}\mathcal{N}\left(w_{g\text{'}\mathrm{,}j}|\mathrm{0,}{\gamma }_{g\text{'}}\right)\mathrm{,}$(16)where $w_{g\text{'}\mathrm{,}j}$ denotes the $j$-th element of ${\mathit{w}}_{g\text{'}}$, $\mathit{\gamma }:={\left[{\mathit{\gamma }}_1\mathrm{,}{\mathit{\gamma }}_2\mathrm{,}\dots \mathrm{,}{\mathit{\gamma }}_g\right]}^{\text{T}}\in {ℝ}^{gn}$, ${\mathit{\gamma }}_{g\text{'}}:=\left[{\gamma }_{g\text{'}}\mathrm{,}\dots \mathrm{,}{\gamma }_{g\text{'}}\right]\in {ℝ}^{1\times n}$. Define an auxiliary variable $\mathit{\varpi }\in {ℝ}^{gn}$, and the likelihood function $p(\mathit{\tau }|\mathit{w}\mathrm{,}\sigma )$ can be written as [28]$p(\mathit{\tau }|\mathit{w}\mathrm{,}\sigma )={\left(2\text{π}\sigma \right)}^{-\frac{H}{2}}\exp \left(-{\left(2\sigma \right)}^{(-1)}\parallel \mathit{\tau }-\mathit{\Psi }\mathit{w}{\parallel }^2\right)=\underset{\varpi }{{\displaystyle \max }}{p}^{\prime }(\mathit{\tau }|\mathit{w}\mathrm{,}\sigma ;\mathit{\varpi })\mathrm{,}$(17)where$\begin{array}{l}{p}^{\prime }(\mathit{\tau }|\mathit{w}\mathrm{,}\sigma ;\mathit{\varpi }):={\left(2\text{π}\sigma \right)}^{-\frac{H}{2}}\exp \left(-{\left(2\sigma \right)}^{(-1)}R(\mathit{w}\mathrm{,}\mathit{\varpi })\right)\mathrm{,}\\ R(\mathit{w}\mathrm{,}\mathit{\varpi }):=\parallel \mathit{\tau }-\mathit{\Psi }\mathit{\varpi }{\parallel }^2+2(\mathit{w}-\mathit{\varpi }{)}^{\text{T}}{\mathit{\Psi }}^{\text{T}}(\mathit{\Psi }\mathit{\varpi }-\mathit{\tau })+\beta \parallel \mathit{w}-\mathit{\varpi }{\parallel }^2\mathrm{,}\end{array}$(18)$\beta =\text{eig}({\mathit{\Psi }}^{\text{T}}\mathit{\Psi })+\zeta$, $\zeta \in {ℝ}^{+}$ is a small positive constant, and $\text{eig}(\cdot )$ denotes eigenvalues. We use the strict lower bound function $p^{\prime }(\mathit{\tau }|\mathit{w}\mathrm{,}\sigma ;\widehat{\varpi })$ of the likelihood function $p(\mathit{\tau }|\mathit{w}\mathrm{,}\sigma )$ to compute the posterior distribution of $\mathit{w}$, as follows:$p(\mathit{w}|\mathit{\tau }\mathrm{,}\sigma \mathrm{,}\mathit{\gamma })\approx \frac{p^{\prime }(\mathit{\tau }|\mathit{w}\mathrm{,}\sigma ;\widehat{\varpi })p(\mathit{w}|\mathit{\gamma })}{{\displaystyle \int }{p}^{\prime }(\mathit{\tau }|\mathit{w}\mathrm{,}\sigma ;\widehat{\varpi })p(\mathit{w}|\mathit{\gamma })\text{d}\mathit{w}}\sim \mathcal{N}\left(\mathit{w}|{\mathit{\mu }}_{\mathit{w}}\mathrm{,}{\mathit{\Sigma }}_{\mathit{w}}\right)\mathrm{,}$(19)where $\widehat{\varpi }$ is the estimated fixed vector, and the posterior covariance ${\mathit{\Sigma }}_{\mathit{w}}\in {ℝ}^{gn\times gn}$ and mean ${\mathit{\mu }}_{\mathit{w}}\in {ℝ}^{gn}$ of $\mathit{w}$ are given by$\begin{array}{l}{\mathit{\Sigma }}_{\mathit{w}}={\left({\mathit{\Gamma }}^{-1}+\beta {\sigma }^{-1}{\mathit{I}}_{gn}\right)}^{-1}\mathrm{,}\mathit{\Gamma }=\text{diag}(\mathit{\gamma })\mathrm{,}\\ {\mathit{\mu }}_{\mathit{w}}={\sigma }^{-1}{\mathit{\Sigma }}_{\mathit{w}}\left(\beta \widehat{\varpi }-{\mathit{\Psi }}^{{\rm T}}\mathit{\Psi }\widehat{\varpi }+{\mathit{\Psi }}^{\text{T\hspace{0.05em}}}\mathit{\tau }\right)\mathrm{.}\end{array}$(20)The purpose is to estimate the unknown parameters $\mathit{\varpi }$, $\mathit{\gamma }$, and $\sigma$ using the evidence maximization method [28], the optimal values of $\mathit{\gamma }$ and $\sigma$ are obtained by maximizing the marginalized probability density function $p(\mathit{\tau }\mathrm{,}\sigma \mathrm{,}\mathit{\gamma })$ as follows:$\begin{array}{c}\left({\mathit{\gamma }}^{*}\mathrm{,}{\sigma }^{*}\right)=\arg \underset{\sigma \mathrm{,}\mathit{\gamma }}{{\displaystyle \max }}p(\mathit{\tau }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )\mathrm{,}\\ p(\mathit{\tau }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )={\displaystyle \int }p(\mathit{\tau }|\mathit{w}\mathrm{,}\sigma )p(\mathit{w}|\mathit{\gamma })\text{d}\mathit{w}={\displaystyle \int }\underset{\mathit{\varpi }}{{\displaystyle \text{max}}}{p}^{\prime }(\mathit{\tau }|\mathit{w}\mathrm{,}\sigma ;\mathit{\varpi })p(\mathit{w}|\mathit{\gamma })\text{d}\mathit{w}\\ \ge \underset{\mathit{\varpi }}{{\displaystyle \text{max}}}{\displaystyle \int }{p}^{\prime }(\mathit{\tau }|\mathit{w}\mathrm{,}\sigma ;\mathit{\varpi })p(\mathit{w}|\mathit{\gamma })\text{d}\mathit{w}\mathrm{,}\end{array}$(21)where the last inequality is obtained by swapping the order of integration and maximization [28]. As a result,$\begin{array}{c}({\mathit{\varpi }}^{*}\mathrm{,}{\mathit{\gamma }}^{*}\mathrm{,}{\sigma }^{*})=\text{arg}\underset{\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma }{{\displaystyle \text{max}}}{\displaystyle \int }{p}^{\prime }(\mathit{\tau }|\mathit{w}\mathrm{,}\sigma ;\mathit{\varpi })p(\mathit{w}|\mathit{\gamma })\text{d}\mathit{w}\\ =\text{arg}\underset{\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma }{{\displaystyle \text{max}}}{\left(2\text{π}\sigma \right)}^{-\frac{H}{2}}{\left|\mathit{\Gamma }{\mathit{\Sigma }}_{\mathit{w}}^{-1}\right|}^{-\frac{12}{}}\exp \left\{-\frac{12}{}J(\mathit{\varpi }\mathrm{,}\sigma \mathrm{,}\mathit{\gamma })\right\}\mathrm{,}\end{array}$(22)with$J(\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )={\sigma }^{-1}\left(\parallel \mathit{\tau }-\mathit{\Psi }\mathit{\varpi }{\parallel }^2-2{\mathit{\varpi }}^{\text{T}}{\mathit{\Psi }}^{\text{T}}(\mathit{\Psi }\mathit{\varpi }-\mathit{\tau })+\beta {\mathit{\varpi }}^{\text{T}}\mathit{\varpi }\right)-{\mathit{\mu }}_{\mathit{w}}^{\text{T}}{\mathit{\Sigma }}_{\mathit{w}}^{-1}{\mathit{\mu }}_{\mathit{w}}\mathrm{.}$(23)Taking $-2\ln (\cdot )$ of Eq. (22), we obtain the following objective function to be minimized:$\mathcal{L}(\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )=n_1\ln {\sigma }^{-1}+{\displaystyle \sum _{g^{\prime }=1}^g}n\ln \left(\sigma +\beta {\gamma }_{g^{\prime }}\right)+J(\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )\mathrm{,}$(24)with $n_1=gn-H$. From Eq. (20), one has $\sigma {\mathbf{\text{Σ}}}_{\mathit{w}}^{-1}{\mathit{\mu }}_{\mathit{w}}=\beta \mathit{\varpi }-{\mathit{\Psi }}^{\text{T}}\mathit{\Psi }\mathit{\varpi }+{\mathit{\Psi }}^{\text{T}}\mathit{\tau }\mathrm{.}$ As a result,$J(\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )={\sigma }^{-1}\left(\parallel \mathit{\tau }-\mathit{\Psi }\mathit{\varpi }{\parallel }^2+2{\left({\mathit{\mu }}_{\mathit{w}}-\mathit{\varpi }\right)}^{\text{T}}{\mathit{\Psi }}^{\text{T}}(\mathit{\Psi }\mathit{\varpi }-\mathit{\tau })+\beta {\Vert {\mathit{\mu }}_{\mathit{w}}-\mathit{\varpi }\Vert }^2\right)+{\mathit{\mu }}_{\mathit{w}}^{\text{T}}{\mathit{\Gamma }}^{-1}{\mathit{\mu }}_{\mathit{w}}\mathrm{.}$(25)Note that ${\mathit{w}}^{\star }={\mathit{\mu }}_{\mathit{w}}$, one has$J(\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )=\underset{\mathit{w}}{{\displaystyle \min }}{\sigma }^{-1}\mathit{R}(\mathit{w}\mathrm{,}\mathit{\varpi })+{\mathit{w}}^{\text{T}}{\mathit{\Gamma }}^{-1}\mathit{w}\mathrm{.}$(26)Therefore, the joint objective function is obtained as follows:$\mathcal{L}(\mathit{w}\mathrm{,}\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )={\sigma }^{-1}\mathit{R}(\mathit{w}\mathrm{,}\mathit{\varpi })-n_1\ln \sigma +{\displaystyle \sum _{g\text{'}=1}^g}\left[n\ln \left(\sigma +\beta {\gamma }_{g\text{'}}\right)+\frac{{\displaystyle {\sum }_{j=1}^n{w}_{g\text{'}\mathrm{,}j}^2}}{{\gamma }_{g\text{'}}}\right]=g(\mathit{w}\mathrm{,}\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )+f(\mathit{\gamma }\mathrm{,}\sigma )\mathrm{,}$(27)with$\begin{array}{c}g(\mathit{w}\mathrm{,}\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma ):={\displaystyle \sum _{g\text{'}=1}^g}\left[\frac{{\displaystyle {\sum }_{j=1}^n{w}_{g\text{'}\mathrm{,}j}^2}}{{\gamma }_{g\text{'}}}\right]-n_1\ln \sigma +{\sigma }^{-1}R(\mathit{w}\mathrm{,}\mathit{\varpi })\mathrm{,}\\ f(\mathit{\gamma }\mathrm{,}\sigma ):={\displaystyle \sum _{g\text{'}=1}^g}n\ln \left(\sigma +\beta {\gamma }_{g\text{'}}\right)\mathrm{.}\end{array}$(28)Note that $g(\mathit{w}\mathrm{,}\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )$ is convex with respect to $\left\{\mathit{w}\mathrm{,}\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma \right\}$, and $f(\mathit{\gamma }\mathrm{,}\sigma )$ is concave with respect to $\left\{\mathit{\gamma }\mathrm{,}\sigma \right\}$. Hence, $\mathcal{L}(\mathit{w}\mathrm{,}\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )$ is a convex-concave procedure problem [29], which can be solved as follows:${\mathit{w}}^{(k+1)}=\text{arg}\underset{\mathit{w}}{{\displaystyle \text{min }}}g(\mathit{w}\mathrm{,}{\mathit{\varpi }}^{(k)}\mathrm{,}{\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)})\mathrm{,}$(29a)${\mathit{\varpi }}^{(k+1)}=\text{arg}\underset{\mathit{\varpi }}{{\displaystyle \text{min }}}g({\mathit{w}}^{(k+1)}\mathrm{,}\mathit{\varpi }\mathrm{,}{\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)})\mathrm{,}$(29b)${\sigma }^{(k+1)}=\text{arg}\underset{\sigma }{{\displaystyle \text{min}}}g({\mathit{w}}^{(k+1)}\mathrm{,}{\mathit{\varpi }}^{(k+1)}\mathrm{,}{\mathit{\gamma }}^{(k)}\mathrm{,}\sigma )+\langle \sigma -{\sigma }^{(k)}\mathrm{,}\nabla f\left({\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\rangle \mathrm{,}$(29c)${\mathit{\gamma }}^{(k+1)}=\text{arg}\underset{\gamma }{{\displaystyle \text{min }}}g({\mathit{w}}^{(k+1)}\mathrm{,}{\mathit{\varpi }}^{(k+1)}\mathrm{,}\mathit{\gamma }\mathrm{,}{\sigma }^{(k+1)})+\langle \mathit{\gamma }-{\mathit{\gamma }}^{(k)}\mathrm{,}\nabla f\left({\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\rangle \mathrm{.}$(29d)Since the objective functions in Eqs. (29a)–(29c) are convex, setting the gradient to zero yields:${\mathit{\Sigma }}_{\mathit{w}}^{(k+1)}={\left({\mathit{\Gamma }}^{{(k)}^{-1}}+\frac{\beta }{{\sigma }^{(k)}}{\mathit{I}}_{gn}\right)}^{-1}\mathrm{,}$(30a)${\mathit{w}}^{(k+1)}=\frac{1}{{\sigma }^{(k)}}{\mathit{\Sigma }}_{\mathit{w}}^{(k+1)}\left(\beta {\mathit{\varpi }}^{(k)}-{\mathit{\Psi }}^{\text{T\hspace{0.05em}}}\mathit{\Psi }{\mathit{\varpi }}^{(k)}+{\mathit{\Psi }}^{\text{T\hspace{0.05em}}}\mathit{\tau }\right)\mathrm{,}$(30b)${\mathit{\varpi }}^{(k+1)}={\mathit{w}}^{(k+1)}\mathrm{,}$(30c)${\varrho }^{(k+1)}=n{\displaystyle \sum _{g\text{'}=1}^g}{\left({\sigma }^{(k)}+\beta {\gamma }_{g\text{'}}^{(k)}\right)}^{-1}\mathrm{,}$(30d)${\sigma }^{(k+1)}=\frac{n_1+\sqrt{n_{12}^{}+4{\varrho }^{(k+1)}\left(\parallel \mathit{\tau }-\mathit{\Psi }{\mathit{\varpi }}^{(k+1)}{\parallel }^2\right)}}{2{\varrho }^{(k+1)}}\mathrm{,}$(30e)${\gamma }_{g\text{'}}^{(k+1)}={\left(\frac{{\sigma }^{(k)}+\beta {\gamma }_{g\text{'}}^{(k)}}{n\beta }{\displaystyle \sum _{j=1}^n}{\left(w_{g\text{'}\mathrm{,}j}^{(k+1)}\right)}^2\right)}^{\frac{12}{}}\mathrm{.}$(30f)The final $\widehat{\theta }\in {ℝ}^g$ is obtained by averaging over each block ${\widehat{w}}_{g\text{'}}$ in $\widehat{w}$ as follows:$\widehat{\theta }=[{\widehat{\theta }}_1\mathrm{,}{\widehat{\theta }}_2\mathrm{,}\dots \mathrm{,}{\widehat{\theta }}_g]\mathrm{,}{\widehat{\theta }}_{g\text{'}}=\frac{1}{n}{\displaystyle \sum _{j=1}^n}{\widehat{w}}_{g\text{'}\mathrm{,}j}\mathrm{.}$(31)

Note that Eq. (30a) only involves the inversion of a diagonal matrix, which has an operation of $\mathcal{O}(n)$. Consider the matrix multiplication ${\mathit{\Psi }}^{\text{T\hspace{0.05em}}}\mathit{\Psi }\mathit{\varpi }$ in Eq. (30b), the proposed EBSBL has a computational complexity of $\mathcal{O}(n^2)$. In contrast, for conventional block sparse Bayesian learning algorithm [30], each iteration requires computing the inverse of a non-diagonal matrix, resulting in $\mathcal{O}(n^3)$ computational complexity. As a result, EBSBL reduces the computational complexity from $\mathcal{O}(n^3)$ to $\mathcal{O}(n^2)$.

Then, the coordinated dynamics ${\mathit{f}}_i({\tilde{p}}_1(t)\mathrm{,}{\tilde{p}}_2(t)\mathrm{,}\dots \mathrm{,}{\tilde{p}}_n(t))$ of the $i$-th USV can be identified as follows:${\mathit{f}}_i({\tilde{p}}_1(t)\mathrm{,}{\tilde{p}}_2(t)\mathrm{,}\dots \mathrm{,}{\tilde{p}}_n(t))={\mathit{\Phi }}_i\left[{\mathit{\theta }}_{i[x]}\mathrm{,}{\mathit{\theta }}_{i[y]}\right]\mathrm{.}$(32)

Theorem 1. The sequence ${\left\{\mathcal{L}\left({\mathit{w}}^{(k)}\mathrm{,}{\mathit{\varpi }}^{(k)}\mathrm{,}{\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\right\}}_{k=0}^{\infty }$ generated using EBSBL is non-increasing and locally convergent. Moreover, $\left\{{\mathit{w}}^{(k)}\mathrm{,}{\mathit{\varpi }}^{(k)}\mathrm{,}{\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right\}$ is bounded.

Proof. Define the surrogate function ${\tilde{\mathcal{L}}}_k(\mathit{w}\mathrm{,}\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )$ as follows:${\tilde{\mathcal{L}}}_k(\mathit{w}\mathrm{,}\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma ):=g(\mathit{w}\mathrm{,}\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )+\langle \left(\mathit{\gamma }\mathrm{,}\sigma \right)-\left({\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\mathrm{,}\nabla f\left({\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\rangle +f\left({\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\mathrm{.}$(33)As $\langle \left(\mathit{\gamma }\mathrm{,}\sigma \right)-\left({\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\mathrm{,}\nabla f\left({\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\rangle +f\left({\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)$ is an affine function and $g(\mathit{w}\mathrm{,}\mathit{\varpi }\mathrm{,}\mathit{\gamma }\mathrm{,}\sigma )$ is convex, one has$\begin{array}{l}{\tilde{\mathcal{L}}}_k\left({\mathit{w}}^{(k+1)}\mathrm{,}{\mathit{\varpi }}^{(k+1)}\mathrm{,}{\mathit{\gamma }}^{(k+1)}\mathrm{,}{\sigma }^{(k+1)}\right)\le {\tilde{\mathcal{L}}}_k\left({\mathit{w}}^{(k+1)}\mathrm{,}{\mathit{\varpi }}^{(k+1)}\mathrm{,}{\mathit{\gamma }}^{(k+1)}\mathrm{,}{\sigma }^{(k)}\right)\le {\tilde{\mathcal{L}}}_k\left({\mathit{w}}^{(k+1)}\mathrm{,}{\mathit{\varpi }}^{(k+1)}\mathrm{,}{\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\\ \le {\tilde{\mathcal{L}}}_k\left({\mathit{w}}^{(k+1)}\mathrm{,}{\mathit{\varpi }}^{(k)}\mathrm{,}{\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\le {\tilde{\mathcal{L}}}_k\left({\mathit{w}}^{(k)}\mathrm{,}{\mathit{\varpi }}^{(k)}\mathrm{,}{\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)=\mathcal{L}\left({\mathit{w}}^{(k)}\mathrm{,}{\mathit{\varpi }}^{(k)}\mathrm{,}{\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\mathrm{.}\end{array}$(34)The concavity of $f(\mathit{\gamma }\mathrm{,}\sigma )$ leads to$\mathcal{L}\left({\mathit{w}}^{(k+1)}\mathrm{,}{\mathit{\varpi }}^{(k+1)}\mathrm{,}{\mathit{\gamma }}^{(k+1)}\mathrm{,}{\sigma }^{(k+1)}\right)\le {\tilde{\mathcal{L}}}_k\left({\mathit{w}}^{(k+1)}\mathrm{,}{\mathit{\varpi }}^{(k+1)}\mathrm{,}{\mathit{\gamma }}^{(k+1)}\mathrm{,}{\sigma }^{(k+1)}\right)\le \mathcal{L}\left({\mathit{w}}^{(k)}\mathrm{,}{\mathit{\varpi }}^{(k)}\mathrm{,}{\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\mathrm{.}$(35)Therefore, the sequence ${\left\{\mathcal{L}\left({\mathit{w}}^{(k)}\mathrm{,}{\mathit{\varpi }}^{(k)}\mathrm{,}{\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\right\}}_{k=0}^{\infty }$ is non-increasing. Since $\mathit{\gamma }>0$, $\sigma >0$, and $R(\mathit{w}\mathrm{,}\mathit{\varpi })>0$, the cost function $\mathcal{L}\left(w\mathrm{,}\varpi \mathrm{,}\gamma \mathrm{,}\sigma \right)$ is lower bounded. By the monotone convergence theorem [31], the nonincreasing sequence ${\left\{\mathcal{L}\left({\mathit{w}}^{(k)}\mathrm{,}{\mathit{\varpi }}^{(k)}\mathrm{,}{\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right)\right\}}_{k=0}^{\infty }$ is locally convergent. Hence, the sequence $\left\{{\mathit{w}}^{(k)}\mathrm{,}{\mathit{\varpi }}^{(k)}\mathrm{,}{\mathit{\gamma }}^{(k)}\mathrm{,}{\sigma }^{(k)}\right\}$ is bounded, which completes the proof.

Cross-domain coordinated tracking

The UKF [32] is employed in this system to predict the position of USVs by handling the nonlinear dynamics. Unlike EKF, the UKF does not require linearization, making it more suitable for complex systems. UKF generates sigma points around the current state estimate and propagates them through the nonlinear model, providing more accurate state and covariance estimates. While direct trajectory estimation based on USV dynamics does not account for uncertainties such as sensor inaccuracies and environmental disturbances, UKF integrates the dynamics model with measurements in a probabilistic framework. It iteratively updates the trajectory estimate using the relative position data observed by UAVs, correcting the estimate at each time step based on the new measurement and the predicted state from the previous step.

Let ${\widehat{p}}_{t|t-\Delta t}$ and ${\widehat{p}}_{t-\Delta t|t-\Delta t}$ denote the posterior estimated position and the estimated position from the previous step of $i$-th USV, respectively, where $\Delta t$ represents the time step. The predicted position ${\widehat{p}}_{t|t-\Delta t}$ of the $i$-th USV at the current time $t$ is obtained from the previous state estimate ${\widehat{p}}_{t-\Delta t|t-\Delta t}$ and the system’s dynamic model ${\mathit{f}}_i$. The UKF predicts the position at the current time step as follows:${\widehat{p}}_{t|t-\Delta t}={\widehat{p}}_{t-\Delta t|t-\Delta t}+{\mathit{f}}_i({\tilde{p}}_1(t)\mathrm{,}{\tilde{p}}_2(t)\mathrm{,}\dots \mathrm{,}{\tilde{p}}_n(t))\Delta t\mathrm{.}$(36)UKF generates a set of sigma points ${\mathcal{X}}_b$ to approximate the probability distribution of the system state, which are derived from the current state estimate ${\widehat{p}}_{t|t}$ and the associated covariance matrix, describing the uncertainty in the current state estimate. The sigma points ${\mathcal{X}}_b$ are generated and propagated through the system model as follows:${\mathcal{X}}_0={\widehat{p}}_{t|t}\mathrm{,}{\mathcal{X}}_b={\widehat{p}}_{t|t}+\sqrt{2+\xi }{\mathit{S}}_b\mathrm{,}{\mathcal{X}}_{b+2}={\widehat{p}}_{t|t}-\sqrt{2+\xi }{\mathit{S}}_b\mathrm{,}b=\mathrm{1,2,}$(37)where $\xi$ is a scaling parameter, and ${\mathit{S}}_b$ the column vectors of the square root of the covariance matrix ${\mathit{P}}_{t|t-\Delta t}$. The predicted states are yielded as follows:${\widehat{p}}_{t+\Delta t|t}={\displaystyle \sum _{b=0}^4}{\mathit{W}}_b^{(m)}{\mathcal{X}}_b^{-}\mathrm{,}$(38)where ${\mathcal{X}}_b^{-}$ denotes the sigma points propagated through ${\mathit{f}}_i$ over $\Delta t$, and ${\mathit{W}}_b^{(m)}$ are the weights associated with each sigma point. The predicted covariance is given by${\mathit{P}}_{t+\Delta t|t}={\displaystyle \sum _{b=0}^4}{\mathit{W}}_b^{(c)}({\mathcal{X}}_b^{-}-{\widehat{p}}_{t+\Delta t|t})({\mathcal{X}}_b^{-}-{\widehat{p}}_{t+\Delta t|t}{)}^{\text{T}}+\mathit{Q}\mathrm{,}$(39)where ${\mathit{W}}_b^{(c)}$ are the covariance weights associated with each sigma point, $\mathit{Q}\in {ℝ}^{2\times 2}$ represents the external noise covariance. The new sigma points ${\widehat{\mathcal{X}}}_b$ are generated using the updated ${\widehat{p}}_{t+\Delta t|t}$ and ${\mathit{P}}_{t+\Delta t|t}$. The measurement update step involves generating predicted measurements ${\mathcal{Z}}_i$ from the predicted sigma points, using the measurement function $h({\tilde{p}}_{t|t})$${\mathcal{Z}}_b=h({\widehat{\mathcal{X}}}_b)\mathrm{,}{\widehat{z}}_{t+\Delta t}={\displaystyle \sum _{b=0}^4}{\mathit{W}}_b^{(m)}{\mathcal{Z}}_b\mathrm{.}$(40)The measurement covariance S and cross-covariance C are given by$\begin{array}{l}\mathit{S}={\displaystyle \sum _{b=0}^4}{\mathit{W}}_b^{(c)}({\mathcal{Z}}_b-{\widehat{z}}_{t+\Delta t})({\mathcal{Z}}_b-{\widehat{z}}_{t+\Delta t}{)}^{\text{T}}+\mathit{D}\mathrm{,}\\ \mathit{C}={\displaystyle \sum _{b=0}^4}{\mathit{W}}_b^{(c)}({\widehat{\mathcal{X}}}_b-{\widehat{p}}_{t+\Delta t|t})({\mathcal{Z}}_b-{\widehat{z}}_{t+\Delta t}{)}^{\text{T}}\mathrm{,}\end{array}$(41)where $\mathit{D}\in {ℝ}^{2\times 2}$ is the measurement noise covariance. The Kalman gain in UKF is computed as ${\mathit{K}}_{t+\Delta t}=\mathit{C}{\mathit{S}}^{-1}$. The current state estimate ${\widehat{p}}_{t+\Delta t|t+\Delta t}$ and updated covariance ${\mathit{P}}_{t+\Delta t|t+\Delta t}$ are updated as follows:$\begin{array}{l}{\widehat{p}}_{t+\Delta t|t+\Delta t}={\widehat{p}}_{t+\Delta t|t}+{\mathit{K}}_{t+\Delta t}({\tilde{p}}_{t+\Delta t}-{\widehat{z}}_{t+\Delta t})\mathrm{,}\\ {\mathit{P}}_{t+\Delta t|t+\Delta t}={\mathit{P}}_{t+\Delta t|t}-{\mathit{K}}_{t+\Delta t}\mathit{S}{\mathit{K}}_{t+\Delta t}^{\text{T}}\mathrm{.}\end{array}$(42)

The controller for the $j$-th UAV to monitor the $i$-th USV is designed as follows [33]:${\mathit{u}}_j^{\text{init}}(t)=k_p{\mathit{e}}_p(t)+k_d{\mathit{e}}_d(t)\mathrm{,}$(43)where$\begin{array}{c}{\mathit{e}}_p(t)={\left[{\widehat{p}}_{i[x]}(t)-q_{j[x]}(t)\mathrm{,}{\widehat{p}}_{i[y]}(t)-q_{j[y]}(t)\mathrm{,0}\right]}^{\text{T}}\mathrm{,}\\ {\mathit{e}}_d(t)={\left[{\mathit{f}}_i({\widehat{p}}_1(t)\mathrm{,}\dots \mathrm{,}{\widehat{p}}_n(t)\mathrm{,0}\right]}^{\text{T}}-{\mathit{\kappa }}_j(t)\mathrm{,}\end{array}$(44)$k_p\in {ℝ}^{+}$ and $k_d\in {ℝ}^{+}$ are the proportional and derivative gains, respectively. To avoid collisions among UAVs, a potential field-based controller is utilized for UAVs. The repulsive force ${\mathit{F}}_{jk}(t)$ exerted on the $j$-th UAV by $k$-th UAV is defined as follows [34]:${\mathit{F}}_{jk}(t)=\{\begin{array}{ll}{k}_r\frac{{\mathit{q}}_j(t)-{\mathit{q}}_k(t)}{\parallel {\mathit{q}}_j(t)-{\mathit{q}}_k(t){\parallel }^2}\mathrm{,}\hfill & \parallel {\mathit{q}}_j(t)-{\mathit{q}}_k(t)\parallel <d_s\mathrm{,}\hfill \\ \mathrm{0,}\hfill & \parallel {\mathit{q}}_j(t)-{\mathit{q}}_k(t)\parallel \ge d_s\mathrm{,}\hfill \end{array}$(45)where $k_r\in {ℝ}^{+}$ is the repulsive force coefficient. According to Ref. [35], the controller ${\mathit{u}}_j^{\text{avoid}}(t)$ of the $j$-th UAV is${\mathit{u}}_j^{\text{avoid}}(t)=-\eta {\displaystyle \sum _{k\ne j}}{\mathit{F}}_{jk}(t)\mathrm{,}$(46)where $\eta \in {ℝ}^{+}$ is a positive gain parameter. As a result, the control input of UAVs is as follows:${\mathit{u}}_j(t)={\mathit{u}}_j^{\text{init}}(t)+{\mathit{u}}_j^{\text{avoid}}(t)\mathrm{.}$(47)

The complete monitoring procedures are summarized in Algorithm 1 with the associated diagram shown in Figure 1.

NUMERICAL AND EXPERIMENTAL RESULTS

Setups

In this section, we demonstrate the effectiveness and superiority of the proposed MUMU-EBSBL by both numerical simulation and real-lake experiments. For comparison, we construct three variants by replacing the EBSBL module with mainstream system identification methods, i.e., block SBL (BSBL) [30], vanilla SBL (VSBL) [36], and LASSO [37], yielding MUMU-BSBL, MUMU-VSBL, and MUMU-LASSO, respectively. The MUMU framework and all the other settings are kept identical across variants. Circular formation [25] and line formation [8] are selected for the coordinated dynamics of USVs.

In numerical simulations, we consider four UAVs monitoring four USVs. In real-lake experiments, we introduce a self-developed cross-domain platform, including three HUSTER-12c USVs and three M-200 UAVs. As shown in Figure 2(a) and (b), the HUSTER-12c USV has a length of 1.2 m and a width of 0.42 m. It is equipped with two CA-6152A GPS antennas, an STM32F407 control module, and a TP-Link TLBS520 Wi-Fi module. Each M-200 UAV measures 0.65 m in both length and width, and is fitted with a DJI Matrice 200 Series GPS module, a Manifold 2 control module, and the same TP-Link TLBS520 Wi-Fi module. Figure 2(c) shows the coordination workflow of the cross-domain platform. UAVs track the positions of USVs through observation and establish communication using a WiFi 5G network, allowing them to generate the required guidance signals for navigation. The base station receives and logs all states, including positions, tracking errors, etc., transmitted over the WiFi 5G network. To quantify the performance of the algorithms, the error metric is defined as $e_i=\sqrt{{(q_{j[x]}(t)-p_{i[x]}(t))}^2+{(q_{j[y]}(t)-p_{i[y]}(t))}^2}$.

Figure 2 Architecture of the real-lake experimental platform. (a) HUSTER-12c USV, (b) M-200 UAV, and the detailed components. (c) Operation procedure of the cross-domain monitoring system, which consists of three HUSTER-12c USVs, three M-200 UAVs, and a WiFi 5G (TP-link TLBS520) wireless communication station. USVs and UAVs have independent communication networks, with USVs not sharing information with UAVs.

Numerical simulation results

Figure 3 shows the tracking errors evolution of numerical simulations with circular formation. Figure 3(a) presents the tracking error results for four USVs under the proposed MUMU-EBSBL algorithm. The errors decrease steadily over time, demonstrating the effectiveness of MUMU-EBSBL. Figure 3(b) compares performance of the four algorithms at $t=\mathrm{5,10,15}$ s. It is observed that all the algorithms exhibit decreasing tracking errors, whereas the proposed MUMU-EBSBL always achieves the best performance with a considerable margin. For $t=15$ s, the error of MUMU-EBSBL is below $0.5$ m, indicating satisfactory performance. For $t=10$ s, the average error of MUMU-EBSBL is reduced by $39.0\%$, $67.8\%$, and $68.6\%$ relative to MUMU-BSBL, MUMU-VSBL, and MUMU-LASSO, respectively.

Figure 3 The tracking errors evolution of numerical simulations with circular formation. (a) The tracking errors of four USVs under MUMU-EBSBL. (b) The errors comparison among MUMU-EBSBL, MUMU-BSBL, MUMU-VSBL, and MUMU-LASSO at $t=\mathrm{5,10,15}$ s.

Figure 4 shows the tracking errors evolution of numerical simulations with line formation, where Figure 4(a) presents the results under MUMU-EBSBL, demonstrating that UAVs successfully track USVs. Figure 4(b) compares performance of the four algorithms at $t=\mathrm{5,10,15}$ s. It is observed that MUMU-EBSBL always achieves the best performance among all the four algorithms. For $t=5$ s, the average error of MUMU-EBSBL is more accurate than all other algorithms with the reduction of $36.6\%$, $48.6\%$, and $54.9\%$, compared with MUMU-BSBL, MUMU-VSBL, and MUMU-LASSO, respectively. The effectiveness and superiority of the proposed MUMU-EBSBL are thus verified.

Figure 4 The tracking errors evolution of numerical simulations with line formation. (a) The tracking errors of four USVs under MUMU-EBSBL. (b) The errors comparison among MUMU-EBSBL, MUMU-BSBL, MUMU-VSBL, and MUMU-LASSO at $t=\mathrm{5,10,15}$ s.

Real-lake experimental results

Figure 5 shows the experimental snapshots and tracking errors evolution of real-lake experiments with circular formation. Figure 5(c) presents the tracking errors evolution for three USVs under MUMU-EBSBL. It is observed that tracking errors gradually decrease. Figure 5(d) shows tracking errors for the four algorithms at $t=\mathrm{10,20,30}$ s. The results indicate that MUMU-EBSBL consistently outperforms the other algorithms. At $t=20$ s, the error of MUMU-EBSBL is reduced by $56.5\%$, $67.3\%$, and $81.2\%$ relative to MUMU-BSBL, MUMU-VSBL, and MUMU-LASSO, respectively. Equivalently, for $t=30$ s, the error of MUMU-EBSBL is $22.6\%$, $14.1\%$, and $7.5\%$ of that of MUMU-BSBL, MUMU-VSBL, and MUMU-LASSO, respectively.

Figure 5 The experimental snapshots and the tracking errors evolution of real-lake experiments with circular formation, where blue circles denote USVs, yellow circles denote UAVs, and red circle denotes trajectory. (a) Initial scene with USVs at their starting positions. (b) UAVs performing real-time tracking and monitoring of the motion of USVs. (c) The tracking errors of three USVs under MUMU-EBSBL. (d) The tracking errors comparison among MUMU-EBSBL, MUMU-BSBL, MUMU-VSBL, and MUMU-LASSO at $t=\mathrm{10,20,30}$ s.

Figure 6 illustrates the experimental snapshots and tracking errors evolution of real-lake experiments with line formation, where Figure 6(c) presents the results under the proposed MUMU-EBSBL. It is observed that tracking errors gradually decrease, demonstrating the effectiveness of the proposed MUMU-EBSBL. Figure 6(d) shows the comparison of tracking errors among MUMU-EBSBL, MUMU-BSBL, MUMU-VSBL, and MUMU-LASSO. The results indicate that MUMU-EBSBL yields the best tracking performance among all four algorithms. For $t=10$ s, the error of MUMU-EBSBL is $65.5\%$, $46.5\%$, and $38.6\%$ of that of MUMU-BSBL, MUMU-VSBL, and MUMU-LASSO, respectively. In addition, the error of MUMU-EBSBL for $t=30$ s is reduced by $60.9\%$, $86.8\%$, and $91.5\%$ relative to MUMU-BSBL, MUMU-VSBL, and MUMU-LASSO, respectively. Both the effectiveness and superiority of MUMU-EBSBL are thus demonstrated.