# Multitarget Parameter Estimation of Monopulse Radar Based on RJ-MCMC Algorithm

Authors: Jian Gong; Yiduo Guo; Hui Yuan; Qun Wan
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

A multiple parameter estimation method based on RJ-MCMC for multiple nondiscernible targets is proposed in this paper. Different from the traditional estimation methods, the proposed method can simultaneously complete the joint estimation of the target number and the target location parameters. More importantly, the method proposed in this chapter is applicable to many situations with different power and nondistinguishable target. The simulation results show that the method proposed in this chapter requires less observation time to obtain similar and even better estimation performance than the ML-MDL method, which is of great significance for real-time processing.


---

## Body

## 1. Introduction

In the process of fighter penetration, as shown in Figure1, a two-tier echelon and a team's wedge formation are usually used [1].Figure 1
The fighting order of fighter plane.
(a)
Two-shipper echelon (b)
“Wedge” team vertical view (c)
“Wedge” team flat viewIt can be seen that two or more than two targets are likely to exist in the same distance-angle resolution unit of the radar when the target is launched with a double “Echelon” or a team ’s “wedge” formation. With the influence of side lobe jamming and false target jamming [2–5], radar’s display on the two-dimensional display of distance angle is shown in Figure 2. The red pentagram represents the target that radar is detecting or tracking, and the dot represents other targets or false targets.Figure 2
Radar range-angle two-dimensional display diagram.From Figure2, we can see that, during the penetration process, many false targets may be observed in the range angle dimension of radar. The black dots represent the false targets in the sidelobe region, which can be suppressed by the side lobe countermeasures. The blue dots represent the false targets in the different range resolution units in the radar main beam resolution unit and can be identified or suppressed by the method of distance dimension false target interference. However, the false targets (green dots) in the same distance and angle resolution units are still possible to cause interference to the radar angle measurement and tracking system. Serious deviation will occur when the angle of the target is estimated by using the traditional single pulse processing.From the spatial domain, the composite jamming described above can deceive monopulse radar from a multisource coherent/incoherent angle. The solution to this problem from the beam is the angle superresolution method based on the improved monopulse processing method, and it mainly includes complex monopulse processing method [6–10], parameter estimation method based on joint bin processing (JBP) [11, 12], and maximum likelihood (ML) estimation. The complex monopulse processing method assumes that the target just falls on the sampling point of matched filter and only uses one matched filter to sample information, so it can only solve the problem of angle estimation of up to two indistinguishable targets. The ML estimation algorithm based on JBP model is computationally intensive and easy to fall into local extremum. When the number of samples is small and the SNR is low, the minimum description length (MDL) criterion is used to estimate the number of targets in the beam and it is very easy to fail. The standard MCMC method cannot “jump” between subspaces of different dimensions, so it does not have the ability to estimate the number of targets. A direct solution is similar to the method proposed by [11, 12]. In order to get the target parameters estimation, the order MCMC sampling algorithm is run independently and the number of targets is fixed at each time. Then the number of targets is judged by MDL criterion. However, there is no substantial improvement in the amount of computation and estimation accuracy. Therefore, in this paper, we use RJ-MCMC method to jointly estimate target number and target position parameters. RJ-MCMC method can “jump” from different dimension space when the dimension of parameter space changes; that is, it can directly sample according to the joint distribution of A and finally get the estimation results of target number and target position parameters.The remainder of this paper is organized as follows. In Section2, the joint bin processing (JBP) model is established. Then, the prior and posterior distributions of Bayesian parameter estimation are derived in Section 3, followed by a multitarget parameter estimation algorithm for monopulse radar Based on RJ-MCMC is proposed in Section 4. Some simulations are conducted to verify the performance of the proposed method in Section 5. Finally, we conclude the paper in Section 6.


## 2. Joint Bin Processing Model

