Parabolic Equation Accumulated Split Error and Its Correction Method for Tropospheric Radio Propagation

An increasingly popular approach of modeling electromagnetic wave propagation in the troposphere is to use the slip-step Fourier transform (SSFT) numerical algorithm to solve a parabolic equation (PE). However, available SSFT does not provide perfect solutions by introducing accumulated split error (ASE) for its step-by-step iteration process with a fixed propagation range step. The common analysis on one step split error is only a simple discussion on the split error of PE. Therefore, the main motivation of this study is to provide improvement on the accuracy of the PE by proposing an error correction method for the first time in the literature for PE. This method is on the basis of defining and deriving the accumulated split error. Its performance is then demonstrated with measured radar sea clutter data from Huanghai Sea, China.


Introduction
Most radio systems have to operate in the lower tropospheric atmosphere.Atmosphere absorbing, scattering, reflection, and refraction could affect radio propagation characteristics.In particular, the occurrence of tropospheric duct can make radio wave beyond-line-of-sight propagation [1], extending communication range and detection distance of radar system.This phenomenon also may cause radio holes and introduce new interfaces, which fundamentally can change the performance of a system that is designed to operate under standard atmospheric conditions [1].This nonstandard atmosphere has a higher demand on effectively describing propagation characteristics.
Over time, radio propagation modeling methods in tropospheric environment mainly focus on ray tracing (RT) [2,3] and parabolic equation (PE) [4,5] methods.RT generally neglects the influence of frequency and fails to describe the distribution of field strength at arbitrary location because of restricting to ray numbers.PE is a full-forward-wave approximation calculation method which is derived from wave equation using paraxial approximation [5].It can deal with electromagnetic propagation problems in arbitrary irregular terrain and inhomogeneous atmosphere with the help of computer numeric calculation.Therefore, PE has been the most popular approach for tropospheric electromagnetic propagation [1][2][3][4][5][6][7][8][9][10].According to propagation angle applied to the Helmholtz wave equation to obtain PE, PE can be divided into wide-angle PE (WPE) and narrowangel PE (or standard PE (SPE)).
At present, there are mainly three algorithms to solve PE: split-step Fourier transform (SSFT) [1], finite difference schemes [11], and finite element method [12], where SSFT has the advantage of higher numerical stability and thus is widely used to solve PE for tropospheric propagation problems.However, PE has approximation error due to paraxial approximation and split error as a result of exponent operator split using SSFT, where approximation error is its inherent drawback.For the past work on the error analysis of PE, approximation error is lower for SPE, and it can be ignored [13].The expression of one-step split error of SPE is described in [14].The two kinds of errors for WPE and SPE are discussed with simple calculation in [4].Some advices are proposed to decrease the split error by studying its relationship with range step in [15].However, an analytical description of the accumulation of split error (called accumulated split error (ASE)) over many range steps has yet to be successfully developed.In summary, these studies are inefficient to degrade split error.For this problem, we conduct the research on ASE and provide its correction method in this paper.In our work, the accumulated split error of SPE is derived for the first time in the literature on account of the definition of split error [4] and the iteration idea of SSFT.The derivation result shows that the ASE is determined by the atmosphere refractivity, the antenna frequency, the propagation elevation angle, the range step, and the field strength with split error.We simulate the sensitivity of ASE resulting from these factors using two statistical methods in a standard and nonstandard atmosphere.With the research conclusions, the ASE can be degraded by taking small range step in a standard atmosphere.However, this solution does not work well for nonstandard atmosphere on some conditions.For solving this problem, an ASE correction method is provided.This method is involved in complex calculation, including higher order accumulation, multiplication, and power operation.In addition, some supplements are provided for this method to analyze the complexity and efficiency.Finally, this method is verified by measured radar sea clutter from Huanghai Sea, China.
The main contributions of this paper are fourfold.Firstly, we derive the accumulated split error of SPE for the first time in the literature.Secondly, the sensitivity analysis of ASE in a standard and nonstandard atmosphere is discussed using two statistical methods.Thirdly, we develop an ASE correction method to improve the accuracy of PE.Lastly, radar sea clutter power calculated by our proposed approach and past PE method is compared with the measured data from Huanghai Sea China, respectively.
The remainder of the paper is organized as follows.Section 2 presents the deviation of the ASE and its iteration model.Section 3 provides simulation results of sensitivity analysis.In Section 4, the ASE correction method is introduced.In addition, experiment research is provided.Lastly, the conclusions are presented in Section 5.

