Algorithms for Designing Unimodular Sequences with High Doppler Tolerance for Simultaneous Fully Polarimetric Radar

Simultaneous fully polarimetric radar uses orthogonal polarization channels to transmit a pair of signals, both of which must have good auto- and cross-correlation characteristics. In this paper, the design of sequences with good correlation properties and Doppler tolerance is investigated. New cyclic algorithms, namely, Cyclic Algorithm-Gradient I (CAGI) and Cyclic Algorithm-Gradient II (CAGII) are proposed to solve the optimization problem. Meanwhile, the sequences designed in this paper have ultra-low auto- and cross-correlation side-lobes in a specified lag interval. Numerical experiments are conducted to demonstrate and validate the superiority of the proposed cyclic algorithms, especially for the measurement of moving targets.


Introduction
In the past decades, polarimetric features, which can be described by a second order backscattering matrix (BSM) of the target, have been widely used in various fields, such as terrain observation, disaster surveillance and atmospheric remote sensing [1][2][3][4]. In order to accurately obtain the BSM, two fully polarimetric measurement schemes, called alternate polarimetric measurement (APM) and simultaneous polarimetric measurement (SPM) schemes have been widely investigated since the 1980s [5][6][7][8].
In alternate mode, the transmitted polarization states are switched alternately between horizontal (H) and vertical (V) polarizations while both polarizations are received simultaneously on reception. The estimations of BSM for targets with high velocity, such as satellites and spacecrafts, are inaccurate in alternate mode due to the change of their polarimetric features during the switch on transmission phase. In contrast, in simultaneous mode, the two orthogonal polarization states are transmitted and received simultaneously. Thus, the BSM of the targets can be retrieved within one pulse recurrent time (PRT). In this case, the limitation caused due to the change of the transmitted polarization states can be overcome. However, the transmitted signals must be orthogonal in order to reduce the interference caused by simultaneous transmission and reception, which is usually evaluated by the cross-correlation properties of the waveforms [9]. Numerous researchers have focused on this problem and attempted to design waveforms with good orthogonal properties. In [10], Babur et al. proposed a method for designing the quasi-orthogonal Frequency Modulated Continuous Waveforms (FMCW) and analyzed its application in multiple-input multiple-output (MIMO) radar. In addition, Liu [11] exploited the genetic algorithm to obtain the coding sequences with good correlation properties for MIMO radar applications. However, the above-mentioned algorithms have not considered the effect of the Doppler shift on the properties of the waveforms.
Designing waveforms with good autocorrelation properties is also of great importance in applications of polarimetric radar. In general, the peak of the side-lobes is related to the probability of false alarm. At the same time, if the peak of the side-lobes is high, the return of the weak targets will be masked by that of the targets with large Radar Cross Section (RCS) in nearby range cells. Therefore, the peak of the side-lobes should be as low as possible in order to improve the detection performance of weak targets. In general, there exist two major merits to measure the autocorrelation property of the waveforms that are Integrated Side-lobes Level (ISL) and Peak Side-lobes Level (PSL). For the minimization of PSL, Deng [12] designed unimodular sequences (in the rest of the article, sequences are used to denote polyphase sequences) by using the simulated annealing algorithm and analyzed their performance in applications of orthogonal netted radar. In [13], Song et al. discussed the problem of minimizing the l p (2 ≤ p < ∞) norm on the side-lobes of autocorrelation function. They approached the minimization of PSL by increasing the value of p. However, their algorithm named "monotonic minimizer (MM) for l p -metric" lacks the ability to suppress a specified part of the autocorrelation function. It should be pointed out that, in some radar applications, the maximum difference between the arrival time of the interest sequences is smaller than the duration of the emitted sequences [13][14][15][16]. Hence, only the partial side-lobes of the correlation function are needed to be low. In addition, Esmaeili-Najafabadi et al. [14] proposed a series of cyclic algorithms (CA), including PSL Minimization Quadratic Approach (PMQA), PSL minimization Algorithm where the smallest Rectangular (PMAR), PSL Optimization Cyclic Algorithm (POCA) and Randomized PSL Optimization Cyclic Algorithm (RPOCA), by using the Chebyshev distance and l ∞ norm. They addressed the problem of optimizing the specified part of the autocorrelation function. Meanwhile, a set of waveforms with low cross-correlation side-lobes can be obtained by applying chaotic waveforms as the initial sequences in their algorithms. However, these four algorithms lack the ability to suppress the specified side-lobes for the cross-correlation function. For the minimization of ISL, P. Stoica and H. He et al. [15,16] proposed a series of cyclic algorithms for unimodular sequences, including CA-pruned (CAP), CA-new (CAN), weighted-CAN (WeCAN) and CA-direct (CAD). In addition to autocorrelation, they minimized the cross-correlation between the generated sequences and optimized the specified part of the correlation (in the rest of the article, correlation is used to denote both auto-and cross-correlation) function. However, the practical convergence rate of these algorithms becomes slow with the increase in the length of the sequences and they lack the ability to suppress all the side-lobes of the correlation function.
Furthermore, all of the above-mentioned algorithms have been proposed for the applications of static or low-velocity targets. For moving targets, Doppler loss occurs at the matched filters on reception [17]. It has been pointed out in [18] that the phase-coded sequences are quite sensitive to the Doppler shift. In other words, even if the velocity of the target is low, the output of the matched filters will still significantly deteriorate compared with the processing results of the static target's echoes. In a recent work, an efficient gradient algorithm has been proposed to optimize the autocorrelation property of the sequences in different Doppler shift [19]. However, this algorithm cannot optimize the cross-correlation properties of the sequences, which is important for the measurement of simultaneous polarimetric radar as mentioned before. Cui et al. [20] proposed an accelerated iterative sequential optimization (AISO) algorithm that can minimize the ISL of the autocorrelation function in different Doppler shift of interest and reduce the computations compared with the algorithm in [19]. However, the AISO algorithm also does not consider the optimization of the cross-correlation properties of the sequences. In addition, Doppler tolerant complementary code sets are also developed these days because of their potential of making all the autocorrelation side-lobes sum to zero, at least in theory, but the orthogonality of the complementary sequences are not considered in most articles [21,22]. In addition, the Doppler frequency in [19,20] is considered to be a constant in one PRT, which means that the target is assumed to be in a uniform rectilinear motion state. However, the instantaneous acceleration of maneuvering targets, such as missile and spacecraft, can reach 10g (g is the gravitational acceleration). At this moment, the Doppler frequency is not a constant in one PRT, and the Doppler phase model should be modified to match the movement of the targets.
The above-mentioned three metrics are usually regarded as the measures of the sequences used in fully polarimetric radars. The low cross-correlation side-lobes improve the estimation accuracy of BSM, the low autocorrelation side-lobes increase the detection performance of weak targets and the high Doppler tolerance enhances the adaptability of the waveforms to the movement of targets. In most studies, this design problem has been usually reformulated as an optimization of a complex matrix [12][13][14][15][16][17]23]. The requirements of Doppler tolerance, autocorrelation and cross-correlation properties make the optimization problem difficult. This paper extends the approaches in [19] and proposes new cyclic algorithms for designing sequences with good correlation properties and Doppler tolerance. The rest of this paper is organized as follows: Section 2 presents the problem formulation; in Section 3, the Cyclic Algorithm-Gradient I (CAGI) and Cyclic Algorithm-Gradient II (CAGII) algorithms are introduced to design the sequences with good properties; the numerical examples are provided in Section 4 to verify the performance of the proposed algorithms, followed by the conclusions in Section 5.
Notation: Boldface uppercase and lowercase letters are used to represent matrices and vectors, respectively. See Table 1 for other notations used in this paper.

