An eLoran Signal Cycle Identiﬁcation Method Based on Joint Time–Frequency Domain

: The eLoran system is an international standardized positioning, navigation, and timing service system, which can complement global navigation satellite systems to cope with navigation and timing warfare. The eLoran receiver measures time-of-arrival (TOA) through cycle identiﬁcation, which is key in determining timing and positioning accuracy. However, noise and skywave interference can cause cycle identiﬁcation errors, resulting in TOA-measurement errors that are integral multiples of 10 µ s. Therefore, this article proposes a cycle identiﬁcation method in the joint time–frequency domain. Based on the spectrum-division method to determine the cycle identiﬁcation range, the time–domain peak-to-peak ratio and waveform matching are used for accurate cycle identiﬁcation. The performance of the method is analyzed via simulation. When the signal-to-noise ratio (SNR) ≥ 0 dB and skywave-to-groundwave ratio (SGR) ≤ 23 dB, the success rate of cycle identiﬁcation is 100%; when SNR ≥ − 13 dB and SGR ≤ 23 dB, the success rate exceeds 75%. To verify its practicability, the method was implemented in the eLoran receiver and tested at three test sites within 1000 km using actual signals emitted by an eLoran system. The results show that the method has a high identiﬁcation probability and can be used in modern eLoran receivers to improve TOA-measurement accuracy.


Introduction
With continuous improvements in informatization, positioning, navigation, and timing (PNT) systems have become indispensable for social development and modern national defense infrastructure, reflecting a country's international status and comprehensive strength [1][2][3]. Although global navigation satellite systems (GNSS) are currently the most widely used type of PNT system, GNSS signals are weak, ineffective at penetration, and susceptible to interference. GNSS systems are thus not suitable for applications in special scenarios, including underground, underwater, and indoors [4][5][6]. Therefore, relying solely on GNSS involves serious risks.
The enhanced long-range navigation (eLoran) system evolved from the Loran-C system. It is an internationally standardized medium and a long-range land-based radio system [7,8]. It can satisfy the requirements of PNT applications in most fields in terms of accuracy, reliability, and continuity. It can also provide a time service better than 100 ns and a position service of 20 m after differential correction [7,9]. When a GNSS signal is unavailable, the eLoran system can be used as a backup, which reduces the risk caused by relying on GNSS; this also allows for coping with timing warfare and navigation Remote Sens. 2022, 14, 250 2 of 20 warfare in modern war [10][11][12]. Therefore, several countries are developing and upgrading their eLoran systems. In August 2019, the US Department of Defense (DoD) released a report titled "Strategy for the DoD Positioning, Navigation, and Timing (PNT) Enterpriseensuring a US military PNT advantage", which stated that Loran-C should be listed as the PNT source, second only to GNSS [13]. China also planned to build three eLoran stations in its western region in 2018 to achieve national coverage of eLoran signals together with the existing eLoran system [14]. The construction of the three stations officially started in 2020. Russia has developed a new land-based radio navigation system called "Scorpion," with the aim of achieving balanced development of satellite-based and land-based PNT systems [15]. South Korea began to develop and deploy its own eLoran system in 2013 due to the deliberate interference of North Korea on its GPS; in 2021, it was announced that the developed eLoran system had commenced operation [16,17]. Although eLoran can be used as to supplement GNSS, the current global eLoran application technology is not sufficiently advanced. Compared with those of GNSS, the accuracy and reliability of the application terminals of eLoran remain inadequate. Therefore, eLoran cannot fully satisfy the requirements of users of PNT services.
The eLoran signal is a long-wave (LW) signal with a frequency of 100 kHz. According to the propagation theory of this frequency band, the eLoran signal is a composite wave consisting of a groundwave signal and a skywave signal [18,19]. The eLoran receiver provides timing and positioning services by measuring the time-of-arrival (TOA) of the groundwave signal in the composite wave. Cycle identification is an important task for the eLoran receiver. It involves correctly identifying the standard zero-crossing (SZC) of the eLoran ground wave signal to measure the TOA [20,21]. Because the SZC point of the eLoran signal has no prominent characteristics, the skywave interference (SWI) signal generated by ionospheric reflection arrives at the receiver a certain cycle after the groundwave signal does, which affects the cycle identification result. The carrier frequency of the eLoran signal is 100 kHz. If the cycle identification is incorrect, it produces errors in the TOA measurement that are integral multiples of 10 µs; this severely affects the positioning and timing accuracy. Moreover, the cycle identification is affected by not only SWI but also noise [19,21]. As the receiving distance increases, the signal amplitude decreases, and the signal-to-noise ratio (SNR) deteriorates, which considerably affects the accuracy of cycle identification.
To solve the problems with eLoran signal cycle identification, in [22], a delayedaddition method was proposed, which uses an eLoran pulse signal to subtract the signal after a certain time delay to achieve phase inversion of the signal at 30 µs. This is achieved by detecting the position of the phase inversion point to realize SZC detection. However, this method exhibits poor anti-noise performance, and SWI causes the result to have multiple phase-reversal points. In [23], the polarity-judgment method was proposed, in which a positive zero crossing is selected, and the signals are accumulated at t − 32.5 µs, t − 22.5 µs, and t − 12.5 µs for a long duration. If the accumulated result exceeds the threshold, then t is the SZC; otherwise, t is moved forward or backward by 10 µs, and the process is repeated. This method is time-consuming, and it is challenging to find a suitable threshold. In [24], a cross-correlation method was proposed, which uses the local standard eLoran signal to correlate with the received signal and calculates the SZC by detecting the maximum value of the correlation result. This method is only applicable to situations where the amplitude of the groundwave is greater than that of the skywave. In [25], the peak-to-peak ratio detection method was proposed, in which the ratio of the signal amplitude at t + 2.5 µs to that at t − 7.5 µs is calculated. If the ratio is close to a standard value, then it is considered as the standard zero-crossing. However, this method is affected by the skywave and noise, which can produce multiple values, resulting in identification errors. In [26], a delay-locked loop (DLL) method was proposed, wherein the lag envelope is subtracted from the advance envelope to obtain the envelope with a peak point of 0, and the SZC is obtained by detecting the zero point of the envelope. This method uses the peak value of the envelope of the ground wave signal; thus, the premise of the Remote Sens. 2022, 14,250 3 of 20 method is to eliminate the skywave. Bian and Last [27] proposed the method of spectrum division to estimate the starting positions of the skywave and groundwave to realize cycle identification. However, this method has a poor estimation resolution under the condition of low SNR. Other researchers have proposed the autoregressive moving average model (ARMA) and multiple signal classification algorithm (MUSIC) methods based on spectrum division to improve the resolution of estimation; however, these methods are not only difficult to implement but are suffer from inaccurate parameter estimation [28][29][30].
Based on the foregoing account of previous research, cycle identification methods can be divided into two types: signal-waveform reconstruction and parameter estimation. In signal-waveform reconstruction methods, cycle identification is realized by constructing the characteristics of SZC. This method can fail under low SNR, and multiple identification results under SWI can be produced; hence, relying solely on this method is hazardous. The parameter-estimation method is mainly based on spectrum-division theory, and the estimation error is high at under SNR. Therefore, a single-cycle identification method is unreliable and has poor anti-noise ability. Hence, this paper proposes a joint timefrequency domain method to realize cycle identification of an eLoran signal. On the basis of improving SNR and repairing the envelope using a digital bandpass filter (BPF) and linear digital averaging (LDA), the range of SZCs is determined via spectral division in the frequency domain, and accurate identification of the SZC is realized through peak-topeak ratio detection and waveform matching. In this study, the performance of the cycle identification method was simulated, and the method was verified via field testing. The methods, findings, and conclusions of this study are presented in the following sections.
The rest of the paper is organized as follows. Section 2 presents the format of the eLoran signal and skywave interference and describes the impact of the skywave interference on cycle identification. The principle and model of the joint cycle identification method in the time-frequency domain are also described. Section 3 describes the verification of the effectiveness of the cycle identification method and the simulation of the performance of the period identification method under different signal-to-noise ratios (SNRs) and skywave interference. The eLoran signals actually received in three places were used to verify the cycle identification method, and the test results were verified. Section 4 summarizes the main conclusions of this study.

