Adaptive Target Detection and DOA Estimation With Uniform Rectangular Arrays in the Presence of Unknown Mutual Coupling

This paper investigates joint adaptive target detection and direction of arrival (DOA) estimation via a uniform rectangular array (URA) affected by mutual coupling. Capitalizing on a bespoke linearization of the array manifold and leveraging the banded symmetric Toeplitz block Toeplitz structure for the coupling matrix description, a vectorial model of the useful target echo return is proposed and used to formulate the detection problem. Two decision rules are designed, i.e., the generalized likelihood ratio (GLR) and multifamily likelihood ratio test (MFLRT), with the latter aimed at handling an unknown number of active mutual coupling coefficients. Both demand the joint maximum likelihood (ML) estimates of the coupling coefficients and the target angular displacement parameters which can be obtained solving a non-convex optimization problem. Toward this goal, an iterative procedure based on the minorization-maximization (MM) algorithm is developed. At the analysis stage, the performance of the proposed methods is assessed in terms of detection probability ( $P_{d}$ ) and DOA root mean square error (RMSE) in comparison with benchmarks and standard strategies that do not account for the mutual coupling phenomenon. The results demonstrate the effectiveness of the proposed approaches to overcome signal mismatches induced by both the DOA uncertainty and mutual coupling.

Abstract-This paper investigates joint adaptive target detection and direction of arrival (DOA) estimation via a uniform rectangular array (URA) affected by mutual coupling. Capitalizing on a bespoke linearization of the array manifold and leveraging the banded symmetric Toeplitz block Toeplitz structure for the coupling matrix description, a vectorial model of the useful target echo return is proposed and used to formulate the detection problem. Two decision rules are designed, i.e., the generalized likelihood ratio (GLR) and multifamily likelihood ratio test (MFLRT), with the latter aimed at handling an unknown number of active mutual coupling coefficients. Both demand the joint maximum likelihood (ML) estimates of the coupling coefficients and the target angular displacement parameters which can be obtained solving a non-convex optimization problem. Toward this goal, an iterative procedure based on the minorization-maximization (MM) algorithm is developed. At the analysis stage, the performance of the proposed methods is assessed in terms of detection probability ( P d ) and DOA root mean square error (RMSE) in comparison with benchmarks and standard strategies that do not account for the mutual coupling phenomenon. The results demonstrate the effectiveness of the proposed approaches to overcome signal mismatches induced by both the DOA uncertainty and mutual coupling.

