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.ObservationData;
27 import org.orekit.files.rinex.observation.ObservationDataSet;
28 import org.orekit.gnss.GnssSignal;
29 import org.orekit.gnss.MeasurementType;
30 import org.orekit.gnss.SatelliteSystem;
31 import org.orekit.time.AbsoluteDate;
32
33 /**
34 * Phase minus code cycle slip detectors.
35 * The detector is based the algorithm given in <a
36 * href="https://gssc.esa.int/navipedia/index.php/Examples_of_single_frequency_Cycle-Slip_Detectors">
37 * Examples of single frequency Cycle-Slip Detectors</a> by Zornoza and M. Hernández-Pajares. Within this class
38 * a 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 (algorithm 1 on Navipedia).
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 * @since 10.2
49 */
50 public class PhaseMinusCodeCycleSlipDetector extends AbstractCycleSlipDetector {
51
52 /** Order of the polynomial used for fitting. */
53 private final int order;
54
55 /** Threshold above which cycle-slip occurs. */
56 private final double threshold;
57
58 /** Polynomial single frequency cycle-slip detector Constructor.
59 * @param dt time gap threshold between two consecutive measurement (if time between two consecutive measurement is greater than dt, a cycle slip is declared)
60 * @param threshold threshold above which cycle-slip occurs
61 * @param n number of measurement before starting
62 * @param order polynomial order
63 */
64 public PhaseMinusCodeCycleSlipDetector(final double dt, final double threshold,
65 final int n, final int order) {
66 super(dt, n);
67 this.threshold = threshold;
68 this.order = order;
69 }
70
71 /** {@inheritDoc} */
72 @Override
73 protected void manageData(final ObservationDataSet observation) {
74
75 // Extract observation data
76 final SatelliteSystem system = observation.getSatellite().getSystem();
77 final int prn = observation.getSatellite().getPRN();
78 final AbsoluteDate date = observation.getDate();
79
80 // Initialize list of measurements
81 final List<ObservationData> pseudoRanges = new ArrayList<>();
82 final List<ObservationData> phases = new ArrayList<>();
83
84 // Loop on observation data to fill lists
85 for (final ObservationData od : observation.getObservationData()) {
86 if (!Double.isNaN(od.getValue())) {
87 if (od.getObservationType().getMeasurementType() == MeasurementType.PSEUDO_RANGE) {
88 pseudoRanges.add(od);
89 } else if (od.getObservationType().getMeasurementType() == MeasurementType.CARRIER_PHASE) {
90 phases.add(od);
91 }
92 }
93 }
94
95 // Loop on phase measurements
96 for (final ObservationData phase : phases) {
97 // Loop on range measurement
98 for (final ObservationData pseudoRange : pseudoRanges) {
99 // Change unit of phase measurement
100 final double wavelength = phase.getObservationType().getSignal(system).getWavelength();
101 final ObservationData phaseInMeters = new ObservationData(phase.getObservationType(),
102 wavelength * phase.getValue(),
103 phase.getLossOfLockIndicator(),
104 phase.getSignalStrength());
105
106 // Check if measurement frequencies are the same
107 if (phase.getObservationType().getSignal(system) == pseudoRange.getObservationType().getSignal(system)) {
108 // Phase minus Code combination
109 final PhaseMinusCodeCombination phaseMinusCode = MeasurementCombinationFactory.getPhaseMinusCodeCombination(system);
110 final CombinedObservationData cod = phaseMinusCode.combine(phaseInMeters, pseudoRange);
111 final String nameSat = setName(prn, observation.getSatellite().getSystem());
112
113 // Check for cycle-slip detection
114 final boolean slip = cycleSlipDetection(nameSat, date, cod.getValue(), phase.getObservationType().getSignal(system));
115 if (!slip) {
116 // Update cycle slip data
117 cycleSlipDataSet(nameSat, date, cod.getValue(), phase.getObservationType().getSignal(system));
118 }
119 }
120 }
121 }
122
123 }
124
125 /**
126 * Compute if there is a cycle slip at a specific date.
127 * @param nameSat name of the satellite, on the predefined format (e.g. GPS - 07 for satellite 7 of GPS constellation)
128 * @param currentDate the date at which we check if a cycle-slip occurs
129 * @param phaseMinusCode phase measurement minus code measurement
130 * @param signal signal used
131 * @return true if a cycle slip has been detected.
132 */
133 private boolean cycleSlipDetection(final String nameSat, final AbsoluteDate currentDate,
134 final double phaseMinusCode, final GnssSignal signal) {
135
136 // Access the cycle slip results to know if a cycle-slip already occurred
137 final List<CycleSlipDetectorResults> data = getResults();
138 final List<Map<GnssSignal, DataForDetection>> stuff = getStuffReference();
139
140 // If a cycle-slip already occurred
141 if (data != null) {
142
143 // Loop on cycle-slip results
144 for (CycleSlipDetectorResults resultPmC : data) {
145
146 // Found the right cycle data
147 if (resultPmC.getSatelliteName().compareTo(nameSat) == 0 && resultPmC.getCycleSlipMap().containsKey(signal)) {
148 final Map<GnssSignal, DataForDetection> values = stuff.get(data.indexOf(resultPmC));
149 final DataForDetection v = values.get(signal);
150
151 // Check the time gap condition
152 if (FastMath.abs(currentDate.durationFrom(v.getFiguresReference()[v.getWrite()].getDate())) > getMaxTimeBeetween2Measurement()) {
153 resultPmC.addCycleSlipDate(signal, currentDate);
154 v.resetFigures( new SlipComputationData[getMinMeasurementNumber()], phaseMinusCode, currentDate);
155 resultPmC.setDate(signal, currentDate);
156 return true;
157 }
158
159 // Compute the fitting polynomial if there are enough measurement since last cycle-slip
160 if (v.getCanBeComputed() >= getMinMeasurementNumber()) {
161 final List<WeightedObservedPoint> xy = new ArrayList<>();
162 for (int i = 0; i < getMinMeasurementNumber(); i++) {
163 final SlipComputationData current = v.getFiguresReference()[i];
164 xy.add(new WeightedObservedPoint(1.0, current.getDate().durationFrom(currentDate),
165 current.getValue()));
166 }
167
168 final PolynomialCurveFitter fitting = PolynomialCurveFitter.create(order);
169 // Check if there is a cycle_slip
170 if (FastMath.abs(fitting.fit(xy)[0] - phaseMinusCode) > threshold) {
171 resultPmC.addCycleSlipDate(signal, currentDate);
172 v.resetFigures( new SlipComputationData[getMinMeasurementNumber()], phaseMinusCode, currentDate);
173 resultPmC.setDate(signal, currentDate);
174 return true;
175 }
176
177 } else {
178 break;
179 }
180
181 }
182
183 }
184
185 }
186
187 // No cycle-slip
188 return false;
189 }
190
191 }