Designing Constant Modulus Sequences with Good Correlation and Doppler Properties for Simultaneous Polarimetric Radar

Abstract

Simultaneous polarimetric radar transmits a pair of orthogonal waveforms both of which must have good auto- and cross-correlation properties. Besides, high Doppler tolerance is also required in measuring the highly maneuvering targets. A new method for the design of sequences with good correlation and Doppler properties is proposed. We formulate a fourth-order polynomial, but unconstrained, minimization problem. An iterative algorithm based on the gradient method on the phases is applied to solve it. Numerical results demonstrate the superiority of the proposed algorithm compared to the previous state-of-the-art method.

1. Introduction

In recent years, the simultaneous polarimetric scheme has been widely used to obtain accurate polarization features, which can be described by a second-order polarization scattering matrix (PSM), of targets [1,2,3,4]. A pair of waveforms transmitted in this mode is required to be stringently orthogonal, which is usually evaluated by the Isolation of the waveforms, to reduce the interference caused by simultaneous transmission and reception [4]. Meanwhile, the autocorrelation properties of the sequences, which usually represent the peak side-lobe level (PSL) and the integrated side-lobe level (ISL), are also of great importance [5]. Generally, the peak of the sidelobes corresponds to falsely detected objects, while high PSL result in masking of the weak targets with small radar cross-section (RCS) next to high signature targets in nearby range cells. Besides, for typical use of the polarimetric radar in meteorology, the ISL is an important metric. In [6,7,8], it has been pointed out that the targets for weather radar are extended volume scatterers and range sidelobes are a major source of error for quantitative applications. Therefore, designing waveforms with good correlation (in the rest of the article, correlation is used to denote both auto- and cross-correlation) properties is of great significance for simultaneous polarimetric radar.

Sequences design with good correlation properties is a traditional problem in radar and communication systems [9,10,11,12]. Focusing on the designing problem, numerous algorithms have been proposed to design the sequences. Among these techniques, a class of typical methods to design the sequences is to use the intelligent algorithms directly. Deng et al. designed unimodular sequences by using Simulated Annealing algorithm and analysed the performance of the sequences in applications of orthogonal netted radar [13]. In addition, Liu et al. used genetic algorithm to obtain sequences with good correlation properties [14]. However, the correlation properties, including PSL and Isolation of the results designed by the above-mentioned algorithms will become worse than the theoretical limiting performance with the sequences length increasing [15,16]. Besides, many researchers have also proposed computationally efficient cyclic optimization algorithms for the design of unimodular sequences. Following a similar line of derivation, Stoica et al. proposed a series of four algorithms containing CAP, CAN, WeCAN and CAD [17], which are based on the minimization of ISL with high computational efficiency. Meanwhile, these procedures can optimize the specified part of the correlation function. It should be pointed out that in some cases, the interest lies in making partial sidelobes small rather than making all sidelobes small [5,15]. However, the practical convergence rate of these four algorithms becomes slow with the sequences length increasing. Based on a customized Limited-Memory Broyden Fletcher Goldfarb and Shanno (LBFGS) algorithm, the problem of minimizing the concerned sidelobes of the correlation function is addressed by Wang et al. in [18]. By defining a new correlation matrix and using the Fast Fourier Transforms (FFT) algorithm, they moderated the computational complexity and improved the convergence rate of the method compared to the WeCAN algorithm. In a recent work [5], Palomar et al. developed the general majorization-minimization (MM) method to tackle the optimization problem arising from the sequences design. By means of the FFT algorithm, they solve the problem of designing sets of very long sequences. However, like LBFGS and MM methods, the Doppler tolerance of the sequences is not considered. In other words, they lack the ability to optimize the Doppler tolerance of the sequence sets.