I. INTRODUCTION
A MONG the several challenges that modern antenna array signal processing algorithms [1] are faced with, the presence of mutual coupling (MC) is definitely a very crucial issue [2], [3], [4], [5], [6], [7]. In fact, elements of an antenna array may electromagnetically interact with their neighbors, altering their electromagnetic characteristics. This unwanted effect depends on a multitude of factors, including the number, type, and relative orientation of each antenna element, as well as their distance and position [5], [8], [9]. As a result, the actual received steering vector could be possibly mismatched with respect to (w.r.t.) the ideal manifold model assumed at the design stage by standard algorithms, eventually causing a significant performance degradation [2], [10], [11], [12], [13]. Despite of the potential measurement campaigns to quantify and compensate accurately the MC, some residual effects due to an unavoidable imperfect measurements are often present. Moreover, the MC effect is not necessarily stationary [2], making it challenging to keep the array properly calibrated over time. In this regard, a frequent calibration process to ensure optimal performance is time-consuming and complex, effectively wasting valuable radar time. Therefore, it is crucial to investigate adaptive signal processing methods capable of realizing an on-the-fly compensation of the MC (or its residual, if a initial calibration stage has been performed).
An interesting aspect of MC is that the magnitude of the coupling coefficient describing the interaction between two radiating elements decreases rapidly with the distance between them, and it is substantially the same between any two elements that are at the same distance [8]. This property, together with special symmetries induced by the array geometrical configuration, yields a special structure in the MC matrix (MCM) which can be capitalized at the design stage of bespoke signal processing schemes [2], [14], [15], [16], [17], [18].
To proceed further, let us now frame the aforementioned issue within a typical radar context, with emphasis on adaptive target detection [19], [20], [21], [22], [23], [24] and direction of arrival (DOA) estimation [25], [26], [27], [28]. In [13], assuming a uniform linear array (ULA) affected by MC, the problem of jointly detecting the target and estimating its bearing is addressed while accounting for MC and the DOA uncertainty along the design phase. This is achieved by modeling the actual steering vector as the product of a MCM and an approximated steering vector depending affinely on the unknown DOA displacement (w.r.t. the looking direction) [29]. The target detection problem, formulated assuming a homogeneous radar interference environment, is addressed by resorting to the generalized likelihood ratio test (GLRT) [30], [31] and the multifamily likelihood ratio test (MFLRT) [32] strategies, with the latter employed when an apriori knowledge on the number of active mutual coupling coefficients is unavailable (also referred to as the unknown model order case).
Following the guidelines of [13], this work aims at simultaneously detecting and estimating the two-dimensional (2D) DOA with a uniform rectangular array (URA) in the presence of MC. Specifically, the main innovative technical achievements as compared to [13] can be summarized as follows: • For URA, the MCM demands a different structural model w.r.t. the ULA case. Specifically, a banded symmetric Toeplitz block Toeplitz MCM [18] is necessary instead of the banded symmetric Toeplitz matrix used for the ULA. From the analytic viewpoint this involves different numbers of parameters to estimate and tailored matrix manipulations to capitalize on the developed matrix model.
• After an appropriate concentration of likelihood functions exploiting the estimates of the interference-plus-noise covariance matrix and the MC coefficients, to handle the joint detection and estimation problem for the URA case it is required the solution of a two-dimensional optimization problem, unlike the ULA sensing scenario [13] involving a one-dimensional problem. The subsequent application of the minorization-maximization (MM) optimization framework for the URA entails solving, at each iteration, a new two-variable box-constrained quadratic optimization problem, wherein its optimal solution is provided in closed form.
• For the URA configuration, the performance of the devised methods is numerically assessed in terms of probability of detection (P d ) and root mean square error (RMSE) of the 2D DOA estimates, via Monte Carlo simulations. The results are also compared with benchmarks and counterparts available in the open literature. In addition, a robustness analysis is conducted to evaluate the effectiveness of the proposed strategies in case of mismatches in the MC model, namely when the actual MCM does not fully exhibit the banded symmetric Toeplitz block Toeplitz structure assumed at the design stage. The paper is organized as follows. The signal model for URA with MC effect is given in Section II. In Section III, the design of decision statistics and parameters estimation processes for known and unknown model order cases are addressed. Numerical simulations are provided in Section IV. Finally, conclusions as well as future research activities are discussed in Section V.

A. Notations
Boldface is used for vectors a (lower case), and matrices A (upper case). The (k, l)-entry (or l-entry) of a generic matrix A (or vector a) is indicated as [ A] k,l (or [a] l ). I and 0 denote respectively the identity matrix and the matrix with zero entries (their size is determined from the context). Moreover, diag(x) indicates the diagonal matrix whose i-th diagonal element is [x] i . The transpose and the conjugate transpose operators are denoted by the symbols (·) T and (·) † , respectively. The determinant and the trace of the matrix A ∈ C N ×N are indicated with | A|, tr{ A}, respectively. ⊗ denotes the Kronecker product. R N and C N are the sets of N -dimensional column vectors of real and complex numbers, respectively. H N and H N ++ represent the set of N × N Hermitian matrices and Hermitian positive definite matrices, respectively. The letter j represents the imaginary unit (i.e., j = √ −1). For any complex number x, |x| indicates the modulus of x and Re{x} denotes its real part. Moreover, for any x ∈ C N , ∥x∥ denotes the Euclidean norm, whereas the Frobenius norm of a matrix A is indicated as ∥ A∥ F . Let f (x) be a real-valued function, ∇ x f (x) denotes the gradient of f (·) w.r.t. x, with the partial derivatives arranged in a column vector. Furthermore, for any x, y ∈ R, max(x, y) returns the maximum between the two argument values. Finally, the symbol O(·) denotes the computational complexity in terms of basic operations.

A. Signal Model With Unknown MC
Let us consider a radar system equipped with a URA consisting of M × N antenna elements on the x y-plane with inter-element spacing equal to d (see Fig. 1).
Assume that a point-like target is located in the far-field with the range R 0 , azimuth φ 0 , and elevation θ 0 . The radar received echo after down-conversion, pulse compression, fasttime sampling, and measurement gathering at the instant of interest, is given by where a ∈ C denotes the unknown complex coefficient accounting for target backscattering as well as the other terms involved in the two-way radar equation, Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. p(u, v) evaluated at (u 0 , v 0 ), with u 0 = cos(φ 0 ) sin(θ 0 ) and v 0 = sin(φ 0 ) sin(θ 0 ) the target location parameters in the space of directional cosines [22], and the spatial steering vectors pertaining to the u and v domains, respectively. Up to now, an ideal array manifold has been considered. In practice, the MC phenomenon within the array elements is almost unavoidable and demands a bespoke description at the modeling stage so as to be appropriately handled at the information extraction level. Therefore, the coupling effects can be accounted for by representing the actual steering vector as where C∈ C M N ×M N is the MCM.
To effectively model C, it is assumed that the MC between any pair of elements at the same distance is essentially the same, 1 whereas it decreases rapidly as their distance increases. Moreover, it can be practically ignored for elements fairly separated in distance from each other [8]. Recently, this standard MCM model, originally developed for linear and circular arrays [2], [15], [16], [17], [33], [34], [35], has been extended to the planar array case [18]. Thereby, assuming that for each radiating element the MC effect is induced only by elements within a specific grid/neighborhood around it [18], P − 1 outgrowing squares could be conceived (see 1 Notice that the assumption is realistic for modern phased array radar systems with plenty of elements. Specifically, it reasonably holds true for internal elements but could not be ensured for those close or exactly on to the edges [5]. Fig. 2 where different colors represent distinct geometric symmetries in the electromagnetic field leakage), resulting in a total of 4P(P − 1) nearest elements mutually coupled with the (m, n)-th (m = 1, . . . , M, n = 1, . . . , N ) element. Accordingly, the MCM is represented accurately as a banded symmetric Toeplitz block Toeplitz matrix [36], whereby each block, denoted by C i ∈ C N ×N , describes the intercolumn and intracolumn MC for i = 1 and i ≥ 2, respectively, experienced by the elements of one column of the URA induced by those lying in the (i − 1)-th adjacent column, if present. Therefore, each matrix block is a banded symmetric Toeplitz matrix. Formally, where is in one-to-one mapping with whereas J m is defined as an N × N matrix having 1s on its m-th (m = 1, 2, . . . , P − 1) upper and lower diagonals, and zeros otherwise, i.e., It is interesting to observe that when the Toeplitz-block-Toeplitz matric C (5) can be factored as C u ⊗C v (with C u and C v banded Toeplitz matrices), multiplying C by the steering vector (4), yields C p = C u p u ⊗ C v p v , which is analytically equivalent to the factored model used in [37] and [38] for Multiple Input Multiple Output (MIMO) systems.
Exploiting the structural features of the overall MCM, i.e., C, it can be further decomposed as with D 0 = I M and D i ∈ C M×M , i = 1, . . . , P − 1, a matrix having 1s on its i-th upper and lower diagonals, and zeros otherwise. To further assess the influence of the MC effect, in the following the angle cosine (cosine similarity) [39] between the ideal and the actual steering vectors is analyzed, by means of As a study case, let us consider a URA with M = 5 and N = 6 experiencing the MC phenomenon, which is modeled assuming P = 3, with the matrices C 1 , C 2 , and C 3 defined as in (6)  18 j] T , respectively. In Fig. 3, (9) is plotted versus θ and φ. Its inspection highlights that a symmetric response is obtained w.r.t. φ = 180 • . Moreover, in some regions the mismatch induced by the MC is quite severe. This is for example the case of ]. In such instances, there is a strong mismatch between the actual and the ideal model as testified by the considerable low values of the cosine similarity.

B. Signal Model via Linearized Array Manifold
Following the guidelines of [29], the ideal received steering vector p(u, v) can be approximated via a first order Taylor expansion around the pointing direction (ū,v) with a resulting functional dependency of the linearized array manifold on the offsets u = u 0 −ū and v = v 0 −v, namely ⊗ṗ v withṗ u andṗ v the vectors whose entries contain the derivative of p u and p v w.r.t. u and v evaluated atū andv, respectively (See Appendix A). Now, letting p u = p u (ū), p v = p v (ū), and including the coupling effects, the echo signal can be expressed as with p am ( u, v) = C p a ( u, v) ∈ C M N the approximated actual steering vector encompassing the MC effect, Before concluding this section, it is worth mentioning that following the same line of reasoning of [13], it can be shown that a sufficient condition for the identifiability of the unknowns in (11) is given by M N ≥ 3P 2 , namely

III. TARGET DETECTION PROBLEM
The target detection problem is formulated as a binary hypothesis test aimed at ascertaining the target presence/absence within the cell under test. Assuming the availability of K secondary data free of any useful target signal (with K > M N ), it can be cast as where • r ∈ C M N and r k ∈ C M N , k = 1, . . . , K , are the vectors of primary and secondary data, respectively; • H( u, v) is functionally dependent on the unknown target DOA displacements w.r.t. the array pointing direction; • b represents the unknown vector accounting for both the complex received target echo return a and the P 2 complex MC coefficients; • n ∈ C M N and n k ∈ C M N , k = 1, . . . , K , denote the interference plus noise components of the received snapshots, modeled as statistically independent, complex, zero-mean, circularly symmetric Gaussian random vectors with unknown positive definite covariance matrix, . . , K . Accordingly, the joint probability density functions (PDFs) of the observations under H 0 and H 1 , can be written respectively as with and The optimal Neyman-Pearson detector for the hypothesis testing problem (12), i.e., the LRT, cannot be implemented without the knowledge of u, v, b, and M. Moreover, another complication connected with the design of a practical detector relies on the fact that the number of unknowns in the array coupling coefficients, namely P 2 − 1, can be either assumed known or unknown. In this respect it is possible to explore GLRT-based detectors when the number of coupling coefficients is known at the design stage, while the case of unknown model order can be studied resorting to the MFLRT criterion.

A. Decision Statistic for Known P
In the following, the target detection problem (12) is addressed under the assumption of P known at the design stage. In this context, resorting to the GLRT criterion, the following decision rule can be defined (16) where γ is the detection threshold set to ensure a desired P f a , A denotes the uncertainty set associated with u, i.e., [−α, α], and B denotes the uncertainty set associated with v, i.e., [−β, β]. In this context, α and β are used to control the uncertainty region related to the unknown target DOA. It is important to emphasize that their values must be carefully chosen to guarantee a high level of accuracy of the array manifold approximation (10). To this end, a viable choice entails setting them to the 3dB single-side beamwidth (u 3d B , v 3d B ) ≜ (0.891/N , 0.891/M) of a planar array when pointing in the boresight direction. 2 Maximizing both the numerator and the denominator of (16) over M and taking their logarithm yields the decision statistic (after standard manipulation), where Hence, letting θ = [ u, v] T ∈ R 2 , the optimization problem in (17) is equivalent to where S = [−α, α]×[−β, β] denotes the nonempty and compact feasible set for θ . Let us now consider the optimization w.r.t. b in (18), where 3 is the Moore-Penrose inverse of H w ( θ ). Thus, by substituting the optimizer of b into (17) leads to where It remains to optimize (21) w.r.t. θ . To this end, let us recast the objective function in (21) as where y ∈ C P 2 , A ∈ H P 2 ++ , andḡ( y, A) is jointly convex w.r.t. A and y.
In order to tackle the challenging optimization problem at hand, an MM-based method [40], [41], [42], [43], [44], [45] is developed, which is an iterative procedure consisting of two distinct steps. The former involves the computation of an appropriate tight minorant (surrogate function) based on the current tentative solution. The latter entails its optimization, generating an updated estimate of the unknowns.
Following the same line of reasoning as in [13], the tangent planeḡ a ( y, A| y 0 , A 0 ) toḡ( y, A) in ( y 0 , A 0 ) is an appropriate choice as surrogate function for the problem at hand. Precisely, it is defined as [13] g a ( y, A| y 0 , where and are the gradients ofḡ( y, A) w.r.t. A and y, respectively. As a result, denoting by θ ⋆(h−1) ] T the output of the MM algorithm at the (h − 1)−th iteration, it yields, with y (h−1) 1) ) while the equality in (26) holds when θ = θ ⋆(h−1) . Then, exploiting the current estimated point θ ⋆(h−1) , the optimization w.r.t. θ at the h-th iteration can be obtained according to the maximization of the right-hand side (RHS) of (26), namely (after some algebra) with θ ⋆(h) being the global (unique) optimum of the objective function at hand (the proof is reported in Appendix C). Moreover, the optimal solution to (27) is obtained according to the following proposition: Proposition 1: The optimal solution θ ⋆(h) to (27) is given by the global optimal point for the unconstrained version of (27), i.e., if this solution is feasible, i.e., θ 1 ∈ S with u and v given in Appendix D. Otherwise, the optimal solution is obtained by restricting the objective function to the boundaries of S and selecting the maximum among the finite set {θ 1 , . . . , θ 5 } of optimal candidate solutions, i.e., where with a ′ , b ′ ±β , a ′′ , and b ′′ ±α are defined in Appendix D. Proof: See Appendix D. A summary of the MM is reported in Algorithm 1, whereby the exit condition is assumed where is the objective function at the h-th iteration and ε 1 > 0 is a user-defined exit threshold.
Obtain θ (h) via Proposition 1; 8. Evaluate Therefore, exploiting the result of Algorithm 1, the GLRT for Linearized Array Manifold (GLRT-LAM) statistic is defined by Furthermore, as long as the number of iterations involved in the MM procedures is limited, the computational complexity is dominated by the evaluation of the sample covariance matrix S, that is O(M N K ). As a final remark, following the same approach as in [13] for the case of ULA, it is straightforward to prove that the derived GLRT-LAM detection architecture verifies the bounded-constant false alarm rate (CFAR) property.

