Multitargets Orientation Technique Based on Reflection Characteristic Analysis Using an Inverse Diffraction Parabolic Equation

Abstract

In this article, the inverse diffraction parabolic equation (IDPE) model based on the finite difference method is proposed, which is first applied in the multiple nonradiation targets orientation technology. In principle, the electromagnetic signal propagating in the transmission path will produce a reflected signal back to the source end while encountering the discontinuous objects. The distribution of the reflection or refraction intensity is directly associated with the distances and heights of the objects, so the location can be determined by means of analyzing the distribution. Here, according to the profile data of field intensity at the source end, the distribution of backward propagating electromagnetic waves are calculated rapidly by the IDPE. Then, the local extreme searching method is applied to search the coordinate of the convergence point of field intensity and the positions of multiple objects are finally determined. The piecewise linear function is used to model the irregular terrain. The influence of discontinuous terrain slopes on the false alarm probability of objects localization is also analyzed. The results show that the localization accuracy of the IDPE algorithm is affected by multiple factors, such as the radio frequency and sampling interval of field intensity. It is proved that the IDPE is a novel and efficient algorithm for multiple nonradiation targets orientation technology in long-range complicated terrain environment.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

In recent years, with the rapid progress of the electromagnetic theory and digital technology, the targets’ orientation based on wireless electromagnetic waves (EWs) is currently widely used in many areas, such as global positioning systems, autonomous driving machines, and radar systems. However, the EW propagation process is very complex; thus, the multiple targets’ location in long-range complicated terrain environment has always been a challenge.

At present, some full-wave algorithms, such as physical optics (PO) [1] and ray tracing (RT) [2], were already quite mature. The shooting and bouncing ray (SBR) [3] method is usually used in radar cross-section computation of cavity scattering. The finite-difference time-domain (FDTD) method is applied to the antenna design and electromagnetic scattering analysis [4]. However, the traditional geometric algorithms are only suitable for simple scenarios. The FDTD method requires to finely mesh the calculation region, which causes the lower computational efficiency. Thus, the above methods cannot satisfy the requirements of the accuracy and efficiency of large-scale EW propagation prediction in complicated scenarios.

The parabolic equation (PE) is a forward iterative algorithm under paraxial approximation, which has the advantages of greater accuracy and efficiency than traditional geometric algorithms [5, 6]. The hybrid method of FDTD and PE has emerged [7]. The conventional PE can be solved by the split-step Fourier transform (SSFT) [8–10]. But SSFT has a lack of flexibility for boundary modeling, which is not suitable for irregular terrain. The finite difference method (FDM) is attractive and convenient because of their flexibility for boundary modeling and greater propagation angle limits [11].

The inverse diffraction parabolic equation (IDPE) is an inverse iterative algorithm based on PE [12]. The main idea is to fast-search radiation sources by backward deduction according to the distribution of the forward-propagating EWs. The studies have shown that the IDPE method is a useful tool to apply in the radiation sources orientation [13]. However, conventional IDPE is based on the SSFT method, which has lower localization accuracy for irregular terrain. Furthermore, there are few studies about the application of the IDPE algorithm in nonradiation targets orientation.

In this article, the IDPE based on the FDM is proposed, which is first applied in the multiple nonradiation targets orientation. In principle, the electromagnetic signal propagating in the transmission path will produce a reflected signal back to the source end while encountering the discontinuous objects. The distribution of the reflection or refraction intensity is directly associated with the distances and heights of the objects, so the position can be determined by means of analyzing the distribution. Here, according to the profile data of field intensity at the source end, the distribution of backward propagating EWs is calculated rapidly by the IDPE. Then, the local extreme searching method is applied to search the coordinate of the convergence point of the field intensity, and the positions of multiple objects are finally determined. The piecewise linear function is used to model the irregular terrain. The influence of discontinuous terrain slopes on the false alarm probability of objects location is also analyzed. The results show that the positioning accuracy is affected by multiple factors, such as EWs frequency and sampling intervals of field intensity. The lower the frequency and sampling intervals, the higher the position accuracy. It is proved that the IDPE is a novel and fast technique for multiple nonradiation targets orientation in an irregular terrain environment.

2. Inverse Diffraction Parabolic Equation Based on the Finite Difference Method

