Improving Detection Performance of Ionospheric Disturbances due to Earthquake by Optimization of Sequential Measurement Combination

. Energy generated from earthquake (EQ) is transferred to the ionosphere and results in co-seismic ionospheric disturbances (CID). CID can be observed in the ionospheric combination using L1, L2 frequency carrier phase. As ionospheric trend due to normal conditions such as elevation angle of satellites is generally larger than disturbances, a proper measure is required to extract disturbance signals. Derivative, or sequential combination, is a simple and effective way to remove the normal trend in the ionospheric delay. When using derivative, however, disturbance signals can often be obscured by noise due to its small amplitude. In order to reduce the noise while preserving the time rate of data, and thus to improve signal-to-noise ratio (SNR), we designed a new derivative method using optimization under a couple of assumptions. With simulation data, it is found that N, the number of epochs used for sequential combination, turned out to be the best when N=160 with maximum SNR. Finally, the proposed algorithm’s SNR was compared to that of the previous study which also used derivative method. 120~260% improvements were observed for the proposed method compared to the conventional method.


Introduction
Global Navigation Satellite System (GNSS) measurement has various errors including satellite orbit error, satellite clock bias, receiver clock bias, ionospheric delay, tropospheric delay, multipath error, etc.As the ionospheric delay is proportional to the total electron content (TEC) from satellite to receiver, electron density in the ionosphere can be estimated using GNSS signal.When earthquakes take place, the energy generated from the event is transferred to the atmosphere by groundatmosphere coupling.Acoustic Gravity Wave (AGW) is responsible for the energy transfer, and with energy conservation and exponential decrease of air density with altitude, disturbance waves are magnified through propagation.This energy disturbs neutral particles in the ionosphere, which in turn disturb ions by momentum transfer [1,2].CID is comprised of several types.First type is induced by Rayleigh wave, the seismic wave on the ground which travels at 3.5km/s [3].This CID is easily distinguishable from other types by its high speed.Second type is induced by direct acoustic wave from the ground source.Acoustic wave is generated from the abrupt vertical movement of the epicentre and moves upward and eventually reaches the ionosphere.The ionospheric disturbance by direct acoustic wave also shows horizontal movement, as the acoustic wave path is curved due to the variable temperature profile by altitude.There are also other types, including tsunami-induced gravity wave and ionospheric gravity wave that incur ionospheric disturbances.Among these various types of waves responsible for disturbances, CID by Rayleigh wave is the main interest in this study as it arrives first with dominant disturbance amplitudes.
In order to detect CID, normal trend needs to be removed.Generally, ionospheric delay measurement shows a U-shape along time series as a satellite rises and falls.This is mainly because ionospheric delay gets larger with low elevation angle, as signal travels longer path inside the ionosphere.There are other factors that contribute to the ionosphere such as diurnal change, seasonal change, solar activity, etc.However, these factors have long period (<1mHz) characteristics that their time rate is insubstantial compared to the main trend.
Several options are available to get rid of the normal trend.Commonly used are bandpass filter, moving average, and derivative.Among them, derivative is a relatively simple and intuitive method using linear combination of sequential measurements [4].When there is CID, the derivative of ionospheric delay will be composed of three parts: the derivative of normal trend, noise, and CID.As ionospheric delay due to normal trend changes little for a short time, its derivative is quite small.Therefore, there remain noise and CID in the derivative of ionospheric delay.
If, however, CID amplitude is not big enough, chances are that CID would not stand out of the noise in the time series.Hence improving CID's SNR is a very important issue.Few studies, however, have focused on reducing the noise part of the derivative and thus improving the SNR of the signal.We propose a new method that deals with the issue which can enhance CID's SNR and improve detection performance.

