Fast method based on CS

The wide-angle EM scattering problem solving by the traditional MoM can be described as a matrix equation with multiple right-hand sides

(1)\begin{equation}\qquad\qquad\qquad\qquad\qquad\qquad {\mathbf{Z}}\left[ {{{\mathbf{I}}_1}{{\mathbf{I}}_2}\cdot \cdot \cdot {{\mathbf{I}}_n}} \right] = \left[ {{{\mathbf{V}}_1}{{\mathbf{V}}_2}\cdot \cdot \cdot {{\mathbf{V}}_n}} \right],\end{equation}

in which, Z is the impedance matrix, V₁ to V_n are the excitation vectors at n different incident angles, and I₁ to I_n represent the n corresponding induced current vectors.

In the fast method, M new excitation vectors based on CS theory are constructed as

(2)\begin{equation}{\mathbf{V}}_i^{{\text{CS}}} = {c_{i1}}{{\mathbf{V}}_1} + {c_{i2}}{{\mathbf{V}}_2} + \cdot \cdot \cdot + {c_{in}}{{\mathbf{V}}_n}\;\;(i = 1,2, \cdot \cdot \cdot ,M),\end{equation}

where c_ij is the element in the measurement matrix Ф. In CS theory, a measurement matrix must satisfy the restricted isometry property (RIP) [Reference Candes, Romberg and Tao8], which ensures the accurate reconstruction of original signal. In general, Gaussian random matrix is often used as Ф.

Substituting equation (2) to equation (1), one can obtain M current vectors under the new excitations by

(3)\begin{equation}{\mathbf{Z}}\left[ {{\mathbf{I}}_1^{{\text{CS}}}{\mathbf{I}}_2^{{\text{CS}}}\cdot \cdot \cdot {\mathbf{I}}_M^{{\text{CS}}}} \right] = \left[ {{\mathbf{V}}_1^{{\text{CS}}}{\mathbf{V}}_2^{{\text{CS}}}\cdot \cdot \cdot {\mathbf{V}}_M^{{\text{CS}}}} \right].\end{equation}

Due to the linearity of the problem, ${\mathbf{I}}_1^{{\text{CS}}}$ to ${\mathbf{I}}_M^{{\text{CS}}}$ can be written as

(4)\begin{equation}{\mathbf{I}}_i^{{\text{CS}}} = {c_{i1}}{{\mathbf{I}}_1} + {c_{i2}}{{\mathbf{I}}_2} + \cdot \cdot \cdot + {c_{in}}{{\mathbf{I}}_n}\;\;(i = 1,2, \cdot \cdot \cdot ,M).\end{equation}

The M current vectors can be regarded as the results of M measurements of the original induced current vectors I₁ to I_n. If the original induced current vectors have a sparse representation, equation (4) can be described as

(5)\begin{equation} {\boldsymbol{\Phi}}\left[\mathbf{I_{1} I_{2} \cdot \cdot \cdot I_{n}}\right]^{\mathbf{T}} = {\boldsymbol{\Phi \Psi}} \left[\boldsymbol{\alpha_{1} \alpha_{2} \cdot \cdot \cdot \alpha_{N}}\right] = \left[\boldsymbol{I_{1}^{CS} I_{2}^{CS} \cdot\cdot\cdot I_{M}^{CS}}\right]^{\boldsymbol{T}},\end{equation}

in which Ψ is the sparse transform, N is the number of the basis functions, and α₁ to α_N are the projections of each column of [I₁ I₂ … I_n]^T in sparse domain. In the wide-angle EM scattering problems, FFT basis, Hermite basis, discrete wavelet transform (DWT) [Reference Chen, Sha and Wu9] basis, or other orthogonal bases are often selected as Ψ [Reference Liu, Qi, Cao, Chen, Deng, Huang and Wu7, Reference Cao, Chen and Wu10].

By utilizing the recovery algorithm (e.g., orthogonal matching pursuit (OMP) [Reference Yang and Hoog de11]), the projections can be approximated by