A detailed description of the inverse diffraction parabolic equation (IDPE) based on the finite difference method will be proposed in this section [14].

Assuming the two-dimensional Cartesian coordinates (x, z), the electric or magnetic fields are independent of the horizontal range y direction, where x represents the horizontal distance, z represents vertical height, and exp^−iωt represents time dependence of the fields.

The electric field E only has one nonzero component $E_y$ for horizontal polarization, the field component ψ is defined by$$\begin{array}{}\text{(1)}& \psi x,z=E_{y}x,z.\end{array}$$

The magnetic field H only has one nonzero component $H_y$ for vertical polarization, and the field component ψ is defined by$$\begin{array}{}\text{(2)}& \psi x,z=H_{y}x,z.\end{array}$$

The two-dimensional vector wave equation is given as follows:$$\begin{array}{}\text{(3)}& \frac{{\partial }^2\psi }{\partial x^2}+\frac{{\partial }^2\psi }{\partial z^2}+k^2{n}^2\psi =0,\end{array}$$where n is the refractive index for free space and k is the radio wave number. Assuming the reduced functions associated with the paraxial direction x is introduced,$$\begin{array}{}\text{(4)}& ux,z=e^{-ikz}\psi x,z.\end{array}$$

Then, the two-dimensional scalar wave equation is given as follows:$$\begin{array}{}\text{(5)}& \frac{{\partial }^{2}u}{\partial x^2}+2ik\frac{\partial u}{\partial x}+\frac{{\partial }^{2}u}{\partial z^2}+k^2{m}^2-1u=0,\end{array}$$where u represents a scalar component of the electric field E for horizontal polarization, or the magnetic field H for vertical polarization, and m is the modified refractive index defined as follows:$$\begin{array}{}\text{(6)}& mx,z=nx,z-1+\frac{z}{a_e}\times {10}^6,\end{array}$$where a_e is the radius of the earth. The scalar wave equation can be split into the two terms, including the forward field u₊ and the backward field u₋. Also, equation (5) can be formally written as follows:$$\begin{array}{}\text{(7a)}& \frac{\partial u_{+}}{\partial x}+ik1-Q{u}_{+}=0,\text{(7b)}& \frac{\partial u_{-}}{\partial x}+ik1+Q{u}_{-}=0.\end{array}$$

Equation (7a) represents the forward PE, and the equation (7b) represents the backward PE.

The pseudodifferential operator Q is given as follows:$$\begin{array}{}\text{(8)}& Q=\sqrt{1+Z}.\end{array}$$

The symbol Z can be defined as follows:$$\begin{array}{}\text{(9)}& Z=\frac{1}{k^2}\frac{{\partial }^2}{\partial z^2}+m^2-1.\end{array}$$

The pseudodifferential operator Q can be approximated using the rational form as follows:$$\begin{array}{}\text{(10)}& Q=\sqrt{1+Z}\approx \frac{{\chi }_1+{\chi }_{2}Z}{{\chi }_3+{\chi }_{4}Z}.\end{array}$$

Here, the Greene approximation coefficients are χ₁ = 0.99987, χ₂ = 0.79624, χ₃ = 1.00, and χ₄ = 0.30102. The Claerbout approximation coefficients are χ₁ = 1, χ₂ = 0.75, χ₃ = 1.00, and χ₄ = 0.25 [15].

Then, the forward PE based on the rational form approximation equation (10) is defined by$$\begin{array}{c}\text{(11)}& \frac{{\partial }^3{u}_{+}}{\partial x\partial z^2}+ik1-\frac{{\chi }_2}{{\chi }_4}\frac{{\partial }^2{u}_{+}}{\partial z^2}+k^2\frac{{\chi }_3}{{\chi }_4}+m^2-1\frac{\partial u_{+}}{\partial x}+\frac{i{k}^3}{{\chi }_4}{\chi }_3-{\chi }_1+{\chi }_4-{\chi }_2{m}^2-1u_{+}=0.\end{array}$$

Also, the backward PE based on the rational form approximation equation (10) is defined by$$\begin{array}{c}\text{(12)}& \frac{{\partial }^3{u}_{-}}{\partial x\partial z^2}+ik1+\frac{{\chi }_2}{{\chi }_4}\frac{{\partial }^2{u}_{-}}{\partial z^2}+k^2\frac{{\chi }_3}{{\chi }_4}+m^2-1\frac{\partial u_{-}}{\partial x}+\frac{i{k}^3}{{\chi }_4}{\chi }_3+{\chi }_1+{\chi }_4+{\chi }_2{m}^2-1u_{-}=0.\end{array}$$

