Radar Scattering Analysis of Multi-Scale Complex Targets by Fast VSBR-MoM Method in Urban Scenes

Abstract

An innovative and efficient hybrid technique, which combines the Method of Moments (MoM) with Volume Meshed Shooting and Bouncing Ray (VSBR), is presented to analyze the scattering of metallic–dielectric mixed multi-scale structures in urban scenes. Additionally, the technique can rapidly generate radar images at different angles, which is useful for the video remote sensing community. By dividing the mixed multi-scale targets into sub-regions, different solvers are employed to compute the scattering contributions based on their varying electrical sizes. For large-scale sub-regions, the VSBR method is employed on both the medium and metal parts, leading to a multilevel electromagnetic field including both the direct induced field and the multi-reflection field. This field contributes to the integral equation in the MoM sub-regions. Additionally, interactions from the MoM region are considered within the VSBR region, followed by a surface current integration to compute the scattered field. When addressing mixed multi-scale electromagnetic problems, this technique proves to be more efficient and easier to implement in general-purpose computer codes. Furthermore, this method is accelerated using the local coupling (LC) theory and fast multipole method (FMM). Using this fast computational method, efficient simulation results of radar images at different angles for the scenes can be obtained and the numerical results demonstrate the efficiency and accuracy of the proposed method.

1. Introduction

The importance of electromagnetic situational awareness in complex urban scenes is increasingly significant, particularly in terms of radar imaging simulation. With the continuous advancement of manufacturing technology, urban scene simulation is facing increasing challenges with multi-scale problems. Thus, accurate modeling and analysis of these multi-scale problems are of significant importance. Scholars around the world have been dedicated to developing innovative methods and techniques to tackle these challenges in engineering practice. The Institute of Electrical and Electronics Engineers Transactions on Antennas and Propagation (IEEE TAP) and the Proceedings of the IEEE (PIEEE) have separately published Special Issues highlighting the latest developments in computational electromagnetics: No. 8 of 2008 [1] and No. 2 of 2013 [2,3]. The former issue concentrates on “Large and Multiscale Computational Electromagnetics”, while the latter predominantly explores multi-scale problems. Electromagnetic scattering modeling for complex large-scale objects is considered essential for applications including vehicle identification, recognition, signal communication, and flight control [4,5]. Due to the prohibitive computational complexity of full-wave methods like the integral equation method, high-frequency approximation methods are favored for their efficiency and accuracy in addressing large-scale scattering problems.

Common high-frequency methods include Physical Optics (PO) [6,7,8], Geometrical Optics (GO) [9,10,11], the Geometrical Theory of Diffraction (GTD) [12,13,14], and the Uniform Theory of Diffraction (UTD) [15,16,17]. Despite the rapid nature of these algorithms, they are prone to significant errors when applied to electromagnetic problems in complex systems, consequently limiting their application. The integration of high-frequency methods with other techniques is frequently observed in current research, particularly with low-frequency numerical approaches. For instance, C.F. Wang accelerated the MoM-PO method via GPU [18], enabling a fast solution for antenna issues on electric power carriers. In addition, J. F. Lee and R. Rernadez-Recio investigated the fusion of FEM and PO [19,20]. Zhang implemented an iterative MoM-PO technique to explore the radiation characteristics of multiple antenna systems on electric power carriers. Liu introduced a hybrid method that merges high-order MOM with PO, facilitating flexible segmentation while maintaining precise modeling. To date, comprehensive methods for computing electromagnetic scattering in complex metallic–dielectric mixed multi-scale targets in urban scenes have been insufficiently explored. Actually, structures such as drones flying at low altitudes, buildings in intricate environments and antenna radomes on targets often involve mixed materials that are non-homogeneous or lossy. Consequently, addressing these challenges efficiently calls for further investigation.

Alternatively, the difficulty in information perception and transmission has increased due to the concentration of the population and the presence of tall buildings in urban scenes [21]. Target recognition and detection also become more difficult in urban scenes in the video remote sensing community. To address this situation, it is essential to start with the accuracy of electromagnetic perception, build a highly realistic urban electromagnetic model, and achieve transparency in urban electromagnetic information. The main content of the research includes constructing the geometric features of urban environments such as trees and buildings. Through the ray tracing method and lossy medium electromagnetic wave propagation theory, an electromagnetic characteristic model of urban scenes and artificial targets is established, enabling the characterization of the scattering characteristics distribution in large-scale scenarios.

In this study, an iterative fast VSBR-MoM method is proposed, is introduced, utilizing both MoM and VSBR. Initially, the complex urban scenes include trees, buildings, rough surfaces, and multi-scale targets. These targets are separated into VSBR and MoM regions according to their structural characteristics. The VSBR region includes smooth metal surfaces and dielectric volumes, while the MoM region comprises more intricate structures. The first-order scattering fields of these regions are determined independently. Subsequently, interactions between the first-order scattering fields of the VSBR and MoM regions occur through current radiation and field reflection, leading to the generation of second-order and third-order coupled fields via iterative processes. This process continues until the interactions become negligible. To enhance efficiency, the local coupling theory and the FMM theory are employed to accelerate current interactions between distant surface elements, significantly reducing the computational demands of simulations. Section 2 outlines the formulation for calculating metallic–dielectric mixed and multi-scale targets. Section 3 presents the proposed fast iterative solution process for coupled fields and discusses its convergence control. The accuracy and efficiency of the proposed method have been confirmed through numerical results. The main contributions of this study are presented as follows:

(1)

Based on a volume meshed ray tracing method for electromagnetic wave propagation in lossy medium, this work analyzes building and tree models in complex urban scenes. Compared with traditional numerical methods, our method significantly enhances computational efficiency.

(2)

The proposed VSBR-MOM hybrid method addresses the resource consumption issues encountered when solving multi-scale problems in urban scenes with the local coupling theory and FMM theory. The local coupling theory is applied to reduce the coupling computation unknowns between high-frequency and low-frequency calculation regions, while the fast multipole method decreases the complexity of the coupled field solution.

(3)

It realizes the efficient computation of electromagnetic scattering from rough surfaces, buildings, trees, and multi-scale targets in complex urban environments, as well as radar imaging simulation of typical scenes from different angles in the video remote sensing community. Compared with traditional high-frequency imaging methods, the images obtained using the method in this paper exhibit more detailed scattering characteristics.

2. Formulations

2.1. Overview

The urban environment includes trees, buildings, rough surface areas, and multi-scale structures. The fast VSBR-MOM method is used to calculate the electromagnetic scattering properties of complex targets in an urban scene. Furthermore, it improves the efficiency of radar imaging simulations. In this section, the complete procedures of the hybrid technique are thoroughly detailed. The buildings are modeled as a non-uniform medium, trees as a combination of non-uniform medium and thin medium, and rough surfaces can be considered reflective surfaces as illustrated in Figure 1a. For these scenes, the volumetric meshed ray-tracing method VSBR is used to compute the scattered field. Complex targets in an urban scene are metallic–dielectric mixed multi-scale models as illustrated in Figure 1b.

In urban scenes, the non-uniform medium, rough surfaces, and metallic–dielectric mixed multi-scale scattered are divided into two regions. This work introduces ray tracing using tetrahedral mesh to compute the scattering from large-sized elements. Then, it describes the process of using the current integration to compute smaller-scale regions. Subsequently, it presents the coupled-field solution method and the acceleration process.

2.2. The VSBR Application in Buildings, Trees, and Large-Scale Regions in an Urban Scene

As shown in Figure 1a, the buildings, trees, and rough surfaces in a complex urban scene are categorized into a large-scale region and analyzed using the VSBR method. A dense grid of ray tubes, representing the incident field, is directed toward the target. For the rough surface part, only ray tracing on the meshed surface is required to calculate the reflection effect. However, for the building and tree parts, volume meshing plays an essential role in solving both the transmission of electromagnetic waves within the medium and the reflection at its surface. In this approach, the VSBR method differentiates between metallic and dielectric regions.

On rough surfaces, when rays are incidentally on a medium, they undergo reflection at the surface and transmission inside it. According to Fresnel’s formula, the reflection and transmission coefficients for a lossless rough surface can be calculated with the following Formulas (1) and (2):

$$R_{\parallel }={\displaystyle \frac{{\epsilon }_2{k}_{1z}-{\epsilon }_1{k}_{2z}}{{\epsilon }_2{k}_{1z}+{\epsilon }_1{k}_{2z}}},\begin{array}{c}{T}_{\parallel }={\displaystyle \frac{2{\epsilon }_2{k}_{1z}}{{\epsilon }_2{k}_{1z}+{\epsilon }_1{k}_{2z}}}\end{array}$$

(1)

$$R_{\perp }={\displaystyle \frac{{\mu }_2{k}_{1z}-{\mu }_1{k}_{2z}}{{\mu }_2{k}_{1z}+{\mu }_1{k}_{2z}}},\begin{array}{c}{T}_{\perp }={\displaystyle \frac{2{\mu }_2{k}_{1z}}{{\mu }_2{k}_{1z}+{\mu }_1{k}_{2z}}}\end{array}$$

(2)

where ${\epsilon }_1$${\epsilon }_2$ and ${\mu }_1$${\mu }_2$ are the dielectric constant and permeability of two contacting mediums (in this case it represents air and ground surface materials). $k_{2z}$ and $k_{1z}$ denote the wave number component perpendicular to the interface. $\perp$ and $\parallel$ represent the vertical and parallel polarizations.

For buildings and trees, the mediums considered here exhibit lossy properties. This results in energy loss as electromagnetic waves propagate through it. As shown in Figure 2a, when a non-uniform plane wave enters such a medium, the refractive index becomes complex, incorporating both real and imaginary components. The imaginary part signifies energy dissipation owing to absorption and scattering within the medium. A non-uniform wave with amplitude vector ${\alpha }_1$ and phase vector ${\beta }_1$ moves from medium 1 to medium 2 across the interface [22]. ${\alpha }_2$ describes the transmitted wave’s amplitude vector, while ${\beta }_2$ represents its phase vector. Therefore, it is assumed that the incident wave has a complex propagation vector ${\gamma }_1={\alpha }_1+j{\beta }_1$. ${\theta }_1$ is the angle between ${\alpha }_1$ and ${\beta }_1$. When the incident wave passes, the angle between the interface normal vector $\widehat{n}$ and phase vector ${\beta }_1$ denoted by ${\phi }_1$. Similarly, in medium 2, the angle between the phase vector ${\beta }_2$ and the normal vector $\widehat{n}$ is ${\phi }_2$. Meanwhile, the angle between ${\alpha }_2$ and ${\beta }_2$ is ${\theta }_2$.

