Bayesian Non-parametric Hidden Markov Model for Agile Radar Pulse Sequences Streaming Analysis

The graphical model of the agile BNP-HMM is shown in Fig. 3, wherein the intercepted pulses $\boldsymbol{P}=\{p_{t}\}_{t\geq 1}$, hidden state assignments $\boldsymbol{S}=\{s_{t}\}_{t\geq 1}$, sufficient statistics $\boldsymbol{\Theta}=({\boldsymbol{\varphi}^{k}})_{k=1}^{\infty}$, initial distribution $\boldsymbol{\pi}=(\pi_{k})_{k=1}^{\infty}$, transition matrix $\boldsymbol{A}=(\boldsymbol{a}_{j})_{j=1}^{\infty},\boldsymbol{a}_{j}=(a_{ji})_{i=1}^{\infty}$, concentration parameter $\alpha_{\pi},\alpha_{A}$, the agile hyper-parameter $\kappa$ and the length of the data batch $T$ are presented. The generative model can be formulated using the Bayesian theory as follows:

Based on the research presented in [36], to ensure the generative model is fully conjugate, the joint distribution of the sufficient statistics are generated using the Gaussian-Gamma distribution, which is given in (9).

wherein $\mu_{k},\sigma_{k}^{2}$ are the mean and variance of the Gaussian distribution, respectively; $\sigma_{k}^{-2}=\frac{1}{\sigma_{k}^{2}}$ represents the precision of the Gaussian distribution; $\xi_{0},\lambda_{0},a_{0},b_{0}$ are hyper-parameters of the Gaussian-Gamma distribution; $N$ refers to Gaussian distribution and $Gamma$ refers to Gamma distribution³³3In MFR applications, each work mode is modeled by a BNP-HMM. The hidden state number in BNP-HMM generally varies across different MFR work modes. The work mode can change over time, and the number of hidden states would also change over time. The infinite assumption in Bayesian-nonparametric theory does not require the number of hidden states be pre-defined during inference.

The Dirichlet distribution to the HMM encourages hidden states to have similar transition distributions. However, it does not differentiate self-transitions from moves between the hidden states. Therefore, when modeling the MFR work mode with the agility characteristic, the self-transition nature of the Dirichlet distribution allows for a high posterior probability of hidden state persistence. Thus it tends to underestimate the hidden state number.

This paper proposed the agile BNP-HMM to determine the number of hidden states. The agile BNP-HMM is governed by placing an agile prior distribution over the state transition probability. The agile prior distribution is used to present the case of a few self-transitions. To implement the agile prior distribution, a novel stick-breaking construction method is designed.

At first, the Dirichlet distribution is reviewed. For any partitions $G=\{G_{1},...,G_{\infty}\}$ of the probability space $\Omega$, the distribution of the base measure’s probability mass on the partition $G$ is defined as follows:

wherein $\mathcal{DIR}$ denotes the Dirichlet process; $H$ is the base measure; $\alpha_{A}$ represents the concentration parameter; $\boldsymbol{a}_{j}$ is a random variable that follows the multi-nominal distribution, and the $\mathcal{DIR}$ can be implemented as a stick-breaking construction⁴⁴4The construction can be briefly explained. Suppose there is a stick with length 1. Let $a_{ji}^{\prime}\sim Beta(1,\alpha_{A})$ for $i=1,2,3,...,\infty$, and regard them as fractions for how much we take away from the remainder of the stick every time. Then $a_{ji}$ can be calculated by the length we take away each time: $a_{j1}=a_{j1}^{\prime},a_{j2}=(1-a_{j1}^{\prime}),...,a_{ji}=a_{ji}^{\prime}\prod_{n=1}^{i-1}(1-a_{jn}^{\prime})$. By the stick-breaking construction, we will have $\sum_{i=1}^{\infty}a_{ji}=1$. [30]:

wherein $Beta$ refers to the beta distribution and $\alpha_{A}$ is the concentration parameter.