How to efficiently calculate the above PE is essential to the EWs prediction precision. The common numerical computing methods are the split-step Fourier transform (SSFT) method and the finite difference method (FDM). The FDM is superior in the accuracy and convenience and has greater angles limits, which is better suited for the irregular terrain environment.

In Figure 1, the three-point FDM scheme is proposed to calculate the PE model. Supposing u^m−1 represents the field at the previous horizontal step, u^m represents the field at the next horizontal step, and the virtual coordinate point ($x_{\kappa }$, $z_j$) is the middle point in the horizontal x direction defined by$$\begin{array}{}\text{(13)}& u{x}_{\kappa },z_j=\frac{u_j^m+u_j^{m-1}}{2}.\end{array}$$

[figure(s) omitted; refer to PDF]

For the forward PE, the first-order difference in the positive x direction is as follows:$$\begin{array}{}\text{(14)}& \frac{\partial u_{+}{x}_{\kappa },z_j}{\partial x}=\frac{u_j^m-u_j^{m-1}}{\Delta x}.\end{array}$$

For the backward PE, the first-order difference in the negative x direction is as follows:$$\begin{array}{}\text{(15)}& \frac{\partial u_{-}{x}_{\kappa },z_j}{\partial x}=\frac{u_j^{m-1}-u_j^m}{\Delta x}.\end{array}$$

The second-order difference forms in the z direction are as follows:$$\begin{array}{c}\text{(16)}& \frac{{\partial }^{2}u{x}_{\kappa },z_j}{\partial z^2}=\frac{1}{2\Delta z^2}{u}_{j+1}^m-2u_j^m+u_{j-1}^m+\frac{1}{2\Delta z^2}{u}_{j+1}^{m-1}-2u_j^{m-1}+u_{j-1}^{m-1}.\end{array}$$

The first-order difference in the positive x direction and second-order difference in the z direction are as follows:$$\begin{array}{c}\text{(17)}& \frac{{\partial }^3{u}_{+}{x}_{\kappa },z_j}{\partial x\partial z^2}=\frac{1}{\Delta x\Delta z^2}{u}_{j+1}^m-2u_j^m+u_{j-1}^m-\frac{1}{\Delta x\Delta z^2}{u}_{j+1}^{m-1}-2u_j^{m-1}+u_{j-1}^{m-1},\end{array}$$and the first-order difference in the negative x direction and second-order difference in the z direction are as follows:$$\begin{array}{c}\text{(18)}& \frac{{\partial }^3{u}_{-}{x}_{\kappa },z_j}{\partial x\partial z^2}=\frac{1}{\Delta x\Delta z^2}{u}_{j+1}^m-2u_j^m+u_{j-1}^m-\frac{1}{\Delta x\Delta z^2}{u}_{j+1}^{m-1}-2u_j^{m-1}+u_{j-1}^{m-1},\end{array}$$where j = 1, …, Z.

Using the above difference scheme in (14), (16), and (17) for (11), the forward rational approximation PE model can be expressed by a matrix as follows:$$\begin{array}{}\text{(19)}& A{u}_{j-1}^m+C_j^m{u}_j^m+A{u}_{j+1}^m=\tilde{A}{u}_{j-1}^{m-1}+{\tilde{C}}_j^m{u}_j^{m-1}+\tilde{A}{u}_{j+1}^{m-1},\end{array}$$where $\tilde{A}=A-x$, $\tilde{B}=B-x$, and ${\tilde{C}}_j^m=C_j^m-x$. The coefficients in equation (19) are shown in the following:$$\begin{array}{}\text{(20)}& \begin{array}{c}A=2{\chi }_4+ik\Delta x{\chi }_4-{\chi }_2,\\ B=2{\chi }_3+ik\Delta x{\chi }_3-{\chi }_1,\\ C_j^m=-2A+k^2\Delta z^{2}A{m_j^m}^2-1+B.\end{array}\end{array}$$

