Patent ID: 12235344

DETAILED DESCRIPTION

The technical scheme of the present invention is further described below with reference to the drawings.

The present invention is applicable to pulse radar system adopting linear frequency modulation waveform, and a target detection processing process is shown asFIG.1and comprises:

Step 1) performing conventional pulse compression and coherent integration on received baseband signal to obtain a range-Doppler map, and performing pre-detection based on the range-Doppler map to obtain interested PDTs.

Since the signal model in spatial domain is similar to the slow time domain, the present invention considers the fast time domain and the slow time domain. Assuming that a number of range cells in one coherent processing interval is L, a number of pulses is K, and the received signal can be represented as a matrix Y∈L×K. Neglecting the range and Doppler migration, the elements thereof can be represented as:

yl,k=∑p=0P-1αp⁢rect⁡(l⁢Ts-2⁢rp/cTpul)⁢exp⁢(j⁢πμ⁡(l⁢Ts-2⁢rp/c)2+j⁢2⁢π⁢fd,p⁢k⁢TI)+nl,k,(1)

wherein yl,kis an lthrange cell and a kthpulse echo signal, 0≤k≤K−1, and 0≤l≤L−1. P is a number of the target, Tsis a sampling time interval, αpis a pthtarget amplitude, rpis a target range, Tpulis a pulse duration, TIis a pulse interval, μ is a frequency modulation slope, and c is light speed. fd,p=2vp/λ is Doppler shift, up is a target speed, and λ is a carrier wavelength. nl,kis additive white Gaussian noise. rect(.) represents a rectangular function.

Pulse compression is performed in fast time domain and coherent integration is performed in slow time domain to obtain a range-Doppler map. Then, based on the range-Doppler map, a conventional CFAR detection method is adopted for target pre-detection to obtain interested PDTs, wherein the corresponding ranges and Doppler frequencies of cells, wherein PDTs are pre-detected, are represented by rζand fζ, respectively.

Step 2) estimating the PDT range and the Doppler parameter obtained by pre-detection to obtain estimated valuesand;

assuming that the number of the PDTs is I, and the corresponding range rζand Doppler frequency fζcan be calculated according to the range cell and the frequency cell wherein the I PDTs are located, respectively.

The PDTs comprise the main components of the received signal, thus, based on I, rζand fζ, the received signal can be approximately represented as:
y≈Aζβ+n,(2)
wherein Y represents the received signal of one coherent processing interval, and is obtained by vectorizing Y through stacking the columns into a vector; β is a dimension-reduction target vector, an ithelement βithereof represents the true complex amplitude of an ithPDT; and if the ithPDT is a false alarm, then βi=0. n is an additive white Gaussian noise vector. Aζ∈LK×Iis an approximate dimension-reduction observation matrix, and an ithcolumn thereof is represented as:
aζ,i=sd(fζ,i)⊗si(rζ,i),  (3)
wherein the symbol ⊗ represents Kronecker product, sd(fζ,i) represents the Doppler domain steering vector corresponding to the PDT with a Doppler frequency of fζ,i, and si(rζ,i) represents the fast-time domain steering vector corresponding to the PDT with a range of rζ,i, sd(fζ,i) and si(rζ,i) are represented respectively as:
sd(fζ,i)=[1, . . . ,exp(j2πfζ,ikTI), . . . ,exp(j2πfζ,i(K−1)TI)]T,  (4)
si(rζ,i)=[qi(0), . . . ,qi(l), . . . ,qi(L−1)]T(5)
wherein the superscript T represents transpose;
wherein,

qi(1)=rect⁡(lTs-2⁢rϛ,i/cTpul)⁢exp⁢(j⁢πμ⁡(l⁢Ts-2⁢rϛ,i/c)2).(6)

Generally, β is still sparse, and reconstruction of β can be performed based on sparse recovery; assuming that the reconstruction result is represented as {circumflex over (β)}, then target detection can be realized based on {circumflex over (β)}. However, in practical applications, the PDTs are not located in integer cells, that is, rζand fζdeviate from the true PDT values, and the reconstruction directly based on the equation (2) is faced with the off-grid problem. In this regard, the present invention firstly estimates the true parametersandof the PDTs based on rζand fζ, and then performs the target reconstruction and detection.

Assuming that the true observation matrix corresponding to the PDTs is, ignoring the effect of missed detection targets during pre-detection, the received signal can be represented as:
y=β+n(7)
wherein=[], the symbol; is used for connecting two vectors to form a vector. An ithcolumn of Ais=sd() ⊗.is unknown, however,is obviously very close to Aζ, that is, the true target parameters are close to parameters corresponding to the integer cells wherein the PDT are located. Let=[rζ: fζ], rζand fζare known. Sinceis very close to Aζ, thenis also very close to θζ, i.e., ∥−θζ∥ is very small. Therefore, in the present invention, it is considered to estimatebased on θζ. Then, the observation matrixcan be obtained. The specific procedures are as follows:

based on a maximum likelihood criterion, estimates ofand β are given as:

{θˆϱ,β^}=argminθ,βy-A⁡(θ)⁢β22,(8)