As established in paper [23], reflection occurs with an angle symmetric to the incident angle, mirroring it precisely. The transmission vector ${\alpha }_2$ and ${\beta }_2$ are derived from boundary conditions in relation to the incidence angle and the dielectric constants. Once these angles are determined, the corresponding reflection and transmission coefficients can also be calculated in (3) and (4).

$$\begin{array}{l}{R}_{\perp }={\displaystyle \frac{{\mu }_2({\alpha }_1\cos {\zeta }_1+j{\beta }_1\cos {\phi }_1)-{\mu }_1({\alpha }_2\cos {\zeta }_2+j{\beta }_2\cos {\phi }_2)}{{\mu }_2({\alpha }_1\cos {\zeta }_1+j{\beta }_1\cos {\phi }_1)+{\mu }_1({\alpha }_2\cos {\zeta }_2+j{\beta }_2\cos {\phi }_2)}}\\ T_{\perp }={\displaystyle \frac{2{\mu }_2({\alpha }_1\cos {\zeta }_1+j{\beta }_1\cos {\phi }_1)}{{\mu }_2({\alpha }_1\cos {\zeta }_1+j{\beta }_1\cos {\phi }_1)+{\mu }_1({\alpha }_2\cos {\zeta }_2+j{\beta }_2\cos {\phi }_2)}}\end{array}$$

(3)

$$\begin{array}{l}{R}_{\parallel }={\displaystyle \frac{-{\epsilon }_2({\alpha }_1\cos {\zeta }_1+j{\beta }_1\cos {\phi }_1)+{\epsilon }_1({\alpha }_2\cos {\zeta }_2+j{\beta }_2\cos {\phi }_2)}{{\epsilon }_2({\alpha }_1\cos {\zeta }_1+j{\beta }_1\cos {\phi }_1)+{\epsilon }_1({\alpha }_2\cos {\zeta }_2+j{\beta }_2\cos {\phi }_2)}}\\ T_{\parallel }={\displaystyle \frac{2{\epsilon }_2({\alpha }_1\cos {\zeta }_1+j{\beta }_1\cos {\phi }_1)}{{\epsilon }_2({\alpha }_1\cos {\zeta }_1+j{\beta }_1\cos {\phi }_1)+{\epsilon }_1({\alpha }_2\cos {\zeta }_2+j{\beta }_2\cos {\phi }_2)}}\end{array}$$

(4)

Eventually, the reflection coefficient and transmission coefficient of horizontal and vertical polarized waves of a non-uniform plane wave in a lossy medium are solved. When calculating the scattering characteristics of leaves, traditional ray methods need to model the leaves with thickness. As shown in Figure 2b, the incident wave is reflected and transmitted when the upper surface of the blade is irradiated. The electromagnetic wave transmitted through the upper surface will continue to be reflected and transmitted through the lower surface. In our theory, part of the energy passes through the blade and propagates outward, while the other part of the energy is still left, spread inside the leaves. In this study, due to the huge number of blades, we simplified the blades. We recorded the transmitted rays of electromagnetic waves penetrating the blades, considering that only one reflection occurred here, and compensated for the phase change $r_{\mathrm{delay}}$ when calculating the transmitted field. $r_{\mathrm{delay}}$ is related to the angle between the incident wave and the blade. The refraction angle is ${\theta }_i^1$ and the direction is ${\widehat{k}}_{\mathrm{i}}^1$. The direction of the reflected rays on the surface 1 of the leaf is calculated according to Formulas (5) and (6). The phase difference $r_{\mathrm{delay}}={\displaystyle \frac{d}{\sin \theta }}$ generated thereby is also considered in the calculations.

$${\widehat{k}}_r={\widehat{k}}_{\mathrm{i}}^1-2 ({\widehat{k}}_{\mathrm{i}}^1\cdot \widehat{n})\widehat{n}$$

(5)

Then, the 1st bouncing transmission direction vector is:

$$s_t^1=v\times \sin ({\theta }_t^1)-\overrightarrow{n}\times \cos ({\theta }_t^1)$$

(6)

Taking the transmission direction of the upper leaf surface as the incident direction of the lower leaf surface, the direction vector of the emitted rays can be obtained as:

$$s_t^1=s_{\mathrm{i}}^2$$

(7)

$$v_{\mathrm{i}}^2=s_{\mathrm{i}}^2\times {\widehat{n}}^2$$

(8)

$$v_t=\widehat{n}\times {\mathbf{v}}_{\mathit{i}}^2$$

(9)

$$s_t^2=v_{\mathrm{t}}\times \sin ({\theta }_t^2)-\widehat{n}\times \cos ({\theta }_t^2)$$

(10)

By adopting this equivalent method, the number of unknowns in the calculation of leaves can be significantly reduced.

2.3. The MOM Application in Multi-Scale Structure in an Urban Scene