(6)\begin{align}& \left[ {{{{{\hat \alpha }}}_1}{{{{\hat \alpha }}}_2}\cdot \cdot \cdot {{{{\hat \alpha }}}_N}} \right]{\text{ = argmin}}{\left\| {\left[ {{{{{\hat \alpha }}}_1}{{{{\hat \alpha }}}_2}\cdot \cdot \cdot {{{{\hat \alpha }}}_N}} \right]} \right\|_L}\;\;s.t.\;\;{{\Theta }}\left[ {{{{\alpha }}_1}{{{\alpha }}_2}\cdot \cdot \cdot {{{\alpha }}_N}} \right]\nonumber \\ & \qquad = {\left[ {{\mathbf{I}}_1^{{\text{CS}}}{\mathbf{I}}_2^{{\text{CS}}}\cdot \cdot \cdot {\mathbf{I}}_M^{{\text{CS}}}} \right]^{T}},\end{align}

where Θ is the sensing matrix and Θ=ФΨ.

Then, the original induced current vectors are reconstructed by

(7)\begin{equation}{\left[ {{{{{\boldsymbol{\hat {\textbf{I}}}}}}_1}{{{{\boldsymbol{\hat {\textbf{I}}}}}}_2}\cdot \cdot \cdot {{{{\boldsymbol{\hat {\textbf{I}}}}}}_n}} \right]^{T} } = {{\Psi }}\left[ {{{{{\boldsymbol{\hat \alpha }}}}_1}{{{{\boldsymbol{\hat \alpha }}}}_2}\cdot \cdot \cdot {{{{\boldsymbol{\hat \alpha }}}}_N}} \right].\end{equation}

The computational complexity for solving equation (6) by OMP is O(nKMN), where K is the sparsity of original induced current vectors in sparse domain, and the inner products of the columns in sensing matrix and the measurement results are the dominant computational cost.

Construction scheme of measurement matrix

To reduce the complexity of recovery, one can make the sensing matrix sparse to decrease the cost of inner products. For the purpose, a novel construction scheme of measurement matrix is proposed as follows:

First, by randomly extracting P columns from the sparse transform Ψ, one can obtain ψ₁ to ψ_P, where ψ represents the column of Ψ.

Afterward, to better satisfy RIP, a linear superposition of these P vectors is implemented as

(8)\begin{equation}{{{\varphi }}_1} = {d_{11}}{{{\Psi }}_1} + {d_{12}}{{{\Psi }}_2} + \cdot \cdot \cdot + {d_{1P}}{{{\Psi }}_P}.\end{equation}

Finally, repeating the above two steps M times and a new measurement matrix [φ₁ φ₂ … φ_M]^T is established, in which

(9)\begin{equation}{{{\varphi }}_i} = {d_{i1}}{{{\Psi }}_1} + {d_{i2}}{{{\Psi }}_2} + \cdot \cdot \cdot + {d_{iP}}{{{\Psi }}_P}\;\;(i = 1,2, \cdot \cdot \cdot ,M).\end{equation}

Obviously, when an orthogonal basis is taken as Ψ, the inner products of φ_i and the (n − P) columns in Ψ that are not extracted in the first step are zero, respectively. In other words, there are only P non-zero elements in the ith row of Θ. A sparse Θ with MP non-zero elements is obtained by using the proposed measurement matrix, and the operations of inner products in recovery algorithm are significantly accelerated. Accordingly, the computational complexity of solution to equation (6) is decreased to O(ηnKMN), where η is the proportion of non-zero elements in Θ and η = P/n. Generally, P is much less than n.

Using the proposed measurement matrix and an orthogonal basis as Ф and Ψ respectively, equation (5) can be transformed as

(10)\begin{align}& {\left[ {{{{\varphi }}_1}{{{\varphi }}_2}\cdot \cdot \cdot {{{\varphi }}_M}} \right]^{T} }{{\Psi }}\left[ {{{{\alpha }}_1}{{{\alpha }}_2}\cdot \cdot \cdot {{{\alpha }}_N}} \right] = {\mathbf{S}}{{{\Psi }}^{T} }{{\Psi }}\left[ {{{{\alpha }}_1}{{{\alpha }}_2}\cdot \cdot \cdot {{{\alpha }}_N}} \right]\nonumber\\ \qquad & = {\mathbf{S}}\left[ {{{{\alpha }}_1}{{{\alpha }}_2}\cdot \cdot \cdot {{{\alpha }}_N}} \right]= {\left[ {{\mathbf{I}}_1^{{\text{CS}}}{\mathbf{I}}_2^{{\text{CS}}}\cdot \cdot \cdot {\mathbf{I}}_M^{{\text{CS}}}} \right]^{T} },\end{align}