We use the Leontovich lower boundary condition on the surface as follows:$$\begin{array}{}\text{(21)}& \frac{\partial x_{\kappa },0}{\partial z}=-ik\eta xu{x}_{\kappa },0,\end{array}$$where the surface impedance η is given by$$\begin{array}{}\text{(22)}& \eta =\begin{array}{cc}\sqrt{\epsilon -1},& \text{horizontal} \text{polarization},\\ \frac{\sqrt{\epsilon -1}}{\epsilon },& \text{vertical} \text{polarization}.\end{array}\end{array}$$

The piecewise linear function is adopted to model the irregular terrain. Assuming the old coordinate system $v,w$, the new coordinate system is defined by (x, z) as follows:$$\begin{array}{}\text{(23)}& x=v,z=w-fv,\end{array}$$where f $v$ represents the terrain function correlated with the terrain slope S.

Then, the tridiagonal structure matrix of forward PE model is obtained as follows:$$\begin{array}{c}\text{(24)}& \begin{array}{ccccccc}G& A& 0& 0& 0& \cdots & 0\\ A& C_1^m& A& 0& 0& \cdots & 0\\ 0& A& C_2^m& A& 0& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \cdots & 0& A& C_{Z-1}^m& A\\ 0& 0& \cdots & 0& 0& A& C_Z^m\end{array}\begin{array}{c}{u}_0^m\\ u_1^m\\ u_2^m\\ ⋮\\ u_{Z-1}^m\\ u_Z^m\end{array}=\begin{array}{ccccccc}\tilde{G}& \tilde{A}& 0& 0& 0& \cdots & 0\\ \tilde{A}& {\tilde{C}}_1^m& \tilde{A}& 0& 0& \cdots & 0\\ 0& \tilde{A}& {\tilde{C}}_2^m& \tilde{A}& 0& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \cdots & 0& \tilde{A}& {\tilde{C}}_{Z-1}^m& \tilde{A}\\ 0& 0& \cdots & 0& 0& \tilde{A}& {\tilde{C}}_Z^m\end{array}\begin{array}{c}{u}_0^{m-1}\\ u_1^{m-1}\\ u_2^{m-1}\\ ⋮\\ u_{Z-1}^{m-1}\\ u_Z^{m-1}\end{array}.\end{array}$$

The coefficient of boundary G is given by$$\begin{array}{}\text{(25)}& G=\frac{C_0^m}{2}+A\xi ,\end{array}$$where ξ = ikηΔz. Equation (24) is the forward PE based on the three-point scheme FDM.

Using the difference scheme in (15), (16), (18), for (12), the inverse rational approximation PE model can be written as follows:$$\begin{array}{}\text{(26)}& D{u}_{j-1}^{m-1}+F_j^m{u}_j^{m-1}+D{u}_{j+1}^{m-1}=\tilde{D}{u}_{j-1}^m+{\tilde{F}}_j^m{u}_j^m+\tilde{D}{u}_{j+1}^m.\end{array}$$

The coefficients in equation (26) are shown in the following:$$\begin{array}{}\text{(27)}& \begin{array}{c}D=2{\chi }_4+ik\Delta x{\chi }_4+{\chi }_2,\\ E=2{\chi }_3+ik\Delta x{\chi }_3+{\chi }_1,\\ F_j^m=-2D+k^2\Delta z^{2}E{m_j^m}^2-1+D.\end{array}\end{array}$$

Also, the tridiagonal structure matrix of backward PE model is obtained as follows:$$\begin{array}{c}\text{(28)}& \begin{array}{ccccccc}H& D& 0& 0& 0& \cdots & 0\\ D& F_1^m& D& 0& 0& \cdots & 0\\ 0& D& F_2^m& D& 0& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \cdots & 0& D& F_{Z-1}^m& D\\ 0& 0& \cdots & 0& 0& D& F_Z^m\end{array}\begin{array}{c}{u}_0^{m-1}\\ u_1^{m-1}\\ u_2^{m-1}\\ ⋮\\ u_{Z-1}^{m-1}\\ u_Z^{m-1}\end{array}=\begin{array}{ccccccc}\tilde{H}& \tilde{D}& 0& 0& 0& \cdots & 0\\ \tilde{D}& {\tilde{F}}_1^m& \tilde{D}& 0& 0& \cdots & 0\\ 0& \tilde{D}& {\tilde{F}}_2^m& \tilde{D}& 0& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \cdots & 0& \tilde{D}& {\tilde{F}}_{Z-1}^m& \tilde{D}\\ 0& 0& \cdots & 0& 0& \tilde{D}& {\tilde{F}}_Z^m\end{array}\begin{array}{c}{u}_0^m\\ u_1^m\\ u_2^m\\ ⋮\\ u_{Z-1}^m\\ u_Z^m\end{array}.\end{array}$$