Parabolic Equation Accumulated Split Error
This section presents the derivation of the ASE on the basis of the one split error of PE.
2.1.Parabolic Approximation.In the rectangular coordinate system shown in Figure 1, the SPE for the field u x, z is expressed as [5].
where k is the wave number and k = 2π/λ (λ is the wavelength.), a is the Earth's radius, and n is the refractive index.The term 2z/a accounts for the spherical shape of the Earth.When it is omitted, (1) describes propagation over a flat Earth.
Only considering the forward wave propagation, (2) is obtained after simple transform of (1) where the expressions of two operator terms A and B are with w = n 2 x, z − 1 + 2z/a .The horizontally inhomogeneous refractivity index n is considered here, shown as n x, z .Integrate over x for (2), then the solution is given as where Δx is the range step.
In order to solve this formula, error is incurred by splitting the exponential with the hypothesis that n is a constant.The simplest split is given by [4].
Then the expression that can be solved is provided by Using (6) and other boundary conditions shown in Figure 1, the field can be described in the propagation domain [5].It is the commonly used method.
The term w is varied with height in the realistic environment, then Therefore, the two terms splitting in the exponent is not exact in nonhomogeneous medium.The split error E error is provided as This formula is the one-step split error at the range x assuming that the field u x, z is the same for the right- 2 International Journal of Antennas and Propagation hand side of (8).However, the field strength values of u x, z are different for the right-hand sides of ( 4) and ( 6), except for the initial field u 0, z at the range x = 0. Therefore, ( 8) is not applied for correcting the ASE between ( 4) and ( 6) over many range steps.
2.2.Accumulated Split Error Modeling.We assume that there is a solution meeting (4), then where F 1 = exp Δx A + B , and u j is the field at the jth (0 ≤ j ≤ N x − 1, N x is the horizontal grid points) range step, which gives the result using PE without ASE.Let u j ′ and u j−1 ′ obey ( 6), thus the field can be solved [4,5], which gives the result using PE with ASE. where In the following derivations, the Hadamard product is introduced to describe the product and division of elements in the vector, where they are expressed as "∘" and ". /," respectively.
The ASE R j at the jth range step is shown as where u j ′ and u j contain N z vertical grid points, respectively.
For the convenience of derivation, we define the correction factor With the detailed derivation found in the Appendix, S j is rewritten as Referring to the iteration process of SSFT, we derive the ASE.
Given the initial field u 0 , the field of next range step x = Δx is shown as Using (12), the field u 1 is given by For the range step x = 2Δx, u 2 satisfies the following relationship Then we get And so on, the field of range x = jΔx is provided by where ∑ k p=0 t p = j − k , 0 ≤ t p ≤ j − k , t p is the integer, and C k j represents the combination, which can be written as C k j = k / j k − j using the factorials.We now get For horizontally homogeneous atmosphere, the refractivity index n only varies with z.Then the term S j in (13) is the same at each range and unified as S 0 .Considering these relationships given in ( 14)-( 19), R j is obtained in homogeneous medium We can conclude from above derivation: with the j increasing, the calculation complexity of u j increases; the ASE R j is affected by environment parameters (n j , ∂n j /∂z, and ∂ 2 n j /∂z 2 ), antenna parameters (the frequency f ), calculation range step Δx, propagation elevation angle α, and u j ′ ; the derivation process is involved in higher order accumulation, multiplication, and power operation.In order to effectively reduce the ASE, we study its sensitivity by varying these parameters.

