Joint Adaptive Sampling Interval and Power Allocation for Maneuvering Target Tracking in a Multiple Opportunistic Array Radar System

In this paper, a joint adaptive sampling interval and power allocation (JASIPA) scheme based on chance-constraint programming (CCP) is proposed for maneuvering target tracking (MTT) in a multiple opportunistic array radar (OAR) system. In order to conveniently predict the maneuvering target state of the next sampling instant, the best-fitting Gaussian (BFG) approximation is introduced and used to replace the multimodal prior target probability density function (PDF) at each time step. Since the mean and covariance of the BFG approximation can be computed by a recursive formula, we can utilize an existing Riccati-like recursion to accomplish effective resource allocation. The prior Cramér-Rao lower boundary (prior CRLB-like) is compared with the upper boundary of the desired tracking error range to determine the adaptive sampling interval, and the Bayesian CRLB-like (BCRLB-like) gives a criterion used for measuring power allocation. In addition, considering the randomness of target radar cross section (RCS), we adopt the CCP to package the deterministic resource management model, which minimizes the total transmitted power by effective resource allocation. Lastly, the stochastic simulation is embedded into a genetic algorithm (GA) to produce a hybrid intelligent optimization algorithm (HIOA) to solve the CCP optimization problem. Simulation results show that the global performance of the radar system can be improved effectively by the resource allocation scheme.


Background and Motivation
The opportunistic array radar (OAR), proposed by the United States Naval Postgraduate School (NPS), is a new system radar for the next generation naval stealth destroyer DD(X) [1][2][3]. In the OAR system, the stealth of the platform is taken as the core and the digital array is regarded as the base, and then the array elements and the transmit/receive (T/R) modules are placed arbitrarily and aperiodically at available open areas over the entire 3-D space of carrier platforms [4,5]. Due to the unique structure Sensors 2020, 20 of the antenna array, the resource management of OAR is flexible. On the other hand, the maneuvering target tracking (MTT) plays a vital part for various commercial and military applications and receives plenty of attention [6][7][8][9]. For example, the application areas include battlefield surveillance, air traffic control, air defense, and fire control. Hence, the resource management for MTT in the OAR system is a significant and worthwhile research direction. So far, many studies has been done to enhance the utilization efficiency of a scarce radar resource in the target tracking process [10][11][12][13][14][15][16]. Reference [10] proposes a sensor selection policy for static target location in distributed multiple-radar architectures. This policy is formulated in a combinatorial optimization framework as a knapsack problem (KP) to obtain a performance level with the lowest cost. Referring to the Bayesian Cramér-Rao lower bound (BCRLB) as a criteria, Reference [11] puts forward a joint power and beams allocation algorithm for multiple target tracking in a co-located multiple input multiple output (MIMO) radar system. Besides, by incorporating an information reduction factor (IRF), the BCRLB is efficiently computed with a measurement origin uncertainty (MOU) and the radar resource is allocated to the targets reasonably in clutter [12]. A joint node selection and power allocation strategy is developed with the objective of tracking multiple targets in Reference [15]. In order to improve the worst-case tracking accuracy with multiple targets, the proposed mechanism implements the optimal resource allocation based on the feedback information in the tracking recursion cycle. In Reference [16], the resource allocation problem concerns the sensor subset, power, and bandwidth. The three-variable optimization problem is simplified by deriving the relationship between the optimal power allocation vector and the bandwidth allocation vector. This algorithm can help achieve better performance within the same resource constraints.
The previously mentioned work provides us with an opportunity, and helps to deal with the resource management problem. However, this work only considers the resource management problems of a uniform moving target [11][12][13][14][15][16]. In the actual battlefield environment, such a favorable circumstance is always an unrealistic assumption. In military applications, the target aircrafts typically fly intelligently in order to achieve their objectives. Even they have electronic equipment that can measure the active sensor's attempts to track them and may respond accordingly. Therefore, an intelligent resource allocation framework is presented herein for tracking a maneuvering target with Markovian switching dynamics (MSD means that the current target state only depends on the previous target states) in multiple OAR systems. This is a hybrid estimation problem where it is required to estimate sequentially both the continuous-valued target state and the discrete-valued switching model variable. There are several existing error bounds that deal with the hybrid system estimation problem. Examples are Bhattacharyya, Bobrovsky-Zakai, and Weiss-Weinstein lower bounds [17]. However, the derivations of these bounds are extremely complex and their implementation is very involved [18,19]. In addition, the direct application of the BCRLB recursive formula for nonlinear filtering [20] would lead to differentiation of terms involving the discrete-valued model variable. Hence, an approximation to the BCRLB for MSD is considered in this case [21][22][23]. The essence of this technique is to replace the multimodal prior target probability density function (PDF) with a best-fitting Gaussian (BFG) distribution at each time step. The mean and covariance of BFG distribution can be computed through a recursive formula.
In a traditional multiple radar system, prior knowledge is not utilized for a resource allocation scheme. In general, the sampling interval is fixed and the transmitted power of spatially separate radars is equal. Such an operation will definitely result in the inadequate use of the limited time and power resources. In terms of time resource, if the sampling interval is too long, it may lead to a large tracking error, or even lose the target. If the sampling interval is too short, the frequent beam illumination is a waste of the radar resource. Through the BFG approximation, we could determine when to transmit beams to acquire data of the tracked targets. The prior Cramér-Rao lower boundary (prior CRLB-like) can be set as a criterion for an adaptive sampling interval [11]. As a consequence, more tracking tasks can be maintained and more radar resources can be applied to track initiation [24]. Regarding the power resource, based on the prior CRLB-like, the BCRLB-like can be calculated as an error measurement contributing to the intelligent power allocation [20].
Generally speaking, in the traditional work of radar resource management, the radar cross section (RCS) of the illuminated target is considered to be a determined value at each sampling instant [11][12][13]15,16]. Due to the effects of the identity, attitude, position, aspect angle, wave length, and polarization, the RCS of the next sampling instant is uncertainty [25], i.e., the RCS to be estimated obeys a certain distribution. To adequately disclose the uncertainty of RCS, we regard it as a random variable. On account of the randomness of constraint conditions, we introduce chance-constraint programming (CCP) to handle the uncertainty by guaranteeing that the stochastic constraint holds at a specific confidence level [26]. When the confidence level is equal to one, we need to provide a radar resource than in other cases to satisfy the constraints. Yet, it can achieve better tracking performance. However, it is not cost-effective to spend superfluous radar resource on the extremely low probability incident that the RCS takes values at the lower boundary of the distribution range. In addition, to obtain an acceptable tracking error, the confidence level cannot be low either. The confidence level is decided by the environments and desired tracking performance. Therefore, the tracking accuracy could be adjusted flexibly conditioned on different confidence levels. Since we abandon the extreme case that satisfies the constraints conditioned on a very low probability, the resource consumption is reduced strikingly to maintain more tracks.
Through the previously mentioned analysis, we propose a joint adaptive sampling interval and power allocation scheme based on CCP for MTT using BFG in multiple OAR systems. The BFG distribution is presented to replace the multimodal prior target PDF at each time step. In order to realize an adaptive sampling interval, the prior CRLB-like is calculated to decide when to illuminate the target. In conjunction with the data Fisher information matrix (FIM) obtained from distributed radars, the approximation of BCRLB is derived for a maneuvering target with MSD, and it is named BCRLB-like. Conditioned on a specified confidence level, we adopt the CCP to package the deterministic resource management model as the final uncertain model. The whole algorithm can be seen as an intelligent response of a radar system to MTT by perceiving environments. The BCRLB-like is to measure the error for target state estimation, and it provides us a criterion to pre-allocate the radar resource. The optimization problem can be solved by a hybrid intelligent optimization algorithm (HIOA) integrating a stochastic simulation and genetic algorithm (GA).

