A Novel Hybrid Algorithm for Transient Near-Field Scattering from Rough Surface in 2-D Case

A novel hybrid algorithm is proposed to reduce the computation cost of the ﬁnite-diﬀerence time-domain (FDTD) method in calculating the transient near-ﬁeld scattering from rough surface. The scattering problem is split into the FDTD calculation of equivalent sources on the contour enclosing rough surface and the calculation of the near-ﬁeld radiation with the two-dimensional (2-D) time-domain Huygens’ principle. The radiation ﬁelds are found from a surface integral of the temporal convolution for which the direct numerical integration of the convolution is computationally expensive. In this paper, the 2-D time-domain Green’s function as the convolution kernel is approximated with a sum of exponential terms by using the Prony’s method. Then, the semianalytical recursive convolution (SARC) approach is applied to complete the update of the near-ﬁeld radiation. Compared with the traditional FDTD, this hybrid algorithm can signiﬁcantly reduce the memory usage and run time, especially for the large distance between the rough surface and observation point.


Introduction
With the extensive application of the wide-band radar, the transient electromagnetic (EM) scattering from the rough surface has gained increasing attention, which plays important role in the field of radar surveillance, EM imagining, target identification, and stealth, etc. [1][2][3][4]. e previous researches are mainly focused on the simulation of the EM scattering from rough surface in the far-field zone [5][6][7]. In fact, the near-field scattering problem is the widespread presence of the target working above background, such as the low-altitude cable over ground or sea. Only by mastering the near-field scattering characteristics from rough surface, we can systematically analyze the EM interference from background and obtain the relatively complete coupling result of target. erefore, it is of critical importance to study an efficient time-domain algorithm on the near-field scattering for the coupling analysis of the low-altitude target and the protection of the sensitive electronic device.
In the past decade, the numerical method has been widely used in the EM analysis owning to the rapid development of the high performance computer [8][9][10][11]. Among them, the finite-difference time-domain (FDTD) method is a powerful and popular technique for numerically solving the transient field problem, and it has achieved successful applications in the far-field scattering from rough surface [12][13][14]. However, when the FDTD is used to simulate the near-field scattering, a brute-force way is to enlarge the computational domain to enclose the rough surface and observation point, which will lead a heavy burden to the computation resource available. Moreover, the farther the near-zone observation point is from rough surface, the heavier the computation burden will be. By now, some hybrid algorithms combining the FDTD and analytical-approximate method have been proposed to reduce the FDTD space lattice and improve the computational efficiency. In [15], the FDTD associated with the Kirchhoff surface integral was proposed to calculate the near-field response of the dipole antenna in the three-dimensional (3-D) case, where the dipole antenna was assigned in the FDTD region, and the near-field response outside the FDTD space was calculated by the Kirchhoff surface integral. For the two-dimensional (2-D) problem, by means of the 2-D timedomain Huygens' principle, the temporal convolution integral needs to be solved to obtain the near-field response exterior to the FDTD region, which complicates the computation even more. e direct integration of the temporal convolution is computationally expensive, and hence it is applicable for the EM analysis of the small model such as the infinitely long line source in [16,17]. Xu et al. in [18] employed the trapezoidal approximation to solve the temporal convolution, and then the near-field scattering exterior to the FDTD region was computed by applying the scheme of treating the radiation field as the equivalent cylindrical wave, which can be accurately used for the target with cylinder shape. However, for the large rough surface model, the characteristics of the structure and scattering wave make the simulation of the near-field scattering with these approaches difficult. In recent years, a fruitful semianalytical recursive convolution (SARC) approach was introduced in the dispersive medium problem, realizing the high-efficient computation of the convolution integral [19,20]. But when the 2-D time-domain Green's function is acted as the convolution kernel, the recursive operation is difficult to realize due to its specific structure. In this paper, this difficulty is overcome by applying the Prony's method to approximate the convolution kernel with a sum of exponential terms, which avoids the complex analysis of the convergence of the series expansion scheme in [21]. e hybrid algorithm combining the FDTD with the 2-D time-domain Huygens' principle is proposed to calculate the transient near-field scattering from the rough surface, in which the normal part of the convolution integral is calculated by the SARC approach and the singular part is done using the linear interpolation. Because a great deal of FDTD meshes between the rough surface and near-zone observation point are removed, the memory usage and run time in the proposed algorithm are dramatically reduced than those in the traditional FDTD method. Meanwhile, this hybrid algorithm retains the unique advantages of the FDTD method. e rest of this paper is organized as follows: in Section 2, the hybrid algorithm combining the FDTD with the 2-D time-domain Huygens' principle is elaborated, including the Prony approximation of the convolution kernel and the SARC approach for solving the near-field radiation of the equivalent sources. In Section 3, the accuracy and efficiency of the proposed hybrid algorithm are verified through the comparison of the numerical results by the hybrid algorithm and by the traditional FDTD method. Section 4 draws the conclusions of this paper and offers the recommendations for future work.

