周跳的探测及修复

前言：

本章节代码均在Gitee中开源：

导航工程: 导航工程及其有关的所有项目 - Gitee.comhttps://gitee.com/Ehundred/navigation-engineering/tree/master/%E5%8D%AB%E6%98%9F%E5%AF%BC%E8%88%AA%E5%8E%9F%E7%90%86/%E5%91%A8%E8%B7%B3%E6%8E%A2%E6%B5%8B其中涉及到载波定位的知识，再另一章有详细讲解，感兴趣可以移步到：

我们如何收到卫星信号？（导航电文，载波与测距码）_导航测距码-CSDN博客https://blog.csdn.net/qq_74260823/article/details/139411651~~因为这章是学校作业，所以稍微正经点.~~

什么是周跳

载波测量观测值

载波相位测量，因为接收机只可以观测到不足一周的相位，但是卫星发射信号到接收机上，中间一定会有很多整周期波：

所以，实际上的总相位为：

而这里的整周数，因为是无法直接测量出来的，所以我们就称他为——整周模糊度

为了计算出整周模糊度，我们还需要对卫星连续观测。在第一次观测的基础上，我们会不断受到电磁波，不足一周的相位会不断增加。当增加到一周时，就将这一周单独提出来。于是，不足一周的相位部分永远小于一周，而满一周时，便将多出的一周加到整周计数上，于是新的表达式为：

换一种方式解释：

电子钟

我们现在有个电子钟，但是这个电子钟不太智能，只能看到秒数的部分：

但是他的实际时间，肯定不是30秒，因为我们要考虑到看不见的分钟部分。他的实际时间是：

我们一直看这个钟。现在是30秒，慢慢他会变成31，32……，一直到，他变成了61秒

但是，61秒就超过了一分钟，会发生进位。虽然我们看不到分钟，但是当61秒变成1秒的时候，我们还是知道，分钟的部分加了1。此时的时间为：

我们就这样一直看，当秒数从60跳到1时，我们就手动记录多加了一分钟。就算我们不知道分钟数，但是我们可以知道，刚开始的分钟数xx，加上我们手动记录的分钟，就是当前的分钟数，我们把xx当成了一个未知的常量。

于是，新的表达式可以写成：

这样有什么意义？当我们想计算两个时间的差，分钟模糊度会被消掉，剩下的秒数和手动记录的分钟都是已知的，可以直接计算出来。

周跳

还是电子钟。我们一直观测这个电子钟，但是突然，发生了意外，电子钟卡住了！

这个时候，我们看到的秒数是不变的，一直为21秒，而因为一直没有秒数的变化，我们也不敢擅自给手动记录的分钟数进行操作。于是，手动记录的分钟一直停在了5分钟。
突然，电子钟又恢复了正常。我们又看到了秒数，但是这个秒数，因为中间卡了一段时间，所以突然跳到了另一个值；而手动记录的分钟数，因为我们并不知道卡住的时候，过了多少分钟，所以这个5分钟就不可信了。

所以之前的数据就没用了吗？当然不是。虽然手动记录的分钟数差了一个值，但是我们把这个差值，归到本就不知道的分钟模糊度中，这样会产生两个结果：

从5分钟开始继续计数，往前的数据可以使用，往后的数据也可以使用
虽然分钟模糊度变成了一个新值，但是新值和旧值中间只差了一个常数，这个常数我们可以粗略估计出来，虽然是新值，但是可以和旧值联系起来。

这样做有什么意义？如果我们想计算卡住之前的时间和卡住之后的时间的差，如果废弃之前的数据，那么两段时间就完全联系不到一起。但是如果用新的表达式，分钟模糊度还是可以被消掉，手动记录的分钟数也是可以被使用的，两个时间还是可以按之前的方法正常计算，只不过多了一个卡住的时间这一常量。

所以，这个卡住的时间，就叫做周跳；而周跳产生的原因，就是电子钟卡了。

周跳的产生

周跳产生的原因有很多，但是根本原因，都来源于观测的中断。卫星和接收机中间隔了一栋楼，信号接收异常，观测就中断了一段时间；因为周围在考试有电磁干扰，卫星信号很弱，观测也中断了一段时间。
中断的时候，不满一周部分无法正确计数，而整周计数部分因为看到不满一周部分没有变化，所以整周计数也不敢擅自变化，于是观测值就一直停在了一个值了。而等信号恢复正常，接收机接收到了正常的信号，此时不满一周部分会突然变化为另一个值，而整周计数因为一直没变，所以还是从原来的值开始继续计数。但是和电子钟一样，在观测失锁的时候，并不知道过去了多少周，于是就把这个常数放在整周模糊度中，而整周计数按照原来的值继续计数，还是一个可以使用的值。

