Ray-Tracing-Assisted SAR Image Simulation under Range Doppler Imaging Geometry

Abstract

In order to achieve an effective balance between SAR image simulation fidelity and efficiency, we proposed a ray-tracing-assisted SAR image simulation method under range doppler (RD) imaging geometry. This method utilizes the spatial traversal mode of RD imaging geometry to transmit discrete electromagnetic (EM) waves into the SAR radiation area and follows the Nyquist sampling law to set the density of transmitted EM waves to effectively identify the beam radiation area. The ray-tracing algorithm is used to obtain the backscatter amplitude and real-time slant range of the transmitted EM wave, which can effectively record the multiple backscattering among the components of the distributed target so that the backscattering subfields of each component can be correlated. According to the RD condition equation, the backscattering amplitude is assigned to the corresponding range gate, and the three-dimensional (3D) target is mapped into the two-dimensional (2D) SAR slant-range coordinate system, and the SAR target simulated image is directly obtained. Finally, the simulation images of the proposed method are compared qualitatively and quantitatively with those obtained by commercial simulation software, and the effectiveness of the proposed method is verified.

1. Introduction

Currently, it is difficult to obtain typical target SAR images under various imaging conditions by actual measurements, which restricts the ability to accurately detect and recognize targets in SAR images under complex scenes. SAR image simulation can effectively reproduce the complex scattering mechanism between EM waves and typical targets and effectively simulate the geometry and radiometric features of the SAR image. It can also generate typical target simulation images under different imaging parameters, which plays an extremely valuable role in the rapid generation of typical target SAR image data sets and assisting the interpretation of SAR images [1,2].

The feature-based SAR image simulation method uses the RD imaging geometry model, which strictly expresses the relationship between objects and SAR imaging pixels. According to the real-time slant range, the amplitude of each scattering is mapped to each corresponding range gate, and the target echo energy at the same range gate is superimposed. There is no definite phase relationship between inter-pulse and intra-pulse, and no echo generation and focusing imaging steps are involved. Some research institutions and scholars have developed some representative SAR image simulation methods. Stefan Auer et al. of the German Aerospace Center (DLR), based on the ray-tracing algorithm, carry out high-resolution SAR image simulation research, which can map 3D targets directly to the SAR slant-range coordinate system and does not involve the process of echo generation and imaging. Capable of quickly simulating the geometry and radiative properties of SAR images, the team developed typical commercial software for SAR image simulation of 3D target models; at present, it is still widely used in the fields of target data set expansion and decision diagram interpretation [3,4,5]. Eric Pottier’s team at the University of Rennes, France, has developed an open-source software program that focuses on Polarimetric Synthetic Aperture Radar (PolSAR) Synthetic Aperture Radar analysis and teaching. The software includes PolSARprosim, a simulation module designed specifically for simulating forest polarization interference data. This module can simulate the polarization scattering characteristics of vegetation and generate simulated interference data sets. The simulation module has been successfully integrated into PolSARpro, the European Space Agency’s European Space Agency PolSAR data analysis education kit, so that more users can benefit from this advanced simulation technology in actual teaching and research work [6]. Horst Hammer et al., of the Institute of Optoelectronics and Pattern Recognition in Germany, have developed a SAR image simulator (CohRaSS), which uses ray-tracing algorithms to accurately simulate the propagation and scattering of EM waves in complex scenes. The SAR simulation image is generated. It is mainly used to simulate high-resolution SAR images in small scenes and generate sample data for machine learning training and auxiliary interpretation [7,8]. Takao Kobayashi of Tohoku University, Japan carried out SAR imaging simulations of the lunar surface using Gauss rough surfaces to model the topography of the lunar sea region. The Gustav Kirchhoff approximation was used to calculate the complex scattering rate of the superimposed meshes, simulating HF radar detection of the lunar surface and subsurface [9,10]. Tsuchida et al. developed a SAR image simulator for land, ocean, and man-made targets. The simulator divides the surface of a target into small panels; then, the backscattering coefficient is calculated according to different panel attributes, and the simulated SAR image is generated by RD imaging geometry mapping [11]. The Balz Timo team developed SARViz, a SAR image simulation software that integrates computer graphics theory and makes full use of GPU technology to speed up the simulation. The complex simulation environment, including the urban scene, the forest region, and the typical target, is constructed by using computer CAD software, and the surface geometry structure of various targets in the scene is accurately expressed by using a small triangle facet model. According to the principle of SAR imaging and the material properties of the target surface, the scattering characteristics of each surface element are determined, and the speckle is simulated by using the Perlin noise algorithm [12,13]. The Ferdinando Nunziata research team at the National Oceanographic Centre in Southampton, UK has developed an educational tool for SAR wave simulation based on MATLAB; it has been applied to the simulation and analysis of SAR sea surface imaging processes in related studies [14,15]. The University of Illinois and the United States Air Force Research Laboratory jointly developed XPATCH software based on the shooting and bouncing ray (SBR) method for high-frequency EM scattering simulation of targets and SAR image simulation of typical military targets [16,17,18]. The Italian University of Trento Adamo Ferro et al. analyzed the effect of the incident angle of the sensor on the multiple scattering intensity of buildings under different material backgrounds and proposed a prediction model for building structure characteristics based on the perturbation model; building SAR simulation images can be generated [19]. Tao et al. developed a new approach for supporting the automatic interpretation of meter-resolution SAR images in complex urban scenarios [20].

According to the research status of SAR image simulation based on features, most of the existing methods calculate the backscattering field by integrating the scattering intensity of the surface micro-facet. It is convenient for matching with the existing high-frequency EM scattering algorithms based on the plane wave assumption, and good simulation results have been achieved in many application scenarios, but it is also limited by the fixed simulation architecture of commercial EM scattering software. It is generally considered that the current intensity at each part of the micro-facet is equal to the current intensity at the center of the facet, while the current intensity and phase produced by the actual incident EM wave at different positions on the target surface are different. In order to ensure the fidelity of the simulation results, both the cross-section and the facet size of the ray tube should be reduced. At the same time, the size of the target model should be controlled at the wavelength level, which will inevitably result in the over-segmentation of the distributed complex targets and the generation of a large number of surface elements. Even with the support of a highly configured hardware acceleration platform, the computational efficiency of this method is still not ideal, which makes it difficult to carry out the SAR image simulation research of complex distributed targets effectively. It is disadvantageous to the practical application performance of SAR image simulation technology. Therefore, the method of calculating the target backscattering field by using the surface integral of the scattering intensity of micro-facets has great limitations and poor expansibility. It is difficult to improve the traditional SAR image simulation architecture and balance the relationship between simulation fidelity and efficiency.

To address the aforementioned limitations, we proposed a ray-tracing-assisted SAR image simulation method under RD imaging geometry. Firstly, the method adopts the spatial traversal mode of SAR beam to transmit the EM waves while maintaining the dual-scale facet size of the distributed target model. Then, combining the SAR mode ray-tracing architecture and the physical illumination model, the target facets in the radiation area are converted into the corresponding discrete EM waves. The integration of the scattering field of the micro-facet is equivalent to the vector superposition of scattering field of discrete EM waves, refining the local incident angle of the target facets. In the process of calculating the local scattering intensity of the target facets, various scattering mechanisms and the material properties on the surface of each target component are fully considered. The amplitude and real-time slant range of multiple scattering of each discrete EM wave are recorded dynamically by using SAR mode ray tracing architecture to ensure the correct matching between the two, the multiple scattering mechanism can be taken into account when recording the propagation path, and the backscattering sub-fields of the components of the distributed target can be correlated.

The main contributions of this paper are as follows:

• (1). A novel multiple-backscattering simulation method based on ray tracing is designed. The method utilizes the spatial traversal mode of RD imaging geometry to transmit discrete EM waves into the SAR radiation area of the distributed target and combines with the designed SAR mode ray-tracing architecture to simulate the multiple backscattering of EM waves.

