I Tomashchuk

I Bendaoud

J-M Jouvard

P Sallamand

The numerical challenges in multiphysical modeling of laser welding with arbitrary Lagrangian-Eulerian method

Keywords: laser welding, metallic alloys, moving mesh, multiphysics

The interaction of high power laser beam with metallic materials produces a number of interconnected phenomena that represent a serious challenge for numerical modeling, especially for creation of auto-consistent models. Additional difficulty consists in lack of data on materials properties at the temperatures superior to their melting point. The present work summarizes the numerical challenges in creation and validation of free-surface models using ALE moving mesh coupled with heat transfer equation and Navier-Stokes fluid flow.

Introduction

ALE method has been successfully used for multiphysical modelling of pulsed [START_REF] Bruyere | Comparison between Phase Field and ALE methods to model the keyhole digging during spot laser welding[END_REF][START_REF] Tomashchuk | Multiphysical modelling of keyhole formation during dissimilar laser welding[END_REF] and continous 2 laser welding, drilling [START_REF] Sharma | A study of hydrodynamics of melt expulsion in pulsed Nd:YAG laser drilling of titanium[END_REF] , direct laser metal deposition [START_REF] Morville | 2D axial-symmetric model for fluid flow and heat transfer in the melting and resolidification of a vertical cylinder[END_REF] and thermal matter ablation [START_REF] Frei | Modeling Thermal Ablation for Material Removal[END_REF] . The displacement of mesh in these models is based on the velocity field calculated with Navier-Stokes equation strongly coupled with heat equation [START_REF] Bruyere | Comparison between Phase Field and ALE methods to model the keyhole digging during spot laser welding[END_REF][START_REF] Tomashchuk | Multiphysical modelling of keyhole formation during dissimilar laser welding[END_REF][START_REF] Sharma | A study of hydrodynamics of melt expulsion in pulsed Nd:YAG laser drilling of titanium[END_REF][START_REF] Morville | 2D axial-symmetric model for fluid flow and heat transfer in the melting and resolidification of a vertical cylinder[END_REF] or on the propertydepending condition [START_REF] Frei | Modeling Thermal Ablation for Material Removal[END_REF] . 2D axisymmetric models with µ-metric mesh size allowed an accurate description of the keyhole formation [START_REF] Bruyere | Comparison between Phase Field and ALE methods to model the keyhole digging during spot laser welding[END_REF][START_REF] Sharma | A study of hydrodynamics of melt expulsion in pulsed Nd:YAG laser drilling of titanium[END_REF] . The main limitation of ALE consists in the impossibility to represent the complex modifications of topology, like the formation of bubbles or porosities; however, it provided quite accurate description of free surface movement and good agreement with experimental dimensions of the impact zone. The liquid phase is modelled as an incompressible [START_REF] Tomashchuk | Multiphysical modelling of keyhole formation during dissimilar laser welding[END_REF][START_REF] Sharma | A study of hydrodynamics of melt expulsion in pulsed Nd:YAG laser drilling of titanium[END_REF] or weakly compressible [START_REF] Bruyere | Comparison between Phase Field and ALE methods to model the keyhole digging during spot laser welding[END_REF] Newtonian liquid that undergoes laminar flow under the influence of Marangoni convection, natural convection and recoil pressure in the keyhole, produced by the metallic vapor. The recoil pressure that is represented by Clausius-Clapeyron equation or in form of adjusted spatial function is an important factor that determines the depth of the keyhole. To suppress the parasite displacement in the solid part of modelled domain, two approaches exist: to stop the displacement by applying an important frictional dissipation depending on the local fraction of liquid phase (Carman-Kozeny apprroximation [START_REF] Frei | Modeling Thermal Ablation for Material Removal[END_REF] ) or by introducing temperaturedependent fictive viscosity (for example 100 Pa•s [START_REF] Sharma | A study of hydrodynamics of melt expulsion in pulsed Nd:YAG laser drilling of titanium[END_REF] or 1000 Pa•s 2 ). The previous models differ in a way to introduce the energy distribution and adsorption coefficient, which is not perfectly known and depends not only on the aggregate state of irradiated matter, but also on local free surface curvature. Another difficulty consists in the lack of data about exact values recoil pressure due to the partial condensation of metallic vapor in the keyhole and about the variation of materials properties in the liquid phase. The motivation for the present work was determined by the interest to represent both keyhole formation and collapse in a fully auto-consistent 3D model, for further application to multimaterial welding. The modelling in 3D accentuated the convergence problems, mainly because of the limitation of minimal mesh size. Our previous model of dissimilar welding 2 suffered from noticeable lack of weld penetration, compared with experimental results. Moreover, the influence of many numerical (inconsistent stabilization for Navier-Stokes equation, type of ALE smoothing, type of used Navier-Stokes equation) and physical parameters (dynamic viscosity, adsorption coefficient, and condensation coefficient) on the convergence and on the resulting behavior of the impact zone was not fully comprehended. The present work attempts to bring some clarity about the influence of these factors on the simple example of standalone laser pulse applied to the plate of Ti6Al4V alloy.

Governing Equations

Heat transfer

Heat transfer equation was used in time-dependent form:

  T T U t T c eq p                     (1)
Laser energy supply is represented by Gauss heat source applied at the top of butt joint. Pulsed laser beam is represented as follows:

  pulse r y x L L t t e r A P q              2 0 2 2 2 0  , ( 2
)
where PLlaser power, r0beam radius and t pulsepulse duration A -absorption coefficient of laser radiation by each individual material, represented as :

A = Asolid +(Aliquid -Asolid)  flc2hs(T-Tm, ∆T), (3)
where Asolid = 0.4 and the absorption coefficient in liquid phase increases when the surface temperature nears vaporization temperature:

Aliquid = Asurf+(Akh-Asurf) flc2hs(T-Tv, ∆T), (4)
where Asurf = 0.25 and Akh = 0.8.

Latent heat of fusion and evaporation are taken in account by means of equivalent enthalpy approach [START_REF] Bonacina | Numerical solutions of phase change problems[END_REF] :

m m p eq p L D c c .   ( 5
)
where Cp is heat capacity as function of temperature, Lflatent heat of fusion and Dm -Gauss function normalized around melting temperature Tm:

  2 2 2 T e D T T T m m       ( 6
)
where ∆T is smoothing interval of 100 K.

Materials properties of Ti6Al4V are considered as interpolation functions of temperature 7 .

Fluid flow