Methodology
GPS L1, L2 carrier phase are denoted in the below.The GPS orbit error, multipath error and inter-frequency bias are assumed to be zero.
Subscript 1 and 2 indicate L1, L2 frequency respectively. is the distance between satellite to receiver,  receiver clock bias,  satellite clock bias,  ionospheric delay error,  tropospheric delay error,  wavelengh of carrier phase,  ambiguity, and  carrier phase noise.Ionospheric combination, or geometry-free combination can be derived in Eq (2).
As the second term on the right side is a bias it disappears by taking derivative.In the end, the n-th derivative of ionopsheric combination is denoted as Eq (3).
CID is included in  1 () .To derive the optimal sequential combination that can reduce the noise part, , while preserving  1 () , two assumptions are adopted.

Assumptions
The first assumption is that the ionospheric delay can be approximated as linear for a short time.This assumption holds true as normal trends in the ionosphere have very long periods (>1000sec).As ionospheric delay by normal trend changes very slowly, it is nearly linear within a short time span.
The other assumption is that the noise in the ionospheric delay is Gaussian random.It means that there is no time correlation in the noise and it follows normal distribution.This assumption was adopted to simplify the calculation.

Minimum Noise Derivative (MND)
Ionospheric combination measurement can be denoted as follows.Notation is changed for concise expression.
where  is ionospheric combination,  true ionospheric delay, and  Gaussian noise with ~(0,  2 ) .Now, Taylor series expansion for n sequential epochs are as follows.
Here, the 2 nd and higher order terms, (ℎ 2 ), are zero in accordance with the first assumption.In order to get the derivative value at the i-th epoch, we will derive   ′ by linear combination of   ,  +1 , … ,  +−1 .Also, as this study deals with 1-sec interval data only, set ℎ equal to 1.
By linear combination, where   ′ is first derivative of ionospheric delay, and   ′ the derivative of noise at i-th epoch.The noise standard deviation (STD) of derivative,   ′ , is root square sum of its coefficients.That is, Therefore, in order to minimize    ′ , the square sum of sequential combination coefficients should be the minimum.Use zero gradient property at extreme values. where, With Eq(9-11) the following equation can be derived.
Also, in the first row of Eq(7) it can be shown that the equation does not change even if the numerator and denominator are both divided by  1 .This means that   ′ is determined by relative size of  1 ,  2 , … ,  −1 and not by their absolute values.Therefore, as a pivot we can set  1 equal to 1.
In the end, the sequential combination which minimizes noise level is, where,   = −6(−1)+12(−1) Then STD of noise derivative,    ′ is Fig. 1 shows relative STD level by N, the number of epochs used for derivative.The y-axis is in log scale.As N increases, the noise level diminishes drastically.With N=100, for example, STD of derivative noise shrinks to 0.35% of its original noise STD, .

Derivative Level to Detect CID
Though normal condition trend changes slowly in the 1 st derivative of ionospheric combination, it eventually makes an overall slant in time series.This is because unlike noise, normal trend accumulates with time for its low frequency.Hence, to completely remove the effect of normal trend, higher order derivatives are needed.We will use the 3 rd derivative to detect CID as it contains no tendency in the normal condition.This is the same derivative level as the previous study [4].

N Selection by Simulation Data
SNR is widely used to indicate the relative amplitude of signal and ambient noise.It is defined as the ratio of signal variance to noise variance.This concept is important in CID detection as well, for SNR is closely related to detection.However, it is hard to calculated the variance of CID by Rayleigh wave as it lasts only for several minutes.Hence, here in this article we define SNR as the ratio of maximum absolute value of CID to STD of noise as shown in Eq(16).Here  indicates the maximum absolute value of CID by Rayleigh wave and   denotes STD of noise without CID.
While the designed derivative ensures the minimum noise level for its derivative, SNR of CID does not monotonically increase with N.This is because CID is not exactly linear for the time span N and induces error, which might curtail the peaks of CID in the derivative.Therefore, we need to find the best N maximizing the SNR of CID.
Simulation data were used to find the best N.For data generation, CID was designed as a sine wave with amplitude 5 and period 225 sec.This period corresponds to the Airy phase of Rayleigh wave [5,6].Noise was modeled as Gaussian normal with STD value 1.For the ionospheric trend by elevation of satellites, 6 hr period sine wave with amplitude 10 was adopted.Here 6 hr period for normal trend is harsher than real trend, as typical normal trend has a 6 hr long U-shape, which can be approximated as one half of 12 hr period sine wave.Fig. 2 shows generated simulation data (black) and its pure CID (top, red).In the bottom plot, MND 3 rd result with N=160 is drawn in red solid line.By using MND 3 rd , CID is extracted with a sufficiently high SNR.Here SNR means the ratio of maximum absolute value in 2.8~3.2hr to STD of 0.5~2.5 hr, as shown in the Fig. 2 bottom plot.To find out the best N for CID detection, a numerical simulation with 100 iterations was performed.With N ranging from 5 to 200 with step size 5, SNR was computed for 100 different simulation data and averaged.Black dots in Fig. 3 represents the result of each data set with a particular N. Red solid line is the average value of 100 iterations.Consequently, N=160 turned out to be the best option to detect CID.