When solving multi-scale targets in an urban scene, we separately use the MoM and the VSBR method to calculate large-scale and small-scale regions. The multi-scale model of a drone in an urban scene is displayed in Figure 3. $E^{in}$ denotes the incident EM field. $E_{VSBR}^{direct}$ and $E_{\mathit{MoM}}^{direct}$ represent the direct scattering field from the VSBR and MOM regions. $E_{VSBR\to \mathit{MoM}}^{couple}$ and $E_{\mathit{MoM}\to VSBR}^{couple}$ are the current radiation and multiple reflection field between the VSBR region and the MoM region. In the composite model, regions with intricate PEC structures, such as the rotor blades in the drone shown in Figure 3, are classified into the MoM region. When the incident wave strikes the model, an initial induced current $J_0^{MoM}$ is generated within the MoM region. By disregarding the coupling field, $J_0^{MoM}$ can be solved by forming a matrix equation derived from the electric field integral equation (EFIE):

$$\left[Z_{mn}^{MoM}\right]\left[I_n^{MoM}\right]=\left[V_m^{inc}\right]$$

(11)

When the scattering field in the VSBR region interacts with the MoM region, a new integral equation is formulated. In this scenario, the surface current $J_0^{VSBR}$ serves as the excitation source to produce the scattered electric field within the MoM region.

2.4. The Iterative VSBR-MoM Application in Two Regions

The iterative VSBR-MoM algorithm includes the process of incident wave irradiation, iterative coupling, and convergence judgment in the overall framework as shown in Figure 4. Ray tracing in the VSBR method is divided into one-bounce reflection and multi-bounce reflection. The tracing of a ray tube can be concluded when rays are emitted from the medium to the air or when energy attenuation falls below a specified threshold.

Subsequently, the scattering field or coupling field produced by the ray tube can be determined. For ideal conductors, surface currents are computed directly. By contrast, for dielectric materials, the calculation of the surface equivalent magnetic current through the ray tube is necessary. This involves determining both the equivalent current density $J_{\mathrm{s}}$ and effective magnetic flux density $M_{\mathrm{s}}$. The far-scattered electric field resulted from aggregating the fields produced by these equivalent electromagnetic flows. Figure 4 illustrates the i-th coupling within the MoM-VSBR process. The total field-induced surface current within the VSBR region acts as the excitation source, facilitating the generation of the couple field in the MoM region. Taking the surface currents as an example,

$$E_1^{MoM\leftarrow VSBR}\left(r\right)={\displaystyle {\int }_s-jk\eta \overline{\overline{G}}\left(r,r^{\prime }\right)\cdot J_0^{VSBR}\left(r\prime \right)dS\prime }$$

(12)

By combining the coupling field with the initial incident field to form a new incident field in both regions, the scattering field for the first coupling is determined.

$$\left[Z_{mn}^{MoM}\right]\left[I_{n,1}^{MoM}\right]=\left[V_m^{inc}\right]+\left[〈E_1^{MoM\leftarrow VSBR},f_m〉\right]$$

(13)

where $I_{n,1}^{MoM}$ is the unknown coefficients for the MoM current, $Z_{mn}^{MoM}$ stands for the impedance matrix in the MoM region, $V_m^{inc}$ represents the excitation matrix by the direct incident field, and $〈E_1^{MoM\leftarrow VSBR},f_m〉$ denotes the coupling interaction from the VSBR region to the MoM region.

The coupling field in the VSBR region is generated by employing the MoM current as an excitation source. This process is expressed using Huygens’ principle through the following Formula (14):

$$H_1^{VSBR}\left(r\right)={\displaystyle {\int }_s\nabla \times \overline{\overline{G}}\left(r,r\prime \right)\cdot J_0^{MoM}\left(r\prime \right)d{S}^{\prime }}$$

(14)

Next, the coupling field is combined with the initial incident field to create a new incident field in the VSBR region. Then, the updated coupling current is recalculated, marking the completion of the first interaction. The process is repeated until the current stabilizes, satisfying the following condition (${\epsilon }_i$ generally is 0.01):

$${\epsilon }_i=‖\left[I_{i+1}^{MoM}\right]-\left[I_i^{MoM}\right]‖/‖\left[I_i^{MoM}\right]‖\le \epsilon$$

(15)

The distinction method is that if the mesh size suitable for VSBR, which is greater than 1/8 wavelength, can fully construct the shape of the target, it is considered a large-scale region. If the shape cannot be reconstructed, it is considered a small-scale region.

2.5. The Acceleration of Fast VSBR-MoM Algorithm

The coupling effect of scattering fields between buildings, trees, rough surfaces, and multi-scale targets in complex urban scenes is limited by computational efficiency. At first, the large area of the ground surface, the large size of buildings, and the numerous leaves on trees contribute to the complexity. Secondly, the multi-scale structure of complex targets leads to fine details of the targets. This results in a very large total computational load. Therefore, the work in this paper requires us to reduce the unknown quantities in electromagnetic scattering calculations and lower the computational complexity, in order to achieve fast multi-view imaging simulation of complex scenes. This necessitates the use of a local coupling iterative approach to modifying the way global variables interact, dividing the regions with strong coupling interactions. Meanwhile, the coupling computation process can also be accelerated using the fast multipole method (FMM).

2.5.1. Local Coupling on the Ray Tracing Path to Reduce the Unknown Quantities

