Modeling long range Loran-C signal propagation numerically is challenging because of its extremely high computational cost. Other analytical/semi- analytical methods are not accurate enough due to unavoidable approximations. In this paper we present a solution using the compute unified device architecture (CUDA) based parallel finite-difference time-domain (FDTD) method with an adaptive moving window (AMW) technique. The moving velocity of the window is varying with the wave speed adaptively. To achieve the AMW technique, the original Loran-C signal is truncated first. A new method used for extracting the electric field amplitude and phase is proposed. The electric field amplitude and phase data of each mesh in the computational space are synchronously obtained from spatial domain as the FDTD updates, without additional storage cost and post processing in time domain. With all these efforts a 400 km propagation path was simulated successfully within 22 minutes. Measurement results between Jinxian and Shangrao in Jiangxi Province, China were taken to validate the numerical approach.
Index Terms'Ground-wave propagation, Loran-C, parallel FDTD, compute unified device architecture, adaptive moving window.
Ground wave propagation along the earth's surface has been widely studied using both analytical solution and numerical approaches such as finite element method (FEM), method of moment (MoM), and finite-difference time-domain (FDTD) method [1-5]. The analytical/semi-analytical approaches are often very fast with little requirement on computational resources. However, they suffer from poor accuracy especially for rapidly changing terrains due to model simplifications. On the other hand, numerical methods model the terrain details accurately but their requirement on computational resources becomes extremely high when long range propagation is of concern.
Global Positioning System (GPS) by US goverment and GLONASS (Globalnaya Navigatsionnay Sputnikovaya Sistema) by the Russian Federations are both satellite based positioning systems. As GNSS (GPS and GLONASS) are inherently vulnerable to interference, disruption, jamming, and spoofing, many government and civilian organizations around the word are devoting to the upgrading research of Long Range Navigation System (Loran-C) [6-10]. Loran-C was a ground-based navigation system operated by the U.S. Coast Guard. The prediction error of wave propagation characteristics is one of the major reasons why the systems suffer from poor positioning/timing accuracy. The FDTD method has proven to be accurate in predicting the electric field strength and secondary time delay of Loran-C ground wave propagation [11-13]. In , the combined IE-FDTD algorithm is used to solve the Loran-C ground wave propagation along long-distance paths with partial complex area. In this work, the original Loran-C signal is truncated first. Then the FDTD algorithm is implemented using the compute unified device architecture (CUDA) based parallel programming  together with an adaptive moving window (AMW) technique. With this approach the simulation of long range complex path, which was never possible using the conventional FDTD method, can be managed in a reasonable short time.
II. PARALLEL AWM FDTD BASED ON CUDA
Fig. 1. Truncated Loran-C current signal
A Loran-C signal is a 100 kHz Gaussian modulated pulse. In conventional Loran-C operations, the negative peak value and positive zero-crossing time of the third carrier cycle of the Loran-C signal are measured as the field strength and the arrival time to prevent contamination by the sky wave. By the nature of the Loran-C tracking, a windowed scheme can be employed in the FDTD simulation. For the duration of a Loran-C signal is about 500??s, if the computation window contains the whole signal, the cell number will be too large to simulate with high precision. As the signal propagation characteristics are not concerned after the third cycles, the original Loran-C is truncated first. Fig. 1 shows the truncated Loran-C current signal, given by
where and The peak value of the current signal envelope in (1) can be obtained by
Fig. 2(a) demonstrates a two-dimensional (2D) rotational symmetric AMW FDTD simulation model. More details are shown in Fig. 2(b). The process of the AMW FDTD of this model is as follows.
1) Initial stage: a dipole source is located at the coordinate axis one grid higher than the earth surface. The window stands still at this stage and the truncated Loran-C signal propagates along the ?? direction in it. The left boundary is an axisymmetric boundary and lossy ground termination is used at the bottom. Ten cells of perfectly matched layer (PML) are set on top. There is no absorbing boundary on the right, as the signal wavefront will not reach the right boundary of the window. The windows of two different time steps (n1 and n2, n1 < n2) are shown in Fig. 2(b).
2) Adaptive moving stage: an adaptive moving window covering the first three cycles of the Loran-C carrier is defined. The window moves synchronously with the propagating signal to have the signal centered in it. The following principle is used to determine whether the signal is in the window center. Assuming the ??-direction mesh number is N??. Analyze the spatial-domain waveform constituted by the electric fields of all meshes one grid higher than the earth surface from right side to the left of the window. If it happens that the zero-crossing point of the third carrier cycle is in the center of the analyzed meshes (??-direction mesh: N??/2), the signal is considered to be in the center of the computational region and the window keeps still. Otherwise, it needs to move forward a grid distance. To absorb the reflected waves, 10 cells of PML are set on the left. The windows of two different time steps (n3 and n4, n3 < n4) are shown in Fig. 2(b).
3) End of the stage: when the signal is about to reach the boundary of the entire simulation area, we stop moving the window, but continue updating the electromagnetic field until the signal pulse completely propagates out of the window. At this stage, a PML boundary is set on the right. The windows of two different time steps (n5 and n6, n5 < n6) are shown in Fig. 2(b).
It should be noted that, for the moving velocity of the window is varying with the wave speed adaptively, the technique used is not exactly the same as that in [16-19], in which the moving velocity is fixed. This technique used is more suitable for solving the pulse wave propagation problems in complex mediums as the wave velocity is constantly changing. The feature makes this method can be applied to other more occasions.
With the AMW technique, a new method used for extracting the Loran-C electric field amplitude and phase is proposed. In order to get the trend of the electric field amplitude and phase changes along the propagation path, in the traditional simulations, we need to store the time-domain data of multiple sampling points. For each sampling point, the time-domain information at least including the time steps of three carrier cycles should be recorded. Then the negative peak value and zero-crossing time of the third carrier cycle can be obtained by post processing in time domain. Through this way, a mass of store resources will be cost especially when we need to know the electric field amplitude and phase of all meshes on the ground in the entire simulation space.
In this paper, the amplitude and phase information are not obtained from time-domain waveform, but from spatial waveform automatically. The electric fields of all meshes one grid higher than the earth surface have been analyzed before to determine whether the window should move forward or not, as shown in Fig 2(b). At each time step, record the time and the third zero-crossing position (point A) of the spatial pulse waveform. The time recorded is exactly the zero-crossing time of the third carrier cycle in time domain at position A, from which the phase information of the radio wave can be obtained. The negative peak value position (point B) nearest to the zero-crossing location and the peak value are recorded too. The peak value is used to approximate the amplitude information needed to be extracted in time domain of position B. By real-time processing, the amplitude and phase information of each mesh on the earth surface will be obtained without additional storage space cost.
Fig.2. 2D AMW FDTD simulation model
In fact, with the AMW technique the analyzed data can be reduced further, for once the window starts to move steadily, the zero-crossing point of the third carrier cycle of the space waveform on the earth surface will always keep nearby the center of the computational window (the position error is no more than 1 grid). Thus, in Fig. 2(b), only the mesh data between the two red points (within about one-quarter wavelength) need to be processed.
The AMW FDTD then is accelerated by using a large number of parallel stream processing on the graphics processing unit (GPU) device. The approach and use of CUDA has been reported in . In this work, the field update of the AMW FDTD is implemented in the GPU. The judgment of adaptive moving and the field extraction at the observation points are performed in the CPU. The flow charts of the traditional algorithm and the new algorithm are shown in Fig.3.
(a) Flow chart of the traditional algorithm
(b) Flow chart of the new algorithm
Fig.3. Flow charts of the traditional algorithm and the new algorithm
III. NUMERICAL VERIFICATION
A. Flat earth model
We first use the flat earth model to validate the parallel AMW-FDTD algorithm by comparing the results to the flat-earth formula [13, 20-21] and full FDTD solution.
The ground parameters are taken as ??r=13 and ??=3??10-3S/m. The total propagation distance is 100 km. In order to limit numerical dispersion, the FDTD cell size is set to be 1/320 of the wavelength, i.e., 9.325m; the time step is 15.625 ns. The full FDTD computational domain is meshed by 11,000??600 cells; while the AMW-FDTD model contains 2,600??600 cells. Both FDTD calculations run 25,000 steps. The total moving number is set as 8,400. In the AMW-FDTD model, the spatial field and topography distribution of six time steps (n=3,000, 4,000, 6,000, 10,000, 22,000, 23,000) are extracted as shown in Fig. 4-6. The steps n=3,000 and 4,000 are in the initial stage, n=6,000 and 10,000 are in the adaptive moving stage, while n=22,000 and 23,000 are in the end of the stage. The moving number (denotes as Mn), third zero-crossing position A (denotes as PA), negative peak value position B (denotes as PB), range of position A (denotes as dA), range of position B (denotes as dB), time of signal arriving (denotes as Ttoa), and the electric field of position B (denotes as ez) of each step spatial pulse waveform are recorded in the table 1, where
The parameter A expressed in (1) is set with the following equation:
where is the radiation power, is the wavelength, and is the length of the dipole source. In the following examples is set to 1 kW.
(a) n=3,000, Mn =0 (b) n=4,000, Mn =0
Fig.4. Spatial distributions of the electric field and topography as time step n=3,000 and 4,000 (initial stage, window stands still)
(a) n=6,000, Mn =799 (b) n=10,000, Mn =2793
Fig.5. Spatial distributions of the electric field and topography as time step n=6000 and 10000 (adaptive moving stage, window moves synchronously with the signal)
(a) n=22,000, Mn =8400 (b) n=23,000, Mn =8400
Fig.6. Spatial distributions of the electric field and topography as time step n=22,000 and 23,000 (end of the stage, window stands still)
DATA OF SIX TIME STEPS FOR THE FLAT EARTH
Time step Moving number Position A Rang A (km) Ttoa (??s) Position B Range B
(km) ez(??V/m) Moving window
3,000 0 604 5.6625 46.875 676 6.3375 -24219 No
4,000 0 1103 10.3406 62.500 1177 11.0344 -13762 No
6,000 799 1300 19.6781 93.750 1375 20.3813 -7323 No
10,000 2793 1300 38.3719 156.250 1376 39.0844 -3724 No
22,000 8400 1680 94.5000 343.750 1756 95.2125 -1439 No
23,000 8400 2180 99.1875 359.375 2255 99.8906 -1365 No
The results of the AMW FDTD are compared with that of the full FDTD as the time step n=10,000. In the full FDTD, standard Loran-C signal in (1) is used. The time-domain waveforms of the corresponding meshes at position A and B are extracted and shown in Fig. 7. The ??-direction mesh numbers are 4093 (PA+Mn) and 4169 (PB+Mn). The electric fields at the time step 10,000 are markd.
(a) ??-direction mesh number is 4093 (b) ??-direction mesh number is 4169
Fig.7. Time-domain waveforms of the electric field as ??-direction mesh numbers are 4093 and 4169
From Fig.7 (a) it can be seen that the positive zero-crossing point of the third carrier cycle is just at the time step 10,000. From Fig.7 (b) it can be seen that when the time step is 10,000, the field value is very close to the negative peak value of the third carrier cycle. The results of the AMW FDTD and full FDTD agree well. It proves that the AMW FDTD method is effective. Furthermore, the AMW FDTD in this paper needs not to store the time-domain data as shown in Fig.7.
The vertical electric field amplitude |Ez| and the secondary time delay tw along the whole propagation path are calculated and compared in Fig. 8. The theoretical results (frequency=100 kHz) with the flat-earth formula are given too. In the two FDTD models
where refers to the time of signal arriving as the propagation path is perfectly conductor.
(a) result comparisons
(b) result comparisons
Fig.8. Comparisons of the FDTD and analytical results over the flat earth surface
As seen in Fig. 8, both FDTD results agree very well the flat-earth formula solution. The maximum error of the electric field amplitude is found to be less than 0.2 dB (??V/m), while the maximum error of secondary time delay is less than 10 ns.
B. Complex terrain
As a second example we test a complex propagation path containing two mountains and one island, as shown in Fig. 9. The electric parameters are defined as: ??r=13, ??=3??10-3 S/m for ground; ??r=80, ??=5S/m for the sea; ??r=6, ??=2??10-4 S/m for the mountain and island, respectively.
In the AMW-FDTD model, the spatial field and topography distribution of six time steps (n=3,000, 4,000, 6,000, 10,000, 22,000, 23,000) are extracted as shown in Fig. 10-12. The parameters Mn , PA, PB, dA, dB, Ttoa, and ez of each step spatial pulse waveform are recorded in the table 2.
Fig. 9. Topography of the propagation path
(a) n=3,000, Mn =0 (b) n=4,000, Mn =0
Fig.10. Spatial distributions of the electric field and topography as time step n=3,000 and 4,000
(a) n=6,000, Mn =797 (b) n=10,000, Mn =2750
Fig.11. Spatial distributions of the electric field and topography as time step n=6,000 and 10,000
(a) n=22,000, Mn =8400 (b) n=23,000, Mn =8400
Fig.12. Spatial distributions of the electric field and topography as time step n=22,000 and 23,000
DATA OF SIX TIME STEPS FOR THE IRREGULAR EARTH
Time step Moving number Position A Rang A (km) te (??s) Position B Range B
(km) ez(??v/m) Moving window
3,000 0 604 5.6625 46.875 676 6.3375 -24219 No
4,000 0 1103 10.3406 62.500 1177 11.0344 -13761 No
6,000 797 1301 19.6687 93.750 1372 20.3344 -7284 Yes
10,000 2750 1300 37.9688 156.250 1376 38.6812 -2494 No
22,000 8400 1646 94.1813 343.750 1724 94.9125 -820 No
23,000 8400 2150 98.9063 359.375 2228 99.6375 -800 No
The results of the AMW FDTD are compared with that of the full FDTD as the time step n=10,000. In the full FDTD, the time-domain waveforms of the corresponding meshes are extracted and shown in Fig. 13. The ??-direction mesh numbers are 4050 (PA+Mn) and 4126 (PB+Mn).
(a) ??-direction mesh number is 4050 (b) ??-direction mesh number is 4126
Fig.13. Time-domain waveforms of the electric field as ??-direction mesh numbers are 4050 and 4126
From Fig.13 (a) it can be seen that the positive zero-crossing point of the third carrier cycle is just in the time step 10,000. From Fig.13 (b) it can be seen that when the time step is 10,000, the field value is very close to the negative peak value of the third carrier cycle. The results of the AMW FDTD and full FDTD agree well. It proves that the AMW FDTD method is effective.
The vertical electric field amplitude |Ez| and the secondary time delay tw along the whole propagation path are calculated and compared in Fig. 14. The theoretical results (frequency=100kHz) with the integral-equation formula [13, 20-21] are given too.
(a) result comparisons
(b) result comparisons
Fig. 14. Comparisons of the FDTD and analytical results over a complex propagation path
From Fig. 13 and 14 it can be seen that in the first mountain and island areas (20-40km and 70-90km), the results of the analytical, Full FDTD and AMW FDTD methods agree very well as the terrains are not very steep. Only in the second mountain area (50km nearby) , where the terrain is very steep, the two kinds of FDTD results are still in good agreement, while the integral-equation results have big difference from the others. In the derivation of the integral-equation method, one of the prerequisites is that the terrain shouldn't be too steep [13, 20-21]. It is concluded that the errors between the FDTD and theoretical results are introduced by the limitation of the integral equation method.
The computational time and memory requirement of the parallel AMW-FDTD and full FDTD method are compared in table 3. The simulation was performed on an Inter (R) Core (TM) i5 CPU 750 and the graphics card was a high-end GeForce GTX 260 with 216 stream processors.
RUN TIME AND MEMORY REQUIREMENT COMPARISON
Method Computation time (s) Memory (Byte)
Full FDTD 107224 CPU: 470M
Parallel AMW FDTD 339 GPU: 112M
C. Long distance verification
To examine the AMW FDTD algorithm over the complex propagation path of long range, the Loran-C signals propagating over an irregular path within 400 km are considered. Three cases are simulated. In the first case, the Loran-C signal expressed in (2) is truncated as , and the MW is meshed by =2,200??1,000. In the second case, the Loran-C signal is truncated as , while the MW is meshed by =4,800??1,000. In the last case, the Loran-C signal is truncated as , while the MW is meshed by =4,800??1,000. The results with larger MW size are thought to contain more scattered/reflected waves. The topography data of the path are shown in Fig. 15. From it we can see the path is extremely complex with higher and steeper mountains than the previous examples. The electrical parameters are set as =13, S/m. The source peak radiated power is set as 1 kW. The AMW FDTD algorithm is run for 85,000 steps. The time-domain waveform comparisons of the observation points at different distances ( 165, 265, and 379.8 km) are shown in Fig. 16 (a)-(c). Other FDTD parameters are set as before.
Fig. 15 Topography data of a long irregular path
Fig. 16 AMW FDTD time-domain result comparisons of three cases as the Loran-C signals propagate over the irregular path
Fig. 16 shows that the AMW FDTD time-domain results in the first three carrier cycles of the three cases all agree very well. This example further verifies the effectiveness of the AMW FDTD algorithm used for long-range irregular propagation paths.
IV. PARALLEL AMW FDTD ALGORITHM FOR REAL LONG-DISTANCE PREDICTIONS
In this section, the parallel AMW-FDTD algorithm based on the CUDA platform is used to predict the Loran-C signal propagating over real paths. The simulated secondary time delays are compared to the measurement results.
(a) Distribution of measurement points
(b) 3D topography of the experiment area
Fig.17. Measurement points distribution and 3D topography of the experimental area
Following the approach in , massive measurements were taken in an area between Jinxian and Shangrao, in Jiangxi Province, China, with the Loran-C transmitter located at Xuancheng station which is 400 km away from Jinxian and 280 km away from Shangrao. Fig. 17 (a) shows the experimental area and the distribution of the measurement points. Fig. 17 (b) shows the three-dimensional (3D) topography of the propagation paths. The area is about 200 km long and 35 km wide.
The parallel AMW FDTD calculated secondary time delay values are compared to the measurement results, as shown in Fig. 18. In the simulation, we take ??r=13, ??=3??10-3 S/m for flat area, ??r=10 and ??=1.4??10-3 S/m for mountain region. Fig. 19 shows the statistics of the errors between the two results.
Fig. 18. Secondary time delay comparison
Fig. 19. Statistics of the secondary time delay errors
Most of the phase errors are less than 400ns. On one hand, the errors between the measured and AMW FDTD results are from the earth electrical parameters (such as dielectric constant and conductivity) used in the simulation for lacking of detailed real parameters. On the other hand, the errors are from the 2D model approximation of the FDTD. The algorithm may introduce large errors, especially when the terrains change greatly along ?? direction. To strictly distinct the errors introduced by the electrical parameters or the terrain, 3D simulation analysis is the only way.
Even though, for a complex path the prediction error of a traditional algorithm is often over 1 microsecond. Obviously, it denotes that the proposed algorithm has an exciting precision. The next step is to expand the algorithm to the three dimensional coordinates for better precision.
Based on the CUDA parallel platform, long-range Loran-C signal propagation problems are solved by the AMW FDTD method. The moving window technique proposed combined with the parallel GPU implementation has made large-scale engineering applications possible. With this algorithm, we can analyze the propagation attenuation, time delay, and dispersion of the LF navigation and timing signals along long complex paths.
Source: Essay UK - http://www.essay.uk.com/free-essays/science/ground-wave-propagation.php