Research on SAR Imaging Simulation Based on Time-Domain Shooting and Bouncing Ray Algorithm

The signal source of synthetic aperture radar (SAR) usually adopts the linear frequency modulation (LFM) signal, which exhibits the characteristics of a wide pulse. Hence, in the case of usage of the LFM signal by time-domain shooting and bouncing ray (TDSBR) to simulate the SAR echo signal, numerous time sampling points are generated, resulting in huge computational efforts; thereby, it is hard to exploit the TDSBR algorithm for simulating the SAR echo. In order to unlock this dilemma, the hybrid approach, the transfer function in conjunction with the range frequency-domain pulse coherence, is developed, in which the transfer function is stated by the radar cross-section. The proposed methodology is capable of enhancing the computational efficiency through avoiding the massive time sampling, so that the suggested TDSBR approach could be more conveniently applied to the SAR echo simulations. Furthermore, because of the efficiency advantage of the TDSBR in the calculation of wideband scattering field, the proposed methodology exhibits higher computational efficiency than the frequency-domain shooting and bouncing ray in the SAR image simulation. Finally, the key equation of the TDSBR for the transient scattering field of dielectric targets based on a closed-form integration formula is analytically derived, which is the basis of SAR imaging simulation for the composite scene of ships above the sea surface.


I. INTRODUCTION
T IME-DOMAIN (TD) method enables us to physically interpret the behaviors of transient waves. Due to the advantages of the TD approach in the computational efficiency of wideband electromagnetic (EM) scattering, examining responses of transient scattering from complex targets has received increased interest in recent years. Up till now, researchers have proposed several transient scattering calculation algorithms classified as low-frequency numerical methods, high-frequency numerical methods, and hybrid methods. Finite-difference timedomain (FDTD) [1], [2], TD integral equation [3], and TD method of moments [4] are some of the low-frequency numerical methods. Time-domain physical optics (TDPO) [5], [6], [7], the TD physical theory of diffraction [8], TD equivalent edge current [9], [10], and time-domain shooting and bouncing ray (TDSBR) [11] are some of the high-frequency numerical methods. In addition, several hybrid algorithms, including FDTD-TDPO [12], [13], have also been established. The high-frequency method, especially the TDSBR, is more suitable for the EM scattering prediction from electrically large targets because of its advantages in computational efficiency. Ling and Bhalla [14] and Jeng et al. [15] introduced the TDSBR method based on the frequency-domain shooting and bouncing ray (FDSBR) solutions. The FDSBR algorithm mainly involves geometric optics (GO) and physical optics (PO) methodologies [16], [17]. The GO is commonly utilized to trace the path and the field intensity of the reflected rays, while the PO is usually implemented to further analyze the scattering field. Accordingly, the TDSBR includes the time-domain version GO (TDGO) and TDPO. Jeng et al. [14] and Ling et al. [15] have derived a closed-form TD ray-tube integration formula for the transient response of perfect electronic conductor (PEC) targets through employing the inverse Fourier transform (IFT) of the corresponding frequency-domain ray-tube integration formula. There also exist some discrepancies between the TDSBR and the FDSBR as follows. The FDSBR is usually employed to evaluate the field scattered from the target at a single frequency, whose incident source is the continuous plane wave. While the time signal sequence of echo with the physical interpretation can be solved by TDSBR, in which the incident plane wave is a pulse instead of a continuous wave. Moreover, the scattering results with multiple frequencies can be obtained by processing the echo signal with fast Fourier transform (FFT), which makes the TDSBR have better efficiency than the FDSBR in predicting a wideband scattering field or radar cross-section (RCS) [18], [19].
Ray tracing used in the GO and TDGO is a time-consuming process for the SBR (both in the TD and FD). The algorithm based on Octree or KD-tree [20] and the parallel algorithm on basis of graphics processing units improve the efficiency of ray tracing [21]. Zhou et al. [11] improved the efficiency of the TDSBR algorithm by introducing the beam tracing technique. The OpenGL and the neighborhood search algorithm [22] are used to accelerate the ray tracing process in this article, which the related research belongs to the previous work of our team, and will not be introduced in detail.
Guan et al. [23] obtained the TDPO closed integral expression for an ideal conductor target using the Radon transform. The developed approach makes the size of the integral This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ region independent of the wavelength, and so large patches could be used to construct the target that enhances the computational efficiency . Based on this efficient approach, the formula for the transient scattering field from medium targets  is derived in this article, suggesting the TDSBR usable for  calculating the transient response from the ship on the sea  surface. At present, the incident signal sources adopted by the TDSBR are almost narrow pulse signals because the use of wide pulse signals, such as the linear frequency modulation (LFM) signals, produces a large number of time sampling points, resulting in a huge increase in computational burden, which restricts the application ability of the TDSBR method in the synthetic aperture radar (SAR) echo signal simulation. In order to solve this problem, we proposed a method combining the transfer function with the range frequency-domain pulse coherence (RFPC) to simulate the SAR echo with the LFM signal as the transmitting signal. In this way, the TDSBR is capable of using the narrow pulse as the signal source, and the characteristics of the LFM signal are reflected in the range frequency domain. Still, the narrow pulse is exploited as an incident source in this methodology, so the amount of calculation is almost not increased, which considerably enhances the computational efficiency. Because the RFPC implements the wideband scattering field to simulate the SAR echo [24], it can be conveniently combined with the transfer function. Hence, the transfer function and the RFPC method are utilized to simulate the SAR echo signal of the composite scene that the ship on the sea surface, and the range Doppler (RD) algorithm [25] was implemented to produce the SAR images. Since the TDSBR exhibits better efficiency than FDSBR in computing wideband scattering field, the proposed method is capable of saving time consumed in SAR image simulation, which is validated in the experimental part.
This article is organized as follows. The section on methods introduces the formulations pertinent to the TDSBR (TDPO and TDGO), which involved a detailed derivation process from the FD to TD. More importantly, the next formulas show the reason that RFPC can be conveniently combined with the transfer function expressed by RCS to simulate SAR echo signal. In Section IV, the efficacy of the TDSBR is validated through numerical examples, and the effectiveness of the proposed method of the transfer function is demonstrated by the simulation of SAR images. Comparing the difference in time consumption for different methods, the efficiency of our method is illustrated. Finally, the completed work is summarized in the conclusions section.