An agile hyper-parameter is introduced to the Dirichlet distribution to control the hidden states self-transitioning tendency of MFR pulse sequences, resulting in agile Dirichlet process (agile $\mathcal{DIR}$), which is given by:

wherein $\kappa$ is the agile hyper-parameter, and a larger $\kappa$ value encourages an HMM to move from the current state $j$ to another state instead of staying in the current state.

The agile $\mathcal{DIR}$ can be implemented via a new stick-breaking construction as follows:

wherein $\sum_{i}^{\infty}a_{ji}=1$, and $\delta$ is the Dirichlet-delta function. When $i=j$ , $a_{jj}^{\prime}$ is sampled from $Beta(1,\alpha_{A}+\kappa)$, the $Beta$ distribution tends to put more mass on smaller values. Thus $a_{jj}^{\prime}$ tends to smaller than $a_{ji}^{\prime}$, leading to break a shorter stick.

In the inference procedure, it is necessary to estimate the posterior distribution from the observed data and the given hyper-parameters $\boldsymbol{\Upsilon}=\{\alpha_{\pi},\alpha_{A},a_{0},b_{0},\lambda_{0},\xi_{0},\kappa\}$. The inference of agile BNP-HMM is illustrated as Algorithm 1.

Variational inference is an efficient way to approximate the intractable posterior, which introduces a variational distribution q to approximate the true posterior. The basic idea of variational inference is to minimize the distance between the variational distribution q and the exact posterior p. The likelihood function is formulated as (14):

wherein $\mathcal{L}$ denotes the well-known evidence lower bound (ELBO), and $\mathcal{KL}$ is the Kullback-Leibler divergence. Minimizing the similarity between the variational distribution and the exact posterior is to maximize the evidence lower bound [36].

This paper considers the partial mean-field variational family [37] and truncation assumption for variational inference. The partial mean-field family assumes that $(\boldsymbol{\pi,A,\Theta,\Upsilon})$ and $(\boldsymbol{S})$ are mutually independent. The truncation assumption considers the probabilities of $L$ states as the infinite hidden states, wherein $L$ is called the truncation level, and $L$ should be large enough to ensure the accuracy. The variational distribution is defined as follows:

Exponential family is chosen for each component of the variational distribution. The variational distributions can then be represented as follows:

Then the evidence lower bound can be derived according to:

This paper uses the coordinate ascent algorithm [38] to maximize the ELBO. The universal update function is given by [36]:

wherein $\boldsymbol{Z}=\{\boldsymbol{\pi,A,\Theta,S}\}$, $n$ refers to the $n$th hidden variable. The details of the derivation and the update function can be found in the supplementary material.

After setting the hyper-parameters $\Upsilon$, pulse PRI values were collected in a buffer called the initial-batch. Since this paper assumes no knowledge about the data in the initial-batch, the Dirichlet Process Mixture Model (DPMM) [39] is used to obtain an initial guess of agile BNP-HMM model parameters of the data in initial-batch. For subsequent PE task, due to the streaming nature, six hyper-parameters that need to be tuned $\boldsymbol{\Upsilon}=\{\alpha_{A},\alpha_{\pi},\xi_{0},\lambda_{0},a_{0},b_{0}\}$, wherein $\boldsymbol{\Upsilon}$ is time variant and thus not suitable for change point detection tasks. In this work, the agile hyper-parameter $\kappa$ represents the time-series characteristics over time, and $\kappa$ is the only parameter need to be tuned for the PE task. The whole framework is illustrated as Algorithm 2.

Since our model is fully conjugate, we only need to initialize the sufficient statistics of (16)-(19). Specifically, $\eta_{1}(\pi_{1}^{\prime}),\eta_{2}(\pi_{2}^{\prime}),\eta_{1}(a_{ji}^{\prime}),\eta_{2}(a_{ji}^{\prime}),a_{j},b_{j}$ are all set to one in the initialization step. The mean $\mu_{j}$ and variance $\sigma_{j}^{-2}/\lambda_{j}$ of each Gaussian distributions are set to 0 and 0.1 in the initialization step.