Materials and Methods
In this section, we first describe the single-pulse signal format and pulse group signal characteristics broadcast by the eLoran system. Second, along with the propagation channel of the eLoran signal, the influence of SWI on cycle identification is analyzed. Finally, the principle of the cycle identification method based on the joint time-frequency domain is presented, and the key parameters of the algorithm are discussed.

eLoran Signal Format
The working frequency allocated to the eLoran system by the International Telecommunication Union (ITU) is 100 kHz, and the bandwidth is 90-110 kHz. It is required that 99% of the transmitted signal energy be concentrated within the specified bandwidth. At the beginning of the system design, the center frequency, signal-waveform, transmission period, and data-modulation mode of the transmission signal are considered comprehensively, and the existing pulse-phase signal system is formed. The single-pulse signal waveform of the eLoran system is defined by the United States Coast Guard (USCG) as [31].
where t is time (µs), A is a normalization constant related to the peak amplitude, f is the carrier frequency of 100 kHz, τ is the envelope-to-cycle difference (ECD) (µs), and P c is the phase-code parameter, which is 0 for positive phase-coded (PC) pulses and π for negative PC pulses. The eLoran pulse signal is generated according to Equation (1). The time-domain and frequency-domain characteristics of a single eLoran signal are shown in Figure 1. It can be seen that the leading edge of the signal first rises rapidly, reaching a peak at 65 µs, and then falls rapidly.
where is time (μs), A is a normalization constant related to the peak amplitude, is the carrier frequency of 100 kHz, is the envelope-to-cycle difference (ECD) (μs), and is the phase-code parameter, which is 0 for positive phase-coded (PC) pulses and for negative PC pulses. The eLoran pulse signal is generated according to Equation (1). The timedomain and frequency-domain characteristics of a single eLoran signal are shown in Figure 1. It can be seen that the leading edge of the signal first rises rapidly, reaching a peak at 65 μs, and then falls rapidly. The eLoran system can be a timing station or a navigation chain. The timing station generally has only one master station, and the navigation chain generally consists of one master station and two to five secondary stations. The master station is generally represented by M, and the secondary station is represented by letters such as W, X, Y, and Z. To distinguish between the signals of the master and secondary stations, each group of the master station transmits nine eLoran pulse signals. The interval between the first eight pulse signals is 1 ms, and that between the eighth and ninth pulse signals is 2 ms. The secondary station transmits only eight pulses. The pulse-group signals of the master and secondary stations are shown in Figure 2. The eLoran system can be a timing station or a navigation chain. The timing station generally has only one master station, and the navigation chain generally consists of one master station and two to five secondary stations. The master station is generally represented by M, and the secondary station is represented by letters such as W, X, Y, and Z. To distinguish between the signals of the master and secondary stations, each group of the master station transmits nine eLoran pulse signals. The interval between the first eight pulse signals is 1 ms, and that between the eighth and ninth pulse signals is 2 ms. The secondary station transmits only eight pulses. The pulse-group signals of the master and secondary stations are shown in Figure 2.
where is time (μs), A is a normalization constant related to the peak amplitude, is the carrier frequency of 100 kHz, is the envelope-to-cycle difference (ECD) (μs), and is the phase-code parameter, which is 0 for positive phase-coded (PC) pulses and for negative PC pulses. The eLoran pulse signal is generated according to Equation (1). The timedomain and frequency-domain characteristics of a single eLoran signal are shown in Figure 1. It can be seen that the leading edge of the signal first rises rapidly, reaching a peak at 65 μs, and then falls rapidly. The eLoran system can be a timing station or a navigation chain. The timing station generally has only one master station, and the navigation chain generally consists of one master station and two to five secondary stations. The master station is generally represented by M, and the secondary station is represented by letters such as W, X, Y, and Z. To distinguish between the signals of the master and secondary stations, each group of the master station transmits nine eLoran pulse signals. The interval between the first eight pulse signals is 1 ms, and that between the eighth and ninth pulse signals is 2 ms. The secondary station transmits only eight pulses. The pulse-group signals of the master and secondary stations are shown in Figure 2. The interval between the eLoran pulse groups is called the group-repetition interval (GRI). The GRI cycle range is defined as 40,000-99,990 µs, with a resolution of 10 µs. The Remote Sens. 2022, 14, 250 5 of 20 eLoran system continuously emits pulse-group signals according to its own GRI. The phase-code cycles are listed in Table 1. The two code cycles are alternately modulated with the signal; "+" indicates a phase code of 0 radians, and "−" indicates a phase code of π radians. Each pulse signal in the eLoran pulse group of the master and secondary stations is phase-coded. The third to eighth pulse signals of each pulse group are used for the Loran data channel (LDC), which is used for transmitting system status, time information, etc. The modulation mode adopts Eurofix data-modulation technology [32,33]. Table 1. Phase codes of eLoran signals.

Ground
Master Secondary In the eLoran system, the sixth zero crossing (30 µs) of the first pulse of the pulse group is designated as the SZC. The SZC has two key functions: (i) it serves as the pulse-time reference (PTR) for measuring the signal specification of the eLoran broadcasting station; (ii) it allows the eLoran receiver to measure the TOA of the signal to realize high-precision positioning and timing.