II. TDSBR FOR EM SCATTERING FROM DIELECTRIC TARGETS
The TDSBR algorithm mainly involves the TDGO and the TDPO methods. The TDGO is used to determine which facets of a target are illuminated by the rays. The TDPO is utilized to evaluate the far-field from those illuminated facets, which can be derived based on the FDPO formulas. These two approaches are introduced in the following subsections.

A. TDPO Method
As shown in Fig. 1, the center position vector of a triangular plane S, with the unit normal vectorn, is signified by r 0 , and r denotes a point on the plane S. The position of the observation point is presented by r , the distance between the original point and the observation point is evaluated by r = | r |. In addition, ishows the direction of the incident wave andŝ presents the scattering direction vector. According to the literature [26], the FD scattering field of the triangular facet can be expressed as Considering a dielectric target, there is a magnetic current component in (1). According to the method of modified equivalent current approximation in [27], the equivalent electric current and equivalent magnetic current in (1) can be expressed as follows: where H tot represents the total magnetic field on the target's surface, E tot is the total electric field, E inc TE signifies the horizontal component of the incident electric field, and E inc TM denotes the vertical component of the incident electric field. Similarly, H inc TE and H inc TM represent the horizontal and vertical components of the incident magnetic field, respectively. The parameters R TE and R TM are the reflection coefficients of the target surface for different polarization conditions, whose corresponding expressions are given in [28]. The factorsê TE andê TM stand for the unit vectors with different polarization directions, as shown in Fig. 2. By substituting (2) and (3) into (1), the scattering field is modified to where A is independent of the integral term, and can be state as follows: (5) Because the wave number is given by k = ω/c, where ω is the circular frequency and c is the speed of light, the FD scattering field can be written as The TD scattering field shown as follows can be obtained by applying the IFT to (6): (7) According to the properties of the impact response and convolution operation, the transient response can be rewritten in the following form: where the symbol " * " represents the convolution operation and δ(·) is the impact response. It is crucial for (8) to solve the integral term. In this article, the Radon transform method is exploited to solve this problem. The integral term can be expressed as follows [23]: in which ΔS denotes the size dimension of the triangular facet and t i is stated as where v i is the vertex of the triangular facet. In order to ensure the calculation accuracy of (7), more fine patches are required for the traditional numerical integration method with the increase of frequency, which in turn is timeconsuming. On the contrary, (9) is essentially analytical, and its integral results are independent of frequency, which can calculate the scattering results of any triangles. Therefore, under the premise of satisfying the accuracy of the model, the targets can be constructed by larger patches to improve the computation efficiency.