Before giving the model, we first make the following assumptions:(I) Assuming that the target model is a point target model, the extended target is not considered.(II) When the noise is not considered, the output of matched filter is triangular waveform; that is, the radar transmits the signal with rectangular envelope.(III) The sampling interval after matched filtering is pulse widthT, which means that a target can only appear in its adjacent matched filtering sampling points.Figure3 shows the “leakage” model of two targets at two adjacent sampling points. ΔT represents the time offset of the real target position relative to the sampling point 1. x=acos(ϕ) (or x=asin(ϕ)) is the peak amplitude of the target echo in the same phase (or quadrature) channel after matched filtering. If α=ΔT/T is defined, the amplitude of the target echo at the first and second sampling points is (1-α)x and αx, respectively.Figure 3
Two targets at two adjacent sampling points “leakage” model.Suppose that there arek targets between the two sampling points. When the radar launches the first pulse, we can get the target’s observed values at two sampling points, as shown in Formula (1):(1)si1m=∑jk1-αjxjm+nsi1msi2m=∑jkαjxjm+nsi2mdai1m=∑jk1-αjηajxjm+ndai1mdai2m=∑jkαjηajxjm+nai2mdei1m=∑jk1-αjηejxjm+nei1mdei2m=∑jkαjηejxjm+nei2mwhere, {αj}j=1k, {ηaj}j=1k, and {ηej}j=1k are the parameters to be estimated, respectively, thek target in the subgate in relative distance, azimuth, and elevation information. si1(m)m=1M, dai1(m)m=1M, and dei1(m)m=1M are the observation data for sum, azimuth, and elevation difference channels for sampling point 1. si2(m)m=1M, dai2(m)m=1M, and dei2(m)m=1M are the observation data for sum, azimuth, and elevation difference channels for sampling point 2. xj(m) is the matching pulse output of them pulse of thej target. xj(m)m=1M obeys the Gauss distribution of zero mean and σtj2 variance.The observation model given by Formula (1) is expressed as vector matrix form:(2)Z:,m=simsqmdaimdaqmdeimdeqmT=DX:,m+W:,mwhere(3)D=1-α11-α2⋯1-αkα1α2⋯αk1-α1ηa11-α2ηa2⋯1-αkηakα1ηa1α2ηa2⋯αkηak1-α1ηe11-α2ηe2⋯1-αkηekα1ηe1α2ηe2⋯αkηek,W=nsi11nsi12⋯nsi1Mnsi21nsi22⋯nsi2Mndai11ndai12⋯ndai1Mndai21ndai22⋯ndai2Mndei11ndei12⋯ndei1Mndei21ndei22⋯ndei2M.Combined with the M group observations, we can get the following:(4)z=Z:=Rskak+vwhere sk≜αT;ηaT;ηeT is a set of unknown parameters, α=α1α2⋯αkT, ηa=ηa1ηa2⋯ηakT, ηe=ηe1ηe2⋯ηekT, and ·T indicates a transposing of a matrix or vector. z=Z(:), ak=X(:), and v=W(:), respectively, represent the vectors of observation matrix Z, amplitude matrix X, and noise matrix W.


## 3. Bayesian Parameter Estimation Method

### 3.1. Prior Distribution Hypothesis

Because the actual location of the target is unknown, it can be assumed that it is evenly distributed between two consecutive sampling points [13]:(5)αj~U0,1,ηaj~U-1,1,ηej~U-1,1j=1,2,⋯kBecause the independence between targets, the parameter space of sk can be expressed as (6)Φk≜0,1k×-1,1k×-1,1kUnder the condition that the number of targets is k, the location parameters of all k targets are(7)psk∣k=∏j=1k11-0·11--1·11--1=14kThe prior distribution of target number k obeys censored Poisson distribution: (8)pk∝Λkk!e-Λ,k∈0,kmaxDocument [14] points out that the number of targets that can be distinguished within the beam cannot exceed 5 at most, that is, kmax=5. Λ can be understood as the number of desired targets, and its influence can be eliminated by calculating Bias factor pz∣k1/pz∣k2. ‎Under the condition ofk,sk is known, the target amplitude ak obeys the Gauss distribution whose mean value is zero and the covariance is σk2Σ. Σ=diag(Σk,⋯,Σk︸M个‎‎), Σk-1=δ-2Ik×k. δ2 denotes the expected signal-to-noise ratio and obeys the conjugate Gama prior distribution, that is,(9)δ2~Igα02,β02The noise power σk2 can also be assumed to obey the conjugate Gama prior distribution:(10)σk2~Igυ02,γ02And α0, β0, υ0, and γ0 stand for superparameters.