In this part, the coupling effects between the VSBR region and the MoM region are calculated using an iterative method. These effects are derived from the integral process of mutual surface currents. A truncation method is employed [24] to concentrate on retaining only the major contribution energy from the main lobe or side lobes with a relatively larger amount of energy. This process is known as the local coupling approach. Then, the process begins with ray tracing in both the VSBR and MoM regions, recording mirror reflection directions and intersection surface elements, such as the red triangle surface element and the ray path shown in Figure 5. Once the path is determined, adjacent boxes within a certain range can be identified based on the Octree box code where the interface elements reside. The current interactions of the interface elements within the box are then calculated to resolve the local coupling field.

Formerly, for VSBR-MoM method, each $VSBR\to M\mathrm{o}M$ iteration requires evaluation of multiple field integrals over $N_{VSBR}$ elements of the source distribution at $N_{MoM}$ observation elements. There is an $O\left(N_{VSBR}\times N_{MoM}\right)$ computational complexity. Each $M\mathrm{o}M\to VSBR$ iteration requires evaluation of multiple field integrals over $N_{MoM}$ elements of the source distribution at $N_{VSBR}$ observation elements. The total number of unknown variables is $N_{MoM}\times N_{VSBR}$. For the enclosed targets, we set the number of iterations as the same as the number of reflections that occur [25], ensuring that the number of current iteration times is consistent with the number of reflections and that an acceptable accuracy is obtained.

Now, the coupling effect calculation process of the proposed acceleration theory consists of two parts: ray tracing and current iteration. The computational complexity of the ray tracing part is consistent with that of the ray tracing method. For the current iteration part, we select $N_{\mathit{sel}}$ elements instead of $N_{VSBR}$ elements to calculate the current re-radiated field to one observation intersected elements on one ray path. Thus, the number of operations that occur in the current interactions is related to $N_{\mathit{sel}}\times N_{MoM}$. We defined the unknowns of the fast VSBR-MoM method as $N_{VSBR-MoM}$ and $N_{VSBR-MoM}=N_{\mathit{sel}}\times N_{\mathit{MoM}}$. $N_{sel}=N_{VSBR }/8^{T_{select}}$. $T_{select}$ is the level of the octree boxes we choose to calculate the triangle element currents in it. We can easily observe that $N_{\mathit{sel}}\times N_{\mathit{MoM}}$ is much smaller than $N_{MoM}\times N_{VSBR}$.

2.5.2. FMM Acceleration

For elements in the MoM region and the VSBR region, when their coupling is directly calculated, each element acts as a scattering center. By employing the FMM [26], interactions within the far zone group are computed using aggregation, translation, and disaggregation operators, as depicted in Figure 6, which significantly reduces the number of radiation current sources. In conjunction with FMM theory, the dyadic Green’s function can be expressed in the following form (16).

$$\begin{array}{l}{\displaystyle \frac{e^{-jk\left|r+d\right|}}{4\pi \left|r+d\right|}}={\displaystyle \frac{-jk}{16{\pi }^2}}{\displaystyle ∯e^{-jkd}{\displaystyle \sum _{l=0}^{\infty }{(-j)}^l}\left(2l+1\right)h_l^{\left(2\right)}\left(kr\right)P_l\left(\widehat{k}\cdot \widehat{r}\right) }{d}^2\widehat{k}\\ \approx {\displaystyle \frac{-jk}{16{\pi }^2}}{\displaystyle ∯e^{-jkd}{\displaystyle \sum _{l=0}^L{(-j)}^l}\left(2l+1\right)h_l^{\left(2\right)}\left(kr\right)P_l\left(\widehat{k}\cdot \widehat{r}\right)} d^2\widehat{k}\end{array}$$

(16)

The schematic diagram illustrating the aggregation, transfer, and configuration of the two elements is presented in Figure 6. The field point is $r_m$, the source point is $r_{\mathrm{n}}$, and the group center points are $r_p$ and $r_q$.

$$\begin{array}{l}\begin{array}{ccc}{r}_{mp}=r_m-r_p& r_{pq}=r_p-r_q& r_{qn}=r_q-r_n\end{array}\\ r_{mn}=r_{mp}+r_{pq}+r_{qn}\end{array}$$

(17)

When $\left|r_{mp}+r_{qn}\right|<\left|r_{pq}\right|$, then,

$$\nabla G({\overrightarrow{r}}_m,{\overrightarrow{r}}_n)={\displaystyle \frac{-k^2}{16{\pi }^2}}{\displaystyle ∯\widehat{k} e^{-jk\cdot (r_{mp}+r_{qn})}{\alpha }_{pq}(\widehat{k}\cdot r_{pq}) d^2\widehat{k}}$$

(18)

$$\overline{\overline{G}}(r_m,r_n)={\displaystyle \frac{-jk}{16{\pi }^2}}{\displaystyle ∯ (\overline{\overline{I}}-\widehat{k}\widehat{k})e^{-jk\cdot (r_{mp}+r_{qn})}{\alpha }_{pq}(\widehat{k}\cdot r_{pq}) d^2\widehat{k}}$$

(19)

The angular spectrum expansion of Green’s function is incorporated into the accelerated method to calculate the composite scattering of rough surfaces and targets. The induced electromagnetic currents on the surfaces of the VSBR and MoM regions are iteratively updated until they reach stability.

3. Numerical Results

In the present section, several numerical examples are simulated to illustrate the accuracy and efficiency of the proposed method. In this study, the computations in this section are carried out on the computer of Intel Xeon E7-4850 CPU equipped with 36GB RAM. This paper provides a validation of the correctness of the simulation methods for building models, tree models, and complex structures in urban scenes, and discusses the improvements made in computational efficiency. Ultimately, it achieves fast imaging simulation of urban scene models.