B. TDGO Method
Similar to the frequency-domain geometric optics method, the TDGO includes ray path tracing and field intensity tracing. There is no difference between the TD and FD methods for ray path tracing. In order to improve the computational efficiency of path tracing, a variety of ray tracing acceleration techniques have been proposed, such as octree and OpenGL. Herein, the algorithm combines both the OpenGL and neighborhood search techniques proposed in [22] to ensure computational efficiency.
The field intensity tracing of the TD method is different from that of the FD approach. The FD method should solve the phase change of the field intensity, while the TD method focuses on the time delay of the field intensity.
As shown in Fig. 3, the intensity of the reflection field after the ray passes through the first facet can be expressed as If the reflection field propagates from facet 1 to facet 2, the time delay of the scattering field increases further. Thereby, the incident field of facet 2 can be outlined as follows: The second-order TD response could be achieved by substituting (12) into (8). If the target exhibits a higher order reflection field, the derivation of the field intensity tracing would be similar to the methodology mentioned above.

A. Transfer Function Expressed by RCS
Assuming that there are two different incident waves whose spectra are specified by E i1 (f ) and E i2 (f ), respectively, and their scattering spectra are denoted by E s1 (f ) and E s2 (f ). Then, the transfer function H(f ) should satisfy Equation (13) shows that the transfer function does not vary in the form of the incident wave at each frequency point, where the RCS could satisfy this property. Further, the wideband RCS could be easily obtained by implementing the FFT to the time signal generated by the TDSBR method given as where F T [·] signifies the Fourier transform. Because the far field RCS of a specific target does not change with the signal form, so that it can be readily solved by any incident signal.

B. SAR Echo Simulation Based on RFPC and Transfer Function
The formula of the RFPC is given in this part. According to this relation, it can be found that the transfer function approach can be easily combined with the RFPC method for SAR echo simulation, which makes TDSBR method can be extended to the SAR image simulation.
The geometry structure of the SAR imaging scene is provided in Fig. 4, where H denotes the height of the radar, t a represents the radar movement time, R(t a ) is the distance between the radar and the point P in the scene, and the SAR is capable of transmitting and receiving the LFM signal at each of the sampling points. The transmit signal could be presented by The echo signal after scattering is where τ = 2R(t a )/c represents the time delay, K r denotes the frequency modulation ratio, f c = c/λ is the center frequency, and σ(t) is the transient response of the scene. According to (16), the baseband echo signal can be expressed as If the number of scattering sources in the scene is N , then the echo signal can be expressed as According to the properties of the impact response and convolution operation, the echo could be written as follows [24]: where s t (t) represents the baseband signal transmitted by the radar and σ i (t) is the transient response of each scattering source. Generally, the pulsewidth of the LFM signal is much larger than the narrow pulse signal. Therefore, if the TDSBR method is directly utilized for evaluating the transient response, many time sampling points are required, resulting in extremely low-computational efficiency. According to the RFPC, the frequency-domain echo could be attained by the Fourier transform of (18) provided as As can be seen from (20), the wideband scattering response (σ i (f ), as the transfer function) can be efficiently calculated by the TDSBR, where the modulated Gaussian pulse is the incident source. The SAR FD echo (s r (f, t a )) is obtainable by multiplying the transmitted LFM signals expressed by s t (f ). The SAR TD echo can be achieved by the IFT, and finally, the scene SAR image can be simulated by the RD imaging algorithm. This method can avoid the massive computational cost resulted by exploiting the LFM signal as the incident source. Fig. 5 displays the SAR image simulation process based on the TDSBR.

