Nomenclature

AOA

angle of arrival

DSSPA

distributed stochastic subgradient projection algorithm

GDOP

geometry dilution of precision

LOS

line of sight

LSSA

line of sight separation angle

MA

mayfly algorithm

PSO

particle swarm optimisation

SSE

stochastic subgradient error

UAV

unmanned aerial vehicle

WSN

wireless sensor network

Variables

$\boldsymbol{{G}}$

graph of network

$\boldsymbol{{N}}$

node set for bearing-only sensors

$\boldsymbol{{R}}$

edge set for bearing-only sensors

$N$

the number of UAVs in the formation

${d_i}(t)$

degree of node $i$ at time $t$

${S_i}$

the $i - {\rm{th}}$ UAV equipped with a pure orientation sensor

${\theta _i}$

elevation angle measured by ${S_i}$

${\varphi _i}$

azimuth angle measured by ${S_i}$

${\gamma _{ij}}$

LOS angle between the two bearing-only sensor nodes and the target

${\odot_{i \to j}}$

locating the spatial circle of the optimal neighbouring node ${S_j}$ based on node ${S_i}$

${\rm{\Delta }}H$

safe flight altitude for UAV formation

$l$

length of baseline

$Q$

number of formation flight levels based on constraints

${x_i},{y_i},{z_i}$

coordinate position of the $i - {\rm{th}}$ UAV

${x_t},{y_t},{z_t}$

coordinate position of the target

${J_i}$

indicator function to be optimised

${\psi _i}$

construction angle to be optimised

${\sigma _i}$

measurement errors inherent in UAVs

2.0 System model

In the collaborative detection model of UAV formation, each UAV is equipped with a bearing-only sensor, which can independently measure the elevation angle and azimuth angle towards the target. For the sake of simple and direct description in this paper, the entire UAV formation is defined as a bearing-only sensor network model with the same topology, in which each UAV is defined as a node in the network.

2.1 Bearing-only UAVs cooperative positioning model

In the AOA method for target positioning, the measurement unit of the system consists of a pair of bearing-only sensors. The intersection point of the two lines of sight from the UAV to the target is the spatial position of the target.

The communication topology of the UAV formation is a time-varying undirected graph $\boldsymbol{{G}}\left( {\boldsymbol{{N}},\boldsymbol{{R}}(t)} \right)$ . $N = \left\{ {1,2, \cdot \cdot \cdot, N} \right\}$ is the sensor node set and $N$ is the number of nodes in the network. $\boldsymbol{{R}}(t)$ is a directed edge set, and edge $\left( {i,j} \right) \in \boldsymbol{{R}}(t)$ means that node $j$ can exchange information with node $i$ as its neighbour. ${d_i}(t)$ represents the degree of node $i$ at time $t$ . For each node $i \in \boldsymbol{{N}}$ , ${\boldsymbol{{N}}_i} = \left\{ {j\left| {\left( {j,i} \right) \in \boldsymbol{{R}}} \right.} \right\}$ , otherwise, ${\boldsymbol{{N}}_i}\backslash \left\{ j \right\}$ .

The geometric situation of the UAV formation and the target is shown in Fig. 1. ${\left( {{x_i},{y_i},{z_i}} \right)^T}$ and ${\left( {{x_j},{y_j},{z_j}} \right)^T}$ are the coordinates of ${S_i}$ and ${S_j}$ , respectively. The elevation angle and azimuth angle measured by ${S_i}$ are ${\theta _i}$ and ${\varphi _i}$ , respectively. The elevation angle and azimuth angle measured by ${S_j}$ are ${\theta _j}$ and ${\varphi _j}$ , respectively.

Figure 1. Bearing-only detection model.

Through the AOA positioning method, the coordinates of the target can be obtained as