The coefficient of boundary H is given by$$\begin{array}{}\text{(29)}& H=\frac{F_0^m}{2}+D\xi .\end{array}$$

The above tridiagonal matrix equation (28) is the IDPE based on the three-point scheme FDM.

In Figure 2, supposing the emitter located at a distance of 0 m, the forward EWs propagate in the positive x direction at step 1 and step 2. It occurs the reflection EWs when meets the first obstacle at a distance of $x_1$, the backward EWs start propagating along with the negative x direction, and the forward EWs continue to propagate on its way at step 3 and step 4. The backward initial field is given as follows [16]:$$\begin{array}{}\text{(30)}& u_{-}{x}_1,z_j=\begin{array}{cc}-u_{+}{x}_1,z_j{e}^{2ik{x}_1},& \text{for} 0<z_j<z_1,\\ 0,& \text{for} z_j>z_1.\end{array}\end{array}$$

[figure(s) omitted; refer to PDF]

When the forward EWs meets the second obstacle at a distance of $x_2$, it produce the reflection EWs again and backward EWs propagate along with the negative x direction at step 5–8. The backward initial field is given as follows:$$\begin{array}{}\text{(31)}& u_{-}{x}_2,z_j=\begin{array}{cc}-u_{+}{x}_2,z_j{e}^{2ik{x}_2},& \text{for} 0<z_j<z_2,\\ 0,& \text{for} z_j>z_2.\end{array}\end{array}$$

The forward diffraction EWs appear and propagate along with the positive x direction at steps 5 and 6.

The total backward EWs is the sum of forward EWs and all the backward EWs as follows:$$\begin{array}{}\text{(32)}& ux,z={\displaystyle \sum }_{n=1}^m{e}^{-ikx}{u}_{-}^{n}x,z+e^{ikx}{u}_{+}x,z.\end{array}$$

The radio propagation prediction can be achieved by the forward PE model. The distribution characteristic of EWs for the entire space can be calculated from step 1 to step 8.

The target location is based on the distribution of reflection EWs. According to the profile sampling data of field intensity at the source end, the reflection distribution of backward EWs can be effectively predicted using the IDPE model from step 9 to step 14.

3. Numerical Results and Discussion

In this section, the comparison between the PE model and the traditional geometric algorithms is firstly carried out in order to verify the efficiency and superiority of the PE model. Then, the numerous simulation experiments for single and multiple rectangles objects localization on flat or irregular terrain are conducted. Finally, the influence of discontinuous terrain slopes on the false alarm probability of objects location is analyzed. Also, the effects of EWs frequency and sampling intervals of field intensity on the position accuracy are investigated.

First, the PE model is compared with the dual-ray model and Miller–Brown model in Figure 3. The simulation conditions are set as follows: In the standard atmosphere, the transmitter is a Gaussian horizontally polarized (HP) antenna at a height of 30 m. The frequency is 400 MHz, and the half-power beam width is 10° with the inclination angle of 0°. The maximum distance and height are set to be 4 km and 1 km, respectively. The horizontal step ∆x and vertical step ∆z are both set to be one-half of the wavelengths. The Tukey window function is designed to be the 30% of total height to eliminate the reflected EWs coming from the upper boundary. The terrain is set to be Sine structure with an angular frequency of 0.0109° and amplitude of 3.0783 m. The relative permittivity is ${\epsilon }_g$ = 30, the conductivity is ${\sigma }_g$ = 0.01 S/m [17].

[figure(s) omitted; refer to PDF]

The results show that the curves of the three models are basically consistent in Figure 3. But the curve of the PE model shows the slight oscillation. This is because the geometric structure of rough ground surface will generate strong reflection and refraction effects on EWs. The dual-ray model and Miller–Brown model cannot simulate accurately the distribution of EWs on complicated terrain. Furthermore, the time for each simulation experiments based on the PE is within 5 min in this section, which is significantly faster than the FDTD method. In general, the accuracy and efficiency of the proposed method based on the PE model are better than those of other algorithms.