Simultaneous Polarimetric Measurement and Doppler Analysis
The simplified signal processing flow chart of the SPM is depicted in Figure 1. Let s H (t) and s V (t) represent the simultaneously transmitted signals via two orthogonal polarization channels. To facilitate the discussion, the transmitted signals can be given in vector form as:  For a point target, the received signals are the Doppler-shifted and time-delayed version of the transmitted signals. Thus, they can be expressed as [24]: where r H (t) and r V (t) represent the received signals from the two orthogonal polarization channels, respectively, τ is the round-trip propagation delay and f d is the Doppler frequency determined by the target velocity and the carrier wavelength. R R R and T T T represent the effect of channel, antenna and propagation path on the BSM during the reception and transmission, respectively, and S S S is the BSM of the target. Processed by the matched filters shown in Figure 1, the output can be expressed as: where The subscripts p, q denote the H and V polarization channels. h H (t) and h V (t) represent the matched filters of H and V polarization channels, respectively. When the transmitted signals satisfy the orthogonal condition, i.e., then, if R R R and T T T are known, all the parameters of the BSM can be estimated accurately. The sequences used for simultaneous polarimetric measurement can be written as: where x p (n) = e jφ p (n) , p = H, V and n = 1, · · · , N are the sequences to be designed, φ p (n) can be an arbitrary value between [−π, π] (in the rest of the article, p and q are used to denote H and V), N is the length of the sequences, τ 0 is the time duration of the sub-pulse and α(t) is the envelope function with unit amplitude. The (aperiodic) cross-correlation function of x p (n) N n=1 and x q (n) N n=1 at lag k is defined as: When p = q, Equation (9) becomes the autocorrelation function of x p (n) N n=1 . It should be pointed out that, in the real radar system, each code interval has more than one point dependent on a frequency of discretization. This leads to the output of the matched filters to be different from the discrete correlation function r pq (k) defined in this paper. However, the authors have proved in [25] that the discrete correlation function r pq (k) has almost the same properties as the continuous correlation r pq (t), where t is the continuous delay. Therefore, it is reasonable to optimize the discrete correlation function instead of the continuous correlation function to obtain the sequences with good correlation properties. Moreover, the correlation functions of the sequences need to be modified because, when the target moves, the Doppler loss occurs at the matched filters on the receiving end. Here, the radial velocity and the acceleration of the target are assumed to be v 0 and a, respectively. Under the condition of uniformly accelerated rectilinear motion, the Doppler frequency of the received signals in different times can be written as: where λ 0 is the carrier wavelength. Based on Equation (10), the Doppler phase shift of the received signals in different times can be written as: Then, the maximum Doppler phase shift can be obtained as φ max = 2π f d (N)Nτ 0 when n is equal to N, and the Doppler frequency is an arithmetic sequence, which means f d (n + 1) − f d (n) is a constant. Based on Equations (10) and (11), Equation (9) needs to be modified due to the existence of the Doppler phase shift as follows:

Cyclic Algorithm for Designing Sequences under One Motion State
The sequences with good correlation properties under a certain motion state can be designed by minimizing the following criterion: The first two terms of Equation (13) represent the energy of the side-lobes of the autocorrelation functions and the last term is the energy of the cross-correlation functions. Obviously, reducing the energy terms will improve the correlation properties of the sequences. According to Equation (12), another thing that should be noticed is: This means that, if the motion state of the target is given, r pp (0) is a constant that is not related to the designed sequences x p (n) N n=1 . Thus, the following criterion can also be optimized to design the sequences with good correlation properties: Before solving the optimization problem of Equation (15), it should be emphasized that, in some radar applications, the maximum difference between the arrival times of the sequences of interest is smaller than the duration of the emitted sequences [13][14][15][16]26]. This means that only the partial side-lobes of the correlation function need to be small, instead of making all the side-lobes small. Therefore, a more proper minimization criterion than Equation (15) is given by: where P is selected based on the prior information about the applications. In this paper, the gradient descent method is used to solve the optimization problem of Equation (16). In order to facilitate the following discussion, the transmitted sequences are denoted by: Since x p is a complex sequence with unit envelope, the gradient can be obtained by differentiating the Equation (16) with respect to its phase. Using the chain rule of the derivative, the derivative with respect to the phase of x p that is written as ∂ψ ∂φ p (n) can be obtained by calculating the derivatives with respect to the real and the imaginary parts of x p : More specifically, next, ∂ψ ∂φ H (n) is calculated and ∂ψ ∂φ V (n) can be obtained by the same steps. It can be observed that where As |b| 2 = Re (b) 2 + Im (b) 2 where b is a complex number, it can be derived that: + Re(r HV (k)) ∂ Re(r HV (k)) ∂(·) + Im(r HV (k)) ∂ Im(r HV (k)) ∂(·) The details for calculating the partial derivatives with respect to Re(x H (n)) and Im(x H (n)) are provided in Appendix A. Here, the results are given by Equations (22) and (23): By substituting Equations (22) and (23) into the chain rule Equation (18), it can be derived that: The ∂ψ ∂φ V (n) can be obtained following the same procedures from Equations (19)- (24). Thus, generally, the gradient of ψ with respect to φ p (n), which is the phase of the sequence, can be written as: In order to facilitate the following discussion, r r r pq = r pq (n)  (25) can be rewritten as: by defining: Equation (26) can be expressed in a vector form: It is well known that the convolution operation in the time domain corresponds to the product operation in the frequency domain. If F (·) denotes the Fast Fourier Transforms (FFT) operation, it can be derived that the convolution of a M-dimensional sequence b b b and a N-dimensional sequence d d d is where and Therefore, Equation (29) provides an efficient computation of the gradient Equation (28), which is expressed as a sum of convolution operations. Then, by using the gradient ∇ p , the following algorithm shown in Table 2, named the Cyclic Algorithm-Gradient I (CAGI), can be performed to obtain the sequences with good correlation properties under one certain motion state. Table 2. Steps for the CAGI Algorithm.
Step 0: Set {x H (n)} N n=1 and {x V (n)} N n=1 to initial sequences (e.g., x p (n) N n=1 can be set to e jφ p (n) N n=1 , where φ p (n) N n=1 are independent random variables distributed in [0, 2π]). Fix the motion state, which means v 0 and a should be given, and compute the vector D D D by Equation (11).
Step 1: Calculate the gradient ∇ p according to Equation (28) Step 2: Renew the phases of the sequences using x i+1 p = x i p · e −jβ i ∇ p , where the step length β i is computed according to the line search algorithm [27].
Step 3: Repeat Steps 1 and 2 until a stop criterion is satisfied, e.g., ψ i+1 − ψ i ≤ ε, where ψ i is the objective function at the ith iteration and ε is a predefined threshold.