After initialization, each pulse in the initial-batch belongs to a particular hidden state. In parameter estimation step, a new pulse is received and combined with previously arrived pulses to create a new data batch. Given the new data batch, the agile BNP-HMM algorithm is used to update the transition matrix and Gaussian distribution parameters. Due to the non-ideal observation conditions, there will be $\Delta TOA$ values caused by missing and spurious pulses. These $\Delta TOA$ values will result in a false hidden state with a relatively low visiting probability and a high transition probability. An acceptance threshold $\gamma_{u}$ is set, when $\mathbb{E}[\boldsymbol{a}_{j}]\leq\gamma_{u}$ satisfies, the corresponding hidden states are considered as $\Delta TOA$ values caused by missing or spurious pulses and are deleted. The remaining hidden states are assumed originating from true PRI values.

Streaming Bayesian updating enables high efficient parameter estimation and reduces the risk of falling into local optima during the PE task. Assume that $\boldsymbol{P}^{n}=\{p_{t}\}_{1}^{m+n}$, and index $m$ and $n$ as the initial-batch length and the $n$th arrived pulse after the initial-batch, respectively. For the initial-batch, applying the Bayesian rule yields $p(\boldsymbol{Z}^{0}|\boldsymbol{P}^{0})=p(\boldsymbol{P}^{0}|\boldsymbol{Z}^{0})/p(\boldsymbol{P}^{0})$. Superscript $0$ denotes batch index $0$. When a new pulse arrives, the batch index increases by one and a random hidden variable is assigned to the new arrived pulse. All underlying state sequences are combined and denoted as $p(\boldsymbol{Z}^{1})$. By applying the Bayesian rule again, instead of placing the prior $p(\boldsymbol{Z}^{1})$, the previous data batch’s variational posterior $q(\boldsymbol{Z}^{0}|\boldsymbol{P}^{0})$ is set as the current data batch’s prior. In variational inference, the true posterior for the $n$th batch $p(\boldsymbol{Z}^{n}|\boldsymbol{P}^{1:n})$ is approximated by the variational distribution $q^{n}(\boldsymbol{Z}^{n})$. Thus, the online variational inference is obtained by substituting $q^{n-1}(\boldsymbol{Z}^{n-1})$ for the prior as follows:

The estimated parameter tuple $(\boldsymbol{\pi},\boldsymbol{A},\boldsymbol{\Theta})^{n}$ for the $n$th batch is sent to a stack. All parameter tuples estimated in the previous time step are pushed in the stack. Then, the newly estimated tuple and the stack are sent to the change point detector.

Detecting the change point of the MFR work mode is equivalent to detecting the change point for BNP-HMM. However, BNP-HMM’s non-ergodicity and infinite hidden state space render general results [28, 27] invalid. Additionally, computing the invariant measure of HMM is computationally demanding. Directly detecting changes from HMM contradicts the basic principle of designing a system in MFR applications - that it should be computationally efficient while clear enough to distinguish various work modes. Therefore, we treat the prior distribution and the transition distributions as nuisance parameters similar to [29], which makes previous research results [32] valid. Our approach resembles Stephen Boyd’s method [33] by transforming learning complex structure problem into multivariate Gaussian inverse covariance matrices to simplify computationally intensive problems.

The change point detector considers the weighted SPRT [32]. Taking modulation parameter level change point as an example, the agile BNP-HMM sequentially returns the parameter set $\boldsymbol{\Theta}_{0},\boldsymbol{\Theta}_{1},\boldsymbol{\Theta}_{2},...,\boldsymbol{\Theta}_{L-m}$. The parameter sets can be described as a multi-variate Gaussian distribution as follows:

wherein $\boldsymbol{\theta}_{0}=(\boldsymbol{\mu}_{0}^{1},...,\boldsymbol{\mu}_{0}^{K})$ and $\boldsymbol{\theta}_{1}=(\boldsymbol{\mu}_{1}^{1},...,\boldsymbol{\mu}_{1}^{K})$ are the pre- and post-change parameters, respectively. $\boldsymbol{\Sigma}$ is the unit covariance matrix.

Following common assumptions in asymptotic performance analysis of the change point detector, we assume that the pre-change parameter $\boldsymbol{\theta}_{0}$ is known, and the post-change parameter $\boldsymbol{\theta}_{1}$ is unknown. The change amplitude $b$ is known and given by:

The stopping time $N$ is expressed by:

wherein $\bar{S}_{k}^{t}=\sum_{i=k}^{t}\left(\boldsymbol{\Theta}_{i}-\boldsymbol{\theta}_{0}\right)$,$G(d,z)=1+\frac{z}{d}+\cdots+\frac{z^{n}}{d(d+1)\ldots(d+n-1)n!}+\cdots$ is generalized hyper-geometric function, wherein $n$ is the index of the $n$ th item of the expansion. The strategy is asymptotic optimal in the sense of mini-max criterion [32]. The result is stand by the following theorem:

The proof of Theorem 1 can be found in supplementary material⁵⁵5Currently, we do not consider the situation of changes only in transition matrix between adjacent work modes. Because this situation is relatively rare in MFR applications. But as stated in the Introduction Section, the change point detection of this situation is non-trivial. We leave this for future study.. We also show that the result of the comparison between the simulation results (green dots) and the asymptotic case (red dashed line) in Fig. 4, wherein each green dot represents the average value of 100 Monte Carlo simulations. The linear relationship between the MDD and the log of MT2FA in asymptotic cases is verified. The simulation results approach the asymptotic case when $\ln{\bar{T}}$ is large (such as $\ln{\bar{T}}\geq 3\times 10^{3}$). When $\ln\bar{T}\to\infty$ is not satisfied, the Berk’s theorem is not satisfied (see supplementary material), and the linear relationship between the $\ln\bar{T}$ and $\bar{\tau}^{*}$ can not be guaranteed. Meanwhile, a smaller value of $\ln\bar{T}$ implies that it is easier for the detector to trigger an alarm. Consequently, a smaller value of $\ln\bar{T}$ leads to a smaller value of MDD. This explanation also clarifies why the MDD of the simulation result (green dots) is lower than the asymptotic case (red dashed line) when the value of $\ln\bar{T}\to\infty$ is small.

Simulation datasets included PRI sequences with multiple segments generated by the G-HMM described in Section III. Given the modulation type, the transition matrix of the single Work Mode (WM) is randomly sampled, and corresponding parameters are sampled based on the G-HMM sufficient statistics. The modulation types and parameters of each simulation category are given in Table II. Square brackets contain modulation parameters. The three elements in the parentheses of sliding PRI and D&S PRI indicate the minimum and maximum values for sliding PRI, as well as the number of PRI levels within a period. The number of one PRI level of D&S PRI is randomly sampled from 5 to 10. Jittered PRI is represented by parentheses indicating its mean value and variance. Considering the signal propagation error and the measurement error, additive Gaussian measuring noise with zero mean and variance 1 is added.

The first-category of the simulation explored the PE performance of agile BNP-HMM for different agile hyper-parameter $\kappa$ values. Each pulse sequence contains pulses from one single work mode with 1000 pulses. The comparison of results and the discussions are given in Section IV-B.

The second and the third categories of simulation examine the CPD performance, including the detection accuracy and timing performance. Each pulse sequence sample contains four work modes, each of which lasts 150 pulses. The work mode change in the Same modulation Type and Variable modulation Parameters (STVP) (D2–D5), and Variable modulation Types and Variable modulation Parameters (VTVP) (D6) are simulated in the second and the third categories, respectively. The results of the STVP, VTVP and different parameter settings are presented in Section IV-C.

