Multi-task hidden Markov modeling of spectrogram feature from radar high-resolution range profiles

Abstract

In radar high-resolution range profile (HRRP)-based statistical target recognition, one of the most challenging task is the feature extraction. This article utilizes spectrogram feature of HRRP data for improving the recognition performance, of which the spectrogram is a two-dimensional feature providing the variation of frequency domain feature with time domain feature. And then, a new radar HRRP target recognition method is presented via a truncated stick-breaking hidden Markov model (TSB-HMM). Moreover, multi-task learning (MTL) is employed, from which a full posterior distribution on the numbers of states associated with the targets can be inferred and the target-dependent states information are shared among multiple target-aspect frames of each target. The framework of TSB-HMM allows efficient variational Bayesian inference, of interest for large-scale problem. Experimental results for measured data show that the spectrogram feature has significant advantages over the time domain sample in both the recognition and rejection performance, and MTL provides a better recognition performance.

1. Introduction

A high-resolution range profile (HRRP) is the amplitude of coherent summations of the complex time return from target scatterers in each range resolution cell, which represents the projection of the complex returned from the target scatting centers onto the line-of-sight (LOS), as shown in Figure 1. Since it contains the target structure signatures, such as target size and scatterer distribution, radar HRRP target recognition has received intensive attention from the radar automatic target recognition (RATR) community [1–16].

Figure 1

Illustration of an HRRP sample from a plane target, where the scatters on the plane target are represented in circles. This figure is cited from [1].

Several studies [8–16] show that statistical recognition is an efficient method for RATR. Figure 2 shows a typical flow chart of radar HRRP statistical recognition. By statistical recognition is meant the feature vector y extracted from test HRRP sample x will be assigned to the class with maximum posterior probability p(c|y), where c ∈ {1, ..., C} denotes the class membership. According to Bayes algorithm, p(c|y) ∝ p(y|c)p(c), where p(y|c) is the class-conditional likelihood and p(c) is the prior class probability. Since the prior class probability is usually assumed to be uniformly distributed, estimation of the posterior probability p(c|y) of each class is turned into estimation of the class-conditional likelihood p(y|c) of each class. There are usually two stages (training and classification) in the statistical recognition procedure. We suppose the class-conditional likelihood p(y|c) can be described via a model with a set of parameters (i.e., a parametric model). In the training phase, these parameters are estimated via training data (known as statistical modeling); and in the classification phase, as discussed above, given a test sample x, we first extract the feature vector y from the test sample, then we calculate the class-conditional likelihood p(y|c) for each target c, finally, the test sample is associated with target c' if $c^{\prime }=arg\underset{c}{max}p\left(c|\mathbf{y}\right)$. The focus of this article is on the feature extraction and statistical modeling.

Figure 2

Typical flow chart of radar HRRP statistical recognition.

Feature extraction from HRRP is a key step of our recognition system. One of the general feature extraction methods is feature dimensionality reduction [5]. This method is generally supervised, and some discriminative information may be lost during the dimensional reduction procedure. Another general feature extraction method is feature transformation. The study given in [6] investigates various higher-order spectra, and shows that the power spectrum has the best recognition performance. Rather than just utilizing the frequency domain feature as in [6], this study exploits the spectrogram feature of HRRP data for combining both time domain and frequency domain features, which is a two-dimensional feature providing the variation of frequency domain feature with time domain feature. Some statistical models [8–16] are developed for HRRP-based RATR, of which [14–16] successfully utilized the hidden Markov model (HMM) for modeling the feature vectors from the HRRP sequence. Since the HRRP sample is a typical high-dimensional distributed signal, it is computationally prohibitive to build such a high-dimensional HMM for describing HRRP sequences directly. Therefore, to avoid this problem, some dimensionality reduction methods are utilized. For example, in [14, 15], the relax algorithm is employed to extract the waveform constituents from the HRRP radar echoes, and then an HMM is utilized to characterize the features; in [16], a nonstationary HMM is utilized to characterize the following two features, i.e., the location information feature of scattering centers which are extracted via the multirelax algorithm, and the moments of HRRP radar echoes. Nevertheless, some information contained in HRRP samples will be inevitably lost during the dimensional reduction procedure. Moreover, since multiple aspect-dependent looks at a single target are utilized in the classification phase, the angular velocity of the target relative to the LOS of the radar is required to remain the same for the training and classification phases, which can hardly be satisfied for a noncooperative target.

The study reported here seeks an alternative way of exploiting HMM, in which we characterize the spectrogram feature from a single HRRP sample via the hidden Markov structure. In our model, the spectrogram feature extracted from each single HRRP is viewed as a d-dimensional sequence (d is the length of spectrogram feature in the frequency dimension), thus only a single HRRP sample is required for the classification phase rather than an aspect-dependent sequence. The main contribution of this study can be summarized as follows.

• (a) Spectrogram feature: The time domain HRRP samples only characterize the time domain feature of the target, which is too simple to obtain good performance. By contrast, the spectrogram feature introduced in this article is a time-frequency representation of HRRP data. The physical meaning of the spectrogram feature extracted from an HRRP sample is that the spectrogram feature in each time bin characterizes the frequency domain property of a fragment of the target, which can reflect the scattering properties of different physical constructions. Therefore, the spectrogram feature should be a better choice for the recognition problem.

• (b) Nonparametric model selection via stick-breaking construction: In the context of target recognition using HMMs, a key issue is to develop a methodology for defining an appropriate set of states to avoid over- or under-fitting. A Bayesian nonparametric infinite HMM (iHMM) which constituted by the hierarchical Dirichlet process (HDP) has proven effective to infer the number of states in acoustic sensing scenarios [17, 18]. However, the lack of conjugacy between the two-level HDP means that a truly variational Bayesian (VB) solution is difficult for HDP-HMM, which makes computationally prohibitive in large data problems. Recent study [19] proposes another way to constitute iHMM, where each row of transition matrix and initial state probability is given with a fully conjugate infinite dimensional stick-breaking prior, which can accommodate an infinite number of states, with the statistical property that only a subset of these states will be used with substantial probabilities, referred to as the stick-breaking HMM (SB-HMM). We utilize the truncated version of such stick-breaking construction in our model to characterize the HMM states, which is referred to as the truncated stick-breaking HMM (TSB-HMM).

• (c) Multi-task learning (MTL): A limitation of statistical model is that it usually requires substantial training data, assumed to be similar to the data on which the model is tested. However, in radar target recognition problems, one may have limited training data. In addition, the test data may be obtained under different motion circumstance. Rather than building models for each data subset associated with different target-aspect individually, it is desirable to appropriately share the information among these related data, thus offering the potential to improve overall recognition performance. If the modeling of one data subset is termed one learning task, the learning of models for all tasks jointly is referred to as MTL [20]. We here extend the TSB-HMM learning in a multi-task setting for spectrogram feature-based radar HRRP recognition.

• (d) Full Bayesian inference: We present a fully conjugate Bayesian model structure, which does have an efficient VB solution.