Main Contributions and Innovations
The main contributions and innovations of this paper are as follows.
(1) The BFG approximation is employed so as to allocate the radar resource conveniently. Due to the diversities of motion states of maneuvering targets, it is difficult to predict the target state, and the radar resource allocation has no referenced criterion at the next sampling instant. Guaranteeing that the state vectors have the same mean and covariance under different models, the BFG approximation is introduced to replace the multimodal prior target PDF at each time step. The target state can be predicted easily by a single motion equation, and then accomplish the resource pre-allocation. (2) A joint adaptive sampling interval and power allocation scheme is presented. The prior CRLB-like as a measurement criteria is compared with the upper boundary of the given tracking error threshold to determine the next optimal sampling interval. The tracking BCRLB-like is computed for power allocation among the distributed radars. The diagonal elements of BCRLB-like provide a referenced boundary on the variances of the estimation of the target's hybrid state. (3) The CCP is brought in to handle the uncertainty of target information in a resource management model. The target RCS is regarded as a random variable. The CCP balances the radar resource and the tracking performance by adjusting the confidence level. If the target environment is simple or the tracking performance requirement is low, the confidence level could be lowered appropriately to save more resources for other tasks. The remainder of this paper is structured as follows. In Section 2, we provide the system model for a linear jump Markov system. In Section 3, the CCP model is formulated. By introducing the BFG approximation, the BCRLB-like is derived as the inverse of the Bayesian information matrix (BIM) for MTT. This model integrates the adaptive sampling interval and power allocation in multiple OAR systems. Section 4 proposes the HIOA to solve the model, and two algorithms of target state estimation are given in this section. The simulation results and corresponding analysis are provided in Section 5 to verify the effectiveness of the proposed resource allocation scheme. Lastly, we conclude this paper in Section 6.

System Model
Suppose a two-dimensional multiple radar system including M (M ≥ 2) distributed monostatic radars. The qth radar is located at (x q , y q ) (q = 1, 2, . . . , M). We assume that: (1) the carrier frequency of each radar is different. (2) There is only a matched filter for the respective transmitted signal. In terms of the two assumptions, each radar receives all the echo signals from the target, but the exclusive matched filter can filter out the echo signals transmitted by the other radars. Hence, each radar operates in a monostatic way. In addition, the centralized tracking is adopted for the above multiple radar systems. All radars generate measurements at each sampling instant and report those measurements to the central fusion center (CFC) [27]. It fuses all the acquired measurements in turn and updates the tracks. The fusion architecture is in Figure 1.
Sensors 2020, 20, x FOR PEER REVIEW 4 of 24 estimation are given in this section. The simulation results and corresponding analysis are provided in Section 5 to verify the effectiveness of the proposed resource allocation scheme. Lastly, we conclude this paper in Section 6.

System Model
Suppose a two-dimensional multiple radar system including M (M ≥ 2) distributed monostatic radars. The qth radar is located at (xq, yq) (q = 1, 2, …, M). We assume that: (1) the carrier frequency of each radar is different. (2) There is only a matched filter for the respective transmitted signal. In terms of the two assumptions, each radar receives all the echo signals from the target, but the exclusive matched filter can filter out the echo signals transmitted by the other radars. Hence, each radar operates in a monostatic way. In addition, the centralized tracking is adopted for the above multiple radar systems. All radars generate measurements at each sampling instant and report those measurements to the central fusion center (CFC) [27]. It fuses all the acquired measurements in turn and updates the tracks. The fusion architecture is in Figure 1. where Ri (i = 1,2,…M) denotes the qth radar. Hence, in accordance with the fusion architecture, the block diagram of a close-loop intelligent tracking system is given in Figure 2.  The tracked target is initially located at (xT0, yT0), and the target is located at (xTk, yTk) at the kth sampling interval. For convenience, the (xTk, yTk) is denoted as (xk, yk).

Signal Model
It is assumed that the transmitted signal of the qth radar to the target at the kth sample interval is sq,k(t): where Pq,k denotes the transmitted power from the qth radar to the target at the kth sample interval, and f c q is the carrier frequency. Sq,k(t) denotes the complex envelope of the transmitted signal. The baseband representation of the echo signal for the qth radar at the kth sample interval is: where R i (i = 1,2, . . . M) denotes the qth radar. Hence, in accordance with the fusion architecture, the block diagram of a close-loop intelligent tracking system is given in Figure 2.
Sensors 2020, 20, x FOR PEER REVIEW 4 of 24 estimation are given in this section. The simulation results and corresponding analysis are provided in Section 5 to verify the effectiveness of the proposed resource allocation scheme. Lastly, we conclude this paper in Section 6.

System Model
Suppose a two-dimensional multiple radar system including M (M ≥ 2) distributed monostatic radars. The qth radar is located at (xq, yq) (q = 1, 2, …, M). We assume that: (1) the carrier frequency of each radar is different. (2) There is only a matched filter for the respective transmitted signal. In terms of the two assumptions, each radar receives all the echo signals from the target, but the exclusive matched filter can filter out the echo signals transmitted by the other radars. Hence, each radar operates in a monostatic way. In addition, the centralized tracking is adopted for the above multiple radar systems. All radars generate measurements at each sampling instant and report those measurements to the central fusion center (CFC) [27]. It fuses all the acquired measurements in turn and updates the tracks. The fusion architecture is in Figure 1. where Ri (i = 1,2,…M) denotes the qth radar. Hence, in accordance with the fusion architecture, the block diagram of a close-loop intelligent tracking system is given in Figure 2.  The tracked target is initially located at (xT0, yT0), and the target is located at (xTk, yTk) at the kth sampling interval. For convenience, the (xTk, yTk) is denoted as (xk, yk).

Signal Model
It is assumed that the transmitted signal of the qth radar to the target at the kth sample interval is sq,k(t): where Pq,k denotes the transmitted power from the qth radar to the target at the kth sample interval, and f c q is the carrier frequency. Sq,k(t) denotes the complex envelope of the transmitted signal. The baseband representation of the echo signal for the qth radar at the kth sample interval is: The tracked target is initially located at (x T0 , y T0 ), and the target is located at (x Tk , y Tk ) at the kth sampling interval. For convenience, the (x Tk , y Tk ) is denoted as (x k , y k ).