wherein θ=[r; f], r and f represent the range and the Doppler frequency, respectively; the minimum over β is attained for:
{circumflex over (β)}=(A(θ)HA(θ))−1A(θ)Hy(9)
then the cost function in the equation (8) can be further represented as:
g(θ)=∥y−A(θ)(A(θ)HA(θ))−1A(θ)Hy∥22(10)

Obviously, a minimum value of the equation (10) is obtained at θ=. A first-order Taylor approximation is performed on a first-order derivative of g(θ) atto obtain:
∇θg(θ)≈∇θg()+∇θ2g()(θ−),  (11)

Obviously, ∇θg()=0, then,
≈θ−(∇θ2g())−1(∇θg(θ)).  (12)

Becauseis also very close to θζ, θ is substituted with θζ, andin the Hessian matrix (∇θ2g()) is substituted with θζto obtain the estimate of:
≈θζ−(∇θ2g(θζ))−1(∇θg(θζ)).  (13)

In the equation (13), the first-order derivative vector and the second-order derivative matrix of g(θ) are required to be calculated, and the calculation is very complicated; therefore, a simplified solution method is further given below.

The following equation can be obtained from the equation (7):
Wiy=Wiβ+Win(14)
wherein Wi=diag(wi)=diag(wd⊗wl,i), wdrepresents a normalized window function in the slow time domain, wl,irepresents a normalized window function in the fast time domain of the ithPDT, wl,i=[qw,0, . . . , qw,l, . . . , qw,L−1]T, and

qw,l=rcct⁢(l⁢Ts-2⁢rϛ,i/cTpul)⁢wc(l⁢Ts-2⁢rϛ,i/c).(15)
wc(t) is a continuous form of the window function in a time Tpul. In the case of normalization, obviously, ∥wi∥22=∥Wi∥F=1, wherein the subscript F represents the Frobenius norm of the matrix. Then, based on the least squares criterion, the estimates ofand β is:

{θˆϱ,βˆ}=argminθ,βWi⁢y-Wi⁢A⁡(θ)⁢β22.(16)

The equation (14) actually weights the data in the fast and slow time domains, and the weighting can be used to decouple the ithPDT from other PDTs in the echo, that is, the ithPDT and other PDTs have approximately no mutual influence. Then, g(θ) can be further approximately represented as:
g(θ)=∥Wiy−βiWiai(ηi)∥22+Wiy−WiA\i(θ)β\i∥22−(Wiy)HWiy(17)
wherein ηi=[ri, fi]T, riand firepresent the range and the Doppler frequency of the ithPDT, respectively; airepresents the steering vector of the ithPDT which is calculated by the equation (3), A\idenotes the matrix obtained from A by deleting the ithcolumn, and β\idenotes the vector obtained from β by deleting the ithentry, to minimizing g(θ), ui(ηi)=∥ Wiy−βiWiai(ηi)∥22should attain its minimum. Then, the estimate of the ithPDT, denoted by=[]T, can be obtained by minimizing ui(ηi).

Minimizing ui(ηi) the estimate of βiis given by:
{circumflex over (β)}i=aiH(ηi)Wi2y.(18)

Inserting equation (18) into ui(ηi), then, minimizing ui(ηi) is equivalent to minimizing the following equation:
zi(ηi)=(aiH(ηi)Wi2y)(yHWi2ai(ηi))  (19)

Referring to equation (13), the estimate ofis given by:

(r^ϱ,if^ϱ,i)≈(rϛ,ifϛ,i)-(∂2zi∂ri200∂2zi∂fi2)(ηi=ηϛ,i)-1⁢(∂zi∂ri∂zi∂fi)(ηi=ηϛ,i),(20)
wherein ηζ,i=[rζ,i, fζ,i]T. In order to obtain higher estimation accuracy, the estimate is iteratively updated using the following equation:
=−(∇η2zi())−1∇ηzi(),  (21)
wherein=nζ, i, and t denotes the tthiteration usually, t=2 can meet the actual demand. Parameter estimation is performed for each PDT, and thencan be obtained.

Step 3) establishing a dimension-reduction observation model of a received signal based on the estimated valuesand;

can be obtained based on(i.e.,and), and the received signal, based on the, may be further represented as:
y≈β+n.(22)

Usually, the number I of the PDTs is much smaller than a number of the cells corresponding to the range-Doppler map; therefore, the equation (22) can greatly reduce the dimension of the vector to be reconstructed, and the equation (22) is the dimension-reduction observation model of the present invention.

Step 4) reconstructing a target vector based on the dimension-reduction observation model and the sparse recovery algorithm;

The present invention reconstructs β based on the generalized approximate message Passing (GAMP) algorithm and the equation (22). Assuming that β follows the i.i.d Bernoulli-Gaussian distribution, the marginal probability density function thereof is:
ρ(βi)=(1−p)δ(βi)+ρ(βi;κ,τq)  (23)
wherein δ is Dirac function, ρ denotes the fraction of nonzero components. κ and τqrepresent the mean and variance of the Gaussian components, respectively. ρ, κ and τqare all unknown, and expectation-maximum (EM) algorithm can incorporated to iteratively learn them.