Skywave Interference and Cycle Identification
The eLoran signal belongs to an LW frequency band with a carrier frequency of 100 kHz. According to the propagation theory of this frequency band, the eLoran signal mainly reaches the receiving point through two paths: (i) as a groundwave signal (surface wave) transmitted along the surface of the Earth, and (ii) as a skywave signal (space wave) produced by one or more reflections between the ionosphere and the ground when the signal is projected to the ionosphere at a certain elevation angle. Generally, only onehop skywave is considered for receiving points with a receiving distance of less than 2000 km [34]. The propagation path of the eLoran signal is shown in Figure 3. The interval between the eLoran pulse groups is called the group-repetition interval (GRI). The GRI cycle range is defined as 40,000-99,990 μs, with a resolution of 10 μs. The eLoran system continuously emits pulse-group signals according to its own GRI. The phase-code cycles are listed in Table 1. The two code cycles are alternately modulated with the signal; "+" indicates a phase code of 0 radians, and "−" indicates a phase code of π radians. Each pulse signal in the eLoran pulse group of the master and secondary stations is phase-coded. The third to eighth pulse signals of each pulse group are used for the Loran data channel (LDC), which is used for transmitting system status, time information, etc. The modulation mode adopts Eurofix data-modulation technology [32,33].

Ground
Master In the eLoran system, the sixth zero crossing (30 μs) of the first pulse of the pulse group is designated as the SZC. The SZC has two key functions: (i) it serves as the pulsetime reference (PTR) for measuring the signal specification of the eLoran broadcasting station; (ii) it allows the eLoran receiver to measure the TOA of the signal to realize highprecision positioning and timing.

Skywave Interference and Cycle Identification
The eLoran signal belongs to an LW frequency band with a carrier frequency of 100 kHz. According to the propagation theory of this frequency band, the eLoran signal mainly reaches the receiving point through two paths: (i) as a groundwave signal (surface wave) transmitted along the surface of the Earth, and (ii) as a skywave signal (space wave) produced by one or more reflections between the ionosphere and the ground when the signal is projected to the ionosphere at a certain elevation angle. Generally, only one-hop skywave is considered for receiving points with a receiving distance of less than 2000 km [34]. The propagation path of the eLoran signal is shown in Figure 3.  The composite wave comprising the eLoran groundwave and skywave signals received in space can be expressed as: where ( ) is the standard eLoran signal with a normalized amplitude, is the amplitude of the groundwave signal, is the amplitude of the skywave signal, ∆ is the delay difference of the skywave relative to the groundwave, and ( ) is noise. The composite wave comprising the eLoran groundwave and skywave signals received in space can be expressed as: where s(t) is the standard eLoran signal with a normalized amplitude, A is the amplitude of the groundwave signal, A s is the amplitude of the skywave signal, ∆T is the delay difference of the skywave relative to the groundwave, and n(t) is noise.
The skywave-to-groundwave ratio (SGR) is the amplitude ratio of the skywave relative to the groundwave, expressed as: Remote Sens. 2022, 14, 250 6 of 20 The time-delay difference ∆T of the skywave relative to the groundwave is expressed as: where T TOA and T s are the propagation delays of the groundwave and skywave, respectively. Figure 4 shows two typical eLoran received-signal waveforms. The first case is where the phase difference is 0 • due to the time delay of the skywave relative to the groundwave, which increases the amplitude of the composite signal. The second case is where the phase difference is 180 • due to the time delay of the skywave relative to the groundwave, which reduces the amplitude of the composite signal. In both cases, the groundwave is distorted, and the signal amplitude exhibits two extremes. This significantly affects the cycle identification result of the eLoran groundwave. The skywave-to-groundwave ratio (SGR) is the amplitude ratio of the skywave relative to the groundwave, expressed as: The time-delay difference ∆ of the skywave relative to the groundwave is expressed as: where and are the propagation delays of the groundwave and skywave, respectively. Figure 4 shows two typical eLoran received-signal waveforms. The first case is where the phase difference is 0° due to the time delay of the skywave relative to the groundwave, which increases the amplitude of the composite signal. The second case is where the phase difference is 180° due to the time delay of the skywave relative to the groundwave, which reduces the amplitude of the composite signal. In both cases, the groundwave is distorted, and the signal amplitude exhibits two extremes. This significantly affects the cycle identification result of the eLoran groundwave.  Cycle identification is the process of finding the SZC of the first pulse of the groundwave after the eLoran pulse acquisition. Ideally, only the pure groundwave signal is received, and the SZC can be detected by measuring the peak value (65 μs) of the eLoran groundwave envelope. However, the skywave is generated due to reflection from the ionosphere, and it always arrives at the receiver a certain cycle after the groundwave does. Therefore, the eLoran groundwave signal is susceptible to interference by the skywave delay [21,35]. According to the minimum receiver performance standards for marine eLoran receiving equipment issued by the Radio Technical Commission for Maritime Services (RTCM), the time delay range ∆t of the skywave relative to the groundwave is 37.5-1500 μs, and the amplitude ratio SGR of the skywave relative to the groundwave is 12-26 dB [36]. Therefore, to consider both the stability and the SNR, the positive zero crossing (30 μs) in the third week is usually selected as the SZC [37].

Time-Frequency Domain Cycle Identification Method
Because the phase of the first pulse signal in the eLoran pulse group is "+" and there is no data modulation, cycle identification is performed using the first pulse in the eLoran pulse group. To overcome the low identification probability and poor anti-interference ability of existing cycle identification methods, we propose a time-frequency domain method to identify the SZC of the eLoran signal. The frequency-domain method uses Cycle identification is the process of finding the SZC of the first pulse of the groundwave after the eLoran pulse acquisition. Ideally, only the pure groundwave signal is received, and the SZC can be detected by measuring the peak value (65 µs) of the eLoran groundwave envelope. However, the skywave is generated due to reflection from the ionosphere, and it always arrives at the receiver a certain cycle after the groundwave does. Therefore, the eLoran groundwave signal is susceptible to interference by the skywave delay [21,35]. According to the minimum receiver performance standards for marine eLoran receiving equipment issued by the Radio Technical Commission for Maritime Services (RTCM), the time delay range ∆t of the skywave relative to the groundwave is 37.5-1500 µs, and the amplitude ratio SGR of the skywave relative to the groundwave is 12-26 dB [36]. Therefore, to consider both the stability and the SNR, the positive zero crossing (30 µs) in the third week is usually selected as the SZC [37].