• (2). A novel SAR image simulation method under RD imaging geometry is proposed. The method integrates the ray-tracing algorithm and the multiple backscattering mechanism into the RD range condition equation. The backscatter amplitude is assigned to the corresponding range gate, and the 3D target is mapped into the 2D SAR slant-range coordinate system to obtain the SAR target simulated image directly.

This paper is organized as follows: Section 2 introduces the multiple backscattering simulation based on ray tracing. Section 3 introduces the workflow of the SAR image simulation under RD imaging geometry, mainly including the parameters assignment of target functional components and the detailed workflow of mapping imaging of multiple backscattering. Section 4 carries out the SAR image simulation test on an airport scene model, and the simulation results are compared and analyzed quantitatively and qualitatively with those obtained by commercial simulation software to verify the effectiveness of the proposed method. Section 5 presents the conclusions.

2. Multiple Backscattering Simulation Based on Ray Tracing

The main step of the SAR image simulation under RD imaging geometry is the multiple backscattering simulation. The ray-tracing algorithm is the key technology for multiple scattering simulations. The traditional ray-tracing architecture is mostly used for 3D scene rendering, and only the scattering amplitude information is recorded. In the process of SAR image simulation, the scattering amplitude and phase information need to be recorded at the same time, and the real-time slant range is the key element to obtain the phase information. In this paper, SAR imaging geometry is introduced into the ray-tracing algorithm, which is only suitable for computer graphics, and the SAR mode ray-tracing architecture is designed. Combined with the physical illumination model, it can provide key technical support for SAR image simulation.

2.1. Design of the SAR Mode Ray Tracing

The ray-tracing algorithm is a general architecture to study the path propagation of an EM wave. It can be used to simulate the global rendering effect of a multi-light source combined with a local illumination model. In the field of computer graphics, a large number of EM waves are transmitted from the camera position of the sensor to the pixel grids of the virtual screen. Based on Fresnel’s law of specular reflection, the specular reflection direction of an EM wave interacting with the target surface is calculated and used as a new incident direction for forward tracking; thus, the whole path of an EM wave is obtained recursively. In practical application, considering that the natural scene or man-made target surface is a complex with certain gloss and roughness, the reflected energy field has a certain width of lobes. The specular reflection direction calculated from the normal vector of the macro-facet at the intersection of an EM wave and the target surface is used as the energy principal axis direction of forward tracing. The gray value of a pixel on the generated image is the superposition of the energy of all incident rays in the case of multiple scattering. It is necessary to ensure that the position of the intersection point when scattering occurs is in view of the position of the light source and the pixel grids; otherwise, it is considered an invalid scattering energy.

As shown in Figure 1, the 3D object can be rendered globally by combining the ray tracing of the global rendering mode with the illumination model, which can be expressed as

(1)$\mathrm{Im}g\left[pix\right]={\displaystyle \sum _{n=1}^N{\displaystyle \sum _{k=1}^{K}Ra{y}_n\left[I_s^k\right]}}$

where $\mathrm{Im}g\left[pix\right]$ represents a pixel unit in the imaging plane; $n$ is a ray index number that is visible to the pixel unit; $k$ indicates the k-th scattering of light; $I_s$ represents the scattered energy and is suitable for the arbitrary illumination model. The ray-tracing architecture of the whole rendering mode only superposes the scattered energy of the incident light, which is the main difference from the SAR mode. In addition, the sensor radiation mode is also different. Commonly used global rendering light source types are divided into point light source, spotlight light source, and directional parallel light source. (1) The point light source is an ideal light source model, which assumes that the light source transmits light rays uniformly from a point in space in all directions. Its illumination range is theoretically infinite in practice; however, there is usually a limited radius of action beyond which the intensity of the light will be sharply reduced to zero. The brightness value of the light transmitted from a point light source is inversely proportional to the square of the distance from a point to the light source. (2) In addition to having a definite position point, the spotlight source has the characteristics of directivity and limiting illumination range. The light from a spotlight source does not spread evenly but is concentrated in a cone-shaped area. The inner angle of the cone identifies the central region with the highest light concentration and brightness, while the outer angle marks the boundary where the light begins to weaken noticeably. Outside the cone, the light is almost invisible or very faint. Therefore, spotlight mode can better simulate the radar beam width of 3 db. (3) Directional parallel light source is a special type of light source in a virtual environment, which is mainly used to simulate the sunlight or other extremely remote and seemingly negligible sizes of light source. This type of light source has only color and orientation attributes and no specific light source location because the directional light source assumes that the light source is infinitely far from the scene so that all the light seems to come from the same point and hit the scene in parallel and the direction of the light remains the same throughout the scene; there is no angular offset or intensity decay due to the distance of propagation. Therefore, no matter how far away the objects in the scene are from the light source, as long as they face or go back to the light source, the received light intensity is constant and only depends on the relationship between the normal of the object surface and the direction of the light angle.

In the field of SAR image simulation, considering the functional relationship between the EM wave propagation range and target backscattering intensity, the limitation of the 3 db radar beam width and the propagation mode of spherical wavefront, the final selection is a single spotlight radiation mode. The SAR mode ray-tracing architecture combined with EM scattering theory or a physical illumination model can be used to simulate the complex process of SAR multiple backscattering. In contrast to light rendering in the field of computer graphics, the SAR sensor, shown in Figure 2, is an active and singular transmitter, producing an EM wave at the same location as the received echo. Therefore, it is not necessary to consider the contribution of the environment and its own radiant energy to the backscattering field; only the self-energy of the active transmission of the EM wave of the SAR sensor and the continuous decay process of its multiple scattering energies are considered. In general, the sampling time intervals along the azimuth direction are within the second level, and the sampling time intervals along the range direction are within the microsecond level [21]; these can be defined as the slow sampling time along the azimuth direction and the fast sampling time along the range direction, respectively. In the process of SAR image simulation for local small scenes, the 2D sampling time difference can be generally met, and the “go-stop” mode can be adopted.

The main function of the SAR mode ray-tracing architecture is to determine the principal axis direction of the forward scattering lobe of the EM wave after each scattering. With the EM wave incident on the target surface, the direction of the principal axis of the EM wave forward scattering lobe can be determined by Fresnel’s law of specular reflection. As shown in Figure 3, the normal vector of the macro-facet, the incidence direction of the EM wave and the specular reflection direction are coplanar. In practical applications, the specular reflection direction calculated by the normal vector of the macro-facet at the intersection point of the EM wave is the main axis direction of the forward scattering lobe.

(2)$\overrightarrow{r}=2\left(\overrightarrow{N}\cdot \overrightarrow{v}\right)\overrightarrow{N}-\overrightarrow{v}$

where $\overrightarrow{r}$ is the principal axis direction of front scattering lobe, $\overrightarrow{N}$ is the macroscopic normal vector of facet, and $\overrightarrow{v}$ is the incident direction of the EM wave. In addition, the information recorded by the SAR mode ray-tracing architecture mainly includes the amplitude of each scattering, the real-time slant range, and the coordinates of the corresponding intersection points and the material parameters of the intersecting facets.

The phase information is obtained by using the ray-tracing algorithm to assist the real-time slant range information conversion when the electromagnetic wave records scattering, mainly introducing the real-time slant range information into the classical RD condition equation. The conversion process will be derived in detail in Section 3.2, where the representation is simplified to show the initial architecture for obtaining the backscatter field of distributed targets. The geometric mapping between the 3D target and the 2D SAR slant-range image is represented by the key elements of the nearest slant range, the farthest slant range, the 2D sampling intervals, and the real-time slant range. Phase information is essentially the main basis for spatial location assignment of backscatter amplitude. The real-time slant range obtained by our proposed method is transformed into the slant-range coordinate matrix by using the classical RD condition equation, which is equivalent phase information. Then we combine the above designed SAR mode ray-tracing architecture with the usual physical illumination model to carry out the SAR image simulation of the 3D target:

(3)$E_T={\displaystyle \sum _{n=1}^N{\displaystyle \sum _{k=1}^{K}E{M}_n\left[{\left|{\delta }_k\right|}_n\cdot e_n^{Pha\left(R_k\right)}\right]}}$

where $E_T$ is the backscattering field of a distributed target. $\left|{\delta }_k\right|$ is the amplitude of the $k$-th scattering. $Pha\left(R_k\right)$ is the real-time equivalent phase of the $k$-th scattering. $Pha\left(·\right)$ is used to simplify the representation by introducing real-time slant rang $R_k$ into the classical RD condition equation, which is transformed into the coordinate matrix of the slant range that is the equivalent phase information; the derivation process will be detailed in Equations (9)–(12). $R_k$ is half of the sum of the propagation paths of the EM waves from the transmission to the reception. The EM waves $E{M}_n$ are transmitted by traversing the radiating area within the radar beam. $N$ is the total number of transmitted EM waves. The transmitted EM wave is tracked by multiple scattering, and the total $K$ times backscatter subfield is obtained. The total backscatter subfield $E_T$ can be obtained by vector superposition.

Section 2.1 designs the SAR mode ray-tracing architecture and presents the method of combining it with the physical illumination model for SAR image simulation, which provides the basic conditions for the simulation of multiple backscattering of the EM wave. A distributed target of an electrically large size or larger can be considered composed of a large number of discrete scattering units, each of which can be represented by a point target or a small facet. The backscattering field of the distributed target can be calculated by vector superposition of the backscattering subfield of its independent scattering unit.

2.2. Multiple Backscattering Simulation of the Transmitted EM Wave

Based on the known coordinates of the SAR sensor, we use the SAR space traversal method to send the discrete EM wave to the corresponding radiation area of the distributed target. It is difficult to track and simulate the multiple backscattering process of the transmitted EM wave, but the traditional physical illumination model still belongs to the local scattering model; the multiple scattering mechanism between the components of distributed target is not effective. We will combine with the designed SAR mode ray-tracing architecture to simulate the multiple backscattering of the EM wave. As shown in Figure 4, taking the transmitted discrete EM wave as an example, the amplitude and slant range of multiple scattering can be recorded in turn, which can provide key technical support for SAR image simulation. Using the map imaging method under RD imaging geometry, the amplitude of each scattering is mapped to the corresponding range gate according to the real-time slant range, which only contains the change in the incoherent scattering energy of the target. There is no definite phase relationship between inter-pulse and intra-pulse, which can be regarded as an incoherent scattering field.

Combined with Figure 2 and Figure 4, it can be seen that when the sensor moves in turn according to the preset sampling interval along the azimuth direction and sends EM waves, it receives the target echo signal at the same position. The scattering amplitude and the instantaneous slant range when a discrete EM wave has its k-th scattering are two key parameters for calculating the multiple backscattering field.

(1) We follow Nyquist’s sampling law to set the sampling interval along the azimuth and range directions, respectively; the EM wave transmitting density is equivalent to the 2D grid density divided by the 2D sampling interval, and the size of the sampling interval is about half of the theoretical resolution. The simulated EM wave is transmitted into the radar beam radiation area by SAR spatial traverse mode in combination with the instantaneous position of the sensor, can be expressed as:

(4)$\left\{\begin{array}{l}\begin{array}{c}\begin{array}{l}\begin{array}{cc}\left\{\begin{array}{l}{R}_J=R_{\min }+j\ast {\rho }_r\\ {\theta }_j=\arccos (H/R_J)\\ \overrightarrow{S_{\mathrm{dir}}}=\left[0,\sin {\theta }_j,-\cos {\theta }_j\right]\end{array}\right.& j\le N={\displaystyle \frac{R_{\max }-R_{\min }}{{\rho }_r}}\end{array}\\ \left\{\begin{array}{l}{S}_0=\left[x_i,-H\cdot \tan \theta ,H\right]\\ x_i=i*{\rho }_a+x_{\min }\end{array}\right.\begin{array}{c}i\le M={\displaystyle \frac{x_{\max }-x_{\min }}{{\rho }_a}}\end{array}\end{array}\end{array}\\ EM=\left[S_0,\overrightarrow{S_{\mathrm{dir}}}\right]\end{array}\right.$

where $M$ is the number of azimuth sampling points, and $N$ is the number of slant-range gates. $x_{\max }$ and $x_{\min }$ are the azimuth coordinates of the scene; $R_{\min }$ and $R_{\max }$ are the near range and far range, respectively; ${\rho }_r$ is the slant-range sampling space; $R_J$ is the slant range within the zero-doppler plane; ${\theta }_j$ is the incidence angle; and $\overrightarrow{S_{\mathrm{dir}}}$ is the incidence direction. $S$ is the sensor position; $x_i$ is the sensor azimuthal coordinate; ${\rho }_a$ is the azimuthal sampling space; $\theta$ is the incidence angle; $EM$ is the simulated incident EM wave; $P$ is the coordinates of any point on the target surface; $S_0$ is the platform coordinates, which can be determined by the discrete time axis along the azimuth direction, the platform velocity, and the scene range.

(2) According to Fresnel’s law of reflection, the forward propagation direction of the EM wave is calculated at the point of intersection of the EM wave and the target surface. In the plane in which the reflected EM wave, the incident EM wave, and the normal vector of the intersection point lie, the reflection angle is equal to the incident angle and can be expressed as:

(5)$\left\{\begin{array}{l}\overrightarrow{V_s^k}=\overrightarrow{V_{in}^k}-2*(\overrightarrow{V_{in}^k}\cdot {\overrightarrow{N}}_k)*{\overrightarrow{N}}_k\\ P_k=S_{k-1}+R_k^{\prime }\ast \overrightarrow{V_{in}^k}\end{array}\right.$

where $\overrightarrow{V_s^k}$ denotes the specular reflection direction of the k-th scattering of the EM wave, and the EM wave continues to propagate forward along the specular reflection direction of the energy principal axis; $\overrightarrow{V_{in}^k}$ is the incident direction of the k-th scattering of the EM wave; $R_k$ is the instantaneous slant range of the k-th scattering; $P_k$ is the intersection coordinate of the k-th scattering. The process of simulating multiple scattering of the EM wave is essentially the process of taking the specular reflection direction of the principal axis of EM wave as the incident direction of EM wave.

(3) The intersection of the EM wave and the target model panel is taken as the starting point of EM wave forward propagation, and the specular reflections are taken as the new incident direction of the EM wave. The recursive depth of multiple scattering is determined by the actual propagation path and energy decay condition of the EM wave. In the process of multiple scattering, the EM wave follows the conservation of energy, which is manifested in the fact that the outgoing energy is not higher than the incident energy, and the propagating energy will decrease continuously with the occurrence of multiple scattering. Usually, the surface material of the target has a certain glossiness and roughness at the same time, and it is easy to occur in the directional scattering dominated by the specular reflection mechanism. The macro-specular reflection direction is the main energy axis, and some diffuse energy including subsurface scattering continues to propagate along the main energy axis. In order to calculate the forward propagation energy of the EM wave after each scattering, the macroscopic energy decay coefficient is set for each component facet of the distributed target to quantify the energy loss value of the EM wave after each scattering; the energy remaining value is used as the main energy of the forward propagation. The EM wave intersects with different parts of the target surface, and the material properties of each part are different, which has a different influence on the propagation energy decay. In the experiment, we set the energy threshold to prevent the EM wave from entering the loop, and we can judge and track the complex multiple scattering path of each discrete EM wave.

(4) Due to the complexity of the geometry structure of distributed targets, the multiple scattering intersection point and the instantaneous sensor position cannot always be intervisible. It is necessary to add the intervisibility condition, that is, a ray is transmitted from the intersection point to the instantaneous sensor position to see whether there is an obstruction. If there is, the k-th scattering intersection point and the sensor are not intervisible, and the backscattering energy and phase of the intersection point cannot be recorded. The specific process can be expressed in the following Equation (6):

(6)$\left\{\begin{array}{l}\overrightarrow{d}=normal(P-S)\\ P^{\prime }=S+t*\overrightarrow{d}\end{array}\right.$

where $\overrightarrow{d}$ is the vector direction of the instantaneous sensor $S$ transmitting rays towards the current intersection point $P$; $t$ is the instantaneous distance from $S$ to $P$; $P^{\prime }$ is the actual intersection point position. If $P^{\prime }$ and $P$ are located at the same position, record the current backscatter amplitude and real-time slant range; Otherwise, this situation will be considered an occlusion state, and the contribution of the current backscatter information will be directly ignored.

(5) The coordinates of the intersection point of the k-th scattering of the EM wave and the target surface are obtained by steps (1) to (4), and the instantaneous slant range of the k-th scattering is calculated by the following equation.

(7)$R_k=\left(\sqrt{{\left(S^{\prime }-P_0^{\prime }\right)}^2}+\sqrt{{\left(P_0-P_1^{\prime }\right)}^2}+\dots +\sqrt{{\left(P_{k-1}^{\prime }-P_k^{\prime }\right)}^2}+\sqrt{{\left(P_k^{\prime }-S^{\prime }\right)}^2}\right)/2$

where $R_k$ denotes the slant of the k-th scattering equal to half of the sum of the distances from the sensor’s instantaneous position $S^{\prime }$ to the point of intersection ${P^{\prime }}_k$ of the k-th scattering and then to $S^{\prime }$. $P^{\prime }$ denotes the coordinates of the intersection of each scattering. $S^{\prime }$ is the instantaneous position of the sensor.

In the process of multiple scattering, the EM wave always intersects the surface facet with different material properties. The normal vector of the facets is divided into the macro-normal vector and the micro-normal vector. The macro-normal vector of the facet is an important variable for recursively determining the forward tracking direction of the EM wave, and it is also a key element for calculating the backscattering field of the facet.

The normal vector of the micro-facet is the reasonable correction of backscattered field calculation value from the perspective of the micro-facet. The index for quantifying the glossiness and roughness of the target surface material is set according to the regularity degree of the normal vector distribution of the micro-facet. The specular reflectivity coefficient, diffuse reflectivity coefficient, specular index, energy decay coefficient, and permittivity are set according to the type of material parameters in the physical illumination model. The normal vector and index number of the intersecting panel are recorded when calculating the coordinates of the intersection point of the k-th scattering of the discrete EM wave and the target surface. The index number is used to determine the permittivity, diffuse reflection coefficient, specular reflection coefficient, specular index, energy decay coefficient of materials, and facet normal vectors.

According to the physical illumination model of SAR mode, the backscatter energy of the target mainly includes diffuse reflection energy and specular reflection energy. The ray-tracing algorithm is used to track the propagation path of each discrete EM wave. When the geometry structure of the selected distributed target is relatively regular, the termination condition of multiple scattering is that there is no collision point to end recursion. Assuming that the recursive depth of path tracing for discrete EM waves is k, the k-th backscattering field can be expressed as

(8)$\left\{\begin{array}{l}\left\{\begin{array}{l}{I}_{s(k)}=\left[K_d{\displaystyle \frac{1}{\pi }}+K_f\ast \max {(0,\overrightarrow{V_s^k}\cdot {\overrightarrow{r}}_k)}^{K_S}\right]\ast \rho \left(\epsilon \right)\ast I_{in(k)}\ast \max (0,\overrightarrow{V_{in}^k}\cdot \overrightarrow{N_k})\\ I_{\mathrm{in}(k)}=I_{in(1)}\ast (1-K_{los1})*(1-K_{los2})\dots \dots (1-K_{losk})(I_{in(k)}<I_{end})\\ {\rho }_{hh}\left(\epsilon \right)={\displaystyle \frac{\cos (\overrightarrow{V_{in}^k}\cdot \overrightarrow{N_k})-\sqrt{\epsilon -{\sin }^2(\overrightarrow{V_{in}^k}\cdot \overrightarrow{N_k})}}{\cos (\overrightarrow{V_{in}^k}\cdot \overrightarrow{N_k})+\sqrt{\epsilon -{\sin }^2(\overrightarrow{V_{in}^k}\cdot \overrightarrow{N_k})}}}\\ {\rho }_{vv}\left(\epsilon \right)={\displaystyle \frac{\left(\epsilon -1\right)\left({\sin }^2(\overrightarrow{V_{in}^k}\cdot \overrightarrow{N_k})-\epsilon \left(1+{\sin }^2(\overrightarrow{V_{in}^k}\cdot \overrightarrow{N_k})\right)\right)}{{\left(\epsilon \cos (\overrightarrow{V_{in}^k}\cdot \overrightarrow{N_k})+\sqrt{\epsilon -{\sin }^2(\overrightarrow{V_{in}^k}\cdot \overrightarrow{N_k})}\right)}^2}}\end{array}\right.\\ \left|\delta (\tau ,\eta ,k)\right|={\displaystyle \frac{I_{s(k)}}{I_{\mathrm{in}1}}}\end{array}\right.$

where $k$ indicates the $k$-th backward sub-scattering. $I_s$ is the backscattering energy, $I_{\mathrm{in}}$ is the incident energy, $I_{in(1)}$ is the initial incident energy of the EM wave, and $I_{in(k)}$ is the residual energy value at the $k$-th scattering. $\delta$ is the backscatter coefficient. $\overrightarrow{N}$ is the normal vector of the surface facet; $K_d{\displaystyle \frac{1}{\pi }}$ indicates the diffuse scattering mechanism to the surrounding area. $K_f$ is the specular reflection coefficient, $K_{los}$ is the energy decay coefficient, and $K_s$ is the specular index. ${\rho }_{vv}$ is the horizontal polarimetric scatter coefficient, and ${\rho }_{vv}$ is the vertical polarimetric scatter coefficient. $\epsilon$ is the relative permittivity. $\overrightarrow{V_{in}}$ is the direction of EM wave incidence, and $\overrightarrow{r}$ is the direction of EM wave specular reflection. $\overrightarrow{V_s}$ is the direction of energy reception, and $I_{end}$ is the energy decay threshold, which can be used to terminate the tracking of the complex multiple scattering path of each discrete EM wave.

The value of $k$ is defined adaptively by the recursive energy decay under the constraint of the law of conservation of energy. In fact, the iterative depth value $k$ of multiple scattering of electromagnetic waves is determined by the complexity of the geometric structure of distributed targets to a greater extent. For example, the geometric structure of the local scene is relatively flat, and the path tracking will stop after only one scattering. If the geometric structure of the local scene is complex, the incident electromagnetic wave will scatter many times. Every time the incident electromagnetic wave scatters on the target surface, its energy decays and then stops tracking when the energy decays below the energy threshold $I_{end}$.

It is concluded that the simulation of multiple backscattering based on SAR mode ray tracing is an important step in SAR image simulation, which directly affects the quality of SAR images. Firstly, the distributed targets within the SAR beam radiation area are transformed into discrete EM waves by using the SAR beam scanning method. Combined with SAR-mode ray-tracing architecture and the physical illumination model, the amplitude and real-time slant range of each EM wave when multiple scattering occurs are recorded, which are the necessary parameters for the calculation of the target multiple backscattering field. Then, the amplitude and slant range of multiple scattering are brought into the range condition equation, and the amplitude is directly mapped to the 2D SAR slant-range image to generate a SAR simulation image.

3. SAR Image Simulation under RD Imaging Geometry

SAR image simulation under RD imaging geometry has certain similarities with ray rendering in the field of computer graphics, both using the ray-tracing algorithm to record the complete propagation path of the EM wave. High-frequency radar has a seemingly light-like scattering mechanism for electrically large-sized targets; the difference is that the former uses RD imaging geometry, which will acquire the amplitude and slant range of the EM wave each time the wave occurs. Scattering is converted into a scattering field, which is vectorially superimposed and mapped into a 2D SAR slant-range image. However, the ray rendering only adds the amplitude of multiple scattering to the central perspective image.

The main step of SAR image simulation under RD imaging geometry is the multiple backscattering simulation, as shown in Figure 4. Firstly, SAR beam scanning is used to transmit the EM wave towards the radiation area of the beam. The starting position of the transmitted EM wave is determined by the real-time position of the sensor. The slant sampling interval of each range gate is taken into the range condition equation to calculate the corresponding slant range, and the real-time incident direction of the EM wave is determined according to the geometry relationship between the real-time slant range and the incident angle within the zero-doppler plane. The ray-tracing algorithm is used to track the multiple scattering path of each EM wave transmitted, and the amplitude of each backscattering is calculated by the physical illumination model of SAR mode designed in literature [21,22,23]. At the same time, the intersection coordinates of each scattering are recorded to calculate the real-time slant range.

It is worth noting that the real-time range is the key variable in mapping the 3D scene to the 2D SAR image. After multiple scattering path tracking, the multiple scattering real-time range of each transmitted EM wave is brought into the range condition equation to calculate the corresponding range gate, and the backscattering amplitude at the same range gate is superimposed. When the radar beam traverses the scene, the corresponding range mapping image can be obtained directly, which has SAR imaging geometry and radiation characteristics.

3.1. Parameters Assignment of Target Functional Components

Considering the correspondence between the function of each component of the distributed target and its surface material, we disassemble each component according to its function and assign the parameters of the surface material to it before assembling it. As shown in Figure 5, taking an airport as an example, according to the types of functional components of the airport, it can be roughly divided into five parts: runway, aircrafts, vehicles, buildings, and bare ground. In the actual experiment, we take the Google satellite image as the base image and build the corresponding airport scene model. It is convenient to arrange the targets in the scene, and the simulation scene can be flexibly changed according to the actual application demand. As shown in Figure 6, the target components in the airport scene are disassembled one by one according to their functions, and the fineness of disassembly is flexibly determined by the amount of target modeling information available.

Several OBJ function parts corresponding to the simulation scene are read in turn, and text files with different material parameters are given. As shown in Figure 7, taking the function part OBJ-1 as an example, its corresponding material parameter file is Txt-1, building the structure to store facets information, including geometry coordinate and material parameters. In the actual simulation process, the geometry data (identification number of vertex coordinates: v and identification number of each facet: f) are used to complete the construction of the scene facets and read the corresponding target material parameters Txt-1. We use the SAR mode physical illumination model designed in literature [21,22], material TXT records the relative permittivity, diffuse coefficient, specular coefficient, specular index, and energy decay coefficient.

Finally, the OBJ model of the disassembled parts is sequentially read and input to complete the splicing, and the facets of different parts have different attribute parameters, which are saved in the form of a list of facets. In the process of program design, the process of reading and assigning parameters to the above target functional components is set up as an independent module; an independent interface is reserved for the continuous iteration and improvement of the backscattering model. In the process of the continuous completion of the backscattering model, the more comprehensive the material parameters considered are, the more likely the material Txt file is always matched with the content of the material parameters considered in the backscattering model.

3.2. Analysis of Range Gate Equivalent Position

It is worth mentioning that the proposed method belongs to the category of feature-based SAR image simulation, which is similar to RaySAR simulation software developed by DLR and POLSAR Pro simulation software [3,4,5,6]. Feature-based SAR image simulation is mostly based on the assumption of a high-frequency far-field plane wave, which only superimposes the target echo energy at the same range gate in the center of the beam; there is no definite phase relationship between inter-pulse and intra-pulse, and it does not involve echo generation and focusing imaging.

The SAR image simulation method under RD imaging geometry essentially regards each pixel position in an SAR target 2D slant-range image as the ideal position of the raw radar signal after 2D focusing. This position is shown in the bright yellow square in Figure 8, corresponding to the Doppler center frequency and the range center frequency of the echo signal. The transmitted EM wave of each scattering unit on the target surface is sampled once by the non-broadening linear frequency modulation (LFM) signal, and the sampling amplitude is assigned to the corresponding range gate only along the range direction according to the slant range. The energy of the sampling points located at the same range gate is superimposed, and there is no definite phase relationship between the azimuth lines and the range lines.

The method of SAR image simulation under RD imaging geometry does not involve the steps of echo generation and focused imaging. In order to better analyze the equivalent relationship of range gates under RD imaging geometry mapping, the original echo baseband signal expression is used as the cut-in point to deeply analyze the SAR-simulated image characteristics:

(9)$s_0(\tau ,\eta )=A_0{w}_r\left(\tau \right)w_a\left(\eta \right)\exp \left\{-\mathrm{j}4\pi f_0{R}_0/c\right\}\exp \left\{-\mathrm{j}2\pi f_{\eta c}\eta \right\}$

where $s_0(\tau ,\eta )$ can be represented as a 2D SAR slant-range image directly mapped by distributed targets and $(\tau ,\eta )$ as a 2D sampling time corresponding to SAR 2D slant-range image pixel coordinates. $A_0$ is the amplitude at the center frequency of the echo signal corresponding to a target point on the target surface, which is always in the 2D envelope $w_r\left(\tau \right)w_a\left(\eta \right)=1$ constructed by $w_r$ and $w_a$.

The SAR simulation image under RD imaging geometry does not consider the change of the position relationship between the sensor and the target during the synthetic aperture time. The instantaneous slant range is always the nearest slant range of the target at each sampling point. There is no range migration at each target point during the sampling process, and the azimuth is always located at the Doppler center frequency $f_{\eta c}$. The amplitude of the 2D slant-range image generated by direct mapping is equivalent to that of the raw echo signal after 2D pulse compression. In the process of SAR image mapping, the target points are not sampled by the 2D broadened LFM signal multiple times, and there is no interference of 3 db beam sidelobe energy. This method can intuitively obtain clear geometry texture and radiation information of a distributed target SAR simulation image, which can be used to assist the SAR image interpretation.

In the side-looking strip mode, the sampling points are usually located at the position of zero-doppler frequency and range zero frequency. In the SAR image simulation method under RD imaging geometry, the sampling points are directly located at the 2D zero-frequency sampling position. The whole SAR 2D slant-range image is essentially a 2D grid divided by the zero-doppler line and the range zero-frequency line. The key quantity that maps the 3D target sampling amplitude directly to the SAR 2D slant-range coordinate system is the selected reference slant range. The starting point of the equivalent range sampling time axis is usually selected at the close end of the slant range, and the sampling amplitude of the EM wave is divided into the range zero-frequency position according to the reference slant range.

3.3. Mapping Imaging of Multiple Backscattering

The number of azimuth sampling points $M$ and the number of range gate $N$ are determined according to the range of scene boundary. The reference Equations (7) and (8) can be used to calculate the amplitude and real-time slant range of each EM wave when the k-th scattering occurs:

(10)$\left\{\begin{array}{l}{R}_k=\left(R_1+R_2+\dots +R_k+R^{″}\right)/2\\ {\sigma }_k\sim R_k\end{array}\right.$

Combining the RD condition equation with the slant-range history of EM wave scattering, the following equation is obtained:

(11)$\left\{\begin{array}{l}{R}^2={\left(X-X_S\right)}^2+{\left(Y-Y_S\right)}^2+{\left(Z-Z_S\right)}^2={\left(R_{\mathrm{near}}+{\rho }_r\cdot y\right)}^2\\ V_X\left(X-X_S\right)+V_Y\left(Y-Y_S\right)+V_Z\left(Z-Z_S\right)=-{\displaystyle \frac{\lambda R}{2}}{f}_{dc}\end{array}\right.$

where ${\rho }_r$ is the slant-sampling interval of the SAR system; $y$ is the position of the range gate; $\left(\begin{array}{ccc}{X}_S& Y_S& Z_S\end{array}\right)$ corresponds to the coordinates of the radiation area of the lattice; $R$ is the slant range when the first scattering occurs; $R_{\mathrm{near}}$ is the nearest range delay and the slant-range reference for SAR image simulation. The real-time slant range of multiple scattering is taken into the RD condition equation to assign the target backscattering amplitude to the corresponding range gate. The position of the range gate for each scattering can be calculated:

(12)$\left\{\begin{array}{l}{y}_k=⌈\left(R_k-R_{\mathrm{near}}\right)/{\rho }_r⌉\\ y_k~{\sigma }_k\\ R_{\mathrm{near}}=\sqrt{{\left(H\cdot \tan \theta +y_{\min }\right)}^2+H^2}\end{array}\right.$

where $y_k$ is the coordinate value of the slant range of the SAR simulation image according to the position of the range gate. The process of the scanning scene with the zero-doppler plane is regarded as the process of sending the pulse train by pulse repetition time. Then the zero-doppler scanning plane set is traversed, and the real-time slant range of the transmitted EM wave when multiple scattering occurs is brought into Equation (11); the corresponding backscattering amplitude is divided into the corresponding range gate. The 2D SAR-simulated image in the slant-range coordinate system can be obtained by vector superposition of the backscatter amplitude in the same range gate:

(13)$\left[\sigma \right]={\displaystyle \sum _{m=1}^M{\displaystyle \sum _{k=1}^K{\sigma }_k^m}}$

In the experiment, the amplitude and the real-time slant range of the EM wave multiple scattering are calculated by using the multiple scattering simulation method of the EM wave designed in Section 2, and they are stored into two matrices. At the same time, the amplitude matrix and the pixel coordinate matrix corresponding to the slant range are formed by introducing the real-time recorded total $k$ scattering amplitude and the real-time slant range into the RD condition equation. The 2D SAR slant-range image is generated by direct mapping, as shown in Figure 9. Taking the 4 × 4 grid recorded in a local scene as an example, the values of each cell of the slant-range matrix are integer and disordered; however, there is a certain matching relationship with the cells and positions of the amplitude matrix. Each real-time scanning plane corresponds to an azimuth line, traversing the same azimuth line of two matrices and superimposing the corresponding amplitude values of the same range gate coordinates, as shown in Figure 9. Along the first row of azimuth lines, the amplitude values at the third range gate are ${\sigma }_2$ and ${\sigma }_3$, respectively.

Then fill the result of ${\sigma }_2+{\sigma }_3$ into the third range gate in the first row of azimuth line. If there is no corresponding range gate coordinate in the slant-range coordinate matrix, for example, there is no amplitude value in the second range gate in the first row of azimuth line; the amplitude value is 0 in the final 2D simulation image matrix, showing black. After traversing the different azimuth lines, according to the above operation steps, different amplitude results are filled into the SAR slant coordinates prefabricated according to the 2D boundary of the scene to generate the SAR target simulated image.

The SAR image simulation method under RD imaging geometry adopts the EM wave transceiver mode of beam center slice to obtain the backscattering amplitude and real-time slant range corresponding to the transmitted EM wave. Taking the nearest slant range as the reference range, the multiple backscattering of the EM wave is divided into the corresponding range gate, and the backscattering amplitude of the 3D target is directly mapped into the 2D SAR slant-range coordinate system.

4. Discussion of Simulation Results

4.1. Test Parameters

The SAR system parameters used in the experiment are given in Table 1, and the computer configuration parameters are given in Table 2. SAR image simulation under RD imaging geometry does not involve the process of echo signal generation and imaging. The sampling interval along azimuth and range direction, incident angle of scene center, platform height, and radar velocity are used. The method we propose strictly expresses the relationship between the object and the image in SAR imaging. According to the real-time slant range, the amplitude of each scattering is mapped to each corresponding range gate, and the target echo energy at the same range gate is superimposed.

We choose the typical airport scene for SAR image simulation. Airport is a typical military target with high value, and it is also an important object for SAR image simulation technology. In order to obtain the property parameters of the material, it is necessary to sample the material of the target surface in the scene and carry out a comprehensive physical test. Because some of the objects selected in the experiment are non-cooperative, the parameters of the surface materials of each part of the object model are reasonably set by referring to the physical properties of the principal components of the material and by objective experience, there is an inevitable deviation from the real material parameters that can be continuously adjusted to approach the actual value, the more realistic the simulation results.

During the experiment, we set the permittivity of the vacuum as unit 1, and the relative permittivity of metal-like materials is larger than that of vacuum. The reflection coefficient of a material surface is related to the statistical probability of the normal vector distribution of its micro-facets: when the diffuse reflection coefficient is 1, the distribution of the normal vectors on the micro-facets is completely irregular. When the specular reflection coefficient is 1, the distribution of the normal vectors on the micro-facets is uniform. The reflection index is used to restrict the projection range of the reflected energy in the receiving direction. The natural target surface is usually a composite with a certain degree of roughness and gloss, and the diffuse reflectivity and specular reflectivity are between 0 and 1. The energy decay coefficient is set to quantify the absorbing ability of different materials, and then the recursion number of multiple scattering is controlled.

In the process of the actual simulation experiment, we first use the Google image of the airport as the base map to collect the relevant data of the surface buildings and residents according to the actual scene, reasonable arrangement of aircraft, vehicles, runways, buildings, bare ground, and other compound scene model. In the process of parameterizing the target part of the scene, considering the relatively smooth surface of the runway and the hardened road, a large specular reflection coefficient is set. While the bare surface of the low vegetation is relatively rough, a larger diffusion reflection coefficient is set. The building surface is made of metal steel structure, the airplane has some metal shells, and a large specular reflection coefficient is set.

4.2. Analysis of Simulation Results

We carry out SAR image simulation experiments according to the above-set SAR system parameters and the airport scene model to be simulated, as shown in Figure 10a,b. The typical airport scene model and its corresponding SAR simulation images are given respectively. To facilitate the qualitative and quantitative evaluation of the proposed method, the SAR simulation images provided by SE-Workbench-RF and RaySAR commercial software are used as reference images, as shown in Figure 10c and Figure 10d, respectively. It is worth noting that SE-Workbench-RF is a professional tool kit for RCS calculation and SAR image generation and analysis by the French OKTAL-SE company. The software main module SE-RAY-SAR can be used for SAR image simulation and has been applied in many important projects in Europe. The RaySAR software developed by the German European Space Agency has also been widely used in the field of SAR image simulation. Therefore, the selected reference images have certain professional authority, which can assist in further qualitative and quantitative analysis and verification of the method proposed in this paper.

As can be seen from Figure 10 above, the simulated images using our proposed method are consistent with those of commercial software, and there is a strong correspondence between the SAR imaging geometry and radiation characteristics of the airport scene. The simulation results of commercial software shown above have been widely used in practical projects, and the fidelity of the simulated image has been widely approved in the field of SAR image simulation. To further highlight the fidelity of the simulation results, we select two local areas of the airport scene for detailed simulation evaluation, as shown in Figure 11, where the runway is relatively flat. The sensor receives a weak backscattering energy, which is dark in SAR images. Metal targets of aircrafts and vehicles have a strong backscattering intensity. Some buildings have a certain height and have the characteristics of overlapping and masking. The above qualitative analysis shows that the SAR imaging characteristics of the targets in the airport scene are consistent with the real SAR imaging mechanism, and the validity of the method is preliminarily proved.

To further validate the effectiveness of our method, we also collected a 3D truck model (8.76 × 5.20 × 2.75 m) in the parking of zoom-2 scene part with 35,928 facets, as shown in Figure 12. The real SAR image of the truck is detailed in Figure 13b with a theoretical resolution of 0.3 m. The model size of the truck was magnified 10 times to match the resolution of the measured image. The simulated image using the proposed method is shown in Figure 13a, and the simulation time consumption is 0.09 h. The simulated images of SE-RAY-SAR and RaySAR are shown in Figure 13c and Figure 13d, respectively. The simulation time consumption levels are 0.81 h and 0.96 h, respectively. We compare the simulated images with the real SAR image in detail, and it can be seen from Figure 13 that the proposed method has better fidelity and the highest simulation efficiency.

In order to verify the validity of the proposed method quantitatively, we use four commonly used indexes to evaluate the similarity of two images besides visual qualitative analysis of the simulated images, the evaluation indexes include cosine similarity, structure similarity, histogram similarity, and mean hash similarity. According to Equations (13)–(16), the simulated SAR image is evaluated comprehensively and quantitatively from the aspects of geometry and radiation characteristics. In fact, a SAR-simulated image is more concerned about its relative reference image structure similarity, focusing on the structure similarity between simulated image and the reference image.

(14)${\beta }_1={\displaystyle \frac{{\displaystyle {\sum }_{x,y}\left(A\left(x,y\right)\times B\left(x,y\right)\right)}}{\sqrt{{\displaystyle {\sum }_{x,y}A{\left(x,y\right)}^2}}\times \sqrt{{\displaystyle {\sum }_{x,y}B{\left(x,y\right)}^2}}}}$

where ${\beta }_1$ is the cosine similarity; $A$ and $B$ are the gray matrix of the simulation image and the reference image, respectively. The gray matrix of the simulation image and the reference image is expressed as two vectors, and the similarity of the two images is expressed by calculating the cosine distance between the vectors.

(15)$O_2={\displaystyle \frac{2\left(\overline{A}\overline{B}+C_1\right)\left(2{\sigma }_{12}+C_2\right)}{\left({\overline{A}}^2+{\overline{B}}^2+C_1\right)\left({\sigma }_1^2+{\sigma }_2^2+C_2\right)}}$

where $O_2$ is the structural similarity; ${\sigma }_1^2$ and ${\sigma }_2^2$ are the variance of gray-scale matrix of $A$ and $B$ images, respectively. This coefficient is used to describe the local contrast between the simulation image and reference image. ${\sigma }_{12}$ is the covariance of the $A$ and $B$ images’ gray-scale matrix, which reflects the structural similarity between the simulation and reference image. $C_1$ and $C_2$ are stability factors that prevent the denominator from dividing by zero when it is close to zero, usually setting a small positive value; the normalized values are set to 1, $k_1$ and $k_2$ are small empirical constants, set to 0.01 and 0.02, respectively.

(16)$S_3={\displaystyle \frac{1}{K}}{\displaystyle {\sum }_{x,y}\left(1-\frac{\left|A\left(x,y\right)-B\left(x,y\right)\right|}{\max \left(A\left(x,y\right),B\left(x,y\right)\right)}\right)}$

where $S_3$ is histogram similarity, $K$ is the total number of pixels, and $\max \left(A\left(x,y\right),B\left(x,y\right)\right)$ is the maximum of corresponding position pixels. This index is based on the histogram to count the number of pixels of different gray values in the image, and then compare the image similarity.

(17)$\left\{\begin{array}{l}ListA=bool\left[resA\left(x,y\right)\ge \overline{A_s}\right]\\ ListB=bool\left[resB\left(x,y\right)\ge \overline{B_s}\right]\\ H_4=1-{\displaystyle \frac{{\displaystyle {\sum }_{x,y}bool\left[ListA\ne ListB\right]}}{N^2}}\end{array}\right.$

where $H_4$ is the mean Hasche similarity, $bool$ is 1 if the criteria are met, $res$ is the uniform scale of the image, $k_s$ is the gray matrix element of the scaled image. $\overline{A_s}$ and $\overline{B_s}$ are the gray mean values of 8 × 8 sub-areas in the scaled image gray matrix, $ListA$ and $ListB$ are the records in the form of a list, $N$ is the scale size of the image, usually set to 32 × 32.

The results from Table 3 show that the simulated image of airport scene has high structural similarity with SE-RAY-SAR’s simulated image and RaySAR’s simulated image. The simulated image of the truck has high structural similarity with real SAR image, and it is consistent with the RaySAR’s simulated image and closer to the real SAR image. It can be further proved that the proposed method has good fidelity. At the same time, Table 4 gives the total simulation time of the three methods respectively, which shows that the proposed method can effectively balance the relationship between SAR image simulation fidelity and efficiency.

5. Conclusions

We propose a ray-tracing-assisted SAR image simulation method under RD imaging geometry, which can achieve an effective balance between SAR image simulation fidelity and efficiency. RD imaging geometry model is used to distinguish the SAR beam radiation area, and the target in the radiation area is transformed into the corresponding discrete EM waves. The integral of the backscattering field of the micro-facet is equivalently transformed into the summation of the scattering field vector of the discrete EM waves so that it is suitable for the dual-scale facets target model. On this basis, we use the ray-tracing algorithm to obtain the backscattering amplitude and real-time slant range of the transmitted EM wave. While recording its propagation path, we can effectively take into account the multiple backscattering mechanism among the components of the distributed target so that the backscattering subfields of each component are correlated. According to the RD condition equation, the backscatter amplitude is assigned to the corresponding range gate, and the 3D target is mapped to the 2D SAR slant-range coordinate system to obtain the SAR-simulated image directly. Finally, the simulation results are compared with those obtained by commercial simulation software, which effectively verifies the effectiveness of the proposed method. Furthermore, we expand the method to SAR image simulation research with more variety and scenario complexity and continue to improve the simulation efficiency under the premise of ensuring the fidelity of SAR image simulation.

Author Contributions

Conceptualization, J.L. and G.Z.; methodology, J.L.; software, C.C. and K.W.; validation, J.L.; formal analysis, G.Z. and J.L.; investigation, J.L. and C.H.; resources, J.L. and G.Z.; data curation, J.L. and C.H.; writing—original draft preparation, W.Z., C.H. and C.C.; writing—review and editing, J.L. and K.D.; visualization, J.L. and K.D.; supervision, J.L.; project administration, J.L.; funding acquisition, J.L. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

Acknowledgments

The authors would like to express their sincere thanks to the members of the editorial team and the anonymous reviewers for their guidance and valuable comments on this article.

Conflicts of Interest

Author Gaohao Zhu and Chen Hou were employed by the company Henan Geo-Vision Information Technology Co., Ltd. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Global light rendering.

Figure 2. The multiple backscattering process of SAR mode ray tracing.

Figure 3. Fresnel specular reflection.

Figure 4. The simulation process of multiple backscattering.

Figure 5. Detail of Airport Scene Model. (a) Airport Scene Model. (b) Zoom-1 part of the scene. (c) Zoom-2 part of the scene.

Figure 6. Component disassembly and parameterization of airports.

Figure 7. The process of reading and assigning parameters to OBJ function parts.

Figure 8. Equivalent position of range gate under RD imaging geometry.

Figure 9. Mapping imaging process in SAR slant-range coordinates.

Figure 10. Comparison of airport scene simulated images. (a) Airport scene. (b) Simulated image of the proposed method. (c) Simulated image of SE-RAY-SAR. (d) Simulated image of RaySAR.

Figure 11. The details of the airport scene simulated image, including the simulated results of each target within zoom-1 and zoom-2 scene parts.

Figure 12. 3D Truck model in the parking of zoom-2 scene part.

Figure 13. Comprehensive comparison of real SAR image and simulated SAR images. (a) Real SAR image. (b) Simulated SAR image using the method we proposed. (c) Simulated image of SE-RAY-SAR. (d) Simulated image of RaySAR.

Table 1

System parameters.

Parameters	Value
Bandwidth	180 MHz
Pulse duration	1.0 µs
Range sample space	0.75 m
Incidence angle	45°
Center frequency	15 GHz
Platform height	2.0 km
Radar velocity	300 m/s
Doppler bandwidth	200 Hz
PRF	250 Hz
Azimuth sample space	1.20 m

Table 2

Configuration parameters.

GPUVersion	Graphics Memory	Compiler Environment	CPU Version	Total Memory	Operation System
NVIDIA GeForce RTX3060	8 G	VS2019	11th Gen Intel(R) Core (TM) i7-11800H	32 G	Windows 10

Table 3

Quantitative evaluation of simulation results.

Reference Images Source	Model	Cosine Similarity	Structural Similarity	Histogram Similarity	Mean Hash Similarity
SE-RAY-SAR	Airport scene	0.782	0.833	0.871	0.805
RaySAR	Airport scene	0.727	0.826	0.735	0.722
Real SAR image	Truck	0.835	0.861	0.895	0.864
RaySAR	Truck	0.762	0.812	0.796	0.857

Table 4

Time consumption.

Method	Model	Number of Facets	SAR Image Size		CPU Time
Proposed method	Airport scene	111,900	Azimuth	1000samples	0.36 h
SE-RAY-SAR			Range	1479samples	3.72 h
RaySAR					5.98 h

References

1. Holtzman, J.C.; Frost, V.S.; Abbott, J.L.; Kaupp, V.H. Radar image simulation. IEEE Trans. Geosci. Electron.; 1978; 16, pp. 296-303. [DOI: https://dx.doi.org/10.1109/TGE.1978.294587]

2. Chen, K.S. Principles of Synthetic Aperture Radar Imaging: A System Simulation Approach; CRC Press: Boca Raton, FL, USA, 2015.

3. Auer, S.; Hinz, S.; Bamler, R. Ray Tracing for Simulating Reflection Phenomena in SAR Images. Proceedings of the Geoscience and Remote Sensing Symposium, 2008 IEEE International-IGARSS; Boston, MA, USA, 6–11 July 2008; [DOI: https://dx.doi.org/10.1109/IGARSS.2008.4780143]

4. Auer, S.; Hinz, S.; Bamler, R. Ray tracing simulation techniques for understanding high-resolution SAR images. IEEE Trans. Geosci. Remote Sens.; 2009; 48, pp. 1445-1456. [DOI: https://dx.doi.org/10.1109/TGRS.2009.2029339]

5. Auer, S.; Bamler, R.; Reinartz, P. RaySAR-3D SAR simulator: Now open source. Proceedings of the IEEE International Geoscience and Remote Sensing Symposium (IGARSS); Beijing, China, 10–15 July 2016; [DOI: https://dx.doi.org/10.1109/IGARSS.2016.7730757]

6. Williams, M.; Pottier, E. Forest Coherent SAR Simulator within PolSARPro: An Edudational Toolbox for PolSAR and Pol In-SAR Data Processing. Proceedings of the Asian Conference on Remote Sensing; Kuala Lumpur, Malaysia, 12 November 2007.

7. Hammer, H.; Schulz, K. SAR Simulation for Urban Scene Analysis. Proceedings of the European Conference on Synthetic Aperture Radar; Nuremberg, Germany, 23–26 April 2012; [DOI: https://dx.doi.org/10.1117/12.897452]

8. Hammer, H.; Kuny, S.; Schulz, K. Amazing SAR imaging effects explained by SAR simulation. Proceedings of the European Conference on Synthetic Aperture Radar; Berlin, Germany, 3–5 June 2014.

9. Kobayashi, T.; Oya, H.; Ono, T. B-scan analysis of subsurface radar sounding of lunar highland region. Earth Planets Space; 2002; 54, pp. 983-991. [DOI: https://dx.doi.org/10.1186/BF03352446]

10. Kobayashi, T.; Oya, H.; Ono, T. A-scope analysis of subsurface radar sounding of lunar mare region. Earth Planets Space; 2002; 54, pp. 973-982. [DOI: https://dx.doi.org/10.1186/BF03352445]

11. Tsuchida, M.; Suwa, K.; Iwamoto, M. A simulator of synthetic aperture radar images; land, ocean surface, and man-made targets. Proceedings of the International Geoscience and Remote Sensing Symposium; Anchorage, AK, USA, 20–24 September 2004; [DOI: https://dx.doi.org/10.1109/IGARSS.2004.1370245]

12. Balz, T. Automated CAD model based geo-referencing for high-resolution SAR data fusion in urban environments. IEEE Proc. Radar Sonar Navig.; 2004; 153, pp. 289-293. [DOI: https://dx.doi.org/10.1049/ip-rsn:20045121]

13. Balz, T.; Haala, N. Improved Real-time SAR Simulation in Urban Areas. Proceedings of the IEEE International Geoscience and Remote Sensing Symposium; Denver, CO, USA, 31 July 2006; [DOI: https://dx.doi.org/10.1109/IGARSS.2006.930]

14. Nunziata, F.; Gambardella, A.; Migliaccio, M. A simulator for SAR sea surface waves imaging. Proceedings of the IEEE International Geoscience & Remote Sensing Symposium; Barcelona, Spain, 23–28 July 2007; [DOI: https://dx.doi.org/10.1109/IGARSS.2007.4422914]

15. Nunziata, F.; Gambardella, A.; Migliaccio, M. An educational SAR sea surface waves simulator. Int. J. Remote Sens.; 2008; 29, pp. 3051-3066. [DOI: https://dx.doi.org/10.1080/01431160701469008]

16. Hazlett, M.; Andersh, D.J.; Lee, S.W.; Ling, H.; Yu, C.L. XPATCH: A high-frequency EM scattering prediction code using shooting and bouncing rays. Proceedings of the Targets & Backgrounds: Characterization & Representation; Orlando, FL, USA, 17–21 April 1995; Volume 2469, pp. 266-275. [DOI: https://dx.doi.org/10.1117/12.210627]

17. Andersh, D.; Moore, J.; Kosanovich, S.; Kapp, D.; Bhalla, R.; Kipp, R.; Courtney, T.; Nolan, A.; German, F.; Cook, J. et al. Xpatch 4: The next generation in high frequency EM modeling and simulation software. Proceedings of the IEEE International Radar Conference; Alexandria, VA, USA, 12 May 2000; [DOI: https://dx.doi.org/10.1109/RADAR.2000.851945]

18. Bhalla, R.; Lin, L.; Andersh, D. A fast algorithm for 3D SAR simulation of target and terrain using Xpatch. Proceedings of the IEEE International Radar Conference; Arlington, VA, USA, 9–12 May 2005; [DOI: https://dx.doi.org/10.1109/RADAR.2005.1435854]

19. Ferro, A.; Brunner, D.; Bruzzone, L.; Lemoine, G. On the Relationship Between Double Bounce and the Orientation of Builing in VHR SAR Images. IEEE Trans. Geosci. Remote Sens.; 2010; 48, pp. 1445-1456. [DOI: https://dx.doi.org/10.1109/LGRS.2010.2097580]

20. Tao, J.; Auer, S.; Palubinskas, G.; Reinartz, P. Automatic SAR Simulation Technique for Object Identification in Complex Urban Scenarios. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2014; 7, pp. 994-1003. [DOI: https://dx.doi.org/10.1109/JSTARS.2013.2275928]

21. Wu, K.; Jin, G.; Xiong, X.; Zhang, H.; Wang, L. SAR Image Simulation Based on Effective View and Ray Tracing. Remote Sens.; 2022; 14, 5754. [DOI: https://dx.doi.org/10.3390/rs14225754]

22. Wu, K.; Jin, G.; Xiong, X.; Zhang, H.; Wang, L. Fast SAR Image Simulation Based on Echo Matrix Cell Algorithm Including Multiple Scattering. Remote Sens.; 2023; 15, 3637. [DOI: https://dx.doi.org/10.3390/rs15143637]

23. Wu, K.; Jin, G.; Xiong, X.; Shi, Q. Bounding Volume Hierarchy-Assisted Fast SAR Image Simulation Based on Spatial Segmentation. Appl. Sci.; 2024; 14, 3340. [DOI: https://dx.doi.org/10.3390/app14083340]