### 3.2. Posteriori Distribution Derivation

According to Bayes principle, the joint posterior distribution of all parameters is(11)pk,θk∣z=pk,θk,zpz∝pz∣k,θkpk,θkwhere pz∣k,θk is a likelihood function, which can be obtained by Formula (4), that is,(12)pz∣k,θk=2πσk2-6M/2exp⁡-z-Rskak22σk2The joint posterior distribution of k,θk can be obtained by chain decomposition method, that is,(13)pk,θk=pk,sk,ak,σk2=pk,sk,ak∣σk2pσk2According to the prior distribution hypothesis, we can get the joint distribution of k,sk,ak as follows:(14)pk,sk,ak∣σk2∝Λke-Λk!·14k·exp⁡-akTΣ-1ak/2σk22πσk2Σ1/2Using the above assumptions and derivation, we can get the following:(15)pk,θk∣z∝pz∣k,θkpk,sk,ak∣σk2pσk2∝2πσk2-6M/2exp⁡-γ0+zTPkz2σk2σk2-υ0/2-1Λkk!·14k·2πσk2Σ-1/2exp⁡-ak-nkTNk-1ak-nk2σk2where Nk-1=Σ-1+RT(sk)R(sk), nk=NkRT(sk)z, and Pk=I-R(sk)NkRT(sk).According to the prior distribution hypothesis ofak and σk2, we can eliminate the ak and σk2 in Formula (15) and finally get the joint posterior distribution of k,sk as(16)pk,sk∣z∝γ0+zTPkz2-6M+υ0/2Λ/1+δ2M/2kk!·4kBecause the location parameters of each target are independent of each other, if k, ski, and z are known, then the location parameter vector sk(:,i) of the i target follows the following distribution:(17)psk:,i∣k,ski,z∝γ0+zTPkz2-6M+υ0/2where ski≜sk∖sk(:,i) represents all the parameters except columni in matrix sk, 1≤i≤k.


## 3.1. Prior Distribution Hypothesis

Because the actual location of the target is unknown, it can be assumed that it is evenly distributed between two consecutive sampling points [13]:(5)αj~U0,1,ηaj~U-1,1,ηej~U-1,1j=1,2,⋯kBecause the independence between targets, the parameter space of sk can be expressed as (6)Φk≜0,1k×-1,1k×-1,1kUnder the condition that the number of targets is k, the location parameters of all k targets are(7)psk∣k=∏j=1k11-0·11--1·11--1=14kThe prior distribution of target number k obeys censored Poisson distribution: (8)pk∝Λkk!e-Λ,k∈0,kmaxDocument [14] points out that the number of targets that can be distinguished within the beam cannot exceed 5 at most, that is, kmax=5. Λ can be understood as the number of desired targets, and its influence can be eliminated by calculating Bias factor pz∣k1/pz∣k2. ‎Under the condition ofk,sk is known, the target amplitude ak obeys the Gauss distribution whose mean value is zero and the covariance is σk2Σ. Σ=diag(Σk,⋯,Σk︸M个‎‎), Σk-1=δ-2Ik×k. δ2 denotes the expected signal-to-noise ratio and obeys the conjugate Gama prior distribution, that is,(9)δ2~Igα02,β02The noise power σk2 can also be assumed to obey the conjugate Gama prior distribution:(10)σk2~Igυ02,γ02And α0, β0, υ0, and γ0 stand for superparameters.


## 3.2. Posteriori Distribution Derivation