in which, the proposed measurement matrix is expressed as the multiplication of a sparse random matrix S and the transposed sparse transform Ψ. The ith row of S has P random coefficients d_{i 1} to d_iP in the columns corresponding to the randomly extracted columns from Ψ in the ith measurement, and other entries in S are all zero.

So, equation (10) can be considered as measuring the sparse signals α₁ to α_N directly with the sparse random matrix S, which has been proved to satisfy a different form of RIP, so-called RIP(p) for p equal (or very close) to 1 [Reference Zhang, Liu and Wang12, Reference Gilbert and Indyk13]. Therefore, equation (6) can produce an accurate solution with high probability by using the proposed measurement matrix. It is interesting to note that, if the random coefficients are all set to 1 in (10), the sparse random matrix is simplified to a sparse binary matrix (SBM) [Reference Berinde, Gilbert, Indyk, Karloff and Strauss14], which consists of only 0 and 1. SBM is often applied as the measurement matrix in wireless sensor networks, since it is easy to be implemented on hardware and has low complexity [Reference Mamaghanian, Khaled, Atienza and Vandergheynst15].

Sphere

A sphere with the radius of 5 m illuminated by the plane waves of 300 MHz is considered, who contains 12,620 Rao-Wilton-Glisson (RWG) basis functions. The waves are set in the xoy plane, and the incident angle is divided into 0.1°, 0.2°, …, 360°.

FFT basis is used as the sparse transform. As is shown in Fig. 1, the similar precision can be achieved at the same number (310) of measurements by applying the Gaussian random matrix and the proposed measurement matrix respectively, while the number of extracting columns P is larger than 1250. It means that one can get a sensing matrix with non-zero elements accounting for 1250/3600 (η) in the best case. Thus, the computational cost for acquiring the projections of original induced currents in sparse domain is cut by about two-thirds (1 − η) by using the proposed measurement matrix rather than Gaussian random matrix; meanwhile, the computational complexity of measurement for both is the same. The comparison of the computing time for acquiring the projections is presented in Table 1, which further proves the high efficiency.

Figure 1. Relationship between the recovery error and the number of measurements in the case of FFT basis.

Table 1. Computing time of recovery for sphere (unit: s)

To show the universality of the proposed technique for orthogonal bases, FFT basis is, respectively, replaced by Hermite basis and DWT basis. The corresponding results are shown in Figs. 2 and 3. It is evident that, for both Hermite basis and DWT basis, much less computational complexity for acquiring the projections is available with the proposed measurement matrix than with Gaussian matrix under the same condition of measurement. The comparisons of computing time in the case of two bases are also provided in Table 1.

Figure 2. Relationship between the recovery error and the number of measurements in the case of Hermite basis.

Figure 3. Relationship between the recovery error and the number of measurements in the case of DWT basis.

Missile model

Then, consider a missile model who contains 1963 RWG basis functions, and the angle increment is set as 1°. The second kind of Chebyshev basis and Legendre basis is chosen as the sparse transforms, respectively. Other experimental parameters are the same with the previous one. Both Figs. 4 and 5 indicate that the proposed scheme is also validity for the complex-shaped objects.

Figure 4. Relationship between the recovery error and the number of measurements in the case of the second kind of Chebyshev basis.

Figure 5. Relationship between the recovery error and the number of measurements in the case of Legendre basis.

By using the second kind of Chebyshev basis and the proposed measurement matrix (M = 90, P = 100), and setting the elements less than 10⁻¹⁴ in the sensing matrix to zero, there are only 9000 non-zero elements in the sensing matrix, which is consistent with the expected number MP. Hence, the solution to equation (6) with the help of OMP is accelerated, which is demonstrated in Table 2. The bistatic radar cross section (RCS) of the missile model illuminated at a random incident angle (take 77° as an example) is also provided in Fig. 6, which agrees well with the result solved by the traditional MoM. From Table 2 and Fig. 6, we can see clearly that the computing time can be significantly reduced while the high accuracy is kept by using the proposed measurement matrix.

Figure 6. Comparison of RCS for the missile model illuminated at a random incident angle.

Table 2. Computing time of recovery for missile model (unit: s)