IV. NUMERICAL RESULTS AND ANALYSIS
First, in order to demonstrate the correctness of the TDSBR algorithm, the resulting TD echo signal obtained should be converted to the broadband RCS, and the predicted results are then compared with those evaluated by the frequency-domain method at different frequency points. In this article, the modulated Gaussian pulse signal as a narrow pulse signal is adopted. The specific form of the modulated Gaussian pulse is given as where ω = 2πf c , in which f c is the center frequency. τ denotes a constant, which determines the pulse base width of the Gaussian pulse. It is worth noting that all sampling points of time signals in this article are selected according to Nyquist sampling law. The factor t 0 denotes the moment of the pulse peak, usually taken as 0.8τ . Besides, an example is used to illustrate the low-computational efficiency of the TDSBR when the LFM signal, shown as (15), is used as the incident source. In the second part, the method of combining transfer function and the RFPC is used to carry out the SAR image simulation, and the obtained results are compared with the measured image to illustrate the correctness of the proposed method. Compared with the calculation time of the FD method, the effectiveness of the proposed method is demonstrated.

A. Validation of TDSBR
As shown in Fig. 6, the transient response and monostatic RCS of a plane with a side length of 1 m are compared under two discrete precision conditions. The longest edges of the patches are 1.414 and 0.042 m, and the corresponding numbers of patches are 2 and 2508, respectively. The geometric model is shown in Fig. 7. The bandwidth is in the interval of 10-20 GHz, the center frequency is 15 GHz, the type of polarization is VV, and the pulsewidth is 0.4 ns. Moreover, the incident and the azimuth angles are set as 30 • and 0 • . The comparison results indicate that the size of the patches does not affect the calculation accuracy for the plane due to the simple structure of the plane with no ray reflection. However, the patch size should ensure the accuracy of ray tracing for the target with multiple scattering structures. In addition, if the target contains curved surfaces, the patch size should ensure that it can fit the curved surface, otherwise, serious errors will be caused. Therefore, it is difficult to determine the standard of discrete accuracy for different types of targets. In the following examples, the number of patches is selected in a way as to guarantee the calculation accuracy and efficiency.
In Fig. 8, the transient response and monostatic RCS of the PEC trihedral reflector are calculated. The edge length of the trihedral reflector is set as 0.3 m, the number of facets in this model is 7106, the longest edge of the patch is thereby equal to 0.0085 m, the bandwidth is considered in the range of 2-10GHz, the center frequency is 6 GHz, the polarization is VV, and the pulsewidth is set as 0.5 ns. In Fig. 8(a) and (b), both the incident and the azimuth angles are set as 45 • . The comparison results in Fig. 8(b) prove the effectiveness of the TDSBR in calculating the scattering from multiple targets.
The dimensions of the complex ship model made of the PEC material are shown in Fig. 9. The number of facets in this model is 16188, and the longest edge of the patch is 0.93 m. Furthermore, the transient response and monostatic RCS of the ship are given in Fig. 10. The bandwidth is 10-20 GHz, the center frequency is set equal to 15 GHz, and the VV polarization is taken into account. The pulsewidth is set as 0.4 ns, and the incident and the azimuth angles are 45 • and 0 • , respectively. The comparison results in Fig. 10(b) also prove the effectiveness of the TDSBR for the EM scattering prediction of complex targets. Since the distances between each pair of parts of the ship are different, the scattering contribution arrives at the observation point at dissimilar times. Therefore, as shown in Fig. 10(a), the total scattering field is separated into multiple parts in time.
Besides, the accuracy of the TD method does not vary with the frequency since the integral results of (9) do not rely on the frequency. Here, the Nyquist sampling law is adopted, hence, the number of time sampling points varies with a variety ranges of frequencies, guaranteeing the calculation accuracy. The demonstrated plots in Fig. 11 indicate that the RCS results   show that the calculation accuracy does not change with the frequency when the sampling law is satisfied.
The composite scene of the complex ship target above the sea surface is illustrated in Fig. 12. The size of the sea surface is set as 150 m × 150 m , and the ship parameters are given in Fig. 9. The number of facets in this composite scene model is 45 390, and  therefore the longest edge of the patch is 0.707 m. The bandwidth is 10-20 GHz, the center frequency is 15 GHz, and the type of polarization is set as VV. The double Debye model is adopted to evaluate the sea surface dielectric parameters at the central frequency. Hence, the relative dielectric parameter of the sea is ε r = (41.75, 39.5) and its relative permeability is μ r = 1.0.
The pulsewidth is 0.4 ns. In Fig. 13(a) and (b), the incident and azimuth angles in order are set as 45 • and 0 • . Fig. 13(a) shows the backscattering field of a composite scene and the single ship, as presented in Fig. 12(a). Due to the existence of the sea surface, the echo delay of the composite scene is longer than that of the single ship. The plotted results in Fig. 13(a) indicate that the ship  has a strong scattering structure, so its echo intensity is much stronger than that of the sea surface. In order to better display the influence of the sea on the echo, two periods of time are amplified in the figure. The upper subgraph in Fig. 13(a) illustrates that the backscatter of the ship in the composite scene represents the main component of the total scattering, and the pulse shape is almost the same as that in the case of the ship alone. The lower subgraph in Fig. 13(a) indicates that the composite scattering echo exhibits a long duration, and the pulse shape is irregular, which is mainly caused by the random fluctuation of the sea surface. Besides, the comparison results in Fig. 13(b) prove the accuracy of the TD algorithm in predicting the composite scene scattering.
In order to explain the difference in the calculation time of the TDSBR method when the narrow pulse signal and the LFM signal are implemented as an incident source, the ship model shown in Fig. 9 are taken as examples. The bandwidth is set as 10-11 GHz, the center frequency is 10.5 GHz, and the polarization is taken as VV polarization. The relative dielectric parameter of the sea surface is ε r = (53. 26, 38.85) and the relative permeability is μ r = 1.0. The pulse duration of the modulated Gaussian pulse signal is 8 ns, and the number of sampling points is 57 372. In this case, for the LFM signal, four different pulse durations are chosen, which are 8 ns , 80 ns , 160 ns , and 1 μs. The corresponding sampling points are 57 372, 63 708, 70 748, and 144 668. The calculation results of the FDSBR method are used as the standard to evaluate the accuracy. The calculation results are shown in Fig. 14. The computation time and root-mean-square error are provided in Table I. According to the comparison data presented in Fig. 14 and Table I, it is difficult to ensure the calculation accuracy of employing the LFM signal as a signal source when the pulse duration is shorter. However, the computational time of the algorithm hugely magnifies with the growth of the pulse duration.
In practical applications, the LFM signal pulsewidth is generally longer to ensure the detection distance, so that it is hard to utilize the LFM signal pulse as the incident source of the TDSBR method for the SAR echo simulation.