GAMP can find both a sparse estimate and a non-sparse, noisy estimate of β, denoted by {circumflex over (β)} andβ, respectively. It can be proved that the modulus of the noise in theβapproximately follows The Gaussian distribution, and the simulation result is shown inFIG.3AtoFIG.3C.

Step 5) designing a generalized likelihood ratio detector based on the reconstruction result for target detection and outputting detection results and their parameters.

The present invention is based onβfor target detection. Determining whether the ithPDT is a target or not can be summarized as the following hypothesis tests:
H0:y=\iβ\i+n
H1:y=+\iβ\i+n.(24)

H0and H1are two hypotheses in the hypothesis tests. H0represents that the ithPDT is not the target, and H1represents that the ithPDT is the target.is the true steering vector of the ithPDT, and is also calculated by the equation (3).

Based on the generalized likelihood ratio criterion, the detector is given by:

❘"\[LeftBracketingBar]"β~i❘"\[RightBracketingBar]"⁢H1><H0⁢γ⁢σ~,(25)

wherein {tilde over (σ)} is the variance of the noise modulus in {tilde over (β)}, and γ is the detection threshold; since the noise modulus follows approximately the Gaussian distribution, the equation (25) is a CFAR detector.

The implementation steps of the radar target detection method based on EBD provided by the present invention are described above, and the effectiveness of the method is verified by both and measured tests. The method of the present invention will be referred to as EBD below.FIG.2AtoFIG.2Cshow the parameter performance of EBD in the present invention. Simulation parameters are shown in Table 1, 4 targets are added to the received signal and three following cases are considered: case 1: the target is located in the integer Doppler cells, and the range cells are 100, 100.2, 140.3, and 191.5 (the fractional part represents the value of deviation from the integer cells); case 2: the target is located in the integer range cells, and the frequency cells are 10, 15.2, 31.3 and 42.5; case 3: the Doppler and range cells of targets are all off-grid, and are set as (15.2,105.3), (31.2,140.4), (42.3,191.1) and (53.5,270.5), respectively; the target parameter estimation accuracy in all three cases is measured by mean square error (MSE), and the results are shown inFIG.2AtoFIG.2C. The results show that the method of estimating the parameters before detection can effectively estimate the range and Doppler of the target.

TABLE 1Radar simulation parametersNos.ParametersValues1Pulse interval0.2 ms2Pulse width25 μs3Bandwidth4 MHz4Sampling rate5 MHz5Number of pulses56Number of FFT points647Window used in the fastHammingtime domain8Window used in the slowChebyshev, −45 dBtime domain

FIG.3AtoFIG.3Cshow the statistical characteristics of the noise modulus in the reconstruction result {tilde over (β)}. The following 2 cases are considered: case 1: the received signal contains no target; case 2: 10 targets are added to the received signal, and the SNRs are 0 dB. When the noise characteristics are counted, the true target samples are removed, that is, only the noise samples are counted.FIG.3AandFIG.3Bshow the quantile-quantile plot (Q-Q plot) curves of the noise samples in both cases, respectively, and are compared with the standard Gaussian distribution.FIG.3Cshows the correlation coefficient of the noise samples in both cases. The results show that the noise approximately follows independent Gaussian distribution.

InFIG.4AandFIG.4B, the performance of the EBD provided by the present invention is compared with3other detectors including the conventional detection method (TSPM, windowed), the ideal matched filtering (IMF, not windowed), the method of off-grid sparse recovery (OGSR) based on the equation (2). In the simulation, the parameters shown in Table 1 are still adopted, 10 targets are added in each simulation, the target range and Doppler are randomly generated. Simulation results show that EBD is approximately a CFAR detector within the entire interval of SNR, and the gain of detection performance is about 1.9 dB compared with the conventional detection method (TSPM) under the condition that the false alarm rate is 10−5.

FIG.5A.FIG.5BandFIG.6verify the method of the present invention based on the measured data. The radar parameters corresponding to the adopted data are as follows: the bandwidth is 8 MHz, the sampling rate is 10 MHz, the pulse width is 12 microseconds, the number of single frame pulses is 128, the number of FFT points is 128, and the repetition frequency is 12.5 KHz. In the test, two unmanned aerial vehicles DJI Phantom 3 are used as cooperative targets, 1 frame is selected for analysis, and the range-Doppler map is shown inFIG.5A. The selected frame contains no other target. Because the SNRs of the two unmanned aerial vehicles are high and are 26.3 dB and 21.3 dB, respectively; in order to verify the performance of the provided algorithm, the white Gaussian noise is added to the original received signal. For example, after increasing the noise level by 6 dB, the range-Doppler map using the conventional processing is as shown inFIG.5B.

The performance comparison between the EBD and TSPM is shown inFIG.6, wherein the abscissa represents the noise increment and is represented by Δσn2. It can be seen from the results that the EBD has a better noise resistance ability, i.e., better target detection performance, than TSPM. For both cooperative targets, the EBD has performance gains of 2.2 dB and 1.5 dB, respectively, at a false alarm rate of 10−5.