Signal Model
It is assumed that the transmitted signal of the qth radar to the target at the kth sample interval is s q,k (t): where P q,k denotes the transmitted power from the qth radar to the target at the kth sample interval, and f c q is the carrier frequency. S q,k (t) denotes the complex envelope of the transmitted signal. The baseband representation of the echo signal for the qth radar at the kth sample interval is: where h q,k denotes the target RCS for the qth radar, and it is a random variable [25]. Furthermore, α q,k ∝1/R 4 q,k denotes the variation of the signal strength due to path loss effects along the path of radar q-target-radar q. τ q,k means time delay and f q,k is the Doppler frequency. q,k (t) represents a zero-mean, complex Gaussian white noise with the autocorrelation function σ 2 δ(τ).

Target Motion Model
As described in previous content, we restrict ourselves to linear jump Markovian systems with an additive Gaussian process noise. Herein, we consider three possible motion models. The three motion models are constant velocity (CV), constant acceleration (CA), and coordinated turn (CT) with a known turn rate, respectively. A larger set of possible motions would be more realistic. Nevertheless, the current set is sufficient to verify the advantage of this algorithm. An overview of other target motion models may be found in References [28,29]. In addition, the dimensionalities of target states for different dynamic models are different, and it may result in a change in the dimensionality of the information matrix, whenever there is a model switch [30]. Hence, to avoid this problem, we set the state vectors with the same dimensionality in different dynamic models. The target state vector consists of the position, velocity, and acceleration in a two-dimensional space. Therefore, the dimensionality of the state vector is six. Details of each motion model used in this paper are given as follows.

Constant Velocity Motion Model
The CV motion model in the xy plane is shown below.
where ξ k = x k .
x k 0 y k . y k 0 T denotes the state vector, and its dimension is N ξ = 6. [x k y k ] T and .
x k . y k T denote the position and velocity of the target, respectively. The acceleration is [0, 0] T . [•] T denotes matrix transposition. F 1,k−1 is the 6 × 6 transition matrix.
where ⊗ is the Kronecker operator. I 2 denotes an identity matrix of order 2. T k is the sampling interval (Let k = 0 when the initial time is t = 0, and the first sampling interval T 1 is between k = 0 and k = 1. Thereby, the kth sampling interval between epochs k−1 and k is denoted as T k in (4). The expression has the same meaning in the remaining paper). w 1,k−1 represents a zero-mean, complex Gaussian white noise, and its covariance is shown below.

Constant Acceleration Motion Model
The CA motion model is prescribed as follows. ..
x k y k . y k .. y k T denotes the state vector. ..
x k .. y k T denotes the constant acceleration. F 2,k−1 is the 6 × 6 transition matrix.
w 2,k−1 is a zero-mean, complex Gaussian white noise with the covariance.
where σ 2 is the process noise intensity of CA [31].

Coordinated Turn Motion Model
In this scenario, we assume that the angular turn rate Ω is known and constant. Hence, the target motion model of CT remains linear. Then, the CT motion model is expressed as: where ξ k = x k .
x k 0 y k . y k 0 T denotes the state vector. F 3,k−1 is the 6 × 6 transition matrix.
where Ω is the angular turn rate, and w 3,k−1 is a zero-mean, complex Gaussian white noise with the covariance given by the equation below.
where σ 3 is the process noise intensity of CT [31].

Measurement Model
The measurement equation does not play a role in the BFG approximation. Hence, although the BFG approximation is restricted to linear switching dynamic models, it places no restriction on the measurement equation, which can be nonlinear. The received signal from the target is an attenuated version of the transmitted signal. The range, azimuth, and Doppler frequency can be extracted from the received signal. The nonlinear measurement equation of the qth radar can be written by the equation below. where corresponding to range, azimuth, and Doppler shift. λ q is the carrier wavelength of the qth radar. The measurement error u q,k is a zero-mean Gaussian white noise with the covariance.
, and σ 2 f q,k are the BCRLBs of the mean square error (MSE) of the range, azimuth, and Doppler shift at high signal-to-noise ratio (SNR). The BCRLBs for σ 2 R q,k and σ 2 f q,k are given in Chapter 10 of Reference [31], and the BCRLB for σ 2 θ q,k which follows from Reference [32]. The BCRLBs can be described as follows.
where B q,k and T q,k are the effective bandwidth and effective time duration [31], respectively. B NN is the null-to-null beam width of the receiver antenna. c is the speed of light. SNR q,k is the SNR denoted as follows [33].

Resource Management Model for Maneuvering Target Trackinĝ
ξ k is an estimation of the target state ξ k . Hence,ξ k is a function of the measurement value Z k , which is represented as follows.
The performance of any estimatorξ k is measured by the MSE matrix. In the Bayesian estimation, the MSE can be bounded by the BCRLB C k BCRLB (defined to be the inverse of BIM J(ξ k )) in estimating the random vector ξ k , under suitable regularity conditions. For the target with a single motion model, the MSE of any estimator cannot be less than the BCRLB C k BCRLB [17].
where E ξ k ,Z k [·] denotes the expectation over ξ k and Z k Sensors 2020, 20, 981 The BIM J(ξ k ) can be described as: The joint PDF p(Z k ,ξ k ) can be factorized as the product of PDF p(ξ k ) and conditional PDF p(Z k |ξ k ).
Hence, the BIM J(ξ k ) can be divided into a prior information matrix J P (ξ k ) and a data information matrix J D (ξ k ) [17].
The prior information matrix J P (ξ k ) is determined by the motion equation of targets, and the transmitted power has no impact on it. On the contrary, the data information matrix J D (ξ k ) is affected by the transmitted power allocated. The larger the transmitted power is, the larger the J D (ξ k ) is.

Best-Fitting Gaussian Approximation
To calculate the BIM for linear jump Markovian systems with additive Gaussian process noise through Equation (22), we need to express the dynamics of the system.
with the BFG approximation: where r k (r k = 1, 2, . . . , N m and N m is the number of the target motion models) specifies the target motion model in effect during the time interval between sampling time indexes k−1 and k. Φ k−1 , W k−1 , and Q p,k−1 are the transition matrix, process noise vector, and covariance matrix of the BFG approximation equation, respectively. For convenient representation in the following process of derivation, the notation F i,k and Q i,k in Section 2.2 (i = 1, 2, . . . , N m ) are changed to F r k k−1 and Q r k p,k−1 . is an additive Gaussian process noise with W At each sampling interval, the regime-dependent motion model (given by Equation (23) and known as "model 1") is replaced with a single BFG distribution (denoted by Equation (24) and referred to as "model 2"). Φ k−1 and Q p,k−1 (Q p,k−1 ≥ 0) are calculated to satisfy the requirements that the distribution of ξ k has the same mean and covariance under different models at each stage.
where E[•] and Cov[•] denote the expectation operator and covariance operator, respectively. The computation procurement of Φ k−1 is given in Appendix A. Then, with the aid of Φ k−1 , we obtain Q p,k−1 by the detailed derivation in Appendix B.
It is assumed that the initial distribution of state vector ξ is Gaussian with mean ξ 0 and covariance J −1 P (ξ 0 ). Generally speaking, initially, the target is in a CV motion. Hence, the ε 1 and Sensors 2020, 20, x FOR PEER REVIEW 8 of 24 Hence, the BIM J(ξk) can be divided into a prior information matrix JP(ξk) and a data information matrix JD(ξk) [17].
The prior information matrix