SNR Comparison with Previous Study
To analyze SNR performance of designed algorithm, 2011 Tohoku Earthquake data were used with National Geographic Information Institute (NGII) stations.The total number of stations are 44, and 3 GPS satellites are chosen for their proximity to the epicenter.In Fig. 4 red star denotes epicenter, black squares NGII stations, a filled square KIMC station of NGII.IPP tracks of prn 15,26,27 are drawn with respect to KIMC station.CID by Rayleigh wave has 3.5km/s speed.It starts at 10 min after EQ from the ionospheric height above epicenter, and propagates as a circular wave [3].CID by Rayleigh wave travels 4000 km for 20 min, and dissipates most of its energy at the end of propagation.Therefore, it is safe to say that CID, if detectable, will arrive within 10~30 min time span since EQ.Therefore, the amplitude of CID will be calculated as the largest absolute value in this region.STD of noise will be calculated from 1hr time window before EQ outbreak.In this way we can minimize the effect of elevation angle on noise and guarantee there is no CID in the region.
Previous study by Park used a derivative method called numerical 3 rd order slant TEC derivatives to detect CID with 30-sec data [4].As this method was based on 30-sec data, SNR with 1-sec data using this method shows less desirable results.Therefore, we will compare the SNR of 30-sec data using the previous study with that of 1-sec data using MND.Table 1 is the statistics on average SNR for prn15, 26, and 27.MND showed 266%, 119%, and 226% improvements respectively.

Conclusion
A new derivative method was proposed with 1-sec interval sequential measurement combination.Minimum noise level is guaranteed under the given assumptions for the algorithm, which enhances the SNR of CID.The 3 rd derivative was chosen to remove the normal trend in the ionospheric delay measurements.To find out the best N that maximizes CID detection performance, simulation data was used and N=160 turned out to be the best option.The SNR performance of the designed algorithm was compared with that of previous study and 120~260% improvements were observed.

Fig. 1 .
Fig. 1.Relative STD by N of the designed derivative method.

Fig. 2 .
Fig. 2. Modeled ionospheric delay with CID (black), pure CID (top, red) and MND 3 rd of modelled data (bottm, red).Shades represents the regions from which STD and maximum absolute value were calculated respectably.

Fig. 3 .
Fig. 3. SNR by N of simulation data.N=160 is the best option.

Fig. 5
shows examples of CID detection by each algorithm.Black line represents the result with the previous study, red with MND.Y axis was normalized by STD of 1 hr region before EQ.Maximum absolute value was extracted from the time span of [EQ+10min, EQ+30min], where CID by Rayleigh wave is expected to arrive.These two time regions are notated as shades in the figure.As the Y axis was normalized, maximum values in the plot represent SNR.These values are marked as hollow circles.As shown in Fig. 5, MND shows much improved SNR performance compared to the previous study results.

Fig. 5 .
Fig. 5. SNR comparison between the previous study (Park, black) and the designed algorithm (red).Vertical magenta lines represent EQ outbreak.Refer to Fig.2 for the shaded regions.Each data was normalized by STD of the left shaded region.