(1) \begin{align}\left( {\begin{array}{*{20}{c}}{{x_t} = {x_i} + \dfrac{{\left( {{y_j} - {y_i}} \right){\rm{cos}}\,{\varphi _j}\,{\rm{cos}}\,{\varphi _i} - \left( {{x_j} - {x_i}} \right){\rm{sin}}\,{\varphi _j}\,{\rm{cos}}\,{\varphi _i}}}{{{\rm{sin}}\left( {{\varphi _i} - {\varphi _j}} \right)}}}\\[9pt]{{y_t} = {y_i} + \dfrac{{\left( {{y_j} - {y_i}} \right){\rm{cos}}\,{\varphi _j}\,{\rm{sin}}\,{\varphi _i} - \left( {{x_j} - {x_i}} \right){\rm{sin}}\,{\varphi _j}\,{\rm{sin}}\,{\varphi _i}}}{{{\rm{sin}}\left( {{\varphi _i} - {\varphi _j}} \right)}}}\\[9pt]{{z_t} = {z_i} + \dfrac{{\left( {{y_j} - {y_i}} \right){\rm{cos}}\,{\varphi _j}\,{\rm{sin}}\,{\theta _i} - \left( {{x_j} - {x_i}} \right){\rm{sin}}\,{\varphi _j}\,{\rm{sin}}\,{\theta _i}}}{{{\rm{cos}}\,{\theta _i}{\rm{sin}}\left( {{\varphi _i} - {\varphi _j}} \right)}}}\end{array}} \right.\end{align}

When ${\varphi _i} = {\varphi _j}$ and ${\theta _i} = {\theta _j}$ , the target line of sight (LOS) of the sensor node pair coincides, the system of equations has no solution, and the triangulation positioning method fails.

In the observation model of the bearing-only sensor network depicted in Fig. 1, two sensor nodes form a triangular geometric relationship with the target for state estimation. ${R_i}$ represents the distance from ${S_i}$ to $T$ . ${S_j}$ represents the distance from ${S_j}$ to $T$ . ${R_{ij}}$ represents the distance from ${S_i}$ to ${S_j}$ . ${\gamma _{ij}}$ represents the LOS angle between the two bearing-only sensor nodes and the target, which is called the line of sight separation angle (LSSA). ${\gamma _i}$ is the angle between ${R_i}$ and ${R_{ij}}$ , which ${\gamma _j}$ can be obtained in the same way.

In the triangle formed by the two bearing-only sensors and the target, there is

(2) \begin{align}\left( {\begin{array}{*{20}{l}}{{R_i}} {}= {\dfrac{{{\rm{sin}}\,{\gamma _j}}}{{{\rm{sin}}\,{\gamma _{ij}}}}{R_{ij}}}\\[12pt]{{R_j}} {}= {\dfrac{{{\rm{sin}}\,{\gamma _i}}}{{{\rm{sin}}\,{\gamma _{ij}}}}{R_{ij}}}\\[12pt]{{\gamma _{ij}}} {}{ = \pi - \left( {{\gamma _i} + {\gamma _j}} \right)}\end{array}} \right.\end{align}

In order to analyse the observable preconditions, finding the partial derivative of ${R_i}$ ,

(3) \begin{align}{\delta {R_i}} & = \dfrac{{{\rm{sin}}\,{\gamma _j}}}{{{\rm{sin}}\,{\gamma _{ij}}}}\delta {R_{ij}} + \dfrac{{{R_{ij}}{\rm{cos}}\,{\gamma _j}}}{{{\rm{sin}}\,{\gamma _{ij}}}}\delta {\gamma _j} - \dfrac{{{R_{ij}}{\rm{cos}}\,{\gamma _{ij}}{\rm{sin}}\,{\gamma _j}}}{{{\rm{si}}{{\rm{n}}^2}{\gamma _{ij}}}}\delta {\gamma _{ij}}\nonumber\\ & = \dfrac{1}{{{\rm{sin}}\,{\gamma _{ij}}}}\left\{ {{\rm{sin}}\,{\gamma _j}\delta {R_{ij}} - {R_{ij}}{\rm{sin}}\,{\gamma _j}{\rm{cot}}\left( {{\gamma _i} + {\gamma _j}} \right)\delta {\gamma _i} + {R_{ij}}\left[ {{\rm{cos}}\,{\gamma _j} - {\rm{sin}}\,{\gamma _j}{\rm{cot}}\left( {{\gamma _i} + {\gamma _j}} \right)} \right]\delta {\gamma _j}} \right\}\end{align}

Analysing the partial derivative results, it can be seen that the estimation error of ${S_i}$ is related to $\delta {R_{ij}}$ , $\delta {\gamma _i},$ and $\delta {\gamma _{ij}}$ . When the LSSA is ${\gamma _{ij}} = \pi /2$ , the estimation error is minimal. When the LSSA reaches 0 or $\pi $ , the estimation error reaches infinity. At this time, the bearing-only sensor nodes are on a straight line, the observation model degenerates into a single-sensor detection problem, the target observability is reduced, and the distance estimation cannot be achieved. This is the same conclusion as drawn in Equation (1).

Therefore, in order to ensure that bearing-only sensors can successfully estimate the target state information, the triangular geometric relationship between pairs of sensor nodes and the target must avoid collinearity. The number of nodes in the bearing-only sensor network must satisfy $N,\left( {N \geqslant 2} \right)$ . When the baseline ${R_{ij}}$ between adjacent sensors in the detection array remains consistent, the LSSA approaches ${\gamma _{ij}} = \pi /2$ , and the distance estimation performance improves. Therefore, during the state estimation process in the bearing-only sensor network, the LSSA of node pairs needs to be considered. If the LSSA is too small, it will affect the observability of the system.

2.2 Constraints on the layered mechanism

Once the spatial position of the first bearing-only sensor ${S_i}$ is determined in 3D space, spatial geometric analysis reveals that for the second adjacent bearing-only sensor ${S_j}$ , there exist numerous optimal positions distributed along the intersection circle ${\odot_{i \to j}}$ of two spheres, as shown in Fig. 2. Considering the constrictions on the safe flight height of the UAV formation, the circle ${\odot_{i \to j}}$ is layered divided into $Q + 1$ layers. The bearing-only sensors of each layer are on the same horizontal plane, with the same height in the inertial coordinate system, and the elevation angles of the detected target are also equal. This is shown in the magnified portion of Fig. 2.

Figure 2. The detection model of bearing-only sensor network.

When the baseline ${R_{ij}}$ of the UAV formation is known and fixed, and the safe flight altitude ${\rm{\Delta }}H$ of the UAV formation is known and fixed, then there are $2Q$ spatial positions for the optimal detection formation planning of adjacent bearing-only sensors.

(4) \begin{align}Q = \frac{{z{{\rm{'}}_j} - {z_j}}}{{{\rm{\Delta }}H}}\end{align}

In the engineering application of mid-to-long-range target tracking and localisation, the practical flight control constraint ${\rm{\Delta H}}\lt\lt {{\rm{R}}_i}$ indicates that the safe flying altitude is much lower than the distance between the sensor and the target. The system model described in this paper involves the deployment of bearing-only sensor nodes positioned at various heights. These sensor nodes exhibit minimal variations in elevation angles when detecting the same target [Reference Han, Yang and Byun31]. The accuracy of detection in different arrays mainly stems from differences in azimuth angles. When primarily considering differences in azimuth angles, adjacent bearing-only sensor nodes are regarded as having the same elevation angle. If the nodes are on the same level, then within a set of $2Q$ spatial positions, there exist two optimal solutions $\left\{ {{S_j},S{{\rm{'}}_j}} \right\}$ that are on the same level as the initial node.

Traditional optimal array planning differs from the scenario studied in this paper. If the constraints of baseline and safety altitude for UAV formation are not considered, then the optimal array planning in such cases cannot be inferred from the conclusions of this paper.

3.0 Optimal neighbour node

Based on the constructed layered detection model of UAV formation, the optimal geometric relationship between adjacent UAVs is firstly obtained through spatial geometric analysis. Then, by employing coordinate transformation, the analytical solution for the optimal positions of adjacent UAVs is derived.

3.1 Formation of optimal observation unit

The position $\left( {{x_i},{y_i},{z_i}} \right)$ of the bearing-only sensor node ${S_i}$ in the inertial coordinate system $OXYZ$ is known. At the same time, the neighbour node ${S_j}$ that forms the optimal formation with the target is on the plane of the LOS coordinate system ${X_1}O{Y_1}$ . Therefore, it can be simplified to calculate the position coordinates of ${S_j}$ .

Taking node ${S_i}$ as the coordinate origin, based on the sight coordinate system ${O_1}{X_1}{Y_1}{Z_1}$ , a new coordinate system ${O_2}{X_2}{Y_2}{Z_2}$ is obtained by converting the sight elevation angle ${q_\theta }$ and the sight azimuth angle ${q_\varphi }$ , which is defined as the observation plane coordinate system, as shown in Fig. 3. The X-axis direction is the projection of the line connecting the sensor node ${S_i}$ and the target $T$ on the horizontal plane. The Y-axis is located on the vertical plane of the X-axis in the observation plane coordinate system, and upward is positive. Both sight angles are independently measured by bearing-only sensor node ${S_i}$ .

Figure 3. Conversion of coordinate systems.

In the observation plane coordinate system, the coordinate system of neighbour node ${S_j}$ is $\left( {{x_j},{y_j}} \right)$ , and the target position coordinate is $\left( {{x_T},{y_T}} \right)$ . If and only if ${S_i}T = {S_j}T$ , the pair of observation nodes composed of ${S_i}$ and ${S_j}$ is optimal, as shown in Fig. 4. The proposal and verification of this condition are detailed in Section 4.

Figure 4. Coordinate system conversion.

The triangular relationship between the node pair and the target is expressed as

(5) \begin{align}\sqrt {{{\left( {{x_T} - 0} \right)}^2} + {{\left( {{k_i}{x_T} - 0} \right)}^2}} = \sqrt {{{\left( {{x_T} - {x_j}} \right)}^2} + {{\left( {{k_i}{x_T} - {y_j}} \right)}^2}} \end{align}

where, ${k_i} = {\rm{tan}}{\beta _i}$ . ${\beta _i}$ is the azimuth angle of the target measured by ${S_i}$ , ${\beta _i} = {q_\varphi }$ .

The baseline length from node ${S_i}$ to node ${S_j}$ is $l$ .

(6) \begin{align}\sqrt {{{\left( {{x_i} - {x_j}} \right)}^2} + {{\left( {{y_i} - {y_j}} \right)}^2}} = \sqrt {{{\left( {0 - {x_j}} \right)}^2} + {{\left( {0 - {y_j}} \right)}^2}} = l\end{align}

The LSSA is recorded as ${\gamma _{ij}}$ ,

(7) \begin{align}{\rm{cos}}\,{\gamma _{ij}} = \frac{{4\left[ {{{\left( {{x_T} - 0} \right)}^2} + {{\left( {{k_i}{x_T} - 0} \right)}^2}} \right] - {l^2}}}{{2\left[ {{{\left( {{x_T} - 0} \right)}^2} + {{\left( {{k_i}{x_T} - 0} \right)}^2}} \right]}} = 1 - \frac{{{l^2}}}{{2\left[ {{{\left( {{x_T} - 0} \right)}^2} + {{\left( {{k_i}{x_T} - 0} \right)}^2}} \right]}}\end{align}

The slope of the straight line passing through ${S_j}$ is recorded as ${k_j}$ , and the following relationship can be obtained.

(8) \begin{align}{k_j} = {\rm{tan}}{\beta _j} = {\rm{tan}}\left( {{\beta _i} + {\gamma _{ij}}} \right)\end{align}

Assume that the equation of the straight line passing through the node ${S_j}$ and the $T$ is $y = {k_j}x + b$ , and the target point $T\left( {{x_T},{k_i}{x_T}} \right)$ of the straight line is solved. The intercept of the straight line is $b = {k_i}{x_T} - {k_j}{x_T}$ . The straight line is expressed as:

(9) \begin{align}y = {k_j}x + {k_i}{x_T} - {k_j}{x_T}\end{align}

Substituting the coordinates of node ${S_j}$ into Equation (9).

(10) \begin{align}{y_j} = {k_j}{x_j} + {k_i}{x_T} - {k_j}{x_T} = {\rm{tan}}{\beta _j} \cdot {x_j} + {\rm{tan}}{\beta _i} \cdot {x_T} - {\rm{tan}}{\beta _j} \cdot {x_T}\end{align}

The simultaneous equations of the unknown parameters ${x_T}$ , ${x_j}$ , ${y_j}$ , ${\gamma _{{\rm{i}}j}}$ , ${k_j}$ are:

(11) \begin{align}\left\{ {\begin{array}{*{20}{l}}{\sqrt {{{\left( {{x_T} - 0} \right)}^2} + {{\left( {{k_i}{x_T} - 0} \right)}^2}} = \sqrt {{{\left( {{x_T} - {x_j}} \right)}^2} + {{\left( {{k_i}{x_T} - {y_j}} \right)}^2}} }\\[9pt]{\sqrt {{{\left( {{x_i} - {x_j}} \right)}^2} + {{\left( {{y_i} - {y_j}} \right)}^2}} = \sqrt {{{\left( {0 - {x_j}} \right)}^2} + {{\left( {0 - {y_j}} \right)}^2}} = l}\\[9pt]{{\rm{cos}}\,{\gamma _{ij}} = = 1 - \dfrac{{{l^2}}}{{2\left[ {{{\left( {{x_T} - 0} \right)}^2} + {{\left( {{k_i}{x_T} - 0} \right)}^2}} \right]}}}\\[12pt]{{\rm{tan}}{\beta _j} = {\rm{tan}}\left( {{\beta _i} + {\gamma _{ij}}} \right)}\\[5pt]{{y_j} = {\rm{tan}}{\beta _j}{ \cdot _j} + {\rm{tan}}{\beta _i} \cdot {x_T} - {\rm{tan}}{\beta _j} \cdot {x_T}}\end{array}} \right.\end{align}

In this way, the coordinates of two neighbouring nodes ${({x_j},{y_j},{z_j})_{{O_2}{X_2}{Y_2}{Z_2}}}$ and ${({x_j},y{{\rm{'}}_j},{z_j})_{{O_2}{X_2}{Y_2}{Z_2}}}$ can be determined, and any neighbouring node can be paired with the initial sensor node ${({x_i},{y_i},{z_i})_{{O_2}{X_2}{Y_2}{Z_2}}}$ to form the optimal sensor observation unit.

3.2 Actual coordinates of neighbouring UAVs

Through coordinate transformation, the coordinate ${({x_j},{y_j})_1}$ on the ${x_1}{o_1}{y_1}$ plane under the LOS coordinate can be obtained.

(12) \begin{align}{\left[ {\begin{array}{*{20}{l}}{{x_j}}\\[3pt]{{y_j}}\end{array}} \right]_1} = \left[ {\begin{array}{l@{\quad}l}{{\rm{cos}}\,{q_\varphi }} {}& {{\rm{sin}}\,{q_\varphi }}\\[3pt]{ - {\rm{sin}}\,{q_\varphi }} {}&{{\rm{cos}}\,{q_\varphi }}\end{array}} \right]{\left[ {\begin{array}{*{20}{l}}{{x_j}}\\[3pt]{{y_j}}\end{array}} \right]_2}\end{align}

The coordinates under the LOS coordinate system are:

(13) \begin{align}{\left[ {\begin{array}{c}{{x_j}}\\[3pt]{{y_j}}\\[3pt]{{z_j}}\end{array}} \right]_1}\left[ {\begin{array}{c}{{x_j} \cdot {\rm{cos}}\,{q_\varphi } + {y_j} \cdot {\rm{sin}}\,{q_\varphi }}\\[3pt]{ - {x_j} \cdot {\rm{sin}}\,{q_\varphi } + {y_j} \cdot {\rm{cos}}\,{q_\varphi }}\\[3pt]0\end{array}} \right]\end{align}

The transformation matrix from the LOS coordinate system to the ground coordinate is:

(14) \begin{align}C = \left[ {\begin{array}{c@{\quad}c@{\quad}c}{{\rm{cos}}\,{q_\varphi }{\rm{cos}}\,{q_\theta }} {}& { - {\rm{cos}}\,{q_\theta }{\rm{sin}}\,{q_\varphi }} {} {}& {{\rm{sin}}\,{q_\theta }}\\[3pt]{{\rm{sin}}\,{q_\varphi }} {}& {}{{\rm{cos}}\,{q_\varphi }} {}& {}0\\[3pt]{ - {\rm{sin}}\,{q_\theta }{\rm{cos}}\,{q_\varphi }} {} {}& {{\rm{sin}}\,{q_\theta }{\rm{sin}}\,{q_\varphi }} {}& {}{{\rm{cos}}\,{q_\theta }}\\[3pt]\end{array}} \right]\end{align}

The position coordinates of the neighbour node ${S_j}$ in inertial coordinates can be obtained:

(15) \begin{align}\begin{array}{l}\left[ {\begin{array}{c@{\quad}c@{\quad}c}{{\rm{cos}}\,{q_\varphi }{\rm{cos}}\,{q_\theta }}& {}{ - {\rm{cos}}\,{q_\theta }{\rm{sin}}\,{q_\varphi }} {}&{{\rm{sin}}\,{q_\theta }}\\[3pt]{{\rm{sin}}\,{q_\varphi }} {}&{{\rm{cos}}\,{q_\varphi }}& {}0\\[3pt]{ - {\rm{sin}}\,{q_\theta }{\rm{cos}}\,{q_\varphi }}& {}{{\rm{sin}}\,{q_\theta }{\rm{sin}}\,{q_\varphi }}& {}{{\rm{cos}}\,{q_\theta }}\\\end{array}} \right]\left[ {\begin{array}{c}{{x_j} \cdot {\rm{cos}}\,{q_\varphi } + {y_j} \cdot {\rm{sin}}\,{q_\varphi }}\\[3pt]{ - {x_j} \cdot {\rm{sin}}\,{q_\varphi } + {y_j} \cdot {\rm{cos}}\,{q_\varphi }}\\[3pt]0\end{array}} \right]\\[26pt] = \left[ {\begin{array}{c} {\left( {{\rm{cos}}\,{q_\varphi }{\rm{cos}}\,{q_\theta }} \right)\left[ {{x_j} \cdot {\rm{cos}}\,{q_\varphi } + {y_j} \cdot {\rm{sin}}\,{q_\varphi }} \right] + \left( { - {\rm{cos}}\,{q_\theta }{\rm{sin}}\,{q_\varphi }} \right)\left[ { - {x_j} \cdot {\rm{sin}}\,{q_\varphi } + {y_j} \cdot {\rm{cos}}\,{q_\varphi }} \right]}\\[7pt]{{\rm{sin}}\,{q_\varphi }\left[ {{x_j} \cdot {\rm{cos}}\,{q_\varphi } + {y_j} \cdot {\rm{sin}}\,{q_\varphi }} \right] + {\rm{cos}}\,{q_\varphi }\left[ { - {x_j} \cdot {\rm{sin}}\,{q_\varphi } + {y_j} \cdot {\rm{cos}}\,{q_\varphi }} \right]}\\[7pt]{\left( { - {\rm{sin}}\,{q_\theta }{\rm{cos}}\,{q_\varphi }} \right)\left[ {{x_j} \cdot {\rm{cos}}\,{q_\varphi } + {y_j} \cdot {\rm{sin}}\,{q_\varphi }} \right] + \left( {{\rm{sin}}\,{q_\theta }{\rm{sin}}\,{q_\varphi }} \right)\left[ { - {x_j} \cdot {\rm{sin}}\,{q_\varphi } + {y_j} \cdot {\rm{cos}}\,{q_\varphi }} \right]}\end{array}} \right]\end{array}\end{align}

The coordinate value of the Z-axis is given by the elevation angles of the layer, ${z_i} = {z_j}$ . By analogy, in the detection array planning of the UAV formation, when the optimal position of the UAVn $\left( {n \geqslant 2} \right)$ is determined, a pending position set ${\mathbb{N}^{{2^n} - 1}}$ with ${2^n} - 1$ elements will be generated. The undetermined position set ${\mathbb{N}^{{2^n} - 1}}$ includes the feasible solution set ${\mathbb{R}^{n + 1}}$ and its complement. The feasible solution set includes the optimal position solution set ${\mathbb{R}^{\rm{*}}}$ and the suboptimal position solution set ${\mathbb{R}^{\backslash {\rm{*}}}}$ , but both of them satisfy the optimal adjacent UAV conditions. How to determine the ${\mathbb{R}^{\rm{*}}}$ that can form the optimal detection formation is a complex combinatorial optimisation problem. Especially when there are a large number of UAV formations, a stable and reliable optimisation algorithm that meets real-time requirements is needed to solve this problem.

4.0 Objective function to be optimise

When the node size of the system network is small, the exhaustive combinatorial optimisation method can be used to solve the problem. However, when the system network is used in application scenarios with large-scale node numbers and high real-time requirements, it is necessary to design a distributed optimisation algorithm suitable for this model. In this way, the optimal geometric configuration of the system network for detecting target under layered constraints can be calculated with better real-time performance and interference resistance. The prerequisite for implementing distributed algorithms to optimise the array of UAV formation is to design a reasonable and efficient objective function.

In the process of target localisation using AOA with bearing-only sensors, different geometric configurations of the sensors can lead to variations in the shape and size of intersecting areas, as shown in Fig. 5.

Figure 5. GDOP for node pair positioning.

The accuracy of the system localisation depends not only on the position error and angle measurement error of the sensor nodes, but also on the geometric configuration between the two UAVs and the target. This type of target localisation error, originating from the position error of the UAVs and the angle measurement error of the sensors, and propagated through the relative geometric relationships between the UAVs and the target, is referred to as GDOP.

$\boldsymbol{{X}} \in {\boldsymbol{{R}}^n}$ is the position vector, $\boldsymbol{{Z}} \in {\boldsymbol{{R}}^m}$ is the measurement vector, and the relationship between the measurement vector $\boldsymbol{{Z}}$ and the position vector $\boldsymbol{{X}}$ can be expressed as:

(16) \begin{align}\boldsymbol{{Z}} = h(\boldsymbol{{X}})\end{align}

By deriving the Equation (16) with respect to the position vector,

(17) \begin{align}d\boldsymbol{{Z}} = \frac{{\partial h}}{{\partial \boldsymbol{{X}}}}d\boldsymbol{{X}} = \boldsymbol{{H}}d\boldsymbol{{X}}\end{align}

where, $\boldsymbol{{H}} \in {R^{m \times n}}$ is the Jacobian matrix,

(18) \begin{align}E\left( {d\boldsymbol{{X}}d{\boldsymbol{{X}}^{\rm{T}}}} \right) = E\left[ {{\boldsymbol{{H}}^{ - 1}}d\boldsymbol{{Z}}d{\boldsymbol{{Z}}^{\rm{T}}}{{\left( {{\boldsymbol{{H}}^{ - 1}}} \right)}^{\rm{T}}}} \right] = {\boldsymbol{{H}}^{ - 1}}E\left( {d\boldsymbol{{Z}}d{\boldsymbol{{Z}}^{\rm{T}}}} \right){\left( {{\boldsymbol{{H}}^{ - 1}}} \right)^{\rm{T}}} = {\boldsymbol{{H}}^{ - 1}}\boldsymbol{{R}}{\left( {{\boldsymbol{{H}}^{ - 1}}} \right)^{\rm{T}}}\end{align}

It can be seen that GDOP is the square root of the sum of the diagonal elements of $E\left\{ {d\boldsymbol{{X}}d{\boldsymbol{{X}}^{\rm{T}}}} \right\}$ .

(19) \begin{align}{\rm{GDOP}} = \sqrt {trace[{\boldsymbol{{H}}^{ - 1}}\boldsymbol{{R}}\left( {{\boldsymbol{{H}}^{ - 1}}{)^{\rm{T}}}} \right]} \end{align}

where, $\boldsymbol{{R}}$ is the measurement error covariance matrix.

GDOP reflects the impact of the geometric configuration of the UAV formation and the target on the localisation accuracy. A larger GDOP value indicates greater localisation errors. However, in 3D space, the GDOP objective function is too complex to be directly used as a metric for optimising the UAV formation array.

In the research scenario of this paper, the two UAVs composing the optimal observation unit are located on the same plane. Based on the geometric analysis in Section 3, the elevation angles measured by the optimal observation unit remains consistent. Furthermore, the height difference between the UAV formation and the target is much smaller than the distance between the UAV formation and the target. The influence of the elevation angle on target localisation accuracy is minimal, the triangle formed by the observation unit and the target can be reasonably abstracted as being on the same plane. Thus, only the optimisation problem of the azimuth angle on the observation plane needs to be considered. Building upon the analysis mentioned above, this paper proposes a GDOP objective function within the observation plane, and under certain conditions, the nonlinear objective function can be simplified to a quadratic form.

In the observation plane ${O_2}{X_2}{Y_2}{Z_2}$ , two UAVs and the target form a certain plane triangle geometric relationship, as shown in Fig. 6.

Figure 6. Geometric relations on the observation plane.

The coordinate of the target is $\left( {{x_t},{y_t}} \right)$ , the coordinates of the node pair are $\left( {{x_i},{y_i}} \right)$ and $\left( {{x_j},{y_j}} \right)$ . The measured azimuth angles are ${\beta _i}$ and ${\beta _j}$ , respectively. Assume that the angle measurement errors of the two sensors are independent of each other, the variances are $\sigma _i^2$ and $\sigma _j^2$ , respectively. $\boldsymbol{{R}} = diag\left( {\sigma _i^2,\sigma _j^2} \right)$ is the covariance matrix of the measurement errors.

(20) \begin{align}\boldsymbol{{H}} = \left[ {\begin{array}{*{20}{l}}{ - \dfrac{{{\rm{sin}}\,{\varphi _i}}}{{{R_i}}}} {}{\dfrac{{{\rm{cos}}\,{\varphi _i}}}{{{R_i}}}}\\[12pt]{ - \dfrac{{{\rm{sin}}\,{\varphi _j}}}{{{R_j}}}} {}{\dfrac{{{\rm{cos}}\,{\varphi _j}}}{{{R_j}}}}\end{array}} \right]\end{align}

Then substitute the expressions of $\boldsymbol{{R}}$ and $\boldsymbol{{H}}$ into Equation (19),

(21) \begin{align}{\rm{GDOP}} = \sqrt {\frac{{R_i^2\sigma _i^2 + R_j^2\sigma _j^2}}{{{\rm{si}}{{\rm{n}}^2}\left( {{\varphi _i} - {\varphi _j}} \right)}}} \end{align}

where, ${R_i} = \sqrt {{{({x_t} - {x_i})}^2} + {{({y_t} - {y_i})}^2}} $ . Because ${R_i}$ can be obtained indirectly using the AOA triangulation method, Equation (21) can be further expressed as:

(22) \begin{align}{\rm{GDOP}} = \sqrt {\frac{{R_{ij}^2\left( {\sigma _i^2{\rm{si}}{{\rm{n}}^2}{\varphi _i} + \sigma _j^2{\rm{si}}{{\rm{n}}^2}{\varphi _j}} \right)}}{{{\rm{si}}{{\rm{n}}^4}\left( {{\varphi _i} - {\varphi _j}} \right)}}} \end{align}

where, ${R_{ij}} = \sqrt {{{({x_2} - {x_1})}^2} + {{({y_2} - {y_1})}^2}} $ is the distance between pairs of nodes.

In the array planning process of UAV formation, baseline ${R_{ij}} = l$ is treated as a constant in each optimisation step, and the angular measurement errors of bearing-only sensors are also considered as a constant ${\sigma _i} = \sigma \left( {i = 1,2, \cdot \cdot \cdot, n} \right)$ . The initial objective function is

(23) \begin{align}{J^{\rm{*}}} = \frac{{\sqrt {{\rm{si}}{{\rm{n}}^2}{\varphi _i}(t) + {\rm{si}}{{\rm{n}}^2}{\varphi _j}(t)} }}{{{\rm{si}}{{\rm{n}}^2}\left[ {{\varphi _i}(t) - {\varphi _j}(t)} \right]}}\end{align}

When the distance between the observation unit and the target is much greater than the distance of the baseline, as shown in Fig. 6, let the convergence angle of the LOS of the bearing-only sensor node pair be $2\delta $ . According to the cosine theorem, $2\delta $ is obtained as a small quantity. Drawing a perpendicular line from the convergence point to the baseline of the node pair, let the angle between the perpendicular line and the bisector of the convergence angle be ${\rm{\Delta }}\delta $ . To simplify the problem, only consider the case where ${\rm{\Delta }}\delta $ is a very small quantity.

(24) \begin{align}\left\{ {\begin{array}{*{20}{l}}{} {}{{\varphi _i} = \pi /2 - \left( {\delta + {\rm{\Delta }}\delta } \right)}\\[4pt]{} {}{{\varphi _j} = \pi /2 + \left( {\delta - {\rm{\Delta }}\delta } \right)}\end{array}} \right.\end{align}

Substituting Equation (24) into the function numerator of Equation (23),

(25) \begin{align}& {\sqrt {{\rm{si}}{{\rm{n}}^2}{\varphi _i} + {\rm{si}}{{\rm{n}}^2}{\varphi _j}} }\nonumber\\& = \sqrt {{\rm{si}}{{\rm{n}}^2}\left[ {\pi /2 - \left( {\delta + {\rm{\Delta }}\delta } \right)} \right] + {\rm{si}}{{\rm{n}}^2}\left[ {\pi /2 + \left( {\delta - {\rm{\Delta }}\delta } \right)} \right]} = \sqrt {{\rm{co}}{{\rm{s}}^2}\left( {\delta + {\rm{\Delta }}\delta } \right) + {\rm{co}}{{\rm{s}}^2}\left( {\delta - {\rm{\Delta }}\delta } \right)} \end{align}

Since both $\delta $ and ${\rm{\Delta }}\delta $ are all very small quantities, their algebraic sum and algebraic product are still very small quantities, then there are ${\rm{co}}{{\rm{s}}^2}\left( {\delta + {\rm{\Delta }}\delta } \right) \approx 1$ , ${\rm{co}}{{\rm{s}}^2}\left( {\delta - {\rm{\Delta }}\delta } \right) \approx 1$ . Because ${\varphi _j} - {\varphi _i} = 2\delta $ , ${\rm{sin}}\left( {{\varphi _i} - {\varphi _j}} \right) \approx \left( {{\varphi _i} - {\varphi _j}} \right)$ . Therefore, Equation (24) can be simplified to

(26) \begin{align}{J^{\rm{*}}} \approx \frac{{\sqrt 2 }}{{{{[{\varphi _i}(t) - {\varphi _j}(t)]}^2}}}\end{align}

The optimal detection array can be obtained by minimising the value ${J^{\rm{*}}}$ in Equation (24), and the proposed method requires that the objective function to be optimised is a convex function. In order to obtain the standard form of a quadratic objective function to simplify the solution process, the problem is solved by maximising the reciprocal of the original objective function.

(27) \begin{align}J{{\rm{'}}^{\rm{*}}} = \frac{{\sqrt 2 {{[{\varphi _i}(t) - {\varphi _j}(t)]}^2}}}{2}\end{align}

As shown in Fig. 6, on the observation plane composed of the optimal observation node pair and target. When ${\varphi _i}(t) = {\varphi _j}(t)$ , the minimum objective function value of Equation (24) can be obtained. At this time, the optimised angle ${\psi _i}$ formed by the angular bisector of the LSSA ${\gamma _{ij}}$ and the baseline ${R_{ij}}$ approaches infinitely to $\pi /2$ , then

(28) \begin{align}{J_i} = \frac{{\sqrt 2 }}{2}{\left[ {{\psi _i}(t) - \frac{\pi }{2}} \right]^2}\end{align}

After the optimised angles ${\psi _i}$ are determined, considering that the optimal formation in the observation plane is a rigid isosceles triangle, the azimuth angle ${\varphi _i}$ of ${S_i}$ can be determined. The elevation angle ${\theta _i}$ of ${S_i}$ is constrained and fixed by the layered mechanism, and finally through the coordinate conversion in Section 3.2, the actual position planning of each UAV in the network is completed.

5.0 The proposed distributed algorithm

In the optimisation problem of locating a single target by UAV formation, considering the constraints of flight safety altitude and fixed baseline, the layered mechanism and the GDOP objective function of the observation plane are introduced. By reasonably simplifying the detection array planning problem for bearing-only sensor network, the complexity and implementation difficulty of the algorithm are reduced, ensuring that it can bring practical advantages in specific problem contexts.

In the communication network of UAV formation, information changes dynamically over time and gradually increases, but the available resources are limited. Therefore, it is extremely important and necessary to design a distributed, real-time optimisation algorithm to solve the detection array planning problem. Distributed algorithm can decompose a problem into multiple subproblems, which are processed in parallel by multiple nodes. This can effectively utilise distributed computing resources, accelerate the entire optimisation process, and ensure the real-time requirements of the system. If one or more system nodes fail, other nodes can still continue to work, ensuring system scalability, stability and efficiency [Reference Ge, Wang and Dong32]. More importantly, it can better adapt to large-scale, distributed and real-time application scenarios with high requirements.

The DSSPA is a distributed algorithm used to solve convex optimisation problems, which combines the idea of stochastic gradient descent (SGD) and the projection method. Since DSSPA uses subgradients instead of gradients, it is more suitable for problems where the GDOP objective function contains the parts of nonsmooth and convex optimisation. Additionally, by applying a projection operation to each feasible solution, it can ensure that the constraints of the problem are satisfied when updating parameters. This is crucial for constraints of fixed baseline and safe flight altitude in planning problems. In the process of collaborative optimisation, each node calculates local subgradients and collaborates through communication to obtain global subgradient information, thereby achieving global optimisation progress.

A pair of UAVs that satisfy the conditions of the optimal observation form a unit. Each UAV has a local cost function. The UAV formation needs to seek a globally consistent optimal solution within the intersection of feasible sets. The sequence of convex cost functions is denoted as $\left\{ {{J_1},{J_2}, \cdot \cdot \cdot, {J_N}} \right\}$ .

Transforming the detection array planning problem of UAV formation into an optimisation model with constraints.

(29) \begin{align}\mathop {{\rm{max}}\ J\left( \psi \right)}\limits_{{\bf{x}} \in {\mathbb{N}^{{{(2Q)}^n} - 1}}} {\rm{\;}} = \mathop \sum \limits_{i = 1}^n {J_i}\left( {{\psi _i}} \right){\rm{\;\;\;\;\;\;\;\;s}}.{\rm{t}}.\psi \in {\mathbb{R}^{n + 1}}\end{align}

Where, ${J_i}$ is the cost function of node ${S_i}$ , and $\psi \in {\mathbb{R}^{n + 1}}$ is the constraint set. Each node only knows its own cost function, and it can share the results of its operations with neighbours.

A doubly stochastic matrix has all nonnegative real elements, with the sum of elements in each row and each column equal to 1. In a time-varying undirected network, the weight matrix represents the probability of connections between nodes, and these connections may change over time. Therefore, in a time-varying network, the weight matrix evolves with time and may not be a doubly stochastic matrix. To address this problem, this paper adopts the weight balancing method.

If for any node $i \in N$ , its weight ${\omega _i}$ satisfies:

(30) \begin{align}{\omega _i}{d_i}(t) = \mathop \sum \nolimits_{j \in {N_i}} {\omega _j}\end{align}

Then the weight ${\omega _i}$ is said to balance a time-varying directed graph $\boldsymbol{{G}}(t)$ . The total weight $\mathop \sum \nolimits_{j \in {N_i}} {\omega _j}$ of the neighbour node is equal to the total weight ${\omega _i}{d_i}(t)$ of the node ${S_i}$ .

Definition 1. Each local cost function ${J_i}\left( \psi \right)$ is convex, for all ${\psi _i},{\psi _j} \in {\mathbb{R}^{ + 1}}$ , the function ${J_i}$ satisfies

(31) \begin{align}{J_i}\left( {{\psi _i}} \right) - {J_i}\left( {{\psi _j}} \right) \geqslant \nabla {J_i}{({\psi _j})^{\rm{T}}}\left( {{\psi _i} - {\psi _j}} \right) + \frac{{{\sigma _i}}}{2}{\psi _i} - {\psi _j}^2\end{align}

where, ${\sigma _i} \geqslant 0$ is a constant and $\nabla {J_i}\left( {{\psi _j}} \right)$ is the subgradient of function ${J_i}\left( {{\psi _i}} \right)$ at ${\psi _i} = {\psi _j}$ . When ${\sigma _i} \gt 0$ , each local cost function is a strongly convex function, otherwise, ${J_i}$ is a general convex function.

Definition 2. For any node ${S_i}\left( {i = 1,2, \cdot \cdot \cdot, n} \right)$ , the subgradient $\nabla {J_i}\left( \psi \right)$ of function ${J_i}\left( \psi \right)$ is consistent and bounded at ${\mathbb{R}^{n + 1}}$ . For all $\psi \in {\mathbb{R}^{n + 1}}$ , there is $\nabla {J_i}\left( \psi \right) \le {L_{max}}$ , where, ${L_{max}} = \mathop {ma{x_{}}}\limits_i \;{L_i},i = 1,2, \cdot \cdot \cdot, n$ . In the DSSPA, let ${\psi _i}(t) \in {\mathbb{R}^{n + 1}}$ be the iteration value of node ${S_i}$ at time $t$ , which is the estimated value of the optimal solution. Therefore, the update rule of the estimated value of node $i \in \boldsymbol{{N}}$ at time $t$ is

(32) \begin{align}{z_i}(t) = \left( {1 - {w_i}(t){d_i}(t)} \right){\psi _i}(t) + \mathop \sum \limits_{j \in {N_i}(t)} {w_j}(t){\psi _j}(t)\end{align}

(33) \begin{align}{\psi _i}\left( {t + 1} \right) = {\Omega _\mathbb{R}}\left[ {{z_i}(t) + \alpha (t){P_i}(t)} \right]\end{align}

where, $\alpha (t)$ is the step sequence, ${P_i}(t)$ is the abbreviation of ${P_i}\left( {{\psi _i}(t)} \right)$ . ${P_i}\left( {{\psi _i}(t)} \right)$ is the noise subgradient of function ${J_i}\left( x \right)$ at $\psi = {\psi _i}(t)$ , and ${\Omega _\mathbb{R}}$ is the Euclidean projection operator.

It can be seen from Equations (32) and (33) that each node first linearly fuses the function values of itself and its neighbours, and then adjusts along the noise subgradient direction of its cost function. Finally, by projecting the adjustment value onto the constraint set, an updated objective function is obtained.

${\mathbb{N}^{{2^n} - 1}}$ is a set of projections with undetermined locations. Correspondingly, ${\hat \psi _i} \in {\mathbb{R}^{n + 1}}$ is the calculated position of the UAV at the optimal observation position, and ${P_C}\left( {{{\hat \psi }_i}} \right)$ is the projection of ${\hat \psi _i}$ on the constraint set ${\mathbb{R}^{\rm{*}}}$ . The optimal set is nonempty.

If exists ${\psi _i} \in {\mathbb{R}^{\rm{*}}}$ , then ${P_C}\left( {{{\hat \psi }_i}} \right) = {\psi _i}$ , if ${\psi _i} \notin {\mathbb{R}^{\rm{*}}}$ , then there is

(34) \begin{align}{P_C}\left( {{{\hat \psi }_i}} \right) = \mathop {{\rm{arg\,min}}}\limits_{{\psi _i} \in {\mathbb{R}^{\rm{*}}}} {\rm{\;}}{\psi _i} - {\hat \psi _i}\end{align}

${P_C}\left( {{{\hat \psi }_i}} \right)$ is an element in ${\mathbb{R}^{\rm{*}}}$ that satisfies the minimisation conditions in Equation (34). ${\hat \psi _i}$ does not satisfy the constraint set conditions is transformed into ${\psi _i}$ that satisfies the constraint set conditions through Euclidean norm projection, and ${\psi _i}$ at this time is the closest element to ${\hat \psi _i}$ .

The noise subgradient is

(35) \begin{align}{\boldsymbol{{P}}_i}\left( {{\psi _i}(t)} \right) = \nabla {J_i}\left( {{\psi _i}(t)} \right) + {\delta _i}(t)\end{align}

where, $\nabla {J_i}\left( {{\psi _i}(t)} \right)$ is the subgradient of ${f_i}\left( \psi \right)$ at $\psi = {\psi _i}(t)$ , and ${\delta _i}(t) \in {\mathbb{R}^{n + 1}}$ is the stochastic subgradient error (SSE). Let ${{\rm{\Gamma }}_t}$ represents all historical information generated by Equations (32) and (33). Therefore, the assumption for the SSE is as follows.

Assumption 2. For any node ${S_i}\left( {i = 1,2, \cdot \cdot \cdot, n} \right)$ , it is assumed that the SSE ${\delta _i}(t)$ is an independent variable with a mean value of 0, $E[{\delta _i}(t)|{\Gamma _{t - 1}}] = 0$ . In addition, it is assumed that the norm ${\delta _i}(t)$ of the SSE is consistent and bounded.

(36) \begin{align}E\left[ {{\delta _i}(t)|{{\rm{\Gamma }}_{t - 1}}} \right] \le {\nu _i}\end{align}

where, ${\nu _i}$ is a positive integer. In order to achieve weight balance, node $i \in N$ updates its own weight according to the following rules:

(37) \begin{align}{w_i}\left( {t + 1} \right) = \frac{1}{2}{w_i}(t) + \frac{1}{{{d_i}(t)}}\mathop \sum \limits_{j \in {N_i}(t)} \frac{1}{2}{w_j}(t)\end{align}

where, ${w_i}(t)$ is the weight of node ${S_i}$ at time $t$ . The optimal solution set of the optimisation problem is ${\mathbb{R}^*} \buildrel \Delta \over = \{ \psi \in {\mathbb{R}^{n + 1}}|J\left( \psi \right) \buildrel \Delta \over = \mathop {ma{x_{}}}\limits_{\psi \in {\mathbb{N}^{{2^n} - 1}}} \;J\left( \psi \right)\} $ , Therefore, this paper assumes that the optimal solution set ${\mathbb{R}^*}$ is nonempty.

Theorem 1. In accordance with the validity of Assumptions 1-4. The optimal solution set ${\mathbb{R}^*}$ is nonempty. For any ${S_i}\left( {i = 1,2, \cdot \cdot \cdot, n} \right)$ and positive step size $\alpha (t)$ , the estimated sequence ${x_i}(t)$ is generated by Equations (32) and (33). If the step size $\alpha (t)$ is positive and satisfies the following attenuation conditions. For $t \gt s \geqslant 1$ ,

(38) \begin{align}\mathop \sum \limits_{t = 1}^\infty \alpha (t) = \infty, \;\;\;\;\mathop \sum \limits_{t = 1}^\infty {\alpha ^2}(t)\lt \infty, \;\;\;\;\alpha (t)\alpha \left( s \right)\end{align}

Then each estimated sequence ${x_i}(t)$ converges to some optimal solution $x \in {\mathbb{R}^*}$ with probability 1. For all ${S_i}\left( {i = 1,2, \cdot \cdot \cdot, n} \right)$ , $\mathop {li{m_{}}}\limits_{t \to \infty } \;{x_i}(t) = {x^*}$ holds with probability 1.

It can be seen from Theorem 1 that in a time-varying topology graph, all iteration values can asymptotically converge to some optimal solution. Therefore, the constrained optimisation problem can be solved by the proposed DSSPA.

For the bearing-only collaborative detection array planning problem of UAV formation, the pseudocode of DSSPA is given in Fig. 7.

Figure 7. The pseudocode of DSSPA.

In the analysis of the convergence rate of Equations (32) and (33), the weighted average variable ${\{ {x_i}(t)\} _{t \geqslant 0}}$ of the estimated value sequence is introduced.

For any $t \geqslant 2$ ,

(39) \begin{align}{\hat \psi _i}(t) = \frac{{\mathop \sum \nolimits_{s = 1}^t \left( {s - 1} \right){\psi _i}\left( s \right)}}{{t\left( {t - 1} \right)/2}}\end{align}

Therefore, the following recursion relationship is

(40) \begin{align}{\hat \psi _i}\left( {t + 1} \right) = \frac{{{\psi _i}\left( {t + 1} \right) + S(t){{\hat \psi }_i}(t)}}{{S\left( {t + 1} \right)}}\end{align}

where, $S(t) = t\left( {t - 1} \right)/2$ , $t \geqslant 2$ , ${\hat \psi _i}\left( 1 \right) = {\psi _i}\left( 0 \right)$ . The local cost function ${J_i}\left( \psi \right)$ is a convex function. The convergence rate of Equations (32) and (33) is as follows. If Assumptions 1-4 hold, let $\alpha (t) = 1/\sqrt {t + 1} $ , where $t \geqslant 0$ . For all $i,i = 1,2, \cdot \cdot \cdot, n$ , then

(41) \begin{align}& {E\left[ {f\left( {{{\tilde \psi }_i}\left( T \right)} \right) - f\left( {{\psi ^{\rm{*}}}} \right)} \right] \le \dfrac{n}{{2\sqrt T }}{{\left\| {\bar \psi \left( 0 \right) - {\psi ^{\rm{*}}}} \right\|}^2} + \dfrac{{C\left( {3L + \mathop \sum \nolimits_{i = 1}^n \left( {{L_i} + {\nu _i}} \right)} \right)\mathop \sum \nolimits_{j = 1}^n {\psi _j}{{\left( 0 \right)}_1}}}{{\left( {1 - \lambda } \right)\sqrt T }}}\nonumber\\[7pt]& { + \dfrac{{C\left( {9L + \mathop \sum \nolimits_{i = 1}^n \left( {{L_i} + {\nu _i}} \right)} \right){\rm{log}}\left( {T + 1} \right)}}{{\left( {1 - \lambda } \right)\sqrt T }}\mathop \sum \limits_{j = 1}^n 2\sqrt d \left( {{L_j} + {\nu _j}} \right) + \dfrac{{3\left( {1 + {\rm{log}}\left( {T + 1} \right)} \right)}}{{\sqrt T }}\mathop \sum \limits_{i = 1}^n {{\left( {{L_i} + {\nu _i}} \right)}^2}}\end{align}

Equation (41) is established with probability 1 [Reference Xu, Zhu and Wu33]. Where, ${L_i}$ , ${\nu _i}$ , $T$ are positive constants. $\lambda \in \left( {0,1} \right)$ , $L = \mathop \sum \limits_{i = 1}^n {L_i}$ . $T$ is the number of iterations. ${\bar \psi _i}\left( T \right)$ is defined as

(42) \begin{align}{\bar \psi _i}\left( T \right) = \frac{1}{{\mathop \sum \nolimits_{t = 0}^T \alpha {{(t)}^{t = 0}}}}\mathop {\mathop \sum \limits_{} }\limits^T {\rm{\;}}\alpha (t){\psi _i}(t)\end{align}

The proposed DSSPA has asymptotic convergence by selecting an appropriate step size, Equations (32) and (33) asymptotically converge to the optimal solution. Through the optimisation of the DSSPA with a layered constraint set, the optimal detection array of the bearing-only sensor network can be calculated, and the convergence rate of the algorithm is $O\left( {{\rm{log}}T/T} \right)$ .

6.0 Simulation experiments

In order to verify the effectiveness and superiority of the proposed DSSPA, three simulation experiments of UAV formation for bearing-only collaborative detection are set up. UAV formations track and locate the single target. As shown in Fig. 8, the communication topology networks for the three sets of simulation experiments are depicted. The numbers represent the IDs of the UAVs, and the connecting lines indicate the communication established between UAVs.

Figure 8. Network topologies of three groups of simulation experiments.

It is assumed that the bearing-only sensors of all UAVs complete the allocation of time and space uniformly. ${S_i}$ corresponds to the spatial position of the UAVn in the bearing-only sensor network. The number of UAV formations in the three groups of simulation experiments increases successively to comprehensively validate the models and optimisation algorithms constructed for small-scale, medium-scale, and large-scale scenarios. ${X_T}$ = (10,000m, 10,000m, 10,000m) is the initial state of the target. The fixed baseline is set to $l = 1{\rm{km}}$ , and the safe flight altitude is set to ${\rm{\Delta }}H = 1.5{\rm{km}}$ .

In the first simulation experiment, the network topology of the simulation experiment is shown in Fig. 8(a). The formation of six UAVs is divided into the same layer to track and locate a single target according to the constraints of safe flight altitude. The initial parameters of the first simulation experiment are shown in Table 1.

Table 1. Initial simulation parameters of the first experiment

As the simulation time increases, each UAV in the formation converges to a stable value, as shown in Fig. 9. Sorting from large to small according to the time spent on convergence, the order is UAV3 (6.82s) $ \gt $ UAV6 (6.68s) $ \gt $ UAV1 (6.00s) $ \gt $ UAV4 (5.35s) $ \gt $ UAV5 (4.76s) $ \gt $ UAV2 (4.32s). The brackets indicate the local correlation function value $J_i^{ - 1}$ of each UAV at that moment when convergence is stable, expressed as the reciprocal of the true value of the local GDOP objective function.

Figure 9. Local correlation function value $J_i^{ - 1}$ in the first scenario.

Each UAV in the formation not only calculates the local objective function, but also conducts distributed network communication to obtain global information interaction and unification. As the simulation time increases, the correlation function values $J_i^{ - 1}$ of the global indicators of each UAV gradually converge, as shown in Fig. 10.

Figure 10. Global correlation function value $J_i^{ - 1}$ in the first scenario.

At time 9.85s for the entire UAV formation, the bearing-only detection system consistently converged to the correlation function value of the global indicator $J_i^{ - 1} = 2.72$ . The optimal detection formation planning of UAV formation are as follows: ${S_1}$ (35,045, 10,512, 10,512), ${S_2}$ (35,045, 11,207, 10,005), ${S_3}$ (35,045, 10,512, 95,025), ${S_4}$ (35,045, 24,000, 95,025), ${S_5}$ (35,045, 8,792, 10,005), ${S_6}$ (35,045, 95,025, 10,512).

In the array planning process of the first experiment, a comparative analysis is conducted under different angular measurement error conditions, as shown in Fig. 11. When the angle measurement errors are 0.01, 0.02 and 0.05rad, respectively, the related objective functions of the system can converge to 5.05, 3.96 and 2.82 in sequence. It can be seen intuitively that in a scenario where the number of UAV formations is small, as the UAV formation angle measurement error increases, the limit of the proposed DSSPA optimal planning also decreases. At the same time, it also verifies the effectiveness of the objective function for optimising the formation array provided in this paper. The fundamental factor that affects the optimal detection array for UAV formation still comes from the propagation of angle measurement errors through geometric configurations.

On the basis of ensuring the objective function and constraint set of the optimised array, the convergence speed and operation processing capabilities of the three methods: DSSPA, PSO and MA are statistically analysed, as shown in Fig. 11. In this case where the number of UAV formations is small, DSSPA (10.2), PSO (9.3s) and MA (7.8s) can all converge to the optimal results. But overall, the convergence speed of PSO is 8.82% higher than the DSSPA, and the convergence speed of PSO is 11.86% higher than the DSSPA. As centralised optimisation algorithms, PSO and MA have faster convergence speed and are more suitable for this scenario. However, it is worth noting that the distributed structure of DSSPA inherently provides excellent resistance to interference and robustness.

In the second simulation experiment, the network topology of the second simulation experiment is shown in Fig. 8(b). The formation of nine UAVs are divided into two layers to track and locate a single target according to the constraints of safe flight altitude. The initial parameters of the second simulation experiment are shown in Table 2.

Table 2. Initial simulation parameters of the second experiment

Figure 11. Comparison of different methods and $\sigma $ in the first scenario.

In the second scenario, as the simulation time increases, the changes in the local correlation function values of each agent in the UAV formation are shown in Fig. 12. In the first layer, sorting from large to small according to the time spent on convergence, the order is: UAV7(5.85s) $ \gt $ UAV8 (5.24s) $ \gt $ UAV9 (4.85s). In the second layer, the order is: UAV5 (7.63s) $ \gt $ UAV3 (7.42s) $ \gt $ UAV2(6.92s) $ \gt $ UAV1 (6.82 s) $ \gt $ UAV6 (5.85s) $ \gt $ UAV4 (4.36s).

Figure 12. Local correlation function value $J_i^{ - 1}$ in the second scenario.

As the simulation time increases, the correlation function values of the global indicators of each UAV gradually converge, as shown in Fig. 13. At 9.85s for the entire UAV formation, the system consistently converged to the correlation function value of the global indicator $J_i^{ - 1} = 2.68$ . The optimal detection formation planning of UAV formation are as follows: ${S_1}$ (36,065, 10,524, 10,524), $S2$ (36,065, 11,018, 10,012), ${S_3}$ (36,065, 10,508, 95,035), ${S_4}$ (36,065, 24,020, 95,035), ${S_5}$ (36,065, 8,792, 10,012), ${S_6}$ (36,065, 95,035, 10,524), ${S_7}$ (34,565, 10,005, 10,433), ${S_8}$ (34,565, 10,500, 9,567), ${S_9}$ (34,565, 9,500, 9,567).

Figure 13. Global correlation function value $J_i^{ - 1}$ in the second scenario.

In the array planning process of the second experiment, the changes in the convergence situation of DSSPA under different angle measurement errors are also analysed, as shown in Fig. 14.

Figure 14. Comparison of different methods and $\sigma $ in the second scenario.

When the measurement errors are 0.01, 0.02 and 0.05rad respectively, the related objective functions of the system can converge to 4.87, 4.12 and 2.75 in sequence. It can also be seen intuitively that as the UAV formation angle measurement error increases, the limit of the proposed DSSPA optimal planning also decreases. In the case of this medium number of UAV formations, the three methods of DSSPA (10.6s), PSO (11.6s) and MA (11.2s) can all converge to the optimal results. As a distributed algorithm, the convergence speed of DSSPA has begun to show its advantages in this scenario, and it can obtain the optimal solution at a faster convergence speed. Overall, the convergence speed of DSSPA is increased by 8.62% compared with PSO. Compared with MA, the convergence speed of DSSPA is increased by 5.35%.

In the third simulation experiment, the network topology of the third simulation experiment is shown in Fig. 8(c). The formation of 19 UAVs was divided into three layers to track and locate a single target according to the constraints of safe flight altitude. The initial parameters of the third simulation experiment are shown in Table 3.

Table 3. Initial simulation parameters of the second experiment

In the third scenario, the changes in the local correlation function values of each agent in the UAV formation are shown in Fig. 15. In the first layer, sorting from large to small according to the time spent on convergence, the order is UAV1 (5.85s) $ \gt $ UAV3 (5.24s) $ \gt $ UAV2 (4.85s). In the second layer, the order is: UAV5 (4.36s) $ \gt $ UAV4 (4.82s) $ \gt $ UAV6 (5.23s) $ \gt $ UAV7 (6.0s) $ \gt $ UAV9 (6.64s) $ \gt $ UAV8(6.87 s). In the third layer, the order is UAV10(4.27s) $ \gt $ UAV11 (4.66s) $ \gt $ UAV12 (5.31s) $ \gt $ UAV13 (6.12s) $ \gt $ UAV14 (6.85s) $ \gt $ UAV15 (6.91s) $ \gt $ UAV16 (7.15s) $ \gt $ UAV17 (7.54s) $ \gt $ UAV18 (7.56s) $ \gt $ UAV19 (7.82s).

Figure 15. Local correlation function value $J_i^{ - 1}$ in the third scenario.

As the simulation time increases, the correlation function values of the global indicators of each UAV gradually converge, as shown in Fig. 16. This paper considers the layered constraints of the UAV formation. Inter-layer communication does not constitute an observation unit, but is only for the convenience of measuring the overall array planning effect of the UAV formation. The array planning of the first layer of UAVs converged at 7.56s, and then the array planning of the second layer of UAVs converged at 9.62s. At this time, the first layer of UAVs and the second layer of UAVs conduct inter-layer communication to obtain the local objective function value $J_i^{ - 1} = 3.46$ . The formation planning of the third layer of UAVs converges at 12.68s, and after the formation is determined, inter-layer communication is performed with the second layer to obtain the global objective function value $J_i^{ - 1} = 2.25$ of the UAV formation. The optimal detection array planning of UAV formation are as follows: ${S_1}$ (36,560, 10,005, 10,433), ${S_2}$ (36,560, 10,500, 9,567), ${S_3}$ (36,560, 9,500, 9,567), ${S_4}$ (38,060, 10,512, 10,512), ${S_5}$ (38,060, 11,207, 10,005), ${S_6}$ (38,060, 10,512, 95,025), ${S_7}$ (38,060, 24,000, 95,025), ${S_8}$ (38,060, 8,792, 10,005), ${S_9}$ (38,060, 95,025, 10,512), ${S_{10}}$ (39,560, 10,707, 11,207), ${S_{11}}$ (39,560, 11,514, 10,516), ${S_{12}}$ (39,560, 11,514, 9,528), ${S_{13}}$ (39,560, 10,707, 8,793), ${S_{14}}$ (39,560, 10,032, 8,086), ${S_{15}}$ (39,560, 9,293, 8,793), ${S_{16}}$ (39,560, 8,586, 9,504), ${S_{17}}$ (39,560, 8,586, 10,500), ${S_{18}}$ (39,560, 9,293, 11,207), ${S_{19}}$ (39,560, 10,005, 11,914).

Figure 16. Global correlation function value $J_i^{ - 1}$ in the third scenario.

When the angle measurement errors are 0.01, 0.02 and 0.05rad, respectively, the related objective functions of the system can converge to 4.36, 3.82 and 2.64 in sequence. It can be seen intuitively that as the UAV formation angle measurement error increases, the limit of the proposed DSSPA optimal planning also decreases, as shown in Fig. 17. In the case of such a large number of UAV formations, the three methods of DSSPA (12.8s), PSO (15.2s) and MA (14.6s) can all converge to the optimal results. DSSPA can obtain the optimal solution with the fastest convergence speed in third scenario. Compared with PSO, the convergence speed of DSSPA is increased by 15.79%. Compared with MA, the convergence speed of DSSPA is increased by 14.06%. The research object of this paper is also medium- and large-scale UAV formations. It is obvious that DSSPA has superior convergence performance.

Figure 17. Comparison of different methods and $\sigma $ in the second scenario.

As the number of UAVs $N$ increases, the system network with different angle measurement errors will steadily approach the optimal observation geometry, and the GDOP at this time is closest to $GDO{P_{{\rm{min}}}}$ . $GDO{P_{{\rm{min}}}}$ represents the optimal formation array, and the position of each bearing-only UAV is the most favourable relative to the geometric distribution of target positioning, thereby providing the best positioning accuracy. Because within a fixed sampling period, more bearing-only UAVs can provide more measurement data, making the measurement accuracy of the system closer and closer to the optimal value.

For the DSSPA convergence curve proposed in this paper, as the angle measurement error increases, the $GDO{P_{min}}$ of the bearing-only UAV formation with the same number of nodes becomes larger. The accuracy of angle measurement directly affects the tracking and positioning effect of the system.

When the local cost function of each UAV is a convex function, the UAV can only obtain its own noise subgradient information. The DSSPA can asymptotically converge to the optimal solution. It is proved that the algorithm can converge asymptotically to the optimal solution at the iteration rate $O\left( {{\rm{log}}T/\sqrt T } \right)$ . In the centralised algorithm, assuming that the cost function of each UAV is convex and differentiable, the variance of the noise subgradient is bounded, the convergence rate of this type of algorithm is $O\left( {1/T} \right)$ [Reference Chen, Fu, Zhang, Fang and Xiao34]. In small-scale UAV formation scenarios, the centralised optimisation algorithm may have a faster convergence rate, but the convergence rate of the proposed DSSPA is already very close to this convergence rate. However, in large-scale UAV formation scenarios, the convergence rate of centralised optimisation algorithms is affected by factors such as communication overhead and computing resources. DSSPA can better adapt to large-scale problems and distributed environments. In addition, in the proposed DSSPA, there is no need to assume that the local cost function is differentiable.