Second, supposing the single rectangle object is located at a distance of 5 km, which is 100 meters wide and 15 meters high. The simulation results based on IDPE are shown in Figures 4–7. The frequency is 0.9 GHz at a height of 20 m, and the half-power beam width is 3°. The maximum distance and height are set to be 10 km and 600 m. The surface ground is set to be the flat ground or irregular piecewise linear terrain. Other simulation conditions are unchanged.

[figure(s) omitted; refer to PDF]

In Figure 4, the distribution of the forward EWs are calculated by the positive PE model and it is shown that the irregular terrain causes more reflection and diffraction EWs in Figure 4(b). The reflection EWs affects the forward EWs distribution for the whole region.

In Figure 5, using the reverse-phase values of forward field intensity at the distance of 5 km as the electromagnetic profile data of the initial field of backward EWs, the backward EWs distribution is also calculated along with the negative x direction by the positive PE model. Because of the difference of initial field of backward EWs, combined with the irregular terrain in front of the obstacle, the backward EWs distribution characteristic appears to have an obvious difference.

[figure(s) omitted; refer to PDF]

In Figure 6, supposing the receiver located at a distance of 0 m, and the sampling intervals of received field intensity are set to be one wavelength. Using the received field intensity values as the electromagnetic profile data of initial field of IDPE model, the inverse diffraction EWs distribution is calculated along with the positive x direction based on the IDPE rapidly. The total computing time is 212. 6 seconds. It is found that the IDPE method can accurately and fast model the reflection characteristic of backward EWs for whole calculation region. Also, the distribution of the inverse diffraction EWs presents the symmetrical distributed situation. The horizontal position of the axis of symmetry is the target location.

[figure(s) omitted; refer to PDF]

In Figure 7, the maximum value of inverse diffraction EWs at every marching step is computed first, and the maximum field values curve along with the horizontal x direction is shown. We found that the normalized field intensity of inverse diffraction EWs on flat ground is greater than random terrain. That is because the irregular surface ground produces the strong reflection effect on the total EWs, which reduce the normalized field intensity. The single extreme point in cyan triangle for flat ground and in red circle for random terrain on maximum field values curves can be seen clearly. Then, the location of single target on the flat ground can be fast determined at 5002.7 m distance and 14.3 m height. The results show that the horizontal position error is 2.7 m and the vertical position error is 0.7 m for flat ground. Furthermore, it is shown that the same target on the random terrain is located at 5098.3 m distance and 10.3 m height. While the horizontal position error is 98.3 m, and the vertical position error is 4.7 m for random ground. It is obvious that the precision of location on flat ground is greater than on random terrain.

[figure(s) omitted; refer to PDF]

Third, the double targets are set on the flat ground or random piecewise linear terrain, and the simulation results based on the positive PE and IDPE are shown in Figures 8–11. Supposing the first rectangle object is located at a distance of 3 km, which are 100 meters wide and 10 meters high. The second rectangle object located at a distance of 7 km, which are 100 meters wide and 30 meters high. Other simulation conditions are the same as those in Figure 4 to 6.

[figure(s) omitted; refer to PDF]

In Figure 8, the forward EWs propagate along with the positive x direction. It is obvious that there are many interference stripes scattered in the whole calculation region. The interference effects are generated from the interaction between the direct EWs and reflected EWs from ground and targets.

In Figure 9, the backward propagating EWs distribution is shown. When the front surface of target receives constructive interference waves, it produces the relative strong backward EWs. Conversely, when the front surface of target receives destructive interference waves, it produces the relative weak backward EWs. Thus, the maximum of backward EWs appear at different positions on the front surface of targets.

[figure(s) omitted; refer to PDF]

In Figure 10, the inverse diffraction EWs distribution based on the IDPE is shown. The results show that IPDE method can accurately model the scattering effect from the irregular terrain and targets, and the position of maximum filed intensity is that of the targets. It is obvious that two vertical axis of symmetry appeared, which means the two targets are located.

[figure(s) omitted; refer to PDF]