Sensitivity Analysis
On account of the derivation results of ASE, its sensitivity involved by various parameters is addressed using computer 3 International Journal of Antennas and Propagation calculation and statistical analysis methods.Standard and nonstandard atmosphere are examined, where the nonstandard atmosphere takes evaporation duct model as an example.
Since the transmission loss is the most commonly used parameter to describe the radio propagation characteristics, we utilize this quantity instead of the field to discuss the sensitivity of ASE.The transmission loss calculated from u j ′ and u j is expressed as L depart and L join , respectively.
Simulations are given to compare the differences between L depart and L join .Antenna parameters are shown in Table 1.The distributions of L depart and L join in a standard atmosphere are described in Figure 2 with the condition of f = 10 GHz, and Δx =500 m (calculation domain: 60 km × 75 m).Besides, the distributions of L depart and L join in an evaporation duct are described in Figure 3 with the same calculation condition with Figure 2, where the simplified Paulus-Jeske (PJ) evaporation duct model [16] is taken with duct height of 20 m.
Comparing Figure 2(a) to Figure 2(b), the effect of ASE on transmission loss is more prominent at the range greater than 20 km.It is because that the cumulative effect of split error increases with the range increasing in a standard atmosphere.Comparing Figure 3(a) to Figure 3(b), the difference of transmission loss is greater inside the duct, especially nearby the ocean.The reason is that the greater gradient varieties of refractivity index along the height direction in (13) lead to greater ASE.
In order to further discuss the effect of environmental parameters, Δx, and f on the ASE, two statistical methods are utilized: one dissimilarity method and one similarity method.
The first method defines the mean value ΔL of relative percentage difference which is calculated from the absolute difference between L join x, z and L depart x, z , and dividing by L depart x, z .Therefore, the mean characteristics of transmission loss difference in the given propagation domain are described.
The second method applies correlation coefficient R between L depart and L join ; thus, the bias from the mean value is considered.

23
where μ join and μ depart are the mean of L depart and L join , respectively.
With the two methods, the simulation is conducted with the frequencies of L band, S band, C band, and X band, which are commonly used for microwave electromagnetic propagation.In addition, the range steps are taken from 50 m to 500 m.

Sensitivity Analysis in a Standard
Atmosphere.We utilize the above-mentioned two methods to analyze the sensitivity in a standard atmosphere.
Figure 4 presents the distributions of ΔL and R varying with f and Δx in a standard atmosphere.In the four frequency bands, ΔL increases with range step, but R decreases.As noticed, the range step has weak effect on the frequency of L band, making ΔL smaller than 1% and R greater than 0.99.Thus, the ASE can be neglected for PE in this band.In order to control ΔL ≤ 1% and R > 0.99, the range step of smaller than 350 m should be taken for the frequency of S band.We should use the range step smaller than 150 m for the frequency of C band and smaller than 100 m for the frequency of X band.

Sensitivity Analysis in an Evaporation
Duct.Two statistical methods are applied for analyzing the sensitivity in an evaporation duct.
Figure 5 provides the distributions of ΔL and R varying with f and Δx in an evaporation duct, which is obviously different from Figure 4. Comparing Figure 5 to Figure 4, the difference between the maximum and the minimum of ΔL in Figure 5(a) is smaller than that in Figure 4(a), but the difference range of R in Figure 5(b) is greater than that in Figure 4(b).Therefore, the effect of ASE on transmission loss in an evaporation duct is relatively greater to that in a standard atmosphere.In order to control ΔL ≤ 1% and R>0.99, we utilize range step smaller than 300 m for the frequency of L band, smaller than 200 m for the frequency of S band, smaller than 150 m for the frequency of C band, and smaller than 150 m for the frequency of X band.
In contrast, examining the distribution of the overall metrics (ΔL and R) for the transmission loss inside the duct is shown in Figure 6.The comparisons of Figure 6, Figure 5, and Figure 4 show that L depart in agreement with L join is the least close inside evaporation duct, which leads to the greater difference between the maximum and the minimum of R than that in an entire domain.This is because the first derivative of refractivity index dramatically varieties with heights and the second derivative is nonzero.Therefore, the ASE cannot be decreased by simply taking the smaller range step in the high frequency, such as X band.

Accumulated Split Error Correction Method
On account of the above-mentioned simulations and analysis, the PE-accumulated split error correction method is proposed.This method is verified by transhorizon tropospheric duct test data.We propose the correction method of ASE contained in u j ′ on the basis of its definition and derivation in Section 3. The method can decrease the distribution domain of greater error values, thus improve the calculation accuracy of PE for electromagnetic propagation in the tropospheric environment.