Non-ideal conditions are added in each datasets. A length $N$ pulse sequence with $n$% non-ideal conditions is consisted with $N\times n/2\%$ number of randomly generated spurious pulses and $N\times n/2\%$ randomly generated missing pulses, respectively.

The annotation error and Hamming distance are used to evaluate the parameter estimation performance.

Annotation Error $\Delta K$: The annotation error reflects the ability of the algorithm to estimate the hidden states number:

wherein $\hat{K}$ and $K$ are the estimated and the true number of hidden states, respectively. “$\Delta K$ close to zero” implied the algorithm achieved accurate estimations of hidden state numbers.

Hamming Distance: This metric indicates the ability to estimate an underlying state sequence, and the smaller its value is, the higher the estimation accuracy is. The Munkres algorithm [41] is used to map randomly selected indices of the estimated state sequence to the set of indices that maximized the overlap with the true sequence.

In the change point detection task, the performance is evaluated by following five metrics: F1 score, MDD, MT2FA, FAR, and MR.

F1 Score: The F1 score is given in [42], which reflects the detection accuracy.

Mean Detection Delay (MDD) and Mean Time to False Alarm (MT2FA): The MDD and MT2FA are introduced by Lordon[35], and they are respectively expressed by (5) and (6).

False Alarm Rate (FAR) and Missing Rate (MR): False alarm rate and missing rate are two typical metrics that represent the probability of incorrect detection [43].

Four change point detection methods are used as baseline methods in this study:

1. 1.

   Agile BNP-HMM-FSS is implemented with the combination of the proposed parameter estimation method (agile BNP-HMM) and the sequential change point strategy of $\chi^{2}$-Fixed-Size-Sample ($\chi^{2}$-FSS)[32]. The stopping time of $\chi^{2}$-FSS is calculated by:

   $$N=\inf\{mj:d_{j}=1\};\ d_{j}=\begin{cases}1,\ $if$~{}|\chi_{(j-1)m+1}^{jm}|\geq h\\ 0,\ $if$~{}|\chi_{(j-1)m+1}^{jm}|<h\end{cases}$$

   (35)

   wherein $\chi_{(j-1)m+1}^{jm}$ is given by (27), $m$ is the fixed-size and $h$ is the threshold.

2. 2.

   U-FSS[44] has been used in our previous work on the MFR work mode online change point detection. It combines the generalized likelihood ratio test for the PE task and the original FSS algorithm[45] for the CPD task.

3. 3.

   U-CUSUM is proposed in [44] to fulfill the requirements of the MFR work mode online change point detection task. It combines the generalized likelihood ratio test for the PE task and the original CUSUM algorithm[45] for the CPD task.

4. 4.

   ChangeFinder is a unifying framework for detecting outliers and change points from time series [46]. It combines the auto-regressive model for the PE task and the logarithmic loss for the CPD task. It is used for the comparison under the non-ideal observations in this paper.

For the simplicity of representation, in the following, ABHC is short for Agile BNP-HMM-CUSUM, ABHF is short for Agile BNP-HMM-FSS, and CF is short for ChangeFinder.

The first simulation examines the influence of parameter $\kappa$ settings for both agile and non-agile inter-pulse modulation types. Agile BNP-HMM with $\kappa=0,0.5,1.0$ is used for five modulation types on D1. The estimated results of the hidden state number under different $\kappa$ settings are presented in Fig. 6. In the simulations, the BNP-HMM with hyper-parameter $\kappa=0.5,1$ outperform that with $\kappa=0$ on the agile inter-pulse modulation types. Noted that this study focuses on the parameter estimation for agile inter-pulse modulation types. In subsequently reported simulations, the agile BNP-HMM with $\kappa=0.5$ is used on agile inter-pulse modulation types, and that with $\kappa=0$ was used for comparison.