Time-Frequency Domain Cycle Identification Method
Because the phase of the first pulse signal in the eLoran pulse group is "+" and there is no data modulation, cycle identification is performed using the first pulse in the eLoran pulse group. To overcome the low identification probability and poor anti-interference ability of existing cycle identification methods, we propose a time-frequency domain method to identify the SZC of the eLoran signal. The frequency-domain method uses spectrum division to estimate the approximate starting positions of the skywave and groundwave, which is also the cycle identification range; the time-domain method uses the peak-to-peak ratio and waveform matching to determine the SZC accurately. The specific implementation process of the method is shown in Figure 5. We assume that the received signal is as follows: where r SGR is the amplitude ratio of the skywave to the groundwave, T 1 is the delay difference of the groundwave relative to the local standard eLoran signal, T 2 is the delay difference of the skywave relative to the local standard eLoran signal, ∆T = T 2 − T 1 is the delay difference of the eLoran skywave relative to groundwave, and n(t) is noise in the signal.
Remote Sens. 2022, 14, x FOR PEER REVIEW 7 of 20 spectrum division to estimate the approximate starting positions of the skywave and groundwave, which is also the cycle identification range; the time-domain method uses the peak-to-peak ratio and waveform matching to determine the SZC accurately. The specific implementation process of the method is shown in Figure 5. We assume that the received signal is as follows: where is the amplitude ratio of the skywave to the groundwave, is the delay difference of the groundwave relative to the local standard eLoran signal, is the delay difference of the skywave relative to the local standard eLoran signal, ∆ = − is the delay difference of the eLoran skywave relative to groundwave, and ( ) is noise in the signal.

Digital Bandpass Filtering and Linear Digital Average
Before the cycle identification, the received signal must be filtered to suppress noise and interference and improve the SNR, to obtain relatively clean composite signals of the groundwave and skywave.
From the eLoran frequency band, the cycle identification results are mainly affected by noise and continuous wave interference (CWI). We used BPFs to suppress out-of-band noise and interference. The BPF includes a Hamming window with a bandwidth of 30 kHz and an order of 128, and its out-of-band suppression capability can exceed 30 dB. The bandpass-filtered signal is: where ( ) is the residual noise in the bandpass-filtered signal.
In the eLoran frequency band, the cycle identification results are mainly affected by noise. After bandpass-filtering, we used LDA to suppress in-band noise, which not only suppresses the noise and improves the SNR of the received signal, but also repairs the eLoran signal and reduces the signal distortion. The LDA is the cumulative average processing of the received signal for M times according to the GRI cycle. Assuming that there are M signals with a length equal to the GRI, the LDA result can be expressed as: where ( ) is the residual noise in the signal after LDA. When the number of LDA is M, the gain of the SNR can be increased by 10 log . We chose the linear digital average of M = 64 to ensure that the SNR of the signal can be improved by 18 dB. A relatively clean composite signal with a length equal to the GRI can be obtained using the linear digital average. The cycle identification method in the timefrequency domain was subsequently used to identify the SZC.

Digital Bandpass Filtering and Linear Digital Average
Before the cycle identification, the received signal must be filtered to suppress noise and interference and improve the SNR, to obtain relatively clean composite signals of the groundwave and skywave.
From the eLoran frequency band, the cycle identification results are mainly affected by noise and continuous wave interference (CWI). We used BPFs to suppress out-of-band noise and interference. The BPF includes a Hamming window with a bandwidth of 30 kHz and an order of 128, and its out-of-band suppression capability can exceed 30 dB. The bandpass-filtered signal is: where n 1 (t) is the residual noise in the bandpass-filtered signal.
In the eLoran frequency band, the cycle identification results are mainly affected by noise. After bandpass-filtering, we used LDA to suppress in-band noise, which not only suppresses the noise and improves the SNR of the received signal, but also repairs the eLoran signal and reduces the signal distortion. The LDA is the cumulative average processing of the received signal for M times according to the GRI cycle. Assuming that there are M signals with a length equal to the GRI, the LDA result can be expressed as: where n 2 (t) is the residual noise in the signal after LDA. When the number of LDA is M, the gain of the SNR can be increased by 10 log 10 M. We chose the linear digital average of M = 64 to ensure that the SNR of the signal can be improved by 18 dB. A relatively clean composite signal with a length equal to the GRI can be obtained using the linear digital average. The cycle identification method in the time-frequency domain was subsequently used to identify the SZC.

Spectrum Division in Frequency Domain
The fast Fourier transform (FFT) was used to obtain the signal spectrum after LDA, and the approximate starting positions of the skywave and groundwave were determined via spectrum division. The FFT of x LDA (t) is: where S( f ) is the FFT result of the eLoran signal and N( f ) is the FFT result of the noise. The spectrum-division method, proposed by Bian and Last, is based on high-resolution time delay estimation under multipath conditions [38]. When Equation (8) is divided by the Fourier transform S 0 ( f ) of the standard eLoran signal, the spectrum-division result is: where k g is the amplitude ratio of the received groundwave signal to the standard signal, and W( f ) is the result of dividing the noise by the standard eLoran signal. According to Euler's formula: The Equation (10) is a sine wave with variable f , and the frequency of sine wave is T. Therefore, the result of the spectrum division of the received signal in Equation (9) becomes an estimate of two continuous wave frequencies T 1 and T 2 . The amplitude ratio r SGR is the ratio of the amplitudes of the two continuous waves. Figure 6 shows the results of spectrum division without noise and with SNR = 5 dB. The waveform of the spectrum division result in the absence of noise was composed of continuous waves formed by different time delays, which is consistent with the theoretical analysis. However, after noise was added, the out-of-band noise was amplified by the spectrum division, which increased the out-of-band noise of the eLoran signal and reduced the SNR sharply. The waveform of the two continuous waves was difficult to distinguish. The Hamming window offered the advantages of low side-lobe attenuation and strong out-of-band suppression. Therefore, a Hamming window with a bandwidth of 50 kHz was selected to suppress out-of-band noise in the spectrum-division results. After windowing, Equation (9) became: where H( f ) is the window function.

. Spectrum Division in Frequency Domain
The fast Fourier transform (FFT) was used to obtain the signal spectrum after LDA, and the approximate starting positions of the skywave and groundwave were determined via spectrum division. The FFT of (t) is: where ( ) is the FFT result of the eLoran signal and ( ) is the FFT result of the noise. The spectrum-division method, proposed by Bian and Last, is based on high-resolution time delay estimation under multipath conditions [38]. When Equation (8) is divided by the Fourier transform ( ) of the standard eLoran signal, the spectrum-division result is: where is the amplitude ratio of the received groundwave signal to the standard signal, and ( ) is the result of dividing the noise by the standard eLoran signal. According to Euler's formula: The Equation (10) is a sine wave with variable , and the frequency of sine wave is . Therefore, the result of the spectrum division of the received signal in Equation (9) becomes an estimate of two continuous wave frequencies and . The amplitude ratio is the ratio of the amplitudes of the two continuous waves. Figure 6 shows the results of spectrum division without noise and with SNR = 5 dB. The waveform of the spectrum division result in the absence of noise was composed of continuous waves formed by different time delays, which is consistent with the theoretical analysis. However, after noise was added, the out-of-band noise was amplified by the spectrum division, which increased the out-of-band noise of the eLoran signal and reduced the SNR sharply. The waveform of the two continuous waves was difficult to distinguish. The Hamming window offered the advantages of low side-lobe attenuation and strong out-of-band suppression. Therefore, a Hamming window with a bandwidth of 50 kHz was selected to suppress out-of-band noise in the spectrum-division results. After windowing, Equation (9) became: where ( ) is the window function.
(a) (b) Finally, the inverse fast Fourier transform (IFFT) was used to estimate the parameters of the skywave and groundwave [39,40]. The IFFT of Equation (11) is: (12) where F −1 represents the IFFT. It can be seen from Equation (12) that δ(t) is the impulse response function resulting from the IFFT of the synthetic signals of the groundwave and skywave. The time delays T 1 and T 2 and amplitude ratio r SGR can be estimated by detecting the time-domain pulse peak value after IFFT. However, the time delay only represents the approximate starting positions of the skywave and groundwave, and the estimated position has a bias; thus, the SZC had to be further determined in the time domain.