Liquid metal is assumed as incompressible Newtonian liquid that undergoes laminar flow:

    F U U T pI U U t U t                                          ) (   0    U  (7)
Metal that has T<Tm is considered as highly viscous fluid, that practically stops all convection movements in solid materials. The transition between solid and liquid material is provided with smoothed where the condensation coefficient β = 0.5 and the coefficients a,b, c for Ti6Al4V are given by Kaidalov 8 .

Free surface

Present study uses ALE moving mesh approach for representing the movements of free surface. The advantage of this method consists in possibility to neglect the vapor phase and thus lighten the calculation. The weak point consists in limitations of free surface movement: the formation of droplets or bubbles in the melted zone is forbidden by the impossibility to break or unite individual meshes.

Movement of free surface is determined by the velocity field in the melted zone.

Numerical implementation

Present model involves

Experimental study

Welding experiments were carried out with pulsed Nd:YAG laser of 3 kW maximal power and focused beam diameter () of 560 µm.

Laser beam was focused on the surface of butt joint between 2 mm thick Ti6Al4V plates. Standalone pulses with duration of 4 and 6 ms and power 1.5 kW were realized. Three standalone pulses per condition were made.

The penetration of pulses at the joint line was measured by breaking the weld and observing Ti6Al4V side of joint plane with optical microscope.

Discussion

Comparison between the different methods of mesh smoothing

The main challenge in use of ALE for modelling of keyhole welding consists in rapid surface deformation, which is associated with quite severe elongation of meshes. In the first place, the performances of different smoothing methods were compared (Figure 1, Table 1). The Laplace smoothing (Figure 1.a) led to the rapid divergence after interaction time of only 1.5 ms, associated with quite high distortion of the elements. The Yeoh smoothing (Figure 1.b) allowed to extend the calculation up to only 3.2 ms, followed by divergence, in spite of low estimated distortion of the elements. Parasite ripple-like surface deformation was observed in the solid around the melted zone. The Winslow smoothing (Figure 1.c) allowed to reach the end of laser pulse (6 ms), however the aspect of free surface is very rough and the width of the keyhole is highly overestimated. Finally, the use of hyperelastic smoothing (Figure 1.d) allowed the most accurate keyhole representation, associated with quite low distortion of the elements. This formulation allows attaining the maximal deformation of free surface, and thus is the most appropriate for this kind of simulation, in accordance with conclusion of Bruyere et al. [START_REF] Bruyere | Comparison between Phase Field and ALE methods to model the keyhole digging during spot laser welding[END_REF] However, the important deformation of the meshes situated at the center of laser pulse leads to the loss of mesh quality at the bottom of the keyhole. These preliminary calculations were performed with inconsistent Navier-Stokes stabilization id = 1 and mean viscosity of solid of 100 Pas.

The influence of numerical stabilization of Navier-Stokes equation

In absence of inconsistent Navier-Stokes stabilization, the model is likely to diverge at the moments corresponding to sudden rise of velocity field, as, for example, during the collapse of the keyhole at the end of the pulse. In general way, the increase of id value promotes the smoothing of velocity field and better convergence, but the excessive values can alter the observed phenomenology. The values of id of 1 (default) and 0.5 were compared (Table 2). The calculations were performed with mean viscosity of solid of 100 Pas. It can be concluded that melt depth is not at all sensitive to the reduction of id value, and makes small effect on melted zone width (W = -124 µm), maximal melt temperature (T = -16K) and velocity (U = +0.2 m/s). The default value of id=1 was conserved for the rest of the calculations.

The effect of viscosity formulations

Another important numerical problem is the parasite velocity field in the solid material. It was observed that the mean viscosity of solid of 100 Pas produces unphysical relaxation of already solidified matter in the ring formed around the zone of the impact (Figure 2). Cutting off the velocity field in the solid area using temperature or viscosity-based condition in ALE module is not enough, because, as it was found, the fictive viscosity of the solid makes important influence on interaction phenomenology. An insufficient value of solid viscosity makes an important influence on temperature evolution of the melt and fluid dynamics. The evolution of the impact zone under standalone laser pulse can be divided in four major stages: the initiation of the keyhole (Figure 3.a), the progression of the keyhole (Figure 3.b), the keyhole collapse (Figure 3.c) and the solidification resulting in final surface profile with the ring around the impact zone and slight depression in its center (Figure 3.d).

The main characteristics of the process and the evolution of maximal temperature, velocity field and keyhole depth for solid viscosity values of 100 and 200 Pas are illustrated in Table 3 and Figure 4 respectively. The increase of solid viscosity has important impact on calculation time. Moreover, higher values of solid viscosity (300 Pa•s) leads to the convergence problems. Two calculations produce the similar evolution of the maximal temperature on the stage of keyhole progression and collapse, however, higher solid viscosity results in little higher speed of collapsing liquid (2.45 m/s compared to 1.55 m/s). average experimental dimensions of the melted zone speaks in advantage of higher solid viscosity (Table 4).

For both solid values, the depth of the melt is underestimated, when the width is correct, however higher solid viscosity results in deeper melt. The experimental weld was produced on two distinct plates, when numerically only the condition of symmetry is applied on the joint line, which can partially explain such high difference between calculated and experimental values. The simplified representation of absorption in the keyhole and the lack of information about recombination coefficient β influencing the recoil pressure (equation 9), along with the difficulties of accurate numerical resolution of keyhole bottom, are likely to be responsible for underestimation of melt depth. The increase of solid viscosity has an immediate effect on calculation time, multiplied roughly by 16 (Table 3). The attempt to apply the gradual decrease of fictive viscosity with temperature (Heaviside function, equation 10) can be used to stabilize the numerical solution and produces the behavior of the melt similar to the use of high (200 Pas) solid viscosity, however, very important calculation time makes this approach unattractive. The effect of the viscosity in the liquid phase also was analyzed: the temperature-dependent formulation (equation 11 considered by default in all calculations) was compared with constant value if T>Tm (Table 5). As the viscosity of liquid decreases with temperature, in the first case the digging of the melted zone is facilitated and consequently the higher depth of the melt is obtained. However, no important effect was noticed on the maximal melt velocity or temperature, as well as on weld width. The behavior of the liquid zone is quite similar (Figure 5).

          K K l 500 2 T if 10 1.5 2500 T T if 0.0118 + T 10 4 - 3 - m -6  (11)

The influence of keyhole absorption coefficient

The absorption coefficient in the keyhole makes a strong influence on maximal surface temperature, and thus on the resulting recoil pressure (equation 9). The parametric study for Akh values from 0.6 to 0.9 was carried out. In these calculations, dynamic viscosity was considered as function of temperature (equation 11) for the liquid and 200 Pas for the solid. It was found that the shift from 0.7 to 0.8 produces an important shift in keyhole depth, associated with rise of maximal temperature and velocity field (Table 6).

For Akh  0.7, the first peak in velocity field associated with rapid digging of the keyhole is absent, and the velocities associated with the collapse of the keyhole are also lower (Figure 6). To reach the realistic values of penetration, it is needed to use high adsorption coefficients 0.9-1, but it this cases the convergence was not reached for the moment. No convergence

The influence of condensation coefficient

The condensation coefficient represents the variation of surface pressure in function of evaporation conditions (equation 9): from  = 0.17 for strong evaporation to  = 1 for full thermodynamic equilibrium [START_REF] Fabbro | Quelques processus physiques présents en régime de soudage par keyhole[END_REF] . Strong evaporation diminishes the surface pressure (Figure 7), which can influence the keyhole formation. By default, the value of  = 0.5 is admitted in our calculations. Lower values of  resulted in convergence problems. The parametric study was performed for  = 0.5 -1.

It was found that for  > 0.5 the intense initiation of keyhole formation is absent, which traduces by lower weld penetration, temperature and melt velocity (Table 7). Melted zone width is not very sensible to the variation of . For  = 0.6 -1, quite similar depths of the melted zone, maximal temperatures and maximal velocities of the melt (associated with keyhole collapse) are observed. It looks like the increase of  (and so the recoil pressure) slows down the digging of the keyhole, which seems unphysical. The observation of relative pressures at weld cross-section at the end of the pulse (Figure 7.b) shows the zones of maximal pressure at the bottom of the keyhole, which seems physical, but also around the ring of ejected matter, where the displacement of the meshes is blocked by the absence of velocity field. This second maximum is a numerical artifact with no physical significance, as well as the zones of low pressure situated at the solid material beneath. This local increase of pressure may impede the ejection of liquid from the keyhole, which explains why low melt depth for  > 0.5.

Comparison of incompressible and weakly compressible Navier-Stokes equations

Axisymmetric ALE-based model of keyhole digging proposed by Bruyere et al. [START_REF] Bruyere | Comparison between Phase Field and ALE methods to model the keyhole digging during spot laser welding[END_REF] showed slight loss of mass of the domain. The same problem was noticed in the present 3D model. In the calculations performed with incompressible Navier-Stokes formulation, the noticeable loss of mass produces during the creation of the keyhole (Figure 8). It is proportional to the maximal deformation of the domain and took place for both considered values of solid viscosity. After the collapse of the keyhole, the relaxation of strained meshes allows to "recover" some missing mass, but not completely, which results in the situation when the depression in the center of solidified melt is not fully compensated by the height of solidified ring around the melt. Weakly compressible Navier-Stokes equation offers good mass conservation (Figure 8), comparable with result of Bruyere et al. [START_REF] Bruyere | Comparison between Phase Field and ALE methods to model the keyhole digging during spot laser welding[END_REF] , however, the convergence of the models is more difficult due to higher temperature and velocity gradients and more rapid deterioration of element quality. µm. For higher solid viscosity, the melt becomes  130 µm less profound, and there is no impact on residual depression.

  Heaviside function:(T) = solid +(liquid-solid) flc2hs(T-Tm, ∆T).



  Marangoni convection with  <0 introduced in form of weak formulation,  recoil pressure of vapor filling the keyhole, represented in form of function

  Heat Transfer in Fluids, Laminar flow and Moving Mesh modules that are strongly coupled and solved with time-depended solver. Multiphysics coupling is ensured by the interrelation between heat transfer and Navier-Stokes equations: velocity field in heat equation is provided by resolution of Navier-Stokes equation, when temperature field defines thermophysical properties of the liquid involved in convection process. Consistent (streamline and crosswind) stabilization and inconsistent stabilization (isotropic diffusion) were applied to Navier-Stokes equation. Consistent stabilization was also applied to the heat transfer equation. Moving Mesh solves surface deformation in function of locally calculated velocity field. Only one plate (3mm x 2mm x 2mm) is modelled for symmetry reasons. Tetrahedral mesh of 250 µm maximal size is applied over calculation domain. The top surface of is meshed with 50 µm. The models were solved on work station with Intel 2.20 (2 processors, 44 cores) and 256 Gb RAM.

Figure 1 .

 1 Figure 1. The comparison of ALE smoothing methods applied to the model of pulsed welding of Ti6Al4V (tpulse = 6 ms, =560 µm, PL= 1.5 kW).

Figure 2 .

 2 Figure 2. Unphysical relaxation of solidified weld between t = 14 ms (end of solidification) and t = 20 ms (wireframe) for solid viscosity of (a)100 Pas and (b) 200 Pas.The main difference between the cases with solid of 100 Pas and of 200 Pas consists in the intensity of keyhole initiation process. The application of high viscosity of the solid impedes the "parasite" lateral evacuation of the heat from impact zone, because the residual convection in solid material is efficiently reduced. This results in sudden progression of the keyhole depth starting from 1.5 ms with maximal temperature of 4480 K at the center of the impact and the velocity of the liquid reaching 7 m/s. The solid viscosity of 100 Pas does not produce such peak of temperature and velocity: during all keyhole life they remain at the plateau of 3380 K and 1.6 m/s. The sudden progression of the keyhole for the case of 200 Pas results in 80 µm more profound keyhole.

Figure 3 .

 3 Figure 3. Nd:YAG laser pulse on Ti-6Al-4V plate (6 ms, 1.5 kW, Ø = 560 µm). solid = 200 Pas.

Figure 4 .

 4 Figure 4. The influence of solid viscosity on the evolution of maximal temperature, melt velocity and keyhole depth for standalone Nd:YAG laser pulse on Ti-6Al-4V plate (1,5 kW, Ø = 560 µm).Higher solid viscosity also produces more tardive completion of the solidification of the melt: 15.5 ms compared to 14 ms for 100 Pas. It also modifies the proportion of the melt, making it more profound and less large (Table3). In spite of quite similar keyhole depth (530 µm), the case with 200 Pas results in 217 µm deeper melt, when the thickness of the melted zone is reduced at 130 µm. The comparison with the

  solid = cold +(hot-cold) flc2hs(T-900, ∆T) (10)

Figure 5 .

 5 Figure 5. The influence of liquid viscosity on the evolution of maximal melt velocity for standalone Nd:YAG laser pulse on Ti-6Al-4V plate (PL= 1.5 kW, Ø = 560 µm).

Figure 6 .

 6 Figure 6. The influence of keyhole absorption coefficient on the evolution of melt velocity for standalone Nd:YAG laser pulse on Ti-6Al-4V plate (1.5 kW, Ø = 560 µm).

Figure 7 .

 7 Figure 7. The effect of condensation coefficient (a) on the rise of recoil pressure with temperature and (b) on the relative pressure field at the end of the pulse ( t= 6 ms).

  For instance, only the association of compressible Navier-Stokes formulation with Heaviside function for solid viscosity (equation 10) allowed to reach the convergence along with good mass conservation. The different values of solid viscosity produce different effect in incompressible and weakly compressible Navier-Stokes equations (Table8). For solid viscosity of 100 Pas, the reduction of penetration depth is negligible and the maximal residual depression in the impact center increased by only 20

Figure 8 .

 8 Figure 8. Mass conservation compared for incompressible and compressible Navier-Stokes formulations.

Table 1 :

 1 The comparison of ALE smoothing methods

	Smoothing	Convergence	Element distortion	zKH a (µm)
	Laplace	No	12.84	120
	Yeoh	No	1.10	301
	Winslow	Yes	7.20	477
	Hyperelastic	Yes	3.2	532
	a maximal keyhole depth

Table 2 :

 2 The effect of inconsistent stabilization

	id	MZ (µm) Depth Width	Tmax (K)	Umax (m/s)
	1	611	1454	3380	1.6
	0.5	611	1330	3364	1.8

Table 3 :

 3 The effect of mean viscosity of the solid

	solid	MZ (µm)	Tmax	Umax	Calculation time (h)
	(Pas)	Depth Width	(K)	(m/s)
	100	611	1454	3380 1.6	4h57
	200	828	1324	4480 7.5	72h24
	Hf a	698	1450	4870 10.0	144h52
	a Heaviside function , equation 10

Table 4 :

 4 The comparison with experimental weld dimensions Experimental dimensions of the melted zone (µm)

	Test 1	1055	1286	1269	1553
	Test 2	1059	1349	1259	1426
	Test 3	1008	1332	1222	1379
	Average	1041	1322	1250	1453
	solid (Pas)	Relative error (%) after 4 ms after 6 ms Depth Width Width Depth
	100	57.0	10.2	8.0	51.0
	200	47.0	9.6	9.0	34.0
	Hf a	48.6	3.2	0.3	44.2
	a Heaviside function

Table 5 :

 5 The effect of viscosity of the liquid

		Melted zone			Calculation
	liquid	dimensions	Tmax	Umax	time (h)
	(Pas)	(µm) Depth Width	(K)	(m/s)
	f(T) a	828	1324	4480 7.5	72h24
	410 -3	725	1358	4400 7.5	82h23
	a equation 11

Table 6 :

 6 The effect of keyhole absorption coefficient

	Akh	T max (K)	U max (m/s)	MZ (µm) Depth Width
	0.6	3445	2.28	587	1308
	0.7	3460	2.28	587	1308
	0.8	4480	5.5	828	1324
	0.9

Table 7 :

 7 The effect of condensation coefficient

		MZ (µm)	Tmax (K)	Umax (m/s)
		Depth	Width
	0 No convergence
	0.4 No convergence
	0.5	828	1324	4480	7.5
	0.6	589	1308	3440	2.34
	0.7	594	1312	3420	2.18
	0.8	597	1308	3400	2.04
	0.9	600	1314	3390	2.04
	1	600	1318	3380	2.04

Table 8 :

 8 The effect of Navier-Stokes formulation Conclusions 1. Hyperelastic deformation seems the most appropriate for the modelling keyhole evolution. 2. The velocity field not sensible to the inconsistent stabilization coefficient for Navier-Stokes at the range 0.5-1. 3. Low viscosity of the solid phase not only produces parasite relaxation of solidified domain, but can completely alter the phenomenology of laser/metal interaction. Thus at least 200 Pa•s value is recommended. The solution is much less sensible to the formulation of viscosity in the liquid phase. However, temperature-dependent viscosity promotes the digging of the keyhole. 4. Keyhole adsorption coefficient is the key parameter that controls surface temperature and melt depth. The shift from 70 to 80 adsorption produces the qualitative change in the dynamics of keyhole progression, with much more intense digging in the beginning of the pulse. Higher value would create even more profound weld, and thus reach experimentally observed penetration, but it is associated with too severe mesh deformation producing the convergence problems. This problem can be treated with different remesing approaches that increase the time of calculation to the high degree. 5. The condensation coefficient of 0.5 seems the most appropriate. Higher values impede the digging of the keyhole because unphysical accumulation of pressure around the ejection ring, and lower values produce convergence problems. 6. The shift from incompressible or weakly compressible Navier-Stokes formulation does not alter the observed phenomenology and allows better mass conservation. However, it is likely to be more sensible to the convergence problems. Agency of Research.

		0.6
		0.5
		0.4
	Mass	0.2 0.3
		0.1
		0
		0	0.005	0.01	0.015	0.02
		200 Pa*s Incompressible
		100 Pa*s Incompressible
		200 Pa*s, Compressible
		100 Pa*s Compressible
		Hf, Compressible
			solid	Controlled dimensions (µm)
	Navier-Stokes	(Pas)	Meted zone depth	Impact center depression
	Incompressible 100	611	125
			200	828	150
	Compressible	100	596	144
			Hf a	698	150
	a Heaviside function

Acknowledgements

This work was carried out as a part of joint laboratory project LabCom FLAMme between Laboratoire Interdisciplinaire Carnot de Bourgogne, University of Bourgogne-Franche Comté and SME Laser Rhone-Alpes. This project is funded by French National


Armaël Guéneau

Procrastination A proof engineering technique

We present a small Coq library for collecting side conditions and deferring their proof.

INTRODUCTION

What are the different ways of proving an assertion of the form "∃ x, P x"? From the perspective of the final proof term, there is only one: providing a witness. However, from a "proof engineering" perspective, Coq enables a wider range of possible approaches thanks to its powerful feature of "evars" [5, §2.11].

We describe existing approaches for proving such a goal, and propose a new one, supported by a small library dubbed Procrastination [START_REF] Guéneau | The procrastination library[END_REF], which itself relies on evars. We believe it captures a useful proof pattern, where one wishes to first continue with the proof of P x, in order to discover which constraints x must satisfy, before providing a value for x and proving the side conditions.

We used this library with success in our work on verifying the asymptotic complexity of programs [START_REF] Guéneau | A fistful of dollars: Formalizing asymptotic complexity claims via deductive program verification[END_REF].

PROVING EXISTENTIAL QUANTIFICATIONS IN COQ

Manually guessing a witness. The most straightforward option to prove a Coq goal "∃ x, P x" is to guess and provide upfront a witness, using the exists tactic. This is often a good way of providing external knowledge about how the proof should continue.

Delaying the instantiation using evars. However this falls short if the only intuition is that x should be "whatever makes the proof work". A trial-and-error process can be attempted, where the user alternatively adjusts the witness and the proof until it passes. This does not lead to maintainable proof scripts: in the end, only a "magic" witness remains; it may be large and duplicate information; it may also have to change in the future, following updates to the rest of the proof.

To allow delaying the instantiation of x, Coq provides a versatile feature: evars. An evar corresponds to a "hole" in the proof, and only exists during the pre-typing phase -it will not appear in the proof term, and only helps elaborating it. Using the eexists tactic on the goal "∃ x, P x" turns it into "P ?x", introducing a fresh evar ?x. Later, it can be discovered that this hole should be (definitionally) equal to some term: the evar then gets instantiated. This is generally performed automatically by Coq tactics through unification.

"Big enough" existentials. This technique again falls short if continuing with the proof of P ?x does not exhibit a term that ?x must be equal to, but e.g. lower bounds on ?x. Cohen [1, §5. 1.4] shows that it is a common situation in analysis, where a proof is established by picking a "big enough x". During the proof, inequalities of the form "a ≤ x", "b ≤ x", ... are all trivially satisfied provided that "x is big enough".

Cohen formalizes this seemingly fuzzy reasoning in a small tactic library. On a "∃ x, P x" goal, a tactic instantiates x with the iterated maximum of a list of terms, this list being initially an evar. It is then possible to progressively populate this list, automatically discharging subgoals of the form "a ≤ x".

THE PROCRASTINATION LIBRARY

We are interested in the situation where arbitrary side conditions about x may be encountered -e.g. inequalities of both the form "a ≤ x" and "x ≤ b". We might also be trying to infer not an arithmetic value but a function, e.g. f satisfying "∀ n, f(n)≥ 1 + f(n-1)". In this situation, it seems difficult to guess the shape of x or f beforehand (e.g. as an iterated maximum).

We introduce a small Coq library named Procrastination, which generalizes the ideas of Cohen, and allows collecting arbitrary side conditions (to be proved later) about zero, one or several variables. It is more general than "big enough" in the sense that it can deal with arbitrary side conditions, but it does not try to solve these side conditions: this is typically done in a second time by a dedicated procedure (e.g. for "big enough" this procedure would take the running maximum of the lower-bounds collected).

In other words, the library enforces a pattern where the proof is separated in two phases: first side conditions are collected about some variables, then they are solved and the variables are instantiated. This is crucial in limiting the proliferation of evars, which tend to produce brittle proofs because they can might get incorrectly instantiated by unrelated tactics.

The implementation provides reusable Ltac components, with the hope that they can help applying this technique even in more exotic scenarios. The code is freely available online [START_REF] Guéneau | The procrastination library[END_REF].

EXAMPLE

Our initial motivation for this work was analyzing the time complexity of programs, where one has to find a solution for a set of inequations about a program's execution time. One approach is using Cormen et al.'s substitution method [2, §4]. The idea is to guess a parameterized shape for the solution; substitute this shape into the goal; gather a set of constraints the parameters must satisfy for the goal to hold; finally, show that these constraints are indeed satisfiable.

We demonstrate this strategy using Procrastination on cost inequations derived from a binary search program. We "guess" that the solution should be a monotonic function f with shape λn. if n = 0 then c else (a log n + b); where a, b and c are parameters, about which we wish to gradually accumulate a set of constraints.

Goal ∃ (f: nat → nat), 1 ≤ f 0 ∧ monotonic f ∧ ∀ n, 0 < n → 2 + f (n/2) ≤ f n.
Proof. begin procrastination assuming a b c.

∃ f: nat → nat, 1 ≤ f 0 ∧ monotonic f ∧ ∀ n, 0 < n → 2 + f (n/2) ≤ f n (1)
a, b, c: nat g: Group (?G a b c)

1 ≤ c (2) a, b, c: nat g: Group (1 ≤ c ∧ ?Goal0) Subgoal solved (3)
∃ a b c, 1 ≤ c ∧ c ≤ b ∧ 2 + c ≤ b ∧ 2 ≤ a (4)
The script starts with begin procrastination assuming a b c which introduces the parameters, as well as a "group" д, using an evar to collect the constraints [START_REF] Cohen | Formalized algebraic numbers: construction and first-order theory[END_REF]. Note that a, b and c are not evars themselves, instead they appear as abstract variables. Note also that the type of these variables is inferred to be nat only further on in the proof script. After instantiating the function with the proposed solution, the first inequality simplifies to a condition on c (2). This condition is stored (in д) for later retrieval, using the procrastinate tactic, which discharges the subgoal (3). Similarly, the other two inequalities yield side conditions about a, b and c that can be procrastinated. Finally, end procrastination retrieves the collected conditions, which the user must then prove. Here, this is done by an ad hoc user-defined tactic, eomega.

Observe that the final proof script never mentions the actual witnesses: we expect the proof to be robust in the face of minor changes in the cost inequations, as long as the shape of the solution remains valid.

( 1 )

 1 exists (fun n ⇒ if zerop n then c else a * log2 n + b). repeat split. { simpl. (2) procrastinate. (3) } { ... (* c ≤ b *) procrastinate. ... } { intros. ... (* 2 + c ≤ b *) procrastinate. ... (* 2 ≤ a *) procrastinate. } end procrastination. (4) eomega. Qed. a: ?x b: ?x0 c: ?x1 g: Group (?G a b c)


Donguk Kim
email: donguk319@snu.ac.kr

Junesol Song
email: junesol.song@recherche.enac.fr

Sunkyoung Yu

Changdon Kee
email: kee@snu.ac.kr

Moonbeom Heo

A New Algorithm for High-Integrity Detection and Compensation of Dual-Frequency Cycle Slip Under Severe Ionospheric Storm Conditions

Keywords: cycle-slip detection, cycle-slip compensation, insensitive cycle-slip pairs, high-integrity detection, realtime kinematic (RTK)

Many strategies for treating dual-frequency cycle slip, which can seriously affect the performance of a carrier-phase-based positioning system, have been studied over the years. However, the legacy method using the Melbourne-Wübbena (MW) combination and ionosphere combination is vulnerable to pseudorange multipath effects and high ionospheric storms. In this paper, we propose a robust algorithm to detect and repair dual-frequency cycle slip for the network-based real-time kinematic (RTK) system which generates high-precision corrections for users. Two independent and complementary carrier-phase combinations, called the ionospheric negative and positive combinations in this paper, are employed for avoiding insensitive pairs. In addition, they are treated as second-order time differences to reduce the impact of ionospheric delay even under severe ionospheric storm. We verified that the actual error distributions of these monitoring values can be sufficiently bounded by the normal Gaussian distribution. Consequently, we demonstrated that the proposed method ensures high-integrity performance with a maximum probability of missed detection of 7.5 × 10 -9 under a desired false-alarm probability of 10 -5 . Furthermore, we introduce a LAMBDA-based cycle slip compensation method, which has a failure rate of 1.4 × 10 -8 . Through an algorithm verification test using data collected under a severe ionospheric storm, we confirmed that artificially inserted cycle slips are successfully detected and compensated for. Thus, the proposed method is confirmed to be effective for handling dual-frequency cycle slips of the network RTK system.

Introduction

Recently, the demand for high-precision navigation systems employing carrier-phase observations has been growing rapidly for applications such as automated vehicle driving and monitoring, collision avoidance, and intelligent transportation systems [START_REF] Song | A Study on GPS/GLONASS Compact Network RTK and Analysis on Temporal Variations of Carrier Phase Corrections for Reducing Broadcast Bandwidth[END_REF]. One of the typical techniques to obtain cm-level position accuracy is the real-time kinematic (RTK) technique, which is widely used for geodesy and surveying, however, the RTK has been constrained to short-baselines under 10 km. Over the past few decades, network-based RTK techniques have been studied to enlarge RTK coverage per station from medium to long-baseline, for use by dynamic users. The compact network RTK method developed by GNSS Laboratory at Seoul National University is considered as a candidate solution for land vehicle navigation because it provides cm-level positioning service with fast ambiguity resolution under an extremely low-rate data link [START_REF] Song | A Study on GPS/GLONASS Compact Network RTK and Analysis on Temporal Variations of Carrier Phase Corrections for Reducing Broadcast Bandwidth[END_REF][START_REF] Park | A Study on Reducing Temporal and Spatial Decorrelation Effect in GNSS Augmentation System: Consideration of the Correction Message Standardization[END_REF][START_REF] Park | The Compact Network RTK Method: An Effective Solution to Reduce GNSS Temporal and Spatial Decorrelation Error[END_REF][START_REF] Song | A Study on Improving Performance of Network RTK through Tropospheric Modeling for Land Vehicle Applications[END_REF].

One of the most important issues in high-precision navigation systems including the network RTK is generating the integrity information for ensuring the safety of life. The safety-critical navigation system for aviation applications such as satellite-based augmentation system (SBAS) provides integrity information of its corrections [START_REF] Pullen | Augmented GNSS: Fundamentals and Keys to Integrity and Continuity[END_REF], however, the carrier-phase-based system, which provides a high-precision correction for users, has a challenging problem that integer ambiguity must be correctly resolved and validated [START_REF] Teunissen | GNSS Integer Ambiguity Validation : Overview of Theory and Methods[END_REF]. Furthermore, cycle slip, which is an instantaneous jump of an integer number of a cycle, is a major integrity threat for carrier-phase observations, even after the ambiguity is correctly fixed and validated. A cycle slip occurs unexpectedly when the receiver's phase-locked loop (PLL) has a loss of lock during a temporary signal blockage or an ionospheric scintillation [START_REF] Bisnath | A New Approach to an Old Problem: Carrier-Phase Cycle Slips[END_REF][START_REF] Seo | Future Dual-Frequency GPS Navigation System for Intelligent Air Transportation Under Strong Ionospheric Scintillation[END_REF][START_REF] Luo | A Study of Multi-GNSS Ionospheric Scintillation and Cycle-Slip over Hong Kong Region for Moderate Solar Flux Conditions[END_REF]. The cycle slip must be handled at the data pre-processing stage since it induces an error with an unpredictable range, which can seriously affect the quality of high-precision corrections and user position solution [START_REF] Bisnath | A New Approach to an Old Problem: Carrier-Phase Cycle Slips[END_REF].

A number of processing methods for cycle slip have been studied over the years. In particular, the cycle-slip detection and repair method using dual-frequency observations was proposed by Blewitt [START_REF] Blewitt | An Automatic Editing Algorithm for GPS data[END_REF], Gao and Li [START_REF] Gao | Cycle Slip Detection and Ambiguity Resolution Algorithms for Dual-Frequency GPS Data Processing[END_REF], Bisnath et al. [START_REF] Bisnath | A New Approach to an Old Problem: Carrier-Phase Cycle Slips[END_REF], Liu [START_REF] Liu | A New Automated Cycle Slip Detection and Repair Method for a Single Dual-Frequency GPS Receiver[END_REF], and Cai et al. [START_REF] Cai | Cycle Slip Detection and Repair for Undifferenced GPS Observations Under High Ionospheric Activity[END_REF]. These popular methods employ two complementary geometry-free linear combinations for dynamic users, the Melbourne-Wübbena (MW) combination and ionosphere combination, since the combined cycle slip can be close to zero and cannot be detected using only one combination. For example, the same size of cycle slips on L1 and L2 frequencies such as [START_REF] Blewitt | An Automatic Editing Algorithm for GPS data[END_REF][START_REF] Blewitt | An Automatic Editing Algorithm for GPS data[END_REF] cannot be noticed on the MW combination, but it has a detectable value on the ionosphere combination. Similarly, the special cycle-slip pair (77, 60), which does not influence the ionosphere combination, is detected by the MW combination. We call these special cycle-slip pairs insensitive pairs [START_REF] Song | A Study on Improving Performance of Network RTK through Tropospheric Modeling for Land Vehicle Applications[END_REF][START_REF] Cai | Cycle Slip Detection and Repair for Undifferenced GPS Observations Under High Ionospheric Activity[END_REF][START_REF] Song | A New Approach for GNSS Frequency Selection Considering Detection of Cycle Slip Insensitive Pairs of Ionospheric Combination for Dual-Frequency Receivers[END_REF].

Most dual-frequency insensitive cycle-slip pairs can be detected by the MW combination and ionosphere combination. In addition, this detection algorithm can be easily used for static and even dynamic users because of the benefit of geometry-free combination; however, it has critical limitations when detecting small cycle slips. First, these techniques are affected by noise and pseudorange multipath effects even if a smoothing or averaging technique such as a Hatch filter and low-pass filter is applied. In other words, these techniques cannot be processed instantaneously and might fail to detect small cycle slips under high-multipath conditions [START_REF] Cai | Cycle Slip Detection and Repair for Undifferenced GPS Observations Under High Ionospheric Activity[END_REF][START_REF] Song | A New Approach for GNSS Frequency Selection Considering Detection of Cycle Slip Insensitive Pairs of Ionospheric Combination for Dual-Frequency Receivers[END_REF]. Furthermore, the ionosphere combination is vulnerable to severe ionospheric storm since the remaining bias of the ionospheric variation makes it difficult to detect small cycle slips [START_REF] Cai | Cycle Slip Detection and Repair for Undifferenced GPS Observations Under High Ionospheric Activity[END_REF][START_REF] Liu | A New Real-Time Cycle Slip Detection and Repair Method Under High Ionospheric Activity for a Triple-frequency GPS/BDS Receiver[END_REF][START_REF] Zeng | Real-Time Triple-Frequency Cycle Slip Detection and Repair Method Under Ionospheric Disturbance Validated with BDS Data[END_REF]. The probability of occurrence of these insensitive pairs is incredibly small [START_REF] Blewitt | An Automatic Editing Algorithm for GPS data[END_REF], but it is still considered an integrity threat to a system requiring high-integrity performance.

To overcome these limitations, triple-frequency signals, which allow many additional combinations such as the extra-wide lane combination and ionospheric reduced combination, are being emphasized for improving performance [START_REF] Zhang | Benefits of the Third Frequency Signal on Cycle Slip Correction[END_REF][START_REF] Dai | Instantaneous Triple-Frequency GPS Cycle-Slip Detection and Repair[END_REF][START_REF] Hu | New-Generation BeiDou (BDS-3) Experimental Satellite Precise Orbit Determination with an Improved Cycle-Slip Detection and Repair Algorithm[END_REF][START_REF] Pu | An Improved Algorithm Based on Combination Observations for Real Time Cycle Slip Processing in Triple Frequency BDS Measurements[END_REF][START_REF] Chang | Adaptive Kalman Filter Based on Variance Component Estimation for the Prediction of Ionospheric Delay in Aiding the Cycle Slip Repair of GNSS Triple-Frequency Signals[END_REF]. Especially, many novel combinations of GPS and BeiDou, which are the currently available triple-frequency signals, have been proposed for cycle slip detection without an insensitive pair even under a high-ionospheric storm [START_REF] Liu | A New Real-Time Cycle Slip Detection and Repair Method Under High Ionospheric Activity for a Triple-frequency GPS/BDS Receiver[END_REF][START_REF] Zeng | Real-Time Triple-Frequency Cycle Slip Detection and Repair Method Under Ionospheric Disturbance Validated with BDS Data[END_REF][START_REF] Gu | Detection and Correction of Cycle Slip in Triple-Frequency GNSS Positioning[END_REF]. Despite the many advantages of triple-frequency signals, most of the receivers and GNSS satellites only provide dual-frequency signals, currently. That is, the demand for a cycle-slip detection method based on dualfrequency signals remains very strong.

Considering the limitations discussed above, we propose a novel algorithm for cycle-slip detection and compensation using only the dual-frequency carrier-phase observations without pseudorange. In order for the network RTK to be used as safety-critical systems in the near future, it is necessary to ensure a high-integrity information as well as a high-precision corrections. As a first step to enhancing the network-based RTK system, we focus on the high-integrity detection algorithm for continuously operating reference stations (CORS). In this paper, we introduce two independent and complementary ionosphere combinations: one is the geometry-free ionosphere combination usually used in cycle-slip detection, and the other is a new geometry-based ionosphere combination to replace the MW combination. We name these combinations the ionospheric negative and positive combination, respectively. Song and Kee demonstrated that there is no additional geometry-free carrier-phase combination to replace the MW combination for dynamic users [START_REF] Song | A New Approach for GNSS Frequency Selection Considering Detection of Cycle Slip Insensitive Pairs of Ionospheric Combination for Dual-Frequency Receivers[END_REF]; however, geometry-based combinations can be employed in static permanent stations for generating reliable high-precision corrections. Because our ionospheric positive combination constructs the combined cycle slips with positive signs, it can detect a small cycle slip more efficiently compared to the popularly used geometry-based ionosphere-free combination [START_REF] Dach | Bernese GNSS Software, version 5.2[END_REF][START_REF] Chen | A Double-Differenced Cycle Slip Detection and Repair Method for GNSS CORS Network[END_REF]. Furthermore, we apply a second-order time-differenced observation (or acceleration) to reduce the influence of ionospheric delay under severe ionospheric storm regardless of baseline lengths [START_REF] Cai | Cycle Slip Detection and Repair for Undifferenced GPS Observations Under High Ionospheric Activity[END_REF][START_REF] Liu | A New Real-Time Cycle Slip Detection and Repair Method Under High Ionospheric Activity for a Triple-frequency GPS/BDS Receiver[END_REF][START_REF] Kee | Quality Control Algorithms on WAAS Wide Area Reference Stations[END_REF].

The remainder of this paper is organized as follows: in Section 2, we first describe the new detection and compensation strategy with two ionospheric acceleration combinations used as monitoring values (MVs). Then, we verify whether the actual error distributions of the MVs could be well assumed as a normal Gaussian distribution during severe ionospheric storms. This analysis of statistical error propagation and bounding is very important because it is impossible to collect an extremely large number of data sets for reliable integrity statistics [START_REF] Pullen | Augmented GNSS: Fundamentals and Keys to Integrity and Continuity[END_REF][START_REF] Decleene | Defining Pseudorange Integrity-Overbounding[END_REF][START_REF] Schempp | An Application of Gaussian Overbounding for the WAAS Fault Free Error Analysis[END_REF]. Finally, we evaluate the probability of missed detection and false alarm, which are the most significant indices of a detection algorithm that indicate the integrity performance [START_REF] Altmayer | Enhancing the integrity of integrated GPS/INS systems by cycle slip detection and correction[END_REF][START_REF] Song | Predictions of Allowable Sensor Error Limit for Cycle-Slip Detection[END_REF][START_REF] Kim | GPS Cycle Slip Detection Considering Satellite Geometry based on TDCP/INS Integrated Navigation[END_REF]. Many former studies overlooked the probability of missed detection and only analyzed whether the size of the cycle slip was greater than the threshold; however, it is necessary to demonstrate the performance of the proposed method quantitatively for conserving integrity. Furthermore, we introduce a reliable cycle-slip compensation method based on the LAMBDA technique [START_REF] Teunissen | GNSS Integer Ambiguity Validation : Overview of Theory and Methods[END_REF][START_REF] Teunissen | Performance of the LAMBDA Method for Fast GPS Ambiguity Resolution[END_REF]. In Section 3, we conduct an algorithm test with actual observations to verify the proposed method under a severe ionospheric storm. We inserted the simulated cycle slips artificially into the raw data and checked whether all of them are successfully detected and compensated for. In Section 4, we present our conclusions.

Cycle-Slip Detection and Compensation Algorithm

Dual-Frequency Cycle-Slip Detection Method

TDSD Carrier-Phase Observations

In order to monitor and detect phase anomalies such as cycle slips and impulse outliers, L1 and L2 carrier-phase observations should be handled as time differences (TDs) and single differences (SDs) between static receivers. By single differencing the observations between receivers, satellite clock biases are eliminated. The TDSD L1 and L2 carrier-phase observations of reference stations are defined as follows:

T T T T T T T T T T CS I L B CS I L B 1 2 1 1 1 2 2 2                                           , (1)
where the symbols  T and  indicate the TD and SD operators, respectively;  is the measured carrier phase; the subscripts "1" and "2" represent the L1 and L2 frequencies, respectively;  is the wavelength; CS is an integer cycle slip that may exist in phase; I is the ionospheric slant delay on the L1 frequency; and  is the square of the ratio of the L1 and L2 frequencies. In this process, we can eliminate the geometric distance using satellite ephemeris since we know the precise position of reference stations. Furthermore, the tropospheric slant delay can be reduced using general model such as the Saastamoinen and Hopfield model [START_REF] Hofmann-Wellenhof | GNSS-Global Navigation Satellite Systems: GPS, GLONASS, Galileo, and More[END_REF]. The symbol L denotes the sum of the geometric distance error and the modeling error of tropospheric slant delay. B is the receiver clock offset, and  is the receiver noise of TDSD observation. The noise level of the TDSD observation     T is twice that of the un-differenced carrier phase, which is typically measured as approximately

      1 2
2 mm at receivers of CORS [START_REF] Park | A Study on Reducing Temporal and Spatial Decorrelation Effect in GNSS Augmentation System: Consideration of the Correction Message Standardization[END_REF][START_REF] Misra | Global Positioning System: Signals, Measurements, and Performance[END_REF][START_REF] Li | GNSS Elevation-Dependent Stochastic Modeling and Its Impacts on the Statistic Testing[END_REF].

The bias components of GNSS error sources, except the receiver clock offset, are known to vary slowly over time. The temporal variabilities of the satellite orbit error and tropospheric slant delay are less than 1 and 0.2 mm/s, respectively. Theoretically, the ionospheric slant delay variation can reached 2 cm/s under ionospheric scintillationfree condition [START_REF] Park | The Compact Network RTK Method: An Effective Solution to Reduce GNSS Temporal and Spatial Decorrelation Error[END_REF][START_REF] Olynik | Temporal Characteristics of GPS Error Sources and Their Impact on Relative Positioning[END_REF][START_REF] Olynik | Temporal Impact of Selected GPS Errors on Point Positioning[END_REF][START_REF] Park | RRC Unnecessary for DGPS Messages[END_REF]. The recent observed ionospheric velocity was less than 0.5 cm/s during a high-ionospheric storm at a low latitude [START_REF] Tiwari | Effect of TEC Variation on GPS Precise Point at Low Latitude[END_REF][START_REF] Zhizhao | Study of the Ionospheric TEC Rate in Hong Kong Region and its GPS/GNSS Application[END_REF]. In other words, the time variation of GNSS error sources is much smaller than the size of cycle slips as long as the sampling time is relatively short (1 s in this paper); however, the change of receiver clock, which is unknown and unbounded, should be eliminated to detect the cycle slips.

Receiver Clock Drift Estimate

The SD drift of receiver clock, which has highly stable oscillator, can be estimated precisely from the average value of each satellite's TDSD ionosphere-free (IF) carrier-phase combination expressed as follows [START_REF] Kim | RTK-GPS Correction Generation Technique for Low-Rate Data-Link[END_REF]:

m j T T IF j B m 1 1         , (2)
T IF T IF T T IF T T a a CS L B
where a a

1 1 2 2 1 2 1 , 1 1                                  , (3)
where the symbols

             . ( 4
)
Assuming the noise levels of the TDSD L1 and L2 carrier phase have the same value (i.e.,

               1 2 T T T
4 mm), the variance of the estimated SD receiver clock drift can be expressed as follows:

                    2 2 2 2 2 1 2 1 ˆT T T IF a a B m m . (5)

Cycle-Slip Detection Using the Ionospheric Acceleration

The TDSD L1 and L2 carrier-phase observations after the receiver clock drift compensation can be expressed as follows:

T T T T T T T T T T T T B CS I L B CS I L 1 2 1 1 1 1 2 2 2 2 ˆˆ                                                      . ( 6
)
We call Equation Error! Reference source not found. the TDSD carrier-phase residual. In this section, we describe two L1/L2 linear combinations of carrier phase for obtaining the MVs of cycle-slip detection. The first observation is the geometry-free combination that only contains ionospheric slant delay, which is widely used to detect dual-frequency cycle slips [START_REF] Bisnath | A New Approach to an Old Problem: Carrier-Phase Cycle Slips[END_REF][START_REF] Blewitt | An Automatic Editing Algorithm for GPS data[END_REF][START_REF] Gao | Cycle Slip Detection and Ambiguity Resolution Algorithms for Dual-Frequency GPS Data Processing[END_REF][START_REF] Liu | A New Automated Cycle Slip Detection and Repair Method for a Single Dual-Frequency GPS Receiver[END_REF][START_REF] Cai | Cycle Slip Detection and Repair for Undifferenced GPS Observations Under High Ionospheric Activity[END_REF][START_REF] Kee | Quality Control Algorithms on WAAS Wide Area Reference Stations[END_REF][START_REF] Kim | Instantaneous Real-Time Cycle-Slip Correction for Quality Control of GPS Carrier-Phase Measurements[END_REF]. In this paper, we named it the ionosphere negative (IN) combination since L1 and L2 cycle slips are combined with negative signs:

  T T T T T I I b b I CS CS where b b 1 1 2 2 1 1 2 2 1 2 1 1 1 1 , 1 1                                           . ( 7
)
In general, the size of the combined cycle slips in the IN combination is much greater than the ionospheric velocity and the observed noise level; however, the size can be almost zero even if the sizes of L1 and L2 cycle slips are not small. These special cycle-slip pairs are defined as insensitive pairs, which cannot be detected. Therefore, we need a complementary linear combination that can detect insensitive cycle-slip pairs in the IN combination [START_REF] Song | A Study on Improving Performance of Network RTK through Tropospheric Modeling for Land Vehicle Applications[END_REF][START_REF] Song | A New Approach for GNSS Frequency Selection Considering Detection of Cycle Slip Insensitive Pairs of Ionospheric Combination for Dual-Frequency Receivers[END_REF].

The second observation we have chosen is another ionosphere combination that couples the L1 and L2 cycle slips with positive signs [START_REF] Dach | Bernese GNSS Software, version 5.2[END_REF]. We named it the ionosphere positive (IP) combination, which can be expressed as follows:

T T T T T T I I b b I L CS CS where b b 1 1 2 2 1 1 2 2 1 2 1 1 1 1 1 2 2 1 1 , 2 2                                                        . ( 8
)
The IP combination is not geometry-free; however, the variations of satellite orbit error and tropospheric slant delay are negligible in case of stationary observations, as described in Section 2.1.1. Since two L1/L2 linear combinations combine cycle slips with different signs and are complementary to each other, all cycle slips can be detected under quiet ionospheric storm, provided the ionospheric velocity is sufficiently small. As discussed in Section 2.1.1, the change of ionospheric slant delay cannot be ignored for detecting small cycle slips, because it can be greater than 2 cm/s under high ionospheric storm. Thus, we employ ionospheric acceleration, which is the second-order derivative of ionospheric slant delay, as MVs. Since the bias of ionospheric acceleration is less than 10 -4 m/s 2 even under high ionospheric disturbance [START_REF] Kang | Analysis on Normal Ionospheric Trend and Detection of Ionospheric Disturbance by Earthquake[END_REF], many studies considered it as the MV for robust cycle-slip detection [START_REF] Cai | Cycle Slip Detection and Repair for Undifferenced GPS Observations Under High Ionospheric Activity[END_REF][START_REF] Liu | A New Real-Time Cycle Slip Detection and Repair Method Under High Ionospheric Activity for a Triple-frequency GPS/BDS Receiver[END_REF][START_REF] Kee | Quality Control Algorithms on WAAS Wide Area Reference Stations[END_REF]. The two linear combinations of the ionospheric acceleration can be expressed as follows:

  T T T I I I CS CS 2 2 2 1 1 2 2 1 1                    , ( 9
) T T T T I I I L CS CS 2 2 2 2 1 1 2 2 1 1 1 1 1 2 2                                     . ( 10
)
Equation Error! Reference source not found. corresponds to the IN combination, and Equation

Error! Reference source not found. corresponds to the IP combination. The symbol  2

T represents the second-order TD operator. These MVs have only receiver noise with negligible biases under the non-cycle-slip hypothesis. Therefore, the cycle slips can be detected when the following condition is satisfied:

    T FA T T FA T MV I T K I or MV I T K I 2 2 max 2 2 max                             , (11)
where


T and  T indicate the thresholds for each MV. The detection thresholds are calculated based on the confidence level FA K and the maximum standard deviation of the MV. We will analyze the noise level of the MV in Section 2.2 and discuss how to determine thresholds for the detection criterion in more detail in Section 2.3.

Error Propagation in Monitoring Values

Theoretical Noise Analysis of Monitoring Values

In this section, we theoretically discuss the noise level in the MVs. We should assume the MVs under the noncycle-slip hypothesis as a normal Gaussian distribution. Unfortunately, the actual data does not have a perfect Gaussian distribution, because of, e.g., the influence of carrier-phase multipath effects, clock estimation errors, and the remaining bias of ionospheric acceleration. Therefore, the standard deviation of the MVs must be carefully determined for sufficiently bounding the tails of the actual error distribution [START_REF] Pullen | Augmented GNSS: Fundamentals and Keys to Integrity and Continuity[END_REF][START_REF] Decleene | Defining Pseudorange Integrity-Overbounding[END_REF][START_REF] Schempp | An Application of Gaussian Overbounding for the WAAS Fault Free Error Analysis[END_REF]. In other words, we should theoretically analyze the statistics of the MVs considering the worst possible scenario.

The general form of the linear combination of the TDSD carrier-phase residual can be expressed by Equation Error! Reference source not found., and the variance of that can be calculated from Equation Error! Reference source not found.. We assume the carrier phase of L1 and L2 have the same noise level:

                                          1 2 ( , ) 1 1 2 2 1 1 2 2 ˆT T T T T T T B B , ( 12
)                                 1 2 2 2 2 2 2 2 ( , ) 1 2 1 2 T T T B . ( 13
)
Assuming the biggest value of the variance of the estimated clock drift from Equation ( 4), the variance of the TDSD carrier-phase residual becomes maximum as follows:

                                                 1 2 2 2 2 2 2 2 2 2 1 2 max ( , ) 1 2 1 2 1 2 , 1 T T T a a func m where m . ( 14
)
Since the carrier-phase measurements can be regarded as uncorrelated in time and each observation, the variance of second-order TD observations is 6 times (i.e.,      ) that of the un-difference observations. Therefore, the maximum variance of the MVs of the IN and IP acceleration can be calculated from Equation Error! Reference source not found. and Equation Error! Reference source not found.:

      T T I func b b func b b 2 2 2 2 2 max 1 2 1 2 , 12 ,                    , (15)
      T T I func b b func b b 2 2 2 2 2 max 1 2 1 2 , 12 ,                    . ( 16
)
Finally, the maximum standard deviation of the MVs is theoretically determined as follows if we set the carrier phase noise as 2 mm:

              2 2 max 7.57 15.1 mm s T MV I (17)               2 2 max 8.52 17.1 mm s T MV I (18)

Actual Error Distribution of Monitoring Values

In this section, we analyze the actual error distribution of MVs by using the observed data under various ionospheric storm conditions. We especially validate that the "fat tails" (i.e., non-Gaussian edge) of the actual error distribution are properly bounded to the theoretical error distribution of the MVs that we determined previous section. That is, we prove the actual error distribution can be replaced the assumption of theoretical normal distribution.

A severe geomagnetic storm occurred on 17 March 2015. During this storm, the geomagnetic Kp index reached 8-, and the Dst index dropped to -223 nT [START_REF] Zhang | Ionospheric Disturbances Following the March 2015 Geomagnetic Storm from GPS Observations in China[END_REF]. Figure 1 shows the time history of the geomagnetic indices Kp and Dst for 16-18 March 2015. The raw data of geomagnetic indices were obtained from the World Data Center (WDC) for Geomagnetism (Kyoto, Japan). The ionosphere was disturbed violently from quiet to severe for these three days in the mid-latitude including the Korean region. We collected GPS carrier-phase measurements with 1-s intervals for the three days from the six reference stations of the National Geographic Information Institute (NGII) in Korea. The GANH station is chosen as a master station, and five other stations are assigned to auxiliary stations to calculate the SD between receivers. All the stations have a Trimble NetR9 receiver connected with a TRM59800 antenna, and their locations are shown in Figure 2. We pre-processed cycle slips and outliers in the collected data by using Bernese GNSS software and calculated the IN acceleration and IP acceleration. We call these values the nominal MVs that represent the non-cycle-slip hypothesis. All baseline data of visible satellites that have elevation angles greater than 5° are used for statistical analysis.

Though we collected the almost 10 6 number of sample data, the assumption of Gaussian distribution cannot be represented by experimental data alone. Therefore, we applied the over-bounding analysis technique that is widely used in aviation area such as SBAS and ground-based augmentation system (GBAS) in order to validate the reliable assumption of the normal error distribution [START_REF] Pullen | Augmented GNSS: Fundamentals and Keys to Integrity and Continuity[END_REF][START_REF] Decleene | Defining Pseudorange Integrity-Overbounding[END_REF][START_REF] Schempp | An Application of Gaussian Overbounding for the WAAS Fault Free Error Analysis[END_REF].

Figure 3 shows the probability density function (PDF) of the two nominal MVs. The blue line represents the actual sample-data distribution normalized by the standard deviation of experimental data, and the red line represent the 1-sigma standard Gaussian distribution. As we can see, the red line appears to be a good statistics of the blue line; however, the both tails of actual error distribution exceed the standard Gaussian distribution. In other words, the probability that the actual nominal MV has a large deviation value (e.g., 4-sigma) under the non-cycle-slip hypothesis exceeds the probability that we expected. On the other hand, the yellow line, which represents the Gaussian distribution based on the theoretically determined standard deviation in Section 2.2.1, properly bounded the both tails of the actual error distribution. Because we calculated the theoretical values conservatively (i.e., worst-case), the probability of potential failure condition or called integrity risk make can be made sufficiently small.

The PDF bounding analysis shown Figure 3 can intuitively check the overall distribution of the sample data; however, it is difficult to identify whether the tail of the distribution is actually bounded under the very small probability level such as 10 -7 since the tail probability is quite sensitive. In addition, it is not possible for bounding all deviation values since its integral over the entire space is equal to one [START_REF] Decleene | Defining Pseudorange Integrity-Overbounding[END_REF]. The cumulative distribution function (CDF) bounding analysis, which overcomes the PDF analysis problem, examines the tails of actual sample distributions to show that they are well-bounded and thinner than the tails of a normal Gaussian distribution [START_REF] Pullen | Augmented GNSS: Fundamentals and Keys to Integrity and Continuity[END_REF][START_REF] Decleene | Defining Pseudorange Integrity-Overbounding[END_REF][START_REF] Schempp | An Application of Gaussian Overbounding for the WAAS Fault Free Error Analysis[END_REF].

Figure 4 represents the folded CDF bounding plots of the two nominal MVs normalized by the standard deviation of experimental data. The plot of a standard CDF has an S-like shape, while the folded CDF plots, which folds the top half of the standard CDF graph over, has mountain shape [START_REF] Monti | Folded Empirical Distribution Function Curves-Mountain Plots[END_REF]. In order to show the tail shape of the distribution more clearly, scale of its vertical axis is logarithmic. The blue curve line represents the actual sample-data CDF, the red line represents the 1-sigma standard CDF, and the yellow line means the conservatively determined CDF in Section 2.2.1. We can clearly see that the blue curve line has non-Gaussian tail, which is fatter than the red line, but thinner than the yellow line. This means that the theoretically determined distribution well bounds the actual sample distributions. For example, as shown in Figure 4a, the probability that the normalized error of actual MVs of IN combination (the blue line) exceeds a 15-sigma value is less than or equal to 10 -7 , under the non-cycle-slip hypothesis; however, the theoretical standard deviation that we conservatively calculated in Section 2.2.1 is designed that the probability of the actual error of MV exceeding a 15-sigma is less than or equal to about 3.0 × 10 -6 . In conclusion, we demonstrate that the theoretically determined standard deviation of the MVs are well-bounded at the 10 -7 -level probability of the normal Gaussian distribution even under severe storms. Table 1 summarizes the statistics of the distribution of the two MVs. The theoretical maximum standard deviation is approximately more three times inflated than the standard deviation of actual data.

Threshold Determination Scheme

Probability of False Alarm and Probability of Missed Detection

In this section, we describe the probability of false alarm and of missed detection, which are the important performance indices of the cycle-slip detection. We already demonstrated that our two independent MVs, the IN acceleration and IP acceleration, could be assumed to have an over-bounded Gaussian distribution with a zero mean and variances of   2 MV and   2 MV , respectively in the cycle-slip-free case (i.e., non-cycle-slip hypotheses H0).

    MV MV MV MV x H MV N f x 2 2 2 0 2 2 1 : ~0, 0, exp 2 2                    , ( 19
)     MV MV MV MV x H MV N f x 2 2 2 0 2 2 1 : ~0, 0, exp 2 2                    . ( 20
)
When a cycle slip occurs, the PDFs of the MVs are shifted by an amount equal to the size of the combined cycle slip. We define the cycle-slip hypothesis Hk corresponding to the cycle-slip (or fault) event k as follows:

      k k k k MV MV MV MV x H MV N f x 2 2 2 2 2 1 : ~, , exp 2 2                             , ( 21
)       k k k k MV MV MV MV x H MV N f x 2 2 2 2 2 1 : ~, , exp 2 2                             . ( 22
)
Figure 5 shows the PDF of two independent MVs under the above hypotheses. Considering a given threshold, there is a possibility that the MV exceeds the threshold when a cycle slip does not occur. Such an event is termed as a false alarm (FA). On the other hand, there is the possibility that the MV does not exceed the threshold when cycle slips are actually present. That is, the shifted bias under the cycle-slip hypothesis can be smaller than the given threshold. Such an event is termed as missed detection (MD) [START_REF] Altmayer | Enhancing the integrity of integrated GPS/INS systems by cycle slip detection and correction[END_REF][START_REF] Song | Predictions of Allowable Sensor Error Limit for Cycle-Slip Detection[END_REF][START_REF] Kim | GPS Cycle Slip Detection Considering Satellite Geometry based on TDCP/INS Integrated Navigation[END_REF]. FA and MD can be derived as follows:

Probability of false alarm for each MV:

  FA MV T P P MV T H 0 2 1                            , ( 23
)   FA MV T P P MV T H 0 2 1                            . ( 24
)
Probability of missed detection for each MV:

  k k MD H k MV T P P MV T H ,                      , ( 25
)   k k MD H k MV T P P MV T H ,                      , ( 26
)             2 1 1 ( ) exp 2 2 x where x z dz . ( 27
)
Since cycle slips are detected using two complementary and independent MVs, the total probability of FA is determined as the sum of the FA rate of each MV. On the other hand, the total probability of MD can be calculated by the product of the MD rate of each MV.

Total probability of false alarm:

 

           . ( 29
)
Figure 5. Probability density function of two monitoring values; the region shaded in orange represents the probability of missed detection, and the sky-blue region represents half of the false-alarm probability.

Detection Threshold Determination

A good detection algorithm has small probabilities of both FA and MD; however, these probabilities are correlated with each other by the detection threshold. If we determine the threshold to minimize FA rate, the probability of MD should be increased and vice versa. Obviously, a threshold of high-integrity and safe detection should be determined to minimize the probability of MD, which is more critical and important than the FA rate [START_REF] Altmayer | Enhancing the integrity of integrated GPS/INS systems by cycle slip detection and correction[END_REF][START_REF] Song | Predictions of Allowable Sensor Error Limit for Cycle-Slip Detection[END_REF][START_REF] Kim | GPS Cycle Slip Detection Considering Satellite Geometry based on TDCP/INS Integrated Navigation[END_REF]. Typically, a detection threshold of integrity monitors or a fault-detection algorithm is determined on the basis of an FA probability allocated from the continuity requirement of the system. The MD probability is computed by a given threshold. Next, we verify whether the MD rate falls within the required level. If the probability of MD is greater than the desired integrity requirement, we need to re-design the detection threshold based on the re-allocated FA rate [START_REF] Pullen | Augmented GNSS: Fundamentals and Keys to Integrity and Continuity[END_REF][START_REF] Joerger | Solution Separation and Chi-Squared ARAIM for Fault Detection and Exclusion[END_REF].

According to the above threshold-determination scheme, we allocated the same probability of FA for each MV as 5.0 × 10 -6 . That is, the desired total FA probability is 10 -5 . The corresponding confidence level KFA is 4.565 that determined by the function "norminv" representing the normal inverse CDF in MATLAB. In conclusion, the proposed cycle-slip detection method introduced in Equation Error! Reference source not found. can be re-expressed with the determined threshold values as follows:

    FA MV
                                      . ( 30
)

Probability of Missed Detection for Insensitive Cycle-Slip Pairs

Figure 6 shows the L1/L2 insensitive cycle-slip pairs for each combination. We assumed that there is only a bias effect k  for each cycle-slip pair in the detection test domain. The size of the bias must at least be greater than the threshold in order to be detectable. The red square points represent the insensitive pairs satisfying Equation Error! Reference source not found. that do not jump significantly on the IN combination. Similarly, the blue square points show the insensitive pairs of cycle slip that satisfy Equation Error! Reference source not found. for the IP combination:

  CS CS T 1 1 2 2 1 1          , ( 31
) CS CS T 1 1 2 2 1 1 2               . ( 32
)
In order to identify the advantage of the proposed method, we have confirmed the insensitive pairs of MW combination satisfying condition as follows:

  WL MW CS CS 1 2 4       , ( 33
)
where the wide-lane (WL) wavelength   0.86 m WL and the threshold for detection is set to 4 sigma of the standard deviation   0.4 m MW that is considering the pseudorange multipath level of CORS receivers [START_REF] Song | A Study on Improving Performance of Network RTK through Tropospheric Modeling for Land Vehicle Applications[END_REF][START_REF] Blewitt | An Automatic Editing Algorithm for GPS data[END_REF][START_REF] Cai | Cycle Slip Detection and Repair for Undifferenced GPS Observations Under High Ionospheric Activity[END_REF]. These insensitive pairs are represented by the green circle points of Figure 6.

As shown in Figure 6, the general method, which employed the IN and MW combination (the red square points and green circle points) as MVs, cannot detect special cycle-slip pairs (4, 3) and [START_REF] Pullen | Augmented GNSS: Fundamentals and Keys to Integrity and Continuity[END_REF][START_REF] Song | A Study on Improving Performance of Network RTK through Tropospheric Modeling for Land Vehicle Applications[END_REF]. The MW combination have a same signs of inclination as the IN combination for insensitive pairs. The pseudorange multipath effect should be smoother for detecting small cycle slips, but there is a limitation because of remaining bias of multipath. On the other hand, the proposed algorithm (the red square points and blue square points) can detect all insensitive cycle-slip pairs since the two complementary combinations have an opposite signs of inclination. Their inclination is completely orthogonal. In other words, they do not share the insensitive pairs with each other. In order to evaluate the detection performance of the proposed method, we calculate the probability of MD for cycle-slip pairs of each linear combination. The evaluated probability of MD is shown in the color map in Figure 7. Figure 7a represents the probability of MD for IN combination by Equation Error! Reference source not found., and Figure 7b represents the probability of MD for IP combination by Equation Error! Reference source not found.. The cycle slips are easily detectable when the color area is blue, while it is difficult to detect when the color changes toward red. The insensitive cycle-slip pairs for each individual combination have a large MD rate close to one; however, the total probability of MD by the product of the MD rate of each combination is very small because the insensitive pairs can surely be detected by the complementary and orthogonal MVs. For example, probability of MD for pair (1, 1) of IN combination is 0.174. That is, the IN combination may fail to detect the pair (1, 1) with a probability of 17.4%, whereas the pair (1, 1) has a MD probability of 4.3 × 10 -8 by the IP combination. Therefore, the total probability of MD is 7.5 × 10 -9 that is the product value of 0.174 and 4.3 × 10 -8 .

The maximum probability of missed detection is 7.5 × 10 -9 when the cycle-slip pair is (1, 1) or (-1, -1). The remaining pairs have a much smaller probability. Consequently, our detection method using the acceleration of the IN and IP combinations ensures high-integrity detection with a maximum MD probability of 7.5 × 10 -9 under the desired total FA probability of 10 -5 . Table 2 summarizes the list of principal insensitive pairs for each combination with their MD probabilities. Even if the cycle-slip pair has an opposite sign, it has the same probability symmetrically.

In addition, we evaluate the probability of MD for general method, which employed the IN and MW combination as MVs, with same strategy. The comparison results of general method and proposed method are summarized in Table 3. As shown in Table 3, the general method may fail to detect some special cycle-slip pairs such as (4, 3), [START_REF] Pullen | Augmented GNSS: Fundamentals and Keys to Integrity and Continuity[END_REF][START_REF] Song | A Study on Improving Performance of Network RTK through Tropospheric Modeling for Land Vehicle Applications[END_REF], and (9, 7) with a probability more than 40%, while the proposed method can be detect all cycle-slip pairs with improved MD probability. Though it can be applied only stationary observations, we verified that the IP combination can be effectively detect cycle-slip with the IN combination over the MW combination. Therefore, the proposed method is suitable for treating dual-frequency cycle slips on network-based RTK technique. Cycle-Slip Compensation Method

Cycle-Slip Identification Method

The detected cycle slip should be removed or repaired in order to maintain the integer ambiguity on the carrier phase. In this section, we discuss the method for cycle-slip compensation (or recovery). The first step is cycle-slip identification for estimating the exact size of the cycle slip. We can construct the observation matrix A for estimating the original L1 and L2 cycle slips with the two independent combinations, IN and IP acceleration:

    T T I I zz zz MV MV CS MV z Ax CS MV where N Q Q diag 2 2 1 2 1 2 1 2 2 2 1 1 1 1 1 1 2 2 ~0, ,                                                                 . ( 34
)
Therefore, the float solution of L1 and L2 cycle slips can be obtained by the weighted least-squares estimation as follows:

    T T zz zz CS x A Q A A Q z CS 1 1 1 1 2 ˆ             , ( 35
)   T xx zz Q A Q A 1 1 ˆˆ   . ( 36
)
These float solutions are needed to fix an integer solution of cycle slip, and we employ the LAMBDA integer mapping method, which is widely used for ambiguity resolution [START_REF] Teunissen | Performance of the LAMBDA Method for Fast GPS Ambiguity Resolution[END_REF]. The LAMBDA method gives the best integer solution through the integer least squares; that is, the success probability of the LAMBDA method is better than that of the popular integer rounding method [START_REF] Teunissen | GNSS Integer Ambiguity Validation : Overview of Theory and Methods[END_REF][START_REF] Teunissen | Success Probability of Integer GPS Ambiguity Rounding and Bootstrapping[END_REF]. Therefore, the LAMBDA method is reasonable for identifying the integer size of cycle slip.

The lower bound (or worst) probability of correct L1/L2 cycle-slip identification can be calculated from the success rate of bootstrapping (or conditional integer rounding) as follows [START_REF] Teunissen | GNSS Integer Ambiguity Validation : Overview of Theory and Methods[END_REF][START_REF] Teunissen | Success Probability of Integer GPS Ambiguity Rounding and Bootstrapping[END_REF]:

  i s i CS P P x x P x x 2 1 1 1 ˆ2 1 2 2                                     , ( 37
)
where the symbol  is the CDF of the standard normal distribution given by Equation

Error! Reference source not found., the conditional standard deviations In order to evaluate the performance of the proposed identification method, we calculate the probability of failure, which refers to the probability that the float solution is pulled to a wrong integer cycle slip. Finally, our identification method ensures high-integrity performance with the upper-bound failure rate f s P P

8 1 1.4 10      .

Cycle-Slip Validation Method

After the cycle-slip identification, we conduct a validation test to ensure whether the actual cycle slip occurred. The outlier, which corresponds to impulse anomaly with an error of an arbitrary size at a specific epoch, appear similar jump as the cycle slip in the time-difference domain. Therefore, we should distinguish them through the validation test.

Once the size of the cycle slip is determined by the identification method, the two MVs used in the cycle-slip detection can be repaired. These repaired MVs cannot exceed the detection threshold if the actual cycle slip is compensated correctly. Therefore, the cycle slips can be validated when the following condition is satisfied:

                  . ( 38
)
The identified cycle slip passing the validation test is compensated and employed as the normal carrier-phase observation, while a phase anomaly that does not satisfy the above condition is considered an outlier that must be removed.

A block diagram of the proposed cycle-slip detection and compensation algorithm is shown in Figure 8.

Algorithm Test Results

Algorithm Test Environment

The proposed algorithm is verified using GPS carrier-phase measurements collected under a severe geomagnetic storm that occurred on 17 March 2015. The test environment is summarized in Table 4. We collected the GPS dualfrequency carrier phase with 1-s intervals from the GANH (37.72° N, 126.39° E) and CHJU (33.51° N, 126.53° E) stations of NGII in Korea. Their locations are shown in Figure 2, and the baseline length between them is 467 km. The maximum ionospheric disturbance was observed between 15:00:00 and 21:00:00 local time in the Korean region (UTC + 9). Figure 9 represents the satellite geometry during the test time in the CHJU station and the doubledifferenced ionospheric delay obtained from the long-baseline RTK system. As can be seen in this figure, doubledifferenced ionospheric delay is dramatically disturbed at a low elevation.

Therefore, we conducted the algorithm test using the low-elevation data of G14, G19, and G27, which are specifically affected by the ionospheric disturbance. Practically, cycle slips frequently occur in the low-elevation data owing to the low signal-to-noise ratio and the possibility of signal blockage [START_REF] Song | A Study on Improving Performance of Network RTK through Tropospheric Modeling for Land Vehicle Applications[END_REF][START_REF] Bisnath | A New Approach to an Old Problem: Carrier-Phase Cycle Slips[END_REF]. First, we confirmed that no real cycle slip occurred in these target data by pre-processing using Bernese GNSS software. Next, we forcibly inserted the simulated cycle slips, which are the principal insensitive pairs listed on Table 2, every 100 epochs into the raw observations.

Data Test Results

We verified the performance of the proposed cycle-slip detection and compensation algorithm through the observation data in which simulated cycle slips were added. The algorithm test results are summarized in Table 5, and the time series of MVs of each satellite are illustrated in Figure 10 As can be seen from the results, all of the cycle slips including the pair (1, 1) for which the probability of MD has the maximum value are correctly detected by the two independent and complementary MVs, IN and IP acceleration. For example, the cycle-slip pair [START_REF] Pullen | Augmented GNSS: Fundamentals and Keys to Integrity and Continuity[END_REF][START_REF] Song | A Study on Improving Performance of Network RTK through Tropospheric Modeling for Land Vehicle Applications[END_REF], one of the most challenging pairs to detect in many previous studies, is not detected by the IN combination but is successfully detected using the IP combination. In contrast, the insensitive pairs in the IP combination are apparently detected under the IN combination.

The proposed algorithm also correctly identifies the original size of L1 and L2 cycle slips. Columns 6 and 7 of Table 5 list the float solutions of the estimated cycle slips. These float solutions are mapped to integer values by the LAMBDA technique and reliably compensated from the validation test.

Conclusions

This paper proposed a cycle-slip detection and compensation algorithm using only dual-frequency carrier-phase stationary observations. In order to detect and repair cycle slips with high-integrity performance, two complementary ionospheric combinations called the IN and IP combinations are employed. They can successfully detect all of the cycle slips since two L1/L2 combinations combine cycle slips with opposite signs for uniquely detecting insensitive pairs. Furthermore, we treat the second-order time-differenced observation as the MVs to reduce the impact of ionospheric delay even for long-baseline distance of CORS network. We verified that the actual error distributions of the MVs under high ionospheric storm are sufficiently bounded and properly assumed at the 10 -7 -level probability of the normal Gaussian distribution from a theoretical analysis. The detection threshold is determined by the total desired probability of FA of 10 -5 . Consequently, our detection method ensures high-integrity performance with a maximum probability of MD of 7.5 × 10 -9 . We also introduced a LAMBDA-based cycle-slip identification method, for which the upper-bound failure rate is 1.4 × 10 -8 . Finally, through a cycle-slip validation test, we can correctly compensate for a real cycle slip and distinguish it from outliers. An algorithm verification test was conducted using actual data collected under a severe ionospheric storm. As a result, all artificially inserted cycle slips, including small cycle slips, were successfully detected and compensated for. In summary, we demonstrated that the proposed method is suitable for monitoring dual-frequency cycle slips on network RTK systems, which should generate high-precision corrections with high-integrity information.

In this paper, although we described an algorithm based on the SD between stationary receivers for applying the network-based RTK technique, we expect that the method can be modified and extended for many dynamic applications through INS integration or Doppler aiding.

  , and SD observations is twice (i.e.,

Figure 1 .

 1 Figure 1. The geomagnetic indices (Kp and Dst) during 16-18 March 2015.

Figure 2 .

 2 Figure 2. Locations of the six reference stations of NGII in Korea.

Figure 3 .

 3 Figure 3. Probability density function of nominal monitoring values: (a) ionospheric negative acceleration; (b) ionospheric positive acceleration.

Figure 4 .

 4 Figure 4. Folded cumulative distribution function of nominal monitoring values (note that the scale of its vertical axis is logarithmic): (a) ionospheric negative acceleration; (b) ionospheric positive acceleration.

Figure 6 .

 6 Figure 6. L1/L2 insensitive cycle-slip pairs; the red square points represent the insensitive pairs of the IN combination, the blue square points represent that of the IP combination, and the green circle points represent that of the MW combination.

Figure 7 .

 7 Figure 7. Probability of missed detection for insensitive cycle-slip pairs: (a) only the IN combination; (b) only the IP combination; (c) total probability of missed detection using the two combinations simultaneously. The maximum probability is 7.5 × 10 -9 when the cycle-slip pair is (1, 1) or (-1, -1).

 ( i 1, 2 

 2 for each frequency) are the square roots of the entries of the diagonal matrix obtained from the LDL decomposition of the covariance matrix xx Q ˆˆ of least-squares estimation.

Figure 8 .

 8 Figure 8. Block diagram of the overall cycle-slip detection and compensation algorithm.

Figure 9 .

 9 Figure 9. (a) Sky plot in the CHJU station; (b) double-differenced ionospheric delay (RTK results).

-Figure 12 .

 12 In these figures, the top subplot represents IN acceleration, and the bottom subplot represents IP acceleration. The blue line shows the MV and the dash-single-dotted magenta line indicates the threshold for cycle-slip detection. The red star points show the missed cycle slip, which correspond to the insensitive pairs for each combination.

Figure 10 .

 10 Figure 10. Cycle-slip detection results for G14.

Figure 11 .

 11 Figure 11. Cycle-slip detection results for G19.

Figure 12 .

 12 Figure 12. Cycle-slip detection results for G27.

  T B   indicate the estimated SD receiver clock drift; m is the number of the TDSD IF observation

	T    ; the superscript j represents satellite index. We must check and screen the float size of cycle slip IF	IF CS , which
	may remain in TDSD IF observation before the clock estimation. In order to remove possible cycle slips without the
	influence of the receiver clock, triple-differenced observations	i	    j T	IF are used;  j i	represents the SD between
	the satellites i and j . Only the TDSD IF observations satisfying the following condition are used to calculate receiver
	clock drift. A threshold of 3 times the standard deviation of triple-differenced IF observations is set. The value
	8		3 2 is due to triple differencing [23]:
			i	j	T	IF	3 8 (	a 1	1	2 ) (	a 2	2	)	2	where j i

Table 1 .

 1 Statistics of actual error distribution of the monitoring values.

	Combination	Number of Sample	Actual Sample Sigma	Theoretical Sigma	Bounding Scale Factor
	Negative (IN)	991,7379	4.6 mm/s 2	15.1 mm/s 2	3.3
	Positive (IP)	991,7379	4.9 mm/s 2	17.1 mm/s 2	3.5

Table 2 .

 2 Principal insensitive pairs of cycle slip with the probability of missed detection of proposed method (any probability less than 10 -100 is entered as zero in the table).

	Cycle Slip	Iono. Negative	Iono. Positive	Total
	(cycle)	Bias (m)	PMD	Bias (m)	PMD	PMD
	(1, 0)	0.294	3.1 × 10 -50	0.095	0.156	4.9 × 10 -51
	(0, 1)	0.378	1.9 × 10 -92	0.074	0.588	1.1 × 10 -92
	(1, 1)	0.083	0.174	0.169	4.3 × 10 -8	7.5 × 10 -9
	(-1, 1)	0.672	0	0.021	1.000	0
	(-1, 2)	1.049	0	0.053	0.928	0
	(-2, 2)	1.343	0	0.042	0.982	0
	(-2, 3)	1.721	0	0.032	0.996	0

Table 3 .

 3 Comparison results of the probability of missed detection (any probability less than 10 -100 is entered as zero in the table).

	Cycle Slip	General Method	Proposed Method
	(cycle)	(IN + MW)	(IN + IP)
	(1, 0)	6.2 × 10 -61	4.9 × 10 -51
	(0, 1)	0	1.1 × 10 -92
	(1, 1)	6.2 × 10 -3	7.5 × 10 -9
	(4, 3)	0.497	0
	(5, 4)	0.613	0
	(8, 6)	9.5 × 10 -4	0
	(9, 7)	0.401	0
	(10, 8)	5.9 × 10 -3	0

Table 4 .

 4 Algorithm test environment.

	Date	17 March 2015
	Time	UTC 06:00:00~11:59:59 (6 h)
	Interval	1 s
	Baseline	GANH-CHJU (467 km)
	Receiver	Trimble NetR9
	Antenna	Trimble Zephyr (TRM59800)
	Kp index	8-(Daily maximum)

Table 5 .

 5 Cycle-slip detection and compensation results.

				Inserted	Monitoring Value	Estimated CS	Validated
	PRN	Epoch	El (°)	Cycle Slip	(meter)	(cycle)	Cycle Slip
				(cycle)	MV -	MV +	L1 CS	L2 CS	(cycle)
		10200	14.81	(10, 8)	-0.083	1.553	10.029	8.029	(10, 8)
		10300	14.09	(5, 4)	-0.043	0.768	4.996	4.006	(5, 4)
		10400	13.37	(-3, 4)	-2.403	0.024	-2.943	4.062	(-3, 4)
	G14	10500	12.65	(-2, 2)	-1.341	-0.039	-1.998	2.011	(-2, 2)
		10600	11.94	(0, 1)	-0.386	0.083	0.030	1.037	(0, 1)
		10700	11.22	(1, 0)	0.289	0.101	1.017	0.020	(1, 0)
		10800	10.52	(1, 1)	-0.084	0.171	1.013	1.015	(1, 1)
		11600	5.99	(1, 1)	-0.088	0.183	1.057	1.072	(1, 1)
		11700	6.67	(1, 0)	0.291	0.103	1.017	0.012	(1, 0)
		11800	7.36	(0, 1)	-0.385	0.070	-0.051	0.979	(0, 1)
	G19	11900	8.05	(-1, 1)	-0.669	-0.041	-1.089	0.935	(-1, 1)
		12000	8.74	(-2, 3)	-1.722	0.027	-2.046	2.942	(-2, 3)
		12100	9.44	(-4, 5)	-3.074	-0.011	-4.020	5.015	(-4, 5)
		12200	10.14	(8, 6)	0.091	1.193	7.956	5.952	(8, 6)
		7900	7.83	(1, 1)	-0.090	0.166	0.941	0.948	(1, 1)
		8000	8.52	(1, 0)	0.296	0.100	1.003	-0.014	(1, 0)
		8100	9.20	(0, 1)	-0.378	0.078	0.007	1.003	(0, 1)
	G27	8200	9.90	(-1, 2)	-1.047	0.059	-0.980	2.010	(-1, 2)
		8300	10.59	(-3, 3)	-2.020	-0.063	-3.011	2.981	(-3, 3)
		8400	11.29	(4, 3)	0.051	0.610	4.019	3.002	(4, 3)
		8500	11.99	(9, 7)	0.000	1.383	9.015	7.014	(9, 7)

Acknowledgments: This work has been supported by a grant (18TLRP-C113269-03) from Transportation logistics research Program funded by Ministry of Land, Infrastructure and Transport (MOLIT) of Korean government, contracted through the SNU-IAMD. This research was supported (in part) by the Institute of Advanced Aerospace Technology at Seoul National University. The Institute of Engineering Research at Seoul National University provided research facilities for this work.