The remainder of this article is organized as follows. We introduce spectrogram feature of HRRP data and analyze its advantages over time domain samples in Section 2. Section 3 briefly reviews the traditional HMMs. In Section 4, the proposed model construction is introduced, and the model learning and classification are implemented based on VB inference. We present experimental results on both single-task and multi-task TSB-HMM with time domain feature and spectrogram feature of measured HRRP data in Section 5. Finally, the conclusions are addressed in Section 6.

2. Spectrogram feature extraction

2.1 Definition of the spectrogram

The spectrogram analysis is a common signal processing procedure in spectral analysis and other fields. It is a view of a signal represented over both time and frequency domains, and has widely been used in the fields of radar signal processing, and speech processing [21, 22], etc.

Spectrograms can readily be created by calculating the short-time Fourier transform (STFT) of the time signal. The STFT transform may be represented as

$$\mathbf{STFT}\left(\tau ,\omega \right)={\int }_{-\infty }^{\infty }x\left(u\right)w\left(u-\tau \right)e^{-j\omega u}du$$

(1)

where x(u) is the signal to be transformed, w(·) is the window function.

The spectrogram is given by the squared magnitude of the STFT function of the signal:

$$\mathbf{Y}\left(\tau ,\omega \right)=|\mathbf{STFT}\left(\tau ,\omega \right){|}^2$$

(2)

From (1) and (2), we can see that the spectrogram function shows how the spectral density of signal varies with time.

In RATR problems, employing some nonlinear transformation (e.g., power transform metric) in feature domain may correct for the departures of samples from normal distribution to some extent and improve average recognition of learning models [2]. The power transform metric is defined as

$$y=x^a;0<a<1$$

(3)

where a is the power parameter.

2.2 Spectrogram feature of HRRP data

The sequential relationship across the range cells within a single HRRP echo can reflect the physical composition of the target. This can be illustrated by Figure 3, which presents the HRRP samples and corresponding spectrogram features from three plane targets, i.e., Yark-42, An-26, and Cessna Citation S/II.

Figure 3

The HRRP sample and corresponding spectrogram feature of plane targets. (a) HRRP echo (time domain feature) of plane targets, (b) corresponding spectrogram feature of original HRRP echoes of plane targets.

The advantages of spectrogram are as follows: (i) HRRP scattering from complex targets are a strong function of the target sensor orientation; and even a slight variation of the target-aspect may yield the scatterers at the edges of the target moving across some range cells [23]. When the target-aspect changes a little, the scatterers within several continuous range cells (referred to as a chunk) are more robust than the scatterers in a single range cell. Therefore, the sequential relationship across the chunks in spectrogram of a single HRRP echo, rather than that across the range cells within a single HRRP, can reflect the target physical composition more robustly. (ii) Spectrogram is a time-frequency representation of a signal. It describes not only time domain feature, but also the spectral density varying with the time. (iii) At each discrete time (each chunk or each range cell), the observation of a spectrogram feature is a vector, while that of time domain feature (HRRP sample) is a scalar. Thus, the high-dimensional feature vector may reflect more details than a single point for discrimination.

3. Review of traditional HMM: finite HMM and infinite HMM

The HMM [24] has widely been used in speech recognition and target recognition. It is a generative representation of sequential data with an underlying markovian process selecting state-dependents distributions, from which observations are drawn. Specially for a sequence of length T, a state sequence s = (s₁, s₂, ..., s _T ) is drawn from P(s _t |s_t-1). Given the observation model f(·), the observation sequence x = (x₁, x₂, ..., x _T ) can then be drawn as $f\left({\theta }_{s_t}\right)$, where ${\theta }_{s_t}$ is a set of parameters for the observation models which is indexed by the state at time t.

An HMM can be modeled as Φ = {w₀, w, θ}, each parameter is defined as

$$\begin{array}{c}{\mathbf{w}}_0=\left\{w_{0,i}\right\},w_{0,i}=Pr\left(s_1=i\right):\mathsf{\text{Initial}}\mathsf{\text{state}}\mathsf{\text{distribution}}\hfill \\ \mathbf{w}=\left\{w_{i,j}\right\},w_{i,j}=Pr\left(s_{t+1}=j|s_t=i\right):\mathsf{\text{State}}\mathsf{\text{transition}}\mathsf{\text{probability}}\hfill \\ \mathbf{\theta }=\left\{{\theta }_i\right\},x_t~f\left({\theta }_{s_t}\right):\mathsf{\text{Emission}}\mathsf{\text{likelihood}}\hfill \end{array}$$

(4)

Given model parameters Φ, the probability of the complete data can be expressed as

$$p\left(\mathbf{x},\mathbf{s}|\mathbf{\Phi }\right)=w_{0,i}\prod _{t=1}^{T-1}{w}_{s_t,s_{t+1}}\prod _{t=1}^{T}p\left(x_t|{\theta }_{s_t}\right)$$

(5)

And the data likelihood p(x|Φ) can be obtained by integrating over the states using the forward algorithm [24]. In a classical HMM [24], the number of the states with the HMM is initialized and fixed. Therefore, we have to specify model structure before learning. However, in many practical applications, it needs an expansive model selection process to obtain a correct model structure. To avoid the model selection process, a fully nonparametric Bayesian approach with countably infinite state spaces is employed, first proposed by Beal, and termed infinite Markov model (iHMM) [25].

Recent study [19] proposes the iHMM with stick-breaking priors (SB-HMM), which can be used to develop an HMM with an unknown number of states. In this model, each row of the infinite state transition matrix w is given a stick-breaking prior. The model is expressed as follows

$$\begin{array}{c}{\mathbf{G}}_i=\sum _{j=1}^{\infty }{w}_{i,j}\delta \left({\theta }_j\right),w_{i,j}=v_{i,j}\prod _{h=1}^{j-1}\left(1-v_{i,h}\right),v_{i,j}|{\beta }_i~\mathbf{Beta}\left(1,{\beta }_i\right)\\ {\beta }_i~\mathbf{G}\mathbf{a}\left(a_{\alpha },b_{\alpha }\right),{\theta }_j~H;i=1,2,\cdots \infty ;j=1,2,\cdots ,\infty \end{array}$$

(6)

where the mixture distribution G _i has weights represented as w _i = [w_{i, 1}, w_{i, 2}, ..., w _{i, j} , ...], δ(θ _j ) is a point measure concentrated at θ _j , Beta(1, β _i ) represents the Beta distribution with hidden variable β _i , the drawn variables ${\left\{v_{i,j}\right\}}_{j=1}^{\infty }$ are independent and identically distributed (i.i.d), Ga(a _α , b _α ) represents the Gamma distribution with preset parameters a _α , b _α , and H denotes a prior distribution from which the set ${\left\{{\theta }_j\right\}}_{j=1}^{\infty }$ is i.i.d drawn. The initial state probability mass function, w₀, is also constructed according to an infinite stick-breaking construction. When w _i terminates at some finite number I - 1 with $w_{i,I}\equiv 1-{\sum }_{J=1}^{I-1}{w}_{i,j}$, this result is a draw from a general Dirichlet distribution [26], which is denoted as w _i ~ GDD(1_{(I-1) × 1}, [β _i ]_{(I-1) × 1}), where 1_{(I-1) × 1}denotes an (I - 1)-length vector of ones, [β _i ]_{(I-1) × 1}is an (I - 1)-length vector of β _i , and I represents the truncation number of states.

