一种基于球谐展开的GNSS精密单点定位方法

一种基于球谐展开的gnss精密单点定位方法
技术领域
1.本发明属于卫星导航定位技术领域，涉及一种基于球谐展开的gnss精密单点定位方法。


背景技术：

2.gnss单点定位方法，是指根据卫星星历给出的观测瞬间卫星在空间的位置和卫星钟差，以及由一台gnss接收机所测定的从卫星到gnss接收机之间的距离，通过距离交会的方法，独立测定该gnss接收机在地球坐标系中的三维坐标的定位方法。
3.利用广播星历提供的卫星位置和卫星钟差，以及伪距观测值进行的定位方式，称为标准单点定位。由于广播星历和伪距观测值精度的限制，且其数学模式中各项改正误差的影响又无法完全消除，因而，标准单点定位的精度通常只能达到米级。因此，标准单点定位主要用于导航定位及资源调查、勘测等低精度定位领域。利用精密星历和精密卫星钟差数据，由相位观测数据进行的定位方式，称为精密单点定位，精密单点定位的精度可以达到厘米量级。
4.在影响单点定位结果的三类误差源中，与信号传播有关的误差(主要是指电离层延迟和对流层延迟误差)无法很好地被模型化，因此，对定位的影响较大。
5.对于电离层延迟，单频接收机需要通过模型加以改正；双频接收机则可以通过线性组合法消除电离层延迟的影响；对于对流层延迟，由于对流层的性质为非弥散性，故不能用双频观测数据线性组合的方法消除对流层延迟误差，其对定位精度的影响不可以忽略。
6.对流层延迟分为干延迟和湿延迟两部分。在通常的定位解算策略中，天顶方向的干延迟由对流层延迟模型计算得出，通过映射函数将其映射到传播路径；而天顶方向的湿延迟则作为未知参数求解，并通过映射函数映射到传播路径。gnss测量中对流层延迟模型或方法有萨斯塔莫宁模型、unb3模型、气象元素法等。
7.以上三种对流层延迟改正方法存在如下问题：萨斯塔莫宁模型的精度相对较高，但是需要用到实测的气象数据，气象数据的获取问题，限制了萨斯塔莫宁模型的应用；unb3模型不需要实测气象数据注入，但是unb3模型的精度有限；气象元素法需要用到外部气象设备采集数据，不同设备的精度存在差异且设备笨重、昂贵，增加了外业采集数据的工作量；同时气象元素法也降低了单点定位的工作效率，很难为全天候的单点定位提供所需数据。
8.以上三种对流层延迟改正方法构建多种模型，以适应不同条件下改正对流层延迟的需要，考虑到实用性和普适性会使得测量外业和内业工作量增大，数据处理效率降低。此外，对于单点定位来说，对流层的延迟误差校正多采用经验模型，此类经验模型的参数估计也多采用经验值，因而，难以符合实际单点定位中的参数估计。


技术实现要素：