Peak-to-Peak Ratio and Waveform Matching in Time Domain
The result of spectrum division determined the starting position range T 1 -T 2 of the groundwave and the skywave in the composite wave. Therefore, we determined the position of the SZC using the peak-to-peak ratio and waveform methods in the time domain.
The peak-to-peak ratio method first determined all positive zero crossing points (t) in the range of T 1 -T 2 and then calculated the ratio of the positive peak value in the next week (t + 2.5 µs) and the previous week (t − 7.5 µs) of each positive zero crossing. If the calculated peak-to-peak ratio was close to the peak-to-peak ratio of the standard eLoran signal a certain range, then it was considered as a possible SZC [23]. The theoretical formula for calculating the peak-to-peak ratio is: The peak-to-peak ratios of positive zero crossings in the theoretical calculation are displayed in Table 2. It can be seen that the peak-to-peak ratio at 30 µs was 1.5338. The peak-to-peak ratios of all positive zero crossings in the range T 1 -T 2 after LDA are calculated first; then, the errors of the peak-to-peak ratios are calculated: where h n (t) is the peak-to-peak ratio of zero crossings, and n is the number of zero crossings. Because the eLoran signal was distorted during propagation and signal processing, which can lead to deviations in the peak-to-peak ratios, when ∆h(n) < 0.3, we determined that the zero crossing n was the closest to the SZC. Thus, the peak-to-peak ratio results contained multiple possible SZCs. Therefore, we used waveform matching to calculate the root mean square (RMS) error between the received signal and the standard waveform over a cycle of time. The minimum value of the RMS error corresponded to the zero crossing to determine the SZC. The waveform-matching calculation formula was as follows: where s(t) is the standard eLoran signal, s n (t) is the received eLoran signal, and t 1 -t 2 is the waveform-matching range. To prevent the influence of SWI, we adopted a waveform matching range of 10-50 µs. Therefore, when σ(n) was the smallest, the corresponding n was the SZC.

Results
In this section, we analyze the effectiveness and performance of the joint timefrequency cycle identification method through simulation. Then, we validate the cycle identification method by receiving actual eLoran signals at different locations.

Validation of Cycle Identification Method
To verify the effectiveness of the proposed method, we performed a simulation, the parameters of which were set as follows. The GRI of the eLoran signal generated by the simulation was 60 ms. After adding noise to the signal, the SNR of groundwave was −5 dB; the amplitude ratio of the skywave to the groundwave was SGR = 10 dB; the time-delay difference between the skywave and the groundwave was ∆T = 62.5 µs; and the sampling rate of the simulation process signal was f s = 2 MHz. The simulation results are shown in Figure 7. where ( ) is the standard eLoran signal, ( ) is the received eLoran signal, andis the waveform-matching range. To prevent the influence of SWI, we adopted a waveform matching range of 10-50 μs. Therefore, when ( ) was the smallest, the corresponding n was the SZC.

Results
In this section, we analyze the effectiveness and performance of the joint time-frequency cycle identification method through simulation. Then, we validate the cycle identification method by receiving actual eLoran signals at different locations.

Validation of Cycle Identification Method
To verify the effectiveness of the proposed method, we performed a simulation, the parameters of which were set as follows. The GRI of the eLoran signal generated by the simulation was 60 ms. After adding noise to the signal, the SNR of groundwave was −5 dB; the amplitude ratio of the skywave to the groundwave was SGR = 10 dB; the timedelay difference between the skywave and the groundwave was ∆ = 62.5 μs; and the sampling rate of the simulation process signal was = 2 MHz. The simulation results are shown in Figure 7. As can be seen from Figure 7a, the composite wave of the skywave and groundwave generated in the simulation is submerged in noise, and the groundwave signal is severely distorted. Figure 7b presents the bandpass-filtered signal. It can be seen that although the noise in the signal is suppressed, the signal distortion is relatively large. Figure 7c shows the signal after LDA. It can be seen that the signal not only became cleaner, but the signal waveform also became complete through repair. The spectrum-division method was then used to obtain the initial positions of the groundwave and skywave on the linearly averaged signal, as shown in Figure 7d. The range of the SZC was observed to be 320-445. The results of calculating the peak-to-peak ratio errors of all zero crossings in the determined range are shown in Figure 7e. As observed, the peak-to-peak ratio errors of zero crossings 3 and 4 were less than 0.3. Waveform matching was performed on the two zero-crossing points, and the results are shown in Figure 7f,g. The waveform-matching coefficients of the two zero-crossings were calculated to be 0.0456621 and 0.1117095, respectively. Thus, zero crossing 3 was the SZC. The SZC determined using the cycle identification method was consistent with that set based on the simulation, which proves the effectiveness of the method.

Performance Analysis of Cycle Identification Method
Noise and SWI are the most important interference sources that affect cycle identification, and they are inevitable. According to the minimum performance standard for  As can be seen from Figure 7a, the composite wave of the skywave and groundwave generated in the simulation is submerged in noise, and the groundwave signal is severely distorted. Figure 7b presents the bandpass-filtered signal. It can be seen that although the noise in the signal is suppressed, the signal distortion is relatively large. Figure 7c shows the signal after LDA. It can be seen that the signal not only became cleaner, but the signal waveform also became complete through repair. The spectrum-division method was then used to obtain the initial positions of the groundwave and skywave on the linearly averaged signal, as shown in Figure 7d. The range of the SZC was observed to be 320-445. The results of calculating the peak-to-peak ratio errors of all zero crossings in the determined range are shown in Figure 7e. As observed, the peak-to-peak ratio errors of zero crossings 3 and 4 were less than 0.3. Waveform matching was performed on the two zero-crossing points, and the results are shown in Figure 7f,g. The waveform-matching coefficients of the two zero-crossings were calculated to be 0.0456621 and 0.1117095, respectively. Thus, zero crossing 3 was the SZC. The SZC determined using the cycle identification method was consistent with that set based on the simulation, which proves the effectiveness of the method.

