## 1.

## Introduction

Inspired by the achievement in multiple-input multiple-output (MIMO) communication,^{1} the MIMO concept was introduced into the radar field, which has then rapidly drawn more and more attention.^{2}3.4.5.6.7.8.9.10.11.12.13.14.15.16.17.18.19.20.21.22.23.24.25.26.27.^{–}^{28} Compared with phased array, MIMO radar has the capability of transmitting arbitrary waveforms, which is regarded as waveform diversity.^{2} In terms of the spacing between its antennas, MIMO radar can be classified into two categories shown as (1) widely separated antennas (e.g., ^{3}) and (2) colocated (e.g., ^{2}). The former employs the widely spaced transmit and receive elements along with diverse transmitted waveforms to view the different aspects of the target thereby improving the target detection performance. In contrast, the latter employs the close-spaced elements in transmit and receive arrays to obtain the identical target radar cross sections observed from all transmit/receive paths, which can increase the spatial resolution for MIMO radar by exploiting the waveform diversity. Extensive researches have illustrated that, by designing the transmitted waveforms elaborately, MIMO radar can significantly outperform its phased radar counterpart such as improved spatial resolution, parameter identifiability,^{4}^{,}^{5} and more flexible design of the transmitting beampattern.^{6}7.8.9.10.11.12.13.14.15.16.17.18.19.20.21.22.^{–}^{23}

As a consequence, waveform design becomes one of the particularly critical issues in MIMO radar, which has received vibrant attentions recently.^{6}7.8.9.10.11.12.13.14.15.16.17.18.19.20.21.22.^{–}^{23} According to the objects needed to be optimized, the existing design methods can be classified as the following two categories: (1) only the transmitter to be designed^{6}7.8.9.10.11.12.13.^{–}^{14} and (2) the transmitter and receiver to be designed jointly.^{15}16.17.18.19.20.21.22.^{–}^{23}

Specific issues that have been considered in first category include the transmit beam pattern design and radar ambiguity function design. In Refs. 6 and 7, the waveform covariance matrix (WCM) was designed to attain a desired beam pattern, while the constant modulus signal design is considered in Refs. 89.–10. Meanwhile, the spatial-range-Dopper domain characteristics of the transmitted waveforms were considered jointly in Ref. 11 to improve the radar performance. In Ref. 12, the transmitted waveforms are designed to acquire the improvement in the detection probability for MIMO space–time adaptive processing (STAP) with perfect target and clutter prior knowledge. However, these parameters must be estimated with error in application. Therefore, the robust waveform optimization for MIMO–STAP with imperfect space–time array steering vector prior knowledge is considered for improving the worst-case (i.e., the point at which the worst value can be obtained over the convex uncertainty set) detection probability in Ref. 13. Besides, focusing on improving the worst-case estimation accuracy in the absence of clutter, a robust waveform design is proposed in Ref. 14.

Some works have been done in the second category approaches to investigate the waveform design problem by jointly optimizing the radar transmitter and receiver. In Ref. 15, the transmit waveforms were optimized for multiple point targets based on several design criteria, e.g., minimizing the trace of the Cramér–Rao Bound (CRB) matrix. The output signal-to-interference-plus-noise ratio (SINR) was maximized in Ref. 16 to acquire the better MIMO radar detection probability in case of extended target by exploiting a gradient-based method. Unfortunately, the method proposed in Ref. 16 cannot guarantee nondecreasing SINR in each step. For the purpose of guarantee convergence, Ref. 17 proposed a new iterative method. The mutual information between the echo and the target radar signatures was employed to design the transmit waveform.^{18}19.^{–}^{20} In Ref. 21, MIMO waveform was devised by minimizing the estimation error of the minimum mean-squared error (MMSE) estimators for uncorrelated and correlated targets. The joint design of the transmit WCM and receive weights is considered for improving the parameter estimation performance in Refs. 22 and 23.

Based on CRB, Li et al.^{15} have studied the waveform design issue for narrowband MIMO radar to acquire the better estimation accuracy, while Wang et al.^{22} have addressed the issue of joint design of the transmitting waveforms and biased estimator in the presence of clutter for the case of known characteristics of targets of interest on the basis of constrained CRB. Obviously, solving the design issues in these literatures must rely on the specification of parameters, e.g., angle of arrival, angle of departure, and reflection coefficients. However, the acquired knowledge of these parameters may be inaccurate due to the fact that they have to be estimated in practice, thus there exists uncertainty as to them. Consequently, the resultant estimation accuracy may be sensitive to these uncertainties in parameters, which has been depicted by numerical examples in Refs. 15 and 22.

The considerations mentioned above raises the problem of robust MIMO radar waveform design in the presence of prior uncertainty as to the target characteristics, which has been investigated in Ref. 27 in the absence of clutter, and studied in Ref. 28 in the presence of clutter. However, it is well-known that the received data are inevitably contaminated by the clutter in many applications (see, e.g., Refs. 8 and 9), and this problem was tackled in Ref. 28 based on the assumption proposed in Ref. 22 that is reasonable only under some certain conditions (see Ref. 22 for more details). Focusing on these issues, partially on the basis of our previous works,^{28} following the min–max approach, the problem of robust WCM optimization in the context of clutter is considered in this paper aiming to systematically ease the sensitivity of parameter estimation performance to the uncertainty in the initial parameter estimates by including the parameter estimation error model into the design issue, which is based on minimizing the trace of the worst-case CRB matrix. With constraints about the transmitted waveforms and uncertainties in the initial estimates, the formulated robust design problem is an NP-hard issue.^{29} To tackle this problem, a diagonal loading (DL)-based^{30} iterative approach is developed to design the WCM for improving the worst-case estimation accuracy, in which the inner optimization problem can first be decomposed to some independent subproblems by using the Hadamard’s inequality,^{31} then these subproblems can be reformulated into convex issues by using DL method, as well as the outer optimization problem also can be relaxed to a convex issue by translating the nonlinear function into a linear one, and, hence, both of them can be solved very effectively by using the CVX toolbox that is MATLAB-based modeling system for convex optimization.^{32} In what following, an optimal solution to the original issue can be obtained via the least-squares (LS) fitting of the solution acquired by the iterative algorithm.

This paper is organized as follows. First, Sec. 2 presents the system model and illustrates the WCM design issue including the convex uncertainty set of parameter estimation errors. Section 3 proposes a DL-based iteration approach and constructs an optimal solution of the original issue. Section 4 shows the efficiency of the proposed approach via numerical simulations. Finally, conclusions are summarized in Sec. 5.

Notation: Vectors and matrices are indicated by boldface lowercase and uppercase letters, respectively. Transpose, complex conjugate, and conjugate transpose are, respectively, denoted by ${\{\xb7\}}^{T}$, ${\{\xb7\}}^{*}$, and ${\{\xb7\}}^{H}$. $\mathrm{Re}\{\xb7\}$ and $\mathrm{Im}\{\xb7\}$ represent the real and image part of a complex-valued argument, respectively. Denote $E\{\xb7\}$ as the expectation operator, $\otimes $ as the Kronecker product, * as the Khatri–Rao product, $\odot $ as the Hadamard product, $\mathbf{I}$ as the identity matrix, ${\Vert \mathbf{A}\Vert}_{F}$ as the Frobenius norm of $\mathbf{A}$, $\mathrm{diag}\{\mathbf{a}\}$ as the diagonal matrix formed via its arguments, and $\mathrm{tr}\{\xb7\}$ as the trace of a square matrix. Given $\mathbf{f}:{\mathbb{R}}^{n}\to {\mathbb{R}}^{k}$, $\frac{\partial \mathbf{f}}{\partial \mathbf{\theta}}$ indicates the $k\times n$ matrix with $ij$’th element being $\partial {\mathbf{f}}_{i}/\partial {\mathbf{\theta}}_{j}$. $\mathbf{M}\u2aaf\mathbf{N}$ indicates $\mathbf{N}-\mathbf{M}$ is positive semidefinite.

## 2.

## Problem Statement

Consider a colocated narrow band MIMO radar equipped with ${M}_{t}$ transmit and ${M}_{r}$ receive elements with arbitrarily spaced distance, where the $i$’th transmit element emits the discrete-time baseband signal ${\mathbf{\phi}}_{i}\in {\mathbb{C}}^{L\times 1}$, $i=\mathrm{1,2},\dots ,{M}_{t}$ with $L$ denoting the number of snapshots. Collecting all the transmitting discrete-time baseband signals, the transmitted waveform matrix can be denoted by $\mathbf{\Phi}={[{\mathbf{\phi}}_{1},{\mathbf{\phi}}_{2},\dots ,{\mathbf{\phi}}_{{M}_{t}}]}^{T}$. If the propagation is nondispersive and the clutter returns in the target range bin can be modeled as a superimposition of a number of independent clutter patches,^{33} the data collected by the receive array can be written as (see, e.g., Refs. 2021.–22)

## (1)

$$\mathbf{Y}=\sum _{k=1}^{K}{\beta}_{k}\mathbf{a}({\theta}_{k}){\mathbf{b}}^{T}({\theta}_{k})\mathbf{\Phi}+{\mathbf{H}}_{c}\mathbf{\Phi}+\mathbf{W},$$Based on the MIMO radar model shown in Eq. (1), the CRB of $\mathbf{\theta}={[{\theta}_{1},{\theta}_{2},\dots ,{\theta}_{K}]}^{T}$ and ${\{{\beta}_{k}\}}_{k=1}^{K}$ can be written as (see Ref. 22 for the detail of the derivation)

## (2)

$$\mathbf{C}=\frac{1}{2}{\left[\begin{array}{ccc}\mathrm{Re}({\mathbf{F}}_{11})& \mathrm{Re}({\mathbf{F}}_{12})& -\mathrm{Im}({\mathbf{F}}_{12})\\ {\mathrm{Re}}^{T}({\mathbf{F}}_{12})& \mathrm{Re}({\mathbf{F}}_{22})& -\mathrm{Im}({\mathbf{F}}_{22})\\ -{\mathrm{Im}}^{T}({\mathbf{F}}_{12})& -{\mathrm{Im}}^{T}({\mathbf{F}}_{22})& \mathrm{Re}({\mathbf{F}}_{22})\end{array}\right]}^{-1},$$## (3)

$${[{\mathbf{F}}_{11}]}_{ij}={\beta}_{i}^{*}{\beta}_{j}{\dot{\mathbf{h}}}_{i}^{H}\{{[\mathbf{I}+({\mathbf{R}}_{\mathrm{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\dot{\mathbf{h}}}_{j}$$## (4)

$${[{\mathbf{F}}_{12}]}_{ij}={\beta}_{i}^{*}{\dot{\mathbf{h}}}_{i}^{H}\{{[\mathbf{I}+({\mathbf{R}}_{\mathrm{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\mathbf{h}}_{j},$$## (5)

$${[{\mathbf{F}}_{22}]}_{ij}={\mathbf{h}}_{i}^{H}\{{[\mathbf{I}+({\mathbf{R}}_{\mathrm{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\mathbf{h}}_{j},$$Observing Eqs. (2)–(5), it is noted that the CRB is a function with respect to (w.r.t.) $\mathbf{\theta}$, ${\{{\beta}_{k}\}}_{k=1}^{K}$, ${\mathbf{H}}_{c}$, and $\mathbf{W}$. In practice, the knowledge of these parameters may be inaccurate due to the fact that they must be estimated with errors. Therefore, the designed waveforms on the basis of the CRB exploiting a certain estimate may provide a rather low accuracy in case of another reasonable one (see the numerical results in Refs. 15 and 22).

In this work, similar to the model developed in Refs. 3435.–36, it can be assumed that the derivation of the $k$’th target channel matrix is modeled as follows:

where ${\tilde{\mathbf{h}}}_{k}$ and ${\mathbf{h}}_{k}$ indicate, respectively, the actual and corresponding presumed the $k$’th target channel vector, and ${\mathbf{\delta}}_{k}$ is the error of ${\tilde{\mathbf{h}}}_{k}$, which lies in the uncertainty set ${\mathcal{U}}_{1}=\{{\mathbf{\delta}}_{k}|{\Vert {\mathbf{\delta}}_{k}\Vert}_{F}\le {\zeta}_{k},k=\mathrm{1,2},\dots ,K\}$, and in which ${\dot{\tilde{\mathbf{h}}}}_{k}$ and ${\dot{\mathbf{h}}}_{k}$ stand for, respectively, the actual and corresponding presumed derivation of ${\mathbf{h}}_{k}$, and ${\dot{\mathbf{\delta}}}_{k}$ is the error of ${\dot{\tilde{\mathbf{h}}}}_{k}$, which lies in the uncertainty set ${\mathcal{U}}_{2}=\{{\dot{\mathbf{\delta}}}_{k}|{\Vert {\dot{\mathbf{\delta}}}_{k}\Vert}_{F}\le {\sigma}_{k},k=\mathrm{1,2},\dots ,K\}$.## Remark:

The uncertainty of ${\tilde{\mathbf{h}}}_{k}$ can be determined via many methods shown in Ref. 36. Due to the fact that ${\dot{\tilde{\mathbf{h}}}}_{k}$ relies on the configuration of MIMO radar, once the system geometry is fixed, similar to that of ${\tilde{\mathbf{h}}}_{k}$, the bound of ${\mathcal{U}}_{2}$ can also be obtained. Consequently, the bounds ${\zeta}_{k}$ and ${\sigma}_{k}$ can be regarded as a prior knowledge. Moreover, it has to be pointed out that ${\zeta}_{k}$ should satisfy ${\zeta}_{k}<{\Vert {\tilde{\mathbf{h}}}_{k}\Vert}_{F}$ and ${\sigma}_{k}<{\Vert {\dot{\tilde{\mathbf{h}}}}_{k}\Vert}_{F}$, i.e., the bound of the uncertainties of ${\mathbf{\delta}}_{k}$ and ${\dot{\mathbf{\delta}}}_{k}$ cannot be too large. Otherwise, ${\tilde{\mathbf{h}}}_{k}$ and ${\dot{\tilde{\mathbf{h}}}}_{k}$ may be zero and then the robust design studied here will make no sense. Moreover, similar to the discussion in Ref. 36, the specification of ${\zeta}_{k}$ and ${\sigma}_{k}$ will affect the performance of the proposed method, which should be illustrated via numerical examples in Sec. 4.

Based on the discussion above, the improvement of the worst-case estimation accuracy can be implemented by minimizing the worst-case CRB over all parameter estimates lying in ${\mathcal{U}}_{1}$ and ${\mathcal{U}}_{2}$ with the transmitted power constraint $\mathrm{tr}({\mathbf{R}}_{\mathbf{\Phi}})=\mathrm{LP}$, where $P$ stands for the transmitted power. By minimizing the trace of the worst-case CRB matrix, this issue can be illustrated as

## (8)

$$\underset{{\mathbf{R}}_{\mathbf{\Phi}}}{\mathrm{min}}\text{\hspace{0.17em}}\underset{{\{{\mathbf{\delta}}_{k}\}}_{k=1}^{K},{\{{\dot{\mathbf{\delta}}}_{k}\}}_{k=1}^{K}}{\mathrm{max}}\text{\hspace{0.17em}}\mathrm{tr}(\mathbf{C})\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{0.17em}\hspace{0.17em}}{\mathbf{\delta}}_{k}\in {\mathcal{U}}_{1},{\dot{\mathbf{\delta}}}_{k}\in {\mathcal{U}}_{2}\phantom{\rule{0ex}{0ex}}\mathrm{tr}({\mathbf{R}}_{\mathbf{\Phi}})=\mathrm{LP}\phantom{\rule{0ex}{0ex}}{\mathbf{R}}_{\mathbf{\Phi}}\u2ab0\mathbf{0},$$^{29}It is not easy to find a method, e.g., the convex optimization approach,

^{29}to acquire a satisfactory solution with accepted computation cost.

## 3.

## Diagonal Loading-Based Iterative Method

To tackle the NP-hard problem shown in Eq. (8), a new DL-based iterative method is developed in this section. To proceed, we first consider the inner maximization problem. To simplify this issue, we present the following lemma:^{31}

## Lemma 1 (Hadamard’s Inequality).

Assume $\mathbf{M}$ be an $N\times N$ positive semidefinite Hermitian matrix, then the following inequality

holds, where the equality is achieved if and only if $\mathbf{M}$ is diagonal.With Lemma 1, the inner optimization in Eq. (8) can be simplified as

## (10)

$$\underset{{\{{\mathbf{\delta}}_{k}\}}_{k=1}^{K},{\{{\dot{\mathbf{\delta}}}_{k}\}}_{k=1}^{K}}{\mathrm{max}}\text{\hspace{0.17em}}\sum _{k=1}^{K}\frac{1}{{[2\text{\hspace{0.17em}}\mathrm{Re}(\mathbf{F})]}_{kk}}\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{0.17em}\hspace{0.17em}}{\mathbf{\delta}}_{k}\in {\mathcal{U}}_{1},{\dot{\mathbf{\delta}}}_{k}\in {\mathcal{U}}_{2}\mathrm{.}$$## (11)

$$\underset{{\{{\mathbf{\delta}}_{k}\}}_{k=1}^{K},{\{{\dot{\mathbf{\delta}}}_{k}\}}_{k=1}^{K}}{\mathrm{max}}\text{\hspace{0.17em}}\sum _{k=1}^{K}\frac{1}{{\beta}_{k}^{*}{\dot{\tilde{\mathbf{h}}}}_{k}^{H}\{{[\mathbf{I}+({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\dot{\tilde{\mathbf{h}}}}_{k}{\beta}_{k}+{\tilde{\mathbf{h}}}_{k}^{H}\{{[\mathbf{I}+({\mathbf{R}}_{\mathrm{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\tilde{\mathbf{h}}}_{k}}\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{1em}}\text{\hspace{1em}}{\mathbf{\delta}}_{k}\in {\mathcal{U}}_{1},{\dot{\mathbf{\delta}}}_{k}\in {\mathcal{U}}_{2}$$Equation (11) shows that the denominator of the $k$’th term in the summation only relies on ${\mathbf{\delta}}_{k}$ and ${\dot{\mathbf{\delta}}}_{k}$. As a consequence, the problem Eq. (11) can be regarded as maximization of every term in the summation subject to the corresponding constraint, i.e., it can be written as $K$ independent problems shown as

## (12)

$$\underset{\{{\mathbf{\delta}}_{k}{\}}_{k=1}^{K},\{{\dot{\mathbf{\delta}}}_{k}{\}}_{k=1}^{K}}{\mathrm{max}}\frac{1}{{\beta}_{k}^{*}{\dot{\tilde{\mathbf{h}}}}_{k}^{H}\{{[\mathbf{I}+({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\dot{\tilde{\mathbf{h}}}}_{k}{\beta}_{k}+{\tilde{\mathbf{h}}}_{k}^{H}\{{[\mathbf{I}+({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\tilde{\mathbf{h}}}_{k}}\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{1em}}\text{\hspace{1em}}{\Vert {\dot{\mathbf{\delta}}}_{k}\Vert}_{F}\le {\sigma}_{k},{\Vert {\mathbf{\delta}}_{k}\Vert}_{F}\le {\zeta}_{k}\phantom{\rule[-0.0ex]{1em}{0.0ex}}k=\mathrm{1,2},\dots ,K.$$As $\mathbf{P}\succ \mathbf{0},{\mathbf{R}}_{\mathbf{\Phi}}\u2ab0\mathbf{0},{\mathbf{R}}_{{\mathbf{H}}_{c}}\u2ab0\mathbf{0}$, then $({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}$ is indefinite matrix,^{31} thus Eq. (12) is not easy to be solved. To tackle it, the DL approach^{30} is employed to ${\mathbf{R}}_{\mathrm{\Phi}}$ such that

## (13)

$${\tilde{\mathbf{R}}}_{\mathbf{\Phi}}={\mathbf{R}}_{\mathbf{\Phi}}+\u03f5\mathbf{I}\succ \mathbf{0},$$It is evident that ${\beta}_{k}^{*}{\dot{\tilde{\mathbf{h}}}}_{k}^{H}\{{[\mathbf{I}+({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})]{\dot{\tilde{\mathbf{h}}}}_{k}{\beta}_{k}$ is convex w.r.t. ${\dot{\mathbf{\delta}}}_{k}$, and ${\tilde{\mathbf{h}}}_{k}^{H}\{{[\mathbf{I}+({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\tilde{\mathbf{h}}}_{k}$ w.r.t. ${\mathbf{\delta}}_{k}$.^{29} Consequently, Eq. (12) can be regarded as minimization of this term w.r.t. ${\dot{\mathbf{\delta}}}_{k}$ and ${\mathbf{\delta}}_{k}$ can be expressed as

## (14)

$$\underset{{\dot{\mathbf{\delta}}}_{k},{\mathbf{\delta}}_{k}}{\mathrm{min}}\text{\hspace{0.17em}}{\beta}_{k}^{*}{\dot{\tilde{\mathbf{h}}}}_{k}^{H}\{{[\mathbf{I}+({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\dot{\tilde{\mathbf{h}}}}_{k}{\beta}_{k}+{\tilde{\mathbf{h}}}_{k}^{H}\{{[\mathbf{I}+({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\tilde{\mathbf{h}}}_{k}\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{0.17em}\hspace{0.17em}}{\Vert {\dot{\mathbf{\delta}}}_{k}\Vert}_{F}\le {\sigma}_{k},{\Vert {\mathbf{\delta}}_{k}\Vert}_{F}\le {\zeta}_{k}.$$Similar to that of Eq. (11), it can be observed that each term in the summation only depends on ${\mathbf{\delta}}_{k}$ or ${\dot{\mathbf{\delta}}}_{k}$ and the uncertainty in ${\{{\mathbf{\delta}}_{k}\}}_{k=1}^{K}$ and ${\{{\dot{\mathbf{\delta}}}_{k}\}}_{k=1}^{K}$ is independent of each other. Following that, the problem above can be illustrated as the following two separate optimization issues:

## (15)

$$\underset{{\dot{\mathbf{\delta}}}_{k}}{\mathrm{min}}\text{\hspace{0.17em}}{\beta}_{k}^{*}{\dot{\tilde{\mathbf{h}}}}_{k}^{H}\{{[\mathbf{I}+({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\dot{\tilde{\mathbf{h}}}}_{k}{\beta}_{k}\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{0.17em}\hspace{0.17em}}{\Vert {\dot{\mathbf{\delta}}}_{k}\Vert}_{F}\le {\sigma}_{k}$$## (16)

$$\underset{{\mathbf{\delta}}_{k}}{\mathrm{min}}\text{\hspace{0.17em}}{\tilde{\mathbf{h}}}_{k}^{H}\{{[\mathbf{I}+({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\tilde{\mathbf{h}}}_{k}\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{0.17em}\hspace{0.17em}}{\Vert {\mathbf{\delta}}_{k}\Vert}_{F}\le {\zeta}_{k}.$$## (17)

$$\underset{{\dot{\mathbf{\delta}}}_{k},t}{\mathrm{min}}\text{\hspace{0.17em}}t\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{0.17em}\hspace{0.17em}}{\beta}_{k}^{*}{\dot{\tilde{\mathbf{h}}}}_{k}^{H}\{{[\mathbf{I}+({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\dot{\tilde{\mathbf{h}}}}_{k}{\beta}_{k}\le t\phantom{\rule{0ex}{0ex}}{\dot{\mathbf{\delta}}}_{k}^{H}{\dot{\mathbf{\delta}}}_{k}\le {\sigma}_{k}^{2}$$## (18)

$$\underset{{\mathbf{\delta}}_{k},t}{\mathrm{min}}\text{\hspace{0.17em}}t\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{0.17em}\hspace{0.17em}}{\tilde{\mathbf{h}}}_{k}^{H}{\{(\mathbf{I}+({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})\}{\tilde{\mathbf{h}}}_{k}\le t\phantom{\rule{0ex}{0ex}}{\mathbf{\delta}}_{k}^{H}{\mathbf{\delta}}_{k}\le {\zeta}_{k}^{2},$$The problems Eqs. (17) and (18) can be reformulated into semidefinite programming (SDP) issues according to the following lemma:^{31}

## Lemma 2 (Schur’s Complement).

Assume $\mathbf{Z}$ is a Hermitian matrix, which can be partitioned as $\mathbf{Z}=\left[\begin{array}{cc}{\mathbf{Z}}_{11}& {\mathbf{Z}}_{12}^{H}\\ {\mathbf{Z}}_{21}& {\mathbf{Z}}_{22}\end{array}\right]$, then $\mathbf{Z}\u2ab0\mathbf{0}$ if and only if ${\mathbf{Z}}_{22}\succ \mathbf{0}$ and ${\mathbf{Z}}_{11}-{\mathbf{Z}}_{12}^{H}{\mathbf{Z}}_{22}^{-1}{\mathbf{Z}}_{21}\u2ab0\mathbf{0}$.

With Lemma 2, the solutions to Eqs. (15) and (16) can be acquired, respectively, by the following SDPs:

## (19)

$$\underset{t,{\dot{\mathbf{\delta}}}_{k}}{\mathrm{min}}\text{\hspace{0.17em}}t\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{0.17em}\hspace{0.17em}}\left[\begin{array}{cc}{\sigma}_{k}^{2}& {\dot{\mathbf{\delta}}}_{k}^{H}\\ {\dot{\mathbf{\delta}}}_{k}& \mathbf{I}\end{array}\right]\u2ab0\mathbf{0}\phantom{\rule{0ex}{0ex}}\left[\begin{array}{cc}t& {\beta}_{k}^{*}{\dot{\tilde{\mathbf{h}}}}_{k}^{H}\\ {\beta}_{k}{\dot{\tilde{\mathbf{h}}}}_{k}& {\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}\otimes \mathbf{P}+{\mathbf{R}}_{{\mathbf{H}}_{c}}\end{array}\right]\u2ab0\mathbf{0}$$## (20)

$$\underset{t,{\mathbf{\delta}}_{k}}{\mathrm{min}}\text{\hspace{0.17em}}t\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{0.17em}\hspace{0.17em}}\left[\begin{array}{cc}{\zeta}_{k}^{2}& {\mathbf{\delta}}_{k}^{H}\\ {\mathbf{\delta}}_{k}& \mathbf{I}\end{array}\right]\u2ab0\mathbf{0}\phantom{\rule{0ex}{0ex}}\left[\begin{array}{cc}t& {\tilde{\mathbf{h}}}_{k}^{H}\\ {\tilde{\mathbf{h}}}_{k}& {\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}\otimes \mathbf{P}+{\mathbf{R}}_{{\mathbf{H}}_{c}}\end{array}\right]\u2ab0\mathbf{0}.$$Substituting ${\{{\dot{\mathbf{\delta}}}_{k}\}}_{k=1}^{K}$ and ${\{{\mathbf{\delta}}_{k}\}}_{k=1}^{K}$ obtained from Eqs. (19) and (20) into Eq. (8), we now consider the outer optimization problem.

With Eq. (13), the following proposition can recast the nonlinear objective w.r.t. ${\mathbf{R}}_{\mathbf{\Phi}}$ in Eq. (8) as a linear issue and reformulate the constraints into a linear matrix inequality (LMI) form, the proof of which can be found in the Appendix.

## Proposition.

Exploiting some matrix operations, the constraints in Eq. (8) can be reshaped as

where## (22)

$$\mathbf{D}={[\mathbf{I}+({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})$$With Lemma 2 as well as the proposition above, similar to the inner optimization problem, the outer optimization issue in Eq. (8) is equivalent to an SDP

## (23)

$$\underset{\mathbf{X},\mathbf{D}}{\mathrm{min}}\text{\hspace{0.17em}}\mathrm{tr}(\mathbf{X})\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{0.17em}\hspace{0.17em}}\alpha \mathbf{I}\u2aaf\mathbf{D}\u2aaf\beta \mathbf{I}\phantom{\rule{0ex}{0ex}}\left[\begin{array}{cc}\mathbf{X}& \mathbf{I}\\ \mathbf{I}& \mathbf{F}\end{array}\right]\u2ab0\mathbf{0},$$After obtaining the optimum $\mathbf{D}$, ${\mathbf{R}}_{\mathbf{\Phi}}$ can be obtained via the LS fitting of it, which can be expressed as

## (24)

$${\mathbf{R}}_{\mathbf{\Phi}}=\mathrm{arg}\underset{{\mathbf{R}}_{\mathbf{\Phi}}}{\mathrm{min}}{\Vert {({\mathbf{D}}^{-1}-{\mathbf{R}}_{{\mathbf{H}}_{c}})}^{-1}-{\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}\Vert}_{F}\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{tr}({\mathbf{R}}_{\mathbf{\Phi}})=\mathrm{LP}\phantom{\rule{0ex}{0ex}}{\mathbf{R}}_{\mathbf{\Phi}}\u2ab0\mathbf{0},$$## (25)

$$\underset{{\mathbf{R}}_{\mathbf{\Phi}},t}{\mathrm{min}}\text{\hspace{0.17em}}t\phantom{\rule{0ex}{0ex}}\mathrm{s}.\mathrm{t}.\text{\hspace{0.17em}\hspace{0.17em}}\left[\begin{array}{cc}t& {\mathrm{vec}}^{H}[{({\mathbf{D}}^{-1}-{\mathbf{R}}_{{\mathbf{H}}_{c}})}^{-1}-{\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}]\\ \mathrm{vec}[{({\mathbf{D}}^{-1}-{\mathbf{R}}_{{\mathbf{H}}_{c}})}^{-1}-{\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}]& \mathbf{I}\end{array}\right]\u2ab0\mathbf{0}\phantom{\rule{0ex}{0ex}}\mathrm{tr}({\mathbf{R}}_{\mathbf{\Phi}})=\mathrm{LP}\phantom{\rule{0ex}{0ex}}{\mathbf{R}}_{\mathbf{\Phi}}\u2ab0\mathbf{0}.$$By solving the inner and outer optimization issue, according to Algorithm 3 given in Ref. 17, a DL-based iterative approach can be proposed to minimize the worst-case CRB, which is illustrated as follows.

**Algorithm.** Given the uncertainty sets, i.e., ${\mathcal{U}}_{1}$ and ${\mathcal{U}}_{2}$, and the initial value of ${\mathbf{R}}_{\mathbf{\Phi}}$ (uncorrelated waveforms can be considered as this initial value in the following numerical examples), then the WCM and the parameter estimation error can be obtained by the following steps.

1. Compute Eqs. (19) and (20) to obtain ${\dot{\mathbf{\delta}}}_{k}$, ${\mathbf{\delta}}_{k}$.

2. Compute Eq. (23) to obtain $\mathbf{D}$.

3. Repeat steps 1 and 2 until the worst-case CRB improvement becomes insignificant, the criterion $\Vert {\mathbf{C}}^{(i)}-{\mathbf{C}}^{(i-1)}\Vert \le {10}^{-9}$ is exploited in the following numerical examples, where $i$ is the iteration number.

Following that, an optimal WCM, i.e., ${\mathbf{R}}_{\mathbf{\Phi}}$, can be obtained by solving Eq. (25).

The solutions to Eqs. (19), (20), (24), and (25) can be acquired very effectively by the CVX optimization toolbox in Ref. 32. It is noted that only the WCM is obtained in this paper. In order to generate the ultimate waveforms based on the optimized WCM, some practical constraints about the transmitted waveforms, e.g., the peak-to-average power ratio, finite-alphabet sets, and so on, should be taken into account (see Ref. 25 and the cited references for more details).

## 4.

## Numerical Simulations

In the following, several numerical simulations are given to demonstrate the benefits of the proposed method, compared to the nonrobust method proposed in Ref. 22, the robust MMSE-based method developed in Ref. 21, and uncorrelated waveforms, which can be illustrated from the following perspectives: the improvement of the worst-case parameter estimation performance, the robustness of the three methods, and the effect of the bound of uncertainties of ${\tilde{\mathbf{h}}}_{k}$ and ${\dot{\tilde{\mathbf{h}}}}_{k}$ on the worst-case CRB.

A linear MIMO radar antenna array of ${M}_{t}=3$ transmit elements and ${M}_{r}=3$ receive elements is exploited in the following numerical simulations with various configurations: half-wavelength spacing between adjacent elements in both transmitter and receiver, denoted by MIMO radar 1, and one-and-half-wavelength spacing between adjacent elements of transmitter and half-wavelength spacing in receiver, denoted by MIMO radar 2. It is assumed that ${N}_{C}=\mathrm{10,000}$ clutter patches are uniformly distributed and homogeneous within the considered range bin and the clutter-to-noise ratio is 30 dB. Due to the limited number of receive element, we assume only one jammer at $-5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$ with jammer-to-noise ratio being 60 dB, and a unit amplitude for the target of interest at $\theta =20\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$. Besides, we set the number of snapshots $L=256$.

From the discussion in Sec. 2, it can be seen that the calculation of CRB must exploit the initial location parameter estimate. These parameter estimates can be obtained by using various methods proposed in Ref. 26. Moreover, for the convenience of comparison, we need to calculate the CRB of the robust MMSE-based method. Similar to the derivation of that in Ref. 22, based on the model shown in Ref. 21, the elements of Fisher information matrix for the robust MMSE-based method can be represented as follows:

## (26)

$${\mathbf{F}}_{11}=(\{{[(\dot{\mathbf{B}}*\mathbf{A})+(\mathbf{B}*\dot{\mathbf{A}})]}^{H}{[\mathbf{I}+({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})[(\dot{\mathbf{B}}*\mathbf{A})+(\mathbf{B}*\dot{\mathbf{A}})]\})\odot ({\mathbf{\beta}}^{*}{\mathbf{\beta}}^{T})$$## (27)

$${\mathbf{F}}_{12}=\mathrm{diag}({\mathbf{\beta}}^{*})\{{(\dot{\mathbf{B}}*\mathbf{A}+\mathbf{B}*\dot{\mathbf{A}})}^{H}{[\mathbf{I}+({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})(\mathbf{B}*\mathbf{A})\}$$## (28)

$${\mathbf{F}}_{22}={(\mathbf{B}*\mathbf{A})}^{H}{[\mathbf{I}+({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\mathbf{R}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})(\mathbf{B}*\mathbf{A})$$## (29)

$$\mathbf{A}=[\begin{array}{cccc}\mathbf{a}({\theta}_{1}),& \mathbf{a}({\theta}_{2}),& \dots ,& \mathbf{a}({\theta}_{K})\end{array}],\phantom{\rule[-0.0ex]{1em}{0.0ex}}\mathbf{B}=[\begin{array}{cccc}\mathbf{b}({\theta}_{1}),& \mathbf{b}({\theta}_{2}),& \dots ,& \mathbf{b}({\theta}_{K})\end{array}],\mathbf{\beta}={[\begin{array}{cccc}{\beta}_{1},& {\beta}_{2},& \dots ,& {\beta}_{K}\end{array}]}^{T}$$## (30)

$$\dot{\mathbf{A}}=\left[\begin{array}{ccc}\frac{\partial \mathbf{a}({\theta}_{1})}{\partial {\theta}_{1}}& \cdots & \frac{\partial \mathbf{a}({\theta}_{K})}{\partial {\theta}_{K}}\end{array}\right],\phantom{\rule[-0.0ex]{1em}{0.0ex}}\dot{\mathbf{B}}=\left[\begin{array}{ccc}\frac{\partial \mathbf{b}({\theta}_{1})}{\partial {\theta}_{1}}& \cdots & \frac{\partial \mathbf{b}({\theta}_{K})}{\partial {\theta}_{K}}\end{array}\right],$$In the following, the performance of the proposed method is demonstrated in two cases, i.e., one is that only the initial angle estimate error is considered, and the other is that only the calibration error in both the transmitting and receiving arrays is considered.

## 4.1.

### Uncertainty in Initial Angle Estimation

In this scenario, we assume that the initial angle estimate lies in an uncertainty $\mathrm{\Delta}\theta =[-3\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg},3\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}]$, i.e., $\tilde{\theta}=[17\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg},23\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}]$, where $\tilde{\theta}$ is the estimate of $\theta $. After calculating, we can obtain $\zeta =5.4382$, $\sigma =7.6593$ for MIMO radar 1 and $\zeta =27.6329$, $\sigma =29.6754$ for MIMO radar 2.

First, we show the transmitting beampatterns optimized by the proposed approach when the array signal-to-noise ratio (ASNR) is equal to 10 dB, where ASNR is defined as ${\mathrm{PM}}_{t}{M}_{r}/{\sigma}_{\mathbf{W}}^{2}$. Note the proposed approach puts a peak around the target location, which implies the maximum power can be emitted toward the location with the worst-case CRB in the estimation error set and then the worst-case accuracy can be improved. Besides, it can be seen that some grating lobes are shown in Fig. 1(b) for MIMO radar 2, which is due to the sparse transmitting array.^{2}

Next, the worst-case CRBs obtained by the proposed method, the robust MMSE-based algorithm proposed in Ref. 21, and that of uncorrelated waveforms against ASNR are compared in Fig. 2 to verify the improvement of the worst-case parameter estimation performance with the CRB acquired by the nonrobust method with perfect knowledge of $\tilde{\mathbf{h}}$ and $\dot{\tilde{\mathbf{h}}}$ as a benchmark. Obviously, the worst-case CRB of the four waveforms decreases with the increase of ASNR. Also, it can be seen the waveforms generated by the proposed method outperform uncorrelated waveforms and that of the robust MMSE-based algorithm significantly while the waveforms obtained by the robust MMSE-based algorithm surpass uncorrelated waveforms considerably, regardless of ASNR. The reason for these can be illustrated as following: the optimized waveforms generated by the proposed method that aims at lowering the worst-case CRB focus the transmitted energy on the uncertainty set of the initial parameter estimates while uncorrelated waveforms emit omnidirectionally; the robust MMSE-based algorithm is to reduce the worst-case MMSE rather than CRB by designing the transmitted waveforms, yet the variance obtained by the MMSE estimator is lower than CRB considerably.^{37} Moreover, it can be seen that the gap between the worst-case CRB generated by the proposed method and the CRB acquired by the nonrobust method with exactly known $\tilde{\mathbf{h}}$ and $\dot{\tilde{\mathbf{h}}}$ is rather small, which means that the worst-case parameter estimation performance can be improved effectively by the proposed method. Besides, we can obtain a better CRB in Fig. 2(b) as compared with that of Fig. 2(a), because the virtual receive array aperture produced by MIMO radar 2 is larger than MIMO radar 1.^{2}

In Fig. 3, the average worst-case CRBs obtained by the proposed method and the robust MMSE-based algorithm versus ASNR are shown, as compared with those of uncorrelated waveforms and the nonrobust method, which is based on 100 Monte-Carlo simulations. One can observe that the transmitted waveforms obtained by the proposed method have a lower worst-case CRB than the others, and the waveforms generated by the robust MMSE-based algorithm outperform the nonrobust method and uncorrelated waveforms. In other words, the proposed method has a better robustness against $\zeta $ and $\sigma $.

To verify the effect of the bound of the uncertainty of $\delta $ and $\dot{\delta}$, i.e., $\zeta $ or $\sigma $, on the performance of the proposed method, the worst-case CRBs obtained by the proposed method and uncorrelated waveforms are compared in Fig. 4, as a function of various size of the considered uncertainty sets, where $\mathrm{ASNR}=10\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$. One can see that the worst-case CRB of these two methods increases as the increasing of $\zeta $ or $\sigma $, which is similar to that of Ref. 36. Moreover, it is noted that the proposed method always outperforms uncorrelated waveforms.

## 4.2.

### Calibration Error in the Transmitting and Receiving Arrays

In this scenario, both the transmitting and receiving arrays are assumed to have calibration errors (the sensor amplitude and phase error as well as position error). Each element of the transmit and receive array steering vectors is perturbed by a CSCG random variable with zero-mean and variance ${\sigma}_{e}^{2}=0.03$. After calculating, we can obtain $\zeta =13.4764$, $\sigma =14.5712$ for MIMO 1, $\zeta =29.8362$, $\sigma =32.6573$ for MIMO 2.

Figure 5 depicts the beampattern generated by the proposed method when $\mathrm{ASNR}=10\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{dB}$. From Fig. 5, we can draw a conclusion similar to that of Fig. 1.

The worst-case CRB obtained by the proposed method versus ASNR for this case is shown in Fig. 6 and that of the robust MMSE-based algorithm and uncorrelated waveforms, as well as the nonrobust method. It is obvious that the results obtained from Fig. 6 is similar to that of Fig. 2.

In the following, four average worst-case CRBs versus ASNR in this case are shown in Fig. 7 to verify the robustness of the proposed method, where 100 Monte-Carlo simulations are used. One can see that the results drawn from Fig. 7 are similar to that of Fig. 3.

Finally, the effect of the bound of the uncertainty of $\mathbf{\delta}$ and $\dot{\mathbf{\delta}}$ on the performance of the proposed method for this case is shown in Fig. 8. Similar conclusions to that of Fig. 4 can be obtained from Fig. 8.

## 5.

## Conclusions

The robust waveform design issue for collocated MIMO radar has been studied in this work aiming to improve the worst-case estimation accuracy in the presence of clutter for the case that the prior knowledge of the receive and transmit array steering vectors of targets of interest is imprecise and lies in a convex uncertainty set. For the purpose of improving the worst-case accuracy, the CRB matrix is considered as an important optimization metric and the cost function of the robust design is established. To tackle the formulated complicated problem, a new DL-based iterative algorithm has been developed to optimize the WCM to acquire the better worst-case accuracy, each step in which can be relaxed as an SDP issue, and therefore can be solved very effectively. An optimal solution to the initial problem has been obtained via the LS fitting of the solution acquired by the iterative approach. Numerical examples have demonstrated the effectiveness and superiority of the proposed algorithm compared to the nonrobust method, the robust MMSE-based algorithm, and uncorrelated waveforms.

## Appendices

## Appendix:

### Proof of Proposition

For the purpose of recasting the function in Eq. (8) as a linear one, note that

## (31)

$${[\mathbf{I}+({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})={[{({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})}^{-1}+{({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1})}^{-1}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}\otimes {\mathbf{P}}^{-1}){\mathbf{R}}_{{\mathbf{H}}_{c}}]}^{-1}\phantom{\rule{0ex}{0ex}}={({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}\otimes \mathbf{P}+{\mathbf{R}}_{{\mathbf{H}}_{c}})}^{-1}$$## (32)

$$\mathbf{D}={({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}\otimes \mathbf{P}+{\mathbf{R}}_{{\mathbf{H}}_{c}})}^{-1}.$$## (34)

$${\mathbf{D}}^{-1}={\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}\otimes \mathbf{P}+{\mathbf{R}}_{{\mathbf{H}}_{c}}.$$According to Weyl’s theorem,^{31} we can obtain

## (35)

$${\lambda}_{\mathrm{max}}({\mathbf{D}}^{-1})\le {\lambda}_{\mathrm{max}}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}\otimes \mathbf{P})+{\lambda}_{\mathrm{max}}({\mathbf{R}}_{{\mathbf{H}}_{c}}),\phantom{\rule{0ex}{0ex}}{\lambda}_{\mathrm{min}}({\mathbf{D}}^{-1})\ge {\lambda}_{\mathrm{min}}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}\otimes \mathbf{P})+{\lambda}_{\mathrm{min}}({\mathbf{R}}_{{\mathbf{H}}_{c}}),$$## (36)

$${\lambda}_{\mathrm{max}}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}\otimes \mathbf{P})\le {\lambda}_{\mathrm{max}}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}){\lambda}_{\mathrm{max}}(\mathbf{P}),\phantom{\rule{0ex}{0ex}}{\lambda}_{\mathrm{min}}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}\otimes \mathbf{P})\ge {\lambda}_{\mathrm{min}}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}){\lambda}_{\mathrm{min}}(\mathbf{P}).$$## (37)

$${\lambda}_{\mathrm{max}}({\mathbf{D}}^{-1})\le {\lambda}_{\mathrm{max}}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}){\lambda}_{\mathrm{max}}(\mathbf{P})+{\lambda}_{\mathrm{max}}({\mathbf{R}}_{{\mathbf{H}}_{c}}),\phantom{\rule{0ex}{0ex}}{\lambda}_{\mathrm{min}}({\mathbf{D}}^{-1})\ge {\lambda}_{\mathrm{min}}({\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}){\lambda}_{\mathrm{min}}(\mathbf{P})+{\lambda}_{\mathrm{min}}({\mathbf{R}}_{{\mathbf{H}}_{c}}).$$## (38)

$$\frac{1}{\mathrm{LP}+\u03f5}\mathbf{I}\u2aaf{\tilde{\mathbf{R}}}_{\mathbf{\Phi}}^{-1}\u2aaf\frac{1}{\u03f5}\mathbf{I}.$$## (39)

$${\lambda}_{\mathrm{max}}({\mathbf{D}}^{-1})\le \frac{1}{\u03f5}{\lambda}_{\mathrm{max}}(\mathbf{P})+{\lambda}_{\mathrm{max}}({\mathbf{R}}_{{\mathbf{H}}_{c}}),\phantom{\rule{0ex}{0ex}}{\lambda}_{\mathrm{min}}({\mathbf{D}}^{-1})\ge \frac{1}{\mathrm{LP}+\u03f5}{\lambda}_{\mathrm{min}}(\mathbf{P})+{\lambda}_{\mathrm{min}}({\mathbf{R}}_{{\mathbf{H}}_{c}}),$$## Acknowledgments

The authors would like to thank the anonymous reviewers for their thoughtful and to-the-point comments and suggestions, which greatly improved the manuscript. This work was supported in part by the National Natural Science Foundation of China under Grant No. 61301258, in part by China Postdoctoral Science Foundation Funded Project under Grant No. 2016M590218, in part by the Education Department of Liaoning Province Natural Science Research Program under Grant No. L2015026, and in part by the Key Scientific and Technological Project of Henan Province under Grant No. 142102210599.

## References

## Biography

**Rongyan Guo** received her BS degree in electronic engineering from Zhengzhou University in 2000 and her MS degree in communication engineering from Huazhong University of Science and Technology (HUST) in 2008. She joined the School of Physics and Electrical Engineering, Zhoukou Normal University in 2000, where she is currently an associate professor. Her research interests include space–time adaptive processing and signal processing for radar.

**Hongyan Wang** received his PhD from Xidian University in 2012. He joined the College of Information Engineering, Dalian University, in 2013, where he is currently a lecturer. His research interests include MIMO radar waveform optimization, MIMO communication, space–time adaptive processing, and array signal processing.