B. SAR Images Simulation
In this part, we choose a cargo ship as the simulation object, and the maritime mobile service identity of the real cargo ship is 353 601 000. Its optical picture and geometric model picture are demonstrated in Fig. 15. The ship's length, width, and height in order are 183, 32, and 51 m, respectively. Real cargo ships have numerous complex and finely detailed structures, so the presented ship CAD model in Fig. 15(b) is only capable of guaranteeing a certain degree of accuracy with inevitable errors. We conduct SAR image simulations for the composite scene of the ship and sea surface. The composite scene is presented in Fig. 15(c), where Φ = 147 • represents the deflection angle of the ship. The sea surface size is 375 m × 375 m. The conditions are listed in Table II, which contains three resolution conditions. The number of facets in this composite scene model is 140 625, and the longest edge of the patch is 0.707 m. First, according to the calculation condition I, the TD method proposed in this article is employed to simulate the SAR image of the composite scene, and the simulation results are compared with those of the FD method and the measured SAR image from the Advanced Land Observation Satellite. The relative dielectric parameter of the sea is set as ε r = (73.57, 65.4) and its relative permeability is μ r = 1.0. The obtained results are shown in Fig. 16. The depicted results reveal that there are discrepancies between the results of the simulation method and those of the measured image, which is mainly due to the difference between the simulation ship model and the actual ship model, and the actual imaging conditions are more complex. However, on the whole, the simulation image of the TD method can reflect the structural characteristics similar to the measured image, indicating the effectiveness of the simulation method.
In order to illustrate the effectiveness of the TD method more specifically, cosine similarity is employed to calculate the similarity between the simulation results and the measured data, the calculation formula can be expressed as follows:  If X i and Y i are centralized (i.e., X = 0 ,Y = 0), then r could be expressed by According to (23), the similarity between the image generated by the TD method and the actual data is 91.3%. And the similarity between the image generated by the FD method and satellite data is 89.1%. Comparison shows that the image simulated by the TD method has a high degree of similarity.
In order to analyze the efficiency of the TD and FD method, further simulations are carried out based on conditions II and III. The relative dielectric parameter of the sea is set as ε r = (54.65, 38.57) and its relative permeability is μ r = 1.0. Fig. 17 shows the imaging results, in which Fig. 17(a) and (c) is simulated by the TD method, and Fig. 17(b) and (d) is generated by the FD method. According to the comparison results in Figs. 16 and 17, the imaging results of the TD method could better reflect most of the structure of the ship, while the calculation results of the FD method are more obvious to reflect the edge of the ship and the structure with strong scattering characteristics.
The times taken by the two methods to simulate the SAR images are listed in Table III. It can be seen from the comparison results that the efficiency advantage of the TD method is not apparent in generating low-resolution images. But the efficiency advantage is more apparent in generating high-resolution images since a greater number of slow-time sampling points are required. At each slow-time sampling point, the FD method should to calculate the scattering results for several frequency     I  COMPUTATION TIME AND RMSE   TABLE II  MAJOR PARAMETERS OF SAR IMAGE SIMULATION   TABLE III  COMPARISON OF SIMULATION TIME points, in turn significantly leads to rising the computational time. It is also worth noting that the time required by the TD method to calculate the broadband scattering inconspicuously magnifies for a certain level of the echo delay. Therefore, the proposed method on the basis of TDSBR to simulate the SAR echo is more suitable for high-resolution SAR image simulations.

V. CONCLUSION
In this article, a method combining transfer function and RFPC is proposed to extend the TDSBR method to SAR echo simulation. First, the correctness of the TDSBR algorithm is verified by the simulation results of the trihedral reflector, simple ship model, and ship-sea composite model. By comparing the computational times of the TDSBR method in the case of employing the modulated Gaussian pulse signal and the LFM signal as the incident source, it is shown that if the LFM signal is directly adopted as the source of the TDSBR method, the corresponding pulse duration would be longer than that of the narrow pulse signal. For actual SAR echoes, the LFM signal usually needs a longer pulse duration, as a result, this method is not suitable for SAR echo simulation. Further, the simulation results based on the LFM signal as an incident source of the TDSBR are not given in the part of the SAR images simulation since the calculation time is unbearable. The effectiveness of the proposed method for the simulation of SAR echo is also demonstrated by comparing the measured image results and those of the TD method. The comparison of the computational times proves that the proposed method exhibits better computational efficiency.