Accumulated Split Error Correction
Method.The steps applied to correct the error are provided as follows, and its diagram is shown in Figure 7. (2) Build propagation environment using refractivity profile.
Using parameters from (1) and a given range step Δx, ASE correction factor S j (j ∈ 1, N x ) is calculated by (13) in the environment of (2).
On the basis of (3) and ( 4), the ASE is corrected using ( 14)-( 19).( 6) Output the field without the ASE.Using (5) yields the field u j j ∈ 1, N x .As noticed in Figure 7, the SSFT to solve PE is involved in complex techniques, including the absorbing boundary, surface boundary modeling, and split-step calculation.Besides, the interaction process from ( 14) to (19) shows that higher order accumulation, multiplication, and power operation are introduced.Therefore, this correction method increases the calculation complexity to improve the accuracy of PE.Taking the standard atmosphere for an example, we provide the complexity and efficiency analysis.The simulation results are shown in Figure 8.These results are based on the two aspects as follows: (1) The statistical relationship between calculation step range and execution time is analyzed for PE with ASE method and PE without ASE method, respectively.These results are simulated using the computer CPU of i7-4270HQ and shown in Figure 8(a).
(2) Given the upper threshold of error ΔL, the statistical relationship between frequency and execution time is analyzed and given in Figure 8(b).
As noticed in Figure 8(a), the execution time decreases with the range step increasing for the two methods.However, the execution time of PE without ASE method (red color) is much greater than that of PE with ASE method (blue color) on the condition of the same range step, especially for smaller   The radar antenna parameters are given in Table 2.
Since there is a clear weather for this measured data, smooth sea surface and range-independent refractivity profile are considered.By filtering the weather echo, ground echo, noise, and so on of sea clutter strength X in Figure 9(a), we can obtain the sea clutter power P obs at range r using where C takes into account radar parameters and so on.Applying the measured modified refractivity M into the radar equation [17], the power can be computed by where L M, r is the transmission loss at range r calculated by PE with ASE and without ASE method; σ 0 is the radar cross section of the sea surface.In order to decrease the effect of range step on calculation results, we take Δx = 150 m.The field u j ′ is computed as the first step result.Then the field u j can be obtained.On the basis of the field, the theoretical power with ASE P c ′ M, r and power without ASE P c M, r is shown in Figure 10, respectively, where the mean of error between P c ′ and P obs is 12.3 dB, but that between P c and P obs is 9.4 dB.Therefore, the parabolic equation ASE correction method can effectively improve the calculation accuracy.

Conclusion
In this paper, we derived the ASE due to exponent split using SSFT numerical algorithm to solve PE.In addition, sensitivity analysis of ASE in standard and nonstandard atmosphere is addressed using two statistical methods.On account of above studies, we propose an ASE correction method for PE.We also verify that PE without ASE method has a greater accuracy by using the measured data from Huanghai, China.Although the ASE of SPE is only derived and simulated, the study idea we developed is also applied to WPE.This is also our future work.

Figure 1 :
Figure 1: The calculation domain of PE.

Figure 2 :Figure 3 :Figure 4 :
Figure 2: The transmission loss distribution in a standard atmosphere.(a) The distribution of L depart .(b) The distribution of L join .

Figure 5 :
Figure 5: The distribution of ΔL and R for the 20 m evaporation duct.(a) ΔL and (b) R.

Figure 6 :
Figure 6: The distribution of ΔL and R for the region of inside duct (z < 20 m).(a) ΔL and (b) R.

6
International Journal of Antennas and Propagation range step.This result agrees with foregoing analysis that the calculation complexity increases by involving in complex correction factor.On the basis of accepted error threshold and execution time in Figure 8(b), the range step can be determined for a given frequency by referring to Figure 8(a).Therefore, Figure 8(b) gives a reference for realistic application.4.2.Experimental Example.The main purpose of this section is to verify the performance of the ASE correction method of PE.The sea clutter data (Figure 9(a)) and corresponding data of a refractivity profile (Figure 9(b)) measured in Huanghai, China, August 1, 2014 are considered here.As noticed in Figure 9(a), the north direction is the azimuth 0 o , and the sea clutter at 58 5 o azimuth is used as an example to analyze the performance of parabolic equation ASE correction method.

Figure 8 :
Figure 8: The complexity and efficiency analysis of PE.(a) The statistical relationship between calculation step range and execution time.(b) The statistical relationship between frequency and execution time.