The key advantage of stick-breaking construction is that the corresponding state-dependent parameters θ _j are drawn separately, effectively detaching the construction θ from the construction initial-state probability w₀ and state transition matrix w. This is contrast with HDP priors [27], where these matrices are linked with two-level construction. Therefore, the stick-breaking construction makes fast variational inference feasible. In addition, the SB-HMM has a good sparse property, which promotes a sparse utilization of the underlying states [19].

4. The SB-HMM for HRRP data

4.1 MTL model construction

According to the scattering center model [4], for a high-resolution radar system, a radar target does not appear as a "point target" any more, but consists of many scatterers distributed in some range cells along radar LOS. For a certain target, the scattering center model varies throughout the whole target-aspect. Therefore, pre-processing techniques should be applied to the raw HRRP data. In our previous work [11–13], we divide the HRRP into frames according to the aspect-sectors without most scatterers' motion through resolution cells (MTRC), and use distinct parametric models for statistical characterization of each HRRP frame, which are referred to as the aspect-frames and corresponding single-task learning (STL) models in our articles.

For the motivating HRRP recognition problems of interested here, we utilize TSB-HMM for analyzing spectrogram features extracted from HRRP data. For a multi-aspect HRRP sequence of target c (c ∈ {1, ..., C} with C denoting the number of targets here), we divide the data into M _c aspect frames, e.g., the m th set (here m ∈ {1, ..., M _c }) is ${\left\{{\mathbf{x}}^{\left(c,m,n\right)}\right\}}_{n=1}^N$ where N denotes the number of samples in the frame, and x^{(c, m, n)}= [x^{(c, m, n)}(1), ..., x^{(c, m, n)}(L _x )]^Trepresents the n th HRRP sample in the m th frame, with L _x denoting the number of range cells in an HRRP sample. Each aspect frame corresponds to a small aspect-sector avoiding scatters' MTRC [13], and the HRRP samples inside each target-aspect frame can be assumed to be i.i.d. We extract the spectrogram feature of each HRRP sample, and Y^{(c, m, n)}= [y^{(c, m, n)}(1), ..., y^{(c, m, n)}(L _y )] denotes the spectrogram feature of x^{(c, m, n)}as defined in (2) with L _y denoting the number of time bins in spectrogram feature.

If learning a separate TSB-HMM for each frame of the target, i.e., ${\left\{{\mathbf{Y}}^{\left(c,m,n\right)}\right\}}_{n=1}^N$, is termed the single-task TSB-HMM (STL TSB-HMM). Here, we wish to learn a TSB-HMM for all the aspect-frames (tasks) of one target jointly, which is referred to as multi-task TSB-HMM (MTL TSB-HMM). MTL is an approach to inductive transfer that improves generalization by using the domain information contained in the training samples of related tasks as an inductive bias [20]. In our learning problems, the aspect-frames of one target may be viewed as a set of related learning tasks. Rather than building models for each aspect-frame individually (due to target-aspect sensitivity), it is desirable to appropriately share the information among these related data. Therefore, the training data for each task are strengthened and overall recognition performance is potentially improved.

The construction of the MTL TSB-HMM with parameters for target c is represented as