Hybrid Algorithm for Near-Field Scattering from Rough Surface
For the FDTD calculation of the transient near-field scattering, a common way is to expand the FDTD region to enclose the rough surface and the near-zone observation point, as shown in Figure 1(a). With the increase in the height of observation point, the memory usage and run time will sharply increase, thus often making it difficult to realize the FDTD calculation of the near-field scattering. For the proposed hybrid algorithm, however, the FDTD region only contains the rough surface as indicated in Figure 1(b); the near-field response at the observation point P is regarded as the radiation of the equivalent sources on the output boundary of the FDTD region. Comparing Figures 1(a) and 1(b), it can be seen that the computational domain in the hybrid algorithm is much smaller than that in the FDTD method, especially in the case with high observation point.

eory of Near-Field Transformation in 2-D Case.
Assume that E and H are respectively the electric and magnetic fields on the output boundary surface of the FDTD domain. According to the equivalence principle, the surface electric and magnetic currents (J, M) and charges (ρ s , ρ m ) can be obtained by where n denotes the unit normal vector on the output boundary. For the 2-D problem, the vector and scalar potentials can be expressed with the 2-D time-domain Green's function as Here, μ and ε represent the permeability and permittivity, R � |r − r ′ | is the distance from the observation point r to source field r ′ , and c is the speed of the EM wave. In accordance with the vector and scalar potentials, the 2-D radiation field exterior to the FDTD region can be written as Substituting equation (2) into equation (3) and shifting the time reference backward by R/c, the formals for calculating the 2-D radiation field of the equivalent sources can be obtained from [21], in which, for the TM case, are 2 International Journal of Antennas and Propagation e corresponding formulas in the TE mode can be derived with the principle of duality, which are not described here. It is obvious that the convolution kernel in equation (4) has a square-root singularity at t ′ � t. An efficient way in addressing the singularity is to approximately calculate the convolution integral in the interval [t − Δt, t] using the linear interpolation. Assuming t − t ′ � τ, the time derivative of the equivalent source, e.g., M(r, t ′ ), in the time-reversed and offset form can be given by For instance, substituting (5) into the second term of (4a), we have A similar procedure can be implemented for other terms in equation (4) to obtain the corresponding approximate results in the interval [t − Δt, t]. e convolution integral in the interval [0, t − Δt] is computed by using the high-efficient SARC approach, as will be elaborated later.

Prony Approximation of Convolution Kernel.
Although equation (4) is applicable for the computer execution, the direct integration of it is such cumbersome as to be impractical, due to the need for storing and processing all past values of the time derivative of the equivalent source. e computational complexity of this scheme scales as O(N 2 max ), where N max is the total time step. In this paper, the SARC approach will be applied to decrease the complexity of the convolution integral. It can be found from equation (4) that the convolution kernel mainly consists of two forms: 1/ �������������� � (t + R/c) 2 − R 2 /c 2 and (t + R/c)/ �������������� � (t + R/c) 2 − R 2 /c 2 , which make it difficult to execute the recursive operation of the convolution. Based on the idea of the signal estimation in the highspeed circuit [22], the Prony's method is used to approximate the convolution kernel with a sum of the exponential terms: International Journal of Antennas and Propagation 3 Here, f(N) is the discrete sample of the convolution kernel at t � NΔt, where N denotes the equally spaced time step, Q represents the total number of poles, and A i and C i represent the coefficient and amplitude of the ith pole. e where N m − Q > Q. e roots α 1 , α 2 , . . . , α Q can be determined by solving this oversampled set of linear equations with the least-squares approach. By virtue of equation (8), the μ can be found as the roots of the following Qth degree equation: If the complex conjugate roots arise in the solution to equation (9), use their magnitude to substitute for the complex pair of μ, and cancel any duplicate. Substituting μ i to equation (7) and introducing N m discrete samples of the convolution kernel, a new set of the linear equations with the Q values of the C i can be established as follows: where N m > Q. e coefficients C 1 , C 2 , . . . , C Q can be obtained by using the least-squares approach to equation (10).

���������������
(τ + R/c) 2 − R 2 /c 2 are respectively expressed as convolution kernels 1 and 2 for convenience. e total numbers of the samples and poles are selected to be 30 and 10, respectively. It can be observed from Figure 2 that the nearly identical results are produced with the two methods, implying the reliability of the Prony's method.

Semianalytical Recursive Convolution Approach.
e RASC approach is a fruitful scheme for the calculation of the convolution integral, which is based on the idea of the efficient analysis of the transfer response in the digital signal processing (DSP) [23]. For a linear time-invariance and dynamical system, as shown in Figure 3, the system response y(t) can be expressed as the convolution of the input signal x(t) with the impulse response h(t). We have In general, the impulse response in an asymptotically stable system is of the form, where q stands for the number of poles, and U(t) is the switch function: Substituting equation (12) into equation (11), the system response at t � nΔt can be obtained as e system response for the qth pole is Dividing the integral of equation (15) into two intervals, [0, (n − 1)Δt] and [(n − 1)Δt, nΔt], we get to the recursive formula of system response as

International Journal of Antennas and Propagation
In this paper, the near-field radiated by the equivalent source is analogous to the system response in the DSP. e convolution kernel can be regard as the impulse response after it is approximated with a sum of exponential terms. e time derivative of the equivalent source corresponds to the input signal. us, the near-field radiation can be updated by applying the SARC approach. Next, focus on the SARC approach for evaluating the radiation field in equation (4), where the representative convolutions are of the forms, For equation (17a) as example, using the Prony's method mentioned above can approximate the convolution kernel as According to equation (16), the temporal convolution in the interval [0, (n − 1)Δt] can be expressed as where W(t) can be approximated in accordance with the linear interpolation approach as follows: Substituting equation (20) into equation (19) and analytically solving the integral over the interval [(n − 2)Δt, (n − 1)Δt], we find the following field-update equation: where From equation (21), it is clear that the SARC approach has a high efficiency in calculating the near-field radiation.
is is because this approach avoids the requirement for storing and processing of the complete history of the equivalent source on the output boundary of the FDTD region, and its computational complexity is reduced as only O(N max ). Beyond that, the SARC approach has a better accuracy than the general recursive approach, due to the analytical calculation of the convolution integral in the discrete time interval. Figure 4 shows the geometry of the near-field scattering from rough surface. e observation point P is located above the rough surface, and θ i denotes the incident angle. Without loss of generality, the one-dimensional Weierestrass-Mandelbrot function is used here to generate the rough surface [24]. In the numerical calculation, the length of the rough surface is taken to be L � 160 m and the spatial increment is taken to be Δ � Δx � Δy � 0.05 m to allow for reasonable resolution of the rough surface model discretized with the staircase approximation, and the time increment is set to be Δt � Δ/2c in accordance with Courant stability condition [22]. e 2-D FDTD space lattice is terminated by the uniaxial perfectly matched layer (UPML) [25], and its thickness is chosen to be 10Δ. Consider the rough surface is exposed to the Gaussian pulse, whose expression is

Results and Discussion
e parameters of the incident wave are adopted as follows: t 0 � 60Δt and τ � 0.8t 0 .
To illustrate the accuracy and efficiency of the proposed hybrid algorithm, Figures 5 and 6 present the comparison of the near-field scattering results by the hybrid algorithm with those by the traditional FDTD method for the TE and TM polarized conditions, respectively. In Figures 5(a)  is result demonstrates the accuracy of the hybrid algorithm proposed in this paper. Beyond that, hybrid algorithm can dramatically lessen the computation cost compared with the traditional FDTD method. For the TE polarized wave incidence and 4000 time-steps updating as example, the computation costs required for two methods are provided in Table 1. When the height of the near-zone observation point is h � 5 m, the number of FDTD meshes, the memory usage, and the run time in the hybrid algorithm are reduced to about 72%, 76%, and 74% of those in the traditional FDTD method. For h � 16 m, the key factors above using the hybrid algorithm are only 39%, 42%, and 33% of those using the traditional method, respectively. It is obvious that the hybrid algorithm becomes more efficient with the increase in the height of the near-zone observation point. is is because that the FDTD region in the traditional method has to be greatly enlarged to enclose the rough surface and observation point, whereas that in the hybrid algorithm remains unchanged when the observation point rises. In addition, we can find from Table 1 that the proposed algorithm differs little from the traditional FDTD method in the run time for h � 2.5 m. is is owing to the fact that, in this case of very low observation point, the advantage of the hybrid algorithm in reducing the FDTD meshes is almost offset by its calculation of the near-field radiation. As the height of the near-zone observation point increases, the proposed hybrid algorithm becomes more favorable in the computation efficiency.

Conclusions
In this paper, a novel hybrid algorithm has been proposed to calculate the transient near-field scattering from rough surface for both TE and TM polarization. is hybrid algorithm is based on the combination of the FDTD and the 2-D time-domain Huygens' principle, and it is the development of the traditional FDTD method. In the hybrid algorithm, the FDTD is only performed in the computational domain containing the rough surface, which makes a large number of FDTD meshes eliminated in the space between the rough surface and observation point. e proposed hybrid algorithm is validated by comparing the near-field scattered waveforms of the rough surface for both the hybrid algorithm and the traditional FDTD method. e numerical results show that the hybrid algorithm is an accurate and efficient technique in the calculation of the transient near-field scattering from the rough surface. Moreover, the relative efficiency of the hybrid algorithm becomes higher with the increase in the height of the observation point. In the future, the hybrid algorithm based on SARC approach will be developed to study the transient composite scattering and coupling of the dispersive target and background.

Data Availability
e simulation data used to support the finding of this study are included within the article.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.