这个失锁导致未知的常数，就叫做周跳。周跳产生的原因，就是信号中断了。

为什么要探测周跳？

废话，因为计算结果需要周跳。

计算机是很笨的。你们可能觉得，只需要看相位的变化，还有总相位的值，不就可以了吗？但是：

首先，计算机的观测并不是我们想象中的完全连续。计算机的观测是一种均匀中断的连续，啥意思？就是比如每隔1秒观测一次，一直连续观测下去，最终观测到的是一个又一个的点，而非一段连续的三角函数
就算相位发生了极大的变化，计算机也不知道。因为接收机只负责接收，而数据的分析要交给另外的机器去做，也就是我们现在应该去处理的，周跳的探测
相位等值发生较大变化，并不百分百因为产生了周跳，还可能因为观测产生误差，或者的的确确真实数据是如此。所以，探测周跳，还应考虑排除误差

周跳探测的方法

周跳探测，一般会采用以下几种方法：

屏幕扫描法，就是看一眼过去，看看相位有没有太大的变化
高次差法，多项式拟合法
MW组合观测值法
电离层残差法
三差法

这里，我们介绍中间三种用的最多的方法：

高次差法

说简单点，就是差分数组。
（补充：一命通关差分-CSDN博客）

如果在某一个时间点，发生了周跳，往后的所有整周模糊度一定会加上一个周跳常数ΔN，所以我们可以写出整周模糊度的函数表达：

此时，我们把整周模糊度写成一个数组：N[n]，然后进行差分：

但是这可能是误差导致的，再把dif1差分：

就这样一直差分下去，周跳附近的点的差分值会越来越规律，一直到呈现出一个很规律的图形：

我们不断差分，其实就是为了排除误差，而最终出现这样的图形，就可以断定，发生了周跳。

线性组合观测值法

电离层残差（GF）和MW组合观测值法，实际上都是采用线性组合观测值的方法来探测周跳，所以就放一起讲了。

线性组合观测

首先，什么是线性组合观测？我们知道，一段卫星信号，有两个频率的载波，波段分别为L1和L2。
换句话说，对同一个卫星，一个接收站会受到这个卫星发来的两个不同频率的载波，而这个不同频率的载波，产生的观测值也是不同的。

假设L1载波的观测值为φ1，L2载波的观测值为φ2，φ为不足一周的相位。

那么组合观测值，φ=n*φ1+m*φ2，当n和m取不同的值时，产生的不同结果，就叫不同的线性组合观测，常见的有：

无电离层，消去电离层延迟
宽巷，n=1，m=-1
窄巷，n=1，m=1
电离层残差，n=λ1，m=-λ2

所以对一个线性组合观测值，我们直接采用一个基类：

//单个载波
struct carrier
{
	double f;//频率，从L1和L2波段中选择
	double Lambda;//波长,f*l=c
	double Fai;//不足一周计数
	double N;//整周模糊度
	double count;//整周计数
	double distance;//伪距
	double TEC;//电子含量
};


//双载波观测值
class Dual_frequency_carrier
{
public:
	Dual_frequency_carrier(carrier c1,carrier c2)
		:_carrier1(c1),
		_carrier2(c2)
	{}

	carrier _carrier1;
	carrier _carrier2;
};

//线性组合观测值（基类）
class LineObservation
{
public:
	LineObservation(const Dual_frequency_carrier& carriers)
		:_carriers(carriers)
	{}

	double Observation()
	{
		return _observation;
	}
protected:
	Dual_frequency_carrier _carriers;
	double _observation;//线性组合观测值

	virtual double _Observation() = 0;//线性组合的方法
};

而对所有的线性组合，只需要去继承这个基类，然后重写一下线性组合的方法，就可以了。

电离层残差组合

//GF观测值
class GeometryFree : public LineObservation
{
public:
	GeometryFree(const Dual_frequency_carrier& carriers)
		:LineObservation(carriers)
	{
		_Observation();
	}
protected:
    //重写虚函数，即线性组合的方法
	virtual double _Observation()
	{
		carrier car1 = _carriers._carrier1;
		carrier car2 = _carriers._carrier2;
		double A1 = -40.3 * car1.TEC;
		double Vion1 = A1 / pow(car1.f, 2);

		_observation = car1.Lambda * car1.N - car2.Lambda * car2.N + (1 - pow(car1.f, 2) / pow(car2.f, 2)) * Vion1;//直接用推导出的公式
		return _observation;//线性组合观测值LGF
	}
};

此时，我们再来看表达式：