According to Bayes principle, the joint posterior distribution of all parameters is(11)pk,θk∣z=pk,θk,zpz∝pz∣k,θkpk,θkwhere pz∣k,θk is a likelihood function, which can be obtained by Formula (4), that is,(12)pz∣k,θk=2πσk2-6M/2exp⁡-z-Rskak22σk2The joint posterior distribution of k,θk can be obtained by chain decomposition method, that is,(13)pk,θk=pk,sk,ak,σk2=pk,sk,ak∣σk2pσk2According to the prior distribution hypothesis, we can get the joint distribution of k,sk,ak as follows:(14)pk,sk,ak∣σk2∝Λke-Λk!·14k·exp⁡-akTΣ-1ak/2σk22πσk2Σ1/2Using the above assumptions and derivation, we can get the following:(15)pk,θk∣z∝pz∣k,θkpk,sk,ak∣σk2pσk2∝2πσk2-6M/2exp⁡-γ0+zTPkz2σk2σk2-υ0/2-1Λkk!·14k·2πσk2Σ-1/2exp⁡-ak-nkTNk-1ak-nk2σk2where Nk-1=Σ-1+RT(sk)R(sk), nk=NkRT(sk)z, and Pk=I-R(sk)NkRT(sk).According to the prior distribution hypothesis ofak and σk2, we can eliminate the ak and σk2 in Formula (15) and finally get the joint posterior distribution of k,sk as(16)pk,sk∣z∝γ0+zTPkz2-6M+υ0/2Λ/1+δ2M/2kk!·4kBecause the location parameters of each target are independent of each other, if k, ski, and z are known, then the location parameter vector sk(:,i) of the i target follows the following distribution:(17)psk:,i∣k,ski,z∝γ0+zTPkz2-6M+υ0/2where ski≜sk∖sk(:,i) represents all the parameters except columni in matrix sk, 1≤i≤k.


## 4. Parameter Estimation Algorithm Based on RJ-MCMC

RJ-MCMC is essentially a generalized state space MH algorithm [15, 16]. It includes three states: the birth process, the destruction process, and the renewal process. The corresponding probabilities are bk, dk, and uk respectively, and for all k∈0,kmax, there are bk+dk+uk=1. bk and dk+1 are as follows:(18)bk≜c·min⁡1,pk+1pkdk+1≜c·min⁡1,pkpk+1where c=0.5 is a jump parameter, which controls the jump probability of the renewal process to the birth and death process.The main process of the RJ-MCMC algorithm can be described as follows [17].Step 1.
Initialization:k(0),θk(0)∈ΘStep 2.
Theith iteration process is as follows: Sampleu~U(0,1);Ifu≤bk(i), execute the birth process;Otherwise, ifu≤bk(i)+dk(i), then execute the destruction process;Otherwise, execute the renewal process.Step 3.
i ← i + 1, and return to Step 2.


### 4.1. The Birth Process

Suppose the current parameter space isk,sk and k<kmax. After randomly generating a new target, the parameter space becomes k+1,sk+1. According to the prior distribution of target parameters, we can get the recommended distribution of birth process and extinction process, that is,(19)qk+1,sk+1∣k,sk=bk·11-0·11--1·11--1=bk4(20)qk,sk∣k+1,sk+1=dk+1k+1Combining formulas (19) and (20), we can get the reception rate and the reception probability are(21)rbirth=pk+1,sk+1∣zpk,sk∣z·qk,sk∣k+1,sk+1qk+1,sk+1∣k,sk=γ0+zTPkzγ0+zTPk+1z6M+ν0/2·11+δ2M/2·1k+1Abirth=min⁡1,rbirthThe specific steps of the birth process can described as below.Step 1.
A new target is generated randomly according to the prior distribution of target location parameters.Step 2.
Calculate theAbirth by the Formula (21); Take sampleu~U(0,1).Step 3.
Ifu≤Abirth, the state of Markov chain is changed to k+1,sk+1; Otherwise, the original statek,sk is kept.


### 4.2. The Destruction Process

Suppose the current parameter space isk+1,sk+1, k>0. Then we remove one target randomly and get the new parameter space k,sk. Similar to the birth process, we can get the reception rate and the reception probabilities of the destruction process are(22)rdeath=pk-1,sk-1∣zpk,sk∣z·qk,sk∣k-1,sk-1qk-1,sk-1∣k,sk=γ0+zTPkzγ0+zTPk-1z6M+ν0/2·1+δ2M/2·kAdeath=min1,rdeathThe specific steps of the destruction process can described as below.Step 1.
Remove one target randomly from the existingk+1 goals.Step 2.
Calculate theAdeath by Formula (22); Take sampleu~U(0,1).Step 3.
Ifu≤Adeath, the state of Markov chain is changed to k,sk; Otherwise, the original statek+1,sk+1 is kept.