Best-Fitting Gaussian Approximation
To calculate the BIM for linear jump Markovian systems with additive Gaussian process noise through Equation (22), we need to express the dynamics of the system.

( )
with the BFG approximation: where rk (rk = 1, 2, …, Nm and Nm is the number of the target motion models) specifies the target motion model in effect during the time interval between sampling time indexes k-1 and k. Φk-1, Wk-1, and Qp,k-1 are the transition matrix, process noise vector, and covariance matrix of the BFG approximation equation, respectively. For convenient representation in the following process of derivation, the notation Fi,k and Qi,k in Section 2.2 (i = 1, 2, …, Nm) are changed to 1 k r k − F and p, 1 At each sampling interval, the regime-dependent motion model (given by Equation (23) and known as "model 1") is replaced with a single BFG distribution (denoted by Equation (24) and referred to as "model 2"). Φk-1 and Qp,k-1 (Qp,k-1 ≥ 0) are calculated to satisfy the requirements that the distribution of k ξ has the same mean and covariance under different models at each stage.
The computation procurement of Φk-1 is given in Appendix A. Then, with the aid of Φk-1, we obtain Qp,k-1 by the detailed derivation in Appendix B.
Φk-1 is given in Appendix A. Then, with the aid of Φk-1, we n Appendix B.
Step 3. By derivation, obtain Φ k−1 of BFG distribution as (see Appendix A): Step 4. With the aid of Φ k−1 , derive the covariance matrix ively. For convenient representation in the following process of derivation, the i,k in Section 2.2 (i = 1, 2, …, Nm) are changed to 1 k r k − F and p, 1 pling interval, the regime-dependent motion model (given by Equation (23) and l 1") is replaced with a single BFG distribution (denoted by Equation (24) and odel 2"). Φk-1 and Qp,k-1 (Qp,k-1 ≥ 0) are calculated to satisfy the requirements that the k has the same mean and covariance under different models at each stage.
denote the expectation operator and covariance operator, respectively.
ation procurement of Φk-1 is given in Appendix A. Then, with the aid of Φk-1, we e detailed derivation in Appendix B.
d that the initial distribution of state vector ξ is Gaussian with mean 0 ξ and . Generally speaking, initially, the target is in a CV motion. Hence, the ε1 and 1 easily. The transition probability is pij  To rsive computation procedure for Φk- 1 and Qp,k-1 is given as follow: 1, and initialize εk-1, k-1, and the mode probabilities pk-1(r), r = 1, 2, …, Nm late the mode probabilities.
dexes k-1 and k. Φk-1, Wk-1, and Qp,kmatrix of the BFG approximation llowing process of derivation, the is an additive odel (given by Equation (23) and on (denoted by Equation (24) and to satisfy the requirements that the erent models at each stage.
A. Then, with the aid of Φk-1, we is Gaussian with mean 0 ξ and CV motion. Hence, the ε1 and 1 To is given as follow: lities pk-1(r), r = 1, 2, …, Nm Step 5. Calculate Q p,k−1 as follows: he number of the target motion models) specifies the target motion rval between sampling time indexes k-1 and k. Φk-1, Wk-1, and Qp,ks noise vector, and covariance matrix of the BFG approximation nient representation in the following process of derivation, the e regime-dependent motion model (given by Equation (23) and with a single BFG distribution (denoted by Equation (24) and Qp,k-1 (Qp,k-1 ≥ 0) are calculated to satisfy the requirements that the ean and covariance under different models at each stage.
expectation operator and covariance operator, respectively.
distribution of state vector ξ is Gaussian with mean 0 ξ and ing, initially, the target is in a CV motion. Hence, the ε1 and 1 sition probability is pij  To procedure for Φk- 1 and Qp,k-1 is given as follow: εk-1, k-1, and the mode probabilities pk-1(r), r = 1, 2, …, Nm babilities.
nd Nm is the number of the target motion models) specifies the target motion e time interval between sampling time indexes k-1 and k. Φk-1, Wk-1, and Qp,kx, process noise vector, and covariance matrix of the BFG approximation or convenient representation in the following process of derivation, the ction 2.2 (i = 1, 2, …, Nm) are changed to 1 k r k − F and p, 1 terval, the regime-dependent motion model (given by Equation (23) and replaced with a single BFG distribution (denoted by Equation (24) and . Φk- 1 and Qp,k-1 (Qp,k-1 ≥ 0) are calculated to satisfy the requirements that the e same mean and covariance under different models at each stage.
enote the expectation operator and covariance operator, respectively.
ocurement of Φk-1 is given in Appendix A. Then, with the aid of Φk-1, we d derivation in Appendix B.
Step 6. Update the mean of the state vector.

Prior Information J P (ξ k )
The prior information matrix J P (ξ k ) can be formulated by the equation below [20].
According the calculated Φ k−1 and Q p,k−1 by the linear BFG approximation, the D 11 k−1 , D 12 k−1 , and D 22 k−1 are deterministic and the expectation operator can be dropped out.

Data Information J D (ξ k )
On the other hand, the data information matrix J D (ξ k ) can be expressed by the equation below [20].
Since the measurements from different radars are independent from each other, J D (ξ k ) can be written as follows.
In practice, the expected value in Equation (39) may be evaluated using Monte Carlo techniques. To shorten the computation time, the predictive BIM is approximated as follows [34].
where ξ k|k−1 denotes the predicted target state for the case of zero process noise.

Predictive Bayesian Cramér-Rao lower bound (BCRLB-like)
The key feature of the resource allocation algorithm is the predictability. The predictive target tracking performance gives us the ability to make decisions in advance based on current knowledge. In the maneuvering target tracking, Φ k−1 and Q p,k−1 are updated through Equations (30) and (32) before each revisit time. They are both functions of the sampling interval T k . Hence, the mathematical symbol of the prior information matrix can be modified from J P (ξ k ) to J P (T k )| ξk (J P (T k ) has nothing to do with ξ k . Nevertheless, to keep the consistency of the format, it is denoted as J P (T k ) ξ k here). The data information matrix is affected by the transmitted power of the radars and the target RCS. The larger the transmitted power P k and the RCS h k are, the larger the data information matrix is. Consequently, the data information matrix can be rewritten as J D (P k ,h k )| ξk . Considering the updated BIM J(ξ k−1 ), sampling interval T k , transmitted power P k and RCS h k , we can calculate the predictive BIM.
where the meaning of ξ k|k−1 refers to Equation (40). The predictive target tracking performance can be calculated as the inverse of BIM. However, because of the BFG approximation technique involved, the predictive target tracking performance no longer provides a lower boundary on error performance. Hence, the predictive target tracking performance in this paper is named BCRLB-like instead of BCRLB.

of 25
The diagonal elements of BCRLB-like provide a referenced boundary on the variances of the estimation of the target's hybrid state. It is sufficient for us to utilize the tracking BCRLB-like factor of the illuminated target as a criterion for power allocation.
where Equation (43) denotes the integrated tracking performance of multiple radars at the kth sampling instant.
Having the tracking BCRLB-like F(T k ,P k ,h k ), we can allocate the power resource reasonably while the beams illuminate the target. However, due to the scarce radar resource, the radars do not have to illuminate the target continually. The beams only need to be transmitted when the tracking error is larger than a given tracking error threshold. As shown in Equations (22) and (36), the prior information matrix J P (T k ) is only related to a target motion equation. Therefore, the inverse of J P (T k ) is utilized to determine when to transmit the beams.
To accomplish an adaptive sampling interval, the prior CRLB-like is compared with the upper boundary of the given tracking error threshold. It can be used as a criterion to determine the optimal sampling period. The representation of the prior CRLB-like factor is given as follows.

Modeling of Chance-Constraint Programming (CCP)
In accordance with Equation (16), it can be seen that the target tracking accuracy is impacted by several parameters. The decision variables involved in this paper are the sampling interval T k and the transmit power P k . The random vector is the RCS h k . For the predetermined tracking error range at each sampling instant, the aim of our work is to minimize the transmitted power by optimally allocating the limited time and power resources.
where 1 T M = [1,1, . . . ,1] 1×M . P min is the lower boundary of the transmitted power of each beam. P total is the available total power. η 1,k and η 2,k are the lower boundary and upper boundary of the tracking error range, respectively. The adaptive sampling interval is included in the resource management model. If the prior CRLB-like F P (T k ) of the next sampling instant is greater than the given tracking error threshold η 2,k , T k is set as the sampling interval and the beams are transmitted to track the target.
In Equation (46), the model can make sure that the tasks can be accomplished by as many as possible conditions on the given tracking error threshold.
In practice, the target RCS is related to the identification, attitude, and position of the target, and it is affected by an aspect angle, wavelength, and polarization, etc. [25], i.e., the target RCS is unknown and uncertain. Consequently, the target RCS is considered as a random variable in this paper. Then the deterministic resource allocation model cannot show the characteristic of the target well. In view of the above situation, the stochastic CCP of the resource management is introduced [26]. According to Equation (46), the resource management model can be reformulated as: where Pr{•} is a probability measure operator. α is a confidence level.

Basic of the Technique
As a type of stochastic programming pioneered by Charnes and Cooper [35], the CCP offers a powerful means of modeling stochastic decision systems with the assumption that the stochastic constraints will hold at least α of time, where α is referred to as the confidence level provided as an appropriate safety margin by the decision-maker. It is convenient and general to deal with them by stochastic simulations. Hence, Liu [26] integrates the stochastic simulation and GA to produce a hybrid intelligent optimization algorithm (HIOA) for solving stochastic programming models.
In recent years, some new filtering algorithms have emerged one after another. The particle filter has the highest filtering accuracy [36,37]. However, since it is based on Monte Carlo methods, the calculation time of the particle filter is longer. The shadowing filter offers a robust methodology to position and track a moving target from limited positional information [38,39]. Nevertheless, the target information of this paper is not incomplete, or else the resource allocation process cannot be performed successfully. Under the condition that the filtering accuracy satisfies the calculation requirement of this paper, the unscented Kalman Filter (UKF) has a faster filtering speed [40]. In addition, it can also solve the filtering problems whose target information is complete. For the previously mentioned reasons, the UKF is selected.
In addition, since the dynamics of the linear jump Markovian system are expressed by the BFG approximation, the radars can track the maneuvering target by a UKF. On the other hand, we also employ the interacting multiple model UKF (IMM-UKF) to track the maneuvering target [40][41][42]. Thus, the tracking performance of BFG-UKF can be compared with the tracking performance of IMM-UKF to verify the efficiency of the BFG approximation.
The resource allocation processing procedure can be detailed as follows.

Stochastic Simulation
For the randomness of the target RCS, the stochastic simulation is used to calculate the stochastic CCP [26]. According to the expert experience and historical measurement data, it is assumed that there exists Nh q,k (Nh q,k is the number of measurements relative to the qth radar.) measurements before sampling, and each measurement is denoted as h q i,k (i = 1, 2, . . . , Nh q,k ). Let N H = Nh 1,k ×Nh 2,k × . . . ×Nh q,k . When the resource will be pre-allocated at the kth sampling time, we can select h j,k (j = 1, 2, . . . , N H ) from the measurement set of RCS to produce F(T k ,P k ,h j,k ). Let N denote the number of occasions on which F(T k ,P k ,h j,k ) ≤ η 1,k (i.e., the number of random vectors satisfy the system of inequalities). Let us define the following.
It follows from the strong law of large numbers that: This converges a.s. to Pr{F(T k ,P k ,h k )≤η 1,k }. Thus, the probability measure can be estimated by N /N H , which provided that N H is sufficiently large. The solving steps is illustrated in Table 1. Table 1. The stochastic simulation.
Step 1. Set N =0; Step 2. Calculate Φ k and Q p,k according to the BFG approximation; Step 3. Select a measurement h j,k (j=1, 2, . . . , N H ) from the measurement set and produce F(T k ,P k ,h j,k ); Step 4. If F(T k ,P k ,h j,k )≤η 1,k , N = N +1; Step 5. Repeat the third and fourth steps N H times; Step 6. Pr{F(T k ,P k ,h k )≤η 1,k }= N /N H .

Hybrid Intelligent Optimization Algorithm
In Equation (47), if F P (T k )≥η 2,k , we will solve the stochastic CCP model to obtain the optimal power allocation P k,opt . Hence, the optimal value T k,opt has been calculated before solving the stochastic CCP model. Then, we just need to solve the stochastic CCP model to obtain the optimal power allocation P k,opt . The stochastic simulation is introduced and embedded into GA to constitute HIOA for solving the stochastic CCP model [26]. The steps are described in Table 2. (1) Initialize pop_size chromosomes, and check the feasibility of the generated chromosomes by the stochastic simulation in Table 1; (2) Update the chromosomes by crossover and mutation operations in which the feasibility of offspring can be checked by the stochastic simulation in Table 1, and, if they do not satisfy the constraint, correct the chromosomes; (6) Repeat the second to fifth steps for a given number of cycles; (7) Report the best chromosome as the optimal solution P k,opt .

Target State Estimation
For a maneuvering target, the key to successful target tracking lies in the effective extraction of useful information about the target's state from observations. Furthermore, a good model of the target will facilitate this information extraction process to a great extent. In this paper, the dynamics of the linear jump Markov system can be replaced by the BFG approximation with a Gaussian distribution. The transition matrix and process noise of the BFG equation need to be updated before each sampling instant. Therefore, the maneuvering target state can be estimated by BFG-UKF (Herein, for the sake of clarity, the UKF including the BFG approximation is renamed as BFG-UKF). In addition, the IMM-UKF is a wildly accepted state estimation method for MTT [41,42]. Hence, we compare the tracking effects between BFG-UKF and IMM-UKF. The target state estimation processes are detailed in the following, respectively.
Step 2. BFG approximation: through the BFG approximation in Section 3.1, determine the mode probability p k (i), and then calculate Φ k−1 and Q p,k−1 (the detailed calculation process is given in Section 3.1).
Step 3. Prediction procedure: predict the target state and covariance.
To implement the sequential updating scheme of centralized tracking of multiple radars [27,43], the target state and covariance are rewritten as: Step 4. Calculation of sigma sampling points and their weights. For M radars, the measurement from the most accurate radar should be updated first so as to reduce subsequent linearization errors. Let q = 1, and we start the recursion from the predicted state and covariance denoted by the equations below.
) is the lth row or column of the matrix square root. N ξ is the dimensionality of the state vector. ω l,k is the weight that is associated with the lth point.
Step 5. Prediction procedure of observations. Transform each sigma point through the measurement equation, and compute the mean z q k|k−1 , covariance C q zz,k , and cross covariance C q ξz,k . Step 6. Target state updating: calculate the gain K q k , and update the state vector ξ q k|k and covariance C q k|k .
where z q,k is the measurement of the qth radar.
Step 7. Let q = q + 1, and if q ≤ M, go to Step 4. Otherwise, go to Step 8.
Step 9. Power optimal allocation: implement the power allocation algorithm in Table 2 and send the optimal allocation result to the multiple OAR system.

Process of Interacting Multiple Model Unscented Kalman Filter (IMM-UKF)
As we all know, IMM-UKF is an organic combination of IMM and UKF. Then, there exist many common steps between the processes of IMM-UKF and BFG-UKF. Hence, in view of the length of this paper and avoiding redundancy in content, we will simplify the steps that have appeared in the process of BFG-UKF.
Step 1. Initialization: let k = 1, and initialize T k,opt = ∆T (∆T denotes a short enough time slot during which the change of the target tracking error can be neglected), P k,opt = P 0 (P 0 denotes equal power allocation), the motion model probability p k−1 (i), the transition probability p ij , the state ξ k−1|k−1 (i), and the covariance C k−1|k−1 (i) = J −1 P (ξ k−1|k−1 (i)) for a mode-matched filter i, where i, j = 1,2, . . . , N m .
Step 2. Interactive input: calculate the mixing probability.
where p k|k−1 (j) is a predicted mode probability.
Thus, the mixing state estimation and the mixing covariance are computed interactively as: Step 3. Filtering of each motion model (j = 1, 2, . . . , N m ): for motion model j, the filtering process is the same as the process from Step 3 to Step 7 of BFG-UKF.
where the model likelihood function is: Step 5. Estimation fusion: fuse the states and covariances of all the mode-matched filters to compute the overall estimate and overall covariance.
Step 6. Adaptive sampling interval: this step is the same as Step 8 of BFG-UKF.
Step 7. Power optimal allocation: this step is the same as Step 9 of BFG-UKF.

Simulation Results and Analysis
To better illustrate the effectiveness of the proposed resource management scheme, the relevant numerical examples are given in this section. First, we demonstrate the advantage of the adaptive sampling algorithm by comparing it to the fixed sampling algorithm. Then, the total transmitted power is calculated to different confidence levels to account for the impacts of confidence levels on power consumption. Third, by comparing the results of the optimal power allocation and the equal power allocation, the optimal power allocation can not only save the total transmitted power strikingly, but also allocate the power among the distributed radars reasonably. Lastly, the tracking performance of BFG-UKF is compared with the tracking performance of IMM-UKF to demonstrate the validity of the BFG approximation. Now, we give the simulation setup first in the following.
Suppose that a radar network with M = 3 distributed radars is considered. The carrier frequency of each radar is set as 10 GHz, and, thus, the carrier wavelength is λ q = 0.03 m. The effective bandwidth and effective time duration of each radar beam are set as B q,k = 5 MHz and T q,k = 1 ms, respectively. The time slot of the tracking process is set to ∆T = 0.5 s. The number of the coherent pulse is 64. The lower bound of the transmitted power of the qth radar is P q,min = 0.01P total . The desired tracking error range is [400, 1000] m. The target is initially located at (40, 50) km. From 0 to 50 s, the target first flies with a constant speed of (−200, 200) m/s, and from 51 to 70 s. A CA motion is taken with an acceleration of 10 m/s 2 in the x and y coordinate directions, and from 71 to 85 s, it takes a CT motion with a known angular turn rate Ω= 0.15 rad/s. Lastly, it keeps moving in a CV motion mode from 86 to 135 s. The maneuvering target trajectory and the deployment of the radars in the tracking scenario is in Figure 3.
To compare the simulation results fairly and conveniently, it is assumed that all the simulation examples have the same target RCS. The total sampling times of the adaptive sampling algorithm cannot be predetermined in the tracking process. Then, we cannot also determine the number of the observed value of RCS. Herein, we give the changing curve of RCS according to the fixed sampling algorithm with the most sampling times. range is [400, 1000] m. The target is initially located at (40, 50) km. From 0 to 50 s, the target first flies with a constant speed of (−200, 200) m/s, and from 51 to 70 s. A CA motion is taken with an acceleration of 10 m/s 2 in the x and y coordinate directions, and from 71 to 85 s, it takes a CT motion with a known angular turn rate Ω = 0.15 rad/s. Lastly, it keeps moving in a CV motion mode from 86 to 135 s. The maneuvering target trajectory and the deployment of the radars in the tracking scenario is in Figure 3. To compare the simulation results fairly and conveniently, it is assumed that all the simulation examples have the same target RCS. The total sampling times of the adaptive sampling algorithm cannot be predetermined in the tracking process. Then, we cannot also determine the number of the observed value of RCS. Herein, we give the changing curve of RCS according to the fixed sampling algorithm with the most sampling times. Where |h i | is the modulus of target RCS relative to the ith radar in Figure 4. Where |h i | is the modulus of target RCS relative to the ith radar in Figure 4. is in Figure 3. To compare the simulation results fairly and conveniently, it is assumed that all the simulation examples have the same target RCS. The total sampling times of the adaptive sampling algorithm cannot be predetermined in the tracking process. Then, we cannot also determine the number of the observed value of RCS. Herein, we give the changing curve of RCS according to the fixed sampling algorithm with the most sampling times. Where |h i | is the modulus of target RCS relative to the ith radar in Figure 4. To proceed, the simulation results and corresponding analysis are presented in the following subsections. In the first three subsections, we all adopt the BFG-UKF for target state estimation. In the last subsection, we compare the tracking effects between BFG-UKF and IMM-UKF.

Adaptive Sampling Interval
This subsection supports the evaluation of the adaptive sampling algorithm. According to the previous analysis in Section 3.4, the adaptive sampling happens before the power allocation, and we predict the tracking error of the target by F P (T k ) to determine the sampling interval. The visual simulation results are shown in Figure 5.
The simulation results of the adaptive sampling algorithm and the fixed sampling algorithm are compared in Figure 5. In Figure 5a, besides the adaptive sampling algorithm, we also give two simulation results of the fixed sampling algorithm: T k = 6 s and T k = 6.5 s. When T k = 6.5 s, it can be seen from Figure 5 that the prior CRLB-like value does not satisfy the error threshold constraint in many moments. When T k = 6 s, the prior CRLB-like basically meets with the error threshold constraint during the whole tracking process. Hence, the fixed sampling interval T k = 6 s is compared with the adaptive sampling interval detailed in Figure 5b,c. During the whole tracking process, the adaptive sampling times 18 is significantly less than the fixed sampling times 22. In Figure 5b, the adaptive sampling interval is approximately 1.5 to 2 times of the fixed sampling interval T k = 6 s. The motion state is simple and the velocity is low from 0 to 70 s, so the target tracking error increases slowly, and then the radar system will select a long sampling interval intelligently. In Figure 5c, with the increase of maneuverability and velocity of the target, the length of the adaptive sampling interval is reduced to the length of the fixed sampling interval.
the last subsection, we compare the tracking effects between BFG-UKF and IMM-UKF.

Adaptive Sampling Interval
This subsection supports the evaluation of the adaptive sampling algorithm. According to the previous analysis in Section 3.4, the adaptive sampling happens before the power allocation, and we predict the tracking error of the target by FP(Tk) to determine the sampling interval. The visual simulation results are shown in Figure 5. The simulation results of the adaptive sampling algorithm and the fixed sampling algorithm are compared in Figure 5. In Figure 5a, besides the adaptive sampling algorithm, we also give two simulation results of the fixed sampling algorithm: Tk = 6 s and Tk = 6.5 s. When Tk = 6.5 s, it can be seen from Figure 5 that the prior CRLB-like value does not satisfy the error threshold constraint in many moments. When Tk = 6 s, the prior CRLB-like basically meets with the error threshold constraint during the whole tracking process. Hence, the fixed sampling interval Tk = 6 s is compared with the adaptive sampling interval detailed in Figures 5b and 5c. During the whole tracking process, the adaptive sampling times 18 is significantly less than the fixed sampling times 22. In Figure 5b, the adaptive sampling interval is approximately 1.5 to 2 times of the fixed sampling interval Tk = 6 s. The motion state is simple and the velocity is low from 0 to 70 s, so the target tracking error increases slowly, and then the radar system will select a long sampling interval intelligently. In Figure 5c, with the increase of maneuverability and velocity of the target, the length of the adaptive sampling interval is reduced to the length of the fixed sampling interval.

Optimal Allocation of Power
In this part, when the adaptive sampling algorithm is used for determining the sampling interval, we contrastively analyze the results of optimal power allocation and equal power allocation conditioned on the desired tracking error threshold. By restricting generality, the confidence level is set as α = 0.95 in this simulation.
The power consumption rate is given in Figure 6. The power consumption rate is defined as the ratio of the total transmitted power P sum and the total power P total of the radar system. Since the adaptive sampling algorithm is used for the two power allocation methods, it can be seen that they have equal sampling times. Under the condition of the same desired tracking error threshold, the optimal power allocation algorithm can save more power resources than the equal power allocation algorithm, especially between 0 s and 70 s. P sum and P total are the total transmitted power and the total power of the radar system, respectively. The power allocation algorithm proposed in this case can also intelligently allocate the power among the radars, and the corresponding simulation results are shown in Figure 7. With the moving of the target, the distance between the target and each radar is time-varying. Attributing to the spatially decentralized radars, the azimuth changing rate and target velocity relative to different radars are different. According to Equations (13) and (14), the measurement values of the range, azimuth, and the Doppler shift can all affect the power allocation, i.e., the power allocation is a comprehensive influence of the three elements. Therefore, the power allocation ratios are changed all the time when the target is moved. ratio of the total transmitted power Psum and the total power Ptotal of the radar system. Since the adaptive sampling algorithm is used for the two power allocation methods, it can be seen that they have equal sampling times. Under the condition of the same desired tracking error threshold, the optimal power allocation algorithm can save more power resources than the equal power allocation algorithm, especially between 0 s and 70 s. Psum and Ptotal are the total transmitted power and the total power of the radar system, respectively.
The power allocation algorithm proposed in this case can also intelligently allocate the power among the radars, and the corresponding simulation results are shown in Figure 7. With the moving of the target, the distance between the target and each radar is time-varying. Attributing to the spatially decentralized radars, the azimuth changing rate and target velocity relative to different radars are different. According to Equations (13) and (14), the measurement values of the range, azimuth, and the Doppler shift can all affect the power allocation, i.e., the power allocation is a comprehensive influence of the three elements. Therefore, the power allocation ratios are changed all the time when the target is moved.  Psum and Ptotal are the total transmitted power and the total power of the radar system, respectively.
The power allocation algorithm proposed in this case can also intelligently allocate the power among the radars, and the corresponding simulation results are shown in Figure 7. With the moving of the target, the distance between the target and each radar is time-varying. Attributing to the spatially decentralized radars, the azimuth changing rate and target velocity relative to different radars are different. According to Equations (13) and (14), the measurement values of the range, azimuth, and the Doppler shift can all affect the power allocation, i.e., the power allocation is a comprehensive influence of the three elements. Therefore, the power allocation ratios are changed all the time when the target is moved. P q,k denotes the transmitted power of the qth radar at the kth sampling instant. To examine the tracking performance of the proposed power allocation method, we compare the root MSEs (RMSE) of the state vectors of the two power allocation methods. In this case, the position RMSE is defined by the equation below.
where N MC is the number of Monte Carlo trials and it is set to 100 in this paper. (x i,k ,ŷ i,k ) is the state estimate of the target at the ith trial.
In Figure 8, although the optimal power allocation consumes less transmitted power, its tracking RMSE is not worse than the equal power allocation. This is because the intelligent increase/decrease of the transmitted power of the radars is certainly better than the simultaneous increase/decrease. In addition, the BCRLB-like no longer provides a lower bound on error performance. Meanwhile, to satisfy the specified confidence level of the CCP, the RCS used in the pre-allocation process is less than the RCS used in the state estimation process. Hence, there is no comparability between the RMSE in Figure 8 and the BCRLB-like in Figure 5a.
of the transmitted power of the radars is certainly better than the simultaneous increase/decrease. In addition, the BCRLB-like no longer provides a lower bound on error performance. Meanwhile, to satisfy the specified confidence level of the CCP, the RCS used in the pre-allocation process is less than the RCS used in the state estimation process. Hence, there is no comparability between the RMSE in Figure 8 and the BCRLB-like in Figure 5a. RMSE is the root mean square error of the target.

Chance-Constraint Programming
To cope with the randomness of RCS, the CCP is utilized in the simulation. Selecting an appropriate confidence level is crucial. Not only the identity, motion state, and number of the targets, but also the target environments and the amount of radar resource are the factors that impact the confidence level. Considering the previously mentioned factors comprehensively, we specify the confidence level α as 0.99, 0.95, and 0.9 in Figure 9, respectively. As seen from the simulation results, the power consumption is increased with the rise of the confidence levels.   RMSE is the root mean square error of the target.

Chance-Constraint Programming
To cope with the randomness of RCS, the CCP is utilized in the simulation. Selecting an appropriate confidence level is crucial. Not only the identity, motion state, and number of the targets, but also the target environments and the amount of radar resource are the factors that impact the confidence level. Considering the previously mentioned factors comprehensively, we specify the confidence level α as 0.99, 0.95, and 0.9 in Figure 9, respectively. As seen from the simulation results, the power consumption is increased with the rise of the confidence levels.
addition, the BCRLB-like no longer provides a lower bound on error performance. Meanwhile, to satisfy the specified confidence level of the CCP, the RCS used in the pre-allocation process is less than the RCS used in the state estimation process. Hence, there is no comparability between the RMSE in Figure 8 and the BCRLB-like in Figure 5a. RMSE is the root mean square error of the target.

Chance-Constraint Programming
To cope with the randomness of RCS, the CCP is utilized in the simulation. Selecting an appropriate confidence level is crucial. Not only the identity, motion state, and number of the targets, but also the target environments and the amount of radar resource are the factors that impact the confidence level. Considering the previously mentioned factors comprehensively, we specify the confidence level α as 0.99, 0.95, and 0.9 in Figure 9, respectively. As seen from the simulation results, the power consumption is increased with the rise of the confidence levels.

Target State Estimation with BFG-UKF and IMM-UKF
The regime-dependent motion model is replaced by a single BFG distribution. Consequently, we can use the BFG-UKF to estimate the target state in the previously mentioned three subsections. In this subsection, we compare the RMSE of BFG-UKF and IMM-UKF in Figures 10 and 11.
Where BFG-UKF denotes the UKF combined with BFG. For BFG-UKF (or IMM-UKF) alone in Figure 10, the RMSEs increase with the decrease of the confidence levels. With the decrease of the confidence levels, the power required to satisfy constraints decreases. Then the tracking performance of BFG-UKF (or IMM-UKF) also gets worse. From 0 s to 60 s and 115 s to 135 s, the tracking effects of BFG-UKF and IMM-UKF are nearly the same. However, between 60 s and 115 s, due to the maneuverability of the target, the estimation performance of the two filters both gets worse. By contrast, the estimation performance of BFG-UKF is better than the estimation performance of IMM-UKF. Now, we give an explanation. Although the traditional UKF is not suitable for tracking the maneuvering target, the BFG is introduced to improve the performance of it and a well simulation is obtained. In addition, in the simulation scenario of this paper, the spatial distribution of radars may result in the inexact calculation of the likelihood function, and then the estimation of the target motion mode is influenced in IMM-UKF. Hence, the simulation results in Figure 10 appear.
Meanwhile, Figure 11 verifies our viewpoint again. In addition, we give a local enlarged drawing to emphasize the simulation results. The confidence level is set as α = 0.95 in Figure 11.

Target State Estimation with BFG-UKF and IMM-UKF
The regime-dependent motion model is replaced by a single BFG distribution. Consequently, we can use the BFG-UKF to estimate the target state in the previously mentioned three subsections. In this subsection, we compare the RMSE of BFG-UKF and IMM-UKF in Figures 10 and 11. Where BFG-UKF denotes the UKF combined with BFG. For BFG-UKF (or IMM-UKF) alone in Figure 10, the RMSEs increase with the decrease of the confidence levels. With the decrease of the confidence levels, the power required to satisfy constraints decreases. Then the tracking performance of BFG-UKF (or IMM-UKF) also gets worse. From 0 s to 60 s and 115 s to 135 s, the tracking effects of BFG-UKF and IMM-UKF are nearly the same. However, between 60 s and 115 s, due to the maneuverability of the target, the estimation performance of the two filters both gets worse. By contrast, the estimation performance of BFG-UKF is better than the estimation performance of IMM-UKF. Now, we give an explanation. Although the traditional UKF is not suitable for tracking the maneuvering target, the BFG is introduced to improve the performance of it and a well simulation is obtained. In addition, in the simulation scenario of this paper, the spatial distribution of radars may result in the inexact calculation of the likelihood function, and then the estimation of the target motion mode is influenced in IMM-UKF. Hence, the simulation results in Figure 10 appear.
Meanwhile, Figure 11 verifies our viewpoint again. In addition, we give a local enlarged drawing to emphasize the simulation results. The confidence level is set as α = 0.95 in Figure 11.   Where BFG-UKF denotes the UKF combined with BFG. For BFG-UKF (or IMM-UKF) alone in Figure 10, the RMSEs increase with the decrease of the confidence levels. With the decrease of the confidence levels, the power required to satisfy constraints decreases. Then the tracking performance of BFG-UKF (or IMM-UKF) also gets worse. From 0 s to 60 s and 115 s to 135 s, the tracking effects of BFG-UKF and IMM-UKF are nearly the same. However, between 60 s and 115 s, due to the maneuverability of the target, the estimation performance of the two filters both gets worse. By contrast, the estimation performance of BFG-UKF is better than the estimation performance of IMM-UKF. Now, we give an explanation. Although the traditional UKF is not suitable for tracking the maneuvering target, the BFG is introduced to improve the performance of it and a well simulation is obtained. In addition, in the simulation scenario of this paper, the spatial distribution of radars may result in the inexact calculation of the likelihood function, and then the estimation of the target motion mode is influenced in IMM-UKF. Hence, the simulation results in Figure 10 appear.
Meanwhile, Figure 11 verifies our viewpoint again. In addition, we give a local enlarged drawing to emphasize the simulation results. The confidence level is set as α = 0.95 in Figure 11.  The joint resource allocation algorithm of this paper can be extendable to a three-dimensional scenario. In the three-dimensional scenario, the z-coordinate is added into target motion equations. In the measurement model, the elevation angle is added. Hence, compared with the two-dimensional scenario, the BCRLB of the elevation angle is added into Equation (16) in the three-dimensional scenario. The predictive BIM of the elevation angle needs to be calculated. Besides the elevation angle, the other predictive BIMs in Equation (40) also need to be derived again. The filter needs to be changed to a three-dimensional form. The joint resource allocation algorithm of the maneuvering target for a three-dimensional scenario is worthy of studying. The formula derivation and simulation are more complicated. For the length of this paper, the content of the three-dimensional scenario is not involved. We will study this part in the future work.

Conclusions
This paper presents a joint adaptive sampling interval and power allocation scheme based on CCP for MTT in a multiple OAR system. Initially, we develop a general approach to replace the MSD of a maneuvering target with a signal BFG distribution. Based on this, the prior CLRB-like is calculated as a criterion for determining the adaptive sampling interval. Meanwhile, the BCRLB-like, which provides a referenced boundary on the tracking error, is exploited as an aid in performing efficient power allocation for MTT. Afterwards, directing at the randomness of target RCS, we package the deterministic resource allocation model as an uncertain model through CCP conditioned on a specified confidence level. Lastly, the resulting optimization problem is solved through HIOA. Simulation results demonstrate the effect of the adaptive sampling algorithm and the power allocation algorithm, the correctness of HIOA for solving the CCP, and the accuracy of BFG-UKF for estimating the maneuvering target state. As can be seen from the derivation, the resource allocation scheme for a single maneuvering target can be easily generalized to a scenario including multiple maneuvering targets, even the distributed OAR case (consists of M 1 transmitting and N 1 receiving radars). In the future work, we will introduce the BFG approximation into the multiple maneuvering target tracking scenario of the distributed OAR.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
It is assumed that the target motion switches between N m models. The evolution of the motion model sequence is modeled by a time-homogenous Markov chain. Then, we define the motion model probabilities.
p k (r) Pr(r k = r) where is a defining operation symbol; r = 1, 2, . . . , N m . Using the total probability theorem, we get: If the initial model probabilities p 1 (r) and the transition probabilities p ij are known, we can determine the p 2 (r), p 3 (r), and so on.
Thus, the mean of ξ k under model 1 is shown below.
Substituting (25) into (A3), we can get the equation below. where