📄 preprocess.cpp
字号:
//i2nd = NOPNT;
calc[j].token = BREAKS * calc[j].token + jb;
n++;
m = 1;
}
else if ( i2nd < 0 )
{
if ( fabs( comb[j].Bw - comb[i1st].Bw ) > TolStep )
{
calc[i1st].token = DEL_TURBOEDIT_OUTLIER;
ncut++;
fprintf( fpOutPut, "%ld\t", obs[i1st].epoch );
i1st = j;
calc[j].token = BREAKS * calc[j].token + jb;
}
else
{
//i1st = i1st;
i2nd = j;
calc[j].token = BREAKS * calc[j].token + jb;
n++;
m++;
mean = ( comb[i1st].Bw + comb[i2nd].Bw ) / 2.0;
sigma = fabs( comb[i1st].Bw - comb[i2nd].Bw ) / sqrt( 2.0 );
}
}
else
{
if ( fabs( comb[j].Bw - mean ) >= MultiFactor * sigma )
{
for ( k = j + 1, l = 0; k < N, l < 1; k++ )
{
if ( calc[k].token == b )
{
nextPnt = k;
l++;
}
else if ( calc[k].token > b )
{
break;
}
}
if ( l < 1 )
{
calc[j].token = DEL_TURBOEDIT_OUTLIER;
ncut++;
fprintf( fpOutPut, "%ld\t", obs[j].epoch );
if ( m < TEMinArc )
{
for ( k = j - 1, l = 0; k >= 0, l < m; k-- )
{
if ( calc[k].token == BREAKS * b + jb )
{
calc[k].token = DEL_TURBOEDIT_TOOSHORT;
ncut++;
fprintf( fpOutPut, "%ld\t", obs[k].epoch );
l++;
}
//else if ( ... ) { ...... }
}
}
}
else
{
if ( fabs( comb[nextPnt].Bw - comb[j].Bw ) > TolStep )
{
calc[j].token = DEL_TURBOEDIT_OUTLIER;
ncut++;
fprintf( fpOutPut, "%ld\t", obs[j].epoch );
}
else
{
if ( m < TEMinArc )
{
for ( k = j - 1, l = 0; k >= 0, l < m; k-- )
{
if ( calc[k].token == (BREAKS * b + jb ))
{
calc[k].token = DEL_TURBOEDIT_TOOSHORT;
ncut++;
fprintf( fpOutPut, "%ld\t", obs[k].epoch );
l++;
}
//else if ( ... ) { ...... }
}
//jb = jb;
}
else
{
jb++;
}
i1st = j;
i2nd = NOPNT;
calc[j].token = BREAKS * calc[j].token + jb;
n++;
m = 1;
}
}
}
else
{
calc[j].token = BREAKS * calc[j].token + jb;
n++;
m++;
sigma = sqrt( pow( sigma, 2.0 ) * ( ( double )n - 2 ) / ( ( double )n - 1 ) + pow( ( comb[j].Bw - mean ), 2.0 ) / ( double )n );
mean = mean + ( comb[j].Bw - mean ) / ( double )n;
}
}
}
else if ( calc[j].token > b )
{
break;
}
}
i = j - 1;
}
}
fprintf( fpOutPut, "\n共剔除野值（太短弧段历元） %d 个\nTurboEdit方法分段信息\n 历元编码\tb1\tb2\tbw\t记号\n", ncut );
for ( i = 0; i < N; i++ )
{
if ( calc[i].token >= 0 )
{
fprintf( fpOutPut, "%5ld\t%16.5lf\t%16.5lf\t%16.5lf\t%d\n", obs[i].epoch, calc[i].cb1, calc[i].cb2, comb[i].Bw, calc[i].token );
}
else
{
fprintf( fpOutPut, "%5ld\t%s\t%s\t%s\t%d\n", obs[i].epoch, strNOVALUE, strNOVALUE, strNOVALUE, calc[i].token );
}
}
fflush( fpOutPut );
return TURBOEDITDETECTOK;
}
static int TurboEdit_PolyFit( const OBSDAT obs[], const CALCDAT calc[],
COMBI comb[], const long int N )
{
fprintf( fpOutPut, "使用最小二乘法对电离层组合伪距进行多项式拟和\n" );
double tempPi, tempLP;
unsigned int m = ( unsigned int )( ( double )N / 100.0 + 1.0 );//truncation of any fractional part
m = m < 6 ? m : 6;
double * a;
double dt[3];
double * pEpoch = NULL;
double * pPi = NULL;
long int i, j;
int value;
if ( ( a = ( double * )malloc( ( m + 1 ) * sizeof( double ) ) ) == NULL )
return TURBOEDIT_ERR_MEMORY;
if ( ( pEpoch = ( double * )malloc( N * sizeof( double ) ) ) == NULL )
{
free( a );
return TURBOEDIT_ERR_MEMORY;
}
if ( ( pPi = ( double * )malloc( N * sizeof( double ) ) ) == NULL )
{
free( a );
free( pEpoch );
return TURBOEDIT_ERR_MEMORY;
}
for( i = 0, j = 0; i < N; i++ )
{
if ( calc[i].token != DEL_PRECALC_P )
{
pEpoch[j] = obs[i].epoch;
pPi[j] = obs[i].P2 - obs[i].P1;
j++;
}
}
if ( ( value = PolyFit( pEpoch, pPi, j, m, a, dt, LEASTSQUARE ) ) != ( int )m )
{
value = TURBOEDIT_ERR_POLYFIT;
}
else
{
fprintf( fpOutPut, " 拟和数据总量： %ld\n", j );
fprintf( fpOutPut, " 拟和多项式次数： %d\n", value );
fprintf( fpOutPut, " 表达式： %.10lf", a[0] );
for ( i = 1; i <= value; i++ )
fprintf( fpOutPut, "+(%.10lf)*x^%d", a[i], i );
fprintf( fpOutPut, "\n 残差向量2－范数（残差平方和开方）： %.12lf\n", dt[0] );
fprintf( fpOutPut, " 残差向量1－范数： %.12lf\n", dt[1] );
fprintf( fpOutPut, " 残差向量∞－范数： %.12lf\n", dt[2] );
fprintf( fpOutPut, "Pi，deltaLi统计分析\n历元编码\tPi原始值\tPi拟和值\tPi残差\tdeltaLi原始值\t拟和Pi后deltaLi值\tdeltaLi残差\t记号\n" );
for ( i = 0; i < N; i++ )
{
if ( calc[i].token != DEL_PRECALC_PHI && calc[i].token != DEL_PRECALC_P &&
calc[i].token != DEL_HIGHERDIFF_OUTLIER && calc[i].token != DEL_TURBOEDIT_OUTLIER )
{
tempLP = comb[i].LP;
comb[i].LP = comb[i].LP + ( obs[i].P2 - obs[i].P1 ) - PolyVal( ( double )obs[i].epoch, a, value );
}
if ( calc[i].token != DEL_PRECALC_P )
{
tempPi = PolyVal( ( double )obs[i].epoch, a, value );
fprintf( fpOutPut, "%5ld\t%16.10lf\t%16.10lf\t%16.10lf\t", obs[i].epoch, obs[i].P2 - obs[i].P1, tempPi, obs[i].P2 - obs[i].P1 - tempPi );
if ( calc[i].token != DEL_PRECALC_PHI &&
calc[i].token != DEL_HIGHERDIFF_OUTLIER && calc[i].token != DEL_TURBOEDIT_OUTLIER )
{
fprintf( fpOutPut, "%16.10lf\t%16.10lf\t%16.10lf\t", tempLP, comb[i].LP, tempLP - comb[i].LP );
}
else
{
fprintf( fpOutPut, "%s\t%s\t%s\t", strNOVALUE, strNOVALUE, strNOVALUE );
}
fprintf( fpOutPut, "%d\n", calc[i].token );
}
else
{
fprintf( fpOutPut, "%5ld\t%s\t%s\t%s\t%s\t%s\t%s\t%d\n", obs[i].epoch, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, calc[i].token );
}
}
value = TURBOEDITPOLYFITOK;
}
fflush( fpOutPut );
free( a );
free( pEpoch );
free( pPi );
return value;
}
static int TurboEdit_Repair( CALCDAT calc[], COMBI comb[],
const OBSDAT obs[],
const long int N )
{
long int i, j, k;
int b;
int jb;
long int m;
long int Pnt[2];
double meanBw0, meanLP0;
double meanBw, meanLP;
double corrBw, corrLP, corrB1, corrB2;
for ( i = 0; i < N; i++ )
{
if ( calc[i].token >= 0 )
{
b = calc[i].token / BREAKS;
jb = calc[i].token % BREAKS;
fprintf( fpOutPut, "预解算分段编码： %d （TurboEdit方法修复周跳以首段为基准段）\n", b );
fprintf( fpOutPut, " TurboEditda方法分段编码\t各段Bw平均值\t各段LP（deltaLi）平均值\t各段Bw改正值\t各段LP（deltaLi）改正值\t各段B1改正值\t各段B2改正值\n" );
meanBw0 = comb[i].Bw;
meanLP0 = comb[i].LP;
m = 1;
Pnt[0] = i;
for ( j = i + 1; j < N; j++ )
{
if ( calc[j].token == calc[i].token )
{
m++;
meanBw0 = meanBw0 + ( comb[j].Bw - meanBw0 ) / m;
meanLP0 = meanLP0 + ( comb[j].LP - meanLP0 ) / m;
}
else if ( calc[j].token > calc[i].token )
{
break;
}
}
Pnt[1] = j - 1;
fprintf( fpOutPut, " %5ld [ %5ld , %5ld ]\t%16.5lf\t%16.5lf\n", jb, obs[Pnt[0]].epoch, obs[Pnt[1]].epoch, meanBw0, meanLP0 );
while ( j < N && calc[j].token / BREAKS == b )
{
meanBw = comb[j].Bw;
meanLP = comb[j].LP;
m = 1;
Pnt[0] = j;
for ( k = j + 1; k < N; k++ )
{
if ( calc[k].token == calc[j].token )
{
m++;
meanBw = meanBw + ( comb[k].Bw - meanBw ) / m;
meanLP = meanLP + ( comb[k].LP - meanLP ) / m;
}
else if ( calc[k].token > calc[j].token )
{
break;
}
}
Pnt[1] = k - 1;
corrBw = meanBw - meanBw0;
corrLP = meanLP - meanLP0;
corrB2 = ( lambda1 * corrBw - corrLP ) / lambdai;
corrB1 = corrBw + corrB2;
j = k;
for ( k = Pnt[0]; k <= Pnt[1]; k++ )
{
comb[k].Bw -= corrBw;
comb[k].LP -= corrLP;
calc[k].cb1 -= corrB1;
calc[k].cb2 -= corrB2;
}
fprintf( fpOutPut, " %5ld [ %5ld , %5ld ]\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\n",
jb, obs[Pnt[0]].epoch, obs[Pnt[1]].epoch, meanBw, meanLP, corrBw, corrLP, corrB1, corrB2 );
}
i = j - 1;
}
}
fprintf( fpOutPut, "\n观测数据预处理结果\n 历元编码\t改正后Bw\t改正后LP\t改正后B1\t改正后B2\n" );
for ( i = 0; i < N; i++ )
{
if ( calc[i].token >= 0 )
{
fprintf( fpOutPut, "%5ld\t%16.5lf\t%16.5lf\t%16.5lf\t%16.5lf\t%d\n", obs[i].epoch, comb[i].Bw, comb[i].LP, calc[i].cb1, calc[i].cb2, calc[i].token );
}
else
{
fprintf( fpOutPut, "%5ld\t%s\t%s\t%s\t%s\t%d\n", obs[i].epoch, strNOVALUE, strNOVALUE, strNOVALUE, strNOVALUE, calc[i].token );
}
}
fflush( fpOutPut );
return TURBOEDITREPAIROK;
}
static void Patch( CALCDAT calc[], const long int N )
{
long int i;
for ( i = 0; i < N; i++ )
{
if ( calc[i].token < 0 )
{
calc[i].cb1 = 0.0;
calc[i].cb2 = 0.0;
}
}
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -