Parametric Electromagnetic Analysis of Radar-Based Advanced Driver Assistant Systems

Efficient and optimal design of radar-based Advanced Driver Assistant Systems (ADAS) needs the evaluation of many different electromagnetic solutions for evaluating the impact of the radome on the electromagnetic wave propagation. Because of the very high frequency at which these devices operate, with the associated extremely small wavelength, very fine meshes are needed to accurately discretize the electromagnetic equations. Thus, the computational cost of each numerical solution for a given choice of the design or operation parameters, is high (CPU time consuming and needing significant computational resources) compromising the efficiency of standard optimization algorithms. In order to alleviate the just referred difficulties the present paper proposes an approach based on the use of reduced order modeling, in particular the construction of a parametric solution by employing a non-intrusive formulation of the Proper Generalized Decomposition, combined with a powerful phase-angle unwrapping strategy for accurately addressing the electric and magnetic fields interpolation, contributing to improve the design, the calibration and the operational use of those systems.


Introduction
Radar is a widely employed technology which relies on wave propagation to detect surrounding objects. It consists in emitting a radio wave from a transmitter and measuring the reflected wave with a receiving antenna. The data collected from this measurement can provide information about a detected object such as its location or its nature. An in-depth analysis of the electromagnetic field may be necessary to ensure the information inferred from the data is valid.
Radio waves are oscillations of the electromagnetic field and can therefore be computed by solving Maxwell's equations. In the present work the electromagnetic problem is solved using the Finite Differences Time Domain (FDTD) method implemented in CEM One R , a commercial software provided by ESI Group, then computing the discrete Fourier transform of the result to obtain the solution in the frequency domain.
The fundamentals of radar technology is nowadays used in many driver assistance systems, needing for a precise quantification of performances and limitations [1]. Its design and use need addressing many topics, raging from waveform design, propagation and pattern recognition [2][3][4], with challenges and opportunities that are driving their evolution [5]. However, its use is not limited to driver assistance, this technology is nowadays largely employed in many domains like short-range localization [6], monitoring worker activity [7], aids for visually impaired people [8], physiological monitoring [9], non-contact identity authentication [10], among many others. Signal pollution is not excluded and sometimes it must be repaired before making use of it [11].
In outdoor applications, the transmitter and antenna are placed behind a radome, a shell made of a dielectric material designed to protect the electric components against the weather. Although dielectric materials are radio-wave transparent, they have an altering effect on the electromagnetic field and must therefore be taken into account when analyzing the data. Simulation of radio wave propagation through the radome provides knowledge of the latter's precise influence on the system accuracy and performances. Thus, an important issue concerns the effect of bumpers on the radar performances, topic that attracted great interest [12][13][14][15][16][17]. In the just referred works the coupling radar-radome effects (e.g., attenuation, signal pollution, . . . ) was investigated but a parametric numerical modeling was not considered.
In our work we address precisely the analysis of that coupling, depending on the radar orientation, however because of the high frequency, the spatial mesh resolution for describing the solution of Maxwell equations leads to extremely fine meshes, with the consequent impact on the computing time. Simulation details are included in Appendix B.
To alleviate that issue, we investigate the use of model order reduction (MOR) techniques that successfully accomplished numerous parametric studies in other engineering domains. Model order reduction has been successfully applied in problems involving dynamics and waves within the so-called projection formulation, needing for a certain degree of intrusiveness when using commercial softwares. Thus, the radial approximation proposed in [18] was extended to mid-frequency dynamics within the so-called variational theory of complex rays [19]. In [20], a parametric solution of the Helmholtz equation was successfully obtained using the usual rank-one greedy PGD constructor. In [21] authors proposed a consistent reduced bases interpolation.
Non-intrusive formulations were proposed for constructing parametric solutions from a number of high-fidelity simulations performed for different choices of the model parameters, while trying to reduce as much as possible the size of the sampling for addressing multi-parametric models [22,23]. In those circumstances, usual surrogate models exhibit limitations when the number of parameters increases, and alternative technologies combining two main ingredients, the separation of variables and sparse sampling and approximations, appeared and proved their performances [22,24].
When these techniques were employed in radar engineering different difficulties appeared. The first related to the fact that the sampling scales with the characteristic length of the solution, that is with the wavelength, extremely small. On the other hand interpolation of complex-valued electric and magnetic fields produces spurious solutions. Thus, for example, the average of 1 + 0i and −1 + 0i results 0 + 0i, even if one is expecting having 0 + i. This limitation was solved by employing and alternative formulation based on the amplitude and phase. For example, in the scenario just described, the average of 1|0 and 1|π results 1|π/2, result that seems more physically consistent.
However, the use of an amplitude/phase description faces the difficulty related to its 2π periodicity, and the associated spurious discontinuities found when combining a phase close to 2π, e.g., 2π − θ (θ > 0, very small), with another very close too, e.g., 2π + θ (θ > 0, very small), but that being higher than 2π reduces to θ , originating the just referred spurious discontinuity: 2π − θ → θ . Even if this issue was addressed in many works that proposed the use of the so-called unwrapping, usual unwrapping algorithms only performs well when the data sampling is dense enough, but they fail in the sparse sampling case [25][26][27] .
The present paper proposes a technique able to conciliate unwrapping and sparse sampling, technique that constitutes the most important scientific achievement of the preset work.
Before introducing the problem statement in Section 2, the proposed strategy in Section 3 and proving its performances in Section 4, we are revisiting the main key features of non-intrusive PGD-based model order reduction in the last part of the present introduction section.

Model Order Reduction
A generic problem in physics consists of a differential operator L(·) acting on the so-called unknown field (here without loss of generality assumed scalar) u(x, t), with x ∈ Ω x ⊂ R 3 and t ∈ Ω t ⊂ R. The problem can involve a series of parameters grouped in vector µ (µ ∈ Ω µ ⊂ R P ) that being part of the optimization process, one would like to retain explicitly in the problem solution, i.e., u(x, t; µ).
The separated representations at the heart of the so-called proper generalized decomposition (PGD), initially proposed in the space-time settings [18] to generalize the proper orthogonal decomposition (POD) [28], were extended to represent parametric solutions of parametrized problems [29] u(x, t, µ 1 , . . . , enablig simulation, optimization, inverse analysis, uncertainty propagation and simulation-based control, all them under real-time constraints [28,29]. However, using such a separated representation becomes very intrusive in the sense that specific algorithms are needed for sequentially compute the different functions involved in (1) by solving 3D partial differential equations (PDE) for computing functions X i (x), 1D ordinary differential equations (ODE) for calculating the time functions T i (t) and a series of algebraic problems for obtaining the functions depending on the parameters M j i (µ j ). To reduce the intrusiveness, different strategies were proposed, among them the so-called Sparse Subspace Learning (SSL) [22] and its sparse counterpart, the so-called sPGD [24] revisited below.
If we consider the parameters µ 1 , . . . , µ P involved in the model grouped in vector µ ∈ Ω µ , the sampling constituting the DoE results in the parameters choices µ k , k = 1, . . . , K, trying to cover as much as possible the parametric domain Ω µ , to limit as much as possible extrapolation. Now, from the computed solutions at µ k , noted by u k (x, t), one could try to infer the solution for any other value of the parameter µ by using a simple interpolation, as usually considered when constructing surrogate models, However, being the parametric space in general high-dimensional standard interpolation, e.g., kriging, radial or polynomial approximation bases, ... represented by the functions N k (µ), fail for accomplishing an accurate enough approximation in the low-data limit, i.e., K P, limit in which even linear regressions do not work.
The sparse-PGD (sPGD) strategy proposed in [24] considers a separated approximation where standard polynomial bases are considered in F i j (µ n ), with the degree of the polynomial considered to represent F i j (µ n ) increasing with the separated representation mode i, i.e., the higher i, the higher the polynomial degree considered to represent F i j (µ n ) in order to avoid overfitting.
In order to calculate the approximation unknowns w i,n j (x, t), we consider the extended collocation form, within the sPGD rationale with δ µ k (µ) the Dirac mass located at µ k .

Parametric Electromagnetic Fields
In the present context, the transmitter emits waves at a constant frequency in a direction denoted by an angle θ between 0 • and 90 • . For each component of the electric and magnetic fields, the data generated by the simulations consists of an array of N × N θ values, N being the number of points in the spatial discretization and N θ the number of points in the angular discretization. The aim of this work is to provide a method to determine through the values of the electromagnetic field in all N discretized points in the geometry for any value of θ between 0 • and 90 • . In what follows, we will consider the problem of constructing the interpolation of the z-component of the electric field, E z , for all θ in one point of the geometry (on a grid of points will be considered later).

Real-Imaginary Interpolation Versus Magnitude-Phase Interpolation
Since the solutions provided by the solver are in the frequency domain, they are complex-valued functions. Complex numbers can be represented either by their real and imaginary parts or by their magnitude and phase, therefore the interpolation of complex-valued functions can be computed in multiple ways which are not equivalent. The first option is to treat the real and imaginary parts as two real-valued functions, interpolate them separately and combine the results. The second option is to treat the magnitude and phase as two real-valued functions, interpolate them separately and combine the results.
The performances of these two methods depend on the considered problem. In our case, the data from the simulation have real and imaginary parts oscillating quite fast compared to the sampling frequency (Figure 1), hence they are not good candidates for interpolation. However the magnitude is much smoother which means it may give trustworthy interpolation results. The phase looks more chaotic but it is not a good indicator because of its natural 2π-periodicity that induced spurious discontinuities: the phase needs to be "unwrapped" before being interpolated.

Phase Unwrapping in Smooth Parametric Settings
Phase unwrapping is the reconstruction of a "physically-meaningful" representation of the phase (as a function of the parameter θ) by adding multiples of 2π to some of its values in order to make it a continuous function. This step is very important because it determines the number of periods between two successive values of θ, regardless of the interpolation method used. The unwrapped phase is defined as the unique representation of the phase which can lead to a correct continuous interpolation. Unwrapping has been extensively addressed, the interested reader can refer to [26,27] and the numerous references therein.
The goal is to find a sequence (k n ) 1≤n≤N θ ∈ Z N θ , such that the unwrapped phase φ ∈ R N θ verifies: where Arg is the principal value of the phase in the interval (−π, π], under some regularity constraint.
The proposed solution consists of assuming the derivative of the phase does not vary too much, or, to put it another way, that the second derivative remains small. This hypothesis discussed later leads to a minimization of the variation of the derivative using a finite differences scheme to compute sequentially the values of the unwrapped phase: Note that θ 2 must be chosen close enough to θ 1 to ensure |φ 2 − φ 1 | < π. This algorithm can also be run in descending n order, for example to avoid initialization in a noisy area.
Once applied to our problem, phase unwrapping produces a very smooth curve, which allows great performance for interpolation ( Figure 2). The functioning interval for the just proposed strategy depends on the difference of the phase derivative between two consecutive sampling points, with respect to the sampling points distance. Thus, one can expect that ∆φ · ∆θ < C, with C = π from our numerical experiments. In order to enlarge the applicability domain, an important enhancement is proposed in Section 2.

Robustness Issues
The phase unwrapping method discussed in the previous section performs well on smooth data but it lacks robustness, mainly for two reasons:

1.
Even with very regular data, the phase can have a chaotic behavior when the magnitude is close to zero, which means our hypothesis on the regularity of the phase is not valid.

2.
Since the different phase values are computed sequentially using the previous ones, local irregularities result in high global error.
However, because a failure of the method is caused by the invalidity of the hypothesis, computing the finite differences approximation of the second derivative of the phase is a good way to detect when the unwrapping fails.

Comparing the Proposed Algorithm with Traditional Unwrapping
To illustrate the advantages of the just introduced methodology, the function having as unwrapped form with k f = 3000, is considered. Figure 3 depicts the wrapped counterpart of the above function.
For unwrapping it, we compare the proposed algorithm with the standard procedure included in the commercial software Matlab, the unwrap() function. The Matlab function applies when the difference between consecutive angles is greater than π, and shifts the angles by adding multiples of ±2π until the difference becomes less than π [30]. As noticed that function achieves a good unwrapping if a fine enough mesh is used (involving more than 35 nodes). However, the Matlab function fails to perform correctly in the case of coarser meshes (less than 30 nodes).
On the other hand, the procedure proposed in this paper achieves good results with only 15 nodes, as it can be seen in Figure 4, with the associated computing time saving.
In order to check the performances when the derivatives involved in the solution increase, we consider the previous function while increasing the factor 1/k f of 20% and 50%, with the same coarse mesh consisting of 15 nodes. The comparison is shown in Figures 5 and 6. As it can be noticed, the proposed procedure seems quite robust when compared with traditional unwrapping.

Improving the Robustness of Phase Unwrapping
So far, we have solved part of the original problem: phase unwrapping can be attempted on all N points in the spatial discretization and accepted or rejected by setting a threshold for the values of the second derivative. Thus, a part of the geometry has therefore been dealt with. In this section, we will discuss how to use these solved points to correct the phase unwrapping in the rest of the geometry.

Reduced Basis of the Search Space
The geometry can be divided into subdomains inside which the variations of the electromagnetic field are small. In each subdomain, the search space is reduced using the proper orthogonalized decomposition (POD) computed with the method of the snapshots where the snapshots are the unwrapped phases from the points in the subdomain which are already solved. The local coherence of the electromagnetic field allows the POD to have a very low dimension, typically 2 or 3, denoted in the following as r.

Phase Unwrapping in the Reduced Search Space
Out of all the possible configurations of the phase, we are looking for the one closest to the reduced search space. This can be expressed as a minimization problem: where (w i ) i=1,...,r are the basis vectors of the search space and (α i ) i=1,...,r are the vector coefficients with respect to this base. The search of k can easily be limited to [k − , k + ] N θ ⊂ Z N θ , k − and k + being two integers which depend on the size of the unwrapping. However, this problem can still be very hard to solve hence we will approximate the resolution.
The proposed algorithm to work around this high complexity follows three main steps: 1.
Using only a few of the raw phase values, generate a discrete subset of the search space which is likely to be close to the optimal solution.

2.
Fit the rest of the phase values to each curve of this subset by adding or subtracting multiples of 2π.

3.
Select the configuration which allowed for the best fit.
Thus, r points (θ n 1 , . . . , θ n r ) are chosen from the θ discretization. For (k n 1 , . . . , k n r ) ∈ [k − , k + ] r , we can find the unique curve in the search space which can be fit to the r points (Arg(E n 1 ) + 2πk n 1 , . . . , Arg(E n r ) + 2πk n r ), by solving for the coefficients α i : In practice, we choose the n r -tuple (θ n 1 , . . . , θ n r ) minimizing the condition number of the linear system represented by Equation (9). Let α i : (k n 1 , . . . , k n r ) → α i (k n 1 , . . . , k n r ) be the function that associates with each configuration the vector coefficients solution to the linear system (9). We can now solve instead of (8) the following problem:

Validation of the Optimization Procedure
We consider the following phase function involving 30 points in the discretization of the θ parametric coordinate.
Once it is artificially wrapped, both the standard unwrapping algorithm and the procedure proposed in this paper based on the second derivative, fail to unwrap it.
Assuming we have computed the reduced basis consisting of the three functions {sin(4θ), sin(8θ), sin(12θ)}, we attempt to unwrap Arg e iφ(θ)+N (0,0.5) using the just described optimization procedure. The unwrapped phase is exactly the sum of φ(θ) and the added noise, as Figure 7 proves.

Convergence Analysis
Since the simulations of radar devices are too expensive, in this section we consider a simpler problem able to provide the required data for performing the convergence analysis.
Thus, we consider the Helmholtz problem with damping in a 2D rectangular domain Ω, that looks for u(x, y) verifying where Ω = [0, 9] × [0, 2], ω = 10, ν = 1 and f (x, y; θ) = e 120i(cos(θ)x+sin(θ)y) . We discretize and solve this problem using the Finite Element Method (FEM) with Lagrange P1 elements for a number of values of θ, and interpolate the solution (complex field) at θ = 3 • . Then, we compare the phase of the interpolated field to the phase of the solution obtained with the FEM, whose difference is depicted in Figure 8.
The previous figure proves the convergence of the proposed strategy when enriching the sampling. It can also be noticed that a quite sparse sampling suffices for reaching an acceptable accuracy.

Results on the RADAR Problem
The electromagnetic field on a patch antenna (8 × 16 rectangular patches) has been computed for 33 different values of θ between 0 • and 90 • , 32 of which are used to apply our method and one as reference to compare the interpolation with the simulation. The reference electromagnetic field corresponds to the parameter value θ re f = 79.8 • . The two closest values of theta used to compute the interpolation are 76.8 • and 82.4 • hence θ re f lies in a gap of size ∆θ = 82.4 • − 76.8 • = 5.6 • . The values of the electromagnetic field below 0.2 V/m for the electric component (E) and 0.5 mA/m for the magnetic component (H) are not considered in this study because the signal is not considered significant enough and the phase cannot be studied because of its singularity (Appendix A). For this reason and due to the polarization of the electromagnetic field, only the z-component of the electric field E z and y-component of the magnetic field H y provide results. We calculate the absolute error between the phase of the reference field and the phase of the interpolated field and plot the mean error over each patch of the antenna using our method (Figures 9 and 10) and the classical unwrapping algorithm for comparison (Figures 11 and 12).   With the proposed method, the average error over the entire antenna is 1.6 • for the electric field and 1.4 • for the magnetic field while with the standard unwrapping algorithm, the average error over the entire antenna is 62.0 • for the electric field and 63.7 • for the magnetic field. In comparison, the average error on the phase when interpolating real and imaginary parts is 166.8 • for the electric field and 166.7 • for the magnetic field.
Our method is accurate and robust as the mean phase error is consistently under 6 • on the entire antenna for both the electric and magnetic fields, whereas real and imaginary parts interpolation completely fails when dealing with such a sparse discretization of the parameter θ.

Conclusions
In this paper, we tackled the computation of parametric solutions for electromagnetic wave propagation in radar applications. These simulations currently require hundreds of hours of computations to obtain suitable accuracy.
We have described an interpolation method for complex-valued fields which allows reducing the computational cost significantly. By exploiting the natural regularity of electromagnetic fields, we are able to retrieve highly accurate parametric solutions from a very sparse set of simulations.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Phase Singularity
The phase of a complex-valued function can have very fast variations when its magnitude is close to zero, even if the function is very regular ( Figure A1). This may induce large error in the interpolation of the phase, which cannot be easily avoided. However since the magnitude is very small, a large error on the phase does not have too much impact on the complex number itself, which can usually simply be regarded as zero. Figure A1. Example of a smooth parametric curve on the complex plane (left) and the phase as a function of the parameter (right).