In [19], Kretschmer et al. have proved that the correlation performance of the phase-coded sequences is sensitive to the target velocity. Even if the target velocity is low, the PSL and Isolation of the waveforms will seriously deteriorate compared with the same metrics of the static target’s echoes. Focusing on this problem, Pezeshki et al. investigated the Doppler tolerant waveforms designing problem for the polarimetric radar [20,21,22]. The constituent waveforms are Golay complementary codes, which have been used in many active sensing and communication systems, for instance, radar pulse compression [23], orthogonal frequency-division multiplexing (OFDM) [24] and channel estimation [25], because of the perfect autocorrelation properties along the zero Doppler shift. Based on the Prouhet-Thue-Morse sequence, they constructed Doppler resilient sequences and the range sidelobes almost vanished at modest Doppler shifts. However, the cross-correlation properties of the sequences cannot be optimized by their methods. Meanwhile, the correlation properties of the sequences keep good only in a small Doppler shifts interval around the zero Doppler shift. Their method lacks the ability to design sequences with good correlation properties in specific range and Doppler bins of the ambiguity function. In [26], Cui et al. considered the local ambiguity function shaping and proposed an accelerate iterative sequential optimization (AISO) algorithm to minimized the average value of the weighted integrated sidelobe level (WISL) over specific Doppler bins and range bins of interest. However, the orthogonality of the sequences that is of great importance for simultaneous polarimetric radar is not considered in AISO algorithm.

In this article, based on the gradient method shown in [18,27], we propose a new cyclic algorithm that can design sequences not only with good correlation properties but also with high Doppler tolerance. The obtained sequences can be used for the space surveillance radar to improve the measurement accuracy of the PSM of the highly maneuvering targets, including satellites, spacecraft, etc. The proposed method can be summarized as follows:

• Suppressing a specific part of the correlation function
• Good correlation properties under motion states of interest
• Better sidelobes both in terms of PSL and ISL

The remainder of this paper is organized as follows: the problem for designing sequences with good correlation and Doppler properties is formulated in Section 2; in Section 3, the Iterative Algorithm–Gradient (IAG) algorithm is proposed to solve the design problem; simulations are presented to validate the method in Section 4; and Section 5 concludes the paper.

$MathematicalNotations$: In this paper, it is assumed that a lower-case letter (e.g., a) denotes a scalar; a boldface lower-case letter (e.g., $\mathbf{a}$) denotes a vector; and a boldface upper-case letter (e.g., $\mathbf{A}$) indicates a matrix. Additionally, the symbols ${\mathbf{A}}^{\mathrm{T}}$ and ${\mathbf{A}}^{\mathrm{H}}$ represent transpose and conjugate transpose of the matrix $\mathbf{A}$, respectively. $\left|a\right|$ and ${\left(a\right)}^{*}$ indicate the absolute value and the conjugate of the scalar a. Besides, ${∥·∥}_F$ denotes the Frobuinous norm of a matrix and ⊙ represents the Hadamard product of two matrices with the same dimension. $Re\left(a\right)$ represents the real part of the scalar a and $\mathrm{Tr}\left(\mathbf{A}\right)$ denotes the trace of the matrix $\mathbf{A}$.

2. Problem Formulation

A pair of constant modulus sequences used for simultaneous polarimetric radar can be written as

$$\left\{\begin{array}{c}{s}_{\mathrm{H}}\left(t\right)=\frac{1}{\sqrt{N\tau }}\sum _{n=1}^N\alpha \left[t-\left(n-1\right)\tau \right]x_{\mathrm{H}}\left(n\right)\hfill \\ s_{\mathrm{V}}\left(t\right)=\frac{1}{\sqrt{N\tau }}\sum _{n=1}^N\alpha \left[t-\left(n-1\right)\tau \right]x_{\mathrm{V}}\left(n\right)\hfill \end{array}\right.$$

(1)

where

$$x_p\left(n\right)=e^{j{\varphi }_p\left(n\right)},p=\mathrm{H},\mathrm{V}\mathrm{and}n=1,2,\cdots ,N$$

(2)

are the sequences to be designed (it is assumed that the phases $\left\{{\varphi }_p\left(n\right)\right\}$ can be arbitrary values from $\left[-\pi ,\pi \right]$), H and V represent the horizontal and vertical polarized channels, $\tau$ is the time duration of one subpulse and $\alpha \left(t\right)$ is the shaping function, e.g., a rectangular pulse with amplitude 1. When the target moves, the Doppler shift of the received signals in different times is

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{d},l}\left(n\right)=2\pi f_{\mathrm{d},l}\left(n\right)n\tau =4\pi \frac{v_l+a_{l}n\tau }{\lambda }n\tau \\ \hfill n=1,2,\cdots ,N\mathrm{and}l=1,2,\cdots ,L,\end{array}$$

