preprocess.cpp

综合利用三次法和blewitt法检测与修复周跳

//-------------------------------------------------------------------------------------------------
//文件名称：	Preprocess.Cpp
//文件版本：	V 1.2
//摘要：		    使用高次差（Higher Difference）方法和 TurboEdit方法实现对双频载波相位和伪距观测
//				值的预处理，其中使用结构体对原始观测数据和解算数据进行了封装。提供了三种编译模式：
//				1、标准运行模式；2、测试但不实时分析模式；3、测试且实时分析模式
//-------------------------------------------------------------------------------------------------



#include "stdafx.h"										//MSVS8编译器自动生成的代码，在其它编译器下不是必需的


#ifndef _PREPROCESS_H_
#include "Preprocess.h"
#endif


#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include "Calc.h"										//自定义的头文件，提供某些数值算法支持


unsigned int BREAKS = 0;								//全局变量，TurboEdit 方法分段所需乘常数
															//将由 Prepare 函数根据( *pbreaks )重新赋值
#define NOPNT			-1								//宏定义，负值，TurboEdit 方法中标识段内第一、二点状态


	//结构体 COMBI 用于存储观测值的两种线性组合
typedef struct Combi
{
	double Bw;											//bw	双频模糊度互差
	double LP;											//ΔLi	电离层组合后相位（转化为距离）和伪距之差
} COMBI;

	//结构体 TREDIFF 用于存储观测值的两种线性组合的三差值
typedef struct TreDiff
{
	double Bw;											//bw	双频模糊度互差的三差值
	double LP;											//ΔLi	电离层组合后相位（转化为距离）和伪距之差的三差值
} TREDIFF;


/*
 *	以下宏定义标识各个具体函数正常执行，由各个具体函数返回
 */
#define PREPAREOK				0X0A00
#define PRECALCOK				0X0A00
#define HIGHERDIFFOK			0X0A00
#define HIGHERDIFFCALCOK		0X0A00
#define HIGHERDIFFDELOK			0X0A00
#define HIGHERDIFFCORRCSOK		HIGHERDIFFOK
#define TURBOEDITOK				PREPROCESSOK
#define TURBOEDITDETECTOK		0X0A00
#define TURBOEDITPOLYFITOK		0X0A00
#define TURBOEDITREPAIROK		TURBOEDITOK


/*
 *	以下两个条件编译命令确保仅在测试状态下才准备输出测试报告文件和/或启动 MATLAB 引擎
 */
//#if !defined(PREPROCESS_DEBUG) && defined(PREPROCESS_DEBUG_strREPORT)
//#undef PREPROCESS_DEBUG_strREPORT
//#endif

//#ifdef PREPROCESS_DEBUG

#include <stdio.h>
#include <time.h>

FILE * fpOutPut = NULL;									//指向测试报告文件的指针
#define strNOVALUE		"      NULL      "

//#endif

/*
 *	函数介绍：	Preprocess 函数调用的第一个子函数，进行观测数据预处理开始前的准备工作，主要包括：
 *				1、检查观测数据总数是否足够
 *				2、检查客户端函数输入的控制参数是否合理
 *				3、在测试模式下，输出控制参数列表
 *
 *	输入参数：	N		原始观测值数目
 *				tol		预处理的控制参数集
 *				pbreaks		指向分段所需乘常数变量的指针
 *
 *	返回值：	int 类型，标识预处理准备工作的执行状态，非 PREPAREOK 的返回值都是执行错误的
 */
static int Prepare( const long int N,	unsigned int * pbreaks );


/*
 *	函数介绍：	Preprocess 函数调用的第二个子函数，进行基于原始观测数据的计算，即首先计算非修正过的
 *				bw、ΔLi、b1、b2 的值
 *
 *	输入参数：	obs			双频载波相位和伪距的原始观测数据 数组的首地址
 *				calc		解算成果（暂时用于存储未修正过的双频模糊度） 数组的首地址
 *				comb		观测值的两种线性组合（bw、ΔLi） 数组的首地址
 *				N			原始观测值数目
 *				Minphi、Maxphi、MinP、MaxP、MaxLossDat、HDMinArc
 *							从预处理控制参数集中提取出来的若干控制参数
 *
 *	返回值：	int 类型，标识预处理起步计算的执行状态，非 PRECALCOK 的返回值都是执行错误的
 */