Performance Analysis of Cycle Identification Method
Noise and SWI are the most important interference sources that affect cycle identification, and they are inevitable. According to the minimum performance standard for eLoran receivers, the requirements for an eLoran receiver to function normally were SNR ≥ −10 dB, SGR < 12 dB, and ∆T > 37.5 µs [36]. Therefore, to verify the performance of the proposed cycle identification method, we analyzed the skywave and groundwave parameter-estimation performance of the spectrum-division method and the probability of cycle identification based on the joint time-frequency domain under different SNRs.
We simulated the parameter-estimation performance of the spectrum-division method under SNR = 5, 0, and -10 dB. For SGR = 5-10 dB and ∆T = 37-150 µs, the amplitude-ratio error and delay-difference error after spectrum estimation are shown in Figure 8.
Remote Sens. 2022, 14, x FOR PEER REVIEW 12 of 20 eLoran receivers, the requirements for an eLoran receiver to function normally were SNR ≥ −10 dB, SGR < 12 dB, and ∆ > 37.5 μs [36]. Therefore, to verify the performance of the proposed cycle identification method, we analyzed the skywave and groundwave parameter-estimation performance of the spectrum-division method and the probability of cycle identification based on the joint time-frequency domain under different SNRs. We simulated the parameter-estimation performance of the spectrum-division method under SNR = 5, 0, and -10 dB. For SGR = 5-10 dB and ∆ = 37-150 μs, the amplitude-ratio error and delay-difference error after spectrum estimation are shown in Figure  8.  Figure 8 shows that with a decrease in SNR, the estimation error of the time-delay difference and the amplitude ratio of the skywave to the groundwave based on the spectrumdivision method increased. When SNR = 5 dB, the amplitude-ratio estimation error was less than 1 dB, and the delay-difference estimation error was less than 2 µs. When SNR = 0 dB, the amplitude-ratio estimation error was less than 1 dB, and the delay-difference estimation error was less than 10 µs. When SNR = −10 dB, the amplitude-ratio estimation error was less than 2 dB, and the delay-difference estimation error was less than 10 µs in most cases; however, when ∆T > 50 µs and SGR < −3 dB, the delay-difference estimation error was greater than 10 µs. In the actual received signal, the SNR may deteriorate, resulting in a greater estimation error in the delay difference. Therefore, based on the simulation results, the spectrum-division method cannot be used for complete cycle identification; more accurate cycle identification based on spectrum division using the time-domain method is necessary. To address the large estimation error of the spectrum-division method under a low SNR, we proposed a cycle identification method based on the combined time-frequency domain. Through spectral division, we further improved the accuracy of cycle identification in the time domain.
Remote Sens. 2022, 14, x FOR PEER REVIEW 13 of 20 Figure 8 shows that with a decrease in SNR, the estimation error of the time-delay difference and the amplitude ratio of the skywave to the groundwave based on the spectrum-division method increased. When SNR = 5 dB, the amplitude-ratio estimation error was less than 1 dB, and the delay-difference estimation error was less than 2 μs. When SNR = 0 dB, the amplitude-ratio estimation error was less than 1 dB, and the delay-difference estimation error was less than 10 μs. When SNR = −10 dB, the amplitude-ratio estimation error was less than 2 dB, and the delay-difference estimation error was less than 10 μs in most cases; however, when ∆ > 50 μs and SGR < −3 dB, the delay-difference estimation error was greater than 10 μs. In the actual received signal, the SNR may deteriorate, resulting in a greater estimation error in the delay difference. Therefore, based on the simulation results, the spectrum-division method cannot be used for complete cycle identification; more accurate cycle identification based on spectrum division using the time-domain method is necessary. To address the large estimation error of the spectrumdivision method under a low SNR, we proposed a cycle identification method based on the combined time-frequency domain. Through spectral division, we further improved the accuracy of cycle identification in the time domain.
We simulated the parameter estimation performance of the cycle identification method under SNR = 0, −5, −10, and −13 dB. Figure 9 shows the error rate of cycle identification when SGR = 5-10 dB and ∆ = 37-150 μs. The simulation results in Figure 9 show that when SGR ≤ 23 dB and SNR ≥ 0 dB, the success rate of cycle identification was 100%. When SGR ≤ 23 dB and SNR ≥ −10 dB, the success rate was better than 75%. When SGR ≤ 23 dB and SNR ≥ −13 dB, the success rate  The simulation results in Figure 9 show that when SGR ≤ 23 dB and SNR ≥ 0 dB, the success rate of cycle identification was 100%. When SGR ≤ 23 dB and SNR ≥ −10 dB, the success rate was better than 75%. When SGR ≤ 23 dB and SNR ≥ −13 dB, the success rate was better than 55%. The cycle identification was mainly affected by the relative amplitudes of the skywave and groundwave, especially under a low SNR. Although the success rate of cycle identification decreased with an increase in SGR, accurate cycle identification could be achieved via multiple cycle identification calculations. Therefore, the-simulation results showed that the proposed cycle identification method is effective and offers a high identification probability.

Experimental Verification of Cycle Identification Method
We validated the cycle identification method with actual signals transmitted by a real eLoran system. The test scheme is shown in Figure 10. The eLoran receiver obtains the signal transmitted by the eLoran station through the antenna, filters and amplifies the signal received by the radiofrequency (RF) unit, and converts it into a digital signal through the A/D chip. Next, the digital-signal-processing unit of the receiver captures the eLoran pulse-group signal to obtain the approximate position of the pulse signal. Finally, the proposed cycle identification method was used to identify the SZC of the first eLoran signal of the pulse group, and the identification result was saved in a computer. The digital signal processing in Figure 10 was completed on an FPGA chip. The test platform built according to the test scheme is shown in Figure 11. was better than 55%. The cycle identification was mainly affected by the relative amplitudes of the skywave and groundwave, especially under a low SNR. Although the success rate of cycle identification decreased with an increase in SGR, accurate cycle identification could be achieved via multiple cycle identification calculations. Therefore, the-simulation results showed that the proposed cycle identification method is effective and offers a high identification probability.

Experimental Verification of Cycle Identification Method
We validated the cycle identification method with actual signals transmitted by a real eLoran system. The test scheme is shown in Figure 10. The eLoran receiver obtains the signal transmitted by the eLoran station through the antenna, filters and amplifies the signal received by the radiofrequency (RF) unit, and converts it into a digital signal through the A/D chip. Next, the digital-signal-processing unit of the receiver captures the eLoran pulse-group signal to obtain the approximate position of the pulse signal. Finally, the proposed cycle identification method was used to identify the SZC of the first eLoran signal of the pulse group, and the identification result was saved in a computer. The digital signal processing in Figure 10 was completed on an FPGA chip. The test platform built according to the test scheme is shown in Figure 11.  The actual received signal is the signal emitted by the BPL timing system in China. The coordinates of the BPL timing system are (034 ∘ 57 N, 109 ∘ 33 E). Its GRI and transmission power are 60 ms and 2 MW, respectively. In July 2021, we selected three test sites within 1000 km of the system to receive actual signals. The coordinates, great circle distances, and electric-field strengths of the test sites are listed in Table 3. was better than 55%. The cycle identification was mainly affected by the relative amplitudes of the skywave and groundwave, especially under a low SNR. Although the success rate of cycle identification decreased with an increase in SGR, accurate cycle identification could be achieved via multiple cycle identification calculations. Therefore, the-simulation results showed that the proposed cycle identification method is effective and offers a high identification probability.