9.本发明的目的之一在于提出一种基于球谐展开的gnss标准单点定位方法，该方法
基于球谐展开表示与测站和卫星间的方位角和高度角有关的误差项，以便能够快速、高效、准确地进行标准单点定位。本发明为了实现上述目的，采用如下技术方案：
10.一种基于球谐展开的gnss标准单点定位方法，包括如下步骤：
11.i.1.建立利用伪距观测值进行的标准单点定位观测方程，如公式(1)所示：
[0012][0013]
式中，i表示第i观测历元，j表示卫星编号；ρ表示伪距观测值，表示卫星j在第i观测历元的伪距观测值；表示接收机的位置与卫星j在第i观测历元的位置之间的几何距离；c表示光速；tr表示测站接收机钟差，(tr)i表示第i观测历元测站的接收机钟差；ts表示卫星钟差，表示卫星j在第i观测历元的钟差；表示卫星j在第i观测历元的信号传播路径上的对流层延迟误差；表示卫星j在第i观测历元的信号传播路径上的电离层延迟误差；ε
ρ
表示伪距观测数据残差；
[0014]
定义卫星j在第i观测历元观测瞬间的空间三维坐标为测站的近似坐标为(x0,y0,z0)，则卫星到测站近似位置的几何距离表示为：
[0015][0016]
i.2.基于球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，误差项包括对流层延迟误差以及电离层延迟误差；基于球谐展开的标准单点定位观测方程表示为：
[0017][0018]
式中，n为球谐展开的阶数，m为球谐展开的次数，nmax为球谐展开的最大阶数；表示n阶m次的缔合勒让德多项式；c
nm
和s
nm
分别表示n阶m次的球谐展开的系数，c
nm
和s
nm
为球谐展开的待求参数；分别表示测站与卫星j在第i观测历元之间的高度角和方位角；为了方便表示球谐展开，令：
[0019][0020]
则公式(2)简化为：
[0021]
由简化后的标准单点定位观测方程(3)，得到标准单点定位误差方程(4)，公式如下：
[0022][0023]
其中，v
ρ
表示伪距观测数据的改正值；
[0024]
表示卫星j在第i观测历元的伪距观测数据的改正值；
[0025]
i.3.基于公式(4)得到标准单点定位误差方程的线性化表达式，并进行简化处理，然后对简化后的标准单点定位误差方程的线性化表达式进行滑动解算；
[0026]
i.3.1.对标准单点定位误差方程(4)在测站的近似坐标(x0,y0,z0)处进行泰勒级数展开并保留一阶项，则标准单点定位误差方程的线性化表达式如公式(5)所示；
[0027][0028]
公式(5)中，令：
[0029][0030][0031][0032]
其中，表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dx的系数，待求参数dx是测站近似坐标x0的改正数；表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dy的系数，待求参数dy是测站近似坐标y0的改正数；表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dz的系数，待求参数dz是测站近似坐标z0的改正数；表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数c
nm
的系数；表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数s
nm
的系数；d(tr)i表示在第i观测历元测站接收机钟差的待求参数；
[0033]
则简化后的标准单点定位误差方程的线性化表达式，如公式(6)所示；
[0034][0035]
公式(6)中，i＝1,2,
…
,epoch，epoch表示最大观测历元数；表示由卫星j在第i观测历元伪距观测值、卫星j在第i观测历元的钟差以及测站近似坐标与卫星j在第i观测历元之间的几何距离计算出的常数项，其中，
[0036]
设定滑动窗口选择i个历元，每个历元最多能够观测到j颗卫星，则在这个滑动窗口下所有可视卫星的伪距观测值组成imax个方程，其中，imax＝i
×
j；
[0037]
定义标准单点定位误差方程的系数矩阵、伪距观测值的改正数向量矩阵、常数项矩阵以及未知参数矩阵分别为矩阵b、v、l、x，则它们的表达式分别如下：
[0038][0039][0040][0041]
x＝[dxdydz(dtr)1…
(dtr)
ic00
…cnmaxnmaxsnmaxnmax
]
t
；
[0042]
则标准单点定位误差方程的线性化表达式(6)用一个矩阵形式表达，如公式(7)所示：
[0043]
v＝bx-l(7)
[0044]
未知参数x的最小二乘估值为：x＝(b
t
pb)-1bt
pl(8)
[0045]
公式(8)中，p为单位权矩阵；
[0046]
i.3.2.判断公式(7)中系数矩阵b是否为病态，经过判断，若系数矩阵b是病态，则执行步骤i.3.3；否则，系数矩阵b是良态，执行步骤i.3.4；
[0047]
i.3.3.首先利用截断奇异值分解的方法，为最小二乘谱修正迭代提供待求解的未知参数x的初值，然后基于最小二乘谱修正迭代计算未知参数x的值；
[0048]
设系数矩阵b∈rn×m，rn×m表示n行m列的实数阵；则系数矩阵b的奇异值分解为：
[0049]
b＝usv
t
(9)
[0050]
公式(9)中，u∈n
×
n，v∈m
×
m，u、v均为正交矩阵，s∈n
×
m为一对角矩阵；
[0051]
系数矩阵b∈rn×m的截断奇异值矩阵bk定义为：
[0052][0053]
公式(10)中，矩阵s中最小的(r-k)个非奇异值用零代替，即被截断，且k≤r；r表示系数矩阵b的秩，k表示矩阵s中保留的奇异值数量；
[0054]
ui表示矩阵u对应的向量，vi表示矩阵v对应的向量，σi表示矩阵s中保留的奇异值；
[0055]
计算矩阵s中奇异值的均值，将该均值作为截断奇异值的阈值，其中，矩阵s中奇异值大于截断奇异值的阈值的保留，小于截断奇异值的阈值的进行零化处理；
[0056]
则公式(7)的截断奇异值解为：
[0057]
公式(11)中，
[0058]
根据最小二乘原理v
t
pv＝min，公式(8)写成以下形式：
[0059]
(b
t
pb)x＝(b
t
pl)(12)
[0060]
将公式(12)的等式左端和右端同时加未知参数kx，并化简得到公式(13)：
[0061]
(b
t
pb ki)x＝b
t
pl kx(13)
[0062]
由公式(13)整理得，最小二乘谱修正迭代公式为：
[0063]
x＝(b
t
pb ki)-1
(b
t
pl kx)(14)
[0064]
公式(14)中，k为任意实数；
[0065]
基于最小二乘谱修正迭代计算未知参数x的解，转到步骤i.3.5；
[0066]
i.3.4.利用最小二乘法求解未知参数x的解，转到步骤i.3.5；
[0067]
i.3.5.根据计算得到的未知参数x的解，得到未知参数x中包含的参数值，即测站的位置参数(dx,dy,dz)、接收机钟差dtr以及球谐展开的系数c
nm
和s
nm
；
[0068]
将未知参数x中包含的位置参数(dx,dy,dz)，分别加到测站的近似坐标(x0,y0,z0)上，即标准单点定位求出的地心地固坐标系下的测站坐标(x,y,z)；
[0069]
i.3.6.将测站坐标由地心地固坐标系转换到站心坐标系，实现gnss标准单点定位。
[0070]
本发明的目的之二在于提出一种基于球谐展开的gnss精密单点定位方法，该方法基于球谐展开表示与测站和卫星间的方位角和高度角有关的误差项，以便能够快速、高效、准确地进行精密单点定位。本发明为了实现上述目的，采用如下技术方案：
[0071]
一种基于球谐展开的gnss精密单点定位方法，包括如下步骤：
[0072]
ii.1.建立利用载波相位观测值进行的精密单点定位观测方程，如公式(1)所示：
[0073][0074]
公式(1)中，i表示第i观测历元；j表示卫星编号；表示载波相位观测值，即等效距离，表示卫星j在第i观测历元的载波相位观测值，即等效距离；表示接收机的位置与卫星j在第i观测历元的位置之间的几何距离；tr表示测站接收机钟差，(tr)i表示第i观测历元测站接收机钟差；ts表示卫星钟差，表示卫星j在第i观测历元的钟差；表示卫星j在第i观测历元的信号传播路径上的对流层延迟误差；表示卫星j在第i观测历元的信号传播路径上的电离层延迟误差；n表示卫星相位观测数据的整周模糊度参数；表示卫星j在第i观测历元的相位观测数据整周模糊度参数；c表示光速；λ表示载波的波长；表示载波相位观测数据残差；
[0075]
定义卫星j在第i历元观测瞬间的空间三维坐标为测站的近似坐标为(x0,y0,z0)，则卫星到测站近似位置的几何距离表示为：
[0076][0077]
ii.2.基于球谐展开表示与测站和卫星间的高度角和方位角有关的误差项，误差项包括对流层延迟误差以及电离层延迟误差；基于球谐展开的精密单点定位观测方程表示为：
[0078][0079]
式(2)中，n为球谐展开的阶数，m为球谐展开的次数，nmax为球谐展开的最大阶数；表示n阶m次的缔合勒让德多项式；c
nm
和s
nm
表示n阶m次的球谐展开的系数，c
nm
和s
nm
为球谐展开待求参数；和分别表示测站与卫星j在第i观测历元之间的高度角和方位角；为了方便表示球谐展开，令：
[0080][0081]
则基于球谐展开的精密单点定位观测方程(2)简化为公式(3)；
[0082][0083]
由公式(3)中精密单点定位观测方程得到精密单点定位误差方程，如公式(4)所示；
[0084][0085]
其中，表示载波相位观测数据的改正值；
[0086]
表示卫星j在第i观测历元的载波相位观测数据的改正值；
[0087]
ii.3.基于公式(4)得到精密单点定位误差方程的线性化表达式，并进行简化处理，然后对简化后的精密单点定位误差方程的线性化表达式进行滑动解算；
[0088]
ii.3.1.对精密单点定位误差方程(4)在测站的近似坐标(x0,y0,z0)处进行泰勒级数展开并保留一阶项，则精密单点定位误差方程的线性化表达式如公式(5)所示；
[0089][0090]
公式(5)中，表示卫星j在第i历元的相位观测数据模糊度初值；表示整周模糊度参数的改正数；d(tr)i表示在第i观测历元测站接收机钟差的待求参数；
[0091]
公式(5)中，令：
[0092][0093][0094][0095]
为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dx的系数，待求
参数dx是测站近似坐标x0的改正数；为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dy的系数，待求参数dy是测站近似坐标y0的改正数；为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dz的系数，待求参数dz是测站近似坐标z0的改正数；a
nm
表示球谐展开待求参数c
nm
的系数，b
nm
表示球谐展开待求参数s
nm
的系数；表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数c
nm
的系数，表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数s
nm
的系数；
[0096]
简化后的精密单点定位误差方程的线性化表达式，如公式(6)所示：
[0097][0098]
公式(6)中，i＝1,2,,
…
,epoch，epoch表示最大观测历元；
[0099]
表示由卫星j在第i观测历元载波相位观测值、卫星j在第i观测历元的钟差、卫星j在第i历元的相位观测数据整周模糊度初始值以及测站近似坐标与卫星j在第i观测历元之间的几何距离计算出的常数项，
[0100]
设定滑动窗口选择i个历元，每个历元最多能够观测到j颗卫星，则在这个滑动窗口下所有可视卫星的载波相位观测值组成imax个方程，其中，imax＝i
×
j；
[0101]
定义精密单点定位误差方程的系数矩阵、载波相位观测值的改正数向量矩阵、常数项矩阵以及未知参数矩阵分别为矩阵b、v、l、x，则它们的表达式分别如下：
[0102][0103][0104][0105][0106]
表示第1颗卫星在第1个观测历元的整周模糊度参数，表示第j颗卫星在第i个观测历元的整周模糊度参数；
[0107]
则精密单点定位误差方程的线性化表达式(6)用一个矩阵形式表达，如公式(7)所示；
[0108]
v＝bx-l(7)
[0109]
未知参数x的最小二乘估值为：x＝(b
t
pb)-1bt
pl(8)
[0110]
公式(8)中，p为单位权矩阵；
[0111]
ii.3.2.判断公式(7)的系数矩阵b是否为病态，若系数矩阵b是病态，则执行步骤ii.3.3；否则，若系数矩阵b为良态，执行步骤ii.3.4；
[0112]
ii.3.3.首先利用截断奇异值分解的方法，为最小二乘谱修正迭代提供待求解的未知参数x的初值，然后基于最小二乘谱修正迭代计算未知参数x的值；
[0113]
设系数矩阵b∈rn×m，rn×m表示n行m列的实数阵；则系数矩阵b的奇异值分解为：
[0114]
b＝usv
t
(9)
[0115]
公式(9)中，u∈n
×
n，v∈m
×
m，u、v均为正交矩阵，s∈n
×
m为一对角矩阵；
[0116]
系数矩阵b∈rn×m的截断奇异值矩阵bk定义为：
[0117][0118]
公式(10)中，矩阵s中最小的(r-k)个非奇异值用零代替，即被截断，且k≤r；r表示系数矩阵b的秩，k表示矩阵s中保留的奇异值数量；
[0119]
ui表示矩阵u对应的向量，vi表示矩阵v对应的向量，σi表示矩阵s中保留的奇异值；
[0120]
计算矩阵s中奇异值的均值，将该均值作为截断奇异值的阈值，其中，矩阵s中奇异值大于截断奇异值的阈值的保留，小于截断奇异值的阈值的进行零化处理；
[0121]
则公式(7)的截断奇异值解为：
[0122]
公式(11)中，
[0123]
根据最小二乘原理v
t
pv＝min，公式(8)写成以下形式：
[0124]
(b
t
pb)x＝(b
t
pl)(12)
[0125]
将公式(12)的等式左端和右端同时加未知参数kx，并化简得到公式(13)：
[0126]
(b
t
pb ki)x＝b
t
pl kx(13)
[0127]
由公式(13)整理得，最小二乘谱修正迭代公式为：
[0128]
x＝(b
t
pb ki)-1
(b
t
pl kx)(14)
[0129]
公式(14)中，k为任意实数；
[0130]
基于最小二乘谱修正迭代计算未知参数x的解，转到步骤ii.3.5；
[0131]
ii.3.4.利用最小二乘法求解未知参数x的解，转到步骤ii.3.5；
[0132]
ii.3.5.根据计算得到的未知参数x的解，得到未知参数x中包含的参数值，即测站的位置参数(dx,dy,dz)、接收机钟差dtr、模糊度参数dn以及球谐展开的系数c
nm
和s
nm
；
[0133]
其中，将未知参数x中包含的位置参数(dx,dy,dz)，分别加到测站的近似坐标(x0,y0,z0)上即精密单点定位求出的地心地固坐标系下的测站坐标；
[0134]
ii.3.6.将测站坐标由地心地固坐标系转换到站心坐标系，实现gnss精密单点定位。
[0135]
本发明具有如下优点：
[0136]
如上所述，本发明述及了一种基于球谐展开的gnss单点定位方法，该方法采用球谐展开描述与测站和卫星之间的高度角和方位角有关的误差项，包括对流层延迟误差以及电离层延迟误差，使用gnss观测数据和卫星星历数据即可实现测站空间位置的计算。本发
明能高效、快速地获得测点的位置信息，具有实际操作简单、数据处理便捷以及计算效率高等优势。
附图说明
[0137]
图1为本发明实施例1中基于球谐展开的gnss标准单点定位的流程示意图；
[0138]
图2为本发明实施例2中基于球谐展开的gnss精密单点定位的流程示意图；
[0139]
图3为本发明实施例基于球谐展开的gnss标准单点定位方法的结果统计图；
[0140]
图4为本发明实施例基于球谐展开的gnss精密单点定位方法的结果统计图。
具体实施方式
[0141]
下面结合附图以及具体实施方式对本发明作进一步详细说明：
[0142]
实施例1
[0143]
本实施例述及了一种基于球谐展开的gnss标准单点定位方法。如图1所示，该基于球谐展开的gnss标准单点定位方法包括如下步骤：
[0144]
i.1.建立利用伪距观测值进行的标准单点定位观测方程。
[0145]
实现标准单点定位，读取观测数据和广播星历数据，处理伪距观测数据，计算与测站和卫星间的高度角和方位角无关或关系不密切的误差，包括卫星星历误差和卫星钟差。由于卫星星历误差的量级与标准单点定位误差的量级基本相同，因此，广播星历可用于低精度的标准单点定位，在标准单点定位中卫星星历误差可以暂时忽略。
[0146]
在标准单点定位中，由广播星历内插计算得到每个历元对应的卫星位置；而卫星钟差使用广播星历提供的卫星钟差参数求得，一般卫星钟差ts由一个多项式拟合得到：ts＝α0 α1(t
i-t0) α2(t
i-t0)2，α0、α1和α2分别为钟差、钟漂和钟的频漂系数，可以在导航电文中读取；t0为多项式计算卫星钟差的参考时刻，ti表示信号接收时刻即数据记录时刻。
[0147]
地面接收机钟大多采用石英钟，地面接收机的精度低于卫星钟，钟差变化无规律性，无法通过多项式较好拟合，因此，将接收机钟差tr作为待求参数进行估计。
[0148]
建立利用伪距观测值进行的标准单点定位观测方程，如公式(1)所示：
[0149][0150]
式中，i表示第i观测历元，j表示卫星编号；ρ表示伪距观测值，表示卫星j在第i观测历元的伪距观测值；表示接收机的位置与卫星j在第i观测历元的位置之间的几何距离；c表示光速；tr表示测站接收机钟差，(tr)i表示第i观测历元测站的接收机钟差；ts表示卫星钟差，表示卫星j在第i观测历元的钟差；表示卫星j在第i观测历元的信号传播路径上的对流层延迟误差；表示卫星j在第i观测历元的信号传播路径上的电离层延迟误差；ε
ρ
表示伪距观测数据残差；定义卫星j在第i观测历元观测瞬间的空间三维坐标为测站的近似坐标为(x0,y0,z0)，则卫星到测站近似位置的几何距离表示为：
[0151][0152]
gnss伪距观测数据类型分为单频观测数据和双频观测数据。若进行标准单点定位的gnss接收机为双频接收机，则电离层延迟误差通过双频消电离层组合来消除，组合形式如下：ρ为伪距消电离层组合观测值；f1、f2为相应观测值的频率；ρ1、ρ2分别表示两种频率对应的伪距观测值。若观测数据为双频观测值时，则公式(1)中，电离层延迟误差的一阶主项通过双频消电离层组合消除，电离层延迟误差项为0。此时，本实施例1中的标准单点定位观测方程中，球谐展开主要用于表示对流层延迟误差项，球谐展开的阶次与单频伪距观测值构建的观测方程有所不同。
[0153]
i.2.在标准单点定位中，对流层延迟误差和电离层延迟误差都是卫星信号穿过大气层时产生的一种信号延迟误差，都与卫星信号的传播有关；而卫星信号的传播路径又与测站和卫星的位置有关，利用测站坐标和卫星位置可以计算出测站和卫星之间的高度角和方位角，应用球谐展开，就能够表示与测站和卫星间的高度角和方位角有关的误差项，包括对流层延迟误差和电离层延迟误差，以消除或最大程度的削弱它们对标准单点定位结果的影响。基于球谐展开的标准单点定位观测方程表示为：
[0154][0155]
式中，n为球谐展开的阶数，m为球谐展开的次数，nmax为球谐展开的最大阶数；表示n阶m次的缔合勒让德多项式；c
nm
和s
nm
分别表示n阶m次的球谐函数模型的系数，c
nm
和s
nm
为球谐展开的待求参数；表示测站与卫星j在第i观测历元之间的高度角，表示测站与卫星j在第i观测历元之间的方位角。
[0156]
经过公式(2)表示的基于球谐展开的标准单点定位观测方程，不涉及对流层延迟误差的改正模型，无需考虑不同模型在不同地区的适用问题，并且本实施例将球谐展开的系数当作待求参数直接求解，来改正与测站和卫星之间的高度角和方位角有关的误差项，主要包括对流层延迟误差，无需考虑气象信息，更加具有普适性。为了方便表示球谐函数展开，令：
[0157][0158]
则公式(2)简化为：
[0159]
由简化后的标准单点定位观测方程(3)，得到标准单点定位误差方程(4)，公式如下：
[0160][0161]
其中，v
ρ
表示伪距观测数据的改正值；
[0162]
表示卫星j在第i观测历元的伪距观测数据的改正值。
[0163]
i.3.基于公式(4)得到标准单点定位误差方程的线性化表达式，并进行简化处理，对简化后的标准单点定位误差方程的线性化表达式进行滑动解算。
[0164]
i.3.1.对标准单点定位误差方程(4)在测站的近似坐标(x0,y0,z0)处进行泰勒级数展开并保留一阶项，则标准单点定位误差方程的线性化表达式为：
[0165][0166]
公式(5)中，令公式(5)中，令
[0167][0168]
表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dx的系数，待求参数dx是测站近似坐标x0的改正数；表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dy的系数，待求参数dy是测站近似坐标y0的改正数；表示测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dz的系数，待求参数dz是测站近似坐标z0的改正数；表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数c
nm
的系数，表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数s
nm
的系数；tr表示测站接收机钟差，因为测站接收机钟差无法很好地被多项式拟合，所以将其作为待求参数进行解算，作为待求参数求解时，dtr的系数为1；d(tr)i表示在第i观测历元测站接收机钟差的待求参数。
[0169]
简化后的标准单点定位误差方程的线性化表达式，如公式(6)所示；
[0170][0171]
公式(6)中，i＝1,2,
…
,epoch，epoch表示最大观测历元数；表示由卫星j在第i观测历元伪距观测值、卫星j在第i观测历元的钟差以及测站近似坐标与卫星j在第i观测历元之间的几何距离计算出的常数项，
[0172]
由公式(6)能够看出，此时，有epoch
×
j个观测方程。而待求参数包括3个位置参数(dx,dy,dz)，epoch个接收机钟差dtr参数以及(nmax 1)2个球谐展开的参数c
nm
、s
nm
。
[0173]
当epoch充分多时，能够通过最小二乘求得所有未知数的最优解。因此只要改正项是可以线性描述的，就能通过增加观测历元的方式来增加观测方程数量，从而解算出这些未知参数。
[0174]
通过以上分析，为了确保未知参数的估计为最优解，本实施例选择合适的滑动窗口，对简化后的标准单点定位误差方程的线性表达式(6)进行滑动解算。
[0175]
设定滑动窗口选择i个历元，每个历元最多能够观测到j颗卫星，则在这个滑动窗口下所有可视卫星的伪距观测值组成imax个方程，其中，imax＝i
×
j。
[0176]
定义标准单点定位误差方程的系数矩阵、伪距观测值的改正数向量矩阵、常数项矩阵以及未知参数矩阵分别为矩阵b、v、l、x，则它们的表达式分别如下：
[0177][0178][0179][0180]
x＝[dxdydz(dtr)1…
(dtr)
ic00
…cnmaxnmaxsnmaxnmax
]
t
。
[0181]
则标准单点定位误差方程的线性化表达式(6)用一个矩阵形式表达，如公式(7)所示：
[0182]
v＝bx-l(7)
[0183]
在构建公式(7)的系数矩阵b时，只需要通过简单的数学计算，计算每一个待求参数(位置参数、测站接收机钟差参数及球谐展开的参数)的系数即可，计算过程简单，构建误差方程的速度迅速，实现标准单点定位的效率高。目前，与测站和卫星之间的高度角和方位角有关的误差项，如对流层延迟误差等还无法用模型精确改正，参考处理测站接收机钟差参数的方式，将其当做待求参数进行求解，是一种很好地方式。基于球谐展开的gnss标准单点定位方法的提出正是基于这一思想，利用球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，在标准单点定位方程中直接求解球谐展开的系数。
[0184]
未知参数x的最小二乘估值为：x＝(b
t
pb)-1bt
pl(8)
[0185]
公式(8)中，p为单位权矩阵。
[0186]
i.3.2.如果标准单点定位误差方程的系数矩阵b是良态，则用公式(8)求出的未知参数的解是最优无偏估计，当系数矩阵b是病态时，则利用最小二乘法求出的不再是最优解。
[0187]
因此，需要首先判断标准单点定位误差方程的系数矩阵b是否为病态，若系数矩阵b是病态，则执行步骤i.3.3；否则，若系数矩阵b为良态，执行步骤i.3.4。
[0188]
本实施例通过计算的系数阵b的条件数，设置条件数的经验阈值，并将条件数与经验阈值作比较，当条件数大于经验阈值时，则表明误差方程(7)为病态，否则为良态。
[0189]
i.3.3.为了解决误差方程的病态问题，本实施例提出了利用截断奇异值分解的方法为最小二乘谱修正迭代提供待求参数x的初值，基于最小二乘谱修正迭代计算待求参数x
的值。截断奇异值法的思想是计算矩阵s中奇异值的均值，将该均值作为截断奇异值的阈值；筛选出小于阈值的奇异值将其零化。由于截断奇异值分解方法是将方程的最小二乘法转换为直接解法，因而避免了病态方程的解过分偏离真实解，有效解决了误差方程系数矩阵b的病态问题。
[0190]
设系数矩阵b∈rn×m，rn×m表示n行m列的实数阵；则系数矩阵b的奇异值分解为：
[0191]
b＝usv
t
(9)
[0192]
公式(9)中，u∈n
×
n，v∈m
×
m，u、v均为正交矩阵，s∈n
×
m为一对角矩阵。
[0193]
系数矩阵b∈rn×m的截断奇异值矩阵bk定义为：
[0194][0195]
公式(10)中，矩阵s中最小的(r-k)个非奇异值用零代替，即被截断，且k≤r；r表示系数矩阵b的秩，k表示矩阵s中保留的奇异值数量；ui表示矩阵u对应的向量，vi表示矩阵v对应的向量；σi表示矩阵s中保留的奇异值。计算矩阵s中奇异值的均值，将该值作为截断奇异值的阈值；矩阵s中奇异值大于该阈值的保留，小于该阈值的进行零化处理，则公式(7)的截断奇异值解为：
[0196]
公式(11)中，
[0197]
根据最小二乘原理v
t
pv＝min，公式(8)写成以下形式：
[0198]
(b
t
pb)x＝(b
t
pl)(12)
[0199]
将公式(12)的等式左端和右端同时加未知参数kx，并化简得到公式(13)：
[0200]
(b
t
pb ki)x＝b
t
pl kx(13)
[0201]
由公式(13)整理得，最小二乘谱修正迭代公式为：
[0202]
x＝(b
t
pb ki)-1
(b
t
pl kx)(14)
[0203]
公式(14)中，k为任意实数。最小二乘谱修正迭代是一种基于最小二乘的迭代计算方法；基于最小二乘谱修正迭代计算未知参数x的解。求解过程如下：
[0204]
基于最小二乘谱修正迭代解算x包括迭代
①
和迭代
②
两个迭代过程；其中，设定用于判断迭代
①
是否收敛的第一比较阈值，设定用于判断迭代
②
是否收敛的第二比较阈值。
[0205]
迭代
①
第一次迭代时，设定k的初始值为1，将公式(11)得到的截断奇异值解作为公式(14)第一次迭代时x的初值，并赋值给公式(14)等式右边的x；
[0206]
每次迭代过程中，将上一次迭代时利用公式(14)得到的未知参数x的值，作为本次迭代过程中x的初值，并赋值给公式(14)等式右边的x；
[0207]
基于公式(14)求解得到每次迭代过程中未知参数x的值；
[0208]
将每次迭代过程中求解得到的x的值，与每次迭代过程中x的初值进行差值运算；
[0209]
若经过差值运算，得到的差值大于第一比较阈值，则表示迭代不收敛，需要改正最小二乘谱修正迭代中的k值，即将k值递增1，继续上述迭代
①
过程；
[0210]
若经过差值运算得到的差值小于或等于第一比较阈值，则跳出迭代
①
，进入迭代
②
过程。
[0211]
迭代
②
将利用公式(14)求解得到的未知参数x的值与第二比较阈值进行比较；
[0212]
若利用公式(14)求解得到的未知参数x的值大于第二比较阈值，则表示迭代
②
未收敛，此时，将位置参数(dx,dy,dz)分别加到测站的近似坐标(x0,y0,z0)上，更新测站的近似坐标，通过公式(5)对更新测站近似坐标后的标准单点定位误差方程进行线性化，通过公式(6)简化线性化后的标准单点定位误差方程；重新进入迭代
①
过程求解出未知参数x的值。
[0213]
若利用公式(14)求解得到的未知参数x的值小于或等于第二比较阈值，则表示迭代
②
收敛，跳出迭代；此时利用公式(14)得到的未知参数x的值，即未知参数x的解。
[0214]
求得未知参数x的解后，转到步骤i.3.5。
[0215]
i.3.4.利用最小二乘法求解未知参数x的解，转到步骤i.3.5。
[0216]
步骤i.3.4中，基于最小二乘法求解未知参数x的解的过程如下：
[0217]
设定用于判断迭代是否收敛的第三比较阈值；每次迭代过程中，利用公式(8)求解出未知参数x的值，将本次迭代得到的未知参数x的值与第三比较阈值进行比较；
[0218]
若利用公式(8)求解得到的未知参数x的值大于第三比较阈值，则表示迭代未收敛，此时，将位置参数(dx,dy,dz)，分别加到测站的近似坐标(x0,y0,z0)上，更新测站的近似坐标，执行公式(5)-公式(8)的计算过程，并重新求解未知参数x的值；
[0219]
若利用公式(8)求解得到的未知参数x的值小于或等于第三比较阈值，则表示迭代收敛，跳出迭代；此时，利用公式(8)求解得到的未知参数x的值，即未知参数x的解。
[0220]
i.3.5.根据计算得到的未知参数x的解，得到未知参数x中包含的参数值，即测站的位置参数(dx,dy,dz)、接收机钟差dtr以及球谐展开的系数c
nm
和s
nm
；
[0221]
其中，未知参数x中的位置参数(dx,dy,dz)分别加到测站的近似坐标(x0,y0,z0)上，即标准单点定位求出的地心地固坐标系下的测站坐标(x,y,z)。
[0222]
通过求解得到的球谐展开的系数c
nm
和s
nm
可直接用于后续标准单点定位。
[0223]
具体应用过程为：读取已有的球谐展开系数，通过球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，主要包括对流层延迟误差和电离层延迟误差；构建单点定位观测方程时，不再考虑与测站和卫星之间的高度角和方位角有关的误差项的改正问题，省略定位中对流层延迟误差和电离层延迟误差的求解过程，使得定位方程仅仅需要解算测站的位置参数、测站接收机钟差参数即可，极大地简化了标准单点定位流程。
[0224]
i.3.6.将测站坐标由地心地固坐标系转换到站心坐标系，实现gnss标准单点定位。
[0225]
将测站坐标由地心地固坐标系转换到站心坐标系的过程如下：
[0226]
地心地固坐标系转经纬高坐标系的公式为：
[0227]
式中，(x,y,z)为在地心地固坐标系的坐标，(lon,lat,alt)为在经纬高坐标系的
坐标；m表示基准椭球体的曲率半径，e表示椭球偏心率；根据测站近似坐标(x0,y0,z0)，利用公式(15)计算得到该测站近似坐标(x0,y0,z0)在经纬高坐标系下的对应坐标，即(lon0,lat0,alt0)。
[0228]
地心地固坐标系转换到站心坐标系的公式为：
[0229][0230]
式中，(e1,n1,u1)为以测站近似坐标(x0,y0,z0)为原点求出的(x,y,z)坐标的位置。
[0231]
s为坐标转换矩阵，
[0232]
根据以上坐标转换关系，将测站坐标由地心地固坐标系转换到站心坐标系。
[0233]
本实施例1中的gnss标准单点定位方法，基于球谐展开表示与测站和卫星间的方位角和高度角有关的误差项，从而能够快速、高效、准确地进行标准单点定位。
[0234]
实施例2
[0235]
本实施例2述及了一种基于球谐展开的gnss精密单点定位方法。如图2所示，该基于球谐展开的gnss精密单点定位方法包括如下步骤：
[0236]
ii.1.建立利用载波相位观测值进行的精密单点定位观测方程。
[0237]
读取观测数据、广播星历数据以及精密星历数据，对载波相位观测数据进行预处理，计算与测站和卫星间的高度角和方位角无关或关系不密切的误差；以上误差主要包括计算卫星钟差、探测周跳以及计算载波相位整周模糊度参数的初始值等。
[0238]
对于精密单点定位，广播星历的精度无法满足精密单点定位的要求，因此必须采用卫星精密星历计算卫星位置；利用二阶多项式估计的卫星钟差无法满足精密单点定位的精度要求；由于不同的卫星钟差均不相同，必须事先确定其大小，再代入到观测方程消除其影响。
[0239]
目前，igs提供的15min、5min采样间隔的事后精密星历和5min、30s采样间隔的卫星钟差产品，能够满足精密定位的精度要求。观测值的采样间隔一般都小于上述文件中数据的采样间隔，因此，需要内插计算得到每个历元对应的卫星位置和卫星钟差。
[0240]
地面接收机钟大多采用石英钟，精度低于卫星钟，钟差变化无规律性，无法通过多项式较好拟合，因此，将接收机钟差tr作为待求参数进行估计。
[0241]
建立的利用载波相位观测值进行的精密单点定位观测方程，如公式(1)所示。
[0242][0243]
公式(1)中，i表示第i观测历元；j表示卫星编号；表示载波相位观测值，将载波相位原始观测值乘上相应的波长λ，得到以距离表示的载波相位观测数据，称为等效距离；表示卫星j在第i观测历元的载波相位观测值，即等效距离；表示接收机的位置与卫星j在第i观测历元位置之间的几何距离；tr表示接收机钟差，(tr)i表示第i观测历元接收机钟
差；ts表示卫星钟差，表示卫星j在第i观测历元的钟差；表示卫星j在第i观测历元的信号传播路径上的对流层延迟误差；表示卫星j在第i观测历元的信号传播路径上的电离层延迟误差；c表示光速；λ表示载波的波长；n表示卫星相位观测数据的整周模糊度参数，表示卫星j在第i观测历元的相位观测数据整周模糊度参数；表示载波相位观测数据残差。定义卫星j在第i历元观测瞬间的空间三维坐标为测站的近似坐标为(x0,y0,z0)，则卫星到测站近似位置的几何距离表示为：
[0244][0245]
gnss载波相位观测数据类型分为单频观测数据和双频观测数据。
[0246]
若进行标准单点定位的gnss接收机为双频接收机，则电离层延迟误差可以通过双频消电离层组合来消除，具体的组合形式如下：式中，表示载波相位消电离层组合观测值(等效距离)；f1、f2表示相应观测值的频率；分别表示两种频率对应的载波相位观测值(等效距离)。若观测数据为双频观测值时，则上述公式(1)中，电离层延迟误差的一阶主项通过双频消电离层组合消除，电离层延迟误差项为0。
[0247]
此时，本实施例2中的精密单点定位观测方程中，球谐展开主要用于表示对流层延迟误差项，球谐展开的阶次与单频载波相位观测值构建的观测方程有所不同。
[0248]
在精密单点定位中，利用载波相位观测值组成gf(几何无关)组合观测值，进行整周模糊度的解算和周跳的探测。具体的组合形式如下：
[0249]
式中，表示gf组合载波相位观测值(等效距离)；分别表示两种频率对应的载波相位观测值；λ1、λ2分别表示两种频率对应的载波相位的波长；f1、f2表示两种载波相位的频率。gf组合观测值与几何距离无关，因此，该组合适合用来探测周跳。
[0250]
ii.2.基于球谐展开表示与测站和卫星间的高度角和方位角有关的误差项，其中，误差项主要包括对流层延迟误差以及电离层延迟误差。在精密单点定位中，对流层延迟误差和电离层延迟误差都是卫星信号穿过大气层时产生的一种信号延迟误差，都与卫星信号的传播有关；而卫星信号的传播路径又与测站和卫星的位置有关，利用测站坐标和卫星位置可以计算出测站和卫星之间的高度角和方位角，应用球谐展开，就可以表示与测站和卫星间的高度角和方位角有关的误差项，主要包括对流层延迟误差和电离层延迟误差，以消除或最大程度的削弱其对精密单点定位结果的影响。经过球谐展开后的精密单点定位观测方程表示为：
[0251][0252]
式(2)中，n为球谐展开的阶数，m为球谐展开的次数，nmax为球谐展开的最大阶数；
p
nm
(cose)表示n阶m次的缔合勒让德多项式；c
nm
和s
nm
表示n阶m次的球谐展开的系数；表示测站与卫星j在第i观测历元之间的高度角，表示测站与卫星j在第i观测历元之间的方位角。公式(2)提出的基于球谐展开的精密单点定位观测方程，不涉及对流层延迟误差的改正模型，无需考虑不同模型在不同地区的适用问题，并且本实施例将球谐展开的系数当作待求参数直接求解，来改正与测站和卫星之间的高度角和方位角有关的误差项，主要包括对流层延迟误差，无需考虑气象信息，更加具有普适性；为了方便表示球谐函数展开，令：
[0253][0254]
则公式(2)简化为：
[0255]
由公式(3)中的精密单点定位观测方程，得到精密单点定位误差方程：
[0256][0257]
其中，表示载波相位观测数据(等效距离)的改正值。
[0258]
表示第i观测历元卫星j的载波相位观测数据(等效距离)的改正值。
[0259]
ii.3.基于公式(4)得到精密单点定位误差方程的线性化表达式，并简化处理，然后对简化后的精密单点定位误差方程的线性化表达式进行滑动解算。
[0260]
ii.3.1.对精密单点定位误差方程(4)在测站的近似坐标(x0,y0,z0)处进行泰勒级数展开并保留一阶项，则精密单点定位误差方程的线性化表达式为：
[0261][0262]
公式(5)中，表示卫星j在第i历元的模糊度初值；表示卫星j在第i历元的相位观测数据整周模糊度参数，dn表示相位观测数据整周模糊度参数的改正数，表示卫星j在第i历元的相位观测数据整周模糊度参数的改正数，在精密单点定位中，由于整周模糊度参数不能被很好地模型化，因此将其改正数作为待求参数进行求解，作为待求参数求解时，dn
ij
的系数为1；d(tr)i表示在第i观测历元测站接收机钟差的待求参数，与整周模糊度的处理方式相似，因为测站接收机钟差无法很好地被多项式拟合，因此将其作为待求参数进行求解，作为待求参数求解时，dtr的系数为1。
[0263]
令令
[0264][0265]
为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dx的系数，待求
参数dx是测站近似坐标x0的改正数；为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dy的系数，待求参数dy是测站近似坐标y0的改正数；为测站近似坐标与卫星j在第i观测历元坐标计算出的待求参数dz的系数，待求参数dz是测站近似坐标z0的改正数；a
nm
表示球谐展开待求参数c
nm
的系数，b
nm
表示球谐展开待求参数s
nm
的系数；表示由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数c
nm
的系数，为由测站近似坐标与卫星j在第i观测历元坐标计算出的球谐展开的待求参数s
nm
的系数。
[0266]
简化后的精密单点定位误差方程的线性化表达式，如公式(6)所示：
[0267][0268]
公式(6)中，i＝1,2,,
…
,epoch，epoch表示最大观测历元；
[0269]
表示由卫星j在第i观测历元载波相位观测值、卫星j在第i观测历元的钟差、卫星j在第i历元的相位观测数据整周模糊度初始值以及测站近似坐标与卫星j在第i观测历元之间的几何距离计算出的常数项，由公式(6)可以看出，此时有epoch
×
j个观测方程。而待求参数包括：3个位置参数(dx,dy,dz)，epoch个接收机钟差参数dtr，j个模糊度参数dn以及(nmax 1)2个球谐展开的系数c
nm
、s
nm
。
[0270]
当epoch充分多时，能够通过最小二乘求得所有未知数的最优解，只要改正项是可以线性描述的，就可以通过增加观测历元的方式来增加观测方程数量，从而解算出这些未知参数。
[0271]
设定滑动窗口选择i个历元，每个历元最多能够观测到j颗卫星，则在这个滑动窗口下所有可视卫星的载波相位观测值组成imax个方程，其中，imax＝i
×
j。
[0272]
定义精密单点定位误差方程的系数矩阵、载波相位观测值的改正数向量矩阵、常数项矩阵以及未知参数矩阵分别为矩阵b、v、l、x，则它们的表达式分别如下：
[0273][0274][0275][0276][0277]
表示第1颗卫星在第1个观测历元的整周模糊度参数，表示第j颗卫
星在第i个观测历元的整周模糊度参数，精密单点定位中求得的是整周模糊度初始值的改正数。
[0278]
则精密单点定位误差方程的线性化表达式(6)用一个矩阵形式表达，如公式(7)所示；
[0279]
v＝bx-l(7)
[0280]
构建精密单点定位误差方程式(7)的系数矩阵b时，只需要通过简单的数学计算，计算每一个待求参数(位置参数、测站接收机钟差参数、整周模糊度参数及球谐展开的参数)的系数即可，计算过程简单，构建误差方程的速度迅速，实现精密单点定位的效率高。目前，与测站和卫星之间的高度角和方位角有关的误差项，如对流层延迟误差等还无法用模型精确改正，参考处理测站接收机钟差参数的方式，将其当做待求参数进行求解，是一种很好地方式；基于球谐展开的gnss精密单点定位方法的提出正是基于这一思想，利用球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，在精密单点定位方程中直接求解球谐展开的系数。未知参数x的最小二乘估值为：x＝(b
t
pb)-1bt
pl(8)
[0281]
公式(8)中，p为单位权矩阵。
[0282]
ii.3.2.如果精密单点定位误差方程的系数矩阵b是良态，则用公式(8)求出的未知参数的解时最优无偏估计，当系数矩阵b是病态时，则利用最小二乘法求出的不再是最优解。
[0283]
因此，需要首先判断精密单点定位误差方程的系数矩阵b是否为病态，若系数矩阵b是病态，则执行步骤ii.3.3；否则，系数矩阵b为良态，执行步骤ii.3.4。
[0284]
本实施例通过计算精密单点定位误差方程的系数阵b的条件数，设置条件数的经验阈值，将条件数与经验阈值作比较，当条件数大于经验阈值时，则误差方程为病态。
[0285]
ii.3.3.为了解决误差方程的病态问题，本实施例提出了利用截断奇异值分解的方法，为最小二乘谱修正迭代提供待求参数x的初值，基于最小二乘谱修正迭代计算待求参数x的值。由于截断奇异值分解方法是将方程的最小二乘法转换为直接解法，因而避免了病态方程的解过分偏离真实解，有效解决了误差方程系数矩阵b的病态问题。
[0286]
设系数矩阵b∈rn×m，rn×m表示n行m列的实数阵，则系数矩阵b的奇异值分解为：
[0287]
b＝usv
t
(9)
[0288]
公式(9)中，u∈n
×
n，v∈m
×
m，u、v均为正交矩阵，s∈n
×
m为一对角矩阵。
[0289]
系数矩阵b∈rn×m的截断奇异值矩阵bk定义为：
[0290][0291]
公式(10)中，矩阵s中最小的(r-k)个非奇异值用零代替，即被截断，且k≤r；r表示系数矩阵b的秩，k表示矩阵u中保留的奇异值数量，ui表示矩阵u对应的向量，vi表示矩阵v对应的向量；σi表示矩阵u中保留的奇异值。计算矩阵s中奇异值的均值，将该值作为截断奇异值的阈值；矩阵s中奇异值大于该截断奇异值的阈值的保留，小于该截断奇异值的阈值的进行零化处理。则公式(7)的截断奇异值解为：
[0292]
[0293]
公式(11)中，
[0294]
最小二乘谱修正迭代是一种基于最小二乘的迭代计算方法。基于最小二乘谱修正迭代计算待求参数x；根据最小二乘原理v
t
pv＝min，公式(8)写成以下形式：
[0295]
(b
t
pb)x＝(b
t
pl)(12)
[0296]
将公式(12)的等式左端和右端同时加未知参数kx，并化简得到公式(13)：
[0297]
(b
t
pb ki)x＝b
t
pl kx(13)
[0298]
由公式(13)整理得，最小二乘谱修正迭代公式为：
[0299]
x＝(b
t
pb ki)-1
(b
t
pl kx)(14)
[0300]
公式(14)中，k为任意实数；基于最小二乘谱修正迭代计算未知参数x的解。
[0301]
基于最小二乘谱修正迭代公式计算x的解的过程如下：
[0302]
基于最小二乘谱修正迭代解算x包括迭代
①
和迭代
②
两个迭代过程；其中，设定用于判断迭代
①
是否收敛的第一比较阈值，设定用于判断迭代
②
是否收敛的第二比较阈值；
[0303]
迭代
①
第一次迭代时，设定k的初始值为1，将公式(11)得到的截断奇异值解作为公式(14)第一次迭代时x的初值，并赋值给公式(14)等式右边的x；
[0304]
每次迭代过程中，将上一次迭代时利用公式(14)得到的未知参数x的值，作为本次迭代过程中x的初值，并赋值给公式(14)等式右边的x；
[0305]
基于公式(14)求解得到每次迭代过程中未知参数x的值；
[0306]
将每次迭代过程中求解得到的x的值，与每次迭代过程中x的初值进行差值运算；
[0307]
若经过差值运算，得到的差值大于第一比较阈值，则表示迭代不收敛，需要改正最小二乘谱修正迭代中的k值，即将k值递增1，继续上述迭代
①
过程；
[0308]
若经过差值运算得到的差值小于或等于第一比较阈值，则跳出迭代
①
，进入迭代
②
过程；
[0309]
迭代
②
将利用公式(14)求解得到的未知参数x的值与第二比较阈值进行比较；
[0310]
若利用公式(14)求解得到的未知参数x的值大于第二比较阈值，则表示迭代
②
未收敛，此时，将位置参数(dx,dy,dz)分别加到测站的近似坐标(x0,y0,z0)上，更新测站近似坐标，通过公式(5)对更新测站近似坐标后的精密单点定位误差方程进行线性化，通过公式(6)简化线性化后的精密单点定位误差方程；重新进入迭代
①
过程求解出未知参数x的值；
[0311]
若利用公式(14)求解得到的未知参数x的值小于或等于第二比较阈值，则表示迭代
②
收敛，跳出迭代；此时利用公式(14)得到的未知参数x的值，即待求解的未知参数x的解。
[0312]
未知参数x的解后，转到步骤ii.3.5。
[0313]
ii.3.4.利用最小二乘法求解未知参数x的解，转到步骤ii.3.5。
[0314]
步骤ii.3.4中，基于最小二乘法求解未知参数x的解的过程如下：
[0315]
设定用于判断迭代是否收敛的第三比较阈值；每次迭代过程中，利用公式(8)求解出未知参数x的值，将本次迭代得到的未知参数x的值与第三比较阈值进行比较；
[0316]
若利用公式(8)求解得到的未知参数x的值大于第三比较阈值，则表示迭代未收敛，此时，将位置参数(dx,dy,dz)分别加到测站的近似坐标(x0,y0,z0)上，更新测站近似坐标；
[0317]
执行公式(5)-公式(8)的计算过程，并重新求解未知参数x的值；
[0318]
若利用公式(8)求解得到的未知参数x的值小于或等于第三比较阈值，则表示迭代收敛，跳出迭代；此时，利用公式(8)求解得到的未知参数x的值，即待求解的未知参数x的解。
[0319]
ii.3.5.根据计算得到的未知参数x的解，得到未知参数x中包含的参数值，即测站的位置参数(dx,dy,dz)、接收机钟差dtr、模糊度参数dn以及球谐展开的系数c
nm
和s
nm
；
[0320]
其中，位置参数(dx,dy,dz)分别加到测站的近似坐标(x0,y0,z0)上，即精密单点定位求出的地心地固坐标系下的测站坐标(x,y,z)。通过求解得到的球谐展开的系数c
nm
和s
nm
，可以直接用于后续精密单点定位，具体应用过程如下：
[0321]
读取已有的球谐展开系数，通过球谐展开表示与测站和卫星之间的高度角和方位角有关的误差项，主要包括对流层延迟误差和电离层延迟误差，构建的精密单点定位观测方程可以不再考虑与测站和卫星之间的高度角和方位角有关的误差项的改正问题，省略定位中的对流层延迟误差和电离层延迟误差的求解过程，这将使得定位方程仅仅需要解算测站的位置参数、测站接收机钟差参数以及整周模糊度参数即可，极大地简化了精密单点定位流程。
[0322]
ii.3.6.将测站坐标由地心地固坐标系转换到站心坐标系，实现gnss精密单点定位。其中，本实施例2中坐标转换过程与上述实施例1中坐标转换过程相同，此处不再赘述。
[0323]
本实施例2中述及的gnss精密单点定位方法，基于球谐展开表示与测站和卫星间的方位角和高度角有关的误差项，从而能够快速、高效、准确地进行精密单点定位。
[0324]
下面对利用本发明中基于球谐展开的gnss单点定位方法所能达到的精度进行说明：
[0325]
以上海佘山站为例：
[0326]
获取2019年1月7日佘山站的gnss观测数据和卫星星历数据。
[0327]
gnss观测数据采样间隔为30s，观测数据选择gps卫星的双频伪距观测值以及双频载波相位观测值，站点的参考坐标为igs提供的站点坐标。
[0328]
1.首先利用本发明中基于球谐展开的gnss标准单点定位方法获取站点位置信息。
[0329]
观测数据选择gps卫星双频伪距观测值，则电离层延迟误差通过伪距线性组合消除。
[0330]
通过球谐展开表示与卫星和测站间的高度角和方位角有关的误差项，主要是对流层延迟误差。球谐展开各阶次的系数、以及测站的接收机钟差需要当作待求参数求解。
[0331]
将球谐函数展开到3阶，滑动窗口选择10个历元。待求参数包括：点位坐标3个方向的改正值，10个接收机钟差(每个历元一个接收机钟差参数)，球谐展开到3阶3次的16个系数：a
00
、a
10
、a
11
、b
11
、a
20
、a
21
、b
21
、a
22
、b
22
、a
30
、a
31
、b
31
、a
32
、b
32
、a
33
、b
33
。
[0332]
卫星的截止高度角设置为5
°
，最小二乘谱修正迭代的系数k设为2。
[0333]
gnss标准单点定位方法的定位结果如图3所示，统计结果见表1。
[0334]
表1标准单点定位精度(单位/m)
[0335][0336][0337]
表2 bernese定位精度(单位/m)
[0338] enumax6.09175.75159.7484min-3.8416-5.8132-11.1756mean0.1268-0.12590.5520std1.27441.60023.0886rms1.28071.60523.1375
[0339]
通过表1统计结果不难发现，本发明标准单点定位方法得到的定位结果在e方向上的rms优于0.39m，在n方向上的rms优于0.40m，在u方向上的rms优于1.43m。
[0340]
由瑞士伯尔尼大学天文研究所开发的国际知名gnss高精度数据处理软件bernese软件解算的标准单点定位结果，如表2所示。
[0341]
可见，gnss高精度数据处理软件bernese软件解算的定位结果在e方向上的rms优于1.28m，在n方向上的rms优于1.61m，在u方向上的rms优于3.14m。
[0342]
从统计表中看出，本发明标准单点定位结果与bernese软件标准单点定位结果精度相当。
[0343]
2.利用本发明中基于球谐展开的gnss标准单点定位方法获取站点位置信息。
[0344]
观测数据选择gps卫星双频载波相位观测值，则电离层延迟误差通过载波相位线性组合消除。通过球谐展开表示与卫星和测站间的高度角和方位角有关的误差项，主要是对流层延迟误差；球谐展开各阶次的系数需要当作待求参数求解。
[0345]
每颗卫星载波相位观测数据数据的整周模糊度以及测站的接收机钟差也需要作为待求参数求解。将球谐函数展开到6阶，滑动窗口选择15个历元。
[0346]
待求参数包括：点位坐标3个方向的改正值，15个接收机钟差(每个历元一个接收机钟差参数)，每个历元下每颗卫星的模糊度参数，球谐展开到6阶6次的49个系数：a
00
、a
10
、a
11
、b
11
、a
20
、a
21
、b
21
、a
22
、b
22
、a
30
、a
31
、b
31
、a
32
、b
32
、a
33
、b
33
、a
40
、a
41
、b
41
、a
42
、b
42
、a
43
、b
43
、a
44
、b
44
、a
50
、a
51
、b
51
、a
52
、b
52
、a
53
、b
53
、a
54
、b
54
、a
55
、b
55
、a
60
、a
61
、b
61
、a
62
、b
62
、a
63
、b
63
、a
64
、b
64
、a
65
、b
65
、a
66
、b
66
。
[0347]
卫星的截止高度角设置为5
°
，最小二乘谱修正迭代的系数k设为3。
[0348]
gnss精密单点定位方法的结果如图4所示，统计结果见表3。
[0349]
表3精密单点定位精度(单位/m)
[0350]
[0351][0352]
表4 bernese定位精度(单位/m)
[0353] enumax0.02740.03010.0461min-0.0178-0.0209-0.0738mean0.00360.0048-0.0082std0.00750.00920.0200rms0.00830.01040.0216
[0354]
通过表3统计结果不难发现，本发明方法得到的定位结果在e方向上的rms优于0.02m，在n方向上的rms优于0.03m，在u方向上的rms优于0.04m。由瑞士伯尔尼大学天文研究所开发的国际知名gnss高精度数据处理软件bernese软件解算的标准单点定位结果如表4所示。可见，gnss高精度数据处理软件bernese软件解算的定位结果在e方向上的rms优于0.009m，在n方向上的rms优于0.02m，在u方向上的rms优于0.03m。
[0355]
从统计表中看出，本发明精密单点定位结果与bernese软件精密单点定位结果精度相当。
[0356]
本发明基于球谐展开的gnss单点定位方法，利用观测数据和卫星星历数据高效快速获得测点的点位信息，相比于现有的利用经验模型进行对流层延迟改正方法，具有解算效率高，无需提供气象文件，操作简便、工作难度小，测量结果可靠，测量精度高等优势。
[0357]
当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。