static int PreCalc( const OBSDAT obs[],		CALCDAT calc[],
					COMBI comb[],			const long int N,
				    const double Minphi,	const double Maxphi,
					const double MinP,		const double MaxP,
					const int MaxLossDat,	const int HDMinArc );


/*
 *	函数介绍：	Preprocess 函数调用的第三个子函数，使用高次差方法剔除野值和探测、修复周跳
 *
 *	输入参数：	同上，TolZero、TolScale 也是从预处理控制参数集中提取出来的若干控制参数
 *
 *	返回值：	int 类型，标识预处理高次差方法的执行状态，非 HIGHERDIFFOK 的返回值都是执行错误的
 */
static int HigherDiff( const OBSDAT obs[],		CALCDAT calc[],
					   COMBI comb[],			const long int N,
					   const double TolZero,	const double TolScale );


/*
 *	函数介绍：	Preprocess 函数调用的第四个子函数，使用 TurboEdit 方法剔除野值和探测、修复周跳
 *
 *	输入参数：	同上，TEMinArc、MultiFactor、TolStep 也是从预处理控制参数集中提取出来的若干控制参数
 *
 *	返回值：	int 类型，标识预处理 TurboEdit 方法的执行状态，非 TUEBOEDITOK 的返回值都是执行错误的
 */
static int TurboEdit( const OBSDAT obs[],		CALCDAT calc[],
					  COMBI comb[],				const long int N,
					  const int TEMinArc,
					  const double MultiFactor,	const double TolStep );


/*
 *	函数介绍：	高次差方法的第一步，计算差值（指三差值，在测试模式下还包括一差值、二差值）
 *				由 HigherDiff 函数调用
 *
 *	输入参数：	同上
 *				tre			观测值的两种线性组合（bw、ΔLi）的三差值 数组的首地址
 *				BwDiff1		bw 的一差值 数组的首地址
 *				BwDiff2		bw 的二差值 数组的首地址
 *
 *	返回值：	int 类型，标识高次差方法计算差值的执行状态，非 HIGHERDIFFCALCOK 的返回值都是执行错误的
 */
static int HigherDiff_Calc( const CALCDAT calc[],	const COMBI comb[],
						    TREDIFF tre[],

							const OBSDAT obs[],		double BwDiff1[],	double BwDiff2[],

							const long int N );


/*
 *	函数介绍：	高次差方法的第二步，剔除野值（在要剔除的观测数值上用“记号”( token )标识）
 *				由 HigherDiff 函数调用
 *
 *	输入参数：	同上，TolZero、TolScale 是经 HigherDiff 函数传递的若干控制参数
 *
 *	返回值：	int 类型，标识高次差方法剔除野值的执行状态，非 HIGHERDIFFDELOK 的返回值都是执行错误的
 */
static int HigherDiff_DelOutlier( CALCDAT calc[],		const COMBI comb[],
								  TREDIFF tre[],

								  const OBSDAT obs[],	double BwDiff1[],		double BwDiff2[],

								  const long int N,		const double TolZero,	const double TolScale );


/*
 *	函数介绍：	高次差方法的第三步，探测和修复周跳，由 HigherDiff 函数调用
 *
 *	输入参数：	同上
 *
 *	返回值：	int 类型，标识高次差方法探测和修复周跳的执行状态，非 HIGHERDIFFCORRCSOK 的返回值都是执行错误的
 */
static int HigherDiff_CorrCycleslip( CALCDAT calc[],		COMBI comb[],
									 TREDIFF tre[],
//#ifdef PREPROCESS_DEBUG
									 const OBSDAT obs[],	double BwDiff1[],	double BwDiff2[],
//#endif
									 const long int N,		const double TolZero );


/*
 *	函数介绍：	TurboEdit 方法的第一步，剔除野值，进行分段，由 TurboEdit 函数调用
 *
 *	输入参数：	同上，TEMinArc、MultiFactor、TolStep 是经 TurboEdit 函数传递的若干控制参数
 *
 *	返回值：	int 类型，标识 TurboEdit 方法剔除野值和分段的执行状态，非 TURBOEDITDETECTOK 的返回值都是执行错误的
 */
static int TurboEdit_Detect( CALCDAT calc[],	const COMBI comb[],
							 const long int N,
//#ifdef PREPROCESS_DEBUG
							 const OBSDAT obs[],
//#endif
							 const int TEMinArc,	const double MultiFactor,	const double TolStep );