Experimental Verification of Cycle Identification Method
We validated the cycle identification method with actual signals transmitted by a real eLoran system. The test scheme is shown in Figure 10. The eLoran receiver obtains the signal transmitted by the eLoran station through the antenna, filters and amplifies the signal received by the radiofrequency (RF) unit, and converts it into a digital signal through the A/D chip. Next, the digital-signal-processing unit of the receiver captures the eLoran pulse-group signal to obtain the approximate position of the pulse signal. Finally, the proposed cycle identification method was used to identify the SZC of the first eLoran signal of the pulse group, and the identification result was saved in a computer. The digital signal processing in Figure 10 was completed on an FPGA chip. The test platform built according to the test scheme is shown in Figure 11.  The actual received signal is the signal emitted by the BPL timing system in China. The coordinates of the BPL timing system are (034 ∘ 57 N, 109 ∘ 33 E). Its GRI and transmission power are 60 ms and 2 MW, respectively. In July 2021, we selected three test sites within 1000 km of the system to receive actual signals. The coordinates, great circle distances, and electric-field strengths of the test sites are listed in Table 3. The actual received signal is the signal emitted by the BPL timing system in China. The coordinates of the BPL timing system are (034 • 57 N, 109 • 33 E). Its GRI and transmission power are 60 ms and 2 MW, respectively. In July 2021, we selected three test sites within 1000 km of the system to receive actual signals. The coordinates, great circle distances, and electric-field strengths of the test sites are listed in Table 3. The environments of the three selected test points were different. Test 1 was in a closed indoor environment, and Figure 12 shows the collected eLoran signal. It can be seen that although the receiving distance was short and there was only a groundwave signal, the signal was full of pulse interference. Test 2 was in the surrounding open environment, and Figure 13 shows the collected eLoran signal. It can be seen that there was SWI in the signal, which made the second half of the groundwave distorted. Test 3 was in a mountainous location at a large circle distance, and Figure 14 shows the collected eLoran signal. Not only was the skywave interference for this test strong, but the SNR was also very small, which made the groundwave signal severely distorted.  The environments of the three selected test points were different. Test 1 was in a closed indoor environment, and Figure 12 shows the collected eLoran signal. It can be seen that although the receiving distance was short and there was only a groundwave signal, the signal was full of pulse interference. Test 2 was in the surrounding open environment, and Figure 13 shows the collected eLoran signal. It can be seen that there was SWI in the signal, which made the second half of the groundwave distorted. Test 3 was in a mountainous location at a large circle distance, and Figure 14 shows the collected eLoran signal. Not only was the skywave interference for this test strong, but the SNR was also very small, which made the groundwave signal severely distorted.   The environments of the three selected test points were different. Test 1 was in a closed indoor environment, and Figure 12 shows the collected eLoran signal. It can be seen that although the receiving distance was short and there was only a groundwave signal, the signal was full of pulse interference. Test 2 was in the surrounding open environment, and Figure 13 shows the collected eLoran signal. It can be seen that there was SWI in the signal, which made the second half of the groundwave distorted. Test 3 was in a mountainous location at a large circle distance, and Figure 14 shows the collected eLoran signal. Not only was the skywave interference for this test strong, but the SNR was also very small, which made the groundwave signal severely distorted.  The results after the received signals of the three test sites were processed via the joint time-frequency domain cycle identification method are shown in Figures 15-17. The number of cycles identified is 700. The peak-to-peak ratio errors of the three test points were 0.0124, 0.0676, and 0.0266, respectively, and the errors were all less than 0.3, meeting the peak-to-peak ratio evaluation conditions. However, as the distance increased, the dispersion of the eLoran signal became significant in the propagation process, leading to a larger standard deviation of the calculated result of the peak-to-peak ratio. The standard deviations of the peak-to-peak ratios were 0.0199, 0.0292, and 0.0938, but these results had no influence on the evaluation of the peak-to-peak ratio. The cycle identification results of the actual signals showed very high cycle identification success rates of tests 1 and 2 and only one identification error at test point 1. Although test point 3 was in a far-off mountainous area, the success rate of cycle identification for this test reached 72.86%. Comparing the experimental results with the simulation results, it can be seen that the cycle identification success rate of the actual measurement results is slightly worse. The main reason is that the actual eLoran signal received contains pulse interference, continuous wave interference and other interference, which would cause the quality of the eLoran signal to deteriorate and affect the success rate of cycle identification. In practical applications, accurate cycle identification can be achieved through multiple identification calculations. Therefore, the experimental test fully verifies the practicality of the proposed joint timefrequency domain eLoran cycle identification method. The results after the received signals of the three test sites were processed via the joint time-frequency domain cycle identification method are shown in Figures 15-17. The number of cycles identified is 700. The peak-to-peak ratio errors of the three test points were 0.0124, 0.0676, and 0.0266, respectively, and the errors were all less than 0.3, meeting the peak-to-peak ratio evaluation conditions. However, as the distance increased, the dispersion of the eLoran signal became significant in the propagation process, leading to a larger standard deviation of the calculated result of the peak-to-peak ratio. The standard deviations of the peak-to-peak ratios were 0.0199, 0.0292, and 0.0938, but these results had no influence on the evaluation of the peak-to-peak ratio. The cycle identification results of the actual signals showed very high cycle identification success rates of tests 1 and 2 and only one identification error at test point 1. Although test point 3 was in a far-off mountainous area, the success rate of cycle identification for this test reached 72.86%. Comparing the experimental results with the simulation results, it can be seen that the cycle identification success rate of the actual measurement results is slightly worse. The main reason is that the actual eLoran signal received contains pulse interference, continuous wave interference and other interference, which would cause the quality of the eLoran signal to deteriorate and affect the success rate of cycle identification. In practical applications, accurate cycle identification can be achieved through multiple identification calculations. Therefore, the experimental test fully verifies the practicality of the proposed joint time-frequency domain eLoran cycle identification method. The results after the received signals of the three test sites were processed via the joint time-frequency domain cycle identification method are shown in Figures 15-17. The number of cycles identified is 700. The peak-to-peak ratio errors of the three test points were 0.0124, 0.0676, and 0.0266, respectively, and the errors were all less than 0.3, meeting the peak-to-peak ratio evaluation conditions. However, as the distance increased, the dispersion of the eLoran signal became significant in the propagation process, leading to a larger standard deviation of the calculated result of the peak-to-peak ratio. The standard deviations of the peak-to-peak ratios were 0.0199, 0.0292, and 0.0938, but these results had no influence on the evaluation of the peak-to-peak ratio. The cycle identification results of the actual signals showed very high cycle identification success rates of tests 1 and 2 and only one identification error at test point 1. Although test point 3 was in a far-off mountainous area, the success rate of cycle identification for this test reached 72.86%. Comparing the experimental results with the simulation results, it can be seen that the cycle identification success rate of the actual measurement results is slightly worse. The main reason is that the actual eLoran signal received contains pulse interference, continuous wave interference and other interference, which would cause the quality of the eLoran signal to deteriorate and affect the success rate of cycle identification. In practical applications, accurate cycle identification can be achieved through multiple identification calculations. Therefore, the experimental test fully verifies the practicality of the proposed joint timefrequency domain eLoran cycle identification method.