3.1. Scattering Simulation for Building Models

For buildings in urban scenes, we should model them as a combination of multiple non-uniform cubes with walls and windows. This work simulates two cases of cube models, as shown in Figure 7a,c below, presenting the Bistatic RCS results based on our method. By comparing the commercial software FEKO with the FDTD solver in Figure 7b,d, it can be found that the method proposed in this work can accurately calculate the scattering fields of the building models. The simulation frequency is set to 3 GHz, with a bistatic incident angle of $\theta ={45}^{\circ }$ and $\phi =0^{\circ }$. The computation times for the methods presented in this work and for the commercial software are summarized in

Table 1. Calculation time and memory use of our method and the FEKO-FDTD solver.

Solver		FEKO FDTD	Our Method
Case 1	Calculation time	21 min	36 s
	Memory	1.62 GB	502 MB
Case 2	Calculation time	24 min	43 s
	Memory	1.78 GB	621 MB

, demonstrating that our method significantly outperforms the FDTD method of commercial software in terms of computational efficiency. In case 1 and case 2, FEKO’s numerical algorithm FDTD takes an average of 34 times longer for computation compared to the method presented in this work. Regarding memory consumption, the method presented in this work uses only 1/3 of the memory required by the precise algorithm.

3.2. Scattering Simulation for Tree Models

For trees, our work adopts the one-bounce ray tracing approach. Firstly, we need to demonstrate that for thin dielectric leaves, using one-bounce reflection and transmission can simulate the scattering results of electromagnetic waves undergoing three reflections and transmissions on the leaves and that this simulation is a valid approximation. Taking the leaf model shown in Figure 8a below as an example with relative permittivity ${\epsilon }_r=3+\mathrm{j}0.1$, the structure of the leaf is depicted, and the simulation results offer one-bounce and three-bounce computations for a bistatic RCS in Figure 8b. The simulation frequency is 1 GHz, and the incident angle is $\theta =-{45}^{\circ }$ and $\phi =0^{\circ }$. From the simulation results for one bounce and three bounces in this Figure, it can be observed that the proposed method with one bounce yields results that are close in accuracy to the method with three bounces. In terms of computation time, the one-bounce calculation takes 1 s, while the three-bounce calculation takes 2 s.

To further analyze the advantages of the method in calculating the scattering from leaves, the other examples compare the computation times for multiple leaves in Figure 8c,e,g. We performed monostatic scattering RCS simulations for multiple-leaf models. As shown in Figure 8d,f,h, we simulated different numbers of leaves. Since the leaves in this method are simplified as thickness-less surface patches, only the scattering field on the meshed surface is computed. The unknowns for the three models are 2270, 4506, and 6078, respectively. The relative permittivity ${\epsilon }_r=3+\mathrm{j}0.1$. The monostatic scattering field is calculated at an angle of $\theta =-{90}^{\circ }\sim {90}^{\circ }$, $\phi =0^{\circ }$ and a frequency of 3 GHz. The FEKO-MoM solver is used in these cases as a reference for accurate results. During the MoM calculation, a volumetric grid discretization is employed, considering the leaf thickness d = 0.003 m. From the comparison of RCS for the three cases shown in the figure, it can be observed that, regardless of the number of facets, the results of the proposed method and the MoM method are consistent. Below are the computational efficiency statistics for the three types of leaves.

Table 2. Calculation time and memory use of our method and the FEKO-MoM solver.

Leaves		FEKO MoM	Our Method
1270 mesh elements	Calculation time	0.18 min for 1 angle	4 s for 1 angle
	Memory	0.17 GB	102 MB
2506 mesh elements	Calculation time	1.06 min for 1 angle	5 s for 1 angle
	Memory	1.97 GB	221 MB
6078 mesh elements	Calculation time	16.2 min for 1 angle	7 s for 1 angle
	Memory	4.55 GB	521 MB

presents the mesh counts, computation time, and memory consumption for three leaf models. The computation time and memory consumption of the MoM method increased exponentially, and its computational complexity is consistent with that of $O(N^2)$. By contrast, the computation time and efficiency of the method proposed in our work increase approximately linearly. Furthermore, the computational efficiency of the method proposed in this work is significantly better than that of numerical MoM methods.

Finally, in this section, we computed a model of the tree, as shown in Figure 9. The tree was subdivided into 9388 tetrahedral elements, and the reflection and transmission of electromagnetic waves were calculated in the trunk. Attenuation was formed inside the tree. The relative permittivity of the trunk is ${\epsilon }_r=8+\mathrm{j}0.01$, and that of the leaves is ${\epsilon }_r=3+\mathrm{j}0.1$. The HH polarization monostatic scattering RCS is calculated at an angle of $\theta =0^{\circ }\sim {90}^{\circ }$ and a frequency of 3 GHz. We compared the MoM algorithm in the FEKO software with the traditional high-frequency algorithm PO. The PO method only considers the scattering effect of surface currents generated by electromagnetic waves incident on the surface, without considering the direct interaction between the target elements. The following presents a comparison of the computational errors from the three methods. From the error comparison, it can be found that the method in this paper significantly improves the calculation accuracy when compared with traditional high-frequency methods. In comparison with precise numerical methods, it also provides a clear advantage in computational efficiency.