The performance under the non-ideal observations is presented in Fig. 7 (a) and 7 (b). The agile BNP-HMM performs better with $\kappa=0.5$ than $\kappa=0$. The $\Delta K$ is close to zero when $\kappa=0.5$ under various non-ideal settings. Note that in Fig. 7 (a), the agile BNP-HMM with $\kappa=0$ tends to underestimate the number of hidden states (the $y$-axis denotes $-\Delta K$), meanwhile the presence of more non-ideal pulses provides a more accurate hidden states assignments (the $\Delta K$ is close to zero when non-ideal ratio increase), which is the advantage of this method.

To analyze the reason of the phenomena, the PE task results of the non-ideal ratio values of 1% and 10% are shown in Fig. 8. The hidden state labels inferred by the agile BNP-HMM are indicated by stars of different colors. The black stars represent $\Delta TOA$ values caused by missing and spurious pulses. The blue lines indicate the temporal relationship between different real PRI values and $\Delta TOA$ values. Theoretically, the original Dirichlet prior $(\kappa=0)$ encourages hidden states to have similar transition distributions $\mathbb{E}[a_{jk}|\alpha_{A}]$. A relatively high self-transition probability indicates that each PRI value belongs to a mutual hidden state, as shown in the left part of Fig. 8. However, as shown in the right part of Fig. 8, the difference between the $\Delta TOA$ values and real PRI values are large, thus the $\Delta TOA$ values could be easily identified as a new external hidden state with extremely low self-transition probability. Since outliers are randomly encountered and would violate the similar state transition assumption. When the detector encounters a $\Delta TOA$ value, the algorithm are more likely to move to another hidden state or create a new hidden state. By accumulating the transition tendency, more hidden states will be created. In the middle part, the detail of estimated hidden states from a pulse segment under 1% and 10% non-ideal ratio are shown. Only one hidden state is inferred under 1% non-ideal ratio, and four hidden states are inferred under 10% non-ideal ratio.

The second category of simulation is performed to evaluate the results of temporal relationship prediction. The normalized Hamming distance is calculated for the three typical agile inter-pulse modulation types at different non-ideal ratios. As shown in Fig. 9 (a), with the increase in the non-ideal observation ratio, the agile BNP-HMM achieves better performance. The explanation is the same for the results in Fig. 8.

From Fig. 9 (a), the proposed model at $\kappa=0.5$ has a smaller normalized Hamming distance than the model at $\kappa=0$. There are two conclusions that can be drawn based on the results presented in Fig. 9 (b). Qualitatively, the agile BNP-HMM at $\kappa=0.5$ is robust to the non-ideal situation with the increase in the non-ideal ratio, and the normalized Hamming distance remains close to zero. The variance of agile PRI is larger than the other two modulation types. The hidden state transition distribution of agile PRI follows a uniform distribution, resulting in some miss-assignments. Quantitatively, the values of normalized Hamming distance are between 0.005 and 0.01, implying that out of 1,000 pulses in the dataset. Five to ten pulses are assigned to wrong hidden states. The normalized Hamming distance is getting lower with the increase in the non-ideal ratio.

Generally, we recommend setting the $\kappa$ to 0.5 when estimating the parameter of agile MFRs, and setting $\kappa$ to 0 when estimating the parameter of non-agile MFRs. In practical systems, we can not always obtain the ground truth to calculate $\Delta K$ and Hamming distance, thus we propose to perform the PE task on two branches with $\kappa$ values of 0.5 and 0. By calculating these branches in parallel, we can select the estimated value from the branch with the higher likelihood.

This section evaluates the performance of the proposed method and baselines under the optimal parameter settings. The parameters of the change detector are optimized based on the F1 score. As shown in Table III, the proposed ABHC method generally outperforms other baseline methods.