### 4.3. The Renewal Process

The renewal process is actually a standard MCMC sampling process. In this article, we use mixed MH sampling to generatesk, in order to achieve optimal search and traversal of target location parameter space, as shown in Table 1.Table 1
The mixed MH algorithm.
Forj=1,⋯,k Take sampleu1~U(0,1); Ifu1<λ (0<λ<1), take the proposal distribution q1sk∗:,j∣sk:,j. And the target distribution of psk:,j∣z,skj as invariant distribution ofπ(sk(:,j)), MH performed a step sampling; Else, take the proposal distributionq2sk∗:,j∣sk:,j, And the target distribution of psk:,j∣z,skj as invariant distribution ofπ(sk(:,j)), MH performed a step sampling; End if End forThe receiving probability of the update process can be obtained as(23)Aupdate=min⁡1,rupdate=min⁡1,γ0+zTPkzγ0+zTPk∗z6M+ν0/2qisk,j∣sk,j∗qisk,j∗∣sk,jTo determine the convergence of the direct determination method, calculate a mean ergodic every 50 iterations, if the absolute value of three consecutive mean calculation results the difference of the control coefficient error is less than the set, to determine the estimate convergence and to stop the iteration.


## 4.1. The Birth Process

Suppose the current parameter space isk,sk and k<kmax. After randomly generating a new target, the parameter space becomes k+1,sk+1. According to the prior distribution of target parameters, we can get the recommended distribution of birth process and extinction process, that is,(19)qk+1,sk+1∣k,sk=bk·11-0·11--1·11--1=bk4(20)qk,sk∣k+1,sk+1=dk+1k+1Combining formulas (19) and (20), we can get the reception rate and the reception probability are(21)rbirth=pk+1,sk+1∣zpk,sk∣z·qk,sk∣k+1,sk+1qk+1,sk+1∣k,sk=γ0+zTPkzγ0+zTPk+1z6M+ν0/2·11+δ2M/2·1k+1Abirth=min⁡1,rbirthThe specific steps of the birth process can described as below.Step 1.
A new target is generated randomly according to the prior distribution of target location parameters.Step 2.
Calculate theAbirth by the Formula (21); Take sampleu~U(0,1).Step 3.
Ifu≤Abirth, the state of Markov chain is changed to k+1,sk+1; Otherwise, the original statek,sk is kept.


## 4.2. The Destruction Process

Suppose the current parameter space isk+1,sk+1, k>0. Then we remove one target randomly and get the new parameter space k,sk. Similar to the birth process, we can get the reception rate and the reception probabilities of the destruction process are(22)rdeath=pk-1,sk-1∣zpk,sk∣z·qk,sk∣k-1,sk-1qk-1,sk-1∣k,sk=γ0+zTPkzγ0+zTPk-1z6M+ν0/2·1+δ2M/2·kAdeath=min1,rdeathThe specific steps of the destruction process can described as below.Step 1.
Remove one target randomly from the existingk+1 goals.Step 2.
Calculate theAdeath by Formula (22); Take sampleu~U(0,1).Step 3.
Ifu≤Adeath, the state of Markov chain is changed to k,sk; Otherwise, the original statek+1,sk+1 is kept.


## 4.3. The Renewal Process

The renewal process is actually a standard MCMC sampling process. In this article, we use mixed MH sampling to generatesk, in order to achieve optimal search and traversal of target location parameter space, as shown in Table 1.Table 1
The mixed MH algorithm.
Forj=1,⋯,k Take sampleu1~U(0,1); Ifu1<λ (0<λ<1), take the proposal distribution q1sk∗:,j∣sk:,j. And the target distribution of psk:,j∣z,skj as invariant distribution ofπ(sk(:,j)), MH performed a step sampling; Else, take the proposal distributionq2sk∗:,j∣sk:,j, And the target distribution of psk:,j∣z,skj as invariant distribution ofπ(sk(:,j)), MH performed a step sampling; End if End forThe receiving probability of the update process can be obtained as(23)Aupdate=min⁡1,rupdate=min⁡1,γ0+zTPkzγ0+zTPk∗z6M+ν0/2qisk,j∣sk,j∗qisk,j∗∣sk,jTo determine the convergence of the direct determination method, calculate a mean ergodic every 50 iterations, if the absolute value of three consecutive mean calculation results the difference of the control coefficient error is less than the set, to determine the estimate convergence and to stop the iteration.