Cyclic Algorithm for Designing Sequences under a Set of Motion States
In this subsection, a new cyclic algorithm is developed to design sequences with good correlation properties for a set of motion states. Here, matrix L L L is used to describe the set of motion states and is defined as: Obviously, each row of the matrix L L L represents one motion state with a velocity and acceleration of v 0,l and a l , respectively, where 1 ≤ l ≤ L. Then, the criterion for designing sequences with good correlation properties for a set of motion states can be written as: and Equation (33) represents the energy of the side-lobes under different motion states of interest. Similarly, the sequences with good correlation properties under the motion states of interest can be obtained by reducing the energy of the side-lobes. According to the above discussion, it is obvious that the gradient ofψ with respect to the phase of the sequences is: where ∂ψ l ∂φ p (n) can be calculated by Equations (18)- (28). Apparently, the gradient ofψ consists on a sum of the gradient of ψ under different motion states. Therefore, Equation (36) can also be efficiently calculated by using the FFT. By defining: the following algorithm shown in Table 3, named the Cyclic Algorithm-Gradient II (CAGII), can be performed to obtain the sequences with good correlation properties under a set of motion states of interest. Table 3. Steps for the CAGII Algorithm.
Step 0: Set {x H (n)} N n=1 and {x V (n)} N n=1 to initial sequences (e.g., x p (n) N n=1 can be set to e jφ p (n) N n=1 , where φ p (n) N n=1 are independent random variables distributed in [0, 2π]). Fix the set of motion states, which means the matrix L L L should be given.
Step 2: Renew the phases of the sequences using x i+1 p = x i p · e −jβ i∇ p , where the step length β i is computed according to the line search algorithm [27].
Step 3: Repeat Steps 1 and 2 until a stop criterion is satisfied, e.g., ψ i+1 −ψ i ≤ ε, whereψ i is the objective function at the ith iteration and ε is a predefined threshold.

Suppressing Specified Parts of the Correlation Functions without Doppler Shift
The performance of the proposed CAGI algorithm and the cyclic algorithm WeCAN in [16] are assessed under a scenario where the target of interest is stationary, which means v 0 and a in Equation (10) are zero. In particular, consider the example given as follows: All of the following simulations are implemented through Matlab 2016b (MathWorks, Natick, State of Massachusetts, the United States) on an i5 3.2 GHz machine (Intel Corporation, Santa Clara, State of California, the United States) with a 4 GB RAM. The comparison between the normalized correlation of CAGI and WeCAN is presented in Figure 2 for parameter ε = 10 −6 . It can be seen from the figure that both algorithms can suppress the side-lobes to almost zero in a specified lag interval. The PCL and I of the sequences designed by the two algorithms are about −100 dB, which is low enough for the applications of simultaneous polarimetric radar [28]. The main difference between the two algorithms is the consumed time. The WeCAN algorithm consumed over 1 h while the proposed CAGI algorithm only consumed 78.12 s to generate this figure. The reason is the efficient computation of the gradient by using the FFT for CAGI algorithm.