/*
 *	函数介绍：	TurboEdit 方法的第二步，对所有有效的电离层组合 Pi = P2 - P1 进行最小二乘拟和，并评价拟和结果
 *				由 TurboEdit 函数调用
 *
 *	输入参数：	同上
 *
 *	返回值：	int 类型，标识 TurboEdit 方法拟和电离层组合的执行状态，非 TURBOEDITPOLYFITOK 的返回值都表示拟和失败
 */
static int TurboEdit_PolyFit( const OBSDAT obs[],	const CALCDAT calc[],
							  COMBI comb[],			const long int N );


/*
 *	函数介绍：	TurboEdit 方法的第三步，计算段内 bw、ΔLi 的平均值，计算段间 bw、ΔLi、b1、b2 的修正值
 *				完成双频模糊度 b1、b2 的解算，由 TurboEdit 函数调用
 *
 *	输入参数：	同上
 *
 *	返回值：	int 类型，标识 TurboEdit 方法修正 b1、b2 的执行状态，非 TURBOEDITREPAIROK 的返回值都是执行错误的
 */
static int TurboEdit_Repair( CALCDAT calc[],		COMBI comb[],
//#ifdef PREPROCESS_DEBUG
							 const OBSDAT obs[],
//#endif
							 const long int N );


/*
 *	函数介绍：	修饰解算成果，将 calc 数组中对应已被剔除的原始观测值历元编码的 b1、b2 设为零
 *
 *	输入参数：	解算成果 calc 数组的首地址
 *
 *	返回值：	无
 */
static void Patch( CALCDAT calc[], const long int N );




int Preprocess( const OBSDAT obs[],		CALCDAT calc[],
				const long int N,		const TOL tol,
				unsigned int * pbreaks )
{

//#ifdef PREPROCESS_DEBUG

//#ifdef PREPROCESS_DEBUG_strREPORT

	if ( ( fpOutPut = fopen( PREPROCESS_DEBUG_strREPORT, "w" ) ) == NULL )
		return DEBUG_ERR_FAILREPORT;

//#endif
	fflush( fpOutPut );
//#endif

	int value;
	COMBI * pcombi = NULL;

	if ( ( pcombi = ( COMBI * )malloc( N * sizeof( COMBI ) ) ) == NULL )
		return PREP_ERR_MEMORY;

	if ( ( value = Prepare( N,pbreaks ) ) != PREPAREOK )
		goto END;

	if ( ( value = PreCalc( obs, calc, pcombi, N, tol.Minphi, tol.Maxphi, tol.MinP, tol.MaxP, tol.MaxLossDat, tol.HDMinArc ) ) != PRECALCOK )
		goto END;

	if ( ( value = HigherDiff( obs, calc, pcombi, N, tol.TolZero, tol.TolScale ) ) != HIGHERDIFFOK )
		goto END;

	if ( ( value = TurboEdit( obs, calc, pcombi, N, tol.TEMinArc, tol.MultiFactor, tol.TolStep ) ) != TURBOEDITOK )
		goto END;

	Patch( calc, N );

END:

	free( pcombi );

//#ifdef PREPROCESS_DEBUG



	fflush( fpOutPut );

//#ifdef PREPROCESS_DEBUG_strREPORT

	fclose( fpOutPut );

//#endif

//#endif

	return value;
}




static int Prepare( const long int N,	unsigned int * pbreaks )
{
	if ( N < MINN )
		return PREPARE_ERR_N;

	BREAKS = ( unsigned int )( * pbreaks );
	if  ( BREAKS >= INT_MAX )
		return PREPARE_ERR_BREAKS;
	else
		* pbreaks = BREAKS;
	return PREPAREOK;
}