B. Decision Statistic for Unknown P
In this subsection, the target detection problem is addressed for the case where the actual number of MC coefficients is unknown. This entails an additional complexity in the development of a tailored processing scheme for target information extraction. To deal with this problem, a multiple composite alternative hypothesis testing problem is formulated, with each alternative hypothesis pertaining to a given number of unknown signal parameters. Formally, whereN ≤ M N /2 is the maximum 4 allowed model order, and with Remarkably, since the alternative hypotheses are nested in (36), i.e., H i ⊂ H j , i < j, a viable solution to handle (36) is to resort to the MFLRT framework [32]. Thus, the target presence can be established according to the decision rule where • 2(i 2 + 1) is the number of unknown parameters under the H i hypothesis; 5 G is the GLRT-like statistic (17) derived assuming P = i; •γ is the threshold guaranteeing the demanded P f a ; • u(t) is the unit step function, i.e., u(t) = 1 as long as t ≥ 0 and zero otherwise. Specifically, denoting by θ (i) MM the estimate of the directional cosines displacements provided by Algorithm 1 assuming P = i, Finally, by following the same methodology used in [13] for ULA, it can be easily proven that the derived MLRT-LAM detector (39) also satisfies the bounded-CFAR property.

IV. SIMULATION RESULTS
In this section, numerical examples are provided to evaluate both the detection and estimation capabilities of the proposed processing schemes for a URA with M = 5 and N = 6 in the presence of MC. The inter-element spacing among the antennas is set to half wavelength and the number of secondary data is K = 4M N = 120. Moreover, the uncertainty set where the unknown target DOA is assumed lying is defined by α = u 3d B ≜ 0.891/M = 0.1782 and β = v 3d B ≜ 0.891/N = 0.1485, where u 3d B and v 3d B are the 3 dB singleside beamwidth of a planar array pointing at the boresight direction. It is also assumed that the pointing direction is (θ, φ) = (33 • , 73 • ). As to the target location parameters, two cases are analyzed. The former considers cosine mismatches Moreover, in order to assess the detection performance, the P d is used as evaluation criterion, computed resorting to N MC = 5000 Monte Carlo runs and assuming a Probability of False Alarm (P f a ) equals to 10 −4 . In this respect, 100 P f a Monte Carlo trails are evaluated to set the detection thresholds of the considered detectors. As to the estimation performance, needless to say evaluated under H 1 , the RMSE is considered as figure of merit. It is computed as where u r and v r are the estimates provided at the r -th trial by a given technique. Furthermore, the SINR is defined as For the evaluation of the MFLRT-based methods,N ∈ {2, 3, 4} is analyzed, where the value ofN is specified as subscript. Furthermore, the two-stage (referred to as "2S") architectures of both the derived GLRT-based and MFLRTbased methods is also implemented, namely, after a first computation of the unknown angular displacement θ ] T (as described in Section III), the ideal steering vector (8) is re-linearized around (ū + u ⋆ ,v + v ⋆ ) and the estimation procedure is executed once again. For comparison, the following detectors are also included: • the GLRT using the actual array manifold with known target DOA and known coupling coefficients • the GLRT using the ideal array manifold (not encompassing the MC phenomenon) with known target DOA • the GLRT using the actual array manifold with known target DOA but with the useful target steering vector affected by the MC phenomenon where • the standard GLRT using the ideal array manifold with the nominal pointing direction (ū,v) (which refers to a fully mismatched case) [20] • the Subspace Detector (SD) [46], namely a mismatched GLRT detector using as useful signal directions those given by the columns of H S D , i.e., with H S D = [ p,p u ,p v ]. In the reported simulations, three different scenarios are examined. First, the useful signal is assumed buried in white Gaussian noise; then, a jamming interference situation is assessed. Finally, within this last context, random perturbations in the elements of the MCM are considered.

A. Gaussian White Noise Case
In this case study, the disturbance covariance matrix is modeled as M = σ 2 n I M N with σ 2 n the internal noise power level, which is assumed, without loss of generality, equal to 0 dB. Not surprisingly, the P d results pertaining to the first analyzed DOAs scenario, depicted in Fig. 4 (a), reveal that the GLRT-ben and the GLRT-ben-DOA yield upper bounds to the performance of all the considered detectors. As to the devised strategies, both the GLRT-LAM and the MFLRT (but for the caseN = 2, where an undersized model order is considered), are able to provide satisfactory detection results, yielding performance close to that of the GLRT-ben-DOA. Specifically, the gap between the corresponding curves is less than 1 dB at P d = 0.9. As already slightly pinpointed, inspection of the figure highlights that usingN = 2 for the computation of the MFLRT entails an inadequate signal model that leads to a SINR loss, at P d = 0.9, in the order of 3 dB w.r.t. the GLRT-LAM. Furthermore, the second stage of processing is unable to provide any improvement in the detection performance, with a negligible loss w.r.t. the corresponding one stage curves. This is not surprising, as the best expansion point has already been used in the first linearization. Finally, the performance of the detectors that do not account for the coupling phenomenon, i.e., ben-GLRT-NC and GLRT-mis, is unsatisfactory, with the only exception of the SD, which is instead capable to deliver acceptable performance, although inferior compared to the proposed strategies. This behavior reinforces the need for the development of tailored signal processing schemes to cope with the challenges posed by the MC. The capabilities of the devised strategies are also corroborated by the analysis of the estimation performance reported in Fig. 4 (b), within the u and v domains, respectively. In particular, looking over the plots highlights that both the GLRT and MFLRT methods (with the exception of the MFLRT 2 ) yield reliable estimates of the DOA displacements. Moreover, the 2S variants of the GLRT and MFLRT architectures deliver almost the same estimation performance, with slight improvement only in the high SINR regime.
Regarding the other DOAs scenario, i.e., [ u, v] T = [0.03, 0.05] T , illustrated in Figs. 4 (c) and (d), the detection performance of the devised methods remains relatively unchanged as compared to the previously analyzed case, indicating that they are able to severely handle the mismatch caused by DOA uncertainty and to endow robustness to the joint detection and estimation process in the presence of MC effects. Furthermore, both the MFLRT 2 (and its 2S variant) and the SD scheme experience a performance loss, compared to the case presented in Fig. 4 (a), of approximately 1 dB at P d = 0.9, reasonably due to the stronger mismatch on the signal model. The highlighted behavior of the devised methods is also pinpointed by the inspection of the RMSE curves, depicted in Fig. 4 (d). At low and medium SINR regimes, the achieved estimation performance is quite similar to that analyzed in the previous DOAs scenario, with some significant differences especially at high SINR. In particular, the estimates of u obtained with the 1S version of the proposed methods saturate, while those obtained with the 2S variants perform considerably better. In fact, the GLRT-LAM-2S provides increasingly accurate estimation performance as the SINR increases. This improvement is difficult to observe in the P d curves since it is attained for SINR values where the detection performance of the aforementioned methods has substantially already achieved saturation.

B. Interference Case
In this section, the radar system is assumed affected by some interference from two narrowband jammers, angularly located at (u 1 , v 1 ) = (0.4698, 0.1710) and (u 2 , v 2 ) = (0.5567, 0.6634), respectively. As a consequence, the interference-plus-noise covariance matrix can be modeled as [47] where p m (u i , v i ) denotes the actual steering vector of the i-th (i = 1, 2) interference source with σ 2 1 and σ 2 2 the corresponding powers. Moreover, it is assumed that the Jammer to Noise Ratio (JNR) of the two emitters is J N R 1 ≜ σ 2 1 /σ 2 n = 30 dB and J N R 2 ≜ σ 2 2 /σ 2 n = 40 dB, respectively. The effectiveness of the methods under consideration are analyzed in Fig. 5, using the same simulation setup considered in Fig. 4 but for the interference scenario which is now modeled according to (47). Specifically, the P d curves evaluated assuming [ u, v] T = [0, 0] T are displayed in Fig. 5 (a), which shows no significant differences as compared to the corresponding white noise case. As a result, all detectors maintain their rankings. Still, a clear performance advantage over the mismatched and SD detectors is present, corroborating the effectiveness of the proposed strategies. Furthermore, the estimation performance, depicted in Fig. 5 (b), is consistent with the white noise disturbance situation as well, with only a slight loss (experienced by all the methods) on the inference of u, especially at the high SINR regime.
On the other hand, regarding the P d behavior for [ u, v] T = [0.03, 0.05] T DOAs case, displayed in Fig. 5 (c), it is possible to observe detection performance similar to the case related to Fig. 4 (c), although with a slight improvement obtained by the 2S variants over the corresponding 1S counterparts. As a matter of fact, the RMSE curves in Fig. 5 (d) indicate a corresponding interesting performance improvement at high values of SINR. Specifically, starting from 22 dB in the u domain and 15 dB in the v domain, some advantages of our refined 2S processing scheme can be observed. Notably, being at this SINR regime the P d not yet achieved saturation, detection gains are experienced too. Remarkably, the 2S strategy exploits the DOA estimation obtained in the first stage and provides, with the second stage, a more accurate approximation of the actual array manifold, resulting in improved performance and yielding accurate DOA estimates without saturation.

C. Perturbations in the Actual MCM
In practice, some perturbations on the nominal banded symmetric Toeplitz block Toeplitz structure of the MCM could be experienced, as for instance induced by radome reflection effects, sensor position misplacements, as well as changes in the environment (e.g., due to the movement of metal objects near the antenna array) [2] that might induce deviations in the nominal conditions of the radar system. To this end, possible random perturbation can be accounted for via the MCM model 6Ĉ where N per t ∈ C N ×N ∼ CN (0, ξ I) is a random perturbation matrix whose complex entries are modeled as independent and identically distributed zero-mean circularly-symmetric complex Gaussian random variables, with variance ξ ≥ 0.
Within the same setup of Figs. 5 (c) and (d) correspondingly, in particular, to [ u, v] T = [0.03, 0.05] T , the aim of the subsequent analysis is a robustness analysis of the devised GLRT-and MFLRT-based signal processing strategies considering, at each Monte Carlo trial, a different realization of the random perturbation term N per t , assuming ξ = 0.1. The results, presented in terms of P d and RMSE versus SINR, are shown in Fig. 6. A comparison of the detection and estimation results with those obtained in the absence of perturbations, analyzed in Figs. 5 (c) and (d), shows that the perturbations induce a detection performance degradation of approximately 1 dB at P d = 0.9 for all considered methods. Moreover, unlike the corresponding ideal case, the MFLRT 4 and its 2S variant exhibit detection capabilities quite close to the ben-GLRT-DOA, with a slight advantage over the GLRT-LAM and GLRT-LAM-2S counterparts. This is attributed to the MFLRT capability to leverage higher model orders, which induces more degrees of freedom for the subsumed array manifold inference, thereby counteracting the effects of perturbations in the structure. This behavior is also reflected in the RMSE curves depicted in Fig. 6 (b). Specifically, at a low SINR regime, the MFLRT 4 yields the best performance in both u and v domains while at high SINR values, it achieves the best performance only for the v estimation, showing a performance saturation in the estimation of u. Still, in the high SINR regime, only the GLRT-LAM and its 2S variant do not experience a performance saturation in both the u and v domains, likely due to their prior knowledge of the number of active MC coefficients. Summarizing, the proposed techniques exhibit a good level of robustness to perturbations in the banded symmetric Toeplitz block Toeplitz structure of the MCM, providing performance levels that are comparable to the nominal case, assumed at the design stage, where the MCM is modeled according to (5).