Suppressing Specified Parts of the Correlation Functions for One Motion State
The performance of the CAGI and the WeCAN algorithms are again assessed under the second scenario, where suppressing a specified part of the correlation functions under a certain motion state is considered. Here, the highly maneuvering targets, such as satellites and spacecrafts, are of concern and the simulation parameters are set in Table 4:  The normalized correlation functions of the sequences designed by the two algorithms for parameter ε = 10 −6 are illustrated in Figures 3 and 4. Obviously, the CAGI sequences provide much lower PCL and I compared with the sequences designed by the WeCAN algorithm. A fundamental reason is that the velocity and the acceleration of the target are considered in the objective function of the CAGI. In addition, the execution time, the iteration number, PCL and I (here, PCL and I are the average values of the sequences x H and x V ) are shown in Table 5. It can be easily observed from the table that the execution time per iteration of CAGI algorithm is shorter than that of the WeCAN algorithm. The reason is that the use of the FFT improves the computation efficiency of the gradient. At the same time, the CAGI algorithm takes much less iterations than that of WeCAN algorithm. Thus, the total execution time is about five percent in comparison with that of the WeCAN algorithm.     Another drawback of the WeCAN is that it cannot easily cope with large length of sequences (i.e., the lengths more than N = 1000). In these situations, the CAGI algorithm can be applied. Table 6 provides the simulation parameters for the CAGI to optimize the sequences with length N = 5000. The normalized correlation for the sequences with parameter ε = 10 −6 are shown in Figure 5. The PCL and I of the sequences are about −93 dB, which are almost equal to that of the sequences with short length shown in Figure 3.

Suppressing Specified Parts of the Correlation Functions for a Set of Motion States
In this subsection, another situation is considered, where suppressing a specified part of the correlation functions under a set of motion states is of concern. Similarly, highly maneuvering targets, such as satellites and spacecrafts, are considered. Assume that the velocity and the acceleration of the target are within v ∈ (1000, 2000) m/s and a ∈ (75, 150) m/s 2 , respectively. In order to facilitate the calculation, the velocity and the acceleration intervals are uniformly discretized with the grid size ∆v = 10 m/s and ∆a = 5 m/s 2 , respectively. Then, the matrix L L L that describes the motion states of interest can be written as: L L L = 1000 1010 · · · 2000 · · · · · · 1000 1010 · · · 2000, 75 75 · · · 75 · · · · · · 150 150 · · · 150, T Moreover, λ 0 = 0.03 m, τ 0 = 5 × 10 −9 s, N = 256 and P = 30 are set for the simulations. Utilizing the proposed CAGII algorithm, the sequences with good correlation properties under a set of motion states can be obtained. Figure 6a-d illustrate the correlation properties of the sequences designed by the CAGII algorithm with parameter ε = 10 −6 . It can be seen clearly that within the given velocity and acceleration interval, the PCL and I of the sequences are under −50 dB, and the fluctuation of the correlation properties is within 5 dB. The reason is that the object function Equation (33) of the CAGII is the weighted sum of correlation functions with different motion states, and the weighting coefficients are all equal to 1. It can also be observed that the correlation properties of the sequences change slightly in this situation with the increasing of the acceleration. This behavior is attributed to the fact that the time length of the sequence, which is equal to the product of N and τ 0 , is short. Thus, the change of the velocity is unapparent in the duration of the sequences. However, when the sequence length is increased or the bandwidth of the sequence is reduced, meaning the duration of the sequence becomes longer, the change of the velocity will be obvious within one PRT of the transmitted sequences, and, accordingly, the Doppler frequency between the sub-pulse of the sequences will be different significantly. Furthermore, in order to verify the performance of the sequences shown in Figure 6, the auto-ambiguity functions (AF) and the cross-ambiguity functions (CF) for the sequences x H and x V with acceleration a = 100 m/s 2 are shown in Figures 7-9, respectively. It can be seen that the PCL and I in the area of interest that refers to the area with k ∈ [−29, 29] and v 0 ∈ (1000, 2000) m/s of the ambiguity function, is lower than −50 dB. In addition, it should be noticed from Figure 9 that the fluctuation of the peak of the normalized autocorrelation function caused by the mismatch of the Doppler loss is within 3 dB when the velocity of the target changes from v = 0 m/s to v = 5000 m/s. This phenomenon further confirms that the sequences designed by the CAGII algorithm have high Doppler tolerance under the given motion states.

Conclusions
In this paper, two new cyclic algorithms, called Cyclic Algorithm-Gradient I (CAGI) and Cyclic Algorithm-Gradient II (CAGII), are proposed for designing a pair of unimodular sequences for simultaneous polarimetric radar. The CAGI approach can design sequences with both good auto-and cross-correlation properties under one motion state. The CAGII is an extension of the CAGI and can solve the optimization problem under a set of motion states. The application of the FFT in the CAGI and the CAGII algorithms enable both approaches to obtain quite long sequences with good correlation properties. Several numerical simulations and comparative analysis are conducted to demonstrate and validate the superior performance of the sequences designed using the proposed methods compared with those obtained using other algorithms in the literature. Furthermore, the authors plan to conduct research on the extension of the proposed algorithms to the finite alphabet phase case and improve the Doppler tolerance and orthogonality of the complementary coded sequences in future.