The results validate the robustness of the proposed framework. As discussed in Section IV-B, the agile BNP-HMM with $\kappa=0.5$ achieves slightly inferior results in estimating the true number of hidden states and hamming distance for non-agile inter-pulse modulation types (such as Jittered PRI). In the change point detection, either the ABHC or the ABHF has an superior performance. Although the agile BNP-HMM may not always accurately identify all hidden states, it can successfully abstract the time series dynamics of different radar work modes, thus the performance of the CPD task is not degraded facing the non-agile work mode.

In practical applications, we have considered two scenarios. In the first scenario, there are some pre-accumulated pulse sequences with change point information analyzed and recorded by experts. Due to the existence of electronic intelligence these datasets are generally available. In this scenario, the “optimal” parameters can be optimized by some pre-defined objectives, such as the F1 score. These parameters are then used for subsequent pulse sequences. In the second scenario, there is no available dataset like in the first scenario. For instance, the MFR is an unknown MFR or a known MFR with unknown or software defined signals. In this case the corresponding radar signals may not be recorded by interceptor. In this scenario, parameters are directly defined and used for all change point detection tasks.

In this section, a detection scatter plot, as shown in Fig. 10, is designed to present the trade-off between MDD and MT2FA by tuning the change detector parameters. The $x$-axis represents the MT2FA, and the $y$-axis denotes the MDD. The feature of a scatter can be represented by a tuple $(\mathcal{T},\mathcal{S})$, wherein $\mathcal{T}$ and $\mathcal{S}$ represent the transparency and the size of a scatter, respectively. Since U-FSS and U-CUSUM is lack of interpretation, they are not simulated in the subsequent simulations. Fig. 10 illustrates the influence of the thresholds selection on the quality metrics of the change point detection approaches. The parameters and their mapping on tuples were defined:

wherein $(r,k)$ represents the decay and lagging factors of the CF; $(m,h)$ is the fixed-size and the threshold of the ABHF; and $\gamma_{p}$ represents the threshold of the ABHC. We first define the parameter range $(\mathcal{T}_{min},\mathcal{T}_{max}),(\mathcal{S}_{min},\mathcal{S}_{max})$. The candidate parameters are generated using equidistant sampling.

From the perspective of the scatter location, the major scatters of CF (green) are located at the bottom left of every sub-figure of Fig. 10, which implies a high false alarm rate. Meanwhile, most scatters of the ABHC framework (red) are located at the bottom right of the each sub-figure of Fig. 10, which demonstrates better timing performances. As for the ABHF framework, the specific design of the parameter estimation task ensures better performance at MT2FA. In some cases (see Fig. 10 (a), 10 (b), and 10 (d)), the MT2FA value reaches 150, which indicated the zero false alarm rate on these data samples. The FSS strategy uses the data window for data points accumulation and decides whether there is a change point at the end of the data window (see (35)), resulting in longer detection delays. The CUSUM strategy makes decisions for each incoming data point (see (25)), resulting in smaller detection delay.

From the perspective of parameter settings, the ABHC framework is more suitable for practical applications than the baseline methods. As shown in Fig. 10 (e), the CF reconstructs the data via Auto-Regressive (AR) model, which is commonly used for modeling the stationary time series, causing the un-interpretable pattern for agile inter-pulse modulation. For ABHF, the scatters can be divided into five groups in terms of MDD as five $m$ were set in the simulations, reflecting the trade-off between the MDD and MT2FA. The transparency of each group from the bottom to the top decreased. The higher the fixed-size $m$ was set, the higher MT2FA and MDD were, achieving a trade-off between them. The size of each scatter increases in terms of the MT2FA value, indicating that a higher threshold $h$ would result in a higher MT2FA value. For ABHC framework, there is only one threshold parameter to be tuned, and the MT2FA increases (the false alarm rate decreased) with the scatter size increases (i.e., by setting a higher threshold). Varying the parameters would result in a stable MDD at a relatively low level, and the concentration of the ABHC scatters are much higher than in other methods.