In Figure 11, the maximum filed intensity curves of inverse diffraction EWs at every horizontal step are shown. The total computing time is 296.9 seconds. The location of the first target on the flat terrain can be fast determined in green triangle on the maximum value curves, which is at 2993.7 m distance and 9 m height. While the first target is located in pink circle on maximum values curves, which is at 3021 m distance and 8 m height on the irregular terrain. The results show that the horizontal position error is 6.3 m, and the vertical position error is 1 m for the first target on flat ground. While the horizontal position error is 21 m and the vertical position error is 2 m for the first target on random fractal terrain. Then, the location of the second target on the flat terrain can also be fast determined in red triangle on maximum values curves, which is at 7009 m distance and 28.7 m height. While the second target is located in a cyan circle on a maximum values curves, which is at 7112.3 m distance and 26 m height on the irregular terrain. The results show that the horizontal position error is 9 m, and the vertical position error is 1.3 m for the second target on flat ground. While the horizontal position error is 112.3 m and the vertical position error is 4 m for the second target on random terrain. It is found that the farther the targets away from the transmitter, the less the backward EWs from the surface targets, and the lower the location precision is.

[figure(s) omitted; refer to PDF]

Furthermore, the analysis of the influence of the sampling interval of backward field intensity at source end on the position accuracy is shown in Figure 12. The EWs frequency is set to be 300 MHz, 1.5 GHz, and 2.5 GHz, respectively. The single target is located at a distance of 5 km with the height of 10 m. The sampling interval of backward field intensity from the receiver at initial position is set to be 0.2 m, 0.5 m, 1 m, 2 m, 3 m, 4 m, 5 m, and 6 m. The position error curves show that the horizontal position errors is increased withing the increasing of sampling interval of the received field intensity. When the sampling interval is less than 0.5, the horizontal position error is less than 15 m at different frequencies. Furthermore, it is found that as the radio frequency decreases, the effect of sampling interval of received field intensity on the position accuracy decreases. Based on a large amount of statistical data, we have also observed that the range of vertical position error is relatively large. This is because the convergence point of inverse diffraction field located in the most radioactive areas of constructive interference waves, which affected by many factors, such as frequency, beam width, inclination angle, and terrain. Thus, the error of height identification for multiple targets is comparatively large.

[figure(s) omitted; refer to PDF]

Finally, the influence of piecewise linear terrain slopes on the false alarm probability of objects localization is also analyzed in Table 1. Supposing the discontinuous terrain slope S is set to be from 0° to 15°. It is found that with the discontinuous terrain slope increasing, the mean value and standard deviation of the terrain heights increases significantly. The results show that the higher the discontinuous terrain slope, the higher the false alarm probability. Furthermore, the false alarm probability increases with increasing the numbers of location objects.

Table 1

False alarm probability for piecewise linear terrain.

Parameters of piecewise linear terrain			False alarm probability
Discontinuous terrain slope (°)	Mean of terrain heights (m)	STD of the terrain heights (m)	Single object	Multiple objects
0∼5	0∼1.8	0∼1.1	<0.01%	<0.01%
5∼10	0.9∼3.5	0.4∼1.6	<1%	<2%
10∼15	1.4∼7.3	0.7∼2.7	<3%	<5%

4. Conclusions

This article presented the inverse diffraction parabolic equation (IDPE) based on FDM. It is first applied to study multiple nonradiation targets orientation technology for large-scale irregular terrain environment. The piecewise linear function is used to model the irregular terrain. According to the sampled profile data of field intensity at source position, the distribution of reflection and refraction waves for whole space can be calculated efficiently by IDPE. Then, the distance and heights of the multiple targets can be determined using the local extreme search method rapidly. The results show that horizontal position errors is increased withing the increasing of the radio frequency and sampling interval of field intensity. Furthermore, with the increasing of the discontinuous terrain slope and the numbers of location objects, the location accuracy decreases. In general, these studies provide a novel and fast orientation method for multiple nonradiation targets in large-scale complicated terrain environment.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant no. 61901532) and the National Defense Basic Research Program (JCKY2020210C614240301).

References

[1] F. Hacivelioglu, L. Sevgi, P. Y. Ufimtsev, "On the modified theory of physical optics," IEEE Transactions on Antennas and Propagation, vol. 61 no. 12, pp. 6115-6119, DOI: 10.1109/tap.2013.2280042, 2013.

