1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
|
/***********************************************/
/**
* @file gnssResiduals2AccuracyDefinition.cpp
*
* @brief Compute antenna accuracies from observation residuals.
*
* @author Torsten Mayer-Guerr
* @date 2012-11-19
*/
/***********************************************/
// Latex documentation
#define DOCSTRING docstring
static const char *docstring = R"(
Compute antenna accuracies from observation \configFile{inputfileResiduals}{instrument}.
The \configFile{inputfileStationInfo}{platform} is needed to assign
the residuals to the equipped antenna at observation times.
The \configFile{outputfileAccuracyDefinition}{gnssAntennaDefinition} contains
at first step the same accuracy information for all antennas as the input file.
Only the azimuth~$A$ and elevation~$E$ dependent grid points of the patterns
where enough residuals are available ($>$ \config{minRedundancy})
are replaced by estimated accuracy
\begin{equation}
\sigma(A,E) = \sqrt{\frac{\sum_i e_i^2(A,E)}{\sum_i r_i(A,E)}},
\end{equation}
where $e_i$ are the azimuth and elevation dependent residuals and $r_i$ the
corresponding redundancies (number of observations minus the contribution to
the estimated parameters).
The \configFile{inputfileAccuracyDefinition}{gnssAntennaDefinition} can be modified
to the demands before with \program{GnssAntennaDefinitionCreate}
(e.g. with \config{antenna:resample}).
To verify the results the \configFile{outputfileAntennaMean}{gnssAntennaDefinition}
and the accumulated \configFile{outputfileAntennaRedundancy}{gnssAntennaDefinition}
of the computed pattern grid points can be written.
Example: Analysis of TerraSAR-X residuals of one month shows that low elevation
GPS satellites are not tracked by the onboard receiver. An estimation of accuracies
for these directions is not possible from the residuals and the apriori accuracies
are left untouched. The other directions show very low phase noise hardly elevation
and azimuth dependent for L2W. A nearly zero mean indicates the use of adequate antennca
center variations in the processing.
\fig{!hb}{0.8}{gnssResiduals2AccuracyDefinition}{fig:gnssResiduals2AccuracyDefinition}{L2W accuracies of TerraSAR-X determined from residuals of one month}
See also \program{GnssResiduals2TransmitterAccuracyDefinition}.
)";
/***********************************************/
#include "programs/program.h"
#include "files/fileInstrument.h"
#include "files/filePlatform.h"
#include "inputOutput/system.h"
/***** CLASS ***********************************/
/** @brief Compute antenna definition from observation residuals.
* @ingroup programsGroup */
class GnssResiduals2AccuracyDefinition
{
public:
void run(Config &config, Parallel::CommunicatorPtr comm);
};
GROOPS_REGISTER_PROGRAM(GnssResiduals2AccuracyDefinition, SINGLEPROCESS, "Compute accuracy definition from observation residuals", Gnss)
GROOPS_RENAMED_PROGRAM(GnssResiduals2AntennaDefinition, GnssResiduals2AccuracyDefinition, date2time(2020, 6, 26))
/***********************************************/
void GnssResiduals2AccuracyDefinition::run(Config &config, Parallel::CommunicatorPtr /*comm*/)
{
try
{
FileName fileNameAntennaMean, fileNameAntennaAccuracy, fileNameAntennaRedundancy;
FileName fileNameStationInfo, fileNameAntenna;
std::vector<FileName> fileNameResiduals;
Bool isTransmitter;
Double thresholdOutlier, minRedundancy;
renameDeprecatedConfig(config, "outputfileAntennaDefinition", "outputfileAntennaMean", date2time(2020, 7, 4));
readConfig(config, "outputfileAccuracyDefinition", fileNameAntennaAccuracy, Config::OPTIONAL, "", "elevation and azimuth dependent accuracy");
readConfig(config, "outputfileAntennaMean", fileNameAntennaMean, Config::OPTIONAL, "", "weighted mean of the residuals");
readConfig(config, "outputfileAntennaRedundancy", fileNameAntennaRedundancy, Config::OPTIONAL, "", "redundancy of adjustment");
readConfig(config, "inputfileAccuracyDefinition", fileNameAntenna, Config::MUSTSET, "", "apriori accuracies");
readConfig(config, "inputfileStationInfo", fileNameStationInfo, Config::MUSTSET, "", "to assign residuals to antennas");
readConfig(config, "isTransmitter", isTransmitter, Config::DEFAULT, "0", "stationInfo is of a transmitter");
readConfig(config, "thresholdOutlier", thresholdOutlier, Config::DEFAULT, "10", "ignore residuals with sigma/sigma0 greater than threshold");
readConfig(config, "minRedundancy", minRedundancy, Config::DEFAULT, "3", "min number of residuals. to estimate sigma");
readConfig(config, "inputfileResiduals", fileNameResiduals, Config::MUSTSET, "", "GNSS receiver residuals");
if(isCreateSchema(config)) return;
// ============================
Platform stationInfo;
if(!fileNameStationInfo.empty())
readFilePlatform(fileNameStationInfo, stationInfo);
GnssType typePRN;
if(isTransmitter)
typePRN = GnssType("***"+stationInfo.markerNumber);
std::vector<GnssAntennaDefinitionPtr> antennaList;
if(!fileNameAntenna.empty())
{
readFileGnssAntennaDefinition(fileNameAntenna, antennaList);
stationInfo.fillGnssAntennaDefinition(antennaList);
}
// ============================
for(const FileName &fileName : fileNameResiduals)
{
logStatus<<"read GNSS residuals <"<<fileName<<">"<<Log::endl;
if(!System::exists(fileName))
{
logWarning<<"file not exist -> continue"<<Log::endl;
continue;
}
InstrumentFile fileReceiver(fileName);
for(UInt arcNo=0; arcNo<fileReceiver.arcCount(); arcNo++)
{
GnssReceiverArc arc = fileReceiver.readArc(arcNo);
for(auto &epoch : arc)
{
// find antenna for epoch
auto ant= stationInfo.findEquipment<PlatformGnssAntenna>(epoch.time);
if(!ant)
continue;
GnssAntennaDefinitionPtr antenna = ant->antennaDef;
if(!antenna)
throw(Exception(epoch.time.dateTimeStr()+": antenna not found: "+ant->str()));
UInt idObs = 0;
for(GnssType satType : epoch.satellite)
{
// find type for the satellite system
UInt idType = 0;
while(epoch.obsType.at(idType) != satType)
idType++;
// azimuth and elevation
if((epoch.obsType.at(idType+0) != (GnssType::AZIMUT + GnssType::L1)) ||
(epoch.obsType.at(idType+1) != (GnssType::ELEVATION + GnssType::L1)) ||
(epoch.obsType.at(idType+2) != (GnssType::AZIMUT + GnssType::L2)) ||
(epoch.obsType.at(idType+3) != (GnssType::ELEVATION + GnssType::L2)))
throw(Exception("azimuth and elevation expected"));
Double azimuth = epoch.observation.at(idObs+0);
Double elevation = epoch.observation.at(idObs+1);
if(isTransmitter)
{
azimuth = epoch.observation.at(idObs+2);
elevation = epoch.observation.at(idObs+3);
}
idObs += 4; // skip azimuth and elevation
idType += 4;
// resiudals, redundancy, sigma/sigma0
while((idType<epoch.obsType.size()) && (idObs<epoch.observation.size()) && (epoch.obsType.at(idType) == satType))
{
GnssType type = epoch.obsType.at(idType++);
Double value = epoch.observation.at(idObs++);
Double redundancy=0, sigma=0;
if((idType < epoch.obsType.size()) && (type == epoch.obsType.at(idType))) // next redundancy?
{
type = epoch.obsType.at(idType++);
redundancy = epoch.observation.at(idObs++);
}
if((idType < epoch.obsType.size()) && (type == epoch.obsType.at(idType))) // next sigma?
{
type = epoch.obsType.at(idType++);
sigma = epoch.observation.at(idObs++);
}
while((idType < epoch.obsType.size()) && (type == epoch.obsType.at(idType))) // other additional information?
idObs++, idType++;
if(value)
for(GnssAntennaPattern &pattern : antenna->patterns)
if(type+satType == pattern.type+typePRN)
{
const UInt idxL = static_cast<UInt>(std::round((Double(azimuth)+2*PI)/(2*PI)*pattern.pattern.rows()))%pattern.pattern.rows();
const UInt idxB = static_cast<UInt>(std::round((PI/2-Double(elevation))/Double(pattern.dZenit)));
if(idxB >= pattern.pattern.columns())
break;
if(!pattern.count.size())
{
pattern.sum = Matrix(pattern.pattern.rows(), pattern.pattern.columns());
pattern.ePe = Matrix(pattern.pattern.rows(), pattern.pattern.columns());
pattern.redundancy = Matrix(pattern.pattern.rows(), pattern.pattern.columns());
pattern.count = Matrix(pattern.pattern.rows(), pattern.pattern.columns());
}
// residuals?
if((redundancy > 0) && (sigma > 0) && (sigma <= thresholdOutlier))
{
const Double p = 1./std::pow(sigma, 2); // weight
pattern.ePe(idxL,idxB) += std::pow(value, 2);
pattern.redundancy(idxL,idxB) += redundancy;
pattern.sum(idxL,idxB) += p * value;
pattern.count(idxL,idxB) += p;
}
// observations?
if(std::isnan(redundancy) || (redundancy <= 0))
{
pattern.ePe(idxL,idxB) += std::pow(value, 2);
pattern.redundancy(idxL,idxB) += 1;
pattern.sum(idxL,idxB) += value;
pattern.count(idxL,idxB) += 1;
}
break;
}
} // while()
} // for(satType)
} // for(epoch)
} // for(arcNo)
} // for(idFile)
// ============================
// only one value at zenith
for(auto &antenna : antennaList)
for(auto &pattern : antenna->patterns)
if(pattern.count.size())
{
copy(Vector(pattern.pattern.rows(), sum(pattern.sum .column(0))), pattern.sum .column(0));
copy(Vector(pattern.pattern.rows(), sum(pattern.ePe .column(0))), pattern.ePe .column(0));
copy(Vector(pattern.pattern.rows(), sum(pattern.redundancy.column(0))), pattern.redundancy.column(0));
copy(Vector(pattern.pattern.rows(), sum(pattern.count .column(0))), pattern.count .column(0));
}
// ============================
if(!fileNameAntennaAccuracy.empty())
{
logStatus<<"write accuracy definition <"<<fileNameAntennaAccuracy<<">"<<Log::endl;
for(auto &antenna : antennaList)
for(auto &pattern : antenna->patterns)
if(pattern.count.size())
{
pattern.offset = Vector3d();
for(UInt i=0; i<pattern.pattern.rows(); i++)
for(UInt k=0; k<pattern.pattern.columns(); k++)
if(pattern.redundancy(i, k) >= minRedundancy)
pattern.pattern(i, k) = std::sqrt(pattern.ePe(i, k)/pattern.redundancy(i, k));
}
writeFileGnssAntennaDefinition(fileNameAntennaAccuracy, antennaList);
}
for(auto &antenna : antennaList)
for(auto &pattern : antenna->patterns)
{
pattern.offset = Vector3d();
pattern.pattern *= NAN_EXPR;
}
if(!fileNameAntennaMean.empty())
{
logStatus<<"write antenna definition <"<<fileNameAntennaMean<<">"<<Log::endl;
for(auto &antenna : antennaList)
for(auto &pattern : antenna->patterns)
if(pattern.count.size())
for(UInt i=0; i<pattern.pattern.rows(); i++)
for(UInt k=0; k<pattern.pattern.columns(); k++)
pattern.pattern(i, k) = pattern.sum(i, k)/pattern.count(i, k);
writeFileGnssAntennaDefinition(fileNameAntennaMean, antennaList);
}
if(!fileNameAntennaRedundancy.empty())
{
logStatus<<"write redundancy <"<<fileNameAntennaRedundancy<<">"<<Log::endl;
for(auto &antenna : antennaList)
for(auto &pattern : antenna->patterns)
if(pattern.count.size())
for(UInt i=0; i<pattern.pattern.rows(); i++)
for(UInt k=0; k<pattern.pattern.columns(); k++)
if(pattern.redundancy(i, k) >= minRedundancy)
pattern.pattern(i, k) = pattern.redundancy(i, k);
writeFileGnssAntennaDefinition(fileNameAntennaRedundancy, antennaList);
}
}
catch(std::exception &e)
{
GROOPS_RETHROW(e)
}
}
/***********************************************/
|