$$\begin{array}{c}{\mathbf{y}}^{\left(c,m,n\right)}\left(l\right)~f\left({\mathbf{\theta }}_{s_l^{\left(c,m,n\right)}}\right);l=1,\dots ,L_y;n=1,\dots ,N;m=1,\dots ,M_c;\hfill \\ s_l^{\left(c,m,n\right)}~\left\{\begin{array}{cc}\hfill {\mathbf{w}}_0^{\left(c,m\right)}\hfill & \hfill \mathsf{\text{if}}l=1\hfill \\ \hfill {\mathbf{w}}_{s_{l-1}^{\left(c,m,n\right)}}^{\left(c,m\right)}\hfill & \hfill \mathsf{\text{if}}l\ge 2\hfill \end{array}\right.\hfill \\ {\mathbf{w}}_i^{\left(c,m\right)}|{\beta }_i~\mathbf{G}\mathbf{D}\mathbf{D}\left({\mathbf{1}}_{\left(I-1\right)\times 1},{\left[{\beta }_i\right]}_{\left(I-1\right)\times 1}\right),{\beta }_i|a_{\alpha },b_{\alpha }~\mathbf{G}\mathbf{a}\left(a_{\alpha },b_{\alpha }\right);i=0,\dots ,I\hfill \\ {\mathbf{\theta }}_i~H;i=1,\dots ,I\hfill \end{array}$$

(7)

where y^{(c, m, n)}(l) is l th time chunk of n th sample's spectrogram in the m th aspect-frame of the c th target, $s_l^{\left(c,m,n\right)}$ denotes the corresponding state indicator, (a _α , b _α ) are the preset hyperparameters. Here, the observation model f(·)is defined as independently normal distribution, and each corresponding element in H(·) is normal-Gamma distribution to preserve conjugacy requirements. Since each time bin of spectrogram feature of a plane corresponds to a fragment of the plane, the HMM states can characterize the frequency domain properties of different fragments of the plane target, i.e., the scattering properties of different physical constructions. A graphical representation of this model is shown in Figure 4a, and Figure 4b depicts that the sequential dependence across time chunks for a given aspect-frame is characterized by an HMM structure.

Figure 4

Graphical representation of the MTL TSB-HMM, the detailed generative process is described in (7). (a) Graph of the MTL TSB-HMM, the detailed generative process is described in (6). Here, ${\mathbf{W}}^{\left(c,m\right)}=\left[{\mathbf{w}}_0^{{\left(c,m\right)}^T};{\mathbf{w}}_1^{{\left(c,m\right)}^T};\dots ;{\mathbf{w}}_I^{{\left(c,m\right)}^T}\right]$, i is the state index, c is the target index, l is the range cell index, y^{(c, m, n)}(l) is l th time chunk of n th sample's spectrogram in the m th aspect-frame of the c th target. For the symbol I, M_c, N, and L _y of a box indicate that there is a stack of such box for the index i, m, n, and l, where I denotes the truncation number of states, M _c denotes the number of aspect-frames of target c, N denotes the number of samples in an aspect frame, and L _y denotes the number of time bins in spectrogram feature. (b) The spatial dependence across time cell of spectrogram feature for a given aspect frame.

The main difference between MTL TSB-HMM and STL TSB-HMM is in the proposed MTL TSB-HMM, all the multi-aspect frames of one target are learned jointly, each of the M _c tasks of target c is assumed to have an independent state-transition statistics, but the state-dependent observation statistics are shared across these tasks, i.e., the observation parameters are learned via all aspect-frames; while in the STL TSB-HMM, each multi-aspect frame of target c is learned separately, therefore, each target-aspect frame builds its own model and the corresponding parameters are learned just via this aspect-frame.

4.2 Model learning

The parameters of proposed MTL TSB-HMM model are treated as variables, and this model can readily be implemented by Markov Chain Monte Carlo (MCMC) [28] method. However, to approximate the posterior distribution over parameters, MCMC requires large computational resources to assess the convergence and reliability of estimates. In this article, we employ VB inference [19, 29, 30], which does not generate a single point estimation of the parameters, but regard all model parameters as possible, with the goal of estimating the posterior density function on the model parameters, as a compromise between accuracy and computational cost for large-scale problems.

The goal of Bayesian inference is to estimate the posterior distribution of model parameters Φ. Given the observation data X and hyper parameters γ, by Bayes' rule, the posterior density for the model parameters may be expressed as

$$p\left(\mathbf{\Phi }|\mathbf{X},\gamma \right)=\frac{p\left(\mathbf{X}|\mathbf{\Phi },\mathbf{\gamma }\right)p\left(\mathbf{\Phi }|\mathbf{\gamma }\right)}{\int p\left(\mathbf{X}|\mathbf{\Phi },\mathbf{\gamma }\right)p\left(\mathbf{\Phi }|\mathbf{\gamma }\right)d\mathbf{\Phi }}$$

(8)

where the denominator ∫p(X|Φ, γ)p(Φ|γ)d Φ = p(X|γ) is the model evidence (marginal likelihood).

VB inference provides a computationally tractable way which seeks a variational distribution q(Φ) to approximate the true posterior distribution of the latent variables p(Φ|X, γ), we obtain the expression

$$logp\left(\mathbf{X}|\mathbf{\gamma }\right)=L\left(q\left(\mathbf{\Phi }\right)\right)+KL\left(q\left(\mathbf{\Phi }\right)\left|\right|p\left(\mathbf{\Phi }|\mathbf{X},\mathbf{\gamma }\right)\right)$$

(9)

where $L\left(q\left(\mathbf{\Phi }\right)\right)=\int q\left(\mathbf{\Phi }\right)log\frac{p\left(\mathbf{X}|\mathbf{\Phi },\mathbf{\gamma }\right)p\left(\mathbf{\Phi }|\mathbf{\gamma }\right)}{q\left(\mathbf{\Phi }\right)}d\mathbf{\Phi }$, and KL(q(Φ)|| p(Φ|X, γ)) is the Kullback-Leibler (KL) divergence between the variational distributions q(Φ) and the true posterior p(Φ|X, γ). Since KL(q(Φ)|| p(Φ|X, γ)) ≥ 0, and it reaches zero when q(Φ) = p(Φ|X, γ), this forms a lower bound for log p(X|γ), so we have log p(X|γ) ≥ L(q(Φ)). The goal of minimizing the KL divergence between the variational distribution and the true posterior is equal to maximize this lower bound, which is known as the negative free energy in statistical physics.

For the computational convenience and intractable of the negative-free energy, we assume a factorized q(Φ), i.e., $q\left(\mathbf{\Phi }\right)=\prod _k{q}_k\left({\phi }_k\right)$, which has same form as employed in p(Φ|X, γ). With this assumption, the mean field approximation of the variational distributions for the proposed MTL TSB-HMM with target c may be expressed as

$$q\left(\mathbf{\Phi }\right)=\prod _{m=1}^{M_c}\prod _{i=0}^{I}q\left({\mathbf{w}}_i^{\left(c,m\right)}\right)\left[\prod _{m=1}^{M_c}\prod _{n=1}^N\prod _{l=1}^{L_y}q\left(s_l^{\left(c,m,n\right)}\right)\right]q\left(\mathbf{\theta }\right)q\left(\mathbf{\beta }\right)$$

(10)

where $\left\{{\left\{{\mathbf{w}}_i^{\left(c,m\right)}\right\}}_{m=1,i=0}^{M_c,I},{\left\{s_l^{\left(c,m,n\right)}\right\}}_{m=1,n=1,l=1}^{M_c,N,L_y},\mathbf{\theta },\mathbf{\beta }\right\}$ are the latent variables in this MTL model.

A general method for performing variational inference for conjugate-exponential Bayesian networks outlined in [17] is as follows: for a given node in a graphic model, write out the posterior as though everything were known, take the logarithm, the expectation with respect to all known parameters and exponentiate the result. We can implement expectation-maximization (EM) algorithm in variational inference. The lower bound is increased in each of iteration until the algorithm converges. In the following experiments, we terminate EM algorithm when the changes of the lower bound can be neglected (the threshold is 10^-6). Since it requires computational resource comparable to EM algorithm, variational inference is faster than MCMC methods. The detailed update equations for the latent variables and hyperparameters of MTL TSB-HMM with HRRP spectrogram feature are summarized in the Appendix.

4.3 Main procedure of radar HRRP target recognition based on the proposed MTL TSB-HMM algorithm

The main procedure of radar HRRP target recognition based on the proposed MTL TSB-HMM algorithm is shown as follows.

4.3.1. Training phase

• (1) Divide the training samples of target c (c = 1, 2, ..., C) into HRRP frames ${\left\{{\mathbf{x}}^{\left(c,m\right)}\right\}}_{m=1}^{M_c}$, where M _c is the number of tasks of target c, ${\mathbf{x}}^{\left(c,m\right)}={\left\{{\mathbf{x}}^{\left(c,m,n\right)}\right\}}_{n=1}^N$ denotes the m th range aligned and amplitude normalized HRRP frame, N is the number of echoes a frame contains.

• (2) Extract the spectrogram feature ${\left\{{\mathbf{Y}}^{\left(c,m,n\right)}\right\}}_{m=1,n=1}^{M_c,N}$ of each HRRP sample with Y^{(c, m, n)}= [y^{(c, m, n)}(1), y^{(c, m, n)}(2), ..., y^{(c, m, n)}(L _y )] denoting the spectrogram feature of x^{(c, m, n)}as defined in (2).

• (3) For each target, we construct an MTL TSB-HMM model, and learn the parameters of $w_{0,s_1^{\left(c,m\right)}}^{\left(c,m\right)}$, $w_{i,j}^{\left(c,m\right)}$, and θ _i for all aspect-frames of the target via using spectrogram feature, where $w_{0,s_1^{\left(c,m\right)}}^{\left(c,m\right)}$ is the initial state probability for the index frame m of target c, $w_{i,j}^{\left(c,m\right)}$ is state transition probability from state i to the j for the index frame m of target c, and θ _i are the parameters of observation model associated with corresponding state i (c ∈ {1, ..., C}, m ∈ {1, ..., M _c }, i, j ∈ {1, ..., I}). The detailed learning procedure of the parameters of MTL TSB-HMM with HRRP spectrogram feature are discussed in Section 4.3 and the Appendix.

• (4) Store the parameters of initial state probability ${\left\{w_{0,s_1^{\left(c,m\right)}}^{\left(c,m\right)}\right\}}_{m=1}^{M_c}$, state transition probability ${\left\{{}_{\mathbf{w}}^{\left(c,m\right)}\right\}}_{m=1}^{M_c}$ and the parameters of observation model ${\left\{{\mathbf{\theta }}_i\right\}}_{i=1}^I$ for each target c with c = 1, 2, ..., C.

4.3.2. Classification phase

• (1) The amplitude normalized HRRP testing sample is time-shift compensated with respect to the averaged HRRP of each frame model via the slide correlation processing [23].

• (2) Extract the spectrogram feature ${\left\{{\mathbf{Y}}_{\mathsf{\text{test}}}^{\left(c,m\right)}\right\}}_{c=1,m=1}^{C,M_c}$of the slide-correlated HRRP testing sample x_test, where ${\mathbf{Y}}_{\mathsf{\text{test}}}^{\left(c,m\right)}=\left\{\left[{\mathbf{y}}_{\mathsf{\text{test}}}^{\left(c,m\right)}\left(1\right),{\mathbf{y}}_{\mathsf{\text{test}}}^{\left(c,m\right)}\left(2\right),\dots ,{\mathbf{y}}_{\mathsf{\text{test}}}^{\left(c,m\right)}\left(L_y\right)\right]\right\}$ denotes the spectrogram feature of HRRP testing sample correlated with m th frame of target c as defined in (2).

• (3) The frame-conditional likelihood of target can be calculated as

  $$p\left({\mathbf{Y}}_{\mathsf{\text{test}}}^{\left(c,m\right)}|c,m\right)=\sum _{\mathbf{s}}〈w_{0,s_1^{\left(c,m\right)}}^{\left(c,m\right)}〉\prod _{l=2}^{L_y}〈w_{s_{l-1}^{\left(c,m\right)},s_l^{\left(c,m\right)}}^{\left(c,m\right)}〉\prod _{l=1}^{L_y}{f}^{\left(c,m\right)}\left({\mathbf{y}}_{\mathsf{\text{test}}}^{\left(c,m\right)}\left(l\right)|〈{\mathbf{\theta }}_{s_l^{\left(c,m\right)}}〉\right)$$

  (11)

where 〈·〉 means the posterior expectation for the latent variable over the corresponding distribution on it, e.g., $〈w_{0,s_1^{\left(c,m\right)}}^{\left(c,m\right)}〉$ denotes the posterior expectations of initial state probability, $〈w_{s_{l-1}^{\left(c,m\right)},s_l^{\left(c,m\right)}}^{\left(c,m\right)}〉$ denotes the posterior expectations of state transition probability from state $s_{l-1}^{\left(c,m\right)}$ to the $s_l^{\left(c,m\right)}$ for the frame m of target c, and $〈{\mathbf{\theta }}_{s_l^{\left(c,m\right)}}〉$ denotes the posterior expectations of the observation model parameters associated with state $s_l^{\left(c,m\right)}$, with the corresponding state indicator for the l th time chunk $s_l^{\left(c,m\right)}\in \left\{1,\cdots ,I\right\}$. Then, $p\left({\mathbf{Y}}_{\mathsf{\text{test}}}^{\left(c,m\right)}|c,m\right)$ can be calculated by forward-backward procedure [24] for each m (m ∈ {1, ..., M _c }) and c (c ∈ {1, ..., C}).

• (4) We calculate the class-conditional likelihood $p\left({\mathbf{Y}}_{\mathsf{\text{test}}}^{\left(c\right)}|c\right)$ for each target c

  $$p\left({\mathbf{Y}}_{\text{test}}^{\left(c\right)}|c\right)=\underset{m}{max}p\left({\mathbf{Y}}_{\text{test}}^{\left(c,m\right)}|c,m\right);m=1,\cdots ,M_c$$

  (12)

• (5) As discussed in Section 1, the testing HRRP sample will be assigned to the class with the maximum class-conditional likelihood, with the assumption that the prior class probabilities are same for all targets of interests,

  $$k=arg\underset{c}{max}p\left({\mathbf{Y}}_{\mathsf{\text{test}}}^{\left(c\right)}|c\right);c=1,\cdots ,C$$

  (13)

5. Experimental results

5.1 Measured data

We examine the performance of the TSB-HMM on the 3-class measured data, including three actual aircrafts (a propeller plane An-26, a small yet plane Cessna Citation S/II, and a big yet plane Yark-42), the radar works on a C band with bandwidth of 400 MHz, the range resolution of the HRRP is about 0.375 m.

The parameters of the targets and radar are shown in Table 1, and the projections of target trajectories onto the ground plane are displayed in Figure 5, from which the aspect angle of the airplane can be estimated according to its relative position to radar. As shown in Figure 5, all aspects of targets were measured repeatedly several times in this dataset. The requirements of choosing training data and test data are that the training data and the test data are from different data segments, and the training data cover almost all of target-aspect angles of the test data, but their elevation angles are different. The second and the fifth segments of Yark-42, the sixth and the seventh segments of Cessna Citation S/II and the fifth and the sixth segments of An-26 are taken as the training samples while the remaining data are left for testing. These training data almost cover all of the target-aspect angles. Also, we need test data from different target to measure the rejection performance of our model. Here, we use 18,000 truck HRRP samples generated by the electromagnetic simulator software, XPATCH, as a confuser target. In addition, the HRRP samples are 128-dimensional vectors.

Table 1 Parameters of planes and radar in the experiment

Figure 5

Projections of target trajectories onto the ground plane: (a) An-26. (b) Cessna Citation S/II. (c) Yark-42.

As discussed in the literature [11, 12], it is a prerequisite for radar target recognition to deal with the target-aspect, time-shift, and amplitude-scale sensitivity. According to radar parameters and the condition of aspect sectors without MTRC, for training data from 3 targets we totally have 135 HRRP frames, of which 35 from Yark-42, 50 from Cessna Citation S/II and 50 from An-26. Similar to the previous study [11–13], HRRP training samples should be aligned by the time-shift compensation techniques in ISAR imaging [23] to avoid the influence of time-shift sensitivity. Each HRRP sample is normalized by L₂ normalization algorithm to avoid the amplitude-scale sensitivity. In the rest of the article, the training HRRPs in each frame are assumed to have been aligned and normalized.

Nine training datasets are considered for training: 1 × 135, 2 × 135, 4 × 135, 8 × 135, 16 × 135, 32 × 135, 64 × 135, 128 × 135 and 1024 × 135, where 2 × 135 means 2 HRRP samples randomly selected from each of the 135 target-aspect frames, i.e., there are totally 2 × 135 = 270 HRRP training samples, similar meaning for other size of HRRP datasets. Since MTL needs load the whole HRRP training dataset of a target to share the information among them, it requires more memory resource than STL. Due to the limited memory resource of our computer, we do not consider 1024 × 135 training dataset in MTL. Since there is no prior knowledge about how many states we should use, and how to set these states, the HMM states are not manually set like [31]. We set a large truncation number in our model to learn the meaningful states automatically. In the following experiments, we set the truncation level I to 40 for both spectrogram feature and time domain feature in STL, and I to 60 for both spectrogram feature and time domain feature in MTL. Similar results were found for lager truncations. In our model, since the parameter β _i controls the prior distribution on the number of states, we set the hyper-parameters a _α = b _α = 10^-6 for each β _i to promote sparseness on states.

5.2 Time domain feature versus spectrogram feature

In this experiment, STL TSB-HMM and MTL TSB-HMM of HRRP training datasets within each frame are learned, respectively, and the two features, i.e., time domain and spectrogram features, are compared. When using the HRRP time domain feature, we can just substitute the scalar x^{(c, m, n)}(l) for the vector y^{(c, m, n)}(l) in (6), where x^{(c, m, n)}(l) represents the n th HRRP sample in the m th frame of target c, with l denoting the corresponding range cell in the HRRP sample.

Figure 6 shows that the performances of STL TSB-HMMs based on time domain feature are better than those based on spectrogram feature when training dataset no more than 32 × 135. The reason is that more parameters need to be estimated for the model with spectrogram feature than that with time domain feature. For example, for a target-aspect frame with 32 training samples, we need to estimate the 40 16-dimensional states for the model with spectrogram feature; while 25 1-dimensional states for the model with time domain feature. When training data size larger than 32 × 135, spectrogram feature obtains obviously better performance. Table 2 further compares the confusion matrices and average recognition rates of using time domain feature and spectrogram feature with 1024 × 135 training samples via STL. We can clearly find that the average recognition rate obtained by spectrogram feature is about 6.6% points larger than that obtained by the time domain feature. The performances of MTL TSB-HMMs are shown in Figure 7. Since MTL sharing states between different tasks of a target, which is better for parameter learning with small training data size, spectrogram domain feature outperforms time domain feature even with few training data.

Figure 6

Average recognition rates of STL TSB-HMMs versus the number of training samples for time domain feature and spectrogram feature.

Table 2 Confusion matrices and average recognition rates of time domain feature and spectrogram feature based on STL TSB-HMMs with 1,024 training samples per target-aspect frame

Figure 7

Average recognition rates of MTL TSB-HMMs versus the number of training samples for time domain feature and spectrogram feature.

The posterior state distributions of MTL TSB-HMM with spectrogram feature for all the three plane targets with 128 × 135 training data are shown in Figure 8. In this example, the state truncation level of I = 60 is employed for each plane. Further for each plane, the 60 hidden states are shared across all aspect-frames, and there are 46, 48, and 49 meaningful states with the posterior state usage larger than zero for Yark-42, Cessna Citation S/II, and An-26, respectively, and those of other 14, 12, and 11 states are zero, which justifies using the truncated version stick-breaking prior for our data.

Figure 8

The hidden state distributions for all the three plane targets with 128 × 135 training data, generated via a run of the MTL TSB-HMM with spectrogram feature. In this example, the state truncation level of I = 60 is employed for each plane. For each plane, 60 hidden states are shared across all aspect-frames, and 46, 48, and 49 meaningful states are inferred for Yark-42, Cessna Citation S/II, and An-26, respectively.

Next, we consider the target rejection problem. Three planes targets are considered as "in-class targets", while the simulated data consists of 18,000 truck HRRP samples are considered as confuser targets. Two examples of confuser targets HRRP samples are shown in Figure 9. Our goal is to test whether a new data is in the family of the in-class targets or not. Figure 10 presents the rejection performance evaluated by the receiver operation characteristic (ROC) curves. The ROC curve depicts the detection probability versus the false alarm probability. For a fixed false alarm probability, a method with the higher detection probability is better. The dataset size is selected as 128 × 135 in the training phase. The ROC curves are shown in Figure 10. The spectrogram feature-based TSB-HMMs outperforms the time domain feature-based TSB-HMMs, especially for Yark-42. Figure 11 shows the test likelihoods of 1,200 Yark-42 samples and 18,000 confuser target samples obtained with STL TSB-HMM. As shown in Figure 11a, when using the time domain feature, the test likelihoods of Yark-42 samples are relative low, and many confuser samples have higher test likelihoods than Yark-42 samples. That is to say, when we set a high discrimination threshold, the detection probability is very low, and the false alarm probability is high. By contrast, from Figure 11b, we can find that when using the spectrogram feature, the test likelihoods of Yark-42 samples are higher than most of the test likelihoods of the confuser samples. Therefore, the detection performance of the spectrogram feature is much better than that of the time domain feature for Yark-42.

Figure 9

Two HRRP samples of simulated scattering data of truck targets.

Figure 10

Rejection results for different sets of targets using STL TSB-HMMs with different features. (a) STL TSB-HMM with time domain feature (average AUC: 0.738). (b) STL TSB-HMM with spectrogram feature (Average AUC: 0.924).

Figure 11

The test likelihoods of 1,200 Yark-42 samples and 18,000 truck samples (confuser target samples) obtained with STL TSB-HMM. The STL TSB-HMM is trained with the corresponding feature extracted from 4480 (128 × 35) Yark-42 training samples. (a) STL TSB-HMM with the time domain feature. (b) STL TSB-HMM with the spectrogram feature.

5.3 STL versus MTL

In order to model spectrogram feature, two parameters need to be set first, i.e., the width of window function and the length of the overlapped window.

In the feature space, the spectrogram varies with the width of the window function and the overlap across the windows. For HRRP data analysis, a wide window function provides better frequency resolution, but worsens the time resolution, and vice versa. Physically, since the width of a window function determines the length of segments of a target, the longer the segment we divide a target into, the more physical composition of the target will be contained in the each component of the observation vector; meanwhile, the overlap across the windows determines the redundancy of the segments.

We build a set of MTL TSB-HMMs to search these two parameters. The width of the window function is chosen from 10 to 40 range cells with an increment of 1 range cells, and the overlap length is fixed as typically half width of the window function. In this experiment, we use 64 × 135 training samples. As demonstrated in Figure 12, the optimal width of the window function is 33 range cells. We then fix the optimal window function width, and set the overlap length from 1 range cell to 29 range cells to determine the optimal overlap length. From Figure 13, the optimal overlap length is 16 range cells. Therefore, we extract the spectrogram feature with the window function width of 33 range cells and the overlap length of 16 range cells for training.

Figure 12

Recognition rates of spectrogram features with different width of the window function via MTL TSB-HMMs. The overlap length is fixed as typically half width of the window function.

Figure 13

Average recognition rate of spectrogram features with different overlap length via MTL TSB-HMMs. The width of the window function is 33 range cells.

We compare two methods for spectrogram feature based TSB-HMM: (i) the proposed MTL TSB-HMMs method, for which we learn target-aspect frames of a target collectively; (ii) the STL TSB-HMMs method, for which each target-aspect frame of targets modeled separately. As shown in Figure 14, the proposed MTL TSB-HMMs method consistently outperforms the STL TSB-HMMs method and the improvement is more significant when there is only a small amount of training data available. This is because MTL exploits the sharing states between different tasks and uses the sharing information to enhance the overall performance. In addition, in the state truncation level of I = 60 is employed for each of the three planes. For each plane, the 60 states are shared across the aspect-frames. Therefore, we only impose 60 × 3 = 180 states in the MTL TSB-HMMs. However, in the training phase of STL model, the state truncation level of I = 40 is employed for each aspect-frame. As discussed in Section 5.1, we have 50 + 50 + 35 = 135 aspect-frames from the three targets. Therefore, we totally impose 40 × 135 = 5400 states in the STL TSB-HMMs. Table 3 summarizes the confusion matrix and average recognition performance of STL TSB-HMMs and MTL TSB-HMMs, with 2 × 135 and 128 × 135 training samples. Note that when we only use 2 × 135 samples for training, the average recognition rate of STL TSB-HMMs is only 52.7%, while the average recognition rate of MTL TSB-HMMs is 88.0%. With the training data increasing, the performance of STL becomes close to that of MTL. The average recognition rate obtained by MTL is about 2.3% higher than that by STL for 128 × 135 training samples.

Figure 14

Average recognition rates of the STL TSB-HMMs and MTL TSB-HMMs with spectrogram feature versus the number of training samples.

Table 3 Confusion matrices and average recognition rates of spectrogram feature based on STL TSB-HMMs and MTL TSB-HMMs with 2 training samples and 128 training samples per target-aspect frame

We also consider the target rejection problem here. The ROC curves of MTL TSB-HMMs are presented in Figure 15. Compare with Figure 10, the area under curve (AUC) of STL are slightly larger than that of MTL.

Figure 15

Rejection results for different sets of targets using MTL TSB-HMMs with different features. (a) MTL TSB-HMM with time domain feature (average AUC: 0.724). (b) MTL TSB-HMM with spectrogram feature (Average AUC: 0.905).

5.4 MTL with power transformed spectrogram feature

We employ the power transform metric for all the observation vectors (the vectors in each fixed time) of the spectrogram. Figure 16 show that the optimal parameter a* = 0.4. Figure 17 shows the performance of MTL TSB-HMMs is better than STL-HMMs for power transformed spectrogram feature. Compared with the average recognition rate of original spectrogram feature in Figure 14, that of power transformed spectrogram feature shown in Figure 17 are much larger, especially for small training data sets.

Figure 16

Average recognition rate of spectrogram features with different power transformed via MTL TSB-HMMs.

Figure 17

Average recognition rates of the STL TSB-HMMs and MTL TSB-HMMs with spectrogram feature transformed by optimal power (0.4) versus the number of training samples.

The confusion matrix and average recognition rates of STL TSB-HMMs and MTL TSB-HMMs with power transformed spectrogram feature are shown in Table 4, where 2 × 135 and 128 × 135 training samples are used for learning the models. Note that, when we consider 2 × 135 training samples for MTL TSB-HMMs based on spectrogram feature with the optimal power transformation, the average of recognition rate is nearly equivalent to the case of considering 128 × 135 training samples for that based on original spectrogram feature. When using 128 training samples per target-aspect frame for STL TSB-HMMs, the average recognition rate is gained by 4.3% via power transformation; while for MTL TSB-HMMs the gain is 2.5%.

Table 4 Confusion matrices and average recognition rates of power transformed spectrogram feature (power parameter is 0.4) based on STL TSB-HMMs and MTL TSB-HMMs with 2 training samples and 128 training samples per target-aspect frame

Similarly, as shown in Figure 18, we obtained the ROC curve for power transformed spectrogram feature in the same experimental environment as that we mentioned in Section 5.3. The AUC of STL and MTL with transformed spectrogram features are gained by 3.3 and 5.0%; therefore, the model of using the transformed spectrogram features can improve the rejection performance.

Figure 18

Rejection results for different sets of targets using STL TSB-HMMs and MTL TSB-HMMs with spectrogram feature transformed via the optimal power. (a) STL TSB-HMMs (average AUC: 0.955). (b) MTL TSB-HMMs (Average AUC 0.955).

5.5 Computation burden

All experiments have been performed in nonoptimized programme written in Matlab, on a Pentium PC with 3.2-GHz CPU and 2 GB RAM. In our VB algorithm, when the relative change of lower bound between two consecutive iterations is less than the threshold 10^-6, we believe our algorithm converges. Generally, in a practical application, the larger training dataset requires the huger computational burden in the training phase. When the training dataset contains 128 × 135 training samples, the VB algorithm of the MTL TSB-HMM with time domain feature and the truncation number I = 60 converges after about 400 iterations and requires about 14 h, and the VB algorithm of MTL TSB-HMMs with spectrogram feature and the truncation number I = 60 converges after about 200 iterations and requires about 2 h. Although the above computation is pretty expensive, we know that the computation cost in the training phase can be ignored for an off-line learning (or training) system, and it is more important to evaluate the computation cost in the classification phase. The MTL TSB-HMMs with time domain feature and spectrogram feature require 0.6680 and 1.5893 s, respectively, to match a test sample with all frame models. The computation time given here is averaged over ten runs.

6. Conclusion

We have utilized spectrogram feature of HRRP data and presented an MTL-based hidden Markov model with truncated stick-breaking prior (MTL TSB-HMM) for radar HRRP target recognition. The construction of this model allows VB inference, which extremely decreases the computational burden.

After resolving the three sensitivity problems, i.e., the target-aspect, time-shift, and amplitude scale sensitivity of HRRP, respectively, we first compare the spectrogram feature of HRRP with the time domain feature of HRRP data via single-task and multi-task learning-based hidden Markov model with truncated stick-breaking prior (STL and MTL TSB-HMM). Second, we measure the performance of STL TSB-HMM and MTL TSB-HMM with spectrogram feature, where in MTL TSB-HMM, the multiple tasks are linked by different target-aspects. Finally, we introduce power transformation metric to improve the recognition performance of spectrogram feature. It is shown that using spectrogram feature not only have a better ROC, but also obtain a better recognition performance than using HRRP time domain feature. MTL shares the underlying state information among different target-aspects, and can provide a better recognition performance compared to STL. In addition, the power transformation metric can enhance both average recognition rate and ROC. It is worth to point out that our MTL model with spectrogram feature can obtain a good recognition performance with much fewer training data compared with the conventional radar HRRP-based statistical recognition methods, which is a good property for RATR.

Appendix

Derivation of update equations in VB approach

1. Update q(θ _i ):

For the MTL TSB-HMM model introduced in Section 4.3 and the corresponding mean-field variational distribution described in (9), q(θ _i ) is defined by the specific application:

$$logq\left({\mathbf{\theta }}_i\right)\propto logH\left({\mathbf{\theta }}_i|{\tilde{\mathbf{\theta }}}_i\right)+\sum _{m=1}^{M_c}\sum _{n=1}^N\sum _{l=1}^{L_y}<s_{l,i}^{\left(c,m,n\right)}>logf\left({}_{\mathbf{y}}^{\left(c,m,n\right)}\left(l\right)|{\mathbf{\theta }}_i\right)i=1,\dots ,I$$

(14)

where q(θ _i ) is defined by the specific application, and ${\tilde{\mathbf{\theta }}}_i$ denotes the parameters in the q(θ _i ), $<s_{l,i}^{\left(c,m,n\right)}>$ represents the expected number of state indicator $s_l^{\left(c,m,n\right)}$ with outcome i. If model is conjugate exponential, that is, H(·) is conjugate to the likelihood f(·), we can readily obtain the update equation for q(θ _i ).

2. Update q(β _i )

In MTL model, we assume $q\left({\beta }_i\right)=\mathbf{G}\mathbf{a}\left({\tilde{a}}_{\alpha }^{\left(i\right)},{\tilde{b}}_{\alpha }^{\left(i\right)}\right)$, the updating equation for q(β _i ) with m = 1, ..., M _c are expressed as follows:

$$\left\{\begin{array}{c}{\tilde{a}}_{\alpha }^{\left(i\right)}=a_{\alpha }+M_c\left(I-1\right)\hfill \\ {\tilde{b}}_{\alpha }^{\left(i\right)}=b_{\alpha }-\sum _{m=1}^{M_c}\sum _{j=1}^{I-1}\left[\psi \left({\tilde{\beta }}_{i,j}^{\left(c,m\right)}{|}_2\right)-\psi \left({\tilde{\beta }}_{i,j}^{\left(c,m\right)}{|}_1+{\tilde{\beta }}_{i,j}^{\left(c,m\right)}{|}_2\right)\right]\hfill \end{array}\right.;i=1,\cdots ,I$$

(15)

where ψ(·) is the digamma function.

3. Update $q\left({\mathbf{w}}_0^{\left(c,m\right)}\right)$ and $q\left({\mathbf{w}}_i^{\left(c,m\right)}\right)$

For the given prior $p\left({\mathbf{w}}_0^{\left(c,m\right)}\right)=\mathbf{G}\mathbf{D}\mathbf{D}\left({\mathbf{1}}_{\left(I-1\right)\times 1},{\left[{\beta }_0\right]}_{\left(I-1\right)\times 1}\right)$ and $p\left({\mathbf{w}}_i^{\left(c,m\right)}\right)=\mathbf{G}\mathbf{D}\mathbf{D}\left({\mathbf{1}}_{\left(I-1\right)\times 1},{\left[{\beta }_i\right]}_{\left(I-1\right)\times 1}\right)$, with m = 1, ..., M _c and i = 1, ..., I, assume $q\left({\mathbf{w}}_0^{\left(c,m\right)}\right)=\mathbf{G}\mathbf{D}\mathbf{D}\left({\tilde{\mathbf{\beta }}}_0^{\left(c,m\right)}{|}_1,{\tilde{\mathbf{\beta }}}_0^{\left(c,m\right)}{|}_2\right)$ and $q\left({\mathbf{w}}_i^{\left(c,m\right)}\right)=\mathbf{G}\mathbf{D}\mathbf{D}\left({\tilde{\mathbf{\beta }}}_i^{\left(c,m\right)}{|}_1,{\tilde{\mathbf{\beta }}}_i^{\left(c,m\right)}{|}_2\right)$. The updating equation for $q\left({\mathbf{w}}_0^{\left(c,m\right)}\right)$ and $q\left({\mathbf{w}}_i^{\left(c,m\right)}\right)$ are given as follows:

$$\begin{array}{c}\left\{\begin{array}{c}\hfill {\tilde{\beta }}_{0,j}^{\left(c,m\right)}{|}_1=1+\sum _{n=1}^N<s_{1,j}^{\left(c,m,n\right)}>\hfill \\ \hfill {\tilde{\beta }}_{0,j}^{\left(c,m\right)}{|}_2=\frac{{\tilde{a}}_{\alpha }^{\left(0\right)}}{{\tilde{b}}_{\alpha }^{\left(0\right)}}+\sum _{n=1}^N\sum _{h=j+1}^I<s_{1,h}^{\left(c,m,n\right)}>\hfill \end{array}\right.;j=1,\dots ,I-1\hfill \\ \left\{\begin{array}{c}\hfill {\tilde{\beta }}_{i,j}^{\left(c,m\right)}{|}_1=1+\sum _{n=1}^N\sum _{l=2}^{L_y}<s_{l-1,i}^{\left(c,m,n\right)}><s_{l,j}^{\left(c,m,n\right)}>\hfill \\ \hfill {\tilde{\beta }}_{i,j}^{\left(c,m\right)}{|}_2=\frac{{\tilde{a}}_{\alpha }^{\left(i\right)}}{{\tilde{b}}_{\alpha }^{\left(i\right)}}+\sum _{n=1}^N\sum _{h=j+1}^I\sum _{l=2}^{L_y}<s_{l-1,i}^{\left(c,m,n\right)}><s_{l,h}^{\left(c,m,n\right)}>\hfill \end{array}\right.;i=1,\cdots ,I;j=1,\dots ,I-1\hfill \end{array}$$

(16)

4. Update log q(s)

Given the approximate distribution of the other variables, the update equation for q(s^{(c, m, n)}) are given as follows:

$$\begin{array}{c}logq\left({\mathbf{s}}^{\left(c,m,n\right)}\right)\propto <log{\mathbf{w}}_{0,s_1^{\left(c,m,n\right)}}^{\left(c,m\right)}>+\sum _{l=1}^{L_y-1}<log{\mathbf{w}}_{s_l^{\left(c,m,n\right)},s_{l+1}^{\left(c,m,n\right)}}^{\left(c,m\right)}>+\sum _{l=1}^{L_y}<logf\left({}_{\mathbf{y}}^{\left(c,m,n\right)}\left(l\right)|{\mathbf{\theta }}_{s_l^{\left(c,m,n\right)}}\right)>;\hfill \\ m=1,\dots ,M_c,n=1,\dots ,N\hfill \end{array}$$

(17)

where < · > denotes the expectation of the associated variables function. One may derive that

$$\begin{array}{c}<log{\mathbf{w}}_{i,j}^{\left(c,m\right)}>=\sum _{h=1}^{j-1}\left[\psi \left({\tilde{\beta }}_{i,h}^{\left(c,m\right)}{|}_2\right)-\psi \left({\tilde{\beta }}_{i,h}^{\left(c,m\right)}{|}_1+{\tilde{\beta }}_{i,h}^{\left(c,m\right)}{|}_2\right)\right]+\left[\psi \left({\tilde{\beta }}_{i,j}^{\left(c,m\right)}{|}_1\right)-\psi \left({\tilde{\beta }}_{i,j}^{\left(c,m\right)}{|}_1+{\tilde{\beta }}_{i,j}^{\left(c,m\right)}{|}_2\right)\right];\\ i=0,\dots ,I,j=1,\dots ,I\end{array}$$

(18)