Discussion
As an important part of the PNT system, the eLoran system can form a powerful supplement and backup to GNSS. Therefore, many countries are developing and improving eLoran systems to improve the service capabilities of the national PNT systems. Compared with GNSS, the accuracy and reliability of the eLoran system needs to be improved. Cycle identification is a key factor in determining the positioning and timing accuracy of eLoran receivers. However, due to noise in the space and the SWI caused by the signal propagation channel, the cycle identification success rate in the eLoran receiver is affected. The existing methods are poor in reliability and anti-interference; thus, this paper proposes a joint time-frequency domain cycle identification method. After the cycle identification range was determined by the spectrum division method, the peak-to-peak ratio and the waveform matching method were used to determine the accurate SZC, which can improve the success rate of cycle identification.
Firstly, based on the analysis of the propagation channel and interference characteristics of the eLoran signal, the principle of the joint time-frequency domain cycle identification method was given. The effectiveness of the method was verified by simulation, and the anti-noise and SWI performance of the method was analyzed. Simulation results showed that this method had a high success rate of cycle identification. When SNR ≥ 0 and SGR ≤ 23 dB, the success rate of cycle identification was 100%; when SNR ≥ −10 and SGR ≤ 23 dB, the success rate of cycle identification was greater than 75%; when SNR ≥ Cycle identification times Peak-to-peak ratio

Discussion
As an important part of the PNT system, the eLoran system can form a powerful supplement and backup to GNSS. Therefore, many countries are developing and improving eLoran systems to improve the service capabilities of the national PNT systems. Compared with GNSS, the accuracy and reliability of the eLoran system needs to be improved. Cycle identification is a key factor in determining the positioning and timing accuracy of eLoran receivers. However, due to noise in the space and the SWI caused by the signal propagation channel, the cycle identification success rate in the eLoran receiver is affected. The existing methods are poor in reliability and anti-interference; thus, this paper proposes a joint time-frequency domain cycle identification method. After the cycle identification range was determined by the spectrum division method, the peak-to-peak ratio and the waveform matching method were used to determine the accurate SZC, which can improve the success rate of cycle identification.
Firstly, based on the analysis of the propagation channel and interference characteristics of the eLoran signal, the principle of the joint time-frequency domain cycle identification method was given. The effectiveness of the method was verified by simulation, and the anti-noise and SWI performance of the method was analyzed. Simulation results showed that this method had a high success rate of cycle identification. When SNR ≥ 0 and SGR ≤ 23 dB, the success rate of cycle identification was 100%; when SNR ≥ −10 and SGR ≤ 23 dB, the success rate of cycle identification was greater than 75%; when SNR ≥ Cycle identification times Peak-to-peak ratio Cycle identification results (μs) Figure 17. Results of cycle identification for Test 3: (a) peak-peak ratio; (b) cycle identification result.

Discussion
As an important part of the PNT system, the eLoran system can form a powerful supplement and backup to GNSS. Therefore, many countries are developing and improving eLoran systems to improve the service capabilities of the national PNT systems. Compared with GNSS, the accuracy and reliability of the eLoran system needs to be improved. Cycle identification is a key factor in determining the positioning and timing accuracy of eLoran receivers. However, due to noise in the space and the SWI caused by the signal propagation channel, the cycle identification success rate in the eLoran receiver is affected. The existing methods are poor in reliability and anti-interference; thus, this paper proposes a joint time-frequency domain cycle identification method. After the cycle identification range was determined by the spectrum division method, the peak-to-peak ratio and the waveform matching method were used to determine the accurate SZC, which can improve the success rate of cycle identification.
Firstly, based on the analysis of the propagation channel and interference characteristics of the eLoran signal, the principle of the joint time-frequency domain cycle identification method was given. The effectiveness of the method was verified by simulation, and the anti-noise and SWI performance of the method was analyzed. Simulation results showed that this method had a high success rate of cycle identification. When SNR ≥ 0 and SGR ≤ 23 dB, the success rate of cycle identification was 100%; when SNR ≥ −10 and SGR ≤ 23 dB, the success rate of cycle identification was greater than 75%; when SNR ≥ −13 and SGR ≤ 23 dB, the success rate of cycle identification was greater than 55%. When the success rate of the current cycle identification method reached 100%, the delay-locked loop method required SNR ≥ 20 dB, the waveform matching method required SNR ≥ 23 dB, and the delayed-addition method required SNR ≥ 24 dB. Although the performance of these methods can be improved by 10-18 dB after filtering, these methods can only meet the requirements under the condition of SGR ≤ 6 dB. Obviously, the cycle identification method proposed in this paper was better than the single cycle identification method for a combination of noise and skywave interference.
Finally, we applied the cycle identification method to the developed eLoran receiver and receive actual signals within 1000 km for testing. The test results show that under the condition of high SNR, the cycle identification success rate can reach 100%. Even in the 862.5-km mountain test environment, the method can achieve a 72.86% cycle recognition success rate. In practical applications, higher identification success rates can also be achieved by counting the results of multiple cycle identifications. Thus, the method has high engineering application value, and it is optimized for the design of the new eLoran navigation and timing receivers. It will be able to improve the accuracy of the application of modern eLoran systems, making it a more reliable backup for GNSS.

Conclusions
In this paper, a cycle identification method based on a joint time-frequency domain was proposed for the standard zero crossing detection of groundwave signals in eLoran positioning and timing receivers. We proved that this method has a good ability of anti-skywave and anti-noise detection and can address the limitation of poor detection performance of the current cycle identification method. This method also has a greater application potential, because it can be used for period identification as well as to obtain the time delay and amplitude of the skywave in the received signal, which is important for the receiver. In the future, we will use this method to separate the skywave and groundwave, in order to reduce the impact of skywave instability on the groundwave measurement. This can improve the TOA measurement accuracy of the eLoran signal and the information demodulation ability of the signal.  Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Restrictions apply to the availability of these data. The ownership of data belongs to the National Time Service Center (NTSC). These data can be available from the corresponding author with the permission of NTSC.