1 /* Copyright 2002-2025 CS GROUP
2 * Licensed to CS GROUP (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.estimation.measurements.gnss;
18
19 import java.util.ArrayList;
20 import java.util.List;
21 import java.util.Map;
22
23 import org.hipparchus.fitting.PolynomialCurveFitter;
24 import org.hipparchus.fitting.WeightedObservedPoint;
25 import org.hipparchus.util.FastMath;
26 import org.orekit.files.rinex.observation.ObservationDataSet;
27 import org.orekit.gnss.GnssSignal;
28 import org.orekit.gnss.MeasurementType;
29 import org.orekit.gnss.SatelliteSystem;
30 import org.orekit.time.AbsoluteDate;
31
32
33 /**
34 * Geometry free cycle slip detectors.
35 * The detector is based the algorithm given in <a
36 * href="https://gssc.esa.int/navipedia/index.php/Detector_based_in_carrier_phase_data:_The_geometry-free_combination">
37 * Detector based in carrier phase data: The geometry-free combination</a> by Zornoza and M. Hernández-Pajares. Within this class
38 * a second order polynomial is used to smooth the data. We consider a cycle-slip occurring if the current measurement is too
39 * far from the one predicted with the polynomial.
40 * <p>
41 * For building the detector, one should give a threshold and a gap time limit.
42 * After construction of the detectors, one can have access to a List of CycleData. Each CycleDate represents
43 * a link between the station (define by the RINEX file) and a satellite at a specific frequency. For each cycle data,
44 * one has access to the begin and end of availability, and a sorted set which contains all the date at which
45 * cycle-slip have been detected
46 * </p>
47 * @author David Soulard
48 * @author Bryan Cazabonne
49 * @since 10.2
50 */
51 public class GeometryFreeCycleSlipDetector extends AbstractCycleSlipDetector {
52
53 /** Threshold above which cycle-slip occurs. */
54 private final double threshold;
55
56 /**
57 * Constructor.
58 * @param dt time gap threshold between two consecutive measurement (if time between two consecutive measurement is greater than dt, a cycle slip is declared)
59 * @param threshold threshold above which cycle-slip occurs
60 * @param n number of measurement before starting
61 */
62 public GeometryFreeCycleSlipDetector(final double dt, final double threshold, final int n) {
63 super(dt, n);
64 this.threshold = threshold;
65 }
66
67
68 /** {@inheritDoc} */
69 @Override
70 protected void manageData(final ObservationDataSet observation) {
71
72 // Extract observation data
73 final int prn = observation.getSatellite().getPRN();
74 final AbsoluteDate date = observation.getDate();
75 final SatelliteSystem system = observation.getSatellite().getSystem();
76
77 // Geometry-free combination of measurements
78 final GeometryFreeCombination geometryFree = MeasurementCombinationFactory.getGeometryFreeCombination(system);
79 final CombinedObservationDataSet cods = geometryFree.combine(observation);
80
81 // Initialize list of measurements
82 final List<CombinedObservationData> phasesGF = new ArrayList<>();
83
84 // Loop on observation data to fill lists
85 for (final CombinedObservationData cod : cods.getObservationData()) {
86 if (!Double.isNaN(cod.getValue()) && cod.getMeasurementType() == MeasurementType.CARRIER_PHASE) {
87 phasesGF.add(cod);
88 }
89 }
90
91 // Loop on Geometry-free phase measurements
92 for (CombinedObservationData cod : phasesGF) {
93 final String nameSat = setName(prn, observation.getSatellite().getSystem());
94 // Check for cycle-slip detection
95 final GnssSignal signal = cod.getUsedObservationData().get(0).getObservationType().getSignal(system);
96 final boolean slip = cycleSlipDetection(nameSat, date, cod.getValue(), signal);
97 if (!slip) {
98 // Update cycle slip data
99 cycleSlipDataSet(nameSat, date, cod.getValue(), cod.getUsedObservationData().get(0).getObservationType().getSignal(system));
100 }
101 }
102
103 }
104
105 /**
106 * Compute if there is a cycle slip at an specific date.
107 * @param nameSat name of the satellite, on the pre-defined format (e.g.: GPS - 07 for satellite 7 of GPS constellation)
108 * @param currentDate the date at which we check if a cycle-slip occurs
109 * @param valueGF geometry free measurement
110 * @param signal frequency used
111 * @return true if a cycle slip has been detected.
112 */
113 private boolean cycleSlipDetection(final String nameSat, final AbsoluteDate currentDate,
114 final double valueGF, final GnssSignal signal) {
115
116 // Access the cycle slip results to know if a cycle-slip already occurred
117 final List<CycleSlipDetectorResults> data = getResults();
118 final List<Map<GnssSignal, DataForDetection>> stuff = getStuffReference();
119
120 // If a cycle-slip already occurred
121 if (data != null) {
122
123 // Loop on cycle-slip results
124 for (CycleSlipDetectorResults resultGF : data) {
125
126 // Found the right cycle data
127 if (resultGF.getSatelliteName().compareTo(nameSat) == 0 && resultGF.getCycleSlipMap().containsKey(signal)) {
128 final Map<GnssSignal, DataForDetection> values = stuff.get(data.indexOf(resultGF));
129 final DataForDetection dataForDetection = values.get(signal);
130
131 // Check the time gap condition
132 final double deltaT = FastMath.abs(currentDate.durationFrom(dataForDetection.getFiguresReference()[dataForDetection.getWrite()].getDate()));
133 if (deltaT > getMaxTimeBeetween2Measurement()) {
134 resultGF.addCycleSlipDate(signal, currentDate);
135 dataForDetection.resetFigures(new SlipComputationData[getMinMeasurementNumber()], valueGF, currentDate);
136 resultGF.setDate(signal, currentDate);
137 return true;
138 }
139
140 // Compute the fitting polynomial if there are enough measurement since last cycle-slip
141 if (dataForDetection.getCanBeComputed() >= getMinMeasurementNumber()) {
142 final List<WeightedObservedPoint> xy = new ArrayList<>();
143 for (int i = 0; i < getMinMeasurementNumber(); i++) {
144 final SlipComputationData current = dataForDetection.getFiguresReference()[i];
145 xy.add(new WeightedObservedPoint(1.0, current.getDate().durationFrom(currentDate),
146 current.getValue()));
147 }
148
149 final PolynomialCurveFitter fitting = PolynomialCurveFitter.create(2);
150 // Check if there is a cycle_slip
151 if (FastMath.abs(fitting.fit(xy)[0] - valueGF) > threshold) {
152 resultGF.addCycleSlipDate(signal, currentDate);
153 dataForDetection.resetFigures(new SlipComputationData[getMinMeasurementNumber()], valueGF, currentDate);
154 resultGF.setDate(signal, currentDate);
155 return true;
156 }
157
158 } else {
159 break;
160 }
161
162 }
163
164 }
165
166 }
167
168 // No cycle-slip
169 return false;
170 }
171
172 }