[2] M. M. Taygur, T. F. Eibert, "A ray-tracing algorithm based on the computation of (exact) ray paths with bidirectional ray-tracing," IEEE Transactions on Antennas and Propagation, vol. 68 no. 8, pp. 6277-6286, DOI: 10.1109/tap.2020.2983775, 2020.

[3] J. Li, W. Meng, L. X. Guo, Y. J. Xi, "Fast solution of scattering from moving target by dynamic-octree-based SBR algorithm," IEEE Antennas and Wireless Propagation Letters, vol. 22 no. 2, pp. 362-366, DOI: 10.1109/lawp.2022.3212416, 2023.

[4] S. Burns, F. Gasdia, J. J. Simpson, R. A. Marshall, "3-D FDTD modeling of long-distance VLF propagation in the earth-ionosphere waveguide," IEEE Transactions on Antennas and Propagation, vol. 69 no. 11, pp. 7743-7752, DOI: 10.1109/tap.2021.3070621, 2021.

[5] Z. El Ahdab, F. Akleman, "Radiowave propagation analysis with a bidirectional 3-D vector parabolic equation method," IEEE Transactions on Antennas and Propagation, vol. 65 no. 4, pp. 1958-1966, DOI: 10.1109/tap.2017.2670321, 2017.

[6] M. S Mikhailov, E. S Malevich, V. A Permyakov, "Modeling of radio–wave propagation in forest by the method of parabolic equation," International Journal of Engineering and Technology, vol. 7 no. 2.23, pp. 111-113, DOI: 10.14419/ijet.v7i2.23.11895, 2018.

[7] L. Huang, X. P. Wu, Z. Li, Y. Lu, M. Wang, Y. Long, "A parallel FDTD/ADI-PE method for ultralarge-scale propagation modeling of ILS signal analysis," IEEE Antennas and Wireless Propagation Letters, vol. 19 no. 12, pp. 2245-2249, DOI: 10.1109/lawp.2020.3029193, 2020.

[8] H. Qin, X. Q. Zhang, "Efficient radio wave propagation modeling in tunnels with a sparse fourier transform-based split-step parabolic equation method," IEEE Antennas and Wireless Propagation Letters, vol. 22 no. 10, pp. 2442-2446, DOI: 10.1109/lawp.2023.3290645, 2023.

[9] J. Li, C. Zhang, K. Sun, C. Zhang, C. Zhou, X. Yan, J. Li, "Application of a 3D parabolic equation model with higher-order operator approximation," Proceedings of the IEEE 6th International Conference on Electronic Information and Communication Technology (ICEICT), .

[10] G. Apaydin, L. Sevgi, "A Novel split-step parabolic-equation package for surface-wave propagation prediction along multiple mixed irregular-terrain paths," IEEE Antennas and Propagation Magazine, vol. 52 no. 4, pp. 90-97, DOI: 10.1109/map.2010.5638238, 2010.

[11] Z. He, H. Zeng, R. S. Chen, "Two-way propagation modeling of expressway with vehicles by using the 3-D ADI-PE method," IEEE Transactions on Antennas and Propagation, vol. 66 no. 4, pp. 2156-2160, DOI: 10.1109/tap.2018.2805018, 2018.

[12] T. A. Spencer, R. A. Walker, R. M. Hawkes, "GNSS Interference localisation method employing inverse diffraction integration with parabolic wave equation propagation," Proceedings of the 17th International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS 2004), pp. 87-98, .

[13] T. A. Spencer, R. A. Walker, R. M. Hawkes, "Inverse diffraction parabolic wave equation localisation system (IDPELS)," Journal of Global Positioning Systems, vol. 4, pp. 245-257, DOI: 10.5081/jgps.4.1.245, 2005.

[14] G. Apaydin, L. Sevgi, Radio Wave Propagation and Parabolic Equation Modeling, 2017.

[15] M. Levy, Parabolic Equation Methods for Electromagnetic Wave Propagation, 2000.

[16] O. Ozgun, "Recursive two-way parabolic equation approach for modeling terrain effects in tropospheric propagation," IEEE Transactions on Antennas and Propagation, vol. 57 no. 9, pp. 2706-2714, DOI: 10.1109/tap.2009.2027166, 2009.

[17] International Radio Consultative Committee, "Recommendations and Reports of the CCIR," Propagation In Non-Ionized Media, vol. 5, 1982.