(3)

where $\lambda$ is the carrier wavelength, $v_l$ and $a_l$ are the radial velocity and acceleration, respectively, the subscript l represents the $l\mathrm{th}$ motion state of interest, and L is the amount of the concerning motion states. The (aperiodic) correlation of ${\left\{x_p\left(n\right)\right\}}_{n=1}^N$ and ${\left\{x_q\left(n\right)\right\}}_{n=1}^N$ (in the rest of the paper, $p,q$ are both used to denote H and V) is defined as

$$r_{pq,l}\left(k\right)=\left\{\begin{array}{c}\sum _{n=k+1}^N{x}_p\left(n\right){\left[x_q\left(n-k\right)e^{j{\varphi }_{\mathrm{d},l}\left(n-k\right)}\right]}^{*},0\le k\le N-1\hfill \\ \sum _{n=1}^{N+k}{x}_p\left(n\right){\left[x_q\left(n-k\right)e^{j{\varphi }_{\mathrm{d},l}\left(n-k\right)}\right]}^{*},-N+1\le k<0\hfill \end{array}\right.$$

(4)

When $p=q$, (4) becomes the autocorrelation of ${\left\{x_p\left(n\right)\right\}}_{n=1}^N$. Denote the transmitted sequences by an N-by-2 matrix

$$\mathbf{X}=\left[{\mathbf{x}}_{\mathrm{H}}{\mathbf{x}}_{\mathrm{V}}\right],$$

(5)

and the Doppler echoes sequences by an N-by-$2L$ matrix

$${\mathbf{X}}_{\mathrm{d}}=\left[{\mathbf{x}}_{\mathrm{H}}\odot {\mathbf{D}}_1{\mathbf{x}}_{\mathrm{V}}\odot {\mathbf{D}}_1\cdots {\mathbf{x}}_{\mathrm{H}}\odot {\mathbf{D}}_L{\mathbf{x}}_{\mathrm{V}}\odot {\mathbf{D}}_L\right]$$

(6)

where

$${\mathbf{x}}_p={\left[x_p\left(1\right),\cdots ,x_p\left(N\right)\right]}^{\mathrm{T}},p=\mathrm{H},\mathrm{V}$$

(7)

and

$${\mathbf{D}}_l={\left[e^{j{\varphi }_{\mathrm{d},l}\left(1\right)},\cdots ,e^{j{\varphi }_{\mathrm{d},l}\left(N\right)}\right]}^{\mathrm{T}},l=1,2,\cdots ,L.$$

(8)

Then, the correlation entries for the kth lag in (4) are given by the elements of the following correlation matrix [18]

$${\mathbf{R}}_k={\mathbf{X}}_{\mathrm{d}}^{\mathrm{H}}{\mathbf{J}}_k\mathbf{X}$$

(9)

where

$${\mathbf{J}}_k={\left[\begin{array}{c}\mathbf{0}{\mathbf{I}}_{N-k}\hfill \\ \mathbf{0}\mathbf{0}\hfill \end{array}\right]}_{N\times N}={\mathbf{J}}_{-k}^{\mathrm{H}}$$

(10)

and $\mathbf{I}$ is the identity matrix. It should be noted that the elements of ${\mathbf{R}}_k$ are made up of the autocorrelations and cross-correlations of the sequences. A compact optimization model for deigning constant modulus sequences with low correlation sidelobes levels can be formulated as

$$\begin{array}{c}\underset{{\mathbf{x}}_{\mathrm{H}},{\mathbf{x}}_{\mathrm{V}}}{min}\sum _{k=-G}^G{w}_k{∥{\mathbf{X}}_{\mathrm{d}}^{\mathrm{H}}{\mathbf{J}}_k\mathbf{X}-N\delta \left(k\right)\mathbf{Y}∥}_F^2\hfill \\ \mathrm{subject}\mathrm{to}\left|x_p\left(n\right)\right|=1,n=1,2,\cdots ,N.\hfill \end{array}$$

(11)

where

$$\mathbf{Y}={\left[{\mathbf{I}}_2{\mathbf{I}}_2\cdots {\mathbf{I}}_2\right]}_{2L\times 2}^{\mathrm{T}},$$

(12)

and $\delta \left(k\right)$ is the impulse function, i.e., $\delta \left(k\right)=1$ for $k=0$ and otherwise $\delta \left(k\right)=0$. In (11), ${\left\{w_k\right\}}_{k=-G}^G$ are positive real weights chosen by the user, and $G<N$ is also a preset positive integer, which represents the lag interval of interest. The model (11) involves minimizing a fourth-order polynomial with nonlinear equality constraints, which is numerically difficult to handle. However, using the phases ${\left\{{\varphi }_p\left(n\right)\right\}}_{n=1}^N$ as new variables, the constant modulus constraints can be dropped and (11) can be formulated as an unconstrained fourth-order polynomial minimization problem

$$\underset{\mathsf{\varphi }}{min}f\left(\mathsf{\varphi }\right)=\sum _{k=-G}^G{w}_k{∥{\mathbf{X}}_{\mathrm{d}}^{\mathrm{H}}\left(\mathsf{\varphi }\right){\mathbf{J}}_k\mathbf{X}\left(\mathsf{\varphi }\right)-N\delta \left(k\right)\mathbf{Y}∥}_F^2$$

(13)

where

$$\mathsf{\varphi }={\left[{\varphi }_{\mathrm{H}}\left(1\right)\cdots {\varphi }_{\mathrm{H}}\left(N\right){\varphi }_{\mathrm{V}}\left(1\right)\cdots {\varphi }_{\mathrm{V}}\left(N\right)\right]}^{\mathrm{T}}.$$

(14)

Compared with (11), obtaining the global optimal solution of (13) is still a NP-hard problem. However, considering the local optimal solutions, it is easier to handle the optimization problem. In this paper, an iterative algorithm based on the gradient method, which can be applied directly, is proposed to solve the problem of (13). In the gradient method, the consuming time of each iteration is determined by the calculation of $f\left(\mathsf{\varphi }\right)$ and $\nabla f\left(\mathsf{\varphi }\right)$. The computational costs of calculating $f\left(\mathsf{\varphi }\right)$ and $\nabla f\left(\mathsf{\varphi }\right)$ directly are $O\left(NGL\right)$ and $O\left(N^{2}GL\right)$, respectively, according to (13) [18]. In the following, an efficient method is proposed to compute them.

3. Solving the Model by IAG Algorithm

Introducing the notation, the matrix $\mathbf{M}$ is defined as

$$\mathbf{M}=\left[{\mathbf{m}}_{-N+1}\cdots \underset{{\mathbf{M}}_{-G}}{\underbrace{{\mathbf{m}}_{-G}\cdots {\mathbf{m}}_{-1}}}{\mathbf{m}}_0\underset{{\mathbf{M}}_G}{\underbrace{{\mathbf{m}}_1\cdots {\mathbf{m}}_G}}\cdots {\mathbf{m}}_{N-1}\right]$$

(15)

where ${\mathbf{m}}_k=\sqrt{w_k}\mathrm{vec}\left({\mathbf{X}}_{\mathrm{d}}^{\mathrm{H}}\left(\varphi \right){\mathbf{J}}_k\mathbf{X}\left(\varphi \right)\right)$, and $\mathrm{vec}\left(·\right)$ vectorizes a matrix by stacking its columns on top of one another. Then the function $f\left(\mathsf{\varphi }\right)$ can be rewritten as

$$f\left(\mathsf{\varphi }\right)={∥{\mathbf{m}}_0-\mathbf{h}∥}_2^2+{∥{\mathbf{M}}_{-G}∥}_F^2+{∥{\mathbf{M}}_G∥}_F^2$$

(16)

where $\mathbf{h}=\sqrt{w_0}\mathrm{vec}\left(N\mathbf{Y}\right)$ and ${∥·∥}_2$ denotes the 2-norm of a vector. According to (16), the gradient $\nabla f\left(\mathsf{\varphi }\right)$ can be computed as follows:

$$\frac{\partial f\left(\mathsf{\varphi }\right)}{\partial {\varphi }_p\left(n\right)}=2Re\left[{\left({\mathbf{m}}_0-\mathbf{h}\right)}^{\mathrm{H}}\frac{\partial {\mathbf{m}}_0}{\partial {\varphi }_p\left(n\right)}+\mathrm{Tr}\left({\mathbf{M}}_{-G}^{\mathrm{H}}\frac{\partial {\mathbf{M}}_{-G}}{\partial {\varphi }_p\left(n\right)}+{\mathbf{M}}_G^{\mathrm{H}}\frac{\partial {\mathbf{M}}_G}{\partial {\varphi }_p\left(n\right)}\right)\right]$$

(17)

Furthermore, it should be noted that each row vector of $\mathbf{M}$ can be computed by the convolution product, ${\mathbf{x}}_p\otimes {\left({\left({\mathbf{x}}_q\odot {\mathbf{D}}_l\right)}^{\mathrm{r}}\right)}^{*}$, where the sequence ${\mathbf{x}}_q^{\mathrm{r}}$ can be obtained by reversing the order of the entries of ${\mathbf{x}}_q$. In other words, the entries in ${\mathbf{m}}_0$, ${\mathbf{M}}_{-G}$ and ${\mathbf{M}}_G$ can be computed by truncating ${\mathbf{x}}_p\otimes {\left({\left({\mathbf{x}}_q\odot {\mathbf{D}}_l\right)}^{\mathrm{r}}\right)}^{*}$ from the indices $-G$ to G. Besides, the convolution operation in the time domain corresponds to the product operation in the frequency domain. Thus, the FFT algorithm can be applied to obtain ${\mathbf{m}}_0$, ${\mathbf{M}}_{-G}$ and ${\mathbf{M}}_G$ by $12L\left(2N-1\right){log}_2\left(2N-1\right)$ complex multiplication operations at most. Computing ${∥{\mathbf{m}}_0-\mathbf{h}∥}_2^2+{∥{\mathbf{M}}_G∥}_F^2+{∥{\mathbf{M}}_{-G}∥}_F^2$ takes $4L\left(2G+1\right)$ complex multiplication operations. So the computation cost of obtaining $f\left(\mathsf{\varphi }\right)$ is $O\left(LG+LN{log}_{2}N\right)$. As for $\nabla f\left(\mathsf{\varphi }\right)$, there are $2N$ entries. According to (4), for $k\ge 0$ (the same as $k<0$), it should be noticed that

$$\begin{array}{cc}\hfill r_{pq,l}\left(k\right)& =\sum _{m=k+1,m\ne n}^N{x}_p\left(m\right){\left[x_q\left(m-k\right)e^{j{\varphi }_{\mathrm{d},l}\left(m-k\right)}\right]}^{*}\hfill \\ \hfill & +x_p\left(n\right){\left[x_q\left(n-k\right)e^{j{\varphi }_{\mathrm{d},l}\left(n-k\right)}\right]}^{*}\hfill \end{array},$$

(18)

$$\begin{array}{cc}\hfill \frac{\partial r_{pq,l}\left(k\right)}{\partial {\varphi }_p\left(n\right)}& ={\left[x_q\left(n-k\right)e^{j{\varphi }_{\mathrm{d},l}\left(n-k\right)}\right]}^{*}\frac{\partial x_p\left(n\right)}{\partial {\varphi }_p\left(n\right)}\hfill \\ \hfill & ={\left[x_q\left(n-k\right)e^{j{\varphi }_{\mathrm{d},l}\left(n-k\right)}\right]}^{*}j{x}_p\left(n\right)\hfill \end{array},$$

(19)

and

$$\begin{array}{c}{r}_{pp,l}\left(k\right)=\sum _{\begin{array}{c}m=k+1\\ m\ne n,m\ne n+k\end{array}}^N{x}_p\left(m\right){\left[x_p\left(m-k\right)e^{j{\varphi }_{\mathrm{d},l}\left(m-k\right)}\right]}^{*}\hfill \\ +x_p\left(n\right){\left[x_p\left(n-k\right)e^{j{\varphi }_{\mathrm{d},l}\left(n-k\right)}\right]}^{*}+x_p\left(n+k\right){\left[x_p\left(n\right)e^{j{\varphi }_{\mathrm{d},l}\left(n\right)}\right]}^{*}\hfill \end{array},$$

(20)

$$\begin{array}{cc}\hfill \frac{\partial r_{pp,l}\left(k\right)}{\partial {\varphi }_p\left(n\right)}& ={\left[x_p\left(n-k\right)e^{j{\varphi }_{\mathrm{d},l}\left(n-k\right)}\right]}^{*}j{x}_p\left(n\right)\hfill \\ \hfill & -x_p\left(n+k\right)e^{-j{\varphi }_{\mathrm{d},l}\left(n\right)}j{x}_p^{*}\left(n\right)\hfill \end{array}.$$

(21)

The entries in ${\mathbf{m}}_0^{\mathrm{H}}\frac{\partial {\mathbf{m}}_0}{\partial {\varphi }_p\left(n\right)}$, ${\mathbf{M}}_{-G}^{\mathrm{H}}\frac{\partial {\mathbf{M}}_{-G}}{\partial {\varphi }_p\left(n\right)}$ and ${\mathbf{M}}_G^{\mathrm{H}}\frac{\partial {\mathbf{M}}_G}{\partial {\varphi }_p\left(n\right)}$ can be obtained according to (19) and (21), and correspondingly, it takes $2L\left(2G+1\right)+1$ complex multiplication operations. Thus, the computational cost of $\nabla f\left(\mathsf{\varphi }\right)$ is $O\left(NGL\right)$. The computational complexities to compute $f\left(\mathsf{\varphi }\right)$ and $\nabla f\left(\mathsf{\varphi }\right)$ using the original expression (13) and new expression (16) are shown in

Table 1. The computational complexity using original expression (13) and new expression (16).

	Original (13)	New (16)
$f\left(\mathsf{\varphi }\right)$	$O\left(NGL\right)$	$O\left(LG+LN{log}_{2}N\right)$
$\nabla f\left(\mathsf{\varphi }\right)$	$O\left(N^{2}GL\right)$	$O\left(NGL\right)$

. Furthermore, using the gradient $\nabla f\left(\mathsf{\varphi }\right)$, the following IAG algorithm shown in

Table 2. Steps for the IAG algorithm.

Step 0:  Set X to initial sequences, and fix the motion states of interest, which means the $v_l$ and $a_l$ should be determined. Step 1:  Compute the gradient $\nabla f\left(\mathsf{\varphi }\right)$ according to (17), (19) and (21). Step 2:  Renew the phases of the sequences using $\mathrm{vec}\left({\mathbf{X}}^{i+1}\right)=\mathrm{vec}\left({\mathbf{X}}^i\right)·e^{-j{\mu }^i\nabla f\left(\mathsf{\varphi }\right)}$, where the step length ${\mu }^i$ can be obtained according to the line search algorithm and the index i represents the $i\mathrm{th}$ iteration [28]. Step 3:  Repeat Step 1 and Step 2 until a termination criterion is satisfied, e.g., $\left|f{\left(\mathsf{\varphi }\right)}^{i+1}-f{\left(\mathsf{\varphi }\right)}^i\right|\le \epsilon$, where $\epsilon$ is a predefined threshold.

can be performed to obtain the sequences with good correlation properties and high Doppler tolerance. One thing should be pointed out is that for the classical gradient method, the line search rule is based on the Wolfe conditions which can ensure the stability of the iterations. In this paper, to improve the efficiency of the algorithm, a low computation complexity line search rule, named Armijo-rule, is used [28].

4. Numerical Results

Here, we provide numerical examples to illustrate the performance of the proposed IAG algorithm. The state-of-the-art WeCAN algorithm proposed in [17] is considered for comparison. The metrics Peak Side-lobe Level ($PSL$), Integrated Side-lobe Level ($ISL$) and the Isolation (I) are used to evaluate the performance of the sequences which are defined as follows [4]

$$PSL=20{log}_{10}\frac{max\left(\left|r_{pp,l}\left(k\right)\right|\right)}{\left|r_{pp,l}\left(0\right)\right|},k\in \left[-G,-1\right]\cup \left[1,G\right],$$

(22)

$$ISL=\frac{2{\displaystyle \sum _{k=1}^G}\left({\left|r_{HH,l}\left(k\right)\right|}^2+{\left|r_{VV,l}\left(k\right)\right|}^2\right)+{\displaystyle \sum _{k=-G}^G}\left({\left|r_{HV,l}\left(k\right)\right|}^2+{\left|r_{VH,l}\left(k\right)\right|}^2\right)}{2{\left|r_{HH,l}\left(0\right)\right|}^2}$$

(23)

$$I=20{log}_{10}\frac{max\left(\left|r_{pq,l}\left(k\right)\right|\right)}{\left|r_{pp,l}\left(0\right)\right|},p\ne q\mathrm{and}k\in \left[-G,G\right].$$

(24)

The running time is obtained using the Matlab 2016b version, running on a standard PC (with a 2.6 GHz Core i7 CPU and 4-GB RAM).

In the first numerical example, we consider the measurement of the stationary target, meaning the Doppler shift in (3) is zero. Simulation parameters are shown in

Table 3. Simulation parameters for the first example.

N	G	${\mathit{w}}_{\mathit{k}}$	$\mathit{\epsilon }$
256	39	${\left\{w_k\right\}}_{k=-G}^G=1$	${10}^{-14}$

. Assume the sequence length is $N=256$ and the lag interval of interest is $G=39$. The weight coefficients are ${\left\{w_k\right\}}_{k=-G}^G=1$ for the IAG method. For the WeCAN algorithm, the parameter $w_0$ should be large enough to guarantee the stability of the algorithm. Thus, the weights are set as 30 while $k=0$ and 1 while $k\ne 0$ for the WeCAN method. The predefined threshold of the predefined threshold is $\epsilon ={10}^{-14}$ for the two methods. To avoid unreasonable solution, the minimum iteration number is set to 1000 in the simulations.

Figure 1a,b show the convergence curves of the two algorithms. The convergence rate of the $PSL$ and I of the IAG method is much faster than that of the WeCAN method. The $PSL$ and I of the sequences designed by the IAG method achieve about $-325\mathrm{dB}$ after 700 iterations, however, the corresponding number of iterations is over ${10}^5$ for the WeCAN method. Further, Figure 2a,b illustrate the correlation level of the sequences designed by the two methods. The $PSL$ and I of the IAG sequences are about $-325\mathrm{dB}$ better than that of the WeCAN sequences with $-265\mathrm{dB}$ under this situation. Meanwhile, the $ISL$ is $8.92\times {10}^{-32}$ for IAG and $8.92\times {10}^{-32}$ for WeCAN, which indicates the energy of the sidelobes is lower for IAG sequences. In other words, the IAG sequences are suitable for observing the volume targets that extend over range areas. The consuming time and iteration number are shown in the

Table 4. Comparison between IAG and WeCAN.

	Iteration Number	Total Execution Time (s)	Execution Time per Iteration (s)
IAG	1000	196.7	0.194
WeCAN	112,432	25,881.5	0.231

. It can be observed that the average execution time of per iteration of the IAG method is shorter than that of the WeCAN method. As mentioned before, by the new expression (16) and the using of FFT, the consuming time for computing the objective function and the corresponding gradient becomes less, leading the average execution time of the IAG method shorter. Meanwhile, the IAG method takes much less iterations than the WeCAN method. Thus, the total execution time is less than one in a hundred in comparison with that of the WeCAN method.

In the second numerical example, designing sequences for measuring the polarization features of the highly maneuvering target is simulated. In [29], Li has pointed out that typical highly maneuvering targets contain satellites, spacecrafts, etc. The velocity of these targets is usually about Mach 10 or even higher. Meanwhile, the acceleration is about $10g$ (g is the acceleration of gravity). Thus, in the simulations, the velocity and the acceleration of the target are assumed to be within $v\in \left(2000,3000\right)\mathrm{m}/\mathrm{s}$ and $a\in (50,200)\mathrm{m}/{\mathrm{s}}^2$, respectively. We uniformly discretize the velocity and acceleration intervals into L bins with the grid size $\mathsf{\Delta }v=50\mathrm{m}/\mathrm{s}$ and $\mathsf{\Delta }a=10\mathrm{m}/{\mathrm{s}}^2$. Then the motion states of interest in (3) can be expressed as the column of the following matrix $\mathbf{S}$.

$$\begin{array}{cc}\hfill \mathbf{S}& ={\left[\begin{array}{c}{v}_1{v}_2\cdots v_L\hfill \\ a_1{a}_2\cdots a_L\hfill \end{array}\right]}_{2\times L}\hfill \\ \hfill & ={\left[\begin{array}{c}20002050\cdots 3000\cdots \cdots 20002050\cdots 3000\hfill \\ 5050\cdots 50\cdots \cdots 200200\cdots 200\hfill \end{array}\right]}_{2\times L}\hfill \end{array}$$

(25)

where the amount of concerning motion states L can be computed as

$$L=\left(\frac{3000-2000}{50}+1\right)·\left(\frac{200-50}{10}+1\right)=336.$$

(26)

Other simulation parameters are shown in

Table 5. Simulation parameters for the second example.

N	G	${\mathit{w}}_{\mathit{k}}$	$\mathit{\epsilon }$	${\mathit{f}}_{\mathbf{0}}$	$\mathit{\tau }$
256	39	${\left\{w_k\right\}}_{k=-G}^G=1$	${10}^{-14}$	$10\mathrm{GHz}$	$5\times {10}^{-9}\mathrm{s}$

. The first four parameters are same as those in Table 3. Besides, taking typical simultaneous polarimetric radars, including the MERIC radar [30] and the CSU-CHILL radar [31], as references, the carrier frequency is set to $f_0=10\mathrm{GHz}$. The carrier wavelength in (3) can be computed by $\lambda =c/f_0$, where c is the speed of light. The time duration of the subpulse is $\tau =5\times {10}^{-9}\mathrm{s}$. Then the IAG algorithm shown in Table 2 can be performed to obtain the sequences with good correlation properties under the motion states of interest. Figure 3 and Figure 4 show the correlation properties of the designed sequences designed by the IAG and WeCAN algorithms through the metrics $PSL$ and I, respectively. It clearly demonstrates that the correlation properties of the sequences designed by the IAG method are better than that of the sequences designed by the WeCAN method under the given motion states. Correspondingly, the $PSL$ and I of the IAG sequences are under $-53\mathrm{dB}$, however, the same metrics of the WeCAN sequences are about $-25\mathrm{dB}$ in this situation. The $ISL$ as a function of velocity and acceleration is shown in Figure 5. It can be observed that under the same conditions, the $ISL$ of the IAG sequences is about one in a thousand in comparison with that of the WeCAN sequences. The reason is that the influence of the velocity and the acceleration of the target on the correlation properties of the sequences, which is not considered in the WeCAN method, is taken into account in the objective function (11) of the IAG method. Besides, it can be observed that the metrics, including $PSL$, $ISL$ and I, change slightly with the acceleration increasing. Since the time duration of the subpulse $\tau$ is set to $5\times {10}^{-9}\mathrm{s}$ in the simulation, the time length of the sequence is equal to $N\tau$, which is quite small in this situation, leading the slight change of the velocity. Thus, the properties of the sequences change slightly with the increasing of the acceleration.

5. Conclusions

In this paper, we propose an iterative algorithm based on the gradient method to design the constant modulus sequences with good correlation and Doppler properties for simultaneous polarimetric radar. By transforming the objective function to a new expression, the computation complexities of the gradient and the objective function are reduced by using the FFT algorithm. Compared with the state-of-the-art WeCAN approach, the proposed IAG method has a better performance in terms of correlation sidelobes levels, Doppler tolerance and execution time.