两个载波的波长是已知且为常量的。
电离层延迟在电离层电子含量TEC变化不大的时候，也是不变的，而在短时间内TEC不会发生太大的变化，所有我们也可以认为他是一个常量。
最后只剩下整周模糊度，因为整周模糊度如果不发生周跳，那么整周模糊度也是一个常量，

所以，只要不发生周跳，那么LGF观测值一定不会发生太大的变化
换句话说，如果LGF发生了太大的变化，那么一定发生了周跳

我们把这个变化的最大上限，定为0.05cm，还是采用差分的方法，如果超过了最大值，就表示发生了周跳：

//GF周跳探测
class GF_detect
{
public:
	GF_detect(const vector<Dual_frequency_carrier>& data)
	{
		//把所有的载波变为GF观测值
		for (const auto& e : data)
		{
			_data.push_back(GeometryFree(e));
		}

		_Detect();
	}

	bool Detect()
	{
		return _detect;
	}

	vector<int> Slip()
	{
		return _slip;
	}

private:
	vector<GeometryFree> _data;
	bool _detect;//是否发生周跳，发生周跳为true
	vector<int> _slip;//发生周跳的具体时间点
	double MAX_GAP = 5E-4;//最大阈值

	bool _Detect()
	{
		//遍历每一个历元差，如果差值超过了阈值，那么一定发生了周跳
		for (int i = 1; i < _data.size(); i++)
		{
			if (_data[i].Observation() - _data[i - 1].Observation() > MAX_GAP)
			{
				_detect = true;//发生了周跳
				_slip.push_back(i);//发生周跳的时间点
			}
		}

		return _detect;
	}
};

MW组合观测

为啥叫MW组合？因为是一个姓M和一个姓W的人提出来的。但是推导方法非常非常简单，如果把我放在那个年代，说不定能叫YB观测值。

因为符号太多，所以只能采取手写的方法推导了，见谅

最后得出的结果为：

在这里，我们以周为单位，来进行代码的示例

class MW :public LineObservation
{
public:
	MW(const Dual_frequency_carrier& carriers)
		:LineObservation(carriers)
	{
		_Observation();
	}
private:
	virtual double _Obersvation()
	{
		carrier car1 = _carriers._carrier1;
		carrier car2 = _carriers._carrier2;

		double f1 = car1.f;
		double f2 = car2.f;
		double L1 = car1.Lambda;
		double L2 = car2.Lambda;
		double P1 = car1.distance;
		double P2 = car2.distance;

		_observation = (f1 - f2) / (f1 + f2) * (P1 / L1 + P2 / L2) - (car1.Fai - car2.Fai);//周数
		return _observation;
	}
};

而还是和电离层残差组合一样，

等式左边都是变化不大的值，而等式右边是整周模糊度的差值。按理说，左边计算出来的值，应该是一个变化不会太大的常量，如果变化很大，那么就表示发生了周跳。

前面的不用管，我们只需要看后面的检测方法——用取平均值和计算方差的方法，来检测到底是周跳还是误差：

class MW_detect
{
	MW_detect(const vector<Dual_frequency_carrier>& data)
	{

		//将载波数据计算为MW组合
		for (const auto& e : data)
		{
			_data.push_back(MW(e));
		}

		_Detect();
	}

	bool Detect()
	{
		return _detect;
	}

	vector<int> Slip()
	{
		return _slip;
	}
private:
	vector<MW> _data;
	bool _detect;//是否发生周跳，发生周跳为true
	vector<int> _slip;//周跳具体时间点
	double MAX_GAP = 0.5;//最大阈值：0.5周

	bool _Detect()
	{
		double sum = _data[0].Observation();
		double average = _data[0].Observation();
		double sigma = 0.5;
		bool record = false;

		for (int i = 1; i < _data.size(); i++)
		{
			//如果可能发生周跳
			if (record)
			{
				//判断是周跳还是误差
				if (abs(_data[i].Observation() - _data[i - 1].Observation()) <= 1)
				{
					_detect = true;
					_slip.push_back(i);
				}
			}

			if (abs(_data[i].Observation() - average) >= sigma * 4)
				record = true;//可能发生周跳，也可能是误差

			sum += _data[i].Observation();
			average = sum / (i + 1);
			sigma = sigma * sigma + (pow(_data[i].Observation() - average, 2) - sigma * sigma) / (i + 1);
		}

		return _detect;
	}
};

最后，给自己叠个甲。因为自己才是导航工程大二的本科生，有些概念理解可能不到位，而又想用最容易理解的方式表达出来，所以可能正确性会稍微有些偏差。但是对初学者来说，应该不会存在太大的错误，如果可以帮到你，真的荣幸之极。还有