static int PreCalc( const OBSDAT obs[],		CALCDAT calc[],
					COMBI comb[],			const long int N,
				    const double Minphi,	const double Maxphi,
					const double MinP,		const double MaxP,
					const int MaxLossDat,	const int HDMinArc )
{

	double k = ( f1 - f2 ) / ( f1 + f2 );
	long int prePntEph	= obs[0].epoch-1;
	long int i;
	long int j;
	long int jn;
	int b;

	for ( i = 0, b = 0, jn = 0; i < N; i++ )
	{
			comb[i].Bw = obs[i].phi1 - obs[i].phi2 - k * ( obs[i].P1 / lambda1 + obs[i]. P2 / lambda2 );
			comb[i].LP = lambda1 * obs[i].phi1 - lambda2 * obs[i].phi2 + obs[i].P1  - obs[i].P2;
			
			calc[i].cb2 = ( obs[i].P2 - obs[i].P1 + lambda1 * comb[i].Bw + lambda2 * obs[i].phi2 - lambda1 * obs[i].phi1 ) / lambdai;
			calc[i].cb1 = calc[i].cb2 + comb[i].Bw;

			if ( ( obs[i].epoch - prePntEph - 1 ) > MaxLossDat )
			{
				if ( jn < HDMinArc )
				{
					for ( j = i - 1; j >= 0, jn > 0; j-- )
					{
						if ( calc[j].token == b )
						{
							calc[j].token = DEL_PRECALC_TOOSHORT;
							jn--;
						}
					}
					b = b;
				}
				else
				{
					b++;
				}
				calc[i].token = b;
				jn = 1;
			}
	
			else
			{
				calc[i].token = b;
				jn++;
			}
			
			prePntEph = obs[i].epoch;
		
	}

	return PRECALCOK;
}




static int HigherDiff( const OBSDAT obs[],		CALCDAT calc[],
					   COMBI comb[],			const long int N,
					   const double TolZero,	const double TolScale )
{
	int value;
	TREDIFF * pTre;

	if ( ( pTre = ( TREDIFF * )malloc( ( N - 3 ) * sizeof( TREDIFF ) ) ) == NULL )
		return HIGHERDIFF_ERR_MEMORY;

	double * pbDiff1 = NULL;
	double * pbDiff2 = NULL;

	if ( ( pbDiff1 = ( double * )malloc( ( N - 1 ) * sizeof( double ) ) ) == NULL )
	{
		free( pTre );
		return HIGHERDIFF_ERR_MEMORY;
	}
	if ( ( pbDiff2 = ( double * )malloc( ( N - 2 ) * sizeof( double ) ) ) == NULL )
	{
		free( pTre );
		free( pbDiff1 );
		return HIGHERDIFF_ERR_MEMORY;
	}
	if ( ( value = HigherDiff_Calc( calc, comb, pTre,obs, pbDiff1, pbDiff2,N ) ) != HIGHERDIFFCALCOK )
		goto FINISH;
	
	if ( ( value = HigherDiff_DelOutlier( calc, comb, pTre,obs, pbDiff1, pbDiff2,N, TolZero, TolScale ) ) != HIGHERDIFFDELOK )
		goto FINISH;
	if ( ( value = HigherDiff_CorrCycleslip( calc, comb, pTre,obs, pbDiff1, pbDiff2,N, TolZero ) ) != HIGHERDIFFCORRCSOK )
		goto FINISH;
FINISH:

	free( pTre );
	free( pbDiff1 );
	free( pbDiff2 );
	return value;
}




static int HigherDiff_Calc( const CALCDAT calc[],	const COMBI comb[],  TREDIFF tre[],
							const OBSDAT obs[],		double BwDiff1[],	double BwDiff2[],
							const long int N )
{

	long int i, j, k;
	long int nextPnt[3];
	int b;

	for ( i = 0; i < N; i++ )
	{
		if ( ( b = calc[i].token ) >= 0 )
		{
			for ( j = i + 1, k = 0; j < N, k < 3; j++ )
			{
				if ( calc[j].token == b )
				{
					nextPnt[k++] = j;
				}
				else if ( calc[j].token > b )
				{
					break;
				}
			}


			for ( j = i, k = 0; j < nextPnt[2]; j++ )
			{
				if ( calc[j].token != b )
				{}
				else
				{				
					if ( k == 1 )
					{	
						BwDiff1[j-1] = comb[j].Bw - comb[i].Bw;
					}

					if ( k == 2 )
					{	
						BwDiff1[j-1] = comb[j].Bw - comb[nextPnt[0]].Bw;
						BwDiff2[j-2] = comb[j].Bw - 2 * comb[nextPnt[0]].Bw + comb[i].Bw;
					}

					k++;

				}
			}

			j = nextPnt[2];
			nextPnt[2] = i;
			for ( k = 1 + 3; j < N; j++ )