Table 3. Calculation time, memory use and RMS error of our method PO method and the FEKO-MoM solver.

Tree	FEKO MoM	PO	Our Method
Calculation time	5.7 h	0.04 h	0.19 h
Memory	12 GB	1011 MB	1650 MB
RMS Error $\theta =0^{\circ }\sim {70}^{\circ }$	Bench mark	8.23 dB	4.71 dB

presents the statistical results for the computation time and efficiency of the entire tree model. It is clear that the method proposed in this paper significantly outperforms in terms of computational efficiency, with a computation time 30 times faster than the numerical method (MoM), and a memory consumption ratio of 7.5:1. Although the PO method also has the advantage of shorter computation time and lower memory consumption, its computational accuracy is clearly lower than that of the method presented in this paper, making it unsuitable for calculating models in complex urban scenarios.

3.3. Scattering Simulation of Metallic–Dielectric Mixed Multi-Scale Models

In this section, a metallic–dielectric mixed SLICY model operating at a frequency of 9.4 GHz is analyzed to verify the validity of our method. The orange cylinder in Figure 10a is PEC and designated as the MoM region. The remaining parts are mediums and belong to the VSBR region. The observation angle is $\phi =0^{\circ }$ and $\theta =0^{\circ }~{360}^{\circ }$. The MLFMM is used as a reference to provide precise numerical simulation results [27], allowing for validation of the correctness of the method proposed in this paper through comparison. The results shown in Figure 10b suggest that, at most angles, the two curves align quite well, particularly near the specular reflection angle. The PO method is used to demonstrate that it is unsuitable for calculating such metal-dielectric hybrid models.

Next, the monostatic scattering field for a drone model (length 1.5 m, width 0.9 m, and height 0.3 m) has been calculated at 3 GHz. Four different cases were adopted for dividing the regions (MoM and VSBR), as illustrated in Figure 11a,b. The drone’s fuselage serves as a large-scale region, utilizing the ray tracing computation of VSBR. The rotors and landing gear, which are made of metal, are divided into different MOM and VSBR regions in the four cases. The monostatic observation angle is $\phi =90^{\circ }$ and $\theta =0^{\circ }~{180}^{\circ }$. The fuselage part is medium with ${\epsilon }_{\mathrm{r}}=1.5-j0.015$. By comparing the VV polarization results from these four different division scenarios, a deeper understanding of the techniques involved in using the hybrid method was achieved. Figure 11b displays the four different cases for dividing the drone model, and Figure 11c presents the calculation results for these cases.

Table 4. Calculation time of our method in different processes.

Case	Region Unknown		Ray Tracing Time in VSBR and Coupling Process	MoM Region Calculating Process	Interacting Time of Coupling Effect	Total Time for 181 Angles	Quantitative Evaluation
	MoM	VSBR
1	3082	150,320	6.2 s	52 s	16 s (Local current + FMM)	3.4 h	1:3.4
					175 s (no local current and no FMM)	11.6 h
2	4157	142,340	6.8 s	97 s	20 s (Local current + FMM)	5.58 h	1:4.2
					366 s (no local current and no FMM)	23.2 h
3	5054	131,062	7.0 s	215 s	23 s (Local current + FMM)	12.4 h	1:3.1
					554 s (no local current and no FMM)	38.6 h
4	5920	121,223	7.2 s	421 s	25 s (Local current + FMM)	22.92 h	1:2.3
					695 s (no local current and no FMM)	53.0 h

presents the computed unknowns, ray tracing time, MoM calculation time, coupling solution time, and the time for accelerated iterative coupling for the four simulation scenarios. In the simulation of four scenarios, a gradual decrease in the size of the VSBR region is observed, while the MoM region expands concurrently. Various partitioning cases were recorded with respect to unknown quantities and computation times, allowing for a comparison of computational complexity as these unknown quantities vary across different regions. The analysis results indicate that the efficiency of calculations is determined by the choice and size of the partitions. Clearly, the computation time of the coupled solving process is greatly reduced after accelerating with LC and FMM theory.

From Figure 11c, it is observed that the RCS results for cases 1, 2, and 3 exhibit RMS errors relative to case 4, with values of 4.2 dB, 2.1 dB, and 1.1 dB, respectively. As the MoM method serves as a more accurate calculation solver, a larger region leads to increased accuracy in the fast MoM-VSBR method. The computational process of the proposed fast MoM-VSBR method can be divided into three components: the ray tracing process within the VSBR-MoM region (as shown in Figure 12a), the matrix-solving process in the MoM region (as shown in Figure 12b) and the iterative solution process of the coupled field accelerated by local current theory and FMM (as shown in Figure 12c). During the computation of the coupled field, the use of the FMM theory and the local current coupling solving approach implies that the unknowns are associated with the number of elements in the two regions. Consequently, the computational complexity is about $O\left((N_{VSBR-M\mathrm{o}M})\log (N_{VSBR-M\mathrm{o}M})\right)$.

3.4. Radar Image Simulation of Urban Scenes

This section performs imaging simulation for a multi-scale drone target over the urban ground. As shown in Figure 13a, the dimensions of the maritime drone model are as follows: length 2.4 m, width 1.6 m, and height 1.2 m. The square ground surface has a side length of 8 m, with a root mean square (RMS) height of 0.1 m. The drone is placed 4 m above the ground surface. In this example, the region division of the drone model is the same as in Figure 11 case 1. The dielectric constant of the ground is set to ${\epsilon }_r=3-j0.5$. The incident wave enters from $\theta {=60}^{\circ },\phi {=45}^{\circ }$, with the center frequency of the incident wave set to 10 GHz and a bandwidth of 1.75 GHz. The scan angles are set such that the resolution in both the range and azimuth directions is 0.16 m, with the imaging window being 8 m × 8 m. The proposed work used our fast scattering calculation method combined with local coupling and FMM theory to image, and the imaging results are shown in Figure 13b,c. Figure 13b gives the imaging result simulated by our hybrid method. Figure 13c gives the imaging result simulated by the traditional ray tracing method. From the comparison of the two figures, it can be seen that our method can characterize the radar features of detailed parts such as the drone’s blades, whereas traditional methods cannot achieve this.

Table 5. Calculation time in different steps.

Solvers	Steps	Conditions	Calculation Time
Traditional ray tracing method	Ray Tracing: calculate the amplitude and phase data of RCS over frequency and angle sweeps	Number of Sampling Points: 51 for Frequency, 51 for Azimuth	13.5 h
	2D IFFT: compute the 2D image		3 s
Our method	Ray Tracing + current calculating + couple effect calculating: calculate the amplitude and phase data of RCS over frequency and angle sweeps		47.2 h
	2D IFFT: compute the 2D image		3 s
FEKO-MoM	Calculate the amplitude and phase data of RCS over frequency and angle sweeps		The calculation time for 1 point is more than 20 h

gives the calculation time for different steps. It can be seen that the imaging time for our method and the traditional high-frequency methods is acceptable, with it being 3.4 times less than the traditional method. When using numerical methods for imaging, the computation time for each sampling point is very long, making it impossible to perform imaging simulations.

Finally, we conducted scattering field simulations using the method described in this work and generated radar imaging results for port scenes from multiple radar observation angles in the video remote sensing community. The port is set as concrete structures, with ships moored on the sea surface. The sea surface is modeled using the PM sea spectrum model through linear filtering, with a wind speed of 3 m/s. The central frequency of the incident radar wave is 6 GHz. The structure of the port is shown in Figure 14a, with ships being multi-scale structures, the sea being a rough surface, the oil pipes and trees at the port being non-uniform mediums, and cranes and containers being uniform mediums. The radar imaging simulation observation angle is set to $\theta =60^{\circ }$ and $\phi ={30}^{\circ }~{170}^{\circ }$, the radar imaging bandwidth is 150 MHz, and the scene size for imaging is 200 m × 200 m. The imaging simulation results with an interval of 10 degrees in azimuth angle are provided below, along with the imaging simulation time. It can be seen that the imaging time using our method is relatively short, allowing for the efficient generation of a large amount of simulated data sets. In this scene simulation, the number of unknowns in small-scale regions is low, and the simulation time is shorter. Therefore, multiple imaging results from different angles can be effectively generated.

4. Discussion

In Section 3.1, the building model is created. Here, the computational efficiency of the method presented in this paper is compared with the FDTD method in FEKO, showing an approximately 30-fold improvement in computation speed. In Section 3.2, the tree model is created. In this section, the computational efficiency of the method in this paper is compared with the MoM. As the number of unknowns increases, the computation time of the MoM method rises exponentially, whereas the computation time of the method in this paper increases slowly in a linear fashion. In Section 3.3, a multi-scale model is created and simulated. The multi-scale model includes the VSBR region and the MoM region, and the different region divisions lead to varying solution accuracies. The corresponding solution times also differ. The paper presents a comparison of solution times for the drone model under four different region divisions. It can be seen that the maximum speedup ratio, depending on whether the acceleration method proposed in this paper is applied, is 1:4.2. From the simulation comparison results of the above three types of models, it can be seen that the method presented in this paper demonstrates significant advantages in accelerating the coupled computation process. This method has also been proven to be an effective approach for computing complex targets in urban scenarios. However, the method presented in this paper has the following drawbacks: (1) The distinction between the VSBR and MoM regions has a significant impact on both the computational accuracy and efficiency of the method. How to choose an appropriate region division strategy will be explored in future work. (2) The method in this paper starts from the condition of a plane wave radar, without considering different radar imaging systems. Therefore, future work needs to incorporate actual radar operating modes for simulation verification.

5. Conclusions

A novel hybrid VSBR-MoM method is proposed in this paper to address the challenges of effectively obtaining the scattering characteristics and radar images of complex targets in urban scenes, particularly when dealing with buildings, trees and muti-scale structures. Using the method in this paper can reduce the computational cost of electromagnetic scattering and enable the process of radar imaging simulation using scattered electric fields in the video remote sensing community.

The RCS of complex multi-scale structures is examined using this iterative fast VSBR-MoM method that leverages ray tracing and the concept of iterative reradiation coupled current interaction. Different solvers are employed for the MoM and VSBR regions, allowing for enhanced computational efficiency. The local coupling theory is utilized to decrease the number of unknowns, while the FMM is applied to simplify computational complexity. This method surpasses numerical methods and traditional high-frequency methods concerning time and resource consumption. J.G contributed equally to this work.