## 5. Computer Simulation Results

Assume there is a maximum of 5 goals between two consecutive sampling assumptions. Location parameters and their settings in the document [11] are the same, as shown in Table 2. Given the hyperparameters and other related parameters as [12], α0=2, β0=10, ν0=γ0=0.1, λ=0.2, and Σq2=I/30M. Monte Carlo simulation 200 times is to get the following results. We will call the method in the document [11, 14] abbreviated as ML-MDL.Table 2
Two consecutive sampling points between the 5 goals of the position.
Target number Sub gate distance Azimuth Pitch 1 0.15 -0.8 -0.9 2 0.35 -0.6 0.7 3 0.55 0.1 -0.4 4 0.75 0.5 0.8 5 0.95 0.9 -0.8Experiment 1 (performance analysis for estimating target location parameters).
When the target has the same power, Figures4–6 show the wavelet door relative distance of target 1, azimuth, and pitch angle estimation accuracy, under different SNR conditions and different real target number conditions.
From Figures4–6, we can see that, in any case, the estimated RMSE will increase as the number of real targets increases. For ML-MDL method, increasing the number of observation pulses can improve the accuracy of parameter estimation to a certain extent. In addition, we can see that the estimation performance of the RJ-MCMC method is between the other two cases under the condition of 10 subpulses. That is to say, compared with the ML-MDL method, the RJ-MCMC method has better estimation performance under the same sub pulse number condition.Figure 4
RMSE ofα estimation under the same power, different target numbers, and different SNR conditions.Figure 5
RMSE ofηa estimation under the same power, different target numbers, and different SNR conditions.Figure 6
RMSE ofηe estimation under the same power, different target numbers, and different SNR conditions.Experiment 2 (RMSE estimation of two targets under different SNR conditions).
Suppose there are two different power targets and the SNR of target 1 changes from 15dB to 30dB and SNR of target 2 varies from 17dB to 32dB. Every 5dB results are achieved across 200 Monte-Carlo simulation, as shown in Figures7–9. In the diagram, the SNR of the abscissa represents the SNR of target 1.
As can be seen from Figures7–9, when there are two targets with different power and when the target SNR is greater than 15dB, the target’s relative distance, azimuth, and pitch angle of the target in the subwave gate can be, respectively, 0.045, 0.08, and 0.09, respectively. And obviously, the parameter estimation accuracy of the target with larger SNR is better than that of the target with smaller SNR.Figure 7
RMSE ofα estimation of target 1 and target 2 under the different power, and different SNR conditions.Figure 8
RMSE ofηa estimation of target 1 and target 2 under the different power, and different SNR conditions.Figure 9
RMSE ofηe estimation of target 1 and target 2 under the different power, and different SNR conditions.Both the proposed algorithm and the ML-MDL algorithm are the concrete implementation of MCMC algorithm. The computational complexity of the two methods corresponds, mainly depending on the number of subpulses M. However, as can be seen from Figures7–9, the estimation performance of RJ-MCMC method under 10 subpulses is between the other two cases. That is to say, compared with ML-MDL method, RJ-MCMC method has better estimation performance under the same number of subpulses. This also means that compared with ML-MDL method, RJ-MCMC method needs less observation time to obtain the same or similar estimation accuracy as ML-MDL method, which is of great significance for real-time processing.


## 6. Conclusions

In this paper, we propose a multiple parameter estimation method based on RJ-MCMC for multiple nondiscernible targets. Different from the traditional estimation methods, we need to estimate the target parameters and then estimate the order of the model (the number of the targets). The proposed method can simultaneously complete the joint estimation of the target number and the target location parameters. More importantly, the method proposed in this chapter is applicable to many situations with different power and nondistinguishable target. The simulation results show that the method proposed in this chapter requires less observation time to obtain similar and even better estimation performance than the ML-MDL method, which is of great significance for real-time processing.



---

*Source: 1504308-2019-01-06.xml*