V. CONCLUSION
Simultaneous target detection and 2D DOA estimation for a URA affected by MC has been addressed. At the modeling stage, a suitable description of the actual echo return has been developed accounting for the unknown target 2D angular displacements w.r.t. the pointing direction and an unknown structured matrix that embeds the effects of MC and target backscattering. Leveraging this bespoke model, the target detection problem in Gaussian interference has been formulated and tackled by resorting to two sub-optimal criteria, i.e., the GLRT and the MFLRT. The former is adopted for the known model order case whereas the latter considers the situation where the number of actual coupling coefficients is unknown. The practical implementation of the detectors has demanded the joint ML estimates of the unknown angular displacements and coupling coefficients which stem as the solution to a non-convex optimization problem. To get a good quality solution, an ad-hoc iterative procedure based on the MM algorithm has been designed. It involves at each iteration the solution of a box-constrained quadratic optimization problem whose optimal point has been obtained in a closed form. Remarkably, both the GLRT and the MFLRT ensure the bounded-CFAR property. At the analysis stage, the performance of the proposed signal processing strategies has been assessed in terms of P d and DOA RMSE in different scenarios, including the presence of jammers as well as random perturbations in the coupling matrix structure assumed at the design stage. Comparisons with benchmark detectors and some counterparts that do not encompass the MC phenomenon have also been conducted. The results have highlighted the effectiveness of the devised architectures showing performance levels comparable to the benchmark that assume known target DOA and model order.
Possible future research avenues might include the analysis on real (or on high-fidelity electromagnetically simulated) phased array data to assess possible discrepancies with the theoretical results, the derivation of additional decision statistics based on different criteria, such as Adaptive Matched Filter (AMF) [48], Rao, and Wald tests [31], as well as the application of the devised framework in a pulse-Doppler scenario accounting, at the design stage, for both the Doppler uncertainty and the MC effect. Finally, the extension to the multi-polarization case [49] could be definitely an interesting complement to this work.
C. Proof of Strict Concavity of the Objective Function in (27) The existence of the global maximizer of (27) is guaranteed by the Weierstrass theorem [50], being the objective function in (27) continuous and the feasible set S non-empty and compact. It is also worth noting that the Hessian matrix of the objective function in (27)   where the first inequality holds for the Cauchy-Schwarz inequality 7 applied to the first term. Hence, H e ≺ 0, indicating that g a ( θ| θ ⋆(h−1) ) is strictly concave in θ .

D. Proof of Proposition 1
Proof: Let us define for simplicity: • g a = g a ( y, A| y (h−1) 0 , A (h−1) 0 ); then according to the first order optimality conditions, the optimal points can be obtained by forcing the gradient of g a w.r.t. u and v to be zero, i.e., ∂g a ∂ u = 0 and ∂g a ∂ v = 0, which yields = y ( θ ⋆(h−1) ) , whose optimal solution is given by v * ± = max(min(−b ′′ ±α /(2a ′′ ), β), −β) with As a consequence, the candidate optimal solutions associated with the left and right edges are and